阅读说明：本技术 一种基于精确Zoeppritz方程的时移地震AVA差异反演方法 (Time-lapse seismic AVA difference inversion method based on accurate Zoeppritz equation ) 是由 廉培庆 段太忠 张文彪 李景叶 赵磊 李蒙 周林 刘彦锋 于 2019-05-05 设计创作，主要内容包括：本申请提供了一种基于精确Zoeppritz方程的时移地震AVA差异反演方法,包括以下步骤：步骤1,基于所述精确Zoeppritz方程建立差异数据正演模型；步骤2,基于贝叶斯理论对所述差异数据正演模型进行时移地震AVA多波联合差异反演,获取油藏弹性参数的变化。通过该反演方法,能够好的避免了近似公式带来的计算误差,直接利用差异数据进行差异反演,摆脱了分别反演做差带来的计算量的问题,有效的提高了运算效率,有利于实际数据的处理。(The application provides a time-lapse seismic AVA difference inversion method based on an accurate Zoeppritz equation, which comprises the following steps: step 1, establishing a difference data forward model based on the accurate Zoeppritz equation; and 2, performing time-lapse seismic AVA multi-wave joint difference inversion on the difference data forward model based on a Bayesian theory to obtain the change of the oil reservoir elastic parameters. By the inversion method, calculation errors caused by an approximate formula can be well avoided, difference inversion is directly carried out by using difference data, the problem of calculation amount caused by difference making of respective inversion is solved, the operation efficiency is effectively improved, and the method is favorable for processing actual data.)

1. A time-lapse seismic AVA difference inversion method based on an accurate Zoeppritz equation is characterized by comprising the following steps:

step 1, establishing a difference data forward model based on the accurate Zoeppritz equation;

and 2, performing time-lapse seismic AVA multi-wave joint difference inversion on the difference data forward model based on a Bayesian theory to obtain the change of the oil reservoir elastic parameters.

2. The method of claim 1, wherein the difference data forward model is represented by the following equation:

wherein, delta d is time-lapse seismic difference data, delta m is the change quantity of the elastic parameters of the oil reservoir, and m₁Is the elastic parameter corresponding to the basic data, and L is the elastic parameter m corresponding to the basic data based on the forward modeling operator of the Zoeppritz equation₁The first partial derivative of (a), e, is the error term directly related to the difference data.

3. A method according to claim 2, characterized in that step 2 may comprise the following sub-steps:

step 21, determining a likelihood function according to the difference data forward model in the step 1;

step 22, constructing a prior model, wherein the prior model comprises a Gaussian distribution term and a block constraint term of a differential Laplace distribution;

step 23, determining an inversion target function according to the likelihood function determined in step 21 and the prior model determined in step 22;

and 24, solving the inversion objective function to obtain the change of the oil reservoir elastic parameter.

4. The method of claim 3, wherein in step 21, the negative logarithm of the likelihood function is represented by:

where F (Δ d | Δ m) is the negative logarithm of the likelihood function, C_DConst is a constant term for the covariance matrix of the noise.

5. The method of claim 4, wherein in step 22, the negative logarithm of the prior model is represented by:

F(Δm)＝const+F₁(Δm)+F₂(Δm)

wherein, F (Δ m), F₁(Δ m) and F₂(Δ m) represents a negative logarithm of the prior model, a negative logarithm of a gaussian distribution term, and a block distribution term, respectivelyNegative logarithm of (C)_ΔmIs a covariance matrix containing the correlation of three difference data, mu is the mean vector of the difference data, N is the length of the model parameter, D is the first order differential operator, k_lIs a scale factor, l ═ 1,2 or 3.

6. The method of claim 5, wherein the inversion objective function is expressed as follows for multi-wave seismic data:

α＝σ_Δpp/σ_Δps

β＝σ_Δpp

wherein J (Δ m) is the inversion objective function, α controls the proportion of PS wave data, β controls the proportion of prior information, σ_ΔPPIs the noise variance, σ, of the PP-wave difference seismic data_ΔPSIs the noise variance of the PS-wave difference seismic data.

7. The method according to any one of claims 3 to 6, wherein in step 24, an iterative solution is performed using an iterative reweighted least squares method to obtain the change in reservoir elasticity parameters.

8. The method of any one of claims 1 to 6, wherein the reservoir elasticity parameters comprise compressional wave velocity, shear wave velocity, and density.

Technical Field

The invention relates to a method for predicting parameters of oil and gas fields, shale gas seismic exploration and reservoirs, in particular to a time-lapse seismic AVA difference inversion method based on an accurate Zoeppritz equation.

Background

The time-lapse seismic oil reservoir monitoring technology can measure changes of seismic response caused by oil and gas exploitation, then quantifies the changes by means of a seismic inversion technology to obtain changes of reservoir elastic parameters, and finally obtains changes of physical parameters such as reservoir oil-gas saturation, pressure, temperature and the like according to rock physical relations, and provides important seismic information for development of middle and later-stage residual oil distribution and reservoir fine description, so that the technology is widely applied to reservoir characterization, monitoring and management.

For time-lapse seismic inversion, Sarkar et al in 2003 (Sarkar S, Gouveia W P, Johnston DH.2003.on the inversion of time-lapse differential data.73th Annual International testing, SEG, Expanded Abstract, 1489-. So if repeatability between the base survey and the survey data is adequate, the difference inversion is the most appropriate time-lapse seismic inversion method. In order to further improve the accuracy of inversion, Theune et al studied a block inversion method for time-lapse seismic AVA (Amplitude Variation with azimuth of an acquisition plane) data on the basis of Bayes theory difference inversion in 2010 (Theune, U.S., Jenses I.O., and Eidsvik, J.2010, Analysis of prior models for a block inversion of semi-seismic AVA data: Geophysics, vol.75, No.3,25-35), and the actual data processing result shows that the method can better reflect the Variation of time-lapse elasticity parameters. The WA waveform difference inversion method based on the Bayesian theory is proposed in 2012 by WA GWA et al (WA GWA, WAVE, time-lapse seismic data Bayes AVO waveform inversion, geophysical science report, 2012,55 vol 7, 2422 and 2431), but the accuracy of an operator is greatly lost due to an approximate formula still adopted by a positive operator in the method. In addition, only depending on the improved Cauchy constraint in the prior constraint, the abrupt change condition of the boundary in the difference inversion is difficult to reflect, so the method has certain defects in the precision and accuracy of the difference inversion, and the block constraint needs to be added in the prior constraint to improve the inversion precision. The inversion method based on the Zoeppritz equation approximation formula limits the improvement of the pre-stack inversion accuracy in many aspects, so that the solution by using the precise Zoeppritz equation is important. In addition, under the condition of possessing good multi-wave data, the joint inversion can further enhance the stability of the inversion process and improve the accuracy of the inversion result.

In the methods described above, difference data is obtained by performing difference after respective inversion, or only inversion of non-difference data is performed, and the difference values of three parameters cannot be directly obtained by directly inverting the difference data. In addition, most of the methods construct a positive operator based on an approximate formula of the Zoeppritz equation, which reduces inversion accuracy to a certain extent.

Disclosure of Invention

Aiming at the problems, the invention provides a time-lapse seismic AVA difference inversion method based on a precise Zoeppritz equation, which utilizes the precise Zoeppritz equation to obtain a forward operator with high accuracy, applies Bayesian theory, adds block constraints which effectively highlight the formation boundary effect into a prior model, and applies multi-wave joint inversion, thereby obtaining a difference inversion result with high precision and high resolution, avoiding artificial errors introduced when the traditional time-lapse seismic inversion method is used for inverting the pre-earthquake and the post-earthquake respectively, providing technical support for direct difference inversion, and finally obtaining a difference value of an elastic parameter.

The time-lapse seismic AVA difference inversion method based on the accurate Zoeppritz equation provided by the invention comprises the following steps: step 1, establishing a difference data forward model based on the accurate Zoeppritz equation; and 2, performing time-lapse seismic AVA multi-wave joint difference inversion on the difference data forward model based on the Bayesian theory to obtain the change of the oil reservoir elastic parameters.

In a preferred embodiment, the difference data forward model is represented by the following equation:

wherein, delta d is time-lapse seismic difference data, delta m is the change quantity of the elastic parameters of the oil reservoir, and m₁Is the elastic parameter corresponding to the basic data, and L is the elastic parameter m corresponding to the basic data based on the forward modeling operator of the Zoeppritz equation₁The first partial derivative of (a), e, is the error term directly related to the difference data.

In a preferred embodiment, step 2 may comprise the following sub-steps: step 21, determining a likelihood function according to the difference data forward model in the step 1; step 22, constructing a prior model, wherein the prior model comprises a Gaussian distribution term and a block constraint term of a differential Laplace distribution; step 23, determining an inversion objective function according to the likelihood function determined in step 21 and the prior model determined in step 22; and 24, solving the inversion objective function to obtain the change of the oil reservoir elastic parameter.

In a preferred embodiment, in step 21, the negative logarithm of the likelihood function is represented by the following equation:

where F (Δ d | Δ m) is the negative logarithm of the likelihood function, C_DConst is a constant term for the covariance matrix of the noise.

In a preferred embodiment, in step 22, the negative logarithm of the prior model is represented by:

F(Δm)＝const+F₁(Δm)+F₂(Δm)

wherein, F (Δ m), F₁(Δ m) and F₂(Δ m) represents the negative logarithm of the prior model, the negative logarithm of the gaussian distribution term, and the negative logarithm of the block distribution term, C, respectively_ΔmIs a covariance matrix containing the correlation of three difference data, mu is the mean vector of the difference data, N is the length of the model parameter, D is the first order differential operator, k_lIs a scale factor, l ═ 1,2 or 3.

In a preferred embodiment, the inverted objective function is expressed for multi-wave seismic data by the following equation:

α＝σ_Δpp/σ_Δps

β＝σ_Δpp

wherein J (Δ m) is the inversion objective function, α controls the proportion of PS wave data, β controls the proportion of prior information, σ_ΔPPIs the noise variance, σ, of the PP-wave difference seismic data_ΔPSIs the noise variance of the PS-wave difference seismic data.

In a preferred embodiment, in step 24, an iterative solution is performed using an iterative reweighed least squares method to obtain the change in the reservoir elasticity parameter.

In a preferred embodiment, the reservoir elastic parameters include compressional wave velocity, shear wave velocity, and density.

The method constructs a forward operator of the time-lapse seismic difference data based on the accurate Zoeppritz equation, constructs an inversion target function from the forward operator, and well avoids calculation errors caused by an approximate formula. In addition, the difference inversion is directly carried out by using the difference data, so that the problem of calculation amount caused by difference making of respective inversion is solved, the operation efficiency is effectively improved, and the processing of actual data is facilitated; in addition, the method can stably and reasonably accurately acquire the change information of two elastic parameters of the longitudinal wave velocity and the transverse wave velocity from the difference seismic data, can also accurately estimate the variation of the density information, verifies the feasibility and the effectiveness of the new method, and provides important data for the fine description of the reservoir stratum and the well position deployment in the middle and later stages of development.

The features mentioned above can be combined in various suitable ways or replaced by equivalent features as long as the object of the invention is achieved.

Drawings

The invention will be described in more detail hereinafter on the basis of embodiments and with reference to the accompanying drawings. Wherein:

FIG. 1 shows a schematic flow diagram of a method for time-lapse seismic AVA difference inversion based on the exact Zoeppritz equation, in accordance with an embodiment of the present invention;

FIG. 2 shows a schematic flow of an inversion based on Bayesian theory according to an embodiment of the invention;

FIG. 3 shows a schematic comparison of well logs before and after single well fluid substitution in accordance with an embodiment of the present invention;

FIG. 4 shows true variance values of reservoir elasticity parameters and an initial model map according to an embodiment of the invention;

FIG. 5 shows a synthetic prestack angle gather from forward modeling using the exact Zoeppritz equation, in accordance with an embodiment of the present invention;

FIG. 6(a) shows the elastic parameter variation obtained by separate inversion and then difference in the prior art;

FIG. 6(b) shows the variation of elastic parameters obtained by the inversion method according to the embodiment of the invention;

FIG. 6(c) shows the variation of elastic parameters obtained by direct inversion of difference data using the inversion method of Aki-Richards approximation formula;

FIG. 7 is a schematic diagram showing the results of noisy time-lapse seismic difference inversion according to an embodiment of the present invention;

FIG. 8 is a graph showing the results of a joint inversion of time lapse seismic difference data based on the exact Zoeppritz equation, in accordance with an embodiment of the present invention.

In the drawings, like parts are provided with like reference numerals. The drawings are not to scale.

Detailed Description

The invention will be further explained with reference to the drawings.

FIG. 1 shows a time-lapse seismic AVA difference inversion method 100 based on the precise Zoeppritz equation provided by the present invention, as shown in FIG. 1, the method 100 comprises the following steps:

s110, establishing a difference data forward model based on the accurate Zoeppritz equation;

and S120, performing time-lapse seismic AVA multi-wave joint difference inversion on the difference data forward model based on the Bayesian theory to obtain the change of the oil reservoir elastic parameters.

Specifically, aiming at the characteristic that the distribution of the parameters to be inverted of the time-lapse earthquake is discontinuous, a more reasonable prior model is introduced, the vertical block constraint item which obeys differential Laplace distribution is added while the elastic parameters obey Gaussian distribution, so that the artificial errors caused by the traditional time-lapse earthquake inversion method when the pre-and post-development earthquakes are inverted respectively are avoided, the accuracy of the inversion result is improved, the change of the elastic parameters of the reservoir is highlighted, the technical support is provided for direct difference inversion, and the difference numerical value of the elastic parameters is finally obtained. In the inversion process, the multi-wave joint inversion is utilized, and the information of the PS wave is added, so that the inversion result is more accurate.

In step S110, the exact Zoeppritz equation describes the reflection and transmission of a plane wave incident at the interface of two semi-infinite homogeneous isotropic elastic media (Zoeppritz, 1919; Aki and Richards, 1980). The expression is as follows:

in time-lapse seismic inversion, for well-matched time-lapse seismic data, the change in seismic response is usually confined to a particular formation, whereas in the non-reservoir portion, the base and monitor data should be the same. The difference data is obtained by subtracting the monitoring data and the basic data, so that necessary coupling can be obtained for the two times of seismic data acquired at different times, the difference data obtained by subtracting is directly inverted, and further the change of the oil deposit elastic parameters caused by mining is obtained.

Forward modeling processes of time-lapse seismic basic data and monitoring data are respectively as follows:

d₁＝G(m₁)+n₁(2)

d₂＝G(m₂)+n₂(3)

wherein G is a positive operator of the Zoeppritz equation. d₁As basic data, d₂To monitor the data. m is₁Based on the corresponding elastic parameter, m₂For monitoring the corresponding elastic parameter of the data, n₁And n₂Respectively, noise data at the time of two data acquisitions.

Corresponding equation (3) to the elasticity parameter m of the basic data₁The position of the stent is subjected to Taylor expansion,

subtracting equation (2) from equation (4),

order to

Then there is

Δd＝LΔm+e (6)

Where Δ d is the difference seismic data, Δ m is the amount of change in the elastic parameter, L (m)₁) Based on the Zoeppritz equationThe elastic parameter m of the positive operator corresponding to the three elastic parameters of longitudinal, transverse and density in the basic data₁The first partial derivative of (a), e, is an error term directly related to the difference data, which may also be understood as noise.

In S120, time-lapse seismic AVA multi-wave joint difference inversion is carried out on the difference data forward model based on Bayesian theory, and the change of the oil reservoir elastic parameters is obtained. Specifically, this step may include the following sub-steps:

s121, determining a likelihood function according to the difference data forward model in the step S110;

s122, constructing a prior model, wherein the prior model comprises a Gaussian distribution term and a block constraint term of a derivative Laplace distribution;

s123, determining an inversion target function according to the likelihood function determined in the step S121 and the prior model determined in the step S122;

and S124, solving the inversion objective function to obtain the change of the oil reservoir elastic parameter.

First, in determining the likelihood function, assuming that the error term e in equation (6) follows a zero-mean gaussian distribution, the likelihood function can be expressed as,

in the above formula, C_DIs the covariance matrix of the noise, N_dIs the length of the observed data. The negative logarithm is taken for the above formula,

time-lapse seismic inversion generally only inverts the variation of the elastic parameter of the target interval, so that the elastic parameter of the small segment can be considered to be normally distributed. Due to the blocking phenomenon of the time-lapse seismic difference data, a block constraint term is added into the prior model, and according to the study of Theune U et al (2010), the differential Laplace distribution is assumed to be obeyed. In summary, the prior model consists of two parts: (1) a gaussian distribution term containing a priori low frequency trends and internal covariances between different model parameters; (2) a block constraint term characterized by a long tail distribution.

Assuming that the prior model follows a gaussian distribution, then:

the negative logarithm is taken for the above formula,

wherein, C_ΔmThe covariance matrix containing the correlation of the three difference data, N is the length of the model parameter, and μ is the mean vector of the difference data (three elastic parameters need to be obtained separately). In order to obtain drastic elastic parameter changes, the study adds vertical block constraints to the prior model, and then the negative logarithm of the prior model is defined as,

F(Δm)＝const+F₁(Δm)+F₂(Δm) (11)

wherein F₂(Δ m) represents a block constraint term, which is expressed as follows:

where D is a first order differential operator, k_lWhere l is 1,2,3 is a scale factor, the three elasticity parameters may differ from one another.

Based on Bayes theory, constant terms are reduced, and an inversion target function shown as follows can be obtained:

in practical processing, it is generally assumed that the noise of the observed data is uncorrelated, and the covariance matrix of the noise is reduced to a diagonal matrix, i.e. C_D＝σ_n ²I, whereinσ_n ²Representing the variance of the noise, I being N_d×N_dIdentity matrix of N_dIs the length of the observed data. For multi-wave seismic data, the objective function can be expanded, and the noise variance of PP wave difference seismic data is assumed to be sigma_ΔPPNoise variance of PS wave difference seismic data is σ_ΔPSThen the above-mentioned objective function can be expressed as,

wherein α ═ σ_Δpp/σ_ΔpsControlling the proportion of PS wave data, beta ═ sigma-_ΔppControlling the proportion of the prior information.

The above objective function is derived with respect to am and the derivative is made zero,

the mixture is obtained by finishing the raw materials,

wherein the content of the first and second substances,

since the right-end element a of the above equation also contains Δ m, an iterative reweighted least squares algorithm is used for iterative solution. Due to the positive operator L (m)₁) The reservoir elastic parameter m corresponding to the three elastic parameters of longitudinal wave velocity and transverse wave velocity and density in basic data based on a positive operator of an accurate Zoeppritz equation is contained₁The first partial derivative, and therefore the pre-stack inversion based on the exact Zoeppritz equation needs to be performed on the matched base data to obtain the positive mathematical submatrix L. This term can be obtained by means of a nonlinear AVO three-parameter inversion method based on the exact Zoeppritz equation, proposed by Zhou et al (2016).

However, those skilled in the art will appreciateCompared with the inversion algorithm based on an approximate formula, the pre-stack AVA inversion algorithm based on the Zoeppritz equation has a slightly larger calculation amount, and according to the above, the time-lapse seismic difference inversion needs to be performed on the basis of the one-time pre-stack inversion, and if the processed actual work area data is larger, the calculation amount can be increased. Through calculation and comparison, the positive operator L can find that the value m is taken by the positive operator L₁And m₀The matrix element difference calculated in the case of (inversion result after independent inversion by using basic data) is small, so that the calculation amount can be saved

To approximate to replace

In order to verify the feasibility and stability of the inversion algorithm, the applicant has carried out practice in exploration work. The solid line in fig. 3 is the log before fluid replacement, i.e., before development, and the dashed line is the log after fluid replacement of the actual reservoir region using the Gassmann equation, where the abscissa represents the wave velocity and the ordinate represents time. In fig. 4, the solid line represents the variation of the real elastic parameter obtained by subtracting the well-log curves before and after development shown in fig. 3, and the dotted line represents the initial model obtained by smoothing the real variation multiple times, wherein the abscissa represents the wave velocity and the ordinate represents the time. A forward evolution with the exact Zoeppritz equation resulted in a synthetic pre-stack angle gather (10 angles between 4 ° and 40 ° at 4 ° intervals) as shown in fig. 5, where the abscissa represents angle and the ordinate represents time. In the forward process, a Rake wavelet with the main frequency of 30Hz is adopted, and the sampling interval is 2 ms. In order to check the noise immunity of the inversion algorithm, random noise with a signal-to-noise ratio (root-mean-square amplitude ratio) of 2 is added to the synthesized basic data and the monitored data, and inversion is performed by using the algorithm.

Fig. 6a shows the elastic parameter variation obtained by subtracting the synthesized basic data and the monitored data after being respectively inverted by using an inversion method based on accurate Zoeppritz, wherein the abscissa represents the wave velocity, the ordinate represents the time (the same below), and fig. 6b shows the elastic parameter variation obtained by using the new inversion method, and it can be seen from the comparison between fig. 6a and fig. 6b that the inversion result of the difference data is obviously better than the inversion result, especially the transverse wave velocity and the density. To further compare the superiority of the new method, the results of conventional time-lapse seismic inversion are presented herein, and fig. 6c is the variation of elastic parameters obtained by direct inversion of difference data using an inversion method based on the Aki-Richards approximation formula. Since the approximate formula has limited calculation accuracy and many assumptions, it can be seen from the comparison of fig. 6b and 6c that the accuracy of the inversion result of the new method is significantly higher than that of the conventional inversion result. Fig. 7 is an inversion result obtained by adding random noise with a signal-to-noise ratio of 2 to the basic data and the monitored data, and performing inversion by using a research method, wherein the abscissa represents the wave velocity and the ordinate represents the time. In addition, if there is well-matched multi-wave seismic data, the method can also be used for performing multi-wave joint inversion, fig. 8 is a longitudinal and transverse wave joint inversion result based on synthetic data, wherein the abscissa represents wave velocity, and the ordinate represents time, and as can be seen from the comparison between fig. 8 and fig. 6b, joint inversion can further enhance the stability of the inversion process and improve the accuracy of the inversion result.

The method constructs a forward operator of the time-lapse seismic difference data based on the accurate Zoeppritz equation, constructs an inversion target function from the forward operator, and well avoids calculation errors caused by an approximate formula. In addition, the difference inversion is directly carried out by using the difference data, the problem of calculation amount caused by difference making of respective inversion is solved, the operation efficiency is effectively improved, and the method is favorable for processing actual data.

In the description of the present invention, it is to be understood that the terms "upper", "lower", "bottom", "top", "front", "rear", "inner", "outer", "left", "right", and the like, indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, are only for convenience in describing the present invention and simplifying the description, and do not indicate or imply that the device or element being referred to must have a particular orientation, be constructed in a particular orientation, and be operated, and thus, should not be construed as limiting the present invention.

Although the invention herein has been described with reference to particular embodiments, it is to be understood that these embodiments are merely illustrative of the principles and applications of the present invention. It is therefore to be understood that numerous modifications may be made to the illustrative embodiments and that other arrangements may be devised without departing from the spirit and scope of the present invention as defined by the appended claims. It should be understood that features described in different dependent claims and herein may be combined in ways different from those described in the original claims. It is also to be understood that features described in connection with individual embodiments may be used in other described embodiments.