Classification of ground deformation using sentinel-1 persistent scatterer interferometry time series

ABSTRACT Displacement time series (TS) provides temporal and spatial information related to ground deformation. This study aims to investigate temporal behavior of ground deformation TS, including classification of displacement trends and periodicity evaluation, which ease the interpretation of movements. To this end, we propose several modifications to an existing automatic classification workflow of Persistent Scatterers Interferometry (PSI) TS using new tests to classify ground deformations into seven main trends: Stable, Linear, Quadratic, Bilinear, Phase Unwrapping Errors (PUE), Discontinuous with constant and different velocities. We illustrate our approach over 1500 km2 of the Granada region and the metropolitan area of Barcelona, which were monitored using Sentinel-1 images and a PSI technique. This study provided the spatial distribution of different ground movement types and was useful to detect several TS anomalies due to PUE. The proposed approach also identified stable targets, which were wrongly classified as moving scatterers by the existing classification method. A periodicity analysis was finally performed using the Welch’s power spectral density estimator to investigate seasonal and yearly fluctuations. The method was validated using simulated data, where the classified TSs characterized by probable phase unwrapping errors were verified by PSI experts. The overall classification accuracy was 77.8%, indicating that the proposed method has a considerable TS classification potential.


Introduction
This paper proposes a new procedure to the existing workflow (Berti et al. 2013) to analyze the temporal behavior of ground deformations by implementing statistical techniques on Differential Interferometric Synthetic Aperture Radar (DInSAR) time series (TS) products. The analysis highlights the trends and behaviors of the deformations, along with the physicaltemporal relationship of the movements. The analysis can be used to classify Ground Motion Services (GMS) TS. These novel services, which are based on Persistent Scatterer Interferometry (PSI), Small BAseline Subset (SBAS), and Distributed Scatterers (DS), provide consistent and reliable information on ground movements related to artificial or natural targets (Crosetto et al. 2020). The GMS application fields are broad and include natural and man-induced geohazard risk assessment, land management, urban planning, infrastructure development and management, mining and other natural resources extraction, dam, groundwater monitoring, etc. (EEA 2020).
The PSI techniques was first introduced by Rocca (2001, 2000), and has been subsequently developed by several other groups (see Crosetto et al. (2016) for a review). The deformation time series and velocity estimation at PS (Persistent Scatterer) points are main outcomes of the PSI technique (Crosetto et al. 2016). The PSI method can measure displacement of the available PSs, which typically corresponds to buildings, structures, infrastructures, outcrops, and exposed rocks.
The PSI TS contains the accumulated deformation in correspondence to the time-ordered SAR images, starting from a reference image. The accumulated displacements are estimated along the satellite Line of Sight (LoS) direction. A TS is defined as a set of quantitative observations arranged in chronological order (Kirchgässner, Wolters, and Hassler 2012), which generally includes three main components (Kirchgässner, Wolters, and Hassler 2012;Puech et al. 2020): (i) the trend component (i.e. a long-term development), (ii) the periodic component (i.e. a cyclical component with periods of one or more than a year or a component that includes the ups and downs within a year, called the seasonal cycle), and (iii) the random component (i.e. there is no long-term trend nor the cyclic and the seasonal component). Most of the published deformation TS studies have been conducted in the amplitude domain (i.e. the accumulated deformation) and in a univariate context (i.e. no correlation with other time series is considered).
Empirical methods have been widely used to categorize TS deformations using statistical criteria. Cigna et al. (2011) labeled each displacement TSs as "Affected targets" and "Unaffected targets." Deformation rates of TS sub-samples were calculated using a simple linear regression to evaluate deformation, velocity, and acceleration patterns. Later, Cigna, Tapete, and Casagli (2012) developed a method to calculate two indexes to depict the deviation using a predefined deformation model. The first index identified precursor movements that can be assigned to significant events, such as accelerations recorded before landslide failures or volcanic eruptions. The second index identified sudden changes in displacement at specific dates, such as those corresponding to tectonic events or structural collapses. Then, Berti et al. (2013) proposed an automatic classification method of PSI TS based on a conditional sequence of statistical tests. The TSs were classified into distinctive predefined classes, which described various types of ground deformation. It was concluded that the method is independent of the dataset. Notti et al. (2015) considered the Berti et al. (2013) and Cigna et al. (2012Cigna et al. ( , 2011 methods to highlight the ability of automatic detection of motion trends in the interpretation of DInSAR data. Recently, ancillary information has been utilized to overcome the limitations of empirical methods. More recent studies Solari et al. 2019) generated continuously updated ground deformation data and highlighted many anomalous trends and/or acceleration affecting the area of interest. Intrieri et al.
(2018) applied a detailed TS analysis on one of the accelerating PS to predict the time of a slope failure. Then in 2019, Raspini et al. (2019) presented the results of continuous monitoring of ground deformation in the Tuscany region based on the difference of deformation rates with respect to the slope and aspect derived from a Digital Elevation Model (DEM).
Several studies have provided periodicities analysis to separate deformation trends from seasonal deformations and understand relationships between triggering factors and ground displacements. For instance, Tomás et al. (2016) identified localized intermittent periodicities using three wavelet tools to interpret InSAR TS information for a landslide. Schlögl, Widhalm, and Avian (2021) have recently performed a decomposition into trend, periodic, and random components on PSI TS to analyze the deformation of a bridge.
In addition to utilizing ancillary information, Machine Learning (ML)-based methods have been recently studied to analyze TSs (Anantrasirichai et al., 2020;Ansari et al. 2021;Hakim, Achmad, and Lee 2020;Tiwari, Narayan, and Dikshit 2020;Valade et al. 2019). For instance, Valade et al. (2019) and Anantrasirichai et al. (2020) have investigated the pros and cons of Convolutional Neural Networks (CNN) in large deformations detection using multisensory satellite-based data. This technique has also been employed for the selection of PS pixels (Tiwari, Narayan, and Dikshit 2020). Hakim, Achmad, and Lee (2020) and Ranjgar et al. (2021) have mapped land subsidence susceptibility using functional and optimized ML algorithms. Finally, an autoencoder clustering model, deep Long Short Term Memory (LSTM), has been proposed to exploit temporal relations for efficient TS mining (Ansari et al. 2021).
According to the state-of-the-art, the TS types have not been classified using large datasets. Furthermore, the Phase Unwrapping Errors (PUE) have not been addressed in any classification studies. This error is due to a wrong reconstruction of the unwrapped phases and cause jumps in TSs that are a multiple of half the wavelength. In this study, we present an automatic approach to classify TSs in seven classes, which is a modification of the method proposed by Berti et al. (2013), and to find those TSs affected by PUE. We also propose a verification step to identify stable (non-moving) targets, which are incorrectly categorized as unstable points. Additionally, the PSI TS are evaluated to identify periodic variations. Section 2 describes the characteristics of the two study areas, along with the Sentinel-1 data information. The method section (Section 3) contains the definition of TS trends, the classification steps, and the periodicity analysis. Classified maps over two study areas and several cases of displacements are provided in Section 4, followed by a discussion about the advantages and limitations of the study. The paper finishes with the conclusion.

Study areas
In this study, two areas located in Spain have been selected (Figure 1(a)), which a PSI technique so far exploited TS products of both study areas. The first study area, called the Granada region, is located south of Granada and Malaga provinces (Figure 1(b)). Three significant movements cases have been recently distinguished (Reyes-Carmona et al. 2020), utilized to validate the methodology. It has an area of 575 km 2 . The second study area lies in the Barcelona metropolitan areas (Figure 1(c)), where around 900 km 2 have been analyzed for potential deformations and the reliability of the proposed methodology in displacements TS classification.

Data
Sentinel-1 globally provides continuous acquisition of C-band SAR images, which are freely accessible. Sentinel-1 consists of a constellation of two satellites, Sentinel-1A and Sentinel-1B, which cover Europe with a revisit time of six days that enables Europe-wide displacement monitoring. In this study, 138 and 249 of Sentinel-1 images were exploited over the Granada region and Barcelona's metropolitan area, respectively. Additionally, using the PSI chain of the Geomatics Division of CTTC (PSIG) (Devanthéry et al. 2014), deformation TSs of each study area have been extracted.
Six thousand six hundred and sixty-four interferograms were generated from 138 SAR images over Granada from 10 March 2015 to 26 September 2018 and 453,315 PS were extracted. The second dataset includes 249 images, which were acquired from 6 March 2015 to 19 June 2020. In this area, 1,623,767 PS points were generated from 1245 interferograms.

Method
The main goal of this study is to provide an advanced procedure to classify land deformation TSs, investigate TS with PUE, and develop an automatic system for processing large TS datasets. To this end, various statistical analyses have been implemented using the C++ language and MATLAB TM .
Ground displacements were categorized by Berti et al. (2013) in six types. In this study, an additional PUE class has been added to those six classes ( Figure 2): • Stable: no deformation occurs. There is no trend (e.g. linear, quadratic, and bilinear) in the analyzed points (Figure 2(a)). The TS consists of small random variations of the state of the ground surface.
• Linear: the TS is characterized by a constant slope (Figure 2(b)). It indicates long timescale movements, such as the steady motion of inactive landslides, natural subsidence, and creep (Berti et al. 2013). • Quadratic: the displacements follow a quadratic law (Figure 2(c)). Nonlinear classes stand for trends characterized by changes in the deformation rate. The Quadratic displacement shows a continuous movement, such as the deformation in mining areas, causing a progressive change in the rate of TS and revealing a continuous and nonlinear movement (Wang, Deng, and Zheng 2020). • Bilinear: the TS is characterized by a breakpoint, which divides the TS into two periods, where each period has a different linear rate (Figure 2 (d)). The different rates point to an abrupt change after the breakpoint, which can indicate a change point in deformations. In the case of a significant diversity after breakpoints, seismic activities can cause abrupt changes in the TS rates. Figure 3 shows the Berti et al. (2013) and the proposed workflows. The basic idea to automatically classify PSI TS based on conditional sequences of statistical tests was proposed by Berti et al. (2013), where the results were finally categorized into three classes (uncorrelated, linear, and non-linear trends), instead of the above listed six classes. It also opened a procedure to enhance the radar interpretation capability using an automatic classifier. In this research, the following improvements are proposed: (i) detecting PU errors in four stages using various statistical tests, (ii) utilizing statistical tests to reclassify stable points correctly, which were wrongly classified as Linear, Quadratic, and Bilinear, and (iii) using more conditions and additional statistical tests after the linearity check. The steps of the proposed approach are as follows: a) First, several intervals and thresholds are examined to explore a practical conditional test for identifying the PSs characterized by phase unwrapping errors (PU Error_1). b) A linear regression is fitted over the remaining PSs to identify Stable points. c) A segmented regression tests the unstable points to check for abrupt changes in TSs. d) If there is no breakpoint inside a TS, a quadratic regression is fitted to classify it as Linear or Quadratic. This primary classification is then tested against the second and third stages of PU errors (PU Error_2 and PU Error_3), along with uncorrelation check (Uncorrelated_1 and Uncorrelated_2) to distinguish PUE and Stable PSs. e) A discontinuity test is applied to separate Bilinear points. The uncorrelation and PU errors tests (Uncorrelated_3 and PU Error_4) are also applied at this level. f) Finally, a test for the equality of slopes is computed to split TSs with constant or varying velocities. The overview of the entire process is given in Table 1 and explained in the following sections comprehensively.

First stage for detecting phase unwrapping errors
The TSs are checked by a condition using the running mean of every seven samples for probable PUE (PU Error_1). A PU error is a multiple of half the Figure 3. Flowchart of the PS TS classification procedure. The left part is the workflow proposed by (Berti et al. 2013), and the right one contains the improvements presented in this paper.
wavelength, approximately 28.3 mm in Sentinel-1 SAR images. It causes jumps in TS values, leading to misclassification. In some cases, this amount is affected by residual Atmospheric Phase Screen (APS) and noise. In this study, values with probable PUE (PU Error_1) were identified based on the average moving technique (Equation 1): where y k , for 7 � k � n À 7 indicates the k th value of the TS and n stands for the number of images in the TS. The statistics, window size, M and the PUE threshold parameters, J were chosen based on various experiments (Table 2) and experts judgment. Accordingly, we found that the mean averaging outperforms the median ones. Moreover, from those eight different windows and seven thresholds considered, window size W ¼ 7 and threshold J ¼ 19mm are the optimal combinations for the first stage of PUE classification. Detailed results on the outputs of the experiment are presented in Table 2. Table 2 consists of Precision (P) and False Alarm Rate (FAR) percentages of all experiments, indicating the portion of correctly classified PUE and incorrectly classified other samples as PUE, respectively.

Linear regression
In this step, we pursued the technique proposed by Berti et al. (2013): every TS is fitted by a linear regression model. We are implementing an ANalysis Of VAriance (ANOVA) Fisher test to calculate the F statistic for each TS. A one-way ANOVA obtains information about the relationship between the dependent (e.g. the rate of deformations) and independent variables (e.g. deformations). The F value also indicates whether the correlation is significant enough to reject the null hypothesis. The null hypothesis stands for an absence of linear correlation in the  TS. In this work, the probability has been chosen 0.01. This value indicates a 1% risk of concluding that a linear correlation exists when there is no actual correlation (Frost 2021), and it is used to find the corresponding critical value of the underlying F distribution (e.g. via F distribution tables) (Christensen 2016; Montgomery and Runger 2018). Therefore, the points where the F statistic is smaller than the F distribution value have a slope close to zero and are classified as Stable (i.e. the linearity test fails).

Breakpoint detection
The TSs with a linear correlation are analyzed for detecting breakpoints using the procedure proposed by Berti et al. (2013). Each TS is split at moving breakpoints into two segments (the minimum number of data points of each segment is five), and each segment is fitted with linear regression. In the first step, the first segment consists of 5 values (data points from 1 to 5) and the second one includes the remaining data points (from 6 th to the last one). A bilinear model is applied over the TS to obtain the breakpoint occurring as an abrupt change in the slope of TS (Berti et al. 2013). A threshold is then defined using an evidence ratio computed by Wagenmakers and Farrell (2004) method to investigate a significant breakpoint. The evidence ratio also describes the evidence about fitted models and how a variable influences the detection probability or occupancy. This ratio is calculated using the Bayesian Information Criterion (BIC) (Equation 2) of a two-line unconstrained model (k ¼ 3), a quadratic model (k ¼ 2), and a linear model (k ¼ 1) over the TS.
where the Residual Sum of Squares (RSS) measures the discrepancy between the real data and estimated values. The BIC establishes selective criteria, where the model with lowest BIC can be preferred as the best fitted model (Claeskens and Hjort 2008). Thus, an evidence value is calculated based on the computed BICs to decide whether TS contains a breakpoint or not. To this end, the evidence ratios (Equation 3) are calculated using the weights of the two-line unconstrained (w 1 ), linear (w 2 ), and quadratic (w 3 ) models.
, À BIC twoÀ line , and Δ 3 ¼ BIC Quadratic À BIC twoÀ line . Because TSs displacements consist of physical and natural effects, a restrictive criterion has been proposed (Berti et al. 2013), using a normalized form of the relative likelihoods. Therefore, a TS with the evidence ratio of weights greater than and equal to the selected threshold reveals a breakpoint in the TS (i.e. a TS with bilinear regression). A PS that does not satisfy the above condition is tested against a quadratic regression in the following section.

Linear and quadratic TS
At this level, there are points with linear and quadratic correlations and with no first stage of PUE TSs. The TS must be fitted by a quadratic regression to test for the significance of the quadratic trend against the linear one (Berti et al. 2013;Davis 1986). Many TSs have been found as misclassified points using Berti et al. (2013) method, including points without a linear or quadratic correlation (i.e., uncorrelated TS with stable trend) and points with PU errors. Therefore, two steps of calibration are carried out ((i): Uncorrelated_1 and PU Error_2; (ii): Uncorrelated_2 and PU Error_3) using the coefficient of determination ( � R 2 , Equation 4), slopes, and F statistic test to label these points accurately.
In Equation 4, ŷ i is the estimated valued by the regression model and � y indicates the average of all temporal samples in each PS. The two terms n À 1 ð Þ and (n À p À 1) indicate the degrees of freedom. Indeed, the R-square statistic is adjusted based on the residual degree of freedom, p, which refers to the number of variables in the regression models; thus, it is 1 and 2 for linear and quadratic regressions, respectively. As discussed in Section 3.2., the calculated F statistic of each TS from a stepwise regression is compared with the corresponding value in the F distribution table to evaluate the significance of quadratic regression (Christensen 2016). If the calculated F is greater than or equal to the table value, the TS is considered as Quadratic; otherwise, as Linear. Afterward, the following steps are implemented to detect (i) actual Quadratic and Linear trends, (ii) Stable TSs misclassified as unstable, and (iii) PUE TSs.
Firstly, PSs are categorized into two groups (pseudo-quadratic and pseudo-linear) using the comparison between F statistic calculated from a stepwise regression. Afterward, the pseudoquadratic TSs that weakly fit the quadratic regression level (Uncorrelated_1) -quadratic correlation determination smaller than 25% and very low values of the quadratic model constants-are excluded and classified as the Stable class (Figure 4(a)). Secondly, the points that agree simultaneously with two conditions -very high values of the F statistic and quadratic correlation determination smaller than 90% (PU Error_2)-are classified as PUE. This choice is justified by the fact that in the one-way ANOVA, a high F statistic value indicates large variations of the data inside the analyzed TS. Figure 4(b) shows a PUE TS classified wrongly as Quadratic. This example contains a vertical jump, which makes the observations deviate from their mean value. Finally, those points that fail both the Uncorrelated_1 and PU Error_2 tests are labeled as Quadratic (Figure 2(c)).
Furthermore, the classification of pseudo-linear TSs is realized similarly to the pseudo-quadratic ones. After excluding stable points (Figure 4(c)) using the Uncorrelated_2 test, the remaining TSs are categorized based on their slopes and linear correlation level ( � R 2 ). In the case that a TS has a high slope and low level of correlation, it is characterized by a vertical jump (PU Error_3). Figure 4(d) illustrates a TS with a high slope classified wrongly as Linear. Finally, the remaining points are classified as Linear after failing Uncorrelated_2 and PU Error_3 tests (Figure 2(b)).

Bilinearity
The remaining points after the above-mentioned steps are characterized by bilinear trends or discontinuities. In order to separate these three remaining classes, a two-sided Student's t discontinuity test with the upper-tail probability of 0.025 (95% confidence level) is first calculated to predict two intervals associated to the two linear segments preceding and following the breakpoint. In fact, the prediction intervals contain a range of values at the breakpoint for evaluating an overlap between two segments (Weisberg 2005). Therefore, if there is an overlap between the two prediction intervals, the TS is Bilinear with 95% probability; on the contrary if the TS has a vertical jump and goes to the next step. Furthermore, points are analyzed for probable Stable or PUE following the logic of the previous section.

Equal slopes test
Points in this level are discontinuous TSs at the breakpoints with vertical jumps. The equality of slope approach tests two or more linear regression lines to see if the lines have equal slopes. The partial F test demonstrates whether the slope (i.e. the linear velocity) of segments preceding and following the breakpoint is equal with a 99% level of significance. If the calculated F statistic is greater than the F percent point function, the hypothesis of equal slopes is failed, then, the TS is classified as DDV. In the alternative hypothesis, the TS is a DCV point (Equation 5) (Liu et al. 2007;Sachs 1982). for equality of slopes before and after a breakpoint (b).

DCV if F stat E � F
In which: where k indicates the number of regressions. Furthermore, subscripts symbols of L and R stand for segments before and after the breakpoint.

Periodicity Analysis
Time series have different variations, including secular trends, irregular variations, and cyclic fluctuations (Saad, Prokhorov, and Wunsch 1998). Secular trends (i.e. the general behavior of the data for a long period) have been discussed comprehensively in previous sections. On the other hand, cyclic fluctuations are defined by the period of their cycle, which may differ depending on the underlying phenomena over time (Puech et al. 2020;Schlögl, Widhalm, and Avian 2021). The spectral analysis assesses the frequency content of TSs. One of the basic spectral measures is the Power Spectral Density (PSD), which contains the spectral intensity of the input data for each analyzed frequency. From the analysis of the PSD, it is possible to study the underlying periodic phenomena by isolating them from noisy signals (Blackman and Tukey 1958;Cooley and Tukey 1965;Gardner 1988). The PSD estimate applies some sort of averaging or smoothing to reduce the effect of random variations. In this work, we employ the Weighted Overlapped Segment Averaging (WOSA) method (Welch, 1961, Welch 1967, also called Welch. This PSD estimate is obtained by first dividing the TS into overlapped segments, then averaging the periodogram values for all the segments. This technique has shown robustness against random variations. This periodogram averaging method enables us to find and evaluate cyclic fluctuations of TSs. It also reduces the variance of the periodogram using a windowed overlapped segmented averaging method in the frequency-domain operations (Solomon 1991). Finally, we implement Welch's PSD estimate on all TSs except those affected by PUE, which contain no periods. The detrending processing is first implemented on Linear, Quadratic, and Bilinear points. Since these three classes contain trends, this process can remove them to improve the accuracy of the periodical frequency. In this study, we estimate the PSD of TSs for seasonal periods. Considering the weather conditions in Spain, a season is considered as a interval of 170 to 190 days. TSs characterized by a dominant peak on each of the mentioned periods are selected as a periodic TS. The result is provided in the result section and Figure 9.

Results
The automatic PS TS classification procedure was applied over 1475 km 2 of Granada and Barcelona datasets with 2,077,082 PS. The aim is to classify accurately all TSs in seven classes. The threshold values utilized in the various classification stages were first validated using synthetic data to test their accuracies. A confusion matrix of 1250 random simulated targets in five classes was analyzed, followed by a validation of 10,000 real TSs, classified as PUE. Then, the TS classification maps of Granada region and Barcelona metropolitan area are provided to visualize the distribution of the land displacement trends generally. Overall accuracy (OA) estimates the overall effectiveness of the classification algorithm using the proportion of the total correctly classified true samples. The producer accuracy (PA) indicates how well the samples for each class were identified. It is computed by dividing the correctly classified class samples by the total reference points. Dividing correctly classified samples by the total number of classified points presents the user accuracy (UA), which illustrates the probability that a sample belongs in its true class. This process is shown with sample cases over two study areas and discussed. Finally, the outcome of the periodicity analysis using Welch's PSD estimate is provided.

Validation of the classifier
In this research, we proposed a synthetic evaluation to assess the proposed method's capability in the classification of ground displacement types. Since the method includes various statistical steps, a set of synthetic data was simulated with five types of classes: Stable, Linear, Quadratic, Bilinear, and PUE ( Figure A). 250 samples of every trend were randomly simulated, classified, and examined by a confusion matrix, including overall, producer, and user accuracies. It should be noted that noise was randomly added using the uniform pseudorandom number generator (Mersenne Twister (Matsumoto and Nishimura 1998)) within the 16-26 mm interval. The average and standard deviation of the jumps at 250 simulated PUE samples were respectively 20.9 mm and 3.2 mm. Figure 5 shows that the method classifies virtually the synthetic data with a 78% OA. It also indicates that the most accurate estimation was over the Stable TS, followed by Linear, Quadratic, and Bilinear. UA (sensitivity) and PA (specificity) represent what portions of classes that are correctly classified by the method. For instance, 89.28% shows that 21 and 10 data of Linear and Bilinear samples, respectively, that were wrongly classified as Stable. The diagonal also shows how many samples of each class were correctly predicted and the other values indicate the incorrect predictions of the method.
Moreover, PSI experts have evaluated 10,000 TSs of Granada and Barcelona datasets among the 83,540 PSs detected as PU Error_1 to _4, to check the validation accuracy. It was concluded that 91% of TSs classified in PU Error_1 was accurately detected, followed by PU Error_2 (73%) and Error_3 (68%) of linear and quadratic regressions. On the other hand, the Bilinear TSs characterized by PU error (PU Error_4) showed the lowest accuracy (62%).

Granada and barcelona PS TS classification
The method was implemented using the C++ language and run on a computer equipped with 32 Gb of RAM, Intel i7-8700 CPU with six 3.2 GHz cores on a Windows 10 64-bit operating system. The two datasets were processed separately. Over the Granada dataset, the processing approximately took 300 seconds for 458,249 PSs with 138 images. The processing of 1,618,293 PSs with 249 images of the Barcelona dataset took 18 minutes to classify all TSs. Figure 6 illustrates the final TS maps of the Granada region and Barcelona. The PS TS classified maps provide the spatial distribution of target trends. They can also highlight the anomalies of targets movements over each period using monthly, seasonal, and yearly maps. Around 70% of both study areas are dominated by Stable PS, that were depicted in green. In both datasets, the Stable class was followed by Bilinear; also, the lowest number of TSs were classified as discontinuous classes. Over Granada region and Barcelona metropolitan, the classes were: 16.2% and 15.7% Bilinear (in blue), 5.5% and 5.4% Quadratic (in red), 4.4% and 4.1% Linear (in yellow), 3% and 5.3% PUE (in white), and less than one percent for DDV and DCV classes. The computational load of classifying ground displacements trends is due to the breakpoints inside four classes (Bilinear, DCV, DDV, and PUE). This causes similar statistical characteristics among TSs. More times checking labels by additional sequential conditions can adjust the complexity. For instance, various metrics can evaluate the performance of models, such as Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), Akaike's Information Criteria (AIC), and Mallows Cp. However, it increases the processing time.
Three cases characterized by Bilinear, Quadratic, and Linear trends are shown in Figure 7. Figure 7 (a) includes Bilinear TSs, where is located over a -5000 m 2 building. More than 30 Bilinear points were detected, where three of them have been shown in Figure 7(a). Figure 7(b) also shows around 25 spots of ground deformations, including 20 Quadratic TSs in the white circle and three samples of Quadratic TSs. As an example of the Linear class (Figure 7(c)), among multiple Linear points distinguished over a 1250 m 2 building and a road, three TSs have been provided.  The study monitored the activity of three landslides (called Lorenzo-1, Rules Viaduct and El Arrecife) using PSIG products and in-situ observations, along with a detailed photo interpretation and a DEM. Accordingly, the PS TS classifier output was compared to the three land deformation cases. Figure 8(a) shows the Cortijo de Lorenzo landslide and three TSs of the PSs. It is a rotational landslide causing several cracks and piping. Over this landslide, the PS points were classified as Bilinear. Furthermore, multiple Quadratic, Linear, and Bilinear PS demonstrate a significant movement in Figure 8(b). Reyes-Carmona et al. (2020) investigated that a hazardous active landslide affected by active piping and opened, fresh cracks (Rules Viaduct landslide).
Finally, Figure  8(c) indicates a deformation surrounding and across a road with multiple Quadratic PS and several Linear trends. The El Arrecife landslide is a translational landslide characterized mostly by slopes. A large number of open cracks and pieces of pavements were found over this area.
As mentioned in section 3.7, a Welch PSD estimate was employed over TSs to indicate seasonal and cyclic fluctuations. Temperature variations, vegetation growth, humidity, and groundwater   changes are the natural causes of yearly fluctuations. In this study, since the PSI product of the Barcelona dataset was carried out for 5 years, we used 2.5-years as the longest cycle to assess a strong calculation. Considering the interval of datasets images, 15-days was the minimum period (2 � TS interval). However, we decided to choose six months to avoid noisy TSs and unreliable data in terms of the minimum number of images per month and maximum period. Figure 9 shows the number of TSs extracted from four classes in 6-month period. While the estimator was applied over more periods (see section 3.7), only periods shown in the Figure 9 were considered because most of the TSs indicated various peaks. On the other hand, only less than 1% of the PSs included periodic fluctuations; however, it can reveal critical movements regarding seasonal or yearly human activities, such as water pumping or injection. Although the estimate could also be calculated for higher or lower periods, the reliability of the output is limited due to the following reasons. First, regarding the shorter periods (e.g. one month), there should be more data over different days. Second, longer periods would provide better performances in terms of yearly periods.

Discussion
Time series classification over large areas improves the understanding of the spatial distribution of the deformations and improves the overall qualitive interpretation of DInSAR products. However, unknown and identified sources of noise and errors can undoubtedly affect the reliability of the output. These error sources are affecting the PS TS products more seriously than the deformation velocity product. In this study, we modified the method illustrated in (Berti et al. 2013) and proposed a workflow to investigate the misclassified TS and a critical noise source (i.e. vertical jumps caused by PUE), as well as significantly reduce the classification commission errors. Following are the discussed features and advantages represented by this research: • The spatial distribution of ground movement trends was shown using seven classes. The classified maps could reasonably demonstrate the spatial propagation of moving and non-moving targets in the study areas.
• Four stages of PUE detection were firstly proposed to classify TSs affected by these errors.
The tests on simulated and real data showed that 84.4% and 73.5% of them were labeled correctly. • A dataset of five movements trends was simulated to validate the method. Although the method was not strong to identify bilinear data as other trends, it could correctly classify 78% of all simulated data. Furthermore, 10,000 of classified PUE TSs were validated by experts, showing the proposed workflow classified approximately with an accuracy of 73.5%. • We also supplied four tests to obtain and separate stable TSs classified incorrectly as unstable. Uncorrelated_1 to 3 tests were applied on the output of the original classification method (Berti et al. 2013), including Linear, Quadratic and Bilinear TSs to discover actual stable points. Around 40% of misclassified TSs were accurately detected and labeled as Stable. • Various anomalies were recognized in the TSs, which were mostly showing similar behavior to Bilinear and PUE trends. Through bilinear anomalies, many significant sudden movements may be investigated. Anomalies detection can develop procedures to find TSs with interesting trends and variations from the noise and temporal analysis perspectives. Figure 10 compares the proposed classification with the one of Berti et al. (2013). Figures 10(d) and (b) indicate that using the method of Berti et al. (2013) Besides, TSs affected by 19 mm or higher PU errors were initially investigated at the beginning of the process. In total, 83,540 PUE from both datasets were classified in the first stage. 16,160 TSs were also detected in the following stages: Linear, Quadratic, and Bilinear. Several PUE TSs were incorrectly in the Bilinear class due to behavior similar to the PUE and Bilinear trends. Although less than 6% of TS were PUE, it can trigger the automated step in detecting PUE PSs in the final DInSAR products using statistical methods to refine PSI products from PU errors and increase the reliability of interpretations.
Furthermore, the relatively high percentage of Bilinear class illustrates another possible variation, including PU errors with different fluctuations unrecognized from four PUE detection stages (Figure 11 (a)), noisy PSs (Figure 11(b)), and TSs with rapid jumps at final images (Figure 11(c)). Mostly, these groups are not indicating significant displacements, such as Figure 11(b) and c TSs, which demonstrate Stable targets. On the other hand, the Bilinear class consists of TSs with significant vertical jumps  showing notable anomalies. For instance, a high slope in the second segment of the Bilinear trend states a significant abrupt displacement (Figure 11 (d)). Since the breakpoint time can be identified using segmented regression (i.e. the start time of the slope's increasing or decreasing), various helpful information can be extracted, including the starting time of a significant displacement (e.g. October of 2019 in Figure 11(d)).
The classification approaches based on statistical analysis do not generally involve spatial information of the neighbor PSs. A group of unstable PSs over a region highly demonstrates the existence of ground deformation, while scattered PSs or an isolated point may draw little attention. Serving the neighborhood characteristics can attach input data to train a supervised method, as well as the land movements inventory, DEM, and land cover maps. Nevertheless, some TSs consist of more than one jump, which makes the interpretation complicated. Indeed, the date of the higher jump is diagnosed as the breakpoint. Although they are correctly classified as Bilinear, the user may be interested in other breakpoints. Finally, when classifying over a large area and a huge number of satellite images, several limitations should be considered with respect to the reliability, the level of accuracy, ancillary data, and the complexity of the classification method. The classifier labeled each PS in 6 � 10 À 4 seconds since more than 60% of points were analyzed more than two times, which increases the processing time. The processing time still needs an improvement to consider the algorithm for carrying out over large datasets. Despite these limitations, this paper has highlighted the value of classifying PS TSs over large datasets using big satellite data.

Conclusion
A PS TS classifier has been implemented to evaluate the distribution of deformations trends, detect significant movements, and test the reliability of the proposed method over large areas. Various statistical analyses have been applied to classify PS TSs into seven stable and unstable trends. Using simulated data, the method accurately classified 77.7% of the random synthetic data, whereas stable points were completely classified. Furthermore, four stages of analyses detected points affected by PU errors for the first time. 10,000 TSs have also been validated qualitatively by experts to evaluate the reliability of the PU errors detection.
Utilizing statistical tests and data processing of wide areas afford accurate and reliable trends classification. 2,077,082 PS points over 1475 km 2 of high slopes and urban areas have been analyzed to achieve an accurate classification. The produced maps can add complementary information of the ground state to available datasets over Barcelona and Granada.