CN104932008B - The method that the multiple time-frequency spectrum of compensation J conversion improves seismic profile resolution ratio - Google Patents

Abstract

The present invention relates to the method that a kind of multiple time-frequency spectrum of compensation J conversion improves seismic profile resolution ratio, it is analyzed by the mean amplitude spectrum to seismic profile, determine the decay frequency range for needing energy compensating, and the amplitude spectrum attenuation function changed with frequency is fitted, time-frequency spectrum compensating factor and penalty function are constructed with this attenuation function.Then, with the multiple time-frequency spectrum of this penalty function regulation.Finally, reconstruct seismic signal is converted by inverse J, realizes improving the purpose of seismic signal resolution ratio.Method has fidelity and good enhancing resolution ratio function.

Description

Compensation J conversion it is multiple when-frequency spectrum improve seismic profile resolution ratio method

Technical field

The invention belongs to geophysics and oil-gas seismic exploration development field, it is related to a kind of by seismic prospecting signal J conversion it is multiple when-frequency spectrum compensates method to improve vertical seismic profiling (VSP) resolution ratio.

Background technology

In oil-gas exploration, the fine description to oil reservoir and geologic structure needs high-resolution geological data.However, by In the viscosity of the earth stratum media so that seismic wave is subject to the earth to absorb and frequency dispersion effect in underground propagation, relative reduction The energy of seismic wave HFS, causes the reduction of seismic signal resolution ratio.

The method that researcher proposes various raising seismic exploration signal resolutions.In early stage, based on simple convolution mould Type, it is assumed that geological data is stationary signal, proposes deconvolution algorithms to improve the resolution ratio of geological data.Then, people's understanding There is attenuation by absorption to act on seismic wave to stratum, cause geological data for non-stationary signal, so that the various non-stationaries for proposing are anti- Convolution method and inverse Q filtering method improve the resolution ratio of seismic signal.Document《Seismic data process application technology》2008 3 In month petroleum industry publishing house page 306, it is believed that non-stationary the Method of Deconvolution is devoted to the research of the Mathematical Modeling of deconvolution in itself, Relax assumed condition, and inverse Q filtering method is by improving the data characteristic before deconvolution, being allowed to better adapt to deconvolution will Ask.

Document Margrave and Lamoureux (2001) exists《CREWES Research Report》The in volume 13 Be used for for Gabor transformation to improve seismic resolution by the article " Gabor deconvolution " of 241-276, propose one it is new Non-stationary the Method of Deconvolution, it is white spectra that method is based on reflectance factor, and is minimum phase by the source wavelet of attenuation by earth absorption Position the two it is assumed that method realize need estimate time-varying seismic wavelet.Then there are many scholars to be improved the method.And lead to Often may be mixed-phase by the source wavelet of attenuation by earth absorption.Hale (1981) proposes inverse Q filtering earliest, equally has Many scholars carry out various improvement to the method.Inverse Q filtering can be compensated at aspects such as frequency, amplitude, phases to signal, its Shortcoming needs to know degree of accuracy Q values.

This invention is to a new method for raising seismic resolution.Method by seismic signal when-frequency spectrum carries out Regulation, to compensate seismic signal energy attenuation, reach improve seismic resolution purpose, wherein, seismic signal when-frequency spectrum lead to J conversion is crossed to obtain.New method need not calculate Q values, not be related to the estimation of source wavelet and minimum phase to assume.Method is only assumed that Reflectance factor is the white spectra sequence of fluctuating change, and source wavelet is a short signal containing abundant frequency content.

The content of the invention

The technical problem to be solved

The existing method for improving seismic exploration signal resolution is mainly inverse Q filtering and non-stationary deconvolution.Inverse Q filtering Method needs to calculate Q values, and the acquisition of Q values is the difficult work of a comparing.Non-stationary deconvolution it is more famous have Gabor anti- Convolution under the hypothesis of minimum phase, it is necessary to need to estimate time-varying seismic wavelet.In order to avoid the deficiencies in the prior art part, this When invention proposes that a kind of J by compensating seismic signal converts multiple-frequency spectrum to be improving the new method of seismic exploration signal resolution. Method assumes that stratum reflection coefficient sequence is white spectra, and the method for nearly all raising seismic exploration signal resolution is false using this If, and assume that source wavelet is a short signal containing abundant frequency content.

Technical scheme

Seismic exploration signal is an image of underground structure, contains abundant geophysical information.But, due to the earth The viscosity of medium, seismic exploration wavelet is in underground propagation, it may occur that energy attenuation and dispersion phenomenon, particularly high-frequency information Energy attenuation is serious.So as to cause seismic signal resolution ratio to decline so that the less earth formation letter related to high-frequency information Breath can not be reflected.Such case is unfavorable for that geologist carries out seismic signal explanation, and further to underground structure Make fine description.

The attenuation by absorption of the earth functions as a wave filter.This wave filter has by force to the radio-frequency component of source wavelet Attenuation by absorption effect so that so that radio-frequency component carry underground structure information can not show in the seismic data.But It is that these information do not disappear, is simply covered by high-energy component or noise, it is impossible to dominant to show.It is appropriate by finding Method of Mathematical Physics the high-frequency energy information of attenuation by earth absorption is compensated, make the structural information entrained by radio-frequency component It is dominant to represent, the purpose of seismic exploration signal resolution is improved to reach.

Assuming that seismic wavelet w (t) is an impact signal in short-term, contain abundant frequency content.Stratum reflectance factor r T (), is a white spectra signal for changing greatly.If without earth filtering effect, the preferable seismic exploration letter of wave detector pickup Number it is

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

In formula (1), symbol " * " is convolution operation symbol.Due to being acted on by attenuation by earth absorption, the ground for actually obtaining Shaking signal 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.Done at formula (2) two ends Fourier is converted, and obtains its Spectrum Relationship

Amplitude spectrum is taken at formula (3) two ends, and is reduced to

In formula (4),It is the Fourier transforming amplitudes spectrum of s (t),It is s₀The Fourier transforming amplitudes of (t) Spectrum, At (f) is a synthesis result of noise and the earth decay, is called amplitude decay function.At (f) can be cutd open by actual seismic Face mean amplitude spectrum is obtained.Thus amplitude decay function constructs compensating factor, and actual seismic section is compensated.

First, compensation causes that actual seismic section amplitude spectrum returns to perfect condition

Secondly, earth formation is complicated and changeable, therefore the attenuation by absorption rule of seismic signal is complicated.Due to earthquake The decay of signal is changed over time, thus, restoration and compensation process should have time-varying feature；Again due on synchronization stratum pair Different frequency content attenuation by absorption is different, thus restoration and compensation should also be with frequency qualitative change.Therefore, this hair The bright attenuation compensation method that is given is realized in the time and frequency zone of seismic signal, when being converted by the J to seismic exploration signal-frequency Spectrum does appropriate regulation, plays a part of attenuation compensation, so as to realize the target of the resolution ratio of raising seismic signal.

A kind of compensation J conversion it is multiple when-the frequency spectrum method that improves seismic profile resolution ratio, it is characterised in that step is as follows：

Step 1：FFT is done to each seismic channel in seismic profile and obtains amplitude spectrum, ask for the mean amplitude spectrum of seismic profile, Mean amplitude spectrum frequency-decible chart is represented, the corresponding frequency F of maximum amplitude point is found out in frequency-decible chart_m, according to Amplitude spectrum variation tendency will be greater than F_mFrequency range part be divided into three sections：First paragraph, adjacent F_mFrequency, its upper and lower ripple of amplitude spectral amplitude ratio It is dynamic, trend is not decreased obviously, frequency range is referred to as shaken, this frequency range is needed not compensate for；Second segment, its amplitude spectral amplitude ratio has one soon Speed declines, and referred to as decay frequency range, and this frequency range needs compensation；3rd section, its amplitude spectral amplitude ratio downward trend is gentle, and amplitude is very Small, referred to as high noisy frequency range, the Signal-to-Noise of this frequency range is relatively low, when should not carry out-frequency spectrum compensation；

Step 2：Decay frequency range is done by linear fit or done sectional linear fitting by different downward trends and is obtained amplitude spectrum and is declined Subtraction function At (f)=af+b, wherein a and b are the parameter of fitting, during construction-the frequency spectrum compensation factor：

Wherein, f₁It is the beginning frequency point of the frequency range that decays, f_nIt is the ending frequency point of the frequency range that decays；

Step 3：By the seismic channel contained by seismic profile by road carry out J conversion obtain multiple when-spectral matrix：

Described J converts expression formula：

Wherein, Tr () be seismic channel, Tr (k Δ T) for seismic channel k-th time sampling point, common N+1 sampled point, Δ T is time sampling interval；Δ F is frequency sampling interval, when n is represented-temporal sampled point sequence number, when m is represented-frequency spectrum frequently Rate sampled point sequence number, always points are M+1 to frequency sampling, and σ is resolution factor；

Step 4：When calculating multiple-mould of spectral matrix each element：

The maximum modulus value of every row is found out, maximum modulus value vector (M (0) M (1) ... M (N)) is obtained^T, during construction-frequency spectrum benefit Repay function：

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

Step 5：When used time-frequency spectrum compensation function pair is multiple-spectral matrix compensates by row, column：

Step 6：Using inverse J conversion to compensation after it is multiple when-frequency spectrum reconfiguration seismic channel；

Described inverse J converts expression formula：

Step 7：With trace equalization standardized method to reconstruct seismic channel amplitude carry out normalization correction, make its order of magnitude with Seismic channel before conversion is consistent：Sampling point amplitude is equal to maximum S in finding out seismic traces_MSampling point number, be designated as K, The preceding K big value to reconstructing seismic channel directly replaces with S_M, the big value of its k-th is designated as T_K, other data standardize as the following formula：

Described ε spans are 0.001 >=ε ＞ 0.

Described p spans are 1 >=p ＞ 0.

Described σ spans are more than 0.4.

Beneficial effect

It is proposed by the present invention it is a kind of by compensate that J converts it is multiple when-the frequency spectrum method that improves seismic profile resolution ratio, can be bright The aobvious resolution ratio for improving seismic exploration signal.There is fidelity and stability by actual seismic exploration data verification this method. Method And Principle understands that realize simply, user's operation is few.

Brief description of the drawings

Fig. 1 is original seismic profile, without by enhancing resolution processes；

Fig. 2 is original seismic profile amplitude spectrum and mean amplitude spectrum；

Fig. 3 is the frequency-decibel curve of original seismic profile mean amplitude spectrum；

Fig. 4 is that the frequency range of original seismic profile mean amplitude spectrum is divided.Decay frequency is determined according to this frequency-decibel curve figure Section, and fit attenuation function；

Fig. 5 is the seismic profile after attenuation compensation.Compared with the original seismic profile shown in Fig. 1, earthquake is cutd open after compensation The resolution ratio in face is significantly improved；

Fig. 6 is the amplitude spectrum of the seismic profile after attenuation compensation

Specific embodiment

In conjunction with embodiment, accompanying drawing, the invention will be further described：

The present invention provides following technical scheme, when converting multiple by J to seismic signal-frequency spectrum compensates, and plays right The effect that seismic signal energy attenuation is compensated, so as to realize improving the purpose of seismic resolution.It is small that method improves Gabor Ripple, its inverse transformation is constructed based on this improved Gabor wavelet.Due to the different small echo inversions with traditional sense of this inverse transformation Change, be a new conversion, the present invention is referred to as J conversion.Its implementation is as follows：

1) Gabor wavelet is unsatisfactory for small echo admissible condition, thus, there is no wavelet reconstruction formula.Gabor wavelet function is

The present invention is transformed formula (6).η=2 π is taken, and changes magnitude parameters and beσ is normal parameter, can be used for adjusting Save the time frequency window rectangular shape of Gabor wavelet, that is, when reconciling-frequency division resolution, the present invention referred to as resolution factor, general value More than 0.4.

Gabor mother wavelet functions are revised as

In formula (7), time parameter t is translated or change of scale can obtain Morlet wavelet functions race { g_{F, τ}(t)}

In formula (8), parameter f is the frequency window center frequency of the time frequency window of Gabor wavelet；τ is the time frequency window of Gabor wavelet When the window center time.If seismic signal is Tr (t), its Gabor wavelet is transformed to

In formula (9),It is g_{F, τ}The complex conjugate function of (t), operator " * " convolution operation.

2) signal reconstruction is realized, must there is inverse transformation, i.e. reconstruction formula.Allow bar because Gabor wavelet is unsatisfactory for small echo Part, is not to allow small echo, hence without the reconstruction formula of traditional small echo.The present invention is given and does not need admissible condition limit The inverse transformation of the Gabor wavelet of system

The wavelet inverse transformation of traditional sense, therefore the present invention are different from due to the Gabor wavelet inverse transformation that formula (10) is given The conversion that formula (9) and formula (10) are given to being called that J is converted, wherein, formula (9) is J direct transforms, and formula (10) is inverse J conversion.For side Just numerical operation, formula (10) can carry out conversion below

In formula (11), note

M (t, f)=JT (t, f) * e^j2πft (12)

When being carried out to seismic signal-frequency division solution and with compensation after when-frequency spectrum reconfiguration seismic signal when, need to be to formula (9), (11) and (12) discretization, to realize numerical operation.

The present invention also provides following technical scheme, and the high-frequency energy that method will greatly be absorbed to seismic signal is mended Repay, the purpose of seismic signal resolution ratio is improved to reach.Rational compensation have to be understood that the attenuation by absorption change rule of seismic signal Rule.Due to the complexity and the otherness of seismic exploration data acquisition and processing (DAP) mode of underground structure, the present invention uses statistics side Method obtains the technical scheme of the rule that seismic signal energy changes with frequency decay.Implementation method is as follows：

1) FFT is to each seismic channel in seismic profile, the amplitude spectrum of seismic channel is obtained；

2) to the amplitude spectrum of all seismic channels, add up, averagely, obtain a mean amplitude spectrum；

3) maximum amplitude of note mean amplitude spectrum is A_E, the corresponding frequency of the amplitude is designated as F_E, it is called maximum amplitude frequently Rate.On the basis of maximum amplitude frequency, the frequency-decible chart of mean amplitude spectrum is made；

4) with maximum amplitude frequency as branch, the frequency of amplitude spectrum is divided into two sections, is referred to as anterior frequency range and rear portion Frequency range.Anterior frequency range, i.e. frequency are less than F_EFrequency range, it is believed that comparatively energy does not decline the frequency information contained in it Subtract, its energy is needed not compensate for.Rear portion frequency range, i.e., more than F_EFrequency range, it is believed that the frequency information contained in it is come relatively Say that energy has greater attenuation, its energy demand has selection to compensate.

5) according to the Changing Pattern of amplitude spectrum, rear portion frequency range can substantially divide three sections.Can be divided into according to amplitude spectrum Changing Pattern Three types, by apart from maximum amplitude frequency F_EDistance be divided into three sections.First paragraph, adjacent F_E, assignment is shaken up and down in the frequency range Swing, do not show attenuation trend, this frequency range need not compensation, referred to as oscillation frequency bands；Second segment, amplitude spectrum significantly decreases Gesture, this frequency range needs compensation, and referred to as decay frequency range；3rd section, away from F_EFarthest, amplitude spectrum general trend is still downward trend, but Downward trend is gentle and amplitude spectral amplitude ratio is very low, and noise element is very high, and this frequency range is without compensation, referred to as high noisy frequency range.In reality In the problem of border, decay frequency range may have the segment of different downward trends, can subdivide.High noisy frequency range compensates it to reduce to be believed Number signal to noise ratio, seismic signal resolution ratio is reduced on the contrary, it is therefore not necessary to compensate, even can reject；

6) the decay frequency range of rear portion frequency range, is the frequency range for needing compensation.Linear fit can be done, or press according to actual conditions Different downward trends, sectional linear fitting obtains amplitude spectrum attenuation function.When then being constructed with amplitude spectrum attenuation function-frequency spectrum Compensating factor and penalty function.

The present invention also provides following technical scheme：By road with J convert, obtain seismic channel it is multiple when-frequency spectrum；When-frequency spectrum Constrained by two variables of time and frequency during regulation, i.e., compensation be both time-varying and frequency become；During multiple after to compensation- Frequency spectrum, reconstruct seismic signal is converted with inverse J.Comprise the following steps：

The first step, J conversion is done to current seismic road, obtain current seismic road it is multiple when-spectral matrix, then calculate multiple When-spectral matrix each element modulus value, its modular matrix is obtained, when referred to as-spectrogram, and the time point is searched out on each time point It is maximum when-frequency modulus value M (t)；

Second step, to it is multiple when-spectral matrix in each element, used time-frequency spectrum compensation function is adjusted, is compensated When afterwards again-spectral matrix；

3rd step, during by inverse J conversion using after compensation again-spectral matrix reconstruct seismic channel；

4th step, the seismic channel of reconstruct of standardizing, obtains high-resolution new seismic channel.

The present invention realizes improving seismic profile resolution ratio using following steps：

The first step, data analysis.

1) the maximum frequency that seismic signal can be showed is determined by the time-domain sampling rate and Nyquist sampling theorems of cross-sectional data Rate scope [0, F_max], wherein, F_maxIt is the highest frequency of signal；

2) for seismic profile, FFT is by road, and counts the mean amplitude spectrum for obtaining seismic channel in section, search for section The maximum of the amplitude spectrum of upper all road FFT and maximum mean amplitude of tide point (A_E, F_E), wherein, A_EIt is maximum average amplitude value, F_E It is the corresponding frequency of maximum average amplitude value；

3) frequency-decible chart of each channel amplitude spectrogram of drafting normalization, and mean amplitude spectrum；

4) with the corresponding frequency F of maximum mean amplitude of tide_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, information energy loss of the frequency that anterior frequency range is included corresponding to it is needed not compensate for. Rear portion frequency range can be divided into three parts, according to amplitude spectrum feature, by F_EDistance from closely to being far designated as respectively：Oscillating part, decay Part and high noisy part, wherein, the adjacent F of oscillating part_E, the corresponding amplitude spectrum of this component frequency vibrates up and down, without obvious Downward trend, its amplitude is often fluctuated up and down between 0~-8dB in -4dB horizontal lines, and the information of this component frequency need not be mended Repay；Attenuation portions, in the middle of the frequency range of rear portion, the amplitude spectrum of this band frequency shows the trend of rapid decrease, span from- 4dB or so drops to -35dB or so, and the information of this component frequency needs compensation；High noisy part, positioned at rear portion frequency range most end, Under -30dB, amplitude spectral amplitude ratio very little, this section of amplitude spectrum downward trend is gentle, and energy ezpenditure is maximum, can almost regard as amplitude Noise.The corresponding information of the band frequency, there is relatively low signal to noise ratio, if compensation can increase signal noise, reduces resolution ratio, this portion The information of crossover rate is also without compensation.

Second step, Amplitude spectrum attenuation function, constructs time-frequency spectrum compensating factor.

Although all frequency contents of seismic signal can all be subject to the attenuation by absorption on stratum to act in communication process, herein Only consider relative attenuation, therefore, with frequency F_ECorresponding amplitude A_EOn the basis of, it is believed that the only first step 4) described in rear portion frequently Section [F_E, F_max] in the information of " attenuation portions " and " high noisy part " this two band frequency there is relative energy to decay.

Energy compensating is only implemented to " attenuation portions ", therefore, this step is first to " attenuation portions " in frequency-decible chart Curve makees linear fit or sectional linear fitting.If curve point is

(f₁, A₁),(f₂, A₂),…,(f_n, A_n) (13)

Wherein, f_iIt is frequency, A_iBe its corresponding decibel amplitude, i=1,2 ..., n. then amplitude spectrum attenuation function it is linear Fitting function is

At (f)=af+b (14)

In formula (14), independent variable is the parameter that frequency f, a and b are fitting.When corresponding-the frequency spectrum compensation factor is

Because " attenuation portions " variation tendency may be inconsistent, sectional linear fitting, now, time-frequency can be done by variation tendency Spectrum compensating factor will be piecewise linear function.By formula (15) it can be seen that when-the frequency spectrum compensation factor be frequency function.To need not The frequency range of compensation, compensating factor is constant 1, that is, have

3rd step, when being carried out by road to seismic profile-frequency spectrum compensation.

Energy attenuation to seismic profile compensates the realization of Shi Zhu roads.First, when being made of J conversion to seismic channel-frequency Decompose, obtain its it is multiple when-frequency spectrum；During construction then ,-frequency spectrum compensation function, with penalty function adjust seismic channel it is multiple when-frequency Spectrum；Finally, seismic channel is reconstructed with inverse J transformation for mula.Implement step as follows：

0) common parameter is set for the entire profile treatment：Smoothing factor p, and 1 >=p ＞ 0；Regularization factors ε, and 0.001 >=ε ＞ 0；Obtain pending geological data maximum S_M, road sampling number (N+1) and time sampling interval Δ T. by Nyquist Sampling Theorems, calculate the largest interval Δ F of frequency sampling_max, and it is arranged as required to the frequency used during treatment Sampling interval Δ_F, and have Δ F≤Δ F_max, frequency sampling points are calculated, it is designated as (M+1)；

1) f=m Δ F, fixed frequency f are taken, by m_{F, 0}T () discretization in time domain, has

Using formula (17), discrete type (9), and noteHave

When can answer-spectral matrix

Wherein, J (n, m) is plural number (n=0,1 ..., N；M=0,1 ..., M）.When calculating multiple-modular matrix of spectral matrix, And when being called-spectrogram matrix, i.e.,

Temporally, i.e., the maximum modulus value at each time point is found out by row, there is 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 construction-frequency spectrum compensation function,

During with formula (23) to multiple shown in formula (20)-matrix that is given of frequency spectrum in each element interaction, can be compensated During afterwards multiple-spectral matrix

4th step, during with after compensation again-frequency spectrum reconfiguration seismic signal, to improve resolution ratio.

In order to realize restructing operation, need to be to formula (11) and formula (12) discretization.First to index tuning signal discretization.Take Frequency f=m Δ F, fixed frequency f, have

By formula (25) can discretization formula (12), have

Can discretization formula (11), the seismic signal rebuild by formula (26)Have

1) seismic signal is reconstructed.During to after compensation again-frequency spectrum, shown in formula (24), provide method by formula (26) and process, have

To formula (28), method is given with formula (27), it is cumulative by row, seismic channel can be reconstructed

2) width normalization of going bail for reconstruct seismic channel is processed.Find out seismic traces sample value and take maximum S_MSampled point Number, is designated as K, and to reconstruct seismic channel data sequence, the big value of its k-th is designated as T_K, the preceding K big value to reconstructing seismic channel is directly Replace with S_M, other data standardize as the following formula

To, per being operated by previous step respectively together, being capable of achieving resolution ratio enhancing effect on section.

Claims (1)

1. a kind of compensation J conversion it is multiple when-the frequency spectrum method that improves seismic profile resolution ratio, it is characterised in that step is as follows：

Step 1：FFT is done to each seismic channel in seismic profile and obtains amplitude spectrum, ask for the mean amplitude spectrum of seismic profile, will be flat Equal amplitude spectrum frequency-decible chart represented, the corresponding frequency F of maximum amplitude point is found out in frequency-decible chart_m, according to amplitude Spectrum variation tendency will be greater than F_mFrequency range part be divided into three sections：First paragraph, adjacent F_mFrequency, its amplitude spectral amplitude ratio fluctuates up and down, does not have Trend is decreased obviously, frequency range is referred to as shaken, this frequency range is needed not compensate for；Second segment, under its amplitude spectral amplitude ratio has one quickly Drop, referred to as decay frequency range, and this frequency range needs compensation；3rd section, its amplitude spectral amplitude ratio downward trend is gentle, and amplitude very little, claims It is high noisy frequency range, the Signal-to-Noise of this frequency range is relatively low, when should not carry out-frequency spectrum compensation；

Step 2：Linear fit is done to decay frequency range or sectional linear fitting is done by different downward trends and obtains amplitude spectrum decay letter Number At (f)=af+b, wherein a and b is the parameter of fitting, during construction-the frequency spectrum compensation factor：

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

CF (f)=1.0；

Wherein, f₁It is the beginning frequency point of the frequency range that decays, f_nIt is the ending frequency point of the frequency range that decays；

Step 3：By the seismic channel contained by seismic profile by road carry out J conversion obtain multiple when-spectral matrix：

Described J converts expression formula：

$\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|}{\mathrm{\σ}\sqrt{2\mathrm{\π}}}\exp \{-\frac{{\left(\mathrm{m\ΔF}\right)}^2{\left[(n-k)\mathrm{\ΔT}\right]}^2}{{2\mathrm{\σ}}^2}\}\exp \left[i2\mathrm{\πm\ΔF}(n-k)\mathrm{\ΔT}\right]\\ n=\mathrm{0,1},...,N\end{array}$

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

Wherein, Tr () is seismic channel, and Tr (k Δ T) is k-th time sampling point of seismic channel, common N+1 sampled point, and Δ T is Time sampling interval；Δ F is frequency sampling interval, when n is represented-temporal sampled point sequence number, when m is represented-spectral frequencies adopt Sample point number, always points are M+1 to frequency sampling, and σ is resolution factor；

Step 4：When calculating multiple-mould of spectral matrix each element：

The maximum modulus value of every row is found out, maximum modulus value vector (M (0) M (1) ... M (N)) is obtained^T, during construction-frequency spectrum compensation letter Number：

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

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

Step 5：When used time-frequency spectrum compensation function pair is multiple-spectral matrix compensates by row, column：

Step 6：Using inverse J conversion to compensation after it is multiple when-frequency spectrum reconfiguration seismic channel；

Described inverse J converts expression formula：

Step 7：Normalization correction is carried out to the amplitude for reconstructing seismic channel with trace equalization standardized method, makes its order of magnitude with conversion Preceding seismic channel is consistent：Sampling point amplitude is equal to maximum S in finding out seismic traces_MSampling point number, be designated as K, counterweight Preceding K big value of structure seismic channel directly replaces with S_M, the big value of its k-th is designated as T_K, other data standardize as the following formula：

Described ε spans are 0.001 >=ε ＞ 0；

Described p spans are 1 >=p ＞ 0；

Described σ spans are more than 0.4.