Machine learning approaches to debris flow susceptibility analyses in the Yunnan section of the Nujiang River Basin

The Random Forest model

The Random Forest model is effective in capturing and simulating complex nonlinear relationships between debris flows and evaluation factors, and can handle large-scale, high-dimensional debris flow datasets without overfitting. Furthermore, RF indicates the relative importance of each evaluation factor, guiding the understanding of which factors have the greatest impact on debris flow susceptibility. It also demonstrates good adaptability to noise and outliers in the data, making it less susceptible to interference. Therefore, RF exhibits excellent applicability to assessments of debris flow susceptibility (Duan et al., 2022; Zhang & Wu, 2019).

RF is an integrated algorithm consisting of multiple unrelated decision trees that determine the DFS, in which the final output is determined by all decision trees together, and is defined as Breiman (2001),

(1) $\begin{array}{c}H\left(x\right)=\arg \begin{array}{c}max\\ z\end{array}\sum _{i=1}^{k}I\left(h_i\left(x\right)=Z\right)\end{array}$where $H\left(x\right)$ denotes the results of the model’s predicted DFS for each watershed unit, $h_i\left(x\right)$ denotes the $i$th decision tree, x denotes attributes, $h_i\left(x\right)=Z$ is the prediction of variable $Z$ using the $i$th tree in variable x, and $I\left(\text{.}\right)$ is the prediction of each decision tree.

Given a database, RF can be interpreted via the following three steps. Firstly, sample subsets are extracted using the Bootstrap resampling method. Specifically, n sample subsets of the same size as the original sample are extracted using the put-back method. Secondly, a decision tree is constructed for each sample subset. Among the attributes of the sample subset, k attributes are randomly selected. Then, the best partitioning attributes of the nodes between decision trees are selected based on the Gini Index, which is calculated as

(2) $\begin{array}{c}Gini\left(p\right)=\sum _{k=1}^k{p}_k\left(1-p_k\right)=1-\sum _{k=1}^k{p_k}^2\end{array}$where $p_k$ indicates the probability that the selected sample belongs to category $k$. Smaller Gini Indexes mean that the probability of a selected sample in the set being misclassified is smaller. Finally, n decision trees are combined to generate a random forest (Fig. 3).

Hyperparameters of the RF model can be divided into two categories: those that determine the sampling method, such as bootstrap and the number of classifiers that determine the sampling method, and the number of decision trees. Parameters that determine the decision tree, such as maximum depth (max_depth), minimum number of samples for a leaf node (min_samples_leaf), minimum number of samples required to split an internal node (min_samples_split), the maximum number of features randomly selected as candidates for splitting (max_features), and a criterion that determines the maximum depth, minimum number of samples for a leaf node, minimum number of samples required to split an internal node, the maximum number of features randomly selected as candidates for splitting, and a criterion for the optimal split attribute.

The support vector machine

Support vector machine (SVM) demonstrates superior classification performance on unseen data due to its outstanding generalization capability, laying a crucial foundation for the credibility and practicality of the model in real-world applications. In practical applications, areas prone to debris flows encompass complex environmental features, and SVM’s excellent handling of high-dimensional data provides robust support for assessing susceptibility considering multiple evaluation factors. Additionally, by employing various kernel functions, SVM exhibits flexibility in modeling nonlinear relationships in high-dimensional space, thereby enhancing its applicability in the assessment of DFS. These attributes position SVM as a powerful tool when facing real and complex datasets, providing a reliable analytical framework for accurately evaluating DFS.

There are two cases of linearly divisible and linearly indivisible sample data in the feature space (Fig. 4). As an example of binary data, a binary classification space $D\left(X_i,Y_i\right)$, $i=1,\dots ,l.X_iϵR_n,Y_iϵ\left\{1,-1\right\},$ where, $l$ represents the number of samples, and n denotes input dimensionality. The hyperplane $\omega x+b=0$ can be found in the original space when the sample data are linearly divisible separating the two classes of samples completely. When the sample data are linearly indivisible, it is necessary to perform nonlinear mapping $\mathrm{\Phi }\left(x\right)$, mapping it from the input space to a certain feature space, the classification hyperplane can be expressed as $\omega \mathrm{\Phi }\left(x\right)+b=0;$ meantime, the optimal hyperplane that requires $2/\omega$ is the largest, and the problem is transformed into a quadratic programming problem, with the application of the Lagrange multiplier method for the solution, namely,

(3) $\begin{array}{c}\{\begin{array}{c}min{\displaystyle \frac{{\omega }^2}{2}+C\sum _{i=1}^l{\epsilon }_i,}\\ s.t.y_i\left(\omega \cdot x_i+b\right)\ge 1-{\epsilon }_i,\\ {\epsilon }_i\ge 0,i=1,2,...,l\end{array}\end{array}$where ${\epsilon }_i$ is the slack variable and C is the penalty factor. According to the Kuhn-Tucker (K-T) condition, the following dyadic problem can be obtained:

(4) $\begin{array}{c}\begin{array}{c}max\sum _{i=1}^l{a}_i-{\displaystyle \frac{1}{2}\sum _{i=1}^l\sum _{j=1}^l{a}_i{a}_j{y}_i{y}_j\phi \left(x_i\right)\cdot \phi \left(x_j\right),}\\ s.t.0\le a_i\le C,\sum _{i=1}^l{a}_i{y}_i=0\end{array}\end{array}$

By solving the dyadic problem of this quadratic programming, the discriminant function is obtained as:

(5) $\begin{array}{c}f\left(x\right)=sign\sum _{i=1}^l{a}_i{y}_i\left[\mathrm{\Phi }\left(x_i\right)\cdot \mathrm{\Phi }\left(x_j\right)\right]+b\end{array}$

According to the relevant theory of generalized functions, as long as a kind of kernel function $K\left(x_i,y_i\right)$ satisfies the Mercer condition, it corresponds to an inner product in a certain transformed space, $K\left(x_i,y_i\right)$ = $\varphi (x_i)\cdot \varphi (x_j)$, and linear classification of a certain nonlinear transformation can be achieved by using a different inner product function in the optimal classification surface (Suykens & Vandewalle, 1999).

Linear kernel, polynomial kernel, sigmoid kernel, and the radial basis function (RBF) are commonly used in support vector machine models, where

Linear kernel

(6) $\begin{array}{c}K\left(y,y^{\mathrm{\prime }}\right)=y^T{y}^{\mathrm{\prime }}\end{array}$

RBF

(7) $\begin{array}{c}K\left(y,y^{\mathrm{\prime }}\right)=\exp \left(-{\displaystyle \frac{1}{2{\sigma }^2}y-{y^{\mathrm{\prime }}}^2}\right)\end{array}$where $y$ and $y^{\mathrm{\prime }}$are both basis vectors in the feature space, and $\sigma$ is a model’s hyperparameter. Compared with the linear kernel, the RBF can transform the features dimensionality for reducing computational complexity, which is extremely suitable for predicting DFS in high-dimensional feature spaces (Lin & Lin, 2003). The penalty parameter C, an empirical parameter in the SVM model, is used to control the tolerance of systematic outliers, allowing for a few outliers to exist in the opponent classification. A higher value of the penalty parameter leads to fewer outliers in the opponent classification. What’s more, the radial basis function kernel has an additional kernel parameter γ i.e., kernel bandwidth to be optimized, where γ = 1/2σ². As γ increases, the fit changes towards non-linear.

Accuracy evaluation metrics

Evaluation of model’s accuracy is critical for decision-makers and relevant institutions. High-precision model predictions contribute to more precise decision-making, ensuring that measures taken are scientifically sound and effective. Comparing the accuracy of different models in practical applications contributes to the selection of the most suitable model for a given region or topography, thereby bolstering the credibility of predictions. Additionally, accuracy evaluation provides feedback on the current model performance, guiding continuous improvement efforts.

Accuracy, Precision, Recall, Kappa, F1-score, receiver operating characteristic (ROC) curve, and area under ROC (AUC) are used as evaluation metrics. Accuracy represents the proportion of correctly classified debris flow samples, serving as a key indicator for overall model performance. Precision refers to the proportion of samples with positive cases that are correctly predicted which are critical for reducing false positives and ensuring the rational use of limited resources. Recall indicates the proportion of positive cases that are correctly predicted in the true sample, and is essential for minimizing the risk of overlooking potential hazard zones. F1-score offers a balanced assessment of precision and recall, guiding the establishment of reasonable warning and management strategies. A higher F1-score indicates greater model accuracy. Kappa measures model consistency, indicating its ability to make similar judgments under different conditions or at different times. In the context of dynamic changes in debris flow risk, model stability is essential for providing continuous and effective risk assessments. A higher Kappa means greater classification accuracy. The ROC curve aids decision-makers in balancing sensitivity and specificity by using true positive rate (TPR) as the vertical axis and false positive rate (FPR) as the horizontal axis. The AUC, obtained by integrating the ROC curve, reflects the model’s classification effectiveness. Even in situations with imbalanced positive and negative samples, a higher AUC value indicates superior model accuracy. Their definitions are as follows,

(8) $\mathrm{A}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}=\frac{\mathrm{T}\mathrm{P}+\mathrm{T}\mathrm{N}}{\mathrm{T}\mathrm{P}+\mathrm{F}\mathrm{P}+\mathrm{T}\mathrm{N}+\mathrm{F}\mathrm{N}}$

(9) $\begin{array}{c}Precission={\displaystyle \frac{TP}{TP+FP}}\end{array}$

(10) $\begin{array}{c}Recall={\displaystyle \frac{TP}{TP+FN}}\end{array}$

(11) $\begin{array}{c}F1-score={\displaystyle \frac{2\ast Precision\ast Recall}{Precision+Recall}}\end{array}$

(12) $\begin{array}{c}Kappa={\displaystyle \frac{Accuracy-p_e}{1-p_e}}\end{array}$

(13) $\begin{array}{c}{p}_e={\displaystyle \frac{\left(TP+FP\right)\ast \left(TP+FN\right)+\left(FN+TN\right)\ast \left(FP+TN\right)}{\left(TP+FP+FN+TN\right)\ast \left(TP+FP+FN+TN\right)}}\end{array}$

The ROC curve is generated using the true positive rate (TPR) as the vertical axis and FPR as the horizontal axis. The AUC is obtained by integrating the ROC curve and it reflects the classification effect of the model. The closer its value to 1, the better and more accurate the model. The calculation is as follows,

(14) $\begin{array}{c}TPR={\displaystyle \frac{TP}{TP+FN}}\end{array}$

(15) $\begin{array}{c}FPR={\displaystyle \frac{FP}{FP+TN}}\end{array}$where in Eqs. (8)–(15), TP stands for the true positive rate, FP represents the false positive rate, TN signifies the true negative rate, and FN denotes the false negative rate.

To facilitate horizontal comparisons of model predictions, many studies use the equal spacing method to classify DFS into five zones and confirm the method’s applicability (Liang et al., 2020; Liu, Miao & Tian, 2017). Therefore, in this study, we divided the predicted debris flow susceptibility calculated by the three models into five classes, from low to high, corresponding to the very low susceptibility zones (0 to <0.2), low susceptibility zones (0.2 to <0.4), medium susceptibility zones (0.4 to <0.6), high susceptibility zones (0.6 to <0.8) and very high susceptibility zones (0.8 to <1), respectively (Zhang & Wu, 2019).

Evaluation units

Raster cells and watershed units in DEM are commonly used for susceptibility assessments (Zou et al., 2017). Raster cells are more convenient for modelling and calculation because of their regular shape and uniform size, while watershed units can represent integrated geomorphological characteristics of hydrological processes, and that helps in obtaining the actual conditions of debris flow (Liu et al., 2018; Qiang et al., 2019; Zhang et al., 2022). Therefore, we adopted watershed units as the basic evaluation units and used ArcMap’s hydrological analysis tools to categorize the 30 m spatial resolution DEM data of the YNR Basin. Finally, the study area was divided into 1,070 watershed units (Esri, 1981).

Evaluation factors

The formation of debris flows is determined by a combination of factors (Huang et al., 2022). After conducting a thorough investigation, data collection, and preliminary analysis of existing data on the historical background, geological structure, topography and geomorphology, hydro-meteorology, soil and vegetation, and human activities in the study area, we selected evaluation factors from five aspects: topographic, rainfall, geological, and vegetation conditions, and human activities (Table 2).

(1) Topographic conditions. Topographic conditions are critical to the formation of debris flows. We chose relative elevation difference and slope to represent potential energy of watersheds and the ability to carry rocky soil, respectively, and watershed area for runoff and sediment yield calculations based on previous research in this area by Liu & Tang (1995), Sun et al. (2021). Therefore, based on the DEM with a spatial resolution of 30 m, the mean relative elevation difference and average slope of each watershed unit were calculated using the ArcMap function of zonal statistics as a table, and the area of each watershed unit was extracted using the “calculate geometry” function (Esri, 1981).

(2) Rainfall conditions. Rainfall is necessary for debris flow incubating and triggering (Cui, Yang & Chen, 2003; Liu, Miao & Tian, 2017; Xu, 2016). Antecedent rainfall serves mainly to wet or soften the soil and reduce the stability of rocky soil. Short-duration heavy rainfalls create a strong mechanical impact on the soil that is almost saturated, and disrupt the equilibrium of the slope causing debris flows (Pan et al., 2012; Tan, Yang & Shi, 1990; Zhang & Guo, 2021). To address the effect of rainfall, we selected three factors to characterize the triggering effect of heavy rainfall on debris flows in 2020: the number of heavy rainstorms, the number of rainstorms, and the number of heavy rains. Based on the distinctive interannual variability in precipitation distribution in the study area between dry and rainy seasons, we chose the average rainfall during the rainy season (April to September) to characterize the effect of antecedent precipitation on developing debris flows. Therefore, a total of 183 daily precipitation data were extracted for dates from 2020/04/01 to 2020/09/30 in the study area, and the field calculator and the regional statistical function of ArcMap was used to calculate the average rainfall for the 2020 rainy season for each watershed unit. Using ArcMap’s model builder (Esri, 1981), a total of 366 daily precipitation data for 2020 were sequentially screened for the number of heavy rainstorms with a cumulative daily precipitation of 100–250 mm, the number of rainstorms with 50–100 mm, and the number of heavy rains of 25–50 mm; then, we summed the number of days in compliance with the raster cells using the field calculator, and calculated the average number of heavy rainstorms, rainstorms, and the heavy rainfall for each watershed unit with the ArcMap function of zonal statistics table (Esri, 1981).

(3) Geological conditions. Fracture zones affect the continuity and stability of rocky soil, and surface soil provides a rich sediment source for debris flows (Pham et al., 2016). Consequently, we used fracture zone density and soil texture to characterize the influence of geological conditions on debris flows. Fracture zone density was calculated by dividing the length of the fracture zone in each watershed unit by the watershed area. Soil texture was calculated separately using the average content of clay, silt, and sand within the watershed unit and the function of zonal statistics.

(4) Vegetation conditions. Plant roots stabilize the rock and soil mass, and increase soil resistance to erosion (Huang et al., 2022). To some extent, plant roots inhibit erosion and hinder the sliding of the topsoil (Zhao, Wu & Wang, 2006). We used the ArcMap function of zonal statistics to calculate the average Normalized Difference Vegetation Index of each watershed unit to characterize vegetation cover (Esri, 1981).

(5) Human activities. Geological structure and surface become unstable with land use changes, which also generate loose deposits providing material sources for debris flows (Huang et al., 2022; Tien Bui et al., 2017; Xu, 2016). Road construction can indicate regional land development intensity. We used the ArcMap function to calculate land use types with the largest proportions for each watershed unit. Densities of highway density, railway, urban primary roads, urban secondary roads, urban tertiary roads, urban quaternary roads, and county and town roads were calculated using corresponding road type lengths divided by the watershed area (Esri, 1981).

Data pre-processing

We used the Random Forest (RF) model to reduce the noise in the data sets and to select evaluation factors following the approach of Kursa & Rudnicki (2011).

First, a RF model was built with all the evaluation factors in Python after digitizing, data formatting, and unifying georeferencing. Second, the model was trained again using the GridSearchCV module (Cournapeau, 2007b), which iterates through all permutations of incoming parameters to find the best hyperparameter. Third, contributions of the evaluation factors were obtained and ranked. Finally, evaluation factors were filtered. Through analyzing the initial factors and their contribution in Table 3, we determined that the contribution of six factors, namely railway, highway density, the density of urban primary roads, the density of urban secondary roads, the density of urban tertiary roads, the density of urban quaternary roads, were less than 0.01%, and they apparently did not in the same order of magnitude as other factors. Therefore, we considered these six factors as noisy data and removed them. From further evaluation. The remaining evaluation factors are shown in Fig. 5.

(1) Data annotation. The debris flow inventory was derived from the nationwide geohazard census done by the Resource and Environment Science and Data Center (www.resdc.cn/). After data verification with local histories, books, reports, statistic data, relevant field surveys, and related literature, a total of 274 debris flow hazards in the study area were found as of the end of 2019, and the attribute information included the field number, geographic location, damage, groundwater grade, and the current degree of stability, and so on. We used this inventory to annotate each watershed unit. Those units that had experienced debris flows were assigned a label of ‘1’, while those that had not were assigned a label of ‘0’. The entire set of watershed units was subsequently divided into two groups, namely ‘debris flow’ and ‘no debris flow’.

(2) Data sampling. To prevent sample imbalance from affecting model accuracy, we used the synthetic minority oversampling technique (SMOTE) to balance the sample size, which analyzed a small number of data and added simulated new data to the dataset when needed (Wu, Yang & Niu, 2020; Zhu & Zhang, 2022). Eventually, the watershed unit ratio with ‘debris flow’ to those with ‘no debris flow’ was 1:1, and the total sample size was 2,140.

(3) Datasets division. The whole dataset was divided into two subsets of 70% to 30% for DFS model training and testing, respectively (Huang et al., 2022).

Analysis of the pre-processing results

The AUC of the RF model using test data increased from 0.73 to 0.97. The ROC curve was closer to the upper left corner (Fig. 6). Further, model training time decreased by 23%, from 105 to 84 s.

Data preprocessing eliminated the impact of redundant data on the model and on the remaining evaluation factors, so the contribution rate of the remaining evaluation factors such as topographic conditions, human activities, and vegetation conditions increased by 0.068, 0.065, and 0.017, respectively, while those of rainfall and geological conditions decreased by 0.115 and 0.035, respectively (Fig. 7). In addition, the contribution rate of the most important evaluation indicator (relative elevation difference) increased by 0.09 (Table 3).

Analysis of evaluation factors

Experimental results indicated that topographic conditions were the decisive factors during the formation of debris flow in the YNR Basin, and geological conditions, rainfall conditions, human activities, and vegetation conditions were important factors, with contribution rate corresponding to 0.425, 0.195, 0.173, 0.133, and 0.074 (Fig. 7), respectively. The top three explanatory conditions were topographic, geological, and rainfall.

The relative elevation difference (contribution rate: 0.274) was the most important evaluation indicator with a key role in the formation of debris flows. This was followed by watershed area, NDVI, and density of county and town roads ranking 3^rd and 4^th, respectively. The others were in the order of average slope, average rainfall during the rainy season, land use, sand content, silt content, clay content, the number of heavy rainstorms, fracture zone density, and the number of heavy rains. The number of rainstorms had a lesser impact on the formation of debris flows (Fig. 8).

Model accuracy analysis

The RF model had higher values of Accuracy, Precision, Recall, F1-score, and Kappa (Table 5), and its ROC curve converged faster than that of RBFSVM and Linear SVM (Fig. 9). This indicated that RF was most suitable of the three models for examination of the spatial correlation between historical debris flows and elevation factors, improving the assessment accuracy of DFS.

Susceptibility analysis

Susceptibility class and density distribution of debris flows were positively correlated with a higher density of debris flows leading to higher susceptibility classes. The difference in density between very high and very low susceptibility zones serves as an indicator of the predictive performance of the model (Li et al., 2022). Our results showed that the prediction performance of RF was better than that of RBFSVM, while that of linear SVM was the lowest (Fig. 10). All three models demonstrated an increase in debris flow density with increasing susceptibility class but the difference in debris flow density between the very high and very low susceptibility zones was greatest for the RF model with 47 debris flows per 1,000 km², followed by RBFSVM with 37 debris flows per 1,000 km², and Linear SVM with 11 debris flows per 1,000 km². These results indicated that the RF model was more adept at discriminating very high and high susceptibility zones and exhibited superior predictive performance when compared to the other two models.

In addition, we used prediction rate curves to visualize the model’s predictive performance. On the X-axis, we sorted the predicted susceptibility of the test dataset in descending order (i.e., from large to small susceptibility), using the method of accumulation of the corresponding watershed area; The cumulative proportion of the number of debris flow in the corresponding watershed unit was represented on the Y-axis. The farther the verification rate curve was from the diagonal from 0 to 1, the better the predictive power of the model. The larger the gradient in the first half of the curve, the stronger the prediction ability (Chung & Fabbri, 1999; Remondo et al., 2003). As can be seen from Fig. 11, the RF model had the strongest prediction ability, followed by RBFSVM, and Linear SVM was the lowest.

Based on the analysis of historical debris flows and the overlap in susceptibility zones (Fig. 12), we determined that RF and RBFSVM modelled susceptibility zones in the YNR Basin better than the linear SVM. The linear SVM predicted many of the very high and high susceptibility zones which had no historical debris flows, therefore the credibility of its predictions was low. The high dimensionality of debris flow data may be a contributing factor, as linear SVM had difficulty effectively capturing the complex nonlinear relationships in the dataset. This resulted in tendency to overfit the training data, reducing its accuracy in predicting susceptibility to debris flow hazards. The Kappa value of the linear SVM model was only 0.37 (Table 5), indicating that the model was unable to make consistent predictions at different conditions or times, indirectly confirming the existence of the overfitting problem. This instability limits the applicability of the model to DFS prediction.

The DFS spatial distribution map obtained with RF showed that the very high susceptibility zones in the upstream section of the YNR Basin were mainly distributed along the mainstream of the Nujiang River. The dominant factors were topographic and geological conditions in the very high susceptibility zones of Gongshan, Fugong, and Lushui counties associated with active neotectonic movements resulting in the development of numerous bulge structures and compressional folds. The swiftly flowing water, driven by the effects of tectonic activity and fluvial erosion, results in the formation of steep riverbanks. After the Holocene, the bulge of the mountains and the deepening of the river valleys led to the formation of alpine-valley landscapes. Furthermore, the deep and large fractures of the Nujiang River, along with numerous tectonic fractures and joint fissures, control the geomorphological development in the area, leading to extreme gravity geomorphology, including debris flows (Liu & Tang, 1995; Tang & Zhu, 2003).

The RF model performed more robustly for susceptibility zoning in the YNR Basin (Fig. 13), as shown by the overlap between historical debris flows and susceptibility areas and by the model accuracy. The modeling results of RBFSVM showed an overall higher susceptibility than that indicated by historical data. In the Linear SVM, many historical debris flow areas were distributed in very low and low susceptibility zones and few historical debris flow areas were distributed in some high susceptibility zones, which rendered predictions less reliable, lead to overfitting of the training data, and resulted in lower accuracy in predicting the susceptibility to debris flow disasters.

Consequently, the RF model was used for susceptibility zoning in the downstream section of the YNR Basin. The results showed that areas with very high and high susceptibility to debris flows were concentrated in northern Longling County, northern Longyang District, eastern and southeastern Zhenkang County, eastern Shidian County, and central Yongde County, and those of susceptibility probability response to 0.86, 0.85, 0.84, and 0.81, respectively.