CN106291700A - Based on the earthquake weighted average instantaneous frequency distilling method synchronizing extruding conversion - Google Patents

Abstract

The invention discloses a kind of based on the earthquake weighted average instantaneous frequency distilling method synchronizing extruding conversion, the present invention is on the basis of Three parameter wavelet converts, by carrying out energy rearrangement after estimating rearrangement criterion (instantaneous frequency), obtain more accurately, more sparse time-frequency representation, more accurately to determine useful signal Energy distribution space；Introduce threshold denoising on this basis, carry out noise compacting；Finally, the method for estimation of signals and associated noises weighting instantaneous frequency is given.The present invention is by calculating without making an uproar and the instantaneous frequency of noisy synthetic seismogram, the noiseproof feature of contrast distinct methods, the result of calculation of context of methods has good noiseproof feature and precision, and the method is applied to real data, it is possible to more clearly from portray the feature of reservoir.

Description

Based on the earthquake weighted average instantaneous frequency distilling method synchronizing extruding conversion

Technical field

The invention belongs to the signal processing field in geophysical exploration, be specifically related to a kind of based on synchronizing extruding conversion Earthquake weighted average instantaneous frequency distilling method.

Background technology

Seismic properties refers to the geometric shape of relevant seismic wave, the kinesiology spy derived by geological data through mathematic(al) manipulation Levy, dynamic characteristic and statistics feature, these features can be used to predict subsurface lithologic, describes the Reservoir Characters within oil reservoir Deng, process of seismic data processing has played the biggest effect.20 century 70s, seismic attributes analysis is directed initially into ground During shake is explained.Since the nineties in 20th century, the needs explained due to layer description and 3D data volume, Seismic attribute analysis technology Develop rapidly.But there is presently no a generally acknowledged seismic properties classification, according to physical significance, when seismic properties can be divided into Between, amplitude, frequency, several big classes such as relevant, decay, according to attribute pick-up method, when every class can be divided into instantaneous attribute, single track again Window attribute when window attribute, multiple tracks.Wherein instantaneous attribute is the attribute of pickup on the position that seismic wave arrives, and shakes including instantaneous Width, instantaneous phase, instantaneous frequency, instant bandwidth etc..In recent years, the quantity of seismic properties is uprushed, but, instantaneous attribute is still It is the pillar of geological data geologic interpretation, is also the technology modules of current commercial processes software indispensability.

In signal processing field, the concept of the signal transient attribute with instantaneous amplitude, instantaneous phase, instantaneous frequency as representative Long-standing, after the work of Gabor and Ville, a large amount of scholars are studied in association area, and Boashash is to these works Summarize.Taner et al. introduced complex seismic trace in 1979, proposed to be asked for the side of instantaneous attribute by complex seismic trace Method, and give the physical significance of these attributes and the application in seismic interpretation.Instantaneous amplitude becomes with the lithology of adjacent layer Change and oil-gas accumulation is relevant.The instantaneous phase reflection discontinuity at interface, tomography, plane of unconformity and sequence boundaries etc..Instantaneous frequency Rate, as one of the most frequently used instantaneous attribute, can effectively portray thickness and the variation of lithological on stratum, the distribution of Indication of Oil-Gas Deng.Zeng utilizes instantaneous frequency extremely to indicate thin layer.Gao et al. utilizes instantaneous frequency to carry out seismic data Q-value estimation.

The extraction of the present invention main Study of Seismic data instantaneous frequency.Instantaneous frequency attribute is a kind of time-varying parameter, definition The spectrum peak position converted over time for input non-stationary signal.According to its original definition, instantaneous frequency provides one To signal at the metric form of frequency domain energy accumulating.The conventional method extracting earthquake instantaneous frequency can be divided into following three Class:

(1) Hilbert alternative approach.Seismic channel doing Hilbert conversion, is converted into complex seismic trace, wherein real part is former Seismic channel data, imaginary part is its Hilbert conversion.After obtaining complex seismic trace, it is possible to calculate instantaneous frequency at each sampled point Rate.After the work of Taner, Hilbert converter technique is widely used in calculating earthquake instantaneous frequency, and major part business is soft so far Part still uses the method.The method is developed and has been improved by a lot of scholars, and Barnes proposed the multiple earthquake of two dimension in 1996 The concept of trace analysis.Luo et al. proposed generalized Hilbert transformation in 2003 and provides its application in terms of geophysics. Barnes proposed the concept of weighting instantaneous frequency in 2007.Hilbert alternative approach is promoted by Lu and Zhang, proposes The method calculating instantaneous frequency based on sef-adapting filter.

(2) method based on time frequency analysis.Time-Frequency Analysis Method utilizes the time-frequency distributions of signal to ask for instantaneous frequency. Boashash et al. proposes self adaptation instantaneous Frequency Estimation method based on time frequency analysis.Stankovic et al. utilizes self adaptation Window time-frequency distributions calculates instantaneous frequency.Gao Jinghuai et al. proposes the method calculating earthquake instantaneous frequency in phase space.Marfurt Et al. propose based on narrow-band spectrum analyze seismic data instantaneous frequency distilling method.Huang et al. proposes to divide based on empirical mode Solve the method that (EMD) calculates instantaneous frequency.Complete overall experience Mode Decomposition (CEEMD) method is used for extracting ground by Han et al. The instantaneous frequency of shake data.

(3) method based on inverting.Fomel et al. propose utilize the method for inverting to obtain local attribute, compared to wink Time frequency, local frequencies physical significance definitely, and application effect obvious.Liu et al. proposes method based on inverting, calculates The instantaneous frequency of seismic data.

Summary of the invention

It is an object of the invention to overcome above-mentioned deficiency, it is provided that a kind of based on the earthquake weighted average wink synchronizing extruding conversion Time frequency extraction method, there is good noiseproof feature and precision, the method be applied to real data, it is possible to more clearly from Portray the feature of reservoir.

In order to achieve the above object, the present invention comprises the following steps:

Step one, carries out Three parameter wavelet conversion according to real seismic trace；

Step 2, calculates and resets criterion, carries out energy by this criterion and resets the time-frequency representation more concentrated；

Step 3, resets energy, the time-frequency distributions after being extruded；

Step 4, carries out the determination in useful signal Energy distribution space, obtains the conversion coefficient that useful signal is corresponding；

Step 5, is weighted average instantaneous frequency distilling, completes average instantaneous frequency.

In described step one, the wavelet transformation of real seismic trace s (t) is defined as

$\begin{array}{c}{W}_s(a,b)=a^{-1}\∫s\left(t\right)\stackrel{\‾}{\mathrm{\ψ}}\left(\frac{t-b}{a}\right)dt\\ =\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}S\left(\mathrm{\ω}\right)\stackrel{\‾}{\widehat{\mathrm{\ψ}}}\left(a\mathrm{\ω}\right)e^{i\mathrm{\ω}b}d\mathrm{\ω},\end{array}---\left(1\right)$

Wherein,For the complex conjugate of ψ (t), a is scale factor, and b is time shift method,It it is wavelet ψ (t) Fourier converts,It is the Fourier conversion of s (t), wavelet square integrable, without DC component, meet and allow below Property condition:

$C_{\mathrm{\&psi;}}={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\frac{|\mathrm{\&psi;}\left(\mathrm{\&omega;}\right){|}^2}{\left|\mathrm{\&omega;}\right|}d\mathrm{\&omega;}<\mathrm{\&infin;}.---\left(2\right)$

The definition of Three parameter wavelet is:

$\mathrm{\&psi;}(t;\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&lsqb;cos\left(\mathrm{\&sigma;}t\right)-k(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&rsqb;+iq(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})sin\left(\mathrm{\&sigma;}t\right)\},---\left(3\right)$

Wherein, τ is the energy attenuation factor, and β is energy relay time, and σ is analysis wavelet modulating frequency, (σ, τ, β ∈ R and σ, τ ＞ 0), remember parameter sets, then ψ (t with vector Λ=(σ, τ, β)；σ, τ, β) ψ (t can be designated as；Λ), other amounts be similar to, with to The form of amount, above formula can be abbreviated as:

$\mathrm{\&psi;}(t;\mathrm{\&Lambda;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p\left(\mathrm{\&Lambda;}\right)\&lsqb;\cos \left(\mathrm{\&sigma;}t\right)-k\left(\mathrm{\&Lambda;}\right)\&rsqb;+iq\left(\mathrm{\&Lambda;}\right)\sin \left(\mathrm{\&sigma;}t\right)\},---\left(4\right)$

Wherein,

$\begin{array}{c}k\left(\mathrm{\&Lambda;}\right)=e^{-\frac{{\mathrm{\&sigma;}}^2}{4\mathrm{\&tau;}}}\&lsqb;\cos \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+i\frac{q\left(\mathrm{\&Lambda;}\right)}{p\left(\mathrm{\&Lambda;}\right)}\sin \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)\&rsqb;,\\ p\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\cos }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}},\\ q\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\sin }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}}.\end{array}---\left(5\right)$

(4) formula is made Fourier conversion, and obtaining its frequency domain form is

$\begin{array}{c}\widehat{\mathrm{\&psi;}}(\mathrm{\&omega;};\mathrm{\&Lambda;})={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\mathrm{\&psi;}(t;\mathrm{\&Lambda;})e^{-i\mathrm{\&omega;}t}dt\\ =\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)+q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}-\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}-\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}+\\ \sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)-q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}+\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}+\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}-\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}p\left(\mathrm{\&Lambda;}\right)k\left(\mathrm{\&Lambda;}\right)e^{-i\mathrm{\&beta;}\mathrm{\&omega;}-\frac{{\mathrm{\&omega;}}^2}{4\mathrm{\&tau;}}}.\end{array}---\left(6\right)$

In described step 2, it is assumed that a single-frequency cosine signal, i.e. s (t)=Acos (ω₀T), the Fourier of s (t) becomes It is changed toAccording to (1) formula, signal s (t) is about Three parameter wavelet ψ (t；Small echo Λ) Transformation results is

$\begin{array}{c}{W}_s(a,b)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}\\ =\frac{A}{2\mathrm{\&pi;}}\&lsqb;\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}+\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(-{\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{-{\mathrm{i\&omega;}}_{0}b}\&rsqb;.\end{array}---\left(7\right)$

Because Three parameter wavelet does not has negative frequency components, i.e. ω ＜ 0,And only consider positive yardstick, then (7) Formula can be reduced to

$W_s(a,b)=\frac{A}{2\mathrm{\&pi;}}\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}.---\left(8\right)$

From (8) if formula is it will be seen that Three parameter wavelet ψ (t；Crest frequency Λ) is ξ_M, then the result of wavelet transformation Will be at yardstick a=ξ_M/ω₀Place gets maximum, and forms a yardstick band centered by the yardstick of this energy maximum, causes The diffusion of energy, for the time-frequency distributions more concentrated, calculates and resets criterion instantaneous frequency:

${\mathrm{\&omega;}}_s(a,b)=\frac{{\&part;}_b{W}_s(a,b)}{2{\mathrm{\&pi;iW}}_s(a,b)},---\left(9\right)$

Wherein, W_s(a, b) ≠ 0, work as W_s(a, when b)=0, signal is in noenergy herein, and (9) formula, by infinity, calculates instantaneous Frequencies omega_s(a, b) nonsensical, so requirement when calculatingWhen noise is less, threshold value can be selectedIf noise is relatively big, needs to suppress noise, choose adaptive threshold as threshold valueWherein σ_η Being relevant with noise level, it is by the front n of wavelet conversion coefficient_vIndividual yardstick is estimated, Median is the built-in function in matlab, represents and takes intermediate value, and 0.6745 is for height The regularization factors of this noise, it is assumed that have N number of moment, it is possible to the σ that N number of moment is tried to achieve_ηTake average, optimum can be obtained Threshold value.

When reference signal is single-frequency cosine signal, its wavelet transformation can form a yardstick band, but can be counted by (9) Calculating its instantaneous frequency is

ω_s(a, b)=ω₀ (10)

From (10) formula it will be seen that the instantaneous frequency of (9) formula definition is exactly the frequency of single-frequency cosine signal, i.e. for one Cosine signal, its Three parameter wavelet conversion obtain time m-scale domain result can near the maximum yardstick of certain energy shape Become a yardstick band, but these yardstick correspondence wavelet coefficients are all cosine signal by the instantaneous frequency that (10) formula is calculated Frequencies omega₀, therefore imagine and the energy of these yardsticks all focused on ω₀On, thus the time-frequency representation more concentrated.

In described step 3, after (9) formula is calculated instantaneous frequency, learn that energy is concentrated, namely in this frequency Wavelet conversion coefficient is added in this frequency content, i.e. (a, b) → (ω_s(a b), b), carries out energy rearrangement.

For given signal s (t) ∈ L²(R) Three parameter wavelet is transformed to W_s(a, b), has a following expression:

$\begin{array}{c}{\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ ={\&Integral;}_0^{+\mathrm{\&infin;}}\frac{\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(\mathrm{\&omega;};\mathrm{\&Lambda;})}{\mathrm{\&omega;}}d\mathrm{\&omega;}\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}.\end{array}---\left(11\right)$

NoteAndIt is the analytic signal of s (t), thus (11) formula can be written as:

${\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=C_{\mathrm{\&psi;}}{s}_a\left(b\right).---\left(12\right)$

Because s_aB () is the analytic signal of s (b), therefore have

It will be seen that be multiplied by factor a-again by wavelet conversion coefficient from (12)¹The result obtained and the solution of original signal Analysis signal only differs from an invariant, and this is readily available its inverse transformation, as shown in (13) formula；If by W_s(a,b)a-¹All divide On dispensing corresponding instantaneous frequency composition mentioned above, then the time-frequency distributions after being extruded, it follows that time m-scale domain arrive The mapping of time-frequency domain is as follows:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_s(a,b)=\mathrm{\&omega;},W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}da.---\left(14\right)$

By (14) formula, at the moment b that certain is fixing, calculate its instantaneous frequency ω_s(a, b), by all instantaneous frequencys all Wavelet coefficient for a certain frequencies omega is accumulated together by (14), this completes the reassignment of energy, is extruded After time-frequency distributions, the equivalent form of value of (14) is:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}\mathrm{\&delta;}\left({\mathrm{\&omega;}}_s\right(a,b)-\mathrm{\&omega;})da.---\left(15\right)$

In described step 4, in order to more effectively determine the Energy distribution space of useful signal, it is proposed that after resetting Temporal frequency domain coefficient uses percentage threshold strategy to carry out noise compacting, i.e. soft-threshold operation, is defined as

$S\left(T_s\right)=\left\{\begin{array}{cc}{T}_s-\mathrm{\&delta;}\frac{T_s}{\left|T_s\right|},& \left|T_s\right|\&GreaterEqual;\mathrm{\&delta;}\\ 0,& \left|T_s\right|<\mathrm{\&delta;}\end{array}\right.,---\left(16\right)$

Wherein threshold parameter δ is calculated by following formula:

δ=p max{ | T_s|, (17) wherein p is percentage ratio, T_sFor the temporal frequency domain coefficient after conversion.

In described step 5, after determining the Energy distribution space that useful signal is corresponding, it is possible to utilize stability higher, more For sparse time-frequency domain conversation coefficient, calculate and weight instantaneous frequency:

$f\left(t\right)=\frac{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}f*T_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df}{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}{T}_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df},---\left(18\right)$

Wherein f is frequency, b and t is the time, and ξ (t) is Hamming window.

Compared with prior art, the present invention is on the basis of Three parameter wavelet converts, by estimating rearrangement criterion (instantaneous frequency Rate) after carry out energy rearrangement, obtain more accurately, more sparse time-frequency representation, more accurately to determine useful signal energy Amount distribution space；Introduce threshold denoising on this basis, carry out noise compacting；Finally, give signals and associated noises and weight instantaneous frequency The method of estimation of rate.The present invention, by calculating without making an uproar and the instantaneous frequency of noisy synthetic seismogram, contrasts the anti-of distinct methods Making an uproar performance, the result of calculation of context of methods has good noiseproof feature and precision, the method is applied to real data, it is possible to More clearly from portray the feature of reservoir.

Accompanying drawing explanation

Fig. 1 is the time-frequency domain spectrogram of the Ricker wavelet of 40Hz of the present invention；Wherein, the Ricker wavelet that (a) nothing is made an uproar, (b), C () is without making an uproar the wavelet transformation of Ricker wavelet and the present invention proposes the time-frequency figure of conversion；D () adds signal to noise ratio is the Gauss of 5dB Ricker wavelet after white noise, the wavelet transformation of (e), (f) noisy Ricker wavelet and the present invention propose the time-frequency figure of conversion；

Fig. 2 is that the present invention is without synthetic seismogram and the instantaneous frequency using distinct methods to calculate of making an uproar；Wherein, (a) The true record of synthesis, the instantaneous frequency that (b) calculates based on HT, (c), (d) are respectively and become based on Fourier conversion and synchronization extruding Change the weighted average instantaneous frequency of calculating；

Fig. 3 is the noisy synthetic seismogram of the present invention and the instantaneous frequency using distinct methods to calculate；Wherein, (a) Noisy synthetic seismogram, the instantaneous frequency that (b) calculates based on HT, (c), (d) are respectively based on Fourier conversion and synchronize The weighted average instantaneous frequency of extruding transformation calculations；

Fig. 4 is 2D real seismic record section of the present invention；Wherein, the position of arrow mark is that delta sandbody is stacked, because of Seismic data resolution is relatively low and influence of noise, can not reflect stacked feature；

Fig. 5 is the instantaneous frequency of the present invention the most calculated 2D actual seismic data；Wherein, (a) based on The instantaneous frequency that HT calculates, (b), (c), (d) are respectively based on Fourier conversion, wavelet transformation and synchronization extruding transformation calculations Weighted average instantaneous frequency.

Detailed description of the invention

The present invention will be further described with embodiment below in conjunction with the accompanying drawings.

(1) wavelet transformation definition

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

$\begin{array}{c}{W}_s(a,b)=a^{-1}\&Integral;s\left(t\right)\stackrel{\&OverBar;}{\mathrm{\&psi;}}\left(\frac{t-b}{a}\right)dt\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}S\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}\left(a\mathrm{\&omega;}\right)e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;},\end{array}---\left(1\right)$

Wherein,For the complex conjugate of ψ (t), a is scale factor, and b is time shift method,It it is wavelet ψ (t) Fourier converts,It it is the Fourier conversion of s (t).Wavelet square integrable, without DC component, meets and allows below Property condition:

$C_{\mathrm{\&psi;}}={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\frac{|\mathrm{\&psi;}\left(\mathrm{\&omega;}\right){|}^2}{\left|\mathrm{\&omega;}\right|}d\mathrm{\&omega;}<\mathrm{\&infin;}.---\left(2\right)$

Contain the feature of fast-changing amplitude and frequency component for thin interbed seismic signal, propose with Gao Jinghuai Based on the small echo of good coupling seismic wavelet (BMSW), traditional Morlet small echo is changed by the work using for reference Harrop et al. Entering, Gao Jinghuai et al. proposes a kind of wavelet function Three parameter wavelet (TPW).Owing to having three parameters, by regulation, It can adapt to our needs flexibly, is applied to it in sequence detection and the resolution of high precision seismic data all achieve very Significantly effect.The present invention chooses the Three parameter wavelet morther wavelet as wavelet transformation.The definition of Three parameter wavelet is:

$\mathrm{\&psi;}(t;\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&lsqb;cos\left(\mathrm{\&sigma;}t\right)-k(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&rsqb;+iq(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})sin\left(\mathrm{\&sigma;}t\right)\},---\left(3\right)$

Wherein, τ is the energy attenuation factor, and β is energy relay time, and σ is analysis wavelet modulating frequency, (σ, τ, β ∈ R and σ, τ ＞ 0).In order to write conveniently, remember parameter sets, then ψ (t with vector Λ=(σ, τ, β)；σ, τ, β) ψ (t can be designated as；Λ), its He measures similar.In vector form, above formula can be abbreviated as:

$\mathrm{\&psi;}(t;\mathrm{\&Lambda;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p\left(\mathrm{\&Lambda;}\right)\&lsqb;\cos \left(\mathrm{\&sigma;}t\right)-k\left(\mathrm{\&Lambda;}\right)\&rsqb;+iq\left(\mathrm{\&Lambda;}\right)\sin \left(\mathrm{\&sigma;}t\right)\},---\left(4\right)$

Wherein,

$\begin{array}{c}k\left(\mathrm{\&Lambda;}\right)=e^{-\frac{{\mathrm{\&sigma;}}^2}{4\mathrm{\&tau;}}}\&lsqb;\cos \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+i\frac{q\left(\mathrm{\&Lambda;}\right)}{p\left(\mathrm{\&Lambda;}\right)}\sin \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)\&rsqb;,\\ p\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\cos }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}},\\ q\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\sin }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}}.\end{array}---\left(5\right)$

(4) formula is made Fourier conversion, and obtaining its frequency domain form is

$\begin{array}{c}\widehat{\mathrm{\&psi;}}(\mathrm{\&omega;};\mathrm{\&Lambda;})={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\mathrm{\&psi;}(t;\mathrm{\&Lambda;})e^{-i\mathrm{\&omega;}t}dt\\ =\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)+q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}-\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}-\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}+\\ \sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)-q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}+\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}+\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}-\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}p\left(\mathrm{\&Lambda;}\right)k\left(\mathrm{\&Lambda;}\right)e^{-i\mathrm{\&beta;}\mathrm{\&omega;}-\frac{{\mathrm{\&omega;}}^2}{4\mathrm{\&tau;}}}.\end{array}---\left(6\right)$

(2) rearrangement criterion (instantaneous frequency) is calculated

Assume a single-frequency cosine signal, i.e. s (t)=Ac ω o₀The Fourier of s (, ts) (t) is transformed toAccording to (1) formula, signal s (t) is about Three parameter wavelet ψ (t；Wavelet transformation Λ) Result is

$\begin{array}{c}{W}_s(a,b)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}\\ =\frac{A}{2\mathrm{\&pi;}}\&lsqb;\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}+\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(-{\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{-{\mathrm{i\&omega;}}_{0}b}\&rsqb;.\end{array}---\left(7\right)$

Because the Three parameter wavelet that we select is almost without negative frequency components, i.e.And only examine Consider positive yardstick, then (7) formula can be reduced to

$W_s(a,b)=\frac{A}{2\mathrm{\&pi;}}\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}.---\left(8\right)$

From (8) if formula is it will be seen that Three parameter wavelet ψ (t；Crest frequency Λ) is ξ_M, then the result of wavelet transformation Will be at yardstick a=ξ_M/ω₀Place gets maximum, and forms a yardstick band centered by the yardstick of this energy maximum, causes The diffusion of energy, for the time-frequency distributions more concentrated, calculates and resets criterion instantaneous frequency:

${\mathrm{\&omega;}}_s(a,b)=\frac{{\&part;}_b{W}_s(a,b)}{2{\mathrm{\&pi;iW}}_s(a,b)},---\left(9\right)$

Wherein, W_s(a,b)≠0.Work as W_s(a, b) during 0=, signal is in noenergy herein, and (9) formula, by infinity, calculates instantaneous Frequencies omega_s(a, b) nonsensical, so requirement when calculatingWhen noise is less, threshold value can be selectedIf noise is relatively big, needing to suppress noise, we choose adaptive threshold as threshold valueIts Middle σ_ηRelevant with noise level, it is by the front n of wavelet conversion coefficient_vIndividual yardstick is estimated, Median is the built-in function in matlab, and representative takes Intermediate value, 0.6745 is the regularization factors for Gaussian noise.Assume there is N number of moment, it is possible to the σ that N number of moment is tried to achieve_η Take average, the threshold value of optimum can be obtained.

When reference signal is single-frequency cosine signal, its wavelet transformation can form a yardstick band, but can be counted by (9) Calculating its instantaneous frequency is

ω_s(a, b)=ω₀ (10)

From (10) formula it will be seen that the instantaneous frequency of (9) formula definition is exactly the frequency of single-frequency cosine signal, i.e. for one Cosine signal, its Three parameter wavelet conversion obtain time m-scale domain result can near the maximum yardstick of certain energy shape Become a yardstick band, but these yardstick correspondence wavelet coefficients are all cosine signal by the instantaneous frequency that (10) formula is calculated Frequencies omega₀, therefore we are it is envisioned that focus on ω by the energy of these yardsticks₀On, thus the frequency schedule more concentrated Show.

(3) energy is reset Time-Scale Domain energy and is mapped to temporal frequency domain

(9) after formula is calculated instantaneous frequency, it is understood that energy should be concentrated, namely by small echo in this frequency Conversion coefficient is added in this frequency content, i.e. (a, b) → (ω_s(a b), b), carries out energy rearrangement.

For given signal s (t) ∈ L²(R) Three parameter wavelet is transformed to W_s(a, b), has a following expression:

$\begin{array}{c}{\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ ={\&Integral;}_0^{+\mathrm{\&infin;}}\frac{\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(\mathrm{\&omega;};\mathrm{\&Lambda;})}{\mathrm{\&omega;}}d\mathrm{\&omega;}\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}.\end{array}---\left(11\right)$

NoteAndIt is the analytic signal of s (t), thus (11) formula can be written as:

${\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=C_{\mathrm{\&psi;}}{s}_a\left(b\right).---\left(12\right)$

Because s_aB () is the analytic signal of s (b), therefore have

It will be seen that be multiplied by factor a again by wavelet conversion coefficient from (12)^-1The result obtained and the parsing of original signal Signal only differs from an invariant, and this is readily available its inverse transformation, as shown in (13) formula.If we are by W_s(a,b)a^-1All Distribute on our corresponding instantaneous frequency composition mentioned above, then can be obtained by extruding after time-frequency distributions, thus we When drawing, m-scale domain is as follows to the mapping of time-frequency domain:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_s(a,b)=\mathrm{\&omega;},W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}da.---\left(14\right)$

By (14) formula, at the moment b that certain is fixing, calculate its instantaneous frequency ω_s(a, b), by all instantaneous frequencys all Wavelet coefficient for a certain frequencies omega is accumulated together by (14), this completes the reassignment of energy, is extruded After time-frequency distributions, the equivalent form of value of (14) is:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}\mathrm{\&delta;}\left({\mathrm{\&omega;}}_s\right(a,b)-\mathrm{\&omega;})da.---\left(15\right)$

(4) determination in useful signal Energy distribution space

Signals and associated noises is carried out wavelet transformation, by it when selecting suitable wavelet function (present invention chooses Three parameter wavelet) The when of projecting to Time-Scale Domain, the energy of useful signal will be distributed in less subspace V, and the energy of noise can spread To bigger subspace V ', m-scale domain time the most whole.In other words, useful signal is corresponding with the coefficient of minority, and Noise is almost distributed on whole coefficients.If we determined that the subspace V that useful signal is corresponding, obtain useful signal pair The coefficient answered, by other coefficient zero setting, then noise will be pressed at transform domain, and signal to noise ratio is improved.The present invention is by same After step extruding conversion, time scale domain coefficient, through resetting, transforms to temporal frequency domain, obtains the most sparse time-frequency knot Really, useful signal Energy distribution space diminishes, and is more beneficial for noise compacting.

In order to more effectively determine the Energy distribution space of useful signal, we have proposed the temporal frequency domain after resetting Coefficient uses percentage threshold strategy to carry out noise compacting, i.e. soft-threshold operation, is defined as

$S\left(T_s\right)=\left\{\begin{array}{cc}{T}_s-\mathrm{\&delta;}\frac{T_s}{\left|T_s\right|},& \left|T_s\right|\&GreaterEqual;\mathrm{\&delta;}\\ 0,& \left|T_s\right|<\mathrm{\&delta;}\end{array}\right.,---\left(16\right)$

Wherein threshold parameter δ is calculated by following formula:

δ=p max{ | T_s|, (17) wherein p is percentage ratio, T_sTemporal frequency domain coefficient after converting for the present invention.

By this filtering strategies, obtain the time-frequency domain coefficient that useful signal is corresponding, these coefficients constitute noise be pressed with After time-frequency domain spectrogram corresponding to signal.Fig. 1 (a) depicts the spectrogram of the Ricker wavelet of the most noisy 40Hz, Fig. 1 B () and (c) is the time-frequency spectrum after wavelet transformation and present invention conversion, it is seen that present invention conversion can more be concentrated, more Add sparse time-frequency representation.Fig. 1 (d) is the spectrogram of the Ricker wavelet (signal to noise ratio snr=5dB) of noisy 40Hz, it is seen that it By sound pollution, Fig. 1 (e) and (f) be use filtering strategies after the time-frequency spectrum that converts of the wavelet transformation that obtains and the present invention Figure, it can be seen that owing to present invention conversion is the most sparse, noise is substantially suppressed, and useful signal is clearly showed.

(5) weighted average instantaneous frequency distilling

After determining the Energy distribution space that useful signal is corresponding, it is possible to the stability utilizing the present invention to obtain is higher, more For sparse time-frequency domain conversation coefficient, calculate and weight instantaneous frequency:

$f\left(t\right)=\frac{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}f*T_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df}{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}{T}_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df},---\left(18\right)$

Wherein f is frequency, b and t is the time, and ξ (t) is Hamming window.

Embodiment:

(1) composite signal example

The synthetic seismogram that the Ricker wavelet of 30Hz is generated by we with reflection coefficient sequence convolution is studied.Figure 2 (a) is the most noisy composite traces, and (b), (c), (d) are respectively the instantaneous frequency of Hilbert transformation calculations, based on Fourier Conversion and the present invention calculated weighting instantaneous frequency.The instantaneous Frequency Estimation result big rise and fall of visible (b) and (c), no Enough stable；The estimated result of the present invention is more stable, instantaneous frequency is mainly distributed near Ricker wavelet dominant frequency 30Hz, In the range of 20Hz-40Hz.

Fig. 3 (a) is the noisy record (signal to noise ratio snr=5dB) plus white Gaussian noise, and (b), (c), (d) are respectively The instantaneous frequency of Hilbert transformation calculations, based on Fourier conversion and the present invention calculated weighting instantaneous frequency.From Fig. 3 (b) it can be seen that in the instantaneous frequency result that obtained by traditional Hilbert alternative approach, the instantaneous frequency of useful signal Because the existence of noise is covered completely.Based on Fourier average weighted instantaneous Frequency Estimation result, although have certain Noise immunity, but estimated result is still not sufficiently stable.The result (Fig. 3 (d)) calculated by the inventive method is remained able to clearly Ground, stably show the instantaneous frequency of useful signal.

(2) actual seismic data example

In order to check the effectiveness of the inventive method further, we use it for actual seismic data.Fig. 4 is three dimensions According to a section of body, this data random noise is extremely serious, has a strong impact on reservoir prediction.Disclose according to well data message, layer It is delta facies between H6 and H7 of position, is that delta sandbody is stacked at red arrow indication, but is limited by seismic resolution and noise System, it is impossible to identify delta sandbody border and be stacked.Fig. 5 (a-d) respectively be use Hilbert transformation calculations instantaneous frequency and The weighting instantaneous frequency calculated based on Fourier conversion, wavelet transformation, the inventive method.Fig. 5 (a) can be seen that delta Sand body is stacked, but causes sand body the most continuous owing to signal to noise ratio is relatively low.Owing to introducing window function, to a certain extent in Fig. 5 (b) Reduce influence of noise, but resolution is the most impacted, can not portray delta sandbody well and be stacked phenomenon simultaneously.Figure 5 (c) has suppressed noise to a certain extent, and accuracy of detection is relatively preferable, but the highest owing to being limited to time frequency resolution, noise Underpressing is obvious, and section signal to noise ratio is the most relatively low.Fig. 5 (d) result noise immunity is remarkably reinforced, and estimates the instantaneous frequency letter obtained Ratio of making an uproar is substantially improved, and the stacked phenomenon of delta sandbody becomes apparent from, it is seen that this method is the effective means of detection lithology.

Claims (6)

1. based on the earthquake weighted average instantaneous frequency distilling method synchronizing extruding conversion, it is characterised in that comprise the following steps:

Step one, carries out Three parameter wavelet conversion according to real seismic trace；

Step 2, calculates and resets criterion, utilize this criterion can carry out resetting the time-frequency representation more concentrated；

Step 3, resets energy, the time-frequency distributions after being extruded；

Step 4, carries out the determination in useful signal Energy distribution space, obtains useful signal correspondence distribution space (i.e. transformation series Number)；

Step 5, is weighted average instantaneous frequency distilling, completes average instantaneous frequency.

The most according to claim 1 based on synchronizing the earthquake weighted average instantaneous frequency distilling method that extruding converts, it is special Levying and be, in described step one, the wavelet transformation of real seismic trace s (t) is defined as

$\begin{array}{c}{W}_s(a,b)=a^{-1}\&Integral;s\left(t\right)\stackrel{\&OverBar;}{\mathrm{\&psi;}}\left(\frac{t-b}{a}\right)dt\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}S\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}\left(a\mathrm{\&omega;}\right)e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;},\end{array}---\left(1\right)$

Wherein,For the complex conjugate of ψ (t), a is scale factor, and b is time shift method,It it is wavelet ψ (t) Fourier converts,It is the Fourier conversion of s (t), wavelet square integrable, without DC component, meet and allow below Property condition:

$C_{\mathrm{\&psi;}}={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\frac{|\mathrm{\&psi;}\left(\mathrm{\&omega;}\right){|}^2}{\left|\mathrm{\&omega;}\right|}d\mathrm{\&omega;}<\mathrm{\&infin;}.---\left(2\right)$

The definition of Three parameter wavelet is:

$\mathrm{\&psi;}(t;\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&lsqb;cos\left(\mathrm{\&sigma;}t\right)-k(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})\&rsqb;+iq(\mathrm{\&sigma;},\mathrm{\&tau;},\mathrm{\&beta;})sin\left(\mathrm{\&sigma;}t\right)\},---\left(3\right)$

Wherein, τ is the energy attenuation factor, and β is energy relay time, and σ is analysis wavelet modulating frequency, (σ, τ, β ∈ R and σ, τ ＞ 0), parameter sets, then ψ (t are remembered with vector Λ=(σ, τ, β)；σ, τ, β) ψ (t can be designated as；Λ), other amounts are similar to, with vector Form, above formula can be abbreviated as:

$\mathrm{\&psi;}(t;\mathrm{\&Lambda;})=e^{-\mathrm{\&tau;}{(t-\mathrm{\&beta;})}^2}\{p\left(\mathrm{\&Lambda;}\right)\&lsqb;\cos \left(\mathrm{\&sigma;}t\right)-k\left(\mathrm{\&Lambda;}\right)\&rsqb;+iq\left(\mathrm{\&Lambda;}\right)\sin \left(\mathrm{\&sigma;}t\right)\},---\left(4\right)$

Wherein,

$\begin{array}{c}k\left(\mathrm{\&Lambda;}\right)=e^{-\frac{{\mathrm{\&sigma;}}^2}{4\mathrm{\&tau;}}}\&lsqb;\cos \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+i\frac{q\left(\mathrm{\&Lambda;}\right)}{p\left(\mathrm{\&Lambda;}\right)}\sin \left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)\&rsqb;,\\ p\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\cos }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}},\\ q\left(\mathrm{\&Lambda;}\right)={\left(\frac{2\mathrm{\&tau;}}{\mathrm{\&pi;}}\right)}^{\frac{1}{4}}{\&lsqb;4(e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}-e^{-\frac{3{\mathrm{\&sigma;}}^2}{8\mathrm{\&tau;}}})\&times;{\sin }^2\left(\mathrm{\&beta;}\mathrm{\&sigma;}\right)+1-e^{-\frac{{\mathrm{\&sigma;}}^2}{2\mathrm{\&tau;}}}\&rsqb;}^{-\frac{1}{2}}.\end{array}---\left(5\right)$

(4) formula is made Fourier conversion, and obtaining its frequency domain form is

$\begin{array}{c}\widehat{\mathrm{\&psi;}}(\mathrm{\&omega;};\mathrm{\&Lambda;})={\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}\mathrm{\&psi;}(t;\mathrm{\&Lambda;})e^{-i\mathrm{\&omega;}t}dt\\ =\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)+q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}-\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}-\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}+\\ \sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}\frac{p\left(\mathrm{\&Lambda;}\right)-q\left(\mathrm{\&Lambda;}\right)}{2}{e}^{-i\mathrm{\&beta;}(\mathrm{\&omega;}+\mathrm{\&sigma;})-\frac{{(\mathrm{\&omega;}+\mathrm{\&sigma;})}^2}{4\mathrm{\&tau;}}}-\sqrt{\frac{\mathrm{\&pi;}}{\mathrm{\&tau;}}}p\left(\mathrm{\&Lambda;}\right)k\left(\mathrm{\&Lambda;}\right)e^{-i\mathrm{\&beta;}\mathrm{\&omega;}-\frac{{\mathrm{\&omega;}}^2}{4\mathrm{\&tau;}}}.\end{array}---\left(6\right)$

The most according to claim 1 based on synchronizing the earthquake weighted average instantaneous frequency distilling method that extruding converts, it is special Levy and be, in described step 2, it is assumed that a single-frequency cosine signal, i.e. s (t)=Acos (ω₀T), the Fourier conversion of s (t) ForAccording to (1) formula, signal s (t) is about Three parameter wavelet ψ (t；Small echo Λ) becomes Changing result is

$\begin{array}{c}{W}_s(a,b)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}\\ =\frac{A}{2\mathrm{\&pi;}}\&lsqb;\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}+\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(-{\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{-{\mathrm{i\&omega;}}_{0}b}\&rsqb;.\end{array}---\left(7\right)$

Because Three parameter wavelet does not has negative frequency components, i.e. ω ＜ 0,And only considering positive yardstick, then (7) formula is permissible It is reduced to

$W_s(a,b)=\frac{A}{2\mathrm{\&pi;}}\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}({\mathrm{a\&omega;}}_0;\mathrm{\&Lambda;})e^{{\mathrm{i\&omega;}}_{0}b}.---\left(8\right)$

From (8) if formula is it will be seen that Three parameter wavelet ψ (t；Crest frequency Λ) is ξ_M, then the result of wavelet transformation will be at chi Degree a=ξ_M/ω₀Place gets maximum, and forms a yardstick band centered by the yardstick of this energy maximum, causes energy Diffusion, for the time-frequency distributions more concentrated, calculates and resets criterion instantaneous frequency:

${\mathrm{\&omega;}}_s(a,b)=\frac{{\&part;}_b{W}_s(a,b)}{2{\mathrm{\&pi;iW}}_s(a,b)},---\left(9\right)$

Wherein, W_s(a, b) ≠ 0, work as W_s(a, when b)=0, signal is in noenergy herein, and (9) formula, by infinity, calculates instantaneous frequency ω_s(a, b) nonsensical, so requirement when calculatingWhen noise is less, threshold value can be selectedIt is 10^-8；As Really noise is relatively big, needs to suppress noise, chooses adaptive threshold as threshold valueWherein σ_ηIt is and noise water Flat relevant, it is by the front n of wavelet conversion coefficient_vIndividual yardstick is estimated, Median is the built-in function in matlab, represents and takes intermediate value, and 0.6745 is for height The regularization factors of this noise, it is assumed that have N number of moment, it is possible to the σ that N number of moment is tried to achieve_ηTake average, optimum can be obtained Threshold value；

When reference signal is single-frequency cosine signal, its wavelet transformation can form a yardstick band, but can be calculated it by (9) Instantaneous frequency is

From (10) formula it will be seen that the instantaneous frequency of (9) formula definition is exactly the frequency of single-frequency cosine signal, i.e. for a cosine Signal, the time m-scale domain result that its Three parameter wavelet conversion obtains can be formed about one at the yardstick that certain energy is maximum Individual yardstick band, but these yardstick correspondence wavelet coefficients are by frequency that the instantaneous frequency that (10) formula is calculated is all cosine signal Rate ω₀, therefore imagine and the energy of these yardsticks all focused on ω₀On, thus the time-frequency representation more concentrated.

The most according to claim 1 based on synchronizing the earthquake weighted average instantaneous frequency distilling method that extruding converts, it is special Levy and be, in described step 3, after (9) formula is calculated instantaneous frequency, learns that energy is concentrated in this frequency, namely will Wavelet conversion coefficient is added in this frequency content, i.e. (a, b) → (ω_s(a b), b), carries out energy rearrangement；

For given signal s (t) ∈ L²(R) Three parameter wavelet is transformed to W_s(a, b), has a following expression:

$\begin{array}{c}{\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ =\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(a\mathrm{\&omega;};\mathrm{\&Lambda;})e^{i\mathrm{\&omega;}b}{a}^{-1}dad\mathrm{\&omega;}\\ ={\&Integral;}_0^{+\mathrm{\&infin;}}\frac{\stackrel{\&OverBar;}{\widehat{\mathrm{\&psi;}}}(\mathrm{\&omega;};\mathrm{\&Lambda;})}{\mathrm{\&omega;}}d\mathrm{\&omega;}\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_0^{+\mathrm{\&infin;}}\hat{s}\left(\mathrm{\&omega;}\right)e^{i\mathrm{\&omega;}b}d\mathrm{\&omega;}.\end{array}---\left(11\right)$

NoteAndIt is the analytic signal of s (t), thus (11) formula Can be written as:

${\&Integral;}_0^{+\mathrm{\&infin;}}{W}_s(a,b)a^{-1}da=C_{\mathrm{\&psi;}}{s}_a\left(b\right).---\left(12\right)$

Because s_aB () is the analytic signal of s (b), therefore have

It will be seen that be multiplied by factor a again by wavelet conversion coefficient from (12)^-1The result obtained and the analytic signal of original signal Only differing from an invariant, this is readily available its inverse transformation, as shown in (13) formula；If by W_s(a,b)a^-1It is distributed on Face is mentioned on corresponding instantaneous frequency composition, then the time-frequency distributions after being extruded, it follows that time m-scale domain the most m- The mapping of frequency domain is as follows:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,{\mathrm{\&omega;}}_s(a,b)=\mathrm{\&omega;},W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}da.---\left(14\right)$

By (14) formula, at the moment b that certain is fixing, calculate its instantaneous frequency ω_s(all instantaneous frequencys b), are all a certain by a The wavelet coefficient of frequencies omega is accumulated together by (14), this completes the reassignment of energy, after being extruded time Frequency division cloth, the equivalent form of value of (14) is:

$T_s(\mathrm{\&omega;},b)={\&Integral;}_{\{a:a>0,W_s(a,b)\&NotEqual;0\}}{W}_s(a,b)a^{-1}\mathrm{\&delta;}\left({\mathrm{\&omega;}}_s\right(a,b)-\mathrm{\&omega;})da.---\left(15\right)$

The most according to claim 1 based on synchronizing the earthquake weighted average instantaneous frequency distilling method that extruding converts, it is special Levy and be, in described step 4, in order to more effectively determine the Energy distribution space of useful signal, it is proposed that to reset after time Between frequency domain coefficient use percentage threshold strategy to carry out noise compacting, i.e. soft-threshold operation, be defined as

$S\left(T_s\right)=\left\{\begin{array}{cc}{T}_s-\mathrm{\&delta;}\frac{T_s}{\left|T_s\right|},& \left|T_s\right|\&GreaterEqual;\mathrm{\&delta;}\\ 0,& \left|T_s\right|<\mathrm{\&delta;}\end{array}\right.,---\left(16\right)$

Wherein threshold parameter δ is calculated by following formula:

δ=p max{ | T_s|, (17)

Wherein p is percentage ratio, T_sFor the temporal frequency domain coefficient after conversion.

The most according to claim 1 based on synchronizing the earthquake weighted average instantaneous frequency distilling method that extruding converts, it is special Levy and be, in described step 5, after determining the Energy distribution space that useful signal is corresponding, it is possible to utilize stability higher, more For sparse time-frequency domain conversation coefficient, calculate and weight instantaneous frequency:

$f\left(t\right)=\frac{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}f*T_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df}{{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}{T}_s(2\mathrm{\&pi;}f,b)\mathrm{\&xi;}(t-b)df},---\left(18\right)$

Wherein f is frequency, b and t is the time, and ξ (t) is Hamming window.