CN114137606A

CN114137606A - Stable spectrum simulation deconvolution method

Info

Publication number: CN114137606A
Application number: CN202111324788.8A
Authority: CN
Inventors: 钱峰; 刘章波; 范昱琪; 胡光岷
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2021-11-10
Filing date: 2021-11-10
Publication date: 2022-03-04

Abstract

The invention discloses a robust spectrum simulation deconvolution method, which is applied to the field of seismic physical exploration and aims at the problem of uncontrollable and unstable traditional spectrum simulation deconvolution.

Description

Stable spectrum simulation deconvolution method

Technical Field

The invention belongs to the field of seismic physical exploration, and particularly relates to a seismic data processing technology.

Background

The natural gas and the petroleum are industrial products which rely on the existence of the society nowadays, relate to various aspects of the national civilization and the living and production, and have self-evident importance. Today global development requires more and more energy, and it is more and more difficult to obtain available underground reservoirs under existing conditions. The accurate acquisition of their distribution in the ground is an important position in energy exploitation.

Seismic physical prospecting is a common method for obtaining information about the specific distribution of energy sources, such as oil and gas, and geological structures, and is one of the most effective and important methods for studying and ascertaining geological structures. The method obtains geological formation information from the elasticity of the rock. FIG. 1 simply depicts the main working principle of seismic physical prospecting: the geophysical prospecting personnel can generate a reflection phenomenon when an artificial seismic source excited on the ground meets an area boundary surface formed between different geological layers in the underground propagation process, and the earth surface detector can obtain seismic data of the area by receiving reflected waves transmitted back to the ground. Since the seismic records transmitted back to the surface are filtered through different earth layers, their characteristics are related to the structure of the properties of the different rock layers in the subsurface. The explorationist can obtain the information such as wave velocity, frequency and energy by processing and analyzing the collected seismic records, so as to deduce the distribution of resources such as underground mineral oil gas and the like, and carry out the next energy exploitation.

Generally, the resolution of original seismic records directly acquired through seismic exploration is generally low, and in order to acquire more geological information, the seismic data needs to be processed in various ways, so that the effect of improving the resolution of seismic data is achieved. Through long-term research, many techniques for improving resolution, such as denoising, inversion, static correction, deconvolution, superposition, offset imaging, etc., have been proposed to increase the signal-to-noise ratio and improve resolution to some extent. Because the resolution ratio of the seismic signal has important guiding significance in lithofacies and geological horizon analysis, the improvement of the resolution ratio of seismic data and the acquisition of reflection coefficients occupy important positions in the field of seismic data processing. Generally speaking, the essence of the high-resolution seismic exploration technology is to recover the weak effective information in the seismic data, so as to achieve the purpose of obtaining the detailed geological structure distribution in a complex area. The method has important value for the development of energy exploration business at the present stage of China.

The convolution model is a common model in the field of seismic exploration. It shows that reflection coefficient is one of the important components constituting seismic data. Since reflection coefficients contain a large amount of information about the geological structure, obtaining a reflection coefficient from seismic data is of great significance in oil exploration work. Deconvolution techniques are commonly used in the field of seismic data processing to improve the resolution of seismic data. The technology can extract the reflection coefficient sequence and enhance the high-frequency information of the seismic data, and finally the time resolution of the seismic data is improved.

The spectral analog deconvolution method is a conventional and common approach. Unlike other methods, it no longer makes assumptions about the form of the distribution of the reflection coefficients, which makes it effective to reduce the influence of the reflection coefficients on the wavelet estimate. In addition, the method assumes that the seismic sub-spectrum is smooth, and utilizes the amplitude spectrum of the seismic record to directly fit a wavelet as the seismic sub-spectrum under the least square criterion, so as to widen the frequency band and effectively improve the seismic data resolution.

In the traditional spectrum simulation deconvolution method, the polynomial order of the fitting wavelet needs to be given before the experiment is carried out, and the accuracy of the selected parameters needs to be tested through a large number of experiments because the fitting result cannot be judged in advance. The selected order is too small, so that under-fitting easily occurs, and the order is too large, so that over-fitting easily occurs due to noise interference. In addition, because the spectral simulation deconvolution is to fit a wavelet from the single-channel seismic data amplitude spectrum during processing, but because the data of different gathers in the same region are affected by factors such as receiving conditions and the like, the fitting parameter value obtained according to a certain piece of data is not necessarily suitable for the whole data, and errors are caused to the deconvolution result at this time. In order to solve the two problems, it is necessary to research a spectral simulation deconvolution method which is convenient to operate and has a wide application range.

Disclosure of Invention

To solve the above technical problem, the present invention provides a robust spectral analog deconvolution method,

the technical scheme adopted by the invention is as follows: a robust spectral modeling deconvolution method, assuming that the spectral waveform of a seismic wavelet is the same as the rake wavelet waveform, on the basis that the seismic wavelet amplitude spectrum is smooth, comprising the steps of:

s1, solving a mean trace gather m of the seismic data D;

s2, solving autocorrelation of m, and judging the dominant frequency f of the wavelet_m1；

S3, generating the main frequency f_m1、f_m1+f_increaseThe combined wavelet of (1); f. of_increaseRepresenting the degree of frequency broadening;

s4 solving cutoff frequency f of dominant frequency combined wavelet_c；

S5, according to the combined wavelets and the dominant frequency f_m1The Rake wavelet calculation operator O (f) and at f_cProcessing it robustly;

s6, acting the operator after the steady processing on D to obtain a deconvolution result D_r。

The specific implementation process of step S2 includes the following sub-steps:

s21, calculating the autocorrelation of the seismic data in the mean trace gather m and the Rake wavelets with different dominant frequencies;

s22, calculating the time interval of the first pair of zero-crossing points of the autocorrelation relative to the central symmetry point;

s23, comparing the time interval obtained by autocorrelation calculation of the seismic data with the time interval obtained by autocorrelation calculation of the Rake wavelets with different dominant frequencies;

and S24, taking the dominant frequency of the Rake wavelet with the minimum difference with the seismic data autocorrelation time interval as the dominant frequency of the seismic wavelet corresponding to the seismic data.

The expression of the combined wavelet described in S3 is:

step S5 specifically includes: at a cut-off frequency f_cIn each frequency interval with a certain length, curve fitting is carried out on the seismic data amplitude spectrum, and the fitted curve in the interval is used for replacing the Rake wavelet in the same interval to obtain a new wavelet R_new(f) Obtaining an updated deconvolution operator O using the following equation_new(f)

The invention has the beneficial effects that: under the assumption that the seismic wavelet waveform is the same as the Rake wavelet, the invention estimates the dominant frequency of the seismic wavelet by utilizing seismic data autocorrelation, and then provides a zero-phase wavelet as expected output on the basis of the Rake wavelet so as to quantitatively and controllably improve the spectrum width of the seismic data; finally, calculating a deconvolution operator and acting the operator on the seismic data to realize spectrum simulation deconvolution with obvious and steady effect; the method of the invention realizes the prediction of high-frequency information by using the medium and low frequency information of the seismic data so as to reduce the amplification of high-frequency noise and realize the stable prediction of the high-frequency information of the seismic data.

Drawings

FIG. 1 is a schematic diagram of the working principle of seismic physical prospecting;

FIG. 2 is a schematic diagram of time domain waveforms of different phase sub-waves according to an embodiment of the present invention;

wherein, (a) is a zero-phase wavelet; (b) is the maximum phase wavelet; (c) is the minimum phase wavelet; (d) is a mixed phase wavelet;

FIG. 3 is a deconvolution filtering process provided by an embodiment of the present invention;

FIG. 4 is a basic flow chart of deconvolution provided by an embodiment of the present invention;

FIG. 5 is a graph of a simulated fit of spectra provided by an embodiment of the present invention;

FIG. 6 is a flow chart of a robust deconvolution method based on spectral modeling according to an embodiment of the present invention;

FIG. 7 is a Rake wavelet form provided by an embodiment of the present invention;

wherein, (a) is a time domain diagram of the Rake wavelet; (b) is the amplitude spectrum of the Rake wavelet;

FIG. 8 is a comparison graph of time domain waveforms of Rake wavelets with different dominant frequencies according to an embodiment of the present invention;

FIG. 9 is a graph comparing the amplitude spectra of Rake waves with different main frequencies according to an embodiment of the present invention;

FIG. 10 is a frequency domain waveform of a combined wavelet provided by an embodiment of the present invention;

FIG. 11 is a time domain waveform diagram of a combined wavelet provided by an embodiment of the present invention;

FIG. 12 is a graph comparing the effects of autocorrelation length provided by embodiments of the present invention;

FIG. 13 is a comparison of waveforms provided by embodiments of the present invention;

FIG. 14 is a comparison of waveforms provided by embodiments of the present invention;

FIG. 15 is a sequence of reflection coefficients provided by an embodiment of the present invention;

FIG. 16 is a frequency domain comparison graph of single-pass data of deconvolution result of synthetic data;

FIG. 17 is a time domain comparison graph of a single trace of data resulting from deconvolution of synthetic data;

FIG. 18 is a time domain comparison graph of a single trace of data of an actual data deconvolution result;

FIG. 19 is a comparison graph of the results of the 3 deconvolution methods provided by the present invention;

wherein, (a) is original data, (b) is a deconvolution result obtained by the method of the invention, (c) is a traditional spectrum simulation deconvolution result, and (d) is a pulse deconvolution result;

FIG. 20 is a comparison graph of local actual data results of 3 deconvolution methods provided by an embodiment of the present invention;

FIG. 21 is a comparison graph of the local effect of the deconvolution method of commercial software according to the present invention;

wherein, (a) is original data, (b) is a deconvolution result obtained by the method of the invention, (c) is a GeoEast system surface consistency deconvolution effect, and (d) is an Omega deconvolution result;

FIG. 22 is a graph comparing spectra of a commercial software deconvolution method according to an embodiment of the present invention.

Detailed Description

In order to facilitate the understanding of the technical contents of the present invention by those skilled in the art, the present invention will be further explained with reference to the accompanying drawings.

The method of the invention involves the following techniques:

1. seismic record convolution model

Seismic wavelets are a signal with a definite starting time, a finite duration and a finite energy, and are the most basic units of composition in seismic records. Generally, it can be divided into non-time-varying wavelets and time-varying wavelets. The waveform of the time-invariant wavelet does not change with time, and the morphology is determined by the type of the seismic source and the rock properties of the wavelet. The time-varying wavelet shape changes over time under the influence of a number of factors, such as formation absorption. If the difference between the stratums is more than half of the duration of the reflected wave, the geological stratum generates stratum wavelets containing unique lithology information of the stratum wavelets. The explorationist can obtain useful information about lithology explanation from wavelets generated by different strata, and judge the distribution condition of underground strata.

In the actual seismic data acquisition, because the lithology of the stratum has higher complexity, how to acquire the seismic wavelets also becomes a complicated problem. Seismic wavelets can generally be considered to be a function consisting of a certain frequency, period, attenuation factor and initial phase:

b(t)＝cp(t)sin(α+ωt) (1)

in formula (1), b (t) represents seismic wavelets, t represents time, c represents amplitude, p (t) represents a specific attenuation function, alpha represents the initial phase of the wavelets, omega represents frequency, and t represents time.

P (t) in the formula (2) may be set in various forms according to the degree of attenuation. E.g. linear piecewise function, gaussian function, exponential function p (t) ═ e^-λt(λ>0) And the like. Alpha depends on the phase delay information of the determined wavelet.

The wavelets can be classified into zero phase, maximum phase, minimum phase, and mixed phase, etc. according to the phase delay property. FIG. 2 shows the time domain waveforms corresponding to four different phase wavelets.

It can be seen from fig. 2 that the phase properties of the wavelets exhibit their delay properties. The zero-phase wavelet is symmetric about the time 0 point. According to the difference of the distribution positions of the signal energy in the time domain, the wavelet phase information of which the energy is concentrated in the front part, the middle part and the rear part corresponds to a minimum phase, a mixed phase and a maximum phase respectively.

Among the seismic wavelets, the zero-phase wavelet has the strongest resolution and is most beneficial to seismic interpretation. The advantages of the zero-phase wavelet are mainly embodied in the following 3 aspects:

(1) the narrower the wavelet shape, the higher the resolution, and the smaller the side lobe of the zero-phase wavelet than other phase wavelets with the same bandwidth, which shows that the energy can be concentrated in a short time, so that the high resolution is realized.

(2) Because the pulse reflection time of the zero-phase wavelet corresponds to the maximum value of the waveform, the time can be calculated more accurately, and the subsequent seismic data interpretation is facilitated.

(3) The zero-phase wavelet can solve the problem of thin layer distinguishing more favorably. For example, a zero-phase wavelet may accurately distinguish two 16ms apart spikes, which is not the capability of wavelets of several other phases.

The seismic convolution model is currently one of two major geophysical models. Convolution models are less computationally intensive and do not require high quality seismic data than other wave equation models, and are therefore heavily used in current practical seismic data processing. The convolution model illustrates the manner in which synthetic seismic data is generated, and to further understand this model, we need to understand how the actual seismic records are formed. In fig. 1, how a seismic physical prospecting method obtains a seismic record is shown, a sharp pulse excited by a seismic source is subjected to filtering action of each geological stratum in the process of underground propagation, and high-frequency components, amplitude and energy are gradually reduced, and finally, the sharp pulse is changed into a seismic wavelet with a certain continuous length.

The proposal of the convolution model lays a solid foundation for the formation of a plurality of important seismic exploration technologies, such as wave impedance inversion, deconvolution and the like. The formula of the convolution model is shown as follows:

d(t)＝b(t)*x(t) (2)

wherein d (t) is seismic record, b (t) is seismic wavelet, x (t) is reflection coefficient, and x represents convolution operation. Equation (2) is a model under ideal conditions, and there is usually some noise influence in the actual seismic trace recording, and now equation (2) is modified to obtain the following equation:

d(t)＝b(t)*x(t)+n(t)＝s(t)+n(t) (3)

where n (t) represents a noise term, s (t) x (t) n (t).

Deconvolution is a seismic signal processing method which is established on the model and improves the resolution of seismic data by extracting reflection coefficients. The reflection coefficient contains information such as rock attributes of all underground geological layers, and therefore, the reflection coefficient extraction has important practical significance in seismic data interpretation. The wave impedance of a rock reflects the ability of a stress wave to penetrate and reflect in the rock and determines the reflection coefficients of different rock layers. The positive and negative values of the reflection coefficient can show the density property of the upper rock layer and the lower rock layer of the reflection surface. The number from sparse to dense is positive, and the number from dense to sparse is negative. Therefore, the prospecting personnel can extract the reflection coefficient from the seismic record through the deconvolution technology and analyze the reflection coefficient, so as to obtain the information which can effectively explain the stratum distribution.

2. Deconvolution enhanced resolution

In order to obtain richer stratum information, the improvement of the resolution of seismic data is a target pursued by seismic prospecting personnel in the seismic data processing. In field acquisition, seismic data with higher resolution cannot be directly acquired. One important reason for this is that the difference between the time of the different reflection layers back to the earth surface is very small, usually tens of milliseconds or even several milliseconds, which causes the interference and overlapping between different waveforms when observing the acquired seismic records, and is difficult to distinguish, so the resolution of the acquired seismic data is very low. If the reflected wave can be expressed as a sharp pulse in the seismic record, the resolving and interpreting capabilities of the seismic record can be effectively improved.

For equation (2), if one can solve for an operator w (t) is exactly the inverse of b (t), i.e.

w(t)＝b^-1(t) (4)

Convolution of w (t) and s (t) to obtain

The reflection coefficient is thus obtained, which is called deconvolution, since this is the inverse of the convolution operation. To improve resolution, deconvolution requires compressing b (t), which has a certain run length, into a sharp pulse form.

The process of deconvolution filtering is shown in fig. 3, and the basic flow of deconvolution is shown in fig. 4.

3. Principle of spectrum simulation deconvolution

Although the above mentioned deconvolution methods are also common methods today, they all have a limitation that the influence of the reflection coefficient cannot be avoided in the solving process, for example, the least square deconvolution cannot distinguish the wavelet from the reflection coefficient in the calculation. In the actual seismic data processing, the spectrum simulation deconvolution is a method capable of reducing the influence of reflection coefficients.

Rosa et al were inspired by Ricker research, and it is believed that wavelet amplitude spectra could be obtained by fitting from seismic traces with mathematical tools on the premise that the seismic wavelet amplitude spectra were smooth.

The hypothesis is first analyzed for plausibility. Generally, the seismic wavelet excited by an explosive source on land is in the form of a spike pulse, and when the seismic wavelet is transmitted underground, high-frequency information is more lost due to absorption and reflection of each rock stratum, so that the rock stratum can be considered to cause loss with low-pass characteristic to the seismic wave. Receivers that receive seismic records at the surface are generally considered to be high-pass in frequency response. While the geophone's coupled response to the ground is considered a bandpass response because it imposes a limit on the high frequency information of the seismic signal. After the seismic wavelet passes through the total seismic response formed by connecting different types of band-pass responses in series and passes through a series of data preprocessing links, the amplitude spectrum of the seismic wavelet at the moment can be considered to have the characteristics of single peak and smoothness.

As verified by extensive experimental development, the mathematical model selected by Rosa et al is as follows:

in the formula

Representing the fitting of wavelet amplitude spectra, α is a constant, and H (f) is a polynomial of order β of frequency f. Spectral modeling deconvolution is performed by fitting the seismic amplitude spectrum with equation (6) under the least squares error criterion to solve for the wavelet amplitude spectrum. In fact, the wavelet spectrum fitting result obtained by solving the formula (6) is a unimodal smooth curve similar to a Rake wavelet.

When alpha and beta are determined, coefficients for each of H (f) are calculated and determined

The specific expression of (1). The traditional spectrum simulation deconvolution is obtained by calculation and solution

Then, can be based on

And designing a deconvolution operator. The deconvolution operator h (f) can be obtained by calculating equation (7):

where λ represents a penalty factor, 0< λ <1, and its value decreases as the effective band center distance d increases.

And multiplying H (f) by the amplitude spectrum of the seismic record under the condition of keeping the phase spectrum unchanged, and obtaining a deconvolution output result of a time domain through inverse Fourier change. Since this process does not change the phase information of the seismic data, spectral modeling deconvolution is also considered a zero-phase deconvolution method.

The spectrum simulation deconvolution can effectively widen the amplitude spectrum of seismic data, reduce the influence of reflection coefficients in the wavelet amplitude spectrum fitting process, and finally realize the improvement of seismic data resolution.

Based on the technology, the method firstly introduces the generation of the confrontation network to better learn and represent the statistical distribution and the potential characteristics of the data so as to ensure the completeness, the accuracy and the like of the solution, and then introduces the seismic data into the network training to solve the serious dependence of the inversion result on the logging data. The method specifically comprises the following steps:

1. zero phase wavelet based on Rake wavelet

The Rake wavelet is a zero-phase wavelet frequently used in forward modeling and seismic record synthesis in the field of seismic exploration. The waveform of the wave form is simple, two sides of the only main peak are respectively provided with a harmonic side lobe, and the convergence time period quickly approaches to zero. Formula (8) is with the dominant frequency f₀The expression of the Rake wavelet in the time domain shows that the Rake wavelet has the characteristic of exponential decay. The shape of the waveform of the Rake wavelet is determined with the choice of dominant frequencies.

r(t)＝[1-(2πf₀t)²]exp[-(πf₀t)²] (8)

Fourier transform is carried out on the formula to obtain an amplitude spectrum expression form of the Rake wavelet:

dominant frequency f₀The waveform and amplitude spectrum of a 20Hz rake wavelet are shown in fig. 7(a) and (b), respectively. From fig. 7(a), it is seen that the rake wavelet is approximated to a sharp pulse signal in the time domain, and the signal is concentrated in a narrow time range. Equation (8) shows that the Rake wavelet has only time t as an argument, and accordingly the positive semi-axis wavelet side lobe in the time domain gradually attenuates with time. As can be seen from FIG. 7(b)The amplitude spectrum of the Rake wave has only one main peak, and the peak frequency and the main frequency f₀And correspondingly.

FIG. 8 shows a group of signals using a main frequency f₀Different Rake wavelet time domain waveform comparison results. As can be seen from fig. 8, the energy of the rake wavelet is strongest at the initial 0 time. The maximum values of the wavelets of different dominant frequencies in the negative amplitude direction are equal and are-0.45. Dominant frequency f of wavelets₀The larger the amplitude, the more concentrated the energy and the smaller the harmonic side lobe. And with the continuous reduction of the main frequency, the attenuation of the wavelet along the time direction gradually becomes slow, and finally the energy approaches to 0 on a negative half shaft.

Since the Rake wavelet is a zero phase wavelet, the phase of the wavelet is always equal to zero, so that the frequency domain formula of formula (9) can be directly used to express the amplitude spectrum of the wavelet. Fig. 9 corresponds to an amplitude spectrum comparison chart of the neutron wave in fig. 8. As is clear from fig. 9, the rake wavelet gradually decreases in peak value of the frequency spectrum with an increase in dominant frequency, while the wavelet amplitude spectrum gradually becomes gradually gentle and wider. f. of₀The smaller the frequency spectrum, the more concentrated the frequency spectrum in the narrower frequency band, and the more energy is dispersed along with the continuous increase of the main frequency.

On the basis of the Rake wavelets, the invention provides a novel zero-phase combined wavelet. The combined wavelet has an analytic form in the frequency domain, and a specific function expression is as follows, wherein f₀<f₁：

From equation (11), it can be seen that the spectrum of the combined wavelets is composed of two dominant frequencies, f₀And f₁The amplitude spectrum of the Rake wavelet is normalized. Frequency from 0 to f₀In time between, the image of S (f) is the main frequency f₀The waveform of the Rake wavelet after normalization in the corresponding frequency interval; at a frequency of f₁The latter part being the main frequency f₁The waveform of the Rake wavelet after normalization in the corresponding interval; at f₀And f₁The section in between is 1.

The comparison of the combined wavelet and two Rake wavelets of different dominant frequencies in the frequency domain is shown in FIG. 10. As can be seen from fig. 10, the combined wavelet is obtained by combining two rake wavelets with different dominant frequencies, and the low-frequency portion of the rake wavelet with a lower dominant frequency and the high-frequency portion of the rake wavelet with a higher dominant frequency are reserved, so that the bandwidth of the combined wavelet in the frequency domain is higher than that of the original two rake wavelets, and more low-frequency and high-frequency information can be contained. Meanwhile, the combined wavelet has a definite expression in a frequency domain and a simple form, and is beneficial to subsequent calculation.

Fig. 11 shows a comparison graph of the combined wavelet and two rake wavelets with different dominant frequencies normalized in the time domain, and it can be seen from fig. 11 that the side lobe of the combined wavelet proposed by the present invention in the time domain is smaller than two rake wavelets, and the energy of the combined wavelet is more concentrated than that of the other two wavelets.

For a discrete trace of seismic data Y_n(N-1, 2, …, N), and Y_nNot equal to 0, order

V_yI.e. signal y_nNormalized Variance Norm (Normalized) capable of characterizing signal y_nThe degree of energy concentration.

The normalized variance mode has 2 basic properties as follows:

properties 1:

and when V_y＝1，Y_nIs a sharp pulse, i.e.

When V is_yWhen 1/N is equal to Y_nThe amplitude value at each instant being the same, i.e. | y₁|＝|y₂|＝…＝|y_n|。

Properties 2: for a discrete spike A_n，

Wherein I is 1,2, …, I; n is more than or equal to 1₁<n₂<…<n_I≤N₁Let the minimum time interval between adjacent spikes be d min { n ═ n_i-n _i-12, …; let B_n＝(b₁,b₂,…,b_N2),N₂D-1 or less, making Y_n＝A_n*B_nThen, there are:

V_y＝V_a*V_b (15)

in the following table 1, the dominant frequencies f with the same data length are calculated respectively₀The normalized variance mode of the different Ricker wavelets and the combined wavelets, as can be seen from Table 1, the dominant frequency parameter proposed by the present invention is f₀,f₁The normalized variance mode of the combined wavelet is better than the dominant frequency of f₀And f₁The standard variance mode of the Rake wavelets shows that under the condition of the same dominant frequency parameter, the energy concentration degree of the combined wavelet provided by the invention on the time domain is better than that of the Rake wavelets.

TABLE 1 different norm Fang-Difference-mode calculation results

2. Robust spectrum simulation deconvolution algorithm based on Rake wavelets

In order to solve the problem of uncontrollable wavelet fitting of spectral simulation deconvolution and improve the deconvolution effect, the invention provides a robust spectral simulation deconvolution method based on a Rake wavelet on the basis of the conventional spectral simulation deconvolution method. The method assumes that the spectrum waveform of the seismic wavelet is the same as the waveform of the Rake wavelet on the basis that the amplitude spectrum of the seismic wavelet is smooth. Thus, conventional spectral modeling deconvolution canesObtained by fitting the amplitude spectrum of the line wavelet

Becoming a waveform of a determined type, i.e. dominant frequency f_m1Unknown Rake wavelets R (f).

The algorithm sets the desired output to the dominant frequency f_m1、f_m2The combined wavelet S (f) of (a) can obtain a frequency expansion operator O (f) by using the following formula to realize that the frequency band expansion of R (f) is changed into S (f)

The operator can realize the frequency band promotion of the quantification of the Rake wavelet, so that f_increase＝f_m2-f_m1I.e. the extent of the bandwidth broadening achievable by the operator. Considering the difference in detail between the morphology of the rake wavelet and the actual seismic data, the boundary robust processing will be performed on o (f), and this specific processing will be described in detail later.

For a seismic data, solving the amplitude spectrum D (f) of the seismic data, and then acting the operator O (f) on the amplitude spectrum D (f), the frequency band broadening of the seismic data quantification can be realized, and the effect of improving the data resolution is achieved.

D_r(f)＝O(f)·D(f) (17)

D in the above formula (17)_r(f) Namely, the amplitude spectrum of the deconvolution result is obtained, and then the seismic data with better resolution improvement can be obtained by performing inverse Fourier transform on the amplitude spectrum.

As shown in fig. 6, the implementation process of the present invention first preprocesses data, generally, performs denoising, and the specific implementation process is shown as algorithm 1:

the algorithm provided by the invention adds the assumption that the amplitude spectrum of the seismic wavelet is the amplitude spectrum of the zero-phase Rake wavelet on the basis of the assumption of traditional spectrum simulation deconvolution, so that the main frequency of the wavelet is only required to be determined when the seismic wavelet is solved. The present invention will describe how to solve for wavelet dominant frequencies by computing the autocorrelation of the seismic data.

The present invention has demonstrated that the autocorrelation of both the seismic record and the seismic wavelet differ in value by only a constant factor. However, if the two are considered to be equal, the influence of the reflection coefficient sequence on the wavelet estimation cannot be avoided, and an inaccurate wavelet estimation result is obtained.

This problem is explained in detail by means of fig. 12. Fig. 12 shows the result of solving the open form of the amplitude spectrum after performing different length truncation on the autocorrelation of a seismic signal with a length of 300 points, wherein the first row of data retains the autocorrelation information of the data, the second row of data retains 1/2 information by truncation of the autocorrelation, and the third row of data retains 1/4 autocorrelation information. As can be seen from fig. 12, as the autocorrelation information is kept less and less, the waveform obtained by squaring the autocorrelation amplitude spectrum also gradually approaches the spectral shape of the wavelet. Illustrating that with this truncation, there is a diminishing decrease in information about the reflection coefficient in the autocorrelation of the seismic signal.

The specific way for solving the wavelet dominant frequency is as follows:

(1) calculating autocorrelation of the seismic data and the Rake wavelets with different dominant frequencies;

(2) calculating a time interval for a first pair of zero crossings of the autocorrelation about the centrosymmetric point;

(3) comparing the time interval obtained by autocorrelation calculation of the seismic data with the time interval obtained by autocorrelation calculation of the Rake wavelets with different dominant frequencies;

(4) and the dominant frequency of the Rake wavelet with the minimum difference with the autocorrelation time interval of the seismic data is the dominant frequency of the seismic wavelet corresponding to the seismic data.

According to the method, on the premise that the seismic wavelet form is set, the wavelet dominant frequency is solved by utilizing the autocorrelation characteristics, and as only the main autocorrelation information is utilized during solving, compared with the case that the autocorrelation of the seismic data is directly used as the autocorrelation of the seismic wavelets, the influence of the reflection coefficient on wavelet extraction can be reduced to a great extent. Meanwhile, compared with wavelet amplitude spectrum fitting of traditional spectrum simulation deconvolution, the method provided by the invention can obtain the seismic wavelets with known forms (Rake wavelets), and the fitting is not carried out like the traditional spectrum simulation method by continuously adjusting parameters.

Due to the detail difference between the Rake wavelets and the actual seismic data, if the operator O (f) is directly applied to the frequency spectrum of the seismic data, or the data is caused to be at the wavelet cut-off frequency f_cA nearby abnormal overflow. Therefore, the seismic data is combined with the operator O (f) at f_cThe vicinity is subjected to a certain correction process.

The algorithm is chosen at the cut-off frequency f_cIn each frequency interval with a certain length, curve fitting is carried out on the seismic data amplitude spectrum, and the fitted curve in the interval is used for replacing the Rake wavelet in the same interval to obtain a new wavelet R_new(f) Obtaining an updated deconvolution operator using the following equation

Operator O obtained after boundary robust processing_new(f) The method can ensure that the result of acting on the seismic data amplitude spectrum does not generate abnormal overflow near the cut-off frequency. The results of the curve fitting are shown in FIG. 13 below, including: amplitude spectra of seismic data mean gathers; dominant frequency f_m1(ii) is considered a seismic wavelet; fitting results of the seismic data; dominant frequency f_m1、f_m2The combined wavelet s (f), i.e., the desired output.

FIG. 14 shows the operator O to be solved_new(f) The results of the amplitude spectrum d (f) of the seismic data are applied. As can be seen from fig. 14, the deconvolution operator after robust processing can achieve quantitative improvement on the high-frequency information of the seismic data, and it is ensured that no abnormal overflow occurs at the medium-high frequency.

In order to verify the correctness of the algorithm provided by the invention, the artificially synthesized data is obtained by utilizing the convolution generation of the Rake wavelets and the random sequence. The main frequency of the Rake wavelet is set to be 20Hz, the length of the reflection coefficient sequence is set to be 1600ms, and the sampling interval is set to be 4 ms. Convolution is carried out on the reflection coefficient and the wavelets to generate 30 channels of data in total, and Gaussian white noise with the signal-to-noise ratio of 30dB is added to form synthetic data. The reflection coefficient sequence selected for the experiment is shown in fig. 15, and is a non-white noise sequence.

For the experiment of synthetic data, the extension degree of the algorithm is set to be 10Hz, parameters alpha and beta of the spectrum simulation deconvolution are determined to be 0.2 and 5 respectively through a certain experiment, the white noise coefficient of the pulse deconvolution is 0.001, and the operator length is 200 ms.

FIG. 16 is a frequency domain amplitude spectrum of the 10 th trace of data from several deconvolution experiments. As can be seen from fig. 16, compared with other methods, the spectral analog deconvolution method provided by the present invention can significantly improve the high frequency information while maintaining the low frequency information of the original data, and the resulting spectrum shape is closer to the original data, instead of changing the spectrum shape into a spectrum similar to white gaussian noise, which excessively generates high frequency noise.

FIG. 17 is a time domain waveform of the data of trace 10 from the experimental results of several deconvolution methods. As can be seen from fig. 17, the time domain results obtained by the method of the present invention contain more information about the reflection coefficient than the original data and the conventional spectral modeling method.

The results of the overall data obtained by several deconvolution algorithms are compared in fig. 18, and it can be seen from the results that the pulse deconvolution method is not suitable for seismic data with a reflection coefficient that is not white, and excessive noise is generated in the solution results, which is not favorable for obtaining effective interpretation information. And the spectrum simulation deconvolution can have good processing effect on the data of the type. Compared with the traditional spectrum simulation deconvolution algorithm, the method provided by the invention obtains more detailed information and has more obvious effect improvement. In a comparison experiment of synthetic data, the method provided by the invention obtains a good resolution improvement effect.

The verification results prove that the method provided by the invention can effectively improve the resolution of the reflection coefficient non-white noise data. In order to show the application effect of the algorithm of the invention in actual seismic data, a single shot data collected in an actual work area is selected to carry out experiments of 3 deconvolution methods, and comparison is carried out. The data is processed by field static correction, surface wave suppression, linear suppression, near-offset high-energy noise suppression, outlier suppression, spherical diffusion compensation, earth surface consistency amplitude compensation and the like in advance, the sampling interval is 4ms, the number of time points is 1000 points, and the data is 288 in total.

For the experiment of actual data, the extension degree of the algorithm is set to be 6Hz, parameters alpha and beta of the spectrum simulation deconvolution are set to be 0.2 and 6 respectively, the white noise coefficient of the pulse deconvolution is 0.001, and the operator length is 200 ms.

FIG. 19 is a time domain waveform diagram of the 100 th data in the experimental results of 3 deconvolution methods. As seen from fig. 19, compared with the conventional spectrum simulation deconvolution method, the method of the present invention has a high cut-off frequency, can limit data of an ultrahigh frequency portion, can adapt to a situation of high ultrahigh frequency noise, and can significantly improve high frequency information of original data while maintaining low frequency information of the original data, so as to obtain more data details and improve time resolution of the data. In addition, the method can improve the high-frequency information through the frequency extension parameter f before the experiment is carried out_increaseThe quantitative control is carried out, and the effect on the deconvolution is controllable and predictable. Meanwhile, compared with the pulse deconvolution method, the method provided by the invention keeps the waveform consistency with the original data on the frequency spectrum, does not have the limit that the reflection coefficient is a white noise sequence similar to the pulse deconvolution, and also shows that the method has stronger adaptability to the non-white-noise seismic data.

FIG. 20 compares the results of the 3 deconvolution methods applied to the actual data as a whole. Compared with the original data, the 2 comparison algorithms improve the resolution of the data to a certain extent, but the method provided by the invention has better effect particularly at some details. Fig. 21 is a diagram obtained by enlarging the experimental result in the grid in fig. 20, and it can be seen from comparison of data of the position indicated by the arrow that the algorithm provided by the present invention can obtain more and more continuous horizon information, so that the resolution enhancement effect is better.

The experimental results of the actual data and the original data are combined, so that the method provided by the invention can effectively and steadily improve the high-frequency information of the seismic data, broaden the frequency spectrum of the seismic data, and finally realize the improvement of the seismic data resolution. Meanwhile, the method is not limited by the fact that the reflection coefficient is a white noise sequence, and has better adaptability to different types of seismic data. Compared with the traditional spectrum simulation deconvolution method, the method can realize the quantification of the data and the predictable resolution improvement without setting sensitive parameters, thereby greatly reducing the time and the cost.

Because the algorithm of the invention depends on the cooperative topic 'deconvolution operator intelligent estimation and robust deconvolution technical research' of the research center of geophysical corporation in east China oil and gas Co., Ltd in laboratories and China, the invention compares the operation result of the algorithm, the GeoEast system surface consistency deconvolution result and the robust deconvolution experimental effect in the commercial software Omega seismic data processing system on the GeoEast seismic interpretation software system of China oil and gas Co., Ltd.

The data selected by the experiment are acquired in a certain work area of Pakistan, the sampling interval is 2ms, the time length is 3001 points, and the number of the data line is 330-380. The algorithm and the comparison algorithm provided by the invention act on the pre-stack seismic data and then compare the deconvolution effect of the post-stack seismic data. The algorithm provided by the invention sets the extension frequency to be 10 Hz.

FIG. 22 below shows the deconvolution results for the actual data on the GeoEast system at line number 330. In fig. 23, (a) is a partial result of the raw data after superposition, (b) is a deconvolution result obtained by the present algorithm, (c) is a GeoEast system surface consistency deconvolution result, and (d) is a deconvolution result in the commercial software Omega. Comparing fig. 23, it can be seen that the three deconvolution methods all achieve the improvement of the seismic data resolution to a certain extent, but the algorithm provided by the present invention can obtain more detailed information, especially at the deeper layer of the data, as shown by the result in the circle in fig. 22, it can be seen that the robust deconvolution based on the rake wavelet can obtain more and clearer information of the event, and has a stronger ability of improving the data resolution. In the position indicated by the arrow in fig. 22, the algorithm proposed by the present invention can obtain more information of small horizons, so as to provide convenience for the subsequent seismic interpretation work.

The spectrum results for several sets of data as a whole are shown in fig. 22, and several sets of data are more clearly compared in decibel form with the amplitude spectrum chosen in fig. 22. It can be seen from fig. 22 that the algorithm of the invention can effectively broaden the spectrum of the seismic data on the basis of retaining low-frequency information in the data, enhance the high-frequency information of the seismic data, and finally realize the resolution improvement of the seismic data. Meanwhile, the frequency band lifting effect of the algorithm is equivalent to the Omega steady deconvolution effect.

It will be appreciated by those of ordinary skill in the art that the embodiments described herein are intended to assist the reader in understanding the principles of the invention and are to be construed as being without limitation to such specifically recited embodiments and examples. Various modifications and alterations to this invention will become apparent to those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the claims of the present invention.

Claims

1. A robust spectral modeling deconvolution method, wherein, on the basis that the seismic wavelet amplitude spectrum is smooth, the seismic wavelet is assumed to have the same spectral waveform as the rake wavelet waveform, said method comprising the steps of:

s1, solving a mean trace gather m of the seismic data D;

s4 solving cutoff frequency f of dominant frequency combined wavelet_c；

2. The robust spectral modeling deconvolution method of claim 1, wherein the step S2 is implemented by the following sub-steps:

3. A robust spectral simulation deconvolution method according to claim 1, wherein the expression of the combined wavelet of S3 is:

4. the robust spectral modeling deconvolution method according to claim 1, wherein step S5 is specifically: at a cut-off frequency f_cIn each frequency interval with a certain length, curve fitting is carried out on the seismic data amplitude spectrum, and the fitting curve in the interval is used for replacing the Rake wavelets in the same interval to obtainA new wavelet R_new(f) Obtaining an updated deconvolution operator O using the following equation_new(f)