A constrained small baseline subsets (CSBAS) InSAR technique for multiple subsets

During the Small Baseline Subsets (SBAS) InSAR processing, the interferograms are often separated into multiple independent subsets, which results inevitably in a rank deficiency problem. Singular value decomposition (SVD) or linear interpolation based least squares (LS) is generally adopted, which causes systematically biased deformation estimation. Therefore, we presents a constrained SBAS method, which takes the period of the time series deformation detected by frequency-spectrum analysis as constraints to solve the rank deficiency problem. This method is illustrated with simulated data in detail. The results of the SVD, LS and our methods are in agreement with the true value in the first subset, but biased in the second subset with the magnitudes of 17.39 cm, 15.90 cm and -0.53 cm, respectively, where our method is the best one. Lastly, the new method is successfully verified using the real SAR data over Southern California from 2003 to 2006. The averaged STD of the differences between our method and GPS observations in four stations are 5.0, 3.62, 6.31 and 5.87 mm/year, respectively, which is much better than those from SVD and LS methods. This outcome confirms the validity of the newly proposed method. ARTICLE HISTORY Received 24 June 2019 Revised 23 October 2019 Accepted 18 December 2019


Introduction
To mitigate the temporal and/or spatial decorrelation and atmospheric effects (Zebker, Rosen, & Hensley, 1997;Zebker & Villasenor, 1992), various types of advanced Interferometric Synthetic Aperture Radar (InSAR) algorithms based on multi-interferograms have been proposed. The typical three categories are Permanent Scatterers InSAR (PS-InSAR) (Ferretti, Prati, & Rocca, 2000, 2001Hooper, Segall, & Zebker, 2007;Kampes, 2006), Small Baseline Subsets (SBAS) InSAR (Berardino, Fornaro, Lanari, & Sansosti, 2002;Lanari et al., 2004a;Lauknes, Zebker, & Larsen, 2010), and SqueeSAR (Ferretti et al., 2011). PS-InSAR deals with single-master interferograms, whereas multiplemaster interferograms are formed in SBAS-InSAR by setting spatial and temporal baseline thresholds, which has enabled wide applications in earth-science studies and engineering-oriented applications such as land subsidence (e.g. Del Ventisette et al., 2013;Deng et al., 2016;Hu, Li, Ding, Zhu, & Sun, 2013;Kim, Lu, Jia, & Shum, 2015) and landslides (e.g. Tong & Schmidt, 2016;Zhao, Kang, Zhang, Lu, & Li, 2018;Zhao, Lu, Zhang, & de La Fuente, 2012;Zhao et al., 2016Zhao et al., , 2013. As for SBAS-InSAR, least squares (LS) are applied to all unwrapped interferograms when only one small baseline subset is available and the cumulative time series deformation can be readily retrieved (Berardino et al., 2002;Usai, 2003). However, if the available SAR acquisitions are distributed within different subsets, an issue of rank deficiency arises in the following cases: (i) No SAR acquisitions are available that cover a certain period. For example, the ERS-1SAR images were available from May 1992 to August 1993, whereas the Envisat ASAR images were available from August 2003 to October 2010 in Taiyuan, Shanxi Province, China. No SAR data were available between August 1993 and August 2003 (Liu, Zhao, Zhang, Yang, & Zhang, 2018). (ii) Some interferometric pairs cannot meet the given spatial and temporal baseline thresholds. For example, the orbit of ALOS PALSAR-1 satellite experienced a jump in May and June 2008, causing the spatial baseline of the interferograms to exceed the threshold; thus the isolated subsets (i.e. no overlap in the time domain) cannot be connected (e.g. Sun, Hu, Zhang, & Ding, 2016;Zhao et al., 2012). (iii) Multi-sensor data such as ERS and Envisat ASAR data are involved (Bonano, Manunta, Marsella, & Lanari, 2012;Del Ventisette et al., 2013;Pepe, Sansosti, Berardino, & Lanari, 2005). Because of the different imaging geometries, the SAR images from different sensors cannot form interferograms. (iv) Low coherence and phase unwrapping errors also cause the disconnection of interferograms in some pixels. Thus, least squares (LS) method cannot be used to derive the entire time series deformation correctly when the subsets of interferograms are large than one. Typically, the singular value decomposition (SVD) method is adopted to connect such interferograms to obtain a unique solution and to increase the deformation time series, which in essence applies the minimum-norm constraint on the corrections of deformation phases or velocities between the adjacent SAR acquisitions. However, the solution is just a mathematical meaning without any physical meaning (Berardino et al., 2002;Xu, 1997;Xu, Shimada, Fujii, & Tanaka, 2000). That is, the bias between each subset is meaningless. To overcome the weakness of the SVD method, Usai (2003) addressed this problem by interpolating the disconnected subsets, which is equivalent to imposing a tight constraint on the deformation between the independent subsets and to obtain the LS solution. But, this solution may lead to large discontinuities in the final result if the constructed "missing" interferogram does not match the reality. Hu et al. (2013) proposed a method to add manually an additional goodquality interferometric pair between the acquisitions from different subsets to obtain a unique solution. However, it seems difficult to find such high-quality interferograms if the SAR images are not dense enough. López-Quiroz, Doin, Tupin, Briole, andNicolas (2009), Deng et al. (2016) and Zhang et al. (2019) proposed some polynomial models in the time domain as constraints to connect the SAR data acquired from single or multiple platforms, where the final time series results are highly dependent on the given polynomial models. The bias can hardly be avoided if the polynomial model does not accord with the surface deformation. Samsonov and d'Oreye (2012) applied the Tikhonov regularization or low-pass filtering to overcome the rank deficiency caused by multi-sensor datasets, which also lead to the biased result due to the nature of the regularization method. Pepe, Bonano, Zhao, Yang, and Wang (2016) presented a method by joint use of multiple satellite SAR data and geotechnical models to study the surface deformation of Shanghai ocean-reclaimed lands, which was just suitable for the reclaimed deformations.
This paper proposes a novel method to solve the abovementioned rank deficiency problem caused by multiple subsets and to overcome the shortcomings of the existing methods. We take the maximum period detected from the time series deformation in each independent subset as constraints to link each independent subset, which allows us to retrieve more straightforward and physically sound long-term deformation time series. This method is an extension of the traditional SBAS method, which can be used in the deformation monitoring with periodical phenomenon.
The paper is organized as follows. After a detailed description of the constraint imposed by the SVD method in multiple subsets in Section 2, the proposed method is then explained in Section 3. Furthermore, the performance of the proposed method is evaluated using both simulated experiment and real SAR data in Section 4. Finally, conclusions on the proposed method and plans for future work are given in Section 5.

SBAS method
We consider N þ 1 SAR images acquired at ordered times t 1 ; t 2 ; Á Á Á ; t Nþ1 ½ ; thus, M interferograms are generated by setting the thresholds for temporal baseline and/or spatial baseline. The unwrapped phase δϕ j at generic pixel ðx; rÞ (x and r are the azimuth and range coordinates, respectively) of the jÀth differential interferogram can be expressed as Equation (1). For the detailed description of the parameters involved in Equation (1), we refer the readers to Berardino et al. (2002).
Furthermore, Berardino et al. (2002) considered the mean deformation velocity v k;kþ1 between the timeadjacent acquisitions as the estimated unknowns to yield a physically reasonable result. Assuming that only the deformation phase component is considered, Equation (1) can be accordingly simplified as follows (Berardino et al., 2002): where δφ disp j denotes the deformation phase of the jÀth differential interferogram; k denotes the acquisition time indexes between t A and t B , as shown in Figure 1.
Then, rewriting Equation (2) in a matrix form yields where v ¼ ½v 1;2 ; Á Á Á ; v N;Nþ1 T denotes the unknown values associated with the deformation rates of the considered pixel. The generic ðj; kÞ element of design matrix B will be Bðj; When M unwrapped interferograms are within a single subset, the design matrix B in Equation (3) is an NÀrank matrix. Then, the deformation rates between the time adjacent acquisitions will be uniquely obtained using the LS method and accordingly the cumulative time series deformation can also be achieved with an additional integration operation by assuming that the deformation at the earliest SAR acquisition date t 1 is zero. However, if L ðL ! 2Þ subsets are available, the design matrix B is rank deficiency with the rank of N À L þ 1. Consequently, ðB T BÞ is a singular matrix and infinite solutions can be obtained from Equation (3). In practice, the SVD method is simply applied to obtain the pseudoinverse of coefficient matrix ðB T BÞ with the constraint ðv T vÞ ¼ min .
Combining Equations (4) and (5), the general solution can then be obtained in this example: where b 1 and b 7 are arbitrary constant.
Moreover, the following datum equations are also imposed as the constraint to give the minimum-norm LS solution in rank-deficiency-free network adjustment (Tao, 2002(Tao, , 2009):

Constrained SBAS (CSBAS) InSAR
To solve the abovementioned rank-deficient problem caused by multiple independent subsets and to overcome the shortcomings of the SVD method, a CSBAS-InSAR method with physically meaningful constraint is proposed, which will better reveal the temporal evolution of surface deformation. The core idea of the proposed method can be reduced to two main key points: the time series deformation in each independent subset is firstly generated using least squares (LS) method; the obtained independent time series results are then connected using the period detected from the aforementioned time series deformation as constraints to obtain the optimal solution that conforms to the law of surface deformation. The following is the way to detect the period from the time series deformation in independent subsets, which is used to connect the isolated subsets: We assume that M unwrapped interferograms belong to L ðL ! 2Þ independent subsets. First, the linear deformation component and topographic residual phase caused by inaccurate DEM are simultaneously estimated using the LS method based on Equation (9) (Berardino et al., 2002): where v is the linear deformation rate and δϕ res;j represents the residual phase mainly caused by nonlinear surface deformation, atmospheric artefacts and noises. By subtracting the estimated linear deformation and topographic error from each differential interferogram, the residual phase δϕ res;j can be obtained. Then, least squares (LS) inversion of the residual phase is conducted for retrieving the corresponding time series deformation by formulating the observation Equation (3) in each independent subset. Subsequently, the Lomb-Scargle periodogram (LSP) method (Lomb, 1976) is applied to carry out the frequency-spectrum analysis of each independent time series result. The corresponding periods T m ðm ¼ 1; 2; Á Á Á ; LÞ can be accordingly derived by taking the inverse of the peak frequency in each independent subset: where f max;m represents the peak frequency of mÀth subset. Generally speaking, the periods T m ðm ¼ 1; 2; Á Á Á ; LÞ detected from different subsets are approximately equal, so the period T of the entire time series can be accordingly determined. Finally, an integer multiplication of the detected period T(e.g. T; 2T; Á Á Á ; Num Ã T), that is larger than the time gap between different subsets, is selected as the constraints to connect the independent residual time series deformation. The constraints can be expressed as follows under an assumption that the cumulative residual deformation within the detected period is zero: where Num is an integer number, v i denotes the deformation rate between the times t iþ1 and t i . It should be noted that v k and v kþmÀ1 must be in the different subset; otherwise, the imposed constraints will be of no use to solve the problem caused by the rank deficiency.
Combining the observation Equation (3) and constraint Equation (11), the design matrix B then becomes full-rank when the number of imposed constraint equations is equal to or greater than the number of independent subsets ðLÞ. Therefore, the deformation rate between successive SAR acquisitions can be simply obtained by the LS method and the unique residual time series deformation are thereby generated by integrating the estimated deformation rate in time domain. Thus, the final time series deformation can be obtained by calculating the sum of subtracted linear component and estimated residual deformation component. The main processing steps of the proposed method are represented in Figure 3.

Simulated data
To demonstrate the performance of the different solutions and obtain the physically meaningful solution from multiple independent subsets, we designed three different schemes on two simulated independent subsets (Figure 4). Similar to the reference (Lauknes et al., 2010), we only simulate the deformations and test the performance on a single point. It is assumed that the simulated InSAR observations are also affected by the decorrelation noise and atmospheric artefacts with zero mean additive Gaussian noise with standard deviation of 0.1 cm and up to 1.8 cm, respectively. Without loss of generality, the unwrapped phase is considered in each interferogram. The time series deformation estimation by SVD, LS with linear interpolation or extrapolation (LS hereafter) and CSBAS-InSAR (New hereafter) method is separately implemented.

Scheme 1: with one overlap in the time domain
In order to make the simulation more realistic, we adopted the same SAR sensor parameters as those from the real Envisat ASAR dataset of the Southern California area used in the next section. For a specific pixel, assuming no deformation in the first epoch, a sinusoidal time series deformation is simulated with the period of 350 days and sampling interval of 35 days. The amplitude of the simulated signal is set to 10 cm. The interferometric combination is shown in Figure 4(a), which includes 16 interferograms. It can be visible that the interferograms with one overlap in the time domain are included between the two independent subsets. Figure 5 shows the time series deformation and corresponding residuals when the standard deviation of atmospheric artefact is set as 1.8 cm. The results of the three methods in the first subset agree well with each other. However, comparison among the results of the SVD, LS, and the true values shows that there exist obvious biases in the second subset, which reach −1.83 cm and −3.11 cm in the SVD and LS method, respectively. In contrast to the aforementioned two methods, the new method can provide a more consistent result with the simulated true values in the entire time series deformation. The bias between the result and true value is −0.32 cm, which is much smaller than those of the SVD and LS methods.
Scheme 2: without overlap in the time domain Another scheme, two independent subsets without overlaps in the time domain, is then considered. The parameters of the simulated time series are the same as those in Scheme 1 and the interferometric combination in this scheme is shown in Figure 4(b), including 16 interferograms.
The estimated time series deformation and corresponding residuals are shown in Figure 6. Here, the results of the SVD, LS and new methods are also in agreement with the simulated true value in the first Scheme 3: with multiple overlaps in the time domain A more realistic scenario is performed, which considers multiple overlaps in the time domain between the two independent subsets. The interferometric combination is shown in Figure 4(c), including 32 interferograms. The parameters of the simulated time series are also the same as those in Scheme 1.
In this scheme, the time series deformation calculation is firstly conducted in each independent subset. Subsequently, Lomb-Scargle periodogram (LSP) method is utilized to perform frequency-spectrum analysis to derive the corresponding period of the obtained residual time series deformation for each  independent subset (Figure 7). Then, the detected period (T % 350 days) is selected as the constraint in Equation (11) to connect these two disconnected subsets.
The estimated time series deformation and the corresponding residuals are shown in Figure 8. It can also be concluded in this example that if t i belongs to the first subset, the corresponding time series deformations derived by these three methods are all consistent with the simulated true values. However, if t i belongs to the second subset, the deformations derived both from the SVD and LS methods show obvious biases with the given values, while the results derived by the new method can match well with the given values.
To better evaluate the performance of these three methods, a Monte Carlo simulation experiment with 100,000 iterations for the atmospheric artefacts with different noise levels is carried out and the averaged root-mean-square errors (RMSEs) of the residuals between the simulated and estimated time series deformation for SVD, LS and new methods are shown in Figure 9. As expected, for all the cases, the results derived from the new method have lowest RMSEs in three schemes, which demonstrates that the proposed method can improve the accuracy of the estimated time series deformation when the interferograms belong to two or more different subsets, especially for the case that there is no or just one overlap between subsets in the time domain. It is reasonable to see the averaged RMSEs increase along with the increase of the noise level of the atmospheric artefacts. Meanwhile, the experiments indicate the period is estimated correctly by LSP.
It can also be summarized as follows from these experiments (Figures 5, 6, and 8): (1) the constraint imposed by the SVD method depends on the deformation rate v k;kþ1 in which t k and t kþ1 are within different subsets (e.g. t k in the first subset; t kþ1 in the second subset); (2) the constraint imposed by the LS method depends on the selected interferogram, which is used to connect the SAR acquisitions in different subsets; (3) the constraint imposed by the new method depends on the detected period of the time series deformation,  which considers the cumulative deformation within the specified period equals to 0.

Real data
Fifteen ascending Envisat ASAR scenes, which covers Southern California from 29 October 2003 to 27 December 2006, are used in this test. Twenty-four interferograms are generated with the perpendicular baseline varying from −384 m to 375 m, which are automatically separated into two independent subsets. One arc-second Shuttle Radar Topography Mission (SRTM) DEM data with the resolution of 30 m is used to eliminate the effect of topographic phase. The preliminary orbit state vectors are then refined using the DORIS data provided by the European Space Agency (ESA). To suppress the effect of noise, a multilooking operation with 2 pixels in range direction and 10 pixels in azimuth direction is performed. Then, adaptive spectral filtering of the interferograms is carried out to further reduce the phase noise (Goldstein & Werner, 1998). Minimum cost flow (MCF) (Costantini, 1998) and Delaunay triangulation algorithms are applied to obtain the unwrapped differential interferograms, in which only the pixels with coherence higher than a given threshold are considered. Moreover, a combined model is adopted to reduce the artefacts of orbital error and atmospheric disturbance, which consist of a 2-D quadratic model for the orbital error (Shirzaei & Walter, 2011) and a linear model for the elevation-dependent error (i.e. stratified atmospheric delay) (Chaabane, Avallone, Tupin, Briole, & Maitre, 2007). Finally, similar to the previous studies, the GPS site, namely, East L.A. Science Center (ELSC) on 29 October 2003 is used as the reference point and reference time for the calculation of the time series deformation (Hu et al., 2013;Lanari, Lundgren, Manzo, & Casu, 2004b). Figure 10 shows the estimated mean deformation rate from October 2003 to December 2006 and its corresponding standard deviation (STD) by the LS inversion based on Equation (9). The low-coherence regions with the average coherence lower than 0.3 are masked out because the interferometric phases in those regions cannot be correctly unwrapped. It should be noted that the obtained deformation rate is in the line  of sight (LOS) direction. The positive values indicate that the ground moves towards the satellite sensor, and vice versa. It can be clearly seen from Figure 10 that an evident increasing rate reaches 33 mm year −1 in the southern part of the investigated area while subsidence occurs in Pomona with a rate of about 29 mm year −1 . The corresponding measurement precision ranges from 0.1 mm year −1 to 6.5 mm year −1 (Figure 10). Compared with the previous researches (Hu et al., 2013), the obtained result in this study shows the similar deformation rate and pattern, indicating that the estimated deformation rate is reasonable and reliable.
After the removal of the linear deformation component and topographic error, the residual time series deformation is calculated by the aforementioned three methods (i.e. SVD, LS and the new methods). To better evaluate these three methods, the time series deformation of four typical SCIGN sites (i.e. BLSA, LBC1, PMHS and SACY) are involved. As for LS method, a "simulated" deformation result (e.g. 20031029_20040630) is calculated by linear interpolation of the pre-existed interferogram (e.g. 20031029_20040908) to connect the two different subsets. While as for the new method, the period of the residual time series deformation derived from SVD method is firstly detected by the LSP method as T % 385 days, which then is taken as the constraint to solve the unique time series deformation. To reasonably compare the InSAR-derived results with GPS  measurements, the coherent points that lie within 200 m of the GPS station are selected and averaged as the deformation value of the GPS station. Figure 11 (a-d) show that the comparison between the results derived from the three methods and those from the GPS measurements. The averaged STD values of the differences between InSAR results from these three methods and GPS observations at points BLSA, LBC1, PMHS and SACY are listed in Table 1. It can be seen that the result obtained from the new method is in better agreement with the GPS measurements than those from SVD and LS methods, which shows the same conclusion as the simulated experiments. The reasons for the discrepancies between InSAR results and GPS measurements are possibly due to the uncompensated atmospheric artefacts, orbit ramps, and even errors in the GPS observations. The results among the SVD, LS, and new methods show fixed deviations, which are the intrinsic unsolved problem of SVD and LS.

Conclusions
A constrained SBAS InSAR method is proposed to address the problem of rank deficiency caused by multiple independent subsets. To solve the deformation datum between two independent subsets, the maximum period is added as the external constraint, which is estimated by using frequency-spectrum analysis with LSP method in each subset after removing the linear deformation, topographic error and atmospheric artefact components from the interferograms, i.e. from the residual time series deformation. Then, one or more interferograms are generated across the independent subsets to link them together under the assumption that the cumulative residual deformation within the detected period is zero. Finally, LS inversion is carried out to derive more reliable time series deformation. To evaluate the performance of the proposed method, the biases among SVD, LS and the proposed methods are quantitatively compared in the simulated and real data experiments. As for real data experiment, in-situ GPS measurements are involved to assess the accuracy of three methods. The time series deformation derived by the new method agree well with the GPS measurements than the other methods, which indicates that the new method can derive more accurate deformation than those with the SVD and LS methods. However, we currently only consider one maximum period as the constraint to solve the deformation datum between two independent datasets. Once any other nonlinear deformation like exponential, logarithmic or sigmoid functions occurred, this method cannot be applicable. In the further research, more flexible physically meaningful constraints based on the hypothesis tests will be considered.