Horizontal-stress vector prediction method and apparatus, and storage medium

Abstract

A horizontal-stress vector prediction method, which relates to the technical field of in situ stress prediction. The method comprises: calculating elastic parameters and anisotropic parameters of a target reservoir; using a horizon interpretation for the target reservoir as a constraint, and using the calculated elastic parameters and anisotropic parameters to calculate logging curves of Young's modulus, Poisson's ratio and the fracture density; on the basis of the calculated logging curves of Young's modulus, Poisson's ratio and the fracture density, constructing an initial azimuthal elastic impedance model; on the basis of the initial azimuthal elastic impedance model, inverting Young's modulus, Poisson's ratio and the fracture density; substituting into a prediction formula the fracture density and Poisson's ratio, which are obtained from the inversion, and the vertical principal stress of the target reservoir, and performing calculation to obtain the maximum and minimum horizontal stresses of the target reservoir; and on the basis of least-squares ellipse fitting, obtaining the directional distribution of the maximum and minimum horizontal stresses within the target reservoir, so as to realize quick, accurate and highly operable horizontal-stress vector prediction. The present invention further relates to a horizontal-stress vector prediction apparatus, a computer device, and a machine-readable storage medium.

Description

Methods, devices and storage media for predicting horizontal stress vectors

[0001] Cross-references to related applications

[0002] This application claims the benefit of Chinese Patent Application No. 202411515355.4, filed on October 29, 2024, the contents of which are incorporated herein by reference. Technical Field

[0003] This application belongs to the field of geophysical exploration and development technology, specifically relating to a horizontal stress vector prediction method, a horizontal stress vector prediction device, a computer device, and a machine-readable storage medium. Background Technology

[0004] In-situ stress is an important evaluation parameter in the field of oil and gas exploration and development, and it is the most important research content in the large-scale and efficient development of shale gas. The magnitude and direction of the maximum and minimum horizontal stress can provide reliable parameter basis for well location deployment, horizontal well trajectory design, and fracturing parameter optimization in shale gas development.

[0005] Methods for measuring geostress direction include core paleomagnetic methods and differential strain analysis. Paleomagnetic methods separate residual magnetism and determine core orientation through thermal demagnetization and alternating demagnetization. Differential strain analysis obtains stress direction by using linear fitting to calculate the slope of the stress-strain curve from the core experiment (Dong Pingchuan, 2004; Han Jun et al., 2005). Liu Haojuan (2021) established a geostress prediction model. Based on refined 3D seismic interpretation and pre-stack 3D seismic inversion, she used well point data to simulate and select regional adaptability parameters for stress calculation, conducted 3D simulation of the geostress field, and predicted the direction of maximum horizontal stress, the direction of minimum horizontal stress, and the horizontal stress difference coefficient.

[0006] Elastic parameters such as λ, μ, and ρ (λ and μ are Lamé coefficients, and ρ is density), vertical and lateral Poisson's ratios obtained from pre-stack seismic inversion can be used to predict stress-related parameters, such as maximum horizontal stress and closure stress gradient (Starr, 2011). Gray et al. (2012) used parameters closely related to rock mechanics characteristics, such as Young's modulus and Poisson's ratio, based on wide-azimuth seismic data inversion, and obtained stress parameters such as maximum horizontal stress, minimum horizontal stress, and horizontal stress difference ratio to quantitatively determine the stress characteristics of reservoirs. Zhang Guangzhi et al. (2015) analyzed the mineral, porosity, fluid, and anisotropic characteristics of shale reservoirs, established a reliable rock physics model, and carried out shear wave velocity prediction. Then, they predicted the minimum horizontal stress of shale reservoirs through the elastic stiffness tensor. Qi Qing (2018) used thin-plate theory and bending theory to calculate the stress parameters of shale reservoirs based on structure, velocity, and density. Through regional stress analysis, they analyzed the induced fractures and production characteristics during fracturing. The direction and difference coefficient of stress played an important role in the design of horizontal wells. Qu Yang (2019) used resistivity imaging logging and shear wave anisotropy logging to explain the direction of the maximum horizontal stress in the area. He used the Newberry model to calculate the minimum principal stress in the study area and calculated the maximum principal stress based on dual-wellbore analysis, clarifying the trajectory direction of horizontal well deployment in the area and providing a basis for reservoir fracturing. Liu Chang et al. (2019) used finite element simulation to simulate the paleotectonic stress field of tight sandstone in the southern Qinshui Basin, reconstructing the paleostress field under compressive stress. They used the stress field simulation results to construct a comprehensive fracture rate, thereby evaluating the degree of fracture development and predicting the sweet spot of fractures in complex tight gas reservoirs. Zhao Xiaolong et al. (2020) established the relationship between fracture weakness in VTI media and rock mechanical parameters such as Young's modulus, Poisson's ratio, and horizontal in-situ stress. Using shale experimental data for numerical simulation, they concluded that the horizontal stress value decreases with increasing fracture weakness. The approximate equation and simulation results showed a high degree of agreement, effectively improving the prediction of in-situ stress. Lin Haiyu et al. (2022) compared the accuracy of the Huang's (or combined spring) geostress calculation model with the transversely isotropic geostress calculation model, concluding that the transversely isotropic model has higher accuracy. They used this model to calculate the well logging geostress profile of the F group shale oil reservoir in the MH Depression of Xinjiang, providing an important basis for the selection of fracturing intervals and safe drilling. Xu Ke et al. (2022) believed that current geostress has a significant impact on the reservoir quality of the Kuqa Depression in the Tarim Basin. By comprehensively utilizing core test data and well logging data, they established stress concentration parameters expressed by maximum horizontal stress, formation pore pressure, and uniaxial compressive strength, thereby evaluating the reservoir quality.

[0007] However, these technical solutions still need improvement in terms of methodological feasibility, computational complexity, and prediction accuracy. Summary of the Invention

[0008] The purpose of this application is to provide a horizontal stress vector prediction method, a horizontal stress vector prediction device, a computer device, and a machine-readable storage medium to overcome at least one of the defects of existing prediction schemes for ground stress distribution, such as poor feasibility, high computational complexity, and poor prediction accuracy.

[0009] To achieve the above objectives, a first aspect of this application provides a method for predicting horizontal stress vectors, comprising: calculating the elastic parameters and anisotropic parameters of the target reservoir using an anisotropic rock physics model of the target reservoir; using the layer interpretation of the target reservoir obtained through seismic data volume interpretation as a constraint, calculating logging curves of Young's modulus, Poisson's ratio, and fracture density using the calculated elastic parameters and anisotropic parameters; constructing an initial azimuth elastic impedance model based on the calculated logging curves of Young's modulus, Poisson's ratio, and fracture density; and using the initial azimuth elastic impedance model... Based on the initial model, Young's modulus, Poisson's ratio, and fracture density are obtained by inversion within a Bayesian framework. The obtained fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir are substituted into the prediction formula to calculate the maximum and minimum horizontal stress of the target reservoir. The prediction formula uses only the vertical principal stress, fracture density, and Poisson's ratio as independent variables and the maximum and minimum horizontal stress as dependent variables. Least square ellipse fitting is performed on the azimuth stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

[0010] In a specific embodiment of this application, the step of constructing an initial azimuth elastic impedance model based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves includes: calculating the azimuth elastic impedance curve using a first formula based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves; and interpolating the azimuth elastic impedance curves in the seismic data volume to obtain the initial azimuth elastic impedance model; the first formula is:

[0011] Wherein, θ represents the incident angle; φ represents the azimuth angle; E, υ, and e represent Young's modulus, Poisson's ratio, and fracture density, respectively; E0 represents the mean of Young's modulus; υ0 represents the mean of Poisson's ratio; a represents a coefficient related to Young's modulus, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; b represents a coefficient related to Poisson's ratio, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; c represents a coefficient related to fracture density, and is expressed by the incident angle, the azimuth angle, and the square of the P-wave velocity ratio.

[0012] In a specific embodiment of this application, the step of obtaining Young's modulus, Poisson's ratio, and fracture density based on the initial azimuth elastic impedance model under a Bayesian framework includes: performing anisotropic elastic impedance inversion using the constructed target reservoir sub-azimuth stacked gathers and the initial azimuth elastic impedance model to obtain azimuth elastic impedance data volume; and extracting Young's modulus, Poisson's ratio, and fracture density based on the azimuth elastic impedance data volume under a Bayesian framework.

[0013] In a specific embodiment of this application, based on the azimuth elastic impedance data volume, Young's modulus, Poisson's ratio, and crack density are extracted within a Bayesian framework. This includes: logarithmizing the first formula to obtain a linearized azimuth elastic impedance equation; using the matrix form of the linearized azimuth elastic impedance equation as the parameter extraction equation, and extracting Young's modulus, Poisson's ratio, and crack density within a Bayesian framework; wherein, the parameter extraction equation is:

[0014] Where θ1, θ2, and θ3 represent different incident angles, φ1, φ2, and φ3 represent different azimuth angles, EI represents elastic impedance, and EI0 represents the mean value of elastic impedance.

[0015] In a specific embodiment of this application, the prediction formula is:

[0016] Where, σ H σ represents the maximum horizontal stress. V Let υ represent the vertical principal stress, ν represent Poisson's ratio, e represent the crack density, C be a constant, and σ represent the vertical principal stress. h This represents the minimum horizontal stress.

[0017] In a specific embodiment of this application, when performing least-squares elliptic fitting on the azimuth-divided superimposed gathers, the following operations are performed: removing singular points present in the fitting; constraining the solution space of the elliptic fitting using prior information on the geostress direction; and calculating the elliptic fitting error using the least squares median criterion.

[0018] In a specific embodiment of this application, least-squares ellipse fitting is performed on the azimuth stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir. This includes: performing least-squares ellipse fitting on the azimuth stacked gathers of the target reservoir to determine the direction of the major axis and the direction of the minor axis of the ellipse; determining the direction of the major axis of the ellipse as the direction of the maximum horizontal stress within the target reservoir; and determining the direction of the minor axis of the ellipse as the direction of the minimum horizontal stress within the target reservoir.

[0019] In a specific embodiment of this application, there are multiple azimuth elastic impedance initial models, and each azimuth elastic impedance initial model is an azimuth elastic impedance initial model for each target reservoir gather partition with a different offset distance.

[0020] In a specific embodiment of this application, constructing an anisotropic rock physics model of the target reservoir includes: interpreting the logging data of the target reservoir to obtain logging curves of sonic transit time, density, porosity and saturation; and constructing an anisotropic rock physics model of the target reservoir using the logging curves of sonic transit time, density, porosity and saturation.

[0021] In a specific embodiment of this application, the horizontal stress vector prediction method further includes: preprocessing the seismic data volume of the target reservoir and jumping to trigger the execution of inversion of Young's modulus, Poisson's ratio and fracture density based on the initial azimuth elastic impedance model within a Bayesian framework.

[0022] In a specific embodiment of this application, the seismic data volume is a pre-stack seismic data volume, and the preprocessing includes: filtering out noise in the pre-stack seismic data volume using singular value decomposition; removing multiples in the pre-stack seismic data volume using pull transform; equalizing the amplitude energy of different seismic traces in the pre-stack seismic data volume using a spectral balance algorithm; and enhancing the gather resolution in the pre-stack seismic data volume using a compressed sensing algorithm.

[0023] A second aspect of this application provides a horizontal stress vector prediction device, comprising: a first calculation module for calculating the elastic parameters and anisotropic parameters of a target reservoir using a constructed or obtained anisotropic rock physics model of the target reservoir; a second calculation module for calculating logging curves of Young's modulus, Poisson's ratio, and fracture density using the calculated elastic parameters and anisotropic parameters, with the layer interpretation of the target reservoir obtained through seismic data volume interpretation as a constraint; an inversion model construction module for constructing an initial azimuth elastic impedance model based on the calculated logging curves of Young's modulus, Poisson's ratio, and fracture density; and an inversion module for using the azimuth elastic... Based on the initial impedance model, Young's modulus, Poisson's ratio, and fracture density are obtained through inversion within a Bayesian framework. The geostress magnitude calculation module is used to substitute the inverted fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir into the prediction formula to calculate the maximum and minimum horizontal stresses of the target reservoir. The prediction formula uses only the vertical principal stress, fracture density, and Poisson's ratio as independent variables and the maximum and minimum horizontal stresses as dependent variables. The geostress direction prediction module is used to perform least-squares ellipse fitting on the azimuth-stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

[0024] A third aspect of this application provides a computer device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the horizontal stress vector prediction method described in the first aspect of this application.

[0025] A fourth aspect of this application provides a machine-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the horizontal stress vector prediction method described in the first aspect of this application.

[0026] The aforementioned technical solution combines the inversion of Young's modulus, Poisson's ratio, and fracture density based on azimuth elastic impedance, with the direct calculation of maximum and minimum horizontal stresses based on the inverted fracture density and Poisson's ratio. First, the azimuth elastic impedance equation is used to invert elastic parameters and anisotropy parameters, improving the accuracy of the obtained Young's modulus, Poisson's ratio, and fracture density. Then, the improved Poisson's ratio and fracture density are used as parameters for calculating the maximum and minimum horizontal stresses. The calculation process for maximum and minimum horizontal stresses involves fewer parameters, reducing computational complexity and minimizing errors caused by parameter errors, thus improving accuracy. Simultaneously, the directions of maximum and minimum horizontal stresses are obtained through least-squares elliptic fitting. This enables rapid and accurate prediction of the magnitude and direction of reservoir stress, with strong feasibility.

[0027] Other features and advantages of the embodiments of this application will be described in detail in the following detailed description section. Attached Figure Description

[0028] Figure 1 schematically illustrates a first flowchart of a horizontal stress vector prediction method according to an embodiment of this application;

[0029] Figure 2 schematically illustrates a second flowchart of the horizontal stress vector prediction method according to an embodiment of this application;

[0030] Figure 3 schematically illustrates a third flowchart of the horizontal stress vector prediction method according to the above embodiments;

[0031] Figure 4 schematically illustrates a post-stack time offset data volume in a specific application example;

[0032] Figure 5 schematically illustrates the pre-stack OVT gather after optimization in a specific application example;

[0033] Figure 6 schematically illustrates the superimposed cross-section obtained by short-path superposition in a specific application example;

[0034] Figure 7 schematically shows a cross-sectional view of the superimposed mid-path in a specific application example;

[0035] Figure 8 schematically illustrates a composite profile obtained by long-range overlay in a specific application example;

[0036] Figure 9 schematically illustrates the well seismic calibration results of well A3 in a specific application example;

[0037] Figure 10 schematically shows the well seismic calibration results of well A5 in a specific application example;

[0038] Figure 11 schematically illustrates the hierarchical interpretation results in a specific application example;

[0039] Figure 12 schematically shows the azimuth elastic impedance curve of well A3 in a specific application example;

[0040] Figure 13 schematically shows the azimuth elastic impedance curve of well A7 in a specific application example;

[0041] Figure 14 schematically shows the Poisson's ratio inversion profile of well A8 obtained through anisotropic elastic impedance inversion in a specific application example;

[0042] Figure 15 schematically shows the fracture density inversion profile of well A8 obtained through anisotropic elastic impedance inversion in a specific application example;

[0043] Figure 16 schematically shows the Poisson's ratio inversion planar diagram obtained by anisotropic elastic impedance inversion in a specific application example;

[0044] Figure 17 schematically shows the crack density inversion planar diagram obtained by anisotropic elastic impedance inversion in a specific application example;

[0045] Figure 18 schematically illustrates the predicted profile of maximum horizontal stress in well A8 in a specific application example;

[0046] Figure 19 schematically illustrates the minimum horizontal stress prediction profile of well A8 in a specific application example;

[0047] Figure 20 schematically illustrates the maximum horizontal stress vector prediction planar diagram in a specific application example;

[0048] Figure 21 schematically illustrates the minimum horizontal stress vector prediction planar diagram in a specific application example;

[0049] Figure 22 schematically illustrates a structural block diagram of a computer device according to an embodiment of this application. Detailed Implementation

[0050] The specific embodiments of this application will be described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are for illustration and explanation only and are not intended to limit the embodiments of this application.

[0051] If the embodiments of this application involve descriptions such as "first" or "second," these descriptions are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined with "first" or "second" may explicitly or implicitly include at least one of those features. Furthermore, the technical solutions of the various embodiments can be combined with each other, but this must be based on the ability of those skilled in the art to implement them. If the combination of technical solutions is contradictory or impossible to implement, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed in this application.

[0052] The following embodiments illustrate the specific implementation process using the horizontal stress vector prediction method provided in this application to shale reservoirs. It should be understood that this application is not limited to shale reservoirs, meaning that the applicability of the horizontal stress vector prediction method provided in this application to other unconventional and conventional reservoirs is not excluded.

[0053] Figure 1 shows a first flowchart of a horizontal stress vector prediction method according to an embodiment of this application. As shown in Figure 1, the horizontal stress vector prediction method provided by this application includes the following steps 202 to 210.

[0054] Step 202: Calculate the elastic parameters and anisotropic parameters of the target reservoir using the anisotropic rock physics model of the target reservoir.

[0055] In this application, the anisotropic rock physics model of the target reservoir can be constructed using well logging data.

[0056] Step 204: Using the layer interpretation of the target reservoir obtained through the seismic data volume interpretation as a constraint, the logging curves of Young's modulus, Poisson's ratio and fracture density are calculated using the elastic parameters and anisotropy parameters obtained in step 202.

[0057] It is known that the seismic data volume of the target reservoir is usually an optimized seismic data volume through a preprocessing procedure. In this application, the seismic data volume may be an optimized pre-stack seismic data volume.

[0058] Step 206: Construct an initial model of azimuth elastic impedance based on the logging curves of the calculated Young's modulus, Poisson's ratio, and fracture density.

[0059] Step 208: Based on the constructed initial model of azimuth elastic impedance, Young's modulus, Poisson's ratio, and crack density are obtained by inversion within a Bayesian framework.

[0060] Step 210: Substitute the fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir obtained from the inversion into the prediction formula to calculate the maximum and minimum horizontal stress of the target reservoir. The prediction formula uses only the vertical principal stress, fracture density, and Poisson's ratio as independent variables, and the maximum and minimum horizontal stress as dependent variables.

[0061] In this application, the vertical principal stress refers to the pressure generated by the weight of the overlying strata borne by the shale reservoir. The calculation of the vertical principal stress can be combined with the process in the general embodiment. For example, the vertical principal stress can be obtained by integrating the density data of the target reservoir obtained by pre-stack inversion and performing density calculations at various stratum depths.

[0062] As in the above embodiment, the formulas for calculating the maximum and minimum horizontal stresses using only the vertical principal stress, fracture density, and Poisson's ratio as independent variables were first derived. To obtain accurate fracture density and Poisson's ratio data, the azimuth elastic impedance formula was then characterized using Young's modulus, fracture density, and Poisson's ratio as parameters. Using the Bayesian framework's inversion method based on elastic parameters and anisotropic parameters of azimuth elastic impedance, Young's modulus, fracture density, and Poisson's ratio were inverted. The more accurate inversion results were substituted into the formulas for calculating the maximum and minimum horizontal stresses to calculate the magnitude of the geostress in the target reservoir. Therefore, the above embodiment made the following three optimizations to improve the feasibility of the prediction process, reduce the computational complexity of the prediction process, and improve the accuracy of the geostress prediction results: deriving the formulas for calculating the maximum and minimum horizontal stresses using only the vertical principal stress, fracture density, and Poisson's ratio as independent variables; inverting Young's modulus, fracture density, and Poisson's ratio based on azimuth elastic impedance; and directly calculating the magnitude of the geostress using Young's modulus, fracture density, and vertical principal stress.

[0063] In one specific embodiment of this application, step 206, constructing an initial azimuth elastic impedance model based on the well logging curves of the calculated Young's modulus, Poisson's ratio, and fracture density, includes the following steps: using the first formula, calculating the azimuth elastic impedance curve based on the well logging curves of the calculated Young's modulus, Poisson's ratio, and fracture density; and interpolating the azimuth elastic impedance curve in the seismic data volume of the target reservoir to obtain the initial azimuth elastic impedance model.

[0064] The first formula can be expressed as:

[0065] Where θ represents the incident angle; φ represents the azimuth angle; E, υ, and e represent Young's modulus, Poisson's ratio, and fracture density, respectively; E0 represents the mean Young's modulus; υ0 represents the mean Poisson's ratio; a represents a coefficient related to Young's modulus, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; b represents a coefficient related to Poisson's ratio, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; c represents a coefficient related to fracture density, and is expressed by the incident angle, the azimuth angle, and the square of the P-wave velocity ratio. It is known that the mean Young's modulus and the mean Poisson's ratio can be calculated from the well logging data of the target reservoir.

[0066] Based on the above embodiments, step 208, using the constructed initial azimuth elastic impedance model as a basis, inverts Young's modulus, Poisson's ratio, and fracture density under a Bayesian framework, including the following steps: using the constructed target reservoir's azimuth stacked gathers and the initial azimuth elastic impedance model to perform anisotropic elastic impedance inversion, respectively, to obtain azimuth elastic impedance data volume; based on the azimuth elastic impedance data volume, extracts Young's modulus, Poisson's ratio, and fracture density under a Bayesian framework.

[0067] In this application, the target reservoir azimuth stacked gather refers to the gather obtained by stacking gathers by azimuth angle.

[0068] For example, based on the azimuth elastic impedance data volume, Young's modulus, Poisson's ratio, and crack density are extracted under a Bayesian framework, including the following steps: logarithmically transforming the first formula to obtain the linearized azimuth elastic impedance equation; using the matrix form of the linearized azimuth elastic impedance equation as the parameter extraction equation, and extracting Young's modulus, Poisson's ratio, and crack density under a Bayesian framework.

[0069] The parameter extraction equation is as follows:

[0070] In the above parameter extraction equation, θ1, θ2, and θ3 represent different incident angles, φ1, φ2, and φ3 represent different azimuth angles, EI represents elastic impedance, and EI0 represents the mean value of elastic impedance.

[0071] As in the above embodiment, by combining the linear logarithmic transformation of the first formula with the inversion under the Bayesian framework, the stability and prediction accuracy of the geostress earthquake prediction process are improved.

[0072] In one specific embodiment of this application, in step 210, the prediction formula used to calculate the maximum and minimum horizontal stress can be expressed as follows:

[0073] In Formulas 2 and 3, σH σ represents the maximum horizontal stress. V Let υ represent the vertical principal stress, ν represent Poisson's ratio, e represent the crack density, C be a constant, and σ represent the vertical principal stress. h This represents the minimum horizontal stress.

[0074] In one specific embodiment of this application, the elastic parameters calculated by the anisotropic rock physics model of the target reservoir include shear wave velocity, longitudinal wave velocity, and density.

[0075] In one specific embodiment of this application, an anisotropic rock physics model of the target reservoir is constructed using well logging data, including the following steps: interpreting the well logging data of the target reservoir to obtain well logging curves of sonic transit time, density, porosity, and saturation; and constructing an anisotropic rock physics model of the target reservoir using the well logging curves of sonic transit time, density, porosity, and saturation.

[0076] In another specific embodiment of this application, an anisotropic rock physics model of the target reservoir is constructed using well logging data, including the following steps: interpreting the well logging data of the target reservoir to obtain well logging curves of sonic transit time, density, TOC content, porosity, and saturation; and constructing an anisotropic rock physics model of the target reservoir using the well logging curves of sonic transit time, density, TOC content, porosity, and saturation.

[0077] For example, when constructing anisotropic petrophysical models of target reservoirs using well logging curves based on sonic transit time, density, TOC content, porosity, and saturation, the focus is on considering the TOC content, organic porosity, horizontal bedding fractures, and reservoir anisotropy of shale reservoirs. The modulus of the rock matrix is ​​obtained by using the VRH model to calculate the TOC content of quartz, carbonate rocks, and other minerals, while simultaneously using Backus averaging to incorporate clay into the content of other minerals. The modulus of dry rock is obtained by adding organic porosity to the KT model and inorganic porosity to the DEM model. Approximate horizontal bedding fractures are introduced using the Pade approximation equation, and approximate vertical structural fractures are added to the Hudson model, yielding the equivalent modulus of the fractured dry rock skeleton. Anisotropic fluid substitution is performed using the Brown-Korringa model to obtain the equivalent modulus of saturated rock. Finally, shear wave velocity, P-wave velocity, density, and anisotropic parameters are calculated using elastic theory and anisotropic theory.

[0078] In one specific embodiment of this application, to obtain an optimized pre-stack seismic data volume, the following preprocessing steps can be performed: using singular value decomposition to filter out noise in the seismic data volume; using pull transform to remove multiples in the seismic data volume; using spectral balance algorithm to equalize the amplitude energy of different seismic traces in the seismic data volume; and using compressed sensing algorithm to enhance the gather resolution in the seismic data volume.

[0079] Through the above preprocessing, random noise and multiple waves in the angle gathers are eliminated, making the amplitude energy of near and far channels consistent and improving the longitudinal resolution of the gathers.

[0080] In one specific embodiment of this application, the stratigraphic interpretation of the target reservoir is obtained through seismic data volume interpretation, which can employ the process described in the general embodiments. For example, it may include the following steps:

[0081] 1) Well-seismic calibration: The wave impedance curve is calculated using sonic and density logging curves, and the reflection coefficient curve is constructed. The reflection coefficient curve is convolved with the seismic wavelet to obtain the well-side synthetic seismic gather. The synthetic seismic record is subjected to drift and stretching processing. Through multiple iterations, the correlation coefficient between the synthetic seismic gather and the well-side seismic gather reaches the first preset value, and the similarity of waveform and wave group characteristics reaches the second preset value. The correspondence between time and depth obtained at this time is the optimal time-depth relationship. The first preset value is a preset high correlation coefficient, and the second preset value is a preset high similarity coefficient. The specific values ​​of the first and second preset values ​​can be determined according to the application requirements.

[0082] 2) Stratigraphic Interpretation: Using time depth as the vertical standard and seismic phase axis as the horizontal constraint, considering the multilayer rationality of shale reservoirs, and combining with existing geological knowledge, stratigraphic interpretation of the top and bottom of shale reservoirs is carried out in three-dimensional space.

[0083] Example 2

[0084] Figure 2 schematically illustrates a second flowchart of the horizontal stress vector prediction method according to an embodiment of this application. As shown in Figure 2, in one embodiment of this application, the difference between the provided horizontal stress vector prediction method and the first embodiment is that the number of initial azimuth elastic impedance models obtained in step 206 is multiple, and each initial azimuth elastic impedance model is an initial azimuth elastic impedance model for each target reservoir gather partition with a different offset. Specifically, the horizontal stress vector prediction method provided in this embodiment may include steps 202 to 210.

[0085] Step 202: Calculate the elastic parameters and anisotropic parameters of the target reservoir using the anisotropic rock physics model of the target reservoir.

[0086] Step 204: Using the layer interpretation of the target reservoir obtained through the seismic data volume interpretation as a constraint, the logging curves of Young's modulus, Poisson's ratio and fracture density are calculated using the elastic parameters and anisotropy parameters obtained in step 202.

[0087] Step 206: Construct an initial model of azimuth elastic impedance for multiple target reservoir gather zones with different offsets based on the logging curves of Young's modulus, Poisson's ratio and fracture density obtained from calculation.

[0088] Step 208: Based on the constructed initial elastic impedance models for each orientation, the Young's modulus, Poisson's ratio, and fracture density of each target reservoir gather partition are obtained by inversion within the Bayesian framework.

[0089] Step 210: Substitute the fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir obtained from the inversion into the prediction formula to calculate the maximum and minimum horizontal stress of the target reservoir. The prediction formula uses only the vertical principal stress, fracture density, and Poisson's ratio as independent variables, and the maximum and minimum horizontal stress as dependent variables.

[0090] For example, in one specific embodiment, step 206, constructing an initial azimuth elastic impedance model for multiple target reservoir gather partitions at different offsets based on the calculated logging curves of Young's modulus, Poisson's ratio, and fracture density, includes the following steps: using a first formula, calculating the azimuth elastic impedance curves for multiple target reservoir gather partitions at different offsets based on the calculated logging curves of Young's modulus, Poisson's ratio, and fracture density; interpolating each azimuth elastic impedance curve one-to-one in the seismic data volume of each target reservoir gather partition to obtain the initial azimuth elastic impedance model for each target reservoir gather partition. The first formula is as shown in Formula 1 above.

[0091] In the above embodiments, by performing partitioned inversion of elastic parameters and anisotropic parameters based on azimuth elastic impedance for the target reservoir, compared with Embodiment 1, the predicted geostress of the target reservoir is closer to the actual geostress distribution, thus improving the accuracy of geostress prediction.

[0092] For example, in one specific embodiment, step 208, based on the constructed initial models of elastic impedance in each azimuth, inverts the Young's modulus, Poisson's ratio, and fracture density of each target reservoir gather partition under a Bayesian framework, including the following steps: for each initial model of elastic impedance, anisotropic elastic impedance inversion is performed using the sub-azimuth superimposed gathers of the target reservoir gather partition corresponding to the initial model of elastic impedance and the initial model of elastic impedance in each azimuth to obtain an azimuth elastic impedance data volume, and then, based on the azimuth elastic impedance data volume, the Young's modulus, Poisson's ratio, and fracture density are extracted under a Bayesian framework.

[0093] Example 3

[0094] In one embodiment of this application, the horizontal stress vector prediction method differs from Embodiment 1 and Embodiment 2 in that it further includes the following steps after step 210:

[0095] Step 212: Perform least-squares ellipse fitting on the azimuth stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

[0096] Figure 3 schematically illustrates a flowchart of a horizontal stress vector prediction method according to the above embodiments.

[0097] In this application, least squares ellipse fitting refers to fitting an ellipse using the least squares method.

[0098] For example, in least-squares ellipse fitting, the five elements of the ellipse are obtained using the following equations: the x-coordinate of the center point x0, the y-coordinate of the center point y0, the major axis p, the minor axis q, and the ellipse orientation θ, as shown below:

[0099] In the above formula, F(a,x)=Ax 2 +Bxy+Cy 2 +Dx+Ey+F=0, where A, B, C', D, and E are different fitting coefficients. If the coordinates of N points are known, the coefficients of the above quadratic curve equation can be obtained through the least squares fitting method, and thus the directional distribution of the horizontal stress can be obtained.

[0100] For example, the construction process of azimuth-based stacked gathers may include the following steps: Offset Vector Tile (OVT) processing is performed on the raw shot records of the target reservoir acquired in the field to obtain pre-stack OVT gathers; azimuth-based stacking is performed on the pre-stack OVT gathers to obtain azimuth-based stacked gathers; based on the velocity field of the target reservoir, the pre-stack common reflection point (CRP) gathers at the corresponding azimuth angles are converted into common reflection angle gathers; and post-stack time migration data is obtained through preprocessing. The preprocessing includes static correction, spherical compensation, noise suppression, deconvolution, and migration imaging.

[0101] As an improvement to the above embodiments, one embodiment of this application proposes an improved least-squares ellipse fitting algorithm. Specifically, when performing least-squares ellipse fitting on the azimuth superimposed gathers, the following operations are performed: removing singular points present in the fitting; constraining the solution space during ellipse fitting using prior information on the geostress direction to determine the major axis direction of the ellipse; and calculating the ellipse fitting error using the least squares median criterion.

[0102] In this application, the prior information on the geostress direction can be the known point geostress direction, or it can be obtained through the analysis of separated shear waves, imaging logging, and core data.

[0103] In the above embodiments, removing singular points can improve the accuracy of ellipse fitting. To address the difficulty in determining the major and minor axes reflecting the direction of geostress during the fitting process, prior information on the geostress direction is introduced to constrain the solution space during ellipse fitting, thereby accurately determining the direction of the major axis. The least squares median criterion is adopted, requiring the median of the squared residuals to be minimized, excluding points with large deviations, which improves the fitting effect. Combining these improvements enhances the stability and accuracy of least-squares ellipse fitting, thus improving the accuracy of geostress direction prediction in shale reservoirs.

[0104] Referring to Figures 4 to 21, the effectiveness of the above embodiment is illustrated by taking a three-dimensional shale gas reservoir seismic work area as an example. Based on pre-stack seismic data, the seismic work area is divided into zones for geostress vector prediction. The zones for each target reservoir gather are near-field seismic zone, mid-field seismic zone, and far-field seismic zone. The specific process includes the following steps B1 to B10.

[0105] Step B1, Pre-stack Seismic Data Processing and Conversion: Offset Vector Tile (OVT) processing is performed on the raw shot records of the shale reservoir acquired in the field to obtain pre-stack OVT gathers. The OVT gathers are then stacked by azimuth angle, and based on the velocity field, the pre-stack common reflection point (CRP) gathers at the corresponding azimuth angles are converted into common reflection angle gathers. Simultaneously, post-stack time migration data is obtained through static correction, spherical compensation, noise suppression, deconvolution, and migration imaging. The post-stack time migration data volume is shown in Figure 4. In Figure 4, the interpretation horizon indicated by the black line between 1600ms and 2100ms represents the actual bottom interpretation result of the shale formation. The shale reservoir shown in the figure exhibits relatively strong peak reflections, stable reflection characteristics, and relatively high signal-to-noise ratio and resolution, indicating that the processed seismic data is of high quality and can provide reliable seismic data for subsequent seismic interpretation. Figures 6 to 8 show the superimposed profiles obtained by superimposing near-track, middle-track, and far-track seismic traces, respectively. The near-track is a seismic trace with an offset greater than or equal to the first value and less than the second value. The middle-track is a seismic trace with an offset greater than or equal to the second value and less than the third value. The far-track is a seismic trace with an offset greater than or equal to the third value. The first value is less than the second value, and the second value is less than the third value. As can be seen from the figures, the superimposed profiles of the near-track, middle-track, and far-track seismic traces have similar structural features, but there are certain differences in amplitude energy.

[0106] Step B2, pre-stack seismic gather optimization: Angle gathers often contain random noise and multiples, and exhibit near-to-far channel energy differences, resulting in low longitudinal resolution. Therefore, we use singular value decomposition to remove random noise, pull transform to remove multiples, spectral balancing to equalize near-to-far channel energy, and compressed sensing to improve gather resolution. The final result is pre-stack seismic data that meets the requirements for geostress vector seismic prediction, as shown in Figure 5. In Figure 5, the phase axis at 1940 ms represents the main waveform characteristic of the target layer. The phase axis of the target layer exhibits a significant AVO effect, indicating better removal of random noise and relatively high resolution. This significantly improves the quality of the gather data, demonstrating the rationality of the seismic data optimization process.

[0107] Step B3, Well-Seismic Calibration: Impedance curves are calculated using acoustic and density logging curves, and reflection coefficient curves are constructed. These curves are then convolved with seismic wavelets to obtain a well-side synthetic seismic gather. The synthetic seismic record is then subjected to drift and stretching processing. Through multiple iterations, a high correlation coefficient is achieved between the synthetic seismic gather and the well-side seismic data, with similar waveforms and wave group characteristics. The time-depth correspondence obtained at this point represents the optimal time-depth relationship. Figures 9 and 10 show the well-seismic calibration results for wells A3 and A5 in this seismic area, respectively. From left to right, the figures show velocity curves, density curves, reflection coefficient curves, logging stratification, synthetic seismic records, and well-side seismic data. From shallow to deep layers, the waveform characteristics of the synthetic seismic record and the well-side seismic data show good similarity and a high correlation coefficient, indicating good well-seismic calibration results.

[0108] Step B4, Stratigraphic Interpretation: Using time depth as the vertical standard and seismic phase axis as the lateral constraint, considering the multi-layered rationality of shale reservoirs and combining existing geological knowledge, stratigraphic interpretation of the top and bottom of shale reservoirs is carried out in three-dimensional space, providing an interpretative basis for the construction of low-frequency inversion models and pre-stack seismic inversion. Figure 11 shows the stratigraphic interpretation of the target shale reservoir layers carried out along the seismic phase axis from the well point, based on well-seismic calibration. According to the stratigraphic interpretation results, the overall structure of the study area shown in the figure exhibits a trend of high in the north and low in the south, with east-west zonation and north-south block division. The uplift structures are relatively narrow and steep, mainly fault anticlines, while the negative structures are characterized by several relatively wide and gentle synclines.

[0109] Step B5: Calculation of shear wave velocity, p-wave velocity, density, and anisotropy parameters: Well logging curves for sonic transit time, density, TOC content, porosity, and saturation are obtained through well logging interpretation. An anisotropic petrophysical model of the shale reservoir is constructed, focusing on TOC content, organic porosity, horizontal bedding fractures, and reservoir anisotropy. Quartz, carbonate rocks, and TOC content are calculated using the VRH model. Simultaneously, the modulus of the rock matrix is ​​obtained by mixing clay into other mineral contents using Backus averaging. Organic porosity is added using the KT model, and inorganic porosity is added using the DEM model to obtain the modulus of dry rock. Approximate horizontal bedding fractures are introduced using the Pade approximation equation, and approximate vertical structural fractures are added using the Hudson model to obtain the equivalent modulus of the fractured dry rock skeleton. Anisotropic fluid substitution is performed using the Brown-Korringa model to obtain the equivalent modulus of saturated rock. Finally, shear wave velocity, p-wave velocity, density, and anisotropy parameters are calculated using elasticity theory and anisotropy theory.

[0110] Step B6, Establishment of the initial model of near, middle and far azimuth elastic impedance: Based on the time depth established in step B3, the layer obtained in step B4 is interpreted as a lateral constraint. Using the P-wave velocity, S-wave velocity, density and anisotropy parameters obtained in step B5, the logging curves of Young's modulus, Poisson's ratio and fracture density are calculated respectively. Then, the azimuth elastic impedance curves of the near, middle and far gathers are calculated according to the equation shown in Formula 1 above. Lateral interpolation is performed on each curve in the three-dimensional data volume space to finally obtain the initial model of near, middle and far azimuth elastic impedance. Figures 12 and 13 show the calculation results of the azimuth elastic impedance curves for the near, middle, and far channels. Figure 12 shows the calculation results of the azimuth elastic impedance curve for well A3, and Figure 13 shows the calculation results of the azimuth elastic impedance curve for well A7. In both figures, from left to right, the curves are P-wave velocity, S-wave velocity, density, near channel azimuth elastic impedance curve, middle channel azimuth elastic impedance curve, far channel azimuth elastic impedance curve, and logging stratification. As can be seen from the figures, the azimuth elastic impedances of the near, middle, and far channels have similar variation characteristics, but there are certain impedance differences.

[0111] Step B7, Seismic parameter inversion based on anisotropic elastic impedance inversion: Using the low-frequency models of near-, mid-, and far-directional elastic impedance obtained in step A6 as initial models, pre-stack inversion equations for near-, mid-, and far-directional elastic impedance under a Bayesian framework are established respectively, yielding inversion results for different angles. The logarithm of both sides of Equation 1 is taken for linearization, expressing the near-, mid-, and far-directional elastic impedance as the product of the Lamé parameters of the coefficient matrix and the density vector (see Equation 2). The inverse problem of Equation 2 is solved, finally obtaining the inversion results for Young's modulus, Poisson's ratio, and fracture density. Figure 14 shows the Poisson's ratio inversion profile of well A8 obtained through anisotropic elastic impedance inversion, and Figure 15 shows the fracture density inversion profile of well A8 obtained through anisotropic elastic impedance inversion. The profile results show that the shale reservoir section has a low Poisson's ratio and a high fracture density. Figure 16 shows the Poisson's ratio inversion planar diagram obtained by anisotropic elastic impedance inversion, and Figure 17 shows the fracture density inversion planar diagram obtained by anisotropic elastic impedance inversion. As can be seen from the figures, the Poisson's ratio is low in the northern and southeastern regions of wells A3 and A7, and fractures are more developed in the northern and southeastern regions of wells A3 and A7.

[0112] Step B8, Vertical Principal Stress Calculation: Combining the velocity field obtained from seismic data processing, the vertical principal stress is calculated by integration using the density data volume obtained from pre-stack inversion. The integration process is as follows:

[0113] In Formula 4, ρ is the density of the stratum, g is the gravitational acceleration, and H is the stratum depth.

[0114] Step B9, Maximum and Minimum Horizontal Stress Prediction: Using the fracture density and Poisson's ratio obtained from the inversion in Step B7 and the vertical principal stress calculated in Step B8, the magnitudes of the maximum and minimum horizontal stresses of the reservoir in this work area are calculated according to Formula 3 above. Figure 18 shows the maximum horizontal stress prediction profile through Well A8, and Figure 19 shows the minimum horizontal stress prediction profile through Well A8. As can be seen from the figures, the profiles show that the maximum and minimum horizontal stress prediction results gradually increase from top to bottom.

[0115] Step B10, Prediction of the directions of maximum and minimum horizontal stress: Based on the improved least-squares ellipse fitting algorithm in the above embodiments, the direction distribution of maximum and minimum horizontal stress is obtained. The direction of maximum horizontal stress is the major axis of the ellipse, and the direction of minimum horizontal stress is the minor axis of the ellipse. Figure 20 shows the predicted planar diagram of the maximum horizontal stress vector in the seismic area, and Figure 21 shows the predicted planar diagram of the minimum horizontal stress vector in the seismic area. The colored fills in the figures represent the magnitude of the horizontal stress, and the black lines represent the direction of the horizontal stress. As can be seen from the figures: Regarding the distribution of principal stress magnitudes, the magnitudes of the maximum and minimum horizontal stresses have similar characteristics, with higher horizontal stresses in the northwest, south-central, and northeast parts of the seismic area; Regarding the distribution of principal stress directions, the main distribution directions of the maximum horizontal stress are near east-west, northwest-southeast, and southwest-northeast, while the main distribution directions of the minimum horizontal stress are near north-south and northwest-southeast. The prediction results have a high degree of consistency with known geological knowledge, and the prediction results confirm the effectiveness and applicability of the embodiments of this application.

[0116] Example 4

[0117] This application embodiment also provides a shale reservoir geostress seismic prediction device, comprising: a first calculation module, used to calculate the elastic parameters and anisotropic parameters of the target reservoir using a constructed or obtained anisotropic rock physics model of the target reservoir; a second calculation module, used to calculate the logging curves of Young's modulus, Poisson's ratio, and fracture density using the calculated elastic parameters and anisotropic parameters, with the layer interpretation of the target reservoir obtained through seismic data volume interpretation as a constraint; an inversion model construction module, used to construct an azimuth elastic impedance initial model based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves; and an inversion module, used to construct an azimuth initial model based on the constructed azimuth initial model. Based on the initial model of elastic impedance, Young's modulus, Poisson's ratio, and fracture density are obtained by inversion within a Bayesian framework. The geostress magnitude calculation module is used to substitute the inverted fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir into the prediction formula to calculate the maximum and minimum horizontal stresses of the target reservoir. The prediction formula uses only the vertical principal stress, fracture density, and Poisson's ratio as independent variables and the maximum and minimum horizontal stresses as dependent variables. The geostress direction prediction module is used to perform least-squares ellipse fitting on the azimuthally stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

[0118] In one specific embodiment, the step of constructing an initial azimuth elastic impedance model based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves includes: calculating the azimuth elastic impedance curve using a first formula based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves; and interpolating the azimuth elastic impedance curves in the seismic data volume to obtain the initial azimuth elastic impedance model; the first formula is:

[0119] Wherein, θ represents the incident angle; φ represents the azimuth angle; E, υ, and e represent Young's modulus, Poisson's ratio, and fracture density, respectively; E0 represents the mean of Young's modulus; υ0 represents the mean of Poisson's ratio; a represents a coefficient related to Young's modulus, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; b represents a coefficient related to Poisson's ratio, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; c represents a coefficient related to fracture density, and is expressed by the incident angle, the azimuth angle, and the square of the P-wave velocity ratio.

[0120] In one specific embodiment, based on the constructed initial azimuth elastic impedance model, Young's modulus, Poisson's ratio, and fracture density are obtained by inversion within a Bayesian framework. This includes: performing anisotropic elastic impedance inversion using the target reservoir's azimuth stacked gathers and the initial azimuth elastic impedance model to obtain azimuth elastic impedance data volume; and extracting Young's modulus, Poisson's ratio, and fracture density within a Bayesian framework based on the azimuth elastic impedance data volume.

[0121] In one specific embodiment, based on the azimuth elastic impedance data volume, Young's modulus, Poisson's ratio, and crack density are extracted under a Bayesian framework, including: performing logarithmic processing on the first formula to obtain a linearized azimuth elastic impedance equation; using the matrix form of the linearized azimuth elastic impedance equation as the parameter extraction equation, and extracting Young's modulus, Poisson's ratio, and crack density under a Bayesian framework.

[0122] The parameter extraction equation is as follows:

[0123] Where θ1, θ2, and θ3 represent different incident angles, φ1, φ2, and φ3 represent different azimuth angles, EI represents elastic impedance, and EI0 represents the mean value of elastic impedance.

[0124] In one specific embodiment, the prediction formula is:

[0125] Where, σ H σ represents the maximum horizontal stress. VLet υ represent the vertical principal stress, ν represent Poisson's ratio, e represent the crack density, C be a constant, and σ represent the vertical principal stress. h This represents the minimum horizontal stress.

[0126] In one specific embodiment, the elastic parameters include transverse wave velocity, longitudinal wave velocity, and density.

[0127] In one specific embodiment, when performing least-squares elliptic fitting on the azimuth-divided superimposed gathers, the following operations are performed: removing singular points present in the fitting; constraining the solution space of the elliptic fitting using prior information on the geostress direction; and calculating the elliptic fitting error using the least squares median criterion.

[0128] In one specific embodiment, there are multiple azimuth elastic impedance initial models, and each azimuth elastic impedance initial model is an azimuth elastic impedance initial model for each target reservoir gather partition with different offset distances.

[0129] In one specific embodiment, constructing an anisotropic rock physics model of the target reservoir includes: interpreting the logging data of the target reservoir to obtain logging curves of sonic transit time, density, porosity, and saturation; and constructing an anisotropic rock physics model of the target reservoir using the logging curves of sonic transit time, density, porosity, and saturation.

[0130] In one specific embodiment, the device further includes: a preprocessing module for preprocessing the seismic data volume of the target reservoir and jumping to trigger the execution of inversion of Young's modulus, Poisson's ratio and fracture density based on the initial azimuth elastic impedance model within a Bayesian framework.

[0131] In one specific embodiment, the seismic data volume is a pre-stack seismic data volume, and the preprocessing includes: filtering out noise in the pre-stack seismic data volume using singular value decomposition; removing multiples in the pre-stack seismic data volume using pull transform; equalizing the amplitude energy of different seismic traces in the pre-stack seismic data volume using a spectral balance algorithm; and enhancing the gather resolution in the pre-stack seismic data volume using a compressed sensing algorithm.

[0132] Figure 22 schematically illustrates a structural block diagram of a computer device according to an embodiment of this application. In one embodiment, a computer device is provided, which may be a terminal, and its internal structure diagram may be as shown in Figure 22. The computer device includes a processor A01, a network interface A02, a display screen A04, an input device A05, and a memory (not shown) connected via a system bus. The processor A01 of the computer device provides computing and control capabilities. The memory of the computer device includes internal memory A03 and a non-volatile storage medium A06. The non-volatile storage medium A06 stores an operating system B01 and a computer program B02. The internal memory A03 provides an environment for the operation of the operating system B01 and the computer program B02 in the non-volatile storage medium A06. The network interface A02 of the computer device is used for communication with an external terminal via a network connection. When the computer program is executed by the processor A01, it implements a horizontal stress vector prediction method. The display screen A04 of the computer device can be an LCD screen or an e-ink screen. The input device A05 of the computer device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the casing of the computer device, or an external keyboard, touchpad, or mouse, etc.

[0133] Those skilled in the art will understand that the structure shown in Figure 22 is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0134] In one embodiment, the shale reservoir in-situ stress seismic prediction device provided in this application can be implemented as a computer program, which can run on a computer device as shown in FIG22. The memory of the computer device can store the various program modules that constitute the shale reservoir in-situ stress seismic prediction device. The computer program, composed of the various program modules, causes the processor to execute the steps in the horizontal stress vector prediction method of the various embodiments of this application described in this specification.

[0135] In one embodiment, this application also provides a machine-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the horizontal stress vector prediction method described in the above embodiments.

[0136] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0137] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0138] It should also be noted that the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0139] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A method for predicting horizontal stress vectors, characterized in that, include: The elastic parameters and anisotropic parameters of the target reservoir are calculated using an anisotropic rock physics model of the target reservoir. Using the layer interpretation of the target reservoir obtained through seismic data volume interpretation as a constraint, the well logging curves of Young's modulus, Poisson's ratio and fracture density are calculated using the calculated elastic parameters and anisotropic parameters. An initial model of azimuth elastic impedance is constructed based on the logging curves obtained from the calculation of Young's modulus, Poisson's ratio, and fracture density. Based on the initial model of azimuth elastic impedance, Young's modulus, Poisson's ratio and crack density are obtained by inversion within a Bayesian framework; Substituting the fracture density, Poisson's ratio obtained from the inversion and the calculated vertical principal stress of the target reservoir into the prediction formula, the maximum and minimum horizontal stress of the target reservoir are calculated. The prediction formula uses only the vertical principal stress, fracture density and Poisson's ratio as independent variables and the maximum and minimum horizontal stress as dependent variables. Least-square ellipse fitting is performed on the azimuth stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

2. The method according to claim 1, characterized in that, The initial model of azimuth elastic impedance is constructed based on the logging curves obtained from the calculated Young's modulus, Poisson's ratio, and fracture density, including: Using the first formula, the azimuth elastic impedance curve is calculated based on the logging curves obtained from Young's modulus, Poisson's ratio, and fracture density. The azimuth elastic impedance curve is used to interpolate in the seismic data volume to obtain the initial azimuth elastic impedance model. The first formula is: Wherein, θ represents the incident angle; φ represents the azimuth angle; E, υ, and e represent Young's modulus, Poisson's ratio, and fracture density, respectively; E0 represents the mean of Young's modulus; υ0 represents the mean of Poisson's ratio; a represents a coefficient related to Young's modulus, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; b represents a coefficient related to Poisson's ratio, and is expressed by the incident angle, the power of the target reservoir density with respect to the P-wave velocity, and the square of the P-wave velocity ratio; c represents a coefficient related to fracture density, and is expressed by the incident angle, the azimuth angle, and the square of the P-wave velocity ratio.

3. The method according to claim 2, characterized in that, Based on the initial model of azimuth elastic impedance, the Young's modulus, Poisson's ratio, and crack density are obtained by inversion within a Bayesian framework, including: Anisotropic elastic impedance inversion is performed using the constructed target reservoir azimuth stacked gathers and the initial azimuth elastic impedance model to obtain the azimuth elastic impedance data volume. Based on the aforementioned azimuth elastic impedance data, Young's modulus, Poisson's ratio, and crack density were extracted within a Bayesian framework.

4. The method according to claim 3, characterized in that, Based on the aforementioned azimuth elastic impedance data, Young's modulus, Poisson's ratio, and crack density are extracted within a Bayesian framework, including: Logarithmizing the first formula yields the linearized azimuth elastic impedance equation. The matrix form of the linearized azimuthal elastic impedance equation is used as the parameter extraction equation, and Young's modulus, Poisson's ratio and crack density are extracted under the Bayesian framework. The parameter extraction equation is as follows: Where θ1, θ2, and θ3 represent different incident angles, φ1, φ2, and φ3 represent different azimuth angles, EI represents elastic impedance, and EI0 represents the mean value of elastic impedance.

5. The method according to claim 1, characterized in that, The prediction formula is: Where, σ H σ represents the maximum horizontal stress. V Let υ represent the vertical principal stress, ν represent Poisson's ratio, e represent the crack density, C be a constant, and σ represent the vertical principal stress. h This represents the minimum horizontal stress.

6. The method according to claim 1, characterized in that, When performing least-squares elliptic fitting on the aforementioned azimuth superimposed gathers, the following operations are performed: Remove singularities present in the fit; The solution space for ellipse fitting is constrained by prior information on the direction of geostress. The ellipse fitting error is calculated using the least squares median criterion.

7. The method according to claim 1, characterized in that, Least-squares ellipse fitting was performed on the azimuth-based stacked gathers of the target reservoir to obtain the directional distributions of the maximum and minimum horizontal stresses within the target reservoir, including: Least square ellipse fitting is performed on the azimuth stacked gathers of the target reservoir to determine the direction of the major axis and the direction of the minor axis of the ellipse. The direction of the major axis of the ellipse is determined as the direction of the maximum horizontal stress within the target reservoir. The direction of the minor axis of the ellipse is determined as the direction of the minimum horizontal stress within the target reservoir.

8. The method according to claim 1, characterized in that, The number of initial azimuth elastic impedance models is multiple, and each initial azimuth elastic impedance model is an initial azimuth elastic impedance model for each target reservoir gather partition with different offset distances.

9. The method according to claim 1, characterized in that, Construct an anisotropic rock physics model of the target reservoir, including: The logging data of the target reservoir is interpreted to obtain logging curves for sonic transit time, density, porosity and saturation; An anisotropic rock physics model of the target reservoir is constructed using well logging curves based on sonic transit time, density, porosity, and saturation.

10. The method according to claim 1, characterized in that, Also includes: The seismic data volume of the target reservoir is preprocessed and then jumps to trigger the execution of Young's modulus, Poisson's ratio and fracture density inversion under a Bayesian framework based on the initial azimuth elastic impedance model.

11. The method according to claim 10, characterized in that, The seismic data volume is a pre-stack seismic data volume, and the preprocessing includes: Singular value decomposition is used to filter out noise in pre-stack seismic data volumes; Multiples in pre-stack seismic data volumes are removed using pull transform; The amplitude energy of different seismic traces in the pre-stack seismic data volume is equalized using the spectral balance algorithm; Compressed sensing algorithms are used to enhance the resolution of gathers in pre-stack seismic data volumes.

12. A horizontal stress vector prediction device, characterized in that, include: The first calculation module is used to calculate the elastic parameters and anisotropic parameters of the target reservoir using the constructed or obtained anisotropic rock physics model of the target reservoir; The second calculation module is used to use the layer interpretation of the target reservoir obtained through the interpretation of seismic data volume as a constraint, and to calculate the logging curves of Young's modulus, Poisson's ratio and fracture density using the calculated elastic parameters and anisotropic parameters. The inversion model construction module is used to construct an initial model of azimuth elastic impedance based on the calculated Young's modulus, Poisson's ratio, and fracture density logging curves. The inversion module is used to invert Young's modulus, Poisson's ratio, and crack density based on the initial model of azimuthal elastic impedance within a Bayesian framework. The geostress magnitude calculation module is used to substitute the fracture density, Poisson's ratio, and calculated vertical principal stress of the target reservoir obtained by inversion into the prediction formula to calculate the maximum and minimum horizontal stress of the target reservoir. The prediction formula only uses the vertical principal stress, fracture density, and Poisson's ratio as independent variables and the maximum and minimum horizontal stress as dependent variables. The geostress direction prediction module is used to perform least-squares ellipse fitting on the azimuth stacked gathers of the target reservoir to obtain the directional distribution of the maximum and minimum horizontal stresses within the target reservoir.

13. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the horizontal stress vector prediction method according to any one of claims 1-12.

14. A machine-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the horizontal stress vector prediction method according to any one of claims 1-11.

Images