Simultaneous direct inversion method, system and equipment based on fluid sensitive factor and brittleness index, and storage medium

Abstract

The invention discloses a simultaneous direct inversion method, system and equipment based on a fluid sensitive factor and a brittleness index, and a storage medium, and the method comprises the steps: obtaining an expression of a stratum reflection coefficient, obtaining an approximate expression of a first seismic record S (theta) according to the expression of the stratum reflection coefficient, simplifying the approximate expression of the first seismic record S (theta), obtaining a forward equation, and obtaining a forward inversion result of the first seismic record S (theta); and based on the forward equation, obtaining inversion results of the fluid sensitive factor F, the brittleness index BI and the density rho through an inversion algorithm. According to the method, high-precision parameters capable of reflecting reservoir fluid properties and brittleness characteristics can be obtained on the basis of a simultaneous direct inversion method of the fluid sensitive factors and the brittleness indexes, fluid identification and brittleness evaluation results of a unified scale can be provided, and a reference can be provided for reservoir fluid identification and reservoir transformation potential; and support is provided for depicting oil and gas reservoir space distribution and developing dessert areas.

Description

Technical Field

[0001] The invention belongs to the technical field of oil and gas field development engineering, and in particular relates to a simultaneous direct inversion method, system, equipment and storage medium based on a fluid sensitivity factor and a brittleness index. Background Art

[0002] Shale gas has become an important research area in the exploration and development of unconventional oil and gas resources. Shale is widely distributed in the world, and most shale gas reservoirs have low permeability and low porosity. It is necessary to create reservoir cracks through fracturing to increase the production capacity of natural gas. Finding sweet spots in shale gas reservoirs is the key factor for its successful exploration.

[0003] In recent years, seismic inversion using prestack seismic data to invert fluid sensitivity factors for evaluating gas content in shale gas target reservoirs and brittleness indices for quantifying reservoir brittleness characteristics has been developing rapidly. Currently, most studies have only conducted a single evaluation of reservoir fluid properties or rock brittleness, and are usually indirect inversion methods, that is, first inverting elastic parameters and then calculating the required fluid sensitivity factors and brittleness indices through rock physics equations. Due to the uncertainty of the inherent "ill-conditioned solutions" of seismic inverse problems, indirect methods will inevitably cause cumulative errors, resulting in strong multi-solutions in the inversion results, which is difficult to meet the needs of actual shale reservoir fluid identification and rock brittleness characteristic evaluation. Summary of the invention

[0004] The purpose of the present invention is to provide a simultaneous direct inversion method, system, equipment and storage medium based on fluid sensitivity factor and brittleness index, so as to solve the problem that the existing single parameter is difficult to meet the demand for accurate evaluation of fluid-containing shale reservoir characteristics, and the indirect inversion method will cause cumulative errors, resulting in strong multi-solution problems in the inversion results.

[0005] The embodiment of the present application is implemented in this way, based on the simultaneous direct inversion method of the fluid sensitivity factor and the brittleness index, including: obtaining an expression of the formation reflection coefficient, as follows:

[0006]

[0007] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0008] According to the expression of formation reflection coefficient, the approximate formula of the first seismic record S(θ) is as follows:

[0009] S(θ)=K F (θ)WDL F +K BI (θ)WDL BI +K ρ (θ)WDL ρ

[0010]

[0011] Where S(θ) is the first seismic record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0012] By simplifying the approximate formula of the first seismic record S(θ), the forward equation is as follows:

[0013] d=Gr

[0014] Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0015] Based on the forward equation, the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ are obtained through the inversion algorithm.

[0016] Optionally, in some embodiments of the present application, an expression for the relationship between the reflection coefficient and the longitudinal wave velocity, the shear wave velocity and the density is provided as follows:

[0017]

[0018] Where R(θ) is the formation reflection coefficient; θ is the incident angle; V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the strata on both sides of the interface; Sis the shear wave velocity; ΔV S is the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is the density; and / or

[0019] The fluid sensitivity factor F is expressed as:

[0020]

[0021] Where F is the fluid sensitivity factor; λ is the first Lame constant; ρ is the density; E is Young's modulus; σ is Poisson's ratio; and / or

[0022] The brittleness index is expressed as:

[0023]

[0024] Where BI is the brittleness index; E is Young's modulus; λ is the first Lame constant; σ is Poisson's ratio.

[0025] Optionally, in some embodiments of the present application, in an isotropic elastic medium, the expressions for the longitudinal wave velocity and the shear wave velocity are respectively as follows:

[0026]

[0027] Where V P is the longitudinal wave velocity; V S is the shear wave velocity; ρ is the density; M is the longitudinal wave modulus; μ is the shear wave modulus;

[0028] According to the chain derivation rule of multivariate calculus, the expressions of longitudinal wave velocity and transverse wave velocity become:

[0029]

[0030] Where V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the strata on both sides of the interface; S is the shear wave velocity; ΔV S is the change in shear wave velocity of the formation on both sides of the interface; Δρ is the change in density of the formation on both sides of the interface; ρ is density; M is the P-wave modulus; ΔM is the change in P-wave modulus of the formation on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the formation on both sides of the interface; and / or

[0031] In isotropic elastic media, the expressions for the longitudinal wave modulus and the transverse wave modulus are as follows:

[0032]

[0033] Where M is the longitudinal wave modulus; μ is the transverse wave modulus; E is Young's modulus; σ is Poisson's ratio.

[0034] Optionally, in some embodiments of the present application, an intermediate variable H is introduced A , H B and H C , H A , H B and H C The expressions are as follows:

[0035]

[0036] Where σ is Poisson's ratio;

[0037] get:

[0038] Where, F is the fluid sensitivity factor; BI is the brittleness index; E is Young's modulus; M is the longitudinal wave modulus; μ is the transverse wave modulus, and ρ is the density;

[0039] because:

[0040]

[0041] Where σ is Poisson's ratio, ΔH A H is the intermediate variable A Change; ΔH B H is the intermediate variable B Change; ΔH C H is the intermediate variable C Change; H A , H B and H C All are intermediate variables;

[0042] get:

[0043]

[0044] Where F is the fluid sensitivity factor; ΔF is the change in the formation fluid sensitivity factor on both sides of the interface; BI is the brittleness index; ΔBI is the change in the brittleness index of the formation on both sides of the interface; E is Young's modulus; ΔE is the change in Young's modulus of the formation on both sides of the interface; M is the P-wave modulus; ΔM is the change in P-wave modulus of the formation on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the formation on both sides of the interface; ρ is the density; Δρ is the change in density of the formation on both sides of the interface; σ is the Poisson's ratio; Δσ is the change in Poisson's ratio;

[0045] because:

[0046]

[0047] Where σ is Poisson’s ratio, and

[0048] Obtain an expression for the formation reflectivity.

[0049] Optionally, in some embodiments of the present application, the method for obtaining the approximate formula of the first seismic record S(θ) is as follows:

[0050] The first seismic record is formed by convolution of the seismic wavelet and the formation reflection coefficient, which is expressed as:

[0051] S=WR+N

[0052] In the formula, S is seismic data; R is seismic reflection coefficient sequence; N is noise; W is wavelet matrix;

[0053] Likewise, the reflection coefficient can be expressed in terms of impedance:

[0054]

[0055] In the formula, R i is the reflection coefficient of the i-th layer and the i+1-th layer; Z i is the wave impedance value of the i-th layer; Z i+1 is the wave impedance value of the i+1th layer; D is the difference matrix, D is

[0056]

[0057] Combined with the expression of the formation reflection coefficient, an approximate expression of the first seismic record S(θ) is obtained.

[0058] Optionally, in some embodiments of the present application, the method for obtaining the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ by the inversion algorithm is as follows:

[0059] According to the forward equation, the first objective function of the inversion is obtained:

[0060]

[0061] Where, J(r) is the first objective function; is the L2 norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0062] The Lp norm is used to construct the reflection sparse constraint term, and the initial model constraints of fluid sensitivity factor F, brittle index BI and density ρ are added to obtain the second objective function:

[0063]

[0064] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter;

[0065] The Lagrangian term R and the dual term C are used to transform the second objective function into the third objective function:

[0066]

[0067] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρis the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of fluid sensitivity factor F; μ2 is the initial model constraint parameter of brittleness index BI; μ3 is the initial model constraint parameter of density ρ; λ is the sparse term constraint parameter; R is the Lagrangian term; C is the dual term; η is the weight coefficient of the dual term.

[0068] Optionally, in some embodiments of the present application, the third objective function is decomposed into single-parameter sub-problems related to r, R and C, thereby obtaining solutions of different sub-functions:

[0069]

[0070] R i+1 =max(|Dr i+1 +C i |-(η / λ) p-2 ·|Dr i+1 +C i | p-1 ,0)·sign(Dr i+1 +C i )

[0071] C i+1 =C i +Dr i+1 -R i+1

[0072] Where E is the unit matrix; are the initial model constraint parameters of elastic parameters; is the logarithmic operation result of the initial model of the parameters to be inverted; r i+1 =[L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; R i+1 is the Lagrangian term; C i+1 is the dual term; i+1 is the number of iterations; d is the seismic record; G is the forward operator; D represents the difference matrix; η is the weight coefficient of the dual term; λ is the sparse term constraint parameter; max is the maximum value; p is the sparse term parameter, ranging from 0 to 1; the sign function is:

[0073]

[0074] By coordinating the solutions of the sub-functions about the elastic parameter r to be inverted, the sub-functions about the Lagrangian term R, and the sub-functions about the dual term C, the optimal solution r of the third objective function can be obtained. The exponential operation of r is performed to obtain the final inversion result:

[0075] [F,BI,ρ] T =exp(r)

[0076] Where r is the optimal solution of the inversion objective function, F is the final inversion result of the fluid sensitivity factor, BI is the final inversion result of the brittleness index, and ρ is the final density inversion result.

[0077] Accordingly, the embodiment of the present application also provides a simultaneous direct inversion system based on a fluid sensitivity factor and a brittleness index, including:

[0078] The formation reflection coefficient construction module is used to obtain the expression of the formation reflection coefficient, as follows:

[0079]

[0080] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0081] The first seismic record S(θ) construction module is used to obtain the approximate formula of the first seismic record S(θ) according to the expression of the formation reflection coefficient as follows:

[0082] S(θ)=K F (θ)WDL F +K BI (θ)WDL BI +K ρ (θ)WDL ρ

[0083]

[0084] Where S(θ) is the first seismic record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0085] The forward equation construction module is used to simplify the approximate formula of the first seismic record S(θ) and obtain the forward equation as follows:

[0086] d=Gr

[0087] Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0088] The inversion algorithm module obtains the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ through the inversion algorithm based on the forward equation.

[0089] Accordingly, an embodiment of the present application also provides a computer device, including a storage device and a processor, wherein the storage device stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the above method.

[0090] Accordingly, an embodiment of the present application also provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the processor executes the steps of the above method.

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

[0092] In order to simultaneously obtain high-precision inversion results of fluid sensitivity factor and brittleness index, the present application derives an expression of formation reflection coefficient related to fluid sensitivity factor, brittleness index and density. According to the expression of formation reflection coefficient, an approximate expression of the first seismic record S(θ) is obtained. The approximate expression of the first seismic record S(θ) is simplified to obtain a forward equation. Based on the forward equation, the inversion results of fluid sensitivity factor, brittleness index and density are obtained through an inversion algorithm, thereby achieving simultaneous direct inversion of fluid sensitivity factor and brittleness index, avoiding the calculation error of indirectly obtaining a single fluid sensitivity factor or a single brittleness index by traditional methods, and improving the accuracy of the inversion results of fluid sensitivity factor and brittleness index. The present application can obtain high-precision parameters that can reflect the properties and brittleness characteristics of reservoir fluids based on the simultaneous direct inversion method of fluid sensitivity factor and brittleness index, can provide unified scale fluid identification and brittleness evaluation results, can provide a reference for the identification of reservoir fluids and the potential for reservoir transformation, and provide support for the characterization of the spatial distribution of oil and gas reservoirs and the development of sweet spots. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] Figure 1 A comparison diagram of the positive wave impedance interface reflection coefficient provided by the present invention;

[0094] Figure 2 A comparison diagram of the negative wave impedance interface reflection coefficient provided by the present invention;

[0095] Figure 3 The multi-angle pre-stack seismic synthetic recording images under different noise conditions provided by the present invention;

[0096] Figure 4 The inversion result diagram of fluid sensitivity factor, brittleness index and density under noise-free condition provided by the present invention;

[0097] Figure 5 The inversion result diagram of fluid sensitivity factor, brittleness index and density provided by the present invention when the signal-to-noise ratio is 5:1

[0098] Figure 6 The inversion result diagram of fluid sensitivity factor, brittleness index and density provided by the present invention when the signal-to-noise ratio is 2:1. DETAILED DESCRIPTION

[0099] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

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

[0101] In a first aspect, the present application provides a simultaneous direct inversion method based on a fluid sensitivity factor and a brittleness index, comprising:

[0102] S01. The expression of formation reflection coefficient is obtained as follows:

[0103]

[0104] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0105] S02. According to the expression of the formation reflection coefficient, the approximate formula of the first seismic record S(θ) is obtained as follows:

[0106] S(θ)=K F (θ)WDL F+K BI (θ)WDL BI +K ρ (θ)WDL ρ

[0107]

[0108] Where S(θ) is the first seismic record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0109] S03. Simplify the approximate formula of the first seismic record S(θ) and obtain the forward equation as follows:

[0110] d=Gr

[0111] Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0112] S04. Based on the forward equation, an inversion algorithm is used to obtain inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ.

[0113] In order to simultaneously obtain high-precision inversion results of fluid sensitivity factor and brittleness index in this application, an expression of formation reflection coefficient related to fluid sensitivity factor, brittleness index and density is derived. According to the expression of formation reflection coefficient, an approximate formula of the first seismic record S(θ) is obtained. The approximate formula of the first seismic record S(θ) is simplified to obtain a forward equation. Based on the forward equation, the inversion results of fluid sensitivity factor, brittleness index and density are obtained through an inversion algorithm, realizing the simultaneous direct inversion of fluid sensitivity factor and brittleness index, avoiding the calculation error of the traditional method of indirectly obtaining a single fluid sensitivity factor or a single brittleness index, and improving the accuracy of the inversion results of fluid sensitivity factor and brittleness index. This application can obtain high-precision parameters that can reflect the properties and brittleness characteristics of reservoir fluids based on the simultaneous direct inversion method of fluid sensitivity factor and brittleness index, can provide unified scale fluid identification and brittleness evaluation results, can provide a reference for the identification of reservoir fluids and the potential for reservoir transformation, and provide support for the characterization of the spatial distribution of oil and gas reservoirs and the development of sweet spots.

[0114] In the S01:

[0115] In some embodiments, the relationship between the reflection coefficient and the longitudinal wave velocity, the shear wave velocity and the density is provided as follows:

[0116]

[0117] Where R(θ) is the formation reflection coefficient; θ is the incident angle; V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the formation on both sides of the interface; S is the shear wave velocity; ΔV S is the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is the density;

[0118] In some embodiments, the fluid sensitivity factor F is expressed as:

[0119]

[0120] Where F is the fluid sensitivity factor; λ is the first Lame constant; ρ is the density; E is the Young's modulus; σ is the Poisson's ratio.

[0121] It can be understood that the fluid sensitivity factor can sensitively characterize rock characteristics and is a commonly used elastic parameter in fluid prediction.

[0122] In some embodiments, the brittleness index is expressed as:

[0123]

[0124] Where BI is the brittleness index; E is Young's modulus; λ is the first Lame constant; σ is Poisson's ratio.

[0125] It can be understood that the brittleness index is E / λ, which is the ratio of Young's modulus to the first Lame constant, and it can better reflect the characteristics of brittle shale.

[0126] In some embodiments, in an isotropic elastic medium, the expressions for the longitudinal wave velocity and the shear wave velocity are respectively as follows:

[0127]

[0128] Where V P is the longitudinal wave velocity; V S is the shear wave velocity; ρ is the density; M is the longitudinal wave modulus; μ is the shear wave modulus.

[0129] It can be understood that an isotropic elastic medium is a medium in which the elastic properties and the propagation speed of seismic waves in the medium do not change depending on the propagation direction.

[0130] Furthermore, according to the chain derivation rule of multivariate calculus, the expressions of longitudinal wave velocity and transverse wave velocity become:

[0131]

[0132] Where V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the formation on both sides of the interface; S is the shear wave velocity; ΔV S is the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is density; M is the P-wave modulus; ΔM is the change in P-wave modulus of the strata on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the strata on both sides of the interface.

[0133] In some embodiments, in an isotropic elastic medium, the expressions for the longitudinal wave modulus and the transverse wave modulus are respectively as follows:

[0134]

[0135] Where M is the longitudinal wave modulus; μ is the transverse wave modulus; E is Young's modulus; σ is Poisson's ratio.

[0136] In some embodiments, an intermediate variable H is introduced A , H B and H C , H A , H B and H C The expressions are as follows:

[0137]

[0138] Where σ is Poisson's ratio.

[0139] Further, we obtain:

[0140] Wherein, F is the fluid sensitivity factor; BI is the brittleness index; E is Young's modulus; M is the longitudinal wave modulus; μ is the transverse wave modulus, and ρ is the density.

[0141] Understandably, The formula is Bring in and get.

[0142] Furthermore, due to:

[0143]

[0144] Where σ is Poisson's ratio, ΔH A H is the intermediate variable A Change; ΔH B H is the intermediate variable B Change; ΔH C H is the intermediate variable C Change; H A , H B and H C All are intermediate variables.

[0145] Understandably, There are the following expressions:

[0146]

[0147] Where, σ is Poisson's ratio; σ1 is Poisson's ratio of the upper stratum; σ2 is Poisson's ratio of the lower stratum; H A is the intermediate variable; ΔH A H is the intermediate variable A Amount of change;

[0148] At the same time, there are the following relationships:

[0149]

[0150] Wherein, σ is Poisson's ratio; Δσ is the change of Poisson's ratio; σ1 is the Poisson's ratio of the upper stratum; σ2 is the Poisson's ratio of the lower stratum;

[0151] When the change of elastic parameters across the elastic interface is relatively weak, the high-order terms can be approximately ignored and the following results can be obtained:

[0152]

[0153] Similarly, we can get:

[0154]

[0155] Where σ is Poisson's ratio; Δσ is the change in Poisson's ratio; ΔH B H is the intermediate variable B Change; ΔH C H is the intermediate variable C Change; H B and H C All are intermediate variables.

[0156] Further, we obtain:

[0157]

[0158] Where F is the fluid sensitivity factor; ΔF is the change in the formation fluid sensitivity factor on both sides of the interface; BI is the brittleness index; ΔBI is the change in the brittleness index of the formation on both sides of the interface; E is Young's modulus; ΔE is the change in Young's modulus of the formation on both sides of the interface; M is the longitudinal wave modulus; ΔM is the change in the longitudinal wave modulus of the formation on both sides of the interface; μ is the shear wave modulus; Δμ is the change in the shear wave modulus of the formation on both sides of the interface; ρ is the density; Δρ is the change in the density of the formation on both sides of the interface; σ is the Poisson's ratio; Δσ is the change in the Poisson's ratio.

[0159] Furthermore, because:

[0160]

[0161] Where σ is Poisson’s ratio, and

[0162] The expression for obtaining the formation reflection coefficient is as follows:

[0163]

[0164]

[0165] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity.

[0166] It can be understood that the expression of the formation reflection coefficient is

[0167] and Obtained by calculation.

[0168] In the S02:

[0169] In some embodiments, the approximate formula of the first seismic record S(θ) is obtained as follows:

[0170] The first seismic record is formed by convolution of the seismic wavelet and the formation reflection coefficient, which is expressed as:

[0171] S=WR+N

[0172] In the formula, S is seismic data; R is seismic reflection coefficient sequence; N is noise; W is wavelet matrix;

[0173] Likewise, the reflection coefficient can be expressed in terms of impedance:

[0174]

[0175] In the formula, R i is the reflection coefficient of the i-th layer and the i+1-th layer; Z i is the wave impedance value of the i-th layer; Z i+1 is the wave impedance value of the i+1th layer; D is the difference matrix, D is

[0176]

[0177] Combined with the expression of the formation reflection coefficient, an approximate expression of the first seismic record S(θ) is obtained.

[0178] In said S03:

[0179] It can be understood that the first seismic record S(θ) is the same as the second seismic record d, and the first and second are added only to distinguish them.

[0180] In said S04:

[0181] In some embodiments, the method of obtaining the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ by the inversion algorithm is as follows:

[0182] According to the forward equation, the first objective function of the inversion is obtained:

[0183]

[0184] Where, J(r) is the first objective function; is the L2 norm; d is the second earthquake record; G is the forward operator; r = [LF ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0185] The Lp norm is used to construct the reflection sparse constraint term, and the initial model constraints of fluid sensitivity factor F, brittle index BI and density ρ are added to obtain the second objective function:

[0186]

[0187] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter;

[0188] The Lagrangian term R and the dual term C are used to transform the second objective function into the third objective function:

[0189]

[0190] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F=lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter; R is the Lagrangian term; C is the dual term; η is the weight coefficient of the dual term;

[0191] The third objective function is decomposed into single-parameter sub-problems related to r, R and C, and then the solutions of different sub-functions are obtained:

[0192]

[0193] R i+1 =max(|Dr i+1 +C i |-(η / λ) p-2 ·|Dr i+1 +C i | p-1 ,0)·sign(Dr i+1 +C i )

[0194] C i+1 =C i +Dr i+1 -R i+1

[0195] Where E is the unit matrix; are the initial model constraint parameters of elastic parameters; is the logarithmic operation result of the initial model of the parameters to be inverted; r i+1 =[L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; R i+1 is the Lagrangian term; C i+1 is the dual term; i+1 is the number of iterations; d is the seismic record; G is the forward operator; D represents the difference matrix; η is the weight coefficient of the dual term; λ is the sparse term constraint parameter; max is the maximum value; p is the sparse term parameter, ranging from 0 to 1; the sign function is:

[0196]

[0197] By coordinating the solutions of the sub-functions about the elastic parameter r to be inverted, the sub-functions about the Lagrangian term R, and the sub-functions about the dual term C, the optimal solution r of the third objective function can be obtained. The exponential operation of r is performed to obtain the final inversion result:

[0198] [F,BI,ρ] T =exp(r)

[0199] Where r is the optimal solution of the inversion objective function, F is the final inversion result of the fluid sensitivity factor, BI is the final inversion result of the brittleness index, and ρ is the final density inversion result.

[0200] It can be understood that the initial model constraints of fluid sensitivity factor, brittleness index and density are added. Since the reflection coefficient of the elastic parameter is also sparse, the Lp norm is used to construct the reflection sparse constraint term to obtain the second objective function.

[0201] It can also be understood that the second objective function includes L p and L2 norm, and the alternating direction multiplier algorithm is used to solve it. The Lagrangian term R is introduced to transform the second objective function into a constrained optimization problem. Then the dual term C is introduced to transform the second objective function into an unconstrained optimization problem, and the third objective function is obtained.

[0202] It can also be understood that the third objective function containing multiple parameters is decomposed into single-parameter sub-problems related to r, R and C, and then solutions of different sub-functions are obtained.

[0203] In a second aspect, the embodiment of the present application provides a simultaneous direct inversion system based on a fluid sensitivity factor and a brittleness index, including: a formation reflection coefficient construction module, which is used to obtain an expression for the formation reflection coefficient, as follows:

[0204]

[0205] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0206] The first seismic record S(θ) construction module is used to obtain the approximate formula of the first seismic record S(θ) according to the expression of the formation reflection coefficient as follows:

[0207] S(θ)=K F (θ)WDL F +K BI (θ)WDL BI +K ρ (θ)WDL ρ

[0208]

[0209] Where S(θ) is the first seismic record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity;

[0210] The forward equation construction module is used to simplify the approximate formula of the first seismic record S(θ) and obtain the forward equation as follows:

[0211] d=Gr

[0212] Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0213] The inversion algorithm module obtains the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ through the inversion algorithm based on the forward equation.

[0214] In the formation reflection coefficient construction module:

[0215] In some embodiments, the relationship between the reflection coefficient and the longitudinal wave velocity, the shear wave velocity and the density is provided as follows:

[0216]

[0217] Where R(θ) is the formation reflection coefficient; θ is the incident angle; V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the strata on both sides of the interface; Sis the shear wave velocity; ΔV S is the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is the density;

[0218] In some embodiments, the fluid sensitivity factor F is expressed as:

[0219]

[0220] Where F is the fluid sensitivity factor; λ is the first Lame constant; ρ is the density; E is the Young's modulus; σ is the Poisson's ratio.

[0221] It can be understood that the fluid sensitivity factor can sensitively characterize rock characteristics and is a commonly used elastic parameter in fluid prediction.

[0222] In some embodiments, the brittleness index is expressed as:

[0223]

[0224] Where BI is the brittleness index; E is Young's modulus; λ is the first Lame constant; σ is Poisson's ratio.

[0225] It can be understood that the brittleness index is E / λ, which is the ratio of Young's modulus to the first Lame constant, and it can better reflect the characteristics of brittle shale.

[0226] In some embodiments, in an isotropic elastic medium, the expressions for the longitudinal wave velocity and the shear wave velocity are respectively as follows:

[0227]

[0228] Where V P is the longitudinal wave velocity; V S is the shear wave velocity; ρ is the density; M is the longitudinal wave modulus; μ is the shear wave modulus.

[0229] It can be understood that an isotropic elastic medium is a medium in which the elastic properties and the propagation speed of seismic waves in the medium do not change depending on the propagation direction.

[0230] Furthermore, according to the chain derivation rule of multivariate calculus, the expressions of longitudinal wave velocity and transverse wave velocity become:

[0231]

[0232] Where V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the formation on both sides of the interface; S is the shear wave velocity; ΔV Sis the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is density; M is the P-wave modulus; ΔM is the change in P-wave modulus of the strata on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the strata on both sides of the interface.

[0233] In some embodiments, in an isotropic elastic medium, the expressions for the longitudinal wave modulus and the transverse wave modulus are respectively as follows:

[0234]

[0235] Where M is the longitudinal wave modulus; μ is the transverse wave modulus; E is Young's modulus; σ is Poisson's ratio.

[0236] In some embodiments, an intermediate variable H is introduced A , H B and H C , H A , H B and H C The expressions are as follows:

[0237]

[0238] Where σ is Poisson's ratio.

[0239] Further, we obtain:

[0240] Wherein, F is the fluid sensitivity factor; BI is the brittleness index; E is Young's modulus; M is the longitudinal wave modulus; μ is the transverse wave modulus, and ρ is the density.

[0241] Understandably, The formula is Bring in and get.

[0242] Furthermore, due to:

[0243]

[0244] Where σ is Poisson's ratio, ΔH A H is the intermediate variable A Change; ΔH B H is the intermediate variable B Change; ΔH C H is the intermediate variable C Change; H A , H B and H C All are intermediate variables.

[0245] Understandably, There are the following expressions:

[0246]

[0247] Where, σ is Poisson's ratio; σ1 is Poisson's ratio of the upper stratum; σ2 is Poisson's ratio of the lower stratum; H A is the intermediate variable; ΔH A H is the intermediate variable A Amount of change;

[0248] At the same time, there are the following relationships:

[0249]

[0250] Wherein, σ is Poisson's ratio; Δσ is the change of Poisson's ratio; σ1 is the Poisson's ratio of the upper stratum; σ2 is the Poisson's ratio of the lower stratum;

[0251] When the change of elastic parameters across the elastic interface is relatively weak, the high-order terms can be approximately ignored and the following results can be obtained:

[0252]

[0253] Similarly, we can get:

[0254]

[0255] Where σ is Poisson's ratio; Δσ is the change in Poisson's ratio; ΔH B H is the intermediate variable B Change; ΔH C H is the intermediate variable C Change; H B and H C All are intermediate variables.

[0256] Further, we obtain:

[0257]

[0258] Where F is the fluid sensitivity factor; ΔF is the change in the formation fluid sensitivity factor on both sides of the interface; BI is the brittleness index; ΔBI is the change in the brittleness index of the formation on both sides of the interface; E is Young's modulus; ΔE is the change in Young's modulus of the formation on both sides of the interface; M is the longitudinal wave modulus; ΔM is the change in the longitudinal wave modulus of the formation on both sides of the interface; μ is the shear wave modulus; Δμ is the change in the shear wave modulus of the formation on both sides of the interface; ρ is the density; Δρ is the change in the density of the formation on both sides of the interface; σ is the Poisson's ratio; Δσ is the change in the Poisson's ratio.

[0259] Furthermore, because:

[0260]

[0261] Where σ is Poisson’s ratio, and

[0262] The expression for obtaining the formation reflection coefficient is as follows:

[0263]

[0264] Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in stratum density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity.

[0265] It can be understood that the expression of the formation reflection coefficient is

[0266] In the first earthquake record S(θ) construction module:

[0267] In some embodiments, the approximate formula of the first seismic record S(θ) is obtained as follows:

[0268] The first seismic record is formed by convolution of the seismic wavelet and the formation reflection coefficient, which is expressed as:

[0269] S=WR+N

[0270] In the formula, S is seismic data; R is seismic reflection coefficient sequence; N is noise; W is wavelet matrix;

[0271] Likewise, the reflection coefficient can be expressed in terms of impedance:

[0272]

[0273] In the formula, R i is the reflection coefficient of the i-th layer and the i+1-th layer; Z i is the wave impedance value of the i-th layer; Z i+1 is the wave impedance value of the i+1th layer; D is the difference matrix, D is

[0274]

[0275] Combined with the expression of the formation reflection coefficient, an approximate expression of the first seismic record S(θ) is obtained.

[0276] In the forward equation construction module:

[0277] It can be understood that the first seismic record S(θ) is the same as the second seismic record d, and the first and second are added only to distinguish them.

[0278] In the inversion algorithm module:

[0279] In some embodiments, the method of obtaining the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ by the inversion algorithm is as follows:

[0280] According to the forward equation, the first objective function of the inversion is obtained:

[0281]

[0282] Where, J(r) is the first objective function; is the L2 norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix;

[0283] The Lp norm is used to construct the reflection sparse constraint term, and the initial model constraints of fluid sensitivity factor F, brittle index BI and density ρ are added to obtain the second objective function:

[0284]

[0285] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρis the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter;

[0286] The Lagrangian term R and the dual term C are used to transform the second objective function into the third objective function:

[0287]

[0288] Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter; R is the Lagrangian term; C is the dual term; η is the weight coefficient of the dual term;

[0289] The third objective function is decomposed into single-parameter sub-problems related to r, R and C, and then the solutions of different sub-functions are obtained:

[0290]

[0291] R i+1 =max(|Dr i+1 +C i |-(η / λ) p-2 ·|Dr i+1 +C i | p-1 ,0)·sign(Dr i+1 +C i )

[0292] C i+1 =C i +Dr i+1 -Ri+1

[0293] Where E is the unit matrix; are the initial model constraint parameters of elastic parameters; is the logarithmic operation result of the initial model of the parameters to be inverted; r i+1 =[L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; R i+1 is the Lagrangian term; C i+1 is the dual term; i+1 is the number of iterations; d is the seismic record; G is the forward operator; D represents the difference matrix; η is the weight coefficient of the dual term; λ is the sparse term constraint parameter; max is the maximum value; p is the sparse term parameter, ranging from 0 to 1; the sign function is:

[0294]

[0295] By coordinating the solutions of the sub-functions about the elastic parameter r to be inverted, the sub-functions about the Lagrangian term R, and the sub-functions about the dual term C, the optimal solution r of the third objective function can be obtained. The exponential operation of r is performed to obtain the final inversion result:

[0296] [F,BI,ρ] T =exp(r)

[0297] Where r is the optimal solution of the inversion objective function, F is the final inversion result of the fluid sensitivity factor, BI is the final inversion result of the brittleness index, and ρ is the final density inversion result.

[0298] It can be understood that the initial model constraints of fluid sensitivity factor, brittleness index and density are added. Since the reflection coefficient of the elastic parameter is also sparse, the Lp norm is used to construct the reflection sparse constraint term to obtain the second objective function.

[0299] It can also be understood that the second objective function includes L p and L2 norm, and the alternating direction multiplier algorithm is used to solve it. The Lagrangian term R is introduced to transform the second objective function into a constrained optimization problem. Then the dual term C is introduced to transform the second objective function into an unconstrained optimization problem, and the third objective function is obtained.

[0300] It can also be understood that the third objective function containing multiple parameters is decomposed into single-parameter sub-problems related to r, R and C, and then solutions of different sub-functions are obtained.

[0301] In a third aspect, the present application provides a computer device comprising a storage device and a processor, wherein the storage device stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the simultaneous direct inversion method based on fluid sensitivity factor and brittleness index as described above.

[0302] The computer device may be a desktop computer, a notebook, a PDA, a cloud server, etc. The computer device may interact with the user through a keyboard, a mouse, a remote control, a touch pad, or a voice control device.

[0303] The memory includes at least one type of readable storage medium, and the readable storage medium includes flash memory, hard disk, multimedia card, card-type memory (for example, SD or D interface display memory, etc.), random access memory (RAM), static random access memory (SRAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), programmable read-only memory (PROM), magnetic memory, disk, optical disk, etc. In some embodiments, the memory can be an internal storage unit of the computer device, such as a hard disk or memory of the computer device. In other embodiments, the memory can also be an external storage device of the computer device, such as a plug-in hard disk equipped on the computer device, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. Of course, the memory can also include both the internal storage unit of the computer device and its external storage device. In this embodiment, the memory is often used to store the operating system and various application software installed on the computer device, such as the program code of the simultaneous direct inversion method based on fluid sensitivity factor and brittleness index, etc. In addition, the memory can also be used to temporarily store various types of data that have been output or are to be output.

[0304] The processor may be a central processing unit (CPU), a controller, a microcontroller, a microprocessor, or other data processing chip in some embodiments. The processor is generally used to control the overall operation of the computer device. In this embodiment, the processor is used to run the program code stored in the memory or process data, such as running the program code of the simultaneous direct inversion method based on the fluid sensitivity factor and the brittleness index.

[0305] In a fourth aspect, the present application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, enables the processor to perform the steps of the simultaneous direct inversion method based on fluid sensitivity factor and brittleness index as described above.

[0306] The computer-readable storage medium stores an interface display program, and the interface display program can be executed by at least one processor to enable the at least one processor to perform the steps of the simultaneous direct inversion method based on fluid sensitivity factor and brittleness index as described above.

[0307] Test Example 1

[0308] Approximate accuracy analysis:

[0309] The three-layer gas-bearing sandstone and shale model given by Goodway et al. based on measured data was used for testing and analysis. The relevant elastic parameters are shown in Table 1.

[0310] Table 1 Related elastic parameters

[0311]

[0312] Based on the above model, the reflection coefficients at different interfaces within the incident angle range of 0 to 40° were calculated using the precise Zoeppritz equation, the Aki-Richard approximate expression and the expression of the formation reflection coefficient proposed in this application, as well as the difference between the reflection coefficient calculated by the expression of the formation reflection coefficient and the reflection coefficient calculated by the precise Zoeppritz equation.

[0313] Figure 1 For the comparison of the positive wave impedance interface reflection coefficient, Figure 1 (a) is the comparison of reflection coefficients under different equations. Figure 1 (b) is the comparison of reflection coefficient differences under different equations. Figure 2 is the comparison of the reflection coefficient of the negative wave impedance interface, Figure 2 (a) is the comparison of reflection coefficients under different equations. Figure 1 (b) Figure 2 (b) is the comparison of reflection coefficient differences under different equations. Figure 1 (a) Figure 2 (a) shows the reflection coefficient response curves of the three expressions at the reflection interface of overlying sandstone and underlying shale and the reflection interface of overlying shale and underlying gas-bearing sandstone, respectively. Figure 1(b) is the difference between the Aki-Richard approximate expression and the expression of the formation reflection coefficient proposed in this application and the reflection coefficient calculated by the precise Zoeppritz equation. As can be seen from the figure, the expression of the formation reflection coefficient proposed in this application is basically consistent with the precise Zoeppritz equation within the range of 40° of the incident angle, and its error is also within an acceptable range, which can be used for AVO elastic parameter seismic inversion.

[0314] Test Example 2

[0315] Model testing:

[0316] In order to verify the reliability and practicality of the direct inversion method combining fluid sensitivity factor and brittleness index in this application, three sets of synthetic seismic data tests were carried out based on actual logging data.

[0317] The reflectivity of the three-layer gas-bearing sandstone and shale model was convolved with a Richer wavelet with a main frequency of 30 Hz. Figure 3 The multi-angle pre-stack seismic data are obtained by adding different degrees of random Gaussian noise to the seismic records. Figure 3 (b) and Figure 3 (c) Noisy seismic records with signal-to-noise ratios (SNRs) of 5 and 2, respectively, are used to test the stability and robustness of the inversion method. Figure 3 (a) is a noise-free seismic record.

[0318] Figure 4 The figure shows the fluid sensitivity factor F, brittleness index BI and density results inverted from the noise-free seismic angle gather. Figure 4 It can be seen from the inversion results that in the noise-free case, all elastic parameters can be well inverted and are consistent with the original well curve trajectory. Figure 5 and Figure 6 These are the inversion results when the SNR (seismic signal noise ratio) is 5 and 2 respectively. Figure 5 and Figure 6 It can be seen that the accuracy of the inversion results is slightly reduced due to the addition of noise, but even when the SNR is 2, the fluid sensitivity factor F and the brittleness index BI can be well inverted, showing good consistency with the logging curve. Overall, the inversion results are acceptable.

[0319] Figures 1 to 3 In, Reflection coefficients: reflection coefficient; Incident angles: incident angle; ER: error; Zoepprizt: exact Zoeppritz equation; Aki: Aki-Richard approximate expression; Proposed: approximate expression proposed in this application; Time (ms): time (milliseconds).

[0320] Figures 4 to 6 In Chinese, Well: well; Initial model: initial model; Inversion: inversion.

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

Claims

1. A simultaneous direct inversion method based on fluid sensitivity factor and brittleness index, characterized in that: include: The expression of formation reflection coefficient is obtained as follows: Where, R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the variation of the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the brittleness index of the formation on both sides of the interface; Δρ is the change in formation density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity; According to the expression of formation reflection coefficient, the approximate formula of the first seismic record S(θ) is as follows: S(θ)=K F (θ)WDL F +K BI (θ)WDL BI +K ρ (θ)WDL ρ Where, S(θ) is the first earthquake record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity; By simplifying the approximate formula of the first seismic record S(θ), the forward equation is as follows: d=Gr Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; Based on the forward equation, the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ are obtained through the inversion algorithm.

2. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 1 is characterized in that: The relationship between the reflection coefficient and the longitudinal wave velocity, shear wave velocity and density is provided as follows: Where R(θ) is the formation reflection coefficient; θ is the incident angle; V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the strata on both sides of the interface; S is the shear wave velocity; ΔV S is the change in shear wave velocity of the strata on both sides of the interface; Δρ is the change in density of the strata on both sides of the interface; ρ is the density; and / or The fluid sensitivity factor F is expressed as: Where F is the fluid sensitivity factor; λ is the first Lame constant; ρ is the density; E is Young's modulus; σ is Poisson's ratio; and / or The brittleness index is expressed as: Where BI is the brittleness index; E is Young's modulus; λ is the first Lame constant; σ is Poisson's ratio.

3. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 2 is characterized in that: In isotropic elastic media, the expressions for the longitudinal wave velocity and the shear wave velocity are as follows: Where V P is the longitudinal wave velocity; V S is the shear wave velocity; ρ is the density; M is the longitudinal wave modulus; μ is the shear wave modulus; According to the chain derivation rule of multivariate calculus, the expressions of longitudinal wave velocity and transverse wave velocity become: Where V P is the longitudinal wave velocity; ΔV P V is the change in the P-wave velocity of the strata on both sides of the interface; S is the shear wave velocity; ΔV S is the change in shear wave velocity of the formation on both sides of the interface; Δρ is the change in density of the formation on both sides of the interface; ρ is density; M is the P-wave modulus; ΔM is the change in P-wave modulus of the formation on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the formation on both sides of the interface; and / or In isotropic elastic media, the expressions for the longitudinal wave modulus and the transverse wave modulus are as follows: Where M is the longitudinal wave modulus; μ is the transverse wave modulus; E is Young's modulus; σ is Poisson's ratio.

4. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 3, characterized in that: Introducing the intermediate variable H A , H B and H C , H A , H B and H C The expressions are as follows: OR B =1+σ, Where σ is Poisson's ratio; get: BI=H A OR B , Where, F is the fluid sensitivity factor; BI is the brittleness index; E is Young's modulus; M is the longitudinal wave modulus; μ is the transverse wave modulus, and ρ is the density; because: Where σ is Poisson's ratio, ΔΗA H is the intermediate variable A Amount of change; ΔΗB H is the intermediate variable B Amount of change; ΔΗC H is the intermediate variable C Change; H A , H B and H C All are intermediate variables; get: Where F is the fluid sensitivity factor; ΔF is the change in the formation fluid sensitivity factor on both sides of the interface; BI is the brittleness index; ΔBI is the change in the brittleness index of the formation on both sides of the interface; E is Young's modulus; ΔE is the change in Young's modulus of the formation on both sides of the interface; M is the P-wave modulus; ΔM is the change in P-wave modulus of the formation on both sides of the interface; μ is the S-wave modulus; Δμ is the change in S-wave modulus of the formation on both sides of the interface; ρ is the density; Δρ is the change in density of the formation on both sides of the interface; σ is the Poisson's ratio; Δσ is the change in Poisson's ratio; because: Where σ is Poisson’s ratio, and Obtain an expression for the formation reflectivity.

5. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 1, characterized in that: The approximate formula of the first seismic record S(θ) is obtained as follows: The first seismic record is formed by convolution of the seismic wavelet and the formation reflection coefficient, which is expressed as: S=WR+N In the formula, S is seismic data; R is seismic reflection coefficient sequence; N is noise; W is wavelet matrix; Likewise, the reflection coefficient can be expressed in terms of impedance: In the formula, R i is the reflection coefficient of the i-th layer and the i+1-th layer; Z i is the wave impedance value of the i-th layer; Z i+1 is the wave impedance value of the i+1th layer; D is the difference matrix, D is Combined with the expression of the formation reflection coefficient, an approximate expression of the first seismic record S(θ) is obtained.

6. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 1, characterized in that: The method of obtaining the inversion results of fluid sensitivity factor F, brittleness index BI and density ρ through the inversion algorithm is as follows: According to the forward equation, the first objective function of the inversion is obtained: Where, J(r) is the first objective function; is the L2 norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; The Lp norm is used to construct the reflection sparse constraint term, and the initial model constraints of fluid sensitivity factor F, brittle index BI and density ρ are added to obtain the second objective function: Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of the fluid sensitivity factor F; μ2 is the initial model constraint parameter of the brittleness index BI; μ3 is the initial model constraint parameter of the density ρ; λ is the sparse term constraint parameter; The Lagrangian term R and the dual term C are used to transform the second objective function into the third objective function: Where J(r) is the objective function; is the L2 norm; is the Lp norm; d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; D is the difference matrix; l F is the logarithmic operation result of the initial model of the fluid sensitivity factor; l BI is the logarithmic operation result of the initial model of the brittleness index; l ρ is the logarithmic operation result of the initial model of density; μ1 is the initial model constraint parameter of fluid sensitivity factor F; μ2 is the initial model constraint parameter of brittleness index BI; μ3 is the initial model constraint parameter of density ρ; λ is the sparse term constraint parameter; R is the Lagrangian term; C is the dual term; η is the weight coefficient of the dual term.

7. The simultaneous direct inversion method based on fluid sensitivity factor and brittleness index according to claim 6, characterized in that: The third objective function is decomposed into single-parameter sub-problems related to r, R and C, and then the solutions of different sub-functions are obtained: R i+1 =max(|D r i+1 +C i |-(η / λ) p-2 ·|D r i+1 +C i | p-1 ,0)·sign(D r i+1 +C i ) C i+1 =C i +D r i+1 -R i+1 Where E is the unit matrix; are the initial model constraint parameters of elastic parameters; is the logarithmic operation result of the initial model of the parameters to be inverted; r i+1 =[L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; R i+1 is the Lagrangian term; C i+1 is the dual term; i+1 is the number of iterations; d is the seismic record; G is the forward operator; D represents the difference matrix; η is the weight coefficient of the dual term; λ is the sparse term constraint parameter; max is the maximum value; p is the sparse term parameter, ranging from 0 to 1; the sign function is: By coordinating the solutions of the sub-functions about the elastic parameter r to be inverted, the sub-functions about the Lagrangian term R, and the sub-functions about the dual term C, the optimal solution r of the third objective function can be obtained. The exponential operation of r is performed to obtain the final inversion result: [F,BI,ρ] T =exp(r) Where r is the optimal solution of the inversion objective function, F is the final inversion result of the fluid sensitivity factor, BI is the final inversion result of the brittleness index, and ρ is the final density inversion result.

8. A simultaneous direct inversion system based on fluid sensitivity factor and brittleness index, characterized in that: include: The formation reflection coefficient construction module is used to obtain the expression of the formation reflection coefficient as follows: Where R(θ) is the formation reflection coefficient; θ is the incident angle; F is the fluid sensitivity factor; BI is the brittleness index; ρ is the density; ΔF is the change in the formation fluid sensitivity factor on both sides of the interface; ΔBI is the change in the formation brittleness index on both sides of the interface; Δρ is the change in the formation density on both sides of the interface; V S is the shear wave velocity, V P is the longitudinal wave velocity; The first seismic record S(θ) construction module is used to obtain the approximate formula of the first seismic record S(θ) according to the expression of the formation reflection coefficient as follows: S(θ)=K F (θ)WDL F +K BI (θ)WDL BI +K ρ (θ)WDL ρ Where, S(θ) is the first earthquake record; θ is the incident angle; W is the wavelet matrix; D is the difference matrix; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; V S is the shear wave velocity, V P is the longitudinal wave velocity; The forward equation construction module is used to simplify the approximate formula of the first seismic record S(θ) and obtain the forward equation as follows: d=Gr Where d is the second earthquake record; G is the forward operator; r = [L F ,L BI ,L ρ ] T is the elastic parameter to be inverted; L F =lnF, which is the logarithmic value of the fluid sensitivity factor F; L BI =lnBI, which is the logarithm of the brittleness index BI; L ρ =lnρ, which is the logarithm of density ρ; T is the transpose operation of the matrix; The inversion algorithm module obtains the inversion results of the fluid sensitivity factor F, the brittleness index BI and the density ρ through the inversion algorithm based on the forward equation.

9. Computer device, characterized in that The invention comprises a storage and a processor, wherein the storage stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that A computer program is stored, and when the computer program is executed by a processor, the processor is caused to perform the steps of the method according to any one of claims 1 to 7.

Images