A time-synchronous compression extraction and transformation method for time-frequency spectrum

Abstract

The invention relates to a time-synchronous compression extraction transformation method for a time-frequency spectrum, and belongs to the technical field of oil and gas field exploration. The method first performs conventional time-frequency transformation on seismic data to calculate a time-frequency spectrum and a group delay operator, then constructs a time-synchronous compression extraction operator through group delay and a Gaussian window, and finally performs synchronous compression extraction transformation on the time-frequency spectrum according to the synchronous compression extraction operator, thereby obtaining an amplitude spectrum with a certain width in the time direction after the compression extraction transformation, and the amplitude spectrum of each frequency has the same time resolution. Compared with the amplitude spectrum before the compression extraction transformation, the amplitude spectrum has improved resolution in the time direction, thereby improving the accuracy of time-frequency analysis. Compared with the amplitude spectrum obtained by the conventional time-synchronous compression method, the amplitude spectrum of each frequency has the same time resolution, and the reliability of comparative analysis of amplitude spectrums of multiple frequencies is improved.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas field exploration, and in particular to a time synchronization compression extraction and transformation method of a time spectrum. Background Art

[0002] Time-frequency analysis is a method that uses the joint density function of time and frequency to describe the energy density or intensity of the signal in the time domain and frequency domain, and analyzes and processes it. Since time-frequency analysis enables us to analyze the energy distribution of seismic signals in the time-frequency domain, time-frequency analysis technology is widely used in the field of geophysical exploration, such as oil and gas identification, thin layer identification, river sand body detection, etc. After years of theoretical expansion, time-frequency analysis has become the mainstream seismic data processing technology. Conventional time-frequency analysis methods can be roughly divided into two categories: (1) linear time-frequency analysis methods: such as short-time Fourier transform (STFT), wavelet transform (WT) and S transform, etc., this type of method obtains time-frequency characteristics by calculating the convolution of the signal and the basis function; (2) bilinear time-frequency analysis methods: such as Wigner-Ville distribution and its variants, etc., this type of method obtains higher time-frequency resolution by calculating the Fourier transform of its local signal correlation.

[0003] However, conventional time-frequency methods have some defects. For example, the linear time-frequency analysis method is limited by the Heisenberg uncertainty principle, that is, it cannot simultaneously obtain high time and frequency resolution (Zhang Jinqiang, 2023), and due to the poor time-frequency focusing ability, the spectrum values ​​will be distributed on a band in the time-frequency spectrum, and its geometry depends on the signal and the window. Therefore, the wrong energy distribution will be observed at the position where energy should not exist on the time-frequency spectrum (Zhonglai Huang, 2016). The bilinear time-frequency analysis method will introduce unexpected cross terms for multi-component signals, and the time-frequency spectrum will be affected by the overlap or cross of energy density, making it difficult to complete the effective identification of the signal (Zhang Hongmao, 2023).

[0004] In order to overcome the defect of insufficient resolution of conventional time-frequency analysis methods, some post-processing methods based on traditional time-frequency analysis methods have been proposed in recent years, such as the reallocation method (RM); synchronous compression methods in the frequency direction, such as the synchronized compression transform (SST), synchronized compression extraction transform (SET), etc.; synchronous compression methods in the time direction, such as the time synchronized compression transform (TSST), time synchronized compression extraction transform (TET), etc.

[0005] The RM method achieves better time-frequency positioning of the signal by redistributing the time-frequency coefficients to the instantaneous frequency trajectory, but it does not allow the reconstruction of the signal (Gang Yu, 2021). Synchronous compression methods in the frequency direction, such as synchronous compression and synchronous compression extraction, only compress the time-frequency coefficients into the instantaneous frequency in the frequency direction, so they have the ability to reconstruct the signal and can provide higher time-frequency resolution than the RM method (Zhen Li, 2020). However, this type of method only performs synchronous compression in the frequency direction and cannot provide good time-frequency resolution in the time direction (F. Auger, 2013). The synchronous compression method in the time direction introduces a group delay operator to synchronize the time-frequency spectrum, which can achieve a relatively highly concentrated time-frequency resolution for pulse-like signals (Dong He, 2019). Synchronous compression methods in the time direction include time synchronous compression transform (TSST) and time synchronous compression extraction transform (TET), among which TET has stronger noise resistance than TSST (Gang Yu, 2021). However, when processing seismic signals, the synchronous compression method in the time direction often compresses the energy into a ridge line in the time-frequency spectrum, which is not conducive to practical analysis and processing, such as comparative analysis of amplitude spectra of multiple frequencies. The present invention proposes a time-frequency synchronous compression extraction transformation method for time-frequency spectrum, which aims to solve the problems existing in conventional technologies (such as TET, etc.) by constructing a time-frequency synchronous compression extraction operator, compressing the amplitude spectrum energy in the time direction within a certain width range (advantage 1: improving resolution), and maintaining the amplitude spectrum of a single frequency with the same time resolution (advantage 2: comparative analysis of amplitude spectra of multiple frequencies can be performed), thereby improving the accuracy and reliability of reservoir prediction based on time-frequency spectrum analysis. Summary of the invention

[0006] In order to solve the problem that the conventional time-frequency analysis method has low resolution and the time synchronization compression method in the high-resolution time-frequency analysis method compresses the energy into a ridge line on the time-frequency spectrum, which is not conducive to the analysis of actual data (such as the inability to perform comparative analysis of the amplitude spectra of multiple frequencies), the present invention proposes a time synchronization compression extraction transformation method for the time spectrum. The method first performs conventional time-frequency transformation (such as short-time Fourier transform, wavelet transform, S transform, etc.) on the data to calculate the time spectrum, then calculates the group delay operator of the data, and constructs a synchronization compression extraction operator based on the group delay operator, and finally performs synchronization compression extraction transformation on the time spectrum based on the synchronization compression extraction operator, thereby obtaining an amplitude spectrum with a certain width in the time direction after the compression extraction transformation, and the amplitude spectrum of each frequency has the same time resolution. Compared with the amplitude spectrum before compression extraction transformation (amplitude spectrum of conventional time-frequency method), the amplitude spectrum has improved resolution in the time direction, thereby improving the accuracy of time-frequency analysis; compared with the amplitude spectrum obtained by conventional time synchronization compression method (such as the amplitude spectrum obtained by TET method), its amplitude spectrum energy is no longer a ridge line on the time-frequency spectrum, but is distributed within a certain width range, and the amplitude spectrum of each frequency has the same time resolution, which improves the reliability of comparative analysis of amplitude spectra of multiple frequencies.

[0007] The specific technical solutions are as follows:

[0008] S1: Perform conventional time-frequency transformation on the data to calculate the time-frequency spectrum. Taking the wavelet transform method in conventional time-frequency transformation as an example, the calculation equation is as follows:

[0009]

[0010] Where x(t) is the time domain representation of the data, W(a,b) is the wavelet transform result, where a is the scale, b is the time, ψ(t) is the wavelet basis function, and () * Represents complex conjugate. In practical applications, complex Morlet wavelet is usually selected as the wavelet basis function. Its expression is as follows:

[0011]

[0012] where f c is the center frequency of the complex Morlet wavelet, f b is the bandwidth frequency of the complex Morlet wavelet;

[0013] S2: Take tψ(t) as the basis function and substitute it into the wavelet transform formula to get W tg (a,b), its expression is as follows,

[0014]

[0015] S3: Calculate the group delay operator, which is expressed as follows,

[0016]

[0017] In the formula, This is the obtained group delay operator, which is in complex form, but only its real part is used in actual use, so the real part of the result needs to be processed;

[0018] S4: The obtained group delay operator in the time scale domain And wavelet transform W(a,b) converted to the time frequency domain to obtain W(ω,t), where the conversion expression between scale and frequency is as follows,

[0019]

[0020] S5: Generate a one-dimensional array of Gaussian window, the expression is as follows,

[0021]

[0022] Assuming that the total length of the Gaussian window is 2*win+1, the value range of x is 1 to 2*win+1, where c is the position of the center point of the Gaussian window, which is win, σ is the Gaussian standard deviation, and the maximum value of the Gaussian window is 1;

[0023] S6: Get the synchronous compression extraction operator TEO(ω,t), which is expressed as follows:

[0024]

[0025] In the formula, Δt is the time interval between adjacent sampling points;

[0026] S7: Get the time-frequency transform result TE(ω,t) of x(t), which is expressed as follows:

[0027] TE(ω,t)=TEO(ω,t).*W(ω,t),

[0028] Where (.*) indicates multiplication of corresponding index elements.

[0029] Innovation description: The core innovation of this invention is to construct a time-synchronous compression extraction operator, the operator is as follows,

[0030]

[0031] The extraction operator is implemented through a decaying Gaussian window and a group delay operator. The amplitude spectrum energy in the time direction can be compressed within a certain width range, thereby improving the resolution of the amplitude spectrum, and further improving the accuracy of time-frequency analysis. The amplitude spectrum of a single frequency can be kept with the same time resolution, so that comparative analysis of the amplitude spectrum of multiple frequencies can be performed, improving the reliability of time-frequency spectrum analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is a technical flow chart of the present invention;

[0033] Figure 2 This is a comparison chart of time-frequency transformation of different methods of simulating signals. Figure 2 -(a) is the analog signal, Figure 2 -(b) Figure 2 -(c) and Figure 2 -(d) are respectively its wavelet transform, time synchronization compression extraction transform and time-frequency transform results obtained by the method of the present invention;

[0034] Figure 3 To generate the simulated 2D seismic profile and its single-channel data graph, Figure 3 -(a) is a two-dimensional seismic data profile generated, Figure 3 -(b) is the time domain diagram of the signal from which one channel is extracted;

[0035] Figure 4 For Figure 3 The data is processed by wavelet transform and the different frequency component slices of the present invention are compared. Figure 4 -(a) to 4-(d) are the processing results of wavelet transform of 20 Hz, 30 Hz, 40 Hz and 60 Hz frequency components, respectively. Figure 4 -(e) to 4-(h) are respectively the processing results of 20 Hz, 30 Hz, 40 Hz and 60 Hz frequency components of the method of the present invention. DETAILED DESCRIPTION

[0036] The following is further described in conjunction with the accompanying drawings and embodiments. The specific steps of the present invention are as follows:

[0037] Step 1: Extract single-channel data based on seismic data, perform conventional time-frequency transform on it and calculate the time-frequency spectrum. In this example, wavelet transform in conventional time-frequency transform is used. The parameters of wavelet transform are set according to the number of sampling points and sampling frequency of the data. In this example, the number of sampling points is 800 and the sampling frequency is 1000. The wavelet transform result is W(a,b).

[0038]

[0039] In the formula, x(t) represents the time domain signal, a is the scale, b is the time, the total number of scales in the example is set to 500, ()* represents the complex conjugate, ψ(t) is the wavelet basis function, and the complex Morlet wavelet is selected as the wavelet basis function in the example, and its expression is,

[0040]

[0041] where f b is the bandwidth frequency of the complex Morlet wavelet, f c is the center frequency of the complex Morlet wavelet, which can be controlled by b and f c The time-frequency resolution of the wavelet transform can be adjusted. In this example, f b is 1.5, f c is 1.7;

[0042] Step 2: Use tψ(t) as the basis function of wavelet transform to calculate the wavelet transform of the single-channel signal and obtain W tg (a,b),

[0043]

[0044] Where ψ(t) and its parameters should be the same as in step 1;

[0045] Step 3, according to the obtained W(a,b),W tg (a,b) and scale a to get the group delay

[0046]

[0047] Among them, rel() means taking the real part of the complex number, and the obtained group delay And the wavelet transform W(a,b) are all in the scale domain, which needs to be converted to the frequency domain. The calculation method is,

[0048]

[0049] In the formula, ω is the frequency, a is the scale, and f c is the center frequency of the selected wavelet basis function. Through conversion, we can get the wavelet transform W(ω,t) and group delay in the time-frequency domain

[0050] Step 4: Select a suitable window according to actual needs and generate a Gaussian window one-dimensional array gauss(x).

[0051]

[0052] The maximum value of the Gaussian window is 1, and the attenuation on both sides is 0. The length of the window can be controlled by setting the window length. c is the position of the window center point. In the example, the total length of the window is 21, the window center point is 10, and σ is the Gaussian standard deviation, which can control the shape of the Gaussian window. The value in the example is 8;

[0053] Step 5: According to the Gaussian window gauss(x) and group delay The synchronous compression extraction operator TEO(ω,t) is calculated.

[0054]

[0055] In the formula, c represents the position of the center point of the Gaussian window, and its value is 10 in the example, and Δt is the time interval between adjacent sampling points, and its value is 0.001 in the example;

[0056] Step 6, multiply the corresponding elements of the wavelet transform result W(ω,t) and the extraction operator TEO(ω,t) to obtain the synchronous compression result TE(ω,t),

[0057] TE(ω,t)=TEO(ω,t).*W(ω,t),

[0058] As attached Figure 2 As shown, 2-(a) is an analog signal, 2-(b), 2-(c) and 2-(d) are its wavelet transform, time-synchronous compression and extraction transform and the time-frequency transform result obtained by the method, respectively. It can be seen that the resolution of the wavelet transform result is poor, and the time-synchronous compression and extraction transform compresses the amplitude spectrum in the time direction into a ridge line. The present invention compresses the amplitude spectrum energy in the time direction within a certain width range while taking into account the high time resolution of both low frequency and high frequency, and keeps the amplitude spectrum of a single frequency with the same time resolution, thereby improving the accuracy and reliability of time-frequency spectrum analysis.

[0059] Step 7, extracting seismic data one by one, looping through steps 1 to 6, and storing the obtained results in the data volume in sequence;

[0060] Step 8: You can select different frequency components for slicing and then conduct comparative analysis, as shown in the attached figure. Figure 3 , 4 As shown, attached Figure 3 -(a) is a generated two-dimensional seismic data profile, 3-(b) is a time domain diagram of the signal extracted from one of the channels, Figure 4 For attachment Figure 3 The data is sliced ​​into different frequency components after time-frequency transformation, where the attached Figure 4-(a) to 4-(d) are the processing results of wavelet transform for 20 Hz, 30 Hz, 40 Hz and 60 Hz frequency components. It can be seen that due to the different window lengths of low-frequency and high-frequency components, the time resolution of the low-frequency component slices of wavelet transform is low and fuzzy in the time direction, which is not conducive to comparative analysis with high-frequency component slices. Figure 4 -(e) to 4-(h) are the processing results of 20 Hz, 30 Hz, 40 Hz and 60 Hz frequency components of this method respectively. It can be seen that this method compresses the amplitude spectrum energy in the time direction within a certain width range, improves the resolution in the time direction, thereby improving the accuracy of time-frequency analysis, and the amplitude spectrum of each frequency has the same time resolution, which improves the reliability of comparative analysis of amplitude spectra of multiple frequencies.

Claims

1. A method for time-synchronous compression extraction and transformation of a time-frequency spectrum comprises the following specific steps: Step 1: extract single-channel data based on actual seismic data, perform conventional time-frequency transform on it to calculate the time-frequency spectrum, use the wavelet transform method in time-frequency transform, set the wavelet transform parameters according to the number of sampling points and sampling frequency of the data, and obtain the wavelet transform result W(a,b). In the formula, x(t) represents the time domain signal, a is the scale of the wavelet transform, b is its time, ()* represents the complex conjugate of the function, ψ(t) is the wavelet basis function, which should be a complex wavelet function. In practical applications, the complex Morlet wavelet is selected as the wavelet basis function, and its expression is, Where i is a complex unit, f b is the bandwidth frequency of the complex Morlet wavelet, f c is the center frequency of the complex Morlet wavelet, which can be controlled by b and f c To adjust the time-frequency resolution of wavelet transform; Step 2: Use tψ(t) as the basis function and substitute it into the wavelet transform formula to obtain W tg (a,b), Where ψ(t) and its parameters should be the same as in step 1; Step 3: Based on the obtained W(a,b), W tg (a,b) and scale a to get the group delay Among them, rel() means taking the real part of the complex number, and the obtained group delay The wavelet transform W(a,b) obtained in step 1 is in the scale domain and needs to be converted to the frequency domain. The calculation method is: Where ω is the frequency, a is the scale, and f c is the center frequency of the selected wavelet basis function. By conversion, we can get the wavelet transform W(ω,t) and group delay in the time-frequency domain. Where t is the b mentioned in the previous step; Step 4: Select a suitable window according to actual needs and generate a Gaussian window one-dimensional array gauss(x). The maximum value of the Gaussian window is 1, and the attenuation on both sides is 0. The length of the window can be controlled by changing c, where c is the position of the center point of the window. The total length of the window is 2c+1, and the value of x is an integer from 1 to 2c+1. σ is the Gaussian standard deviation, which can control the shape of the Gaussian window. Step 5: According to the Gaussian window gauss(x) and group delay The synchronous compression extraction operator TEO(ω,t) is calculated, and its discrete calculation method is: In the formula, Δt is the time interval between adjacent sampling points; Step 6, multiply the corresponding elements of the wavelet transform result W(ω,t) and the extraction operator TEO(ω,t) to obtain the synchronous compression result TE(ω,t), TE(ω,t)=TEO(ω,t).*W(ω,t), In the formula, (.*) indicates the multiplication of corresponding elements; Step 7, extracting seismic data one by one, looping through steps 1 to 6, and storing the obtained results in the data volume in sequence; Step 8: Select different frequency components for slicing and then conduct comparative analysis.

Images