Prediction of Safety Factor of Newly Failed Red Clay Slope Based on XGBoost-PSO-SVR Model

1. Introduction

Among various geological disasters, landslides are characterized by their wide distribution, high frequency of occurrence, and significant damage. Extensive research has been conducted on slope stability analysis, yielding substantial results. However, most studies focus on slopes before sliding occurs. Even though a few scholars have explored the stability of already failed slopes, their research primarily focus on ancient landslides. Such as Yang et al., Hu et al., Zhu et al. and Frink et al. studied the stability of ancient landslides by field survey and monitoring, satellite remote sensing, optical remote sensing dynamic monitoring, drone aerial survey and numerical simulation [1-8]. While these studies have contributed to understanding ancient landslide revival, the shear strength of the sliding zone in newly failed slopes has not been restored, and surface cracks remain unblocked. Consequently, findings from ancient landslide studies cannot be directly applied to newly failed slopes. The stability of the newly failed slope at Xintang Gao Kanzi was analyzed by using the transfer coefficient method [9]. However, only the weight of the slope and the weight under a once-in-50-years heavy rain scenario were considered. The parameters considered were also limited to the natural unit weight, internal friction angle, and cohesion of both the slope and the sliding zone. The hydraulic coupling numerical simulation was used to simulate its stability by Wang et al [10]. But only three kinds of rainfall conditions, such as light rain, moderate rain and heavy rainstorm were considered, without considering the influence of crack state. It is suggested that abundant loose material sources and dominant joint structures could provide fundamental conditions for the transition from shallow to deep sliding of the slope [11]. However, only natural conditions and heavy rainfall scenarios were considered, which was insufficient. The reverse analysis method and cloud model method were used to evaluate the stability of a newly failed slope [12]. But the influence degree of crack length, width, depth and position on the stability of the slope was not considered.

Although the aforementioned analytical methods have their merits, scholars have rarely considered both rainfall conditions and crack states simultaneously, despite their interaction and mutual constraints [13-18]. At present, machine learning has been widely used to predict in various research fields. Models such as Decision Tree (DT) [19], K-Nearest Neighbors (KNN) [20], Naive Bayes (NB) [21], Random Forest (RF) [22], eXtreme Gradient Boosting (XGBoost) [23], and Support Vector Regression (SVR) [24] each have their advantages and disadvantages. For example, the XGBoost model has the advantages of non-linear data processing, low computational load, faster operation speed, and better prevention of overfitting [25-26]. However, it has the disadvantages of easy overfitting and sensitivity to outliers. The SVR model can effectively solve practical issues such as small sample sizes, nonlinearity, high dimensionality, and local minima, demonstrating excellent generalization performance [27-28]. However, it suffers from slow training and difficult parameter selection. It was believed that combining multiple single machine learning models can yield a predictive model with superior performance [29-31].

The study selected 11 parameters commonly used for red clay slope evaluation, including slope height, slope angle, cohesion, internal friction angle, rainfall conditions, and crack status. The safety factor of the slope under different parameter values was simulated and analyzed using GeoStudio software. Subsequently, an adaptive weighted XGBoost-PSO-SVR hybrid model was trained with the simulation results to establish a prediction model. The model's effectiveness of this model was verified by comparing its prediction results with those of single machine learning models like XGBoost, PSO-SVR, and DT. Finally, the model's accuracy was further validated through a case study of a recently failed slope in Yongchun County, Fujian Province. This study provided a new approach for the stability prediction of recently failed slopes under comprehensive consideration of various factors, such as rainfall conditions and crack status.

2. Materials and Methods

2.1Acquisition of Research Data

To obtain the safety factor of recently failed red clay slopes under rainfall conditions, a finite element model was established using GeoStudio software for numerical simulation. The slope soil assumed to be a single layer of red clay. Since recently failed slopes have numerous and scattered cracks, it would be a huge task to represent all the cracks in the model. Therefore, the cracks were simplified into three main types: the main crack in the slip zone, the crack at the top of the slope, and the crack on the slope surface, respectively. The length and width of the main crack in the slip zone could be equivalently processed using Equation 1 and 2. The length of the cracks at the top and surface of slope was calculated as half of the main crack length, and their width could be equivalently processed using Equation 3.

$L_m={\displaystyle \frac{{\displaystyle \sum _{i=1}^n{l}_i{x}_i}}{{\displaystyle \sum _{i=1}^n{x}_i}}}$ （1）

$D_m={\displaystyle \frac{{\displaystyle \sum _{i=1}^n{d}_i{x}_i}}{{\displaystyle \sum _{i=1}^n{x}_i}}}$ （2）

$D_n={\displaystyle \frac{{\displaystyle \sum _{i=1}^n{S}_i}}{S}}\times L$ （3）

Where L_m is the equivalent length of the main crack; l_i is the length of each crack in the slip zone; x_i is the horizontal distance of each crack from the center of the slip zone; D_m​ is the equivalent width of the main crack; d_i is the width of each crack. D_n​ is the equivalent width of the cracks at the top of the slope or on the surface of slope; S_i is the area of each crack at the top of the slope; S is the total area of the top of the slope; L is the length of the top of the slope.

Since the study focused on recently failed slopes, a slip zone was assumed to exist. Therefore, the most dangerous slip surface of the single-layer homogeneous soil in the model was set as the slip zone. However, whether the final slip surface passed through this slip zone was determined by the numerical simulation of GeoStudio software, without forcing the slip surface to pass through the slip zone in this simulation. Taking a model with slope height of 10 meters and slope angle of 45° as an example, the final model was shown in Figure 1. In the figure, N1 represented the red clay layer base; N2 was the recently failed red clay layer; H1was the slip zone; L1 was the main crack in the slip zone; L2 was the crack at the top of the slope; L3 was the crack on the surface of slope; W1 was the groundwater level. The constitutive model of this model was selected as the Mohr-Coulomb model. Since the rainfall intensity in this simulation was relatively high, for the good infiltration channel cracks and slip zone soil, the rainfall intensity was set as a flow boundary. But for the red clay layer without cracks, the rainfall boundary condition was set as a zero-pressure water head boundary. The seepage in this model was transient seepage, and the influence of groundwater was not considered. Therefore, the groundwater level was set as close to the bottom of the slope as possible, as shown in the position of W1 in Figure 1.

The parameters selected for simulation included slope height (H), slope angle (β), cohesion (c), internal friction angle (φ), unit weight (γ), rainfall intensity (I_r), rainfall duration (T_r), main crack width (D_m), main crack depth (L_m), crack area ratio at the top of the slope (S_t), and crack area ratio on the surface of slope (S_f). According to the Engineering Geology Manual, the minimum unit weight of red clay is 16.5 kN/m³, and the maximum is 18.5 kN/m³. The benchmark values for the three groups were selected as 17.0 kN/m³, 17.5 kN/m³, and 18.0 kN/m³ using the equal division method. The grouping benchmark values for internal friction angle and unit weight were selected in the same way. The meteorological department defines rainfall less than 10 mm in 24 hours as light rain, between 10 mm and 25 mm as moderate rain, and between 25 mm and 50 mm as heavy rain. Therefore, 10 mm, 25 mm, and 50 mm were selected as the benchmark values for each group. For parameters such as slope height, slope angle, and rainfall duration, the benchmark values were selected based on the common classification standards used by scholars. The final benchmark values for each group of parameters were shown in Table 1. The safety factor of the slope was then calculated for different combinations of conditions in each group, with each parameter varying by ±5%, ±10%, ±15%, ±20%, and ±25% of its group benchmark value. For example, the combinations of conditions for Group II were shown in Table 2.

The density, cohesion, and internal friction angle of the cracks in this numerical simulation were all set to zero. The density of the slip zone soil was taken to be consistent with the density of the red clay being simulated. The cohesion and internal friction angle of the slip zone soil were referenced from the research results of Tang [32] and Ren [33], and were set at 19.5 kPa and 10.73°, respectively. The slope model was set as unsaturated, and the sample material was selected from the built-in clay material of GeoStudio software. The saturated and residual water content of the soil were both set at 45% and 10%, respectively. The permeability coefficients of red clay, slip zone soil, and cracks were taken as 5×10⁻¹⁰ m/s, 5×10⁻⁶ m/s, and 1 m/s, respectively. The relationship curves of matric suction with volumetric water content and water X-conductivity of the slip zone soil were shown in Figure 2 and Figure 3, respectively.

According to the above simulation scheme and parameter values, numerical simulations were conducted using GeoStudio software. Since the safety factor of the slope dynamically changes during rainfall, the safety factor obtained in this simulation was the one at the last moment of the rainfall duration. A total of 363 sets of simulation results were obtained, and a part of the results from Group II were shown in Table 3.

2.2. XGBoost

XGBoost is an algorithm based on decision trees, with decision trees being its fundamental components. During the decision tree process, subsequent trees are trained based on the residuals of the previous tree. Through continuous iterative optimization, the residuals are minimized, ultimately enhancing the overall model's prediction accuracy. The objective function of the XGBoost model consists of a loss function and a regularization term, calculated according to Equation 4 [34].

$O={\displaystyle \sum _{i=1}^{n}l\left(y_i,\stackrel{\wedge }{y_i}\right)}+{\displaystyle \sum _{k=1}^k\Omega \left(f_k\right)}$ （4）

where O is the objective function; y_i​ is the measured value of the i target; $\stackrel{\wedge }{y_i}$is the predicted value of the i target; l(y_i, $\stackrel{\wedge }{y_i}$) is the difference between y_i and $\stackrel{\wedge }{y_i}$​; n is the number of samples; Ω(f_k) is the complexity of the tree model for the k sample feature parameter f_k​; k is the number of sample feature parameters.

The objective function was approximated by performing a second-order Taylor expansion on it, thus transforming Equation 4 into Equation 5.

$O^{(t)}\approx {\displaystyle \sum _{i=1}^n\left[l\left(y_i,y_i^{t-1}\right)+g_i{f}_t(x_i)+0.5h_i{f}_t^2(x_i)\right]}+\Omega \left(f_t\right)+C$ （5）

where g_i and h_i are the first and second derivatives of $l\left(y_i,y_i^{t-1}\right)$, respectively.

Since the goal of the model was to minimize the objective function, the constant term was temporarily disregarded. After removing the constant term $l\left(y_i,y_i^{t-1}\right)$and C from Equation 5 and summing the objective function in the form of leaf nodes, Equation 6 was obtained [35].

$O^{(t)}={\displaystyle \sum _{j=1}^T\left[w_j{\displaystyle \sum _{i\in I}{g}_i}+0.5({\displaystyle \sum _{i\in I}{h}_i}+\lambda )w_j^2\right]}+\gamma T$ （6）

where I is the set of samples on each leaf; w_j​ is the output score of each tree leaf node; T is the number of leaf nodes of the split tree; λ and γ are weight factors, controlling the weights of the corresponding parts.

2.3. SVR

SVR is a small-sample creative machine learning method based on statistical learning theory, aiming to minimize the model's structural risk. When dealing with nonlinear problems, this learning method maps the original data x to a high-dimensional feature space to obtain φ(ｘ), thereby transforming it into a linear problem for solution. It has strong generalization performance and is effective for regression problems. Suppose there exists a training set {(x_i,y_i)|i=1,2,3, ∙∙∙,n}, where x_i​ is the input vector, y_i​ is the output target, and n is the number of samples. The input vector and output target can be described by Equation 7 [36].

$f(x)=w^T\phi (x)+b$ （7）

where f(x)is the predicted value; w^T is the weight vector; φ(x) is the mapping of the input variable x in the high-dimensional feature space; b is the threshold.

After a series of transformations and the introduction of the Lagrange function and kernel function, the objective function and kernel function of SVR were shown in Equations 8 and 9 [37].

$f(x)={\displaystyle \sum _{i=1}^n\left({\alpha }_i-{\alpha }_i^{*}\right)}K(x_i,x_j)+b$ （8）

$K(x_i,x_j)=exp\left(-{\displaystyle \frac{\parallel x_i-x_j{\parallel }_2^2}{2{\sigma }^2}}\right)=exp\left(-g\parallel x_i-x_j{\parallel }_2^2\right)$ （9）

where α_i and α_i^*​ are Lagrange multipliers; K(x_i,x_j) is the kernel function; ∥x_i​−x_j​∥² is the squared Euclidean distance between two feature vectors; σ is the width of the kernel function; g is the parameter of the kernel function.

C is the penalty coefficient of SVR, mainly used to control the error range of the model to avoid underfitting or overfitting; g is the kernel parameter, mainly used to control the distribution of data in the new feature space, determining the number of support vectors, and thereby affecting the speed of training and prediction. Therefore, determining the optimal parameters is a crucial part of the SVR algorithm [38].

2.4. PSO Algorithm

The PSO algorithm, proposed by Kennedy [39] and Eberhart [40], is a search algorithm inspired by the foraging behavior of birds. It is characterized by its high efficiency and fast search speed. The algorithm mainly iteratively calculates the initial position and velocity of a group of random particles to find the optimal solution [41]. Before the algorithm runs, a group of particles with vector dimension n is initialized. The position of a particle can be denoted as a point in an n-dimensional search space, with its coordinates represented as x_i=(x_i,1,x_i,2,∙∙∙,x_i,D), which is also considered a solution in the n-dimensional optimization space. Its flight velocity is denoted as v_i=(v_i,1,v_i,2,∙∙∙,v_i,D); the historical optimal coordinates of the i particle are P_i=(P_i,1,P_i,2,∙∙∙,P_i,D); and the optimal coordinates experienced by each particle are P_g=(P_g,1,P_g,2,∙∙∙,P_g,D). During the flight process, the particle swarm can iteratively calculate, as shown in Equations 10 and 11 [42].

$v_{id}^{k+1}=\omega v_{id}^k+c_1{r}_1(p_{best}-x_{id}^k)+c_2{r}_2(g_{best}-x_{id}^k)$ （10）

$x_{id}^{k+1}=x_{id}^k+v_{id}^{k+1},i=1,2,\cdots ,m;d=1,2,\cdots D$ （11）

where m is the size of the particle swarm; D is the dimension of the particle swarm; v^k_id​ is the velocity; x^k_id​ is the position; k is the iteration number; c₁​ and c₂​ are acceleration factors, controlling the state of particles maintaining p_best​ and g_best​; r₁ and r₂ are random numbers between [0,1]. ω is the inertia weight, used to control the influence of the original speed on the new speed. When ω is large, the algorithm has strong global search capability, and vice versa, it has strong local search capability, which can be expressed by Equation 12.

$\omega ={\omega }_{\max }-{\displaystyle \frac{k\left({\omega }_{\max }-{\omega }_{\min }\right)}{k_{\max }}}$ （12）

where k is the current iteration step; k_max​ is the maximum iteration step; ω_max and ω_min​ are the maximum and minimum values of ω, respectively.

2.5. Adaptive Weighting Combination Model

The adaptive weighting combination model is an improvement based on the residual weighting method. Its main approach is to assign weights to the current sample model based on the average weight of the previous m samples [19]. The optimal m needs to be determined through trial calculations, which can be done using Equations 13 to 15.

$\stackrel{\wedge }{y_i}={\displaystyle \sum _{j=1}^n\left[{\omega }_j(i)\stackrel{\wedge }{y_i^j}\right],i\ge 2}$ （13）

${\omega }_j(i)={\displaystyle \frac{\frac{1}{{\overline{\epsilon }}_j(i)}}{{\displaystyle \sum _{j=1}^n\frac{1}{{\overline{\epsilon }}_j(i)}}}}$ （14）

${{\omega }^{\prime }}_j(i)={\displaystyle \frac{1}{m}}{\displaystyle \sum _{k=1}^m{\omega }_j(j-k)}$ （15）

where $\stackrel{\wedge }{y_i}$ is the predicted value of the i sample of the combination model; $\stackrel{\wedge }{y_i^j}$ is the predicted value of the i sample of the j model; ω_j(i)​ is the residual weighting combination model weight of the j model for the i sample; ${\overline{\epsilon }}_j(i)$ is the sum of squared prediction errors of the j model for the k sample; ${{\omega }^{\prime }}_j(i)$ is the adaptive weighting combination model weight.

3. Results and Discussion

3.1. Comparison of Prediction Results of PSO-SVR, XGBoost, and XGBoost-PSO-SVR

The 363 simulation results were input into the PSO-SVR and XGBoost models for training. Both models had 290 training samples and 73 testing samples. The kernel functions of SVR mainly include linear kernel, polynomial kernel, Gaussian kernel, and Sigmoid kernel. The Gaussian kernel was chosen for SVR due to its ability to handle complex non-linear problems. The PSO algorithm was used to optimize the penalty coefficient C and parameter g with 5-fold cross-validation. According to the experience of scholars, the particle swarm size N in the PSO algorithm was set to 50, the inertia weight ω to 1.2, and the learning factors c₁​ and c₂​ to 2, with the maximum iteration number G_kset to 60. The parameter settings for the XGBoost model were shown in Table 4. Based on the above parameter settings, the prediction results for PSO-SVR and XGBoost were obtained. The adaptive weighting combination model was then used to combine these two results, yielding the prediction results for XGBoost-PSO-SVR. The results were shown in Figure 4.

From Figure 4, it could be intuitively seen that when trained using the XGBoost method, two samples had significant deviations between predicted and actual values, while the predicted values of other samples were quite close to the actual values. The mean squared error (MSE) of the test set for this method was 0.0056979, with a corresponding R² of 0.98378, indicating that the model trained using the XGBoost method could explain 98.4% of the variance in the dependent variable. This demonstrated that the training effect of this model was good. From Figure 5, it could be seen that when trained using the PSO-SVR method, three samples had significant deviations between predicted and actual values. The corresponding MSE and R² were 0.0037367 and 0.98515, respectively, indicating that this method reduced the MSE by 34.4% compared to the XGBoost method, resulting in better training performance. From Figure 6, it can be seen that when using the XGBoost-PSO-SVR combined algorithm, only one sample had a significant deviation between predicted and actual values, and the fit of the predicted values to the actual values was higher than the above two methods. The corresponding MSE and R² were 0.0016001 and 0.9919, respectively. It representing a 71.9% reduction in MSE and a 0.83% increase in the accuracy of the dependent variable explanation compared to the XGBoost method; and a 57.8% reduction in MSE and a 0.69% increase in the accuracy of the dependent variable explanation compared to the PSO-SVR method. Therefore, the XGBoost-PSO-SVR combined algorithm could further reduce the model's MSE and improve accuracy, resulting in the best fit. In MATLAB software, the tic and toc functions were used to calculate the start and end times of the algorithm, respectively, to determine the algorithm's execution time. The calculated execution times for PSO-SVR, XGBoost, and XGBoost-PSO-SVR were 8.14 s, 6.21 s, and 13.32 s, respectively. It can be seen that XGBoost has the fastest computation time, while XGBoost-PSO-SVR has the slowest. Since the XGBoost-PSO-SVR algorithm includes the running time of the XGBoost algorithm, the optimization time of the PSO algorithm, and the prediction time of SVR, its running time was longer than that of the individual XGBoost or PSO-SVR algorithms. However, since the data volume analyzed in this case was relatively small and the running times were all relatively short, they were within an acceptable range. If the data volume to be analyzed is larger, the appropriate algorithm should be selected by considering factors such as accuracy and running time.

3.2. Comparison of Prediction Results of XGBoost-PSO-SVR with Other Machine Learning Models

To further verify the prediction accuracy of the XGBoost-PSO-SVR combination model, four single machine learning models, namely DT, KNN, NB, and RF, were used to predict the safety factor. The comparison of the predicted safety factor values with the actual values obtained from these four machine learning models was shown in Figure 7. The MSE and R² of the safety factor prediction results of the XGBoost-PSO-SVR combination model and the four single machine learning models were shown in Table 5.

From Figure 7 and Table 5, it could be seen that the XGBoost-PSO-SVR model had the smallest MSE, followed by the DT model, and the KNN model had the largest; the XGBoost-PSO-SVR model had the largest R², followed by the DT model, and the KNN model had the smallest. Therefore, the training effects of the five models, from highest to lowest, were: XGBoost-PSO-SVR, DT, NB, RF, and KNN. Considering the characteristics of each machine learning model, the main reasons were as follows: The DT model usually assumed independence between attributes during construction, but in reality, the factors affecting the slope safety factor were intercoupled; the performance of the KNN model was easily affected when the number of samples of different categories varies greatly, and the samples in this case include parameters from three different groups; the NB model also assumed independence between attributes; the RF model might be affected by the majority class samples, leading to a decrease in the prediction performance of the model for minority class samples, which might occur when the model randomly assigned training and testing samples. SVR was good at handling high-dimensional data and could effectively solve nonlinear classification problems through kernel trick techniques, and the sample dimension in this case was 11, which was suitable for this model. Meanwhile, XGBoost constructed a strong learner by integrating multiple decision trees, and this integration method allowed XGBoost to significantly improve prediction accuracy. Therefore, the training effect of the XGBoost-PSO-SVR model on this sample was the best.

4. Justification

Taking a recently failed slope in Lengshuicun, Yongchun County in China, as reported in [12], the XGBoost-PSO-SVR model established in this paper was used for comparison. Since the opening width at the rear edge of the slide mass of this slope was approximately 0.3 to 1 m, the main crack width was taken as 1.0 m. Also, since the rear edge of the slope body has already moved down as a whole by about 2 to 2.5 m, the main crack depth was taken as 2.5 m. Since there was rainfall from November 4 to 7, 2016, with a total rainfall of 126.1 mm, the rainfall duration was taken as 4 days, and the rainfall intensity was taken as 31.525 mm/day. According to the field investigation, as shown in Figure 8 and Figure 9, the crack area ratio at the top of the slope was estimated to be 20%, and the crack area ratio on the surface of slope was taken as 10%. The values of the other parameters were shown in Table 6.

The parameters listed in Table 6 were input into the XGBoost-PSO-SVR model for prediction, and the results were shown in Figure 10.

From Figure 10, it could be seen that the model predicted the 74th sample for this slope, with a value of 0.966. According to the Technical Code for Building Slope Engineering of Chia (GB50330-2013), the stability state of this slope was unstable. According to the reference [12], the soil at the front edge of the slope moved forward by about 0.5 m from November 4 to 7, verifying the accuracy of the XGBoost-PSO-SVR model. The safety factors obtained by using XGBoost, PSO-SVR, DT, KNN, NB, and RF for this slope were shown in Table 7.

The numerical analysis of this slope using GeoStudio software yielded a safety factor of 0.988. From Table 7, it could be seen that the deviation between the predicted value obtained by XGBoost-PSO-SVR and the numerical analysis value was 0.022. Although this deviation was larger than that of XGBoost and DT, the prediction values of XGBoost and DT were both greater than 1. According to the Technical Code for Building Slope Engineering of Chia, the stability state of this slope should be judged as under stable, which did not match the actual situation of the slope sliding again. Therefore, in general, compared with other methods, the prediction of XGBoost-PSO-SVR was closer to the actual situation.

5. Conclusions

(1) By comparing the MSE and R2 indicators of the models, the XGBoost-PSO-SVR combined algorithm reduced the mean squared error by 71.9% and 57.8% compared to the single XGBoost and PSO-SVR methods, respectively, and increased the accuracy of the dependent variable explanation by 0.83% and 0.69%, respectively. This model had significant advantages in training accuracy and fitting effect. Moreover, compared with the four single models DT, NB, RF, and KNN, the XGBoost-PSO-SVR combined algorithm also achieved the best training effect.

(2) Using the XGBoost-PSO-SVR combined algorithm to predict the stability of a recently failed slope in Lengshuicun, Yongchun County, from November 4 to 7, 2016, the predicted safety factor was 0.966, indicating that the slope was in an unstable state. This prediction result was consistent with the actual situation.

(3) Due to the limitation of workload, this study simplified the cracks on the slope in the numerical simulation, considering only the main crack in the slip zone, the crack at the top of the slope, and the crack on the surface of slope. At the same time, the influence of groundwater and the spatial variability of slope parameters were not considered, which may cause some differences from the actual situation. Future research will focus on more reasonably considering the impact of the above factors on the stability of recently failed slopes.