Fuyuan County is located in the eastern part of Yunnan Province, China, between 103°58′and 104°49′E longitude and 25°00′to 25°58′N latitude. It covers a total land area of approximately 3,348 square kilometers and has a population of around 675,200. The region’s climate transitions from a south temperate semi-humid zone in the south to a north subtropical mountain monsoon climate zone in the north. It experiences an average annual temperature of 13.8 °C, annual precipitation of approximately 1,421.1 mm, and approximately 1,773.9 h of sunshine per year. Fuyuan County is one of the most severely affected areas by landslides in Yunnan Province, characterized by their frequency, wide distribution, sudden onset, and destructive nature. According to historical data, approximately 239 landslides have occurred in the study area, with their spatial distribution shown in Figure 1.

Fuyuan County lies at the boundary between the Eastern Yunnan Fold Belt of the Yangzi Paraplatform and the Southeastern Yunnan Fold Belt of the South China Fold System, delineated by the Neocathaysian Fuyuan-Mile Fault Zone. The northwest pertains to the Eastern Yunnan Fold Belt, and the southeast to the Southeastern Yunnan Fold Belt. This fault zone, a regional controlling structure, exhibits compressional-torsional characteristics with significant displacement and extensive length. Dominant tectonic systems include Shan-shaped, lotus-shaped, northeast-trending, northwest-trending, and arcuate structures. Fuyuan County displays diverse geomorphological types, including karst, erosional-tectonic, tectonic-erosional, denudational-tectonic, and planation-terrace landforms. Karst landforms, covering approximately one-third of the study area, are distributed from eastern Mohong to Laochang Heike, with intense karstification. Erosional-tectonic landforms (820 km²) occur along the Mohong-Dahe corridor. Tectonic-erosional landforms (740 km²), located in areas like Zhuyuan-Yingshang, are heavily incised by rivers and tectonic activity. Denudational-tectonic landforms (18 km²) are found in Zhong’an-Housuo. Planation-terrace landforms (85 km²), distributed along Xiaoyangchang-Yidehe, exhibit gentle topography.

The landslide disaster inventory is a critical component of Landslide Susceptibility Mapping, providing essential data support for model training and prediction. Through field surveys and historical data analysis, it was found that landslides in Fuyuan County predominantly occur in rock masses, Quaternary soil layers, and residual slope deposits, controlled by weak surfaces or joint fractures, exhibiting characteristics of “numerous points, wide coverage, and uneven distribution.” As shown in Figure 2, among the 12 townships in the county, Zhongan Township, Shengjing Township, and Gudan Township have fewer disaster points, while the other townships show more concentrated distributions. Landslides mostly occurred between 2000 and 2023, with June to September being the peak period. Approximately 50% of landslides are located in or near fault zones or stratigraphic contact zones, influenced by tectonic activities, resulting in poor rock mass integrity, fragmented rocks, well-developed joints, and thick weathered soil layers, which provide favorable conditions for landslide formation. Fuyuan County is rich in coal resources, with about one-quarter of its land area being coal-rich zones. A total of 186 coal mine shafts with a history of mining or still in production were documented, and mining activities have further exacerbated landslide risks.

This study focuses on investigating the geological conditions of historical landslides and potential geological hazard points, combined with GF-1 (Gaofen-1) remote sensing imagery for county-wide geological disaster and environmental geology interpretation. In key areas, higher-precision QuickBird data was used to supplement the analysis, providing detailed baseline data for LSM. Ultimately, by integrating official records from the Yunnan Geological Survey Bureau (https://ynddj.org.cn/), optical remote sensing analysis, and field verification, a total of 239 historical landslide events across the 12 townships of Fuyuan County from 2000 to 2023 were compiled. Each event includes precise geographic coordinates, landslide scale, and movement type. Typical landslide site photos are shown in Figure 3.

To effectively address the issues of class imbalance and spatial autocorrelation, this study employed a constrained random sampling method based on a GIS platform to systematically construct a balanced training dataset. Following established practices in landslide susceptibility mapping (Zhao et al., 2021; Guo et al., 2023), 239 non-landslide points were selected from geologically stable units in the ArcGIS environment to match the quantity of landslide points. The specific sampling constraints included: 1) spatial distribution uniformity constraint, ensuring sufficient dispersion of sample points in geographical space to reduce spatial autocorrelation effects; 2) geological environment matching constraint, maintaining comparability between non-landslide and landslide points in the distribution of key environmental factors (e.g., elevation, slope, lithology); and 3) stability verification constraint, confirming through remote sensing time-series monitoring and verification with geological reports that all selected areas have remained stable with no landslide records since official documentation began. This systematic sampling strategy not only effectively balanced the class distribution but also significantly reduced spatial dependency interference during model training, establishing a solid data foundation for developing a robust landslide susceptibility model. The spatial distribution of the ultimately adopted landslide and non-landslide points is shown in Figure 4.

Modern machine learning algorithms (e.g., decision trees, neural networks, kernel methods) inherently capture nonlinear feature interactions, with L1/L2 regularization and Dropout effectively suppressing linearly redundant feature adverse effects. Since Pearson correlation coefficients only detect linear associations, this study uses VIF detection instead. In applications like landslide susceptibility mapping, VIF is a critical diagnostic tool for model robustness and interpretability (Guoliang et al., 2021). In our modeling, this approach ensured parsimonious geomorphological variable selection while maintaining computational efficiency when integrating multi-source topographic, geological, and hydrological data. The VIF <5 threshold was adopted based on empirical machine learning studies and geohazard modeling domain considerations.

Figure 7 shows the model features’ multicollinearity analysis results. VIF analysis indicates all features have values below the critical threshold of 5 (significantly lower than 10), meaning no severe multicollinearity exists. This ensures subsequent model training stability and reliability by avoiding redundant feature interactions, making the feature set suitable for input into the machine learning framework.

The FR results provide an intuitive and physically meaningful explanation of the landslide occurrence mechanisms in the study area. This establishes a geoscientifically grounded expectation and reference point for subsequent machine learning model predictions. Before conducting model training, we already knew which factors were considered important in traditional analysis. As a pre-modeling analysis, FR provides prior knowledge for feature selection and plausibility assessment of model results.

The FR method was employed to evaluate the statistical relationship between landslide occurrence and contributing factors (Lee and Pradhan, 2007). This technique computes the ratio of landslide density within each factor class (e.g., slope intervals or lithology types) to the overall landslide density across the study area. FR values >1 indicate elevated landslide propensity, while values <1 suggest reduced susceptibility. The method effectively processes both continuous variables (e.g., slope angle) and categorical variables (e.g., land use), making it particularly suitable for preliminary susceptibility assessment in our study region. Prior to FR calculation, all factor classes were discretized using natural breaks classification (Jenks optimization). The derived ratios were subsequently aggregated to generate a composite susceptibility index. This robust analytical approach has been extensively validated in geohazard research.

This study developed two ensemble learning methods (utilizing XGboost and LightGBM) for landslide susceptibility prediction, while also selecting several traditionally used evaluation models for comparison. The selected models encompass diverse paradigms to capture different data patterns. The characteristics of these models are presented in Table 2.

XGBoost is an efficient ensemble algorithm that models complex nonlinear relationships, handles missing data, and minimizes prediction errors through boosted trees in landslide susceptibility assessment (Matougui et al., 2025), as formalized in Equation 2: Obj = ∑ i = 1 n l y i , y ^ i + ∑ k = 1 K Ω f k (2) where l is the loss function (e.g., logistic loss), y i and y ^ i are observed and predicted values, f k denotes the kth regression tree, and Ω f k is a regularization term defined by: Ω f = γ T + 1 2 λ ∑ j = 1 T w j 2 (3) where T is the number of leaf nodes, w j is the score on leaf j, and γ , λ are regularization parameters, as given in Equation 3. Compared to traditional methods, XGBoost offers higher accuracy, reduces overfitting, and provides insights into feature importance, making it a valuable tool for landslide risk prediction.

The study used Randomized Search instead of Grid Search for hyperparameter optimization, as Randomized Search is more efficient in high-dimensional parameter spaces and can explore a broader range of parameter combinations with limited computational resources. We adopted a 5-fold stratified cross-validation approach to ensure the ratio of positive to negative samples in each fold is consistent with that of the original dataset. Through this randomized search hyperparameter optimization, the optimal configuration was determined with the following parameters: colsample_bytree = 0.721, gamma = 0.041, learning_rate = 0.0045, max_depth = 10, n_estimators = 552, reg_alpha = 0.065, reg_lambda = 7.866, and subsample = 0.783. This parameter combination effectively controls model complexity while maintaining strong predictive capability for landslide susceptibility assessment.

This study employs LightGBM, a fast and accurate GBDT-based algorithm, to model landslide susceptibility using standardized environmental factors as input features for large-scale, high-dimensional, and imbalanced data (Danlu et al., 2023), as formulated in Equation 4: L = ∑ i = 1 n l y i , y ^ i t − 1 + f t x i + Ω f t (4) where l is the loss function, y i and y ^ i are the true and predicted values, f t is the new tree, and Ω is the regularization term.

LightGBM uses a leaf-wise tree growth strategy that prioritizes splits with the highest loss reduction, with the gain of a split computed as Equation 5: G a i n = 1 2 G L 2 H L + λ + G R 2 H R + λ − G L + G R 2 H L + H R + λ − γ (5) where G and H are the sum of first- and second-order gradients, λ and γ are regularization parameters. Its histogram-based decision tree construction improves training speed and memory efficiency, while embedded feature selection enhances interpretability and generalization.

A similar strategy was adopted for the LightGBM model. The LightGBM model was optimized using histogram-based algorithms for efficient training. The hyperparameter optimization yielded the following optimal configuration: colsample_bytree = 0.704, learning_rate = 0.0098, max_depth = 7, min_child_samples = 49, n_estimators = 274, num_leaves = 48, reg_alpha = 0.013, reg_lambda = 0.179, and subsample = 0.987. This parameter set ensures model robustness while preventing overfitting, making it suitable for landslide susceptibility prediction applications.

The ANN, originally for handwritten digit recognition, is a specialized neural network architecture particularly suited for processing grid-structured data (Haykin, 1994). In this study, we directly employ the ANN model as a classifier for landslide susceptibility prediction. The implementation involves three main steps: 1) transforming raw landslide data into three-dimensional data structures; 2) constructing corresponding ANN architectures tailored to different data dimensions; and 3) conducting susceptibility analysis based on these dimension-specific data and ANN structures. For clarity, we illustrate the analysis process using a 1D-ANN architecture (Figure 8). The 1D-ANN takes 15 × 1 factor vectors as input. Through convolutional and pooling operations in the hidden layers, it extracts high-dimensional landslide features. These features are then mapped to a lower-dimensional space via fully connected layers and finally transformed into the sample label space through nonlinear activation functions, outputting both landslide/non-landslide classifications and their corresponding probability values (Landslide Susceptibility Index, LSI). By applying this 1D-ANN model to predict susceptibility indices for each grid cell, we generate the final landslide susceptibility map.

SHAP (SHapley Additive exPlanations) is an interpretability analysis method grounded in game theory, designed to explain the predictions of any machine learning model (Li and Tian, 2025). Its core principle involves calculating the contribution of each feature to an individual prediction (i.e., the SHAP value), thereby quantitatively revealing how that feature influences the model’s final decision. For instance, in a landslide prediction model, SHAP analysis can explicitly indicate which factors contribute most to a specific high-risk prediction. This approach transforms the “black box” decision-making process of the model into a transparent one, enhancing trust in the model, clarifying the basis of its decisions, and identifying key driving factors. SHAP (SHapley Additive exPlanations) not only explains machine learning model predictions but also effectively conducts feature contribution analysis, quantifying the relative importance of each evaluation factor in landslide susceptibility prediction. Therefore, in landslide susceptibility assessment, there is no need to eliminate any evaluation factors, as SHAP can clearly reveal their contribution to model outputs through calculated SHAP values. This approach retains all factors (e.g., topographic, geological, hydrological, environmental, and anthropogenic) and identifies key drivers while avoiding information loss from subjective elimination, thereby enhancing the scientific rigor and predictive accuracy of the model.

In this study, model performance was assessed through multiple complementary metrics. The Receiver Operating Characteristic (ROC) curve and Area Under the Curve (AUC) served as primary evaluation measures, with AUC values closer to one indicating superior discrimination ability between landslide and non-landslide areas (Fawcett, 2006). For comprehensive assessment, four additional key metrics were employed: Accuracy (%), Precision (%), Recall (%), and F1 score (%). Accuracy quantifies the overall proportion of correct classifications, while Precision reflects the model’s capability to minimize false positives in predicted susceptibility zones. Recall (Sensitivity) evaluates the detection rate of actual landslide locations, particularly crucial for risk mitigation purposes. The F1 score provides a harmonic mean of Precision and Recall, offering a balanced performance measure especially valuable for landslide studies where class imbalance between landslide and stable areas is typically observed. Together, this multi-metric approach (AUC-ROC plus Accuracy/Precision/Recall/F1) ensures robust evaluation of both the model’s predictive power and its operational reliability for landslide risk management applications.

To evaluate the reliability of landslide susceptibility predictions in complex geographical environments, this study employed a bootstrap resampling approach to quantify prediction uncertainty. The bootstrap method generates a large number (500 iterations in this study) of resampled training sets by randomly sampling the original training dataset with replacement (Hesterberg, 2011). For each bootstrap sample, an ensemble learning model was independently trained, resulting in a model ensemble. This process involved two key steps: Model Stability Assessment: The distribution of performance metrics (e.g., AUC values) of this model ensemble on an independent test set was calculated. The central tendency (e.g., mean) of this distribution reflects the model’s average performance, while its dispersion (e.g., standard deviation or confidence interval) directly quantifies the stability uncertainty of the model arising from training data variability.

Figure 9 presents the FR analysis results of 15 landslide-influencing factors, quantifying the correlation between each factor’s categories and landslide occurrence through three key metrics: landslide proportion, area percentage, and FR value. Among all factors, lithology exhibits the most pronounced discriminative power, with relatively soft mudstone (Figure 9g) showing the highest FR value of 4.56, accompanied by a 31.38% landslide proportion despite accounting for only 6.88% of the total area, indicating a strong positive correlation with landslide susceptibility. Proximity to rivers (Figure 9e) emerges as a critical topographic factor, where the <200 m category displays an FR of 2.44, with 37.24% of landslides concentrated in this zone—more than twice its proportional area (15.26%). Mine density (Figure 9j) also demonstrates a significant influence, with the highest density category (0.098–0.149) yielding an FR of 3.48, reflecting enhanced landslide risk in areas with intensive mining activities. In contrast, hard limestone (Figure 9g) and high-elevation zones (>2,260 m, Figure 9a) exhibit FR values of 0.16 and 0.03, respectively, indicating their inhibitory effects on landslide occurrence due to their stable geological properties or limited human activity. For topographic factors, mid-slope gradients (15.39°–22.79°, Figure 9b) and moderate precipitation ranges (13,322–14,447 mm, Figure 9i) show moderate positive correlations (FR = 1.34 and 2.33, respectively), while flat terrain (Aspect: Plane, Figure 9c) and high vegetation coverage (>0.869, Figure 9k) display low FR values (<0.65), consistent with their role in reducing slope instability. These quantitative results provide a robust basis for identifying high-risk zones and understanding the multi-factor synergistic effects on landslide occurrence.

Table 3 presents performance metrics (accuracy, precision, recall, F1 - score) of ML models. LightGBM achieved the highest accuracy (84.10%) and precision (80.25%), with recall at 90.50% and F1 at 83.17%. XGBoost followed, with accuracy of 82.30%, recall reaching 91.67%, and F1 at 83.81%. RF (81.70% accuracy) and KNN (81.25% accuracy) showed competitive results. Conversely, SVM, LR, NB, and ANN exhibited lower performance across metrics. Overall, gradient - boosting (LightGBM, XGBoost) and ensemble (RF) models outperformed traditional (SVM, LR) and ANN approaches, with LightGBM demonstrating the most balanced predictive capability.

Figure 10 illustrates the ROC curves of different models for landslide susceptibility assessment on the test set. The LightGBM achieves the highest AUC of 0.89, outperforming other models, while RF, XGBoost, and KNN show comparable performance with an AUC of 0.87. SVM, ANN, and NB have AUC values of 0.85, 0.85, and 0.83, respectively. These results indicate that ensemble learning models, particularly LightGBM, are more effective for landslide susceptibility assessment.

The bootstrap method is used for quantifying the LightGBM model uncertainty, and the results are shown in Figure 11.

Figure 11a illustrates the distribution of prediction standard deviation (used as a quantitative indicator of uncertainty). The results show that the mean value of prediction uncertainty is 0.141, and the overall distribution exhibits a left-skewed pattern—most samples have low prediction uncertainty (with standard deviations concentrated in the 0.0–0.1 range), while only a small number of samples have uncertainty higher than 0.2. This indicates that the model has high certainty in predicting most samples, with significant prediction fluctuations observed only in a few cases.

Figure 11b presents the distribution characteristics of the 95% confidence interval width, with a mean value of 0.361. Similar to the prediction uncertainty distribution, the confidence interval width also shows a bimodal feature: on the one hand, a large number of samples have narrow confidence intervals (concentrated in the 0.0–0.2 range); on the other hand, some samples have confidence interval widths close to 1.0. This distribution suggests that the model provides narrow credible intervals for the prediction results of most samples, with a wide range of uncertainty only for a small number of samples.

Figure 11c displays the relationship between the average prediction probability and prediction standard deviation (uncertainty) in the form of a scatter plot, where the color of the scatter points represents the true label (0/1) of the samples. An obvious “U-shaped” trend can be observed: when the prediction probability is close to 0 or one (i.e., the model is highly certain about the class prediction), the prediction uncertainty (standard deviation) is low; when the prediction probability is close to 0.5 (i.e., the model is ambiguous about the class judgment), the uncertainty increases significantly. This pattern holds for samples of both labels, indicating that the model’s uncertainty quantification can effectively reflect its “degree of confidence” in the prediction results.

Figure 11d evaluates the model’s stability using the Bootstrap method, with AUC (Area Under the Curve) as the indicator. The results show that the mean value of Bootstrap AUC is 0.964, and the distribution is concentrated in the range of 0.955–0.975. This demonstrates that the model maintains high discriminative ability in multiple resampling validations, featuring good stability and further supporting the reliability of the model’s prediction results.

In summary, the LightGBM model demonstrates overall excellent performance and stability. Its prediction results are reliable and stable in most scenarios, with high uncertainty observed only in a small number of samples, indicating outstanding overall performance.

The landslide susceptibility map generated by the LightGBM model (Figure 12) illustrates the spatial distribution of landslide potential across the study area, which covers 11 townships within Fuyuan county. The map, categorized into five susceptibility levels (Very low, Low, Moderate, High, Very high), reveals a non - uniform distribution of landslide risks.

High and Very high susceptibility zones (depicted in red and dark red) are concentrated in certain townships, such as Yingshang Town and parts of Dahe Town and Zhuyuan Town. These zones coincide with a significant number of recorded landslides (marked as red dots), indicating the model’s ability to identify areas prone to landslide occurrences. In contrast, the southern and some western townships (e.g., Mohong Town, Zhongan town, Shibalian mountain Town) are dominated by Low and Very low susceptibility zones (shown in green), where fewer landslides are observed.

Table 4 presents the landslide susceptibility zoning results generated by the LightGBM model, which classified the study area into five susceptibility levels: Very low, low, Moderate, High, and Very high. The distribution characteristics of landslides across these zones revealed a clear positive correlation between susceptibility levels and actual landslide occurrence. Specifically, the Very high susceptibility zone, accounting for 16.55% of the total area (535.59 km²), contained the largest proportion of landslides, comprising 60.67% (145 landslides) of the total landslide count, with a landslide density of 0.271 number/km². The High susceptibility zone, occupying 13.28% of the area (429.53 km²), hosted 26.36% of the landslides (63 landslides) and had a relatively high density of 0.147 number/km². In contrast, the lower susceptibility zones showed significantly fewer landslides: the Moderate, low, and Very low zones accounted for 7.53% (18 landslides), 5.02% (12 landslides), and 0.42% (1 landslide) of the total landslides, respectively, with corresponding densities of 0.039, 0.020, and 0.001 number/km². These results indicate that the LightGBM model effectively concentrated most landslides in higher susceptibility zones, where the landslide density increased progressively with rising susceptibility levels, demonstrating the model’s ability to capture the spatial association between susceptibility and actual landslide distribution.

Figure 13 shows landslide density rises with higher susceptibility zones in (a). In (b), very high - susceptibility zones have a far larger proportion of landslides than their area percentage, highlighting high landslide - prone characteristics.

The SHAP beeswarm plot for the LightGBM model (Figure 14) visualizes feature impacts on predictions. Each dot represents a feature’s SHAP value in one instance, with the x-axis for SHAP values and color gradient (blue - red) for feature values (low - high).

Mine density shows a wide SHAP value spread, indicating a critical role. Population density and distance to rivers also have notable influences, seen from their SHAP value ranges. Features like FVC, distance to faults, and elevation contribute variably, per their distinct SHAP distributions. Overall, the plot reveals how mine-related and geospatial factors (e.g., mine density, population density) drive the LightGBM model’s susceptibility assessments, clarifying key feature-prediction relationships.

The SHAP dependence plot for mine density in the LightGBM model (Figure 15) illustrates the relationship between standardized mine density, SHAP values, and the interacting feature of distance to rivers (color - coded). As standardized mine density increases, SHAP values generally show an upward trend, indicating a positive influence on the model’s prediction of landslide susceptibility. Notably, at higher mine density values (especially >0.5), SHAP values are predominantly positive and exhibit greater variability, suggesting a stronger and more complex impact on predictions. The color gradient reveals that instances with shorter distances to rivers (warmer colors, e.g., yellow) often coincide with higher SHAP values at given mine densities, implying a synergistic effect: areas with both high mine density and proximity to rivers may face elevated hazard potential, as captured by the LightGBM model. This interaction highlights the importance of considering combined feature effects in interpreting the model’s output.

The SHAP feature importance plot for the LightGBM model (Figure 16) quantifies the contribution of each feature to the model’s predictions using the mean absolute SHAP value. Among the evaluated features, mine density emerges as the most influential, exhibiting the highest mean absolute SHAP value. This indicates that variations in mine density have a substantial impact on predicting landslide susceptibility. Population density and distance to rivers follow, with relatively high mean absolute SHAP values, signifying their considerable roles in shaping the model’s output. In contrast, features such as land use, plane curvature, and TWI show lower mean absolute SHAP values, suggesting a lesser influence on the model’s predictions. Overall, the plot highlights the dominance of anthropogenic and proximity-related factors (mine density, population density, distance to rivers) in driving the LightGBM model’s landslide susceptibility assessments, while terrain and land cover features (e.g., slope, land use) have a more muted impact.

The FR analysis revealed significant patterns in the influence of geological, topographic, and anthropogenic factors on landslide susceptibility. Soft mudstone in lithology exhibited a notably high FR value of 4.56, attributed to its tendency to soften and argillize upon water contact, poor permeability, rapid reduction in rock strength, susceptibility to weathering, low shear strength, and well-developed bedding planes and fractures that facilitate sliding. This aligns with the observation that mudstone areas account for 31.38% of landslides despite occupying only 6.88% of the spatial distribution, consistent with previous findings in Yunnan Province (Liu et al., 2023; Liu et al., 2025). Proximity to rivers (<200 m) showed a strong positive correlation (FR = 2.44), as lateral river erosion continuously undercuts slope toes, disrupting mechanical equilibrium and triggering instability. The FR gradually decreased with increasing distance from rivers. High mine density (0.098–0.149) elevated risk (FR = 3.48) through slope disturbance and hydrological damage caused by mining activities. In contrast, stable factors such as hard limestone (FR = 0.16) and dense vegetation cover (>0.869, FR < 0.65) mitigated instability. Field surveys in Fuyuan County confirmed that areas with high vegetation coverage, characterized by tall trees and deep-rooted plants, effectively inhibited landslide occurrence. Their well-developed root systems penetrated weak layers and anchored deep rock and soil, significantly enhancing shear strength and overall slope stability, creating a “biological anchor” reinforcement effect.

The comparative evaluation of machine learning models highlighted the superiority of gradient boosting methods, among which the LightGBM model achieved the highest accuracy (84.10%) and balanced performance (F1 score = 83.17%), outperforming traditional statistical methods (SVM, LR) and artificial neural networks (ANN) — a finding consistent with previous studies (Danlu et al., 2023). This advantage stems from LightGBM’s ability to capture nonlinear relationships and higher-order interactions, as evidenced by its effectiveness in concentrating 60.67% of landslide points within 16.55% of the “very high” susceptibility zones. Although deep learning methods such as ANN are suitable for large-scale assessments, their potential is difficult to be fully realized in regions like Fuyuan County. The uncertainty quantification and stability results of the LightGBM model reflect the synergistic effect of data features, algorithm characteristics, and quantification methods: most samples exhibit low prediction standard deviation and narrow confidence intervals, benefiting from the strong discriminative power of the mainstream features in the dataset and the model’s effective capture of nonlinear relationships; the “U-shaped” correlation between average prediction probability and uncertainty reflects the consistency between the model’s decision-making logic and quantification indicators; the high AUC mean (0.964) and concentrated distribution (0.955–0.975) in Bootstrap validation are attributed to the reasonable data distribution and the model’s strong generalization ability. SHAP analysis clearly elucidated the model’s interpretability, showing that each influencing factor contributes to landslide occurrence but to varying degrees, with mining density identified as the most critical factor. The color gradient in Figure 15 indicates that, at a given mining density, areas closer to rivers typically correspond to higher SHAP values — a result of the combined effect of geomorphic factors (e.g., steep slopes and eroded landforms) that exacerbate slope instability, and geophysical processes (e.g., water-induced soil saturation and reduced shear strength near rivers) that amplify risks. Additionally, anthropogenic factors, particularly intensive mining activities, exacerbate surface disruption and weaken slope cohesion, forming a synergistic effect with natural factors to significantly increase landslide risk.

The high-resolution LSM exerts a profound practical impact on regional development, with its core value lying in transforming geological hazard prevention and control from passive post-disaster response to proactive pre-disaster prevention. Classifying regions into five susceptibility levels (very low, low, moderate, high, very high), the LSM serves as a key guideline for land-use planning, addressing local governments’ critical needs in urban development (Mandal and Mandal, 2018; Yeganeh et al., 2025; Cemiloglu et al., 2023). Specifically, governments can achieve dynamic zoning updates by integrating LSM results into mandatory planning revisions (e.g., designating ecological red lines in very high-risk areas to prohibit large-scale infrastructure and urban expansion, and setting differentiated construction standards involving anti-slide piles or drainage facilities in moderate-risk areas) and regularly updating zones using dynamic data such as seasonal precipitation and geological changes. Additionally, the model can be linked with GIS technology and real-time monitoring devices to build a full-process intelligent early warning system, generating multi-level alerts and supporting targeted disaster response plans. For mining regulation, the LSM acts as a core tool for license approval (rejecting applications in very high-risk areas and restricting mining scope in moderate-risk areas) and dynamic supervision via remote sensing, forming a “warning-inspection-penalty” closed-loop management system to enforce slope reinforcement and vegetation restoration.

Furthermore, if end-users (e.g., civil defense departments) have limited time or financial resources, they may prioritize focusing on truly critical influencing factors while temporarily setting aside secondary ones—a strategy supported by existing research findings in the region.

This study still has certain limitations. Due to constraints of surveying and mapping conditions as well as regional accessibility, the landslide inventory may have missed a small number of small-scale or hidden landslides. The 30 m spatial resolution makes it difficult to accurately characterize fine-grained features such as microtopography, imposing certain restrictions on the representation of local small-scale landslides. Precipitation and population density data are obtained through interpolation methods, whose accuracy depends on the distribution of original observation stations and interpolation algorithms, potentially leading to slight errors in areas with complex terrain or sparse data.

These limitations have a limited impact on the overall evaluation results. Future research can be further improved by optimizing data collection methods, enhancing spatial resolution, and refining interpolation techniques.