CN103995289A - Time-varying mixed-phase seismic wavelet extraction method based on time-frequency spectrum simulation - Google Patents

Abstract

The invention discloses a time-varying mixed-phase seismic wavelet extraction method based on time-frequency spectrum simulation. The method is characterized in that on the basis that improved generalized S transformation is carried out on a non-stationary seismic record, a time-frequency filter is introduced for the first time to filter a time-frequency spectrum of the seismic record, then a wavelet amplitude spectrum at each moment is estimated according to the spectrum simulation method, a mixed-phase spectrum of the wavelets is reconstructed in the seismic record with the time-varying wavelet amplitude spectrums removed according to the bispectrum method of high-order cumulants under the assumed condition that the phases of the wavelets are time-invariant, and the amplitude spectrums and the phase spectrum are combined to achieve extraction of the time-varying mixed-phase wavelets. Compared with a conventional stage extraction method, the time-varying wavelet extraction method has the advantages that the assumption that stages of the wavelets are stable is not needed, the time-varying wavelets can be accurately estimated even though the attributes of the wavelets at adjacent stages vary greatly, the mixed phase of the extracted wavelets is more consistent with the reality and seismic data can be analyzed and processed more finely and more accurately. The method has great application value in improving the resolution of the seismic data.

Description

Based on time-frequency spectrum simulation time become method of mixed phase wavelet extraction

Technical field:

The invention belongs to seismic data processing field.

Background technology:

Seismic wavelet estimation is one of key issue in seismic data processing, is the basis of the seismic processing techniques such as Seismic forward, seismic inversion and deconvolution.The increase of oil-gas exploration difficulty makes the requirement of earthquake data process and interpretation harsher.In order to obtain high-resolution seismic data, the non-stationary composition in seismologic record can not be ignored.This is because under the absorption and filter action on stratum, wavelet its amplitude and energy in underground medium communication process are absorbed gradually, radio-frequency component is decayed gradually, therefore wavelet to be estimated has dynamic attenuation, and in the time of in recent years, the extraction of varitron ripple becomes oil-gas seismic exploration field high resolution processing technique important subject urgently to be resolved hurrily.

In fields such as signal analysis and processings, use mainspring and the starting point of time frequency analysis to be summed up as: time frequency analysis is the powerful of research non-stationary signal, its basic thought is the Copula of design time and frequency, signal is described simultaneously in energy density or the intensity of different time and frequency, the characteristic that represents signal with the form of joint time-frequency distribution, can locate comparatively exactly a certain moment and occur which frequency component.Consider when seismic wavelet is propagated in underground medium because its energy attenuation of effect, the frequency bands such as wavefront diffusion and attenuation by absorption narrow down, seismic wavelet is time dependent in communication process, the seismic data of record is a non-stationary signal, this seismologic record is carried out to time frequency analysis and can effectively portray its non-stationary property.

S conversion is in conjunction with Short Time Fourier Transform and wavelet transformation, overcome Short Time Fourier Transform time frequency resolution constant and the scale factor of wavelet transformation and the shortcoming of frequency-independent when inheriting their localization advantages.Yet the variation tendency of S conversion window function is changeless, can not adjust flexibly according to actual needs, and a lot of scholars are improved and proposed multiple generalized S-transform the window function of S conversion for this reason.When conventional generalized S-transform adopts, window width is with frequency f the Gauss function that inverse proportion changes, and when high band, window is narrower, can obtain higher temporal resolution, and when low frequency place, window is wider, can obtain higher frequency resolution.These generalized S-transforms can obtain yupin effect when good for most of non-stationary signals, but still be not well positioned to meet non-stationary seismologic record, do not process desired resolution.The window function of improved generalized S-transform has been introduced two regulatory factors, and the proportional example of its width and frequency changes, and at low frequency place, obtains higher temporal resolution, at high frequency treatment, obtains higher frequency resolution, and this character meets seismologic record dynamic attenuation characteristic.Seismologic record is carried out to improved generalized S-transform and there is good time-frequency focusing, can obtain good Time-Frequency Information, at compacting noise, the judgement position in reflection horizon and the aspects such as direct demonstration of oil gas, can obtain good application.

Consider that the time-frequency spectrum that time-frequency conversion obtains exists dispersion phenomenon, between adjacent reflection horizon, will certainly produce interference, this will affect the accuracy of latter earthquake data processing.In order better to adapt to the seismic data of thin interbed, process and explain, eliminating the impact between adjacent reflection horizon, introducing a kind of Time frequency Filter.This Time frequency Filter is the smooth function of time and frequency, and spatial point is had to time-frequency focusing.The time-frequency spectrum of seismologic record is carried out after filtering at time-frequency domain, can carry out varitron ripple when meticulous and estimate and seismic data processing.

Spectrum analog is generally used in deconvolution processing, is a kind of relatively more effective and stronger seismic wavelet analogy method of robustness.The spectral amplitude of its hypothesis wavelet is similar to the unimodal smooth curve of Ricker wavelet spectrum, and it is carried out to mathematical modeling in frequency field, adopts the method for least square fitting from seismologic record frequency spectrum, to obtain the frequency spectrum of wavelet.Compare with conventional autocorrelation estimation wavelet amplitude method, it is the restriction of white noise that spectrum analog method has been abandoned reflection coefficient, and can carry out at time-frequency domain.At time-frequency domain, carry out spectrum analog and fully retained the non-stationary composition in seismologic record, reacted truly the time variation of wavelet, can extract more subtly the wavelet of corresponding interval.

The feature that Higher Order Cumulants is applicable to seismic wavelet estimation problem can be summed up as: the complete information that (1) Higher Order Cumulants comprises signal; (2) because two are added up the cumulative amount sum that the semi-invariant of stochastic process independently equals two stochastic processes, so the Higher Order Cumulants of seismic convolution model can be decomposed into the semi-invariant of wavelet and reflection coefficient convolution and the form that noise storage amount is added; (3) for Gauss's coloured noise, Higher Order Cumulants perseverance is zero, can avoid gaussian additive noise on the impact of estimating.Therefore, Higher Order Cumulants has retained the abundanter information (especially phase place) of signal, can suppress Gaussian noise, can in coloured Gaussian noise, extract non-Gaussian signal, aspect the analysis of non-linear, non-gaussian sum non-minimum phase system and processing, demonstrate unique advantage.Consider that real seismic record neutron deficiency is mixed-phase, the wavelet extraction method based on Higher Order Cumulants is widely used in seismic data processing.

Summary of the invention:

The object of the invention is to propose a kind of based on time-frequency spectrum simulation time become method of mixed phase wavelet extraction.

The invention is characterized in, fully retained the non-stationary composition in seismologic record, and do not need reflection coefficient sequence to carry out white noise hypothesis, solved conventional stage extraction method in the situation that adjacent layer cross-talk ripple attribute exists larger variation, wavelet is estimated inaccurate problem, realizes and becomes seismic wavelet estimation when meticulousr; Adopt the mixed-phase of the two spectrometries extraction of high-order wavelet, widened the hypothesis of traditional wavelet extraction method neutron deficiency zero phase or minimum phase, the sub-wave phase attribute of extraction is more realistic; Aspect raising seismic data resolution, the method does not need the dynamic attenuation mechanism of Study of Seismic record, directly seismic data is processed, and has avoided error accumulation; The introducing of Time frequency Filter, has improved temporal resolution greatly, while having eliminated time frequency analysis, due to the impact existing between the adjacent reflection horizon of dispersion phenomenon, carries out more accurately on this basis follow-up work for the treatment of.The method comprises following steps successively:

Step (1) is carried out improved generalized S-transform time frequency analysis to non-stationary seismologic record, obtains time-frequency spectrum NGST (τ, f);

Step (2) is introduced Time frequency Filter the time-frequency spectrum of seismologic record is carried out to filtering processing, obtains the time-frequency spectrum S (τ, f) after processing;

The value of step (3) set time τ obtains putting sometime the frequency spectrum S (T of the seismologic record at place, f), the method of employing spectrum analog simulates the wavelet amplitude W (T in this moment, f), the position occurring according to reflection horizon changes the value of time τ, obtains the wavelet amplitude at each reflection spot place;

Step (4) becomes the impact that becomes wavelet amplitude when Wiener Filter Method is removed from seismologic record time-frequency spectrum while utilizing, residual sub-wave phase in the seismologic record after processing;

Step (5) supposes that sub-wave phase is time-independent, utilizes the mixed-phase spectrum φ (f) of two spectrometry reconstruct wavelets of Higher Order Cumulants;

Step (6) by time varitron wave amplitude spectrum and phase spectrum combination, obtain becoming when complete mixed phase wavelet.

Experiment simulation explanation:

The wavelet of the present invention by theogram estimate experiment check based on time-frequency domain spectrum analog time become the effect of mixed phase wavelet extracting method.First adopt the sparse model with five reflection coefficient pulses to test, its length is 1000ms, and sample frequency is 1ms, as shown in Fig. 2 (a).Position and the size of each reflection coefficient pulse is respectively (200ms, 1.0), (240ms, 0.8), (520ms,-1.0), (570ms, 0.9) and (840ms ,-0.8), comprising two adjacent reflection coefficients in the same way and two adjacent incorgruous reflection coefficients.Adopt dominant frequency to be respectively 60HZ, 50HZ, 40HZ, 35HZ, 30HZ, and amplitude energy the zero phase Ricker wavelet and the synthetic seismologic record of decaying of reflection coefficient sequence convolution of decaying gradually, as shown in Fig. 2 (b).

Synthetic seismologic record is carried out respectively to S conversion (Fig. 3 (a)), generalized S-transform (Fig. 3 (b)) and improved generalized S-transform (Fig. 3 (c)).Contrast can find out, improved generalized S-transform has been inherited the advantage of S conversion and generalized S-transform, has good time-frequency focusing.Seismologic record is carried out to improved generalized S-transform, can obviously find out the decay of amplitude energy, and can reflect exactly frequency content and the difference thereof on adjacent two stratum, under the prerequisite not reducing in retention time resolution, improved the precision of frequency resolution.Yet, the time point that former reflection coefficient occurs is 200ms, 240ms, 520ms, 570ms, 840ms, from improved generalized S-transform spectrum, although the energy at these reflection spot places is maximum, but in neighborhood, have frequency dispersion, this has a certain impact to subsequent estimation wavelet and then removal wavelet composition, and Fig. 3 (d) is to the result after improved generalized S-transform spectral filter, dispersion phenomenon weakens greatly, can obtain better temporal resolution.

Obtain after the time-frequency spectrum after seismologic record is processed, we extract respectively each seismic wavelet constantly.Original wavelet amplitude when Fig. 4 is respectively 200ms, 520ms, 570ms, 650ms, the spectral amplitude of seismologic record and the wavelet amplitude comparison diagram simulating.Fig. 5 is respectively original wavelet and the abdominal muscle wavelet time-frequency comparison diagram at 200ms, 520ms, 570ms place.From Fig. 4 and Fig. 5, can find out: the wavelet amplitude that (1) simulates and original wavelet amplitude are substantially identical, inverse transformation is also so to time domain, thereby proof spectrum analog extracts the accuracy of wavelet amplitude; (2) contrast 4 (b) and Fig. 4 (c), although exist and influence each other between adjacent stratum, method in this paper still can accurately extract the wavelet amplitude of every one deck; (3), along with the increase of time, its dominant frequency of the wavelet simulating reduces gradually and energy reduces, and can reflect the dynamic attenuation characteristic of wavelet in communication process; (4) extract the wavelet at reflection spot time point place in addition, although can obtain corresponding wavelet amplitude, as shown in Fig. 4 (d), the more original wavelet of its frequency has certain error, and energy is lower.

After obtaining the spectral amplitude of zero-phase wavelet, inverse transformation, to time domain, has just completed the extraction of wavelet.When seismologic record is carried out with the time varitron ripple obtaining, become Wiener filtering, dynamic deconvolution is processed, and result as shown in Figure 6.The time varitron ripple that used time frequency domain spectra simulation is extracted carries out the seismologic record after deconvolution processing, has slackened the impact of wavelet.The treatment effect of (a) rear (b) before filtering in comparison diagram 6, clearly filtering is processed deconvolution result afterwards and is more approached pulse at reflection spot place.If adopt the method for staging treating, extract wavelet, 200ms and 240ms and 520ms and 570ms were probably assigned in the period, within each period, extract like this wavelet under average meaning, the actual wavelet at place, this wavelet stratum adjacent with each all has certain error, and it must be inaccurate by this wavelet, carrying out deconvolution processing.

In order to verify the necessity of accurate extraction wavelet mixed-phase, select a mixed phase wavelet (Fig. 7 (a)) and above-mentioned reflection coefficient sequence to synthesize non-stationary seismologic record, as shown in Fig. 7 (b).

While first adopting the method for time-frequency spectrum in this paper simulation to extract after varitron wave amplitude, by time become Wiener Filter Method from removing (result is as shown in Fig. 8 (b)) seismologic record, the seismologic record of the two spectrometries that then adopt Higher Order Cumulants varitron wave-amplitude when removing, extract sub-wave phase, thereby while having completed, become the extraction of mixed phase wavelet, and further from removing the mixed-phase (as shown in Fig. 8 (c)) of wavelet seismologic record, obtain reflection coefficient sequence.From Fig. 8 (b), can find out, remove the seismologic record of wavelet amplitude, the impact of its wavelet obviously reduces, and frequency band broadens, and so, due to the existence of sub-wave phase, still has certain gap with reflected impulse, can not well judge the polarity at reflection spot place.Fig. 8 (c) further removes the result of wavelet mixed-phase on the basis of Fig. 8 (b), clearly at reflection spot place, more approach pulse, and polarity is clear and legible.If this explanation phase estimation is inaccurate, residual phase place can have a strong impact on seismic data treatment effect, becomes wavelet amplitude and phase spectrum is the basis of latter earthquake data processing and explanation while therefore extracting exactly.

In order to further illustrate problem, finally synthetic one relatively approaches actual non-stationary seismologic record.The random sparse reflectivity (Fig. 9 (a)) that generates, produces a non-stationary seismologic record (Fig. 9 (b)) with the mixed phase wavelet convolution of (2) part.By method in this paper, this seismologic record is processed, be directly used in deconvolution processing while obtaining after varitron ripple, result is as shown in Fig. 9 (c).Before and after processing, the contrast of geological data, can it is evident that, the seismic data of processing by method in this paper, more approaches pulse train at reflection spot place, thereby varitron ripple when explanation this method can accurately be extracted effectively improves seismologic record resolution.

In sum, the estimated accuracy of varitron ripple when the method that the present invention proposes can effectively improve, varitron ripple in the time of particularly also accurately estimating when adjacent layer cross-talk ripple attribute exists larger variation, for follow-up meticulousr, analyzing and processing seismic data provides assurance exactly.

Accompanying drawing explanation:

The process flow diagram of Fig. 1, accompanying method of the present invention

Fig. 2, reflection coefficient sequence (a) and non-stationary theogram (b)

Improved generalized S-transform (b) result after Fig. 3, S conversion (a), generalized S-transform (b), improved generalized S-transform (c) and filtering are processed

The spectral contrast figure that Fig. 4,200ms (a), 520ms (b), 570ms (c), 650ms (d) locate

The time domain wavelet comparison diagram that Fig. 5,200ms (a), 520ms (b), 570ms (c) locate

Before Fig. 6, filtering, after (a), the seismologic record after removing is composed in (b) pointwise

Fig. 7, mixed phase wavelet (a) and non-stationary theogram (b)

Fig. 8, theogram (a), become the result (c) of wavelet amplitude and phase spectrum when becoming the result (b) after wavelet amplitude and removing while removing

Result (c) after Fig. 9, reflection coefficient sequence (a), theogram (b) and deconvolution

Embodiment:

The present invention proposes based on time-frequency spectrum simulation time become method of mixed phase wavelet extraction.This method has taken into full account the actual attribute of seismologic record, and its embodiment has:

(1) with arma modeling mixed phase wavelet to be estimated and the synthetic non-stationary seismic wavelet of the random reflection coefficient sequence generating;

(2) seismologic record is carried out to improved generalized S-transform, theoretical analysis and the simulation experiment result show that this time-frequency conversion can access the good time-frequency spectrum of time-frequency focusing, the introducing of Time frequency Filter, has further improved the temporal resolution of time-frequency spectrum, the meticulous extraction of varitron ripple while being conducive to;

(3) adopt the method for spectrum analog from the seismologic record time-frequency spectrum after processing, to simulate the wavelet amplitude of each reflection spot, extracting wavelet amplitude with correlation method compares, it is the restriction of white noise that the method has been abandoned reflection coefficient, at time-frequency domain, carry out spectrum analog, react truly the time variation of wavelet, can estimate subtly the wavelet amplitude of corresponding interval;

(4) by time become the impact that becomes wavelet amplitude when Wiener Filter Method is removed from seismologic record, adopt on this basis the mixed-phase of two spectrometry reconstruct wavelets of Higher Order Cumulants, the hypothesis of having widened conventional wavelet extraction method wavelet zero phase or minimum phase, the sub-wave phase of extraction more tallies with the actual situation;

(5), by the spectral amplitude of said extracted and phase spectrum combination, obtain becoming when complete mixed phase wavelet.

Concrete principle is as follows:

Varitron wave amplitude spectrum when the present invention utilizes the method for time-frequency spectrum simulation to extract under constant condition, adopts two spectrometry reconstruct wavelet mixed-phases of Higher Order Cumulants to compose when hypothesis sub-wave phase, thereby obtains becoming when complete mixed-phase seismic wavelet.

The time frequency analysis of seismologic record

S conversion is when inheriting Short Time Fourier Transform and wavelet transformation advantage, has overcome Short Time Fourier Transform time frequency resolution constant and the scale factor of wavelet transformation and the shortcoming of frequency-independent.S conversion adopts the variable Gauss function relevant with frequency, and its window function is defined as follows:

$w\left(t\right)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{(-\frac{f^2{t}^2}{2})}---\left(1\right)$

Wherein, f is frequency.The variation tendency of S conversion window function is changeless, can not adjust flexibly according to actual needs.Conventional generalized S-transform adopts window function width to be with frequency f the Gauss function that inverse proportion changes, and when high band, window is narrower, can obtain higher temporal resolution, and when low-frequency range, window is wider, can obtain higher frequency resolution.These generalized S-transforms can obtain yupin effect when good for most of non-stationary signals, but still be not well positioned to meet non-stationary seismologic record, do not process desired resolution.Adopt a kind of improved generalized S-transform, its window function is as follows for this reason

$G(t,f)=\frac{1}{\sqrt{2\mathrm{\π}}q{\left|f\right|}^p}\exp (-\frac{t^2}{2q^2{f}^{2p}})---\left(2\right)$

Wherein, q, p is greater than zero regulatory factor.Width and the proportional example of frequency that from formula (2), can find out window function change, and at low frequency place, can obtain higher temporal resolution, at high frequency treatment, obtain higher frequency resolution, and this character meets seismologic record dynamic attenuation characteristic.The time-frequency spectrum of signal s (t) is:

$\mathrm{NGST}(\mathrm{\τ},f)={\∫}_{-\∞}^{\∞}s\left(t\right)\{\frac{1}{\sqrt{2\mathrm{\π}}q{\left|f\right|}^p}\exp \left(\frac{-{(\mathrm{\τ}-t)}^2}{2q^2{f}^{2p}}\right)\×\exp (-i2\mathrm{\πft})\}\mathrm{dt},f\≠0---\left(3\right)$

Utilize the method to carry out time frequency analysis to seismologic record, can obtain higher time frequency resolution, and there is good time-frequency focusing, can effectively tell the variation of seismic data medium frequency composition, be conducive to analyze more subtly stratal configuration.The method also can be carried out the inverse transformation of noenergy loss, accurately reconstitution time territory seismic signal.

The time-frequency filtering of seismologic record time-frequency spectrum

All inevitably there is dispersion phenomenon to a certain degree in time-frequency conversion, the time-frequency spectrum of obtained non-stationary seismologic record is directly processed, and has certain effect, but be not very good.In order better to adapt to seismic analysis and the explanation of thin interbed, eliminate the impact between adjacent reflection stratum, introduce a kind of Time frequency Filter the time-frequency spectrum of seismologic record is carried out to filtering processing.This wave filter is the smooth function of time and frequency, and has following characteristic:

(1) to τ sometime, when t=τ, H (t, ω)=1; When | t-τ | during → ∞, H (t, ω)=0;

(2) to ξ sometime, when ω=ξ, H (t, ω)=1; When | ω-ξ | during → ∞, H (t, ω)=0,

Be that H (t, ω) has time frequency localization character to spatial point (τ, ξ).Be defined as follows accordingly three parameters (λ, p, the v) Time frequency Filter of form

$H_{\mathrm{\τ}}(t,\mathrm{\ω})={[1+\frac{{(t-\mathrm{\τ})}^2{\mathrm{\ω}}^2}{v}]}^{\frac{v+1}{2}}\×[1+\frac{2}{\mathrm{\π}}\arctan \left(\mathrm{\λ}{\left|(t-\mathrm{\τ})\mathrm{\ω}\right|}^p\right)]---\left(4\right)$

The rear portion of function plays decay control action to front portion, different λ, the decay that p is corresponding different.λ is less, and p is larger, and H (t, ω) is more sharp-pointed, and support is less, and in practical application, general λ gets negative value, p get on the occasion of.

Utilize above-mentioned wave filter to carry out filtering to the time-frequency spectrum of seismologic record

S(τ,f)＝∫NGST(t _r,f)H _T(τ-t _r,f)dt _r,f≠0 (5)

Wherein S (τ, f) is filtered seismologic record time-frequency spectrum, the time window width that T is Time frequency Filter, t _rfor the time point of reflection coefficient appearance, and T/2 is less than adjacent reflection coefficient minimum interval half.

Based on time-frequency spectrum simulation time become wavelet amplitude and extract

The spectral amplitude of spectrum analog hypothesis seismic wavelet is level and smooth, thinks that wavelet amplitude is similar to the unimodal smooth curve of Ricker wavelet spectrum, and it is carried out to mathematical modeling in frequency field, and mathematic(al) representation is as follows:

$W\left(f\right)={\left|f\right|}^k\exp \left(\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{a}_n{f}^n\right)---\left(6\right)$

Wherein, k is a constant, and N is exponent number, a _nit is the multinomial coefficient about f.Remember that a certain moment wavelet amplitude is W (τ, f), adopt the method for least square fitting from seismologic record time-frequency spectrum, to obtain its estimated value.

The reconstruct of mixed-phase spectrum

Suppose that in the poststack seismologic record after processing, time-varying characteristics only show on wavelet amplitude, constant when phase place is relative.The application of considering high-order statistic is to be based upon on the basis of seismologic record stationarity hypothesis, before extracting sub-wave phase by time become Wiener Filter Method and from seismologic record, remove wavelet amplitude, seismologic record after processing is approximately stably, then adopts two spectrometry reconstruct sub-wave phases of Higher Order Cumulants.If the nonstationary random process that x (k) is zero-mean, its three rank semi-invariant is

c _x(τ ₁,τ ₂)＝E[x(k)x(k+τ ₁)x(k+τ ₂)] (7)

Pass through Fourier transform

$\begin{array}{c}{B}_x({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)=\underset{{\mathrm{\τ}}_1=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}\underset{{\mathrm{\τ}}_2=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}{c}_x({\mathrm{\τ}}_1,{\mathrm{\τ}}_2)e^{-j({\mathrm{\ω}}_1{\mathrm{\τ}}_1+{\mathrm{\ω}}_2{\mathrm{\τ}}_2)}=X\left({\mathrm{\ω}}_1\right)X\left({\mathrm{\ω}}_1\right)X^{*}({\mathrm{\ω}}_1+{\mathrm{\ω}}_2)\\ =\left|B_x({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)\right|\exp \left(\mathrm{j\ψ}({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)\right)\end{array}---\left(8\right)$

Obtain the three rank spectrums of x (k), i.e. two spectrums.The phase spectrum of this pair of spectrum is:

ψ(ω ₁,ω ₂)＝φ(ω ₁)+φ(ω ₂)-φ(ω ₁+ω ₂) (9)

Here φ (ω) is the phase spectrum of x (k).

In seismic prospecting, seismic reflection is recorded x (t) and can be expressed as

x(t)＝r(t)*w(t)(10)

Wherein r (t) is reflection coefficient, and w (t) is seismic wavelet.According to the character of two spectrums, have:

B _x(ω ₁,ω ₂)＝B _r(ω ₁,ω ₂)B _w(ω ₁,ω ₂) (11)

If reflection coefficient is Gauss's white noise sequence, B _r(ω ₁, ω ₂)=σ ², σ is variance.Can be determined by two spectrums of seismologic record the frequency spectrum of wavelet.The frequency spectrum of wavelet is designated as

Here with | W (ω) | be respectively phase spectrum and the spectral amplitude of w (t).From (9) formula

According to the symmetry of two spectrums, know | ψ (ω ₁, ω ₂) | symmetrical, only need N/2 independent formula.For convenient, calculate, make ω ₁=2k π/N, ω ₂=2l π/N, N is sampling number,

Get k=1,2 ..., N/2, corresponding l=k, k+1 ..., N-k, thus following equation obtained

Above-mentioned equation is write as to the form of matrix

$\widetilde{\mathrm{\ψ}}=A\widetilde{\mathrm{\φ}}---\left(16\right)$

Wherein,

$\widetilde{\mathrm{\ψ}}={(\mathrm{\ψ}\left(\mathrm{1,1}\right),\mathrm{\ψ}\left(\mathrm{1,2}\right),\·\·\·,\mathrm{\ψ}\left(\mathrm{2,2}\right),\mathrm{\ψ}\left(\mathrm{2,3}\right),\·\·\·,\mathrm{\ψ}(N/2,N/2))}^T$

$A=\left[\begin{array}{ccccccccc}2& -1& 0& 0& & 0& \·\·\·& & 0\\ 1& 1& -1& 0& & 0& \·\·\·& & 0\\ 1& 0& 1& -1& & 0& \·\·\·& & 0\\ & & \·& & & & & & \·\\ & & \·& & & & & & \·\\ & & \·& & & & & & \·\\ 1& 0& 0& 0& & 0& \·\·\·& 1& -1\\ 0& 2& 0& -1& & 0& \·\·\·& & 0\\ 0& 1& 1& 0& & -1& \·\·\·& & \\ & & \·& & & & & & \·\\ & & \·& & & & & & \·\\ & & \·& & & & & & \·\\ & & & & \·\·\·& & & 0& -1\end{array}\right]$

A is (N/2) ²the non-singular matrix of * N, can obtain by least square method solution be

$\widetilde{\mathrm{\φ}}={\left(A^{T}A\right)}^{-1}{A}^T\widetilde{\mathrm{\ψ}}---\left(17\right)$

Obtain the mixed-phase spectrum of wavelet.

So far, varitron wave amplitude spectrum and phase spectrum when we estimate respectively from seismologic record, become mixed phase wavelet by the two when complete in conjunction with just obtaining.

Claims (2)

1. Based on time-frequency spectrum simulation time become the thought of method of mixed phase wavelet extraction, its feature one is, the method contains following steps successively:

   Step (1) is carried out improved generalized S-transform time frequency analysis to non-stationary seismologic record:

   $\mathrm{NGST}(\mathrm{\τ},f)=\underset{-\∞}{\overset{\∞}{\∫}}s\left(t\right)\{\frac{1}{\sqrt{2\mathrm{\π}}q{\left|f\right|}^p}\exp \left(\frac{-{(\mathrm{\τ}-t)}^2}{{2q}^2{f}^{2p}}\right)\×\exp (-i2\mathrm{\πft})\}\mathrm{dt},f\≠0$

   Wherein s (t) is non-stationary seismologic record, window function for improved generalized S-transform;

   Step (2) is introduced Time frequency Filter the time-frequency spectrum of seismologic record is carried out to filtering processing:

   $S(\mathrm{\τ},f)={\∫}_{-\∞}^{\∞}\mathrm{NGST}(t,f)H_T(t-t_r,f){\mathrm{dt}}_r,f\≠0$

   Wherein H (t, f) is Time frequency Filter, the time window width that T is Time frequency Filter, t _rtime point for reflection coefficient appearance.Time frequency Filter is expressed as follows:

   $H_{\mathrm{\τ}}(t,\mathrm{\ω})={[1+\frac{{(t-\mathrm{\τ})}^2{\mathrm{\ω}}^2}{v}]}^{\frac{v+1}{2}}\×[1+\frac{2}{\mathrm{\π}}\arctan \left(\mathrm{\λ}{\left|(t-\mathrm{\τ})\mathrm{\ω}\right|}^p\right)]$

   This wave filter is the time of three parameters (λ, p, v) and the smooth function of frequency;

   The value of step (3) set time τ obtains putting sometime the frequency spectrum S (T of the seismologic record at place, f), the method of employing spectrum analog simulates the wavelet amplitude W (T in this moment, f), the position occurring according to reflection horizon changes the value of time τ, obtains the wavelet amplitude at corresponding reflection spot place;

   Step (4) becomes the impact that becomes wavelet amplitude when Wiener Filter Method is removed from seismologic record time-frequency spectrum while utilizing, residual sub-wave phase in the seismologic record after processing;

   Step (5) is supposed when sub-wave phase is constant, utilizes the mixed-phase spectrum φ (f) of two spectrometries estimation wavelets of Higher Order Cumulants;

   Step (6) is varitron wave amplitude spectrum and phase spectrum combination when above-mentioned, obtain becoming when complete mixed phase wavelet, and the impact that the mixed-phase that step (5) is obtained is further removed sub-wave phase for step (4) just obtains reflection coefficient sequence.
2. Based on time-frequency spectrum simulation time become method of mixed phase wavelet extraction, its feature two is that the method has fully retained the non-stationary composition in seismologic record, and do not need reflection coefficient sequence to carry out white noise hypothesis, while having solved conventional stage extraction, varitron wave method is in the situation that there is larger variation in adjacent layer cross-talk ripple attribute, wavelet is estimated inaccurate problem, realizes and becomes seismic wavelet estimation when meticulousr; Adopt the mixed-phase of the two spectrometries extraction of high-order wavelet, widened the hypothesis of conventional wavelet extraction method neutron deficiency zero phase or minimum phase, the sub-wave phase attribute of extraction is more realistic; Aspect raising seismic data resolution, the method does not need the dynamic attenuation mechanism of Study of Seismic record, directly seismic data is processed, and has avoided error accumulation; The introducing of Time frequency Filter, has improved temporal resolution greatly, while having eliminated time frequency analysis, due to the impact existing between the adjacent reflection horizon of dispersion phenomenon, carries out more accurately on this basis follow-up work for the treatment of.