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

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.


Introduction
e aging stability of the tunnel is influenced by the creep characteristics of surrounding rock, which attracts more and more attention.
e 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. erefore, the study of rock creep characteristics has guiding significance for the stability evaluation of the tunnel. e 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][3][4][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][10][11][12][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. e back analysis of tunnel displacement by using field monitoring data provides a new idea for the study of creep parameters. e 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. e 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.
e 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.

The Creep Parameter Inversion Method
Based on GP-DE

e Problem of Creep Parameter
Inversion. e inversion of creep parameters of tunnel surrounding rock is essential to optimize the parameters and find the optimal solution. e upper and lower limits of model parameters are determined by the actual physical meaning. Assuming that there are m observations, the inversion optimization equation is as follows: where m is the number of observations, Y k 0 is the measured displacement, Y k is the calculated displacement, x k is the creep parameter, and N represents the number of parameters.
e creep parameters of the tunnel are selected, and FLAC 3D 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 FLAC 3D is used for numerical simulation, the selection of creep parameters requires a lot of forward calculation. erefore, this paper uses the GP model to construct the response surface of creep parameters and tunnel displacement.

Gaussian Process.
In Gaussian process, the joint probability distribution of random variable X 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: 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: where ε is an independent random variable, ε ∼ N(0, σ 2 n ), and σ 2 n is the variance. e prior distribution of the output value of the training sample is where K is the n-order symmetric positive definite covariance matrix; K ij of K � K (X, X) measures the correlation between x i and x j . Gaussian process (GP) is used to compute the predictive distribution of the function values y * at test points x * . e joint distribution of target value and function value can be written as follows: 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 value x * 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: where μ y * and σ y * are the expectation and variance of y * ; α � (K + σ n 2 I) − 1 y; and I is the identity matrix. e covariance function can represent the kernel function in the GP model. is article uses the squared exponential covariance function (SE): where x p and x q are the learning samples or prediction samples, l is the distance between x p and x q , σ 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: 2 Advances in Materials Science and Engineering where GP h (θ) is the estimated output data, Y h is the actual output of the test sample, and θ � (σ f , σ n ) represents the hyperparametric vector.

Differential Evolutionary
Algorithm. e 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.
e DE algorithm was proposed by Storm and Price in 1995. e algorithm has excellent performance in dealing with nonlinear and complex problems. e 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: where r and ∈ [0, 1] and x ij U and x ij L are the upper and lower limits of a single variable.

Mutation Operation.
ree individuals x r1,G , x r2,G , and x r3,G are randomly selected from the population for mutation operation, as shown in Figure 1. e variation vector is as follows: where r 1 , r 2 , and r 3 are unequal integers between [1,NP] and F∈[0,1] is the variation factor.

Crossover Operation.
A new vector u i,G+1 is generated by crossing the target vector x i,G with the variation vector v i,G+1 as follows: 1], and n i is a random integer.

Select Operation.
e new generation of the population obtained by the selection operation is expressed as follows: where f (·) represents the fitness function that corresponds to equation (8).

Creep Parameters Inversion
Process. According to equation (1), Y k can be calculated using the GP model, which can be expressed by the following equation: where m is the number of measuring points, YJ is the monitoring value of the j-th measuring point, and Y j is the monitoring value of the j-th measuring point. Firstly, DE is used to search the optimal hyperparameters of GP to improve the nonlinear regression ability of Gaussian process. en, 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. e creep parameter

Advances in Materials Science and Engineering
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. (5) Mutation, crossover, and selection are carried out by DE to find the optimal super parameters of GP, improving the prediction ability. (6) e mapping relationship between creep parameters and displacement is established through Gaussian training process. (7) en, the DE algorithm is used to inverse the creep parameters. e 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.

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. e terrain inclines from southeast to northwest with an elevation of 700-1000 m. It is a low mountain landform. e 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°. e 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 Figure 3. e 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. e 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

Creep Model of Tunnel
e 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 Figure 5.
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. e specimen produces instantaneous elastic deformation under stress loading, and then the deformation slowly     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. e whole creep process of gneiss presents four stages: instantaneous elastic deformation, attenuation creep, decay creep, and sudden failure.
erefore, the Cvsic model is selected to describe the creep characteristics of gneiss, as shown in Figure 7. In the Cvsic model, viscoelastic constitutive relation corresponds to the Burger model and plastic constitutive relation corresponds to the Mohr-Coulomb model. For the Cvsic model, the strain tensor can be written as where e K ij is the Kelvin body strain, e M ij is the Maxwell body strain, and e P ij is the Mohr-Coulomb body strain. For the viscoelastic unit, the constitutive laws for the Kelvin and Maxwell units are formulated by where S ij is the stress tensors, G K and η K are the shear modulus and the viscosity of Kelvin, and G M 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. Table 1 shows the creep parameters of the Cvsic model that are identified by Istopt software.
To verify the accuracy of the model and its parameters, the test curve is compared with the model curve, as shown in

Numerical Simulation Model.
According to the geological characteristics of tunnel LK49 + 885∼LK49 + 915, a numerical model is established. e 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. e model has a total of 40541 cells and 28230 nodes, as shown in Figure 9.
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. ree displacement monitoring points are arranged at arch foot, arch waist, and arch crown of each section, as shown in Figure 9.
According to the test analysis results, the Cvsic model can be used to study tunnel creep. e creep parameters G K , G M , η K , and η M of the Cvsic model are inversed intelligently. Table 2 shows the range of creep parameters to be inversed. e orthogonal parameter scheme and uniform design parameters are shown in Tables 3 and 4. Table 3 is the  training data, and Table 4 is the verification data. Using FLAC 3D 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. e monitoring data of section LK49 + 900 are compared with the inversion results.

Back Analysis Results.
e creep parameters of tunnel surrounding rock are inversed by the GP-DE algorithm. Setting relevant parameters as follows: the population size NP � 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. e results of the uniformly parameters schemes are used for verification.
According to the inversion results, FLAC 3D software is used again for forwarding analysis, and the calculated displacement of section LK49 + 900 is compared with the measured displacement, in Figure 10. e 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.
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 Table 7.
As shown in Table 7, 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. e fitness value of the GP-DE algorithm group under different iteration times is shown in Figure 11. It can be seen Table 3: e orthogonal parameter sample.  Table 4: e uniform design parameter sample.  from the figure that the fitness value tends to converge as the number of iterations increases, and the distribution in the space gradually decreases.

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 Figure 12. It can be seen from the figure that the prestressed anchor rod plays a small role, while the grouting reinforcement plays an obvious role. erefore, it is recommended to adopt the scheme of initial lining + pipe shed + advanced grouting anchor rod. Figure 13 shows the creep deformation curve at the tunnel crown. e 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.
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. When t⟶∞, 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. erefore, it is necessary to speed up the construction of secondary lining or reduce the thickness of primary lining. Figure 14 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. is shows that it is very necessary to select reasonable super parameters to ensure the inversion accuracy of the GP model.

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. is paper analyzes the influence of mutation factor F, crossover factor CR, population size NP, and different difference strategies on the DE algorithm and selects reasonable parameters. Figure 15 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. e results show that the convergence performance of the DE algorithm is better when CR � 0.9 and F � 0.6 are selected.    Figure 16 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.

Influence of Differentiation Strategy on Inversion
Accuracy. In the DE algorithm, the different differential strategies are as follows:

DE best/2
: v i,g � x best,g + F x r1,g − x r2,g + x r3,g − x r4,g , DE rand − to − best/2 : v i,g � x r1,g + F x best,g − x r2,g + x r3,g − x r4,g , where x best,g is the best individual and x r1 , x r2 , x r3 , x r4 , and x r5 are random individuals. e iterative curves of different difference strategies are shown in Figure 17. 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.

Conclusions
(1) According to the creep test of gneiss, the creep model and parameter range of the tunnel are determined. e results show that the whole creep process of gneiss presents four stages: instantaneous elastic deformation, attenuation creep, stable creep, and sudden failure. e 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. e algorithm performance test and parameter analysis are carried out. Selecting CR � 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) e stability analysis of the tunnel and the selection of a reasonable construction plan are carried out. e 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.
Data Availability e datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Advances in Materials Science and Engineering