Back Analysis of Geotechnical Engineering Based on Data-Driven Model and Grey Wolf Optimization

Abstract

Geomaterial mechanical parameters are critical to implementing construction design and evaluating stability through feedback analysis in geotechnical engineering. The back analysis is widely utilized to identify and calibrate the geomaterial mechanical properties in geotechnical engineering. This study developed a novel back-analysis framework by combining a reduced-order model (ROM), grey wolf optimization (GWO), and numerical technology. The ROM was adopted to evaluate the response of the geotechnical structure based on a numerical model. GWO was used to search and identify the geomaterials properties based on the ROM. The developed back analysis framework was applied to a circular tunnel and a practical tunnel for determining the mechanical property of the surrounding rock mass. The results showed that the ROM could be an excellent surrogated model and replaced it with the numerical model. The obtained geomaterial properties were in excellent agreement with the actual properties. The deformation behavior captured by the developed framework was consistent with the theoretical solution in a circular rock tunnel. The developed framework provides a practical, accurate, and convenient approach for calibrating the geomaterial properties based on field monitoring data in practical geotechnical engineering applications.

1. Introduction

The geomaterial property is critical to guiding construction and assuring safety during the feedback analysis and dynamic design for geotechnical engineering [1,2,3]. The back analysis based on field monitoring data has been widely used in geotechnical engineering. The back analysis aims to obtain the geomaterial properties by minimizing the objective function, which represents the response difference of geotechnical structure between the monitored field and the predicted by the physical model, based on the optimization technology. The core of the back analysis technique is the optimization technique and physical model. Generally, the numerical method is widely used to calculate the structure response field in practical geotechnical engineering [4,5]. However, a numerical method is time-consuming due to the repetitive computing in back analysis, especially for the practical large-scale geotechnical structure. Meanwhile, optimal techniques often have problems in back analysis, such as high dimensionality, complexity, and local minima. In order to solve the above two problems, various surrogate models and optimal technology have been proposed in the back analysis [2,3,6]. In this study, the reduced order model (ROM) and the grey wolf algorithm (GWO) were combined to improve the practicability of back analysis.

Last century, displacement back analysis was proposed by Sakurai and Takeuchi to identify the rock mechanical parameters based on the monitoring information of the tunnel [7]. However, the numerical method is time-consuming and limits the engineering application of back analysis. In order to overcome the limitations of the numerical method, various machine learning-based surrogate models were a focus of attention to approximate the response of the geotechnical structure in the past decades [8,9,10,11]. The neural network method was utilized to construct the intelligent displacement back analysis model for identifying the mechanical property of the surrounding rock mass [12,13,14,15]. The support vector machine and the relevance vector machine were selected to build a displacement back analysis model to recognize the geomaterials parameters [16,17,18]. Machine learning provides an excellent tool for predicting the structural response and is selected as the surrogate model in the geotechnical back analysis. However, some limits of artificial intelligence, such as overfitting, trapping local minimum, etc., hinder the practical application in back analysis for geotechnical engineering. Meanwhile, the traditional surrogate model does not reflect the physical mechanism of geological engineering and only obtains information about the discrete response at some monitoring points. Therefore, obtaining the total response fields using this classical approach is challenging due to the learning efficiency. The ROM was developed that contains some knowledge about the engineering structure under consideration and can overcome the limits of machine learning and the traditional surrogate model. So, this study adopted the ROM to establish the surrogate model for back analysis.

In order to acquire the appropriate geomaterial property, various intelligent optimization methods, such as genetic algorithm [16,19,20], artificial bee colony [18], particle swarm optimization [21,22], etc., have been widely utilized to seek the geomechanical property in the back analysis. Due to the complexity of geomaterials, trapping the local minimum solution is the main drawback of the back analysis method, which hinders its application in practical engineering. Grey wolf optimization (GWO) is an efficient metaheuristic method developed recently [23]. The GWO algorithm is suitable for solving nonlinear and complex problems due to its simple concept, small number of adjustable parameters, fast convergence, and strong global optimization capacity. It only considers the function evaluation and does not need the derivative information, which is suitable for black box global optimization problems. So, GWO was selected as an optimal technology for the back analysis in this study.

This study proposed a novel back analysis approach by combining numerical models, ROM, and GWO to identify the geomaterial mechanical property in geotechnical engineering. ROM was utilized to construct a surrogate model to approximate and capture the response of the geotechnical structure for replacing the numerical model in the back analysis. The GWO algorithm was regarded as an optimal technology for seeking the unknown geomaterial property based on the idea of the back analysis. The developed framework was applied to a circular tunnel and an actual tunnel project. The remainder of this study is stated as follows. First, ROM and GWO algorithms are briefly introduced in Section 2 and Section 3, respectively. Section 4 introduces the main ideas and procedures of the developed back analysis in detail. In Section 5, a circular tunnel and an actual tunnel, i.e., the experimental tunnel in the Goupitan Water conservancy project, China, are used to verify and investigate the developed back analysis framework. Finally, some conclusions are drawn in Section 6.

2. Reduced-Order Model

2.1. Constructing the ROM Model

The ROM was used to predict the system response using a low-order model based on numerical methods and the proper orthogonal decomposition in engineering. For any x_i, i = 1, 2, …, I, and θ_j, j = 1, 2, …, J, the proper orthogonal decomposition was utilized to obtain the following equation [24].

$${\tilde{u}}^h\left(x_i,{\theta }_j\right)={{\displaystyle \sum }}_{k=1}^K{\beta }_k\left({\theta }_j\right){\phi }^k\left(x_i\right)+\tilde{g}\left(x_i,{\theta }_j\right)$$

(1)

where ${\tilde{u}}^h$ denots the solution of field variables for geotechnical structure, θ_j and x_i denote the parameter and design variables of a geotechnical model, $\phi$ and $\beta$ note the unknown coefficient of ROM, and $\tilde{g}\left(x,\theta \right)$ is an extended boundary condition in the entire domain.

$$\tilde{g}\left(x,\theta \right)=\{\begin{array}{l}g\left(x,\theta \right) on \partial \Omega \\ 0 elsewhere\end{array}$$

(2)

Equation (1) can be rewritten as follows:

$${\tilde{u}}^h=\phi \beta +\tilde{g}$$

(3)

By utilizing the Latin hypercube sampling (LHS), the set of design variables ${\theta }_j, j=1,2,\dots ,J$. was constructed for determining the unknown coefficient $\phi$. Then, the corresponding discrete solutions (snapshots) of the numerical model $w_j=u^h\left({\theta }_j\right)-\tilde{g}\left({\theta }_j\right), j=1,2,\dots ,J,$ were acquired based on numerical techniques such as finite element, discrete element, boundary element, etc. The spatial Gram matrix by M^x can be obtained as follows:

$$M_{ij}^x=\left(w_i·w_j\right), i,j=1,2,\dots ,J$$

(4)

where $\left(w_i·w_j\right)$ notes the inner product between w_i and w_j.

The descending order of the positive eigenvalues of M^x is listed in the following form.

$${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_J\ge 0$$

(5)

The first K eigenfunctions ${\phi }^k\left(x\right), k=1,2,\dots K$ corresponds to the first K eigenvalues, providing the orthogonal principal direction of snapshots. If $r^k={\left(r_j^k\right)}_{j=i,i,\dots J}$ denotes the kth eigenvector of M^x, then its dual kth eigenfunctions ${\phi }^k\left(x\right)$ can be determined according to the following form.

$${\phi }^k\left(x\right)={{\displaystyle \sum }}_{j=1}^K{r}_j^k{w}_j\left(x\right)$$

(6)

where K notes the basis size of the proper orthogonal decomposition and can be solved in the following inequation.

$$\frac{{{\displaystyle \sum }}_{i=1}^K{\lambda }_i^{}}{{{\displaystyle \sum }}_{i=1}^J{\lambda }_i^{}}>k$$

(7)

where k is the user-specified tolerance and equals 0.9999 in this study.

The following penalized minimization problem can solve the unknown coefficient β.

$$\underset{{\beta }_j\in R^K}{\min }\Vert u^{h,j}-\phi {\beta }_j-{\tilde{g}}_j{\Vert }^2+\mu \Vert {\beta }_j{\Vert }^2$$

(8)

${\beta }_j$ can be solved by the following normal equation.

$$\left({\phi }^T\phi +\mu I_K\right){\beta }_j={\phi }^T\left(u^{h,j}-{\tilde{g}}_j\right), j=1,2,\dots ,J$$

(9)

where $\mu$ notes a small regularization parameter.

2.2. Predicting the Field Variables

In order to determine the field variables for unknown variables θ and x, the radial basis function (RBF) was adopted to expand the coefficient ${\beta }_k\left(\theta \right)$ in the following form.

$${\beta }_k\left(\theta \right)={{\displaystyle \sum }}_{j=1}^J{\alpha }_{jk}\psi \left(\frac{\left|\theta -{\theta }_j\right|}{\sigma }\right)$$

(10)

For any ${\theta }_{j^{\prime }}, j^{\prime }=1,2,\dots ,J$, Equation (10) can be presented as follows.

$${{\displaystyle \sum }}_{j=1}^J{\alpha }_{jk}\psi \left(\frac{\left|{\theta }_{j^{\prime }}-{\theta }_j\right|}{\sigma }\right)={\beta }_{k{j}^{\prime }}$$

(11)

where ${\beta }_{k{j}^{\prime }}$ are determined by Equation (9). The above equations can be rewritten in the compacted form:

$$A{\alpha }_k={\beta }_k$$

(12)

The following equation can solve the unknown coefficient ${\alpha }_k$.

$$\left(A^{T}A+\mu I_J\right){\alpha }_k=A^T{\beta }_k, k=1,2,\dots ,K$$

(13)

2.3. Procedure of the ROM Model

This study developed a ROM by combining the numerical method (including finite element, discrete element and boundary element, etc.) and the proper orthogonal decomposition. Using the LHS, the set of design variables was constructed for the snapshot. Then the numerical method was adopted to calculate the corresponding solution of design variable in the above set. Based on the obtained snapshots, the proper orthogonal decomposition basis vector and its coefficient were acquired using the proper orthogonal decomposition algorithm for geotechnical engineering. In order to acquire the unknown field of the new design variable, RBF functions were adopted to expand the coefficient of the proper orthogonal decomposition basis and then to determine the coefficients of the orthogonal decomposition ROM. The unknown field variables corresponding to the new design were determined using the ROM. Figure 1 shows the main flowchart and procedure of the ROM. In what follows, the procedure of the ROM model is presented in detail.

• Step 1: Collect the data of the geotechnical engineering, including project property, geo-stress, boundary conditions, etc.;
• Step 2: Establish the numerical model (FEM) based on the above engineering information;
• Step 3: Construct the design variables set θ for the numerical model using LHS;
• Step 4: Calculate the field variables w_i (displacement or stress field) at space domain X using a numerical method for each design variable. Collect all the field variables and acquire the snapshots;
• Step 5: Build the spatial Gram matrix M^x based on the above snapshots;
• Step 6: Solve the eigenvalues λ and eigenvectors r based on the spatial Gram matrix;
• Step 7: Determine the rank number K of M^x and the first K eigenfunction vector φ;
• Step 8: Determine the undetermined coefficient β based on eigenfunction vector φ and snapshots;
• Step 9: To a new design variable θ, construct element ɸ based on the design variables θ generated by LHS using the RBF function;
• Step 10: Determine the interpolation matrix A of elements ɸ;
• Step 11: Determine the vector of element α using the penalized linear systems;
• Step 12: Solve the coefficients β(θ) based on the RBF function;
• Step 13: Calculate the unknown field variables ${\tilde{u}}^h\left(\theta \right)$ based on coefficients β(θ) and eigenfunction vector φ using the ROM.

3. Grey Wolf Optimization (GWO)

Grey Wolf Optimization (GWO) is a heuristic optimization strategy inspired by the social hierarchy and hunting techniques of grey wolves. GWO mimics the leadership hierarchy and hunting technology of grey wolves. The hierarchy of the grey wolves is divided into four levels (alpha (α), beta (β), delta (δ), and omega (ω)), which present the optimal solution, the suboptimal solution, the third optimal solution, and the remaining candidate solutions, respectively [23]. There are three main stages of grey Wolf hunting, namely searching, encircling, and attacking prey. The GWO is the mathematical model of the hunting strategy and social hierarchy of grey wolf. A grey wolf can determine the prey by randomly changing its position based on the GWO algorithm. The hunting process is guided by α, β, δ, and ω according to the above three kinds of wolves in GWO.

In the encircling phase, the encircling behavior of wolves can be expressed in the following mathematical model:

$$\overrightarrow{D}=\left|\overrightarrow{D}\cdot {\overrightarrow{X}}_p(t)-\overrightarrow{X}(t)\right|$$

(14)

$$\overrightarrow{X}(t+1)={\overrightarrow{X}}_p(t)-\overrightarrow{A}\cdot \overrightarrow{D}$$

(15)

where t denotes the step of the iteration, $\overrightarrow{D}$ denotes the searching vector, $\overrightarrow{X}$ and ${\overrightarrow{X}}_p$ are the vector and denote the position of a grey wolf and the prey, respectively. $\overrightarrow{A}$ and $\overrightarrow{C}$ denote the coefficient vectors and can be determined according to the following form:

$$\overrightarrow{A}=2\overrightarrow{a}\cdot {\overrightarrow{r}}_1-\overrightarrow{a}$$

(16)

$$\overrightarrow{C}=2{\overrightarrow{r}}_2$$

(17)

where ${\overrightarrow{r}}_1$ and ${\overrightarrow{r}}_2$ denote the vector and selected randomly in the range of zero to unity, the component $\overrightarrow{a}$ decreases linearly from 2 to 0 with the iterations.

In the hunting phase, the locations of other search agents (including omega) were updated according to alpha, beta, and delta knowledge based on the following equations.

$${\overrightarrow{D}}_{\alpha }=\left|{\overrightarrow{C}}_1\cdot {\overrightarrow{X}}_{\alpha }-\overrightarrow{X}\right|$$

(18)

$${\overrightarrow{D}}_{\beta }=\left|{\overrightarrow{C}}_2\cdot {\overrightarrow{X}}_{\beta }-\overrightarrow{X}\right|$$

(19)

$${\overrightarrow{D}}_{\delta }=\left|{\overrightarrow{C}}_3\cdot {\overrightarrow{X}}_{\delta }-\overrightarrow{X}\right|$$

(20)

$${\overrightarrow{X}}_1={\overrightarrow{X}}_{\alpha }-{\overrightarrow{A}}_1\cdot {\overrightarrow{D}}_{\alpha }$$

(21)

$${\overrightarrow{X}}_2={\overrightarrow{X}}_{\beta }-{\overrightarrow{A}}_2\cdot {\overrightarrow{D}}_{\beta }$$

(22)

$${\overrightarrow{X}}_3={\overrightarrow{X}}_{\delta }-{\overrightarrow{A}}_3\cdot {\overrightarrow{D}}_{\delta }$$

(23)

$$\overrightarrow{X}(t+1)=\frac{{\overrightarrow{X}}_1+{\overrightarrow{X}}_2+{\overrightarrow{X}}_3}{3}$$

(24)

where the subscripts α, β, and δ represent the alpha wolf, beta wolf, and delta wolf, respectively.

In the attacking prey phase, the final attack is determined by decreasing the $\overrightarrow{a}$ from 2 to 0 with the iterations $\overrightarrow{A}$ selected randomly in the range [$-2\overrightarrow{a}$, $2\overrightarrow{a}$]. $\overrightarrow{A}$ will decrease with the reduction to the $\overrightarrow{a}$ and force the wolfs to approach the prey while $\left|\overrightarrow{A}\right|$ is less than 1.

In search of prey, grey wolves follow the leader, dispersing from one another in search of prey and gathering to attack. The number of wolves N_w and the generation N_G are the two essential parameters of GWO. Each generation represents the decision movement of a wolf. The number of wolves represents the function computational times in each generation. Figure 2 shows the flowchart of GWO and the main procedure. A detailed introduction of the GWO can be found in the literature [23].

4. ROM-Based Back Analysis Using GWO

Back analysis technology has been commonly utilized to identify the geomaterial property in geotechnical engineering. The numerical method and the optimal technology are the two critical elements of back analysis. This study developed a ROM-based back analysis combining numerical methods, ROM, and GWO. The ROM model was utilized to predict the nonlinear response of the geotechnical structure based on the numerical method. GWO was selected as an optimal technology to seek the geomaterial properties.

4.1. Back Analysis

In the 1970s, back analysis was proposed to identify the rock mass properties in rock engineering [25]. The back analysis provides a simple but effective way to identify the geomaterial properties based on the field data and numerical analysis. It also provides a helpful tool for guiding the dynamic design, reinforcement of surrounding rock mass, and safe construction of geotechnical engineering. Figure 3 shows the main parts of the back analysis and its basic idea. Field measurements provide basic information for back analysis. The physical model is the heart of the back analysis. Due to the complex geological conditions, it is not easy to determine the closed-form solution of geotechnical engineering. Meanwhile, optimal technology is essential to back analysis due to multi-extremum and multi-constrained optimization problems. This study adopted the ROM model to capture the physical model. GWO was selected as the optimal technology due to its excellent global optimizing capability.

4.2. ROM-Based Surrogate Model

The physical model characterizes the nonlinear implicit function mapping between geomaterial properties and their response during construction. In this study, a ROM-based physical model was established to capture the nonlinear mapping function between the geomaterial properties (Deformation modulus, Poisson’s ratio, strength property, and in-situ stress) and corresponding structural response (displacement, stress, strain, etc.). The following equations define the physical model ROM(X).

$$ROM\left(\mathit{X}\right):R^N\to R^Q,$$

(25)

Y = ROM (X),

(26)

where X = (x₁, x₂, …, x_N) is a vector and x_i (i = 1, 2, …, N) is the ith geomaterial properties and Y = (y₁, y₂, …, y_N) is a Q dimension vector and denotes the response induced by construction.

Some known training samples are necessary for establishing the surrogate model ROM(X) of the physical model in the back analysis. It is necessary to obtain the training samples for ROM based on a numerical method and the design of experiment.

4.3. Objection Function

This study constructed the objective function based on the geotechnical structure response difference between the field value and ROM prediction. The objection function forced the optimal technology to seek the optimal variable. The following root means square defines the objective function:

$$fitness=\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(y_{pi}-y_i\right)}^2/n}$$

(27)

where n denotes the number of monitoring points, y_i and y_pi denote the predicted by ROM and monitoring response of the geotechnical structure in ith measurement point.

4.4. Procedure of the Developed Framework

This study developed a novel back-analysis framework combining ROM, GWO, and numerical technology. The ROM was adopted to capture the nonlinear mapping between the geomaterial properties and the corresponding response during excavation in combination with the ROM and numerical model. The experimental design was adopted to construct the combination of the unknown properties, and then the numerical method was utilized to calculate the structural response at each combination. The snapshots consist of a combination of the unknown parameters and the corresponding response. ROM was built based on the above snapshots. GWO is an optimal approach to seeking geomaterial properties based on the ROM. The detailed procedures of the proposed method are listed as follows (as seen in Figure 4):

• Step 1: Collect the engineering data, such as the unknown (need to determine by back analysis) and known geomaterial mechanical and physical properties, boundary conditions, and the range of unknown geomaterial properties;
• Step 2: Generate the combination of the unknown properties based on experimental design and calculate the structural response at each training sample. The snapshots consist of the combination of the unknown parameters and the corresponding response;
• Step 3: Based on the determined snapshots, generate the ROM to capture the nonlinear function mapping between the geomaterial properties and the corresponding structural response in geotechnical engineering;
• Step 4: Establish the objective function and call the GWO to seek the geomaterial properties based on the monitored data during the construction.

5. Numerical Example and Application

5.1. Numerical Example

A circular tunnel is excavated in a continuous, homogeneous, and isotropic rock mass. The hydrostatic far-field stress p₀ and uniform support pressure p_i are shown in Figure 5. When support pressure p_i is not enough to meet critical pressure p_cr, a plastic zone will exist. The values of p_cr could be computed as follows:

$$p_{cr}=\frac{2p_0-{\sigma }_c}{k+1}$$

(28)

where σ_c notes the uniaxial compression strength. It could be obtained using the following equation.

$${\sigma }_c=\frac{c(k-1)}{\tan \phi }$$

(29)

where φ and c denote the cohesion and the friction angle, respectively. k could be determined as follows:

$$k=\frac{1+\sin \phi }{1-\sin \phi }$$

(30)

According to the Mohr–Coulomb criterion, Duncan (1993) analyzed and inferred the inward displacement of tunnel wall u_ip and the plastic zone radius r_p [26].

$$\frac{r_p}{r_0}={\left[\frac{2(p_0+s)}{(k+1)(p_i+s)}\right]}^{\frac{1}{k-1}}$$

(31)

$$\frac{u_{ip}}{r_0}=\left[\frac{(1+\mu )}{E}\right]\left[2(1-\mu )(p_0-p_{cr}){(\frac{r_p}{r_0})}^2-2(1-2\mu )(p_0-p_i)\right]$$

(32)

where E notes the elastic modulus and μ notes Poisson’s ratio. The values are computed as follows:

$$s=\frac{{\sigma }_c}{k-1}$$

(33)

In this study, back analysis was utilized to identify the far-field stress p_0, cohesion c, and friction angle φ based on the deformation of the surrounding rock mass in the tunnel. Five horizontal direction monitoring points were placed to record the deformation of surrounding rock mass in the circular tunnel. The distance between the center of the tunnel and the five monitored points are 1.0 m, 1.2 m, 1.6 m, 1.8 m, and 2.0 m, respectively. The displacements of monitored points could be calculated using the above formula (Equation (32)). The tunnel radius is 1.0 m. The value of far-field stress p_0, cohesion c, and friction angle φ are 32 MPa, 6.8 GPa, 3.2 MPa, and 32°, respectively. The displacements of 5 monitored points were calculated by the analytical solution and adopted as field measurement to back-calculate the unknown parameters of the surrounding rock mass using the proposed method. The snapshots were constructed and generated based on the experimental design and a numerical method.

Once the snapshots are obtained, the surrogate model could be established according to the ROM algorithm. Figure 6 shows the calculated displacement comparison between the ROM and the analytical solution (Equation (32)). The predicted displacement using the ROM surrogate model is in good agreement with the analytical solution. It shows that the ROM captured well the nonlinear function mapping between unknown properties and the tunnel deformation. The ROM-based surrogate model provides a feasible way to replace the analytical solution in back analysis.

ROM-based back analysis was utilized to identify the unknown properties of the surrounding rock mass using the above ROM surrogate model and to predict displacement in the tunnel. The far-field stress p₀, cohesion c, and friction angle φ are 32.01 MPa, 6.73 GPa, 3.40 MPa, and 31.03°, respectively (

Table 1. The results and comparison.

	Actual	This Study	Relative Error (%)
p0/MPa	32.00	32.01	−0.03
E/MPa	6800.00	6730.87	1.02
c/MPa	3.20	3.40	−6.25
φ/°	32.00	31.03	3.03

). The relative error is −0.03%, 1.02%, −6.25% and 3.03%, respectively. The maximum relative error is less than 7%. It shows that the identified parameters agree with the actual parameters using the proposed back analysis framework. The displacements comparison between the predicted by back analysis, ROM surrogate model, and calculated by the analytical by actual parameters are shown in Figure 7. Figure 8 shows the stress and displacement of the surrounding rock mass in the tunnel calculated based on the different methods. The results show that the proposed framework can be utilized to identify the mechanical property of the surrounding rock mass in the tunnel.

Figure 9 shows the convergence of the unknown mechanical property of the surrounding rock mass using the developed back analysis method. The convergence property of the developed back analysis is shown in Figure 10. The unknown property can converge to the final value quickly. The developed method has excellent convergence and global optimization performance.

5.2. Application: Goupitan Experimental Tunnel

Goupitan Hydropower Station is a landmark west-to-east transmission project located on the Wujiang River in Guizhou Province, China [27]. According to the preliminary design of the underground powerhouse, the right bank tailwater tunnel and the construction diversion tunnel pass through the soft clay rock mass. An experiment tunnel with a buried depth of 70 m was excavated to understand the rheological characteristics of clay rock. The tunnel size was 2 m in width and 2 m in height, respectively. Some monitoring points are set up to obtain the deformation of surrounding rock during tunnel excavation. Rock stratum $S_{2h}^{1-2}$ and $S_{2h}^{1-1}$ are located approximately 3 m below and 30 m above the tunnel, respectively (Figure 11). The 4# and 6# borehole of 7 m depth is set at the position of 11.6 m in the tunnel, where the 4# borehole is horizontal, the 5# borehole is 45° inclined, and the 6# borehole is vertical. Five monitoring points numbered 1–5 were arranged at a depth of 0, 1, 2, 4, and 6 m (Figure 11).

Table 2. The monitoring displacement in the monitored borehole.

Time (Day)	Displacement (mm)
	4#	6#
3	2.558	1.778
5	3.789	2.377
11	4.531	2.685

lists the deformation of the 3-day, 5-day, and 11-day at each monitoring point.

The developed back analysis framework was utilized to identify the rheological properties of the rock mass based on the 3rd, 5th, and 11th day monitored displacements at 4# and 6# boreholes. For the developed back analysis, a physical model is critical to identify the rheological mechanical properties of the rock mass. In this study, the rheological properties of the rock mass include the shear modulus G₁ and viscosity η₁ for the Kelvin model, the shear modulus G₂ and viscosity η₂ for the Maxwell model. In geotechnical engineering, it is not easy to determine an analytical solution for use. Although the numerical method is commonly used to understand the rheological mechanism and deformations behavior of the rock mass, it is time-consuming in practical and large-scale geotechnical engineering.

In this study, a ROM surrogate model, which replaced the numerical method, was used to improve the efficiency of back analysis. A uniform design was utilized to construct a group of 42 samples, and Fast Lagrangian Analysis of Continua (FLAC) software was utilized to solve the displacement of the tunnel wall in the rock mass. Rock masses $S_{2h}^{1-1}$ and $S_{2h}^{1-2}$ were regarded as Burger‘s material [28,29], and their rheological mechanical properties were identified by the developed method.

Table 3. Ranges of unknown rheological properties.

Clay-Green Clay Rock ${\mathit{S}}_{2\mathit{h}}^{1-1}$				Purple Clay Rock ${\mathit{S}}_{2\mathit{h}}^{1-2}$
G1h (GPa)	G2h (GPa)	η2h (GPa·d)	η1h (GPa·d)	G1z (GPa)	G2z (GPa)	η2z (GPa·d)	η1z (103 GPa·d)
0.5–4.5	0.1–3.5	0.1–3.5	15–35	1–15	5–20	1–15	1.5–4.5

lists the range of unknown rheological mechanical parameters based on laboratory tests and a field survey. Figure 12 indicates the predicted displacement comparison by the ROM and the numerical model along 4# and 6# boreholes, respectively. The calculated deformations by the ROM were found very close to the deformation computed by the numerical solution. It proves that the ROM surrogate model characterized well the nonlinear mapping relationship between unknown rheological properties and deformation of the tunnel wall in the rock mass. Hence, it could replace the numerical method in the back analysis.

According to the back analysis theory and the proposed method, the rheological parameters of $S_{2h}^{1-1}$ and $S_{2h}^{1-2}$ are determined based on the displacement of the surrounding rock mass on the 3rd, 5th and 11th days of the 4# and 6# boreholes (

Table 4. Obtained rheological parameters.

Number of Monitored Day		3rd, 5th and 11th
Clay-green clay rock S2h1−1	G1h (GPa)	1.39
	G2h (GPa)	0.20
	η2h (GPa·d)	0.12
	η1h (GPa·d)	35.00
Purple clay rock S2h1−2	G1z (GPa)	1.00
	G2z (GPa)	20.00
	η2z (GPa·d)	8.68
	η1z (103 GPa·d)	1.96

). The developed method identified the rheological mechanical properties in a reasonable way. Compared with other methods, the rheological mechanical properties obtained by the proposed framework are closer to the actual properties of rock mass, and the rheological properties could be determined dynamically by rationally using the on-site displacement data monitored during construction. The developed back analysis framework costs approximately 73.65 s in PC with Intel(R) Xeon(R) Gold 5218 CPU @ 2.30 GHz &2.29 GHz to obtain the rheological properties. However, it takes approximately 2 min for a single tunnel stability analysis using the numerical simulation method. It is obvious that the ROM-based surrogate model could dramatically improve the efficiency of back analysis.

It is critical to identify and understand the deformation and failure mechanism of the surrounding rock mass. The developed back analysis identifies the rheological properties based on the displacements monitored. The deformation of the rock mass was investigated based on the rheological properties identified by the developed back analysis framework using the ROM and GWO. The ROM then predicted the displacements of the monitored borehole. Their comparisons are shown in Figure 13 in which it is evident that the displacement determined by the ROM is in good agreement with the monitored displacement during the excavation. It also shows that the identified rheological properties well characterized the rheological behavior of the rock mass during the tunnel construction. This confirms that the developed method can be used for determining rheological properties and evaluating the time-dependent behavior of the rock mass. The rheological properties obtained by the developed method can be used for stability analysis, design, and safety construction during excavation in rock engineering. In addition, the predicted deformation by ROM diverges from the monitoring data on the 34th day due to the complexity of the construction site, which brings errors and uncertainty to the monitoring data. With the increasing monitoring data, the rheological properties of the surrounding rock mass will be updated dynamically to capture the trend.

The relationship between the rheological property (clay-green clay rock) and the number of iterations is shown in Figure 14 for the GWO. The variation process objective function is plotted in Figure 15. The results are similar to the above numerical example. GWO could seek the appropriate mechanical property of the rock mass quickly. It proved again that the GWO is an excellent optimal technology and has a good performance of global optimization.

6. Conclusions

This study developed a novel back analysis framework to identify the geomaterial property by combining a numerical model, ROM, and GWO. The ROM was utilized to establish the surrogate model for capturing deformation during excavation in geotechnical engineering. The numerical method was adopted to construct snapshots for the ROM based on the design of the experiment. The geomaterial properties were identified based on the monitored displacement data by the developed framework. Meanwhile, GWO was selected as an optimal technology for back analysis due to its global optimization performance. The developed back analysis framework was illustrated successfully by a numerical example and Goupitan experimental tunnel. The geomaterial properties identified by this study were compared with the ROM prediction. The predicted displacement of the surrounding rock mass was also close to the actual monitored displacement. The results showed that the developed back analysis framework provides a convenient, practical, and accurate way to understand the geomaterial properties based on the monitored response during construction.

(1)

The ROM model was utilized to construct a low-order surrogate model for capturing the response-induced excavation in geotechnical engineering and replacing the numerical model in the back analysis. It is critical to practical engineering due to the difficulties in obtaining the analytical solution for geotechnical engineering;

(2)

Back analysis is a scientific and practical tool widely used in geotechnical engineering. The numerical model and optimal technology are the two critical components of back analysis. The developed back analysis framework takes full advantage of the merits of ROM and GWO and provides a feasible way for determining the property of the surrounding rock mass in geotechnical engineering;

(3)

ROM is an excellent physics-based data-driven surrogate model that can capture the mechanism of surrounding rock mass. GWO is an efficient metaheuristic method developed recently and is suitable for solving the black-box problem. However, ROM depends on the numerical fidelity model, and the parameters of the GWO algorithm influence the optimal performance. In a future study, the authors will further improve the developed framework by absorbing and combining the advantages and merits of various methods.