1. Introduction
Synthetic Aperture Radar Interference Strategy (InSAR) has been used to detect geophysical deformation signals such as landslides [
1], volcanoes [
2], and land subsidence [
3], as well as other geological disasters. For InSAR data, the accuracy of its products strongly depends on the accuracy of the interferometric phase [
4]. However, owing to several decorrelation factors, the SAR interferometric phase is often degraded. Many advanced algorithms have been developed to improve the reliability of InSAR products, such as the persistent scatterers algorithm in single-primary interferograms [
5,
6,
7] and the distributed scatterers (DS) algorithm in small baseline (SB) interferograms [
8,
9]. DS targets are characterized by varying coherence and phase stability [
10], especially when the interferometric SAR data pairs have long spatial and temporal baselines.
In multi-temporal InSAR technology, two strategies are usually employed to improve the density of monitoring points. The first is to use spatiotemporal filter technology to improve the phase stability [
11], and the second is to select appropriate interferometric SAR pairs to mitigate the effects of phase decorrelation [
8].
Within spatiotemporal filter strategies, the averaging of statistically homogeneous pixels (SHPs) can improve the phase stability of DS targets, in which SHPs can be determined by SAR amplitude information [
12]. Further,
N × (
N − 1)/2 non-redundant interferograms (after averaging SHPs) can be applied to achieve temporal phase optimization using the phase triangle inconsistency information [
9,
13,
14,
15]. These methods are more effective in the case where the surface features are relatively distinguishable or the phase gradient is small. However, for small-scale or larger gradient deformation targets, such as landslides and mining collapses, such advanced InSAR approaches have two disadvantages.
First, in the deformation region with dense interferometric fringes, averaging SHPs obscures the phase continuity in the spatial domain, especially for large-gradient deformation monitoring, and in the case of using low- and medium-resolution SAR data, such as Sentinel-1A/B images [
16]. Second, interferograms with longer temporal intervals and larger deformation are usually unusable, even with the SHPs filter and temporal phase optimization methods. Additionally, dense fringes of interferograms or sparse spatial pixels make phase unwrap** challenging. Pepe et al. proposed SB interferograms to perform temporal phase optimization by nonlinear optimization in order to reconstruct wrapped phase time series [
17]. Meanwhile, a sequential temporal filter method was proposed by the subspace clustering of block diagonal technology [
15]. Today, SAR satellites provide large amounts of SAR images at unprecedented temporal resolutions. Further, stepwise temporal phase optimization can promote InSAR dynamical processing and deformation monitoring, which is greatly beneficial to disaster mitigation. Therefore, a new spatiotemporal filtering strategy is required.
From appropriate interferometric SAR data pairs strategies, selecting appropriate SB interferograms can mitigate the effects of decorrelation in what is known as small baseline subset (SBAS)-InSAR processing. Some approaches for selecting optimal SB interferograms have been previously explored, including the minimum spanning tree [
18], simulated annealing [
17], semi-automatic [
19], and error bound algorithms [
20], as well as the graph theory and variance-covariance matrix [
21]. In fact, each pixel has different appropriate SB interferograms. Therefore, the intermittent SBAS (ISBAS) technique [
22] and variable length deformation time series [
23] have been proposed. However, the accuracy of ISBAS-derived deformation products is lower than those derived by the original SBAS method [
24]. Previous studies have explored the accuracy of InSAR deformation parameters with external observations, such as leveling and continuous GPS measurement data, residual noise, network configuration, and data gaps [
25], as well as using the relative error upper bound approach according to the error perturbation theory [
26]. The presence of adaptive SB interferograms in ISBAS processing directly affects its ability to resist errors. Therefore, a quality assessment method for the ISBAS method is required.
More recently, in the big SAR era, sequential InSAR deformation time series theories have been proposed, such as the complex sequential least squares (SLS) method for pixel offset SBAS [
27], as well as the SLS [
28], robust SLS [
29], modified SLS [
30], and Kalman filter [
30,
31,
32] methods for InSAR time series. With the increase in SAR data, ground targets may be subject to significant change over the long term, meaning that many monitoring point targets could become decorrelated. Therefore, it is necessary to extend the sequential InSAR methods from the SBAS frame to the ISBAS frame.
How to use the big SAR data to dynamically promote hazard monitoring? To improve the monitoring capability of InSAR in serious decorrelation, we propose a sequential DS-ISBAS InSAR deformation parameter dynamic estimation and quality evaluation for serious decorrelation scenarios algorithm. The proposed method can obtain reliable monitoring point targets in serious decorrelation scenarios and dynamically estimate deformation time series.
Firstly, to improve computational efficiency and promote InSAR dynamic data processing, truncated SB interferograms are used to dynamically perform step-by-step temporal phase optimization after multi-look and Goldstein filtering. Secondly, to obtain an increased number of monitoring point targets, a sequential ISBAS algorithm is proposed to dynamically update the deformation time series, in which the variance-covariance matrix in the Gauss-Markov model is used to evaluate the quality of the deformation parameters pixel-by-pixel. The proposed approach is then applied to monitor the post-failure deformation of the Baige landslide at the **sha River in China, using a stack of 72 Sentinel-1 SAR data images collected from 8 November 2018 to 2 April 2021.
2. Methodology
Figure 1 shows the flow chart of the proposed approach, including four steps. First, we initialize the interferograms, deformation parameters, and the quality file. Second, when we obtain a new SAR image, new filtered SB interferograms are added. Third, we execute stepwise SB interferograms filtering. Finally, we update the ISBAS deformation time series and its quality.
2.1. Review Phase Triangulation Algorithms
Any filtering with each interferogram method can be expected to lead to the disclosure phase within the temporal phase triangle. The mathematical framework of phase triangle algorithms can be expressed as in Equation (1) [
33]:
where
is the filtered interferometric phase, and
is the optimal phase,
m and
n correspond to SAR acquisition dates, and
w is the weight. The difference between these algorithms (temporal filter or phase triangulation algorithms) is that they rely on different weights
w [
33], such as the weight
in the SqueeSAR algorithm [
9];
in the EVD algorithm [
14], where
u is the eigenvector; and
in the nonlinear optimization algorithm [
17]; as well as the coherence weight
[
33].
Phase linking (PL) is an iterative strategy to perform temporal phase optimization [
13] and can be expressed as Equation (2):
in which
and
k is the number of iterations. To avoid the uncertainties and ill-conditioned problems caused by the coherence matrix inversion operation, we use the coherence weight
in the PL algorithm [
11,
33,
34].
2.2. Stepwise SB Phase Optimization
The increasing number of SAR satellites provide large amounts of SAR images at unprecedented temporal resolutions, which require a large amount of interference phase and make temporal phase optimization very time-consuming. Although strategies involving SB interferograms [
17] and block interferograms [
15] have been proposed, they cannot dynamically implement the temporal phase filter to improve the phase stability step-by-step in cases when a new SAR image is obtained and new interferograms are generated. It makes temporal phase optimization very time-consuming.
Inspired by the truncated SB interferograms to dynamically deal with real-valued unwrapped triangular error using the SLS method [
30], the temporal phase filter can be considered as a network adjustment with complex value. This strategy is suitable for existing phase triangulation algorithms, such as PL, SqueeSAR, and EVD. Therefore, we use the truncated SB interferograms algorithm to dynamically implement the temporal phase under the PL. To maintain the integrity of interferometric fringes in interferograms with dense interferometric fringes, we use the small multi-look number and implement the Goldstein filter [
11,
35] in the frequency domain to improve the phase stability, instead of averaging the SHPs. Phase triangle algorithms are dependent on a redundant phase. Conventional SB interferograms are used to perform temporal phase optimization, as shown in
Figure 2A [
17]; thus, we use adequate redundant SB interferograms to perform PL, as shown in
Figure 2B.
In addition, we use a vector calculation strategy for the PL algorithm, which is very fast and does not need to loop for each pixel; however, this strategy consumes more memory. The proposed strategy includes spatial domain filtering (i.e., multi-look), frequency domain filtering (i.e., Goldstein filter), and time domain filtering (i.e., PL).
Redundant interferograms make the optimization processing of temporal SB interferograms meaningful. For the archived SAR data, to obtain more redundant SB interferograms, we set the
K nearest SAR to generate SB interferograms
K × (
K − 1)/2 to perform PL. When we obtain a new SAR image, new
K − 1 interferograms are generated with
K-1 nearest SAR. Then, the appropriate multi-look and Goldstein filters are used for each new SB interferogram. Next, we use the optimized SB interferograms (
K − 1 SAR related) and new
K − 1 interferograms to update the related optimized interferograms and optimize new SB interferograms by the PL method with the sliding window (truncated SB interferograms network) with one SAR data step, as shown in
Figure 3. Note that when we obtain a new SAR image and related interferograms, in the block interferograms, the first column and first row interferograms are deleted and the new interferograms are added in the end column and end row.
2.3. Review Sequential SBAS InSAR Deformation Parameter Dynamic Estimation
For the archived SAR data, by setting the spatiotemporal baseline thresholds,
N archived SAR images can generate
M1 interferograms. After obtaining the spatiotemporal optimized phase, we perform two-dimensional phase unwrap** in the space domain; for example, using minimum cost flow (MCF) [
36]. The topography-related atmospheric phase is removed from the best-fitting linear relation between the phase delay and topography [
37]. The digital elevation model error-caused phase can then be estimated using the relationship between the perpendicular baseline and phase [
38]. We model
M1 redundant unwrapped interferograms using the Gauss-Markov model as
, and estimate the deformation time series as
. The subscript
1,
,
,
,
, and
represent the related archived data, design matrix, observations data, weight matrix, estimated deformation time series, and its cofactor matrix, respectively. The superscript
T represents the transposition of a matrix.
When we obtain a new SAR image,
M2 new interferograms are generated by combining the new SAR image with recent acquisitions using the spatiotemporal baseline thresholds. Then, after two-dimensional phase unwrap** and error correction for the new
M2 interferograms, we model the new
M2 unwrap** interferograms and dynamically estimate the deformation time series using Equation (3), where the weight matrix is
and the design matrix is extended to
A2 and
B [
28]:
where
is the updated archived deformation time series,
is the new cumulative deformation,
is the new cofactor matrix, and
is the gain matrix.
2.4. Sequential ISBAS Deformation Time Series and Quality Assessment
For one generic pixel, conventional sequential SBAS InSAR techniques require each SB interferogram to contain qualified pixels to retrieve the deformation time series. If an interferogram does not meet the requirement, the pixel will be discarded, reducing the number of available monitoring points. To obtain more monitoring targets and avoid the selection of interferograms, the adaptive SB interferogram (i.e., ISBAS) is introduced into the sequential InSAR processing. It should be noted that, unlike in the conventional sequential SBAS InSAR method, in sequential ISBAS, the archived unwrap** interferogram and its design matrix , as well as the new unwrap** interferogram and its design matrix , are adaptively processed for each pixel. This means that the design matrix and observation data are variable for each pixel.
The condition number of the design matrix (SB interferograms) for each pixel is used to identify network connectivity and reliable monitoring point pixels. A very large number of conditions means that the SB interferograms is so ill-conditioned that the pixel will be discarded, while more than one SB subset will lead to the lack of deformation datum, resulting in bias of deformation time series.
In the concept of adaptive interferograms, i.e., the ISBAS method proposed by [
22,
23], the accuracy of the ISBAS-derived deformation results is lower than that of the SBAS-derived results by statistical analysis, as was performed by [
24]. In InSAR data processing, phase decorrelation, mis-selected pixels, and unwrap** error will inevitably lead to the uncertainty of the deformation parameters.
Next, adaptive SB interferograms are employed for each pixel to obtain deformation parameters. Redundant interferograms can reduce the influence of the error phase and improve the stability of deformation parameters, which is a key step to determine the number of final monitoring points and the accuracy of deformation parameters. In geodetic network adjustment, residual noise/error can be expressed as Equation (4) [
39,
40]:
where the
R matrix is determined by the SB interferograms (design matrix
A),
is the noise, and
I is the identity matrix.
represents the number of redundant observations for SB interferograms. The residual error
is determined by the matrix
and noise
.
is the redundant observation component (
). More redundant observations result in a higher
value. When
= 0, the number of redundant observations is 0, which is a chain SB interferograms. All noise of observations is transferred to the deformation parameter, resulting in the fluctuation of the deformation time series. When
= 1, the observation error is totally transferred to the residual error, which means the observation is unnecessary.
We performed 1000 simulations of the relationship between the residual noise percentage of the deformation parameter and the number of redundant observations, and the average of the results is shown in
Figure 4. The blue line is the residual noise percentage in the estimated deformation parameter and the green line is the redundant observation component.
Figure 4 illustrates that the more redundant the SB interferograms are, the smaller the noise content in the deformation time series is.
The temporal coherence of least squares residual was proposed to measure the quality of the phase unwrap** [
41]. The accuracy of the estimated parameters is not only affected by residual error, but also by the SB interferograms. Inspired by geodetic network adjustment, we assessed the accuracy of the estimated parameters using the variance-covariance matrix in the Gauss-Markov model. It contains the number of redundant interferograms, the residual error, and the SB interferograms network structure. The variance-covariance matrix was used to describe the accuracy of geodetic network adjustment. Therefore, it is appropriate to assess the accuracy of sequential ISBAS InSAR-derived deformation parameters. For archived SAR data, the variance-covariance matrix of the deformation time series
is written as Equation (5):
where
is variance of unit weight. When we obtain a new SAR image in Equation (3), the variance-covariance matrix of the deformation time series
is written as Equation (6):
The accuracy of the estimated deformation parameters can be assessed through the variance-covariance matrix, which consists of two parts: and . In unit weight variance, the sum of the squares of the residual phase is an absolute precision index of the deformation parameters, which determines their reliability. The weight or cofactor matrices of the estimated deformation parameters can then be obtained by error propagation, which determines the relative accuracy among the deformation time series vector from the observation (unwrap** interferograms).
We use five indices to describe the quality of the deformation parameters in sequential ISBAS processing,
described as follows:
- (1)
The number of redundant interferograms is used to describe adaptive SB networks pixel-by-pixel;
- (2)
The sum of the residual is used to describe the phase unwrap** error, which is the absolute precision index of deformation parameters;
- (3)
The average trace value of the cofactor matrix is used to describe the relative accuracy of the deformation parameters;
- (4)
The average of the standard deviation (STD) is used to describe the quality of the deformation rate;
- (5)
The STD is used to describe the accuracy of the deformation time series.
3. Study Area SAR Dataset
InSAR monitoring of cree** landslides is a challenging task in mountainous regions because of decorrelation [
42,
43]. To test the proposed approach, we take the Baige landslide as an example, which occurred on the left bank of the **sha River in China over two separate days, on 10 October and 3 November 2018. The failed material (approximately 25.13 × 10
6 m
3) rushed into the **sha River and a total of 67,449 inhabitants were affected [
44,
45]. The post-failure of the Baige landslide was mapped from 8 November 2018 to 9 December 2019 using multiple radars; however, the resulting spatial-temporal deformation pattern is incomplete due to the large deformation gradient and decorrelation [
46]. We employ 72 Sentinel-1 SAR data images acquired from 8 November 2018 to 2 April 2021 to demonstrate the applicability of the proposed approach and investigate the deformation pattern of the Baige landslide.
Multi-look SB interferograms have bias [
47,
48,
49,
50]. Redundant interferograms can facilitate the identification and correction of the phase inconsistency and bias. In this study, we use block-shaped SB interferograms to mitigate the bias by the PL algorithm, in which full interferograms are generated for each block SAR data. Redundant interferograms characterized by a not so small temporal baseline can mitigate the effect of the bias, to some extent. The mitigation of the bias is proportional to the number of redundant interferograms. In fact, the dense fringes or sparse spatial pixels make phase unwrap** challenging. How to select the size of block-shaped SB interferograms for stepwise phase optimization? It is a compromise strategy, in which we should consider both the availability of the optimized phase and computational efficiency.
In the experiment, because of the large deformation gradient and decorrelation in the post-failure of the Baige landslide area, we set the
K = 4 nearest SAR (i.e., 40 days in Sentinel-1) to generate 210 redundant SB interferograms, as seen in
Figure 5. In which, each SAR data connect the nearest six SAR data and generate six interferograms.
5. Discussions
Figure 11 shows the enlarged Baige landslide deformation from the black rectangular marks the area of
Figure 9D, in which
Figure 11A,B illustrates the line-of-sight direction deformation rate and quality evaluation index (i.e., STD) of the proposed approach. In
Figure 11, the blue polygon indicates the collapsed area in November 2018, and shows that the south side of the landslide is being deformed.
The comparison of InSAR- (red dot) and GPS LOS-projected (blue dot) deformation time series is shown in
Figure 12, which is located at pixel P3 in
Figure 11A. To some extent, the results demonstrate the reliability of InSAR to obtain deformation in the Baige landslide.
Figure 11.
Post-failure Baige landslide deformation. (
A,
B) are the linear deformation rate and average STD of deformation time series, respectively. The blue polygon indicates the collapsed area in November 2018 and shows that the south side of the landslide is being deformed. The pixels P1 and P2 indicated by the arrows are shown in
Figure 13.
Figure 11.
Post-failure Baige landslide deformation. (
A,
B) are the linear deformation rate and average STD of deformation time series, respectively. The blue polygon indicates the collapsed area in November 2018 and shows that the south side of the landslide is being deformed. The pixels P1 and P2 indicated by the arrows are shown in
Figure 13.
The deformation time series of pixels P1 and P2 are shown in
Figure 13, in which
Figure 13A,B is the deformation time series with different STD (51 and 3 mm, respectively), and
Figure 13C,D shows the redundant observation components with different number interferograms (61 and 136, respectively).
The STD of deformation time series is similar to the method described by [
52]. The deformation time series of point P1 fluctuated, while that of P2 was smooth. P1 had sparse SB interferograms, as can be seen in
Figure 13C,D. If the redundant observation component is zero, any noise or error phase may lead to the deviation of the subsequent deformation time series [
39,
40].
The proposed method includes two key steps, which are to improve the big InSAR dynamic processing and obtain the deformation time series with dense coherent targets in areas of serious decorrelation. Firstly, the frequency domain filtering method [
11,
35] can maintain the phase integrity, avoiding obscuration of the phase continuity by averaging spatial SHPs [
51], which are shown in
Figure 7. In the stepwise SB phase optimization, we generate redundant candidate SB interferograms to improve the quality of the interferometric phase through the PL method, in which we should set a compromised threshold
K ≥ 3 to generate block interferograms. As the amount of SAR data increases, compared to SBAS and full-based interferograms methods [
9,
17], the number of interferograms in the proposed method remains consistent, which avoids having large amounts of interferograms and improves the computational efficiency in temporal phase optimization. In addition, we use vector calculation for the PL algorithm, in which all pixels can be processed without iteration for each pixel. Although this process consumes more memory, it still improves the efficiency of large-scale SAR data processing.
Secondly, in the big SAR era, the stepwise SB phase optimization method facilitates sequential ISBAS processing. Compared with the conventional sequential SBAS-InSAR method [
28], the proposed sequential DS-ISBAS method avoids the selection of interferograms. In this experiment, we set one subset of adaptive interferograms to update the deformation time series. In special cases, such as missing SAR data, we set multiple subsets of adaptive interferograms. Although the conventional ISBAS method can improve the density of monitoring points, its reliability needs to be evaluated. Meanwhile, the quality assessment method is applicable to the SBAS processing process. For the large gradient deformation, Massonet proposed the maximum deformation gradient that InSAR can monitor, which is the ratio of half the wavelength to the pixel size [
53].