A New Approach to Spatial Landslide Susceptibility Prediction in Karst Mining Areas Based on Explainable Artificial Intelligence

Abstract

Landslides are a common and costly geological hazard, with regular occurrences leading to significant damage and losses. To effectively manage land use and reduce the risk of landslides, it is crucial to conduct susceptibility assessments. To date, many machine-learning methods have been applied to the landslide susceptibility map (LSM). However, as a risk prediction, landslide susceptibility without good interpretability would be a risky approach to apply these methods to real life. This study aimed to assess the LSM in the region of Nayong in Guizhou, China, and conduct a comprehensive assessment and evaluation of landslide susceptibility maps utilizing an explainable artificial intelligence. This study incorporates remote sensing data, field surveys, geographic information system techniques, and interpretable machine-learning techniques to analyze the sensitivity to landslides and to contrast it with other conventional models. As an interpretable machine-learning method, generalized additive models with structured interactions (GAMI-net) could be used to understand how LSM models make decisions. The results showed that the GAMI-net model was valid and had an area under curve (AUC) value of 0.91 on the receiver operating characteristic (ROC) curve, which is better than the values of 0.85 and 0.81 for the random forest and SVM models, respectively. The coal mining, rock desertification, and rainfall greater than 1300 mm were more susceptible to landslides in the study area. Additionally, the pairwise interaction factors, such as rainfall and mining, lithology and rainfall, and rainfall and elevation, also increased the landslide susceptibility. The results showed that interpretable models could accurately predict landslide susceptibility and reveal the causes of landslide occurrence. The GAMI-net-based model exhibited good predictive capability and significantly increased model interpretability to inform landslide management and decision making, which suggests its great potential for application in LSM.

1. Introduction

In Karst regions, landslides are a significant natural hazard distinguished by their widespread occurrence, high frequency, and tremendous damage [1,2]. The frequency, size, and destructiveness of landslides have recently increased due to enhanced human activity and frequent catastrophic weather occurrences [3,4]. Landslides provide a hazard to 74 million people in China alone, with yearly direct economic damages typically running into the hundreds of millions [5]. While landslides may not be entirely avoidable, managers may create effective mitigation measures to lessen mortality and damage by researching landslide condition factors, modeling landslide risk, and mapping and interpreting hazard-prone locations. Landslide susceptibility assessments are valuable for disaster prevention, land resource planning, and risk assessment in landslide-prone areas.

Landslide susceptibility assessment is a method of quantitatively and spatially estimating the likelihood of a landslide occurring in a specific area [6]. It considers various factors contributing to landslides, such as topographic features including slope, aspect, and material properties such as soil type and lithology [7]. In addition, environmental factors such as climate, hydrology, ecology, human activities, and seismic activity could also play a role in the occurrence of landslides [8]. Recent advances in remote sensing technology and spatial datasets have provided valuable data for geological hazard monitoring and landslide susceptibility assessment [9]. The use of remote sensing allows for the collection of a wide range of data, including optical and microwave observations, which can provide information on normalized difference vegetation index (NDVI), soil moisture content, and surface deformation [10,11]. Zhang and Shen et al. discuss the effect of vegetation distribution on LSM implementation and the use of Interferometric Synthetic Aperture Radar (InSAR) surface deformation information to achieve dynamic refinement of LSM [12,13]. However, there are differences in geography, topography, and climatic conditions in different areas, and landslide patterns, mechanisms, and drivers vary in these areas. Therefore, landslide influencing factors need to be selected accordingly to the characteristics of the study area. This study considered the impact of frequent mining activities on landslides in the Nayong area.

Researchers have used landslide influencing factors for landslide susceptibility assessment through physical models, statistical analysis methods, or machine learning techniques [14,15,16,17,18,19,20,21]. Researchers used finite element models and limit equilibrium methods in the numerical simulation methods to assess slope stability [14]. This method simplifies assumptions about some environmental conditions (e.g., soil moisture content and rock characteristics) and landslide structures, and therefore requires costly field exploration and is not conducive to large-scale regional monitoring [15]. Statistical methods are reproducible and objective, and typical methods include logistic regression (LR), frequency ratio methods, and index of entropy [16,17,18]. Machine learning is the current mainstream method for LSM, which has significantly improved the accuracy and precision of landslide susceptibility assessment models. Common machine learning methods include support vector machines (SVM), random forest (RF), and deep neural networks [19,20,21].

Previous studies have shown that machine learning (ML) approaches and deep learning (DL) methods can assess landslide susceptibility in various area [22]. However, researchers have tended to neglect the explanatory aspects of models in seeking higher accuracy. Because of the ‘black box model’, these methods have difficulty explaining how landslide factors contribute to and dominate the occurrence of all or individual landslides in the study area. The lack of interpretation is a significant drawback of the above machine models and a fundamental flaw in their application to high-risk events (e.g., landslides, earthquakes, and lives) [23]. In these applications, decisions can affect lives, property, and generate huge costs. Based on the above analysis, the study provided an explainable artificial intelligence (XAI) method, called generalized additive models with structured interactions (GAMI-net), for generating LSM by integrating landslide influence factors and explained the effects of variables (features) on predictions and their responses in predictors. The GAMI-net model is an interpretable model proposed by Professor Zhang Aijun in 2021, consisting of a generalized additive model and several artificial neural networks (ANN) and sub-networks [24]. The interpretable model results help decision makers better understand the factors contributing to landslide risk and identify areas that may require interventions such as slope stabilization or land management.

The study extracted the corresponding landslide impact factors from satellite remote sensing data and geological information. Specifically, the DEM data derived terrain elevation, slope, aspect, and curvature. The NDVI is a measure of vegetation cover that indicates dense vegetation, soil moisture content, and soil type and directly impacts soil quality movement. Lithology is an essential factor indicating the strength of the soil. Heavy rainfall may lead to soil saturation, resulting in rainfall-induced landslides. In addition, mining, roads, and rivers are also contributing factors to landslides. The GAMI-net model constructed LSM by these factors and provided a global interpretation of landslide susceptibility. Next, the performance of traditional models (LR, SVM, and RF) and GAMI-net models are evaluated and analyzed. Finally, the Shapley adaptive interpretation (SHAP) method was used to analyze the key factors leading to landslides in a typical landslide, taking the town of Zongling in Nayong County as an example.

2. Materials and Methods

2.1. Study Area and Landslide Inventory

The study area is located in Nayong County, southeast of Bijie City, Guizhou Province, at longitude 104°55′–105°38′ East and latitude 26°30′–27°05′ North (Figure 1). The Karst is widely developed in the region, the geological environment is complex, karst action is obvious, and long-term mining activities have caused frequent geological disasters in the area with strong concealment [25]. On 28 August 2017, a large-scale landslide occurred in the Qiaobian group of the Pusha community in Zhangjiawan Town, which caused more than 30 deaths [26]. Since 2010, two hundred ninety-three geological disaster events have been recorded, impacting 1716 people in Nayong [27]. As a result, landslide hazards seriously threaten the safety of people’s lives and property in the study area.

The Third Institute of Surveying and Mapping of Guizhou Province provided the landslides inventory, including a total of 293 historical landslides (2010–2020) in Nayong County. In addition, the same number of non-landslides were used as a reference group. According to the geological hazard survey results, there were 194 small landslides, 88 medium landslides, eight large landslides, and three massive landslides.

2.2. Data Preparation

Landslides are natural hazards influenced by various factors, including geometric conditions such as elevation, slope, and slope orientation, as well as hydrological factors such as the topographic moisture index (TWI) and distance from rivers [28,29]. In addition, geological and environmental factors, such as lithology, land use, and human activity conditions can also impact the likelihood of landslides [30,31]. These factors may vary depending on the region in which the hazard occurs. Accurate detection of landslides requires the consideration of multiple factors, including those related to both natural phenomena and human activities [32,33].

Based on the above considerations and relevant literature, 14 influencing factors were selected for the study, including topography and geology, hydrology, land cover, and human activities [8,9,11,28,29,30,31,32,33]. The spatial distribution of these influencing factors is shown in Figure 2. Table 1 illustrated the factors selected, resolution and data source. The main landslide factors are the following categories:

Topography: The steepness and orientation of a slope can also play a role in its susceptibility to landslides. Steep slopes are generally more prone to landslides than gentle ones, and slopes facing specific directions may be more susceptible depending on the direction of predominant weather patterns and the location of water sources.

Geology: The type of rock and soil present on a slope can affect its stability. For example, slopes composed of weak, poorly consolidated soils or those with a high clay content may be more prone to landslides.

Hydrology: Water’s presence can significantly affect a slope’s stability. Water can seep into the ground and weaken the soil, making it more prone to landslides. Additionally, heavy rain or rapid snowmelt can cause landslides by adding weight to the slope or increasing the soil’s water content.

Land use: Human activities can also contribute to landslide susceptibility. For example, the construction of buildings or roads on steep slopes can increase the weight on the slope and make it more prone to landslides. Additionally, vegetation removal can destabilize a slope by reducing its ability to absorb water.

Climate: Changes in temperature and precipitation patterns can also influence landslide susceptibility. For example, prolonged drought can cause soil to become dry and brittle, making it more prone to landslides, while heavy rain or rapid snowmelt can increase the water content of the soil and lead to landslides.

To gather data on these factors, the study employed a variety of sources, including DEM, ground surveys, Landsat-8, and Global Precipitation Measurement (GPM). The data processing implemented in this study consisted of the following steps.

(1) Based on Landsat-8 OIL images obtained from the Google Earth Engine, the NDVI of the study area was produced;

(2) Obtaining SRTM DEM data of the study area from USGS to calculate elevation, slope, aspect, plan curvature, and profile curvature (https://urs.earthdata.nasa.gov/ accessed on 12 August 2022);

(3) The GPM data was obtained from the USGS and converted to spatial raster data by Kriging interpolation and meteorological station monitoring data;

(4) Soil type, geological data, and land use and land cover (LULC) were all rasterized, and only resampling was required.

Table 1. List of landslide factors and data scale and source.

Class of Factors	Factors Selected	Scale	Data Source
Topography	Elevation	30 m	Calculated from DEM (NASA)
	Slope angle	30 m
	Aspect	30 m
	Topographic wetness index	30 m
	Profile curvature	30 m
	Plan curvature	30 m
Geolithology	Distance to an active fault	30 m	The Third Survey and Mapping Institute of Guizhou Province
	Lithology	30 m
Land use and land cover	Land use and land cover	30 m	Ministry of Natural Resources Data Center
Morphology	Distance to rivers and road	30 m	Openstreetmap
Climate	Rainfall	Vector/0.1°	Meteorological station monitoring data and GPM
Vegetation	NDVI	30 m	Landsat
Human	Coal mining	Vector	The Third Survey and Mapping Institute of Guizhou Province

Table 1. List of landslide factors and data scale and source.

Class of Factors	Factors Selected	Scale	Data Source
Topography	Elevation	30 m	Calculated from DEM (NASA)
	Slope angle	30 m
	Aspect	30 m
	Topographic wetness index	30 m
	Profile curvature	30 m
	Plan curvature	30 m
Geolithology	Distance to an active fault	30 m	The Third Survey and Mapping Institute of Guizhou Province
	Lithology	30 m
Land use and land cover	Land use and land cover	30 m	Ministry of Natural Resources Data Center
Morphology	Distance to rivers and road	30 m	Openstreetmap
Climate	Rainfall	Vector/0.1°	Meteorological station monitoring data and GPM
Vegetation	NDVI	30 m	Landsat
Human	Coal mining	Vector	The Third Survey and Mapping Institute of Guizhou Province

Figure 2. Influencing factors in Nayong. (a) Elevation; (b) slope; (c) profile curvature; (d) plan curvature; (e) aspect; (f) fault and coal mining; (g) rainfall; (h) road and river; (i) land use and land cover; (j) NDVI; (k) TWI; and (l) lithology.

Figure 2. Influencing factors in Nayong. (a) Elevation; (b) slope; (c) profile curvature; (d) plan curvature; (e) aspect; (f) fault and coal mining; (g) rainfall; (h) road and river; (i) land use and land cover; (j) NDVI; (k) TWI; and (l) lithology.

The NDVI time series was calculated from 2016 to 2021, as shown in Figure 3. The 2018 NDVI mean was the lowest. Therefore, the study chose the 2018–2021 NDVI mean as an input factor in this study. Figure 4 depicts a time series graph of NDVI variations from 2016 to 2021, correlating to landslides (purple) and non-landslides (red). The NDVI trends associated with non-landslides showed greater values than the NDVI associated with landslides. This may imply that places with less vegetation are more vulnerable to landslides. For precipitation data, the Global Precipitation Measurement (GPM) mission (IMERG) satellite precipitation products have a spatial resolution of 0.1° (approximately 10 km), which is difficult for regional landslide susceptibility analysis. Therefore, the topographic data (elevation, slope, slope direction) to establish a mapping between precipitation and topographic factors in the study area using the geo-weighted regression approach [34]. Following, descending sampling was used to generate precipitation data for Nayong Town at a resolution of 1 km. Finally, five meteorological station data were used for validation in the study area.

3. Methods

3.1. Overall Methodology

As shown in Figure 5, the study began by collecting and extracting various landslide features from various data sources, including topography, lithology, land use and land cover, morphology, climate, vegetation, and human activity. These features were essential for understanding the underlying causes of landslides and identifying areas at high risk. The extracted landslide features were then converted into a stacked dataset using band synthesis and resampling to ensure data consistency and accuracy. The study split the landslides dataset (50% for the landslide sample and 50% for the non-landslide sample) into training (80%) and test (20%) sets upon random permutation. Next, the GAMI-net and other algorithms (LR, SVM, and RF) were applied to the stack dataset for model calculation and parameter optimization. The grid search was utilized to optimize the parameters by adjusting the weights and biases of the model to achieve the best performance. After the model calculation and parameter optimization, the model’s performance was evaluated using various metrics, such as accuracy, precision, ROC, and recall. These metrics determine how well the model can predict landslides based on the input data. This evaluation helps to identify any issues or limitations of the model and to improve it if necessary. Finally, the model was interpreted and used to create a landslide susceptibility map (LSM) that can be used to identify areas at high risk of landslides.

3.2. Machine Learning Methods

3.2.1. Logistic Regression (LR)

Multiple logistic regression was widely used to examine LSM [35]. The basic idea behind logistic regression is to use Equations (1) and (2) to predict the likelihood of landslides based on a set of input factors. By training the model on a dataset of landslides and non-landslides, it is possible to use the trained model to predict the likelihood of a slope failure in a new, unseen location. The probability values of landslide susceptibility were calculated by the model for all landslide factors in the study area. In general, the probability value indicated the landslide susceptibility of the area and the closer it was to 1, the greater the landslide susceptibility. The logistic function is described as:

$$f\left(x\right)=a_0+{{\displaystyle \sum }}_{k=1}^n\left(a_k{x}_k\right)$$

(1)

$$\mathrm{P} =\frac{1}{1+\left(ex{p}^{-f\left(x\right)}\right)}$$

(2)

where $\mathrm{P}$ is the probability of landslide occurrence and $f\left(x\right)$ a linear combination of casual $x_k$ factors.

3.2.2. Random Forest (RF)

Random forests are a machine-learning algorithm that can predict the likelihood of a landslide occurring in a particular area [36,37]. The algorithm creates many decision trees, each of which is trained on a random subset of the data. The trees are then combined to form a “forest”, and the predictions of the individual trees are averaged to produce a final prediction. To use random forests for landslide susceptibility analysis, a dataset containing information about the characteristics of an area, such as soil type, slope angle, and vegetation cover, is needed. The algorithm is then trained on this dataset to learn the relationship between the predictor variables and the occurrence of landslides in the area. Once the model has been trained, it can predict the likelihood of a landslide occurring in a new area by inputting the relevant predictor variables. The model’s output is a probability, indicating the likelihood of a landslide in the given area.

3.2.3. Support Vector Machines (SVM)

Support Vector Machines (SVM) is a supervised learning algorithm that can be adopted for classification and regression tasks [38,39]. In the context of landslide susceptibility, SVM can classify areas into different susceptibility classes (e.g., low, moderate, and high susceptibility) based on a set of input features or predictors.

SVM works by locating the hyperplane in a high-dimensional space that best divides the various classes. The hyperplane is chosen to maximize the distance between the hyperplane and the nearest data points from each class (called the support vectors). This distance is known as the margin. The larger the margin, the lower the generalization error and the better the classifier will perform on unseen data.

The SVM algorithm involves solving a convex optimization problem to find the hyperplane that maximally separates the classes. The optimization problem can be written as:

minimize w, b, ξ

subject to:

$$\begin{gathered} y_i\left({\mathrm{W}}^{\mathrm{T}}{x}_i+b\right)\ge 1-{ \mathsf{\xi }}_i ,\mathrm{i} =1, 2, \dots , \mathrm{n} \\ {\mathsf{\xi }}_i\ge 0, \mathrm{i} =1, 2, \dots , \mathrm{n} \end{gathered}$$

(3)

where w and b are the weight and bias of the hyperplane, respectively, $x_i$ is the input feature vector, $y_i$ is the class label (0 or 1), and ${\mathsf{\xi }}_i$ is the slack variable that allows for misclassification. The optimization problem is subject to the constraints that the distance between the hyperplane and the nearest data points from each class (the support vectors) is at least 1 and the slack variables are non-negative.

In the case of non-linearly separable data, SVM can still find a good hyperplane by introducing additional dimensions through the use of kernels. A kernel is a function that transforms the input data into a higher-dimensional space where the data may be more easily separable.

3.2.4. GAMI-Net

The generalized additive model (GAM) is a semi-parametric smoothing model that provides a valuable tool for investigating factor interaction by Equation (4) [40].

$$g(E(y|x))=\mu +{\displaystyle \sum }_{j\in S_1}^M{h}_j\left(x_j\right)+\epsilon$$

(4)

where $h_j\left(x_j\right)$ is a smooth function, and j denotes various factors. The $\mu$ indicates the predictor variable category relationship that is unaffected by nonlinear transformation, while the $\epsilon$ represents the residual.

In a traditional GAM network, a feature was made to have monotonically increasing or decreasing properties on the target value by setting constraints and penalty functions or adjusting the number of spline functions [41]. The GAM model are more explanatory than the other black box ML models; they could not be represented as a single function describing the estimated relationship between independent and dependent variables. The GAM model does not consider the factors’ interactions, and the addable spline function does not have sufficient predictive accuracy [42].

To solve the shortcomings of the GAM model, GAMI-net considered factor interactions and replaced the spline function in GAM with an ANN network to improve prediction accuracy. Under the additive model, it was found that it explored and explained a particular variable by being able to control for other variables. As related research has demonstrated, the explanatory nature of ANNs can be achieved under appropriate constraints. This enhancement can explain neural networks with an excellent balance between prediction performance and model interpretation.

$$g(E(y|x))=\mathsf{\mu }+{\displaystyle \sum }_{j\in S_1}^M{h}_j\left(x_j\right)+{\displaystyle \sum }_{\left(j,k\right)\in S_2}^N{f}_{jk}\left(x_j,x_k\right)$$

(5)

According to Figure 6, the GAMI-net consisted of the individual sub-networks, interaction sub-networks, and three constraints. The all sub-networks estimated the individual and interaction factors. At the same time, the three constraints ensured interpretability. The sparsity constraint ranks the features according to their importance, retaining only the more critical feature ridge functions. The heredity constraint ensures that the interaction feature pair must have at least one feature from a vital feature retained in the individual feature sub-network module. The marginal clarity constraint prevents individual features and their corresponding interaction features from absorbing each other, making the model more stable. According to the GAMI-net structure, the hyperparameters are learning rate, batch size, epoch, the maximum number of features, and crossover features. In this study, the learning rate and batch size were set by default, with learning rates of 0.001 and the batch size of 200. The training process was divided into seven stages as follows (

Table 2. The iterative training process of GAMI-net.

Step	Training Process
1	Training all individual factors sub-networks;
2	Selecting the most important $x_i$ individual factors sub-networks based on factors importance;
3	Fine-tuning the individual factors sub-network modules;
4	Selecting interaction pairwise factors based on the ranking of interaction factors;
5	Training the interaction factors sub-network modules;
6	Selecting the top $interactio{n}_i$ significant interaction factors sub-networks based on feature importance;
7	Fine-tuning the interaction factors sub-network modules.

). First, the individual factors sub-networks were trained iteratively for 1000 epochs, and then the interaction factors subnetwork module was trained for 1000 epochs. In each module, the features were ranked according to their importance. Then, the unimportant tail features were removed according to the magnitude of the loss function value and then trained iteratively for 500 epochs. In this paper, the maximum number of individual features was (4, 8, 12, 16, 20), and the maximum number of interaction features was (0, 6, 9, 12, 15, 20) in the grid search.

3.3. Explainability of GAMI-Net

3.3.1. Importance Ratio

The importance ratio (IR) values are used to measure each factor contribution to the overall forecast. The definition of IR referenced the Sobol index, the main difference being that the Sobol index was generated under the premise that all variables were independent and uniformly distributed, whereas IR was based on the empirical distribution of the predicted variables [43,44]. Each factor had an importance score based on the model, and the important ratio $IR\left(j\right)$ of that factor can be calculated using Equation (6). Correspondingly, the importance of each pairwise factor $IR\left(j,k\right)$ could be obtained from Equation (7).

$$IR\left(j\right)=D\left(h_j\right)/T$$

(6)

$$IR\left(j,k\right)=D\left(f_{jk}\right)/T$$

(7)

where $D\left(h_j\right)$ and $D\left(f_{jk}\right)$

$$T={{\displaystyle \sum }}_{j\in S_1}D\left(h_j\right)+{{\displaystyle \sum }}_{\left(j,k\right)\in S_2}D\left(f_{jk}\right)$$

(8)

In addition to the IR factors, the ridge function was used to explain the relationship between one or two specific factors and the landslides. The $h_j\left(x_j\right)$ value drew one-dimensional linear graphs for continuous factors and the bar charts for categorical factors to show the factor-landslides relationship. Each of the pairwise factors $f_{jk}\left(x_j,x_k\right)$ can be visualized using a 2D heat map showing the effect of the factor on the slippage.

3.3.2. Shapley Additive Explanations (SHAP)

For a trained complex machine learning model, SHAP interprets the model by calculating the Shapley values of the features to measure the importance of the features [45,46]. The Shapley value can be calculated for each feature component of each sample, and the Shapley values of the sample feature components are additive. The magnitude of the Shapley value of a sample feature component tells how much the value taken by the sample on each feature affects the prediction result for that sample [47]. The overall feature importance can be obtained by calculating the absolute value of the Shapley values of all samples on each feature to find out the magnitude of each feature’s influence on the prediction results of the complex model by averaging the overall importance.

With the SHAP method, GAMI-net may directly extract the value of each additive component, such as the primary factor or the pairwise interaction factor, for every sample x. The sequence of these marginal effects can help to clarify how the input x’s choices should be understood. The 1D line plots (or bar charts) and 2D heat maps can be used to statistically analyze how sensitive the forecast is to modest changes in the explanatory factors.

4. Results

4.1. Model Assessments for RF, SVM, and GAMI-Net

The variation of model accuracy (ROC) for the different maximum number of feature parameters (individual and interaction) in the parametric grid search was shown in Figure 7. As the maximum number of features increased, the model accuracy improved. The highest accuracy was obtained when the maximum number of features was 16 and 15 for single and interactive features, respectively. After this point, the model accuracy decreased, probably due to the increase in error in additive models with too many parameters, indicating the need for fine-tuning.

Table 3. Evaluation of each model.

Algorithm	Accuracy	Precision	F1 Score	Recall	AUC ROC
LR	0.6189	0.6160	0.5754	0.5400	0.6482
SVM	0.6708	0.7013	0.6129	0.5448	0.7530
RF	0.8110	0.8222	0.7897	0.7581	0.9027
GAMI-net	0.8730	0.8648	0.8683	0.8747	0.9442

illustrated the evaluation metrics obtained by predicting 510 samples (20% of the total sample was the test set) within the study area using each algorithm, with LR generally having the weakest algorithm performance and the SVM, RF, and GAMI-net algorithms having better performance. Based on the ROC curves, the algorithms were ranked from smallest to most significant in terms of accuracy: LR, SVM, RF, and GAMI-net.

As shown in Figure 8, the GAMI-net method provided comparable performance to the RF methods, with prediction accuracies all greater than 0.8, F1 scores of 0.7897 and 0.8683, respectively, and ROC values all greater than 0.9. The GAMI-net model outperformed LR and SVM and slightly outperformed RF in predicting landslide susceptibility in the study area, demonstrating that the algorithm could achieve good results.

Figure 9 shows the training and validation losses of the GAMI-net with 25 interaction effects considered. In the first stage, the loss rate does not change significantly when only the main effects are considered, indicating that landslide susceptibility cannot be effectively predicted using only the main effects. According to Table 3, as paired interactions were added to the network, the loss rate reduced dramatically, demonstrating the necessity to add pairwise interactions to the GAMI-net. There was a significant surge in training and validation loss at the start of the fine-tuning phase, which coincided with a fine-tuning operation on the paired factors. Figure 10 also depicts the validation loss for establishing the ideal number of main effects and pairwise interactions. The number of main effects and pairwise interactions included is indicated on the left and right image x-axes, respectively. The red star symbols represent the number of main features or pairwise interactions that should be used. According to the findings, the GAMI-net had 15 main features and 13 pairwise interactions. The rise in the primary and pair effects beyond this period had no significant influence on the loss rate.

4.2. Comparing Landslides Susceptibility Map

Landslide susceptibility maps can be generated using a variety of methods, including support vector machines (SVM), random forests, and GAMI-net. Each of these methods has its own advantages and disadvantages, and the most appropriate method will depend on the specific needs of the application.

The map divides Nayong into several categories of landslide susceptibility, including high, moderate, low, and very low. High-susceptibility areas are those with a high probability of landslides occurring, while low-susceptibility areas have a low probability of landslides.

In Figure 11, the map showed that the highest susceptibility areas are found in the coal mining area, particularly in the central and southern part of the study area. As coal is extracted from underground mines, the void left behind can cause the overlying rock and soil to collapse, leading to landslides. This is especially true in Karst areas, where the rock is naturally porous and prone to landslides. In addition, the southern region is characterized by a highly weathered layer of loose dolomite, which is subject to intense erosion due to the impact of precipitation. This erosion process significantly disrupts the mechanical equilibrium of the otherwise fragile grass and irrigation-dominated system in a short period, leading to an increase in the regional susceptibility to landslides. This underscores the need for careful management of erosion and landslides in this region. Moderate susceptibility areas are found in the northern and central parts of the state, including the Mazong Hill. These areas are characterized by a mix of steep and gentle slopes, as well naturally porous rock and prone-to-collapse slope, along with some coal mining operations. They are also prone to mining and heavy rainfall, which can increase the risk of landslides. Low-susceptibility areas are found in the western part of the state. These areas are characterized by gentle slopes and sandy soils. They are generally less prone to landslides, although they can still be affected by earthquakes and other natural disasters.

4.3. Interpretative Results for GAMI-Net and Shapley Additive Explanations

Figure 12 showed the importance estimates for all impact factor variables for the GAMI-net model. The results for the study area indicated that rainfall, mining areas, lithology, elevation, rivers, and faults are the main factors. In brief, human activity and climate had a strong influence on landslides.

The global interpretation of the GAMI-net model factors was shown in Figure 13. With the introduction of 13 main effects and 25 interaction effects, it is found that rainfall, mine, river, failure, and lithology are the factors that account for 6.8%, 5.8%, 2.8%, 1.8%, and 1.7% of the landslide susceptibility factors, respectively. In terms of interaction effects, lithology rainfall, mining, rainfall, and elevation-rainfall account for 32.8% and 20.7% of the landslide susceptibility factors. The other susceptibility factors and their interaction effects were below 1%. In relative terms, the importance of the main effect is 18.9%, and the importance of the pairwise effect is 81.1%.

Firstly, rainfall was the most dominant influence factor. Rainfall causes soil destabilization and is highly susceptible to landslides. Extreme rainfall (class 1) significantly impacts landslides, with subsequent impacts on landslides decreasing as rainfall values decrease. Rainfall can cause landslides in several ways. For the distance factor of the mining area, the closer the mining area (the smaller the classification value in the graph), the higher the model score, indicating a higher landslide susceptibility factor, which suggests that mining activities are the main factor for landslide susceptibility in the area. Certain types of rock, such as shale (class 13) and clay (class 7), have high water retention capabilities due to their high porosity and low permeability. When these rocks become saturated with water, they can become unstable and prone to landslides. In contrast, rocks with low water retention capabilities, such as sandstone (class 3) and granite (class 11), are less prone to landslides because they are able to drain excess water more effectively. Earthquake-induced landslides are less common in the region due to few earthquakes and low seismic intensities. In addition, faults can create cracks and other types of pressure in the ground, which can increase the likelihood of landslides, and this probably will increase under the influence of mining and rainfall.

Secondly, in terms of interaction effects, the interaction between rainfall and lithology, mining, and elevation are the main factors, and the vulnerability of fragile rocks to landslides under rainfall conditions is evident. For the Nayong area, long-term mining activities have increased the fragile surface of the mountain. Landslide susceptibility was then significantly increased under external conditions such as rainfall.

Overall, the research measured the influence of various factors on landslide susceptibility using 1D/2D heat maps. For the study area, rainfall, mining, and lithology are the most significant influences.

4.4. Factorial Interpretability Analysis of the Zongling Landslides

The location of the Karst landslide in Zongling Town, Guizhou Province, China, was shown in Figure 14. It is a typical karst erosion landscape in the Guizhou plateau. Figure 14 is an optical image and InSAR deformation of Zongling, with the black polygon showing the landslide area. The InSAR surface deformation monitoring results were from our previous study [13,48]. Since the 1980s, Zongling Town has gradually become a central coal mining area, and a large amount of coal mining has created many landslide hazards. In particular, there is a cluster of landslides on the southern side of the Mazongling Mountains, running from west to east for a length of 5.6 km, and these landslides seriously threaten the people’s safety. The rear edge of the slope is an almost vertical escarpment over 30 m on top and below. There are numerous ground fissures on the summit, ranging from several hundred meters to several kilometers in length, with fissures ranging from nearly 2 to 6 m in width. The study used Zongling landslides as an example for interpretable analysis based on the landslide susceptibility results above. As shown in Figure 14, points A and B are located in the more deformed areas of the Zongling landslide. The results of the interpretable model analysis in Figure 15 reveal that rainfall, mining, lithology, and their coupling were the key factors in the landslides at these two places, which is congruent with the ground research results. As shown in Figure 16, mining causes stress damage to the Karst hills in the area; lack of vegetation, extreme climatic background, and rock desertification further exacerbate landslide susceptibility and threaten the safety of the roads and people below. It suggests that our interpretable results are credible and can provide some basis for disaster mitigation and resilience.

5. Discussion

In this study, we used a combination of GIS, remote sensing, and machine learning to analyze landslide susceptibility in Nayong County, Guizhou Province. It analyzed complex functional relationships through ANN networks of main effect factors and pairwise interaction factors and considered the coupling between the landslide factors, increasing the model’s interpretability. In addition, the relationship between the influence of landslide factors on landslide susceptibility can be easily assessed using 1D line/bar diagrams and 2D heat maps. Our results showed that landslide susceptibility in the region is significantly correlated with several landslide influence factors, including rainfall, slope, and lithology. By introducing interaction factors, we also found that coupling several factors had significant effects on landslide susceptibility, such as rainfall and mining, lithology and rainfall, rainfall and elevation, and elevation and mining.

Validation of the four models showed that GAMI-net was the most reliable method, with the best model performance on the landslide dataset, followed by RF and SVM, while LR had the worst. The result suggests that GAMI-net can provide a strong interpretation while guaranteeing better accuracy. Compared with the four machine-learning models, LR has fast classification speed and easy interpretation of classification results but has a specific error rate in classification decisions; RF has high training efficiency and high practicality but ignores the correlation between data. SVM can avoid overfitting but depends on the choice of kernel function; GAMI-net has the highest accuracy, considers factors’ coupling effect, has good interpretability, and quantifies the influence of landslide features.

In addition, our study uses interpretable algorithms to account for the changing influence of these factors on landslide susceptibility in the region. This method adds credibility to the landslide susceptibility model and provides a method for quantitatively assessing the influence of natural and human factors on landslides. Ranking rainfall and mining as more influential features in landslide prediction was in line with a series of regional studies [25,49,50]. Tao and Shi et al. used the Jinhaihu landslide event in Bijie as an example of a physical experiment and field investigation to show that lithology, construction, and rainfall (snow) were the main factors contributing to the landslide [49]. Wang et al. used time-series InSAR techniques to illustrate the persistent influence of mining activities on the Zongling landslide [50]. Unusually, the slope factor was classified as a less influential feature in Nayong, contrary to the findings of Zhao and Masoumi et al. in Iran [51], which may reflect the characteristics and regional effects of landslide events, landslide causal factors, and the dependence of the algorithm from one study area to another. In the present study, although the slope was positively correlated, its effect was small (Figure 13). The Nayong region was Karst landform with steep slopes at high elevations. The slopes in most parts of the region were eligible for landslide occurrence, making slope a relatively low-impact feature in the internal modeling process.

Our results showed that targeted measures, such as rational mining and rock desertification reduction, can effectively reduce landslide risks. In addition, our study could be used to inform land resource managers and policymakers. There have several limitations in the study that should be acknowledged. First, our analysis is based on remotely sensed data, which may only capture some relevant variables affecting landslides, such as cohesion and angle of internal friction. In addition, our study area is limited to mining areas and karst landscapes, so extension to other areas with different climatic and land use patterns will require corresponding changes.

6. Conclusions

LSM is an essential tool for landslide hazard reduction. Existing LSM research has mainly focused on machine-learning algorithms, with less interpretable research on LSM. Therefore, the GAMI-net method was applied as an interpretable algorithm and three traditional machine-learning models (LR, SVM, RF) to generate landslide sensitivity maps for Nayong region of Guizhou Province.

Firstly, the performance of the GAMI-net model and the three traditional ML models were examined and compared against several evaluation metrics. The results showed that GAMI-net outperformed the other three models (accuracy = 87.3%, ROC = 0.944). GAMI-net shows excellent potential for susceptibility assessment and interpretation, and it would be advantageous to explore its predictive capabilities in landslide susceptibility assessment.

Second, the LSM shows that most of the high-susceptibility study areas are located in mining and heavy rainfall areas, which are affected by erosion processes and rock desertification. The low susceptibility areas are the main areas of gentle slopes and lithological stability.

Thirdly, GAMI-net allows for an assessment of the significance of the characteristics of the sample (global and local interpretation). The findings suggest rainfall, mining, and lithology are the main causal factors for local landslides. In addition, some interaction factors such as rainfall and mining, lithology and rainfall, rainfall and elevation, and elevation and mining are also worth noting.

Overall, GAMI-net could be used to develop LSM and provide clear explanations for why the model made that prediction. This could help decision-makers better understand the risks of landslides in their area and take appropriate actions to prevent or mitigate them. The results of this work could be useful for decision-makers and planners in planning land use in areas prone to landslides.