Predicting Factor of Safety of Slope Using an Improved Support Vector Machine Regression Model

Abstract

To reduce the disasters caused by slope instability, this paper proposes a new machine learning (ML) model for slope stability prediction. This improved SVR model uses support vector machine regression (SVR) as the basic prediction tool and the grid search method with 5-fold cross-validation to optimize the hyperparameters to improve the prediction performance. Six features, namely, unit weight, cohesion, friction angle, slope angle, slope height, and pore pressure ratio, were taken as the input of the model, and the factor of safety was taken as the model output. Four statistical indicators, namely, the coefficient of determination (R²), mean absolute percentage error (MAPE), mean absolute error (MAE), and root mean squared error (RMSE), were introduced to assess the generalization performance of the model. Finally, the feature importance score of the features was clarified by calculating the importance of the six features and visualizing them. The results show that the model can well describe the nonlinear relationship between features and the factor of safety. The R², MAPE, MAE, and RMSE of the testing dataset were 0.901, 7.41%, 0.082, and 0.133, respectively. Compared with other ML models, the improved SVR model had a better effect. The most sensitive feature was unit weight.

1. Introduction

The mechanical balance of slopes is often destroyed due to natural and human processes, inducing landslides and other geological disasters and further causing great economic loss and significant casualties [1,2,3]. As shown in Figure 1, slope is often encountered in geotechnical engineering practice, which is of great significance for safe construction and operation, scientific design, and long-life operation of geotechnical engineering. Landslides are the main type of geological hazard in most mountainous areas, and accurate prediction of slope stability allows for the timely relocation of people and property to minimize losses [4,5,6]. Therefore, rapid and accurate slope stability prediction is always one of the most important research topics in the field of geological engineering [7,8,9].

Since Terzaghi [10] published the classic work entitled “Mechanism of Landslides” in 1950, numerous studies have focused on analyzing slope stability using various methods such as theoretical, analytical, experimental, numerical, and statistical approaches. The stability of the slope is generally stated by the term “factor of safety” (FS), and a slope is considered safe when its FS is greater than 1 [11]. In general, traditional approaches to solving FS can be categorized into two main groups: analytical methods and numerical modeling methods. The analytical methods are mainly based on the slope displacement model and are used to identify the most possible sliding surface. Among them, the limit equilibrium method and the circular/non-circular failure surface method are the most common. Wang et al. [12] improved the traditional slope stability analysis method (i.e., the simple wedge failure mode) and used the slope stability analysis of the channel bank based on the circular failure surface method to simulate the lateral enlargement in dam breaches. Utilizing the limit equilibrium method as a foundation, Faramarzi et al. [13] analyzed the stability of the Zhanshiba power station pit in Kohgiluyeh Va Boyer-Ahmad Province, southwest of Iran. Duncan et al. [14,15] provided an admirable summary of the limit equilibrium method and its importance in understanding the various applications of geotechnical engineering. However, owing to its inherent characteristics, such as simplification and the requirement for input features across the study area, the analytical methods fall short in providing a comprehensive understanding of slope behavior [16]. For example, the rigid body limit equilibrium method (RBLEM) cannot determine the overall deformation and stress field distribution of the structure, particularly in the case of sliding surfaces [13]. Therefore, Song et al. [17] pointed out that these analytical methods are only suitable for slopes with simple geometry.

Compared with analytical methods, numerical modeling methods are a more realistic and strict method in theory [18]. At present, a large number of numerical methods, such as the discrete element method (DEM), particle flow code (PFC), boundary element method (BEM), finite element method (FEM), and numerical manifold method (NMM), have been used for factor of safety, displacement, and failure behavior analysis [19,20,21]. Duncan was the first to apply the FEM to study slope stability, ushering in a new era of numerical modeling for slope stability problems. Subsequently, Tsiampousi et al. [22] used FEM to investigate the effects of pore water pressures on the stability and serviceability of a typical cut slope in England. Wu and Hsieh [23] adopted the 3D discrete element method to analyze the seismically induced Chiu-fen-erh-shan landslide in central Taiwan. More recently, Chand and Koner [24] combined an unmanned aerial vehicle with the numerical modeling method to create a 3D numerical model in FLAC3D to evaluate dump slope stability under static and dynamic loads. The effectiveness of the numerical method depends on the computational model used for modeling and the input parameters. The input parameters require back analysis using a large number of in situ measurements, which is limited in some cases [25,26].

Scholars have developed a number of models to assess landslide susceptibility, which can generally be divided into four categories: heuristic models, physically based models, statistically based models, and machine learning models [27,28,29]. Reference [30] provides a valuable introduction to these models. Among them, machine learning (ML) methodologies have garnered increasing attention and adoption within the domain of civil engineering due to their high performance in characterizing the nonlinear relationships between high-dimensional data [31,32,33,34,35,36,37,38,39,40,41,42]. In the case of specific field or laboratory parameters, ML-based methods can be an alternative to traditional models [43]. Feng [44] pioneered the definition of “intelligent rock mechanics” and developed an intelligent slope stability prediction model using neural networks. Zhao et al. [45] utilized the relevance vector machine (RVM) to investigate the correlation between slope stability and various influencing parameters, offering valuable insights into the field. Qi and Tang [25] proposed and compared the generalization ability of six integrated ML models in slope stability prediction: the neural network, decision tree, support vector machine (SVM), gradient boosting machine, random forest, and logistic regression. Manouchehrian et al. [46] collected data from the published literature to build a database of 121 cases. Then, the regression model for slope stability analysis was found by means of a genetic algorithm, and the corresponding FS prediction model of circular mode failure was established. For slope units, Chang et al. [47] proposed a new multi-scale segmentation (MSS) method using the random forest (RF) algorithm and the multi-layer perceptron (MLP) algorithm to develop different ML models. Recently, Bardhan and Samui [48] gave a review of the application of machine learning models in slope stability analysis, discussing the latest applications of various ML models to solve slope stability problems. On the basis of the literature review,

Table 1. Previous works about slope stability using ML-based methods.

Source	Technique	Input	Number of Data
Lu and Rosenbaum [53]	Back propagation Neural network (BPNN)	γ, c, φ, α, H, ru	32
Huang et al. [54]	Convolutional neural network (CNN)	γ, c, φ, α, H, ru	64
Sakellariou and Ferentinou [55]	BPNN	γ, c, φ, α, H, ru	46
Wang et al. [56]	BPNN	γ, c, φ, α, H	27
Samui [57]	SVM	γ, c, φ, α, H, ru	46
Zhao [58]	SVM	γ, c, φ	10
Choobbasti et al. [59]	Artificial neural network (ANN)	γ, c, φ, α, H, ru	36
Zhang and Luo [60]	K-nearest neighbor (KNN)	γ, c, φ, α, H, ru	39
Das et al. [61]	ANN	γ, c, φ, α, H, ru	46
Zhao et al. [45]	RVM	γ, c, φ, α, H, ru	80
Liu et al. [62]	Extreme learning machine (ELM)	γ, c, φ, α, H, ru	97
Hoang and Pham [63]	Least squares support vector classification (LS-SVC)	γ, c, φ, α, H, ru	168
Suman et al. [64]	Functional networks (FNs), multivariate adaptive regression splines (MARS), and multigene genetic programming (MGGP)	γ, c, φ, α, H, ru	103
Verma et al. [65]	ANN	c, φ, α, Pore pressure	100
Rukhaiyar et al. [11]	Particle swarm optimization (PSO)–ANN	γ, c, φ, α, H, ru	83
Xue [66]	PSO–LSSVM	γ, c, φ, α, H, ru	46
Feng et al. [67]	NBC naive bayes classifier (NBC)	γ, c, φ, α, H, ru	82
Qi and Tang [25]	Six integrated AI approaches	γ, c, φ, α, H, ru	148
Zhou et al. [68]	Gradient boosting machine (GBM)	γ, c, φ, α, H, ru	221
Niu et al. [69]	Genetic algorithm (GA)–SVM	γ, c, φ, α, H, ru	40
Manouchehrian et al. [46]	GA	γ, c, φ, α, H, ru	121
Abdollahi et al. [70]	PyCaret AutoML environment	/	/

Note: unit weight γ, cohesion c, friction angle φ, slope angle α, slope height H, and pore pressure ratio r_u.

summarizes previously published ML-based methods on slope stability. Despite many efforts in ML modeling, there are still some issues that need to be further investigated. First, the algorithm’s hyperparameters are not determined by standard methods, and more efficient hyperparameters may be ignored [49]. Second, the existing studies have not been cross-validated, and the accuracy of the predicted results may not be scientific [50]. Last but not least, most studies validate the applicability of their algorithms with only a few cases, which can lead to unreliable conclusions [51]. A predictive model that performs well on a site-specific landslide may perform poorly in others [52].

To address the above problems, this study proposes a new ML model for predicting FS. Taking support vector machine regression (SVR) as the main prediction tool, the hyperparameters of the model were optimized by combining the grid search method and 5-fold cross-validation. A total of 166 samples were collected to test the predictive ability of the established improved SVR model. Four statistical indicators, namely, the coefficient of determination (R²), mean absolute percentage error (MAPE), mean absolute error (MAE), and root mean squared error (RMSE), were adopted to assess the prediction performance. Subsequently, it was compared with other ML models. Finally, sensitivity analysis was conducted to ascertain the sensitivity of each input feature to FS.

2. Dataset Construction and Model Evaluation Metrics

2.1. Input and Output Features

Natural geomaterials such as rock and soil have a complex formation history, inherent heterogeneity in composition, and unique configuration in structure [71,72,73]. The heterogeneity and complexity make it difficult to truly characterize the slope. There are many features that affect the slope stability, such as geological structure and stress, properties of geotechnical materials, the role of water, and the geometric shape and surface morphology of the slope. Referring to the various influencing factors listed in Table 1, as shown in Figure 2, the most representative six features were selected to evaluate the FS. The six key features were unit weight γ, cohesion c, friction angle φ, slope angle α, slope height H, and pore pressure ratio r_u. H and α reflect the geometric features that affect the slope stability. γ, c, and φ reflect the basic physical and mechanical properties of geotechnical materials. The smaller the γ and the larger the c and φ, the better the slope stability. r_u reflects the magnitude of pore water pressure. An increase in r_u will directly lead to the reduction in slope stability.

As shown in Equation (1), the mapping f was established according to these six features to solve the FS value.

$$FS=f\left(\gamma ,c,\phi ,\alpha ,\mathrm{H},r_u\right)$$

(1)

A total of 166 samples of slope were collected from Refs. [74,75,76], and they included the above six input features and the output feature FS. The statistical information for all these features is given in

Table 2. Statistical characteristics of various features.

Feature	N Total	Unit	Mean	Median Absolute Deviation	Skewness	Kurtosis	Standard Deviation	Maximum
γ	166	kN/m3	22.1108	3.16	−0.1522	−0.3978	4.9362	31.3
c	166	kPa	26.4072	14.84	1.8900	4.4461	30.1166	150.05
φ	166	°	30.1155	3	−1.4863	2.2215	10.0662	45
α	166	°	36.7319	6	−0.6007	−0.7066	10.0275	53
H	166	m	127.9484	42	1.1404	−0.0887	150.5761	511
ru	166	-	0.2602	0.06	−0.5997	0.3467	0.1292	0.5
FS	166	-	1.2561	0.18	0.6031	0.1997	0.3256	2.05

. This table provides a thorough array of statistical information, including the mean, standard deviation, and kurtosis of the observations. Using the entire dataset to train the ML model can lead to an overfitting problem. Therefore, before dividing the entire dataset, the order of the samples was shuffled, and the data set was randomly split into a training dataset and a testing dataset at a division ratio of 4:1. The training dataset was used for training, and the testing dataset- was used for testing the ML model. It is important to note that the training dataset and testing dataset can be enriched with new data in the future to improve the generalizability and predictive accuracy of the ML model.

The distribution characteristics of the dataset were visualized using violin plots. As depicted in Figure 2, the violin plot integrates the characteristics of density estimation and box plots to simultaneously show the dataset’s maximum, minimum, mean, and distribution density. Specifically, the width of the outer contour of Figure 3 denotes the density of the distribution of the dataset. The rectangle in the center of Figure 3 is a box plot, where the upper and lower edges are the upper and lower quartiles, respectively, and the middle is the median [77]. A cursory inspection of this violin plot reveals that the dataset is evenly distributed across the input features, with no data points that deviate significantly.

Figure 4 illustrates the correlations of the six features considered and presents the corresponding statistical distribution for each feature. The interaction of the input features is represented by the correlation coefficient. As shown, the correlation coefficients between the individual features are poor (Pearson’s r < 0.7). Specifically, there is a moderate correlation between γ and H but almost no correlation between γ and r_u (Pearson’s r = 0.046). From Figure 4, it can be seen that this dataset is very widely distributed, and most of the features have asymmetric distributions. Therefore, before the actual operation, all input features need to be scaled to the range [0, 1] according to the minimum and maximum values of each input feature. This can greatly improve the operational efficiency and computation speed of the ML model.

2.2. Model Evaluation Metrics

In this study, different statistical indicators, including R², MAPE, MAE, and RMSE, were adopted to assess the prediction ability of the model. The definitions of these statistical indices are given by Equations (2)–(5) [78,79]. These statistical indicators were adopted to evaluate the relationship between the true and predicted values of FS.

Generally speaking, the R² reflects the degree of agreement between the predicted and measured values, while the statistical indicators (MAPE, MAE, and RMSE) reflect the deviation between the predicted and true values. As suggested by [80,81], the closer R² is to 1 and the closer the statistical indicators (MAPE, MAE, and RMSE) are to 0, the better the model performance.

$$MAE={\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^M\left|Y_i^{mea}-Y_i^{pre}\right|}$$

(2)

$$RMSE=\sqrt{{\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^M{\left(Y_i^{mea}-Y_i^{pre}\right)}^2}}$$

(3)

$$MAPE={\displaystyle \frac{1}{M}}{\displaystyle \sum _{i=1}^M\left|\frac{Y_i^{mea}-Y_i^{pre}}{Y_i^{mea}}\right|}\times 100\%$$

(4)

$$R^2=1-{\displaystyle \frac{{\displaystyle \sum _{i=1}^M{\left(Y_i^{mea}-Y_i^{pre}\right)}^2}}{{\displaystyle \sum _{i=1}^M{\left(Y_i^{mea}-\overline{Y^{mean}}\right)}^2}}}$$

(5)

where Y_i^mea and Y_i^pre are the ith original and predicted FS values, respectively; $\overline{Y^{mean}}$ is the average of Y_i^mea; and M represents the number of observations in the dataset.

3. Methodology

This section first introduces the basic principle of the SVR model. Then, it introduces the process of introducing the grid search method and the k-fold cross-validation method to modify the original SVR model. On this basis, an improved SVR model is established.

3.1. Support Vector Machine Regression

SVR is an extension of SVM for regression analysis and is widely used by earth and computer scientists, especially in landslide evaluation [82,83,84]. SVR is a common model of SVM to solve regression problems. This section gives a brief introduction to SVR, and the reader is referred to Smola and Schölkopf [85] for more information.

The basic principle of SVR can be simply described as follows. Let the training dataset D contain m training samples, which can be represented as shown in Equation (6):

$$D=\left\{\left(X_1,Y_1\right),\left(X_2,Y_2\right),\dots ,\left(X_m,Y_m\right)\right\}\subset R^n\times R$$

(6)

where X_i denotes the input feature vector, Y_i denotes the regression value, and n denotes the input feature dimension.

The input feature vector is mapped into the high-dimensional space R^m by nonlinear mapping. The decision function in high-dimensional space is shown as follows:

$$\begin{array}{ccc}f\left(X\right)={\omega }^T\phi \left(X\right)+b& with& \omega \in R^m,\end{array}b\in R$$

(7)

where f(X) is the optimal function, ω denotes the adjustable weight vector, b denotes the threshold value, and φ(X) represents a set of nonlinear mapping functions.

Subsequently, b and ω need to be optimized as the following objective functions:

$$\begin{array}{cc}\mathit{min}& {\displaystyle \frac{{‖\omega ‖}^2}{2}}\\ s.t.& \left\{\begin{array}{c}{Y}_m-{\omega }^T\phi \left(X_m\right)-b\le \epsilon \\ {\omega }^T\phi \left(X_m\right)+b-Y_m\le \epsilon \end{array}\right.\end{array}$$

(8)

where ε denotes the maximum error allowed by regression analysis.

Similar to SVM, the relaxation factors ξ_i and ξ_i^* are introduced to transform Equation (8) into an optimization problem as follows.

$$\begin{array}{cc}\mathit{min}i& {\displaystyle \frac{1}{2}}{‖\omega ‖}^2+C{\displaystyle \sum _i^n\left({\xi }_i+{\xi }_i^{*}\right)}\\ s.t.& \left\{\begin{array}{c}{y}_i-{\omega }^T\phi \left(x_i\right)-b\le \epsilon +{\xi }_i\\ {\omega }^T\phi \left(x_i\right)+b-y_i\le \epsilon +{\xi }_i^{*}\\ {\xi }_i,{\xi }_i^{*}\ge 0\end{array}\right.\end{array}$$

(9)

where C denotes the penalty parameter. C denotes the complexity of the model. A high value of C means that the model will try to fit as much of the data as possible, which may cause the model to overfit [86].

When solving Equation (9), we introduce the Lagrange function and convert it into the dual form:

$$\left\{\begin{array}{c}\begin{array}{cc}\mathit{max}& -{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{i=j}^m\left(a_i-a_i^{*}\right)\left(a_j-a_j^{*}\right)}}K\left(X_i,X_j\right)-{\displaystyle \sum _{i=1}^m\left(a_i+a_i^{*}\right)\epsilon }+{\displaystyle \sum _i^m\left(a_i-a_i^{*}\right)Y_i}\end{array}\\ \begin{array}{cc}s.t.& \left\{\begin{array}{c}{\displaystyle \sum _i^m\left(a_i-a_i^{*}\right)=0}\\ 0\le a_i\le C\\ 0\le a_i^{*}\le C\end{array}\right.\end{array}\end{array}\right.$$

(10)

For the solution of Equation (10), the kernel function K(X_i, X_j) can be introduced, and the SVR model under nonlinear conditions can be used. The role of K(X_i, X_j) is to represent the training dataset from the original space to a higher-dimensional space so that it becomes linearly separable data. The feature spaces mapped by different K(X_i, X_j) are highly different. Referring to Yin et al. [87], the Gaussian radial basis function with higher adaptability is selected as the kernel function in this paper.

Finally, the corresponding regression equation can be derived:

$$K\left(X_i,X_j\right)=\mathit{exp}\left(-{\displaystyle \frac{‖X_i-X_j‖}{2g^2}}\right)$$

(11)

where g denotes the kernel parameter governing the complexity of the solution within the Gaussian radial basis function.

3.2. Hyperparameter Optimization

To improve the generalization ability of the ML model, appropriate data preprocessing and model hyperparameter optimization are essential. Because the numerical attributes of the datasets have large proportional differences, such as H and r_u not being the same order of magnitude. Therefore, as shown in Equation (12), the original data are normalized and transformed into dimensionless pure numerical values before model training. This preprocessing process speeds up the training process of the ML model and improves the training efficiency of the model [88].

$$x_i^{\prime }={\displaystyle \frac{x_i-x_{min}}{x_{max}-x_{\mathit{min}}}}$$

(12)

where x_max and x_min represent the original maximum and minimum values before normalization, and x_i and x_i^′ represent the values before and after normalization, respectively.

Each ML model has its hyperparameters, and the setting of the hyperparameters greatly affects the actual effect of the algorithm. For the SVR model with Gaussian radial basis function, the performance of the model largely depends on the hyperparameters g and C [89]. In the past, scholars have often used the genetic algorithm, ant colony algorithm, particle swarm optimization algorithm, and so on to investigate optimal hyperparameter values. However, these methods often face the dilemma of easily falling into local optima and a slow convergence rate.

In this study, the grid search method is adopted to search and optimize the hyperparameters of the original SVR model, and a new improved SVR model is established. The possible values of each hyperparameter are listed one by one, and the permutation and combination of different values of different hyperparameters will eventually form a parameter space. In this parameter space, the parameters are adjusted sequentially by step size. The learner is trained with the adjusted parameters, and the best-performing hyperparameter combination (highest accuracy or minimum error) is found while traversing all grids in the parameter space. The goal of parameter optimization is to find the combination of hyperparameters that makes the ML model generalize the most, which is also a training and comparison process. Therefore, in this process, k-fold cross-validation is used to calculate evaluation indicators such as RMSE, mean square error (MSE), and R² to reflect the generalization ability of the proposed model. The k-fold cross-validation is an effective way to evaluate the robustness and generalization of a model, which can avoid overfitting and underfitting to some extent. The training dataset is divided into k subsets: one subset is set aside as a new validation set for performance evaluation, and the union of the remaining k-1 subsets is used as a new training dataset for training the model. Then, the search range and step size are adjusted reversely based on evaluation index results so as to determine the optimal values of the hyperparameters. The grid search method avoids the shortcomings of local optimality and the slow convergence of traditional optimization algorithms. Compared with other parameter optimization techniques, the grid search method has good versatility and can ensure the global optimal solution in the given parameter space is found.

A generalized operational framework of the proposed improved SVR model is illustrated in Figure 5, including dataset normalization, model building, hyperparameter optimization, and performance evaluation.

4. Results and Discussion

4.1. Results of the Hyperparameters

In this study, the grid search method was used for hyperparameter tuning, the 5-fold cross-validation method was adopted to evaluate the model performance, and the best hyperparameter combination with the best prediction performance was selected. The search range of the hyperparameters was set as (2⁻⁵,2⁵), and the step size was 2^0.2. A three-dimensional grid space was created with the possible ranges of C, g, and MSE. A grid node in Figure 6 represents a kind of hyperparameter combination, and the vertical axis denotes the corresponding MSE value. As the grid color changes from yellow to dark blue, lower MSE values are achieved. By optimizing the hyperparameters C and g, the optimal solution of the hyperparameters in this search range is obtained, and the improved SVR model also achieves the best overall performance.

4.2. Results of the Improved SVR Model

Applying the proposed improved SVR model to the dataset in Section 2, the comparison between the actual value of FS and the model prediction results is shown in Figure 7. In order to better highlight the generalization ability of the improved SVR model, the original SVR model was used as a comparison for the same training and testing datasets.

To provide a clearer reflection of the prediction results of each ML model, the calculated results of the statistical indicators of the two ML models are shown in

Table 3. Performance indices of the improved SVR and SVR models.

Model	Data	R2	MAE	RMSE	MAPE
Improved SVR	Training dataset	0.9159	0.043	0.0888	3.793%
	Testing dataset	0.9010	0.082	0.1325	7.41%
	All set	0.9129	0.051	0.099	4.49%
SVR	Training dataset	0.7465	0.104	0.1617	8.12%
	Testing dataset	0.7459	0.188	0.2386	13.91%
	All set	0.7464	0.120	0.179	9.24%

. Both models have better statistical indicators on the training dataset compared to the testing dataset. As shown in Table 3, the improved SVR model has a higher R² and smaller statistical indicators (MAPE, MAE, and RMSE), and the generalization ability of the improved SVR model is significantly better than that of the SVR model in both the training and testing datasets. The statistical indicators of the improved SVR model on the training set are R² = 0.9129, MAE = 0.051, RMSE = 0.099, MAPE = 4.49%, and the corresponding statistical indicators of the SVR model are R² = 0.7464, MAE = 0.120, RMSE = 0.179, and MAPE = 9.24%. Compared with the original SVR model, the R² of the improved SVR model increases by 22.3% on average, the RMSE reduces by 44.693% on average, the MAPE reduces by 51.41% on average, and the MAE reduces by 57.50% on average. Compared with the SVR model, the improved SVR model has better robustness and fault tolerance due to the use of the best-performing hyperparameter combination.

4.3. Validation of the Performance

In order to substantiate the superiority of the proposed improved SVR model, it is meaningful to compare the proposed model with existing ML models. Part of the data in this dataset was selected, and

Table 4. Comparison of the prediction results of each model.

No.	FS
	True Values	BP Model [44]		GA-BP Model [44]		GP Model [76]		Proposed Model
		Predicted Results	Errors (%)	Predicted Results	Errors (%)	Predicted Results	Errors (%)	Predicted Results	Errors (%)
1	0.96	1.03	7.2917	1.0928	13.8333	1.0258	6.8542	1.05678	10.0813
2	1.15	1.16	0.8696	1.2543	9.0696	1.3147	14.3217	1.23214	7.1426
3	1.34	1.44	7.4627	1.2218	−8.8209	1.3696	2.2089	1.35395	1.0410
4	1.20	1.14	−5	1.2022	0.1833	1.1997	−0.025	1.20637	0.5308
5	1.55	1.62	4.5161	1.5005	−3.1936	1.4207	−8.3419	1.52003	−1.9336
6	1.45	1.41	−2.7586	1.4570	0.4828	1.3880	−4.2759	1.43555	−0.9966
7	1.31	1.47	12.2137	1.2678	−3.2214	1.3262	1.2366	1.31395	0.3015
8	1.49	1.6	7.3826	1.5594	4.6577	1.4674	−1.5168	1.48754	−0.1651
9	1.20	1.18	−1.6667	1.1977	−0.1917	1.2107	0.8917	1.18327	−1.3942
10	1.52	1.64	7.8947	1.4883	−2.0855	1.4670	−3.4868	1.50523	−0.9717
11	1.2	1.14	−5	1.2063	0.525	1.2654	5.45	1.21410	1.175
R2		0.891		0.851		0.866		0.974
RMSE		0.086		0.058		0.046		0.023

lists the FS prediction results and errors of the BP [44], GA-BP [44], and Gaussian process (GP) models [76]. From the prediction error, we can see that the proposed model has the best performance, with the largest R² (R² = 0.974) and the smallest RMSE (RMSE = 0.023). In summary, these findings imply that the proposed model can be regarded as a suitable predictive model for forecasting FS.

4.4. Feature Importance Score

A feature importance score is used to determine feature relevance and model interpretability, which is conducive to feature selection and improves model robustness to some extent [90]. In addition, for the problem of slope stability assessment, the feature importance score can also reflect the relative influence degree of various factors on slope stability from the side, so as to provide an effective reference for the formulation of slope safety protection programs. The feature importance score of each input is visualized as shown in Figure 8.

As can be seen, each of these six input features has a non-negligible feature importance score. Obviously, the unit weight γ is the most sensitive feature affecting slope stability, accounting for 29% of the feature importance scores of input features. The feature importance scores of the remaining input features are ranked from highest to lowest as follows: c (0.22) > φ (0.18) > α (0.16) > r_u (0.13) > H (0.09). Although c, φ, and r_u are all features from the perspective of mechanical properties, their effects on slope stability are quite different. These features form the basic modeling parameters for engineering projects and contribute to the study and prediction of FS. The analysis results can provide a reference for slope design with similar conditions. It is noteworthy that varying feature importance scores may arise when employing different ML models and datasets. With the emergence of more representative slope examples in the future, more valuable results can be obtained.

4.5. Contribution and Limitations

The primary innovation of this paper is the proposal of the improved SVR model applied to slope stability prediction. Slope cases from published literature are adopted to verify the proposed model. The new model relies on the powerful learning and reasoning functions of machine learning to directly estimate the FS and does not need to establish a complex mechanical calculation model. The input data exists objectively or is easy to measure, and the FS can be effectively estimated. The new model is simple, practical and economical, and easy to master by engineers and technicians and can provide reliable results for similar research. The proposed improved SVR model has great potential to be more widely used in other geotechnical engineering practices, where nonlinear problems are widely encountered.

It is worth noting that external factors such as seismic effects or other human factors that induce slope instability failure were not considered in the current study. Due to statistical data collection limitations, it was not possible to consider all features comprehensively [91,92]. Ignoring other features that affect slope stability, such as engineering disturbances, rainfall, vegetation cover, and geological structure, is an obvious limitation of this paper. Other features, such as dynamic load and scale effect, can be added in subsequent studies to further enhance the prediction ability of the model. In addition, the dataset is still relatively small. When more data become available, the performance of the proposed improved SVR model will improve. An important aspect that needs to be further explored is the scope of application of the proposed improved SVR model, particularly in relation to landslide types. This method has been shown to be suitable for the slope cases collected in Section 2. Whether this method might be applicable to other types of landslides remains to be explored. Another direction that may be extended is the potential application of the proposed method in dynamic analysis, obtaining a slope safety factor and changing some values of the input predictors along time. Last but not least, despite the growing familiarity with machine learning methods, the application of these methods to solve practical engineering problems is still relatively limited. Therefore, it is of great practical value to develop and utilize related intelligent systems in the future.

5. Conclusions

Slope stability plays an important role in engineering design and excavation. In this paper, a grid search method is introduced to optimize the hyperparameters of the model, and an improved SVR model is proposed to predict the FS. Learning and testing was conducted on 166 groups of slope cases, and an R² of 0.9 was achieved, which demonstrates that the proposed improved SVR model has strong practicability in FS prediction. The model was compared with the original SVR model and three widely used ML models to evaluate its predictive performance. Performance evaluation results demonstrate that this new model is superior to other models in R², RMSE, MAPE, and MAE. There is no doubt that once the proposed model is further well developed (such as using embedded software for secondary development), it will be more beneficial to engineering.

In this study, six input features were considered to predict FS, including γ, c, φ, α, H, and r_u. Through the feature importance score, it was found that the most important feature was γ. There was little difference in the importance scores of each input feature, and they all had a non-negligible effect on slope stability. Therefore, all of them should be taken into account when predicting slope stability. It is worth mentioning that the quantitative results obtained from the feature importance score are only applicable to the dataset adopted in this study.