CN104880731A - Seismic instantaneous attribute extraction method based on generalized Morse frame - Google Patents

Abstract

The invention discloses a seismic instantaneous attribute extraction method based on a generalized Morse frame. The generalized Morse frame is used for converting the determination of an effective signal energy distribution space into an optimization problem, which can be solved by using an iterative contraction threshold algorithm and a fast iterative shrink threshold algorithm. On the basis that the effective signal energy distribution space is determined, a method for analyzing noisy signal instantaneous attributes is provided by using the relationship between wavelet transform and Hilbert transform. The method provided by the invention has good anti-noise performance and high accuracy, and the obtained instantaneous frequency cross section can more clearly reflect the change in the instantaneous frequency and indicate an abnormal area.

Description

Based on the earthquake instantaneous attribute extracting method of broad sense Morse frame

Technical field

The invention belongs to the signal transacting field in geophysical survey, relate to the instantaneous attribute extracting method of seismic data, particularly relate to a kind of earthquake instantaneous attribute extracting method based on broad sense Morse frame.

Background technology

Seismic properties refers to the geometric shape of the relevant seismic event of being derived through mathematic(al) manipulation by geological data, kinematics character, dynamic characteristic and statistics feature, and these features visually can reflect form and the oil-gas possibility thereof of reservoir.20 century 70s, seismic attributes analysis starts to be introduced in seismic interpretation.Since the nineties in 20th century, due to the needs that layer description and 3-D data volume are explained, Seismic attribute analysis technology develops rapidly, and seismic properties is widely applied in stratigraphic structure explanation, reservoir lithology and physical property characteristic description, prediction of oil-gas reserve and dynamic surveillance.Go back the seismic properties classification that neither one is generally acknowledged at present, according to physical significance, seismic properties can be divided into several large classes such as time, amplitude, frequency, relevant, decay, according to attribute pick-up method, and window attribute when window attribute, multiple tracks when every class can be divided into again instantaneous attribute, single track.Wherein instantaneous attribute is the attribute picked up on the position that seismic event arrives, and comprises instantaneous amplitude, instantaneous phase, instantaneous frequency, instant bandwidth etc.In recent years, the quantity of seismic properties is uprushed, but instantaneous attribute remains the pillar of geological data geologic interpretation, is also the technology modules of current commercial processes software indispensability.

In signal transacting field, the concept of the signal transient attribute being representative with instantaneous amplitude, instantaneous phase, instantaneous frequency is long-standing, after the work of Gabor and Ville, a large amount of scholar is studied in association area, and Boashash summarizes these work.The people such as Taner introduced complex seismic trace in 1979, proposed the method being asked for instantaneous attribute by complex seismic trace, and gave the physical significance of these attributes and the application in seismic interpretation.Lithological change and the oil-gas accumulation of instantaneous amplitude and adjacent layer are relevant.The instantaneous phase reflection uncontinuity at interface, tomography, unconformity surface and sequence boundaries etc.The change of instantaneous frequency effectively can portray thickness and the lithological change on stratum, the distribution etc. of Indication of Oil-Gas.Instantaneous attribute is used for hole sandstone thickness of thin layer and estimates by Robertson and Nogami, Chopra and Marfurt utilizes instantaneous attribute to carry out uncontinuity, tomography and horizontal discontinuity detecting technique.Liu and Marfurt utilizes instantaneous attribute detect and portray the distribution of meandering river and determine its thickness.Zeng utilizes instantaneous frequency extremely to indicate thin layer.The people such as Gao utilize instantaneous frequency to carry out the estimation of seismic data Q value.

The method of conventional extraction earthquake instantaneous attribute can be divided into following three classes:

(1) Hilbert transform method.Do Hilbert conversion to seismic trace, be converted into complex seismic trace, wherein real part is former seismic channel data, and imaginary part converts for its Hilbert.After obtaining complex seismic trace, just can at attributes such as each sampled point calculated amplitude, phase place and frequencies, i.e. instantaneous attribute.After the work of Taner, Hilbert converter technique is widely used in calculating earthquake instantaneous attribute, and most of business software still adopts the method so far.A lot of scholar is developed the method and is improved, and Barnes proposed the concept of two-dimentional complex seismic trace in 1996.The people such as Luo proposed generalized Hilbert transformation in 2003 and provide its application in geophysics.Barnes proposed the concept of weighting instantaneous frequency in 2007.Hilbert transform method is promoted by Lu and Zhang, proposes the method calculating instantaneous frequency based on sef-adapting filter.

(2) based on the method for time frequency analysis.Time-Frequency Analysis Method utilizes the time-frequency distributions of signal to ask for instantaneous attribute.The people such as Boashash propose the self-adaptation instantaneous Frequency Estimation method based on time frequency analysis.The people such as Stankovic utilize adaptive windows time-frequency distributions to calculate instantaneous frequency.The people such as Gao Jinghuai propose the method calculating earthquake instantaneous attribute in phase space.The people such as Marfurt propose the seismic data instantaneous attribute extracting method analyzed based on narrow-band spectrum.Steeghs and Drijkoningen proposes based on the seismic sequence analysis of Bilinear TFD and attributes extraction method.The people such as Huang propose the method calculating instantaneous frequency based on empirical mode decomposition (EMD).The people such as Han by complete overall experience Mode Decomposition (CEEMD) method for extracting the instantaneous frequency of seismic data.

(3) based on the method for inverting.The people such as Fomel proposes to utilize the method for inverting to obtain local attribute, and compared to instantaneous frequency, definitely, and effect is obvious for local frequencies physical significance.The people such as Liu propose the method based on inverting, calculate the instantaneous frequency of seismic data.

But based on Hilbert conversion, conventional estimates that the method for instantaneous parameters is very sensitive to noise, and due to the truncation effect of wave filter, make the instantaneous attribute precision that calculates low.

Summary of the invention

The object of the invention is to overcome the deficiencies in the prior art, provide a kind of earthquake instantaneous attribute extracting method based on broad sense Morse frame, have good noiseproof feature and accuracy, the instantaneous frequency section obtained clearly can reflect the change of instantaneous frequency.

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

Based on an earthquake instantaneous attribute extracting method for broad sense Morse frame, comprise the following steps:

Step 1: acquiring seismic data body;

Step 2: calculate the broad sense Morse frame coefficient that the noisy geological data in each road is corresponding;

The wavelet transformation of real seismic trace s (t) is defined as:

$W_{\mathrm{\ψ}}(t,a)=\frac{1}{a}{\∫}_{-\∞}^{\∞}s\left(\mathrm{\τ}\right){\mathrm{\ψ}}^{*}\left(\frac{\mathrm{\τ}-t}{a}\right)\mathrm{d\τ},---\left(22\right)$

In formula, wavelet ψ (t) selects broad sense Morse small echo;

Broad sense Morse wavelet frequency domain expression formula is:

${\mathrm{\Ψ}}_{\mathrm{\β},\mathrm{\γ}}\left(\mathrm{\ω}\right)=U\left(\mathrm{\ω}\right){\mathrm{\α}}_{\mathrm{\β},\mathrm{\γ}}{\mathrm{\ω}}^{\mathrm{\β}}{e}^{-{\mathrm{\ω}}^{\mathrm{\γ}}},---\left(23\right)$

Wherein U (ω) is unit step function, β and γ is the parameter of small echo, and β > 0, γ > 0, α _{beta, gamma}for normalization constant, α _{beta, gamma}≡ 2 (e γ/β) ^β/γ;

Little wave system after discretize can be expressed as:

${\mathrm{\ψ}}_{m,n}\left(x\right)=a_0^{-m}\mathrm{\ψ}\left[a_0^{-m}(x-n{t}_0)\right]---\left(24\right)$

Wherein a ₀the discretize step-length of scale factor a and a ₀> 1, t ₀for the discretize step-length of shift factor t;

With the wavelet representation after the discretize that formula (1) is corresponding be:

C _m,n＝<s,ψ _m,n>, (25)

Wherein C _m,nfor frame coefficient;

Step 3: the broad sense Morse frame coefficient iteration for each road obtains coefficient corresponding to useful signal;

Definition K be from be mapped to operator, by frame coefficient be mapped as on signal, namely

The adjoint operator K of operator K ^*for from be mapped to operator, will on signal project on frame coefficient:

K ^*s＝<s,ψ′ _m,n>, (27)

Then s=KK ^*s, (28)

Operator K is composite operator, K ^*for analyze operator and Wavelet Frames ψ ' _m,ncorresponding;

Signals and associated noises is expressed as y=s+n=Kx+n, (29)

Wherein y represents signals and associated noises, and s is not noisy useful signal, and n is white Gaussian noise, and K is the composite operator in formula (22), and x represents the coefficient of transform domain, obtains coefficient corresponding to useful signal by solving following optimization problem

$\tilde{x}=\arg \underset{x}{\min }{\left|\right|x\left|\right|}_{1}s.t.\left|\right|y-\mathrm{Kx}{\left|\right|}_2\&le;\mathrm{\&epsiv;},---\left(30\right)$

Wherein ε is relevant with the noise level of signal to be analyzed, with Lagrange multiplier λ, this problem is converted into following unconstrained problem:

${\tilde{x}}_{\mathrm{\&lambda;}}=\arg \underset{x}{\min }\frac{1}{2}{\left|\right|y-\mathrm{Kx}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x\left|\right|}_1,---\left(31\right)$

Wherein λ is called as regularization parameter, adopts the alternative manner shown in following formula to upgrade coefficient x

x ^(k+1)＝T _λ[x ^(k)+K ^*(y-Kx ^(k))],k＝1,2,…,N,(32)

Wherein T _λfor threshold function table, be defined as

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&lambda;}\frac{x}{\left|x\right|}& \left|x\right|\&GreaterEqual;\mathrm{\&lambda;},\\ 0& \left|x\right|<\mathrm{\&lambda;}.\end{array}\right.---\left(33\right);$

Step 4: calculate analytic signal by formula (2);

$c\left(t\right)\&equiv;\frac{1}{K_{\mathrm{\&psi;}}}{\&Integral;}_0^{\&infin;}{W}_{\mathrm{\&psi;}}(t,a)\frac{\mathrm{da}}{a}=s\left(t\right)+\mathrm{ih}\left(t\right),---\left(34\right)$

Wherein W _ψ(t, a) is the coefficient that useful signal is corresponding, and h (t) is the Hilbert conversion of s (t), and c (t) is the analytic signal that s (t) is corresponding;

Step 5: utilize c (t) to calculate instantaneous amplitude, instantaneous phase and instantaneous frequency:

$\hat{e}\left(t\right)=\sqrt{{\mathrm{Re}}^2\left[c\left(t\right)\right]+{\mathrm{Im}}^2\left[c\left(t\right)\right]},---\left(35\right)$

$\widehat{\mathrm{\&theta;}}\left(t\right)=\arctan \left[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\right],---\left(36\right)$

$\hat{f}\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d}{\mathrm{dt}}\left[\arctan \right[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\left]\right],---\left(37\right)$

Wherein Re [c (t)] and Im [c (t)] represents real part and the imaginary part of c (t) respectively.

As the further preferred version of the present invention, described step 3 formula (26) is each time in iterative process, adopts index threshold to decline strategy, as shown in formula (28):

${\mathrm{\&lambda;}}_k={\left(\frac{{\mathrm{\&lambda;}}_{\min }}{{\mathrm{\&lambda;}}_{\max }}\right)}^{(k-1)/(N-1)}\&CenterDot;{\mathrm{\&lambda;}}_{\max },k=\mathrm{1,2},...,N.---\left(38\right)$

Wherein

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\max }=p_{\max }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\\ {\mathrm{\&lambda;}}_{\min }=p_{\min }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\end{array}\right.---\left(39\right)$

λ in formula (29) _maxand λ _minbe respectively minimum value and the maximal value of regularization parameter, p _maxand p _minfor minimum and maximum number percent.

As the further preferred version of the present invention, in order to stop iterative process in time in described step 3, define a dynamic stopping criterion

$\frac{{\left|\right|x^{(k+1)}-x^{\left(k\right)}\left|\right|}^2}{{\left|\right|x^{(k+1)}\left|\right|}^2}<\mathrm{tolerance},---\left(40\right)$

Wherein tolerance is allowable value given before iteration, and under this criterion, when continuing iteration and can not obtaining better result, iterative process just stops automatically.

As the further preferred version of the present invention, in order to improve speed of convergence in described step 3, iteratively faster atrophy thresholding algorithm is adopted to see formula (31);

$\left\{\begin{array}{c}{\tilde{x}}^{\left(k\right)}=x^{\left(k\right)}+\frac{t^{\left(k\right)}-1}{t^{(k+1)}}(x^{\left(k\right)}-x^{(k-1)}),\\ x^{(k+1)}=T_{\mathrm{\&lambda;}}[{\tilde{x}}^{\left(k\right)}+K^{*}(y-K{\tilde{x}}^{\left(k\right)})],\end{array}\right.---\left(41\right)$

Wherein x ⁽⁰⁾=0, t ⁽⁰⁾=1, t ^(k)meet following recursion formula:

$t^{(k+1)}=(1+\sqrt{1+4{\left(t^{\left(k\right)}\right)}^2})/2---\left(42\right).$

The present invention is based on the earthquake instantaneous attribute extracting method of broad sense Morse frame, adopt broad sense Morse frame, the determination in useful signal energy distribution space is changed into an optimization problem, utilizes iteration atrophy thresholding algorithm and iteratively faster atrophy thresholding algorithm to solve.On the basis determining useful signal energy distribution space, the relation utilizing wavelet transformation and Hilbert to convert, proposes the method for signals and associated noises instantaneous attribution analysis.Method of the present invention has good noiseproof feature and accuracy, and the instantaneous frequency section obtained clearly can reflect the change of instantaneous frequency, instruction abnormal area.

Further, in order to improve convergence of algorithm speed, the present invention adopts index threshold strategy and dynamic stopping criterion, makes by less iterations, just can in the hope of coefficient corresponding to useful signal.

Accompanying drawing explanation

The time scale spectrogram of the Ricker wavelet of Figure 150 Hz;

The wavelet transformation spectrogram of (a) not noisy ricker wavelet

The wavelet transformation spectrogram of (b) noisy ricker wavelet

C spectrogram that the coefficient after () iteration is formed

The instantaneous frequency that Fig. 2 test signal and distinct methods calculate;

The instantaneous frequency that the noisy test signal of Fig. 3 and distinct methods calculate;

The instantaneous frequency of the noisy test signal that Fig. 4 distinct methods calculates;

(a) ricker wavelet instantaneous frequency

(b) theogram instantaneous frequency

(c) real data instantaneous frequency

Fig. 5 is the instantaneous frequency signal to noise ratio (S/N ratio) of three test signals;

The FSNR of (a) ricker wavelet

The FSNR of (b) theogram

The FSNR in (c) actual seismic road

Fig. 6 is real data example;

(a) noisy seismic section

B method instantaneous frequency section that () adopts Hibert to obtain

C noisy seismic section that () adopts the present invention to obtain

Fig. 7 is the instantaneous frequency horizon slice of certain three dimensional seismic data

A instantaneous frequency horizon slice that () utilizes business software to calculate

B instantaneous frequency horizon slice that () calculates based on the inventive method

The instant bandwidth horizon slice of Fig. 8 three dimensional seismic data

A instant bandwidth horizon slice that () utilizes business software to calculate

B instant bandwidth horizon slice that () calculates based on the inventive method

The bound of the frame that the little wave system after table 1 discretize is formed and comfort level

Table 2 asks for the flow process (false code) of useful signal coefficient of correspondence by iteration

The instantaneous frequency signal to noise ratio (S/N ratio) FSNR of table 3 distinct methods.

Embodiment

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

Wavelet Transform calculates instantaneous attribute

The wavelet transformation of real seismic trace s (t) is defined as

$W_{\mathrm{\&psi;}}(t,a)=\frac{1}{a}{\&Integral;}_{-\&infin;}^{\&infin;}s\left(\mathrm{\&tau;}\right){\mathrm{\&psi;}}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)\mathrm{d\&tau;}.---\left(43\right)$

When wavelet ψ (t) is for resolving small echo, Gao Jinghuai etc. demonstrate following theorem:

If ψ (t) is for resolving small echo, its real part ψ _rt () is for even function and its Fourier transform Ψ _r(ω) following condition is met: then to arbitrary signal following equalities is had to set up

$c\left(t\right)\&equiv;\frac{1}{K_{\mathrm{\&psi;}}}{\&Integral;}_0^{\&infin;}{W}_{\mathrm{\&psi;}}(t,a)\frac{\mathrm{da}}{a}=s\left(t\right)+\mathrm{ih}\left(t\right),---\left(44\right)$

Wherein h (t) is the Hilbert conversion of s (t), and c (t) is the analytic signal that s (t) is corresponding.This theorem establishes the relation between wavelet transformation and Hilbert conversion.Therefore, the present invention utilizes c (t) to calculate instantaneous amplitude, instantaneous phase and instantaneous frequency:

$\hat{e}\left(t\right)=\sqrt{{\mathrm{Re}}^2\left[c\left(t\right)\right]+{\mathrm{Im}}^2\left[c\left(t\right)\right]},---\left(45\right)$

$\widehat{\mathrm{\&theta;}}\left(t\right)=\arctan \left[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\right],---\left(46\right)$

$\hat{f}\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d}{\mathrm{dt}}\left[\arctan \right[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\left]\right],---\left(47\right)$

Wherein Re [c (t)] and Im [c (t)] represents real part and the imaginary part of c (t) respectively.

Choosing of wavelet function

Broad sense Morse small echo is the little wave system of the two-parameter parsing of a class, derives when studying time frequency localization operator by people such as Daubechies the earliest, is the solution of time-frequency combination localization problem.Its frequency-domain expression is

${\mathrm{\&Psi;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}\left(\mathrm{\&omega;}\right)=U\left(\mathrm{\&omega;}\right){\mathrm{\&alpha;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}{\mathrm{\&omega;}}^{\mathrm{\&beta;}}{e}^{-{\mathrm{\&omega;}}^{\mathrm{\&gamma;}}},---\left(48\right)$

Wherein U (ω) is unit step function, β and γ is the parameter of small echo, and β > 0, γ > 0.α _{beta, gamma}for normalization constant:

α _β,γ≡2(eγ/β) ^β/γ.(49)

The present invention selects broad sense Morse small echo to calculate the instantaneous attribute of signal, and main cause has:

(1) broad sense Morse small echo is rigorous analytic, and the conventional Morlet small echo just approximate analysis when modulating frequency is enough large, when the time, localization requirement was higher, Morlet small echo there will be negative frequency to be revealed.The more important thing is, above require that wavelet is resolve small echo in theorem, small echo departs from analyticity and this theorem can be caused no longer to set up, the accuracy of the required instantaneous attribute of impact.And broad sense Morse small echo rigorous analytic, meet theorem requirement completely;

(2) broad sense Morse small echo has higher degree of freedom, by adjustment parameter, can show different forms and feature, to mate different seismic signals;

(3) the broad sense Morse small echo after discrete more easily forms tight frame, and the operand of wavelet transformation is reduced greatly, also ensure that below solving-optimizing problem time convergence.

The discretization method of wavelet function

Small echo atom carries out flexible and translation by wavelet and obtains, namely

${\mathrm{\&psi;}}_{t,a}\left(\mathrm{\&tau;}\right)=\frac{1}{a}\mathrm{\&psi;}\left(\frac{\mathrm{\&tau;}-t}{a}\right).---\left(50\right)$

When carrying out numerical evaluation, need to carry out discretize to scale factor a and translation factor t, it is discrete that the present invention carries out exponentiate to scale factor, is designated as wherein a ₀for discretize step-length and a ₀the discretize step-length of > 1, shift factor t is t ₀.Like this, the little wave system after discretize can be expressed as

${\mathrm{\&psi;}}_{m,n}\left(x\right)=a_0^{-m}\mathrm{\&psi;}\left[a_0^{-m}(x-n{t}_0)\right]---\left(51\right)$

In discretization method of the present invention, t ₀identical with the sampling interval of signal to be analyzed, therefore the little wave system of discretize meets translation invariance, can regard as by dictionary translation produces.For the dictionary possessing translation invariance, there is following theorem:

If there is B>=A > 0, make for have

$A\&le;\underset{m}{\mathrm{\&Sigma;}}{\left|{\mathrm{\&Phi;}}_m\left(\mathrm{\&omega;}\right)\right|}^2\&le;B,---\left(52\right)$

Then this translation invariant dictionary forms a frame, A and B is respectively the upper bound and the lower bound of frame.The condition that dictionary is formed frame by this theorem is converted into φ _mthe Fourier transform Φ of (x) _m(ω) condition that should meet.

For wavelet function, based on discretization method of the present invention, have then the condition equivalence of above-mentioned theorem is

$A\&le;\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{\left|{\mathrm{\&Phi;}}_m\left(a_0^m\mathrm{\&omega;}\right)\right|}^2\&le;B.---\left(53\right)$

The present invention adopts comfort level δ to portray the tight degree of frame:

$\mathrm{\&delta;}=\frac{B-A}{B+A}.---\left(54\right)$

Comfort level is less, illustrates that the lower bound A of frame and upper bound B is more close, and frame is more close to tight frame.When δ=0, bound is equal, and frame is tight frame.

Note $F_{a_0}\left(\mathrm{\&omega;}\right)\&equiv;\underset{m=-\&infin;}{\overset{+\&infin;}{\mathrm{\&Sigma;}}}{\left|{\mathrm{\&Psi;}}_m\left(a_0^m\mathrm{\&omega;}\right)\right|}^2,---\left(55\right)$

Then calculate at different Scale Discreteness step-length a ₀under, the bound of the frame that the Morlet small echo after discretize and broad sense Morse small echo (GMW) are formed and comfort level, as shown in table 1.The parameter that the present invention chooses broad sense Morse small echo is β=1, and γ=3, in order to ensure the approximate analysis of Morlet small echo, parameter gets σ=6.As can be seen from the table, based on discretize mode of the present invention, the little wave system of Morlet and the little wave system of broad sense Morse all can form frame, but the comfort level of broad sense Morse small echo is less, more easily forms tight frame, particularly at sampling step length a ₀time larger, the little wave system of broad sense Morse can be seen as tight frame, and this can reduce operand greatly.

The bound of the frame that the little wave system after table 1 discretize is formed and comfort level

For the little wave system that can form tight frame, the wavelet transformation after the discretize corresponding with formula and signal reconstruction can be expressed as:

C _m,n＝<s,ψ _m,n>,(56)

$s=\frac{1}{A}\underset{m,n}{\mathrm{\&Sigma;}}{C}_{m,n}{\mathrm{\&psi;}}_{m,n},---\left(57\right)$

Wherein C _m,nfor frame coefficient, s is the vector that the signal of discretize is formed, { Ψ _m,nbe the tight frame of small echo, <s, ψ _m,n> represents the two inner product, and A is boundary's (upper bound is equal with lower bound) of frame.The little wave system of discretize forms tight frame, calculate wavelet transformation coefficient and become very convenient from wavelet coefficient reconstruction signal, also ensure that equally below solving-optimizing problem time convergence.

The operator representation of Wavelet Frames

Order

ψ′ _m,n＝A ^-1/2ψ _m,n,(58)

Then formula (14) and (15) can be write as

C _m,n＝<s,ψ′ _m,n>,(59)

$s=\underset{m,n}{\mathrm{\&Sigma;}}{C}_{m,n}{\mathrm{\&psi;}}_{m,n}^{\&prime;},---\left(60\right)$

Namely

$s=\underset{m,n}{\mathrm{\&Sigma;}}<s,{\mathrm{\&psi;}}_{m,n}^{\&prime;}>{\mathrm{\&psi;}}_{m,n}^{\&prime;}.---\left(61\right)$

The present invention define K for from be mapped to operator, by frame coefficient be mapped as on signal, namely

The adjoint operator K of operator K ^*for from be mapped to operator, will on signal project on frame coefficient:

K ^*s＝<s,ψ′ _m,n>,(63)

Then formula (19) can be write as

s＝KK ^*s,(64)

Operator K is composite operator, K ^*for analyzing operator, and Wavelet Frames ψ ' _m,ncorresponding.

The determination in useful signal energy distribution space

When selecting suitable wavelet function, wavelet transformation is carried out to signals and associated noises, when being projected to Time-Scale Domain, the energy of useful signal will be distributed in less subspace V, and the energy of noise can be diffused into larger subspace V ', m-scale domain time even whole.In other words, useful signal is corresponding with the coefficient of minority, and noise is almost distributed on whole coefficients.If determine the subspace V that useful signal is corresponding, obtain the coefficient that useful signal is corresponding, by other coefficient zero setting, then noise will be pressed at transform domain, and signal to noise ratio (S/N ratio) is improved.

In order to determine the energy distribution of useful signal, wishing that useful signal projects on the least possible Wavelet Frames atom, that is wishing that signal obtains expression sparse as far as possible at transform domain, so just can suppress more noise, improve signal to noise ratio (S/N ratio).Specifically, a signals and associated noises can be expressed as

y＝s+n＝Kx+n,(65)

Wherein y represents signals and associated noises, and s is not noisy useful signal, and n is white Gaussian noise, and K is the composite operator in formula (22), and x represents the coefficient of transform domain.Wish to obtain an optimum coefficient make (1) sparse as far as possible; (2) value little.Openness and l ₀norm minimum problem is corresponding, is non-convex problem, solves difficulty very large, is usually converted into l under certain condition ₁norm minimum problem.Therefore, can obtain by solving following optimization problem

$\tilde{x}=\arg \underset{x}{\min }{\left|\right|x\left|\right|}_{1}s.t.\left|\right|y-\mathrm{Kx}{\left|\right|}_2\&le;\mathrm{\&epsiv;},---\left(66\right)$

Wherein ε is relevant with the noise level of signal to be analyzed.According to the method that the people such as Elad propose, with Lagrange multiplier λ, this problem can be converted into following unconstrained problem:

${\tilde{x}}_{\mathrm{\&lambda;}}=\arg \underset{x}{\min }\frac{1}{2}{\left|\right|y-\mathrm{Kx}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x\left|\right|}_1,---\left(67\right)$

Wherein λ is called as regularization parameter, so, can directly adopt iteration atrophy thresholding algorithm (IST) to solve this problem.The people such as Daubechies point out, adopt the alternative manner shown in following formula to upgrade coefficient x, and when iterations is enough large time, the result of iteration converges on the solution of optimization problem shown in formula (24).

x ^(k+1)＝T _λ[x ^(k)+K ^*(y-Kx ^(k))],k＝1,2,…,N,(68)

Wherein T _λfor threshold function table, be defined as

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&lambda;}\frac{x}{\left|x\right|}& \left|x\right|\&GreaterEqual;\mathrm{\&lambda;},\\ 0& \left|x\right|<\mathrm{\&lambda;}.\end{array}\right.---\left(69\right)$

It is pointed out that due to what adopt in iteration atrophy thresholding algorithm (IST) it is soft-threshold (as Suo Shi formula (27)), now threshold value is equal with regularization parameter, represents with λ.

In solving-optimizing problem process, need to take threshold value to decline tactful, that is: along with the increase of iterations, the value of λ constantly reduces.Conventional threshold value decline tactful linear decline, index decreased etc., the people such as Gao study discovery, adopt index threshold descending method greatly can improve the speed of algorithm convergence.

In iterative process each time, index threshold is adopted to decline strategy, as shown in formula (28):

${\mathrm{\&lambda;}}_k={\left(\frac{{\mathrm{\&lambda;}}_{\min }}{{\mathrm{\&lambda;}}_{\max }}\right)}^{(k-1)/(N-1)}\&CenterDot;{\mathrm{\&lambda;}}_{\max },k=\mathrm{1,2},...,N.---\left(70\right)$

Wherein

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\max }=p_{\max }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\\ {\mathrm{\&lambda;}}_{\min }=p_{\min }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\end{array}\right.---\left(71\right)$

λ in formula (29) _maxand λ _minbe respectively minimum value and the maximal value of regularization parameter, p _maxand p _minfor minimum and maximum number percent.

Just threshold strategies above only can stop when iterations reaches maximal value N, in order to iterative process can be being stopped in time, invention defines a dynamic stopping criterion

$\frac{{\left|\right|x^{(k+1)}-x^{\left(k\right)}\left|\right|}^2}{{\left|\right|x^{(k+1)}\left|\right|}^2}<\mathrm{tolerance},---\left(72\right)$

Wherein tolerance is allowable value given before iteration.Under this criterion, when continuing iteration and can not obtaining better result, iterative process will stop automatically, to reduce unnecessary iterations.

The speed of convergence of iteration atrophy thresholding algorithm (IST) is very slow, is only O (1/k), namely has single order speed of convergence.Beck and Teboulle proposes iteratively faster atrophy thresholding algorithm (FIST), sees formula (31), and the method makes algorithm the convergence speed rise to O (1/k ²), namely there is second order convergence speed.

$\left\{\begin{array}{c}{\tilde{x}}^{\left(k\right)}=x^{\left(k\right)}+\frac{t^{\left(k\right)}-1}{t^{(k+1)}}(x^{\left(k\right)}-x^{(k-1)}),\\ x^{(k+1)}=T_{\mathrm{\&lambda;}}[{\tilde{x}}^{\left(k\right)}+K^{*}(y-K{\tilde{x}}^{\left(k\right)})],\end{array}\right.---\left(73\right)$

Wherein x ⁽⁰⁾=0, t ⁽⁰⁾=1, t ^(k)meet following recursion formula:

$t^{(k+1)}=(1+\sqrt{1+4{\left(t^{\left(k\right)}\right)}^2})/2.---\left(74\right)$

Table 2 gives the flow process asking for useful signal coefficient of correspondence (false code), and by iteration, the present invention obtains the solution of optimization problem, the coefficient that the signal after being namely pressed with noise is corresponding fig. 1 (a)-(b) depicts the transform domain spectrogram that not noisy and signal to noise ratio (S/N ratio) is the Ricker wavelet of the 50Hz of 5dB respectively, can see the impact that the distribution of noise on useful signal causes.Fig. 1 (c) depicts the spectrogram of the coefficient formation obtained after iteration, and in figure, noise is effectively neutralized, and accurately features the energy distribution of useful signal.

Table 2 asks for the flow process (false code) of useful signal coefficient of correspondence by iteration

The time scale spectrogram of the Ricker wavelet of Figure 150 Hz, can find out that useful signal energy distribution is subject to the impact of noise in the wavelet transformation spectrogram of signals and associated noises, after iteration, noise is suppressed, and spectrogram clearly reflects the energy distribution of useful signal, has had coefficient afterwards, the present invention can obtain the analytic signal corresponding with useful signal, calculates instantaneous attribute further.

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

Step 1: calculate broad sense Morse frame coefficient corresponding to each road geological data according to formula (17);

Step 2: the flow process of utilization table 2 carries out iteration, obtains the coefficient that useful signal is corresponding;

Step 3: calculate analytic signal by formula (2);

Step 4: calculate instantaneous attribute according to formula (3)-(5).

Effect analysis

Composite signal example

The present invention's three test signals check distinct methods to calculate the effect of instantaneous frequency.Fig. 2 (a)-(c) is the Ricker wavelet of 50Hz respectively, the seismologic record synthesized by Ricker wavelet and the reflection coefficient sequence convolution of 50Hz and actual seismic track data.First, the traditional Hilbert method of the present invention calculates the instantaneous frequency of three signals, as shown in Fig. 2 (d)-(f), then, adopt the method that the present invention proposes, select β=1, the broad sense Morse Wavelet Frames of γ=3 calculates instantaneous frequency, as shown in Fig. 2 (g)-(i), the result that visible two kinds of methods obtain is almost identical.

Fig. 2 is the instantaneous frequency that test signal and distinct methods calculate, and as seen in not noisy situation, two kinds of methods all can obtain instantaneous frequency result of calculation accurately.

But when signal is by noise, the instantaneous frequency accuracy that Hilbert converter technique calculates reduces greatly.The present invention adds noise to three test signals, its signal to noise ratio (S/N ratio) is 5dB (Fig. 3 (a)-(c)), then again instantaneous attribute is calculated, as shown in Fig. 3 (d)-(f) He Fig. 3 (g)-(i) by two kinds of methods.

The instantaneous frequency that the noisy test signal of Fig. 3 and distinct methods calculate, visible in the result of calculation of Hilbert converter technique, the instantaneous frequency of useful signal is completely by noise takeover, and adopt the result of calculation of the inventive method, even if under very noisy condition, the instantaneous frequency of useful signal still clearly can be demonstrated

In the method as proposed in the present invention, by solving-optimizing problem, obtain the coefficient corresponding with useful signal, then calculate analytic signal and instantaneous frequency.Visible, even if when noise ratio is more serious, the method still accurately can obtain the instantaneous frequency of test signal.And the result that traditional Hilbert transform method obtains, the instantaneous frequency of useful signal, cannot identification completely by noise effect.

Below, together with the method propose the present invention and the instantaneous frequency of the noisy test signal calculated based on Hilbert transform method are drawn in the instantaneous frequency of not signals and associated noises, and select different Wavelet Frames to calculate in the method as proposed in the present invention, result as shown in Figure 4.

The instantaneous frequency of the noisy test signal that Fig. 4 distinct methods calculates, wherein blue line is the instantaneous frequency of noisy test signal, and red line is the instantaneous frequency of not signals and associated noises.The method adopting Hilbert method and the present invention to propose respectively calculates, and in the method as proposed in the present invention, selects broad sense Morse frame (GMW) and Morlet Wavelet Frames respectively.Can find out that the method that the present invention proposes has better noiseproof feature and precision, meanwhile, adopt the result of broad sense Morse frame to be better than Morlet frame

Fig. 4 (a)-(c) depicts the instantaneous frequency calculated with Hilbert converter technique, Fig. 4 (d)-(f) depicts the method proposed with the present invention, select the instantaneous frequency that broad sense Morse frame calculates, the method that Fig. 4 (g)-(i) proposes for adopting the present invention, selects the instantaneous frequency that Morlet frame calculates.

In order to the noiseproof feature of the instantaneous frequency that quantitative measurement distinct methods calculates, invention defines instantaneous frequency signal to noise ratio (S/N ratio):

$\mathrm{FSNR}=10{\log }_{10}\left(\frac{{\left|\right|f_0\left|\right|}^2}{{\left|\right|f_n-f_0\left|\right|}^2}\right)\left(\mathrm{dB}\right),---\left(75\right)$

Wherein f ₀the instantaneous frequency of not signals and associated noises, f _nfor the instantaneous frequency of signals and associated noises.

The instantaneous frequency signal to noise ratio (S/N ratio) FSNR of each method is listed in table 3.

The instantaneous frequency signal to noise ratio (S/N ratio) FSNR of table 3 distinct methods

Can find out, due to the impact of noise, the instantaneous frequency accuracy based on Hilbert transformation calculations reduces greatly, even cannot pick out the instantaneous frequency of useful signal.And the instantaneous frequency that the method that the present invention proposes calculates has good noiseproof feature and accuracy.Can find out, the frame formed after adopting the precision of broad sense Morse frame to be better than conventional Morlet wavelet discrete simultaneously.

Below, the method proposed in the present invention is asked in the process of instantaneous frequency, adopt iteration atrophy thresholding algorithm (IST) and iteratively faster atrophy thresholding algorithm (FIST) solving-optimizing problem respectively, after every single-step iteration, all calculate an instantaneous frequency, obtain instantaneous frequency signal to noise ratio (S/N ratio) FSNR, as shown in Figure 5.

Result shows, the speed of convergence of iteratively faster atrophy algorithm (FIST) is better than iteration atrophy thresholding algorithm (IST), the solution of the problem that can be optimized after less iterations.In actual applications, due to the requirement of computing velocity, too much iterations is unrealistic, usually needs to arrange a maximum iteration time.In this experiment, the present invention arranges maximum iteration time N=24, can see that the method can obtain convergence solution within 10 times.Due to the effect of dynamic stopping criterion, make the actual iterations used little, this can save computing time greatly, is convenient to the attributive analysis carrying out large data volume.

The instantaneous frequency signal to noise ratio (S/N ratio) of Fig. 5 tri-test signals, in figure, red line is for adopting iteration atrophy threshold value (IST) algorithm, blue line adopts iteratively faster threshold shrink (FIST) algorithm, can find out that the latter can obtain higher instantaneous frequency signal to noise ratio (S/N ratio) after less iterations

Actual seismic data example

The present invention have chosen a section of certain oil field poststack 3-D data volume, calculates instantaneous frequency respectively, as shown in Figure 6 by the method that Hilbert converter technique and the present invention propose.Can see that method of the present invention has good noiseproof feature and accuracy, the instantaneous frequency section obtained clearly can reflect the change of instantaneous frequency, instruction abnormal area.

Fig. 6 real data example.The method adopting Hilbert converter technique and the present invention to propose respectively calculates instantaneous frequency section, and the latter has good noiseproof feature, more clearly can reflect the change of instantaneous frequency, instruction abnormal area

Below, the method extracting instantaneous attribute based on broad sense Morse frame is used for the D seismic data processing of certain oil field Sandstone Gas Reservoir by the present invention, utilize business software respectively and extract based on broad sense Morse frame instantaneous frequency and the instant bandwidth that the method for instantaneous attribute calculates this three dimensional seismic data, then cut into slices along extraction rock stratum, destination layer position.The instantaneous frequency section utilizing business software to extract is as shown in Fig. 7 (a), and because the business software of routine calculates instantaneous frequency based on traditional Hilbert transform method, very responsive to noise, result is serious by noise.Based on broad sense Morse frame extract instantaneous frequency horizon slice as shown in Fig. 7 (b).Can be seen by contrast, the method that the present invention proposes has very high anti-noise ability, portrays more obvious to the space distribution of destination layer sand body.

Adopt the instant bandwidth horizon slice of routine business software and the inventive method calculating respectively as shown in Figure 8.Can see that the former is affected by noise comparatively large equally, bring a lot of inconvenience to seismic interpretation, and the latter be affected by noise less, clearly can portray geologic structure.

A instantaneous frequency horizon slice that () utilizes business software to calculate

B instantaneous frequency horizon slice that () calculates based on the inventive method

The instantaneous frequency horizon slice of Fig. 7 three dimensional seismic data

A instant bandwidth horizon slice that () utilizes business software to calculate

B instant bandwidth horizon slice that () calculates based on the inventive method

The instant bandwidth horizon slice of Fig. 8 three dimensional seismic data.

Claims (4)

1., based on the earthquake instantaneous attribute extracting method of broad sense Morse frame, it is characterized in that comprising the following steps:

Step 1: acquiring seismic data body;

Step 2: calculate the broad sense Morse frame coefficient that the noisy geological data in each road is corresponding;

The wavelet transformation of real seismic trace s (t) is defined as:

$W_{\mathrm{\&psi;}}(t,a)\&equiv;\frac{1}{a}{\&Integral;}_{-\&infin;}^{+\&infin;}s\left(\mathrm{\&tau;}\right){\mathrm{\&psi;}}^{*}\left(\frac{\mathrm{\&tau;}-t}{a}\right)\mathrm{d\&tau;},\left(1\right)$

In formula, wavelet ψ (t) selects broad sense Morse small echo;

Broad sense Morse wavelet frequency domain expression formula is:

${\mathrm{\&psi;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}\left(\mathrm{\&omega;}\right)=U\left(\mathrm{\&omega;}\right){\mathrm{\&alpha;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}{\mathrm{\&omega;}}^{\mathrm{\&beta;}}{e}^{-{\mathrm{\&omega;}}^{\mathrm{\&gamma;}}},\left(2\right)$

Wherein U (ω) is unit step function, β and γ is the parameter of small echo, and β > 0, γ > 0, α _{beta, gamma}for normalization constant, ${\mathrm{\&alpha;}}_{\mathrm{\&beta;},\mathrm{\&gamma;}}\&equiv;2{(\mathrm{e\&gamma;}/\mathrm{\&beta;})}^{\mathrm{\&beta;}/\mathrm{\&gamma;}};$

Little wave system after discretize can be expressed as:

${\mathrm{\&psi;}}_{m,n}\left(x\right)=a_0^{-m}\mathrm{\&psi;}\left[a_0^{-m}(x-{\mathrm{nt}}_0)\right]---\left(3\right)$

Wherein a ₀the discretize step-length of scale factor a and a ₀> 1, t ₀for the discretize step-length of shift factor t;

With the wavelet representation after the discretize that formula (1) is corresponding be:

C _m,n＝<s,ψ _m,n>, (4)

Wherein C _m,nfor frame coefficient;

Step 3: the broad sense Morse frame coefficient iteration for each road obtains coefficient corresponding to useful signal;

Definition K be from be mapped to operator, by frame coefficient be mapped as on signal, namely

The adjoint operator K of operator K ^*for from be mapped to operator, will on signal project on frame coefficient:

K ^*s＝<s,ψ′ _m,n>, (6)

Then s=KK ^*s, (7)

Operator K is composite operator, K ^*for analyze operator and Wavelet Frames ψ ' _m,ncorresponding;

Signals and associated noises is expressed as y=s+n=Kx+n, (8)

Wherein y represents signals and associated noises, and s is not noisy useful signal, and n is white Gaussian noise, and K is the composite operator in formula (22), and x represents the coefficient of transform domain, obtains coefficient corresponding to useful signal by solving following optimization problem

$\begin{array}{ccc}\tilde{x}=\arg \underset{x}{\min }{\left|\right|x\left|\right|}_1& s.t.& {\left|\right|y-\mathrm{Kx}\left|\right|}_2\&le;\mathrm{\&epsiv;}\end{array},,---\left(9\right)$

Wherein ε is relevant with the noise level of signal to be analyzed, with Lagrange multiplier λ, this problem is converted into following unconstrained problem:

${\tilde{x}}_{\mathrm{\&lambda;}}=\arg \underset{x}{\min }\frac{1}{2}{\left|\right|y-\mathrm{Kx}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x\left|\right|}_1,---\left(10\right)$

Wherein λ is called as regularization parameter, adopts the alternative manner shown in following formula to upgrade coefficient x

x ^(k+1)＝T _λ[x ^(k)+K ^*(y-Kx ^(k))],k＝1,2,…,N, (11)

Wherein T _λfor threshold function table, be defined as

$T_{\mathrm{\&lambda;}}\left(x\right)=\left\{\begin{array}{cc}x-\mathrm{\&lambda;}\frac{x}{\left|x\right|}& \left|x\right|\&GreaterEqual;\mathrm{\&lambda;},\\ 0& \left|x\right|<\mathrm{\&lambda;}.\end{array}\right.---\left(12\right);$

Step 4: calculate analytic signal by formula (2);

$c\left(t\right)\&equiv;\frac{1}{K_{\mathrm{\&psi;}}}{\&Integral;}_0^{\&infin;}{W}_{\mathrm{\&psi;}}(t,a)\frac{\mathrm{da}}{a}=s\left(t\right)+\mathrm{ih}\left(t\right),---\left(13\right)$

Wherein W _ψ(t, a) is the coefficient that useful signal is corresponding, and h (t) is the Hilbert conversion of s (t), and c (t) is the analytic signal that s (t) is corresponding;

Step 5: utilize c (t) to calculate instantaneous amplitude, instantaneous phase and instantaneous frequency:

$\hat{e}\left(t\right)=\sqrt{{\mathrm{Re}}^2\left[c\right]\left(t\right)+{\mathrm{Im}}^2\left[c\right]\left(t\right)},---\left(14\right)$

$\widehat{\mathrm{\&theta;}}\left(t\right)=\arctan \left[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\right],---\left(15\right)$

$\hat{f}\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d}{\mathrm{dt}}\left[\arctan \right[\frac{\mathrm{Im}\left[c\left(t\right)\right]}{\mathrm{Re}\left[c\left(t\right)\right]}\left]\right],---\left(16\right)$

Wherein Re [c (t)] and Im [c (t)] represents real part and the imaginary part of c (t) respectively.

2. the earthquake instantaneous attribute extracting method based on broad sense Morse frame according to claim 1, it is characterized in that: described step 3 formula (26) is each time in iterative process, index threshold is adopted to decline tactful, as shown in formula (28):

${\mathrm{\&lambda;}}_k={\left(\frac{{\mathrm{\&lambda;}}_{\min }}{{\mathrm{\&lambda;}}_{\max }}\right)}^{(k-1)/(N-1)}\&CenterDot;{\mathrm{\&lambda;}}_{\max },k=\mathrm{1,2},...,N.---\left(17\right)$

Wherein

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\max }=p_{\max }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\\ {\mathrm{\&lambda;}}_{\min }=p_{\min }\&CenterDot;\underset{i}{\max }\left\{\right|{\left(K^{*}y\right)}_i\left|\right\},\end{array}\right.---\left(18\right)$

λ in formula (29) _maxand λ _minbe respectively minimum value and the maximal value of regularization parameter, p _maxand p _minfor minimum and maximum number percent.

3. the earthquake instantaneous attribute extracting method based on broad sense Morse frame according to claim 1, is characterized in that: in order to stop iterative process in time in described step 3, defines a dynamic stopping criterion

$\frac{{\left|\right|x^{(k+1)}-x^{\left(k\right)}\left|\right|}^2}{{\left|\right|x^{(k+1)}\left|\right|}^2}<\mathrm{tolerance},---\left(19\right)$

Wherein tolerance is allowable value given before iteration, and under this criterion, when continuing iteration and can not obtaining better result, iterative process just stops automatically.

4. the earthquake instantaneous attribute extracting method based on broad sense Morse frame according to claim 1, is characterized in that: in order to improve speed of convergence in described step 3, adopts iteratively faster atrophy thresholding algorithm to see formula (31);

$\left\{\begin{array}{c}{\stackrel{!}{x}}^{\left(k\right)}=x^{\left(k\right)}+\frac{t^{\left(k\right)}-1}{t^{(k+1)}}(x^{\left(k\right)}-x^{(k-1)}),\\ x^{(k+1)}=T_{\mathrm{\&lambda;}}[{\stackrel{!}{x}}^{\left(k\right)}+K^{*}(y-K{\stackrel{!}{x}}^{\left(k\right)})],\end{array}\right.---\left(20\right)$

Wherein x ⁽⁰⁾=0, t ⁽⁰⁾=1, t ^(k)meet following recursion formula:

$t^{(k+1)}=(1+\sqrt{1+{4\left(t^{\left(k\right)}\right)}^2})/2---\left(21\right).$