# Local Stress Field Correction Method Based on a Genetic Algorithm and a BP Neural Network for In Situ Stress Field Inversion

Authors: Tianzhi Yao; Zuguo Mo; Li Qian; Jianhua He; Jianhai Zhang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The in situ stress field is the fundamental factor causing deformation and damage in geotechnical engineering, so it is the main basis for underground engineering design and excavation. However, it is difficult to accurately obtain the in situ stress through most existing inversion methods in areas with complex geological conditions. For the problem of a relatively discrete and nonlinear relationship of measured stress in the Yebatan Hydropower Station area, a new in situ stress inversion method called the local stress field correction (LSFC) method combining a genetic algorithm (GA), backpropagation (BP) neural network, and submodel method is proposed. The inverted in situ stress results produced by this method show that the distribution of in situ stress is greatly influenced by tectonic movements in the Yebatan area, there is no obvious linear relationship with depth, and the stress release phenomenon occurs at the faults. By comparison with the multiple regression method, it is found that the method still has high inversion accuracy under complex geological conditions, and the average relative error of LSFC inversion results is 17.05%, which is much lower than the value of 43.58% via the multiple regression method. Therefore, the LSFC method can be used for the inversion of in situ stress in complex geological regions and provide a reference for engineering design and construction.


---

## Body

## 1. Introduction

In situ stress is the most basic engineering load in geotechnical engineering and is the main basis for underground engineering design, excavation, and stability analysis [1, 2]: this stress exerts a great influence on the failure and deformation of the rock mass [3]. Therefore, it is essential to simulate the distribution of the in situ stress field, which is influenced by a number of factors, including topography, lithology, geologic structure, gravity, and groundwater [4, 5].The Yebatan Hydropower Station is located on the upper reaches of the Jinsha River in the southeastern Qinghai-Tibet Plateau. The area has undulating topography, frequent geological movements, complex geological structure, and many faults. Complex in situ stress fields can cause large deformation and damage in underground caverns during excavation, presenting unpredictable safety risks. Therefore, grasping the distribution law of the in situ stress field is of great significance to the design and construction of engineering projects.In situ measurement is a direct method of obtaining the in situ stress field, and the main methods of measurement include stress relief, acoustic emission, overcoring, hydraulic fracturing test, and AE-based ground stress methods [6–9]. However, due to site and financial constraints, the scope and number of measurement points are limited, especially in complex and steep mountainous regions. In addition, the causes of the in situ stress field are complex, and there are many influencing factors. The actual measurement results often have a certain degree of dispersion and exhibit certain limitations [10–13]. The actual distribution of the in situ stress field cannot be assessed solely on the basis of the measured results. Therefore, to obtain a more reasonable and applicable range of in situ stress fields to meet the needs of the project, in situ stress inversion should be carried out using effective analysis methods based on actual measurement data and geological structure conditions [14–19].At present, the main inversion methods used by scholars are as follows:(1)
Multiple linear regression inversion method [20–23]. This method is based on the employment of topography and lithology to build a numerical model, taking into account the effects of gravity and tectonic movement, and constructing linear regression equations, for which the regression coefficients of each stress field caused by each factor individually can be solved. The actual in situ stress field is a linear combination of these regression coefficients multiplied by the corresponding stress fields [24]. The method is convenient, efficient, and widely used for situations where the distribution of in situ stress is smooth and simple [18]. However, the results obtained by this method differ significantly from the measured stress values under complex geological conditions [25].(2)
Displacement back analysis method [26–28]. First, the method assumes an initial stress field, yields the mechanical parameters of the rock mass, and obtains the excavation disturbance displacement field by numerical calculation. Then, the mechanical parameters of the rock are adjusted until the calculated values correspond to the actual monitoring displacement values. Finally, the inverse in situ stress is obtained. However, this method can only be used during the construction period.(3)
Lateral pressure coefficient method [29, 30]. The lateral pressure coefficients are calculated according to the relationship between the different horizontal stress components and the vertical stresses at the measured points. The initial stress field is obtained by multiplying the lateral pressure coefficient under the action of gravity.(4)
Submodel secondary inversion method [31, 32]. The numerical model is divided into two models: a large model for obtaining the approximate distribution law of the stress field and a submodel for refining the analysis. The region including the critical research subjects is removed from the large-scale model as the submodel. Finally, the stress field in the large model is mapped into the small model to complete the secondary inversion of the submodel and obtain a more refined stress field. The method is suitable for the study of critical areas containing complex geological structures.(5)
Artificial intelligence combining the above methods. Several artificial intelligence algorithms, such as genetic algorithms (GA) [25, 33, 34] and neural networks [35–37], have been developed and are widely used for in situ stress inversion. These methods can find nonlinear relationships between stress values, boundary conditions, and relative positions, so they can be used to invert the in situ stress fields under complex geological conditions.In summary, there are a number of methods that can be used for the inversion of the in situ stress field, but few of them can solve the problem of inversion where measured stress shows no linear relationship with burial depth or topography. The multiple regression method requires a linear relationship between the stresses and the depth of burial and between the stresses at the measured points. The lateral pressure coefficient method requires that the stress tends to increase with increasing depth. Current artificial intelligence methods invert the in situ stress field mainly by changing the structure of the neural network. In addition to changing the structure of the neural network, multiple regression and lateral pressure methods are often implemented through the neural network or directly fit the relationship between location and stress, which suits simple geological conditions [38, 39].The author attempted to apply the multiple regression method and lateral pressure coefficient method to the Yebatan Hydropower Station area, but it was found that the results were unsatisfactory for extremely complicated geological conditions and discrete measured stresses. Therefore, this paper proposes a new method called the local stress field correction (LSFC) method to solve the problem of in situ stress inversion in complex geological conditions such as the Yebatan area. The LSFC method combines GA, a BP neural network, and a submodel method. This method can directly invert the stresses in the Yebatan area based on the relationship between the stresses and the relative position of the measured points without any assumptions about geological conditions. The inverted in situ stress in the Yebatan area is compared with the inversion results of the multiple regression method, which justifies the validity of the proposed LSFC method.


## 2. In Situ Stress Inversion Method Based on Local Stress Field Correction by a GA and BP Neural Network

### 2.1. Approximate Stress Field Obtained by the Genetic Algorithm

The genetic algorithm is a global search optimization algorithm inspired by Darwin’s theory of evolution. GA achieves satisfactory performance even in noisy environments and does not require any linearizing assumption. The process of the GA is as follows: (1) randomly generate an initial population; (2) determine the fitness function; (3) select individuals with high fitness from the current population; (4) generate new individuals through crossover and mutation; (5) terminate the algorithm if the termination conditions are met and the global optimum is obtained; otherwise, repeat the above steps.Different fitness functions are used to solve different problems, so a fitness function that is suitable for revealing the distribution law of the initial stress field must be determined. The rock mass ground stress field is mainly caused by gravity and geological tectonism. The influences of temperature and groundwater can be ignored [40]. Therefore, gravity and tectonic movement have been chosen as the fundamental influencing factors for the inversion, which is simulated by applying the force of gravity and the unit distributed loads to the numerical model boundaries [41]. Based on the idea of fitting the measured stress and the calculated stress values of the numerical model in the multiple linear regression method [42, 43], an inversion method based on GA was proposed. The fitness function in GA can be defined as shown in equation (1):(1)fitness=min1M∑k=1MγkFuk−Fmk2,where M is the number of stress measurement points, k indicates the six stress components (k = 1–6), Fmk is the measured value of the kth stress component, γk is the coefficient of the kth stress field caused by the unit displacement loading on the boundary, Fuk denotes the kth stress components at the measured points caused by loading unit displacement at the boundary, as shown in Figure 1, Fu1 is the stress caused by a 1-meter uniform extrusion displacement Ux along the x-axis applied in the y-z plane, as shown in Figure 1(a), Fu2 is the stress caused by a 1-meter uniform extrusion displacement Uy along the y-axis applied to the x-z plane, as shown in Figure 1(b), Fu3 is the stress caused by gravity G, as shown in Figure 1(c), Fu4 is the stress caused by a 1-meter tangential distribution of displacement Uyz along the z-axis applied to the x-z plane, as shown in Figure 1(d), Fu5 is the stress caused by a 1-meter tangential distribution of displacement Uxz along the z-axis applied to the y-z plane, as shown in Figure 1(e), and Fu6 is the stress caused by a 1-meter tangential distribution of displacement Uxy along the y-axis applied to the y-z plane, as shown in Figure 1(f). The above loading and calculation are implemented by FLAC3D.Figure 1
Applying boundary displacement and gravity. (a)X-direction extrusion tectonic movement. (b) Y-direction extrusion tectonic movement. (c) Z-direction gravity movement. (d) YZ-direction shearing tectonic movement. (e) XZ-direction shearing tectonic movement. (f) XY-direction shearing tectonic movement.
(a)(b)(c)(d)(e)(f)The optimalγk value can be obtained by selection, crossover, and mutation in GA. The stress components of each element in the numerical model can be expressed as equation (2); an approximate stress field may be thus obtained.(2)σxxn=γ1σxxun,σyyn=γ2σyyun,σzzn=γ3σzzun,τyzn=γ4τyzun,τxzn=γ5τxzun,τxyn=γ6τxyun.The results obtained from equation (2) can approximately show the distribution of the stress field, but there is still a difference between the calculated stress and the actual stress values at the measurement points.


### 2.2. Local Stress Field Correction Based on BP Neural Network Optimization in a Large Model

The BP neural network achieves good performance in solving nonlinear problems. It is composed of an input layer, hidden layer, and output layer, with different layers containing multiple parallel neuronal signals. The BP neural network produced by forward propagation is combined with the backward adjustment of the error to find the optimal weights between the nodes and ultimately obtain the prediction model. Forward propagation and backward propagation are the two main steps of the BP neural network, and these two steps proceed together until the results of the model are close to the expected output.


#### 2.2.1. Design of Input Layer Parameters

It is important to determine the appropriate input and output layers in the BP neural network because they affect the performance of the inversion. The BP neural network is applied to establish the relationship between the depth, the distance between measured points and other points in the numerical model, and the stress components:(3)Fσxx,σyy,σzz,τyz,τxz,τxy=dxnk,dynk,hnT,where F (σxx, σyy, σzz, τyz, τxz, τxy) is the calculated value of the six stress components, n is the total number of nodes in the numerical model, k is the number of measured points, hn is the burial depth of point n in the numerical model, and dxnk and dynk are the horizontal distances between point n in the numerical model and a measured point k in the x-direction and in the y-direction, respectively, as shown in Figure 2.Figure 2
Distance between thenth point in the numerical model and measured points in the horizontal x-y plane.To study the distribution law of the local stress, the relative positions of then points must first be determined. The stress at the surface is 0, so the burial depth is used to determine the relative position in the vertical direction. Both the x-y coordinates and dxnk and dynk can reflect the relative positions of points in the x-y plane. However, setting dxnk and dynk as the input parameters is more conducive to the training of the neural network than the x-y coordinates, so dxnk and dynk are used to determine the relative positions of points in the horizontal direction of the numerical model. The input layer matrix is expressed as follows:(4)X=dx11dy11dx12dy12⋯dx1kdy1kh1dx21dy21dy22dy22⋯dx2kdy2kh2dx31dy31dx32dy32⋯dx3kdx3kh3⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯dxn1dyn1dxn1dxn2⋯dxnkdynkhnn×2k+1,where X is the input matrix of the BP neural network, the number of rows of X is determined by the number of data points in the numerical model, and the number of columns is determined by the number of measured points, which can be expressed as (2k + 1), in which the last column is the burial depth hn of each point. The 2k and (2k − 1)th columns are the horizontal distances of the points from measurement point k. The structure of the BP neural network is shown in Figure 3.Figure 3
The structure of the BP neural network.


#### 2.2.2. Local Stress Field Correction Method

Under complex geological conditions, the ground stress calculated by equation (2) is often very coarse and cannot be used for further excavation analysis, so a method called local stress field correction (LSFC) has been proposed to correct the stress according to the measured values. The LSFC method is shown in Figure 4.Figure 4
Local stress field correction method data points. (a) Horizontalx-y plane data points. (b) Vertical x-z plane data points.
(a)(b)In Figure4, the red point is the measured point, the blue plane in Figure 4(a) is a section of the horizontal X–Y plane through the measured point, the yellow plane in Figure 4(b) is a section of the vertical X-Z plane through the measured point, and the other black points indicate the relative position and stress components extracted at that point as training samples for the BP neural network. The black points are selected at regular intervals from the measured points up to the boundary of the model.The local stress field correction process is as follows:Step 1. The stress fields caused by different influencing factors shown in Figure 1 are calculated separately, and an approximate stress field is obtained by equation (2).Step 2. The stress component values calculated by Step 1 are applied to the numerical model built in FLAC3D and the program is run to equilibrate the stress field.Step 3. The stress component values obtained in Step 2 at the measured points are replaced by the measured values.Step 4. Along the X, Y, and Z directions, points are taken every 50 m from the measured point, that is, dxnk = 0, 50, 100, …, dynk = 0, 50, 100, …, h = 0, 50, 100, …, up to the model boundary, as shown in Figure 4. The data for each point consist of position information and stress components.Step 5. A neural network model following the mapping relationship in equation (3) is built using the data in Steps 3 and 4. This neural network model can be trained to predict the stress components after inputting hn, dxnk, and dynk for points in the numerical model. Therefore, we obtain the predicted six stress components for the whole numerical model.Step 6. The corrected stress components are applied to the numerical model, the FLAC3D program is run again to equilibrate the stress field, and the output stress field is taken as the inverse stress field of the large model.It should be pointed out that, after the training of the BP neural network in Step 5, the stress values at the measured points approach the real stress values and are in accordance with the actual stress field distribution law. At Step 6, the stress values at the measured points in the numerical model are corrected after the stress equilibrium calculation.


### 2.3. Secondary Inversion Method for Submodels with Faults

The secondary inversion method selects the key area in the large model as the submodel and further refines the submodel grid for stress field inversion [44]. Its main theoretical basis is Saint-Venant’s principle, which states that when the actual distributed load is replaced by an equivalent load, the stress and strain are only changed around the area where the load is applied. Therefore, as long as the submodel boundary is far from the stress concentration area, a more accurate model can be obtained [45].The implementation process is as follows:Step 1. Map the stress values calculated in Section 2.2 from the large model to the submodel by the interpolation procedure and run the stress equilibrium calculation in FLAC3D.Step 2. The stress component values in Step 1 at the measured points are replaced by the measured values.Step 3. Along the X, Y, and Z directions, points are taken every 50 m from the measured point, that is, dxnk = 0, 50, 100, …, dynk = 0, 50, 100, …, h = 0, 50, 100, …, up to the submodel boundary, as shown in Figure 4. The data for each point consist of position information and stress components.Step 4. A neural network model following the mapping relationship in equation (3) is built using the data in Step 2 and Step 3. This neural network model is trained to predict the stress components after inputting hn, dxnk, and dynk for points in the numerical model. Therefore, we obtain the corrected six stress components for the submodel.Step 5. The corrected stress components are applied to the numerical model, the FLAC3D program is run again to equilibrate the stress field, and the output stress field is taken as the final inverted stress field of the submodel with faults.


## 2.1. Approximate Stress Field Obtained by the Genetic Algorithm

The genetic algorithm is a global search optimization algorithm inspired by Darwin’s theory of evolution. GA achieves satisfactory performance even in noisy environments and does not require any linearizing assumption. The process of the GA is as follows: (1) randomly generate an initial population; (2) determine the fitness function; (3) select individuals with high fitness from the current population; (4) generate new individuals through crossover and mutation; (5) terminate the algorithm if the termination conditions are met and the global optimum is obtained; otherwise, repeat the above steps.Different fitness functions are used to solve different problems, so a fitness function that is suitable for revealing the distribution law of the initial stress field must be determined. The rock mass ground stress field is mainly caused by gravity and geological tectonism. The influences of temperature and groundwater can be ignored [40]. Therefore, gravity and tectonic movement have been chosen as the fundamental influencing factors for the inversion, which is simulated by applying the force of gravity and the unit distributed loads to the numerical model boundaries [41]. Based on the idea of fitting the measured stress and the calculated stress values of the numerical model in the multiple linear regression method [42, 43], an inversion method based on GA was proposed. The fitness function in GA can be defined as shown in equation (1):(1)fitness=min1M∑k=1MγkFuk−Fmk2,where M is the number of stress measurement points, k indicates the six stress components (k = 1–6), Fmk is the measured value of the kth stress component, γk is the coefficient of the kth stress field caused by the unit displacement loading on the boundary, Fuk denotes the kth stress components at the measured points caused by loading unit displacement at the boundary, as shown in Figure 1, Fu1 is the stress caused by a 1-meter uniform extrusion displacement Ux along the x-axis applied in the y-z plane, as shown in Figure 1(a), Fu2 is the stress caused by a 1-meter uniform extrusion displacement Uy along the y-axis applied to the x-z plane, as shown in Figure 1(b), Fu3 is the stress caused by gravity G, as shown in Figure 1(c), Fu4 is the stress caused by a 1-meter tangential distribution of displacement Uyz along the z-axis applied to the x-z plane, as shown in Figure 1(d), Fu5 is the stress caused by a 1-meter tangential distribution of displacement Uxz along the z-axis applied to the y-z plane, as shown in Figure 1(e), and Fu6 is the stress caused by a 1-meter tangential distribution of displacement Uxy along the y-axis applied to the y-z plane, as shown in Figure 1(f). The above loading and calculation are implemented by FLAC3D.Figure 1
Applying boundary displacement and gravity. (a)X-direction extrusion tectonic movement. (b) Y-direction extrusion tectonic movement. (c) Z-direction gravity movement. (d) YZ-direction shearing tectonic movement. (e) XZ-direction shearing tectonic movement. (f) XY-direction shearing tectonic movement.
(a)(b)(c)(d)(e)(f)The optimalγk value can be obtained by selection, crossover, and mutation in GA. The stress components of each element in the numerical model can be expressed as equation (2); an approximate stress field may be thus obtained.(2)σxxn=γ1σxxun,σyyn=γ2σyyun,σzzn=γ3σzzun,τyzn=γ4τyzun,τxzn=γ5τxzun,τxyn=γ6τxyun.The results obtained from equation (2) can approximately show the distribution of the stress field, but there is still a difference between the calculated stress and the actual stress values at the measurement points.


## 2.2. Local Stress Field Correction Based on BP Neural Network Optimization in a Large Model

The BP neural network achieves good performance in solving nonlinear problems. It is composed of an input layer, hidden layer, and output layer, with different layers containing multiple parallel neuronal signals. The BP neural network produced by forward propagation is combined with the backward adjustment of the error to find the optimal weights between the nodes and ultimately obtain the prediction model. Forward propagation and backward propagation are the two main steps of the BP neural network, and these two steps proceed together until the results of the model are close to the expected output.


### 2.2.1. Design of Input Layer Parameters

It is important to determine the appropriate input and output layers in the BP neural network because they affect the performance of the inversion. The BP neural network is applied to establish the relationship between the depth, the distance between measured points and other points in the numerical model, and the stress components:(3)Fσxx,σyy,σzz,τyz,τxz,τxy=dxnk,dynk,hnT,where F (σxx, σyy, σzz, τyz, τxz, τxy) is the calculated value of the six stress components, n is the total number of nodes in the numerical model, k is the number of measured points, hn is the burial depth of point n in the numerical model, and dxnk and dynk are the horizontal distances between point n in the numerical model and a measured point k in the x-direction and in the y-direction, respectively, as shown in Figure 2.Figure 2
Distance between thenth point in the numerical model and measured points in the horizontal x-y plane.To study the distribution law of the local stress, the relative positions of then points must first be determined. The stress at the surface is 0, so the burial depth is used to determine the relative position in the vertical direction. Both the x-y coordinates and dxnk and dynk can reflect the relative positions of points in the x-y plane. However, setting dxnk and dynk as the input parameters is more conducive to the training of the neural network than the x-y coordinates, so dxnk and dynk are used to determine the relative positions of points in the horizontal direction of the numerical model. The input layer matrix is expressed as follows:(4)X=dx11dy11dx12dy12⋯dx1kdy1kh1dx21dy21dy22dy22⋯dx2kdy2kh2dx31dy31dx32dy32⋯dx3kdx3kh3⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯dxn1dyn1dxn1dxn2⋯dxnkdynkhnn×2k+1,where X is the input matrix of the BP neural network, the number of rows of X is determined by the number of data points in the numerical model, and the number of columns is determined by the number of measured points, which can be expressed as (2k + 1), in which the last column is the burial depth hn of each point. The 2k and (2k − 1)th columns are the horizontal distances of the points from measurement point k. The structure of the BP neural network is shown in Figure 3.Figure 3
The structure of the BP neural network.


### 2.2.2. Local Stress Field Correction Method

Under complex geological conditions, the ground stress calculated by equation (2) is often very coarse and cannot be used for further excavation analysis, so a method called local stress field correction (LSFC) has been proposed to correct the stress according to the measured values. The LSFC method is shown in Figure 4.Figure 4
Local stress field correction method data points. (a) Horizontalx-y plane data points. (b) Vertical x-z plane data points.
(a)(b)In Figure4, the red point is the measured point, the blue plane in Figure 4(a) is a section of the horizontal X–Y plane through the measured point, the yellow plane in Figure 4(b) is a section of the vertical X-Z plane through the measured point, and the other black points indicate the relative position and stress components extracted at that point as training samples for the BP neural network. The black points are selected at regular intervals from the measured points up to the boundary of the model.The local stress field correction process is as follows:Step 1. The stress fields caused by different influencing factors shown in Figure 1 are calculated separately, and an approximate stress field is obtained by equation (2).Step 2. The stress component values calculated by Step 1 are applied to the numerical model built in FLAC3D and the program is run to equilibrate the stress field.Step 3. The stress component values obtained in Step 2 at the measured points are replaced by the measured values.Step 4. Along the X, Y, and Z directions, points are taken every 50 m from the measured point, that is, dxnk = 0, 50, 100, …, dynk = 0, 50, 100, …, h = 0, 50, 100, …, up to the model boundary, as shown in Figure 4. The data for each point consist of position information and stress components.Step 5. A neural network model following the mapping relationship in equation (3) is built using the data in Steps 3 and 4. This neural network model can be trained to predict the stress components after inputting hn, dxnk, and dynk for points in the numerical model. Therefore, we obtain the predicted six stress components for the whole numerical model.Step 6. The corrected stress components are applied to the numerical model, the FLAC3D program is run again to equilibrate the stress field, and the output stress field is taken as the inverse stress field of the large model.It should be pointed out that, after the training of the BP neural network in Step 5, the stress values at the measured points approach the real stress values and are in accordance with the actual stress field distribution law. At Step 6, the stress values at the measured points in the numerical model are corrected after the stress equilibrium calculation.


## 2.2.1. Design of Input Layer Parameters

It is important to determine the appropriate input and output layers in the BP neural network because they affect the performance of the inversion. The BP neural network is applied to establish the relationship between the depth, the distance between measured points and other points in the numerical model, and the stress components:(3)Fσxx,σyy,σzz,τyz,τxz,τxy=dxnk,dynk,hnT,where F (σxx, σyy, σzz, τyz, τxz, τxy) is the calculated value of the six stress components, n is the total number of nodes in the numerical model, k is the number of measured points, hn is the burial depth of point n in the numerical model, and dxnk and dynk are the horizontal distances between point n in the numerical model and a measured point k in the x-direction and in the y-direction, respectively, as shown in Figure 2.Figure 2
Distance between thenth point in the numerical model and measured points in the horizontal x-y plane.To study the distribution law of the local stress, the relative positions of then points must first be determined. The stress at the surface is 0, so the burial depth is used to determine the relative position in the vertical direction. Both the x-y coordinates and dxnk and dynk can reflect the relative positions of points in the x-y plane. However, setting dxnk and dynk as the input parameters is more conducive to the training of the neural network than the x-y coordinates, so dxnk and dynk are used to determine the relative positions of points in the horizontal direction of the numerical model. The input layer matrix is expressed as follows:(4)X=dx11dy11dx12dy12⋯dx1kdy1kh1dx21dy21dy22dy22⋯dx2kdy2kh2dx31dy31dx32dy32⋯dx3kdx3kh3⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯dxn1dyn1dxn1dxn2⋯dxnkdynkhnn×2k+1,where X is the input matrix of the BP neural network, the number of rows of X is determined by the number of data points in the numerical model, and the number of columns is determined by the number of measured points, which can be expressed as (2k + 1), in which the last column is the burial depth hn of each point. The 2k and (2k − 1)th columns are the horizontal distances of the points from measurement point k. The structure of the BP neural network is shown in Figure 3.Figure 3
The structure of the BP neural network.


## 2.2.2. Local Stress Field Correction Method

Under complex geological conditions, the ground stress calculated by equation (2) is often very coarse and cannot be used for further excavation analysis, so a method called local stress field correction (LSFC) has been proposed to correct the stress according to the measured values. The LSFC method is shown in Figure 4.Figure 4
Local stress field correction method data points. (a) Horizontalx-y plane data points. (b) Vertical x-z plane data points.
(a)(b)In Figure4, the red point is the measured point, the blue plane in Figure 4(a) is a section of the horizontal X–Y plane through the measured point, the yellow plane in Figure 4(b) is a section of the vertical X-Z plane through the measured point, and the other black points indicate the relative position and stress components extracted at that point as training samples for the BP neural network. The black points are selected at regular intervals from the measured points up to the boundary of the model.The local stress field correction process is as follows:Step 1. The stress fields caused by different influencing factors shown in Figure 1 are calculated separately, and an approximate stress field is obtained by equation (2).Step 2. The stress component values calculated by Step 1 are applied to the numerical model built in FLAC3D and the program is run to equilibrate the stress field.Step 3. The stress component values obtained in Step 2 at the measured points are replaced by the measured values.Step 4. Along the X, Y, and Z directions, points are taken every 50 m from the measured point, that is, dxnk = 0, 50, 100, …, dynk = 0, 50, 100, …, h = 0, 50, 100, …, up to the model boundary, as shown in Figure 4. The data for each point consist of position information and stress components.Step 5. A neural network model following the mapping relationship in equation (3) is built using the data in Steps 3 and 4. This neural network model can be trained to predict the stress components after inputting hn, dxnk, and dynk for points in the numerical model. Therefore, we obtain the predicted six stress components for the whole numerical model.Step 6. The corrected stress components are applied to the numerical model, the FLAC3D program is run again to equilibrate the stress field, and the output stress field is taken as the inverse stress field of the large model.It should be pointed out that, after the training of the BP neural network in Step 5, the stress values at the measured points approach the real stress values and are in accordance with the actual stress field distribution law. At Step 6, the stress values at the measured points in the numerical model are corrected after the stress equilibrium calculation.


## 2.3. Secondary Inversion Method for Submodels with Faults

The secondary inversion method selects the key area in the large model as the submodel and further refines the submodel grid for stress field inversion [44]. Its main theoretical basis is Saint-Venant’s principle, which states that when the actual distributed load is replaced by an equivalent load, the stress and strain are only changed around the area where the load is applied. Therefore, as long as the submodel boundary is far from the stress concentration area, a more accurate model can be obtained [45].The implementation process is as follows:Step 1. Map the stress values calculated in Section 2.2 from the large model to the submodel by the interpolation procedure and run the stress equilibrium calculation in FLAC3D.Step 2. The stress component values in Step 1 at the measured points are replaced by the measured values.Step 3. Along the X, Y, and Z directions, points are taken every 50 m from the measured point, that is, dxnk = 0, 50, 100, …, dynk = 0, 50, 100, …, h = 0, 50, 100, …, up to the submodel boundary, as shown in Figure 4. The data for each point consist of position information and stress components.Step 4. A neural network model following the mapping relationship in equation (3) is built using the data in Step 2 and Step 3. This neural network model is trained to predict the stress components after inputting hn, dxnk, and dynk for points in the numerical model. Therefore, we obtain the corrected six stress components for the submodel.Step 5. The corrected stress components are applied to the numerical model, the FLAC3D program is run again to equilibrate the stress field, and the output stress field is taken as the final inverted stress field of the submodel with faults.


## 3. Engineering Profile and In Situ Stress Measurement Analysis

### 3.1. Engineering Profile

The Yebatan surface elevation of the underground plant area above 2,870 m is steep bedrock with a topographic slope of 40–50°, and the slope locally reaches up to 55–60°; the area below the elevation of 2,870 m to the edge of the Jinsha River is a steep slope zone of overburden with a general slope of 32–40°.The underground plant cavern group of Yebatan Hydropower Station is large in scale, the main powerhouse is 250.9 m long with a span of 25.9 m, and the maximum excavation height for the plant is 71.54 m. The main transformer room is 176.00 m long, with an excavation span of 19.40 m and an excavation height of 26.85 m. The total length of the tailwater pressure chamber is 155.00 m, and the section size of the pressure chamber is 25.0 m × 70.0 m × 87.1 m (width × length × height). The numerical model parameters are shown in Table1.Table 1
Rock mass parameters related to the numerical model.
CategoryModulus of deformationE (GPa)Poisson’s ratioμSpecific weightt (m3)Factor of frictionfCohesionC (MPa)Unweathered18.50.212.721.281.70Weakly weathered9.00.262.701.081.10Strongly weathered2.50.332.670.700.60Fault0.20.362.230.350.05


### 3.2. Regional Geological Formations

The basement rock stratigraphy exposed in the dam site area is a single, complex rock strain composed of moderately acidic rocks that intruded during the Variscan period, mainly reflecting quartz diorite and granodiorite. There are no major regional fractures passing through the site area, and the tectonics are mainly faults, joints, and fissures. Faults F2, F3, f80, f26, f23, f9, f85, and f28 developed in the plant area, as shown in Table2 and Figure 5.Table 2
Fault information.
Fault numberAttitudeWidth (m)F2NE80°/S∠60°1.0F3NW70°/SW∠40°1.4f80NE85°/S∠50°0.3f26EW/S∠55°0.4f23NE73°/NW∠75°0.6f9NW85°S∠60°1.0f85NE85°S∠55°0.3f28NE65°NW∠80°0.6Figure 5
Schematic diagram of model range, faults, and locations of measured points. The red points in the figure are the measured points.


### 3.3. In Situ Stress Measurement Analysis

The measured in situ stress results and the stress components are listed in Tables3 and 4. Figures 6 and 7 show the relationship between stress, stress components, and depth, respectively, for which it can be found that in situ stresses cannot be expressed by a simple linear relationship. The measurement stress in Yebatan has the following distribution characteristics:(1)
The maximum principal stressσ1 is 16.51–26.30 MPa, the intermediate principal stress σ2 is 8.06–19.51 MPa, and the minimum principal stress σ3 is 3.85–8.51. The fitted curves in Figures 6 and 7 show that the stress values tend to increase with burial depth. However, the curve fit performance is poor, and the stress distribution is not only related to the burial depth.(2)
The average values ofσxx/σzz and σyy/σxx are 1.30 and 1.73, respectively, indicating that σyy > σxx > σzz in Yebatan. This indicates that the distribution of the stress field is significantly influenced by tectonic movements.(3)
The average Poisson’s ratioμ of the rock mass is 0.27, and the lateral pressure coefficient λ = μ/(1 − μ) is 0.37. The average values of σxx/λσzz and σyy/λσzz are 3.51 and 4.67, respectively, which are larger than 1.0. This indicates that the horizontal stress component of the actual stress field is greater than the stress components created by the Poisson ratio effect of the rock mass gravity. Therefore, the Yebatan area is subject to the action of both tectonic movement and gravity.Table 3
Measured in situ stress results and stress components.
Measured pointDepth (m)σ1σ2σ3Magnitude (MPa)Direction (°)Dip angle (°)Magnitude (MPa)Direction (°)Dip angleMagnitude (MPa)Direction (°)Dip angle (°)118316.51288.230.012.3625.312.16.32134.757.2216318.51305.712.09.4640.822.84.349.9−63.9329226.30283.0−6.719.51199.845.13.856.544.1431223.30262.128.48.0610.830.85.21138.545.7534724.54278.012.712.18153.368.58.5112.017.1σ1 is the maximum principal stress, σ2 is the intermediate principal stress, and σ3 is the minimum principal stress.Table 4
Coordinate stress components for measured in situ stress (MPa).
Measured pointσxxσyyσzzτyzτxzτxy112.7813.289.14−4.410.44−1.2228.7717.825.72−2.85−1.61−1.07315.2922.3612.015.026.17−6.43414.5112.0010.06−5.694.52−6.06512.8819.8912.47−1.772.32−6.68σxx, σyy, σzz, τyz, τxz, and τxy are the stress components, where positive normal stress indicates tension and negative normal stress indicates compression.Figure 6
Relationship between measured stress and depth.Figure 7
Relationship between measured stress components and depth.In summary, the law of stress field distribution and the geological conditions in the Yebatan are relatively complex.


## 3.1. Engineering Profile

The Yebatan surface elevation of the underground plant area above 2,870 m is steep bedrock with a topographic slope of 40–50°, and the slope locally reaches up to 55–60°; the area below the elevation of 2,870 m to the edge of the Jinsha River is a steep slope zone of overburden with a general slope of 32–40°.The underground plant cavern group of Yebatan Hydropower Station is large in scale, the main powerhouse is 250.9 m long with a span of 25.9 m, and the maximum excavation height for the plant is 71.54 m. The main transformer room is 176.00 m long, with an excavation span of 19.40 m and an excavation height of 26.85 m. The total length of the tailwater pressure chamber is 155.00 m, and the section size of the pressure chamber is 25.0 m × 70.0 m × 87.1 m (width × length × height). The numerical model parameters are shown in Table1.Table 1
Rock mass parameters related to the numerical model.
CategoryModulus of deformationE (GPa)Poisson’s ratioμSpecific weightt (m3)Factor of frictionfCohesionC (MPa)Unweathered18.50.212.721.281.70Weakly weathered9.00.262.701.081.10Strongly weathered2.50.332.670.700.60Fault0.20.362.230.350.05


## 3.2. Regional Geological Formations

The basement rock stratigraphy exposed in the dam site area is a single, complex rock strain composed of moderately acidic rocks that intruded during the Variscan period, mainly reflecting quartz diorite and granodiorite. There are no major regional fractures passing through the site area, and the tectonics are mainly faults, joints, and fissures. Faults F2, F3, f80, f26, f23, f9, f85, and f28 developed in the plant area, as shown in Table2 and Figure 5.Table 2
Fault information.
Fault numberAttitudeWidth (m)F2NE80°/S∠60°1.0F3NW70°/SW∠40°1.4f80NE85°/S∠50°0.3f26EW/S∠55°0.4f23NE73°/NW∠75°0.6f9NW85°S∠60°1.0f85NE85°S∠55°0.3f28NE65°NW∠80°0.6Figure 5
Schematic diagram of model range, faults, and locations of measured points. The red points in the figure are the measured points.


## 3.3. In Situ Stress Measurement Analysis

The measured in situ stress results and the stress components are listed in Tables3 and 4. Figures 6 and 7 show the relationship between stress, stress components, and depth, respectively, for which it can be found that in situ stresses cannot be expressed by a simple linear relationship. The measurement stress in Yebatan has the following distribution characteristics:(1)
The maximum principal stressσ1 is 16.51–26.30 MPa, the intermediate principal stress σ2 is 8.06–19.51 MPa, and the minimum principal stress σ3 is 3.85–8.51. The fitted curves in Figures 6 and 7 show that the stress values tend to increase with burial depth. However, the curve fit performance is poor, and the stress distribution is not only related to the burial depth.(2)
The average values ofσxx/σzz and σyy/σxx are 1.30 and 1.73, respectively, indicating that σyy > σxx > σzz in Yebatan. This indicates that the distribution of the stress field is significantly influenced by tectonic movements.(3)
The average Poisson’s ratioμ of the rock mass is 0.27, and the lateral pressure coefficient λ = μ/(1 − μ) is 0.37. The average values of σxx/λσzz and σyy/λσzz are 3.51 and 4.67, respectively, which are larger than 1.0. This indicates that the horizontal stress component of the actual stress field is greater than the stress components created by the Poisson ratio effect of the rock mass gravity. Therefore, the Yebatan area is subject to the action of both tectonic movement and gravity.Table 3
Measured in situ stress results and stress components.
Measured pointDepth (m)σ1σ2σ3Magnitude (MPa)Direction (°)Dip angle (°)Magnitude (MPa)Direction (°)Dip angleMagnitude (MPa)Direction (°)Dip angle (°)118316.51288.230.012.3625.312.16.32134.757.2216318.51305.712.09.4640.822.84.349.9−63.9329226.30283.0−6.719.51199.845.13.856.544.1431223.30262.128.48.0610.830.85.21138.545.7534724.54278.012.712.18153.368.58.5112.017.1σ1 is the maximum principal stress, σ2 is the intermediate principal stress, and σ3 is the minimum principal stress.Table 4
Coordinate stress components for measured in situ stress (MPa).
Measured pointσxxσyyσzzτyzτxzτxy112.7813.289.14−4.410.44−1.2228.7717.825.72−2.85−1.61−1.07315.2922.3612.015.026.17−6.43414.5112.0010.06−5.694.52−6.06512.8819.8912.47−1.772.32−6.68σxx, σyy, σzz, τyz, τxz, and τxy are the stress components, where positive normal stress indicates tension and negative normal stress indicates compression.Figure 6
Relationship between measured stress and depth.Figure 7
Relationship between measured stress components and depth.In summary, the law of stress field distribution and the geological conditions in the Yebatan are relatively complex.


## 4. Engineering Applications

### 4.1. Large Model Stress Field Inversion Analysis by the LSFC Method

The size of the large model in the horizontal plane is 1100 m × 1100 m, and the coordinate position and orientation are shown in Figure5. The numerical model in FLAC3D is shown in Figure 8. The z-axis is consistent with the geodetic coordinates along the vertical direction.Figure 8
Meshes of the large model.According to Step 1 of Section 2.2.2, first, the unit displacement is loaded to the boundary, as shown in Figure1, to obtain the stress distribution under the basic influencing factors. Using equation (1) as the fitness function, γk of each stress component was obtained by the GA algorithm in MATLAB version 2016a. The calculated results are γ1 = 0.951, γ2 = 1.019, γ3 = 1.214, γ4 = −8.66, γ5 = 14.47, and γ6 = 32.86. The stress components of each element in the numerical model can be obtained from equation (2). Then, through Step 2 of 2.2.2, the stress field of the initial inversion is obtained after the stress equilibrium calculation. The results are shown in Table 5, and the relative error is defined as (5)Δ=Calculated value−Measured value2Measured value2×100%,where 2 is the 2-norm.Table 5
Comparison of measured stress component values calculated by GA and LSFC.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—GA12.4716.477.35−1.902.25−4.3127.09LSFC10.8812.009.05−3.67−0.951.3417.922Measured stress8.7717.825.72−2.85−1.61−1.07—GA12.5715.406.63−1.762.09−3.9731.84LSFC8.9113.606.84−3.870.810.8125.953Measured stress15.2922.3612.015.026.17−6.43—GA12.1215.1911.38−1.902.82−1.7638.13LSFC13.8919.299.47−2.244.70−2.2330.34Measured stress14.5112.0010.06−5.694.52−6.06—GA10.7517.6510.71−2.612.99−2.4036.24LSFC11.6014.007.86−2.153.50−2.2228.925Measured stress12.8819.8912.47−1.772.32−6.68—GA13.2912.8811.00−1.803.25−1.2932.51LSFC13.7016.3710.63−1.701.74−1.9422.57An approximate in situ stress field was calculated by the GA method with a relative error of 27.09% to 38.13%, which can be further optimized by the LSFC method.Through Step 3 of Section2.2.2, the stress component values at the measured points are replaced by the measured values.Through Steps 4 and 5 of Section2.2.2, each of the 5 measured points is used as the center point, and the data are extracted for BP neural network training. Since there are 5 measured points, the column number of equation (4) is 11. There are 11 input data points (hn, dxn1, dyn1, dxn2, dyn2, dxn3, dyn3, dxn4, dyn4, dxn5, and dyn5) and 6 target data points (σxx, σyy, σzz, τyz, τxz, and τxy), and 15 nodes are set in the hidden layer of the BP neural network. After training, a corrected stress component prediction model can be obtained.Finally, through Step 6 of Section2.2.2, the relative position of all elements is converted into the form of equation (4) and imported into the trained BP neural network to calculate the stress components of all elements. After performing the stress equilibrium calculation in FLAC3D, a more accurate stress field was calculated.The results of the LSFC method are shown in Table5. In comparison with the GA method, the LSFC method was more accurate. The x = 135 m section is in the direction of the main powerhouse axis and passes through measured points 1 and 2, as shown in Figure 5. Figures 9(a) and 9(b) show the maximum and minimum principal stress contours for the x = 135 m section, respectively. The maximum principal stresses are more densely distributed and have a large gradient in the river valley. The distribution of the maximum and the minimum principal stresses in the central area of the model away from the boundary increases with increasing burial depth and changes with topography. The gradient of the minimum principal stress changes relatively slowly near the surface. The inversion results of the large model reflect the actual distribution law of the stress field and the influences of topography and geomorphology.Figure 9
Principal stress contour of the sectionx = 135 m for the large model. (a) Maximum principal stress contour of section x = 135 m for the large model. (b) Minimum principal stress contour of section x = 135 m for the large model.
(a)(b)


### 4.2. Submodel Stress Field Inversion Analysis

The key area with faults was taken as the submodel within the larger model area, as shown in Figure5. The submodel was calculated in the range x = 0–500 m and y = −100–400 m, and the elevation of the bottom of the submodel was 2,600 m. The submodel has a total of 8 faults (F2, F3, f80, f26, f23, f9, f85, and f28), as shown in Figures 10(a) and 10(b).Figure 10
Submodel with faults. (a) Meshes of the submodel. (b) Perspective view of submodel with faults.
(a)(b)According to Step 1 of Section2.3, the stress field obtained from the large model is loaded into the submodel by means of the interpolation procedure. A new submodel stress field is obtained after the stress equilibrium calculation.Through Step 2 of Section2.3, the stress components at the measured points are replaced by the measured values.Through Steps 3 and 4, each of the 5 measured points is used as the center point, and the data are extracted for BP neural network training. There are 11 input data points (hn, dxn1, dyn1, dxn2, dyn2, dxn3, dyn3, dxn4, dyn4, dxn5, and dyn5) and 6 target data points (σxx, σyy, σzz, τyz, τxz, and τxy), and 15 nodes are set in the hidden layer of the BP neural network. After training, a corrected stress component prediction model can be obtained.Finally, through Step 5 of Section2.3, the relative position of all elements is converted into the form of equation (4) and imported into the trained BP neural network to calculate the stress components of all elements. After the stress equilibrium calculation in FLAC3D, the final inversion result of the in situ stress field for the submodel is obtained.The results of the LSFC method of the submodel are shown in Table6. According to statistical data, the error of measurement results can reach 25%–30% [46]. Therefore, as long as the inversion error of the measured points is less than 30%, the inversion results are considered to be in line with the actual stress field distribution. The relative errors of the submodel in situ stress field obtained by the LSFC method are all less than 30%. In comparison with the large model, the in situ stress field of the submodel is more accurate.Table 6
Comparison of the measured stress component values and those calculated by the large model and the submodel.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—Large model10.8812.009.05−3.67−0.951.3417.92Submodel10.3313.178.84−3.25−0.13−1.4813.282Measured stress8.7717.825.72−2.85−1.61−1.07—Large model8.9113.606.84−3.870.810.8125.95Submodel9.9315.665.28−1.52−1.210.4815.513Measured stress15.2922.3612.015.026.17−6.43—Large model13.8919.299.47−2.244.70−2.2330.30Submodel15.5618.726.951.833.04−3.7025.994Measured stress14.5112.0010.06−5.694.52−6.06—Large model11.6014.007.86−2.153.50−2.2228.92Submodel13.7813.409.42−4.613.84−4.8710.445Measured stress12.8819.8912.47−1.772.32−6.68—Large model13.7016.3710.63−1.701.74−1.9422.57Submodel13.5416.0311.25−1.752.36−2.9320.02Figures11(a) and 11(b) show the maximum and minimum principal stress contours for the x = 135 m section, respectively. The figure shows the following:(1)
The range of the maximum principal stress in the main powerhouse is 17.5–25 MPa, and the range of minimum principal stress in the main powerhouse is 7.5–15 MPa.(2)
Because of the weak mechanical parameters of the faults, stress release occurs within the faults. Faults have a large effect on the local stress values, causing local stress elevation at the fault edge due to extrusion.(3)
The maximum and minimum principal stress values tend to increase with depth, but the distribution of the stress field is also influenced by tectonic movement, gravity, and topographic relief, in agreement with the results analyzed in Section3.3. The results of this inversion can be used to guide underground engineering design and construction.Figure 11
Principal stress contour of sectionx = 135 m for submodel. (a) Maximum principal stress contour of section x = 135 m for submodel. (b) Minimum principal stress contour of section x = 135 m for submodel.
(a)(b)


## 4.1. Large Model Stress Field Inversion Analysis by the LSFC Method

The size of the large model in the horizontal plane is 1100 m × 1100 m, and the coordinate position and orientation are shown in Figure5. The numerical model in FLAC3D is shown in Figure 8. The z-axis is consistent with the geodetic coordinates along the vertical direction.Figure 8
Meshes of the large model.According to Step 1 of Section 2.2.2, first, the unit displacement is loaded to the boundary, as shown in Figure1, to obtain the stress distribution under the basic influencing factors. Using equation (1) as the fitness function, γk of each stress component was obtained by the GA algorithm in MATLAB version 2016a. The calculated results are γ1 = 0.951, γ2 = 1.019, γ3 = 1.214, γ4 = −8.66, γ5 = 14.47, and γ6 = 32.86. The stress components of each element in the numerical model can be obtained from equation (2). Then, through Step 2 of 2.2.2, the stress field of the initial inversion is obtained after the stress equilibrium calculation. The results are shown in Table 5, and the relative error is defined as (5)Δ=Calculated value−Measured value2Measured value2×100%,where 2 is the 2-norm.Table 5
Comparison of measured stress component values calculated by GA and LSFC.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—GA12.4716.477.35−1.902.25−4.3127.09LSFC10.8812.009.05−3.67−0.951.3417.922Measured stress8.7717.825.72−2.85−1.61−1.07—GA12.5715.406.63−1.762.09−3.9731.84LSFC8.9113.606.84−3.870.810.8125.953Measured stress15.2922.3612.015.026.17−6.43—GA12.1215.1911.38−1.902.82−1.7638.13LSFC13.8919.299.47−2.244.70−2.2330.34Measured stress14.5112.0010.06−5.694.52−6.06—GA10.7517.6510.71−2.612.99−2.4036.24LSFC11.6014.007.86−2.153.50−2.2228.925Measured stress12.8819.8912.47−1.772.32−6.68—GA13.2912.8811.00−1.803.25−1.2932.51LSFC13.7016.3710.63−1.701.74−1.9422.57An approximate in situ stress field was calculated by the GA method with a relative error of 27.09% to 38.13%, which can be further optimized by the LSFC method.Through Step 3 of Section2.2.2, the stress component values at the measured points are replaced by the measured values.Through Steps 4 and 5 of Section2.2.2, each of the 5 measured points is used as the center point, and the data are extracted for BP neural network training. Since there are 5 measured points, the column number of equation (4) is 11. There are 11 input data points (hn, dxn1, dyn1, dxn2, dyn2, dxn3, dyn3, dxn4, dyn4, dxn5, and dyn5) and 6 target data points (σxx, σyy, σzz, τyz, τxz, and τxy), and 15 nodes are set in the hidden layer of the BP neural network. After training, a corrected stress component prediction model can be obtained.Finally, through Step 6 of Section2.2.2, the relative position of all elements is converted into the form of equation (4) and imported into the trained BP neural network to calculate the stress components of all elements. After performing the stress equilibrium calculation in FLAC3D, a more accurate stress field was calculated.The results of the LSFC method are shown in Table5. In comparison with the GA method, the LSFC method was more accurate. The x = 135 m section is in the direction of the main powerhouse axis and passes through measured points 1 and 2, as shown in Figure 5. Figures 9(a) and 9(b) show the maximum and minimum principal stress contours for the x = 135 m section, respectively. The maximum principal stresses are more densely distributed and have a large gradient in the river valley. The distribution of the maximum and the minimum principal stresses in the central area of the model away from the boundary increases with increasing burial depth and changes with topography. The gradient of the minimum principal stress changes relatively slowly near the surface. The inversion results of the large model reflect the actual distribution law of the stress field and the influences of topography and geomorphology.Figure 9
Principal stress contour of the sectionx = 135 m for the large model. (a) Maximum principal stress contour of section x = 135 m for the large model. (b) Minimum principal stress contour of section x = 135 m for the large model.
(a)(b)


## 4.2. Submodel Stress Field Inversion Analysis

The key area with faults was taken as the submodel within the larger model area, as shown in Figure5. The submodel was calculated in the range x = 0–500 m and y = −100–400 m, and the elevation of the bottom of the submodel was 2,600 m. The submodel has a total of 8 faults (F2, F3, f80, f26, f23, f9, f85, and f28), as shown in Figures 10(a) and 10(b).Figure 10
Submodel with faults. (a) Meshes of the submodel. (b) Perspective view of submodel with faults.
(a)(b)According to Step 1 of Section2.3, the stress field obtained from the large model is loaded into the submodel by means of the interpolation procedure. A new submodel stress field is obtained after the stress equilibrium calculation.Through Step 2 of Section2.3, the stress components at the measured points are replaced by the measured values.Through Steps 3 and 4, each of the 5 measured points is used as the center point, and the data are extracted for BP neural network training. There are 11 input data points (hn, dxn1, dyn1, dxn2, dyn2, dxn3, dyn3, dxn4, dyn4, dxn5, and dyn5) and 6 target data points (σxx, σyy, σzz, τyz, τxz, and τxy), and 15 nodes are set in the hidden layer of the BP neural network. After training, a corrected stress component prediction model can be obtained.Finally, through Step 5 of Section2.3, the relative position of all elements is converted into the form of equation (4) and imported into the trained BP neural network to calculate the stress components of all elements. After the stress equilibrium calculation in FLAC3D, the final inversion result of the in situ stress field for the submodel is obtained.The results of the LSFC method of the submodel are shown in Table6. According to statistical data, the error of measurement results can reach 25%–30% [46]. Therefore, as long as the inversion error of the measured points is less than 30%, the inversion results are considered to be in line with the actual stress field distribution. The relative errors of the submodel in situ stress field obtained by the LSFC method are all less than 30%. In comparison with the large model, the in situ stress field of the submodel is more accurate.Table 6
Comparison of the measured stress component values and those calculated by the large model and the submodel.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—Large model10.8812.009.05−3.67−0.951.3417.92Submodel10.3313.178.84−3.25−0.13−1.4813.282Measured stress8.7717.825.72−2.85−1.61−1.07—Large model8.9113.606.84−3.870.810.8125.95Submodel9.9315.665.28−1.52−1.210.4815.513Measured stress15.2922.3612.015.026.17−6.43—Large model13.8919.299.47−2.244.70−2.2330.30Submodel15.5618.726.951.833.04−3.7025.994Measured stress14.5112.0010.06−5.694.52−6.06—Large model11.6014.007.86−2.153.50−2.2228.92Submodel13.7813.409.42−4.613.84−4.8710.445Measured stress12.8819.8912.47−1.772.32−6.68—Large model13.7016.3710.63−1.701.74−1.9422.57Submodel13.5416.0311.25−1.752.36−2.9320.02Figures11(a) and 11(b) show the maximum and minimum principal stress contours for the x = 135 m section, respectively. The figure shows the following:(1)
The range of the maximum principal stress in the main powerhouse is 17.5–25 MPa, and the range of minimum principal stress in the main powerhouse is 7.5–15 MPa.(2)
Because of the weak mechanical parameters of the faults, stress release occurs within the faults. Faults have a large effect on the local stress values, causing local stress elevation at the fault edge due to extrusion.(3)
The maximum and minimum principal stress values tend to increase with depth, but the distribution of the stress field is also influenced by tectonic movement, gravity, and topographic relief, in agreement with the results analyzed in Section3.3. The results of this inversion can be used to guide underground engineering design and construction.Figure 11
Principal stress contour of sectionx = 135 m for submodel. (a) Maximum principal stress contour of section x = 135 m for submodel. (b) Minimum principal stress contour of section x = 135 m for submodel.
(a)(b)


## 5. Multiple Regression Method for Inversion of in Situ Stress Field and Comparison with LSFC Method

### 5.1. Multiple Regression Method for Inversion of In Situ Stress Field

Multiple regression methods have been widely used in different engineering projects in the past. This type of method mainly simulates the tectonic stresses by applying unit loads on the boundary and constructs the regression equation through the principle of superimposed stresses. The in situ stress field inversion of the large model in Section2.1 is based on the idea of the multiple regression method. As shown in Figure 1, the stress field is caused by gravity, extrusion, and shearing tectonic movement in the multiple regression method. The multiple regression method of in situ stress can be expressed as(6)σCS=a1σTS1+a2σTS2+a3σTS3+a4σTS4+a5σTS5+a6σTS6+a7σTS7+b,where σCS is the calculated in situ stress, a1–a7 are the coefficients, σTS1 is the stress caused by a 1-meter uniform extrusion displacement Ux along the x-axis applied in the y-z plane, as shown in Figure 1(a), σTS2 is the stress caused by a 1-meter uniform extrusion displacement Uy along the y-axis applied to the x-z plane, as shown in Figure 1(b), σTS3 is the stress caused by gravity G, as shown in Figure 1(c), σTS4 is the stress caused by a 1-meter tangential distribution of displacement Uyz along the z-axis applied to the x-z plane, as shown in Figure 1(d), σTS5 is the stress caused by a 1-meter tangential distribution of displacement Uxz along the z-axis applied to the y-z plane, as shown in Figure 1(e), and σTS6 and σTS7 are the stresses caused by a 1-meter tangential distribution of displacement Uxy along the y-axis applied to the y-z plane and along the x-axis applied to the x-axis plane, respectively. The above loading is applied to the submodel in Section 4.2 in FLAC3D.After calculating the tectonic stress with respect to the different individual factors, the measured stress value is used as the target of the regression equation. The coefficients of the regression equation can be solved by the least square method, and the coefficients are imported into equation (6) to obtain the regression equation for the in situ stress field of Yebatan:(7)σCS=26.7σTS1+47.1σTS2−72.2σTS3+7.9σTS4−0.6σTS5+10.4σTS6−11.3σTS7−12.2.The regression residual sum of squares of equation (7) is Q = 356.975, the multicorrelation coefficient is R = 0.91, the regression sum of squares is U = 1748.41, and the F-statistic shows F = 15.39, which exceeds the critical value F (7, 30–7-1) = 2.46 at the significance level of 0.05. Therefore, it can be proven that the linear relationship of the variables of the regression equation is significant. R = 0.91 indicates that the regression equation (7) has a good correlation.After the stress equilibrium calculation in FLAC3D, the inversion result calculated by the multiple regression method is obtained. The comparison between the stresses measured and calculated by multiple regression methods and the LSFC method is shown in Table 7.Table 7
Comparison between measured stress and stress calculated by multiple regression methods and the LSFC method.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—Multiple regression13.3315.795.900.472.192.0428.62LSFC10.3313.178.84−3.25−0.13−1.4813.282Measured stress8.7717.825.72−2.85−1.61−1.07—Multiple regression15.9016.997.930.393.380.5945.80LSFC9.9315.665.28−1.52−1.210.4815.513Measured stress15.2922.3612.015.026.17−6.43—Multiple regression19.4818.2617.593.197.572.9840.27LSFC15.5618.726.951.833.04−3.7025.994Measured stress14.5112.0010.06−5.694.52−6.06—Multiple regression12.0310.004.681.411.321.5447.96LSFC13.7813.409.42−4.613.84−4.8710.445Measured stress12.8819.8912.47−1.772.32−6.68—Multiple regression21.3319.8521.562.373.092.2755.27LSFC13.5416.0311.25−1.752.36−2.9320.02


### 5.2. Comparison and Analysis of the Multiple Regression and LSFC Methods

The comparison of the multiple regression and LSFC methods is shown in Table8 and Figure 12. The relative error of multiple regression is greater than 30% for all 4 measured points except point 1, and all values are greater than the relative error of the LSFC method. The average relative error of LSFC is 17.05%, but the average relative error of the multiple regression method reaches as high as 43.58%. Therefore, the multiple regression method is very ineffective under complex geological conditions and is not applicable to Yebatan engineering. By comparison with the multiple regression method and through analyzing the stress field distribution of the inversion results in Section 4.2, the advantages of LSFC compared to the existing method are determined to be as follows:(1)
The LSFC method does not require a linear relationship between stress and depth, so it has much broader applicability. Although there is no obvious linear relationship between the stress and the buried depth in the Yebatan area (Figures6 and 7), the relative errors of the LSFC method are below 30% and much less than the relative errors of the multiple regression method.(2)
The LSFC method can more effectively reflect the influence of tectonic movements on the stress distribution. As shown in Figure11, the stress distribution in the Yebatan area does not strictly increase with depth and is mainly influenced by tectonic movements.Table 8
Comparison of the average relative errors of the multiple regression and LSFC methods.
MethodAverage relative error (%)Multiple regression43.58LSFC17.05Figure 12
Comparison of the relative errors of multiple regression and the LSFC method.


## 5.1. Multiple Regression Method for Inversion of In Situ Stress Field

Multiple regression methods have been widely used in different engineering projects in the past. This type of method mainly simulates the tectonic stresses by applying unit loads on the boundary and constructs the regression equation through the principle of superimposed stresses. The in situ stress field inversion of the large model in Section2.1 is based on the idea of the multiple regression method. As shown in Figure 1, the stress field is caused by gravity, extrusion, and shearing tectonic movement in the multiple regression method. The multiple regression method of in situ stress can be expressed as(6)σCS=a1σTS1+a2σTS2+a3σTS3+a4σTS4+a5σTS5+a6σTS6+a7σTS7+b,where σCS is the calculated in situ stress, a1–a7 are the coefficients, σTS1 is the stress caused by a 1-meter uniform extrusion displacement Ux along the x-axis applied in the y-z plane, as shown in Figure 1(a), σTS2 is the stress caused by a 1-meter uniform extrusion displacement Uy along the y-axis applied to the x-z plane, as shown in Figure 1(b), σTS3 is the stress caused by gravity G, as shown in Figure 1(c), σTS4 is the stress caused by a 1-meter tangential distribution of displacement Uyz along the z-axis applied to the x-z plane, as shown in Figure 1(d), σTS5 is the stress caused by a 1-meter tangential distribution of displacement Uxz along the z-axis applied to the y-z plane, as shown in Figure 1(e), and σTS6 and σTS7 are the stresses caused by a 1-meter tangential distribution of displacement Uxy along the y-axis applied to the y-z plane and along the x-axis applied to the x-axis plane, respectively. The above loading is applied to the submodel in Section 4.2 in FLAC3D.After calculating the tectonic stress with respect to the different individual factors, the measured stress value is used as the target of the regression equation. The coefficients of the regression equation can be solved by the least square method, and the coefficients are imported into equation (6) to obtain the regression equation for the in situ stress field of Yebatan:(7)σCS=26.7σTS1+47.1σTS2−72.2σTS3+7.9σTS4−0.6σTS5+10.4σTS6−11.3σTS7−12.2.The regression residual sum of squares of equation (7) is Q = 356.975, the multicorrelation coefficient is R = 0.91, the regression sum of squares is U = 1748.41, and the F-statistic shows F = 15.39, which exceeds the critical value F (7, 30–7-1) = 2.46 at the significance level of 0.05. Therefore, it can be proven that the linear relationship of the variables of the regression equation is significant. R = 0.91 indicates that the regression equation (7) has a good correlation.After the stress equilibrium calculation in FLAC3D, the inversion result calculated by the multiple regression method is obtained. The comparison between the stresses measured and calculated by multiple regression methods and the LSFC method is shown in Table 7.Table 7
Comparison between measured stress and stress calculated by multiple regression methods and the LSFC method.
Measured pointMethodStress components (MPa)Relative error (%)σxxσyyσzzτyzτxzτxy1Measured stress12.7813.289.14−4.410.44−1.22—Multiple regression13.3315.795.900.472.192.0428.62LSFC10.3313.178.84−3.25−0.13−1.4813.282Measured stress8.7717.825.72−2.85−1.61−1.07—Multiple regression15.9016.997.930.393.380.5945.80LSFC9.9315.665.28−1.52−1.210.4815.513Measured stress15.2922.3612.015.026.17−6.43—Multiple regression19.4818.2617.593.197.572.9840.27LSFC15.5618.726.951.833.04−3.7025.994Measured stress14.5112.0010.06−5.694.52−6.06—Multiple regression12.0310.004.681.411.321.5447.96LSFC13.7813.409.42−4.613.84−4.8710.445Measured stress12.8819.8912.47−1.772.32−6.68—Multiple regression21.3319.8521.562.373.092.2755.27LSFC13.5416.0311.25−1.752.36−2.9320.02


## 5.2. Comparison and Analysis of the Multiple Regression and LSFC Methods

The comparison of the multiple regression and LSFC methods is shown in Table8 and Figure 12. The relative error of multiple regression is greater than 30% for all 4 measured points except point 1, and all values are greater than the relative error of the LSFC method. The average relative error of LSFC is 17.05%, but the average relative error of the multiple regression method reaches as high as 43.58%. Therefore, the multiple regression method is very ineffective under complex geological conditions and is not applicable to Yebatan engineering. By comparison with the multiple regression method and through analyzing the stress field distribution of the inversion results in Section 4.2, the advantages of LSFC compared to the existing method are determined to be as follows:(1)
The LSFC method does not require a linear relationship between stress and depth, so it has much broader applicability. Although there is no obvious linear relationship between the stress and the buried depth in the Yebatan area (Figures6 and 7), the relative errors of the LSFC method are below 30% and much less than the relative errors of the multiple regression method.(2)
The LSFC method can more effectively reflect the influence of tectonic movements on the stress distribution. As shown in Figure11, the stress distribution in the Yebatan area does not strictly increase with depth and is mainly influenced by tectonic movements.Table 8
Comparison of the average relative errors of the multiple regression and LSFC methods.
MethodAverage relative error (%)Multiple regression43.58LSFC17.05Figure 12
Comparison of the relative errors of multiple regression and the LSFC method.


## 6. Conclusion

From the inversion results, it is clear that the in situ stress field in the Yebatan Hydropower Station area is affected by tectonic stress, gravity, and topography, which is consistent with the results analyzed in Section3.3. Based on the genetic algorithm, BP neural network, and FLAC3D, an LSFC method was proposed for in situ stress field inversion. The effect of the fault on the stress field in the key area has also been analyzed by means of the submodel secondary inversion method. The following conclusions can be drawn:(1)
The proposed local stress field correction method (LSFC) can effectively set up the mapping relationship between the stress components and the relative positions. By a refined submodel with faults, the influence of faults is considered, and the relative errors of the inversion results are further decreased versus those of the large model without faults. The final errors of Yebatan in situ stress are in the range of 10.44%–25.99%, which is less than the measurement error, indicating that the above method is feasible and can satisfy the inversion requirement of the in situ stress field.(2)
The burial depth is generally considered to be the main influencing factor of the stress distribution in past studies, so it is a common practice to invert the in situ stress based on the burial depth. However, the results show that the distribution of in situ stress in the Yebatan area is mainly influenced by tectonic movements. The newly proposed LSFC method does not require a linear relationship between stress and depth, so it has much wider applicability.(3)
By comparison with the multiple regression method, it is found that the average relative error of LSFC inversion results is 17.05%, which is much lower than that of 43.58% for the multiple regression method, indicating that the LSFC method can better reflect the influence of tectonic movements on the in situ stress distribution and has higher inversion accuracy under complex geological conditions.(4)
Due to the weak mechanical parameters of the fault, a stress release occurs within the fault. The local stress values at the edges of the fault increase to some extent due to extrusion. This indicates that the fault has some influence on the local stress values.



---

*Source: 4396168-2021-10-13.xml*