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

Method for removing wave noise interferences on offshore earthquake data Download PDF

Info

Publication number
CN103245973A
CN103245973A CN2012100266843A CN201210026684A CN103245973A CN 103245973 A CN103245973 A CN 103245973A CN 2012100266843 A CN2012100266843 A CN 2012100266843A CN 201210026684 A CN201210026684 A CN 201210026684A CN 103245973 A CN103245973 A CN 103245973A
Authority
CN
China
Prior art keywords
air gun
wave noise
sampling
sequence
expression
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN2012100266843A
Other languages
Chinese (zh)
Other versions
CN103245973B (en
Inventor
高少武
祝宽海
陈继红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Petroleum Corp
BGP Inc
Original Assignee
China National Petroleum Corp
BGP Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China National Petroleum Corp, BGP Inc filed Critical China National Petroleum Corp
Priority to CN201210026684.3A priority Critical patent/CN103245973B/en
Publication of CN103245973A publication Critical patent/CN103245973A/en
Application granted granted Critical
Publication of CN103245973B publication Critical patent/CN103245973B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)

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 [ nΔτ ] = sin π ( f 2 + f 1 ) nΔτ sin π ( f 2 - f 1 ) nΔτ ( f 2 - f 1 ) π 2 n 2 Δ τ 2 , -M h≤n≤M h (1)
In the formula, f 1And f 2Be 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 [ nΔt ] = Σ l = 0 M w - 1 w [ lΔτ ] h [ nΔt - 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 [ k ] = Σ n = 0 M b - 1 b [ n ] W N 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 2=-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 = max k max ∈ [ 0 , N / 2 - 1 ] { A [ k ] } - - - ( 6 )
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 [ k ] = A [ k ] A max k = 0,1,2 , Λ , k max - 1 1 k = k max , k max + 1 , Λ , N / 2 - 1 - - - ( 7 )
In the formula, H 1[k] is the high-pass filtering factor;
8) calculate the high-pass filtering decay factor;
H 2 [ k ] = 10 - β 20 log 2 ( M β k + 1 ) k = 0,1,2 , Λ , M β - 1 1 k = M β , M β + 1 , Λ , N / 2 - 1 - - - ( 8 )
In the formula, H 2[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log 2() expression is the logarithm at the end with 2; M βExpression cut-off frequency number of samples, and
M β = [ f 0 Δf ] - - - ( 9 )
In the formula, [] expression rounding operation, f 0The low cut-off frequency rate of expression, Δ f represents that frequency sampling at interval.
9) calculate the wave noise reduction factor;
H 3[k]=H 1[k]H 2[k],k=0,1,2,Λ,N/2-1 (10)
H 3 [ k ] = H ‾ 3 [ N - k - 1 ] , k=N/2,N/2+1,N/2+2,Λ,N-1 (11)
In the formula, H 3[k] is frequency field wave noise reduction factor,
Figure BDA0000134335600000051
Be H 3The complex conjugate of [k];
10) territory wave noise attentuation computing time operator:
h 3 [ n ] = Σ k = 0 N - 1 H 3 [ k ] W N - kn , n=0,1,2,Λ,N-1 (12)
In the formula, h 3[n] is time domain wave noise attentuation operator, H 3[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 [ n ] = Σ k = 0 N - 1 h 3 [ k ] 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 3[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 [ k ] = Σ n = 0 N x - 1 x [ n ] W N 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 3[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 3[k] is frequency field wave noise reduction factor;
(3) Fourier inversion:
y [ n ] = Σ k = 0 N - 1 Y [ k ] W N - 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 [ nΔτ ] = sin π ( f 2 + f 1 ) nΔτ sin π ( f 2 - f 1 ) nΔτ ( f 2 - f 1 ) π 2 n 2 Δ τ 2 , -M h≤n≤M h (1)
In the formula, f 1And f 2Be 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 [ nΔt ] = Σ l = 0 M w - 1 w [ lΔτ ] h [ nΔt - 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 [ k ] = Σ n = 0 M b - 1 b [ n ] W N 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 2=-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 = max k max ∈ [ 0 , N / 2 - 1 ] { A [ k ] } - - - ( 6 )
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 [ k ] = A [ k ] A max k = 0,1,2 , Λ , k max - 1 1 k = k max , k max + 1 , Λ , N / 2 - 1 - - - ( 7 )
In the formula, H 1[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 [ k ] = 10 - β 20 log 2 ( M β k + 1 ) k = 0,1,2 , Λ , M β - 1 1 k = M β , M β + 1 , Λ , N / 2 - 1 - - - ( 8 )
In the formula, H 2[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log 2() expression is the logarithm at the end with 2; M βExpression cut-off frequency number of samples, and
M β = [ f 0 Δf ] - - - ( 9 )
In the formula, [] expression rounding operation, f 0The 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 3[k]=H 1[k]H 2[k],k=0,1,2,Λ,N/2-1 (10)
H 3 [ k ] = H ‾ 3 [ N - k - 1 ] , k=N/2,N/2+1,N/2+2,Λ,N-1 (11)
In the formula, H 3[k] is frequency field wave noise reduction factor, Be H 3The 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 [ n ] = Σ k = 0 N - 1 H 3 [ k ] W N - kn , n=0,1,2,Λ,N-1 (12)
In the formula, h 3[n] is time domain wave noise attentuation operator, H 3[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 [ n ] = Σ k = 0 N - 1 h 3 [ k ] 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 3[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 [ k ] = Σ n = 0 N x - 1 x [ n ] W N 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 3[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 3[k] is frequency field wave noise reduction factor;
3) Fourier inversion:
y [ n ] = Σ k = 0 N - 1 Y [ k ] W N - 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 [ nΔt ] = Σ l = 0 M w - 1 w [ lΔτ ] h [ nΔt - 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 [ k ] = Σ n = 0 M b - 1 b [ n ] W N 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 2=-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 = max k max ∈ [ 0 , N / 2 - 1 ] { A [ k ] } - - - ( 6 )
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 [ k ] = A [ k ] A max k = 0,1,2 , Λ , k max - 1 1 k = k max , k max + 1 , Λ , N / 2 - 1 - - - ( 7 )
In the formula, H 1[k] is the high-pass filtering factor;
8) calculate the high-pass filtering decay factor;
H 2 [ k ] = 10 - β 20 log 2 ( M β k + 1 ) k = 0,1,2 , Λ , M β - 1 1 k = M β , M β + 1 , Λ , N / 2 - 1 - - - ( 8 )
In the formula, H 2[k] is the high-pass filtering decay factor; β is Hi-pass filter die-away curve slope, and unit is dB/octave (decibel/octave), log 2() expression is the logarithm at the end with 2; M βExpression cut-off frequency number of samples, and
M β = [ f 0 Δf ] - - - ( 9 )
In the formula, [] expression rounding operation, f 0The low cut-off frequency rate of expression, Δ f represents that frequency sampling at interval.
9) calculate the wave noise reduction factor;
H 3[k]=H 1[k]H 2[k],k=0,1,2,Λ,N/2-1 (10)
H 3 [ k ] = H ‾ 3 [ N - k - 1 ] , k=N/2,N/2+1,N/2+2,Λ,N-1 (11)
In the formula, H 3[k] is frequency field wave noise reduction factor,
Figure FDA0000134335590000032
Be H 3The complex conjugate of [k];
10) territory wave noise attentuation computing time operator:
h 3 [ n ] = Σ k = 0 N - 1 H 3 [ k ] W N - kn , n=0,1,2,Λ,N-1 (12)
In the formula, h 3[n] is time domain wave noise attentuation operator, H 3[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 [ n ] = Σ k = 0 N - 1 h 3 [ k ] 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 3[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 [ nΔτ ] = sin π ( f 2 + f 1 ) nΔτ sin π ( f 2 - f 1 ) nΔτ ( f 2 - f 1 ) π 2 n 2 Δ τ 2 , -M h≤n≤M h (1)
In the formula, f 1And f 2Be 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 [ k ] = Σ n = 0 N x - 1 x [ n ] W N 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 3[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 3[k] is frequency field wave noise reduction factor;
(3) Fourier inversion:
y [ n ] = Σ k = 0 N - 1 Y [ k ] W N - 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.
CN201210026684.3A 2012-02-07 2012-02-07 A kind of method of eliminating marine seismic data wave noise jamming Active CN103245973B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201210026684.3A CN103245973B (en) 2012-02-07 2012-02-07 A kind of method of eliminating marine seismic data wave noise jamming

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201210026684.3A CN103245973B (en) 2012-02-07 2012-02-07 A kind of method of eliminating marine seismic data wave noise jamming

Publications (2)

Publication Number Publication Date
CN103245973A true CN103245973A (en) 2013-08-14
CN103245973B CN103245973B (en) 2016-05-25

Family

ID=48925606

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201210026684.3A Active CN103245973B (en) 2012-02-07 2012-02-07 A kind of method of eliminating marine seismic data wave noise jamming

Country Status (1)

Country Link
CN (1) CN103245973B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105758795A (en) * 2016-04-11 2016-07-13 山东省科学院海洋仪器仪表研究所 Seawater detection pre-treatment method
CN108780159A (en) * 2015-12-16 2018-11-09 Pgs 地球物理公司 Being operated alone in the subarray of source
CN113655519A (en) * 2021-08-23 2021-11-16 中海石油(中国)有限公司 Method and system for acquiring throttling action coefficient and gas release efficiency parameters of air gun

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2082771A (en) * 1980-08-20 1982-03-10 Mobil Oil Corp F-K Geophysical operations including filtering of seismic records
CN1837859A (en) * 2005-03-25 2006-09-27 中国石油天然气集团公司 Three-dimensional seismic data processing quality monitoring technology
WO2009048683A2 (en) * 2007-10-08 2009-04-16 Schlumberger Canada Limited Controlling seismic source elements based on determining a three-dimensional geometry of the seismic source elements
CN102141634A (en) * 2010-12-17 2011-08-03 西南交通大学 Method for suppressing interference of neutral line of prestack seismic signal based on curvelet transform
CN102262243A (en) * 2010-05-31 2011-11-30 中国石油天然气集团公司 Method for suppressing harmonic interference in seismic data of controlled source by filtering

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2082771A (en) * 1980-08-20 1982-03-10 Mobil Oil Corp F-K Geophysical operations including filtering of seismic records
CN1837859A (en) * 2005-03-25 2006-09-27 中国石油天然气集团公司 Three-dimensional seismic data processing quality monitoring technology
WO2009048683A2 (en) * 2007-10-08 2009-04-16 Schlumberger Canada Limited Controlling seismic source elements based on determining a three-dimensional geometry of the seismic source elements
CN102262243A (en) * 2010-05-31 2011-11-30 中国石油天然气集团公司 Method for suppressing harmonic interference in seismic data of controlled source by filtering
CN102141634A (en) * 2010-12-17 2011-08-03 西南交通大学 Method for suppressing interference of neutral line of prestack seismic signal based on curvelet transform

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108780159A (en) * 2015-12-16 2018-11-09 Pgs 地球物理公司 Being operated alone in the subarray of source
CN105758795A (en) * 2016-04-11 2016-07-13 山东省科学院海洋仪器仪表研究所 Seawater detection pre-treatment method
CN105758795B (en) * 2016-04-11 2018-10-19 山东省科学院海洋仪器仪表研究所 A kind of seawater detection pre-treating method
CN113655519A (en) * 2021-08-23 2021-11-16 中海石油(中国)有限公司 Method and system for acquiring throttling action coefficient and gas release efficiency parameters of air gun
CN113655519B (en) * 2021-08-23 2023-10-13 中海石油(中国)有限公司 Air gun throttling action coefficient and gas release efficiency parameter acquisition method and system

Also Published As

Publication number Publication date
CN103245973B (en) 2016-05-25

Similar Documents

Publication Publication Date Title
CN101598809A (en) A kind of self-adaptation is eliminated the method for linear programming noise and multiple reflection interference
CN100349008C (en) Method for carrying out inversion for wave impedance of earthquake wave
US7243029B2 (en) Systems and methods of hydrocarbon detection using wavelet energy absorption analysis
US9575196B2 (en) Coherent noise attenuation
CN107678062B (en) The integrated forecasting deconvolution of the domain hyperbolic Radon and feedback loop methodology multiple suppression model building method
US20080270033A1 (en) Methods of hydrocarbon detection using spectral energy analysis
US8995223B2 (en) Method for removing Scholte waves and similar ground roll type waves from seismic sea bottom data shallow waters
CN106896409B (en) A kind of varying depth cable ghost reflection drawing method based on wave equation boundary values inverting
CN101620276A (en) Method for attenuation of multiple reflections in seismic data
CN101923176B (en) Method for oil and gas detection by utilizing seismic data instantaneous frequency attribute
WO2012044480A2 (en) Interferometric seismic data processing for a towed marine survey
CN103926622A (en) Method for suppressing multiple waves based on L1 norm multichannel matched filtering
CN112285767B (en) Ocean bottom seismograph four-component ocean surface wave multi-order frequency dispersion energy imaging device and method
CN111290017B (en) Surface wave exploration method for jointly extracting Rayleigh wave frequency dispersion characteristics through seismic electric wave field
CN105116445B (en) A kind of method and device of land and water detector seismic data merging treatment
US20090154291A1 (en) Attenuating Noise in Seismic Data
EP2075597A2 (en) Spectral conditioning for surface seismic data
CN104330826A (en) A method for removing various noises under the condition of complex surface
CN103675910A (en) Amphibious detector seismic data scaling factor retrieval method
CN100349012C (en) Method for pressing random noise in seismological record with low SNR
CN106526678A (en) Reflection acoustic logging wave field separation method and device
CN107436451A (en) A kind of automatic amplitude spectral method for calculating geological data optical cable coupled noise degree of strength
CN102841380B (en) A kind of surface relief complex structure Earthquakes grouping of data is concerned with noise attenuation method
CN102103215B (en) Method for suppressing surface waves of three-dimensional high-density seismic prospecting records before stack
CN103135133A (en) Method and device of vector noise reduction of multi-component seismic data

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant