Reservoir physical property parameter seismic probabilistic inversion method for inhomogeneous medium

Abstract

The invention provides a reservoir physical parameter seismic probabilistic inversion method for a non-uniform medium, and the method is combined with the posterior expression of Bayesian linear inversion to form a semi-analytical method for inverting the posterior probability of the reservoir physical parameter from an elastic parameter estimation result. In order to describe heterogeneity caused by a complex pore structure in a reservoir, a dual-pore medium Biot-Rayleigh model is introduced to establish a rock physics positive operator, an analytical expression of posterior distribution is formed under the assumption of a Gaussian mixture model, and meanwhile, in order to improve the reliability of probability distribution estimation, the probability distribution of the reservoir is estimated. More sufficient data samples are provided by adopting Monte Carlo simulation, and finally inversion of physical property parameters such as porosity, water saturation and soft pore content is realized. Actual well-side seismic trace data tests show the effectiveness of the method, and compared with a traditional method based on a uniform hypothesis model, the method has higher inversion precision.

Description

Technical Field

[0001] The invention belongs to the field of unconventional anisotropic reservoir seismic exploration, and in particular relates to a seismic probabilistic inversion method for reservoir physical property parameters of inhomogeneous media. Background Art

[0002] Rock physics builds a bridge between the physical parameters of rock skeleton, pores and elastic properties, while seismic inversion is an important technology for mining reservoir elastic information from seismic observation data. Rock physics-driven seismic inversion quantitatively extracts reservoir physical property information from seismic data and is gaining more and more attention in oil and gas exploration and evaluation.

[0003] According to the adopted rock physics model, rock physics inversion can be performed based on the theoretical rock elastic relationship. For example, the effective medium theory evaluates the overall elastic properties of the rock skeleton and pore space as a whole, including the inclusion model and the contact model. Another type is the two-phase medium theory that describes the problem of seismic wave propagation in fluid-containing porous rocks, among which the Gassmann equation is most commonly used to simulate the elastic response of fluid-saturated rocks at the seismic scale. Due to the feasibility of solving related inversion problems, commonly used rock physics models are applied under the assumption of rock homogeneity. In order to improve the effect of seismic rock physics inversion of complex reservoirs, many studies have tried and improved rock physics models. For example, in the rock physics modeling of carbonate rocks, pores are defined as hard pores, fractures and reference pores, and rock physics inversion is carried out to predict the distribution of pore types; the analytical form of the linear rock physics model is used to form the posterior probability expression of rock physics inversion under the Bayesian framework, and the recursive exact solution of the approximate posterior probability of seismic rock physics inversion. For example, some scholars have quantitatively calculated the linear forward operators of carbonate physical properties and elastic parameters, and realized the linear inversion of physical property parameters. These studies all adopted the Gassmann equation and its extended model based on the homogeneous assumption, which are not applicable to heterogeneous reservoirs with complex pore structures such as carbonate rocks. Therefore, the reservoir heterogeneity problem has not been well solved in the seismic rock physics inversion problem. Summary of the invention

[0004] In view of the problems existing in the prior art, the present invention provides a seismic probabilistic inversion method for reservoir physical parameters in heterogeneous media, which is combined with the posterior expression of Bayesian linear inversion to form a semi-analytical method for estimating the posterior probability of reservoir physical parameters from the inversion results of elastic parameters. In order to describe the heterogeneity caused by the complex pore structure inside the reservoir, the dual-porosity medium Biot-Rayleigh model is introduced to establish a rock physics forward operator, and an analytical expression of the posterior distribution is formed under the assumption of a Gaussian mixture model. At the same time, in order to improve the reliability of the probability distribution estimation, Monte Carlo simulation is used to provide more sufficient data samples, which can solve the reservoir heterogeneity problem in the seismic rock physics inversion problem.

[0005] In order to solve the above technical problems, the present invention provides the following technical solutions: a seismic probabilistic inversion method for reservoir physical property parameters of heterogeneous media, comprising the following steps:

[0006] S1. Based on the heterogeneous dual-porosity Biot-Rayleigh equation, a Biot-Rayleigh rock physics model is established to estimate the elastic parameters of the reservoir;

[0007] S2, using prestack AVO inversion to estimate elastic parameters based on the angle gathers of seismic data;

[0008] S3, calculate the relative content of soft pores in the well, and combine the measured saturation and porosity to form the physical property parameter y, and use the expectation maximization algorithm to calculate the statistical parameters of the prior distribution of the physical property parameter y;

[0009] S4. Establish a priori probability distribution of physical property parameters, perform Monte Carlo simulation to obtain expanded data of physical property parameters, estimate corresponding elastic parameters by establishing Biot-Rayleigh rock physics model, and establish joint data samples;

[0010] S5. Based on the joint samples, the expectation maximization algorithm is used to calculate the statistical parameters of the joint probability distribution, the error-related parameters are initialized, and the iterative cycle of physical property parameter inversion is performed;

[0011] S6. Based on the estimation results of the statistical parameters of the joint probability distribution, calculate the statistical parameters and the posterior mean of the posterior probability distribution;

[0012] S7, using the posterior mean, the corresponding elastic parameters are calculated by establishing the Biot-Rayleigh rock physics model, and the error between the measured elastic data and the predicted elastic data is calculated;

[0013] S8. Determine whether the error is higher than the error threshold. If yes, return to step S5. Otherwise, output the posterior mean as the inversion result.

[0014] Furthermore, the aforementioned step S1 specifically comprises: using the Voigt-Reuss-Hil average to calculate the elastic modulus of the mineral mixture, then using the isotropic differential effective medium model to add soft pores and hard pores to the rock matrix of the mineral mixture, and finally calculating the velocity of the saturated fluid rock based on the Biot-Rayleigh equation, including the longitudinal wave velocity V P and shear wave velocity V S ;

[0015] Among them, the Biot-Rayleigh equation is a wave propagation equation describing the dual porosity model, and its plane wave equation is:

[0016]

[0017] Where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 , b 11 , b 12 , b 13 , b 21 , b 22 , b 23 , b 31 , b 32 , b 33 is the wave equation coefficient;

[0018] The wave equation coefficients are as follows:

[0019] a 11 =A+2N+i(Q2φ1-Q1φ2)x1,a 12 =Q1+i(Q2φ1-Q1φ2)x2

[0020] a 13 =Q2+i(Q2φ1-Q1φ2)x3,a 21 =Q2-iR1φ2x3

[0021] a 22 =R1-iR1φ2x2, a 23 =-iR1φ2x3

[0022] a 31 =Q2+iR2φ1x1,a 32 =iR2φ1x2,a 33 =R2+iR2φ1x3

[0023] b 11=iω(b1+b2)-ρ 11 ω 2 , b 12 =-iωb1-ρ 11 ω 2

[0024] b 13 =-iωb2-ρ 13 ω 2 , b 21 =-iωb1-ρ 12 ω 2

[0025] b 22 =iωb1-ρ 22 ω 2 , b 23 =0, b 32 =0

[0026] b 31 =-iωb2-ρ 13 ω 2 , b 33 =-ρ 33 ω 2 +iωb2

[0027] Among them, ρ 11 , 12 , 13 is the density parameter of the main phase medium, ρ 22 and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary number sign, x1, x2, x3 are attenuation factor parameters, A, N, Q1, Q2, R1, R2 are Biot stiffness coefficients.

[0028] Furthermore, the inversion objective function of the aforementioned step S2 for estimating elastic parameters using prestack AVO inversion is as follows:

[0029] J(x)=||d obs -G(x)||2+λ(x-μ x ) Τ ·(Σ x )·(x-μ x )

[0030] Among them, d obs is the observed seismic angle gather; x is the elastic parameter vector, including V P 、V S and ρ; G is the forward operator from elastic parameters to seismic data; μ x and Σ x are the expected matrix and covariance matrix of x respectively, and λ is the adjustment parameter.

[0031] Furthermore, the aforementioned step S3 includes the following sub-steps:

[0032] S3.1. Calculate the relative content of soft pores in the well as follows:

[0033] r sp =φ sp / φ

[0034] Among them, r sp is the relative content of soft pores, φ sp is the porosity of soft pores, and φ represents the porosity.

[0035] S3.2. Calculate the statistical parameters of the prior distribution of the physical property parameter y using the expectation maximization algorithm. Specifically, set the prior distribution function of the physical property parameter to a two-component Gaussian mixture model, and use the expectation maximization method to estimate the statistical parameters of the prior distribution, including the weight parameters, expectations and covariance matrices of different Gaussian components;

[0036] The prior distribution function of the physical property parameters is set to a two-component Gaussian mixture model, as follows:

[0037] P prior (y) = β1N1(y; μ 1|y ,Σ 1|y )+β2Ν2(y;μ 2|y ,Σ 2|y )

[0038] Among them, P prior (y) is the prior distribution function of the physical property parameters, β1 and β2 are the weight parameters in the prior distribution function, N1 and N2 represent the first and second Gaussian components, μ 1|y and Σ 1|y are the expectation and covariance of the first Gaussian component in the prior distribution function, μ 2|y and Σ 2|y are the expectation and covariance of the second Gaussian component in the prior distribution function, y is the physical parameter vector; β1, β2, μ 1|y , μ 2|y ,Σ 1|y and Σ 2|y All of them are estimated using the expectation maximization method.

[0039] Furthermore, the aforementioned step S4 specifically includes: obtaining the physical property parameter expansion data y by Monte Carlo simulation based on the prior distribution estimation result of the physical property parameter mc , and expand the physical property parameter data y mc As a known value, the Biot-Rayleigh rock physics model is used to calculate the value that is consistent with the data y mc Corresponding elastic parameter expansion data xmc , establish a joint data sample {x mc ,y mc}.

[0040] Furthermore, the aforementioned step S5 includes the following sub-steps:

[0041] S5.1. Assuming the joint probability distribution to be a two-component Gaussian mixture model, the joint probability distribution function expression is as follows:

[0042]

[0043] Among them, x mc and mc represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, P joint (x mc ,y mc ) is the joint probability distribution of elastic-physical parameters, β1 and β2 are weight parameters in the joint distribution function, N1 and N2 represent the first and second Gaussian components, and are the expectation and covariance of the first Gaussian component in the joint distribution function, and are the expectation and covariance of the second Gaussian component in the joint distribution function, β1, β2, and All are estimated by the expectation maximization method;

[0044] and The mathematical forms are as follows:

[0045]

[0046] Among them, μ 1|x and μ 1|y are the prior expectations of the elastic parameters and physical parameters of the first Gaussian component, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and is the cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data; the cross-covariance matrix is ​​as follows:

[0047]

[0048] Among them, δ1 represents the covariance between different parameters in the statistical parameters of the first Gaussian component, V P is the longitudinal wave velocity, V Sis the shear wave velocity, ρ is the density, φ is the porosity, S W is water saturation, r sp is the soft hole content; at the same time, the cross covariance matrix satisfies the following relationship:

[0049] and The mathematical forms are as follows:

[0050]

[0051] Among them, μ 2|x and μ 2|y are the prior expectations of the elastic parameters and physical parameters of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the elastic and physical parameter expansion data of the second Gaussian component is as follows:

[0052]

[0053] Among them, δ2 represents the covariance between different parameters in the statistical parameters of the second Gaussian component, and the cross covariances satisfy the following relationship:

[0054] S5.2. Initialize error-related parameters, including data error ε and data error variance Σ ε , starting the iterative cycle of physical property parameter inversion.

[0055] Furthermore, the aforementioned step S6 includes the following sub-steps:

[0056] S6.1. The posterior probability distribution of the physical property parameters is described by a two-component Gaussian mixture model, expressed as follows:

[0057] P post (y|x)=α1N1(y; μ 1|(y|x) ,Σ 1|(y|x) )+α2N2(y;μ 2|(y|x) ,Σ 2|(y|x) )

[0058] Among them, x is the elastic parameter, y is the physical parameter, P post (y|x) is the posterior probability function of the physical parameter y when the elastic parameter x is known, α1 and α2 are the weights of different Gaussian distributions, N1 and N2 represent two Gaussian components, μ 1|(y|x) and Σ1|(y|x) represents the expectation and covariance of the first Gaussian component in the posterior probability, μ 2|(y|x) and Σ 2|(y|x) represents the expectation and covariance of the second Gaussian component in the posterior probability;

[0059] S6.2. Calculate the expectation and covariance corresponding to the first Gaussian component in the posterior probability function as follows:

[0060]

[0061] Among them, μ 1|y is the prior expectation of the physical parameters in the statistical parameters of the first Gaussian component, μ 1|x is the prior expectation of the elastic parameter in the statistical parameter of the first Gaussian component, x is the elastic parameter result, x mc and mc They represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, Σ ε is the preset data error covariance, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data;

[0062] The expectation and covariance corresponding to the second Gaussian component in the posterior probability function are calculated as follows:

[0063]

[0064] Among them, μ 2|y is the prior expectation of the physical parameters in the statistical parameters of the second Gaussian component, μ 2|x is the prior expectation of the elastic parameter in the statistical parameter of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data;

[0065] S6.3. After estimating the expectation of different Gaussian components in the posterior probability distribution, the posterior expectation of the physical property parameters is calculated, and its mathematical expression is as follows:

[0066] μ PE =α1μ 1|(y|x) +α2μ 2|(y|x)

[0067] Among them, μ PE is the posterior expectation of the physical property parameters.

[0068] Furthermore, the aforementioned step S7 includes the following sub-steps:

[0069] S7.1. Using the posterior expectation μ of physical property parameters PE , the Biot-Rayleigh rock physics model is used to predict the elastic parameters as follows:

[0070]

[0071] Among them, x pred Represents the predicted elastic parameters, including the predicted longitudinal wave velocity Predicted shear wave velocity and the predicted density ρ pred , F BR It is a forward operator based on the Biot-Rayleigh rock physics model;

[0072] S7.2. Calculate the error between the predicted speed and the input data as follows:

[0073]

[0074] Where S represents the error function, and are the input data of P-wave velocity and S-wave velocity, κ1, κ2 and κ3 are the weight parameters of P-wave velocity, S-wave velocity and density data items.

[0075] Furthermore, in the aforementioned step S8, the posterior mean is as follows:

[0076]

[0077] Among them, y inv is the inversion result of physical property parameters, i represents the i-th iteration, Represents the posterior expectation of the physical property parameters estimated in the i-th iteration.

[0078] Compared with the prior art, the beneficial technical effects of the above technical solutions adopted by the present invention are as follows: the previous reservoir physical property parameter inversion adopts the rock physics theory of uniform medium assumption, such as Gassmann model, etc., which cannot describe the heterogeneity caused by complex pore structure or patchy fluid in complex reservoirs, which will lead to limited accuracy of reservoir prediction. The present invention introduces the Biot-Rayleigh model of heterogeneous medium into the reservoir physical property parameter inversion to improve the simulation accuracy of rock physics forward modeling. In addition, the rock physics model in the previous reservoir physical property parameter inversion is relatively simple and nonlinear, so Bayesian or least squares linear inversion methods are mostly used. However, since the rock physics forward modeling operator based on the Biot-Rayleigh model is very nonlinear, in order to improve efficiency and ensure accuracy at the same time, the present invention proposes a semi-analytical probabilistic inversion method based on Biot-Rayleigh rock physics forward modeling, and in order to ensure the accuracy of statistical parameter estimation, it is proposed to combine Biot-Rayleigh rock physics modeling and Monte Carlo simulation to provide more sufficient data samples. Compared with conventional inversion of uniform medium physical property parameters, the present invention has higher inversion accuracy and is more suitable for prediction of complex heterogeneous reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 It is a flow chart of the method of the present invention.

[0080] Figure 2 It is a prior distribution diagram of physical property parameters based on a two-component Gaussian mixture model. (a) is the distribution result diagram of water saturation and porosity, and (b) is the distribution result diagram of porosity and soft pore content.

[0081] Figure 3 It is a joint probability distribution diagram based on a two-component Gaussian mixture model and Monte Carlo simulation. (a) is the joint distribution diagram of P-wave velocity and porosity, and (b) is the joint distribution diagram of S-wave velocity and porosity.

[0082] Figure 4 It is a schematic diagram of the elastic parameter inversion results of the wellbore seismic data. In the figure, (a) is a schematic diagram of the longitudinal wave velocity inversion results, (b) is a schematic diagram of the shear wave velocity inversion results, and (c) is a schematic diagram of the density inversion results.

[0083] Figure 5 These are the probabilistic inversion diagrams of physical parameters based on the Gassmann rock physics model. (a) is the probabilistic inversion diagram of porosity, and (b) is the probabilistic inversion diagram of shear wave velocity.

[0084] Figure 6These are the probabilistic inversion diagrams of physical property parameters based on the Biot-Rayleigh rock physics model. (a) is the probabilistic inversion diagram of porosity, (b) is the probabilistic inversion diagram of water saturation, and (c) is the probabilistic inversion diagram of soft pore content. DETAILED DESCRIPTION

[0085] In order to better understand the technical content of the present invention, specific embodiments are given and described as follows in conjunction with the accompanying drawings.

[0086] Various aspects of the invention are described herein with reference to the accompanying drawings, in which many illustrative embodiments are shown. The embodiments of the invention are not limited to those described in the accompanying drawings. It should be understood that the invention is implemented by any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed in the invention are not limited to any implementation. In addition, some aspects disclosed in the invention may be used alone or in any appropriate combination with other aspects disclosed in the invention.

[0087] refer to Figure 1 The present invention provides a seismic probabilistic inversion method for reservoir physical property parameters of heterogeneous media, comprising the following steps:

[0088] S1. Based on the heterogeneous dual-porosity Biot-Rayleigh equation, a Biot-Rayleigh rock physics model is established to estimate the elastic parameters of the reservoir;

[0089] S2, using prestack AVO inversion to estimate elastic parameters based on the angle gathers of seismic data;

[0090] S3, calculate the relative content of soft pores in the well, and combine the measured saturation and porosity to form the physical property parameter y, and use the expectation maximization algorithm to calculate the statistical parameters of the prior distribution of the physical property parameter y;

[0091] S4. Establish a priori probability distribution of physical property parameters, perform Monte Carlo simulation to obtain expanded data of physical property parameters, use Biot-Rayleigh rock physics model to estimate corresponding elastic parameters, and establish joint data samples;

[0092] S5. Based on the joint samples, the expectation maximization algorithm is used to calculate the statistical parameters of the joint probability distribution, the error-related parameters are initialized, and the iterative cycle of physical property parameter inversion is performed;

[0093] S6. Based on the estimation results of the statistical parameters of the joint probability distribution, calculate the statistical parameters and the posterior mean of the posterior probability distribution;

[0094] S7, using the posterior mean, the Biot-Rayleigh rock physics model is used to calculate the corresponding elastic parameters, and the error between the measured elastic data and the predicted elastic data is calculated;

[0095] S8. Determine whether the error is higher than the error threshold. If yes, return to step S5. Otherwise, output the posterior mean as the inversion result.

[0096] As a preferred embodiment of the present invention, in step S1, first, the elastic modulus of the mineral mixture is calculated using the Voigt-Reuss-Hill (VRH) average method. The VRH average method is used to establish the mineral mixture, and the mineral mixture stiffness matrix C is calculated by the following formula:

[0097]

[0098] in,

[0099]

[0100] Where μ is the shear modulus of the mineral mixture, K is the bulk modulus of the brittle mineral mixture, C is the stiffness matrix of the mineral mixture, and μ is V and μ R is the VRH average upper limit shear modulus and VRH average lower limit shear modulus of the brittle mineral mixture, K V and K R are the VRH average upper limit bulk modulus and VRH average lower limit bulk modulus of the brittle mineral mixture.

[0101] Then, the soft and hard pores are added to the rock matrix of the mineral mixture using an isotropic differential effective medium (DEM) model.

[0102] Finally, the velocity of the fluid-saturated rock, including the P-wave velocity and S-wave velocity, is calculated based on the Biot-Rayleigh equation.

[0103] The Biot-Rayleigh equation is a wave propagation equation that describes the dual-porosity model. Its plane wave equation is:

[0104]

[0105] Where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 、b 11 、b 12 、b 13 、b 21 、b 22 、b 23, b 31 , b 32 , b 33 are the coefficients of the wave equation. The specific expressions of the wave equation coefficients are as follows:

[0106] a 11 = A + 2N + i(Q2φ1 - Q1φ2)x1, a 12 = Q1 + i(Q2φ1 - Q1φ2)x2

[0107] a 13 = Q2 + i(Q2φ1 - Q1φ2)x3, a 21 = Q2 - iR1φ2x3

[0108] a 22 = R1 - iR1φ2x2, a 23 = -iR1φ2x3

[0109] a 31 = Q2 + iR2φ1x1, a 32 = iR2φ1x2, a 33 = R2 + iR2φ1x3

[0110] b 11 = iω(b1 + b2) - ρ 11 ω 2 , b 12 = -iωb1 - ρ 11 ω 2

[0111] b 13 = -iωb2 - ρ 13 ω 2 , b 21 = -iωb1 - ρ 12 ω 2

[0112] b 22 = iωb1 - ρ 22 ω 2 , b 23 = 0, b 32 = 0

[0113] b 31 = -iωb2 - ρ 13 ω 2 , b 33 = -ρ 33 ω 2 + iωb2

[0114] where ρ 11 , ρ 12 , ρ 13 are the density parameters of the main-phase medium, ρ 22and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary number sign, x1, x2, x3 are attenuation factor parameters, A, N, Q1, Q2, R1, R2 are Biot stiffness coefficients.

[0115] As a preferred embodiment of the present invention, the probabilistic inversion method for physical parameters proposed in step S2 is an indirect inversion strategy, which requires first obtaining elastic parameter results using pre-stack seismic gather inversion, and then estimating the physical parameters from the elastic parameter results using probabilistic inversion based on the Biot-Rayleigh rock physics model.

[0116] Here, the prestack AVA inversion based on the Zoeppritz equation is used to obtain the inversion results of the elastic parameters. The forward expression of the Zoeppritz equation is as follows:

[0117] G(x)=W·R PP

[0118] Among them, G is the forward operator from elastic parameters to seismic data, W is the wavelet matrix, and R PP is the reflection coefficient. The expression of reflection coefficient is as follows:

[0119]

[0120] Among them, VP1, VS1, ρ1 represent the longitudinal wave velocity, shear wave velocity and density of the upper medium respectively; VP2, VS2, ρ2 represent the longitudinal wave velocity, shear wave velocity and density of the lower medium respectively; They represent the longitudinal wave incident angle, transverse wave incident angle, longitudinal wave reflection angle, and transverse wave reflection angle, respectively. PP , R PS 、T PP 、T PS They represent the PP wave reflection coefficient, PS wave reflection coefficient, PP wave transmission coefficient and PS wave transmission coefficient respectively.

[0121] Based on the seismic observation data, conventional pre-stack AVA inversion is used to obtain the elastic parameters x of the rock formation. These parameters include the longitudinal wave velocity V P , shear wave velocity V S and density ρ. The inversion objective function is as follows:

[0122] J(x)=||d obs -G(x)||2+λ(x-μ x ) Τ ·(Σ x )·(x-μ x )

[0123] Among them, d obsis the observed seismic angle gather; x is the elastic parameter vector, including V P 、V S and ρ; G is the forward operator from elastic parameters to seismic data; where μ x and Σ x are the expected matrix and covariance matrix of x respectively, and λ is the adjustment parameter.

[0124] Figure 4 The elastic parameter results obtained by prestack AVA inversion using wellside seismic gathers. In the figure, (a) shows the P-wave velocity inversion results, (b) shows the S-wave velocity inversion results, and (c) shows the density inversion results. The black is the logging data, the blue is the initial model, and the red is the inversion result. The subsequent physical parameter estimation uses the elastic parameter inversion results (red) as known input data.

[0125] As a preferred embodiment of the present invention, step S3 includes the following sub-steps:

[0126] S3.1. Calculate the relative content of soft pores in the well as follows:

[0127] r sp =φ sp / φ

[0128] Among them, r sp is the relative content of soft pores, φ sp is the porosity of soft pores, and φ represents the porosity.

[0129] S3.2. Calculate the statistical parameters of the prior distribution of the physical property parameter y using the expectation maximization algorithm. Specifically, set the prior distribution function of the physical property parameter to a two-component Gaussian mixture model, and use the expectation maximization method to estimate the statistical parameters of the prior distribution, including the weight parameters, expectations and covariance matrices of different Gaussian components;

[0130] The prior distribution function of the physical property parameters is set to a two-component Gaussian mixture model, as follows:

[0131] P prior (y) = β1N1(y; μ 1|y ,Σ 1|y )+β2Ν2(y;μ 2|y ,Σ 2|y )

[0132] Among them, P prior (y) is the prior distribution function of the physical property parameters, β1 and β2 are the weight parameters in the prior distribution function, N1 and N2 represent the first and second Gaussian components, μ 1|y and Σ 1|y are the expectation and covariance of the first Gaussian component in the prior distribution function, μ 2|y and Σ2|y are the expectation and covariance of the second Gaussian component in the prior distribution function, y is the physical parameter vector; β1, β2, μ 1|y , μ 2|y ,Σ 1|y and Σ 2|y As a preferred embodiment of the present invention, step S4 is specifically: based on the prior distribution estimation result of the physical property parameters, the physical property parameter expansion data y is obtained by Monte Carlo simulation. mc , and expand the physical property parameter data y mc As a known value, the Biot-Rayleigh rock physics model is used to calculate the value that is consistent with the data y mc Corresponding elastic parameter expansion data x mc , thereby establishing a joint data sample {x mc ,y mc}. Figure 2 The results of the two-dimensional prior distribution of physical property parameters are shown. (a) is the distribution result of water saturation and porosity, and (b) is the distribution result of porosity and soft pore content. The prior distribution of physical property parameters is set as a two-component Gaussian mixture distribution. The scattered points are the physical property parameter expansion data obtained by Monte Carlo simulation based on well logging data. The colored curve is the estimated result of the prior distribution. It can be seen that the distribution estimation result is reasonable and can better describe the data characteristics.

[0133] As a preferred embodiment of the present invention, based on the obtained joint data sample, it can be used to calculate the joint probability distribution of elasticity-physical property parameters. Step S5 includes the following sub-steps:

[0134] S5.1. Assuming the joint probability distribution to be a two-component Gaussian mixture model, the joint probability distribution function expression is as follows:

[0135]

[0136] Among them, x mc and mc represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, P joint (x mc ,y mc ) is the joint probability distribution of elastic-physical parameters, β1 and β2 are weight parameters in the joint distribution function, N1 and N2 represent the first and second Gaussian components, and are the expectation and covariance of the first Gaussian component in the joint distribution function, and are the expectation and covariance of the second Gaussian component in the joint distribution function, β1, β2, and All are estimated by the expectation maximization method;

[0137] and The mathematical forms are as follows:

[0138]

[0139] Among them, μ 1|x and μ 1|y are the prior expectations of the elastic parameters and physical parameters of the first Gaussian component, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and is the cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data; the cross-covariance matrix is ​​as follows:

[0140]

[0141] Among them, δ1 represents the covariance between different parameters in the statistical parameters of the first Gaussian component, V P is the longitudinal wave velocity, V S is the shear wave velocity, ρ is the density, φ is the porosity, S W is water saturation, r sp is the soft hole content; at the same time, the cross covariance matrix satisfies the following relationship:

[0142] and The mathematical forms are as follows:

[0143]

[0144] Among them, μ 2|x and μ 2|y are the prior expectations of the elastic parameters and physical parameters of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the elastic and physical parameter expansion data of the second Gaussian component is as follows:

[0145]

[0146]

[0147] Among them, δ2 represents the covariance between different parameters in the statistical parameters of the second Gaussian component, and the cross covariances satisfy the following relationship:

[0148] S5.2. Initialize error-related parameters, including data error ε and data error variance Σ ε , starting the iterative cycle of physical property parameter inversion.

[0149] Figure 3 The results of the two-dimensional joint probability distribution of elastic-physical parameters are shown. (a) shows the joint distribution of P-wave velocity and porosity, and (b) shows the joint distribution of S-wave velocity and porosity. The joint probability is set as a two-component Gaussian mixture distribution, and the scatter points are the joint sample data of elastic-physical parameters obtained by Monte Carlo simulation based on logging data and the Biot-Rayleigh rock physics model. The colored curve is the estimated result of the joint distribution. It can be seen that the distribution estimation result is reasonable and can better describe the data characteristics.

[0150] As a preferred embodiment of the present invention, step S6 includes the following sub-steps:

[0151] S6.1. The posterior probability distribution of the physical property parameters is described by a two-component Gaussian mixture model, expressed as follows:

[0152] P post (y|x)=α1N1(y; μ 1|(y|x) ,Σ 1|(y|x) )+α2N2(y;μ 2|(y|x) ,Σ 2|(y|x) )

[0153] Among them, x is the elastic parameter, y is the physical parameter, P post (y|x) is the posterior probability function of the physical parameter y when the elastic parameter x is known, α1 and α2 are the weights of different Gaussian distributions, N1 and N2 represent two Gaussian components, μ 1|(y|x) and Σ 1|(y|x) represents the expectation and covariance of the first Gaussian component in the posterior probability, μ 2|(y|x) and Σ 2|(y|x) represents the expectation and covariance of the second Gaussian component in the posterior probability;

[0154] The weight parameter α1 is as follows:

[0155]

[0156] Among them, β1 and β2 are the weights of different Gaussian components in the prior distribution of elastic parameters, μ 1|x and Σ 1|x represents the expectation and covariance of the first Gaussian component in the prior distribution of the elastic parameters, μ2|x and Σ 2|x represents the expectation and covariance of the second Gaussian component in the prior distribution of the elasticity parameter;

[0157] The weight parameter α2 is as follows:

[0158]

[0159] S6.2. Calculate the expectation and covariance corresponding to the first Gaussian component in the posterior probability function as follows:

[0160]

[0161] Among them, μ 1|y is the prior expectation of the physical parameters in the statistical parameters of the first Gaussian component, μ 1|x is the prior expectation of the elastic parameter in the statistical parameter of the first Gaussian component, x is the elastic parameter result, x mc and mc They represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, Σ ε is the preset data error covariance, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data;

[0162] The expectation and covariance corresponding to the second Gaussian component in the posterior probability function are calculated as follows:

[0163]

[0164] Among them, μ 2|y is the prior expectation of the physical parameters in the statistical parameters of the second Gaussian component, μ 2|x is the prior expectation of the elastic parameter in the statistical parameter of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data;

[0165] S6.3. After estimating the expectations of different Gaussian components in the posterior probability distribution, the posterior expectations of the physical property parameters are calculated. The expression is as follows:

[0166] μ PE =α1μ1|(y|x) +α2μ 2|(y|x)

[0167] Among them, μ PE is the posterior expectation of the physical property parameters.

[0168] As a preferred embodiment of the present invention, step S7 includes the following sub-steps:

[0169] S7.1. Using the posterior expectation μ of physical property parameters PE , the Biot-Rayleigh rock physics model is used to predict the elastic parameters as follows:

[0170]

[0171] Among them, x pred represents the predicted elastic parameters, including the predicted longitudinal wave velocity V P pred , predict the shear wave velocity V S pred and the predicted density ρ pred , F BR It is a forward operator based on the Biot-Rayleigh rock physics model;

[0172] S7.2. Using the calculated elastic parameters and the measured elastic parameters in the well, the error of the elastic parameters can be obtained, thereby quality controlling the estimation accuracy of the physical property parameters. The error between the predicted velocity and the input data is calculated as follows:

[0173]

[0174] Where S represents the error function, and are the input data of P-wave velocity and S-wave velocity, κ1, κ2 and κ3 are the weight parameters of P-wave velocity, S-wave velocity and density data items.

[0175] The error value is calculated by the error function. If the error value is higher than the preset error threshold, the iterative process will be repeated until the error is lower than the error threshold, and the posterior mean is output as the inversion result:

[0176] As a preferred embodiment of the present invention, in step S8, the posterior mean is as follows:

[0177]

[0178] Among them, y inv is the inversion result of physical property parameters, i represents the i-th iteration, Represents the posterior expectation of the physical property parameters estimated in the i-th iteration.

[0179] Figure 5This is the probabilistic inversion result of the physical parameters of the conventional Gassmann rock physics modeling. Since Gassmann does not consider different pore types, the inversion results only include porosity and saturation. (a) shows the probabilistic inversion result of porosity, and (b) shows the probabilistic inversion result of shear wave velocity. Black is the measured data, green is the uncertainty interval, the red curve is the inversion result, the yellow area corresponds to the high probability value, and the blue area corresponds to the low probability value. It can be seen that the results based on the Gassmann model have large errors. Figure 6 The results of the probabilistic inversion of physical parameters based on the Biot-Rayleigh model proposed in the present invention are shown, (a) is the probabilistic inversion diagram of porosity, (b) is the probabilistic inversion diagram of water saturation, and (c) is the probabilistic inversion diagram of soft pore content. It can be seen that the results of this method have greatly improved the accuracy compared with the Gassmann method, and are highly consistent with the measured data, which illustrates the effectiveness of the method of the present invention.

[0180] Although the present invention has been described above with preferred embodiments, it is not intended to limit the present invention. A person skilled in the art of the present invention may make various modifications and improvements without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be determined by the definition of the claims.

Claims

1. A seismic probabilistic inversion method for reservoir physical property parameters in heterogeneous media, characterized in that: The following steps are involved: S1. Based on the heterogeneous dual-porosity Biot-Rayleigh equation, a Biot-Rayleigh rock physics model is established to estimate the elastic parameters of the reservoir; S2, using prestack AVO inversion to estimate elastic parameters based on the angle gathers of seismic data; S3, calculate the relative content of soft pores in the well, and combine the measured saturation and porosity to form the physical property parameter y, and use the expectation maximization algorithm to calculate the statistical parameters of the prior distribution of the physical property parameter y; S4. Establish a priori probability distribution of physical property parameters, perform Monte Carlo simulation to obtain expanded data of physical property parameters, estimate corresponding elastic parameters by establishing Biot-Rayleigh rock physics model, and establish joint data samples; S5. Based on the joint samples, the expectation maximization algorithm is used to calculate the statistical parameters of the joint probability distribution, the error-related parameters are initialized, and the iterative cycle of physical property parameter inversion is performed; S6. Based on the estimation results of the statistical parameters of the joint probability distribution, calculate the statistical parameters and the posterior mean of the posterior probability distribution; S7, using the posterior mean, the corresponding elastic parameters are calculated by establishing the Biot-Rayleigh rock physics model, and the error between the measured elastic data and the predicted elastic data is calculated; S8. Determine whether the error is higher than the error threshold. If yes, return to step S5. Otherwise, output the posterior mean as the inversion result.

2. A seismic probabilistic inversion method for reservoir physical property parameters of heterogeneous media according to claim 1, characterized in that: Step S1 specifically includes: using Voigt-Reuss-Hil average to calculate the elastic modulus of the mineral mixture, then using the isotropic differential effective medium model to add soft pores and hard pores to the rock matrix of the mineral mixture, and finally, calculating the velocity of the saturated fluid rock based on the Biot-Rayleigh equation, including the P-wave velocity V P and shear wave velocity V S ; Among them, the Biot-Rayleigh equation is a wave propagation equation describing the dual porosity model, and its plane wave equation is: Where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 、b 11 、b 12 、b 13 、b 21 、b 22 、b 23 、b 31 、b 32 、b 33 is the wave equation coefficient; The wave equation coefficients are as follows: a 11 =A+2N+i(Q2φ1-Q1φ2)x1,a 12 =Q1+i(Q2φ1-Q1φ2)x2 a 13 =Q2+i(Q2φ1-Q1φ2)x3,a 21 =Q2-iR1φ2x3 a 22 =R1-iR1φ2x2,a 23 =-iR1φ2x3 a 31 =Q2+iR2φ1x1,a 32 =iR2φ1x2,a 33 =R2+iR2φ1x3 b 11 =iω(b1+b2)-ρ 11 oh 2 ,b 12 =-iωb1-ρ 11 oh 2 b 13 =-iωb2-ρ 13 oh 2 ,b 21 =-iωb1-ρ 12 oh 2 b 22 =iωb1-ρ 22 oh 2 ,b 23 =0,b 32 =0 b 31 =-iωb2-ρ 13 oh 2 ,b 33 =-ρ 33 oh 2 +iωb2 Among them, ρ 11 , 12 , 13 is the density parameter of the main phase medium, ρ 22 and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary number sign, x1, x2, x3 are attenuation factor parameters, A, N, Q1, Q2, R1, R2 are Biot stiffness coefficients.

3. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: In step S2, the inversion objective function of estimating elastic parameters using prestack AVO inversion is as follows: J(x)=||d obs -G(x)||2+λ(x-μ x ) Τ ·(S x )·(x-μ x ) Among them, d obs is the observed seismic angle gather; x is the elastic parameter vector, including V P 、V S and ρ; G is the forward operator from elastic parameters to seismic data; μ x and Σ x are the expected matrix and covariance matrix of x respectively, and λ is the adjustment parameter.

4. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: Step S3 includes the following sub-steps: S3.

1. Calculate the relative content of soft pores in the well as follows: r sp =φ sp / f Among them, r sp is the relative content of soft pores, φ sp is the porosity of soft pores, and φ represents the porosity. S3.

2. Calculate the statistical parameters of the prior distribution of the physical property parameter y using the expectation maximization algorithm. Specifically, set the prior distribution function of the physical property parameter to a two-component Gaussian mixture model, and use the expectation maximization method to estimate the statistical parameters of the prior distribution, including the weight parameters, expectations and covariance matrices of different Gaussian components; The prior distribution function of the physical property parameters is set to a two-component Gaussian mixture model, as follows: P prior (y)=β1Ν1(y;μ 1|y ,S 1|y )+β2Ν2(y;μ 2|y ,S 2|y ) Among them, P prior (y) is the prior distribution function of the physical property parameters, β1 and β2 are the weight parameters in the prior distribution function, N1 and N2 represent the first and second Gaussian components, μ 1|y and Σ 1|y are the expectation and covariance of the first Gaussian component in the prior distribution function, μ 2|y and Σ 2|y are the expectation and covariance of the second Gaussian component in the prior distribution function, y is the physical parameter vector; β1, β2, μ 1|y , μ 2|y ,Σ 1|y and Σ 2|y All of them are estimated using the expectation maximization method.

5. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: Step S4 specifically includes: obtaining the physical property parameter expansion data y based on the prior distribution estimation result of the physical property parameter through Monte Carlo simulation mc , and expand the physical property parameter data y mc As a known value, the Biot-Rayleigh rock physics model is used to calculate the value that is consistent with the data y mc Corresponding elastic parameter expansion data x mc , establish a joint data sample {x mc ,y mc }.

6. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: Step S5 includes the following sub-steps: S5.

1. Assuming the joint probability distribution to be a two-component Gaussian mixture model, the joint probability distribution function expression is as follows: Among them, x mc and mc represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, P joint (x mc ,y mc ) is the joint probability distribution of elastic-physical parameters, β1 and β2 are weight parameters in the joint distribution function, N1 and N2 represent the first and second Gaussian components, and are the expectation and covariance of the first Gaussian component in the joint distribution function, and are the expectation and covariance of the second Gaussian component in the joint distribution function, β1, β2, and All are estimated by the expectation maximization method; and The mathematical forms are as follows: Among them, μ 1|x and μ 1|y are the prior expectations of the elastic parameters and physical parameters of the first Gaussian component, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and is the cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data; the cross-covariance matrix is ​​as follows: Among them, δ1 represents the covariance between different parameters in the statistical parameters of the first Gaussian component, V P is the longitudinal wave velocity, V S is the shear wave velocity, ρ is the density, φ is the porosity, S W is water saturation, r sp is the soft hole content; at the same time, the cross covariance matrix satisfies the following relationship: and The mathematical forms are as follows: Among them, μ 2|x and μ 2|y are the prior expectations of the elastic parameters and physical parameters of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the elastic and physical parameter expansion data of the second Gaussian component is as follows: Among them, δ2 represents the covariance between different parameters in the statistical parameters of the second Gaussian component, and the cross covariances satisfy the following relationship: S5.

2. Initialize error-related parameters, including data error ε and data error variance Σ ε , starting the iterative cycle of physical property parameter inversion.

7. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: Step S6 includes the following sub-steps: S6.

1. The posterior probability distribution of the physical property parameters is described by a two-component Gaussian mixture model, expressed as follows: P post (y|x)=α1N1(y;μ 1|(y|x) ,S 1|(y|x) )+α2N2(y;μ 2|(y|x) ,S 2|(y|x) ) Among them, x is the elastic parameter, y is the physical parameter, P post (y|x) is the posterior probability function of the physical parameter y when the elastic parameter x is known, α1 and α2 are the weights of different Gaussian distributions, N1 and N2 represent two Gaussian components, μ 1|(y|x) and Σ 1|(y|x) represents the expectation and covariance of the first Gaussian component in the posterior probability, μ 2|(y|x) and Σ 2|(y|x) represents the expectation and covariance of the second Gaussian component in the posterior probability; S6.

2. Calculate the expectation and covariance corresponding to the first Gaussian component in the posterior probability function as follows: Among them, μ 1|y is the prior expectation of the physical parameters in the statistical parameters of the first Gaussian component, μ 1|x is the prior expectation of the elastic parameter in the statistical parameter of the first Gaussian component, x is the elastic parameter result, x mc and mc They represent the elastic parameters and physical property parameter expansion data obtained by Monte Carlo simulation, Σ ε is the preset data error covariance, represents the conditional covariance matrix of the first Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the first Gaussian component physical property parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data; The expectation and covariance corresponding to the second Gaussian component in the posterior probability function are calculated as follows: Among them, μ 2|y is the prior expectation of the physical parameters in the statistical parameters of the second Gaussian component, μ 2|x is the prior expectation of the elastic parameter in the statistical parameter of the second Gaussian component, represents the conditional covariance matrix of the second Gaussian component elastic parameter augmented data, represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data; S6.

3. After estimating the expectation of different Gaussian components in the posterior probability distribution, the posterior expectation of the physical property parameters is calculated, and its mathematical expression is as follows: m PE =α1μ 1|(y|x) +a2m 2|(y|x) Among them, μ PE is the posterior expectation of the physical property parameters.

8. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: Step S7 includes the following sub-steps: S7.

1. Using the posterior expectation μ of physical property parameters PE , the Biot-Rayleigh rock physics model is used to predict the elastic parameters as follows: Among them, x pred Represents the predicted elastic parameters, including the predicted longitudinal wave velocity Predicted shear wave velocity and the predicted density ρ pred , F BR It is a forward operator based on the Biot-Rayleigh rock physics model; S7.

2. Calculate the error between the predicted speed and the input data as follows: Where S represents the error function, and are the input data of P-wave velocity and S-wave velocity, κ1, κ2 and κ3 are the weight parameters of P-wave velocity, S-wave velocity and density data items.

9. The method for seismic probabilistic inversion of reservoir physical property parameters for heterogeneous media according to claim 1, characterized in that: In step S8, the posterior mean is as follows: Among them, y inv is the inversion result of physical property parameters, i represents the i-th iteration, Represents the posterior expectation of the physical property parameters estimated in the i-th iteration.

Images