CN106772604B - Prestack seismic inversion method based on the fluid volume compressed coefficient - Google Patents

Abstract

The invention discloses the prestack seismic inversion methods based on the fluid volume compressed coefficient, comprising: and it is theoretical based on pore media, using theoretical and experience petrophysical model, establish contacting for the fluid volume compressed coefficient and other elastic parameters and physical parameter；On the basis of this relationship, the reflectance signature equation based on the fluid volume compressed coefficient is derived；Using contacting between reflectance signature establishing equation seismic data and the fluid volume compressed coefficient, and then realize the inverting to fluid coefficient of bulk compressibility.Air water can be separated well on the fluid volume compressed coefficient inverting section obtained by means of the present invention, improve the precision of fluid identification.

Description

Pre-stack seismic inversion method based on fluid volume compression coefficient

Technical Field

The invention relates to a prestack seismic data inversion method in the process of predicting an oil and gas reservoir, in particular to a prestack seismic inversion method based on a fluid volume compression coefficient.

Background

Reservoir fluid identification is an important link in reservoir exploration and reservoir evaluation. Today, where human oil and gas exploration activities have gone through two or more centuries, we are faced with increasingly complex subterranean conditions, highly cryptic residual oil and gas resource distribution, identification and description of difficult lithologic formation reservoirs, and high exploration and development costs. Meanwhile, with the continuous progress and development of seismic acquisition, processing and interpretation technologies, the identification of the fluid contained in the reservoir based on the information such as amplitude, frequency and the like contained in seismic data becomes possible under the guidance of geological knowledge and oil-gas enrichment rules. The reservoir fluid identification with the seismic data as the main body can increase the exploration success rate to a certain extent, improve the production efficiency and reduce the exploration and development cost, and is one of the research hotspots in the field of oil and gas exploration at the present stage.

In order to reduce the risk of exploration and development, improving the exploration accuracy becomes a problem of intensive research in recent years by geophysical workers, and therefore, reservoir fluid identification technology based on seismic data becomes a hot point of research. Smith and Gidlow (1987) proposed that lithology and fluids could be predicted by stacking different weighting functions using prestack data to obtain fluid factors and pseudo-Poisson's ratio profiles, and first proposed the concept of fluid factors that prompted the development of AVO intercept-gradient intersection techniques for fluid anomaly identification (Verm and Hilterman, 1995; Castagna, 1998). Goodway et al (1997) proposed a lambda-mu-rho technique for reservoir fluid identification using Lamei elasticity parameters. Hilterman (2001) summarized AVO-based fluid identification techniques based on the results of Goodway and Hedlin et al. Batzle (2001) compared the fluid indicator factor, proposed that the Lame parameter combination property is most sensitive to the fluid type for clastic rock, and particularly emphasizes the selection of the fluid factor sensitivity according to the regional characteristics in practical application. George (2003) provides concepts of fluid factor angles and intersection angle according to prestack AVO analysis, and finds that the two attributes have strong identification capability on reservoir fluid types through model trial calculation and practical application. Ningzhonghua et al (2006) proposed the concept of a high sensitivity fluid factor based on a summary of the pre-analysis methods. Mark et al (2006) propose the concept of Poisson impedance. Russell et al (2003, 2006) summarize the previous view taking into account porous saturated elastic media, and use the Biot-Gassmann equation to rewrite the compressional velocity equation under saturated fluid conditions, propose it as a fluid indicator, and point out that the fluid term can directly participate in fluid detection and reservoir prediction as a fluid factor. Hilterman (2009) combined with example application focuses on studying the fluid factor sensitivity of unconsolidated sandstone and indicates that the key to determining the fluid identification sensitivity of unconsolidated sandstone reservoir is the selection of the longitudinal and transverse wave velocity.

The calculation mode of the fluid factor at the present stage is based on indirect combination operation of elastic parameters, and the quality of fluid identification mainly depends on two aspects, namely whether the inversion of the basic elastic parameters is reliable or not; the second is whether the constructed fluid factor is sensitive to the type of the pore fluid. The current common fluid factor has the defect of poor inversion accuracy in the inversion process. Considering that pre-stack seismic inversion is a main means for extracting basic elastic parameters, the reliability of the elastic parameters can be better improved by improving an inversion method; on the other hand, the fluid factor based on indirect combination inevitably causes accumulated errors, the direct inversion of the fluid factor can better solve the problem, the sensitive fluid factor is constructed based on the double-phase medium rock physics theory, the direct inversion of the fluid factor can be realized by researching the internal relation between the sensitive fluid factor and the seismic motion reflection characteristic and utilizing pre-stack seismic data, and therefore the indication sensitivity of the fluid factor and the reliability of the estimation of the fluid factor are improved.

Disclosure of Invention

The invention aims to provide a prestack seismic inversion method based on fluid volume compression coefficients aiming at the defects of the prior art.

The technical scheme adopted by the invention is as follows.

The prestack seismic inversion method based on the fluid volume compression coefficient comprises the following steps: establishing a relation between the volume compression coefficient of the fluid and other elastic parameters and physical parameters by using a theory and an empirical rock physical model based on a pore medium theory; on the basis of the relation, a reflection characteristic equation based on the volume compressibility of the fluid is deduced; and establishing a relation between the seismic data and the fluid volume compression coefficient by using a reflection characteristic equation so as to further realize inversion of the fluid volume compression coefficient.

Further, the prestack seismic inversion method based on the fluid volume compression coefficient is characterized by comprising the following steps:

step 1: establishing a relation between the volume compression coefficient of the fluid and other elastic parameters and physical parameters by using a theory and an empirical rock physical model based on a pore medium theory; other elastic parameters comprise Gassmann fluid factor, shear modulus and density, and the physical parameters comprise water saturation, porosity and argillaceous content;

step 2: under the assumption of plane waves, according to the relationship between the fluid volume compression coefficient obtained in the step 1 and other elastic parameters or physical parameters, an AVO (amplitude of seismic waves changes along with offset) seismic reflection characteristic approximation equation and an elastic impedance equation based on the fluid volume compression coefficient are deduced;

and 3, establishing a seismic data-Bayes elastic impedance inversion-fluid volume compression coefficient inversion process by using the elastic impedance equation obtained in the step 2, and realizing direct inversion of the fluid volume compression coefficient under a Bayes theory framework to form a relatively complete and stable pre-stack seismic inversion method based on the fluid volume compression coefficient.

Further, after the inversion is completed, the method also comprises the step of identifying the reservoir fluid by using the fluid volume compression coefficient obtained by the inversion.

And further, during reservoir fluid identification, the method also comprises the steps of carrying out sensitivity analysis on the fluid volume compressibility in the step 1 based on rock physics theory and well logging information, and verifying the fluid identification precision of the fluid volume compressibility.

Further, the method of the sensitivity analysis is as follows:

combining logging information, and comparing the fluid volume compression coefficient with the fluid factor type based on the conventional single-phase medium theory to obtain fluid sensitivity; defining a certain fluid factor sensitivity as:

wherein mean is_gas,mean_waterThe average values, std, of the gas layer and the water layer corresponding to the fluid factor of the target interval on the well are_gasThen it is the standard deviation.

Further, the fluid factor is one or a combination of several of longitudinal and transverse wave velocity, longitudinal and transverse wave impedance, poisson ratio, lame constant density, shear modulus density, Gassmann fluid factor and poisson ratio.

Further, in step 1, the pore media theory is a two-phase media theory.

Relative reduction rate of volume per unit pressure at a certain temperature

And the compression coefficient C of the two-phase medium is phi C_f+(1-φ)C_sIn which C is_f＝S_oC_o+S_wC_w+S_gC_g，C_g,C_oAnd C_wRespectively representing the fluid volume compression coefficients of gas, oil and water; s_g、S_oAnd S_wRespectively represent the gas, oil and water saturations, and S_g+S_w+S_o＝1。

The Batzle and Han find that the main factors influencing the value of Gassmann fluid item f are the volume modulus of the pore fluid and the porosity of the solid skeleton; dehua Han and Batzle et al studied the coupling effect of the solid effect (porosity, mineral modulus, etc.) of the pore fluid and the rock skeleton in the Biot-Gassmann theory on the rock modulus information by performing rock physics statistics on clastic rock, and proposed that the factor C for removing the rock skeleton fluid is highlighted_fThe empirical formula for the Gassmann fluid term of (1):

f＝G(φ)/C_f (2)。

wherein,wherein the gain function G (phi) represents the combined effect of the rock skeleton minerals and the porosity.

By combining an empirical relation obtained by rock physics experiments, the fluid volume compressibility coefficient C can be obtained by utilizing logging information_fCarrying out estimation; therefore, the fluid volume compressibility C is extracted from the seismic data by a geophysical method_fCompressing the volume of the fluid by a factor C_fThe solid skeleton and the fluid elastic effect can be decoupled by taking the fluid factor as a fluid factor to participate in fluid identification, so that the reliability of reservoir pore fluid identification is effectively improved.

Further, in step 2, an AVO seismic reflection characteristic approximation equation and an elastic impedance equation based on the fluid volume compressibility are derived by using the relationship between the fluid volume compressibility and the Gassmann fluid term.

The approximate equation for the reflection coefficient containing the Gassmann fluid term is shown below:

wherein f, mu and rho respectively represent the average values of Gassmann fluid terms, shear modulus and density of media on two sides of the interface; Δ f, Δ μ, and Δ ρ represent the difference in shear modulus and density, respectively, of the Gassmann fluid terms on either side of the interface.

Substituting the formula (2) into the Russell approximation formula (4), and performing corresponding transformation to obtain:

the dry rock shear modulus μ is used herein in view of its immunity to the pore fluid_dryInstead of μ, equation (5) can be further simplified to:

nur has shown through extensive research that for rocks less than the critical porosity, it isThe bulk modulus and shear modulus of the dry rock can be adjusted to the critical porosity phi_cThe relevant linear function is expressed, namely a critical porosity model, and the expression is shown as the following formula:

wherein phi is_cRepresents the critical porosity, K_dryDenotes the bulk modulus, μ, of the dry rock_dryDenotes the shear modulus, K, of the dry rock_mDenotes the bulk modulus, μ, of the solid mineral matrix_mRepresenting the shear modulus of the mineral matrix; further expansion of equation (6) with the Nur model as the tie can result in:

will be provided withSubstituting G (φ), further expansion can result in:

suppose thatThen there is F_poro＝φ_cPhi mu; the formula can be further simplified as:

if f is set_mAnd phi mu, the fluid volume compressibility approximation formula is finally obtained as follows:

wherein f is_mPhi mu, which is called the solid stiffness parameter;

by using the idea of Connolly to derive elastic impedance, the reflection coefficient is expressed by the elastic impedance, and the following is obtained:

substituting equation (12) into (11) yields:

expressing the relative change of the elastic parameter in logarithmic form, we can get:

further variations, one can obtain:

namely, it is

And (3) integrating and indexing two sides of the above equation, eliminating differential terms and logarithmic terms on two sides of the equation, and further removing an integral constant of 0 to obtain:

EI(θ)＝(C_f)^a(θ)(f_m)^b(θ)(ρ)^c(θ)(φ)^d(θ) (17)。

wherein,

similar to the conventional elastic impedance formula, formula (17) also has a problem that the numerical dimension changes with angle; four reference constants, namely A, are introduced here₀，C_f0，f_m0，ρ₀And phi₀(ii) a The formula (17) is subjected to standardization processing, so that a standardized elastic impedance equation based on the fluid factor of the rock-removing framework can be obtained;

EI(θ)＝A₀(C_f/C_f0)^a(θ)(f_m/f_m0)^b(θ)(ρ/ρ₀)^c(θ)(φ/φ₀)^d(θ) (18)。

C_f0，f_m0，ρ₀and phi₀Are respectively defined as C_f0，f_mρ and φ; a. the₀For the normalization factor, the specific expression is:

further, in step 3, a Bayesian inversion framework is utilized, a likelihood function reflecting seismic information is combined with a prior geological constraint reflecting a parameter to be inverted, an inversion target function is established in a mode of solving a maximum posterior probability density function, conventional Bayesian AVO inversion is conducted on an algorithm equation, correlation characteristics among parameters are not considered in the prior constraint, four-variable Cauchy distribution is adopted as a prior regular constraint to improve four-parameter inversion of Bayesian AVO on the basis of the decorrelation processing of the algorithm equation in order to further improve inversion quality, and model trial calculation and practical application show that the reliability of the four-parameter inversion is improved to a certain extent.

The formula (11) is expressed as a matrix form according to the difference of the incident angles:

wherein, a_i(i＝1,2…,m)，b_i(i＝1,2…,m)，c_i(i-1, 2 …, m) and d_i(i ═ 1,2 …, m) respectively represent the corresponding coefficients for the ith angle of incidence; the method is popularized to the condition that m incident angles and n interfaces exist, and the matrix is subjected to blocking processing to obtain the following results:

wherein R is_i(i ═ 1,2 …, m) represents the reflection coefficient vector for the ith angle of incidence, consisting of n elements; a. the_i(i＝1,2…,m)，B_i(i＝1,2…,m)，C_i(i-1, 2 …, m) and D_i(i ═ 1,2 …, m) respectively represent forward modeling coefficient matrices corresponding to the i-th incident angle, and are respectively diagonal matrices of n × n dimensions;R_ρand R_φRespectively representing the fluid factor of the rock-removing framework, the solid rigidity parameter, the density and the relative change rate vector of the porosity, and respectively consisting of n elements.

With the wavelet matrix W introduced based on the assumption that the seismic records conform to the convolution model, equation (21) becomes further:

wherein d is_i(i ═ 1,2 …, m) column vectors consisting of seismic data for the ith angle of incidence, all containing n elements;

generating a covariance matrix according to a sample statistical method of logging data, namely:

wherein,n represents the number of samples;

covariance matrix C of parameters to be inverted after decorrelation transformation_xThe non-diagonal elements of (1) are zero, which indicates that the parameters after transformation become mutually independent, and is beneficial to improving the reliability of parameter inversion.

Based on a Bayesian inversion framework, an inversion target function is constructed by solving a maximum posterior probability density function, specifically to the prestack AVO inversion problem, the posterior probability density function can be expressed as:

P(R|d,I)＝const₀×P(d|R,I)P(R|I) (24)。

wherein P (d | R, I) is a likelihood function, P (R | I) is a prior distribution function, d represents prestack seismic data changing along with an incident angle, I represents basic geological information, R represents a model parameter to be inverted, const₀Is a probability normalization constant;

since the inversion idea based on Bayes theory only concerns the shape of the posterior probability density function, const₀Can be ignoredHowever, equation (24) can be further simplified as:

P(R|d,I)＝P(d|R,I)P(R|I) (25)。

the value of the maximum posterior probability density function is the maximum value of the product of the two random functions, and the uncertainty of the parameter to be estimated is represented by the width of posterior distribution; the likelihood function is mainly used for expressing the relation between observed seismic data and parameters to be inverted, forward records are obtained by means of a forward equation, and the likelihood function is constructed by researching the characteristics of errors (namely noise) between the observed data and the forward records; assuming that the seismic noise follows a gaussian distribution and the noises of different measurement conditions satisfy mutually independent conditions, the likelihood function is expressed as:

wherein σ_mRepresenting the standard deviation of the noise signal.

The prior function mainly represents the statistical characteristics of the parameters to be inverted, and Gaussian distribution and Cauchy distribution are distributed in a common prior way; in general, the Cauchy prior constraint is used under the condition that different elastic parameter distribution characteristics of different interfaces are assumed to be the same and are mutually independent, as discussed in the decorrelation processing of an equation of evolution, certain statistical correlation characteristics exist between the relative change rate of longitudinal and transverse wave impedance and the relative change rate of density parameters in an actual condition, so that the Cauchy prior regular constraint of a single variable is adopted, and the inversion quality of elastic parameters is influenced by neglecting the correlation characteristics of the parameters; the distribution characteristics of the relative change rate of the impedance of the longitudinal wave and the transverse wave and the relative change rate of the density are described by adopting four-variable Cauchy distribution on the assumption that the distribution of different interface parameters accords with the independent characteristics, so that the correlation characteristics among the four parameters are fully considered.

The prior function can be expressed as follows:

wherein,C_xis a covariance matrix, D_iIs a 4 × 4n matrix, and the values of the constituent elements are defined as follows:

therefore, substituting the likelihood function and the prior function into the posterior probability density function can obtain:

the constants are omitted, the maximum posterior probability is solved, and an inversion target function can be obtained, and the specific form is as follows:

(G'^TG'+2Q)R'＝G'^Td (30)。

wherein,the first item on the left side of the formula is mainly used for constraining the similarity degree between the forward record and the actual seismic record, the second item is a four-variable Cauchy regular constraint item and is mainly used for constraining the sparsity degree of inversion parameters, and then the target optimization is carried out on the inversion equation by adopting an iterative reweighted least square algorithm.

The calculation mode of the fluid factor is based on indirect combination operation of elastic parameters, and the quality of fluid identification mainly depends on two aspects, namely whether the inversion of the basic elastic parameters is reliable or not; the second is whether the constructed fluid factor is sensitive to the type of the pore fluid. Considering that pre-stack seismic inversion is a main means for extracting basic elastic parameters, the reliability of the elastic parameters can be better improved by improving an inversion method; on the other hand, the fluid factor based on indirect combination inevitably causes accumulated errors, the direct inversion of the fluid factor can better solve the problem, the sensitive fluid factor is constructed based on the double-phase medium rock physics theory, the direct inversion of the fluid factor can be realized by researching the internal relation between the sensitive fluid factor and the seismic motion reflection characteristic and utilizing pre-stack seismic data, and therefore the indication sensitivity of the fluid factor and the reliability of the estimation of the fluid factor are improved. The compressibility of the fluid and the solid can be quantitatively characterized, wherein the compressibility of the gas, the liquid and the solid has different degrees of compressibility, and the compressibility is defined as a relative reduction rate of pressure increase by one unit volume at a certain temperature, and the compressibility of the fluid and the solid can be quantitatively characterized, and generally, the compressibility of the gas, the liquid and the solid has the following relation: gas > liquid > solid. The method starts from the construction of pore fluid sensitive parameters, considers the difference between the compression coefficients of pore fluids (gas and liquid) and a rock framework, takes the fluid volume compression coefficient as a fluid factor, and further researches the inversion method of the fluid volume compression coefficient.

Drawings

Fig. 1 is a graph of the trend of ramet constant density with porosity and water saturation.

Fig. 2 is a graphical representation of the change in rem constant density with porosity and water saturation within the range of the dashed line in fig. 1 (water saturation between 10% and 40%).

FIG. 3 shows the fluid volume compressibility factor C_fTrend plots with porosity and water saturation.

FIG. 4 is a graph showing the volumetric compressibility C of the fluid within the range of the dotted line in FIG. 3 (water saturation of 10% -40%)_fSchematic of the changes with porosity and water saturation.

FIG. 5 is a graph of different fluid factor sensitivity analyses; wherein, 1-longitudinal wave velocity; 2-longitudinal wave impedance; 3-ramet parameter density; 4-shear modulus x density; 5-fluid factor; 6-Poisson's ratio; 7-transverse wave impedance; 8, transverse wave speed; 9-fluid volume compressibility.

FIG. 6 shows the reflection coefficient calculation results for model one.

Fig. 7 shows the calculation result of the reflection coefficient of model two.

FIG. 8 is a post-stack seismic section for use with the present invention.

Fig. 9 is a section of the inversion result of the fluid volume compressibility corresponding to fig. 7.

FIG. 10 is a detailed technical roadmap of the prestack seismic inversion method based on fluid volume compressibility of the present invention.

FIG. 11 is a flow chart of the present invention.

Detailed Description

The invention is further illustrated by the following figures and examples.

Example 1. The prestack seismic inversion method based on the fluid volume compression coefficient is characterized by comprising the following steps of:

step 1, establishing a relation between a fluid volume compression coefficient and other elastic parameters and physical parameters by using a theory and an empirical rock physical model based on a pore medium theory.

And 2, deriving an AVO approximate equation and an elastic impedance equation based on the fluid volume compression coefficient.

And 3, realizing direct inversion of the fluid volume compression coefficient under a Bayesian theory framework, and forming a relatively complete and stable pre-stack seismic inversion method based on the fluid volume compression coefficient.

And 4, applying the method to actual data.

In step 1, the technology firstly needs to establish the relation between the fluid volume compression coefficient and other elastic parameters and physical parameters by using theoretical and empirical rock physical models according to the pore medium theory.

Relative reduction rate of volume per unit pressure at a certain temperature

The dual-phase medium theory fully considers the relationship between the rock skeleton structure and the pore fluid property of the medium and the local characteristics and the overall effect, expresses the fluid-containing reservoir as a composite of a solid phase and a fluid phase, and respectively considers the influence of the solid and the fluid and the mutual coupling of the solid and the fluid on the seismic wave propagation. And the compression coefficient C of the two-phase medium is phi C_f+(1-φ)C_sIn which C is_f＝S_oC_o+S_wC_w+S_gC_g，C_g,C_oAnd C_wRespectively representing the fluid volume compression coefficients of gas, oil and water; s_g、S_oAnd S_wRespectively represent the gas, oil and water saturations, and S_g+S_w+S_o＝1。

Batzle and Han discovered through research that the main factors influencing the value of the Gassmann fluid item f are the bulk modulus of the pore fluid and the porosity of the solid skeleton. Dehua Han and Batzle et al studied the coupling effect of the solid effect (porosity, mineral modulus, etc.) of the pore fluid and the rock skeleton in the Biot-Gassmann theory on the rock modulus information by performing rock physics statistics on clastic rock, and proposed that the factor C for removing the rock skeleton fluid is highlighted_fThe empirical formula for the Gassmann fluid term of (1):

f＝G(φ)/C_f (2)。

wherein,wherein the gain function G (phi) represents the combined effect of the rock skeleton minerals and the porosity.

By combining an empirical relation obtained by rock physics experiments, the fluid volume compressibility coefficient C can be obtained by utilizing logging information_fAn estimation is performed. FIG. 1 shows C for different porosities and water saturations_fThe fluid volume compressibility C can be seen from the graph_fHas a completely linear change trend with the water saturation and is completely free from the influence of porosity. And as the gas saturation increases, the fluid volume compressibility factor C_fThe change is particularly obvious and indicates the fluid volume compressibility factor C_fThe fluid is very sensitive, and the accuracy of fluid identification is improved. Therefore, if we use a particular geophysical method to extract the fluid volume compressibility C from the seismic data_fCompressing the volume of the fluid by a factor C_fThe solid skeleton and the fluid elastic effect can be decoupled by taking the fluid factor as a fluid factor to participate in fluid identification, so that the reliability of reservoir pore fluid identification is effectively improved.

And (3) comparing the fluid factor of the detritus skeleton with the conventional fluid factor types (longitudinal and transverse wave impedance, Poisson ratio, Lambda [ rho ] and the like) and Gassmann fluid items based on the single-phase medium theory by combining logging data. Defining a certain fluid factor sensitivity as:

wherein mean is_gas,mean_waterThe average values, std, of the gas layer and the water layer corresponding to the fluid factor of the target interval on the well are_gasThen it is the standard deviation. Combining well data of a work area, performing sensitivity comparative analysis on conventional 9 fluid factors to obtain a graph 5It is seen that the fluid volume compressibility has the highest sensitivity to reservoir fluids. When the sensitivity is analyzed by using the logging data, an AVO approximation equation is not needed, and the general method is to use (or derive) the corresponding AVO approximation equation to perform subsequent work such as fluid identification after the fluid factor with the highest sensitivity is obtained.

The sensitivity analysis of rock is not used when the fluid volume compressibility is inverted, but used before. The sensitivity analysis of the rock guides the construction of the fluid volume compression coefficient, and aims to verify the practical feasibility of the constructed fluid volume compression coefficient. The fluid volume compressibility written in fig. 5 refers to the reversed fluid volume compressibility (which is different from the fluid volume compressibility used in the rock sensitivity analysis in terms of time), and the rock sensitivity analysis can see that the sensitivity of the fluid volume compressibility is high, so that the reversed fluid volume compressibility can be used as an indicator factor of fluid identification.

In step 2, an AVO seismic reflection characteristic approximation equation and an elastic impedance equation based on the fluid volume compressibility are derived by using the relationship between the fluid volume compressibility and the Gassmann fluid term.

Based on the Biot-Gassmann theory, Russell et al studied the AVO theory of saturated fluid porous media, and proposed the approximate formula of reflection coefficient including Gassmann fluid terms for the first time in the research report in 2006, and discussed it in the official papers published in Geophysics in 2011, and the Russell approximate formula is shown as follows:

wherein f, mu and rho respectively represent the average values of Gassmann fluid terms, shear modulus and density of media on two sides of the interface; Δ f, Δ μ, and Δ ρ represent the difference in shear modulus and density, respectively, of the Gassmann fluid terms on either side of the interface. Russell et al show through model trial calculation that the approximation formula meets the approximation accuracy requirement in the range of incidence angles less than 50 degrees.

Substituting the formula (2) into the Russell approximation formula (4), and performing corresponding transformation to obtain:

the dry rock shear modulus μ is used herein in view of its immunity to the pore fluid_dryInstead of μ, equation (5) can be further simplified to:

nur has shown through a great deal of research that for rocks with a porosity less than the critical porosity, the bulk modulus and shear modulus of the dry rock can be used in combination with the critical porosity phi_cThe relevant linear function is expressed, namely a critical porosity model, and the expression is shown as the following formula:

wherein phi is_cRepresents the critical porosity, K_dryDenotes the bulk modulus, μ, of the dry rock_dryDenotes the shear modulus, K, of the dry rock_mDenotes the bulk modulus, μ, of the solid mineral matrix_mRepresenting the shear modulus of the mineral matrix. Further expansion of equation (6) with the Nur model as the tie can result in:

will be provided withSubstituting G (φ), further expansion can result in:

suppose thatThen there is F_poro＝φ_cPhi mu is measured. The formula can be further simplified as:

if f is set_mAnd phi mu, the fluid volume compressibility approximation formula is finally obtained as follows:

wherein f is_mIs called the solid stiffness parameter. In order to verify the accuracy of the formula, two-layer sandstone models are designed by combining actual data of a clastic rock work area and using a fluid substitution method to test the clastic rock work area. The model parameters are shown in Table 2-1, where the critical porosity φ_cAll are 27%. The first model is that the upper sandstone contains gas, the lower sandstone contains water, and the porosity of the two layers is the same; and the two layers of sandstone of the second model are different in fluid and porosity, the upper layer of sandstone contains gas and has the porosity of 20%, and the lower layer of sandstone contains water and has the porosity of 25%. Calculating the reflection coefficient of the interface of the two models by respectively adopting an accurate Zoeppritz equation, an Aki-Richard approximation formula and an AVO approximation formula of the rock-removing skeleton, wherein the calculation results of the two models are respectively as errors! No reference source is found. 6. As shown in fig. 7.

TABLE 1 two sandstone model parameters

From fig. 6 and 7, it can be found that when the porosities of two sides of the interface are not changed and only the pore fluids are different, the accuracy of the AVO approximation of the rock-removing skeleton is better matched with that of Aki-Richard approximation, and both can be better matched with the accurate result of Zoeppritz; it can be seen from fig. 3 that when there is a difference in porosity and pore fluid across the interface, the derelike framework AVO approximation is still nearly identical to the Aki-Richard approximation, but both again deviate from the exact Zoeppritz equation calculations as the angle of incidence increases. The reason for this analysis is that the derricking skeleton AVO approximation is derived from the Aki-Richard approximation, so the consistency of the approximation accuracy of the two is not difficult to understand. Furthermore, the main reasons that both will generate large errors in case of porosity variation are: the AVO linear approximation formula is established on the assumption that small elastic disturbance exists on two sides of an interface, and the change of the porosity has large influence on rock modulus information, so that large errors are caused under the condition of large-angle incidence. However, based on the analysis result, it can be seen that the derricking framework AVO approximation can satisfy the approximation accuracy of the reflection coefficient under the condition of small-angle incidence, and considering that the angle gathers used in practical application generally do not exceed 30 degrees, the derricking framework AVO approximation can be considered to completely satisfy the requirement of the reflection coefficient accuracy within the error allowable range.

By using the idea of Connolly to derive elastic impedance, the reflection coefficient is expressed by the elastic impedance, and the following is obtained:

substituting equation (12) into (11) yields:

expressing the relative change of the elastic parameter in logarithmic form, we can get:

further variations, one can obtain:

namely, it is

And (3) integrating and indexing two sides of the above equation, eliminating differential terms and logarithmic terms on two sides of the equation, and further removing an integral constant of 0 to obtain:

EI(θ)＝(C_f)^a(θ)(f_m)^b(θ)(ρ)^c(θ)(φ)^d(θ) (17)。

wherein,

similar to the conventional elastic impedance formula, formula (17) also has a problem that the numerical dimension changes with angle. Four reference constants, namely A, are introduced here₀，C_f0，f_m0，ρ₀And phi₀. By normalizing equation (17), a normalized elastic impedance equation based on the detritus skeleton fluid factor can be obtained.

EI(θ)＝A₀(C_f/C_f0)^a(θ)(f_m/f_m0)^b(θ)(ρ/ρ₀)^c(θ)(φ/φ₀)^d(θ) (18)。

C_f0，f_m0，ρ₀And phi₀Are respectively defined as C_f0，f_mρ and φ. A. the₀For the normalization factor, the specific expression is:

in step 3, a Bayesian inversion framework is utilized, a likelihood function reflecting seismic information is combined with a prior geological constraint reflecting a parameter to be inverted, an inversion target function is established in a mode of solving a maximum posterior probability density function, conventional Bayesian AVO inversion is used for decorrelating an evolution equation without considering correlation characteristics among parameters in the prior constraint, in order to further improve inversion quality, four-variable Cauchy distribution is adopted as a prior regular constraint to improve Bayesian AVO four-parameter inversion on the basis of decorrelation processing of the evolution equation, and model trial calculation and practical application show that the reliability of four-parameter inversion is improved to a certain extent by the method. The bayesian inversion is an inversion algorithm, and the boxes under the black arrows mentioned in fig. 10 refer to some data that needs to be input in the bayesian inversion.

The reflection coefficient linear approximation formula is a bridge for establishing the relation between elastic parameters and seismic data, and the types and the information quantity of the elastic parameters contained in different reflection coefficient linear approximation formulas are different. The research adopts an AVO reflection coefficient linear approximation formula of the rock-removing framework to establish a forward model. The formula (11) is expressed as a matrix form according to the difference of the incident angles:

wherein, a_i(i＝1,2…,m)，b_i(i＝1,2…,m)，c_i(i-1, 2 …, m) and d_i(i ═ 1,2 …, m) denotes the corresponding coefficient for the ith angle of incidence, respectively. The method is popularized to the condition that m incident angles and n interfaces exist, and the matrix is subjected to blocking processing to obtain the following results:

wherein R is_i(i ═ 1,2 …, m) represents the reflection coefficient vector for the ith angle of incidence, consisting of n elements; a. the_i(i＝1,2…,m)，B_i(i＝1,2…,m)，C_i(i-1, 2 …, m) and D_i(i ═ 1,2 …, m) respectively represent forward modeling coefficient matrices corresponding to the i-th incident angle, and are respectively diagonal matrices of n × n dimensions;R_ρand R_φRespectively representing the fluid factor of the rock-removing framework, the solid rigidity parameter, the density and the relative change rate vector of the porosity, and respectively consisting of n elements.

With the wavelet matrix W introduced based on the assumption that the seismic records conform to the convolution model, equation (21) becomes further:

wherein d is_i(i ═ 1,2 …, m) column vectors consisting of seismic data for the i-th angle of incidence, all containing n elements.

In order to improve the solving stability of the inversion problem, some scholars introduce covariance matrixes in the process of carrying out mathematical characterization on the forward problem to carry out decorrelation processing on parameters to be inverted. Generating a covariance matrix according to a sample statistical method of logging data, namely:

wherein,n represents the number of samples.

Covariance matrix C of parameters to be inverted after decorrelation transformation_xThe non-diagonal elements of (1) are zero, which indicates that the parameters after transformation become mutually independent, and is beneficial to improving the reliability of parameter inversion.

Based on a Bayesian inversion framework, an inversion target function is constructed by solving a maximum posterior probability density function, specifically to the prestack AVO inversion problem, the posterior probability density function can be expressed as:

P(R|d,I)＝const₀×P(d|R,I)P(R|I) (24)。

wherein P (d | R, I) is a likelihood function, P (R | I) is a prior distribution function, d represents prestack seismic data changing along with an incident angle, I represents basic geological information, R represents a model parameter to be inverted, const₀Is a probability normalization constant.

Since the inversion idea based on Bayes theory only concerns the shape of the posterior probability density function, const₀Can be omitted, then equation (24) can be further simplified to:

P(R|d,I)＝P(d|R,I)P(R|I) (25)。

the value of the maximum posterior probability density function is the maximum value of the product of the two random functions, and the width of posterior distribution characterizes the uncertainty of the parameter to be estimated. The likelihood function is mainly used for expressing the relation between observed seismic data and parameters to be inverted, forward records are obtained by means of a forward equation, and the likelihood function is constructed by researching the characteristics of errors (namely noise) between the observed data and the forward records. Assuming that the seismic noise follows a gaussian distribution and the noises of different measurement conditions satisfy mutually independent conditions, the likelihood function is expressed as:

wherein σ_mRepresenting the standard deviation of the noise signal.

The prior function mainly represents the statistical characteristics of the parameters to be inverted, and Gaussian distribution and Cauchy distribution are distributed in a relatively common prior distribution. In general, the cauchy prior constraint used is that different elastic parameter distribution characteristics of different interfaces are assumed to be the same and are mutually independent, as discussed in the decorrelation processing of an equation of evolution, a certain statistical correlation characteristic exists between the relative change rate of the longitudinal and transverse wave impedance and the relative change rate of the density parameter in an actual situation, and therefore, the cauchy prior regular constraint of a single variable is adopted, and the inversion quality of the elastic parameter is influenced by ignoring the correlation characteristic of the parameter. The distribution characteristics of the relative change rate of the impedance of the longitudinal wave and the transverse wave and the relative change rate of the density are described by adopting four-variable Cauchy distribution on the assumption that the distribution of different interface parameters accords with the independent characteristics, so that the correlation characteristics among the four parameters are fully considered.

The prior function can be expressed as follows:

wherein,C_xis a covariance matrix, D_iIs a 4 × 4n matrix, and the values of the constituent elements are defined as follows:

therefore, substituting the likelihood function and the prior function into the posterior probability density function can obtain:

the constants are omitted, the maximum posterior probability is solved, and an inversion target function can be obtained, and the specific form is as follows:

(G'^TG'+2Q)R'＝G'^Td (30)。

wherein,the first item on the left side of the formula is mainly used for constraining the similarity between the forward record and the actual seismic record, the second item is a four-variable Cauchy regular constraint item and is mainly used for constraining the sparsity of inversion parameters, and then an Iterative weighted least square algorithm (IRLS) is adopted to carry out target optimization on the inversion equation.

The compressibility of the fluid and the solid can be quantitatively characterized, wherein the compressibility of the gas, the liquid and the solid has different degrees of compressibility, and the compressibility is defined as a relative reduction rate of pressure increase by one unit volume at a certain temperature, and the compressibility of the fluid and the solid can be quantitatively characterized, and generally, the compressibility of the gas, the liquid and the solid has the following relation: a gaseous liquid solid. The method is based on the construction of pore fluid sensitive parameters, considers the difference between the compression coefficients of pore fluid (gas and liquid) and a rock framework, takes the fluid volume compression coefficient as a fluid factor, further researches an inversion method of the fluid volume compression coefficient, and has practical significance for improving the reservoir fluid identification precision.

The method comprises the steps of firstly, carrying out sensitivity analysis on the fluid volume compressibility based on the rock physics theory and well logging information, and verifying the fluid identification accuracy of the fluid volume compressibility. Then, based on the pore elasticity theory, by utilizing theoretical and empirical rock physics models, the influence of rock modulus, pore fluid and pore size on reservoir fluid is fully considered, a plane wave reflection characteristic approximation equation based on a fluid volume compressibility is deduced, and the application conditions and the sensitivity of the equation are analyzed, wherein the application conditions and the sensitivity of the equation correspond to the conditions from the left to the right in a two-phase medium theory → a fluid volume compressibility → an AVO approximation equation → an elastic impedance equation in fig. 10, the two-phase medium theory is one of the pore medium theories, and the plane wave reflection characteristic approximation equation is the deduced AVO approximation equation.

Example 2. FIG. 10 is a technical roadmap for a pre-stack seismic inversion method based on fluid volume compressibility of the present invention. The main process of this embodiment:

constructing a fluid volume compression coefficient based on a pore elasticity theory;

deducing a seismic reflection characteristic equation based on the fluid volume compression coefficient;

a prestack seismic inversion method of fluid volume compression coefficients under a Bayesian framework;

the actual data application.

Step 101: fluid volume compression coefficient construction based on pore elasticity theory

Relative reduction rate of volume per unit pressure at a certain temperature

Batzle and Han discovered through research that the main factors influencing the value of the Gassmann fluid item f are the bulk modulus of the pore fluid and the porosity of the solid skeleton. Dehua Han and Batzle et al studied the coupling effect of the solid effect (porosity, mineral modulus, etc.) of the pore fluid and the rock skeleton in the Biot-Gassmann theory on rock modulus information by performing rock physics statistics on clastic rock, and proposed that the volume compressibility coefficient C of the fluid is highlighted_fThe empirical formula for the Gassmann fluid term of (1):

f＝G(φ)/C_f (2)。

wherein,wherein the gain function G (phi) represents the combined effect of the rock skeleton minerals and the porosity.

By comparing the fluid volume compressibility with the conventional fluid factor type based on single-phase medium theory in combination with well log data, as shown in fig. 5, it can be seen that the fluid volume compressibility has a very high sensitivity in the work area.

Step 102: derivation of seismic reflection characteristic equation based on fluid volume compression coefficient

Based on the Biot-Gassmann theory, Russell et al studied the AVO theory of saturated fluid porous media, and proposed the approximate formula of reflection coefficient including Gassmann fluid terms for the first time in the research report in 2006, and discussed it in the official papers published in Geophysics in 2011, and the Russell approximate formula is shown as follows:

according to the related rock physical model and the calculation method, the fluid volume compression coefficient approximation formula is finally obtained as shown in the following formula:

by using the idea of Connolly to derive elastic impedance, the reflection coefficient is expressed by the elastic impedance, and the following is obtained:

substituting equation (12) into (11) yields:

expressing the relative change of the elastic parameter in logarithmic form, we can get:

further variations, one can obtain:

namely, it is

And (3) integrating and indexing two sides of the above equation, eliminating differential terms and logarithmic terms on two sides of the equation, and further removing an integral constant of 0 to obtain:

EI(θ)＝(C_f)^a(θ)(f_m)^b(θ)(ρ)^c(θ)(φ)^d(θ) (17)。

wherein,

similar to the conventional elastic impedance formula, formula (17) also has a problem that the numerical dimension changes with angle. Four reference constants, namely A, are introduced here₀，C_f0，f_m0，ρ₀And phi₀. By normalizing equation (17), a normalized elastic impedance equation based on the volume compressibility of the fluid can be obtained.

EI(θ)＝A₀(C_f/C_f0)^a(θ)(f_m/f_m0)^b(θ)(ρ/ρ₀)^c(θ)(φ/φ₀)^d(θ) (18)。

C_f0，f_m0，ρ₀And phi₀Are respectively defined as C_f0，f_mρ and φ. A. the₀For the normalization factor, the specific expression is:

the accuracy of equation (19) can be tested by different models, as shown in fig. 6 and 7, the offset/incidence angle of the acquired seismic data is within a certain range, and the maximum incidence angle is generally less than 35 degrees. The thick dotted line in fig. 6 and 7 is a reflection coefficient curve of the exact equation, which describes a theoretical AVO characteristic, the thinner dotted line represents a reflection coefficient curve obtained by the AKI approximation equation, which describes a theoretical AVO characteristic under linear approximation, and the solid line is a reflection coefficient curve obtained by the newly derived approximation equation in this embodiment. In the range of the incident angle less than 30 °, it can be seen from fig. 6 that the 3 curves are almost consistent, and in fig. 7, although there is a difference between the 3 curves, and the main reason for the difference is that the difference between the two physical properties of the model is significant, so it can be considered that the error is within the allowable range, and therefore the new process has high accuracy in a certain angle range, so it can be applied to practical production.

Step 103: pre-stack seismic inversion method for fluid volume compression coefficient under Bayes framework

And finally, the reflection coefficient linear approximation formula is a bridge for establishing the relation between the elastic parameters and the seismic data, and the types and the information quantity of the elastic parameters contained in different reflection coefficient linear approximation formulas are different. The invention adopts a fluid volume compression coefficient AVO reflection coefficient linear approximation formula to establish a forward model. The formula (11) is expressed as a matrix form according to the difference of the incident angles:

wherein, a_i(i＝1,2…,m)，b_i(i＝1,2…,m)，c_i(i-1, 2 …, m) and d_i(i ═ 1,2 …, m) denotes the corresponding coefficient for the ith angle of incidence, respectively. The method is popularized to the condition that m incident angles and n interfaces exist, and the matrix is subjected to blocking processing to obtain the following results:

wherein R is_i(i ═ 1,2 …, m) represents the reflection coefficient vector for the ith angle of incidence, consisting of n elements; a. the_i(i＝1,2…,m)，B_i(i＝1,2…,m)，C_i(i-1, 2 …, m) and D_i(i ═ 1,2 …, m) respectively represent forward modeling coefficient matrices corresponding to the i-th incident angle, and are respectively diagonal matrices of n × n dimensions;R_ρand R_φRespectively representing the fluid factor of the rock-removing framework, the solid rigidity parameter, the density and the relative change rate vector of the porosity, and respectively consisting of n elements.

With the wavelet matrix W introduced based on the assumption that the seismic records conform to the convolution model, equation (21) becomes further:

wherein d is_i(i ═ 1,2 …, m) column vectors consisting of seismic data for the i-th angle of incidence, all containing n elements.

Step 104 actual data application

Under the Bayes theory frame, an elastic impedance inversion method is used for obtaining an elastic impedance inversion result of the fluid volume compression coefficient, and a direct relation between the elastic impedance and the fluid volume compression coefficient is established by using a formula (22), so that a final result of the fluid volume compression coefficient can be obtained. In a specific solving process, based on a Bayesian theory framework, a likelihood function and a prior function are substituted into a posterior probability density function, so that the following can be obtained:

the constants are omitted, the maximum posterior probability is solved, and an inversion target function can be obtained, and the specific form is as follows:

(G'^TG'+2Q)R'＝G'^Td (30)。

wherein,the first item on the left side of the formula is mainly used for constraining the similarity between the forward record and the actual seismic record, the second item is a four-variable Cauchy regular constraint item and is mainly used for constraining the sparsity of inversion parameters, and then an Iterative weighted least square algorithm (IRLS) is adopted to carry out target optimization on the inversion equation.

The resulting elastic impedance results were used to extract the volumetric compressibility of the fluid, as shown in fig. 8 and 9, A, B, C being 3 wells, one of which was a water well and the other two wells were gas wells. The gas and the water are well separated on the inversion section of the fluid volume compression coefficient obtained by the method, and the accuracy of fluid identification is improved.

The embodiment relates to a prestack seismic inversion method based on fluid volume compression coefficients. The compressibility of the fluid and the solid can be quantitatively characterized, and generally, the compressibility of the gas, the liquid and the solid has the following relation: gas > liquid > solid, so from the construction of pore fluid sensitive parameter, consider the difference between the compression coefficients of pore fluid (gas, liquid) and rock skeleton, regard fluid volume compression coefficient as the fluid factor, have practical meaning to improve reservoir fluid identification accuracy. Firstly, sensitivity analysis is carried out on the fluid volume compressibility based on the rock physics theory and well logging information, and the sensitivity analysis is used for verifying the fluid identification precision of the fluid volume compressibility. And then based on a pore elasticity theory, utilizing theoretical and empirical rock physical models, fully considering the influence of rock modulus, pore fluid and pore size on reservoir fluid, deducing a plane wave reflection characteristic approximation equation based on a fluid volume compression coefficient, and analyzing the application condition and sensitivity of the equation so as to provide input data meeting the condition and select proper inversion parameters during inversion. In this embodiment, the AVO approximation equation is derived based on the plane wave assumption, and an accurate plane wave equation is difficult to be applied in practice, so that the AVO approximation equation needs to be derived under a certain assumption, so that the AVO approximation equation is a special form of the plane wave equation. And finally, deducing an elastic impedance equation based on the fluid volume compression coefficient, and establishing a seismic-elastic impedance-fluid volume compression coefficient inversion process under a Bayes theory framework to form a prestack seismic inversion method based on the fluid volume compression coefficient.

Finally, it is noted that the above-mentioned preferred embodiments illustrate rather than limit the invention, and that, although the invention has been described in detail with reference to the above-mentioned preferred embodiments, it will be understood by those skilled in the art that various changes in form and detail may be made therein without departing from the scope of the invention as defined by the appended claims.

Claims (1)

1. A prestack seismic inversion method based on fluid volume compression coefficients is characterized by comprising the following steps: establishing a relation between the volume compression coefficient of the fluid and other elastic parameters and physical parameters by using a theory and an empirical rock physical model based on a pore medium theory; on the basis of the relation, a reflection characteristic equation based on the volume compressibility of the fluid is deduced; establishing a relation between the seismic data and the fluid volume compression coefficient by using a reflection characteristic equation so as to realize inversion of the fluid volume compression coefficient;

the method comprises the following steps:

step 1: establishing a relation between the volume compression coefficient of the fluid and other elastic parameters and physical parameters by using a theory and an empirical rock physical model based on a pore medium theory; other elastic parameters comprise Gassmann fluid factor, shear modulus and density, and the physical parameters comprise water saturation, porosity and argillaceous content;

step 2: under the assumption of plane waves, according to the relationship between the fluid volume compression coefficient obtained in the step 1 and other elastic parameters or physical parameters, an AVO seismic reflection characteristic approximation equation and an elastic impedance equation based on the fluid volume compression coefficient are deduced;

step 3, establishing a seismic data-Bayes elastic impedance inversion-fluid volume compression coefficient inversion process by using the elastic impedance equation obtained in the step 2, and realizing direct inversion of the fluid volume compression coefficient under a Bayes theory framework;

after the inversion is finished, the method also comprises the step of identifying reservoir fluid by using the fluid volume compression coefficient obtained by the inversion;

during reservoir fluid identification, the method also comprises the steps of carrying out sensitivity analysis on the fluid volume compressibility in the step 1 based on the rock physics theory and logging information, and verifying the fluid identification precision of the fluid volume compressibility;

the method of the sensitivity analysis is as follows:

combining logging information, and comparing the fluid volume compression coefficient with the fluid factor type based on the conventional single-phase medium theory to obtain fluid sensitivity; defining a certain fluid factor sensitivity as:

wherein mean is_gas,mean_waterThe average values, std, of the gas layer and the water layer corresponding to the fluid factor of the target interval on the well are_gasThen is the standard deviation;

the fluid factor is one or a combination of more of longitudinal and transverse wave velocity, longitudinal and transverse wave impedance, Poisson ratio, Lame constant density, shear modulus density, Gassmann fluid factor and Poisson ratio;

in step 1, the pore media theory is a biphasic media theory;

relative reduction rate of volume per unit pressure at a certain temperature

And the compression coefficient C of the two-phase medium is phi C_f+(1-φ)C_sIn which C is_f＝S_oC_o+S_wC_w+S_gC_g，C_g,C_oAnd C_wRespectively representing the fluid volume compression coefficients of gas, oil and water; s_g、S_oAnd S_wRespectively represent the gas, oil and water saturations, and S_g+S_w+S_o＝1；

Coefficient of compression of fluid volume C_fThe empirical formula for the Gassmann fluid term of (1):

f＝G(φ)/C_f (2)

wherein,wherein the gain function G (phi) represents the comprehensive action of rock skeleton minerals and porosity;

by combining an empirical relation obtained by rock physics experiments, the fluid volume compressibility coefficient C can be obtained by utilizing logging information_fCarrying out estimation; therefore, the fluid volume compressibility C is extracted from the seismic data by a geophysical method_fCompressing the volume of the fluid by a factor C_fThe solid skeleton and the fluid elastic effect can be decoupled by taking the fluid factor as a fluid factor to participate in fluid identification, so that the reliability of reservoir pore fluid identification is effectively improved;

in step 2, an AVO seismic reflection characteristic approximation equation and an elastic impedance equation based on the fluid volume compression coefficient are deduced by utilizing the relationship between the fluid volume compression coefficient and the Gassmann fluid term;

the approximate equation for the reflection coefficient containing the Gassmann fluid term is shown below:

wherein f, mu and rho respectively represent the average values of Gassmann fluid terms, shear modulus and density of media on two sides of the interface; Δ f, Δ μ, and Δ ρ represent differences in shear modulus and density of the Gassmann fluid terms on both sides of the interface, respectively;

substituting the formula (2) into the Russell approximation formula (4), and performing corresponding transformation to obtain:

the dry rock shear modulus μ is used herein in view of its immunity to the pore fluid_dryInstead of μ, equation (5) can be further simplified to:

the critical porosity model expression is shown as follows:

wherein phi is_cRepresents the critical porosity, K_dryDenotes the bulk modulus, μ, of the dry rock_dryDenotes the shear modulus, K, of the dry rock_mDenotes the bulk modulus, μ, of the solid mineral matrix_mRepresenting the shear modulus of the mineral matrix; further expansion of equation (6) with the Nur model as the tie can result in:

will be provided withSubstituting G (φ), further expansion can result in:

suppose thatThen there is F_poro＝φ_cPhi mu; the formula can be further simplified as:

if f is set_mAnd phi mu, the fluid volume compressibility approximation formula is finally obtained as follows:

wherein f is_mPhi mu, which is called the solid stiffness parameter;

by using the idea of Connolly to derive elastic impedance, the reflection coefficient is expressed by the elastic impedance, and the following is obtained:

substituting equation (12) into (11) yields:

expressing the relative change of the elastic parameter in logarithmic form, we can get:

further variations, one can obtain:

namely, it is

And (3) integrating and indexing two sides of the above equation, eliminating differential terms and logarithmic terms on two sides of the equation, and further removing an integral constant of 0 to obtain:

EI(θ)＝(C_f)^a(θ)(f_m)^b(θ)(ρ)^c(θ)(φ)^d(θ) (17)；

wherein,

similar to the conventional elastic impedance formula, formula (17) also has a problem that the numerical dimension changes with angle; four reference constants, namely A, are introduced here₀，C_f0，f_m0，ρ₀And phi₀(ii) a The formula (17) is subjected to standardization processing, so that a standardized elastic impedance equation based on the fluid factor of the rock-removing framework can be obtained;

EI(θ)＝A₀(C_f/C_f0)^a(θ)(f_m/f_m0)^b(θ)(ρ/ρ₀)^c(θ)(φ/φ₀)^d(θ) (18)；

C_f0，f_m0，ρ₀and phi₀Are respectively defined as C_f0，f_mρ and φ; a. the₀For standardizing factors, specially expressingThe formula is as follows:

in step 3, a Bayesian inversion framework is utilized, a likelihood function reflecting seismic information is combined with a prior geological constraint reflecting a parameter to be inverted, and an inversion target function is established by solving a maximum posterior probability density function;

the formula (11) is expressed as a matrix form according to the difference of the incident angles:

a_iwherein i is 1,2 …, m; b_iWherein i is 1,2 …, m; c. C_iWherein i is 1,2 …, m; d_iWherein i is 1,2 …, m; respective coefficients representing the ith incident angle; the method is popularized to the condition that m incident angles and n interfaces exist, and the matrix is subjected to blocking processing to obtain the following results:

R_iwherein i is 1,2 …, m; a reflection coefficient vector representing an i-th incident angle, composed of n elements; a. the_iWherein i is 1,2 …, m; b is_iWherein i is 1,2 …, m; c_iWherein i is 1,2 …, m; d_iWherein i is 1,2 …, m; respectively representing forward modeling coefficient matrixes corresponding to the ith incidence angle, wherein the forward modeling coefficient matrixes are respectively oblique diagonal matrixes with dimensions of n multiplied by n;R_ρand R_φRespectively representing the fluid factor of the rock skeleton, the solid rigidity parameter, the density and the relative change rate vector of the porosity, and respectively consisting of n elements;

with the wavelet matrix W introduced based on the assumption that the seismic records conform to the convolution model, equation (21) becomes further:

d_iwherein i is 1,2 …, m; the column vectors composed of the seismic data expressed as the ith incidence angle comprise n elements;

generating a covariance matrix according to a sample statistical method of logging data, namely:

wherein,n represents the number of samples;

covariance matrix C of parameters to be inverted after decorrelation transformation_xThe non-diagonal elements of the data are zero, which shows that the parameters after transformation become mutually independent, and the reliability of parameter inversion is favorably improved;

based on a Bayesian inversion framework, an inversion target function is constructed by solving a maximum posterior probability density function, specifically to the prestack AVO inversion problem, the posterior probability density function can be expressed as:

P(R|d,I)＝const₀×P(d|R,I)P(R|I) (24)

wherein P (d | R, I) is a likelihood function, P (R | I) is a prior distribution function, d represents prestack seismic data changing along with an incident angle, I represents basic geological information, R represents a model parameter to be inverted, const₀Is a probability normalization constant;

since the inversion idea based on Bayes theory only concerns the shape of the posterior probability density function, const₀Can be omitted, then equation (24) can be further simplified to:

P(R|d,I)＝P(d|R,I)P(R|I) (25)

the value of the maximum posterior probability density function is the maximum value of the product of the two random functions, and the uncertainty of the parameter to be estimated is represented by the width of posterior distribution; assuming that the seismic noise follows a gaussian distribution and the noises of different measurement conditions satisfy mutually independent conditions, the likelihood function is expressed as:

wherein σ_mRepresents the standard deviation of the noise signal;

assuming that different interface parameter distributions conform to independent characteristics, describing the distribution characteristics of the relative change rate of the longitudinal and transverse wave impedance and the relative change rate of the density by adopting four-variable Cauchy distribution, thereby fully considering the relevant characteristics among the four parameters;

the prior function can be expressed as follows:

wherein,C_xis a covariance matrix, D_iIs a 4 × 4n matrix, and the values of the constituent elements are defined as follows:

therefore, substituting the likelihood function and the prior function into the posterior probability density function can obtain:

the constants are omitted, the maximum posterior probability is solved, and an inversion target function can be obtained, and the specific form is as follows:

(G'^TG'+2Q)R'＝G'^Td (30)

wherein,

and then performing target optimization on the inversion equation by adopting an iterative reweighted least square algorithm.