Ground-penetrating radar (GPR) is a subsurface exploration method that uses high-frequency to microwave electromagnetic waves to detect reflections. The system works by emitting electromagnetic waves from a transmitter antenna directed towards the surrounding rock. The antenna is positioned on the rock’s surface, and as the waves propagate through the subsurface, they are reflected when encountering interfaces with varying electrical properties. The reflected waves are captured by the receiver antenna. The initial signal received from the transmitter to the receiver serves as the system’s time zero reference point. After processing the reflected signal, the depth of the reflection target is determined by halving the round-trip time and multiplying by the radar wave velocity in the corresponding medium. The characteristics of the reflection target, such as fractured rock masses, cavities, loose zones, weak interlayers, and cracks in unfavorable geological formations, are identified based on the intensity, shape, and variations of the reflected wave in both horizontal and vertical directions. The detection principle of geological radar is illustrated in Figure 2.

The Ground Penetrating Radar (GPR) provides critical input data for the 3DEC model, facilitating accurate simulation of slope stability. Through radar detection, structural information of the rock mass, such as fractures, joints, weak interlayers, and fault zones, can be obtained, which is essential for defining the mechanical properties in the 3DEC model. Additionally, GPR offers depth-related data, enabling precise determination of layer thickness, distribution, and the strength and stiffness of rock masses at different depths. Variations in the intensity and frequency of the reflected radar waves reveal the characteristics of weak interlayers and fracture zones, thereby aiding in the accurate representation of the rock mass heterogeneity in the 3DEC model. The collected data, along with field conditions, surrounding geological survey reports, and geological investigation data, are analyzed to provide a comprehensive interpretation of the geological conditions in the detected area, with particular focus on axial direction, amplitude, and frequency.

In this study, radar lines 1# and 2# were positioned on the safety platform of the open-pit slope. The combined results of the field geological survey and radar detection indicate the following: (1) The surrounding rock joints and fractures are generally well to moderately developed, and the rock mass exhibits a blocky structure; (2) At a depth of 30–35 m below the surface, significant variations in the frequency and amplitude of the reflected waves suggest poor rock integrity, well-developed joint fractures, and the presence of weak interlayers and fault zones. The field survey shows that these weak interlayers consist of chlorite mineral layers, where the minerals have fragmented into flaky structures. Some of these weak interlayers exhibit wide openings, and the minerals have been eroded by external forces, leading to the formation of fractures of varying widths. These observations are consistent with the findings from the geological radar survey. The detection results are presented in Figures 3a,b.

After determining the Geological Strength Index (GSI) value, rock mass mechanical parameters were analyzed based on the Hoek-Brown strength criterion. The final rock mass mechanical parameters for each zone are shown in Table 3. In the 3DEC calculation process, the rock mechanics parameters from Table 3 are selected as input variables for numerical calculations. These parameters, such as compressive strength and shear strength of the rock mass, are chosen based on field geological surveys and laboratory test data, accurately reflecting the physical and mechanical properties of the rock mass. This provides a reliable basis for subsequent slope stability analysis.

3DEC (3-Dimensional Distinct Element Code), developed by Itasca Consulting Group, Inc., is a three-dimensional numerical analysis software based on the Distinct Element Method (DEM). This program is specifically designed to simulate the mechanical behavior of discontinuous media and is equipped with an extensive library of material models.

During the model process, 3DEC can realistically represent the geometric characteristics of rock masses with structural planes in three-dimensional space. It offers distinct advantages in analysing rock mass deformation and failure that are governed by structural features. As a result, it has been widely applied in studies of the mechanical behavior of slopes, surrounding rock of underground excavations, and other jointed rock masses.

3DEC provides powerful and comprehensive pre-processing capabilities, including slope model construction, rock mass structural division, and structural plane net-work simulation. These features allow for a more accurate representation of the geo-metric characteristics of engineering slopes, making 3DEC an excellent tool for de-tailed investigations of three-dimensional stress responses in rock slopes. Lattice element-based modeling of cemented geomaterials provides further validation for using discontinuum-based tools, such as 3DEC, in simulating slope failure mechanisms with weak interlayers (Rizvi et al., 2020).

The basic principle of the strength reduction method is to divide the strength parameters of the slope soil-rock mass, namely, cohesion c and internal friction angle φ, by a reduction factor F simultaneously, obtaining a new set of values. These reduced values are then used as the new parameters for the calculation. The reduced shear strength parameters are as follows Equations 1, 2: C ′ = C F (1) φ ′ = arctan tan ⁡ φ F (2)

As shown in Figure 7, the value of F when Intensity line after discounting is tangent to the Mohr’s stress circle is the slope safety factor required for the calculation.

This calculation only presents a portion of the results. As shown in Figure 8, the schematic diagrams of the models with a slope angle of 60°, bench heights of 15 m and 17 m, and safety platform widths of 5 m, 4 m, and 3 m, respectively.

As shown in Figure 9a, the slope safety factor varies with different bench widths and heights. Specifically, Figure 9a presents the safety factor for a slope with a 60° angle, a bench height of 15 m, and a bench width of 5 m. Figure 9b shows the safety factor for the same slope angle and height but with a bench width of 4 m, while Figure 9c illustrates the safety factor for a bench width of 3 m. Figures 9d–f correspond to a slope with a height of 17 m and bench widths of 5 m, 4 m, and 3 m, respectively. According to the safety factor results in Table 5, as the bench width decreases, the safety factor gradually declines. Although increasing the bench height slightly reduces the safety factor, the change is minimal, suggesting that higher benches still maintain relatively good stability.

Support Vector Machine (SVM) is a widely used machine learning algorithm. Its main principle is to construct a decision surface—a classification hyperplane—that maximizes the separation boundary between positive and negative examples. There-fore, the SVM algorithm exhibits excellent generalization performance in classification problems. The algorithm’s principles and system structure are shown in Figures 10, 11.

If the training set T of the SVM is given in the following form are as follows Equations 3–6: T = x 1 , y 1 , x 2 , y 2 , . . . . . . , x l , y l ∈ X × Y ′ (3)

In the equation: x i ∈ X = R n , y i ∈ Y = 1 , − 1 , Where, i = 1,2,…,l.

Select an appropriate kernel function K (without loss of generality, we choose the Radial Basis Function (RBF) as the kernel) and parameter C to construct the optimization problem, as shown below. min α 1 2 ∑ i = 1 l ∑ j = 1 l α i α j y i y j K x i , x j − ∑ j = 1 l α j (4) K x , x i = exp − gx − x i 2   g > 0 (5) s . t . ∑ i = 1 l α i y i = 0 0 ≤ α i ≪ C   i , j = 1 , 2 , . . . , l (6)

In the equation, α represents the Lagrange multiplier.

Thus, the optimal solution α=(α1,α2,…,αl)T can be obtained. Select a positive component 0<αl < C to calculate the threshold are as follows Equations 7, 8. b * = y j − ∑ i = 1 l α i * K x i − x j (7)

The final decision function can be constructed as follows: f x = sgn ∑ i = 1 l α i * y i K x , x i + b * (8)

Overall, the SVM algorithm has the following advantages: it is well-structured according to the framework of Vapnik-Chervonenkis generalization theory; it possesses good robustness; and it can be applied to a wide range of functions.

When using the machine learning SVM algorithm for model training, the penalty coefficient ccc and the kernel function parameter g have a significant impact on the model’s accuracy. The penalty coefficient c represents the tolerance for errors. A larger c indicates a smaller tolerance for errors, which may lead to overfitting, while a smaller c may lead to underfitting. On the other hand, both excessively large and small values of c can weaken the model’s generalization ability. The kernel function parameter g affects the number of support vectors: as g increases, the number of support vectors decreases, and as g decreases, the number of support vectors increases. However, in practical applications of the SVM algorithm, the parameters c and g are often set based on empirical experience or system default values, lacking theoretical foundation, and are not necessarily the optimal solution for a specific problem.

Bayesian Optimization is a global optimization method based on Bayes' theorem and Gaussian Processes (GP), primarily used to optimize black-box functions, i.e., functions that do not have an analytical form, are computationally expensive, and are noisy. In general, Bayesian Optimization allows us to find the best solution with fewer trials in a smarter way, making it particularly suitable for problems where the “trial-and-error cost” is high.

Problem Definition:

Suppose we need to optimize an unknown objective function is as follows Equation 9: f : x → R (9)

Here, χ is the domain (typically a high-dimensional continuous space), and f x is the objective function to be optimized (the black-box function).

The goal is to find the optimal solution are as follows Equations 10, 11: x * = arg max f x x ∈ x (10) or x * = arg min f x x ∈ x (11)

Since f x is unknown and computationally expensive, ordinary grid search or gradient descent cannot be used, so Bayesian Optimization is employed.

Bayes Theorem.

The core idea of Bayesian Optimization is to use existing observational data, combined with Bayes' Theorem, to update the estimation of the objective function is as follows Equation 12: p θ / D = P D / θ P θ P D (12)

Where:

θ represents the parameters of the objective function (i.e., the variables we want to estimate); D = x i , y i i = 1 n is the existing observational data, with y i = f x i + ε being the noisy observations; P θ is the prior probability, representing our prior assumption about the objective function; P D / θ is the likelihood function, representing the likelihood of the data given the parameters; P D is the normalization factor, ensuring that the sum of all probabilities equals 1.

The core idea of Bayesian Optimization is to continuously update the posterior distribution of the objective function after each sampling, and then use the posterior distribution to select the next optimal sampling point.

Gaussian Process (GP) for Modeling the Objective Function.

Since the objective function f x is unknown, we need to represent it using a probabilistic model. Bayesian optimization typically uses a Gaussian process (GP) to model f x .

Definition of Gaussian Process.

The Gaussian process is a distribution over functions, meaning that for any set of input points x 1 , x 2 , . . . . , x n , the corresponding function values f x follow a multivariate Gaussian distribution is as follows Equation 13. f x ∼ G P m x , k x , x ′ (13) m x is the mean function (usually set to 0) is as follows Equation 14: m x = E f x (14) k x , x ′ is the covariance function (or kernel function), which represents the correlation between different input point is as follows Equation 15. k x , x ′ = C o v f x , f x ′ (15)

Commonly used kernel functions include the RBF kernel (Radial Basis Function kernel), Matern kernel, polynomial kernel, and so on. The definition of the RBF kernel is as follows Equation 16: k x , x ′ = σ 2 ⁡ exp − x − x ′ 2 2 l 2 (16) Where:

σ 2 : Controls the variance.

Prediction with Gaussian Processes

Given n observed points D n = x i , y i i = 1 n , we want to predict the distribution of f x * at a new point x * . Assume the existing data y = y , 1 . . . y n , with mean vector and covariance matrix is as follows Equation 17:

μ n = m X = m x 1 , . . . . , m x n

For the new point x * , define are as follows Equations 18, 19: k * = [ k ( x * , x 1 , [ k ( x * , x 2 . . . . [ k ( x * , x n (18) k * * = k x * , x * (19)

The predicted distribution at x * is as follows Equation 20: μ * = k * T K n − 1 y σ * 2 = k * * − k * T K n − 1 k * (20) the new point x * is as follows Equation 21: f x * D n ∼ Ν μ * , σ * 2 (21)

This gives the posterior distribution of f x , which is used to guide the selection of the next sampling point.

Acquisition function.

The acquisition function (AF) is used to determine the next sampling point, aiming to efficiently explore the optimal solution while minimizing the number of trials.

Common acquisition functions include:

Expected Improvement (EI) is as follows Equation 22. E I x = E max ⁡ f x − f x + , 0 (22)

Probability of Improvement, PI is as follows Equation 23 P I x = P f x > f x + + ξ (23)

Where f x + is the current optimal value. Where ξ is the exploration parameter.

Upper Confidence Bound, UCB is as follows Equation 24. U C B x = μ x + κ σ x (24)

Where κ controls the trade-off between exploration and exploitation.

The next step selects the point that maximizes the acquisition function is as follows Equation 25. x n + 1 = arg max A F x x ∈ x (25)

NSGA-II is a genetic algorithm-based multi-objective optimization method introduced by Deb et al., in 2002. It is designed to solve multi-objective optimization problems (MOPs) and can simultaneously optimize multiple conflicting objectives, generating a set of uniformly distributed Pareto optimal solutions. The NSGA-II algorithm mainly involves steps such as nondominated sorting, crowding distance calculation, selection, crossover, and mutation.

1. Definition of Multi-Objective Optimization Problem is as follows Equation 26

f m x represents the m-the objective function, m = 1 , 2 , . . . , M ; The definition of a Pareto optimal solution is as follows: If there exists a solution x * such that:

∃ f m x * ≤ f m x f m x * < f m x and there exists at least one m such that f m x * < f m x , then x * is called a Pareto optimal solution.

2. Non-dominated Sorting

In NSGA-II, individuals are ranked in layers based on the non-dominance relationship. For a given solution a and b, the dominance relationship is defined as follows Equation 27: ∀ m ∈ 1 , 2 , . . . , M , f m a ≤ f m b (27) and there exists at least 1 m such that f m a < f m b , this is called “a dominates b”.

In the non-dominated sorting process, all individuals that are not dominated by any other solution form the first front F 1 is as follows Equation 28. F 1 = x i ∀ x j ∈ S , x j ⊄ x i (28)

After removing F 1 , the Pareto front F 2 of the remaining population is determined by finding all individuals that are not dominated by any solution in the remaining population is as follows Equation 29. F 2 = x i ∀ x j ∈ S − F 1 , x j ⊄ x i (29)

This process is iteratively repeated until all individuals are assigned to a specific level or front.

3. Crowding Distance Calculation

To improve the uniformity of the Pareto solutions' distribution, NSGA-II introduces the crowding distance (Crowding Distance) metrics CCC and DDD to evaluate the sparsity of an individual’s neighborhood.

For the sorted individual iii, its crowding distance in the m objective direction is calculated is as follows Equation 30: C D i m = f m i + 1 − f m i − 1 f m max − f m min (30) f m i + 1 and f m i − 1 are the function values of the neighboring individuals in the m -the objective dimension.

f m max and f m min are the maximum and minimum values of the objective in that dimension.

The total crowding distance calculation is as follows Equation 31: C D i = ∑ m = 1 M C D i m (31)

The crowding distance of boundary solutions is set to infinity is as follows Equation 32. C D min = ∞ , C D max = ∞ (32)

4. Selection mechanism

The individual with the lower level is given priority; if the levels are the same, the individual with higher congestion is chosen.

5. Crossover and Mutation Operations

Simulated Binary Crossover:

Crossover between parent individuals x i and x j to produce offspring are as follows Equations 34, 35: x i ′ = 1 + β x i + 1 − β x j 2 (34) x j ′ = 1 − β x i + 1 + β x j 2 (35)

Where,

β = 2 r 1 η + 1 if r ≤ 0.5

β = 2 1 − r − 1 η + 1 if r > 0.5

Polynomial Mutation are as follows Equations 36, 37 x i ′ = x i + Δ i x max − x min (36) Δ i 2 r 1 η + ! r < 0.5 2 1 − r − 1 η + 1 r ≥ 0.5 (37)

6. Population Update Mechanism

Q t is the new population generated by selection, crossover, and mutation with size N; after merging, R will have a size of 2N.

2. Perform Non-Dominated Sorting

Perform fast non-dominated sorting on R t and assign individuals to different Pareto layers is as follows Equation 39. F 1 , F 2 , . . . F k (39)

Perform non-dominated sorting on R t and select the top N individuals based on hierarchy. F 1 is the Pareto-optimal front (not dominated by any individual); F 2 is the next-best front (only dominated by individuals in F1); and so on, forming multiple ranks.

3. Select the top N individuals is as follows Equation 40

If a certain front F_k cannot fully fill P t + 1 , select the best solutions based on crowding distance ranking.

The safety factor obtained by model the predicted slope step parameters using 3DEC is shown in the As shown in Figure 16a the result calculated by SVM-NSGAII is 1.211, Figure 16b, the result obtained from the SVM-BO predicted parameters is 1.254. In Table 9 The error for SVM-NSGAII is 4.0%, and the error for SVM-BO is 6.0%. By comparing the results with machine learning metrics such as MAE, MSE, MAPE, and R², it is evident that the optimization result obtained using SVM-NSGAII is more accurate.

Based on the relative error verification calculation shown in Table 10, the error rate for the SVM-NSGAII algorithm is 3.8%, while the error rate for the SVM-BO algorithm is 5.6%. From the analysis, it can be concluded that using SVM-NSGAII is more accurate.