Mountain Landslide Monitoring Using a DS-InSAR Method Incorporating a Spatio-Temporal Atmospheric Phase Screen Correction Model

Abstract

The detection of potential rural mountain landslide displacements using time-series interferometric Synthetic Aperture Radar has been challenged by both atmospheric phase screens and decoherence noise. In this study, we propose the use of a combined distributed scatterer (DS) and the Prophet_ZTD-NEF model to rapidly map the landslide surface displacements in Diqing Tibetan Autonomous Prefecture, China. We conducted tests on 28 full-resolution SENTINEL-1A images to validate the effectiveness of our methods. The conclusions are as follows: (1) Under the same sample conditions, confidence interval estimation demonstrated higher performance in identifying SHPs compared to generalized likelihood ratio test. The density of DS points was approximately eight times and five times higher than persistent scatterer interferometry and small baseline subset methods, respectively. (2) The proposed Prophet_ZTD-NEF model considers the spatial and temporal variability properties of tropospheric delays, and the root mean square error of measured values was approximately 1.19 cm instead of 1.58 cm (PZTD-NEF). (3) The proposed Prophet_ZTD-NEF method reduced the mean standard deviation of the corrected interferograms from 1.88 to 1.62 cm and improved the accuracy of the deformation velocity solution by approximately 8.27% compared to Global Position System (GPS) measurements. Finally, we summarized the driving factors contributing to landslide instability.

1. Introduction

Landslides generally occur in geographically complex, weather-prone and ecologically fragile areas in most regions worldwide [1,2,3,4,5]. The developed differential interferometric synthetic aperture radar (D-InSAR) directly uses radar phase differences measured from satellite or airborne SAR images to measure displacement. This not only overcomes the limitations of Global Position System (GPS), precise leveling measurements, and high management costs associated with installing instruments in inaccessible or hazardous areas but also avoids the destruction of fragile surface environments, making monitoring more convenient and efficient [6,7,8,9]. However, the practical application and research of D-InSAR are still limited by inherent surface scattering noise [10,11], inaccurate orbital errors [12], digital elevation model (DEM) errors [13], and APS factors [14,15,16,17,18], which limit the applicability of correctly represented deformation signals. Over the past 25 years, a series of time-series analysis methods based on point targets rather than image pixels have been proposed. Among them, is the persistent scatterer interferometry (PSI) method, which uses phase variations over time of PS [19].The small baseline subset (SBAS) method based on spatial correlation, are two pioneering time-series InSAR (TS-InSAR) methods [20].

Due to the strong scattering of roads, buildings, and bridges, PSI offers high precision in urban areas [21,22,23]. Unfortunately, when applied to areas with vegetation cover or low reflectivity signals, the spatial distribution of measurement points (MPs) becomes uneven. Unlike the PSI method, the SBAS method typically focuses on processing point targets where backward scatterers are not dominant in the radar resolution pixel [24]. These slowly decorrelating filtered phase (SDFP) pixels are typically correlated within short time intervals, thus allowing sufficiently scatterers to be selected even in bare soil, vegetation, and arable land. This has made the SBAS method one of the most used techniques for analyzing multi-temporal interferometric stacks at low resolution. However, SBAS methods typically result in loss of image resolution and contamination from different scattering characteristics [25]. Especially in texture-rich areas, the scattering characteristics and phase variations of individual pixels vary substantially, posing greater challenges for accurate phase recovery in deformation analysis [26]. In 2011, a new TS-InSAR technique named SqueeSAR™ was proposed, which is based on DS and is designed to apply a Phase Triangulation (PT) algorithm to accurately recover phase in the filtered interferogram [27].

The DS interferometry (DSI) method has been widely applied in detecting potentially unstable slopes, landslides, and surface motion change [28,29,30,31]. However, the Kolmogorov–Smirnov (KS) hypothesis test used in the DSI method has low power, especially under small sample conditions (N ≤ 20). This may result in the inclusion of many heterogeneous points in the selected homogeneous set [2]. The generalized likelihood ratio (GLR) test is used to calculate and determine whether there is a significant difference in variance between the two sample data, thereby identifying homogeneous pixels and more accurately evaluating the significance difference between models. However, the GLR test still has low power issues when dealing with large sample sizes [32]. Recently, based on the principles of DSI, some extended homogeneous selection algorithms have been developed, such as the Anderson–Darling (AD) test, adaptive test, and confidence interval (CI) estimation [26,33,34].

Another limiting factor for TS-InSAR landslide monitoring comes from the atmospheric phase screen (APS). This is mainly caused by spatial and temporal variations in pressure, temperature, and relative humidity in the lower troposphere [14]. Terrain-related stratified delays may introduce seasonal oscillation bias to the deformation time-series measured by InSAR in steep terrain [35,36]. Since landslides often occur in steep mountainous areas with harsh environments, established correction methods based on the interferometric phase, such as empirical modeling and variance covariance estimation, are usually susceptible to interference from other contaminants in the interferometric phase including topographic errors, decorrelation noise, and deformation [37,38,39]. Recent studies have shown the near real-time availability of atmospheric parameters output from global atmospheric models (GAMs), reflecting the higher development potential of GAMs in mitigating InSAR tropospheric delays [16,40,41,42]. However, modeling for atmospheric spatial and temporal variability is often neglected in the current GAMs-based correction methods, resulting in GAMs methods being limited by spatial and temporal accuracy issues. Guo et al. [43] used the PZTD-NEF model to describe the spatial and temporal variability of the troposphere. However, the nonlinear properties of the tropopause delay in the time dimension may be more complex than the expected seasonal oscillatory properties.

In this study, we proposed a combined APS correction model with DSI for deformation monitoring of mountain landslides to overcome the effects of tropospheric delay and decoherence noise. First, we used the computationally efficient fast SHP selection (FaSHPS) test on the full-resolution dataset. Furthermore, we used a temporal phase optimization method to estimate the optimal interferometric phase from the complex coherence matrix of each DS target. Subsequently, we constrained the tropospheric delay reference point to the ground, modeled the negative exponential atmospheric vertical profile using the elevation decline factor and the surface zenith total delay (ZTD), and then introduced the Prophet time-series prediction model to fit the time-dimensional variations of the model factors to model the time-series negative exponential function ZTD correction (Prophet-ZTD-NEF), aiming to mitigate the bias in tropospheric delay estimation caused by atmospheric spatial and temporal variability. Finally, we integrated these methods into the standard StaMPS framework for 3D unwrapping, atmospheric delay estimation, and time-series estimation. We used archived SENTINEL-1A single look complex (SLC) data covering from 2022 to 2023 to study the landslide of the Lancang River Basin in Deqin Tibetan Autonomous Prefecture, China to verify the feasibility of the proposed method.

2. Materials and Methods

In this study, we proposed a combined APS correction model and DSI method for deformation monitoring of the Riwagong landslide. The technical flowchart is illustrated in Figure 1. The open source and interoperable of StaMPS allowed us to extend Hooper by adding homogeneous point selection and temporal phase optimization steps [44]. Considering the complex weather conditions in mountainous areas, we also integrated the constructed Prophet_ZTD-NEF model into the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) to mitigate APS error from each interferogram. The entire data processing involved the combined use of GAMMA, StaMPS, and TRAIN. We will explain the steps shown in Figure 1 in detail in the following section.

2.1. Statistically Homogeneous Pixels and Temporal Phase Optimization

We then used the GAMMA to StaMPS program to obtain a total of N SAR stacks registered and N − 1 interferogram stacks. The interferometric phase of the $i\mathrm{t}\mathrm{h}$ image at pixel x can be written as

$${\mathsf{\phi }}_{x,i}={\mathsf{\phi }}_{t_B,x,i}-{\mathsf{\phi }}_{t_A,x,i}={\mathsf{\Delta }\mathsf{\phi }}_{def,x,i}+\mathsf{\Delta }{\mathsf{\phi }}_{atm,x,i}+\mathsf{\Delta }{\mathsf{\phi }}_{orb,x,i}+\mathsf{\Delta }{\mathsf{\phi }}_{dem,x,i}+{\mathsf{\Delta }\mathsf{\phi }}_{noi,x,i}$$

(1)

where $\mathsf{\Delta }{\mathsf{\phi }}_{def,x,i}$ denotes the phase change due to the movement of the image element in the satellite line-of-sight (LOS) direction. $\mathsf{\Delta }{\mathsf{\phi }}_{atm,x,i}$ denotes the atmospheric delay component from $t_B$ to $t_A$, $\mathsf{\Delta }{\mathsf{\phi }}_{orb,x,i}$ is the phase contribution due to orbital inaccuracies. ${\mathsf{\Delta }\mathsf{\phi }}_{dem,x,i}$ denotes residual terrain phase due to DEM errors. $\mathsf{\Delta }{\mathsf{\phi }}_{noi,x,i}$ denotes the noise term due to the variability of pixel scattering and is small enough that it does not completely obscure the signal.

We assume that the first three terms in Equation (1) are correlated on a particular spatial scale. Meanwhile, the last two terms are uncorrelated, so that low-frequency signals can be removed using low-pass adaptive filtering. The residual phase of the $ith$ image at pixel x can be written as

$$W\left\{{\mathsf{\phi }}_{x,i}-{\overline{\mathsf{\phi }}}_{x,i}\right\}=W\left\{\mathsf{\Delta }{\mathsf{\phi }}_{dem,x,i}+{\mathsf{\Delta }\mathsf{\phi }}_{noi,x,i}+{\mathsf{\Delta }\mathsf{\phi }}_{\epsilon ,x,i}\right\}$$

(2)

where $W$ denote the wrap operator. ${\overline{\mathsf{\phi }}}_{x,i}$ denotes the phase after low pass filtering of ${\mathsf{\phi }}_{x,i}$. $\mathsf{\Delta }{\mathsf{\phi }}_{dem,x,i}$ and ${\mathsf{\Delta }\mathsf{\phi }}_{noi,x,i}$ denote the non-space related parts, respectively. ${\mathsf{\Delta }\mathsf{\phi }}_{\epsilon ,x,i}$ can be largely excluded, since only the non-spatially relevant parts of the first three terms of Equation (1) are represented.

For the selection of PS targets, we first identify PS candidate points based on the amplitude dispersion index $D_A$ [19].

$$D_A={\displaystyle \frac{{\sigma }_A}{{\mu }_A}}\cong {\sigma }_{\phi }$$

(3)

where $A$ denotes amplitude, ${\sigma }_A$ and ${\mu }_A$ are the standard deviation and mean of a series of amplitude values, respectively. ${\sigma }_{\phi }$ indicates phase standard deviation (STD).

The measure of temporal coherence $r_x^{PS}$ is used to analyze the stability of the phase of each PSCs and serves as an indicator of whether a pixel is a PS or not [44,45].

$$r_x^{PS}={\displaystyle \frac{1}{N}}\left|{\sum }_{i=1}^{N}exp\left\{j\left({\mathsf{\Phi }}_{x,i}-{\overline{\mathsf{\Phi }}}_{x,i}-\mathsf{\Delta }{\widehat{\mathsf{\Phi }}}_{h,x,i}\right)\right\}\right|$$

(4)

where $\mathsf{\Delta }{\widehat{\mathsf{\Phi }}}_{h,x,i}$ denotes the interferogram corrected for flat terrain. $j$ denotes an imaginary number. Finally, we consider both ${\mathrm{D}}_A$ and $r_x^{PS}$ to determine the final PS targets.

Further, we selected DS based on confidence interval (CI) estimation. For the reference pixel $p$ and the target pixel $q$, according to the definition of the standard normal distribution, the average amplitude ${\overline{A}}_q={\displaystyle \frac{1}{N}}\sum _{i=1}^N{A}_{i\left(q\right)}$ tends to be normally distributed as the number of samples $N$ increases, and then the CI of ${\overline{A}}_p$ is expressed as [26]:

$$\mathrm{P}\left\{u_p-z_{1-{\displaystyle \frac{\alpha }{2}}}·\sqrt{{\displaystyle \frac{Var\left(A_p\right)}{N}}}<{\overline{A}}_q<u_p+z_{1-{\displaystyle \frac{\alpha }{2}}}·\sqrt{{\displaystyle \frac{Var\left(A_p\right)}{N}}}\right\}=1-\mathsf{\alpha }$$

(5)

where $\mathrm{P}\left\{·\right\}$ denotes the probability, $z_{1-{\displaystyle \frac{\alpha }{2}}}$ denotes the distribution point corresponding to $\mathsf{\alpha }$ in the normal distribution,$A_p$ denotes the amplitude value of $p$, and $u_p$ denotes the expected value of $p$.

Then, the average amplitude $A_p$ of the image pixel to be detected is compared with the interval by setting a small $\alpha$ in set Ω. If it is within the interval, it means that the target pixel $q$ is a homogeneous image element of the reference pixel $p$.

Since DS can be affected by decorrelation noise, it is necessary to optimize the phase of DS. This is unlike conventional spatial and temporal dimensional filtering, which improves the signal-to-noise ratio at the cost of reduced spatial resolution [46]. Here, we used a temporal phase optimization algorithm based on the principle of PT to find the phase optimal solution under the condition of satisfying PT [47]. Three covariant single look interferograms (L = 1) can be written as

$${\phi }_{nm}=W\left\{{\phi }_{om-}{\phi }_{on}\right\}=W\left\{\left[\begin{array}{cc}1& \left.-1\right]\end{array}\right.\left[\left.\begin{array}{c}{\phi }_{om}\\ {\phi }_{on}\end{array}\right]\right.\right\}$$

(6)

where ${\phi }_{nm}$, ${\phi }_{om}$, and ${\phi }_{on}$ correspond to the interferometric phases of the three interferograms, respectively. Based on the set Ω selected, a sample covariance matrix of pixel $p$ was designed for the FaSHPs test. For the intra-temporal N-dimensional scattering vectors $\mathrm{y}\left(p\right)=\left[\left.{\mathrm{y}}_1\left(p\right),{\mathrm{y}}_2\left(p\right),\cdots ,{\mathrm{y}}_n\left(p\right)\right]\right.$ covariance matrix $\mathrm{C}$ can be estimated using the coherence matrix $\mathrm{T}$:

$$\mathrm{C}\approx \mathrm{T}=\mathrm{E}\left[\left.y{y}^H\right]\right.\approx {\displaystyle \frac{1}{N_p}}{\sum }_{y\in \mathsf{\Omega }}y{y}^H$$

(7)

where $\mathrm{E}\left[·\right]$ denotes the expectation, $N_p$ denotes the number of N homogeneous points, and $\mathsf{\Omega }$ denotes the set consisting of $N_p$ homogeneous targets. $H$ refers to the Ermitian conjugate.

We estimated the optimal phase sequence $\mathsf{\theta }=\left[{\mathsf{\theta }}_1,{\mathsf{\theta }}_2,\cdots ,{\mathsf{\theta }}_n\right]$ for each pixel through the coherence matrix. The core theory of temporal phase optimization is a maximum likelihood estimation (MLE) under a multidimensional complex circular Gaussian distribution and considers $N\left(N-1\right)/2$ interferometric phase information. Therefore, to reduce the computational complexity, the first value of $\mathsf{\theta }$ is usually set to zero (${\mathsf{\theta }}_1=0$). Then, we can estimate the objective function $\widehat{\mathsf{\theta }}$ of the PT algorithm through MLE.

$$\widehat{\mathsf{\theta }}=\underset{\mathsf{\theta }}{\arg \mathit{max}}\left\{{\mathsf{\Lambda }}^H\left(-{\left|T\right|}^{-1}\circ T\right)\mathsf{\Lambda }\right\}$$

(8)

where $\mathsf{\Lambda }=\left[0,e^{j{\mathsf{\theta }}_2}\cdots ,e^{j{\mathsf{\theta }}_n}\right]$, $\circ$ represents the mathematic operator of the Hadamard product.

We evaluated the quality of DS phase optimization and the extent to which DS is affected by time decoherence by calculating the fit ${\gamma }_{DS}$ between before and after phase optimization

$${\gamma }_{DS}={\displaystyle \frac{2}{N(N-1)}}Re\left({\sum }_{m=1}^N{\sum }_{n=m+1}^N{e}^{j\left({\phi }_{n,m}-\left({\phi }_m-{\phi }_n\right)\right)}\right)$$

(9)

where ${\gamma }_{DS}$ represents the goodness of fit index. ${\phi }_{n,m}$ represents the interferometric phase before optimization. ${\phi }_m$, ${\phi }_n$ refer to the optimized interferometric phase, respectively. $Re\left(·\right)$ denotes the real part of the complex data.

Finally, we selected those candidates with ${\gamma }_{DS}$ values above a predefined threshold as the final optimized DS targets.

2.2. Prophet_ZTD-NEF APS Correction Model

The APS in interferograms is caused by temporal and spatial variations in temperature, pressure and relative humidity of the troposphere between two repeated SAR observations. After removing orbital errors, topographic residuals, and decoherence noise, the InSAR phase difference observation ${\mathsf{\Delta }\phi }_{obs}$ can be expressed as

$${\mathsf{\phi }}_{x,i}={\mathsf{\phi }}_{t_2,x,i}-{\mathsf{\phi }}_{t_1,x,i}={\mathsf{\Delta }\mathsf{\phi }}_{def,x,i}+{\mathsf{\phi }}_{atm,x,i}^{t_2}-{\mathsf{\phi }}_{atm,x,i}^{t_1}$$

(10)

where ${\mathsf{\Delta }\phi }_{def,x,i}$ represents the displacement component of the radar LOS direction from $t_2$ to $t_1$; ${\mathsf{\phi }}_{atm,x,i}^{t_1}$ and ${\mathsf{\phi }}_{atm,x,i}^{t_2}$ denotes the atmospheric delay component at $t_1$ and $t_2$, respectively.

We calculated the zenith total delay (ZTD) on different layers by integrating atmospheric parameters based on the ERA5 reanalysis outputs.

$$ZTD={\int }_{h_0}^{h_{top}}{k}_1{\displaystyle \frac{P}{T}}+k_2^{\prime }{\displaystyle \frac{e}{T}}+k_3{\displaystyle \frac{\mathrm{e}}{T^2}}dh$$

(11)

where ${\rm T}$ denotes the absolute temperature; $k_1=0.776KP{a}^{-1}$, $k_2^{\prime }=0.233KP{a}^{-1}$, $k_3=3.75\times {10}^3{\mathrm{K}}^2\mathrm{P}{\mathrm{a}}^{-1}$; $P$ and $e$ denotes the extracted stratified total pressure and water vapor pressure at each grid point, respectively.

We modeled the zenith total delay negative exponential function (ZTD-NEF) in the vertical space dimension to describe the spatial variations in atmospheric delay. Then, the fitted β was calculated based on the weighted least squares algorithm in steps of one-hour intervals [48].

$${ZTD}_{APS\_model}=\alpha ·\mathrm{e}\mathrm{x}\mathrm{p}(\beta ·h)$$

(12)

where $\alpha$ is the formula coefficient; $\beta$ represents the elevation decline factor in cm/km; $h$ represents the elevation in km.

We constrain the reference point for the tropospheric delay to the surface, therefore, ZTD at different heights at any grid can be formulated as

$${ZTD}_{APS\_model}={ZTD}_r·exp(\left(h-h_r\right)·\beta )$$

(13)

where $h_r$ denotes the reference attitude; ${ZTD}_r$ denotes the reference ZTD at reference attitude.

Considering the deviation of the atmospheric parameters output by the GAMs from the moment of SAR acquisition, we implemented a time-dimensional compensation scheme for the ZTD calculated by integration instead of directly using the closest ERA5 outputs. The Prophet time-series prediction model was used to fit hourly intervals of ${ZTD}_r$ and $\beta$.

$$\left\{\begin{array}{c}{ZTD}_r\left(t\right)=g_{{ZTD}_r}\left(t\right)+s_{{ZTD}_r}\left(t\right)+h_{{ZTD}_r}\left(t\right)+{\epsilon }_{{ZTD}_r}\left(t\right)\\ \beta \left(t\right)=g_{\beta }\left(t\right)+s_{\beta }\left(t\right)+h_{\beta }\left(t\right){+\epsilon }_{\beta }\left(t\right)\\ g(t)={\displaystyle \frac{C}{1+exp(-k(t-m\left)\right)}}\\ s\left(t\right)=\sum _{n=1}^N\left(a_n\mathit{cos}{\displaystyle \frac{2\pi nt}{365.25}}+b_n\mathit{sin}{\displaystyle \frac{2\pi nt}{365.25}}\right)\\ h\left(t\right)=\sum _j\left(k_j·I\left(t\in {holiday}_j\right)\right)\end{array}\right.$$

(14)

where ${\mathrm{Z}\mathrm{T}\mathrm{D}}_{\mathrm{r}}\left(t\right)$ and $\beta \left(t\right)$ denote ${ZTD}_r$ and $\beta$ time series forecasting models, respectively. $g\left(t\right)$ denotes the saturated growth function, $C$ is the maximum asymptotic value, $k$ represents the growth rate, and $m$ denotes the offset. $s\left(t\right)$ denotes the Fourier periodic function, and $a_n$ and $b_n$ are the coefficients of $s\left(t\right)$. When n = 2, it is the PZTD-NEF method [43]. $h\left(t\right)$ denotes holidays and events, $k_j$ is the intensity parameter of the holiday effect, and $I(·)$ is the indicator function; and $\epsilon \left(t\right)$ is the normal distribution of the error term. The final Prophet-ZTD_NEF model can be combined as

$$\left\{\begin{array}{c}{ZTD}_{AP{S}_{model}}(t)={ZTD}_r\left(t\right)·exp(\left(h-h_r\right)·\beta \left(t\right))\\ {ZTD}_r\left(t\right)=g_{{ZTD}_r}\left(t\right)+s_{{ZTD}_r}\left(t\right)+h_{{ZTD}_r}\left(t\right)+{\epsilon }_{{ZTD}_r}(t)\\ \beta \left(t\right)=g_{\beta }\left(t\right)+s_{\beta }\left(t\right)+h_{\beta }\left(t\right){+\epsilon }_{\beta }(t)\end{array}\right.$$

(15)

We used the L-BFGS algorithm for model fitting, which has a higher convergence rate in solving unconstrained nonlinear programming problems [49].

Then, we established the prediction framework of tropospheric delay according to the 3D coordinates and time of each coherence point. In the same altitude plane, we used the ZTD outputted by Prophet-ZTD_NEF and inverse distance weighting strategy to realize the accurate calculation of the ZTD at any point to be measured [43].

Further, we converted the ZTD to an atmospheric delay phase based on the pseudorange conversion factor and the local radar incidence angle, which is used to correct the APS effect in the interferograms

$${\phi }_{AP{S}_{model}}(t)={\displaystyle \frac{-4\pi }{\lambda }}{10}^{-6}{\left(\cos \theta \right)}^{-1}·{ZTD}_{AP{S}_{model}}(t)$$

(16)

where $\theta$ represents the radar incidence angle; $\lambda$ represent the radar wavelength; ${\displaystyle \frac{-4\pi }{\lambda }}$ represent a conversion factor to convert from pseudo-range increase to phase delay [14].

Finally, we implemented the time-series analysis with the StaMPS toolkit. The DEM errors were first corrected by subtracting the estimated values of $\mathsf{\Delta }{\widehat{\mathsf{\Phi }}}_{h,x,i}$, and after correcting DS points for spatially uncorrelated look angle errors. Then, the 3D unwrapping algorithm was used to unwrap the phase of each combined interferogram to solve for the time-series phase of each pixel [50]. A summary of the applied StaMPS processing parameters is shown in

Table 1. Key parameters set in DS time-series analysis.

Parameter	Value	Parameter	Value
Max topo error	35	Weed time win	365
Select method	Percent	Unwrap method	3D
Percent rand	20	GAMMA iterations	6
Weed max noise	2	Unwrap grid size	25

. By separating the above error terms, the displacement phase on the radar LOS is obtained.

We compared the PZTD-NEF and Prophet_ZTD-NEF model performance using the ${ZTD}_r$ of grid point 28°N, 99°E as a reference value for validation, and plotted the correlation diagram as shown in Figure 2. The disadvantage of the PZTD- NEF model is that the model-estimated ${ZTD}_r$ may be underfitted where the atmosphere varies by a large magnitude over a short period of time. However, the Prophet_ZTD-NEF method has a better fitting performance for the non-linear variation of the model coefficients ${ZTD}_r$, reducing the root mean square error (RMSE) to the real value by approximately 24.7% compared to the PZTD-NEF method. In total, 100 and 500 samples were randomly selected and the RMSE of Prophet_ZTD-NEF with respect to the real value was 1.19 and 1.40 cm instead of 1.58 and 1.72 cm (PZTD-NEF), respectively.

3. Test Sites and Datasets

3.1. Geological Setting

Diqing Tibetan Autonomous Prefecture is located in the transverse valley region on the southeastern edge of the Qinghai–Tibet Plateau. The region is traversed by three major rivers, namely the Nujiang River, Lancang River, and Jinsha River, along with their main and tributary streams, making it the central area of the world-renowned landscape known as the Three Parallel Rivers as shown in Figure 3a. One of the deepest canyons in the world, Tiger Leaping Gorge, and the highest peak in northwest Yunnan, Meili Snow Mountain, are both located in this region. The terrain exhibits a stepped distribution, with the elevation gradually decreasing from northwest to southeast with the high mountain and canyon landforms are the main topographic features of the area.

The surrounding area has complex geological structures with a high concentration of active faults, particularly due to the uplifting effect of the geotectonic movements in the Himalayan mountain range, which has resulted in frequent seismic activity ranging from moderate to strong earthquakes (4~5.9 magnitude). The geological formations are primarily composed of Mesozoic and Cenozoic rocks, including sedimentary, metamorphic, and igneous rocks. Due to the strong development of tectonics at all levels, coupled with the erosive effect of rivers, resulting in broken and loose rock layers. This makes the region highly susceptible to landslide disasters under the influence of geotectonic movements and seismic events.

The Riwagong landslide is located approximately 20 km from the urban area of Deqin County in Diqing Tibetan Autonomous Prefecture, as shown in Figure 3b. The local climate is characterized as a cold temperate mountainous monsoon climate, influenced by the alternating Southwest Monsoon and the southern branch of the Westerlies. The average annual temperature ranges from 6.3 to 11.8 °C, and the region experiences distinct wet and dry seasons. Due to the topographical difference of being lower in the south and higher in the north, the southern part experiences relatively abundant precipitation and a humid climate. Meanwhile, the northern part is relatively dry, creating a unique climate pattern known as “one mountain with four seasons, ten miles with different weather.” We collected daily variations of PWV in the European Centre for Medium-Range Weather Forecasts (ECMWF) from 2013 to 2023, as shown in Figure 3c. We found that PWV exhibits seasonal variation patterns, with peak values occurring between May and October. We anticipate that continuous heavy rainfall during the period from May to October is the main factor triggering landslides.

3.2. Datasets

We collected SENTINEL-1A SLC data using 28 descending mode acquisitions from January 2022 to January 2023, obtained from the European Space Agency. The imaging time was approximately 23:14 UTC, covering the area indicated by the green box in Figure 3a. We used a 1-arcsecond SRTM DEM for input to simulate the topography phase. Based on PSI, SBAS, and DSI methods, we combined both spatial and temporal baselines, as shown in Figure 4. The ERA5 raw output grids corresponding to all SAR acquisition times were downloaded from ECMWF and used for modeling and estimating the atmospheric conditions with the ERA5 grid locations shown as black circles in Figure 3a. Additionally, we collected precise orbital data to correct inaccurate orbital errors.

In response to the call for disaster prevention and control, more than 200 GPS monitoring stations were installed in the potential landslide areas and on the surfaces of residential buildings in the northwest region of Yunnan Province. We collected GPS monitoring data covering the research area from 1 May to 31 October 2022, as shown in Figure 3b. Two GPS base stations were installed on a stable bedrock near the GPS monitoring stations at approximately 10 m. All the GPS base stations used the L1 frequency band with a center frequency of 1575.42 MHz. The GPS monitoring data were processed using GAMIT/GLOBK v10.71 software to obtain three-dimensional information, with a sampling interval of 30–60 min.

4. Results

4.1. DSI Method Performance Evaluation

For SHP identification, the traditional DeSpecKS may be unreliable for small sample sizes. Meanwhile, it increases the processing cost for large sample data, which poses a challenge to achieving fast landslide point retrieval applications [35]. The GLR test was shown to give satisfactory results even for small samples ($10<\mathrm{N}<20$). However, for large samples, the computational efficiency is still not as high as parametric testing methods that require only logical operations. In contrast, our DSI method uses FaSHPs to test for homogeneous points. This significantly improves the computational efficiency compared to the GLR test and has a greater application potential for achieving rapid monitoring of short-term landslides. Another outstanding advantage of FaSHPs is that all temporal pixels at the same location can be estimated by the same spatial neighbors. This is crucial for analyzing the change in properties on the temporal pixels [26]. To verify the reasonableness of the FaSHPs test method for SHP identification, we used 28 SENTINELl-1A images as experimental data and tested GLR and FaSHPs at different stacks N = 10, 16, 22, and 28, respectively. We set an estimation window of 11 × 21 and a significance level of 5% and selected the riverbanks as homogeneous regions, with blue indicating the reference pixel and green indicating the homogeneous pixels. The background shows the mean SAR intensity image, as shown in Figure 5.

From Figure 5, we can observe that the GLR test is more sensitive to the number of stacks and cannot distinguish the heterogeneous image elements well in small samples (e.g., N = 10 and 16). The accuracy of homogeneous point recognition can only be guaranteed when the number of samples is sufficient. Meanwhile, the homogeneous points identified by the two tests tend to be similar with the increase in the number of stacks. However, the FaSHPs test still produced satisfactory results even in small samples.

The basic idea of temporal phase optimization is based on the PT algorithm with a MLE based on spectral estimation theory. Compared with the traditional Boxcar spatial filtering, the DSI method considers the temporal dimension and achieves an optimal estimation of the phase using the full information of the N(N − 1)/2 interferometric pairs. This is not limited to the combination of the short baseline subset of interferograms. Figure 6 shows the original, adaptively filtered (Goldstein filter) and optimized interferometric phases (SAR images from 10 May 2022 and 28 April 2022). In the low-coherence region with rugged terrain and sparse vegetation cover, the original phase contains a large amount of decoherence noise, which almost completely masks the interferometric fringes. Both temporal phase optimization and adaptive filtering can mitigate the noise. However, adaptive filtering performs better, with significant suppression of the decoherence noise, resulting in more coherent interferometric fringes and improved accuracy of the subsequent phase unwrapping.

4.2. APS Corrected Time-Series InSAR

After effective suppression of decoherence noise in the original interferograms, we estimated the APS for each SAR acquisition using the Prophet_ZTD-NEF model. We compared it with other GAMs-based atmospheric delay correction methods. The GAMs correction method has the advantage of being independent of the interferometric phase and is, therefore, suitable for use in PSI, SBAS, or combined methods.

4.2.1. Overall Assessment of APS Improvement

We first assessed the accuracy of the APS predicted by the Prophet_ZTD-NEF model. We simulated the surface atmospheric delay distribution for 10 January 2022 and 21 July 2022 at UTC = 23:00 using the DEM covering the study area as input data, as shown in Figure 7. The Prophet_ZTD-NEF method is in strong agreement with the atmospheric delays calculated by the generic atmospheric correction online service (GACOS) and PZTD-NEF models. For simulated atmospheric delays under different seasonal conditions in winter and summer, the RMSE of the Prophet_ZTD-NEF calculated values versus the GACOS and PZTD-NEF methods were 0.98 cm, 0.15 cm, and 3.21 cm, 1.96 cm, respectively. These three GAMs-based methods have in common that they all use a negative exponential function to depict the vertical profile in tropospheric delay. However, both our method and PZTD-NEF constrain the model reference to the land surface. Therefore, there is little difference between the two methods. Meanwhile, GACOS constrains the model reference to the top of the troposphere where the change in atmospheric delay is essentially zero. The reason we do this is to consider the fact that changes are more pronounced in the lower troposphere. Therefore, constraining the model reference to the surface provides a more accurate estimate of atmospheric delay changes above the land surface. The consistency of the atmospheric delays estimated by the three methods is also higher because Diqing is situated in a high mountain valley area. Here, the water vapor content of the air is very low in winter, it generally does not exceed about 5 mm, and has a small variability within a day. However, the water vapor content in the air is at its summer peak in July, and the water vapor varies more substantially throughout the day. The black box shows the difference in tropospheric delays estimated by the different methods. The GACOS method uses subjective linear interpolation in the time dimension to estimate the atmospheric delays at the moment of SAR acquisition. This may be ineffective in estimating atmospheric delays under some extreme weather conditions such as heavy rainfall. Prophet considers the trending, seasonality, and cyclicality of the data for modeling, which is more capable of capturing the time-varying characteristics of the atmosphere. Therefore, the RMSE of the atmospheric delay calculated by GACOS and Prophet_ZTD-NEF methods is also larger, reaching 3.21 cm. This also suggests that interferograms under summer conditions may be more significantly affected by APS.

We also applied three atmospheric models to correct the APS in each interferogram and systematically evaluated the phase STD changes of the interferograms before and after atmospheric delay corrections. Considering the effects of deformation signals and topographic residuals, we manually excluded interferograms with temporal baselines over 100 days and spatial baselines greater than 300 m, thus ignoring the interference of surface deformation in the short term. Finally, 169 interferogram samples were selected from the 28 collected SAR images from 10 January 2022 to 24 December 2022 and the mean phase STD for every 10 interferograms was counted (Figure 8). Overall, after three GAMs-based (GACOS, PZTD-NEF and Prophet_ZTD-NEF) methods of correction, all of them mitigated the APS in the interferograms, where our new method reduced the interferogram mean phase STD from 1.88 cm to 1.62 cm instead of 1.69 cm (PZTD-NEF) and 1.72 cm (GACOS), indicating the advantage of our new method in mitigating the spatial variability of tropospheric delays.

4.2.2. Mitigation of Topography De-Correlation Effects

In this study, all three GAMs-based atmospheric delay correction methods were modeled in the vertical spatial dimension and are, therefore, well suited to be used to assess the mitigation performance of the vertical stratified signal. In general, the larger absolute value of the correlation coefficient ($\left|\mathrm{k}\right|$) between the original interferometric phase and the local topography indicating a more severe stratified delay signal. The lower the absolute value of the corrected correlation coefficient, the more significant the mitigation of the stratified signal. We show two cases of different methods to correct the vertical stratified delay and plot the correlation between elevation and phase before and after the atmospheric delay correction, as shown in Figure 9.

As can be observed in Figure 9, due to the fact that the study area is located in a high mountain valley region, there are vertically stratified delays associated with topography in the original interferogram, especially on the mountain ranges on both sides of the canyon. The APS associated with topography is significantly mitigated after correction by the three GAM-based methods. The correlation coefficients $\left|\mathrm{k}\right|$ of the interferometric phases with topography are reduced from 5.7881 to 1.6936, 2.7917, and 0.0208 rad/km, respectively, which suggests that the three atmospheric delay correction models are effective in correcting for the stratified delays. We also analyzed the correlation coefficients of the mitigation magnitude with |k| for all interferometric phases, as shown in Figure 10. The atmospheric delay correction model is more capable of correcting interferograms the more dominant the vertical stratified delay is. Moreover, our method achieves better results with the least negative correction cases.

4.3. Accuracy Assessment of DSI Measurements

4.3.1. PSI, SBAS, and DSI Cross-Comparison

We performed TS-InSAR processing on the images covering the slopes of alpine valleys along the Lancang River Basin and calculated the deformation velocity in the LOS direction using PSI, SBAS, and DSI, respectively, as shown in Figure 11. In terms of the spatial density of observation points, the DSI method determines more MPs (>5 times) than PSI or SBAS. The spatial densities of the PSI and SBAS methods are only 51 and 92/km², respectively. This also illustrates the limitation of using conventional TS-InSAR methods applied to the densely vegetated alpine valley region, and our proposed method can overcome the decorrelation effect. The number of MPs detected by DSI is approximately eight and five times higher than that of PSI and SBAS, respectively. The total number of MPs detected by the three methods and the corresponding spatial densities are counted in

Table 2. Number of MPs and their spatial density in the results of the three methods.

Methods	Number of MPs	Spatial Density (MPs/km2)
PSI	44,307	51
SBAS	79,734	92
DSI	367,512	424

.

To assess the validity of the DSI method, we verified the consistency between the PSI, SBAS, and DSI calculated the deformation velocity. We counted all the PS and SDFP points within a 3 m radius of each other from the DS and further computed the average velocity of each PS point and SDFP point that fell within the radius. A total of 6492 PS points and 7597 SDFP were selected as the proxy for DS. Figure 12 shows the scatterplot of the cross-validation of the three methods and reports the corresponding correlation coefficients and RMSE. The highest displacement velocity correlation was monitored by PSI and DSI, with a correlation coefficient of 0.8459. The displacement velocity RMSE of approximately 6.3 mm/a, and an error between ±5 mm/a for 76.9% of the MPs. The correlation between SBAS and the DSI method was approximately 0.8371. The lowest displacement velocity correlation was monitored by PS and SBAS, with a correlation coefficient of 0.6315, and an RMSE of approximately 8.1 mm/a, with only 55.5% of the MPs having an error between ±5 mm/year. From Figure 11 and Figure 12, we can obtain the following conclusions. First, there is a discrepancy between the PSI and SBAS methods applied to the area with dense vegetation cover. The C-wavelength radar signals have a weak penetration performance to the vegetation, and we retained a part of the low-coherence points during the data processing to ensure the reliability of the unwrapping. This causes the deformation velocity solving to introduce an error. Meanwhile, it is caused by the different point selection strategies. Second, the DSIs are in high agreement with the PSI and SBAS methods. However, there are also slight differences, which are mainly due to slight differences in MP location and resolution units, and decoherence errors caused by external environments.

4.3.2. Validation of DSI Measured Deformation Velocity Using GPS Measurements

We also compared the results from DSI with GPS measurements. We selected all DS pixels in the vicinity of a horizontal distance ≤ 25 m from the GPS station, calculated the average LOS displacement velocity, and projected the LOS-directed displacements to the vertical direction by using $U_V=D_{LOS}/\cos \theta$ [41,51]. Considering the higher temporal sampling of GPS, we calculated a one-day average value to be used as a proxy for comparison with GPS measurements value, with GPS and InSAR displacements calibrated with reference to 10 May 2022 for each period, as shown in Figure 13 and

Table 3. RMSE of GPS and other atmospheric delay correction methods.

GPS Number	Original mm/yr	GACOS mm/yr	PZTD-NEF mm/yr	Prophet_ZTD-NEF mm/yr
1	5.39	5.18	4.33	4.52
2	7.69	7.59	7.87	7.57
3	3.48	3.49	3.87	3.15
4	10.03	9.90	9.91	9.15

. The RMSE of the displacement time-series of DSI (with Prophet-ZTD_NEF APS correction) measurements with G1, G2, G3, and G4 were 4.52, 7.57, 3.15, and 9.15 mm/a, instead of 5.18, 7.59, 3.49 and 9.90 mm/a by GACOS and 4.43, 7.87, 3.87 and 9.91 mm/a by PZTD_NEF. The displacements computed by the different methods had a similar time evolution trend, and the new method achieved the largest improvement in performance compared to the without APS correction, resulting in an improvement in the accuracy of the deformation velocity solution of about 8.27%.

The differences between DSI and GPS measurements are attributed to two main factors. One is the spatial resolution mismatch between DS and in situ GPS measurement. DS represents the ensemble value of MPs within the surroundings, whereas GPS is the accurate value of a single location. As shown in Figure 11, the terrain around Riwagong has a high degree of undulation due to the high mountain valley terrain. A distance of a few meters may lead to significant changes in the displacement measurements. Another key factor is the asynchronous observation of the two techniques in the time dimension. Compared with the GPS near real-time observation, the InSAR time sampling interval is 12 days, and different observation periods and different time sampling schemes will result in differences in displacement measurements. However, compared with the sparse MPs of GPS, the DSI method detects denser observation points, which is valuable for depicting the spatial morphology of landslides.

Figure 14 shows the deformation velocity estimated by different atmospheric delay correction methods. We observe a residual APS signal in the original deformation velocity map, which appears as a slight uplift deformation. After the correction of three methods, the atmospheric delay error related to terrain is alleviated to some extent. The new method yields substantially better results (Figure 14d) than the other methods, especially on the sides of the canyon.

5. Discussion

Landslide prevention and control in Deqin County has been a key concern for the government of Yunnan Province, China. Considering that almost all villages are located in the middle of narrow, deep-cut longitudinal valleys and on steep slopes, and are affected by both the Lancang River and fault zones running through the area, the houses of the residents and other infrastructures are always under threat from unstable slopes. Although the government departments are currently ready to relocate the entire county town, some remote mountain villages, such as Riwagong village in this study, still have limited government control over the area due to their dispersed population and remote mountainous areas. In the following, we focus on the Riwagong landslide detected by the DSI method, with a maximum deformation velocity of more than 110 mm/a. We analyzed the spatial distribution of the landslide deformation and discussed the influencing factors.

5.1. Spatial Distribution of the Landslide Deformation

Using DSI methods, we can detect the deformation contours of landslides more intuitively after reducing DEM errors and APS, and the spatial distribution of surface deformation matches geomorphological features well, as shown in Figure 15. However, the PSI technique only detected some PS points in some artificial structures and rocks located at the bottom of the canyon. Meanwhile, the detection deformation velocity of PS was extremely low in the high mountains to the east, as shown in Figure 11a. Therefore, the PSI is obviously not successful when applied to mountainous areas covered with thick vegetation. The small baseline strategy forms interferograms only between images with short intervals of less than 24 d and small differences in appearance and slant angle. This to some extent preserves a very low decorrelation signal and increases the density of coherent points. However, this also makes a large number of poor-quality signals in the SDFP, which is also susceptible to incorrect estimation of deformation velocity. By comparing the PSI and the SBAS method, we can observe the superiority of DSI in small-scale landslide investigations with vegetation cover. Figure 15a shows that the landslide moves more strongly in the northwest direction, with a maximum movement velocity of approximately 110 mm/a. Although the terrain is higher in the northeast direction, with an elevation of more than 3800 m, the deformation velocity is instead slower, at approximately −20 mm/a. The main reason for this is that the vegetation cover in the northeast is more luxuriant and the soil is more resistant to erosion. This curtails the occurrence of landslide creep and, therefore, the velocity of surface displacement is smaller. However, high vegetation cover may also increase slope loading, which can lead to landslide creep. The slight subsidence in the northeast is, therefore, equally worthy of attention from those living at the foot of the hillside.

We plotted the deformation velocity profiles along A–B and C–D, as shown in Figure 14. There is a response relationship between the intensity of the landslide deformation and the slope of the landslide surface, as shown in Figure 15b, c. The deformation of the landslide was more affected by the topographic relief, and the sliding surface in the northwestern part of the area was steeper than that in the northeastern part of the area. The steeper sliding surface is more likely to produce greater displacement. Therefore, the spatial distribution pattern of the surface displacement of this landslide is reasonable. Figure 16 shows the time-series of landslide LOS displacements measured using the DSI method. The temporal evolution of the landslide shows an approximately linear growth trend and reached its maximum deformation in October, with a maximum cumulative displacement of approximately 150 mm over the one-year cycle. We expect that there is a response of deformation to seasonal rainfall, which is usually concentrated in northwestern Yunnan from June to October, while the PWV in November–December 2022 is significantly lower, and the decrease in rainfall may lead to the rebound of the deformed area.

5.2. Analysis of Factors Affecting Landslide Instability

5.2.1. Effects of Geological Tectonic Movements

The Deqin–Weixi area is located in the south–central section of the Hengduan Mountains on the southeastern edge of the Tibetan Plateau. The region has a complex tectonic environment, with widely distributed canyons, mountains, rivers, and deforestation, coupled with strong exergonic denudation. This has resulted in the formation of a large area of exposed bedrock. The Riwagong landslide is located between the Lancangjiang fault and the Wixi–Qiaohou fault, and the Deqin–Zhongdian fault and the Jinshajiang fault are also distributed in the south.

Influenced by the southeastern extrusion of the Sichuan–Yunnan rhombic block, the Weixi–Qiaohou fault exhibits right-handed strike-slip characteristics, with the active nature of tensional components. Since 2013, many earthquakes of Mw5.0 or higher have occurred around this fault zone, including the Yangbi Mw6.4 earthquake that occurred on 21 May 2021 at the western boundary of the Sichuan–Yunnan block. This is also thought to have been caused by an NW-directed sub-fracture zone on the western side of the Weixi–Qiaohou fault. Figure 15 shows that the Weixi–Qiaohou fault runs through the landslide surface from north to south and is consistent with the landslide boundary defined by InSAR. The strong deformations are located on the west side of the fault zone. Therefore, the landslide may have been caused by the action of the strike-slip activity of the fracture, resulting in an anomalous stress field within the slope body, which makes it easy to form large-scale sliding.

Another fracture zone is the Lancang River Fault. When the Lancang River flows through the Hengduan Mountains, it forms a V-shaped canyon with a relative height difference of 1000–2000 m due to river erosion. River erosion is also one of the factors triggering landslide destabilization [2]. The Riwagong landslide is located in the deep-cut alpine valley area, where the topographic slope varies more than 25°. The valley area developed by the Lancang river fault is one of the most active geological and tectonic movements in China. The fragmented rock bodies formed by the long-term action of an active fault, combined with factors such as gravity, precipitation, and groundwater result in a decrease in the mechanics of the rock. Consequently, under the influence of external forces such as fault activity, earthquakes, and river erosion, slope instability is more likely to occur.

5.2.2. Response Analysis with Rainfall

Heavy precipitation can lead to an increase in pore–water pressure and reduce the effective positive stress in the basement shear zone, thus triggering landslides [52]. To understand the response of rainfall to Riwagong landslide drivers, we selected P1 and P2 MPs, located at the back and central edges of the landslide, for analysis. We compared the deformation time-series with PWV. We extracted the daily PWV value at 23:00 closest to the SAR acquisition of UTC 23:14 using the meteorological parameters provided by ERA5 and plotted the correlation between the deformation time-series of P1 and P2 and the amount of daily PWV, as shown in Figure 17.

From the deformation time-series, the first stage of the landslide occurred in April 2022, although the rainfall increased month by month from January to October. However, P1, which is located at the back edge of the slope, already showed signs of sliding in April and slid more than 50 mm in one month, which can easily cause a landslide. According to local reports, on 5 April 2022, due to continuous rainfall, the middle bridge at the exit of the Jiabi Tunnel in Yunling Township, Deqin County, was struck by falling rocks from above, causing severe damage to the bridge. The Riwagong slope, 18 km from Yunling Township may have also slid. The second phase of sliding occurred from April to October 2022. We found that the correlation between P2 deformation and precipitation was stronger, which indicates that sliding in the center was significantly stronger than at the back edge, with an average deformation velocity of 22.8 mm/month. Meanwhile, deformation at the back edge showed a seasonal oscillation pattern similar to that of precipitation. The third stage of sliding occurred from October 2022 to January 2023, which was affected by the continuous high-density rainfall in the previous period. This resulted in an increase in the water content of the soil body, dissipation of matrix suction, and a reduction in the strength of the soil body. When the shear strength of the soil body was insufficient to counteract the downward sliding force, the soil body was deformed under the action of shear stress. This, in turn, leads to the displacement of the slope body.

The seasonal snowmelt at the top of the hill may have exacerbated the landslide instability. Due to Diqing Autonomous Prefecture being located in the high alpine zone of northwestern Yunnan, the average altitude is more than 3000 m. The highest altitude of the Riwagong landslide reaches 3848 m. The top layer of the soil body of the landslide mostly freezes from November every year to March next year, and then gradually melts and thaws in March. Snow and ice melting and water seepage gradually increase the weight of the slope body, reducing the strength of the rock and soil body, and causing slope instability.

6. Conclusions

In this study, to map landslide surface displacements in a densely vegetated rural mountain environment, we extend the StaMPS/PS method aiming to jointly use both PS and DS targets to more reliably estimate ambiguous phases by increasing the spatial sampling density of observation points. In addition, we proposed the Prophet_ZTD-NEF model based on GAMs aimed at modeling APS at coherent points by considering the temporal and spatial variability of tropospheric delays. To assess the effectiveness of the method, we mapped the displacement of the Riwagong landslide in Deqin County, Diqing Autonomous Prefecture, China, using a data stack of 28 SENTINEL-1A SLC images. The following conclusions can be summarized.

First, the DSI method in this study is in strong agreement with PS and SBAS methods in landslide displacement monitoring, with correlations of 0.84 and 0.83, respectively. DSI maps the spatial displacement of the Riwagong landslide in more detail than the traditional PSI and SBAS methods and detects significantly more MPs (>5 times) than the other two methods. This effectively overcomes the key problem of untangling errors due to the sparsity of observation points.

Second, we counted the changes in SHPs under different stack samples. Overall, the GLR test is more sensitive to stack samples, and the FaSHPs test based on CI in this study has a greater advantage than the GLR test when the number of samples is limited.

Third, the Prophet_ZTD-NEF model can effectively overcome the estimation bias caused by the spatial and temporal variability of the tropospheric delays. The RMSE of the model estimates to the real values is reduced by 24.7% compared to the PZTD-NEF model. The STD of the 169 interferometric phases counted is also reduced from 1.88 cm to 1.62 cm instead of 1.69 cm (PZTD-NEF) and 1.72 cm (GACOS). The GAMs-based correction method has the advantage of being independent of the interferograms and can be applied to other TS-InSAR methods.

Fourth, by comparing the DSI results with the in situ GPS measurement, we found that the accuracy of atmospheric corrected DSI results to the SENTINEL-1A data is approximately 6.01 mm/a and improved the accuracy of the deformation velocity solution by approximately 8.27%. This suggests that the proposed Prophet_ZTD-NEF method has considerable potential to become a reliable tool to complement traditional geodetic techniques for landslide deformation.

Finally, the spatial and temporal characteristics of the surface displacement of the Riwagong landslide are analyzed in detail. The DSI results of the SENTINEL-1A data show the displacement time-series of the landslide, with a maximum LOS displacement velocity of approximately 120 mm/a. The northwest part of the landslide moves significantly faster than the southeast part, and there is a response between the landslide boundary and the distribution of the fracture zone. On this basis, a qualitative analysis of the various factors affecting the destabilization of the Riwagong landslide was performed. Heavy rainfall was the main driver of landslide deformation. Meanwhile, the difference in deformation velocity between the north and south parts was largely due to differences in the degree of vegetation cover and the distribution of the fracture zone.