Constructing Adaptive Deformation Models for Estimating DEM Error in SBAS-InSAR Based on Hypothesis Testing

: The Interferometric Synthetic Aperture Radar (InSAR) technique has been widely used to obtain the ground surface deformation of geohazards (e.g., mining subsidence and landslides). As one of the inherent errors in the interferometric phase, the digital elevation model (DEM) error is usually estimated with the help of an a priori deformation model. However, it is difﬁcult to determine an a priori deformation model that can ﬁt the deformation time series well, leading to possible bias in the estimation of DEM error and the deformation time series. In this paper, we propose a method that can construct an adaptive deformation model, based on a set of predeﬁned functions and the hypothesis testing theory in the framework of the small baseline subset InSAR (SBAS-InSAR) method. Since it is difﬁcult to ﬁt the deformation time series over a long time span by using only one function, the phase time series is ﬁrst divided into several groups with overlapping regions. In each group, the hypothesis testing theory is employed to adaptively select the optimal deformation model from the predeﬁned functions. The parameters of adaptive deformation models and the DEM error can be modeled with the phase time series and solved by a least square method. Simulations and real data experiments in the Pingchuan mining area, Gaunsu Province, China, demonstrate that, compared to the state-of-the-art deformation modeling strategy (e.g., the linear deformation model and the function group deformation model), the proposed method can signiﬁcantly improve the accuracy of DEM error estimation and can beneﬁt the estimation of deformation time series.

It is acknowledged that the interferometric phase includes not only the interested deformation component, but also undesirable noise components (e.g., decorrelation noise, atmospheric delay, and digital elevation model (DEM) error) [18]. Although the specific implementation process of different TS-InSAR methods is different, they should all first correct or weaken these noises before obtaining reliable deformation results. For example, decorrelation noise can be suppressed by multi-looking operation or the spatial filter [15]; stratified atmospheric delay can be mitigated by an elevation-dependent model [19,20] and turbulence atmospheric delay is generally reduced by a spatial-temporal filter or external meteorological datasets [14,21,22]; long-wavelength orbit error can be effectively alleviated by fitting second-order polynomials [23,24]. As for DEM error, it is an undesired component for deformation estimation. However, if the DEM error can be precisely estimated, it can be added back to the existing DEM data to generate a more accurate version of DEM data [25,26]. The DEM error-related interferometric phase is proportional to the spatial baseline of the interferogram [27][28][29][30]. In this case, the DEM error can be theoretically estimated by mathematical modeling. Generally, the deformation time series can be assumed to satisfy an a priori time-dependent deformation model, and the DEM error and model parameters can be simultaneously solved. In state-of-the-art deformation modeling strategies, a linear deformation model [13,31,32] is one of the most used a priori deformation models, which has been employed to obtain the deformation parameters associated with tectonic movement, ground subsidence, and so on. For specific deforming situations (e.g., permafrost and active caldera), we can employ the period model [33], cubic model [15], or even the function group deformation model [34]. Note that the accuracies of the DEM error and the final deformation time-series are both dependent on the reliability of the a priori deformation model. However, in reality it is impracticable to understand the real deformation evolution beforehand, and it is difficult to find a reliable model that can fit the long time-span deformations. Although the DEM error can also be estimated by a nonparametric estimator (e.g., independent component analysis [35]), the performance of this estimator may be degraded in some situations when only a small number of observations in space are available or high-level noise and outliers exist in observations [36].
To tackle the aforementioned problem, relative to the deformation model, we propose a method to construct adaptive deformation models for estimating DEM error based on hypothesis testing in the framework of the SBAS-InSAR method. The basis of this method is that the natural ground deformation is usually a time-dependent smooth process [37], and the increasing temporal sampling of SAR images (e.g., Sentinel-1 SAR images) will make the deformations at adjacent acquisitions more consistent with this smooth process. Under these circumstances, if the whole time span is divided into several small time spans, a set of time-dependent functions (e.g., cubic polynomials and periodic function) can be used to fit the deformation time series in each small time span. However, too many time-dependent functions easily induce the overfitting of deformation time series, resulting in possible biases in the estimation of parameters. Hypothesis testing is an inferential statistical process that uses sample data to assess the plausibility of established null hypotheses, which has been used to adaptively determine the influencing factors of the deformation model [38][39][40]. Therefore, we employed hypothesis testing here to statistically test the parameter significance of the predefined functions, and then those significant parameters were adaptively selected to construct the final deformation model used to assist the estimation of DEM error. This self-adaptive deformation model not only avoids the error induced by the simple model that cannot effectively fit the true deformation, but also reduces the possible overfitting error of the too complex function group.
To validate the superiority of the proposed method, both simulations and real data experiments over the Pingchuan mining area were conducted. The Pingchuan mining area belongs to the Jingyuan coal field, which is one of the most important coal bases in the Gansu Province, China. Coal mining began in this area in as early as the 1970s, and the ground deformations time series are significant and nonlinear, as a result of the continuous mining activities.
The rest of this paper is structured as follows: Section 2 presents the proposed method of constructing adaptive deformation models based on hypothesis testing; then, the validation of the proposed method with simulated experiments and real data experiments in the Pingchuan mining area, Gaunsu, China, is shown in Section 3. Finally, some discussions are given and conclusions drawn in Sections 4 and 5, respectively.

Methodology
Since the proposed method is based on the framework of the SBAS-InSAR method, we first introduce the basic idea of the SBAS-InSAR method in Section 2.1; then, details about the proposed method for constructing adaptive deformation models, based on hypothesis testing, are subsequentially presented in Sections 2.2-2.4.

Basic Idea of the SBAS-InSAR Method
In the SBAS-InSAR method, multi-prime interferograms with short spatial-temporal baselines are used for the estimation of deformations, which can significantly suppress the decorrelation noise in the vegetation area, compared to single-prime interferograms. Generally, the interferometric phase between two SAR images for a pixel can be expressed as [18]: where δϕ represents the observed interferometric phase, δϕ de f represents the surface deformation phase along the line-of-sight (LOS) direction during two acquisitions, δϕ topo represents the topographic residuals phase due to the DEM error, δϕ orb represents the orbit error-related phase and can be modeled by second-order polynomials, δϕ atm represents the difference in atmospheric delay between two acquisitions, and δϕ noise represents the random noise.
Assuming that the orbit error phase is mostly removed by existing methods and the atmospheric delay can be regarded as random noise in the temporal domain, the interferometric phase can be written as: where the topographic residuals phase can be modeled as [27]: where B ⊥ is the perpendicular spatial baseline between the two acquisitions, r is the slant range between the SAR antenna and ground surface, θ is the looking angle, and dz represents the DEM error. Since the observation system that includes both the deformation time series and the DEM error is underdetermined, an a priori deformation model is commonly employed in the SBAS-InSAR method before the estimation of the final deformation time series. In the SBAS-InSAR method, a linear model is generally employed to fit the deformation time series [15,41]: where λ is the radar wavelength, ∆t is the temporal baseline of the interferogram, and v is the unknown deformation model parameter.
Given that there are M interferograms obtained from N SAR images based on the spatial-temporal baseline thresholds, the DEM error and the deformation model parameters can be solved with Equation (2) in a least-squares sense. Then, the DEM error phase components are removed from the M interferograms, after which the interferogram residuals are used to estimate the N deformation time series by taking the first acquisition as the reference.
In this paper, the innovation lies in the adaptive construction of the deformation model (i.e., Equation (4)). To this aim, phase time series were firstly estimated from the multi-prime unwrapped interferograms, which is more appropriate for the conduction of hypothesis testing compared to interferograms (Section 2.2). Secondly, hypothesis testing was employed to adaptively choose the deformation models from a set of time-dependent Remote Sens. 2021, 13, 2006 4 of 19 functions (Section 2.3). Finally, the adaptive deformation model parameters and the DEM error were simultaneously estimated based on a least-squares method (Section 2.4).

Derivation of Phase Time Series from Multi-Prime Interferograms
Before constructing the adaptive deformation model based on hypothesis testing, multi-prime unwrapped interferograms are used to obtain the phase time series, which can decrease the sensitivity of parameter estimation to the interferogram network [30]. The relationship between the M unwrapped phases δϕ = [δϕ 1 , δϕ 2 , . . . , δϕ M ] T and the N phase time series ϕ = [ϕ 1 , ϕ 2 , . . . , ϕ N ] T can be modeled as where B denotes the design matrix between the interferograms and the phase time series. If the M interferograms belong to a single baseline subset, Equation (5) can be directly solved by the least-squares method (i.e., Equation (6)).
However, if not, the singular value decomposition (SVD) or the iteratively reweighted least-squares (IRLS) methods can used to estimate the unknowns.
Since the characteristics of real ground deformation are complex and changeable, it is difficult to fit the deformations by only one model for long time series. To choose a model that is closer to the real deformation, the observations are divided into different groups according to a fixed time span. Moreover, the adjacent groups have 20% overlapping observations to allow different groups to constrain one another. Supposing there are J groups, the number of observations in each group is n j , and the number of overlapping observations is l. Generally, the time series in one year are taken as a group, and details about the determination of the time span for one group are presented in Section 4.1.

Construction of the Adaptive Deformation Model Based on Hypothesis Testing
In this paper, the functions in Equation (7) are used as the original deformation models: where k j is the constant term, v j is the mean velocity, a j is the acceleration, ∆a j is the acceleration variation, s j and c j are the coefficients of the sine and cosine functions, T is defined as 365 days, ϕ j = [ϕ  (7) can be written as: where k j , v j , a j , ∆a j , s j , c j are the unknown deformation model parameters; Moreover, two variables are defined for better expression of the following process.
Then, hypothesis testing is conducted for each group to test the significance of the coefficient of time-dependent variables in Equation (8) [38,39]. For generality, the process of hypothesis testing for the jth group is illuminated as follows.
First, the significance of the original model (i.e., Equation (8)) should be tested. The null hypothesis is: If H 0 can be accepted, the original deformation model is considered non-significant, and vice versa. To test the significance of the deformation model, the following statistical variable is established [38,39]: where F(n1, n2) is the F distribution whose degree of freedom is n1 and n2, , ϕ j is the mean value of ϕ j , andφ j i is the reestimated value of ϕ j i based on the unknown estimation. The deformation model is significant when F > F α p j 0 , n j − p j 0 − 1 , where α is the significance level, and can usually be set as 0.01.
However, it may happen that the model is insignificant due to the existence of severe noise (e.g., atmospheric delay). In this case, the following hypothesis testing process would not be conducted, and the null hypothesis H 0 would be adopted, which means that a constant model is used to fit the deformations.
Second, if the deformation model is significant, the significance of each deformation model parameter, except the constant term, will be tested. The null hypothesis is: If H 0 is accepted, x j u is considered insignificant and the parameter x j u should be removed from the model, and vice versa. Similarly, the following statistical variable is established [38,39]: , and t(n) represents the t distribution whose degree of free- dom is n. The corresponding parameter is significant when |T u | > t α 2 n j − p j − 1 , where α is the significance level, and can usually be set as 0.01.
After hypothesis testing of all model parameters, a new deformation model can be obtained by removing the insignificant parameters from Equation (8): where B j and X j can be obtained by removing the corresponding parts of the insignificant parameters from B j 0 and X j 0 , respectively. The deformation model of each group of phase time series can be adaptively constructed based on the aforementioned process (i.e., Equations (8)- (14)), and are used in the following subsection to simultaneously estimate the deformation model parameters and the DEM error.

Estimation of the Deformation Model Parameters and DEM Error
The DEM error-related phase is usually small in phase time series compared to the deformation signal, and the natural ground deformation is usually a time-dependent smooth process. In this paper, the phase differences of the adjacent acquisitions are taken as the observations to increase the proportion of the DEM error-related phase among the observations. Since the phase time series are divided into several overlapping groups, there are two kinds of observation equations that can be established, including (1) the observation equations for each group related to the DEM error and their individual model parameters of each group, and (2) the observation equations in the overlapping time span between two adjacent groups. The deformation time series should be equal for the two sets of model parameters of the two adjacent groups, based on which the latter kind of observation equation can be established.
As for the first kind of observation equation, by referring to Equation (2), we can model the jth group of observations as: Based on the adaptive deformation model constructed in the last subsection and the relationship between the DEM error dz and the related phase ∆ϕ j topo , Equation (15) can be rewritten as: where: . .  (7)) representing the relationship between the deformation model parameters and the ith observation of the jth T is the deformation model parameter, B ⊥ t j i is the perpendicular baseline of the SAR image at the ith acquisition of the jth group t j i relative to the initial reference SAR image, with i = 1, 2, . . . , n j , dz is the DEM error, λ is the wavelength of the SAR signal, r is the slant distance between the sensor and the ground point, and θ is the incident angle. In summary, the number of this first observation equation is In the overlapping region between two adjacent groups, the following pseudoobservation equations can be established based on the constraint that the deformation time series should be equal for the two sets of model parameters: where: and l is the number of overlapping observations between two adjacent groups. In summary, the number of this second observation equation is (l − 1) × (J − 1). Combining Equations (16) and (17) the following system can be obtained, from which the deformation model parameters and the DEM error can be simultaneously estimated based on the least-squares method.
After the correction of the DEM error from the phase time series, deformation time series can be obtained. In addition, if the deformation residuals are still severely affected by atmospheric delays, a step to mitigate the atmospheric delays is necessary [14,22]. The overall processing chain of the proposed method for constructing an adaptive deformation model based on hypothesis testing is shown in Figure 1.

Simulated Experiments
In the simulation, 159 interferograms were generated from 52 SLCs real spatial-temporal baselines of the Sentinel-1A SAR data acquired ov mining area, Gansu Province, China, which were used in the real data spatial pattern of the deformation field (200 × 200 pixels) was simulated tional volumetric change of the subsurface fluid [42] (see Figure 2). Th mation patterns were simulated by four kinds of functions, i.e., linear, and complex functions, and the simulated deformation time series at th Figure 2 are shown in Figure 3. The simulated deformation time series be obtained by multiplying the pixel value ( Figure 2) and the temporal terns ( Figure 3). The DEM error is the difference between Shuttle Radar sion (SRTM) DEM and TanDEM-X DEM (see Figure 4). The atmospheric ulated using a fractal surface with the fractal dimension being 2.2 and the being 1.0 rad for each SAR image [43].

Simulated Experiments
In the simulation, 159 interferograms were generated from 52 SLCs according to the real spatial-temporal baselines of the Sentinel-1A SAR data acquired over the Pingchuan mining area, Gansu Province, China, which were used in the real data experiment. The spatial pattern of the deformation field (200 × 200 pixels) was simulated based on the fractional volumetric change of the subsurface fluid [42] (see Figure 2). The temporal deformation patterns were simulated by four kinds of functions, i.e., linear, periodic, logistic, and complex functions, and the simulated deformation time series at the central point of Figure 2 are shown in Figure 3. The simulated deformation time series at each point can be obtained by multiplying the pixel value ( Figure 2) and the temporal deformation patterns (Figure 3). The DEM error is the difference between Shuttle Radar Topography Mission (SRTM) DEM and TanDEM-X DEM (see Figure 4). The atmospheric delays were simulated using a fractal surface with the fractal dimension being 2.2 and the maximum value being 1.0 rad for each SAR image [43]. The decorrelation noise was modeled as a zero-mean Gaussian random process with the standard deviation σ = (1 − γ 2 )/2/γ 2 , where γ = exp(−∆t/T) is the coherence simulation based on the temporal baseline ∆t and a time constant T (T = 60 days in this paper to let γ > 0.5) [37]. Figure 5 shows these simulated components of one interferogram.         As shown in Figure 6, the DEM errors estimated by Models 1 and 2 are serious affected by deformation signals, which is expected, since the deformation model used these two methods cannot well describe the deformation evolution. On the contrary, t proposed method can achieve a higher accuracy of DEM error estimation due to the ada tive deformation model and the grouping strategy. For linear deformations (Figure 6aall three methods can obtain very similar DEM error results, since the linear deformati can be well fitted by the linear model, the function group model, and the proposed ada tive deformation model. For periodic (Figure 6d-f) and logistic deformations (Figure 6 i), the proposed method can obtain more reasonable DEM error results compared to t other two methods. For complex deformations (Figure 6j-l), although it is difficult to the deformation evolution with these common functions, the proposed method still pe forms much better than the other two methods. It should be noted that, for the compl model, the location of the discontinuity concerns the accuracy of the results. If the disco tinuity is located in the overlapping area of two adjacent groups, the fitting reliability both of these groups would be influenced. Otherwise, only one group would be infl enced. In real situations, the moment of deformation jump can usually be known befo the data processing. In this case, the grouping strategy can be modified based on defo mation jump moment to suppress this influence on the result. For comparison, the traditional SBAS-InSAR method with the linear deformation model (Model 1), the traditional SBAS-InSAR method with function group (i.e., Equation (7)) (Model 2), and the proposed method (Model 3) were used to estimate the DEM error and the final deformation time series in the four temporal patterns of deformations. Note that the Model 1 and Model 2 are applied to the whole time series. The accuracy of the estimation is quantitatively described by the root-mean-square error (RMSE), which can be calculated by: whereê i and e i are the parameter estimation and simulated true value, respectively, and NU M is the number of samples. Note that NUM equals the number of pixels in the calculation of the RMSE of DEM error, and equals the production of the number of pixels and the number of SAR images in the calculation the RMSE of deformation. As shown in Figure 6, the DEM errors estimated by Models 1 and 2 are seriously affected by deformation signals, which is expected, since the deformation model used by these two methods cannot well describe the deformation evolution. On the contrary, the proposed method can achieve a higher accuracy of DEM error estimation due to the adaptive deformation model and the grouping strategy. For linear deformations (Figure 6a-c), all three methods can obtain very similar DEM error results, since the linear deformation can be well fitted by the linear model, the function group model, and the proposed adaptive deformation model. For periodic (Figure 6d-f) and logistic deformations (Figure 6g-i), the proposed method can obtain more reasonable DEM error results compared to the other two methods. For complex deformations (Figure 6j-l), although it is difficult to fit the deformation evolution with these common functions, the proposed method still performs much better than the other two methods. It should be noted that, for the complex model, the location of the discontinuity concerns the accuracy of the results. If the discontinuity is located in the overlapping area of two adjacent groups, the fitting reliability of both of these groups would be influenced. Otherwise, only one group would be influenced. In real situations, the moment of deformation jump can usually be known before the data processing. In this case, the grouping strategy can be modified based on deformation jump moment to suppress this influence on the result. Furthermore, the accuracy of the final deformation time series was also in and the residual histograms of the deformation time series with respect to diffe ods are shown in Figure 7. As can be seen, the accuracy of the deformation tim consistent with the accuracy of the DEM error estimation, i.e., the proposed m achieve better deformation results than the other two methods. This indicates formation model involved in the SBAS-InSAR process has a nonnegligible eff deformation time series. Furthermore, the accuracy of the final deformation time series was also investigated, and the residual histograms of the deformation time series with respect to different methods are shown in Figure 7. As can be seen, the accuracy of the deformation time series is consistent with the accuracy of the DEM error estimation, i.e., the proposed method can achieve better deformation results than the other two methods. This indicates that the deformation model involved in the SBAS-InSAR process has a nonnegligible effect on final deformation time series.
For determining the adaptive deformation model, the phase time series without the correction of atmospheric delay were used as the observation. In this case, the existence of atmospheric delay was not conducive to the construction of adaptive deformation model. Therefore, we conducted a series of simulated experiments with the maximum atmospheric delay value ranging from 0.0 to 2.0 rad to investigate the effect of atmospheric delay on the estimation of DEM error. As shown in Figure 8, the RMSEs of the DEM error estimation increase with an increase in the atmospheric delay magnitude for all kinds of deformations and all three methods. Nevertheless, the proposed adaptive deformation model achieved higher accuracy for different magnitudes of atmospheric delay compared to the other two methods, demonstrating the superiority of the proposed adaptive deformation model. determining the adaptive deformation model, the phase time series without the correction of atmospheric delay were used as the observation. In this case, the existence of atmospheric delay was not conducive to the construction of adaptive deformation model. Therefore, we conducted a series of simulated experiments with the maximum atmospheric delay value ranging from 0.0 to 2.0 rad to investigate the effect of atmospheric delay on the estimation of DEM error. As shown in Figure 8, the RMSEs of the DEM error estimation increase with an increase in the atmospheric delay magnitude for all kinds of deformations and all three methods. Nevertheless, the proposed adaptive deformation model achieved higher accuracy for different magnitudes of atmospheric delay compared to the other two methods, demonstrating the superiority of the proposed adaptive deformation model.

Real Data Experiments in the Pingchuan Mining Area
The Pingchuan mining area belongs to the Jingyuan coal field, which is one of the important coal bases in Gansu Province, China (see Figure 9). The Pingchuan mining area is rich in coal reserves, and coal mining began in this area as early as the 1970s [44]. With the continuous mining of underground coal, the ground surface suffers serious defor-

Real Data Experiments in the Pingchuan Mining Area
The Pingchuan mining area belongs to the Jingyuan coal field, which is one of the important coal bases in Gansu Province, China (see Figure 9). The Pingchuan mining area is rich in coal reserves, and coal mining began in this area as early as the 1970s [44]. With the continuous mining of underground coal, the ground surface suffers serious deformation, which has a negative impact on the safety of people's lives and property. The SBAS-InSAR method can be used to monitor such deformation, therefore providing a reliable database for policy decisions. However, the classical SBAS-InSAR method is vulnerable to the unreliable a priori deformation model, which is especially important for the case of highly nonlinear deformation time series in the mining area. Therefore, in this paper, the proposed method was employed to obtain more reliable estimations of the DEM error and deformation for the Pingchuan mining area.
Fifty-two Sentinel-1A SAR images were acquired over the Pingchuan mining area from 27 March 2017 to 15 July 2019 (see Figure 9a), and 158 interferograms were generated with the maximum temporal and spatial baselines of 36 days and 200 m, respectively (Figure 9b). To ensure all interferograms belong to a single set, we manually connected the SAR images acquired on 5 October 2017 and 22 November 2017 to generate the 159th interferogram (see the red line in Figure 9b). During the InSAR processing, the topographic phase was removed based on the Shuttle Radar Topography Mission (SRTM) 1 arc-second DEM [45].

Real Data Experiments in the Pingchuan Mining Area
The Pingchuan mining area belongs to the Jingyuan coal field, which is one of the important coal bases in Gansu Province, China (see Figure 9). The Pingchuan mining area is rich in coal reserves, and coal mining began in this area as early as the 1970s [44]. With the continuous mining of underground coal, the ground surface suffers serious deformation, which has a negative impact on the safety of people's lives and property. The SBAS-InSAR method can be used to monitor such deformation, therefore providing a reliable database for policy decisions. However, the classical SBAS-InSAR method is vulnerable to the unreliable a priori deformation model, which is especially important for the case of highly nonlinear deformation time series in the mining area. Therefore, in this paper, the proposed method was employed to obtain more reliable estimations of the DEM error and deformation for the Pingchuan mining area.   (Figure 10d), it is easy to see that the spatial pattern of the DEM error estimation by Models 1 and 2 is somewhat similar to the spatial pattern of the deformation rate, indicating that the deformation and the DEM error are not well distinguished in these two methods. On the contrary, the DEM error estimated by Model 3 seems more reasonable and shows less correlation with the deformation. Furthermore, the deformation time series of eight points (i.e., P1-P8 in Figure 10) are presented in Figure 11. The deformations at P1, P3, P5, P6, and P7 show an obvious nonlinear pattern. Compared to Figure 10a,b, it can be found that the DEM error signals at these points are all very large, indicating that these DEM error signals may be caused by the inaccurate deformation models used in Models 1 and 2. As for the points P2, P4, and P8, at which the deformation time series are not as nonlinear as the deformation at P1, P3, P5, P6, and P7, the magnitude of the DEM error estimation is not so correlated with the deformation rate.
Compared to the deformation rate map (Figure 10d), it is easy to see that the spatial pattern of the DEM error estimation by Models 1 and 2 is somewhat similar to the spatial pattern of the deformation rate, indicating that the deformation and the DEM error are not well distinguished in these two methods. On the contrary, the DEM error estimated by Model 3 seems more reasonable and shows less correlation with the deformation. Furthermore, the deformation time series of eight points (i.e., P1-P8 in Figure 10) are presented in Figure 11. The deformations at P1, P3, P5, P6, and P7 show an obvious nonlinear pattern. Compared to Figure 10a,b, it can be found that the DEM error signals at these points are all very large, indicating that these DEM error signals may be caused by the inaccurate deformation models used in Models 1 and 2. As for the points P2, P4, and P8, at which the deformation time series are not as nonlinear as the deformation at P1, P3, P5, P6, and P7, the magnitude of the DEM error estimation is not so correlated with the deformation rate. P8 are the selected points.The deformation time series at P3 were taken as an example to show the difference between the different methods (see Figure 11b). As can be seen, the differences in the deformation time series obtained by the different methods are much smaller compared to the long-term accumulated deformation. Nevertheless, the deformation varies as much as several millimeters in some moments between the different  (Figure 11c). This magnitude of deformation cannot be negligible when millimeter-level precision is expected in the monitoring of ground deformation. Figure 12 shows the deformation time series retrieved from the proposed method, in which only one deformation map per month is illustrated. The maximum cumulate deformation along the LOS direction, which occurred around the vicinity of P6, was approximately -88.5 cm, where the negative value indicates the ground moving away from the satellite.  The deformation time series at P3 were taken as an example to show the difference between the different methods (see Figure 11b). As can be seen, the differences in the deformation time series obtained by the different methods are much smaller compared to the long-term accumulated deformation. Nevertheless, the deformation varies as much as several millimeters in some moments between the different methods ( Figure 11c). This magnitude of deformation cannot be negligible when millimeter-level precision is expected in the monitoring of ground deformation. Figure 12 shows the deformation time series retrieved from the proposed method, in which only one deformation map per month is illustrated. The maximum cumulate deformation along the LOS direction, which occurred around the vicinity of P6, was approximately −88.5 cm, where the negative value indicates the ground moving away from the satellite.

Determining the Time Span of Each Group in the Grouping Process
Since one of the key steps of the proposed method is to group observations, we verified the rationality of this grouping strategy based on the simulation of four kinds of mining subsidence deformations (see Figure 13a) and simultaneously determined the time span of each group in the grouping process. The mining subsidence was simulated based on the logistic function [46]:

Determining the Time Span of Each Group in the Grouping Process
Since one of the key steps of the proposed method is to group observations, we verified the rationality of this grouping strategy based on the simulation of four kinds of mining subsidence deformations (see Figure 13a) and simultaneously determined the time span of each group in the grouping process. The mining subsidence was simulated based on the logistic function [46]: where d(t) denotes the cumulative subsidence at time t, d 0 represents the maximum subsidence value, a and b are the shape parameters of the logistic function. Here, d 0 = −0.8 m, and the values of (a, b) are (2400, 0.015), (2400, 0.030), (1200, 0.015), and (1200, 0.030), respectively. Similar to the Sentinel-1 SAR satellite, we simulated 70 acquisitions with a 12-day sampling interval across two and half years, and the time span of the observations in each group varied from six months to two and a half years. As shown in Figure 13b, the RMSEs of the fitting deformations by grouping increase with the increase in the observation time span in a group for the four kinds of deformations, indicating that the shorter the time span in each group (i.e., the more groups), the higher the fitting accuracy of the deformation model. However, the short time span decreases the number of independent time series observations, and would increase the risk of model overfitting. Therefore, a one-year time span in a group was preferred in this paper to obtain a tradeoff between the accuracy of deformation and the possibility of overfitting, and this grouping strategy was adopted in both the simulated and the real data experiments. (1200, 0.030), respectively. Similar to the Sentinel-1 SAR satellite, we simulated 70 acquisitions with a 12-day sampling interval across two and half years, and the time span of the observations in each group varied from six months to two and a half years. As shown in Figure 13b, the RMSEs of the fitting deformations by grouping increase with the increase in the observation time span in a group for the four kinds of deformations, indicating that the shorter the time span in each group (i.e., the more groups), the higher the fitting accuracy of the deformation model. However, the short time span decreases the number of independent time series observations, and would increase the risk of model overfitting. Therefore, a one-year time span in a group was preferred in this paper to obtain a tradeoff between the accuracy of deformation and the possibility of overfitting, and this grouping strategy was adopted in both the simulated and the real data experiments.

Decreasing Correlation between the Deformation Rate and the DEM Error Estimation Based on the Proposed Method
In order to further demonstrate the superiority of the proposed method, we calculated the correlation coefficients between the deformation rate and the DEM error estimations of the different methods at the vicinity of P1-P8 (see Table 1).

Decreasing Correlation between the Deformation Rate and the DEM Error Estimation Based on the Proposed Method
In order to further demonstrate the superiority of the proposed method, we calculated the correlation coefficients between the deformation rate and the DEM error estimations of the different methods at the vicinity of P1-P8 (see Table 1). As shown in Table 1, the DEM error estimations at the points whose deformation patterns are roughly linear (P2, P4, and P8) show low correlation with the deformation rate for all three deformation models. However, at those points, with obviously nonlinear deformations (P1, P3, P5, P6, and P7), the DEM error estimation shows a strong correlation with the deformation rate for Model 1, with correlation coefficients greater than 0.5; this correlation is decreased for Model 2 compared to Model 1, since the function group in Model 2 can better fit nonlinear deformation than the linear function in Model 1. For Model 3, the DEM error estimations show the slightest correlation with the deformation rate compared to the results of Models 1 and 2, indicating that the proposed method can construct a more reliable deformation model and then obtain more reasonable DEM error estimation.

Conclusions
DEM error is one of the components in the InSAR phase and is usually estimated with the help of an a priori deformation model. However, since the temporal deformation evolution is usually unknown, it is difficult to use a definite function to fit temporal deformation, which is not conducive to the estimation of DEM error and deformation time series. In this paper, we proposed a method to construct adaptive deformation models in the SBAS-InSAR framework based on the hypothesis testing theory. In particular, the phase time series were first divided into several groups with overlapping regions. In each group, we used a set of predefined functions to fit the deformation time series. In order to prevent the occurrence of overfitting, the hypothesis testing theory was introduced to adaptively select the optimal deformation model from the predefined functions, after which the parameters of the adaptive deformation models in each group and the DEM error were simultaneously solved by a least-squares method. Both simulated and real data experiments in the Pingchuan mining area were conducted. The results show that the proposed method can effectively improve the accuracy of DEM error estimation and can also benefit the accuracy of deformation time series.
It should be noted that the proposed method can also be embedded in other state-ofthe-art TS-InSAR techniques to construct more suitable deformation models and to obtain more reliable DEM error, as well as deformation. Moreover, the proposed method can be further improved by adjusting the original deformation model functions (i.e., Equation (7)) based on a priori information about the deformation evolution process. For example, a time-related jump function (e.g., the Heaviside function [47]) can be used if there is an abrupt event.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.