Abstract
The efficiency and accuracy in determining mechanical parameters of joint of rock mass have significant influences on the safety of construction. In this paper, a new back-analysis method of joint parameters of the ubiquitous-joint model is proposed. This study combines the ubiquitous-joint model, random forest algorithm (RF), and cuckoo search algorithm (CS) to construct the parameters identification method of a jointed rock mass. The parameter determination is transformed into a global optimization problem. In this method, the mean square error between the measured displacement (stress) and the calculated displacement (stress) is taken as the objective function, and the joint parameters are taken as the decision variables. The method is applied to the project of Qingniwaqiao Metro Station of Dalian Metro Line 5. This method can be used for the back-analysis of parameters of a jointed rock mass in underground engineering construction.
1. Introduction
Accurate determination of rock constitutive model and model parameters is the key to geotechnical engineering design and application. To determine the mechanical parameters of rock mass, the back-analysis method is proposed based on field measurement information. Using different monitoring information, the back-analysis method can be divided into stress back-analysis method, displacement back-analysis method, and mixed back-analysis method. Among these approaches, the mixed back-analysis method combines the information of stress and displacement, so the results from the back-analysis are more accurate.
In the process of back-analysis, different constitutive models reflect different mechanical behaviors. The late 1970s to early 1980s is the starting period of back-analysis. The research object of this period is a two-dimensional elastic linear medium [1, 2]. With the development of computing science and computer technology, the research object has changed into three-dimensional media, nonlinear or elastic-plastic media, and viscoelastic media. Zhang et al. [3] proposed a displacement back-analysis method for underground rock mass engineering considering displacement loss. Gao et al. [4] investigated the impact of rock mass in a time-dependent manner using the displacement back-analysis method. Zuo et al. [5] proposed a reliability back-analysis method to determine these shear strength parameters of landslide based on nonlinear Mohr-Coulomb (M-C) failure criterion. The previously mentioned back-analysis methods are mostly limited to anisotropic homogeneous model, and it is difficult to accurately describe the mechanical properties of the jointed rock mass. During the excavation of tunnels, the existence of joint of rock mass is one of the main factors affecting the deformation of supporting structures [6, 7]. Ubiquitous-joint (U-J) model can consider the properties of joints and rock blocks at the same time, and it can better reflect the characteristics of jointed rock and engineering practice [8–10]. There are fruitful results obtained by analyzing tunnel surrounding rock parameters, as a comparison, the ubiquitous-joint model is rarely studied. It is majorly ascribed to the time-consuming process of complex constitutive modelling within numerical back-analysis method and theoretical methods, which renders a more difficult process of parameter acquisition [11, 12]. At present, some mature intelligent algorithms have been widely used in data analysis in various fields. Many scholars have studied the application of intelligent algorithms in geotechnical engineering and made good progress [13–16]. For an improved calculation speed, parameters back-analysis based on intelligent algorithm has become an important development direction in the research of geotechnical engineering back-analysis method. Feng et al. [17] proposed using a neural network to predict the deformation value of surrounding rock instead of finite element calculation. Gao et al. [18] and Sun et al. [19] combined particle swarm optimization (PSO) algorithm with a numerical method to carry out back-analysis of tunnel rock mass mechanical parameters. The intelligent algorithm has been widely used in geotechnical engineering. However, the convergence speed and prediction accuracy of the neural network are easily affected by the initial weights and thresholds. In addition, the PSO algorithm also has the problem of premature convergence, which can not guarantee the global optimal solution. Therefore, it is necessary to explore a better intelligent algorithm for parameters back-analysis of ubiquitous-joint model.
Random forest algorithm (RF) is a new machine learning method which combines ensemble learning and decision tree [20–24]. RF has the advantages of few parameters, high training efficiency, and not easy to overfit. It has a good effect on solving multivariable mapping and is honored as one of the best machine learning algorithms at present. Cuckoo search algorithm (CS) is simple and easy to implement, with few parameters, strong search ability and fast convergence speed, which has obvious advantages in solving optimization problems [25, 26].
Therefore, in this paper, RF and CS are applied to the parameter identification of U-J model. The three-dimensional numerical simulation, combined with the orthogonal design method, was used to generate training data samples. Besides, the joint parameters were correlated with the displacement and stress of surrounding rock with the help of RF. Moreover, the CS was used to optimize the parameters of the RF model hence to improve the mapping effect. Subsequently, the CS algorithm integrated with well-trained RF model was used to determine the parameters of joints. Additionally, the previously mentioned method is successfully applied to the construction of Qingniwaqiao Metro Station of Dalian Metro Line 5, and it proves valid and quite useful for parameters back-analysis of joints. Finally, the prediction accuracy of RF and the selection process of parameters of RF and CS are analyzed.
2. The Back-Analysis of Parameters Based on RF-CS
2.1. U-J Model
The U-J model is an extension of the M-C model; that is, a weak joint surface is added to the M-C model, and the weak joint surface also obeys the M-C yield criterion. The stress of rock mass is the sum of the stress calculated by M-C model and the local stress calculated by universal joint model converted to the global coordinate system. The failure may first occur in rock mass or joint weak plane. The joint parameters include surface cohesion (Jcoh), surface friction angle (Jfric), surface tensile strength (Jten), surface dilatation angle (Jdila), surface dip (Jdip), and surface dip angle (Jdd). The U-J model can simulate the influence of joints on the stability of caverns without adding extra finite element mesh, so it has good operability and high computational efficiency in the stability analysis of caverns.
2.2. RF Regression Algorithm (RF)
The RF is a fusion algorithm based on decision tree classifier. Its basic idea is based on statistical theory, using bootstrap resampling method to extract multiple samples from the original samples, building a decision tree for each bootstrap sample, and then taking the average value of all decision trees as the final prediction result. RF can be regarded as a strong predictor integrated by many weak predictors (decision trees). Set the training set to be extracted from the independent random vector (X, Y); the input vector is X and the output vector is Y. Then, the mean square generalization error of the predicted output h(X) is
The prediction output of RF is obtained by taking the average value of k decision trees {h(θ, Xk)}, which has the following definition:(1)when k ⟶∞, Let the right side of expression [2] be PE, which is the generalization error of RF. The average generalization error PE of each decision tree can be defined as(2)For all θ,where is the weighted correlation coefficient of residuals and , and and are independent of each other.
RF reduces the average error of decision tree by .
The steps of RF can be summarized as follows.
Set up θ as a random parameter vector, and the corresponding decision tree is T(θ). Note that B is the field of X, that is, , where p ∈ N is the dimension of the independent variable. Each leaf node of the decision tree corresponds to a rectangular space of B. Note that the rectangular space of each leaf node is (l = 1,2,…,L). For every x ∈ B, if and only if a leaf node l satisfies x ∈ Rl, let the leaf node of the decision tree T(θ) be l (x,θ).
(1)The bootstrap method is used to resample and randomly generate K training sets θ1,θ2,…,θk; the corresponding decision trees generated by each training set are {T(x,θ1)},{T(x, θ2) } ,…,{T(x, θk)}.(2)Assuming that the feature has m dimension, m features are randomly selected from the M dimension features as the split feature set of the current node, and the node is split in the best way among the m features.(3)Each decision tree can grow to the maximum without pruning.(4)For the new data, the prediction of a single decision tree T(θ) can be obtained by averaging the observed values of leaf nodes l(x, θ). If an observed value Xi belongs to the leaf node l(x, θ) and is not 0, let the weight ωi(x, θ) be where the sum of weights is equal to 1.(5)The prediction of a single decision tree is based on the observed values Yi (i = 1,2, ..., n) of the dependent variable and is obtained by the weighted average. The predicted value of a single decision tree can be obtained by the following formula: )For decision tree weights ω i(x, θ t) (t = 1,2,…,k), take the average value to obtain each observed value Yi∈(1, 2…,n) weight ω i(x):
The predicted value of RF is as follows:
2.3. The Process of RF Optimized by CS
There are two important parameters of RF: one is the number of preselected variables of tree nodes (Mtry), and the other is the number of random forest trees (Ntree). The accuracy of RF predictions is affected by these two parameters. If Mtry is too small, the classifier will be overfitted and the accuracy of prediction classification will be reduced, and if it is too large, it will lead to the slowdown of running speed. With the decrease of Mtry, the correlation between trees also decreases gradually, and the classification accuracy also decreases. If Ntree is too small, it will lead to insufficient training, and if it is too large, it will increase the amount of model computation. Because it is difficult to select these parameters manually, this paper uses intelligent optimization algorithm CS to select the optimal parameters of RF.
The steps of optimizing RF parameters by CS algorithm are as follows.(1)Set the number of Mtry and Ntree initially, and determine the value range of parameters Mtry and Ntree.(2)A population with n nests is initialized, including the number of nests n, the probability of being found by the nest owner Pa, and the number of iterations K. The position of nests is generated randomly, and each nest position represents a pair of parameters (Mtry, Ntree).(3)Using the mean square error () as the fitness value, the fitness value of each nest is calculated, and the optimal nest position is reserved for the next generation. The definition of mean square error is as follows: where yi is the calculated value; is the predicted value; n is the number of prediction samples.(4)Follow Levy’s flight process to update the nest: where and are the positions of bird nests in the t+1 and t iterations; α is the scale factor of step size, α > 0; s is the step length; ⊕ is point-to-point multiplication. L(s, λ) is the path of random walk and obeys the exponential distribution of u = t − λ, where λ is a constant, λ∈(1,3].(5)Calculate the fitness of the updated nest, compare the fitness with that of the previous generation nest, and update the optimal position of the nest. Judge whether the termination condition is met. If so, stop the iteration. Otherwise, return to step 3.(6)According to the optimal nest position, the corresponding Mtry and Ntree are output and assigned to RF model. Thus, the optimal parameters of RF algorithm are selected.
2.4. Back-Analysis Process of Joint Parameters
The joint parameters were correlated with the displacement and stress of surrounding rock with the help of RF. Subsequently, the CS algorithm integrated with well-trained RF model was used to determine the parameters of joints.
The back-analysis calculation steps of joint parameters are listed as follows.(1)With the orthogonal design and uniform design methods, the different combination schemes of six parameters of jointed rock mass (Jcoh, Jfrc, Jten, Jdila, Jdip, and Jdd) are designed. The constructed training data sets are corresponsive to the joint parameters, the displacement, and stress of key points.(2)The RF model is used to validate the nonlinear responses of displacement and stress against the joint parameters. The orthogonal design scheme is employed as the training sample, and the uniform design scheme is employed as the prediction test sample.(3)According to the condition of training samples, set initial parameters, including the number of nests n, the probability of being found by the nest owner Pa, the number of iterations K, and the range of parameters Mtry and Ntree.(4)CS algorithm uses the mean square error as the fitness value and uses Levy’s flight process to update the nest.(5)If the number of iterations is met, the optimal parameters will be output assigned to RF model, and then proceed to the next step. Otherwise, step (4) will be returned.(6)The parameters of RF-CS are initialized, and the initial population in the solution space is a group of randomly generated parameters.(7)The nest was updated following the mean square error and the Levy’s flight process. A set of parameters is input into the above RF for displacement and stress measurements of the key points. Equation (9) is utilized to evaluate the fitness of the current individual, and the nest replacement operation is adopted.(8)If the number of iterations or the minimum error requirements were matched, the optimal parameters of jointed rock mass will be output and the calculation will end; otherwise, return to step 7.
Figure 1 shows the process of back-analysis of joint parameters.

3. Engineering Application
3.1. Engineering Overview
Qingniwaqiao Metro Station is the station of Dalian Metro Line 5, which is arranged along Youhao street from north to south. The terrain in the site fluctuates greatly. The rock mass where the main body of the station is located has developed joints and fissures. The location and geological conditions of the station are shown in Figure 2.

3.2. Numerical Simulation Model
The geotechnical engineering software FLAC3D is used to simulate the construction process of the station. According to the excavation size of the station, the height of the calculation model is 98 m, the width is 212 m, and the longitudinal length is 20 m (Figure 3(a)). The representative geological sections in the construction site are selected. The main structure of the station is in the moderately weathered slate. The U-J model is selected as constitutive model. This paper mainly studies the excavation of pilot tunnels. The excavation sequence of pilot tunnel is ①, ②, ③, and ④. The layout of monitoring points in the excavation stage of pilot tunnel ③ is shown in Figure 3(b).

(a)

(b)
The surrounding rock parameters of moderately weathered slate are listed: modulus of elasticity is 1300 MPa, Poisson’s ratio is determined at 0.24, the value of cohesion is 3.7 MPa, and internal friction angle is 38°. The range of joint parameters is determined as shown in Table 1.
25 orthogonal design schemes and 5 uniform design schemes are formed, respectively. The corresponding parameters of each test scheme are brought into FLAC3D for numerical calculation. The displacement and stress values of the supporting structure during the excavation of pilot tunnel ③ are recorded. The excavation footage of pilot tunnel ③ is 1.0 m. The vault settlement is expressed as AZ and CZ, respectively, the bottom uplift is BZ, and the horizontal convergence is DZ. Use the data in Table 2 as training samples and the data in Table 3 as test samples.
To verify the influence of selected joint parameters on the stability of rock mass around the station, the values of other parameters were controlled as fixed values, and select Jdd as 146° and 186° and Jdip as 19° and 27° for analysis. Figure 4 shows the distribution of plastic zones with different Jdd and Jdip.

(a)

(b)

(c)

(d)
It can be seen that the change of Jdd and Jdip has certain influence on the stability of surrounding rock around the station. The plastic zone distribution of rock mass caused by station construction is different with different Jdd and Jdip. The influence of the existence of joints on the stability of surrounding rock during construction cannot be ignored.
3.3. Back-Analysis of Joint Parameters
After the excavation of pilot tunnel 3, the joint parameters are back-analyzed according to the field monitoring results. Set the maximum number of iterations K = 100 of CS, the probability of nest discovery Pa = 0.25, and the number of nests N = 50. The field monitoring data were as follows: AZ was 14.34 mm, BZ was 7.54 mm, CZ was 11.24 mm, DZ was 5.12 mm, P1 was 146.36 kPa, and P2 was 124.53 kpa. Table 4 shows the results of back-analysis.
Using the joint parameters obtained from back-analysis, the excavation simulation analysis of pilot tunnel ③ is carried out, and the monitoring values of each monitoring point are obtained, as shown in Table 5. Compared with the field monitoring data, the largest relative error is determined at only 5.78%, verifying the accuracy and feasibility of the back-analysis method.
Figure 5 shows the fitness values of RF-CS algorithm population in the 1st, 5th, 15th, and 20th iterations. The results show that, with the increase of evolution times, the distribution of solution vector gradually decreases in space, exhibiting a converging tendency. For instance, a relatively scattered population was observed in the 1st evolution, and the overall fitness value was close to 1.6. In contrast, the overall fitness decreased in the 5th evolution, with an overall fitness value of about 0.6. It further decreased to 0.25 in the 15th evolution. Finally, the overall fitness value was less than 0.02 in the 20th evolution.

(a)

(b)

(c)

(d)
Figure 6 shows the change of the six parameters searched by the RF-CS as a function of the number of iterations. It can be seen that, upon initial iteration, these parameters are relatively scattered away from the optimal solution and fluctuate. With increasing iteration times, when the number of iterations reaches 22, these parameters tend to be the optimal solution and keep stable. The results show that when the number of iterations reaches 22 times, the optimal solution is obtained, which shows that the optimization process has good convergence.

3.4. Optimization Analysis of Excavation Footage Based on the Back-Analysis of Parameters
In the construction process, the selection of reasonable excavation footage of pilot tunnel can effectively improve the construction efficiency. To seek the optimal scheme in line with the project, the joint parameters obtained by back-analysis are used to simulate the construction process. The construction simulation is carried out with different cyclic excavation footage of pilot tunnel ④.
Combined with the actual construction situation, four kinds of cyclic footage schemes of 0.5 m, 1.0 m, 1.5 m, and 2.0 m are selected for three-dimensional numerical simulation. During the excavation process, pilot tunnel ④ is excavated after pilot tunnel ③ is excavated.
The variation trend of the LDP curve of the vault under the four working conditions is similar (Figure 7). When the footage is selected to be 0.5–1.5 m for excavation, the maximum difference of the settlement of the vault is only 0.3 mm. However, when the excavation footage of 2.0 m is selected, the settlement increases greatly, reaching 9.63 mm, which is 1.89 mm higher than that in the condition of the excavation footage of 0.5 m. The average excavation rate of pilot tunnel is 4.5 m/d, and the maximum settlement rates of vault under four excavation footage are 1.65 mm/d, 1.77 mm/d, 1.81 mm/d, and 2.16 mm/d. During the construction, the monitoring control value of surface settlement rate is 2.0 mm/d, and when the cyclic excavation footage is 2.0 m, the settlement rate has exceeded the monitoring control value. The excavation footage shall be controlled within 2.0 m.

Figure 8 shows the distribution of plastic zone of pilot tunnel under different excavation footage. When the excavation footage is 0.5 m, the plastic zone is distributed at the top and bottom of the pilot tunnel, and the distribution range is small, and the degree of rock yield failure is low. When the excavation footage is 1.0 m, the plastic zone increases in a small range at the arch bottom and arch foot. When the excavation footage is 1.5 m, the plastic zone is similar to that of 1.0 m, and the increase of excavation footage does not cause great disturbance to the rock mass. When the excavation is 2.0 m, the plastic zone of the vault extends to both sides, and the plastic zone of the bottom extends to the lower part. Almost all the rock mass around the pilot tunnel has plastic deformation.

(a)

(b)

(c)

(d)
In the process of construction, the width of each steel grill is set as 0.5 m. Therefore, the excavation footage of pilot tunnel ④ can be set as 1.5 m.
4. Discussion
4.1. The Accuracy of Algorithm
Figure 9 shows the influence of RF model on prediction accuracy under different parameters. The change of Mtry and Ntree will affect the prediction error of RF. The accuracy of prediction results is closely related to the selection of parameters. To ensure the accuracy of the prediction results, the optimal parameters of RF need to be selected.

Table 6 shows the comparison of CS-RF model with general RF model and ANN model. The previously mentioned algorithm has a good prediction effect on the displacement and stress of the station during construction. The prediction results of CS-RF, RF, and ANN are basically consistent with the field monitoring results.
The R2 values of displacement and stress predicted by CS-RF, RF, and ANN and monitoring results are shown in Figure 10. R2 value shows that the three techniques perform well in displacement and stress prediction. Comparing the R2 values of these methods, we can conclude that CS-RF has the best performance, and RF performs is better than ANN.

Table 7 shows the comparison of RF-CS model with particle swarm optimization (PSO) and genetic algorithm (GA). The previously mentioned algorithm has a good application effect in the process of back-analysis. Comparing the error results of these methods, we can conclude that RF-CS has the best performance. In the actual calculation, RF-CS algorithm is easy to implement and has fewer parameters, strong search ability, and fast convergence.
4.2. The Influence of CS Parameters on Optimization Results
In RF-CS algorithm, CS algorithm is more complex with varying influencing factors, playing a crucial role in RF-CS algorithm. Taking back-analysis of joint parameters as an example, the influence of parameter selection of CS algorithm on optimization results is analyzed in construction site.
In the CS, the probability of being discovered by the host (Pa) and the number of nests (N) have an effect on the convergence rate. Pa and N were settled, respectively. As shown in Figure 11, the iterative convergence curves under different Pa and N are presented, where the influences of Pa and N values on the search process were analyzed by comparing the convergence curves of different Pa and N values in the iteration process.

First, the number of nests is selected as N = 30. As can be seen, Pa is ranged from 0.15 to 0.45, with a stable iterative process. Notably, the convergence effect is promising but at varying convergence speeds. The number of iteration steps is small to reach convergence when Pa = 0.25. Then, Pa is fixed at 0.25 and N is set to 30∼100. When N = 50, the convergence rate is relatively good. The results show that the reasonable selection of CS parameters can improve the convergence speed. Pa = 0.25 and N = 50 are recommended in this study.
5. Conclusions
Aiming at the problem that the parameters of jointed of rock mass are difficult to obtain, a hybrid algorithm of RF and CS is established, and its calculation principle, flowchart, and application are studied. The main conclusions are as follows:(1)In this paper, U-J model, RF, and CS are integrated to establish the back-analysis method of joint parameters. The advantages of both RF and CS are fully utilized via this method, improving the calculation speed as well as avoiding the problem of local optimal solution.(2)This method has been successfully applied to the construction of metro station with arch cover method, and the back-analysis of joint parameters of surrounding rock has been realized. The back-analyzed results are highly consistent with the field monitoring data, with a relative error of 5.78%. Importantly, this method fulfills the engineering guidance since its calculation speed is fast, and its convergence is good with high accuracy.(3)The construction scheme of station pilot tunnel is optimized following the back-analysis results of joint parameters. The excavation footage of the middle tunnel is determined to be 1.5 m, which can ensure the safety of construction and effectively improve the construction efficiency. This conclusion has guiding significance for the study of arch cover method construction of metro station. It offers a feasible and effective approach for informational construction and feedback allocation of jointed surrounding rock tunnel.(4)The RF algorithm is optimized by CS algorithm. Compared with the general machine learning algorithms RF and ANN, the prediction error of CS-RF is relatively small. The parameter of CS plays an important role on the convergence rate. Lastly, Pa = 0.25 and N = 50 are recommended through calculation and analysis.
Data Availability
The data used in this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this article.
.
Acknowledgments
This work was supported by the Liaoning Revitalization Talents Program (no. XLYC1905015) and National Natural Science Foundation of China (nos. 52078093and 52108363).