Comparative Performance of Machine Learning Models for Landslide Susceptibility Assessment: Impact of Sampling Strategies in Highway Buffer Zone

1. Introduction

China experiences some of the highest frequencies of natural disasters. In 2024 alone, these disasters resulted in direct economic losses amounting to 2.5 billion. Landslides, which constitute 96.96% of these natural disasters, are one of the most prevalent forms of disaster in China [1]. Minimizing the consequential impacts of landslides continues to be a significant challenge to policy makers in various fields.

Landslide susceptibility assessment [2], is often associated with landslide susceptibility mapping, takes into account the location of previous landslides and the contributing factors that caused them (such as slope, precipitation, and land use, etc.) within the Areas of Interest (AOIs). It then predicts the distribution of areas that are prone to landslides and draws a zoning map based on the probability of landslide initiation when triggering events occur. Landslide susceptibility mapping has proven to be a highly useful tool in landslide management over the years. The approaches employed to conduct landslide susceptibility mapping can be divided into 3 main categories: inventory mapping, knowledge-based models and data-driven models [3,4].

Prior to the maturity of GIS technology, the assessment was mainly done using inventory mapping. The accuracy of the knowledge-based models, including Analytical Hierarchy Process (AHP) and Weighted Overlay Analysis, reflected on the experts’ knowledge regarding the landslide mechanisms, thus lacking objectivity. The studies based on data-driven models have proliferated in recent years with booming machine learning techniques and more accessible public aerial imagery databases. In landslide susceptibility assessment, machine learning applications are predominantly focused on supervised learning approaches [5], with the most widely adopted models including logistic regression (LR), support vector machines (SVM) [6-8], random forest (RF) [7,9], convolutional neural network (CNN) [10], deep belief network [11], and recurrent neural network (RNN) [6,10].

Developing a universally applicable landslide susceptibility model is challenging for it is observed that the weights of the contributing indicators showed some consistency across models within a single AOI, but significant variability across AOIs [12]. Also, the predictive performance of a specific model can vary under entirely different geological and geomorphological conditions. In light of this, the primary goal of this research is to identify the most suitable model for landslide susceptibility assessment within G3033 Duku Expressway buffer zones. This research aims to examine a range of approaches, from basic models to hybrid methods and more advanced techniques, while also assessing the impact of different sampling strategies on model performance. The findings of this research are expected to contribute not only to improving landslide susceptibility modeling within road buffer zones but also to providing valuable insights for landslide mitigation throughout the entire lifecycle of road construction projects.

Although machine learning and deep learning techniques have been widely applied in landslide susceptibility modeling [13], the influence of non-landslide sampling strategies on model performance and robustness remains underexplored [9,14]. Moreover, comparative studies integrating multiple models under consistent sampling frameworks are limited, which restricts understanding of optimal model selection and interpretability. Addressing these gaps is critical for developing reliable and interpretable susceptibility maps that better support disaster risk management and land-use planning.

2. Study Area

The G3033 Duku Expressway, located in Xinjiang Uyghur Autonomous Region, spans approximately 460 kilometers in length and has a ratio of bridges and tunnels exceeding 50%. It functions as the replacement of the original G217 Duku Highway which only opens from July to October for safety considerations due to heavy snowfall. It enables all-weather traffic and is an important component for enhancing the resilience of Xinjiang’s road network. Serving as a crucial line connecting southern and northern Xinjiang, the G3033 Duku Expressway holds significant importance for the exploitation of mineral resources, tourism development, and economic growth in the Tianshan Mountains area.

The G3033 Duku Expressway crosses the Tianshan Mountains, which was characterized by significant topographical variations, complex landforms, densely distributed faults and strong tectonic movement [15]. The combined effects of the area’s lithology, lineaments, topography, meteorology, and hydrology have resulted in frequent occurrence of high-speed and long-run-out landslides ranging from 200m-2km [16,17]. Thus, the 2km buffer zone of the route was chosen as our AOI (Figure 1). It is situated within the following geographic coordinates: 44.5465° N for the northern boundary, 82.8549° W for the western boundary, 85.1121° W for the eastern boundary, and 41.5999° N for the southern boundary, and covers area of 2215.00 km². Figure 2 shows the significant elevation changes along the Duku Expressway, showing that it traverses through diverse terrains with varying altitudes. Specifically, the elevation of the AOI ranges from 358 meters (green) to 4678 meters (red).

The annual average temperature is 9.6°C, with a typical January average of -15°C and an extreme minimum of -35°C. In July, the average temperature reaches 18°C, with an extreme maximum of 32°C. This area is classified as semi-arid, with precipitation being low in all seasons except summer. Rainfall primarily occurs as heavy rain or showers, with the regions along the route experiencing heavy rain almost every year. The heavy rainfall typically concentrates between June and August, accounting for over 50% of the total annual precipitation. The main lithologies are limestone, sandstone, conglomerates, and sand-gravel.

The surrounding population density is quite low, and the availability of resources for emergency road recovery is limited. Geological hazards along the route exhibit characteristics such as frequent occurrence, recurrence, and hidden nature. Once geological disasters occur, ensuring timely rescue operations becomes challenging. To achieve the overall goal of a safe and all-weather accessible expressway, it is imperative to conduct a landslide susceptibility assessment along the G3033 Duku Expressway, establish an accurate and efficient landslide susceptibility model, and provide quantitative data for strengthening disaster prevention and mitigation capabilities.

3. Material and Method

3.1. Landslide Inventory

Building a landslide inventory is a crucial preliminary step in susceptibility evaluation [18,19]. This process involves identifying the existing landslides, including information of their location, scale, type, and other features. A landslide inventory database can be developed through various methods such as field surveys, existing landslide distribution maps, satellite image interpretation, Synthetic Aperture Radar (SAR) image shift tracking, or a combination of these approaches [20]. Paul, et al. combined field geologic investigations, remote sensing impact interpretation, and historical landslide reports to locate 816 historical landslides in the study area [21]. Faming, et al. cross-referenced the local Department of Natural Resources (DNR) historical field survey data with Landsat 8 Land Imager (OLI) remote sensing image interpretation results in order to build the landslide inventory database [7]. In this study, 101 landslide events were identified by field survey and Google Earth satellite image interpretation. Figure 3 shows two landslides identified during the field survey. In this study, the centroid of each landslide deposit was used to construct the landslide inventory dataset.

Figure 4 presents the workflow for landslide susceptibility assessment within the 2 km buffer zone alongside the G3033 Duku Expressway.

3.2. Contributing Factors

Landslide susceptibility evaluation indicators can be categorized into two primary types: internal and external factors [22]. Internal indicators refer to topographic and geomorphological data of the AOI, such as elevation, slope, and aspect, etc. External factors, often referred to as triggering factors, include variables like rainfall and seismic activity. Depending on the characteristics of the dataset, these indicators can further be classified into continuous indicators (elevation, slope, etc.) and nominal indicators (stratigraphic lithology, land use type, etc.).

Currently, there is no unified standard or specification for constructing landslide susceptibility evaluation indicator systems [19,23,24]. The selected evaluation indicators vary significantly across different regions. Table 1 compiles landslide susceptibility evaluation indicators from various studies, showing substantial differences in the selection and number of indicators. Most studies’ evaluation systems consist solely of internal indicators (topographic and geomorphological indicators), with commonly used indicators being elevation, slope, aspect, and plan and profile curvature. These indicators can be computed from the Digital Elevation Model (DEM) of the study area using Geographic Information Systems (GIS) tools. The distance to rivers, distance to roads, and distance to faults is also frequently considered. Some studies use the distance to linear structures as an indicator. The concept of linear structures refers to significant landscape lines formed by joints, fractures, shear zones, and faults, which can more comprehensively reveal the topographic and geomorphological features of the study area [25].

The selection of external (triggering) indicators was based on prominent disaster-inducing factors specific to the local area. Considering the meteorological and topographical characteristics as well as landslide mechanisms observed in the Tinshan region and supported by a comprehensive literature review, we initially selected 16 indicators: plain curvature, profile curvature, roughness, STI, TWI, elevation (DEM), aspect, slope, NDVI, land cover, distance to roads, distance to rivers, distance to faults, lithology, and precipitation (including rainfall and snowfall) as well as snowmelt. Precipitation and snowmelt were chosen as triggering indicators since they are recognized as primary landslide-inducing factors during the warm seasons in Xinjiang [30,31]. The selection is also consistent with frequently used factors identified in previous studies, including a review by Md. Sharafat et al. [5], which summarized 45 landslide susceptibility indicators commonly used between 2008 and 2020, highlighting slope, land use/cover, aspect, elevation, geology, NDVI, distance to roads, rainfall, and distance to rivers as the most frequent.

The digital elevation model (DEM) was acquired through the ASTER Global Digital Elevation Model V003 from https://search.earthdata.nasa.gov/ with a resolution of 30 meters. Plain curvature, profile curvature, roughness, STI, TWI, elevation (DEM), aspect, and slope were calculated and extracted in ArcMap 10.8 based on the DEM of our area of interest (AOI). The road and river datasets were obtained through the OpenStreetMap program (https://www.openstreetmap.org). The faults dataset was sourced from the results published in July 2024 by the research team led by Xu et al. at the Institute of Geology [32], China Earthquake Administration. The multiple ring buffer function was employed to construct the distances to roads, rivers, and faults. The normalized difference vegetation index (NDVI) was calculated based on Landsat 8 satellite imagery (https://www.usgs.gov/). Land cover and lithology data were acquired from the global SOTER database, developed by the Food and Agriculture Organization (FAO), the United Nations Environment Programme (UNEP), and ISRIC. Snowmelt data were obtained from the monthly-scale snowmelt dataset (1951-2020) provided by Yong, et al. [33]. Additionally, the precipitation dataset was sourced from ERA5 monthly averaged data on pressure levels from 1940 to the present, available at ECMWF (https://www.ecmwf.int), the point rainfall data were processed using the Kriging method [34]. All datasets were converted to 30-meter resolution raster format in ArcMap 10.8 to facilitate future calculations (Figure 5).

3.3. Landslide Susceptibility Assessment

3.3.1. Sampling the training & validation dataset

This study used the common 1:1 landslide/non-landslide sampling strategy [35], while the non-landslide points were sampled using two different methods (Figure 6): Group 1: sample randomly within the AOI. Group 2: sample within the interpolated section of the AOI and the Low and Medium hazard level area of The Global Landslide Hazard Map (GLHM) to avoid accidentally including points located in actual landslide-prone area into the non-landslide dataset. The application of Group 2 (GLHM) has resulted in a notable reduction in the concentration of elevation within the dataset. Specifically, the median elevation for Group 2 is 2619 m, compared to 2756 m for the Group 1 (Random) (Figure 7).

Global landslide hazard map is a quantitative landslide hazard map for the whole world produced by Arup under the joint request of the World Bank and the Global Facility for Disaster Reduction and Recovery (GFDRR). The dataset comprises gridded maps of estimated annual frequency of significant landslides per square kilometer, categorized as follows: Very low: <0.0001 per square kilometer; Low: 0.0001 - 0.001; Medium: 0.001-0.01; High: >0.01. This project was conducted based on global-level datasets to show spatial and temporal variation of landslide occurrence internationally, however lacking accuracy on regional level. Though only three events in Xinjiang are listed in the global landslide hazard map inventory, to the extent that the map reflects the homogeneous features of landslides reported worldwide.

The values of contributing factors for all landslide and non-landslide points were extracted from their respective raster grids to form the dataset. The dataset—comprising 101 landslide points (assigned a value of 1) and 101 non-landslide points (assigned 0)—was then randomly split into training (70%) and validation (30%) groups for the modeling procedure. Prior to training, the input features were normalized using L2 normalization to ensure consistent scaling across variables, thereby reducing the risk of overfitting and enhancing the stability and efficiency of the learning process.

3.3.2. Modeling techniques

To address the critical concern of overfitting, especially given our limited dataset of 204 instances (102 landslide, 102 non-landslide) and 12 geological indicators, we implemented tailored regularization strategies for each model: Naive Bayes (NB), Support Vector Machine (SVM), Random Forest (RF), and XGBoost. We employed 5-fold stratified cross-validation to conduct both hyperparameter tuning and evaluate overall model performance. A brief overview of each technique is provided below.

• Ÿ Naïve Bayes (NB)

Naïve Bayes is a probabilistic classifier based on Bayes' theorem with the assumption of feature independence. Given a feature vector $\mathbf{x}=(x_1,x_2,\dots ,x_n)$ and class label $C_k$, the posterior probability is calculated as [36]:

$$P\left(C_k∣x\right)={\displaystyle \frac{P\left(C_k\right){\prod }_{i=1}^{n}P\left(x_i∣C_k\right)}{P\left(x\right)}}$$

(1)

Where $P\left(C_k\right)$ is the prior probability of class $C_k$, and $P\left(x_i∣C_k\right)$ is the likelihood of feature $x_i$, given class $C_k$. The classifier assigns the class with the highest posterior probability.

For Naive Bayes, traditional regularization techniques like L2 regularization or Dropout were not applicable, as it is a probabilistic model without iterative training or adjustable weights. Additive smoothing was used to prevent zero probabilities for unseen feature combinations, which was crucial for robust performance on small datasets. We tuned the alpha parameter for MultinomialNB or BernoulliNB using cross-validation. The optimal alpha value was 1.0.

• Support Vector Machine (SVM)

Support Vector Machine is a supervised learning algorithm that seeks to find the optimal hyperplane that maximizes the margin between two classes [37]. Given a training dataset ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^n$, where $x_i\in R^d$ and $y_i\in \{-\mathrm{1,1}\}$, SVM solves the following optimization problem:

$$\underset{w,b}{\min }{\displaystyle \frac{1}{2}}|\mathit{w}{|}^2\hspace{1em}\mathrm{subject}\mathrm{to}\hspace{1em}{y}_i\left(w^{\mathrm{T}}{x}_i+b\right)\ge 1$$

(2)

SVM inherently regularize by maximizing the classification margin. This regularization was controlled by the C parameter and, when using the RBF kernel, the gamma parameter. A smaller C value led to stronger regularization by allowing more misclassifications to achieve a larger margin, beneficial for our small dataset. Similarly, a smaller gamma with the RBF kernel resulted in a simpler, more generalized decision boundary. We used a logarithmic grid search for both parameters. The optimal C value was 0.1, and the optimal gamma (RBF kernel) was 0.001. Dropout and Early Stopping were not applicable to SVMs.

• SVM–RF Hybrid Modeling Framework (SVM-RF)

The SVM-RF model is a hybrid approach that combines the strengths of Support Vector Machines (SVM) and Random Forests (RF). While SVM excels at finding the optimal decision boundary in the feature space, Random Forests provide robust feature selection and reduce overfitting by aggregating predictions from multiple decision trees [38]. The integration of these methods has been shown to enhance classification accuracy in landslide susceptibility studies by leveraging the strengths of both techniques.

Each model was individually optimized using distinct strategies. For SVM, we employed the same optimization method as with our standalone SVM model, tuning its box constraint (C) and kernel scale (gamma) within the ranges [0.1, 100] and [0.001, 10], respectively. For random forest, it naturally resisted overfitting through their ensemble nature, employing bootstrap aggregating and random feature selection during tree construction. Our regularization strategy constrained the complexity of individual trees and enhanced ensemble diversity. Key parameters that were tuned included max_depth, min_samples_split, min_samples_leaf, and max_features. The optimal parameters found were max_depth=7, min_samples_split=5, min_samples_leaf=3, and max_features='sqrt'. For our small dataset, we favored smaller max_depth values and larger min_samples_split/min_samples_leaf to build more generalized trees.

After tuning, models were evaluated on a validation set. If one model significantly outperformed the other (difference in accuracy > 1%), the better model was chosen. If performance was similar, a weighted ensemble was used, where the final prediction was a weighted average of SVM and RF outputs, with weights inversely proportional to validation error.

• Extreme Gradient Boosting (XGBoost)

XGBoost is a tree-based ensemble learning method that builds models sequentially, where each new tree corrects the residuals of the previous ones. It minimizes a regularized objective function [39]:

$$\mathcal{L}=\sum _{i=1}^{n}l(y_i,{\hat{y}}_i^{\left(t\right)})+\sum _{k=1}^t\mathsf{\Omega }\left(f_k\right)$$

(3)

Where $l$ is a differentiable loss function (e.g., log loss), and $\mathsf{\Omega }$ is a regularization term that penalizes model complexity.

We employed both L1 (reg_alpha) and L2 (reg_lambda) regularization to penalize the weights of the leaf nodes, along with Early Stopping to prevent excessive boosting iterations. The optimal reg_lambda (L2) was 1.0, and the optimal reg_alpha (L1) was 0.01. We used Early Stopping by monitoring performance on a validation set; training halted if validation performance did not improve for 20 early_stopping_rounds. Other crucial regularization parameters for XGBoost were also tuned, with optimal values including max_depth=4, min_child_weight=1, subsample=0.8, and colsample_bytree=0.7.

3.3.3. Assessing predictive performance

The Receiver Operating Characteristic (ROC) curve is a graphical representation used to evaluate the performance of a binary classification model. At each step, the True Positive Rate (TPR) and the False Positive Rate (FPR) were calculated, and they were plotted as the vertical and horizontal coordinates, respectively, to obtain the “ROC curve.” It uses the True Positive Rate (TPR) on the vertical axis and the False Positive Rate (FPR) on the horizontal axis. Other performance metrics derived from the confusion matrix are presented in Equations 8-11.

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(8)

$$Recall={\displaystyle \frac{TP}{TP+FN}}$$

(9)

$$Accuracy={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(10)

$$F_1=2·{\displaystyle \frac{Precision·Recall}{Precision+Recall}}$$

(11)

Five-fold cross-validation was used to evaluate model performance. The dataset was partitioned into five equal or nearly equal subsets, with each subset used once as the test set while the remaining subsets were used for training [40]. To ensure representative fold composition, K-means clustering was applied to group data into distinct clusters based on feature similarity prior to partitioning [41,42]. K-means-based 5-fold cross-validation is well-suited for landslide susceptibility assessment because it accounts for the inherent spatial heterogeneity of landslide occurrences, for validating the model on spatially disjoint clusters tests its ability to generalize across different terrain characteristics, reflecting real-world scenarios where predictions are applied to new, unseen areas [12].

Model performance was evaluated using two metrics: area under the receiver operating characteristic curve (AUROC) and accuracy. AUROC measures the ability of a classifier to distinguish between classes, with a value of 1 indicating perfect classification and 0.5 indicating random guessing. Accuracy represents the proportion of correctly classified instances among all samples and was chosen as a complementary metric to provide an intuitive measure of overall correctness. The calculation of accuracy is given in Equation 10.

The cross-validation procedure was repeated 20 times to obtain stable performance estimates. The resulting data were used for subsequent statistical significance analyses. Statistical significance of performance differences was assessed to determine if observed variations were attributable to factors beyond random variability. The Mann–Whitney U test was used to compare performance between different sampling methods, while the Kruskal–Wallis one-way ANOVA test was applied to compare multiple susceptibility models. To control for false discoveries arising from multiple comparisons, the Benjamini–Hochberg procedure was used to control the false discovery rate (FDR).

This evaluation framework is critical not only for quantifying model performance but also for ensuring that model improvements are statistically meaningful. It supports the reliability and reproducibility of the susceptibility assessment, making the findings applicable to real-world risk management and decision-making scenarios.

The modeling and performance metric calculations were conducted using MATLAB R2023a. All statistical significance analyses were performed using IBM SPSS Statistics 27.

3.3.4. Assessing factors’ importance in models

To quantify the contribution of each contributing factor to the predictive performance of the employed machine learning models, this study adopts the SHAP (SHapley Additive exPlanations) method [43]. SHAP is a game-theoretic approach rooted in cooperative game theory, which assigns an importance value to each feature by calculating its marginal contribution across all possible combinations of features.

For a given model, the SHAP value (${\varphi }_i$) of a feature $i$ represents the average change in the model’s output when i is included in a prediction, weighted across all possible feature subsets. This is formalized as:

$${\varphi }_i=\sum _{S\subseteq F\setminus \left\{i\right\}}{\displaystyle \frac{\left|S\right|!\left(\right|F|-|S|-1)!}{\left|F\right|!}}\left(f\right(S\cup \left\{i\right\})-f(S\left)\right)$$

(12)

where $F$ is the set of all features, $S$ is a subset of features excluding $i$, and $f\left(S\right)$ is the model’s prediction for subset $S$.

SHAP provides locally interpretable explanations of these landslide susceptibility prediction models, thereby addressing the limitation raised by a few studies - the poor interpretability of most machine learning approaches [6,12,49]. The SHAP analysis was performed using the shapley function available in the Statistics and Machine Learning Toolbox of MATLAB R2023a.

4. Results and Discussion

4.1. Multicollinearity assessment of the contributing indicators

To minimize multicollinearity among input variables, we conducted Pearson correlation analysis between all pairs of the initial 16 contributing factors. Indicators with an absolute correlation coefficient exceeding 0.5 were considered highly correlated. The Pearson correlation coefficients between each pair of the 16 contributing indicators are shown in Figure 8. As a result, plain curvature, roughness, rainfall, and snowmelt were excluded from the final modeling process. Notably, rainfall and snowmelt are the only triggering factors in our dataset.

The correlation coefficient between elevation and rainfall reached 0.618, while that between elevation and snowmelt was 0.725, indicating strong linear relationships. These findings are consistent with previous research highlighting the influence of topography on precipitation patterns in mid- and high-latitude mountainous regions [45]. For example, Bingqi et al. [46] reported a significant linear relationship between precipitation and elevation in Northern Xinjiang. Similarly, Zheng et al. [47] identified a positive elevation–precipitation relationship during the warm season in Xinjiang. In our study area, warm-season precipitation contributes substantially to the annual total, which further justifies the strong observed correlation.

Therefore, despite rainfall and snowmelt being triggering indicators, their high correlation with elevation suggests that excluding them will not significantly reduce the model’s ability to capture triggering conditions. This approach also helps avoid redundancy and improves model stability.

Consequently, the final set of 12 contributing factors included in the model were profile curvature, sediment transport index (STI), topographic wetness index (TWI), elevation, aspect, slope, normalized difference vegetation index (NDVI), land cover, distance to roads, distance to rivers, distance to faults, and lithology.

Given this evidence and consistent with prior studies, excluding the highly correlated indicators does not compromise the comprehensiveness or robustness of the dataset. Similar correlation-based indicator selection procedures have been widely adopted in landslide susceptibility modeling to ensure model stability and interpretability [18,48].

4.2. Landslide susceptibility maps

Maps of landslide susceptibility were produced for visual comparison based on the probability of landslide initiation at each raster point, derived from the tuned aforementioned models that provided predictions along with probabilities. The Jenks natural breaks method [49] was employed to reclassify the probability of the original map into five classes: very low, low, medium, high, and severe, representing the susceptibility of 3,093,091 raster points in total, namely the potential of landslide initiation when triggering events happened (Figure 9).

By comparing the susceptibility maps of each model, it is evident that while the area proportions of the predicted risk levels vary, the distribution of high-risk and severe-risk zones shows a high degree of consistency. These zones are consistently concentrated along the section of the G3033 Duku Expressway that traverses the Tianshan Mountains. Although the models differ in their underlying mechanisms, they all identified similar high-risk areas along this segment, indicating that topographic and geological factors played a dominant role in model determination. The consistency in the distribution of high-risk zones also suggests that, despite differences in algorithms and non-landslide unit sampling strategies, the models exhibit a certain degree of robustness in landslide susceptibility assessment.

Table 2 shows the zonal statistics of landslide susceptibility maps, providing a detailed breakdown of the spatial distribution of landslide susceptibility areas across different models and sampling strategies. The spatial distribution, number of landslide points, and landslide density within each zone are used to assess the effectiveness of each model.

Among all models, the SVM-GLHM and SVM-RF-GLHM models exhibit the highest landslide point densities in the “Severe” zones, reaching 0.1679 point/km² and 0.1791 point/km², respectively. This indicates a strong clustering of actual landslide occurrences in areas predicted as high risk, demonstrating strong model performance under the GLHM strategy.

Conversely, models utilizing random sampling, such as SVM-Random and SVM-RF-Random, also captured a substantial proportion of landslides in their "Severe" zones, achieving point densities of 0.2010 point/km² and 0.1525 point/km², respectively. While these models showed a more dispersed distribution of landslide points across zones, their high-risk categories still effectively identified landslide locations.

The XGBoost-GLHM model achieved the highest landslide point density (0.2110 point/km²) in the “Severe” zone, while predicting 0 points in the “Very Low” zone, suggesting effective discrimination.

In contrast, models like NB-GLHM and NB-Random present lower densities across all classes, with 0.0980 point/km² and 0.2139 point/km² in the “Severe” zone, respectively, indicating relatively weaker spatial prediction power compared to other models.

Overall, models trained with the GLHM sampling strategy tend to result in a more concentrated landslide distribution in the “High” and “Severe” zones, effectively delineating high-risk areas. This highlights the advantages of the GLHM sampling method in capturing spatial heterogeneity and improving susceptibility map accuracy.

4.3. The predictive performance

Figure 10 presents the ROC curves of all models under both GLHM and Random sampling strategies, along with their corresponding Area Under the Curve (AUC) values. All models showed acceptable to excellent classification ability, with AUC values above 0.80, indicating their effectiveness in distinguishing landslide-prone areas.

Among the GLHM-based models, XGBoost-GLHM achieved the highest AUC value (0.9861), followed by SVM-RF-GLHM (0.9683), SVM-GLHM (0.9107), and NB-GLHM (0.8963). The ensemble learning models (XGBoost, SVM-RF) clearly outperformed the traditional models (e.g., Naive Bayes and SVM), highlighting their superior capability in capturing complex feature interactions.

Under the Random sampling strategy, XGBoost-Random still achieved the highest AUC (0.9746), followed by SVM-RF-Random (0.9362), SVM-Random (0.8440), and NB-Random (0.8190). While the performance of simpler models such as NB and SVM declined under Random sampling, ensemble models like XGBoost and SVM-RF remained robust and continued to deliver high classification accuracy.

Among all the evaluated models, ensemble learning approaches—particularly XGBoost and SVM-RF—consistently achieved the highest AUC values across both sampling strategies, highlighting their strong overall performance. These models’ ability to capture nonlinear relationships and effectively integrate multiple predictive features likely contributed to their superior classification accuracy. Notably, XGBoost maintained top performance under both GLHM (AUC = 0.9861) and Random (AUC = 0.9746) sampling, while SVM-RF also delivered high AUC scores (0.9683 and 0.9362, respectively), underscoring their robustness even in the absence of spatially structured data. In contrast, traditional models such as Naive Bayes and SVM exhibited more pronounced performance declines under Random sampling, suggesting greater sensitivity to sampling variability. Overall, the consistent and high performance of XGBoost and SVM-RF indicates their strong generalization capacity and practical utility for landslide susceptibility mapping in diverse scenarios.

To further explore differences between sampling methods and models, we conducted statistical significance tests. The data used for these tests were obtained from 20 repetitions of 5-fold spatial cross-validation. Figure 11 and Figure 12 show boxplots of AUROC and accuracy for models under different sampling methods. Table 3 presents Mann–Whitney U test results comparing performance across sampling strategies, while Table 4 shows Kruskal–Wallis test results for model comparisons. These tests evaluate whether the differences observed in model performance are statistically significant.

The comparison between sampling methods reveals significant differences in model performance across all tested algorithms. Models trained with GLHM-informed sampling consistently outperformed those using random sampling within the AOI in both AUROC and accuracy metrics.

For AUROC, the NB model showed a notable improvement, increasing from a median of 82.05% (IQR 4.29) with random sampling to 90.49% (IQR 3.44) under GLHM sampling, representing a significant gain of +8.44% (p < 0.001). Similarly, the SVM model improved by +7.11% (p < 0.001), moving from 83.71% (IQR 5.75) to 90.82% (IQR 2.93). The SVM-RF model also saw an improvement of +3.45% (p = 0.023), and XGBoost improved by +3.04% (p = 0.029). Accuracy gains followed the same trend. The NB model increased by +11.30%, from 71.31% (IQR 2.68) to 82.62% (IQR 0.76), which was highly significant (p < 0.001). The SVM model showed an +8.33% increase (from 74.51% to 82.84%, p < 0.001), while SVM-RF improved by +7.40% (p = 0.002), and XGBoost increased by +8.31% (p < 0.001).

These results indicate that the GLHM-based sampling strategy leads to statistically significant enhancements in model predictive performance compared to random sampling. This suggests that sampling guided by hazard-level information improves the quality of the training data, enabling models to better distinguish between landslide and non-landslide areas.

Table 4 presents the results of the pairwise statistical comparisons among the evaluated models for both AUROC and Accuracy metrics, under different sampling strategies. These comparisons rigorously assessed the statistical significance of observed performance differences. The null hypothesis for each test stated that the distributions of the compared model pairs were identical. All significance values were adjusted using the Bonferroni correction to account for multiple comparisons, with a significance level of 0.050.

Under the GLHM sampling strategy, the pairwise comparisons for AUROC revealed a statistically significant difference only between NB and XGBoost (Adj. Significance = 0.044)*, indicating that XGBoost significantly outperformed Naive Bayes in terms of AUROC. No other model pairs (NB-SVM, NB-SVM-RF, SVM-SVM-RF, SVM-XGBoost, SVM-RF-XGBoost) showed statistically significant differences in AUROC, suggesting their performance distributions were comparable under this sampling approach.

For Accuracy under GLHM sampling, a consistent pattern of statistically significant differences emerged. SVM-RF notably outperformed NB (Adj. Significance = 0.044), and similarly, SVM-RF also showed a significant difference when compared to XGBoost (Adj. Significance = 0.044). Furthermore, NB was significantly outperformed by XGBoost (Adj. Significance = 0.044)*. These results indicate clearer distinctions in accuracy performance among models when trained with GLHM sampling.

Under the Random sampling strategy, the pairwise comparisons demonstrated more pronounced statistically significant differences. For AUROC, NB was significantly outperformed by SVM-RF (Adj. Significance = 0.044)* and even more strongly by XGBoost (Adj. Significance < 0.001*). Similarly, SVM also showed highly significant differences compared to XGBoost (Adj. Significance < 0.001*) and a significant difference compared to SVM-RF (Adj. Significance = 0.044)*. This suggests that XGBoost and SVM-RF, as well as SVM, were generally more distinct from their counterparts in AUROC under random sampling.

For Accuracy under Random sampling, the differences were even more striking. NB was significantly outperformed by SVM (Adj. Significance < 0.001*) and by XGBoost (Adj. Significance = 0.044)*. This highlights the substantial improvement these models offered over Naive Bayes in terms of accuracy when random sampling was employed. Other pairs, such as NB-SVM-RF, SVM-RF-SVM, SVM-RF-XGBoost, and SVM-XGBoost, did not show statistically significant differences in accuracy.

These pairwise comparison results provide crucial statistical insights into the relative performance of the models and the impact of sampling strategies. While some models (e.g., SVM, SVM-RF, XGBoost) often performed numerically better, these tables confirm when such differences were statistically significant after Bonferroni correction. Notably, XGBoost consistently showed statistically significant superiority over Naive Bayes in both AUROC and Accuracy across both sampling strategies, underscoring its robust performance. The GLHM sampling strategy generally clarified differences in accuracy, while Random sampling revealed more profound distinctions in AUROC among model pairs.

Under the Random sampling strategy, the analysis indicated that more complex models, specifically SVM-RF and XGBoost, generally demonstrated significantly better predictive performance and greater robustness compared to simpler models such as NB and SVM. This was evident as NB was significantly outperformed by SVM-RF and XGBoost in AUROC, and by SVM and XGBoost in Accuracy. Similarly, SVM was significantly outperformed by SVM-RF and XGBoost in AUROC. This highlights the stronger capacity of SVM-RF and XGBoost to handle the higher levels of label noise and variability inherent in randomly sampled datasets, which disproportionately impacted the accuracy and stability of NB and SVM. Despite the noisier data, SVM-RF and XGBoost often exhibited statistically comparable predictive power among themselves under random sampling.

Conversely, under the GLHM sampling strategy, the performance landscape shifted. While XGBoost still demonstrated a statistically significant superiority over NB in terms of AUROC, and SVM-RF, XGBoost, and NB all showed significant differences in Accuracy for certain pairs, the overall trend was towards more consistent and robust performances, often reflected by narrower IQRs (Table 3). This suggests that hazard-informed sampling effectively reduces noise and label ambiguity, leading to a higher baseline performance across models, even if some statistically significant differences persist.

Based on the results, model selection should consider the quality of the sampling strategy and the associated data noise:

• When using hazard-informed sampling methods such as GLHM, which reduce label noise and improve data quality, multiple models including: SVM-RF, XGBoost, SVM, and NB demonstrate comparable and robust performance. In this scenario, model choice can be guided by other factors such as computational cost, interpretability, and ease of deployment, since performance differences are minimal and models exhibit consistent robustness.
• When data are sampled randomly without hazard-level constraints, resulting in noisier and less reliable labels, more complex models like SVM-RF and XGBoost should be preferred. These models show significantly better predictive performance and greater robustness compared to simpler models, which are more sensitive to label noise and exhibit lower accuracy and stability.

4.4. The importance ranking of contributing factors

Table 5 summarizes feature importance for each model under the two sampling methods, quantified by mean absolute SHAP values. These feature importance patterns are further illustrated by the bi-modal violin plots of SHAP value distributions for each model under both sampling methods (Figure 13).

Under the GLHM sampling, STI exhibited the highest importance in SVM (0.294) and NB (0.065), while NDVI was the most influential feature for SVM-RF (0.320) and XGBoost (0.172). Slope was notably important in NB (0.137). Other features such as TWI and profile curvature showed moderate influence (e.g., TWI: 0.203 for SVM, 0.165 for SVM-RF).

In the Random sampling group, STI remained highly important for SVM (0.392) and SVM-RF (0.345). NDVI was also among the top features for these models (SVM: 0.227; SVM-RF: 0.254). Lithology and elevation gained relatively more importance compared to GLHM, with lithology at 0.090 in NB and elevation up to 0.126 in SVM. Generally, feature importance values were lower and exhibited greater variability under Random sampling, which likely reflects the noisier nature of the data acquired through this method.

The consistent prominence of STI, TWI and NDVI across models and sampling strategies highlights their critical roles in landslide susceptibility. These findings underscore the critical role of hydro-dynamic processes in landslide initiation within the Tianshan region. Specifically, the combined effects of spring snowmelt and intense summer rainfall create dynamic hydrological conditions that are pivotal in triggering landslide events.

This analysis revealed nuances in feature reliance based on model complexity and data quality. For instance, slope consistently emerged as a key factor in simpler models like Naive Bayes, especially under GLHM sampling, affirming the significance of terrain steepness. The increased importance of lithology and elevation under Random sampling suggests that models may depend more heavily on fundamental geological features when confronted with higher levels of label noise. This observed variability in feature importance between the two sampling methods powerfully highlights the direct impact of data quality and noise levels on model interpretability and the features models prioritize for prediction.

4.5. Possible Use

This study enables rapid identification of areas highly prone to landslides along the route. By analyzing slope stability and integrating key hazard indicators—such as displacement, rainfall, and crack development—into a monitoring system, it efficiently locates typical high-risk slopes and supports road safety during operation. The approach emphasizes efficiency in early detection and risk mitigation, which can also further support prioritization of slope reinforcement and drainage upgrades along the route.

It also aids in strategically positioning emergency supply depots (e.g., rescue materials) by concentrating resources in high-risk road sections characterized by steep terrain, significant elevation fluctuations, and intense snowmelt/rainfall during warm seasons—are prioritized due to their heightened vulnerability.

5. Conclusions

Landslides remain the most significant challenge in disaster prevention and mitigation, particularly along highway buffer zones, where multiple fatal accidents last year have drawn increased attention. In this study, NB, SVM, SVM-RF, and XGBoost models, along with two different non-landslide sampling strategies (GLHM and Random), were employed to assess landslide susceptibility in the 2 km buffer zone of the G3033 Duku Expressway. 16 contributing factors were preliminary selected. To avoid multicollinearity, we excluded plain curvature, roughness, rainfall, and snowmelt due to their high correlation (|r| > 0.5) with other indicators, a pattern consistent with known topographic-precipitation relationships in alpine regions.

The key findings are summarized below:

• XGBoost consistently achieved the best performance under both sampling strategies. Under GLHM sampling, XGBoost obtained the highest AUROC (94.61%, IQR: 3.58) and accuracy (84.30%, IQR: 2.81), outperforming the other models. Pairwise comparisons further confirmed the statistical superiority of XGBoost, with significant differences observed especially when compared to NB and SVM models (adjusted p-values < 0.05 or < 0.001 in multiple comparisons). These results highlight the strong generalization capability of XGBoost in complex landslide susceptibility assessments, benefiting from its ability to model nonlinear interactions and handle high-dimensional feature spaces efficiently.
• GLHM sampling method demonstrated a clear advantage over random sampling across all models. Both AUROC and accuracy were consistently improved under GLHM. For AUROC, the improvements for NB, SVM, SVM-RF, and XGBoost under GLHM reached +8.44 (p<0.001), +7.11 (p<0.001), +3.45 (p=0.023), and +3.04 (p=0.029), respectively. A similar trend was observed for accuracy, with increases of +11.30% (p<0.001) for NB, +8.33% (p<0.001) for SVM, +7.40% (p=0.002) for SVM-RF, and +8.31% (p<0.001) for XGBoost. These results indicate that GLHM effectively enhances model performance, likely by improving the representativeness of the training data and better capturing the underlying distribution of landslide and non-landslide units.
• Interpretability analysis using SHAP values further demonstrated that the choice of sampling method not only affects model performance but also the attribution of contributing factors. Under GLHM, top-ranked features were consistent across models, with STI (e.g., 0.2936 in SVM, 0.2947 in SVM-RF), NDVI (e.g., 0.3202 in SVM-RF, 0.2358 in SVM), and slope (e.g., 0.1374 in NB) appearing most frequently among the top three features. In contrast, under random sampling, feature rankings varied more widely, and models exhibited greater reliance on features such as elevation (e.g., 0.1260 in SVM) and lithology, which may reflect artifacts introduced by sampling from spatially mixed hazard contexts.

6. Limitations and Future Work

• Simplified Landslide Representation: This study used the centroids of historical landslide deposits as representative points. While practical, this may not fully capture landslide morphology. Future work could explore alternative point selections, such as the headscarp, or adopt polygon-based landslide datasets to better reflect their spatial extent and improve model fidelity.
• Sampling Strategy Constraints: The GLHM-based method for selecting non-landslide samples helped reduce bias, but its reliance on coarse, global-scale hazard data limits regional applicability. Future research could explore alternative non-landslide sampling strategies to improve the spatial representativeness and robustness of the dataset.
• Static Input Variables: The current models use static environmental factors, ignoring temporal triggers like rainfall or land use change. Future studies should integrate dynamic variables and time-series data to improve the predictive capability and adaptability of susceptibility models.