Landslide Susceptibility Mapping Based on Ensemble Learning in the Jiuzhaigou Region, Sichuan, China

Abstract

Accurate landslide susceptibility mapping is vital for disaster forecasting and risk management. To address the problem of limited accuracy of individual classifiers and lack of model interpretability in machine learning-based models, a coupled multi-model framework for landslide susceptibility mapping is proposed. Using Jiuzhaigou County, Sichuan Province, as a case study, we developed an evaluation index system incorporating 14 factors. We employed three base models—logistic regression, support vector machine, and Gaussian Naive Bayes—assessed through four ensemble methods: Stacking, Voting, Bagging, and Boosting. The decision mechanisms of these models were explained via a SHAP (SHapley Additive exPlanations) analysis. Results demonstrate that integrating machine learning with ensemble learning and SHAP yields more reliable landslide susceptibility mapping and enhances model interpretability. This approach effectively addresses the challenges of unreliable landslide susceptibility mapping in complex environments.

1. Introduction

Landslides are a widespread geological hazard, causing significant casualties and property damage annually, thus posing a serious threat to public safety and infrastructure. Factors such as accelerated glacial melting, erratic rainfall patterns, increasing extreme weather events, and intensified human activity are contributing to the rising incidence and severity of landslides [1]. Accurately identifying potential landslide locations is essential for mitigating their impacts and safeguarding communities.

Landslide susceptibility mapping (LSM) research typically falls into two categories: knowledge-driven and data-driven methods [2,3]. Knowledge-driven methods for LSM encompass methods like Analytic Hierarchy Process (AHP) [4,5], Fuzzy AHP [6], Fuzzy Relation AHP [7], and Fuzzy Logic [8,9], among others. Conversely, data-driven methods primarily include the frequency ratio [10], evidence weight [11], and various machine learning techniques such as logistic regression [12], support vector machines [13], random forests [14], and deep learning methods such as recurrent neural networks [15], deep neural networks [16], and convolutional neural networks [17].

Despite their efficacy, these models often suffer from a lack of interpretability, complicating the understanding of their decision-making processes. The complexity and scale of machine learning models render them “black boxes”, making it challenging to discern how inputs are transformed into outputs. Each geographical area exhibits unique geomorphological, topographical, hydrological, and anthropogenic characteristics, limiting the effectiveness of traditional single-model approaches. Thus, developing new methodologies that enhance model interpretability and improve prediction accuracy has become a critical research focus [18,19,20,21].

To thoroughly investigate the potential of ensemble learning and improve model interpretability, this study combines diverse datasets and employs both traditional and ensemble learning techniques for landslide susceptibility mapping in Jiuzhaigou County, Sichuan Province. A SHAP analysis is utilized on the trained models to identify key factors and their values influencing each landslide, thereby enhancing model interpretability.

2. Study Area and Data

2.1. Study Area

Jiuzhaigou County, located in the northwestern Aba Tibetan and Qiang Autonomous Prefecture of Sichuan Province, China, lies between approximately 32°54′ to 33°19′N latitude and 103°46′ to 104°4′E longitude (Figure 1). This area, situated at the juncture of the Qinghai–Tibet Plateau and the Sichuan Basin, features a monsoon-dominated and semi-humid climate with significant topographical variation—elevations range from low in the east to high in the west, creating a relative height difference of 3712 m [22].

The region exhibits folded fractures due to intense neotectonic activity, primarily in the Minjiang, Huya, Xueshanliangzi, and Tazang Fracture Zones, with ground shaking peak accelerations of 0.04 g to 0.26 g. Jiuzhaigou County is prone to geological disasters, having experienced over 50 earthquakes of magnitude 5.0 or higher in the past century, according to Fan et al. [23]. The average annual rainfall is 552.3 mm, contributing to its susceptibility to landslides and other geological hazards.

2.2. Database

2.2.1. Historical Landslides

In the initial phase of the LSM process, we collected high-precision field observations from local geological surveys. By the end of 2020, we documented 164 significant slope failure events through field investigations, historical records, and satellite imagery interpretation. These data include medium-sized landslides, avalanches, debris flows, and unstable slopes, forming the basis of our landslide inventory (Figure 1).

Given the scarcity of detailed data on landslide boundaries, we represented each event by its center point. To mitigate the challenge of imbalanced training samples, we performed random sampling to match an equal number of non-landslide points. A 100 m buffer zone was established around landslide points, and non-landslide points were sampled outside these zones [24].

2.2.2. Landslide Conditioning Factors

The success of LSM hinges on the careful selection of conditioning factors [25]. Among these, altitude plays a critical role in the formation and distribution of landslides [26,27]. Using ArcGIS 10.8, we analyzed Digital Elevation Model (DEM) data to compute six key conditioning factors: the slope, aspect, plan curvature, profile curvature, topographic wetness index (TWI), and stream power index (SPI). A buffer analysis was employed to derive distance factors related to roads, historical earthquakes, faults, and human activities [28]. Detailed information about the data sources and nomenclature for these factors is provided in

Table 1. Description of lithology of study area.

Legend
Number	Lithology	Number	Lithology
a	Holocene Alluvial	j	Ice Peak Group
b	Pleistocene Glaciers Accumulate	k	Qiongyi Group
c	Oligocene Geru Group	l	Late Carboniferous Hollowa Group
d	Baoding Group or Baiguowan Group	m	Gangou Group
e	Sanzhushan Group–Maso Mountain Group	n	Posongchong Group–Qujing Group
f	Qingtianbao Group–Yanyuan Group	o	Early Devonian Segala Group
g	Dashibao Group–Pinecigou Group	p	Lou Shanguan Group
h	Permian Mafic Rocks	q	Wooden Seat Group–Crystal Group
i	Gunda Overview Group	l	Yang Tianba Group

and

Table 2. Factors utilized in the present study.

Source	Time	Scale	Format	Factor
Landslides		-	Shapefile	-
ASTER DEM (https://earthexplorer.usgs.gov/)	2011	30 m	TIFF	Slope
				Aspect
				Profile curvature
				Plan curvature
				SPI
				TWI
Underlying geographic information data (https://www.webmap.cn/)	2021	1:1,000,000	Shapefile	Distance to road
				Distance to human active point
History earthquake point (https://data.earthquake.cn/)	2010–2022	-	Shapefile	Distance to history earthquake point
Geological map	2017	1:1,000,000	TIFF	Distance to faults
				Lithology
LULC map (http://data.ess.tsinghua.edu.cn/)	2020	30 m	TIFF	LULC
Landsat 8	2022	30 m	TIFF	NDVI
Chinese Academy of Sciences, Center for Resource and Environmental Data and Sciences	2002–2022	1000 m	TIFF	Annual precipitation

, as well as Figure 2.

Slope influences pressure distribution and sediment accumulation [2], while variations in the solar radiation of different slopes impact landslide stability [29]. The profile and plan curvature elucidate the spatial geometry of unstable slopes; profile curvature affects flow velocity, whereas plan curvature influences flow direction [30,31].

The TWI is associated with slope and water flow, indicating areas with potentially higher soil moisture. The SPI indicates river erosive capacity, contributing to slope instability through erosion [32]. Lithological variations in the Jiuzhaigou region affect slope stability due to differences in strength, weathering resistance, and permeability [33]. The lithology map for the Jiuzhaigou region is presented in Table 1.

The normalized difference vegetation index (NDVI) serves as an indicator of vegetation growth, which indirectly affects slope stability by influencing rock and soil erosion [34]. The proximity to faults and historical seismic sites provides valuable insights into regional tectonic activity and fault instability, both critical factors in landslide formation [35]. Additionally, distances to roads, human activity points, and land use/land cover (LULC) data reflect the effects of anthropogenic activities on slope stability [36]. Meteorological factors, especially annual precipitation, play a key role in triggering landslides [27].

After a comprehensive analysis, 14 landslide conditioning factors were selected. These factors were converted to raster format and resampled to a 30 m resolution. Each factor was categorized using established classifications for lithology and LULC, while other factors were classified based on expert knowledge to optimize value arrangements (Figure 2). Continuous variables were manually classified according to expert insights and previous studies [3,26,27,28]. To normalize each factor for the analysis, the maximum–minimum method was applied, scaling values from 0 to 1, ensuring a zero mean and unit variance for each feature.

3. Methods

The study comprises three stages, as illustrated in Figure 3. First, we conducted comprehensive data collection, including conditioning factors and historical landslide inventories. Next, we developed three conventional machine learning models and four ensemble models to create an extensive LSM. In the final stage, we thoroughly evaluated and discussed the resulting susceptibility maps. This systematic approach facilitated a detailed assessment and interpretation of our findings, contributing valuable insights to the field of landslide susceptibility research.

3.1. Evaluation of Conditioning Factors

Incorporating additional environmental variables and triggers enhances the understanding of landslide susceptibility but also increases computational burden of the model. Excessive redundant information can obscure the importance of key factors, potentially distorting results [37]. Therefore, optimizing the selected assessment factors is essential, which includes diagnosing multicollinearity.

To evaluate multicollinearity, we employed the Variance Inflation Factor (VIF) and tolerance (TOL) metrics. This analysis examined the linear correlations among selected factors, identifying redundancies that could compromise model accuracy.

3.2. Base Classifiers

3.2.1. Support Vector Machines

A support vector machine (SVM) is used for binary classification through supervised learning, leveraging a kernel trick for classification [38]. The fundamental idea is to identify the optimal hyperplane that maximizes the distance between data points of the two classes, thereby ensuring better generalization to unseen data. The effectiveness of SVMs is highly dependent on the choice of kernel and regularization parameters, which are typically tuned through cross-validation to optimize performance.

Assuming a training set $x$ = {$x_1,x_2,···,x_n$} ∈ R^B×N, the hyperplane can be minimized as follows:

$${\displaystyle \frac{1}{2}}{\Vert w\Vert }^2+C{\displaystyle \sum }_i{\xi }_i,s.t.y_i\left(w·x_i+b\right)\ge 1-{\xi }_i$$

where $\Vert w\Vert$ represents the norm of the hyperplane normal vector, $w$ is the weight vector, $b$ denotes the bias, ${\xi }_i$ are the slack variables for non-separable data, and $C$ is the penalty parameter. The cost function is defined as follows:

$${\displaystyle \sum }_{i=1}^n{a}_i-{\displaystyle \frac{1}{2}}{\displaystyle \sum }_{i=1,j=1}^n{a}_i{a}_j{y}_i{y}_{j}K\left(x,x_i\right),s.t.0\le a_i\le C,{\displaystyle \sum }_{i=1}^l{a}_i{y}_i=0$$

where $a_i$ are Lagrange multipliers and $K\left(x,x_i\right)$ is the kernel function.

3.2.2. Logistic Regression

Logistic regression (LR) is a generalized linear regression analysis model that falls under supervised learning in machine learning. The LR model is constructed as follows:

$$p={\displaystyle \frac{1}{1+e^{-y}}}$$

where $p$ is the probability of a landslide occurring and $y$ denotes a linear combination function as follows:

$$y=b_0+b_1{x}_1+b_2{x}_2+···+b_n{x}_n$$

where ($x_1,x_2,···,x_n$) denote the explanatory variables affecting the occurrence of landslides, $b_0$ is a constant value, and ($b_1,b_2,···,b_n$) denote the regression coefficients.

3.2.3. Gaussian Naive Bayes

Gaussian Naive Bayes (GNB) is a supervised learning algorithm that directly estimates the probabilistic relationship between labels and features. By applying Bayes’ theorem, GNB calculates the posterior probability of a sample belonging to each class based on its feature distribution. GNB is characterized by its stable classification efficiency, strong performance even with small datasets, and reduced sensitivity to missing data. The formula for Gaussian Naive Bayes with each feature’s values is as follows:

$$P\left(x_i|y_c\right)={\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{ci}^2}}}\exp \left(-{\displaystyle \frac{{\left(x_i-{\mu }_{ci}\right)}^2}{2{\sigma }_{c_i}^2}}\right)$$

where $x_i$ denotes the ith feature dimension, and ${\sigma }_{ci}$ and ${\mu }_{ci}$ denote the standard deviation and expectation corresponding to the features under the category.

3.3. Ensemble Learning Methods

3.3.1. Bagging

Bagging was initially proposed for integrating classification and regression trees [39]. It employs bootstrap sampling to create multiple subsets from the original training set, with each subset training an individual base model. This results in several base classifiers, each operating on different data. During prediction, the outputs from these models are aggregated using a Voting mechanism, enhancing model stability and accuracy by reducing variance and mitigating overfitting.

3.3.2. Voting

Voting is an ensemble learning technique for classification, utilizing two strategies: simple averaging and weighted averaging [40]. Simple averaging combines predicted probabilities from all base classifiers to generate a final prediction. In contrast, weighted averaging assigns greater importance to models with higher AUC values, which reflects their effectiveness in distinguishing classes, potentially leading to more accurate predictions.

3.3.3. Boosting

Boosting improves predictive performance by combining homogeneous base learning algorithms through dataset reweighting and updates across multiple models [41,42]. It begins by training the first classifier on the initial dataset and evaluating its error rate. The dataset is then adjusted based on this error, and the reweighted data are used to train the next classifier. This iterative process continues, ultimately aiming to construct a strong ensemble model. Final predictions are made by aggregating the outputs of all classifiers, typically using a weighted combination reflecting each model’s performance.

3.3.4. Stacking

Stacking, introduced by Wolpert (1992), is a heterogeneous ensemble learning method that employs meta-learning to integrate multiple base classifiers [43]. The Stacking framework consists of two components: a Basic Learner (Level 0) and Meta-Learner (Level 1). During training, the dataset is split into two subsets to develop the base and meta-learners. In the Base Level Training stage, various base models are trained independently on the original dataset, with their predictions serving as inputs for the Meta-Level Training stage [39]. Here, a meta-model is trained using the outputs from the base models as features.

3.4. Model Evaluation Measures

Accurately assessing model performance is crucial for evaluation. The receiver operating characteristic curve (ROC) is a widely used method for this purpose, with the area under the ROC (AUC) indicating the model’s predictive capability for unseen data; higher AUC values signify better predictions.

By setting a threshold (commonly 0.5), continuous predictions can be converted into binary classifications. Evaluation metrics such as the accuracy, Matthews Correlation Coefficient (MCC), and confusion matrix are typically employed to gauge predictive performance [44,45]. For binary classification tasks, outcomes can be categorized into four scenarios, as shown in

Table 3. Four cases of prediction results of binary classification model.

Confusion Matrix		Predicted Value
		Positive Class (1)	Negative Class (0)
True value	Positive Class (1)	True Positive (TP)	False Negative (FN)
	Negative Class (0)	False Positive (FP)	True Negative (TN)

.

Kappa measures categorization accuracy and consistency. While precision and recall focus on the positive class, neglecting the negative class, the MCC offers a balanced assessment by considering both classes equally, making it particularly useful for imbalanced datasets. The evaluation metrics are computed based on the confusion matrix as follows:

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

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

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

$$MCC={\displaystyle \frac{TP\times TN-TP\times FN}{\sqrt{\left(TP+FP\right)\times \left(TP+FN\right)\times \left(TN+FP\right)\times \left(TN+FN\right)}}}$$

$$Kappa={\displaystyle \frac{P_o-P_e}{1-P_e}}$$

Among

$$P_o=Accuracy$$

$$P_e={\displaystyle \frac{\left(TP+FP\right)\times \left(TP+FN\right)+\left(FP+TN\right)\times \left(FN+TN\right)}{{\left(TP+FP+FN+TN\right)}^2}}$$

3.5. SHapley Additive exPlanations

SHAP is a unified framework for interpreting machine learning models by providing consistent and interpretable feature importance scores [46]. Formally, the SHAP value ${\varphi }_i$ for a feature $i$ is computed using the formula

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

where $N$ denotes the set of all features, $S$ represents a subset of features excluding $i$, and $f\left(S\right)$ is the model’s prediction with the feature subset $S$. This formula computes the average marginal contribution of feature $i$ across all possible feature subsets, ensuring a fair and comprehensive assessment of its impact.

4. Results

4.1. Multicollinearity

We assessed multicollinearity among the 14 selected assessment factors using the VIF and TOL methods. Research indicates that factors are considered independent when VIF is below 10 and TOL is above 0.1 [18,47]. All 14 factors exhibited VIF values below 10 and TOL values exceeding 0.1, confirming their suitability for landslide susceptibility modeling and analyses (Figure 4).

4.2. Model Performance and Evaluation

The performance of each model was assessed using the training set, as illustrated in

Table 4. Evaluation of landslide models using training dataset.

	Model	Accuracy	Precision	AUC	Recall	Kappa	Matthews
Base classifiers	SVM	0.8672	0.8952	0.9267	0.8318	0.7345	0.7363
	LR	0.8230	0.8119	0.9020	0.8407	0.6560	0.6464
	GNB	0.8761	0.9126	0.9174	0.8318	0.7522	0.7551
Ensemble learning methods	Voting	0.9070	0.9259	0.9789	0.8849	0.8141	0.8149
	Bagging	0.9380	0.9380	0.9847	0.9380	0.8761	0.8762
	Stacking	0.9161	0.9126	0.9516	0.8718	0.7922	0.8081
	Boosting	0.9336	0.9437	0.9789	0.9115	0.8672	0.8681

and Figure 5. SVM served as the base model for Bagging and Boosting due to its highest verification accuracy. For Voting and Stacking, three models were employed as base classifiers. The results show that all models exhibit strong performance, with SVM achieving the highest AUC of 0.9267. Among the ensemble methods, Bagging outperformed Boosting, Voting, and Stacking. Integrating SVM with Bagging improved accuracy by 7.08%, precision by 4.28%, AUC by 5.80%, recall by 10.62%, Kappa by 14.16%, and Matthews by 13.99%. Boosting yielded improvements of 6.64%, 4.85%, 5.22%, 7.97%, 13.27%, and 13.18%, respectively. Stacking and Voting also showed enhancements, achieving AUC values of 0.9516 and 0.9789.

Subsequently, the generalization ability of the model was evaluated (

Table 5. Evaluation of landslide models using validation dataset.

	Model	Accuracy	Precision	Auc	Recall	Kappa	Matthews
Base classifiers	SVM	0.8571	0.9069	0.9192	0.7959	0.7142	0.7197
	LR	0.8265	0.8200	0.9042	0.8367	0.6530	0.6531
	GNB	0.8571	0.9069	0.9017	0.7959	0.7142	0.7197
Ensemble learning methods	Voting	0.8775	0.9092	0.9296	0.8571	0.7551	0.7557
	Bagging	0.8775	0.9243	0.9340	0.8571	0.7551	0.7557
	Stacking	0.8673	0.9285	0.9291	0.8124	0.7522	0.7551
	Boosting	0.8857	0.8869	0.9180	0.8163	0.7628	0.7629

and Figure 5). Among ensemble models, the optimal values for accuracy, precision, AUC, recall, Kappa, and Matthews were 0.8857, 0.9285, 0.9340, 0.8571, 0.7628, and 0.7629, respectively, with Bagging demonstrating the best predictive performance. Voting followed closely, while Stacking ranked moderately, and Boosting provided the least enhancement to predictive capabilities.

Each ensemble strategy—Bagging, Boosting, Stacking, and Voting—offers unique advantages and limitations for landslide susceptibility mapping. Bagging reduces variance by averaging predictions from correlated models, as shown by the 5.80% AUC improvement over SVM. Boosting, such as XGBoost, improves accuracy by sequentially minimizing the loss function, although it may risk overfitting. Stacking effectively utilizes outputs from multiple base learners through a secondary model for predictions, demonstrating improved accuracy. Voting aggregates predictions via majority or weighted voting, enhancing robustness against overfitting. Overall, the findings indicate that while diverse independent models minimize errors, the careful consideration of each ensemble strategy is essential for optimizing predictive performance in complex environments.

4.3. Landslide Susceptibility Mapping

Using well-trained landslide susceptibility models, we calculated a susceptibility index ranging from 0 to 1, with higher values indicating greater landslide susceptibility. To enhance interpretability, we applied the natural breaks method to classify this index into five categories, providing a clearer depiction of susceptibility distribution for informed decision making (Figure 5) [33,48,49]. The generated susceptibility maps reveal consistent spatial patterns across methodologies, with areas classified as very high susceptibility aligning with regions historically prone to landslides (Figure 6).

The results demonstrate that both machine learning and ensemble learning techniques exhibit similar spatial distribution characteristics in landslide susceptibility mapping. High- and very-high-susceptibility zones correlate with significant landslide occurrences, indicating robust predictive capabilities aligned with terrain features. We also computed landslide density within these five susceptibility zones (

Table 6. Calculation results of statistical results (%) and landslide density (no./km²) in each landslide sensitivity zone of LSM models.

Level	Very Low		Low		Moderate		High		Very High
	Density	Ratio	Density	Ratio	Density	Ratio	Density	Ratio	Density	Ratio
SVM	0.0249	46.23	0.0867	17.76	0.0992	18.56	0.1577	10.22	1.0223	7.23
LR	0.0071	40.93	0.0323	27.23	0.1235	16.52	0.2438	9.22	1.5450	6.10
GNB	0.0035	60.30	0.0218	27.12	0.1280	12.22	0.2022	6.32	1.4233	6.26
Voting	0.0013	41.72	0.0316	30.24	0.1273	12.83	0.2012	7.12	1.2011	6.63
Boosting	0.0025	52.23	0.0225	19.60	0.1125	11.29	0.2532	7.26	1.4324	9.62
Bagging	0.0004	42.57	0.0023	18.31	0.0183	21.95	0.1943	11.29	1.6792	5.68
Stacking	0.0020	56.07	0.0113	17.48	0.0204	15.22	0.1822	5.21	1.6172	6.02

).

High-sensitivity areas predominantly cluster along roads, rivers, and valleys. Notably, landslide-prone areas (very-high- and high-susceptibility zones) encompass only 16% of the study area, yet they account for nearly all historical landslides. In Jiuzhaigou County, intensive human activities, particularly highway construction, exacerbate landslide risk by increasing slope exposure. Previous studies highlight the need for extensive infrastructure due to rising tourism, with 138 geological disasters recorded in Zhangzha Town, 85 of which occurred in scenic areas [50]. Thus, future research should focus on balancing sustainable tourism development in Jiuzhaigou while implementing effective strategies to mitigate landslide risks and ensure local resident safety.

5. Discussion

5.1. SHAP Model Interpretation

This study utilized the SHAP interpreter to elucidate the Bagging model and analyze how various factors influence landslide spatial distribution (Figure 7). Figure 7a presents the importance ranking of factors based on average absolute SHAP values, identifying distance to roads, rainfall, and lithology as the most significant contributors to landslides. Figure 7b illustrates the SHAP values for each feature.

The results indicate that proximity to roads correlates with lower slope stability, while areas further from roads exhibit greater stability. Annual average rainfall initially promotes landslide occurrences, but this effect diminishes at higher rainfall levels. Human activities, including deforestation, urbanization, and mining, significantly destabilize slopes. Deforestation reduces root systems that anchor soil, increasing erosion and landslide risk. Urbanization leads to increased runoff and soil saturation, further destabilizing slopes. Mining creates steep, unsupported slopes and disrupts natural drainage, heightening failure risks.

Excessive rainfall saturates soil, weakening its stability and increasing susceptibility to sliding. Meanwhile, moderate rainfall can trigger landslides, and excessive saturation results in reduced cohesion. Heavy rainfall raises pore water pressure, destabilizing slopes, while erosion can undermine surface integrity. Beyond a certain threshold, increased rainfall is associated with a decreased likelihood of landslides [51].

Lithology also significantly impacts landslide development by influencing slope stability, water retention, and erosion. Different soil types affect water retention; for example, clay and silt can increase pore pressure, reducing stability, while sandy soils may become loose when wet. Rock types further influence stability: igneous rocks tend to be stable, while sedimentary and metamorphic rocks may exhibit weaknesses due to layered structures or foliation. Weathering processes differ among lithologies, producing loose materials that impact slope integrity.

Local interpretations often provide more precise insights than global explanations for specific regions. The waterfall chart in Figure 8 illustrates the magnitude and trend of each variable within a selected landslide prediction unit. The vertical axis represents the landslide impact factors and their values, while the horizontal axis shows the corresponding SHAP values.

The chart indicates that distance to roads and NDVI are the most critical variables affecting landslide susceptibility. Notably, the rankings for rain and lithology differ when considering individual landslides compared to the overall context.

5.2. Limitation and Future Work

This study recognizes the influence of expert knowledge in the Holdout method for data partitioning. Although we used a one-to-one ratio for landslide and non-landslide data, we did not adequately address the potential imbalance between these datasets, which could affect their contributions to the model.

Additionally, the analysis utilized grid cells, which may inadequately represent geological, topographical, and hydrogeological features. Future research should consider using slope units derived from terrain segmentation, as they more effectively capture the factors influencing landslides [52].

Model interpretability is crucial in landslide research due to the significant safety and economic implications of decisions based on these models. Beyond SHAP, techniques like Class Activation Mapping (CAM) can further enhance interpretability [53]. Future studies should explore the physical mechanisms underpinning the relationships identified by these models to strengthen their practical applications.

6. Conclusions

This study leverages ensemble learning with SHAP to develop an interpretable LSM model based on 164 landslides and 14 factors in Jiuzhaigou County, Sichuan Province, achieving optimal landslide susceptibility assessment with clear insights into factor importance. Our findings reveal that ensemble learning enhances the accuracy and robustness of individual models, with SHAP clarifying the roles of precipitation, lithology, and distance to roads in landslide occurrence. This integration of interpretability and ensemble learning offers a novel approach to regional LSM. Future work could explore more advanced, interpretable machine learning models to further advance the field.