Refined InSAR Mapping Based on Improved Tropospheric Delay Correction Method for Automatic Identification of Wide-Area Potential Landslides

Abstract

Slow-moving landslides often occur in areas of high relief, which are significantly affected by tropospheric delay. In general, tropospheric delay correction methods in the synthetic-aperture radar interferometry (InSAR) field can be broadly divided into those based on external auxiliary information and those based on traditional empirical models. External auxiliary information is hindered by the low spatial–temporal resolution. Traditional empirical models can be adaptable for the spatial heterogeneity of tropospheric delay, but are limited by preset window sizes and models. In this regard, this paper proposes an improved tropospheric delay correction method based on the multivariable move-window variation model (MMVM) to adaptively determine the window size and the empirical model. Considering topography and surface deformation, the MMVM uses multivariate variogram models with iterative weight to determine the window size and model, and uses the Levenberg–Marquardt (LM) algorithm to enhance convergence speed and robustness. The high-precision surface deformation is then derived. Combined with hotspot analysis (HSA), wide-area potential landslides can be automatically identified. The reservoir area of the Baihetan hydropower station in the lower reaches of the Jinsha River was selected as the study area, using 118 Sentinel-1A images to compare with four methods in three aspects: corrected interferograms, derived deformation rate, and stability of time-series deformation. In terms of mean standard deviation, the MMVM achieved the lowest value for the unwrapped phase in the non-deformed areas, representing a reduction of 56.4% compared to the original value. Finally, 32 landslides were identified, 16 of which posed a threat to nearby villages. The experimental results demonstrate the superiority of the proposed method and provide support to disaster investigation departments.

1. Introduction

Synthetic-aperture radar interferometry (InSAR) is a crucial tool for identifying potential geological hazards over wide areas, especially slow-moving landslides (<60 mm/year) [1]. Compared to ground surveys and leveling measurements, it offers advantages such as extensive coverage, high efficiency, and cost-effectiveness. Landslides are usually identified on the basis of the surface deformation rate at measurement points (MPs) [2]. After terrain and baseline errors have been mitigated for, two factors affect the identification of landslides, geometric distortions and atmospheric delay, especially in high-relief and low-coherence regions [3]. The influence of geometric distortions can be mitigated by evaluating the terrain visibility of SAR images [4,5,6]. However, the latter, characterized by spatial heterogeneity, remains the primary source of error in InSAR measurement and affects the identification of geological hazards [7]. Advanced time-series InSAR techniques such as PS-InSAR [8] and SBAS-InSAR [9] use either short spatial–temporal baselines or spatial–temporal filtering to mitigate the atmospheric delay [10]. Due to the non-uniform distribution characteristics, atmospheric delay still impacts the InSAR results [11].

The electromagnetic waves emitted by the sensor can cause phase delays during atmospheric propagation, and a 20% change in the relative humidity of water vapor in the atmosphere can lead to deformations on a centimeter scale [12]. The near-surface atmosphere, where water vapor is present, has a relatively mixed composition, including a terrain-dependent atmosphere in the vertical direction and a terrain-independent atmosphere in the horizontal direction, i.e., atmospheric turbulence [13]. Atmospheric delays are characterized as temporally uncorrelated but spatially correlated. Therefore, on a spatial scale of a few kilometers, they can be mitigated by temporal high-pass and spatial low-pass filters [8]. In addition, the atmospheric delay can also be reduced by stacking the unwrapped phase [14]. Stratified tropospheric delay can introduce significant bias into InSAR-derived deformation, especially in regions of high relief [15]. The path-dependent nature of the hydrostatic delay in the troposphere results in a strong correlation between the unwrapped phase and local terrain [16].

For decades, numerous tropospheric delay correction methods have been proposed, which can be broadly divided into two categories: those based on external auxiliary information and those based on phase information of the interferogram. The methods based on external auxiliary information can be further divided into methods based on meteorological data and methods based on measured data. The use of meteorological reanalysis data, including NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) in combination with the fifth generation of the European Centre for Medium-Range Weather Forecasts Reanalysis (ERA5) and the incorporation of meteorological models, has shown superior performance in many cases [17,18]. The release of the Generic Atmospheric Correction Online Service for InSAR (GACOS) facilitated the process of tropospheric delay correction [19,20] and it shows excellent performance in mining areas, seismic zones, and other challenging cases [21,22]. However, the low spatial–temporal resolution and sensitivity to cloud cover limit the applicability of these methods. In addition, these methods do not analyze the spatial–temporal characteristics of the atmosphere, making it difficult to model an atmosphere characterized by randomly distributed atmospheric turbulence. Measurement-based methods allow for the retrieval of atmospheric water vapor in a given region, offering the potential to correct the tropospheric delay [23]. However, the sparse spatial distribution, such as from the Global Navigation Satellite System (GNSS) network [24], affects the effectiveness of the method.

The phase-based correction methods, known as the empirical model approach, use the relationship between the unwrapped phase and the local terrain [25,26]. Typically, this method involves window segmentation and blockwise processing of the image. Through a simulation experiment and the actual case, a window size of 5 km was chosen to estimate the tropospheric delay with a linear model [27]. Adaptive partitioning based on the elevation variation in the covered area, and the experimental results in the Hawaiian Islands demonstrated the effectiveness of the method [12]. Based on a quadtree method, the elevation gradient was used to adaptively determine the window size [26,28]. The adaptive method offers advantages in mining areas and regions with crustal activity. However, the manually set threshold for elevation gradient and the heavy computational burden need to be improved. In addition, at scales smaller than 10 km, the degree of correlation between phase and topography becomes highly heterogeneous, influenced by local topographic relief such as over a single valley or ridge in mountainous areas [29]. Therefore, the window size should be set at a relatively large spatial scale [27]. For example, at small scales, windows may be completely occupied by atmospheric turbulence, reducing the reliability of the estimated results, and increasing the computational burden [30]. On the contrary, if the scale is too large, the topography within the window may have a small elevation gradient, which is insufficient to establish the relationship model [26].

Furthermore, several empirical models have been developed to describe the relationship between phase and terrain in non-deformed areas. A linear model [31], an exponential model [32], and a variable power law model have been adopted to account for the spatial variability in atmospheric properties [30]. Advanced empirical models have become a trend [33]. The above models offer great effects in some cases but are not suitable for all regions [34]. The correction method involving a single empirical relationship is not robust for effective correction [27,28]. And spatial clustering methods, such as hotspot analysis (HSA), automatically identify areas of anomalous deformation [35,36]. Given that the better the correction quality, the more accurately the at-risk zone can be identified, it can be employed to evaluate the performance of the correction model. In the meantime, the HSA facilitates the identification of potential at-risk zones over wide areas.

Considering the existing challenges and previous research, we propose an improved tropospheric delay correction method based on the multivariable move-window variation model (MMVM), combined with HSA, for automatic identification of wide-area potential landslides. Taking topographic factors and surface deformation into consideration, the MMVM uses multivariate variogram models with iterative weight derived from Stacking-InSAR to adaptively determine the window size and the model based on the distribution characteristics between the unwrapped phase and the local terrain, and uses the Levenberg–Marquardt (LM) algorithm to improve both the convergence speed and robustness. The proposed method was applied to the lower reaches of the Jinsha River above the Baihetan Hydropower Reservoir, using a total of 118 Sentinel-1A images. A comparative analysis was carried out with the exponential model (Exp), ERA5, GACOS, and the traditional SBAS method in terms of corrected interferograms, derived deformation rate, and stability of time-series deformation to demonstrate the superiority of the proposed method.

The remainder of the article is organized as follows. Section 2 introduces the study area and datasets. Section 3 describes the principle and details of the proposed method. Section 4 presents experimental results and analysis. The superiority and limitations of the proposed method are discussed in Section 5. Section 6 concludes the paper.

2. Study Area and Datasets

2.1. Study Area

The topography of a typical alpine valley, characterized by steep ravines and widespread hazardous rocks, on both sides of the Jinsha River makes it susceptible to geological hazards, as shown in Figure 1. Other factors, such as steep slopes, sparse vegetation, and erosion from rivers, have made the situation even more critical. It is also an excellent place to compare the performance of tropospheric delay correction methods, because the temperatures in the region have a distinct vertical distribution with the local terrain [37]. Numerous geological hazards have occurred along the Jinsha River, exemplified by cases such as Wangjiashan [38], Tuandigou [39], Wulipo [40], and others [41,42], and have attracted widespread attention. The study area is located at the junction of Sichuan, Yunnan, and Guizhou provinces in China, in the lower reaches of the Jinsha River, covering an area of 32,818 km². The reservoir area extends from Leibo, Sichuan province, in the north to Dongchuan, Yunnan province, in the south, covering a length of 270.24 km along the Jinsha River. The area includes the Baihetan hydropower station, which is the second-largest hydropower station in China after the Three Gorges hydropower station [40].

2.2. Datasets

A total of 118 Sentinel-1A images were used to cover the study area with a resolution of 2.33 m in range and 13.99 m in azimuth. The incidence angle is 33.92°. The coverage period of the data is from 11 January 2020 to 12 May 2022. The reference image of the dataset is selected as 17 January 2021. The spatial–temporal baselines of interferogram pairs are shown in Figure 2, with an average temporal baseline of 24.65 days and an average spatial baseline of 64.04 m. The AW3D30 DEM is used in the standard differential InSAR processing [43].

3. Methodology

In this study, we utilize the SBAS-InSAR, incorporating an improved tropospheric delay correction method, to derive high-precision surface deformation information of the study area. Based on the hotspot analysis and spatial clustering, the potential landslides over a large area can be automatically identified. The workflow of the proposed method is shown in Figure 3, which can be divided into the following three parts.

3.1. InSAR Processing

In this study, two types of InSAR time-series techniques are involved: SBAS-InSAR and Stacking-InSAR. Based on the consideration of the characteristics of the phase components, the deformation derived from SBAS-InSAR using iterative regression analysis is more accurate than that from Stacking-InSAR [4]. On the one hand, the deformed signal within the segmentation window affects the relationship modeling. On the other hand, to avoid atmospheric overcorrection in the deformed areas, it is essential to distinguish between deformed and non-deformed areas. Stacking-InSAR is a fast deformation inversion technique [44]. By setting an appropriate threshold, the deformation rate derived from Stacking-InSAR can be used to distinguish between the deformed and non-deformed areas. The mask of the deformed areas can be updated synchronously with tropospheric-delay-corrected interferograms for each optimization iteration. Therefore, this study primarily utilizes the SBAS-InSAR workflow to obtain the accurate deformation information, including the deformation rate and cumulative deformation. The interferograms generated by SBAS-InSAR are further corrected by the proposed correction method. The mask of the deformed areas generated by Stacking-InSAR is added to the correction process and updated synchronously with the iteration. As a result, the MMVM-based method produces a higher quantity and quality of MPs that refine the InSAR deformation map.

3.1.1. SBAS-InSAR Processing

The standard differential InSAR processing includes the registration and resampling of SAR images and multi-look processing [8]. Based on the multiple reference images, interferograms with short spatial–temporal baselines are generated to reduce the decorrelation effect. External DEM can be used to remove the flat earth and topographic phases. Baseline error can be mitigated by using a quadratic polynomial model [45]. Singular value decomposition (SVD) is applied to obtain the linear deformation vector based on the minimum solution criterion. After the SVD operation, the residual unmodeled phase still contains phase terms related to the atmospheric phase delay, baseline errors, nonlinear deformation, and signal noise. Given the different characteristics in the spatial–temporal domain, operations such as temporal high-pass filtering, spatial low-pass filtering, and iterative optimization are performed to generate the deformation information [4]. Finally, high-precision deformation information of the study area is obtained, including the annual deformation rate and the time-series of cumulative deformation.

3.1.2. Stacking-InSAR Processing

Stacking-InSAR becomes a fast deformation inversion technique by assuming that the atmospheric effect is independent for each observation and randomly distributed over time [46]. Assuming there are m interferograms, the annual deformation rate $v_{stg}$ derived from Stacking-InSAR [47] can be expressed as

$$v_{stg}=\sum _{i=1}^m\left[-\frac{\lambda }{4\pi }({\varphi }_{diff}^i-{\varphi }_{diff,r}^i)\Delta T_i\right]/\sum _{i=1}^m\Delta T_i^2$$

(1)

where $\lambda$ is the radar wavelength, $\Delta T_i$ is the acquired time interval of the ith interferogram, ${\varphi }_{diff}^i$ is the unwrapped phase of the arbitrary point in the ith interferogram, and ${\varphi }_{diff,r}^i$ is the unwrapped phase of the reference point in the ith interferogram. The reference point can be selected in areas with high and stable coherence over time, such as the central urban region [48,49].

3.2. Improved Tropospheric Delay Correction Method

A fixed window size and single empirical relationship are not suitable for the spatial variation in the heterogeneous atmosphere [32]; their selection should be based on the characteristics of the atmospheric conditions. To address these issues, multivariate variogram models are introduced to adaptively determine the optimal window size. And, based on the optimal window size, taking into account the topographic factor and surface deformation, the MMVM uses multivariate variogram models with iterative weight to estimate the tropospheric delay phase.

3.2.1. Adaptive Window Selection

A window size of 5 km was chosen for block partitioning, which performed well in the experimental case [27]. However, it is challenging to adapt it to other cases, as well as the relationship model. Adaptive partitioning based on elevation changes is computationally intensive, although it has performed well in mining areas [28]. Multivariate variogram models can be used to solve this problem. Note that deformed and low-quality signals, such as those in non-visible and low-coherence areas, can affect the model. A coherence threshold is set to exclude low-quality signals. And, the annual deformation rate $v_{stg}$ derived from Stacking-InSAR [44], is used to distinguish between the deformed and non-deformed areas with a threshold, typically set at 5 mm/year [27]. $R_{index}$, a terrain visibility evaluation model based on the cosine of the angle between the local terrain surface and the radar beam, is calculated to analyze the visibility of the study area [5]:

$$R_{index}=sin\{{\alpha }_{inc}+arctan[tan{\alpha }_{slp}cos({\alpha }_{azm}-{\alpha }_{asp})]\}$$

(2)

where ${\alpha }_{inc}$ is the radar incidence angle, ${\alpha }_{slp}$ is the slope of the terrain, ${\alpha }_{azm}$ is the radar azimuth angle, and ${\alpha }_{asp}$ is the aspect of the terrain. After masking the unreliable area, the baseline errors for the reliable point $(x,y)$ are corrected using precise orbit information and a quadratic polynomial model [45]. The residual unwrapped phase $\tilde{\varphi }$ can be expressed as

$$\tilde{\varphi }=\varphi \left(X\right|x,y,h)+{\varphi }_e$$

(3)

where h is the elevation, ${\varphi }_e$ includes Gaussian white noise, decorrelation noise, system noise, etc., $\varphi \left(X\right|x,y,h)$ is a terrain-dependent phase, and A requires that the segmentation window be adapted to account for the spatial variability in the tropospheric signal. Assuming that c represents the overlap of the window and $\Delta {\mathrm{D}}_{pixel}$ represents the actual range of one pixel on the interferogram, the relationship between the number of windows n and the window size $S\left(n\right)$ can be given by

$$\begin{array}{c}S\left(n\right)=\frac{A}{n-nc-c}\hfill \\ n\sim [2,floor(\frac{A+{10}^{4}c/\Delta {\mathrm{D}}_{pixel}}{{10}^4(1-c)/\Delta {\mathrm{D}}_{pixel}})]\hfill \end{array}$$

(4)

where A is the number of pixels in the azimuth and range directions of the interferogram. When n is 1, overlap is not considered, hence c is 0 and $S\left(1\right)$ is A. And, if the number of reliable points within the window $S\left(n\right)$ is k, Pearson’s correlation coefficient R between the observation ${\varphi }_{S\left(n\right)}$ and the fitting value ${\stackrel{⌢}{\varphi }}_{S\left(n\right)}$ can be expressed as

$$R=\frac{1}{k-1}\sum _{i=1}^k\left(\frac{{\varphi }_{S\left(n\right)}^i-\overline{{\varphi }_{S\left(n\right)}}}{\sigma \left({\varphi }_{S\left(n\right)}\right)}\right)\left(\frac{{\stackrel{⌢}{\varphi }}_{S\left(n\right)}^i-\overline{{\stackrel{⌢}{\varphi }}_{S\left(n\right)}}}{\sigma ({\stackrel{⌢}{\varphi }}_{S\left(n\right)})}\right)$$

(5)

where $\overline{(\ast )}$ represents the average, and $\sigma (\ast )$ represents the standard deviation. Under the condition that $S\left(n\right)$ is greater than 10 km, the optimal scale at which the correlation reaches its maximum is the optimal size of the window:

$$S{\left(n\right)}_{opt}=argmax\left(R\right)$$

(6)

where $argmax(\ast )$ represents the value of the argument at which the given function reaches its maximum. Multivariable variogram models, e.g., linear, exponential, spherical models, etc., are used to assess the degree of fit between the unwrapped phase and the local terrain at different spatial scales within the reliable areas.

$$\begin{array}{cc}\hfill X& =argmin\left\{\mathcal{F}({\varphi }_{S\left(n\right)}-{\stackrel{⌢}{\varphi }}_{S\left(n\right)})\right\}\hfill \\ & \hfill =argmin\left\{{\sum _{i=1}^k\left[\begin{array}{c}{\varphi }_{S\left(n\right)}^i\\ {\varphi }_{S\left(n\right)}^i\\ {\varphi }_{S\left(n\right)}^i\\ ⋮\end{array}-\left(\begin{array}{c}{x}_1·h_{S\left(n\right)}+x_2\\ x_3·{x_4}^{h_{S\left(n\right)}·x_5}+x_6\\ x_7·h_{S\left(n\right)}^2+x_8·h_{S\left(n\right)}+x_9\\ ⋮\end{array}\right)\right]}^2\right\}\end{array}$$

(7)

where $argmin(\ast )$ represents the value of the argument at which the given function reaches its minimum, $\mathcal{F}(\ast )$ represents the error function, x represents the unknown parameters of the variogram model, and $h_{S\left(n\right)}$ represents the elevation within the window $S\left(n\right)$. The LM algorithm [50], an optimization algorithm with the advantages of both gradient descent and Newton’s method, is used to solve for the parameters required by the model, based on the criterion of minimizing the following objective function:

$${\left(\begin{array}{c}J\left(X\right)\\ \mu I\end{array}\right)}^T\left(\begin{array}{c}J\left(X\right)\\ \mu I\end{array}\right)\Delta X={\left(\begin{array}{c}J\left(X\right)\\ \mu I\end{array}\right)}^T\left(\begin{array}{c}J{\left(X\right)}^T\mathcal{F}\left(X\right)\\ 0\end{array}\right)$$

(8)

$$X_{k+1}=X_k-{(J{\left(X\right)}^{T}J\left(X\right)+\mu I)}^{-1}J{\left(X\right)}^T\mathcal{F}\left(X\right)$$

(9)

where $J\left(X\right)$ is the Jacobian matrix of $E\left(X\right)$ with respect to X. When $\mu$ is small, the step is equal to the step of Newton’s method, and when $\mu$ is large, the step is equal to the step of gradient descent.

Using the correlation coefficient as the weight, a weighted averaging operation is used to determine the optimal window size:

$$S_{opt}=\frac{\sum \left[S{\left(n\right)}_{opt}·R({\varphi }_{S\left(n\right)},{\stackrel{⌢}{\varphi }}_{S\left(n\right)})\right.}{\sum R({\varphi }_{S\left(n\right)},{\stackrel{⌢}{\varphi }}_{S\left(n\right)})}$$

(10)

3.2.2. The MMVM-Based Tropospheric Delay Correction Method

Based on the optimal window size, the interferograms with short spatial–temporal baselines are block-partitioned with an overlap, as well as the visibility map, coherence map, and the annual deformation rate derived from Stacking-InSAR. After masking the phases in deformed, non-visible, and low-quality areas, the MMVM is used to establish the relationship between the residual phase and the local terrain for estimating tropospheric delay phases. The deformed signal affects the correction modeling in a localized region. In addition, the weights derived from the deformation rate can also be used for the atmospheric delay correction in non-deformed regions. The detailed workflow of the MMVM is shown in Figure 4.

The corrected phase is obtained by weighted subtraction of the estimated tropospheric delay phase from the original unwrapped phase. For the ith point of the jth unwrapped image, the MMVM of the kth iteration can be expressed as

$$w_{i,k}^j=\left\{\begin{array}{cc}1\hfill & if:abs\left(v_{i,k}^j\right)\le v_{thr}\&i\in reliableareas\hfill \\ abs{\left[abs\left(v_{i,k}^j\right)-v_{thr}+1\right]}^{\gamma }\hfill & if:abs\left(v_{i,k}^j\right)>v_{thr}\&i\in reliableareas\hfill \\ 0\hfill & else\hfill \end{array}\right.$$

(11)

where $\gamma$ is the weight, $abs(\ast )$ represents the absolute value, and $v_{thr}$ is the velocity threshold used to distinguish between the deformed and non-deformed areas.

$$\begin{array}{cc}\hfill {\phi }_{i,k+1}^j& ={\phi }_{i,k}^j-w_{i,k}^j·\stackrel{⌢}{\phi }\hfill \\ & ={\phi }_{i,k}^j-w_{i,k}^j·\underset{\stackrel{⌢}{\phi }}{\mathrm{argmax}}[R({\stackrel{⌢}{\phi }}_{i,k}^j,{\phi }_{i,k}^j)]\hfill \\ & ={\phi }_{i,k}^j-w_{i,k}^j·\underset{\stackrel{⌢}{\phi }}{\mathrm{argmax}}(h_i^j\otimes X_{i,k}^j)\hfill \\ & ={\phi }_{i,k}^j-w_{i,k}^j·\underset{\stackrel{⌢}{\phi }}{\mathrm{argmax}}\{h_i^j\otimes \underset{X}{\mathrm{argmax}}[\mathcal{F}(X)]\}\hfill \end{array}$$

(12)

It should be noted that the estimated tropospheric delay phase of each overlapping subwindow should be interpolated to cover the entire image. And the adaptive weighting contributes to preserving the signals from the deformed areas. In each sliding subwindow, the value corresponding to the maximum of the correlation function between the unwrapped phase and the estimated tropospheric delay phase is the estimated tropospheric delay phase. The above function can be further transformed into a set of variogram models with unknown parameters X to describe the relationship between the unwrapped phase and the local topography. Multivariate variogram models are used to accurately estimate the tropospheric delay phase, and the model with the highest correlation is selected as the relationship function. The optimal parameters required by the relationship function are calculated using the LM algorithm.

Before establishing the MMVM for the next iteration, it is necessary to recalculate the deformation rate derived from Stacking-InSAR and update the corresponding weights. If the standard deviation (SD) of the phase within the stable areas is greater than the previous one, this image is skipped. The above steps are repeated until the termination condition is reached. After reducing small-scale atmospheric turbulence using temporal high-pass and spatial low-pass filters, the refined high-precision surface deformation information of the study area can be obtained.

3.3. Hotspot Analysis and Spatial Clustering

Hotspot analysis, which describes the clustering degree of MPs within a search distance, can be used to automatically identify areas of anomalous deformation, including potential landslides. The spatial search distance is the key parameter for hotspot analysis. By analyzing the incremental spatial autocorrelation, the optimal search distance can be determined. The Getis–Ord $G_i$ statistic [35] is used to evaluate the spatial variability in each monitoring point. For the ith point, the $G_i^{\ast }\left(d\right)$ within a distance d can be represented as

$$G_i^{\ast }\left(d\right)=\frac{\sum v_i+v_i-n_{ij}\times \overline{v}}{\sqrt{S^{\ast }\left\{\left[\left(n\times n_{ij}\right)-{n_{ij}}^2\right]/\phantom{\left[\left(n\times n_{ij}\right)-{n_{ij}}^2\right]\left(n-1\right)}\left(n-1\right)\right\}}}$$

(13)

where j is the jth surrounding point, $n_{ij}$ is the number of MPs within the distance d, v is the deformation rate, $\overline{v}$ is the average of the deformation rate, S is the SD of the deformation rate, and n is the total number of MPs. Based on the Z-score threshold and p-value threshold, the points with high aggregation can be extracted as hotspots [36]. However, there are still a few discrete points. Based on the point clustering method, the clustered hotspots within a certain distance can be delineated as polygons, i.e., potential geologic hazard zones.

4. Results

Following the workflow, based on the annual deformation rate derived from SBAS-InSAR, the HSA and spatial clustering methods were used to automatically identify potential geological disaster risk zones over the reservoir area of Baihetan hydropower station. There are several keys that should be noted for the correction model. A threshold of 0 for $R_{index}$ was set to exclude the signals from non-visible areas [5]. A threshold of 0.3 for coherence was used to identify coherent pixels in each coherence map. Pixels with values above the threshold were considered as high-coherence areas containing useful phase information [51]. A deformation rate threshold of 5 mm/year for Stacking-InSAR was set to extract signals from non-deformed areas [27]. Given the advantages of focusing on the local features of the data while minimizing global errors, in this paper, cubic spline interpolation was used to restore the unknown signals [52]. The trend of the overall corrected unwrapped phase was eliminated by a least-squares fit of a bilinear equation. For hotspot analysis, a Z-score threshold of ±1.96 and a p-value threshold of 95% were used to identify hotspots [35]. Hotspots were automatically delineated into at-risk zones using a distance of 200 m [4].

The original identification results of at-risk zones and those after correction by the MMVM are illustrated in Figure 5. The accuracy of the identified at-risk zones is determined through the analysis of the Stacking-InSAR deformation rate, SBAS-InSAR deformation rate, cumulative deformation, optical imagery, visual interpretation, and factors such as the undulations of the local terrain, distance from the river, and whether it is historically a site of geological hazards. First of all, based on the deformation rate derived from SBAS-InSAR, hotspot analysis and the point aggregation method are used to automatically identify the risk zones in the study area. The SBAS-InSAR deformation rates of the identified risk zones are then superimposed on the optical image and compared with the maximum cumulative deformation and the Stacking-InSAR deformation rate to comprehensively determine whether they are potential risk zones by visual interpretation. Based on their geographic location in the optical image, these potential geologic hazard zones can be further classified into different types of geologic disasters, such as landslides, urban subsidence, ground cracks, fissures, and others. Finally, the identified at-risk areas are further verified, if possible, by field investigation. The accuracy statistics of the identification results before and after the MMVM correction are shown in

Table 1. Comparison of the identification results before and after MMVM correction.

Method	At-Risk Points	At-Risk Zones	Accuracy
Original	3301	176	47.159%
After Correction	2297	133	96.241%

.

By applying the HSA to the original deformation rate, a total of 3301 MPs were identified as at-risk points, and 176 at-risk zones were formed. After correction by the MMVM method, a total of 2297 MPs were identified as at-risk points, representing a reduction of 30.415%, forming 133 at-risk zones, representing a reduction of 24.432%. The accuracy rate of at-risk zones identified by hotspot analysis increased from 47.159% to 96.241%. A total of 32 active landslides were identified, 6 of which were along the Jinsha River. And, 16 identified active landslides pose a threat to nearby villages, with detailed information summarized in

Table 2. Detailed information about detected active landslides threatening nearby villages.

No.	Name	Longitude	Latitude	Area	Maximum Rate	Maximum Deformation	Aspect	Threat Object
		(°)	(°)	(km2)	(mm/year)	(mm)
1	Luotianba	103.58	27.98	0.031	−24.098	−36.274	W	River, villages
2	Sunjialiangzi	103.59	27.94	0.017	−28.744	−40.041	SW	Villages, roads
3	Lijiaping	103.52	27.93	0.168	−25.042	−57.358	NW	Villages, river, and roads
4	Wozitou	103.52	27.71	0.039	−57.750	−113.923	W	Villages, farmlands, and roads
5	Qianligou	103.32	27.69	0.022	−24.243	−34.381	S	Villages, farmlands, and roads
6	Gongshan	103.24	27.59	0.055	−27.401	−48.558	SW	Villages, farmlands
7	Liangshanjing	103.23	27.48	0.025	−41.725	−75.183	SW	Villages, farmlands
8	Yujiapingzi	103.25	27.46	0.022	−25.286	−66.031	W	Villages, farmlands, and roads
9	Dayadong	103.25	27.40	0.032	−26.053	−38.562	S	Villages, roads
10	Galuo	102.95	27.45	0.336	−31.396	−54.545	W	Villages, roads
11	Shanshu	103.25	27.38	0.001	−29.191	−58.825	S	Villages, roads
12	Niupingyan	103.12	27.39	0.034	−35.932	−72.767	SW	Villages, roads
13	Ertun	103.00	27.40	0.098	−32.063	−72.047	SE	Villages, roads
14	Youyicun	102.85	27.14	0.130	−34.710	−51.933	SE	Villages, farmlands
15	Huodi	103.06	26.90	0.570	−30.339	−52.472	SW	Villages, farmlands
16	Xintian	103.12	26.67	0.015	−22.794	−48.206	E	River, villages

. The largest landslide covers an area of 0.57 km², while the smallest covers an area of 8642 m². The Wozitou landslide exhibits the maximum deformation rate of −57.75 mm/year, which corresponds to a cumulative deformation of −113.923 mm and covers an area of 0.039 km². This landslide required additional attention as it threatened villages, farmland, and roads.

Overall, the original results exhibit sparser MPs and numerous outliers compared to the MMVM-corrected results, particularly in the southern part of the study area. The high-relief regions exhibit spatial variability in atmospheric properties, which affects the retrieval of deformation information. After the atmospheric correction, both the quantity and quality of MPs have been improved to some extent. Both successfully identified the subsidence zones of area S1, providing evidence of the identification capability of the HSA.

We select two areas, area S1 and area S2, for further analysis of the identification results from HSA. The original identification results of the at-risk zones and those after correction by the MMVM method, covering areas S1 and S2, are illustrated in Figure 6. It should be noted that the boundaries of the at-risk zones vary due to differences in the identified hotspots, which are extracted based on the statistical significance of the spatial clustering for each MP. Figure 6a,b demonstrate an improvement in both the quantity and quality of MPs after tropospheric delay correction. The points in stable areas have more than doubled, with a lower SD of 2.756 mm/year compared to the original value of 2.813 mm/year. Area S3 is characterized by severe subsidence with sparse MPs. The deformed areas, Mahao and Hongshan, become apparent after correction by the MMVM method. In area S4 of Figure 6b, MPs cover almost the whole area, and after tropospheric delay correction, six new at-risk zones have been identified, resulting in a reduction of 14 at-risk zones from the original identification results. Comparing Figure 6c,d, the number of MPs increased from 6147 to 11,704, with the SD being 2.818 mm/year in contrast to the initial 2.849 mm/year. The performance is particularly striking in the area along the western side of the Jinsha River. There were initially 22 at-risk zones, and after correction using the MMVM method, the number of at-risk zones increased to 30, with 20 common zones. Especially in the Gantianba area, the increase in identified hotspots contributed to the expansion of the at-risk zone. In area S5 of Figure 6d, the ground subsidence of the Siguxi area is noticeable, as the overall atmospheric noise is suppressed, leading to an improvement in the quality of the deformation results. Area S6 shows a situation similar to that of area S5, where several common hotspots have been identified, such as Miaozhai.

5. Discussion

5.1. Performance Evaluation of MMVM-Based Correction Method

Numerous tropospheric delay correction methods have been proposed by establishing an exponential model over the entire area to remove the tropospheric delay component associated with the terrain, known as the Exp method [32]. Based on the meteorological data from ECMWF, the ERA5 model is employed to mitigate the atmospheric phase of interferograms [17]. The GACOS [19] is an online data service based on GNSS and ECMWF data, which facilitates the InSAR atmospheric correction. In addition, the traditional SBAS method follows the standard processing flow [5]. The above methods perform well in many cases, and are used to assess the effectiveness and performance of the proposed atmospheric correction method. Since the proposed method involves iterative optimization, we present results for 5 iterations, denoted as MMVM(5), and 10 iterations, denoted as MMVM(10). Referring to the metrics of [27,53], the performance comparison of the different correction methods is conducted based on three aspects: statistical evaluation of corrected interferograms, improvement in derived deformation rate, and stability of time-series deformation.

5.1.1. Statistical Evaluation of Corrected Interferograms

Without taking into account the phase errors such as height error, orbit error, etc., that exist in the unwrapped phases, both the mean and the SD are used to evaluate the correction effect. The interferometric pair of 1–13 November 2021 is chosen as an illustrative example to show the correction effects of different tropospheric delay methods, as shown in Figure 7. The mean and the SD within the stable areas are shown in the lower-right corner of the subplot. To provide an intuitive comparison between the unwrapped phase and the elevation, the terrain of the study area is shown in Figure 7h.

After the correction by ERA5 and GACOS, the unwrapped phase still shows terrain-related atmospheric components when compared to the corresponding terrain. By using an exponential model to establish the relationship between the atmospheric delay phase and the local terrain, Exp yields fewer terrain-related atmospheric components compared to ERA5 and GACOS. The SBAS method shows a remarkable correction performance and its SD of 1.268 rad is lower than those of ERA5 of 2.234 rad and GACOS of 1.785 rad. The SBAS method incorporates a linear regression function to estimate the terrain-related atmospheric components and a polynomial least-squares model to eliminate the residual phase slope. Compared to MMVM(5), MMVM(10) shows a phenomenon of overcorrection in some areas, with a corresponding decrease in the SD value. It is therefore recommended to set a reasonable threshold to avoid overcorrection. The advantage of MMVM lies in its ability to eliminate atmospheric components while preserving deformation information, as can be seen in the red anomaly area in the center of Figure 7f.

The mean and SD values of unwrapped phases before and after correction by different methods within stable areas are shown in Figure 8. The statistical characteristics of the corrected phases within stable areas were carried out for six tropospheric delay correction methods, as listed in

Table 3. Statistical evaluation of original unwrapped phases and corrected unwrapped phases.

Method	Number with Reduced SD	SD (rad)	Number with Mean Closer to 0	Mean (rad)
Original	110	2.3228	110	−0.0327
Exp	110	2.0449	106	0.0101
ERA5	107	1.8854	61	−0.0606
GACOS	43	2.1990	42	0.0326
SBAS	39	2.2711	36	0.0103
MMVM(5)	110	1.1893	106	0.0040
MMVM(10)	110	1.0123	107	0.0045

.

The correction effect of terrain-related tropospheric delay varies among the six methods. The original unwrapped phase has a mean SD of 2.3228 rad, while the alternative methods give the following values: 2.0449, 1.8854, 2.199, 2.271, 1.1893, and 1.0123 rad. Among the first four correction methods, ERA5 has the maximum SD reduction, reaching 1.8854 rad, and the corresponding number of images is 107. Both the Exp and MMVM methods have a constant number of images in the SD reduction, at 100%. However, the MMVM outperforms the others in both the quantity and magnitude of the SD reduction. Given that the closer the mean is to 0, the better the correction effect will be, the Exp method performs well among the first four correction methods, and the number of corresponding maps is 107, consistent with the MMVM(10) method. In summary, the MMVM method has the minimum SD value and the mean value closest to 0 of all the methods, with corresponding correction ratios of 100% and 97.27%, respectively.

5.1.2. Improvement in Derived Deformation Rate

The identification of geological hazards, such as landslides, depends on the deformation rate derived from InSAR at MPs. The quality of the derived deformation can be assessed by analyzing the statistical characteristics of the MPs within both non-deformed and deformed areas of the study area. In this study, two types of deformation rate were obtained using SBAS-InSAR and Stacking-InSAR. The former uses the SVD to extract the atmospheric and nonlinear deformation information, followed by a spatial–temporal filter, resulting in reliable results. The deformation rate derived from the Stacking-InSAR assumes stationary atmospheric statistics for the set of interferograms. Given that potential landslides often occurred on both sides of the Jinsha River, we established a 10 km buffer zone over the reservoir area of the Baihetan hydropower station to extract the corresponding deformation information. The coverage of the reservoir area is shown in Figure 9a. The original deformation rate derived from SBAS-InSAR and the deformation rate after correction by various methods are shown in Figure 9b–g. The distributions of the LOS deformation rate for MPs are shown in Figure 10. The statistical characteristics of the SBAS-InSAR deformation rate using various correction methods are shown in

Table 4. Statistical characteristics of SBAS-InSAR deformation rate using various correction methods.

Method	MPs	Non-Deformed Area		Deformed Area
		MPs	SD	Mean	SD
			(mm/year)	(mm/year)	(mm/year)
Original	173,349	61,530	2.8320	−2.9009	13.3601
Exp	218,352	81,333	2.8310	−1.8157	13.4583
ERA5	170,706	62,201	2.8274	−2.7907	13.3704
GACOS	152,046	56,331	2.8346	−2.5942	13.3697
SBAS	245,442	95,205	2.8307	−2.0812	12.8820
MMVM	296,813	150,071	2.7633	−1.3111	11.0191

.

Compared to the original results, the MMVM provides a higher density of MPs, as seen in the enlarged views of area S1. Severe subsidence is observed in area S1, attributed to open-pit mining, with a maximum deformation rate of −75.495 mm/year. The original results have a total of 173,349 MPs, with 61,530 MPs (35.495%) within non-deformed areas. The GACOS method shows poor performance in atmospheric phase correction. The SBAS method gives favorable results, which can be attributed to the removal of the residual phase tendency. Compared to the original results, the Exp method yields more MPs, increasing by 25.961%, with the number of MPs within stable areas reaching 81,333. For SBAS, the number of MPs within stable areas increases to 95,205. However, with MMVM this increases to 150,071, a growth of 143.899%, with the low SD of 2.7633 mm/year. Compared to the other correction methods in the deformed areas, the mean of the MMVM is closer to 0 and the SD is the minimum at 11.0191 mm/year. The SBAS method has the second-lowest SD among all the methods. For the statistical characteristics of all MPs, the mean values are −1.93, −1.18, −1.83, −1.68, −1.32, and −0.74 mm/year, while the SD values are 10.94, 10.83, 10.87, 10.81, 10.28, and 8.01 mm/year. The ERA5 and GACOS methods exhibit poor correction results with SDs of 10.87 mm/year and 10.81 mm/year due to significant atmospheric delay variations during image acquisition. In contrast, the SBAS method shows favorable correction results with an SD of 10.28 mm/year. The MMVM method stands out with a mean of −0.74 mm/year and SD of 8.01 mm/year, as shown in Figure 10f.

From the enlarged views of area S1 in Figure 9, the deformation rates derived by the ERA5 and GACOS methods have a higher degree of displacement deviation. The derived deformation of the majority of the surface monitoring points should tend toward stability, which is shown in green on the deformation rate map. There are numerous erroneous deformation signals in non-deformed areas, especially the positive deformation signals in canyon areas. Due to the severe terrain undulations and the presence of the atmospheric delay, the phase unwrapping fails and the erroneous deformation signals appear. The effectiveness of different correction methods varies. The GACOS products introduce uncertainties to a certain extent, particularly in the Qiaoia area of area S1, which makes MPs sparse and does not provide good coverage of the deformed area. However, compared to the original result, there is some improvement. The Exp method, through the establishment of an exponential model, gives significantly better results than the ERA5 and GACOS methods. The deformation results of the method proposed in this paper are purer, with a higher density and more pronounced deformation area.

The statistical characteristics of Stacking-InSAR deformation rate for the entire area and area S1 are listed in

Table 5. Statistical characteristics of Stacking-InSAR deformation rate within the entire study area.

Method	Non-Deformed Areas		Deformed Areas
	MPs	SD	Mean	SD
		(mm/year)	(mm/year)	(mm/year)
Original	4,329,987	0.6469	−0.0020	2.5300
Exp	5,601,942	0.6464	−0.0248	3.3740
ERA5	4,772,330	0.6444	0.0006	2.3842
GACOS	4,657,488	0.6478	0.0054	2.8266
SBAS	4,644,076	0.6464	−0.0017	2.5607
MMVM	8,897,582	0.6236	−0.0482	2.4970

and

Table 6. Statistical characteristics of Stacking-InSAR deformation rate within area S1.

Method	Non-Deformed Areas		Deformed Areas
	MPs	SD	Mean	SD
		(mm/year)	(mm/year)	(mm/year)
Original	98,255	0.6431	0.0008	2.2065
Exp	107,778	0.6426	−0.0004	2.7814
ERA5	118,150	0.6374	0.0007	2.0490
GACOS	86,101	0.6347	0.0011	2.4283
SBAS	97,495	0.6443	0.0008	2.2063
MMVM	179,376	0.5882	−0.0005	1.7713

, respectively. In Table 5, the number of stable MPs for the MMVM more than doubles compared to the original result, accompanied by a smaller SD. This is attributed to the iterative weighting operation, which protects the deformed zone and stabilizes the non-deformed zone. In terms of the number of non-deformed MPs, the Exp method achieves a larger number of MPs, an increase of 29.375%. However, the GACOS method exhibits the smallest increase, with 10.216%. Note that negative values mean that the surface moves away from the satellite along the line of sight (LOS). The Exp method has a negative mean and a small SD of 2.3842 mm/year. Compared to the other methods, the MMVM has the smallest mean and SD of −0.0482 mm/year and 2.4970 mm/year, respectively. The superiority of the MMVM method can also be seen in Table 6. Within the non-deformed areas, the MMVM method has the maximum number of MPs and the lowest SD value of 0.5882 mm/year. The ERA5 method performs well compared to the other four methods, with a large number of MPs and a low SD value within area S1.

5.1.3. Stability of Time-Series Deformation

Atmospheric distortion can affect the deformation trend in the study area. Assuming no residual atmospheric phase and decorrelation noise, for stable MPs, the mean of the time-series deformation and the standard deviation from fit (SDF) are theoretically close to 0 [27]. SDF can be used for the assessment of the stability of the evolutionary trend in the time-series deformation for MPs. Two feature points are selected to evaluate the correction effect: the maximum subsidence point (FP1) in area S1 and the stable point (FP2) in area S2. The parameters using various correction methods at FP1 and FP2 are listed in

Table 7. Comparison of parameters using various correction methods at FP1 and FP2.

Method	Standard Deviation from Fit
	FP1 (mm)	FP2 (mm)
Original	6.278	4.213
Exp	5.965	4.267
ERA5	6.259	4.321
GACOS	6.303	4.863
SBAS	6.684	5.306
MMVM	4.773	3.022

.

FP1 corresponds to the point with the maximum subsidence in area S1, located at 103.209°E, 27.462°N, with an annual subsidence rate of −75.495 mm/year. According to optical imagery and ground surveys, the open-pit mining of mineral resources has resulted in severe subsidence, as shown in Figure 11a. A cubic polynomial method was used to fit the time-series deformations, with a minimum SDF of 4.773 mm for the MMVM method. Figure 11c shows the time-series cumulative deformations and the corresponding fit curve at FP1 during the observation period. FP2 represents the stable point in the Gantianba at-risk zone of area S1, located at 103.543°E, 27.942°N, with an annual subsidence rate of 1.427 mm/year. The location of FP2 in the optical imagery is shown in Figure 11b. The SDF of the MMLM at FP2 is the lowest of all the methods at 3.022 mm. The results demonstrate the superiority of the MMVM over other methods and highlight the importance of tropospheric delay correction.

5.2. Advantages and Limitations of the Proposed Method

Based on the deformation rate derived from SBAS-InSAR the Getis–Ord $G_i$ statistic is applied to evaluate the spatial variability in each MP within the study area. Then, using the Z-score and p-value thresholds for hotspot analysis, all the feature points with high aggregation, i.e., the deformation points, can be identified. The point aggregation method was used to automatically delineate these hotspots into risk zones with a certain distance threshold. However, the spatially heterogeneous tropospheric delay biases the InSAR-derived surface deformation, especially in regions of high relief. Considering topographic factors and surface deformation, the MMVM method is introduced to correct the stratified tropospheric delay for generating high-accuracy displacement maps. We then developed an automatic identification method to delineate wide-area geological hazard risk areas, including potential landslides. The proposed workflow can be extended to other experimental scenarios, such as risk assessment of geological hazards like the catastrophic landslides in the Himalayas. The proposed MMVM-based tropospheric delay correction method is free from the limitations of the traditional methods [25], which require the presetting of the segmentation window size and the priori relationship model, and greatly improves the computational efficiency. Compared to the method of Liang et al. [24], which is based on elevation gradient for window segmentation, the MMVM-based method adopts an iterative weight updating strategy to overcome the effect of the heterogeneous distribution of atmospheric turbulence on the establishment of relationship models. Unlike the GACOS and ERA5 methods, the MMVM-based method does not need to take into account meteorological conditions and avoids uncertainties in the estimated atmospheric delay phase introduced by phase errors due to the low spatial resolution of meteorological data. Furthermore, compared to a single empirical model like the Exp method [30], the proposed method based on the multivariate variogram models is able to account for the spatial heterogeneity in the atmospheric delay, thus improving the adaptability and robustness of the algorithm. A comparison of the performance of the different correction methods is shown in

Table 8. Comparison of the performance of different correction methods.

Method	Plain	Steep Terrain	Computational Efficiency	External Data	Advantages	Disadvantages
Exp	Good	Good	Medium	No	Satisfies the spatial heterogeneity	Sensitive to deformed or turbulent signals, small-scale
ERA5	Good	Poor	Low	Yes	Wide-scale	Uncertainty in estimated tropospheric delay phases
GACOS	Good	Poor	Medium	Yes	Wide-scale	Uncertainty in estimated tropospheric delay phases
SBAS	Good	Good	High	No	Various monitoring scenarios	Unable to satisfy the spatial heterogeneity
MMVM	Excellent	Excellent	Low	No	Wide-scale, various monitoring scenarios, and satisfies the spatial heterogeneity	Phase overcorrection in some areas

. As a result, the MMVM-based method partially produces MPs of a higher quantity and quality for the refinement of the InSAR deformation map.

However, the proposed method still has some limitations. The weight derived from Stacking-InSAR deformation information affects the correction performance of the MMVM model. When the weights are small, it maintains the actual deformation information of the MPs, but increases the running time. Therefore, the appropriateness of the formula used to derive the weights needs to be further investigated. The low coherence of the pixels in steep mountainous areas leads to unwrapped phase errors. As a result, the delay phase estimated by the correction model is inaccurate. To solve this problem, a robust phase unwrapping method is required. In addition, hotspot analysis for automatic landslide identification has certain limitations, especially in mountainous areas. Hotspot analysis extracts the risk points based only on the degree of clustering of the deformation rate for each MP. The annual deformation rate provides an important indication of where geological hazards occur. On the one hand, the deformation rate ignores some areas of nonlinear deformation. On the other hand, hotspot analysis does not take into account other factors such as lithology, dip, etc., although it has been successfully applied and validated in many cases. Therefore, the method can only delineate the areas where potential geological hazards may occur and further investigation by ground investigation is required.

6. Conclusions

Stratified tropospheric delay can introduce significant biases in the deformation derived from InSAR, significantly affecting the identification of slow-moving landslides, especially in regions of high relief. Therefore, it is crucial to study the tropospheric delay correction.

In this study, we propose an improved tropospheric delay correction method based on the MMVM, and combined with hotspot analysis for the automatic identification of wide-area landslides. Considering topographic factors and surface deformation, the MMVM uses multivariate variogram models with iterative weight derived from Stacking-InSAR to adaptively determine the window size and the model based on the distribution characteristics between the unwrapped phase and the local terrain, and uses the LM algorithm to improve both the convergence speed and robustness. The lower reaches of the Jinsha River above the reservoir area of the Baihetan hydropower station were selected as the study area. Using a total of 118 Sentinel-1A images, an exponential model, ERA5, GACOS, and the traditional SBAS methods were applied to demonstrate the performance of the proposed method in terms of corrected interferograms, derived deformation rate, and stability of time-series deformation. In the corrected interferograms, the MMVM method has the minimum SD and the mean closest to 0 among all the methods. Regarding the derived deformation, we made a comparison of the statistical characteristics of MPs, and the MMVM-corrected results provide a higher quantity and quality of MPs in the study area. To compare the stability of the evolutionary trend in time-series deformation for MPs, two feature points are selected: the maximum subsidence point (FP1) in area S1 and a stable point (FP2) in area S2. Both of them show the minimum SDF of 4.773 mm and 3.022 mm, compared to the original values of 6.278 mm and 4.213 mm, respectively. In addition, the accuracy of the at-risk zones identified by the HSA has improved from an initial of 47.159% to 96.241% after the correction by MMVM. Finally, a total of 32 active landslides were identified, 6 of which occurred along the Jinsha River. The proposed method can provide support for the geological hazard investigation department.