An Improved SBAS-InSAR Processing Method Considering Phase Consistency: Application to Landslide Monitoring in Hualong County, Qinghai Province, China

TL;DR

This paper introduces a new InSAR method to improve phase consistency for better landslide monitoring in Hualong County, China.

Contribution

The novel NCP-SBAS workflow enforces phase consistency to correct fading biases and unwrapping errors in deformation estimation.

Findings

Fading biases affect deformation estimates based on magnitude and direction.

Unwrapping errors lead to underestimation of deformation.

NCP-SBAS shows significant advantages in landslide monitoring.

Abstract

Phase consistency is a critical prerequisite for achieving high-precision time-series InSAR deformation retrieval. However, conventional SBAS-InSAR methods provide only limited consideration of phase consistency during the inversion process. Within the SBAS-InSAR workflow, two principal categories of error sources are primarily responsible for phase inconsistency, manifested as non-zero closure phase (NCP): (1) fading biases introduced during multilooking and filtering prior to phase unwrapping; and (2) unwrapping errors caused by large deformation gradients, low coherence, or inappropriate selection of unwrapping algorithms. To address these issues, this study introduces an improved SBAS-InSAR processing workflow, termed NCP-SBAS, designed to improve the accuracy of deformation field estimation and thereby enhance its applicability to geological hazard monitoring. The key idea of the method is to enforce phase consistency as a constraint, jointly accounting for the spatiotemporal characteristics of fading biases and the valid deformation signals, thereby enabling effective correction of NCP. To evaluate the effectiveness of NCP-SBAS, this study conducted a detailed analysis of deformation differences in Hualong County, Qinghai Province, before and after NCP correction, highlighting the significant advantages of the proposed approach. The results indicate that the influence of fading biases on deformation estimates depends on both the magnitude and direction of deformation, while unwrapping errors primarily lead to an underestimation of deformation. In addition, the study provides an in-depth discussion of how fading biases and unwrapping errors affect landslide monitoring and identification.

Figures

Full text

1. Introduction

Multi-temporal Interferometric Synthetic Aperture Radar (MT-InSAR) has become one of the most powerful tools for retrieving high-precision deformation fields and has been widely employed in the monitoring of geological hazards, such as landslides, land subsidence, and mining deformation [1,2,3,4,5]. By exploiting stacks of interferograms, MT-InSAR enables the extraction of long-term, spatially dense deformation time series, providing valuable insights into hazard evolution. Commonly used MT-InSAR approaches include Persistent Scatterer InSAR (PS-InSAR), Small Baseline Subset InSAR (SBAS-InSAR), SqueeSAR, and StaMPS [6,7,8,9,10]. Among these, SBAS-InSAR, which relies on distributed scatterers (DS) and time series inversion of unwrapped interferograms, has become one of the most widely applied methods for geological hazard monitoring [11,12,13,14].

Effective error correction is a prerequisite for obtaining reliable deformation fields. In conventional SBAS-InSAR workflows, well-established methods exist for mitigating error components with distinct spatial low-pass characteristics, such as topographic residuals, orbital inaccuracies, and atmospheric delays. However, temporal error sources that compromise phase consistency have received comparatively little attention. Phase consistency is fundamental for robust phase estimation and accurate time series inversion in MT-InSAR [10,15,16]. According to the phase closure principle, the theoretical sum of phase differences around a closed loop of three or more interferograms should equal zero. Yet, in practice, two major error sources disrupt this principle: (1) fading bias introduced by multilooking and filtering operations; and (2) phase unwrapping errors caused by large deformation gradients, decorrelation, and algorithmic limitations. These errors give rise to nonzero closure phase (NCP) both before and after unwrapping, which accumulates along the temporal baseline and ultimately degrades deformation accuracy.

Recent studies have proposed a variety of strategies to mitigate fading bias. SBAS-InSAR commonly applies multilooking and filtering to enhance interferometric phase quality, but this inevitably sacrifices spatial resolution and mixes signals from scatterers with inconsistent temporal behavior, introducing fading bias [17]. To reduce such effects, closure triplets have been used to build observation matrices, with singular value decomposition (SVD) employed to estimate interferogram-specific biases [18]. However, when the number of closure loops is smaller than that of interferograms, the inversion becomes underdetermined, leading to large uncertainties. To overcome this, some approaches introduce linear dependence assumptions among closure loops with different temporal baselines [19], while others apply averaging strategies based on the law of large numbers [20]. More advanced techniques model amplitude evolution across long-span closure stacks [21] or exploit time-invariant and time-variant fading models with temporal windowing [22]. Despite these developments, existing methods either rely on oversimplified assumptions or fail to capture localized and spatially heterogeneous bias characteristics.

Similarly, unwrapping errors represent another critical challenge. Phase unwrapping aims to reconstruct continuous phase values from wrapped observations in the interval [−π, π], enabling deformation retrieval [23]. In theory, correctly unwrapped interferograms forming closed loops should yield zero closure phase. Thus, NCP serves as a sensitive indicator of unwrapping errors [24,25], particularly in two-dimensional algorithms, where NCP reflects temporal–spatial consistency. To correct such errors, two main categories of approaches have been developed. The first focuses on global correction via observation matrices and least-squares inversion, with regularization strategies such as L1-norm or Tikhonov frameworks applied to improve stability [26,27]. The second emphasizes localized corrections, integrating phase closure with morphological and re-estimation techniques [28]. However, these methods become ineffective when multiple interferograms within a closure loop simultaneously contain errors, leading to compounded ambiguities. To address this issue, iterative correction frameworks combining phase gradients, closure phases, and spatial coherence have been introduced [29]. While effective for large-scale regions, their applicability to small-scale hazards, such as localized landslides, remains limited due to mismatched detection window sizes.

To address the two aforementioned sources of error, this study proposes an improved SBAS-InSAR processing approach that explicitly incorporates the principle of phase consistency. By leveraging NCP derived from closure loops both before and after unwrapping, we aim to generate more accurate deformation fields and quantitatively assess the impacts of phase errors on geological hazard monitoring.

The remainder of this paper is organized as follows. Section 2 introduces the methodology for fading bias and unwrapping error correction. Section 3 presents the datasets and experimental results. Section 4 provides further discussion. Finally, Section 5 concludes the paper.

2. Methods and Data Processing

The proposed SBAS-InSAR workflow in GAMMA (version 2.9) includes data pre-processing, multilooking, spatiotemporal network construction, differential interferometry, phase unwrapping via the Minimum Cost Flow algorithm, orbital and atmospheric corrections, and time-series inversion using SVD. Two additional modules: (1) fading bias correction and (2) pixel-wise unwrapping error correction in deformation zones, are integrated into the standard pipeline, as detailed below.

2.1. Triplet Phase Closure

Assuming a dataset of N SAR acquisitions, a spatiotemporal baseline network can be constructed to generate M interferograms, with the interferometric phases before and after unwrapping denoted as $[eqn]$ and $[eqn]$ , respectively. Let $[eqn]$ denote the acquisition times of the master and slave images. Select three interferograms that form a closed loop in time: $[eqn]$ . Ideally, the triplet phase closure (TPC) should be zero both before and after phase unwrapping [10,15,16], as shown in Equations (1) and (2).

$$[eqn]$$ $$[eqn]$$

Fading bias and unwrapping errors both disrupt phase consistency, producing nonzero closure phases. In the wrapped phase, fading bias $[eqn]$ dominates the deviation, while in the unwrapped phase, both fading bias and unwrapping error, highlighting the importance of correcting fading bias. Hereafter, closure phases due to multilooking and filtering are denoted $[eqn]$ , and those from unwrapping errors as $[eqn]$ .

2.2. Fading Bias Correction

Fading bias refers to phase deviations introduced by temporally inconsistent scattering behavior of multiple scatterers within a pixel, typically after multilooking and filtering. Its spatiotemporal characteristics are closely related to changes in surface scattering mechanisms. Intensity, a direct indicator of backscattering strength, effectively reflects regional scattering variations. For example, urban areas with stable scatters usually exhibit lower fading bias than agricultural areas, where soil moisture fluctuations result in higher temporal decorrelation [30,31]. Several studies have validated the effectiveness of amplitude-based fading bias correction methods [16,21]. Coherence, a key parameter for evaluating interferometric phase quality, is also closely associated with fading bias, determining its spatial distribution [30]. These observations suggest that both amplitude variation and coherence are important indicators of fading bias behavior. Moreover, fading bias demonstrates temporal decay characteristics, with stronger effects typically observed at shorter temporal baselines [19]. In this study, we account for the spatiotemporal behavior of fading bias in the correction process. The detailed correction steps are presented below.

2.2.1. Calculation of ∆NCPw, ∆ins, and ∆coh

Fading bias is sensitive to both spatial averaging and temporal decorrelation effects introduced by multilooking and filtering. All possible triangle loops are iterated for each pixel, and the $[eqn]$ , $[eqn]$ , and $[eqn]$ are computed accordingly. For both multilook ratios of 10:2 and 5:1, the wrapped phase, amplitude values, and coherence values are extracted. The difference between multilook ratios 5:1 and 10:2 captures the sensitivity of $[eqn]$ to fading bias while preserving fine-scale deformation. Using Equation (1), the $[eqn]$ is calculated. To avoid the suppression of subtle deformation signals that may arise from excessively large multilooking windows, this study employs multilook ratios of 5:1 and 10:2. These configurations achieve a balanced trade-off between noise reduction and the preservation of fine-scale deformation information, ensuring that the subsequent closure-phase analysis retains both spatial detail and phase stability.

To address the attenuation of fading bias, temporal baselines of the interferograms were incorporated into the calculation of intensity and coherence variations. The changes in fading bias, intensity, and coherence were then computed using Equations (3)–(5). Consider a set of three interferograms composed of master-slave pairs $[eqn]$ , $[eqn]$ , and $[eqn]$ , forming a triangle loop. The key parameters are defined as follows: the closure phase $[eqn]$ and the change in intensity $[eqn]$ and coherence $[eqn]$ of each interferogram $[eqn]$ , where m denotes the multilook ratio. The temporal baseline for interferogram $[eqn]$ is denoted by $[eqn]$ . The relative contributions of intensity ( $[eqn]$ ) and coherence variations ( $[eqn]$ ) to fading bias were quantified. Due to heterogeneous deformation patterns and differential ground target composition across multilook configurations, the sensitivity of $[eqn]$ to these variations exhibited notable regional differences.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Following the extraction of $[eqn]$ , $[eqn]$ , and $[eqn]$ , it is necessary to establish a quantitative relationship between these indicators and the fading bias itself. Since both intensity and coherence exhibit spatially and temporally varying sensitivities to fading bias, their relative contributions cannot be assumed to be constant across the scene or between different multilooking configurations. Therefore, instead of assigning empirical or fixed weights, we develop a data-driven weighting model that characterizes the joint influence of intensity and coherence variations on fading bias.

2.2.2. Construction of a Data-Driven Pixel Weighting Function

To quantify the contribution of each pixel to the observed closure-phase variations, we model the empirical relationship between fading-bias-related phase deviations and two key indicators $[eqn]$ and $[eqn]$ using a multivariate linear formulation. This strategy allows the weighting process to adapt to spatially heterogeneous scattering environments and naturally incorporate the temporal characteristics of fading bias.

(1) Data normalization. Because the dynamic ranges of $[eqn]$ , $[eqn]$ , and $[eqn]$ differ significantly, all variables are normalized to the interval $[eqn]$ using Min-Max scaling to ensure numerical stability and to avoid parameter bias.

(2) To capture the potentially coupled effects of $[eqn]$ and $[eqn]$ on fading bias-induced phase variations, we employ a multivariate linear regression model with an intercept and interaction term:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the primary contributions of amplitude and coherence variations; $[eqn]$ captures their interaction, reflecting the fact that fading bias is often jointly influenced by scattering strength and temporal decorrelation; $[eqn]$ is the intercept.

(3) Mapping regression output to a physically meaningful weight. Since the regression output is not restricted to $[eqn]$ , its predicted value is transformed into a physically interpretable pixel weight using a sigmoid function:

$$[eqn]$$

This transformation ensures that the resulting pixel-wise weight satisfies $[eqn]$ . $[eqn]$ represents the predicted normalized closure-phase deviation obtained from the multivariate linear regression model. It reflects the expected contribution of each pixel to the closure-phase noise given its intensity and coherence values, and serves as the input to the sigmoid function to generate the physically interpretable pixel weight $[eqn]$ . For each interferogram, the same normalization and regression model can be applied to obtain pixel-wise weights, which are then used in the weighted closure-phase SVD estimation.

2.2.3. Bias Estimation Based on Weighted SVD

Given a dataset of $[eqn]$ SAR images with spatial dimensions $[eqn]$ , $[eqn]$ interferograms are generated. From these, $[eqn]$ triangle loops are formed, with each $[eqn]$ computed using Equation (1). Based on the learned pixel-wise weights, a diagonal weighting matrix $[eqn]$ is constructed, and the corresponding weight matrix $[eqn]$ is defined as:

$$[eqn]$$

A design matrix $[eqn]$ is formulated based on triangle phase closure (TPC) relationships, and the unknown vector $[eqn]$ represents the fading bias for each interferogram at the pixel level. By incorporation the pixel-wise weights, the weighted design matrix and weighted observation vector are defined as:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ denotes the $[eqn]$ observations at pixel location $[eqn]$ .

The weighted design matrix is then decomposed via SVD:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are orthogonal matrices, and $[eqn]$ is a diagonal matrix containing the singular values of $[eqn]$ .

The inversion is performed independently for each pixel $[eqn]$ using the Moore–Penrose pseudoinverse of $[eqn]$ . The pseudoinverse $[eqn]$ id obtained by inverting the non-zero singular values in $[eqn]$ . To improve numerical stability and suppress the influence of ill-conditioned components, singular values smaller than a predefined threshold are discarded during the inversion. The fading bias vector is then estimated as:

$$[eqn]$$

The SVD-based inversion is performed independently for each pixel $[eqn]$ . The resulting $[eqn]$ values are then assembled to form a two-dimensional fading bias field for each interferogram.

Finally, the estimated fading bias is removed from each interferogram to obtain the corrected interferogram stack:

$$[eqn]$$

where $[eqn]$ denotes the corrected phase of interferogram $[eqn]$ .

2.3. Unwrapping Error Correction in Deformation Areas

After excluding low coherence points, rapid spatial variations caused by severe atmospheric delays and surface deformation are the primary factors leading to phase unwrapping errors [24,28]. $[eqn]$ based correction methods require a sufficient number of correctly unwrapped interferograms as reliable references. For triangle loops without error compounding, single interferogram errors often associated with the longest temporal baselines allow the integer phase offset to be directly inferred from $[eqn]$ . Correction proceeds by prioritizing such non-compounded loops, rapidly expanding the set of correctly unwrapped pixels and enabling subsequent correction of more complex cases. Moreover, these corrected interferograms provide reliable deformation-phase constraints, serving as a priori knowledge for precise correction in deforming regions. In this study, we integrate a sequential iterative strategy with deformation phase constraints to enable pixel-wise unwrapping error correction while addressing the impact of phase compounding. The detailed implementation is presented in the following sections.

2.3.1. Initial High-Quality Phase Acquisition

The $[eqn]$ of each triangle loop serves as the criterion for the selection of an initial high-quality phase set. A threshold $[eqn]$ of 1 rad is applied to the $[eqn]$ values calculated from $[eqn]$ triangle loops. The threshold was determined empirically based on statistical analysis of correctly unwrapped pixels. We found that, under normal phase gradients, the closure phase consistently remains below 1 rad, whereas pixels affected by unwrapping errors exhibit significantly larger deviations. This empirically derived threshold therefore provides a reliable criterion for separating valid closure phases from those contaminated by unwrapping errors. Pixels with $[eqn]$ exceeding λ are masked to obtain the initial high-quality phase set, denoted as $[eqn]$ . Although this initial set may include pixels affected by error compounding, the high quality phase set is iteratively refined, allowing such pixels to be progressively identified and corrected during subsequent iterations.

2.3.2. Deformation Zone Localization and Reliable Deformation Information Extraction

Examination of pixels known to exhibit stable deformation behavior shows that their sign consistency defined as the proportion of interferograms sharing the same deformation sign typically exceeds 0.75, whereas pixels dominated by noise, atmospheric artifacts, or unwrapping uncertainties exhibit substantially lower consistency. Thus, $[eqn]$ provides a robust and conservative criterion that effectively separates deformation dominated pixels from non-deformation pixels. This threshold ensures high reliability of the deformation mask while retaining sufficient spatial coverage for subsequent correction procedures.

$$[eqn]$$

$[eqn]$ denotes the final set of deformation pixels, and $[eqn]$ represents the initial high-quality phase list for a given pixel. Using the deformation area mask, the deformation phase series of each deformation pixel across all interferograms is extracted. Incorporating the temporal baselines $[eqn]$ , phase sets following the deformation trend ( $[eqn]$ ) are selected. $[eqn]$ is the union of $[eqn]$ . The mean of $[eqn]$ , calculated by Equation (15) as $[eqn]$ , is used to replace the median value.

$$[eqn]$$

Using the median absolute deviation (MAD) method, the absolute deviation of each phase from $[eqn]$ calculated to obtain the set of absolute deviations $[eqn]$ (Equation (16)). Based on a threshold criterion ( $[eqn]$ ), outliers are identified and the phase values affected by compounding errors ( $[eqn]$ ) are further removed (Equation (17)). The resulting valid deformation phase set $[eqn]$ is then used to reestimate integer phase ambiguities (Equation (18)).

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

2.3.3. Re-Estimation of Integer Phase Ambiguities in Deformation Erroneous Zones

For $[eqn]$ in triangle loops with unwrapping errors, two scenarios can be distinguished. In the first scenario, errors occur in only one interferogram, allowing direct correction using the corresponding triangle loop. In the second scenario, multiple interferograms contain simultaneous unwrapping errors, generating complex interactions that obscure true phase discontinuities. To address this, a sequential iterative correction strategy is applied. Closure triplets are ranked by coherence, inferred from temporal baselines and $[eqn]$ values, and corrected in descending order. Each iteration updates the deformation estimates, gradually incorporating corrected pixels into the reliable phase set, which then contribute to subsequent corrections. For sufficiently redundant networks, this iterative procedure progressively converts compounded-error loops into effectively non-compounded cases, enabling robust unwrapping error correction in complex regions. The workflow comprises three main components:

(1) Calculate the sum of temporal baselines for each triangle loop, and use the closure phase residuals of triangle loop to derive the unwrapping nonzero closure phase set $[eqn]$ . These two parameters serve as sequential sorting references: $[eqn]$ as the primary (Equation (19)), $[eqn]$ as the secondary (Equation (20)). The triangle loops are then sorted in ascending order to obtain a ranked sequence $[eqn]$ (Equation (21)), which defines the correction order.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a weighting factor used to balance the relative contributions of the temporal baseline and the unwrapping-related closure phase residual. Both terms are normalized prior to combination to remove scale effects. Triangle loops are then sorted in ascending order of $[eqn]$ . Loops with smaller scores are considered more reliable and are corrected first, whereas loops with larger scores, typically associated with longer temporal baselines and stronger unwrapping inconsistencies are deferred to later iterations. This ordering ensures that reliable deformation information is progressively accumulated and used to resolve more complex unwrapping error patterns.

(2) Identify a set of valid deformation phases with the same temporal baseline to estimate the integer phase ambiguity. If the interferogram to be corrected has a temporal span of $[eqn]$ and the erroneous phase is denoted as $[eqn]$ , the valid deformation phases with the same temporal span, $[eqn]$ , can be extracted from $[eqn]$ (Equation (22)). The mean value $[eqn]$ of $[eqn]$ (Equation (23)) is taken as the reference value for estimating the integer ambiguity $[eqn]$ (Equation (24)). The resulting $[eqn]$ is then applied to correct the errors in the original phase.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

(3) Recalculate the $[eqn]$ of the triangle loops that include the corrected interferograms to assess the effectiveness of the correction. The successfully corrected pixels are continuously incorporated into the valid deformation phase set, enabling iterative refinement. This process is repeated until all abrupt phase errors are corrected, resulting in the final corrected interferogram set $[eqn]$ (Equation (25)).

The NCP-SBAS data processing workflow was established, as shown in Figure 1.

3. Study Area and Datasets

The study area is situated in the southern part of Hualong County, Qinghai Province, in northwestern China, within the transitional zone between the Qinghai–Tibet Plateau and the Loess Plateau (Figure 2). The region has an average elevation exceeding 2800 m and exhibits a typical plateau continental climate, with annual precipitation ranging from 300 to 500 mm, primarily occurring between June and September. The terrain is highly dissected, characterized by steep slopes, deeply incised valleys, and intense fluvial erosion. Quaternary loess deposits are present in portions of the area, which are prone to instability and exhibit high susceptibility to water-induced softening and collapse. In addition, extensive bedrock weathering has formed fractured zones. Due to these geomorphological and geological conditions, the region is highly susceptible to landslides, rendering it an ideal site for landslide monitoring and related geohazard studies. The fundamental data processing steps for the study area—including data preprocessing, multilook filtering, spatiotemporal network formation, differential interferometry, phase unwrapping, orbital error correction, atmospheric delay compensation, and SVD-based time-series inversion—were all performed using the GAMMA software (version 2.9).

A total of 115 Sentinel-1 IW-mode single-look complex (SLC) images from 7 March 2015 to 23 December 2019 were acquired and organized into networks of consecutive four image sets (Figure 3), yielding 450 interferograms. A coherence threshold of 0.4 was applied, and phase unwrapping was conducted using the Minimum Cost Flow (MCF) algorithm. A multilook ratio of 10:2 was used for deformation inversion, while bias sensitivity estimation employed multilook ratios of 5:1 and 10:2. Goldstein filtering was applied for noise reduction. Monthly precipitation data for Haidong City were obtained from the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/, assessed on 15 July 2024).

4. Results and Analysis

4.1. Bias Sensitivity Analysis

Using three wrapped interferograms (denoted as 20150307_20150331, 20150331_20150424, and 20150307_20150424) as an illustrative example, we analyze the relative contributions of amplitude and coherence variations to the three-interferogram $[eqn]$ at the pixel level. A pixel-wise sensitivity index $[eqn]$ is defined as:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the regression coefficients corresponding to $[eqn]$ and $[eqn]$ , respectively. By construction, $[eqn]$ indicates that the closure phase is predominantly influenced by coherence variations, while $[eqn]$ indicates that amplitude variations are the dominant factor. Finally, the bias sensitivity map was computed using Equation (26), as illustrated in Figure 4d. According to Figure 4, $[eqn]$ (Figure 4a) and $[eqn]$ (Figure 4c) exhibit a generally positive correlation at the regional scale. Larger $[eqn]$ values are typically associated with more pronounced variations in $[eqn]$ . Since $[eqn]$ is negatively correlated with interferometric coherence, fading bias tends to increase as coherence decreases. This finding is consistent with the results reported by [30]. Distinct sensitivity differences among various surface features can be observed in Figure 4d, allowing for clear delineation of land-cover boundaries.

Three representative regions were selected for further analysis, as indicated by the black solid boxes (A, B and C) in Figure 4d. Region A corresponds to a largely unvegetated mountainside. Despite noticeable amplitude variations, it exhibits relatively low fading bias and minor changes in coherence. Appearing in blue in Figure 4d, this region is minimally affected by fading bias, indicative of high pixel quality and a strong correlation between bias and coherence. Region B is an agricultural area, where the spatial distributions of fading bias and amplitude variation are closely aligned, as clearly seen in Figure 4a,b. Coherence changes are moderate. The region appears primarily in red in Figure 4d, suggesting a strong dependence of fading bias on amplitude variation. Region C is located near a watershed with intensive human activity. It shows large amplitude fluctuations and significant coherence degradation, along with the highest bias levels. This highly unstable region appears light purple in Figure 4d, near the median sensitivity level, indicating that fading bias here is influenced by both amplitude variation and coherence.

4.2. InSAR Deformation Results and Analysis

We obtained the mean deformation rates after (Figure 5a) NCP correction. In addition, we quantified the errors arising from fading bias (Figure 5b) and from unwrapping errors (Figure 5c). The combined effect of these error components equals the total correction impact.

The spatial distribution of deformation signals remains largely consistent before and after correction. However, significant deformation zones demonstrate systematic underestimation in uncorrected results. Figure 5c further reveals that unwrapping errors typically lead to underestimated displacement rates. The maximum difference before and after correcting unwrapping errors reaches up to 8.5 cm/yr, with an average difference of 34 mm/yr. As shown in Figure 6b, the fading bias exhibits significant spatial variability. The maximum difference before and after bias correction reaches 7.12 mm/yr, with an average difference of 1.02 mm/yr.

Two main regions are notably affected by fading bias: the area within the solid rectangular box in Figure 6b, and the deformation zone. The area in the box shows pronounced systematic bias, manifesting as deformation signal overestimation. This results from steep mountainous terrain and complex scattering mechanisms that amplify spatiotemporal backscattering inconsistencies. In the deformation zone, fading biases are concentrated in regions with large deformation magnitudes and those along deformation boundaries. The bias magnitude is proportional to the deformation level. Four representative landslide areas, labeled 1, 2, 3 and 4 were selected for detailed analysis (Figure 5a). For each area, we extracted the average deformation velocity along with the corresponding corrections for fading bias and unwrapping errors.

As no unwrapping errors were detected in Area 4, its unwrapping errors are not shown in Figure 6. Both fading bias and unwrapping errors are predominantly distributed in regions with larger deformation magnitudes (Figure 6a–j). Additionally, fading bias are often concentrated along the deformation zone edges. For areas with minor overall bias trends (outside the solid-line rectangle), the bias inside deformation zones aligned with the displacement direction, resulting in fading bias underestimation. As shown in Figure 6k, in regions with significant overall bias trends (outlined by the solid-line rectangle in Figure 5b), the bias magnitude correlates with deformation magnitude. In one deformation zone, large-deformation regions exhibit biases align with the deformation direction, whereas small-deformation regions have biases aligned with the overall bias trend.

Based on the cumulative displacement statistics within the deformation zones, the average difference in cumulative displacement is 5.1 cm, with a maximum difference reaching 15 cm. To further investigate the time series, we selected three representative ponit1, point2 and point3 in Figure 6a,d,g. For each point, we analyzed three displacement time series: (1) uncorrected raw data, (2) unwrapping error corrected only, and (3) fully corrected (both unwrapping errors and fading bias). As shown in Figure 7, unwrapping errors result in deformation underestimation. The influence of fading bias varies across regions. At point3 (negligible fading bias), the corrected and uncorrected time series exhibit minimal fluctuations, and the final cumulative displacements are nearly identical. This indicates that unwrapping errors are the primary error source at this point. Both point1 and point2 are affected by fading bias, with before and after cumulative displacement differences of 14 mm and 8 mm, respectively. After correction, the fading bias leads to underestimation. In addition, we obtained the monthly average precipitation data for the study area in March 2015 from to December 2019. A sharp increase in rainfall occurred in July and August 2018, which triggered accelerated deformation rates (Figure 7). Notably, both unwrapping errors and fading biases became more pronounced following the rainfall event.

5. Discussion

5.1. In Effectiveness of Sequential Iteration for Unwrapping Errors Correction in Error Compounding Regions

When correct, unwrapping errors in deformation zones, $[eqn]$ serve two key functions: identifying areas of discontinuity prior to correction, and validating the correction effectiveness. As shown in Figure 8, in regions with either corrected or inherently absent unwrapping errors, the three interferograms forming a triangle loop maintain consistent deformation, resulting $[eqn]$ approaching zero.

In triangle loops affected by phase unwrapping errors, longer temporal baselines are generally associated with lower coherence and larger deformation magnitudes. In a triangle loop with relatively short average temporal baselines, unwrapping errors are most likely to occur in the interferogram with the longest baseline [29]. However, longer temporal baselines and reduced coherence lead to increasingly complex unwrapping errors patterns. When two of the three images contain same-sign errors, unwrapping error compounding may occur, whereby the phase errors cancel out during the triangle loop computation. In such cases, the $[eqn]$ fails to capture the actual unwrapping error, leading to incorrect localization of phase errors and significantly compromising the correction accuracy. Therefore, sequential and iterative correction strategies are essential for effective unwrapping errors mitigation. To investigate the necessity of sequential iteration for correcting unwrapping discontinuities in deformation zones, we analyzed a representative area using two triangle loops: ① 20161120_20170107, 20170107_20170131, and 20161120_20170131; and ② 20161120_20170131, 20170131_20170212, and 20161120_20170212. For clarity, letters are used to denote interferogram IDs. Group ① includes $[eqn]$ , $[eqn]$ , and $[eqn]$ ; Group ② includes $[eqn]$ , $[eqn]$ , and $[eqn]$ .

We extracted unwrapped interferograms from Groups ① and ②, with orbital and atmospheric errors corrected. Figure 9a show group ① and Figure 9d show group ②. In Group ①, no compounding occurs. The phase discontinuity appears in $[eqn]$ , which has the longest temporal baseline. The computed $[eqn]$ (Figure 9b) enables accurate localization of the error region. In contrast, Group ② exhibits aliasing, as both $[eqn]$ and $[eqn]$ experience unwrapping error with the same sign. These errors cancel out during $[eqn]$ calculation, as shown in Figure 9e, preventing effective identification of the unwrapping errors. Specifically, the phase unwrapping errors in $[eqn]$ can be corrected using the method described in Section 2.3, yielding a corrected interferogram $[eqn]$ (Figure 9c). This corrected interferogram is then included in the high-quality dataset and contributes to the correction of Group ②, as illustrated in Figure 9f. In this case, compounding is eliminated, and the $[eqn]$ result (Figure 9g) successfully identifies the unwrapping error region in $[eqn]$ . This iterative correction continues in a similar manner, effectively mitigating the impact of compounding on unwrapping error correction. For the landslide in Area 3 of Figure 6a, the deformation zone was treated as an integrated unit for correction. Both non-sequential and sequential iterative strategies were applied. Among the 450 interferograms, 181 phase discontinuities were detected. The non-sequential correction strategy addressed only 88 errors, while the sequential iterative strategy successfully corrected 157.

5.2. Comparison with Other Fading Bias Correction Method and Validation of Effectiveness

To further validate the effectiveness of the proposed method, we compared the fading bias correction performance of $[eqn]$ with the $[eqn]$ (Time Variant model) [22]. The $[eqn]$ adopts a 72 day temporal window for bias estimation and correction. Under error-free conditions, deformation time series derived from the same set of SAR images, regardless of different spatiotemporal baseline networks, should be consistent. Therefore, any discrepancy in the results obtained from different baseline configurations can serve as an indicator of correction performance.

In this study, we conducted deformation time series analysis using two baseline configurations: the full baseline set, constructed with connections to the subsequent four images, and a long baseline subset, obtained by removing interferometric pairs with temporal baselines shorter than 36 days from the full set. Differences in the estimated average deformation rates between the two configurations are shown in Figure 10. To isolate the specific contribution of fading bias correction, this comparative analysis intentionally excluded phase unwrapping error correction. Overall, the mean deformation rate differences are 1.09 mm/yr, 0.91 mm/yr, and 0.62 mm/yr, respectively, demonstrating the superior performance of the $[eqn]$ approach.

To evaluate temporal consistency, we selected two representative points ( $[eqn]$ and $[eqn]$ ) where both $[eqn]$ and $[eqn]$ methods demonstrated effective correction. Figure 11a,b show the fading bias time series, derived by differencing with the original deformation data. While both methods exhibit strong temporal agreement, the $[eqn]$ method shows reduced effectiveness in active deformation zones, and areas with pronounced local bias trends. To further assess the correction performance in such areas, we conducted a statistical analysis on the region outlined by the black rectangle in Figure 10, where significant bias trends were observed. The histograms of deformation rate differences for the raw, $[eqn]$ -corrected, and $[eqn]$ -corrected results are shown in Figure 11c, Figure 11d, and Figure 11e, respectively. The corresponding mean differences are 0.59 mm/yr, 0.51 mm/yr, and 0.30 mm/yr, confirmed that $[eqn]$ achieves the best correction performance. This improvement is likely due to the underlying mechanisms of the two methods. The $[eqn]$ model estimates systematic bias primarily based on temporal baselines, with limited consideration of spatial characteristics. In contrast, $[eqn]$ incorporates both the spatial sensitivity of fading bias and its temporal fading behavior, enabling a more comprehensive spatiotemporal correction.

In addition, the effectiveness of the proposed method is further validated and evaluated by comparing deformation estimates obtained from different baseline combinations. A deformation mask threshold is determined based on the standard deviation of the deformation results, and deformation areas are subsequently masked out using this threshold. The effectiveness of the correction is then assessed by analyzing the statistical characteristics of the non-deforming areas. Specifically, the statistical analyses of the non-deformation regions in Figure 10a,c are presented, with the corresponding histogram distributions and RMSE values shown in Figure 12. The results indicate that the RMSE of the deformation estimates obtained using the proposed method is 0.7 mm, compared to 1.5 mm for the original deformation results, corresponding to an overall accuracy improvement of approximately 53%.

5.3. Impacts of Fading Bias on Landslide Detection

Fading bias introduces systematic errors with pronounced spatial heterogeneity, which may lead to misestimation of the spatial distribution and magnitude of landslide deformation. For slowly moving or low-magnitude landslides, the deformation signal may be obscured by the bias. Therefore, this section conducts landslide identification within the study area and assesses the impacts of fading bias on detection. Landslide identification was conducted through visual interpretation of the deformation fields before and after $[eqn]$ correction, assisted by Google Earth.

A total of 108 landslides were detected before correction, primarily distributed along mountain ridges and riverbanks. After correction, 118 landslides were identified, including 10 newly recognized landslides (Figure 13a). These newly detected sites exhibit relatively small deformation extent and magnitude, with an average area of 13,775 m^2^ and a mean deformation rate of 1.61 mm/yr. To illustrate the improvement in landslide detection, a representative area outlined by the rectangle in Figure 13a is selected for detailed analysis (Figure 13b). Within this area, five new landslides were detected. Among them, three landslides (labeled L-1, L-2 and L-3 in Figure 13b) were selected for further analysis. The deformation fields before and after correction, along with corresponding optical imagery, are shown in Figure 13c–k, respectively. All the three landslides exhibit clear linear cracks along the rear or lateral margins of the hillslopes within the delineated boundaries in the optical images, consistent with typical visual features of landslides. Before fading bias correction, the deformation magnitudes are often comparable to the background noise level. These scattered and spatially limited signals are difficult to confidently identified as landslides. After bias correction, deformation signals exhibit significantly enhanced spatial coherence and recognizability.

6. Conclusions

This study proposes an improved SBAS-InSAR processing method that accounts for NCP. Correction strategies for fading biases and unwrapping errors in deformation areas are introduced to enhance phase consistency throughout the inversion process. This approach enables retrieve more accurate deformation fields and allows quantitatively assess the effects of fading biases on landslide detection and analysis. The key contributions of this work are as follows: (1) A weighted SVD-based method is proposed for fading bias correction. The approach first constructs a pixel-wise weighting model using amplitude and coherence information. Finally, $[eqn]$ is calculated based on triangle loops to achieve bias correction considering both spatial and temporal variability. (2) An iterative optimization based algorithm is developed to correct unwrapping errors in deformation areas. The method first identifies potential unwrapping errors and selects an initial set of high-quality interferograms by $[eqn]$ , on which basis deformation zones are accurately localized. Leveraging reliable deformation signals and a sequential iterative strategy, a pixel-wise correction model is established to address unwrapping errors, including cases affected by compounding. After error correction, the number of identified landslides increased from 108 to 118, corresponding to an improvement of approximately 9.3% in landslide detection capability. Compared with the uncorrected result and TVm, the proposed method reduces the mean difference from 0.59 mm/yr and 0.51 mm/yr to 0.30 mm/yr, corresponding to error reductions of approximately 49.2% and 41.2%, respectively. In terms of the differences in deformation results under different baseline combinations, the proposed method improves the accuracy by 53%. One limitation of this study lies in the $[eqn]$ estimation model, which does not yet fully account for seasonal variations, potentially introducing uncertainties in long-term monitoring. Moreover, the correction of unwrapping errors in deformation areas relies on the temporal consistency of deformation trends. When the deformation exhibits alternating uplift and subsidence, a more comprehensive deformation model is required to support accurate correction.