CN104880730A - Seismic data time-frequency analysis and attenuation estimation method based on Synchrosqueezing transform - Google Patents

Abstract

The invention discloses a seismic data time-frequency analysis and attenuation estimation method based on Synchrosqueezing transform. A new time-frequency analysis tool Synchrosqueezing transform is first used for seismic data time-frequency analysis, more aggregated time-frequency representation is obtained through the rearrangement of transform domain coefficients, the time-frequency resolution is greatly improved, the obtained time-frequency representation is used for actual seismic data analysis and tight sandstone model gas bearing detection, the location of a reservoir can be accurately defined, the geological structure like a river channel, a fault and the like can be indicated, and further material interpretation and well location determination are facilitated. The invention provides the seismic attenuation estimation method based on the Synchrosqueezing transform and gives a specific implementation process, the attenuation estimation results of three-dimensional data of an oilfield tight sandstone reservoir and drilling results have a good consistency, and the method can be used for helping geologists to indicate a gas bearing reservoir and determine a drilling location.

Description

Based on seismic data time frequency analysis and the decay behavior method of Synchrosqueezing conversion

Technical field

The invention belongs to the signal transacting field in geophysical survey, the time frequency analysis and the seismic attenuation that relate to seismic data are estimated, are specifically related to a kind of seismic data time frequency analysis based on Synchrosqueezing conversion and decay behavior method.

Background technology

Seismic signal is complicated non-stationary signal, and time frequency analysis can be used for describing seismic signal rule over time, portrays the local feature of seismic signal, thus reflects the geologic structure corresponding to these features and reservoir information.Therefore, time frequency analysis is the important means of seism processing and explanation.Along with the development of signal processing theory in modern age, emerge a large amount of time frequency analyzing tool and Time-Frequency Analysis Method, be wherein much also used in seismic signal analysis.

Seismic event is after gas-bearing reservoir, and rapidly, cause seismic event to reduce in the local dominant frequency in this region, this exception can be used to refer to the distribution of hydrocarbon reservoir in radio-frequency component decay.But this is abnormal not obvious on original seismic data, but can be enhanced on the section of frequency resolution.Therefore, different time frequency analyzing tool is used to detect this exception.In this course, the T/F resolution of time frequency analyzing tool just becomes the key of this problem.

Fourier conversion as classics analysis of spectrum instrument, seismic data processing a lot of in be applied, such as: spectrum analysis, noise compacting etc., it is also the indispensable module in current geophysical software.After whole track data is carried out Fourier conversion, the spectrum distribution of data can be obtained, but the localization characteristics of frequency cannot be reflected, namely cannot portray frequency over time, and a key character of this seismic signal exactly.Therefore, a large amount of time-frequency domain Conjoint Analysis instruments is used in seismic signal analysis, and the seismic signal in time domain is extended to time-frequency domain after conversion, and the time-frequency figure obtained reflects the situation of frequency content along with time variations.

Enlarge section (WFT) the window function intercept signal of movement, then does Fourier conversion to the signal inside each time window, obtains the T/F distribution of signal.Enlarge section is used for the spectral factorization of Gulfian block 3-d seismic data set by the people such as Partyka, determines position and the thickness of thin layer, portrays the uncontinuity of geologic body.Marfurt and Kirlin enlarge section analyzes thin layer tuning effect.

In enlarge section, the selection of window function is the compromise of temporal resolution and frequency resolution, once window function is selected, temporal resolution and the frequency resolution of time-frequency figure are also determined thereupon.As typical non-stationary signal, the frequency of seismic signal constantly changed along with the time, need to carry out multiscale analysis, namely the ramped composition of wider window analytic signal is used, by the mutagenic components of narrower window analytic signal, the window width of enlarge section is fixed, and obviously cannot be competent at the requirement of multiscale analysis.In order to address this problem, wavelet transformation arises at the historic moment.The earliest using wavelet transformation as time m-yardstick localization instrument proposes is geophysicist Morlet, after this, physicist Grossmann cooperates with Morlet, gives the strict difinition of wavelet transformation.Wavelet transformation controls the width of window function by yardstick, and obtain wide window function by large scale, small scale obtains narrow window function, has different T/F resolution, achieve the multiresolution analysis to signal in low-and high-frequency.Because window width is controlled by yardstick, thus signal through wavelet transformation obtain be time m-size distribution, be called wavelet scale spectrum (scalogram), the corresponding relation that yardstick and frequency are not determined, depends on the selection of morther wavelet.Wavelet transformation is used for seismic signal analysis by Chakraborty and Okaya, and compares with enlarge section, illustrates its superiority.The people such as Siha propose the concept of time-frequency continuous wavelet transform, the time m-scale spectrum that wavelet transformation produces are converted into T/F spectrum, for detecting the visuality of the low frequency shadow relevant with hydrocarbon reservoir and the tiny geologic structure of enhancing.The people such as Gao Jinghuai have studied the On The Choice of seism processing Wavelets, propose, with mating the function of seismic wavelet as wavelet, seismic data to be carried out to the method for denoising and frequency division explanation.The people such as Gao Jinghuai propose Three parameter wavelet, and analyze for thin interbed seismic data.Wavelet transformation is used for the 3-D data volume analysis of certain block of Germany by the people such as Kazemeini, and result can reflect the deposition in river course and the migration of oil gas.

Above-mentioned enlarge section, wavelet transformation are linear transformation, and its T/F associating resolution is subject to the restriction of uncertainty principle.Wigner-Ville distribution has the highest T/F associating resolution, is also widely used in seismic signal analysis.Wigner-Ville distribution is used for the spectral factorization of seismic data by Li and Zheng, portrays Carbonate Reservoir.

In recent years, existing time frequency analyzing tool, for the feature of seismic signal analysis, carried out promoting and developing by a lot of scholar, or new phase space analysis instrument is used for seismic signal.2013, a two-dimentional deconvolution operator is acted on enlarge section time-frequency spectrum by Lu and Liu, obtain deconvolution enlarge section time-frequency spectrum, improve time frequency resolution, the spectral factorization based on the method is used to the Carbonate Reservoir of the cave type portraying In A Certain Area of Tarim Basin block.2013, Han and van der Baan utilized empirical mode decomposition (EMD) to carry out the time frequency analysis of seismic signal, and for the Cretaceous period meandering river detection in certain basin of Canada and portraying of tiny geologic structure.

But, traditional Time-Frequency Analysis Method, as Short Time Fourier Transform and wavelet transformation etc. have limitation for during seismic signal time-frequency analysis, be mainly reflected in the impact due to time-frequency atom, cause the energy dispersal of signal on time-frequency plane and the decline of time frequency resolution.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, provides the high seismic data time frequency analysis based on Synchrosqueezing conversion of a kind of time frequency resolution and decay behavior method.

For achieving the above object, the present invention by the following technical solutions:

Seismic data Time-Frequency Analysis Method based on Synchrosqueezing conversion:

Synchrosqueezing conversion specifically comprises the following steps:

Step a: continuous wavelet transform

Signal wavelet transformation time domain and frequency domain be expressed as:

$W_f(a,b)\begin{array}{c}=\frac{1}{\sqrt{a}}{\∫}_{-\∞}^{+\∞}f\left(t\right){\mathrm{\ψ}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}\\ =\frac{\sqrt{a}}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}F\left(\mathrm{\ω}\right){\mathrm{\Ψ}}^{*}\left(\mathrm{a\ω}\right)e^{\mathrm{ib\ω}}\mathrm{d\ω},\end{array}---\left(7\right)$

Wherein ψ (t) is wavelet, a and b is respectively scale factor and shift factor, and F (ω) and Ψ (ω) converts for f (t) and the Fourier of ψ (t); Suppose that wavelet function does not almost have negative frequency components, namely as ω < 0, Ψ (ω) ≈ 0;

The calculating of step b:FM frequency, demodulation frequency

${\mathrm{\&omega;}}_f(a,b)=\left\{\begin{array}{cc}\frac{-i{\&PartialD;}_b{W}_f(a,b)}{W_f(a,b)}& \left|W_f(a,b)\right|>0,\\ \&infin;& \left|W_f(a,b)\right|=0.\end{array}\right.---\left(8\right);$

Step c: time m-scale domain to the mapping of time-frequency domain

(1) conitnuous forms

$T_f(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_f(a,b)=\mathrm{\&omega;},W_f(a,b)\&NotEqual;0\}}{W}_f(a,b)a^{-3/2}\mathrm{da}---\left(9\right)$

By formula (7), time wavelet coefficient corresponding to m-scale domain all and frequencies omega combine, at time-frequency domain again by the position of energy " extruding " to frequencies omega place;

(2) discrete form

When carrying out numerical evaluation, need to carry out discrete to the yardstick a in formula (7) and frequencies omega.Yardstick after discretize is designated as { a _k, wherein a _k> a _k-1, yardstick is spaced apart a _k-a _k-1=(Δ is a) _k; Frequency is divided, is designated as { ω _l, wherein ω _l> ω _l-1, frequency interval is ω _l-ω _l-1=Δ ω; The discrete form of Synchrosqueezing conversion is expressed as:

$T_f({\mathrm{\&omega;}}_l,b)=\underset{a_k:|{\mathrm{\&omega;}}_f(a_k,b)-{\mathrm{\&omega;}}_l|\&le;\mathrm{\&Delta;\&omega;}/2}{\mathrm{\&Sigma;}}{W}_f(a_k,b){a_k}^{-3/2}{\left(\mathrm{\&Delta;a}\right)}_k,---\left(10\right);$

Seismic data Time-Frequency Analysis Method specifically comprises the following steps:

Step 1: choose typical road in 3-D data volume seismic section, carries out Synchrosqueezing conversion to this track data, finds out abnormal area respective frequencies;

Step 2: carry out Synchrosqueezing conversion to whole seismic section, extracts the section of abnormal area respective frequencies;

Step 3: carry out Synchrosqueezing conversion to whole 3-D data volume, obtains frequency data body, then extracts a horizon slice, carries out seismic data interpretation for geological personnel.

Carry out the method for seismic attenuation estimation based on Synchrosqueezing conversion, comprise the following steps:

Step 1: determine destination layer scope, carries out spectrum analysis to the 3-d seismic data set near destination layer, chooses suitable high frequency f _hwith low frequency f _l;

Step 2: utilize Synchrosqueezing to convert the single-frequency data volume T (x, y, t, the f that obtain high frequency _h) and single-frequency data volume T (x, y, t, the f of low frequency _l);

Step 3: the layer position H above destination layer _anear (x, y), the amplitude difference of high and low frequency is eliminated in advance;

Step 4: the decay near estimating target layer.

Further, in described step 1, high frequency f _hamplitude and low frequency f _lamplitude roughly the same.

Further, in described step 3, utilize formula (14) to calculate modifying factor α (x, y) and do level and smooth;

$\mathrm{\&alpha;}(x,y)=\frac{T(x,y,H_A(x,y),f_L)}{T(x,y,H_A(x,y),f_H)}.---\left(11\right).$

Further, in described step 4, by the decay AS (x, y, t) near formula (15) estimating target layer;

AS(x,y,t)＝T(x,y,t,f _L)-α(x,y)T(x,y,t,f _H)(12)。

New time frequency analyzing tool Synchrosqueezing conversion is used for seismic data time frequency analysis by the present invention first, and this conversion is by the rearrangement to coefficient in transform domain, and obtain a time-frequency representation more assembled, time frequency resolution improves greatly.Use it for actual seismic analysis and tight sand model Gas potential detection, accurately can define the position of reservoir, the instruction geologic structure such as river course and tomography, and then be conducive to further data interpretation and well location is determined.

The present invention also proposes the seismic attenuation method of estimation based on Synchrosqueezing conversion, and provide specific implementation flow process, good consistance is had to the decay behavior result of certain oil field Sandstone Gas Reservoir 3-D data volume and drilling well result, the method can help geological personnel to indicate gas-bearing reservoir, determines drilling well position.

Accompanying drawing explanation

Fig. 1 is wavelet transformation and the Synchrosqueezing Transformation Graphs of cosine signal;

(a) cosine signal; (b) wavelet transformation; C () Synchrosqueezing converts;

Fig. 2 is that Synchrosqueezing converts schematic diagram;

Fig. 3 is test signal and the result adopting different Time-Frequency Analysis Method;

(a) test signal; (b) enlarge section; (c) wavelet transformation; D () Synchrosqueezing converts;

Fig. 4 is seismic section;

Fig. 5 is the time frequency analysis of seismic trace;

(a) seismic channel data; (b) enlarge section; (c) wavelet transformation; D () Synchrosqueezing converts;

Fig. 6 is the 30Hz frequency slice schematic diagram of seismic section in Fig. 2;

(a) wavelet transformation frequency slice; B () Synchrosqueezing conversion frequency is cut into slices;

Fig. 7 is the horizon slice schematic diagram of 30Hz data volume;

(a) enlarge section horizon slice; B () Synchrosqueezing converts horizon slice;

Fig. 8 is that the seismic attenuation of certain Sandstone Gas Reservoir 3-D data volume is estimated;

(a) work area base map; (b) zone of interest position; C () converts the section of decay behavior result along zone of interest based on Synchrosqueezing.

Embodiment

Below by way of specific embodiments and the drawings, the present invention program is illustrated:

Used tool Synchrosqueezing conversion of the present invention specifically comprises the following steps:

Step a: continuous wavelet transform

Signal wavelet transformation time domain and frequency domain be expressed as

$W_f(a,b)\begin{array}{c}=\frac{1}{\sqrt{a}}{\&Integral;}_{-\&infin;}^{+\&infin;}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}\\ =\frac{\sqrt{a}}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}F\left(\mathrm{\&omega;}\right){\mathrm{\&Psi;}}^{*}\left(\mathrm{a\&omega;}\right)e^{\mathrm{ib\&omega;}}\mathrm{d\&omega;},\end{array}---\left(13\right)$

Wherein ψ (t) is wavelet, a and b is respectively scale factor and shift factor, and F (ω) and Ψ (ω) converts for f (t) and the Fourier of ψ (t).

Here, suppose that wavelet function does not almost have negative frequency components, namely as ω < 0, Ψ (ω) ≈ 0, the concentration of energy of small echo is at positive frequency ω=ω ₀near.As shown in Fig. 1 (a), with cosine signal h (t) of a single-frequency=Acos (Ω t), wherein A=1, Ω=100 π are example, and the result of its wavelet transformation is

$W_h(a,b)=\frac{A}{2}\sqrt{a}{\mathrm{\&Psi;}}^{*}\left(\mathrm{a\&Omega;}\right)e^{\mathrm{ib\&Omega;}}---\left(14\right)$

Through wavelet transformation, the energy of cosine signal h (t) time m-yardstick plane around yardstick a=ω ₀the straight line that/Ω is corresponding spreads, and is distributed on the yardstick band centered by this straight line, as shown in Fig. 1 (b).Can find out, h (t) was a simple signal originally, owing to being subject to the impact of wavelet function, spread at wavelet transformed domain energy, obtained the distribution of " obfuscation ".

The present invention attempts the impact removing wavelet function, from the result of wavelet transformation by frequency omega " demodulation " out.In other words, the energy that yardstick band distributes is all caused by a single-frequency component, and the present invention is converted by Synchrosqueezing, by the yardstick a=ω of the energy on yardstick band again corresponding to " extruding " to this single-frequency component ₀on the straight line at/Ω place.

The calculating of step b:FM frequency, demodulation frequency

In order to find out frequency content corresponding to wavelet coefficient, the result of wavelet transformation is proceeded as follows:

${\mathrm{\&omega;}}_f(a,b)=\left\{\begin{array}{cc}\frac{-i{\&PartialD;}_b{W}_f(a,b)}{W_f(a,b)}& \left|W_f(a,b)\right|>0,\\ \&infin;& \left|W_f(a,b)\right|=0.\end{array}\right.---\left(15\right)$

For cosine signal h (t), formula (2) is substituted into formula (3), can obtain working as | W _h(a, b) | when ≠ 0,

ω _h(a,b)＝Ω,(16)

That is, through the operation of formula (3), from the result of wavelet transformation, again frequency omega can be obtained.Through this step, the frequency content corresponding to each wavelet coefficient can be obtained.

Step c: time m-scale domain to the mapping of time-frequency domain

After obtaining frequency content corresponding to each wavelet coefficient, time m-scale domain carry out " extruding " operation, carry out energy rearrangement, by the coefficient sets of corresponding same frequency composition altogether, by time m-yardstick Planar Mapping then m-frequency plane, i.e. (a, b) → (ω (a, b), b).

For signal note f _at () is analytic signal corresponding to this signal.When wavelet ψ (t) is for resolving small echo,

$\begin{array}{c}{\&Integral;}_0^{\&infin;}{W}_f(a,b)a^{-3/2}\mathrm{da}=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}{\&Integral;}_0^{\&infin;}f\left(\mathrm{\&xi;}\right){\mathrm{\&Psi;}}^{*}\left(\mathrm{a\&xi;}\right)e^{\mathrm{ib\&xi;}}{a}^{-1}\mathrm{d\&xi;da}\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}F\left(\mathrm{\&xi;}\right)e^{\mathrm{ib\&xi;}}\mathrm{d\&xi;}{\&Integral;}_0^{\&infin;}{\mathrm{\&Psi;}}^{*}\left(\mathrm{a\&xi;}\right)a^{-1}\mathrm{da}\\ =\frac{1}{2}{\&Integral;}_0^{\&infin;}{\mathrm{\&Psi;}}^{*}\left(\mathrm{\&xi;}\right){\mathrm{\&xi;}}^{-1}\mathrm{d\&xi;}\&CenterDot;\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{\&infin;}2F\left(\mathrm{\&xi;}\right)e^{\mathrm{ib\&xi;}}\mathrm{d\&xi;}\\ =\frac{1}{2}{\&Integral;}_0^{\&infin;}{\mathrm{\&Psi;}}^{*}\left(\mathrm{\&xi;}\right){\mathrm{\&xi;}}^{-1}\mathrm{d\&xi;}\&CenterDot;f_a\left(b\right),\end{array}---\left(17\right)$

Order $C_{\mathrm{\&psi;}}=\frac{1}{2}{\&Integral;}_0^{\&infin;}{\mathrm{\&Psi;}}^{*}\left(\mathrm{\&xi;}\right){\mathrm{\&xi;}}^{-1}\mathrm{d\&xi;},$ Can obtain

$f\left(b\right)=\mathrm{Re}\left[C_{\mathrm{\&psi;}}^{-1}{\&Integral;}_0^{\&infin;}{W}_f(a,b)a^{-3/2}\mathrm{da}\right].---\left(18\right)$

In formula (6), Re () represents real part.Thus, can be following several form by Synchrosqueezing transform definition:

(1) conitnuous forms

$T_f(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_f(a,b)=\mathrm{\&omega;},W_f(a,b)\&NotEqual;0\}}{W}_f(a,b)a^{-3/2}\mathrm{da}---\left(19\right)$

By formula (7), time m-scale domain, wavelet coefficient corresponding to all and frequencies omega combines, at time-frequency domain again by the position of energy " extruding " to frequencies omega place.

(2) distribution form

$T_f(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,W_f(a,b)\&NotEqual;0\}}{W}_f(a,b)\mathrm{\&delta;}({\mathrm{\&omega;}}_f(a,b)-\mathrm{\&omega;})a^{-3/2}\mathrm{da},---\left(20\right)$

Formula (8) adopts unit impulse function to represent, and the people such as Synchrosqueezing converts, Daubechies explain δ (ω from the angle of distribution _f(a, b)-ω), this form can be called as distribution form.

(3) approximate form

$S_{f,\mathrm{\&epsiv;}}^{\mathrm{\&delta;}}(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,W_f(a,b)>\mathrm{\&epsiv;}\}}{W}_f(a,b)\frac{1}{\mathrm{\&delta;}}h\left(\frac{\mathrm{\&omega;}-{\mathrm{\&omega;}}_f(a,b)}{\mathrm{\&delta;}}\right)a^{-3/2}\mathrm{da},---\left(21\right)$

In formula (9), the unlimited smooth function that h (t) is local support, namely and ∫ h (t) dt=1.When δ → 0, the limit of this method for expressing unit pulse is similar to and replaces unit impulse function.Notice and work as W _fwhen (a, b) is very little, when using formula (3) calculates frequency content corresponding to wavelet coefficient, there will be the situation of numerical value instability, therefore only exist | W _f(a, b) | calculate during > ε.Formula (9) is called the Synchrosqueezing conversion with thresholding ε and precision δ by the present invention.

(4) discrete form

When carrying out numerical evaluation, need to carry out discrete to the yardstick a in formula (7) and frequencies omega.Yardstick after discretize is designated as { a _k, wherein a _k> a _k-1, yardstick is spaced apart a _k-a _k-1=(Δ is a) _k.Frequency is divided, is designated as { ω _l, wherein ω _l> ω _l-1, frequency interval is ω _l-ω _l-1=Δ ω.The discrete form of Synchrosqueezing conversion can be expressed as

$T_f({\mathrm{\&omega;}}_l,b)=\underset{a_k:|{\mathrm{\&omega;}}_f(a_k,b)-{\mathrm{\&omega;}}_l|\&le;\mathrm{\&Delta;\&omega;}/2}{\mathrm{\&Sigma;}}{W}_f(a_k,b){a_k}^{-3/2}{\left(\mathrm{\&Delta;a}\right)}_k,---\left(22\right)$

Visible, by formula (10), corresponding frequency content is at ω _lthe wavelet coefficient of neighbouring (part in Fig. 2 between dotted line) is reconfigured, and energy is " squeezed. " ω in T/F plane _lposition.

If Fig. 1 is that the wavelet transformation of cosine signal and Synchrosqueezing convert, can finds out that the latter has better time-frequency locality, feature the time-frequency distributions of 50Hz cosine signal more exactly.

If Fig. 2 is that Synchrosqueezing converts schematic diagram, the wavelet coefficient of respective frequencies between dotted line is reconfigured.

Fig. 1 (c) depicts the Synchrosqueezing transformation results of cosine signal h (t), and its energy reassembles on the straight line at ω=Ω place at time-frequency plane as seen.

Based on the seismic attenuation method of estimation of Synchrosqueezing conversion

When seismic wavelet is through Isotropic Decay medium, its spectral amplitude is

S(f)＝S ₀(f)e ^-αz(23)

Wherein S ₀f spectral amplitude that () is seismic wavelet, z is propagation distance, and α is attenuation coefficient, and S (f) is the later spectral amplitude of seismic wavelet propagation distance z.Aki and Richards supposes that quality factor q is not with frequency change, gives the relation between attenuation coefficient and quality factor:

$\mathrm{\&alpha;}=\frac{\mathrm{\&pi;f}}{\mathrm{QV}},---\left(24\right)$

In formula, V is phase velocity, and f is frequency.(12) are substituted into (11) obtain

$S\left(f\right)=S_0\left(f\right)e^{-\frac{\mathrm{\&pi;f}}{\mathrm{QV}}z}.---\left(25\right)$

Visible, the high fdrequency component attenuation ratio low frequency component in seismic wavelet wants fast, when seismic event is through the lower hydrocarbon reservoir of Q value, decays more obvious.The different attenuation characteristic of high-low frequency weight can be used to refer to hydrocarbon reservoir, and in view of the good nature of Synchrosqueezing conversion, the present invention is used for portraying the decay of seismic event through tight sand gas-bearing reservoir, carries out the instruction of hydrocarbon reservoir.

The attenuation characteristic difference of research high-low frequency weight, needs to select suitable high frequency f _hwith low frequency f _l, first spectrum analysis is done to the geological data near zone of interest, then selects f _hand f _l, make it meet the following conditions:

(1) high frequency f on spectral amplitude _hamplitude and low frequency f _lamplitude should be roughly the same;

(2) high frequency f _hwith low frequency f _lbetween frequency interval should be enough large, to ensure that high fdrequency component and low frequency component have enough difference in attenuation;

(3) high frequency f _hwith low frequency f _lall should in the frequency band range of seismic wavelet.

Note destination layer position is H _t(x, y), the high fdrequency component utilizing Synchrosqueezing to convert to obtain is T (x, y, t, f _h), low frequency component is T (x, y, t, f _l), below the difference of high fdrequency component and low frequency component just can be adopted to portray layer position H _tdecay near (x, y).Due to the complicacy of underground structure and layer position medium, seismic event there occurs decay before arrival destination layer, produce the difference of high-low frequency weight, in order to portray destination layer to neighbouring decay, need above destination layer position, the amplitude difference of high and low frequency to be eliminated in advance, therefore, remember that the layer position above destination layer is H _a(x, y), is defined as follows modifying factor:

$\mathrm{\&alpha;}(x,y)=\frac{T(x,y,H_A(x,y),f_L)}{T(x,y,H_A(x,y),f_H)}.---\left(26\right)$

Due to the complicacy of underground medium, α (x, y) change is violent, easily makes result of calculation unstable, therefore needs in actual applications the smoothing operation of α (x, y).

The destination layer position H obtained is converted by Synchrosqueezing _tdecay near (x, y) is defined as

AS(x,y,t)＝T(x,y,t,f _L)-α(x,y)T(x,y,t,f _H) (27)。

Material base of the present invention is seismic data volume, employing by road treating method.

The present invention is based on Synchrosqueezing conversion to carry out seismic data time frequency analysis concrete steps and be:

Step 1: choose typical road in 3-D data volume seismic section, carries out Synchrosqueezing conversion to this track data, finds out abnormal area respective frequencies;

Step 2: carry out in the result of Synchrosqueezing conversion to whole seismic section, extracts frequency slice;

Step 3: carry out Synchrosqueezing conversion to whole 3-D data volume, obtains frequency data body, then extracts a horizon slice.

The present invention carry that to carry out the method summary of seismic attenuation estimation based on Synchrosqueezing conversion as follows:

Step 1: determine destination layer scope, carries out spectrum analysis to the 3-d seismic data set near destination layer, chooses suitable high frequency f _hwith low frequency f _l;

Step 2: utilize Synchrosqueezing to convert and obtain the single-frequency data volume of high frequency for T (x, y, t, f _h) and the single-frequency data volume of low frequency be T (x, y, t, f _l);

Step 3: the layer position H above destination layer _anear (x, y), utilize formula (14) to calculate modifying factor α (x, y) and do level and smooth;

Step 4: by the decay AS (x, y, t) near formula (16) estimating target layer.

Effect analysis

One, Numerical Simulation Results

First, with enlarge section and wavelet transformation, the effect of time frequency analysis is carried out to composite signal by contrast Synchrosqueezing conversion.Composite signal s (t) is as shown in Fig. 3 (a), and it is made up of following two components:

s ₁(t)＝sin(3(140πt+30sin(3πt))),(28)

s ₂(t)＝sin(3(80πt+20sin(3πt))).(29)

The present invention adopts enlarge section, wavelet transformation and Synchrosqueezing to convert, as shown in Fig. 3 (b)-(d) to s (t) respectively.Wherein enlarge section adopts the Hamming window of 128, and wavelet transformation adopts the Morlet small echo of σ=6.Although the result of three kinds of conversion can both distinguish two components, the resolution of enlarge section and wavelet transformation is starkly lower than Synchrosqueezing conversion.Due to the impact of window function or wavelet function, the energy of signal spreads at transform domain, makes to occur between two components overlapping.And the result of Synchrosqueezing conversion can distinguish two component of signals, accurate description frequency rule over time well.

As Fig. 3, test signal and the result adopting different Time-Frequency Analysis Method, with traditional enlarge section compared with wavelet transformation, Synchrosqueezing conversion can obtain meticulousr time-frequency distributions, clearly distinguish two components, depict the frequency process over time of each component exactly

Two, actual seismic data

Below, different time frequency analyzing tool is used for actual seismic data.Fig. 4 illustrates the seismic section of a 3-D data volume.This section has 400 roads, per pass 500 sampled points, and sampling interval is 2ms, wherein between the 60th road to the 120th road, and about 1.22s, and between the 255th road to the 305th road, about 1.25s has typical river feature, respectively with oval mark.

As shown in Figure 4, seismic section, ellipse indicates the position in river course.

Carry out time frequency analysis by extracting the 90th road, this road is through the river course on the left side.This track data and enlarge section, wavelet transformation and Synchrosqueezing conversion thereof are as shown in Figure 5.Therefrom can find out the frequency decay trend that the absorption due to stratum is brought.The abnormal area being wherein positioned at the 30Hz of about 1.22s correspond to the position in this river course.Due to the energy dispersal of enlarge section and wavelet transformation, its transform domain spectrogram is comparatively fuzzy, conceals the frequency change feature of unlike signal component.And Synchrosqueezing conversion creates a more sparse time-frequency distributions, there is higher time frequency resolution, disclose the localized variation feature of signal frequency, indicate the abnormal area corresponding with reservoir.

If Fig. 5 is the time frequency analysis of seismic trace.The result of Synchrosqueezing conversion has more sparse distribution and higher time frequency resolution, reflects the localized variation feature of the frequency of the different component of signal, clearly indicates the abnormal area of about the 1.25s 30Hz relevant with river course.

Next, time frequency analysis is carried out to whole seismic section.The frequency slice of 30Hz is extracted, as shown in Figure 6 from the result that wavelet transformation and Synchrosqueezing convert.The frequency slice of Synchrosqueezing conversion is level and smooth unlike the frequency slice of wavelet transformation, some are had to be similar to the feature of " noise ", but, there is more sparse distribution just because of it, there is higher time frequency resolution, more clearly can portray the feature of geologic body.Indicate the position in river course in figure with ellipse, the result of visible Synchrosqueezing conversion very clearly reflects position and the border in river course.And wavelet transformation is due to the impact of wavelet function, energy dispersal on time-frequency figure causes that the component of different frequency intersects, aliasing, the resolution showing as single-frequency section declines, and in figure, join together in the energy of position, river course and both sides, cannot differentiate the border in river course.

As the 30Hz frequency slice that Fig. 6 is seismic section in Fig. 2.The result of Synchrosqueezing change has higher time frequency resolution, indicates position and the border in river course exactly.

Finally, time frequency analysis is carried out to whole 3-D data volume.Adopt enlarge section and Synchrosqueezing conversion respectively, obtain the frequency data body of 30Hz, then extract a horizon slice, as shown in Figure 7.Can find out that the feature of river course (red arrow indicating section) and tomography (green arrow indicating section) becomes very fuzzy on the former figure, and due to very high time frequency resolution, the result of Synchrosqueezing conversion clearly can reflect the feature of river course and tiny tomography.

As the horizon slice that Fig. 7 is 30Hz data volume, wherein red arrow and green arrow indicate the position of river course and tomography respectively, and the result of Synchrosqueezing conversion more clearly reflects the feature in tomography and river course.

The seismic attenuation method of estimation based on Synchrosqueezing conversion the present invention proposed is used for the decay behavior of the three-dimensional real data of certain tight sand.This Sandstone Gas Reservoir main manifestations is the Sandbody Reservoirs feature that factor of porosity is less than the low-permeability of 10%, and net sandstone distribution compares dispersion, small scale, poor continuity, and wave impedance difference is little, brings difficulty to reservoir prediction.This work area base map is as shown in Fig. 8 (a), and destination layer position is as shown in (b).This areal distribution six mouthfuls of wells, represent with WELL1 to WELL6 respectively from north orientation south, wherein the WELL3 of northeastward and WELL4, WELL5, WELL6 of southwestward is I class and II class gas-producing well, I class is prolific well, II class is taken second place, in the drawings with blue mark, northwest to WELL1 and WELL2 be III class well, be dry-well or gas production rate is extremely low, in the drawings with yellow mark.

By the seismic attenuation method of estimation that will convert based on the Synchrosqueezing seismic data for the treatment of this region.After spectrum analysis, select 10Hz as low frequency respectively, 40Hz, as high frequency, after utilizing Synchrosqueezing conversion to obtain attenuation results, cuts into slices, as shown in (c) along the zone of interest shown in (b).Show overdamp characteristic near visible WELL3, WELL4, WELL5 and WELL6, this is that I class and II class gas-producing well meet completely with it, and near WELL1 and WELL2, substantially do not observe obvious decay, this is III class wells with them is identical.In a word, the decay behavior result converted based on Synchrosqueezing and drilling well result have reasonable consistance, as the indicator of hydrocarbon reservoir, geological personnel can be helped to carry out well location and determine to estimate with reservoir gas-bearing amount under certain condition.

As the seismic attenuation drawing for estimate that Fig. 8 is certain Sandstone Gas Reservoir 3-D data volume.

Performance summary of the present invention:

(1) the present invention proposes the method for carrying out seismic signal time-frequency analysis based on Synchrosqueezing conversion, the example of test signal shows, with traditional time frequency analyzing tool, as enlarge section, wavelet transformation etc. are compared, Synchrosqueezing conversion is by the rearrangement to coefficient in transform domain, the time-frequency distributions that can be more assembled, depicts the frequency process over time of each component exactly.

(2) on time-frequency figure, Synchrosqueezing conversion creates a more sparse time-frequency distributions, has higher time frequency resolution, discloses the localized variation feature of signal frequency, indicate the abnormal area corresponding with reservoir;

(3) in the single-frequency section of seismic section, the result of Synchrosqueezing conversion more clearly can reflect position and the border in river course;

(4) on the horizon slice of data volume, the result of Synchrosqueezing conversion more clearly reflects the feature in tomography and river course;

(5) the present invention proposes the seismic attenuation method of estimation based on Synchrosqueezing conversion, provide its realization flow, and the method is used for the decay behavior of certain oil field Sandstone Gas Reservoir 3-D data volume, estimated result and drilling well result have good consistance, can as the direct indicator of gas-bearing reservoir, carry out reservoir gas-bearing amount to estimate and well location is determined with helping geological personnel.

Claims (5)

1., based on the seismic data Time-Frequency Analysis Method of Synchrosqueezing conversion, it is characterized in that:

Synchrosqueezing conversion specifically comprises the following steps:

Step a: continuous wavelet transform

Signal wavelet transformation time domain and frequency domain be expressed as:

$\begin{array}{c}{W}_f(a,b)=\frac{1}{\sqrt{a}}{\&Integral;}_{-\&infin;}^{+\&infin;}f\left(t\right)\mathrm{\&psi;}*\left(\frac{t-b}{a}\right)\mathrm{dt}\\ =\frac{\sqrt{a}}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}F\left(\mathrm{\&omega;}\right){\mathrm{\&psi;}}^{*}\left(\mathrm{a\&omega;}\right)e^{\mathrm{ib\&omega;}}\mathrm{d\&omega;}\end{array},---\left(1\right)$

Wherein ψ (t) is wavelet, a and b is respectively scale factor and shift factor, and F (ω) and Ψ (ω) converts for f (t) and the Fourier of ψ (t); Suppose that wavelet function does not almost have negative frequency components, namely as ω < 0, Ψ (ω) ≈ 0;

The calculating of step b:FM frequency, demodulation frequency

${\mathrm{\&omega;}}_f(a,b)=\left\{\begin{array}{cc}\frac{-i{\&PartialD;}_b{W}_f(a,b)}{W_f(a,b)}& \left|W_f(a,b)\right|>0,\\ \&infin;& \left|W_f(a,b)\right|=0.\end{array}\right.---\left(2\right);$

Step c: time m-scale domain to the mapping of time-frequency domain

(1) conitnuous forms

$T_f(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_f(a,b)=\mathrm{\&omega;},W_f(a,b)\&NotEqual;0\}}{W}_f(a,b)a^{-3/2}\mathrm{da}---\left(3\right)$

By formula (7), time wavelet coefficient corresponding to m-scale domain all and frequencies omega combine, at time-frequency domain again by the position of energy " extruding " to frequencies omega place;

(2) discrete form

When carrying out numerical evaluation, need to carry out discrete to the yardstick a in formula (7) and frequencies omega; Yardstick after discretize is designated as { a _k, wherein a _k> a _k-1, yardstick is spaced apart a _k-a _k-1=(Δ is a) _k; Frequency is divided, is designated as { ω _l, wherein ω _l> ω _l-1, frequency interval is ω _l-ω _l-1=Δ ω; The discrete form of Synchrosqueezing conversion is expressed as:

$T_f({\mathrm{\&omega;}}_l,b)=\underset{a_k:|{\mathrm{\&omega;}}_f(a_k,b)-{\mathrm{\&omega;}}_l|\&le;\mathrm{\&Delta;\&omega;}/2}{\mathrm{\&Sigma;}}{W}_f(a_k,b)a_k^{-3/2}{\left(\mathrm{\&Delta;a}\right)}_k,---\left(4\right);$

Seismic data Time-Frequency Analysis Method specifically comprises the following steps:

Step 1: choose typical road in 3-D data volume seismic section, carries out Synchrosqueezing conversion to this track data, finds out abnormal area respective frequencies;

Step 2: carry out Synchrosqueezing conversion to whole seismic section, extracts the section of abnormal area respective frequencies;

Step 3: carry out Synchrosqueezing conversion to whole 3-D data volume, obtains frequency data body, then extracts a horizon slice, carries out seismic data interpretation for geological personnel.

2. carry out the method for seismic attenuation estimation based on Synchrosqueezing conversion, it is characterized in that comprising the following steps:

Step 1: determine destination layer scope, carries out spectrum analysis to the 3-d seismic data set near destination layer, chooses suitable high frequency f _hwith low frequency f _l;

Step 2: utilize Synchrosqueezing to convert the single-frequency data volume T (x, y, t, the f that obtain high frequency _h) and single-frequency data volume T (x, y, t, the f of low frequency _l);

Step 3: the layer position H above destination layer _anear (x, y), the amplitude difference of high and low frequency is eliminated in advance;

Step 4: the decay near estimating target layer.

3. method of carrying out seismic attenuation estimation based on Synchrosqueezing conversion according to claim 2, is characterized in that: in described step 1, high frequency f _hamplitude and low frequency f _lamplitude roughly the same.

4. method of carrying out seismic attenuation estimation based on Synchrosqueezing conversion according to claim 2, is characterized in that: in described step 3, utilizes formula (14) to calculate modifying factor α (x, y) and does level and smooth;

$\mathrm{\&alpha;}(x,y)=\frac{T(x,y,H_A(x,y),f_L)}{T(x,y,H_A(x,y),f_H)}.---\left(5\right).$

5. method of carrying out seismic attenuation estimation based on Synchrosqueezing conversion according to claim 2, is characterized in that: in described step 4, by the decay AS (x, y, t) near formula (15) estimating target layer;

AS(x,y,t)＝T(x,y,t,f _L)-α(x,y)T(x,y,t,f _H) (6)。