A landslide susceptibility prediction method combining terrain combination features and partition explainable analysis

Abstract

The application discloses a landslide susceptibility prediction method combining terrain combined features and partition explainable analysis, and relates to the technical field of geological disaster evaluation. Slope, plane curvature and profile curvature factors are extracted from a digital elevation model of a research area, and after standardization, Gaussian mixture model is used for unsupervised clustering. The class labels obtained by clustering are taken as terrain combined features and are coded, and are combined with landslide influence factors to form an input feature set. The input feature set is input into an XGBoost model for training and the hyperparameters are optimized through grid search. Based on the trained model, the SHAP method is used to calculate the global feature contribution degree, and the contribution degrees of various factors are calculated according to terrain combination types, and the spatial differences of the factors under different terrain units are analyzed. The terrain types obtained by clustering are used as explicit features to participate in model training, and the spatial heterogeneity of the factors is revealed through grouped contribution degree analysis, which provides a new technical means for landslide susceptibility evaluation.

Description

Technical Field

[0001] This invention belongs to the field of geological disaster prediction and evaluation technology, specifically involving a landslide susceptibility prediction method that integrates topographic combination features and zoning interpretable analysis. Background Technology

[0002] Landslide susceptibility assessment is a fundamental task in regional geological disaster prevention and control. Its task is to identify the spatial distribution patterns of potential landslides, providing a basis for land use planning and engineering site selection. With the development of remote sensing technology and geographic information systems, machine learning methods have been widely used in landslide susceptibility prediction.

[0003] In existing methods for predicting landslide susceptibility, there is room for improvement in the way topographic factors are used. Slope, planar curvature, and profile curvature are three basic factors describing topographic features. Slope reflects the degree of surface inclination, planar curvature reflects the morphological changes of the slope along contour lines, and profile curvature reflects the morphological changes of the slope along the direction of maximum slope. From the perspective of geomorphic evolution and slope stability, there is a certain physical relationship between these three factors: in areas with steep slopes, the slope curvature often determines the stress distribution pattern and hydrological response characteristics within the slope. On convex slopes, under gravity, a tensile stress zone is generated at the top and a shear stress concentration zone at the toe; this stress distribution is conducive to slope instability. Concave slopes, on the other hand, have water-collecting characteristics, with surface runoff and groundwater flowing towards the middle and lower parts of the slope, leading to increased water content and pore water pressure in the soil and rock, and a corresponding decrease in effective stress. At the same time, planar curvature affects the degree of convergence or divergence of water flow; in areas with negative planar curvature, water flow converges, further amplifying the impact of hydrological processes on slope stability. Therefore, specific combinations of slope and curvature correspond to significantly different slope mechanical states and hydrogeological conditions, which together control the spatial development pattern of landslides.

[0004] However, existing technologies typically input the three topographic factors mentioned above as independent variables into the model, leading to the failure to distinguish their synergistic relationships. Some studies have attempted to classify topographic factors by setting thresholds, but the selection of thresholds lacks objective and unified standards, and the differences in topographic characteristics across different study areas make it difficult to generalize fixed thresholds. Regarding clustering analysis methods, existing studies have used K-means clustering or Gaussian mixture models for geomorphic zoning, but the clustering results are usually only used for sample stratification or zonal modeling. The topographic type information obtained from clustering is not used in model training as explicit features, and its multi-factor combination information is not fully utilized. In terms of model interpretability analysis, the SHAP method has been used to calculate the global importance ranking of various influencing factors. However, landslides, as complex geological phenomena influenced by multiple factors, have dominant factors that vary with topographic conditions. Global ranking reflects the average situation across the entire study area and cannot reflect the differences in the effects of factors within different topographic units. For example, on steep and convex ridges, slope geometry plays a more prominent role in controlling stability; while on steep and concave valley slopes, hydrological conditions reflected by water system distance may become the dominant factor. A global analysis alone would mask this spatial heterogeneity, making it difficult to provide targeted decision-making basis for regional prevention and control.

[0005] Based on the above analysis, there is a need in the existing technology for a landslide susceptibility prediction method that can objectively integrate multi-factor topographic information, explicitly incorporate it into the prediction model training, and reveal the spatial differentiation of factor contributions under different topographic conditions. Summary of the Invention

[0006] The purpose of this invention is to provide a landslide susceptibility prediction method that integrates topographic combination features and zonal interpretable analysis. This method integrates slope, plane curvature and profile curvature into topographic combination features through a Gaussian mixture model. After encoding, these features are used as input variables in the training of the XGBoost model. Furthermore, SHAP analysis grouped by topographic type reveals the spatial differences in the contribution of factors.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] A landslide susceptibility prediction method integrating topographic features and zonal interpretable analysis includes the following steps:

[0009] S1, Data Acquisition

[0010] Acquire digital elevation model data, landslide logging data, non-landslide sample data, and landslide influencing factor data for the study area. The non-landslide sample data is selected from non-landslide areas within the study area.

[0011] The landslide influencing factor data are selected from at least one of geomorphological factors, geological structure factors, hydrological environmental factors, vegetation cover factors, and human engineering activity factors.

[0012] Specifically:

[0013] Digital elevation model (DEM) data of the study area is acquired, which is used to extract basic topographic factors. Historical landslide logging data of the study area is acquired, including the spatial coordinates and boundary ranges of landslide locations. Non-landslide sample data is selected from non-landslide areas within the study area, with the number of selected samples balanced with the number of landslide samples. Landslide influencing factor data is acquired, which is selected from at least one of geomorphological factors, geological structure factors, hydrological environmental factors, vegetation cover factors, and human engineering activity factors.

[0014] S2. Extraction and Standardization of Basic Topographic Factors

[0015] Based on the digital elevation model data, slope factor, plane curvature factor, and profile curvature factor are extracted, and the extracted factors are standardized; specifically:

[0016] Three basic terrain factors are extracted based on the digital elevation model data.

[0017] The slope factor reflects the degree of inclination of a surface unit relative to the horizontal plane. It is calculated by the local window fitting method. A moving window is constructed with the target grid as the center, and the local surface equation is fitted. The maximum slope value is obtained by the first-order partial derivative, and the unit is degrees.

[0018] The plane curvature factor reflects the change in slope morphology along the contour line. A positive value indicates that the slope has a divergent shape, while a negative value indicates that the slope has a convergent shape.

[0019] The profile curvature factor reflects the change in slope morphology along the direction of maximum slope. A positive value indicates that the slope is convex, while a negative value indicates that the slope is concave.

[0020] The extracted slope factor, plane curvature factor, and profile curvature factor are standardized to eliminate the influence of dimensional differences on subsequent clustering. Standardization can be achieved using the min-max standardization method, mapping the values ​​of each factor to a uniform interval.

[0021] S3, Construction of Terrain Combination Features

[0022] Unsupervised clustering of the standardized slope factor, plane curvature factor, and profile curvature factor was performed using a Gaussian mixture model to divide the study area into multiple terrain combination types. Each grid cell was assigned a terrain combination category label as a terrain combination feature. The number of clusters was adaptively determined by the information criterion.

[0023] The number of clusters is determined by the Bayesian information criterion, and the number of clusters with the smallest Bayesian information criterion value is selected as the number of the multiple terrain combination types.

[0024] The multiple terrain combination types are automatically clustered by the Gaussian mixture model based on the statistical distribution characteristics of slope factor, plane curvature factor and profile curvature factor. Each type has significant differences in slope value range and positive and negative curvature combination form, and corresponds to different topographic units in the study area.

[0025] The clustering process of the Gaussian mixture model includes: assuming that the terrain factor data is generated by mixing multiple Gaussian distributions, estimating the mean vector, covariance matrix and mixing weight of each Gaussian component through the expectation-maximization algorithm, calculating the posterior probability of each grid cell belonging to each cluster, and taking the cluster corresponding to the maximum posterior probability as the terrain combination category of the grid cell.

[0026] Specifically:

[0027] Unsupervised clustering of standardized slope, planar curvature, and profile curvature factors was performed using a Gaussian mixture model. The Gaussian mixture model assumes that the observed data are generated by a mixture of multiple Gaussian distributions, and its probability density function is expressed as a weighted sum of the probability densities of each Gaussian component. Model parameters include the mean vector, covariance matrix, and mixture weights of each Gaussian component. These parameters are iteratively estimated using an expectation-maximization algorithm: in the expectation step, the posterior probability of each sample belonging to each Gaussian component is calculated based on the current parameters; in the maximization step, the parameter estimates of each component are updated based on the posterior probabilities. This iteration continues until the parameters converge. Finally, the posterior probability of each grid cell belonging to each cluster is calculated, and the cluster with the highest posterior probability is taken as the terrain combination category of that grid cell.

[0028] The number of clusters was determined using the Bayesian information criterion. The Bayesian information criterion value strikes a balance between model goodness of fit and model complexity; a smaller value indicates better overall model performance on given data. Bayesian information criterion values ​​were calculated for different preset numbers of clusters, and the number that minimized this value was selected as the final number of clusters. In practical calculations, for situations with a large amount of raster data across the entire area, spatial sampling can be used to calculate the Bayesian information criterion value to improve computational efficiency. The sampling method should ensure spatial representativeness and statistical unbiasedness of the samples. The adaptive determination of the number of clusters ensures the applicability of the method under different terrain conditions.

[0029] After clustering, each raster cell is assigned a corresponding terrain combination category label, which constitutes the terrain combination feature. Each cluster type is automatically formed by Gaussian mixture model based on the statistical distribution characteristics of slope, planar curvature, and profile curvature, taking into account the spatial heterogeneity of the data itself. The types have significant differences in the range of slope values ​​and the positive and negative combination forms of curvature, and correspond to different topographic units in the study area.

[0030] S4. Input Feature Set Construction

[0031] The terrain features are encoded, and the encoded features are combined with landslide impact factor data to construct an input feature set. The input feature set is then subjected to collinearity factor removal.

[0032] The encoding process employs at least one of one-hot encoding, label encoding, or target encoding to convert terrain combination category labels into numerical features, enabling categorical features and continuous influencing factors to participate in model training within a unified feature space.

[0033] The collinearity factor removal process employs variance inflation factor analysis or Pearson correlation coefficient analysis to remove factors exhibiting strong collinearity.

[0034] Specifically:

[0035] Terrain combination features are categorical variables; their label values ​​only indicate category distinctions and do not have numerical magnitude. To enable this feature to be trained in the same feature space as the continuous landslide influencing factor, the terrain combination feature is encoded. Encoding can employ at least one of one-hot encoding, label encoding, or target encoding to convert the terrain combination category labels into numerical features. After encoding, the categorical features are transformed into numerical features with dimensions corresponding to the original number of categories, which, together with the continuous influencing factor, construct the input feature matrix.

[0036] The encoded terrain features are fused with various landslide influencing factors to construct an initial input feature set. Collinearity testing and elimination are performed on the initial feature set: Pearson correlation coefficients are calculated between factors; for factor pairs with absolute correlation coefficients exceeding a preset threshold, one factor is retained based on its geological correlation with landslide occurrence—when the correlation coefficients of two factors exceed the preset threshold, factors with higher data acquisition accuracy, clearer physical mechanisms, or a more direct causal relationship with landslide development are prioritized for retention; the variance inflation factor value of each factor is calculated, and factors with variance inflation factor values ​​exceeding a preset threshold are eliminated. After collinearity elimination, the retained factors constitute the final input feature set used for model training.

[0037] S5 and XGBoost model training

[0038] The input feature set is input into the XGBoost model for training. The hyperparameters of the XGBoost model are optimized by a grid search method to obtain a trained landslide susceptibility prediction model. The hyperparameter optimization of the XGBoost model includes grid search optimization of at least one of the learning rate, maximum tree depth, subsample ratio, and regularization parameter.

[0039] Specifically:

[0040] The XGBoost model was used to train the input feature set. XGBoost is an ensemble learning algorithm based on the gradient boosting framework, which iteratively constructs multiple weak learners and optimizes the objective function using gradient descent. This algorithm features regularization support, automatic handling of missing values, and the ability to capture high-order interactions between features. Since terrain features essentially reflect the synergistic relationship between slope, planar curvature, and profile curvature, their informational value requires the model to capture nonlinear interactions between features to be realized; the XGBoost model effectively meets this requirement.

[0041] The training process includes a hyperparameter optimization phase. A grid search method is used to traverse various parameter combinations within a preset range of values. Cross-validation is employed to evaluate the model performance under each parameter group, with the area under the receiver operating characteristic (ROC) curve as the primary evaluation metric to select the optimal hyperparameter combination. Hyperparameters requiring optimization may include one or more of the following: learning rate, maximum tree depth, subsample ratio, L1 regularization parameter, and L2 regularization parameter.

[0042] After the model is trained, it is applied to the grid cells of the entire study area to output the predicted probability of landslide occurrence for each grid cell.

[0043] S6. Grouped SHAP Interpretability Analysis

[0044] Based on the landslide susceptibility prediction model that has been trained, the SHAP method is used to calculate the global feature contribution. The grid cells of the study area are grouped according to the terrain combination type, and the SHAP value of each landslide influencing factor in each group is calculated. The spatial differentiation pattern of factor contribution under different terrain combination types is analyzed.

[0045] The SHAP analysis grouped by terrain combination type includes: dividing all raster cells in the study area into multiple subsets according to terrain combination category labels, calculating the mean SHAP contribution of each landslide influencing factor for each subset, and determining the dominant influencing factor under each type by comparing the contribution differences of the same factor among different subsets; wherein the SHAP value of each subset is calculated based on the same global model trained in step S5, so that the comparison of factor contributions between different terrain types is carried out at the same prediction logical scale.

[0046] Specifically:

[0047] Based on the trained XGBoost model, the SHAP method is used for interpretability analysis. The SHAP method originates from the Shapley value theory in game theory, and its core idea is to decompose the model's prediction for a single sample into the sum of the contributions of each feature. For tree models, the TreeSHAP algorithm can efficiently calculate the marginal contribution of each feature to the prediction result for each sample.

[0048] It is important to clarify that the SHAP values ​​of each subset in the group analysis are calculated based on the same globally trained model, rather than training independent models for each subset separately. This approach ensures that the comparison of factor contributions between different terrain units is conducted on the same predictive logistic scale, directly reflecting the sensitivity differences of the same model under different terrain conditions, thereby revealing the spatial heterogeneity of factor contributions.

[0049] SHAP analysis is divided into two levels.

[0050] The first level is the global feature contribution analysis. The mean of the absolute values ​​of each feature's SAP in all test samples is calculated, and the samples are sorted in descending order to obtain the global feature contribution ranking, reflecting the average importance of each factor across the entire study area.

[0051] The second level involves group feature contribution analysis. All raster cells in the study area are divided into multiple subsets based on the terrain combination category labels obtained in step S3, with each subset corresponding to a specific terrain combination type. The mean SHAP value of each landslide influencing factor is calculated for each subset. By comparing the contribution differences of the same factor in different terrain type subsets, the spatial variation of the factor's intensity under different terrain conditions can be quantitatively revealed; by comparing the contribution ranking of factors within the same terrain type subset, the dominant influencing factor under that terrain condition can be identified.

[0052] S7 Landslide Susceptibility Mapping

[0053] Based on the output of the landslide susceptibility prediction model, a spatial distribution map of landslide susceptibility in the study area is generated. Specifically, the predicted landslide occurrence probability value of each grid cell output by the model in step S5 is spatially visualized to generate a spatial distribution map of landslide susceptibility in the study area. The probability values ​​can be divided into different susceptibility levels using the natural breakpoint method or the geometric interval method. Combined with the grouped SHAP analysis results in step S6, the dominant influencing factor types and their corresponding terrain combination characteristics within each susceptibility level area can be identified.

[0054] The present invention has the following beneficial effects:

[0055] By organically integrating slope, planar curvature, and profile curvature into a combined terrain feature using a Gaussian mixture model, the method achieves objective condensation of terrain information in a data-driven manner, effectively overcoming the shortcomings of traditional methods that use terrain factors independently or subjective thresholding. The number of clusters is adaptively determined using the Bayesian information criterion, ensuring the applicability of the method under different terrain conditions.

[0056] By encoding the combined terrain features and incorporating them as explicit input variables into the XGBoost model training, the model is forced to focus on the comprehensive geomorphological information reflected by the combination of multiple terrain factors, which is different from the indirect use of clustering results for partition modeling or sample selection in existing technologies.

[0057] By employing the SHAP analysis method, which groups factors according to terrain combination type, and calculating the contribution of each factor within the same globally trained model, the spatial differentiation patterns of landslide influencing factors are revealed. This grouped interpretation method allows the local dominant factors that were previously masked by the global average effect to emerge, providing a quantitative analytical tool for understanding the differences in controlling factors of landslide occurrence in different geomorphic locations.

[0058] Furthermore, since the XGBoost model can effectively capture the interaction relationships between features, the multi-factor synergistic information contained in the terrain combination features is fully utilized in this model, providing a practical basis for model selection in landslide susceptibility prediction. Attached Figure Description

[0059] Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention.

[0060] Figure 2 This is a schematic diagram of the spatial distribution of terrain combination types in an embodiment of the present invention.

[0061] Figure 3 This is a comparison chart of the subject operating characteristic curves in two scenarios: one with terrain combination features and one without.

[0062] Figure 4 This is a ranking diagram of global SHAP feature contributions in an embodiment of the present invention.

[0063] Figure 5 This is one of the comparison charts of SHAP contribution calculated by grouping according to terrain combination type in an embodiment of the present invention.

[0064] Figure 6 This is the second comparison chart of SHAP contribution calculated by grouping according to terrain combination type in this embodiment of the invention.

[0065] Figure 7 This is the third comparison chart of SHAP contribution calculated by grouping according to terrain combination type in this embodiment of the invention.

[0066] Figure 8 This is the fourth comparison chart of SHAP contribution calculated by grouping according to terrain combination type in the embodiments of the present invention.

[0067] Figure 9 This is a spatial distribution map of landslide susceptibility generated in an embodiment of the present invention. Detailed Implementation

[0068] The present invention will now be described in detail with reference to embodiments and accompanying drawings. However, it should be understood that the embodiments and drawings are for illustrative purposes only and do not constitute any limitation on the scope of protection of the present invention. All reasonable modifications and combinations included within the inventive spirit of the present invention fall within the scope of protection of the present invention.

[0069] This embodiment uses a specific location as the study area. This area is located in the southern segment of the Longmenshan fault zone, characterized by significant topographic relief, complex geological structures, and well-developed landslide hazards. It possesses the geological and geomorphological conditions necessary for verifying the method of this invention. According to... Figure 1 The process shown includes the following steps:

[0070] Step S1: Data Acquisition

[0071] Digital elevation model data for the region was acquired, with a spatial resolution of 12.5m. This data was derived from ALOS PALSAR satellite observation data products.

[0072] Historical landslide data for the study area were obtained from local geological hazard survey data, including the spatial coordinates and boundary ranges of landslide locations. In this embodiment, a total of 203 landslide sample points were collected. Outside the buffer zone surrounding the landslide points, 203 non-landslide sample points were randomly selected, equal in number to the number of landslide sample points, to form a balanced training sample set.

[0073] Landslide influencing factors data were obtained, specifically including eight factors: elevation, aspect, slope, distance to water system, distance to fault, distance to road, vegetation cover NDVI, and average annual rainfall. Elevation was directly obtained from the Digital Elevation Model (DEM). Slope aspect and slope were calculated using ArcGIS software from the DEM. Distance to water system was calculated using Euclidean distance from regional water system linear vector data. Distance to fault was calculated using Euclidean distance from regional fault structure linear vector data. Distance to road was calculated using Euclidean distance from regional road linear vector data. Vegetation cover NDVI was calculated from the near-infrared and red band reflectance of Landsat remote sensing imagery. Average annual rainfall was obtained from regional meteorological station observation data through spatial interpolation. All factor data were uniformly resampled to a spatial resolution of 12.5m, consistent with the DEM.

[0074] Step S2: Extraction and standardization of basic terrain factors

[0075] Based on digital elevation model data, three basic terrain factors are extracted.

[0076] When extracting the slope factor, a 3×3 moving window is constructed with the target grid as the center. The elevation values ​​of each grid within the window are fitted with a quadratic surface equation. The slope value of the central grid is calculated from the first-order partial derivative of the surface equation, in degrees.

[0077] The plane curvature factor is obtained by calculating the rate of change of slope aspect along the contour lines. Positive values ​​indicate that the slope surface diverges along the contour lines, while negative values ​​indicate a convergent shape.

[0078] The profile curvature factor is obtained by calculating the rate of change of the slope along the direction of maximum slope. A positive value indicates that the slope is convex upwards along the direction of maximum slope, while a negative value indicates that it is concave downwards.

[0079] The three factors were subjected to min-max standardization to map their values ​​to the [0,1] interval. All of these operations were performed in ArcGIS software.

[0080] Step S3: Construction of Terrain Combination Features

[0081] Unsupervised clustering of standardized slope, plane curvature, and profile curvature factors was performed using a Gaussian mixture model.

[0082] The Gaussian mixture model assumes that the observed data are generated by a mixture of K Gaussian distributions, and its probability density function is:

[0083] (1);

[0084] Where x is a three-dimensional vector composed of slope, plane curvature, and profile curvature; K is the number of clusters; The mixing weights for the k-th Gaussian component satisfy the following condition: =1; Let be the mean vector of the k-th component; Let be the covariance matrix of the k-th component; Θ represents all model parameters; N(·) is the multivariate Gaussian distribution density function.

[0085] The model parameters are estimated using the expectation-maximization algorithm. In the expectation step, the parameters for each sample are calculated. Posterior probability of belonging to the k-th Gaussian component :

[0086] (2);

[0087] In equation (2), Let be the posterior probability that the i-th sample belongs to the k-th Gaussian component; The mixing weight of the k-th component at the t-th iteration; Let be the mean vector of the k-th component at the t-th iteration; Let be the covariance matrix of the k-th component at the t-th iteration; the denominator is the normalization factor, representing the sample... The weighted sum of probability densities for all K Gaussian components.

[0088] In the maximization step, the parameters are updated based on the posterior probability:

[0089] (3);

[0090] In equation (3), is the mixed weight of the updated k-th component; m is the total number of samples.

[0091] (4);

[0092] In equation (4), Let be the mean vector of the updated k-th component; the numerator is the mean vector of all samples. This is the weighted sum of the weights, and the denominator is the number of valid samples corresponding to that component.

[0093] (5);

[0094] In equation (5), Let be the covariance matrix of the updated k-th component; This represents the transpose of the bias vector; the numerator is the sum of all samples... It is the sum of the outer products of the weighted deviations.

[0095] The expected step and the maximization step are executed alternately, that is, formulas (2) to (5) are repeated until the parameter change is less than the preset threshold or the maximum number of iterations is reached. Finally, each sample is assigned to the cluster with the highest posterior probability.

[0096] The number of clusters K is determined using the Bayesian information criterion. The formula for calculating the Bayesian information criterion is:

[0097] (6);

[0098] Where L is the model likelihood function value, p is the number of model free parameters, and m is the sample size. The BIC values ​​for different values ​​of K are calculated, and the K with the smallest BIC value is selected. Considering the large amount of raster data across the entire area, this embodiment uses spatial stratified sampling to select a portion of the raster data for BIC value calculation, ensuring that the sample proportion of each terrain location is consistent with the overall distribution. In this embodiment, the BIC value is smallest when K=4, therefore the study area is divided into four terrain combination types.

[0099] Spatial distribution and characteristics of four terrain combination types (see appendix) Figure 2 )as follows:

[0100] Type 1: The average slope is approximately 28.3°, the average profile curvature is positive, and the average planar curvature is close to zero. It is mainly distributed on both sides of ridgelines, the upper parts of steep cliffs, and the turning points of convex slopes. Under the influence of gravity, the convex slope surface experiences tensile stress at the top and shear stress concentration at the toe, creating stress conditions conducive to landslides.

[0101] Type 2: The average slope is approximately 27.6°, with a negative average profile curvature and a negative average planar curvature. It is mainly distributed in the lower and middle parts of steep slopes, the sidewalls of gullies, and the catchment areas of concave slopes. The concave shape facilitates surface runoff and groundwater accumulation. Water infiltration leads to increased pore water pressure and reduced effective stress in the soil and rock, weakening the slope's shear strength.

[0102] Type 3: The average slope is approximately 19.2°, and the average profile curvature is negative. It is mainly distributed in the lower and middle parts of moderately sloping surfaces.

[0103] Type 4: The average slope is about 9.8°, and the average profile curvature is close to zero or slightly positive. It is mainly distributed in river valley terraces, plateaus, and gentle slope areas, where the terrain conditions are relatively stable.

[0104] Each grid cell is assigned a terrain combination category label of 1 to 4 as a terrain combination feature.

[0105] Step S4: Input Feature Set Construction

[0106] The terrain combination feature is a categorical variable, with labels 1, 2, 3, and 4 representing only category distinctions. To ensure this feature participates in model training within the same feature space as the continuous influence factors, it undergoes one-hot encoding. The four terrain combination labels are converted into four binary feature columns: if a sample belongs to type one, the corresponding column has a value of 1, and the rest have a value of 0; if it belongs to type two, the corresponding column has a value of 1, and the rest have a value of 0; and so on. After encoding, the terrain combination feature is represented by four binary variables, which, together with the eight continuous influence factors, constitute a 12-dimensional initial input feature matrix.

[0107] Collinearity tests and eliminations were performed on the initial feature matrix. Pearson correlation coefficients were calculated among the factors. Factor pairs with absolute correlation coefficients greater than 0.7 were further screened based on their geological significance in relation to landslide occurrence. When the correlation coefficients between two factors were too high, factors with a more direct physical connection to landslide development were retained. The variance inflation factor of each factor was also calculated, and factors with a variance inflation factor greater than 10 were eliminated. After screening, the final input feature set retained 6 factors.

[0108] Simultaneously, a comparative input feature set excluding terrain combination features was constructed, containing only the aforementioned eight landslide influencing factors, and subjected to the same collinearity removal process. This comparative feature set was used to verify the effect of terrain combination features on improving the model's predictive performance.

[0109] Step S5: XGBoost Model Training

[0110] The training samples were randomly divided into a training set and a test set in a 7:3 ratio.

[0111] An XGBoost model was used to train the input feature set containing terrain features. Hyperparameter optimization employed a grid search combined with five-fold cross-validation. The grid search parameter space was set as follows: learning rate in {0.01, 0.05, 0.1, 0.3}, maximum tree depth in {3, 5, 7, 9}, subsample ratio in {0.6, 0.8, 1.0}, L1 regularization parameter in {0, 0.1, 1}, and L2 regularization parameter in {1, 2, 3}. The optimal hyperparameter combination was determined with the objective of maximizing the area under the receiver operating characteristic (ROC) curve.

[0112] The test set evaluation results show (see appendix) Figure 3 The area under the receiver operating characteristic (ROC) curve of the model with terrain features was 0.913, while it was 0.900 without terrain features. Introducing terrain features increased the ROC curve by 0.013. Simultaneously, accuracy, recall, and F1 score all improved accordingly. These results indicate that terrain features positively contribute to the predictive performance of the XGBoost model.

[0113] The trained model is applied to the entire study area, and the predicted probability of landslide occurrence is calculated grid by grid.

[0114] Step S6: Grouped SHAP Interpretability Analysis

[0115] Based on the trained XGBoost model, the SHAP value of each sample is calculated.

[0116] First, a global feature contribution analysis was performed. The mean of the absolute values ​​of the SAPs of each factor in the test set was calculated and ranked. The results are shown in the appendix. Figure 4 As shown, distance from roads, elevation, and distance from faults rank among the top three in global contribution. The global SHAP mean of terrain combination features is relatively low, and their actual value may be underestimated if only the global ranking is considered.

[0117] Further group contribution analysis was conducted. All raster cells in the study area were divided into four subsets based on terrain combination category labels, and the mean SHAP value of each factor within each subset was calculated. The SHAP values ​​of each subset were calculated based on the same globally trained model, rather than training independent models for each subset. This approach ensures that the comparison of factor contributions across different terrain types is conducted under a completely consistent prediction logic. The group analysis results are attached. Figures 5 to 8 As shown.

[0118] like Figures 5 to 8 As shown, among the four terrain combination types, the SHAP contribution of the three factors—distance from the road, elevation, and distance from the fault—ranks in the top three, but the specific contribution values ​​vary significantly between different types.

[0119] In the Type 1 subset, the mean SHAP values ​​for distance from road, elevation, and distance from fault are 1.35, 0.89, and 0.34, respectively. In the gentle slopes and terraces where this type is distributed, the contributions of each factor are relatively low compared to other types, but the distance from road remains the highest, reflecting that even in relatively flat terrain, road disturbance is a significant contributing factor to landslides.

[0120] In the Type II subset, the mean SHAP values ​​for distance from the road, elevation, and distance from the fault are 1.82, 0.95, and 0.28, respectively. In the medium-slope area of ​​this type distribution, the contribution of the distance from the road factor is significantly higher than that of Type I, and the contribution of the elevation factor is also slightly higher, indicating that as the slope increases, the impact of road construction on slope stability is enhanced.

[0121] In the third subset of type, the mean SHAP values ​​for distance from road, elevation, and distance from fault were 2.13, 1.49, and 0.23, respectively. In the lower and middle parts of steep slopes and the sidewalls of valleys in this type of distribution, the contributions of both distance from road and elevation reached the highest values ​​among the four types, reflecting that the superposition of steep terrain and road engineering activities significantly enhanced the sensitivity to landslides. The landslide susceptibility in this type of area is most prominently controlled by the combined effects of human engineering activities and topographical conditions.

[0122] In the Type 4 subset, the mean SHAP values ​​for distance from the road, elevation, and distance from the fault were 1.77, 0.96, and 0.24, respectively. In the upper regions of ridges and steep cliffs in this type of distribution, the contribution of each factor decreased compared to Type 3, but remained at a high level, indicating that the impact of road cutting disturbance on slope stability in this terrain location is also significant.

[0123] Based on the combined analysis of the four types, the distance to the road factor consistently maintained the highest contribution across all types, reaching a peak of 2.13 in type three. The elevation factor also showed the most significant contribution in type three, reaching 1.49. The contribution of the distance to the fault factor remained relatively stable across all types, fluctuating between 0.23 and 0.34. Using the spatial division framework provided by topographic features, the differences in the intensity of factor effects in different geomorphic units were quantitatively presented. Type three, exhibiting the highest contribution from all factors, should be considered a key area for landslide prevention and control.

[0124] Step S7: Landslide susceptibility mapping

[0125] The predicted landslide occurrence probabilities output by the model in step S5 are spatially visualized. Using the natural discontinuity method, the probability values ​​are divided into five levels: extremely high susceptibility, high susceptibility, moderate susceptibility, low susceptibility, and extremely low susceptibility, generating a spatial distribution map of landslide susceptibility in the study area (see appendix). Figure 9 The spatial distribution of extremely high-risk and high-risk areas overlaps significantly with that of type II and type III terrain combination units, which corroborates the conclusion in the grouped SHAP analysis that the above-mentioned type factors contribute significantly.

[0126] This embodiment verifies the feasibility and effectiveness of the method of the present invention. By incorporating combined terrain features into the XGBoost model training, a measurable improvement in prediction accuracy is achieved; through grouped SHAP analysis based on the same global model, the spatial differentiation of factor intensity under different terrain conditions is revealed. The method's workflow framework is universal and can be extended to landslide susceptibility assessment work in other regions.

[0127] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, and improvements made within the scope of the technical concept and claims of the present invention should be included within the scope of protection of the present invention.

Claims

1. A landslide susceptibility prediction method integrating topographic features and zoning interpretable analysis, characterized in that, Includes the following steps: S1. Obtain digital elevation model data, landslide logging data, non-landslide sample data, and landslide influencing factor data for the study area. The non-landslide sample data is selected from non-landslide areas within the study area. S2. Extract slope factor, plane curvature factor and profile curvature factor based on the digital elevation model data, and standardize the extracted factors. S3. Unsupervised clustering of the standardized slope factor, plane curvature factor and profile curvature factor is performed using a Gaussian mixture model to divide the study area into multiple terrain combination types. Each grid cell is assigned a terrain combination category label as a terrain combination feature. The number of clusters is adaptively determined by the information criterion. S4. Encode the terrain combination features, construct an input feature set by combining the encoded features with landslide impact factor data, and remove collinearity factors from the input feature set. S5. Input the input feature set into the XGBoost model for training, and optimize the hyperparameters of the XGBoost model using the grid search method to obtain the trained landslide susceptibility prediction model. S6. Based on the landslide susceptibility prediction model that has been trained, the SHAP method is used to calculate the global feature contribution. The grid cells of the study area are grouped according to the terrain combination type. The SHAP value of each landslide influencing factor in each group is calculated respectively. The spatial differentiation law of factor contribution under different terrain combination types is analyzed. S7. Based on the output of the landslide susceptibility prediction model, generate a spatial distribution map of landslide susceptibility in the study area.

2. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The number of clusters mentioned in step S3 is determined by the Bayesian information criterion, and the number of clusters with the smallest Bayesian information criterion value is selected as the number of the multiple terrain combination types.

3. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The multiple terrain combination types mentioned in step S3 are automatically clustered by the Gaussian mixture model based on the statistical distribution characteristics of slope factor, plane curvature factor and profile curvature factor. Each type has significant differences in slope value range and positive and negative curvature combination form, and corresponds to different topographic units in the study area.

4. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The clustering process of the Gaussian mixture model in step S3 includes: assuming that the terrain factor data is generated by mixing multiple Gaussian distributions, estimating the mean vector, covariance matrix and mixing weight of each Gaussian component through the expectation-maximization algorithm, calculating the posterior probability of each grid cell belonging to each cluster, and taking the cluster corresponding to the maximum posterior probability as the terrain combination category of the grid cell.

5. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The encoding process described in step S4 employs at least one of one-hot encoding, label encoding, or target encoding to convert terrain combination category labels into numerical features, enabling categorical features and continuous influencing factors to participate in model training within a unified feature space.

6. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The collinearity factor removal process in step S4 uses variance inflation factor analysis or Pearson correlation coefficient analysis to remove factors with strong collinearity.

7. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The hyperparameter optimization of the XGBoost model in step S5 includes performing grid search optimization on at least one of the learning rate, maximum tree depth, subsample ratio, and regularization parameter.

8. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... Step S6, the SHAP analysis grouped by terrain combination type, includes: dividing all raster cells in the study area into multiple subsets according to the terrain combination category label; calculating the mean SHAP contribution of each landslide influencing factor for each subset; and determining the dominant influencing factor under each type by comparing the contribution differences of the same factor among different subsets. The SHAP values ​​of each subset are calculated based on the same global model trained in step S5, so that the comparison of factor contributions between different terrain types is carried out on the same prediction logical scale.

9. The landslide susceptibility prediction method according to claim 1, which integrates topographic features and zoning interpretable analysis, is characterized in that... The landslide influencing factor data are selected from at least one of geomorphological factors, geological structure factors, hydrological environmental factors, vegetation cover factors, and human engineering activity factors.

Images