Slope Stability Prediction Using k-NN-Based Optimum-Path Forest Approach

Abstract

Slope instability can lead to catastrophic consequences. However, predicting slope stability effectively is still challenging because of the complex mechanisms and multiple influencing factors. In recent years, machine learning (ML) has received great attention in slope stability prediction due to its strong nonlinear prediction ability. In this study, an optimum-path forest algorithm based on k-nearest neighbor (OPF_k-NN) was used to predict the stability of slopes. First, 404 historical slopes with failure risk were collected. Subsequently, the dataset was used to train and test the algorithm based on randomly divided training and test sets, respectively. The hyperparameter values were tuned by combining ten-fold cross-validation and grid search methods. Finally, the performance of the proposed approach was evaluated based on accuracy, F₁-score, area under the curve (AUC), and computational burden. In addition, the prediction results were compared with the other six ML algorithms. The results showed that the OPF_k-NN algorithm had a better performance, and the values of accuracy, F₁-score, AUC, and computational burden were 0.901, 0.902, 0.901, and 0.957 s, respectively. Moreover, the failed slope cases can be accurately identified, which is highly critical in slope stability prediction. The slope angle had the most important influence on prediction results. Furthermore, the engineering application results showed that the overall predictive performance of the OPF_k-NN model was consistent with the factor of safety value of engineering slopes. This study can provide valuable guidance for slope stability analysis and risk management.

1. Introduction

Slope instability is a global geological problem, which is one of the three major geological problems in nature besides earthquakes and volcanoes. Many geotechnical projects, such as open-pit mining, mountain roads, tailings dams, and landfills, are seriously threatened by slope instability. A serious slope instability disaster can cause casualties, building damages, and huge economic losses. For example, on 20 December 2015, a catastrophic landslide occurred at the Hong’ao landfill in Shenzhen, China, resulting in 77 deaths, 33 buildings buried, and direct economic losses of more than 880 million RMB [1]. On the evening of 11 March 2017, a landslide at the Koshe landfill in Ethiopia’s capital, Yah, caused 113 deaths and more than 80 people missing [2]. Due to heavy rainfall on 18 October 2020, a landslide occurred in Vietnam’s Quang Tri province, claiming the lives of 22 soldiers [3]. Because of its serious consequences, predicting the risk of slope instability is crucial and plays a significant role in disaster prevention.

The prediction methods of slope stability can be classified into four categories. The first one is instrumental monitoring technology. Currently, many on-site monitoring techniques of slope deformation have been applied to monitor the early warning signs of slope instability. For example, Zhang et al. [4] used distributed fiber optic strain sensors to monitor the shear displacement in the Three Gorges Reservoir region in China, and two potential circular sliding surfaces were successfully identified. Dixon et al. [5], Shiotani [6], and Codeglia et al. [7] adopted the acoustic emissions (AE) technology to monitor the signals generated by the fracture of soil and rock materials in the slope. By analyzing the relationship between AE characteristics and slope deformation, AE-based criteria were used to evaluate the long-term stability of slopes. In addition, some other techniques, such as remote sensing [8], terrestrial laser scanning [9], synthetic aperture radar [10], and time domain reflectometry [11], were applied to slope stability monitoring. These technologies have relatively high prediction accuracy because the precursor information of slope instability can be obtained directly, but the installation process is complicated, and the cost is high.

The second one is the theoretical analysis method. It is proposed from the view of mechanical mechanisms. Many theoretical and analytical approaches have been used to analyze slope stability, such as the limit equilibrium method (LEM) [12], the strength reduction method (SRM) [13], and the limit analysis method [14]. The factor of safety (FOS), calculated by the ratio of resisting force to driving force, is used to evaluate the stability of the slope. When the value of FOS is larger than 1, the slope is stable; otherwise, it is unstable [15]. Faramarzi et al. [16] employed LEM to calculate the FOS and analyzed the rock slope stability of the Chamshir dam pit. Liu [17] adopted the SRM to obtain the FOS of the established slope model. Mbarka et al. [18] combined the Monte Carlo approach, LEM, and SRM for the reliability analysis of homogeneous slopes with circular-type failure. Although the theoretical and analytical methods are simple, they are unsuitable for slopes with complex conditions due to the simplified formulas and assumptions.

The third one is the numerical simulation technique. With the rapid development of numerical simulation methods, finite element method (FEM) [19], boundary element method [20], discrete element method [21], numerical manifold method [22], and other methods have been widely used in slope stability analysis. Sun et al. [23] simulate the progressive failure process of jointed rock slopes based on the combined finite-discrete element method. Ma et al. [24] analyzed the slope stability under a complex stress state with saturated and unsaturated seepage using the fast Lagrangian analysis of continua. Wei et al. [25] investigated the kinetic features of slope instability based on particle flow code. Haghnejad et al. [26] analyzed the effect of blast-induced vibration on slope stability using dynamic pressure in three dimensions distinct element codes. Song et al. [27] adopted an improved smoothed-particle hydrodynamics method to calculate the slope safety factor. Zhang et al. [28] adopted a realistic failure process analysis to evaluate the stability and investigated the failure mode of the high rock slope during excavations. In addition, some researchers have integrated numerical simulation and mathematical methods to analyze the slope stability. For example, Dyson and Tolooiyan [29] adopted FEM and Monte Carlo to determine the FOS and damage probability of slopes. Although the numerical simulation methods are convenient to operate, the accuracy strongly depends on constitutive models and mechanical parameters [30].

The fourth one is the machine learning (ML) algorithm. With the accumulation of slope cases, some researchers attempted to develop slope stability prediction models using ML algorithms. There are two types of predicted outputs: FOS and stability status. Lu and Rosenbaum [31] adopted an artificial neural network to estimate the FOS and SS on 46 slope cases collected by Sah et al. [32]. Based on the same database, Samui [33] and Yang et al. [34] used a support vector machine (SVM) and genetic programming to determine FOS, respectively. Amirkiyaei and Ghasemi [35] constructed two tree-based models to assess circular-type failure slopes based on 87 cases. Zhou et al. [36] collected 221 slope cases and employed the gradient-boosting machine to predict the SS. Wang et al. [37] hybridized a genetic algorithm with a multi-layer perceptron to predict FOS using 630 cases. In addition, several researchers performed a comparative analysis of multiple ML algorithms. Hoang and Tien Bui [38] carried out a comparative study of SS prediction using a radial basis function neural network, an extreme learning machine, and least squares SVM. Mahmoodzadeh et al. [39] adopted Gaussian process regression, support vector regression, decision trees (DT), long-short-term memory, deep neural networks, and k-nearest neighbors (k-NN) to determine FOS. All the above ML algorithms performed well on slope stability prediction. However, a large number of slope stability cases are required to improve its credibility.

Compared with other approaches, ML algorithms can obtain reliable prediction results by establishing the nonlinear relationship between input and output. It is a promising method for predicting slope stability. But to date, there is no one ML algorithm that can be applied to all slope engineering conditions under the consensus of the geotechnical industry. Accordingly, it is meaningful to investigate more robust ML algorithms to achieve better prediction results. Recently, the optimum-path forest (OPF) algorithm has been successfully applied in many fields, such as face recognition [40], Parkinson’s disease identification [41], laryngeal cancer pathology detection [42], land use classification [43], and network intrusion detection [44]. However, the OPF algorithm is susceptible to outliers. In response to this deficiency, Papa et al. [45] proposed the OPF algorithm based on k-NN (OPF_k-NN), and the discriminative performance of the OPF model was improved. In combination with the k-NN algorithm, the OPF_k-NN algorithm can provide better performance for classification tasks by leveraging the topological properties of the data [46]. Compared to other classification algorithms, the OPF_k-NN algorithm has several advantages, including (1) it is free of hyperparameters, (2) it does not assume separability of the feature space, (3) it has a unique feature selection and classification mechanism that can effectively handle the high-dimensional and nonlinear data with outliers, (4) and its training step is usually much faster than traditional ML approaches.

Considering that the OPF_k-NN has great predictive performance and has not yet been employed to predict the stability of slopes, this study aims to investigate the feasibility of OPF_k-NN for predicting slope stability. In addition, a comparison against OPF, radial basis function support vector machine (RBF-SVM), random forest (RF), DT, k-NN, and logistic regression (LR) classifiers is performed.

2. Methodology

2.1. k-NN Based OPF Classifier

The OPF is a graph-based classifier [47,48]. Its classification principle is to denote the training samples as nodes and connect them by path. Then, the optimal path tree (OPT) is constructed by executing the shortest path algorithm on the graph. Finally, the test sample is mapped onto the OPT, and its class is determined. Figure 1 shows the schematic diagram of the OPF-based classifiers. The nodes with different colors in the set S represent different classes, and the nodes outside the set S are the samples to be classified. A series of adjacent nodes are defined as path π. Among all paths, the one with the maximum path-cost function f(π_t) is called OPT, and all OPTs constitute OPF. There are three different classes in Figure 1; the blue sample s is the root node of the OPT where sample t is located, so sample t is classified as blue.

The OPF_k-NN is a variant of the OPF algorithm, and the main difference between them is the adjacency of the samples in the training set. The latter is to construct a complete graph, while the former is to construct a k-NN graph [45]. The OPF_k-NN algorithm is divided into training and classification phases.

2.1.1. Training Phase

The first step is to construct a k-NN graph G_k based on the training set Z₁. The sample s is weighed by a probability density function ρ(s):

$$\rho \left(s\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}\left|G_k^{\ast }\left(s\right)\right|}{\displaystyle \sum _{\forall t\in G_k^{\ast }\left(s\right)}\exp \left(\frac{-d^2\left(s,t\right)}{2{\sigma }^2}\right),}$$

(1)

where $\sigma =\frac{d_f}{3}$, d_f is the maximum arc weight in G_k, and d(s, t) is the distance between sample s and sample t.

The second step is to calculate the path cost function f_min, which is defined as:

$$\begin{array}{l}{f}_{\min }\left(\langle t\rangle \right)=\{\begin{array}{cc}\rho \left(t\right)& \mathrm{if} t\in S\hfill \\ \rho \left(t\right)-1& \mathrm{otherwise}\hfill \end{array}\\ f_{\min }\left({\pi }_s\cdot \langle s,t\rangle \right)=\min \left\{f_{\min }\left({\pi }_s\right),\rho \left(t\right)\right\}\end{array},$$

(2)

According to the method proposed by Papa et al. [50], the k value of k-NN is determined by maximizing the accuracy of the training set in the range [1, k_max]. The value of k_max defaults to 5. After determining the value of k, the algorithm is applied to retrain the classifier. The function f_min is replaced by $f_{\min }^{\prime }$, which is defined as:

$$\begin{array}{l}{f}_{\min }^{\prime }\left(\langle t\rangle \right)=\{\begin{array}{cc}\rho \left(t\right)& \mathrm{if} t\in S\hfill \\ \rho \left(t\right)-1& \mathrm{otherwise}\hfill \end{array}\\ f_{\min }^{\prime }\left({\pi }_s\cdot \langle s,t\rangle \right)=\{\begin{array}{cc}-\infty \hfill & \mathrm{if} \lambda \left(t\right)\ne \lambda \left(s\right)\\ \min \left\{f_{\min }^{\prime }\left({\pi }_s\right),\rho \left(t\right)\right\}\hfill & \mathrm{otherwise}\end{array}\end{array}.$$

(3)

Figure 2 is the schematic diagram of the training phase, where Figure 2a indicates the k-NN graph generated from the training set, Figure 2b represents the minimum spanning tree calculated by the k-NN graph, Figure 2c denotes the two samples of different colors labeled as prototype samples (marked by black dashed circles), and Figure 2d signifies the OPF_k-NN classifier composed by all the OPTs. The red squares and green circles represent different classes, respectively.

2.1.2. Classification Phase

After training the OPF_k-NN classifier, the sample t in the test set Z₂ is classified. The k-NN is first calculated from Z₁ to a testing sample t. Then, it is verified which sample s ∊ Z₁ satisfies the equation below:

$$V\left(t\right)=\max \left\{\min \left[V\left(s\right),\rho \left(t\right)\right]\right\}\forall s\in Z_1$$

(4)

Figure 3 indicates the classification process of OPF_k-NN. The blue triangle is the sample to be classified. Figure 3a shows that the blue triangle is connected to the k-nearest training samples in the generated OPF, and Figure 3b illustrates that the triangle is conquered by the samples of the red squares class and labeled as red.

2.2. Proposed Approach

Figure 4 depicts the flowchart of the proposed approach. First, due to the different units of indicators and the diversity of data distribution, the raw data is pre-processed. The dataset is standardized using a Gaussian distribution with zero mean and unit standard deviation. Subsequently, 80% of samples are used for training, and the remaining 20% are adopted for testing [51,52]. For the k-NN, RBF-SVM, RF, DT, and LR algorithms, the grid search and ten-fold cross-validation (CV) methods are used to select the optimal hyperparameters. Finally, the test set is predicted, and the optimal classifier is determined according to the evaluation metrics.

The OPF_k-NN and OPF are implemented based on the Python library “opfython” [53], and the k-NN, RBF-SVM, RF, DT, and LR are conducted on the Python library “scikit-learn” [54]. All experiments are conducted using a Windows1064 bits computer with 8Gb of RAM running an Intel^® Core™ i7-9700F CPU @ 3.00 GHz × 2.

If the predictive performance of our proposed approach is acceptable, it can be used for engineering applications in several ways. For example, it can be integrated into slope monitoring systems to provide real-time alerts for potential instability. The model can also be used to evaluate slope stability during the design phase of construction projects to ensure the safety and stability of the slope. Additionally, the model can be applied to slope stability analysis and risk management, which can be used by geotechnical engineers in various projects related to slope instability.

2.3. Performance Evaluation Metrics

In this study, several metrics are used to evaluate the performance of classifiers and figure out the optimal classifier for slope stability prediction [30,55].

A confusion matrix, which can also be called a likelihood table or error matrix, is used to visually represent whether the performance is ideal or not.

Table 1. Confusion matrix for slope stability prediction.

Actual Condition	Predicted Condition
	Stable	Failed
Stable	True positive	False negative
Failed	False positive	True negative

shows the confusion matrix for the slope-stability prediction, where true positive (TP) means the number of stable cases predicted correctly, false positive (FP) means the number of stable cases predicted incorrectly, true negative (TN) means the number of failed cases predicted correctly, and false negative (FN) means the number of failed cases predicted incorrectly. According to Table 1, true negative rate (TN/(TN + FP)) and false positive rate (FP/(FP + TN)) can be defined.

Accuracy indicates the ratio of the cases correctly predicted to the total cases, which can be calculated by: accuracy = (TP + TN)/(TP + TN + FP + FN).

F₁-score indicates the harmonic mean of precision and recall, which can be calculated by: F₁-score = 2precision · recall/(precision + recall), where precision = TP/(TP + FP), recall = TP/(TP + FN).

The area under the curve (AUC) is defined as the area under the receiver operating characteristic (ROC) curve, which is commonly used to evaluate the performance of classifiers. Bradley [56] proposed classification criteria of AUC as follows: not discriminating (0.5–0.6), poor (0.6–0.7), fair (0.7–0.8), good (0.8–0.9), and excellent (0.9–1).

Computational burden is used to evaluate the computational efficiency of algorithms. The mean and standard deviation of computation time are used as the evaluation metrics in this study.

2.4. Hyperparameter Tuning

In general, the performance of most ML algorithms is highly dependent on the hyperparameters. There are several hyperparameter tuning methods, such as manual search, grid search, random search, Bayesian optimization, gradient-based optimization, and evolutionary optimization [57]. In this study, the grid search algorithm is combined with the k-fold CV method to select the optimal hyperparameters.

The grid search algorithm is to grid the hyperparameters in a fixed range in equal steps, compare all hyperparameter combinations exhaustively, and then select the optimal hyperparameters. To avoid the risk of overfitting or selection bias in the model, the k-fold CV method is used in the hyperparameter tuning process, illustrated in Figure 5. The original training set is randomly split into k folds, of which k − 1 folds are used as the training sub-set, and the remaining fold is used as the validation set in turn. Then, the average accuracy of k rounds is calculated to evaluate the performance and determine the optimal hyperparameters [58]. In this study, k was selected as 10 after considering the calculation time and variance [59].

3. Data Collection

3.1. Dataset Description

The failure surfaces of the slopes are prone to occur near the potential slip surface. Because of the excavation at the foot of the slope or water seepage at the top of the slope, the shear stress on the potential slip surface exceeds the shear strength, causing the local slope instability, as shown in Figure 6. A large number of engineering cases and theoretical analyses indicate that there are three main aspects that affect slope stability: the physical–mechanical properties of the potential slip surface, the basic geometrical parameters, and the external triggering factors. [12,18,60,61,62]. Considering the independent correlation between indicators and the easy availability of indicator values, six indicators were selected in this study, including the unit weight (γ), the cohesion (c), the internal friction angle (φ), the slope angle (β), the slope height (H), and the pore pressure ratio (r_u). The detailed descriptions of these indicators are displayed in

Table 2. Descriptions of input indicators.

Indicator	Description	Measurement Method
γ (kN/m3)	It indicates the weight of soil/rock per unit volume.	It can be measured by performing the standard mass volume method, mercury displacement method, or gravimeter method in the laboratory.
c (kPa)	It indicates the attraction between molecules on the surface of adjacent material particles within the soil/rock.	It can be determined by performing direct shear tests and triaxial compression tests in the laboratory.
φ (°)	It indicates a measure of the ability of a unit of soil/rock to withstand shear stress.	It can be determined by performing direct shear tests and triaxial compression tests in the laboratory.
β (°)	It indicates the angle between the slope plane and the slope bottom.	It can be measured in the field by an inclinometer.
H (m)	It indicates the vertical distance from the slope bottom to the slope top.	It can be measured in the field using a surveying instrument such as a total station.
ru	It is defined as the ratio of the pore pressure and normal stress at a certain point within a slope.	It can be measured by installing pore water piezometers on-site or by performing immersion tests or infiltration tests in the laboratory.

.

In this study, a database of 404 slopes with failure risk from various countries was collected (available in “Appendix A”) [32,36,57,63,64,65,66,67,68,69,70,71,72]. There are two statuses of slope stability: stable (207 cases) and failed (197 cases). Among them, most of the failed slopes were circular-type failures. The distribution of slope SS on the overall dataset is given in Figure 7, and the statistical values of data samples are illustrated in

Table 3. Statistical values of slope stability dataset.

Value Type	γ (kN/m3)	c (kPa)	φ (°)	β (°)	H (m)	ru
Minimum	10.06	0	0	4.24	3.45	0
Median	21.38	29.7	28.27	34.03	51	0.2
Maximum	31.3	300	57.36	59.35	565	0.75
Mean	21.69	39.38	27.74	34.19	84.26	0.18
Standard	3.84	40.54	9.63	10.86	94.97	0.17

.

3.2. Dataset Analysis

The violin plots of six indicators are shown in Figure 8. They were a combination of box plots and density plots and indicated the overall distribution of the dataset. For each violin plot, the white dot in the center was the median of the samples, the top and bottom of the thick black line represented the third and first quartile of the samples, and the top and bottom of the thin black line indicated the upper and lower adjacent value. From Figure 8, it can be seen that the distribution of γ, φ, β was relatively balanced, and the medians were basically in the middle of the violin plots. While for c, H, r_u, there were some individual outliers.

The heatmap of the Pearson correlation coefficient between each indicator is shown in Figure 9. According to Figure 9, all correlation coefficients were less than 0.5, and the highest correlation was only 0.41, which indicated that the correlation between indicators was poor. Therefore, all indicators were relatively independent and important for predicting slope stability.

To visualize the distribution of the dataset, the correlation pair plots of the two slope SS were displayed in Figure 10. The distribution plots of these six indicators were shown on the diagonal line, and the correlation scatter plots between indicators were shown on the non-diagonal line. It can be seen that the differences in the distribution of indicators for both slope statuses were slight, and there was no apparent correlation among the indicators. Therefore, it was difficult to classify the slope SS only using one indicator, and the effect of all indicators should be incorporated for better accuracy.

4. Results and Analysis

4.1. Results of Hyperparameters Tuning

The average accuracy of ten-fold CV corresponding to different hyperparameters for k-NN, LR, DT, RF, and RBF-SVM algorithms is shown in Figure 11. According to Figure 11, the overall performance can be observed. With the increase of hyperparameter values, the average accuracy of LR decreased, but the other models had several peaks. Compared with other models, the results of RF were more stable. Based on the best average accuracy of ten-fold CV, the optimal hyperparameter values were determined. The scope, interval, and final optimization results of hyperparameter values are indicated in

Table 4. Results of hyperparameters tuning.

ML Algorithms	Hyperparameters	Scope of Values	Interval of Values	Optimal Values
k-NN	n_neighbors	(1, 31)	1	7
LR	Inverse of regularization strength C1	(0.1, 10)	0.1	0.1
DT	max_depth	(1, 10)	1	7
	min_samples_leaf	(1, 10)	1	3
RF	n_estimators	(1, 101)	10	31
	max_depth	(1, 20)	1	11
RBF-SVM	gamma	(0.01, 0.6)	0.01	0.55
	Penalty coefficient C2	(3, 4)	0.1	3.3

.

4.2. Models Comparison and Evaluation

After the hyperparameters were tuned, these seven ML algorithms were used to predict slope stability based on the test set. The confusion matrix, accuracy, and F₁-score were calculated to compare the performance of each algorithm, which were illustrated in

Table 5. Confusion matrix, accuracy, and F₁-score of the classifiers.

Classifiers	Actual Condition	Predicted Condition		Accuracy	F1-Score
		Stable	Failed
OPFk-NN	Stable	37	4	0.901	0.902
	Failed	4	36
OPF	Stable	38	3	0.876	0.884
	Failed	7	33
RF	Stable	37	7	0.827	0.841
	Failed	7	30
k-NN	Stable	36	8	0.815	0.828
	Failed	7	30
RBF-SVM	Stable	35	9	0.802	0.814
	Failed	7	30
DT	Stable	33	11	0.765	0.776
	Failed	8	29
LR	Stable	28	16	0.679	0.683
	Failed	10	27

. It can be observed that OPF_k-NN performed best with the highest accuracy of 0.901, followed by OPF, RF, k-NN, RBF-SVM, and DT with an accuracy of 0.876, 0.827, 0.815, 0.802, and 0.765, respectively. LR performed worst with an accuracy of 0.679. Furthermore, the rank was the same when using the F₁-score. Therefore, based on the overall prediction performance, the rank was OPF_k-NN > OPF > RF > k-NN > RBF-SVM > DT > LR.

The ROC curves and AUC values of these seven classifiers are presented in Figure 12. It can be seen that the ROC curve of the OPF_k-NN classifier was closer to the left and upper axes than others, indicating better performance. The AUC values of OPF_k-NN, RBF-SVM, RF, OPF, k-NN, DT, and LR were 0.895, 0.885, 0.876, 0.870, 0.783, and 0.720, respectively. According to the AUC classification criterion mentioned in Section 2.3, only OPF_k-NN performed excellently, RBF-SVM, RF, OPF, and k-NN performed well, while DT and LR performed fair.

The average computation time for each classifier during the training and testing phases over 20 runs was calculated, as listed in

Table 6. Average computation time of classifiers over 20 runs.

Time	OPFk-NN	OPF	RF	RBF-SVM	DT	LR	k-NN
Train	0.915 ± 0.027	0.322 ± 0.034	117.751 ± 1.163	82.326 ± 1.149	2.162 ± 0.067	0.402 ± 0.035	0.187 ± 0.011
Test	0.042 ± 0.002	0.097 ± 0.004	0.055 ± 0.003	0.171 ± 0.005	0.008 ± 0.001	0.008 ± 0.001	0.011 ± 0.002
Total	0.957 ± 0.026	0.419 ± 0.033	117.806 ± 1.164	82.497 ± 1.149	2.170 ± 0.067	0.410 ± 0.036	0.198 ± 0.011

. The results were presented in the following format: x ± y, where x and y indicated the average time and standard deviation, respectively—noted that the values in bold indicated the minimum time consumed. It can be observed that the k-NN took the least time in the training phase, followed by LR, OPF, OPF_k-NN, RBF-SVM, and RF. In the testing phase, the time consumed by each classifier was not significantly different, and the difference between the maximum and minimum values was less than 0.2 s. For the total time, the rank was k-NN > LR > OPF > OPF_k-NN > DT > RBF-SVM > RF. The total computation time of the OPF_k-NN classifier was less than 1 s.

4.3. Relative Importance of Indicators

The relative importance of indicators was significant for the design of support structures in slope engineering. In this study, the relative importance of each indicator was calculated by combining the OPF_k-NN model with the permutation feature importance technique [73]. The permutation feature importance is a model inspection technique available in the Python library “scikit-learn” [54]. Values of indicators were shuffled in turn within the test set, the slope stability prediction results were generated by the OPF_k-NN model, and the accuracy changes were recorded. Then, the prediction accuracy changes of indicators were ranked, and the relative importance was derived. As shown in Figure 13, the slope angle was the most important indicator with an importance value of 30.5%, followed by internal friction angle (22%), cohesion (19.7%), unit weight (12.3%), slope height (7.93%), and pore pressure ratio (7.63%).

5. Discussions

The prediction of failed cases is particularly important, which may lead to the development of slope instability if predicted incorrectly [74]. Therefore, the false positive rate and true negative rate were presented together in Figure 14. It can be seen that the false positive rate and true negative rate of RBF-SVM, k-NN, and RF were the same, and the OPF_k-NN had the largest true negative rate and the lowest false positive rate. From this view, the OPF_k-NN classifier performed better. The reason is that the OPF_k-NN algorithm can effectively process high-dimensional and nonlinear slope data with outliers, improve the data quality of the model in the training phase, and predict the failed slope cases more accurately.

When the trade-off between AUC and the computational burden was considered, the OPF_k-NN classifier was the most prominent because it demonstrated the optimal performance (AUC = 0.901) in less computation time (total time < 1 s) among the seven classifiers. It is worth noting that the OPF_k-NN classifier was pretty much faster than RBF-SVM (86.2 times faster) and RF (123.1 times faster), although the difference in their performance was not significant. Therefore, the OPF_k-NN classifier achieved the best trade-off between performance and efficiency.

According to the importance scores, all indicators were non-negligible for slope stability prediction. The physical–mechanical properties had the greatest influence on the slope stability (φ = 22%, c = 19.7%, γ = 12.3%), followed by the geometrical parameters (β = 30.5%, H = 7.93%). Some measures can be adopted to improve the slope stability from two directions. One is to optimize the slope geometry parameters, especially the slope angle. Another is to improve the physical–mechanical properties by using grouting-reinforcement techniques.

Although the OPF_k-NN approach obtained excellent results in the slope stability prediction, there are also some limitations:

(1) More indicators should be considered. Although the six indicators in this study affect the slope stability significantly, other factors such as excavation, the properties of clay minerals, vegetation coverage, earthquake, and rainfall also have an effect on the slope stability. It is significant to analyze the influences of these indicators on the prediction results;

(2) The dataset is relatively small. The performance of ML algorithms greatly depends on the quantity and quality of data. Although the OPF_k-NN algorithm performs well on this dataset, a better dataset might further improve the predictive performance. Therefore, it is necessary to build a larger slope database;

(3) Slopes are typically composed of multiple layers of various geotechnical materials whose properties and spatial distribution can significantly affect slope stability. As the number of slope failure cases increases, a comprehensive and diverse slope dataset should be expanded in future work. Such efforts are crucial for advancing the field of geotechnical engineering and ensuring the safety of human lives and infrastructure;

(4) The safety factor of slope stability can reflect the percentage of slope instability, and the slope stability analysis can be better considered a regression problem. Therefore, it is necessary to compile relevant data and develop relevant ML models for slope FOS value estimation in future work.

6. Engineering Application

In order to further verify the reliability of the proposed OPF_k-NN model, it was necessary to apply it to evaluate the stability of engineering slopes. For this, eight typical slopes were collected from the Jing-xin expressway in Hebei Province, China, where landslides frequently occurred [75].

The FOS values of these eight slopes and the estimation results of the OPF_k-NN model were recorded in

Table 7. Predictive results of OPF_k-NN model on engineering slopes.

Slopes	γ (kN/m3)	c (kPa)	φ (°)	β (°)	H (m)	ru	FOS	Status
1	22.4	20.0	27.0	30.0	54.0	0.12	1.208	Stable
2	21.4	31.5	42.0	34.0	18.0	0.23	2.448	Stable
3	19.0	50.0	32.0	42.0	26.0	0.17	1.786	Stable
4	19.6	17.8	29.2	41.2	50.0	0.31	0.979	Failed
5	20.2	16.7	22.3	42.4	26.6	0.47	0.869	Failed
6	20.4	25.0	20.4	35.0	65.9	0.42	0.833	Failed
7	20.0	20.0	36.0	45.0	50.0	0.14	1.102	Stable
8	23.0	18.3	25.2	39.6	61.2	0.30	0.824	Failed

. It can be seen that the overall prediction performance of the OPF_k-NN model was consistent with the FOS values of the slopes.

7. Conclusions

Slope stability prediction is a crucial task in geotechnical engineering. This study investigated the performance of the OPF_k-NN algorithm for the stability prediction of slopes. A total of 404 historical slope cases with failure risk from various countries were collected after considering the slope damage mechanism and geological conditions simultaneously. The OPF_k-NN, OPF, RBF-SVM, RF, k-NN, DT, and LR were used to evaluate and compare the predictive performance. To avoid the risk of overfitting or selection bias, ten-fold CV and grid search methods were selected to tune the hyperparameters. Overall, the prediction results of the OPF_k-NN algorithm were better and more reliable, and its prediction accuracy and F₁-score were 0.901 and 0.902, respectively. According to the ROC curves and AUC values, the performance rank of the seven classifiers was OPF_k-NN > RBF-SVM > RF > OPF > k-NN > DT > LR. In addition, the OPF_k-NN achieved the highest TNR and the lowest FPR, which indicated that it could predict failed slope cases better. After considering the total calculation time, the OPF_k-NN classifier achieved the optimal trade-off between performance and efficiency. Based on the importance scores of indicators, the slope angle was the most influential indicator on prediction results. Furthermore, the engineering application results showed that the overall predictive performance of the OPF_k-NN model was consistent with the FOS value of engineering slopes.

In the future, more parameters such as excavation, the properties of clay minerals, geological formation, vegetation coverage, earthquake, and rainfall can be considered so that the feasibility of the OPF_k-NN classifier can be further validated using more comprehensive and diverse slope datasets. In addition, the proposed methodology can be recommended for the application of other mining and geotechnical engineering projects, such as rockburst risk prediction and pillar stability prediction.