A novel method for identifying geomechanical parameters of rock masses based on a PSO and improved GPR hybrid algorithm

In view of the shortcomings of existing artificial neural network (ANN) and support vector regression (SVR) in the application of three-dimensional displacement back analysis, Gaussian process regression (GPR) algorithm is introduced to make up for the shortcomings of existing intelligent inversion methods. In order to improve the generality of the standard GPR algorithm with single kernel function, an improved Gaussian process regression (IGPR) algorithm with combined kernel function is proposed by adding two single kernel functions. In addition, in the training process of IGPR model, the particle swarm optimization (PSO) is combined with the IGPR model (PSO-IGPR) to optimize the parameters of the IGPR model. After the IGPR model can accurately map the relationship between geomechanical parameters and rock mass deformation, the PSO algorithm is directly used to search the best geomechanical parameters to match the deformation calculated by igpr model with the measured deformation of rock mass. The application case of Beikou tunnel shows that the combined kernel function GPR has higher identification accuracy than the single kernel function GPR and SVR model, the IGPR model with automatic correlation determination (ARD) kernel function can obtain higher identification accuracy than the IGPR model with isotropic (ISO) kernel function, and the PSO-IGPR hybrid model based on ARD kernel function has the highest identification accuracy. Therefore, this paper proposes a displacement back analysis method of the PSO-IGPR hybrid algorithm based on ARD kernel function, which can be used to identify the geomechanical parameters of rock mass and solve other engineering problems.

A novel method for identifying geomechanical parameters of rock masses based on a PSO and improved GPR hybrid algorithm Hanghang Yan 1 , Kaiyun Liu 1* , Chong Xu 2 & Wenbo Zheng 3 In view of the shortcomings of existing artificial neural network (ANN) and support vector regression (SVR) in the application of three-dimensional displacement back analysis, Gaussian process regression (GPR) algorithm is introduced to make up for the shortcomings of existing intelligent inversion methods. In order to improve the generality of the standard GPR algorithm with single kernel function, an improved Gaussian process regression (IGPR) algorithm with combined kernel function is proposed by adding two single kernel functions. In addition, in the training process of IGPR model, the particle swarm optimization (PSO) is combined with the IGPR model (PSO-IGPR) to optimize the parameters of the IGPR model. After the IGPR model can accurately map the relationship between geomechanical parameters and rock mass deformation, the PSO algorithm is directly used to search the best geomechanical parameters to match the deformation calculated by igpr model with the measured deformation of rock mass. The application case of Beikou tunnel shows that the combined kernel function GPR has higher identification accuracy than the single kernel function GPR and SVR model, the IGPR model with automatic correlation determination (ARD) kernel function can obtain higher identification accuracy than the IGPR model with isotropic (ISO) kernel function, and the PSO-IGPR hybrid model based on ARD kernel function has the highest identification accuracy. Therefore, this paper proposes a displacement back analysis method of the PSO-IGPR hybrid algorithm based on ARD kernel function, which can be used to identify the geomechanical parameters of rock mass and solve other engineering problems.
Back analysis techniques are often used in geotechnical engineering to determine unknown geomechanical parameters of rock masses by using field measurements of displacements, strains, and stresses during excavation or construction phases [1][2][3][4][5] . The displacement of rock masses induced by excavation can be measured reliably with relative ease. Therefore, displacement back analysis is the most common and effective method for identifying geomechanical parameters in rock engineering [6][7][8][9][10][11][12][13] . There are two types of displacement back analysis, namely, the inverse solving method and the direct optimization method. Compared with the inverse solving method, the direct optimization inversion method is not limited by the constitutive relationship of materials and has better applicability.
Two-dimensional numerical analysis is typically used in displacement back analysis to shorten the calculation time to meet the needs of engineering construction. However, in many cases, two-dimensional numerical simulations cannot account for the time and spatial effects of engineering construction. Several studies have introduced three-dimensional (3D) numerical simulations into displacement back analysis [12][13][14] , however, the computation time of 3D displacement back analysis is costly. Moreover, the parameters of rock masses exhibit considerable uncertainty, variability, and nonlinearity. Based on the advantages of fast calculation, self-fault tolerance, self-adaptation and considering high uncertainty, artificial intelligence technology is also widely used in the field of geotechnical engineering [15][16][17][18][19][20] .Thus, the intelligent back analysis method was introduced by incorporating

Improved Gaussian process regression (IGPR) algorithms
The GP is a type of machine learning method based on the Gaussian stochastic process and Bayesian learning theory. In the theory of statistics, the GP is a random process, and the distribution of any finite set of variables is a Gaussian distribution. That is, the joint probability distribution of any set of random variables {x i ∈ X, i = 1, …, n} and its corresponding process state {Y(x 1 ), …, Y(x n )} follows an n-dimensional Gaussian distribution. All statistical characteristics of a GP are determined by its mean value and covariance functions 35,40 where x, x ′ ∈ X are both arbitrary random variables and x i is a multi-dimensional input vector. The GP can be defined as f (x) ∼ GP(µ(x), C(x, x ′ )).
A set of n observation data D = {(x i , t i ), i = 1, · · · , n} is considered the training set of the Gaussian model. The observed target t is corrupted by noise, and the difference in the actual output value is ε; thus, the Gaussian noise model can be written as follows: where t i denotes the output scalar, ε denotes independent random variables with ε ∼ N(0, σ 2 n ) (N is the normal distribution) and σ n denotes the variance of the noise. Under the framework of Bayesian linear regression f(x) = Ψ(x) T w, the prior distribution of the target value t of Eq. (1) is obtained as: (2) t i = f (x i ) + ε i , i = 1, ..., n www.nature.com/scientificreports/ where I denotes the unit matrix. According to Eq. (3), the joint Gaussian prior distribution of the training sample output t and a test sample output t* can be deduced as follows: where C(X, X) is a n × n positive definite covariance matrix. Its arbitrary c ij measures the correlation between x i and x j . C(X, x * ) is a n × 1 covariance matrix of the tested x * and the input points of X, and C(x * , x * ) denotes the covariance of test point x * . For a given test point x * and training set D, the goal of Bayesian probability prediction is to calculate p(t * |x * , D ) . According to Bayesian posterior probability, Eq. (5) can be obtained: The expectation and variance of t * are as expressed by Williams 35 : Because the covariance function of the GP method in the finite set of input requirements, which is a symmetric function satisfying the Mercer condition, is positive, the covariance function is equivalent to the kernel function. Equations (6) and (7) can be rewritten as: The mean value of the prediction is a linear combination of kernel function K. The nonlinear relational data can be mapped into the feature space and transformed into a linear relation; therefore, the complex nonlinear problem can be transformed into a linear problem.
Kernel functions can be divided into two categories: isotropic (ISO) kernel functions and automatic relevance determination (ARD) kernel functions. In this paper, the following two types of ARD kernel functions are applied: (1) Square exponential (SE) kernel function (2) Rational quadratic (RQ) kernel function where {M} = diag(ℓ −2 ) denotes the diagonal matrix of the hyperparameters and ℓ represents the hyperparameter, which determines the relevance between the input variable and the output variable. A larger value of ℓ indicates a lower degree of relevance between the input and output variables. In addition, σ 2 f denotes the signal variance of the kernel function and can be used to control the degree of local relevance, α represents the shape parameter of the kernel function, σ 2 n is the noise variance, and δ ij is the Kronecker symbol: The ARD-type kernel function and the ISO-type kernel function are the same in form; the difference between the two types of functions is the dimension of hyperparameter ℓ . The dimension of ℓ is the same as that of the input variable x in the ARD-type kernel function. In other words, if the input variable x is a d-dimensional vector, ℓ is also a d-dimensional vector. However, for the ISO-type kernel function, the dimension of ℓ cannot be changed with the input variable dimension. Instead, it is always a one-dimensional scalar, which means that the relevance between all components of the input variable x and the output variable is the same 35,40 .
The hyperparameters of the kernel function have a considerable influence on the learning and prediction results. In the GPR algorithm, the optimal hyperparameters can be adaptively obtained by using the maximum likelihood method to establish a log marginal likelihood function of training samples and then obtaining the www.nature.com/scientificreports/ partial derivative of the hyperparameters. Finally, the optimum solution of the hyperparameters can be calculated using the conjugate gradient optimization method as follows 35,40 : where θ denotes the vector that contains all of the hyperparameters.
The following combination kernel (CK) function can be obtained by adding Eqs. (10) and (11) to improve the generalizability of GPR 40 : The improved Gaussian process regression (IGPR) algorithm is a GPR algorithm that uses the CK function shown in Eq. (15).

PSO-IGPR hybrid algorithm for geomechanical parameter identification
IGPR model for the non-linear relationship between geomechanical parameters and rock displacement. A non-linear relationship of geomechanical parameters X (x 1 ,x 2 ,…,x n ) with the displacement of the rock mass at key location points Y (y 1 , y 2 …, y m ) is represented by IGPR(X) as.
where X = (x 1 ,x 2 ,…,x n ) represents the input variables, including the geomechanical parameters of the rock mass that need to be identified, and Y = (y 1 , y 2 …, y m ) is the displacement of key location points.
To obtain IGPR(X), a training process based on the known data set is needed. The best way is to train the IGPR algorithm and establish an optimal IGPR model to fit the non-linear relation between geomechanical parameters and displacement of rock mass by using the displacement data and the corresponding geomechanical parameters of rock mass measured in situ. However, such data are generally not available. To solve this problem, numerical simulations of rock masses are usually used. Under the premise of the given geomechanical parameter range of the rock mass, the corresponding rock mass displacements under different geomechanical parameters are obtained by numerical simulations. Numerical simulations in this paper are conducted with FLAC 3D . The IGPR algorithm is trained with these data, and the intelligent non-linear mapping model IGPR(X) between geomechanical parameters and displacements of rock mass can be established. The intelligent IGPR(X) model can replace numerical calculations for displacement back analysis.
To establish the IGPR model, the objective function,g(θ ) , is defined as follows: where θ are the hyperparameters of the IGPR algorithm. FLAC(X, d) i represents the displacement of the ith testing sample calculated by FLAC simulations under the given geomechanical parameters X of the rock mass and the value of distance d between the tunnel face and the monitoring section. IGPR(θ |X,d) i indicates the displacement of the ith testing sample calculated by the IGPR model under the given geomechanical parameters X of the rock mass and the value of distance d but with unknown hyperparameters θ of the IGPR model. S denotes the number of testing samples. Equation (10) shows that the SE kernel function has three hyperparameters: σ f , ℓ , and σ n . As shown in Eq. (11), the RQ kernel function has four hyperparameters: σ f , ℓ, α , and σ n . The combination kernel function shown in Eq. (15) has 6 hyperparameters: σ SE f , ℓ SE , σ RQ f , ℓ RQ , α , and σ n . Because the generalization performance of the GPR algorithm is directly affected by the selection of each hyperparameter value, the value of each hyperparameter of GPR must be optimized, which is a typical multiparameter combination optimization problem. The conjugate gradient method can be used to obtain the optimal GPR hyperparameters (as shown in Eqs. 13 and 14), but it is challenging to determine the initial value and number of iterations for finding a globally optimal solution.
In recent years, swarm intelligence bionic optimization algorithm has been highly valued and widely used 52,53 . As a bionic global optimization algorithm, the objective function of PSO does not need to be differentiable, and the optimization results do not depend on the initial values and implicit parallel searches. PSO has been widely used to solve multiparameter, multivariable, and multiobjective optimization problems. These advantages of PSO compensate for the shortcomings of the conjugate gradient method. Compared with the genetic algorithm, the PSO algorithm is simpler in theory and easier to implement in programs. In this paper, PSO replaces the conjugate gradient method and automatically searches for the optimal hyperparameters during the training process www.nature.com/scientificreports/ of the IGPR algorithm to establish the non-linear intelligent mapping model IGPR(X) between geomechanical parameters and displacements of rock mass.

PSO coupled with the IGPR algorithm (PSO-IGPR). PSO is based on the following model and is used
to solve optimization problems [54][55][56] . In PSO, each solution of the optimization problem is a bird, or "particle", in the search space. Each particle has a fitness value determined by the optimized function and a speed that determines the direction and distance of flight. The particles will follow the current optimal particle to search the solution space.
In PSO, a group of random particles (random solutions) is initialized. Then, the optimal solution is iteratively determined. In each iteration, particles are updated by tracking two extremes. The first is the optimal solution found by the particle itself. This solution is called the individual extremum pbest, and the other extremum is the optimal solution found by the entire population. This extremum is the global extremum gbest. After finding these two optimal values, the speed and position of the particle are updated according to the following equations: where P k i and P k+1 i are the n-dimensional vectors that represent the position of the i-th particle in the search space at iterations k and k + 1, respectively;v k i and v k+1 i denote its speed at iterations k and k + 1, respectively; pbest k i and gbest k are defined as the individual extremum of the i-th particle and the global extremum at iteration k, respectively; r 1 and r 2 are random numbers between 0 and 1; and c 1 and c 2 represent learning factors. In most cases, c 1 = c 2 = 2.
The individual extremum pbest k+1 i of the i-th particle and the global extremum gbest k+1 at iteration k + 1 are updated according to the following equation: where f represents the fitness function of the PSO algorithm.
IGPR training is a process of continual readjustment between the hyper-parameters to reduce the prediction error of the IGPR model to a pre-set minimum or stop at a pre-set iteration step. Then, the testing samples are input to the trained IGPR model, and the forecasting results are obtained.
There, we use the PSO algorithm, which is characterized by fast global optimization, to search for the optimal hyper-parameters of the IGPR model and avoid trial calculations. The steps involved in implementing the PSO and IGPR (PSO-IGPR) coupled algorithm are as follows.
(1) Initialization of PSO network parameters. The particle population size, number of iterations, initial particle velocities and initial locations of particles are determined. Each particle vector represents an IGPR model. Different models correspond to different parameters of the IGPR algorithm. (2) The IGPR algorithm learns knowledge from the learning samples and predicts the displacement value of each testing sample. Then, the PSO algorithm calculates the individual fitness value f i of each particle. (3) Compare the fitness f i calculated in step (2) with the best fitness f(pbest i ) of the particle in the iteration history. If f(pbest i ) > f i , f i is used to replace f(pbest i ) from the previous round, and new particles replace the particles from the previous round. (4) Compare the best fitness f(pbest i ) of each particle with the best fitness f(gbest) of all the particles. If f(pbest i ) < f(gbest), then replace the original global best fitness value f(gbest) with the best fitness value f(pbest i ) of the particle and save the current state of the particle. (5) When the calculation reaches the pre-set number of iterations, the calculation is finished, and the particle with the minimum fitness value is returned to obtain the optimal solution. Otherwise, a new round of iterations is needed to update the positions and speeds of the particles and generate new particles. Then, the process returns to step (2) until the maximum number of iterations is reached and the calculation is finished. Figure 1 shows the whole flow chart of the back-analysis of displacement based on the PSO-IGPR hybrid algorithm. Part I of the flow chart represents the modelling stage, i.e., the establishment stage of an IGPR intelligent model that can best fit the non-linear relationship between geomechanical parameters and displacement of rock mass. As shown in Fig. 1, the calculation steps of this stage are as follows: Step 1: Different geomechanical parameters of the rock mass are used for numerical simulation of tunnel construction. The displacements of key points under different geomechanical parameters and different distances between the tunnel face and the monitoring section are recorded, which are used as the training samples of the IGPR algorithm. The training samples (X i , Y i )(i = 1, 2, · · · , k; X ∈ R n , Y ∈ R m ) are divided into two parts: the learning samples for the IGPR algorithm to obtain knowledge and the testing samples used to examine the generalizability of the IGPR model. www.nature.com/scientificreports/ Step 2: The PSO is initialized, and the initial population of the hyperparameters θ of the IGPR algorithm is randomly generated with N p individuals. Each individual in the population determines a set of hyperparameters θ of the IGPR algorithm, that is, an IGPR model; the generation counter g is set to 1.
Step 3: The IGPR algorithm reads (X i , Y i )(i = 1, 2, · · · , k − s; X ∈ R n , Y ∈ R m ) in the learning samples and the input variables X i (i = k − 1 + s, · · · , k) of the testing samples. The hyperparameters θ of each IGPR model www.nature.com/scientificreports/ in the initial population are read, and then, the learning samples are trained by the IGPR algorithm to obtain knowledge and predict the testing samples.
Step 4: The prediction result of each individual of the testing samples is delivered to the following fitness function of the PSO to calculate the fitness of each individual.
i indicates the minimal relative error between the IGPR model calculation displacement and the FLAC 3D calculation displacement of the s testing samples.
According to Eq. (22), f (θ ) ≥ 1 . As the fitness value f (θ ) decreases, the calculated results of the IGPR model more closely approach the simulated values by FLAC. When f (θ ) = 1 , the calculated results of the IGPR model are equal to the simulated values by FLAC.
Step 5: Compare the fitness f i calculated in previous step with the best fitness f(pbest i ) of the particle in the iteration history. If f(pbest i ) > f i , f i is used to replace f(pbest i ) from the previous round, and new particles replace the particles from the previous round. Then, compare the best fitness f(pbest i ) of each particle with the best fitness f(gbest) of all the particles. If f(pbest i ) < f(gbest), then replace the original global best fitness value f(gbest) with the best fitness value f(pbest i ) of the particle and save the current state of the particle.
Step 6: Whether the prespecified number of evolution generations of the PSO was reached is determined. If reached, then the calculation stops and the individual IGPR model with the minimal fitness f (θ * ) is returned; the optimal hyperparameters θ * of the IGPR model have been obtained. Otherwise, a new round of iterations is needed to update the positions and speeds of the particles and generate new particles. Then, the process returns to step 3 until the maximum number of iterations is reached and the calculation is finished. So far, the optimal hyperparameters θ * of the IGPR model have been obtained.
To date, the best non-linear IGPR intelligent model with the optimal hyperparameters θ * mapping the relationship between geomechanical parameters and displacement of rock mass has been established. We can use this best IGPR model to replace the traditional 3D numerical simulation to carry out displacement back analysis, which can greatly reduce the calculation time to meet the construction needs.

Back analysis for identifying geomechanical parameters of a rock mass using PSO-IGPR.
The best IGPR intelligent model with optimal hyperparameters θ * that can map the nonlinear relationship between the geomechanical parameters and the displacement of the rock mass is established in Sect. PSO coupled with the IGPR algorithm (PSO-IGPR). For intelligent displacement back analysis, the objective function of the geomechanical parameter identification stage for the rock mass is shown in the following equation: where u i represents the measured displacement values of the ith monitoring point (crown subsidence or horizontal convergence) and IGPR(X θ * ,d) i indicates the displacement values of the ith monitoring point calculated by the best IGPR model under the optimal hyperparameters θ * and the value of distance d but with unknown geomechanical parameters X of the rock mass. In addition, p denotes the number of monitoring points (p = 1 here).
As shown in Eq. (23), iteration optimization is used in the proposed displacement back analysis method to search the optimal geomechanical parameters X of the rock mass by using PSO, which can ensure that the displacement calculated by the best IGPR model is as close to the measured displacement as possible; these parameters are considered the "real geomechanical parameters" of the rock mass, thus completing parameter identification.
Part II of Fig. 1 represents the identification stage, i.e., the identification stage of geomechanical parameters of rock mass based on PSO-IGPR hybrid algorithm. The calculation steps of this stage are as follows: Step 1: The PSO step is initialized, and then, the initial population of the geomechanical parameters X of the rock mass is randomly generated with N p individuals. Each individual in the population determines a set of geomechanical parameters X of the rock mass; the generation counter g is set to 1.
Step 2: The best IGPR model reads in each rock mass geomechanics parameter individual of the initial population and the specified distance d between the tunnel face and the monitoring section. The key point displacements can be predicted by the excellent generalization performance of the IGPR algorithm.
Step 3: The PSO step reads in the corresponding measured displacements of each key point, calculates the fitness values of each rock mass geomechanical parameter according to the fitness function defined below, and evaluates the individual fitness of each rock mass geomechanical parameter.
where arg min[ |ui−IGPR(X|θ * ,d) i | u i × 100%] p i indicates the minimal relative error between the displacements predicted by the best IGPR model and the measured displacements of the p monitoring points.
According to Eq. (24), f (X) ≥ 1 . The smaller the fitness value is, the closer the displacements predicted by the best IGPR model are to the measured displacements. When f (X) = 1 , the calculated displacements of the best IGPR model are equal to the measured displacements. www.nature.com/scientificreports/ Step 4: Is the prespecified generation number of the PSO reached? If reached, then the calculation stops and the geomechanical parameter individual X * with the minimal fitness f (X * ) is returned; the optimal geomechanical parameters X * of the rock mass have been obtained. If not, the calculation procedure proceeds to the next step.
Step 5: Perform similar operations similar to steps 5 and 6 of the above training process until the prespecified number of evolution generations of the PSO is reached and the optimal particle individual (G best ) is returned, which is the geomechanical parameter identification value of rock mass.
The geomechanical parameter identification of rock masses based on the PSO-IGPR hybrid algorithm has been completed. Generally, this method consists of two stages: STAGE 1: This stage is the modelling stage, i.e., training the IGPR algorithm with samples obtained from numerical simulation. In this training process, the PSO algorithm is used to automatically search the optimal IGPR hyperparameters θ * , which can bring the displacement predicted by the IGPR algorithm closest to the displacement calculated by numerical simulation, and to establish the IGPR model, which can best fit the nonlinear relationship between geomechanical parameters and the displacement of the rock mass.

STAGE 2:
The second stage is the identification stage. That is, PSO is used to automatically search the optimal geomechanical parameters X * within the geomechanical parameter range of the rock mass, which can bring the displacement predicted by the best IGPR model established in the previous stage closest to the measured displacement, and complete the geomechanical parameter identification of the rock mass. The quality rating result of the rock masses around the tunnel based on the basic quality (BQ) classification method of China is shown in Table 1 based on the initial geological survey.

IGPR model of rock mass deformation around the Beikou Tunnel. Three-dimensional numerical
simulations. According to the flow chart described in the previous section, first, numerical simulations should be used to form samples for IGPR training. In this paper, a numerical tool, FLAC3D, is used to carry out 3D numerical simulations to obtain the deformation of rock masses in Beikou tunnel construction.
The chainage numbers at the entrance of the Beikou tunnel for the left and right lines are LK59 + 130 and LK59 + 255, respectively. The entrances of the left and right lines are 125 m apart in the direction of excavation. In addition, the distance from the axis of the left tunnel to the axis of the right tunnel is approximately 71 m, and the maximum excavation width of the tunnel is approximately 11 m. Therefore, the construction of the left tunnel has limited influence on the right tunnel and vice versa. This study modelled only the construction of the left line to shorten the computation time for numerical simulations.    Because there are no in situ stress data, the ground stress in the vertical direction is taken as the weight of the overlying strata, and the horizontal stress in two directions is assumed to be equal to the vertical stress multiplied by the lateral pressure coefficient.
The excavation was performed by the upper and lower bench excavation method; the excavation of the lower bench of the tunnel face was divided into two parts: left and right. Compared to the right part of the lower bench, the excavation footage of the left part was 15 m ahead. Compared to the left part of the lower bench, the excavation footage of the upper bench was 15 m ahead. Each cyclic excavation was 1.5 m. In this method, preliminary bracing was used for the tunnel face. The initial support parameters of the tunnel excavation are shown in Table 2. The cable structural element and the shell element in the FLAC software are used to simulate the anchor bolt and primary shotcrete, respectively.
Thirty test schemes are designed using the uniform test method for the mechanical parameters of the rock mass and the lateral pressure coefficient. The construction of the left line entrance is simulated by FLAC3D software. The distances d between the tunnel face and monitoring section LK59 + 150 are 6, 7.5, and 9 m, respectively. The calculated results are extracted randomly to compose 36 learning samples and 5 testing samples for training the IGPR algorithm, as shown in Table 3.
IGPR model for mapping geomechanical parameters to displacement of rock mass. By training on the samples shown in Table 3, PSO is used to search for the optimal hyperparameters θ * of the IGPR model that can make the   www.nature.com/scientificreports/ displacement calculated by the IGPR model closest to the displacement calculated by the FLAC3D simulations during the sample training process. Therefore, the best IGPR intelligent model with the optimal hyperparameters θ * that maps the nonlinear relation between the geomechanical parameters and the displacement of the rock mass is established. Three different algorithms with the same inversion process are used to execute the displacement back analysis for comparison purposes. These algorithms are the PSO and standard GPR coupled (PSO-GPR) algorithm with a single kernel function, the PSO-IGPR algorithm and the PSO and SVR coupled (PSO-SVR) algorithm. The PSO parameters of the above three identification methods are identical. The CK function shown in Eq. (15) Table 3. Training sample set. Here and elsewhere, the symbol "d" represents the distance from the tunnel face to monitoring section LK59 + 150. www.nature.com/scientificreports/ where r is a positive integer that represents the order of the polynomial. There are two other parameters of the SVR algorithm: the penalty parameter C and fitting error ε 56,57 . The search intervals of SVR parameters C and ε are [0, 5000] and [0, 1], respectively, and the kernel parameter k is 1, 2 or 3. The population size of the PSO is 20, and the number of evolution generations is 500. The output variable of the standard GPR and standard SVR algorithms can only be a one-dimensional scalar 32,33 ; therefore, two intelligent models mapping the relationship of the crown subsidence and the horizontal convergence to the geomechanical parameters of the rock mass are established.
Rasmussen 34 developed the open source program of the GPR algorithm based on MATLAB language. On this basis, the code of the PSO-IGPR hybrid algorithm is programmed in MATLAB according to the flowchart shown in Fig. 1 and Eqs. (22) and (24) and then applied to the training sample set shown in Table 3. The best GPR, IGPR and SVR models with the optimal hyperparameters and parameters searched by the PSO on the basis of the steps described in Section "PSO coupled with the IGPR algorithm (PSO-IGPR)" are shown in Tables 4,  5 and 6, respectively.
Rock mass parameter identification for the Beikou tunnel using PSO. The geomechanical parameters of the rock mass to be identified are the deformation modulus E, cohesion C, internal friction angle Φ, horizontal lateral pressure coefficient λ and Poisson's ratio μ. According to the value range of the geomechanical parameters of the different grades of rock mass specified in the Standard for engineering classification of rock masses 58 , the search intervals of the geomechanical parameters are provided in Table 7.  www.nature.com/scientificreports/ The best intelligent GPR, IGPR and SVR models (Tables 4, 5 and 6) that can map the nonlinear relationship between the displacement and geomechanical parameters of the rock mass around the Beikou tunnel are established in Section "IGPR model of rock mass deformation around the Beikou tunnel". The measured horizontal convergence and crown subsidence values are input into these intelligent models. The distance between the tunnel face and monitoring section LK59 + 150 is 6 m. The number of evolution generations and the population size of the PSO are 500 and 20, respectively. The other parameters of the PSO remain unchanged. The fitness function of the PSO is adopted the Eq. (24) during the geomechanical parameter identification stage.
The identification results obtained using the PSO are shown in Table 8. Then, the identified parameter values shown in Table 8 are input into the intelligent models shown in Tables 4, 5 and 6; the horizontal convergence and crown subsidence values in monitoring section LK59 + 150 can be forecast for the following excavation. The results are shown in Tables 9 and 10.
Analysis of calculation results. Table 9 illustrates that the maximum prediction error of GPR is slightly larger than that of SVR, but the average relative error is less than that of the SVR algorithm when using identification parameters to predict the surrounding rock deformation for five successive excavation steps. This result Table 7. Searching intervals of the geomechanical parameters to be identified.  Table 8. Geomechanical parameters identified by the displacement of the rock mass around the tunnel of section LK59 + 150 (d = 6.0 m).  www.nature.com/scientificreports/ means that all of the GPR algorithms have higher identification accuracy than the SVR algorithm. A similar result can be obtained from the crown subsidence results forecasted using the identified parameters for subsequent excavation in Table 10. Overall, the prediction precisions of the three types of GPR algorithms are better than that of the SVR algorithm regardless of whether an ARD-type or ISO-type kernel function is used. Moreover, the GPR algorithm with an ARD-type kernel function is better than the GPR algorithm with an ISOtype kernel function for the same form of the kernel. When the same type of kernel function is used, the IGPR algorithm has higher identification accuracy than the GPR algorithm with a single kernel function. The SVR algorithm has the lowest identification accuracy. The prediction accuracy of the IGPR algorithm with an ARD-type combination kernel function is exceptionally high. The maximal prediction relative error of the surrounding rock deformation using the identification parameters for five successive excavation steps is only 10.26%; the average prediction relative errors of the horizontal convergence and the crown subsidence are 4.3% and 7.39%, respectively, which fully satisfies the needs of tunnel construction. For the IGPR algorithm with an ISO-type combination kernel function, the maximal prediction relative error of the surrounding rock deformation using the identification parameters for five successive excavation steps is 38.46%, and the average prediction relative errors of the horizontal convergence and crown subsidence are 6.34% and 25.97%, respectively. Therefore, the IGPR algorithm with an ARD-type combination kernel function has higher identification accuracy than the IGPR algorithm with an ISO-type combination kernel function.

Identification value of parameters E (GPa) C (MPa) Φ (°) λ μ
Discussion. Robustness of CKGPR model. All three machine learning algorithms, namely, SVR, standard GPR and IGPR, are comparatively presented briefly for intelligent displacement back analysis in this paper. In contrast, SVR does not perform as well as standard GPR and IGPR under less raw data conditions. Since GPR derives the mapping from probability theory with prior and covariance functions, the prediction of GPR is a probabilistic distribution; GPR can apply the mean of the distribution as point predictions to avoid robust point predictions like that in SVR, which leads to GPR performing slightly better than SVR in the practical application of displacement back analysis.
In order to test the robustness of the IGPR (ARD type) model established in this paper, the second test sample in Table 3 is taken as the model input benchmark sample and adopted the super parameters of IGPR model shown in Table 5. The input parameters of this sample are changed step by step, the corresponding prediction output of the model is checked, and compared with the sample output value (benchmark value) to observe the robustness of the model. The change range of each input parameter is increased or decreased by 10%, the calculation results are shown in the following table: It can be seen from the calculation results in Table 11 that when the input parameters of IGPR model increase by 10% one by one, compared with the benchmark value, the maximum relative error of horizontal convergence prediction of the model is 6.5%, and the average relative error of prediction is only 4.03%; The maximum relative error of the model is 11.5%, and the average relative error is only 5.12% When the input parameters of IGPR model are reduced by 10% one by one, compared with the benchmark value, the maximum relative error of horizontal convergence prediction of the model is 16.8%, and the average relative error of prediction is only 5.13%; The maximum relative error of the model is 15.4%, and the average relative error is only 8.96% It can be seen that the ard type igpr model established in this paper has good robustness.
Computational efficiency analysis. For the GPR algorithm, the relevant measure hyperparameter ℓ of the ARDtype kernel function is more than that of the ISO-type kernel function. Therefore, the generalization performance of the GPR algorithm increases with an increasing number of parameters. By comparison, the number of hyperparameters of the IGPR algorithm with an ARD-type combination kernel function is the largest, and its generalization performance is the best among the three algorithms.
Computation time is monitored by defining a clock function in the program; the computation times for different algorithms in the sample training and parameter identification phases are shown in Table 12 (personal computer configuration for a dual-core processor, 3.1 GB and 4 GB of memory). www.nature.com/scientificreports/ Table 12 demonstrates that in the sample training stage, the computation time of the SVR algorithm is approximately 2-5 times that of GPR; in the parameter identification stage, the computation time of SVR is approximately 4-8 times that of GPR. The GPR algorithm is far more efficient than the SVR algorithm.
Limitations and future studies. For multi-parameter direct optimization inversion, the biggest problem is the non-uniqueness of the solution, which brings great difficulties to practical application. The feasible solution is to limit the search interval of each parameter to be inversed as much as possible, that is, to enhance the priori of the parameters, which is the next research focus of multi-parameter displacement back analysis.
The amount of data is a critical factor that affects the analysis results for any machine learning algorithm. With limited data, it is difficult to obtain an accurate and robust model. Therefore, the trained model should be used carefully. In addition, the prediction result of the learning machine also depends on the quality of the data and its generalization performance, which is directly affected by the value of the parameters. In this paper, PSO is a coupled learning machine to automatically search for the best parameters of the learning machine in the training process, and the effect is excellent. However, it is necessary to further study the parameter selection theory of the GP.

Conclusions
The following conclusions are obtained through the above comparative study: (1) The GPR algorithm can be used to identify the geomechanical parameters of rock masses. This can improve the accuracy of the identification compared to that when using the SVR algorithm. (2) The IGPR algorithm can improve the generalization performance and the identification precision of the GPR algorithm with a single kernel function. The average relative prediction error of the IGPR algorithm with an ARD-type kernel function is approximately 27-74% lower than that of the GPR algorithm with a single kernel function; that of the IGPR algorithm with an ISO-type kernel function is approximately  www.nature.com/scientificreports/ 11-52% lower than that of the GPR algorithm with a single kernel function. In comparison, the ARD-type combination kernel GPR has the largest number of hyperparameters, but the identification accuracy is significantly improved without a considerable increase in computation time. (3) The maximal relative error of each excavation step and the maximal average relative error of all excavation steps forecasted by the IGPR algorithm based on the ARD-type kernel function are only 10.26% and 7.39%, respectively, over five excavation steps. This result shows that the novel 3D displacement back analysis method based on the PSO-IPGR hybrid algorithm with the ARD-type kernel function proposed in this paper has very high accuracy when identifying the geomechanical parameters of the rock mass and can be applied for parameter identification in geotechnical engineering applications.

Method
The geomechanical parameters of rock mass are the key parameters directly related to the stress and strain calculation of rock mass engineering, and the geomechanical parameters of rock mass are usually "inaccurate" whether in indoor experiment or in-situ test. Therefore, the displacement back analysis method is usually used to identify the geomechanics of rock mass in engineering. In order to meet the needs of construction, artificial intelligence technology is widely used in displacement back analysis to identify the geomechanical parameters of rock mass. Aiming at the shortcomings of existing artificial neural network (ANN) and support vector regression (SVR) in the application of three-dimensional displacement back analysis, Gaussian process regression (GPR) algorithm is introduced to make up for the shortcomings of existing intelligent identification methods. In order to improve the generalization performance of single kernel GPR algorithm, two single kernel functions are added to form a combined kernel function. An improved Gaussian process regression (IGPR) algorithm based on combined kernel function is proposed. In addition, in the training process of IGPR model, in order to quickly find the super parameters of IGPR model with the best training effect, particle swarm optimization (PSO) and IGPR model are combined to automatically search the parameters of IGPR model to form PSO-IGPR coupling model. FLAC3D software is used for three-dimensional numerical simulation of tunnel engineering construction to form a training sample set with rock mass geomechanical parameters and the distance between tunnel face and monitoring section as input parameters, and the calculated crown subsidence and horizontal convergence as output parameters, which is used for the training of IGPR model. The training steps of the IGPR model are as follows: (1) Initialization of PSO parameters. The particle population size, number of iterations, initial particle velocities and initial locations of particles are determined. Each particle vector represents an IGPR model. Different models correspond to different parameters of the IGPR algorithm. (2) The IGPR algorithm learns knowledge from the learning samples and predicts the displacement value of each testing sample. Then, the PSO algorithm calculates the individual fitness value f i of each particle. (3) Compare the fitness f i calculated in step (2) with the best fitness f(pbest i ) of the particle in the iteration history. If f(pbest i ) > f i , f i is used to replace f(pbest i ) from the previous round, and new particles replace the particles from the previous round. (4) Compare the best fitness f(pbest i ) of each particle with the best fitness f(gbest) of all the particles. If f(pbest i ) < f(gbest), then replace the original global best fitness value f(gbest) with the best fitness value f(pbest i ) of the particle and save the current state of the particle. (5) When the calculation reaches the preset number of evolution generations of PSO, the calculation is finished, and the particle with the minimum fitness value is returned to obtain the optimal solution. Otherwise, a new round of iterations is needed to update the positions and speeds of the particles and generate new particles. Then, the process returns to step (2) until the maximum number of iterations is reached and the calculation is finished.
After the training process is completed, the IGPR model has been able to accurately map the relationship between geomechanical parameters and rock mass deformation, then directly use PSO algorithm to search the best geomechanical parameters of rock mass, so as to make the deformation calculated by IGPR model closest to the measured deformation of rock mass, and complete the identification of rock mass geomechanical parameters. The steps of training process are as follows: 1. The PSO is initialized, and then, the initial population of the geomechanical parameters X of the rock mass is randomly generated with N p individuals. Each individual in the population determines a set of geomechanical parameters X of the rock mass; the generation counter g is set to 1. 2. The best IGPR model reads in each rock mass geomechanics parameter individual of the initial population and the specified distance d between the tunnel face and the monitoring section. The key point displacements can be predicted by the excellent generalization performance of the IGPR algorithm. 3. The PSO reads in the corresponding measured displacements of each key point, calculates the fitness values of each rock mass geomechanical parameter according to the fitness function of PSO, and evaluates the individual fitness of each rock mass geomechanical parameter. 4. Is the prespecified generation number of the PSO reached? If reached, then the calculation stops and the geomechanical parameter individual with the minimal fitness is returned; the optimal geomechanical parameters of the rock mass have been obtained. If not, the calculation procedure proceeds to the next step. 5. Perform similar operations similar to step 3 to 5 of the above training process until the prespecified number of evolution generations of the PSO is reached and the optimal particle individual (G best ) is returned, which is the geomechanical parameter identification value of rock mass. www.nature.com/scientificreports/ The application example of Beikou tunnel shows that the combined kernel function IGPR has higher identification accuracy than the single kernel function GPR and SVR model, the IGPR model with automatic correlation deterministic (ARD) kernel function has higher identification accuracy than the IGPR model with isotropic (ISO) kernel function, and the PSO-IGPR hybrid model based on ARD kernel function has the highest identification accuracy. Therefore, a displacement back analysis method of PSO-IGPR hybrid algorithm based on ARD kernel function is proposed in this paper, which can be used to identify the geomechanical parameters of rock mass and solve other engineering problems.