Geological Engineering Dual Sweet Spot Seismic Identification Method for Shale Reservoirs

Abstract

The present invention discloses a seismic identification method for dual sweet spots in shale reservoir geological engineering, including: obtaining a new equation for the longitudinal wave reflection coefficient of VTTI medium containing fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameter; obtaining the inclined fracture density e<supgt;TTI< / supgt>, fluid bulk modulus K<subgt;f< / subgt>, brittleness index BI, pressure correlation parameter σ<subgt;E< / subgt>, and horizontal fracture density e<supgt;VTI< / supgt> through an inversion method. This application realizes a direct seismic prediction method for multiple sweet spot parameters covering shale gas fluid parameters, in-situ stress parameters, brittleness index, and horizontal and inclined fracture densities based on horizontal and inclined fracture media (VTTI medium).

Description

Technical Field

[0001] The present invention belongs to the technical field of oil and gas field development engineering, and particularly relates to a method for seismic identification of double sweet spots in shale reservoir geological engineering. Background Art

[0002] At present, for the sweet spot seismic identification technology of unconventional fractured gas-bearing shale reservoirs, the existing traditional seismic sweet spot inversion methods mainly indirectly calculate sweet spot parameters such as brittleness index and fracture density after inverting elastic parameters from seismic data. Their limitations are mainly reflected in: (1) It is necessary to use the inversion results of elastic parameters to indirectly calculate shale gas sweet spot parameters such as brittleness index. There are cumulative errors in this calculation process, resulting in inaccurate prediction of sweet spot parameters; (2) It is impossible to directly identify multiple sweet spot parameters by seismic at the same time, and there is an urgent need to improve the simultaneous prediction ability of multiple complex sweet spot parameters; (3) When predicting the sweet spot parameter of fracture density, only simple vertical fracture media or orthogonal fracture (horizontal fracture and vertical fracture) media are considered, and the situation where horizontal fractures and inclined fractures may coexist in actual fractured shale reservoirs is not considered. Summary of the Invention

[0003] The purpose of the present invention is: a method for seismic identification of double sweet spots in shale reservoir geological engineering, so as to solve the geophysical problem that it is currently impossible to directly and simultaneously predict multiple sweet spot parameters of geological engineering integration for unconventional fractured gas-bearing shale reservoirs.

[0004] The embodiment of the present application is implemented as follows. A method for seismic identification of double sweet spots in shale reservoir geological engineering is provided, including: obtaining a new equation for the longitudinal wave reflection coefficient of VTTI medium containing fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameters. The new equation is:

[0005]

[0006] Wherein,

[0007] In the formula: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient sum of the horizontal fracture density is the term coefficient of the inclined fracture density; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the perturbation of the fluid bulk modulus between the upper and lower formations, ΔBI is the perturbation of the brittleness index between the upper and lower formations, Δσ E is the perturbation of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation of the rock density between the upper and lower formations, Δe VTI is the perturbation of the horizontal fracture density between the upper and lower formations, Δe TTI is the perturbation of the inclined fracture density between the upper and lower formations; is the average value of the fluid bulk modulus between the upper and lower formations, is the average value of the brittleness index between the upper and lower formations, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton, is the square of the P-wave velocity to S-wave velocity ratio in the saturated rock;

[0008] The inclined fracture density e is obtained by the inversion method TTI 、fluid bulk modulus K f 、brittleness index BI, pressure correlation parameter σ E and horizontal fracture density e VTI .

[0009] In some embodiments, for the VTTI medium containing inclined fractures and horizontal fractures, the elastic stiffness coefficient matrix is obtained as:

[0010]

[0011] Among them, the constituent elements of the matrix in the elastic stiffness coefficient matrix are:

[0012]

[0013] In the formula, is the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, represents the stiffness coefficient at different positions in the stiffness coefficient matrix of the dry rock skeleton e VTIis the horizontal fracture density; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; e TTI is the inclined fracture density; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture; is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; δ N is the normal fracture weakness, δ T is the tangential fracture weakness.

[0014] In some embodiments, for dry or gas-saturated fractures, the fracture weakness is:

[0015]

[0016] where, is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; the superscript A is VTI or TTI. When A is VTI, e VTI is the horizontal fracture density, is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; when A is TTI, e TTI is the inclined fracture density, is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture.

[0017] In some embodiments, the calculation formula for the elastic stiffness elements of any anisotropic fluid-saturated rock is expressed as follows:

[0018]

[0019] where,

[0020]

[0021] )]]where, C sat is the stiffness coefficient of the saturated rock, C dry is the stiffness coefficient of the dry rock skeleton; φ is the total porosity of the rock; K m is the matrix bulk modulus, K f is the fluid bulk modulus; is the anisotropic Gassmann pore modulus; β and are all intermediate variable parameters with no special meaning; the subscripts I and J are the row and column positions of the stiffness coefficient elements in the stiffness coefficient matrix, and ∑ is the summation symbol.

[0022] In some embodiments, the elastic stiffness coefficient matrix and the formula of the matrix components in the elastic stiffness coefficient matrix are substituted into get:

[0023]

[0024] Where, is the P-wave modulus of the dry rock skeleton, K dry is the bulk modulus of dry rock skeleton; It is an intermediate variable parameter with no special meaning; is the normal weakness of horizontal cracks, is the normal weakness of the inclined crack;

[0025] Will Substitution get:

[0026]

[0027] Where, is the anisotropic Gassmann pore modulus; φ is the total porosity of the rock; K m is the matrix bulk modulus, K f is the bulk modulus of the fluid, K dry is the bulk modulus of dry rock skeleton; is the P-wave modulus of the dry rock skeleton; is the normal weakness of horizontal cracks, is the normal weakness of the inclined crack;

[0028] Will Substitution get:

[0029]

[0030] β4=β6=0

[0031] Where, β is the anisotropic Biot coefficient; is the P-wave modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; β0=1-K dry / K m is the isotropic rock parameter; is the normal weakness of horizontal cracks, is the normal weakness of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; ζ is the inclination angle of the inclined fracture.

[0032] In some embodiments, based on the weak anisotropic medium theory and the pore rock physics theory, it is considered that the fracture weakness parameter is close to 0 and the fluid bulk modulus is much smaller than the matrix bulk modulus, and the following is obtained:

[0033]

[0034] In the formula, is the P-wave modulus of the dry rock skeleton, K dry is the bulk modulus of the dry rock skeleton, K m is the matrix bulk modulus; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture; φ is the total porosity of the rock; K f is the fluid bulk modulus;

[0035] Based on is approximately equal to:

[0036]

[0037] In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of the rock; K f is the fluid bulk modulus;

[0038] Since φ≈φ p , therefore, according to the relationship between porosity and vertical effective stress, the quantitative relationship between porosity and vertical effective stress is obtained:

[0039]

[0040] In the formula, φ0 is the initial empirical porosity, σ V is the vertical effective stress, and is the pressure correlation parameter, denoted as σ E ; β l is the effective stress coefficient; e is the natural constant;

[0041] According to the gain function related to the bulk modulus of the dry rock skeleton, the matrix bulk modulus and the porosity, the following is obtained:

[0042] G n (φ p )=(1 - K dry / K m ) 2 / φ p =φ p / φ c 2 =σ E φ0 / φ c 2

[0043] where G n is the gain function; φ p is the background porosity; K dry is the dry rock frame bulk modulus, K m is the matrix bulk modulus; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter;

[0044] Substitute the formula

[0045]

[0046] β4 = β6 = 0 and the formula into the formula During the intermediate process of this substitution and simplification, all terms related to and are ignored, and an approximate formula for the stiffness coefficient of saturated rock is obtained:

[0047]

[0048]

[0049] where represents the stiffness coefficients at different positions in the stiffness coefficient matrix of saturated VTTI rock; K f is the fluid bulk modulus; δ ij and δ kl are Kronecker deltas; C dry is the stiffness coefficient matrix of dry VTTI medium; is the P-wave modulus in the dry rock frame, is the first Lamé coefficient in the dry rock frame; μ b is the shear modulus of saturated rock, μ b is equal to the shear modulus of the dry rock frame ζ is the dip angle of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock frame; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter.

[0050] In some embodiments, based on the approximate formula for the stiffness coefficient of saturated rock, according to the principle of weak anisotropy and weak contrast or small perturbation of the background elastic modulus at the interface, the formulas for the perturbation quantities of the corresponding stiffness coefficients in the approximate formula for the stiffness coefficient of saturated rock are respectively:

[0051]

[0052] In the formula, and are the perturbation quantities of the stiffness coefficients of saturated rock at different positions in the stiffness coefficient matrix between the upper and lower formations; is the P-wave modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b is equal to the shear modulus of the dry rock skeleton ζ is the dip angle of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter;

[0053] Based on the Born stationary phase method linearization function, that is:

[0054]

[0055] In the formula, R PP is the longitudinal wave reflection coefficient; θ is the longitudinal wave incident angle; ρ is the rock density; S(r0) is the position vector on the interface of two weakly isotropic media; Δρ represents the change in rock density between the upper and lower formations; is a known coefficient representing the relationship with the incident angle and azimuth angle; is the perturbation quantity of the stiffness coefficient at the position of the I-th row and J-th column in the saturated rock stiffness coefficient matrix, which specifically includes and

[0056] Substitute the formula for the perturbation quantity of the corresponding stiffness coefficient in the approximate formula for the stiffness coefficient of saturated rock into the formula to obtain the longitudinal wave reflection coefficient of VTTI medium as:

[0057] Among them,

[0058]

[0059] In the formula, is the longitudinal wave reflection coefficient of the VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal fracture reflection coefficient, is the anisotropic TTI inclined fracture reflection coefficient; θ is the longitudinal wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, μ b is the dry rock frame shear modulus, σ E is the pressure correlation parameter, ρ is the rock density; is the fluid bulk modulus, is the shear modulus, is the pressure correlation parameter, k ρ is the term coefficient of the rock density; ΔK f represents the change in fluid bulk modulus between the upper and lower formations; Δμ b represents the change in shear modulus between the upper and lower formations; Δσ E represents the change in pressure correlation parameter between the upper and lower formations; Δρ represents the change in rock density between the upper and lower formations; represents the change in the normal weakness of the horizontal fracture between the upper and lower formations; represents the change in the tangential weakness of the horizontal fracture between the upper and lower formations; represents the change in the normal weakness of the inclined fracture between the upper and lower formations; represents the change in the tangential weakness of the inclined fracture between the upper and lower formations; represents the average value of the fluid bulk modulus between the upper and lower formations; represents the average value of the shear modulus between the upper and lower formations; represents the average value of the pressure correlation parameter between the upper and lower formations; represents the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock frame, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack.

[0060] In some embodiments, based on the formula G n (φ p ) = (1 - K dry / K m ) 2 / φ p = φ p / φ c 2 = σ E φ0 / φ c 2 , let get:

[0061]

[0062] In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the fluid bulk modulus; μ b is the shear modulus of the saturated rock, and μ b is equal to the shear modulus of the dry rock skeleton G n is the gain function; φ p is the background porosity; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock;

[0063] For the formula take the total differential of μ b to get

[0064]

[0065] In the formula, μ b is the shear modulus of the saturated rock; Δμ b represents the change in the shear modulus between the upper and lower strata; ΔK f represents the change in the fluid bulk modulus between the upper and lower strata; Δσ E represents the change in the pressure correlation parameter between the upper and lower strata; represents the change in the parameter between the upper and lower strata; represents the average value of the shear modulus; represents the average value of the fluid bulk modulus between the upper and lower strata; represents the average value of the pressure correlation parameter between the upper and lower formations; represents between the upper and lower formations the average value of the parameter; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; K f is the fluid bulk modulus, σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock;

[0066] Introduce the brittleness index where BI is the brittleness index, and are the Young's modulus and the first Lamé coefficient of the saturated rock respectively; because this brittleness index can be obtained from expressed as: Therefore, the total differential of BI with respect to can be obtained, resulting in:

[0067]

[0068] In the formula, BI is the brittleness index; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; is the ratio of the P-wave modulus to the shear modulus in the saturated rock; ΔBI represents the change in the brittleness index between the upper and lower formations; represents between the upper and lower formations the change in the parameter; represents the average value of the brittleness index between the upper and lower formations; represents between the upper and lower formations the average value of the parameter.

[0069] Substitute formula and formula into formula Substitute formula into formula and formula to finally obtain a new equation for the P-wave reflection coefficient of the VTTI medium containing the fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameters. The new equation is:

[0070]

[0071] where,

[0072] In the formula: is the longitudinal wave reflection coefficient of the VTTI medium; θ is the incident angle of the longitudinal wave; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density and is the term coefficient of the inclined fracture density; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the perturbation of the fluid bulk modulus between the upper and lower formations, ΔBI is the perturbation of the brittleness index between the upper and lower formations, Δσ E is the perturbation of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation of the rock density between the upper and lower formations, Δe VTI is the perturbation of the horizontal fracture density between the upper and lower formations, Δe TTI is the perturbation of the inclined fracture density between the upper and lower formations; is the average value of the fluid bulk modulus between the upper and lower formations, is the average value of the brittleness index between the upper and lower formations, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock.

[0073] In some embodiments, the inversion of the inclined fracture density e TTI , includes: subtracting the seismic data of two vertically orthogonal azimuths to obtain the seismic amplitude difference Δd as the seismic data d input in the formula d = Gm + e', and then, in order to obtain the coefficient matrix A in the formula G = WAD, substituting Φ1 = 0° and Φ1 = 90° into the formula and subtracting the results to obtain:

[0074]

[0075] In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the P-wave reflection coefficient of the VTTI medium; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle;

[0076] Simplify the formula to obtain:

[0077]

[0078] In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the azimuth perturbation coefficient of the inclined fracture density; e TTI is the inclined fracture density; Δe TTI represents the perturbation quantity of the inclined fracture density between the upper and lower formations; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; represents the term coefficient of the inclined fracture density;

[0079] Invert the inclined fracture density e TTI The coefficient matrix used is expressed as:

[0080]

[0081] Among them,

[0082] In the formula, M is the number of incident angles; the subscript i is the serial number of the incident angle; N is the number of sampling points of the inclined fracture density; is the azimuth perturbation coefficient of the inclined fracture density; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; K is the coefficient matrix composed of the azimuth perturbation coefficients of the inclined fracture density under a single incident angle; the superscript "T" is the transpose operation symbol of the matrix; A1 is the coefficient matrix composed of the azimuth perturbation coefficients for inverting the inclined fracture density; the symbol "diag" is the diagonal matrix that arranges the elements in a diagonal pattern;

[0083] Take A1 in the formula as the coefficient matrix A of the parameter to be inverted in the formula G = WAD and input it. Furthermore, the parameter to be inverted is expressed as:

[0084]

[0085] In the formula, m1 is the matrix of the parameter to be inverted composed of the inclined fracture density; e TTIis the density of inclined fractures; the subscripts 1, 2, 3, …, N represent the sampling point numbers; N is the number of sampling points for the density of inclined fractures; the superscript “T” represents the transpose operation symbol of the matrix;

[0086] Finally, take m1 in the formula as the parameter m to be inverted in the formula d = Gm + e′, and then the density of inclined fractures e can be inverted TTI ; and / or

[0087] Invert the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density parameter:

[0088] Use the processed remaining seismic data d re as the seismic data d input in the formula d = Gm + e′

[0089] and then represent the coefficient matrix of the parameter to be inverted as:

[0090]

[0091] where

[0092]

[0093] X(θ i ) = diag[k BI (θ i ) k BI (θ i ) … k BI (θ i )] (N-1)×(N-1)

[0094]

[0095] U(θ i ) = diag[k ρ (θ i ·) k ρ (θ i ) … k ρ (θ i )] (N-1)×(N-1)

[0096]

[0097] In the formula, A2 is the coefficient matrix composed of the term coefficients for inverting the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density; M is the number of incident angles; the subscript i is the sequence number of the incident angle; N is the number of sampling points for each parameter; the superscript “T” is the transpose operator of the matrix; θ is the incident angle of the P-wave; is the term coefficient of the fluid bulk modulus, k BIis the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density, is the term coefficient of the inclined fracture density; the symbol "diag" is a diagonal matrix that arranges elements diagonally; F is a coefficient matrix composed of the term coefficients of the fluid bulk modulus at a single incident angle; X is a coefficient matrix composed of the term coefficients of the brittleness index at a single incident angle; V is a coefficient matrix composed of the term coefficients of the pressure correlation parameter at a single incident angle; U is a coefficient matrix composed of the term coefficients of the rock density at a single incident angle; Y is a coefficient matrix composed of the term coefficients of the inclined fracture density at a single incident angle;

[0098] Input A2 in the formula as the coefficient matrix A in the formula G = WAD, and the parameters to be inverted in this step are expressed as:

[0099]

[0100] where,

[0101]

[0102]

[0103] In the formula, is the matrix of parameters to be inverted for the fluid bulk modulus, m BI is the matrix of parameters to be inverted for the brittleness index, is the matrix of parameters to be inverted for the pressure correlation parameter, m ρ is the matrix of parameters to be inverted for the rock density, the matrix of parameters to be inverted for the horizontal fracture density;

[0104] K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; the subscripts 1, 2, 3,..., N are the sampling point numbers; the symbol "ln" is the natural logarithm operation; N is the number of sampling points for each parameter; the superscript "T" is the matrix transpose operator;

[0105] Finally, take m2 as the parameter m to be inverted in the formula d = Gm + e′, and the fluid bulk modulus K f 、the brittleness index BI, the pressure correlation parameter σ E and the horizontal fracture density e VTI can be inverted.

[0106] In some embodiments, according to Bayesian theory, the posterior probability distribution of the parameter to be inverted is expressed as follows:

[0107]

[0108] where p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, representing the fitting degree between the observed data and the model parameters; p(d) is the normalization constant; p(m|d) is the posterior probability distribution function of the parameter to be inverted; m is the parameter to be inverted;

[0109] Then, the linear forward operator is denoted by G, and the calculation formula of G is: G = WAD

[0110] where G is the linear forward operator; A is the coefficient matrix; D is the known first-order difference matrix, and W is the known wavelet matrix;

[0111] The forward model of the seismic data is given as: d = Gm + e′

[0112] where d is the input seismic data; G is the linear forward operator; m is the parameter to be inverted; e′ is the error term, following a Gaussian distribution with an expectation of 0 and a covariance matrix of C d ;

[0113] Since the forward operator G is linear and the likelihood function p(d|m) follows a Gaussian distribution, assuming that the prior distribution of the parameter m follows a Gaussian distribution, which satisfies an expectation of μ m and a prior covariance matrix C m , satisfying m ∼ N(μ m , C m ),

[0114]

[0115] where μ m is the expectation satisfied by the parameter m to be inverted; C m is the prior covariance matrix; C0 is the covariance matrix between the parameters to be inverted; the symbol is the Kronecker operation; C t is the time-domain spatial correlation matrix. Assuming that the length varying with time is t l , and the function varying with time is v(τ), then C t is expressed as:

[0116]

[0117] where C t is the time-domain spatial correlation matrix; t l is the length varying with time; v is the function varying with time;

[0118] The prior probability distribution function and the likelihood function are expressed as follows in sequence:

[0119]

[0120] where p(m) is the prior probability distribution function of the parameter to be inverted; m is the parameter to be inverted; p(d|m) represents the likelihood function; C m is the prior covariance matrix; m is the parameter to be inverted; N is the number of sampling points of the parameter to be inverted; μ m is the expectation satisfied by the parameter m to be inverted; the superscript "-1" is the symbol for matrix inversion operation; the superscript T represents the symbol for matrix transpose operation; d is the input seismic data; C m represents the prior covariance matrix satisfied by the parameter m to be inverted; C d represents the covariance matrix for which the error term e′ conforms to a Gaussian distribution; e is the natural constant; G is the linear forward operator; the symbol || represents the absolute value;

[0121] When the derivative of is zero, the obtained m is the posterior mean solution, and the analytical forms of the posterior mean and variance are given respectively as:

[0122] μ m|d = μ m +(GC m ) T (GC m G T + C d ) -1 (d - Gμ m )

[0123] C m|d = μ m -(GC m ) T (GC m G T + C d ) -1 (GC m )

[0124] where μ m|d is the posterior mean solution constrained by the seismic data, C m|d is the posterior covariance constrained by the seismic data; C m is the prior covariance matrix; the superscript "-1" represents the symbol for matrix inversion operation; d is the input seismic data; G is the linear forward operator; the superscript "T" represents the symbol for matrix transpose operation; μ m represents the expectation satisfied by the parameter m to be inverted; C d represents the covariance matrix for which the error term e′ conforms to a Gaussian distribution.

[0125] In this application:

[0126] OA: Orthogonal anisotropy; VTI: Vertical transverse isotropy; TTI: Tilted transverse isotropy; VTTI medium: A monoclinic fractured medium considering the VTI background, that is, a tilted transverse isotropy (TTI) fractured formation medium embedded in vertical transverse isotropy (VTI). Further, it can be understood that: In order to distinguish from ordinary TTI monoclinic media, the present invention refers to the monoclinic fractured medium considering the VTI background as the VTTI medium.

[0127] In summary, due to the adoption of the above technical solutions, the beneficial effects of the present invention are:

[0128] Based on the Gassmann anisotropic fluid substitution equation, this application determines the stiffness coefficient matrix of the saturated VTTI medium with high precision. In this process, by combining the fluid bulk modulus, brittleness index, vertical effective stress-related parameters, and horizontal and tilted fracture densities, a new equation for the linear PP-wave reflection coefficient is derived, providing an equation basis for directly predicting multiple sweet spot parameters in the integration of geological engineering simultaneously; based on the Bayesian inversion strategy, these multiple sweet spot parameters in the new equation are inverted from seismic data in the offset vector tile (OVT) domain. This application realizes a direct seismic prediction method for multiple sweet spot parameters covering shale gas fluid parameters, in-situ stress parameters, brittleness index, and horizontal and tilted fracture densities based on horizontal fracture and tilted fracture media (VTTI medium). Description of the Drawings

[0129] Figure 1 It is a comparison diagram of the approximate values (red dashed lines) and accurate values (blue solid lines) under different stiffness coefficients provided by the embodiments of the present invention;

[0130] Figure 2 It is another comparison diagram of the approximate values (red dashed lines) and accurate values (blue solid lines) under different stiffness coefficients provided by the embodiments of the present invention;

[0131] Figure 3 It is a well logging interpretation data diagram of vertical well A in the shale reservoir provided by the embodiments of the present invention;

[0132] Figure 4 It is a curve diagram of elastic parameters and predicted fracturing parameters of vertical well A provided by the embodiments of the present invention;

[0133] Figure 5 It is a seismic stacking data diagram of any survey line passing through well A provided by the embodiments of the present invention;

[0134] Figure 6 It is a data diagram of azimuth amplitude difference at different central incident angles provided by the embodiments of the present invention;

[0135] Figure 7 It is the inversion result diagram of the elasticity and fracture parameters of any survey line provided by the embodiments of the present invention. Specific embodiments

[0136] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0137] The technical solution of this application is as follows:

[0138] The embodiments of this application provide a method for seismic identification of double sweet spots in shale reservoir geological engineering, including:

[0139] S01. Obtain a new equation for the longitudinal wave reflection coefficient of VTTI medium containing fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density and horizontal fracture density parameters. The new equation is:

[0140]

[0141] Among them,

[0142] In the formula: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density and is the term coefficient of the inclined fracture density; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, [[ID=]56] is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the perturbation amount of the fluid bulk modulus between the upper and lower strata, ΔBI is the perturbation amount of the brittleness index between the upper and lower strata, Δσ Eis the perturbation amount of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation amount of the rock density between the upper and lower formations, Δe VTI is the perturbation amount of the horizontal fracture density between the upper and lower formations, Δe TTI is the perturbation amount of the inclined fracture density between the upper and lower formations; is the average value of the fluid bulk modulus between the upper and lower formations, is the average value of the brittleness index between the upper and lower formations, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton, is the square of the P-wave velocity to S-wave velocity ratio in the saturated rock;

[0143] S02. The inclined fracture density e is obtained by the inversion method TTI , the fluid bulk modulus K f , the brittleness index BI, the pressure correlation parameter σ E and the horizontal fracture density e VTI .

[0144] It can be understood that the double sweet spots specifically refer to the double types of sweet spots in the field of shale gas sweet spot oil and gas exploration, namely engineering sweet spots and geological sweet spots, which can be subdivided into multiple sweet spot parameters.

[0145] Based on the Gassmann anisotropic fluid substitution equation, this application determines the stiffness coefficient matrix of the saturated VTTI medium with high precision. In this process, by combining the fluid bulk modulus, brittleness index, vertical effective stress related parameters, and horizontal and inclined fracture densities, a new equation for the linear PP-wave reflection coefficient is derived, providing an equation basis for directly predicting multiple sweet spot parameters of the integration of geology and engineering simultaneously; based on the Bayesian inversion strategy, these multiple sweet spot parameters in the new equation are inverted from the seismic data in the offset vector tile (OVT) domain. This application realizes a direct seismic prediction method for multiple sweet spot parameters covering shale gas fluid parameters, in-situ stress parameters, brittleness index, horizontal and inclined fracture densities based on the horizontal fracture and inclined fracture medium (VTTI medium).

[0146] Applying the method of this application to the example of the shale fracture reservoir in the Sichuan Basin, China shows reliable effects.

[0147] In the S01:

[0148] In some embodiments, for the VTTI medium containing inclined fractures and horizontal fractures, the elastic stiffness coefficient matrix obtained is:

[0149]

[0150] Among them, the constituent elements of the matrix in the elastic stiffness coefficient matrix are:

[0151]

[0152] In the formula, is the stiffness coefficient matrix of the dry rock skeleton of VTTI medium, represents the stiffness coefficients at different positions in the stiffness coefficient matrix of the dry rock skeleton ; e VTI is the horizontal fracture density; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; e TTI is the inclined fracture density; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture; is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; δ N is the normal fracture weakness, δ T is the tangential fracture weakness.

[0153] It can be understood that is the stiffness coefficient matrix of the dry rock skeleton of VTTI medium, represents the stiffness coefficients at different positions in the stiffness coefficient matrix of the dry rock skeleton Specifically, the subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, 66 all represent the row and column positions of different stiffness coefficient components in the stiffness coefficient matrix. Taking the subscript 11 as an example, represents the stiffness coefficient in the 1st row and 1st column of the stiffness coefficient matrix of the dry rock skeleton of VTTI medium.

[0154] In some embodiments, for dry or gas-saturated fractures, the fracture weakness is:

[0155]

[0156] In the formula, is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; the superscript A is VTI or TTI. When A is VTI, e VTI is the horizontal fracture density, is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; when A is TTI, e TTI is the inclined crack density, is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack.

[0157] It can be understood that is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton.

[0158] It can be understood that for anisotropic fractured rocks, based on the low-frequency Gassmann poroelastic theory, the elastic stiffness elements of any anisotropic fluid-saturated rock are obtained, and the calculation formula for the elastic stiffness elements of any anisotropic fluid-saturated rock is:

[0159]

[0160] In the formula, C sat is the stiffness coefficient matrix of the saturated VTTI medium; K m is the matrix bulk modulus, K f is the fluid bulk modulus; δ ij and δ kl are the Kronecker symbols; C dry is the stiffness coefficient matrix of the dry VTTI medium; φ is the total porosity of the rock, and φ is equal to the sum of the fracture porosity φ f and the background porosity φ p ; the subscripts i, j, k, l, a, b, c, d all represent the order.

[0161] In some embodiments, the calculation formula for the elastic stiffness elements of any anisotropic fluid-saturated rock is expressed as follows:

[0162]

[0163] Among them,

[0164]

[0165] In the formula, C sat is the stiffness coefficient of the saturated rock, C dry is the stiffness coefficient of the dry rock skeleton; φ is the total porosity of the rock; K m is the matrix bulk modulus, K f is the fluid bulk modulus; is the anisotropic Gassmann pore modulus; β and They are all intermediate variable parameters without special meaning; the subscripts I and J are the row and column positions of the stiffness coefficient elements in the stiffness coefficient matrix, and ∑ is the summation symbol.

[0166] It can be understood that the subscripts I and J are the row and column positions of the stiffness coefficient elements in the stiffness coefficient matrix. For example, when I = 1 and J = 1, it represents the stiffness coefficient in the first row and first column of the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, and so on; ∑ is the summation symbol.

[0167] Furthermore, substitute the formula of the elastic stiffness coefficient matrix and the constituent elements of the matrix in the elastic stiffness coefficient matrix into to obtain:

[0168]

[0169] In the formula, is the P-wave modulus of the dry rock skeleton, and K dry is the bulk modulus of the dry rock skeleton; is an intermediate variable parameter without special meaning; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture;

[0170] Substitute into to obtain:

[0171]

[0172] In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of the rock; K m is the matrix bulk modulus, K f is the fluid bulk modulus, K dry is the bulk modulus of the dry rock skeleton; is the P-wave modulus of the dry rock skeleton; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture;

[0173] Substitute into to obtain:

[0174]

[0175] β4 = β6 = 0

[0176] In the formula, β is the anisotropic Biot coefficient; is the P-wave modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; β0 = 1 - Kdry / K m is the isotropic rock parameter; is the normal weakness of horizontal fractures, is the normal weakness of inclined fractures; is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton; ζ is the dip angle of the inclined fracture.

[0177] Furthermore, based on the weak anisotropy medium theory and pore rock physics theory, it is considered that the fracture weakness parameter is close to 0 and the fluid bulk modulus is much smaller than the matrix bulk modulus, and we get:

[0178]

[0179] In the formula, is the P-wave modulus of the dry rock skeleton, K dry is the bulk modulus of the dry rock skeleton, K m is the matrix bulk modulus; is the normal weakness of horizontal fractures, is the normal weakness of inclined fractures; φ is the total rock porosity; K f is the fluid bulk modulus;

[0180] Based on is approximately equal to:

[0181]

[0182] In the formula, is the anisotropic Gassmann pore modulus; φ is the total rock porosity; K f is the fluid bulk modulus;

[0183] Since φ≈φ p , so according to the relationship between porosity and vertical effective stress, the quantitative relationship between porosity and vertical effective stress is obtained:

[0184]

[0185] In the formula, φ0 is the initial empirical porosity, σ V is the vertical effective stress, and is the pressure correlation parameter, denoted as σ E ; β l is the effective stress coefficient; e is the natural constant;

[0186] It can be understood that based on the classical thin coin-shaped fracture model theory, the fracture aspect ratio is close to 0, then it is considered that the fracture porosity φ f is almost negligible compared to the background porosity φ p , so the total porosity is approximately equal to the background porosity, that is, φ≈φ p .

[0187] Based on the gain function related to the bulk modulus of the dry rock skeleton, the matrix bulk modulus, and the porosity, we obtain:

[0188] G n (φ p ) = (1 - K dry / K m ) 2 / φ p = φ p / φ c 2 = σ E φ0 / φ c 2

[0189] Where G n is the gain function; φ p is the background porosity; K dry is the bulk modulus of the dry rock skeleton, K m is the matrix bulk modulus; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter;

[0190] Substitute the formula

[0191]

[0192] β4 = β6 = 0 and the formula into the formula During the intermediate process of this substitution and simplification, ignore all terms related to and to obtain an approximate formula for the stiffness coefficient of the saturated rock:

[0193]

[0194] Where represents the stiffness coefficients at different positions in the stiffness coefficient matrix of the saturated VTTI rock; K f is the fluid bulk modulus; δ ij and δ kl are the Kronecker deltas; C dry is the stiffness coefficient matrix of the dry VTTI medium; is the P-wave modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; μ b is the shear modulus of the saturated rock, μ b is equal to the shear modulus of the dry rock skeleton ζ is the dip angle of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter.

[0195] It can be understood that represents the stiffness coefficient matrix of the saturated VTTI rock The stiffness coefficients at different positions in the matrix. The specific explanations of the subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, 66 all represent the row and column positions of different stiffness coefficient components in the stiffness coefficient matrix. Taking the subscript 11 as an example, represents the stiffness coefficient matrix of the saturated VTTI rock The stiffness coefficient in the first row and the first column of the matrix.

[0196] It can be understood that there are approximate operations in the derivation process of the approximate formula for the stiffness coefficient of the saturated rock. Therefore, in order to verify the accuracy of the derivation result in the approximate formula for the stiffness coefficient of the saturated rock, the approximate stiffness coefficient calculated by the formula:

[0197]

[0198] is compared with the accurate stiffness coefficient calculated by the formula for accuracy analysis. The matrix minerals used in the rock physics model experiment include calcite (its bulk modulus is 76.8 GPa, shear modulus is 32 GPa, density is 2.71 g / cm 3 ), quartz (its bulk modulus is 37 GPa, shear modulus is 44 GPa, density is 2.65 g / cm 3 ), and clay (bulk modulus is 21 GPa, shear modulus is 7 GPa, density is 2.6 cm 3 ). The background pores and fracture pores are filled with a mixed fluid composed of different proportions of water (its bulk modulus is 2.2 GPa, shear modulus is 0, density is 1.02 g / cm 3 ) and gas (gas bulk modulus is 0.002 GPa, shear modulus is 0, density is 0.00065 g / cm 3 ). When the gas saturation is 80%, the horizontal and inclined crack densities are 0.25φ, and when the dip angle of the inclined crack is 60° (see Figure 1) and 90° (see Figure 2 ), in the case of Figure 1 and Figure 2 it is known that from the formula:

[0199]

[0200] The approximate stiffness coefficients calculated separately are represented by red dashed lines, and the accurate stiffness coefficients calculated by the formula are represented by blue solid lines. Under different total rock porosities φ, Figure 1 (a)-(f) show that the approximate values and accurate values of different stiffness coefficients C 11 , C 12 , C 13 , C 22 , C 23 , C 33 all match very well, proving that in gas-bearing shale reservoirs with high gas content, low porosity, and high quartz content, the approximate stiffness coefficients calculated by the formula

[0201]

[0202] have high precision and meet the requirements;

[0203] Furthermore, based on the derivation of the approximate formula for the stiffness coefficient of saturated rocks, according to the principle of weak anisotropy and weak contrast or small perturbation of the interfacial background elastic modulus, the formulas for the perturbation amounts of the corresponding stiffness coefficients in the approximate formula for the stiffness coefficient of saturated rocks are respectively:

[0204]

[0205] where and are the perturbation amounts of the stiffness coefficients of saturated rocks at different positions in the stiffness coefficient matrix between the upper and lower formations; is the P-wave modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; μ b is the shear modulus of the saturated rock, μ b is equal to the shear modulus of the dry rock skeleton ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture; φ cis the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter.

[0206] It can be understood that and are the perturbation amounts of the stiffness coefficients of saturated rocks at different positions in the stiffness coefficient matrix between the upper and lower formations. Specifically, the subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, 66 all represent the row and column positions of different stiffness coefficient components in the stiffness coefficient matrix. Taking the subscript 11 as an example, represents the stiffness coefficient matrix of saturated rocks the perturbation amount of the stiffness coefficient at the position of the 1st row and 1st column between the upper and lower formations.

[0207] Furthermore, based on the Born stationary phase method linearization function, that is:

[0208]

[0209] In the formula, R PP is the longitudinal wave reflection coefficient; θ is the longitudinal wave incident angle; ρ is the rock density; S(r0) is the position vector on the interface of two weakly isotropic media; Δρ represents the change in rock density between the upper and lower formations; is a known coefficient related to the incident angle and azimuth angle; is the perturbation amount of the stiffness coefficient at the position of the Ith row and Jth column in the saturated rock stiffness coefficient matrix, which specifically includes and

[0210] Substitute the formula of the perturbation amount of the corresponding stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rocks into the formula to obtain the longitudinal wave reflection coefficient of the VTTI medium in this application as:

[0211] Among them,

[0212]

[0213] In the formula, is the longitudinal wave reflection coefficient of the VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal fracture reflection coefficient, is the anisotropic TTI inclined fracture reflection coefficient; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, μ b is the dry rock frame shear modulus, σ Eis the pressure correlation parameter, and ρ is the rock density; is the fluid bulk modulus, is the shear modulus, is the pressure correlation parameter, k ρ is the term coefficient of the rock density; ΔK f represents the change in the fluid bulk modulus between the upper and lower formations; Δμ b represents the change in the shear modulus between the upper and lower formations; Δσ E represents the change in the pressure correlation parameter between the upper and lower formations; Δρ represents the change in the rock density between the upper and lower formations; represents the change in the normal weakness of the horizontal fracture between the upper and lower formations; represents the change in the tangential weakness of the horizontal fracture between the upper and lower formations; represents the change in the normal weakness of the inclined fracture between the upper and lower formations; represents the change in the tangential weakness of the inclined fracture between the upper and lower formations; represents the average value of the fluid bulk modulus between the upper and lower formations; represents the average value of the shear modulus between the upper and lower formations; represents the average value of the pressure correlation parameter between the upper and lower formations; represents the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture.

[0214] Further, based on the formula G n (φ p ) = (1 - K dry / K m ) 2 / φ p = φ p / φ c 2 = σ E φ0 / φ c 2 , let Obtained:

[0215]

[0216] In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the fluid bulk modulus; μ b is the shear modulus of the saturated rock, μ b is equal to the shear modulus of the dry rock skeleton G n is the gain function; φ p is the background porosity; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock.

[0217] Furthermore, for the formula obtain the total differential of μ b to get

[0218] ]>

[0219] In the formula, μ b is the shear modulus of the saturated rock; Δμ b represents the change in the shear modulus between the upper and lower formations; ΔK f represents the change in the fluid bulk modulus between the upper and lower formations; Δσ E represents the change in the pressure correlation parameter between the upper and lower formations; represents the change in the parameter between the upper and lower formations; represents the average value of the shear modulus; represents the average value of the fluid bulk modulus between the upper and lower formations; represents the average value of the pressure correlation parameter between the upper and lower formations; represents the change in the parameter between the upper and lower formations; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, [[ID=6८]]is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; K f is the fluid bulk modulus, σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock;

[0220] Introduce the brittleness index where BI is the brittleness index, and are the Young's modulus and the first Lamé coefficient of the saturated rock, respectively. Since the brittleness index can be expressed by as: Therefore, the total differential of BI with respect to can be obtained as:

[0221]

[0222] In the formula, BI is the brittleness index; is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton, is the square of the P-wave velocity to S-wave velocity ratio in the saturated rock; is the ratio of the P-wave modulus to the shear modulus in the saturated rock; ΔBI represents the change in the brittleness index between the upper and lower formations; represents the change in the parameter between the upper and lower formations; represents the average value of the brittleness index between the upper and lower formations; represents the average value of the parameter between the upper and lower formations.

[0223] Substitute formula and formula into formula

[0224] , and substitute formula into formula and formula . Finally, a new equation for the P-wave reflection coefficient of the VTTI medium including the fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameters is obtained. The new equation is:

[0225]

[0226] where

[0227] In the formula: is the P-wave reflection coefficient of the VTTI medium; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, The term coefficient for the pressure correlation parameter, k ρ The term coefficient for the rock density, The term coefficient for the horizontal fracture density and The term coefficient for the inclined fracture density; The term coefficient for the normal weakness of the horizontal fracture, The term coefficient for the tangential weakness of the horizontal fracture; The term coefficient for the normal weakness of the inclined fracture, The term coefficient for the tangential weakness of the inclined fracture; ΔK f The perturbation of the fluid bulk modulus between the upper and lower formations, ΔBI is the perturbation of the brittleness index between the upper and lower formations, Δσ E The perturbation of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation of the rock density between the upper and lower formations, Δe VTI The perturbation of the horizontal fracture density between the upper and lower formations, Δe TTI The perturbation of the inclined fracture density between the upper and lower formations; The average value of the fluid bulk modulus between the upper and lower formations, The average value of the brittleness index between the upper and lower formations, The average value of the pressure correlation parameter between the upper and lower formations, The average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; The square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, The square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock.

[0228] In the S02:

[0229] In some embodiments, the inversion of the inclined fracture density e TTI , including: taking the difference of seismic data in two perpendicular azimuths to obtain the seismic amplitude difference Δd as the input of the seismic data d in the formula d = Gm + e', and then, in order to obtain the coefficient matrix A in the formula G = WAD, substituting Φ1 = 0° and Φ1 = 90° into the formula and subtracting the results to obtain:

[0230]

[0231] In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the P-wave reflection coefficient of the VTTI medium; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle;

[0232] Simplifying the formula , we get:

[0233]

[0234] In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the azimuth perturbation coefficient of the inclined fracture density; e TTI is the inclined fracture density; Δe TTI represents the perturbation amount of the inclined fracture density between the upper and lower formations; θ is the incident angle of the longitudinal wave; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; represents the term coefficient of the inclined fracture density;

[0235] Inverting the inclined fracture density e TTI The coefficient matrix used is expressed as:

[0236]

[0237] Among them,

[0238] In the formula, M is the number of incident angles; the subscript i is the serial number of the incident angle; N is the number of sampling points of the inclined fracture density; is the azimuth perturbation coefficient of the inclined fracture density; θ is the incident angle of the longitudinal wave; ζ is the dip angle of the inclined fracture; K is the coefficient matrix composed of the azimuth perturbation coefficients of the inclined fracture density under a single incident angle; the superscript "T" is the transpose operation symbol of the matrix; A1 is the coefficient matrix composed of the azimuth perturbation coefficients for inverting the inclined fracture density; the symbol "diag" is the diagonal matrix that arranges the elements in a diagonal pattern;

[0239] Substituting A1 in the formula as the coefficient matrix A of the parameter to be inverted in the formula G = WAD, and furthermore, the parameter to be inverted is expressed as:

[0240]

[0241] In the formula, m1 is the matrix of the parameter to be inverted composed of the inclined fracture density; e TTI is the inclined fracture density; the subscripts 1, 2, 3,..., N represent the sampling point serial numbers; N is the number of sampling points of the inclined fracture density; the superscript "T" represents the transpose operation symbol of the matrix;

[0242] Finally, substituting m1 in the formula as the parameter m to be inverted in the formula d = Gm + e′, the inclined fracture density e TTI can be inverted.

[0243] Furthermore, according to the Bayesian theory, the posterior probability distribution of the parameter to be inverted is expressed as follows:

[0244]

[0245] Wherein, p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, representing the fitting degree between the observed data and the model parameters; p(d) is the normalization constant; p(m|d) is the posterior probability distribution function of the parameter to be inverted; m is the parameter to be inverted;

[0246] Then, the linear forward operator is represented by G, and the calculation formula of G is: G = WAD

[0247] Wherein, G is the linear forward operator; A is the coefficient matrix; D is the known first-order difference matrix, and W is the known wavelet matrix;

[0248] Furthermore, the forward model of the seismic data is given as: d = Gm + e′

[0249] Wherein, d is the input seismic data; G is the linear forward operator; m is the parameter to be inverted; e′ is the error term, which conforms to the Gaussian distribution with an expectation of 0 and a covariance matrix of C d of.

[0250] It can be understood that:

[0251] Since the forward operator G is linear and the likelihood function p(d|m) follows the Gaussian distribution, assuming that the prior distribution of the parameter m follows the Gaussian distribution, which satisfies the expectation of μ m and the prior covariance matrix C m , satisfying m ~ N(μ m , C m ),

[0252]

[0253] Wherein, μ m is the expectation satisfied by the parameter m to be inverted; C m is the prior covariance matrix; C0 is the covariance matrix between the parameters to be inverted; the symbol is the Kronecker operation; C t is the time-domain spatial correlation matrix. Assuming that the length varying with time is t l , and the function varying with time is v(τ), then C t is expressed as:

[0254]

[0255] Wherein, C t is the time-domain spatial correlation matrix; t l is the length varying with time; v is the function varying with time;

[0256] The prior probability distribution function and the likelihood function are respectively expressed as follows:

[0257]

[0258] Wherein, p(m) is the prior probability distribution function of the parameter to be inverted; m is the parameter to be inverted; p(d|m) represents the likelihood function; C m is the prior covariance matrix; m is the parameter to be inverted; N is the number of sampling points of the parameter to be inverted; μ m is the expectation satisfied by the parameter m to be inverted; the superscript "-1" is the symbol for matrix inversion operation; the superscript T represents the transpose operation symbol of the matrix; d is the input seismic data; C m represents the prior covariance matrix satisfied by the parameter m to be inverted; C d represents the covariance matrix of the error term e' conforming to the Gaussian distribution; e is the natural constant; G is the linear forward operator; the symbol || represents the absolute value;

[0259] When the derivative of is zero, the obtained m is the posterior mean solution, and the analytical forms of the posterior mean and variance are respectively given as:

[0260] μ m|d = μ m +(GC m ) T (GC m G T + C d ) -1 (d - Gμ m )

[0261] C m|d = μ m -(GC m ) T (GC m G T + C d ) -1 (GC m )

[0262] Wherein, μ m|d is the posterior mean solution constrained by the seismic data, C m|d is the posterior covariance constrained by the seismic data; C m is the prior covariance matrix; the superscript "-1" represents the matrix inversion operation symbol; d is the input seismic data; G is the linear forward operator; the superscript "T" represents the transpose operation symbol of the matrix; μ m represents the expectation satisfied by the parameter m to be inverted; C d represents the covariance matrix of the error term e' conforming to the Gaussian distribution.

[0263] In some embodiments, the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density parameter are inverted:

[0264] Using the remaining processed seismic data d re As the seismic data d in the formula d = Gm + e′

[0265] is input, and then the coefficient matrix of the parameters to be inverted is expressed as:

[0266]

[0267] Wherein,

[0268]

[0269] X(θ i ) = diag[k BI (θ i ) k BI (θ i ) … k BI (θ i )] (N-1)×(N-1)

[0270]

[0271] U(θ i ) = diag[k ρ (θ i ) k ρ (θ i ) … k ρ (θ i )] (N-1)×(N-1)

[0272]

[0273] In the formula, A2 is the coefficient matrix composed of the term coefficients for inverting the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density; M is the number of incident angles; the subscript i is the sequence number of the incident angle; N is the number of sampling points for each parameter; the superscript "T" is the matrix transpose operator; θ is the P-wave incident angle; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density, The term coefficient for the inclined fracture density; the symbol "diag" is a diagonal matrix that arranges elements diagonally; F is a coefficient matrix composed of the term coefficients of the fluid bulk modulus at a single incident angle; X is a coefficient matrix composed of the term coefficients of the brittleness index at a single incident angle; V is a coefficient matrix composed of the term coefficients of the pressure correlation parameter at a single incident angle; U is a coefficient matrix composed of the term coefficients of the rock density at a single incident angle; Y is a coefficient matrix composed of the term coefficients of the inclined fracture density at a single incident angle;

[0274] Substitute A2 in the formula as the coefficient matrix A in the formula G = WAD, and the parameters to be inverted in this step are expressed as:

[0275]

[0276] where,

[0277]

[0278]

[0279] In the formula, is the matrix of parameters to be inverted for the fluid bulk modulus, m BI is the matrix of parameters to be inverted for the brittleness index, is the matrix of parameters to be inverted for the pressure correlation parameter, m ρ is the matrix of parameters to be inverted for the rock density, the matrix of parameters to be inverted for the horizontal fracture density;

[0280] K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; the subscripts 1, 2, 3,..., N are the sampling point numbers; the symbol "ln" is the natural logarithm operation; N is the number of sampling points for each parameter; the superscript "T" is the matrix transpose operator;

[0281] Finally, take m2 as the parameter m to be inverted in the formula d = Gm + e′, and then the fluid bulk modulus K f 、the brittleness index BI, the pressure correlation parameter σ E and the horizontal fracture density e VTI can be inverted.

[0282] Application Example

[0283] The actual working area is located in the Sichuan Basin, China, and is an unconventional fractured shale gas reservoir in the Longmaxi-Wufeng formation. A vertical well, Well A, has been drilled on this fracturing platform, which can provide imaging logging information (see Figure 3 ), and conventional logging elastic parameters and fracture parameter curves (see Figure 4 ). Figure 3 In Figure 3 , (a) FMI imaging, (b) core photos, (c) dip fracture azimuth rose diagram, and (d) fracture dip statistical histogram. It can be seen from Figure 4 that the logging imaging data of vertical Well A shows that natural fractures in the reservoir are developed, including bedding horizontal fractures and high-angle inclined fractures. The strike of the fractures is mainly southwest-northeast, and the dip angle of the inclined fractures mainly concentrates on 70°. Therefore, it is appropriate to regard the gas-bearing shale reservoir in this working area as a VTTI medium to verify the effectiveness of the proposed method. Figure 4 In

[0284] , azimuthal seismic data are first extracted from the OVT domain, along the strike of natural fractures and its perpendicular direction respectively, that is, the azimuth angles are 45° and 135°. Then, with 3°, 9°, 15°, 21°, and 27° as the central incident angles, the azimuthal seismic data are stacked to improve the signal-to-noise ratio of the seismic data. Finally, the stacked seismic data of partial central incident angles at these two azimuth angles are obtained, see Figure 5 . Figure 5 contains different azimuth angles and central incident angles, where (a)-(e) the observation azimuth Φ abs1 is 45°, (f)-(j) the observation azimuth Φ abs2 is 135°, and is divided into different central incident angles: the central incident angle θ1 of (a) and (f) is 3°, the central incident angle θ2 of (b) and (g) is 9°, the central incident angle θ3 of (c) and (h) is 15°, the central incident angle θ3 of (d) and (i) is 21°; the central incident angle θ3 of (e) and (j) is 27°. It can be seen fromFigure 5 It can be seen that the azimuthal seismic data used in this example application is divided into two azimuth angles: 45 degrees and 135 degrees, and the two azimuths are stacked respectively according to the central incident angles of 3°, 9°, 15°, 21°, and 27° to obtain seismic stacked data with different azimuths and different central incident angles, which serves as the basis for the subsequent inversion seismic data. Among the five central incident angles, we preferentially select three relatively large central incident angles, namely 15°, 21°, and 27°, and generate azimuth amplitude difference seismic data between the above two azimuth angles, as shown in Figure 6 , Figure 6 In (a), the central incident angle θ3 is 15°, in (b), the central incident angle θ3 is 21°, and in (c), the central incident angle θ3 is 27°. From Figure 6 it can be seen that by subtracting the seismic profile data of different azimuths with central incident angles of 15°, 21° Figure 5 and 27° under two different azimuths in

[0285] , the amplitude difference seismic data of azimuthal seismic with different central incident angles can be obtained. It can be seen at the stratification position of Well T3 that the amplitude difference seismic data shows obvious anomalies, which is closely related to the development of natural fractures in the reservoir. To invert the anisotropic inclined fracture density. In Figure 6 , at the reservoir position of the well stratification T2 - T3, obvious amplitude differences can be observed between the two azimuth angles, which is related to the development of natural fractures in the reservoir. Based on the multi - sweet - spot parameter seismic inversion method based on the Bayesian inversion strategy introduced in the present invention, the inversion result of this survey line corresponding to the seismic data is obtained (see Figure 7 ).

[0286] From Figure 7 it can be seen that Figure 7 the fluid bulk modulus in (a) shows a low - value anomaly, indicating a relatively high shale gas content; Figure 7 ] the brittleness index in (b) shows a high - value anomaly, indicating good brittleness of the shale rock; Figure 7 the stress correlation parameter in (c) shows a high - value anomaly, which must indicate that due to the relatively high shale gas content and large pore pressure, the vertical effective pressure is relatively small; Figure 7 the shale gas density in (d) shows a high - value anomaly. Figure 7 The horizontal fracture density in (e) and Figure 7 the inclined fracture density in (f) show high - value anomalies, indicating the development of horizontal and inclined natural fractures in the reservoir. Therefore, this reservoir has the characteristics of low fluid bulk modulus, high brittleness index, high effective stress - related parameters, and high horizontal and inclined fracture densities, and the prediction is consistent with the actual geological background and objective understanding.

[0287] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for seismic identification of double sweet spots in shale reservoir geological engineering, characterized in that, Including: Obtaining a new equation for the longitudinal wave reflection coefficient of a VTTI medium including fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameter, and the new equation is: Wherein, In the formula: is the longitudinal wave reflection coefficient of the VTTI medium; θ is the incident angle of the longitudinal wave; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density and is the term coefficient of the inclined fracture density; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the perturbation of the fluid bulk modulus between the upper and lower formations, ΔBI is the perturbation of the brittleness index between the upper and lower formations, Δσ E is the perturbation of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation of the rock density between the upper and lower formations, Δe VTI is the perturbation of the horizontal fracture density between the upper and lower formations, Δe TTI is the perturbation of the inclined fracture density between the upper and lower formations; is the average value of the fluid bulk modulus between the upper and lower formations, is the average value of the brittleness index between the upper and lower formations, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; Obtain the inclined fracture density e through the inversion method TTI , fluid bulk modulus K f , brittleness index BI, pressure correlation parameter σ E and horizontal fracture density e VTI .

2. The method for seismic identification of dual sweet spots in shale reservoir geological engineering according to claim 1, wherein For a VTTI medium containing inclined fractures and horizontal fractures, the elastic stiffness coefficient matrix is obtained as: Wherein, the constituent elements of the matrix in the elastic stiffness coefficient matrix are: In the formula, is the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, represents the stiffness coefficients at different positions in the stiffness coefficient matrix of the dry rock skeleton ; e VTI is the horizontal fracture density; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; e TTI is the inclined fracture density; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture; is the P-wave modulus in the dry rock skeleton, μ b is the shear modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; δ N is the normal fracture weakness, δ T is the tangential fracture weakness.

3. The method for seismic identification of double sweet spots in shale reservoir geological engineering according to claim 2, characterized in that, For dry or gas-saturated fractures, the fracture weakness is: In the formula, is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton; the superscript A′ is VTI or TTI. When A′ is VTI, e VTI is the horizontal fracture density, is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; when A′ is TTI, e TTI is the inclined fracture density, is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture.

4. The method for seismic identification of double sweet spots in shale reservoir geological engineering according to claim 2, wherein The calculation formula for the elastic stiffness elements of any anisotropic fluid-saturated rock is expressed as follows: Wherein, where C sat is the stiffness coefficient of saturated rock, and C dry is the stiffness coefficient of dry rock skeleton; φ is the total porosity of rock; Km is the matrix bulk modulus, and K f is the fluid bulk modulus; is the anisotropic Gassmann pore modulus; β is the anisotropic Biot coefficient; is an intermediate variable parameter with no special meaning; the subscripts I and J are the row and column positions of the stiffness coefficient elements in the stiffness coefficient matrix, and ∑ is the summation symbol.

5. The method for seismic identification of double sweet spots in shale reservoir geological engineering according to claim 4, wherein Substitute the elastic stiffness coefficient matrix and the formula for the elements that make up the matrix in the elastic stiffness coefficient matrix into to obtain: In the formula, is the P-wave modulus of the dry rock skeleton, and K dry is the bulk modulus of the dry rock skeleton; is an intermediate variable parameter with no special meaning; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture; Substitute into to obtain: wherein, is the anisotropic Gassmann pore modulus; φ is the total porosity of the rock; K m is the matrix bulk modulus, K f is the fluid bulk modulus, K dry is the dry rock frame bulk modulus; is the P-wave modulus of the dry rock frame; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture; Substitute into to obtain: β4＝β6＝0 In the formula, β is the anisotropic Biot coefficient; is the P-wave modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; β0 = 1 - K dry / K m is the isotropic rock parameter; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; ζ is the dip angle of the inclined fracture.

6. The method for seismic identification of double sweet spots in shale reservoir geological engineering according to claim 5, characterized in that Based on the weak anisotropy medium theory and pore rock physics theory, it is considered that the fracture weakness parameter is close to 0 and the fluid bulk modulus is much smaller than the matrix bulk modulus, and it is obtained that: In the formula, is the P-wave modulus of the dry rock skeleton, K dry is the bulk modulus of the dry rock skeleton, K m is the matrix bulk modulus; is the normal weakness of the horizontal fracture, is the normal weakness of the inclined fracture; φ is the total porosity of the rock; K f is the fluid bulk modulus; Based on Approximately equal to: In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of the rock; K f is the fluid bulk modulus; Since φ≈φ p , a quantitative relationship between porosity and vertical effective stress is obtained according to the relationship between porosity and vertical effective stress: where φ0 is the initial empirical porosity, and σ V is the vertical effective stress, and is denoted as σ E ; β l is the effective stress coefficient; e is the natural constant; According to the gain function related to the bulk modulus of the dry rock skeleton, matrix bulk modulus, and porosity, it is obtained that: G n (φ p ) = (1 - K dry / K m ) 2 / φ p = φ p / φ c 2 = σ E φ0 / φ c 2 where G n is the gain function; φ p is the background porosity; K dry is the dry rock frame bulk modulus, K m is the matrix bulk modulus; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter; The formula β4 = β6 = 0 and the formula Substitute into the formula During the intermediate process of this substitution and simplification, all terms related to and are ignored, and an approximate formula for the stiffness coefficient of saturated rock is obtained: In the formula, represents the stiffness coefficient at different positions in the stiffness coefficient matrix of saturated VTTI rock ; K f is the fluid bulk modulus; C dry is the stiffness coefficient of the dry rock skeleton; is the P-wave modulus in the dry rock skeleton, is the first Lamé coefficient in the dry rock skeleton; μ b is the shear modulus in the dry rock skeleton; ζ is the dip angle of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter.

7. The method for seismic identification of double sweet spots in shale reservoir geological engineering according to claim 6, wherein Based on the derived approximate formula for the stiffness coefficient of the saturated rock, according to the principle of weak anisotropy and weak contrast or small perturbation of the background elastic modulus at the interface, the formulas for the perturbation amounts of the corresponding stiffness coefficients in the approximate formula for the stiffness coefficient of the saturated rock are respectively: In the formula, and are the perturbation amounts of the stiffness coefficients of saturated rocks at different positions in the stiffness coefficient matrix between the upper and lower formations; is the P-wave modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; μ b is the shear modulus in the dry rock skeleton; ζ is the dip angle of the inclined crack; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter; Based on the Born stationary phase method linearization function, that is: Wherein, R PP is the longitudinal wave reflection coefficient; θ is the incident angle of the longitudinal wave; ρ is the rock density; Δρ represents the perturbation of the rock density between the upper and lower formations; is a known coefficient related to the incident angle and azimuth angle; is the perturbation of the stiffness coefficient at the position of the I-th row and J-th column in the stiffness coefficient matrix of the saturated rock, which specifically includes and Substitute the formula of the perturbation of the corresponding stiffness coefficient in the approximate formula of the stiffness coefficient of the saturated rock into the formula to obtain the longitudinal wave reflection coefficient of the VTTI medium as: Wherein, Wherein, is the longitudinal wave reflection coefficient of the VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal fracture reflection coefficient, is the anisotropic TTI inclined fracture reflection coefficient; θ is the longitudinal wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, μ b is the dry rock frame shear modulus, σ E is the pressure correlation parameter, ρ is the rock density; is the term coefficient of the fluid bulk modulus, is the term coefficient of the shear modulus, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density; ΔK f represents the perturbation of the fluid bulk modulus between the upper and lower formations; Δμ b represents the change in the shear modulus between the upper and lower formations; Δσ E represents the perturbation of the pressure correlation parameter between the upper and lower formations; Δρ represents the perturbation of the rock density between the upper and lower formations; represents the change in the normal weakness of the horizontal fracture between the upper and lower formations; represents the change in the tangential weakness of the horizontal fracture between the upper and lower formations; represents the change in the normal weakness of the inclined fracture between the upper and lower formations; represents the change in the tangential weakness of the inclined fracture between the upper and lower formations; represents the average value of the fluid bulk modulus between the upper and lower formations; represents the average value of the shear modulus between the upper and lower formations; represents the average value of the pressure correlation parameter between the upper and lower formations; represents the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock frame, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the normal weakness of the horizontal fracture, is the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; is the normal weakness of the inclined fracture, is the tangential weakness of the inclined fracture.

8. The method for seismic identification of dual sweet spots in shale reservoir geological engineering according to claim 7, characterized in that, Based on formula G n (φ p ) = (1 - K dry / K m ) 2 / φ p = φ p / φ c 2 = σ E φ0 / φ c 2 , assume Get: In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the fluid bulk modulus; μ b is the shear modulus in the dry rock skeleton; G n is the gain function; φ p is the background porosity; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock; For the formula Find μ b The total differential of, we get where μ b is the shear modulus in the dry rock skeleton; Δμ b represents the change in shear modulus between the upper and lower formations; ΔK f represents the perturbation of the fluid bulk modulus between the upper and lower formations; Δσ E represents the perturbation of the pressure correlation parameter between the upper and lower formations; represents between the upper and lower formations the change in the parameter; represents the average value of the shear modulus; represents the average value of the fluid bulk modulus between the upper and lower formations; represents the average value of the pressure correlation parameter between the upper and lower formations; represents between the upper and lower formations the average value of the parameter; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock; K f is the fluid bulk modulus, σ E is the pressure correlation parameter; is the ratio of the P-wave modulus to the shear modulus in the saturated rock; Introduction of brittleness index where BI is the brittleness index, and are the Young's modulus and the first Lame coefficient of the saturated rock, respectively; since this brittleness index can be represented by as: Therefore, the total differential of BI with respect to can be obtained as: Where BI is the brittleness index; is the square of the P-wave velocity to S-wave velocity ratio in the dry rock skeleton, is the square of the P-wave velocity to S-wave velocity ratio in the saturated rock; is the ratio of the P-wave modulus to the shear modulus in the saturated rock; ΔBI represents the perturbation of the brittleness index between the upper and lower formations; represents between the upper and lower formations the change in the parameter; represents the average value of the brittleness index between the upper and lower formations; represents between the upper and lower formations the average value of the parameter; Substitute the formula and the formula into the formula Substitute the formula into the formula and the formula to finally obtain a new equation for the P-wave reflection coefficient of the VTTI medium that includes the fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density, and horizontal fracture density parameter. The new equation is as follows: Wherein, In the formula: is the longitudinal wave reflection coefficient of the VTTI medium; θ is the longitudinal wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density and is the term coefficient of the inclined fracture density; is the term coefficient of the normal weakness of the horizontal fracture, is the term coefficient of the tangential weakness of the horizontal fracture; is the term coefficient of the normal weakness of the inclined fracture, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the perturbation of the fluid bulk modulus between the upper and lower formations, ΔBI is the perturbation of the brittleness index between the upper and lower formations, Δσ E is the perturbation of the pressure correlation parameter between the upper and lower formations, Δρ is the perturbation of the rock density between the upper and lower formations, Δe VTI is the perturbation of the horizontal fracture density between the upper and lower formations, Δe TTI is the perturbation of the inclined fracture density between the upper and lower formations; is the average value of the fluid bulk modulus between the upper and lower formations, is the average value of the brittleness index between the upper and lower formations, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of the rock density between the upper and lower formations; ζ is the dip angle of the inclined fracture; is the square of the ratio of the P-wave velocity to the S-wave velocity in the dry rock skeleton, is the square of the ratio of the P-wave velocity to the S-wave velocity in the saturated rock.

9. The method for seismic identification of dual sweet spots in shale reservoir geological engineering according to claim 1, wherein Inversion of inclined fracture density e TTI , including: taking the difference between seismic data in two perpendicular orthogonal azimuths to obtain the seismic amplitude difference Δd as the input of seismic data d in the formula d = Gm + e′, and then, in order to obtain the coefficient matrix A in the formula G = WAD, substituting Φ1 = 0° and Φ1 = 90° into the formula , subtracting the results to obtain: In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the P-wave reflection coefficient of the VTTI medium; θ is the P-wave incident angle; ζ is the dip angle of the inclined fracture; Φ is the azimuth angle; Simplify the formula to obtain: In the formula, is the azimuth perturbation reflection coefficient caused by the TTI inclined fracture; is the azimuth perturbation coefficient of the inclined fracture density; e TTI is the inclined fracture density; Δe TTI represents the perturbation amount of the inclined fracture density between the upper and lower formations; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined fracture; Φ is the azimuth angle; represents the term coefficient of the inclined fracture density; Inversion of the inclined fracture density e TTI The coefficient matrix used is expressed as: Among them, where M is the number of incident angles; the subscript i is the serial number of the incident angle; is the azimuth perturbation coefficient of the inclined fracture density; θ is the incident angle of the P-wave; ζ is the dip angle of the inclined fracture; K is the coefficient matrix composed of the azimuth perturbation coefficients of the inclined fracture density under a single incident angle; the superscript "T" is the transpose operation symbol of the matrix; A1 is the coefficient matrix composed of the azimuth perturbation coefficients for inverting the inclined fracture density; the symbol "diag" is a diagonal matrix that arranges the elements in a diagonal pattern; Input A1 in the formula as the coefficient matrix A of the parameter to be inverted in the formula G = WAD. Furthermore, the parameter to be inverted is expressed as: where m1 is the matrix of parameters to be inverted composed of the inclined crack density; e TTI is the inclined crack density; the subscripts 1, 2, 3,..., N represent the sampling point numbers; the superscript "T" represents the transpose operation symbol of the matrix; Finally, taking m1 in the formula as the parameter m to be inverted in the formula d = Gm + e′, the inclined crack density e can be inverted TTI ; and / or Inverting the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density parameter: Using the processed residual seismic data d re as the formula d = Gm + e′ Input the seismic data d in, and then express the coefficient matrix of the parameter to be inverted as: Wherein, In the formula, A2 is a coefficient matrix composed of term coefficients for inverting the fluid bulk modulus, brittleness index, pressure correlation parameter, and horizontal fracture density; M is the number of incident angles; the subscript i is the sequence number of the incident angle; N is the number of sampling points for each parameter; the superscript "T" is the transpose operator of the matrix; θ is the P-wave incident angle; is the term coefficient of the fluid bulk modulus, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure correlation parameter, k ρ is the term coefficient of the rock density, is the term coefficient of the horizontal fracture density, is the term coefficient of the inclined fracture density; the symbol "diag" is a diagonal matrix that arranges elements in a diagonal pattern; F is a coefficient matrix composed of term coefficients of the fluid bulk modulus under a single incident angle; X is a coefficient matrix composed of term coefficients of the brittleness index under a single incident angle; V is a coefficient matrix composed of term coefficients of the pressure correlation parameter under a single incident angle; U is a coefficient matrix composed of term coefficients of the rock density under a single incident angle; Y is a coefficient matrix composed of term coefficients of the inclined fracture density under a single incident angle; Input A2 in the formula as the coefficient matrix A in the formula G = WAD, and the parameters to be inverted in this step are expressed as: Wherein, In the formula, is the matrix of parameters to be inverted for the fluid bulk modulus, m BI is the matrix of parameters to be inverted for the brittleness index, is the matrix of parameters to be inverted for the pressure correlation parameter, m ρ is the matrix of parameters to be inverted for the rock density, is the matrix of parameters to be inverted for the horizontal fracture density; K f is the fluid bulk modulus, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal fracture density, e TTI is the inclined fracture density; the subscripts 1, 2, 3, ..., N are the sampling point numbers; the symbol "ln” represents the natural logarithm operation; N is the number of sampling points for each parameter; the superscript "T” is the matrix transpose operator; Finally, taking m2 as the parameter m to be inverted in the formula d = Gm + e′, the bulk modulus K of the fluid can be inverted f , brittleness index BI, pressure correlation parameter σ E and horizontal fracture density e VTI .

10. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 1, wherein According to the Bayesian theory, the posterior probability distribution of the parameter to be inverted is expressed as follows: In the formula, p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, representing the fitting degree between the observed data and the model parameters; p(d) is the normalization constant; p(m|d) is the posterior probability distribution function of the parameter to be inverted; m is the parameter to be inverted; Then, the linear forward operator is represented by G, and the calculation formula of G is: G = WAD In the formula, G is the linear forward operator; A is the coefficient matrix; D is the known first-order difference matrix, and W is the known wavelet matrix; The forward model of the seismic data is given as: d = Gm + e′ where d is the input seismic data; G is the linear forward operator; m is the parameter to be inverted; e′ is the error term, which follows a Gaussian distribution with an expectation of 0 and a covariance matrix of C d ; Since the forward operator G is linear and the likelihood function p(d|m) follows a Gaussian distribution, assuming that the prior distribution of the parameter m follows a Gaussian distribution with an expectation of μ m and a prior covariance matrix C m , such that m ~ N(μ m , C m ), where μ m is the expectation that the inversion parameter m satisfies; C m is the prior covariance matrix; C0 is the covariance matrix between the inversion parameters; the symbol is the Kronecker operation; C t is the time-domain spatial correlation matrix, assuming that the length varying with time is t l , and the function varying with time is v(τ), then C t is expressed as: where C t is the time-domain spatial correlation matrix; t l is the length varying with time; v is the function varying with time; The prior probability distribution function and the likelihood function are respectively expressed as follows: where p(m) is the prior probability distribution function of the parameter to be inverted; m is the parameter to be inverted; p(d|m) represents the likelihood function; C m is the prior covariance matrix; m is the parameter to be inverted; N is the number of sampling points of the parameter to be inverted; μ m is the expectation that the parameter m to be inverted satisfies; the superscript "-1” is the symbol for matrix inversion operation; the superscript T represents the symbol for matrix transpose operation; d is the input seismic data; C m represents the prior covariance matrix that the parameter m to be inverted satisfies; C d represents the covariance matrix that the error term e′ conforms to the Gaussian distribution; e is the natural constant; G Is the linear forward operator; the symbol || represents the absolute value; When has a derivative of zero, the resulting m is the posterior mean solution, and the analytical forms of the posterior mean and variance are given respectively as: μ m|d = μ m + (GC m ) T (GC m G T + C d ) -1 (d - Gμ m ) C m|d = μ m -(GC m ) T (GC m G T + C d ) -1 (GC m ) where μ m|d is the posterior mean solution constrained by seismic data, C m|d is the posterior covariance constrained by seismic data; C m is the prior covariance matrix; the superscript "-1" represents the inverse operation symbol of the matrix; d is the input seismic data; G is the linear forward operator; the superscript "T" represents the transpose operation symbol of the matrix; μ m represents the expectation that the parameter m to be inverted satisfies; C d represents the covariance matrix of the error term e′ conforming to the Gaussian distribution.

Images