An interpretable model for landslide susceptibility assessment based on Optuna hyperparameter optimization and Random Forest

Abstract

This study proposed an interpretable model that combines Random Forest (RF), Optuna hyperparameter optimization, and SHapley Additive exPlanations (SHAP) to achieve optimal landslide susceptibility evaluation and provide explanations in the northwest region of Yunnan Province in China. First, an inventory of 4447 landslides and 23 related factors was considered for the landslide susceptibility assessment. Subsequently, a hyperparameter-optimized RF model was developed using the Optuna framework and the training dataset to generate landslide susceptibility maps. The performance of the models were evaluated using accuracy (ACC), precision (PPV), recall (TPR), F1-score (F1), and the Area Under the Curve (AUC) based on the Receiver Operating Characteristic. Furthermore, the interpretability of the model was enhanced through the implementation of SHAP. The proposed model demonstrated outstanding performance on the test set, achieving an ACC of 0.7792, PPV of 0.7448, TPR of 0.8769, F1 of 0.8055, and an AUC of 0.8387. The interpretability analysis revealed that elevation, population density, distance from roads, and normalized difference vegetation index were the primary factors influencing landslide occurrences in the study area. This study provides a comprehensive framework for evaluating landslide susceptibility in specific regions and offers invaluable insights for the prevention and management of landslide disasters.

Keywords:

• Landslide susceptibility
• Random Forest
• hyperparameter optimization
• explainability
• influence factors

1. Introduction

Landslides, a prevalent type of geological disaster, are characterized by unpredictability, extensive distribution, and destructive power, often posing significant threats and causing substantial losses to human lives and the natural environment (Xu et al. 2013; Froude and Petley 2018). Between 1998 and 2017, 378 major landslides occurred globally, resulting in the deaths of 18,414 individuals, injuries to 4.8 million people, and approximately 8 billion dollars in economic damages (Liao et al. 2022). Therefore, the precise assessment of regional landslide susceptibility is of critical importance for mitigating disaster risks, safeguarding public safety and property, and guiding the rational utilization of regional land resources.

In recent years, landslide susceptibility research has emerged as a focal point of the study, being widely recognized as one of the most effective methods to mitigate landslide disasters (Saha et al. 2021). Various methodologies have been employed for landslide susceptibility analysis, broadly categorized as follows: (1) traditional statistical methods, including frequency ratio model (Li et al. 2016; Zhang et al. 2020), weight of evidence (Ilia and Tsangaratos 2016), and analytic hierarchy process (Wei et al. 2021); (2) index-based methods, such as entropy index and statistical index (Wang et al. 2016; Wu et al. 2016; Shirani et al. 2018); (3) machine learning models, including Random Forest (RF) (Stumpf and Kerle 2011), Boosted Regression Trees (Saha et al. 2021), Support Vector Machine (SVM) (Chen et al. 2022), Logistic Regression (Zhao et al. 2019), and eXtreme Gradient Boosting (XGBoost) (Zhou et al. 2022). Although numerous approaches have demonstrated viability in assessing landslide susceptibility, the optimal method remains undetermined. Machine learning techniques have exhibited distinct advantages when dealing with large datasets and nonlinear relationship (Prodhan et al. 2022), and has gained widespread application in landslide susceptibility assessment (Liu et al. 2023; Long et al. 2023; Ma et al. 2023). In comparative studies of various machine learning methods, RF has found widespread adoption in the field of geological hazard susceptibility assessment. This is due to its strengths in handling high-dimensional data, evaluating variable importance, and its accuracy and stability in model construction (Trigila et al. 2015; Hong et al. 2017; Wang et al. 2020).

The accuracy of machine learning models is not solely determined by the algorithms themselves but is also significantly influenced by a predetermined set of hyperparameters before model training (Wang et al. 2021). Complex machine learning models often involve a multitude of hyperparameters, and optimizing them is a critical step in the model-tuning process. However, due to the laborious nature of hyperparameter optimization during model development, researchers consider various alternatives to manual tuning. Techniques such as random search, grid search, and genetic algorithms have been employed for hyperparameter optimization (Li and Kong 2014; Ribeiro da Silva et al. 2019; Tang et al. 2019). Nevertheless, random and grid search suffer from low optimization efficiency and high computational costs, while genetic algorithms tend to converge to local optima, potentially overlooking global performance. In recent years, traditional Bayesian algorithms have been increasingly utilized for hyperparameter optimization in machine learning models. For example, Xie et al. (2021) employed Bayesian optimization with SVM to assess landslide hazards in Nanping City, China. Sun et al. (2020) applied Bayesian optimization to RF model for assessing landslide susceptibility in Fengjie, China. However, traditional Bayesian algorithms, which employ Gaussian processes for the optimization process of their probabilistic models, are computationally intensive and better suited to continuous hyperparameters, making them less efficient for discrete ones. Optuna, a novel hyperparameter optimization framework, utilizes Tree-structured Parzen Estimator (TPE) approaches to model the objective function, providing more flexibility and efficiency than traditional Gaussian Process Bayesian optimization (Akiba et al. 2019). Its suitability for various types of hyperparameters (Lin et al. 2023). However, limited research has been conducted on the application of Optuna in the field of landslide susceptibility to date.

While machine learning models have demonstrated substantial advantages in handling complex data and exhibited commendable performance, the non-interpretability of the majority of these models remains a challenging issue. This lack of interpretability hinders a comprehensive assessment or trust in these models based on a single metric (Chen and Fan 2023a; Youssef et al. 2023). In recent times, the advent of a post-hoc explanatory method called SHapley Additive exPlanations (SHAP) has presented a novel approach to interpreting the outcomes of landslide susceptibility assessments. By iteratively hypothesizing and testing, SHAP illuminates how data influences model functioning, attracting the attention of researchers from various disciplines. For instance, Guliyev and Mustafayev (2022) utilized SHAP in conjunction with RF and XGBoost to interpret factors influencing the dynamics of crude oil price. Abdollahi and Pradhan (2023) utilized SHAP to interpret the influencing factors in wildfire susceptibility prediction models. Dikshit and Pradhan (2021) applied SHAP to investigate how deep-learning models arrive at specific outcomes under various drought conditions. In summary, SHAP has the potential to clarify model predictions and provide transparency to ‘black-box’ models; however, research specifically focusing on the application of SHAP in the domain of landslide susceptibility is still limited.

To tackle the aforementioned issues, an innovative interpretable model is proposed in this study by combining RF, Optuna hyperparameter optimization, and SHAP. This amalgamation aims to achieve an optimal assessment of landslide susceptibility while providing transparent explanations. The main objectives of this study are twofold: first, to develop a hyperparameter-optimized RF model utilizing the Optuna framework for the identification of landslide susceptibility; and second, to employ the SHAP method to interpret the model and further evaluate the factors influencing landslide susceptibility. This study will bridge the existing gap in machine learning for landslide susceptibility research and provide a more reliable foundation for practical applications in related fields.

2. Study area and data

2.1. Study area

Yunnan Province is situated on the edge of southwestern China, encompassing approximately 394,100 km² of land area. Its topography is characterized by highlands and mountains, featuring steep peaks and deep river valleys in the western regions. Influenced by a subtropical monsoon climate, Yunnan experiences distinct dry and wet seasons, with abundant and concentrated precipitation during May to October (Wang and Yu 2021). These complex topographical conditions and diverse climatic patterns render Yunnan as one of the provinces in China most susceptible to landslides. Notably, the northwestern part of Yunnan (23°77′ − 29°15′N, 97°31′ − 99°83′E, covering an area of about 63,478 km²) is marked by densely distributed mountain ranges and valleys (Figure 1) (Li et al. 2015). The area is traversed by the Lancang river and the Nujiang river, both of which are characterized as typical mountain rivers known for their rapid flow and significant drop (Liu et al. 2023). The river valley slopes on either side are relatively steep (He et al. 2022). The turbulent runoff within the area, coupled with concentrated rainfall, provides favourable hydrodynamic conditions for the frequent occurrence of landslides. Combined with the formidable topography, these factors make the area a focal point for landslide prevention and control efforts.

Figure 1. Location of the study area (a) and kernel density map of landslides (b).

2.2. Data

2.2.1. Landslide inventory

In this study, a landslide inventory was compiled using survey database, consisting of 4447 landslides that occurred in the northwestern part of Yunnan Province. These landslides were derived from spatial distribution data of geological hazard points in China, provided by the Resources and Environmental Science Data Center (RESDC) (http://www.resdc.cn), as well as detailed survey database on geological landslide disasters in Yunnan Province. The kernel density map of landslide events is presented in Figure (1b). The density of landslides reflects the number of occurrences within a unit area. Overall, landslide events are primarily concentrated in the southern part of the study area, with a maximum density of 0.4442 events/km². Fewer events are distributed in the northern part, with most clusters observed around river valleys.

2.2.2. Conditioning factors

This study, building upon previous studies and the characteristics of landslide development in the study area, collected foundational data for landslide factors from three aspects: precipitation, surface environment, and human activities. Detailed data used in this study can be found in Table 1. Daily precipitation data from 1986 to 2018 were obtained from the Yunnan Provincial Climate Center, comprising records from 125 meteorological stations. For the surface environment category, the digital elevation model (DEM) was obtained from the Shuttle Radar Topography Mission (SRTM) (Rodríguez et al. 2006). Soil types (Soiltype) and land use/land cover (LULC) data were obtained from RESDC. Geological lithology data with a spatial resolution of 250 m were sourced from Sayre et al. (2014). To calculate the Normalized Difference Vegetation Index (NDVI), Landsat 5/8 satellite imagery spanning from 1986 to 2020 was retrieved from the Google Earth Engine platform (https://earthengine.google.com). River data were obtained from the National Catalogue Service for Geographic Information (https://www.webmap.cn/), while fault zone data were sourced from the Comprehensive karst hydrogeological map of Yunnan Province. Within the human activities category, population density (POP) data were selected from the widely used WorldPop dataset (Tatem 2017) (https://hub.worldpop.org/). Roads data, available in shapefile format, were obtained from the National Catalogue Service for Geographic Information.

Table 1. Detailed information of the data.

Data category	Data name	Time	Spatial resolution	Data source
Precipitation	Daily precipitation	1986–2018	–	Yunnan Province Climate Center
Surface environment	Elevation	2009	30 m	SRTMv3.0
	Lithology	2012	250 m	Sayre et al. (2014)
	Soil types	1995	1 km	RESDC
	NDVI	1986–2020	30 m	Landsat 5/8
	LULC	2005	1 km	RESDC
	Rivers	2015	–	National Catalogue Service for Geographic Information
	Fault zone	–	–	Comprehensive karst hydrogeological map of Yunnan Province
Human activities	Population density	2000–2020	1 km	Tatem (2017)
	Roads	2015	–	National Catalogue Service for Geographic Information

3. Methodology

Figure 2 depicts the methodological roadmap employed in this study for the assessment of landslide susceptibility. The assessment process comprises multiple steps. First, an extensive inventory consisting of 4447 landslides and 23 associated factors was prepared for the evaluation of landslide susceptibility. Next, three models were developed: RF, Light Gradient Boosting Machine (LGBM), and Long Short-Term Memory (LSTM), each optimized with hyperparameters using the Optuna framework. The models’ performance was assessed using accuracy (ACC), precision (PPV), recall (TPR), F1-score (F1), and Area Under the Curve (AUC) based on the Receiver Operating Characteristic (ROC). Additionally, RF model utilized a training dataset to generate maps depicting the susceptibility of landslides. Finally, the interpretability of RF model was enhanced through the implementation of SHAP.

Figure 2. Methodological framework of this study.

3.1. Selection and data processing of landslide influencing factors

3.1.1. Precipitation factors

Precipitation is a particularly crucial factor in the analysis of landslide susceptibility due to its significant impact on slope stability and soil cohesion. In this study, we selected nine key precipitation index (Table 2) as factors that influence landslides (Zhang et al. 2011). Indices such as annual total wet-day precipitation (PRCPTOT) and simple daily intensity index (SDII) were specifically chosen to reflect the spatial distribution characteristics of precipitation. To accurately compute these precipitation indices, we employed the RClimDex software. Subsequently, we conducted spatial interpolation of the nine precipitation indices, collected from 125 meteorological sites, using the ordinary Kriging interpolation method in ArcGIS 10.4, at a spatial resolution of 1 km (Figure 3). This step aimed to generate raster data for precipitation factors, thereby providing spatial continuity and detailed richness for our landslide analysis.

Figure 3. Spatial distribution of precipitation factors. (a) PRCPTOT; (b) SDII; (c) RX1D; (d) RX5D; (e) R95P; (f) R99P; (g) R10MM; (h) R20MM; (i) R50MM.

Table 2. Definition of precipitation index.

Precipitation index	Abbreviation	Definition	Unit
Annual total wet-day precipitation	PRCPTOT	Annual total precipitation from days ≥ 1 mm	mm
Simple daily intensity index	SDII	The ratio of annual total precipitation to the number of wet days (≥1 mm)	mm/days
Max 1-day precipitation amount	RX1D	Monthly maximum 1-day precipitation	mm
Max 5-day precipitation amount	RX5D	Monthly maximum consecutive 5-day precipitation	mm
Very wet days	R95P	Annual total precipitation from days > 95th percentile	mm
Extremely wet days	R99P	Annual total precipitation from days > 99th percentile	mm
Number of heavy precipitation days	R10MM	Annual count when precipitation ≥ 10 mm	Days
Number of very heavy precipitation days	R20MM	Annual count when precipitation ≥ 20 mm	Days
Number of extreme heavy precipitation days	R50MM	Annual count when precipitation ≥ 50 mm	Days

3.1.2. Surface environment factors

In the study of landslide environmental impacts, the role of topography and land surface characteristics is pivotal. Elevation has consistently been identified as one of the most common influencing factors in landslide research (Reichenbach et al. 2018). Terrain parameters derived from the DEM, such as Slope and Aspect (SA), directly impact surface water flow and soil permeability. An increase in slope gradient signifies reduced soil infiltration and increased surface runoff, thereby contributing to landslide occurrence. SA influences precipitation and sunlight exposure, subsequently affecting soil moisture and landslide risk (Van den Eeckhaut et al. 2006). The Topographic Wetness Index (TWI) highlights the influence of terrain on the process of moisture accumulation, while Surface Roughness (SR) reflects the irregularity of the terrain (Hong et al. 2018). Topographic Relief (TR) indicates the elevation difference between the highest and lowest points within a specified area, establishing a direct correlation between landform morphology and potential landslide occurrence (Schmidt and Montgomery 1995; Qiu et al. 2018). Soiltype impacts water discharge based on its texture, permeability, and structure, while Lithology controls the water percolation process, both exerting direct effects on landslide occurrence (Trigila et al. 2015; Mojaddadi et al. 2017). LULC modify topsoil properties, affecting water storage and surface runoff, thereby indirectly influencing landslide susceptibility (Shu et al. 2019). Distance from rivers (DistRivers) influences landslide occurrence as water bodies can undermine slope stability through erosion and increased saturation, which weakens soil and rock strength near riverbanks (Gokceoglu et al. 2005). Distance from faults (DistFaults) is a critical factor in landslide susceptibility, as seismic activities associated with faults can trigger landslides by inducing ground shaking and fracturing the earth’s surface. Lastly, NDVI reflects the area’s plant biomass, indicating vegetation coverage’s role in regulating runoff processes and indirectly influencing landslide risk (Chu et al. 2009; Khosravi et al. 2016). All the aforementioned data were resampled to a spatial resolution of 1 km, and the spatial distribution of all factors is depicted in Figure 4.

Figure 4. Spatial distribution of the surface environment factors. (a) Elevation; (b) SA; (c) Slope; (d) LULC; (e) SR; (f) TR; (g) TWI; (h) Lithology; (i) NDVI; (j) DistRivers; (k) DistFaults; (l) Soiltype.

3.1.3. Human activity factors

Human activities play a pivotal role in influencing landslide occurrences, as they often significantly alter the natural landscape and its stability. In this study, we examined two critical human activity factors: POP and distance from roads (DistRoads). These factors were used to assess the intensity and extent of human activities. POP indicates the degree of population concentration within the area, while DistRoads demonstrates the development of the transportation network. Together, they contribute to landslide risk. The spatial distribution of these factors is meticulously illustrated in Figure 5.

Figure 5. Spatial distribution of human activity factors. (a) POP; (b) DistRoads.

3.2. Methods

3.2.1. Spearman correlation coefficient

In this study, the Spearman correlation coefficient (SCC) was used to assess the relationship between various landslide-contributing factors and landslide occurrences. Introduced by Charles Spearman in 1904 (de Winter et al. 2016), the SCC is a nonparametric measure of dependency, suitable for evaluating the correlation between rank-based variables. Its application allows for the precise quantification of the association between these influencing factors and landslide events. The SCC is defined as follows: (1) $$\rho =\frac{\mathrm{\text{Cov}}\left(R\left(A_i\right),R\left(B_i\right)\right)}{\sqrt{\mathrm{v}\mathrm{\text{ar}}\left(R\left(A_i\right)\right)\hspace{0.17em}\bullet \hspace{0.17em}\mathrm{\text{var}}\left(\mathit{R}\left({\mathit{B}}_{\mathit{i}}\right)\right)}}$$(1) where $A_i$ represents the 23 landslide-contributing factors, and $B_i$ denotes the total number of landslide occurrences within the corresponding grid cells.

3.2.2. Interaction detector

The geodetector, a spatial analysis method designed to detect spatial heterogeneity and its driving forces, was proposed by Wang in 2010 (Wang et al. 2010). It consists of four primary types of detectors: factor detector, interaction detector, risk detector, and ecological detector. In this study, the interaction detector is predominantly utilized to examine the interplay between two factors in relation to landslides, and the extent of their relationship is presented in Table 3. This degree of relationship is measured by the q-value, calculated using following formula: (2) $$q=1-\frac{\sum _{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}=1-\frac{\mathrm{\text{SSW}}}{\mathrm{\text{SST}}}$$(2) (3) $$\mathrm{\text{SSW}}=\sum _{h=1}^L{N}_h{\sigma }_h^2$$(3) (4) $$\mathrm{\text{SST}}=N{\sigma }^2$$(4)

Table 3. Types of interaction between two covariates.

Description	Interaction
q(X1∩X2) < Min(q(X1), q(X2))	Weakened, nonlinear
Min(q(X1), q(X2)) < q(X1∩X2) < Max(q(X1), q(X2))	Weakened, single factor nonlinear
q(X1∩X2) > Max(q(X1), q(X2))	Enhanced, bi
q(X1∩X2) = q(X1) + q(X2)	Independent
q(X1∩X2) > q(X1) + q(X2)	Enhanced, nonlinear

In formula (2)(2) $$q=1-\frac{\sum _{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}=1-\frac{\mathrm{\text{SSW}}}{\mathrm{\text{SST}}}$$(2) : $h$ = 1,…, L is the stratum of variable Y or factor X; $N_h$ and $N$ are the number of units in layer $h$ and the whole area respectively; ${\sigma }_h^2$ and ${\sigma }^2$ are the variances of the Y value of layer h and the whole area respectively. $\mathrm{\text{SSW}}$ and $\mathrm{\text{SST}}$ are within the sum of squares and the total sum of squares.

3.2.3. Model

3.2.3.1. Random Forest

The RF algorithm, initially proposed by Breiman (2001), is an ensemble machine-learning approach that operates based on multiple decision trees for both data training and prediction. The core components of a RF consist of individual decision trees and the bagging technique. As a hierarchically structured classifier, a decision tree classifies data through its levels. Bagging contributes by randomly drawing subsets from the training data to construct each tree, while the remaining data (out-of-bag data) are utilized for cross-validation to assess the performance of the RF classifier. RF is known for its superior performance compared to other machine learning algorithms in aspects such as parameter tuning, training speed, noise reduction, and overfitting prevention. Moreover, it efficiently manages and processes large-scale and high-dimensional data, making the RF an ideal choice for handling remote sensing information.

3.2.3.2. Light Gradient Boosting Machine

LGBM is an efficient algorithm based on Gradient Boosting Decision Tree frameworks, employing gradient-based one-sided sampling, histogram algorithms, exclusive feature bundling, and leaf-wise tree growth strategies to enhance training speed, computational efficiency, and reduce memory consumption (Song et al. 2024). The histogram algorithm discretizes continuous feature values and constructs histograms to identify optimal split points, while the leaf-wise growth strategy, which allows splitting only for the leaf with the maximum information gain, reduces the number of tree nodes. This strategy not only boosts splitting efficiency but also prevents overfitting, offering LGBM a distinct advantage in processing large-scale data.

3.2.3.3. Long Short-Term Memory

LSTM networks represent a specialized subclass of Recurrent Neural Networks (RNNs) that are designed to learn long-term dependencies within data sequences. Distinct from conventional RNNs, LSTMs are equipped with integrated mechanisms that facilitate the selective retention and discarding of information during the learning process. This architecture enhances the model’s ability to manage data over extended time intervals, thereby mitigating issues related to vanishing or exploding gradients commonly encountered in standard RNNs (Chen and Fan 2023b).

3.2.4. Optuna-based hyperparameter optimization framework

To enhance the performance of the RF model in landslide identification applications, the Optuna optimization framework for hyperparameter tuning was adopted in this study. Optuna (Akiba et al. 2019) is an advanced hyperparameter optimization framework that selects the optimal combination of hyperparameters by optimizing a predefined objective function. The framework’s core components include the objective function, individual trials, and the optimization process. In contrast to traditional grid search methods, Optuna employs sampling-based approaches and pruning strategies to select hyperparameters, effectively shortening the model fitting time while improving performance. Integration of the Optuna framework with the RF algorithm enabled more precise and efficient identification of landslides in this study. The process for optimizing hyperparameters with Optuna involves four main steps: (1) define the objective function, which is to maximize the F1 score, and specify the range of hyperparameters for the classification model; (2) in each trial, train the classification model using the given hyperparameters, predict on the validation data, and calculate the F1 score; (3) perform multiple trials to find the hyperparameters that yield the highest F1 score, thus optimizing the classifier’s performance; (4) output the best set of hyperparameters and the corresponding F1 score. The hyperparameters search spaces for all classifier models and the best parameters will be displayed in Table 4.

Table 4. Hyperparameter search spaces of three models.

Classification model	Hyperparameters search spaces	Best parameters
RF	n_estimators: 10-1000; max_depth: 2-32; min_samples_split: 2-20; max_feature: [0.5, ’sqrt’, ’log2’].	n_estimators: 483; max_depth: 17; min_samples_split: 3; max_feature: 0.5.
LGBM	n_estimators: 200-1000; max_depth: 3-8; num_leaves, 8-256; reg_alpha: 1e-3-10.0; reg_lambda 1e-3-10.0; colsample_bytree: 0.4-1.0; learning_rate: 0.01-0.3; subsample: 0.4-1.0; min_child_samples: 5-100.	n_estimators: 535; max_depth: 5; num_leaves, 109; reg_alpha: 0.3764; reg_lambda 0.0018; colsample_bytree: 0.533; learning_rate: 0.024; subsample: 0.562; min_child_samples: 12.
LSTM	batch_size: [32, 64, 128, 256, 512]; hidden_size: [32, 64, 128, 256, 512]; learning_rate: 0.001-0.01.	batch_size: 32; hidden_size: 256; learning_rate: 0.0079.

3.2.5. Model evaluation

Landslide susceptibility models are commonly evaluated using confusion matrices and ROC curves. The ROC curve, which plots false positive rates (1-specificity) on the x-axis and true positive rates (sensitivity) on the y-axis, provides a visualization of the model’s performance. The AUC serves as a critical indicator of the predictive capabilities of the model (Fawcett 2006). Confusion matrices offer an intuitive understanding of model predictions, including instances of misclassification. Key metrics such as ACC, PPV, TPR, and F1 are used to assess the overall model performance. These metrics are defined as follows: (5) $$\mathrm{\text{Accuracy}}\mathrm{ }\left(\mathrm{\text{ACC}}\right)=\frac{\mathrm{\text{TP}}+\mathrm{\text{TN}}}{\mathrm{\text{TP}}+\mathrm{\text{FP}}+\mathrm{\text{TN}}+\mathrm{\text{FN}}}$$(5) (6) $$\mathrm{\text{Precision}} \left(\mathrm{\text{PPV}}\right)=\frac{\mathrm{\text{TP}}}{\mathrm{\text{TP}}+\mathrm{\text{FP}}}$$(6) (7) $$\mathrm{\text{Recall}}\mathrm{ }\left(\mathrm{\text{TPR}}\right)=\frac{\mathrm{\text{TP}}}{\mathrm{\text{TP}}+\mathrm{\text{FN}}}$$(7) (8) $$\mathrm{F}1‐\mathrm{\text{score}}\left(\mathrm{F}1\right)=\frac{2\text{Precision*Recall}}{\mathrm{\text{Precision}}+\mathrm{\text{Recall}}}$$(8)

where $\mathrm{\text{TP}}$ means that the actual value and predicted value are both positive; $\mathrm{\text{FP}}$ means that the actual value is negative and the predicted value is positive; $\mathrm{\text{FN}}$ means that the actual value is positive and the predicted value is negative; $\mathrm{\text{TN}}$ means that both the actual value and the predicted value are negative.

3.2.6. SHAP interpreter

SHAP is a game-theoretic approach used to interpret machine learning models, as described by Lundberg and Lee (2017). Its primary advantage lies in the ability to quantify the contribution of each feature to every individual prediction, including the direction of these contributions, whether positive or negative. This holistic approach quantifies feature importance and also accounts for the complex interactions that underlie model predictions. The essence of the SHAP method involves calculating the marginal contribution of each feature to the model’s output, known as Shapley values. These values facilitate explanations of the model’s ‘black box’, both globally and locally. The mathematical representation for computing these marginal contributions is as follows: (9) $${\varphi }_i=\sum _{S\subseteq N\left\{i\right\}}\frac{\left|s\right|!\left(n-\left|s\right|-1\right)!}{n!}\left[f\left(s\cup \left\{i\right\}\right)-f\left(s\right)\right]$$(9) where ${\varphi }_i$ is the contribution of feature $\mathrm{ }i,$ $N$ is the set containing all features, $n$ is the number of features in $N,$ $s$ is the subset of $N$ that contains feature $i,$ $f\left(s\cup \left\{i\right\}\right)$ and $f\left(s\right)$ denote the model results with or without the feature $i,$ respectively (Lundberg et al. 2019).

SHAP quantifies the contribution of each feature to model predictions and interprets the predictions as the sum of Shapley values assigned to each input feature: (10) $$g={\varphi }_0+\sum _{i=1}^n{\varphi }_i$$(10) where $g$ is the value of the model, ${\varphi }_0$ is the constant that explains the model (i.e. the predicted mean of all training samples), and ${\varphi }_i$ is the imputed value (Shapley value) of each feature.

SHAP provides a quantitative measure of each feature’s contribution to the model’s predictions, revealing not only the significance of each feature but also how these features interact to influence the prediction outcomes. This insight is invaluable for decision-makers, offering a clear understanding of which features are pivotal and how they shape the predictive results. Such knowledge is crucial for developing targeted, data-informed strategies, enhancing the decision-making process.

3.2.7. Computational efficiency

The duration of model training is significantly influenced by their complexity and the extent of their hyperparameter search space. In this investigation, the RF model necessitates the construction of numerous decision trees to attain optimal performance, resulting in a training time of approximately 40 min. Conversely, the LGBM demonstrates enhanced efficiency in processing large datasets, requiring merely 2 min for training. Both models underwent training on a CPU. Given the substantial computational demands associated with deep learning models, the LSTM’s training was facilitated using a GPU, culminating in a duration of roughly 6 min. All models in this study were developed using Python 3.9. The RF was implemented with the scikit-learn 1.2.1 library, LGBM was facilitated using the lightgbm 4.3.0 library, and LSTM was constructed with the torch 1.7.1 library, utilizing the PyCharm Community Edition 2022.3.2 as the integrated development environment. Computations were performed on a system equipped with an Intel(R) Xeon(R) W-2295 CPU @ 3.00 GHz, 64.0 GB of RAM, and an NVIDIA Quadro RTX 4000 GPU.

4. Results

4.1. Correlation analysis between influencing factors and landslides

Rational selection of influencing factors is of paramount importance in constructing geological disaster models. Regional landslides result from the combined effects of various influencing factors, which often exhibit a certain degree of correlation. This study focuses on three major aspects: precipitation, human activities, and surface environment. A total of 23 key influencing factors were selected, and their correlations were analyzed using the SCC. The results, depicted in Figure 6, revealed negative correlations between landslides and factors such as Elevation, DistRoads, DistRivers, LULC, Slope, SR, and TR. Other influencing factors showed positive correlations. It is noteworthy that POP demonstrated the most significant correlation with landslide events, with a SCC of 0.182, emphasizing the substantial impact of human activities on landslide occurrence. Furthermore, Elevation exhibited a SCC of −0.151, suggesting that landslides are more likely to occur at lower elevations. DistRoads, with a SCC of −0.133, was also a significant factor that cannot be overlooked. Among these 23 influencing factors, all except DistFaults and SA reached a significance level of 0.05 for their SCC, indicating their important connections to the occurrence of landslides.

Figure 6. Spearman correlation coefficient between landslides and influence factors.

4.2. Bivariate interaction effects

The occurrence of landslides results from the combined influence of multiple factors. To investigate the interrelationships among these factors, the interaction detector from the geodetector toolkit was utilized (Song et al. 2020). Analysis of Figure 7 revealed two primary interaction modes between pairs of factors: bivariate enhancement and nonlinear enhancement. This discovery suggests that the interpretive power of landslides is amplified when considering the combined effects of multiple factors. In other words, the interaction between two factors has a more substantial impact on landslides than their independent effects. Larger points in Figure 7 represent higher degrees of landslide impact by the two interacting factors. Among all the bivariate interactions, the most significant influences on landslide occurrence were identified as the interaction between Elevation and LULC, LULC with Soiltype, and Elevation and RX5D. In summary, the cumulative effect of multiple factors increases the likelihood of landslide occurrence. Moreover, Elevation, POP, and Soiltype exhibited strong explanatory power for landslides when combined with various other factors. The combinations involving DistFaults and SA with other factors resulted in a nonlinear enhancement of the explanatory power for landslides.

Figure 7. Interaction effects of influence factors of landslide based on geodetector (red indicates enhanced, bi; blue indicates enhanced, nonlinear).

4.3. Landslide susceptibility assessment

4.3.1. Model construction and evaluation

In this study, a significant dataset imbalance was observed, with the number of non-landslide instances far exceeding the number of landslide instances. To address this issue, all landslide samples were initially extracted, followed by the extraction of an equal number of non-landslide samples from the study area. This process adhered to two primary principles: ensuring a uniform distribution of the non-landslide samples throughout the study area, and halting sampling when a 1:1 ratio of landslide to non-landslide samples was achieved to attain sample balance. As a result, the dataset comprised a total of 7378 records, with an equal representation of landslide and non-landslide data, which were then divided into a training set and a test set in an 8:2 ratio.

The model training and refinement process involved the integration of RF, LGBM, and LSTM models with the Optuna framework for hyperparameter optimization. Optimal hyperparameters were searched within the predefined hyperparameter space. The F1 score was utilized as the benchmark for performance valuation, while the TPE was adopted as the searching technique within the realm of hyperparameter optimization. As shown in Table 5, RF model excelled during the training phase, registering laudable outcomes on various evaluative metrics, including ACC, PPV, TPR, and the F1, which were recorded at 0.9615, 0.9278, 1, and 0.9625, respectively. These metrics substantiate the model’s exceptional fidelity to the data. On the test dataset, the model sustained commendable performance, with an ACC of 0.7792 and an F1 of 0.8055. While the LGBM and LSTM were marginally outperformed by the RF model in the training dataset, the LGBM showcased superior PPV on the test set, and the LSTM model was notably distinguished by a TPR of 0.9124 on the same dataset. Figure 8 delineates the AUC for the trio of models across both training and test datasets, with the RF model achieving AUC values of 0.9997 and 0.8387, respectively, thereby outperforming the LGBM and LSTM models. In conclusion, the RF model stands out for its paramount overall performance, underlining its potent predictive capacity for real-world applications, characterized by a high degree of predictive accuracy and reliability.

Figure 8. (a) Receiver operating characteristic curves for the training set; (b) Receiver operating characteristic curves for the test set.

Table 5. Results of model accuracy evaluation.

Data set	Model	ACC	PPV	TPR	F1
Training set	RF	0.9615	0.9278	1	0.9625
	LGBM	0.8684	0.8155	0.9485	0.8770
	LSTM	0.7304	0.6965	0.8067	0.7475
Test set	RF	0.7792	0.7448	0.8769	0.8055
	LGBM	0.7771	0.7488	0.8613	0.8012
	LSTM	0.7459	0.6945	0.9124	0.7887

4.3.2. Distribution pattern of landslide susceptibility

In this study, landslide susceptibility was predicted using the hyperparameter-optimized RF model implemented with the Optuna framework. The susceptibility levels were classified into five categories: very low (0–0.2), low (0.2–0.4), moderate (0.4–0.6), high (0.6–0.8), and very high (0.8–1) (He et al. 2021). Analysis of the susceptibility map (Figure 9) reveals that landslide-prone areas are predominantly concentrated in the southern part of the study area, showing a patchy distribution pattern. It is preliminarily inferred that this distribution pattern may be significantly influenced by human activities. Additionally, in the northern part of the region, particularly around the Nujiang and Lancang Rivers, areas prone to landslides display noticeable signs of development along river courses, presenting a linear and densely populated distribution that could be attributed to topography.

Figure 9. Landslide susceptibility map.

A detailed statistical analysis of the landslide susceptibility zones (Table 6) revealed that very low risk areas constitute 35.23% of the total study area, yet historically, only 0.24% of landslide events occurred in these zones. In contrast, despite the high and very high-risk zones collectively occupying just 28.17% of the study area, they accounted for 93.79% of the landslide occurrences, with densities of 0.0913/km² and 0.5329/km² respectively in these areas. Additionally, the results indicate a correlation between the level of landslide susceptibility and the density of landslide occurrences, consistent with typical characteristics of geohazard-prone areas. Importantly, the model’s predictions demonstrated a high degree of alignment with observed data, thus affirming the model’s accuracy and practical utility.

Table 6. Landslide susceptibility zoning statistics.

Landslide susceptibility	Number of grids	Percentage of area/%	Number of landslides	Percentage of landslides/%	Landslide density/(Pcs·km^2)
Very low	23344	35.23	9	0.24	0.0004
Low	9708	15.48	43	1.17	0.0044
Moderate	11987	19.12	177	4.80	0.0148
High	13482	21.50	1231	33.37	0.0913
Very high	4183	6.67	2229	60.42	0.5329

4.4. Interpretability analysis of landslide influence factors

In this study, we employed the SHAP interpreter to provide a detailed interpretation of the RF model for investigating the influence mechanisms of different factors on the spatial distribution of landslides. Figure 10(a) presents the ranking of factor importance obtained from the mean absolute SHAP values for each sample, highlighting Elevation, POP, DistRoads, NDVI, and LULC as the most significant factors contributing to landslides. In Figure 10(b), each feature is represented by a separate row, with the horizontal axis reflecting SHAP values. Each point corresponds to a sample, and the color scale ranging from redder to bluer represents higher to lower feature values, respectively. The analysis reveals a significant negative correlation between Elevation and landslide occurrence within the study area, indicating that lower elevations have a more pronounced positive impact on the occurrence of landslides. Higher NDVI values, reflecting greater vegetation cover, are associated with a reduced probability of landslide events. As anticipated, factors strongly associated with human activities, namely POP and DistRoads, have been found to significantly amplify the risk of landslide incidents. Additionally, LULC types comprising water bodies, urban areas, rural settlements, and other constructed areas have been identified as playing a crucial role in landslide risk. These observations further validate the prominence of human activities as essential factors influencing landslide occurrences within the study area.

Figure 10. SHAP-based feature importance. (a) Factor importance plot; (b) summary plot of SHAP.

5. Discussion

5.1. Advantages of the proposed model

This study presents an interpretable model that integrates RF, Optuna framework, and SHAP. This unified model demonstrates significant advantages in optimal landslide susceptibility assessment and explanatory capabilities. First, the RF model, as an ensemble classifier using decision trees as base models, has been widely utilized in various domains, including landslide susceptibility analysis, due to its noise reduction capability, ability to address overfitting, and apt handling of large, high-dimensional datasets (Kim et al. 2018; Liu et al. 2018; Wu et al. 2023). Furthermore, Optuna is acknowledged as a leading hyperparameter optimization framework tailored for machine and deep learning applications, offering seamless integration with the RF model (Imani and Arabnia 2023). Lastly, SHAP, a machine learning explanation method rooted in game theory and local interpretability, is employed to elucidate model outputs, with a particular focus on the well-established application of its tree-based models to the RF model (Lundberg et al. 2020; Sun et al. 2023). The application of the proposed model in landslide susceptibility assessment, as evidenced in Table 5 and Figure 8, demonstrates remarkable performance on both training (ACC = 0.9615, PPV = 0.9278, TPR = 1, F1 = 0.9625, AUC = 0.9997) and test datasets (ACC = 0.7792, PPV = 0.7448, TPR = 0.8769, F1 = 0.8055, AUC = 0.8387).

The effects of hyperparameter optimization on RF were compared using different methods: manual tuning, random search, grid search, Bayesian optimization, and the Optuna framework. Table 7 presents the performance metrics for RF model under different hyperparameter optimization methods. Regardless of dataset, the utilization of the Optuna framework for hyperparameter optimization yielded improved model performance for the RF model. Furthermore, previous research has demonstrated the effectiveness of the Optuna framework in mitigating model overfitting (Zhang et al. 2023).

Table 7. Performance metrics for RF model under different hyperparameter optimization methods.

Data set	Methods	ACC	PPV	TPR	F1	AUC
Training set	Manual tuning	0.8648	0.7936	0.9821	0.8778	0.9423
	Grid search	0.8651	0.7926	0.9852	0.8784	0.9438
	Random search	0.8651	0.7926	0.9852	0.8784	0.9438
	Bayesian optimization	0.8853	0.8143	0.9952	0.8957	0.9664
	Optuna	0.9615	0.9278	1	0.9625	0.9997
Test set	Manual tuning	0.7650	0.7274	0.8782	0.7957	0.8308
	Grid search	0.7663	0.7295	0.8769	0.7964	0.8292
	Random search	0.7663	0.7295	0.8769	0.7964	0.8292
	Bayesian optimization	0.7697	0.7294	0.8873	0.8007	0.8313
	Optuna	0.7792	0.7448	0.8769	0.8055	0.8387

The importance of features in RF model is determined by averaging the contributions of features in splits (Figure 11), while SHAP values calculate the impact of features on the predictive outcome. When Figures 10(a) and 11 are compared, the rankings of most features remain consistent across both methods, especially for factors that significantly influence the model, such as Elevation, POP and NDVI. However, discrepancies exist in the ranking of less important features between the two methods, such as LULC, Soiltype, SA, and TWI. These inconsistencies may arise from the fact that the average contributions of features evaluated by RF model do not entirely correspond to the marginal contributions measured by SHAP values. The SHAP offers a more detailed perspective for elucidating model decisions, beyond solely ranking the importance of features; it also reveals the positive and negative effects of features and their distribution, providing a more comprehensive and stable explanation for the model’s predictive outputs.

Figure 11. Importance of features built into the RF model.

5.2. Impacts of input drought factors on model results

To assess the model’s robustness to input data variability and potential biases, the leave-one-out method was utilized in this study for a systematic sensitivity analysis involving 23 landslide factors (Xiao et al. 2024), as shown in Table 8. Each factor was sequentially excluded from the model, which then utilized the remaining 22 factors to forecast landslides, allowing for an evaluation of each factor’s impact on the prediction outcomes. It was observed that model performance consistently declined upon the exclusion of any factor. Notably, sensitivity was highest with respect to the POP factor; absence of POP resulted in the lowest values across all performance indicators. Nonetheless, while the omission of any factor invariably reduced performance, the magnitude of decline was generally modest. Performance metrics such as ACC, PPV, and F1 remained above 0.72, with TPR and AUC exceeding 0.83. These results highlight the critical role of all input landslide factors in sustaining model performance and demonstrate that the model effectively maintains data sensitivity without compromising predictive accuracy.

Table 8. Model prediction accuracy of the RF model using 22 or 23 drought factors.

	ACC	PPV	TPR	F1	AUC
All	0.7792	0.7448	0.8769	0.8055	0.8387
Elevation	0.7603	0.7377	0.8381	0.7847	0.8335
DistFaults	0.7738	0.7425	0.8666	0.7998	0.8373
DistRoads	0.7549	0.7250	0.8536	0.7841	0.8307
DistRivers	0.7725	0.7425	0.8627	0.7981	0.8380
Lithology	0.7752	0.7436	0.8679	0.8010	0.8386
LULC	0.7738	0.7414	0.8692	0.8002	0.8345
NDVI	0.7664	0.7346	0.8640	0.7940	0.8305
POP	0.7535	0.7294	0.8381	0.7800	0.8308
PRCPTOT	0.7738	0.7420	0.8679	0.8000	0.8396
R10MM	0.7677	0.7357	0.8653	0.7952	0.8385
R20MM	0.7677	0.7362	0.8640	0.7950	0.8379
R50MM	0.7738	0.7436	0.8640	0.7993	0.8385
R95P	0.7758	0.7450	0.8666	0.8012	0.8394
R99P	0.7758	0.7455	0.8653	0.8010	0.8390
RX1D	0.7670	0.7364	0.8614	0.7940	0.8371
RX5D	0.7725	0.7409	0.8666	0.7988	0.8368
SA	0.7765	0.7458	0.8666	0.8017	0.8391
SDII	0.7711	0.7398	0.8653	0.7976	0.8367
Slope	0.7752	0.7453	0.8640	0.8002	0.8390
Soiltype	0.7704	0.7384	0.8666	0.7974	0.8378
SR	0.7725	0.7409	0.8666	0.7988	0.8374
TR	0.7758	0.7418	0.8744	0.8026	0.8385
TWI	0.7731	0.7412	0.8679	0.7995	0.8372

5.3. Uncertainty in this study

The number of factors influencing landslide susceptibility is extensive, and the primary factors vary by region (Carrara et al. 2008; Xie et al. 2017; Wang et al. 2020; Berhane et al. 2021; Sun et al. 2021). Although this study has considered 23 landslide factors, whether these factors are sufficient to meet the needs for landslide susceptibility prediction in all regions remains a question worth considering.

Additionally, considering the diversity and varying spatial scales of the data used in this study, resampling was employed to standardize the spatial scales for computational convenience. However, this resampling process may inadvertently omit critical information corresponding to actual landslide factors, potentially diminishing the model’s predictive accuracy.

Furthermore, previous research indicates that the selection of sampling areas and the ratios of non-landslide to landslide samples significantly impact susceptibility assessment results (Hong et al. 2019). In this study, non-landslide and landslide samples were randomly sampled at a 1:1 ratio, which might adversely affect the evaluation outcomes. There is a potential risk during the random sampling process of misclassifying areas prone to landslides as non-landslide samples, possibly leading to a misinterpretation of their true susceptibility. Additionally, the 1:1 ratio may not reflect the actual proportion in the real world, potentially affecting the model’s accuracy and sensitivity.

In light of these considerations, future work will focus on developing a method that both unifies the scale of landslide factors and minimizes information loss, as well as identifying the optimal sampling regions and proportions for non-landslide samples to enhance the model’s predictive accuracy.

5.4. Model generalizability

The model proposed in this study has been applied to the northwestern Yunnan region, demonstrating potential for generalizability across areas with diverse geological and environmental characteristics. Despite robust performance in the specified region (F1 = 0.8055 and AUC = 0.8387), the model’s generalizability may require adjustments when deployed in varying environments. High-quality meteorological, terrestrial, and human activity data are crucial for model construction. In regions where such data are lacking, the model may not achieve the expected results. This study identified 23 landslide factors from precipitation, surface environments, and human activities, insufficient to meet the landslide susceptibility forecasting needs of all regions. However, for areas heavily influenced by specific conditions, it may be necessary to incorporate additional landslide factors to enhance the predictive adaptability and accuracy.

6. Conclusions

This study presents an interpretable model that integrates RF model, Optuna hyperparameter optimization, and SHAP to achieve optimal evaluation of landslide susceptibility and provide explanations in Northwest Yunnan Province, China. An inventory of 4447 landslides and 23 related factors was used for the landslide susceptibility assessment. The relevance of these landslide factors was ensured through the utilization of the SCC and geodetector. The results indicated that 21 out of the 23 landslide factors exhibited a significant correlation with landslides, and the interaction of any two different factors enhanced the explanatory power for landslides. By leveraging the hyperparameter-optimized RF model using the Optuna framework, in conjunction with landslide factors and historical landslides, spatial distribution maps of landslide susceptibility were generated. The model’s performance analysis revealed outstanding results on the test set, with an ACC of 0.7792, F1 of 0.8055, and an AUC of 0.8387. Furthermore, the interpretation of the model was facilitated through the utilization of SHAP. The results emphasized the paramount importance of factors such as Elevation, POP, DistRoads, NDVI, and LULC in determining landslide susceptibility. Areas at lower elevations, especially those that are aquatic, heavily impacted by human activities, or characterized by sparse vegetation coverage, are more susceptible to landslides. The proposed new interpretable model offers a fresh perspective for evaluating landslide susceptibility at a regional scale. Moving forward, further exploration of advanced machine learning models and comprehensive landslide factors holds promise for advancing this field.

Data availability statement

The data presented in this study are available on request from the corresponding author.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Key Research and Development Program of China [Grant No. 2022YFF1302405].