Interpretable machine learning incorporating major lithology for regional landslide warning in northern and eastern Guangdong

Exploratory data analysis

The data analysis focuses on investigating the statistical distributions and relationships of geotechnical and geological parameters with lithology to establish its significance as an independent factor in modeling.

To examine the data comprehensively, we began by analyzing the distribution of each feature (Fig. 1). Histograms revealed that most features exhibit unimodal distributions, although their symmetry and spread vary. For example, features such as water content and internal friction angle demonstrate near-symmetric distributions, suggesting quasi-normal behavior. However, certain parameters, including compression coefficient and slope length, exhibit pronounced skewness, indicating potential non-linear relationships or the presence of outliers. This variability highlights the need for more advanced approaches to ensure robust modeling.

The distributions of each feature under two major lithologies are shown in Fig. 2. A comparative analysis of these distributions highlights significant variability in geotechnical properties and geological metrics between igneous and sedimentary lithologies. To quantify lithology-dependent contrasts, we applied two-sample Wilcoxon rank-sum tests (for medians) and Student’s t-tests (for means) to each non-rainfall feature, comparing igneous and sedimentary groups. These formal tests yielded large test statistics with p-values < 0.05 for all features (most ≪ 0.001), supporting that the differences observed in Fig. 2 are statistically significant and reflect robust lithological effects rather than random variation.

Sedimentary rocks exhibit higher water content due to secondary porosity from fractures and weathering, which enhances their permeability and water retention. Their lower compression coefficient and cohesion are linked to weaker cementation and a more granular structure, leading to lower shear strength and compressive resistance³³. In contrast, igneous rocks typically have lower water content, not only due to their dense crystalline structure but also because of fewer fractures, although weathering can create void spaces that affect permeability. The data also show higher plasticity index and internal friction angle, reflecting the strong interlocking of mineral crystals in granitic rocks, which enhances their mechanical strength and resistance to deformation³⁴. These findings also suggest the reason why sedimentary rock landslides are more frequent than igneous rock events.

Slope geometry further reflects lithological influences, with sedimentary formations associated with longer, more confined slopes and igneous formations forming steeper, broader slopes due to differences in material strength and erosion resistance. Aspect distributions, however, appear less affected by lithology, likely influenced more by regional tectonics or climate. These findings underline the necessity of using lithology as an additional prerequisite to differentiate models. By accounting for the specific mechanical behaviors and morphological patterns associated with each lithology, such models can more accurately capture the interactions between material properties and failure mechanisms.

To assess the independence and interrelationship of the features, we assumed all variables follow independent continuous Gaussian distributions and then calculated Pearson correlation coefficients for each pair of features (Fig. 3). Pearson correlation measures the strength and direction of linear relationships between variables, with values ranging from −1 (perfect negative correlation) to 1 (perfect positive correlation). The results indicated that most features exhibited weak correlations (absolute values below 0.3), while a few pairs showed moderate correlations (absolute values between 0.3 and 0.5).

We also calculated the variance Inflation Factor (VIF) for each input variable to ensure the independence of the predictors. VIF is a statistical measure used to assess the degree of multicollinearity in regression models. As a rule of thumb, a VIF value greater than 10 suggests significant multicollinearity. We found that all variables had a VIF value below 10, indicating that multicollinearity is not a significant issue in our work. While some moderate multicollinearity (5 < VIF < 10) was detected among certain features, such as slope and water content, this is expected in terrain modeling and does not compromise the model’s validity or performance, especially given the robustness of tree-based algorithms. The findings suggest that these features can be considered independent, providing a robust foundation for subsequent modeling and offering insights into potential mechanisms underlying landslide occurrences.

Comparison of machine learning models

To select the most appropriate method for our landslide prediction model, we compared multiple ML models, including Logistic Regression (LR), Support Vector Machine (SVM), Random Forest (RF), and Artificial Neural Network (ANN). Here, we adopt a broad definition of ML, under which LR is included as a baseline classifier for binary prediction tasks. At this stage, lithology was not differentiated but treated as a categorical feature encoded into the dataset for training. Since some models require data normalization, we applied Z-score standardization to ensure all features were on a comparable scale.

We use the ROC curve to evaluate the models’ performance (Fig. 4). The RF model outperformed all others with an AUC of 0.98, indicating its superior performance. The ANN model followed closely behind, with the SVM model showing competitive results as well. However, the LR model performed the weakest across all metrics. This outcome suggests that while LR is simple and easy to use, it struggles with more complex problems compared to advanced algorithms. The RF model’s robust handling of complex data and high-dimensional features makes it the most suitable choice for this task. Given the promising results of SVM and ANN, these models could also be useful for highly demanding classification tasks, but based on our comparison, we chose RF to enhance the accuracy and reliability of landslide early warning predictions.

Model training and results

For igneous lithology, we randomly divided the dataset of 2184 entries into an 80% training set and a 20% testing set. Utilizing the Bayesian optimization algorithm combined with ten-fold cross-validation, we determined the optimal set of hyperparameters. The default output threshold for the RF classifier is set at 0.5. Considering the specific policies of our study area, it is prudent to issue warnings and take preventive measures when the relative probability of a landslide occurring exceeds 0.4. Therefore, we adjusted the threshold to 0.4 to match the threshold specified by the warning center. This adjustment means that when the model predicts the probability of an event occurring to be over 0.4, it is classified as a positive instance (indicating a prediction of occurrence). Conversely, predictions falling below this threshold are considered negative instances (indicating a prediction of non-occurrence). This adjustment enhances the model’s relevance and responsiveness to the specific risk tolerance and safety requirements of the study area.

The efficacy of the RF model in predicting landslides based on an integrated analysis of precipitation, geological environment, and geotechnical mechanics features is underscored by the confusion matrix (Fig. 5a). This matrix visualizes the predictive accuracy of machine learning models, clearly depicting their performance in classifying events correctly. Based on the confusion matrix, we derived several key numerical metrics (Fig. 5b). On the testing dataset, the model achieved a f1 score of 0.9071, an accuracy of 96.11%, a hit rate of 94.3%, a miss rate of 5.7%, and a false alarm rate of only 3.4%. This indicates that the model performs well, with only three misses and four false alarms in numerical terms.

In addition to the confusion matrix, the Receiver Operating Characteristic (ROC) curve (Fig. 5c) and the Precision-Recall (PR) curve (Fig. 5d) provide further insights into the model’s predictive capabilities. With an Area Under the ROC Curve (AUC) of 0.989 and an Average Precision (AP) of 0.965, the model demonstrates excellent generalization ability. In summary, the model has been well-trained and exhibits robust performance.

For sedimentary lithology, we also randomly divided the dataset of 4899 entries into an 80% training set and a 20% testing set. Using the same approach, we determined the optimal set of hyperparameters. On the testing dataset, the model achieved a f1 score of 0.9114, an accuracy of 97.14%, a hit rate of 90.6%, a miss rate of 9.4%, and a false alarm rate of 1.6%, resulting in only three missed and twelve false alarms numerically. Additionally, the model performed well with an AUC of 0.992 and an AP of 0.974 (Fig. 6). Both models, for igneous and sedimentary lithologies, demonstrated satisfactory performance, indicating their readiness for practical application in landslide prediction and early warning systems.

Global Interpretability of RF model

The Random Forest (RF) algorithm includes a built-in feature importance measure³⁵, which quantifies the contribution of each feature to the model’s predictive power. This measure ranks predictive variables based on their usefulness in making accurate predictions, thereby providing valuable insights into the key factors contributing to landslides (see Fig. 7). This global explanation helps to understand the overall contribution of each feature across all predictions. Feature importance is calculated based on the decrease in impurity caused by each feature across all trees in the forest. Features that result in larger decreases in impurity are considered more important. This metric not only aids in understanding the model but also in feature selection and optimization by highlighting which features are most influential in predicting the target output.

SHAP also provides a similar feature importance ranking through the use of summary beeswarm plots (Fig. 8). These plots offer an intuitive visualization of how each feature contributes to the model’s predictions. In a beeswarm plot, the horizontal axis represents the SHAP values of each sample for the corresponding feature. A SHAP value greater than 0 indicates that the feature has a positive contribution to predicting a landslide, while a SHAP value less than 0 suggests a negative contribution to the prediction. The color encoding in the plot indicates the feature value, with colors ranging from low to high values, providing insight into how the magnitude of the feature value influences the prediction. This visualization helps to identify not only which features are important but also how their specific values affect the model’s predictions in different samples.

The discrepancy between the feature importance rankings in Fig. 7 (RF model) and Fig. 8 (SHAP analysis) arises from the different mechanisms used by these two techniques. RF importance provides a global view of feature significance based on overall model performance, while SHAP values offer a local, instance-level explanation, revealing how individual feature values contribute to specific predictions. Despite slight differences in the ranking, the most important 8 to 10 features identified by both methods are largely consistent, reflecting a similar overall mechanism. The use of SHAP for generating PDPs in the next section was specifically chosen because it provides a more detailed, localized understanding of feature contributions, which is crucial for interpreting the varying influence of features under different conditions.

SHAP decision heatmaps are another powerful global explanation tool that helps visualize the relationship between feature contributions and model predictions across the entire dataset (Fig. 9). In a SHAP decision heatmap, the horizontal axis represents the test set samples for landslide data after hierarchical clustering, grouping similar samples together. The color encoding indicates the magnitude of the SHAP values for each feature, with warmer colors representing higher SHAP values and cooler colors indicating lower contributions. The vertical axis displays the various features, while the baseline, denoted as f(x), represents the model’s default prediction before considering the specific feature contributions. If the resulting output curve is above the baseline, the model predicts a landslide, otherwise the model predicts non-landslide.

The heatmap shows how each feature’s contribution varies across the testing set, with warmer colors indicating higher contributions to landslide prediction. The figure highlights the importance of specific features, such as precipitation and geotechnical factors, and shows how these features interact with one another in influencing the model’s output. While precipitation is a dominant factor in predicting landslide occurrences, non-precipitation factors serve as important modifiers that can adjust the risk levels based on the unique characteristics of the terrain, geotechnical properties, and slope conditions. This will highlight the complexity of the model and the interplay between various factors influencing landslide predictions.

Precipitation factors unsurprisingly emerged as the most influential features in our model, with rainfall on the day of the event being significantly more important than other variables. The combined importance of precipitation features over a three-day span account for more than half of the total importance. This underscores the dominance of rainfall in triggering landslides in the studied region. However, it is equally important to note that non-precipitation factors also played a crucial supportive role. These factors were able to significantly adjust the influence of precipitation in certain cases, correcting for situations where the rainfalls alone overstated or understated the landslide risk.

Sedimentary rocks typically exhibit higher porosity and permeability compared to igneous rocks, making them highly susceptible to infiltration during rainfall events. This infiltration can elevate pore water pressure and significantly reduce effective stress, leading to a decrease in shear strength and eventual slope failure³⁶. Moreover, the sensitivity of sedimentary rocks to precipitation can be linked to their composition, which often includes clay-rich materials prone to swelling and softening under wet conditions³⁷.

In contrast, igneous rock landslides show a slightly lower cumulative importance for precipitation factors. The relatively low porosity and higher intact strength of igneous rocks diminish the direct influence of rainfall on slope stability. Instead, precipitation’s effect is more indirect, such as through increased slope weight or localized saturation zones. Nevertheless, rainfall continues to play a significant role in igneous rock landslides. We propose that the higher degree of weathering observed in the igneous rocks within the study area have enhanced their sensitivity to rainfall, resulting in a greater vulnerability to precipitation-induced instability than would typically be observed in less weathered, more intact igneous materials³⁸.

The mechanical stability of igneous slopes is predominantly governed by their strength and structural integrity. Water content is particularly influential, as localized saturation can weaken discrete planes of weakness or fault zones³⁹. Additionally, terrain attributes such as slope height and aspect directly affect the gravitational driving forces, making them critical in high-steep terrain typically associated with igneous rocks. The role of internal friction angle and cohesion further reflects the strength-dependent nature of igneous rock landslides, where failure occurs primarily through structural or planar mechanisms. The relatively stable structure of igneous rocks contributes to the confidence and stability of model predictions. This structural consistency results in less variability in the model’s output, as the predictable nature of the mechanical properties in igneous slopes reduces uncertainty in landslide susceptibility.

For sedimentary rocks, slope length, in particular, is a key factor as it determines the extent of water infiltration and accumulation along the slope, which can destabilize the material over time. The plasticity index’s high importance is indicative of the clay-rich composition of many sedimentary rocks, which exhibit significant shear strength reduction when wet. The compression coefficient reflects the compressibility of sedimentary materials under loading or wetting, further contributing to slope failure. These factors collectively suggest that sedimentary landslides are strongly influenced by hydrological processes and material properties that control their response to water.

These findings reveal distinct landslide mechanisms in our study region: igneous rock landslides are predominantly terrain- and strength-controlled, with failure modes involving structural instability, while sedimentary rock landslides are more water-driven, with mechanisms influenced by hydrological processes and material weakening.

Local Interpretability from features

To further interpret the model, we perform local Interpretability at the feature levels. The three most influential non-precipitation features will be focused on using PDP, which illustrate the relationship between each feature and landslide occurrence, as well as their interactions.

For igneous model, the three key features are water content, cohesion, and slope height. A decrease in water content below 26% is associated with a marked increase in landslide dependency, suggesting that low water content increases the material’s apparent cohesion. Higher cohesion (above 29 KPa) and increased slope height (above 12 meters) are linked to a significantly higher likelihood of landslides, as shown in Fig. 10. These results highlight the importance of the degree of weathering in the igneous rocks of the study area. The higher apparent cohesion in weathered igneous rocks, which is primarily due to matrix suction or weak cementation, can appear substantial under dry conditions but significantly decreases upon wetting or disturbance, making the slope more susceptible to failure. Additionally, the combination of slope height and cohesion plays a critical role in the landslide risk, especially when the material’s strength is reduced or the slope is very steep⁴⁰.

The sedimentary rock PDP reveals a different set of relationships (Fig. 11). The key features here are slope length, aspect, and slope. Slope length and aspect are both negatively correlated with landslide dependency. Specifically, when slope length exceeds ~25 meters, the aspect of the slope becomes the dominant factor influencing landslide risk. Aspect is primarily determined by the terrain’s orientation, which governs water drainage and sun exposure—two factors that significantly affect slope stability. Slope angle also plays a crucial role in landslide prediction with relatively high dependency throughout the entire range. For steep slopes, the risk of failure increases a bit due to potential sliding along bedding planes or jointed structures, where the potential slip surface is less steep than the slope itself, facilitating sliding. For gentler slopes, the risk slightly increases due to the formation of perched water tables on low-permeability layers, which leads to an increase in pore pressure and material weakening. This can result in translational sliding along weak or residual failure surfaces, contributing to the observed landslide dependency⁴¹.

These findings clearly demonstrate that non-precipitation features have a greater influence on landslide prediction in igneous rocks than in sedimentary rocks. This further underscores the idea that igneous rock slopes exhibit a lower dependency on precipitation compared to sedimentary rock slopes, highlighting the distinct behaviors of these two rock types in landslide susceptibility.

Graded warning model performance comparison

Since 2003, China has progressively implemented regional landslide meteorological warning services, achieving notable success in disaster prevention and mitigation⁴². In alignment with the Technical Guidelines for Meteorological Risk Early Warning of Geological Hazards at the Municipal Level in Guangdong Province, we classify the landslide warning probabilities based on model predictions. For an output probability P ≤ 20%, no warning is issued (indicated as a gray warning); for P > 20% and P ≤ 40%, a blue warning is issued; for P > 40% and P ≤ 60%, a yellow warning is dispatched; for P > 60% and P ≤ 80%, an orange warning is released; and for P > 80%, a red warning is proclaimed, signaling the highest level of alert.

To demonstrate the advantages of our approach, we compared it with two traditional statistical models: the implicit statistical warning model and the explicit statistical warning model. The implicit model categorizes precipitation into different ranges by calculating early effective precipitation amounts:

where \({R}_{p}\) is the amount of precipitation that impacts landslides during the precipitation process prior to their occurrence, \({R}_{n}\) is the daily precipitation on the nth day before, \(n\) is the number of effective precipitation days, and \(k\) is the effective precipitation coefficient set to 0.84. Corresponding warnings are issued based on the range within which the precipitation falls⁴³. Building on this, the explicit model incorporates geological environmental factors and employs Principal Component Analysis (PCA) to extract main components that account for a total contribution rate of 90%. It calculates a potentiality score for each landslide point^44,45. These scores, combined with effective precipitation, are then used to delineate warning levels:

where \(G\) is the potentiality score from geological environment, \({R}_{d}\) is the daily precipitation on the day of the landslide occurrence, and \({R}_{p}\) is the early effective precipitation amount.

In evaluating our model’s effectiveness, we conducted a comprehensive comparison across the testing dataset against traditional statistical methods, which utilize all landslide points for analysis. To facilitate a nuanced examination, tests were conducted separately on positive and negative cases across different lithologies, calculating hit rates and false alarm rates for each. It is noteworthy that blue warnings, typically issued online without subsequent action, are classified as non-landslide indicators in our analysis. The results (see Fig. 12) reveal that our RF model significantly outperforms the traditional statistical models in both igneous and sedimentary rock regions, delivering higher hit rates (94.32% and 90.57%) and lower false alarm rates (3.44% and 1.58%). Furthermore, the accuracy in tier classification is markedly improved; for positive cases, our model predominantly issues red or orange warnings, indicative of an imminent landslide. Conversely, for negative cases, it almost exclusively issues gray or blue warnings, denoting a negligible risk of occurrence.

Application on rainfall-induced landslides in Pingyuan County in June 2024

Pingyuan County, located in northeastern Guangdong, China, experienced severe rainfall-induced landslides on June 16, 2024. This event led to significant fatalities and economic losses, with rainfall reaching up to 360 mm in a single day. The region’s steep terrain, granite lithology, and prolonged rainfall contributed to the widespread slope failures, highlighting the vulnerability of the area to monsoon-driven landslides⁴⁶.

Pingyuan County is situated within our study area, and there are three available rainfall gauge stations in the region (Fig. 13a). For this application, we utilized daily rainfall data from June 1 to June 16, 2024, from these stations. Using the corresponding lithology-based RF model, we conducted landslide predictions for this period by a 5-day window to assess the practical feasibility of our method in a real-world setting. The model issued daily warnings from June 5, confirming its potential to predict landslides and assist in early disaster prevention.

During the period leading up to the landslides, when rainfall had not yet reached extreme levels, the model did show a small number of false alarms (Fig. 13b). These occurred due to slight increases in rainfall on the days when the 5-day sliding window was updated. However, it is important to note that these false alarms were all within the yellow warning level, meaning they were relatively low-risk predictions and were quickly corrected as the rainfall data progressed.

On June 16, as rainfall levels surged, all three stations issued warnings at or above the orange level, which closely aligned with the significant increase in precipitation on that day. Notably, the final warning levels and probability values varied across the stations, and upon closer examination, we found that these differences corresponded well with the number of landslides within the warning zone of each station. This suggests that the model was effectively accounting for both the intensity of the rainfall and the local topographic factors, enhancing its ability to provide tailored and reliable early warnings for landslides. This case study serves as a validation of our approach and demonstrates its applicability in real-world landslide forecasting.