Rapid Landslide Extraction from High-Resolution Remote Sensing Images Using SHAP-OPT-XGBoost

Abstract

Landslides, the second largest geological hazard after earthquakes, result in significant loss of life and property. Extracting landslide information quickly and accurately is the basis of landslide disaster prevention. Fengjie County, Chongqing, China, is a typical landslide-prone area in the Three Gorges Reservoir Area. In this study, we newly integrate Shapley Additive Explanation (SHAP) and Optuna (OPT) hyperparameter tuning into four basic machine learning algorithms: Gradient Boosting Decision Tree (GBDT), Extreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (LightGBM), and Additive Boosting (AdaBoost). We construct four new models (SHAP-OPT-GBDT, SHAP-OPT-XGBoost, SHAP-OPT-LightGBM, and SHAP-OPT-AdaBoost) and apply the four new models to landslide extraction for the first time. Firstly, high-resolution remote sensing images were preprocessed, landslide and non-landslide samples were constructed, and an initial feature set with 48 features was built. Secondly, SHAP was used to select features with significant contributions, and the important features were selected. Finally, Optuna, the Bayesian optimization technique, was utilized to automatically select the basic models’ best hyperparameters. The experimental results show that the accuracy (ACC) of these four SHAP-OPT models was above 92% and the training time was less than 1.3 s using mediocre computational hardware. Furthermore, SHAP-OPT-XGBoost achieved the highest accuracy (96.26%). Landslide distribution information in Fengjie County from 2013 to 2020 can be extracted by SHAP-OPT-XGBoost accurately and quickly.

1. Introduction

Landslides are a type of slope instability phenomenon caused by the combined effects of internal factors such as terrain, geological structure, lithology, and external factors such as meteorology, hydrology, and weathering, as well as human activities like engineering. As one of the world’s most hazardous geological disasters, landslides are characterized by their suddenness, high concealment, and destructive power, and can cause casualties, infrastructure damage, farmland destruction, building collapse, and other significant losses [1,2]. Traditional methods for landslide detection and mapping mainly involve field surveys and remote sensing visual interpretation [3]. However, field surveys have limitations, including a low efficiency, high risk, and limited view, particularly in areas with a large number and widespread distribution of landslides, which require more time, workforce, and financial costs [4]. Although visual interpretation based on remote sensing technology can overcome unfavorable conditions such as terrains, transportation, and weather, most landslide identification techniques still rely on human–computer interactions. The interpretation expertise of the personnel involved is critical, and the quality of the interpretation result is greatly influenced by the experience and subjectivity of experts [5]. Furthermore, when dealing with large-scale regions, there are still many drawbacks, such as a low identification efficiency, long duration, weak specificity, and susceptibility to omission [6]. Therefore, traditional methods are easily constrained in landslide investigations. By combining machine learning methods with high-resolution remote sensing images, the efficiency of landslide identification can be improved while reducing subjective bias caused by reliance on human experts [7]. The sample application and spatial generalization ability of machine learning methods are outstanding, making them one of the current mainstream technologies for automatic or semi-automatic extraction of landslide spatial information [8].

Most scholars use remote sensing satellite images combined with classical algorithms such as random forests (RFs) [9,10], support vector machines (SVMs) [11,12], convolutional neural networks (CNNs) [13,14], and decision trees (DTs) [15] in machine learning to build classifiers for landslide identification and extraction. Although these algorithms have solved the problem of expressing high-dimensional nonlinear relationships, they still have some shortcomings, such as poor interpretability of RF, difficulty controlling the operation of the model’s internals, and poor handling of missing data by DT, which can lead to overfitting. Boosting ensemble algorithms, such as GBDT, XGBoost, LightGBM, and AdaBoost, have been proven to outperform other commonly used machine learning methods, including RF and SVM, in terms of classification performance, robustness to training data quality, running speed, sample size requirements, and prediction accuracy. XGBoost, especially, has emerged as a powerful and efficient algorithm that is popular for both regression and classification tasks. This ensemble learning framework is based on gradient-boosting trees and was designed by Dr. Tianqi Chen from the University of Washington in 2016 [16]. XGBoost boasts a strong robustness, superior performance, and good processing speed and accuracy for low-dimensional data. It excels when there are limited training samples, a short training time, and a lack of tuning knowledge, making it ideal for extracting landslides rapidly.

SHAP and Optuna are two powerful analytical tools. SHAP adopts game-theoretic methods to calculate the marginal contribution of each feature to the model’s output and assign a specific predictive importance value to each feature, ensuring good global and local interpretability [17]. Therefore, SHAP has a wide range of applications in various fields [18]. Hyperparameter optimization is the process of finding the optimal set of hyperparameters for a machine learning model. Hyperparameters can greatly impact the performance of the model, and finding the optimal values can be a time-consuming and complex task. Optuna is a Python library used for hyperparameter optimization, which uses Bayesian optimization and approximate target algorithms to rapidly and reliably identify the optimal hyperparameter configuration for improved model performance [19]. By incorporating the SHAP algorithm for feature selection, the speed and quality of feature selection can be enhanced. Moreover, the introduction of the Optuna learning framework for Bayesian hyperparameter automatic optimization can mitigate the base model’s overfitting and low generalization ability, achieving a faster and more cost-effective way of building models with better predictive performance.

Fengjie County is located in the hinterland of the Three Gorges Reservoir Area and is a region prone to landslide disasters in southwest China [20]. We focused on Fengjie, Chongqing, and employed multiple data sources including high-resolution remote sensing satellite image data such as GF-1, GF-2, GF-6, and ZY3-01, an advanced land observing satellite (ALOS) digital elevation model (DEM) with a 12.5 m resolution, and a historical landslide disaster inventory. We newly integrated the explanatory framework SHAP and Optuna hyperparameter tuning into GBDT, XGBoost, LightGBM, and AdaBoost, and constructed SHAP-OPT-GBDT, SHAP-OPT-XGBoost, SHAP-OPT-LightGBM, and SHAP-OPT-AdaBoost rapid landslide extraction models for the first time (hereafter referred to as SHAP-OPT-Models throughout the rest of the paper). We compared the performance and efficiency of different SHAP-OPT-Models and quantitatively analyzed their temporal and spatial distribution in the study area for the years from 2013 to 2020.

The remaining sections of this paper are organized as follows: Section 2 presents an overview of the study area and the data sources used. Section 3 provides a comprehensive description of the experimental methods and procedures, encompassing high-resolution image preprocessing, the conception of the core algorithms, and the development of landslide extraction models based on SHAP and Optuna. Section 4 and Section 5 presents the experimental results and discussions. Finally, Section 6 summarizes the key findings and conclusions derived from the study.

2. Study Area and Datasets

2.1. Study Area

Fengjie County (Figure 1) is situated in Chongqing City, China. It is located in the eastern Sichuan Basin and at the heart of the Three Gorges Reservoir Area (between latitude 30.02°N to 30.93°N and longitude 110.30°E to 110.87°E), covering an area of 4098 km² and belonging to the subtropical humid monsoon climate, with an average annual precipitation of 1132 mm.

The geological environment in the research area is extremely complex, mainly consisting of anticlines and syncline valleys, which roughly correspond to the structural lines. The relief degree of the land surface (RDLS) is large, with the lowest altitude at the Qutang Gorge, which was 86 m before the Three Gorges Project reservoir was filled, and the highest altitude at Maoerliang, which is 2123 m above sea level. The area is symmetrical surrounding the Yangtze River, with high terrain in the north and south, and gentle terrain in the central and western regions. The Yangtze River flows through Fengjie County for 43 km, with a dense distribution of tributaries. According to Fengjie’s hydrological data, there are a total of 17 river basins in the research area, with an area exceeding 50 km². The highest water level can reach 125.13 m, with a runoff volume of 27.95 × 108 m³ flowing into the Yangtze River. The average annual runoff of the tributaries is 420.8 mm. During the rainy season, when the Yangtze River floods, the flow rate can reach a maximum of 69,500 m³/s, and heavy summer rainfall can trigger a large number of landslides, making it a typical landslide-prone area.

2.2. Data Collection and Preprocessing

Due to the long time span and broad scope of this research, it was necessary to consider the hardware constraints of computers while achieving high-precision extraction of landslides. The high-resolution satellite imagery from the Gaofen series satellite was found to be well-matched to the study’s needs. However, due to incomplete coverage of Fengjie County during the study years, ZY3-01 satellite imagery was additionally selected to construct a complete long time-series image data source ranging from 2013 to 2020. In total, 71 images with a raw data size of approximately 98.5 GB were selected. Specifically, this imagery came from the Chongqing Data and Application Center of the Chinese High-Resolution Earth Observation System, and most of the imaging took place during winter when precipitation was low and vegetation was in a non-growing state. For specific details, please refer to Appendix A,

Table A1. List of high-resolution historical remote sensing image data.

Year	Satellite	Sensor	Spatial Resolution/m		Cloud Cover/%	Date of Acquisition	Precipitation Situation	Vegetation Growth Status
			Panchromatic	Multispectral
2013	ZY3-01	MUX NAD	2.1	5.8	2	2013/2/10	less rain	non-growing
			2.1	5.8	9	2013/3/26	less rain	growing
			2.1	5.8	1	2013/8/11	rainy	growing
			2.1	5.8	9	2013/8/11	rainy	growing
			2.1	5.8	0	2013/12/2	less rain	non-growing
			2.1	5.8	0	2013/12/2	less rain	non-growing
2014	GF-1	PMS1	2	8	0	2014/3/26	less rain	growing
		PMS1	2	8	0	2014/3/26	less rain	growing
		PMS1	2	8	7	2014/7/27	rainy	growing
		PMS1	2	8	2	2014/7/27	rainy	growing
		PMS1	2	8	5	2014/7/27	rainy	growing
		PMS1	2	8	14	2014/7/27	rainy	growing
		PMS2	2	8	2	2014/7/27	rainy	growing
		PMS2	2	8	18	2014/7/27	rainy	growing
		PMS2	2	8	12	2014/7/27	rainy	growing
		PMS1	2	8	5	2014/7/31	rainy	growing
		PMS2	2	8	2	2014/7/31	rainy	growing
2015	GF-1	PMS2	2	8	4	2015/2/17	less rain	non-growing
		PMS1	2	8	0	2015/3/30	less rain	growing
		PMS1	2	8	0	2015/3/30	less rain	growing
		PMS2	2	8	0	2015/3/30	less rain	growing
		PMS2	2	8	0	2015/3/30	less rain	growing
		PMS2	2	8	31	2015/3/30	less rain	growing
		PMS1	2	8	1	2015/5/14	rainy	growing
		PMS1	2	8	3	2015/5/14	rainy	growing
		PMS1	2	8	22	2015/5/14	rainy	growing
		PMS1	2	8	9	2015/8/16	rainy	growing
		PMS1	2	8	22	2015/8/16	rainy	growing
		PMS1	2	8	33	2015/8/16	rainy	growing
2016	GF-1	PMS1	2	8	0	2016/8/19	rainy	growing
		PMS1	2	8	1	2016/8/19	rainy	growing
		PMS1	2	8	0	2016/9/17	rainy	growing
		PMS1	2	8	2	2016/9/17	rainy	growing
		PMS1	2	8	2	2016/9/17	rainy	growing
		PMS1	2	8	15	2016/9/17	rainy	growing
		PMS2	2	8	0	2016/9/17	rainy	growing
		PMS2	2	8	1	2016/9/17	rainy	growing
		PMS2	2	8	0	2016/12/4	less rain	non-growing
2017	ZY3-01	MUX NAD	2.1	5.8	3	2017/10/28	less rain	non-growing
			2.1	5.8	19	2017/10/28	less rain	non-growing
	GF-1	PMS1	2	8	1	2017/5/13	rainy	growing
		PMS1	2	8	1	2017/5/13	rainy	growing
		PMS2	2	8	0	2017/5/13	rainy	growing
		PMS1	2	8	1	2017/11/5	less rain	non-growing
		PMS2	2	8	0	2017/11/5	less rain	non-growing
		PMS2	2	8	0	2017/11/5	less rain	non-growing
		PMS1	2	8	0	2017/11/9	less rain	non-growing
		PMS2	2	8	0	2017/11/9	less rain	non-growing
		PMS2	2	8	5	2017/11/9	less rain	non-growing
2018	ZY3-01	MUX NAD	2.1	5.8	1	2018/8/24	rainy	growing
			2.1	5.8	2	2018/8/24	rainy	growing
			2.1	5.8	0	2018/8/29	rainy	growing
			2.1	5.8	0	2018/8/29	rainy	growing
	GF-1	PMS1	2	8	0	2018/1/14	less rain	non-growing
	GF-2	PMS2	2	8	0	2018/1/22	less rain	non-growing
		PMS2	2	8	0	2018/1/22	less rain	non-growing
2019	ZY3-01	MUX NAD	2.1	5.8	0	2019/8/13	rainy	growing
			2.1	5.8	3	2019/8/13	rainy	growing
			2.1	5.8	32	2019/8/13	rainy	growing
	GF-6	PMS	2	8	1	2019/11/3	less rain	non-growing
			2	8	7	2019/11/3	less rain	non-growing
2020	GF-1	PMS1	2	8	1	2020/1/30	less rain	non-growing
		PMS1	2	8	1	2020/1/30	less rain	non-growing
		PMS1	2	8	1	2020/1/30	less rain	non-growing
		PMS1	2	8	5	2020/1/30	less rain	non-growing
		PMS2	2	8	3	2020/1/30	less rain	non-growing
		PMS2	2	8	13	2020/1/30	less rain	non-growing
		PMS2	2	8	0	2020/11/8	less rain	non-growing
		PMS2	2	8	0	2020/11/8	less rain	non-growing
		PMS2	2	8	0	2020/11/8	less rain	non-growing
		PMS2	2	8	0	2020/11/8	less rain	non-growing

for the data list. In addition, the ALOS digital elevation model (DEM) data and a landslide inventory map were also used. The DEM data were from Chongqing Public Service Platform of Geographic Information. The landslide inventory map was obtained through historical landslide hazard inventories and field surveys, providing the location information of the landslide.

Preprocessing of the raw remote sensing images is an important prerequisite for the extraction of landslide temporal and spatial distribution information. Preprocessing mainly includes radiometric calibration, atmospheric correction, geometric correction, orthorectification, image fusion, unification of the coordinate system, image mosaic, image cropping, and color balance, which can ensure the normativity, uniformity, and validity of the input data. In addition, it is necessary to extract the stratigraphic lithology and structure information from the geological map through artificial vectorization.

3. Methodology

In this paper, eCognition, ArcGIS, PyCharm, and Anaconda software were used as the key technology platforms. Firstly, comprehensive consideration was given to the geometry, texture, topography, and spectral factors, resulting in the construction of datasets comprising landslide and non-landslide samples, along with an initial feature set consisting of 48 features. High-resolution remote sensing images were utilized as the primary data, integrating feature selection based on the SHAP and automatic hyperparameter optimization utilizing the Optuna framework into foundational machine learning algorithms including GBDT, XGBoost, LightGBM, and AdaBoost. The optimal hyperparameters were obtained iteratively to create SHAP-OPT-Models. By comparing their performance, a rapid landslide extraction model was constructed to extract the spatial distribution information of landslides from 2013 to 2020 in the study area. The creation and optimization of SHAP-OPT-Models were completed through Python programming and secondary development. The primary third-party libraries involved in this process included XGBoost, LightGBM, AdaBoost, SHAP, and Optuna. The experiment was completed on the Windows 10 platform, with Intel(R) Core(TM) i9-10900F CPU, 32 GB running memory and a NVIDIA GeForce RTX 2060 graphics card. An overview of this study is shown in Figure 2.

3.1. Image Segmentation

This paper utilized high-resolution satellite imagery data with a spatial resolution of 2 m, resulting in a substantial data volume. However, computational constraints of the hardware led to extended processing times for the landslide extraction model [21]. To address this challenge, a regular grid approach was implemented to partition the image data of the study area into multiple tile datasets. These tiles were processed in parallel using synchronized image segmentation, enabling the generation of geographic objects. Subsequently, the resulting datasets were integrated to facilitate a comprehensive analysis of the area. The fractal net evolution approach (FNEA) is a widely adopted image multiresolution segmentation algorithm extensively utilized in earth sciences and implemented using the commercial software eCognition [22,23]. Its objective is to generate non-overlapping geographical object units by dividing images based on the heterogeneity index between adjacent pixels, adhering to the principle of minimum heterogeneity [24,25]. To overcome the challenge of varying scales among different surface covers and segmentation scales and to mitigate the subjectivity in selecting the optimal segmentation scale for landslides, this study used the method of estimating the scale parameter 2 (ESP2), which could automatically optimize the scale parameter for landslide segmentation [26,27]. During the multi-scale segmentation process, ESP2 computed the rates of change of the local variance (ROC-LV) for image objects under different segmentation scale parameters. The identification of the optimal segmentation scale for delineating object features was accomplished by observing the peak in the ROC-LV curve [28].

3.2. Construction Samples and Initial Features

3.2.1. Constructing Samples

A rational and representative sample set is crucial for optimizing and evaluating the accuracy of landslide extraction models. Due to the inherent spectral differences between multi-temporal images, the features of the same landslide may not be entirely consistent across different periods. Therefore, it is necessary to select samples from different regions as well as from the same location at different times, including both landslide and non-landslide samples [29]. By combining actual land cover and landslide distribution data supplemented by field survey results with a primary focus on manual interpretation, a total of 1844 segmented objects were selected from images of various years to form a sample set, as shown in Figure 3. Among them, there were 460 positive samples representing landslides and 1384 negative samples representing different land cover types, including cultivated land, bare land, buildings, shadows, water bodies, woodlands, and roads.

3.2.2. Building Initial Feature Set

The spectral, geometric, and textural characteristics of high-resolution remote sensing images provide crucial information sources and act as technical foundations in the process of landslide extraction [30,31,32]. In landslide-prone regions, vegetation degradation and exposed bedrock phenomena are commonly observed, demonstrating a strong correlation between the spatial distribution of landslides and terrain features. Furthermore, terrain features, as well as normalized difference vegetation index (NDVI) and normalized difference soil index (NDSI), also serve as significant reference indicators [33]. The spectral, geometric, and textural features were calculated based on the objects after image segmentation, and the terrain features were extracted from the DEM [34]. By comprehensively considering them, a set of 48 initial features was constructed, as illustrated in Appendix B 

Table A2. A comprehensive set of 48 initial features.

Type	Feature	Feature Meaning
spectrum	Mean I (i = Red, Green, Blue, Nir)	Band means, mean for red, green, blue, near-infrared bands.
	Standard deviation i (i = Red, Green, Blue, Nir)	Standard deviation of the object in the red, green, blue, and near-infrared bands.
	Brightness	Average brightness value of all bands in the image.
	Max. diff.	Maximum spectral difference value among all image bands.
Texture	GLCM Homogeneity (all dir.)	Homogeneity of grey-level co-occurrence Matrix (GLCM): Measures the local gray-level uniformity of the image.
	GLCM Contrast (all dir.)	Contrast of GLCM: Measures the total amount of local variation in the image.
	GLCM Dissimilarity (all dir.)	Dissimilarity of GLCM: Similar to contrast, measures the amount of local variation in the image.
	GLCM Entropy (all dir.)	Entropy of GLCM: Measures the amount of information in the image.
	GLCM Ang. 2nd moment (all dir.)	Second moment of GLCM: Measures the uniformity of the gray-level distribution in the image.
	GLCM Mean (all dir.)	Mean of GLCM: Reflects the regularity and uniformity of gray levels in the image.
	GLCM StdDev (all dir.)	Standard deviation of GLCM: Reflects the deviation between gray-level values and their mean in the image.
	GLCM Correlation (all dir.)	Correlation of GLCM: Reflects the length of extension of certain gray-level values along a certain direction in the image.
	GLDV Ang. 2nd moment (all dir.)	Second moment of GLDV: Measures the local homogeneity of the image.
	GLDV Entropy (all dir.)	Entropy of GLDV: Measures the complexity of the image.
	GLDV Mean(all dir.)	Mean of GLDV: Reflects the regularity and uniformity of gray levels in the image.
	GLDV Contrast (all dir.)	Contrast of GLDV: Measures the total amount of local variation in the image.
Geometric	Area (Pxl)	Area: Number of pixels in the object.
	Border length (Pxl)	Boundary length: Total number of edge pixels in objects shared with other objects.
	Length (Pxl)	Length: Product of the total number of pixels in the object and the aspect ratio of length to width.
	Length/Width	Aspect ratio: Ratio of length to width of the object.
	Volume (Pxl)	Volume: Volume of the object in the image.
	Width (Pxl)	Width: Ratio of the total number of pixels in the object and the aspect ratio of length to width.
	Asymmetry	Asymmetry: Relative length of the object.
	Border index	Boundary index: Indicates the degree of irregularity of the object.
	Compactness	Compactness: Describes the compactness of the object.
	Radius of smallest enclosing ellipse	Minimum radius of the external ellipse: Describes the similarity between the object’s shape and an ellipse.
	Elliptic Fit	Fitting degree of the ellipse: Describes the degree of approximation between the object and a similar-sized ellipse.
	Density	Density: Spatial distribution of pixels in the object.
	Rectangular Fit	Fitting degree of the rectangle: Degree of approximation between the object and a similar-sized rectangle.
	Radius of largest enclosing ellipse	Maximum radius of the internal ellipse: Describes the similarity between the object and an ellipse.
	Roundness	Roundness: Degree of similarity between the object and an ellipse.
	Shape index	Shape index: Smoothness of the object boundary.
Index	NDVI	Normalized difference vegetation index (NDVI): Calculated as (NIR − R)/(NIR + R), where NIR is the near-infrared band and R is the red band.
	NDSI	Normalized difference soil index (NDSI): Calculated as (R − G)/(R + G), where R is the red band and G is the green band.
Terrain	Mean i (i = DEM, Slope, Aspect, Relief)	Mean of terrain features: Average value of elevation, slope, aspect, and relief bands in the image object.
	Standard deviation i (i = DEM, Slope, Aspect, Relief)	Standard deviation of terrain features: Standard deviation of elevation, slope, aspect, and relief bands in the image object.

.

3.3. Basic Machine Learning Model

3.3.1. GBDT

GBDT is a machine-learning model originally proposed by Jerome Friedman in 1999 [35]. This model proved to be a cornerstone for the development of algorithms such as XGBoost and LightGBM. GBDT performs classification tasks using classification and regression trees (CART) and applies softmax functions to map the output classification structure [36]. It trains a set of regression trees serially, and the final prediction is obtained by summing the predictions of all regression trees, resulting in a robust learner.

3.3.2. XGBoost

XGBoost uses CART as the base learner and optimizes the traditional GBDT to achieve ensemble learning of multiple trees, which can be used to solve machine learning problems such as classification and regression [37,38]. The construction process of the XGBoost model is illustrated in Figure 4. Firstly, an initial tree is constructed using the training set for model training, and the residuals of the model prediction and the actual values are obtained. Then, in each iterative process, a tree is added to fit the residuals predicted by the model in the previous iteration until the model learning process is terminated, resulting in an iterative residual tree ensemble composed of multiple tree models [39,40].

The predicted value ${\hat{y}}_i$ is defined as:

$${\hat{y}}_i={{\displaystyle \sum }}_{k=1}^K{f}_k(x_i)$$

(1)

where ${\hat{y}}_i$ represents the final model prediction, K represents the number of classification regression trees, $x_i$ represents the input variables, and $f_k$ represents the model of the kth tree.

The objective function of XGBoost can be written as:

$$O_{\mathit{bj}}={{\displaystyle \sum }}_{i=1}^m{l(y}_i,{\widehat{ y}}_i)+{{\displaystyle \sum }}_{k=1}^K\mathsf{\Omega }(f_k)$$

(2)

where the first term is the training loss function, which measures the difference between the true value $y_i$ and the predicted value. The second term is the regularization term ${\mathsf{\Omega }(f}_k)$, which can be defined as $\mathsf{\Omega }\left(f\right)=\mathsf{\gamma }T+\frac{1}{2}\mathsf{\lambda }\Vert w{\Vert }^2$ to control the model complexity and prevents overfitting.

3.3.3. LightGBM

LightGBM is a machine-learning model introduced in late 2016 by Guolin Ke [41]. Although both use the negative gradient of the loss function and residuals to fit a new decision tree, LightGBM improves upon GBDT in several ways: (1) Using a histogram-based algorithm to discretize data for better decision tree splitting instead of scanning the entire dataset for each feature at every node; (2) Optimizing the decision tree using a leaf-wise strategy to control model complexity; (3) Employing a gradient-based one-sided sampling (GOSS) algorithm to remove low-gradient data, minimizing time and memory consumption; (4) Using mutual exclusive feature (MEF) binding to reduce dimensionality by merging mutually exclusive features into one. The practical application of LightGBM has greatly expanded the scope and variety of GBDT’s applications, addressing the low training rate challenge of GBDT in large-scale data training [42,43,44].

3.3.4. AdaBoost

AdaBoost was proposed by Yoav Freund and Robert Schapire in 1995 [45]. Unlike GBDT, AdaBoost adjusts the data structure and the sample weights before building weak learners [46,47]. In classification problems, AdaBoost’s fundamental concept has two main aspects. Firstly, if a sample is misclassified by a previous weak learner, its weight will be increased in the next round of training, and if it is correctly classified, its weight will be reduced. An increased weight means that the misclassified samples in the previous round received more attention in subsequent weak learners’ training. Secondly, AdaBoost combines weak learners using a weighted majority vote mechanism. In each round, the classification error rate is calculated, and the weight of weak learners is adjusted relative to the error rate, with lower error rates leading to higher weights and vice versa. Ultimately, the weak learners with lower error rates significantly influence the voting decision, resulting in improved classification accuracy [48,49].

3.4. Rapid Landslide Extraction Models

3.4.1. Feature Selection Using SHAP

SHAP is a unified framework proposed by Lundberg in 2017 for explaining machine learning models [19]. The framework considers each feature to be a ‘contributor’ to the model and calculates its average marginal contribution to measuring its importance for the model’s classification result [50]. A high contribution indicates a feature’s high importance ranking. By utilizing this framework, high-dimensional features can be normalized and their dimensionality reduced without compromising the model’s accuracy [51]. When conducting a contribution analysis for landslide extraction, identifying the features that have the greatest impact can be challenging. Therefore, all 48 features listed in Appendix B, Table A2 were input into the basic model for analysis. Then, the features with a high contribution rank in landslide extraction were selected and retained in the final data column to produce optimal results.

3.4.2. Optuna Hyperparameter-Tuning

Optuna is a hyperparameter tuning tool that uses Bayesian optimization to determine the most promising hyperparameter selection and iteratively adjusts the search using a Bayesian optimization algorithm known as the Tree-structured Parzen Estimator [52]. It compensates for deficiencies in tuning approaches such as grid search and random search, which are computationally inefficient, costly, and susceptible to local optimum values. Furthermore, by optimizing the hyperparameters, it can adjust the model’s complexity to improve its generalization ability and learning effectiveness on independent datasets, enhancing its overall system performance [53].

Optuna has two main advantages: it allows users to dynamically construct the search space by expressing hyperparameter optimization as the process of a group of hyperparameters as input and solving for the minimum/maximum value of the objective function, and it uses efficient search and pruning strategies to improve its optimization efficiency [54]. Furthermore, the Optuna code is simple and easy to understand, and the optimization process only requires three core concepts: objective, trial, and study. The process for optimizing hyperparameters with Optuna involves four main steps: (1) Initialize the objective and define the type and search space of the parameters; (2) Initialize the study and set the optimization method and the maximum number of iterations for a single trial; (3) Begin tuning the hyperparameters; (4) Output the optimal hyperparameters and objective function value when the maximum number of iterations is reached [55].

3.4.3. Construction SHAP-OPT-XGBoost Landslide Extraction Model

To extract the landslide information, we developed SHAP-OPT-Models, as is shown in Figure 5, which involved various steps such as data preprocessing, selection of training, validation, testing sets, feature selection, automatic hyperparameter optimization, and model accuracy comparison. Initially, the multidimensional input features were processed during data preprocessing. Next, SHAP was used to explain the initial model established using the training set, and several high-contribution features were retained to generate an optimal feature data column. This column was then used to select the best feature data column from the training set to form the validation set, and the Optuna framework based on Bayesian optimization was used to automatically perform hyperparameter tuning. During the optimization iterations, the accuracy metric values obtained from 10-fold cross-validation on the validation set were used as references to evaluate the automatically generated parameter combinations and select the best hyperparameter configuration. Then OPT-SHAP-Models were evaluated for their performance on the testing set. To obtain the best strategy for landslide information extraction, we trained and compared SHAP-OPT-Models with a consistent set of sample data, initial feature data, feature selection criteria, and hyperparameter optimization. Finally, SHAP-OPT-XGBoost, which had the best performance, was selected to extract landslides. This approach allowed us to determine the most effective method for extracting landslide information.

3.5. Accuracy Evaluation

This paper evaluated the accuracy of SHAP-OPT-Models based on their accuracy, precision, recall, Kappa coefficient, and receiver operating characteristic (ROC) curve, combined with the area under curve (AUC) [56,57,58,59]. The indicators are expressed as follows, where TP, FP, FN, and TN are shown in

Table 1. Confusion matrix.

Prediction Situation	Actual Situation
	Landslide	Non-Landslide
Landslide	True positive (TP)	False positive (FP)
Non-landslide	False negative (FN)	True negative TN

. Specifically, TP indicates the total number of correct positive extractions, while FP represents the number of negative examples wrongly identified as positive. FN is the number of positive examples that were missed during extraction, and TN indicates the number of negative examples that were not extracted [60].

(1) Accuracy refers to the proportion of correctly classified samples over the total number of samples in a classification model. It is commonly used for balanced datasets and is calculated using the following formula:

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

(3)

where TP represents the number of true positives, FP represents false positives, TN represents true negatives, and FN represents false negatives.

(2) Precision refers to the proportion of correctly classified positive samples over the total number of samples classified as positive in a classification model. It is commonly used for imbalanced datasets and is calculated using the following formula:

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

(4)

(3) Recall refers to the proportion of true positive samples over the total number of true positive samples in a classification model. It is calculated using the following formula:

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

(5)

(4) Kappa is a measure of the classification model accuracy that takes into account random chance and random guessing. It is considered meaningful only if the prediction results of the classifier exceed random guessing. Kappa is calculated using the following formula:

$$\mathit{Kappa}=\frac{P_0-P_e}{1-P_e}$$

(6)

where P₀ represents the probability of correctly classified samples and P_e represents the probability of random guessing.

(5) ROC is a graphical plot that illustrates the performance of a binary classifier system as its discrimination threshold is varied. It is used to evaluate the performance of a classification model and to determine the optimal threshold value for binary classification. The x-axis of the ROC curve represents the false positive rate (FPR), while the y-axis represents the true positive rate (TPR). The FPR and TPR can be calculated as follows:

$$\mathit{FPR}=\frac{\mathit{FP}}{\mathit{TN}+\mathit{FP}}$$

(7)

$$\mathit{TPR}=\frac{\mathit{TP}}{\mathit{TP}+\mathit{FN}}$$

(8)

To obtain a ROC curve, the prediction results of a classifier are first sorted in descending order of their probabilities. The classification threshold is then gradually adjusted, and the corresponding FPR and TPR points are calculated, resulting in a complete ROC curve. The classifier performance can then be measured by computing the AUC value.

4. Results

4.1. Preprocessing of High-Resolution Images

Radiometric, atmospheric, and geometric corrections were carried out on the remote sensing imagery of Fengjie based on ENVI 5.6 software to reduce sensor noise, eliminate atmospheric scattering and absorption, and remove geometric distortions. Preprocessing also involved stitching together multiple images to create a seamless and complete picture of the area. Additionally, image fusion techniques were implemented to integrate complementary information from panchromatic and multispectral data. Part of the preprocessing results are shown in Figure 6 and Figure 7.

4.2. Segmentation

The utilization of high-resolution imagery datasets with a 2 m spatial resolution presents a significant challenge due to their massive size and computational hardware limitations. Consequently, the operational period of the landslide extraction model was extended, making it difficult to fulfill the demands of model training and iterative calculations. To tackle this issue, this paper used a regular grid method for partitioning the dataset into tiled data blocks.

In the multiscale segmentation process of ESP2, ROC-LV was calculated to measure the local variations in image object homogeneity at different segmentation scale parameters. When the rate of change exhibited a peak value, the corresponding scale indicated the optimal segmentation scale for the objects, as shown in Figure 8.

Several prominent peaks can be observed in the ROC-LV image, and the top four peak values, sorted in descending order, correspond to the scales illustrated in Figure 9b. To further confirm the optimal segmentation scale for matching landslide objects, a compactness parameter of 0.6 and a shape parameter of 0.2 were set. The image was then segmented using the segmentation scales corresponding to the top four peak values, resulting in the segmented images depicted in Figure 9c ①–④. Referring to the pre-segmentation image, as shown in Figure 9a, significant over-spitting of landslides occured at scale parameters of 118 and 160. At a scale parameter of 281, the landslides were completely segmented with clear contours, but the surrounding vegetation and cultivated land were overly segmented. At a scale parameter of 202, it effectively avoided under-segmentation of small-scale landslides and over-spitting of large-scale landslides. The integrity and clarity of the landslide objects were high, and the segmentation results for other land cover classes were also satisfactory. Therefore, the final determined image segmentation parameters were as follows: scale = 202, compactness = 0.6, shape factor = 0.2.

4.3. SHAP Feature Selection

As shown in Figure 10, SHAP was used to explain the features of GBDT, XGBoost, LightGBM, and AdaBoost. The features were ranked in descending order of importance based on the SHAP value’s absolute average, which reflects the feature’s impact on the model output. Positive and negative SHAP values indicate positive and negative effects, respectively. The color band in the figure shows that the redder the band, the larger the feature value, while the bluer the band, the smaller the feature value. Additionally, wider color bands correspond to more significant feature influence, indicating that such features are more critical.

For instance, for XGBoost, the mean slope of the landslide objects has the most significant impact. The SHAP values of the mean slope demonstrate a wide distribution with two clusters; the blue cluster indicating that a lower mean slope corresponds to a lower probability of landslide objects, while the red cluster indicates a higher probability of landslide objects within a certain range of mean slope. However, there is a decreasing trend in the probability of landslide objects outside of this range. Similarly, NDSI exhibits a feature performance similar to the mean slope. In contrast, the distribution pattern of the NDVI feature is less clear, with red scatter points and a small number of blue scatter points in the range of SHAP values less than 0. This indicates that high values of the NDVI feature hurt landslide detection, and the lower the value, the less significant probability of identifying the object as a landslide. This may be due to the low separability of landslide from non-landslide types, such as buildings and water bodies, in the NDVI feature. In GBDT, LightGBM, and AdaBoost, features like the mean slope, NDSI, and mean relief are important features for landslide extraction, as they exhibit a similar “blue on the left and red on the right” characteristic with high a contribution across these models. Specifically, the mean slope has a higher contribution ranking in BGDT, LightGBM, and AdaBoost, with the top positions in LightGBM and AdaBoost. NDSI is ranked first, third, and sixth in BGDT, LightGBM, and AdaBoost, respectively. Mean relief ranks 3rd in all three models. On the other hand, NDVI has a significant difference, and its contribution ranking in these three models is not as high as it is in XGBoost, and even does not rank in the top 15 in AdaBoost. Therefore, while reducing dimensionality, it is still necessary to construct a set of optimal multi-dimensional feature combinations for landslide extraction.

In order to optimize the feature selection process, all the features of the four models were ranked according to their importance, and the number of feature combinations was chosen with a step size of 2. The sample extraction accuracy was evaluated using AUC. As illustrated in Figure 11, XGBoost exhibited a distinct inflection point on the curve at eight features, achieving an extraction accuracy of 0.932. GBDT displayed an inflection point on the curve at 14 features, with a sample extraction accuracy of 0.820. AdaBoost achieved a sample extraction accuracy of 0.783, highlighting a clear inflection point on the curve with the use of 10 features. LightGBM utilizes six features, with a clear inflection point on the curve and a sample extraction accuracy of 0.902. It can be observed that the sample extraction accuracy did not continue to rise with an increase in the number of features, but gradually shifted from a balanced state to a slow decline. Excessive features not only increase the computational complexity but also generate information interference within the feature space, leading to a decrease in extraction accuracy. Therefore, based on a comprehensive analysis of the experimental results, GBDT selected NDSI, GLCM entropy (all dir.), mean relief, mean slope, and mean blue, along with the first 14 features, as the preferred feature set. XGBoost selected the mean slope, NDSI, NDVI, GLDV contrast (all dir.), GLCM entropy (all dir.), standard deviation green, mean relief, and standard deviation NIR as the preferred feature set. LightGBM selected the mean slope, GLCM Entropy (all dir.), mean relief, NDSI, mean red, and NDVI as the preferred feature set. AdaBoost preferred the mean slope, length/width, mean relief, GLCM mean (all dir.), and GLCM entropy (all dir.), along with the first 10 features.

4.4. Optuna Hyperparameter Tuning

When using Optuna for hyperparameter tuning, AUC was employed as the performance metric for the validation set. Specifically, the optimizer returned the average AUC value gathered using a 10-fold cross-validation mechanism as the objective value to the main function. Iterating the optimization process multiple times resulted in several optimization objective values. The parameter combination corresponding to the maximum value is the optimal hyperparameter scheme of the model.

By observing the interaction between the hyperparameters and the objective function during the automatic optimization process, as shown in Figure 12, it was found that the AUC values improved significantly, with the GBDT AUC increased from 0.8417 to 0.9167. Similarly, XGBoost demonstrated an increase in AUC from 0.8409 to 0.9335. LightGBM also showed an impressive performance gain, with its AUC rising from 0.8401 to 0.8964. The AUC for AdaBoost increased from 0.8131 to 0.9343. The start and end values of each column in the figure represent the search range of the corresponding hyperparameter. For example, the search range of the iteration number of parameters for XGBoost mainly concentrated on [300, 700]. The search range of the learning rate varied extensively, but was mainly concentrated below 0.3. The intensity of the blue line represents the average evaluation metric value of different parameter combinations during the optimization process, with darker blue lines indicating a better model performance.

As shown in Figure 13, the experiment was designed to execute the automated parameter tuning command 50 times. The orange line in the figure depicts the maximum achieved accuracy value during the hyperparameter optimization process, while the blue dots represent the current experiment’s AUC value. For GBDT, the cross-validated AUC values were mostly above 0.885, with the highest value of 0.9163 achieved in the 11th round of hyperparameter tuning. The cross-validated AUC values for XGBoost were predominantly above 0.9, with the highest value of 0.9324 obtained in the 21st round of hyperparameter tuning. Similarly, the cross-validated AUC values for LightGBM were primarily above 0.84, with the highest value of 0.8961 obtained in the 21st round of hyperparameter tuning. Finally, for AdaBoost, the cross-validated AUC values were generally above 0.91, with the highest value of 0.9342 obtained in the 35th round of hyperparameter tuning.

We retrained the model using the optimal hyperparameter combination, which was obtained when XGBoost achieved the top AUC score, as shown in

Table 2. XGBoost hyperparameter tuning results.

Hyperparameter	Defaults	Optimal Value	Parameter Meaning
learning_rate	0.3	0.25	learning rate
max_depth	6	10	the maximum depth of the tree
n_estimators	500	700	Number of estimators
min_child_weight	1	2	Min leaf weight
subsample	1	0.4	Subsample of training instances
colsample_bytree	1	0.5	Feature subsampling
reg_alpha	0	7	L1 regularization of weights
reg_lambda	1	4	L2 regularization of weights
gamma	0	0.3	Minimum loss reduction

.

Comparing the training process curves before and after hyperparameter optimization, as shown in Figure 14a, the performance of GBDT before optimization demonstrates subpar accuracy curves for both the training and test sets. Training with default iteration parameters, the accuracy curve for the training set failed to reach a stable state, resulting in a test set accuracy of 0.818. As shown in Figure 15, the experimental results reveal the performance of XGBoost before optimization to be excellent on the training set. The accuracy curve of the training set exhibited a turning point after the 19th iteration, and the AUC value gradually approached 1. However, its performance on the test set was poor. The model briefly improved its AUC value at the 11th iteration, but then gradually decreased until ultimately stabilizing at only 0.934. As shown in Figure 16, LightGBM demonstrated outstanding performance on the training set before optimization. The accuracy curve of the training set shows a turning point after the 35th iteration, and the AUC value gradually approached 1. However, it performed poorly on the test set, with a brief increase in the AUC value at the fifth iteration, followed by a gradual decrease, ultimately stabilizing at only 0.90. AdaBoost, before optimization, as shown in Figure 17a, exhibited excellent performance on the training set. The accuracy curve of the training set experienced a turning point after the 17th iteration, and the AUC value gradually approached 1. However, it performed poorly on the test set, with the accuracy curve fluctuating throughout, resulting in a test set accuracy of 0.83. This indicates that the models have weak generalization capabilities and lower system performance before optimization when applied to independent datasets.

In contrast, the model training curve based on the optimal parameter combination shows that the fitting effect of the training set and the test set was good. Compared to the pre-optimized models, GBDT improved its test set AUC value to 0.967. XGBoost showed a significant improvement in the AUC value on the test set, steadily stabilizing at 0.978. LightGBM improved to 0.974, and AdaBoost improved to 0.957.

4.5. Comparison of Model Accuracy

We compared the performance of SHAP-OPT-Models, which were trained using consistent initial feature datasets, feature selection rules, and hyperparameter optimization strategies. The optimal features were selected using the Optuna hyperparameter-tuning steps, resulting in fourteen features for GBDT, eight features for XGBoost, six features for LightGBM, and ten features for AdaBoost. The accuracy of the four SHAP-OPT-Models was evaluated on the test set, as shown in

Table 3. Accuracy comparison of SHAP-OPT-Models.

Algorithm	Accuracy/%	Precision/%	Recall/%	Kappa	Training Time/s
SHAP-OPT-XGBoost	96.26	90.91	85.71	0.8602	1.16
SHAP-OPT-GBDT	93.93	82.35	80.00	0.7754	1.28
SHAP-OPT-LightGBM	95.79	86.11	88.57	0.8480	1.06
SHAP-OPT-AdaBoost	92.99	83.33	71.43	0.7282	0.97

and Figure 18. All the SHAP-OPT-Models achieved high accuracy on the test set, demonstrating good performance. Among them, SHAP-OPT-XGBoost outperformed the others, with an accuracy of 96.26%, precision of 90.91%, Kappa coefficient of 0.8602, training time of 1.16 s, and AUC value of 0.9705. However, its recall rate was slightly lower than that of SHAP-OPT-LightGBM, with 85.71% compared to SHAP-OPT-LightGBM’s 88.57%. Moreover, all the SHAP-OPT-Models demonstrated rapid and efficient performance, with short landslide extraction times. Overall, SHAP-OPT-XGBoost showed superior performance, achieving a higher accuracy in landslide extraction while maintaining the extraction speed.

4.6. Landslides Information Extraction

In this paper, a SHAP-OPT-XGBoost rapid landslide extraction model was built that successfully extracted landslide spatial information in the whole Fengjie County from 2013 to 2020. SHAP-OPT-XGBoost’s accuracy was verified, and the results indicated that it effectively and rapidly extracted the landslide information. Our findings provide a valuable tool for disaster management and mitigation efforts in the region, helping to reduce the impact of landslides on local communities and infrastructure.

4.6.1. Typical Individual Landslide Analysis

Figure 19 shows the temporal and spatial changes in a single landslide that were extracted using the SHAP-OPT-XGBoost model. The topography of the landslide area is generally north-high and south-low. The landslide crest sits at a high point of approximately 1385 m, while the foot of the landslide is a small tributary of the Meixi River at an elevation of approximately 1000 m. The slope gradient is approximately 25° to 45° with a height difference of 385 m. Detailed field investigations and landslide data records indicate that the Miaowan landslide, located in Changping Village, Ping’an Township, Fengjie County, which occurred in September 2002, is an old landslide. In 2013, the landslide’s boundaries were difficult to distinguish, and its shape was extremely irregular, with an area of approximately 4.4 × 10⁴ m². The soil in the landslide body was dense, lacking discernible slide marks, and it had a large exposed bedrock area. The front slope of the landslide was steep and accompanied by high vegetation coverage. It can be inferred that the landslide underwent a series of transformations, such as erosion and leveling after its initial occurrence, and its basic elements had gradually weakened. However, traditional remote sensing interpretation methods may not accurately and quickly make decisions. In 2015, the area of the landslide was approximately 8.0 × 10⁴ m², and the difference between the landslide and the surrounding objects significantly increased. The landslide body was lighter in color, with a grayish-white hue, shaped like a tongue, with a clear boundary. The middle part of the slope was lower than the sides, with no significant slide markings, but with evident piles of pebbles and debris on the gentle slope below. In 2018, the extracted landslide almost completely covered the surface that was present in 2013. The area of the landslide expanded to 8.7 × 10⁴ m² compared to 2015, with the right wing extending and a small-scale landslide appearing. The vegetation coverage of the landslide remained low, and the surface markings were smoother and more even compared to 2015. In 2020, the area of the landslide was the smallest: approximately 1.7 × 10⁴ m². The surrounding vegetation had recovered well, mainly concentrated in the upper-middle part of the landslide, with the color of the landslide gradually darkening and the slide marks gradually weakening. The small landslide added in 2018 on the right side of the landslide was almost completely covered by vegetation.

According to the investigation and prevention proposal for the Miaowan landslide, the stability of the landslide is affected by both intrinsic and external factors. Intrinsic factors include the geological structure, topography, material composition, and vegetation of the landslide zone. External factors mainly consist of atmospheric precipitation, changes in the water level of the reservoir, and human engineering activities. Among them, the storage of water in the reservoir and rainfall are important factors influencing the stability of the landslide. Changes in the water level of the reservoir can cause the surface water to infiltrate into the landslide, resulting in soil softening and decreased mechanical strength, as well as an increase in hydrostatic pressure that is unfavorable for the stability of the landslide. Rainfall can saturate the landslide soil, increase the load, and reduce the mechanical strength, negatively affecting the stability of the landslide. Other factors such as topography, geological conditions, material composition, and human engineering activities can also interact and cause changes in the stability of the landslide.

An analysis of true-color remote sensing images in 2020 and multi-temporal landslide vector data revealed that the landslide had a period of strong activity from 2013 to 2018, followed by weakened activity from 2018 to 2020. Despite vegetation cover showing signs of recovery in 2020, the landslide is situated on a hidden fault zone where the developed fault represents a surface rupture or stress concentration zone, increasing the likelihood of the landslide recurring. As such, this area should remain a priority focus for landslide monitoring.

4.6.2. Regional Landslides Analysis

To quantitatively evaluate the accuracy of SHAP-OPT-XGBoost, several performance metrics were employed, including accuracy, misclassification rate, and omission rate. The evaluation process involved the utilization of historical landslide disaster data, along with three-dimensional interactive manual interpretation and field investigation photos to rectify the hidden danger points. The rectification included adding unrecorded landslide point attributes, supplementing and updating missing values of landslide occurrence time, excluding individual abnormal landslide hidden danger points, and dividing the reference data for each period of landslide accuracy evaluation. Subsequently, the accuracy verification reference data was used in ArcGIS to analyze the spatial correlation between landslide disaster points and extraction results using spatial connection technology. In the northwest of Fengjie, eight townships had relatively dense landslide distributions, which were selected as sub-areas. Taking the sub-area as an example, the calculated landslide extraction accuracy indicators for four periods are presented in

Table 4. Extraction accuracy of landslides in sub-areas in different years.

Year	Accuracy/%	Misclassification Rate/%	Omission Rate/%
2013	74.31	25.69	16.92
2015	86.76	13.24	8.34
2018	89.77	10.23	18.56
2020	82.14	17.86	11.33

.

A comparative analysis of the landslide extraction results is shown in Figure 20. A statistical analysis revealed that SHAP-OPT-XGBoost correctly extracted 52 landslide objects while erroneously extracting or omitting 11 landslide objects. Of the misclassified objects, five were falsely identified as non-landslides, and six were falsely identified as landslides. The accuracy of the landslide information extraction was 91.23%, and the misclassification rate and omission rate were 8.77% and 10%, respectively. Figure 20b visually demonstrates the clear and identifiable contour of the correctly extracted landslide, showcasing its prominent texture and distinct geometric shape. The omission of landslides can be attributed to natural erosion and human transformation that occurred over several years, resulting in less obvious characteristic elements. On the other hand, some landslides experienced significant vegetation recovery within their surfaces, gradually concealing the exposed features present during the initial stage of occurrence. This vegetation cover makes it challenging for SHAP-OPT-XGBoost to recognize these landslides. Misclassifying non-landslide surfaces as landslides is a common issue caused by the similarity in spectral and terrain information between cultivated or bare lands and actual landslide objects. However, upon verification, these surfaces lack the necessary environmental conditions and significant characteristic features for landslide occurrence.

Additionally, a three-dimensional scene view, generated by overlaying the correctly extracted landslide vector data, remote sensing image, and DEM, provides a more comprehensive representation, as shown in Figure 21. This view offers clear insights into the landslide’s occurrence location, range, shape, as well as the surrounding terrain and landform environment. Based on Figure 21, it is evident that the slope exhibits a high degree of landslide development, with severe destruction of the original vegetation cover and debris accumulation at the lower edge of the landslide body. Some of the debris has been carried into the river by rainwater, resulting in sedimentation and blockage.

Using SHAP-OPT-XGBoost, we extracted the landslide spatial distribution data of Fengjie County from 2013 to 2020. The landslide extraction results of eight townships in the northern part of Fengjie are shown in Figure 22.

In the eight townships, the total area of landslides was 8.44 km² in 2013, accounting for only 0.21% of Fengjie’s total area. The landslide scale was mainly small to medium-sized, and they were scattered and sparse in spatial distribution, with denser distributions near river junctions and valleys. Their overall development degree was low. By 2015, the total area of the landslides had increased to 20.89 km², accounting for 0.51% of Fengjie’s total area. This was a 148% increase compared to 2013, with its density showing a significant increase during 2013–2015. The landslide expansion was mainly distributed in a belt shape along various water systems, with the Kangle Township, Qinglian Township, and Dashu Township having the largest areas. After consulting the relevant books and data, it was concluded that the increase in landslides in 2015 was due to the impact of extremely heavy rains in 2014, which induced multiple landslides of different sizes and quantities. In 2018, the landslide area was decreased to 15.25 km², accounting for 0.37% of Fengjie’s total area; a 0.14% decrease compared to 2015. During 2015–2018, the landslides were mainly concentrated in the southeast and northwest of the sub-area, and the new landslides mainly appeared in the northwest, east and southeast of the sub-area. Overall, the number of newly formed landslides did not increase significantly, and the landslide area decreased compared to 2015. By 2020, the extracted landslide area was 17.29 km², accounting for 0.42% of Fengjie’s total area: an increase of 0.05% compared to 2018. During 2018–2020, the growth rate of the landslide area was relatively small, mainly concentrated on the banks of the mainstream and some tributaries. Additionally, during this period, the degree of development of bank landslides was high, and the phenomenon of resurrecting old landslides was significant.

In summary, the landslide area had significant interannual changes. The landslides remained at a high level of development throughout the year, and if extreme weather was encountered such as heavy rains, they would show explosive growth. The spatial distribution of landslides between townships was uneven, but they were distributed in a concentrated manner in certain areas. Due to the close relationship between landslides and water systems, they were distributed in a strip shape along valleys and both riverbanks.

4.6.3. Time-Series Landslides Extraction in Fengjie from 2013 to 2020

The landslide contour boundaries extracted by SHAP-OPT-XGBoost exhibited a high level of stability and reliability with clear edges. The landslide area statistics for each year from 2013 to 2020 are presented in

Table 5. Landslide area extracted based on SHAP-OPT-XGBoost from 2013 to 2020.

Year	2013	2014	2015	2016	2017	2018	2019	2020
landslide area (km2)	19.26	32.92	55.42	42.74	39.65	45.76	35.02	39.35

. To highlight the landslide information, the other land cover types are not represented in the spatial distribution map, as illustrated in Figure 23. The results show that landslides in Fengjie County mainly occur in the north and southwest regions, with the northern area exhibiting the most significant development, with a larger area and more concentrated distribution. Conversely, the eastern and southern regions exhibit lower levels of development. The areas of high landslide development overlap with high-risk geological disaster zones, indicating a close relationship between landslide hazard and development distribution. The spatial distribution of landslides in this region underwent significant changes from 2013 to 2020, and the average level of development fluctuated but has remained relatively high. Therefore, effective monitoring and prevention strategies are necessary, and a reliable landslide disaster warning system must be established. Further optimization is required to enhance the accuracy and precision of landslide detection and analysis in this region.

5. Discussion

(1) Through comparing the performance of the four basic models before and after optimization, it is demonstrated that the optimal parameter combinations established by Optuna effectively reduce the computational complexity and prevent overfitting, thus enhancing the models’ performance. However, the generalizability of the SHAP-OPT-Models to different datasets and environmental contexts remains unclear. Further studies will be needed to address this.

(2) When evaluating the accuracy of the SHAP-OPT-XGBoost rapid landslide extraction model, factors such as human activities, vegetation restoration, and similarities in spectral and terrain information between cultivated or bare land and landslides can lead to the misclassification of non-landslide surfaces as landslides. In future work, a multi-source data fusion approach that combines optical remote sensing data with other relevant data types, including radar remote sensing data, can be used to improve the accuracy of the SHAP-OPT-XGBoost landslide extraction model.

(3) While extracting landslide information, it has been observed that individual landslides generally exhibit limited spatial variability over the long term. The location is closely related to topographic factors such as slope gradient. Variation primarily manifests in terms of color, texture, geometric shape, and vegetation coverage of the landslide. However, the spatio temporal distribution pattern is complex and susceptible to external factors, such as sudden heavy rainfall, which can lead to a sudden increase in landslide magnitude. Furthermore, fluctuations in water levels in reservoir areas can also increase the likelihood of landslides. Therefore, it is crucial to pay attention to climate changes and establish monitoring and early warning systems. During the rainy season, areas prone to landslide hazards should implement warning mechanisms for prompt evacuation and relocation efforts.

6. Conclusions

In this paper, SHAP-OPT-XGBoost was proposed by introducing SHAP for feature selection and the Optuna framework for Bayesian hyperparameter tuning, which enables rapid extraction of spatial distribution information and quantitative analysis of landslide spatiotemporal distribution.

(1) SHAP-OPT-XGBoost achieved higher accuracy, precision, Kappa coefficient, and AUC values (96.26%, 90.91%, 0.8602, and 0.9705, respectively) than SHAP-OPT-GBDT, SHAP-OPT-LightGBM, and SHAP-OPT-Adaboost. Meanwhile, the SHAP-OPT-Models’ training time was less than 1.3 s using mediocre computational hardware, demonstrating a rapid and efficient performance. Although there were some cases of misclassification and omission, SHAP-OPT-XGBoost’s overall performance was significant and met the accuracy requirements for landslide spatial information extraction in the study area.

(2) From 2013 to 2020, the overall developmental degree of landslides in Fengjie was high, and the distribution density of landslides increased significantly from 2013 to 2015. The expansion trend of landslide spatial distribution was vigorous. In terms of the time dimension, the landslides in the study area maintained a high level of development throughout the year since 2013. The landslide area in 2015 increased by approximately 55.42 km²; a 188% increase from 2013, showing an explosive growth pattern. In terms of their spatial dimensions, the landslides were unevenly distributed between towns and exhibited a belt-shaped concentrated distribution along the valleys and on both sides of the river, with characteristics of regional concentration.

(3) The spatial and temporal distribution of landslides is the result of the comprehensive effects of various factors, including terrain environment, geological conditions, meteorology, hydrology, human engineering, and social economy. In the future, it will be necessary to combine multi-source and multi-temporal data to explore the essential relationships between the developmental laws of landslides and disaster-prone factors such as terrain environment, as well as to analyze the disaster-prone mode of landslide hazards and further improve disaster prevention and relief plans. To achieve this, we can further develop and optimize remote sensing methods and techniques to obtain more comprehensive and accurate information on landslides.