CN104932009A - Method for enhancing resolution of seismic section through compensating Morlet wavelet transform complex time-frequency spectrum - Google Patents

Abstract

The invention relates to a method for enhancing resolution of a seismic section through compensating Morlet wavelet transform complex time-frequency spectrum, and aims at realizing seismic signal time-frequency decomposition and reconstruction. Morlet wavelet is improved by the method, based on which a Morlet wavelet inverse transform formula is constructed. The method comprises the steps that mean amplitude spectrum of the seismic section is analyzed, an attenuation frequency band requiring energy compensation is confirmed, and an amplitude spectrum attenuation function changing with frequency is obtained through fitting so that a time-frequency spectrum compensation factor and a compensation function are constructed on the basis of the attenuation function. Then the complex time-frequency spectrum is adjusted by the compensation function. Finally, a seismic signal is reconstructed through Morlet wavelet inverse transform so that an objective of enhancing resolution of the seismic signal is realized. The method has fidelity and a great resolution enhancement function.

Description

When compensating Morlet wavelet transformation multiple-frequency spectrum improves the method for seismic section resolution

Technical field

The invention belongs to geophysics and oil-gas seismic exploration development field, be specifically related to a kind of by the Morlet wavelet transformation to seismic trace multiple time-frequency spectrum compensates to improve the method for vertical seismic profiling (VSP) resolution.

Background technology

In oil-gas exploration, high-resolution geological data is needed to the fine description of oil reservoir and geologic structure.But, due to the glutinousness of the earth stratum media, make seismic event be subject to the earth when underground propagation and absorb and frequency dispersion effect, relatively reduce the energy of seismic event HFS, cause seismic signal resolution to reduce.

Researchist proposes the method for various raising seismic exploration signal resolution.In early days, based on simple convolution model, suppose that geological data is stationary signal, propose deconvolution algorithms to improve the resolution of geological data.Subsequently, people recognize that there is attenuation by absorption effect on stratum to seismic event, cause geological data to be non-stationary signal, thus the various non-stationary the Method of Deconvolution proposed and inverse Q filtering method improve the resolution of seismic signal.In document " seismic data process application technology " in March, 2008 petroleum industry publishing house 306 pages, think that non-stationary the Method of Deconvolution is devoted to the research of the mathematical model of deconvolution itself, relax assumed condition, and inverse Q filtering method is by improving the data characteristic before deconvolution, make it to adapt to deconvolution requirement better.

Gabor transformation is used for improving seismic resolution by the article " Gabor deconvolution " that document Margrave and Lamoureux (2001) rolls up 241-276 in 13 at " CREWES Research Report ", a new non-stationary the Method of Deconvolution is proposed, method is white spectrum based on reflection coefficient, and through the source wavelet of attenuation by earth absorption be minimum phase these two hypothesis, method realize need estimate time become seismic wavelet.Many scholars are had to improve the method subsequently.And usually through the source wavelet of attenuation by earth absorption may be mixed-phase.Hale (1981) proposes inverse Q filtering the earliest, has many scholars to carry out various improvement to the method equally.Inverse Q filtering can compensate signal in frequency, amplitude, phase place etc., and its shortcoming needs to know accuracy Q value.

This invention improves the new method of seismic resolution to one.Method by seismic signal time-frequency spectrum regulates, and to compensate seismic signal energy attenuation, reaches the object improving seismic resolution, wherein, seismic signal time-frequency spectrum obtained by Morlet wavelet transformation.New method does not need to calculate Q value, does not relate to estimation and the minimum phase hypothesis of source wavelet.Method only suppose reflection coefficient be fluctuation change white spectral sequence, and source wavelet be one containing the short signal enriching frequency content.

Summary of the invention

The technical matters solved

Method mainly inverse Q filtering and the non-stationary deconvolution of existing raising seismic exploration signal resolution.Inverse Q filtering method needs to calculate Q value, and the acquisition of Q value is a more difficult job.What non-stationary deconvolution was more famous has Gabor deconvolution, needs to become seismic wavelet when needing to estimate under the hypothesis of minimum phase.In order to avoid the deficiencies in the prior art part, the present invention propose a kind of Morlet wavelet transformation by compensating seismic signal multiple time-frequency spectrum to be to improve the new method of seismic exploration signal resolution.Method hypothesis stratum reflection coefficient sequence be white spectrum, and the method for nearly all raising seismic exploration signal resolution all adopts this to suppose, and suppose source wavelet be one containing the short signal enriching frequency content.

Technical scheme

Seismic exploration signal is a reflection of underground structure, containing abundant geophysical information.But due to the glutinousness of earth medium, when underground propagation, can there is energy attenuation and dispersion phenomenon in seismic exploration wavelet, particularly high-frequency information energy attenuation is serious.Thus, cause seismic signal resolution to decline, the less stratal configuration information relevant to high-frequency information can not be reflected.This situation is unfavorable for that geologist carries out seismic signal explanation, and does fine description to underground structure further.

The attenuation by absorption effect of the earth is equivalent to a wave filter.The radio-frequency component of this wave filter to source wavelet has strong attenuation by absorption effect, thus the underground structure information that radio-frequency component is carried can not show in the seismic data.But these information does not disappear, and is just covered by high-energy component or noise, can not dominantly show.Compensating by finding the high-frequency energy information of appropriate Method of Mathematical Physics to attenuation by earth absorption, can make that the structural information entrained by radio-frequency component is dominant to be represented, to reach the object improving seismic exploration signal resolution.

Suppose that seismic wavelet w (t) is an impact signal in short-term, containing abundant frequency content.Stratum reflection coefficient r (t) is a white spectrum signal changed greatly.If do not have earth filtering effect, the desirable seismic exploration signal of wave detector pickup is

s ₀(t)＝w(t)*r(t) (1)

In formula (1), symbol " * " is convolution operation symbol.Owing to being subject to the effect of attenuation by earth absorption, the actual seismic signal obtained is

s(t)＝α(t)*s ₀(t)+n(t) (2)

In formula (2), function alpha (t) embodies the attenuation by absorption effect of the earth, and n (t) is noise.Do Fourier conversion at formula (2) two ends, obtain its Spectrum Relationship

$\hat{s}\left(f\right)=\widehat{\mathrm{\α}}\left(f\right){\hat{s}}_0\left(f\right)+\hat{n}\left(f\right)---\left(3\right)$

Get spectral amplitude at formula (3) two ends, and be reduced to

$A_{\hat{s}}\left(f\right)=\mathrm{At}\left(f\right)A_{{\hat{s}}_0}\left(f\right)---\left(4\right)$

In formula (4), the Fourier transforming amplitudes spectrum of s (t), s ₀t the Fourier transforming amplitudes spectrum of (), At (f) is a synthesis result of noise and the earth decay, is called amplitude decay function.Because underground structure is complicated and changeable, amplitude decay function At (f) is difficult to determine, the present invention adopts statistical method, is obtained by actual seismic section mean amplitude spectrum.Amplitude decay function structure compensating factor function thus, and actual seismic section is compensated.

First, compensation makes actual seismic section spectral amplitude return to perfect condition

$A_{{\hat{s}}_0}\left(f\right)=A_{\hat{s}}\left(f\right)/\mathrm{At}\left(f\right)---\left(5\right)$

Secondly, stratal configuration is complicated and changeable, and therefore the attenuation by absorption rule of seismic signal is complicated.Because the decay of seismic signal changes in time, when thus restoration and compensation process should have, become feature; Again due to different to different frequency content attenuation by absorption effects on synchronization stratum, thus restoration and compensation also should have with frequency change character.Therefore, the attenuation compensation method that the present invention provides realizes at the time and frequency zone of seismic signal, by Morlet wavelet transformation to seismic exploration signal time-frequency spectrum does suitable adjustment, to play the effect of attenuation compensation, thus reaches the object of the resolution improving seismic signal.

When compensating Morlet wavelet transformation multiple-frequency spectrum improves the method for seismic section resolution, it is characterized in that step is as follows:

Step 1: be FFT to each seismic trace in seismic section and obtain spectral amplitude, ask for the mean amplitude spectrum of seismic section, represents mean amplitude spectrum frequency-decible chart, finds out the frequency F that maximum amplitude point is corresponding in frequency-decible chart _m, will F be greater than according to spectral amplitude variation tendency _mfrequency range part be divided into three sections: first paragraph, adjacent F _mfrequency, its spectral amplitude amplitude fluctuates up and down, does not have obvious downtrending, is called oscillation frequency bands, and this frequency range does not need to compensate; Second segment, its spectral amplitude amplitude has one to decline fast, is called decay frequency range, and this frequency range needs to compensate; 3rd section, its spectral amplitude amplitude downtrending is flat slow, and amplitude is very little, is called high noisy frequency range, and the Signal-to-Noise of this frequency range is lower, when should not carry out-and frequency spectrum compensation;

Step 2: do linear fit to decay frequency range and obtain spectral amplitude attenuation function At (f)=af+b, wherein a and b is the parameter of matching, during structure-the frequency spectrum compensation factor:

$\mathrm{CF}\left(f\right)=\exp (-\frac{\mathrm{At}\left(f\right)}{10})=\exp (-\frac{\mathrm{af}+b}{10});f\∈[f_1,f_n]$

$\mathrm{CF}\left(f\right)=1.0;f\∉[f_1,f_n]$

Wherein, f ₁for the beginning frequency point of the frequency range that decays, f _nfor the ending frequency point of the frequency range that decays;

Step 3: by the seismic trace contained by seismic section by road carry out Morlet wavelet transformation obtain seismic trace multiple time-spectral matrix:

Described Morlet wavelet transformation expression formula is:

$\begin{array}{c}J(n,m)\stackrel{\mathrm{def}}{=}\mathrm{JT}(\mathrm{n\ΔT},\mathrm{m\ΔF})\\ ={\mathrm{\Σ}}_{k=0,-K\≤n-k\≤K}^N\mathrm{Tr}\left(\mathrm{k\ΔT}\right)\frac{\left|\mathrm{m\ΔF}\right|}{\sqrt{2\mathrm{\π}}}\exp \{-\frac{{\left(\mathrm{m\ΔF}\right)}^2{\left[(n-k)\mathrm{\ΔT}\right]}^2}{2}\}\exp \left[i2\mathrm{\πm\ΔF}(n-k)\mathrm{\ΔT}\right]\end{array}$

$n=\mathrm{0,1},...,N;m=\mathrm{0,1},...,M$

$J(n,0)\stackrel{\mathrm{def}}{=}\mathrm{JT}(\mathrm{n\ΔT},0)=0,n=\mathrm{0,1},...,N$

In formula, Tr () is seismic trace, the kth time-sampling point that Tr (k Δ T) is seismic trace, altogether N+1 sampled point, and Δ T is time sampling interval; Δ F is frequency sampling interval, when n represents-and temporal sampled point sequence number, when m represents-spectral frequencies sampled point sequence number, frequency sampling is always counted as M+1;

Step 4: when calculating multiple-mould of each element of spectral matrix:

Find out the maximum modulus value of often row, obtain maximum modulus value vector (M (0) M (1) ... M (N)) ^t, during structure-frequency spectrum compensation function:

$\mathrm{CFu}(n,m)=\frac{{\left[\right|\left|J(n,m)\right||\×\mathrm{CF}\left(m\right)]}^p}{\left|\right|J(n,m)\left|\right|+\mathrm{\ϵ}\×M\left(n\right)},n=\mathrm{0,1},...,N;m=\mathrm{0,1},...,M;$

Wherein, ε is regularization factors, and p is smoothing factor,

Step 5: used time-frequency spectrum compensation function to time multiple-spectral matrix compensates by row, column:

Step 6: adopt inverse Morlet wavelet transformation to multiple after compensating time-frequency spectrum reconfiguration seismic trace;

Described inverse Morlet wavelet transformation expression formula is:

$\begin{array}{c}\widetilde{\mathrm{Tr}}\left(n\right)\stackrel{\mathrm{def}}{=}\mathrm{Tr}\left(\mathrm{n\ΔT}\right)\\ {={\mathrm{\Σ}}_{m=0}^M\mathrm{\Σ}}_{k=0,-N\≤n-k\≤N}^N\mathrm{CJ}(k,m)\exp \left[j2\mathrm{\π}(n-k)\mathrm{\ΔTm\ΔT}\right]\\ n=\mathrm{0,1},...,N\end{array}$

Step 7: carry out normalization by the amplitude of trace equalization standardized method to reconstruct seismic trace and correct, the seismic trace before making its order of magnitude and converting is consistent: find out seismic traces sampling point amplitude and equal maximal value S _msampling point number, be designated as K, before finding out the seismic channel data of reconstruct, the large value of K directly replaces with S _m, its K large value is designated as T _k, to other data, standardize by following formula:

$\mathrm{NT}\left(n\right)=S_M\×\widetilde{\mathrm{Tr}}\left(n\right)/T_K.$

Described ε span 0.001 >=ε > 0.

Described p span 1 >=p > 0.

Beneficial effect

The present invention propose a kind of by compensating Morlet wavelet transformation multiple time-frequency spectrum improves the method for seismic section resolution, can significantly improve the resolution of seismic exploration signal.By actual seismic exploration data verification this method, there is fidelity and stability.Method And Principle is clear, and realize simple, user operation is few.

Accompanying drawing explanation

Fig. 1 is original seismic section, through strengthening resolution processes;

Fig. 2 is original seismic section spectral amplitude and mean amplitude spectrum;

Fig. 3 is the frequency-decibel curve of original seismic section mean amplitude spectrum;

Fig. 4 is that the frequency range of original seismic section mean amplitude spectrum divides.Determine according to this frequency-decibel curve figure the frequency range that decays, and simulate attenuation function;

Fig. 5 is the seismic section after attenuation compensation.Compared with the original seismic section shown in Fig. 1, after compensating, the resolution of seismic section is significantly improved

Fig. 6 is the spectral amplitude of the seismic section after attenuation compensation

Embodiment

Now in conjunction with the embodiments, the invention will be further described for accompanying drawing:

The invention provides following technical scheme, by time multiple to the Morlet wavelet transformation of seismic signal-frequency spectrum compensates, play and earthquake signal energy decay the effect compensated, thus realize the object of raising seismic resolution.Method improves Morlet small echo, and the Morlet Construction of Wavelets improved based on this goes out its inverse transformation.Implementation is as follows:

1) Morlet small echo does not meet small echo admissible condition, thus, does not have wavelet reconstruction formula.Morlet small echo is a special case of Gabor wavelet, and Gabor wavelet function is

$\mathrm{\ψ}\left(t\right)=\frac{1}{{\left({\mathrm{\σ}}^2\mathrm{\π}\right)}^{1/4}}{e}^{-\frac{t^2}{2{\mathrm{\σ}}^2}}{e}^{\mathrm{j\ηt}}---\left(6\right)$

In formula (6), when σ=1, Gabor wavelet is Morlet small echo.The present invention gets η=2 π, and amendment Morlet mother wavelet function is

$m\left(t\right)=\frac{1}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{t^z}{z}}{e}^{j2\mathrm{\πt}}---\left(7\right)$

In formula (7), translation is carried out to time parameter t or change of scale can obtain Morlet wavelet function race { m _{f, τ}(t) }

$m_{f,\mathrm{\τ}}\left(t\right)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^z{(t-\mathrm{\τ})}^z}{z}}{e}^{j2\mathrm{\πf}(t-\mathrm{\τ})}---\left(8\right)$

In formula (8), parameter f is the frequency window center frequency of the time frequency window of Morlet small echo; τ be the time frequency window of Morlet small echo time the window center time.If seismic signal is Tr (t), its Morlet wavelet transformation is

$\begin{array}{c}\mathrm{JT}(\mathrm{\τ},f)={\∫}_{-\∞}^{+\∞}\mathrm{Tr}\left(t\right)\stackrel{\‾}{m_{f,\mathrm{\τ}}}\left(t\right)\mathrm{dt}\\ ={\∫}_{-\∞}^{+\∞}\mathrm{Tr}\left(t\right)\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^z{(t-\mathrm{\τ})}^z}{z}}{e}^{-j2\mathrm{\πf}(t-\mathrm{\τ})}\mathrm{dt}\\ =\mathrm{Tr}\left(\mathrm{\τ}\right)*m_{f,0}\left(\mathrm{\τ}\right)\end{array}---\left(9\right)$

In formula (9), m _{f, τ}the complex conjugate function of (t), operational symbol " * " convolution operation.

2) to signal reconstruction be realized, must inverse transformation be had, i.e. reconstruction formula.Because Morlet small echo does not meet small echo admissible condition, not allow small echo, thus there is no the reconstruction formula of traditional small echo.The present invention provides the inverse transformation being only suitable for Morlet small echo

$\mathrm{Tr}\left(t\right)={\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}\mathrm{JT}(\mathrm{\τ},f)e^{-j2\mathrm{\πf}(t-\mathrm{\τ})}\mathrm{dfd\τ}---\left(10\right)$

Conveniently numerical operation, formula (10) is changed below can carrying out

$\begin{array}{c}\mathrm{Tr}\left(t\right)={\∫}_{-\∞}^{+\∞}[\mathrm{JT}(t,f)*e^{j2\mathrm{\πft}}]\mathrm{dt}\\ ={\∫}_{-\∞}^{+\∞}M(t,f)\mathrm{df}\end{array}---\left(11\right)$

In formula (11), note

M(t，f)＝JT(t，f)*e ^j2πft(12)

When carrying out seismic signal-frequency division solution and with after compensating time-frequency spectrum reconfiguration seismic signal time, need to formula (9), (11) and (12) discretize, to realize numerical operation.

The present invention also provides following technical scheme, and method is compensated by the high-frequency energy greatly absorbed seismic signal, to reach the object improving seismic signal resolution.Rational compensation must understand the attenuation by absorption Changing Pattern of seismic signal.Due to the complicacy of underground structure and the otherness of seismic exploration data acquisition and processing (DAP) mode, the present invention adopts statistical method to obtain the technical scheme of the rule that seismic signal energy changes with frequency decay.Implementation method is as follows:

1) FFT is done to each seismic trace in seismic section, obtain the spectral amplitude of seismic trace;

2) to the spectral amplitude of all seismic traces, cumulative, on average, obtain a mean amplitude spectrum;

3) maximum amplitude of note mean amplitude spectrum is A _e, frequency corresponding to this amplitude is designated as F _e, be called maximum amplitude frequency.With maximum amplitude frequency for benchmark, make the frequency-decible chart of mean amplitude spectrum;

4) with maximum amplitude frequency for branch, the frequency of spectral amplitude is divided into two sections, is called anterior frequency range and rear portion frequency range.Anterior frequency range, namely frequency is less than F _efrequency range, information contained by the frequency can thinking in it comparatively speaking energy does not decay, and the signal energy of this frequency range does not need to compensate.Rear portion frequency range, is namely greater than F _efrequency range, information contained by the frequency can thinking in it comparatively speaking energy has comparatively high attenuation, and its energy demand has to be selected to compensate.

5) according to the Changing Pattern of spectral amplitude, rear portion frequency range roughly can divide three sections.Three types can be divided into, by distance maximum amplitude frequency F according to spectral amplitude Changing Pattern _edistance be divided into three sections.First paragraph, adjacent F _e, in this frequency range, amplitude is vibrated up and down, does not show attenuation trend, and this frequency range need not compensate, and is called oscillation frequency bands; Second segment, spectral amplitude significantly decreases trend, and this frequency range needs to compensate, and is called decay frequency range; 3rd section, apart from F _efarthest, spectral amplitude general trend is still for downtrending, but downtrending is mild and spectral amplitude amplitude is very low, and noise element is very high, and this frequency range, without the need to compensating, is called high noisy frequency range.In practical problems, decay frequency range may have the segment of different downtrending, can segment again.High noisy frequency range compensates the signal to noise ratio (S/N ratio) by reducing signal to it, reduces seismic signal resolution on the contrary, therefore, without the need to compensating, even can reject;

6) the decay frequency range of rear portion frequency range is the frequency range needing to compensate.According to actual conditions, linear fit can be done, or by different downtrending, sectional linear fitting, obtains spectral amplitude attenuation function.When then constructing with spectral amplitude attenuation function-frequency spectrum compensation Summing Factor penalty function.

The present invention also provides following technical scheme: by road Morlet wavelet transformation, when obtaining seismic trace multiple-and frequency spectrum; Time-adjustment of frequency spectrum is subject to the constraint of time and frequency Two Variables, that is: what become when compensation is is also frequently become; To multiple after compensating time-frequency spectrum, with the inverse transformation of Morlet small echo reconstruct seismic signal.Comprise the following steps:

The first step, Morlet wavelet transformation is done to current seismic road, when obtaining current seismic road multiple-spectral matrix, when then calculating multiple-each element modulus value of spectral matrix, obtain its modular matrix, when being called-spectrogram, and when searching out this time point maximum on each time point-frequency modulus value M (t);

Second step, to time multiple-spectral matrix in each element, used time-frequency spectrum compensation function regulates, be compensated rear multiple time-spectral matrix;

3rd step, when after being compensated by the conversion of inverse Morlet small echo, seismic trace is multiple-spectral matrix reconstruct seismic trace;

4th step, the seismic trace of normalization reconstruct, obtains high-resolution new seismic trace.

The present invention adopts following steps to realize improving seismic section resolution:

The first step, data analysis.

1) maximum frequency range [0, the F that can be showed by time-domain sampling rate and the Nyquist sampling theorem determination seismic signal of cross-sectional data _max], wherein, F _maxfor the highest frequency of signal;

2) for seismic section, be FFT by road, and statistics obtains the mean amplitude spectrum of seismic trace in section, the maximal value of the spectral amplitude of search section Shang Ge road FFT and mean amplitude spectrum maximum amplitude point (A _e, F _e), wherein, A _emaximum average amplitude value, F _eit is the frequency that maximum average amplitude value is corresponding;

3) each channel amplitude spectrogram of normalization is drawn, and the frequency-decible chart of mean amplitude spectrum;

4) with the frequency F that maximum average amplitude is corresponding _ethe frequency of seismic signal is divided into two sections: anterior frequency range [0, F _e], and rear portion frequency range [F _e, F _max].In invention, the information energy loss of the frequency that anterior frequency range comprises corresponding to it does not need to compensate.Rear portion frequency range can be divided into three parts, by arriving F _edistance is designated as according to spectral amplitude feature: oscillating part, attenuation portions and high noisy part near respectively to far, wherein, and the adjacent F of oscillating part _e, spectral amplitude corresponding to this part frequency vibrates up and down, and do not have obvious downtrending, its amplitude is between 0 ~-8dB, and-4dB horizontal line of being everlasting fluctuates up and down, and the information of this part frequency does not need to compensate; Attenuation portions, be positioned in the middle of the frequency range of rear portion, the spectral amplitude of this band frequency shows quick downward trend, span drops to about-35dB from about-4dB, the information of this part frequency needs to compensate, this frequency range middle ware or have little fluctuation or frequency range two ends scope slightly to change to have no significant effect method effect; High noisy part, is positioned at rear portion frequency range most end, and amplitude is under-30dB, and spectral amplitude amplitude is very little, and this section of spectral amplitude downtrending is mild, and energy ezpenditure is maximum, almost can regard noise as.The information that this band frequency is corresponding, has lower signal to noise ratio (S/N ratio), if compensation can increase signal noise, reduces resolution, and the information of this part frequency does not need to compensate yet.

Second step, Amplitude spectrum attenuation function, during structure-the frequency spectrum compensation factor.

Although all frequency contents of seismic signal all can be subject to the attenuation by absorption effect on stratum in communication process, only consider relative attenuation herein, therefore, with frequency F _ecorresponding amplitude A _efor benchmark, think only there is the first step 4) described in rear portion frequency range [F _e, F _max] in the information of " attenuation portions " and " high noisy part " this two band frequency have relative energy to decay.

In the present invention, energy compensating is only implemented " attenuation portions ", and therefore, this step only makes linear fit or sectional linear fitting to the curve of " attenuation portions " in frequency-decible chart.If curve point is

(f ₁，A ₁),(f ₂，A ₂),…,(f _n，A _n) (13)

Wherein, f _ifor frequency, A _ifor the decibel amplitude of its correspondence, i=1,2 ..., the linear fit function of n. then spectral amplitude attenuation function is

At(f)＝af+b (14)

In formula (14), independent variable is frequency f, a and b is the parameter of matching.Time corresponding-the frequency spectrum compensation factor is

$\mathrm{CF}\left(f\right)=\exp (-\frac{\mathrm{At}\left(f\right)}{10})=\exp (-\frac{\mathrm{af}+b}{10});f\∈[f_1,f_n]---\left(15\right)$

Because " attenuation portions " variation tendency may be inconsistent, sectional linear fitting can be done by variation tendency, now, time-the frequency spectrum compensation factor will be piecewise linear function.When can be found out by formula (15)-the frequency spectrum compensation factor is the function of frequency.To the frequency range without the need to compensating, compensating factor is constant 1, namely has

$\mathrm{CF}\left(f\right)=1.0;f\∉[f_1,f_n]---\left(16\right)$

3rd step, when seismic section is undertaken by road-frequency spectrum compensation.

To the energy attenuation of seismic section compensate Shi Zhu road realize.First, when seismic trace Morlet wavelet transformation is done-frequency division solution, obtain its multiple time-frequency spectrum; Then, during structure-frequency spectrum compensation function, when regulating seismic trace multiple with penalty function-frequency spectrum; Finally, with the inverse transformation expression re-formation seismic trace of Morlet small echo.Specific implementation step is as follows:

0) common parameter is set for the entire profile process: smoothing factor p, and 1>=p > 0; Regularization factors ε, and 0.001>=ε > 0; Pending geological data maximal value S is obtained by geological data header file _m, road sampling number (N+1) and time sampling interval Δ T.By Nyquist Sampling Theorem, the largest interval Δ F of calculated rate sampling _max, and as required set handling time the frequency sampling interval delta F that uses, and have Δ F≤Δ F _max, calculated rate sampling number, is designated as (M+1);

1) get f=m Δ F, fixed frequency f, by m _{f, 0}t () be discretize in time domain, have

$\begin{array}{c}m(n,m)\stackrel{\mathrm{def}}{=}{m}_{\mathrm{m\ΔF},0}\left(\mathrm{n\ΔT}\right)=\frac{\left|\mathrm{m\ΔF}\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{{\left(\mathrm{m\ΔF}\right)}^2{\left(\mathrm{n\ΔT}\right)}^2}{2}]\exp \left(i2\mathrm{\πm\ΔFn\ΔT}\right)\\ n=-K,-K+1,...,-\mathrm{1,0,1},...,K;m=\mathrm{0,1},...,M\end{array}---\left(17\right)$

Utilize formula (17), discrete type (9), and note have

$\begin{array}{c}J(n,m)\stackrel{\mathrm{def}}{=}\mathrm{JT}(\mathrm{n\ΔT},\mathrm{m\ΔF})\\ ={\mathrm{\Σ}}_{k=0,-K\≤n-k\≤K}^N\mathrm{tr}\left(k\right)m(n-k,m)\\ n=\mathrm{0,1},...,N\end{array}---\left(18\right)$

$J(n,0)\stackrel{\mathrm{def}}{=}\mathrm{JT}(\mathrm{n\ΔT},0)=0,n=\mathrm{0,1},...,N---\left(19\right)$

When can answer-spectral matrix

Wherein, J (n, m) be plural number (n=0,1 ..., N; M=0,1 ..., M).When calculating multiple-modular matrix of spectral matrix, and when being called-spectrogram matrix, namely

Temporally, namely by row, find out the maximum modulus value of each time point, have maximum modulus value vector

(M(0) M(1) … M(N) ^T(22)

2) according to the frequency spectrum compensation factor, formula (15) and (16), during structure-frequency spectrum compensation function,

$\begin{array}{c}\mathrm{CF}\left(m\right)\stackrel{\mathrm{def}}{=}\mathrm{CF}\left(\mathrm{m\ΔF}\right)\\ \mathrm{CFu}(n,m)=\frac{{\left[\right|\left|J(n,m)\right||\×\mathrm{CF}\left(m\right)]}^p}{\left|\right|J(n,m)\left|\right|+\mathrm{\ϵ}\×M\left(n\right)}\\ n=\mathrm{0,1},...,N;m=\mathrm{0,1},...,M\end{array},---\left(23\right)$

With formula (23) to multiple shown in (20) time-spectral matrix in each element act on, during multiple after can being compensated-spectral matrix

Time after 4th step compensates again-and frequency spectrum reconfiguration seismic signal, to improve resolution.

In order to realize restructing operation, need to formula (11) and formula (12) discretize.First to index harmonic ringing discretize.Get frequency f=m Δ F, fixed frequency f, has

$\begin{array}{c}e(k,m)\stackrel{\mathrm{def}}{=}\exp \left(j2\mathrm{\πk\ΔTm\ΔF}\right)\\ k=-N,-N+1,...,-\mathrm{1,0,1},...,N;m=\mathrm{0,1},...,M\end{array}---\left(25\right)$

Can discretize formula (12) by formula (25), have

$\begin{array}{c}M(n,m)={\mathrm{\Σ}}_{k=0,-N\≤n-k\≤N}^N\mathrm{CJ}(k,m)e(n-k,m)\\ n=\mathrm{0,1},..,N;m=\mathrm{0,1},...,M\end{array}---\left(26\right)$

Can discretize formula (11) by formula (26), obtain the seismic signal rebuild have

$\begin{array}{c}\widetilde{\mathrm{Tr}}\left(n\right)\stackrel{\mathrm{def}}{=}\mathrm{Tr}\left(\mathrm{n\ΔT}\right)\\ \begin{array}{cc}={\mathrm{\Σ}}_{m=0}^{M}M(n,m)& n=\mathrm{0,1},...,N\end{array}\end{array}---\left(27\right)$

1) seismic signal is reconstructed.When answering compensation is rear-and frequency spectrum, shown in formula (24), provide method process by formula (26), have

To formula (28), with formula (27) give method, add up by row, can seismic trace be reconstructed

2) to reconstruct seismic trace go bail for width normalization process.Processing mode is: find out seismic traces sample value and get maximal value S _msampled point number, be designated as K, to the sequence of reconstruct seismic channel data, its K is large, and value is designated as T _k, directly S is replaced with to K large value before reconstruct seismic trace _m, other data, standardize by following formula

$\mathrm{NT}\left(n\right)=S_M\×\widetilde{\mathrm{Tr}}\left(n\right)/T_K---\left(29\right)$

Often press previous step operation together respectively on section, resolution can be realized and strengthen effect.

Claims (3)

1. when compensating Morlet wavelet transformation multiple-frequency spectrum improves the method for seismic section resolution, it is characterized in that step is as follows:

Step 1: be FFT to each seismic trace in seismic section and obtain spectral amplitude, ask for the mean amplitude spectrum of seismic section, represents mean amplitude spectrum frequency-decible chart, finds out the frequency F that maximum amplitude point is corresponding in frequency-decible chart _m, will F be greater than according to spectral amplitude variation tendency _mfrequency range part be divided into three sections: first paragraph, adjacent F _mfrequency, its spectral amplitude amplitude fluctuates up and down, does not have obvious downtrending, is called oscillation frequency bands, and this frequency range does not need to compensate; Second segment, its spectral amplitude amplitude has one to decline fast, is called decay frequency range, and this frequency range needs to compensate; 3rd section, its spectral amplitude amplitude downtrending is flat slow, and amplitude is very little, is called high noisy frequency range, and the Signal-to-Noise of this frequency range is lower, when should not carry out-and frequency spectrum compensation;

Step 2: do linear fit to decay frequency range and obtain spectral amplitude attenuation function At (f)=af+b, wherein a and b is the parameter of matching, during structure-the frequency spectrum compensation factor:

$\mathrm{CF}\left(f\right)=\exp (-\frac{\mathrm{At}\left(f\right)}{10})=\exp (-\frac{\mathrm{af}+b}{10});f\∈[f_1,f_n]$

$\mathrm{CF}\left(f\right)=1.0;f\∉[f_1,f_n]$

Wherein, f ₁for the beginning frequency point of the frequency range that decays, f _nfor the ending frequency point of the frequency range that decays;

Step 3: by the seismic trace contained by seismic section by road carry out Morlet wavelet transformation obtain seismic trace multiple time-spectral matrix:

Described Morlet wavelet transformation expression formula is:

$\begin{array}{c}J(n,m)\stackrel{\det }{=}\mathrm{JT}(\mathrm{n\ΔT},\mathrm{m\ΔF})\\ ={\mathrm{\Σ}}_{k=0,-K\≤n-k\≤K}^N\mathrm{Tr}\left(\mathrm{k\ΔT}\right)\frac{\left|\mathrm{m\ΔF}\right|}{\sqrt{2\mathrm{\π}}}\exp \{-\frac{{\left(\mathrm{m\ΔF}\right)}^2{\left[(n-k)\mathrm{\ΔT}\right]}^2}{2}\}\exp \left[i2\mathrm{\πm\ΔF}(n-k)\mathrm{\ΔT}\right]\\ n=\mathrm{0,1},...,N;m=\mathrm{0,1},...,M\end{array}$

$J(n,0)\stackrel{\det }{=}\mathrm{JT}(\mathrm{n\ΔT},0)=0,n=\mathrm{0,1},...,N$

In formula, Tr () is seismic trace, the kth time-sampling point that Tr (k Δ T) is seismic trace, altogether N+1 sampled point, and Δ T is time sampling interval; Δ F is frequency sampling interval, when n represents-and temporal sampled point sequence number, when m represents-spectral frequencies sampled point sequence number, frequency sampling is always counted as M+1;

Step 4: when calculating multiple-mould of each element of spectral matrix:

Find out the maximum modulus value of often row, obtain maximum modulus value vector (M (0) M (1) ... M (N)) ^t, during structure-frequency spectrum compensation function:

$\mathrm{CFu}(n,m)=\frac{\left[\right|\left|J(n,m)\right||\×\mathrm{CF}\left(m\right){]}^p}{\left|\right|J(n,m)\left|\right|+\mathrm{\ϵ}\×M\left(n\right)},n=\mathrm{0,1},...,N;m=\mathrm{0,1},...,M;$

Wherein, ε is regularization factors, and p is smoothing factor,

Step 5: used time-frequency spectrum compensation function to time multiple-spectral matrix compensates by row, column:

Step 6: adopt inverse Morlet wavelet transformation to multiple after compensating time-frequency spectrum reconfiguration seismic trace;

Described inverse Morlet wavelet transformation expression formula is:

Step 7: carry out normalization by the amplitude of trace equalization standardized method to reconstruct seismic trace and correct, the seismic trace before making its order of magnitude and converting is consistent: find out seismic traces sampling point amplitude and equal maximal value S _msampling point number, be designated as K, before finding out the seismic channel data of reconstruct, the large value of K directly replaces with S _m, its K large value is designated as T _k, to other data, standardize by following formula:

2. during compensation Morlet wavelet transformation according to claim 1 multiple-frequency spectrum improves the method for seismic section resolution, it is characterized in that described ε span 0.001 >=ε > 0.

3. during compensation Morlet wavelet transformation according to claim 1 multiple-frequency spectrum improves the method for seismic section resolution, it is characterized in that described p span 1 >=p > 0.