A Novel Method of Generating Deformation Time-Series Using Interferometric Synthetic Aperture Radar and Its Application in Mexico City

As a result of rapid societal development and urbanization, the pumping of groundwater has gradually increased. Land subsidence has thus become a common geological disaster, which can result in huge economic losses. Interferometric synthetic aperture radar (InSAR), with its large-scale and high-accuracy monitoring characteristics, can attain information on Earth surface deformation using the interferometric phase between couples of SAR images acquired at different times. Time-series results for the ground surface are the key information required to understand the deformation pattern and further study the reason for the subsidence. However, in recent research, most methods for resolving time-series deformation—like the Berardino method—that use residuals in functional model solving and distinguish high-pass displacement and the atmospheric component by filtering do not generally work well and functional models focusing on prior information in the time-series solution process are not always available. In this paper, to solve the above problems, 34 Sentinel-1A descending mode scenes of Mexico City captured between 2015/04/13 and 2016/09/10 are used as experimental data. Firstly, a new functional model is provided to obtain the deformation time-series. The nonlinear deformation and atmospheric phase are combined as an unknown parameter and the method of singular value decomposition (SVD) is used to solve this variable. The nonlinear displacement and atmospheric phase are then separated by the singular spectrum analysis (SSA) method. Finally, the total land subsidence time-series is obtained by adding together the linear displacement and nonlinear displacement. Two typical methods and the proposed method were compared using both unit weights and adaptive weights. The experimental results show that the proposed method can obtain a more accurate time-series deformation result. Moreover, the different weights do not result in significant differences and the solved atmospheric and nonlinear phases have good consistency with the interferogram phase.


Introduction
As one of the most important disaster types in regions of dense population, land subsidence is common and can be induced by either natural or human activity.Land subsidence has the characteristics of a long duration, a large spatial range and a slow rate of deformation [1][2][3].Ground deformation can not only result in economic losses but also poses a great threat to the safety and survival of humans [4].As a result, continuous monitoring of surface deformation is necessary to understand the changing process and minimize damage.
Compared with conventional observation technologies, such as leveling and global positioning systems (GPS), interferometric synthetic aperture radar (InSAR) does not need manual intervention and can monitor the surface changes both day and night and in all weather conditions [5][6][7].Simultaneously, the successful launch of multiple sensors has provided us with more data sources and has shortened the repeat period of single look complex (SLC) images [8].In particular, the free availability of Sentinel-1 data has further boosted the application of InSAR [9].As a result, InSAR has great advantages in deformation monitoring and has been successfully applied in the monitoring of various kinds of geological disasters, such as volcanic eruptions, earthquakes, landslides and subsidence [10][11][12][13][14][15][16].
Despite the many obvious advantages of InSAR, errors are an inevitable part of the interferometric phase.Phase decorrelation, such as spatial decorrelation and temporal decorrelation [17], usually causes low coherence and a low signal-to-noise ratio (SNR), which hinders accurate measurement and the retrieval of reliable ground deformation information.When microwave signals are transported in air and electromagnetic waves experience different atmospheric phase screens at different time acquisition, hence the atmospheric components are variable.The phases of interferograms are the difference of signals at two different times, as a result, the interferograms constructed in different atmospheric conditions will not be exactly the same, because of the atmospheric error.When this error is serious, it can reach the decimeter level [18,19].In addition, orbital inaccuracy and digital elevation model (DEM) deviation are also potential error sources [20].Therefore, how to reduce the effect of errors and thus attain an accurate deformation time-series is the key issue for InSAR technology [21].
Deformation time-series acquisition methods using InSAR technology have been developed over the last 15 years, the most representative method is small baseline subset (SBAS) technique proposed by Berardino et al. [22].The increased availability of images has further accelerated the development of different methods.According to the different observational objects, time-series acquisition methods can be divided into three categories.The first category takes the average phase and amplitude of all the targets in a resolution element as the pixel observation [22][23][24][25].The second category is persistent scatterer (PS) interferometry [26].The third category is known as temporally coherent point (TCP) InSAR [27].Compared with PS interferometry, this approach is not limited by requiring high coherence throughout all the time period, which increases the target density.Other methods are hybrids of the above three approaches.
The atmospheric error is a non-negligible part in the generation of a displacement time-series [20].Generally speaking, the atmosphere is made up of the troposphere and ionosphere.When dealing with C-band data, we are mainly concerned with the troposphere, which includes both an elevation-related component and a turbulent component.The elevation-related part can be diminished by a fitting function and the turbulent atmosphere resolution methods generally fall into two categories.One category is filter processing, where it is assumed that the deformation is made up of high-pass and low-pass components.The low-pass deformation is solved in the functional model and the high-pass atmospheric component is obtained from the residual phase by filtering technology.The other approach is to ignore the atmospheric effect by building a functional model between adjacent targets.
Mexico City suffers from severe land subsidence, which has received wide attention.Since the deformation phenomenon was first observed in the city in 1925 [28], there have been lots of studies using InSAR technology, which have all assumed that the deformation behaves in a linear manner and the nonlinear component has usually been ignored [29][30][31][32][33][34].However, the fact is that not all areas demonstrate a linear variation tendency and some regions display nonlinear displacement [35][36][37].
The Sentinel-1 mission [9], which followed the ERS-1, ERS-2 and Envisat satellites in the European radar observation system, is made up of a constellation of two C-band satellites-Sentinel-1A and Sentinel-1B-which were launched in April 2014 and April 2016, respectively.The Sentinel-1 constellation has a revisit time of six days, contains multiple imaging patterns and polarization modes and provides abundant information for Earth environmental monitoring [38,39].
According to the characteristics of nonlinear deformation and the atmospheric phase, this paper presents a new time-series generation method, which takes the nonlinear deformation and atmospheric phase as an independent unknown parameter related to the acquisition time in the functional model.Instead of extracting the nonlinear deformation and atmospheric phase from the residual phase, the proposed method separates these components with the singular spectrum analysis (SSA) method [40].
The rest of this paper is arranged as follows.Section 2 describes the proposed method, including the data preprocessing and the mathematical model.Section 3 details the study area and the InSAR data.Section 4 provides the experimental results.Sections 5 and 6 are the discussion and conclusions, respectively.

Data Preprocessing
The interferometric phase represents the distance change between two moments in the radar signal direction: where ϕ ij is the interferometric phase during time i and j; is the ground deformation in the line of sight (LOS); ϕ res,orb ij is the residual orbital error due to an inaccurate sensor position; ϕ res,dem ij is the DEM error; ϕ res,unw ij is the unwrapping error; ϕ atm ij is the atmospheric error; and ϕ noise ij is the summary effect of time decorrelation and random noise, such as sensor thermal noise.The resolution process involves separating ϕ disp ij from ϕ ij in the above formula and diminishing the effect of the other components.ϕ res,orb ij is the trend error in the InSAR interferogram, which is due to the imprecise orbit production and often behaves as systematic error.This is generally removed by masking the deformation area and applying polynomial fitting, the form of which is The interferometric phase induced by residual topography can be represented as , where λ is the microwave wavelength, R is the distance between the satellite and ground target in the LOS direction, θ is the incidence angle, B ⊥ is the spatial perpendicular baseline and ∆z is the difference between the real elevation and the applied DEM.When having a great influence, ϕ res,dem ij needs to be estimated in the functional model; otherwise, it can be neglected when it has less of a magnitude than the ground deformation.
A cycle phase of ϕ res,unw ij represents 2.8 cm for C-band sensors.Thus, the unwrapping error behaves as a gross error and needs to be preprocessed for more accurate deformation results.The detection and correction of the unwrapping error is based on the closed-loop residual phase [41].If the closed loop is constructed by ϕ ij , ϕ ja , ϕ ai , then the residual is: whereϕ disp ija and ϕ atm ija are the deformation and atmospheric phase in the closed-loop residual, which are relative to the attained times and should be zero.As a result, after dealing with the error of ϕ res,orb ija and ϕ res,dem ija , the unwrapping error of the interferometry ϕ res,unw ij can be checked when combined with multiple ϕ ija .A similar approach was explored in Reference [42] where the temporal coherence factor was introduced as a quality index of the unwrapped data in the framework of a minimum cost flow phase unwrapping algorithm.
The likelihood of unwrapping error increases with a larger interferometric phase gradient.Thus, as is shown in Figure 1, we can correct the unwrapping error based on the multiple closed-loop interference pairs by removing the orbital error and low-pass displacement from the original interferometric phase ϕ ij , unwrapping the residual phase and finally checking if the unwrapping error exists in ϕ ij .For further details, we refer the reader to Biggs et al. [41]

Mathematical Model
If the N interferograms,  ,  , ⋯ ,  are generated by conjugating multiplication of M+1 SLC images, the functional model is then: where  is the vector of the unknown parameters,  is the interferometric phase vector,  is the coefficient matrix and  is the residual vector.If the weighted matrix  is taken into account, according to the least-squares criterion, the final estimated parameters are: The above formula is the general mathematical model between the subsidence time-series and the interferogram phase, where the unwrapping error and orbital residual phase have been removed.Thus, how to build a proper functional relationship, define the weight values of the observations and solve the deformation parameters are the most important questions.
In this paper, according to different hypotheses of the ground surface change process, we test three schemes-the Usai method, the Berardino method and our proposed method-to generate the deformation time-series (Table 1).
Usai [23] computed firstly deformation time-series with InSAR technique, he did not consider the effect of other interference errors in the functional model and set the deformation of each time point as the unknown parameter: where  , is the interferogram phase,  is the deformation at time  and Δ is the residual.Although the time-series deformation is easily attained and the functional model is simple, the solution result  tends to show discontinuity in the least-squares adjustment.
To avoid discontinuity of the deformation time-series, Berardino et al. [22] set the deformation velocity instead of each time point displacement as the unknown parameter.In addition, the atmospheric error was filtered based on its spatio-temporal feature.
Equation ( 6) increases the restrictive condition and takes the deformation process as the polynomial form.The computation process is as follows: 1. Estimate the deformation parameters: ̅ ,  , ∆ .2. Remove the low-pass deformation and topography error from the original interferogram and unwrap the residual phase.3. Filter the atmospheric phase.4. Subtract the atmospheric phase from the original interferogram and then repeat steps 1-3 to refine the solved results.

Mathematical Model
generated by conjugating multiplication of M+1 SLC images, the functional model is then: where X is the vector of the unknown parameters, L is the interferometric phase vector, A is the coefficient matrix and V is the residual vector.If the weighted matrix P is taken into account, according to the least-squares criterion, the final estimated parameters are: The above formula is the general mathematical model between the subsidence time-series and the interferogram phase, where the unwrapping error and orbital residual phase have been removed.Thus, how to build a proper functional relationship, define the weight values of the observations and solve the deformation parameters are the most important questions.
In this paper, according to different hypotheses of the ground surface change process, we test three schemes-the Usai method, the Berardino method and our proposed method-to generate the deformation time-series (Table 1).
Usai [23] computed firstly deformation time-series with InSAR technique, he did not consider the effect of other interference errors in the functional model and set the deformation of each time point as the unknown parameter: where δϕ t B ,t A is the interferogram phase, d t B is the deformation at time t B and ∆ϕ res is the residual.Although the time-series deformation is easily attained and the functional model is simple, the solution result x tends to show discontinuity in the least-squares adjustment.
To avoid discontinuity of the deformation time-series, Berardino et al. [22] set the deformation velocity instead of each time point displacement as the unknown parameter.In addition, the atmospheric error was filtered based on its spatio-temporal feature.
Equation ( 6) increases the restrictive condition and takes the deformation process as the polynomial form.The computation process is as follows: 1. Estimate the deformation parameters: v, a, ∆a. 2. Remove the low-pass deformation and topography error from the original interferogram and unwrap the residual phase.3. Filter the atmospheric phase.4. Subtract the atmospheric phase from the original interferogram and then repeat steps 1-3 to refine the solved results.
The atmospheric phase in the interferogram is related to the atmospheric conditions when the signal goes through the air and the nonlinear displacement shows a low-pass characteristic in space and time.The new method sets the unknown parameters as the linear deformation v and the sum of the nonlinear deformation and atmospheric phase S.
where S t i is the sum of the nonlinear displacement and atmospheric phase at time t i .The atmospheric phase is obtained by the SSA method from S according to the features of the atmospheric phase and nonlinear displacement.For pixel i, the following processing steps: where, ϕ non disp is the time-series nonlinear displacement in pixel i, ϕ atm is the time-series atmospheric production in the pixel i, SSA(S i ) is the SSA processing of combination of atmosphere and nonlinear deformation, ϕ non−disp B and ϕ non−disp A are nonlinear displacement of pixel i at the time t B and t A respectively.
SSA, which was developed on the basis of Karhunen-Loève decomposition theory [43] and is not limited by sine wave characteristics, extracts the signal components that show significant oscillatory behavior and rebuilds the time-series by recognizing and strengthening the periodic signal by the use of the time-domain spectral characteristics method.For the detailed mathematical process, we refer the reader to Montillet [44].Due to its practicality, SSA has been extensively applied in oceanography and nonlinear dynamics fields.

Parameter estimation, filter
In general, the weights of the observations are considered to be equal when estimating the deformation and elevation error.Lauknes et al. [25] first employed the unit array as the initial value and then defined the weights iteratively, on the basis of the observation residuals vector.To avoid the residual value n i from being too small and the corresponding weight overlarge, they introduced a minor positive factor ε and rectified the weight matrix: The value of P is unit matrix, while k is equal to 1.Among the many different weighting functions, the IGG-III scheme divides the observation by quality into three categories and makes full use of the effective information, while limiting the suspicious information and eliminating the harmful information [45]: where ε k i is the standardization residual after k iterations, c 0 = 1 ∼ 1.5 and c 1 = 3.0 ∼ 8.0.In the experiment undertaken in this study, we applied the IGG-III scheme and set c 0 = 1.5 and c 1 = 4.

Geological Background
Mexico City, which has an average elevation of ~2200 m, is located on the southern plain of the Mexican Basin.Mexico City is surrounded by a number of active volcanoes [46].The soil of Mexico City is composed of highly compressible, weakly permeable lacustrine clay [47], the distribution coverage is displayed in Figure 2. Industrial and agricultural water consumption has led to excessive groundwater extraction in the alluvial volcanic debris aquifer, strata compression and depression and consolidation of the aquitard [48].
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 17 where  is the standardization residual after  iterations,  = 1~1.5 and  = 3.0~8.0.In the experiment undertaken in this study, we applied the IGG-III scheme and set  = 1.5 and  = 4.

Geological Background
Mexico City, which has an average elevation of ~2200 m, is located on the southern plain of the Mexican Basin.Mexico City is surrounded by a number of active volcanoes [46].The soil of Mexico City is composed of highly compressible, weakly permeable lacustrine clay [47], the distribution coverage is displayed in Figure 2. Industrial and agricultural water consumption has led to excessive groundwater extraction in the alluvial volcanic debris aquifer, strata compression and depression and consolidation of the aquitard [48].

Data Processing
The 34 Sentinel-1 scenes of Mexico City, from 2015/04/13 to 2016/09/10, were obtained from the European Space Agency (ESA) website.For the interferometric application of TOPSAR mode data,

Data Processing
The 34 Sentinel-1 scenes of Mexico City, from 2015/04/13 to 2016/09/10, were obtained from the European Space Agency (ESA) website.For the interferometric application of TOPSAR mode data, which is different from the strip scanning mode, the phase ramps of the SLC pair require registration accuracy of the milli-pixel order.In this study, we utilized initial registration based on intensity and accurate registration based on spectral diversity to remove the inconsistent Doppler centers between bursts and reduce the registration error to a level of one-thousandth of a pixel.An adaptive filter was then employed to reduce the phase noise, the topographic phase was removed by external SRTM DEM simulation and minimum cost flow phase unwrapping was used to recover the real displacement.Based on the GAMMA platform, 287 interferograms were selected from the interferogram pairs (Figure 3), with the limit of 180 days of temporal baseline and 120 m of spatial perpendicular baseline, obtaining pixels with a resolution of 150 × 150 m after 31 × 6 multi-look processing.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 17 which is different from the strip scanning mode, the phase ramps of the SLC pair require registration accuracy of the milli-pixel order.In this study, we utilized initial registration based on intensity and accurate registration based on spectral diversity to remove the inconsistent Doppler centers between bursts and reduce the registration error to a level of one-thousandth of a pixel.An adaptive filter was then employed to reduce the phase noise, the topographic phase was removed by external SRTM DEM simulation and minimum cost flow phase unwrapping was used to recover the real displacement.Based on the GAMMA platform, 287 interferograms were selected from the interferogram pairs (Figure 3), with the limit of 180 days of temporal baseline and 120 m of spatial perpendicular baseline, obtaining pixels with a resolution of 150 × 150 m after 31 × 6 multi-look processing.The influence of topography is defined as  ∆ =  ∆ ( • sin ) ⁄ , according to the parameters of the Sentinel-1 data.The distance between the satellite sensor and ground target  = 9.14 × 10 5 m; the center incidence angle in the study area  = 44.5°; the interferometric perpendicular baseline  < 120 m; and the maximum relative error of the Shuttle Radar Topography Mission (SRTM) DEM is 10 m.If we set ∆ = 10 m, the computed max( ∆ ) = 1.8 mm.Due to the real annual rate being larger than 10 cm, the impact of elevation error can be ignored in the latter solution.
Satellite positioning precision in orbit is a key factor for an accurate interferometric phase.ESA provides precise orbit ephemerides (POE) data for precise orbit determination (POD).In this study, we utilized the accurate orbit data to diminish the orbital effect.Satellite positioning precision in orbit is a key factor for an accurate interferometric phase.ESA provides precise orbit ephemerides (POE) data for precise orbit determination (POD).In this study, we utilized the accurate orbit data to diminish the orbital effect.

Annual Velocity
To diminish the effect of the atmosphere and the annual deformation rate could be obtained with the stacking method.Figure 4 shows the mean land subsidence between 2015/04/13 and 2016/09/10 generated with the 34 Sentinel-1 scenes, where the reference point have been selected in stable regions.
The results show that the main deformation areas are concentrated in the west.Therefore, the shortened study area not only reduces the workload but also displays only necessary deformation information.The land subsidence rate reaches a maximum value of −22 cm/year in the State of Mexico but other errors, such as atmospheric error, are still present.

Annual Velocity
To diminish the effect of the atmosphere and the annual deformation rate could be obtained with the stacking method.Figure 4 shows the mean land subsidence between 2015/04/13 and 2016/09/10 generated with the 34 Sentinel-1 scenes, where the reference point have been selected in stable regions.
The results show that the main deformation areas are concentrated in the west.Therefore, the shortened study area not only reduces the workload but also displays only necessary deformation information.The land subsidence rate reaches a maximum value of −22 cm/year in the State of Mexico but other errors, such as atmospheric error, are still present.

Mathematical Model
In this paper, three kinds of functional models are selected.Method I directly solves the deformation of each time by ignoring the effect of other errors.Method II supposes that the deformation process obeys a polynomial rule, calculates the polynomial coefficients and separates the atmospheric phase by filter processing.Method III computes the annual deformation rate and takes the combination of the atmospheric phase and nonlinear displacement as an independent parameter related to the SLC image acquisition time.
According to previous studies [31,32,36], the regions of boxes a and b in Figure 4 are stable districts, which means that the real-time deformation can be considered to be zero.The time-series

Mathematical Model
In this paper, three kinds of functional models are selected.Method I directly solves the deformation of each time by ignoring the effect of other errors.Method II supposes that the deformation process obeys a polynomial rule, calculates the polynomial coefficients and separates the atmospheric phase by filter processing.Method III computes the annual deformation rate and takes the combination of the atmospheric phase and nonlinear displacement as an independent parameter related to the SLC image acquisition time.
According to previous studies [31,32,36], the regions of boxes a and b in Figure 4 are stable districts, which means that the real-time deformation can be considered to be zero.The time-series deformation was obtained with the different functional models and weighting functions.The statistical results obtained for the stable areas represent the accuracy of the corresponding methods.Figure 5 show the deformation time-series of the different functional models and weighting functions in the stable regions.When using the equal-weighting function, the result of Method III is the closest to zero and therefore has the highest accuracy.This is followed by Methods I and II, which shows that the selected functional model must conform to the deformation characteristics of the actual study area.When using the same functional model, the results obtained using unit weights and adaptive weights are depicted in different colors in Figure 5.The different weighting functions in Methods I and III obtain roughly the same results.However, for Method II, the difference between the two weighting functions can reach approximately 1 cm.Based on the above analysis, Method III has the best precision and is not affected by the weighting function and the equal-weight result of Method III can be taken as the final time-series deformation.The results displayed in Figure 6 are based on Method III with equal weights.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 17 deformation was obtained with the different functional models and weighting functions.The statistical results obtained for the stable areas represent the accuracy of the corresponding methods.Figure 5 show the deformation time-series of the different functional models and weighting functions in the stable regions.When using the equal-weighting function, the result of Method III is the closest to zero and therefore has the highest accuracy.This is followed by Methods I and II, which shows that the selected functional model must conform to the deformation characteristics of the actual study area.When using the same functional model, the results obtained using unit weights and adaptive weights are depicted in different colors in Figure 5.The different weighting functions in Methods I and III obtain roughly the same results.However, for Method II, the difference between the two weighting functions can reach approximately 1 cm.Based on the above analysis, Method III has the best precision and is not affected by the weighting function and the equal-weight result of Method III can be taken as the final time-series deformation.The results displayed in Figure 6 are based on Method III with equal weights.

Time-Series Deformation
In Figure 6, we display the surface displacement time-series between 2015/04/13 and 2016/09/10 in the SAR coordinate system.Two typical subsidence zones can be observed.The maximum cumulative displacement in the study area exceeds −30 cm and some displacement areas have no results due to decorrelation of the agricultural region.To view the details of some specific points, we chose three points-P1, P2 and P3-in the cumulative deformation map (See Figures 2 and 7 for the location and deformation time-series).P1 and P2 are located in the first deformation area and P3 is chosen in the second deformation region (see locations on Figure 4).Analysis of the time-series shows that the movement of the deformation region is very close to linear motion.

Time-Series Deformation
In Figure 6, we display the surface displacement time-series between 2015/04/13 and 2016/09/10 in the SAR coordinate system.Two typical subsidence zones can be observed.The maximum cumulative displacement in the study area exceeds −30 cm and some displacement areas have no results due to decorrelation of the agricultural region.To view the details of some specific points, we chose three points-P1, P2 and P3-in the cumulative deformation map (See Figures 2 and 7 for the location and deformation time-series).P1 and P2 are located in the first deformation area and P3 is chosen in the second deformation region (see locations on Figure 4).Analysis of the time-series shows that the movement of the deformation region is very close to linear motion.

Atmospheric and Nonlinear Results
Figure 8 shows the atmospheric results relative to the first SLC image acquisition time and the reference point calculated by Method III.Due to the lack of atmospheric products with a sufficient resolution, we chose four interferograms that are not contaminated by other error, for checking and validating the precision of the atmospheric phase and nonlinear deformation The selected interferograms in Figure 9 were not involved in the calculation of the deformation time-series and the linear deformation has been removed.Thus, in theory, only nonlinear deformation and atmospheric phase are left.However, interferogram 20160326_20160817 demonstrates that residual nonlinear displacement is still present, although the atmospheric phase has been effectively removed.The other atmospheric and interferometric phases show a similar spatial distribution and the maximum value reaches 5 cm, with the nonlinear deformation showing a smaller value.The results show that the atmospheric phase can be correctly computed and separated from the nonlinear deformation using the SSA method, to some extent.
Table 2 lists the mean and standard deviation values of the atmospheric phase and nonlinear displacement in stable areas a and b.The change of the residual interference phase represents the effectiveness of the nonlinear displacement and atmospheric phase results obtained with Method III.If the phase of IFG-NONL is larger than that of IFG, then the NONL result can be considered suspicious.In area a, except for the precision of interferogram 20160326_20160817 being reduced by a few millimeters, the mean accuracies of the other interferograms are improved by 96%, 67% and 99%, respectively.In area b, the IFG and IFG-ATM-NONL values show roughly the same precision and the mean accuracies of the other interferograms are improved by 87%, 79% and 61%.Further analysis of the results of the interferograms of 20160326_20160817 and 20160501_20160513 shows that the differences are of the millimeter magnitude.The topography error may have introduced contamination and the results can be considered as acceptable.Therefore, the results of Method III can be deemed as credible, although some interferograms show a slight reduction in precision.

Atmospheric and Nonlinear Results
Figure 8 shows the atmospheric results relative to the first SLC image acquisition time and the reference point calculated by Method III.Due to the lack of atmospheric products with a sufficient resolution, we chose four interferograms that are not contaminated by other error, for checking and validating the precision of the atmospheric phase and nonlinear deformation.
The selected interferograms in Figure 9 were not involved in the calculation of the deformation time-series and the linear deformation has been removed.Thus, in theory, only nonlinear deformation and atmospheric phase are left.However, interferogram 20160326_20160817 demonstrates that residual nonlinear displacement is still present, although the atmospheric phase has been effectively removed.The other atmospheric and interferometric phases show a similar spatial distribution and the maximum value reaches 5 cm, with the nonlinear deformation showing a smaller value.The results show that the atmospheric phase can be correctly computed and separated from the nonlinear deformation using the SSA method, to some extent.
Table 2 lists the mean and standard deviation values of the atmospheric phase and nonlinear displacement in stable areas a and b.The change of the residual interference phase represents the effectiveness of the nonlinear displacement and atmospheric phase results obtained with Method III.If the phase of IFG-NONL is larger than that of IFG, then the NONL result can be considered suspicious.In area a, except for the precision of interferogram 20160326_20160817 being reduced by a few millimeters, the mean accuracies of the other interferograms are improved by 96%, 67% and 99%, respectively.In area b, the IFG and IFG-ATM-NONL values show roughly the same precision and the mean accuracies of the other interferograms are improved by 87%, 79% and 61%.Further analysis of the results of the interferograms of 20160326_20160817 and 20160501_20160513 shows that the differences are of the millimeter magnitude.The topography error may have introduced contamination and the results can be considered as acceptable.Therefore, the results of Method III can be deemed as credible, although some interferograms show a slight reduction in precision.

Mathematical Models
Method III can not only obtain the land subsidence time-series but also independently computes the nonlinear displacement and atmospheric phase as an unknown parameter in the functional model.Furthermore, the comparison with the other two methods in stable areas shows that the proposed method performs better than the others.
The assumption of low-pass deformation is important for the deformation solution.The general deformation trend is linear in the study area.Method II computing the low-pass deformation with polynomial form can introduce other errors into the solution.
When the functional model roughly fits the deformation process, the impact of the weighting function is not great.We used unit weights and adaptive weights in the three methods to check the influence of the weighting function.Except for Method II, whose functional is not linear, the results of the other methods using different weighting functions do not demonstrate significant differences.

Spatio-Temporal Changes of the Deformation in Mexico City
The temporal evolution of the deformation identified in this study is roughly linear and is consistent with that reported in former studies [30,31,34,35].The results of our study indicate that the maximum cumulative deformation exceeds −30 cm in the metropolitan area of Mexico City (see Figure 6).
By studying the relationship between the deformation characteristics and the geological structure of Mexico City, we can find that the deformation areas are mostly distributed in the lacustrine clay region.In addition to the geological environment, human activity is another important factor.Increased population and rapid urbanization has led to excessive pumping of groundwater and the non-resilient compression in the lacustrine aquitard has caused continuous subsidence of the Earth's surface.

Conclusions
The free availability of the Sentinel data offers an opportunity for large-area and slow surface deformation research.The design of shorter temporal and spatial baselines can also help to reduce the incoherence problem of InSAR technology.In the meantime, reasonable data preprocessing provides a good basis for time-series acquisition.
The spectral registration overcomes the problem of different Doppler rates between bursts.The correction of the unwrapping error in the data preprocessing avoids gross error and provides good conditions for time-series calculation.At the same time, the negligible topography error reduces the number of unknown parameters and helps with the deformation acquisition.
A suitable mathematical model is the most important part in deformation time-series generation.The selected functional model needs to be consistent with the actual deformation characteristics.Once the functional model is correctly selected, unit weight and adaptive weight models make no difference to the final results.Compared with Methods I and II, Method III is more reliable, not only for obtaining nonlinear deformation but also for obtaining the phase information of the atmosphere, which can play an important role.It can be seen from the graphs in this paper that the atmospheric phase in some time is close to 5 cm, which can be ignored in large deformation regions but in smaller deformation areas, calculation of the atmospheric phase becomes more important.Therefore, Method III is beneficial, not only to the study of the time-series but also to the analysis of the atmospheric conditions in the study area.
Consistent with the former annual deformation rates reported for Mexico City, the maximum rate obtained in this study reached −20/year cm in the LOS direction between 2015/04/13 and 2016/09/10.The deformation characteristics of points P1, P2 and P3 were also approximately linear.

Figure 1 .
Figure 1.Framework of interferogram unwrapping, error detection and correction.The boxes represent the detailed operations, the parallelograms indicate the processed results and the arrows point to the next process.

Figure 1 .
Figure 1.Framework of interferogram unwrapping, error detection and correction.The boxes represent the detailed operations, the parallelograms indicate the processed results and the arrows point to the next process.

Figure 2 .
Figure 2. The coverage of Mexico City soil composition (plot according to [31]) and study area.The black box in the upper right figure indicates the scope of the Sentinel-1 dataset and the blue rectangle is the study area.In the main figure, the purple and red lines represent the scope of the transitional zone and the lacustrine clay border and the blue box denotes the study area.The lacustrine clay border defines the scope of clayed layer overlaying sand and gravel aquifer, transitional zone distinguishes the different geotechnical units.

Figure 2 .
Figure 2. The coverage of Mexico City soil composition (plot according to [31]) and study area.The black box in the upper right figure indicates the scope of the Sentinel-1 dataset and the blue rectangle is the study area.In the main figure, the purple and red lines represent the scope of the transitional zone and the lacustrine clay border and the blue box denotes the study area.The lacustrine clay border defines the scope of clayed layer overlaying sand and gravel aquifer, transitional zone distinguishes the different geotechnical units.

Figure 3 .
Figure 3.The baseline distribution of the interferograms, where the blue triangles show the times of the SLC image acquisitions and the purple lines display the formed interferograms.

Figure 3 .
Figure 3.The baseline distribution of the interferograms, where the blue triangles show the times of the SLC image acquisitions and the purple lines display the formed interferograms.The influence of topography is defined as d ∆z = B ⊥ •z/(r∆ sin θ), according to the parameters of the Sentinel-1 data.The distance between the satellite sensor and ground target r = 9.14 × 10 5 m; the center incidence angle in the study area θ = 44.5 • ; the interferometric perpendicular baseline B ⊥ < 120 m; and the maximum relative error of the Shuttle Radar Topography Mission (SRTM) DEM is 10 m.If we set ∆z = 10 m, the computed max(d ∆z ) = 1.8 mm.Due to the real annual rate being larger than 10 cm, the impact of elevation error can be ignored in the latter solution.Satellite positioning precision in orbit is a key factor for an accurate interferometric phase.ESA provides precise orbit ephemerides (POE) data for precise orbit determination (POD).In this study, we utilized the accurate orbit data to diminish the orbital effect.

Figure 4 .
Figure 4. Subsidence velocity derived by the stacking method (unit: cm/year), where the black rectangles represent stable areas, the purple star indicates the reference point, the blue triangles are the feature points for the time-series deformation analysis and the blue box is the study area, the grey line indicates the administrative boundaries.

Figure 4 .
Figure 4. Subsidence velocity derived by the stacking method (unit: cm/year), where the black rectangles represent stable areas, the purple star indicates the reference point, the blue triangles are the feature points for the time-series deformation analysis and the blue box is the study area, the grey line indicates the administrative boundaries.

Figure 5 .
Figure 5. Time-series statistical chart in the stable areas by different functional models and weighting functions.Unit weight and adp weight represent equal and adaptive weight distributions, respectively.The left column shows the results of area a and the right column shows the results of area b.The figures in the three lines show the time-series results of the Usai method, the Berardino method and the proposed method, respectively.

Figure 5 .
Figure 5. Time-series statistical chart in the stable areas by different functional models and weighting functions.Unit weight and adp weight represent equal and adaptive weight distributions, respectively.The left column shows the results of area a and the right column shows the results of area b.The figures in the three lines show the time-series results of the Usai method, the Berardino method and the proposed method, respectively.

Figure 6 .
Figure 6.Displacement time-series along LOS direction with unit weights by Method (unit: cm).The negative value represents the increased distance between ground target and satellite antenna direction.

Figure 6 .
Figure 6.Displacement time-series along LOS direction with unit weights by Method (unit: cm).The negative value represents the increased distance between ground target and satellite antenna direction.

Figure 7 .
Figure 7. Time-series results of points P1, P2 and P3 (see Figure 4) between 2015/04/13 and 2016/09/10, where the horizontal axis represents time and the vertical axis indicates the deformation in the LOS direction.

Figure 7 .
Figure 7. Time-series results of points P1, P2 and P3 (see Figure 4) between 2015/04/13 and 2016/09/10, where the horizontal axis represents time and the vertical axis indicates the deformation in the LOS direction.

Figure 8 .
Figure 8. Atmospheric time-series generated by Method III in the LOS direction (unit: cm).Figure 8. Atmospheric time-series generated by Method III in the LOS direction (unit: cm).

Figure 8 .
Figure 8. Atmospheric time-series generated by Method III in the LOS direction (unit: cm).Figure 8. Atmospheric time-series generated by Method III in the LOS direction (unit: cm). .

Table 1 .
Summary of the methods used to obtain the deformation time-series.