CN105467441A - Device using average incidence angle trace gathers to carry out PP save and PS wave combined AVO inversion - Google Patents

Abstract

The invention provides a device using average incidence angle trace gathers to carry out PP save and PS wave combined AVO inversion. The device comprises a first calculation module used for calculating a compression coefficient for compressing PS waves to a time domain of PP waves; a compression module used for applying the compression coefficient to PS wave AVA trace gathers and obtaining PS wave AVA trace gathers in the time domain of the PP waves; a selecting module used for selecting PP wave AVA trace gathers and the compressed PS wave AVA trace gathers as AVA trace gathers respectively of a first incidence angle range and a second incidence angle range; a superposing and extracting module used for carrying out weighted superposing on respectively AVA trace gathers of first and second incidence angle ranges of the PP waves and the PS waves, forming the average incidence angle trace gathers, and respectively extracting seismic sub-waves of the PP waves and the PS waves in the same time window; and a second calculation module used for extrapolating a P wave speed, an S wave speed and a density obtained by logging into an initial model so as to calculate reflective coefficients of the PP waves and the PS waves, and forming the average incidence angle trace gathers by combining the reflective coefficients with the seismic sub-waves. The device provided by the invention solves the problem that an existing inversion method has a relatively large error.

Description

Device for performing PP wave and PS wave combined AVO inversion by using average incidence angle gather

Technical Field

The invention relates to the technical field of seismic exploration, in particular to a device for performing PP wave and PS wave combined AVO inversion by using an average incidence angle gather.

Background

Along with the continuous deepening of the exploration and development of oil and gas fields in China, the constructed oil and gas reservoirs are less and less, and unconventional oil and gas resources such as compact oil, shale gas, coal bed gas and the like become important directions for exploration. In the process, technical means of seismic exploration are also continuously developed to deal with the increasingly complex oil and gas exploration problems. Among many new seismic exploration technologies, the multi-component seismic technology has obvious advantages in solving the exploration problem of complex reservoirs. A number of theoretical studies have shown that: the PP wave and the PS wave obtained by the multi-component earthquake are utilized to carry out joint AVO inversion, so that reliable longitudinal and transverse wave speed and density information can be provided, more attribute parameters reflecting the lithology and the fluid-containing property of the stratum are derived, the multi-solution of the earthquake explanation is reduced to a great extent, and a new thought is provided for the prediction of a complex oil and gas reservoir.

The basis of the multi-wave AVO inversion is to solve the seismic reflection coefficient by a Zoeppritz equation, but the equation is complex in form and does not give an "explicit" relationship between the reflection coefficient and the formation elasticity parameters. Therefore, an important aspect of research on AVO inversion at home and abroad is to approximate the Zoeppritz equation to obtain a linear relationship between the reflection coefficient and the elastic parameters of the stratum. Among these, the most commonly used approximation formula is the Aki-Richards approximation formula, which is based on the assumption that the elastic parameters of both sides of the formation interface are weak contrasts, and can be used to directly reverse the longitudinal and transverse wave velocities and density variations of both sides of the single interface from the reflection coefficient. The theory is firstly broken through on longitudinal waves and is fully applied to actual exploration.

Although the AVO inversion combining the PP wave and the PS wave has obvious advancement in theory, the effect is not obvious in practical application. The reason mainly exists in the following two aspects: (1) firstly, the signal-to-noise ratio of the PS wave incident angle gather is poor, and the reliability of the joint inversion is reduced. AVO inversion needs to obtain a high-quality prestack gather as a guarantee, but at present, the processing level is difficult to ensure that the AVA gather of the PS wave has a wider incidence angle range, and the covering times of each incidence angle are uniform and the signal-to-noise ratio is higher. (2) Another reason is that the reflection coefficient approximation formula used in the inversion is not suitable for the formation interface with strong contrast on the assumption of weak contrast. The difference between the elastic parameters of a plurality of reservoirs, such as coal-bearing formations, unconsolidated sandstone formations, igneous rock and shale formations and the like, and the surrounding rock is large, so that the AVO inversion method based on the Aki-Richards reflection coefficient approximation formula has large errors.

Disclosure of Invention

The invention mainly aims to provide a device for performing PP wave and PS wave combined AVO inversion by using an average incident angle gather so as to overcome the problem of larger error of the conventional inversion method.

In order to solve the above problems, the apparatus for performing a PP wave and PS wave joint AVO inversion by using an average incident angle gather according to an embodiment of the present invention includes a first calculation module, a compression module, a selection module, a superposition and extraction module, and a second calculation module, wherein the first calculation module calculates a compression coefficient of a time domain when a PS wave is compressed to a PP wave; the compression module is connected with the first calculation module to apply the compression coefficient to a PS wave AVA gather to obtain a PS wave AVA gather of a PP wave time domain; the selection module is connected with the compression membrane block to sort the PP wave and the compressed PS wave AVA gather into AVA gathers in a first incidence angle range and a second incidence angle range; the stacking and extracting module is connected with the selecting module to perform weighted stacking on AVA gathers of the first and second incidence angle ranges of the PP wave and the PS wave respectively to form an average incidence angle gather, and seismic wavelets of the PP wave and the PS wave are extracted respectively in the same time window, and the second calculating module is connected with the stacking and extracting module to extrapolate the P wave velocity, the S wave velocity and the density obtained by well logging into an initial model to calculate the reflection coefficients of the PP wave and the PS wave and synthesize the average incidence angle gather with the seismic wavelets.

The input data used by the method is an average incident angle gather, and a uniform objective function is utilized to combine residual errors between a forward modeling gather of PP waves and PS waves and an actual data average incident angle gather together so as to calculate accurate reflection coefficients of the PP waves and the PS waves and perform least square inversion on longitudinal wave velocity, transverse wave velocity and density at the same time, so that the parameters obtained by the method are accurate.

Drawings

The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this application, illustrate embodiment(s) of the invention and together with the description serve to explain the invention without limiting the invention. In the drawings:

FIG. 1a and FIG. 1b are AVO curves of PP wave reflection coefficient and PS wave reflection coefficient calculated by Zoeppritz equation, respectively;

FIG. 2 is a block diagram of an apparatus for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the present invention;

FIG. 3 is another block diagram of an apparatus for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the present invention;

FIG. 4 is a flow chart of a method for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the present invention;

FIG. 5 is another flow diagram of a method for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the present invention;

FIGS. 6a and 6c are waveform diagrams of AVA gathers of PP and PS waves, respectively;

FIGS. 6b and 6d are waveform diagrams of PP and PS channel sets composed of incident angles of 5 and 20 degrees, respectively;

FIG. 7a, FIG. 7b, and FIG. 7c are schematic diagrams of the results of inversion of 5 degree and 20 degree angle gathers of PP and PS waves, respectively, using the inversion apparatus and method of the present invention;

FIGS. 8a, 8b, and 8c are inversion results based on approximate reflection coefficients, respectively;

FIGS. 9a and 9c are waveform diagrams of a PP wave and a PS wave channel set after 20% random noise is added, respectively;

FIGS. 9b and 9d are waveform diagrams of PP and PS channel sets composed of incident angles of 5 degrees and 20 degrees, respectively, and added with 20% of random noise;

fig. 10a, 10b, and 10c are schematic diagrams of the results obtained by inversion of 5-degree and 20-degree angle gathers of PP and PS waves with 20% random noise added by the inversion apparatus and method of the present invention, respectively.

Detailed Description

The main idea of the invention is that the used input data is an average incident angle gather, a uniform target function is utilized to combine residual errors between a forward gather of PP waves and PS waves and an actual data average incident angle gather to calculate accurate reflection coefficients of the PP waves and the PS waves, and simultaneous least square inversion of longitudinal wave velocity, transverse wave velocity and density is carried out, and further, in the inversion process, an initial model is modified through multiple iterations to enable the target function to reach the minimum, so that the parameters obtained by the invention are accurate.

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings and specific embodiments.

First, in an isotropic medium, when a planar longitudinal wave is obliquely incident to a boundary surface of two media, four kinds of waves, i.e., a reflected longitudinal wave, a reflected transverse wave, a transmitted longitudinal wave, and a transmitted transverse wave, are generated, and Snell's law (Snell's law) is satisfied, as shown in equation (1.1):

$\frac{V_{p1}}{\sin \mathrm{\α}}=\frac{V_{s1}}{\sin \mathrm{\β}}=\frac{V_{p2}}{{\mathrm{sin\α}}^{\′}}=\frac{V_{s2}}{{\mathrm{sin\β}}^{\′}}.---\left(1.1\right)$

wherein p is a ray parameter; v_P1And V_S1Longitudinal wave speed and transverse wave speed of the upper medium respectively; v_P2And V_S2Respectively the longitudinal wave speed and the transverse wave speed of the lower medium; rho₁And ρ₂Density of upper and lower layers of medium, α respectively₁And α₂Incident angle (reflection angle) and transmission angle of longitudinal wave β₁And β₂The shear wave reflection angle and the transmission angle.

Let the longitudinal wave reflection coefficient be R_PPConverting the transverse wave reflection coefficient to R_PSLongitudinal wave transmission coefficient of T_PPAnd converting the transverse wave transmission to T_PSA zeppartz (Zoeppritz) equation system can be derived that the displacement amplitudes of the four waves should satisfy, as shown in equation (1.2):

$\left(\begin{array}{cccc}-{\mathrm{sin\α}}_1& -{\mathrm{cos\β}}_1& {\mathrm{sin\α}}_2& -{\mathrm{cos\β}}_2\\ {\mathrm{cos\α}}_1& -{\mathrm{sin\β}}_1& {\mathrm{cos\α}}_2& {\mathrm{sin\β}}_2\\ -\cos 2{\mathrm{\β}}_1& \frac{V_{S1}}{V_{P1}}{\mathrm{sin\β}}_1& \frac{{\mathrm{\ρ}}_2}{{\mathrm{\ρ}}_1}\frac{V_{P2}}{V_{P1}}{\mathrm{cos\β}}_2& \frac{{\mathrm{\ρ}}_2}{{\mathrm{\ρ}}_1}\frac{V_{S2}}{V_{P1}}{\mathrm{sin\β}}_2\\ \sin 2{\mathrm{\α}}_1& \frac{V_{P1}}{V_{S1}}{\mathrm{cos\β}}_1& \frac{{\mathrm{\ρ}}_2}{{\mathrm{\ρ}}_1}\frac{V_{S2}^2}{V_{S1}^2}\frac{V_{P1}}{V_{P2}}\sin 2{\mathrm{\α}}_2& -\frac{{\mathrm{\ρ}}_2}{{\mathrm{\ρ}}_1}\frac{V_{S2}{V}_{P1}}{V_{S1}^2}\cos 2{\mathrm{\β}}_2\end{array}\right)\left(\begin{array}{c}{R}_{PP}\\ R_{PS}\\ T_{PP}\\ T_{PS}\end{array}\right)=\left(\begin{array}{c}{\mathrm{sin\α}}_1\\ {\mathrm{cos\α}}_1\\ \cos 2{\mathrm{\β}}_1\\ \sin 2{\mathrm{\α}}_1\end{array}\right)---\left(1.2\right)$

based on the formula (1.2), the elastic parameters of the media on two sides of the interface are given, and the accurate reflection coefficients of the PP wave and the PS wave can be calculated through a numerical solution of a linear equation system.

In order to verify whether lithology and oil-gas-bearing property of the stratum can be reversely deduced by combining the reflection coefficients of a small incidence angle and a large incidence angle on AVO curves of a PP wave and a PS wave, three equivalent stratum models are established, model parameters are shown in a table 1, and the AVO curves calculated accurately through a formula (1.2) are shown in a figure 1a and a figure 1 b. Taking a small incident angle of 5 degrees and a large incident angle of 20 degrees, corresponding points on the reflectance curves of the PP wave and the PS wave are picked up in fig. 1a and 1b, respectively, as shown by the dotted lines in fig. 1a and 1b, and the "two points" are connected in a straight line, and it can be seen that: the combination mode of the stratum cannot be distinguished by using the connection line of the two-point reflection coefficients of the PP wave (the dotted line in figure 1 a) alone, and particularly, the connection line of the two-point reflection coefficients of the shale/water saturated sandstone interface and the shale/gas saturated sandstone interface is very close to each other and cannot be distinguished; only the difference between the PP wave two-point reflection coefficient connecting line of the limestone/gas saturated sandstone interface and other two connecting lines is larger. Therefore, the reliability of AVO inversion cannot be guaranteed by singly using the average incident angle gathers of two angles of the PP wave.

Comparing fig. 1b with fig. 1a, it can be seen that: the "two-point" reflection coefficient connection of the PS wave alone (the dashed line in fig. 1 b) is difficult to distinguish between the shale/gas saturated sandstone interface and the limestone/gas saturated sandstone interface, but the "two-point" reflection coefficient connection of the PS wave of the shale/water saturated interface shows a greater difference from the other two connections. Therefore, the true condition of the oil and gas bearing property of the stratum cannot be inverted by using the average incidence angle gather of two angles of the PS wave alone. However, if the PP wave and the PS wave "two-point" reflection coefficient link are combined together, corresponding to an incident angle gather forming four angles, the three interface types can be completely distinguished. The conclusion also provides a theoretical basis for the joint inversion of the average incidence angle gathers of the PP wave and the PS wave.

The average incident angle gathers are equivalent to limiting the angles of the AVA gathers of the PP and PS waves to two angles each. Wherein the smaller average incident angle seismic traces are obtained by weighted superposition of seismic data of a small incident angle range in the AVA trace set; similarly, a larger average angle of incidence seismic trace is obtained by weighted superposition of seismic data for a relatively larger range of angles of incidence. The multi-wave local superposition data of the average incidence angle can be used for inverting the elastic parameters. In order to make the amplitude of the average incidence angle seismic trace approximate to the true seismic amplitude of the angle, the prestack data processing should adopt amplitude preserving algorithm as much as possible. The average incidence angle is selected according to the actual situation of multi-component seismic data acquisition, but in order to keep the variation trend of the reflection coefficients of the PP wave and the PS wave, the large and small average incidence angles are selected to be different as much as possible, and the average incidence angles are operable in actual data processing.

Under the condition that the relation between the longitudinal wave velocity and the density meets a Gardner formula and the longitudinal wave velocity change rate and the transverse wave velocity change rate of the stratums on two sides of the interface are close, the optimal average incident angle of the independent inversion of the PP wave is 0 degree; the far-offset incidence angle of the PS wave is selected, and the 30-degree average incidence angle is suitable for the single inversion of the PS wave in consideration of the factors of absorption attenuation and dynamic correction stretching ratio of the stratum. However, in the joint inversion of the PP wave and the PS wave, the dominant amplitude of the PP wave superposition is concentrated in medium and small offset distances, and the dominant amplitude of the PS wave superposition is concentrated in medium and large offset distances; the 0 and 30 degree average incidence angle seismic traces are not suitable for both PP and PS wave inversion. Since the reflection coefficient of the PS wave at an incident angle of 0 degree is 0, the joint inversion should select a small incident angle of non-zero degree, and therefore the incident angle is preferably 5 degrees. For large incidence angles, the trade-off is that the co-inversion is preferably chosen to have a large incidence angle of 20 degrees, given that the reflection coefficient of the PP wave may be small at large incidence angles.

TABLE 1 elastic parameters of three fluid-bearing formation models

In addition, in order to enable the inversion method to be applicable to wider stratum conditions, accurate reflection coefficients are calculated based on the Zoeppritz equation to carry out joint inversion of the PP wave and the PS wave, so that the inversion accuracy is improved. The joint inversion objective function Q, established by the least squares method, is shown in equation (1.3):

Q(V)＝||W_PP*R_PP-D_PP||²+||W_PS*R_PS-D_PS||²(1.3)

wherein, $D_{PP}=\left(D_{PP}^{\left({\mathrm{\θ}}_1\right)},D_{PP}^{\left({\mathrm{\θ}}_2\right)}\right),D_{PS}=\left(D_{PS}^{\left({\mathrm{\θ}}_1\right)},D_{PS}^{\left({\mathrm{\θ}}_2\right)}\right)$ the average incident angles of the actual PP wave and the PS wave are theta₁And theta₂The seismic record vector of (a); w_PP、W_PSSeismic wavelet vectors of PP waves and PS waves respectively,the average incident angles of the PP wave and the PS wave are respectively theta₁And theta₂The reflection coefficient of (2). V ═ V (V)_P,V_Sρ) is a vector of formation elastic parameters on both sides of the interface and has V_P＝(V_P1,V_P2)、V_S＝(V_S1,V_S2)、ρ＝(ρ₁,ρ₂)。

Taylor expansion is carried out on the reflection coefficient vectors of the PP wave and the PS wave near the initial model, and the unified matrix expression form is shown as a formula (1.4):

$R(V_0+\mathrm{\Δ}V)=R\left(V_0\right)+\frac{\∂R}{\∂V}\·\mathrm{\Δ}V---\left(1.4\right)$

wherein, V₀Vectors formed from initial guesses of compressional velocity, shear velocity and density△ V is the initial model correction quantity, let Jacobian matrixSubstituting equation (1.4) into equation (1.3) yields equation (1.5), as follows:

Q(V)＝||W_PP*(R_PP0+G_PP△V)-D_PP||²+||W_PS*(R_PS0+G_PS△V)-D_PS||²(1.5)

writing equation (1.5) as an expanded form of the matrix, as shown in equation (1.6):

Q(V)＝(W_PP*R_PP0-D_PP)^T(W_PP*R_PP0-D_PP)+(W_PP*R_PP0-D_PP)^T(W_PP*G_PP△V)

+(W_PP*G_PP△V)^T(W_PP*R_PP0-D_PP)+(W_PP*G_PP△V)^T(W_PP*G_PP△V)(1.6)

+(W_PS*R_PS0-D_PS)^T(W_PS*R_PS0-D_PS)+(W_PS*R_PS0-D_PS)^T(W_PS*G_PS△V)

+(W_PS*G_PS△V)^T(W_PS*R_PS0-D_PS)+(W_PS*G_PS△V)^T(W_PS*G_PS△V)

in order to minimize the value of the objective function Q (V), the requirement existsEquation (1.6) can therefore be converted to equation (1.7) as follows:

$\begin{array}{c}\frac{\∂Q\left(V\right)}{\∂\mathrm{\Δ}V}=\[{(W_{PP}*R_{PP0}-D_{PP})}^T(W_{PP}*G_{PP})+{(W_{PP}*G_{PP}\mathrm{\Δ}V)}^T(W_{PP}*G_{PP})\]\\ +\[{\left(W_{PS}*R_{PS0}-D_{PS}\right)}^T\left(W_{PS}*G_{PS}\right)+{\left(W_{PS}*G_{PS}\mathrm{\Δ}V\right)}^T\left(W_{PS}*G_{PS}\right)\]=0\end{array}---\left(1.7\right)$

further, Δ V can be solved as shown in equation (1.8):

△V＝[(W_PP*G_PP)^T(W_PP*G_PP)+(W_PS*G_PS)^T(W_PS*G_PS)]^-1(1.8)

·[(W_PP*G_PP)^T(D_PP-W_PP*R_PP0)+(W_PS*G_PS)^T(D_PS-W_PS*R_PS0)]

wherein G is_PP、G_PSRespectively are Jacobian matrixes of the PP wave and the PS wave; r_PP0、R_PS0Each represents R_PP、R_PSThe zero order term for Taylor expansion at the initial model △ V in the above equation is compared to the initial elastic parameter vector V₀Carrying out vector summation to obtain an updated initial model, and carrying out iteration of the process; when the value of the objective function Q (V) is minimum, the inversion process is finished. Updated elastic parameter vector V₀The inversion result of the model is obtained.

In practical joint inversion, a time series consisting of a series of sampling points is input, and the vector in the formula (1.8) needs to be expanded to be suitable for inversion of a time series model. Extended jacobian matrix G_PPAnd G_PSHas a uniform form, as shown in equation (1.9):

$G=\left[\begin{array}{cccc}\frac{\∂R_1}{\∂V_1}& \frac{\∂R_1}{\∂V_2}& L& \frac{\∂R_1}{\∂V_n}\\ \frac{\∂R_2}{\∂V_1}& \frac{\∂R_2}{\∂V_2}& L& \frac{\∂R_2}{\∂V_n}\\ M& M& O& M\\ \frac{\∂R_n}{\∂V_1}& \frac{\∂R_n}{\∂V_2}& L& \frac{\∂R_n}{\∂V_n}\end{array}\right]---\left(1.9\right)$

the subscripts 1, 2 … n in equation (1.9) represent the various time sample points, R_iThe reflection coefficient of the PP wave or the PS wave representing the position of the ith sampling point,the reflection coefficient of the PP wave or the PS wave at the ith sampling point is represented by the variation of the elastic parameter at the jth sampling point, and has the following formula (1.10):

$\frac{\∂R_i}{\∂V_j}=\left(\frac{\∂R_i}{\∂V_{Pj}},\frac{\∂R_i}{\∂V_{Sj}},\frac{\∂R_i}{\∂{\mathrm{\ρ}}_j}\right)---\left(1.10\right)$

in the above equation, the solution of the reflection coefficient must be numerically solved by the equation (1.2) for the system of linear equations; if Aki-Richards reflection coefficient approximation formula is used to establish the Jacobian matrix of formula (1.9), the inversion method of the present invention degenerates to an inversion method based on Aki-Richards reflection coefficient approximation formula.

In order to modify the elastic parameters of the time sampling points simultaneously in one iteration, the average incident angle gather of the PS wave must be on the same time scale as the PP wave, so the PS wave gather must be compressed to the time of the PP wave. In the process, the comparison and interpretation of the horizon can be carried out through the superposition section of the PP wave and the PS wave so as to obtain the compression coefficients of different horizons, thereby realizing the compression of the PS wave. If the difference of the same-phase axis appearances of the PP wave and the PS wave reflected waves of the same geological interface is too large, enough logging synthetic records in a work area are needed to assist in completing the matching of the layer positions. In addition, the seismic wavelet needs to be extracted after the PS wave average incident angle gather is compressed to PP wave time, so that the consistency of the time scale of the seismic wavelet is ensured. Because of the difference in wavelet types between PP and PS waves, the PP and PS waves must be separately wavelet extracted for inversion.

While the above description has been provided to illustrate the related equations needed for implementing the embodiments of the present invention, the following description will be provided to illustrate the corresponding embodiments. According to the embodiment of the invention, a device for performing combined AVO inversion of PP waves and PS waves by using average incidence angle gathers is provided.

FIG. 2 is a block diagram of an apparatus for performing a PP-wave and PS-wave joint AVO inversion with an average incident angle gather according to an embodiment of the present invention. The device 200 for performing PP wave and PS wave joint AVO inversion by using an average incident angle gather includes a first calculating module 210, a compressing module 220, a selecting module 230, a stacking and extracting module 240, and a second calculating module 250.

The first calculation module 210 calculates a compression coefficient for compressing the PS wave into the PP wave time domain.

The compression module 220 is connected to the first calculation module 210, so as to apply the compression coefficient to the PS-wave AVA gather to obtain the PS-wave AVA gather of the PP-wave time domain.

The selecting module 230 is connected to the compressing module 220 to select the PP wave and the compressed PS wave AVA gathers into AVA gathers of a first incident angle range and a second incident angle range.

The stacking and extracting module 240 is connected to the selecting module 130 to weight and stack the AVA gathers of the first and second incident angle ranges of the PP wave and the PS wave, respectively, to form an average incident angle gather, and extract the seismic wavelets of the PP wave and the PS wave in the same time window, respectively.

The second computation module 250 is connected to the stacking and extracting module 240 to extrapolate the P-wave velocity, S-wave velocity and density obtained by logging into an initial model to compute the reflection coefficients of the PP wave and the PS wave, and synthesize an average incident angle gather with the seismic wavelets.

FIG. 3 is another block diagram of an apparatus for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the present invention.

The apparatus 300 for performing PP wave and PS wave joint AVO inversion by using an average incident angle gather includes a first calculating module 210, a compressing module 220, a selecting module 230, a stacking and extracting module 240, a second calculating module 250, a third calculating module 310, a determining module 320, and a control module 330. The connection relationship and operation of the first calculating module 210, the compressing module 220, the selecting module 230, the superimposing and extracting module 240, and the second calculating module 250 can refer to the embodiment of fig. 2, and therefore, the description thereof is omitted.

The third calculation module 310 is coupled to the second calculation module 250 to calculate an objective function based on the average angle of incidence gathers of the synthetic seismic records with actual data.

The determining module 320 determines whether the objective function is smaller than a predetermined value to generate a determination result.

The control module 330 is connected to the determining module 320, the first calculating module 210, the compressing module 220, the selecting module 230, the superimposing and extracting module 240, and the second calculating module 250, so as to output a model modifier vector when the determination result is that the objective function is smaller than a predetermined value. When the determination result is that the objective function is not less than the preset value, the control module 330 generates an initial model correction amount to correct the initial model, and controls the first calculation module 210, the compression module 220, the selection module 230, the superposition and extraction module 240, and the second calculation module 250 to restart, so that the initial model is modified through multiple iterations until the objective function reaches a minimum value.

In addition, according to an embodiment of the invention, a method for performing PP wave and PS wave joint AVO inversion by using an average incident angle gather is provided.

FIG. 4 is a flow chart of a method for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the invention.

In step S402, a compression coefficient of the PS wave compressed to the PP wave time domain is calculated.

Step S404, the compression coefficient is applied to the PS wave AVA gather to obtain the PS wave AVA gather of the PP wave time domain, so that the average incidence angle gather of the PS wave and the PP wave are on the same time scale.

Step S406, the PP wave and the compressed PS wave AVA gathers are sorted into AVA gathers of a first incidence angle range and a second incidence angle range. The first incident angle range is different from the second incident angle range, that is, the first incident angle range is smaller than the second incident angle range, or the first incident angle range is larger than the second incident angle range.

Step S408, the AVA gathers of the first and second incidence angle ranges of the PP wave and the PS wave are weighted and superposed respectively to form an average incidence angle gather, and the seismic wavelets of the PP wave and the PS wave are extracted respectively in the same time window.

And S410, extrapolating the P-wave velocity, the S-wave velocity and the density obtained by logging into an initial model to calculate the reflection coefficients of the PP wave and the PS wave, and synthesizing an average incident angle gather with the seismic wavelets. The reflection coefficients of the PP wave and the PS wave are obtained according to, for example, formula (1.2).

FIG. 5 is another flow chart of a method for performing a PP-wave and PS-wave joint AVO inversion with an average angle of incidence gather according to an embodiment of the invention.

In step S402, a compression coefficient of the PS wave compressed to the PP wave time domain is calculated.

Step S404, the compression coefficient is applied to the PS wave AVA gather to obtain the PS wave AVA gather of the PP wave time domain, so that the average incidence angle gather of the PS wave and the PP wave are on the same time scale.

Step S406, the PP wave and the compressed PS wave AVA gathers are sorted into AVA gathers of a first incidence angle range and a second incidence angle range. The first incident angle range is different from the second incident angle range, i.e. the first incident angle range is smaller than the second incident angle range, or the first incident angle range is larger than the second incident angle range

Step S408, the AVA gathers of the first and second incidence angle ranges of the PP wave and the PS wave are weighted and superposed respectively to form an average incidence angle gather, and the seismic wavelets of the PP wave and the PS wave are extracted respectively in the same time window.

And S410, extrapolating the P-wave velocity, the S-wave velocity and the density obtained by logging into an initial model to calculate the reflection coefficients of the PP wave and the PS wave, and synthesizing an average incident angle gather with the seismic wavelets. The reflection coefficients of the PP wave and the PS wave are obtained according to, for example, formula (1.2).

Step S502, calculating an objective function based on the average incident angle gather of the synthetic seismic record and the actual data. Wherein the objective function can be solved using equation (1.3).

Step S504, judging whether the objective function is smaller than a preset value, and generating a judgment result.

Step S506, when the judgment result is that the objective function is smaller than a preset value, outputting a model modifier vector.

Step S508, when the determination result is that the objective function is not smaller than the preset value, generating an initial model correction amount to correct the initial model, and returning to step S410, repeating the above process, and stopping when the objective function reaches the minimum. The initial model correction amount is formula (1.8), and can be solved by formulas (1.3) to (1.7).

In order to verify the inversion device and the inversion method, the invention provides a stratum model containing a thin layer, wherein model parameters are shown in table 2, the thickness of the 3 rd layer is thinner, the difference between the thickness of the 3 rd layer and the elastic parameters of surrounding rocks is larger, and both top and bottom interfaces are strong contrast interfaces; other bed boundaries are weak contrast interfaces. Synthetic recordings were made using ray tracing based on exact reflection coefficient simulations. The seismic wavelets of the PP waves and the PS waves are respectively set as 40Hz and 30Hz Rake wavelets, and the sampling rate is 1 ms; the reflection time of the PS wave has been compressed to the PP wave time.

Fig. 6a and 6c are waveform diagrams of AVA gathers of PP and PS waves, respectively, and fig. 6b and 6d are waveform diagrams of PP and PS gathers composed of 5 and 20 degree incident angles, respectively. Fig. 7a, 7b, and 7c are schematic diagrams of results obtained by inversion of 5-degree and 20-degree angle gathers of PP waves and PS waves, respectively, by using the inversion apparatus and method of the present invention. Wherein, reference numeral 701 is an initial model curve, reference numeral 702 is a real model curve, and reference numeral 703 is an average incident angle gather joint inversion curve.

Fig. 7a, 7b, 7c show that: under the condition that the initial model has no interface change, the inversion device and the inversion method based on the accurate reflection coefficient can obtain an inversion result which is very close to the real model; especially for the thin layer of the 3 rd layer, the sudden change of the inversion curve at the interface position is better matched with the real model. In general, compared with a real model, the inversion device and the inversion method have the advantages that the model curve obtained by the inversion device and the inversion method is slightly smooth at the position where the elastic parameter of the stratum suddenly changes, but the error is small.

Fig. 8a, 8b, and 8c are inversion results based on the approximate reflection coefficient. Wherein, reference numeral 801 is an initial model curve, reference numeral 802 is a real model curve, and reference numeral 802 is an average incident angle gather joint inversion curve. Fig. 8a, 8b, 8c compare fig. 7a, 7b, 7c to see that: the inversion based on the approximate reflection coefficient cannot well reflect the real change rule of the elastic parameters of the stratum; even on the interface with weak contrast, the result error is large, and particularly the inversion result of the longitudinal wave velocity of the 2 nd layer is obviously small and the density is obviously large. The results of the above numerical analysis showed that: for some stratums with large elastic parameter contrast and thin thickness, AVO inversion must be carried out by adopting an accurate reflection coefficient; the accuracy of AVO inversion based on Zoeppritz' equation exact reflection coefficients is high even for weakly contrasted formations.

After 20% of random noise is added to the seismic record, the influence of the noise on the inversion device and the inversion method are tested. As shown in FIGS. 9b and 9d, 5-degree and 20-degree average incident angle seismic traces are obtained by weighted stacking of seismic traces with incident angles of 0-10 degrees and seismic traces with incident angles of 15-25 degrees, respectively, on the basis of noisy AVA gathers (FIGS. 9a and 9 c). The inversion device and the inversion method of the invention are used for respectively carrying out inversion of precise reflection coefficients based on a Zoeppritz equation on an AVA gather and an average incident angle gather, as shown in fig. 10a, 10b and 10c, an average incident angle gather inversion node 1001 and an AVA gather inversion result 1003 are both close to a real model 1002; however, the noise amplitude of the AVA gather inversion result 1003 is stronger than the average incident angle gather inversion junction 1001. Therefore, direct inversion of an AVA gather is more robust to noise than inversion with an average incident angle gather.

In summary, according to the technical solution of the present invention, the input data is an average incident angle gather, and a uniform objective function is used to combine residuals between a forward gather of PP and PS waves and an actual average incident angle gather of data, so as to calculate accurate reflection coefficients of PP and PS waves, perform simultaneous least square inversion of longitudinal wave velocity, transverse wave velocity and density, and further modify the initial model through multiple iterations in the inversion process, so that the objective function is minimized, and the parameters obtained by the present invention are more accurate.

The above description is only an example of the present invention, and is not intended to limit the present invention, and it is obvious to those skilled in the art that various modifications and variations can be made in the present invention. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the claims of the present invention.

Claims (2)

1. A device for performing PP wave and PS wave combined AVO inversion by using an average incident angle gather comprises a first calculation module, a compression module, a selection module, a superposition and extraction module and a second calculation module, wherein the first calculation module calculates a compression coefficient of a time domain from the compression of a PS wave to the PP wave; the compression module is connected with the first calculation module to apply the compression coefficient to a PS wave AVA gather to obtain a PS wave AVA gather of a PP wave time domain; the selection module is connected with the compression membrane block to sort the PP wave and the compressed PS wave AVA gather into AVA gathers in a first incidence angle range and a second incidence angle range; the stacking and extracting module is connected with the selecting module to perform weighted stacking on AVA gathers of the first and second incidence angle ranges of the PP wave and the PS wave respectively to form an average incidence angle gather, and seismic wavelets of the PP wave and the PS wave are extracted respectively in the same time window, and the second calculating module is connected with the stacking and extracting module to extrapolate the P wave velocity, the S wave velocity and the density obtained by well logging into an initial model to calculate the reflection coefficients of the PP wave and the PS wave and synthesize the average incidence angle gather with the seismic wavelets.

2. The apparatus of claim 1, further comprising a third computing module, a determining module and a control module, wherein the third computing module is connected to the second computing module to compute an objective function based on the average incident angle gather of the synthetic seismic record and the actual data, the determining module is connected to the third computing module to determine whether the objective function is smaller than a predetermined value and generate a determination result, the control module is connected to the determining module, the first computing module, the compressing module, the selecting module, the stacking and extracting module and the second computing module to output a model modifier vector when the determination result is that the objective function is smaller than a predetermined value, and generate an initial model modifier when the determination result is that the objective function is not smaller than the predetermined value, so as to correct the initial model and control the first calculation module, the compression module, the selection module, the superposition and extraction module and the second calculation module to restart.