CN103245973A - Method for removing wave noise interferences on offshore earthquake data - Google Patents

Abstract

The invention provides a method for removing wave noise interferences on offshore earthquake data, and belongs to oil field exploring, developing and exploiting technology. An offshore exploration ship air gun source provokes and collects earthquake data, re-sampling filters are calculated, air gun wavelet re-sampling enables the sampling rate of the air gun wavelet to be equal to the sampling rate of earthquake data, Fourier transform is performed on the air gun wavelet data, air gun wavelet amplitude spectrum is calculated, the maximum value of the air gun wavelet amplitude spectrum is determined, high-pass filter coefficient is calculated, wave noise attenuation operator within a time domain is calculated, and wave noise interferences are removed within the time domain and the frequency domain. According to the invention, that optimal wave noise interference removing operator is determined by using air gun wavelet provoked by using an offshore air gun source directly can be achieved, not only treatment is performed on wave noise interferences to earthquake data, but also low frequency components of effective signals of earthquake data can be effectively protected, and the method has the advantages of small calculated amount, high computing speed, excellent stability and high computation accuracy.

Description

A kind of method of eliminating marine seismic data wave noise

Technical field

The present invention relates to exploration, exploitation, the production technique in oil field, specifically is to provide the earthquake figure of high resolving power, high s/n ratio and a kind of method of eliminating marine seismic data wave noise of data for reflection subsurface formations layer position, reservoir description.

Background technology

The process of offshore seismic exploration is exactly that the offshore survey ship is pulling focus (air gun source) and record towing cable (wave detector is positioned at wherein), along the navigation route that designs in advance, navigates by water with even speed across the sea.Lay according to designing requirement source location and wave detector position, excites and receive seismic event continuously in the survey vessel navigation.On a series of time points, utilize air gun source earthquake-wave-exciting in water, seismic event is to underground propagation, when running into wave impedance (seismic event in stratum media to the speed of underground propagation and the product of Media density) interface (namely the unequal face of stratum wave impedance) up and down, seismic event produces reflex on the wave impedance interface, the seismic wave propagation direction changes, seismic event begins upwards to propagate, a series of acceptance points in the record towing cable are being settled receiver, receive the seismic data of upwards propagating.Constantly excite and reception work by the offshore survey ship, finish offshore seismic exploration work.Yet the actual geological data that receives is also comprising the information of shot point and acceptance point locus and arrangement position and various noise etc.It is exactly that the seismic data recording of propagating that makes progress in the ground observation process is handled that geological data is handled, and keeps the information at reflection subsurface formations wave impedance interface, and eliminates other information, and this information is exactly post-stack seismic data.Post-stack seismic data only reflects structure and the structure of subsurface formations.

Along with the offshore seismic exploration develop rapidly, the zone of exploratory engineering of off-shore petroleum/gas reservoir and area are also increasing, and be also more and more higher to resolution and the signal to noise ratio (S/N ratio) requirement of seismic data.Because conditions such as inclement weather cause the wave on the sea to spring up, during offshore seismic exploration, seismograph can be recorded to this wave and spring up, and this wave springs up concerning the seismic prospecting signal, is a kind of noise, is called the wave noise.Wave noise in the seismologic record has reduced the signal to noise ratio (S/N ratio) of geological data.Particularly in the shoal water zone, the wave noise is strong especially sometimes.When therefore marine seismic data is handled, must at first remove this wave noise.

The wave noise shows as the low frequency vertical band at geological data big gun collection record.In common geological data is handled, use conventional low-resistance filtering method (being filtering method), remove wave noise (the Yilmaz work on the geological data, Liu Huaishan etc. translate, seismic data is analyzed the processing of----seismic data, inverting and explanation (first volume), Beijing: petroleum industry publishing house, 2006, P628).Conventional low-resistance filtering method is exactly according to the low-resistance frequency, designs a low-cut filter, the frequency content below the low-resistance frequency is filtered out, and the frequency content more than the low-resistance frequency is remained.Like this in the following frequency content of low-resistance frequency that filters out, the useful signal frequency content is also filtered out, because this method is not considered useful signal frequency content and wave noise frequency content, but they is filtered out fully, can not keep the low frequency useful signal, and counting yield is low.

Summary of the invention

The object of the present invention is to provide and a kind ofly improve the geological data signal to noise ratio (S/N ratio), eliminate the method for marine seismic data wave noise fast.

The technical solution used in the present invention may further comprise the steps:

1) excites with acquiring seismic data with offshore survey ship air gun source and do pre-service;

The described pre-service of step 1) refers to geological data is put label, definition recording geometry.

2) calculate resample filter;

Step 2) described resample filter h[n Δ τ] computing formula is:

$h\left[\mathrm{n\Δ\τ}\right]=\frac{\sin \mathrm{\π}(f_2+f_1)\mathrm{n\Δ\τ}\sin \mathrm{\π}(f_2-f_1)\mathrm{n\Δ\τ}}{(f_2-f_1){\mathrm{\π}}^2{n}^2\mathrm{\Δ}{\mathrm{\τ}}^2},$ -M _h≤n≤M _h (1)

In the formula, f ₁And f ₂Be two frequency parameters of ideal low-pass filter, Δ τ is wave filter time-sampling rate, and n is wave filter time-sampling sequence number, M _hBe just half time-sampling sampling point number of wave filter, 2M _hThe+1st, wave filter time-sampling sampling point number.

3) carry out the resampling of air gun wavelet and make the sampling rate of air gun wavelet equate with the sampling rate of geological data, air gun wavelet resampling computing formula is:

$b\left[\mathrm{n\Δt}\right]=\underset{l=0}{\overset{M_w-1}{\mathrm{\Σ}}}w\left[\mathrm{l\Δ\τ}\right]h[\mathrm{n\Δt}-\mathrm{l\Δ\τ}],$ 0≤n≤M _b-1 (2)

In the formula, w[l Δ τ] be air gun wavelet sequence before resampling, b[n Δ t] be the back air gun wavelet sequence that resamples, h[n Δ t-l Δ τ] be the resample filter sequence, Δ τ is air gun wavelet sequence time sampling rate before resampling, Δ t is the back air gun wavelet sequence time sampling rate that resamples, it also is the time-sampling rate of geological data, l is air gun wavelet sequence time sampling sequence number before resampling, and n is the back air gun wavelet sequence time sampling sequence number that resamples, M _wBe air gun wavelet sequence time sampling sampling point number before resampling, M _bIt is the back air gun wavelet sequence time sampling sampling point number that resamples;

4) air gun wavelet data are made Fourier transform:

$B\left[k\right]=\underset{n=0}{\overset{M_b-1}{\mathrm{\Σ}}}b\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (3)

Wherein,

W _N＝e ^-j2π/N (4)

In the formula, b[n] expression air gun wavelet data sequence, B[k] the Fourier transform sequence of expression air gun wavelet data correspondence; K is Fourier transform sequence frequency sampling sequence number, and n is air gun wavelet sequence time sampling sequence number, and j represents imaginary unit, and j ²=-1; W _NThe expression N point Fourier transform factor, N represents Fourier transform sequence number of samples, and N=2 ^m〉=(M _b+ N _x), m is a suitable positive integer; N _xExpression geological data number of samples, M _bIt is air gun wavelet sequence time sampling number of samples;

5) calculate the air gun wavelet amplitude:

A[k]＝|B[k]|，k＝0，1，2，Λ，N-1 (5)

In the formula, A[k] expression air gun wavelet amplitude;

6) determine air gun wavelet amplitude maximal value;

Seek k=0,1,2, Λ, the peak swing spectrum value of the interval air gun wavelet amplitude of N/2-1 with and corresponding frequency sampling sequence number, its formula is:

$A_{\max }=\underset{k_{\max }\∈[0,N/2-1]}{\max }\left\{A\right[k\left]\right\}---\left(6\right)$

In the formula, A _MaxExpression peak swing spectrum value, k _MaxThe frequency sampling sequence number of expression peak swing spectrum value correspondence;

7) calculate the high-pass filtering factor;

$H_1\left[k\right]=\left\{\begin{array}{cc}\frac{A\left[k\right]}{A_{\max }}& k=\mathrm{0,1,2},\mathrm{\Λ},k_{\max }-1\\ 1& k=k_{\max },k_{\max }+1,\mathrm{\Λ},N/2-1\end{array}\right.---\left(7\right)$

In the formula, H ₁[k] is the high-pass filtering factor;

8) calculate the high-pass filtering decay factor;

$H_2\left[k\right]=\left\{\begin{array}{cc}{10}^{-\frac{\mathrm{\β}}{20}{\log }_2\left(\frac{M_{\mathrm{\β}}}{k+1}\right)}& k=\mathrm{0,1,2},\mathrm{\Λ},M_{\mathrm{\β}}-1\\ 1& k=M_{\mathrm{\β}},M_{\mathrm{\β}}+1,\mathrm{\Λ},N/2-1\end{array}\right.---\left(8\right)$

In the formula, H ₂[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log ₂() expression is the logarithm at the end with 2; M _βExpression cut-off frequency number of samples, and

$M_{\mathrm{\β}}=\left[\frac{f_0}{\mathrm{\Δf}}\right]---\left(9\right)$

In the formula, [] expression rounding operation, f ₀The low cut-off frequency rate of expression, Δ f represents that frequency sampling at interval.

9) calculate the wave noise reduction factor;

H ₃[k]＝H ₁[k]H ₂[k]，k＝0，1，2，Λ，N/2-1 (10)

$H_3\left[k\right]={\stackrel{\‾}{H}}_3[N-k-1],$ k＝N/2，N/2+1，N/2+2，Λ，N-1 (11)

In the formula, H ₃[k] is frequency field wave noise reduction factor,

Be H ₃The complex conjugate of [k];

10) territory wave noise attentuation computing time operator:

$h_3\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{H}_3\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，Λ，N-1 (12)

In the formula, h ₃[n] is time domain wave noise attentuation operator, H ₃[k] is frequency field wave noise reduction factor, W _NThe expression N point Fourier transform factor is calculated definite by formula (4);

11) eliminate the wave noise in time domain, computing formula is:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{h}_3\left[k\right]x[n-k],$ n＝0，1，2，Λ，N _x-1 (13)

In the formula, x[n] be the geological data that comprises marine wave noise, y[n] be the geological data after the marine wave noise of elimination, h ₃[n] is time domain wave noise attentuation operator; N represents time domain wave noise attentuation sequence of operators number of samples, N _xExpression geological data number of samples; K is time domain wave noise attentuation sequence of operators time-sampling sequence number, and n is the geological data time-sampling sequence number after the marine wave noise of elimination;

12) eliminate the wave noise in frequency field;

Describedly eliminate the wave noise in frequency field and adopt following method:

(1) geological data Fourier transform:

$X\left[k\right]=\underset{n=0}{\overset{N_x-1}{\mathrm{\Σ}}}x\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (14)

In the formula, x[n] be the geological data that comprises marine wave noise, X[k] expression geological data Fourier transform sequence, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number, and n is geological data time-sampling sequence number;

(2) the wave noise is eliminated and is handled:

Y[k]＝H ₃[k]X[k]，k＝0，1，2，Λ，N-1 (15)

In the formula, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, H ₃[k] is frequency field wave noise reduction factor;

(3) Fourier inversion:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}Y\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，Λ，N _x-1 (16)

In the formula, y[n] expression eliminates geological data sequence after the wave noise, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number after the elimination wave noise, and n is the geological data time-sampling sequence number after the marine wave noise of elimination.

The present invention is specially adapted to the earthquake data before superposition that marine air gun source excites, and handles by eliminating marine seismic data wave noise, reaches and eliminates the wave noise influence that wave produces in the geological data, improves geological data signal to noise ratio (S/N ratio) purpose.

A kind of method of eliminating marine seismic data wave noise of the present invention has realized the air gun wavelet that the marine air gun source of direct use excites, and determines best wave noise elimination operator.And in the conventional wave noise removal method, adopt filtering method, it has not only removed the wave noise, and filtering the low-frequency component of geological data useful signal.A kind of method of eliminating marine seismic data wave noise of the present invention is only handled geological data wave noise, and has effectively protected the low-frequency component of geological data useful signal.

A kind of method of eliminating marine seismic data wave noise of the present invention has that calculated amount is little, computing velocity is fast, good stability and the high characteristics of computational accuracy.

Description of drawings

Air gun source wavelet and logarithmic spectrum thereof contrast before and after Fig. 1 resamples

(a) 2ms sampling wavelet;

(b) 2ms sampling wavelet logarithmic spectrum;

(c) 4ms sampling wavelet;

(d) 4ms sampling wavelet logarithmic spectrum

Fig. 2 wave noise attentuation operator and logarithmic spectrum thereof

(a) wave noise attentuation operator;

(b) wave noise attentuation operator logarithmic spectrum

The actual shot gather data result contrast of Fig. 3

(a) raw data;

(b) high-pass filtering HP (4,8);

(c) high energy disturbs;

(d) high-pass filtering+high energy disturbs;

(e) the inventive method

The contrast of the actual shot gather data result of Fig. 4 gain effect

(a) raw data;

(b) high-pass filtering HP (4,8);

(c) high energy disturbs;

(d) high-pass filtering+high energy disturbs;

(e) the inventive method

The contrast of the actual shot gather data result of Fig. 5 logarithmic spectrum

(a) raw data;

(b) high-pass filtering HP (4,8);

(c) high energy disturbs;

(d) high-pass filtering+high energy disturbs;

(e) the inventive method.

Embodiment

A kind of method of eliminating marine seismic data wave noise of the present invention is exactly the marine air gun source wavelet according to simulation, designs a low-cut filter, and filtering marine seismic data wave noise is with effective raising geological data signal to noise ratio (S/N ratio).A kind of method of eliminating marine seismic data wave noise of the present invention, by calculating the air gun source wavelet amplitude to determine wave filter, handle by this filter filtering, can effectively remove the wave noise, keep the low frequency useful signal, and calculate and to save time, fast the counting yield height.

A kind of method of eliminating marine seismic data wave noise of the present invention adopts following technical scheme to realize, may further comprise the steps:

(1) excites with acquiring seismic data with offshore survey ship air gun source and do pre-service;

The described pre-service of step (1) refers to geological data is put label, definition recording geometry.

(2) calculate resample filter;

When the time-sampling rate of the time-sampling rate of air gun source wavelet and geological data not simultaneously, must at first resample to the air gun source wavelet, make the time-sampling rate of air gun source wavelet identical with the time-sampling rate of geological data.In the geological data processing procedure, usually the air gun source wavelet time-sampling rate that provides is little, therefore the sampling precision height must resample to air gun source wavelet sequence, makes the time-sampling rate of air gun source wavelet sequence identical with the time-sampling rate of geological data sequence.

Resample filter is an ideal low-pass filter.Resample filter h[n Δ τ] computing formula is:

$h\left[\mathrm{n\Δ\τ}\right]=\frac{\sin \mathrm{\π}(f_2+f_1)\mathrm{n\Δ\τ}\sin \mathrm{\π}(f_2-f_1)\mathrm{n\Δ\τ}}{(f_2-f_1){\mathrm{\π}}^2{n}^2\mathrm{\Δ}{\mathrm{\τ}}^2},$ -M _h≤n≤M _h (1)

In the formula, f ₁And f ₂Be two frequency parameters of ideal low-pass filter, Δ τ is wave filter time-sampling rate, and n is wave filter time-sampling sequence number, M _hBe just half time-sampling sampling point number of wave filter, 2M _hThe+1st, wave filter time-sampling sampling point number.

(3) the air gun wavelet resamples;

The air gun wavelet resamples, and refers to the air gun wavelet is carried out resampling, makes the sampling rate of air gun wavelet equate with the sampling rate of geological data.Air gun wavelet resampling computing formula is:

$b\left[\mathrm{n\Δt}\right]=\underset{l=0}{\overset{M_w-1}{\mathrm{\Σ}}}w\left[\mathrm{l\Δ\τ}\right]h[\mathrm{n\Δt}-\mathrm{l\Δ\τ}],$ 0≤n≤M _b-1 (2)

In the formula, w[l Δ τ] be air gun wavelet sequence before resampling, b[n Δ t] be the back air gun wavelet sequence that resamples, h[n Δ t-l Δ τ] be the resample filter sequence, Δ τ is air gun wavelet sequence time sampling rate before resampling, Δ t is the back air gun wavelet sequence time sampling rate that resamples, it also is the time-sampling rate of geological data, l is air gun wavelet sequence time sampling sequence number before resampling, and n is the back air gun wavelet sequence time sampling sequence number that resamples, M _wBe air gun wavelet sequence time sampling sampling point number before resampling, M _bIt is the back air gun wavelet sequence time sampling sampling point number that resamples;

(4) air gun wavelet data are made Fourier transform:

$B\left[k\right]=\underset{n=0}{\overset{M_b-1}{\mathrm{\Σ}}}b\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (3)

Wherein,

W _N＝e ^-j2π/N (4)

In the formula, b[n] expression air gun wavelet data sequence, B[k] the Fourier transform sequence of expression air gun wavelet data correspondence; K is Fourier transform sequence frequency sampling sequence number, and n is air gun wavelet sequence time sampling sequence number, and j represents imaginary unit, and j ²=-1; W _NThe expression N point Fourier transform factor, N represents Fourier transform sequence number of samples, and N=2 ^m〉=(M _b+ N _x), m is a suitable positive integer; N _xExpression geological data number of samples, M _bIt is air gun wavelet sequence time sampling number of samples;

(5) calculate the air gun wavelet amplitude:

The air gun wavelet amplitude is exactly the absolute value of air gun wavelet data Fourier transform sequence.Air gun wavelet amplitude computing formula is

A[k]＝|B[k]|，k＝0，1，2，Λ，N-1 (5)

In the formula, A[k] expression air gun wavelet amplitude;

(6) determine air gun wavelet amplitude maximal value;

Seek k=0,1,2, Λ, the peak swing spectrum value of the interval air gun wavelet amplitude of N/2-1 with and corresponding frequency sampling sequence number, its formula is:

$A_{\max }=\underset{k_{\max }\∈[0,N/2-1]}{\max }\left\{A\right[k\left]\right\}---\left(6\right)$

In the formula, A _MaxExpression peak swing spectrum value, k _MaxThe frequency sampling sequence number of expression peak swing spectrum value correspondence;

(7) calculate the high-pass filtering factor;

Design a high-pass filtering factor, its computing formula is

$H_1\left[k\right]=\left\{\begin{array}{cc}\frac{A\left[k\right]}{A_{\max }}& k=\mathrm{0,1,2},\mathrm{\&Lambda;},k_{\max }-1\\ 1& k=k_{\max },{\mathrm{<figure-callout\; id="1"\; label="k\; max\; +"\; filenames="HDA0000134335610000021.png,HDA0000134335610000041.png"\; state="\{\{state\}\}">k</figure-callout>}}_{\max }+1,\mathrm{\&Lambda;},N/2-1\end{array}\right.---\left(7\right)$

In the formula, H ₁[k] is the high-pass filtering factor;

(8) calculate the high-pass filtering decay factor;

Design a high-pass filtering decay factor, its computing formula is

$H_2\left[k\right]=\left\{\begin{array}{cc}{10}^{-\frac{\mathrm{\&beta;}}{20}{\log }_2\left(\frac{M_{\mathrm{\&beta;}}}{k+1}\right)}& k=\mathrm{0,1,2},\mathrm{\&Lambda;},M_{\mathrm{\&beta;}}-1\\ 1& k=M_{\mathrm{\&beta;}},M_{\mathrm{\&beta;}}+1,\mathrm{\&Lambda;},N/2-1\end{array}\right.---\left(8\right)$

In the formula, H ₂[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log ₂() expression is the logarithm at the end with 2; M _βExpression cut-off frequency number of samples, and

$M_{\mathrm{\&beta;}}=\left[\frac{{\mathrm{<figure-callout\; id="0"\; label="f"\; filenames="HDA0000134335610000011.png,HDA0000134335610000012.png"\; state="\{\{state\}\}">f</figure-callout>}}_0}{\mathrm{\&Delta;f}}\right]---\left(9\right)$

In the formula, [] expression rounding operation, f ₀The low cut-off frequency rate of expression, Δ f represents that frequency sampling at interval.

(9) calculate the wave noise reduction factor;

The wave noise reduction factor is exactly the product of the high-pass filtering factor and high-pass filtering decay factor, and its computing formula is

H ₃[k]＝H ₁[k]H ₂[k]，k＝0，1，2，Λ，N/2-1 (10)

$H_3\left[k\right]={\stackrel{\&OverBar;}{H}}_3[N-k-1],$ k＝N/2，N/2+1，N/2+2，Λ，N-1 (11)

In the formula, H ₃[k] is frequency field wave noise reduction factor, Be H ₃The complex conjugate of [k];

(10) territory wave noise attentuation computing time operator:

Time domain wave noise attentuation operator is exactly the Fourier inversion of frequency field wave noise reduction factor, and its computing formula is

$h_3\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{H}_3\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，Λ，N-1 (12)

In the formula, h ₃[n] is time domain wave noise attentuation operator, H ₃[k] is frequency field wave noise reduction factor, W _NThe expression N point Fourier transform factor is calculated definite by formula (4);

(11) eliminate the wave noise in time domain, computing formula is:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{h}_3\left[k\right]x[n-k],$ n＝0，1，2，Λ，N _x-1 (13)

In the formula, x[n] be the geological data that comprises marine wave noise, y[n] be the geological data after the marine wave noise of elimination, h ₃[n] is time domain wave noise attentuation operator; N represents time domain wave noise attentuation sequence of operators number of samples, N _xExpression geological data number of samples; K is time domain wave noise attentuation sequence of operators time-sampling sequence number, and n is the geological data time-sampling sequence number after the marine wave noise of elimination;

(12) eliminate the wave noise in frequency field;

Describedly eliminate the wave noise in frequency field and adopt following method:

1) geological data Fourier transform:

$X\left[k\right]=\underset{n=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}x\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (14)

In the formula, x[n] be the geological data that comprises marine wave noise, X[k] expression geological data Fourier transform sequence, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number, and n is geological data time-sampling sequence number;

2) the wave noise is eliminated and is handled:

Y[k]＝H ₃[k]X[k]，k＝0，1，2，Λ，N-1 (15)

In the formula, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, H ₃[k] is frequency field wave noise reduction factor;

3) Fourier inversion:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}Y\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，A，N _x-1 (16)

In the formula, y[n] expression eliminates geological data sequence after the wave noise, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number after the elimination wave noise, and n is the geological data time-sampling sequence number after the marine wave noise of elimination.

(13) the wave noise is eliminated in geological data section and the storage after the drafting elimination wave noise

The present invention has carried out elimination wave noise and has handled checking.

Fig. 1 is air gun source wavelet and logarithmic spectrum contrast thereof before and after resampling, the 2ms sampling wavelet that air gun source excites when (a) being the simulation offshore acquisition; (b) be 2ms sampling wavelet logarithmic spectrum; (c) be the back 4ms sampling wavelet that resamples; (d) be 4ms sampling wavelet logarithmic spectrum.

Fig. 2 is wave noise attentuation operator and logarithmic spectrum thereof, (a) is wave noise attentuation operator; (b) be wave noise attentuation operator logarithmic spectrum.

Fig. 3 is an actual shot gather data result contrast, (a) is raw data; (b) be high-pass filtering HP (4,8); (c) be that high energy disturbs; (d) be that high-pass filtering+high energy disturbs; (e) be the inventive method.Fig. 4 is the contrast of the actual shot gather data result of Fig. 3 gain effect, (a) is raw data; (b) be high-pass filtering HP (4,8); (c) be that high energy disturbs; (d) be that high-pass filtering+high energy disturbs; (e) be the inventive method.

Fig. 5 is actual shot gather data result logarithmic spectrum contrast, (a) is raw data; (b) be high-pass filtering HP (4,8); (c) be that high energy disturbs; (d) be that high-pass filtering+high energy disturbs; (e) be the inventive method.Logarithmic spectrum has been chosen five seismic traces, and Taoist monastic name is respectively 1,36,71,106 and 141.By Fig. 3 and Fig. 4 as can be known, exist very strong wave noise in the original earthquake data, flooded seismic signal fully.Though high-pass filtering and high energy disturb two kinds of methods can eliminate most of wave noise, but also there is part wave noise on the data, and the inventive method has been eliminated the wave noise effectively, simultaneously from the logarithmic spectrum of Fig. 5 as can be seen, the present invention has not only eliminated the wave noise fully, also effectively keeping the low frequency useful signal, and high-pass filtering and high energy most of wave noise energy that disturbed two kinds of method filterings, but keeping stronger wave noise energy, high-pass filtering+high energy interference method then fully filtering wave noise energy and low frequency useful signal energy.

Claims (4)

1. method of eliminating marine seismic data wave noise, characteristics are may further comprise the steps:

1) excites with acquiring seismic data with offshore survey ship air gun source and do pre-service;

2) calculate resample filter;

3) carry out the resampling of air gun wavelet and make the sampling rate of air gun wavelet equate with the sampling rate of geological data, air gun wavelet resampling computing formula is:

$b\left[\mathrm{n\&Delta;t}\right]=\underset{l=0}{\overset{M_w-1}{\mathrm{\&Sigma;}}}w\left[\mathrm{l\&Delta;\&tau;}\right]h[\mathrm{n\&Delta;t}-\mathrm{l\&Delta;\&tau;}],$ 0≤n≤M _b-1 (2)

In the formula, w[l Δ τ] be air gun wavelet sequence before resampling, b[n Δ t] be the back air gun wavelet sequence that resamples, h[n Δ t-l Δ τ] be the resample filter sequence, Δ τ is air gun wavelet sequence time sampling rate before resampling, Δ t is the back air gun wavelet sequence time sampling rate that resamples, it also is the time-sampling rate of geological data, l is air gun wavelet sequence time sampling sequence number before resampling, and n is the back air gun wavelet sequence time sampling sequence number that resamples, M _wBe air gun wavelet sequence time sampling sampling point number before resampling, M _bIt is the back air gun wavelet sequence time sampling sampling point number that resamples;

4) air gun wavelet data are made Fourier transform:

$B\left[k\right]=\underset{n=0}{\overset{M_b-1}{\mathrm{\&Sigma;}}}b\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (3)

Wherein,

W _N＝e ^-j2π/N (4)

In the formula, b[n] expression air gun wavelet data sequence, B[k] the Fourier transform sequence of expression air gun wavelet data correspondence; K is Fourier transform sequence frequency sampling sequence number, and n is air gun wavelet sequence time sampling sequence number, and j represents imaginary unit, and j ²=-1; W _NThe expression N point Fourier transform factor, N represents Fourier transform sequence number of samples, and N=2 ^m〉=(M _b+ N _x), m is a suitable positive integer; N _xExpression geological data number of samples, M _bIt is air gun wavelet sequence time sampling number of samples;

5) calculate the air gun wavelet amplitude:

A[k]＝|B[k]|，k＝0，1，2，Λ，N-1 (5)

In the formula, A[k] expression air gun wavelet amplitude;

6) determine air gun wavelet amplitude maximal value;

Seek k=0,1,2, Λ, the peak swing spectrum value of the interval air gun wavelet amplitude of N/2-1 with and corresponding frequency sampling sequence number, its formula is:

$A_{\max }=\underset{k_{\max }\&Element;[0,N/2-1]}{\max }\left\{A\right[k\left]\right\}---\left(6\right)$

In the formula, A _MaxExpression peak swing spectrum value, k _MaxThe frequency sampling sequence number of expression peak swing spectrum value correspondence;

7) calculate the high-pass filtering factor;

$H_1\left[k\right]=\left\{\begin{array}{cc}\frac{A\left[k\right]}{A_{\max }}& k=\mathrm{0,1,2},\mathrm{\&Lambda;},k_{\max }-1\\ 1& k=k_{\max },k_{\max }+1,\mathrm{\&Lambda;},N/2-1\end{array}\right.---\left(7\right)$

In the formula, H ₁[k] is the high-pass filtering factor;

8) calculate the high-pass filtering decay factor;

$H_2\left[k\right]=\left\{\begin{array}{cc}{10}^{-\frac{\mathrm{\&beta;}}{20}{\log }_2\left(\frac{M_{\mathrm{\&beta;}}}{k+1}\right)}& k=\mathrm{0,1,2},\mathrm{\&Lambda;},M_{\mathrm{\&beta;}}-1\\ 1& k=M_{\mathrm{\&beta;}},M_{\mathrm{\&beta;}}+1,\mathrm{\&Lambda;},N/2-1\end{array}\right.---\left(8\right)$

In the formula, H ₂[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log ₂() expression is the logarithm at the end with 2; M _βExpression cut-off frequency number of samples, and

$M_{\mathrm{\&beta;}}=\left[\frac{f_0}{\mathrm{\&Delta;f}}\right]---\left(9\right)$

In the formula, [] expression rounding operation, f ₀The low cut-off frequency rate of expression, Δ f represents that frequency sampling at interval.

9) calculate the wave noise reduction factor;

H ₃[k]＝H ₁[k]H ₂[k]，k＝0，1，2，Λ，N/2-1 (10)

$H_3\left[k\right]={\stackrel{\&OverBar;}{H}}_3[N-k-1],$ k＝N/2，N/2+1，N/2+2，Λ，N-1 (11)

In the formula, H ₃[k] is frequency field wave noise reduction factor,

Be H ₃The complex conjugate of [k];

10) territory wave noise attentuation computing time operator:

$h_3\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{H}_3\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，Λ，N-1 (12)

In the formula, h ₃[n] is time domain wave noise attentuation operator, H ₃[k] is frequency field wave noise reduction factor, W _NThe expression N point Fourier transform factor is calculated definite by formula (4);

11) eliminate the wave noise in time domain, computing formula is:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{h}_3\left[k\right]x[n-k],$ n＝0，1，2，Λ，N _x-1 (13)

In the formula, x[n] be the geological data that comprises marine wave noise, y[n] be the geological data after the marine wave noise of elimination, h ₃[n] is time domain wave noise attentuation operator; N represents time domain wave noise attentuation sequence of operators number of samples, N _xExpression geological data number of samples; K is time domain wave noise attentuation sequence of operators time-sampling sequence number, and n is the geological data time-sampling sequence number after the marine wave noise of elimination;

12) eliminate the wave noise in frequency field.

2. according to the method for claim 1, characteristics are that the described pre-service of step 1) refers to geological data is put label, definition recording geometry.

3. according to the method for claim 1, characteristics are steps 2) described resample filter h[n Δ τ] computing formula is:

$h\left[\mathrm{n\&Delta;\&tau;}\right]=\frac{\sin \mathrm{\&pi;}(f_2+f_1)\mathrm{n\&Delta;\&tau;}\sin \mathrm{\&pi;}(f_2-f_1)\mathrm{n\&Delta;\&tau;}}{(f_2-f_1){\mathrm{\&pi;}}^2{n}^2\mathrm{\&Delta;}{\mathrm{\&tau;}}^2},$ -M _h≤n≤M _h (1)

In the formula, f ₁And f ₂Be two frequency parameters of ideal low-pass filter, Δ τ is wave filter time-sampling rate, and n is wave filter time-sampling sequence number, M _hBe just half time-sampling sampling point number of wave filter, 2M _hThe+1st, wave filter time-sampling sampling point number.

4. according to the method for claim 1, to be that step 12) is described eliminate the wave noise in frequency field and adopt following method characteristics:

(1) geological data Fourier transform:

$X\left[k\right]=\underset{n=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}x\left[n\right]W_N^{\mathrm{kn}},$ k＝0，1，2，Λ，N-1 (14)

In the formula, x[n] be the geological data that comprises marine wave noise, X[k] expression geological data Fourier transform sequence, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number, and n is geological data time-sampling sequence number;

(2) the wave noise is eliminated and is handled:

Y[k]＝H ₃[k]X[k]，k＝0，1，2，Λ，N-1 (15)

In the formula, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, H ₃[k] is frequency field wave noise reduction factor;

(3) Fourier inversion:

$y\left[n\right]=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}Y\left[k\right]W_N^{-\mathrm{kn}},$ n＝0，1，2，Λ，N _x-1 (16)

In the formula, y[n] expression eliminates geological data sequence after the wave noise, Y[k] expression eliminates geological data Fourier transform sequence after the wave noise, W _NThe expression N point Fourier transform factor is calculated definite by formula (4); N represents Fourier transform factor sequence number of samples, N _xExpression geological data number of samples; K is geological data Fourier transform sequence frequency sampling sequence number after the elimination wave noise, and n is the geological data time-sampling sequence number after the marine wave noise of elimination.

Images