Inversion of horizontal in situ stress field based on wide-azimuth seismic data

In situ stress is a significant rock mechanics parameter, which is widely used in many aspects of the petroleum industry. Horizontal in situ stress is an important part of in situ stress. Hence, how to obtain accurate information as to horizontal stress can vary with azimuth has been widely explored. A new prediction method for in situ stress on a horizontal plane is proposed here. This new method analyses the relationship between seismic data, fracture orientation and a horizontal in situ stress field. Through a series of actual data, it is found that azimuth seismic data and azimuth elastic modulus can effectively point out the direction of the horizontal stress. Meanwhile, the magnitude of horizontal stress is shown to be closely related to a combination of rock elastic modulus and seismic amplitude according to a further analysis of log and wide-azimuth seismic data. Under this premise, an objective function of horizontal in situ stress field inversion using azimuth seismic data is established. Finally, the inversion of a horizontal in situ stress field based on wide-azimuth seismic data is realized. This inversion method can simultaneously invert the maximum horizontal stress, minimum horizontal stress and their orientations, which optimizes the prediction process of a horizontal in situ stress field and establish a new method for the prediction of horizontal in situ stress fields. The application of actual data also verifies the prediction accuracy of this method.


Introduction
The in situ tress prediction of a subsurface medium is an important step in the process of petroleum exploration and development. From the geological point of view, in situ stress provides the driving force for the formation, migration and accumulation of petroleum, and ultimately controls the location of reservoirs. In the process of logging, in situ stress provides a basis for calculating wellbore stability (Raoof et al. 2015). During the development stages, in situ stress determines the fracturing ability of reservoir and indicates the engineering 'sweet spot ' . In general, in situ stress can be expressed on the basis of the magnitude and orientation of three mutually perpendicular stresses. Among them, vertical stress ( v ) is equal to the sum of the gravity above it, which can be predicted by the gravity integral of the overlying rock (Lei & Sinha 2012). The magnitude and orientation of horizontal in situ stress (HISS, including maximum horizontal in situ stress H , minimum horizontal in situ stress h and their orientations) are affected by many factors such as rock gravity, pore pressure, tectonic stress and thermal stress (Mallick et al. 1998;Wu et al. 2009): so, the prediction is more complex and less accurate. Therefore, how to obtain accurate HISS information is an urgent problem to be solved.
Many scholars have investigated the prediction methods of HISS . These methods can be roughly divided into four categories: laboratory measurement, welllog analysis, numerical simulation and seismic data prediction. Among them, laboratory measurement methods can only obtain small amounts of data. Well-log analysis can only reflect the in situ stress distribution near the well locations, which is not suitable for stress prediction in a large-scale survey (Campbell-Stone et al. 2012). The numerical simulation method has high requirements for the construction of the model and the application of loads and boundary conditions, leading to poor feasibility. Compared with the other three methods, seismic data are less difficult to obtain and cover more information. Thus, the method of using seismic data to predict stress in reservoir develops. Dillen (2000) investigated the seismic characteristic change with stress field and researched the earthquake mechanisms by observing timelapse seismic data. This thesis set up the foundation for predicting in situ stress with seismic data. Hunt et al. (2011) obtained elastic modulus and curvature characteristics from seismic data, and built relevant models to realize an HISS prediction. Gray et al. (2012) used 3D seismic data to calculate the magnitude and orientation of H and h . Based on a transversely isotropic medium, a new stress prediction model was established by calculating the medium stiffness matrix (Zhang et al. 2015). Information about elastic parameters and anisotropic parameters of underground media was obtained by using the elastic impedance inversion of prestack seismic data . Then they built the HISS prediction model of relevant parameters, and realized the HISS prediction.
The HISS prediction methods based on seismic data are mainly model-driven prediction methods, such as the Nick empirical equation, Andersen model, Huang's model etc. (Ge et al. 1998). These prediction methods invert various parameters related to stress through seismic data, and then substitute these parameters into the prediction model to realize the prediction of HISS. The limitations of these methods include, on the one hand, according to the structure of stress model, restrictions to the applicable scope. For example, some models do not distinguish between maximum and minimum horizontal stress, so they are unsuitable for an area where tectonic movement is obvious. On the other hand, due to the uncertainty of the inversion method itself, the error generated in the inversion process will be transferred and multiplied through the model (Zhang et al. 2020), and ultimately affect the prediction accuracy of HISS. To get rid of the limitation of the prediction model and reduce the cumu- lative error of the inversion process, a new HISS inversion method driven by wide-azimuth seismic data is proposed. First, according to the statistical analysis, a stress-related attribute Q is constructed to indicate the characteristics of HISS changing with azimuth. Then, based on the analysis of rock mechanics, the objective equation of direct inversion of HISS field (including H and h and their orientations) are derived. The equation establishes the direct relationship between seismic elastic information and rock mechanics information, which makes it possible to directly inverse the HISS field of underground medium. Finally, the inversion of HISS field is realized by this new inversion method. Working area data application illustrates that accurate information of horizontal stress field can be obtained from this inversion method.

Orientation of horizontal stress indicated by wide-azimuth seismic data
In the last few years, acquisition and processing technology of wide-azimuth seismic data has been greatly developed, and wide-azimuth seismic data have been gradually applied to many aspects of seismic inversion (Pan et al. 2018). Wideazimuth data expand the azimuth of the observation system, and obtain rich information about the reflection waves of a subsurface medium varying with the azimuth, which is infeasible for narrow-azimuth seismic data. Therefore, wideazimuth seismic data have a significant advantage in identifying geological information with horizontal anisotropy (Pan et al. 2017). Fracture detection technology based on azimuth seismic data is a good example, which reflects its advantages (Mallick & Frazer 1991). As figure 1 shows, fracture detection technology uses an ellipse fitting method (Mitchell & Van 2016) to analyze the seismic amplitude attributes varying with azimuth, and obtains information about major and minor axes of the amplitude attributes ellipse, which can be used to indicate the azimuth of fractures. Fracture occurs when the stress of subsurface medium exceeds its own bearing capacity. As a result, there is a close 132 Downloaded from https://academic.oup.com/jge/article/19/2/131/6551130 by guest on 24 March 2022 relationship between the extension direction of natural fractures in subsurface media and the orientation of horizontal stress (Hubbert & Willis 1957;Mallick et al. 2017;Zhang et al. 2021). As figure 2 shows, the direction of fracture extension can indicate the orientation of H : in the two main types of rock fracture mechanism, the tensile fracture is parallel to the orientation of the H and the acute angle bisector of a group of shear fractures is also consistent with the orientation of H (Haimson 2006). Logically, the long axis direction of azimuth seismic amplitude attribute can also indicate the orientation of H as figures 1 and 2 show.
Formation micro-scanner image (FMI) gets the distribution of horizontal stress around the wellbore by scanning and analyzing the surrounding resistivity (Huang et al. 2012). Wellbore breakout and induced fractures are common in drilling, which can be observed by FMI. These events occur because of the imbalance between the confining pressure and the pressure of drilling fluid during drilling. The orientation of horizontal stress at the drilling position can be indicated by analyzing the azimuth of these events (Zoback 2007). At the wellbore breakout, the H is perpendicular to the long axis of the wellbore, and the induced fractures are parallel to the H (Nian et al. 2016).
Figures 3a and 4a are two imaging logging profiles with obvious wellbore breakout and induced fractures, respectively, in the survey. Figures 3b and 4b are the rose diagrams of the wellbore breakout (or induced fracture) azimuth of FMI profiles. Figures 3c and 4c show the ellipse fitting results of azimuth amplitude attributes at corresponding well locations. Through analysis, it can be found that the wellbore breakout in figure 3a is about 160-180°and the orientation of H is perpendicular to it, about 70-90°, which is consistent with the long axis orientation of ellipse fitting in figure 3c. In figure 4b, the orientation of induced fracture is about 90-110°, indicating the orientation of H , about 90-110°, which matches the long axis of the amplitude ellipse in figure 4c well. These figures demonstrate the ratiocination that the amplitude attribute of azimuth seismic data can indicate the orientation of horizontal stress mentioned before.
It is universal to describe formation anisotropy by using the ellipse characteristics of azimuth seismic data (Downton   & Gray 2006). However, the different combination of the properties of the media above and below the stratigraphic interface leads to the difference of reflection characteristics. Ignoring this difference and considering that the major and minor axes indicate the orientations of H and h , respectively, increases the uncertainty of HISS field prediction, that is, the so-called 90°ambiguity (Li 2013). Therefore, it is controversial to obtain the orientation of horizontal stress only by using azimuthal amplitude attribute. The elastic modulus is an inherent characteristic of rock that describes its deformation resistance. It depends on the internal structure of the rock, mineral composition, density and arrangement of pores. Therefore, the elastic modulus has obvious azimuth characteristics. Normally, the azimuth elastic modulus parallel to the orientation of H is larger than the azimuth elastic modulus parallel to the orientation of h (Sayers 2013). Obviously, in the process of HISS prediction, adding azimuth elastic modulus can effectively avoid the uncertainty and make the prediction result more reasonable.

Magnitude of horizontal stress predict by wide-azimuth seismic data
According to Hooke's law, the magnitude of stress is related to the elastic modulus and strain. Meanwhile, the amplitude attribute of seismic data reflects the displacement of a medium under stress. In the process of energy transfer of seismic wave, rock displacement is caused by overcoming the confining pressure. Therefore, the information of seismic amplitude and elastic modulus can be used to represent the magnitude of horizontal stress indirectly.
To avoid the interference of positive and negative transformation of seismic amplitude, the root mean square (RMS) amplitude attribute is used in this study, and the RMS amplitude can be obtained by equation (1). All descriptions of amplitude in the following paper refer to RMS amplitude.
where A represents RMS amplitude, n represents the number of signal and a i represents amplitude of the ith signal.
To further clarify the relationship between amplitude, elastic modulus and in situ stress, some log and seismic data in the working area were analyzed. Figure 5a and b shows the regression analysis cross plot of wells 13 and 133 in the survey. The x-and y-axes in the figures represent the elastic modu-lus and the measured stress, respectively. The scatter points of different colors in the figure represent the data of different layers in the well logs. The straight lines in different colors are the statistical regression results of the corresponding color scatters. It can be seen that the scatter points of different colors are distributed near the straight line.
From figure 5a and b, it is clear that there is a linear relation between normal stress and elastic modulus. However, different lines (represent different layers) have different intercepts and gradients. Taking seismic amplitude into consideration, we get a more unified result which is shown in figure 6a and b. In these figures, the x-axis is replaced by a product of elastic modulus and amplitude, the y-axis remains unchanged. Here, 10 -8 makes the x-and y-axes values have the same units, MPa. Figure 6 reflects that the gradients of different regression lines become generally consistent. Also, through analyses of the intercept P and the travel time of corresponding layers, it is found that P increases linearly with the travel time. Table 1 shows the travel time interval of different layers and P.
Therefore, it is assumed that stress can be expressed by seismic and rock elastic information as shown in equation (2): where is the normal stress; E, A represent elastic modulus and amplitude, respectively; k is the scale factor to unified the demotions, which in this study is 10 -8 , and P is the intercept, which increases as the layers become deeper. Through equation (2), we can use the information of seismic amplitude and elastic modulus to indirectly express the magnitude of horizontal stress.

Inversion of horizontal in situ stress field from wide-azimuth seismic data
Under equation (2), a new stress attribute Q varying with azimuth data is defined here: where is the azimuth of the data. In equation (3), P can be regarded as a kind of lowfrequency model that increases linearly with the increase in layer travel time. It can be obtained through well data interpolation over the interpreted horizons. Azimuth amplitude attribute can be obtained from wide-azimuth seismic data according to equation (1). Meanwhile, the azimuthal elastic modulus is obtained by inversion of prestack azimuthal seismic data. Ruger (1996) deduces the approximate equation of the relationship between P-wave azimuth reflection coefficient and P-and S-wave velocity, density and three anisotropy coefficients. Based on this, AVAZ technology develops rapidly and is widely used in anisotropic analysis (Xu & Li 2001;Shaw & Sen 2006). According to the connection between seismic elastic parameters and rock mechanics parameters, the reflection coefficient equation related to Poisson's ratio and elastic modulus is established (Zong et al. 2012). Then the theory is extended to azimuth seismic data and constructed the inversion equation to obtain azimuthal elastic modulus (Wang et al. 2021): d( , ) = 1 2 sin 2 cos 2 (1 + sin 2 )(1 − sin 2 cos 2 ), e( , ) = 1 2 sin 4 cos 4 (1 + sin 2 ), where ΔE E , Δ , Δ are the reflection coefficients of elastic modulus, Poisson's ratio and density, , , are three anisotropic parameters (Thomsen 1986), p is square of the ratio of P-and S-wave velocity and , represent the angle of incidence and azimuth, respectively. According to equation (4), at least six azimuthal P-wave reflection coefficient data volumes are required for inversion, and finally the azimuthal elastic modulus information is obtained. After obtaining the amplitude and elastic modulus, the azimuthal Q attribute volume is constructed using equation (3)  In Mohr's circle theory, the normal stress and shear stress in any direction can be expressed by the first and second principal stresses on the plane, as figure 7 shows. is the angle between the stress and normal phase of principal plane. In terms of horizontal stress distribution, the first and second principal stresses are H , h . Then, the normal stress of subsurface medium along any horizontal direction can be expressed as (King et al. 2008;Shahri et al. 2016): According to equations (3) and (5), the relationship between HISS and Q attribute of any azimuth is In equation (6), a nonlinear equation is constructed that contains all the parameters used in HISS expression. Then a linearization process is applied to optimize the nonlinear equation as follows. Define Then equation (6) can be written as In figure 8, there is a certain azimuthal seismic data, and the included angle in the x direction is . In the same coordinate system, we assume the angle between the H and the x direction is . Thus, the angle between the data and the stress orientation is defined as , and = ∅ − according to figure 8. Therefore, Q ( ) = + + − ⋅ cos 2( − ).
(8) Figure 8. Relationship between azimuth seismic data and orientation of horizontal stress.
According to trigonometric transformation, we can get Substitute in to equation (9), then In equation (11), there are three unknown parameters + , M, N. At least three different sets of azimuth data are required to solve this problem.
Finally, the HISS can be expressed as: Equation (11) establishes the potential relationships between azimuth seismic amplitude and H , h and their orientations. The direct inversion of HISS field can be realized to benefit from this equation.
When applied to a data volume, assuming there are m azimuth seismic data, the equation (11) is written in matrix form: 1 cos 2 1 sin 2 1 1 cos 2 2 sin 2 2 . . . Furthermore, for n samples of each seismic trace, there is where Q 1 = (Q 1 1 , Q 2 1 , ...Q n 1 ) T , Equation (14) can be simplified as where The final inversion objective function is established by the damping least-square method to realize the inversion of HISS field (Wang 2016;Zhang & Dai 2016): R is the solution, and is the damping factor or weighting factor, which is used to enhance the accuracy of the inversion result.  The direct inversion process of HISS field is as follows: (i) Obtain the amplitude attribute information of each azimuth from the wide-azimuth seismic data. Meanwhile, the azimuthal elastic modulus is inverted from the same seismic data. (ii) According to equation (3), the amplitude attribute with azimuth information is combined with the azimuthal elastic modulus to obtain the Q attribute of seismic data. (iii) The relationship equation between the azimuthal Q attribute and HISS field is established by equation (11). An objective function with at least three seismic azimuth-stack gathers is found to predict the unknown parameters + , M, N. Finally, the HISS field is inversed by the damping least-square method.

Actual data application
To verify the feasibility and accuracy of the direct inversion method, actual data are applied. Six different azimuth angles of actual seismic data are selected to test the HISS field inversion method in this paper. Figure 9 shows the Q attribute slice of different azimuths at the T2 horizon.  Overall, the value of the Q attribute in figure 9a and f is relatively larger than other slices, which indicates that the in situ stress orientation of the actual data is mainly east to west. However, in some local positions, figure 9c and d show larger Q values, indicating that the stress orientation at these positions is mainly south to north, and the value of Q in the north of the work area is larger than other areas, indicating that stress has an obvious high value here. The slices of Q attribute with different azimuth show pre-inversion information that can provide us a general understanding of HISS. Figure 10a and b shows the maximum and minimum horizontal stress slices of the T2 horizon obtained by the horizontal stress inversion method mentioned previously. These figures show that the stress distribution is consistent with the Q attribute slices, proving the reliability of the inversion method in predicting HISS. Figure 11 shows the orientation of H at T2 horizon obtained by inversion. The direction of the short line in the figure represents the orientation of H , and the different colors of the lines represent the magnitude of the maximum horizontal stress values. Figure 11 parts b and c are enlarged views of the position of the red rectangle in figure 11a. It is clear that the stress orientation in figure 11b is mainly east to west, while the Q values of the 0-30°and 150-180°azimuths in the corresponding positions in figure 9 are larger than other azimuths, so they are consistent. The maximum horizontal stress orientation in figure 11c is mainly in the north to south orientation, while the Q values in the corresponding position of azimuths 60-90°and 90-120°in 139 figure 9 are larger than other azimuth data, which is consistent with the results in figure 11c. These results further prove the accuracy of the horizontal stress inversion method. Figures 12 and 13 are the profiles of the H and the h in the same working area, respectively. The stress logs of well A represent H rmand h are overlapped on the corresponding stress prediction profile separately. Through observing the profile shape and numerical distribution, it is clear that the inversion results of H and h reflect a certain degree of heterogeneity in the transverse direction. In the vertical direction, the inversion results are in good agreement with the actual well logs, which reflects the accuracy of the direct inversion method in horizontal stress prediction.

Conclusion
The prediction of an in situ stress field is a difficult problem in the integration of geological engineering. At present, the prediction methods of in situ stress are still limited to core experiments, logging and other methods. These methods not only have high cost, but also obtain less in situ stress information, which has difficulty in meeting the needs of exploration and development. Wide-azimuth seismic data integrate different pieces of azimuth information and contain sufficient information for analyzing and prospecting. The attribute parameters varying with azimuth can be predicted by using the differences between different azimuth information. Based on this characteristic, a new prediction process is proposed, which sets up a new method for the prediction of HISS fields.
Through the study, we found that: (i) Seismic amplitude ellipse fitting is an effective way to exhibit the orientation of stress subsurface, but the results are not always accurate and may deviate by 90°. The azimuth elastic modulus information is introduced into the orientation prediction of HISS. Through a combination of azimuthal amplitude and azimuthal elastic modulus, the uncertainty of using only the seismic amplitude to predict the orientation of stress is avoided. (ii) Through the statistical analyses of logging data and azimuthal seismic data, the potential associations between HISS, azimuthal seismic amplitude and azimuthal elastic modulus are obtained and presented in the form of equation (3), and then a new stress-related attribute Q is constructed. The azimuthal Q attribute not only indicates the orientation of stress, but also forms a reliable substitute for the magnitude of stress, which can be obtained from wide-azimuth seismic data. (iii) On the basis of rock mechanics analysis, stress in any azimuth can be expressed by HISS. And the stress in any azimuth can be represented as the azimuthal Q attribute. Then, the inversion objective function of pre-dicting HISS field is derived by using the relationship equation. (iv) The H , h and their orientations of actual data are predicted by the direct inversion method. The slices of H , h exhibit a tight connection with Q both in magnitude and orientation, which meet theoretical expectations. Meanwhile, the prediction results of H , h . profiles demonstrate great corresponding to the actual well logs, proving the feasibility of this method in predicting an HISS field.
The method of direct inversion of a horizontal stress field from wide-azimuth seismic data overcomes the limitation of the stress prediction model and avoids the generation of accumulated errors in the process of indirect inversion. Applied to the actual data, the expected prediction effect is achieved.