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.
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 are m observations, the inversion optimization equation is as follows: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 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, , and is the variance.
The prior distribution of the output value of the training sample iswhere 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 values y∗ at test points x∗. The joint distribution of target value and function value can be written as follows:where K (, ) is the covariance matrix of and K (X, ) is the order n × 1 covariance matrix of X and .
When the input value of the test sample and the training set D are known, the output value of the test sample can be calculated through the Gaussian process:where are the expectation and variance of ; ; 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):where xp and xq are the learning samples or prediction samples, is the distance between xp and xq, is the standard deviation of the noise, is the local correlation, and is a sign function.
Hyperparameters and affect the GP training effect, which is expressed as follows:where is the estimated output data, Yh is the actual output of the test sample, and 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:where and and are the upper and lower limits of a single variable.
2.3.2. Mutation Operation
Three individuals , , and are randomly selected from the population for mutation operation, as shown in Figure 1. The variation vector is as follows:where , , and are unequal integers between [1, NP] and F∈[0,1] is the variation factor.
2.3.3. Crossover Operation
A new vector is generated by crossing the target vector with the variation vector as follows: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: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: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 of GP 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 by DE 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, the DE 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.
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 Figure 3.
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 Figure 4.
(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 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, , until the specimen is damaged. Figure 6 shows the results of creep tests.
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 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 aswhere is the Kelvin body strain, is the Maxwell body strain, and is the Mohr–Coulomb body strain.
For the viscoelastic unit, the constitutive laws for the Kelvin and Maxwell units are formulated bywhere is the stress tensors, and are the shear modulus and the viscosity of Kelvin, and and 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 Figure 8. 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.
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.
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 Figure 9.
According to the test analysis results, the Cvsic model can be used to study tunnel creep. The creep parameters GK, GM, ηK, and ηM of the Cvsic model are inversed intelligently. Table 2 shows the range of creep parameters to be inversed.
The 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 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.
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 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. The results of the uniformly parameters schemes are used for verification. Table 6 shows the creep parameters that are obtained by inversion.
According 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.
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.
The fitness value of the GP-DE algorithm group under different iteration times is shown in Figure 11. 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.
(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 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. Therefore, it is recommended to adopt the scheme of initial lining + pipe shed + advanced grouting anchor rod.
(a)
(b)
(c)
(d)
Figure 13 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.
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. 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
Figure 14 shows the influence of super parameters on the inversion accuracy of the GP model. It can be seen from the figure that and play a key role in the inversion accuracy. When and , 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.
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 factor F, crossover factor CR, population size NP, and different difference strategies on the DE algorithm and selects reasonable parameters.
6.2.1. Influence of F and CR on Inversion Accuracy
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. The results show that the convergence performance of the DE algorithm is better when CR = 0.9 and F = 0.6 are selected.
(a)
(b)
6.2.2. Influence of NP on Inversion Accuracy
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.
6.2.3. Influence of Differentiation Strategy on Inversion Accuracy
In the DE algorithm, the different differential strategies are as follows:where is the best individual and are random individuals.
The 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.
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. 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)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.Data Availability
The datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Nos. 51678101 and 52078093), the Basal Research Funds for the Central Universities (No. 3132014326), and the LiaoNing Revitalization Talents Program (No. XLYC1905015).