Shale reservoir geological engineering double-dessert earthquake identification method

Abstract

The invention discloses a shale reservoir geological engineering double-dessert earthquake identification method. The method comprises the steps of obtaining a new longitudinal wave reflection coefficient equation of a VTTI medium containing fluid bulk modulus, brittleness index, pressure correlation parameters, inclined fracture density and horizontal fracture density parameters; the inclined fracture density eTTI, the fluid volume modulus Kf, the brittleness index BI, the pressure correlation parameter sigma E and the horizontal fracture density eVTI are obtained through an inversion method. The direct earthquake prediction method based on the horizontal fracture and inclined fracture media (VTTI media) realizes a plurality of dessert parameters covering shale gas fluid parameters, crustal stress parameters, brittleness indexes and horizontal fracture and inclined fracture densities.

Description

Technical Field

[0001] The invention belongs to the technical field of oil and gas field development engineering, and in particular relates to a double sweet spot seismic identification method for 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 method mainly calculates the sweet spot parameters such as brittleness index and fracture density indirectly after inverting the elastic parameters from seismic data. Its limitations are mainly reflected in: (1) It is necessary to use the inversion results of elastic parameters to indirectly calculate the shale gas sweet spot parameters, such as the brittleness index. There are accumulated errors in this calculation process, resulting in inaccurate prediction of sweet spot parameters; (2) It is impossible to simultaneously realize the direct seismic identification of multiple sweet spot parameters, and it is urgent to improve the simultaneous prediction capability of multiple complex sweet spot parameters; (3) When predicting the fracture density sweet spot parameters, only simple vertical fracture media or orthogonal fracture (horizontal fractures and vertical fractures) media are considered, without considering the situation that horizontal fractures and inclined fractures may develop simultaneously in actual fractured shale reservoirs. Summary of the invention

[0003] The purpose of the present invention is to provide a shale reservoir geological engineering double sweet spot seismic identification method 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 by providing a shale reservoir geological engineering double sweet spot seismic identification method, including: obtaining a new equation of the longitudinal wave reflection coefficient of the VTTI medium including fluid bulk modulus, brittleness index, pressure correlation parameter, inclined fracture density and horizontal fracture density parameters, the new equation is:

[0005]

[0006] in,

[0007] Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock;

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

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

[0010]

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

[0012]

[0013]

[0014]

[0015]

[0016]

[0017]

[0018]

[0019]

[0020]

[0021]

[0022]

[0023]

[0024]

[0025] In the formula, is the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, The stiffness coefficient matrix of the dry rock skeleton is represented by Stiffness coefficient at different positions in the VTI is the horizontal crack density; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; e TTI is the inclined crack density; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; ζ is the inclination 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 The normal crack is weak, δ T is the tangential crack weakness.

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

[0027]

[0028] 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 crack 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.

[0029] In some embodiments, the elastic stiffness element of any anisotropic fluid-saturated rock is calculated as follows:

[0030]

[0031] in,

[0032]

[0033]

[0034]

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

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

[0037]

[0038] In the formula, 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;

[0039] Will Substitution get:

[0040]

[0041] In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of rock; K m is the bulk modulus of the matrix, 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;

[0042] Will Substitution get:

[0043]

[0044]

[0045]

[0046]

[0047] β4=β6=0

[0048] 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.

[0049] In some embodiments, based on the weak anisotropic medium theory and the porous 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:

[0050]

[0051] 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 cracks, is the normal weakness of the inclined fracture; φ is the total porosity of the rock; K f is the bulk modulus of the fluid;

[0052] based on Approximately equal to:

[0053]

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

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

[0056]

[0057] Where φ0 is the initial empirical porosity, σ V is the vertical effective stress, is the pressure-related parameter, denoted by σ E β l is the effective stress coefficient; e is the natural constant;

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

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

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

[0061] The formula

[0062]

[0063]

[0064]

[0065] β4=β6=0 and the formula Substitute into the formula In the intermediate process of this substitution simplification, ignore all and The relevant terms give an approximate formula for the stiffness coefficient of saturated rock:

[0066]

[0067]

[0068]

[0069]

[0070]

[0071]

[0072]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078]

[0079] In the formula, The stiffness coefficient matrix of the saturated VTTI rock is represented by Stiffness coefficient at different positions in the f is the bulk modulus of the fluid; δ ij and δ kl is the Kronecker symbol; C dry is the stiffness coefficient matrix of dry VTTI medium; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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 related parameter.

[0080] In some embodiments, based on the derivation of the approximate formula of the stiffness coefficient of saturated rock, according to the principle of weak anisotropy and weak contrast or small disturbance of the background elastic modulus at the interface, the formulas of the disturbance amount corresponding to the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock are obtained respectively:

[0081]

[0082]

[0083]

[0084]

[0085]

[0086]

[0087]

[0088]

[0089]

[0090]

[0091]

[0092]

[0093]

[0094] In the formula, and is the disturbance of the saturated rock stiffness coefficient between the upper and lower strata at different positions in the stiffness coefficient matrix; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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-related parameter;

[0095] Based on the Born steady phase method, the linearization function is:

[0096]

[0097] In the formula, R PP is the P-wave reflection coefficient; θ is the P-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 strata; represents the known coefficients related to the angle of incidence and azimuth; is the disturbance of the stiffness coefficient at the Ith row and Jth column in the saturated rock stiffness coefficient matrix, which specifically includes and

[0098] Substitute the formula for the disturbance of the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock into the formula The longitudinal wave reflection coefficient of VTTI medium is obtained as:

[0099] in,

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106] In the formula, is the longitudinal wave reflection coefficient of VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal crack reflection coefficient, is the reflection coefficient of the anisotropic TTI inclined crack; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined crack; Φ is the azimuth; K f is the bulk modulus of the fluid, μ b is the shear modulus of dry rock skeleton, σ E is the pressure correlation parameter, ρ is the rock density; is the bulk modulus of the fluid, is the shear modulus, is the pressure correlation parameter, k ρ is the term coefficient of rock density; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δμ b Indicates the change in shear modulus between the upper and lower strata; Δσ E It represents the change of the pressure correlation parameter between the upper and lower strata; Δρ represents the change of the rock density between the upper and lower strata; It indicates the change in the normal weakness of horizontal fractures between the upper and lower strata; Indicates the change in tangential weakness of horizontal fractures between upper and lower strata; It indicates the change in the normal weakness of the inclined fracture between the upper and lower strata; It indicates the change in the tangential weakness of the inclined fracture between the upper and lower strata; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of shear modulus between upper and lower strata; It represents the average value of the pressure correlation parameter between the upper and lower formations; It represents the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the term coefficient of the normal weakness of the inclined crack, 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.

[0107] In some embodiments, based on Formula G n (φ p )=(1-K dry / K m ) 2 / φp =φ p / φ c 2 =σ E φ0 / φ c 2 ,set up get:

[0108]

[0109] In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the bulk modulus of the fluid; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus 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-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock;

[0110] Formula Find μ b The total differential of

[0111]

[0112] In the formula, μ b is the shear modulus of saturated rock; Δμ b Indicates the change in shear modulus between the upper and lower strata; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δσ E Indicates the change in the pressure correlation parameters between the upper and lower formations; Indicates the space between upper and lower strata The amount of change in the parameter; represents the average value of shear modulus; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of the pressure correlation parameter between the upper and lower formations; Indicates the space between upper and lower strata 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 saturated rock; K f is the bulk modulus of the fluid, σ Eis the pressure-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock;

[0113] Introducing the brittleness index BI is the brittleness index, and are Young's modulus and first Lame coefficient of saturated rock respectively; because the brittleness index can be obtained from It is expressed as: Therefore, we can obtain the BI The total differential of , we get:

[0114]

[0115] Where 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 saturated rock; is the ratio of the P-wave modulus to the shear modulus in saturated rock; ΔBI represents the change in the brittleness index between the upper and lower strata; Indicates the space between upper and lower strata The amount of change in the parameter; It represents the average value of brittleness index between upper and lower strata; Indicates the space between upper and lower strata The mean value of the parameter.

[0116] The formula and formula Substitute into the formula The formula Substitute into the formula and formula Finally, a new equation for the longitudinal 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:

[0117]

[0118] in,

[0119] Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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, It is the square of the ratio of the P-wave velocity to the S-wave velocity in saturated rock.

[0120] In some embodiments, the inverted inclined fracture density e TTI , including: using two vertical orthogonal seismic data to make a difference, obtaining the seismic amplitude difference Δd as the seismic data d input in the formula d=Gm+e′, and then substituting Φ1=0° and Φ1=90° into the formula to obtain the coefficient matrix A in the formula G=WAD Subtract the result from , and we get:

[0121]

[0122] In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle;

[0123] The formula Simplifying, we get:

[0124]

[0125] In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the azimuthal perturbation coefficient of the inclined fracture density; e TTI is the inclined crack density; Δe TTI It represents the disturbance of the density of inclined fractures between the upper and lower strata; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; Φ is the azimuth; The term coefficient representing the density of inclined fractures;

[0126] Inversion of inclined fracture density e TTI The coefficient matrix used is expressed as:

[0127]

[0128] in,

[0129] Where, M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points of the inclined crack density; is the azimuthal perturbation coefficient of the inclined fracture density; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; K is the coefficient matrix composed of the azimuthal perturbation coefficient of the inclined fracture density under a single incident angle; the superscript "T" is the transposition operator of the matrix; A1 is the coefficient matrix composed of the azimuthal perturbation coefficient for the inverted inclined fracture density; the symbol "diag" is a diagonal matrix in which the elements are arranged diagonally;

[0130] The formula A1 in the formula G=WAD is used as the coefficient matrix A of the parameters to be inverted, and then the parameters to be inverted are expressed as:

[0131]

[0132] Where m1 is the parameter matrix to be inverted consisting of the inclined fracture density; e TTI is the inclined fracture density; subscripts 1, 2, 3, ..., N represent the sampling point numbers; N is the number of sampling points of the inclined fracture density; the superscript "T" represents the transposition operator of the matrix;

[0133] Finally, the formula m1 in the formula d = Gm + e' is used as the parameter m to be inverted, and the inclined fracture density e can be inverted. TTI ; and / or

[0134] Invert fluid bulk modulus, brittleness index, pressure correlation parameters, and horizontal fracture density parameters:

[0135] Using the processed residual seismic data d re As the formula d = Gm + e'

[0136] The seismic data d in is input, and then the coefficient matrix of the parameters to be inverted is expressed as:

[0137]

[0138] in,

[0139]

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

[0141]

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

[0143]

[0144] Where A2 is the coefficient matrix composed of the coefficients of the inverted fluid bulk modulus, brittleness index, pressure correlation parameter and horizontal fracture density; M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points for each parameter; superscript "T" is the transpose operator of the matrix; θ is the longitudinal wave incident angle; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the term coefficient of rock density, is the term coefficient of horizontal crack density, is the term coefficient of inclined fracture density; the symbol "diag" is a diagonal matrix in which the elements are arranged diagonally; F is the coefficient matrix composed of the term coefficients of the fluid bulk modulus at a single incident angle; X is the coefficient matrix composed of the term coefficients of the brittleness index at a single incident angle; V is the coefficient matrix composed of the term coefficients of the pressure correlation parameter at a single incident angle; U is the coefficient matrix composed of the term coefficients of the rock density at a single incident angle; Y is the coefficient matrix composed of the term coefficients of the inclined fracture density at a single incident angle;

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

[0146]

[0147] in,

[0148]

[0149]

[0150]

[0151]

[0152]

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

[0154] K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined fracture density; 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 transposition operator;

[0155] Finally, m2 is used as the parameter m to be inverted in the formula d = Gm + e' to invert the fluid bulk modulus K. f , brittleness index BI, pressure correlation parameter σ E and horizontal crack density e VTI .

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

[0157]

[0158] Where p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, which indicates the degree of fit 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;

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

[0160] 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;

[0161] The forward model of seismic data is given as: d = Gm + e'

[0162] 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 is expected to be 0, and the covariance matrix is ​​C d Gaussian distribution of

[0163] Since the forward operator G is linear, the likelihood function p(d|m) obeys a Gaussian distribution. Assuming that the prior distribution of the parameter m obeys a Gaussian distribution, it satisfies the expectation μ m and the prior covariance matrix C m , satisfying m~N(μ m ,C m ),

[0164]

[0165] In the formula, μ m is the expectation satisfied by the inverted parameter m; 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 of the time variation is t l , the function that changes with time is v(τ), then C t It is expressed as:

[0166]

[0167] In the formula, C t is the time-domain spatial correlation matrix; t lis the length that changes with time; v is the function that changes with time;

[0168] The prior probability distribution function and likelihood function are expressed as follows:

[0169]

[0170]

[0171] 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 inversion parameter m; the superscript “-1” is the inversion operator of the matrix; the superscript T represents the transposition operator of the matrix; d is the input seismic data; C m represents the prior covariance matrix satisfied by the inversion parameter m; C d Indicates that the error term e′ conforms to the covariance matrix of Gaussian distribution; e is a natural constant; G is a linear forward operator; the symbol || indicates the absolute value;

[0172] when 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 as:

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

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

[0175] In the formula, μ 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” indicates the matrix inversion operator; d is the input seismic data; G is the linear forward operator; the superscript “T” indicates the matrix transposition operator; μm represents the expectation that the inversion parameter m satisfies; C d The covariance matrix representing the error term e′ conforms to the Gaussian distribution.

[0176] In this application:

[0177] OA: orthotropic anisotropy; VTI: vertical transverse isotropy; TTI: inclined transverse isotropy; VTTI medium: monoclinic fracture medium taking into account the VTI background, that is, inclined transverse isotropic (TTI) fracture formation medium embedded in vertical transverse isotropy (VTI). Further, it can be understood that: in order to distinguish the monoclinic medium from the ordinary TTI monoclinic medium, the present invention refers to the monoclinic fracture medium taking into account the VTI background as VTTI medium.

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

[0179] Based on Gassmann anisotropic fluid substitution equation, this application determines the stiffness coefficient matrix of saturated VTTI medium with high precision. In this process, the fluid bulk modulus, brittleness index, vertical effective stress related parameters, and horizontal and inclined fracture density are combined to derive a new linear PP wave reflection coefficient equation, which provides an equation basis for directly and simultaneously predicting multiple sweet spot parameters of geological engineering integration; 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. Based on horizontal fracture and inclined fracture media (VTTI media), this application realizes a direct seismic prediction method covering multiple sweet spot parameters of shale gas fluid parameters, ground stress parameters, brittleness index, horizontal fracture and inclined fracture density. BRIEF DESCRIPTION OF THE DRAWINGS

[0180] Figure 1 A comparison diagram of the approximate value (red dashed line) and the accurate value (blue solid line) under different stiffness coefficients provided in an embodiment of the present invention;

[0181] Figure 2 A comparison diagram of the approximate value (red dashed line) and the accurate value (blue solid line) under different stiffness coefficients provided in another embodiment of the present invention;

[0182] Figure 3 A well logging interpretation data diagram of a vertical well A in a shale reservoir provided by an embodiment of the present invention;

[0183] Figure 4 A graph showing elastic parameters and predicted fracturing parameters of a vertical well A provided in an embodiment of the present invention;

[0184] Figure 5 A seismic stacking data diagram of any survey line passing through Well A provided in an embodiment of the present invention;

[0185] Figure 6 A diagram of azimuth amplitude difference data at different central incident angles provided by an embodiment of the present invention;

[0186] Figure 7 This is a diagram of the inversion results of elasticity and crack parameters of any measuring line provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0187] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with 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.

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

[0189] The present application embodiment provides a shale reservoir geological engineering double sweet spot seismic identification method, including:

[0190] S01. Obtaining a new equation for the longitudinal 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, wherein the new equation is:

[0191]

[0192] in,

[0193] Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK fis the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock;

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

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

[0196] Based on Gassmann anisotropic fluid substitution equation, this application determines the stiffness coefficient matrix of saturated VTTI medium with high precision. In this process, the fluid bulk modulus, brittleness index, vertical effective stress related parameters, and horizontal and inclined fracture density are combined to derive a new linear PP wave reflection coefficient equation, which provides an equation basis for directly and simultaneously predicting multiple sweet spot parameters of geological engineering integration; 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. Based on horizontal fracture and inclined fracture media (VTTI media), this application realizes a direct seismic prediction method covering multiple sweet spot parameters of shale gas fluid parameters, ground stress parameters, brittleness index, horizontal fracture and inclined fracture density.

[0197] The method of the present application was applied to the example of shale fracture reservoirs in the Sichuan Basin of China, showing reliable results.

[0198] In the S01:

[0199] In some embodiments, for a VTTI medium containing inclined cracks and horizontal cracks, the elastic stiffness coefficient matrix is ​​obtained as follows:

[0200]

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

[0202]

[0203]

[0204]

[0205]

[0206]

[0207]

[0208]

[0209]

[0210]

[0211]

[0212]

[0213]

[0214]

[0215] In the formula, is the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, The stiffness coefficient matrix of the dry rock skeleton is represented by Stiffness coefficient at different positions in the VTI is the horizontal crack density; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; e TTI is the inclined crack density; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; ζ is the inclination 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 The normal crack is weak, δ T is the tangential crack weakness.

[0216] Understandably, is the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, The stiffness coefficient matrix of the dry rock skeleton is represented by The stiffness coefficients at different positions in the matrix are explained in detail. Subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, and 66 all represent the row and column positions of the elements of different stiffness coefficients in the stiffness coefficient matrix. Taking subscript 11 as an example, Represents the stiffness coefficient of the first row and first column in the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium.

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

[0218]

[0219] 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 crack 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.

[0220] Understandably, is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton.

[0221] It can be understood that for anisotropic fractured rocks, based on the low-frequency Gassmann poroelasticity theory, the elastic stiffness element of any anisotropic fluid-saturated rock is obtained. The calculation formula of the elastic stiffness element of any anisotropic fluid-saturated rock is:

[0222]

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

[0224] In some embodiments, the elastic stiffness element of any anisotropic fluid-saturated rock is calculated as follows:

[0225]

[0226] in,

[0227]

[0228]

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

[0230] 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, J = 1, It represents the stiffness coefficient of the first row and first column in the stiffness coefficient matrix of the dry rock skeleton of the VTTI medium, and so on; ∑ is the summation symbol.

[0231] Furthermore, the formulas of the elastic stiffness coefficient matrix and the matrix elements in the elastic stiffness coefficient matrix are substituted into get:

[0232]

[0233] In the formula, 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;

[0234] Will Substitution get:

[0235]

[0236] In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of rock; K m is the bulk modulus of the matrix, 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;

[0237] Will Substitution get:

[0238]

[0239]

[0240]

[0241]

[0242] β4=β6=0

[0243] 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.

[0244] Furthermore, based on the weak anisotropic medium theory and porous rock physics theory, it is believed that the fracture weakness parameter is close to 0 and the fluid bulk modulus is much smaller than the matrix bulk modulus, and we can obtain:

[0245]

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

[0247] based on Approximately equal to:

[0248]

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

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

[0251]

[0252] Where φ0 is the initial empirical porosity, σ V is the vertical effective stress, is the pressure-related parameter, denoted by σ E β l is the effective stress coefficient; e is the natural constant;

[0253] It can be understood that based on the classic coin-shaped fracture model theory, when the fracture aspect ratio is close to 0, the fracture porosity φ f Relative to the background porosity φ p It is almost negligible, so the total porosity is approximately equal to the background porosity, that is, φ≈φ p .

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

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

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

[0257] The formula

[0258]

[0259]

[0260]

[0261] β4=β6=0 and the formula Substitute into the formula In the intermediate process of this substitution simplification, ignore all and The relevant terms give an approximate formula for the stiffness coefficient of saturated rock:

[0262]

[0263]

[0264]

[0265]

[0266]

[0267]

[0268]

[0269]

[0270]

[0271]

[0272]

[0273]

[0274]

[0275] In the formula, The stiffness coefficient matrix of the saturated VTTI rock is represented by Stiffness coefficient at different positions in the f is the bulk modulus of the fluid; δ ij and δ kl is the Kronecker symbol; C dry is the stiffness coefficient matrix of dry VTTI medium; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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 related parameter.

[0276] Understandably, The stiffness coefficient matrix of the saturated VTTI rock is represented by The stiffness coefficients at different positions in the matrix are explained in detail. Subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, and 66 all represent the row and column positions of the elements of different stiffness coefficients in the stiffness coefficient matrix. Taking subscript 11 as an example, The stiffness coefficient matrix representing the saturated VTTI rock is The stiffness coefficient of the first row and first column in .

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

[0278]

[0279]

[0280]

[0281]

[0282]

[0283] The approximate stiffness coefficient calculated is the same as that given by the formula The matrix minerals used in the rock physics model experiment include calcite (whose bulk modulus is 76.8 GPa, shear modulus is 32 GPa, and density is 2.71 g / cm 3 ), quartz (whose bulk modulus is 37 GPa, shear modulus is 44 GPa, and density is 2.65 g / cm 3 ) and clay (bulk modulus 21 GPa, shear modulus 7 GPa, density 2.6 cm 3 The background pores and crack pores are filled with different proportions of water (whose bulk modulus is 2.2 GPa, shear modulus is 0, and density is 1.02 g / cm 3 ) and gas (gas bulk modulus is 0.002GPa, shear modulus is 0, density is 0.00065g / cm 3 ) is a mixed fluid. When the gas saturation is 80%, the density of horizontal and inclined fractures is 0.25φ, and the inclination angles of the inclined fractures are 60° (see Figure 1 ) and 90° (see Figure 2 ) Figure 1 and Figure 2 Known, by the formula:

[0284]

[0285]

[0286]

[0287]

[0288]

[0289] The approximate stiffness coefficients calculated are shown as red dashed lines, which are consistent with the equation The calculated accurate stiffness coefficient is shown as the blue solid line, under different rock total porosity φ, Figure 1 Different stiffness coefficients C shown in (a)-(f) 11 , C 12 , C 13 , C 22 , C 23 , C 33 The approximate and exact values ​​of are in good agreement, proving that in gas-bearing shale reservoirs with high gas content, low porosity, and high quartz content, the formula

[0290]

[0291]

[0292]

[0293]

[0294] The calculated approximate stiffness coefficient is highly accurate and meets the requirements;

[0295] Furthermore, on the basis of the derived approximate formula of the stiffness coefficient of saturated rock, according to the principle of weak anisotropy and weak contrast or small disturbance of the background elastic modulus at the interface, the formulas of the disturbance amount corresponding to the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock are obtained:

[0296]

[0297]

[0298]

[0299]

[0300]

[0301]

[0302]

[0303]

[0304]

[0305]

[0306]

[0307]

[0308]

[0309] In the formula, and is the disturbance of the saturated rock stiffness coefficient between the upper and lower strata at different positions in the stiffness coefficient matrix; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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 related parameter.

[0310] Understandably, and is the disturbance of the saturated rock stiffness coefficient at different positions in the stiffness coefficient matrix between the upper and lower strata. The subscripts 11, 12, 13, 15, 22, 23, 25, 33, 35, 44, 46, 55, and 66 represent the row and column positions of the components of different stiffness coefficients in the stiffness coefficient matrix. Taking subscript 11 as an example, The stiffness coefficient matrix representing saturated rock is The disturbance amount of the stiffness coefficient between the upper and lower strata at the 1st row and 1st column in .

[0311] Furthermore, the linearization function is based on the Born phase stability method, namely:

[0312]

[0313] In the formula, R PP is the P-wave reflection coefficient; θ is the P-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 strata; represents the known coefficients related to the angle of incidence and azimuth; is the disturbance of the stiffness coefficient at the Ith row and Jth column in the saturated rock stiffness coefficient matrix, which specifically includes and

[0314] Substitute the formula for the disturbance of the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock into the formula The longitudinal wave reflection coefficient of the VTTI medium of the present application is obtained as follows:

[0315] in,

[0316]

[0317]

[0318]

[0319]

[0320]

[0321]

[0322] In the formula, is the longitudinal wave reflection coefficient of VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal crack reflection coefficient, is the reflection coefficient of the anisotropic TTI inclined crack; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined crack; Φ is the azimuth; K f is the bulk modulus of the fluid, μ b is the shear modulus of dry rock skeleton, σ E is the pressure correlation parameter, ρ is the rock density; is the bulk modulus of the fluid, is the shear modulus, is the pressure correlation parameter, k ρ is the term coefficient of rock density; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δμ b Indicates the change in shear modulus between the upper and lower strata; Δσ E It represents the change of the pressure correlation parameter between the upper and lower strata; Δρ represents the change of the rock density between the upper and lower strata; It indicates the change in the normal weakness of horizontal fractures between the upper and lower strata; Indicates the change in tangential weakness of horizontal fractures between upper and lower strata; It indicates the change in normal weakness of inclined fractures between upper and lower strata; It indicates the change in the tangential weakness of the inclined fracture between the upper and lower strata; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of shear modulus between upper and lower strata; It represents the average value of the pressure correlation parameter between the upper and lower formations; It represents the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the term coefficient of the normal weakness of the inclined crack, 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.

[0323] Furthermore, based on formula G n (φ p )=(1-K dry / K m ) 2 / φ p =φ p / φ c 2 =σ E φ0 / φ c 2 ,set up get:

[0324]

[0325] In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the bulk modulus of the fluid; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus 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-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock.

[0326] Furthermore, the formula Find μ b The total differential of

[0327]

[0328] In the formula, μ b is the shear modulus of saturated rock; Δμ b Indicates the change in shear modulus between the upper and lower strata; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δσ E Indicates the change in the pressure correlation parameters between the upper and lower formations; Indicates the space between upper and lower strata The amount of change in the parameter; represents the average value of shear modulus; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of the pressure correlation parameter between the upper and lower formations; Indicates the space between upper and lower strata 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 saturated rock; K f is the bulk modulus of the fluid, σ E is the pressure-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock;

[0329] Introducing the brittleness index BI is the brittleness index, and are Young's modulus and first Lame coefficient of saturated rock respectively; because the brittleness index can be obtained from It is expressed as: Therefore, we can obtain the BI The total differential of , we get:

[0330]

[0331] Where 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 saturated rock; is the ratio of the P-wave modulus to the shear modulus in saturated rock; ΔBI represents the change in the brittleness index between the upper and lower strata; Indicates the space between upper and lower strata The amount of change in the parameter; It represents the average value of brittleness index between upper and lower strata; Indicates the space between upper and lower strata The mean value of the parameter.

[0332] The formula and formula Substitute into the formula

[0333] , the formula Substitute into the formula and formula Finally, a new equation for the longitudinal 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:

[0334]

[0335] in,

[0336] Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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, It is the square of the ratio of the P-wave velocity to the S-wave velocity in saturated rock.

[0337] In the S02:

[0338] In some embodiments, the inverted inclined fracture density e TTI, including: using two vertical orthogonal seismic data to make a difference, obtaining the seismic amplitude difference Δd as the seismic data d input in the formula d=Gm+e′, and then substituting Φ1=0° and Φ1=90° into the formula to obtain the coefficient matrix A in the formula G=WAD Subtract the result from the above to get:

[0339]

[0340] In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle;

[0341] The formula Simplifying, we get:

[0342]

[0343]

[0344] In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the azimuthal perturbation coefficient of the inclined fracture density; e TTI is the inclined crack density; Δe TTI represents the disturbance of the density of inclined fractures between the upper and lower strata; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; Φ is the azimuth; k eTTI The term coefficient representing the density of inclined fractures;

[0345] Inversion of inclined fracture density e TTI The coefficient matrix used is expressed as:

[0346]

[0347] in,

[0348] Where, M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points of the inclined crack density; is the azimuthal perturbation coefficient of the inclined fracture density; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; K is the coefficient matrix composed of the azimuthal perturbation coefficient of the inclined fracture density under a single incident angle; the superscript "T" is the transposition operator of the matrix; A1 is the coefficient matrix composed of the azimuthal perturbation coefficient for the inverted inclined fracture density; the symbol "diag" is a diagonal matrix in which the elements are arranged diagonally;

[0349] The formula A1 in the formula G=WAD is used as the coefficient matrix A of the parameters to be inverted, and then the parameters to be inverted are expressed as:

[0350]

[0351] Where m1 is the parameter matrix to be inverted consisting of the inclined fracture density; e TTI is the inclined fracture density; subscripts 1, 2, 3, ..., N represent the sampling point numbers; N is the number of sampling points of the inclined fracture density; the superscript "T" represents the transposition operator of the matrix;

[0352] Finally, the formula m1 in the formula d = Gm + e' is used as the parameter m to be inverted, and the inclined fracture density e can be inverted. TTI .

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

[0354]

[0355] Where p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, which indicates the degree of fit 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;

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

[0357] 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;

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

[0359] 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 is expected to be 0, and the covariance matrix is ​​C d Gaussian distribution.

[0360] It is understandable that:

[0361] Since the forward operator G is linear, the likelihood function p(d|m) obeys a Gaussian distribution. Assuming that the prior distribution of the parameter m obeys a Gaussian distribution, it satisfies the expectation μ m and the prior covariance matrix C m , satisfying m~N(μ m ,C m ),

[0362]

[0363] In the formula, μ m is the expectation satisfied by the inverted parameter m; 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 of the time variation is t l , the function that changes with time is v(τ), then C t It is expressed as:

[0364]

[0365] In the formula, C t is the time-domain spatial correlation matrix; t l is the length that changes with time; v is the function that changes with time;

[0366] The prior probability distribution function and likelihood function are expressed as follows:

[0367]

[0368]

[0369] 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 inversion parameter m; the superscript “-1” is the inversion operator of the matrix; the superscript T represents the transposition operator of the matrix; d is the input seismic data; C m represents the prior covariance matrix satisfied by the inversion parameter m; C d Indicates that the error term e′ conforms to the covariance matrix of Gaussian distribution; e is a natural constant; G is a linear forward operator; the symbol || indicates the absolute value;

[0370] when 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 as:

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

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

[0373] In the formula, μ 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” indicates the matrix inversion operator; d is the input seismic data; G is the linear forward operator; the superscript “T” indicates the matrix transposition operator; μ m represents the expectation that the inversion parameter m satisfies; C d The covariance matrix representing the error term e′ conforms to the Gaussian distribution.

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

[0375] Using the processed residual seismic data d re As the formula d = Gm + e'

[0376] The seismic data d in is input, and then the coefficient matrix of the parameters to be inverted is expressed as:

[0377]

[0378] in,

[0379]

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

[0381]

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

[0383]

[0384] Where A2 is the coefficient matrix composed of the coefficients of the inverted fluid bulk modulus, brittleness index, pressure correlation parameter and horizontal fracture density; M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points for each parameter; superscript "T" is the transpose operator of the matrix; θ is the longitudinal wave incident angle; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the term coefficient of rock density, is the term coefficient of horizontal crack density, is the term coefficient of inclined fracture density; the symbol "diag" is a diagonal matrix in which the elements are arranged diagonally; F is the coefficient matrix composed of the term coefficients of the fluid bulk modulus at a single incident angle; X is the coefficient matrix composed of the term coefficients of the brittleness index at a single incident angle; V is the coefficient matrix composed of the term coefficients of the pressure correlation parameter at a single incident angle; U is the coefficient matrix composed of the term coefficients of the rock density at a single incident angle; Y is the coefficient matrix composed of the term coefficients of the inclined fracture density at a single incident angle;

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

[0386]

[0387] in,

[0388]

[0389]

[0390]

[0391]

[0392]

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

[0394] K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined fracture density; 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 transposition operator;

[0395] Finally, m2 is used as the parameter m to be inverted in the formula d = Gm + e' to invert the fluid bulk modulus K. f , brittleness index BI, pressure correlation parameter σ E and horizontal crack density e VTI .

[0396] Application Examples

[0397] The actual work area is located in the Sichuan Basin of China, which is an unconventional fractured shale gas reservoir in the Longmaxi-Wufeng formation. A vertical well A has been drilled on the 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 the figure, (a) FMI imaging, (b) core photo, (c) inclined fracture orientation rose diagram and (d) fracture dip statistical histogram are shown in Figure 2. Figure 3 It can be seen that the logging imaging data of the vertical well A shows that natural fractures are developed in the reservoir, including horizontal fractures in the bedding and high-angle inclined fractures. The fractures are mainly oriented from southwest to northeast, and the inclination angles of the inclined fractures are mainly concentrated at 70°. Therefore, it is appropriate to regard the gas-bearing shale reservoir in this work area as a VTTI medium to verify the effectiveness of the proposed method. Figure 4 In the figure, (a) mineral component content; (b) P wave velocity; (c) S wave velocity; (d) rock density; (e) brittleness index; (f) horizontal fracture density; (g) inclined fracture density; (h) TOC content; (i) fluid bulk modulus; (j) pressure correlation parameter; Figure 4 It can be seen that the mineral components in the reservoir mainly include calcite, quartz, dolomite, clay, pyrite and kerogen, etc. The elastic parameters and fracture parameters can describe the reservoir position (red curve part) as low P-wave velocity, low S-wave velocity, low density, high brittleness index, high horizontal fracture density, low fluid bulk modulus, high pressure correlation parameters and high inclined fracture density. There are two groups of horizontal and high-angle fractures in the shale reservoir. The natural fractures are mainly SW-NE, and the inclination of the inclined fractures is mainly concentrated at 70°. Therefore, it is appropriate to regard the gas-bearing shale reservoir in this work area as a VTTI medium to verify the effectiveness of the proposed method.

[0398] First, azimuth seismic data are extracted from the OVT domain, along the natural fracture strike and its vertical direction, i.e., azimuths of 45° and 135°. Then, the azimuth seismic data are stacked at 3°, 9°, 15°, 21°, and 27° as the central incident angles to improve the signal-to-noise ratio of the seismic data. Finally, partial central incident angle stacked seismic data under these two azimuths are obtained, see Figure 5 , Figure 5 contains different azimuths and central incident angles, where (a)-(e) observation azimuth Φ abs1 is 45°, (f)-(j) observation direction Φ abs2 The central incident angle is 135°, which 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°, Figure 5 It can be seen that the azimuth seismic data used in this example application is divided into two azimuths: 45 degrees and 135 degrees, and the two azimuths are superimposed at the central incidence angles of 3°, 9°, 15°, 21°, and 27°, respectively, to obtain seismic superposition data of different azimuths and different central incidence angles, which serve as the seismic data basis for subsequent inversion. Among the five central incidence angles, we give priority to three relatively large central incidence angles, namely 15°, 21°, and 27°, and use the azimuth amplitude difference seismic data generated between the above two azimuths, see Figure 6 , Figure 6 In (a), the central incident angle θ3 is 15°, (b) the central incident angle θ3 is 21°, and (c) the central incident angle θ3 is 27°. Figure 6 It can be seen that the two Figure 5 At two different orientations, the central incident angles are 15° and 21° respectively.

[0399] By subtracting the seismic profile data at different azimuths of 27°, we can obtain the amplitude difference seismic data of azimuths at different central incident angles. It can be seen at the stratified position of the T3 well that the amplitude difference seismic data shows obvious anomalies, which is closely related to the development of natural fractures in the reservoir. In order to invert the density of anisotropic inclined fractures. Figure 6 In the reservoir position of well layer T2-T3, it can be observed that there is an obvious amplitude difference between the two azimuths, 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 ).

[0400] from Figure 7 It can be seen that Figure 7 The fluid bulk modulus in (a) shows an abnormally low value, indicating a high shale gas content; Figure 7 The brittleness index in (b) shows an abnormally high value, indicating that the shale rock has good brittleness; Figure 7 The stress correlation parameters in (c) show abnormally high values, which must indicate that the vertical effective pressure is small due to the high shale gas content and high pore pressure; Figure 7 The shale gas density in (d) is abnormally high. Figure 7 The horizontal crack density and Figure 7 The high value abnormality of the inclined fracture density in (f) indicates the development of horizontal and inclined natural fractures in the reservoir. Therefore, the reservoir has the characteristics of low fluid bulk modulus, high brittleness index, high effective stress-related parameters, and high horizontal and inclined fracture density, and the prediction is consistent with the actual geological background and objective understanding.

[0401] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A double sweet spot seismic identification method for shale reservoir geological engineering, characterized in that: include: A new equation for the longitudinal 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, and the new equation is: in, Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock; Obtain the inclined fracture density e by the inversion method TTI , fluid bulk modulus K f , brittleness index BI, pressure correlation parameter σ E and horizontal crack density e VTI .

2. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 1 is characterized in that: For the VTTI medium containing inclined cracks and horizontal cracks, the elastic stiffness coefficient matrix is ​​obtained as follows: Among them, 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, The stiffness coefficient matrix of the dry rock skeleton is represented by Stiffness coefficient at different positions in the VTI is the horizontal crack density; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; e TTI is the inclined crack density; is the normal weakness of the inclined crack, is the tangential weakness of the inclined crack; is the P-wave modulus in the dry rock skeleton, is the shear modulus in the dry rock skeleton, is the first Lame coefficient in the dry rock skeleton; ζ is the inclination 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 The normal crack is weak, δ T is the tangential crack weakness.

3. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 2 is characterized in that: For dry or gas-saturated cracks, the crack weakness is: 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 crack 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.

4. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 2 is characterized in that: The calculation formula of the elastic stiffness element of arbitrary anisotropic fluid-saturated rock is expressed as follows: in, In the formula, C sat is the stiffness coefficient of saturated rock, C dry is the stiffness coefficient of the dry rock skeleton; φ is the total porosity of the rock; K m is the bulk modulus of the matrix, K f is the bulk modulus of the fluid; is the anisotropic Gassmann pore modulus; β and All are intermediate variable parameters without special meanings; subscripts I, J are the row and column positions of the stiffness coefficient elements in the stiffness coefficient matrix, and Σ is the summation symbol.

5. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 4 is characterized in that: Substitute the formula of the elastic stiffness coefficient matrix and the matrix elements in the elastic stiffness coefficient matrix into get: In the formula, 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; Will Substitution get: In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of rock; K m is the bulk modulus of the matrix, 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; Will Substitution get: β4＝β6＝0 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.

6. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 5, characterized in that: Based on the weak anisotropic medium theory and porous rock physics theory, it is believed 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: 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 cracks, is the normal weakness of the inclined fracture; φ is the total porosity of the rock; K f is the bulk modulus of the fluid; based on Approximately equal to: In the formula, is the anisotropic Gassmann pore modulus; φ is the total porosity of rock; K f is the bulk modulus of the fluid; Since φ≈φ p , so according to the relationship between porosity and vertical effective stress, the quantitative relationship between porosity and vertical effective stress is obtained: Where φ0 is the initial empirical porosity, σ V is the vertical effective stress, is the pressure-related parameter, denoted by σ E β l is the effective stress coefficient; e is the natural constant; Based on the gain function related to the bulk modulus of the dry rock skeleton, the matrix bulk modulus and the porosity, we obtain: G n (f p )=(1-K dry / K m ) 2 / f p =φ p / f c 2 =s E φ0 / φ c 2 In the formula, G n is the gain function; φ p is the background porosity; K dry is the bulk modulus of dry rock skeleton, K m is the bulk modulus of the matrix; φ c is the critical porosity; φ0 is the initial empirical porosity; σ E is the pressure-related parameter; The formula β4=β6=0 and the formula Substitute into the formula I, J = 1, 2, ..., 6. In the intermediate process of this substitution simplification, ignore all and The relevant terms give an approximate formula for the stiffness coefficient of saturated rock: In the formula, The stiffness coefficient matrix of the saturated VTTI rock is represented by Stiffness coefficient at different positions in the f is the bulk modulus of the fluid; δ ij and δ kl is the Kronecker symbol; C dry is the stiffness coefficient matrix of dry VTTI medium; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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 related parameter.

7. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 6, characterized in that: On the basis of the derived approximate formula of the stiffness coefficient of saturated rock, according to the principle of weak anisotropy and weak contrast or small disturbance of the background elastic modulus at the interface, the formulas of the disturbance amount corresponding to the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock are obtained respectively: In the formula, and is the disturbance of the saturated rock stiffness coefficient between the upper and lower strata at different positions in the stiffness coefficient matrix; is the P-wave modulus in the dry rock skeleton, is the first Latin coefficient in the dry rock skeleton; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus ζ is the inclination 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 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-related parameter; Based on the Born steady phase method, the linearization function is: In the formula, R PP is the P-wave reflection coefficient; θ is the P-wave incident angle; ρ is the rock density; S(r0) is the position vector on the interface between two weakly isotropic media; Δρ Indicates the change in rock density between upper and lower strata; represents the known coefficients related to the angle of incidence and azimuth; is the disturbance of the stiffness coefficient at the Ith row and Jth column in the saturated rock stiffness coefficient matrix, which specifically includes and Substitute the formula for the disturbance of the stiffness coefficient in the approximate formula of the stiffness coefficient of saturated rock into the formula The longitudinal wave reflection coefficient of VTTI medium is obtained as: in, In the formula, is the longitudinal wave reflection coefficient of VTTI medium; R iso is the isotropic reflection coefficient, is the anisotropic VTI horizontal crack reflection coefficient, is the reflection coefficient of the anisotropic TTI inclined crack; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined crack; Φ is the azimuth; K f is the bulk modulus of the fluid, μ b is the shear modulus of dry rock skeleton, σ E is the pressure correlation parameter, ρ is the rock density; is the bulk modulus of the fluid, is the shear modulus, is the pressure correlation parameter, k ρ is the term coefficient of rock density; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δμ b Indicates the change in shear modulus between the upper and lower strata; Δσ E It represents the change of the pressure correlation parameter between the upper and lower strata; Δρ represents the change of the rock density between the upper and lower strata; It indicates the change in the normal weakness of horizontal fractures between the upper and lower strata; Indicates the change in tangential weakness of horizontal fractures between upper and lower strata; It indicates the change in normal weakness of inclined fractures between upper and lower strata; It indicates the change in the tangential weakness of the inclined fracture between the upper and lower strata; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of shear modulus between upper and lower strata; It represents the average value of the pressure correlation parameter between the upper and lower formations; It represents the average value of rock density between the upper and lower strata; ζ is the inclination 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 saturated rock; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the normal weakness of the horizontal crack, is the tangential weakness of the horizontal crack; is the term coefficient of the normal weakness of the inclined crack, 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.

8. The shale reservoir geological engineering double sweet spot seismic identification method 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 ,set up get: In the formula, is the P-wave modulus in the saturated rock skeleton; K f is the bulk modulus of the fluid; μ b is the shear modulus of saturated rock, μ b Equal to the dry rock skeleton shear modulus 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-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock; Formula Find μ b The total differential of In the formula, μ b is the shear modulus of saturated rock; Δμ b Indicates the change in shear modulus between the upper and lower strata; ΔK f Indicates the change in bulk modulus of the fluid between the upper and lower strata; Δσ E Indicates the change in the pressure correlation parameters between the upper and lower formations; Indicates the space between upper and lower strata The amount of change in the parameter; represents the average value of shear modulus; It represents the average value of the bulk modulus of the fluid between the upper and lower formations; It represents the average value of the pressure correlation parameter between the upper and lower formations; Indicates the F between the upper and lower strata ω 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 saturated rock; K f is the bulk modulus of the fluid, σ E is the pressure-related parameter; is the ratio of the P-wave modulus to the shear modulus in saturated rock; Introducing the brittleness index BI is the brittleness index, and are Young's modulus and first Lame coefficient of saturated rock respectively; because the brittleness index can be obtained from It is expressed as: Therefore, we can obtain the BI The total differential of , we get: Where 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 saturated rock; is the ratio of the P-wave modulus to the shear modulus in saturated rock; ΔBI It indicates the variation of brittleness index between upper and lower strata; Indicates the space between upper and lower strata The amount of change in the parameter; It represents the average value of brittleness index between upper and lower strata; Indicates the space between upper and lower strata The mean value of the parameter. The formula and formula Substitute into the formula The formula Substitute into the formula and formula Finally, a new equation for the longitudinal 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: in, Where: is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined crack density; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the coefficient of rock density, is the term coefficient of horizontal crack density and is the term coefficient of inclined crack density; is the term coefficient of the normal weakness of horizontal cracks, is the term coefficient of tangential weakness of horizontal cracks; is the term coefficient of the normal weakness of the inclined crack, is the term coefficient of the tangential weakness of the inclined fracture; ΔK f is the disturbance of the bulk modulus of the fluid between the upper and lower strata, ΔBI is the disturbance of the brittleness index between the upper and lower strata, Δσ E is the disturbance of the pressure correlation parameter between the upper and lower strata, Δρ is the disturbance of the rock density between the upper and lower strata, Δe VTI is the disturbance of horizontal fracture density between the upper and lower strata, Δe TTI is the disturbance of the inclined fracture density between the upper and lower strata; is the average value of the bulk modulus of the fluid between the upper and lower strata, is the average value of the brittleness index between the upper and lower strata, is the average value of the pressure correlation parameter between the upper and lower formations, is the average value of rock density between the upper and lower strata; ζ is the inclination 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, It is the square of the ratio of the P-wave velocity to the S-wave velocity in saturated rock.

9. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 1, characterized in that: Inversion of inclined fracture density e TTI , including: using two vertical orthogonal seismic data to make a difference, obtaining the seismic amplitude difference Δd as the seismic data d input in the formula d=Gm+e′, and then substituting Φ1=0° and Φ1=90° into the formula to obtain the coefficient matrix A in the formula G=WAD Subtract the result from , and we get: In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the longitudinal wave reflection coefficient of VTTI medium; θ is the longitudinal wave incident angle; ζ is the inclination angle of the inclined crack; Φ is the azimuth angle; The formula Simplifying, we get: In the formula, is the azimuthally disturbed reflection coefficient caused by TTI inclined cracks; is the azimuthal perturbation coefficient of the inclined fracture density; e TTI is the inclined crack density; Δe TTI It represents the disturbance of the density of inclined fractures between the upper and lower strata; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; Φ is the azimuth; The term coefficient representing the density of inclined fractures; Inversion of inclined fracture density e TTI The coefficient matrix used is expressed as: in, Where, M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points of the inclined crack density; is the azimuthal perturbation coefficient of the inclined fracture density; θ is the incident angle of the longitudinal wave; ζ is the inclination angle of the inclined fracture; K is the coefficient matrix composed of the azimuthal perturbation coefficient of the inclined fracture density under a single incident angle; the superscript "T" is the transposition operator of the matrix; A1 is the coefficient matrix composed of the azimuthal perturbation coefficient for the inverted inclined fracture density; the symbol "diag" is a diagonal matrix in which the elements are arranged diagonally; The formula A1 in the formula G=WAD is used as the coefficient matrix A of the parameters to be inverted, and then the parameters to be inverted are expressed as: Where m1 is the parameter matrix to be inverted consisting of the inclined fracture density; e TTI is the inclined fracture density; subscripts 1, 2, 3, ..., N represent the sampling point numbers; N is the number of sampling points of the inclined fracture density; superscript "T" represents the transposition operator of the matrix; Finally, the formula The m1 in the formula d = Gm + e' is used as the parameter m to be inverted, and the inclined fracture density e can be inverted. TTI ; and / or Invert fluid bulk modulus, brittleness index, pressure correlation parameters, and horizontal fracture density parameters: Using the processed residual seismic data d re As the formula d = Gm + e' The seismic data d in is input, and then the coefficient matrix of the parameters to be inverted is expressed as: in, X(θ i )=diag[k BI (i i )k BI (i i )…k BI (i i )] (N-1)×(N-1) U(θ i )=diag[k ρ (i i )k ρ (i i )…k ρ (i i )] (N-1)×(N-1) Where A2 is the coefficient matrix composed of the coefficients of the inverted fluid bulk modulus, brittleness index, pressure correlation parameter and horizontal fracture density; M is the number of incident angles; subscript i is the serial number of the incident angle; N is the number of sampling points for each parameter; superscript "T" is the transposition operator of the matrix; θ is the longitudinal wave incident angle; is the coefficient of the bulk modulus of the fluid, k BI is the term coefficient of the brittleness index, is the term coefficient of the pressure-related parameter, k ρ is the term coefficient of rock density, is the term coefficient of horizontal crack density, is the term coefficient of inclined fracture density; the symbol "diag" is a diagonal matrix with elements arranged diagonally; F is the coefficient matrix composed of the term coefficients of the fluid bulk modulus at a single incident angle; X is the coefficient matrix composed of the term coefficients of the brittleness index at a single incident angle; V is the coefficient matrix composed of the term coefficients of the pressure correlation parameter at a single incident angle; U is the coefficient matrix composed of the term coefficients of the rock density at a single incident angle; Y is the coefficient matrix composed of the term coefficients of the inclined fracture density at a single incident angle; The formula A2 in the formula G = the coefficient matrix A input in WAD, and the parameters to be inverted in this step are expressed as: in, In the formula, is the parameter matrix to be inverted of the bulk modulus of the fluid, m BI is the parameter matrix to be inverted of the brittleness index, is the parameter matrix to be inverted of the pressure-related parameters, mρ is the parameter matrix to be inverted of rock density, The parameter matrix to be inverted for horizontal fracture density; K f is the bulk modulus of the fluid, BI is the brittleness index, σ E is the pressure correlation parameter, ρ is the rock density, e VTI is the horizontal crack density, e TTI is the inclined fracture density; 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 transposition operator; Finally, m2 is used as the parameter m to be inverted in the formula d = Gm + e' to invert the fluid bulk modulus K. f , brittleness index BI, pressure correlation parameter σ E and horizontal crack density e VTI .

10. The shale reservoir geological engineering double sweet spot seismic identification method according to claim 1, characterized in that: According to Bayesian theory, the posterior probability distribution of the parameters to be inverted is expressed as follows: Where p(m) is the prior Gaussian distribution function of the parameter to be inverted; p(d|m) is the likelihood function, which indicates the degree of fit 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 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; The forward model of 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 is expected to be 0, and the covariance matrix is ​​C d Gaussian distribution of Since the forward operator G is linear, the likelihood function p(d|m) obeys a Gaussian distribution. Assuming that the prior distribution of the parameter m obeys a Gaussian distribution, it satisfies the expectation μ m and the prior covariance matrix C m , satisfying m~N(μ m ,C m ) In the formula, μ m is the expectation satisfied by the inverted parameter m; 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 of the time variation is t l , the function that changes with time is v(τ), then C t It is expressed as: In the formula, C t is the time-domain spatial correlation matrix; t l is the length that changes with time; v is the function that changes with time; The prior probability distribution function and likelihood function are 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 satisfied by the inversion parameter m; the superscript "-1" is the symbol of the matrix inversion operation; the superscript T represents the symbol of the matrix transposition operation; d is the input seismic data; C m represents the prior covariance matrix satisfied by the inversion parameter m; C d The covariance matrix of the error term e′ conforms to the Gaussian distribution; e is a natural constant; G is a linear forward operator; the symbol || represents the absolute value; when 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 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 ) In the formula, μ 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" indicates the matrix inversion operator; d is the input seismic data; G is the linear forward operator; the superscript "T" indicates the matrix transposition operator; μ m represents the expectation satisfied by the inversion parameter m; C d The covariance matrix representing the error term e′ conforms to the Gaussian distribution.

Images