Landslide Susceptibility Assessment in Yulong County Using Contribution Degree Clustering Method and Stacking Ensemble Coupled Model Based on Certainty Factor

Abstract

To address the subjectivity of traditional factor attribute grading methods and the weak predictive capabilities of single-model classifications, this study focused on Yulong County; the Contribution Degree Clustering Method (CDCM) utilizes the Certainty Factor (CF) as the contribution index to partition continuous factor attribute intervals. Additionally, the Sparrow Search Optimization Algorithm (SSA) is employed for hyperparameter tuning. The CF is incorporated into Support Vector Machine (SVM), Back Propagation Neural Network (BPNN), and Random Forest (RF) models to form the CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF coupling models, respectively. These basic coupling models are further integrated using the Stacking algorithm to create the CF-SSA-Stacking integrated coupling model for constructing a landslide susceptibility assessment system. The results indicate that the CF-SSA-Stacking integrated coupling model achieves the highest accuracy, F1 score, Kappa coefficient, and AUC value, with values of 0.89375, 0.89172, 0.787500, and 0.9522, respectively. These metrics are significantly superior to those of the three basic coupling models, demonstrating better generalization capability and reliability. This suggests that the model can identify more historical landslide occurrences using fewer grid areas classified as extremely-high- or high-susceptibility zones. It is suitable as an effective regional landslide susceptibility assessment method for practical disaster prevention and mitigation applications. Further studies could explore the model’s performance across varying geological settings or with different datasets, providing a roadmap for future research and development in landslide susceptibility assessment.

1. Introduction

Landslides, a frequent geological disaster, have increased in frequency and scale in recent years, primarily influenced by seismic activity, extreme weather events induced by climate change, and accelerated human engineering activities. These factors not only reactivate old landslides but also trigger new ones, posing severe threats to lives, property, and critical infrastructure safety [1]. Given these increasing threats, landslide susceptibility assessment, as an early warning measure, not only aids in understanding regional landslide activity patterns but also predicts the spatial locations of potential landslides. This approach can effectively mitigate the risks to human life and property, providing valuable insights for disaster prevention, mitigation, and informed government decision-making [2].

The results of landslide susceptibility assessment are influenced by multiple factors, with the choice of evaluation method being a crucial component. Current evaluation methods predominantly include empirical models, mathematical statistical analysis models, and machine learning models [3]. Previous studies have indicated that empirical models rely on expert knowledge and geological insights to identify landslide locations and potential risk areas. While these models offer practical and intuitive insights, they are often highly subjective and time-consuming and may lack accuracy and consistency when applied to complex geological conditions [4]. Mathematical statistical analysis models, such as the frequency ratio (FR), Information Value (IV), and Certainty Factor (CF) models, quantify the relationship between landslide influencing factors and landslide occurrences, providing strong interpretability. However, these models may struggle to fully capture the nonlinear relationships among multiple influencing factors when dealing with complex multidimensional data [5]. Machine learning models are better equipped to reveal the intricate nonlinear relationships between landslide development and triggering factors, often yielding superior performance in landslide susceptibility assessment. Nevertheless, as many machine learning models are “black-box” models, their results tend to be less interpretable and are highly dependent on the quality of the dataset [1,6]. Li et al. (2022) demonstrated that by using optimization algorithms such as Simulated Annealing (SA), Particle Swarm Optimization (PSO), and Sparrow Search Algorithm (SSA) to enhance a Deep Belief Network (DBN); the SSA-DBN model achieved the highest assessment accuracy [7]. Recently, the application of model coupling and ensemble methods in landslide susceptibility assessment has increased [1,8]. Ensemble algorithms, which integrate the predictive outcomes of multiple models, exhibit stronger robustness against noisy data and outliers, making them more broadly applicable [8]. Among these, the Stacking ensemble model has shown greater generalization ability in landslide susceptibility assessment compared to individual base models [9]. For example, studies comparing the predictive accuracy of Stacking, Bagging, and Boosting models when integrating Support Vector Machines (SVMs) for susceptibility assessment have found that the Stacking ensemble model exhibits higher accuracy when combining different models. However, when all base models are homogeneous, its predictive performance may not surpass that of Bagging and Boosting models [10]. Zhou (2023) compared the accuracy of Stacking ensemble models with that of Boosting ensemble models, concluding that Stacking achieves higher accuracy when integrating heterogeneous models [11]. Model coupling, which combines traditional statistical models with machine learning models, leverages the ability of statistical models to assess the contribution of different factor intervals to landslide occurrence during attribute interval grading. This information can then be used to interpret feature importance or correlations within machine learning models, thus enhancing predictive accuracy while balancing model interpretability and prediction capability. For instance, Yuan et al. (2022) demonstrated that coupling machine learning models with the Certainty Factor model yields higher accuracy than using standalone machine learning models [12]. Current research on coupling models with Stacking ensembles has primarily focused on the frequency ratio and Information Value models, with little exploration into the effects of coupling the Certainty Factor (CF) model with Stacking ensemble models for susceptibility assessment [9,11].

In addition, the method of grading attribute intervals for evaluation factors is another factor influencing the results of landslide susceptibility assessment. Evaluation factors can be classified into two categories: discrete and continuous. Previous studies have typically graded discrete factors based on their natural attribute intervals, while various methods have been developed for continuous factors. Methods such as expert judgment [13], natural breaks [14], and equal interval [15] approaches grade attribute intervals based on predetermined steps. However, these methods may introduce biases due to subjective judgment errors and the particularities of data distribution. On the other hand, methods like the Contribution Degree Clustering Method (CDCM) [16], which relies on the results of mathematical statistical analysis models to cluster and merge similar intervals, can mitigate the impact of unreasonable interval grading on accuracy. For instance, methods such as the CDCM that combines the frequency ratio method and equal interval method for attribute interval classification have been shown to yield more accurate results than the natural breaks method for attribute interval classification [16]. However, this method relies on constructing histograms to perform fuzzy interval clustering and merging, which can be influenced by subjective factors to some extent. To minimize errors arising from subjective judgments, one could attempt to discretize each evaluation factor into as many equally spaced small intervals as possible and then perform equal interval clustering and merging based on the contribution results, rather than relying on fuzzy interval clustering. Moreover, no studies have explored the effectiveness of using the Contribution Degree Clustering Method as a grading approach for attribute intervals within Stacking models in landslide susceptibility assessments.

In summary, this study takes Yulong Naxi Autonomous County as an example, selecting 14 indicators including slope and lithology, to perform attribute interval classification for continuous factors using the Contribution Degree Clustering Method, with the Certainty Factor (CF) as the measure of contribution. For discrete factors, natural attribute intervals are used for classification. This study then employs the Sparrow Search Algorithm (SSA) for parameter optimization and adopts a Stacking ensemble model for landslide susceptibility assessment. This ensemble model consists of three heterogeneous machine learning models—Random Forest (RF), Support Vector Machine (SVM), and Back Propagation Neural Network (BPNN)—as base classifiers, with a Logistic Regression (LR) model serving as the meta-classifier. By comparing the landslide susceptibility prediction results and model accuracy of the coupled models based on CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF, with those of the CF-SSA-Stacking ensemble model, this study explores the applicability of the CF-SSA-Stacking ensemble model for landslide susceptibility assessment when the Contribution Degree Clustering Method is involved in attribute interval classification. The objective is to enhance the accuracy and practicality of landslide susceptibility assessment, thereby providing a new perspective and method for such evaluations.

2. Study Area and Datasets

2.1. Study Area

Yulong Naxi Autonomous County (hereinafter referred to as “Yulong County”) is located in the northwest of Yunnan Province, China (Figure 1), with geographical coordinates ranging from 99°23′ to 100°32′E and 26°34′ to 27°46′N. It covers a total area of 6198.76 km², with mountainous areas accounting for 94.4% of the total area. Situated at the transitional zone between the Hengduan Mountains on the southeastern edge of the Qinghai–Tibet Plateau and the Yunnan–Guizhou Plateau, Yulong County lies at the junction of the Tibetan–Yunnan Geosyncline and the Yangtze Paraplatform, featuring geomorphological characteristics of both the Hengduan Mountain Gorges and the Western Yunnan Plateau. The terrain is high in the northwest and low in the southeast, belonging to a high mountain gorge region. The county exhibits significant vertical relief and has typical three-dimensional climate characteristics, with an average annual rainfall of 968 mm and notable regional differences in rainfall distribution. An analysis of historical geological disaster records reveals 240 landslide hazard points, many of which are in the creeping or low-speed sliding stages. These landslides are among the main geological disasters threatening the safety of lives and property in the region. The geological structure within Yulong County is well developed, with frequent tectonic movements, widespread faults and folds, frequent geological disasters, and severe activity hazards [17].

2.2. Data Source

The data sources are listed below:

• Historical Landslide Inventory Data were obtained from the 2021 Geological Hazard Risk Survey Project (1:50,000 scale) in Yulong Naxi Autonomous County, Yunnan Province (Figure 1 blue dots);
• DEM data were sourced from the ALOS satellite’s AW3D30 (ALOS World 3D) 30 m dataset, available from JAXA, Japan Aerospace Exploration Agency (https://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm, accessed on 11 January 2024), used for extracting factors such as the slope, aspect, planar curvature, profile curvature, topographic wetness index (TWI), stream power index (SPI), and topographic position index (TPI);
• The Road and River Data from the 2021 Geological Hazard Risk Survey Project of Yulong Naxi Autonomous County, Yunnan Province (1:50,000), and the fracture data from the 1:200,000 geologic data were processed by Euclidean distance to extract the evaluation indexes of distance to a road, river, and fault;
• Annual Average Rainfall factors were obtained through the Kriging interpolation of precipitation data from Yulong County Meteorological Bureau;
• Vegetation cover assessment factors were extracted from the 30 m Annual China Land Cover Dataset (CLCD) [18];
• Slope Length Factor (LR) was derived from DEM data using the Terrain Factors Calculation Tool (Version 2.0) of the Soil Erosion Model [19];
• Lithological Data from 1:200,000 geological data were categorized into eight types based on regional lithological characteristics: dense, massive, hard extrusive rock units; moderately thick to thick bedded weak to moderately karstified, relatively hard carbonate rock units; thin to moderately thick bedded, relatively soft sandstone–mudstone to relatively hard conglomerate; thin to moderately thick bedded, relatively hard sedimentary metamorphic rock units; loose to medium-dense mixed soil; massive to blocky, hard intrusive rock units; thin to moderately thick bedded, weak mudstone, sandstone, and conglomerate units; thin to moderately thick bedded, relatively soft sedimentary metamorphic rock units.

3. Methods

3.1. Research Methodology

The research methodology of this study involves the following steps: (1) selecting evaluation factors and classifying attribute intervals; (2) constructing the dataset and screening evaluation index factors; (3) building the evaluation model and partitioning the model’s prediction results; (4) conducting an accuracy evaluation of the partitioned results, analyzing the distribution pattern of indices, and performing a SHAP-based interpretability analysis (Figure 2). Specifically, this study employs a contribution clustering method, using the Certainty Factor (CF) as the contribution metric to classify the continuous factor attributes. Subsequently, the Sparrow Search Algorithm (SSA) is applied for hyperparameter optimization, and the Certainty Factor is coupled with a Stacking ensemble model to construct an integrated coupling model for landslide susceptibility evaluation in Yulong County. Finally, the model’s applicability is assessed, and an interpretability analysis is performed based on the prediction results.

3.2. Evaluation Factor Selection and Attribute Interval Grading

3.2.1. Evaluation Factor Selection

The selection of basic evaluation units includes grid cells, slope units, watershed units, etc. Among them, grid cells are chosen as the fundamental evaluation unit due to their convenient handling and simple data structure. Compared to other methods, they are more suitable for research involving machine learning models. Therefore, this study selects 30 m × 30 m resolution grid cells as the basic evaluation unit [3].

Based on previous research experience, 11 continuous evaluation factors have been preliminarily selected, including the slope, planar curvature, profile curvature, topographic wetness index (TWI), stream power index (SPI), distance to river, distance to roads, distance to fault, vegetation coverage, slope length steepness factor (LR), and annual average precipitation, along with 3 discrete evaluation factors, the aspect, topographic position index (TPI), and lithology, totaling 14 evaluation factors (Figure 3) [20,21].

3.2.2. Evaluation Factor Attribute Interval Grading

For the preliminary selection of 14 evaluation factors, the Certainty Factor model is used to calculate the certainty results of each factor in small intervals as a contribution indicator for grading the intervals of continuous factors.

• Certainty Factor Model:

The Certainty Factor (CF) model is widely applied in landslide susceptibility assessment studies, similar to frequency ratios. It is used to handle the heterogeneity and uncertainty of combining different data layers and input data [22]. The formula is

$$CF=\left\{\begin{array}{cc}{\displaystyle \frac{p{p}_a-p{p}_s}{p{p}_a\left(1-p{p}_s\right)}}& p{p}_a\ge p{p}_s\\ {\displaystyle \frac{p{p}_a-p{p}_s}{p{p}_s\left(1-p{p}_a\right)}}& p{p}_a\le p{p}_s\end{array}\right.$$

(1)

In Equation (1), CF represents the Certainty Factor of a specific class of an evaluation factor, which is the conditional probability of a landslide event occurring under a certain classification of the factor. $p{p}_a$ is the ratio of the number of landslide points in a certain classification a of the evaluation factor to the area of that classification unit; $p{p}_s$ is the ratio of the total number of landslides in the study area to the total area of the study area. By using the CF model, each class of each evaluation factor is assigned a value that varies within the interval [−1, 1]. A positive value closer to 1 indicates an increase in the certainty of landslide occurrence, while a negative value closer to −1 indicates a decrease in certainty. Values close to 0 suggest that it is difficult to provide information on the certainty of landslide occurrence.

To eliminate the impact of the inconsistent dimensions among factors and the errors caused by the presence of both positive and negative Certainty Factors, the CF values of each category of all conditional factors are normalized to determine the factor category weights ($f_{cf}$). The calculation formula is

$$f_{cf}={\displaystyle \frac{\raisebox{1ex}{$\left(S_i+1\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{{\displaystyle {\sum }_{i=1}^n\left(\raisebox{1ex}{$\left(S_i+1\right)$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}}}$$

(2)

In Equation (2), $S_i$ represents the CF value of each level of the evaluation factor.

2.

Factor Attribute Interval Grading:

Due to the assumption that geological conditions are similar in similar regions, the contribution results of mathematical statistical models are used as the similarity measure criteria for factor interval clustering. There is currently no unified naming standard for the method of recalculating the contribution of clustered intervals after merging similar intervals of continuous-type factors to determine the final grading results of continuous-type factor attribute intervals. Therefore, this paper uses the Contribution Degree Clustering Method (CDCM) as a general term for such methods.

In contrast to the approach proposed by Ke et al. (2023), which involves interval grading of factor attributes using equidistant intervals determined by expert judgment and applies the frequency ratio as a contribution index in contribution-based clustering, this study divides all factors into 20 discrete attribute intervals using an equal interval method [23]. Subsequently, the normalized Certainty Factor for each attribute interval is calculated as the contribution index, with the contribution degree being precise to six decimal places. The grading is further refined by clustering and merging equal indices, after which the contribution degrees for the merged attribute intervals are recalculated. The final results provide the data foundation for constructing the sensitivity evaluation index system (Figure 4). Discrete evaluation factors are graded according to their inherent natural attributes. Some results of grading for both continuous and discrete factors are shown in

Table 1. Certainty Factor Values for Some Landslide Susceptibility Assessment Indicators ¹.

Evaluation Factor	Type	Classification Range	CF	Evaluation Factor	Type	Classification Range	CF
Slope (°)	Continuous	0.1–4.1	0.01709	Distance to Fault (m)	Continuous	1–539	0.10062
		4.2–8.2	0.05351			540–1078	0.08177
		8.3–12.3	0.12942			1079–1617	0.06828
		12.4–16.4	0.13847			1618–2155	0.08447
		16.5–20.4	0.11670			2156–2694	0.09169
		20.5–24.5	0.14384			2695–3233	0.02720
		24.6–28.6	0.13701			3234–3772	0.03523
		28.7–32.7	0.08774			3773–4311	0.11018
		32.8–36.8	0.09173			4312–4850	0.03856
		36.9–40.9	0.02880			4851–5389	0.06814
		41–45	0.01777			5390–5927	0.06520
		45.1–49.1	0.03792			5928–6466	0.06236
		49.2–53.1	0			6467–7005	0.03790
		53.2–57.2	0			7006–7544	0.04672
		57.3 –61.3	0			7545–8083	0.03000
		61.4–65.4	0			8084–8622	0
		65.5 –69.5	0			8623–9160	0.05166
		69.6 –73.6	0			9161–9699	0
		73.7–77.7	0			9700–10,238	0
		77.8–81.8	0			10,239–10,777	0
Aspect	Discrete	Flat	0	TPI	Discrete	Flat area	0
		North	0.1145			Lower slope	0.1839
		Northeast	0.1447			Middle slope	0.1936
		East	0.1313			Ridge	0.2516
		Southeast	0.1432			Upper slope	0.3069
		South	0.1276			Valley	0.0641
		Southwest	0.1197
		West	0.1140

¹ Continuous factors are exemplified by slope and distance to fault; discrete factors are exemplified by aspect and TPI.

.

3.3. Dataset Construction and Factor Selection

Due to the scarcity of landslide samples in certain regions, large datasets can often introduce biases or overlook critical local patterns. Therefore, we adopt a carefully constructed small-sample dataset approach [24]. By selecting and using real landslide sample data, this method effectively captures the local characteristics while minimizing prediction bias caused by overall data inadequacy. Additionally, this approach enhances the stability and reliability of the model, especially in regions with limited data. Following the dataset construction, we optimize it using a factor importance evaluation, multicollinearity analysis, and correlation analysis. These methods help eliminate less important factors, those with high linear correlation, and multicollinearity, improving the explanatory power and predictive accuracy of the model [25].

3.3.1. Dataset Construction

Typically, landslide susceptibility assessment is a binary classification problem aimed at predicting whether grid cells in a region are prone to landslides. Thus, the dataset input into the machine learning model consists of two types of samples: positive and negative. This study selects 240 landslide sample points from the landslide inventory as positive samples. Using the Improved Target Space Exteriorization Sampling (ITSES) method proposed by Miao Yamin et al. (2016), an equal number of 240 non-landslide points were selected as negative samples [26]. The labels “1” and “0” were used for the positive and negative samples, respectively. Subsequently, the contribution indexes of all evaluation factors are extracted to the positive and negative sample sets corresponding to the geospatial location, as the a priori knowledge that can effectively represent the local features to construct the dataset, and finally the constructed dataset is randomly divided into a training set and a test set according to the ratio of 7:3 to be inputted into the model for predicting landslide susceptibility of all raster cells in the study area [14].

3.3.2. Evaluation Factor Screening

Factor importance evaluates the influence of each evaluation factor on landslide occurrence. In this study, the average impurity reduction in Random Forests is used to calculate factor importance values. By filtering based on feature importance, we can prevent the introduction of excessive noise and redundant features, thereby reducing model complexity and improving interpretability.

Only factors with an importance score greater than 0.01 were included in the model training, while those with an importance score lower than 0.01 were considered redundant and removed. As shown in Figure 5, among the 12 evaluation factors, slope has the highest importance, making it the most significant controlling factor for landslides in Yulong County. In contrast, profile curvature has the lowest importance score of 0.012, indicating it has the least influence on landslides in the study area. However, since its importance score is above the 0.01 threshold, it was not removed in the importance analysis [27].

A multicollinearity analysis evaluates whether there is excessively high multicollinearity among factors, as high collinearity can reduce the model’s generalization ability due to interdependencies between factors (

Table 2. Results of multicollinearity analysis.

Evaluation Factor	Multicollinearity Analysis
	VIF	TOL
Distance to Road	1.166118	0.857546
Distance to Fault	1.120423	0.892520
Distance to River	1.144622	0.873651
Vegetation Coverage	1.085449	0.921278
Rainfall	1.162597	0.860144
Profile Curvature	1.309268	0.763786
Planar Curvature	1.551096	0.644705
TWI	1.541271	0.648815
SPI	1.831806	0.545989
LR	1.507743	0.663243
Slope	1.749590	0.571563
Aspect	1.095354	0.912947
Lithology	1.217160	0.821585
TPI	1.804903	0.554046

). In this study, multicollinearity of the evaluation factors is quantified using the variance inflation factor (VIF) and tolerance (TOL), where tolerance is the reciprocal of the variance inflation factor. When VIF is less than 2 and TOL is close to 1, it indicates low collinearity among the data. As shown in Table 2, the VIF and TOL values for all 11 evaluation factors fall within the range of low collinearity; thus, no factor elimination is necessary [28].

A correlation analysis can identify factors that provide similar information through the strength of correlations between pairs of factors, allowing for the selection of less correlated factors for classification and predictive modeling. According to the widely accepted Pearson correlation coefficient classification standards proposed by Cohen, factors with a correlation of less than 0.3 are considered to have a weak correlation, those with correlations between 0.3 and 0.5 are considered to have a moderate correlation, and factors with correlations greater than 0.5 are considered to have a strong correlation [29]. Therefore, this study selects factors with a correlation less than 0.3 as the most representative factors for classification and predictive modeling. Pearson correlation coefficients were calculated for each evaluation factor, resulting in the Pearson correlation coefficient chart (Figure 6) [30]. As shown in Figure 6, TPI has correlations of −0.4 with planar curvature and 0.52 with profile curvature, LR has a correlation of 0.54 with slope, and SPI has correlations of 0.56 with TWI and 0.39 with slope. Consequently, the relatively less commonly used factors TPI, LR, and SPI are excluded from the combinations of factors with correlations greater than 0.3.

Synthesizing the results of the three factor selection methods, this study ultimately selects 11 highly representative and independent evaluation factors from the initial 14: the slope, aspect, planar curvature, profile curvature, distance to road, distance to river, distance to fault, rainfall, TWI, vegetation cover, and lithology. These factors are utilized for susceptibility classification and predictive modeling.

3.4. Evaluation Model Construction and Result Partitioning Method

The dataset constructed with positive and negative samples, containing the contribution indices of the selected evaluation factors, serves as the input for the model algorithm, with the probability of landslide occurrence as the output. Susceptibility calculations are performed using the following algorithm model, and the model’s predicted results are used to create susceptibility maps.

3.4.1. Evaluation Model Construction

Sparrow Search Algorithm

The Sparrow Search Algorithm (SSA), proposed by Xue et al. in 2020, is a swarm intelligence optimization algorithm known for its strong optimization capabilities and fast convergence speed [31]. In machine learning, many models require the adjustment of hyperparameters—settings that control the model’s learning process (e.g., learning rate, regularization strength, tree depth). Unlike parameters learned from data, hyperparameters must be set externally and optimized for better performance. SSA is often used to assist classification algorithms by optimizing these hyperparameters, tuning the model for higher accuracy.

The main principle involves assigning roles of discoverers, joiners, and scouters in the hyperparameter search space, with individuals within the swarm switching among these roles to search for the optimal settings. Discoverers are responsible for the global search to find positions with better fitness (better-performing hyperparameter combinations), and joiners perform local searches to refine the solutions in these positions. By dynamically adjusting the number of discoverers and joiners, SSA balances exploration and exploitation, preventing the algorithm from getting trapped in local optima while improving solution accuracy. SSA also includes a danger avoidance mechanism: if scouters detect areas with poor fitness values, all roles rapidly relocate to better areas, continuing the search for optimal hyperparameters. This mechanism effectively enhances the algorithm’s robustness and convergence capabilities (Figure 7).

Support Vector Machine Model

The Support Vector Machine (SVM) model, proposed by Cortes and Vapnik in 1995, is a supervised learning algorithm that excels in handling binary classification problems, making it suitable for landslide susceptibility assessment. The algorithm primarily uses kernel functions to define the geometric margin of sample data within the hyperplane, mapping nonlinearly separable data into a high-dimensional space. It then maximizes the decision boundary margin and seeks the optimal separating hyperplane in the constructed high-dimensional feature space, thereby achieving efficient data classification and prediction [32].

Back Propagation Neural Network Model

The Back Propagation Neural Network (BPNN) algorithm, proposed by Rumelhart and McClelland in 1986, is a widely used error backpropagation algorithm for training multilayer feedforward neural network models [33]. It is well suited for solving complex, nonlinearly separable problems. Data and weights are input into hidden layers containing several neurons, where activation functions such as sigmoid compute the output. The discrepancy between the predicted output and the actual value determines whether to adjust the weights and thresholds before re-inputting the data. The algorithm iteratively adjusts the network’s weights and thresholds through forward propagation, error calculation, backpropagation, and parameter updates to minimize the sum of squared errors, thereby training and predicting with the model. The neuron’s expression is as follows:

$$y=\log sig\left(x\right)={\displaystyle \frac{1}{1+e^{t-\alpha x}}}$$

(3)

In Equation (3), $x$ represents the input data for the neuron; $\alpha$ denotes the corresponding input weight; $t$ is the threshold of the neuron; and $y$ is the output of the neuron.

Random Forest Model

The Random Forest (RF) model, introduced by Breiman in 2001, is an ensemble learning method widely applied in landslide susceptibility assessment. The main principle involves constructing multiple decision tree models to make predictions. Random Forest models consist of multiple decision trees that are built by Bagging, which randomly selects samples with replacement to create multiple training subsets. Each tree is trained on a randomly selected subset of data, introducing diversity among the trees. During classification prediction, the Random Forest aggregates the predictions of all decision trees using voting or averaging to produce the final prediction [34].

Stacking Ensemble Coupling Model

The Stacking ensemble learning model, initially proposed by Wolpert in 1992, has evolved into a robust method for combining diverse machine learning models, thanks to Breiman’s modern Stacking framework utilizing internal cross-validation [35]. The core idea behind Stacking is to select multiple independent and diverse machine learning models as base classifiers. The dataset is input into these base classifiers through cross-validation to prevent overfitting. The predictions generated by these base classifiers are then combined and used as the dataset for a second-layer meta-classifier, which learns the optimal strategy for combining the predictions of the base classifiers. This approach leverages the strengths of each base classifier, improving the overall predictive accuracy of the landslide susceptibility model.

In this study, we selected three different coupled models—CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF—as the base classifiers. The rationale behind choosing SVM, BPNN, and RF models for coupling as base models lies in their well-established performance in classification tasks and their complementary strengths, which can enhance the model’s overall performance and robustness [36]. Other models, such as decision trees or K-Nearest Neighbors, were not chosen due to their tendencies towards overfitting or scalability issues in the context of large datasets, which might have negatively impacted the ensemble’s predictive accuracy [37]. A Logistic Regression model (LR) was employed as the meta-classifier to construct a two-layer ensemble coupled learning model (Figure 8).

3.4.2. Landslide Susceptibility Mapping Zones

The susceptibility mapping method primarily categorizes susceptibility levels based on either geometric interval proportions or grid area proportions [14,38]. Building upon previous research, this study divides the predicted results of four coupled models (CF-SSA-SVM, CF-SSA-BPNN, CF-SSA-RF, and CF-SSA-Stacking) according to their respective grid area proportions, ranging from highest to lowest, into five levels: very high susceptibility (10%), high susceptibility (20%), moderate susceptibility (40%), low susceptibility (20%), and very low susceptibility (10%). This approach is used to generate a landslide susceptibility zoning map [38].

3.5. Modeling Evaluation Metrics

For ease of comparing the performance differences in susceptibility prediction among the CF-SSA-Stacking ensemble model and the three basic coupled models, this study employs four metrics as dimensional measures to evaluate model effectiveness: the accuracy, F1 score, Kappa coefficient, and ROC curve. Accuracy refers to the proportion of correctly predicted positive samples. Precision indicates the proportion of correctly predicted positive samples out of all predicted positive samples. Recall is the proportion of correctly predicted positive samples out of all actual positive samples. Higher values of these metrics indicate stronger capability of the model to predict positive samples. The F1 score is the harmonic mean of precision and recall, providing an overall reflection of algorithm performance; higher values indicate better classification performance [8]. An ROC curve closer to the upper-left corner signifies a larger AUC value, indicating better classification predictive performance. The Kappa coefficient measures the consistency between model predictions and actual values in classification accuracy evaluation. A Kappa coefficient in the range of [0.6, 1] indicates high reliability of the model [14]. The computation formulas are shown in Equations (4)–(11):

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

(4)

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

(5)

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

(6)

$$Specificity={\displaystyle \frac{\mathrm{TN}}{FP+TN}}$$

(7)

$$F1={\displaystyle \frac{2\times Precision\times Recall}{Precision+Recall}}$$

(8)

$$Kappa={\displaystyle \frac{Accuracy-p_e}{1-p_e}}$$

(9)

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

(10)

$$\mathrm{AUC}={\displaystyle \sum _{i=1}^{\mathrm{n}-1}\frac{\left(\left(Sensitivit{y}_{i+1}+Sensitivit{y}_i\right)\times \left(\left(1-Specificit{y}_{i+1}\right)-\left(1-Specificit{y}_i\right)\right)\right)}{2}}$$

(11)

In Equations (4) and (11), TP (True Positive) and FN (false negative) represent the number of landslide points in the predicted samples correctly identified as landslides and the number of landslide points incorrectly identified as non-landslides, respectively. FP (false positive) and TN (True Negative) represent the number of non-landslide points incorrectly identified as landslides and the number of non-landslide points in the predicted samples correctly identified as non-landslides, respectively.

3.6. SHAP Value

The SHAP (Shapley Additive Explanations) method, proposed by Lundberg and Lee in 2017, is based on Shapley values from game theory [39]. It is designed to explain the predictions of complex machine learning models, particularly those with “black-box” characteristics. SHAP quantifies the contribution of each feature to the model’s output in a way that does not depend on the internal structure of the model. The core idea is to compute the marginal contribution of a single feature across all possible subsets of features, known as Shapley values, to provide both local and global explanations of the model. Due to its powerful data visualization capabilities, SHAP has been widely used to enhance model transparency and interpretability [40]. The formula for Shapley values is provided as follows (12). In this paper, we utilize the SHAP method to visualize the decision-making process of the model, thereby improving its transparency and interpretability.

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

(12)

In Equation (12), the term ${\phi }_i$ on the left-hand side represents the Shapley value of feature $i$, which quantifies the contribution of feature $i$ to the prediction output.$\mathrm{N}$ represents the set of all features in the mode. Meanwhile, $S\subseteq N\backslash \left\{i\right\}$ denotes any subset of $\mathrm{N}$ that excludes the feature $i$. The term $\left|S\right|$ represents the number of features in subset $S$. $f\left(S{\displaystyle \cup \left\{i\right\}}\right)$ describes the prediction made when the feature $i$ is added to the subset $S$. ${\displaystyle \frac{\left|S\right|!\left(\left|N\right|-\left|\mathrm{S}\right|-1\right)!}{\left|\mathrm{N}\right|!}}$ serves as a weighting factor. This factor accounts for the different possible orderings of feature subsets, ensuring that the contribution of feature $i$ is fairly averaged across all possible combinations of feature subsets.

4. Results and Analysis

4.1. Spatial Consistency Analysis of Landslide Susceptibility Mapping Results

Figure 9 presents the susceptibility mapping results of the four coupled models. To visually compare the susceptibility mapping results of various models, this paper conducts a comparative analysis of the landslide susceptibility prediction performance by examining the ratio of raster cell areas in each susceptibility class (A) to the number of landslide points in each class (B), as well as the landslide frequency ratio (B/A) [32]. As can be seen from

Table 3. Historical Landslide Distribution Statistics for Each Susceptibility Mapping.

Model		CF-SSA-SVM	CF-SSA-BPNN	CF-SSA-RF	CF-SSA-Stacking
Very low	A/%	10	10	10	10
	B/%	0.42	0	0.42	0
	B/A	0.04	0	0.04	0
Low	A/%	20	20	20	20
	B/%	3.33	3.33	3.75	3.75
	B/A	0.17	0.17	0.19	0.19
Moderate	A/%	40	40	40	40
	B/%	20.83	22.08	20.83	20.42
	B/A	0.52	0.55	0.52	0.51
High	A/%	20	20	20	20
	B/%	30.42	27.50	27.92	24.17
	B/A	1.52	1.38	1.39	1.21
Very high	A/%	10	10	10	10
	B/%	45.00	47.08	47.08	51.67
	B/A	4.50	4.71	4.72	5.17

, the CF-SSA-Stacking model exhibits the highest landslide frequency ratio (2.528) and the highest proportion of landslide points (75.84%) in the high/very high susceptibility grades compared to other models. Simultaneously, it shows the lowest proportion of landslide points (3.75%) and the lowest landslide frequency ratio (0.125) in the very low/low susceptibility grades. The results indicate that the CF-SSA-Stacking model outperforms the other three models in classifying historical landslides into very-low-, low-, high-, and very-high-susceptibility categories, signifying that its predictive capability is the best among the four well-constructed models [3].

4.2. Modeling Evaluation Metrics Results

Table 4. Accuracy Statistics of Evaluation Metrics.

Model	Accuracy	F1	Kappa Coefficient
CF-SSA-SVM	0.889583	0.887473	0.779167
CF-SSA-BPNN	0.870833	0.870293	0.741667
CF-SSA-RF	0.877083	0.876310	0.754167
CF-SSA-Stacking	0.893750	0.891720	0.787500

and Figure 10, respectively, show the results of the three evaluation metrics and the AUC (Area Under the ROC Curve) values. Generally, if a model’s evaluation metrics show accuracy and F1 values greater than 0.8, a Kappa coefficient greater than 0.6, and a larger AUC value, it indicates high accuracy of the model [23]. The results of accuracy and F1 for the four coupled models all exceed 0.85, and the Kappa coefficients are all above 0.7 and the AUC values of all four coupled models surpass 0.9. This implies that all validation models have met or exceeded the predetermined minimum accuracy criteria, ensuring the reliability and validity of the models.

Among them, the performance of the CF-SSA-Stacking model is particularly outstanding. Compared with the other three base coupling models, the CF-SSA-Stacking model shows improvements in the accuracy, F1 score, Kappa coefficient, and AUC by approximately 2.3%, 2.1%, 4.6%, and 1.1%, respectively. These results indicate that the CF-SSA-Stacking model exhibits significant advantages across all evaluation metrics, making it the most accurate model.

However, although the CF-SSA-Stacking model improves the accuracy of landslide susceptibility assessment, it also introduces greater computational complexity due to the involvement of multiple base models and a meta-learning layer [35]. For example, the CF-SSA-BPNN base model has a lengthy processing time, which increases the overall computational complexity of the CF-SSA-Stacking model, thereby affecting the efficiency of the final prediction results [41]. Consequently, in resource-constrained environments (such as real-time monitoring systems or devices with limited processing capabilities), this complexity may limit its practical application.

Therefore, a balance must be struck between model performance and computational efficiency. Future research could focus on optimizing the base models to reduce the computational burden without significantly decreasing prediction accuracy. Possible optimization methods include using simpler and more efficient base models or improving feature selection criteria, thereby reducing the overall computational complexity and enhancing the model’s feasibility in practical applications.

4.3. Distribution Patterns of Landslide Susceptibility Index

To analyze the differences in the distribution patterns of susceptibility indices among different models, histograms, box plots, and probability density plots of the susceptibility index distribution were further produced (Figure 11). From Figure 11a–d, it can be seen that the distribution patterns of the susceptibility indices differ significantly. The CF-SSA-RF model exhibits a near-normal distribution trend with “high in the middle and low on both ends.” The other models exhibit a bimodal distribution trend with “low in the middle and high on both ends”, with the low-susceptibility area on the left peak being significantly higher than the high-susceptibility area on the right peak. This aligns with the actual situation in the study area, where most landslide developments are in low- and very-low-susceptibility zones [42].

Generally, models with smaller mean susceptibility indices and larger standard deviations have lower uncertainty, and the distribution pattern of susceptibility indices is more reasonable [43,44]. From Figure 11, it is also apparent that the order of mean values of the prediction results of the four coupled models is CF-SSA-BPNN < CF-SSA-Stacking < CF-SSA-SVM < CF-SSA-RF; and the order of standard deviations is CF-SSA-BPNN > CF-SSA-Stacking > CF-SSA-SVM > CF-SSA-RF. It can be observed that the CF-SSA-BPNN model has the smallest mean and the largest standard deviation. This indicates that the average level of its predictions is the lowest and more dispersed compared to other models, meaning the prediction results are closer to the high-value and low-value regions at both ends. However, the mean and standard deviation of the CF-SSA-BPNN model do not align with the order of model accuracy evaluation results, which is contrary to the usual pattern observed in previous studies. This discrepancy arises because the predictions of the CF-SSA-BPNN model are more conservative compared to other models, tending towards lower values rather than higher ones. This leads to a significant increase in the number of false negatives.

From the box plot in Figure 11e, it is observed that the Interquartile Range (IQR) of the CF-SSA-Stacking model is greater than that of the CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF models. This indicates that the prediction results of CF-SSA-Stacking are more dispersed within the middle 50% of the data distribution compared to the other three models, suggesting a stronger capability in handling diverse types of data. Additionally, it can be noted that the median of the CF-SSA-Stacking model, along with CF-SSA-SVM and CF-SSA-BPNN models, is less than 0.5, and the overall IQR is closer to low susceptibility values. This tendency indicates that these three models tend to output lower predicted values, thereby reducing the occurrence of false positives but potentially increasing false negatives, similar to the analysis of the CF-SSA-BPNN model mentioned earlier [45].

From the probability density plots in Figure 11f, it is evident that the curve of the CF-SSA-Stacking model’s predictions is smoother compared to the results of the other three models. This smoothness indicates that the CF-SSA-Stacking model, by integrating the advantages and disadvantages of three base models through ensemble learning, achieves more stable predictions and is less susceptible to abrupt changes caused by local noise compared to the other three models. Moreover, from Figure 11e,f, it can also be seen that there are differences in the susceptibility index ranges of the predicted results among the models. Both the CF-SSA-RF and CF-SSA-Stacking ensemble models have ranges smaller than the interval (0, 1). Additionally, it is evident that the CF-SSA-Stacking model’s predicted probability range lies between the ranges of the other three models. This indicates that the prediction range of the CF-SSA-Stacking model is relatively narrow, suggesting that it is less likely to generate extreme values compared to other models and reduces the impact of outliers.

5. Discussion

5.1. The Advantages of the Proposed Method over Traditional Methods

This study employs the coefficient of determination as a measure of contribution, using contribution-based clustering to classify continuous factor attributes. The Sparrow Search Algorithm (SSA) is then applied to optimize parameters for constructing a coupled ensemble model for landslide susceptibility evaluation. Experimental results demonstrate that all four models achieved satisfactory performance, validating the feasibility of the entire methodology on the current dataset.

Moreover, the optimized CF-SSA-Stacking coupled ensemble model outperformed the CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF base coupled models across multiple performance metrics. Specifically, the CF-SSA-Stacking model exhibited improvements of approximately 2.3%, 2.1%, 4.6%, and 1.1% in the accuracy, F1 score, Kappa coefficient, and AUC, respectively. This indicates that the Stacking model significantly enhances overall performance by integrating the strengths of multiple base models.

Compared to single models, the Stacking method can more effectively combine the advantages of each model, thereby mitigating potential issues of underfitting or overfitting. For instance, SVM might generate unstable decision boundaries in high-dimensional data, BPNN neural networks could suffer from inadequate training due to gradient vanishing, and while RF performs well in certain scenarios, it may exhibit bias when dealing with highly correlated factors. By integrating these models through Stacking, these issues can be effectively addressed, thereby enhancing the model’s overall generalization ability.

Additionally, when predicting with the trained model, the CF-SSA-Stacking model achieved the highest landslide frequency ratio (2.528) and landslide count ratio (75.84%), further demonstrating its robust capability in handling unseen data. This implies that the model not only performs well on known data but also maintains high predictive accuracy when applied to unseen data, indicating strong generalization ability.

In conclusion, the CF-SSA-Stacking coupled ensemble model exhibits significant advantages in landslide susceptibility evaluation. Its improvements across multiple performance metrics and its exceptional performance on unseen data make it an ideal choice for analyzing the current dataset. In future applications, this model is expected to provide more accurate and reliable support for landslide prediction and geological disaster prevention.

5.2. Distribution Patterns of Landslide Susceptibility Index Analysis

The CF-SSA-Stacking and CF-SSA-RF models show a prediction probability range smaller than the interval [0, 1]. Single models like SVM and BPNN usually directly output probability predictions, typically covering the entire [0, 1] interval. In contrast, the prediction results of ensemble models like Stacking and RF are obtained by combining the predictions of multiple base models through averaging and integration methods. This secondary learning process significantly reduces the occurrence of extreme values, resulting in a more concentrated and stable prediction probability range for ensemble models [46,47]. Additionally, the Stacking algorithm can identify appropriate base model fusion rules through a meta-learner, integrating the strengths and balancing the prediction errors of each base model. This characteristic leads to a similarity between the susceptibility index distribution patterns of the CF-SSA-Stacking model and the base models’ prediction results, explaining why the prediction probability range of the CF-SSA-Stacking model falls between that of the CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF models. Previous studies have indicated that generally, higher model accuracy is associated with smaller mean prediction values and larger standard deviations. In this study, although the CF-SSA-Stacking ensemble model achieves the highest accuracy, its susceptibility index distribution pattern is not as ideal as that of the CF-SSA-BPNN model. The analysis in the previous section shows that the CF-SSA-BPNN model’s prediction results are predominantly in the low-value area, resulting in the smallest mean and largest standard deviation. This leads to a higher number of false negatives compared to other models, contributing to the overall lower accuracy of the CF-SSA-BPNN model [43].

Overall, the mean and standard deviation of prediction results reflect the statistical characteristics of the prediction distribution, aiding in understanding the stability and reliability of the predictions. However, they do not directly indicate the model’s performance. Model performance is typically assessed using metrics such as the accuracy, F1 score, Kappa coefficient, and AUC value [48]. Therefore, despite the less-than-ideal mean and standard deviation of the CF-SSA-Stacking model’s prediction results, the model achieves the highest prediction accuracy due to its ability to integrate the strengths of multiple models and mitigate the weaknesses of individual models. This demonstrates its superiority in susceptibility prediction performance compared to the three basic coupled models.

5.3. SHAP-Based Interpretability Analysis of Model Results

In this study, the SHAP (Shapley Additive Explanations) analysis method is employed to comprehensively evaluate the impact of each feature factor on the final landslide susceptibility prediction results and the marginal contribution of each base model to the ensemble model. This approach not only elucidates how the model leverages input features to make decisions but also clarifies the importance of each base model within the ensemble. This is valuable for validating whether the model’s predictions align with the characteristics of known landslide events and for analyzing the extent of participation of each base model in the ensemble. The application of SHAP in this context helps to address the “black-box” issue in ensemble models, providing a more interpretable and trustworthy framework for practitioners.

Figure 12 presents the SHAP heatmaps of four coupled models: CF-SSA-SVM, CF-SSA-BPNN, CF-SSA-RF, and CF-SSA-Stacking. The left y-axis shows the feature importance rankings determined by SHAP values, ordered from highest to lowest contributions. The right side provides the visual representation of these rankings. Each column represents an instance from the test set, where the color intensity indicates the absolute magnitude of the SHAP values—the darker the color, the more significant the feature’s impact on the model’s prediction. The top part visualizes the models’ prediction results under different feature values [49]. From Figure 12, it can be observed that all four models rank the top four features—the slope, planar curvature, distance to river, and distance to road—consistently. This consistency indicates that these models effectively identify key features closely related to geological processes in landslide prediction. Specifically, the slope, as a topographic variable, is a core feature for the models due to its direct influence on landslide occurrences. Additionally, features like the planar curvature, distance to river, and distance to road further assist the models in understanding the terrain and hydrological conditions, thereby enhancing their capability to assess landslide risk. However, beyond these key features, the importance rankings of other features vary between models. This suggests that while the models exhibit consistency in identifying the major risk factors, differences in model structure and learning mechanisms lead to discrepancies in how they evaluate the contribution of secondary features. These variations reflect the models’ distinct preferences in feature selection and sensitivity, highlighting the significant influence of model architecture on feature extraction.

Figure 13 presents the SHAP dependence plot for the slope feature, illustrating its contribution to the model’s output. The x-axis represents the normalized Certainty Factor value, which was previously calculated and reflects the certainty of landslide occurrence, while the y-axis indicates the influence of the slope feature on the model’s landslide susceptibility prediction (y values greater than 0 indicate a positive contribution, while y values less than 0 indicate a negative contribution) [50]. From Figure 13, it can be observed that the contribution trends of the slope feature are similar across all four models, with most of the positively contributing sample points concentrated in the region where the normalized Certainty Factor is greater than 0.08. Based on the analysis of Table 1, the slope range of 8.3–36.8° aligns closely with the slope range of landslide-prone areas (10–45°) reported in the literature [51]. This indicates that the model effectively captures the relevant geological characteristics when predicting landslide risks. Notably, the Stacking model, despite its more complex structure, continues to produce results that are consistent with geological principles. This demonstrates that the model’s complexity does not compromise its reasonable assessment of the slope feature. Therefore, it can be further concluded that the model exhibits strong interpretability and robustness in predicting landslide risk.

Figure 14 illustrates the marginal contributions of each base model to the ensemble model (Stacking), calculated through the SHAP analysis using the mean absolute SHAP values [39]. As shown in Figure 14, CF-SSA-SVM has the highest mean absolute SHAP value (0.7968), indicating that it plays the most critical role in the Stacking model and contributes significantly to the final prediction outcome. In comparison, CF-SSA-RF and CF-SSA-BPNN contribute less, with SHAP values of 0.6851 and 0.6717, respectively, but they still play an important role in capturing different patterns and features within the data. The differences in marginal contributions among the base models suggest that they effectively complement each other in the Stacking structure, compensating for the weaknesses of the other models and significantly enhancing the overall model performance. Although CF-SSA-SVM has the highest marginal contribution, CF-SSA-RF and CF-SSA-BPNN also contribute to improving the model’s robustness and generalization capability. This further demonstrates that, despite the complexity of the ensemble model, the complementary nature of the base models enhances the overall predictive performance and generalizability of the model.

In summary, by combining the SHAP value analysis of feature importance and the marginal contribution analysis of the base models, we not only improve the interpretability of the model but also provide practitioners with more intuitive and transparent decision-making insights for real-world applications. This approach effectively mitigates the “black-box” effect of ensemble models, making the prediction process more transparent and understandable. These analyses demonstrate how the model synthesizes information from both features and base models to enhance overall predictive performance and reliability.

5.4. Implications, Limitations, and Future Directions

Although the model employed in this study demonstrated strong performance in the Yulong County region, challenges may arise when extending its application to other geographic areas. This is primarily due to significant variations in environmental characteristics across different regions, particularly in terms of topography, vegetation cover, and soil stability [52]. Specifically, the model may require adjustments in input data types, feature weights, and evaluation metrics when applied to different geomorphological types. For example, in complex mountainous canyon regions, where the terrain is steep, vegetation cover is sparse, and soil erosion risk is high; these factors should be given higher weight in the model [53]. Conversely, in regions like the Loess Plateau or low mountainous hills, soil properties and rainfall erosion may dominate, necessitating adjustments to the input variables and weighting scheme. Therefore, the model developed for the complex mountainous canyon landscape in this study should primarily be applied to regions with similar geological and geomorphological characteristics. However, for areas with distinctly different geomorphological features, such as the Loess Plateau or low mountainous hills, new evaluation metrics may need to be introduced, and the model recalibrated to ensure its applicability and predictive accuracy in these new geographic contexts. This implies that future research should focus on optimizing and adjusting the model according to the specific conditions of different regions to ensure its generalizability and robustness.

In this study, most of the evaluation factor data were derived from publicly available high-resolution datasets, and such data were accessible for most of the study areas. However, in some regions where high-resolution data are unavailable, it may be necessary to reduce the number of features or use alternative datasets that are more readily accessible. This, however, may affect the model’s predictive accuracy, necessitating retesting and fine-tuning the model when alternative data are used to ensure it adapts to new data characteristics. Additionally, the landslide sample points used in this study are primarily based on the results of field geological disaster investigations, providing a reliable foundation for model training and feature learning. However, the predictive performance of the model is somewhat constrained by the geographic coverage of the current sample dataset. The limited distribution of samples may result in reduced predictive accuracy in areas that are not well represented by the dataset [10]. To further improve the model’s predictive accuracy, future research should focus on enriching the dataset and enhancing the model’s robustness by incorporating multi-source data. The integration of remote sensing data, meteorological data, and soil monitoring data, for instance, can not only increase the diversity of positive and negative samples but also effectively improve the model’s generalization capability, particularly in complex terrain conditions. Existing studies have demonstrated that high-quality multi-source data significantly enhance the model’s predictive performance across different geographic environments.

In summary, the future development of this study should emphasize three main aspects: first, adjusting the input data types and model weights for different geographic environments to ensure model adaptability; second, exploring the use of lower-resolution or alternative data in areas with data scarcity and adjusting the model accordingly to accommodate new data characteristics; and finally, focusing on improving the model’s overall predictive performance and applicability in broader geographic areas by integrating multi-source data. Future research should validate the effectiveness of these strategies to ensure the model’s generalizability and practicality.

6. Conclusions

This study uses Yulong County as a case example, employing the Contribution Degree Clustering Method and the natural attribute interval classification method to grade continuous and discrete factors, respectively. Subsequently, four coupled ensemble models—CF-SSA-SVM, CF-SSA-BPNN, CF-SSA-RF, and CF-SSA-Stacking—were chosen to evaluate landslide susceptibility assessment. A comparative analysis of the index distribution patterns of the susceptibility prediction results and the accuracy of each model revealed that using the Contribution Degree Clustering Method with the determinacy coefficient as the contribution index for grading the attribute intervals of continuous factors yields predictive results consistent with actual conditions, demonstrating high applicability. Furthermore, the accuracy evaluation results of the CF-SSA-Stacking ensemble coupled model show the highest values of the accuracy, F1, Kappa coefficient, and AUC at 0.89375, 0.89172, 0.7875, and 0.9522, respectively, significantly outperforming the CF-SSA-SVM, CF-SSA-BPNN, and CF-SSA-RF basic coupled models. This indicates that the Stacking model achieves better generalization ability and reliability compared to the basic coupled models. Additionally, the landslide frequency ratio (2.528) and the landslide number ratio (75.84%) in the susceptibility mapping results of the CF-SSA-Stacking ensemble coupled model are the highest, suggesting that this model can identify the most historical landslides with the least grid area of very-high-/high-susceptibility zones compared to the basic coupled models.

Overall, the CF-SSA-Stacking ensemble coupled model proves to be an effective method for regional landslide susceptibility assessment, particularly in practical disaster prevention and mitigation applications in complex mountainous terrains and canyon landscapes. However, the model’s applicability to different geographical regions may be affected by variations in environmental characteristics, especially in terms of topography, vegetation cover, and soil stability. Therefore, future research should focus on adjusting the model’s input data types, feature weights, and evaluation criteria according to specific geomorphological conditions to ensure its broader applicability and prediction accuracy. Additionally, by integrating multi-source data, such as remote sensing, meteorological, and soil monitoring data, the model’s generalization capability and robustness could be further enhanced, providing more reliable support for landslide disaster prevention and mitigation in various geographical settings.