Predicting the Thickness of an Excavation Damaged Zone around the Roadway Using the DA-RF Hybrid Model

1  Introduction

After the excavation of the roadway, the original stress equilibrium state of the surrounding rock is broken, resulting in stress redistribution. On the one hand, the radial stress decreases and the tangential stress increases; on the other hand, the stress state of the surrounding rock changes from a three-dimensional stress state to an approximately two-dimensional stress state. Therefore, the surrounding rock will appear strength reduction and stress concentration phenomenon. When the induced stress in the surrounding rock exceeds the strength of the rock mass, the rock mass will be broken, and a ring-shaped broken rock zone (BRZ) will be formed within a certain range around the roadway, which is called the excavation damaged zone (EDZ) or loose circle [1–3]. The thickness of the excavation damaged zone is not only an important basis for judging the stability of the roadway but also crucial for the construction of underground space engineering and the design of supporting structures [4–6]. Therefore, accurate prediction of the thickness of the excavation damaged zone (EDZ) around the roadway is of great importance in engineering practice.

At present, there are several methods for predicting EDZ thickness, including the empirical formula method, on-site measurement method, numerical simulation, and machine learning (ML). Based on onsite tests and theoretical analysis, various empirical formulas for calculating EDZ thickness have been proposed [7–10]. Shemyakin et al. [11,12] proposed the concept of a discontinuous zone and derived an empirical formula for the thickness of the discontinuous zone. Serrano et al. [13] obtained a mathematical expression for determining the radius of the loose circle based on the convergence-confinement method. Wang et al. [14] revised the H-B strength criterion and derived a theoretical formula for the loose circle radius. Qiao et al. [15] used the Drucker-Prager yield criterion to derive a formula for calculating the thickness of EDZ. Although these empirical formulas can predict the thickness of EDZ, they have been validated for only a few engineering sites and are difficult to be applied beyond the original environmental features, and are not applicable to complex engineering environments. With the development of technology, there are a variety of on-site measurement methods to determine the thickness of EDZ, including the ultrasonic method [16,17], geological radar method [18], microseismic monitoring [1,19], complex resistivity method [20], borehole imaging method [21–23], multi-point extensometer method [24], etc. However, the on-site measurement method has disadvantages such as cumbersome operation, high measurement cost, and is easily restricted by on-site conditions. Based on the development of computer technology, numerical simulation methods have been widely used to study the distribution of EDZ around roadways [25–30]. Golshani et al. [31] performed numerical simulations of the EDZ around openings in brittle rocks. Xia et al. [32] used the finite element method to obtain the loose range of the surrounding rock of the roadway through numerical simulation, and the simulation results were very close to the actual results. Yang et al. [33] performed a 3D numerical simulation of the evolution of the EDZ around the tunnel using LS-DYNA software. An et al. [34] proposed a hybrid finite element-discrete element method to simulate the EDZ formation process around the tunnel. Some scholars have used finite elements, discrete elements, marker point processing techniques, and material point dormancy method to build numerical models for estimating EDZ and obtained the formation pattern and distribution characteristics of the loose zone [35–38]. Although the numerical simulation method is low-cost and easy to operate, there are many assumptions in the simulation process, and it is often necessary to simplify the boundary conditions and material properties, resulting in model results that are inconsistent with reality.

Due to the advantages of ML methods in dealing with nonlinear problems, ML has been widely used to solve complex problems in the engineering field [39–48]. In particular, there have been many studies using ML to predict the thickness of EDZ and its related problems. Zhou et al. [3] developed a support vector machine (SVM) model for estimating the broken rock zone thickness (BRZT), and its predictive performance was more accurate than the regression model. Liu et al. [49] established a roadway surrounding rock stability prediction model based on concept lattice theory and improved probabilistic neural network. Xie et al. [2] developed a random forest regression model to predict EDZ, and the results proved that the random forest model has excellent prediction performance. Yu et al. [50] developed four prediction models of the multiple linear regression, artificial neural network, Gaussian process (GP), and particle swarm optimization (PSO)-GP to predict EDZ thickness. Zhao et al. [51] used the sparrow search algorithm (SSA) to optimize the backpropagation neural network (BPNN), Elman neural network (ENN), and support vector regression (SVR), and constructed SSA-BPNN, SSA-ENN, and SSA -SVR three models to predict the thickness of EDZ around the roadway, of which the SSA-BPNN model performed the best. According to the above, ML has the advantages of low cost, less time, high reliability of results, and shows excellent prediction performance in predicting the thickness of EDZ around the roadway.

Random forest (RF), as an efficient machine learning (ML) technique, has been applied in many engineering practices and has shown good predictive ability [52–57]. Dragonfly algorithm (DA) is a new type of intelligent optimization algorithm, which has the characteristics of strong optimization ability, simple principle, and few parameters [58]. Based on the above, this paper constructs a DA-RF hybrid model to predict the thickness of EDZ around the roadway. The model uses the dragonfly algorithm to optimize two hyperparameters in RF, namely the number of features used to grow each tree (mtry) and the number of trees (ntree), to enhance its predictive power. In addition, the DA-RF hybrid model is compared with three classical models to evaluate its predictive performance. Finally, the degree of influence of each output parameter on the prediction of EDZ thickness was determined by sensitivity analysis.

2  Methods

2.1 Dragonfly Algorithm (DA)

The dragonfly algorithm (DA) is a metaheuristic algorithm proposed by Mirjalili [58] based on the hunting and migration behavior of dragonflies, which has been applied in many fields [59–66]. The algorithm simulates two swarming behaviors of dragonflies, namely static swarm and dynamic swarm. In static swarming behaviors, dragonflies divide into small groups and fly back and forth in a small area in search of prey, which is similar to the exploration phase. In dynamic swarming behaviors, a large number of dragonflies form swarms for long-distance migration in the same direction, which is similar to the development stage [67,68]. During the movement of dragonflies, the updating of individual positions is related to five behaviors, namely separation, alignment, cohesion, foraging, and avoiding the enemy, as shown in Fig. 1. The mathematical expressions of the five behaviors are as follows.

Figure 1: Five behavioral patterns of the DA

Separation refers to the behavior of dragonflies to separate from the population in order to avoid static collisions with neighbors during flight, its mathematical model is as follows:

Si=−∑j=1NX−Xj(1)

where X is the position of the current individual, Xj is the position of the j-th neighboring individual, and N is the number of neighboring individuals in the population.

Alignment refers to the speed at which a dragonfly individual maintains the same speed as its neighbors when flying, and its mathematical model is as follows:

Ai=∑j=1NVjN(2)

where Vj denotes the flight speed of the j-th neighboring individual.

Cohesion refers to the behavior that dragonfly individuals tend to gather toward the center of neighboring individuals, and its mathematical model is as follows:

Ci=∑j=1NXjN−X(3)

where X represents the current position of the dragonfly individual.

Foraging is the behavior of a dragonfly that is attracted to food and flies towards it, and its mathematical model is as follows:

Fi=X+−X(4)

where X+ is the location of the food source.

Avoiding the enemy refers to distracting the enemy's attention and escaping from the natural enemy, and its mathematical model is as follows:

Ei=X−−X(5)

where X− represents the position of the enemy.

The position of the food source is the current optimal position of the algorithm, and the fitness value is the best, and the position of the enemy is the current worst position, and the fitness value is the worst. DA is developed based on the particle swarm algorithm, and to imitate the flight direction of the dragonfly in space as well as the step length, a step vector ΔX and a position vector X were introduced. The step vector represents the direction of motion of the dragonfly and is mathematically defined as follows:

ΔXt+1=(sSi+aAi+cCi+fFi+eEi)+wΔXt(6)

where t denotes the number of current iterations, s is the separation weight, a is the alignment weight, c is the cohesion weight, f is the foraging weight, e is the enemy avoidance weight, w is the inertia weight, and ΔXt+1 is the next population location update step.

The position vector represents the position of the individual dragonfly movement. When the former individual has neighboring dragonflies (i.e., there are other dragonflies within his radius (r)), the mathematical model is as follows:

Xt+1=Xt+ΔXt+1(7)

When there is no neighboring dragonfly, to improve the randomness of dragonfly searching, the Lévy flight method is used to fly around the search space to update the position of the dragonfly. The mathematical model is as follows:

Xt+1=Xt+Levy(d)Xt(8)

Levy(x)=0.01×r1×σ|r2|1/β(9)

σ=(Γ(1+β)×sin⁡(πβ2)Γ(1+β2)×β×2(β−12))1β,Γ(x)=(x−1)!(10)

where d is the position vector dimension, r1 and r2 are random numbers in [0,1], and β is a constant.

To adjust the search performance of the DA algorithm, the parameter values (s, a, c, f, e) will be adaptively adjusted during the search process. The algorithm flow of DA is shown in Fig. 2.

Figure 2: The flow chart of DA

2.2 Random Forest (RF) Algorithm

Random Forest (RF) algorithm is a machine learning model based on the decision tree proposed by Breiman [69]. It essentially uses the idea of ensemble learning to integrate multiple decision trees by bagging methods (also known as bootstrap aggregating) and random feature selection to increase the diversity of classification and regression trees (CART) [70], the mean value of each regression decision tree is taken as the prediction result. The RF algorithm has good generalization ability, is simple to implement, has high precision, and is not prone to overfitting. Fig. 3 shows the structure of the random forest. The basic steps to construct the RF are as follows.

Figure 3: Structure of RF algorithm

First, the bagging method is used to repeatedly randomly select k samples from the original training sample set N with replacement to form k training subsets with a sample size of n. Second, a single CART regression tree (decision tree) is constructed with the training subset. In the training process, the splitting rule for each node is to randomly select m (m ≪ M) feature attributes from the set M of all feature attributes, then select the optimal splitting point from these m features, and then divide the left and right subtrees until it cannot be split, and there is no pruning in the formation of the regression decision tree. Repeat the above steps to construct a large number of regression decision trees and maximize the growth of each tree, and these regression decision trees will form a random forest. Finally, based on the idea of ensemble learning, the prediction results of all regression decision trees are averaged as the final prediction results [71–73].

3  Materials

3.1 Data Preparation and Description

To construct a DA-RF hybrid model for predicting EDZ thickness, we built a database by referring to previous literature [51,74,75]. The database has collected a total of 217 sets of data from various mines in different places. In this database, each sample contains embedding depth (ED), drift span (DS), surrounding rock mass strength (RMS), joint index (JI), and the excavation damaged zone thickness (EDZT), where ED, DS, RMS, and JI are used as input parameters and EDZT is used as an output parameter. ED refers to the depth of the roadway from the ground. Since the thickness of the EDZ increases with the increase of the original rock stress, ED is selected to describe the original rock stress. DS refers to the width of the roadway, the size and shape of the roadway section have a certain influence on the thickness of the EDZ, and the roadway span is generally selected as the analysis index. RMS refers to the uniaxial compressive strength of the surrounding rock, and JI refers to the degree of joint development in the surrounding rock. The strength of the surrounding rock reflects the difficulty of the surrounding rock resisting damage and is inversely proportional to the thickness of the EDZ. In addition, the presence of joints and fissures in the rock mass also affects the strength of the rock mass, so RMS and JI are chosen to jointly describe the strength of the rock mass. These indicators are the main influencing factors of EDZ formation, are easily accessible and the information described is independent [51], so these four parameters were chosen as inputs. Pearson correlation coefficients were calculated to describe the correlation between the input parameters and between the input and output parameters, as shown in Fig. 4. As can be seen from the figure, the correlations among the parameters are low, indicating that the input variables are relatively independent, and the input and output parameters are not simply linearly related to each other. For better understanding, Table 1 shows some statistical information about the dataset, and a schematic of the structure of the EDZ around the roadway is shown in Fig. 5.

Figure 4: Correlation plot of EDZ dataset

Figure 5: Schematic diagram of the structure of EDZ around the roadway

3.2 Model Verification and Evaluation

In order to verify and evaluate the prediction performance of the DA-RF hybrid model and ensure that the model can reliably predict the thickness of the EDZ around the roadway, four evaluation indicators were selected in this study, namely mean absolute error (MAE) and root mean square error (RMSE), determination coefficient (R2), and variance accounted for (VAF) [76–78]. The smaller the RMSE and MAE values, the closer the R2 is to 1, and the closer the VAF value is to 100, indicating that the prediction quality of the model is higher. The formulas of the evaluation indexes are as follows [79–85]:

MAE=1N∑i=1N|Yi−Yi^|(11)

RMSE=1N∑i=1N(Yi−Yi^)2(12)

R2=1−∑iN(Yi−Yi^)2∑iN(Yi−Yi¯)2(13)

VAF=[1−var(Yi−Yi^)var(Yi)]×100(14)

where N is the number of samples, Yi is the actual value of EDZ thickness, Yi^ is the predicted value of EDZ thickness, and Yi¯ is the average value of EDZ thickness.

4  Results and Discussion

4.1 DA-RF Prediction Model

In this paper, the RF algorithm is used to predict the thickness of the EDZ around the roadway, and DA is used to optimize the two hyperparameters mtry and ntree in RF to improve the prediction performance of the RF model. From this, a DA-RF hybrid prediction model is constructed. To more intuitively illustrate the architecture of the proposed hybrid prediction model, the construction process of DA-RF model development and evaluation is shown in Fig. 6. As can be seen from Fig. 6, the model can be divided into four stages. First, the dataset was randomly divided into the training set and the test set, of which the training set accounts for 80% and the test set accounts for 20%. After that, all datasets were normalized to the range [−1, 1], which was to attenuate the influence of input parameters of different scales on the results and reduce the computational cost. Next, initialize the relevant parameters, optimize the two hyperparameters mtry and ntree in RF through DA, and the RMSE was used as the fitness value to determine whether the termination condition was satisfied and if so, the best RF model was determined. Finally, the prediction performance of the proposed hybrid models was evaluated by calculating the evaluation indicators (R2, VAF, RMSE, MAE) for the training and test sets.

Figure 6: Flowchart of hybrid model based on RF optimized by DA technique

The correlation between the predicted and actual values of the training and test sets of the DA-RF model is shown in Fig. 7. The confidence and prediction bands in the figure define the upper and lower bounds and width of the correlation interval, with the 95% confidence band indicating a 95% probability of containing the true best-fit linear regression line within the interval and the 95% prediction band indicating the region expected to contain 95% of the new data points. As can be seen from the figure, the confidence band and prediction band of both the training and test sets are narrow, indicating that the prediction performance of the DA-RF hybrid model is relatively stable. It is worth noting that several hybrid models have been developed to predict the thickness of the EDZ. Zhao et al. [51] developed an SSA-BPNN model with an MAE, RMSE, and R2 of 0.1246,0.1636, and 0.9277, respectively, during the test phase. Yu et al [50] constructed the PSO-GP with RMSE, R2, and VAF of 0.1591, 0.9153, and 89.2192, respectively, during the testing phase. The MAE, RMSE, R2, and VAF of the DA-RF model constructed in this paper were 0.1115, 0.1417, 0.9423, and 94.0836, respectively, in the testing phase, which was significantly better than the above two models and possessed higher prediction performance.

Figure 7: Predictive performance of DA-RF model

4.2 Comparison with Classical Models

Several classical models and risk assessment procedures have been applied to high-risk operations and projects such as blasting and quarrying [86–90]. To verify and evaluate the superiority of the DA-RF model prediction performance, we introduced four classical models to predict the thickness of EDZ, which were back propagation neural network (BPNN), extreme learning machine (ELM), radial basis function network (RBF) and RF, and the prediction results of these classical models were compared with the DA-RF model, as shown in Table 2. In the table, we comprehensively score the evaluation indicators of the training set and test set of the models. The higher the score, the better the prediction performance of the model [91,92].

The prediction results of the test set of DA-RF and classic models are shown in Fig. 8. The larger the area of the bubble in the figure represents the larger the absolute error (AE), and the AE is calculated as follows:

AE=|Y−Y^|(15)

Figure 8: Prediction error analysis of the five models test sets

As can be seen from Fig. 6, the AE value of the DA-RF model is generally lower than that of the other four models. To compare the predictive performance of all models more intuitively and accurately, we rank the comprehensive scores of the five models, and the results are shown in Fig. 9. It can be seen that the DA-RF model has the best performance in predicting the thickness of the EDZ.

Figure 9: Comprehensive score ranking of EDZ thickness prediction models

4.3 Sensitivity Analysis

The prediction of the thickness of the EDZ around the roadway is related to the roadway excavation and the design of the support structure, which plays an important role in ensuring the stability of the roadway. To accurately predict the EDZ thickness, the influence of various factors must be comprehensively considered. It is known that the four input parameters selected in this paper have a certain influence on the prediction of the thickness of the EDZ, but the extent of their influence is not clear. Therefore, sensitivity analysis is needed to judge the importance of each input variable in the DA-RF model to the output EDZ thickness.

The RF algorithm has built-in feature importance to calculate the sensitivity of the input parameters [69,93]. The basic idea is to calculate the mean square error (MSE) of each decision tree in the RF using the corresponding out-of-bag (OOB) data sample for it. Then keep the values at other features unchanged, perform random up-down permutation on the value at feature j in OOB, and the MSE of its OOB is calculated again. The importance of feature j can be represented by the difference between the MSEs calculated before and after [73]. Assuming there are N trees, the importance of feature X can be calculated by the formula (16).

Importancej=1N∑n=1N(MSE2nj−MSE1nj)(16)

where MSE1nj and MSE2nj denote the MSE of the OOB before and after the alteration. if feature j is randomly perturbed and has a significant effect on the error of the OOB, it means that this feature is important. The relative importance of the calculated input parameters is shown in Fig. 10, where the most sensitive parameter is the joint index (JI).

Figure 10: Relative importance of input variables in predicting EDZ thickness

5  Conclusions

The thickness of the excavation damaged zone (EDZ) is an important basis for determining and controlling the stability of roadway engineering construction, and it is of great significance to accurately predict the thickness of the EDZ in engineering practice. Therefore, in this paper, a hybrid RF-based prediction model is developed to accurately predict the thickness of the EDZ around the roadway. The model utilized the dragonfly algorithm (DA) to optimize two hyperparameters in RF, namely mtry and ntree. We established a database containing 217 sets of data, of which 174 sets were used as the training set and the remaining 43 sets as the test set. The embedding depth (ED), drift span (DS), surrounding rock mass strength (RMS), and joint index (JI) were used as input variables, and the excavation damage zone thickness (EDZT) was used as output variables. To demonstrate the superiority of the DA-RF model, four models, BPNN, ELM, RBF, and RF, are cited for comparison with the DA-RF model. The results show that the DA-RF model is significantly better than the unoptimized classical model and has the best prediction performance (training set: MAE = 0.1036, RMSE = 0.1514, R2 = 0.9577, VAF = 94.2645; test set: MAE = 0.1115, RMSE = 0.1417, R2 = 0.9423, VAF = 94.0836), so it is recommended to predict the thickness of the EDZ around the roadway. Finally, the relative importance of each input variable was determined by sensitivity analysis, and the order from low to high was DS, ED, RMS, and JI. The hybrid model developed in this paper provides a simple and accurate method for predicting the thickness of the excavation damaged zone, and the algorithm used can be applied to other similar engineering problems. In the future, it is necessary to collect more data samples, consider more influencing factors, and optimize the algorithm to further improve the prediction performance of the hybrid model.

Funding Statement: This research was funded by the National Science Foundation of China (42177164), the Distinguished Youth Science Foundation of Hunan Province of China (2022JJ10073) and the Innovation-Driven Project of Central South University (2020CX040).

Conflicts of Interest: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.