A seismic probabilistic inversion method for reservoir property parameters of non-uniform medium

Abstract

The application provides a reservoir physical parameter seismic probabilistic inversion method for non-uniform medium, which is combined with a posteriori expression of Bayesian linear inversion to form a semi-analytical method for inverting reservoir physical parameter posterior probability from elastic parameter estimation results. In order to describe the non-uniformity caused by complex pore structure inside the reservoir, a Biot-Rayleigh model of double-pore medium is introduced to establish a rock physics forward operator, and an analytical expression of the posterior distribution is formed under the assumption of Gaussian mixture model. Meanwhile, in order to improve the reliability of the probability distribution estimation, Monte Carlo simulation is adopted to provide more sufficient data samples, and finally the inversion of physical parameters such as porosity, water saturation and soft pore content is realized. Actual well seismic trace data test shows that the method has higher inversion accuracy than the traditional method based on the uniform assumption model.

Description

Technical Field

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

[0002] Rock physics bridges the gap between physical parameters like rock structure and pore size and elastic properties, while seismic inversion is a key technique for extracting reservoir elastic information from seismic observation data. Rock physics-driven seismic inversion, which quantitatively extracts reservoir physical property information from seismic data, is gaining increasing attention in oil and gas exploration and evaluation.

[0003] Depending on the rock physics model employed, rock physics inversion can be performed based on theoretical rock-elastic relationships. For example, effective medium theory considers the rock skeleton and pore space as a whole to evaluate its overall elastic properties, including inclusion models and contact models. Another type of theory describes seismic wave propagation in fluid-filled porous rocks. 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. To improve the performance of seismic rock physics inversion for complex reservoirs, numerous studies have attempted to refine rock physics models. For example, in carbonate rock physics modeling, pores are defined as hard pores, fractures, and reference pores, allowing rock physics inversion to predict the distribution of pore types. Analytical forms of linear rock physics models are used to develop posterior probability expressions for rock physics inversion within a Bayesian framework, as well as recursive exact solutions to approximate posterior probabilities for seismic rock physics inversion. For example, some researchers have developed linear forward operators that quantify the physical and elastic parameters of carbonate rocks, achieving linear inversion of physical parameters. These studies all used the Gassmann equation and its extended model based on the homogeneous assumption, which is not applicable to heterogeneous reservoirs with complex pore structures such as carbonate rocks. Therefore, the issue of reservoir heterogeneity has not been well addressed in seismic rock physics inversion. Summary of the Invention

[0004] This invention addresses the challenges of the prior art by providing a probabilistic seismic inversion method for reservoir physical parameters in heterogeneous media. This method, combined with the posterior expression of Bayesian linear inversion, forms a semi-analytical method for estimating the posterior probability of reservoir physical parameters from elastic parameter inversion results. To describe the heterogeneity caused by the complex pore structure within the reservoir, a dual-porosity Biot-Rayleigh model is introduced to establish a rock physics forward operator. An analytical expression for the posterior distribution is then formed under the assumption of a Gaussian mixture model. To improve the reliability of the probability distribution estimate, Monte Carlo simulation is employed to provide a more comprehensive data sample, enabling a more robust solution to reservoir heterogeneity in seismic rock physics inversion.

[0005] To solve the above technical problems, the present invention provides the following technical solution: a seismic probabilistic inversion method for reservoir physical property parameters in 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, estimate elastic parameters using prestack AVO inversion based on seismic data angle gathers;

[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 sample;

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

[0011] S6. Calculate the statistical parameters and the posterior mean of the posterior probability distribution based on the estimation results of the statistical parameters of the joint probability distribution;

[0012] S7. Using the posterior mean, calculate the corresponding elastic parameters by establishing a Biot-Rayleigh rock physics model, and calculate the error between the measured elastic data and the predicted elastic data;

[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] Further, the aforementioned step S1 is specifically: using Voigt-Reuss-Hil average to calculate the elastic modulus of the mineral mixture, then adding soft pores and hard pores into the rock matrix of the mineral mixture by using the isotropic differential effective medium model, and finally calculating the velocity of the saturated fluid rock based on the Biot-Rayleigh equation, including the longitudinal wave velocity V P and the transverse wave velocity V S ;

[0015] wherein the Biot-Rayleigh equation is a wave propagation equation describing a double-porosity model, and the plane wave equation thereof is:

[0016]

[0017] wherein k is a 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 wave equation coefficients;

[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 elastic parameters estimated by prestack AVO inversion in step S2 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, assume the prior distribution function of the physical property parameter is 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, expectation, and covariance matrices of the 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 )+β2N2(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 are estimated using the expectation maximization method.

[0039] Furthermore, the aforementioned 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 of 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 is a two-component Gaussian mixture model, the joint probability distribution function is expressed as follows:

[0042]

[0043] Among them, x mc and y 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 elasticity and physical property 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 using 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 parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data 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 the 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 process of physical property parameter inversion.

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

[0056] S6.1. Assume that 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 respectively, μ 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 y 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 parameter expansion data, and 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 Cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data;

[0065] S6.3. After estimating the expectations of different Gaussian components in the posterior probability distribution, the posterior expectation of the physical property parameters is calculated. The 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 at the i-th iteration.

[0078] Compared to the prior art, the beneficial technical effects of the above technical solutions adopted by the present invention are as follows: Previous reservoir property parameter inversions used rock physics theories based on the homogeneous medium assumption, such as the Gassmann model, which are unable to describe the heterogeneity caused by complex pore structures or patchy fluids in complex reservoirs, resulting in limited reservoir prediction accuracy. The present invention introduces the Biot-Rayleigh model of heterogeneous media into the inversion of reservoir property parameters to improve the simulation accuracy of rock physics forward modeling. In addition, the rock physics models used in previous reservoir property parameter inversions were relatively simple and not highly nonlinear, and therefore Bayesian or least squares linear inversion methods were often used. However, since the rock physics forward modeling operators based on the Biot-Rayleigh model are highly nonlinear, in order to improve efficiency while ensuring accuracy, the present invention proposes a semi-analytical probabilistic inversion method based on the Biot-Rayleigh rock physics forward modeling. At the same time, 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 physical property parameters of uniform media, 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 priori distribution diagram of physical property parameters based on a two-component Gaussian mixture model. In the figure, (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 the two-component Gaussian mixture model and Monte Carlo simulation. In the figure, (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 trace 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 maps of physical property parameters based on the Biot-Rayleigh rock physics model. (a) is the probabilistic inversion map of porosity, (b) is the probabilistic inversion map of water saturation, and (c) is the probabilistic inversion map of soft pore content. DETAILED DESCRIPTION

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

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

[0087] refer to Figure 1 The present invention provides a seismic probabilistic inversion method for reservoir physical property parameters in 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, estimate elastic parameters using prestack AVO inversion based on seismic data angle gathers;

[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 the Biot-Rayleigh rock physics model to estimate the corresponding elastic parameters, and establish a joint data sample;

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

[0093] S6. Calculate the statistical parameters and the posterior mean of the posterior probability distribution based on the estimation results of the statistical parameters of the joint probability distribution;

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

[0095] S8, determining whether the error is higher than an error threshold, if yes, returning to execute step S5, otherwise outputting the posterior mean as the inversion result.

[0096] As a preferred embodiment of the present application, in step S1, first, the elastic modulus of the mineral mixture is calculated using Voigt-Reuss-Hill (V-R-H) average. Using the V-R-H average method, the mineral mixture is established, and the stiffness matrix C of the mineral mixture is calculated by the following formula:

[0097]

[0098] wherein,

[0099]

[0100] In the formula, μ 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, μ V and μ R are the V-R-H average upper limit shear modulus and the V-R-H average lower limit shear modulus of the brittle mineral mixture, K V and K R are the V-R-H average upper limit bulk modulus and the V-R-H average lower limit bulk modulus of the brittle mineral mixture.

[0101] Then, the isotropic differential effective medium (DEM) model is used to add soft pores and hard pores to the rock matrix of the mineral mixture.

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

[0103] The Biot-Rayleigh equation is a wave propagation equation describing the double-porosity model, and its plane wave equation is:

[0104]

[0105] wherein, 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 自 电影发动商发。 The fluctuation equation is the coefficient of the specific expression 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 oh 2 ,b 12 =-iωb1-ρ 11 oh 2

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

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

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

[0114] Among them, p 11 ,r 12 ,r 13 下载相电影 density parameter,r 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 physical parameters from the elastic parameter results using probabilistic inversion based on the Biot-Rayleigh rock physics model.

[0116] Here, we use pre-stack AVA inversion based on the Zoeppritz equation to obtain the inversion results of elastic parameters. The forward expression of the Zoeppritz equation is as follows:

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

[0118] Where 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 the reflection coefficient is as follows:

[0119]

[0120] Where VP1, VS1, and ρ1 represent the longitudinal wave velocity, shear wave velocity, and density of the upper layer, respectively; VP2, VS2, and ρ2 represent the longitudinal wave velocity, shear wave velocity, and density of the lower layer, 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 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 Elastic parameter results obtained using prestack AVA inversion on near-well seismic gathers. (a) shows the P-wave velocity inversion results, (b) shows the S-wave velocity inversion results, and (c) shows the density inversion results. Black represents the well data, blue represents the initial model, and red represents the inversion results. Subsequent physical property 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, assume the prior distribution function of the physical property parameter is 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, expectation, and covariance matrices of the 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 )+β2N2(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 results of the physical property parameters, the physical property parameter expansion data y is obtained 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 of mc Corresponding elastic parameter expansion data x mc , thereby establishing a joint data sample {x mc ,y mc}. Figure 2 The two-dimensional prior distribution results for physical property parameters are shown. (a) shows the distribution results for water saturation and porosity, and (b) shows the distribution results for porosity and soft pore content. The prior distribution for the physical property parameters is assumed to be a two-component Gaussian mixture distribution. The scattered points represent the physical property parameter augmentation data obtained through Monte Carlo simulation based on well logging data. The colored curves represent the estimated prior distribution results. It can be seen that the distribution estimates are reasonable and well describe the data characteristics.

[0133] As a preferred embodiment of the present invention, the obtained joint data sample 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 is a two-component Gaussian mixture model, the joint probability distribution function is expressed as follows:

[0135]

[0136] Among them, x mc and y 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 elasticity and physical property 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 using 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 parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data 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 the 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 process of physical property parameter inversion.

[0149] Figure 3 The two-dimensional joint probability distribution of elastic and physical parameters is 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 assumed to be a two-component Gaussian mixture distribution. The scatter plots represent sample data of the joint elastic and physical parameters obtained using Monte Carlo simulations based on well logging data and the Biot-Rayleigh rock physics model. The colored curves represent the estimated joint distribution, demonstrating that the distribution estimates are reasonable and well describe the data characteristics.

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

[0151] S6.1. Assume that 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 respectively, μ 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 parameter, μ2|x and∑ 2|x denote the mean and covariance of the second Gaussian component in the elastic parameter prior distribution;

[0157] the weight parameter a2 is given by:

[0158]

[0159] S6.2, the mean and covariance corresponding to the first Gaussian component in the posterior probability function are calculated as follows:

[0160]

[0161] where μ 1|y is the prior mean of the property parameter in the first Gaussian component statistical parameter, μ 1|x is the prior mean of the elastic parameter in the first Gaussian component statistical parameter, x is the elastic parameter result, x mc and y mc denote the elastic parameter and property parameter augmented data obtained by Monte Carlo simulation,∑ ε is the preset data error covariance, denotes the conditional covariance matrix of the elastic parameter augmented data of the first Gaussian component, denotes the conditional covariance matrix of the property parameter augmented data of the first Gaussian component, and is the cross covariance matrix of the elastic and property parameter augmented data of the first Gaussian component;

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

[0163]

[0164] where μ 2|y is the prior mean of the property parameter in the second Gaussian component statistical parameter, μ 2|x is the prior mean of the elastic parameter in the second Gaussian component statistical parameter, denotes the conditional covariance matrix of the elastic parameter augmented data of the second Gaussian component, denotes the conditional covariance matrix of the property parameter augmented data of the second Gaussian component, and is the cross covariance matrix of the elastic and property parameter augmented data of the first Gaussian component;

[0165] S6.3, after estimating the mean of different Gaussian components in the posterior probability distribution, the posterior mean of the property parameter is calculated, and the expression is as follows:

[0166] μ PE = a1 μ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 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 continued 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 at the i-th iteration.

[0179] Figure 5These are the probabilistic inversion results of physical parameters using conventional Gassmann rock physics modeling. Because Gassmann does not account for different pore types, the inversion results only include porosity and saturation. (a) shows the probabilistic inversion results for porosity, and (b) shows the probabilistic inversion results for shear wave velocity. Black represents measured data, green represents the uncertainty interval, and the red curve is the inversion result. Yellow areas correspond to high probabilities, while blue areas correspond to low probabilities. As can be seen, the results based on the Gassmann model have significant errors. Figure 6 The results of the probabilistic inversion of physical parameters based on the Biot-Rayleigh model proposed in this paper are shown. (a) is the probabilistic inversion plot of porosity, (b) is the probabilistic inversion plot of water saturation, and (c) is the probabilistic inversion plot of soft pore content. It can be seen that the results of this method have significantly improved accuracy compared to the Gassmann method and are highly consistent with the measured data, demonstrating the effectiveness of the method proposed in this paper.

[0180] While the present invention has been described above with reference to preferred embodiments, this is not intended to limit the present invention. Persons skilled in the art will readily appreciate that various modifications and variations can be made 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 appended claims.

Claims

1. A seismic probabilistic inversion method for reservoir physical parameters in heterogeneous media, characterized by: 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, estimate elastic parameters using prestack AVO inversion based on seismic data angle gathers; 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; The following sub-steps are included: 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, φ 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, assume the prior distribution function of the physical property parameter is 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, expectation, and covariance matrices of the 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 are estimated using the expectation maximization method; 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 sample; S5. Based on the joint samples, the expectation maximization algorithm is used to calculate the statistical parameters of the joint probability distribution, initialize the error-related parameters, and perform an iterative cycle of physical property parameter inversion; S6. Calculate the statistical parameters and the posterior mean of the posterior probability distribution based on the estimation results of the statistical parameters of the joint probability distribution; S7. Using the posterior mean, calculate the corresponding elastic parameters by establishing a Biot-Rayleigh rock physics model, and calculate the error between the measured elastic data and the predicted elastic data; 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. The seismic probabilistic inversion method for reservoir physical property parameters in heterogeneous media according to claim 1, characterized in that: Step S1 specifically involves calculating the elastic modulus of the mineral mixture using the Voigt-Reuss-Hil average, then adding soft pores and hard pores to the rock matrix of the mineral mixture using the isotropic differential effective medium model, and finally calculating the velocity of the fluid-saturated rock based on the Biot-Rayleigh equation, including the longitudinal wave velocity V P and shear wave velocity V S ; The Biot-Rayleigh equation is a wave propagation equation that describes 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 seismic probabilistic inversion method for reservoir physical property parameters in heterogeneous media according to claim 1, characterized in that: In step S2, the inversion objective function for 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 seismic probabilistic inversion method for reservoir physical property parameters in 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 results of the physical property parameters 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 data y mc Corresponding elastic parameter expansion data x mc , establish a joint data sample {x mc ,y mc }.

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

1. Assuming the joint probability distribution is a two-component Gaussian mixture model, the joint probability distribution function is expressed as follows: Among them, x mc and y 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 elasticity and physical property 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 using 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 parameter expansion data, and The cross-covariance matrix of the first Gaussian component elastic and physical parameter expansion data 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 the 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 expanded data of the elastic parameter of the second Gaussian component, 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 process of physical property parameter inversion.

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

1. Assume that 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 respectively, μ 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 y 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 parameter expansion data, and 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 expanded data of the elastic parameter of the second Gaussian component, Represents the conditional covariance matrix of the second Gaussian component physical parameter expansion data, and Cross-covariance matrix of the first Gaussian component elastic and physical parameter augmented data; S6.

3. After estimating the expectations of different Gaussian components in the posterior probability distribution, the posterior expectation of the physical property parameters is calculated. The 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.

7. The seismic probabilistic inversion method for reservoir physical property parameters in 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.

8. The seismic probabilistic inversion method for reservoir physical property parameters in 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 at the i-th iteration.

Images