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

Abstract

The embodiment of the invention provides a shale gas reservoir maximum horizontal principal stress prediction method and system, and belongs to the technical field of shale gas development. The method comprises the following steps: performing cross arrangement trace gather extraction on an acquired original shot gather to obtain a corresponding OVT trace gather, and preprocessing the OVT trace gather; logging data are collected, well-seismic calibration and seismic horizon interpretation are carried out based on the logging data, and an interpretation horizon, meeting crustal stress prediction, of a current reservoir is obtained and serves as a target horizon; constructing a low-frequency model based on logging data and the target horizon, and calculating a Poisson's ratio and fracture density based on the low-frequency model; calculating the vertical principal stress of the target horizon based on the velocity field, the stratum density and the density inversion body; 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 solving model. According to the scheme, the problem that an existing shale gas reservoir maximum horizontal principal stress prediction scheme is not high in precision is solved.

Description

Technical Field

[0001] The present invention relates to the technical field of shale gas development, and particularly 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

[0002] The maximum horizontal principal stress plays an important role in the exploration and development of shale gas in shale reservoirs. In geotechnical engineering, the maximum horizontal principal stress refers to the maximum horizontal stress in the formation, which is the maximum stable stress in the formation. In shale oil and gas exploration, the magnitude of the maximum horizontal principal stress determines the fractures in the formation and the failure of rocks, as well as the migration and storage of oil and gas in the formation. In shale oil and gas development, the maximum horizontal principal stress determines the location of well placement, the direction of well trajectory, and the parameters of fracturing design.

[0003] In recent years, a large number of studies have been carried out on the prediction of the maximum horizontal principal stress, and quantitative prediction studies of in-situ stress have been carried out using seismic and logging data. For example, Liu Haojuan (in 2021) established an in-situ stress prediction model. Based on fine three-dimensional seismic interpretation and three-dimensional prestack seismic inversion, well-point data were used to simulate and select regional adaptability parameters for stress calculation, and three-dimensional simulation of the in-situ stress field was carried out to predict the direction of the maximum horizontal principal stress, the direction of the minimum horizontal principal stress, and the horizontal stress difference coefficient. Qu Yang (in 2019) used resistivity imaging logging and shear wave anisotropy logging methods to interpret the direction of the maximum horizontal principal stress in this area, calculated the minimum principal stress in the study area using the Newberry model, calculated the maximum principal stress based on dual caliper analysis, and clarified the trajectory direction of horizontal well deployment in this area, providing a basis for reservoir fracturing. Qi Dunke (in 2017) proposed a theory and method for extracting a more reliable stress field direction by analyzing microseismic monitoring data of horizontal wells, and obtained the stress field direction of the monitored well area based on the microseismic monitoring data of 3 horizontal wells in the Daqing Oilfield. Li Guangquan (in 2012) studied the method of obtaining the maximum horizontal principal stress using borehole collapse information, and proposed a method of inversely calculating the magnitude of the maximum horizontal principal stress by combining logging data analysis with numerical simulation methods, which can truly reflect the stress environment of deep formations and the mechanical properties of formation rocks.

[0004] However, in the actual application process, there are certain accuracy problems in these solutions, and the calculation results of the maximum horizontal principal stress cannot be directly used in industrial production. Based on this, it is necessary to create a prediction solution for the maximum horizontal principal stress of a shale gas reservoir with guaranteed accuracy. Summary of the Invention

[0005] 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 existing in the existing prediction solutions for the maximum horizontal principal stress of a shale gas reservoir.

[0006] To achieve the above object, a 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 cross-line binning on the collected original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; collecting logging data, and performing well-seismic calibration and seismic horizon interpretation based on the logging data to obtain an interpreted horizon that meets the requirements of in-situ stress prediction for the current reservoir as the target horizon; constructing a low-frequency model based on the logging data and the target horizon, and calculating the Poisson's ratio and fracture density based on the low-frequency model; calculating the vertical principal stress of the target horizon based on the velocity field, formation density, and density inversion body; 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 calculation model.

[0007] Optionally, the performing cross-line binning on the collected original shot gathers to obtain corresponding OVT gathers includes: performing cross-line binning on the collected original shot gathers to obtain seismic gathers of the same geophone point under the same shot line; performing OVT unit division on the gathers of each cross-line, and performing binning processing on the divided gathers of each cross-line to obtain OVT gathers.

[0008] Optionally, the preprocessing of the OVT gathers includes: performing five-dimensional regularization processing on the OVT gathers to obtain a first OVT gather; performing migration 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 from the third OVT gather and converting the third OVT gather into an incident angle gather with azimuth information; performing azimuth angle stacking processing on the incident angle gathers of different azimuth angle segments to obtain an azimuth angle gather; performing one or more of Radon transform, wavelet threshold method, and spectral decomposition on the azimuth angle gather to complete the preprocessing of the OVT gathers.

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

[0010] Optionally, calibrating well-seismic data and interpreting seismic horizons based on the logging data to obtain an interpreted horizon of the current reservoir that meets the requirements for in-situ stress prediction as the target horizon includes: dynamically extracting seismic wavelets based on the logging data, continuously synthesizing seismic records based on the seismic wavelets, and using the seismic wavelets corresponding to the seismic records with a correlation coefficient less than a preset correlation coefficient threshold with the actual seismic records of the well-side trace as the calibration result; based on the calibration result, interpreting the structure of the current seismic horizon based on a preset knowledge graph to obtain an interpreted horizon of the current reservoir that meets the requirements for in-situ stress prediction as the target horizon.

[0011] Optionally, constructing a low-frequency model based on the logging data and the target horizon includes: using the Young's modulus curve, the Poisson's ratio curve, and the fracture density curve as the longitudinal source data, and using the structural interpretation of the target horizon as the lateral constraint for extrapolation and interpolation to obtain a low-frequency model data volume of Young's modulus, Poisson's ratio data volume, and fracture density in three-dimensional space as the low-frequency model.

[0012] Optionally, calculating Poisson's ratio and fracture density based on the low-frequency model includes: calculating seismic reflection coefficients at various incident angles and azimuth angles based on the low-frequency model to obtain a reflection coefficient matrix; constructing a prestack 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 the Bayesian framework and the seismic reflection amplitude matrix, and calculating Poisson's ratio and fracture density based on the objective function.

[0013] Optionally, the rule for constructing the seismic reflection amplitude matrix is:

[0014] d NMK×1 =G NMK×3K m 3K×1 .

[0015] where 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.

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

[0017]

[0018] where σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well-side trace; d is the seismic reflection amplitude matrix; G is the seismic wavelet matrix; T is the sampling period; Q is a diagonal matrix expressed as:

[0019]

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

[0021] Optionally, calculating the vertical principal stress of the target horizon based on the logging data includes: reading the signal velocity field and formation density from the logging data; calculating the vertical principal stress based on the velocity field, the formation density, and a preset density inversion body, and the calculation rule is:

[0022]

[0023] where, σ 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 two-way travel time.

[0024] Optionally, the preset principal stress calculation model is:

[0025]

[0026] where, σ H is the maximum principal stress; ν is Poisson's ratio; g is a preset variable; e is the fracture density.

[0027] In a second aspect of the present invention, a system for predicting the maximum horizontal principal stress of a shale gas reservoir is provided. The system includes: an acquisition unit for performing cross-array trace gathering on the acquired original shot gather to obtain a corresponding OVT gather, and preprocessing the OVT gather; a processing unit for acquiring logging data, and performing well-seismic calibration and seismic horizon interpretation based on the logging data to obtain an interpreted horizon that meets the in-situ stress prediction of the current reservoir as the target horizon; a training unit for: constructing a low-frequency model based on the logging data and the target horizon, and calculating Poisson's ratio and fracture density based on the low-frequency model; calculating the vertical principal stress of the target horizon based on the velocity field, formation density, and density inversion body; 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 calculation model.

[0028] On the other hand, the present invention provides a computer-readable storage medium, which stores instructions thereon that, when run on a computer, cause the computer to execute the above-mentioned method for predicting the maximum horizontal principal stress of a shale gas reservoir.

[0029] Through the above technical solutions, the solution of the present invention first derives a prediction formula for the maximum horizontal principal stress expressed in terms of fracture density and Poisson's ratio. At the same time, starting from the theory of prestack anisotropic media, a longitudinal wave azimuth AVO approximate formula expressed in terms of Young's modulus, Poisson's ratio, and fracture density is established. Under the framework of Bayesian theory, an inversion equation is established using the new approximate equation to achieve prestack inversion of seismic parameters. This method improves the prediction accuracy of the maximum horizontal principal stress through direct inversion.

[0030] Other features and advantages of the embodiments of the present invention will be described in detail in the subsequent specific embodiments section. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] The 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 specific embodiments, they are used to explain the embodiments of the present invention, but do not constitute a limitation to the embodiments of the present invention. In the drawings:

[0032] Figure 1 is a flow chart of the steps for predicting the maximum horizontal principal stress of a shale gas reservoir provided by an embodiment of the present invention;

[0033] Figure 2 is the well-seismic calibration result of an experimental well provided by an embodiment of the present invention;

[0034] Figure 3 is the OVT domain spiral trace gather of Inline 19325 processed by OVT provided by an embodiment of the present invention;

[0035] Figure 4 is the azimuthal angle trace gather obtained by azimuthal processing of the OVT spiral trace gather provided by an embodiment of the present invention;

[0036] Figure 5 is the azimuthal angle trace gather after random noise removal of the azimuthal angle trace gather provided by an embodiment of the present invention;

[0037] Figure 6 is the azimuthal angle trace gather after flattening the trace gather provided by an embodiment of the present invention;

[0038] Figure 7 is the horizon interpretation result of the target layer of the shale gas reservoir provided by an embodiment of the present invention;

[0039] Figure 8 is the seismic profile of the experimental well provided by an embodiment of the present invention;

[0040] Figure 9 is the Poisson's ratio inversion profile passing through the experimental well provided by an embodiment of the present invention;

[0041] Figure 10It is the fracture density inversion profile across the experimental well provided by an embodiment of the present invention;

[0042] Figure 11 It is the Poisson's ratio inversion plan view of the study area provided by an embodiment of the present invention;

[0043] Figure 12 It is the fracture density inversion plan view of the study area provided by an embodiment of the present invention;

[0044] Figure 13 It is the maximum horizontal principal stress profile result across the experimental well provided by an embodiment of the present invention;

[0045] Figure 14 It is the maximum horizontal principal stress magnitude prediction plan view provided by an embodiment of the present invention;

[0046] Figure 15 It is the system structure diagram of the maximum horizontal principal stress prediction system for shale gas reservoirs provided by an embodiment of the present invention. Specific Embodiments

[0047] The following further elaborates on the specific embodiments of the present invention with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are only for the purpose of illustrating and explaining the present invention, and are not intended to limit the present invention.

[0048] Figure 1 It is the method flow chart of the maximum horizontal principal stress prediction method for shale gas reservoirs provided by an embodiment of the present invention. As Figure 1 shown, an embodiment of the present invention provides a method for predicting the maximum horizontal principal stress of a shale gas reservoir, and the method includes:

[0049] Step S10: Perform cross-line trace gathering on the collected original shot gathers to obtain the corresponding OVT gathers, and preprocess the OVT gathers.

[0050] Specifically, the performing cross-line trace gathering on the collected original shot gathers to obtain the corresponding OVT gathers includes: performing cross-line trace gathering on the collected original shot gathers to obtain the seismic trace gathers of the same geophone point under the same shot line; performing OVT unit division on the trace gathers of each cross-line, and performing trace gathering processing on the divided trace gathers of each cross-line to obtain the OVT gathers.

[0051] In the embodiments of the present invention, cross-line sorting is performed on the original shot gathers collected in the field to obtain seismic gathers of the same shot line and the same geophone point. Then, the gathers obtained by cross-line sorting are divided into OVT units. Based on the division, shot gathers are sorted from the OVT units to obtain OVT gathers. The OVT gathers are subjected to five-dimensional regularization processing, omnidirectional migration processing in the OVT domain, and anisotropic correction of the gathers to obtain corrected OVT gathers, providing a data basis for pre-stack fracture and in-situ stress prediction.

[0052] Preferably, the preprocessing of the OVT gathers includes: performing five-dimensional regularization processing on the OVT gathers to obtain first OVT gathers; performing migration processing on the first OVT gathers to obtain second OVT gathers; performing anisotropic correction on the second OVT gathers to obtain third OVT gathers; picking up the velocity field of the third OVT gathers and converting the third OVT gathers into incident angle gathers containing azimuth information; performing azimuth-angle stacking processing on the incident angle gathers of different azimuth angle segments to obtain azimuth gathers; and performing one or more of Radon transform, wavelet threshold method, and spectral decomposition on the azimuth gathers to complete the preprocessing of the OVT gathers.

[0053] Specifically, the velocity field is picked up and utilized for the obtained OVT gathers, and they are converted into incident angle gathers containing azimuth information. According to the incident angle gathers of different azimuth angle segments, azimuth-angle stacking processing is performed on them to obtain incident angle gathers of different azimuths (abbreviation: azimuth-separated gathers). Aiming at the problems of multiple wave interference, low signal-to-noise ratio, and low resolution existing in the azimuth gathers, methods such as high-precision Radon transform, wavelet threshold method, and spectral decomposition are used to improve the quality of the gather data.

[0054] Step S20: Collect logging data, and perform well-seismic calibration and seismic horizon interpretation based on the logging data to obtain the interpreted horizons of the current reservoir that meet the in-situ stress prediction as the target horizons.

[0055] Preferably, after collecting the logging data, the method further includes: performing preprocessing on the logging curves in the logging data, including: performing standardization processing on each logging curve and identifying the logging lacking shear wave curves; wherein, each logging curve includes: acoustic wave curve, density curve, spontaneous potential curve, gamma curve, and resistivity curve; for the logging lacking shear wave curves, constructing a model relationship between the shear wave velocity of the corresponding logging and preset sensitive parameters, and obtaining the shear wave curve based on the neural network and the model relationship; obtaining Young's modulus curve and Poisson's ratio curve based on the shear wave curves of each standardized logging curve; 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.

[0056] Specifically, standardize the logging curves in the work area, such as acoustic curves, density curves, spontaneous potential curves, gamma curves, resistivity curves, etc. For wells lacking shear wave curves, a model relationship between shear wave velocity and sensitive parameters can be constructed, and the shear wave velocity curve can be obtained through a neural network algorithm, the Young's modulus and Poisson's ratio curves can be obtained, and an anisotropic theoretical model can be established to calculate the fracture density curve.

[0057] Further, perform well-seismic calibration and seismic horizon interpretation based on the logging data to obtain the interpreted horizon that meets the in-situ stress prediction for the current reservoir as the target horizon, including: dynamically extracting seismic wavelets based on the logging data, continuously synthesizing seismic records based on the seismic wavelets, and using the seismic wavelets corresponding to the seismic records with a correlation coefficient less than the preset correlation coefficient threshold with the actual seismic records of the well side trace as the calibration result; based on the calibration result, perform structural interpretation on the current seismic horizon based on a preset knowledge graph to obtain the interpreted horizon that meets the in-situ stress prediction for the current reservoir as the target horizon.

[0058] Specifically, through dynamic seismic wavelet extraction, synthetic seismic record production and drift, and comparison between the synthetic seismic record and the well side trace, when the correlation coefficient between the two reaches the preset requirement, the seismic wavelet and the time-depth relationship at this time are the final calibration results. Based on the calibration results, clarify the target layer as the seismic response characteristics, and combine the existing structural understanding to interpret the target horizon to obtain the interpreted horizon that meets the in-situ stress prediction, preventing the phenomenon of the horizon jumping up and down in three-dimensional space.

[0059] Step S30: Construct a low-frequency model based on the logging data and the target horizon, and calculate the Poisson's ratio and fracture density based on the low-frequency model.

[0060] Specifically, using the Young's modulus curve, the Poisson's ratio curve, and the fracture density curve as the longitudinal source data, and the structural interpretation of the target horizon as the transverse constraint, perform extrapolation and interpolation to obtain the Young's modulus, Poisson's ratio data volume, and the low-frequency model data volume of the fracture density in three-dimensional space as the low-frequency model.

[0061] In the embodiment of the present invention, under the control of the time-depth curve obtained in step S20, using the Young's modulus curve, the Poisson's ratio curve, and the fracture density curve as the longitudinal source data, and under the transverse constraint of the target horizon data, perform extrapolation and interpolation to obtain the Young's modulus, Poisson's ratio data volume, and the low-frequency model data volume of the fracture density in three-dimensional space.

[0062] Further, calculating the Poisson's ratio and fracture density based on the low-frequency model includes: calculating the seismic reflection coefficients at various incident angles and azimuth angles based on the low-frequency model to obtain a reflection coefficient matrix; constructing a prestack 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 the Bayesian framework and the seismic reflection amplitude matrix, and calculating the Poisson's ratio and fracture density based on the objective function.

[0063] Specifically, use the Young's modulus, Poisson's ratio data volume and fracture density low-frequency model data volume. Calculate the seismic reflection coefficient, and the calculation rule is:

[0064]

[0065] where θ is the incident angle, φ is the azimuth angle, a is the power exponent of density with respect to the longitudinal wave velocity, g is the square of the transverse-to-longitudinal wave velocity ratio, here a and g are constants, and E, ν, and e represent Young's modulus, Poisson's ratio, and fracture density respectively.

[0066] Further, for the seismic reflection amplitudes at different incident angles and azimuth angles, the coefficient matrix included multiplies the vector of the parameters to be inverted to obtain the seismic reflection coefficient, and the seismic reflection amplitude multiplied by the wavelet matrix constructs a prestack inversion equation:

[0067] d NMK×1 =G NMK×3K m 3K×1 .

[0068] where 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. Among them:

[0069]

[0070] where N is the number of incident angles, M is the number of azimuth angles, K is the number of sampling points, wvlt is the wavelet matrix. S(θ i ,φ j ) is the seismic reflection amplitude with an incident angle of θ i and an azimuth angle of φ j .

[0071] Further:

[0072]

[0073]

[0074]

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

[0076]

[0077] where σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well-side trace; d is the seismic reflection amplitude matrix; G is the seismic wavelet matrix; T is the sampling period; Q is a diagonal matrix, expressed as:

[0078]

[0079] where are the variances of Young's modulus, Poisson's ratio, and fracture density parameter, respectively.

[0080] Step S40: Calculate the vertical principal stress of the target horizon based on the velocity field, formation density, and density inversion body.

[0081] Specifically, read the signal velocity field and formation density based on the well logging data; calculate the vertical principal stress based on the velocity field, the formation density, and a preset density inversion body, and the calculation rule is:

[0082]

[0083] where σ 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 two-way travel time.

[0084] Step S50: Calculate 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 calculation model.

[0085] Specifically, the preset principal stress calculation model is:

[0086]

[0087] where σ H is the maximum principal stress; ν is the Poisson's ratio; g is a preset variable; e is the fracture density.

[0088] In the embodiments of the present invention, the solution of the present invention constructs a prestack inversion equation and an objective function within the framework of Bayesian theory through a newly derived prestack azimuth AVO approximate equation. Starting from the prestack OVT gather, Poisson's ratio and fracture density are directly obtained through prestack anisotropic AVO inversion, which has the advantages of stable algorithm and high prediction accuracy. The solution of the present invention overcomes the problem that it is difficult to obtain the vertical stress in conventional in-situ stress calculation. By using the velocity volume, the pressure generated by the overlying strata is obtained through the integration of the density inversion result as the vertical stress of the reservoir, providing a basis for the calculation of the maximum horizontal principal stress, and having high calculation efficiency and being simple and easy to implement. The solution of the present invention derives a new formula for calculating the maximum horizontal principal stress expressed in terms of 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 calculation storage during large data volume calculations and improve the calculation efficiency of in-situ stress evaluation.

[0089] Embodiment:

[0090] Perform well-seismic calibration for a certain exploration well, and the calibration results are as Figure 2 shown. From left to right are the acoustic travel time curve, density curve, velocity curve, impedance curve, reflection coefficient curve, well logging stratification, synthetic seismogram (red) and seismic data of the well-side trace (black). From shallow to deep, the synthetic seismogram and the seismic waveform characteristics of the well-side trace have good similarity and a high correlation coefficient, indicating that the well-seismic calibration results are good.

[0091] Perform OVT processing on the gather of the corresponding exploration well, and the processing results are as Figure 3 shown. The in-phase axis at 1700 ms in time in the figure is the waveform characteristic of the main target layer. The in-phase axis of the gather fluctuates up and down with travel time and azimuth angle, indicating that the in-phase axis has obvious anisotropic characteristics.

[0092] Perform optimization processing on the gather of the corresponding exploration well, and the processing results are as Figure 4 、 Figure 5 and Figure 6 shown. Figure 4 is the azimuth angle gather after azimuth processing of the OVT spiral gather, Figure 5 is the azimuth angle gather after removing random noise from the azimuth angle gather, Figure 6 is the azimuth angle gather after flattening the gather. The optimized OVT gather has a high signal-to-noise ratio and resolution, and the AVO characteristics of the in-phase axis are obvious, providing a data basis for prestack seismic inversion.

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

[0094] AsFigure 8 The seismic profile of an exploration well implementing the solution of the present invention. It can be seen from the seismic profile that the target layer of the shale reservoir has relatively strong wave peak reflection characteristics. Figure 9 and Figure 10 are the inversion profiles of Poisson's ratio and fracture density passing through the experimental well obtained by pre-stack anisotropic inversion. Figure 9 is the inversion profile of Poisson's ratio of the experimental well. Figure 10 is the inversion profile of fracture density of the experimental well. The profile results show that the shale reservoir section has a relatively low Poisson's ratio and a relatively high fracture density.

[0095] Figure 11 and Figure 12 are the inversion plane maps of Poisson's ratio and fracture density passing through the experimental well obtained by pre-stack AVO anisotropic inversion. Figure 11 is the inversion plane map of Poisson's ratio in the study area. The Poisson's ratio in the northern part of the study area is relatively low, and the Poisson's ratio in the southwestern and eastern regions is relatively high. Figure 12 is the inversion plane map of fracture density in the study area. The fractures are relatively developed in the northwestern part of the study area and less developed in the southeastern part.

[0096] Figure 13 is the profile result of the maximum horizontal principal stress passing through the experimental well. The profile shows that the maximum horizontal principal stress gradually increases from top to bottom. Figure 14 is the predicted plane map of the maximum horizontal principal stress in the study area. The maximum horizontal principal stress is relatively high in the northern and southeastern parts of the study area, and relatively low in the central and southern regions. The prediction results are highly consistent with the known geological understanding, which verifies the effectiveness and applicability of this method.

[0097] Figure 15 is the system structure diagram of the maximum horizontal principal stress prediction system for shale gas reservoirs provided by an embodiment of the present invention. As Figure 15 shown, the embodiment of the present invention provides a maximum horizontal principal stress prediction system for shale gas reservoirs. The system includes: an acquisition unit, which is used to perform cross-line extraction of shot gathers on the collected original shot gathers to obtain corresponding OVT gathers, and preprocess the OVT gathers; a processing unit, which is used to collect logging data, and perform well-seismic calibration and seismic horizon interpretation based on the logging data to obtain the interpreted horizons that meet the in-situ stress prediction for the current reservoir as the target horizons; a training unit, which is used to: construct a low-frequency model based on the logging data and the target horizons, and calculate Poisson's ratio and fracture density based on the low-frequency model; calculate the vertical principal stress of the target horizons based on the logging data; an output unit, which is used to calculate 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 calculation model.

[0098] An embodiment of the present invention also provides a computer-readable storage medium, on which instructions are stored, and when the instructions run on a computer, the computer is caused to execute the above-mentioned method for predicting the maximum horizontal principal stress of a shale gas reservoir.

[0099] Those skilled in the art can understand that all or part of the steps in the methods of the above embodiments can be completed by instructing relevant hardware through a program. The program is stored in a storage medium, including several instructions to cause a single-chip microcomputer, a chip or a 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 that can store program codes 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 disc.

[0100] The optional embodiments of the present invention have been described in detail above 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 scope of the technical concept of the embodiments of the present invention, various simple modifications can be made to the technical solutions of the embodiments of the present invention, and these simple modifications all fall within the protection scope of the embodiments of the present invention. Additionally, it should be noted that, among the various specific technical features described in the above specific embodiments, they can be combined in any appropriate manner without conflict. To avoid unnecessary repetition, the embodiments of the present invention do not separately describe various possible combination methods.

[0101] In addition, any combination can be made among the various different embodiments of the present invention as long as it does not violate the idea of the embodiments of the present invention, and it should also be regarded as the content disclosed by 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 includes: Performing cross-line binning on the collected original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; Collecting logging data, and performing well-seismic calibration and seismic horizon interpretation based on the logging data to obtain the interpreted horizons that meet the in-situ stress prediction for the current reservoir as the target horizons; Constructing a low-frequency model based on the logging data and the target horizons, and calculating the Poisson's ratio and fracture density based on the low-frequency model; Calculating the vertical principal stress of the target horizons based on the velocity field, formation density, and density inversion volume; 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 calculation model.

2. The method according to claim 1, characterized in that, the performing cross-line binning on the collected original shot gathers to obtain corresponding OVT gathers includes: Performing cross-line binning on the collected original shot gathers to obtain the seismic gathers of the same geophone point under the same shot line; Performing OVT unit division on the gathers of each cross-line, and performing binning processing on the divided gathers of each cross-line to obtain OVT gathers.

3. The method according to claim 1, characterized in that, the preprocessing of the OVT gathers includes: Performing five-dimensional regularization processing on the OVT gathers to obtain the first OVT gathers; Performing migration processing on the first OVT gathers to obtain the second OVT gathers; Performing anisotropy correction on the second OVT gathers to obtain the third OVT gathers; Picking up the velocity field from the third OVT gathers and converting the third OVT gathers into incident angle gathers with azimuth information; Performing azimuthal angle stacking processing on the incident angle gathers of different azimuth segments to obtain azimuth gathers; Performing one or more of Radon transform, wavelet threshold method, and spectral decomposition on the azimuth gathers to complete the preprocessing of the OVT gathers.

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: Performing standardization processing on each logging curve and identifying the logging lacking shear wave curves; where the each logging curve includes: acoustic curve, density curve, spontaneous potential curve, gamma curve, and resistivity curve; For the logging lacking shear wave curves, constructing the model relationship between the shear wave velocity of the corresponding logging and the 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 Poisson's ratio curve based on the shear wave curves of each standardized logging curve; Constructing an anisotropic theoretical model based on the Young's modulus curve and the Poisson's ratio curve to calculate the fracture density curve.

5. The method according to claim 4, characterized in that, the performing well-seismic calibration and seismic horizon interpretation based on the logging data to obtain the interpreted horizons that meet the in-situ stress prediction for the current reservoir as the target horizons includes: Dynamically extract seismic wavelets based on logging data, continuously synthesize seismic records based on the seismic wavelets, and use the seismic wavelets corresponding to the synthesized seismic records with a correlation coefficient less than a preset correlation coefficient threshold with the actual seismic records of the well-side trace as the calibration result; Based on the calibration result, perform structural interpretation on the current seismic horizon based on a preset knowledge graph, and obtain an interpreted horizon where the current reservoir meets the in-situ stress prediction as the target horizon.

6. The method according to claim 4, wherein, the constructing a low-frequency model based on the logging data and the target horizon includes: Using the Young's modulus curve, the Poisson's ratio curve, and the fracture density curve as the longitudinal source data, and using the structural interpretation of the target horizon as the lateral constraint, perform extrapolation interpolation to obtain a low-frequency model data volume of Young's modulus, Poisson's ratio data volume, and fracture density in three-dimensional space as the low-frequency model.

7. The method according to claim 6, wherein, the calculating Poisson's ratio and fracture density based on the low-frequency model includes: Based on the low-frequency model, calculate the seismic reflection coefficients at various incident angles and azimuth angles to obtain a reflection coefficient matrix; Construct a prestack inversion model 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, construct an objective function, and calculate Poisson's ratio and fracture density based on the objective function.

8. The method according to claim 7, wherein, the rule for constructing the seismic reflection amplitude matrix is: d NMK×1 = G NMK×3K m 3K×1 . where 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, wherein, the objective function is expressed as: where σ n is the noise variance between the synthetic seismogram and the actual seismic record of the well-side trace; 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: where are the variances of Young's modulus, Poisson's ratio, and fracture density parameter respectively.

10. The method according to claim 9, wherein, the calculating the vertical principal stress of the target horizon based on the logging data includes: Read the signal velocity field and formation density based on the logging data; Calculate the vertical principal stress based on the velocity field, the formation density, and a 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 of gravity; v f is the velocity field; t is the two-way travel time.

11. The method according to claim 10, wherein, the preset principal stress calculation model is: Among them, σ H is the maximum principal stress; ν is Poisson's ratio; g is a preset variable; e is the fracture density.

12. A system for predicting the maximum horizontal principal stress of a shale gas reservoir, wherein, the system includes: An acquisition unit for performing cross-array trace gathering on the acquired original shot gathers to obtain corresponding OVT gathers, and preprocessing the OVT gathers; A processing unit for acquiring logging data, and performing well-seismic calibration and seismic horizon interpretation based on the logging data to obtain an interpreted horizon where the current reservoir meets the in-situ stress prediction as the target horizon; A training unit for: Constructing a low-frequency model based on the logging data and the target horizon, and calculating Poisson's ratio and fracture density based on the low-frequency model; Calculating the vertical principal stress of the target horizon based on the velocity field, formation density, and density inversion volume; The output unit is used to calculate 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, It is 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