Separating and Unwrapping of Residual Topographical and Time-Series Deformation Phases in PS-InSAR Based on Zero-Temporal Baseline Combination and Time-Domain Differentiation

To separate and unwrap the residual topographical and deformation components from the wrapped phases are the key issues for persistent scatterer interferometric synthetic aperture radar (PS-InSAR) technique. The current methods, such as grid-search method based on the assumption that the deformation in time-domain is linear, or the least squares ambiguity decorrelation adjustment method that needs priori information to construct the pseudo-observations, etc., are still difficult to tradeoff accuracy against efficiency. In this article, we propose a new method to separate and unwrap the differential phases for PS-InSAR by dividing them into residual topographical phase and differential deformation phases and performing the phase unwrapping separately. First, we calculate phase ambiguities caused by residual topographical error by zero-temporal baseline combination and single parameter grid-search method. The calculated residual topographical error is used to separate the differential phases dominated by deformations from the original differential phases of PS-InSAR. Then, we proposed time-domain differentiation to unwrap differential phases dominated by deformations without any assumption of deformation model in time-domain. The performance of the proposed method is tested by both simulated and real datasets. The experimental results show that, compared with classic PS-InSAR, the accuracy of residual topographical error estimation is comparable; however, there is an accuracy improvement of 34.4% for the deformation parameter estimation and 35.7% for nonlinear deformation time-series extraction, as well as an improvement of 20–40 times in computational efficiency.


I. INTRODUCTION
A S AN advanced interferometric synthetic aperture radar (InSAR) technique, persistent scatterer (PS) InSAR is widely used to detect ground deformation with millimetric accuracy in time series on a high density of PS points [1], [2], [3].However, since it was proposed, PS-InSAR has been affected by the problem of components separating from the wrapped phase.Phase separation resolution is a critical step in PS-InSAR [4], [5].The performance of phase separation can directly determine the accuracy and reliability of PS-InSAR products.
As we know, the differential interferometric phase can be seen as combinations of the deformation phase, residual topographical phase, orbital error phase, atmospheric delay phase, and other decorrelation noise phases.Phase separation means extracting the parts of interest (such as the deformation phase and residual topographical phase) from the mixed phase.Because the interferometric phase is wrapped, a phase unwrapping process is performed first or simultaneously during the phase separation.Phase unwrapping is a process of recovering unambiguous phase cycles from the phase data that are measured in modulo 2π.Traditional phase unwrapping methods, such as the Branch-and-Cut algorithm [6] and MCF method [7], are not performed well for PS-InSAR because PS-InSAR is worked at sparse and irregularly sampled PS points.In order to achieve the phase separating and unwrapping, four kinds of parameterized methods have been proposed.First, the wrapped phase was transformed by Euler's formula to avoid the effect of phase ambiguities, and deformation parameters could be searched directly in the solution space, which is called the bidimensional period-diagram method [3] or grid-search method [8].However, this kind of method can only calculate deformation rate parameters and residual topographical parameters for arc phase and does not really unwrap the phase and phase ambiguities may still exist when we extract the nonlinear deformation time series from residual phase because underfitted deformation model may cause large residual arc phase, which may affect the accuracy of parameter calculation and result in the existence of phase ambiguities in residual arc phase.In addition, the accuracy of the deformation parameters maybe affected by the upper and lower bounds of the solution space, and partial nonlinear deformation may be mistakenly removed during the atmospheric phase correction process by spatial-temporal filtering.Second, the least squares ambiguity decorrelation adjustment [9], [10] and its improved methods [11], [12], [13] were applied to PS-InSAR to estimate both ambiguities and deformation parameters simultaneously by using the integer least square method and pseudo-observations.Because of the lack of a priori information, it is difficult to construct the pseudo-observations and weight matrices.Third, three-dimensional phase unwrapping methods were applied to PS-InSAR with a representative method proposed by [14] and [15].This method introduced temporal-dimension information and improved the accuracy of phase unwrapping.However, these methods require assumptions that are not always valid.Moreover, the efficiency of these methods is low, and the reduceresolution process can result in new errors during parameter estimation.Fourth, the TCP-InSAR technique [16] was proposed to estimate deformation parameters by eliminating PS arcs that contain phase ambiguities.The process of elimination may break the stability of the PS network and create some individual subnetworks, which can cause the failure of parameter estimation.The above four methods also have some common problems, such as low efficiency and excessive time and computer memory consumption.
The main principle of these parameterized methods is to model phase components so that we can divide the deformation time series into modeled and nonmodeled parts and calculate them separately, which can be affected by the inapplicability of the model [17].In fact, real surface deformation is very complicated in both temporal and spatial domains and hard to be fit with a unified model [18].Nonmodeled deformation often partially lose during the temporal-spatial filtering process, and the real deformation cannot be fully restored by these parameterized methods [19].To avoid the effect of badly underfitted deformation models, some nonparametric methods have been proposed [20], [21], [22], but these methods are complex and can exhibit low robustness.Therefore, separating, unwrapping the phase, and extracting deformation time series accurately and efficiently remain a key issue in further advancing PS-InSAR applications.
Limited by the number of available SAR satellites and their revisit periods, previous articles have focused more on the advantages of high space resolution and wide ground coverage of SAR images.However, with the rapid development of SAR satellite technologies, some SAR satellites were launched in constellation mode, and the revisit period of SAR satellites has been shortened dramatically to several days.Intensive time sampling can provide abundant observational information for SAR datasets and improve the ability to depict actual ground displacement, which can provide a practical basis for phase separating, unwrapping, and the recovery of deformation sequences from the time dimension.
In the general phase model for PS-InSAR, the differential interferometric phase consists of the deformation phase, residual topographical phase, atmospheric phase, orbital error phase, and other noise phases.Accordingly, ambiguities of the differential interference phase are contributed by the above constituent.When we construct PS arcs for PS-InSAR, common space phase errors can be counteracted, and arc phase ambiguities can be considered to be only caused by deformation and residual topographical error.The residual topographical phase is irregular in the time series and can be functionalized with the vertical baseline, but the deformation phase accumulates gradually in the time series and is wrapped regularly.The differential expression of the deformation phase and residual topographical phase in the time dimension offers a new approach to separating and unwrapping phase for the PS-InSAR technique.In this article, we proposed a new phase separating, unwrapping, and deformation time-series extraction method for PS-InSAR by using intensive SAR datasets.This method follows the classic PS-InSAR framework and separates the deformation phase and residual topographical phase from the PS arc phase.Phase ambiguities caused by residual topographical error and ground deformation are calculated, respectively.Then, the weighted least squares method is used to estimate the deformation parameters and thus extract the deformation time series for the PS points.
The rest of this article is organized as follows.Section II introduces the proposed method in detail.The performance of our new method is tested by using synthetic datasets and real SAR datasets in Section III.Some factors that affect the estimation accuracy and success rate of this new method are discussed in Section IV.Finally, Section V concludes this article.

II. METHODOLOGY
Different from the classic methods, we deal with the phase separating and unwrapping problem by considering phase compositions and their change rules directly.The PS arc phase contains differential deformation phase, differential residual topographical phase, and differential noise phase.In our method, the zero-temporal baseline combination method and single parameter grid-search method are used to separate the differential residual topographical phases and calculate phase ambiguities caused by residual topographical error for the arc.Then, a successive difference for the arc phase in the time dimension is engaged to detect and correct the wrapped phase caused by ground deformation.Finally, the weighted least squares method is used to estimate deformation parameters and extract deformation time series for PS points.Section II-A-D describes these innovative steps in detail, and Section II-E shows the improved PS-InSAR processing chain compared with classical PS-InSAR technique.

A. Separating and Unwrapping of Residual Topographical Phase Based on Zero-Temporal Baseline Combination
Given N differential interferograms formed by N SAR images (containing a master-master differential interferogram), the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.interferometric phase of PS point x in the ith differential interferogram can be written as a combination of different contributions as follows: where Φ represents the wrapped phase, ϕ represents the unwrapped phase, W [•] is the wrapping operator, and ϕ defo , ϕ topo , ϕ orb , ϕ atm , and ϕ noise are the deformation phase, residual topographical phase, orbital error phase, atmospheric delay phase, and noise phase, respectively.
After modeling, (1) can be expressed as where λ, R, and θ are the wavelength, satellite-target distance, and incidence angle of PS point x for the sensor, respectively; t i and B ⊥ i are the temporal baseline and perpendicular baseline of the ith differential interferogram, respectively; v(x) and ε h x are the linear deformation rate and residual topographical error, respectively, due to inaccurate DEM data of PS point x; and ϕ res x,t i is the residual phase contributed by atmospheric delay, nonlinear deformation, noise error, orbital error phase, and other obstacles.
After constructing the PS network, the differential phase for an arc that connects two PSs x and y can be derived as (3) where ΔΦ is the phase difference between two neighboring pixels x and y; Δv and Δε h are the differential residual topographical error and differential linear deformation rate of PS pair x and y, respectively; and ϕ res x−y,t i is the differential residual phase.Some error phases, such as atmospheric delay and orbital error phases, can be counteracted [23] because PSs x and y are so nearby in space that those error phases are nearly uniform.Thus, ϕ res x−y,t i is contributed by differential nonlinear deformation and noise.
The differential phase for arc ΔΦ(x, y, t i ) is still a wrapped phase.The noise phase of PSs is very small, which cannot lead to phase wrapping.Phase ambiguities of the PS arc phase are caused by differential residual topographical error and differential deformations.The unwrapped differential phase for the arc can be described by where k def x−y and k topo x−y are the differential phase ambiguities caused by differential deformation and differential residual topographical error of PS pair x and y, respectively.
The differential deformation phase and differential residual topographical phase show different characteristics in the time dimension, so we can separate them from the arc phase.The values of differential residual topographical phase are related only to the geometric parameters of the SAR sensor.PS-InSAR is suitable for continuous and slow-moving ground deformation, so we can assume that the deformation values within adjacent equal time intervals are uniform.Under this assumption, the zero-temporal baseline combination method and single parameter grid-search method are used to separate the differential residual topographical phase and calculate phase ambiguities caused by residual topographical error for the arc.
If the master image is M th in N SAR images, the wrapped differential phase dataset for an arc can be represented as We make differentiation again for this dataset in the time dimension, and the differential dataset is Then, a zero-temporal baseline combination process is performed and a diagrammatic map of pseudophase formation processing is shown in Fig. 1.More specifically, we search within a certain search range (e.g., 30 days forward and backward) in for each time-dimension differential phase, and then we select those that are expected to have the same temporal baseline or multiples temporal baseline in search range to center phase.Then, a pseudophase can be formed for this arc using where ΔΦ a and ΔΦ b are two time-dimensional differential phases with equal or multiple temporal baselines, which means n a Δ t a = n b Δt b and n a , n b ∈ ±1, ±2.The small combinatorial number (n a , n b ) is used to avoid the overexpansion of phase errors, such as residual APS or orbit errors.According to the assumption above, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where ΔΦ a−b is the pseudophase with a "zero" temporal baseline and n a ΔB ⊥ a − n b ΔB ⊥ b is the pseudoperpendicular baseline.
After traversing for all time-dimensional differential phases, a series of pseudophases are obtained and used as observations to calculate a single parameter Δε h by the single-parameter gridsearch method.Then, we can calculate the differential residual topographical phase for the arc by Phase ambiguities k topo x−y caused by residual topographical error can also be calculated easily by measuring in modulo 2π red for Δϕ topo .Therefore, we can complete phase unwrapping of the differential residual topographical phase.

B. Separating and Unwrapping of Deformation Phases Based on Time-Domain Differentiation
In this section, the differential residual topographical phase is first subtracted from the wrapped differential phase dataset for an arc.The differential deformation phase (with noise phase) can be derived as As previously stated, the deformation phase continuously varies in the time dimension due to the intensive acquisition of SAR images, and the phase change is dramatic when wrapping occurs (see Fig. 2); therefore, we can detect and recover the phase wrapping according to the continuous change process of the differential deformation phase in the time dimension.
Concretely, the wrapped differential deformation phase dataset for an arc is We make a difference again for this dataset in the time dimension, and the differential dataset is A phase change threshold (e.g., 1.5π) is used to judge if phase wrapping occurred.Each time-dimensional differential phase is judged, and the corresponding wrap value is added or subtracted from the differential deformation phase.After this process, the differential deformation phase can be unwrapped preliminarily.
Because of the noise phase that is contained in the differential deformation phase, or the absence of SAR images which cause unevenly distributed time intervals, phase unwrapping for the differential deformation phase may contain mistakes.Considering that phase wrapping is always countertrending for deformation time series [see Fig. 2(a)], but the noise phase or unevenly distributed time intervals usually do not change the deformation trend [see Fig. 2(b)], a correction process is performed in our method.First, a low-pass filter is used separately for several differential deformation phases before and after the moment of phase wrapping, as judged by the phase change threshold.Then, slopes are calculated before, during and after the moment of the phase wrapping for the differential deformation phase, and we denote them by L 1 , L 2 , and L 3 .If L 1 and L 3 have the same plus/minus sign but L 2 is different, phase wrapping occurred at this moment.Otherwise, L 1 , L 2 , and L 3 have the same sign, phase wrapping is misjudged, and we need to correct it.
After this correction process, the phase ambiguities caused by deformation can be calculated precisely.

C. Parameters Estimation for Arcs
After subsections B and C, k def x−y and k topo x−y are calculated precisely, and we can restore the unwrapped differential phase for the arc by ( 4), which can be derived as follows: The least squares method is used by taking the unwrapped differential phase for the arc as an observation, and the differential residual topographical error Δε h and differential linear deformation rate Δv of PS pair x and y can be estimated exactly.The multi-image phase coherence proposed by [2] and [3] can also be obtained by where N is the number of interferograms; j is an imaginary unit, and j = √ −1 ; Δϕ res i can be derived using (10).

D. Nonparametric Estimation of Time-Series Deformation
Once we obtain the exact solution of arc parameters, the unwrapped differential residual topographical phase can be separated accurately, and the rest of the unwrapped arc differential phase contains the unwrapped differential deformation phase and noise phase.Some arcs with low γ (<0.7) are discarded.Then, the weighted least squares network adjustment method is performed to obtain the residual topographical error and deformation rate for all PS points.
For the PS network with P points and Q arcs, error equations can be formed as follows: where V is the correction matrix; B, X, and L are the design matrix, parameter matrix, and observation matrix, respectively; and B, X, and L have the following form: Furthermore, we define the weight matrix as follows: According to the methodology of the weighted least squares method, given the known starting data of the global reference point, the parameter can be calculated by We regard the differential residual topographical error and differential linear deformation rate calculated in subsection C as observations, respectively, and the residual topographical error and linear deformation rate can be estimated for each PS point.
In the same way, we can calculate the deformation phase for each PS point by regarding the rest of the unwrapped arc differential phase as observations and perform for every differential interferogram.Because master-master interferograms do not have phase wrapping and ground deformation, we can revise the deformation phase time series by subtracting the master-master deformation phase value for each PS point.After low-pass filtering to eliminate the noise phase and phase transformation, we can extract the deformation time series successfully.

E. Improved PS-InSAR Processing Chain
Compared with classic PS-InSAR, the workflow of the proposed approach can be summarized as follows (see Fig. 3).
1) We register SAR datasets and construct interferograms with a single-master model.Moreover, PSs are selected, and the network is constructed.

TABLE I BASIC PARAMETERS OF THE DATASETS CAPTURED FROM SENTINEL-1 SENSORS
2) We separate the differential residual topographical phase from the arc differential phase and calculate phase ambiguities caused by residual topographical error.3) We separate the differential deformation phase from the arc differential phase and calculate phase ambiguities caused by ground deformation.4) The differential residual topographical error and differential linear deformation rate of arcs can be estimated exactly.5) We estimate residual topographical error and linear deformation rates and extract deformation time series for each PS point by using the weighted least squares network adjustment method.

A. Synthetic Dataset Test
The test is first carried out on 69 synthetic interferograms with 512 × 512 pixels simulated from 69 C-band SAR images from May 2015 to April 2018.The geometric parameters of the SAR images are set with reference to Sentinel-1, as shown in Table I.In these interferograms, the maximum perpendicular baseline is 147 m, and the maximum adjacent time interval of the SAR images is 60 d.A total of 9968 PSs are selected that were evenly distributed in images, and 29 805 arcs are constructed in a Delaunay network.During the simulation, we simulate four types of phase signals: deformation, residual topographical error, atmospheric artifacts, and phase noise.The linear deformation pattern is simulated by a peak function with a range of deformation rates from −3 to 1 cm/a, as shown in Fig. 4(a).Moreover, a cyclic deformation is added in the northeast corner of images with a fixed period of 1 year and variable amplitude in space ranging from 0 to 2 cm, as shown in Fig. 4(b).The residual topographical errors are randomly distributed and range from −40 to 40 m for all PS points, as shown in Fig. 4(c).Then, turbulent atmospheric signal based on Kolmogorov turbulence theory [24] is simulated for each image, and random phase noise is added to all images with a mean of 15°.
Under these conditions, the new method is used to estimate the linear deformation rate and residual topographical error for each PS point, as shown in Fig. 5. Based on Fig. 5, we find that the calculated results are highly consistent with the simulated Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.values.More than 98% of PSs have a deformation rate bias within 1 mm/a, and more than 86% of PSs have a residual topographical error bias within 5 m.The RMSEs of the deformation rate calculated value and residual topographical error calculated value are 0.43 mm/a and 3.66 m, respectively.Fig. 5(d) shows that residual topographical error biases in the nonlinear deformation region are larger than those in other regions because residual topographical error estimation is sensitive to the residual phase, and we use a linear model to estimate nonlinear deformation, which can lead to a large residual phase.
Furthermore, deformation time series are extracted by both classical PS-InSAR and new method with same time series filter method (triangular filtering method) and same filter window, and the results of three pixels [P1, P2, and P3, which are shown in Fig. 5(a)] are shown in Fig. 6.The RMSEs of three points are 4.5, 5.3, and 2.4 mm for new method, but 7.0, 6.1, and 2.9 mm for classical method.The results show that compared with the classical method, the new method is competent for linear deformation but an improvement of 35.7% for nonlinear deformation.This may because linear model is used in classical PS-InSAR technique but new method does not relay on deformation model during deformation time-series estimation.We can also find that the deformation time-series results of new method may be slightly inconsistent in P1.The inconsistence appeared at a time point of 1 year after the initial observation in which about 8 mm displacement occurred during 28 days.Both the classic PS-InSAR technique and the new method were utilized in this experiment, and the results are shown in Fig. 8(a) and (b).The difference between the two methods is shown in Fig. 8(c).The distribution and magnitude were consistent for the two results, and 99.89% of the differences in the deformation rate were within ±2 mm/year.The intercomparison results [Fig.8(d)] between the two methods show that the correlation is 0.97 and the RMSE of the difference is 0.33 mm/year, revealing that the deformation rate mapped by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the new method is in good agreement with the classic method and that the new method is an alternative option in terms of monitoring accuracy.

A. Superiority in Accuracy Compared With Classic PS-InSAR
When we use classic PS-InSAR, the parameters are estimated by the grid-search method for the wrapped arc phase.The estimation may be affected by the wrapped phase, and the differential deformation rate and differential residual topographical error interact during searching.Parameter accuracy depends on the degree of agreement between the model and the actual situation.In the new method, we separate the differential deformation phase and differential residual topographical phase from the arc phase and unwrap them.Parameters are estimated directly by the unwrapped phase, and thus, parameter accuracy is improved for arc calculation and point calculation.
To show the superiority of our method in terms of accuracy, comparisons of classic PS-InSAR [25] and the new method have been performed for both arc parameter calculations and point parameter calculations by using the synthetic datasets described in Section III-A.Bias histograms of the parameters are shown for both arcs [Fig.9(a) and (b)] and points [Fig.9(c) and (d)].The means and STDs of these results are given in Fig. 9. Fig. 9 shows that residual topographical error biases are similar for both arcs and points, which may be because the influence of model errors is similar for the grid-search method in classic PS-InSAR and the least squares method in the new method when we calculate residual topographical error.Compared with the classic method, the STDs of the deformation rate parameter are reduced by 29.17% and 34.38% for arcs and points when we use the new method, respectively, revealing that the new method can achieve a more accurate deformation rate.We also find that the point parameter accuracy is lower than that of the arc results, which may be because of the misestimation of multi-image phase coherence, which is likely to increase random errors for phase observable and deformation products due to inappropriate weighting [26].

B. Superiority in Efficiency Compared With Classic PS-InSAR
Low efficiency is a major limitation in classic PS-InSAR.Although a novel block method [27] was proposed to solve this problem, it pays more attention to the space scale rather than the phase resolution.In classic PS-InSAR, the grid-search method is used to calculate arc parameters.The solution space is a 2-D space, and each parameter combination is tried during the search process.In this article, after simple phase combination and difference, 2-D space can be compressed into 1-D, which can reduce time consumption significantly.Table II lists the computing time of the arc calculation process for the classic Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 10.Effects of noise level and DEM error on success rate for five situations of residual topographical error and nine situations of noise level.

TABLE II COMPARISON OF TIME CONSUMPTION
method and the new method in both the synthetic dataset test the real dataset test.

C. Effects of Noise Level and DEM Error on Phase Unwrapping and Parameters Estimation
The accurate estimation of the differential residual topographical error is the first and most important step for our method.We can find from (3) that the model coefficient of Δε h is very small, so differential residual topographical error estimation is extremely susceptible to noise [28].In addition, low-resolution and imprecise DEM data may cause large residual topographical error and reduce parameter estimation accuracy.
Compared to the true values, two types of phase unwrapping success rates are counted for the arc phase: the network success rate and the interferogram success rate.Network success means that an arc phase has been unwrapped successfully in all interferograms, and interferogram success means that an arc phase has been unwrapped successfully in one interferogram.
The statistical results are given in Fig. 10.We find that the size of the residual topographical error hardly affects the success rate of phase unwrapping, implying that the new method is suitable for most external DEM datasets and scenarios.On the other hand, as the noise level increases, the success rate decreases quickly, especially when the mean noise is larger than 20°.Although arc elimination is used in the new method, too many failing-unwrapping phases will break the PS network and lead to some individual subnetworks, which reminds us that strict PS point selection should be performed in the new method and the noise level of PSs should less than 25°.
Furthermore, the parameter estimation accuracy is obtained for multiple simulations, and the results are shown in Fig. 11(a) and (b).The accuracy of the differential residual topographical error and differential deformation rate slowly changes and then greatly decreases with increasing noise level.Under a higher noise level, the more imprecise the external DEM is, the lower the accuracy of the differential residual topographical error for the arc.After we discard some arcs with the multi-image phase coherence threshold of 0.7, the effect of the noise level is obviously dominant, and the accuracy of the differential deformation rate in all situations is less than 1 mm/a.It should be noted that parameter errors still exist when the noise level is 0°because there are still some atmospheric delay phases present in the differential arc phase.

D. Effects of Image Time Intervals on Deformation Time-Series Extraction
As described in the INTRODUCTION, the new method in this article is proposed benefit from intensive time sampling of SAR images with the rapid development of SAR satellite technologies.Image time interval is an important parameter for new method, which can affect the construction of pseudophase with "zero" temporal baseline and phase wrapping detection of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.To assess the effect of image time intervals, we select part of images from synthetic datasets described in Section III-A and construct four datasets with average time interval of 12, 24, 36, and 48 days, respectively.Deformation time series in two PSs are estimated and the results are shown in Fig. 12.
We can find from Fig. 12 that compared with the estimation of linear deformation time series, the nonlinear deformation time series is more susceptible to the increase of image time interval.With the increase of image time interval, the assumption in Section II-A has become imprecise for nonlinear deformation and the number of pseudophases is reduced, which may lead to inaccurately calculation of residual topographical error and the reduce of multi-image phase coherence (γ) for arc phase.Then, the network is changed after discarding some arcs with low γ and then estimation accuracy of time-series deformation is reduced.In order to improve the accuracy and reliability of the new method, it is suggested that the average time interval of SAR images should be within 30 days for nonlinear deformation monitoring.

E. Effects of Inaccurately Residual Topographical Error Estimation on Deformation Time-Series Extraction
It should be noted that the differential residual topographical error estimation in Section II-A is an approximate value because the assumption is not rigorous and the residual deformation phase may still exist in pseudophase.We can find from Section IV-A and C that, after we set multi-image phase coherence threshold of 0.7 for are calculation, the parameter error of the differential residual topographical error is small (small than 5 m).Cause of the small model coefficient of Δε h in (3), the inaccurately residual topographical error estimation is not enough to cause phase ambiguities.Phase ambiguities caused by residual topographical error can be calculated precisely.
Furthermore, in Section II-B, not only the phase ambiguities caused by residual topographical error was subtracted, but the residual topographical phase.Therefore, the inaccurately estimation of residual topographical dose have impact on the deformation time-series estimation.For the common SAR satellites, such as PALSAR2, Sentinel-1, and TerraSAR-X, the orbit  positioning was controlled steadily.According to the basic geometric parameters [29], [30], [31], the effect of inaccurately residual topographical error estimation (per meter) on deformation time series can be calculated and the results are shown in Table III.We can find from Table III that the effect is slight.As the perpendicular baseline of differential interferograms are random in positive and negative numbers, the inaccurately estimation of residual topographical error cannot cause a trend error for deformation time-series estimation which can be weaken by a low-pass filtering processing.

V. CONCLUSION
In this article, we proposed a novel phase separating, unwrapping, and nonparametric deformation time series extracting method for PS-InSAR.Unlike classic methods, the PS arc phase is divided into differential residual topographical phase and differential deformation phase by phase differentiation in time dimension so that we can unwrap them separately.We use the zero-temporal baseline combination method and single parameter grid-search method to calculate phase ambiguities caused by residual topographical error, and direct detection and correction are used to unwrap the differential deformation phase.Hence, the arc phase can be unwrapped successfully, and then deformation parameters can be calculated for each PS point.Furthermore, the deformation time series can be extracted by integrating the arc differential deformation phase in space using the least squares network adjustment method.The new method is tested by using both simulated datasets and real TerraSAR-X datasets, and the experimental results show that the new method can conveniently achieve high-precision phase unwrapping for PS-InSAR and calculated deformation time series in a nonparameterized form.The RMSE of the topographical residual results and linear deformation rate results of the PS points are 3.66 m and 0.43 mm/year, respectively.Compared with the classic methods, our new method provides insight for using timedimensional information to improve PS-InSAR.The accuracy of topographical residual estimation is comparable with that of classic PS-InSAR, but there is a 34.38% improvement in the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.deformation rate accuracy, as well as 20-40 times improvement of computational efficiency.
There are three major advantages of the proposed method as follows.
1) Time-dimensional phase differentiation is added to our method, which gives us more constraints for phase unwrapping and improves the precision of deformation products.2) Classic methods divide the deformation time series into modeled parts and nonmodeled parts, and a temporalspatial filter is used to extract nonmodeled parts, which likely causes deformation loss.Compared with classic methods, the deformation time series is obtained by a nonparametric process without the assumption of a deformation model and can be restored more completely in our method.
3) The calculation efficiency can be greatly improved in our method.However, it should be noted that the new method is more suitable for SAR datasets with intensive time sampling so that fewer errors would be introduced during the phase difference in the time dimension.Moreover, the noise level can seriously affect the success rate of phase unwrapping; hence, we should select PS points carefully when using the new method.It should also be pointed out that it is difficult to give a quantitative analysis for the maximum deformation that the new method can handle because the threshold is affected by a combination of multiple factors, such as wavelength of SAR sensor, time interval for SAR data sampling, coherence of the observed signal, and density of PSs.But it is definitely that some abrupt deformations, such as earthquakes or ongoing landslides, may cause large displacement and result in misjudgment when we separate and unwrap of deformation phases based on timedomain differentiation, so like traditional PS-InSAR technique, our method is not suitable for abrupt deformations.

Fig. 2 .
Fig. 2. Diagrammatic map of deformation phase wrapping.Dramatic phase changes occur due to phase wrapping (A) or long-time intervals (B).

,
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 3 .
Fig. 3. Workflow comparison between classic PS-InSAR and new methods.The common processes are labeled with black dotted-line boxes.Classic processes are labeled with blue dotted boxes, and new processes are labeled with red dotted boxes.Five portions are labeled with gray symbols.

Fig. 4 .
Fig. 4. (a) Simulated linear deformation velocities for all PSs by the peaks function.(b) Simulated amplitude of cyclic deformation for all PSs.(c) Simulated residual topographical error for all PSs.

Fig. 5 .
Fig. 5. (a) Calculated linear deformation velocities for all PSs.(b) Calculated residual topographical error for all PSs.(c), (d) Differences between calculated results and simulated values of linear deformation rate and residual topographical error.

Fig. 6 .
Fig. 6.Deformation time series of P1, P2, and P3.Blue dots show the simulated results, red triangles show the results of new method, and green rectangle shows the results of classical method.

Fig. 8 .
Fig. 8. (a), (b) Estimated deformation rate by using classic PS-InSAR and the new method.(c) Difference between the above two methods.(d) Intercomparison results of the above two methods.

Fig. 9 .
Fig. 9. Bias histograms of differential residual topographical error (a) and differential deformation rates (b) for arcs and residual topographical (c) and deformation rates (d) for PS points compared with true values.Red indicates production by the calculated method, and blue indicates the new method.

Fig. 11 .
Fig. 11.Effects of noise level and DEM error on parameter accuracy for five situations of residual topographical error and nine situations of noise level.(a), (b) are the results of all arcs; (c), (d) are the results of arcs with multi-image phase coherence larger than 0.7.

Fig. 12 .
Fig. 12.Comparison of deformation time-series estimation in two PSs (nonlinear deformation time series and linear deformation time series) by using four datasets with average time interval of 12, 24, 36, and 48 days.

TABLE III EFFECTS
OF INACCURATELY ESTIMATED TOPOGRAPHIC ERROR ON DEFORMATION TIME SERIES