GAN inversion method of an initial in situ stress field based on the lateral stress coefficient

Abstract

The initial in situ stress field influences underground engineering design and construction. Since the limited measured data, it is necessary to obtain an optimized stress field. Although the present stress field can be obtained by valley evolution simulation, the accuracy of the ancient stress field has a remarkable influence. This paper proposed a method using the generative adversarial network (GAN) to obtain optimized lateral stress coefficients of the ancient stress field. A numerical model with flat ancient terrain surfaces is established. Utilizing the nonlinear relationship between measured stress components and present burial depth, lateral stress coefficients of ancient times are estimated to obtain the approximate ancient stress field. Uniform designed numerical tests are carried out to simulate the valley evolution by excavation. Coordinates, present burial depth, present lateral stress coefficients and ancient regression factors of lateral stress coefficients are input to GAN as real samples for training, and optimized ancient regression factors can be predicted. The present stress field is obtained by excavating strata layers. Numerical results show the magnitude and distribution law of the present stress field match well with measured points, thus the proposed method for the stress field inversion is effective.

Introduction

As a foundation for underground engineering design and stability analysis¹, the initial in situ stress field of a rock mass refers to the natural stress state of the rock mass before engineering construction², which has a great effect on the deformation and failure of the rock mass. To ensure the construction safety of engineering, it is necessary to obtain an optimized initial in situ stress field for deep-buried underground engineering³. The formation of an in situ stress field is a complex process that is affected by topography, geology, rock behavior, tectonic evolution and other factors^4,5,6. Therefore, it is nearly impossible to obtain an optimized in situ stress field for every space point in underground engineering, but it is still important to obtain a more optimized initial stress field by various inversion methods.

In situ stress can be obtained directly by field in situ stress measurement methods, such as strain relief, acoustic emission, borehole cores, and hydraulic fracturing tests^7,8,9,10,11. However, due to limited funding and testing technology, as well as complex geological tectonic and steep terrain conditions, it is hard to obtain a large number of in situ stress measurement results^12,13. Moreover, the measured stress can only reveal the distribution characteristics of the stress field in a local area. Therefore, it is necessary to use an efficient back analysis method combined with advanced numerical simulation to derive a reasonable and applicable initial stress field based on limited measurement results and engineering geological data to meet the needs of engineering^{14,15,16,17,18,19,20}.

• (1) In the multiple linear regression method²¹, a numerical model based on the distribution of rock layers and topographic features is established. By considering gravity and several tectonic stress fields, a multiple linear regression equation is constructed. The solution of the equation gives regression factors for stress fields. The initial stress field is a linear combination of these stress fields multiplied by the corresponding regression factors^22,23,24. This method is convenient and efficient, and the solution is unique, suitable for rock masses with simple geological conditions and has been widely used²⁵. However, for engineering projects with complex geological conditions, in situ stress values obtained by the superposition of regression factors often differ greatly from the measured values. In addition, this method linearly superimposes different in situ stress fields, but in deep-buried underground engineering, the relationship between the in situ stress and burial depth is not linear.

• (2) The boundary load adjustment method^26,27 repeatedly adjusts the calculated boundary load of the numerical model so that the calculated stress values under certain working conditions are close to the measured values. The initial in situ stress obtained by this method can be highly optimized for a specific space point but can hardly fit all measured points. In addition, there is no unique solution for this method, and the trial solution process has no rules to follow and can be time-consuming.

• (3) The displacement back analysis method²⁸ uses the assumed in situ stress field and rock mass mechanical parameters to calculate the displacement of the engineering excavation through numerical simulation. The mechanical parameters of the rock mass are adjusted until the calculated displacement is consistent with the measured displacement. Then, the secondary stress field is obtained^29,30. However, this method is only applicable to underground engineering during construction and cannot be applied in the engineering design stage³¹.

• (4) The valley evolution method^32,33 uses layered excavation of ancient strata to simulate the formation process of the present topography and uses the lateral stress coefficient to indicate the variation in the stress field, thus obtaining the present in situ stress field. However, since large underground engineering stresses are affected by tectonic stress, it is difficult to determine the lateral stress coefficient, and because the lateral stress coefficient is not a constant value at various depths, this method can result in a large error in the initial in situ stress field inversion.

• (5) Artificial intelligence combined with the above methods^34,35,36 can improve the efficiency and accuracy of in situ stress inversion. It can simulate the physiological mechanism of the human brain through a large number of neurons to realize intelligent information processing³⁷. Artificial intelligence has remarkable advantages in nonlinear analysis and fuzzy recognition^38,39,40. Commonly used artificial intelligence methods include neural networks³⁴, genetic algorithms⁴¹, and gray theories². These artificial intelligence methods usually aim at inversion of stress at measured points, not the global stress distribution, which results in good agreement at the measured points but uncertain accuracy elsewhere.

In this paper, a 3D numerical model of the Shuangjiangkou underground hydropower station in ancient times is established, which is used to invert initial in situ stress field. According to measured point data, the lateral stress coefficients of stress components in ancient times are estimated, and multiple combinations of stress components with different lateral stress coefficients are established by a uniform design test. Generative adversarial network (GAN) can automatically establish high-dimensional and global probability distribution of variables by confront training, and can generate new data samples according to the distribution of real data, which can improve the accuracy of rock mechanical parameters and models. Based on the nonlinear analysis of the GAN, two regression factors governing the distribution law of optimized lateral stress coefficients of each stress component in ancient times were determined. After GAN prediction of 12 regression factors for six stress components, optimized lateral stress coefficients are determined. By built model of ancient planation surface with no undulations, the ancient stress field can be reconstructed using the optimized lateral stress coefficients. Then, strata are excavated by layers to simulate valley evolution, and the final stress result represents the present in situ stress field and distribution law. Through comparison and verification with measurement points, this inversion method of initial in situ stress is found to be feasible and suitable for deep buried large underground engineering.

GAN inversion method for initial in situ stress field

The stress field in deep-cut valley area evolves from an ancient regional stress field, accompanied by river valley erosion, surface denudation and other geological effects to shape the valley slope in a long-term unloading process, prompting the rock mass to constantly adjust stress, strain and energy to form a new local stress field⁴². Therefore, top-down excavation of the ancient strata is used to simulate the gradual formation of the river valley geomorphic characteristics and the release of stress in the inversion of the present in situ stress field. The ancient strata denudation process is shown in Fig. 1, and the following assumptions are made about ancient surface and geological formations^43,44.

• (1) The ancient surface is a planation surface with no obvious undulations.

• (2) The ancient stress field consists mainly of a self-gravity stress field and a tectonic stress field, and tectonic movement was completed in the ancient period.

• (3) The present in situ stress field evolved from an ancient stress field and was gradually formed by geological effects such as surface denudation and river erosion.

Figure 1

Schematic diagram of ancient strata denudation process.

Inversion of the ancient stress field based on the lateral stress coefficient

In recent years, a method using lateral stress coefficients to inverse in situ stress fields has been proposed and developed^32,33,42 and is expressed as

$$\sigma_{i}^{{}} (x,y,z) = k_{i} \gamma H$$

(1)

where σ₁–σ₆ are stress components σ_x, σ_y, σ_z, τ_yz, τ_zx and τ_xy. k_i is the lateral stress coefficient for the six stress components. H is burial depth. γ is bulk density.

However, the present stress field is affected by topographical fluctuations and rock mass mechanical properties, leading to the in situ stress coefficient varying greatly in different places³¹. In ancient times, the topography and geomorphology had small fluctuations, and the geological structure was simple. Therefore, utilizing lateral stress coefficients to invert the ancient stress field is more suitable.

Many scholars^25,37 have considered that lateral stress coefficients are a certain value at each place in a rock mass, which means that stress components increase linearly as burial depth increases. However, for deep buried engineering, the lateral stress coefficient is found not to be linear with burial depth⁴⁵. Based on global measured in situ stress data, Brown and Hoek⁴⁵ found that the relationship between the lateral stress coefficient and the burial depth is

$$\frac{100}{H} + 0.3 \le \frac{{\sigma_{H} + \sigma_{h} }}{{2\sigma_{v} }} \le \frac{1500}{H} + 0.5$$

(2)

where σ_H, σ_h and σ_v are the maximum principal stress, minimum principal stress on the horizontal plane and vertical stress, respectively.

Based on the statistical analysis of lateral stress coefficients of different rock masses in mainland China, Jing Feng⁴⁶ reached conclusions similar to the ones of Eq. (2). Therefore, the relationship between the lateral stress coefficients and burial depth is

$$k_{i} = \frac{{a_{i} }}{H} + b_{i}$$

(3)

where a_i and b_i are regression factors.

Based on Eq. (3) and measured stresses in engineering, approximate lateral stress coefficients k_i of different stress components can be obtained. Introducing k_i to Eq. (1), an approximate ancient stress field can be obtained.

Inversion of the present in situ stress field based on GAN

Based on the approximate ancient stress field, numerical simulation can be used to excavate the ancient strata layer by layer to obtain the present stress field. Since the ancient stress field is approximately estimated, there may be a large error in obtaining the present stress field. To solve this problem, this paper introduces the generative adversarial network (GAN) to fine-tune lateral stress coefficients in ancient times to optimize the present in situ stress field to meet the demand of engineering analysis.

GAN is a generative model proposed by Goodfellow et al.⁴⁷, and further developed by Mirza et al.⁴⁸, Odena et al.⁴⁹ and Reed et al.⁵⁰. Its basic thought is inspired by a minimax two-player game, and is widely used in artificial intelligence models, which show high learning and association capabilities and distributed processing capabilities to solve nonlinear problems. GAN consists of a generator and a discriminator trained on confrontational learning mechanism, as shown in Fig. 2. The purpose of GAN is to estimate the potential distribution of existing data and generate new data samples from the same distribution by a generator G which transforms the random variables z into the data that can be faked by continuously learning the probability distribution of real data. The discriminator D is a binary classifier that discriminates whether the input data is real or generated data samples. The two networks are enhanced simultaneously by competing with each other during training, so that the two networks constitute a dynamic game process, until Nash equilibrium⁵¹ is achieved. Both the generator and the discriminator can be designed in conjunction with the current deep neural networks⁵². The calculated process of GAN can be summarized as a binary minimax game, and the objective function can be defined as:

$$\mathop {\min }\limits_{G} \mathop {\max }\limits_{D} V(D,G) = E_{{x\sim P_{data(x)} }} [\lg D(x)] + E_{{z\sim P_{g(z)} }} [\lg (1 - D(G(z)))]$$

(4)

where V(D,G) is the cross-entropy loss of two classifications, P_data(x) is the real data distribution, P_g(z) is the random variables distribution, G(z) is samples generated by the generator based on the random variables, and E(⋅) means the calculated expected value. The global optimal solution is reached when P_data = P_g. For GAN, the loss function of generator and discriminator are expressed as log(D(G(z))) and log(D(x)) + log(1 − D(G(z))).

Figure 2

Calculation process of GAN.

In this paper, deep learning is carried out on the distribution of the lateral stress coefficient based on the superiority of GAN data enhancement. The basic process of using GAN to analyze the lateral stress coefficients is as follows.

• (1) From field survey data, the obtained stress values of the measured points are statistically analyzed, and the effective measured points are selected to determine the approximate lateral stress coefficients by Eq. (3).

• (2) Based on a uniform design test, multiple ancient stress fields with different values of a_i and b_i can be designed for a given floating range of a_i and b_i in Eq. (3).

• (3) According to multiple ancient stress fields, strata excavation is carried out using FLAC3D to obtain multiple present initial stress fields, from which lateral stress coefficients under the present burial depth are determined.

• (4) The coordinates of the measurement point, present burial depth, present lateral stress coefficient and regression factors in ancient times are regarded as real data sample. After data normalization, the data sample are input into the GAN for training.

• (5) After training in GAN, regression factors a_i and b_i of the optimized ancient lateral stress coefficients can be obtained. Then, the optimized ancient stress field is constructed by Eq. (1).

• (6) After excavating the optimized ancient stress field, the optimized present initial stress field can be obtained.

For the selection of a_i and b_i values, many scholars^{25,53,54,55,56} use orthogonal experiments. There are six pairs of regression factors a_i and b_i for six stress components. According to the principle of orthogonal experiments, (6 × 2)² = 144 combinations need to be designed for numerical calculation, which can lead to an excessive number of orthogonal tests and is rather time-consuming. Therefore, it is more appropriate to use a uniform design test instead of an orthogonal test^22,57. The advantage of the uniform design can greatly reduce the number of tests and keep training samples uniformly dispersed. In the uniform design test, multiple combinations with different values of a_i and b_i can be designed. The parameters of the measured points under different combinations of a_i and b_i are

$${\mathbf{x}} = \left[ {x_{1} ,x_{2} , \ldots x_{n} } \right]^{T}$$

(5)

$${\mathbf{y}} = \left[ {y_{1} ,y_{2} , \ldots y_{1n} } \right]^{T}$$

(6)

$${\mathbf{H}}^{{\mathbf{p}}} = \left[ {H_{1}^{p} ,H_{2}^{p} , \ldots H_{n}^{p} } \right]^{T}$$

(7)

$${\mathbf{k}}_{{\mathbf{j}}} = \left[ {\begin{array}{*{20}c} {k_{11} ,k_{12} ,k_{13} ,k_{14} ,k_{15} ,k_{16} } \\ {k_{21} ,k_{22} ,k_{23} ,k_{24} ,k_{25} ,k_{26} } \\ \vdots \\ {k_{{{\text{n}}1}} ,k_{{{\text{n}}2}} ,k_{{{\text{n}}3}} ,k_{{{\text{n}}4}} ,k_{{{\text{n}}5}} ,k_{{{\text{n}}6}} } \\ \end{array} } \right]$$

(8)

$${\mathbf{a}}_{{\mathbf{j}}} = \left[ {\begin{array}{*{20}c} {a_{11} ,a_{12} ,a_{13} ,a_{14} ,a_{15} ,a_{16} } \\ {a_{21} ,a_{22} ,a_{23} ,a_{24} ,a_{25} ,a_{26} } \\ \vdots \\ {a_{n1} ,a_{n2} ,a_{n3} ,a_{n4} ,a_{n5} ,a_{n6} } \\ \end{array} } \right]$$

(9)

$${\mathbf{b}}_{{\mathbf{j}}} = \left[ {\begin{array}{*{20}c} {b_{11} ,b_{12} ,b_{13} ,b_{14} ,b_{15} ,b_{16} } \\ {b_{21} ,b_{22} ,b_{23} ,b_{24} ,b_{25} ,b_{26} } \\ \vdots \\ {b_{n1} ,b_{n2} ,b_{n3} ,b_{n4} ,b_{n5} ,b_{n6} } \\ \end{array} } \right]$$

(10)

where x and y are the coordinate vectors in the horizontal plane of the measured points. n is the number of measured points. H^p is the burial depth vector of the measured points at present. k_j is the lateral stress coefficient matrix of the measured points at present and j is the number of uniform design tests. a_j and b_j are regression factors matrices of stress components σ_x, σ_y, σ_z, τ_yz, τ_zx and τ_xy of measured points in ancient times.

After each uniform design test, the training of the GAN can be carried out. The real data sample input to GAN can be expressed as:

$${\mathbf{P}} = \left[ {\begin{array}{*{20}c} {{\mathbf{x}},{\mathbf{y}},{\mathbf{H}}^{{\mathbf{p}}} ,{\mathbf{k}}_{{\mathbf{1}}} ,{\mathbf{a}}_{{\mathbf{1}}} ,{\mathbf{b}}_{{\mathbf{1}}} } \\ {{\mathbf{x}},{\mathbf{y}},{\mathbf{H}}^{{\mathbf{p}}} ,{\mathbf{k}}_{{\mathbf{2}}} ,{\mathbf{a}}_{{\mathbf{2}}} ,{\mathbf{b}}_{{\mathbf{2}}} } \\ \vdots \\ {{\mathbf{x}},{\mathbf{y}},{\mathbf{H}}^{{\mathbf{p}}} ,{\mathbf{k}}_{{\mathbf{j}}} ,{\mathbf{a}}_{{\mathbf{j}}} ,{\mathbf{b}}_{{\mathbf{j}}} } \\ \end{array} } \right]$$

(11)

After training, input the data of measured points to GAN, the optimized a_i and b_i values can be obtained. The optimized ancient stress field can be calculated by Eq. (3), and the present stress field can then be obtained after a nonlinear elastic–plastic excavation simulation.

Equivalent mechanical parameters of faults

The parameters between the fault and the rock mass are quite different. If the influence of the fault is not considered in the inversion of the present stress field, it causes a large error in the result⁵⁸. Therefore, the parameters of faults should be assigned to the corresponding element before simulating the excavation of ancient strata. Faults and underground caverns do not usually intersect in an orthogonal form and are often simplified into thin elements in numerical simulation analysis. This method of treating faults may bring difficulties in the mesh generation of models due to the limited thickness and complicated intersection relationship of faults.

In this paper, the rock mass is divided into hexahedron elements. Due to the existence of faults, some rock mass elements are cut by faults. These elements contain both rock masses and faults, forming composite elements, as shown in Fig. 3a. For one element, only a set of mechanical parameters can be assigned, so the mechanical parameters of composite elements with both rock masses and faults need to be equalized⁵⁹. It is assumed that the composite element is a transversely isotropic material in the local coordinate system (x′y′z′) established on the fault plane. The composite element is simplified into an equivalent element distributed in layers, and the parameters are shown in Fig. 3b. H_k is the layer thickness, which can be expressed as

$$H_{k} = V_{k} /A$$

(12)

where subscript k equals 1 or 2, representing the parameters of the rock mass and fault, respectively. V_k is volume of rock mass or fault. A is the area of the contact surface between the rock mass and the fault.

Figure 3

Schematic diagram of fault cutting rock mass element.

For the z′ direction, the deformation and stress of the equivalent element are

$$\left\{ {\begin{array}{*{20}l} {\overline{\sigma }_{v} = \sigma_{v1} = \sigma_{v2} } \hfill \\ {\Delta l = \frac{{\sigma_{v1} }}{{E_{1} }}H_{1} + \frac{{\sigma_{v2} }}{{E_{2} }}H_{2} = \frac{{\sigma_{v} }}{{\overline{E}_{v} }}(H_{1} + H_{2} )} \hfill \\ \end{array} } \right.$$

(13)

where ${\tilde{\sigma }}_{v}$, σ_v1 and σ_v2 are the stress applied in the equivalent element, rock mass and fault, respectively, in the z′ direction. ${\overline{E} }_{v}$ is the equivalent elastic modulus in z′ direction.

Based on Eq. (13), the equivalent elastic modulus in the z′ direction is

$$\overline{E}_{v} = {{(H_{1} + H_{2} )} \mathord{\left/ {\vphantom {{(H_{1} + H_{2} )} {\left( {\frac{{H_{1} }}{{E_{1} }} + \frac{{H_{2} }}{{E_{2} }}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {\frac{{H_{1} }}{{E_{1} }} + \frac{{H_{2} }}{{E_{2} }}} \right)}}$$

(14)

For the x′ and y′ directions, assuming that the elongation of the rock mass and fault in the equivalent element are equal, then

$$\left\{ {\begin{array}{*{20}l} {\overline{\varepsilon }_{h} = \varepsilon_{h1} = \varepsilon_{h2} } \hfill \\ {\overline{\sigma }_{h} (H_{1} + H_{2} ) = \sigma_{h1} H_{1} + \sigma_{h2} H_{2} } \hfill \\ \end{array} } \right.$$

(15)

where ${\overline{\varepsilon }}_{h}$, ε_h1 and ε_h2 are the strain in the equivalent element, rock mass and fault in the x′ and y′ directions, respectively. ${\overline{\sigma }}_{h}$, σ_h1 and σ_h2 are the stresses applied in the equivalent element, rock mass and fault in the x′ and y′ directions, respectively.

Equation (15) can be written as

$$\overline{E}_{h} \overline{\varepsilon }_{h} (H_{1} + H_{2} ) = E_{1} \varepsilon_{h1} H_{1} + E_{2} \varepsilon_{h2} H_{2}$$

(16)

where ${\overline{E} }_{h}$ is the equivalent elastic modulus in the x′ and y′ directions.

Based on Eqs. (15) and (16), the equivalent elastic modulus is

$$\overline{E}_{h} = \frac{{E_{1} H_{1} + E_{2} H_{2} }}{{H_{1} + H_{2} }}$$

(17)

Combining Eqs. (15) and (17), the equivalent Poisson’s ratio $\overline{\mu }$ can be expressed as

$$\begin{aligned} \overline{\mu } & = \frac{{\overline{\varepsilon }_{h} }}{{\overline{\varepsilon }_{v} }} = \frac{{\overline{\sigma }_{h} }}{{\overline{E}_{h} }} \cdot \frac{1}{{\overline{\varepsilon }_{v} }} = \frac{{\sigma_{h1} H_{1} + \sigma_{h2} H_{2} }}{{\overline{E}_{h} (H_{1} + H_{2} )}} \cdot \frac{1}{{\overline{\varepsilon }_{v} }} \\ & = \frac{{\mu_{1} E_{1} \varepsilon_{v1} H_{1} + \mu_{2} E_{2} \varepsilon_{v2} H_{2} }}{{\overline{E}_{h} (H_{1} + H_{2} )}} \cdot \frac{1}{{\overline{\varepsilon }_{v} }} \\ & = \frac{{\mu_{1} \sigma_{v1} H_{1} + \mu_{2} \sigma_{v2} H_{2} }}{{\overline{E}_{h} (H_{1} + H_{2} )}} \cdot \frac{{\overline{E}_{v} }}{{\overline{\sigma }_{v} }} \\ & = \frac{{\mu_{1} H_{1} + \mu_{2} H_{2} }}{{H_{1} + H_{2} }} \cdot \frac{{\overline{E}_{v} }}{{\overline{E}_{h} }} \\ \end{aligned}$$

(18)

Project overview and in situ stress measurement analysis

The Shuangjiangkou hydropower station is located on the Dadu River in western Sichuan Province, China. The diversion and power generation structures of the hydropower station are arranged on the left bank of the mountain. The length of the 3D numerical model is 878.14 m in the X direction and 950 m in the Y direction, as shown in Fig. 4. The X direction is N170° E along the axis of the main powerhouse. The length of the main powerhouse in the underground powerhouse is 217.5 m, the arch span is 28.3 m, the maximum excavation depth is 68.3 m, the horizontal burial depth is 400–640 m, and the vertical burial depth is 320–500 m. Two faults F1 and F2 and a lamprophyre L1 are found in the engineering site, as shown in Fig. 5.

Figure 4

3D model of the mountain where the Shuangjiangkou hydropower station located.

Figure 5

Schematic diagram of the underground powerhouses, faults, and measured points.

The distribution of the ten measured points is shown in Fig. 5, and the measured in situ stress data are shown in Table 1. As the measured in situ stress is obtained in the geodetic coordinate system, it needs to be converted to a model coordinate system, as shown in Table 2. The measured in situ stress in the underground powerhouse has the following distribution characteristics.

• (1) The maximum principal stress σ₁ is 14.88–28.96 MPa, the intermediate principal stress σ₂ is 8.53–20.37 MPa, and the minimum principal stress σ₃ is 3.14–10.88 MPa, indicating that the rock mass is in a triaxially unequal pressure state.

• (2) The average values of σ_x/σ_z and σ_y/σ_z are 1.44 and 1.22, respectively, indicating that the measured values of in situ stress show σ_x > σ_y > σ_z. However, σ_y/σ_z of point 5 is significantly less than average.

• (3) The average density ρ and Poisson's ratio μ of the rock mass are 2600 kg/m³ and 0.3, respectively, and the lateral stress coefficient λ = μ/(1 − μ) is 0.4286. The average values of σ_x/λσ_z and σ_y/λσ_z are 3.3559 and 2.8433, respectively, which are greater than 1.0. This indicates that the horizontal stress component in the measured stress is larger than the horizontal stress formed by Poisson’s ratio effect of the rock mass gravity. Therefore, there is a horizontal tectonic stress field in this area.

• (4) σ_x/σ_y is greater than 1.00, indicating that the horizontal tectonic stress in the Y direction is smaller than that in the X direction.

• (5) σ_z/ρgh is the ratio of σ_z at a point to the vertical stress formed by gravity of the overlying rock mass. Except for at the measured point 3, σ_z/ρgh is 1.08–3.06, and the average value is 1.97.

Table 1 In situ stress parameters measured.

Table 2 In situ stress components (MPa).

Point 3 and point 5 are relatively near Fault 2 and the lamprophyre, respectively, which may cause undulation of the measured value¹⁰. Therefore, in the inversion of the initial in situ stress, points 3 and 5 are not used.

Inversion and analysis of stress field

Approximate stress field

By the analysis of the measured in situ stress in “Project overview and in situ stress measurement analysis”, except for measured points 3 and 5, other measured in situ stress points are considered to be representative sample points. In Table 3, the actual lateral stress coefficients k_i for the six stress components can be calculated. Based on Eq. (3), the relationship between the present burial depth and k_i is shown in Fig. 6. The fitted stress function can be imported into FLAC3D to form an approximate ancient stress field. As shown in Fig. 7, the Z direction range of the ancient model is from the bottom elevation of 1970 m to the planation plane. The upper six layers in the model are elements for excavation, which are used to simulate the denudation process of the valley. The parameters input to the ancient model in FLAC3D are listed in Table 4. The coordinates of elements containing both rock masses and faults are exported by the FISH language. Then, the parameters of equivalent elements are calculated based on “Equivalent mechanical parameters of faults” and assigned to the model. The equivalent element is calculated using a transversely isotropic elastic constitutive model.

Table 3 Lateral stress coefficient k_i for six stress components.

Figure 6

Relationship between burial depth and k_i.

Figure 7

The ancient model.

Table 4 The parameters for numerical simulations.

In FLAC3D, each element zone stores six stress components, and each element grid point stores unbalanced forces. If unbalanced forces in the X, Y, Z directions are not applied to model element grid points, stresses in element zones change, and the calculated results of the model greatly deviate from the real results. To optimize the model in subsequent calculations, based on guidelines of the FLAC3D user manual, the model is performed as follows⁶⁰ to perform ground stress balance.

• (1) Based on Eq. (3), ancient stress field can be calculated and imported to the ancient model in FLAC3D, then the parameters of the rock mass are input to the model.

• (2) The boundary normal velocities of the ancient model periphery and bottom surfaces are fixed to zero.

• (3) Calculate the model for one step, that is, run the step 1 command to solve the model, and then record the unbalanced forces on element grid points.

• (4) Apply forces opposite to unbalanced forces recorded in (2) to each element grid point in the ancient model.

• (5) Repeat (3) and (4) several times to make the stress components close enough to the stress components in the ancient model.

After the above process, the balance calculation of the approximate ancient field is completed. Excavating ancient strata layer by layer can obtain the present in situ stress field. Six stress components and lateral stress coefficients at the measured points can be obtained. However, these lateral stress coefficients are obtained by approximating the ancient stress field, so there may be some unacceptable errors, and further optimized regression factors need to be determined.

Optimized stress field

By using a GAN to inversely optimize the present stress field, it is necessary to construct adequate real samples for training. Since the initial in situ stress inversion of this project is to provide accurate ground stress for the underground cavern excavation, the measured points in and around the underground cavern need to be selected as samples. In addition, boundary effects occur in numerical excavation simulations, so measured points near the boundary of the numerical model need to be selected as samples. Therefore, the measured point in Table 3 can be as samples to GAN training. Based on a large number of trial calculations²⁵, the floating range of regression factors in Eq. (3) is determined. The value range of each regression factor a_i or b_i is between 0.5 and 1.5 times of itself. Twelve regression factors are taken as factors of the uniform design test, and different values of regression factors are regarded as different levels. Based on the principle of uniform design, there are 13 levels for regression factors. Then, levels 1–13 correspond to regression factors 0.52, 0.60, 0.68, 0.76, 0.84, 0.92, 1.00, 1.08, 1.16, 1.24, 1.32, 1.40, and 1.48 times, respectively, so the U₁₃ uniform table is selected⁵⁶, as shown in Table 5. Regression factors of each test can be taken into Eqs. (1) and (3) to calculate the stress components of each element, which can be imported into FLAC3D for calculation. In each test, the calculated stress components and lateral stress coefficients of the measured points can be obtained. There are eight effective measured points, so eight training samples can be created in each test, for a total of 13 × 8 = 104 samples.

Table 5 U₁₃ uniform table of levels combination of different factors in each test.

After training with real data samples, input x and y coordinates, present burial depth and lateral stress coefficient of the measured points in the present stress field of measured points, the regression factors of the ancient stress field can be obtained. Then, the optimized ancient stress field can be obtained by Eq. (3). Figure 8 shows the relationship between the loss function of generator and discriminator in GAN and iterations. The generator G and the discriminator network D are continuously confronted, so that the loss function of the two network reduces to a minimum, indicating that the output sample data from generator is stable. After GAN training, the regression factors of the optimized ancient stress field can be predicted, so the optimized ancient stress field can be obtained by Eq. (3). After excavating the optimized ancient stress field, the present stress field can be obtained. The regression factors and calculated stress values of the measured points are listed in Tables 6 and 7. The relative error Δ of the measured points can be calculated:

$$\Delta = \frac{{\left\| {\hat{\sigma }_{i} - \sigma_{i} } \right\|_{2} }}{{\left\| {\sigma_{i} } \right\|_{2} }} \times 100\%$$

(19)

where $\hat{\sigma }_{i}$ is the calculated stress component and $\sigma_{i}$ is the measured stress component. $\left\| \cdots \right\|_{2}$ is the 2-Norm.

Figure 8

The loss function of GAN.

Table 6 Regression factors in ancient stress field obtained by GAN and BP neutral network.

Table 7 Stress component values of each measured point based on GAN and BP neutral network.

For comparison, regression factor of stress components obtained by GAN and BP neural network are also listed in Tables 6 and 7. For these method, the relative errors of σ_x, σ_y and σ_z are less than 15%, which is more optimized than τ_yz, τ_xz and τ_xy. According to statistical data, the error of in situ stress measurement results can reach 25–30%⁶¹. The maximum relative error of the measured points with GAN is 17.46%, so the distribution law of the in situ stress field is reasonable. For BP method, the maximum relative error for Point 8 reaches 28.4%, so the average relative error of GAN is much smaller than that of BP neural network, especially for the inversion of τ_yz, τ_xz and τ_xy.

Analysis of present in situ stress field based on GAN inversion method

Based on the calculated present stress field, the contour maps of y′ = 0 in Fig. 3 including some measured points are shown in Fig. 9. The topography has a significant influence on the in situ stress distribution in the engineering area, while the lithology has a small influence on the stress field. From the top to the bottom of the Earth’s surface, σ_x, σ_y and σ_z increase with increasing depth. Due to the weak mechanical parameters at faults and lamprophyre, there are stress release effects in these areas. The stress components of the faults and lamprophyre all drop sharply, and the stress components of nearby rock masses also decrease. Figure 10 shows the principal stress contour map of Unit 1 in the main powerhouse, indicating that.

• (1) The principal stress changes with the topography, and the surface unloading area has a small deformation modulus due to the relaxation of the rock mass structure, so there is a stress relaxation.

• (2) When the horizontal and vertical burial depths increase, the values of σ₁ and σ₃ gradually increase, and the contour lines are roughly parallel to the hillside.

• (3) With the increase in the depth of the caves, the stress value of the cave surrounding rock also increases, and the stress of the tailwater surge chamber is higher than that of the main powerhouse. The excavation and unloading stresses at the lower part of each cave are higher than those at the top arch, indicating that the excavation unloading effect of the lower part of the cave is more pronounced.

• (4) The stress field presents distribution characteristics of σ_x > σ_y > σ_z, which is advantageous to excavation of the cave sidewalls, but the stability of the end wall is unfavorable.

Figure 9

Stress contour of σ_x, σ_y, σ_z at y′ = 0.

Figure 10

Principal stress contour at x = 0.

The lateral stress coefficients k₁–k₆ along the perpendicular line through Points A, B and C in Fig. 5 of the center of the three caverns are shown in Fig. 11. The calculated values of the lateral stress coefficients are close to the measured points but disperse from the main curve when the location is at the fault or lamprophyre. The lateral stress coefficients decrease with increasing burial depth and gradually tend to a constant value. The fitting curve shows the relationship between burial depth and lateral pressure coefficient Points A, B and C. The relationship is in accordance with Eq. (3), which can be directly imported into the numerical model with a cavern group for excavation simulation.

Figure 11

Relationship between k and burial depth.

Conclusion

The in situ stress field of the Shuangjiangkou Hydropower Station is mainly composed of tectonic stress and gravity stress, and its distribution is affected by the tectonic stress, geological structure and topography. In this paper, the ancient stratum excavation method combined with the lateral stress coefficient and GAN is utilized to invert the in situ stress field. By comparing with measured points, the feasibility of this method to invert the in situ stress field is verified.

• (1) The in situ stress value is positively correlated with the burial depth and the lateral stress coefficient. For deep buried engineering, the lateral stress coefficient is inversely proportional to the burial depth and tends to a constant value, rather than following the commonly used linear relationship.

• (2) Instead of inverting the stress at measured points with the artificial intelligence method, the regression factor of the lateral stress coefficient in ancient times is determined by GAN. Then, the in situ stress field is obtained by strata excavation to simulate the geological effects of surface erosion and weathering to inverse the present in situ stress field.

• (3) Elements containing both rock masses and faults are regarded as equivalent elements of transverse isotropy, thereby reducing the difficulty of numerical model modeling and calculation time.

• (4) The inverted in situ stress components are close to the measured values, so the obtained in situ stress field is reasonable; thus, the inverted in situ stress field can reflect the distribution law of the actual in situ stress field and can provide reasonable guidance for the excavation simulation and stability analysis of underground caverns.

• (5) Based on the in situ stress field obtained by the inversion analysis, the equation of the lateral stress coefficient and the burial depth can be obtained. In subsequent cavern excavation simulations, a refined model with a cavern group can be established, and the reversed geostress equation can be directly imported into the model.

In summary, this paper proposes a new inversion method of the in situ stress field, which can provide a more optimized in situ stress field for deep-buried engineering.