Automated classification of Persistent Scatterers Interferometry time series

We present a new method for the automatic classification of Persistent Scatters Interferometry (PSI) time series based on a conditional sequence of statistical tests. Time series are classified into distinctive predefined target trends, such as uncorrelated, linear, quadratic, bilinear and discontinuous, that describe different styles of ground deformation. Our automatic analysis overcomes limits related to the visual classification of PSI time series, which cannot be carried out systematically for large datasets. The method has been tested with reference to landslides using PSI datasets covering the northern Apennines of Italy. The clear distinction between the relative frequency of uncorrelated, linear and non-linear time series with respect to mean velocity distribution suggests that different target trends are related to different physical processes that are likely to control slope movements. The spatial distribution of classified time series is also consistent with respect the known distribution of flat areas, slopes and landslides in the tests area. Classified time series enhances the radar interpretation of slope movements at the site scale, pointing out significant advantages in comparison with the conventional analysis based solely on the mean velocity. The test application also warns against potentially misleading classification outputs in case of datasets affected by systematic errors. Although the method was developed and tested to investigate landslides, it should be also useful for the analysis of other ground deformation processes such as subsidence, swelling/shrinkage of soils, or uplifts due to deep injections in reservoirs.

deformation. Our automatic analysis overcomes limits related to the visual classification of PSI time series, which cannot be carried out systematically for large datasets. The method has been tested with reference to landslides using PSI datasets covering the northern Apennines of Italy. The clear distinction between the relative frequency of uncorrelated, linear and non-linear time series with respect to mean velocity distribu- 10 tion suggests that different target trends are related to different physical processes that are likely to control slope movements. The spatial distribution of classified time series is also consistent with respect the known distribution of flat areas, slopes and landslides in the tests area. Classified time series enhances the radar interpretation of slope movements at the site scale, pointing out significant advantages in comparison with

Introduction
The detection of ground displacements by space borne synthetic aperture radar interferometry has progressed, over the last two decades, from the use of single interferograms, as in the pioneer work of Massonnet and Feigl (1998), to the use of advanced 25 persistent scatterers multi-interferometry (PSI) techniques, among which Permanent 208 sion model is generally fitted to the data and the average displacement rate is used to describe the entire time series. On the other hand, several studies have shown that even the analysis of few relevant time series may provide useful information on slope dynamics (Meisina et al., 2008;Cigna et al., 2011). Cigna et al. (2011), in particular, manually classify the time series 20 of several radar targets to identify the change in deformation rate caused by a specific slope instability event of known date of occurrence. The visual analysis and supervised manual classification of time series, however, is a tedious, time-consuming and subjective work which cannot be carried out for large datasets. Hence, criteria and methods for an objective, automated classification of time series included in large datasets are pose a PSI time series clustering approach is that of Milone and Scepi (2011). They used the partition-based clustering algorithms CLARA (Clustering for Large Applications) which defines "k" clusters in an entire dataset on the basis of the identification of "k" representative objects in a sub-dataset. The type and the number of clusters are therefore specifically dependant on the dataset itself; if a different dataset is pro- 10 cessed, different clusters are possibly generated. This is a relevant problem in radar interpretation since PSI information from different orbits, tracks and satellites datasets have usually to be used and compared. Our approach aims to overcome such a limitation, by allowing PSI time series to be clustered into fixed predefined target trends (uncorrelated, linear, quadratic, bilinear, discontinuous) which can be potentially recog- 15 nized in any PSI dataset and that are believed to be interpretable in terms of physical processes related to slope instability. Inevitably, time series classification is affected by the quality of the dataset. In this paper, the proposed automated classification algorithm (hereafter referred to as PS-Time) has been tested with reference to landslides in the northern Apennines of Italy by using ENVISAT datasets available in the frame 20 of EPRS-E project (Extraordinary Plan of Remote Sensing of the Environment) of the Italian Ministry of Environment. Datasets of the EPRS-E have been generated using standard PSI processing and are notoriously affected by a high noise-to-signal ratio. Nevertheless, results obtained with PS-Time show that even with such rough quality datasets, the method can provide results that are generally consistent with the nature 25 of the analysed landslide phenomena and that, as such, can potentially improve radar interpretation of mass movements at regional to site-specific scale.

Identification of distinctive target trends in PSI time series
As the main goal of the research is to develop an automatic procedure that classifies the PSI time series according to their peculiar trends, the a-priori identification of trends that might indicate changes in time of physical processes related to slope instability (or 5 other terrain movements) is crucial in order to tailor the selection of suitable statistical techniques.
To this purpose we randomly selected 1000 time series from our sample dataset (see Sect. 3) and performed a visual analysis to find distinct patterns of displacement. Six recurrent patterns were observed and consequently identified as target trends in PSI 10 time series (Fig. 1 characterize ground deformation processes acting over long time scales (creep, natural subsidence, steady motion of dormant landslides). Non-linear trends (type 2 to 5) indicate a change of the displacement rate during the monitoring period. Regardless of whether this change is progressive (type 2), abrupt (type 3), or discontinuous (type 4-5), non-linear trends catch our attention because indicate a variation of the displace-5 ment field.

Statistical tests for automated classification of PSI time series
An automatic procedure based on the sequential application of a number of statistical tests ( Fig. 2) was developed to classify each time series into one of the six target trends shown in Fig. 1. Each testing method is hereafter briefly described.

(A) Linear regression
The first step in the analysis (test A in Fig. 2) is an ANOVA F-test for the significance of the linear regression (Davies, 1986). The PS displacements are plotted against time and fitted by a linear regression model. The F-statistic is then used to compute the probability value p 1 that the regression coefficient β 1 (the slope of the regression line) is 15 equal to zero: if p 1 is less than the selected level of significance α 1 , the null hypothesis of no-correlation (H 0 : β 1 = 0) is rejected and we accept that there is a significant linear relationship between the two variables; conversely, if p 1 > α 1 the null hypothesis cannot be rejected and we conclude that the two variables are not linearly correlated. Obviously, the absence of a linear correlation does not imply that the displacements 20 are randomly distributed with time because the time series could be described by a higher order model (Draper and Smith, 1981;Sen and Srivastava, 1990). Based on our experience, however, a time series with an average slope close to zero (β 1 ≈ 0) typically indicates a stable PS with no appreciable movements. Therefore, if the linearity test fails (p 1 > α 1 ) the time series is directly classified as "uncorrelated" (type = 0; Fig. 2

(B) Segmented regression
If linear regression is significant (p 1 ≤ α 1 ) the time series is tested against a bilinear model (test B, Fig. 2). The segmented regression, also known as piecewise or changepoint regression, is a regression technique in which the independent variable is divided into intervals and a separate line segment is to fit each interval (Main et al., 1999;5 Steven, 2001). The segments are introduced to see if there is an abrupt change in slope in the data and to determine where the change occurs (breakpoint). For this analysis we follow the method proposed by Main et al. (1999): i. the time series t 1 , ...., t n is divided into two parts separated by a breakpoint t b , which is moved along the series from b = 5 to b = n − 5 (we assume that a minimum of 5 data points is required to define a segment, in order to avoid very short segments at the beginning or at the end of the time series) ii. for each breakpoint, a two-line unconstrained model is fitted for the segments t 1 , ...., t b and t b+1 , ...., t n and the goodness of fit is evaluated using the Bayesian Information Criterion BIC (Main et al., 1999): where RSS is the residual sum of squares and k is the number of parameters in the model (k = 3 for a two-line regression); iii. an overall linear and quadratic fit is computed for the entire series t 1 , ...., t n , and the corresponding values BIC L and BIC Q are computed by Eq. (1) using k = 1 for 20 the linear model and k = 2 for a quadratic model.
The BIC criterion (Schwarz, 1978)  way the BIC resolves the problem of overfitting and identifies the best model as the one that strike a balance between fitting the data well (low RSS) using only a few parameters (low k). As an example, Fig. 3 shows the values of BIC(t b ), BIC L and BIC Q computed for the time series shown in Fig. 1d. By comparing the BIC values we infer cases where a double-slope assumption is 5 statistically better than a linear or quadratic model. Specifically, the time series justifies a breakpoint if the minimum value of BIC for the segmented model (BIC min ) is lower than both BIC L and BIC Q : if BIC min < (BIC L and BIC Q ) ⇒ significant breakpoint exists in the time series (2) In the example of Fig. 3 the segmented regression outperforms both the quadratic 10 and the linear fit. There is then statistical evidence of a breakpoint at the beginning of 2006 (compare Fig. 3 with Fig. 1d).
Although condition (2) suffices to establish the existence of a breakpoint, we used a more restrictive criterion that also allows the user to calibrate the test outcome. The breakpoint is considered significant when the so called "evidence ratio" B w of the break-15 point exceeds a predefined threshold B th (with B th ≥ 1), that is: The evidence ratio B w is computed as proposed by Wagenmakers and Farrell (2004): where w 1,2,3 indicate the weights of the segmented, linear and quadratic model respectively. The weights are obtained by the normalization of the relative likelihoods: where the subscript i = 1, 2, 3 indicates the regression model and ∆ i is the difference of BIC between the bilinear model and the other two (∆ 1 = 0; ∆ 2 = BIC L − BIC min ; ∆ 3 = BIC Q − BIC min ). In the example of Fig. 3 the evidence ratio is B w =1.06. Whether or not this value is significant depends on the selected threshold B th , thus the higher is the value of B th the more restrictive is the test. Note that the limit value B w =1 indicates the 5 same fit for the bilinear and the other models (BIC min = BIC Q or BIC min = BIC L ).

(C) Quadratic regression
If the time series is linearly correlated (test A passed) and does not have a breakpoint (test B failed) a quadratic fit is performed to test for the significance of quadratic over linear fit (Davis, 1986). An ANOVA F-test evaluates the probability p 12 that the 10 quadric term is not contributing to the regression. If p 12 ≤ α 12 (the selected level of significance) the quadratic term is making a significant contribution to the regression and should be retained. The time series is therefore classified as "quadratic" (type=2). On the contrary, if p 12 > α 12 the additional term is not contributing significantly to the regression and the time series is classified as "linear" (type = 1).

(D) Discontinuity test
The time series characterized by a significant breakpoint (test B passed) are tested to evaluate if there is a vertical jump in the data. The discontinuity test is based on a simple comparison of the prediction intervals at the 95 % level of significance computed for the two linear segments before and after the breakpoint (Weisberg, 1985), as shown 20 in the example of Fig. 1e. The prediction intervals provide the range of expected values at the breakpoint: if the two intervals overlap each other, the time series is said to be continuous and it is classified as "bilinear" (type=3); conversely, if the two intervals do not overlap, the series is said to be discontinuous at the breakpoint ( Fig. 1e) and the two segments are tested for the equality of slopes.

(E) Equality of slopes
The F-test for the equality of slopes (Quinn and Keough, 2002) is applied to a discontinuous time series (test D passed) to evaluate if there is a significant difference of velocity before and after the breakpoint. The null hypothesis is that the two segments have the same slope (H 0 : V 1 = V 2 ). If the computed probability p V is greater than the 5 selected level of significance (set to 0.05) the null hypothesis cannot be rejected and we conclude that the difference in velocity is not significant (V 1 ≈ V 2 ). The time series is then classified as "discontinuous with constant velocity" (type = 4). If p V ≤ 0.05 the null hypothesis is rejected (V 1 = V 2 ) and the time series is classified as "discontinuous with different velocity" (type = 5). It must be noted that the test for equality of slopes is adversely affected by data scattering and may have not enough power to detect differences that do in fact exist. The result of the test must be then carefully evaluated.

Calibration of thresholds in statistical tests
In hypothesis testing, the significance level α is the criterion used for rejecting the null hypothesis H 0 (Davis, 1986). The value of α is always small because it represents the 15 probability of a Type I error, that is the probability of rejecting H 0 when it is actually true. A typical value of α is 0.05, but the choice of any cutoff significance level is arbitrary. For instance, a stiffer standard α = 0.01 can be adopted if a stronger evidence is needed to reject the null hypothesis. PSI time series have a limited number of data points and are quite noisy. Therefore, 20 statistical tests need to be stringent in identifying non-random data and the significance levels must be properly calibrated. A way to perform this calibration is to run the analysis using different values of the significance levels, and to compare the results with that obtained by an expert visual classification. The suitable significance levels are those providing the best correlation with the expert judgment. This technique has been used Introduction

Annual periodicity of time series
Several time series are characterized by cyclic fluctuations of the displacement values over the monitoring period. The amplitude of these fluctuations appears to be quite variable in time, but the period is relatively constant, with a typical periodicity of about one year. Such a cyclic behavior is usually observed in slow-moving PS and tends to 5 disappear when the velocity increase. In most cases, periodic time series are classified as "uncorrelated" (type = 0) because the linear velocity is close to zero. The detection of periodic time series can be of practical interest. For instance, cyclic movements of the ground surface may indicate significant shrink/swell phenomena in expansive soils, bradyseism, or man's activity such as injection of gas into a reser- A convenient way to evaluate the periodicity of a time series is to compute its power spectral density P (Priestley, 1981). The power spectrum is generated by the basic 15 FFT analysis and provides a representation of the magnitude of the various frequency components (f ) of the series. Figure 4 shows the power spectrum P (f ) computed for two sample PSI time series, one of which clearly characterized by periodic fluctuations (Fig. 4c). The x-axis covers the frequencies from zero to the Nyquist frequency (half the sampling rate) while the y-axis indicates the spectral power. As expected, the 20 periodic series is characterized by a well-defined peak at f ≈1 yr −1 (Fig. 4d) while the random series does not show any peak in that frequency range (Fig. 4b). Both the series, however, show a peak in the low-frequency band at f < 0.5 which correspond to the fundamental frequency f 0 (that is the frequency associated to the longer waveform fitting the data); when the series follows an half-sine trend (as in Fig. 4a) f 0 is where P 0 and P 1 are the spectral peaks in the two frequency bands (Fig. 4b-d). The index AP ranges from 0 (no annual periodicity) to 1 (very strong annual periodicity). In the example of Fig. 4, AP is 0.02 for the random series and 0.91 for the periodic series.

5
Periodicities other than annual are not captured by AP.

Graphical user interface
The classification method is integrated into a Graphical User Interface called PS-Time (Fig. 5). PS-Time is written in Matlab 7.11 and it is freely available as a standalone application that can run on Windows and Unix systems without requiring Matlab or its 10 auxiliary toolboxes (http://www.bigea.unibo.it/it/ricerca/pstime). The input file is an excel spreadsheet containing the PS data (a sample input file can be found in the web page). PS-Time plots the time series of the selected PS and visualizes the result of the statistical analyses (Fig. 5). Beside the classification of the time series and the analysis of annual periodicity, the program analyzes the 15 residuals of the linear model and computes several descriptive statistical indexes such as the coefficient of determination r 2 , the residual mean square error RMSE, and the Standard Deviation of Slope (an index of the scatter of the time series based on the algorithm proposed by Frankel and Dolan (2007) to measure the roughness of natural surfaces).

20
The button "Analyze all" performs a sequential analysis of all the dataset. The output is an Excel spreadsheet ( Table 1) that is ready to be joined to the PS shapefile by the Code field. All the relevant information obtained from the statistical analyses can then viewed into a geo-referenced GIS software to support the interpretation of PSI data.

Test area and sample dataset
PS-Time has been calibrated and tested using a PSI dataset covering about 2000 km 2 of the northern Apennines of Italy (Fig. 6) visat data and to the relatively high coherence threshold adopted for regional scale PSP processing, the dataset contains "only" 63 707 persistent scatterers, which generally correspond to buildings. Average density is 35 PS km −2 , but it varies significantly NHESSD

Calibration of the statistical thresholds
Before applying the automatic classification procedure, suitable levels of significance must be set (see Sect. 2.2). To this purpose three of the authors performed an expertbased classification of 1000 PS randomly selected within the study area. The results were then compared and discussed to reach a shared classification. In general we 10 found good agreement between the subjective assessment of uncorrelated (type = 0) and linear (type = 1) time series, while it was more difficult to agree on the quadratic (type = 2), bilinear (type = 3) and discontinuous (type = 4 and 5) trends. We then decided to group all the "non-linear" time series (types 2 to 5) into a single class and to restrict the calibration analysis on the statistical thresholds α 1 , α 12 , and B th . These 15 are the significance levels used to discriminate between uncorrelated, linear, and nonlinear trends in our automatic procedure (tests A, B, C; Fig. 2).
The optimal values of α 1 , α 12 , and B th were selected using the ROC curve method (Green and Swets, 1966). The three statistical thresholds were varied in a wide range (α 1 , α 12 = 1·10 −5 ÷0.4, 57 log-spaced values; B th = 1.0÷1.5, 11 equally-spaced 20 values) and, for each combination, the performance of the model was evaluated by comparing the automatic with the expert classification. The results of the analysis are shown Fig. 7. Each point represents the capability of the model to predict a certain type of time series (for instance linear) for a given combination of the statistical thresholds. The more distant is a point from the 1:1 line (which 25 indicates a random prediction) the higher is the predictive capability of the model. Therefore, the best combination of α 1 , α 12 , and B th is that providing the more distant NHESSD points for all the three time series types. A good combination seems to be α 1 = 0.01, α 12 = 0.01, and B th = 1.0, that provides a prediction accuracy of 84%, 82%, and 90 % for the uncorrelated, linear, and non-linear time series respectively (see the black dots in Fig. 7). As expected (see Sect. 2.2) our calibrated significance levels are lower than those usually adopted (0.01 instead of the conventional 0.05) because data are noisy 5 and we therefore need a strong evidence (that is a lower level of significance) the null hypothesis.

Automated time series classification
Results of PS-Time can be analyzed without considering the spatial location of PS, in order to provide a general analysis of the time series in the dataset, or they can be 10 displayed on a map to investigate the spatial distribution of target trends with respect to topographic and geologic attributes (such as the known distribution of flat areas, slopes and landslides), or to perform site-specific analyses.

Relative frequency and velocity distribution of classified time series
The pie chart in Fig. 8a shows the percentage distribution of the target trends in the 15 study area. The majority of the time series is classified as "uncorrelated" (59 %) or "linear" (33 %). Non-linear trends (type 2 to 5) account for the 8 % only and mostly consist of bilinear functions. It seems reasonable that in a large area most of the PS are stable or characterized by slow movements without significant variation of velocity, and that non-linear displacements are a small minority. The frequency distributions of the 20 mean displacement rate of the different target trends further support the correctness of the analysis (Fig. 8b). The velocity distribution of uncorrelated time series is centered on zero and nearly all the PS have a velocity between +2 and −2 mm yr −1 . Interestingly, this is the velocity threshold typically used to cutoff significant ground movements in practical applications (e.g. Catani et al., 2012 around +1.5 mm yr −1 . The distribution is skewed to the left and it is clearly separated from that of the uncorrelated PS. Non-linear time series are even more separated and remarkably skewed toward higher displacement rates (both positive or negative). Such a clear distinction between the three frequency distributions suggests that different physical processes are controlling the displacements of uncorrelated, linear, and non-5 linear PS. The latter, in particular, should be of particular interest since they highlight a radical change of behavior during the monitoring period.

Spatial clustering
A simple way to establish the significance of the method is to see whether the classified PS are clustered or randomly distributed over space. Spatial clustering of the target 10 trends would suggest the existence of a common process influencing the temporal behavior of a group of PS (such as a landslide or local subsidence), while a random distribution would indicate that time series classification not add any significant spatial information.
Most of the methods for cluster analysis compare the pattern of the data to that 15 produced by a homogeneous Poisson process as shown in Haining (2003). In the case of PS, however, data are inherently clustered since the reflecting elements (building, pylons, outcrops) are usually grouped in space. Therefore, to evaluate if classified PS are spatially aggregated, we have to see if they are more clustered than the parent distribution.

20
The analysis was performed using the univariate Ripley's K function (Dixon, 2002) and results are shown in Fig. 9. Each curve represents the difference between the "variance stabilized" Ripley's K function (L) and the value expected for a complete spatial randomness distribution (CSR). The difference L − CSR is plotted against the distance d , which is a scale parameter indicating the size of the cluster: the more the curve 25 is away from the x-axis (random distribution) the more clustered are the data. As can be seen, all the classified PS (uncorrelated, linear, non-linear curves) exhibit a strong pattern of spatial clustering and their degree of aggregation is higher than that computed for the parent unclassified PS (overall data curve). Moreover, classified PS fall outside the 99 % confidence intervals of a random distribution (computed by performing 1000 Montecarlo simulations in which the target trends were randomly distributed within the overall clustered domain) thus indicating that the observed differences are 5 statistically significant. On the other hand, spatial clustering can be clearly recognized on the map. Uncorrelated PS (type 0) are generally concentrated in flat areas such as fluvial terraces and valley bottoms, and along stable watershed divides; linear PS (type 1) are mainly located on slopes (both inside or outside mapped landslides) or near the edge of scarps 10 or steep slopes; non-linear PS (types 2 to 5) typically fall inside landslide deposits or in the surrounding areas. The histograms in Fig. 10 illustrate this general tendency showing that the relative percentage of linear and non-linear PS (with respect uncorrelated PS) increases from flat areas to slopes to landslides. Of course, the spatial distribution of classified PS is more complex when analyzed at the local scale. Ancient dormant 15 landslides, for instance, may show uncorrelated time series (type 0) that indicate the absence of detectable movements, while non-linear PS can be found in flat areas or outside the landslides as a consequence of deformation processes not related to slope instability such as ground subsidence, soil swelling, or structural deformations of buildings. Examples site-specific distributions are discussed in the next section.

Enhanced radar interpretation of slope movements
The main strength of the method is to enhance the interpretation of PS data. Different temporal trends likely correspond to different deformation processes, and their spatial distribution can greatly improve our understanding. A first example is given in Fig. 11. the south of Modena. The village is surrounded by several landslide bodies, principally dormant earthflows, and threatened by rockfalls from the steep slopes to the south. The vast majority of the 389 PS identified in the area show a very low mean velocity (from −3.0 to 3.0 mm yr −1 ; Fig. 11a) which lead to consider the area as uniformly stable. Time series analysis, however, provides a more accurate picture of the stability conditions.

5
Although most of the PS are characterized by uncorrelated time series (type 0, stable points) several significant clusters of linear (type 1) and non-linear (types 2 to 5) PS appear in the northern sector of the village (Fig. 11b). More attention must be devoted to these areas in spite of the similarity of the mean displacement rate. Another illustrative example is given in Fig. 12. The map shows the confluence between the Reno and the Silla river, about 40 km to the south of Bologna. All the area is deeply affected by earthflows of variable extent and degree of activity. Nearly all the landslides reach the valley bottom, and in recent years several unsafe slopes have been occupied by urban development due to the lack of flat areas. The general distribution of classified PS follows the one previously described. Most of the uncorrelated 15 (stable) PS are located on fluvial terraces or on gently sloping foot slopes, while linear and non-linear PS typically fall inside the landslides. In particular, a large cluster of non-linear PS characterizes the toe of the big earthflow to the east. The extent of the cluster and the fact that most breakpoints occur in the same period (end of 2006) reveal a potentially critical condition. Results are more difficult to interpret when PS 20 with different trends are mixed together without any clear clustering, as in the case of the landslide located in the upper-left corner of Fig. 12. In these cases a close comparison of the time series and a detailed geomorphological analysis must be performed to judge the significance of the detected differences. Time series analysis may provide useful information even when a critical areal is 25 already clearly visible by the PS velocity alone. The Borgo village, for example, is located about 50 km to the south of Bologna and it is threatened by two dormant earthflows along the southwestern side (Fig. 13). PS velocities clearly indicate on-going deformations behind the landslide scarps (Fig. 13a) the useful information that the displacement rate is everywhere constant (linear trends have constant velocity) or increasing (bilinear trends with V 2 > V 1 ) (Fig. 13b). The lack of decelerating PS and the fact that all the breakpoints fall in the same period (first semester of 2008) are a strong indication of the critical site conditions. 5 The reliability of the results highly depend on the quality of the dataset. Because of the indirect nature of the radar measurement, and of the many sources of errors involved in the data generation, PSI time series typically show a quite high noise-to-signal ratio. While these errors have a limited impact on the overall linear velocity, they may change the classification outcome.

Influence of dataset's quality
In general, random errors add spatial noise to the results but do not invalidate the analysis. All the statistical tests, in fact, compare a pre-defined regression model (linear, bilinear, quadratic) with a random data distribution and are therefore quite robust in this respect. On the other hand, an erratic drift of several consecutive values at the beginning or at the end of the time series might lead to classify a time series as "bi- 15 linear" while it is actually "linear" or "uncorrelated". This problem is quite common in our datasets, and is particularly relevant for slow-moving PS, because random errors have a detrimental effect when the displacement values fluctuate around zero. For this reason a specific function to exclude initial and final data from the time series has been implemented in PS-Time.

20
A more subtle problem occurs when PSI datasets are affected by a systematic error. For instance, a bias in the displacement values will cumulate over time and will not only give an erroneous value of the average velocity, but it will also add a drift to the time series that might lead to a misleading classification.
We have experienced this problem while testing PS-Time using another dataset of  (Fig. 14a, solid line). If this dataset is processed without any correction, the spatial distribution of target trends with respect to flat areas, slopes and mapped landslides is not as expected, since we see no difference in the relative abundance of non-linear time series between the different morphological 5 units (Fig. 14b). However, when the dataset is corrected by applying a velocity offset of −1.15 mm yr −1 (in order to adjust the frequency peak to zero; see the dashed line in Fig. 14a) the analysis provides the expected results and a clear clustering of non-linear trends appears within mapped landslides (Fig. 14b). The applied correction changes the classification of 4932 time series over 8737 (56 %), indicating that even a small bias of few mm yr −1 may have a deep impact on the results. To further highlight the problem let consider the example in Fig. 15. The maps show the slope of Montecreto, which is classified by existing landslide inventory maps as a potentially unstable slope. Time series classification based on the original dataset shows 24 time series as bilinear (type = 3) and only 3 as uncorrelated 15 (type = 0) (Fig. 15a). Moreover, all the bilinear PS are stable until first semester of 2007 and then start moving at average velocities of 6 mm yr −1 (see the example in Fig. 15c).
This might induce to speculate that slope processes are progressing quite rapidly at present and that some cause of reactivation affected the entire slope in 2007. However, time series classification based on the unbiased dataset provides a different, and 20 possibly more realistic, picture. Among the 27 PS only 1 is now classified as "bilinear" and 26 become "uncorrelated" (Fig. 15b). The mean velocity is close to zero for almost all PS, which is consistent to the fact that most of the buildings show no damages and there are no evidences of accelerated slope movements. Still a break in the time series is in 2007 (Fig. 15d)

Conclusions
The new approach for the automatic classification of the Permanent Scatterers time series allows to enhance the radar interpretation of slope movements, providing a better understanding of the deformation phenomena than that based on the conventional analysis of mean velocities alone. The statistical algorithms adopted in the analysis 5 use standard statistical technique to determine whether a measured time series differs significantly from six pre-defined trends recognized by expert visual inspection of 1000 time series. The identified target trends (uncorrelated, linear, quadratic, bilinear, discontinuous with constant velocity, discontinuous with variable velocity) describe different styles of ground deformation and likely indicates the response to different de-10 formation processes. The application in a study area of the northern Apennines of Italy shows the potential of the method. PS with similar trends are spatially clustered and fairly distributed within the main morphological units, with most of the uncorrelated time series (stable PS) located on the fluvial terraces and in the flat areas, linear PS on slopes, and non-linear 15 PS on landslides. At the local scale, the time series analysis has proven to be an effective integration of the conventional velocity analysis allowing to distinguish between areas characterised by different temporal evolution (constant velocity, acceleration, deceleration, sudden changes of the deformation rate). The test application also shows that the reliability of the results highly depend on the quality of the dataset, and that 20 even a small bias in the data may provide misleading classification of the time series. Therefore, it is important to correct a biased dataset for any systematic error.
Although the method was developed and tested to investigate slope stability, it should also be useful for the analysis of other ground deformation processes such as subsidence, swelling/shrinkage of soils, uplift due to deep injections. The ever increasing 25 availability of space borne radar data and their increasing quality will make the time series analysis even more powerful. STDS Standard Deviation of Slope (mm yr −1 ) (Frankel and Dolan, 2007) AP Index of Annual Periodicity computed by Eq. (6) P1 F-statistic for the significance of the linear regression (to be compared with the selected level of significance α 1 ) P2 F-statistic for the significance of the quadratic regression (not used in the classification procedure) P12 F-statistic for the significance of a quadratic term added to a linear regression (to be compared with the selected level of significance α 12 )

BL
Result of the bilinear regression analysis: 0 = the bilinear model does not provide a better fit to the data than the linear and quadratic models; 1 = the bilinear model provides a better fit BICW Evidence ratio of the breakpoint B w computed by Eq. (4) Type Time series classification: 0=uncorrelated; 1 = linear; 2 = quadratic; 3 = bilinear; 4 = discontinuous with constant velocity; 5 = discontinuous with variable velocity.

V1
Mean velocity before the breakpoint (V 1 , mm yr −1 ) computed by bilinear regression for non-linear time series (types 2 to 5)

V2
Mean velocity after the breakpoint (V 2 , mm yr (1) as a function of the breakpoint position t b . BIC L and BIC Q indicate the linear and quadratic fit computed for the entire time series. The inset shows the "evidence ratio" of the breakpoint computed by Eqs. (4) and (5).