Interferometric Synthetic Aperture Radar (InSAR)-Based Absence Sampling for Machine-Learning-Based Landslide Susceptibility Map**: The Three Gorges Reservoir Area, China

Abstract

The accurate prediction of landslide susceptibility relies on effectively handling landslide absence samples in machine learning (ML) models. However, existing research tends to generate these samples in feature space, posing challenges in field validation, or using physics-informed models, thereby limiting their applicability. The rapid progress of interferometric synthetic aperture radar (InSAR) technology may bridge this gap by offering satellite images with extensive area coverage and precise surface deformation measurements at millimeter scales. Here, we propose an InSAR-based sampling strategy to generate absence samples for landslide susceptibility map** in the Badong–Zigui area near the Three Gorges Reservoir, China. We achieve this by employing a Small Baseline Subset (SBAS) InSAR to generate the annual average ground deformation. Subsequently, we select absence samples from slopes with very slow deformation. Logistic regression, support vector machine, and random forest models demonstrate improvement when using InSAR-based absence samples, indicating enhanced accuracy in reflecting non-landslide conditions. Furthermore, we compare different integration methods to integrate InSAR into ML models, including absence sampling, joint training, overlay weights, and their combination, finding that utilizing all three methods simultaneously optimally improves landslide susceptibility models.

1. Introduction

Landslide susceptibility is defined as the likelihood of landslide occurrences under given geo-environmental conditions. It represents a fundamental step toward assessing landslide hazards and develo** mitigation strategies, especially in mountainous areas, where geological hazards are often concealed, harmful, and ubiquitous due to the high altitude, steep terrain, and dense vegetation [1]. China has complex geological environments, with a large proportion (>65%) of mountain areas prone to landslides [2]. Seismic activities, climatic extremes, and the increasingly acute human disturbances from unplanned urban expansions and deforestation frequently occur in the mountains, triggering substantial major landslides and subsequent casualties and economic losses [3,4]. Such socioeconomic burdens have prompted a significant surge in research papers on landslide susceptibility map** (LSM) [5]. Among these, machine learning (ML)-based LSM has gained widespread popularity due to its high performance in addressing challenges and improving accuracies for capturing the spatial locations of landslides and delineating areas prone to potential landslides [6,7].

ML-based landslide susceptibility maps are premised on the reliability of two types of samples: samples describing conditions where landslides have occurred (landslide presence data) and samples representing conditions where no landslides have been observed (landslide absence data) [8]. Landslide presence data are usually generated from landslide inventories produced by field investigations and the interpretation of aerial photos or satellite remote sensing images. Landslide absence data are sampled from areas other than those described in landslide presence data using various strategies [9]. Incorporating landslide absence data properly is crucial because it can appropriately control the expansion of high-susceptibility landslide areas [10]. Various absence sampling methods have been proposed to improve the accuracy and quality of LSM. They can be broadly classified into two types: geographic and feature-based generations [8]. Classic geographic generation methods, such as random sampling and buffer-controlled analysis, collect landside absence data from areas distant from known landslides. However, they cannot ensure that the selected areas are entirely free from landslides, since the known landslide inventory may not cover all locations where landslides have occurred. Feature-based generation assumes that future landslides are more likely to occur under environmental conditions similar to those causing past landslide events. Following this principle, landslide presence and absence data can be separated in the feature space. Target space exteriorization sampling and utilizing Mahalanobis distances are two examples of this approach [11,12]. However, the generated absence data acquired with these methods are virtual and exist only in the feature space, making field validation difficult.

Recent studies have indicated that involving physical laws and expert knowledge in ML methods enhances LSM accuracy and interpretability. Several investigations have utilized physical models to better select landslide absence samples. For instance, Wei et al. (2021) applied the TRIGRS model to convert the spatially discrete lithology information into the continuous safety factor (Fs) within a specified range and pre-selected landslide absence samples by limiting Fs [13]. Then, a convolutional neural network (CNN) model was trained using historical landslide presence samples and physics-informed landslide absence samples. They concluded that compared to the traditional CNN model, this hybrid approach performs better, enhancing CNN’s adaptability to a new study area. Liu et al. (2023) employed a Scoops 3D model for qualitative slope stability analysis, extracting landslide absence samples from areas with Fs greater than 1.5 [9,13]. The generated landslide absence samples are more rational than those derived from random sampling. However, one major drawback of physics-informed methods is that they require a deep understanding of the physical processes governing slope stability and various geotechnical parameters, such as soil depth, shear strength parameters, and groundwater configurations [9]. Obtaining such comprehensive knowledge is challenging and costly, and these limitations force researchers to rely on approximate data, which introduces uncertainties, especially in complex and regional settings. Therefore, finding a more straightforward and interpretable approach that generates landslide absence samples at larger scales is imperative.

The rapid progress of remote sensing technologies, particularly interferometric synthetic aperture radar (InSAR), presents a promising opportunity to bridge the aforementioned gap. InSAR deformation data, known for their extensive area coverage and precise surface deformation measurements at millimeter scales, can effectively capture slope activities such as ground uplift and subsidence [1,14]. This capability facilitates the identification of potential landslide-prone areas [15]. A comparative analysis between ML-based LSM and InSAR deformation results reveals spatial congruence, with higher deformation points coinciding with elevated susceptibility levels and more stable points in regions with lower susceptibility levels [16]. In this context, we may assume InSAR’s ability to distinguish between stable and unstable slopes helps select landslide absence samples. Indeed, previous studies have highlighted the enhanced accuracy of landslide susceptibility models by integrating InSAR technology. Standard integration methods based on InSAR data involve updating landslide inventories [17], cross-validating traditional LSM [16], joint training for landslide conditioning factors [18], and weighted overlays of results either through a contingency matrix [19,20] or an empirical formula [3,15]. However, there is a gap in the literature regarding the exploration of InSAR deformation data for selecting landslide absence samples.

The Three Gorges Reservoir Area (TGRA) in China is recognized as landslide-prone due to its fragile geological conditions, seasonal rainfall, and periodic reservoir filling [19]. Given that more than 5000 landslides have occurred in the TGRA, there is an urgent need for a prompt and accurate LSM assessment. The escalating volume and frequency of remote sensing data, along with advancements in computing platforms for large datasets, create opportunities to explore more accurate LSM assessments with more advanced absence sampling strategies. This study proposes an InSAR-based sampling strategy for generating landslide absence samples. Three commonly used ML methods, including logistic regression (LR), random forest (RF), and support vector machine (SVM), are applied to construct the LSM models, mitigating the randomness associated with a single ML type [21]. To validate the effectiveness of our proposed sampling method, the ML models using landslide absence samples from InSAR-based sampling are compared against those employing the widely used random sampling. The novelty of our study lies in expanding our understanding of integrating InSAR deformation data with LSM and elucidating InSAR’s role in sampling landslide absence data.

2. Materials and Methods

2.1. Study Area and Landslide Inventory

This study focused on the Badong–Zigui section of the TGRA, China, along the Yangtze River (Figure 1). The boundary of the study area is defined by merging a cluster of sub-watersheds, with geographical coordinates spanning from 110°17′ to 110°40′E and 30°56′ to 31°05′N. This section is situated within a low-to-medium mountainous area characterized by deep canyons [22]. The geological structure is complex, with faults and fragmented rock mass. The faults include the ** Fault [23] (Figure 2). The Quaternary Sinian strata are widely exposed, along with a few much older strata, such as Devonian and Silurian ones [24]. The lithology is mostly carbonate, sand shale, marlstone, and mudstone, which is sensitive to landsliding [25].

The study area has a humid subtropical climate. The average annual precipitation accounts for 1100 mm, with approximately 70% occurring from May to September. The reservoir water level of TGRA was initially impounded from 69 to 135 m a.s.l. in 2003. It has been undergoing cyclic hydraulic operations since 2009, fluctuating between 145 and 175 m a.s.l. The rapid rise and drawdown of water levels in the TGRA keep modifying the hydrogeological environment of the reservoir banks and promoting bank instability [26]. As a result, many landslides have occurred on slopes adjacent to the reservoir. Based on high-resolution remote sensing imagery, field investigations, and landslide reports, we identified 134 landslides covering an area of 17.77 km², accounting for 5.06% of the entire study area. Most of them are colluvial landslides, with seasonal rainfall, periodic fluctuation of water level, and anthropogenic activities as the main triggering factors of landslides. These landslides exhibit slow movement, ranging from millimeters to several meters per year, persisting over extended years [19]. Figure 3 illustrates two typical landslides and their micro-geomorphic features, indicating ongoing deformation in the study area.

2.2. Methodology

This study was conducted in four stages (Figure 4). Firstly, layers of 13 landslide conditioning factors were prepared for the study area, and their multi-collinearity was checked. Secondly, the Small Baseline Subset (SBAS) InSAR technology was applied to produce annual average ground deformation and select landslide absence samples. Thirdly, landslide-conditioning datasets were constructed, considering the uncertainty with the locations of landslide absence samples. Finally, three ML models were trained, their performance was evaluated using metrics, and the study area was classified accordingly into five subareas with different levels of landslide susceptibility using the natural breaks method.

2.2.1. Landslide Conditioning Factors

Landslide conditioning factors (LCFs) include the internal geology and the external environment. There is no uniform criterion for selecting LCFs. Considering data availability and the relevance to the landslide mechanism, we selected 13 conditioning factors: elevation, slope, aspect, plan curvature (Plan_C), profile curvature (Profile_C), topographic wetness index (TWI), distance to faults (D_Fault), distance to rivers (D_River), distance to roads (D_Road), land use, the Normalized Difference Vegetation Index (NDVI), average annual precipitation (AAP), and lithology. All LCFs were discretized according to a raster DEM (digital elevation model) with a spatial resolution of 30 m (Figure 5). The data source is presented in

Table 1. Information on the 13 LCFs.

LCFs	Resolution	Variable Type	Source
Elevation	30 m	Continuous	ASTER GDEM (30-m DEM)
Slope	30 m	Continuous	Derived from the DEM
Aspect
Plan_C
Profile_C
TWI
D_Fault	30 m	Continuous	Adapted and digitized from [23,25]
D_River
D_Road	30 m	Continuous	OpenStreetMap https://www.91weitu.com/ (accessed on 11 October 2023)
Land use	30 m	Discrete	GlobaLand30 dataset https://www.resdc.cn/ (accessed on 11 October 2023)
NDVI	30 m	Continuous	Derived from Landsat 8 images
AAP	1 km	Continuous	1-km annual average precipitation dataset for China http://www.gis5g.com/ (accessed on 11 October 2023)
Lithology	30 m	Discrete	Adapted and digitized from [23,25]

.

Multi-collinearity, where two or more covariates are linearly related, can adversely affect the accuracy of a statistical model by increasing the variance of estimates and impede the actual convergence toward the most suitable solution. This paper applied the variance inflation factor (VIF) and Person’s correlation to check the correlations and multi-collinearity of LCFs (Figure 6). VIF is defined as follows:

VIF = 1/(1 − R_j²)

(1)

where R_j is the negative correlation coefficient between the j-th independent variable (X_j) and the other independent variables (X_i). It is widely accepted that VIF values below 5 indicate the absence of multicollinearity, while values between 5 and 10 suggest weak multicollinearity, and values exceeding 10 indicate moderate or higher multicollinearity [5,12]. On the other hand, there is no universal standard for the extent to which correlation represents high collinearity. We follow a previous study by Fu and Wang et al. (2023) and set the correlation values over 0.8 to indicate a strong correlation between two factors [12]. This study identified two LCFs with weak multicollinearity, elevation and AAP, with VIF values of 5.65 and 6.17, respectively. The correlation analysis also revealed that AAP was highly correlated to elevation (0.89). Therefore, AAP was considered responsible for multi-collinearity and was left out of the subsequent analysis.

2.2.2. InSAR-Based Sampling Strategy

Sentinel-1 SAR C-band data (5.5 cm wavelength) were obtained from the European Space Agency (ESA). We compiled a dense stack of ascending Sentinel-1 images from 6 January 2018 to 17 March 2023 to derive time-series deformation histories. This study used an Interferometric Wide swath mode with VV polarization.

The SBAS InSAR, proposed by Berardino et al. (2002), applies SAR interferograms with small baselines to minimize inaccuracies in topographic data and the impacts of baseline decorrelation [1,2,27,28,29,30,31]. The general processing of SBAS-InSAR starts from SAR image co-registration (Figure 4). Then, precise orbit files were downloaded from the ESA and preprocessed to reduce the orbital error. With a perpendicular baseline threshold of 200 m and a temporal baseline of 90 days, we generated 441 interferograms. From this set, we carefully selected 291 interferograms with good coherence for subsequent analysis. Next, adaptive filtering and phase unwrap** were performed to acquire a series of unwrapped differential interferograms. Then, we applied a simple linear model to eliminate the topography-dependent atmospheric phase. This procedure is commonly referred to as atmospheric correction [29]. Finally, the deformation rates were obtained through the singular value decomposition (SVD) and geocoded using the DEM.

This study evaluated slope deformation rates along the Long-of-Sight (LOS) direction of Sentinel-1 images. Cruden et al. (1996) suggested that landslides could be classified into seven types based on their deformation rates: extremely slow, very slow, slow, moderate, rapid, very rapid, and extremely rapid [32]. In this respect, landslide absence samples from areas with lower surface deformation rates are theoretically more stable than their neighboring regions. The LOS deformation measurements in this study were projections of three-dimensional (3D) deformation in the actual world onto the 1D LOS direction. Therefore, negative deformation values indicated movements away from the satellite, while positive deformation values indicated those toward the satellite. While comprehending 3D deformation offers a better understanding of surface processes, obtaining the 3D deformation requires multi-orbit satellite SAR observations from both ascending and descending passes. Because this study only had access to ascending Sentinel-1 images, we relied on LOS InSAR measurements, following methodologies employed in relevant studies [29,33,34,35].

Predefined deformation thresholds were assigned to differentiate landslide-prone areas and safe areas. The choices of thresholds are case-specific and depend on the mechanical properties of the failed material, failure mechanism, sensor measurement precision, and investigation objectives [36,37]. This study addressed this issue statistically, considering the in-phase response between landslide susceptibility and deformation rates identified in previous studies [16]. The Kriging tool in ArcGIS 10.8 was employed to interpolate the SBAS measurement points, fill in the no-data areas associated with the decorrelation of SAR images due to the steep terrain and dense vegetation, and generate the deformation rates for the entire study area. Then, the absolute values of deformation rates were reclassified into five levels using the natural breaks method: very low, low, moderate, high, and very high. Areas with very low deformation rates are considered safe areas from which landslide absence samples can be selected.

2.2.3. Construction of Datasets for LS Modeling

We adopted a cross-validation strategy to partition landslides for training (n = 67) and those for testing (n = 67), as indicated by Figure 1b. Landslide records were discretized into 30 m grids corresponding to the DEM used. Then, all grid cells on the landslides were used to construct landslide presence samples, resulting in 9818 cells for training and 9962 cells for testing.

Landslide absence samples equal to landslide presence samples were randomly extracted from the predefined safe area. The pool of available samples for selection significantly exceeded the required quantity. The uncertainty inherent in selecting landslide absence samples leads to diverse training datasets and the development of varied trained models [38]. To assess uncertainty in selected sample locations, we generated 100 sets of diverse landslide absence samples for training, each derived from partially distinct spatial positions. After combining these absence samples with the landslide presence samples, we constructed 100 training datasets and, thus, 100 training models. Then, the trained models were evaluated using various metrics against a consistent testing dataset whose landslide presence samples were from 67 testing landslides, and land-slide absence samples were from the InSAR-defined safe area, all with fixed locations.

2.2.4. Landslide Susceptibility Prediction Modeling

Logistic regression (LR) is one of the most widely utilized methods in multivariate statistical analysis. It predicts the probability of event occurrences or classifies a dependent variable by establishing the log-odds for the event as a linear combination of predictor variables. The maximum likelihood method is used to calculate the regression coefficients, followed by a gradient iteration to estimate the parameters. Logistic regression is mathematically expressed as

P = 1/(1 + e^−z)

(2)

Z = β₀ + β₁X₁ + β₂X₂ + … + β_nX_n

(3)

where P is the probability of the event occurrence ranging from 0 to 1, Z is the sum of all input variables based on weight values, β_i (i = 1, 2, 3…n) is the regression coefficient based on the training samples, n is the number of independent variables, and X_i (i = 1, 2, 3…n) is the independent variable.

Support vector machine (SVM) can solve problems with high-dimensional, small samples, and nonlinear features, and is therefore frequently used for LSM [39]. The core idea of SVM is the use of kernel functions to search for the best classification hyperplane in the feature space to maximize the interval between positive and negative samples on the training dataset [12]. Once no hyperplane in the current feature space can linearly divide the sample categories, the original feature space can be mapped to a higher-dimensional feature space by a kernel function to make the samples linearly divisible. The selection of the kernel function strongly influences the classification accuracy of SVM. Commonly used kernel functions include linear (LN), polynomial (PL), radial basis function (RBF), and sigmoid (SIG) kernels. This study applies RBF-SVM to calculate the landslide susceptibility index because it can perform better than other kernel functions in non-linear classification [39].

The basic concept of random forest (RF) is to develop multiple decision trees on random subsets [40]. It uses bootstrap aggregation to select samples from the original dataset to train a tree. Each tree votes independently. Then, the voting results of all decision trees are integrated, and the category that appears most in all decision trees is considered the final prediction result. This bootstrap** procedure leads to better model performance because it integrates several decision tree models to obtain the optimal prediction results, decreasing the variance of the model without increasing the bias [41,42]. In general, the number of decision trees (n_estimators) and the maximum ratio of features (max_features) are the most critical parameters influencing the accuracy of RF models [8]. A large value of n_estimators will improve the modeling accuracy, but it also increases the complexity and modeling time [43]. The max_features parameter limits the number of features to be branched. In this study, we applied the grid search method to find the optimal hyperparameters for the RF model. After many trials, the value of n_estimators was set to 100, and that of max_features was set to 4.

2.2.5. Model Performance Evaluation

We selected six metrics to assess the predictive performance of the model. Firstly, the area under the curve (AUC) was derived from the receiver operator curve (ROC) to evaluate the comprehensive performance of the models [43]. Then, a confusion matrix was used to assess the performance of a classification model by comparing the predicted class labels with the actual class labels. The matrix contains four terms: true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). From the confusion matrix, five metrics can be calculated, including accuracy, precision, recall, specificity, and F1-score, as follows:

accuracy = (TP + TN)/(TP + TN + FP + FN)

(4)

precision = TP/(TP + FP)

(5)

recall = TP/(FN + TP)

(6)

specificity = TN/(TN + FP)

(7)

F1 score = 2(precision × recall)/(precision + recall)

(8)

3. Results

3.1. InSAR-Defined Safe Areas

The maximum annual average deformation rate, determined through SBAS-InSAR technology, reached 67.7 mm/yr (Figure 7a). After conducting the Kriging interpolation (Figure 7b), the absolute values of annual deformation rates were classified into five grades using the natural breaks method (Figure 7c). These grades represented 0.527%, 5.47%, 18.86%, 37.22%, and 37.92%, respectively, from high to slow deformation. We designated areas with the slowest deformation grade (i.e., deformation rate < 2.4 mm/yr) as safe areas for selecting non-landslide samples. However, within some areas designated as landslides, cells with the slowest deformation also exist. To avoid contradictions and overlaps between landslide and non-landslide samples, we excluded cells with the slowest deformation within landslides from the designated safe areas (Figure 7d).

3.2. Enhanced Landslide Susceptibility Maps with InSAR-Based Absence Sampling

Figure 8 contrasts landslide susceptibility maps from different ML models employing InSAR-based sampling and random sampling, respectively. The study area, divided into five subareas based on distinct levels of landslide susceptibility using the natural breaks method, exhibits similar spatial patterns with higher-susceptibility areas consistently appearing in a strip along both banks of the Yangtze River. This spatial pattern indicates that the fluctuation of reservoir water elevation significantly influences bank-slope stability [30]. Moreover, these areas are interconnected with densely populated regions, frequent engineering activities, and vegetation depletion, which collectively exacerbate landslide susceptibility [15].

A subsequent examination of the ROC curves reveals that the susceptibility maps using the InSAR-based sampling strategy outperform those using random sampling by approximately 1% in AUC in Figure 9a. This indicates that our proposed sampling strategy performs better for all models. In addition, the InSAR-derived landslide absence sampling improved model performance by increasing the percentages in the ‘high’ and ‘very high’ classes or decreasing those in the ‘low’ and ‘very low’ classes (Figure 10a). Figure 10b further illustrates that after applying the InSAR-based sampling strategy, there was an increase in landslide percentages within the ‘high’ and ‘very high’ classes. However, it is important to note that the landslide density remains lower for the LR and SVM models in these classes (Figure 10c), irrespective of the sampling strategy used. This may imply that compared to RF, the LR and SVM models tend to incorporate more landslide grids at the cost of overestimating the percentages of ‘High’ and ‘Very High’ landslide susceptibility classes, decreasing their landslide density.

3.3. Comparison of Different InSAR Integration Methods and Their Combined Use

We compared the effectiveness of the InSAR-based sampling strategy against two classic InSAR integration methods: joint training and overlay weights methods. The joint training method incorporates InSAR deformation data as one of the LCFs into the base dataset [15,44]. The overlay weights method assigns weights of 0.61 and 0.39 to the RF model and SBAS-InSAR results, respectively, and then combines them through the weighted overlay operation [3,19,20,30]. Here, we chose the RF model for integration with InSAR deformation data due to RF’s superior performance, as indicated in Section 3.2.

Figure 11a–c provide a comparative analysis of landslide susceptibility maps generated by RF models that utilize different InSAR integration methods. The AUC metric reveals a range of values between 0.914 and 0.927 for these susceptibility maps. The RF models that integrate InSAR data through joint training and the proposed InSAR-based sampling strategy exhibit comparable model performance. However, the RF model utilizing the overlay weights method stands out, showing the best performance based on the AUC value. We assume that the outstanding performance of the overlay weights method may stem from its ability to fully utilize InSAR deformation information, encompassing both high- and low-deformation regions. In contrast, the InSAR-based sampling strategy only utilizes information from low-deformation areas for selecting non-landslide samples. Despite the joint training method fully using InSAR deformation information, this advantage may be compromised by the unreasonable use of landslide absence samples from random sampling.

We further explore the potential enhancement of LSM through the combined use of multiple InSAR integration methods. The susceptibility maps in Figure 11d–f, generated by RF models combining different InSAR integration methods, demonstrate noteworthy improvements. Combining InSAR-based sampling and overlay weights increases the model’s AUC value to 0.936. Similarly, InSAR-based sampling and joint training raised the AUC value to 0.952. Further, combining InSAR-based sampling, joint training, and overlay weights yields an AUC value of 0.959. A comparison of ROC curves in Figure 9b intuitively demonstrates the improvement achieved by combining three InSAR integration methods. In addition, Figure 12 shows that this combined approach of the three InSAR integration methods increased the percentages of landslide susceptibility classes and landslides in the high and very high categories while maintaining a rational level of landslide density. These results highlight the effectiveness of integrating multiple InSAR methods in enhancing LSM performance and providing more precise and reliable predictions of landslide-prone areas.

4. Discussion

4.1. Uncertainty Analysis of InSAR-Enhanced LSM

ML models are inherently sensitive to data distributions, with their performance potentially inflated by coincidental factors such as data shift, stochastic variance, overfitting, data quality issues, or inadequate validation strategies [5,45]. To address these challenges and quantify the uncertainty of the InSAR-based sampling strategy, this study trained 100 different models with spatially distinct non-landslide samples. The distributions of evaluation metrics are summarized in Figure 13. The results indicate a positive impact of our InSAR-based sampling strategy on the three models compared to random sampling, with an average improvement of ca. 1% across the six selected evaluation metrics, despite the uncertainty in the chosen locations of landslide absence samples. Thus, we conclude that non-landslide samples derived from InSAR-defined safe areas more accurately represent non-landslide conditions than those obtained through random sampling. Additionally, the RF models show superior performance compared to the LR and SVM models. Similar insights can also be found in previous studies, indicating that RF models have more advanced capabilities in terms of accuracy and robustness [15,45,46]. Furthermore, we also assessed the uncertainty associated with the combined use of multiple InSAR integration methods. The results confirm the consistently optimal performance of ML models when combining three InSAR integration methods with evaluation metrics exhibiting a higher mean value, and a lower standard deviation (Figure 14). This suggests that combining multiple InSAR integration methods could further reduce the uncertainty of the InSAR-based absence sampling strategy and facilitate more stable and superior model performance of LSM.

4.2. Empirical Cross-Checking via InSAR Dynamics at the Da** Landslide Group

Accurately identifying credible landslide-prone areas is crucial for managing geological hazards and enables the efficient allocation of public resources and planning by government agencies [47]. Within this context, cross-validations of landslide susceptibility results must be carried out at localities with suspiciously high susceptibility levels. Some studies applied fieldwork to validate their assessment results [3,23,48]. However, a tricky issue associated with slow-moving landslides is that their deformation might span an exceptionally long period, and close monitoring based on field observations would be labor-intensive and resource-consuming. To solve this dilemma, this study explores the use of InSAR to validate high-susceptibility areas generated by LSM. LSM, inferred from machine learning methods, primarily captures the relationship between landslide occurrences and features [45]. In contrast, InSAR-derived deformation intuitively reflects slope dynamics, making it a popular way to identify potential landslides [37,49,50,51]. Thus, we hypothesize that InSAR-derived deformation can validate LSM-identified landslide-prone areas.

The Da** landslide group data were excluded when training the landslide susceptibility models. Thus, the derived insights regarding the LSM accuracy can be validated with InSAR deformation dynamics in this locality. From Figure 15a, we observe prominent deformation features with substantial deformation primarily manifesting in the central sections of the left block (Figure 15b). In particular, this deformation exhibits significant spatial-temporal variations corresponding to the diversified landslide evolution on the landslide surface. These findings highlight the instability of the Da** landslide group and thus support the credibility of the escalated susceptibility level in this locality, as shown in Figure 11a–f.

4.3. Advantages and Future Recommendations

The rapid progress of InSAR technology has opened new opportunities for more accurate LSM. With the free and open dissemination strategy of Sentinel-1 satellite images at high temporal resolution and global coverage, InSAR has become a low-cost landslide monitoring technique, and high-precision Deformation maps can be quickly generated and regularly updated [34,52,53]. Various time-series InSAR technology, including stacking InSAR, SBAS InSAR, and persistent scatterer/distributed scatter (PS/DS)-InSAR, were developed based on multiple-period SAR images with millimeter accuracy [2,19,35,54,55,56,57,58]. In this study, we apply the SBAS InSAR technique to derive deformation information from a dense Sentinel-1 time series because SBAS-InSAR exhibits superior performance in effectively mitigating errors caused by decoherence and atmospheric delay in the spatial baseline [18,31,59,60,61].

The performance of LSM predictive models has been enhanced by integrating InSAR-derived deformation through absence sampling strategy, joint training, overlay weights, or a combination of these methods. However, there is limited research on using InSAR deformation data to select landslide absence samples. This study highlights InSAR’s unique role in absence sampling for LSM by explicitly associating absence samples with low deformation values. While physics-informed models like TRIGRS and Scoops 3D have been explored to select absence samples in areas with higher safety factor values, their application is restricted in smaller areas due to difficulties in obtaining detailed geotechnical parameters such as soil depth, shear strength, and groundwater configurations [9,13,38,62]. In contrast, our proposed sampling strategy can be applied to larger areas, primarily because of unrestricted access to space-born SAR satellite images, especially from Sentinel-1. This broader application makes the proposed strategy more versatile and practical for large-area landslide monitoring and prediction.

In this study, we ensured that the landslide absence samples had clear physical significance by selecting areas with low deformations as defined by InSAR measurements. Our proposed absence sampling strategy achieved better performance than traditional random sampling. However, we acknowledge that our proposed method suffers some uncertainties due to a limited temporal span of SAR data compared to landslide evolution processes. This may result in inconsistent patterns between landslide susceptibility and deformation [3], although this issue is quite common for slow-moving landslides because their lifetime ranges from days to hundreds of years [63]. In addition, integrating SAR images with other remote sensing data, such as high-resolution optical images and Light Detection and Ranging (LiDAR), may aid in selecting landslide absence samples [1,33,50]. High-resolution optical images and LiDAR can intuitively capture the micro-geomorphological features, including landslide fissures, ravines, and tension cracks. Such features play a vital role in understanding landslide dynamics and assessing stability [14,51,64], thereby eliminating false negatives in selecting landslide absence samples. The collaborative use of multi-source data from space-, sky-, and ground-based platforms provides high-quality and complementary information for landslide prevention and warning [1,2,3,30,65,66,67]. Although this method has not been widely applied in the field of LSM, we believe that in the future, it could help optimize the credibility of landslide absence samples and generate more accurate ML-based algorithms for LSM. This direction is critical to advancing our understanding and predictive capabilities in LSM.

5. Conclusions

This study introduced an InSAR-based sampling strategy for generating absence samples in LSM within the Badong–Zigui area in the TGRA, China. Utilizing the SBAS InSAR technology, we computed annual average ground deformation and selected absence samples from areas displaying very slow deformation. We employed three predictive models (LR, SVM, and RF) to verify the efficacy of our sampling strategy for susceptibility zoning. Our findings highlight a positive impact of the InSAR-based sampling on LR, RF, and SVM models, indicating that samples from InSAR-defined safe areas more accurately represent non-landslide conditions than randomly selected negatives. Additionally, we compared the effectiveness of different InSAR integration methods, including absence sampling, joint training, overlay weights, and their combined usage. The results demonstrated the optimal enhancement of landslide susceptibility models from the use of these three InSAR integration methods. Finally, we cross-checked the improved landslide susceptibility maps using InSAR-derived deformation dynamics at the Da** landslide group, consolidating the effectiveness of incorporating InSAR technology into the LSM predictive models.