Reservoir parameter prediction method based on Bayesian rock physics inversion based on linear Gaussian distribution

Abstract

The present invention relates to a reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution. A new Bayesian rock physics inversion technology, using seismic data as input, combines the Gray reflection coefficient approximation formula with the Gassmann equation and the Nur critical porosity model, and constructs the linear relationship between the bulk modulus, shear modulus and density of fluid-saturated rock and the matrix bulk modulus, matrix shear modulus, matrix density, fluid bulk modulus, fluid density and porosity under the first-order Taylor series approximation. Under the constraints of the Bayesian framework, an iterative algorithm is used to carry out a stable linear Gaussian distribution of rock physics inversion, and then combined with the Voigt-Reuss-Hill average formula, the density-weighted average formula and the Wood formula, to achieve stable prediction of reservoir rock physical parameters such as clay content and water saturation. The method has a low dependence on the initial model, and the posterior mean is assigned to the initial model through continuous iteration until the model converges.

Description

Technical Field

[0001] The invention relates to the technical field of seismic reservoir prediction and seismic rock physics, and in particular to a reservoir parameter prediction technology based on Bayesian rock physics inversion of linear Gaussian distribution. Background Art

[0002] The reservoir parameters of the formation mainly include elastic parameters, fluid types and rock physical properties, etc. Elastic parameters include the longitudinal wave velocity, shear wave velocity, density, bulk modulus, shear modulus, etc. of the formation, fluid types include oil, gas and water in the reservoir pores, etc., and rock physical properties include porosity, permeability, fluid saturation, mineral composition, etc.

[0003] In general, the geophysical inversion problem is to estimate the properties of the subsurface model given a set of geophysical observations, assuming that the physical relationship between the model properties and the observations is known. The most common inversion in reservoir geophysics is seismic inversion, which is widely used in reservoir parameter prediction. Starting from seismic data, the velocity or impedance of the subsurface medium is predicted by seismic amplitude and travel time information, which is called elastic parameter inversion. Similarly, the prediction of rock and fluid properties directly from seismic data, such as rock physical properties such as porosity, mineral composition and fluid saturation, can also be expressed as an inversion, usually called rock physics inversion.

[0004] Seismic wave propagation theory and rock physics models are the bridges that link rock physical properties with seismic observation data. Due to the measurement errors and approximate uncertainties of seismic observation data, the inversion solution is usually not unique. Geophysicists have proposed many methods to solve geophysical inversion problems, including deterministic methods and probabilistic methods. Most of these methods are applied to the inversion solution process of seismic, and then extended to the inversion solution based on rock physical properties. Summary of the invention

[0005] In view of the above problems, the purpose of the present invention is to provide a Bayesian rock physics inversion method based on linear Gaussian distribution. The invention is based on a new linearized forward model, which is based on a linear approximate formula of amplitude changing with offset distance and a linearized rock physics model, and uses a convolution forward model of seismic wavelet and linearized approximate equation to solve the inversion.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] In order to solve this inversion, the seismic forward model must first be linearized. The forward model in the present invention is a convolution model based on seismic wavelets and reflection coefficients. The seismic reflection coefficient is mainly based on the Zoeppritz equation. Many experts and scholars have made many different linear approximations of elastic parameter characterizations, including the Aki-Richards approximation, the Fatti approximation, the Gray approximation, etc. Since the present invention conducts a linearized rock physics inversion, in order to establish a connection with the rock physics model, the reflection coefficient and the bulk modulus K of fluid-saturated rock are used in the present invention. sat , shear modulus μ sat Gray reflectance approximation formula for the linear characteristics between and density ρ.

[0008] Secondly, since linearization of rock physics models is less common in reservoir characterization, the relationship between reservoir elastic parameters and porosity and mineral volume is usually close to linear, but the relationship with water saturation and saturation is usually nonlinear, especially when the pore fluid is uniformly mixed in the pore space. The Gassmann equation gives the elastic modulus of fluid-saturated rock (K sat and μ sat ) and the elastic modulus of dry rock skeleton (K dry and μ dry ), so a model is needed in the Gassmann equation to calculate the elastic modulus of dry rock (K dry and μ dry ), while the Nur critical porosity model is a model that characterizes the elastic modulus of dry rock (K dry and μ dry ) and matrix elastic modulus (K m and μ m ) shows a linear change. In the process of rock physics modeling in the present invention, the Gassmann equation is combined with the Nur critical porosity model to obtain a rock physics model, and the rock physics model is linearized using the first-order expansion approximation of the Taylor series. The linearized rock physics model shows that the bulk modulus K of fluid-saturated rock is sat , shear modulus μ sat and density ρ and bulk modulus K of the matrix m , matrix shear modulus μ m , matrix density ρ m , the bulk modulus K of the fluid f , fluid density ρ f It is linearly related to the porosity φ.

[0009] Again, Gray's approximation gives the reflection coefficient as a function of the bulk modulus K of the fluid-saturated rock: sat , shear modulus K satand density ρ, while the Gassmann equation and Nur critical porosity model give the bulk modulus K of fluid-saturated rock. sat , shear modulus μ sat and density ρ and matrix bulk modulus K m , shear modulus μ m , density ρ m , the bulk modulus K of the fluid f , fluid density ρ f Therefore, the reflection coefficient formula is reparameterized by combining Gray's approximation formula, Gassmann equation and Nur's critical porosity model, and the reflection coefficient is obtained by changing the matrix bulk modulus K. m , shear modulus μ m , density ρ m , the bulk modulus K of the fluid f , fluid density ρ f Expressions for the linear variation of reservoir parameters such as φ and porosity.

[0010] Finally, using the linearized expressions of reservoir parameters and reflection coefficients, and under the constraints of the Bayesian framework, an iterative algorithm is used to carry out rock physics inversion of a stable linear Gaussian distribution to obtain the reservoir parameters (K m , μ m , m , K f , f , φ). The Voigt-Reuss-Hill average formula, density weighted average formula and Wood formula respectively give the reservoir parameters (K m , μ m , m , K f , f , φ) and water saturation s w and clay volume fraction v c The relationship between the rock porosity φ and water saturation s is finally realized by combining the Bayesian rock physics inversion formula. w , clay content v c The reservoir rock physical parameters are predicted.

[0011] The present invention adopts the above technical solution, which has the following advantages:

[0012] 1. The Gassmann equation gives the elastic modulus of fluid-saturated rock (K sat and μ sat ) and the elastic modulus of dry rock skeleton (K dry and μ dry), which itself shows weak nonlinear characteristics. The Nur critical porosity model gives the elastic modulus of dry rock (K dry and μ dry ) and matrix elastic modulus (K m and μ m ) is still a weak nonlinear relationship after combining the two. In order to solve this problem and make the rock physics model have better linear characteristics, the first-order Taylor series expansion is used in the present invention to linearize the weak nonlinear rock physics model and obtain the elastic parameters (K) of the fluid-saturated rock. sat and μ sat ) and reservoir parameters (K m , μ m , m , K f , f , φ).

[0013] 2. The present invention uses seismic data as input. In the Bayesian rock physics inversion process, it has a strong dependence on the rock physics model, and the rock physics model depends on the initial value of the linearized model of the Taylor series expansion. Since the prior value of the reservoir parameter is unknown during the inversion process, the initial value is assumed to be a constant and equal to the mean of the reservoir parameter, but the initial model may often differ greatly from the true value. In order to solve this problem, this patent proposes an iterative method, which has a low dependence on the initial model. The posterior mean is assigned to the initial model through continuous iteration until the model converges. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Various other advantages and benefits will become apparent to those of ordinary skill in the art by reading the detailed description of the preferred embodiments below. The accompanying drawings are only for the purpose of illustrating the preferred embodiments and are not to be considered as limiting the present invention. Throughout the accompanying drawings, the same reference numerals are used to represent the same components. In the accompanying drawings:

[0015] Figure 1 This is a schematic diagram of the logging curve of Well A in a certain work area in the east;

[0016] Figure 2 This is a schematic diagram of the synthetic seismic record using the three angles of the well A data and the wavelet convolution;

[0017] Figure 3 Schematic diagram of the Bayesian rock physics inversion results of the linear Gaussian distribution of Well A;

[0018] Figure 4 This is a schematic diagram of partial superposition of seismic sections at three angles in a certain actual work area in the east;

[0019] Figure 5A schematic diagram of the Bayesian reservoir parameter inversion profile of a linear Gaussian distribution in an actual work area in the east; and

[0020] Figure 6 This is a schematic diagram of the Bayesian rock physical parameter inversion profile of a linear Gaussian distribution in an actual work area in the east.

[0021] The symbols in the accompanying drawings are as follows: DETAILED DESCRIPTION

[0022] The exemplary embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments described herein. On the contrary, these embodiments are provided in order to enable a more thorough understanding of the present invention and to fully convey the scope of the present invention to those skilled in the art.

[0023] A reasonable inversion method should aim to estimate the most likely model while quantifying the uncertainty in the model predictions. Probabilistic approaches to inversion problems provide a natural framework for seismic and rock physics inversion. In probabilistic approaches, the solution to the inversion can be expressed as a probability density function (PDF) of the model variables or a set of model realizations sampled from the PDF, capturing the uncertainty in the model.

[0024] Therefore, the present invention introduces the Bayesian inversion method to evaluate the posterior distribution of model variables given a set of measured geophysical data, and proposes an analytical solution to the seismic inversion problem based on the linearization of the seismic forward model under the Gaussian distribution assumption and the prior distribution of model variables and measurement errors. The Bayesian linearization method was subsequently extended to seismic reservoir prediction and fluid identification. Conventional reservoir prediction and fluid identification can only predict parameters such as reservoir thickness, reservoir fluid type and reservoir scale, and cannot accurately characterize reservoir parameters such as reservoir porosity, reservoir fluid saturation and reservoir mineral volume. Seismic rock physics modeling is to establish the relationship between rock physics and elastic variables through rock physics models.

[0025] Therefore, a Bayesian rock physical reservoir parameter inversion method based on linear Gaussian distribution is developed in the framework of Bayesian distribution of linear Gaussian distribution and combined with seismic rock physics modeling to realize rock physics inversion technology for direct prediction of reservoir parameters.

[0026] The present invention firstly parameters the reflection coefficient model based on the previous research results, and characterizes it into a linear reflection coefficient approximate formula corresponding to the rock physics model. Then, a linear relationship between elastic parameters and differential parameters is constructed through rock physics modeling. Then, the linear reflection coefficient approximate formula is combined with the rock physics model, and the reflection coefficient formula is re-parameterized into a linear representation of reservoir parameters. Finally, Bayesian rock physics inversion prediction of linear Gaussian distribution is carried out under the constraints of the Bayesian framework.

[0027] The entire calculation and implementation process can be divided into four steps:

[0028] Step 1: Derivation of the linearized reflection coefficient approximation formula based on elastic parameters

[0029] Seismic data can be approximately expressed as the convolution of seismic wavelets and a series of reflection coefficients. The magnitude of the reflection coefficient is the difference in elastic parameters of the interfaces between the layers. Therefore, the seismic reflection coefficient depends on the elastic properties of the layers above and below the underground interface. The elastic properties of an isotropic medium can be expressed by three elastic parameters: the longitudinal wave velocity V P , shear wave velocity V S and density ρ. Therefore, seismic data can be characterized as:

[0030] d(t,θ)=w(t,θ)*r PP (t,θ)=∫w(u,θ)r PP (tu,θ)du (1)

[0031] Where d(t,θ) is the synthetic seismic record, w(t,θ) is the seismic wavelet, * is the convolution operator, and r PP (t,θ) is the PP reflection coefficient, t is the two-way travel time of the earthquake, and θ is the incident angle of the earthquake. PP (t,θ) as a function of time t and incident angle θ can be calculated exactly using the Zoeppritz equation or obtained using the Aki-Richards linear approximation.

[0032] In the present invention, the reflection coefficient is approximated by Aki-Richards linear approximation, that is, the intercept R, gradient G and curvature F can be expressed by longitudinal wave velocity, shear wave velocity and density as follows:

[0033]

[0034] in:

[0035]

[0036] Among them, R represents the intercept, G represents the gradient, F represents the curvature, △V P It represents the difference in longitudinal wave velocity between the upper and lower layers. Represents the average value of the longitudinal wave velocity of the upper and lower layers, △V S represents the difference in shear wave velocity between the upper and lower layers of the medium. represents the mean value of the shear wave velocity of the upper and lower layers of the medium, △ρ represents the difference in density between the upper and lower layers of the medium, represents the mean value of the density of the upper and lower layers, t i and t i+1 represents the travel time of two consecutive crossings of the interface. Therefore, the reflection coefficient r PP The approximation of (θ) is:

[0037]

[0038] θ is the seismic incident angle. To simplify the notation, the travel time variable t is omitted. The single interface reflection coefficient in equation (4) can be expanded to a continuous-time reflectivity formula, which can be expressed as:

[0039]

[0040] By expanding equation (5), the reflection coefficient and the saturated rock elastic modulus (K sat and μ sat ) is related to the density ρ and expressed as the AVO approximation based on Gray et al. as follows:

[0041]

[0042] Where K sat , μ sat and ρ represent the bulk modulus, shear modulus and density of saturated rock, respectively.

[0043] Therefore, if the elastic properties of the formation are known, the seismic response can be calculated:

[0044] d(t,θ)=w(t,θ)*r PP (θ) (7)

[0045] The function r PP (θ) can be calculated according to equations (5) and (6).

[0046] Step 2: Construct a rock physics model to obtain the relationship between elastic parameters and reservoir parameters

[0047] Seismic rock physics studies the relationship between rock physical properties (such as porosity, mineral composition and fluid saturation) and elastic properties (such as elastic modulus, velocity, density and impedance). In this patent, a linear rock physics model is mainly constructed based on the Gassmann equation and Nur's critical porosity model. The Gassmann equation is based on the assumption that the underground medium is isotropic and that the pore space is connected and in a state of pressure equilibrium, which is used to estimate the elastic modulus of fluid-saturated rock:

[0048]

[0049] Where K sat is the bulk modulus of fluid-saturated rock, K dry is the bulk modulus of dry rock skeleton, K m is the rock matrix bulk modulus, K f is the bulk modulus of the fluid, φ is the porosity, μ sat is the shear modulus of fluid-saturated rock, μ dry is the shear modulus of the dry rock skeleton. Since the presence of fluid has no effect on the shear modulus of rock, the shear modulus of fluid-saturated rock is μ sat The shear modulus μ of the dry rock skeleton dry equal.

[0050] The elastic modulus of the dry rock skeleton is either measured in the laboratory or calculated by rock physics models. Nur (1992) proposed a critical porosity model to estimate the relationship between the dry rock skeleton modulus and the matrix modulus, where the dry rock elastic modulus is expressed as a linear function of the porosity φ, and the linear regression coefficient depends on the matrix modulus and the critical porosity φ0, and its expression is:

[0051]

[0052] The rock density ρ is usually expressed as a linear combination of matrix density and fluid density, expressed as:

[0053] ρ=ρ m (1-φ)+ρ f φ (10)

[0054] where ρ m is the density of the rock matrix, ρ f is the density of the pore fluid.

[0055] By combining the Gassmann equation of equation (8), the critical porosity model of equation (9) and the density model of equation (10), a rock physics model is obtained. This rock physics model converts the elastic property (K sat , μ sat and ρ) and reservoir properties (Km , μ m , m , K f , f and φ), where the bulk modulus K of the matrix m , shear modulus μ m and density ρ m Depending on the mineral composition of the porous rock, the bulk modulus K of the pore fluid is f and density ρ f Depends on the composition of the pore fluid.

[0056] For fluid-saturated rock shear modulus μ sat and density ρ, Equations (8) and (9) show that μ sat Relative to μ m is linear, and equation (10) shows that ρ is linear with respect to ρ m is linear. Therefore, the elastic model of saturated rock varies linearly with the matrix elastic model. The main nonlinear changes in this rock physics model are usually due to the fluid bulk modulus K f Caused by, where the fluid volume model K f The expression is:

[0057]

[0058] Where K w represents the bulk modulus of water in the pores, K g represents the bulk modulus of the gas in the pores, s w Indicates water saturation. Equations (8)-(11) show that in the presence of multiple pore fluids, the elastic modulus of saturated rock (K sat and μ sat ) and water saturation w and degree show a nonlinear relationship, but are related to the fluid bulk modulus K f It is a linear relationship.

[0059] For fluid-saturated rock, the bulk modulus K sat Since the Gassmann equation exhibits weak nonlinear characteristics, a mathematical method for linearizing weak nonlinear rock physics models using a first-order Taylor series approximation is proposed in the present invention. The purpose of the Taylor series is to approximate an arbitrary function using a polynomial under given values ​​of the independent variable, where the coefficients depend on the derivative of the function. For simple functions, a good approximation can be obtained using a finite number of terms. The Taylor series is used in geophysics to approximate complex nonlinear functions. In rock physics inversion, we are interested in the linear approximation of the rock physics model. Therefore, the Taylor series expansion of the rock physics equation only retains the first term.

[0060] Therefore, it is necessary to linearize the reservoir parameters (matrix bulk modulus K m , matrix shear modulus μ m , matrix density ρ m , fluid bulk modulus K f , fluid density ρ f These reservoir parameters are used as input to calculate the fluid-saturated rock bulk modulus K sat , saturated rock shear modulus μ sat and elastic parameters such as rock density ρ.

[0061] Step 3: Model parameterization reflection coefficient to obtain the linear relationship between reservoir parameters and reflection coefficient

[0062] The main goal of seismic reservoir characterization is to predict models of rock and fluid properties based on well logging and seismic data. Formation elastic parameters depend on rock and fluid properties. When porosity increases, velocity generally decreases; when gas saturation increases and the corresponding water saturation decreases, velocity generally decreases. If d represents the observed seismic data and q represents the elastic parameters that need to be characterized in the real underground, we can relate the elastic parameters q to the observed seismic data d through a set of geophysical operators f:

[0063] d=f(q)+ε (12)

[0064] Where ε represents random noise in the data. The geophysical operator f can include rock physics models and seismic models. The rock physics model includes equivalent medium theory, inclusion model, and empirical model that links reservoir parameters with elastic parameters. The seismic model mainly refers to the geophysical equation that links elastic parameters with seismic response, such as Zoeppritz equation and Aki-Richards linear approximation equation. The elastic parameter q is solved as an inversion.

[0065] Based on the present invention, a rock physics inversion is carried out, q represents elastic parameters, such as P-wave velocity, S-wave velocity and density, m represents reservoir parameters, such as porosity φ, clay content v c and water saturation s w Etc. Therefore, the rock physics model is expressed by function g:

[0066] q=g(m)+e (13)

[0067] Where e represents the reservoir parameter data error. Substituting equation (13) into equation (12), the rock physics inversion can be written as follows:

[0068] d=f(q)+ε=f[g(m)]+ε (14)

[0069] In step 2, the relationship between the elastic parameter q and the reservoir parameter m is constructed through rock physics modeling. Let the reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T , elastic parameter q=[K sat ,μ sat ,ρ] T .

[0070] After taking the logarithm of both sides of formula (13) and expanding it by Taylor series, we can get the linearized expression as follows:

[0071]

[0072] in is the average value of the reservoir parameter, is a constant at m0, is the Jacobian determinant of the rock physics model, and its value at m0 can be written as follows:

[0073]

[0074] Substituting equations (15) and (16) into equation (6), we can obtain the reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The linearized expression of the reflection coefficient is as follows:

[0075]

[0076] in

[0077]

[0078] in

[0079]

[0080] in represents the mean value of the matrix bulk modulus, represents the mean value of matrix shear modulus, represents the mean value of matrix density, represents the mean value of the bulk modulus of the fluid, represents the mean value of the fluid density, represents the mean porosity, and φ0 represents the critical porosity value.

[0081] Step 4: Rock physics inversion based on linear Gaussian distribution in Bayesian framework

[0082] In the Bayesian framework, the posterior distribution P(m|d) of the reservoir parameter m under the condition of the observed seismic data d can be expressed as:

[0083]

[0084] Where P(m) is the prior distribution of reservoir parameters, P(d|m) is the likelihood function of observed seismic data, and P(d) is the probability of observing seismic data, also known as the normalization constant. The formula in equation (20) can be applied to discrete random variables and continuous random variables. In geophysical inversion applications, the inversion problem of continuous random variables includes the inversion of elastic parameters and rock physics of seismic data. The variables of elastic parameter inversion mainly include P-wave velocity, S-wave velocity and density, and the variables of rock physics inversion mainly include porosity, mineral volume and fluid saturation. Discrete random variable inversion mainly includes lithology or rock type classification of seismic data.

[0085] In the case where the reservoir parameter variable m is continuous, the most common PDF for a continuous variable is the Gaussian distribution, often called the normal distribution. We say that the random variable X follows a Gaussian distribution. distribution, if its PDFf X (x) can be written as:

[0086]

[0087] Gaussian distribution The mean is 0 and the variance is 1, which is also called the standard Gaussian distribution (or normal distribution). The Gaussian distribution is symmetrical and unimodal, and can be used to describe many phenomena in nature. The Gaussian distribution is defined over the entire set R of real numbers, and it is always positive.

[0088] In order to solve the problem of reservoir parameter m inversion, the present invention proposes a linear Gaussian Baye inversion method. Assume that the reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The prior distribution of m and the prior covariance matrix is ​​∑ m Gaussian distribution, m:N(m;ξ m ,∑ m ), the error term ε has a mean of 0 and a covariance matrix of ∑ ε Gaussian distribution, ε:N(ε; 0,∑ ε ), and is independent of the reservoir parameter m; the posterior distribution m|d also obeys the conditional mean ξm|d Gaussian distribution, m|d:N(m;ξ m|d ,∑ m|d ), the posterior mean ξ m|d It can be expressed as follows:

[0089] ξ m|d =ξ m +∑ m G T (G∑ m G T +∑ ε ) -1 (d-Gξ m ) (twenty two)

[0090] Conditional covariance matrix∑ m|d as follows:

[0091] ∑ m|d =∑ m -∑ m G T (G∑ m G T +∑ ε ) -1 G∑ m (twenty three)

[0092] Where G represents the matrix associated with the linear rock physics model, which depends on the value of the linearized initial model m0 of the Taylor series expansion. In equations (17) and (18), we assume that the initial model is a constant and equal to the average value of the reservoir parameters, because the prior values ​​of the reservoir parameters are unknown. If the initial value is significantly different from the actual inverted reservoir parameter value, it means that the rock physics model has a strong nonlinear characteristic and the linearization may be inaccurate. In order to solve this problem, this patent proposes an iterative method, assuming that the initial value of the model is Calculate the solution of inversion (22) and (23), and then replace m0 = ξ m|d , and iterate until convergence.

[0093] The main steps of the iterative algorithm are as follows:

[0094] Step 1: Set the initial model parameters

[0095]

[0096] Step 4: If yes, end the calculation; if no, let m0 = ξ m|d ;

[0097] Step 5: Update the matrix G related to the linear rock physics model and jump to step 2;

[0098] The rock reservoir parameter m is calculated by Bayesian inversion algorithm: m ,μ m ,ρ m ,K f ,ρ f ,φ] T , combined with rock physics models such as the Voigt-Reuss-Hill average formula, weighted average formula, and Wood formula, other rock physics parameters such as mineral volume fraction and water saturation are calculated. Assuming that the rock components contain two minerals, quartz and clay, and two fluids, gas and water, the Voigt-Reuss-Hill average formula, density weighted average formula, and Wood formula are characterized as follows:

[0099]

[0100]

[0101] ρ m =v c ρ c +(1-v c )ρ q (26)

[0102] ρ f =s w ρ w +(1-s w )ρ g (27)

[0103]

[0104] where v c Indicates clay content, K c represents the bulk modulus of clay, K q represents the bulk modulus of quartz, μ c represents the shear modulus of clay, μ q represents the shear modulus of quartz, ρ c represents the density of clay, ρ q Indicates the density of quartz, s w represents water saturation, ρ w represents the density of water, ρ g Indicates the density of gas, K w represents the bulk modulus of water, K g represents the bulk modulus of gas.

[0105] The modulus and density (K) of all minerals and fluids c , K q , K w , K g , μ c, μ q , c , q , w , g ) are constants and known, then according to equations (24) to (28), the clay content v c K in the main formulas (24)-(26) m , μ m and ρ m To determine the water saturation s w Mainly from formula (27)-(28) f and K f According to the Bayesian linear rock physics inversion formulas (22) and (23), the clay content v c and water saturation s w The posterior distribution of can have reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The calculation is performed as follows:

[0106]

[0107]

[0108] in Indicates clay content v c The cumulative density function of Indicates water saturation s w The cumulative density function of and Marginal probability density functions of rock and fluid obtained from equations (22) and (23).

[0109] The present invention also provides a practical application embodiment, which is as follows:

[0110] It is known that the logging curve of a certain area in the east is as follows Figure 1 As shown, there are six curves including P-wave velocity, S-wave velocity, density, porosity, clay content and water saturation. Assume that there are three angles partially superimposed, namely small, medium and large angles, corresponding to incident angles of 10, 20 and 30 degrees respectively. Then use 30Hz Ricker wavelet to calculate the corresponding seismic response according to formula (1). First, use equation (5) to calculate the linearized approximate reflection coefficient, and then use the corresponding convolution model in formula (1) to synthesize the corresponding seismic records for different incident angles. The final seismic data set is as follows: Figure 2 As shown, the vertical axis is the round-trip travel time.

[0111] right Figure 2 The Bayesian rock physics inversion technique based on linear Gaussian distribution is carried out on the partial stacked synthetic seismic records in the image, and the corresponding rock physical parameters can be obtained, such as Figure 3 As shown. The black is the measured value of the well logging, the gray is the Bayesian rock physics inversion result of the linear Gaussian distribution of Well A predicted by the patented technology of the present invention, and the filling between the gray and black represents the probability of the inversion result falling in this interval. It can be seen from the figure that the inversion result is in good agreement with the measured value of the well logging, which proves the feasibility of the inversion technology of the present invention from a single well.

[0112] In order to further verify the reliability of the inversion technology in the present invention, the technology was applied to an actual work area of ​​an oil field in the east. Figure 4 Partially superimposed seismic profiles at three angles of an actual work area in the east are given. It can be seen from the figure that there are obvious differences in the amplitude of the seismic data at different incident angles. The present invention utilizes the difference in amplitude variation with incident angle, combines it with the Bayesian framework, and carries out a Bayesian rock physics inversion of linear Gaussian distribution when the prior information satisfies the Gaussian distribution. Figure 5 The Bayesian reservoir parameter (K) of the linear Gaussian distribution is given using formula (17) in the Bayesian framework. m , μ m , m , K f , f , φ) inversion profile. It can be seen from the figure that the matrix bulk modulus K m , matrix shear modulus μ m , matrix density ρ m , fluid bulk modulus K f , fluid density ρ f The inversion profiles of porosity and porosity φ are in good agreement with the measured values ​​of well logging. Figure 5 The inversion results of the Bayesian inversion formula are combined with the Voigt-Reuss-Hill average formula, the density weighted average formula, the Wood formula and the Bayesian inversion formula to finally obtain the Bayesian rock physical parameter inversion results based on the linear Gaussian distribution as follows: Figure 6 As shown in the figure, it can be seen that the porosity φ and clay content v c and water saturation s w All of them are in good agreement with the measured values ​​on the well, further proving that the prediction method in the present invention is a reliable inversion prediction technology.

[0113] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A reservoir parameter prediction method based on Bayesian rock physics inversion with linear Gaussian distribution, characterized in that: include: Gray's approximate formula is used to give the linear relationship between the reflection coefficient and the bulk modulus, shear modulus and density of fluid-saturated rock. The relationship between the elastic modulus of fluid-saturated rock and the elastic modulus of dry rock skeleton is characterized by the Gassmann equation, and the linear relationship between the elastic modulus of dry rock skeleton and the elastic modulus of matrix is ​​characterized by the Nur critical porosity model. The rock physics model is obtained by combining the Gassmann equation with the Nur critical porosity model. The rock physics model is linearized using a first-order expansion approximation of a Taylor series, and the linearized rock physics model indicates that the bulk modulus of the fluid-saturated rock, the shear modulus of the fluid-saturated rock, and the density of the fluid-saturated rock are linearly related to the bulk modulus of the matrix, the shear modulus of the matrix, the density of the matrix, the bulk modulus of the fluid, the density of the fluid, and the porosity of the rock; The Gray approximation formula, Gassmann equation and Nur critical porosity model are combined to parameterize the reflection coefficient model, and an expression is obtained that the reflection coefficient changes linearly with the bulk modulus of the matrix, the shear modulus of the matrix, the density of the matrix, the bulk modulus of the fluid, the density of the fluid and the porosity of the fluid. and Under the constraints of the Bayesian framework, an iterative algorithm is used to carry out the inversion of the rock physics model of the stable linear Gaussian distribution, and then the Voigt-Reuss-Hill average formula, the density weighted average formula and the Wood formula are combined to realize the prediction of reservoir parameters. The reflection coefficient uses the Aki-Richards linear approximation, and the intercept, gradient, and curvature are expressed in terms of longitudinal wave velocity, shear wave velocity, and density as follows: in: Among them, R represents the intercept, G represents the gradient, F represents the curvature, ΔV P It represents the difference in longitudinal wave velocity between the upper and lower layers. Represents the average value of the longitudinal wave velocity of the upper and lower layers, ΔV S represents the difference in shear wave velocity between the upper and lower layers of the medium. represents the mean value of the shear wave velocity of the upper and lower layers of the medium, Δρ represents the difference in density between the upper and lower layers of the medium, represents the mean value of the density of the upper and lower layers, t i and t i+1 They represent two consecutive travel times across the interface respectively; Reflection coefficient r PP The approximation of (θ) is: Among them, r PP (t,θ) represents the PP reflection coefficient, t represents the earthquake two-way travel time, and θ represents the earthquake incident angle; The single interface reflection coefficient in equation (4) is expanded into the continuous-time reflectivity formula: By expanding equation (5) and linking the reflection coefficient with the elastic modulus and density of saturated rock, it can be expressed as follows: Where K sa t represents the bulk modulus of fluid-saturated rock, μ sat represents the shear modulus of fluid-saturated rock, ρ represents the density of fluid-saturated rock; If the elastic properties of the formation are known, the seismic response can be calculated: d(t,θ)=w(t,θ)*r PP (i) (7) Among them, the function r PP (θ) is calculated according to equations (5) and (6); The iterative algorithm steps are as follows: Step 1: Let the initial model parameters Step 2: Calculate the posterior mean ξ md and the conditional covariance matrix Σ md ; Step 3: Determine ||dG(ξ m|d )|| is it less than the set error? Step 4: If yes, end the calculation; if no, let m0 = ξ md ; Step 5: Update the matrix G related to the linear rock physics model and jump to step 2; The rock reservoir parameter m is calculated by Bayesian inversion algorithm. m ,μ m ,ρ m ,K f ,ρ f ,φ] T , assuming that the rock components contain two minerals, quartz and clay, and two fluids, gas and water, the Voigt-Reuss-Hill average formula, the density weighted average formula, and the Wood formula are respectively expressed as follows: r m =v c r c +(1-v c )r q (26) r f =s w r w +(1-s w )r g (27) where v c Indicates clay content, K c represents the bulk modulus of clay, K q represents the bulk modulus of quartz, μ c represents the shear modulus of clay, μ q represents the shear modulus of quartz, ρ c represents the density of clay, ρ q Indicates the density of quartz, s w represents water saturation, ρ w represents the density of water, ρ g Indicates the density of gas, K w represents the bulk modulus of water, K g represents the bulk modulus of gas; The modulus and density (K) of all minerals and fluids c , K q , K w , K g , μ c , μ q , c , q , w , g ) are constants and known, then according to equations (24) to (28), the clay content v c K in the main formulas (24)-(26) m , μ m and ρ m To determine the water saturation s w Mainly from formula (27)-(28) f and K f to decide; According to the Bayesian linear rock physics inversion formulas (22) and (23), the clay content v c and water saturation s w The posterior distribution of is given by the reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The calculation is performed as follows: in Indicates clay content v c The cumulative density function of Indicates water saturation s w The cumulative density function of and Marginal probability density functions of rock and fluid obtained from equations (22) and (23).

2. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 1 is characterized in that: The Gassmann equation assumes that the subsurface medium is isotropic and that the pore space is connected and in pressure equilibrium, and is used to estimate the elastic modulus of fluid-saturated rock: Where K sat is the bulk modulus of fluid-saturated rock, K dry is the bulk modulus of dry rock skeleton, K m is the rock matrix bulk modulus, K f is the bulk modulus of the fluid, φ is the porosity, μ sat is the shear modulus of fluid-saturated rock, μ dry is the shear modulus of dry rock skeleton, and the shear modulus of fluid-saturated rock is μ sat The shear modulus μ of the dry rock skeleton dry equal.

3. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 2 is characterized in that: The Nur critical porosity model is used to estimate the relationship between the dry rock skeleton modulus and the matrix modulus, where the dry rock elastic modulus is expressed as a linear function of porosity, and the linear regression coefficient depends on the matrix modulus and the critical porosity, and its expression is: Where φ0 represents the critical porosity; Rock density is usually expressed as a linear combination of matrix density and fluid density, expressed as: p=p m (1-φ)+ρ f f (10) where ρ m is the density of the rock matrix, ρ f is the density of the pore fluid; By combining equations (8), (9) and (10), the rock physics model is obtained.

4. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 3 is characterized in that: The rock physics model is obtained by: For the shear modulus and density of fluid-saturated rock, equations (8) and (9) show that the shear modulus of fluid-saturated rock is linear with respect to the matrix shear modulus, and equation (10) shows that the rock density is linear with respect to the matrix density. The elastic model of saturated rock changes linearly with the matrix elastic model. The nonlinear changes in the rock physics model are caused by the fluid bulk modulus, where the expression of the fluid bulk model is: in: K f represents the bulk modulus of the fluid, K w represents the bulk modulus of water in the pores, K g represents the bulk modulus of the gas in the pores, s w Indicates water saturation.

5. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 1 is characterized in that: Equation (12) is used to relate the elastic parameters to the observed seismic data: d=f(q)+ε (12) in, d represents earthquake data; f represents a geophysical operator; q represents the elasticity parameter; ε represents the random noise in the data; The rock physics inversion includes: The rock physics model is expressed by the function g: q=g(m)+e (13) Among them, e represents the reservoir parameter data error, q represents the elastic parameter, and m represents the reservoir parameter; Substituting equation (13) into equation (12), we get the rock physics inversion: d=f(q)+ε=f[g(m)]+ε (14) The relationship between the elastic parameter q and the reservoir parameter m is constructed through rock physics modeling. Let reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T , elastic parameter q=[K sat ,μ sat ,ρ] T ; After taking the logarithm of both sides of formula (13) and expanding it by Taylor series, we can get the linearized expression as follows: in is the average value of the reservoir parameter, is a constant at m0, is the Jacobian of the rock physics model; The value at m0 is written as follows: Substituting equations (15) and (16) into equation (6), we obtain Reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The linearized expression of the reflection coefficient is as follows: in Where K m represents the mean value of the matrix bulk modulus, μ m represents the mean value of matrix shear modulus, ρ m represents the mean value of matrix density, K f represents the mean value of the bulk modulus of the fluid, ρ f represents the mean value of the fluid density, represents the mean porosity, and φ0 represents the critical porosity value.

6. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 1, characterized in that: In the Bayesian framework, the posterior distribution P(m|d) of the reservoir parameter m under the condition of the observed seismic data d is expressed as: Among them, P(m) is the prior distribution of reservoir parameters, P(d|m) is the likelihood function of observed seismic data, and P(d) is the probability of observed seismic data; The formula in (20) applies to both discrete random variables and continuous random variables.

7. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 6, characterized in that: In the case where the reservoir parameter variable m is continuous, the most common PDF for a continuous variable is the Gaussian distribution, often called the normal distribution; The random variable X is distributed according to Gaussian distribution. distribution, if its PDF f X (x) is written as: Gaussian distribution A distribution with a mean of 0 and a variance of 1 is called a standard Gaussian distribution.

8. The reservoir parameter prediction method based on Bayesian rock physics inversion of linear Gaussian distribution according to claim 7, characterized in that: The inversion of reservoir parameter m is solved by using the linear Gaussian Baye inversion method. Assuming that reservoir parameter m = [K m ,μ m ,ρ m ,K f ,ρ f ,φ] T The prior distribution of m and the prior covariance matrix is ​​∑ m Gaussian distribution, m N(m;ξ m ,∑ m ), the error term ε has a mean of 0 and a covariance matrix of ∑ ε Gaussian distribution, εN(ε; 0,∑ ε ), and is independent of the reservoir parameter m; the posterior distribution m|d also obeys the conditional mean ξ m|d Gaussian distribution, m|d N(m;ξ m|d ,∑ m|d ), the posterior mean ξ m|d It is expressed as follows: x m|d =ξ m +∑ m G T (GΣ m G T +S ε ) -1 (d-Gξ m ) (22) Conditional covariance matrix Σ m|d as follows: S m|d =∑ m -S m G T (GΣ m G T +∑ ε ) -1 G∑ m (23) Where G represents the matrix associated with the linear rock physics model, which depends on the value of the linearized initial model m0 of the Taylor series expansion.

Images