CN102478671A - Method for suppressing controllable seismic-source harmonic-wave interference - Google Patents

Abstract

The invention discloses a method for suppressing harmonic-wave interference existed in an original seismic record of controllable land seismic source excitation acquisition. Method is characterized by: taking a related previous earthquake data to carry out short-time Fourier transform and transforming a time domain signal to a time-frequency domain; taking a time sampling sequence number and a frequency sampling sequence number in a related previous time-frequency domain signal as two coordinate axes of a rectangular coordinate, carrying out coordinate transformation so that a fundamental wave is parallel to a time axis; separating the harmonic wave interference from the signal, carrying out coordinate inverse transformation of the time-frequency domain and carrying out short-time Fourier inverse transform so as to acquire time domain data which can remove the harmonic wave interference; after all unrelated seismic data collected from a filed is processed, acquiring the related seismic data after the harmonic wave interference is suppressed. By using the method of the invention, disadvantages of a filtering method in a frequency wave number domain can be avoided. Through the controllable seismic source, the seismic data with high qualities can be collected.

Description

A kind of method of suppressing the vibroseis harmonic interference

Technical field

The present invention relates to the oil seismic exploration data acquisition technology, is the method that excites the harmonic interference that exists in the original seismic data of collection to suppress to land vibroseis.

Background technology

The first step is to use the manual method earthquake-wave-exciting in the oil seismic exploration data collection task.Some basic requirements that excite to seismic event.But in order to adapt to the characteristics of various surface conditions and concrete method of work, used focus and exciting method are again diversified.To earthquake-wave-exciting enough strong energy to be arranged, because, in petroleum prospecting, find out the structural feature on a whole set of stratum of underground a few km depth rangees with seismic reflection, seismic event passes to underground reflected back again ground from focus.Propagate long like this distance, it is not all right not having enough energy.The signal that produces enough energy has two kinds of methods, and the one, the amplitude of increase signal, the 2nd, the duration of prolongation signal.Present land earthquake data acquisition mainly uses explosive and two kinds of epicenter excitations of vibroseis to the signal of underground propagation, for explosive, can adopt a series of measures to increase the amplitude of excitation signal, but can not prolong the duration of signal.For vibroseis, can prolong the duration of signal, but be difficult to increase amplitude.

Be at the near surface formation fired charge when explosive excites earthquake data acquisition, the very strong force signal that blast produces, it can be regarded as a pulse that amplitude is very big; Has very wide frequency spectrum; This pulse partly is attenuated absorption in the earth communication process medium-high frequency, becomes a band-limited wavelet, and it carries the geological information that lands down from the underground ground of returning; The wave detector that is laid in ground is noted, the seismologic record of Here it is explosive excites earthquake data acquisition.Facts have proved, in the oil gas field seismic prospecting, use explosive source, good effect is arranged, outstanding advantage is to have stronger energy.But explosive source also exposes a lot of shortcomings, and the most of energy loss that produces such as exploding is on the destruction to rock; Generation does not fit into the frequency content of earth-layer propagation; It is bigger with the required expense of use explosive to bore shot hole; Lack of water, top layer are that drilling well area in hardship constructions such as gravel, quicksand, marsh are inconvenient; The manufacturing district, the densely populated area should not use explosive; The most important thing is that the blasting explosive secure risk is high, the feature of environmental protection is poor.In order to overcome above-mentioned unfavorable factor, vibroseis becomes the main focus of earthquake data acquisition gradually.Vibroseis is compared with the explosive source of routine, and it has a lot of advantages.At first, the frequency band that vibroseis can be selected to be adapted to earth-layer propagation most according to formation characteristics is as scan band, thereby do not produce the vibration frequency that does not propagate on the stratum; Secondly, during vibroseis work, consumed energy does not destroy rock; At last, because it is the correlation analysis principle that adopts, so can avoid many interference, these outstanding advantages make it in earthquake data acquisition, be widely used.

Vibroseis excites the method for earthquake data acquisition different with explosive.Earthquake controllable earthquake focus system is made up of a sweep generator, a focus and correlator three parts.When exciting with vibroseis when carrying out earthquake data acquisition, sweep generator at first produces the electric signal of the length that uniform amplitude, a video frequency gradually changes, and is referred to as sweep signal, reference scan, reference signal or guide's scanning.Focus changes into force signal to sweep signal, and the vibration dull and stereotyped through focus is transmitted into force signal underground.The ground vibration that wave detector records when carrying out earthquake data acquisition through vibroseis is not the convolution result of a band-limited wavelet and underground reflection coefficient sequence; It is not a seismologic record truly; Also need be through relevant way the band-limited wavelet that becomes of propagating sweep signal through the earth, this process is accomplished through correlator.Correlator receives geological data (not relative recording or relevant preceding record) and reference signal simple crosscorrelation to wave detector; To shorten short finite bandwidth wavelet and reflection coefficient sequence convolution result (relative recording) from the long scan signal pressure of underground reflection into, the peak amplitude of signal has strengthened in this process.Three ingredients of earthquake controllable earthquake focus system have constituted three links that vibroseis carries out earthquake data acquisition, accomplish when wherein the 3rd link can be gathered in the open air, also can accomplish in the indoor Data Processing stage.

In three links of vibroseis earthquake data acquisition, whichsoever link all be unable to do without sweep signal, and it is the wavelet to underground propagation that vibroseis produces, so its characteristic will influence the quality of single shot record to a great extent.Mainly confirm to sweep sweep signal through designing several important parameters when vibroseis is gathered, they are sweep signal length L _s, the start-stop frequency f _sAnd f _e, frequency function f (t), the F that exerts oneself of focus.Sweep signal can be expressed as s (t)=Fsin (f (t) * t), and the t in the formula represents the time, and f (t) has determined the type of sweep signal.When f (t) was a constant, s (t) was a sinusoidal signal; As f (t)=f ₁During ± at, s (t) is the linear sweep signal, a=(f _e-f _s)/L is a sweep speed, f ₁Initial frequency for linear sweep; When f (t) was the caret form, s (t) was the nonlinear sweep signal.Linear sweep is to use maximum scan modes in producing, and it divides raising frequency and frequency reducing dual mode, as f (t)=f ₁Be called linear upsweep during+at, as f (t)=f ₁Be called linear down sweep during-at.The energy size that produces when knowing that by sweep signal vibroseis is gathered this shows for : can be through increasing the purpose that sweep time or enhancing signal intensity reach increases seismic source energy.It is suitable that but the size of exerting oneself must be selected, and it can not be too big, can not be too little.If exert oneself above the restriction of focus maximum static lotus ballast, will constitute a threat to the data quality.Because in this case, the vibroseis hydraulic system will be operated in nonlinear area, thereby cause sweep signal to distort, and then influence the quality of data.Simultaneously, the selection of time can not be oversize, because this will have influence on field construction efficient.

There is special interference in vibroseis with the machine-processed different records that make that also it produces of explosive source.The power valve core that reference signal promotes in the electrohydraulic servo valve during work of vibroseis is done the sine sweep campaign; The power valve is input to the focus Vib. to the high pressure liquid from hydraulic system according to the Changing Pattern of reference signal, thereby produces huge mechanical vibration force and the flat board through close proximity to ground acts on the earth.Import in the underground process in reference signal; Because making with resonance, the delay vibration between the earth vibration and the earthquake controllable earthquake focus system imports the signal that underground reference signal produces nonlinear distortion into; It is exactly the not relevant geological data that earthquake-capturing obtains that reference signal and distorted signal are propagated through being laid in behind the earth that ground wave detector notes; Or being called relevant preceding geological data, relevant preceding seismologic record, it is exactly original earthquake data that these data and sweep signal are carried out the relevant relative recording that obtains.Usually call harmonic wave, harmonic interference or resonance to the distorted signal in the geological data before relevant, call first-harmonic to the signal that does not have distortion.Theoretical according to correlation analysis; Geological data is with after sweep signal is relevant before relevant; First-harmonic is relevant with sweep signal, and what produce is the needed useful signal of earthquake data acquisition, and the harmonic interference of the relevant generation with sweep signal of harmonic wave is at " the distortion tail " that show as on the relative recording after appearing at the correlator wavefront.When " distortion tail " was constantly identical with the appearance of relevant wavelet, both mutual superposition made relevant wavelet distortion, and then influenced the continuity of lineups.If the energy level of " distortion tail " is greater than weak reflection correlator wave energy level, lineups just might be submerged in the harmonic interference so.Harmonic interference shows as " pinniform " or " herring-bone form " in big gun collection record.These harmonic interference have a strong impact on the quality of section, and then influence seismic resolution.

Harmonic interference is inevitable phenomenon in the vibroseis work progress, can only hope the generation harmonic distortion of minimum degree in the open air when we construct.The kind of harmonic interference is very complicated, in real work, usually it is divided into four big types: higher hamonic wave, and promptly the integral multiple of fundamental frequency generally is present in the frequency range of the whole acting force of focus; Low-order harmonic, particularly 1/2 subharmonic as fundamental frequency ratio component, appear in the limited frequency range of source signal usually; Ultra subharmonic is the higher harmonic wave of low-order harmonic composition medium frequency, generally in real data, occurs seldom; Irregular or vibration at random in the focal force signal, this hunting may influence the relative recording of vibroseis in whole focus force signals frequency range.

The kind of harmonic wave is complicated, and the reason that harmonic wave produces also is diversified.Can the reason that harmonic interference produces be divided into two types simply: the coupling reason on vibroseis mechanical reason and focus and the face of land.The harmonic interference that the failure factor of vibroseis self causes can be got rid of through maintenance, and the harmonic interference that the coupling factor on the vibroseis and the face of land causes is inevitable phenomenon in the vibroseis work progress.The method of eliminating the vibroseis harmonic interference probably is divided into two big types at present: the first kind is to design rational vibroseis field construction parameter; Second type is the property difference according to harmonic interference and significant wave, adopts the method for digital filtering to remove.Before a kind of method paid attention to by people because it directly has influence on field construction efficient.And because field construction parameter reasonable in design is to start with from the root that produces harmonic interference, so a kind of construction method of science can access more correct field data.Yet from producing the root of harmonic interference, no matter which kind of field acquisition parameter also will produce certain harmonic interference.This just makes a kind of effective digital processing method of research necessitate.People (1992) such as David A.Okaya when the earth's crust data of San Joaquin Valley of research California with vibroseis not related data to transform to temporal frequency domain discovery harmonic interference be linearity with significant wave and have Different Slope.According to this characteristic, at first be employed in the method that temporal frequency excises and suppress harmonic noise.This method can be removed harmonic noise to a certain extent, but has kept harmonic energy greatly.Further carry out the wavenumber domain of two-dimension fourier transform to the time-frequency space, i.e. K to the data of temporal frequency domain _f-K _TThe territory finds that significant wave harmonic feature of noise is more obvious.At K _f-K _TFiltering and inverse transformation are carried out after time domain to harmonic wave in the territory, and harmonic energy has been eliminated basically.This method has certain practicality, but also exists not enough: because at K _f-K _TWhat harmonic interference and significant wave leaned in the territory is very near, so be difficult to confirm a more rational excision extension; And resonance and useful signal are superimposed near initial point, and therefore this method also is difficult to complete harmonic carcellation and disturbs.

Summary of the invention

The object of the invention is to overcome directly at K _f-K _TThe shortcoming that cutting method exists in the territory causes harmonic interference to the coupling factor on the vibroseis and the face of land, and a kind of method of suppressing the vibroseis harmonic interference is provided.

The present invention provides following technical step:

1) field acquisition geological data is got relevant together preceding geological data and is carried out Short Time Fourier Transform, transforms to time-frequency domain to time-domain signal;

The described Short Time Fourier Transform of step 1) is:

$X[n,m]=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[k\right]w[k-n]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(1\right)$

In the formula: n and k are the time-sampling sequence number, n=1, and 2,3 ..., N, k=1,2,3 ..., N;

M is the frequency sampling sequence number, m=1, and 2,3 ..., M;

L is the length of geological data before relevant;

Δ t is the time sampling interval;

Δ f is the frequency sampling interval;

N is the number of samples of time-sampling, N=L/ Δ t;

M is the number of samples of frequency sampling, M=1/2 Δ t Δ f;

W [k-n]=w (k Δ t-n Δ t) is the window function of symmetry, the expression formula of window function when window function adopts Blacknam:

$w\left(t\right)=\left\{\begin{array}{c}0.42+0.5\cos (\mathrm{\&pi;t}/T)+0.08\cos (2\mathrm{\&pi;t}/T),t\&Element;[-T,T]\\ 0,t\&NotElement;(-T,T)\end{array}\right.---\left(2\right)$

In the formula: t=k Δ t-n Δ t is the time, and T is the length of half window function;

X [n]=x (n Δ t) is relevant together arbitrarily preceding discrete geological data;

X [n, m]=X (n Δ t, m Δ f) is the Short Time Fourier Transform of discrete geological data before relevant:

The length of described window function is 0.8 second to 1.2 seconds.

2) will be correlated with before n and m among the time-frequency domain signal X [n, m] as two coordinate axis of rectangular coordinate system, carry out coordinate transform, make first-harmonic parallel: X with time shaft _z[n _z, m _z]=X [n, m];

Step 2) said coordinate transform is at first calculated rotation angle by following (3) formula, then by following (4) and (5) formula to all coordinate transforms among the X [n, m] preceding point (n, m) point (n behind the calculating coordinate change _z, m _z), just accomplished X [n, m] to X _z[n _z, m _z] coordinate transform;

The coordinate transform anglec of rotation is:

${\mathrm{\&theta;}}_0=\arctan \left(\frac{L_s}{(f_e-f_s)(2T+1)\mathrm{\&Delta;t\&Delta;t}}\right)---\left(3\right)$

In the formula: θ ₀Be the coordinate system rotation angle; L _sBe sweep signal length; f _sAnd f _eBe respectively the initial frequency and termination frequency of sweep signal;

Rectangular coordinates transformation becomes polar coordinates to be:

$\left\{\begin{array}{c}r=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">m</figure-callout>}}^2}\\ \mathrm{\&theta;}=\arctan (m/n)\end{array}\right.---\left(4\right)$

The polar coordinate transform coordinate that meets at right angles is:

$\left\{\begin{array}{c}{n}_z=r\cos (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\\ m_z=r\sin (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\end{array}\right.---\left(5\right)$

In the formula: r is a pole axis length; θ is a polar angle.

3) harmonic interference is separated with signal;

Described harmonic interference of step 3) and Signal Separation adopt Ka Hunan-Laue husband conversion.

Described harmonic interference of step 3) and Signal Separation are accomplished as follows:

1. utilize the Householder conversion with X _zBecome the triple diagonal matrix A of symmetry, computing formula is:

A＝U ^TX _zV (6)

In the formula: X _zBe time-frequency domain geological data not relevant after the coordinate transform, U ^TBe the transposition of U, U=u ₁U ₂K u _k, k=min{N, M-1}, V=v ₁V ₂K v _l, l=min{M, N-2}, each conversion u of U _j(j=1,2K k) is with X _zIn element below the j row principal diagonal become 0; Each conversion v of V _i(i=1,2K l) is with X _zIn i become 0 with the element on right minor diagonal element the right of principal diagonal next-door neighbour in capable, calculate following symmetric triple-diagonal matrix A through (6) formula:

$A=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1& {\mathrm{\&beta;}}_1& & \\ {\mathrm{\&beta;}}_1& {\mathrm{\&alpha;}}_2& {\mathrm{\&beta;}}_1& \\ O& O& O& \\ & O& {\mathrm{\&alpha;}}_{N-1}& {\mathrm{\&beta;}}_{N-1}\\ & & {\mathrm{\&beta;}}_{N-1}& {\mathrm{\&alpha;}}_N\end{array}\right],{\mathrm{\&beta;}}_i>0,i=\mathrm{1,2},\mathrm{\&Lambda;},M-1---\left(7\right)$

2. the QL transform method of the implicit displacement of utilization becomes diagonal matrix with the tridiagonal matrix A of symmetry, carries out according to the following steps:

(I) ask the afterbody second order submatrix of A Proper vector be p ₁And p ₂, get ε=p ₁

(II) make r _M=α _M, C _M=S _M=1, G _M=(α _M-ε) ², η _M=α _M-ε, to i=M, M-1, Λ, 2 by formula (8) calculating

With

With

$\left\{\begin{array}{c}{B}_i=C_i{\mathrm{\&beta;}}_{i-1}^2,F_i=S_i{\mathrm{\&beta;}}_{i-1}^2,\\ {\widehat{\mathrm{\&beta;}}}_i^2=F_i+G_i,\\ S_{i-1}=F_i/{\widehat{\mathrm{\&beta;}}}_i^2,C_{i-1}=1-S_{i-1},\\ {\mathrm{\&tau;}}_{i-1}={\mathrm{\&eta;}}_i{\mathrm{\&beta;}}_{i-1}^2/{\widehat{\mathrm{\&beta;}}}_i^2,\\ {\mathrm{\&rho;}}_{i-1}={\mathrm{\&alpha;}}_{i-1}-r_i,{\mathrm{\&kappa;}}_{i-1}=C_{i-1}-S_{i-1},{\mathrm{\&xi;}}_{i-1}={\mathrm{\&rho;}}_{i-1}{S}_{i-1}\\ {\mathrm{\&eta;}}_{i-1}=C_{i-1}{\mathrm{\&xi;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{\mathrm{\&tau;}}_{i-1},\\ G_{i-1}={\mathrm{\&rho;}}_{i-1}{\mathrm{\&eta;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}({\mathrm{\&tau;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{B}_i),\\ d_{i-1}={\mathrm{\&xi;}}_{i-1}+2{\mathrm{\&tau;}}_{i-1},\\ {\widehat{\mathrm{\&alpha;}}}_i=r_i+d_{i-1}\\ r_{i-1}={\mathrm{\&alpha;}}_{i-1}-d_{i-1}\end{array}\right.---\left(8\right)$

${\widehat{\mathrm{\&alpha;}}}_1=r_1,{\widehat{\mathrm{\&beta;}}}_1^2=G_1,$

(III) to i=M; M-1; Λ; 1, make up new triple diagonal matrix

with

and

$\hat{A}=\left[\begin{array}{cccc}{\widehat{\mathrm{\&alpha;}}}_1& {\widehat{\mathrm{\&beta;}}}_1& & \\ {\widehat{\mathrm{\&beta;}}}_1& {\widehat{\mathrm{\&alpha;}}}_2& {\widehat{\mathrm{\&beta;}}}_1& \\ O& O& O& \\ & O& {\widehat{\mathrm{\&alpha;}}}_{M-1}& {\widehat{\mathrm{\&beta;}}}_{M-1}\\ & & {\widehat{\mathrm{\&beta;}}}_{M-1}& {\widehat{\mathrm{\&alpha;}}}_M\end{array}\right]---\left(9\right)$

(IV) if

then

be exactly desired diagonal matrix, otherwise

press the method calculating of I to III again;

3. calculate eigenwert and corresponding left eigenvector and the right proper vector of eigenwert of tridiagonal matrix A;

Described tridiagonal matrix

Diagonal entry be exactly the eigenvalue of tridiagonal matrix A _i,

λ _iTo deserved left eigenvector μ _iAccording to matrix Ask for λ _iTo deserved right proper vector v _iAccording to matrix

Ask for, wherein

Be X _zTransposition.

4. use the singular value of computes tridiagonal matrix A:

is to the descending ordering of singular value;

5. harmonic interference and Signal Separation;

Said harmonic interference and Signal Separation or stick signal perhaps keep harmonic interference;

If stick signal then keeps a preceding S singular value and characteristic of correspondence vector, be 0 with remaining M-S singular value and the tax of characteristic of correspondence vector thereof then.

Then a preceding S singular value and characteristic of correspondence vector are composed 0 if keep harmonic interference, then remaining M-S singular value and characteristic of correspondence vector thereof are kept.

Said harmonic interference and Signal Separation singular value number s are 1 or 2 or 3.

6. with the singular value and the proper vector reconstruct data thereof that remain.

The corresponding data of each singular value of described reconstruct are pressed following formula:

$I_i={\mathrm{\&sigma;}}_i{\mathrm{\&mu;}}_i{v}_i^T,(i=\mathrm{2,3},\mathrm{\&Lambda;},M)---\left(10\right)$

X _z＝I ₁+I ₂+Λ+I _M (11)

(11) X in the formula is exactly the data of removing the data of harmonic interference after the reconstruct or keeping harmonic interference.

4) the coordinate inverse transformation of time-frequency domain is with (the n that has a few after the coordinate transform _z, m _z) get back to former coordinate in inverse transformation, the anglec of rotation of coordinate inverse transformation is the same with direct transform;

The said coordinate inverse transformation of step 4) is pressed following formula to X _z[n _z, m _z] middle (n that has a few _z, m _z) (n m), has just accomplished X to calculation level _z[n _z, m _z] to the coordinate inverse transformation of X [n, m];

$\left\{\begin{array}{c}{r}_z=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n\; z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}_z^{}+m_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}\\ {\mathrm{\&theta;}}_z=\arctan (m_z/n_z)\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}n=r_z\cos ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\\ m=r_z\sin ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\end{array}\right.---\left(13\right)$

(12) in the formula: r _zBe pole axis length; θ _zBe polar angle, (13) formula is the method that transforms polar coordinate system rectangular coordinate system.

5) do Fourier inversion in short-term;

The described Fourier inversion in short-term of step 5) is accomplished by following formula:

$x^{\&prime;}\left[k\right]=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}X[n,m]w[k-m]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(14\right)$

In the formula: the implication of parameter is with (1) formula.

6) obtain removing the time domain data of harmonic interference;

The time domain data of the said removal harmonic interference of step 6) is accomplished by following (15) formula or (16) formula:

If what step 4) kept is the data of removing harmonic interference, then:

x″[k]＝x′[k] (15)

If what step 4) kept is the data of harmonic interference, then:

x″[k]＝x[k]-x′[k] (16)

(15) in formula and (16) formula: " [k] is exactly the time domain geological data of removing harmonic interference to x.

7) get the relevant preceding geological data of other one field acquisition; Repeating step 2) to the process of step 5); Until handling the relevant geological data of all field acquisitions, obtain the downtrodden geological data of harmonic interference, obtain the relevant geological data after harmonic interference is suppressed after being correlated with.

The present invention has avoided the deficiency of filtering method in the frequency-wavenumber domain, has reached the purpose of gathering high-quality geological data through vibroseis better.

Description of drawings

Fig. 1 is that the present invention invents process flow diagram;

Fig. 2 is a relevant preceding wherein big gun geological data of field acquisition;

Fig. 3 shakes the 97th road earthquake data that extract the data relatively from the former beginning and end of Fig. 2;

Fig. 4 is that the single track signal carries out Short Time Fourier Transform and obtains time-frequency figure;

Fig. 5 is the time-frequency figure after Fig. 4 time-frequency domain data are carried out coordinate transform;

Fig. 6 is the singular value distribution plan after Ka Hunan-Laue husband conversion;

Fig. 7 is the time-frequency domain fundamental signal that after Ka Hunan-Laue husband conversion, obtains;

Fig. 8 is the time-frequency domain fundamental signal after the coordinate inverse transformation;

Fig. 9 is the time-domain signal that extracts after fundamental signal is done the compacting harmonic interference that obtains after the Fourier inversion in short-term;

Figure 10 is the time-frequency domain harmonic interference that after Ka Hunan-Laue husband conversion, obtains;

Figure 11 is the time-frequency domain harmonic interference after the coordinate inverse transformation;

Figure 12 extracts harmonic interference, the time-domain signal after time domain is subtracted each other the compacting harmonic interference that obtains then at time-frequency domain;

Figure 13 does not use the record that the present invention suppresses harmonic interference;

Figure 14 is the results of same data after the present invention suppresses harmonic interference.

Embodiment

The present invention can extract fundamental signal from time-frequency domain, can also extract harmonic interference at time-frequency domain, deducts harmonic interference in time domain with original signal then and obtains useful signal.The compacting flow process is as shown in Figure 1, and the rectangle frame among the figure has shown the residing state of data, and oval frame is represented processing that data are done, and arrow is indicated the flow direction of data.

Below in conjunction with actual geological data the specific embodiment of the invention is described respectively.

The enforcement that time-frequency domain extracts fundamental signal divides 6 steps.

1) field acquisition geological data is got relevant together preceding geological data and is carried out Short Time Fourier Transform, transforms to time-frequency domain to time-domain signal;

The described Short Time Fourier Transform of step 1) is:

$X[n,m]=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[k\right]w[k-n]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(1\right)$

In the formula: n and k are the time-sampling sequence number, n=1, and 2,3 ..., N, k=1,2,3 ..., N;

M is the frequency sampling sequence number, m=1, and 2,3 ..., M;

L is the length of geological data before relevant;

Δ t is the time sampling interval;

Δ f is the frequency sampling interval;

N is the number of samples of time-sampling, N=L/ Δ t;

M is the number of samples of frequency sampling, M=1/2 Δ t Δ f;

W [k-n]=w (k Δ t-n Δ t) is the window function of symmetry, the expression formula of window function when window function adopts Blacknam:

$w\left(t\right)=\left\{\begin{array}{c}0.42+0.5\cos (\mathrm{\&pi;t}/T)+0.08\cos (2\mathrm{\&pi;t}/T),t\&Element;[-T,T]\\ 0,t\&NotElement;(-T,T)\end{array}\right.---\left(2\right)$

In the formula: t=k Δ t-n Δ t is the time, and T is the length of half window function;

X [n]=x (n Δ t) is relevant together arbitrarily preceding discrete geological data;

X [n, m]=X (n Δ t, m Δ f) is the Short Time Fourier Transform of discrete geological data before relevant:

The length of described window function is 0.8 second to 1.2 seconds.

Fig. 2 is a relevant preceding wherein big gun geological data of field acquisition, the road number of the relevant preceding geological data of horizontal ordinate numeral among Fig. 2, and ordinate is a sampling sequence number, Fig. 2 is the display result of 2 milliseconds of samplings, therefore maximum sampling sequence number 11000.From this big gun geological data, extract together, such as the 97th road, Fig. 3 is exactly the 97th road earthquake data that extract, the amplitude of horizontal ordinate Digital Telemetry Seismic Wave among Fig. 3, and ordinate is a sampling sequence number.97 road earthquake data are carried out Short Time Fourier Transform, time-domain signal is transformed to time-frequency domain, transformation results is as shown in Figure 4; Horizontal ordinate express time among Fig. 4, ordinate are the expression frequencies, and the signal of many group Different Slope is arranged among the figure; Having only a kind of slope signal is first-harmonic, is useful signal.

2) will be correlated with before n and m among the time-frequency domain signal X [n, m] as two coordinate axis of rectangular coordinate system, carry out coordinate transform, make first-harmonic parallel: X with time shaft _z[n _zm _z]=X [n, m];

Step 2) said coordinate transform is at first calculated rotation angle by following (3) formula, then by following (4) and (5) formula to all coordinate transforms among the X [n, m] preceding point (n, m) point (n behind the calculating coordinate change _z, m _z), just accomplished X [n, m] to X _z[n _z, m _z] coordinate transform;

The coordinate transform anglec of rotation is:

${\mathrm{\&theta;}}_0=\arctan \left(\frac{L_s}{(f_e-f_s)(2T+1)\mathrm{\&Delta;t\&Delta;t}}\right)---\left(3\right)$

In the formula: θ ₀Be the coordinate system rotation angle; L _sBe sweep signal length; f _sAnd f _eBe respectively the initial frequency and termination frequency of sweep signal;

Rectangular coordinates transformation becomes polar coordinates to be:

$\left\{\begin{array}{c}r=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">m</figure-callout>}}^2}\\ \mathrm{\&theta;}=\arctan (m/n)\end{array}\right.---\left(4\right)$

The polar coordinate transform coordinate that meets at right angles is:

$\left\{\begin{array}{c}{n}_z=r\cos (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\\ m_z=r\sin (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\end{array}\right.---\left(5\right)$

In the formula: r is a pole axis length; θ is a polar angle.

Not relevant seismologic record is made up of first-harmonic harmonic two parts; Show as the clinographic curve of Different Slope at time-frequency plane; The purpose of coordinate transform is to do a rotation to these clinographic curves, in postrotational coordinate system, makes first-harmonic parallel with time shaft, and harmonic wave still is a clinographic curve.Fig. 5 is the result who time-frequency domain data shown in Figure 4 is carried out conversion according to coordinate transformation method, the time-frequency figure after the coordinate transform.The horizontal ordinate of data after the coordinate transform does not have the pure time and the notion of frequency, so here we represent the transverse axis and the longitudinal axis to count.Because coordinate transform is the rotation of doing to first-harmonic, the horizontal direction signal is exactly a first-harmonic after the rotation, and the harmonic interference that is of slope is arranged.

3) harmonic interference is separated with signal;

Described harmonic interference of step 3) and Signal Separation adopt Ka Hunan-Laue husband conversion.Ka Hunan-Laue husband conversion is also referred to as the He Tailing conversion, proper vector conversion or svd etc.Because first-harmonic is very strong in time-frequency domain internal linear correlativity, become next its linear dependence that further increases of level through carrying out it coordinate transform, its singular value just is far longer than non-first-harmonic composition on numerical values recited so.

Described harmonic interference of step 3) and Signal Separation are accomplished as follows:

1. utilize the Householder conversion with X _zBecome the triple diagonal matrix A of symmetry, computing formula is:

A＝U ^TX _zV (6)

In the formula: X _zBe time-frequency domain geological data not relevant after the coordinate transform, U ^TBe the transposition of U, U=u ₁U ₂K u _k, k=min{N, M-1}, V=v ₁V ₂K v _l, l=min{M, N-2}, each conversion u of U _j(j=1,2K k) is with X _zIn element below the j row principal diagonal become 0; Each conversion v of V _i(i=1,2K l) is with X _zIn i become 0 with the element on right minor diagonal element the right of principal diagonal next-door neighbour in capable, calculate following symmetric triple-diagonal matrix A through (6) formula:

$A=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1& {\mathrm{\&beta;}}_1& & \\ {\mathrm{\&beta;}}_1& {\mathrm{\&alpha;}}_2& {\mathrm{\&beta;}}_1& \\ O& O& O& \\ & O& {\mathrm{\&alpha;}}_{N-1}& {\mathrm{\&beta;}}_{N-1}\\ & & {\mathrm{\&beta;}}_{N-1}& {\mathrm{\&alpha;}}_N\end{array}\right],{\mathrm{\&beta;}}_i>0,i=\mathrm{1,2},\mathrm{\&Lambda;},M-1---\left(7\right)$

2. the QL transform method of the implicit displacement of utilization becomes diagonal matrix with the tridiagonal matrix A of symmetry, carries out according to the following steps:

(I) ask the afterbody second order submatrix of A

Proper vector be p ₁And p ₂, get ε=p ₁

(II) make r _M=α _M, C _M=S _M=1, G _M=(α _M-ε) ², η _M=α _M-ε, to i=M, M-1, Λ, 2 by formula (8) calculating

With

With

$\left\{\begin{array}{c}{B}_i=C_i{\mathrm{\&beta;}}_{i-1}^2,F_i=S_i{\mathrm{\&beta;}}_{i-1}^2,\\ {\widehat{\mathrm{\&beta;}}}_i^2=F_i+G_i,\\ S_{i-1}=F_i/{\widehat{\mathrm{\&beta;}}}_{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,C_{i-1}=1-S_{i-1},\\ {\mathrm{\&tau;}}_{i-1}={\mathrm{\&eta;}}_i{\mathrm{\&beta;}}_{i-1}^2/{\widehat{\mathrm{\&beta;}}}_{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,\\ {\mathrm{\&rho;}}_{i-1}={\mathrm{\&alpha;}}_{i-1}-r_i,{\mathrm{\&kappa;}}_{i-1}=C_{i-1}-S_{i-1},{\mathrm{\&xi;}}_{i-1}={\mathrm{\&rho;}}_{i-1}{S}_{i-1}\\ {\mathrm{\&eta;}}_{i-1}=C_{i-1}{\mathrm{\&xi;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{\mathrm{\&tau;}}_{i-1},\\ G_{i-1}={\mathrm{\&rho;}}_{i-1}{\mathrm{\&eta;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}({\mathrm{\&tau;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{B}_i),\\ d_{i-1}={\mathrm{\&xi;}}_{i-1}+2{\mathrm{\&tau;}}_{i-1},\\ {\widehat{\mathrm{\&alpha;}}}_i=r_i+d_{i-1}\\ r_{i-1}={\mathrm{\&alpha;}}_{i-1}-d_{i-1}\end{array}\right.---\left(8\right)$

${\widehat{\mathrm{\&alpha;}}}_1=r_1,{\widehat{\mathrm{\&beta;}}}_1^2=G_1,$

(III) to i=M; M-1; Λ; 1, make up new triple diagonal matrix

with and

$\hat{A}=\left[\begin{array}{cccc}{\widehat{\mathrm{\&alpha;}}}_1& {\widehat{\mathrm{\&beta;}}}_1& & \\ {\widehat{\mathrm{\&beta;}}}_1& {\widehat{\mathrm{\&alpha;}}}_2& {\widehat{\mathrm{\&beta;}}}_1& \\ O& O& O& \\ & O& {\widehat{\mathrm{\&alpha;}}}_{M-1}& {\widehat{\mathrm{\&beta;}}}_{M-1}\\ & & {\widehat{\mathrm{\&beta;}}}_{M-1}& {\widehat{\mathrm{\&alpha;}}}_M\end{array}\right]---\left(9\right)$

(IV) if

then

be exactly desired diagonal matrix, otherwise press the method calculating of I to III again;

3. calculate eigenwert and corresponding left eigenvector and the right proper vector of eigenwert of tridiagonal matrix A;

Described tridiagonal matrix

Diagonal entry be exactly the eigenvalue of tridiagonal matrix A _i,

λ _iTo deserved left eigenvector μ _iAccording to matrix

Ask for λ _iTo deserved right proper vector v _iAccording to matrix

Ask for, wherein Be X _zTransposition.

4. use the singular value of computes tridiagonal matrix A:

is to the descending ordering of singular value;

Fig. 6 is the singular value distribution plan after Ka Hunan-Laue husband conversion, and horizontal ordinate is the size of singular value among the figure, and ordinate is the number of singular value.As can be seen from the figure, singular value is big more, and number is more little, and singular value is more little, and number is big more, and the part that singular value is big is a fundamental signal, and the singular value fraction is a harmonic interference.

5. separate harmonic interference and signal;

Said separation harmonic interference and signal stick signal keep a preceding S singular value and characteristic of correspondence vector, are 0 with remaining M-S singular value and the tax of characteristic of correspondence vector thereof then,

The singular value number s of said harmonic interference and Signal Separation is 1 or 2 or 3.

6. with the singular value and the proper vector reconstruct data thereof that remain.

The corresponding data of each singular value of described reconstruct are pressed following formula:

$I_i={\mathrm{\&sigma;}}_i{\mathrm{\&mu;}}_i{v}_i^T,(i=\mathrm{2,3},\mathrm{\&Lambda;},M)---\left(10\right)$

X _z＝I ₁+I ₂+Λ+I _M (11)

(11) X in the formula is exactly data of removing harmonic interference after the reconstruct

Fig. 7 is that the singular value number is the time-frequency domain fundamental signal that after Ka Hunan-Laue husband conversion, obtained in 1 o'clock.As can be seen from the figure, after Ka Hunan-Laue husband conversion, there is the harmonic interference of slope to be removed, only kept the fundamental signal of horizontal direction.

4) the coordinate inverse transformation of time-frequency domain is with (the n that has a few after the coordinate transform _z, m _z) get back to former coordinate in inverse transformation, the anglec of rotation of coordinate inverse transformation is the same with direct transform;

The said coordinate inverse transformation of step 4) is pressed following formula to X _z[n _z, m _z] middle (n that has a few _z, m _z) (n m), has just accomplished X to calculation level _z[n _z, m _z] to the coordinate inverse transformation of X [n, m];

$\left\{\begin{array}{c}{r}_z=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n\; z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}_z^{}+m_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}\\ {\mathrm{\&theta;}}_z=\arctan (m_z/n_z)\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}n=r_z\cos ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\\ m=r_z\sin ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\end{array}\right.---\left(13\right)$

(12) in the formula: r _zBe pole axis length; θ _zBe polar angle, (13) formula is the method that transforms polar coordinate system rectangular coordinate system.

Fig. 8 is the time-frequency domain fundamental signal after the coordinate inverse transformation, can find out after the coordinate inverse transformation that coordinate has in length and breadth had identical meaning with coordinate in length and breadth among Fig. 4, and have only a kind of fundamental signal of slope this moment among the figure, and harmonic interference is removed.

5) do Fourier inversion in short-term;

The described Fourier inversion in short-term of step 5) is accomplished by following formula:

$x^{\&prime;}\left[k\right]=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}X[n,m]w[k-m]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(14\right)$

In the formula: the implication of parameter is with (1) formula.

(14) x ' in the formula [k] is exactly the time domain data of removing after the harmonic interference, and Fig. 9 is exactly the time-domain signal of doing after the compacting harmonic interference that obtains after the Fourier inversion in short-term.

6) get the relevant preceding geological data of other one field acquisition; Repeating step 2) to the process of step 5); Until handling the not relevant geological data of all field acquisitions; Obtain the downtrodden geological data of harmonic interference, the geological data input correlator after the harmonic interference compacting is correlated with, just obtain the relevant geological data after harmonic interference is suppressed.

The present invention can also adopt time-frequency domain to extract the embodiment of harmonic interference, divides 7 steps to accomplish.

1) field acquisition geological data is got relevant together preceding geological data and is carried out Short Time Fourier Transform, transforms to time-frequency domain to time-domain signal;

The described Short Time Fourier Transform of step 1) is:

$X[n,m]=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[k\right]w[k-n]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(1\right)$

In the formula: n and k are the time-sampling sequence number, n=1, and 2,3 ..., N, k=1,2,3 ..., N;

M is the frequency sampling sequence number, m=1, and 2,3 ..., M;

L is the length of geological data before relevant;

Δ t is the time sampling interval;

Δ f is the frequency sampling interval;

N is the number of samples of time-sampling, N=L/ Δ t;

M is the number of samples of frequency sampling, M=1/2 Δ t Δ f;

W [k-n]=w (k Δ t-n Δ t) is the window function of symmetry, the expression formula of window function when window function adopts Blacknam:

$w\left(t\right)=\left\{\begin{array}{c}0.42+0.5\cos (\mathrm{\&pi;t}/T)+0.08\cos (2\mathrm{\&pi;t}/T),t\&Element;[-T,T]\\ 0,t\&NotElement;(-T,T)\end{array}\right.---\left(2\right)$

In the formula: t=k Δ t-n Δ t is the time, and T is the length of half window function;

X [n]=x (n Δ t) is relevant together arbitrarily preceding discrete geological data;

X [n, m]=X (n Δ t, m Δ f) is the Short Time Fourier Transform of discrete geological data before relevant:

The length of described window function is 0.8 second to 1.2 seconds.

Fig. 2 is a relevant preceding wherein big gun geological data of field acquisition, the road number of the relevant preceding geological data of horizontal ordinate numeral among Fig. 2, and ordinate is a sampling sequence number, Fig. 2 is the display result of 2 milliseconds of samplings, therefore maximum sampling sequence number 11000.From this big gun geological data, extract together, such as the 97th road, Fig. 3 is exactly the 97th road earthquake data that extract, the amplitude of horizontal ordinate Digital Telemetry Seismic Wave among Fig. 3, and ordinate is a sampling sequence number.97 road earthquake data are carried out Short Time Fourier Transform, time-domain signal is transformed to time-frequency domain, transformation results is as shown in Figure 4; Horizontal ordinate express time among Fig. 4, ordinate are the expression frequencies, and the signal of many group Different Slope is arranged among the figure; Having only a kind of slope signal is first-harmonic, is useful signal.

2) will be correlated with before n and m among the time-frequency domain signal X [n, m] as two coordinate axis of rectangular coordinate system, carry out coordinate transform, make first-harmonic parallel: X with time shaft _z[n _z, m _z]=X [n, m];

Step 2) said coordinate transform is at first calculated rotation angle by following (3) formula, then by following (4) and (5) formula to all coordinate transforms among the X [n, m] preceding point (n, m) point (n behind the calculating coordinate change _z, m _z), just accomplished X [n, m] to X _z[n _z, m _z] coordinate transform;

The coordinate transform anglec of rotation is:

${\mathrm{\&theta;}}_0=\arctan \left(\frac{L_s}{(f_e-f_s)(2T+1)\mathrm{\&Delta;t\&Delta;t}}\right)---\left(3\right)$

In the formula: θ ₀Be the coordinate system rotation angle; L _sBe sweep signal length; f _sAnd f _eBe respectively the initial frequency and termination frequency of sweep signal;

Rectangular coordinates transformation becomes polar coordinates to be:

$\left\{\begin{array}{c}r=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}^2+{\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">m</figure-callout>}}^2}\\ \mathrm{\&theta;}=\arctan (m/n)\end{array}\right.---\left(4\right)$

The polar coordinate transform coordinate that meets at right angles is:

$\left\{\begin{array}{c}{n}_z=r\cos (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\\ m_z=r\sin (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\end{array}\right.---\left(5\right)$

In the formula: r is a pole axis length; θ is a polar angle.

Not relevant seismologic record is made up of first-harmonic harmonic two parts; Show as the clinographic curve of Different Slope at time-frequency plane; The purpose of coordinate transform is to do a rotation to these clinographic curves, in postrotational coordinate system, makes first-harmonic parallel with time shaft, and harmonic wave still is a clinographic curve.Fig. 5 is the result who time-frequency domain data shown in Figure 4 is carried out conversion according to coordinate transformation method, the time-frequency figure after the coordinate transform.The horizontal ordinate of data after the coordinate transform does not have the pure time and the notion of frequency, so here we represent the transverse axis and the longitudinal axis to count.Because coordinate transform is the rotation of doing to first-harmonic, the horizontal direction signal is exactly a first-harmonic after the rotation, and the harmonic interference that is of slope is arranged.

3) harmonic interference is separated with signal;

The described removal harmonic interference of step 3) adopts Ka Hunan-Laue husband conversion.Ka Hunan-Laue husband conversion is also referred to as the He Tailing conversion, proper vector conversion or svd etc.Because first-harmonic is very strong in time-frequency domain internal linear correlativity, become next its linear dependence that further increases of level through carrying out it coordinate transform, its singular value just is far longer than non-first-harmonic composition on numerical values recited so.

The described removal harmonic interference of step 3) is accomplished as follows:

1. utilize the Householder conversion with X _zBecome the triple diagonal matrix A of symmetry, computing formula is:

A＝U ^TX _zV (6)

In the formula: X _zBe time-frequency domain geological data not relevant after the coordinate transform, U ^TBe the transposition of U, U=u ₁U ₂K u _k, k=min{N, M-1}, V=v ₁V ₂K v _l, l=min{M, N-2}, each conversion u of U _j(j=1 is with X 2Kk) _zIn element below the j row principal diagonal become 0; Each conversion v of V _i(i=1,2K l) is with X _zIn i become 0 with the element on right minor diagonal element the right of principal diagonal next-door neighbour in capable, calculate following symmetric triple-diagonal matrix A through (6) formula:

$A=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1& {\mathrm{\&beta;}}_1& & \\ {\mathrm{\&beta;}}_1& {\mathrm{\&alpha;}}_2& {\mathrm{\&beta;}}_1& \\ O& O& O& \\ & O& {\mathrm{\&alpha;}}_{N-1}& {\mathrm{\&beta;}}_{N-1}\\ & & {\mathrm{\&beta;}}_{N-1}& {\mathrm{\&alpha;}}_N\end{array}\right],{\mathrm{\&beta;}}_i>0,i=\mathrm{1,2},\mathrm{\&Lambda;},M-1---\left(7\right)$

2. the QL transform method of the implicit displacement of utilization becomes diagonal matrix with the tridiagonal matrix A of symmetry, carries out according to the following steps:

(I) ask the afterbody second order submatrix of A

Proper vector be p ₁And p ₂, get ε=p ₁

(II) make r _M=α _M, C _M=S _M=1, G _M=(α _M-ε) ², η _M=α _M-ε, to i=M, M-1, Λ, 2 by formula (8) calculating

With

With

$\left\{\begin{array}{c}{B}_i=C_i{\mathrm{\&beta;}}_{i-1}^2,F_i=S_i{\mathrm{\&beta;}}_{i-1}^2,\\ {\widehat{\mathrm{\&beta;}}}_i^2=F_i+G_i,\\ S_{i-1}=F_i/{\widehat{\mathrm{\&beta;}}}_{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,C_{i-1}=1-S_{i-1},\\ {\mathrm{\&tau;}}_{i-1}={\mathrm{\&eta;}}_i{\mathrm{\&beta;}}_{i-1}^2/{\widehat{\mathrm{\&beta;}}}_{\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">i</figure-callout>}}^2,\\ {\mathrm{\&rho;}}_{i-1}={\mathrm{\&alpha;}}_{i-1}-r_i,{\mathrm{\&kappa;}}_{i-1}=C_{i-1}-S_{i-1},{\mathrm{\&xi;}}_{i-1}={\mathrm{\&rho;}}_{i-1}{S}_{i-1}\\ {\mathrm{\&eta;}}_{i-1}=C_{i-1}{\mathrm{\&xi;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{\mathrm{\&tau;}}_{i-1},\\ G_{i-1}={\mathrm{\&rho;}}_{i-1}{\mathrm{\&eta;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}({\mathrm{\&tau;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{B}_i),\\ d_{i-1}={\mathrm{\&xi;}}_{i-1}+2{\mathrm{\&tau;}}_{i-1},\\ {\widehat{\mathrm{\&alpha;}}}_i=r_i+d_{i-1}\\ r_{i-1}={\mathrm{\&alpha;}}_{i-1}-d_{i-1}\end{array}\right.---\left(8\right)$

${\widehat{\mathrm{\&alpha;}}}_1=r_1,{\widehat{\mathrm{\&beta;}}}_1^2=G_1,$

(III) to i=M; M-1; Λ; 1, make up new triple diagonal matrix

with

and

$\hat{A}=\left[\begin{array}{cccc}{\widehat{\mathrm{\&alpha;}}}_1& {\widehat{\mathrm{\&beta;}}}_1& & \\ {\widehat{\mathrm{\&beta;}}}_1& {\widehat{\mathrm{\&alpha;}}}_2& {\widehat{\mathrm{\&beta;}}}_1& \\ O& O& O& \\ & O& {\widehat{\mathrm{\&alpha;}}}_{M-1}& {\widehat{\mathrm{\&beta;}}}_{M-1}\\ & & {\widehat{\mathrm{\&beta;}}}_{M-1}& {\widehat{\mathrm{\&alpha;}}}_M\end{array}\right]---\left(9\right)$

(IV) if

then

be exactly desired diagonal matrix, otherwise press the method calculating of I to III again;

3. calculate eigenwert and corresponding left eigenvector and the right proper vector of eigenwert of tridiagonal matrix A;

Described tridiagonal matrix

Diagonal entry be exactly the eigenvalue of tridiagonal matrix A _i,

λ _iTo deserved left eigenvector μ _iAccording to matrix

Ask for λ _iTo deserved right proper vector v _iAccording to matrix

Ask for, wherein

Be X _zTransposition.

4. use the singular value of computes tridiagonal matrix A:

; (i=1; 2; Λ, M); To the descending ordering of singular value;

Fig. 6 is the singular value distribution plan after Ka Hunan-Laue husband conversion, and horizontal ordinate is the size of singular value among the figure, and ordinate is the number of singular value.As can be seen from the figure, singular value is big more, and number is more little, and singular value is more little, and number is big more, and the part that singular value is big is a fundamental signal, and the part that singular value is little is a harmonic interference.

5. separate harmonic interference and signal;

Said separation harmonic interference and signal keep harmonic interference, and a preceding S singular value and characteristic of correspondence vector are composed 0, then remaining M-S singular value and characteristic of correspondence vector thereof are kept.

Said harmonic interference and Signal Separation singular value number s are 1 or 2 or 3.

6. with the singular value and the proper vector reconstruct data thereof that remain.

The corresponding data of each singular value of described reconstruct are pressed following formula:

$I_i={\mathrm{\&sigma;}}_i{\mathrm{\&mu;}}_i{v}_i^T,(i=\mathrm{2,3},\mathrm{\&Lambda;},M)---\left(10\right)$

X _z＝I ₁+I ₂+Λ+I _M (11)

(11) X in the formula is exactly the data of harmonic interference after the reconstruct.

Figure 10 is the time-frequency domain harmonic interference that after Ka Hunan-Laue husband conversion, obtains.As can be seen from the figure, after Ka Hunan-Laue husband conversion, only kept the harmonic interference that slope is arranged.

4) the coordinate inverse transformation of time-frequency domain is with (the n that has a few after the coordinate transform _z, m _z) get back to former coordinate in inverse transformation, the anglec of rotation of coordinate inverse transformation is the same with direct transform;

The said coordinate inverse transformation of step 4) is pressed following formula to X _z[n _z, m _z] middle (n that has a few _z, m _z) (n m), has just accomplished X to calculation level _z[n _z, m _z] to the coordinate inverse transformation of X [n, m];

$\left\{\begin{array}{c}{r}_z=\sqrt{{\mathrm{<figure-callout\; id="2"\; label="n\; z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">n</figure-callout>}}_z^{}+m_{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}\\ {\mathrm{\&theta;}}_z=\arctan (m_z/n_z)\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}n=r_z\cos ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\\ m=r_z\sin ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\end{array}\right.---\left(13\right)$

(12) in the formula: r _zBe pole axis length; θ _zBe polar angle, (13) formula is the method that transforms polar coordinate system rectangular coordinate system.

Figure 11 is the time-frequency domain harmonic interference after the coordinate inverse transformation, can find out after the coordinate inverse transformation that coordinate has in length and breadth had identical meaning with coordinate in length and breadth among Fig. 4, has only a kind of harmonic interference of slope this moment among the figure.

5) do Fourier inversion in short-term;

The described Fourier inversion in short-term of step 5) is accomplished by following formula:

$x\left[k\right]=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}X[n,m]w[k-m]e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000361558100091.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(14\right)$

In the formula: the implication of parameter is with (1) formula.

6) obtain removing the time domain data of harmonic interference;

The time domain data of the said removal harmonic interference of step 6) is accomplished by following (16) formula:

If what step 4) kept is the data of harmonic interference, then:

x″[k]＝x[k]-x′[k] (16)

(16) " [k] is exactly the time domain geological data of removing harmonic interference to x in the formula.Figure 12 extracts harmonic interference, the time-domain signal after time domain is subtracted each other the compacting harmonic interference that obtains then at time-frequency domain.

7) get the relevant preceding geological data of other one field acquisition; Repeating step 2) to the process of step 5); Until handling the not relevant geological data of all field acquisitions; Obtain the downtrodden geological data of harmonic interference, the geological data input correlator after the harmonic interference compacting is correlated with, just obtain the relevant geological data after harmonic interference is suppressed.

The present invention utilizes Short Time Fourier Transform to transform to time-frequency domain to not relevant seismologic record; Using coordinate transform makes first-harmonic be tending towards level; The singular value of horizontal lineups after Ka Hunan-Laue husband conversion is maximum; And the linear harmonic interference of some other inclination is very little with the singular value of noise that some do not have linear behavio(u)r, disturbs according to unusual size separation signal harmonic, and harmonic interference just can well be pressed like this.

The present invention has avoided the deficiency of filtering method in the frequency-wavenumber domain, thereby has reached the purpose of gathering high-quality geological data through vibroseis better.

Figure 13 is the relevant back seismologic record that vibroseis is gathered; Be not use the record that the present invention suppresses harmonic interference; The thing of similar among the figure " herring-bone form " is exactly a harmonic interference, and these interference have reduced the signal to noise ratio (S/N ratio) and the resolution of seismic data, has had a strong impact on the quality of geological data.Figure 14 is the results of same data after the present invention suppresses harmonic interference, and as can be seen from the figure, the harmonic interference of similar " fish-bone " shape thing has been removed, and the quality of seismic data has obtained obvious improvement.

Claims (11)

1. method of suppressing the vibroseis harmonic interference, characteristic provides following technical step:

1) field acquisition geological data is got relevant together preceding geological data and is carried out Short Time Fourier Transform, transforms to time-frequency domain to time-domain signal;

2) will be correlated with before time-sampling sequence number n and frequency sampling sequence number m among the time-frequency domain signal X [n, m] as two coordinate axis of rectangular coordinate system, carry out coordinate transform, make first-harmonic parallel: X with time shaft _z[n _z, m _z]=X [n, m];

3) harmonic interference is separated with signal;

4) the coordinate inverse transformation of time-frequency domain is with (the n that has a few after the coordinate transform _z, m _z) get back to former coordinate in inverse transformation, the anglec of rotation of coordinate inverse transformation is the same with direct transform;

5) do Fourier inversion in short-term;

6) obtain removing the time domain data of harmonic interference;

7) get the relevant preceding geological data of other one field acquisition; Repeating step 2) to the process of step 6); Until handling the relevant geological data of all field acquisitions, obtain the downtrodden geological data of harmonic interference, obtain the relevant geological data after harmonic interference is suppressed after being correlated with.

2. method according to claim 1, characteristics are that the described Short Time Fourier Transform of step 1) is:

$X[n,m]=\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[k\right]w[k-n]e^{-j2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(1\right)$

In the formula:

N and k are the time-sampling sequence number, n=1, and 2,3 ..., N, k=1,2,3 ..., N;

M is the frequency sampling sequence number, m=1, and 2,3 ..., M;

L is the length of geological data before relevant;

Δ t is the time sampling interval;

Δ f is the frequency sampling interval;

N is the number of samples of time-sampling, N=L/ Δ t;

M is the number of samples of frequency sampling, M=1/2 Δ t Δ f;

W [k-n]=w (k Δ t-n Δ t) is the window function of symmetry, the expression formula of window function when window function adopts Blacknam:

$w\left(t\right)=\left\{\begin{array}{c}0.42+0.5\cos (\mathrm{\&pi;t}/T)+0.08\cos (2\mathrm{\&pi;t}/T),t\&Element;[-T,T]\\ 0,t\&NotElement;(-T,T)\end{array}\right.---\left(2\right)$

In the formula: t=k Δ t-n Δ t is the time, and T is the length of half window function;

X [n]=x (n Δ t) is relevant together arbitrarily preceding discrete geological data;

X [n, m]=X (n Δ t, m Δ f) is the Short Time Fourier Transform of discrete geological data before relevant:

The length of described window function is 0.8 second to 1.2 seconds.

3. method according to claim 1, characteristics are steps 2) said coordinate transform at first calculates rotation angle by following (3) formula, then by following (4) and (5) formula to all coordinate transforms among the X [n, m] preceding point (n, m) point (n behind the calculating coordinate change _z, m _z), just accomplished X [n, m] to X _z[n _z, m _z] coordinate transform;

The coordinate transform anglec of rotation is:

${\mathrm{\&theta;}}_0=\arctan \left(\frac{L_s}{(f_e-f_s)(2T+1)\mathrm{\&Delta;t\&Delta;t}}\right)---\left(3\right)$

In the formula: θ ₀Be the coordinate system rotation angle; L _sBe sweep signal length; f _sAnd f _eBe respectively the initial frequency and termination frequency of sweep signal;

Rectangular coordinates transformation becomes polar coordinates to be:

$\left\{\begin{array}{c}r=\sqrt{n^2+m^2}\\ \mathrm{\&theta;}=\arctan (m/n)\end{array}\right.---\left(4\right)$

The polar coordinate transform coordinate that meets at right angles is:

$\left\{\begin{array}{c}{n}_z=r\cos (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\\ m_z=r\sin (\mathrm{\&theta;}-{\mathrm{\&theta;}}_0)\end{array}\right.---\left(5\right)$

In the formula: r is a pole axis length; θ is a polar angle.

4. method according to claim 1, characteristics are that described harmonic interference of step 3) and Signal Separation are accomplished as follows:

1. utilize the Householder conversion with X _zBecome the triple diagonal matrix A of symmetry, computing formula is:

A＝U ^TX _zV (6)

In the formula: X _zBe time-frequency domain geological data not relevant after the coordinate transform, U ^TBe the transposition of U, U=u ₁U ₂K u _k, k=min{N, M-1}, V=v ₁V ₂K vl, l=min{M, N-2}, each conversion u of U _j(j=1,2K k) is with X _zIn element below the j row principal diagonal become 0; Each conversion v of V _i(i=1,2K l) is with X _zIn i become 0 with the element on right minor diagonal element the right of principal diagonal next-door neighbour in capable, calculate following symmetric triple-diagonal matrix A through (6) formula:

$A=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1& {\mathrm{\&beta;}}_1& & \\ {\mathrm{\&beta;}}_1& {\mathrm{\&alpha;}}_2& {\mathrm{\&beta;}}_1& \\ O& O& O& \\ & O& {\mathrm{\&alpha;}}_{N-1}& {\mathrm{\&beta;}}_{N-1}\\ & & {\mathrm{\&beta;}}_{N-1}& {\mathrm{\&alpha;}}_N\end{array}\right],{\mathrm{\&beta;}}_i>0,i=\mathrm{1,2},\mathrm{\&Lambda;},M-1---\left(7\right)$

2. the QL transform method of the implicit displacement of utilization becomes diagonal matrix with the tridiagonal matrix A of symmetry, carries out according to the following steps:

(I) ask the afterbody second order submatrix of A

Proper vector be p ₁And p ₂, get ε=p ₁

(II) make r _M=α _M, C _M=S _M=1, G _M=(α _M-ε) ², η _M=α _M-ε, to i=M, M-1, Λ, 2 by formula (8) calculating

With

With

$\left\{\begin{array}{c}{B}_i=C_i{\mathrm{\&beta;}}_{i-1}^2,F_i=S_i{\mathrm{\&beta;}}_{i-1}^2,\\ {\widehat{\mathrm{\&beta;}}}_i^2=F_i+G_i,\\ S_{i-1}=F_i/{\widehat{\mathrm{\&beta;}}}_i^2,C_{i-1}=1-S_{i-1},\\ {\mathrm{\&tau;}}_{i-1}={\mathrm{\&eta;}}_i{\mathrm{\&beta;}}_{i-1}^2/{\widehat{\mathrm{\&beta;}}}_i^2,\\ {\mathrm{\&rho;}}_{i-1}={\mathrm{\&alpha;}}_{i-1}-r_i,{\mathrm{\&kappa;}}_{i-1}=C_{i-1}-S_{i-1},{\mathrm{\&xi;}}_{i-1}={\mathrm{\&rho;}}_{i-1}{S}_{i-1}\\ {\mathrm{\&eta;}}_{i-1}=C_{i-1}{\mathrm{\&xi;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{\mathrm{\&tau;}}_{i-1},\\ G_{i-1}={\mathrm{\&rho;}}_{i-1}{\mathrm{\&eta;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}({\mathrm{\&tau;}}_{i-1}+{\mathrm{\&kappa;}}_{i-1}{B}_i),\\ d_{i-1}={\mathrm{\&xi;}}_{i-1}+2{\mathrm{\&tau;}}_{i-1},\\ {\widehat{\mathrm{\&alpha;}}}_i=r_i+d_{i-1}\\ r_{i-1}={\mathrm{\&alpha;}}_{i-1}-d_{i-1}\end{array}\right.---\left(8\right)$

${\widehat{\mathrm{\&alpha;}}}_1=r_1,{\widehat{\mathrm{\&beta;}}}_1^2=G_1,$

(III) to i=M; M-1; Λ; 1, make up new triple diagonal matrix

with

and

$\hat{A}=\left[\begin{array}{cccc}{\widehat{\mathrm{\&alpha;}}}_1& {\widehat{\mathrm{\&beta;}}}_1& & \\ {\widehat{\mathrm{\&beta;}}}_1& {\widehat{\mathrm{\&alpha;}}}_2& {\widehat{\mathrm{\&beta;}}}_1& \\ O& O& O& \\ & O& {\widehat{\mathrm{\&alpha;}}}_{M-1}& {\widehat{\mathrm{\&beta;}}}_{M-1}\\ & & {\widehat{\mathrm{\&beta;}}}_{M-1}& {\widehat{\mathrm{\&alpha;}}}_M\end{array}\right]---\left(9\right)$

(IV) if

then

be exactly desired diagonal matrix, otherwise press the method calculating of I to III again;

3. calculate eigenwert and corresponding left eigenvector and the right proper vector of eigenwert of tridiagonal matrix A;

4. use the singular value of computes tridiagonal matrix A:

is to the descending ordering of singular value;

5. harmonic interference and Signal Separation;

6. with the singular value and the proper vector reconstruct data thereof that remain.

5. method according to claim 1, characteristics be the said coordinate inverse transformation of step 4) by following formula to X _z[n _z, m _z] middle (n that has a few _z, m _z) (n m), has just accomplished X to calculation level _z[n _z, m _z] to the coordinate inverse transformation of X [n, m];

$\left\{\begin{array}{c}{r}_z=\sqrt{n_z^2+m_z^2}\\ {\mathrm{\&theta;}}_z=\arctan (m_z/n_z)\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}n=r_z\cos ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\\ m=r_z\sin ({\mathrm{\&theta;}}_z+{\mathrm{\&theta;}}_0)\end{array}\right.---\left(13\right)$

(12) in the formula: r _zBe pole axis length; θ _zBe polar angle, (13) formula is the method that transforms polar coordinate system rectangular coordinate system.

6. method according to claim 1, characteristics are that the described Fourier inversion in short-term of step 5) is accomplished by following formula:

$x^{\&prime;}\left[k\right]=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}X[n,m]w[k-m]e^{-j2\mathrm{\&pi;}\left(\mathrm{m\&Delta;f}\right)\left(\mathrm{k\&Delta;t}\right)}---\left(14\right)$

In the formula: the implication of parameter is with (1) formula.

7. method according to claim 1, characteristics are that the time domain data of the said removal harmonic interference of step 6) is accomplished by following (15) formula or (16) formula:

If what step 4) kept is the data of removing harmonic interference, then:

x″[k]＝x′[k] (15)

If what step 4) kept is the data of harmonic interference, then:

x″[k]＝x[k]-x′[k] (16)

(15) in formula and (16) formula: " [k] is exactly the time domain geological data of removing harmonic interference to x.

8. according to claim 1 or 4 described methods, characteristics are described tridiagonal matrixs

Diagonal entry be exactly the eigenvalue of tridiagonal matrix A _i,

λ _iTo deserved left eigenvector μ _iAccording to matrix

Ask for λ _iTo deserved right proper vector v _iAccording to matrix

Ask for, wherein

Be X _zTransposition.

9. according to claim 1 or 4 described methods, characteristics are said harmonic interference and Signal Separation or stick signal, perhaps keep harmonic interference;

If stick signal then keeps a preceding S singular value and characteristic of correspondence vector, be 0 with remaining M-S singular value and the tax of characteristic of correspondence vector thereof then;

Then a preceding S singular value and characteristic of correspondence vector are composed 0 if keep harmonic interference, then remaining M-S singular value and characteristic of correspondence vector thereof are kept.

10. harmonic interference according to claim 9 and Signal Separation singular value number s are 1 or 2 or 3.

11. according to claim 1 or 4 described methods, characteristics are the corresponding data of each singular value of described reconstruct by following formula:

$I_i={\mathrm{\&sigma;}}_i{\mathrm{\&mu;}}_i{v}_i^T,(i=\mathrm{2,3},\mathrm{\&Lambda;},M)---\left(10\right)$

X _z＝I ₁+I ₂+Λ+I _M (11)

(11) X in the formula is exactly the data of removing the data of harmonic interference after the reconstruct or keeping harmonic interference.

Images