Shale solid skeleton and organic matter elastic parameter synchronous prediction method

Abstract

The present application is suitable for oil and gas geophysical exploration field, and provides a shale solid skeleton and organic matter elastic parameter synchronous prediction method, the method comprises the following steps: constructing a petrophysical model representing the elastic modulus of shale reservoir solid skeleton and the elastic modulus of organic matter; developing an inversion method in the well based on the petrophysical model; combining the petrophysical model and the scattering function theory, establishing a reflection coefficient formula linearly characterized by the elastic modulus of shale solid skeleton, the elastic modulus of organic matter and the density parameter; combining the reflection coefficient formula and the Bayesian linear inversion theory, developing a seismic petrophysical inversion method; applying synthetic seismic data, testing and verifying the effectiveness of the seismic petrophysical inversion method constructed in step E; then applied to actual seismic data, predicting the elastic modulus of shale gas reservoir solid skeleton and organic matter in the study area, providing an effective method for the brittleness evaluation of shale gas reservoir engineering sweet spot and the organic matter enrichment prediction of geological sweet spot.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil and gas geophysical exploration, and particularly relates to a shale solid skeleton and organic matter elastic parameter synchronous prediction method. BACKGROUND

[0002] As an important unconventional oil and gas resource, shale gas has been widely concerned by the industry and academia. In order to effectively evaluate the shale gas reservoir, it is necessary to carry out shale gas reservoir engineering sweet spot and geological sweet spot prediction. The main content of engineering sweet spot prediction is the evaluation of the mechanical properties of shale solid skeleton, and the main content of geological sweet spot prediction is to identify shale gas reservoirs with high total organic carbon content (TOC). At present, in the prediction of shale gas reservoir engineering sweet spot, the conventional method mainly calculates the mechanical parameters on the basis of seismic elastic parameter inversion, and the calculation results reflect the mechanical properties of the whole shale. However, due to the influence of organic matter and fluid in shale and other factors, this method cannot completely reflect the mechanical properties of the solid skeleton itself. In the prediction of geological sweet spot, the current method mainly applies the empirical relationship between TOC and conventional seismic inversion elastic parameters (such as density parameters) to predict the TOC content of shale. The main problem of this method is that the universality of the constructed empirical relationship is limited in different regions, and at the same time, the seismic inversion method of density parameters has inherent problems such as poor stability and low calculation accuracy. Therefore, on the basis of rock physics research of shale gas reservoir, it is necessary to develop a seismic inversion method directly for the elastic parameters of shale solid skeleton and organic matter, and to provide an effective method for the evaluation of brittleness of shale gas reservoir engineering sweet spot and the prediction of organic matter enrichment of geological sweet spot. SUMMARY

[0003] In view of the deficiencies of the prior art, the purpose of the present application is to provide a shale solid skeleton and organic matter elastic parameter synchronous prediction method to solve the technical problems mentioned in the background art.

[0004] A shale solid skeleton and organic matter elastic parameter synchronous prediction method, the method comprises:

[0005] Step A, on the basis of microstructure analysis of shale gas reservoir, a rock physics model for characterizing the elastic modulus of shale reservoir solid skeleton ( M frame , μ frame ) and the elastic modulus of organic matter ( M om , μ om ) is constructed;

[0006] Step B: Based on the rock physics model constructed in Step A, develop a well inversion method to predict the elastic modulus of the solid skeleton from conventional logging data (mineral composition, porosity, fluid saturation, P-wave and S-wave velocities). M frame , μ frame ) and the elastic modulus of organic matter ( M om , μ om This provides a physically meaningful initial model for the inversion of actual pre-stack seismic data;

[0007] Step C: Combining the rock physics model and scattering function theory constructed in Step A, establish the elastic modulus of the shale solid skeleton ( M frame , μ frame ), organic matter elastic modulus ( M om , μ om ) and density parameters ( ρ sat The formula for the linearized reflection coefficient;

[0008] Step D: Conduct an accuracy evaluation of the reflection coefficient formula and a sensitivity analysis of the parameters in the reflection coefficient formula (reflectivity of the elastic modulus of the solid skeleton, reflectivity of the elastic modulus of organic matter, and density reflectivity of organic-rich shale) to verify the applicability of the reflection coefficient formula;

[0009] Step E: Combining the reflection coefficient formula established in Step C with Bayesian linear inversion theory, a seismic rock physics inversion method is developed to achieve direct synchronous prediction of the elastic modulus of the shale solid skeleton from pre-stack seismic data. M frame , μ frame ), organic matter elastic modulus ( M om , μ om ) and other parameters;

[0010] Step F: Apply the predicted elastic modulus of the solid skeleton and the elastic modulus of organic matter from Step B to the actual seismic wavelet to calculate and synthesize seismic data, and test and verify the effectiveness of the seismic rock physics inversion method constructed in Step E.

[0011] Step G involves applying the seismic rock physics inversion method developed in Step E and the initial model constructed in Step B to actual pre-stack seismic data to calculate the elastic modulus of the solid skeleton and organic matter of shale gas reservoirs in the study area. This provides a new and effective method for evaluating the brittleness of engineering sweet spots in shale gas reservoirs and predicting the enrichment of organic matter in geological sweet spots.

[0012] Preferably, in step A, the step of constructing a petrophysical model characterizing the horizontal fracture density of the shale gas reservoir comprises:

[0013] Applying Hashin-Shtrikman boundary (HSB) theory to calculate the bulk modulus and shear modulus of the solid matrix composed of complex minerals K 0, μ 0);

[0014] Using self-consistent approximation (SCA) theory to calculate the bulk modulus and shear modulus of the porous matrix composed of the solid matrix and water-saturated inorganic pores K matrix , μ matrix );

[0015] Using Krief formula to calculate the bulk modulus and shear modulus of the solid skeleton K frame , μ frame );

[0016] Using Kuster-Toksöz (KT) formula to calculate the bulk modulus and shear modulus of the organic matter composed of kerogen and gas-saturated organic pores K mix , μ mix );

[0017] The overall bulk modulus and shear modulus of the organic-rich shale composed of the solid skeleton and the organic matter K sat , μ sat ) can be obtained by using Ciz-Shapiro solid replacement theory.

[0018] Preferably, in step B, the well inversion method comprises:

[0019] Two key variables are applied in the petrophysical inversion process: the volume fraction of organic matter f om and the aspect ratio of inorganic pores α ; by fitting the measured and petrophysical model calculated compressional and shear wave velocities, the optimal volume fraction of organic matter f om and the aspect ratio of inorganic pores α are obtained by minimizing the objective function, and the objective function F obj is as follows:

[0020] Equation 18: ;

[0021] where, V P and V S represent the P- and S-wave velocities, respectively; the subscripts "modeled" and "measured" refer to the velocities calculated from the petrophysical modeling and measured, respectively, and the total porosity is set to φ the inorganic and organic pore porosities, respectively, i.e. φ inorganic and φ organic = φ φ organic + φ inorganic ;

[0022] The predicted f om and α values are used as input to the petrophysical model to calculate the elastic modulus of the solid matrix in the well ( M frame and μ frame ) and of the organic matter ( M om and μ om ), and the petrophysical inversion results provide an initial model with physical meaning for the subsequent seismic petrophysical inversion.

[0023] Preferably, in step C, the content of the reflectivity formulation linearized in terms of the elastic modulus of the shale solid matrix ( M frame , μ frame ), the elastic modulus of the organic matter ( M om , μ om ) and the density parameter ( ρ sat ) includes:

[0024] Shaw and Sen (2004, 2006) derived the expression for the PP-wave reflection coefficient at an interface between any weakly anisotropic medium R PP and the scattering function S by combining the perturbation of the stiffness matrix of the saturated rock and the Born approximation:

[0025] Equation 19: ;

[0026] where θ and​ψ represent the incidence and azimuth angles, respectively, ρ sat represent the density of the saturated rock, and r0is the point on the horizontal interface that satisfies Snell's law; the scattering function S is related to the perturbation of the elastic parameters, the horizontal slowness, and the polarization direction, and its expression is:

[0027] Equation 20: ;

[0028] where S (r0) is the scattering function, and ρ sat is the perturbation of the density of the saturated rock, and C ijkl is the perturbation of the stiffness tensor of the saturated rock, and the variable ξ is related to the polarization vector of the PP wave, and the variable η ijkl is related to the polarization vector and the slowness vector of the PP wave;

[0029] According to Shaw and Sen (2006), the dummy index "i" in Equation 20 is simplified to the dummy index "j" according to the following formula: ijkl IJ

[0030] Equation 21: ;

[0031] Equation 22: ;

[0032] where ijkl represents the subscript value of the fourth-order stiffness tensor, i , j , k , l take values of 1 to 3, respectively; IJ represents the simplified subscript value of the fourth-order stiffness tensor, I , J take values of 1 to 6, respectively.

[0033] where δ ij and δ kl are the Kronecker symbols, so Equation 20 can be simplified as:

[0034] Equation 23: ;

[0035] where S (r0) is the scattering function, and ρ sat is the perturbation of the density of the saturated rock, and C ​​ijkl is the perturbation of the stiffness tensor of the saturated rock, and the variable ξ is a physical quantity related to the polarization vector of the PP wave, and the variable η ijkl is a physical quantity related to the polarization vector of the PP wave and the slowness vector

[0036] is the polarization vector of the PP wave t and t may be expressed by the incident and reflected waves, respectively, as

[0037] Equation 24: ;

[0038] Equation 25: ;

[0039] where θ and ψ represent the incidence angle and the azimuth angle, respectively.

[0040] its corresponding slowness vector p and p may be expressed by the polarization vector as

[0041] Equation 26: ;

[0042] Equation 27: ;

[0043] where α sat is the P-wave velocity of the saturated rock;

[0044] When the PP wave is incident and reflected, the expression of the variable ξ is written as

[0045] Equation 28: ;

[0046] t i and represent the first t and t components of the polarization vector i and

[0047] and the variable η IJ is written as

[0048] Equation 29: ;

[0049] where represents the first t component of the polarization vector i , represents the slowness vectorp the jth component of the polarization vector t, tkrepresents the kth component of the polarization vector t, p l represents the lth component of the slowness vector p; η IJ The expression for each component is:

[0050] Equation 30: ;

[0051] η 11 , η 12 , η 13 , η 22 , η 23 , η 33 , η 44 , η 55 , η 66 , η 21 , η 31 , η 32 represent the variables η IJ in different directions.

[0052] For a weakly isotropic medium, the perturbation Δ C IJ of the stiffness matrix in Equation 23 can be simplified to the form:

[0053] Equation 31: ;

[0054] where Δ M sat , Δ λ sat and Δ μ sat are the perturbations of the saturated rock's longitudinal wave modulus M sat , the first Lamé constant λ sat and the second Lamé constant μ sat , respectively;

[0055] Substituting Equations 28, 30 and 31 into Equation 23, the scattering function S for a weakly isotropic medium is obtained, which is expressed as:

[0056] Equation 32: ;

[0057] where Δ M sat , Δ λ sat and Δ μ sat are the perturbation of the longitudinal wave modulus M sat , the first Lamé constant λ sat and the second Lamé constant μ sat of the saturated rock, respectively;

[0058] Finally, substituting equation 32 into equation 19, we obtain the PP wave reflection coefficient approximation formula for any weakly isotropic medium:

[0059] For the petrophysical model, if the media on both sides of the interface satisfy the weak reflection approximation condition, then the perturbation of the longitudinal wave modulus and the shear modulus of the organic-rich shale are:

[0060] Equation 33: ;

[0061] Equation 34: ;

[0062] where Δ M frame and Δ μ frame represent the perturbation of the longitudinal wave modulus and the shear modulus of the solid skeleton, respectively, and Δ M om and Δ μ om represent the perturbation of the corresponding elastic parameters of the organic matter, respectively. The perturbations satisfy the following conditions: |Δ M frame | M frame , |Δ μ frame | μ frame , |Δ M om | M om , |Δ μ om | μ om ;

[0063] Therefore, substituting equations 32, 33 and 34 into equation 19, we obtain the following PP wave linearized reflection coefficient approximation formula:

[0064] Equation 35: ;

[0065] where R PP is the reflection coefficient, θ is the incident angle, Δ ρ sat is the perturbation of the saturated rock density.

[0066] Substituting equations 13 and 17 into equation 35, we have

[0067] Equation 36: ;

[0068] Then, define γ 2 sat = M sat / μ sat and γ 2 frame = M frame / μ frame , which represent the square of the ratio of the P-wave and S-wave velocities of the organic-rich shale and the solid matrix, respectively. Therefore, the above equation can be simplified as

[0069] Equation 37: ;

[0070] Finally, define the parameter p = μ frame / μ sat as the ratio of the shear modulus of the solid matrix and the organic-rich shale. Therefore, equation 37 can be simplified, and the final linearized expression of the reflection coefficient is

[0071] Equation 38: ;

[0072] For simplicity of expression and convenience of subsequent inversion studies, equation 38 can be simplified as

[0073] Equation 39: ;

[0074] The coefficient term in the equation is defined as

[0075] Equation 40: ;

[0076] The above equation represents the case of the organic matter with non-zero shear modulus as a solid filler in the organic-rich shale. For the case of fluid filler, i.e. Figure 2when the remaining pores are filled with fluid, it satisfies the condition p = 1 and the Gassmann equation; under this condition, the shear modulus of the solid matrix is equal to the shear modulus of the whole rock, i.e. μ frame = μ sat Therefore, p = 1 and p = 1 are brought into formula 38 as follows: μ frame = μ sat Formula 41:

[0077] Formula 41: ;

[0078] According to the poroelastic medium theory proposed by Russell et al. (2011), the Gassmann equation describes the relationship between the longitudinal wave modulus of the solid matrix M frame and the shear modulus of the saturated rock μ sat , i.e.

[0079] Formula 42: ;

[0080] Then, the reflectivity of the variable Mframe is calculated, and there is:

[0081] Formula 43: ;

[0082] Therefore, formula 41 is simplified as:

[0083] Formula 44: ;

[0084] According to the reflectivity formula derived by Russell et al.:

[0085] Formula 45: ;

[0086] Formula 44 and formula 45 have the same form, where Δ M om / M om is equivalent to Δ f / f in formula 45.

[0087] Preferably, in step E, the content of the seismic rock physics inversion method includes:

[0088] According to Buland and Omre (2003) and Hampson and Russell (2005), the solid matrix elastic modulus reflectivity Δ Mframe / M frame and Δ μ frame / μ frame , reflectivity of organic matter elastic modulus Δ M om / M om and Δ μ om / μ om , reflectivity of organic matter elastic modulus Δ ρ sat / ρ sat The corresponding natural logarithm value can be used to approximate, i.e.:

[0089] Equation 46: ;

[0090] In the formula, the symbol "ln(·)" represents the natural logarithm of the elastic parameter, and the symbol "Δ" represents the difference between the values on both sides of the reflection interface;

[0091] In addition, the elastic properties of the solid matrix, the density of the organic matter-rich shale, and the elastic properties related to the organic matter are extended to continuous time variables, as follows:

[0092] Equation 47: ;

[0093] Therefore, the reflection coefficient of a single interface represented by Equation 37 is extended to a continuous-time reflection coefficient equation by Equations 45 and 46:

[0094] Equation 48: ;

[0095] In the formula, the symbol "d" represents the difference between the elastic parameters at time t and time t+1 in continuous time; for example, dln[ M frame ( t )] = ln[ M frame ( t +1)] - ln[ M frame ( t )];

[0096] Subsequently, the length of the seismic data in time is set to n+1, and the number of incidence angles is h, so there are n reflection interfaces in the time direction of the seismic data, and the angle-dependent seismic wavelet matrix W is represented as:

[0097] Equation 49: ;

[0098] Block matrix W in matrix W i Including the angle of incidence θ i The relevant seismic wavelet is expressed as follows:

[0099] Formula 50: ;

[0100] in w i Indicates the angle of incidence as θ i Seismic wavelet at time, i =1, 2, ..., h. s 1, s 2, ..., s l The time sampling points represent the seismic wavelet;

[0101] The seismic wavelet matrix W is convolved with the reflection coefficient sequence to obtain the seismic forward modeling equation with data error terms, as shown below:

[0102] Formula 51: ;

[0103] In the formula, S represents the seismic response, W is the seismic wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, and e represents the data residual term, which satisfies the condition that the mean is 0 and the variance is Σ. d The Gaussian distribution is given by m, which is the model parameter vector, and its expression is:

[0104] Formula 52: ;

[0105] Where the superscript T denotes the transpose of the vector, and the length of vector m is... q = 5( n +1), this is because each of the n+1 strata contains five elastic parameters, which ultimately form n reflection points;

[0106] The coefficient matrix Λ has a dimension of (h×n)×5n, and its expression is:

[0107] Formula 53: ;

[0108] in A ( θ i ), B ( θ i ), C ( θ i ), D ( θi ) and E ( θ i ) is a diagonal matrix with the expression:

[0109] Equation 54: ;

[0110] where a ( t j , θ i ), b ( t j , θ i ), c ( t j , θ i ), d ( t j , θ i ) and e ( t j , θ i ) are coefficients varying with angle and time, where j = 1, 2,..., n and i = 1, 2,..., h ;

[0111] The difference matrix D has a dimension of 5n x 5(n+1) and the expression is:

[0112] Equation 55: ;

[0113] where the block matrix Ddhas a dimension of n x (n+1) and the expression is:

[0114] Equation 56: ;

[0115] According to the Bayesian theory, the posterior probability density function P(m|d) is obtained from the prior probability density function P(d|m) and the likelihood function P(m):

[0116] Equation 57: ;

[0117] In the present application, the prior model parameter vector m is set to satisfy a multi-dimensional Gaussian distribution with a mean of μ m and a covariance of Σ m , and the data residual term satisfies a Gaussian distribution with a mean of 0 and a covariance of Σ dGaussian distribution, therefore, the posterior probability density function P(m|d) satisfies Gaussian distribution, the mean μ m|d and covariance Σ m|d The expression is:

[0118] Equation 58: ;

[0119] Equation 59: ;

[0120] Where the matrix G = WΛD is the forward operator;

[0121] Therefore, the maximum likelihood solution m opt The expression is:

[0122] Equation 60: ;

[0123] Where μ m|d And Σ m|d The mean and covariance of the posterior probability density function P(m|d).

[0124] The beneficial effects realized by the present application are:

[0125] The present application realizes direct synchronous prediction of shale solid skeleton elastic modulus, organic matter elastic modulus and other parameters from prestack seismic data. Synthetic seismic data is applied to test and verify the effectiveness of the developed seismic rock physics inversion method. Finally, the developed seismic rock physics inversion method and the constructed initial model are applied to actual prestack seismic data to calculate the elastic modulus of the solid skeleton and organic matter of the shale gas reservoir in the study area, which provides a new effective method for the brittleness evaluation of shale gas reservoir engineering sweet spots and the organic matter enrichment prediction of geological sweet spots. BRIEF DESCRIPTION OF DRAWINGS

[0126] Figure 1 It is a technical flow chart of a shale solid skeleton and organic matter elastic parameter synchronous prediction method;

[0127] Figure 2 It is a rock physics modeling flow chart constructed according to geological and core data;

[0128] Figure 3 The change of the longitudinal wave modulus and shear modulus of the organic matter-rich shale with the organic matter content f om And the inorganic pore aspect ratio α The longitudinal wave modulus and shear modulus of the organic matter-rich shale M sat And μ sat With the organic matter content f oma and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir α M frame a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir μ frame a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir f om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir α

[0129] Figure 4 a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir f om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir α M om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir μ om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir f om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir α f om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir M om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir μ om a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir

[0130] Figure 5

[0131] Figure 6

[0132] Figure 7 a and b are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir

[0133] Figure 8

[0134] Figure 9 ​​​​​​​is the actual well log of Well A in Example 2, a is the natural gamma ray curve, b is the P-wave velocity curve, c is the S-wave velocity curve, d is the density curve, e is the total porosity curve, and f is the mineral volume fraction curve;

[0135] Figure 10 is the petrophysical inversion result of Well A in Example 2, a is the P-wave velocity curve, b is the S-wave velocity curve, where the red curve and the black curve represent the petrophysical predicted and measured velocities, respectively, c is the volume fraction of organic matter, i.e., the organic matter content f om , d is the inorganic pore aspect ratio, and e is the porosity curve, where the red curve and the black curve represent the organic pore porosity and the total porosity, respectively;

[0136] Figure 11 is the target curve of petrophysical modeling based on the petrophysical inversion result of Well A in Example 2, where a is the P-wave modulus curve, b is the shear modulus curve, c is the density curve, the red curve and the black curve represent the properties of shale solid skeleton and organic-rich shale, respectively, d and e represent the organic matter-related P-wave modulus and shear modulus curves, respectively;

[0137] Figure 12 is the synthetic seismic data constructed based on Figure 10 in Example 3, a is the solid skeleton P-wave modulus curve, b is the solid skeleton shear modulus curve, c is the density curve of organic-rich shale, d is the organic matter-related P-wave modulus curve, e is the organic matter-related S-wave modulus curve, and f is the synthetic seismic data;

[0138] Figure 13 is the inversion result of the synthetic seismic data in Example 3, a is the solid skeleton P-wave modulus curve, b is the solid skeleton shear modulus curve, c is the density curve of organic-rich shale, d is the organic matter-related P-wave modulus curve, e is the organic matter-related S-wave modulus curve, where the black curve, the red curve, and the blue curve represent the actual model, the inversion result, and the initial model, respectively, and f is the synthetic seismic data, where the red curve represents the seismic data simulated by forward modeling using the inversion result, and the black curve represents the original synthetic seismic data;

[0139] Figure 14 is the post-stack seismic profile of Well A in Example 4;

[0140] Figure 15 is the initial model of Well A in Example 4, a is the solid skeleton P-wave modulus, b is the solid skeleton shear modulus, c is the density of organic-rich shale, d is the organic matter-related P-wave modulus, and e is the organic matter-related S-wave modulus;

[0141] Figure 16These are the inversion results from well A in Experiment Example 4. a is the P-wave modulus of the solid skeleton, b is the shear modulus of the solid skeleton, c is the density of organic-rich shale, d is the P-wave modulus related to organic matter, and e is the S-wave modulus related to organic matter. Detailed Implementation

[0142] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0143] Please see Figure 1 This invention provides a method for simultaneously predicting the elastic parameters of shale solid skeleton and organic matter, the method comprising:

[0144] Step A: Based on the microstructure analysis of shale gas reservoirs, construct a characterizing elastic modulus of the solid skeleton of shale reservoirs. M frame , μ frame ) and the elastic modulus of organic matter ( M om , μ om Rock physics model;

[0145] Step B: Based on the rock physics model constructed in Step A, develop a well inversion method to predict the elastic modulus of the solid skeleton from conventional logging data (mineral composition, porosity, fluid saturation, P-wave and S-wave velocities). M frame , μ frame ) and the elastic modulus of organic matter ( M om , μ om This provides a physically meaningful initial model for the inversion of actual pre-stack seismic data;

[0146] Step C: Combining the rock physics model and scattering function theory constructed in Step A, establish the elastic modulus of the shale solid skeleton ( M frame , μ frame ), organic matter elastic modulus ( M om , μ om ) and density parameters ( ρ sat The formula for the linearized reflection coefficient;

[0147] Step D, carry out precision evaluation and sensitivity analysis of the reflection coefficient formula, verify the applicability of the reflection coefficient formula;

[0148] Step E, further combine the reflection coefficient formula established in step C and the Bayesian linear inversion theory to develop a seismic rock physics inversion method to realize direct and synchronous prediction of shale solid skeleton elastic modulus (G0) M frame , μ frame ), organic matter elastic modulus (Gorg) M om , μ om , etc. from prestack seismic data;

[0149] Step F, apply the predicted solid skeleton elastic modulus and organic matter elastic modulus in step B to calculate the synthetic seismic data combined with the actual seismic wavelet, test and verify the effectiveness of the seismic rock physics inversion method constructed in step E;

[0150] Step G, apply the seismic rock physics inversion method developed in step E and the initial model constructed in step B to the actual prestack seismic data to calculate the elastic modulus of the solid skeleton and organic matter of the shale gas reservoir in the study area, providing a new effective method for the brittleness evaluation of shale gas reservoir engineering sweet spots and the organic matter enrichment prediction of geological sweet spots.

[0151] Please refer to Figure 2 In step A of the embodiment, the step of constructing a rock physics model representing the horizontal fracture density of the shale gas reservoir includes:

[0152] Applying Hashin-Shtrikman boundary (HSB) theory to calculate the bulk modulus and shear modulus (K0 K 0, μ 0) of the solid matrix composed of complex minerals;

[0153] Using self-consistent approximation (SCA) theory to calculate the bulk modulus and shear modulus (Kp K matrix , μ matrix ) of the pore matrix composed of solid matrix and water-saturated inorganic pores;

[0154] Using Krief formula to calculate the bulk modulus and shear modulus (K0 K frame , μ frame ) of the solid skeleton;

[0155] The Kuster-Toksöz (KT) formula was used to calculate the bulk modulus and shear modulus of organic matter and gas-saturated organic pores from kerogen. K mix , μ mix );

[0156] The overall bulk modulus and shear modulus of organic-rich shale composed of a solid framework and organic matter ( K sat , μ sat It can be obtained using the Ciz-Shapiro solid substitution theory.

[0157] In step A of this embodiment, the Hashin-Shtrikman boundary (HSB) theory includes the following:

[0158] For a medium composed of two different components, its equivalent bulk modulus and shear modulus can be calculated using the HSB theory:

[0159] Formula 1: ;

[0160] in K HS± and μ HS± These are represented as the upper boundary + and lower boundary – of the equivalent medium bulk modulus, respectively. K i and μ i Representing the first i The elastic modulus of the component; f i Representing the i The volume fraction of each component; Makvo et al. (2009) pointed out that for isotropic mixtures containing two-phase media, the theoretical boundary calculated by HSB theory is the narrowest.

[0161] For an equivalent medium composed of multiple components, its equivalent elastic modulus can be calculated using the following formula:

[0162] Formula 2: ;

[0163] in K HS+ and K HS- These represent the upper and lower boundaries of the equivalent medium's bulk modulus, respectively. μ HS+ and μ HS- These represent the upper and lower boundaries of the equivalent medium shear modulus, respectively. K maxand μ max These represent the maximum values ​​of the bulk modulus and shear modulus of each component constituting the equivalent medium, respectively. K min and μ min These represent the minimum values ​​of the bulk modulus and shear modulus of each component that make up the equivalent medium, respectively.

[0164] Formula 3: ;

[0165] in, z The independent variable in the function, when calculating the upper limit, z for K max or μ max When calculating the lower limit, z is K min or μ min Λ( z ), Γ( z )and ζ ( K , μ () represents a function in the calculation process. K Represents bulk modulus and shear modulus; r Representing the first in the rock r The composition has several values, which are related to the rock physics modeling process. The symbol <·> represents calculating a volume-weighted average for each composition.

[0166] Bulk modulus of equivalent medium K 0 and shear modulus μ 0 can be calculated from the average of the theoretical lower and lower limits of HSB, that is:

[0167] Formula 4: .

[0168] In step A of this embodiment, the Kuster-Toksöz (KT) formula includes the following:

[0169] For an equivalent medium containing inclusions of various material properties or different shapes within a background medium, Kuster and Toksöz (1974) provided expressions for its bulk modulus and shear modulus. In this invention, kerogen is assumed to be filled with gas-saturated organic pores to form an organic mixture; therefore, its elastic modulus is expressed as:

[0170] Formula 5: ;

[0171] in K * mix and μ *mix These represent the bulk modulus and shear modulus of the organic mixture, respectively. K k and μ k These represent the bulk modulus and shear modulus of kerogen, respectively. K g and μ g The bulk modulus and shear modulus of gases; the shape of organic pores is determined by the structure factor. P kg and G kg For its specific expression, please refer to Mavko et al. (2009). x g The proportion of organic pores in the organic matter mixture, a variable. ζ It has the following forms of expression:

[0172] Formula 6: .

[0173] In step A of this embodiment, the self-compatible approximation (SCA) theory includes the following:

[0174] Given the elastic properties and geometry of the components of the equivalent medium, the self-consistent approximation (SCA) theory can accurately calculate the elastic modulus of the equivalent medium. In this invention, the bulk modulus and shear modulus of the porous matrix can be obtained by adding water-saturated inorganic pores to the solid matrix. For two-phase mixtures, Wu (1966) gave an estimation formula for its elastic modulus:

[0175] Formula 7: ;

[0176] in, K matrix and μ matrix These represent the bulk modulus and shear modulus of the porous matrix, respectively. K 0 and μ 0 represents the bulk modulus and shear modulus of the solid matrix, respectively; K w and μ w These are the bulk modulus and shear modulus of water, respectively. x w The inorganic pores constitute the volume fraction of the porous matrix; the shape of the inorganic pores can be determined by a factor. P 0w and G 0wwhere its specific expression can be found in Mavko et al., (2009); for more complex N-phase mixtures, Berryman (1980) gives a general form of self-consistent approximation for it in future studies:

[0177] Equation 8: ;

[0178] where i is the volume modulus of the i-th material, i K i and μ i are the bulk modulus and shear modulus of the i-th material, respectively, i x i is its volume content, P 0w and G 0w denote the shape of the inclusion, and the superscript i refers to the geometric factor of the inclusion i with respect to the background medium in the calculation process.

[0179] In step A of the embodiment, the Krief formula includes the following content:

[0180] Krief gives the relationship between the elastic modulus of dry rock and solid matrix. In the present application, the Krief formula is used to calculate the bulk modulus K frame and shear modulus μ frame of the solid skeleton:

[0181] Equation 9: ;

[0182] where K matrix and μ matrix are the bulk modulus and shear modulus of the pore matrix, respectively; f om is the volume fraction of the organic matter mixture composed of kerogen and organic pores in the entire rock; m is an empirical parameter related to the rock facies, and the present application takes a value of 5.

[0183] In step A of the embodiment, the Ciz-Shapiro solid replacement theory includes the following content:

[0184] ​​An organic mixture consisting of kerogen and organic pores developed inside it has a non-zero bulk and shear modulus (Carcione, 2000; Zhao et al., 2016), so the Gassmann equation, which is valid for fluids (zero shear modulus), is no longer applicable. For the case of pores filled with solid material, Ciz and Shapiro (2007) proposed a solid substitution theory to calculate the elastic moduli of saturated rocks. For the isotropic case, the expressions for the bulk modulus K sat and shear modulus μ sat of a saturated rock, i.e., an organic-rich shale, are:

[0185] Equation 10: ;

[0186] Equation 11: ;

[0187] where f om is the volume fraction occupied by the organic mixture in the entire rock; K mix and μ mix are the bulk and shear moduli of the organic mixture;

[0188] Equations 10 and 11 are further simplified in the form:

[0189] Equation 12: ;

[0190] Equation 13: ;

[0191] where, K om and μ om are the bulk and shear moduli associated with the organic matter;

[0192] Thus, the P-wave modulus of an organic-rich shale and its solid matrix are:

[0193] Equation 14: ;

[0194] Equation 15: ;

[0195] where, M sat and M frame are the P-wave moduli of the organic-rich shale and the solid matrix, respectively, μ sat and μframe Gshear for the organic-rich shale and the solid matrix, respectively;

[0196] Substituting equations 12 and 13 into equation 14 gives:

[0197] Equation 16: ;

[0198] Further substituting equations 14 and 15 into equation 16 gives:

[0199] Equation 17: ;

[0200] where K om and μ om are the bulk modulus and shear modulus associated with the organic matter, M om = K om + 4 / 3 μ om denotes the P-wave modulus associated with the organic matter.

[0201] In step B of the present embodiment, the inversion method in the well includes:

[0202] Two key variables are applied in the petrophysical inversion process: the volume fraction of the organic matter (fom) and the aspect ratio of the inorganic pores (a); the optimal values of fom and a are obtained by fitting the measured and petrophysically calculated P- and S-wave velocities, and minimizing the objective function. f om α f om and the aspect ratio of the inorganic pores α The objective function F obj is as follows:

[0203] Equation 18: ;

[0204] where, V P and V S denote the P- and S-wave velocities, respectively; the subscripts “modeled” and “measured” denote the petrophysically calculated and measured velocity values, respectively; in the petrophysical modeling process, the total porosity φ is set to be the sum of the inorganic pore porosity φ inorganic and the organic pore porosity φ organic , i.e. φ = φ ​​organic + φ inorganic ;

[0205] And, the predicted f om And α The numerical values are taken as inputs of the rock physics model to calculate the elastic modulus of the solid skeleton in the well (E M frame With μ frame ) and the elastic modulus of the organic matter (E M om With μ om ) in the well, the rock physics inversion results provide an initial model with physical meaning for subsequent seismic rock physics inversion.

[0206] In step C of the embodiment, the content of the reflection coefficient formula linearly characterized by the shale solid skeleton elastic modulus (E M frame , μ frame ), the organic matter elastic modulus (E M om , μ om ) and the density parameter (D ρ sat ) includes:

[0207] Shaw and Sen (2004, 2006) derived the PP wave reflection coefficient on the interface of any weak anisotropic medium by combining the perturbation of the saturated rock stiffness matrix and the Born approximation: R PP The relationship with the scattering function S is:

[0208] Formula 19: ;

[0209] Wherein θ and ψ represent the incidence angle and the azimuth angle, ρ sat respectively, and D represents the density of the saturated rock, and r0 is the point satisfying the Snell law on the horizontal interface; the scattering function S is related to the perturbation of the elastic parameters, the horizontal slowness and the polarization direction, and its expression is:

[0210] Formula 20: ;

[0211] Wherein S (r0) is the scattering function, and D ρ sat is the perturbation of the saturated rock density, and D Cijkl is the perturbation of the stiffness tensor of the saturated rock, variable ξ is a physical quantity related to the polarization vector of the PP wave, while variable η ijkl is a physical quantity related to both the polarization vector of the PP wave and the slowness vector;

[0212] According to Shaw and Sen (2006), the dummy index “ ijkl ” in equation 20 is simplified to the dummy index “ IJ ” according to the following equation:

[0213] Equation 21: ;

[0214] Equation 22: ;

[0215] where ijkl represents the subscript value of the fourth-order stiffness tensor, i , j , k , l , and IJ , respectively, take values from 1 to 3; I represents the subscript value of the simplified fourth-order stiffness tensor, J , , respectively, take values from 1 to 6.

[0216] δ where ij and δ kl are the Kronecker symbols, therefore, equation 20 can be simplified as:

[0217] Equation 23: ;

[0218] where S (r0) is the scattering function, Δ ρ sat is the perturbation of the density of the saturated rock, Δ C ijkl is the perturbation of the stiffness tensor of the saturated rock, variable ξ is a physical quantity related to the polarization vector of the PP wave, while variable η ijkl is a physical quantity related to both the polarization vector of the PP wave and the slowness vector

[0219] the polarization vector of the PP wave t and t ’ can be represented by the incident wave and the reflected wave, respectively:

[0220] Equation 24: ;

[0221] Equation 25: ;

[0222] in θ and ψ These represent the angle of incidence and the azimuth, respectively.

[0223] Its corresponding slowness vector p and p Then, represented by the polarization vector, it is:

[0224] Formula 26: ;

[0225] Formula 27: ;

[0226] in α sat The longitudinal wave velocity of saturated rock;

[0227] When the PP wave is incident and reflected, the variable ξ The expression is written as:

[0228] Formula 28: ;

[0229] t i and They represent polarization vectors respectively. t and t 'the i Each component.

[0230] And variables η IJ The expression is written as:

[0231] Formula 29: ;

[0232] In the formula, Represents polarization vector t 'the i One portion, Represents the slowness vector p The j-th component of ', tk represents the k-th component of the polarization vector t, p l The l-th component variable represents the slowness vector p; η IJ The expression for each component is:

[0233] Formula 30: ;

[0234] η 11 , η 12 , η 13 ,η 22 , η 23 , η 33 , η 44 , η 55 , η 66 , η 21 , η 31 , η 32 represent the perturbations of the variable η IJ in different directions.

[0235] For a weakly isotropic medium, the perturbations Δ C IJ in the stiffness matrix in equation 23 can be simplified to the form

[0236] Equation 31: ;

[0237] where Δ M sat , Δ λ sat and Δ μ sat are the perturbations of the saturated rock's longitudinal wave modulus M sat , the first Lamé constant λ sat and the second Lamé constant μ sat respectively;

[0238] Substituting equations 28, 30 and 31 into equation 23, we obtain the scattering function S for a weakly isotropic medium, which is given by

[0239] Equation 32: ;

[0240] where Δ M sat , Δ λ sat and Δ μ sat are the perturbations of the saturated rock's longitudinal wave modulus M sat , the first Lamé constant λ sat and the second Lamé constant μ sat respectively;

[0241] Finally, the PP-wave reflection coefficient approximation formula for any weakly anisotropic medium can be obtained by substituting equation 32 into equation 19.

[0242] For the petrophysical model, if the media on both sides of the interface satisfy the weak reflection approximation condition, then the perturbation of the P-wave modulus and shear modulus of the organic-rich shale is:

[0243] Equation 33:

[0244] Equation 34:

[0245] where Δ M frame and Δ μ frame represent the perturbation of the P-wave modulus and shear modulus of the solid skeleton, respectively, and Δ M om and Δ μ om represent the perturbation of the corresponding elastic parameters of the organic matter, respectively, and the perturbations satisfy the following conditions: |Δ M frame | M frame , |Δ μ frame | μ frame , |Δ M om | M om , |Δ μ om | μ om

[0246] Therefore, the PP-wave linearized reflection coefficient approximation formula is obtained by substituting equations 32, 33, and 34 into equation 19 as follows:

[0247] Equation 35:

[0248] where R PP is the reflection coefficient, θ is the incidence angle, Δ ρ sat is the perturbation of the saturated rock density.

[0249] Subsequently, equations 13 and 17 are substituted into equation 35 to obtain:

[0250] Equation 36:

[0251] After that, define γ ​​​​​2 sat = M sat / μ sat and γ 2 frame = M frame / μ frame , respectively, represent the square of the ratio of the P- and S-wave velocities of the organic-rich shale and the solid matrix, thus, the above equation can be simplified as

[0252] Equation 37: ;

[0253] Finally, define the parameter p = μ frame / μ sat as the ratio of the shear modulus of the solid matrix and the organic-rich shale, thus, Equation 37 can be further simplified and the final linearized reflectivity equation is obtained as:

[0254] Equation 38: ;

[0255] For simplicity of the equation expression and convenience of the subsequent inversion study, Equation 38 is simplified as:

[0256] Equation 39: ;

[0257] The coefficient term in the equation is defined as:

[0258] Equation 40: ;

[0259] The above equation represents the case of the organic matter with non-zero shear modulus as the solid filler in the organic-rich shale; for the case of fluid filling, i.e., the remaining pore filling fluid in Figure 2 , it satisfies the condition of p = 1 and satisfies the Gassmann equation; under this condition, the shear modulus of the solid matrix is equal to the shear modulus of the whole rock, i.e., μ frame = μ sat , thus p = 1 and μ frame = μ sat are brought into Equation 38 to obtain the following equation:

[0260] Equation 41: ;

[0261] According to the poroelastic medium theory proposed by Russell et al. (2011), the Gassmann equation describes the relationship between the longitudinal wave modulus of the solid matrix and the shear modulus of the saturated rock, i.e.: M frame The relationship between the shear modulus of the saturated rock and the bulk modulus of the solid matrix is: μ sat i.e.:

[0262] Equation 42: ;

[0263] After calculating the reflectivity of the variable Mframe, we have:

[0264] Equation 43: ;

[0265] Therefore, Equation 41 is further simplified as:

[0266] Equation 44: ;

[0267] According to the reflectivity formula derived by Russell et al.:

[0268] Equation 45: ;

[0269] Equation 44 and Equation 45 have the same form, where Δ M om / M om is equivalent to Δ f / f in Equation 45.

[0270] In step E of this embodiment, the content of the seismic rock physics inversion method includes:

[0271] The reflectivity of the elastic modulus of the solid matrix Δ M frame / M frame and Δ μ frame / μ frame , the reflectivity of the elastic modulus of the organic matter Δ M om / M om and Δ μ om / μ om , and the reflectivity of the density of the organic-rich shale Δ ρ sat / ρ sat can be approximately expressed by their corresponding natural logarithm values, i.e.:

[0272] Equation 46: ;

[0273] where the symbol "ln( )" denotes the natural logarithm of the elastic parameter, and the symbol "Δ" denotes the difference in the value on either side of the reflecting interface;

[0274] In addition, the elastic properties of the solid matrix, the density of the organic-rich shale, and the organic-matter-related elastic properties are extended to continuous time variables in logarithmic form as follows:

[0275] Equation 47: ;

[0276] Thus, the reflection coefficient of a single interface represented by Equation 37 is extended to a continuous-time reflection coefficient equation using Equations 45 and 46:

[0277] Equation 48: ;

[0278] where the symbol "d" denotes the difference in the elastic parameter between the tth time and the (t+1)th time in continuous time; for example, dln[ M frame ( t )] = ln[ M frame ( t +1)] - ln[ M frame ( t )].

[0279] Subsequently, the length of the seismic data in time is set to n+1, and the number of incidence angles is h, so there are n reflecting interfaces in the time direction of the seismic data, and thus the angle-dependent seismic wavelet matrix W is represented as:

[0280] Equation 49: ;

[0281] The block matrix W i in the matrix W contains the seismic wavelet related to the incidence angle θ i , and is expressed as:

[0282] Equation 50: ;

[0283] where w i denotes the seismic wavelet when the incidence angle is θ i , i = 1, 2,..., h. s 1, s 2,..., s lThe time sampling points represent the seismic wavelet;

[0284] The seismic wavelet matrix W is convolved with the reflection coefficient sequence to obtain the seismic forward modeling equation with data error terms, as shown below:

[0285] Formula 51: ;

[0286] In the formula, S represents the seismic response, W is the seismic wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, and e represents the data residual term, which satisfies the condition that the mean is 0 and the variance is Σ. d The Gaussian distribution is given by m, which is the model parameter vector, and its expression is:

[0287] Formula 52: ;

[0288] Where the superscript T denotes the transpose of the vector, and the length of vector m is... q = 5( n +1), this is because each of the n+1 strata contains five elastic parameters, which ultimately form n reflection points;

[0289] The coefficient matrix Λ has a dimension of (h×n)×5n, and its expression is:

[0290] Formula 53: ;

[0291] in A ( θ i ), B ( θ i ), C ( θ i ), D ( θ i )and E ( θ i ) is a diagonal matrix, and its expression is:

[0292] Formula 54: ;

[0293] in a ( t j , θ i ), b ( t j , θ i ), c ( t j ,θ i ), d t j , θ i ) and e t j , θ i ) are angle and time dependent coefficients, where j = 1, 2, …, n and i = 1, 2, …, h ;

[0294] The difference matrix D has a dimension of 5n x 5(n+1) and its expression is:

[0295] Equation 55: ;

[0296] where the block matrix Dd has a dimension of n x (n+1) and its expression is:

[0297] Equation 56: ;

[0298] According to the Bayesian theory, the posterior probability density function P(m|d) is obtained from the prior probability density function P(d|m) and the likelihood function P(m):

[0299] Equation 57: ;

[0300] In the present application, the prior model parameter vector m is set to satisfy a multi-dimensional Gaussian distribution with a mean of μ m and a covariance of Σ m , and the data residual term satisfies a Gaussian distribution with a mean of 0 and a covariance of Σ d , therefore, the posterior probability density function P(m|d) satisfies a Gaussian distribution, and the expression of its mean μ m|d and covariance Σ m|d is:

[0301] Equation 58: ;

[0302] Equation 59: ;

[0303] where the matrix G = WΛD is the forward operator;

[0304] Therefore, the expression of the maximum likelihood solution m opt of the inversion parameter is:

[0305] Equation 60: ​​;

[0306] Where μ m|d and Σ m|d Let P(m|d) be the mean and covariance of the posterior probability density function.

[0307] Experimental Example 1, Rock Physical Model Test:

[0308] Figure 3 The simulated saturated P-wave modulus and saturated shear modulus of organic-rich shale are shown. M sat , μ sat ), and the longitudinal wave modulus and shear modulus of its solid skeleton ( M frame , μ frame ) f om and α The simulation results of the changes, among which f om The range of variation is 0-0.1. α The range of variation is 10 -1.5 – 1. For ease of representation, f om The x-axis is... α Take the common logarithm, i.e., lg( α i () is the ordinate. Figure 3 a and b show that as the organic matter content increases or the aspect ratio of inorganic pores decreases, M sat and μ sat It exhibits a significant non-linear decreasing trend. This suggests that higher organic matter content or narrower inorganic pores reduce the stiffness of organic-rich shale. Figure 3 c and d demonstrate M frame and μ frame A similar trend was observed. However, under the same organic matter content and inorganic porosity ratio, the elastic modulus of the solid skeleton was consistently lower than that of organic-rich shale. These results confirm the influence of organic matter on shear modulus and P-wave modulus.

[0309] Using the same simulation results, the organic matter-related longitudinal wave modulus and shear modulus were calculated using formulas 17 and 13. M om , μ om The result is as follows Figure 4 As shown in figures a and b. It can be seen from the figure that... M om and μom Both increase with the increase of f om and α i , and are very sensitive to the change of f om . Therefore, in order to highlight this dependence, Figure 4 c and d are converted to be plotted with f om as the horizontal coordinate, and the elastic modulus M om and μ om as the vertical coordinate, respectively, while lg( α ) as the color bar. The simulation results show that M om and μ om present a very significant approximate linear positive correlation with f om , which proves their potential as indicators directly representing the enrichment degree of organic matter in shale gas reservoirs.

[0310] In order to verify the accuracy and precision of the proposed reflection coefficient formula 38, a seismic reflection model containing three shale layers is designed according to the geological background and logging data of the study area, as shown in Figure 5 .

[0311] Figure 6 The reflection coefficients of the top and bottom interfaces calculated by the Zoeppritz equation (1919), the Aki-Richards approximate formula (2002) and the reflection coefficient approximate formula proposed in the present application are shown. The numerical results show that the reflection coefficients calculated by the proposed reflection coefficient approximate formula are highly consistent with the calculation results of the traditional Zoeppritz equation and the Aki-Richards approximate formula, and only when the incident angle is greater than about 35°, the reflection coefficients calculated by the present method appear a small deviation compared with the more accurate Zoeppritz equation.

[0312] The above research results confirm the reliability of the proposed reflection coefficient approximate formula in simulating the PP wave reflection characteristics, and verify its applicability in predicting the elastic properties of solid skeleton and organic matter respectively in seismic rock physics inversion.

[0313] Figure 7 The contribution of each elastic parameter term to the reflection coefficient is shown, and the parameter change range is set to-0.3-0.3. As shown in Figure 7 a- Figure 7 c, the longitudinal wave modulus reflectivity (Δ M frame / of the solid skeleton Mframe ) and the density reflectance (Δ) of organic-rich shale ρ sat / ρ sat The shear modulus reflectance of the solid skeleton contributes to the reflection coefficient of the PP wave across all incident angles. Conversely, the shear modulus reflectance of the solid skeleton (Δ...) contributes to the reflection coefficient of the PP wave across all incident angles. μ frame / μ frame Then, the reflection coefficient begins to make a significant contribution when the incident angle is about 10°. Figure 7 d and e show the longitudinal wave modulus reflectivity (Δ) related to organic matter. M om / M om ) and shear modulus reflectance (Δ μ om / μ om The contribution of organic matter-related elastic modulus reflectance to the reflection coefficient was investigated. The results showed that the contribution of organic matter-related elastic modulus reflectance to the reflection coefficient exhibited a similar trend to that of the corresponding solid framework elastic parameters, but the contribution of organic matter-related elastic modulus reflectance to the reflection coefficient was relatively small. This phenomenon is consistent with the prediction of rock physics theory, namely that organic matter, as inclusions, has a weaker influence on the overall elastic parameters of shale compared to the solid framework, and consequently a relatively weaker influence on the reflection coefficient.

[0314] Experimental Example 2, Rock Physical Inversion in Wells:

[0315] Figure 8 The flowchart of forward and inverse rock physics modeling in the well is shown.

[0316] Figure 9 The logging curves of well A in the study area are shown, including natural gamma (GR) and P-wave velocity (P-wave velocity). V P ), transverse wave velocity ( V S ),density( ρ ), total porosity ( φ (The image shows) and mineral volume fraction curves. The shale layer ranges from approximately 2670 to 2825 meters, with the two black lines representing the depth range of the target layer.

[0317] Subsequently, the solid skeleton and organic matter elastic modulus were inverted from the logging data of Well A using the in-well inversion method. Figure 10 The model-driven rock physics inversion results are presented. For example... Figure 10 The calculated P-wave and S-wave velocities (red lines) from rock physics calculations (a and b) are in high agreement with the measured P-wave and S-wave velocities (black lines), which verifies the effectiveness of the rock physics model for organic-rich shale and the corresponding rock physics inversion method proposed in this study.Figure 10 c and d show the inversion results in the well for parameters f om and α respectively. Based on the inversion results, the shale organic porosity (Porg) is calculated using the formula φ organic = f om - V kerogen , as shown in φ organic e. Figure 10

[0318] According to the petrophysical forward and inversion procedure in Figure 8 , the inversion results of f om and α will be used as the input of the petrophysical model to calculate the elastic modulus (E and μ) from the well logging data, as shown in M frame , μ frame , M om , μ om . Figure 11 Figure 11 a and b show the P-wave modulus and shear modulus of the organic-rich shale and the solid matrix, respectively, represented by the black and red lines. Figure 11 c shows the density of the organic-rich shale (black line) and the corresponding density of the solid matrix (red line). Finally, the organic matter related P-wave modulus and shear modulus are calculated, as shown in Figure 11 d and e.

[0319] Example 3, synthetic data test:

[0320] Based on the petrophysical forward and inversion results in Figure 12 , first, the time-depth conversion is performed to convert the well logging curves in the depth domain to the time domain. Then, the time domain well logging curves are resampled at 2 ms interval to match the time sampling interval of the actual seismic data. The real model in the time domain is shown in Figure 12 a- Figure 12 e. Figure 12 f shows the noise-free synthetic seismic record generated using the formula 50 and the actual seismic wavelet.

[0321] Subsequently, the seismic petrophysical inversion method proposed by the present application will be applied to the synthetic seismic data as shown in Figure 12 f. The inversion results are shown in Figure 13 a Figure 13 ​​e, where the black solid line represents the true model, the blue solid line represents the initial model, and the red dashed line represents the inversion result. Figure 13 f shows the seismic record (red) obtained by using the inversion result for seismic forward modeling, compared with Figure 12 the noise-free synthetic seismic record in f. The high agreement between the inversion result and the true model, and the consistency of the corresponding seismic response, verify the effectiveness and robustness of the inversion method and result.

[0322] Experimental Example 4, application of actual data:

[0323] Figure 14 is the post-stack seismic data of Well A.

[0324] Figure 15 is the actual well logging data (Vsh, GR, and DEN) of Well A. Figure 11 ) is the initial model obtained by using the Kriging interpolation method to interpolate and extrapolate the actual well logging data (Vsh, GR, and DEN) of Well A, and performing low-pass filtering on the interpolation result.

[0325] Based on the initial model of Figure 15 , the post-stack seismic data, and the actual seismic wavelet, the elastic parameters M frame , μ frame , ρ sat , M om and μ om , as shown in Figure 16 . According to the petrophysical meaning of the proposed petrophysical model, Figure 16 a and b, which represent the elastic modulus of the solid skeleton (Ea and Eb) M frame , μ frame , can directly reflect the elastic modulus of the shale pore matrix and eliminate the influence of organic matter. These moduli are mainly affected by the mineral composition of shale, especially the relative content of clay, quartz, and carbonate minerals in the shale in the study area. This indicates that, through appropriate calibration of comprehensive well logging data, M frame and μ frame can be used as quantitative indicators for evaluating the brittleness index, which is crucial for optimizing the hydraulic fracturing effect of shale gas reservoirs. At the same time, through the petrophysical and well logging data analysis in the above, the elastic modulus related to organic matter (Eo) M om , μ om ​) can be used as a key indicator of the quantitative characterization of the occurrence of organic matter. The inverted parameters can provide an effective method for the evaluation of the brittleness of shale gas reservoir engineering sweet spots and the prediction of the enrichment of organic matter of geological sweet spots.

[0326] The above merely describes the preferred embodiments of the present application, and does not limit the patent scope of the present application, and any equivalent structure or equivalent process transformation made by the present application, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present application.

Claims

1. A method for simultaneous prediction of shale solid matrix and organic matter elastic parameters, characterized in that, The method comprises: Step A, on the basis of shale gas reservoir microstructure analysis, a rock physics model is constructed to characterize the elastic modulus of the solid skeleton and the elastic modulus of the organic matter in the shale reservoir; Step B, based on the rock physics model constructed in step A, a well inversion method is developed to predict the elastic modulus of the solid skeleton and the elastic modulus of the organic matter from conventional logging data, and to provide an initial model with physical meaning for actual pre-stack seismic data inversion; Step C, in combination with the rock physics model constructed in step A and the scattering function theory, a reflection coefficient formula is established, which is linearly characterized by the elastic modulus of the shale solid skeleton, the elastic modulus of the organic matter and the density parameter; Step D, precision evaluation and sensitivity analysis of the parameters in the reflection coefficient formula are carried out to verify the applicability of the reflection coefficient formula; Step E, in combination with the reflection coefficient formula established in step C and the Bayesian linear inversion theory, a seismic rock physics inversion method is developed to directly and synchronously predict the elastic modulus of the shale solid skeleton and the elastic modulus of the organic matter from pre-stack seismic data; Step F, the solid skeleton elastic modulus and the organic matter elastic modulus predicted in step B are applied to calculate the synthetic seismic data in combination with the actual seismic wavelet, and the effectiveness of the seismic rock physics inversion method constructed in step E is tested and verified; Step G, the seismic rock physics inversion method developed in step E and the initial model constructed in step B are applied to actual pre-stack seismic data to calculate the elastic modulus of the solid skeleton and the organic matter in the shale gas reservoir in the study area, and a new effective method is provided for the brittle evaluation of shale gas reservoir engineering sweet spots and the prediction of organic matter enrichment in geological sweet spots.

2. The method of claim 1, wherein, In step A, the steps of constructing the rock physics model to characterize the elastic modulus of the solid skeleton and the elastic modulus of the organic matter in the shale reservoir include: The Hashin-Shtrikman boundary theory is applied to calculate the bulk modulus and shear modulus of the solid matrix composed of complex minerals; The self-consistent approximation theory is used to calculate the bulk modulus and shear modulus of the pore matrix composed of the solid matrix and water-saturated inorganic pores; The Krief formula is used to calculate the bulk modulus and shear modulus of the solid skeleton; The Kuster-Toksöz (KT) formula is used to calculate the bulk modulus and shear modulus of the organic matter composed of kerogen and gas-saturated organic pores; The overall bulk modulus and shear modulus of the organic-rich shale composed of the solid skeleton and the organic matter can be obtained by using the Ciz-Shapiro solid replacement theory.

3. The method of claim 2, wherein, In step A, the Hashin-Shtrikman boundary theory includes the following contents: For a medium containing two different components, the equivalent bulk modulus and shear modulus are calculated by the HSB theory as follows: Formula 1: ; where K HS± and μ HS± represent the upper and lower bounds, respectively, of the equivalent medium bulk modulus; K i and μ i represent the elastic modulus of the i-th constituent; f i represents the volume fraction of the i-th constituent; For an equivalent medium composed of multiple components, the equivalent elastic modulus is calculated by the following formula: Formula 2: ; where K HS+ and K HS- represent the upper and lower bounds of the effective medium bulk modulus, μ HS+ and μ HS- represent the upper and lower bounds of the effective medium shear modulus, K max and μ max represent the maximum values of the bulk modulus and shear modulus of the components of the effective medium, K min and μ min represent the minimum values of the bulk modulus and shear modulus of the components of the effective medium; Formula 3: ; where z represents the independent variable in the function, z = K when calculating the upper limit max or μ max , z = K when calculating the lower limit min or μ min ; Λ(z), Γ(z) and ζ(K, μ) represent functions in the calculation process, K represents the bulk modulus and shear modulus; r represents the rth component in the rock, the value is related to the rock physics modeling process; the symbol <·> represents the volume weighted average of each component. The bulk modulus K0 and shear modulus μ0 of the equivalent medium are calculated by the average value of the upper and lower limits of the HSB theory, i.e. Formula 4: .

4. The method of claim 2, wherein, In step A, the Kuster-Toksöz formula includes the following contents: The gas-saturated organic pores in the kerogen form an organic mixture, and the expression of its elastic modulus is as follows: Formula 5: ; where K mix and μ mix represent the bulk and shear moduli of the organic mixture; K k and μ k represent the bulk and shear moduli of the kerogen; K g and μ g represent the bulk and shear moduli of the gas; the shape of the organic pores is represented by the structure factors P kg and Q kg ; x g is the fraction of organic pores in the organic mixture, and the variable ζ k has the following expression: Equation 6: .

5. The method of claim 2, wherein, In step A, the self-consistent approximation theory includes the following contents: The bulk and shear moduli of the porous matrix are obtained by adding water-saturated inorganic pores to the solid matrix, and the formula for estimating the elastic modulus of the two-phase mixture is: Equation 7: ; where K matrix and μ matrix represent the bulk and shear moduli of the pore matrix; K0and μ0are the bulk and shear moduli of the solid matrix; K w and μ w are the bulk and shear moduli of water; x w is the volume fraction of inorganic pores in the pore matrix; the shape of the inorganic pores is represented by the structure factors P 0w and Q 0w ; for more complex N-phase mixtures, the self-consistent approximation takes the form: Equation 8: ; where i is the i-th material, K i and μ i are the bulk and shear moduli of the i-th material, x i is its volume content, P 0w and Q 0w denote the shape of the inclusion, and the superscript *i refers to the geometric factor of inclusion i with respect to the background medium in the calculation.

6. The method of claim 2, wherein, In Step A, the Krief formula includes the following: The bulk modulus K of the solid skeleton is calculated using the Krief formula frame and the shear modulus μ frame : Equation 9: ; where K matrix and μ matrix represent the bulk modulus and the shear modulus of the pore matrix, respectively; f om is the volume fraction of the organic matter mixture consisting of kerogen and organic pores in the entire rock; and m is an empirical parameter related to the rock facies.

7. The method of claim 2, wherein, In Step A, the Ciz-Shapiro solid substitution theory includes the following: For the isotropic case, the expressions for the bulk modulus K sat and shear modulus μ sat of the saturated rock, i.e. the organic-rich shale, are given by Equation 10: ; Equation 11: ; where f om is the volume fraction of the organic-matter mixture in the entire rock; K mix and μ mix are the bulk modulus and shear modulus of the organic-matter mixture; Formulas 10 and 11 are simplified as follows: Equation 12: ; Equation 13: ; where K om and μ om are the bulk modulus and shear modulus associated with the organic matter; Therefore, the longitudinal wave modulus of the organic-rich shale and its solid skeleton is: Equation 14: ; Equation 15: ; where M sat and M frame are the longitudinal wave modulus of the organic-rich shale and the solid matrix, respectively, μ sat and μ frame are the shear modulus of the organic-rich shale and the solid matrix, respectively; Bringing formulas 12 and 13 into formula 14 gives: Equation 16: ; Bringing formulas 14 and 15 into formula 16 gives: Equation 17: ; where K om and μ om are the bulk modulus and shear modulus associated with the organic matter, M om = K om + 4 / 3 μ om represents the longitudinal wave modulus associated with the organic matter.

8. The method of claim 1, wherein, In Step B, the well inversion method includes: By fitting the measured and petrophysical model calculated P-wave and S-wave velocities, the objective function is minimized, and the optimal volume fraction of organic matter f is obtained om and the aspect ratio a value of inorganic pores, the objective function F obj As follows: Equation 18: ; where V P and V S represent the longitudinal and transverse wave velocities, respectively; the subscripts "modeled" and "measured" represent the velocities calculated from the petrophysical model and the measured velocities, respectively. In the petrophysical modeling process, the total porosity φ is set to be composed of the inorganic porosity φ inorganic and the organic porosity φ organic , i.e., φ = φ organic + φ inorganic . And, the predicted f om And the values of a and a are used as input of the rock physics model to calculate the elastic modulus of the solid matrix and the elastic modulus of the organic matter in the well. The rock physics inversion results provide an initial model with physical meaning for the subsequent seismic rock physics inversion.

9. The method of claim 1, wherein, In Step C, the content of establishing the reflectance formula characterized by the linearization of the elastic modulus of the shale solid skeleton, the elastic modulus of the organic matter, and the density parameter includes: PP-wave reflection coefficient R on an arbitrary weakly anisotropic medium interface PP Relationship with scattering function S: Equation 19: ; where θ and ψ represent the incidence angle and azimuth, respectively, and p sat represents the density of the saturated rock, and r0 is the point on the horizontal interface that satisfies Snell's law; the scattering function is related to the perturbation of the elastic parameters, the horizontal slowness, and the polarization direction, and its expression is: Equation 20: ; where S(r0) is the scattering function, Δρ sat is the perturbation of the saturated rock density, ΔC ijkl is the perturbation of the saturated rock stiffness tensor, the variable ξ is a physical quantity related to the polarization vector of the PP wave, and the variable η ijkl is a physical quantity related to the polarization vector of the PP wave and the slowness vector; The dummy index "ijkl" in formula 20 is simplified to the dummy index "IJ" according to the following formula: Equation 21: ; Equation 22: ; Whereijklrepresents the index values of the fourth-order stiffness tensor, and the values of i, j, k, and l are 1 to 3; IJrepresents the index values of the simplified fourth-order stiffness tensor, and the values of I and J are 1 to 6; where δ ij and δ kl is the Kronecker symbol, so that equation 20 simplifies to: Equation 23: ; where S(r0) is the scattering function, Δρ sat is the perturbation of the saturated rock density, ΔC ijkl is the perturbation of the saturated rock stiffness tensor, the variable ξ is a physical quantity related to the polarization vector of the PP wave, and the variable η ijkl is a physical quantity related to the polarization vector of the PP wave and the slowness vector; The polarization vectors t and t' of the PP wave can be represented by the incident wave and the reflected wave, respectively, as follows: Equation 24: ; Equation 25: ; Where θ and ψ represent the incident angle and the azimuth angle, respectively; The corresponding slowness vectors p and p' are represented by the polarization vectors as follows: Equation 26: ; Equation 27: ; where a sat is the P-wave velocity of the saturated rock; When the PP wave is incident and reflected, the expression of the variable ξ is written as: Equation 28: ; t i and ti and t'j represent the i-th component of the polarization vectors t and t', respectively; and the variable η IJ The expression for η is written as: Equation 29: ; wherein represents the i-th component of the polarization vector t', represents the j-th component of the slowness vector p', tkrepresents the k-th component of the polarization vector t, p l represents the i-th component of the polarization vector t', IJ The expression for each component is: Equation 30: ; η 11 , η 12 , η 13 , η 22 , η 23 , η 33 , η 44 , η 55 , η 66 , η 21 , η 31 , η 32 represent the components of the variable η IJ in different directions; For weakly isotropic media, the perturbation of the stiffness matrix ΔC in equation 23 IJ is expressed in the form of simplification as follows: Equation 31: ; where ΔM sat , Δλ sat and Δμ sat are the perturbations of the saturated rock's longitudinal wave modulus M sat , the first Lamé constant λ sat and the second Lamé constant μ sat , respectively; Bringing formulas 28, 30, and 31 into formula 23 gives the scattering function of the weakly isotropic medium, which is expressed as: Equation 32: ; where ΔM sat , Δλ sat , and Δμ sat are the perturbations in the saturated rock's longitudinal wave modulus M sat , the first Lamé constant λ sat , and the second Lamé constant μ sat , respectively; Finally, bringing formula 32 into formula 19 gives the PP wave reflectance approximation formula for any weakly isotropic medium; For the rock physics model, the media on both sides of the interface satisfy the weak reflection approximation condition, and the perturbation of the longitudinal wave modulus and the shear modulus of the organic-rich shale is: Equation 33: ; Equation 34: ; where ΔM frame and Δμ frame represent the perturbation of the longitudinal modulus and shear modulus of the solid skeleton, respectively, ΔM om and Δμ om represent the perturbation of the corresponding elastic parameters of the organic matter, respectively, and the perturbations satisfy the following conditions: |ΔM frame |M frame , |Δμ frame |μ frame , |ΔM om |M om , |Δμ om |μ om ; Therefore, bringing formulas 32, 33, and 34 into formula 19 gives the PP wave linearization reflectance approximation formula as follows: Equation 35: ; where R PP is the reflection coefficient, θ is the angle of incidence, Δρ sat is the perturbation of the saturated rock density; Subsequently, bringing formulas 13 and 17 into formula 35 gives: Equation 36: ; After that, define γ2 sat = M sat / μ sat and γ2 frame = M frame / μ frame , which represent the square of the ratio of the P- and S-wave velocities of the organic-rich shale and the solid matrix, respectively, so the above equation can be simplified as Equation 37: ; Finally, define parameter p = μ frame / μ sat is the ratio of the shear modulus of the solid matrix and the organic-rich shale, so formula 37 can be simplified and the final linearized characterization of the reflectivity formula is obtained as follows: Equation 38: ; In order to simplify the formula expression and facilitate subsequent inversion research, formula 38 is simplified as follows: Equation 39: ; The coefficient term in the formula is defined as: Equation 40: ; Taking p = 1 and μ frame = μ sat into equation 38 gives the following equation: Formula 41: ; The Gassmann equation describes the longitudinal wave modulus M of a solid skeleton. frame With saturated rock shear modulus μ sat The relationship is as follows: Equation 42: ; After that, the reflectivity of the variable Mframe is calculated, and the following is obtained: Equation 43: ; Therefore, formula 41 is simplified as follows: Equation 44: ; The reflectance formula: Equation 45: ; Equation 44 has the same form as Equation 45, where AM om / M om is equivalent to Δf / f in Equation 45.

10. The method of claim 1, wherein, In Step E, the content of the seismic rock physics inversion method includes: The relative quantity of the linearization characterized reflectance formula is approximately represented by its corresponding natural logarithm value, that is: Equation 46: ; In the formula, the symbol "ln(·)" represents the natural logarithm of the elastic parameter, and the symbol "Δ" represents the difference between the values on both sides of the reflection interface; In addition, the elastic properties of the solid skeleton, the density of the organic-rich shale, and the elastic properties related to the organic matter expressed in the logarithmic form are expanded to continuous time variables as follows: Equation 47: ; Therefore, the reflectance of a single interface represented by formula 37 is expanded to the reflectance equation of continuous time by formulas 45 and 46: Equation 48: ; In the formula, the symbol "d" represents the difference between the elastic parameters at the t time and the t+1 time on the continuous time; The polarization vectors t and t' of the PP wave can be represented by the incident wave and the reflected wave, respectively, as follows: Where θ and ψ represent the incident angle and the azimuth angle, respectively; The corresponding slowness vectors p and p' are represented by the polarization vectors as follows: When the PP wave is incident and reflected, the expression of the variable ξ is written as: Bringing formulas 28, 30, and 31 into formula 23 gives the scattering function of the weakly isotropic medium, which is expressed as: Finally, bringing formula 32 into formula 19 gives the PP wave reflectance approximation formula for any weakly isotropic medium; For the rock physics model, the media on both sides of the interface satisfy the weak reflection approximation condition, and the perturbation of the longitudinal wave modulus and the shear modulus of the organic-rich shale is: Therefore, bringing formulas 32, 33, and 34 into formula 19 gives the PP wave linearization reflectance approximation formula as follows: Subsequently, bringing formulas 13 and 17 into formula 35 gives: In order to simplify the formula expression and facilitate subsequent inversion research, formula 38 is simplified as follows: The coefficient term in the formula is defined as: After that, the reflectivity of the variable Mframe is calculated, and the following is obtained: Therefore, formula 41 is simplified as follows: The reflectance formula: In Step E, the content of the seismic rock physics inversion method includes: The relative quantity of the linearization characterized reflectance formula is approximately represented by its corresponding natural logarithm value, that is: In the formula, the symbol "ln(·)" represents the natural logarithm of the elastic parameter, and the symbol "Δ" represents the difference between the values on both sides of the reflection interface; In addition, the elastic properties of the solid skeleton, the density of the organic-rich shale, and the elastic properties related to the organic matter expressed in the logarithmic form are expanded to continuous time variables as follows: Therefore, the reflectance of a single interface represented by formula 37 is expanded to the reflectance equation of continuous time by formulas 45 and 46: In the formula, the symbol "d" represents the difference between the elastic parameters at the t time and the t+1 time on the continuous time; Then, the length of seismic data is set as n+1 and the number of incident angles is h, so there are n reflection interfaces in the time direction of seismic data, and the angle-dependent wavelet matrix W is expressed as: Equation 49: ; The block matrix W in the matrix W i The seismic wavelet associated with the incidence angle θ i is expressed by Equation 50: ; where w i represents the seismic wavelet when the incidence angle is θ i l h; s1, s2, …, s l represent the time sampling points of the seismic wavelet; The seismic forward equation with data error term is obtained by the convolution of the wavelet matrix W and the reflection coefficient sequence, as shown below: Equation 51: ; where S represents the seismic response, W is the seismic wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, e represents the data residual term, which satisfies the mean of 0 and the variance of Σ d Gaussian distribution, and m is the model parameter vector, which is expressed as: Equation 52: ; where the superscript T represents the transpose of the vector, the length of the vector m is q = 5(n+1), because each of the n+1 layers contains five elastic parameters, and finally n reflection points are formed; The dimension of the coefficient matrix Λ is (h×n)×5n, and its expression is: Equation 53: ; where A(θ i ), B(θ i ), C(θ i ), D(θ i ), and E(θ i ) are diagonal matrices whose expressions are: Equation 54: ; Where a(t) j , θ i b(t) j , θ i ), c(t) j , θ i ), d(t) j , θ i ) and e(t j , θ i ) is a coefficient that varies with angle and time, where j = 1, 2, ..., n and i = 1, 2, ..., h; The dimension of the difference matrix D is 5n×5(n+1), and its expression is: Equation 55: ; where the dimension of the block matrix Dd is n×(n+1), and its expression is: Equation 56: ; According to the Bayesian theory, the posterior probability density function P(m|d) is obtained from the prior probability density function P(d|m) and the likelihood function P(m): Equation 57: ; The prior model parameter vector m is set to satisfy a multi-dimensional Gaussian distribution with a mean μ m and a covariance Σ m , and the data residual term satisfies a Gaussian distribution with a mean of 0 and a covariance Σ d , so the posterior probability density function P(m|d) satisfies a Gaussian distribution with a mean μ m|d and a covariance Σ m|d , and the expression is: Equation 58: ; Equation 59: ; where the matrix G = WΛD is the forward operator; Thus, the maximum likelihood solution m opt for the parameters is given by the expression Equation 60: ; where μ m|d and Σ m|d are the mean and covariance of the posterior probability density function P(m|d).

Images