CN103728662A - Method for estimating stratum medium quality factors based on seismic signal envelope peak - Google Patents

Method for estimating stratum medium quality factors based on seismic signal envelope peak Download PDF

Info

Publication number
CN103728662A
CN103728662A CN201410003193.6A CN201410003193A CN103728662A CN 103728662 A CN103728662 A CN 103728662A CN 201410003193 A CN201410003193 A CN 201410003193A CN 103728662 A CN103728662 A CN 103728662A
Authority
CN
China
Prior art keywords
mrow
msub
msubsup
mfrac
seismic
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201410003193.6A
Other languages
Chinese (zh)
Inventor
赵伟
高静怀
赵静
杨森林
姜秀娣
翁斌
朱振宇
陈剑军
印海燕
王清振
刘志鹏
丁继才
江南森
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Offshore Oil Corp CNOOC
Xian Jiaotong University
CNOOC Research Institute Co Ltd
Original Assignee
China National Offshore Oil Corp CNOOC
Xian Jiaotong University
CNOOC Research Institute Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China National Offshore Oil Corp CNOOC, Xian Jiaotong University, CNOOC Research Institute Co Ltd filed Critical China National Offshore Oil Corp CNOOC
Priority to CN201410003193.6A priority Critical patent/CN103728662A/en
Publication of CN103728662A publication Critical patent/CN103728662A/en
Pending legal-status Critical Current

Links

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to a method for estimating stratum medium quality factors based on a seismic signal envelope peak. The method comprises the steps that 1) media are divided into a plurality of small sheets with the distance between two adjacent detectors as a thickness; 2) constant phase wavelets are used for approaching direct waves received by the top of the ith small sheet, and seismic wavelet parameters are estimated through a high-order cumulant matching method based on MSMG; 3) the parameters of a wavelet transform kernel function are estimated through the high-order cumulant matching method based on the MSMG; 4) the signal receiving instantaneous frequency of the upper interface and the lower interface of the ith small sheet is calculated; 5) the changes of envelope peak instantaneous frequency in the ith small sheet are calculated; 6) the quality factor Q value of the ith small sheet is calculated; 7) the steps 2)-6) are repeated, the quality factors of the other N-1 small sheets except the ith small sheet are calculated in sequence; 8) the oil-gas possibility of a reservoir stratum is forecast through the quality factor Q values, obtained in the step 6) and the step 7), of thee N small sheets, and then an estimated attenuation curve is obtained.

Description

Stratum medium quality factor estimation method based on seismic signal envelope peak
Technical Field
The invention relates to the field of seismic surveying, in particular to a method for estimating a stratum medium quality factor based on an envelope peak of a seismic signal.
Background
During the propagation of seismic waves in the subsurface, the seismic waves are absorbed due to the viscoelasticity of the formation, resulting in amplitude attenuation and velocity dispersion. The magnitude of seismic attenuation is typically measured by the quality factor Q of the formation. Laboratory and actual data measurements show that the quality factor Q is related to rock properties, fluid saturation and other factors. Therefore, the quality factor Q value is an effective tool for reservoir identification and hydrocarbon detection. In addition, the Q value of the quality factor has important significance in better explaining the effect of AVO (Amplitude Versus Offset, the variation of Amplitude along with Offset), improving the seismic imaging resolution, detecting and monitoring fluid reservoirs in time-lapse seismic, improving the research of formation physical properties and the like.
Various methods have been proposed for formation attenuation parameter estimation. The quality factor Q value is typically estimated using the amplitude of the seismic signal. In the time domain, the quality factor Q is generally calculated by using the pulse amplitude attenuation, pulse rise time, pulse stretching and the like. These methods all require the use of pulse amplitudes, but since the amplitude information of seismic pulses is often affected by scattering, geometric diffusion and other factors, the accuracy of the quality factor Q values estimated by these time domain methods is low. Attenuation estimation methods in the frequency domain typically include the log-ratio method (LSR), center frequency shift method (CFS), and peak frequency shift method. Firstly, intercepting a section of seismic record by using a time window, and then calculating a Fourier spectrum of the intercepted seismic record; however, in practical use, the spectrum estimation may be inaccurate once the type and length of the time window are not properly selected, and the attenuation estimation accuracy is affected. In the prior art, a method for estimating the quality factor Q value by using the change of the peak scale in a wavelet domain is provided by assuming that a source wavelet is an ideal pulse, but the actual source wavelet and pulse signals have large difference, so the method is limited in practical application.
The quality factor Q value may also be estimated by the instantaneous frequency variation of the seismic wavelets. Gabor teaches the concept of instantaneous frequency, which Taner uses for seismic interpretation. Tonn, Barnes et al give the relationship between different seismic instantaneous spectral measurements and seismic attenuation, respectively. Barnes gives a relation between instantaneous frequency and quality factor Q value and transmission time on the assumption that the power spectrum of the source wavelet is an ideal band-pass wavelet, but the method needs to be improved because the power spectrum of the actual source is greatly different from the ideal band-pass wavelet. Mathney and Nowack propose an instantaneous frequency matching method, namely an iterative process is adopted to match a reference pulse acted by a causal attenuation operator and a weighted instantaneous frequency at a target pulse envelope peak value, and the attenuation of shot gather data is estimated by the method; dasios et al estimated the attenuation of full-wave sonic logs using an instantaneous frequency matching method. This instantaneous frequency matching method overcomes some of the disadvantages of the log-ratio method, such as the need to select a variable frequency band range, but it requires the computation of the instantaneous frequency using the Hilbert transform and a complex iterative process to match the instantaneous frequency. It is well known that the Hilbert transform is susceptible to noise, and therefore the application of the instantaneous frequency matching method to noisy seismic signals is limited. High quiet eye et al propose a WEPIF (instantaneous frequency at wavelet envelope peak) method that calculates the instantaneous frequency in the wavelet domain and estimates the attenuation using the variation of the instantaneous frequency. However, this method still does not solve the following problems: 1. how to estimate wavelet parameters with high precision in WEPIF; 2. and (3) optimal selection of parameters in the three-parameter wavelet.
Disclosure of Invention
In view of the above problems, an object of the present invention is to provide a medium quality factor estimation method based on seismic signal envelope peaks, which can effectively improve quality factor estimation accuracy and has better random noise resistance.
In order to achieve the purpose, the invention adopts the following technical scheme: a stratum medium quality factor estimation method based on seismic signal envelope peaks comprises the following steps: 1) dividing the medium into a plurality of small thin plates by taking the distance between two adjacent detectors in the well as the thickness; 2) for the ith small thin plate, a common phase wavelet is used for approaching a direct wave received at the top of the small thin plate, and the seismic wavelet parameters are estimated by a high-order cumulant matching method based on MSMG; wherein MSMG is a multivariate Gaussian mixture model; the method for approaching the top of the small thin plate to receive the direct wave by using the constant phase wavelet comprises the following steps: (1) in the horizontal layered viscoelastic medium, the quality factor Q value of each layer is taken as a constant, and only the single-pass wave propagation of the plane wave is considered, so that the frequency expression of the seismic source wavelet at the earth surface when the seismic source wavelet reaches the depth z is as follows:
<math> <mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>)</mo> </mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mo>[</mo> <mo>-</mo> <mfrac> <mi>i&omega;&Delta;z</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mfrac> <mi>&omega;&Delta;z</mi> <mrow> <mn>2</mn> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>Q</mi> </mrow> </mfrac> <mo>]</mo> <mo>,</mo> </mrow> </math>
where Δ z is the transmission distance of the source wavelet, ω is the angular frequency of the source wavelet, G (Δ z) is a factor independent of frequency and absorption,frequency domain expression for source waveletExp is an exponential function based on the natural logarithm e,
Figure BDA0000452740090000023
c (omega) is the phase velocity of the seismic source wavelet; (2) adopting a constant phase wavelet to approximate a seismic source wavelet, wherein the function expression of the constant phase wavelet is as follows:
Figure BDA0000452740090000024
wherein,
Figure BDA0000452740090000025
the phase constant of the seismic source wavelet is used, sigma is a scale parameter of the seismic source wavelet, the width of the wavelet is controlled, delta is a frequency parameter of the seismic source wavelet, and t is the propagation time of the seismic source wavelet; the process of estimating the parameters of the seismic wavelets by using the MSMG-based high-order cumulant matching method comprises the following steps: (1) estimation of phase parameters of seismic wavelets using existing phase rotation methods(2) Estimating a scale parameter sigma and a frequency parameter delta by using an cumulant matching method; 3) estimating parameters of a kernel function of wavelet transformation by a high-order cumulant matching method based on MSMG; 4) calculating the instantaneous frequency of the received signal at the upper and lower interfaces of the ith small thin plate; 5) calculating the change of the instantaneous frequency of the envelope peak value in the ith small thin plate; 6) calculating the Q value of the quality factor of the ith small sheet: three parameters sigma of the matched seismic wavelet obtained in the step 3)iiAnd τiThe value of (d) and the change Δ f of the instantaneous frequency of the envelope peak in the ith platelet obtained in step 5)p(i)Substituting into the analytic relational expression between the envelope peak instantaneous frequency and the quality factor Q value to obtain the estimated value of the quality factor Q value. The analytic relation between the envelope peak instantaneous frequency and the quality factor Q value is as follows:
<math> <mrow> <mi>Q</mi> <mo>&ap;</mo> <mfrac> <mrow> <mi>&tau;</mi> <msup> <mi>&delta;</mi> <mn>2</mn> </msup> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mn>4</mn> <mi>&pi;&Delta;</mi> <msub> <mi>f</mi> <mi>p</mi> </msub> </mrow> </mfrac> <mo>,</mo> </mrow> </math>
wherein, <math> <mrow> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> <mi>&eta;</mi> <msup> <mi>&Phi;</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mn>2</mn> <mi>&pi;&eta;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <msup> <mrow> <mo>-</mo> <mn>2</mn> <mi>&pi;</mi> </mrow> <mn>2</mn> </msup> <msup> <mi>&eta;</mi> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </math> <math> <mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mo>&infin;</mo> </mrow> <mi>x</mi> </msubsup> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>t</mi> <mn>2</mn> </msup> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mi>dt</mi> </mrow> </math> is a probability integral function of a standard normal distribution, t is an integral variable,
Figure BDA0000452740090000034
7) repeating the steps 2) to 6), and sequentially calculating the quality factors of the N-1 small thin plates except the ith small thin plate; 8) predicting the oil-gas-bearing performance of the reservoir by using the Q values of the quality factors of the N small sheets obtained in the steps 6) and 7), and further obtaining an estimated attenuation curve.
In the step 2), in the process (2) of estimating the seismic wavelet parameters by using the MSMG-based high-order cumulant matching method, the estimation method comprises the following steps: calculating the high-order cumulant of the seismic record: in the convolution model of the seismic record, a reflection coefficient sequence is assumed to be a random process with zero mean, mutual independence, same distribution and non-Gaussian, a seismic wavelet is a stable linear time-invariant system, a corresponding seismic record is a random process with zero mean and stability, the k-order stability of the process is further assumed, and according to a BBR formula, for the high-order cumulant of the seismic record, the high-order cumulant of the seismic record is as follows:
<math> <mrow> <msub> <mi>c</mi> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&gamma;</mi> <mi>k</mi> </msub> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>
where k denotes the order, y denotes the seismic record, ckyRepresenting the higher order cumulant of the seismic record, cky1,…,τk-1) Representing the higher order cumulant as a function of the amount of time shift, τ representing the time shift, τ1,…τk-1Representing different amounts of time shift, n representing a certain sampling point, ykIs the k-order cumulant of the reflection coefficient; u (n) is seismic wavelet, and u (n + tau) represents seismic wavelet after time shift; secondly, matching the seismic wavelets by calculating the ratio of the high-order cumulants: matching seismic wavelets by calculating the ratio of cumulants, as given by equation (1)
Figure BDA0000452740090000036
Wherein r iskyRepresenting the ratio of k-order cumulants for the seismic records; thirdly, estimating the parameters of the seismic wavelets by using an objective function: the physical meaning of the objective function yields an expression of the objective function defined as:
<math> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </munder> <mo>|</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msup> <mi>u</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> </mrow> </math>
by using an objective function
Figure BDA0000452740090000038
Obtaining minimum parameter estimation of seismic wavelets; fourthly, estimating the high-order cumulant of the seismic signal based on the MSMG model: the calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>k</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>k</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>=</mo> <msub> <mi>&omega;</mi> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω1,…,ωnA variable representing a fourier transform of the characteristic function,
Figure BDA0000452740090000042
viis the derivative order, k ═ v1+v2+…vn,ω=[ω1,…,ωn]T(ii) a Calculating the ratio of one-dimensional cumulant according to 2, 4 and 6 orders cumulant in the actual estimation of the seismic wavelet to match the seismic wavelet, and expressing the ratio of various cumulants:
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>3</mn> </mfrac> <mo>+</mo> <mfrac> <mn>2</mn> <mn>3</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>b</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>5</mn> </mfrac> <mo>+</mo> <mfrac> <mn>4</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>d</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>3</mn> <mn>5</mn> </mfrac> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>2</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>3</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>e</mi> <mo>)</mo> </mrow> </mrow> </math>
where ρ iss(m) is the correlation coefficient, ρs(0) R denotes the cumulative quantity ratio, m is the time delay, ynRepresenting seismic records, r(6)And (4) expressing the ratio of the 6-order cumulant, substituting the calculation result of the formula (b) into the target function in the step (c), and performing inversion to obtain the parameters of the seismic source wavelet.
In the step 4), the instantaneous frequency of the received signals at the upper and lower interfaces of the small thin plate is as follows:
<math> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </mfrac> <mfrac> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <mi>dH</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>]</mo> </mrow> <mi>dt</mi> </mfrac> <mo>-</mo> <mi>H</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>]</mo> <mfrac> <mrow> <mi>dx</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mi>dt</mi> </mfrac> </mrow> <mrow> <mi>e</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&epsiv;</mi> <mi>d</mi> </msub> <msub> <mi>e</mi> <mi>m</mi> </msub> </mrow> </mfrac> <mo>,</mo> </mrow> </math>
wherein,
Figure BDA0000452740090000049
for the analysis part of the seismic signal x (t) calculated by wavelet transform, Im (-) denotes the imaginary part, CgFor the permissive condition, Wf (t, s) is the reversible integral transform of the continuous wavelet transform of seismic signals x (t), the real number s is a scale factor and s > 0; epsilondIs damping coefficient and 0 < epsilond<<1;e(t)=x2(t)+(H[x(t)])2Is the square of the instantaneous amplitude of the direct wave of the top and bottom received signals of the ith plateletmMax (e (t)) is the maximum square of the instantaneous amplitude.
In the step 8), the process of predicting the oil-gas content of the reservoir by using the quality factor Q value is as follows: firstly, defining a target area range to be surveyed according to logging information; finding out the position corresponding to the target area defined in the step one on the estimated attenuation curve; checking whether the estimated attenuation curve is matched with the amplitude and the trend of the priori logging information or not, thereby checking whether the parameters of the seismic source wavelet and the parameters of the three-parameter wavelet determined by a high-order cumulant matching method based on the MSMG model are correct or not; if the target area on the attenuation profile is matched with the prior information, the parameter selection is proper, and the oil-gas content of other blocks is predicted by using the strong absorption area on the attenuation profile; and otherwise, continuing the steps, and re-estimating the parameters of the source wavelet and the matched seismic wavelet by using an iterative inversion method.
In the third step, the calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>r</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>r</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>=</mo> <msub> <mi>&omega;</mi> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> <mo>,</mo> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω ═ ω1,…,ωn]TThe MSMG model was introduced to simulate the multidimensional probability density of seismic recordings.
Due to the adoption of the technical scheme, the invention has the following advantages: 1. the invention provides an analytic relation between the instantaneous frequency of the envelope peak value and the Q value of the quality factor, the Q value of the medium quality factor can be conveniently estimated by using the analytic relation, the oil-gas content of a reservoir is predicted by using the Q value of the quality factor, and the estimation precision of the quality factor is effectively improved. 2. The method estimates the Q value of the medium quality factor by using the instantaneous frequency at the envelope peak of the seismic signal, avoids the problem of time-adding window, and solves the problem of low attenuation estimation precision caused by inaccurate spectrum estimation possibly introduced by time-adding window. 3. The wavelet domain envelope peak value instantaneous frequency method provided by the invention has strong anti-reflection wave interference capability at the interface, and the longitudinal resolution is 10-20 meters higher than that of the common method. 4. The invention has better random noise resistance by calculating the instantaneous frequency in the effective signal energy distribution space of the wavelet domain. 5. The method has high precision in estimating the parameters of the seismic source wavelet and the parameters of the kernel function by using a high-order cumulant matching method based on the MSMG model. 6. The wavelet domain envelope peak value instantaneous frequency method provided by the invention is used for models and practical examples, and the obtained medium has good corresponding relation between the absorption strength of the medium to seismic waves and the oil-gas content of a reservoir. The invention can be widely applied in the field of seismic survey.
Drawings
FIG. 1 is a schematic sheet of the present invention for zero-offset VSP (vertical seismic profile) data; wherein the thick solid line is the recording time, the thin solid line is the instantaneous frequency, and the dashed line is the instantaneous amplitude.
Detailed Description
The invention is described in detail below with reference to the figures and examples.
The invention provides a medium quality factor estimation method based on a seismic signal envelope peak, and particularly provides a seismic surveying method for estimating a medium quality factor by using VSP (vertical seismic profiling) data based on the change of instantaneous frequency at the seismic signal envelope peak. Which comprises the following steps:
1) the distance between two adjacent detectors in the well is taken as the thickness to divide the medium into a plurality of small thin plates.
As shown in FIG. 1, the top and the bottom of the rectangle are two adjacent detectors, and the medium is divided into N small thin plates according to the depth indication direction by taking the distance between the two adjacent detectors as the thickness, wherein N is a natural number.
2) For the ith thin plate, a constant phase wavelet is used for approximating a direct wave received at the top of the thin plate, and the seismic wavelet parameters are estimated by a high-order cumulant matching method based on MSMG (multivariate Gaussian mixture model).
The method for receiving the direct wave by approximating the top of the small thin plate by the constant phase wavelet comprises the following steps:
(1) in the horizontal layered viscoelastic medium, the quality factor Q value of each layer is set as a constant (i.e. Q is independent of frequency), only one-way wave propagation of a plane wave is considered, and a frequency expression of a seismic source wavelet located at the earth surface when the seismic source wavelet reaches a depth z is given out:
<math> <mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>)</mo> </mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mo>[</mo> <mo>-</mo> <mfrac> <mi>i&omega;&Delta;z</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mfrac> <mi>&omega;&Delta;z</mi> <mrow> <mn>2</mn> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>Q</mi> </mrow> </mfrac> <mo>]</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>
where Δ z is the transmission distance of the source wavelet, ω is the angular frequency of the source wavelet, G (Δ z) is a factor independent of frequency and absorption (including geometric dispersion),
Figure BDA0000452740090000062
in the frequency domain expression of the source wavelet, exp is an exponential function with the base of the natural logarithm e,
Figure BDA0000452740090000063
c (ω) is the phase velocity of the source wavelet.
(2) The following ordinary phase wavelets are adopted to approximate the seismic source wavelet, and the function expression of the ordinary phase wavelets is as follows:
Figure BDA0000452740090000064
wherein,the phase constant of the seismic source wavelet is shown, sigma is a scale parameter of the seismic source wavelet, the width of the wavelet is controlled, delta is a frequency parameter of the seismic source wavelet, and t is the propagation time of the seismic source wavelet. Since the wavelet in equation (2) has 3 undetermined parameters
Figure BDA0000452740090000066
Therefore, the method can better approximate the actual source wavelet than Ricker wavelet, band-pass wavelet or ideal pulse signal and the like.
The process of estimating the parameters of the seismic wavelets by using the MSMG-based high-order cumulant matching method comprises the following steps:
(1) estimation of phase parameters of seismic wavelets using existing phase rotation methods
Figure BDA0000452740090000067
(2) Estimating a scale parameter sigma and a frequency parameter delta by using a Cumulant Matching (CM) method, which mainly comprises the following steps:
calculating the high-order cumulant of the earthquake record
In the convolution model of seismic record, assuming that the reflection coefficient sequence is a random process with zero mean, mutual independence and same distribution and non-Gaussian, the seismic wavelet is a stable linear time-invariant system, the corresponding seismic record is a random process with zero mean and stability, further assuming that the k-th order of the process is stable, according to the BBR (Bartlett-Brillnger-Rosenblatt) formula, for the high-order cumulant of the seismic record, there are:
<math> <mrow> <msub> <mi>c</mi> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mrow> <mo>,</mo> <mi>&tau;</mi> </mrow> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&gamma;</mi> <mi>k</mi> </msub> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>
in the formula: k is the order, k is the stationary signal, y is the seismic record, ckyRepresenting the higher order cumulant of the seismic record, cky1,…,τk-1) Representing the higher order cumulant as a function of the amount of time shift, τ representing the time shift, τ1,…τk-1Representing different amounts of time shift, n representing a certain sampling point, ykIs the k-order cumulant of the reflection coefficient; u (n) is a seismic wavelet, and u (n + τ) represents a seismic wavelet after time shift.
② matching the seismic wavelets by calculating the ratio of higher order cumulants
Since the distribution of the reflection coefficient is unknown, γkIt is difficult to obtain, and the seismic wavelets can be matched by calculating the ratio of the cumulant, as can be obtained from equation (3)
Figure BDA0000452740090000071
Wherein r iskyThe ratio of k-order cumulants for seismic records is expressed.
Third, estimating the parameters of the seismic wavelet by the target function
The physical meaning of the objective function is to match the seismic wavelets by minimizing the difference between the estimated cumulants and those computed from the seismic records, and according to its physical meaning, the expression of the objective function is defined as:
<math> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </munder> <mo>|</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <munder> <mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> <mi>n</mi> </munder> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msup> <mi>u</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> </mrow> </math>
by using an objective function
Figure BDA0000452740090000079
The minimum is found to obtain a parameter estimate for the seismic wavelet.
Fourthly, estimating high-order cumulant of seismic signal based on MSMG model
The MSMG model can simulate the multidimensional probability density of the seismic record, so that the estimation of the accumulated quantity can be obtained more accurately. The calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>k</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>k</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>=</mo> <msub> <mi>&omega;</mi> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω1,…,ωnA variable representing a fourier transform of the characteristic function,
Figure BDA0000452740090000075
viis the derivative order, k ═ v1+v2+…vn,ω=[ω1,…,ωn]T
Fifthly, calculating the ratio of one-dimensional cumulant to match the seismic wavelets
When the seismic sub-wave is actually estimated, only one-dimensional cumulant ratio (ODCR) estimation needs to be calculated according to 2, 4 and 6 orders cumulant, and the cumulant ratio is calculated as the following form:
<math> <mrow> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>3</mn> </msubsup> <msub> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>4</mn> </msubsup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>2</mn> </msubsup> <msubsup> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mn>2</mn> </msubsup> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>4</mn> </msubsup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>5</mn> </msubsup> <msub> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>6</mn> </msubsup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>4</mn> </msubsup> <msubsup> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mn>2</mn> </msubsup> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>6</mn> </msubsup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msubsup> <mi>r</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>r</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>3</mn> </msubsup> <msubsup> <mi>u</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mn>3</mn> </msubsup> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msubsup> <mi>u</mi> <mi>n</mi> <mn>6</mn> </msubsup> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow> </math>
in the formula: m is a time delay; u represents the windowed intercepted seismic wavelet, n represents the sampling point, r represents the cumulant ratio, u3Denotes the power of u to the 3, ynRepresenting seismic records, r(6)The ratio of the 6 th order cumulant is expressed.
Based on the formulas (5) to (9), and ρs(0) 1, various ODCR expressions can be obtained:
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>3</mn> </mfrac> <mo>+</mo> <mfrac> <mn>2</mn> <mn>3</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mi>b</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>5</mn> </mfrac> <mo>+</mo> <mfrac> <mn>4</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mi>d</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>3</mn> <mn>5</mn> </mfrac> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>2</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>3</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mi>e</mi> <mo>)</mo> </mrow> </mrow> </math>
where ρ iss(m) is a correlation coefficient, which can be calculated from the actual seismic record. Experimental research shows that after multiple experiments, the parameters estimated by the formula (10 b) are most accurate, the error is minimum, the effect is best, the matching performance under the condition of low signal-to-noise ratio is best, and the estimation error fluctuation is minimum, so that the calculation result of the formula (10 b) is substituted into the target function in the step (c) to obtain the parameters of the seismic source wavelet through inversion.
3) And estimating parameters of a kernel function of the wavelet transformation by using a high-order cumulant matching method based on MSMG.
For the noise-containing seismic signals, a good denoising effect can be obtained by calculating the instantaneous frequency in the wavelet domain. The key to obtain higher time-frequency resolution is to select a kernel function which can be well matched with the seismic signal to be analyzed. The invention selects the kernel function of wavelet transform as three-parameter wavelet, and the expression is as follows:
g(t,Γ)=exp[-ξ(t-β)2]{p(Γ)[cos(mt)-k(Γ)]+iq(Γ)sin(mt)}
where the vector Γ ═ (m, ξ, β) is a set of parameters m, ξ, β, m is the modulation frequency, ξ is the energy attenuation factor, β is the energy delay factor.
The parameters of the wavelet kernel function are as follows:
<math> <mrow> <mi>k</mi> <mrow> <mo>(</mo> <mi>&Gamma;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>m</mi> <mn>2</mn> </msup> <mrow> <mn>4</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>[</mo> <mi>cos</mi> <mrow> <mo>(</mo> <mi>m&beta;</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>i</mi> <mi>sin</mi> <mrow> <mo>(</mo> <mi>m&beta;</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <mi>q</mi> <mrow> <mo>(</mo> <mi>&Gamma;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>&Gamma;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>]</mo> </mrow> </math>
<math> <mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>&Gamma;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> <mi>&pi;</mi> </mfrac> <mo>)</mo> </mrow> <mrow> <mn>1</mn> <mo>/</mo> <mn>4</mn> </mrow> </msup> <msup> <mrow> <mo>{</mo> <mn>4</mn> <msup> <mi>cos</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>m&beta;</mi> <mo>)</mo> </mrow> <mo>[</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>m</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mrow> <mn>3</mn> <mi>m</mi> </mrow> <mn>2</mn> </msup> <mrow> <mn>8</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>]</mo> <mo>-</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>m</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>+</mo> <mn>1</mn> <mo>}</mo> </mrow> <mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msup> </mrow> </math>
<math> <mrow> <mi>q</mi> <mrow> <mo>(</mo> <mi>&Gamma;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> <mi>&pi;</mi> </mfrac> <mo>)</mo> </mrow> <mrow> <mn>1</mn> <mo>/</mo> <mn>4</mn> </mrow> </msup> <msup> <mrow> <mo>{</mo> <msup> <mrow> <mn>4</mn> <mi>sin</mi> </mrow> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>m&beta;</mi> <mo>)</mo> </mrow> <mo>[</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>m</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mrow> <mn>3</mn> <mi>m</mi> </mrow> <mn>2</mn> </msup> <mrow> <mn>8</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>]</mo> <mo>-</mo> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>m</mi> <mn>2</mn> </msup> <mrow> <mn>2</mn> <mi>&xi;</mi> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>+</mo> <mn>1</mn> <mo>}</mo> </mrow> <mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> </mrow> </msup> </mrow> </math>
wherein, m, xi, beta are three parameters, if the parameters are selected reasonably, the three-parameter wavelet g (t, Γ) ═ exp [ -xi (t-beta) ]2]{p(Γ)[cos(mt)-k(Γ)]+ iq (Γ) sin (mt) } is the matching seismic wavelet. The estimation of the three parameters can be obtained by a high-order cumulative method based on the MSMG model, namely, the estimation is obtained by utilizing a formula for an objective function
Figure BDA0000452740090000093
Three parameters sigma of minimum obtained matching seismic wavelet are solvediiAnd τiThe value of (c).
Wherein the expression of the objective function is as follows:
<math> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>&xi;</mi> <mo>,</mo> <mi>&beta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </munder> <mo>|</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>m</mi> <mo>,</mo> <mi>&xi;</mi> <mo>,</mo> <mi>&beta;</mi> <mo>)</mo> </mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>;</mo> <mi>m</mi> <mo>,</mo> <mi>&xi;</mi> <mo>,</mo> <mi>&beta;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <mi>m</mi> <mo>,</mo> <mi>&xi;</mi> <mo>,</mo> <mi>&beta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msup> <mi>g</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>m</mi> <mo>,</mo> <mi>&xi;</mi> <mo>,</mo> <mi>&beta;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> <mo>,</mo> </mrow> </math>
in the formula:
4) calculating the instantaneous frequency of the received signal at the upper and lower interfaces of the ith small sheet
For seismic signals x (t) with high signal-to-noise ratio, the instantaneous frequency can be directly calculated by using the definition of the instantaneous frequency, namely the derivative of the instantaneous phase, so that the instantaneous frequency of the received signal at the upper and lower interfaces of the small sheet is as follows:
<math> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </mfrac> <mfrac> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <mi>dH</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mi>dt</mi> </mfrac> <mo>-</mo> <mi>H</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>]</mo> <mfrac> <mrow> <mi>dx</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mi>dt</mi> </mfrac> </mrow> <mrow> <mi>e</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&epsiv;</mi> <mi>d</mi> </msub> <msub> <mi>e</mi> <mi>m</mi> </msub> </mrow> </mfrac> <mo>,</mo> </mrow> </math>
wherein,
Figure BDA0000452740090000097
for the analysis part of the seismic signal x (t) calculated by wavelet transform, Im (-) denotes the imaginary part, CgFor the permissive condition, Wf (t, s) is the reversible integral transform of the continuous wavelet transform of seismic signals x (t), the real number s is a scale factor and s > 0; epsilondIs damping coefficient and 0 < epsilond<<1;e(t)=x2(t)+(H[x(t)])2Is the square of the instantaneous amplitude of the direct wave of the top and bottom received signals of the ith plateletmMax (e (t)) is the maximum square of the instantaneous amplitude. Respectively picking up the arrival time of the envelope peak of the top and bottom direct waves according to the position of the instantaneous frequency of the envelope peak as
Figure BDA0000452740090000098
Andthen calculating the travel time of the seismic waves in the ith small thin plate <math> <mrow> <msub> <mi>&tau;</mi> <mi>i</mi> </msub> <mo>=</mo> <msubsup> <mi>t</mi> <mi>i</mi> <mi>b</mi> </msubsup> <mo>-</mo> <msubsup> <mi>t</mi> <mi>i</mi> <mi>t</mi> </msubsup> <mo>.</mo> </mrow> </math>
5) And calculating the change of the instantaneous frequency of the envelope peak in the ith small thin plate.
Using formulasCalculating the instantaneous frequency of the received signal at the upper and lower interfaces of the ith small thin plate, and respectively recording as:and
Figure BDA0000452740090000104
respectively extracting the envelope peak instantaneous frequency of the direct wave at the upper and lower interfaces of the ith small thin plate, and respectively recording as:
Figure BDA0000452740090000105
then using the formula Δ fp(i)=fp(0)-fpi) Calculating to obtain the change of the instantaneous frequency of the envelope peak value in the ith small thin plate;
6) calculating the Q value of the quality factor of the ith small sheet
Three parameters sigma of the matched seismic wavelet obtained in the step 3)iiAnd τiThe value of (d) and the change Δ f of the instantaneous frequency of the envelope peak in the ith platelet obtained in step 5)p(i)Substituting into the analytic relational expression between the envelope peak instantaneous frequency and the quality factor Q value to obtain the estimated value of the quality factor Q value. The analytic relation between the envelope peak instantaneous frequency and the quality factor Q value is as follows:
<math> <mrow> <mi>Q</mi> <mo>&ap;</mo> <mfrac> <mrow> <mi>&tau;</mi> <msup> <mi>&delta;</mi> <mn>2</mn> </msup> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mn>4</mn> <mi>&pi;&Delta;</mi> <msub> <mi>f</mi> <mi>p</mi> </msub> </mrow> </mfrac> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow> </math>
wherein: <math> <mrow> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> <mi>&eta;</mi> <msup> <mi>&Phi;</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mn>2</mn> <mi>&pi;&eta;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mn>2</mn> <msup> <mi>&pi;</mi> <mn>2</mn> </msup> <msup> <mi>&eta;</mi> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </math> <math> <mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mo>&infin;</mo> </mrow> <mi>x</mi> </msubsup> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>t</mi> <mn>2</mn> </msup> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mi>dt</mi> </mrow> </math> is a probability integral function of a standard normal distribution, t is an integral variable,
7) and repeating the steps 2) to 6), and sequentially calculating the quality factors of the N-1 small thin plates except the ith small thin plate.
8) Predicting the oil-gas-bearing performance of the reservoir by using the Q values of the quality factors of the N small sheets obtained in the steps 6) and 7), and further obtaining an estimated attenuation curve. A large amount of data show that the quality factor Q value of soil and a surface layer is minimum, the quality factor Q value of sandstone is large, the quality factor Q value of shale is large, and the quality factor Q value of limestone is maximum. The sandstone has obviously enhanced absorptivity when containing oil gas, and the Q value of the quality factor is reduced. In summary, the better the elasticity of the medium, the larger the quality factor Q value. Therefore, the oil-gas-bearing property of the reservoir can be predicted by using the quality factor Q value, and the prediction implementation process is as follows:
firstly, defining a target area range to be surveyed according to logging information;
finding out the position corresponding to the target area defined in the step one on the estimated attenuation curve;
checking whether the estimated attenuation curve is matched with the amplitude and the trend of the priori logging information or not, and checking whether the parameters of the seismic source wavelet and the parameters of the three-parameter wavelet determined by a high-order cumulant matching method based on the MSMG model are correct or not;
if the target area on the attenuation profile is matched with the prior information, the parameter selection is proper, and the oil-gas content of other blocks can be predicted by using a strong absorption area on the attenuation profile; if the seismic source wavelet and the matched seismic wavelet are not matched, the steps are continued, and the parameters of the seismic source wavelet and the matched seismic wavelet are re-estimated by using an iterative inversion method.
In the above embodiment, the calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>r</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>r</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>&omega;</mi> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω ═ ω1,…,ωn]TThe MSMG model is introduced to simulate the multidimensional probability density of the seismic record, so that the estimation of the accumulated amount can be obtained more accurately. By modeling the joint probability density of the seismic records in the MSMG model, the high-order cumulant of the seismic records can be expressed as a polynomial of second-order statistic (correlation function), which brings great advantages to cumulant estimation, overcomes the problem that the work can be finished only by a large amount of data in the past, improves the estimation precision and has extremely high noise immunity (simulation experiments show that the lowest noise-to-noise ratio can reach 6 dB). The cumulant is used for estimating the parameters of the matched seismic wavelets, so that the satisfactory estimation effect is achieved in a shorter data set under a lower signal-to-noise ratio, and the problem of the matched wavelet structure in the high-precision frequency spectrum imaging of the seismic data is solved.
In each of the above embodiments, the damping coefficient εdThe introduction of (2) can eliminate the peak appearing at the instantaneous frequency at the smaller amplitude of the signal, and has little influence on the instantaneous frequency at the larger amplitude of the signal.
In the above embodiments, the analytic relational expression between the envelope peak instantaneous frequency and the quality factor Q value is obtained by analyzing the characteristics of seismic wavelets extracted from the zero-offset VSP data, and performing approximate expansion by using a first-order taylor series in combination with the quality factor Q values of different media and the range of travel time in the thin plate. It shows that the attenuation of the seismic wave is related to the variation of the envelope peak instantaneous frequency, the propagation time and the wavelet parameters. The invention estimates the attenuation parameter by the change of the instantaneous frequency of the envelope peak in a certain transmission time. If the change of the instantaneous frequency of the envelope peak value in a specific transmission time is large, the attenuation or the absorption is strong; otherwise, the attenuation or absorption is weak.
The above embodiments are only used for illustrating the present invention, and the structure, connection mode, manufacturing process, etc. of the components may be changed, and all equivalent changes and modifications performed on the basis of the technical solution of the present invention should not be excluded from the protection scope of the present invention.

Claims (5)

1. A stratum medium quality factor estimation method based on seismic signal envelope peaks comprises the following steps:
1) dividing the medium into a plurality of small thin plates by taking the distance between two adjacent detectors in the well as the thickness;
2) for the ith small thin plate, a common phase wavelet is used for approaching a direct wave received at the top of the small thin plate, and the seismic wavelet parameters are estimated by a high-order cumulant matching method based on MSMG; wherein MSMG is a multivariate Gaussian mixture model; the method for approaching the top of the small thin plate to receive the direct wave by using the constant phase wavelet comprises the following steps:
(1) in the horizontal layered viscoelastic medium, the quality factor Q value of each layer is taken as a constant, and only the single-pass wave propagation of the plane wave is considered, so that the frequency expression of the seismic source wavelet at the earth surface when the seismic source wavelet reaches the depth z is as follows:
<math> <mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>&Delta;z</mi> <mo>)</mo> </mrow> <mover> <mi>U</mi> <mo>^</mo> </mover> <mrow> <mo>(</mo> <mn>0</mn> <mo>,</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mo>[</mo> <mo>-</mo> <mfrac> <mi>i&omega;&Delta;z</mi> <mrow> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mfrac> <mi>&omega;&Delta;z</mi> <mrow> <mn>2</mn> <mi>c</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> <mi>Q</mi> </mrow> </mfrac> <mo>]</mo> <mo>,</mo> </mrow> </math>
where Δ z is the transmission distance of the source wavelet, ω is the angular frequency of the source wavelet, G (Δ z) is a factor independent of frequency and absorption,
Figure FDA0000452740080000012
in the frequency domain expression of the source wavelet, exp is an exponential function with the base of the natural logarithm e,
Figure FDA0000452740080000013
c (omega) is the phase velocity of the seismic source wavelet;
(2) adopting a constant phase wavelet to approximate a seismic source wavelet, wherein the function expression of the constant phase wavelet is as follows:
Figure FDA0000452740080000014
wherein,
Figure FDA0000452740080000015
the phase constant of the seismic source wavelet is used, sigma is a scale parameter of the seismic source wavelet, the width of the wavelet is controlled, delta is a frequency parameter of the seismic source wavelet, and t is the propagation time of the seismic source wavelet;
the process of estimating the parameters of the seismic wavelets by using the MSMG-based high-order cumulant matching method comprises the following steps:
(1) estimation of phase parameters of seismic wavelets using existing phase rotation methods
Figure FDA00004527400800000110
(2) Estimating a scale parameter sigma and a frequency parameter delta by using an cumulant matching method;
3) estimating parameters of a kernel function of wavelet transformation by a high-order cumulant matching method based on MSMG;
4) calculating the instantaneous frequency of the received signal at the upper and lower interfaces of the ith small thin plate;
5) calculating the change of the instantaneous frequency of the envelope peak value in the ith small thin plate;
6) calculating the Q value of the quality factor of the ith small sheet: three parameters sigma of the matched seismic wavelet obtained in the step 3)iiAnd τiThe value of (d) and the change Δ f of the instantaneous frequency of the envelope peak in the ith platelet obtained in step 5)p(i)Substituting into the analytic relational expression between the envelope peak instantaneous frequency and the quality factor Q value to obtain the estimated value of the quality factor Q value. The analytic relation between the envelope peak instantaneous frequency and the quality factor Q value is as follows:
<math> <mrow> <mi>Q</mi> <mo>&ap;</mo> <mfrac> <mrow> <mi>&tau;</mi> <msup> <mi>&delta;</mi> <mn>2</mn> </msup> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mn>4</mn> <mi>&pi;&Delta;</mi> <msub> <mi>f</mi> <mi>p</mi> </msub> </mrow> </mfrac> <mo>,</mo> </mrow> </math>
wherein, <math> <mrow> <mi>&kappa;</mi> <mrow> <mo>(</mo> <mi>&eta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> <mi>&eta;</mi> <msup> <mi>&Phi;</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mn>2</mn> <mi>&pi;&eta;</mi> <mo>)</mo> </mrow> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mn>2</mn> <msup> <mi>&pi;</mi> <mn>2</mn> </msup> <msup> <mi>&eta;</mi> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </math> <math> <mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msqrt> <mn>2</mn> <mi>&pi;</mi> </msqrt> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mo>&infin;</mo> </mrow> <mi>x</mi> </msubsup> <mi>exp</mi> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>t</mi> <mn>2</mn> </msup> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> <mi>dt</mi> </mrow> </math> is a probability integral function of a standard normal distribution,t is the variable of the integral and is,
Figure FDA0000452740080000021
7) repeating the steps 2) to 6), and sequentially calculating the quality factors of the N-1 small thin plates except the ith small thin plate;
8) predicting the oil-gas-bearing performance of the reservoir by using the Q values of the quality factors of the N small sheets obtained in the steps 6) and 7), and further obtaining an estimated attenuation curve.
2. The method of claim 1, wherein the method comprises the steps of: in the step 2), in the process (2) of estimating the seismic wavelet parameters by using the MSMG-based high-order cumulant matching method, the estimation method comprises the following steps:
calculating the high-order cumulant of the seismic record:
in the convolution model of the seismic record, a reflection coefficient sequence is assumed to be a random process with zero mean, mutual independence, same distribution and non-Gaussian, a seismic wavelet is a stable linear time-invariant system, a corresponding seismic record is a random process with zero mean and stability, the k-order stability of the process is further assumed, and according to a BBR formula, for the high-order cumulant of the seismic record, the high-order cumulant of the seismic record is as follows:
<math> <mrow> <msub> <mi>c</mi> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&gamma;</mi> <mi>k</mi> </msub> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </math>
where k denotes the order, y denotes the seismic record, ckyRepresenting the higher order cumulant of the seismic record, cky1,…,τk-1) Representing the higher order cumulant as a function of the amount of time shift, τ representing the time shift, τ1,…τk-1Representing different amounts of time shift, n representing a certain sampling point, ykIs the k-order cumulant of the reflection coefficient; u (n) is seismic wavelet, and u (n + tau) represents seismic wavelet after time shift;
secondly, matching the seismic wavelets by calculating the ratio of the high-order cumulants:
matching seismic wavelets by calculating the ratio of cumulants, as given by equation (1)
Figure FDA0000452740080000023
Wherein r iskyRepresenting the ratio of k-order cumulants for the seismic records;
thirdly, estimating the parameters of the seismic wavelets by using an objective function:
the physical meaning of the objective function yields an expression of the objective function defined as:
<math> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>&Sigma;</mi> <mrow> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> </mrow> </munder> <mo>|</mo> <msub> <mover> <mi>r</mi> <mo>^</mo> </mover> <mi>ky</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>,</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mn>1</mn> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <msub> <mi>&tau;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mi>&Sigma;</mi> <mi>n</mi> </munder> <msup> <mi>u</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>n</mi> <mo>;</mo> <mi>&sigma;</mi> <mo>,</mo> <mi>&delta;</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>|</mo> </mrow> </math>
by using an objective function
Figure FDA0000452740080000025
Obtaining minimum parameter estimation of seismic wavelets;
fourthly, estimating the high-order cumulant of the seismic signal based on the MSMG model:
the calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo></mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>k</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>k</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mrow> <mo>=</mo> <mi>&omega;</mi> </mrow> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω1,…,ωnA variable representing a fourier transform of the characteristic function,
Figure FDA0000452740080000032
viis the derivative order, k ═ v1+v2+…vn,ω=[ω1,…,ωn]T
Calculating the ratio of one-dimensional cumulant according to 2, 4 and 6 orders cumulant in the actual estimation of the seismic wavelet to match the seismic wavelet, and expressing the ratio of various cumulants:
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>3</mn> </mfrac> <mo>+</mo> <mfrac> <mn>2</mn> <mn>3</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>b</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>c</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>5</mn> </mfrac> <mo>+</mo> <mfrac> <mn>4</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>d</mi> <mo>)</mo> </mrow> </mrow> </math>
<math> <mrow> <msubsup> <mi>r</mi> <mi>y</mi> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mi>y</mi> <mi>n</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>y</mi> <mrow> <mi>n</mi> <mo>+</mo> <mi>m</mi> </mrow> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>3</mn> <mn>5</mn> </mfrac> <msub> <mi>&rho;</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>2</mn> <mn>5</mn> </mfrac> <msubsup> <mi>&rho;</mi> <mi>s</mi> <mn>3</mn> </msubsup> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>e</mi> <mo>)</mo> </mrow> </mrow> </math>
where ρ iss(m) is the correlation coefficient, ρs(0) R denotes the cumulative quantity ratio, m is the time delay, ynRepresenting seismic records, r(6)And (4) expressing the ratio of the 6-order cumulant, substituting the calculation result of the formula (b) into the target function in the step (c), and performing inversion to obtain the parameters of the seismic source wavelet.
3. The method for estimating the formation medium quality factor based on the envelope peak of the seismic signal as claimed in claim 1 or 2, wherein: in the step 4), the instantaneous frequency of the received signals at the upper and lower interfaces of the small thin plate is as follows:
<math> <mrow> <mi>f</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> </mfrac> <mfrac> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mfrac> <mrow> <mi>dH</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>]</mo> </mrow> <mi>dt</mi> </mfrac> <mo>-</mo> <mi>H</mi> <mo>[</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>]</mo> <mfrac> <mrow> <mi>dx</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mi>dt</mi> </mfrac> </mrow> <mrow> <mi>e</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&epsiv;</mi> <mi>d</mi> </msub> <msub> <mi>e</mi> <mi>m</mi> </msub> </mrow> </mfrac> <mo>,</mo> </mrow> </math>
wherein,
Figure FDA0000452740080000039
for the analysis part of the seismic signal x (t) calculated by wavelet transform, Im (-) denotes the imaginary part, CgFor the permissive condition, Wf (t, s) is the reversible integral transform of the continuous wavelet transform of seismic signals x (t), the real number s is a scale factor and s > 0; epsilondIs damping coefficient and 0 < epsilond<<1;e(t)=x2(t)+(H[x(t)])2Is the ith small filmSquare of instantaneous amplitude of direct wave of top and bottom received signals of board, emMax (e (t)) is the maximum square of the instantaneous amplitude.
4. The method for estimating the formation medium quality factor based on the envelope peak of the seismic signal as claimed in claim 1 or 2, wherein: in the step 8), the process of predicting the oil-gas content of the reservoir by using the quality factor Q value is as follows:
firstly, defining a target area range to be surveyed according to logging information;
finding out the position corresponding to the target area defined in the step one on the estimated attenuation curve;
checking whether the estimated attenuation curve is matched with the amplitude and the trend of the priori logging information or not, thereby checking whether the parameters of the seismic source wavelet and the parameters of the three-parameter wavelet determined by a high-order cumulant matching method based on the MSMG model are correct or not;
if the target area on the attenuation profile is matched with the prior information, the parameter selection is proper, and the oil-gas content of other blocks is predicted by using the strong absorption area on the attenuation profile; and otherwise, continuing the steps, and re-estimating the parameters of the source wavelet and the matched seismic wavelet by using an iterative inversion method.
5. The method of claim 4, wherein the method comprises the steps of: in the third step, the calculation formula of the seismic signal high-order cumulant based on the MSMG model is as follows:
<math> <mrow> <msub> <mi>c</mi> <mrow> <msub> <mi>v</mi> <mn>1</mn> </msub> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mi>i</mi> <mi>r</mi> </msup> </mfrac> <mfrac> <mrow> <msup> <mo>&PartialD;</mo> <mi>r</mi> </msup> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>&omega;</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mn>1</mn> <msub> <mi>v</mi> <mn>1</mn> </msub> </msubsup> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&PartialD;</mo> <msubsup> <mi>&omega;</mi> <mi>n</mi> <msub> <mi>v</mi> <mi>n</mi> </msub> </msubsup> </mrow> </mfrac> <msub> <mo>|</mo> <mrow> <msub> <mi>&omega;</mi> <mn>1</mn> </msub> <mo>=</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <mo>&CenterDot;</mo> <msub> <mi>&omega;</mi> <mi>n</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </msub> <mo>,</mo> </mrow> </math>
where φ is a characteristic function of the seismic signal, ω ═ ω1,…,ωn]TThe MSMG model was introduced to simulate the multidimensional probability density of seismic recordings.
CN201410003193.6A 2014-01-03 2014-01-03 Method for estimating stratum medium quality factors based on seismic signal envelope peak Pending CN103728662A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410003193.6A CN103728662A (en) 2014-01-03 2014-01-03 Method for estimating stratum medium quality factors based on seismic signal envelope peak

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410003193.6A CN103728662A (en) 2014-01-03 2014-01-03 Method for estimating stratum medium quality factors based on seismic signal envelope peak

Publications (1)

Publication Number Publication Date
CN103728662A true CN103728662A (en) 2014-04-16

Family

ID=50452813

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410003193.6A Pending CN103728662A (en) 2014-01-03 2014-01-03 Method for estimating stratum medium quality factors based on seismic signal envelope peak

Country Status (1)

Country Link
CN (1) CN103728662A (en)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103984013A (en) * 2014-04-24 2014-08-13 浪潮电子信息产业股份有限公司 Wavelet domain pre-stack seismic trace set absorption attenuation parameter estimation algorithm
CN105388523A (en) * 2014-09-04 2016-03-09 中国石油化工股份有限公司 High-precision quality factor extraction method
CN106324663A (en) * 2015-06-17 2017-01-11 中国石油化工股份有限公司 Method for obtaining quality factor
CN106547019A (en) * 2015-09-17 2017-03-29 中国石油化工股份有限公司 A kind of method of definitely interval quality factors
CN106650609A (en) * 2016-10-26 2017-05-10 太原理工大学 J-wave detection and classification method based on tunable Q-factor wavelet transform and higher-order cumulant
CN108614295A (en) * 2018-05-15 2018-10-02 中国海洋石油集团有限公司 A kind of stratum Q value calculating methods based on broad sense seismic wavelet
CN109459788A (en) * 2017-09-06 2019-03-12 中国石油化工股份有限公司 Ground interval quality factors calculation method and system
CN110007340A (en) * 2019-02-01 2019-07-12 西安理工大学 Salt dome speed density estimation method based on the direct envelope inverting of angle domain
CN113009569A (en) * 2019-12-20 2021-06-22 中国石油天然气集团有限公司 Seismic migration imaging method and device
CN113805232A (en) * 2020-06-17 2021-12-17 中国石油化工股份有限公司 Method, system and storage medium for estimating quality factor of shallow earth surface

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102288992A (en) * 2011-04-26 2011-12-21 中国海洋石油总公司 Method for estimating quality factor of medium by using peak envelope instantaneous frequency of seismic signal

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102288992A (en) * 2011-04-26 2011-12-21 中国海洋石油总公司 Method for estimating quality factor of medium by using peak envelope instantaneous frequency of seismic signal

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
JINGHUAI GAO ET AL.: "《Estimation of quality factor Q from the instantaneous frequency at the envelope peak of a seismic signal》", 《JOURNAL OF COMPUTATIONAL ACOUSTICS》, vol. 19, no. 2, 31 December 2011 (2011-12-31) *
JING-HUAI GAO ET AL.: "《Estimation of Seismic Wavelets Based on the Multivariate Scale Mixture of Gaussians Model》", 《ENTROPY》, 28 December 2009 (2009-12-28), pages 14 - 26 *
高静怀等: "《利用 VSP 资料直达波的包络峰值处瞬时频率提取介质品质因子》", 《地球物理学报》, vol. 51, no. 3, 31 May 2008 (2008-05-31), pages 853 - 856 *
高静怀等: "《利用地震信号的包络峰值瞬时频率估计介质品质因子》", 《中国应用地球物理论文选集》, 30 September 2009 (2009-09-30), pages 374 - 379 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103984013A (en) * 2014-04-24 2014-08-13 浪潮电子信息产业股份有限公司 Wavelet domain pre-stack seismic trace set absorption attenuation parameter estimation algorithm
CN105388523A (en) * 2014-09-04 2016-03-09 中国石油化工股份有限公司 High-precision quality factor extraction method
CN106324663A (en) * 2015-06-17 2017-01-11 中国石油化工股份有限公司 Method for obtaining quality factor
CN106324663B (en) * 2015-06-17 2018-10-02 中国石油化工股份有限公司 A kind of acquisition methods of quality factor
CN106547019A (en) * 2015-09-17 2017-03-29 中国石油化工股份有限公司 A kind of method of definitely interval quality factors
CN106650609A (en) * 2016-10-26 2017-05-10 太原理工大学 J-wave detection and classification method based on tunable Q-factor wavelet transform and higher-order cumulant
CN106650609B (en) * 2016-10-26 2018-07-31 太原理工大学 Based on the detection of J waves and sorting technique for adjusting Q wavelet transformations and Higher Order Cumulants
CN109459788A (en) * 2017-09-06 2019-03-12 中国石油化工股份有限公司 Ground interval quality factors calculation method and system
CN108614295A (en) * 2018-05-15 2018-10-02 中国海洋石油集团有限公司 A kind of stratum Q value calculating methods based on broad sense seismic wavelet
CN110007340A (en) * 2019-02-01 2019-07-12 西安理工大学 Salt dome speed density estimation method based on the direct envelope inverting of angle domain
CN113009569A (en) * 2019-12-20 2021-06-22 中国石油天然气集团有限公司 Seismic migration imaging method and device
CN113805232A (en) * 2020-06-17 2021-12-17 中国石油化工股份有限公司 Method, system and storage medium for estimating quality factor of shallow earth surface
CN113805232B (en) * 2020-06-17 2024-04-09 中国石油化工股份有限公司 Estimation method, system and storage medium for quality factors of shallow earth surface

Similar Documents

Publication Publication Date Title
CN103728662A (en) Method for estimating stratum medium quality factors based on seismic signal envelope peak
RU2579164C1 (en) Handling method for determining quality of geologic environment
US6931324B2 (en) Method for determining formation quality factor from seismic data
CN105093294B (en) Attenuation of seismic wave gradient method of estimation based on variable mode decomposition
CN104502997B (en) A kind of method of utilization fracture spacing curve prediction fracture spacing body
CN104237945B (en) A kind of seismic data self adaptation high resolution processing method
CN109669212B (en) Seismic data processing method, stratum quality factor estimation method and device
CN102288992A (en) Method for estimating quality factor of medium by using peak envelope instantaneous frequency of seismic signal
CN102478670B (en) Method for forecasting reservoir fluid property through earthquake attenuation attribute
Yanhu et al. A method of seismic meme inversion and its application
CN103257363B (en) A kind of method of fracture dip in Underground fracture-type reservoir
CN102305941A (en) Method for determining stratum stack quality factor by direct scanning of prestack time migration
Gao et al. Estimation of quality factor Q from the instantaneous frequency at the envelope peak of a seismic signal
CN108020863A (en) A kind of thin and interbedded reservoir porosity prediction method based on earthquake parity function
CN103645499B (en) Based on the earth surface consistency vibration amplitude compensation method of poststack reflected energy statistics
CN105044777A (en) Method for detecting strong reflection amplitude elimination of seismic marker layer based on empirical mode decomposition
CN102169188A (en) Method for surveying oil and gas based on Morlet spectrum
CN112946752B (en) Method for predicting fracture probability body based on fracture density model
CN115327632A (en) Method for analyzing short-time-window spectral attenuation of dereflection coefficient
CN112711070B (en) Oil gas detection method and device based on seismic signal decomposition
CN103257362B (en) Carbonatite efficient well forecasting method based on pressure noise density difference inversion
CN114152985B (en) Method for determining boundary of underground ancient river channel and thickness of thin sand body in boundary
CN111505707B (en) Method for extracting dispersion curve from vertical seismic profile data
CN103984013B (en) A kind of wavelet field prestack seismic gather attenuation by absorption parameter estimation algorithm
CN110568491B (en) Quality factor Q estimation method

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20140416

RJ01 Rejection of invention patent application after publication