Longitudinal and transverse wave combined prestack AVO inversion method and related device

Abstract

The invention discloses a longitudinal and transverse wave combined pre-stack AVO inversion method and a related device. The method comprises the steps that a work area elastic parameter prediction model is established according to work area logging, well drilling, vertical seismic profile VSP and geological data, and elastic parameters comprise P-wave velocity of longitudinal waves, S-wave velocity of transverse waves and medium density; according to the elastic parameter prediction model, the pre-stack AVO reflection coefficient equation and the angle gather seismic wavelets, establishing an angle gather formation record; according to the angle gather forming record and the angle gather actual record, establishing a longitudinal and transverse wave joint prestack AVO inversion iteration equation; and performing joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correcting the elastic parameter prediction model until a preset termination condition is met, and determining the elastic parameter of the pre-stack AVO inversion. The problem that existing pre-stack elastic parameter inversion is inaccurate is solved, higher accuracy is achieved in the application of directly locking the oil-gas-containing abnormal target area, and technical support is provided for oil-gas exploration and development.

Description

Technical Field

[0001] The present invention relates to the field of vector elastic wave exploration, and in particular to a method and related device for joint pre-stack AVO inversion of longitudinal and transverse waves. Background Art

[0002] Seismic inversion is to inversely deduce the distribution of rock geophysical parameters such as wave impedance, velocity, and density of underground media based on seismic data, so as to estimate reservoir parameters and conduct reservoir prediction, in order to provide reliable basic data for the exploration and development of oil and gas fields. The variation of seismic wave amplitude with offset (AVO) inversion belongs to pre-stack inversion, which analyzes the variation characteristics of amplitude with offset on the pre-stack CDP gather to identify lithology and predict oil and gas, and is an important technology for studying lithology and identifying oil and gas using amplitude information. Currently, seismic wave AVO inversion mainly adopts P-wave AVO inversion, converted-wave AVO inversion, or joint AVO inversion of P-wave and converted-wave, but these methods cannot obtain accurate shear wave velocity parameters, causing errors in oil and gas prediction work. Summary of the Invention

[0003] In view of the above problems, the purpose of the present invention is to provide a method and related device for joint pre-stack AVO inversion of seismic longitudinal and transverse waves.

[0004] In a first aspect, an embodiment of the present invention provides a method for joint pre-stack AVO inversion of seismic longitudinal and transverse waves, including the following steps:

[0005] According to well logging, drilling, vertical seismic profile (VSP) and geological data of the work area, establish an elastic parameter prediction model for the work area, where the elastic parameters include longitudinal wave P-wave velocity, shear wave S-wave velocity, and medium density;

[0006] According to the elastic parameter prediction model, the pre-stack AVO reflection coefficient equations of longitudinal wave P-wave, pure shear wave SV-wave, and pure shear wave SH-wave, and the angle gather seismic wavelet, establish an angle gather synthetic record;

[0007] According to the angle gather synthetic record and the actual angle gather record, establish a joint pre-stack AVO inversion iteration equation;

[0008] Perform joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correct the elastic parameter prediction model until a preset termination condition is met, and determine the elastic parameters of the pre-stack AVO inversion.

[0009] In one embodiment, the performing joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correcting the elastic parameter prediction model until a preset termination condition is met, and determining the elastic parameters of the pre-stack AVO inversion includes:

[0010] Determine whether the pre-stack AVO inversion objective function corresponding to the pre-stack AVO inversion iteration equation converges. If it does not converge, replace the elastic parameter prediction model of the work area with the elastic parameters determined by the pre-stack AVO inversion, and perform iterative joint inversion again until the pre-stack AVO inversion objective function converges, and determine the elastic parameters of the pre-stack AVO inversion under the convergence condition of the objective function.

[0011] In one embodiment, the pre-stack AVO reflection coefficient equation of the longitudinal P wave is:

[0012] R PP (θ P )≈A + B sin 2 θ P ;

[0013]

[0014] Δρ = ρ2 - ρ1,

[0015] ΔV P =V P2 -V P1 , ΔV S =V S2 -V S1 ;

[0016] Where: R pp (θ p ) is the reflection coefficient of the P wave, θ p is the average of the P wave incident angle and the P wave transmission angle, θ p1 and θ p2 are the incident angle and transmission angle of the P wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V p1 and V p2 are the longitudinal wave velocities of the media above and below the reflection interface respectively, V s1 and V s2 are the shear wave velocities of the media above and below the reflection interface respectively.

[0017] In one embodiment, the pre-stack AVO reflection coefficient equation of the pure shear SV wave is established by the following method:

[0018] Obtain the approximate expression of the reflection coefficient of the SV wave;

[0019] Simplify the approximate expression of the reflection coefficient of the SV wave to obtain the pre-stack AVO reflection coefficient equation of the SV wave;

[0020] The approximate expression of the reflection coefficient of the SV wave is:

[0021]

[0022]

[0023] Δρ = ρ2 - ρ1,

[0024]

[0025] where: R SV (θ s ) is the reflection coefficient of the SV wave, p is the component of the SV wave ray slowness parallel to the interface, θ s is the average of the SV wave incident angle and the SV wave transmitted angle, θ s1 and θ s2 are the incident angle and the transmitted angle of the SV wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the SV wave velocities of the media above and below the reflection interface respectively;

[0026] The pre-stack AVO reflection coefficient equation of the SV wave is:

[0027] R SV (θ S ) ≈ A + B sin 2 θ S ;

[0028]

[0029] In one embodiment, the pre-stack AVO reflection coefficient equation of the longitudinal wave SH wave is established by the following method:

[0030] Obtain the approximate expression of the reflection coefficient of the SH wave;

[0031] Simplify the approximate expression of the reflection coefficient of the SH wave to obtain the pre-stack AVO reflection coefficient equation of the SH wave;

[0032] The approximate expression of the reflection coefficient of the SH wave is:

[0033]

[0034] Δρ = ρ2 - ρ1,

[0035] ΔV S = V S2 - V S1 ,

[0036] where: R SH(θ s ) is the reflection coefficient of the SH wave, p is the component of the SH wave ray slowness parallel to the interface, θ s is the average of the incident angle and the transmitted angle of the SH wave, θ s1 and θ s2 are the incident angle and the transmitted angle of the SH wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the SH wave velocities of the media above and below the reflection interface respectively;

[0037] The pre-stack AVO reflection coefficient equation of the SH wave is:

[0038] R SH (θ S )≈A+B sin 2 θ S ;

[0039]

[0040] In one embodiment, extracting the angular gather seismic wavelets of the P wave, SV wave and SH wave includes:

[0041] Determining the logarithmic spectrum of the angular gather wavelet according to the actual angular gather seismic record;

[0042] The expression of the logarithmic spectrum of the angular gather wavelet is:

[0043]

[0044] Where: is the logarithmic spectrum of the wavelet, ω is the angular frequency, |X(ω)| is the amplitude spectrum of the seismic record, |W(ω)| is the amplitude spectrum of the wavelet, is the imaginary part of the logarithmic spectrum;

[0045] Writing the logarithmic spectrum sequence of the angular gather wavelet as the sum of the odd part and the even part sequence

[0046] For the Do Fourier transform to get

[0047] For Do inverse Fourier transform to obtain the angular gather seismic wavelet.

[0048] In one embodiment, substituting the elastic parameter prediction model into the pre-stack AVO reflection coefficient equation and convolving it with the angular gather seismic wavelet to obtain the angular gather synthetic record.

[0049] In one embodiment, establishing the P-wave and S-wave joint pre-stack AVO inversion iterative equation includes:

[0050] Perform Taylor expansion on the pre-stack AVO inversion objective function to obtain the Taylor expansion formula of the pre-stack AVO inversion objective function;

[0051] Find the extreme value of the Taylor expansion formula to obtain the pre-stack AVO inversion iteration equation.

[0052] In one embodiment, the pre-stack AVO inversion objective function is:

[0053]

[0054]

[0055] n - 1

[0056] f SH (Vs, ρ) = ||S SH (Vs, ρ) Δ - D|| = ∑(S SH (Vs, ρ) i Δ - D i ) 2 ;

[0057] i = 0

[0058] where: f PP (Vp, Vs, ρ) is the AVO inversion objective function of the PP longitudinal wave, f SV (Vs, ρ) is the AVO inversion objective function of the SV shear wave, f SH (Vs, ρ) is the AVO inversion objective function of the SH shear wave, is the PP longitudinal wave angle gather seismic record of the i-th layer of the formation, is the SV shear wave angle gather seismic record of the i-th layer of the formation, is the SH shear wave angle gather seismic record of the i-th layer of the formation, D i is the value of the i-th sampling point of the actual seismic angle gather record.

[0059] Second, the embodiment of the present invention provides a method for determining an oil and gas anomaly target area, including:

[0060] Determine the elastic parameters of the joint pre-stack AVO inversion of the seismic longitudinal and transverse waves;

[0061] Determine the longitudinal and transverse wave velocity ratio and Poisson's ratio according to the elastic parameters of the pre-stack AVO inversion;

[0062] Compare the longitudinal and transverse wave velocity ratio and the Poisson's ratio with a preset threshold, and determine the oil and gas anomaly target area according to the comparison result;

[0063] The steps of determining the elastic parameters of the prestack AVO inversion are obtained by using the joint P-wave and S-wave prestack AVO inversion method as described above.

[0064] In one embodiment, the P-wave to S-wave velocity ratio and Poisson's ratio are calculated according to the following formulas:

[0065]

[0066]

[0067] where: λ is the P-wave to S-wave velocity ratio, and ν is Poisson's ratio.

[0068] In a third aspect, an embodiment of the present invention provides a joint P-wave and S-wave prestack AVO inversion device, including:

[0069] An elastic parameter prediction model establishment module, configured to establish an elastic parameter prediction model for the work area according to well logging, drilling, vertical seismic profile (VSP) and geological data of the work area, where the elastic parameters include P-wave velocity, S-wave velocity and medium density;

[0070] An angle gather synthetic record establishment module, configured to establish an angle gather synthetic record according to the elastic parameter prediction model, the prestack AVO reflection coefficient equations of P-wave, pure S-wave SV and pure S-wave SH, and the angle gather seismic wavelet;

[0071] An inversion iteration equation establishment module, configured to establish a joint P-wave and S-wave prestack AVO inversion iteration equation according to the angle gather synthetic record and the actual angle gather record;

[0072] An inversion elastic parameter determination module, configured to perform joint inversion on the prestack AVO inversion iteration equation, repeatedly correct the elastic parameter prediction model until a preset termination condition is met, and determine the elastic parameters of the prestack AVO inversion.

[0073] In a fourth aspect, an embodiment of the present invention provides a computing device, including: a memory, a processor, and a computer program stored on the memory and executable on the processor, where when the processor executes the program, it implements the foregoing joint P-wave and S-wave prestack AVO inversion method or the foregoing method for determining an oil and gas anomaly target area.

[0074] In a fifth aspect, an embodiment of the present invention provides a computer-readable storage medium, where the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, it implements the foregoing joint P-wave and S-wave prestack AVO inversion method or the foregoing method for determining an oil and gas anomaly target area.

[0075] The beneficial effects of the above technical solutions provided by the embodiments of the present invention at least include:

[0076] The above-mentioned P-S wave joint pre-stack AVO inversion method and related devices provided by the embodiments of the present invention solve the problem of inaccurate pre-stack elastic parameter inversion in the prior art. In the prior art, elastic parameters, namely, P-wave velocity, S-wave velocity, and medium density, are mainly inverted based on the P-wave AVO equation. During the determination of the S-wave parameter term, the approximate P-S wave velocity ratio is directly assumed to be a fixed constant value. As a result, the inverted S-wave parameters, including S-wave velocity and medium density, are inaccurate. The embodiments of the present invention establish a complete set of elastic parameter inversion methods for P-wave, pure S-wave (SV-wave), and pure S-wave (SH-wave). An elastic parameter prediction model is established according to the actual working conditions, and then the AVO pre-stack elastic parameter inversion iterative equations for P-wave, pure S-wave (SV-wave), and pure S-wave (SH-wave) are established respectively. Through joint iterative inversion, the accurate P-S wave elastic parameters under actual working conditions are finally determined accurately. The elastic parameter inversion method provided by the embodiments of the present invention has accurate and reliable results, can effectively combine theory and practice (different working conditions), and has good applicability.

[0077] Furthermore, by adopting the seismic P-S wave joint pre-stack AVO inversion method in the embodiments of the present invention, not only can the P-S wave elastic parameters in different work areas be accurately determined, but also the key elastic parameters, namely, P-S wave velocity ratio and Poisson's ratio, can be further determined. It has higher accuracy in the application of directly locking hydrocarbon-bearing anomaly target areas, providing technical support for oil and gas exploration and development.

[0078] The technical solutions of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Description of the Drawings

[0079] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. They are used together with the embodiments of the present invention to explain the present invention, but do not constitute a limitation to the present invention. In the accompanying drawings:

[0080] Figure 1 is the flow chart of the P-S wave joint pre-stack AVO inversion method in the embodiments of the present invention;

[0081] Figure 2 is the schematic diagram of seismic wave incidence, reflection, and refraction in the embodiments of the present invention;

[0082] Figure 3 is the method for determining hydrocarbon-bearing anomaly target areas in the embodiments of the present invention;

[0083] Figure 4 is the flow schematic diagram of determining hydrocarbon-bearing anomaly target areas by P-S wave joint pre-stack AVO inversion in the embodiments of the present invention;

[0084] Figure 5This is the prediction map of the joint pre-stack AVO inversion of P-wave and S-wave and the oil and gas anomaly target area in a certain work area in the embodiment of the present invention;

[0085] Figure 6 This is the device for joint pre-stack AVO inversion of P-wave and S-wave in the embodiment of the present invention. Specific implementation manners

[0086] The embodiment of the present invention provides a method for joint pre-stack AVO inversion of P-wave and S-wave. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be fully conveyed to those skilled in the art. Specifically, it includes the following steps:

[0087] The embodiment of the present invention provides a method for joint pre-stack AVO inversion of P-wave and S-wave. Referring to Figure 1 as shown, this method includes the following steps:

[0088] S11. According to the logging, drilling, vertical seismic profile (VSP) and geological data of the work area, establish an elastic parameter prediction model for the work area. The elastic parameters include the P-wave velocity, S-wave velocity and medium density;

[0089] S12. According to the elastic parameter prediction model, the pre-stack AVO reflection coefficient equations of the P-wave, pure SV-wave and pure SH-wave, and the angle gather seismic wavelet, establish an angle gather synthetic record;

[0090] S13. According to the angle gather synthetic record and the actual angle gather record, establish a joint pre-stack AVO inversion iteration equation;

[0091] S14. Perform joint inversion on the pre-stack AVO inversion iteration equation, and repeatedly correct the elastic parameter prediction model until the preset termination condition is met, and determine the elastic parameters of the pre-stack AVO inversion.

[0092] In step S11, collect and sort out the existing logging, drilling and vertical seismic profile (VSP) data of the work area, and obtain the elastic parameters of each single well in the work area, including: P-wave velocity, S-wave velocity and medium density. The interpolation method is used to determine the elastic parameter values within the entire work area, and then an elastic prediction model is established, that is, a spatial function of the P-wave velocity, S-wave velocity and medium density.

[0093] In step S12, in order to establish an angle gather synthetic record, it is necessary to establish the pre-stack AVO reflection coefficient equations of the P-wave, pure SV-wave and pure SH-wave respectively. The specific implementation steps are as follows:

[0094] (1) Establish the pre-stack AVO reflection coefficient equation of the P-wave.

[0095] Refer to Figure 2 As shown, according to the Zoeppritz equation, assuming a small perturbation layered medium, a simplified formula for the pre-stack AVO reflection coefficient of P-wave with respect to the P-wave velocity, S-wave velocity, and medium density can be obtained. The expression for the pre-stack AVO reflection coefficient of P-wave is as follows:

[0096] R PP (θ P )≈A + Bsin 2 θ P ;

[0097]

[0098] Δρ = ρ2 - ρ1,

[0099] ΔV P = V P2 - V P1 , ΔV S = V S2 - V S1 ;

[0100] Where: R pp (θ p ) is the reflection coefficient of P-wave, θ p is the average of the incident angle and transmission angle of P-wave, θ p1 and θ p2 are the incident angle and transmission angle of P-wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V p1 and V p2 are the P-wave velocities of the media above and below the reflection interface respectively, V s1 and V s2 are the S-wave velocities of the media above and below the reflection interface respectively.

[0101] (2) Establish the pre-stack AVO reflection coefficient equation for pure shear wave SV-wave.

[0102] Refer to Figure 2 As shown, according to the Rüger equation, assuming a small perturbation layered medium and finite wave angles, a simplified formula for the pre-stack AVO reflection coefficient of SV-wave can be obtained, which can be expressed in terms of the pure shear wave SV-wave velocity and medium density. The Rüger equation for the pre-stack AVO reflection coefficient of SV-wave is as follows:

[0103]

[0104]

[0105] Δρ = ρ2 - ρ1,

[0106]

[0107] where: R SV (θ s ) is the reflection coefficient of the SV wave, p is the component of the SV wave ray slowness parallel to the interface, θ s is the average value of the SV wave incident angle and the SV wave transmission angle, θ s1 and θ s2 are the incident angle and the transmission angle of the SV wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the SV wave velocities of the media above and below the reflection interface respectively;

[0108] Simplify the Rüger equation for the SV wave to obtain the deformed form of the SV wave reflection coefficient. The approximate deformed form of the prestack AVO reflection coefficient of the SV wave is:

[0109]

[0110] Assume θ s < 35°, omit the monomial (the third term) that approaches zero in the above deformed form of the SV wave reflection coefficient to obtain the prestack AVO reflection coefficient equation of the SV wave. The prestack AVO reflection coefficient equation of the SV wave is as follows:

[0111] R SV (θ S ) ≈ A + Bsin 2 θ S ;

[0112]

[0113] where: A is the intercept of the prestack AVO reflection coefficient equation of the SV wave, that is, when the angle is zero, the zero-angle reflection coefficient is the reflection coefficient of the zero-offset profile, and B is the gradient (slope) of the prestack AVO reflection coefficient equation of the SV wave.

[0114] (3) Establish the prestack AVO reflection coefficient equation of the pure shear wave SH wave.

[0115] Referring to Figure 2 as shown, Rüger (2001) obtained the simplified formula of the prestack AVO reflection coefficient of the SH wave by assuming a small perturbation layered medium and a finite wave angle. The Rüger equation of the prestack AVO reflection coefficient of the SH wave is as follows:

[0116]

[0117] Δρ = ρ2 - ρ1,

[0118]

[0119] where: R SH (θ s ) is the reflection coefficient of the SH wave, p is the component of the slowness of the SH wave ray parallel to the interface, θ s is the average of the incident angle and the transmitted angle of the SH wave, θ s1 and θ s2 are the incident angle and the transmitted angle of the SH wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the shear wave SH velocities of the media above and below the reflection interface respectively.

[0120] Assume θ s < 35°, omit the monomial (the third term) that approaches zero in the above SH wave Rüger equation, and obtain the pre-stack AVO reflection coefficient equation of the SH wave. The pre-stack AVO reflection coefficient equation of the SH wave is as follows:

[0121] The pre-stack AVO reflection coefficient equation of the SH wave is:

[0122] R SH (θ S ) ≈ A + Bsin 2 θ S ;

[0123]

[0124] where: A is the intercept of the pre-stack AVO reflection coefficient equation of the SH wave, that is, when the angle is zero, the zero-angle reflection coefficient is the reflection coefficient of the zero-offset profile, and B is the gradient (slope) of the pre-stack AVO reflection coefficient equation of the SH wave.

[0125] In step S12, in order to establish the angle gather synthesis record, it is necessary to extract the angle gather seismic wavelet.

[0126] Assume that the reflection system is a random white noise sequence. Since the autocorrelation of the seismic record is equal to the autocorrelation of the wavelet, that is, the amplitude spectrum |X(ω)| of the seismic record is equal to the amplitude spectrum |W(ω)| of the wavelet, and their logarithmic spectra are also equal, ln|X(ω)| = ln|W(ω)|, as follows:

[0127] X(t) = W(t) * R(t)

[0128] X(ω) = W(ω) × R(ω)

[0129] ln X(ω) = ln W(ω) + ln R(ω)

[0130] Determining the logarithmic spectrum of a wavelet from the geometric mean of several amplitude spectra:

[0131]

[0132] Wherein, is the logarithmic spectrum of the wavelet, ω is the angular frequency, |X(ω)| is the amplitude spectrum of the seismic record, |W(ω)| is the amplitude spectrum of the wavelet, is the imaginary part of the logarithmic spectrum and also the phase of the corresponding wavelet logarithmic spectrum, where |X(ω)| can be directly extracted from the common angle gather record.

[0133] When the wavelet is minimum phase, its logarithmic spectrum sequence is a real causal sequence, and any real sequence can be written as the sum of an odd part and an even part sequence, as follows:

[0134]

[0135] For perform Fourier transform to obtain

[0136] For perform inverse Fourier transform to obtain the common angle gather seismic wavelet W(t), and the common angle gather seismic wavelet is discretely represented as:

[0137] W = {0, W 2-m , … W -1 , W0, W1 …, W m-2 , 0};

[0138] Wherein: W is the symmetric discretized wavelet array, the length of the wavelet array is 2m + 1, W0 is the extreme value term of the wavelet, W m-2 and W 2-m are equal, and m is the number of sample points on the positive or negative half axis of the wavelet signal.

[0139] According to the pre-stack AVO reflection coefficient equations of the longitudinal P-wave, pure transverse SV-wave and pure transverse SH-wave established above, the common angle gather seismic wavelet extracted, and the elastic parameter prediction model established in step S11, using the convolution model, the common angle gather synthetic seismogram can be established. The specific method is as follows:

[0140] First, represent the pre-stack AVO reflection coefficients of the P-wave, SV-wave and SH-wave as abstract functions of the wave velocity v, medium density ρ, and ray angle θ, and the respective abstract functions are:

[0141] R PP (θ) i = R1(Vp i , Vs i , ρ i , Vp i+1 , Vs i+1, ρ i+1 );

[0142] R SV (θ) i = R2(Vs i , ρ i , Vs i+1 , ρ i+1 );

[0143] R SH (θ) i = R3(Vs i , ρ i , Vs i+1 , ρ i+1 );

[0144] Where: R PP (θ) i is the PP longitudinal wave reflection coefficient of the i-th layer of the formation, R SV (θ) i is the SV shear wave reflection coefficient of the i-th layer of the formation, R SH (θ) i is the SH shear wave reflection coefficient of the i-th layer of the formation, Vp i is the PP longitudinal wave velocity of the i-th layer of the formation, Vs i is the shear wave velocity of the i-th layer of the formation, ρ i is the density of the i-th layer of the formation.

[0145] Substitute each elastic parameter prediction model in step S11 into the above abstract functions, and convolve with the angle gather seismic wavelet extracted above to obtain the angle gather synthetic record. The angle gather synthetic record is as follows:

[0146]

[0147]

[0148]

[0149] Where: S is the seismic record, R is the formation reflection coefficient, W is the angle gather seismic wavelet, S PP (θ) i is the PP longitudinal wave angle gather seismic record of the i-th layer of the formation, S SV (θ) i is the SV shear wave angle gather seismic record of the i-th layer of the formation, S SH (θ) i is the SH shear wave angle gather seismic record of the i-th layer of the formation, R PP (θ) k is the PP longitudinal wave reflection coefficient of the k-th layer of the formation, R SV (θ) kThe SV shear wave reflection coefficient of the k-th layer of the formation, R SH (θ) k is the SH shear wave reflection coefficient of the k-th layer of the formation.

[0150] In step S13, the pre-stack migration gather data is preprocessed and converted into an angle gather record, that is, the actual angle gather record is obtained; the synthetic angle gather record has been obtained in the above step S12.

[0151] The actual angle gather seismic record is the seismic record of the average value of the pre-stack migration and stacking angles of the actual seismic data. The actual angle gather record is discretely expressed as:

[0152] D(θ) = {D(θ)0, D(θ)1 …… D(θ) n-1};

[0153] where: D(θ) is the seismic record with the average value of the actual seismic stacking angle being θ, θ is the angle, and D(θ) i is the value of the i-th formation sampling point of the seismic record corresponding to the θ angle, i = 0, 1, 2 …… n - 1, and n is the maximum number of sampling points.

[0154] The angle gather seismic wavelet is discretely expressed as:

[0155] W = {0, W 2-m , … W -1 , W0, W1 …, W m-2 , 0};

[0156] The wave velocity of the longitudinal P wave is discretely expressed as:

[0157] Vp = {Vp0, Vp1, ……, Vp n-1 , Vp n};

[0158] The wave velocity of the shear S wave is discretely expressed as:

[0159] Vs = {Vs0, Vs1, ……, Vs n-1 , Vs n};

[0160] The medium density ρ is discretely expressed as:

[0161] ρ = {ρ0, ρ1, ……, ρ n-1 , ρ n};

[0162] In order to establish the joint pre-stack AVO inversion iterative equation for P and S waves, first establish the pre-stack AVO inversion objective functions for P wave, SV wave and SH wave.

[0163] The pre-stack AVO inversion objective functions for different seismic waves are:

[0164]

[0165]

[0166]

[0167] where: f PP (Vp, Vs, ρ) is the AVO inversion objective function for PP longitudinal waves, f SV (Vs, ρ) is the AVO inversion objective function for SV shear waves, f SH (Vs, ρ) is the AVO inversion objective function for SH shear waves, is the PP longitudinal wave angle gather seismic record of the i-th layer of the formation, is the SV shear wave angle gather seismic record of the i-th layer of the formation, is the SH shear wave angle gather seismic record of the i-th layer of the formation, D i is the value of the i-th sampling point of the actual seismic angle gather record.

[0168] Then, perform a Taylor expansion on the pre-stack AVO inversion objective function established above to transform the non-linear objective function into a linear function. The linear objective functions corresponding to different seismic waves are:

[0169] The linear objective function for pre-stack AVO inversion of PP waves is:

[0170]

[0171] The linear objective function for pre-stack AVO inversion of SV waves is:

[0172]

[0173] The linear objective function for pre-stack AVO inversion of SH waves is:

[0174]

[0175] Next, find the extreme values of the linear objective functions for pre-stack AVO inversion corresponding to different seismic waves above, that is, take the derivatives of ΔVp, ΔVs, and Δρ respectively at the same time to obtain the corresponding pre-stack AVO inversion iteration equations. The pre-stack AVO inversion iteration equations corresponding to different seismic waves are:

[0176] The inversion iteration equation corresponding to the longitudinal P wave is:

[0177]

[0178] The inversion iteration equation corresponding to the pure shear SV wave is:

[0179]

[0180] The inversion iteration equation corresponding to the pure shear wave SH wave is as follows:

[0181]

[0182] Furthermore, rewrite the inversion iteration equations corresponding to the above different seismic waves into matrix form, that is, obtain the inversion iteration matrix equation. The pre-stack AVO inversion iteration matrix equations corresponding to different seismic waves are as follows:

[0183] The inversion iteration matrix equation corresponding to the longitudinal wave P wave is:

[0184]

[0185] The inversion iteration matrix equation corresponding to the pure shear wave SV wave is:

[0186]

[0187] The inversion iteration matrix equation corresponding to the pure shear wave SH wave is:

[0188]

[0189] Where: is the partial derivative term of the longitudinal wave velocity V PP of the linear function after the Taylor expansion of S p (θ), is the partial derivative of the shear wave velocity V PP of the linear function after the Taylor expansion of S s (θ), is the partial derivative of the density ρ PP of the linear function after the Taylor expansion of S is the partial derivative of the shear wave velocity V SV of the linear function after the Taylor expansion of S s (θ), is the partial derivative of the medium density ρ SV of the linear function after the Taylor expansion of S is the partial derivative of the shear wave velocity V SH of the linear function after the Taylor expansion of S s (θ), is the partial derivative of the medium density ρ SH of the linear function after the Taylor expansion of S

[0190] In step S14, it is necessary to perform joint inversion on the established pre-stack AVO inversion iteration equation to determine the elastic parameters of the pre-stack AVO inversion.

[0191] The inventor of the present invention found that there is a problem with the pre-stack inversion formula of the longitudinal wave P wave in the prior art Item, it is necessary to approximate this item The accuracy of the shear wave velocity and medium density calculated in this way is relatively low.

[0192] In this step, first, solve the established pre-stack AVO inversion iterative equation to obtain a set of elastic parameters for pre-stack AVO inversion;

[0193] Then, determine whether the pre-stack AVO inversion objective function converges. For example: Substitute the predicted values of the elastic parameters in the elastic parameter prediction model of the work area established in step S11 into the pre-stack AVO inversion objective function established in step S13, and determine whether the value of the above pre-stack AVO inversion objective function is zero. If the result is zero, it converges; if it is not zero, it does not converge;

[0194] Next, if the above result does not converge, replace the elastic parameter prediction model of the work area with the elastic parameters determined by the pre-stack AVO inversion iterative equation, and perform iterative joint inversion again until the pre-stack AVO inversion objective function converges, and determine the elastic parameters of the pre-stack AVO inversion under the convergence condition of the objective function. The process of solving the pre-stack AVO inversion iterative equation is also the process in which the inversion objective function gradually converges to zero.

[0195] The embodiment of the present invention also provides a method for determining an oil and gas anomaly target area. Refer to Figure 3 As shown, this method includes the steps:

[0196] S31. Determine the elastic parameters of seismic P-S wave joint pre-stack AVO inversion;

[0197] S32. Determine the P-S wave velocity ratio and Poisson's ratio according to the elastic parameters of pre-stack AVO inversion;

[0198] S33. Compare the P-S wave velocity ratio and Poisson's ratio with a preset threshold, and determine the oil and gas anomaly target area according to the comparison result.

[0199] In step S31, the steps of determining the elastic parameters of pre-stack AVO inversion are obtained by using the methods of the foregoing steps S11-S14.

[0200] In step S32, the P-S wave velocity ratio and Poisson's ratio are calculated according to the following formula:

[0201]

[0202]

[0203] Where: λ is the P-S wave velocity ratio, and ν is Poisson's ratio.

[0204] Since there are two types of shear waves, including pure shear wave SV wave and pure shear wave SH wave, two sets of P-S wave velocity ratios and Poisson's ratios in different shear wave directions can be obtained respectively. In specific engineering applications, different parameter values should be selected according to the actual working conditions and experience, and the embodiments of the present invention do not limit this.

[0205] In step S33, thresholds can be preset for the P-S wave velocity ratio and Poisson's ratio respectively. It is required that both the P-S wave velocity ratio and Poisson's ratio are less than their respective thresholds, or it is required that one of the P-S wave velocity ratio or Poisson's ratio is less than a specific threshold to determine the hydrocarbon-bearing abnormal area in the work area. The setting method of the thresholds in the embodiments of the present invention is not limited.

[0206] In summary, combining the P-S wave joint pre-stack AVO inversion method provided by the embodiments of the present invention and the method for determining the hydrocarbon-bearing abnormal target area, a set of processes for determining the hydrocarbon-bearing abnormal target area by P-S wave joint pre-stack AVO inversion is obtained. Refer to Figure 4 As shown, first establish a prediction model of elastic parameters (V p 、V s and ρ1) under actual working conditions, and calculate the reflection coefficient under the elastic parameter prediction model; then extract the angle gather wavelet and use the convolution model to determine the angle gather synthetic record; establish an inversion iteration equation according to the foregoing work for joint inversion to determine the elastic parameters of pre-stack AVO inversion; determine the P-S wave velocity ratio and Poisson's ratio, and identify the reservoir and determine the hydrocarbon-bearing abnormal target area.

[0207] Taking a regional project as an example, according to the above process method, refer to Figure 5 As shown, through P-S wave joint pre-stack AVO inversion, the elastic parameters are inverted: the P-wave velocity and S-wave velocity. Further, the key references are determined: the P-S wave velocity ratio and Poisson's ratio. According to experience, the area where Poisson's ratio in the target formation is less than 0.22 is determined as the gas-bearing layer. It can be easily seen from Figure 5 that there are 5 gas-bearing layers that meet the requirements in this work area.

[0208] The embodiments of the present invention provide a P-S wave joint pre-stack AVO inversion device. Refer to Figure 6 As shown, it includes:

[0209] An elastic parameter prediction model establishment module 61, configured to establish an elastic parameter prediction model for the work area according to well logging, drilling, vertical seismic profile VSP and geological data of the work area. The elastic parameters include P-wave velocity of longitudinal wave, S-wave velocity of transverse wave and medium density;

[0210] An angle gather synthetic record establishment module 62, configured to establish an angle gather synthetic record according to the elastic parameter prediction model, the pre-stack AVO reflection coefficient equations of longitudinal wave P-wave, pure shear wave SV wave and pure shear wave SH wave, and the angle gather seismic wavelet;

[0211] An inversion iteration equation establishing module 63, configured to establish a combined P-wave and S-wave pre-stack AVO inversion iteration equation according to the synthetic record of the angle gather and the actual record of the angle gather;

[0212] An inversion elastic parameter determining module 64, configured to perform combined inversion on the pre-stack AVO inversion iteration equation, repeatedly correct the elastic parameter prediction model until a preset termination condition is met, and determine the elastic parameters of the pre-stack AVO inversion.

[0213] An embodiment of the present invention provides a computing device, including: a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the foregoing combined P-wave and S-wave pre-stack AVO inversion method or the foregoing method for determining an oil and gas anomaly target area is implemented.

[0214] An embodiment of the present invention provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the foregoing combined P-wave and S-wave pre-stack AVO inversion method or the foregoing method for determining an oil and gas anomaly target area is implemented.

[0215] Obviously, those skilled in the art can make various changes to the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications of the present invention fall within the scope of the claims of the present invention and their equivalent technologies, the present invention also intends to include these modifications.

Claims

1. A combined pre-stack AVO inversion method for P-wave and S-wave of earthquake, characterized in that, Including: Based on the logging, drilling, vertical seismic profile (VSP) and geological data of the work area, an elastic parameter prediction model for the work area is established, and the elastic parameters include the longitudinal wave P-wave velocity, the transverse wave S-wave velocity, and the medium density; Based on the elastic parameter prediction model, the pre-stack AVO reflection coefficient equations of the longitudinal wave P-wave, the pure transverse wave SV-wave, and the pure transverse wave SH-wave, and the angle gather seismic wavelet, an angle gather synthetic record is established; Based on the angle gather synthetic record and the actual angle gather record, a joint pre-stack AVO inversion iteration equation for P-wave and S-wave is established; Perform joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correct the elastic parameter prediction model until a preset termination condition is met, and determine the elastic parameters of the pre-stack AVO inversion.

2. The method according to claim 1, characterized in that, The performing joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correcting the elastic parameter prediction model until a preset termination condition is met, and determining the elastic parameters of the pre-stack AVO inversion includes: Judge whether the pre-stack AVO inversion objective function corresponding to the pre-stack AVO inversion iteration equation converges. If it does not converge, replace the elastic parameter prediction model of the work area with the determined elastic parameters of the pre-stack AVO inversion, and perform iterative joint inversion again until the pre-stack AVO inversion objective function converges, and determine the elastic parameters of the pre-stack AVO inversion under the condition of the convergence of the objective function.

3. The method according to claim 1, characterized in that, The pre-stack AVO reflection coefficient equation of the longitudinal wave P-wave is: R PP (θ P )≈A + Bsin 2 θ P ; Δρ = ρ2 - ρ1, Where: R pp (θ p ) is the reflection coefficient of the P-wave, θ p is the average of the incident angle and the transmitted angle of the P-wave, θ p1 and θ p2 are the incident angle and the transmitted angle of the P-wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V p1 and V p2 are the longitudinal wave velocities of the media above and below the reflection interface respectively, V s1 and V s2 are the shear wave velocities of the media above and below the reflection interface respectively.

4. The method according to claim 1, wherein The pre-stack AVO reflection coefficient equation of the pure transverse wave SV-wave is established by the following method: Obtain the approximate expression of the reflection coefficient of the SV-wave; Simplify the approximate expression of the reflection coefficient of the SV-wave to obtain the pre-stack AVO reflection coefficient equation of the SV-wave; The approximate expression of the reflection coefficient of the SV-wave is: Δρ = ρ2 - ρ1, Where: R SV (θ s ) is the reflection coefficient of the SV wave, p is the component of the SV wave ray slowness parallel to the interface, θ s is the average of the incident angle and the transmitted angle of the SV wave, θ s1 and θ s2 are the incident angle and the transmitted angle of the SV wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the SV wave velocities of the media above and below the reflection interface respectively; The pre-stack AVO reflection coefficient equation of the SV-wave is: R SV (θ S )≈A + Bsin 2 θ S ; 5. The method according to claim 1, characterized in that, The pre-stack AVO reflection coefficient equation of the longitudinal wave SH-wave is established by the following method: Obtain the approximate expression of the reflection coefficient of the SH-wave; Simplify the approximate expression of the reflection coefficient of the SH-wave to obtain the pre-stack AVO reflection coefficient equation of the SH-wave; The approximate expression of the reflection coefficient of the SH-wave is: Δρ = ρ2 - ρ1, Where: R SH (θ s ) is the reflection coefficient of the SH wave, p is the component of the SH wave ray slowness parallel to the interface, θ s is the average of the incident angle and the transmitted angle of the SH wave, θ s1 and θ s2 are the incident angle and the transmitted angle of the SH wave respectively, ρ1 and ρ2 are the densities of the media above and below the reflection interface respectively, V s1 and V s2 are the SH wave velocities of the media above and below the reflection interface respectively; The pre-stack AVO reflection coefficient equation of the SH-wave is: R SH (θ S )≈A + Bsin 2 θ S ; 6. The method according to claim 1, wherein The extracting the angle gather seismic wavelets of the P-wave, SV-wave, and SH-wave includes: Based on the actual angle gather seismic record, determine the logarithmic spectrum of the angle gather wavelet; The logarithmic spectrum expression of the angle gather wavelet is: Wherein: is the wavelet logarithmic spectrum, ω is the angular frequency, |X(ω)| is the amplitude spectrum of the seismic record, and |W(ω)| is the amplitude spectrum of the wavelet, is the imaginary part of the logarithmic spectrum; Write the logarithmic spectrum sequence of the corner gather wavelet as the sum of the sequences of the odd part and the even part Perform Fourier transform on the to obtain Perform an inverse Fourier transform on to obtain the angular gather seismic wavelet.

7. The method according to claim 1, characterized in that, Substitute the elastic parameter prediction model into the pre-stack AVO reflection coefficient equation and perform convolution with the angle gather seismic wavelet to obtain the angle gather synthetic record.

8. The method according to claim 7, characterized in that, The establishing the joint pre-stack AVO inversion iteration equation for P-wave and S-wave includes: Perform Taylor expansion on the pre-stack AVO inversion objective function to obtain the Taylor expansion of the pre-stack AVO inversion objective function; Find the extreme value of the Taylor expansion to obtain the pre-stack AVO inversion iteration equation.

9. The method according to claim 8, characterized in that The pre-stack AVO inversion objective function is: where: f PP (Vp, Vs, ρ) is the AVO inversion objective function for PP longitudinal waves, f SV (Vs, ρ) is the AVO inversion objective function for SV shear waves, f SH (Vs, ρ) is the AVO inversion objective function for SH shear waves, is the PP longitudinal wave common-azimuth seismic record of the i-th layer of the formation, is the SV shear wave common-azimuth seismic record of the i-th layer of the formation, is the SH shear wave common-azimuth seismic record of the i-th layer of the formation, D i is the value of the i-th sampling point of the actual seismic common-azimuth record.

10. A method for determining an abnormal target area containing oil and gas, characterized in that, Including: Determine the elastic parameters of the joint pre-stack AVO inversion of P-wave and S-wave in the seismic data; Determine the P-wave to S-wave velocity ratio and Poisson's ratio based on the elastic parameters obtained from pre-stack AVO inversion; Compare the P-wave to S-wave velocity ratio and Poisson's ratio with preset thresholds, and determine the hydrocarbon-bearing anomaly target area according to the comparison results; The steps of determining the elastic parameters of the pre-stack AVO inversion are obtained by using the method described in claims 1-9.

11. The method according to claim 10, characterized in that, The P-wave to S-wave velocity ratio and Poisson's ratio are calculated according to the following formula: Where: λ is the P-wave to S-wave velocity ratio, and ν is Poisson's ratio.

12. A combined P-wave and S-wave pre-stack AVO inversion device, characterized in that, It includes: An elastic parameter prediction model establishment module, configured to establish an elastic parameter prediction model for the work area according to well logging, drilling, vertical seismic profile (VSP) and geological data of the work area, and the elastic parameters include P-wave velocity, S-wave velocity and medium density; An angle gather composite record establishment module, configured to establish an angle gather composite record according to the elastic parameter prediction model, the pre-stack AVO reflection coefficient equations of P-wave, pure S-wave (SV-wave) and pure S-wave (SH-wave), and the angle gather seismic wavelet; An inversion iteration equation establishment module, configured to establish a joint pre-stack AVO inversion iteration equation for P-wave and S-wave according to the angle gather composite record and the actual angle gather record; An inversion elastic parameter determination module, configured to perform joint inversion on the pre-stack AVO inversion iteration equation, repeatedly correct the elastic parameter prediction model until a preset termination condition is met, and determine the elastic parameters of the pre-stack AVO inversion.

13. A computing device, characterized in that, It includes: A memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, it implements the joint pre-stack AVO inversion method for P-wave and S-wave described in any one of claims 1-9, or the method for determining the hydrocarbon-bearing anomaly target area described in any one of claims 10-11.

14. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program, and when the computer program is executed by the processor, it implements the joint pre-stack AVO inversion method for P-wave and S-wave described in any one of claims 1-9, or the method for determining the hydrocarbon-bearing anomaly target area described in any one of claims 10-11.

Images