CN120143231A - Method and system for predicting maximum horizontal principal stress of shale gas reservoir - Google Patents
Method and system for predicting maximum horizontal principal stress of shale gas reservoir Download PDFInfo
- Publication number
- CN120143231A CN120143231A CN202311713818.3A CN202311713818A CN120143231A CN 120143231 A CN120143231 A CN 120143231A CN 202311713818 A CN202311713818 A CN 202311713818A CN 120143231 A CN120143231 A CN 120143231A
- Authority
- CN
- China
- Prior art keywords
- seismic
- ovt
- principal stress
- density
- gather
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明实施例提供一种页岩气储层最大水平主应力预测方法及系统，属于页岩气开发技术领域。所述方法包括：对采集的原始炮集进行十字排列抽道集，获得对应的OVT道集，并对所述OVT道集进行预处理；采集测井数据，并基于所述测井数据进行井震标定和进行地震层位解释，获得当前储层满足地应力预测的解释层位，作为目标层位；基于测井数据和所述目标层位，构建低频模型，并基于低频模型计算泊松比和裂缝密度；基于速度场、地层密度和密度反演体计算目标层位的垂向主应力；基于泊松比、裂缝密度、垂向主应力和预设主应力求取模型计算获得当前储层的最大水平主应力。本发明方案解决了现有页岩气储层最大水平主应力预测方案存在的精准度不高的问题。
The embodiment of the present invention provides a method and system for predicting the maximum horizontal principal stress of a shale gas reservoir, and belongs to the technical field of shale gas development. The method includes: performing cross-arranged extraction of the collected original shot gathers to obtain the corresponding OVT gathers, 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 the interpretation layer of the current reservoir that meets the ground stress prediction as the target layer; based on the well logging data and the target layer, constructing a low-frequency model, 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 calculating the maximum horizontal principal stress of the current reservoir based on the Poisson's ratio, fracture density, vertical principal stress and preset principal stress model. The scheme of the present invention solves the problem of low accuracy of the existing maximum horizontal principal stress prediction scheme for shale gas reservoirs.
Description
Technical Field
The invention relates to the technical field of shale gas development, in particular to a shale gas reservoir maximum horizontal principal stress prediction method and a shale gas reservoir maximum horizontal principal stress prediction system.
Background
The maximum horizontal principal stress plays an important role in shale gas exploration and development of shale reservoirs. In geological engineering, the maximum horizontal principal stress refers to the maximum horizontal stress in the stratum, which is the maximum stable stress in the stratum, and in shale oil and gas exploration, the magnitude of the maximum horizontal principal stress determines fracture and rock damage in the stratum and migration and storage of oil and gas in the stratum, and in shale oil and gas development, the maximum horizontal principal stress determines parameters for well distribution position, well track direction and fracturing design.
In recent years, a great deal of research has been conducted on maximum level principal stress prediction, and quantitative prediction research of ground stress is conducted by using seismic and logging data prediction. For example, liujuan (2021) establishes a ground stress prediction model, based on fine three-dimensional seismic interpretation and three-dimensional seismic prestack inversion, utilizes well point data simulation to select regional adaptability parameters of stress calculation, develops three-dimensional simulation of a ground stress field, and predicts a maximum horizontal main stress direction, a minimum horizontal main stress direction and a horizontal stress difference coefficient. Qu Yang (2019) application of resistivity imaging logging and transverse wave anisotropic logging methods explains the direction of the maximum horizontal principal stress of the area, a Newberry model is adopted to calculate the minimum principal stress of a research area, the maximum principal stress is calculated based on double-well diameter analysis, the track direction of the horizontal well deployment of the area is defined, and the basis is provided for reservoir fracturing. qi Duke (2017) proposes a theory and a method for extracting a more reliable stress field direction by analyzing microseism monitoring data of a horizontal well, and obtaining the stress field direction of the monitoring well region according to microseism monitoring data of 3 horizontal wells of Daqing oilfield. Li Guangquan (2012) researches a method for solving the maximum horizontal main stress by using borehole collapse information, and provides a method for inverting the maximum horizontal main stress by adopting a numerical simulation method in combination with well logging data analysis, so that the stress environment of deep stratum and the mechanical properties of stratum rock can be truly reflected.
However, in the practical application process, the schemes have certain accuracy problems, and the calculation result of the maximum level main stress cannot be directly used for industrial production. Based on the method, a shale gas reservoir maximum level principal stress prediction scheme capable of ensuring accuracy needs to be created.
Disclosure of Invention
The embodiment of the invention aims to provide a method and a system for predicting the maximum horizontal principal stress of a shale gas reservoir, which at least solve the problem of low accuracy in the existing scheme for predicting the maximum horizontal principal stress of the shale gas reservoir.
In order to achieve the aim, the first aspect of the invention provides a shale gas reservoir maximum horizontal principal stress prediction method, which comprises the steps of carrying out cross arrangement gather on an acquired original shot gather to obtain a corresponding OVT gather, carrying out preprocessing on the OVT gather, acquiring well logging data, carrying out well shock calibration and seismic horizon interpretation based on the well logging data to obtain an interpretation horizon of a current reservoir meeting the earth stress prediction as a target horizon, constructing a low-frequency model based on the well logging data and the target horizon, calculating poisson ratio and crack density based on the low-frequency model, calculating vertical principal stress of the target horizon based on a velocity field, the stratum density and a density inversion body, and obtaining the maximum horizontal principal stress of the current reservoir based on the poisson ratio, the crack density, the vertical principal stress and a preset principal stress calculation model.
Optionally, performing cross arrangement channel extraction on the collected original shot set to obtain a corresponding OVT channel set, wherein the cross arrangement channel extraction is performed on the collected original shot set to obtain a seismic channel set of the same detection point under the same shot line, the OVT unit division is performed on each cross arrangement channel set, and channel extraction processing is performed on each cross arrangement channel set after division to obtain the OVT channel set.
Optionally, preprocessing the OVT gathers, wherein the preprocessing comprises five-dimensional regularization processing is performed on the OVT gathers to obtain a first OVT gather, offset processing is performed on the first OVT gather to obtain a second OVT gather, anisotropic correction is performed on the second OVT gather to obtain a third OVT gather, a speed field is picked up on the third OVT gather, the third OVT gather is converted into an incidence angle gather containing azimuth information, azimuth angle gather processing is performed on different azimuth angle section incidence angle gathers to obtain azimuth angle gathers, and one or more of radon transform, wavelet thresholding and spectral decomposition are performed on the azimuth angle gathers to complete the preprocessing of the OVT gathers.
Optionally, after logging data are acquired, the method further comprises the steps of performing preprocessing on logging curves in the logging data, performing standardization processing on the logging curves, identifying logs lacking a transverse wave curve, wherein the logging curves comprise acoustic curves, density curves, natural potential curves, gamma curves and resistivity curves, constructing model relations between corresponding logging transverse wave speeds and preset sensitive parameters for the logs lacking the transverse wave curve, solving the transverse wave curves based on a neural network and the model relations, solving Young modulus curves and Poisson ratio curves based on the transverse wave curves of the logging curves after standardization processing, and constructing an anisotropic theoretical model based on the Young modulus curves and the Poisson ratio curves to calculate a crack density curve.
Optionally, the method comprises the steps of carrying out well seismic calibration and seismic horizon interpretation based on the logging data to obtain an interpretation horizon of the current reservoir meeting the ground stress prediction, wherein the interpretation horizon is used as a target horizon, and comprises the steps of dynamically extracting seismic wavelets based on the logging data, continuously synthesizing seismic records based on the seismic wavelets, taking the seismic wavelets corresponding to the seismic records of which the actual seismic record correlation coefficient of the well side channel is smaller than a preset correlation coefficient threshold value as a calibration result, and carrying out structural interpretation on the current seismic horizon based on a preset knowledge graph based on the calibration result to obtain the interpretation horizon of the current reservoir meeting the ground stress prediction as the target horizon.
Optionally, the constructing a low-frequency model based on the logging data and the target horizon includes using the young modulus curve, the poisson ratio curve and the fracture density curve as longitudinal source data, using the construction of the target horizon as transverse constraint, and performing extrapolation interpolation to obtain a young modulus, a poisson ratio data body and a low-frequency model data body of the fracture density in a three-dimensional space as a low-frequency model.
Optionally, the Poisson's ratio and the crack density are calculated based on the low frequency model, and the Poisson's ratio and the crack density are calculated based on the low frequency model by performing seismic reflection coefficient calculation under each incident angle and azimuth angle to obtain a reflection coefficient matrix, constructing a pre-stack inversion model based on a 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 the crack density based on the objective function.
Optionally, the seismic reflection amplitude matrix construction rule is:
dNMK×1＝GNMK×3Km3K×1.
Wherein d NMK×1 is an earthquake reflection amplitude matrix, G NMK×3K is an earthquake wavelet matrix, and m 3K×1 is a reflection coefficient matrix.
Optionally, the objective function is expressed as:
wherein σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well side channel, d is the seismic reflection amplitude matrix, G is the seismic wavelet matrix, T is the sampling period, Q is the diagonal matrix, expressed as:
wherein the variance of Young's modulus, the variance of Poisson's ratio and the variance of fracture density parameters are respectively.
Optionally, the method for calculating the vertical principal stress of the target horizon based on the logging data comprises the steps of reading a signal speed field and a stratum density based on the logging data, and calculating the vertical principal stress based on the speed field, the stratum density and a preset density inversion body, wherein the calculation rule is as follows:
Wherein sigma V is vertical main stress, ρ is stratum density, g w is gravitational acceleration, v f is velocity field, and t is during double journey.
Optionally, the preset principal stress solving model is:
Wherein sigma H is the maximum principal stress, v is poisson's ratio, g is a preset variable, and e is crack density.
The invention provides a shale gas reservoir maximum horizontal principal stress prediction system which comprises an acquisition unit, a processing unit, a training unit, a speed field, a stratum density and a density inversion body, wherein the acquisition unit is used for carrying out cross arrangement of an acquired original shot gather to obtain a corresponding OVT gather and preprocessing the OVT gather, the processing unit is used for acquiring logging data, carrying out well shock calibration and seismic horizon interpretation based on the logging data to obtain an interpretation horizon of a current reservoir meeting ground stress prediction as a target horizon, the training unit is used for constructing a low-frequency model based on the logging data and the target horizon and calculating a Poisson ratio and a crack density based on the low-frequency model, calculating the vertical principal stress of the target horizon based on the speed field, the stratum density and the density inversion body, and the output unit is used for calculating and obtaining the maximum horizontal principal stress of the current reservoir based on the Poisson ratio, the crack density, the vertical principal stress and a preset principal stress calculation model.
In another aspect, the present invention provides a computer readable storage medium having instructions stored thereon that when run on a computer cause the computer to perform the shale gas reservoir maximum level principal stress prediction method described above.
Through the technical scheme, the maximum horizontal principal stress prediction formula based on the crack density and the poisson ratio is firstly deduced, meanwhile, from the prestack anisotropic medium theory, a longitudinal wave azimuth AVO approximate formula based on the representation of Young modulus, poisson ratio and crack density is established, and an inversion equation is established by utilizing a new approximate equation under the framework of Bayesian theory, so that prestack inversion of seismic parameters is realized. The method improves the prediction precision of the maximum horizontal principal stress through direct inversion.
Additional features and advantages of embodiments of the invention will be set forth in the detailed description which follows.
Drawings
The accompanying drawings are included to provide a further understanding of embodiments of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain, without limitation, the embodiments of the invention. In the drawings:
FIG. 1 is a flowchart of the steps for maximum horizontal principal stress prediction of a shale gas reservoir provided by one embodiment of the present invention;
FIG. 2 is a graph showing the results of a well shock calibration of an experimental well provided by one embodiment of the present invention;
FIG. 3 is an OVT domain spiral gather over an OVT treated line19325 line provided in accordance with one embodiment of the present invention;
FIG. 4 is a split azimuth gather of an OVT spiral gather with split azimuth processing provided by one embodiment of the present invention;
FIG. 5 is a plot of the sub-azimuth gathers after random noise removal according to one embodiment of the present invention;
FIG. 6 is a split azimuth gather flattened through a gather provided by one embodiment of the present invention;
FIG. 7 is a horizon interpretation of a shale gas reservoir destination provided by one embodiment of the present invention;
FIG. 8 is a seismic section of an experimental well provided by one embodiment of the invention;
FIG. 9 is a Poisson's ratio inversion profile of a well over experiment provided by one embodiment of the present invention;
FIG. 10 is a fracture density inversion profile of a well over-tested provided in one embodiment of the present invention;
FIG. 11 is a plot of the Poisson's ratio inversion of a region of interest provided by one embodiment of the present invention;
FIG. 12 is a plan view of an inversion of fracture density in a region of interest provided by an embodiment of the present invention;
FIG. 13 is a graph of maximum horizontal principal stress profile results for an experimental well provided by one embodiment of the invention;
FIG. 14 is a plan view of a predicted maximum horizontal principal stress magnitude provided by an embodiment of the present invention;
fig. 15 is a system structural diagram of a shale gas reservoir maximum horizontal principal stress prediction system provided by an embodiment of the invention.
Detailed Description
The following describes specific embodiments of the present invention in detail with reference to the drawings. It should be understood that the detailed description and specific examples, while indicating and illustrating the invention, are not intended to limit the invention.
Fig. 1 is a flow chart of a method for predicting maximum horizontal principal stress of a shale gas reservoir according to an embodiment of the invention. As shown in fig. 1, an embodiment of the present invention provides a method for predicting maximum horizontal principal stress of a shale gas reservoir, the method comprising:
and S10, carrying out cross arrangement channel extraction on the acquired original gun set to obtain a corresponding OVT channel set, and preprocessing the OVT channel set.
The method comprises the steps of carrying out cross arrangement gather on an acquired original shot set to obtain a corresponding OVT gather, carrying out cross arrangement gather on the acquired original shot set to obtain a seismic gather of the same detection point under the same shot line, carrying out OVT unit division on each cross arrangement gather, and carrying out gather extraction processing on each cross arrangement gather after division to obtain the OVT gather.
In the embodiment of the invention, the original shot set collected in the field is subjected to cross arrangement and gather to obtain the seismic gather of the same shot line and the same detection point, then the gather which are cross arranged is subjected to OVT unit division, the gather of the OVT unit is processed on the basis of division to obtain the gather of the OVT, the gather of the OVT is subjected to five-dimensional regularization, the gather of the OVT domain is subjected to omnibearing offset processing, and the gather is subjected to anisotropic correction to obtain the corrected gather of the OVT, so that a data basis is provided for pre-stack crack and ground stress prediction.
Preferably, the pre-processing of the OVT gathers comprises the steps of performing five-dimensional regularization processing on the OVT gathers to obtain a first OVT gather, performing offset processing on the first OVT gather to obtain a second OVT gather, performing anisotropic correction on the second OVT gather to obtain a third OVT gather, picking up a speed field on the third OVT gather, converting the third OVT gather into an incidence angle gather containing azimuth information, performing azimuth angle superposition processing on incidence angle gathers of different azimuth sections to obtain an azimuth gather, and performing one or more of radon transform, wavelet thresholding and spectral decomposition on the azimuth gather to complete the pre-processing of the OVT gathers.
Specifically, the obtained OVT gather is processed to pick up a velocity field, the velocity field is converted into an incidence angle gather containing azimuth information, azimuth angle division superposition processing is carried out on the velocity field by the incidence angle gather according to different azimuth angle segments to obtain incidence angle gathers (called azimuth angle gather for short) in different azimuth angles, and aiming at the problems of multiple interference, low signal to noise ratio and low resolution of the azimuth angle gather, the quality of gather data is improved by methods such as high-precision radon transform, wavelet thresholding, spectral decomposition and the like.
And S20, acquiring logging data, carrying out well shock calibration and seismic horizon interpretation based on the logging data, and obtaining an interpretation horizon of the current reservoir meeting the ground stress prediction as a target horizon.
The method comprises the steps of carrying out preprocessing on logging curves in logging data, carrying out standardization processing on the logging curves, identifying logging lacking a transverse wave curve, constructing model relations between corresponding logging transverse wave speeds and preset sensitive parameters according to the logging lacking the transverse wave curve, solving transverse wave curves based on a neural network and the model relations, solving Young modulus curves and Poisson ratio curves based on the transverse wave curves of the logging curves after standardization processing, and constructing an anisotropic theoretical model based on the Young modulus curves and the Poisson ratio curves to calculate a crack density curve.
Specifically, the well logging curves in the work area, such as acoustic wave curves, density curves, natural potential curves, gamma curves, resistivity curves and the like, are subjected to standardized treatment. For a well with partial lack of a transverse wave curve, a model relation between transverse wave speed and sensitive parameters can be constructed, a transverse wave speed curve is obtained through a neural network algorithm, young modulus and Poisson ratio curves are obtained, and an anisotropic theoretical model is established to calculate a crack density curve.
Further, the method comprises the steps of carrying out well earthquake calibration and seismic horizon interpretation based on the well logging data to obtain an interpretation horizon of the current reservoir meeting the ground stress prediction, wherein the interpretation horizon is used as a target horizon, the method comprises the steps of dynamically extracting seismic wavelets based on the well logging data, continuously synthesizing seismic records based on the seismic wavelets, taking the seismic wavelets corresponding to the seismic records with the actual seismic record correlation coefficient of the well side channel smaller than a preset correlation coefficient threshold value as a calibration result, and carrying out structural interpretation on the current seismic horizon based on a preset knowledge graph based on the calibration result to obtain the interpretation horizon of the current reservoir meeting the ground stress prediction as the target horizon.
Specifically, through dynamic seismic wavelet extraction, synthetic seismic record production and drift, and comparison of the synthetic seismic record and the well side channel, when the correlation coefficient of the two reaches the preset requirement, the timely depth relationship of the seismic wavelet is the final calibration result. According to the calibration result, the target layer is defined as the seismic response characteristic, and the target layer is interpreted by combining the existing structural knowledge, so that an interpretation layer meeting the ground stress prediction is obtained, and the phenomenon of up-down jumping of the layer in the three-dimensional space is prevented.
And step S30, constructing a low-frequency model based on the logging data and the target horizon, and calculating the Poisson ratio and the crack density based on the low-frequency model.
Specifically, the Young modulus curve, the Poisson ratio curve and the crack density curve are used as longitudinal source data, the construction of a target horizon is interpreted as transverse constraint, extrapolation interpolation is carried out, and a Young modulus data body, a Poisson ratio data body and a low-frequency model data body of the crack density in a three-dimensional space are obtained and used as a low-frequency model.
In the embodiment of the invention, under the control of the time-depth curve obtained in the step S20, the Young modulus curve, the Poisson ratio curve and the crack density curve are used as source data in the longitudinal direction, and extrapolation interpolation is performed under the transverse constraint of the target horizon data, so that the Young modulus, the Poisson ratio data body and the crack density low-frequency model data body in the three-dimensional space are obtained.
Further, the Poisson's ratio and the crack density are calculated based on the low frequency model, and the Poisson's ratio and the crack density are calculated based on the low frequency model, wherein the Poisson's ratio and the crack density are calculated based on the low frequency model, a reflection coefficient matrix is obtained, a seismic reflection amplitude matrix is obtained based on a seismic wavelet matrix and a reflection coefficient matrix component pre-stack inversion model, an objective function is constructed based on a Bayesian framework and the seismic reflection amplitude matrix, and the Poisson's ratio and the crack density are calculated based on the objective function.
Specifically, the Young's modulus, poisson's ratio data volume and crack density low frequency model data volume are utilized. Calculating the seismic reflection coefficient, wherein the calculation rule is as follows:
where θ is the angle of incidence, φ is the azimuth, a is the power of density versus longitudinal wave velocity, g is the square of the transverse-to-longitudinal wave velocity ratio, where a and g are constants, E, ν, E represent Young's modulus, poisson's ratio, and fracture density, respectively.
Further, for the seismic reflection amplitudes under different incidence angles and azimuth angles, the seismic reflection coefficients are obtained by multiplying the coefficient matrix contained by vectors of parameters to be inverted, and the seismic reflection amplitudes multiplied by the wavelet matrix are used for constructing a pre-stack inversion equation:
dNMK×1＝GNMK×3Km3K×1.
Wherein d NMK×1 is an earthquake reflection amplitude matrix, G NMK×3K is an earthquake wavelet matrix, and m 3K×1 is a reflection coefficient matrix. Wherein:
Wherein N is the number of incidence angles, M is the number of azimuth angles, K is the number of sampling points, wvlt is the sub-blog matrix S (theta i,φj) is the amplitude of the seismic reflection with incidence angle theta i and azimuth angle phi j.
Further:
based on the rules, constructing an objective function under a Bayesian framework to obtain final inversion parameter estimation:
wherein σ n is the noise variance between the synthetic seismic record and the actual seismic record of the well side channel, d is the seismic reflection amplitude matrix, G is the seismic wavelet matrix, T is the sampling period, Q is the diagonal matrix, expressed as:
wherein the variance of Young's modulus, the variance of Poisson's ratio and the variance of fracture density parameters are respectively.
And S40, calculating the vertical principal stress of the target horizon based on the velocity field, the stratum density and the density inversion body.
The method comprises the steps of acquiring logging data, reading a signal speed field and stratum density based on the logging data, and calculating vertical main stress based on the speed field, the stratum density and a preset density inversion body, wherein the calculation rule is as follows:
Wherein sigma V is vertical main stress, ρ is stratum density, g w is gravitational acceleration, v f is velocity field, and t is during double journey.
And S50, calculating and obtaining the maximum horizontal principal stress of the current reservoir based on the Poisson ratio, the fracture density, the vertical principal stress and a preset principal stress calculation model.
Specifically, the preset principal stress solving model is as follows:
Wherein sigma H is the maximum principal stress, v is poisson's ratio, g is a preset variable, and e is crack density.
In the embodiment of the invention, the scheme constructs a prestack inversion equation and an objective function under a Bayesian theory framework through a newly deduced prestack azimuth AVO approximation equation, and the poisson ratio and the crack density are obtained through direct inversion of prestack anisotropy AVO from a prestack OVT gather, so that the method has the advantages of stable algorithm and high prediction precision. The scheme of the invention solves the problem that the vertical stress is not easy to be obtained in conventional ground stress calculation, utilizes the velocity body to obtain the pressure intensity generated by the overlying strata through the integral of the density inversion result, and is used as the vertical stress of the reservoir, provides a basis for the calculation of the maximum horizontal main stress, and has high calculation efficiency and is simple and easy to implement. The scheme of the invention derives a new formula for calculating the maximum horizontal principal stress based on crack density and Poisson ratio, and the new formula has the advantages of definite physical meaning, less calculation parameters and high calculation accuracy, so that the consumption of calculation storage can be reduced during large data bulk calculation, and the calculation efficiency of ground stress evaluation can be improved.
Examples:
The well earthquake calibration of a detection well is carried out, the calibration results are shown in figure 2, and the sound wave time difference curve, the density curve, the speed curve, the impedance curve, the reflection coefficient curve, the logging layering, the synthetic earthquake record (red) and the well side channel earthquake data (black) are respectively from left to right, and from shallow layers to deep layers, the synthetic earthquake record and the well side channel earthquake waveform characteristics have good similarity, and the correlation coefficient is high, so that the well earthquake calibration results are indicated to be good.
And carrying out OVT processing on the corresponding detection well set, wherein the processing result is shown in fig. 3, a phase axis at the time 1700ms in the drawing is a main objective layer waveform characteristic, and the phase axis of the gather fluctuates up and down along with traveling and azimuth angle, so that the phase axis has obvious anisotropic characteristics.
The corresponding detection well set is optimized, the processing results are shown in fig. 4, 5 and 6, fig. 4 is a azimuth-dividing gather of the OVT spiral gather after azimuth-dividing processing, fig. 5 is a azimuth-dividing gather after random noise removal of the azimuth-dividing gather, fig. 6 is a azimuth-dividing gather after the aisle-gathering leveling, the OVT gather after optimization processing has higher signal-to-noise ratio and resolution, and the AVO features of the same phase axis are obvious, so that the method can provide a data base for pre-stack seismic inversion.
As shown in fig. 7, the horizon interpretation of the shale reservoir target layer developed along the same phase axis of the earthquake from the well point is performed under the time-depth relation of well earthquake calibration, and the construction in the research area shows the characteristic of high east west.
As shown in fig. 8, a seismic section of a detection well embodying the scheme of the invention shows that the target layer of the shale reservoir is a relatively strong peak reflection characteristic. Fig. 9 and 10 are poisson's ratio and fracture density inversion sections of experimental wells obtained by pre-stack anisotropy inversion. Fig. 9 is a poisson's ratio inversion section of an experimental well. FIG. 10 is a fracture density inversion profile of an experimental well, with profile results showing lower Poisson's ratio and higher fracture density in shale reservoir sections.
Fig. 11 and 12 are poisson's ratio and fracture density inversion plan views of experimental wells obtained by pre-stack AVO anisotropic inversion. Fig. 11 is a plan view of poisson's ratio inversion in the research area, the poisson's ratio in the north of the research area is low, the poisson's ratio in the southwest and eastern areas is high, and fig. 12 is a plan view of crack density inversion in the research area, the crack in the northwest of the research area is relatively developed, and the crack in the southeast of the research area is underdeveloped.
Fig. 13 shows the result of the cross section of the maximum horizontal principal stress of the experimental well, the cross section shows that the maximum horizontal principal stress gradually increases from top to bottom, fig. 14 shows the plan view of the prediction of the maximum horizontal principal stress of the research area, the maximum horizontal principal stress of the north and south regions of the research area is higher, the maximum horizontal principal stress of the middle and south regions is lower, the prediction result has higher consistency with the known geological awareness, and the prediction result proves the effectiveness and applicability of the method.
Fig. 15 is a system structural diagram of a shale gas reservoir maximum horizontal principal stress prediction system provided by an embodiment of the invention. As shown in FIG. 15, the embodiment of the invention provides a shale gas reservoir maximum horizontal principal stress prediction system, which comprises an acquisition unit, a processing unit, a training unit and an output unit, wherein the acquisition unit is used for carrying out cross arrangement and channel extraction on an acquired original shot set to obtain a corresponding OVT channel set and carrying out preprocessing on the OVT channel set, the processing unit is used for acquiring well logging data, carrying out well seismic calibration and seismic horizon interpretation on the basis of the well logging data to obtain an interpretation horizon of the current reservoir meeting the ground stress prediction as a target horizon, the training unit is used for constructing a low-frequency model on the basis of the well logging data and the target horizon and calculating Poisson ratio and crack density on the basis of the low-frequency model, and calculating the vertical principal stress of the target horizon on the basis of the well logging data, and obtaining the maximum horizontal principal stress of the current reservoir on the basis of the Poisson ratio, the crack density, the vertical principal stress and a preset principal stress obtaining model.
The embodiment of the invention also provides a computer readable storage medium, wherein the computer readable storage medium is stored with instructions, and when the computer is run on the computer, the computer is caused to execute the shale gas reservoir maximum horizontal main stress prediction method.
Those skilled in the art will appreciate that all or part of the steps in a method for implementing the above embodiments may be implemented by a program stored in a storage medium, where the program includes several instructions for causing a single-chip microcomputer, chip or processor (processor) to perform all or part of the steps in a method according to the embodiments of the invention. The storage medium includes a U disk, a removable hard disk, a Read-Only Memory (ROM), a random access Memory (RAM, random Access Memory), a magnetic disk, an optical disk, or other various media capable of storing program codes.
The alternative embodiments of the present invention have been described in detail above with reference to the accompanying drawings, but the embodiments of the present invention are not limited to the specific details of the above embodiments, and various simple modifications may be made to the technical solutions of the embodiments of the present invention within the scope of the technical concept of the embodiments of the present invention, and all the simple modifications belong to the protection scope of the embodiments of the present invention. In addition, the specific features described in the above embodiments may be combined in any suitable manner without contradiction. In order to avoid unnecessary repetition, the various possible combinations of embodiments of the invention are not described in detail.
In addition, any combination of the various embodiments of the present invention may be made, so long as it does not deviate from the idea of the embodiments of the present invention, and it should also be regarded as what is disclosed in the embodiments of the present invention.
Claims (13)
1. A method for predicting maximum horizontal principal stress of a shale gas reservoir, the method comprising:
performing cross arrangement channel extraction on the collected original gun set to obtain a corresponding OVT channel set, and preprocessing the OVT channel set;
acquiring logging data, carrying out well shock calibration and seismic horizon interpretation based on the logging data, and obtaining an interpretation horizon of the current reservoir meeting the ground stress prediction as a target horizon;
Constructing a low-frequency model based on the logging data and the target horizon, and calculating poisson's ratio and crack density based on the low-frequency model;
calculating the vertical principal stress of the target horizon based on the velocity field, the formation density and the density inversion body;
And calculating and obtaining the maximum horizontal principal stress of the current reservoir based on the poisson ratio, the fracture density, the vertical principal stress and a preset principal stress calculation model.
2. The method of claim 1, wherein the cross-arranging the gather of the acquired raw shot gather to obtain a corresponding OVT gather comprises:
carrying out cross arrangement channel extraction on the collected original gun set to obtain a seismic channel set of the same detector under the same gun line;
and carrying out OVT unit division on each cross-arranged channel set, and carrying out channel set extraction processing on each cross-arranged channel set after division to obtain the OVT channel set.
3. The method of claim 1, wherein the preprocessing the OVT trace set comprises:
five-dimensional regularization processing is carried out on the OVT channel set, and a first OVT channel set is obtained;
Performing offset processing on the first OVT channel set to obtain a second OVT channel set;
Performing anisotropic correction on the second OVT gather to obtain a third OVT gather;
picking up a velocity field for a third OVT gather, converting the third OVT gather to an incident angle gather containing azimuth information;
Carrying out azimuth angle superposition processing on different azimuth angle section incident angle gathers to obtain azimuth angle gathers;
And performing one or more of radon transform, wavelet thresholding and spectral decomposition on the azimuth gather to complete OVT gather preprocessing.
4. The method of claim 1, wherein after collecting logging data, the method further comprises:
Performing a pretreatment of a log in log data, comprising:
Performing standardization processing on each logging curve, and identifying logging lacking a transverse wave curve, wherein,
The logging curves include:
acoustic curves, density curves, natural potential curves, gamma curves, and resistivity curves;
aiming at a well logging lacking a transverse wave curve, constructing a model relation between a corresponding well logging transverse wave speed and a preset sensitive parameter, and solving the transverse wave curve based on a neural network and the model relation;
based on the transverse wave curves of all logging curves after standardized treatment, a Young modulus curve and a Poisson ratio curve are obtained;
And constructing an anisotropic theoretical model based on the Young modulus curve and the Poisson ratio curve to calculate a crack density curve.
5. The method of claim 4, wherein the performing well-shock calibration and performing seismic horizon interpretation based on the well-logging data to obtain an interpretation horizon for which the current reservoir meets the geostress prediction as a target horizon comprises:
Dynamically extracting the seismic wavelets based on logging data, continuously synthesizing the seismic records based on the seismic wavelets, and taking the seismic wavelets corresponding to the synthesized seismic records with the correlation coefficient smaller than a preset correlation coefficient threshold value of the actual seismic records of the well side channels as calibration results;
based on the calibration result, performing construction explanation on the current seismic horizon based on a preset knowledge graph to obtain an explanation horizon of the current reservoir meeting the ground stress prediction as a target horizon.
6. The method of claim 4, wherein the constructing a low frequency model based on the log data and the target horizon comprises:
And taking the Young modulus curve, the Poisson ratio curve and the crack density curve as longitudinal source data, taking the construction of a target horizon as transverse constraint, and performing extrapolation interpolation to obtain a Young modulus data body, a Poisson ratio data body and a low-frequency model data body of the crack density in a three-dimensional space as a low-frequency model.
7. The method of claim 6, wherein the calculating poisson's ratio and fracture density based on the low frequency model comprises:
Based on the low-frequency model, calculating the seismic reflection coefficient under each incident angle and azimuth angle 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;
and constructing an objective function based on the Bayesian framework and the seismic reflection amplitude matrix, and calculating the Poisson ratio and the crack density based on the objective function.
8. The method of claim 7, wherein the seismic reflection amplitude matrix construction rule is:
dNMK×1＝GNMK×3Km3K×1.
wherein d NMK×1 is the seismic reflection amplitude matrix;
G NMK×3K is a seismic wavelet matrix;
m 3K×1 is a reflection coefficient matrix.
9. The method of claim 7, wherein the objective function is expressed as:
Wherein σ 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 a seismic wavelet matrix;
t is the transpose of the wavelet matrix;
Q is a diagonal matrix, expressed as:
wherein the variance of Young's modulus, the variance of Poisson's ratio and the variance of fracture density parameters are respectively.
10. The method of claim 9, wherein calculating a vertical principal stress of a target horizon based on the logging data comprises:
reading a signal velocity field and formation density based on the log data;
and calculating vertical main stress based on the velocity field, the stratum density and a preset density inversion body, wherein the calculation rule is as follows:
Wherein σ V is the vertical principal stress;
ρ is the formation density;
g w is gravity acceleration;
v f is the velocity field;
t is the travel time of double journey.
11. The method of claim 10, wherein the predetermined principal stress derivation model is:
Wherein σ H is the maximum principal stress;
V is poisson's ratio;
g is a preset variable;
And e is the crack density.
12. A shale gas reservoir maximum horizontal principal stress prediction system, the system comprising:
The acquisition unit is used for carrying out cross arrangement channel extraction on the acquired original gun set to obtain a corresponding OVT channel set and carrying out pretreatment on the OVT channel set;
The processing unit is used for acquiring logging data, carrying out well shock calibration and seismic horizon interpretation based on the logging data, and obtaining an interpretation horizon of the current reservoir meeting the ground stress prediction as a target horizon;
Training unit for:
Constructing a low-frequency model based on the logging data and the target horizon, and calculating poisson's ratio and crack density based on the low-frequency model;
calculating the vertical principal stress of the target horizon based on the velocity field, the formation density and the density inversion body;
and the output unit is used for calculating and obtaining the maximum horizontal principal stress of the current reservoir based on the poisson ratio, the fracture density, the vertical principal stress and a preset principal stress calculation model.
13. A computer readable storage medium having instructions stored thereon, which when run on a computer causes the computer to perform the shale gas reservoir maximum level principal stress prediction method of any of claims 1-11.
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311713818.3A CN120143231A (en) | 2023-12-13 | 2023-12-13 | Method and system for predicting maximum horizontal principal stress of shale gas reservoir |
| PCT/CN2024/139074 WO2025124529A1 (en) | 2023-12-13 | 2024-12-13 | Method and system for predicting maximum horizontal principal stress of shale gas reservoirs |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311713818.3A CN120143231A (en) | 2023-12-13 | 2023-12-13 | Method and system for predicting maximum horizontal principal stress of shale gas reservoir |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN120143231A true CN120143231A (en) | 2025-06-13 |
Family
ID=95945937
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311713818.3A Pending CN120143231A (en) | 2023-12-13 | 2023-12-13 | Method and system for predicting maximum horizontal principal stress of shale gas reservoir |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN120143231A (en) |
| WO (1) | WO2025124529A1 (en) |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2529254A4 (en) * | 2010-01-25 | 2017-08-02 | CGG Veritas Services (U.S.) Inc. | Methods and systems for estimating stress using seismic data |
| CN103258091B (en) * | 2013-04-28 | 2016-09-14 | 北京大学 | Unconventionaloil pool hides the method and device that net horizontal section three-dimensional mechanical models for rock mass is set up |
| CN111239833B (en) * | 2020-03-06 | 2022-02-18 | 中海石油(中国)有限公司 | K-value robust YPD (pre-stack simultaneous inversion) method based on Poisson ratio decomposition |
| CN113552621B (en) * | 2020-04-23 | 2024-05-28 | 中国石油天然气集团有限公司 | Shale gas ground stress determination method and device |
| CN113919196B (en) * | 2021-09-26 | 2024-09-17 | 中国石油大学(华东) | Reservoir three-dimensional stress field simulation method, simulation system, terminal, and storage medium |
-
2023
- 2023-12-13 CN CN202311713818.3A patent/CN120143231A/en active Pending
-
2024
- 2024-12-13 WO PCT/CN2024/139074 patent/WO2025124529A1/en active Pending
Also Published As
| Publication number | Publication date |
|---|---|
| WO2025124529A1 (en) | 2025-06-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113759425B (en) | Method and system for combined well-seismic evaluation of filling characteristics of deep paleokarst reservoirs | |
| CN113759424B (en) | Karst reservoir filling analysis method and system based on spectral decomposition and machine learning | |
| CN108802812B (en) | Well-seismic fusion stratum lithology inversion method | |
| CN111257926B (en) | Method for predicting ancient valley uranium reservoir by using old seismic data | |
| CN114994759A (en) | Intelligent carbon seal storage box identification method and system based on GAN network | |
| CN111983683B (en) | Prediction method and system for lake-facies limestone reservoir under low-well condition | |
| WO2016008105A1 (en) | Post-stack wave impedance inversion method based on cauchy distribution | |
| CN110737018B (en) | Method for modeling anisotropy of VSP seismic data | |
| CN110095811B (en) | Method and device for constructing and processing velocity model of paste rock stratum | |
| CN110568491B (en) | Quality factor Q estimation method | |
| CN113848593A (en) | A method for quantitative prediction of magmatic rock erosion area in coal-bearing strata | |
| CN113514884A (en) | Compact sandstone reservoir prediction method | |
| CN107703548B (en) | Shallow stratigraphic boundary division method based on sediment quality factor and echo loss grade curve peaks and valleys | |
| CN112888970B (en) | Method and device for determining slowness of sound waves | |
| CN113109875B (en) | Inversion method of carbonate rock reservoir under full waveform velocity field constraint | |
| CN112965109B (en) | Method and device for predicting reservoir under strong reflection shielding | |
| CN112835098B (en) | Method and device for predicting energy storage coefficient of weathered-crust karst reservoir | |
| CN115079256B (en) | High-resolution inversion method, device, medium and equipment for deep water turbid ponding channel gas reservoir | |
| CN118295024A (en) | Method for predicting all-stratum three-pressure body by utilizing post-stack earthquake | |
| CN120143231A (en) | Method and system for predicting maximum horizontal principal stress of shale gas reservoir | |
| US12265192B2 (en) | Reconstructing carbonate sediment transport and pathways in the geological record | |
| CN118033740A (en) | Carbonate rock oil-gas reservoir quantitative characterization method and system | |
| CN113534263A (en) | A method for predicting oil and gas saturation without relying on logging data | |
| RU2253885C1 (en) | Method for determining oil productiveness of porous collectors in three-dimensional inter-well space | |
| CN114137607B (en) | Layer sequence stratum division method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |