Slow-Moving Landslide Hazard Assessment Using LS-Unilab Deep Learning Model with Highlighted InSAR Deformation Signal

Abstract

Slow-moving landslides are often precursors of catastrophic failure, posing a major threat to human life and property safety. Interferometric synthetic aperture radar (InSAR) has become a crucial tool for investigating slow-moving landslides hazard because of its high-precision detection capability for slow surface deformation. However, landslides usually occur in alpine canyon areas and vegetation coverage areas where InSAR measurements are still limited by temporal and spatial decorrelation and atmospheric influences. In addition, there are several difficulties in monitoring the multiscale characterization of landslides from the InSAR results. To address this issue, this paper proposes a novel method for slow-moving landslide hazard assessment in low-coherence regions. A window-based atmosphere correction method is designed to highlight the surface deformation signals of InSAR results in low-coherence regions and reduce false alarms in landslide hazard assessment. Then, the deformation annual velocity rate map, coherence map and DEM are used to construct the InSAR sample set. A landslide hazard assessment model named Landslide-SE-Unilab is subsequently proposed. The global–local relationship aggregation structure is designed to capture the spatial relationship between local pixel-level deformation features and global landslides, which can reduce the number of missed assessments and false assessments of small-scale landslides. Additionally, a squeeze-and-excitation network is embedded to adjust the weight relationship between the features of each channel in order to enhance the performance of network evaluation. The method was evaluated in Kangding city and the Jinsha River Valley in the Hengduan Mountains, where a total of 778 potential landslides with slow deformation were identified. The effectiveness and accuracy of this approach for low-coherence landslide hazard assessment are demonstrated through comparisons with optical images and previous research findings, as well as evaluations via time-series deformation results.

1. Introduction

Although slow moving landslides rarely cause catastrophic losses, they can accelerate rapidly, and even catastrophically, due to a variety of internal and external factors [1,2,3,4], causing widespread destruction and casualties. With global changes in climate, the hazard of slow-moving landslides leading to catastrophic collapse is increasing. On the basis of the back analysis of various landslide cases, slow ground motions have often been identified as precursory deformation signals before failure. Accordingly, slow-landslide monitoring is still regarded as an effective way to evaluate hazards, even for early-warning.

Interferometric synthetic aperture radar (InSAR) is a technology that processes the phases of SAR images to obtain surface elevation and deformation information. With the advantages of all-weather, all-time, high resolution and large-scale surface deformation monitoring, InSAR is widely used in geological hazard monitoring. Multitemporal InSAR (MT-InSAR) technology analyzes surface time-series deformation by using high-coherence targets in time-series SAR images to overcome the limitations of traditional DInSAR. MT-InSAR techniques include a series of methods, such as stacking-InSAR [5], PS-InSAR [6], SBAS-InSAR [7], StaMPS [8], TCPInSAR [9,10], and Tomo-PSInSAR [11]. The time-series deformation information provided by MT-InSAR can be used to study the mechanisms and kinematics of landslides, analyze landslide influencing factors, and conduct susceptibility and hazard assessments. Scholars have conducted a great deal of research on landslides via MT-InSAR technology, including detecting and identifying landslides [12,13,14,15], classifying landslide processes and styles [16,17,18], capturing the precursory signals of slope rupture [19], inferring subsurface geometry [20], modeling kinematic behaviors [21], and assessing hazard and risk [22]. Owing to the slow-motion rate of slow-moving landslides, which can last from months to decades or even centuries, with movement rates ranging from a few millimeters to a few meters per year [2,23], the deformation signal obtained via InSAR techniques can be masked by the various noises caused by temporal and spatial decoherence and atmospheric effects, especially in the low coherence area, such as mountainous areas with trees where landslides often occur. Empirical models, which are commonly used to calculate atmospheric errors, encounter challenges in high mountains because of the spatial variability of tropospheric properties [24].

Extraction of the deformation area and evaluation of landslide likelihood are two primary tasks for landslide hazard assessment. For deformation area extraction, techniques such as visual interpretation, threshold segmentation and deep learning are often applied. Visual interpretation relies on personal experience and is less applicable to large areas. Threshold segmentation can lead to the following issues: applying the same threshold to areas with different geological structures, lithologies, and vegetation cover can introduce errors and inaccuracies in deformation area extraction; the determination of the threshold value may depend on subjective judgment by different researchers, while some scholars have made efforts to automate its selection on the basis of statistical characteristics of deformation and have achieved certain advancements [25,26]. Deep learning methods depend on the quality of the training datasets and, in previous studies, deformation results obtained via InSAR have rarely been used as the dominant feature in datasets [27]. The multiscale feature of the InSAR sample set also makes its assessment by existing models unsatisfactory. In terms of landslide-prone zone assessment, environment factors such as slope, aspect, distance to rivers, and geological information are applied to evaluate landslide susceptibility [28]. Susceptibility can only describe the regional probability of landslide occurrence and cannot be used to evaluate individual landslides. These indicators are often regional and cannot be used to determine the hazard level of a single landslide.

On the basis of the issues described and scholars’ research, an automatic slow-moving landslide hazard assessment method applicable to InSAR results is proposed in this paper for low-coherence areas, such as alpine canyon areas and vegetation coverage areas. First, adaptive windows are used to effectively suppress atmospheric signals in InSAR results, thereby highlighting surface deformation feature signals. Then, the deformation annual velocity rate map, coherence map and DEM are superimposed, and the slow-moving landslides are labeled by combining optical images, field investigations and historical landslide inventory data to produce the InSAR deformation sample set. Finally, the Landslide-SE-Unilab (LS-Unilab) semantic segmentation model based on the Uniformer backbone network and the squeeze-and-excitation network is constructed for automated landslide hazard assessment.

2. Methods

The methodology adopted in this study can be summarized in three key steps, as depicted in Figure 1. (1) The SBAS-InSAR deformation features are highlighted by using a window-based atmospheric correction. (2) The deformation annual velocity rate map, coherence map and DEM are superimposed, and slow-moving landslides are labeled with multisource data to produce an InSAR deformation sample set. (3) The LS-Unilab semantic segmentation model is employed to train and perform quantitative analysis on the InSAR sample set to assess landslide hazards.

2.1. SBAS-InSAR Surface Deformation Velocity Inversion

Unlike persistent scatterer interferometry (PSI) techniques, which are more suitable for the monitoring of point scatterers, such as man-made structures and, in general, the built environment, the SBAS technique is well suited to the monitoring of distributed scatterers, i.e., ground textures characterized by a spatially uniform backscattering within the resolution cell, which is typically associated with the natural environment where landslides often occur [29]. However, the deformation results are often contaminated by complex atmospheric conditions in mountainous areas. To highlight the deformation signals, a window-based atmospheric correction method is proposed to minimize the influence of the atmosphere.

The specific procedure is illustrated in Figure 2. First, the optimal interferometric pairs are selected on the basis of temporal and spatial baselines to produce the aligned SLC image stack [30], in which multi-looks of 2 and 8 are applied in the azimuth and range directions, respectively. Next, we calculate the mean spatial coherence of each interferogram and automatically remove those interferograms below the coherence threshold to reduce errors in the inversion of the time-series deformation. Subsequently, the SNAPHU algorithm is employed to produce a stack of unwrapped interferograms [31]. High-coherence targets are then selected, and initial time-series deformation is inverted via either the least squares (LS) method or singular value decomposition (SVD). Finally, a window-based model for atmospheric correction is applied to derive the time-series deformation and residual topographic height error.

In the initial phase estimation stage of InSAR processing, the estimated original phase can be described as:

$$\widehat{\varphi }={\widehat{\varphi }}_{dis}+{\widehat{\varphi }}_{tropo}+{\widehat{\varphi }}_{geom}+{\widehat{\varphi }}_{resid}$$

(1)

where ${\widehat{\varphi }}_{dis}$ is the surface deformation phase, ${\widehat{\varphi }}_{tropo}$ is the atmospheric phase, ${\widehat{\varphi }}_{geom}$ is the topographic phase, and ${\widehat{\varphi }}_{resid}$ is the residual phase. The topographic phase can be generated from the DEM. Finally, the surface displacement is estimated as:

$${dis}_i^{\phi }={\phi }^i-{tropo}_i^{\phi }-{resid}_i^{\phi }$$

(2)

$${dis}_i^{\phi }=\nu \cdot t_i{\displaystyle \frac{-4\pi }{\lambda }}+\xi +{tropo}_i^{\phi }$$

(3)

where ${dis}_i^{\phi }$ is the displacement phase, ${\phi }^i$ is the original phase, ${tropo}_i^{\phi }$ is the atmospheric phase, $\xi$ is the residual phase, and $\nu$ is the linear surface displacement rate obtained via the least squares method. When calculating the annual average deformation rate, a time-weighted approach is typically used:

$$v_{ avg}=f\left(v_i\right)={\displaystyle \frac{\sum _{i=1}^N \left(v_i\cdot T_i\right)}{\sum _{i=1}^N T_i}}\cdot {\displaystyle \frac{365}{T_i}}$$

(4)

where $v_i$ is the surface displacement rate of the interferometric pair, and $T_i$ is the time interval. Then, we remove the errors caused by the atmospheric factor in $v_{avg}$. Phase-based techniques and auxiliary data products, such as the GACOS and ERA5 products, are commonly employed to remove atmospheric interference over an entire image [32]. However, when slow-moving landslide detection is conducted in mountainous regions with substantial differences in elevation, several challenges emerge. First, the spatial heterogeneity of the atmosphere in mountainous areas is relatively strong, and more local correction methods are needed. Second, the effectiveness of using auxiliary data to mitigate atmospheric delay is strongly influenced by the time gap between SAR acquisition and data collection. Finally inadvertently removing the deformation signal associated with slow-moving landslides that are linked to elevation is avoided [33]. Therefore, we used a window-based atmospheric correction method to eliminate atmospheric phase residues in the annual deformation rate. The relationship can be represented as:

$$\widehat{v}=v_0^i-{tropo}_i^{\epsilon }$$

(5)

$${tropo}_i^{\epsilon }=k·h+{\epsilon }_0$$

(6)

where $v_0^i$ is the deformation rate of the $i$ pixel, ${tropo}_i^{\epsilon }$ is the deformation rate error component influenced by the atmospheric factor for the $i$ pixel, $h$ is the elevation value, ${\epsilon }_0$ is the constant term of atmospheric error, representing the component of atmospheric error independent of elevation, and $k$ is the unknown parameter. In our method, the deformation rate map is divided into several windows. Equation (6) is used to estimate the atmospheric error of each window. Finally, the atmospheric errors of all windows are spliced and subtracted to get the corrected rate map.

The entire scene of the deformation rate map is divided into multiple windows (Figure 2), each of which is d × d pixels in azimuth and range directions. In addition, there is an overlapping region of columns and rows in the range and azimuth directions between the neighboring windows. Different window sizes were used for each study area. The window size is determined on the basis of the magnitude of the heterogeneous atmospheric delay and the scale of the maximum deforming area. Barnhart and Lohman’s evidence suggests that the range of terrain-dependent atmospheric delay should be less than 10-km in the alpine region [34]. Within this constraint, the window size is then determined based on the basis of the scale of maximum deformation in each study area. For each window, the atmospheric factor error is estimated with Equation (7) on the basis of the selected M pixels, which can be mathematically written as follows:

$$\left[\begin{array}{cc}{h}_1& 1\\ h_2& 1\\ ⋮& ⋮\\ h_N& 1\end{array}\right]\left[\begin{array}{c}k\\ {\epsilon }_0\end{array}\right]=\left[\begin{array}{c}{tropo}_1^{\epsilon }\\ {tropo}_2^{\epsilon }\\ ⋮\\ {tropo}_N^{\epsilon }\end{array}\right]$$

(7)

The explanation for each term is consistent with Equation (6). Least-squares is typically used to estimate the unknown parameter. Because the tropospheric delay for each window is calculated independently, the parameters of h and ${\epsilon }_0$ inevitably differ between adjacent windows due to varying atmospheric conditions. To address this misalignment, we calculate the mean delay value within overlapping regions to ensure a smoother transition between windows. After estimating the atmospheric phase contribution for each window, we generate an atmospheric artifact map for the entire scene by mosaicking the tropospheric delays from all windows. Finally, the atmospheric factor error is subtracted to obtain the refined annual deformation velocity rate map.

2.2. InSAR Sample Dataset

The InSAR dataset for landslide assessment was obtained by superimposing the deformation annual velocity rate map, coherence map and DEM. The deformation rate results reflect the process of slow landslide change; the coherence indicates the reliability of the deformation results; and the DEM refines the landslide units and provides more accurate landslide boundaries. The process and sample examples are shown in Figure 3.

When the deformation rate maps are obtained by the improved SBAS-InSAR technique, coherence is used as a key metric to assess the noise intensity. The regions below the mean value of the temporal coherence (µ) are masked and severely suffer from decoherence noise. The time-series average coherence Coh is defined by Equation (8)

$${\displaystyle \frac{1}{{Coh}^2}}={\displaystyle \frac{1}{M}}\sum _{k=1}^M \left({\displaystyle \frac{1}{{Coh}_k^2}}\right)$$

(8)

where ${Coh}_k$ is the coherence map of each interferogram pair, and $M$ is the number of interferograms.

The initial deformation regions are determined on the basis of the statistics of the displacement rate. If the absolute value of the deformation rate exceeds twice the standard deviation of the deformation rate (2-σ), the region is initially classified as a deformation region on the basis of Equation (9); otherwise, it is considered a stable region.

$$deformation region=\left\{\begin{array}{c}Coh>\mu \\ \left|V_{los}\right|>2\sigma \end{array}\right.$$

(9)

To ensure the accuracy of the labeling, the initial deformed regions are then further validated by combining optical images, a priori landslide inventories (

Table 1. Landslides already recorded within the study area in NASA’s Global Landslide Catalog (https://gpm.nasa.gov/applications/landslides (accessed on 16 May 2024)).

Latitude	Longitude	Country	Province	Landslide Size
30.0553	101.9649	China	Sichuan	large
30.0694	101.44	China	Sichuan	medium
30.0091	101.9374	China	Sichuan	medium
30.07755	102.1447	China	Sichuan	large

), and field survey data to exclude the deformation regions caused by other factors.

Subsequently, the boundaries of the landslide area are refined on the basis of the topographic features of the landslide with a digital elevation model (DEM).

Finally, the hazard level of each area of a slow-moving landslide is assessed for each landslide pixel on the basis of the rate of deformation. The safety level of the non-landslide area is designated as “safe”, while three intervals are established by adding or subtracting 3, 4, and 5 standard deviations from the average deformation rate. Disaster levels are defined as “low”, “medium”, “high”, and “very high”.

2.3. LS-Unilab Model

Landslide features in the InSAR sample set have multiscale spatial characteristics, so it is necessary to simultaneously capture the spatial relationships between large-scale landslides and single-pixel features and effectively distinguish small-scale landslides from noise. The classical semantic segmentation model Deeplabv3+ [35] uses only dilated convolution, which does not express shallow information adequately and is not capable of capturing multiscale features of landslides. In addition, the model has poor flexibility in terms of feature fusion and dynamic weight adjustment, making it difficult to adapt to multichannel feature inputs and accurately recognize complex landslide areas. In response to these problems, the LS-Unilab model was proposed, which adds a multi-head relation aggregator that integrates multiscale features and channel attention that adjusts the dependency between different feature channels.

Figure 4 illustrates the overall structure of LS-Unilab proposed in this paper. The model is an encoding–decoding structure. The encoder is used to extract the semantic information of slow-moving landslides, and the decoder is used to locate the target and identify the boundary.

2.3.1. Uniformer Block

On the basis of the multi-head relation aggregator, local and global relation aggregators are used to extract multiscale features of landslides. Uniformer [36,37] is a visual model based on the Transformer architecture, which includes key modules, such as dynamic position embedding (DPE), the multi-head relation aggregator (MHRA), and the feed-forward network (FFN).

Given an input deformation rate sample feature $X_{in}\in R^{4\times 3\times 256\times 256}$, where 4 is the batch size, 3 represents the number of input channels (e.g., annual velocity rate, coherence and DEM), and 256 × 256 is the spatial resolution of the input images the processing involves the following:

$$\begin{array}{l}X=DPE\left(X_{in}\right)+X_{in},\\ Y=MHRA\left(Norm\left(X\right)\right)+X,\\ Z=FFN\left(Norm\left(Y\right)\right)+Y.\end{array}$$

(10)

InSAR data usually have high spatial resolution and rich surface information. MHRA is able to efficiently extract landslide features at different scales in the InSAR sample set and adapt to multiscale information by designing local and global labeling affinity mechanisms at shallow and deep levels, and by organically integrating the convolution operation with a self-attention mechanism. The MHRA can be described as:

$$\begin{array}{c}{R}_n\left(X_{in}\right)=A_n{V}_n\left(X_{in}\right)\\ MHRA\left(X_{in}\right)=Concat\left(R_1\right(X_{in});R_2(X_{in});\cdots ;R_N(X_{in}\left)\right)U\end{array}$$

(11)

where the input feature $X_{in}\in R^{4\times 3\times 256\times 256}$ is first reshaped into a token sequence $X_{in}\in R^{L\times C}$, where $L=B\times H\times W$, $R_n\left(\cdot \right)$ denotes the aggregation relation of the nth head, and $U\in R^{C\times C}$ is a learnable parameter matrix for all heads. Each aggregation relation includes token context encoding and token affinity learning. The original token is encoded into a contextual token via a linear transformation $V_n\left(X_{in}\right)\in R^{L\times C/N}$, and the token affinity $A_n\in R^{L\times L}$ is then used to summarize the contextual information and the spatial and logical relationships among the features of each channel of the InSAR sample.

In shallow processing, we use local relational aggregation (LRA) to capture small-scale landslide features in InSAR data by learning local relationships between tokens. Specifically, for a given input feature $X_{in}$, LRA computes the affinity matrix of the token with other tokens in the small domain:

$$A_{local}\left(X_i,X_j\right)=a_{i-j}$$

(12)

where $a_n\in R^{B\times H\times W}$ is a learnable parameter, $X_j$ represents other tokens in the neighborhood, and $(i-j)$ denotes the relative position of the target token to other tokens in the domain. The token in the shallow layer has a smaller neighborhood receptive field, which is particularly suitable for capturing small-scale landslide features and improving the ability to discriminate landslides from noise.

In the deeper network, we introduce global relational aggregation (GRA) for processing larger scale landslide features in the deeper layer. The GRA calculates the global affinity matrix between markers via the following formula:

$$A_{global}\left(X_i,X_j\right)={\displaystyle \frac{e^{Q_n\left(X_i{)}^T{K}_n\right(X_j)}}{\sum _{j^{\prime }} e^{Q_n\left(X_i{)}^T{K}_n\right(X_{j^{\prime }})}}}$$

(13)

Unlike local relational aggregation, $X_j$ can denote any token in the global scope, and $Q_n(\cdot )$ and $K(\cdot )$ are linear transformations of the query and key, respectively, to measure the similarity between any two tokens. With this formulation, MHRA is able to learn large-scale dependencies in InSAR data and identify large-scale landslide features. By combining LRA and GRA, Uniformer can skillfully handle multiscale features in InSAR data, which improves the accuracy of feature extraction and maintains high computational efficiency.

The DPE block encodes the position of the input InSAR sample features via deep convolution, giving the model the ability to handle arbitrarily resolved inputs. The FFN, similar to traditional ViT (Vision Transformer) [38], consists of two linear layers and a nonlinear function (e.g., GELU), enhancing feature expression in the InSAR sample set.

In this network, the input feature layer is transformed by four embedding processes (Figure 4), and the initial three channel numbers are converted to [64, 128, 320, 512]. Each embedding layer is processed by a Uniformer module with depth $L_i$, where the $L_i$ values are [3,4,8]. The first stage converts the number of input channels from 3 to 64 and processes them via 3 LRA modules. The image size changes from 256 × 256 to 64 × 64. The second stage increases the number of channels to 128, reduces resolution to 32 × 32, and uses 4 LRA modules. The third stage increases the number of channels to 320, lowers the resolution to 16 × 16, and uses 8 GRA modules. The final stage increases the number of channels to 512, reduces the resolution to 8 × 8, and uses 3 GRA modules.

Each stage progressively enhances the depth, channel and resolution to extract features at different levels. The shallow level of the network extracts small-scale features of landslides through LRA, especially to capture details and suppress noise, whereas the deeper level extracts large-scale features of landslides through GRA, which is capable of recognizing and dealing with a wider range of terrain variations and complex landslide features.

2.3.2. Channel Attention

In our InSAR sample set, the deformation annual velocity rate map is the main feature of slow-moving landslides, and the coherence map and DEM are used to make landslide boundaries more accurate. To optimize the information representation of different channel features, the weights occupied by the features are adjusted by adding the squeeze-and-excite (SE) channel attention module [39]:

$$z={\displaystyle \frac{1}{H\times W}}\sum _{i=1}^H \sum _{j=1}^W X\left(i,j\right)$$

(14)

$$s=\sigma \left(W_2\cdot ReLU(W_1\cdot z)\right)$$

(15)

$$\tilde{X}=s\cdot X$$

(16)

where $X$ denotes the channel of the input feature map and, after generating the descriptor $z$, the weight $s\in R^{1\times 1\times C}$ for each channel is generated through a two-layer fully connected network; $W_1$and $W_2$ are the weight matrices of the outputs of each layer, and $\sigma (\cdot )$ is a sigmoid function used to normalize the output weight $s$. The generated channel weights $s$ act on the input feature maps to obtain $\tilde{X}$. After multiscale feature extraction via the Uniformer block, the decoder layer uses the multibranch ASPP (atrous spatial pyramid pooling) module to process the deep features obtained from the fourth stage of the feature extraction module through multiscale dilation convolution and global pooling [35]. The semantic information at different scales in the InSAR sample set can be effectively captured by adopting expansion rates of 6, 12 and 18 for the 3 × 3 dilation convolution. The low-level detail features obtained from the second stage Uniformer block are processed by a convolutional layer to reduce the number of channels, and then fused with the up-sampled deep semantic information features to preserve the details and global semantic information. Finally, the fused feature maps are restored to their original resolution by categorizing the convolutional layers and up-sampling operations to generate accurate semantic segmentation results.

2.4. Accuracy Evaluation

Three commonly used metrics are selected to quantitatively evaluate the performance of the surface deformation anomaly assessment model: recall, intersection over union (IoU) and precision. These three quantitative evaluation metrics are defined as follows.

Recall is the proportion of all actual positive instances that are correctly recognized as positive instances. Intersection over union is used to measure the degree of overlap between the predicted boundary and the true boundary, and precision is usually defined as the ratio of the number of correctly classified pixels to the total number of positive pixels. The formulas for the recall, intersection ratio and pixel accuracy are shown below:

$$Recall={\displaystyle \frac{TP}{TP+FN}},IoU={\displaystyle \frac{TP}{FP+TP+FN}},Precision={\displaystyle \frac{TP}{TP+FP}}$$

(17)

where $TP$ represents the number of true positive pixels correctly classified as landslides. $TN$ represents the number of true negative pixels correctly classified as background. $FP$ denotes the count of background pixels mistakenly classified as landslides, and $FN$ signifies the count of landslide pixels incorrectly categorized as background. The detection results are dichotomized, and any region outside the “safe” region is identified as a detected slow-moving landslide. All the metrics are calculated on the basis of the number of pixels recognized and the total number of pixels in the test dataset.

3. Study Area and Data

3.1. Study Area

The area around Kangding city and the middle reaches of Jinsha River Gorge were selected as study areas (Figure 5). Kangding is a city frequently affected by earthquakes because of its proximity to the Xianshuihe fault, which makes it prone to earthquake-triggered landslides. On 5 September 2022, an Ms. 6.8 earthquake struck Luing County near Kangding, which triggered 5336 landslides with a total area of 28.53 km² [40]. The Jinsha River valley is located on the southeastern edge of the Tibetan Plateau at coordinates 29°58′N–31°1′N and 98°31′E–99°14′E. The geology is particularly complex, with high mountains and deep valleys in the form of a V-shape [41]. Considering the complex geology, strong crustal movements, frequent summer rainstorms and human activities (deforestation), the Jinsha River is an area where geohazards, such as mudslides and landslides, are widely distributed and highly active [42,43]. Figure 6 shows field photographs of slow-moving landslides along the Jinsha River, where slow-moving landslides are marked with red rectangular boxes.

3.2. Data

In this study, the multisource data used for hazard assessment of slow-moving landslides are presented in

Table 2. Remote sensing images and geological data.

Data Type	Source
Remote sensing	Sentinel-1
	Sentinel-2
	Google Earth Images
SRTM1 DEM	http://step.esa.int/auxdata/dem/SRTMGL1/ (accessed on 16 May 2024)
Fault distribution	https://docs.gmt-china.org/latest/dataset-CN/CN-faults/ (accessed on 16 May 2024)
NASA landslides list	https://gpm.nasa.gov/applications/landslides(accessed on 16 May 2024)
Field investigation	-

.

The Sentinel-1 SAR data are used for InSAR processing. Sentinel-1 is an all-weather, all-day radar imaging system operating in the C-band, with a 12-day repeat cycle. In its interferometric wide (IW) mode, it offers a resolution of 5 min range and 20 m in azimuth. The Sentinel-1 data used are listed in

Table 3. Sentinel-1 data parameters used in the study.

Path	Number of SAR Images	Ascending/Descending	Acquisition Time
99	20	Ascending	2022.1–2023.9
26	91	Ascending	2022.1–2023.9

. A total of 111 SAR images from three tracks are used for the experiment, covering the period from January 2022 to September 2023.

For InSAR processing, the SRTM1 DEM with 30-m resolution is used for InSAR topographic phase removal, atmospheric phase correction, geocoding and sample set input.

In accordance with the InSAR sample set production method described in Section 2.2, the initial deformation regions are evaluated via Sentinel-2 optical imagery and Google Earth optical imagery from September 2023, landslide inventory data from NASA, and field surveys. The SRTM1 DEM is then resampled to the same 40-m resolution as the InSAR results and used for refining the deformation region boundaries. The slices are sliced according to a size of 256 × 256 with 50% overlap, and the slices with positive samples are retained. To solve the problem of imbalance in the proportion of positive and negative samples due to the small number of real landslide areas, we used rotate, pan, flip and mirror operations to extend the dataset. Eventually, 5660 landslide assessment samples were obtained. The training and test sets were divided at a ratio of 8:2.

The experiments are based on a single NVIDIA RTX 3090 GPU. In the experiments, the batch size was 8 and 350 epochs were trained. The optimizer is adaptive moment estimation (Adam). In this setup, the initial learning rate is 0.0005 and the minimum learning rate is 1% of the initial learning rate. The learning rate is decreased via cosine annealing (‘cos’), which means that the learning rate gradually decreases according to the cosine function during the training process, resulting in a smooth learning rate decay. We used the same parameters in the ablation experiments as in the comparison experiments.

4. Results

4.1. Displacement Maps Derived from the Refined InSAR Method

On the basis of the refined SBAS-InSAR method, we obtain the deformation velocity map of the PATH 26 track of the Sentinel-1 image. The images before and after atmospheric correction are shown in Figure 7. The red dot in the lower right of the image indicates the area affected by the Luding earthquake in Ganzi Prefecture, Sichuan Province, on 5 September 2022. In this image, we set the window size of the azimuth and ranging directions to 100 × 100 pixels with an overlap rate of 30% (the calculated decimal part is rounded to the nearest integer), which is due to the size of the maximum deformation in the Luding earthquake disaster area (the field distance is 4-km). For Sentinel 1 orbital path 99, we used 89 × 89 pixels as the window size of the azimuth and ranging directions with an overlap rate of 30%, which is due to the size of the maximum deformation (the field distance is 3.6-km). Elevation correction is performed directly on the entire image. Figure 7b shows the results of elevation correction, revealing significant residual atmospheric errors in region A. Figure 7c displays the results of adaptive window correction, highlighting numerous small deformation features in region B. In region D, there is less atmospheric error in the original image (Figure 7a), but the result of elevation correction results in more uplift deformation (Figure 7b). An analysis of the overall results reveals that the adaptive window method achieves better results. Yu et al. used the standard deviation (std) of the correction result as an indicator to evaluate atmospheric correction [44]. In those results (Figure 7d), the adaptive window correction reduces the standard deviation from 14.96 mm/y to 13.38 mm/y, making the distribution of deformation results more consistent with the standard normal distribution, whereas the elevation correction increases to 16.82 mm/y.

A total of 111 SAR images were processed to derive the annual deformation velocity map, as depicted in Figure 8. On the basis of the analysis of the annual deformation characteristics, the main subsidence areas of Kangding city are distributed in eastern, northeastern and southern parts of the city. Among them, I, II and III exhibit significant landslide deformation. Optical image data and research results show that landslides have been confirmed to exist in region I [45]. The snow cover area in region II is greater, accompanied by an increased occurrence of false alarm signals for landslides. The Luding earthquake occurred close to region III where the landslide distribution is dense and concentrated. The major subsidence areas in the Jinsha River Gorge are located on both sides of the river. In areas IV, V and VI, many areas of landslides may hinder the Jinsha River channel. Figure 8c is a zoomed-in view of area IV in Figure 8b, where the pink dots indicate the geographical location of the field survey, and its corresponding field photographs are displayed in Figure 6 as P1P6.

4.2. Ablation Experiments and Comparison Experiments of Models

We used the classical semantic segmentation models Unet [46], Swin Transformer [47], and SegFormer [48] as model comparison experiments, and performed ablation experiments on the two improvement points proposed in LS-Unilab. The results of the ablation and comparison experiments are shown in

Table 4. Quantitative evaluation results of landslide hazard assessment.

	Model	IoU	Recall	Precision
1	Deeplabv3+	74.12%	94.52%	77.45%
	SE-DeeplabV3+	77.65%	94.96%	80.99%
	Unilab	82.86%	94.62%	86.96%
	LS-Unilab	84.10%	96.16%	87.02%
2	Unet	79.83%	94.45%	83.75%
	SwinTransformer	74.92%	86.60%	84.73%
	SegFormer	80.53%	92.17%	86.45%

and Figure 9.

As shown in Figure 9 for the results of the comparison experiments, the existing model can meet the basic requirements for slow-moving landslide assessment. Areas (1) and (2) are located along the Jinsha River, areas (3) and (4) are in Kangding city, and the two study areas are different in terms of topography and geomorphology. For the multiscale problem of landslides in the InSAR results, as shown by the white solid circles in Figure 9 (1) and (2), the small-scale slow-moving landslides in area (1), which are misdiagnosed as “safe” by existing models, are effectively evaluated by the LS-Unilab model, and the boundaries of the large-scale slow-moving landslides in area (2), which are beyond or smaller than the true boundaries, can be accurately explained by our model. In the white circular area in region (3), our model has a rare missed detection. In region (4), the evaluation in the white circular box outperforms the other models.

According to the quantitative evaluation results in Table 4, LS-Unilab has the highest score for each metric. Among them, the IoU of the Swin transformer in the test set is less than 75%, whereas that of SegFormer is rather than 80%, and we can see from Figure 9 (1) and (2) that the former outperforms the latter at small scales, but falls short in landslide assessment at larger scales, resulting in a lower IoU. Ablation experiments for the LS-Unilab model designed in this paper were conducted to demonstrate the effectiveness of the Uniformer and SE block components, The comparison results are shown in Part 1 of Table 4. The ablation of the SE module resulted in a decrease in the performance of the model, with a decrease in the IoU of 1.24%, which proves the effectiveness of the channel attention mechanism in controlling the weight of each channel. The ablation of the Uniformer module leads to a decrease in the model’s IoU by 6.45% and recall by1.20%, which demonstrates that the combined use of local and global information by the Uniformer module improves the model’s understanding of landslide features.

4.3. Landslide Hazard Assessment Results

We superimposed the annual deformation velocity map, coherence map, and DEM, and then cropped them into a series of 256 × 256 overlapping region images to prevent information loss, similar to how the training sample was prepared. The model assessment results are shown in Figure 10 with a total of 778 slow-moving landslides. The hazard level assessment results are summarized in Figure 10c,d. The Kangding city area is prone to frequent earthquakes, resulting in relatively high hazard levels. Previous studies have shown that many slow-moving ancient landslides have occurred along the Jinsha River, leading to large “low” areas.

We further analyzed the morphological characteristics of these identified slow-moving landslides via optical images to verify their accuracy. These morphological features are the result of many years of movement of active landslides and are evidence of long-term slow movement [49]. We randomly selected two areas as validation zones.

Figure 11 and Figure 12 depict the verification areas referenced in Figure 10. Panel (a) presents the annual deformation rate results, with red-framed regions marking areas evaluated as slow-moving landslides. Panel (b) displays the Sentinel-2 optical images, where black areas represent cloud cover. Panel (c) shows the model evaluation results, with “safe areas” depicted as transparent and red-framed regions indicating the InSAR dataset. Panel (d) provides the evaluation results when threshold segmentation is used.

In regions A–D in Figure 11c and Figure 12c, the boundaries of our landslide assessment results exhibit a relatively consistent alignment with the boundaries of the ground truth data. These findings indicate that the model can accurately predict the landslide boundary. We used optical images to verify the authenticity of the landslides in the assessed area (Figure 12b). Although there are clear soil flows in the optical image in region C, the area in Figure 12a falls within 2 to 3 times the standard deviation of the deformation rate, leading us to classify it as an old landslide with a “low” hazard level. In contrast, although there are fewer landslide features visible in the optical image of region D, certain areas exhibit deformation rates that are less than three times the standard deviation, indicating their classification as old landslides with significant sliding tendencies.

To show the evaluation accuracy of the model, we used the same threshold value as the sample set for threshold segmentation. In Figure 11d and Figure 12d, the red line is the boundary of the landslide evaluation results of the model. The results show that the threshold segmentation method has a limited ability to discriminate “low” hazard levels and fails to accurately determine the boundaries of landslides. In addition, in Figure 12d, the threshold segmentation method identifies “moderate” to “very high” landslide hazard levels on the left side of region C, whereas our model does not yield any results. The utilization of a lower threshold to mitigate the impact of low deformation rate signals may inadvertently result in the omission of certain areas with low landslide hazards. In region E of Figure 12c, the landslide assessment region is smaller than the ground truth boundary and it identifies two separate landslides as a single complete landslide. This discrepancy may be attributed to the prevalence of zero values in the deformation velocity data, necessitating a more comprehensive assessment using optical images and geomorphic information.

5. Discussion

5.1. Validation of Landslide Assessment Based on Time-Series Results

To verify the accuracy of landslide hazard assessment, we apply the time-series deformation results for displacement process analysis. Figure 13 presents the time-series deformation analysis of selected landslide regions containing validation points (marked in Figure 10), with several monitoring points established for detailed investigation. The selection of these points was based on the central axis of the overall landslide structure. By choosing points along this axis, we aimed to ensure that they represented the most critical areas of deformation, particularly where landslide activity was most prominent. Points in area (a) showed an overall downward trend from January 2022 to September 2023, with a maximum deformation of −120 mm. Optical imagery reveals distinct landslide features at point J04; however, the deformation rate indicates a slow-moving slide, suggesting an ancient landslide that still exhibited sliding tendencies. The linear distribution of monitoring points in region (b) exhibited a pronounced trend of rapid subsidence, indicative of major landslide occurrence within the area. This observation aligns with the “very high” classification obtained from the hazard assessment conducted for this region. The solifluction flows can be clearly observed in area (c) from the optical image. On the basis of the findings from the time-series analysis, it is determined that the old landslide is still undergoing active sliding. The monitoring points in area (d) showed significant fluctuations in deformation, especially at points K03 and K04, with deformations ranging from −120 mm to 50 mm. The scatter plot shows a sinusoidal variation in deformation, whereas the overall trend indicates landslide activity. The analysis of deformation changes at monitoring points verified the accuracy of the deformation and model identification results.

5.2. Validation of Landslide Assessment Based on Existing Work

The results obtained by our model confirm that we can obtain a more comprehensive landslide hazard assessment in low-coherence regions. Regions A and B in Figure 11 and region D in Figure 12 were chosen as the research objects for comparison with the results obtained by other researchers. Figure 14a corresponds to region A and B, and Figure 14c,e correspond to region D. The field survey images obtained in August 2022 are depicted in Figure 14, where the presence of a substantial number of landslides within the evaluated area agrees with the outcomes derived from the model. Notably, the areas marked by red circles in Figure 14b correspond to regions A, B, and C in Figure 14a, further confirming the validity and accuracy of our method. In the black rectangular box of Figure 14f, some regions do not spatially overlap with reported results in the literature, as indicated by the red circles in (d–f). Due to the temporal differences in the study periods, a detailed analysis of these regions was conducted using optical images (g–k). Among these, regions 1 and 2, corresponding to the circular areas in (d,e), exhibited no significant surface changes, suggesting that they stabilized during the monitoring period. For the newly detected regions 3 to 5, optical imagery revealed the presence of developing mountain fractures, confirming that these areas are in an active phase of landslide development. By comparing the results of landslide detection across different time periods, we can also conduct a comprehensive analysis and make informed judgments regarding the geological movement in this area.

6. Conclusions

In this study, we aimed to assess the hazard of slow-moving landslides in alpine canyon areas and vegetation coverage areas. Therefore, a window-based atmosphere correction method and LS-Unilab landslide hazard assessment model were proposed.

First, we processed 111 scenes of Sentinel-1 images to obtain surface deformation information via the SBAS-InSAR technique. The adaptive window atmospheric correction method was applied to acquire the refined annual deformation rate. Second, we labeled slow-moving landslides on InSAR results with the refinement of landslide boundaries using optical images. The standard deviation of the deformation rate was used as the hazard classification standard of the sliding area to produce the InSAR hazard assessment sample set. Finally, the LS-Unilab semantic segmentation model based on Transformer was constructed, achieving a recall rate of 96.16%, an intersection over union ratio of 84.10%, and a precision rate of 87.02% for slow-moving landslide hazard assessment in the testing dataset.

Our method was applied to Kangding city and Jinsha River Gorge. A total of 778 slow landslide hazard areas were detected, 299 of which were in Jinsha River Gorge, accounting for 11.3% of the “very high” areas and 10.4% of the “high” areas, and seven were likely to block large-scale landslides in the river. There are 479 landslides in Kangding city, concentrated east of Kangding city because of the Luding earthquake. Among them, “very high” areas accounted for 27.4%, and “high” areas accounted for 14.9%. Geological disaster investigation agencies and local departments need to conduct real-time monitoring of large landslides along the Jinsha River to prevent major disasters.