Method and system for predicting maximum horizontal principal stress of shale gas reservoirs

Abstract

Embodiments of the present invention relate to the technical field of shale gas development, and provide a method and system for predicting the maximum horizontal principal stress of shale gas reservoirs. The method comprises: performing cross-arranged gather extraction on acquired original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; acquiring well logging data, performing well-to-seismic calibration and seismic horizon interpretation on the basis of the well logging data, and obtaining an interpreted horizon of the current reservoir that meets in-situ stress prediction as a target horizon; constructing a low-frequency model on the basis of the well logging data and the target horizon, and calculating a Poisson's ratio and a fracture density on the basis of the low-frequency model; calculating a vertical principal stress of the target horizon on the basis of a velocity field, a formation density and a density inversion body; and calculating the maximum horizontal principal stress of the current reservoir on the basis of the Poisson's ratio, the fracture density, the vertical principal stress and a preset principal stress solving model. The solution of the present invention solves the problem of low accuracy in existing solutions for predicting the maximum horizontal principal stress of shale gas reservoirs.

Description

Method and system for predicting maximum horizontal principal stress in shale gas reservoirs

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims the benefit of Chinese patent application 202311713818.3 filed on December 13, 2023, the contents of which are incorporated herein by reference. Technical Field

[0003] The present invention relates to the technical field of shale gas development, and in particular to a method for predicting the maximum horizontal principal stress of a shale gas reservoir and a system for predicting the maximum horizontal principal stress of a shale gas reservoir. Background Art

[0004] Maximum horizontal principal stress plays a crucial role in shale gas exploration and development in shale reservoirs. In geological engineering, maximum horizontal principal stress refers to the maximum horizontal stress in a formation. It is the maximum stable stress in the formation. In shale oil and gas exploration, the magnitude of maximum horizontal principal stress determines fractures and rock damage in the formation, as well as the migration and storage of oil and gas in the formation. In shale oil and gas development, maximum horizontal principal stress determines the location of wells, the direction of well trajectories, and the parameters of fracturing design.

[0005] In recent years, extensive research has been conducted on the prediction of maximum horizontal principal stress, using seismic and well logging data to quantitatively predict geostress. For example, Liu Haojuan (2021) established a geostress prediction model. Based on detailed 3D seismic interpretation and 3D seismic prestack inversion, they used wellpoint data to simulate and select regional adaptive parameters for stress calculation. They then conducted 3D simulations of the geostress field and predicted the maximum and minimum horizontal principal stress directions, as well as the horizontal stress difference coefficient. Qu Yang (2019) applied resistivity imaging logging and shear wave anisotropy logging to interpret the direction of the maximum horizontal principal stress in the region. They calculated the minimum principal stress in the study area using the Newberry model and the maximum principal stress based on dual-caliper analysis. This helped clarify the trajectory direction of horizontal well deployment in the region, providing a basis for reservoir fracturing. Qi Dunke (2017) proposed a theory and method for extracting more reliable stress field directions by analyzing microseismic monitoring data from horizontal wells. Based on microseismic monitoring data from three horizontal wells in the Daqing Oilfield, they derived the stress field direction in the monitoring well area. Li Guangquan (2012) studied the method of obtaining the maximum horizontal principal stress using wellbore collapse information, and proposed a method of inverting the maximum horizontal principal stress using numerical simulation methods combined with logging data analysis, which can truly reflect the stress environment of deep formations and the mechanical properties of formation rocks.

[0006] However, in actual application, these schemes all have certain accuracy issues, and the calculation results of maximum horizontal principal stress cannot be directly used in industrial production. Based on this, it is necessary to create a shale gas reservoir maximum horizontal principal stress prediction scheme with guaranteed accuracy. Summary of the Invention

[0007] The purpose of the embodiments of the present invention is to provide a method and system for predicting the maximum horizontal principal stress of a shale gas reservoir, so as to at least solve the problem of low accuracy of existing maximum horizontal principal stress prediction schemes for shale gas reservoirs.

[0008] In order to achieve the above-mentioned objectives, the first aspect of the present invention provides a method for predicting the maximum horizontal principal stress of a shale gas reservoir, the method comprising: performing a cross-arranged extraction of the collected original shot gathers to obtain a corresponding OVT gather, and preprocessing the OVT gathers; collecting well logging data, and performing well-seismic calibration and seismic layer interpretation based on the well logging data to obtain an interpreted layer of the current reservoir that meets the ground stress prediction as a target layer; constructing a low-frequency model based on the well logging data and the target layer, and calculating the Poisson's ratio and fracture density based on the low-frequency model; calculating the vertical principal stress of the target layer based on the velocity field, formation density and density inversion volume; and calculating the maximum horizontal principal stress of the current reservoir based on the Poisson's ratio, the fracture density, the vertical principal stress and a preset principal stress obtaining model.

[0009] Optionally, the cross-arranged channel extraction process is performed on the collected original shot gathers to obtain corresponding OVT channel gathers, including: cross-arranged channel extraction process is performed on the collected original shot gathers to obtain seismic channel gathers at the same detection point under the same shot line; each cross-arranged channel gather is divided into OVT units, and each cross-arranged channel gather after division is subjected to channel extraction process to obtain OVT channel gathers.

[0010] Optionally, the preprocessing of the OVT gather includes: performing five-dimensional regularization processing on the OVT gather to obtain a first OVT gather; performing offset processing on the first OVT gather to obtain a second OVT gather; performing anisotropy correction on the second OVT gather to obtain a third OVT gather; picking up the velocity field on the third OVT gather to convert the third OVT gather into an incident angle gather containing azimuth information; performing azimuth angle superposition processing on incident angle gathers in different azimuth segments to obtain azimuth gathers; performing one or more of Radon transform, wavelet threshold method and spectral decomposition on the azimuth gather to complete the OVT gather preprocessing.

[0011] Optionally, after collecting the logging data, the method further includes: performing preprocessing of the logging curves in the logging data, including: standardizing each logging curve and identifying the logging curves that lack shear wave curves; wherein the various logging curves include: acoustic wave curves, density curves, natural potential curves, gamma curves and resistivity curves; for the logging curves that lack shear wave curves, constructing a model relationship between the corresponding logging shear wave velocity and preset sensitive parameters, and obtaining the shear wave curve based on the neural network and the model relationship; obtaining the Young's modulus curve and the Poisson's ratio curve based on the shear wave curves of each logging curve after standardization; and constructing an anisotropic theoretical model based on the Young's modulus curve and the Poisson's ratio curve to calculate the fracture density curve.

[0012] Optionally, the well seismic calibration and seismic layer interpretation based on the well logging data are performed to obtain the interpreted layer of the current reservoir that meets the ground stress prediction as the target layer, including: dynamically extracting seismic wavelets based on the well logging data, and continuously synthesizing seismic records based on the seismic wavelets, and taking the seismic wavelets corresponding to the seismic records whose correlation coefficient with the actual seismic records of the well bypass is less than a preset correlation coefficient threshold as the calibration result; based on the calibration result, structural interpretation of the current seismic layer is performed based on a preset knowledge graph to obtain the interpreted layer of the current reservoir that meets the ground stress prediction as the target layer.

[0013] Optionally, constructing a low-frequency model based on the logging data and the target layer includes: using the Young's modulus curve, the Poisson's ratio curve and the fracture density curve as vertical source data, and the structural interpretation of the target layer as a horizontal constraint, performing extrapolation and interpolation to obtain a low-frequency model data body of Young's modulus, Poisson's ratio data body and fracture density in three-dimensional space as a low-frequency model.

[0014] Optionally, the calculation of Poisson's ratio and fracture density based on the low-frequency model includes: calculating the seismic reflection coefficient at each incident angle and azimuth angle based on the low-frequency model to obtain a reflection coefficient matrix; constructing a pre-stack inversion model based on the seismic wavelet matrix and the reflection coefficient matrix to obtain a seismic reflection amplitude matrix; constructing an objective function based on a Bayesian framework and the seismic reflection amplitude matrix, and calculating the Poisson's ratio and fracture density based on the objective function.

[0015] Optionally, the seismic reflection amplitude matrix construction rule is: d NMK×1 =G NMK×3K m 3K×1 .

[0016] Among them, d NMK×1 is the seismic reflection amplitude matrix; G NMK×3K is the seismic wavelet matrix; m 3K×1 is the reflection coefficient matrix.

[0017] Optionally, the objective function is expressed as:

[0018] Among them, σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well bypass; d is the seismic reflection amplitude matrix; G is the seismic wavelet matrix; T is the sampling period; Q is a diagonal matrix, which can be expressed as:

[0019] Among them, are the variance of Young's modulus, the variance of Poisson's ratio and the variance of crack density parameter, respectively.

[0020] Optionally, calculating the vertical principal stress of the target layer based on the well logging data includes: reading a signal velocity field and formation density based on the well logging data; and calculating the vertical principal stress based on the velocity field, the formation density, and a preset density inversion volume, wherein the calculation rule is:

[0021] Among them, σ V is the vertical principal stress; ρ is the formation density; g w is the acceleration due to gravity; v f is the velocity field; t is the round-trip travel time.

[0022] Optionally, the preset principal stress obtaining model is:

[0023] Among them, σ H is the maximum principal stress; ν is Poisson's ratio; g is the preset variable; and e is the crack density.

[0024] The second aspect of the present invention provides a shale gas reservoir maximum horizontal principal stress prediction system, the system comprising: an acquisition unit, for performing a cross-arranged extraction of the acquired original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; a processing unit, for acquiring well logging data, and performing well-seismic calibration and seismic layer interpretation based on the well logging data, to obtain an interpreted layer of the current reservoir that meets the ground stress prediction as a target layer; a training unit, for: constructing a low-frequency model based on the well logging data and the target layer, and calculating the Poisson's ratio and fracture density based on the low-frequency model; calculating the vertical principal stress of the target layer based on the velocity field, formation density and density inversion body; and an output unit, for obtaining the maximum horizontal principal stress of the current reservoir by calculating a model based on the Poisson's ratio, the fracture density, the vertical principal stress and a preset principal stress.

[0025] On the other hand, the present invention provides a computer-readable storage medium having instructions stored thereon, which, when executed on a computer, enables the computer to execute the above-mentioned method for predicting the maximum horizontal principal stress of a shale gas reservoir.

[0026] Through the above technical solution, the present invention first derives a maximum horizontal principal stress prediction formula based on fracture density and Poisson's ratio. Simultaneously, based on prestack anisotropic media theory, an approximate formula for the longitudinal wave orientation (AVO) is established based on Young's modulus, Poisson's ratio, and fracture density. Using this new approximate equation within the framework of Bayesian theory, an inversion equation is established to achieve prestack inversion of seismic parameters. This method improves the accuracy of maximum horizontal principal stress prediction through direct inversion.

[0027] Other features and advantages of the embodiments of the present invention will be described in detail in the subsequent detailed description. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] The accompanying drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification. Together with the following detailed description, they are used to explain the embodiments of the present invention, but do not constitute a limitation of the embodiments of the present invention. In the accompanying drawings:

[0029] FIG1 is a flowchart of the steps for predicting the maximum horizontal principal stress of a shale gas reservoir provided by one embodiment of the present invention;

[0030] FIG2 is a well seismic calibration result of an experimental well provided by an embodiment of the present invention;

[0031] FIG3 is an OVT domain spiral gather of Inline19325 line after OVT processing provided by one embodiment of the present invention;

[0032] FIG4 is an azimuth angle gather of an OVT spiral gather provided by one embodiment of the present invention after azimuth processing;

[0033] FIG5 is a diagram of a sub-azimuth gather provided by one embodiment of the present invention after random noise removal;

[0034] FIG6 is a gather of sub-azimuth angles provided by one embodiment of the present invention after gather flattening;

[0035] FIG7 is a result of horizon interpretation of a target layer of a shale gas reservoir provided by one embodiment of the present invention;

[0036] FIG8 is a seismic profile of an experimental well provided by one embodiment of the present invention;

[0037] FIG9 is a Poisson's ratio inversion profile of a test well provided by one embodiment of the present invention;

[0038] FIG10 is a fracture density inversion profile of a test well provided by one embodiment of the present invention;

[0039] FIG11 is a plan view of Poisson's ratio inversion of a study area provided by one embodiment of the present invention;

[0040] FIG12 is a plan view of fracture density inversion in a study area provided by one embodiment of the present invention;

[0041] FIG13 is a maximum horizontal principal stress profile result of a test well provided by one embodiment of the present invention;

[0042] FIG14 is a plan view showing the maximum horizontal principal stress magnitude prediction provided by one embodiment of the present invention;

[0043] FIG15 is a system structure diagram of a shale gas reservoir maximum horizontal principal stress prediction system provided by one embodiment of the present invention. DETAILED DESCRIPTION

[0044] The following describes the specific embodiments of the present invention in detail with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only used to illustrate and explain the present invention and are not intended to limit the present invention.

[0045] FIG1 is a flow chart of a method for predicting the maximum horizontal principal stress of a shale gas reservoir provided by an embodiment of the present invention. As shown in FIG1 , an embodiment of the present invention provides a method for predicting the maximum horizontal principal stress of a shale gas reservoir, the method comprising:

[0046] Step S10: performing cross-arrangement extraction on the collected original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers.

[0047] Specifically, the cross-arranged channel extraction process for the collected original shot gathers to obtain corresponding OVT channel gathers includes: cross-arranged channel extraction process for the collected original shot gathers to obtain seismic channel gathers at the same detection point under the same shot line; OVT unit division of each cross-arranged channel gather, and channel extraction processing of each divided cross-arranged channel gather to obtain OVT channel gathers.

[0048] In an embodiment of the present invention, original shot gathers collected in the field are cross-arranged and extracted to obtain seismic gathers of the same shot line and the same receiver point. The cross-arranged gathers are then divided into OVT units. Based on the division, the OVT unit extracted gathers are processed to obtain OVT gathers. The OVT gathers are subjected to five-dimensional regularization, the OVT domain gathers are subjected to omnidirectional migration, and anisotropic correction is performed on the gathers to obtain corrected OVT gathers, thereby providing a data basis for pre-stack fracture and ground stress prediction.

[0049] Preferably, the preprocessing of the OVT gather includes: performing five-dimensional regularization processing on the OVT gather to obtain a first OVT gather; performing offset processing on the first OVT gather to obtain a second OVT gather; performing anisotropy correction on the second OVT gather to obtain a third OVT gather; picking up the velocity field on the third OVT gather to convert the third OVT gather into an incident angle gather containing azimuth information; performing azimuth angle superposition processing on incident angle gathers in different azimuth segments to obtain azimuth gathers; performing one or more of Radon transform, wavelet threshold method and spectral decomposition on the azimuth gather to complete the OVT gather preprocessing.

[0050] Specifically, the velocity field of the OVT gathers obtained above is processed and picked up to convert it into an incident angle gather containing azimuth information. The incident angle gathers of different azimuth segments are superimposed on them in different azimuth angles to obtain incident angle gathers of different azimuths (abbreviated as azimuth gathers). In view of the problems of multiple wave interference, signal-to-noise ratio and low resolution in the azimuth gathers, the quality of the gather data is improved through high-precision Radon transform, wavelet threshold method and spectral decomposition.

[0051] Step S20: collecting well logging data, and performing well seismic calibration and seismic horizon interpretation based on the well logging data to obtain the interpretation horizon of the current reservoir that meets the in-situ stress prediction as the target horizon.

[0052] Preferably, after collecting the logging data, the method further includes: performing preprocessing of the logging curves in the logging data, including: standardizing each logging curve and identifying the logging curves that lack shear wave curves; wherein the various logging curves include: acoustic wave curves, density curves, natural potential curves, gamma curves and resistivity curves; for the logging curves that lack shear wave curves, constructing a model relationship between the corresponding logging shear wave velocity and preset sensitive parameters, and obtaining the shear wave curve based on the neural network and the model relationship; obtaining the Young's modulus curve and the Poisson's ratio curve based on the shear wave curves of each logging curve after standardization; and constructing an anisotropic theoretical model based on the Young's modulus curve and the Poisson's ratio curve to calculate the fracture density curve.

[0053] Specifically, well logging curves within the work area, such as acoustic wave curves, density curves, natural potential curves, gamma curves, and resistivity curves, are standardized. For wells lacking shear wave curves, a model relationship between shear wave velocity and sensitive parameters is constructed. A neural network algorithm is used to calculate the shear wave velocity curve, Young's modulus, and Poisson's ratio curves. An anisotropic theoretical model is then established to calculate the fracture density curve.

[0054] Furthermore, the well seismic calibration and seismic layer interpretation based on the logging data are performed to obtain the interpretation layer of the current reservoir that meets the ground stress prediction as the target layer, including: dynamically extracting seismic wavelets based on the logging data, and continuously synthesizing seismic records based on the seismic wavelets, and taking the seismic wavelets corresponding to the seismic records whose correlation coefficient with the actual seismic records of the well bypass is less than a preset correlation coefficient threshold as the calibration result; based on the calibration result, the structural interpretation of the current seismic layer is performed based on a preset knowledge graph to obtain the interpretation layer of the current reservoir that meets the ground stress prediction as the target layer.

[0055] Specifically, through dynamic seismic wavelet extraction, synthetic seismic record creation and drift, and comparison of synthetic seismic records with well bypasses, the final calibration result is obtained when the correlation coefficient between the two reaches the preset requirement. Based on this calibration result, the target layer is identified as a seismic response characteristic. Combined with existing structural knowledge, the target layer is interpreted to obtain an interpreted layer that meets the requirements for ground stress prediction and prevents the layer from jumping up and down in three-dimensional space.

[0056] Step S30: constructing a low-frequency model based on the well logging data and the target layer, and calculating the Poisson's ratio and fracture density based on the low-frequency model.

[0057] Specifically, the Young's modulus curve, the Poisson's ratio curve and the fracture density curve are used as vertical source data, and the structural interpretation of the target layer is used as a horizontal constraint to perform extrapolation and interpolation to obtain a low-frequency model data volume of Young's modulus, Poisson's ratio and fracture density in three-dimensional space as a low-frequency model.

[0058] In an embodiment of the present invention, under the control of the time-depth curve obtained in step S20, the Young's modulus curve, the Poisson's ratio curve and the fracture density curve are used as the source data in the vertical direction, and under the lateral constraints of the target layer data, extrapolation interpolation is performed to obtain the Young's modulus, Poisson's ratio data volume and the fracture density low-frequency model data volume in three-dimensional space.

[0059] Furthermore, the calculation of Poisson's ratio and fracture density based on the low-frequency model includes: calculating the seismic reflection coefficient at each incident angle and azimuth angle based on the low-frequency model to obtain a reflection coefficient matrix; obtaining a seismic reflection amplitude matrix based on a pre-stack inversion model of a seismic wavelet matrix and the reflection coefficient matrix component; constructing an objective function based on a Bayesian framework and the seismic reflection amplitude matrix, and calculating the Poisson's ratio and fracture density based on the objective function.

[0060] Specifically, the Young's modulus, Poisson's ratio data volume and fracture density low-frequency model data volume are used to calculate the seismic reflection coefficient. The calculation rule is:

[0061] Where θ is the incident angle, φ is the azimuth angle, a is the power exponent of density with respect to the P-wave velocity, and g is the square of the ratio of the shear to the P-wave velocity. Here a and g are constants, and E, ν, and e represent Young's modulus, Poisson's ratio, and crack density, respectively.

[0062] Furthermore, for the seismic reflection amplitudes at different incident angles and azimuths, the seismic reflection coefficients are obtained by multiplying the coefficient matrix with the vector of the parameters to be inverted, and the seismic reflection amplitudes multiplied by the wavelet matrix are used to construct the prestack inversion equation: NMK×1 =G NMK×3K m 3K×1 .

[0063] Among them, d NMK×1 is the seismic reflection amplitude matrix; G NMK×3K is the seismic wavelet matrix; m 3K×1 is the reflection coefficient matrix.

[0064] Where N is the number of incident angles, M is the number of azimuth angles, K is the number of sampling points, and wvlt is the sub-bo matrix. S(θ i ,φ j ) is the incident angle θ i and azimuth is φ j seismic reflection amplitude.

[0065] Further:

[0066] Based on the above rules, in the Bayesian framework, the objective function is constructed to obtain the final inversion parameter estimation:

[0067] Among them, σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well bypass; d is the seismic reflection amplitude matrix; G is the seismic wavelet matrix; T is the sampling period; Q is a diagonal matrix, which can be expressed as:

[0068] Among them, are the variance of Young's modulus, the variance of Poisson's ratio and the variance of crack density parameter, respectively.

[0069] Step S40: Calculate the vertical principal stress of the target layer based on the velocity field, formation density and density inversion volume.

[0070] Specifically, the signal velocity field and formation density are read based on the logging data; the vertical principal stress is calculated based on the velocity field, the formation density and a preset density inversion volume, and the calculation rule is:

[0071] Among them, σ Vis the vertical principal stress; ρ is the formation density; g w is the acceleration due to gravity; v f is the velocity field; t is the round-trip travel time.

[0072] Step S50: Calculating the maximum horizontal principal stress of the current reservoir based on the Poisson's ratio, the fracture density, the vertical principal stress and a preset principal stress obtaining model.

[0073] Specifically, the preset principal stress obtaining model is:

[0074] Among them, σ H is the maximum principal stress; ν is Poisson's ratio; g is the preset variable; and e is the crack density.

[0075] In an embodiment of the present invention, the scheme of the present invention constructs a prestack inversion equation and an objective function under the framework of Bayesian theory through a newly derived prestack orientation AVO approximate equation, and starting from the prestack OVT data set, directly inverts the Poisson's ratio and fracture density through prestack anisotropic AVO, which has the advantages of stable algorithm and high prediction accuracy. The scheme of the present invention overcomes the difficulty of obtaining vertical stress in conventional geostress calculations. By using the velocity body, the pressure generated by the overlying strata is obtained by integrating the density inversion results as the vertical stress of the reservoir, providing a basis for the calculation of the maximum horizontal principal stress, and has high computational efficiency and is simple and easy to implement. The scheme of the present invention derives a new formula for calculating the maximum horizontal principal stress based on fracture density and Poisson's ratio. The new calculation formula has the advantages of clear physical meaning, few calculation parameters, and high calculation accuracy. This can reduce the consumption of computational storage when calculating large data bodies, and can improve the computational efficiency of geostress evaluation.

[0076] Example:

[0077] A well-seismic calibration was performed on a certain detection well. The calibration results are shown in Figure 2. From left to right, they are the acoustic wave time difference curve, density curve, velocity curve, impedance curve, reflection coefficient curve, well logging layer, synthetic seismic record (red) and well-side channel seismic data (black). From shallow to deep layers, the synthetic seismic record and the well-side channel seismic waveform characteristics have good similarity and high correlation coefficient, indicating that the well-seismic calibration results are good.

[0078] The OVT processing of the corresponding detection well gather is performed, and the processing results are shown in Figure 3. In the figure, the phase axis at time 1700ms is the main waveform feature of the target layer. The phase axis of the gather fluctuates up and down with the travel time and azimuth, indicating that the phase axis has obvious anisotropy characteristics.

[0079] The corresponding detection well gathers are optimized, and the processing results are shown in Figures 4, 5 and 6. Figure 4 is the azimuth gather of the OVT spiral gather after azimuth processing, Figure 5 is the azimuth gather after random noise removal, and Figure 6 is the azimuth gather after gather flattening. The optimized OVT gathers have high signal-to-noise ratio and resolution, and the AVO characteristics of the phase axis are obvious, which can provide a data basis for pre-stack seismic inversion.

[0080] As shown in Figure 7, the stratigraphic interpretation of the shale reservoir target layer is carried out along the seismic phase axis starting from the well point under the time-depth relationship of well-seismic calibration. The structure in the study area shows the characteristics of high in the west and low in the east.

[0081] Figure 8 shows a seismic profile of a well using the present invention. The target shale reservoir layer exhibits a relatively strong peak reflection characteristic. Figures 9 and 10 show the Poisson's ratio and fracture density inversion profiles of the experimental well obtained through prestack anisotropy inversion. Figure 9 shows the Poisson's ratio inversion profile of the experimental well. Figure 10 shows the fracture density inversion profile of the experimental well. The profile results indicate a low Poisson's ratio and a high fracture density in the shale reservoir section.

[0082] Figures 11 and 12 are inversion plots of the Poisson's ratio and fracture density obtained from prestack AVO anisotropy inversion for the experimental wells. Figure 11 shows the Poisson's ratio inversion plot for the study area, which shows low Poisson's ratios in the northern part and high Poisson's ratios in the southwestern and eastern parts. Figure 12 shows the fracture density inversion plot for the study area, which shows well-developed fractures in the northwest and less developed fractures in the southeast.

[0083] Figure 13 is the maximum horizontal principal stress profile result of the experimental well. The profile shows that the maximum horizontal principal stress gradually increases from top to bottom. Figure 14 is the maximum horizontal principal stress prediction plan of the study area. The maximum horizontal principal stress in the northern and southeastern parts of the study area is higher, and the maximum horizontal principal stress in the central and southern regions is lower. The predicted results are highly consistent with the known geological knowledge, and the prediction results confirm the effectiveness and applicability of this method.

[0084] Figure 15 is a system structure diagram of a shale gas reservoir maximum horizontal principal stress prediction system provided by one embodiment of the present invention. As shown in Figure 15, an embodiment of the present invention provides a shale gas reservoir maximum horizontal principal stress prediction system, the system comprising: an acquisition unit for performing cross-arranged extraction of acquired original shot gathers to obtain corresponding OVT gathers and preprocessing the OVT gathers; a processing unit for acquiring well logging data and performing well-seismic calibration and seismic horizon interpretation based on the well logging data to obtain an interpreted horizon in the current reservoir that satisfies the in-situ stress prediction as a target horizon; a training unit for constructing a low-frequency model based on the well logging data and the target horizon, and calculating the Poisson's ratio and fracture density based on the low-frequency model; and calculating the vertical principal stress of the target horizon based on the well logging data; and an output unit for calculating the maximum horizontal principal stress of the current reservoir based on the Poisson's ratio, the fracture density, the vertical principal stress, and a preset principal stress model.

[0085] An embodiment of the present invention further provides a computer-readable storage medium having instructions stored thereon, which, when executed on a computer, enables the computer to execute the above-mentioned method for predicting the maximum horizontal principal stress of a shale gas reservoir.

[0086] Those skilled in the art will appreciate that all or part of the steps in the methods of the aforementioned embodiments can be accomplished by instructing the relevant hardware through a program, which is stored in a storage medium and includes a number of instructions for causing a single-chip microcomputer, chip, or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0087] The above describes in detail the optional embodiments of the present invention in conjunction with the accompanying drawings. However, the embodiments of the present invention are not limited to the specific details in the above embodiments. Within the technical concept of the embodiments of the present invention, a variety of simple modifications can be made to the technical solutions of the embodiments of the present invention, and these simple modifications all fall within the scope of protection of the embodiments of the present invention. It should also be noted that the various specific technical features described in the above specific embodiments can be combined in any suitable manner unless there is any contradiction. In order to avoid unnecessary repetition, the embodiments of the present invention will no longer describe the various possible combinations separately.

[0088] In addition, the various embodiments of the present invention may be arbitrarily combined, and as long as they do not violate the concept of the embodiments of the present invention, they should also be regarded as the contents disclosed in the embodiments of the present invention.

Claims

1. A method for predicting the maximum horizontal principal stress of a shale gas reservoir, characterized in that: The method comprises: Perform cross-arrangement extraction on the collected original shot gathers to obtain corresponding OVT gathers, and preprocess the OVT gathers; Collect well logging data, and perform well seismic calibration and seismic horizon interpretation based on the well logging data to obtain the interpretation horizon of the current reservoir that meets the in-situ stress prediction as the target horizon; Based on the well logging data and the target layer, a low-frequency model is constructed, and Poisson's ratio and fracture density are calculated based on the low-frequency model; Calculate the vertical principal stress of the target layer based on the velocity field, formation density and density inversion volume; The maximum horizontal principal stress of the current reservoir is obtained by calculation based on the Poisson's ratio, the fracture density, the vertical principal stress and the preset principal stress obtaining model.

2. The method according to claim 1, characterized in that The method of extracting the collected original shot gathers by cross arrangement to obtain the corresponding OVT gathers includes: The collected original shot gathers are cross-arranged and extracted to obtain seismic gathers at the same detection point under the same shot line; The cross-arranged gathers are divided into OVT units, and the divided cross-arranged gathers are extracted to obtain OVT gathers.

3. The method according to claim 1, characterized in that The preprocessing of the OVT gathers includes: Perform five-dimensional regularization processing on the OVT gather to obtain the first OVT gather; Performing migration processing on the first OVT gather to obtain a second OVT gather; Perform anisotropy correction on the second OVT gather to obtain a third OVT gather; Pick up the velocity field of the third OVT gather, and transform the third OVT gather into an incident angle gather containing azimuth information; Perform azimuth angle superposition processing on incident angle gathers in different azimuth segments to obtain azimuth gathers; One or more of Radon transform, wavelet threshold method and spectral decomposition are performed on the azimuth gather to complete OVT gather preprocessing.

4. The method according to claim 1, characterized in that After collecting the logging data, the method further includes: performing preprocessing on the logging curves in the logging data, including: Standardize each well logging curve and identify the well logging curves that lack shear wave curves; The various logging curves include: Acoustic curve, density curve, natural potential curve, gamma curve and resistivity curve; For well logging lacking shear wave curves, a model relationship between the corresponding well logging shear wave velocity and preset sensitive parameters is constructed, and the shear wave curve is obtained based on the neural network and the model relationship; Based on the shear wave curves of various well logging curves after standardization, Young's modulus curve and Poisson's ratio curve are obtained; An anisotropic theoretical model is constructed based on the Young's modulus curve and the Poisson's ratio curve to calculate the crack density curve.

5. The method according to claim 4, characterized in that The method of performing well seismic calibration and seismic horizon interpretation based on the well logging data to obtain an interpretation horizon of the current reservoir that satisfies the in-situ stress prediction as a target horizon includes: Dynamically extracting seismic wavelets based on well logging data, and continuously synthesizing seismic records based on the seismic wavelets, and taking the seismic wavelets corresponding to the synthesized seismic records whose correlation coefficients with the actual seismic records of the well bypass track are less than a preset correlation coefficient threshold as calibration results; Based on the calibration results, the current seismic horizon is structurally interpreted based on the preset knowledge graph to obtain the interpretation horizon of the current reservoir that meets the in-situ stress prediction as the target horizon.

6. The method according to claim 4, characterized in that The step of constructing a low-frequency model based on the well logging data and the target layer comprises: Taking the Young's modulus curve, the Poisson's ratio curve and the fracture density curve as the vertical source data and the structural interpretation of the target layer as the horizontal constraint, extrapolation and interpolation are performed to obtain a low-frequency model data volume of Young's modulus, Poisson's ratio and fracture density in three-dimensional space as a low-frequency model.

7. The method according to claim 6, characterized in that The calculating Poisson's ratio and crack density based on the low-frequency model comprises: Based on the low-frequency model, the seismic reflection coefficient is calculated at each incident angle and azimuth to obtain a reflection coefficient matrix; a prestack inversion model is constructed based on the seismic wavelet matrix and the reflection coefficient matrix to obtain a seismic reflection amplitude matrix; based on the Bayesian framework and the seismic reflection amplitude matrix, an objective function is constructed, and the Poisson's ratio and fracture density are calculated based on the objective function.

8. The method according to claim 7, characterized in that The seismic reflection amplitude matrix construction rule is: NMK×1 =G NMK×3K m 3K×1 . Among them, d NMK×1 is the seismic reflection amplitude matrix; G NMK×3K is the seismic wavelet matrix; m 3K×1 is the reflection coefficient matrix.

9. The method according to claim 7, characterized in that: The objective function is expressed as: Among them, σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well bypass; d is the seismic reflection amplitude matrix; G is the seismic wavelet matrix; T is the transpose of the wavelet matrix; Q is a diagonal matrix, expressed as: Among them, are the variance of Young's modulus, the variance of Poisson's ratio and the variance of crack density parameter, respectively.

10. The method according to claim 9, characterized in that Calculating the vertical principal stress of the target layer based on the well logging data includes: Reading signal velocity field and formation density based on the logging data; The vertical principal stress is calculated based on the velocity field, the formation density and the preset density inversion volume, and the calculation rule is: Among them, σ V is the vertical principal stress; ρ is the formation density; g w is the acceleration due to gravity; v f is the velocity field; t is the round-trip travel time.

11. The method according to claim 10, characterized in that The preset principal stress obtaining model is: Among them, σ H is the maximum principal stress; ν is Poisson's ratio; g is the preset variable; e is the crack density.

12. A shale gas reservoir maximum horizontal principal stress prediction system, characterized in that: The system comprises: A collection unit, used for performing cross-arrangement extraction of the collected original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; A processing unit is used to collect well logging data, and perform well seismic calibration and seismic horizon interpretation based on the well logging data to obtain an interpretation horizon of the current reservoir that satisfies the prediction of ground stress as a target horizon; Training unit for: Based on the well logging data and the target layer, a low-frequency model is constructed, and Poisson's ratio and fracture density are calculated based on the low-frequency model; Calculate the vertical principal stress of the target layer based on the velocity field, formation density and density inversion volume; The output unit is used to calculate and obtain the maximum horizontal principal stress of the current reservoir based on the Poisson's ratio, the fracture density, the vertical principal stress and a preset principal stress obtaining model.

13. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores instructions, which, when executed on a computer, enable the computer to execute the method for predicting the maximum horizontal principal stress of a shale gas reservoir as described in any one of claims 1 to 11.

Images