A seismic inversion method for tilted crack parameters and brittle indicator factors

Abstract

This invention discloses a seismic inversion method for inclined fracture parameters and brittleness indicator factors, comprising the following steps: (1) calculating Young's modulus and brittleness indicator factors; calculating two inclined fracture parameters using fracture density and fracture dip angle; performing interpolation and extrapolation based on seismic stratigraphic data and well logging data to obtain an initial model of Young's modulus, brittleness indicator factors, density, fracture density, and two inclined fracture parameters; (2) extracting azimuth angle seismic wavelets; (3) combining azimuth partial angle superimposed seismic data, azimuth angle seismic wavelets, and the initial model, obtaining Young's modulus, brittleness indicator factors, density, fracture density, and two inclined fracture parameters through Bayesian inversion; calculating the fracture dip angle based on the inverted two inclined fracture parameters. This invention can provide stable and reliable prediction results for brittleness indicator factors, fracture density, and fracture dip angle, which is helpful for conducting seismic characterization of brittleness and inclined fractures in shale gas reservoirs.

Description

Technical Field

[0001] This invention belongs to the field of geological exploration technology, and specifically relates to a seismic inversion method for tilted fracture parameters and brittleness indicator factors. Background Technology

[0002] Shale gas reservoirs are characterized by low porosity and low permeability, typically requiring horizontal drilling and hydraulic fracturing techniques to enhance reservoir recovery. The presence of natural fractures significantly improves the permeability of shale gas reservoirs. Geomechanical brittleness reflects the ability of rocks to fracture and maintain fracture open under stress conditions. Generally, rocks with higher brittleness are more easily fractured. Therefore, seismic characterization of natural fractures and rock brittleness is crucial for the exploration and development of shale gas reservoirs. Traditional methods typically equate the reservoir to a set of rotationally invariant vertical fractures embedded in an isotropic background rock-induced horizontal transverse isotropic (HTI) medium. Seismic fracture prediction is achieved through pre-stack inversion under the HTI medium assumption, which is reasonable for reservoirs with a set of vertical fractures. However, core observation data and imaging logging interpretation results indicate that fracture sets developed in reservoir rocks are more inclined than completely vertical. A set of rotationally invariant inclined fractures embedded in an isotropic background rock is equivalent to an inclined transverse isotropic (TTI) medium. In such cases, using a simple HTI medium assumption for seismic inversion may introduce prediction errors. Compared to HTI media, TTI media have an additional model parameter (i.e., fracture dip angle), making TTI medium parameter inversion more complex. Regarding brittleness prediction, previous studies have found that the ratio of Young's modulus (E) to Lamé constant (λ) (i.e., E / λ) is more sensitive to the brittleness of organic-rich, high-porosity shale gas reservoirs than Young's modulus, and therefore can be used as a brittleness indicator. Although previous studies have investigated brittleness indicator factors and fracture density prediction inversion methods, these techniques are only applicable to HTI media containing a set of vertical fractures and cannot solve the brittleness and tilt fracture prediction problems of reservoirs with inclined fractures (TTI media). Directly applying this technique to TTI media will inevitably introduce inversion errors. Summary of the Invention

[0003] The purpose of this invention is to overcome the shortcomings of the prior art and provide a seismic inversion method for tilted fracture parameters and brittleness indicator factors that can provide stable and reliable prediction results of brittleness indicator factors, fracture density and fracture dip angle, which is helpful for carrying out seismic characterization of brittleness and tilted fractures in shale gas reservoirs.

[0004] The objective of this invention is achieved through the following technical solution: a seismic inversion method for tilted crack parameters and brittleness indicator factors, comprising the following steps:

[0005] (1) Constructing an initial model: Young's modulus and brittleness indicator factor are calculated using the P-wave and S-wave velocities and densities from well logging data; two inclined fracture parameters are calculated using fracture density and fracture dip angle; and an initial model of Young's modulus, brittleness indicator factor, density, fracture density and two inclined fracture parameters is obtained by interpolation and extrapolation using seismic stratigraphic data and well logging data.

[0006] (2) Extraction of azimuth angle seismic wavelet: Perform azimuth and partial angle stacking processing on the pre-stack gather of the azimuth to obtain azimuth partial angle stacked seismic data, and extract the azimuth angle seismic wavelet by combining it with well logging data;

[0007] (3) Perform Bayesian seismic inversion: Combine the azimuth part of the angle superimposed seismic data, the azimuth angle seismic wavelet and the initial model, and obtain Young's modulus, brittleness indicator factor, density, crack density and two tilted crack parameters through Bayesian inversion; calculate the crack dip angle based on the two tilted crack parameters obtained by inversion.

[0008] The beneficial effects of this invention are as follows: This invention proposes a seismic inversion method for tilted fracture parameters and brittleness indicator factors, which can predict brittleness indicator factors, fracture density, and fracture dip angle using azimuth seismic data, and then be used for the seismic description of tilted fractures and brittleness in shale gas reservoirs. First, by combining the anisotropic Gassmann fluid substitution equation and the theoretical model of the fracture equivalent medium, the elastic stiffness of the fluid-saturated TTI medium containing fracture density and fracture dip angle is derived. Based on seismic scattering theory, the equation for the P-wave reflection coefficient of the TTI medium, expressed by Young's modulus, brittleness indicator factor, density, fracture density, and two tilted fracture parameters, is further derived. Finally, a linearized Bayesian seismic inversion method is proposed to predict tilted fracture parameters and brittleness indicator factors. Practical applications show that this invention can provide stable and reliable prediction results for brittleness indicator factors, fracture density, and fracture dip angle, which is helpful for the seismic characterization of brittleness and tilted fractures in shale gas reservoirs. Attached Figure Description

[0009] Figure 1 This is a flowchart of the seismic inversion process of the present invention;

[0010] Figure 2 This is a seismic profile diagram with azimuth angle superimposed in this embodiment;

[0011] Figure 3 This is an estimated profile of the elastic parameters and an estimated profile of the crack parameters in this embodiment. Detailed Implementation

[0012] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0013] like Figure 1As shown, the seismic inversion method for inclined crack parameters and brittleness indicator factors of the present invention includes the following steps:

[0014] (1) Constructing an initial model: Young's modulus and brittleness indicator factor are calculated using the P-wave and S-wave velocities and densities from well logging data; two inclined fracture parameters are calculated using fracture density and fracture dip angle; and an initial model of Young's modulus, brittleness indicator factor, density, fracture density and two inclined fracture parameters is obtained by interpolation and extrapolation using seismic stratigraphic data and well logging data.

[0015] The specific implementation method of step (1) is as follows: the linear slip theory simulates the crack as an imperfect bonded interface with linear slip boundary conditions, and assumes that the traction force across the interface is continuous, but the displacement is discontinuous; a set of rotationally invariant horizontal cracks are embedded in the isotropic background rock to form an equivalent vertical and transverse isotropic medium, referred to as VTI medium; the dry rock elastic stiffness matrix of VTI medium is shown below:

[0016]

[0017] In the formula, M dry =λ dry +2μ,τ dry =λ dry / M dry , λ dry and μ represent the first and second Lamé constants of isotropic background dry rocks, respectively, and M dry This represents the longitudinal wave modulus of dry rocks with an isotropic background. and Let represent the normal and tangential weaknesses of dry cracks, respectively, and their relationship with crack density be expressed as follows:

[0018]

[0019]

[0020] In the formula, k dry =μ / M dry , e represents the crack density in the coin-shaped crack model;

[0021] The anisotropic Gassmann fluid substitution theory establishes the relationship between the saturated elastic stiffness of the fluid and the elastic stiffness of dry rock, as expressed below:

[0022]

[0023] in, and Let K represent the element in the i-th row and j-th column of the fluid saturated elastic stiffness matrix and the dry rock elastic stiffness matrix, respectively; m and K f These represent the mineral bulk modulus and the fluid bulk modulus, respectively, where φ is the total porosity and α is the total porosity. i and α j All are generalized Biot coefficients;

[0024] The fluid phases in shale gas reservoirs include gas and water. Since the bulk modulus of fluids is much smaller than that of minerals, the elastic stiffness matrix of fluid-saturated VTI media with fracture density is derived by combining equations (1) to (4):

[0025]

[0026] In the formula, the values ​​of each element are calculated by substituting formulas (1) to (3) into (4), specifically:

[0027]

[0028]

[0029] f represents the fluid / pore term, and α0 represents the isotropic Biot coefficient;

[0030] The Bond transformation can be used for rotational stiffness or elastic brittleness matrices. Applying the Bond transformation to the elastic stiffness matrix of a fluid-saturated VTI medium yields the elastic stiffness matrix of a fluid-saturated TTI medium, i.e.:

[0031]

[0032] In the formula, The crack dip angle represents the angle between the crack normal and the vertical axis; the superscript T indicates transpose.

[0033] Substituting (5) into (6), we derive the elastic stiffness matrix of the fluid-saturated TTI medium containing crack density and crack inclination angle:

[0034]

[0035] In the formula,

[0036]

[0037]

[0038]

[0039]

[0040]

[0041]

[0042]

[0043]

[0044]

[0045]

[0046] When a longitudinal wave is incident on a horizontal interface between two arbitrarily anisotropic media, the relationship between the longitudinal wave reflection coefficient and the scattering function is expressed as:

[0047]

[0048] Where ρ represents density, and θ represents the angle of incidence measured relative to the vertical axis. Indicates the azimuth of the observation in the plane With the crack axis of symmetry The included angle Δ represents the difference in elastic properties between the lower and upper halves of the space (at the underground medium reflection interface); C ij Represents the element in the i-th row and j-th column of the elastic stiffness matrix of a fluid-saturated TTI medium; η ij It is a function related to the incident angle, azimuth angle, and P-wave velocity, see Shaw and Sen (Use of AVOA data to estimate fluid indicator in a vertically fractured medium, 2006);

[0049] Substituting equation (7) into equation (8), and ignoring eΔM dry 、eΔλ dry 、eΔK f The linearized longitudinal wave reflection coefficient of the TTI medium is derived from small quantities related to eΔf and eΔμ, as shown below:

[0050]

[0051] In the formula,

[0052]

[0053]

[0054] k sat =μ / M, and These are the fluid / pore term reflection coefficient, shear wave modulus reflection coefficient, and density reflection coefficient, respectively, with superscripts. - The average characteristic between the two half-spaces is represented by M, which represents the saturated rock P-wave modulus.

[0055] Previous analyses based on shale petrophysical equivalent models have shown that the ratio of Young's modulus (E) to Lamé constant (λ), E / λ, is highly sensitive to brittle shale rich in organic matter, quartz, porosity, and gas saturation. Therefore, the ratio E / λ can be used as an indicator of brittleness, with a higher E / λ value indicating better rock brittleness. The specific calculation method for E / λ will not be elaborated here. Li et al. (2022) derived the linearized P-wave reflection coefficient of HTI media containing Young's modulus, the brittleness indicator E / λ, density, and fracture density. The reflection coefficient containing the brittleness indicator E / λ can be expressed as:

[0056] Using E / λ as a brittleness indicator; the isotropic part of the reflection coefficient is expressed as:

[0057]

[0058] In the formula,

[0059]

[0060] and These represent Young's modulus reflectance and brittleness indicator reflectance, respectively;

[0061] Replacing the isotropic part in equation (9) with equation (10), we derive the linearized longitudinal wave reflection coefficient of the TTI medium, which includes Young's modulus, brittleness indicator factor E / λ, density, crack density, and crack dip angle, i.e.:

[0062]

[0063] To estimate Young's modulus, brittleness indicator E / λ, density, fracture density, and fracture dip angle from azimuth seismic data, equation (11) is further rewritten as:

[0064]

[0065] in,

[0066]

[0067] and It is related to crack density and crack dip angle, and is therefore called the inclined crack parameter.

[0068] (2) Extraction of azimuth angle seismic wavelet: Perform azimuth and partial angle stacking processing on the pre-stack gather of the azimuth to obtain azimuth partial angle stacked seismic data, and extract the azimuth angle seismic wavelet by combining it with well logging data;

[0069] (3) Perform Bayesian seismic inversion: Combine the azimuth part angle superimposed seismic data, azimuth angle seismic wavelet and initial model, and obtain Young's modulus, brittleness indicator factor, density, crack density and two tilted crack parameters through Bayesian inversion; based on the two tilted crack parameters obtained by inversion, the crack dip angle is calculated by formula (20).

[0070] The specific implementation method of step (3) is as follows: Based on the seismic convolution model, azimuth seismic data is generated by convolving the P-wave reflection coefficient with the angular seismic wavelet of each azimuth; therefore, the seismic forward model is represented in the following matrix form:

[0071] d=Gm+n (13)

[0072] Among them, G=WA, m=[R E R BI R ρ Δe Δe2 Δe4] T W is the matrix formed by the azimuth angle seismic wavelet, and A is the weighting coefficient a in equation (12). E (θ), a BI (θ), a ρ (θ), a e (θ) The resulting matrix, where n is the noise term;

[0073] Bayes' theorem is widely used to solve geophysical inversion problems. The main goal of the inversion framework is to infer unknown model parameters m based on a given set of observation data d. In implementing seismic inversion within the Bayesian framework, the goal is to evaluate the posterior model of the unknown model parameters m given the observation data d, i.e.:

[0074] P(m|d)∝P(m)P(d|m) (14)

[0075] Where P(d|m) represents the likelihood function and P(m) represents the prior model;

[0076] In seismic inversion, it is typically assumed that the likelihood function follows a zero-mean, covariance matrix of Σ. d The Gaussian distribution of is in the form of:

[0077] P(d|m)=N(d;Gm,Σ d (15)

[0078] If we assume that the prior model follows a mean of μ m The covariance matrix is ​​Σ m The Gaussian distribution, that is:

[0079] P(m) = N(m; μ) m ,Σ m (16)

[0080] According to Bayes' theorem, the posterior distribution of the model parameters is also a Gaussian distribution, that is:

[0081] P(m|d)=N(m;μ m|d ,Σ m|d (17)

[0082] The mean and covariance matrix of the posterior distribution are as follows:

[0083] μ md =μ m +Σ m G T (GΣ m G T +Σ d ) -1 (d-Gμ m (18)

[0084] Σ md =Σ m -Σ m G T (GΣ m G T +Σ d ) -1 GΣ m (19)

[0085] After obtaining Young's modulus, brittleness indicator factor, density, crack density, and two inclined crack parameters through inversion, the crack dip angle was further calculated:

[0086]

[0087] This embodiment utilizes 2D seismic survey data collected in the southern Sichuan Basin of southwestern China to verify the feasibility of the proposed method. The target layer in the study area is a black organic-rich shale gas reservoir located in the Lower Silurian Longmaxi Formation, with a high content of brittle minerals, approximately 55%–70%. Core observation and imaging logging interpretation show the presence of high-dipping fractures within the reservoir. Therefore, the TTI medium model assumption is adopted, and Bayesian linearized seismic inversion is performed to estimate the brittleness indicator factor, fracture density, and fracture dip angle. The azimuth pre-stack gathers are divided into four azimuths (20°, 65°, 110°, and 155°), and then partial incident angle stacking is performed on the pre-stack gathers of each azimuth, generating post-stack seismic data at near incident angles (10°), medium incident angles (20°), and far incident angles (30°), resulting in a total of 12 azimuth partial angle stacking seismic data sets, such as... Figure 2 As shown, the black solid line represents the P-wave impedance; (a) is the stacked seismic profile at a near incident angle (10°), (b) is the stacked seismic profile at a medium incident angle (20°), and (c) is the stacked seismic profile at a far incident angle (30°). It can be seen that the P-wave impedance in the shale gas reservoir region exhibits relatively low values, with a strong amplitude anomaly between the Longmaxi Formation and the Wufeng Formation (below the Longmaxi Formation). Figure 3 (a) and (b) show the estimated profiles of elastic parameters and fracture parameters, respectively. It can be observed that the estimated brittleness indicator and fracture density exhibit relatively high values ​​near the reservoir. Simultaneously, in the region of high fracture density, the estimated fracture dip angle exceeds 80°, indicating the presence of high-dipping fractures and relatively brittle reservoir properties. Furthermore, the inversion profiles of the brittleness indicator, fracture density, and fracture dip angle demonstrate good lateral continuity, which is helpful for the seismic characterization of rock brittleness and dipping fractures, and also facilitates the identification of areas more suitable for hydraulic fracturing.

[0088] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A seismic inversion method for tilted crack parameters and brittleness indicator factors, characterized in that, Includes the following steps: (1) Initial model construction: Young's modulus and brittleness indicator factor are calculated using the P-wave and S-wave velocities and density from well logging data; two inclined fracture parameters are calculated using fracture density and fracture dip angle; interpolation and extrapolation are performed using seismic stratigraphic data and well logging data to obtain the initial model of Young's modulus, brittleness indicator factor, density, fracture density and two inclined fracture parameters; linear slip theory simulates the fracture as an imperfect bonded interface with linear slip boundary conditions, and assumes that the traction force across the interface is continuous, but the displacement is discontinuous; a set of rotationally invariant horizontal fractures are embedded in the isotropic background rock to form an equivalent vertical and transverse isotropic medium, referred to as VTI medium; The linearized longitudinal wave reflection coefficient of TTI medium is: （12）； in, , ， ， and It is related to crack density and crack dip angle, so it is called the inclined crack parameter, where e represents the crack density in the coin-shaped crack model; , M represents the longitudinal wave modulus of saturated rock. and Let the first and second Lamé constants of the isotropic background dry rock be represented respectively. This represents the longitudinal wave modulus of dry rocks with an isotropic background. , , ; This represents the angle of incidence measured relative to the vertical axis. Indicates the azimuth of the observation in the plane With the crack axis of symmetry The included angle Δ represents the difference in elastic properties between the lower and upper half-spaces; The crack inclination angle represents the angle between the crack normal and the vertical axis. (2) Extraction of azimuth angle seismic wavelet: Perform azimuth and partial angle stacking processing on the pre-stack gather of the azimuth to obtain azimuth partial angle stacked seismic data, and extract the azimuth angle seismic wavelet by combining it with well logging data; (3) Perform Bayesian seismic inversion: Combine the azimuth part of the angle superimposed seismic data, the azimuth angle seismic wavelet and the initial model, and obtain Young's modulus, brittleness indicator factor, density, crack density and two tilted crack parameters through Bayesian inversion; calculate the crack dip angle based on the two tilted crack parameters obtained by inversion.

2. The seismic inversion method for inclined crack parameters and brittleness indicator factors according to claim 1, characterized in that, The specific implementation method of step (1) is as follows: The dry rock elastic stiffness matrix of the VTI medium is shown below: （1）； In the formula, , , , , , , , , and Let the first and second Lamé constants of the isotropic background dry rock be represented respectively. This represents the longitudinal wave modulus of dry rocks with an isotropic background. and Let represent the normal and tangential weaknesses of dry cracks, respectively, and their relationship with crack density be expressed as follows: （2）； （3）； In the formula, , , , e represents the crack density in the coin-shaped crack model; The anisotropic Gassmann fluid substitution theory establishes the relationship between the saturated elastic stiffness of the fluid and the elastic stiffness of dry rock, as expressed below: （4）； in, ; and These represent the elements in the i-th row and j-th column of the fluid saturated elastic stiffness matrix and the dry rock elastic stiffness matrix, respectively. and These represent the bulk modulus of minerals and the bulk modulus of fluids, respectively. Total porosity and All are generalized Biot coefficients; Combining equations (1) to (4), the elastic stiffness matrix of the fluid-saturated VTI medium containing crack density is derived: （5）； In the formula, , , , , , , Represents the fluid / pore term. Represents the isotropic Biot coefficient; By performing a Bond transformation on the elastic stiffness matrix of the fluid-saturated VTI medium, the elastic stiffness matrix of the fluid-saturated TTI medium is obtained, i.e.: （6）； In the formula, ; The crack dip angle represents the angle between the crack normal and the vertical axis; the superscript T indicates transpose. Substituting (5) into (6), the elastic stiffness matrix of the fluid-saturated TTI medium containing crack density and crack inclination angle is derived: （7）； In the formula, , , ， ， ， ， ， ， ， ， ， ， ； When a longitudinal wave is incident on a horizontal interface between two arbitrarily anisotropic media, the relationship between the longitudinal wave reflection coefficient and the scattering function is expressed as: （8）； in, Indicates density, This represents the angle of incidence measured relative to the vertical axis. Indicates the azimuth of the observation in the plane With the crack axis of symmetry The included angle Δ represents the difference in elastic properties between the lower and upper half-spaces; This represents the element in the i-th row and j-th column of the elastic stiffness matrix of the fluid-saturated TTI medium; It is a function related to the incident angle, azimuth angle, and P-wave velocity; Substituting equation (7) into equation (8), ignoring the... , , , and With the relevant small quantities, the linearized longitudinal wave reflection coefficient of the TTI medium is derived as follows: （9）； In the formula, ， ， , , and These are the fluid / pore term reflection coefficient, shear wave modulus reflection coefficient, and density reflection coefficient, respectively, with superscripts. The average characteristic between the two half-spaces is represented by M, which represents the saturated rock P-wave modulus. The ratio of Young's modulus E to Lamé constant λ, E / λ, is used as the brittleness indicator; the reflection coefficient containing the brittleness indicator E / λ is expressed as: （10）； In the formula, , , , and These represent Young's modulus reflectance and brittleness indicator reflectance, respectively; Replacing the isotropic part in equation (9) with equation (10), we derive the linearized longitudinal wave reflection coefficient of the TTI medium, which includes Young's modulus, brittleness indicator factor E / λ, density, crack density, and crack dip angle, i.e.: （11）。 3. The seismic inversion method for inclined crack parameters and brittleness indicator factors according to claim 1, characterized in that, The specific implementation method of step (3) is as follows: Based on the seismic convolution model, azimuth seismic data is generated by convolving the P-wave reflection coefficient with the angular seismic wavelet of each azimuth; therefore, the seismic forward model is represented in the following matrix form: （13）； in, , , It is the matrix formed by the azimuth angle seismic wavelet. It is a matrix composed of the weight coefficients in equation (12). It is a noise term; The goal of performing seismic inversion within a Bayesian framework is to derive seismic inversion from observed values. Evaluate unknown model parameters for given observation data The posterior model, namely: （14）； in, Represents the likelihood function. Represents the prior model; In seismic inversion, it is typically assumed that the likelihood function follows a zero-mean, covariance matrix of... The Gaussian distribution of is in the form of: （15）； If we assume that the prior model follows the mean... The covariance matrix is The Gaussian distribution, that is: （16）； According to Bayes' theorem, the posterior distribution of the model parameters is also a Gaussian distribution, that is: （17）； The mean and covariance matrix of the posterior distribution are as follows: （18）； （19）； After obtaining Young's modulus, brittleness indicator factor, density, crack density, and two inclined crack parameters through inversion, the crack dip angle was further calculated: （20）。

Images