CN103245973B - A kind of method of eliminating marine seismic data wave noise jamming - Google Patents

Abstract

The present invention is exploration, the exploitation in oil field, a kind of method of eliminating marine seismic data wave noise jamming of production technique. Offshore exploration ship air gun source excites and acquiring seismic data, calculate resample filter, carrying out air gun waveform resampling equates the sample rate of air gun waveform and the sample rate of geological data, air gun waveform data are made to Fourier transformation, calculate air gun waveform amplitude spectrum, determine air gun waveform amplitude spectrum maximum, calculate the high-pass filtering factor, territory wave noise attentuation computing time operator, eliminates wave noise jamming in time-domain, frequency domain. The present invention has realized air gun waveform that the marine air gun source of direct use excites and has determined that best wave noise jamming eliminates operator; only geological data wave noise jamming is processed; and effectively protected the low-frequency component of geological data useful signal, have that amount of calculation is little, computational speed fast, good stability and a high feature of computational accuracy.

Description

Method for eliminating sea wave noise interference of marine seismic data

Technical Field

The invention relates to the exploration, development and exploitation technology of oil fields, in particular to a method for eliminating sea wave noise interference of marine seismic data, which provides seismic graphs and data with high resolution and high signal-to-noise ratio for reflecting underground stratum horizons and oil reservoir description.

Background

The marine seismic exploration process is that a marine exploration vessel tows a seismic source (an air gun source) and a recording streamer (in which a geophone is located), and travels at a uniform speed on the sea surface along a pre-designed navigation route. The seismic source position and the geophone position are arranged according to the design requirement, and seismic waves are continuously excited and received in the navigation of the exploration ship. At a series of time points, seismic waves are excited in water by using an air gun seismic source, the seismic waves propagate to the underground, when a wave impedance (the product of the speed of the seismic waves propagating to the underground in a stratum medium and the density of the medium) interface (namely an upper stratum wave impedance and a lower stratum wave impedance are unequal) is met, the seismic waves generate a reflection phenomenon on the wave impedance interface, the propagation direction of the seismic waves is changed, the seismic waves start to propagate upwards, and receivers are arranged at a series of receiving points in a recording streamer to receive the seismic wave data propagating upwards. The offshore seismic exploration work is completed by continuously exciting and receiving the offshore exploration ship. However, the actually received seismic data also contains information of the spatial positions and arrangement positions of the excitation points and the receiving points, various noise interferences, and the like. Seismic data processing is to process the upward-propagating seismic data records in the field exploration process, reserve the information reflecting the wave impedance interface of the underground stratum, and eliminate other information, which is the post-stack seismic data. The post-stack seismic data only reflects the structure and architecture of the subsurface formations.

With the rapid development of offshore seismic exploration, the area and the area of offshore oil and gas exploration are larger and larger, and the requirements on the resolution and the signal-to-noise ratio of seismic data are higher and higher. Sea surging on the sea surface is caused by conditions such as severe weather, and during marine seismic exploration, the seismic recorder can record the sea surging, and the sea surging is noise interference for seismic exploration signals and is called sea wave noise interference. The signal-to-noise ratio of the seismic data is reduced due to the interference of sea wave noise in the seismic records. Especially in shallow water, sometimes the wave noise disturbance is particularly strong. Therefore, such sea wave noise interference must be removed first when processing the marine seismic data.

Sea wave noise interference appears as low frequency vertical banding on seismic data shot gathers. In general seismic data processing, a conventional low-resistance filtering method (i.e., filtering method) is used to remove sea wave noise interference on seismic data (Yilmaz, translation of liuhuashan, etc., seismic data analysis — seismic data processing, inversion and interpretation (upper book), beijing: oil industry press, 2006, P628). The conventional low-resistance filtering method is to design a low-resistance filter according to the low-resistance frequency, filter out the frequency components below the low-resistance frequency, and reserve the frequency components above the low-resistance frequency. In this way, the effective signal frequency components are also filtered out from the filtered-out frequency components below the low-resistance frequency, because the method does not consider the effective signal frequency components and the wave noise interference frequency components, but completely filters them out, the low-frequency effective signal cannot be reserved, and the calculation efficiency is low.

Disclosure of Invention

The invention aims to provide a method for improving the signal-to-noise ratio of seismic data and quickly eliminating sea wave noise interference of marine seismic data.

The technical scheme adopted by the invention comprises the following steps:

1) exciting and acquiring seismic data by using an air gun seismic source of the offshore exploration ship and preprocessing the seismic data;

the preprocessing in the step 1) is to label the seismic data and define an observation system.

2) Calculating a resampling filter;

the resampling filter h [ n Δ τ ] calculation formula in step 2) 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₂Is two frequency parameters of an ideal low-pass filter, delta tau is the time sampling rate of the filter, n is the time sampling sequence number of the filter, M_hThe number of sampling points of positive half-branch time sampling of the filter is 2M_h+1 is the number of filter time sample points.

3) Resampling the airgun wavelet to make the sampling rate of the airgun wavelet equal to that of the seismic data, wherein the calculation formula of the resampling of the airgun wavelet is as follows:

$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 Δ τ ]]For resampling of the air gun wavelet sequence, b [ n Δ t ]]For the resampled air gun wavelet sequence, h [ n Δ t-l Δ τ]For the resampled filter sequence, Δ τ is the air gun wavelet sequence time sampling rate before resampling, Δ t is the air gun wavelet sequence time sampling rate after resampling, and also the time sampling rate of the seismic data, l is the air gun wavelet sequence time sampling sequence number before resampling, n is the air gun wavelet sequence time sampling sequence number after resampling, M is the sampling sequence number of the air gun wavelet sequence after resampling_wIs the number of sampling points of the air gun wavelet sequence time before resampling, M_bThe number of sampling points of the air gun wavelet sequence time after resampling;

4) fourier transform of airgun wavelet data:

$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 ]]Representing air gun wavelet data sequences, B k]Representing a Fourier transform sequence corresponding to the airgun wavelet data; k is the Fourier transform sequence frequency sample number, n is the airgun wavelet sequence time sample number, j represents the imaginary unit, and j²＝-1；W_NDenotes N-point fourier transform factor, N denotes the number of fourier transform sequence samples, and N is 2^m≥(M_b+N_x) M is a suitable positive integer; n is a radical of_xRepresenting the number of seismic data samples, M_bThe number of sampling points of the airgun wavelet sequence time;

5) calculating the amplitude spectrum of the airgun wavelet:

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

in the formula, A [ k ] represents the amplitude spectrum of the airgun wavelet;

6) determining the maximum value of the amplitude spectrum of the airgun wavelet;

the maximum amplitude spectrum value of the amplitude spectrum of the airgun wavelet with k being 0,1,2, Λ, N/2-1 interval and the corresponding frequency sampling serial number are searched, and the formula is as follows:

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

in the formula, A_maxRepresenting the maximum amplitude spectral value, k_maxThe frequency sampling sequence number corresponding to the maximum amplitude spectrum value is represented;

7) calculating a 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 a high pass filter factor;

8) calculating a high-pass filtering attenuation 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]High pass filter attenuation factor, β high pass filter attenuation curve slope, unit is dB/octave, log₂(. cndot.) represents the base 2 logarithm; m_βRepresents the number of truncated samples, an

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

In the formula [ ·]Representing a rounding operation, f₀Representing the low cutoff frequency and af the frequency sampling interval.

9) Calculating a sea wave noise attenuation 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 the wave noise attenuation factor in the frequency domain,is H₃[k]Complex conjugation of (a);

10) calculating a time domain sea wave noise attenuation 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 a time domain sea wave noise attenuation operator, H₃[k]Is the wave noise attenuation factor in the frequency domain, W_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4);

11) wave noise interference is eliminated in a time domain, and the calculation formula is as follows:

$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 ]]Is seismic data containing sea wave noise interference, yn]Is seismic data h after eliminating sea wave noise interference₃[n]Is a time domain sea wave noise attenuation operator; n represents the number of sample points of sea wave noise attenuation operator sequence in time domain, N_xRepresenting the number of seismic data samples; k is a time sampling sequence number of a time domain sea wave noise attenuation operator sequence, and n is a seismic data time sampling sequence number after sea wave noise interference is eliminated;

12) eliminating wave noise interference in a frequency domain;

the method for eliminating the sea wave noise interference in the frequency domain adopts the following steps:

(1) fourier transform of seismic data:

$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 ]]Is seismic data containing sea wave noise interference, X k]Representing Fourier transform sequences, W, of seismic data_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4); n denotes the number of samples of the Fourier transform factor sequence, N_xRepresenting the number of seismic data samples; k is a seismic data Fourier transform sequence frequency sampling sequence number, and n is a seismic data time sampling sequence number;

(2) wave noise interference elimination treatment:

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

in the formula, Y [ k ]]Representing a Fourier transform sequence, H, of the seismic data after the elimination of the wave noise interference₃[k]Is the frequency domain wave noise attenuation factor;

(3) fourier inverse transformation:

$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 ]]Representing the sequence of seismic data after the elimination of the interference of sea wave noise, Yk]Representing a Fourier transform sequence, W, of the seismic data after the elimination of the wave noise interference_NRepresenting N-point Fourier transform factorsCalculating and determining the formula (4); n denotes the number of samples of the Fourier transform factor sequence, N_xRepresenting the number of seismic data samples; k is a seismic data Fourier transform sequence frequency sampling sequence number after sea wave noise interference is eliminated, and n is a seismic data time sampling sequence number after sea wave noise interference is eliminated.

The method is particularly suitable for pre-stack seismic data excited by an offshore air gun seismic source, and achieves the purposes of eliminating the wave noise interference influence generated by waves in the seismic data and improving the signal-to-noise ratio of the seismic data by eliminating the wave noise interference processing of the offshore seismic data.

The method for eliminating sea wave noise interference of marine seismic data realizes that airgun wavelets excited by a marine airgun seismic source are directly used to determine the optimal sea wave noise interference elimination operator. In the conventional wave noise interference removing method, a filtering method is adopted, so that not only wave noise interference is removed, but also low-frequency components of effective signals of seismic data are filtered. According to the method for eliminating the sea wave noise interference of the marine seismic data, only the sea wave noise interference of the seismic data is processed, and the low-frequency component of the effective signal of the seismic data is effectively protected.

The method for eliminating sea wave noise interference of marine seismic data has the characteristics of small calculated amount, high calculating speed, good stability and high calculating precision.

Drawings

FIG. 1 illustrates air gun source wavelets before and after resampling and their log-spectrum comparison

(a) Sampling the wavelets for 2 ms;

(b)2ms sampling wavelet log spectrum;

(c) sampling wavelets for 4 ms;

(d)4ms sampled wavelet log spectrum

FIG. 2 sea wave noise attenuation operator and log spectrum thereof

(a) A sea wave noise attenuation operator;

(b) sea wave noise attenuation operator log spectrum

FIG. 3 comparison of actual shot gather data processing results

(a) Original data;

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

(c) high energy interference;

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

(e) the method of the invention

FIG. 4 comparison of gain effects of actual shot gather data processing results

(a) Original data;

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

(c) high energy interference;

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

(e) the method of the invention

FIG. 5 comparison of log spectra of actual shot gather data processing results

(a) Original data;

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

(c) high energy interference;

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

(e) the method of the invention.

Detailed Description

The invention discloses a method for eliminating sea wave noise interference of marine seismic data, which is to design a low-resistance filter according to a simulated sea air gun seismic source wavelet, filter the sea wave noise interference of the marine seismic data and effectively improve the signal-to-noise ratio of the seismic data. According to the method for eliminating sea wave noise interference of marine seismic data, the filter is determined by calculating the wavelet amplitude spectrum of the air gun seismic source, the filter is used for filtering, the sea wave noise interference can be effectively eliminated, low-frequency effective signals are reserved, the calculation is time-saving and fast, and the calculation efficiency is high.

The invention discloses a method for eliminating sea wave noise interference of marine seismic data, which is realized by adopting the following technical scheme and comprises the following steps:

(1) exciting and acquiring seismic data by using an air gun seismic source of the offshore exploration ship and preprocessing the seismic data;

the preprocessing in the step (1) is to label the seismic data and define an observation system.

(2) Calculating a resampling filter;

when the time sampling rate of the air gun source wavelet is different from the time sampling rate of the seismic data, the air gun source wavelet must be resampled first, so that the time sampling rate of the air gun source wavelet is the same as the time sampling rate of the seismic data. In the process of seismic data processing, the time sampling rate of the air gun source wavelet provided usually is small, and the sampling precision is high, so that the air gun source wavelet sequence needs to be resampled, and the time sampling rate of the air gun source wavelet sequence is the same as the time sampling rate of the seismic data sequence.

The resampling filter is an ideal low-pass filter. The resampling filter h [ n Δ τ ] is calculated by the formula:

$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₂Is two frequency parameters of an ideal low-pass filter, delta tau is the time sampling rate of the filter, n is the time sampling sequence number of the filter, M_hThe number of sampling points of positive half-branch time sampling of the filter is 2M_h+1 is the number of filter time sample points.

(3) Resampling the airgun wavelet;

the air gun wavelet resampling means that the air gun wavelet is resampled, so that the sampling rate of the air gun wavelet is equal to the sampling rate of seismic data. The air gun wavelet resampling calculation formula is as follows:

$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 Δ τ ]]For resampling of the air gun wavelet sequence, b [ n Δ t ]]For the resampled air gun wavelet sequence, h [ n Δ t-l Δ τ]For the resampled filter sequence, Δ τ is the air gun wavelet sequence time sampling rate before resampling, Δ t is the air gun wavelet sequence time sampling rate after resampling, and also the time sampling rate of the seismic data, l is the air gun wavelet sequence time sampling sequence number before resampling, n is the air gun wavelet sequence time sampling sequence number after resampling, M is the sampling sequence number of the air gun wavelet sequence after resampling_wIs the number of sampling points of the air gun wavelet sequence time before resampling, M_bThe number of sampling points of the air gun wavelet sequence time after resampling;

(4) fourier transform of airgun wavelet data:

$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 ]]Representing air gun wavelet data sequences, B k]Representing a Fourier transform sequence corresponding to the airgun wavelet data; k is the Fourier transform sequence frequency sample number, n is the airgun wavelet sequence time sample number, j represents the imaginary unit, and j²＝-1；W_NDenotes N-point fourier transform factor, N denotes the number of fourier transform sequence samples, and N is 2^m≥(M_b+N_x) M is a suitable positive integer; n is a radical of_xRepresenting the number of seismic data samples, M_bThe number of sampling points of the airgun wavelet sequence time;

(5) calculating the amplitude spectrum of the airgun wavelet:

the airgun wavelet amplitude spectrum is the absolute value of the airgun wavelet data Fourier transform sequence. The air gun wavelet amplitude spectrum calculation formula is

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

In the formula, A [ k ] represents the amplitude spectrum of the airgun wavelet;

(6) determining the maximum value of the amplitude spectrum of the airgun wavelet;

the maximum amplitude spectrum value of the amplitude spectrum of the airgun wavelet with k being 0,1,2, Λ, N/2-1 interval and the corresponding frequency sampling serial number are searched, and the formula is as follows:

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

in the formula, A_maxRepresenting the maximum amplitude spectral value, k_maxThe frequency sampling sequence number corresponding to the maximum amplitude spectrum value is represented;

(7) calculating a high-pass filtering factor;

designing a high-pass filter factor, which is calculated by the formula

$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 a high pass filter factor;

(8) calculating a high-pass filtering attenuation factor;

designing a high-pass filter attenuation factor, and calculating the attenuation factor by the formula

$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]High pass filter attenuation factor, β high pass filter attenuation curve slope, unit is dB/octave, log₂(. cndot.) represents the base 2 logarithm; m_βRepresents the number of truncated samples, an

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

In the formula [ ·]Representing a rounding operation, f₀Representing the low cutoff frequency and af the frequency sampling interval.

(9) Calculating a sea wave noise attenuation factor;

the sea wave noise attenuation factor is the product of the high-pass filter factor and the high-pass filter attenuation factor, and the calculation formula is

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 the wave noise attenuation factor in the frequency domain,is H₃[k]Complex conjugation of (a);

(10) calculating a time domain sea wave noise attenuation operator:

the time domain sea wave noise attenuation operator is the Fourier inverse transformation of the frequency domain sea wave noise attenuation factor, and the calculation formula is

$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 a time domain sea wave noise attenuation operator, H₃[k]Is the wave noise attenuation factor in the frequency domain, W_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4);

(11) wave noise interference is eliminated in a time domain, and the calculation formula is as follows:

$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 ]]Is seismic data containing sea wave noise interference, yn]Is seismic data h after eliminating sea wave noise interference₃[n]Is a time domain sea wave noise attenuation operator; n represents the number of sample points of sea wave noise attenuation operator sequence in time domain, N_xRepresenting the number of seismic data samples; k is a time sampling sequence number of a time domain sea wave noise attenuation operator sequence, and n is a seismic data time sampling sequence number after sea wave noise interference is eliminated;

(12) eliminating wave noise interference in a frequency domain;

the method for eliminating the sea wave noise interference in the frequency domain adopts the following steps:

1) fourier transform of seismic data:

$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 ]]Is seismic data containing sea wave noise interference, X k]Representing Fourier transform sequences, W, of seismic data_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4); n denotes the number of samples of the Fourier transform factor sequence, N_xRepresenting the number of seismic data samples; k is a seismic data Fourier transform sequence frequency sampling sequence number, and n is a seismic data time sampling sequence number;

2) wave noise interference elimination treatment:

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

in the formula, Y [ k ]]Representing a Fourier transform sequence, H, of the seismic data after the elimination of the wave noise interference₃[k]Is the frequency domain wave noise attenuation factor;

3) fourier inverse transformation:

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

in the formula, y [ n ]]Representing the sequence of seismic data after the elimination of the interference of sea wave noise, Yk]Representing a Fourier transform sequence, W, of the seismic data after the elimination of the wave noise interference_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4); n represents the number of samples of the fourier transform factor sequence,N_xrepresenting the number of seismic data samples; k is a seismic data Fourier transform sequence frequency sampling sequence number after sea wave noise interference is eliminated, and n is a seismic data time sampling sequence number after sea wave noise interference is eliminated.

(13) Drawing seismic data profile after eliminating wave noise interference and storing and eliminating wave noise interference

The invention carries out processing verification for eliminating wave noise interference.

FIG. 1 is a comparison of air gun source wavelets and their log spectra before and after resampling, (a) is a 2ms sampling wavelet excited by the air gun source during the simulated marine data acquisition; (b) is a 2ms sampled wavelet log spectrum; (c) sampling wavelets 4ms after resampling; (d) is a 4ms sampled wavelet log spectrum.

FIG. 2 is a wave noise attenuation operator and its log spectrum, (a) is a wave noise attenuation operator; (b) is the logarithmic spectrum of the sea wave noise attenuation operator.

FIG. 3 is a comparison of actual shot gather data processing results, (a) raw data; (b) is a high-pass filtering HP (4, 8); (c) is a high energy interference; (d) is high pass filtering + high energy interference; (e) is the process of the present invention. FIG. 4 is a comparison of the gain effect of the actual shot gather data processing results of FIG. 3, (a) is the raw data; (b) is a high-pass filtering HP (4, 8); (c) is a high energy interference; (d) is high pass filtering + high energy interference; (e) is the process of the present invention.

FIG. 5 is a log spectrum comparison of actual shot gather data processing results, (a) is the raw data; (b) is a high-pass filtering HP (4, 8); (c) is a high energy interference; (d) is high pass filtering + high energy interference; (e) is the process of the present invention. Five seismic traces were selected for the log spectrum, with trace numbers 1, 36, 71, 106, and 141, respectively. As can be seen from fig. 3 and 4, strong wave noise interference exists in the original seismic data, and the seismic signals are completely submerged. Although the two methods of high-pass filtering and high-energy interference can eliminate most of the sea wave noise interference, partial sea wave noise interference exists in data, the method effectively eliminates the sea wave noise interference, and meanwhile, as can be seen from a logarithmic spectrum of fig. 5, the method not only completely eliminates the sea wave noise interference, but also effectively retains low-frequency effective signals, the two methods of high-pass filtering and high-energy interference filter most of the sea wave noise interference energy, but retains stronger sea wave noise interference energy, and the method of high-pass filtering and high-energy interference completely filters the sea wave noise interference energy and the low-frequency effective signal energy.

Claims (4)

1. A method for eliminating sea wave noise interference of marine seismic data is characterized by comprising the following steps:

1) exciting and acquiring seismic data by using an air gun seismic source of the offshore exploration ship and preprocessing the seismic data;

2) calculating a resampling filter;

3) resampling the airgun wavelet to make the sampling rate of the airgun wavelet equal to that of the seismic data, wherein the calculation formula of the resampling of the airgun wavelet is as follows:

in the formula, w [ l Δ τ ]]For resampling of the air gun wavelet sequence, b [ n Δ t ]]For the resampled air gun wavelet sequence, h [ n Δ t-l Δ τ]For the resampled filter sequence, Δ τ is the air gun wavelet sequence time sampling rate before resampling, Δ t is the air gun wavelet sequence time sampling rate after resampling, and also the time sampling rate of the seismic data, l is the air gun wavelet sequence time sampling sequence number before resampling, n is the air gun wavelet sequence time sampling sequence number after resampling, M is the sampling sequence number of the air gun wavelet sequence after resampling_wIs the number of sampling points of the air gun wavelet sequence time before resampling, M_bThe number of sampling points of the air gun wavelet sequence time after resampling;

4) fourier transform of airgun wavelet data:

wherein,

W_N＝e^-j2π/N(4)

in the formula, b [ n ]]Representing air gun wavelet data sequences, B k]Representing a Fourier transform sequence corresponding to the airgun wavelet data; k is the Fourier transform sequence frequency sample number, n is the airgun wavelet sequence time sample number, j represents the imaginary unit, and j²＝-1；W_NDenotes N-point fourier transform factor, N denotes the number of fourier transform sequence samples, and N is 2^m≥(M_b+N_x) M is a suitable positive integer; n is a radical of_xRepresenting the number of seismic data samples, M_bThe number of sampling points of the air gun wavelet sequence time after resampling;

5) calculating the amplitude spectrum of the airgun wavelet:

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

in the formula, A [ k ] represents the amplitude spectrum of the airgun wavelet;

6) determining the maximum value of the amplitude spectrum of the airgun wavelet;

the maximum amplitude spectrum value of the amplitude spectrum of the N/2-1 interval airgun wavelet with k being 0,1,2, … and the corresponding frequency sampling serial number are searched, and the formula is as follows:

in the formula, A_maxRepresenting the maximum amplitude spectral value, k_maxThe frequency sampling sequence number corresponding to the maximum amplitude spectrum value is represented;

7) calculating a high-pass filtering factor;

in the formula, H₁[k]Is a high pass filter factor;

8) calculating a high-pass filtering attenuation factor;

in the formula, H₂[k]High pass filter attenuation factor, β high pass filter attenuation curve slope, unit is dB/octave, log₂(. cndot.) represents the base 2 logarithm; m_βRepresents the number of truncated samples, an

In the formula [ ·]Representing a rounding operation, f₀Represents the low cutoff frequency, Δ f represents the frequency sampling interval;

9) calculating a sea wave noise attenuation factor;

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

in the formula, H₃[k]Is the wave noise attenuation factor in the frequency domain,is H₃[k]Complex conjugation of (a);

10) calculating a time domain sea wave noise attenuation operator:

in the formula, h₃[n]Is a time domain sea wave noise attenuation operator, H₃[k]Is the wave noise attenuation factor in the frequency domain, W_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4);

11) wave noise interference is eliminated in a time domain, and the calculation formula is as follows:

in the formula, x [ n ]]Is seismic data containing sea wave noise interference, yn]Is seismic data h after eliminating sea wave noise interference₃[n]Is a time domain sea wave noise attenuation operator; n represents the number of sample points of sea wave noise attenuation operator sequence in time domain, N_xRepresenting the number of seismic data samples; k is a time sampling sequence number of a time domain sea wave noise attenuation operator sequence, and n is a seismic data time sampling sequence number after sea wave noise interference is eliminated;

12) wave noise interference is eliminated in the frequency domain.

2. The method of claim 1, wherein said preprocessing of step 1) is performed by tagging the seismic data to define observation systems.

3. Method according to claim 1, characterized in that said resampling filter h [ n Δ τ ] calculation formula of step 2) is:

in the formula (f)₁And f₂Two frequencies being ideal low pass filtersParameter, Δ τ is the filter time sample rate, n is the filter time sample number, M_hThe number of sampling points of positive half-branch time sampling of the filter is 2M_h+1 is the number of filter time sample points.

4. Method according to claim 1, characterized in that the elimination of the sea wave noise disturbance in the frequency domain according to step 12) employs the following method:

(1) fourier transform of seismic data:

in the formula, x [ n ]]Is seismic data containing sea wave noise interference, X k]Representing Fourier transform sequences, W, of seismic data_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4); n denotes the number of samples of the Fourier transform factor sequence, N_xRepresenting the number of seismic data samples; k is a seismic data Fourier transform sequence frequency sampling sequence number, and n is a seismic data time sampling sequence number;

(2) wave noise interference elimination treatment:

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

in the formula, Y [ k ]]Representing a Fourier transform sequence, H, of the seismic data after the elimination of the wave noise interference₃[k]Is the frequency domain wave noise attenuation factor;

(3) fourier inverse transformation:

in the formula, y [ n ]]Representing seismic data after the elimination of sea wave noise interference, Yk]Representing a Fourier transform sequence, W, of the seismic data after the elimination of the wave noise interference_NExpressing N-point Fourier transform factors, and determining by calculation of formula (4); n denotes the number of samples of the Fourier transform factor sequence, N_xRepresenting the number of seismic data samples; k is the sampling sequence number of the Fourier transform sequence frequency of the seismic data after the wave noise interference is eliminated, n is the elimination of the seaAnd (4) sampling sequence numbers of seismic data after wave noise interference.