阅读说明：本技术 自适应去鬼波与宽频准零相位反褶积联合处理方法及系统 (Adaptive ghost wave removing and broadband quasi-zero phase deconvolution combined processing method and system ) 是由 徐洪斌 周云和 何跃明 于 2020-12-31 设计创作，主要内容包括：本申请公开了一种自适应去鬼波与宽频准零相位反褶积联合处理方法及系统。该方法可以包括：输入原始地震数据；针对原始地震数据进行自适应去鬼波处理,获得鬼波压制后的地震数据；针对鬼波压制后的地震数据进行宽频准零相位反褶积处理,获得宽频处理后的地震数据。本发明解决海洋地震勘探资料子波复杂、宽频处理中可靠拓频及准零相位化的技术问题,提高海洋地震资料处理成果的分辨率和品质,为资料解释提供可靠的基础数据,进而提高勘探开发的成功率,有很好的应用前景。(The application discloses a method and a system for adaptively removing ghost waves and performing wideband quasi-zero phase deconvolution. The method can comprise the following steps: inputting original seismic data; carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data after ghost wave suppression; and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to ghost wave suppression to obtain the seismic data subjected to broadband processing. The invention solves the technical problems of complex wavelets, reliable frequency broadening and quasi-zero phasing in broadband processing of marine seismic exploration data, improves the resolution and quality of marine seismic data processing results, provides reliable basic data for data interpretation, further improves the success rate of exploration and development, and has good application prospect.)

1. A method for adaptive ghost wave removal and wideband quasi-zero phase deconvolution combined processing, comprising:

inputting original seismic data;

carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data subjected to ghost wave suppression;

and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to the ghost suppression to obtain broadband processed seismic data.

2. The method of claim 1, wherein the adaptive de-ghost processing comprises:

performing one-dimensional Fourier transform on the original seismic data to obtain original frequency domain seismic data;

respectively calculating the delay time of the shot-geophone point relative to the primary reflected wave, and further calculating a frequency domain ghost wave removing operator;

calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator to obtain ghost wave removing seismic data;

and performing one-dimensional Fourier inversion on the ghost wave removed seismic data to a t-x domain to obtain the ghost wave suppressed seismic data.

3. The adaptive de-ghosting and wideband quasi-zero phase deconvolution method according to claim 2, wherein the delay time of the shot relative to the primary reflection is calculated by equation (1):

the delay time of the detection point with respect to the primary reflected wave is calculated by equation (2):

wherein, t_sDelay time of shot point relative to primary reflection, t_rIs the delay time of the point of detection relative to the primary reflection, x is the offset of the seismic data trace, d_wIndicating depth of sea floor, d_sDenotes the depth of the excitation point, d_rRepresenting the depth of the wave detection point and v is the speed of the seawater.

4. The adaptive de-ghosting and wideband quasi-zero phase deconvolution method of claim 2, wherein the frequency domain de-ghost operator is calculated by equation (3):

wherein A (f) is a frequency domain ghost wave removing operator, R₀For the sea surface reflection coefficient, ω is the angular frequency, ω ═ 2 π f, and i denotes the complex exponential.

5. The method of claim 1, wherein the wideband quasi-zero phase deconvolution comprises:

for the seismic wavelets of each time window of each channel of data in the shot gather of the seismic data after ghost wave suppression, pulse deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, and then shot gather data of f-x domain is obtained;

calculating the frequency bandwidth of the reflected signal by an energy concentration method;

determining a broadband quasi-zero phase filter operator according to the frequency bandwidth;

and calculating the seismic data after the broadband processing according to the shot gather data of the f-x domain and the broadband quasi-zero phase filtering operator.

6. The adaptive de-ghost and wideband quasi-zero phase deconvolution method of claim 5, wherein calculating the frequency bandwidth of the reflected signal by energy concentration comprises:

calculating an average power spectrum of shot records, and further calculating an extreme point on the average power spectrum;

the frequency corresponding to the extreme point is the central frequency, and the deconvolution expected frequency is determined;

and respectively accumulating power spectrums towards two sides by taking the deconvolution expected frequency as a center, setting an energy threshold value, and calculating the frequency bandwidth.

7. The method of claim 5, wherein the frequency bandwidth is determined by calculating the low-frequency end frequency and the high-frequency end frequency of the effective signal for each time window according to formula (4):

wherein σ is the mean square error of the frequency of the Gaussian-shaped power spectrum, f₀For the center frequency, P (0) is the corresponding power spectrum at zero frequency, and P (f) is the corresponding power spectrum at frequency f.

8. The method of claim 7, wherein the wideband quasi-zero phase filter operator is:

wherein f is_nIs the folding frequency obtained from the data sampling rate.

9. The adaptive de-ghost and wideband quasi-zero phase deconvolution method of claim 5, wherein computing wideband processed seismic data from the f-x domain shot gather data and the wideband quasi-zero phase filter operator comprises:

and calculating the product of the broadband quasi-zero phase filtering operator and the shot gather data of the f-x domain in the frequency domain, and performing one-dimensional Fourier inversion on the calculation result to the t-x domain to obtain the seismic data after broadband processing.

10. An adaptive ghost-canceling and wideband quasi-zero phase deconvolution combined processing system, comprising:

a memory storing computer-executable instructions;

a processor executing computer executable instructions in the memory to perform the steps of:

inputting original seismic data;

carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data subjected to ghost wave suppression;

and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to the ghost suppression to obtain broadband processed seismic data.

Technical Field

The invention relates to the field of processing of geophysical exploration seismic data, in particular to a method and a system for processing self-adaptive ghost wave removing and broadband quasi-zero phase deconvolution in a combined mode.

Background

Ghost waves bring serious adverse effects on the processing quality of marine seismic exploration data, and ghost wave removal is one of the main key points of marine seismic exploration data processing. However, due to the fact that ghost wave cycles are complex and changeable, ghost wave compression difficulty is high, when ghost waves are not completely removed, extracted wavelets are far away from real seismic wavelets and cannot meet wavelet minimum phase hypothesis required by deconvolution, fluctuation characteristics of effective seismic waves are damaged by subsequent prediction deconvolution, pulse deconvolution treatment can cause obvious reduction of treatment quality, deconvolution treatment cannot easily exert due effects of the wavelet deconvolution treatment, and the purpose of high-resolution treatment is difficult to achieve.

Therefore, there is a need to develop a method and a system for adaptive ghost wave elimination and wideband quasi-zero phase deconvolution combined processing, which solve the problem of wideband processing in marine seismic exploration.

The information disclosed in this background section is only for enhancement of understanding of the general background of the invention and should not be taken as an acknowledgement or any form of suggestion that this information forms the prior art already known to a person skilled in the art.

Disclosure of Invention

The invention provides a self-adaptive ghost wave removing and broadband quasi-zero phase deconvolution combined processing method and system, which can solve the technical problems of complex wavelets of marine seismic exploration data, reliable frequency broadening and quasi-zero phasing in broadband processing, improve the resolution and quality of marine seismic data processing results, provide reliable basic data for data interpretation, further improve the success rate of exploration and development and have good application prospect.

In a first aspect, an embodiment of the present disclosure provides a method for joint processing of adaptive ghost wave removal and wideband quasi-zero phase deconvolution, including:

inputting original seismic data;

carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data subjected to ghost wave suppression;

and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to the ghost suppression to obtain broadband processed seismic data.

Preferably, the adaptive de-ghost processing comprises:

performing one-dimensional Fourier transform on the original seismic data to obtain original frequency domain seismic data;

respectively calculating the delay time of the shot-geophone point relative to the primary reflected wave, and further calculating a frequency domain ghost wave removing operator;

calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator to obtain ghost wave removing seismic data;

and performing one-dimensional Fourier inversion on the ghost wave removed seismic data to a t-x domain to obtain the ghost wave suppressed seismic data.

Preferably, the delay time of the shot point with respect to the primary reflection is calculated by the formula (1):

the delay time of the detection point with respect to the primary reflected wave is calculated by equation (2):

wherein, t_sDelay time of shot point relative to primary reflection, t_rIs the delay time of the point of detection relative to the primary reflection, x is the offset of the seismic data trace, d_wIndicating depth of sea floor, d_sDenotes the depth of the excitation point, d_rRepresenting the depth of the wave detection point and v is the speed of the seawater.

Preferably, the frequency domain de-ghost operator is calculated by equation (3):

wherein A (f) is a frequency domain ghost wave removing operator, R₀For the sea surface reflection coefficient, ω is the angular frequency, ω ═ 2 π f, and i denotes the complex exponential.

Preferably, the wideband quasi-zero phase deconvolution process comprises:

for the seismic wavelets of each time window of each channel of data in the shot gather of the seismic data after ghost wave suppression, pulse deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, and then shot gather data of f-x domain is obtained;

calculating the frequency bandwidth of the reflected signal by an energy concentration method;

determining a broadband quasi-zero phase filter operator according to the frequency bandwidth;

and calculating the seismic data after the broadband processing according to the shot gather data of the f-x domain and the broadband quasi-zero phase filtering operator.

Preferably, the calculating the frequency bandwidth of the reflected signal by the energy concentration method includes:

calculating an average power spectrum of shot records, and further calculating an extreme point on the average power spectrum;

the frequency corresponding to the extreme point is the central frequency, and the deconvolution expected frequency is determined;

and respectively accumulating power spectrums towards two sides by taking the deconvolution expected frequency as a center, setting an energy threshold value, and calculating the frequency bandwidth.

Preferably, the low-frequency end frequency and the high-frequency end frequency of the effective signal of each time window are calculated by formula (4), and the frequency bandwidth is determined:

wherein σ is the mean square error of the frequency of the Gaussian-shaped power spectrum, f₀For the center frequency, P (0) is the corresponding power spectrum at zero frequency, and P (f) is the corresponding power spectrum at frequency f.

Preferably, the wideband quasi-zero phase filter operator is:

wherein f is_nIs the folding frequency obtained from the data sampling rate.

Preferably, calculating the broadband processed seismic data according to the f-x domain shot gather data and the broadband quasi-zero phase filter operator includes:

and calculating the product of the broadband quasi-zero phase filtering operator and the shot gather data of the f-x domain in the frequency domain, and performing one-dimensional Fourier inversion on the calculation result to the t-x domain to obtain the seismic data after broadband processing.

In a second aspect, an embodiment of the present disclosure further provides an adaptive ghost wave removing and wideband quasi-zero phase deconvolution combined processing system, including:

inputting original seismic data;

carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data subjected to ghost wave suppression;

and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to the ghost suppression to obtain broadband processed seismic data.

Preferably, the adaptive de-ghost processing comprises:

performing one-dimensional Fourier transform on the original seismic data to obtain original frequency domain seismic data;

respectively calculating the delay time of the shot-geophone point relative to the primary reflected wave, and further calculating a frequency domain ghost wave removing operator;

calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator to obtain ghost wave removing seismic data;

and performing one-dimensional Fourier inversion on the ghost wave removed seismic data to a t-x domain to obtain the ghost wave suppressed seismic data.

Preferably, the delay time of the shot point with respect to the primary reflection is calculated by the formula (1):

the delay time of the detection point with respect to the primary reflected wave is calculated by equation (2):

wherein, t_sDelay time of shot point relative to primary reflection, t_rIs the delay time of the point of detection relative to the primary reflection, x is the offset of the seismic data trace, d_wIndicating depth of sea floor, d_sDenotes the depth of the excitation point, d_rRepresenting the depth of the wave detection point and v is the speed of the seawater.

Preferably, the frequency domain de-ghost operator is calculated by equation (3):

wherein A (f) is a frequency domain ghost wave removing operator, R₀For the sea surface reflection coefficient, ω is the angular frequency, ω ═ 2 π f, and i denotes the complex exponential.

Preferably, the wideband quasi-zero phase deconvolution process comprises:

for the seismic wavelets of each time window of each channel of data in the shot gather of the seismic data after ghost wave suppression, pulse deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, and then shot gather data of f-x domain is obtained;

calculating the frequency bandwidth of the reflected signal by an energy concentration method;

determining a broadband quasi-zero phase filter operator according to the frequency bandwidth;

and calculating the seismic data after the broadband processing according to the shot gather data of the f-x domain and the broadband quasi-zero phase filtering operator.

Preferably, the calculating the frequency bandwidth of the reflected signal by the energy concentration method includes:

calculating an average power spectrum of shot records, and further calculating an extreme point on the average power spectrum;

the frequency corresponding to the extreme point is the central frequency, and the deconvolution expected frequency is determined;

and respectively accumulating power spectrums towards two sides by taking the deconvolution expected frequency as a center, setting an energy threshold value, and calculating the frequency bandwidth.

Preferably, the low-frequency end frequency and the high-frequency end frequency of the effective signal of each time window are calculated by formula (4), and the frequency bandwidth is determined:

wherein σ is the mean square error of the frequency of the Gaussian-shaped power spectrum, f₀For the center frequency, P (0) is the corresponding power spectrum at zero frequency, and P (f) is the corresponding power spectrum at frequency f.

Preferably, the wideband quasi-zero phase filter operator is:

wherein f is_nIs the folding frequency obtained from the data sampling rate.

Preferably, calculating the broadband processed seismic data according to the f-x domain shot gather data and the broadband quasi-zero phase filter operator includes:

and calculating the product of the broadband quasi-zero phase filtering operator and the shot gather data of the f-x domain in the frequency domain, and performing one-dimensional Fourier inversion on the calculation result to the t-x domain to obtain the seismic data after broadband processing.

The beneficial effects are that:

aiming at solving the technical difficulty of ocean wavelet frequency broadening, the invention provides a ocean seismic data broadband processing method based on the combination of self-adaptive ghost wave removing and broadband quasi-zero phase deconvolution, according to the characteristic that the self-adaptive ghost wave removing technology has good ghost wave removing and broadband quasi-zero phase deconvolution and has obvious frequency broadening. The invention can better solve the technical problems of complex wavelet of marine seismic exploration data, reliable frequency broadening and quasi-zero phasing in broadband processing, can improve the resolution and quality of marine seismic data processing results, provides reliable basic data for data interpretation, further improves the success rate of exploration and development, and has good application prospect.

The method and apparatus of the present invention have other features and advantages which will be apparent from or are set forth in detail in the accompanying drawings and the following detailed description, which are incorporated herein, and which together serve to explain certain principles of the invention.

Drawings

The above and other objects, features and advantages of the present invention will become more apparent by describing in more detail exemplary embodiments thereof with reference to the attached drawings, in which like reference numerals generally represent like parts.

FIG. 1 is a flow diagram illustrating the steps of a method for joint processing of adaptive de-ghost and wideband quasi-zero phase deconvolution, according to one embodiment of the invention.

FIGS. 2a and 2b show a schematic of wavelets and spectra, respectively, of raw seismic data, according to one embodiment of the invention.

FIGS. 3a and 3b show a schematic of wavelets and spectra, respectively, of ghost-suppressed seismic data, according to one embodiment of the invention.

FIGS. 4a and 4b show a schematic of wavelets and spectra, respectively, of broadband processed seismic data according to one embodiment of the invention.

Fig. 5a, 5b, and 5c are schematic diagrams illustrating comparison of the monochroom recordings before processing, after de-ghost only, and after de-ghost and wideband quasi-zero phase deconvolution processing, respectively, according to an embodiment of the present invention.

Fig. 6a, 6b, and 6c are schematic diagrams showing PSTM profiles before processing, after depoisoning only, and after depoisoning and wideband quasi-zero phase deconvolution processing, respectively, according to an embodiment of the present invention.

Detailed Description

Preferred embodiments of the present invention will be described in more detail below. While the following describes preferred embodiments of the present invention, it should be understood that the present invention may be embodied in various forms and should not be limited by the embodiments set forth herein.

FIG. 1 is a flow diagram illustrating the steps of a method for joint processing of adaptive de-ghost and wideband quasi-zero phase deconvolution, according to one embodiment of the invention.

The invention provides a self-adaptive ghost wave removing and broadband quasi-zero phase deconvolution combined processing method, which comprises the following steps:

step 101, inputting original seismic data;

specifically, the seismic data loaded with the observation system information channel header is input, the original seismic data includes ghost waves, and the wavelets and the frequency spectrums of the ghost waves are respectively shown in fig. 2a and fig. 2 b.

102, carrying out self-adaptive ghost wave removing processing on original seismic data to obtain seismic data subjected to ghost wave suppression; in one example, adaptive de-ghosting processing includes:

performing one-dimensional Fourier transform on the original seismic data to obtain original frequency domain seismic data;

respectively calculating the delay time of the shot-geophone point relative to the primary reflected wave, and further calculating a frequency domain ghost wave removing operator;

calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator to obtain ghost wave removing seismic data;

and performing one-dimensional Fourier inversion on the seismic data without the ghost waves to a t-x domain to obtain the seismic data after ghost wave suppression.

In one example, the delay time of the shot point with respect to the primary reflection is calculated by equation (1):

the delay time of the detection point with respect to the primary reflected wave is calculated by equation (2):

wherein, t_sDelay time of shot point relative to primary reflection, t_rFor the detection point to be opposite to the primaryDelay time of the seismic wave, x being the offset of the seismic data trace, d_wIndicating depth of sea floor, d_sDenotes the depth of the excitation point, d_rRepresenting the depth of the wave detection point and v is the speed of the seawater.

In one example, the frequency domain de-ghost operator is computed by equation (3):

wherein A (f) is a frequency domain ghost wave removing operator, R₀For the sea surface reflection coefficient, ω is the angular frequency, ω ═ 2 π f, and i denotes the complex exponential.

Specifically, assuming that the seismic record x (t) is formed by stacking primary waves and various types of ghost waves, and s (t) is a primary reflection wave in the seismic record, the seismic record can be represented as:

wherein R is₀Is the sea surface reflection coefficient, approximately-1, t_sAnd t_rFor the delay time of the shot point relative to the primary reflected wave, performing one-dimensional Fourier transform on the formula (6) to obtain the seismic data of the original frequency domain as follows:

the delay time of the shot point with respect to the primary reflected wave is calculated by equations (1) and (2), respectively.

The ghost can be obtained according to equation (7):

then equation (7) can be written as:

X(f)＝S(f)G(f) (9)

as can be seen from equation (8), the original frequency domain seismic data with ghost interference can be regarded as the product of the primary reflection and a filter operator g (f). Therefore, an inverse filter operator a (f) is designed according to the formula (9) as a formula (3), namely, a frequency domain ghost wave removing operator is designed, and ghost wave interference in the original frequency domain seismic data can be removed.

FIGS. 3a and 3b show a schematic of wavelets and spectra, respectively, of ghost-suppressed seismic data, according to one embodiment of the invention.

And calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator in the frequency domain to obtain ghost wave removing seismic data. And performing one-dimensional Fourier inversion on the ghost wave-removed seismic data without the influence of the ghost wave notch to a t-x domain to obtain the seismic data after ghost wave suppression, wherein wavelets and frequency spectrums of the seismic data are respectively shown in fig. 3a and 3 b.

And 103, performing broadband quasi-zero phase deconvolution processing on the seismic data subjected to ghost suppression to obtain broadband processed seismic data. In one example, the wideband quasi-zero phase deconvolution process includes:

aiming at the seismic wavelets of each time window of each channel of data in the shot gather of the seismic data after ghost wave suppression, pulse deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, and then shot gather data of f-x domain is obtained;

calculating the frequency bandwidth of the reflected signal by an energy concentration method;

determining a broadband quasi-zero phase filter operator according to the frequency bandwidth;

and calculating the seismic data after the broadband processing according to the shot gather data of the f-x domain and the broadband quasi-zero phase filtering operator.

In one example, calculating the frequency bandwidth of the reflected signal by energy concentration comprises:

calculating the average power spectrum of the shot record, and further calculating an extreme point on the average power spectrum;

determining the deconvolution expected frequency by taking the frequency corresponding to the extreme point as a central frequency;

and accumulating the power spectrums towards two sides by taking the deconvolution expected frequency as a center, setting an energy threshold value and calculating the bandwidth.

In one example, the low-frequency end frequency and the high-frequency end frequency of the effective signal of each time window are calculated by formula (4), and the frequency bandwidth is determined:

wherein σ is the mean square error of the frequency of the Gaussian-shaped power spectrum, f₀For the center frequency, P (0) is the corresponding power spectrum at zero frequency, and P (f) is the corresponding power spectrum at frequency f.

In one example, the wideband quasi-zero phase filter operator is:

wherein f is_nIs the folding frequency obtained from the data sampling rate.

In one example, computing the broadband processed seismic data from the f-x domain shot gather data and the broadband quasi-zero phase filter operator comprises:

and calculating the product of the broadband quasi-zero phase filtering operator and shot gather data of the f-x domain in the frequency domain, and performing one-dimensional Fourier inversion on the calculation result to the t-x domain to obtain the seismic data after broadband processing.

Specifically, if the seismic wavelet is used as an input for the inverse filtering, the desired output is d (t) spike. The basic idea behind pulse deconvolution is to design an inverse wavelet a (t) operator that transforms a known input seismic signal into a given desired output spike signal, which is pulse deconvolution. The basic equation for pulse deconvolution is:

in general, the seismic wavelet is unknown, and in order to find the inverse filter factor in the case of unknown wavelet, certain assumption conditions must be added to the seismic wavelet and the reflection coefficient sequence, including:

1. let the reflection coefficient sequence r (t) be a random white noise sequence, i.e. its autocorrelation is:

2. the seismic wavelets are assumed to be of minimum phase.

According to hypothesis 1, the autocorrelation R of seismic wavelets_rrAutocorrelation R that can be recorded with seismographs_xxInstead of this. According to the hypothesis 2, it is known that zeros of the Z transform b (Z) of the seismic wavelet are all outside the unit circle, that is, zeros of the denominator polynomial of the Z transform a (Z) 1/b (Z) of the inverse filter factor a (t) are all outside the unit circle, so a (t) is stable and physically realizable. The equation for this pulse deconvolution becomes:

this is the fundamental equation of pulse deconvolution, where each element in the coefficient matrix can be directly derived from the seismic record. After the deconvolution factor a (t) is found, it is convolved with the seismic record x (t), that is:

s(t)＝a(t)*x(t) (13)

s (t) is the new seismic record output after pulse deconvolution. As mentioned above, the pulse deconvolution compresses the seismic wavelet to improve the resolution, and simultaneously enhances the noise energy of the low-frequency end and the high-frequency part except the effective signal to reduce the signal-to-noise ratio of the effective reflected signal.

Inputting seismic shot gather data loaded with an observation system information channel head; on the input shot gather record, the time-sharing window calculates the seismic wavelet of each data local time window in the shot gather by utilizing the autocorrelation function; aiming at the seismic wavelets of each time window of each channel of data in the shot gather, a pulse deconvolution method is utilized, the expected output is a sharp pulse, a deconvolution wavelet operator is obtained by solving a DebBetz matrix, deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, one-dimensional Fourier transform is carried out on the operation result, and the shot gather data of the f-x domain is obtained.

In order to overcome the defect of pulse deconvolution, an attenuation function is designed in a frequency domain according to the frequency band range of an effective signal of seismic data, so that the attenuation function is in an effective frequency band f₁And f₂Is close to 1, and is in f₁And f₂The effective signal is attenuated by 6 decibels rapidly, the effective signal is attenuated to a minimum value rapidly outside an effective frequency band, noise signals of a low-frequency end and a high-frequency end of the effective signal are subjected to suppression attenuation processing, and the signal-to-noise ratio of the effective signal is not reduced while the resolution of the effective signal is improved. The effect is equivalent to designing a broadband zero-phase filter m (t) in the time domain and carrying out broadband quasi-zero-phase convolution operation on the seismic signals after pulse deconvolution processing.

The frequency bandwidth of the reflected signal is estimated by using an energy concentration method, and the basic principle is that aiming at the frequency spectrum of a Gaussian-shaped effective signal, a limited range is set after the central frequency of the signal is estimated in a self-adaptive mode, and then the signal energy is concentrated to the range to estimate the frequency bandwidth, namely f₁And f₂。

The Gaussian-shaped power spectrum is formula (4), and the amplitude is attenuated by 3 dB of variable quantity delta f according to the definition of a half-power point_3dbThe relationship between σ and σ is as follows:

△f_3db＝2.335σ (14)

for any normal distribution, mu is expected, the standard deviation is sigma, and values on each side of mu are expected to be 50%, with 68% between mu-sigma and mu + sigma and 95% between mu-2 sigma and mu +2 sigma. If the variable has n samples, μ and σ are defined by:

in the process of estimating the bandwidth, the average power spectrum of the shot record is firstly calculated as:

wherein E (f) is the amplitude of each channel with frequency f in the shot record, N is the channel number of each shot, and then an extreme point is calculated on the average power spectrum, and the frequency corresponding to the extreme point is the central frequency f₀The average of the next extreme points on both sides of the extreme point is defined as the deconvolution expected frequency f_m. Desired output dominant frequency f by deconvolution_mRespectively accumulating power spectrums from the center to two sides, setting the energy threshold value to be 70%, and automatically calculating f in each processing time window according to a formula (4)₁And f₂。

FIGS. 4a and 4b show a schematic of wavelets and spectra, respectively, of broadband processed seismic data according to one embodiment of the invention.

According to the frequency band range of the effective wave on the shot record, at the low frequency end f of the frequency domain₁And a high frequency end f₂The frequency domain attenuation operator is designed as equation (5). Calculating the product of the broadband quasi-zero phase filtering operator and shot gather data of the f-x domain in the frequency domain; and performing one-dimensional Fourier inversion on the calculation result to a t-x domain to obtain the seismic data after the broadband quasi-zero phase deconvolution, wherein wavelets and frequency spectrums of the seismic data are respectively shown in fig. 4a and 4 b.

Example 1

To facilitate understanding of the solution of the embodiments of the present invention and the effects thereof, a specific application example is given below. It will be understood by those skilled in the art that this example is merely for the purpose of facilitating an understanding of the present invention and that any specific details thereof are not intended to limit the invention in any way.

Fig. 5a, 5b, and 5c are schematic diagrams illustrating comparison of the monochroom recordings before processing, after de-ghost only, and after de-ghost and wideband quasi-zero phase deconvolution processing, respectively, according to an embodiment of the present invention.

For example, as shown in fig. 5a, 5b and 5c, the single shot records before and after being processed by the method of the present invention have good ghost wave suppression and good wave group characteristics.

Fig. 6a, 6b, and 6c are schematic diagrams showing PSTM profiles before processing, after depoisoning only, and after depoisoning and wideband quasi-zero phase deconvolution processing, respectively, according to an embodiment of the present invention.

The PSTM cross-sectional surfaces before and after being processed by the method of the invention are shown in FIG. 6a, FIG. 6b and FIG. 6c, so that the visible resolution is high and the wave group characteristics are good; the method is beneficial to explaining horizon tracking and seismic inversion, and the quality is obviously improved.

The invention provides a self-adaptive ghost wave removing and broadband quasi-zero phase deconvolution combined processing system, which is characterized by comprising: a memory storing computer-executable instructions; a processor executing computer executable instructions in the memory to perform the steps of:

inputting original seismic data;

carrying out self-adaptive ghost wave removing processing on the original seismic data to obtain seismic data after ghost wave suppression;

and carrying out broadband quasi-zero phase deconvolution processing on the seismic data subjected to ghost wave suppression to obtain the seismic data subjected to broadband processing.

In one example, adaptive de-ghosting processing includes:

performing one-dimensional Fourier transform on the original seismic data to obtain original frequency domain seismic data;

respectively calculating the delay time of the shot-geophone point relative to the primary reflected wave, and further calculating a frequency domain ghost wave removing operator;

calculating the product of the original frequency domain seismic data and the frequency domain ghost wave removing operator to obtain ghost wave removing seismic data;

and performing one-dimensional Fourier inversion on the seismic data without the ghost waves to a t-x domain to obtain the seismic data after ghost wave suppression.

In one example, the delay time of the shot point with respect to the primary reflection is calculated by equation (1):

the delay time of the detection point with respect to the primary reflected wave is calculated by equation (2):

wherein, t_sDelay time of shot point relative to primary reflection, t_rIs the delay time of the point of detection relative to the primary reflection, x is the offset of the seismic data trace, d_wIndicating depth of sea floor, d_sDenotes the depth of the excitation point, d_rRepresenting the depth of the wave detection point and v is the speed of the seawater.

In one example, the frequency domain de-ghost operator is computed by equation (3):

wherein A (f) is a frequency domain ghost wave removing operator, R₀For the sea surface reflection coefficient, ω is the angular frequency, ω ═ 2 π f, and i denotes the complex exponential.

In one example, the wideband quasi-zero phase deconvolution process includes:

aiming at the seismic wavelets of each time window of each channel of data in the shot gather of the seismic data after ghost wave suppression, pulse deconvolution operation in each time window of each channel is completed in the t-x domain shot gather, and then shot gather data of f-x domain is obtained;

calculating the frequency bandwidth of the reflected signal by an energy concentration method;

determining a broadband quasi-zero phase filter operator according to the frequency bandwidth;

and calculating the seismic data after the broadband processing according to the shot gather data of the f-x domain and the broadband quasi-zero phase filtering operator.

In one example, calculating the frequency bandwidth of the reflected signal by energy concentration comprises:

calculating the average power spectrum of the shot record, and further calculating an extreme point on the average power spectrum;

determining the deconvolution expected frequency by taking the frequency corresponding to the extreme point as a central frequency;

and accumulating the power spectrums towards two sides by taking the deconvolution expected frequency as a center, setting an energy threshold value and calculating the bandwidth.

In one example, the low-frequency end frequency and the high-frequency end frequency of the effective signal of each time window are calculated by formula (4), and the frequency bandwidth is determined:

wherein σ is the mean square error of the frequency of the Gaussian-shaped power spectrum, f₀For the center frequency, P (0) is the corresponding power spectrum at zero frequency, and P (f) is the corresponding power spectrum at frequency f.

In one example, the wideband quasi-zero phase filter operator is:

wherein f is_nIs the folding frequency obtained from the data sampling rate.

In one example, computing the broadband processed seismic data from the f-x domain shot gather data and the broadband quasi-zero phase filter operator comprises:

and calculating the product of the broadband quasi-zero phase filtering operator and shot gather data of the f-x domain in the frequency domain, and performing one-dimensional Fourier inversion on the calculation result to the t-x domain to obtain the seismic data after broadband processing.

It will be appreciated by persons skilled in the art that the above description of embodiments of the invention is intended only to illustrate the benefits of embodiments of the invention and is not intended to limit embodiments of the invention to any examples given.

Having described embodiments of the present invention, the foregoing description is intended to be exemplary, not exhaustive, and not limited to the embodiments disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments.