阅读说明：本技术 海洋地震数据震源子波提取方法、装置和存储介质 (Method and device for extracting seismic source wavelet of marine seismic data and storage medium ) 是由 成谷 于 2021-09-15 设计创作，主要内容包括：本发明公开了海洋地震数据震源子波提取方法、装置和存储介质,该方法包括：从观测海洋地震数据中截取包含鬼波的直达波窗口,以获得各地震道的叠加信号W-(dg)(t),所述叠加信号W-(dg)(t)包含有当前地震道的直达波信号W-(d)(t)和震源鬼波信号W-(g)(t)；对任一地震道的叠加信号W-(dg)(t)进行频域转换,获得当前地震道的频域叠加信号计算当前地震道上鬼波与直达波的到时差Δt；依据当前地震道上鬼波与直达波的到时差Δt,计算不同圆频率ω时当前地震道对应的1-e～(-iωΔt)的数值；依据当前地震道的频域叠加信号和不同圆频率ω时当前地震道对应的1-e～(-iωΔt)的数值,计算当前地震道的频域直达波信号对当前地震道的频域直达波信号进行时域转换,获得当前地震道的地震震源子波W-(d)(t)；本发明能够从直达波和鬼波叠加信号中直接求取地震震源子波。(The invention discloses a method, a device and a storage medium for extracting seismic source wavelets of marine seismic data, wherein the method comprises the following steps: intercepting direct wave window containing ghost waves from observation marine seismic data to obtain superposition signal W of each seismic channel dg (t), the superimposed signal W dg (t) direct wave signal W containing current seismic trace d (t) and seismic ghost signals W g (t); superimposed signal W for any seismic trace dg (t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel Calculating the arrival time difference delta t of ghost waves and direct waves on the current seismic channel; calculating 1-e corresponding to the current seismic channel at different circular frequencies omega according to the arrival time difference delta t of ghost waves and direct waves on the current seismic channel ‑iωΔt The value of (d); frequency domain according to current seismic traceSuperimposed signal 1-e corresponding to current seismic trace at different circular frequencies omega ‑iωΔt Calculating the frequency domain direct wave signal of the current seismic channel For frequency domain direct wave signal of current seismic channel Performing time domain conversion to obtain seismic source wavelet W of current seismic channel d (t); the invention can directly obtain the seismic source wavelet from the direct wave and ghost wave superposed signals.)

1. A marine seismic source wavelet extraction method is characterized by comprising the following steps:

intercepting direct wave window containing ghost waves from observation marine seismic data to obtain superposition signal W of each seismic channel_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)；

Superimposed signal W for any seismic trace_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

Calculating the arrival time difference delta t of ghost waves and direct waves on the current seismic channel;

calculating 1-e corresponding to the current seismic channel at different circular frequencies omega according to the arrival time difference delta t of ghost waves and direct waves on the current seismic channel^-iωΔtThe value of (d);

frequency domain stack signal according to current seismic trace1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

For frequency domain direct wave signal of current seismic channelPerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t)。

2. The method for extracting source wavelets of marine seismic data as claimed in claim 1, wherein said intercepting a direct wave window containing ghost waves from the observed marine seismic data to obtain a stack signal W for each seismic trace_dg(t), further comprising:

and leveling the take-off time on each seismic channel to a standard time value by taking the take-off time of the direct wave on any seismic channel as the standard time value.

3. The method of extracting a source wavelet of marine seismic data as claimed in claim 2, wherein a direct arrival rise time on a first seismic trace is used as said standard time value.

4. The marine earthquake of claim 1The data source wavelet extraction method is characterized in that the superposition signal W of any seismic channel_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channelThe method specifically comprises the following steps:

computing a stack signal W of any seismic trace using Fourier transform_dg(t) obtaining the frequency domain superposition signal of the current seismic channel according to the corresponding frequency domain

5. The method for extracting seismic source wavelets of marine seismic data according to claim 1, wherein the calculating of the arrival time difference Δ t between the ghost and the direct on the current seismic trace specifically comprises:

obtaining the vertical distance z between the seismic source and the water surface from the observation system_sThe vertical distance z between the towing cable and the water surface_gHorizontal distance x between seismic source and hydrophone_sgLet the distance from the seismic source to the hydrophone be d_sgThe distance between the virtual seismic source of the seismic source relative to the water surface and the hydrophone is d_s′gThe propagation velocity of seismic waves in water is V_wThe arrival time of the direct wave on the current seismic channel is t_dThe arrival time of ghost waves on the current seismic channel is t_g；

According to the formulaCalculating the arrival time t of the direct wave on the current seismic channel_dAnd according to the formulaCalculating the arrival time t of ghost waves on the current seismic channel_g；

According to the formulaAnd calculating the arrival time difference delta t of the ghost waves and the direct waves on the current seismic channel.

6. The method of extracting source wavelets for marine seismic data of claim 1, wherein the frequency domain stacking signal is dependent on the current seismic trace1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channelThe method specifically comprises the following steps:

according to the formulaCalculating frequency domain direct wave signal of current seismic channel

7. The method of extracting seismic source wavelets of marine seismic data of claim 1, wherein the frequency domain direct arrival signal for the current seismic tracePerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t), specifically including:

calculating frequency domain direct wave signal of current seismic channel by utilizing Fourier inverse transformationObtaining the seismic source wavelet W of the current seismic channel_d(t)。

8. An extraction device for seismic source wavelets of marine seismic data, comprising:

an acquisition module configured to intercept a direct wave window containing ghost waves from the observed marine seismic data to obtain a stack signal W for each seismic trace_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)；

A first conversion module configured to add the signal W to any seismic trace_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

The first processing module is configured to calculate the arrival time difference delta t of ghost waves and direct waves on the current seismic channel;

a second processing module configured to calculate 1-e corresponding to the current seismic channel at different circular frequencies omega according to the arrival time difference delta t of the ghost wave and the direct wave on the current seismic channel^-iωΔtThe value of (d);

a third processing module configured to superimpose the signal in the frequency domain according to the current seismic trace1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

A second conversion module configured to convert the frequency domain direct wave signal of the current seismic tracePerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t)。

9. The marine seismic source wavelet extraction device of claim 8, wherein the first conversion module is further configured to level the rise time on each seismic trace to a standard time value with the direct arrival rise time on any seismic trace as the standard time value.

10. A storage medium for storing a computer program, characterized in that: the program when executed by a processor implements a method of marine seismic data source wavelet extraction as claimed in any one of claims 1 to 7.

Technical Field

The invention relates to the technical field of seismic data processing, in particular to a method and a device for extracting marine seismic source wavelets and a storage medium.

Background

The wavelet extraction plays an important role in marine seismic data processing. When the seismic wave is excited and transmitted in a stratum medium, the sharp excitation pulse waveform is changed into a seismic source wavelet form with certain duration and limited energy due to the earth filtering effect. In order to improve the resolution of seismic data, wavelet shaping and other processing are generally performed on the form of wavelets, and wavelets are also important parameters in the marine seismic data processing stages such as data filtering, forward modeling, and waveform inversion. Therefore, extracting the seismic source wavelet from the observed marine seismic data has been an important and difficult point in seismic data processing.

As shown in fig. 1, the direct wave in the seismic data is the waveform recorded by the excitation source propagating directly to the receiving cable (streamer) after the seismic source is excited. Because reflection and refraction do not occur at the underground interface, only the speed information of the shallowest layer is carried, and the information of the underground stratum interface is not carried. The direct wave waveform intuitively represents the form of the seismic source signal wavelet, so that the method for extracting the seismic source wavelet by using the direct wave has certain advantages compared with other seismic wave forms such as reflected waves and the like. In addition, the direct wave propagation path is shortest, usually appears earliest on most seismic traces, is not interfered by other waveforms such as reflected waves and refracted waves, and is relatively easy to identify. Therefore, the extraction of seismic source wavelets by direct waves is one of the main research directions in seismic data processing.

In 2019, Lifuyuan and the like disclose an article of 'extracting seismic source wavelets from ocean seismic data direct waves' in petroleum geophysical exploration journals. The method obtains a relation between the direct wave and a seismic source signal by deducing a time distance curve of the direct wave, considering the combined effect of the seismic source and a receiving system and combining a bubble oscillation theory, and provides an analytic solution of calculating a far-field wavelet of the seismic source signal by the direct wave in a frequency domain. When the seismic data receiving position meets the far field condition of the seismic source combined array, the direct wave can be used for calculating the far field wavelet of the seismic source signal.

The direct arrival wave is utilized to extract the seismic source wavelet, which has certain advantages, but the seismic source ghost wave usually exists in marine seismic exploration (in marine seismic exploration, because an excitation source and a receiving cable are both sunk to a certain depth below a water surface, after seismic wave excitation, the direct arrival wave directly propagates to a hydrophone position on the receiving cable, and the direct arrival wave directly propagates to the water surface after excitation and is reflected once, and then propagates to the hydrophone position on the receiving cable, namely the 'seismic source ghost wave' existing along with the direct arrival wave), the difference between the arrival time of the seismic source ghost wave and the direct arrival wave is not large, and the two waves are overlapped together, so that the seismic source wavelet is difficult to directly extract the seismic source wavelet from the direct arrival wave.

Fig. 2 shows a schematic diagram of propagation paths of direct waves and seismic source ghost waves, wherein a downward ray is a direct wave directly propagated to a hydrophone by a seismic source, and a dotted ray is a seismic source ghost wave path directly propagated to the hydrophone after being excited by a virtual seismic source (the path is identical to a path which is reflected once on a water surface and then reaches the hydrophone after being excited by the seismic source) obtained according to the virtual seismic source principle. The depths of subsidence of the excitation source and the receiving cable are small relative to the total propagation path of the seismic waves, so the arrival time of the ghost wave signals of the seismic source is very close to the arrival time of the direct waves. However, ghost waves have a process of reflecting at the water surface, and the reflection coefficient of seismic waves at the water surface is usually-1. Thus, source ghost waves are typically characterized on the seismic traces by the presence of a small time interval after the direct arrival in the opposite sign takeoff direction. And because the wavelet usually has a certain time duration, the direct wave and the ghost wave of the seismic source immediately behind the direct wave are superposed together and cannot be distinguished in a time domain.

Disclosure of Invention

The invention aims to provide a method, a device and a storage medium for extracting seismic source wavelets of marine seismic data, which can directly obtain seismic source wavelets from direct waves and seismic source ghost wave superposed signals so as to serve the processing requirements of subsequent wavelet shaping, forward modeling, waveform inversion and the like.

In order to achieve the purpose, the invention discloses a method for extracting marine seismic source wavelets, which comprises the following steps:

s1, intercepting direct wave window containing ghost waves of the seismic source from the observed marine seismic data to obtain the superposed signal W of each seismic channel_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)；

S2 superposition signal W for any seismic channel_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

S3, calculating the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel;

s4, calculating 1-e corresponding to the current seismic channel at different circular frequencies omega according to the arrival time difference delta t of the seismic source ghost wave and the direct wave on the current seismic channel^-iωΔtThe value of (d);

s5, stacking signals according to frequency domain of current seismic channel1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

S6, frequency domain direct wave signal of current seismic channelPerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t)。

Compared with the prior art, the invention has the advantages that the direct wave signal W containing the current seismic channel_d(t) and seismic ghost signals W_g(t) superimposed signal W_dg(t) on the basis of the current seismic source on the seismic channel by calculatingThe arrival time difference delta t of ghost waves and direct waves and the 1-e corresponding to the current seismic channel at different circular frequencies omega^-iωΔtDirectly obtaining seismic source wavelet W from direct wave and seismic source ghost wave superposed signals_d(t) for serving the processing requirements of subsequent wavelet shaping, forward modeling and waveform inversion, and the seismic source wavelet W_dAnd (t) is simple and quick in extraction mode and is suitable for being realized through an algorithm.

Preferably, in step S1, a direct wave window including ghost waves of the seismic source is intercepted from the observed marine seismic data to obtain a superposition signal W of each seismic trace_dg(t), further comprising:

and S11, leveling the take-off time on each seismic channel to a standard time value by taking the take-off time of the direct wave on any seismic channel as the standard time value.

Specifically, the time of the direct wave takeoff on the first seismic trace is taken as the standard time value.

Preferably, the step S2 specifically includes:

s21, calculating the superposition signal W of any seismic channel by utilizing Fourier transform_dg(t) obtaining the frequency domain superposition signal of the current seismic channel according to the corresponding frequency domain

Preferably, the S3 specifically includes:

s31, acquiring the vertical distance z between the seismic source and the water surface from the observation system_sThe vertical distance z between the towing cable and the water surface_gHorizontal distance x between seismic source and hydrophone_sgLet the distance from the seismic source to the hydrophone be d_sgThe distance between the virtual seismic source of the seismic source relative to the water surface and the hydrophone is d_s′gThe propagation velocity of seismic waves in water is V_wThe arrival time of the direct wave on the current seismic channel is t_dThe arrival time of the ghost wave of the seismic source on the current seismic channel is t_g；

S32, according to the formulaCalculating the arrival time t of the direct wave on the current seismic channel_dAnd according to the formulaCalculating the arrival time t of the ghost wave of the seismic source on the current seismic channel_g；

S33, according to the formulaAnd calculating the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel.

Preferably, the step S5 specifically includes:

s51, according to the formulaCalculating frequency domain direct wave signal of current seismic channel

Preferably, the step S6 specifically includes:

s61, calculating the frequency domain direct wave signal of the current seismic channel by utilizing Fourier inverse transformationObtaining the seismic source wavelet W of the current seismic channel_d(t)。

Correspondingly, the invention also discloses a device for extracting the seismic source wavelet of the marine seismic data, which comprises the following components:

an acquisition module configured to intercept a direct wave window containing seismic source ghost waves from the observed marine seismic data to obtain a superposition signal W of each seismic trace_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)；

A first conversion module configured to add the signal W to any seismic trace_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

The first processing module is configured to calculate the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel;

a second processing module configured to calculate 1-e corresponding to the current seismic trace at different circular frequencies omega according to the arrival time difference delta t of the seismic source ghost wave and the direct wave on the current seismic trace^-iωΔtThe value of (d);

a third processing module configured to superimpose the signal in the frequency domain according to the current seismic trace1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

A second conversion module configured to convert the frequency domain direct wave signal of the current seismic tracePerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t)。

Preferably, the first conversion module is further configured to level the take-off time on each seismic channel to a standard time value by taking the direct wave take-off time on any seismic channel as the standard time value.

Accordingly, the present invention also discloses a storage medium for storing a computer program, which when executed by a processor implements the marine seismic data source wavelet extraction method as described above.

Drawings

FIG. 1 is a schematic diagram of the propagation paths and waveforms of a direct wave and a reflected wave;

FIG. 2 is a schematic diagram of direct and seismic ghost propagation paths;

FIG. 3 is a block flow diagram of a marine seismic data wavelet extraction method of the present invention;

FIG. 4 is a seismic section of the invention in the form of an impulse obtained using numerical simulation;

FIG. 5 is the arrival time difference between the ghost wave and the direct wave of the seismic source obtained from each seismic trace according to the present invention;

FIG. 6 is a seismic section in the form of wavelets obtained using numerical simulation in accordance with the present invention;

FIG. 7 is the flattened seismic section of FIG. 6;

FIG. 8 is an amplitude spectrum of the flattened seismic section of FIG. 7;

FIG. 9 is a unwrapped phase spectrum of the flattened seismic section of FIG. 7;

FIG. 10 is a comparison graph of the processing effectiveness time domain signals of the present invention;

FIG. 11 is a comparison graph of the amplitude spectra of the treatment effect of the present invention;

FIG. 12 is a phase spectrum comparison graph of the processing effect of the present invention;

FIG. 13 is a wavelet morphology map obtained using the present invention for each of the seismic traces in FIG. 7;

FIG. 14 is an amplitude spectrum of a wavelet obtained using the present invention for each seismic trace of the present invention;

FIG. 15 is a phase spectrum of wavelets obtained by the present invention for each seismic trace of the present invention;

FIG. 16 is a block diagram of the structure of the marine seismic data wavelet extraction device of the invention.

Detailed Description

In order to explain technical contents, structural features, and objects and effects of the present invention in detail, the following detailed description is given with reference to the accompanying drawings in conjunction with the embodiments.

Referring to fig. 3, the method for extracting a seismic source wavelet of marine seismic data according to the present embodiment is a method for extracting a marine seismic source wavelet based on a direct wave, and is suitable for directly obtaining a seismic source wavelet from a direct wave and a seismic source ghost superimposed signal.

The seismic source wavelet is characterized in that when seismic waves are excited and then propagate in a stratum medium, a sharp excitation pulse waveform is changed into a form with certain duration and limited energy due to the action of earth filtering, the earth is equivalent to a low-pass filter, and after a broad-band sharp pulse passes through the earth filter, high-frequency components are more lost relatively to low frequency, so that the frequency bandwidth of signals is narrowed, and the duration length of the signals in a time domain is increased.

Referring to fig. 1 to fig. 3, the method for extracting the marine seismic source wavelet of the marine seismic data of the embodiment includes the following steps:

s1, intercepting direct wave window containing ghost waves of the seismic source from the observed marine seismic data to obtain the superposed signal W of each seismic channel_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)。

In fact, since the direct wave take-off time on different seismic traces is different, in order to unify the take-off time and simplify the subsequent calculation, it is preferable that in step S1, a direct wave window containing the ghost wave of the seismic source is intercepted from the observed marine seismic data to obtain the stacking signal W of each seismic trace_dg(t), further comprising:

and S11, leveling the take-off time on each seismic channel to a standard time value by taking the take-off time of the direct wave on any seismic channel as the standard time value.

Specifically, the time of the direct wave takeoff on the first seismic trace is taken as the standard time value. The first seismic channel is the seismic channel with the earliest direct wave rise time in all the seismic channels.

The purpose of leveling the takeoff time of each seismic trace to the standard time value is to facilitate subsequent processing and comparison of seismic trace waveforms, and if the leveling step is not performed, the subsequent superposed signal W may be affected_dg(t) phase spectrum.

S2 superposition signal W for any seismic channel_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

Preferably, the step S2 specifically includes:

s21, Utility foucaultCalculating the superposition signal W of any seismic channel by inner leaf transformation_dg(t) obtaining the frequency domain superposition signal of the current seismic channel according to the corresponding frequency domain

The superposed signal W of the time domain_dg(t) Fourier transform into a frequency domain superimposed signalSo as to be convenient for subsequent calculation on the frequency domain and simplify the calculation process.

And S3, calculating the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel.

Preferably, the S3 specifically includes:

s31, acquiring the vertical distance z between the seismic source and the water surface from the observation system_sThe vertical distance z between the streamer (or hydrophone) and the water surface_gHorizontal distance x between seismic source and hydrophone_sgThe data can be directly acquired from an observation system for marine seismic exploration data acquisition.

The embodiment can directly obtain the seismic source wavelet W of the current seismic channel through calculation on the basis of the existing data_d(t)。

S32, according to the formulaCalculating the arrival time t of the direct wave on the current seismic channel_dAnd according to the formulaCalculating the arrival time t of the ghost wave of the seismic source on the current seismic channel_g。

Wherein d is_sgDistance of seismic source to hydrophone, d_s′gIs the virtual source-to-hydrophone distance of the source relative to the water surface, V_wIs the propagation velocity of seismic waves in water, t_dIs the arrival time, t, of the direct wave on the current seismic trace_gThe time of arrival of the ghost wave of the seismic source on the current seismic trace is shown.

It can be understood that, in marine seismic exploration, the arrival time of the direct wave only carries information of the velocity of the water layer, that is, the arrival time of the direct wave on each seismic channel is the distance from the seismic source to the hydrophone (the distance can be directly calculated according to the relative coordinates of the seismic source and the hydrophone) divided by the propagation velocity of the seismic wave on the water layer (the normal velocity is usually 1500 m/s), so the arrival time of the direct wave and the seismic source ghost wave on each seismic channel can be directly calculated by the above formula.

S33, according to the formulaAnd calculating the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel.

As can be seen from FIG. 2, for the same seismic trace, the arrival time of the seismic source ghost always lags behind the arrival time of the direct wave, so that the arrival time t of the seismic source ghost of the current seismic trace is obtained_gAnd the arrival time t of the direct wave_dAnd performing difference processing to obtain the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel.

S4, calculating 1-e corresponding to the current seismic channel at different circular frequencies omega according to the arrival time difference delta t of the seismic source ghost wave and the direct wave on the current seismic channel^-iωΔtThe numerical value of (c).

S5, stacking signals according to frequency domain of current seismic channel1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

Preferably, the step S5 specifically includes:

s51, according to the formulaCalculating frequency domain direct wave signal of current seismic channel

S6, frequency domain direct wave signal of current seismic channelPerforming time domain conversion to obtain seismic source wavelet W of current seismic channel_d(t)。

Preferably, the step S6 specifically includes:

s61, calculating the frequency domain direct wave signal of the current seismic channel by utilizing Fourier inverse transformationObtaining the seismic source wavelet W of the current seismic channel_d(t)。

For the current seismic trace, the reflection coefficient of the seismic source ghost wave reflected on the water surface is about-1, and the time lag is a direct wave Δ t, so that the seismic source ghost wave signal and the direct wave signal have the following relationship:

W_g(t)＝-W_d(t-Δt)。

at this time, the time domain superposition signal of the direct wave and the seismic source ghost wave can be expressed by the following formula:

W_dg(t)＝W_d(t)+W_g(t)＝W_d(t)-W_d(t-Δt)，

after the time domain signal is Fourier transformed to the frequency domain, according to the property of Fourier transform, formula W_dg(t)＝W_d(t)+W_g(t)＝W_d(t)-W_d(t-Δt)

Obtaining a formula after Fourier transform

WhereinCorresponds to W_dg(t) a Fourier transform of the (t),corresponds to W_d(t) a Fourier transform of the (t),corresponds to W_g(t) Fourier transform, from which can be derived

That is, the frequency spectrum of the direct wave (seismic source wavelet when there is no seismic source ghost interference) is the frequency spectrum of the superposition wave of the direct wave and the seismic source ghost divided by 1-e^-iωΔtWhere ω is the circle frequency. And transforming the frequency spectrum of the direct wave into a time domain by utilizing Fourier inverse transformation to obtain the seismic source wavelet.

Referring to fig. 1-15, in order to verify the effect of the present embodiment, the following embodiment is verified by numerical simulation:

fig. 4 shows a pulse-type seismic section obtained by numerical simulation in this embodiment, where a total of 48 hydrophones (a time-varying vibration pattern received by each hydrophone is referred to as a seismic trace) on a receiving cable is assumed, a vertical distance between a seismic source and a water surface is 12 meters, a vertical distance between a hydrophone and a water surface is 18 meters, a horizontal distance between the seismic source and a first hydrophone is 150 meters, a horizontal distance between hydrophones is 12.5 meters, a propagation speed of a seismic wave in water is 1500 meters per second, a total of 6000 sampling points are provided for each seismic trace of a simulated time-domain seismic section, and a time sampling interval is 0.0001 seconds. In FIG. 4, each seismic trace has a positive pulse and a negative pulse, which respectively represent the arrival time of the direct wave and the ghost wave on each seismic trace. As can be seen from fig. 4, the arrival time difference between the ghost and the direct on the seismic trace becomes smaller as the horizontal distance between the hydrophone and the seismic source increases. Fig. 5 shows the arrival time difference between the ghost wave and the direct wave of the seismic source obtained in each seismic channel, and the horizontal axis shows the horizontal distance between the hydrophone and the seismic source corresponding to each seismic channel.

Fig. 6 is a seismic section of a wavelet form obtained by numerical simulation in this embodiment, and is a result of convolution of fig. 4 and a Ricker wavelet with a dominant frequency of 30Hz, that is, the Ricker wavelet with a dominant frequency of 30Hz (the time duration of the seismic source wavelet is a period value corresponding to the dominant frequency) is a seismic source wavelet form finally required by the numerical simulation in this embodiment, and it indicates that the processing effect of this embodiment is better if the wavelet form obtained by the subsequent processing in this embodiment is consistent with the wavelet form. FIG. 6 shows observation data simulating the effect of direct arrival waves by seismic source ghost waves. As can be seen from fig. 6, under the influence of the seismic source ghost, the Ricker wavelet form with the main frequency of 30Hz cannot be directly seen from the direct arrival wave, and the arrival time difference between the direct arrival wave and the seismic source ghost decreases as the horizontal distance between the hydrophone and the seismic source increases, and the superposed waveform amplitude value is smaller and deviates from the Ricker wavelet form.

The waveform in each seismic channel in fig. 6 is subjected to flattening processing, that is, the take-off time of the direct waves on all the seismic channels is corrected to the take-off time on the 1 st seismic channel, and the seismic profile after flattening is shown in fig. 7. The purpose of flattening is to facilitate subsequent processing and comparison of the waveforms, e.g., not flattening may affect the phase spectrum of the signal. FIGS. 8 and 9 are amplitude and phase spectra of the flattened seismic section, obtained using Fourier transform from FIG. 7. As can be seen from fig. 8, the amplitude spectrum amplitude of the signal on each seismic trace gradually decreases and the frequency band becomes narrower under the influence of the ghost wave of the seismic source. As can be seen from fig. 9, the phase spectra of the signals on the seismic traces are substantially identical due to the flattening process.

FIG. 10 is a graph of the effect time domain signal contrast processed using the example. The solid line in the figure is a Ricker wavelet waveform adopted in numerical simulation and is a time domain wavelet signal form required by the embodiment, and the consistency of the two forms indicates that the effect of the embodiment is better. The broken line is the wavelet obtained by the processing of the present embodiment calculated using the data of the 10 th seismic trace in fig. 7 (the line having the smaller amplitude in fig. 10). Compared with the solid line, the wavelet obtained by the embodiment and the expected wavelet have high consistency in the signal form in the time domain.

Fig. 11 is a comparison graph of amplitude spectra obtained by processing effects of the example, and is obtained by obtaining a frequency spectrum by fourier transform and then further obtaining an amplitude spectrum in fig. 10. In the figure, the solid line represents the amplitude spectrum of a Ricker wavelet with a main frequency of 30Hz used in numerical simulation. The broken line is an amplitude spectrum of the wavelet obtained by the processing of the present embodiment calculated using the data of the 10 th seismic trace in fig. 7 (the line having a smaller amplitude in fig. 10). The dashed line versus the solid line shows that the wavelet obtained by the present embodiment and the desired wavelet have a high degree of consistency in the frequency domain amplitude spectrum.

Fig. 12 is a comparison graph of phase spectra (both phase spectra are obtained after unwrapping) obtained by the processing effect of the example, and fig. 10 is obtained by obtaining a frequency spectrum by fourier transform and then further obtaining a phase spectrum. The solid line in the figure is the phase spectrum of the Ricker wavelet waveform used in the numerical simulation. The dashed line is the phase spectrum of the wavelet obtained by the embodiment process computed using the data for the 10 th seismic trace in FIG. 7 (the line with the smaller amplitude in FIG. 10). Comparing the solid and dashed lines, it can be seen that the wavelet obtained by the embodiment and the desired wavelet have a high degree of conformity in the frequency domain phase spectrum with a constant phase shift from the phase spectrum of the data shown by the green line.

10-13 are only results obtained by using the embodiment of the 10 th trace in FIG. 7, and in order to compare results obtained by different traces, the same processing is performed on each trace in FIG. 7, and the waveform of the seismic source wavelet is shown in FIG. 13. It can be seen that the wavelet morphology obtained by the embodiment has high consistency in the time domain, although the arrival time difference of the direct wave and the ghost wave of the seismic source is different and the waveform of the superposed wave is different in each seismic trace, which is processed by using different seismic traces in fig. 7.

Fig. 14 shows a spectrum obtained by fourier transform in fig. 13, and then further amplitude spectrum is calculated. It can be seen from the figure that although the arrival time differences of the direct wave and the ghost wave of the seismic source on different seismic channels are different, and the waveforms of the superposed waves are different, the wavelets obtained by using the embodiment on different seismic channels have good consistency in terms of main frequency, amplitude, bandwidth and the like on an amplitude spectrum.

Fig. 15 is a diagram obtained by obtaining a frequency spectrum by fourier transform and then further calculating a phase spectrum after unwrapping in fig. 13. It can be seen from the figure that although the arrival time differences of the direct wave and the ghost wave of the seismic source on different seismic channels are different, and the waveforms of the superposed waves are different, the wavelets obtained by using the embodiment on different seismic channels have good consistency on the phase spectrum.

In summary, as can be seen from fig. 10 to fig. 15, the effectiveness of the embodiment is demonstrated by means of numerical simulation.

Referring to fig. 16, correspondingly, the present invention further discloses a device for extracting a seismic source wavelet of marine seismic data, which includes:

an acquisition module 10 configured to intercept a direct wave window containing seismic source ghost waves from the observed marine seismic data to obtain a superposition signal W for each seismic trace_dg(t), the superimposed signal W_dg(t) direct wave signal W containing current seismic trace_d(t) and seismic ghost signals W_g(t)；

A first conversion module 20 configured to add the signal W to any seismic trace_dg(t) performing frequency domain conversion to obtain a frequency domain superposition signal of the current seismic channel

A first processing module 30 configured to calculate a difference Δ t between arrival time of a seismic source ghost and a direct arrival on a current seismic trace;

a second processing module 40 configured to calculate 1-e corresponding to the current seismic trace at different circular frequencies ω according to the arrival time difference Δ t between the ghost wave of the seismic source and the direct wave on the current seismic trace^-iωΔtThe value of (d);

a third processing module 50 configured to superimpose the signals in the frequency domain according to the current seismic trace1-e corresponding to current seismic trace at different circular frequencies omega^-iωΔtCalculating the frequency domain direct wave signal of the current seismic channel

A second conversion module 60 configured to convert the frequency domain direct wave signal of the current seismic tracePerforming time domain conversionObtaining the seismic source wavelet W of the current seismic channel_d(t)。

Preferably, the first conversion module 20 is further configured to level the take-off time on each seismic trace to a standard time value with the direct wave take-off time on any seismic trace as the standard time value.

Accordingly, the present invention also discloses a storage medium for storing a computer program, which when executed by a processor implements the marine seismic data source wavelet extraction method as described above.

With reference to FIGS. 1-16, the present invention applies to direct wave signals W containing current seismic traces_d(t) and seismic ghost signals W_g(t) superimposed signal W_dgOn the basis of (t), calculating the arrival time difference delta t of the ghost wave and the direct wave of the seismic source on the current seismic channel and the 1-e corresponding to the current seismic channel at different circular frequencies omega^-iωΔtDirectly obtaining seismic source wavelet W from direct wave and seismic source ghost wave superposed signals_d(t) for serving the processing requirements of subsequent wavelet shaping, forward modeling and waveform inversion, and the seismic source wavelet W_dAnd (t) is simple and quick in extraction mode and is suitable for being realized through an algorithm.

The above disclosure is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the scope of the present invention, therefore, the present invention is not limited by the appended claims.