Time domain compression seismic spectrum decomposition method based on matching pursuit

Abstract

The invention relates to the technical field of oil and gas geophysical exploration, in particular to a time domain compression seismic spectrum decomposition method based on matching pursuit, which comprises the following steps of: performing sparse decomposition on seismic signals, namely searching to obtain an optimal matching wavelet, and further obtaining a reflection wavelet; fourier transform and time compression transform are carried out on the reflection wavelets to obtain time domain compression frequency spectrums of the reflection wavelets, and superposition summation is further carried out on the time domain compression frequency spectrums of the reflection wavelets, so that the time resolution of the frequency spectrums during seismic recording is improved, false frequency spectrum abnormity is avoided, and the accuracy of an output result is guaranteed. The time domain compression method is used for improving the time-frequency resolution of the data, fine features in the signals and the characteristic of changing along with time can be highlighted, the precision and reliability of signal processing are improved, and when the method is applied to other fields except for seismic exploration, signal feature recognition can be enhanced, and the signal detection capacity can be improved.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas geophysical exploration, and more specifically, to a time-domain compressed seismic spectrum decomposition method based on matching pursuit. Background Art

[0002] Seismic spectrum decomposition is an important seismic data processing and interpretation technology, which has been widely used in many fields such as description of seismic phases, formation absorption compensation, oil and gas detection, noise suppression, thin layer identification, etc. Applying spectrum decomposition technology to pre-stack seismic data and combining it with AVO (reflection amplitude varies with offset) technology to perform frequency-variable AVO dispersion parameter inversion can realize hydrocarbon detection.

[0003] Traditional spectral decomposition methods, such as Fourier transform or wavelet transform, may have insufficient resolution when processing complex signals; Wigner-Ville distribution (WVD) is a commonly used time-frequency analysis method with high time-frequency resolution, but it has the problem of cross-term interference when analyzing multi-component signals. As an adaptive signal decomposition method, the matching pursuit algorithm can effectively sparsely represent the signal and improve the time-frequency resolution of the signal. The matching pursuit algorithm needs to be combined with other time-frequency analysis methods to calculate the time-frequency spectrum of the signal. The WVD spectrum decomposition method based on matching pursuit aims to solve the cross-term interference problem of WVD in multi-component signal analysis and improve the time-frequency resolution of the signal. However, the WVD based on matching pursuit is an energy spectrum, which is prone to produce false spectrum anomalies and affect the detection effect of hydrocarbons. The instantaneous spectrum based on matching pursuit uses the decomposed wavelet amplitude envelope to calculate the wavelet time-frequency spectrum, and the obtained seismic record time-frequency spectrum has a low time resolution. Summary of the invention

[0004] In order to overcome the defect of low time resolution of seismic record time spectrum obtained by calculating wavelet time spectrum by decomposing wavelet amplitude envelope based on instantaneous spectrum of matching pursuit in the above-mentioned prior art, the present invention provides a time domain compressed seismic spectrum decomposition method based on matching pursuit to improve the time resolution of seismic record time spectrum.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is: a time domain compressed seismic spectrum decomposition method based on matching pursuit, comprising the following steps:

[0006] Step 1: Input the seismic signal and select the appropriate wavelet atom library;

[0007] Step 2: Based on the matching pursuit algorithm, the seismic signal is sparsely decomposed, the best matching wavelet is searched, the corresponding wavelet coefficients are calculated, and a series of reflection wavelets are obtained;

[0008] Step 3: Given a time compression factor, perform Fourier transform and time compression transform on each matching wavelet obtained by decomposition, and calculate the time-domain compressed time-frequency spectrum of the reflected wavelet;

[0009] Step 4: Superimpose and sum the time-domain compressed time-frequency spectra of the reflected wavelets to obtain the time-frequency decomposition result of the entire original seismic signal.

[0010] Preferably, in step 2, the seismic signal is used as the residual signal, the number of iterations is set, the best matching wavelet with the largest correlation coefficient with the current residual signal is searched within the number of iterations, the wavelet coefficient corresponding to the best matching wavelet is calculated, and the residual signal is updated using the wavelet coefficient. When the number of iterations or the residual signal energy is lower than a given threshold, the iteration is stopped to obtain a series of reflection wavelets.

[0011] Preferably, in step 2, for the kth iteration, the following formula is used to search for the best matching wavelet with the largest correlation coefficient with the current residual signal:

[0012]

[0013] Among them, s(t) represents the seismic signal, r 0 (t) represents the residual signal, w(t) represents the wavelet atom, r k-1 (t),w(t)> represents the inner product of the residual signal and the wavelet atom, is the normalization operator of the wavelet atom, and k represents the kth iteration.

[0014] Preferably, in step 2, the formula for calculating the wavelet coefficient corresponding to the best matching wavelet is:

[0015] C k represents the wavelet coefficients at the kth iteration;

[0016] The formula for updating the residual signal is: k (t) = r k-1 (t)-C k w k (t).

[0017] Preferably, in step 2, a series of reflected wavelets are obtained by the following formula:

[0018] Among them, C k w k (t) represents the reflected wavelet.

[0019] Preferably, in step 3, the formula for performing Fourier transform on each matching wavelet is: Where f represents frequency; each matching wavelet is compressed in the time direction, and the formula is: w k (t)→wk (at); where a is the time compression factor.

[0020] Preferably, in step 3, the time-domain compressed time-frequency spectrum calculation formula of each reflected wavelet is:

[0021] S k (t,f)=C k ·W k (f)env[w k (at)];

[0022] Among them, env[w k (at)] represents the amplitude envelope of the time-domain compressed wavelet. By compressing the wavelet in the time direction, the purpose of reducing the wavelet duration and improving the wavelet resolution is achieved.

[0023] Preferably, in step 4, the time-domain compression frequency spectrum superposition summation formula of the reflected wavelet is:

[0024]

[0025] Among them, S(t,f) is the time-frequency spectrum of the entire original seismic signal. Compared with the conventional time-frequency spectrum, this time-frequency spectrum has a higher time resolution and can more clearly describe the change of signal spectrum energy over time.

[0026] Preferably, in step 1, the input seismic signal is a post-stack seismic signal that has undergone a series of data processing such as extraction, dynamic correction, migration and stacking.

[0027] Preferably, in step 1, the wavelet atom library contains a series of Ricker wavelets w(τ) with different main frequencies, and the calculation formula is:

[0028] w(τ)={1-2[πf d τ] 2}exp{-[πf d τ] 2};

[0029] Among them, f d is the main frequency of the Ricker wavelet, the frequency is 5Hz~50Hz, and τ is the sampling time. The Ricker wavelet can accurately describe the propagation and reflection characteristics of seismic waves.

[0030] Compared with the prior art, the invention has the following beneficial effects: by sparsely decomposing the seismic signal, that is, searching for the best matching wavelet, and then further obtaining the reflected wavelet, and then obtaining the time domain compressed spectrum of the reflected wavelet by Fourier transform and time compression transform of the reflected wavelet, and then further superimposing and summing the time domain compressed spectrum of the reflected wavelet, the time resolution of the seismic record time spectrum is improved, false spectrum anomalies are avoided, and the accuracy of the output result is guaranteed. Using the time domain compression method to improve the time and frequency resolution of data can highlight the subtle features in the signal and its characteristics of changing over time, improve the accuracy and reliability of signal processing, and apply it to other fields besides seismic exploration, which can enhance signal feature recognition and improve signal detection capabilities. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is a flow chart of a time domain compression seismic spectrum decomposition method based on matching pursuit of the present invention;

[0032] Figure 2 It is a schematic diagram of seismic signals of a time-domain compression seismic spectrum decomposition method based on matching pursuit of the present invention;

[0033] Figure 3 It is a schematic diagram of time-frequency spectrum decomposition results of a time-domain compression seismic spectrum decomposition method based on matching pursuit of the present invention;

[0034] Figure 4 It is a comparison diagram of time-frequency spectrum decomposition results of a time-domain compression seismic spectrum decomposition method based on matching pursuit of the present invention. DETAILED DESCRIPTION

[0035] The drawings are only for illustrative purposes and cannot be construed as limiting the present invention. To better illustrate the present embodiment, some parts of the drawings may be omitted, enlarged, or reduced, and do not represent the size of the actual product. For those skilled in the art, it is understandable that some well-known structures and their descriptions may be omitted in the drawings. The positional relationships described in the drawings are only for illustrative purposes and cannot be construed as limiting the present invention.

[0036] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right", "long", "short" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limitations on this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0037] The technical solution of the present invention is further described in detail below through specific embodiments and in conjunction with the accompanying drawings:

[0038] Example 1

[0039] like Figure 1 As shown, a time domain compressed seismic spectrum decomposition method based on matching pursuit includes the following steps:

[0040] Step 1: Input the seismic signal and select the appropriate wavelet atom library;

[0041] Step 2: Based on the matching pursuit algorithm, the seismic signal is sparsely decomposed, the best matching wavelet is searched, the corresponding wavelet coefficients are calculated, and a series of reflection wavelets are obtained;

[0042] Step 3: Given a time compression factor, perform Fourier transform and time compression transform on each matching wavelet obtained by decomposition, and calculate the time-domain compressed time-frequency spectrum of the reflected wavelet;

[0043] Step 4: Superimpose and sum the time-domain compressed time-frequency spectra of the reflected wavelets to obtain the time-frequency decomposition result of the entire original seismic signal.

[0044] The beneficial effects of this embodiment are as follows: by sparsely decomposing the seismic signal, that is, searching for the best matching wavelet, and then further obtaining the reflected wavelet, and then obtaining the time domain compressed spectrum of the reflected wavelet by Fourier transform and time compression transform of the reflected wavelet, and then further superimposing and summing the time domain compressed spectrum of the reflected wavelet, the time resolution of the seismic record time spectrum is improved, and false spectrum anomalies are avoided, and the accuracy of the output result is guaranteed. Using the time domain compression method to improve the time and frequency resolution of data can highlight the subtle features in the signal and its characteristics of changing over time, improve the accuracy and reliability of signal processing, and apply it to other fields besides seismic exploration, which can enhance signal feature recognition and improve signal detection capabilities.

[0045] Example 2

[0046] Based on Example 1, the difference from Example 1 is that:

[0047] Furthermore, in step 2, the seismic signal is used as the residual signal, the number of iterations is set, the best matching wavelet with the largest correlation coefficient with the current residual signal is searched within the number of iterations, the wavelet coefficient corresponding to the best matching wavelet is calculated, and the residual signal is updated using the wavelet coefficient. When the number of iterations or the residual signal energy is lower than a given threshold, the iteration is stopped to obtain a series of reflection wavelets.

[0048] For the kth iteration, the following formula is used to search for the best matching wavelet with the largest correlation coefficient with the current residual signal:

[0049]

[0050] Among them, s(t) represents the seismic signal, r 0 (t) represents the residual signal, w(t) represents the wavelet atom, r k-1 (t),w(t)> represents the inner product of the residual signal and the wavelet atom, is the normalization operator of the wavelet atom, and k represents the kth iteration.

[0051] The formula for calculating the wavelet coefficients corresponding to the best matching wavelet is:

[0052] C k represents the wavelet coefficients at the kth iteration;

[0053] The formula for updating the residual signal is: k (t) = r k-1 (t)-C k w k (t).

[0054] A series of reflected wavelets are obtained by the following formula:

[0055] Among them, C k w k (t) represents the reflected wavelet.

[0056] The remaining features and working principles of this embodiment are consistent with those of Embodiment 1.

[0057] Example 3

[0058] Based on Example 1 or Example 2, Example 1 or Example 2 is further limited, and the difference is that:

[0059] Furthermore, in step 3, the formula for Fourier transform of each matching wavelet is: Where f represents frequency; each matching wavelet is compressed in the time direction, and the formula is: w k (t)→w k (at); where a is the time compression factor.

[0060] Among them, the time-domain compressed time-frequency spectrum calculation formula of each reflected wavelet is:

[0061] S k (t,f)=C k ·W k (f)env[w k (at)];

[0062] Among them, env[w k(at)] represents the amplitude envelope of the time-domain compressed wavelet. By compressing the wavelet in the time direction, the purpose of reducing the wavelet duration and improving the wavelet resolution is achieved.

[0063] Furthermore, in step 4, the time-domain compression frequency spectrum superposition summation formula of the reflected wavelet is:

[0064]

[0065] Among them, S(t,f) is the time-frequency spectrum of the entire original seismic signal. Compared with the conventional time-frequency spectrum, this time-frequency spectrum has a higher time resolution and can more clearly describe the change of signal spectrum energy over time.

[0066] Furthermore, in step 1, if Figure 2 As shown in FIG. 1 , the input seismic signal is a post-stack seismic signal after a series of data processing such as extraction, dynamic correction, migration and stacking.

[0067] Furthermore, in step 1, the wavelet atom library contains a series of Ricker wavelets w(τ) with different main frequencies, and the calculation formula is:

[0068] w(τ)={1-2[πf d τ] 2}exp{-[πf d τ] 2};

[0069] Among them, f d is the main frequency of the Ricker wavelet, the frequency is 5Hz~50Hz, and τ is the sampling time. The Ricker wavelet can accurately describe the propagation and reflection characteristics of seismic waves.

[0070] Example 4

[0071] In this embodiment, given a time compression factor a=2, spectral decomposition is performed by the following steps:

[0072] Step 1: If Figure 2 As shown in the figure, the seismic signal after a series of data processing such as extraction, dynamic correction, migration and stacking is input, and a suitable wavelet atom library is selected. The wavelet atom library contains a series of Ricker wavelets w(τ) with different main frequencies. The calculation formula is:

[0073] w(τ)={1-2[πf d τ] 2}exp{-[πf d τ] 2};

[0074] Among them, f d is the main frequency of the Ricker wavelet, the frequency is 5Hz~50Hz, τ is the sampling time;

[0075] Step 2: Based on the matching pursuit algorithm, the seismic signal is sparsely decomposed, the best matching wavelet is searched, the corresponding wavelet coefficients are calculated, and a series of reflection wavelets are obtained;

[0076] The input seismic signal s(t) is regarded as the residual signal r 0 (t):

[0077] r 0 (t) = s(t)

[0078] Here, t represents time.

[0079] For the kth iteration, search for the best matching wavelet with the largest correlation coefficient with the current residual signal,

[0080]

[0081] Among them, w(t) represents the wavelet atom, <r k-1 (t),w(t)> represents the inner product of the residual signal and the wavelet atom, is the normalization operator of the wavelet atom.

[0082] Calculate the wavelet coefficients corresponding to the best matching wavelet:

[0083]

[0084] Update the residual signal,

[0085] r k (t) = r k-1 (t)-C k w k (t)

[0086] When the number of iterations k or the residual signal energy is lower than a given threshold, the iteration stops. The seismic signal is represented by a series of reflected wavelets:

[0087]

[0088] Step 3: Perform Fourier transform and time compression transform on each matching wavelet, and calculate the time-domain compressed time-frequency spectrum of the reflected wavelet;

[0089] For each matching wavelet w k (t) Perform Fourier transform:

[0090]

[0091] Wherein, f represents frequency.

[0092] Given a time compression factor, each matching wavelet is compressed in the time direction:

[0093] w k (t)→w k (at)

[0094] The time-domain compressed time-frequency spectrum of each reflected wavelet is:

[0095] S k (t,f)=C k ·W k (f)env[w k (at)]

[0096] Among them, env[w k (at)] represents the amplitude envelope of the time-domain compressed wavelet.

[0097] Step 4: Superimpose and sum the time-frequency spectrum of the reflected wavelet in the time domain to obtain the time-frequency spectrum decomposition result of the entire original seismic signal:

[0098]

[0099] Among them, S(t,f) is the time-frequency spectrum of the entire original seismic signal.

[0100] Get as Figure 3 The time-frequency spectrum decomposition result is shown. Figure 4 This is the time-frequency spectrum decomposition result obtained by the ordinary method. In comparison, the time-frequency spectrum decomposition result calculated by the present application has a higher resolution.

[0101] The rest of the working principle and working process of this embodiment are consistent with those of Embodiment 1 or Embodiment 2.

[0102] In the specific contents of the above-mentioned specific implementation methods, the various technical features can be combined in any non-contradictory manner. In order to make the description concise, not all possible combinations of the above-mentioned technical features are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0103] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A time-domain compressed seismic spectrum decomposition method based on matching pursuit, characterized in that: The following steps are involved: Step 1: Input the seismic signal and select the appropriate wavelet atom library; Step 2: Based on the matching pursuit algorithm, the seismic signal is sparsely decomposed, the best matching wavelet is searched, the corresponding wavelet coefficients are calculated, and a series of reflection wavelets are obtained; Step 3: Given a time compression factor, perform Fourier transform and time compression transform on each matching wavelet obtained by decomposition, and calculate the time-domain compressed time-frequency spectrum of the reflected wavelet; Step 4: Superimpose and sum the time-domain compressed time-frequency spectra of the reflected wavelets to obtain the time-frequency decomposition result of the entire original seismic signal.

2. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 1, characterized in that: In step 2, the seismic signal is used as the residual signal, the number of iterations is set, and the best matching wavelet with the largest correlation coefficient with the current residual signal is searched within the number of iterations. The wavelet coefficient corresponding to the best matching wavelet is calculated, and the residual signal is updated using the wavelet coefficient. When the number of iterations or the residual signal energy is lower than a given threshold, the iteration is stopped to obtain a series of reflection wavelets.

3. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 2, characterized in that: In step 2, for the kth iteration, the following formula is used to search for the best matching wavelet with the largest correlation coefficient with the current residual signal: r 0 (t)=s(t), Among them, s(t) represents the seismic signal, r 0 (t) represents the residual signal, w(t) represents the wavelet atom, r k-1 (t),w(t)> represents the inner product of the residual signal and the wavelet atom, is the normalization operator of the wavelet atom, and k represents the kth iteration.

4. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 3 is characterized by: In step 2, the formula for calculating the wavelet coefficient corresponding to the best matching wavelet is: C k represents the wavelet coefficients at the kth iteration; The formula for updating the residual signal is: k (t) = r k-1 (t)-C k w k (t).

5. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 4, characterized in that: In step 2, a series of reflected wavelets are obtained by the following formula: Among them, C k w k (t) represents the reflected wavelet.

6. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 1, characterized in that: In step 3, the formula for Fourier transform of each matching wavelet is: Where f represents frequency; each matching wavelet is compressed in the time direction, and the formula is: w k (t)→w k (at); where a is the time compression factor.

7. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 6, characterized in that: In step 3, the time-domain compressed time-frequency spectrum calculation formula of each reflected wavelet is: S k (t,f)=C k ·W k (f)env[w k (at)]; Among them, env[w k (at)] represents the amplitude envelope of the time-domain compressed wavelet.

8. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 1, characterized in that: In step 4, the time-domain compression and frequency spectrum superposition summation formula of the reflected wavelet is: Among them, S(t,f) is the time-frequency spectrum of the entire original seismic signal.

9. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 1, characterized in that: In step 1, the input seismic signal is a post-stack seismic signal that has undergone a series of data processing such as extraction, dynamic correction, migration and stacking.

10. The time-domain compressed seismic spectrum decomposition method based on matching pursuit according to claim 1, characterized in that: In step 1, the wavelet atom library contains a series of Ricker wavelets w(τ) with different main frequencies, and the calculation formula is: w(τ)={1-2[πf d [t] 2 }exp{-[πf d [t] 2 }; Among them, f d is the main frequency of the Ricker wavelet, the frequency is 5Hz~50Hz, and τ is the sampling time.

Images