A method for processing seismic signals

Abstract

The present application relates to a kind of seismic signal processing method, belong to seismic interpretation processing technical field.Processing method includes: obtaining poststack data;The poststack data is converted to time-frequency, and the transformed spectrum is obtained;The spectrum of each seismic trace signal is segmented, and the segmented interval and spectral boundary point are obtained;By the interval and spectral boundary point after segmentation, an adaptive wavelet filter bank is constructed, and the original seismic signal is filtered, and different empirical mode components are obtained;Different empirical mode components are optimized, and the empirical mode component reflecting seismic wavelet information is obtained;The energy balance of the empirical mode component reflecting seismic wavelet information is carried out;The energy balanced empirical mode component and the non-optimized empirical mode component are reconstructed, and high-resolution seismic data body is obtained.The seismic data body obtained in the present application is well protected in low-frequency component, the frequency band width is expanded, has good fidelity, can satisfy the need of thin layer fine portrayal.

Description

TECHNICAL FIELD

[0001] The present application relates to a seismic signal processing method, belonging to the technical field of seismic interpretation processing. BACKGROUND

[0002] The post-stack seismic data (vertical) high resolution processing can improve the vertical identification accuracy of seismic data, highlight the seismic response characteristics of small scale geological bodies (reservoirs), and is an important link of seismic interpretation and reservoir prediction. The traditional post-stack seismic data high resolution processing method mainly takes deconvolution as the basic algorithm, and its principle is to restore the reflection coefficient from the given observation results (seismic signals), and then convolve with the reflection pulse (seismic wavelet) of the compression length to eliminate the earth filter effect of the seismic wave in the underground running process. In the theoretical model, this algorithm can eliminate the multiple wave influence and improve the vertical resolution of seismic data, but in actual application, there are limitations. The primary problem is that if the noise component of the seismic signal is more, improving the data resolution will lead to a large decrease in signal-to-noise ratio. Therefore, removing the noise of the seismic signal becomes the key problem of the subsequent research after high resolution processing, and the noise distribution rule of the seismic signal is that the higher the frequency, the greater the noise, so the time-frequency analysis of the seismic signal has become an important research field of high resolution processing in recent years.

[0003] The processing method commonly used in the prior art to remove the noise of the seismic signal is to perform empirical mode decomposition on the seismic signal, for example: the Chinese patent application file with the application publication number CN 108072894 A discloses a seismic signal processing method, which is a Hilbert spectrum whitening seismic data high resolution processing method, specifically: performing empirical mode decomposition (EMD) on the seismic signal to obtain a plurality of intrinsic mode functions (IMF) with different frequency ranges, then using spectral whitening to reasonably equalize the amplitude of each IMF component, and finally reconstructing the processed IMF component to improve the resolution of the seismic data.

[0004] However, the empirical mode decomposition method has mode aliasing phenomenon, which easily makes the generated intrinsic mode function lose physical meaning and geological meaning, thereby affecting the final high resolution processing result. Therefore, some people propose to decompose the seismic signal based on the variational mode decomposition method, for example: the author Zhang Lipeng, the name of "seismic signal high resolution processing method based on variational mode decomposition" of Chengdu University of Information Engineering Master's degree thesis, first obtains IMF by using variational mode decomposition algorithm (VMD), then performs Hilbert transform on each IMF to obtain instantaneous amplitude and frequency, then uses whitening filter to perform whitening processing on the Hilbert spectrum of each IMF, and finally restores the time domain signal after the amplitude and phase information after whitening processing.

[0005] However, the above method only obtains IMF by VMD, and the VMD method has the problem that the appropriate decomposition order is difficult to be adaptively selected, and cannot adaptively reflect the weak inherent characteristic difference of different seismic traces; in addition, whitening processing is performed on all IMFs, which is not conducive to highlighting weak hydrocarbon response, resulting in inaccurate high resolution processing result of seismic signal. SUMMARY

[0006] The purpose of the present application is to provide a seismic signal processing method to solve the problem of inaccurate existing processing method.

[0007] To achieve the above purpose, the present application provides a technical scheme of a seismic signal processing method, including the following steps:

[0008] 1) obtaining post-stack data, the post-stack data including each seismic trace signal;

[0009] 2) performing FFT transform on the post-stack data to obtain the transformed frequency spectrum;

[0010] 3) adaptively segmenting the frequency spectrum of each seismic trace signal to obtain segmented intervals and spectral boundary points;

[0011] 4) constructing an adaptive wavelet filter bank through the segmented intervals and spectral boundary points to filter the original seismic signal to obtain different empirical mode components;

[0012] 5) optimizing the different empirical mode components in step 4) through mutual information calculation to obtain an empirical mode component reflecting seismic wavelet information;

[0013] 6) energy equalizing the empirical mode component reflecting seismic wavelet information;

[0014] 7) reconstructing the energy equalized empirical mode component and the unoptimized empirical mode component to obtain a high resolution seismic data volume.

[0015] The technical scheme of the seismic signal processing method has the beneficial effects that: the adaptive wavelet filter bank is constructed after the spectrum of the seismic signal is segmented, the empirical mode components are obtained by filtering the original signal through the adaptive wavelet filter bank, then the empirical mode components are optimized to obtain the empirical mode components reflecting the seismic wavelet information, the amplitude spectrum of the signal is widened after the energy of the empirical mode components reflecting the seismic wavelet information is balanced, and then the seismic data volume is reconstructed.

[0016] Further, the time-frequency conversion method in the step 2) is FFT transformation.

[0017] Further, in order to avoid the influence of modal aliasing, the adaptive segmentation method in the step 3) is an alternating direction multiplier adaptive segmentation method.

[0018] Further, the adaptive wavelet filter bank is constructed through Littlewood-Paley and Meyers wavelet in the step 4).

[0019] Further, the calculation process of the mutual information in the step 5) is as follows:

[0020] I(c k ,c k+1 )=H(c k )+H(c k+1 )-H(c k ,c k+1 );

[0021] Wherein, I(c k ,c k+1 ) is the mutual information of the empirical mode component c k and the empirical mode component c k+1 ; H(·) is the calculation of Shannon entropy.

[0022] Further, the energy balancing process in the step 6) is as follows:

[0023]

[0024] Wherein, is the empirical mode component after energy balancing; e k (ω) is the envelope of the empirical mode component reflecting the seismic wavelet information; ε is a noise factor; max(·) is the maximum value operation. BRIEF DESCRIPTION OF DRAWINGS

[0025] Figure 1is a flow chart of the seismic signal processing method of the present application;

[0026] Figure 2 is a certain well profile of a certain area in a certain basin of the present application;

[0027] Figure 3 is a frequency spectrum analysis chart of 2500-2900 ms of a certain well profile in a certain area of the present application;

[0028] Figure 4 is a certain well profile of the seismic data after resolution enhancement of the present application;

[0029] Figure 5 is a frequency spectrum analysis chart of 2500-2900 ms of a certain well profile after resolution enhancement of the seismic data of the present application. DETAILED DESCRIPTION

[0030] Embodiment of the seismic signal processing method:

[0031] The main idea of the seismic signal processing method is that, based on the problem that the existing high resolution processing result of the seismic signal is not accurate, the present application constructs an adaptive wavelet filter set in combination with the spectrum boundary in the segmented interval after adaptive segmentation of the seismic signal, filters the original signal to obtain an empirical mode component; then each empirical mode component is optimized, and the empirical mode component reflecting the seismic wavelet information is selected for energy equalization, and then the equalized empirical mode component and the unequilibrated empirical mode component are reconstructed together to generate a data body after improvement of the signal-to-noise ratio, realizing high resolution processing of the seismic signal.

[0032] Specifically, taking the seismic signal of a certain reservoir in a certain area of a certain basin as an example, the seismic signal processing method of the present application is described, as shown in Figure 1 , comprising the following steps:

[0033] 1) Obtain the stacked data, and determine the analysis time window of the stacked data, wherein each seismic trace signal is included in the stacked data.

[0034] Obtain the three-dimensional stacked pure wave data (i.e. stacked data) of a certain area in a certain basin, accurately mark the horizon of the target layer by using the field geological information, drilling, well logging and fine well-seismic calibration, and determine the analysis time window of the stacked pure wave data. A certain well profile in the data is shown in Figure 2 , which is a two-dimensional stacked migration profile, the analysis time window of the profile is 2500-2900 ms, and the reservoir type is a carbonate reef bank reservoir.

[0035] 2) In the analysis time window range determined in step 1), perform FFT transformation on the stacked data to obtain the transformed frequency spectrum.

[0036] The process of FFT transformation is:

[0037]

[0038] where, is the Fourier transform spectrum of the time domain stacked data; ω is the frequency; s(t) is the time domain stacked data; t is the time, and j is the imaginary unit.

[0039] The spectrum obtained after FFT transform of the target layer is shown in FIG. 1, and the frequency band range is 6-66 Hz, and the main frequency is 25 Hz. Figure 3

[0040] 3) The spectrum of each seismic trace signal is adaptively segmented by using the alternating direction multiplier method to obtain the segmented interval (here, the segmented interval represents the frequency band) and the spectral boundary point.

[0041] The bandwidth objective function for spectrum segmentation by using the alternating direction multiplier method is:

[0042]

[0043]

[0044] where, u k is the kth narrowband function after segmentation; ω κ is the center frequency of the kth narrowband function; K represents the number of segmented intervals (generally selected as 3 or 4); δ(t) is the impulse function; j is the complex unit; represents the Hilbert transform of the original signal; is the partial derivative with respect to time t; ‖‖2 is the second order norm.

[0045] To solve the above formula, a quadratic penalty factor and a Lagrange multiplier are introduced, and the above constraint problem can be changed into an unconstrained problem:

[0046]

[0047] where, λ is the Lagrange multiplier; α is the balance parameter of data fidelity; L(u k ,ω k ,λ) is the objective function; δ t (.) is the impulse function.

[0048] Further, the different narrowband functions in the spectral domain are:

[0049]

[0050] where, is the Fourier transform of the kth narrowband function after the n+1th iteration; is the Fourier transform of the Lagrange multiplier after the nth iteration;​ Here, K is the center frequency of the k-th narrowband function after the n-th iteration; n represents the number of iterations; K is the number of narrowband functions, which is the same as the number of intervals divided, and is generally set to 3 or 4. Based on the actual data, K = 3 in this case.

[0051] The center frequency ω of the k-th narrowband function κ The iterative update is as follows:

[0052]

[0053] Based on the center frequency ω of the kth narrowband function κ The transition zone β between the adjacent intervals of the k-th narrowband function κ The spectral boundary point Ω of the k-th narrowband function can be calculated using the following formula. κ :

[0054]

[0055] Among them, Ω O =0,Ω κ =π.

[0056] 4) An adaptive wavelet filter bank is constructed by dividing the intervals and spectral boundaries to filter the original seismic signal and obtain different empirical mode components.

[0057] In this step, the Littlewood-Paley and Meyers wavelets are used to construct an adaptive wavelet filter bank. The construction process is as follows:

[0058]

[0059]

[0060]

[0061] in, It is an empirical scaling function; Ω is the empirical wavelet function; ω is the frequency; k It is the k-th boundary frequency; γ is a parameter to ensure that there is no overlap between two consecutive intervals; β (χ) Defined as the following function:

[0062]

[0063] The original seismic signal was filtered using a constructed adaptive wavelet filter, and the approximation coefficients W were obtained. s (0,t) is:

[0064]

[0065] Detail coefficient W s (k,t) is obtained from the inner product of the signal and the wavelet function:

[0066]

[0067] Thus, the empirical mode component c k is:

[0068] c0(t) = W s (0,t)*φ1(t);

[0069] c k (t) = W s (k,t)*ψ k (t);

[0070] where c0(t) is the mode component with the spectral boundary of 0; c k (t) is the kth empirical mode component; k = 1,..., n, n is the number of empirical mode components, each empirical mode component is a narrowband function; Φ1 is the time domain scaling function; F -1 is the inverse Fourier transform; is the complex conjugate of the corresponding function; t is time; is the complex conjugate of the frequency domain empirical scaling function; f is frequency; is the empirical wavelet function; is the complex conjugate of the frequency domain wavelet function; φ1(t) is the empirical scaling function corresponding to the initial boundary point; ψ k (t) is the wavelet function corresponding to the kth boundary point.

[0071] 5) The different empirical mode components in step 4) are optimized using the mutual information principle to obtain the empirical mode component reflecting the information of the seismic wavelet.

[0072] The optimization of the empirical mode component is performed by calculating the mutual information of each empirical mode component, and the calculation process of the mutual information is as follows:

[0073] I(c k ,c k+1 ) = H(c k ) + H(c k+1 ) - H(c k ,c k+1 );

[0074] where I(c k ,c k+1 ) is the mutual information of the empirical mode component c k and the empirical mode component c k+1 ; H(·) is the calculation of Shannon entropy.

[0075] The final preferred set of empirical mode components reflecting seismic wavelet information can be expressed as: (c1, c2, …, cN) G , where G is the number of empirical mode components in the set,

[0076] 6) Energy equalization of the empirical mode components reflecting seismic wavelet information in step 5) to broaden the amplitude spectrum of the empirical mode components reflecting seismic wavelet information.

[0077] The process of energy equalization is as follows:

[0078]

[0079] wherein, is the energy equalized empirical mode component; e k (ω) is the envelope (envelope is calculated according to the amplitude spectrum) of the empirical mode component c k reflecting seismic wavelet information; ε is a noise factor, the value of which will affect the resolution and signal-to-noise ratio of the final processing result, and thus needs to be determined according to the actual seismic data; max(·) is the maximum value operation.

[0080] 7) Reverse transformation of the non-preferred empirical mode components in step 5) and the energy equalized empirical mode components in step 6) to the time domain to complete the reconstruction of the seismic data s d (t) together, and the reconstructed seismic data improves the resolution while maintaining the signal-to-noise ratio of the original signal.

[0081] The transformation process of the energy equalized empirical mode components reverse transformed to the time domain is as follows:

[0082]

[0083] wherein, is the reverse transformed energy equalized empirical mode component; F -1 is the reverse transformation function.

[0084] Here, the non-preferred empirical mode components are reconstructed together to ensure the integrity of the signal and not to discard the original component signal.

[0085] 8) The above completes the data reconstruction at a certain time point in a certain seismic trace, and then steps 3) to 7) are repeated point by point (time point) and trace by trace (seismic trace) to obtain a high-resolution seismic data volume as shown in Figure 4 , and further complete the identification of the target layer.

[0086] Comparison Figure 2 and Figure 4It can be seen that the geological regularity reflected by the seismic data volume after improving the resolution is consistent, and the response characteristics of the reservoir inside the reef are more obvious. It is beneficial to the division and identification of the reef development stages.

[0087] The seismic data volume obtained in step 8) is subjected to FFT transformation to obtain a frequency spectrum analysis diagram as shown in Figure 5 , the frequency band range is 6-74 Hz, the main frequency is 30 Hz, and the frequency spectrum analysis diagram is obtained by Figure 3 and Figure 5 It can be seen from the comparison that the data volume after the resolution improvement processing maintains the low-frequency information and widens the high-frequency information, so that the resolution of the seismic data is effectively improved, and the fidelity of the seismic data is also ensured.

Claims

1. A method of seismic signal processing, characterized by, The method comprises the following steps: 1) obtaining post-stack data, determining an analysis time window of the post-stack data, wherein the post-stack data comprises each seismic trace signal; 2) performing FFT transform on the post-stack data in the analysis time window to obtain transformed frequency spectrum; 3) performing adaptive segmentation on the frequency spectrum of each seismic trace signal to obtain segmented intervals and spectral boundary points; 4) constructing an adaptive wavelet filter set through the segmented intervals and spectral boundary points to filter the original seismic signal to obtain different empirical mode components; 5) performing optimization on the different empirical mode components in step 4) through calculation of mutual information to obtain an empirical mode component reflecting seismic wavelet information; 6) performing energy equalization on the empirical mode component reflecting seismic wavelet information to expand the amplitude spectrum thereof, wherein the energy equalization process is as follows: wherein, is the energy-balanced empirical mode component; e k (ω) is the envelope of the empirical mode component reflecting the seismic wavelet information; ε is the noise factor; max(·) is the maximum value operation; 7) reconstructing the energy-equalized empirical mode component and the non-optimized empirical mode component to obtain a high-resolution seismic data volume.

2. The seismic signal processing method of claim 1, wherein, The time-frequency conversion method in step 2) is FFT transform.

3. The seismic signal processing method of claim 1, wherein, The adaptive segmentation method in step 3) is an alternating direction multiplier adaptive segmentation method.

4. The seismic signal processing method of claim 1, wherein, The adaptive wavelet filter set in step 4) is constructed through Littlewood-Paley and Meyers wavelets.

5. The seismic signal processing method of claim 1, wherein, The calculation process of mutual information in step 5) is as follows: I(c k ,c k+1 ) = H(c k ) + H(c k+1 ) - H(c k ,c k+1 ); where I(c k ,c k+1 ) is the mutual information of empirical mode component c k and empirical mode component c k+1 ; H(·) is the calculation of shannon entropy.

Images