Shale solid skeleton and organic matter elastic parameter synchronous prediction method

Abstract

The invention is applicable to the field of oil and gas geophysical exploration, and provides a shale solid skeleton and organic matter elastic parameter synchronous prediction method, which comprises the following steps of: constructing a rock physical model for representing a shale reservoir solid skeleton elastic modulus and an organic matter elastic modulus; developing an in-well inversion method based on the rock physical model; combining a rock physical model and a scattering function theory, and establishing a reflection coefficient formula represented by shale solid skeleton elasticity modulus, organic matter elasticity modulus and density parameter linearization; a reflection coefficient formula and a Bayesian linear inversion theory are combined to develop an earthquake rock physical inversion method; applying synthetic seismic data, and testing and verifying the effectiveness of the seismic rock physical inversion method constructed in the step E; and then the method is applied to actual seismic data to predict the elasticity modulus of a shale gas reservoir solid skeleton and organic matters in a research area, and an effective method is provided for brittleness evaluation of a shale gas reservoir engineering dessert and organic matter enrichment prediction of a geological dessert.

Description

Technical Field

[0001] The present invention relates to the field of oil and gas geophysical exploration, and in particular to a method for synchronously predicting shale solid skeleton and organic matter elastic parameters. Background Art

[0002] Shale gas, as a significant unconventional oil and gas resource, is attracting increasing attention from both industry and academia. Effective assessment of shale gas reservoirs requires the prediction of both engineering and geological sweet spots. The primary focus of engineering sweet spot prediction is the evaluation of the mechanical properties of the shale's solid skeleton, while the primary focus of geological sweet spot prediction is the identification of shale gas reservoirs with high total organic carbon (TOC) content. Currently, conventional methods for predicting engineering sweet spots in shale gas reservoirs primarily rely on the calculation of mechanical parameters based on seismic elastic parameter inversion. These calculations reflect the overall mechanical properties of the shale. However, due to the influence of factors such as organic matter and fluids in the shale, these methods cannot fully reflect the mechanical properties of the solid skeleton itself. Regarding geological sweet spot prediction, current methods primarily use empirical relationships between TOC and elastic parameters (such as density) derived from conventional seismic inversion to predict shale TOC content. However, these methods face the following challenges: The generalizability of these empirical relationships across different regions is limited; Furthermore, seismic inversion methods for density parameters suffer from inherent issues such as poor stability and low computational accuracy. Therefore, it is necessary to develop a seismic inversion method that directly targets the elastic parameters of the shale solid skeleton and organic matter based on the rock physics research of shale gas reservoirs, so as to provide an effective method for the brittleness evaluation of shale gas reservoir engineering sweet spots and the prediction of organic matter enrichment in geological sweet spots. Summary of the Invention

[0003] In view of the shortcomings of the existing technology, the purpose of the present invention is to propose a method for synchronously predicting the elastic parameters of shale solid skeleton and organic matter, so as to solve the technical problems mentioned in the background technology.

[0004] A method for synchronously predicting shale solid skeleton and organic matter elastic parameters, the method comprising: Step A: Based on the microstructure analysis of shale gas reservoir, construct a solid skeleton elastic modulus to characterize the shale reservoir ( M frame , m frame ) and elastic modulus of organic matter ( M om , m om ) rock physics model; Step B: Develop an in-well inversion method based on the rock physics model constructed in step A to predict the solid skeleton elastic modulus ( Mframe , m frame ) and elastic modulus of organic matter ( M om , m om ), providing an initial model with physical meaning for the inversion of actual pre-stack seismic data; Step C: Combine the rock physics model constructed in step A with the scattering function theory to establish the elastic modulus of the shale solid skeleton ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and the density parameter ( r sat ) Linearized reflection coefficient formula; Step D: conducting an accuracy evaluation of the reflection coefficient formula and a sensitivity analysis of the parameters in the reflection coefficient formula (reflectivity of the solid skeleton elastic modulus, reflectivity of the organic matter elastic modulus, and reflectivity of the organic-rich shale density) to verify the applicability of the reflection coefficient formula; Step E: Combine the reflection coefficient formula established in step C with the Bayesian linear inversion theory to develop a seismic rock physics inversion method to achieve direct and synchronous prediction of the shale solid skeleton elastic modulus ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and other parameters; Step F, applying the solid skeleton elastic modulus and organic matter elastic modulus predicted in step B in combination with actual seismic wavelets to calculate synthetic seismic data, and testing and verifying the effectiveness of the seismic rock physics inversion method constructed in step E; Step G: Apply 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 the shale gas reservoir in the study area, providing a new and effective method for evaluating the brittleness of shale gas reservoir engineering sweet spots and predicting the organic matter enrichment of geological sweet spots.

[0005] Preferably, in step A, the step of constructing a rock physics model characterizing the horizontal fracture density of the shale gas reservoir includes: The Hashin-Shtrikman boundary (HSB) theory is applied to calculate the bulk modulus and shear modulus of solid matrices composed of complex minerals ( K 0, m 0); The bulk modulus and shear modulus of the porous matrix consisting of a solid matrix and water-saturated inorganic pores were calculated using the self-consistent approximation (SCA) theory. K matrix , m matrix ); The bulk modulus and shear modulus of the solid skeleton are calculated using the Krief formula ( K frame , m frame ); The Kuster-Toksöz (KT) formula is used to calculate the bulk modulus and shear modulus of organic matter composed of kerogen and gas-saturated organic pores ( K mix , m mix ); The bulk modulus and shear modulus of organic-rich shale composed of solid skeleton and organic matter ( K sat , m sat ) can be calculated using the Ciz-Shapiro solid substitution theory.

[0006] Preferably, in step B, the well inversion method includes: Two key variables are used in the rock physics inversion process: the volume fraction of organic matter ( f om ) and the aspect ratio of the inorganic pores ( α By fitting the measured P- and S-wave velocities with those calculated by the rock physics model, the objective function is minimized to obtain the optimal volume fraction of organic matter. f om and the aspect ratio of inorganic pores α Numerical, objective function F obj As shown below: Formula 18: ; in, V P and V S Represents the P-wave velocity and S-wave velocity respectively; the subscripts “modeled” and “measured” represent the velocity calculated by the rock physics model and the measured velocity value respectively. In the process of rock physics modeling, the total porosity is set to f Porosity of inorganic pores f inorganic and organic pore porosity f organic Composition, that is f = f organic + f inorganic ; And, the predicted f om and α The numerical values ​​are used as input to the rock physics model to calculate the elastic modulus of the solid skeleton in the well ( M frame and m frame ) and the elastic modulus of organic matter ( M om and m om ), the rock physics inversion results provide an initial model with physical meaning for subsequent seismic rock physics inversion.

[0007] Preferably, in step C, the elastic modulus of the shale solid skeleton ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and the density parameter ( r sat )The contents of the linearized reflection coefficient formula include: Shaw and Sen (2004, 2006) derived the PP wave reflection coefficient on the interface of any weakly 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: Formula 19: ; in i and ψ represent the incident angle and azimuth, respectively. r sat represents the density of saturated rock, and r0 is the point on the horizontal interface that satisfies Snell's law. The scattering function S is related to the perturbation of the elastic parameters, horizontal slowness, and polarization direction, and its expression is: Formula 20: ; in S (r0) is the scattering function, Δ r sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable xThe physical quantity related to the polarization vector of the PP wave, and the variable or ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave; According to Shaw and Sen (2006), the dummy indicator in Equation 20 is ijkl "Simplified to a dummy indicator according to the following formula" IJ ”: Formula 21: ; Formula 22: ; in ijkl represents the subscript value of the fourth-order stiffness tensor, i , j , k , l The values ​​of are 1 to 3; IJ represents the simplified subscript value of the fourth-order stiffness tensor, I , J The values ​​are 1 to 6.

[0008] In the formula d ij and d kl is the Kronecker symbol, so Formula 20 can be simplified as: Formula 23: ; in S (r0) is the scattering function, Δ r sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable x The physical quantity related to the polarization vector of the PP wave, and the variable or ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave Polarization vector of PP wave t and t ' can be expressed by the incident wave and the reflected wave respectively as: Formula 24: ; Formula 25: ; in i and ψ represent the incident angle and azimuth, respectively.

[0009] Its corresponding slowness vector p and p ' is expressed by the polarization vector as: Formula 26: ; Formula 27: ; in α sat is the longitudinal wave velocity of saturated rock; When the PP wave is incident and reflected, the variable x The expression is written as: Formula 28: ; t i and Represents the polarization vector t and t ' i A portion.

[0010] And the variable or IJ The expression is written as: Formula 29: ; Where, Represents the polarization vector t ' i A quantity, Represents the slowness vector p 'jth component, tk represents the kth component of the polarization vector t, p l represents the lth component variable of the slowness vector p; or IJ The expression for each component is: Formula 30: ; or 11 , or 12 , or 13 , or 22 , or 23 , or 33 , or 44 , or 55 , or 66 , or 21 , or 31 , or 32 Represents variables or IJ Components in different directions.

[0011] For weakly isotropic media, the perturbation Δ of the stiffness matrix in Equation 23 is C IJ The representation is simplified to the following form: Formula 31: ; where Δ M sat , Δ l sat and Δ m sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant l sat and the second Lamé constant m sat The amount of disturbance; Substituting Equations 28, 30, and 31 into Equation 23, we can obtain the scattering function S of the weakly isotropic medium, which is expressed as follows: Formula 32: ; where Δ M sat , Δ l sat and Δ m sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant l sat and the second Lamé constant m sat The amount of disturbance; Finally, by substituting Equation 32 into Equation 19, we can obtain the approximate formula for the PP wave reflection coefficient for any weakly isotropic medium; For the rock physics model, assuming that the media on both sides of the interface meet the weak reflection approximation condition, the perturbation of the longitudinal wave modulus and shear modulus of the organic-rich shale is: Formula 33: ; Formula 34: ; Where, Δ M frame and Δ m frame Represent the disturbance of the longitudinal wave modulus and shear modulus of the solid skeleton, Δ M om and Δ m om They represent the disturbance of the corresponding elastic parameters of organic matter, and the disturbance satisfies the following conditions: |Δ M frame |<M frame ,|Δ m frame |< m frame ,|Δ M om |< M om ,|Δ m om |< m om ; Therefore, by substituting Equations 32, 33, and 34 into Equation 19, we obtain the following approximate formula for the linearized reflection coefficient of the PP wave: Formula 35: ; in R PP is the reflection coefficient, i is the angle of incidence, Δ r sat is the disturbance of saturated rock density.

[0012] Then, substituting Equations 13 and 17 into Equation 35, we obtain: Formula 36: ; Afterwards, define c 2 sat = M sat / m sat and c 2 frame = M frame / m frame , respectively represent the square of the ratio of the longitudinal and shear wave velocities of the organic-rich shale and the solid skeleton. Therefore, the above formula can be simplified to Formula 37: ; Finally, define the parameters p = m frame / m sat is the ratio of the shear modulus of the solid skeleton and the organic-rich shale, so Equation 37 can be simplified and the final linearized reflection coefficient formula is obtained as: Formula 38: ; In order to simplify the formula expression and facilitate subsequent inversion research, Formula 38 is simplified to: Formula 39: ; The coefficient term in the formula is defined as: Formula 40: ; The above formula represents the case in which organic matter with non-zero shear modulus is used as solid filling in organic-rich shale; for the case of fluid filling, that is, Figure 2 When the remaining pores in the rock are filled with fluid, the p = 1 condition is met and the Gassmann equation is satisfied. Under this condition, the shear modulus of the solid skeleton is equal to the shear modulus of the rock as a whole, that is, m frame = m sat , so the p = 1 and m frame = m sat Substitute them into formula 38 to obtain the following formula: Formula 41: ; According to the poroelastic medium theory proposed by Russell et al. (2011), the Gassmann equation describes the longitudinal wave modulus of the solid skeleton: M frame and saturated rock shear modulus m sat The relationship is: Formula 42: ; After that, the reflectivity of the variable Mframe is calculated as follows: Formula 43: ; Therefore, Equation 41 simplifies to: Formula 44: ; According to the reflection coefficient formula derived by Russell et al.: Formula 45: ; Equation 44 has the same form as Equation 45, where Δ M om / M om Equivalent to Δ in Formula 45 f / f .

[0013] Preferably, in step E, the contents of the seismic rock physics inversion method include: According to Buland and Omre (2003) and Hampson and Russell (2005), the elastic modulus of the solid skeleton reflectivity Δ M frame / M frame and Δ mframe / m frame , reflectivity of elastic modulus of organic matter Δ M om / M om and Δ m om / m om , density reflectivity of organic-rich shale Δ r sat / r sat It can be approximated by its corresponding natural logarithm, namely: Formula 46: ; Where, 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 organic matter expressed in logarithmic form are expanded into continuous time variables as follows: Formula 47: ; Therefore, the reflection coefficient of a single interface expressed in Equation 37 is expanded into a continuous-time reflection coefficient equation using Equations 45 and 46: Formula 48: ; Where the symbol “d” represents the difference in elastic parameters between the tth moment and the t+1th moment in continuous time; for example, dln[ M frame ( t )] = ln[ M frame ( t +1)] - ln[ M frame ( t )]; Then, the time length of the seismic data is set to n+1, and the number of incident angles is set to h. Therefore, there are n reflection interfaces in the time direction of the seismic data. Therefore, the seismic wavelet matrix W that changes with angle is expressed as: Formula 49: ; The block matrix W in the matrix W i Inclusion and angle of incidence i i The related seismic wavelet is expressed as: Formula 50: ; in w i The incident angle is i i The seismic wavelet at i =1, 2…, h. s 1, s 2,…, s l represents the time sampling point of the seismic wavelet; The seismic wavelet matrix W is convolved with the reflection coefficient sequence to obtain the seismic forward equation with data error terms, as shown below: Formula 51: ; Where S represents the earthquake response, W is the earthquake wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, and e represents the data residual term, which satisfies the mean of 0 and the variance of Σ d Gaussian distribution, m is the model parameter vector, and its expression is: Formula 52: ; The superscript T represents the transpose of the vector, and the length of the vector m is q = 5( n +1), this is because each of the n+1 layers contains five elastic parameters, ultimately forming n reflection points; The dimension of the coefficient matrix Λ is (h×n)×5n, and its expression is: Formula 53: ; in A ( i i ), B ( i i ), C ( i i ), D ( i i )and E ( i i ) is a diagonal matrix, and its expression is: Formula 54: ; in a ( t j , i i ), b ( t j , i i ), c ( t j , ii ), d ( t j , i i )and e ( t j , i 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: Formula 55: ; The dimension of the block matrix Dd is n×(n+1), and the expression is: Formula 56: ; According to 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): Formula 57: ; In the present invention, the prior model parameter vector m is set to satisfy the mean μ m , the covariance is Σ m The multidimensional Gaussian distribution of the data residual term satisfies the mean of 0 and the covariance of Σ d Gaussian distribution, therefore, the posterior probability density function P(m|d) satisfies the Gaussian distribution, and its mean μ m|d and covariance Σ m|d The expression is: Formula 58: ; Formula 59: ; The matrix G = WΛD is the forward operator; Therefore, the maximum likelihood solution m for the inversion parameters opt The expression is: Formula 60: ; where μ m|d and Σ m|d is the mean and covariance of the posterior probability density function P(m|d).

[0014] Beneficial effects achieved by the present invention: This invention enables direct and simultaneous prediction of parameters such as the shale rock skeleton elastic modulus and organic matter elastic modulus from prestack seismic data. The effectiveness of the developed seismic rock physics inversion method was tested and verified using synthetic seismic data. Finally, the developed seismic rock physics inversion method and the constructed initial model were applied to actual prestack seismic data to calculate the elastic moduli of the shale rock skeleton and organic matter in the study area. This provides a new and effective method for evaluating the brittleness of shale gas reservoir engineering sweet spots and predicting organic matter enrichment in geological sweet spots. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a technical flow chart of a method for synchronous prediction of shale solid skeleton and organic matter elastic parameters; Figure 2 It is a flow chart of rock physics modeling constructed based on geological and core data; Figure 3 The longitudinal wave and shear moduli of the organic-rich shale and solid skeleton in Example 1 are as a function of organic matter content. f om and inorganic pore aspect ratio α a and b are the longitudinal wave modulus and shear modulus of organic-rich shale. M sat and m sat With organic matter content f om and inorganic pore aspect ratio α c and d are the longitudinal wave modulus and shear modulus of the solid skeleton of the shale gas reservoir M frame and m frame With organic matter content f om and inorganic pore aspect ratio α Graph of changes in Figure 4 The organic matter-related longitudinal wave modulus and shear modulus in Example 1 are related to the organic matter content. f om and inorganic pore aspect ratio α a and b are the organic matter related longitudinal wave modulus and shear modulus M om and m om With organic matter content f om and inorganic pore aspect ratio α The relationship diagram of the change of organic matter content is shown in Figure 2. f om The horizontal axis is the longitudinal wave modulus and shear modulus related to the organic matter in the shale gas reservoir. M om and m om is the intersection diagram of the ordinate; Figure 5 is a schematic diagram of the seismic reflection model constructed based on geological and logging data in Experimental Example 1; Figure 6 is the reflection coefficient curve of the top and bottom interfaces calculated using different reflection coefficient formulas in Experimental Example 1; Figure 7 is the sensitivity analysis result of the linearization parameters of Experimental Example 1, where a is the contribution of the solid skeleton longitudinal wave modulus reflectivity to the reflection coefficient, b is the contribution of the solid skeleton shear modulus reflectivity to the reflection coefficient, c is the contribution of the organic-rich shale density reflectivity to the reflection coefficient, d is the contribution of the organic-related longitudinal wave modulus to the reflection coefficient, and e is the contribution of the organic-related shear modulus reflectivity to the reflection coefficient; Figure 8 This is the forward and inversion flow chart of the well rock physics in Experimental Example 2; Figure 9 is the actual logging curve of Well A in Experimental Example 2, a is the natural gamma 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 content curve; Figure 10 This is the inversion result of the rock physics in Well A in Experiment 2. a is the P-wave velocity curve, b is the S-wave velocity curve, where the red and black curves represent the rock physics predicted and measured velocities, respectively, and c is the organic matter volume fraction, i.e., the organic matter content. f om , d is the inorganic pore aspect ratio, e is the porosity curve, where the red curve and the black curve represent the organic pore porosity and the total porosity, respectively; Figure 11 These are the target curves for rock physics modeling based on the in-well rock physics inversion results for Well A in Experimental Example 2. a is the P-wave modulus curve, b is the shear modulus curve, and c is the density curve. The red and black curves represent the properties of the shale solid skeleton and organic-rich shale, respectively. d and e represent the organic-matter-related P-wave modulus and shear modulus curves, respectively. Figure 12 Based on the experimental example 3 Figure 10 The constructed synthetic seismic data, a is the solid skeleton longitudinal wave modulus curve, b is the solid skeleton shear modulus curve, c is the organic-rich shale density curve, d is the organic matter-related longitudinal wave modulus curve, e is the organic matter-related shear wave modulus curve, and f is the synthetic seismic data; Figure 13is the inversion result of synthetic seismic data in Experimental Example 3, a is the solid skeleton P-wave modulus curve, b is the solid skeleton shear modulus curve, c is the organic-rich shale density curve, d is the organic matter-related P-wave modulus curve, and e is the organic matter-related S-wave modulus curve. The black curve, red curve, and blue curve represent the actual model, inversion result, and initial model, respectively. f is the synthetic seismic data, where the red curve represents the seismic data forward simulated using the inversion result, and the black line represents the original synthetic seismic data. Figure 14 is the post-stack seismic profile through Well A in Experimental Example 4; Figure 15 is the initial model of Well A in Experimental 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-related P-wave modulus, and e is the organic-related S-wave modulus. Figure 16 is the inversion result of Well A in Experimental Example 4, where 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-related P-wave modulus, and e is the organic-related S-wave modulus. DETAILED DESCRIPTION

[0016] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0017] See also Figure 1 The embodiment of the present invention provides a method for synchronously predicting shale solid skeleton and organic matter elastic parameters, the method comprising: Step A: Based on the microstructure analysis of shale gas reservoir, construct a solid skeleton elastic modulus to characterize the shale reservoir ( M frame , m frame ) and elastic modulus of organic matter ( M om , m om ) rock physics model; Step B: Develop an in-well inversion method based on the rock physics model constructed in step A to predict the solid skeleton elastic modulus ( M frame , m frame ) and elastic modulus of organic matter ( Mom , m om ), providing an initial model with physical meaning for the inversion of actual pre-stack seismic data; Step C: Combine the rock physics model constructed in step A with the scattering function theory to establish the elastic modulus of the shale solid skeleton ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and the density parameter ( r sat ) Linearized reflection coefficient formula; Step D: conducting an accuracy evaluation of the reflection coefficient formula and a sensitivity analysis of the parameters in the reflection coefficient formula (reflectivity of the solid skeleton elastic modulus, reflectivity of the organic matter elastic modulus, and reflectivity of the organic-rich shale density) to verify the applicability of the reflection coefficient formula; Step E: Further combine the reflection coefficient formula established in step C with the Bayesian linear inversion theory to develop a seismic rock physics inversion method to achieve direct and synchronous prediction of the shale solid skeleton elastic modulus ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and other parameters; Step F, applying the solid skeleton elastic modulus and organic matter elastic modulus predicted in step B in combination with actual seismic wavelets to calculate synthetic seismic data, and testing and verifying the effectiveness of the seismic rock physics inversion method constructed in step E; Step G: Apply 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 the shale gas reservoir in the study area, providing a new and effective method for evaluating the brittleness of shale gas reservoir engineering sweet spots and predicting the organic matter enrichment of geological sweet spots.

[0018] See also Figure 2 In step A of this embodiment, the step of constructing a rock physics model to characterize the horizontal fracture density of a shale gas reservoir includes: The Hashin-Shtrikman boundary (HSB) theory is applied to calculate the bulk modulus and shear modulus of solid matrices composed of complex minerals ( K 0, m 0); The bulk modulus and shear modulus of the porous matrix consisting of a solid matrix and water-saturated inorganic pores were calculated using the self-consistent approximation (SCA) theory. K matrix , m matrix ); The bulk modulus and shear modulus of the solid skeleton are calculated using the Krief formula ( K frame , m frame ); The Kuster-Toksöz (KT) formula is used to calculate the bulk modulus and shear modulus of organic matter composed of kerogen and gas-saturated organic pores ( K mix , m mix ); The bulk modulus and shear modulus of organic-rich shale composed of solid skeleton and organic matter ( K sat , m sat ) can be calculated using the Ciz-Shapiro solid substitution theory.

[0019] In step A of this embodiment, the Hashin-Shtrikman bound (HSB) theory includes the following: For a medium containing two different components, its equivalent bulk modulus and shear modulus can be calculated by HSB theory: Formula 1: ; in K HS± and m HS± They are represented as the upper bound + and lower bound – of the equivalent medium bulk modulus; K i and m i Representing the i Elastic modulus of the components; f i Representative i The volume fraction of the components; Makvo et al. (2009) pointed out that for isotropic mixtures containing two-phase media, the theoretical boundary calculated by HSB theory has the narrowest range.

[0020] For equivalent media composed of multiple components, the equivalent elastic modulus can be calculated by the following formula: Formula 2: ; in K HS+ andK HS- represent the upper and lower boundaries of the equivalent medium bulk modulus, m HS+ and m HS- Represent the upper and lower boundaries of the equivalent medium shear modulus, respectively. K max and m max Respectively represent the maximum values ​​of bulk modulus and shear modulus of each component constituting the equivalent medium, K min and m min They represent the minimum bulk modulus and shear modulus of each component that constitutes the equivalent medium.

[0021] Formula 3: ; in, z Represents the independent variable in the function. When calculating the upper limit, z for K max or m max , when calculating the lower limit, z is K min or m min ; Λ( z ),Γ( z )and g ( K , m ) represents the function in the calculation process, K represents bulk modulus and shear modulus; r Represents the rock r The value of each component is related to the rock physics modeling process. The symbol <·> represents the volume weighted average of each component; Bulk modulus of the equivalent medium K 0 and shear modulus m 0 can be calculated from the average of the HSB theoretical upper and lower limits, that is: Formula 4: .

[0022] In step A of this embodiment, the Kuster-Toksöz (KT) formula includes the following: For an equivalent medium containing inclusions of various material properties or shapes within a background medium, Kuster and Toksöz (1974) provided expressions for its bulk modulus and shear modulus. In the present invention, kerogen is assumed to be filled with gas-saturated organic pores to form an organic matter mixture, so the expression for its elastic modulus is: Formula 5: ; in K * mix and m * mix represent the bulk modulus and shear modulus of the organic matter mixture, respectively; K k and m k represent the bulk modulus and shear modulus of kerogen, respectively; K g and m 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 The specific expression can refer to Mavko et al., (2009); x g is the proportion of organic pores in the organic matter mixture, and the variable g There are the following expressions: Formula 6: .

[0023] In step A of this embodiment, the self-consistent approximation (SCA) theory includes the following: When the elastic properties and geometric shapes of the components of the equivalent medium are known, the self-consistent approximation (SCA) theory can accurately calculate the elastic modulus of the equivalent medium. In the present 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 a two-phase mixture, Wu (1966) gave an estimation formula for its elastic modulus: Formula 7: ; in, K matrix and m matrix represent the bulk modulus and shear modulus of the pore matrix, respectively; K 0 and m 0 are the bulk modulus and shear modulus of the solid matrix, respectively; K w and m w are the bulk modulus and shear modulus of water, respectively; x w is the volume fraction of inorganic pores in the porous matrix; the shape of inorganic pores can be determined by the factor P 0w and G 0wThe specific expression can be referred to Mavko et al., (2009). In future studies, for more complex N-phase mixtures, Berryman (1980) gave the general form of its self-consistency approximation: Formula 8: ; in i For the i kind of materials, K i and m i Respectively represent i The bulk modulus and shear modulus of the material, x i is its volume content, P 0w and G 0w Indicates the shape of the inclusion, superscript * i Refers to the inclusion of objects in the calculation process i Geometric factor for the background medium.

[0024] In step A of this embodiment, the Krief formula includes the following contents: Krief gave a relationship describing the elastic modulus of dry rock and solid matrix. In the present invention, the bulk modulus of the solid skeleton is calculated using Krief's formula K frame and shear modulus m frame : Formula 9: ; In the formula K matrix and m matrix represent 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 It is an empirical parameter related to lithofacies, and the value of this invention is 5.

[0025] In step A of this embodiment, the Ciz-Shapiro solid replacement theory includes the following: The organic matter mixture composed of kerogen and its internal organic pores has non-zero bulk and shear moduli (Carcione, 2000; Zhao et al., 2016). Therefore, the Gassmann equation applicable to fluids (shear modulus is zero) is no longer applicable. For the case of pores filled with solid matter, Ciz and Shapiro (2007) proposed the solid replacement theory to calculate the elastic modulus of saturated rocks. For the isotropic case, the bulk modulus of organic-rich shale, i.e., saturated rock, is K sat and shear modulus m sat The expression is: Formula 10: ; Formula 11: ; in f om is the volume fraction of the whole rock composed of organic matter mixture; K mix and m mix are the bulk modulus and shear modulus of the organic mixture; Formulas 10 and 11 are further simplified to the following form: Formula 12: ; Formula 13: ; in, K om and m om are the bulk modulus and shear modulus associated with organic matter; Therefore, the longitudinal wave modulus of organic-rich shale and its solid skeleton is: Formula 14: ; Formula 15: ; in, M sat and M frame are the longitudinal wave moduli of organic-rich shale and solid skeleton, m sat and m frame are the shear moduli of organic-rich shale and solid skeleton respectively; Substituting formulas 12 and 13 into formula 14, we get: Formula 16: ; Further substituting formulas 14 and 15 into formula 16, we can obtain: Formula 17: ; in K om and m om are the bulk modulus and shear modulus associated with organic matter, M om = K om +4 / 3 m om Represents the longitudinal wave modulus associated with organic matter.

[0026] In step B of this embodiment, the well inversion method includes: Two key variables are used in the rock physics inversion process: the volume fraction of organic matter ( f om ) and the aspect ratio of the inorganic pores ( α By fitting the measured P- and S-wave velocities with those calculated by the rock physics model, the objective function is minimized to obtain the optimal volume fraction of organic matter. f om and the aspect ratio of inorganic pores α Numerical, objective function F obj As shown below: Formula 18: ; in, V P and V S Represents the P-wave velocity and S-wave velocity respectively; the subscripts “modeled” and “measured” represent the velocity calculated by the rock physics model and the measured velocity value respectively. In the process of rock physics modeling, the total porosity is set to f Porosity of inorganic pores f inorganic and organic pore porosity f organic Composition, that is f = f organic + f inorganic ; And, the predicted f om and α The numerical values ​​are used as input to the rock physics model to calculate the elastic modulus of the solid skeleton in the well ( M frame and m frame ) and the elastic modulus of organic matter ( M om and m om), the rock physics inversion results provide an initial model with physical meaning for subsequent seismic rock physics inversion.

[0027] In step C of this embodiment, the elastic modulus of the shale solid skeleton ( M frame , m frame ), organic matter elastic modulus ( M om , m om ) and the density parameter ( r sat )The contents of the linearized reflection coefficient formula include: Shaw and Sen (2004, 2006) derived the PP wave reflection coefficient on the interface of any weakly 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: Formula 19: ; in i and ψ represent the incident angle and azimuth, respectively. r sat represents the density of saturated rock, and r0 is the point on the horizontal interface that satisfies Snell's law. The scattering function S is related to the perturbation of the elastic parameters, horizontal slowness, and polarization direction, and its expression is: Formula 20: ; in S (r0) is the scattering function, Δ r sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable x The physical quantity related to the polarization vector of the PP wave, and the variable or ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave; According to Shaw and Sen (2006), the dummy indicator in Equation 20 is ijkl "Simplified to a dummy indicator according to the following formula" IJ ”: Formula 21: ; Formula 22: ; in ijkl represents the subscript value of the fourth-order stiffness tensor, i , j ,k , l The values ​​of are 1 to 3; IJ represents the simplified subscript value of the fourth-order stiffness tensor, I , J The values ​​are 1 to 6.

[0028] In the formula d ij and d kl is the Kronecker symbol, so Formula 20 can be simplified as: Formula 23: ; in S (r0) is the scattering function, Δ r sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable x The physical quantity related to the polarization vector of the PP wave, and the variable or ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave Polarization vector of PP wave t and t ' can be expressed by the incident wave and the reflected wave respectively as: Formula 24: ; Formula 25: ; in i and ψ represent the incident angle and azimuth, respectively.

[0029] Its corresponding slowness vector p and p ' is expressed by the polarization vector as: Formula 26: ; Formula 27: ; in α sat is the longitudinal wave velocity of saturated rock; When the PP wave is incident and reflected, the variable x The expression is written as: Formula 28: ; t i and Represents the polarization vector t and t ' i A portion.

[0030] And the variable or IJ The expression is written as: Formula 29: ; Where, Represents the polarization vector t ' i A quantity, Represents the slowness vector p 'jth component, tk represents the kth component of the polarization vector t, p l represents the lth component variable of the slowness vector p; or IJ The expression for each component is: Formula 30: ; or 11 , or 12 , or 13 , or 22 , or 23 , or 33 , or 44 , or 55 , or 66 , or 21 , or 31 , or 32 Represents variables or IJ Components in different directions.

[0031] For weakly isotropic media, the perturbation Δ of the stiffness matrix in Equation 23 is C IJ The representation is simplified to the following form: Formula 31: ; where Δ M sat , Δ l sat and Δ m sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant l sat and the second Lamé constant m sat The amount of disturbance; Substituting Equations 28, 30, and 31 into Equation 23, we can obtain the scattering function S of the weakly isotropic medium, which is expressed as follows: Formula 32: ; where Δ M sat , Δ l sat and Δ m sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant l sat and the second Lamé constant m sat The amount of disturbance; Finally, by substituting Equation 32 into Equation 19, we can obtain the approximate formula for the PP wave reflection coefficient for any weakly isotropic medium; For the rock physics model, assuming that the media on both sides of the interface meet the weak reflection approximation condition, the perturbation of the longitudinal wave modulus and shear modulus of the organic-rich shale is: Formula 33: ; Formula 34: ; Where, Δ M frame and Δ m frame Represent the disturbance of the longitudinal wave modulus and shear modulus of the solid skeleton, Δ M om and Δ m om They represent the disturbance of the corresponding elastic parameters of organic matter, and the disturbance satisfies the following conditions: |Δ M frame |< M frame ,|Δ m frame |< m frame ,|Δ M om |< M om ,|Δ m om |< m om ; Therefore, by substituting Equations 32, 33, and 34 into Equation 19, we obtain the following approximate formula for the linearized reflection coefficient of the PP wave: Formula 35: ; in R PP is the reflection coefficient, i is the angle of incidence, Δ r sat is the disturbance of saturated rock density.

[0032] Then, substituting Equations 13 and 17 into Equation 35, we obtain: Formula 36: ; Afterwards, define c 2 sat = M sat / m sat and c 2 frame = M frame / m frame , respectively represent the square of the ratio of the longitudinal and shear wave velocities of the organic-rich shale and the solid skeleton. Therefore, the above formula can be simplified to Formula 37: ; Finally, define the parameters p = m frame / m sat is the ratio of the shear modulus of the solid skeleton and the organic-rich shale, so Equation 37 can be further simplified, and the final linearized reflection coefficient formula is obtained as follows: Formula 38: ; In order to simplify the formula expression and facilitate subsequent inversion research, Formula 38 is simplified to: Formula 39: ; The coefficient term in the formula is defined as: Formula 40: ; The above formula represents the case in which organic matter with non-zero shear modulus is used as solid filling in organic-rich shale; for the case of fluid filling, that is, Figure 2 When the remaining pores in the rock are filled with fluid, the p = 1 condition is met and the Gassmann equation is satisfied. Under this condition, the shear modulus of the solid skeleton is equal to the shear modulus of the rock as a whole, that is, m frame = m sat , so the p = 1 and m frame = msat Substitute them into formula 38 to obtain the following formula: Formula 41: ; According to the poroelastic medium theory proposed by Russell et al. (2011), the Gassmann equation describes the longitudinal wave modulus of the solid skeleton: M frame and saturated rock shear modulus m sat The relationship is: Formula 42: ; After that, the reflectivity of the variable Mframe is calculated as follows: Formula 43: ; Therefore, Equation 41 is further simplified to: Formula 44: ; According to the reflection coefficient formula derived by Russell et al.: Formula 45: ; Equation 44 has the same form as Equation 45, where Δ M om / M om Equivalent to Δ in Formula 45 f / f .

[0033] In step E of this embodiment, the contents of the seismic rock physics inversion method include: Solid skeleton elastic modulus reflectivity Δ M frame / M frame and Δ m frame / m frame , reflectivity of elastic modulus of organic matter Δ M om / M om and Δ m om / m om , density reflectivity of organic-rich shale Δ r sat / r sat It can be approximated by its corresponding natural logarithm, namely: Formula 46: ; Where, 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 organic matter expressed in logarithmic form are expanded into continuous time variables as follows: Formula 47: ; Therefore, the reflection coefficient of a single interface expressed in Equation 37 is expanded into a continuous-time reflection coefficient equation using Equations 45 and 46: Formula 48: ; Where the symbol “d” represents the difference in elastic parameters between the tth moment and the t+1th moment in continuous time; for example, dln[ M frame ( t )] = ln[ M frame ( t +1)] - ln[ M frame ( t )]; Then, the time length of the seismic data is set to n+1, and the number of incident angles is set to h. Therefore, there are n reflection interfaces in the time direction of the seismic data. Therefore, the seismic wavelet matrix W that changes with angle is expressed as: Formula 49: ; The block matrix W in the matrix W i Inclusion and angle of incidence i i The related seismic wavelet is expressed as: Formula 50: ; in w i The incident angle is i i The seismic wavelet at i =1, 2…, h. s 1, s 2,…, s l represents the time sampling point of the seismic wavelet; The seismic wavelet matrix W is convolved with the reflection coefficient sequence to obtain the seismic forward equation with data error terms, as shown below: Formula 51: ; Where S represents the earthquake response, W is the earthquake wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, and e represents the data residual term, which satisfies the mean of 0 and the variance of Σ dGaussian distribution, m is the model parameter vector, and its expression is: Formula 52: ; The superscript T represents the transpose of the vector, and the length of the vector m is q = 5( n +1), this is because each of the n+1 layers contains five elastic parameters, ultimately forming n reflection points; The dimension of the coefficient matrix Λ is (h×n)×5n, and its expression is: Formula 53: ; in A ( i i ), B ( i i ), C ( i i ), D ( i i )and E ( i i ) is a diagonal matrix, and its expression is: Formula 54: ; in a ( t j , i i ), b ( t j , i i ), c ( t j , i i ), d ( t j , i i )and e ( t j , i 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: Formula 55: ; The dimension of the block matrix Dd is n×(n+1), and the expression is: Formula 56: ; According to 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): Formula 57: ; In the present invention, the prior model parameter vector m is set to satisfy the mean μ m , the covariance is Σ m The multidimensional Gaussian distribution of the data residual term satisfies the mean of 0 and the covariance of Σ d Gaussian distribution, therefore, the posterior probability density function P(m|d) satisfies the Gaussian distribution, and its mean μ m|d and covariance Σ m|d The expression is: Formula 58: ; Formula 59: ; The matrix G = WΛD is the forward operator; Therefore, the maximum likelihood solution m for the inversion parameters opt The expression is: Formula 60: ; where μ m|d and Σ m|d is the mean and covariance of the posterior probability density function P(m|d).

[0034] Experimental Example 1, rock physics model test: Figure 3 The saturated longitudinal modulus and saturated shear modulus of the simulated organic-rich shale are shown ( M sat , m sat ), and the longitudinal wave modulus and shear modulus of its solid skeleton ( M frame , m frame ) f om and α The simulation results of the changes in f om The range of variation is 0-0.1, α The range is 10 -1.5 – 1. For ease of presentation, f om is the horizontal axis, αTake the common logarithm, lg( α i ) is the vertical axis. Figure 3 a and b show that as the organic matter content increases or the inorganic pore aspect ratio decreases, M sat and m sat The results show a significant nonlinear downward trend, which indicates that higher organic matter content or narrower inorganic pores will reduce the stiffness of organic-rich shales. Figure 3 c and d show M frame and m frame The results show a similar trend. However, at the same organic matter content and inorganic pore aspect ratio, the elastic modulus of the solid skeleton is always lower than that of the organic-rich shale. These results confirm the influence of organic matter on the shear modulus and longitudinal wave modulus.

[0035] Using the same simulation results, the organic matter-related longitudinal wave modulus and shear modulus are calculated using Equations 17 and 13 ( M om , m om ), the result is as follows Figure 4 a and b. As can be seen from the figure, M om and m om All follow f om and α i increases with the increase of f om Therefore, in order to highlight this dependency, Figure 4 c and d are converted to f om is the horizontal axis, elastic modulus M om and m om are the vertical coordinates, and lg( α ) is the color column. The simulation results show that M om and m om and f om It shows a very significant approximate linear positive correlation, proving its potential as an indicator for directly characterizing the organic matter enrichment degree of shale gas reservoirs.

[0036] In order to verify the precision and accuracy of the proposed reflection coefficient formula 38, the following seismic reflection model containing three shale layers was designed based on the geological background and logging data of the study area, as follows: Figure 5 shown.

[0037] Figure 6 The reflection coefficients of the top and bottom interfaces are calculated using the Zoeppritz equation (1919), the Aki-Richards approximation (2002), and the proposed reflection coefficient approximation. Numerical results show that the reflection coefficients calculated using the proposed reflection coefficient approximation agree closely with those calculated using the traditional Zoeppritz equation and the Aki-Richards approximation. Only at angles of incidence greater than approximately 35° does the reflection coefficient calculated using the proposed method deviate slightly from the more accurate Zoeppritz equation.

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

[0039] Figure 7 The contribution of each elastic parameter to the reflection coefficient is shown, and the parameter range is set to -0.3–0.3. Figure 7 a- Figure 7 c shows the longitudinal wave modulus reflectivity of the solid skeleton (Δ M frame / M frame ) and the density reflectivity of organic-rich shale (Δ r sat / r sat ) contributes to the PP wave reflection coefficient in all incident angles. In contrast, the shear modulus reflectivity of the solid skeleton (Δ m frame / m frame ) begins to make a significant contribution to the reflection coefficient when the incident angle is greater than about 10°. Figure 7 d and e show the organic matter related longitudinal wave modulus reflectivity (Δ M om / M om ) and shear modulus reflectivity (Δ m om / m om) to the reflectivity. The results show that the contribution of organic matter-related elastic modulus reflectivity to the reflectivity follows a similar trend to the corresponding solid skeleton elastic parameters, but the contribution of organic matter-related elastic modulus reflectivity to the reflectivity is smaller. This phenomenon is consistent with rock physics theory predictions that organic matter, as an inclusion, has a weaker influence on the overall elastic parameters of shale than the solid skeleton, and thus a relatively weaker impact on the reflectivity.

[0040] Experimental Example 2, Well Rock Physics Inversion: Figure 8 The forward and inverse flow chart of wellbore rock physics is shown.

[0041] Figure 9 The well logging curves of Well A in the study area are shown, including natural gamma (GR), compressional wave velocity ( V P ), shear wave velocity ( V S ),density( r ), total porosity ( f ), and mineral volume fraction curves. The shale interval ranges from approximately 2670–2825 m, and the two black lines represent the depth range of the target layer.

[0042] Afterwards, the elastic moduli of the solid skeleton and organic matter were inverted from the logging data of Well A using the well inversion method. Figure 10 The results of model-driven rock physics inversion are shown. Figure 10 a and b, the P- and S-wave velocities calculated by rock physics (red line) are highly consistent with the measured P- and S-wave velocities (black line), 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 parameters respectively. f om and α On this basis, using the formula f organic = f om – V kerogen Calculation of shale organic pore porosity ( f organic ),like Figure 10 As shown in e.

[0043] according to Figure 8 The rock physics forward and inversion process, the above inversion results f om and α will be used as input to the rock physics model to calculate the elastic modulus from the well log data ( M frame , m frame , M om , m om ), the calculation results are as follows Figure 11 shown. Figure 11 a and b show the longitudinal wave modulus and shear modulus of organic-rich shale and solid skeleton, represented by black and red lines, respectively. Figure 11 c shows the density of organic-rich shale (black line) and the corresponding density of the solid skeleton (red line). Finally, the organic matter-related longitudinal wave modulus and shear modulus are calculated, as shown in Figure 11 As shown in d and e.

[0044] Experimental Example 3, synthetic data test: based on Figure 12 The rock physics forward and inversion results are first converted to time-depth, converting the depth domain logging curve to the time domain. Then the time domain logging curve is resampled at 2ms intervals to match the time sampling interval of the actual seismic data. The real model in the time domain is generated as follows Figure 12 a- Figure 12 As shown in e. Figure 12 f shows the noise-free synthetic seismic record generated using Equation 50 and the actual seismic wavelet.

[0045] Subsequently, the seismic rock physics inversion method proposed in this invention will be applied to Figure 12 The synthetic seismic data of f. The inversion results are as follows Figure 13 a Figure 13 e, where the black solid line represents the true model, the blue realization represents the initial model, and the red dashed line represents the inversion result. Figure 13 f shows the earthquake record (red) obtained by using the inversion results to perform earthquake forward modeling. Figure 12 Comparison results of the noise-free synthetic seismic record in f. The high agreement between the inversion results and the real model and the consistency of the corresponding seismic responses verify the effectiveness and robustness of the inversion method and results.

[0046] Experimental Example 4, actual data application: Figure 14 is the post-stack seismic data of Well A.

[0047] Figure 15 In order to use Kriging interpolation method to calculate the actual well logging data ( Figure 11 ) to perform interpolation and extrapolation, and perform low-pass filtering on the interpolation results to obtain the actual initial model.

[0048] based on Figure 15 The elastic parameters in the actual work area were inverted using the initial model, post-stack seismic data and actual seismic wavelets.M frame , m frame , r sat , M om and m om ,like Figure 16 According to the petrophysical meaning of the proposed petrophysical model, Figure 16 a and b show the elastic modulus of the solid skeleton ( M frame , m 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 the shale, especially the relative content of clay, quartz and carbonate minerals in the shale of the study area. This shows that through appropriate calibration of comprehensive logging data, M frame and m frame It can be used as a quantitative indicator to evaluate the brittleness index, which is crucial for optimizing the hydraulic fracturing effect of shale gas reservoirs. At the same time, through the analysis of rock physics and well logging data mentioned above, the elastic modulus related to organic matter ( M om , m om ) can serve as a key indicator for quantitatively characterizing organic matter occurrence. The inverted parameters can provide an effective method for evaluating the brittleness of engineering sweet spots in shale gas reservoirs and predicting organic matter enrichment in geological sweet spots.

[0049] The above are only preferred embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the description of the present invention, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A method for synchronously predicting shale solid skeleton and organic matter elastic parameters, characterized in that: The method comprises: Step A: Based on the microstructural analysis of the shale gas reservoir, a rock physics model is constructed to characterize the elastic modulus of the solid skeleton and the elastic modulus of organic matter in the shale reservoir; Step B: Developing an in-well inversion method based on the rock physics model constructed in Step A, predicting the elastic modulus of the solid skeleton and the elastic modulus of organic matter from conventional well logging data, and providing an initial model with physical meaning for inversion of actual pre-stack seismic data; Step C: combining the rock physics model constructed in step A with the scattering function theory to establish a reflection coefficient formula characterized by the linearization of the shale solid skeleton elastic modulus, the organic matter elastic modulus, and the density parameter; Step D: Conducting an accuracy evaluation of the reflection coefficient formula and a sensitivity analysis of the parameters in the reflection coefficient formula to verify the applicability of the reflection coefficient formula; Step E: Combining the reflection coefficient formula established in step C with the Bayesian linear inversion theory, a seismic rock physics inversion method is developed to achieve direct and synchronous prediction of the shale solid skeleton elastic modulus and organic matter elastic modulus from pre-stack seismic data; Step F, applying the solid skeleton elastic modulus and organic matter elastic modulus predicted in step B in combination with actual seismic wavelets to calculate synthetic seismic data, and testing and verifying the effectiveness of the seismic rock physics inversion method constructed in step E; Step G: Apply 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 the shale gas reservoir in the study area, providing a new and effective method for evaluating the brittleness of shale gas reservoir engineering sweet spots and predicting the organic matter enrichment of geological sweet spots.

2. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step A, the steps of constructing a rock physics model to characterize the horizontal fracture density of a shale gas reservoir include: The Hashin-Shtrikman boundary theory is applied to calculate the bulk modulus and shear modulus of solid matrices composed of complex minerals; The bulk modulus and shear modulus of the porous matrix consisting of a solid matrix and water-saturated inorganic pores are calculated using the self-compatibility approximation theory. The bulk modulus and shear modulus of the solid skeleton are calculated using Krief's formula; The Kuster-Toksöz (KT) formula is used to calculate the bulk modulus and shear modulus of organic matter composed of kerogen and gas-saturated organic pores; The overall bulk modulus and shear modulus of organic-rich shale, which is composed of a solid skeleton and organic matter, can be calculated using the Ciz-Shapiro solid substitution theory.

3. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 2, characterized in that: In step A, the Hashin-Shtrikman boundary theory includes the following: For a medium containing two different components, the equivalent bulk modulus and shear modulus are calculated using the HSB theory: Formula 1: ; in K HS± and μ HS± They are represented as the upper bound + and lower bound – of the equivalent medium bulk modulus; K i and μ i Representing the i Elastic modulus of the components; f i Representative i Volume fraction of each component; For equivalent media composed of multiple components, the equivalent elastic modulus is calculated by the following formula: Formula 2: ; in K HS+ and K HS- represent the upper and lower boundaries of the equivalent medium bulk modulus, μ HS+ and μ HS- represent the upper and lower boundaries of the equivalent medium shear modulus, K max and μ max Respectively represent the maximum values ​​of bulk modulus and shear modulus of each component constituting the equivalent medium, K min and μ min Respectively represent the minimum values ​​of bulk modulus and shear modulus of each component constituting the equivalent medium; Formula 3: ; in, z Represents 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 the function in the calculation process, K represents bulk modulus and shear modulus; r Represents the rock r The value of each component is related to the rock physics modeling process; the symbol <·> represents the volume weighted average of each component; Bulk modulus of the equivalent medium K 0 and shear modulus μ 0 is calculated from the average of the HSB theoretical upper and lower limits, that is: Formula 4: .

4. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step A, the Kuster-Toksöz formula includes the following: Assuming that the kerogen is filled with gas-saturated organic pores to form an organic matter mixture, the expression of its elastic modulus is: Formula 5: ; in K * mix and μ * mix represent the bulk modulus and shear modulus of the organic matter mixture, respectively; K k and μ k 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 express; x g is the proportion of organic pores in the organic matter mixture, and the variable ζ k There are the following expressions: Formula 6: .

5. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step A, the self-consistent approximate theory includes the following: The bulk modulus and shear modulus of the porous matrix are obtained by adding water-saturated inorganic pores to the solid matrix. For a two-phase mixture, the elastic modulus is estimated by: Formula 7: ; in, K matrix and μ matrix represent the bulk modulus and shear modulus of the pore matrix, respectively; K 0 and μ 0 are the bulk modulus and shear modulus of the solid matrix, respectively; K w and μ w are the bulk modulus and shear modulus of water, respectively; x w is the volume fraction of inorganic pores in the pore matrix; the shape of inorganic pores is determined by the structure factor P 0w and G 0w It means that for more complex N-phase mixtures, the self-consistency approximation takes the form: Formula 8: ; in i For the i kind of materials, K i and μ i Respectively represent i The bulk modulus and shear modulus of the material, x i is its volume content, P 0w and G 0w Indicates the shape of the inclusion, superscript * i Refers to the inclusion of objects in the calculation process i Geometric factor for the background medium.

6. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step A, the Krief formula includes the following: Calculating the bulk modulus of a solid skeleton using the Krief formula K frame and shear modulus μ frame : Formula 9: ; In the formula K matrix and μ matrix represent 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 It is an empirical parameter related to lithofacies.

7. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step A, the Ciz-Shapiro solid substitution theory includes the following: For the isotropic case, the bulk modulus of organic-rich shale, i.e., saturated rock, is K sat and shear modulus μ sat The expression is: Formula 10: ; Formula 11: ; in f om is the volume fraction of the whole rock composed of organic matter mixture; K mix and μ mix are the bulk modulus and shear modulus of the organic mixture; Formulas 10 and 11 are simplified to the following form: Formula 12: ; Formula 13: ; in, K om and μ om are the bulk modulus and shear modulus associated with organic matter; Therefore, the longitudinal wave modulus of organic-rich shale and its solid skeleton is: Formula 14: ; Formula 15: ; in, M sat and M frame are the longitudinal wave moduli of organic-rich shale and solid skeleton, μ sat and μ frame are the shear moduli of organic-rich shale and solid skeleton respectively; Substituting formulas 12 and 13 into formula 14, we get: Formula 16: ; Substituting formulas 14 and 15 into formula 16, we get: Formula 17: ; in K om and μ om are the bulk modulus and shear modulus associated with organic matter, M om = K om +4 / 3 μ om Represents the longitudinal wave modulus associated with organic matter.

8. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step B, the well inversion method includes: By fitting the measured P- and S-wave velocities with those calculated by the rock physics model, the objective function is minimized to obtain the optimal volume fraction of organic matter. f om and the aspect ratio of inorganic pores α Numerical, objective function F obj As shown below: Formula 18: ; in, V P and V S Represents the P-wave velocity and S-wave velocity respectively; the subscripts "modeled" and "measured" represent the velocity calculated by the rock physics model and the measured velocity value respectively. In the process of rock physics modeling, the total porosity is set φ Porosity of inorganic pores φ inorganic and organic pore porosity φ organic Composition, that is φ = φ organic + φ inorganic ; And, the predicted f om and α The numerical values ​​are used as inputs to the rock physics model to calculate the elastic modulus of the solid skeleton in the well and the elastic modulus of the organic matter. The rock physics inversion results provide an initial model with physical meaning for subsequent seismic rock physics inversion.

9. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step C, the contents of establishing the reflection coefficient formula represented by the linearization of the shale solid skeleton elastic modulus, organic matter elastic modulus, and density parameters include: PP wave reflection coefficient on the interface of arbitrary weakly anisotropic media R PP The relationship with the scattering function S is: Formula 19: ; in θ and ψ represent the incident angle and azimuth, respectively. ρ sat represents the density of 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, horizontal slowness, and polarization direction, and its expression is: Formula 20: ; in S (r0) is the scattering function, Δ ρ sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable ξ The physical quantity related to the polarization vector of the PP wave, and the variable η ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave; Dummy indicator in formula 20" ijkl "Simplified to a dummy indicator according to the following formula" IJ ”: Formula 21: ; Formula 22: ; in ijkl represents the subscript value of the fourth-order stiffness tensor, i , j , k , l The values ​​of are 1 to 3; IJ represents the simplified subscript value of the fourth-order stiffness tensor, I , J The values ​​of are 1 to 6; In the formula δ ij and δ kl is the Kronecker symbol, so Formula 20 is simplified to: Formula 23: ; in S (r0) is the scattering function, Δ ρ sat is the disturbance of saturated rock density, Δ C ijkl is the perturbation of the saturated rock stiffness tensor, and the variable ξ The physical quantity related to the polarization vector of the PP wave, and the variable η ijkl It is a physical quantity related to the polarization vector and slowness vector of the PP wave; Polarization vector of PP wave t and t ' can be expressed by the incident wave and the reflected wave respectively as: Formula 24: ; Formula 25: ; in θ and ψ represent the incident angle and azimuth respectively; Its corresponding slowness vector p and p ' is expressed by the polarization vector as: Formula 26: ; Formula 27: ; in α sat is the longitudinal wave velocity of saturated rock; When the PP wave is incident and reflected, the variable ξ The expression is written as: Formula 28: ; t i and Represents the polarization vector t and t ' i Quantity And the variable η IJ The expression is written as: Formula 29: ; Where, Represents the polarization vector t ' i A quantity, Represents the slowness vector p 'jth component, tk represents the kth component of the polarization vector t, p l represents the lth component variable of the slowness vector p; η IJ The expression for each component is: Formula 30: ; η 11 , η 12 , η 13 , η 22 , η 23 , η 33 , η 44 , η 55 , η 66 , η 21 , η 31 , η 32 Represents variables η IJ Components in different directions; For weakly isotropic media, the perturbation Δ of the stiffness matrix in Equation 23 is C IJ The representation is simplified to the following form: Formula 31: ; where Δ M sat , Δ λ sat and Δ μ sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant λ sat and the second Lamé constant μ sat The amount of disturbance; Substituting Equations 28, 30, and 31 into Equation 23, we can obtain the scattering function of the weakly isotropic medium, which is expressed as follows: Formula 32: ; where Δ M sat , Δ λ sat and Δ μ sat are the longitudinal wave modulus of saturated rock M sat , the first Lamé constant λ sat and the second Lamé constant μ sat The amount of disturbance; Finally, by substituting Equation 32 into Equation 19, we can obtain the approximate formula for the PP wave reflection coefficient for any weakly isotropic medium; For the rock physics model, assuming that the media on both sides of the interface meet the weak reflection approximation condition, the perturbation of the longitudinal wave modulus and shear modulus of the organic-rich shale is: Formula 33: ; Formula 34: ; Where, Δ M frame and Δ μ frame Represent the disturbance of the longitudinal wave modulus and shear modulus of the solid skeleton, Δ M om and Δ μ om They represent the disturbance of the corresponding elastic parameters of organic matter, and the disturbance satisfies the following conditions: |Δ M frame |< M frame ,|Δ μ frame |< μ frame ,|Δ M om |< M om ,|Δ μ om |< μ om ; Therefore, by substituting Equations 32, 33, and 34 into Equation 19, we obtain the following approximate formula for the linearized reflection coefficient of the PP wave: Formula 35: ; in R PP is the reflection coefficient, θ is the angle of incidence, Δ ρ sat is the disturbance of saturated rock density; Then, substituting Equations 13 and 17 into Equation 35, we obtain: Formula 36: ; Afterwards, define γ 2 sat = M sat / μ sat and γ 2 frame = M frame / μ frame , respectively represent the square of the ratio of the longitudinal and shear wave velocities of the organic-rich shale and the solid skeleton. Therefore, the above formula is simplified to Formula 37: ; Finally, define the parameters p = μ frame / μ sat is the ratio of the shear modulus of the solid skeleton and the organic-rich shale, so Equation 37 can be simplified and the final linearized reflection coefficient formula is obtained as: Formula 38: ; In order to simplify the formula expression and facilitate subsequent inversion research, Formula 38 is simplified to: Formula 39: ; The coefficient term in the formula is defined as: Formula 40: ; Will p = 1 and μ frame = μ sat Substitute them into formula 38 to obtain the following formula: Formula 41: ; Gassmann equation describes the longitudinal wave modulus of solid skeleton M frame and saturated rock shear modulus μ sat The relationship is: Formula 42: ; After that, the reflectivity of the variable Mframe is calculated as follows: Formula 43: ; Therefore, Equation 41 simplifies to: Formula 44: ; Reflection coefficient formula: Formula 45: ; Equation 44 has the same form as Equation 45, where Δ M om / M om Equivalent to Δ in Formula 45 f / f .

10. The method for synchronously predicting shale solid skeleton and organic matter elastic parameters according to claim 1, characterized in that: In step E, the contents of the seismic rock physics inversion method include: The relative amount of the linearized reflection coefficient formula is approximated by its corresponding natural logarithm, namely: Formula 46: ; Where, 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 organic matter expressed in logarithmic form are expanded into continuous time variables as follows: Formula 47: ; Therefore, the reflection coefficient of a single interface expressed in Equation 37 is expanded into a continuous-time reflection coefficient equation using Equations 45 and 46: Formula 48: ; Wherein, the symbol "d" represents the difference in elastic parameters between the tth moment and the t+1th moment in continuous time; Then, the time length of the seismic data is set to n+1, and the number of incident angles is set to h. Therefore, there are n reflection interfaces in the time direction of the seismic data. Therefore, the seismic wavelet matrix W that changes with angle is expressed as: Formula 49: ; The block matrix W in the matrix W i Inclusion and angle of incidence θ i The related seismic wavelet is expressed as: Formula 50: ; in w i The incident angle is θ i The seismic wavelet at i =1, 2…, h; s 1, s 2,…, s l represents the time sampling point of the seismic wavelet; The seismic wavelet matrix W is convolved with the reflection coefficient sequence to obtain the seismic forward equation with data error terms, as shown below: Formula 51: ; Where S represents the earthquake response, W is the earthquake wavelet matrix, Λ is the coefficient matrix, D is the difference matrix, and e represents the data residual term, which satisfies the mean of 0 and the variance of Σ d Gaussian distribution, m is the model parameter vector, and its expression is: Formula 52: ; The superscript T represents the transpose of the vector, and the length of the vector m is q = 5( n +1), this is because each of the n+1 layers contains five elastic parameters, ultimately forming n reflection points; The dimension of the coefficient matrix Λ is (h×n)×5n, and its expression is: Formula 53: ; in A ( θ i ), B ( θ i ), C ( θ i ), D ( θ i )and E ( θ i ) is a diagonal matrix, and its expression is: Formula 54: ; in 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: Formula 55: ; The dimension of the block matrix Dd is n×(n+1), and the expression is: Formula 56: ; According to 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): Formula 57: ; Set the prior model parameter vector m to satisfy the mean μ m , the covariance is Σ m The multidimensional Gaussian distribution of the data residual term satisfies the mean of 0 and the covariance of Σ d Gaussian distribution, therefore, the posterior probability density function P(m|d) satisfies the Gaussian distribution, and its mean μ m|d and covariance Σ m|d The expression is: Formula 58: ; Formula 59: ; The matrix G = WΛD is the forward operator; Therefore, the maximum likelihood solution m for the inversion parameters opt The expression is: Formula 60: ; where μ m|d and Σ m|d is the mean and covariance of the posterior probability density function P(m|d).

Images