# Creep Parameter Inversion and Long-Term Stability Analysis of Tunnel Based on GP-DE Intelligent Algorithm

Authors: Fengrui Zhang; Annan Jiang; Xinping Guo; Xiurong Yang
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

Aiming at the creep problem of Banshi Tunnel in Jilin province, the creep laws of rock are analyzed by the creep test, and the Cvsic model describing the creep characteristics of the tunnel is established. To obtain the creep parameters accurately, considering the advantages of Gauss process and differential evolution algorithm, coupling the two methods, a Gauss process-differential evolution intelligent inversion method is proposed. According to on-site monitoring data, the the creep parameters of the tunnel are accurately inverted. On this basis, the stability analysis of the tunnel and the selection of a reasonable construction plan are carried out. The research results show that to ensure the stability of the tunnel, the construction scheme of initial lining + pipe shed + advanced grouting anchor rod should be adopted. The research results have guiding significance for the long-term stability evaluation of the tunnel.


---

## Body

## 1. Introduction

The aging stability of the tunnel is influenced by the creep characteristics of surrounding rock, which attracts more and more attention. The practical application shows that the deformation of surrounding rock increases with time, which may lead to large deformation and even damage of tunnel support structure [1]. If the influence of creep effect is ignored, it will cause great difficulty in engineering construction. Therefore, the study of rock creep characteristics has guiding significance for the stability evaluation of the tunnel.The primary work of studying creep characteristics of the tunnel is how to determine the reasonable creep constitutive model and creep parameters. Many scholars have studied the rock creep characteristics through creep tests [2–5]. Zhang et al. [6] studied the creep damage characteristics of rocks under multilevel loads and obtained the creep damage evolution equation. Rassouli and Zoback [7] conducted creep tests on horizontal and vertical shale samples. Hu et al. [8] carried out the creep test on artificial layered specimens and obtained the creep deformation law of the test.According to the laboratory test results, a variety of rock creep models have been proposed [9–13]. However, for tunnel engineering, it is impossible to accurately obtain the creep parameters restricted by sampling disturbance, size effect, and test technology. Lin et al. [14] established a nonlinear viscoplastic element, connecting the plastic element and the viscous element in parallel. Zhang et al. [15] used discrete element software to establish serrated specimens and analyzed the creep characteristics of serrated specimens.The back analysis of tunnel displacement by using field monitoring data provides a new idea for the study of creep parameters. The numerical back analysis method and the analytical back analysis method are the traditional methods to study creep parameters [16, 17]. However, these methods have many problems, such as large computational workload and poor stability of solutions. In recent years, intelligent methods including Gaussian process (GP) and differential evolution (DE) algorithm have been widely used in geotechnical parameter inversion. Guan et al. [18] proposed a creep parameter inversion algorithm based on BN and GA. The author in [19] used the face mapping data to predict ground properties in a tunnel back analysis by an artificial neural network. Li et al. [20] used the Gaussian process model to predict tunnel water inrush. Liu and Liu [21] introduced a GA-CCGPR intelligent algorithm for tunnel construction.The GP algorithm can adaptively obtain parameters in the process of model construction and can well establish nonlinear mapping relationship, which is suitable for solving complex problems [22, 23]. However, using the conjugate gradient method to obtain the optimal hyperparameters of GP, the conjugate gradient method is dependent on the initial value, easily falling into local minimum and iterations may not converge. To overcome these shortcomings, DE is used to find the optimal hyperparameters of the GP. DE is a multiobjective optimization algorithm to solve the overall optimal solution in multidimensional space [24, 25].In this paper, considering the advantages of the two algorithms, a Gaussian process-differential evolution (GP-DE) intelligent inversion method for creep parameters is established. Based on the creep problem during the construction of Banshi Tunnel in Jilin Province, the creep constitutive model and parameter range are determined by the triaxial creep test. According to the field monitoring data, the creep parameters are inversed intelligently, and the reasonable creep parameters are obtained. On this basis, the stability of the tunnel is studied, and a reasonable construction scheme is selected.


## 2. The Creep Parameter Inversion Method Based on GP-DE

### 2.1. The Problem of Creep Parameter Inversion

The inversion of creep parameters of tunnel surrounding rock is essential to optimize the parameters and find the optimal solution. The upper and lower limits of model parameters are determined by the actual physical meaning. Assuming that there arem observations, the inversion optimization equation is as follows:(1)minEx1,x2,…,xn=min1m∑k=1mYk0−Yk,xka≤xk≤xkb,k=1,2,…,N,where m is the number of observations, Yk0 is the measured displacement, Yk is the calculated displacement, xk is the creep parameter, and N represents the number of parameters.The creep parameters of the tunnel are selected, and FLAC3D is used to numerically simulate the tunnel engineering and compare the calculated value with the field monitoring displacement value. If the difference is large, the creep parameters are reselected until the difference is small enough.When FLAC3D is used for numerical simulation, the selection of creep parameters requires a lot of forward calculation. Therefore, this paper uses the GP model to construct the response surface of creep parameters and tunnel displacement.


### 2.2. Gaussian Process

In Gaussian process, the joint probability distribution of random variableX and its corresponding process state f(X) obey n-dimensional Gaussian distribution. Statistical characteristics are determined by means of mean and covariance, and the expression is as follows:(2)ft∼GPμt,Ct,t′.D = (X, y) is taken as the training set of the Gaussian model, and learning the training set, the nonlinear mapping relationship between input variables and output vectors is established:(3)y=fx+ε,where ε is an independent random variable, ε∼N0,σn2, and σn2 is the variance.The prior distribution of the output value of the training sample is(4)y∼N0,K+σn2I,where K is the n-order symmetric positive definite covariance matrix; Kij of K = K (X, X) measures the correlation between xi and xj.Gaussian process (GP) is used to compute the predictive distribution of the function valuesy∗ at test points x∗. The joint distribution of target value and function value can be written as follows:(5)yy∗∼N0,KX,X+σn2I,KX,x∗KX,XT,KX∗,X∗,where K (X∗, X∗) is the covariance matrix of X∗ and K (X, x∗) is the order n × 1 covariance matrix of X and x∗.When the input valuex∗ of the test sample and the training set D are known, the output value y∗ of the test sample can be calculated through the Gaussian process:(6)y∗|x∗,D∼Nμy∗,σy∗2,μy∗=Kx∗,XKX,X+σn2I−1y=∑inαiKxi,x∗,where μy∗andσy∗ are the expectation and variance of y∗; α=K+σn2I−1y; and I is the identity matrix.The covariance function can represent the kernel function in the GP model. This article uses the squared exponential covariance function (SE):(7)ksexp,xq=σf2exp−12l2xp−xq2+σn2δpq,where xp and xq are the learning samples or prediction samples, l is the distance between xp and xq, σn is the standard deviation of the noise, σf is the local correlation, and δpq is a sign function.Hyperparametersσf and σn affect the GP training effect, which is expressed as follows:(8)minEθ=min∑h=1KGPhθ−YhYh,h=1,2,…,K,where GPhθ is the estimated output data, Yh is the actual output of the test sample, and θ=σf,σn represents the hyperparametric vector.


### 2.3. Differential Evolutionary Algorithm

The traditional method uses the conjugate gradient method to solve the optimal hyperparameters of GP, and then the conjugate gradient method depends on the initial value, which is easy to fall into the local minimum. To overcome these shortcomings, DE is used to search the optimal hyperparameters of GP.The DE algorithm was proposed by Storm and Price in 1995. The algorithm has excellent performance in dealing with nonlinear and complex problems. The specific optimization steps of DE are as follows.


#### 2.3.1. Generating the Initial Population

NP individuals are randomly generated in D-dimensional space, and the formula is as follows:(9)xi,j,0=xi,jL+randxi,jU−xi,jL,i=1,2,…,NP;j=1,2,…,D,where rand∈0,1 and xijU and xijL are the upper and lower limits of a single variable.


#### 2.3.2. Mutation Operation

Three individualsxr1,G, xr2,G, and xr3,G are randomly selected from the population for mutation operation, as shown in Figure 1. The variation vector is as follows:(10)Vi,G+1=Xr1,G+FXr2,G−Xr3,G,where r1, r2, and r3 are unequal integers between [1, NP] and F∈[0,1] is the variation factor.Figure 1
Schematic diagram of variation operation.


#### 2.3.3. Crossover Operation

A new vectorui,G+1 is generated by crossing the target vector xi,G with the variation vector vi,G+1 as follows:(11)ui,G+1=Vi,G+1,rj≤CRj=ni,Xi,G,rj>CR&j≠ni,where rj∈[0,1], CR∈[0,1], and ni is a random integer.


#### 2.3.4. Select Operation

The new generation of the population obtained by the selection operation is expressed as follows:(12)xi,G+1=ui,G+1,iffui,G+1≤fxi,G,xi,G,ffui,G+1>fxi,G,where f (·) represents the fitness function that corresponds to equation (8).


### 2.4. Creep Parameters Inversion Process

According to equation (1), Yk can be calculated using the GP model, which can be expressed by the following equation:(13)minEx1,x2,…,xN=min1m∑j=1mGPx1,x2,…,xNj0−Yj,xka≤xk≤xkb,k=1,2,...,N,where m is the number of measuring points, YJ is the monitoring value of the j-th measuring point, and Yj is the monitoring value of the j-th measuring point.Firstly, DE is used to search the optimal hyperparameters ofGP to improve the nonlinear regression ability of Gaussian process. Then, the mapping relationship between creep parameters and displacement is established through Gaussian training process. Finally, the creep parameters are obtained by searching the global space of the solution by the differential evolution algorithm. The creep parameter inversion process is shown in Figure 2, and the inversion steps are as follows:(1)
Based on the field data and laboratory tests, the creep characteristics of rock are analyzed, and the reasonable creep constitutive model and the range of creep parameters to be retrieved are determined.(2)
The inverse creep parameters are designed according to the orthogonal design principle.(3)
The numerical calculation model is established. The combination of creep parameters is calculated by FLAC3D software, obtaining the displacement and establishing the learning sample.(4)
Firstly, the kernel function of GP is determined, then the parameters of DE algorithm are set, and finally the super parameters of the kernel function are optimized.(5)
Mutation, crossover, and selection are carried out byDE to find the optimal super parameters of GP, improving the prediction ability.(6)
The mapping relationship between creep parameters and displacement is established through Gaussian training process.(7)
Then, theDE algorithm is used to inverse the creep parameters. The creep parameter to be randomly generated is the initial population, taking fitness function as evaluation index. When the fitness meets the requirements, creep parameters are considered to be found. If the fitness does not meet the requirements, mutation, crossover, selection, and other operations will be carried out until the maximum population iteration or fitness reaches the preset value.Figure 2
Flow chart of creep parameter inversion.


## 2.1. The Problem of Creep Parameter Inversion

The inversion of creep parameters of tunnel surrounding rock is essential to optimize the parameters and find the optimal solution. The upper and lower limits of model parameters are determined by the actual physical meaning. Assuming that there arem observations, the inversion optimization equation is as follows:(1)minEx1,x2,…,xn=min1m∑k=1mYk0−Yk,xka≤xk≤xkb,k=1,2,…,N,where m is the number of observations, Yk0 is the measured displacement, Yk is the calculated displacement, xk is the creep parameter, and N represents the number of parameters.The creep parameters of the tunnel are selected, and FLAC3D is used to numerically simulate the tunnel engineering and compare the calculated value with the field monitoring displacement value. If the difference is large, the creep parameters are reselected until the difference is small enough.When FLAC3D is used for numerical simulation, the selection of creep parameters requires a lot of forward calculation. Therefore, this paper uses the GP model to construct the response surface of creep parameters and tunnel displacement.


## 2.2. Gaussian Process

In Gaussian process, the joint probability distribution of random variableX and its corresponding process state f(X) obey n-dimensional Gaussian distribution. Statistical characteristics are determined by means of mean and covariance, and the expression is as follows:(2)ft∼GPμt,Ct,t′.D = (X, y) is taken as the training set of the Gaussian model, and learning the training set, the nonlinear mapping relationship between input variables and output vectors is established:(3)y=fx+ε,where ε is an independent random variable, ε∼N0,σn2, and σn2 is the variance.The prior distribution of the output value of the training sample is(4)y∼N0,K+σn2I,where K is the n-order symmetric positive definite covariance matrix; Kij of K = K (X, X) measures the correlation between xi and xj.Gaussian process (GP) is used to compute the predictive distribution of the function valuesy∗ at test points x∗. The joint distribution of target value and function value can be written as follows:(5)yy∗∼N0,KX,X+σn2I,KX,x∗KX,XT,KX∗,X∗,where K (X∗, X∗) is the covariance matrix of X∗ and K (X, x∗) is the order n × 1 covariance matrix of X and x∗.When the input valuex∗ of the test sample and the training set D are known, the output value y∗ of the test sample can be calculated through the Gaussian process:(6)y∗|x∗,D∼Nμy∗,σy∗2,μy∗=Kx∗,XKX,X+σn2I−1y=∑inαiKxi,x∗,where μy∗andσy∗ are the expectation and variance of y∗; α=K+σn2I−1y; and I is the identity matrix.The covariance function can represent the kernel function in the GP model. This article uses the squared exponential covariance function (SE):(7)ksexp,xq=σf2exp−12l2xp−xq2+σn2δpq,where xp and xq are the learning samples or prediction samples, l is the distance between xp and xq, σn is the standard deviation of the noise, σf is the local correlation, and δpq is a sign function.Hyperparametersσf and σn affect the GP training effect, which is expressed as follows:(8)minEθ=min∑h=1KGPhθ−YhYh,h=1,2,…,K,where GPhθ is the estimated output data, Yh is the actual output of the test sample, and θ=σf,σn represents the hyperparametric vector.


## 2.3. Differential Evolutionary Algorithm

The traditional method uses the conjugate gradient method to solve the optimal hyperparameters of GP, and then the conjugate gradient method depends on the initial value, which is easy to fall into the local minimum. To overcome these shortcomings, DE is used to search the optimal hyperparameters of GP.The DE algorithm was proposed by Storm and Price in 1995. The algorithm has excellent performance in dealing with nonlinear and complex problems. The specific optimization steps of DE are as follows.


### 2.3.1. Generating the Initial Population

NP individuals are randomly generated in D-dimensional space, and the formula is as follows:(9)xi,j,0=xi,jL+randxi,jU−xi,jL,i=1,2,…,NP;j=1,2,…,D,where rand∈0,1 and xijU and xijL are the upper and lower limits of a single variable.


### 2.3.2. Mutation Operation

Three individualsxr1,G, xr2,G, and xr3,G are randomly selected from the population for mutation operation, as shown in Figure 1. The variation vector is as follows:(10)Vi,G+1=Xr1,G+FXr2,G−Xr3,G,where r1, r2, and r3 are unequal integers between [1, NP] and F∈[0,1] is the variation factor.Figure 1
Schematic diagram of variation operation.


### 2.3.3. Crossover Operation

A new vectorui,G+1 is generated by crossing the target vector xi,G with the variation vector vi,G+1 as follows:(11)ui,G+1=Vi,G+1,rj≤CRj=ni,Xi,G,rj>CR&j≠ni,where rj∈[0,1], CR∈[0,1], and ni is a random integer.


### 2.3.4. Select Operation

The new generation of the population obtained by the selection operation is expressed as follows:(12)xi,G+1=ui,G+1,iffui,G+1≤fxi,G,xi,G,ffui,G+1>fxi,G,where f (·) represents the fitness function that corresponds to equation (8).


## 2.3.1. Generating the Initial Population

NP individuals are randomly generated in D-dimensional space, and the formula is as follows:(9)xi,j,0=xi,jL+randxi,jU−xi,jL,i=1,2,…,NP;j=1,2,…,D,where rand∈0,1 and xijU and xijL are the upper and lower limits of a single variable.


## 2.3.2. Mutation Operation

Three individualsxr1,G, xr2,G, and xr3,G are randomly selected from the population for mutation operation, as shown in Figure 1. The variation vector is as follows:(10)Vi,G+1=Xr1,G+FXr2,G−Xr3,G,where r1, r2, and r3 are unequal integers between [1, NP] and F∈[0,1] is the variation factor.Figure 1
Schematic diagram of variation operation.


## 2.3.3. Crossover Operation

A new vectorui,G+1 is generated by crossing the target vector xi,G with the variation vector vi,G+1 as follows:(11)ui,G+1=Vi,G+1,rj≤CRj=ni,Xi,G,rj>CR&j≠ni,where rj∈[0,1], CR∈[0,1], and ni is a random integer.


## 2.3.4. Select Operation

The new generation of the population obtained by the selection operation is expressed as follows:(12)xi,G+1=ui,G+1,iffui,G+1≤fxi,G,xi,G,ffui,G+1>fxi,G,where f (·) represents the fitness function that corresponds to equation (8).


## 2.4. Creep Parameters Inversion Process

According to equation (1), Yk can be calculated using the GP model, which can be expressed by the following equation:(13)minEx1,x2,…,xN=min1m∑j=1mGPx1,x2,…,xNj0−Yj,xka≤xk≤xkb,k=1,2,...,N,where m is the number of measuring points, YJ is the monitoring value of the j-th measuring point, and Yj is the monitoring value of the j-th measuring point.Firstly, DE is used to search the optimal hyperparameters ofGP to improve the nonlinear regression ability of Gaussian process. Then, the mapping relationship between creep parameters and displacement is established through Gaussian training process. Finally, the creep parameters are obtained by searching the global space of the solution by the differential evolution algorithm. The creep parameter inversion process is shown in Figure 2, and the inversion steps are as follows:(1)
Based on the field data and laboratory tests, the creep characteristics of rock are analyzed, and the reasonable creep constitutive model and the range of creep parameters to be retrieved are determined.(2)
The inverse creep parameters are designed according to the orthogonal design principle.(3)
The numerical calculation model is established. The combination of creep parameters is calculated by FLAC3D software, obtaining the displacement and establishing the learning sample.(4)
Firstly, the kernel function of GP is determined, then the parameters of DE algorithm are set, and finally the super parameters of the kernel function are optimized.(5)
Mutation, crossover, and selection are carried out byDE to find the optimal super parameters of GP, improving the prediction ability.(6)
The mapping relationship between creep parameters and displacement is established through Gaussian training process.(7)
Then, theDE algorithm is used to inverse the creep parameters. The creep parameter to be randomly generated is the initial population, taking fitness function as evaluation index. When the fitness meets the requirements, creep parameters are considered to be found. If the fitness does not meet the requirements, mutation, crossover, selection, and other operations will be carried out until the maximum population iteration or fitness reaches the preset value.Figure 2
Flow chart of creep parameter inversion.


## 3. Study Area

Banshi Tunnel, located in Jilin Province, is divided into left and right sides, with a left side of 1711 m and a right side of 1717 m. The terrain inclines from southeast to northwest with an elevation of 700–1000 m. It is a low mountain landform. The slope angle of the mountain body is 10–30°, the entrance and exit slope of the south side is 12–16°, and the entrance and exit slope of the north side is 20–25°. The surrounding rocks of the tunnel are mainly mixed gneiss with locally developed faults. Rock mass is broken, and its overall stability is poor, as shown in Figure3.Figure 3
Tunnel and surrounding rock.The LK49 + 885-LK49 + 915 area at the left entrance of the tunnel is characterized by weak rock belts, developed joints and fissures, poor joint planes, and loose and fractured rock mass. The tunnel is excavated by the CRD method, and the supporting method is advanced bolt + steel arch truss + shotcrete. After excavation and support, the creep characteristics of surrounding rocks are obvious, resulting in local collapse, support failure, lining cracking, and other phenomena, in Figure4.Figure 4
Tunnel accident. (a) Collapse. (b) Varch failure. (c) Support damage. (d) Lining rupture.
(a)(b)(c)(d)


## 4. Creep Model of Tunnel

The sample is taken from the left portal area of Banshi Tunnel. It is mixed gneiss, mainly composed of feldspar, quartz, and various dark minerals. Intact rock from the same tunnel face in the tunnel is selected and transported back to the laboratory and then processed into the cylindrical standard specimen with a diameter of 50 mm and a height of 100 mm, as shown in Figure5.Figure 5
Gneiss specimen.According to the in situ stress measurement of the tunnel, the confining pressure of the creep test is determined to be 1 MPa. Firstly, the confining pressure is loaded to 1 MPa; then, the axial pressure is loaded in stages, the loading rate is 0.5 MPa/min, and the initial value is 10 MPa. When the creep deformation rate is less than 0.001 mm/d, the next loading is carried out,Δσ=20MPa, until the specimen is damaged. Figure 6 shows the results of creep tests.Figure 6
Creep test curve of gneiss.The specimen produces instantaneous elastic deformation under stress loading, and then the deformation slowly increases with time under the constant stress, but the creep rate gradually decreases and tends to be stable, showing the characteristics of decay creep and stable creep. With the increase in axial stress, the deformation of sample increases gradually. Loading to the last load, the deformation of the specimen increases rapidly, the rock will be suddenly destroyed, and the failure process is extremely short, and there is no obvious accelerated creep stage.The whole creep process of gneiss presents four stages: instantaneous elastic deformation, attenuation creep, decay creep, and sudden failure. Therefore, the Cvsic model is selected to describe the creep characteristics of gneiss, as shown in Figure7. In the Cvsic model, viscoelastic constitutive relation corresponds to the Burger model and plastic constitutive relation corresponds to the Mohr–Coulomb model.Figure 7
The creep model of Cvsic.For the Cvsic model, the strain tensor can be written as(14)e˙ij=e˙ijK+e˙ijM+e˙ijP,where eijK is the Kelvin body strain, eijM is the Maxwell body strain, and eijP is the Mohr–Coulomb body strain.For the viscoelastic unit, the constitutive laws for the Kelvin and Maxwell units are formulated by(15)Sij=2ηKe˙ijK+2GKe˙ijK,(16)e˙ijM=S˙ij2GM+Sij2ηM,where Sij is the stress tensors, GK and ηK are the shear modulus and the viscosity of Kelvin, and GM and ηM are the shear modulus and the viscosity of Maxwell.According to the Boltzmann superposition principle, the creep curves of gneiss under step load are transformed and the creep curves under different loading conditions are obtained. Table1 shows the creep parameters of the Cvsic model that are identified by Istopt software.Table 1
Creep parameters of the Cvsic model.
SampleAxial stress (MPa)Gk (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)Gneiss50.743.900.9266.5100.533.562.71192.61150.398.584.33284.15200.289.727.06442.98250.267.478.80232.15To verify the accuracy of the model and its parameters, the test curve is compared with the model curve, as shown in Figure8. The model curve of the Cvsic model is in good agreement with the test curve, which shows that the Cvsic model can accurately describe the creep characteristics of gneiss. In addition, on the basis of the creep parameters identified by laboratory tests, combined with the field investigation and the reference of similar engineering experience, the range of creep parameters of the tunnel can be determined.Figure 8
Comparison of test curve and fitting curve of the Cvsic model.


## 5. Back Analysis Using GP-DE Algorithm

### 5.1. Numerical Simulation Model

According to the geological characteristics of tunnel LK49 + 885∼LK49 + 915, a numerical model is established. The tunnel section is horseshoe shaped with a width of 10.5 m and a height of 7.8 m.X, Y, and Z direction constraints are applied to the bottom of the model, and normal displacement constraints are applied to the side. The model has a total of 40541 cells and 28230 nodes, as shown in Figure 9.Figure 9
(a) Tunnel numerical model. (b) Layout of monitoring section.To monitor the tunnel displacement, one monitoring section is arranged every 15 m, and three sections (LK49 + 885, LK49 + 900, and LK49 + 915) are installed in total. Three displacement monitoring points are arranged at arch foot, arch waist, and arch crown of each section, as shown in Figure9.According to the test analysis results, the Cvsic model can be used to study tunnel creep. The creep parametersGK, GM, ηK, and ηM of the Cvsic model are inversed intelligently. Table 2 shows the range of creep parameters to be inversed.Table 2
Range of parameters to be back analyzed.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)0.1∼2.51∼50.5∼4.510∼50The orthogonal parameter scheme and uniform design parameters are shown in Tables3 and 4. Table 3 is the training data, and Table 4 is the verification data. Using FLAC3D software to solve the data, the displacement of the tunnel monitoring section is calculated. Among them, the monitoring data of section LK49 + 885 are used for the inversion of creep parameters. The monitoring data of section LK49 + 900 are compared with the inversion results.Table 3
The orthogonal parameter sample.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d10.110.5104.269.6215.3621.636.5815.1419.4723.789.4117.4823.3129.5220.121.5203.979.2314.8420.846.1914.7718.7523.338.9517.0322.7829.2530.132.5303.798.9114.6120.525.7413.917.3621.998.4816.6722.2128.8240.143.5403.538.8414.2620.065.3713.4816.8421.648.0116.2821.8228.0250.154.5503.238.5113.7419.635.0813.0416.2121.037.4515.4721.0927.2660.711.5303.589.3114.4320.785.9414.2117.6321.778.4619.721.9228.4970.722.5403.378.6614.0320.095.3113.5416.4021.467.7615.7821.5527.5180.733.5503.118.3213.4819.374.7412.7215.9420.787.2115.0520.9826.7990.744.5103.178.4713.6119.494.9312.8516.1520.977.315.2321.3226.96100.750.5203.248.6213.8519.805.131316.3021.057.7215.4422.5327.52111.312.5502.767.9813.2319.314.6312.6415.7020.667.1314.9920.6326.45121.323.5102.828.2213.5319.474.8112.9515.3120.967.3715.5122.0526.69131.334.5202.677.6213.0318.504.2512.0314.9419.776.9314.920.6425.13141.340.5302.547.2012.7218.193.7411.2314.7118.126.6215.6119.6224.87151.351.5402.236.8512.4517.773.4010.8714.0317.726.2714.3923.3123.01161.913.5201.716.6411.4917.083.1210.2913.6517.305.8013.8618.6422.78171.924.5301.596.3610.9516.742.909.4812.8317.015.6213.4317.8021.76181.930.5401.536.1910.6416.372.849.3912.3915.745.5112.8217.2821.32191.941.5501.445.7310.0215.792.619.0411.7814.775.2811.7316.4320.20201.952.5101.515.8910.4016.152.719.2612.2915.255.3212.3116.8220.62212.514.5401.365.4310.1414.452.418.7512.7014.324.3811.5714.5819.17222.520.5501.255.249.8713.972.058.5312.1213.784.0910.7313.9418.65232.531.5101.195.169.6613.451.988.1811.9413.114.0210.1213.2918.13242.542.5201.105.049.2812.781.817.9311.3512.683.789.6112.3117.10252.553.5301.024.678.8112.121.937.0910.5312.013.259.1311.4816.43Table 4
The uniform design parameter sample.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d10.121.5303.859.0114.7320.615.9214.3518.0422.838.6516.9022.5429.0120.733.5203.248.5213.8819.875.1413.1216.2421.187.5215.4518.8822.1931.310.5502.306.2711.2517.393.5311.6215.5819.046.2414.0721.4127.2341.952.5101.515.8910.4016.152.719.2612.2915.255.3212.3116.8220.6252.544.5401.064.829.0312.451.877.6510.8012.323.489.3711.8816.80


### 5.2. Back Analysis Results

The creep parameters of tunnel surrounding rock are inversed by the GP-DE algorithm. Setting relevant parameters as follows: the population sizeNP = 100, the variation factor F = 0.6, the cross factor CR = 0.9, the maximum evolution algebra itermax = 100, and the kernel function is square exponential covariance function. Table 5 shows the engineering monitoring data. Based on the calculation results of the orthogonal parameters schemes, the nonlinear implicit relation between creep parameters and displacement is established. The results of the uniformly parameters schemes are used for verification. Table 6 shows the creep parameters that are obtained by inversion.Table 5
The tunnel monitoring data.
Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d1.516.0310.3516.172.829.4713.1817.735.4712.7017.0520.41Table 6
Inversion of creep parameters.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)2.172.104.5016.65According to the inversion results, FLAC3D software is used again for forwarding analysis, and the calculated displacement of section LK49 + 900 is compared with the measured displacement, in Figure 10. The calculated curves are in good agreement with the monitoring curves, and the values are basically the same, indicating that the creep constitutive model and the creep parameters obtained by inversion are basically reasonable.Figure 10
Comparison between the calculated value and the measured value.To verify the accuracy of the GP-DE model, the GP-DE model is compared with GP and ANN models, and the results are shown in Table7.Table 7
Comparison results of different algorithms.
Monitoring pointTime (d)Monitoring displacement (mm)GP-DE (mm)Relative error (%)GP (mm)Relative error (%)ANN (mm)Relative error (%)131.511.463.41.454.21.455.6126.035.862.85.813.75.966.42410.3510.151.99.993.510.645.84816.1715.772.515.772.915.756.1232.822.752.62.713.82.995.9129.479.183.18.995.110.025.82413.1812.961.712.366.214.026.44817.7317.342.217.073.718.326.7335.475.371.85.185.35.816.11212.7012.372.612.273.413.516.42417.0516.562.916.185.118.065.94820.4119.703.518.446.421.756.6As shown in Table7, the GP-DE model has the highest inversion accuracy of creep parameters, while the inversion ability of the GP algorithm and ANN algorithm is general, and the inversion accuracy of the GP-DE model is significantly improved.The fitness value of the GP-DE algorithm group under different iteration times is shown in Figure11. It can be seen from the figure that the fitness value tends to converge as the number of iterations increases, and the distribution in the space gradually decreases.Figure 11
Fitness values with different evolution. (a) 1stevolution. (b) 2th evolution. (c) 10th evolution. (d) 20th evolution.
(a)(b)(c)(d)


### 5.3. Reasonable Support Scheme of Tunnel

To ensure the stability of the tunnel, a reasonable support scheme is selected and the plastic zone of surrounding rock is compared under different support schemes, as shown in Figure12. It can be seen from the figure that the prestressed anchor rod plays a small role, while the grouting reinforcement plays an obvious role. Therefore, it is recommended to adopt the scheme of initial lining + pipe shed + advanced grouting anchor rod.Figure 12
Plastic zone of the tunnel with different support modes. (a) Initial lining. (b) Initial lining + pipe shed. (c) Initial lining + pipe shed + prestressed anchor rod. (d) Initial lining + pipe shed + advanced grouting anchor rod.
(a)(b)(c)(d)Figure13 shows the creep deformation curve at the tunnel crown. The numerical calculation value is close to the measured monitoring value, and the deformation trend is basically the same, which can reflect the creep characteristics of the surrounding rock of the tunnel. In the early stage of creep, the creep deformation of the tunnel vault increases rapidly. With the increase in time, the deformation gradually decreases and tends to be stable, which is basically consistent with the conclusion of the indoor test.Figure 13
Creep deformation curve of the tunnel crown.According to JTG F60-2009, the reasonable supporting time of the tunnel secondary lining is taken as the time when the deformation value of the tunnel reaches 80% of the final deformation. Whent⟶∞, the vault settlement is 24.15 mm, so the time 24 d corresponding to 80% of the final deformation is the reasonable time of secondary lining support. In the actual construction of the tunnel, the time of the second lining support is longer than this value, which shows that the support of the second lining is relatively lagging behind and the design of the primary lining of the tunnel is conservative. Therefore, it is necessary to speed up the construction of secondary lining or reduce the thickness of primary lining.


## 5.1. Numerical Simulation Model

According to the geological characteristics of tunnel LK49 + 885∼LK49 + 915, a numerical model is established. The tunnel section is horseshoe shaped with a width of 10.5 m and a height of 7.8 m.X, Y, and Z direction constraints are applied to the bottom of the model, and normal displacement constraints are applied to the side. The model has a total of 40541 cells and 28230 nodes, as shown in Figure 9.Figure 9
(a) Tunnel numerical model. (b) Layout of monitoring section.To monitor the tunnel displacement, one monitoring section is arranged every 15 m, and three sections (LK49 + 885, LK49 + 900, and LK49 + 915) are installed in total. Three displacement monitoring points are arranged at arch foot, arch waist, and arch crown of each section, as shown in Figure9.According to the test analysis results, the Cvsic model can be used to study tunnel creep. The creep parametersGK, GM, ηK, and ηM of the Cvsic model are inversed intelligently. Table 2 shows the range of creep parameters to be inversed.Table 2
Range of parameters to be back analyzed.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)0.1∼2.51∼50.5∼4.510∼50The orthogonal parameter scheme and uniform design parameters are shown in Tables3 and 4. Table 3 is the training data, and Table 4 is the verification data. Using FLAC3D software to solve the data, the displacement of the tunnel monitoring section is calculated. Among them, the monitoring data of section LK49 + 885 are used for the inversion of creep parameters. The monitoring data of section LK49 + 900 are compared with the inversion results.Table 3
The orthogonal parameter sample.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d10.110.5104.269.6215.3621.636.5815.1419.4723.789.4117.4823.3129.5220.121.5203.979.2314.8420.846.1914.7718.7523.338.9517.0322.7829.2530.132.5303.798.9114.6120.525.7413.917.3621.998.4816.6722.2128.8240.143.5403.538.8414.2620.065.3713.4816.8421.648.0116.2821.8228.0250.154.5503.238.5113.7419.635.0813.0416.2121.037.4515.4721.0927.2660.711.5303.589.3114.4320.785.9414.2117.6321.778.4619.721.9228.4970.722.5403.378.6614.0320.095.3113.5416.4021.467.7615.7821.5527.5180.733.5503.118.3213.4819.374.7412.7215.9420.787.2115.0520.9826.7990.744.5103.178.4713.6119.494.9312.8516.1520.977.315.2321.3226.96100.750.5203.248.6213.8519.805.131316.3021.057.7215.4422.5327.52111.312.5502.767.9813.2319.314.6312.6415.7020.667.1314.9920.6326.45121.323.5102.828.2213.5319.474.8112.9515.3120.967.3715.5122.0526.69131.334.5202.677.6213.0318.504.2512.0314.9419.776.9314.920.6425.13141.340.5302.547.2012.7218.193.7411.2314.7118.126.6215.6119.6224.87151.351.5402.236.8512.4517.773.4010.8714.0317.726.2714.3923.3123.01161.913.5201.716.6411.4917.083.1210.2913.6517.305.8013.8618.6422.78171.924.5301.596.3610.9516.742.909.4812.8317.015.6213.4317.8021.76181.930.5401.536.1910.6416.372.849.3912.3915.745.5112.8217.2821.32191.941.5501.445.7310.0215.792.619.0411.7814.775.2811.7316.4320.20201.952.5101.515.8910.4016.152.719.2612.2915.255.3212.3116.8220.62212.514.5401.365.4310.1414.452.418.7512.7014.324.3811.5714.5819.17222.520.5501.255.249.8713.972.058.5312.1213.784.0910.7313.9418.65232.531.5101.195.169.6613.451.988.1811.9413.114.0210.1213.2918.13242.542.5201.105.049.2812.781.817.9311.3512.683.789.6112.3117.10252.553.5301.024.678.8112.121.937.0910.5312.013.259.1311.4816.43Table 4
The uniform design parameter sample.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d10.121.5303.859.0114.7320.615.9214.3518.0422.838.6516.9022.5429.0120.733.5203.248.5213.8819.875.1413.1216.2421.187.5215.4518.8822.1931.310.5502.306.2711.2517.393.5311.6215.5819.046.2414.0721.4127.2341.952.5101.515.8910.4016.152.719.2612.2915.255.3212.3116.8220.6252.544.5401.064.829.0312.451.877.6510.8012.323.489.3711.8816.80


## 5.2. Back Analysis Results

The creep parameters of tunnel surrounding rock are inversed by the GP-DE algorithm. Setting relevant parameters as follows: the population sizeNP = 100, the variation factor F = 0.6, the cross factor CR = 0.9, the maximum evolution algebra itermax = 100, and the kernel function is square exponential covariance function. Table 5 shows the engineering monitoring data. Based on the calculation results of the orthogonal parameters schemes, the nonlinear implicit relation between creep parameters and displacement is established. The results of the uniformly parameters schemes are used for verification. Table 6 shows the creep parameters that are obtained by inversion.Table 5
The tunnel monitoring data.
Displacement value of 1 (mm)Displacement value of 2 (mm)Displacement value of 3 (mm)3 d12 d24 d48 d3 d12 d24 d48 d3 d12 d24 d48 d1.516.0310.3516.172.829.4713.1817.735.4712.7017.0520.41Table 6
Inversion of creep parameters.
GK (GPa)GM (GPa)ηK (GPa·d)ηM (GPa·d)2.172.104.5016.65According to the inversion results, FLAC3D software is used again for forwarding analysis, and the calculated displacement of section LK49 + 900 is compared with the measured displacement, in Figure 10. The calculated curves are in good agreement with the monitoring curves, and the values are basically the same, indicating that the creep constitutive model and the creep parameters obtained by inversion are basically reasonable.Figure 10
Comparison between the calculated value and the measured value.To verify the accuracy of the GP-DE model, the GP-DE model is compared with GP and ANN models, and the results are shown in Table7.Table 7
Comparison results of different algorithms.
Monitoring pointTime (d)Monitoring displacement (mm)GP-DE (mm)Relative error (%)GP (mm)Relative error (%)ANN (mm)Relative error (%)131.511.463.41.454.21.455.6126.035.862.85.813.75.966.42410.3510.151.99.993.510.645.84816.1715.772.515.772.915.756.1232.822.752.62.713.82.995.9129.479.183.18.995.110.025.82413.1812.961.712.366.214.026.44817.7317.342.217.073.718.326.7335.475.371.85.185.35.816.11212.7012.372.612.273.413.516.42417.0516.562.916.185.118.065.94820.4119.703.518.446.421.756.6As shown in Table7, the GP-DE model has the highest inversion accuracy of creep parameters, while the inversion ability of the GP algorithm and ANN algorithm is general, and the inversion accuracy of the GP-DE model is significantly improved.The fitness value of the GP-DE algorithm group under different iteration times is shown in Figure11. It can be seen from the figure that the fitness value tends to converge as the number of iterations increases, and the distribution in the space gradually decreases.Figure 11
Fitness values with different evolution. (a) 1stevolution. (b) 2th evolution. (c) 10th evolution. (d) 20th evolution.
(a)(b)(c)(d)


## 5.3. Reasonable Support Scheme of Tunnel

To ensure the stability of the tunnel, a reasonable support scheme is selected and the plastic zone of surrounding rock is compared under different support schemes, as shown in Figure12. It can be seen from the figure that the prestressed anchor rod plays a small role, while the grouting reinforcement plays an obvious role. Therefore, it is recommended to adopt the scheme of initial lining + pipe shed + advanced grouting anchor rod.Figure 12
Plastic zone of the tunnel with different support modes. (a) Initial lining. (b) Initial lining + pipe shed. (c) Initial lining + pipe shed + prestressed anchor rod. (d) Initial lining + pipe shed + advanced grouting anchor rod.
(a)(b)(c)(d)Figure13 shows the creep deformation curve at the tunnel crown. The numerical calculation value is close to the measured monitoring value, and the deformation trend is basically the same, which can reflect the creep characteristics of the surrounding rock of the tunnel. In the early stage of creep, the creep deformation of the tunnel vault increases rapidly. With the increase in time, the deformation gradually decreases and tends to be stable, which is basically consistent with the conclusion of the indoor test.Figure 13
Creep deformation curve of the tunnel crown.According to JTG F60-2009, the reasonable supporting time of the tunnel secondary lining is taken as the time when the deformation value of the tunnel reaches 80% of the final deformation. Whent⟶∞, the vault settlement is 24.15 mm, so the time 24 d corresponding to 80% of the final deformation is the reasonable time of secondary lining support. In the actual construction of the tunnel, the time of the second lining support is longer than this value, which shows that the support of the second lining is relatively lagging behind and the design of the primary lining of the tunnel is conservative. Therefore, it is necessary to speed up the construction of secondary lining or reduce the thickness of primary lining.


## 6. Discussion

### 6.1. The Prediction Accuracy of GP Model

Figure14 shows the influence of super parameters on the inversion accuracy of the GP model. It can be seen from the figure that σf and σn play a key role in the inversion accuracy. When lnσf=5⋅42 and lnσn=4⋅83, the inversion accuracy error is the lowest. This shows that it is very necessary to select reasonable super parameters to ensure the inversion accuracy of the GP model.Figure 14
The influence ofσf and σn on the accuracy of GP inversion.


### 6.2. Influence of Differential Evolution Parameters on Optimization Results

When the GP-DE algorithm is used for the inversion of tunnel creep parameters, the DE algorithm plays a key role. This paper analyzes the influence of mutation factorF, crossover factor CR, population size NP, and different difference strategies on the DE algorithm and selects reasonable parameters.


#### 6.2.1. Influence ofF and CR on Inversion Accuracy

Figure15 shows the effects of different F and CR on the iterative curve. It can be seen from the figure that the variation factor F and cross factor CR has an important impact on the convergence speed. When F = 0.6 and CR = 0.5∼0.9, the iterative process converges well and CR 0.9 converges fastest; when CR = 0.9 and F = 0.5∼0.9, F = 0.6 converges faster. The results show that the convergence performance of the DE algorithm is better when CR = 0.9 and F = 0.6 are selected.Figure 15
Iteration curve with change ofCR and F: (a) F = 0.6. (b) CR = 0.9.
(a)(b)


#### 6.2.2. Influence ofNP on Inversion Accuracy

Figure16 shows the influence of different population sizes on the iterative convergence curve. When NP = 1∼100, the optimization accuracy of the model increases gradually with the increase in population size, but the iteration rate decreases; when NP = 100∼200, the optimization accuracy of the model does not change significantly with the increase in population size, so the population size NP = 100 is comprehensively considered.Figure 16
Iteration curve of different population.


#### 6.2.3. Influence of Differentiation Strategy on Inversion Accuracy

In the DE algorithm, the different differential strategies are as follows:(17)DErand/1:vi,g=xr1,g+Fxr2,g−xr3,g,DEbest/1:vi,g=xbest,g+Fxr2,g−xr3,g,DErand/2:vi,g=xr1,g+Fxr2,g−xr3,g+xr4,g−xr5,g,DEbest/2:vi,g=xbest,g+Fxr1,g−xr2,g+xr3,g−xr4,g,DErand−to−best/2:vi,g=xr1,g+Fxbest,g−xr2,g+xr3,g−xr4,g,where xbest,g is the best individual and xr1,xr2,xr3,xr4,andxr5 are random individuals.The iterative curves of different difference strategies are shown in Figure17. It can be seen from the figure that the difference strategy has an important impact on the convergence speed and optimization accuracy. Compared with other difference strategies, DE/Best/1 has faster convergence speed and higher optimization accuracy.Figure 17
Different differential strategies.


## 6.1. The Prediction Accuracy of GP Model

Figure14 shows the influence of super parameters on the inversion accuracy of the GP model. It can be seen from the figure that σf and σn play a key role in the inversion accuracy. When lnσf=5⋅42 and lnσn=4⋅83, the inversion accuracy error is the lowest. This shows that it is very necessary to select reasonable super parameters to ensure the inversion accuracy of the GP model.Figure 14
The influence ofσf and σn on the accuracy of GP inversion.


## 6.2. Influence of Differential Evolution Parameters on Optimization Results

When the GP-DE algorithm is used for the inversion of tunnel creep parameters, the DE algorithm plays a key role. This paper analyzes the influence of mutation factorF, crossover factor CR, population size NP, and different difference strategies on the DE algorithm and selects reasonable parameters.


### 6.2.1. Influence ofF and CR on Inversion Accuracy

Figure15 shows the effects of different F and CR on the iterative curve. It can be seen from the figure that the variation factor F and cross factor CR has an important impact on the convergence speed. When F = 0.6 and CR = 0.5∼0.9, the iterative process converges well and CR 0.9 converges fastest; when CR = 0.9 and F = 0.5∼0.9, F = 0.6 converges faster. The results show that the convergence performance of the DE algorithm is better when CR = 0.9 and F = 0.6 are selected.Figure 15
Iteration curve with change ofCR and F: (a) F = 0.6. (b) CR = 0.9.
(a)(b)


### 6.2.2. Influence ofNP on Inversion Accuracy

Figure16 shows the influence of different population sizes on the iterative convergence curve. When NP = 1∼100, the optimization accuracy of the model increases gradually with the increase in population size, but the iteration rate decreases; when NP = 100∼200, the optimization accuracy of the model does not change significantly with the increase in population size, so the population size NP = 100 is comprehensively considered.Figure 16
Iteration curve of different population.


### 6.2.3. Influence of Differentiation Strategy on Inversion Accuracy

In the DE algorithm, the different differential strategies are as follows:(17)DErand/1:vi,g=xr1,g+Fxr2,g−xr3,g,DEbest/1:vi,g=xbest,g+Fxr2,g−xr3,g,DErand/2:vi,g=xr1,g+Fxr2,g−xr3,g+xr4,g−xr5,g,DEbest/2:vi,g=xbest,g+Fxr1,g−xr2,g+xr3,g−xr4,g,DErand−to−best/2:vi,g=xr1,g+Fxbest,g−xr2,g+xr3,g−xr4,g,where xbest,g is the best individual and xr1,xr2,xr3,xr4,andxr5 are random individuals.The iterative curves of different difference strategies are shown in Figure17. It can be seen from the figure that the difference strategy has an important impact on the convergence speed and optimization accuracy. Compared with other difference strategies, DE/Best/1 has faster convergence speed and higher optimization accuracy.Figure 17
Different differential strategies.


## 6.2.1. Influence ofF and CR on Inversion Accuracy

Figure15 shows the effects of different F and CR on the iterative curve. It can be seen from the figure that the variation factor F and cross factor CR has an important impact on the convergence speed. When F = 0.6 and CR = 0.5∼0.9, the iterative process converges well and CR 0.9 converges fastest; when CR = 0.9 and F = 0.5∼0.9, F = 0.6 converges faster. The results show that the convergence performance of the DE algorithm is better when CR = 0.9 and F = 0.6 are selected.Figure 15
Iteration curve with change ofCR and F: (a) F = 0.6. (b) CR = 0.9.
(a)(b)


## 6.2.2. Influence ofNP on Inversion Accuracy

Figure16 shows the influence of different population sizes on the iterative convergence curve. When NP = 1∼100, the optimization accuracy of the model increases gradually with the increase in population size, but the iteration rate decreases; when NP = 100∼200, the optimization accuracy of the model does not change significantly with the increase in population size, so the population size NP = 100 is comprehensively considered.Figure 16
Iteration curve of different population.


## 6.2.3. Influence of Differentiation Strategy on Inversion Accuracy

In the DE algorithm, the different differential strategies are as follows:(17)DErand/1:vi,g=xr1,g+Fxr2,g−xr3,g,DEbest/1:vi,g=xbest,g+Fxr2,g−xr3,g,DErand/2:vi,g=xr1,g+Fxr2,g−xr3,g+xr4,g−xr5,g,DEbest/2:vi,g=xbest,g+Fxr1,g−xr2,g+xr3,g−xr4,g,DErand−to−best/2:vi,g=xr1,g+Fxbest,g−xr2,g+xr3,g−xr4,g,where xbest,g is the best individual and xr1,xr2,xr3,xr4,andxr5 are random individuals.The iterative curves of different difference strategies are shown in Figure17. It can be seen from the figure that the difference strategy has an important impact on the convergence speed and optimization accuracy. Compared with other difference strategies, DE/Best/1 has faster convergence speed and higher optimization accuracy.Figure 17
Different differential strategies.


## 7. Conclusions

(1)
According to the creep test of gneiss, the creep model and parameter range of the tunnel are determined. The results show that the whole creep process of gneiss presents four stages: instantaneous elastic deformation, attenuation creep, stable creep, and sudden failure. The Cvsic model can reflect the creep characteristics of the tunnel.(2)
Taking full advantage of Gauss process and differential evolution algorithm and coupling the two methods, a Gauss process-differential evolution intelligent inversion method is proposed. The algorithm performance test and parameter analysis are carried out. SelectingCR = 0.9, F = 0.6, and NP = 100 can improve the convergence speed of the algorithm, save the calculation time, and avoid falling into the local optimal solution.(3)
The stability analysis of the tunnel and the selection of a reasonable construction plan are carried out. The prestressed anchor rod plays a small role, while the grouting reinforcement plays an obvious role. It is recommended to adopt the scheme of initial lining + pipe shed + advanced grouting anchor rod. In addition, the reasonable support time for the secondary lining of the tunnel is 24 days. In the actual construction of the tunnel, the secondary lining support is relatively lagging, which can speed up the construction of the secondary lining.



---

*Source: 3769474-2021-09-21.xml*