Seismic data time-frequency analysis method and system based on time sparse S transformation

Abstract

The invention discloses a seismic data time-frequency analysis method and system based on time sparse S transformation, and the method comprises the steps: collecting seismic data, carrying out the time-frequency analysis of the collected seismic data through S transformation, and obtaining a time-frequency spectrum containing the time-frequency characteristics of the seismic data; based on sampling duration and sampling rate characteristics of the seismic data, obtaining a window function of each frequency of the seismic data by using S transformation time-frequency analysis; and constructing a Toeplitz matrix based on a window function, carrying out reverse convolution along the time direction by using the time-frequency spectrum and the Toeplitz matrix, and traversing all frequencies to obtain a seismic data analysis result. The method is applied to time-frequency analysis of seismic data for the first time, a good effect is achieved in multi-scale analysis of two-dimensional actual seismic data, the time resolution of time-frequency analysis of the seismic data is remarkably improved through de-windowing processing, and thin-layer identification of transient signals and seismic data is more accurate especially in a high-frequency part.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic analysis, and in particular relates to a seismic data time-frequency analysis method and system based on time sparse S transform. Background Art

[0002] Time-frequency analysis is one of the important methods for processing and analyzing non-steady-state signals in seismic exploration. It plays an important role in thin-layer reservoir identification, transient signal analysis, and frequency decomposition of seismic data. In the time-frequency spectrum, the frequency components of seismic data at each moment can be determined, which can be used for thin-layer analysis and hydrocarbon detection in seismic data. The resolution of existing window-based time-frequency analysis methods is limited by the uncertainty principle, resulting in the mutual restriction of time resolution and frequency resolution, which cannot reach the highest resolution. Time-frequency analysis relies on Fourier transform, such as windowing the time domain seismic data and then performing Fourier transform, i.e. short-time Fourier transform, or windowing the frequency domain and then performing inverse Fourier transform, i.e. short-frequency Fourier transform (SFFT). However, no matter which transform implementation strategy is used, the window size is fixed, and it has poor adaptability to seismic data with rich low and high frequencies. Therefore, considering the different periods of high and low frequency time signals, wavelet transform (WT) and S transform (ST) with multi-resolution analysis are proposed. Among them, the wavelet basis function used by wavelet transform is used to decompose seismic data, and for different wavelet basis functions, the corresponding relationship between scale and frequency is not fixed. Therefore, the fixed basis function is Morlet wavelet (i.e. fixed center frequency), and the scale of wavelet transform is replaced by frequency to obtain S transform. Similar to wavelet transform, S transform has low time resolution and high frequency resolution at low frequencies, high time resolution and low frequency resolution at high frequencies, and has multi-resolution characteristics.

[0003] Because the window time-frequency analysis method is affected by the uncertainty principle, the time resolution and frequency resolution cannot reach the highest at the same time. Therefore, many scholars assume that the seismic data to be analyzed is a multi-component signal and propose a synchronous squeezing time-frequency analysis method (SST). This method uses the time-frequency spectrum of the window time-frequency analysis method to obtain the instantaneous frequency, and then redistributes the signal within a certain frequency band to the vicinity of the instantaneous frequency to obtain a time-frequency analysis method with high frequency resolution. In order to improve the accuracy of the instantaneous frequency method, high-order synchronous squeezing transform (HighOrder-SST) and multi-step squeezing transform (MSST) are proposed; and in order to obtain a clearer time-frequency spectrum, synchronous extraction transform (SET) is proposed; similar to the above methods, sparse inversion time-frequency analysis methods use part of the Fourier transform matrix as the observation matrix, assume the sparsity of the time-frequency spectrum, and use the sparse inversion strategy to obtain a high-resolution time-frequency spectrum. This type of method is similar to the synchronous squeezing method, and the frequency of the obtained time-frequency spectrum is between the average frequency and the instantaneous frequency. When the seismic data fully meets the assumption of multi-component signals, the adaptive seismic data decomposition strategy and Hilbert-Huang transform can also be used to obtain a high-frequency resolution time-frequency spectrum.

[0004] Although the existing S-transform or window-based time-frequency analysis technology has certain advantages in processing the time-frequency characteristics of seismic data, it also has certain defects and problems that need to be solved. The time resolution and frequency resolution of the S-transform are limited by the uncertainty principle, which means that it is impossible to obtain extremely high time resolution and frequency resolution at the same time; in the S-transform, the seismic data of each frequency component is processed using a fixed-size window, resulting in low time resolution at low frequencies and low frequency resolution at high frequencies; because the S-transform uses multiple filters for multi-band filtering, the window effect may cause distortion of the time-frequency spectrum, affecting the accurate analysis of seismic data. Summary of the invention

[0005] The present invention aims to solve the deficiencies of the prior art and provides the following solutions:

[0006] A method for time-frequency analysis of seismic data based on time sparse S transform comprises the following steps:

[0007] Collecting seismic data, applying S transform to the collected seismic data for time-frequency analysis, and obtaining a time-frequency spectrum containing the time-frequency characteristics of the seismic data;

[0008] Based on the sampling duration and sampling rate characteristics of the seismic data, using S-transform time-frequency analysis to obtain a window function for each frequency of the seismic data;

[0009] A Toeplitz matrix is ​​constructed based on the window function, and the time-frequency spectrum and the Toeplitz matrix are used for reverse convolution along the time direction. After traversing all frequencies, the seismic data analysis result is obtained.

[0010] Preferably, the method for obtaining the time-frequency spectrum comprises:

[0011] The seismic data Windowing is performed to obtain the windowed sub-signal

[0012] According to the windowed sub-signal Get the spectrum

[0013] x i ≈ΓF -1 y i

[0014] Γ=[I|O]

[0015]

[0016] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, WN represents the N-1th row and N-1th column; represents the real number space; represents complex space;

[0017] Let A=ΓF -1 , and obtain the sparse inversion equation at each moment:

[0018]

[0019] Among them, min means minimization, J represents the objective function, G represents the window function, s i represents the i-th sub-signal, λ represents the regularization operator;

[0020] The sparse inversion equation is solved to obtain the time-frequency spectrum.

[0021] Preferably, the method for obtaining the window function includes:

[0022] For each frequency index k, the fixed time window is:

[0023]

[0024] Where w(t) represents the time window, f krepresents the frequency corresponding to the index k, f represents the frequency, and t represents the sampling duration of the seismic data;

[0025] The envelope of the time window is:

[0026]

[0027] in, represents the time window envelope;

[0028] The window function is calculated through the envelope according to the sampling duration and the sampling rate characteristics.

[0029] Preferably, the method for obtaining the seismic data analysis result includes:

[0030] First, a Toeplitz matrix is ​​constructed based on the window function: the time-frequency spectrum is approximately the convolution of the seismic data amplitude with the envelope of the window function. For each frequency index k of the time-frequency spectrum, there exists:

[0031]

[0032] Among them, S(f k , τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k , τ) takes the absolute value; S(f k , τ) matrix, we get:

[0033]

[0034] Among them, S k represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data,

[0035] Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum.

[0036]

[0037] in, Express minimize;

[0038] Solve the high-resolution time-frequency spectrum The seismic data analysis result is obtained.

[0039] The present invention also provides a time-frequency analysis system for seismic data based on time sparse S transform, wherein the analysis system applies any of the above-mentioned analysis methods and comprises: a time-frequency spectrum calculation module, a window function calculation module and a data analysis module;

[0040] The time-frequency spectrum calculation module is used to collect seismic data, apply S transformation to the collected seismic data for time-frequency analysis, and obtain a time-frequency spectrum containing the time-frequency characteristics of the seismic data;

[0041] The window function calculation module obtains the window function of each frequency of the seismic data by using S-transform time-frequency analysis based on the sampling duration and sampling rate characteristics of the seismic data;

[0042] The data analysis module constructs a Toeplitz matrix based on the window function, and uses the time-frequency spectrum and the Toeplitz matrix to perform reverse convolution along the time direction, and obtains seismic data analysis results after traversing all frequencies.

[0043] Preferably, the workflow of the time-frequency spectrum calculation module includes:

[0044] The seismic data Windowing is performed to obtain the windowed sub-signal

[0045] According to the windowed sub-signal Get the spectrum

[0046] x i ≈ΓF -1 y i

[0047] Γ=[I|O]

[0048]

[0049] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, WN represents the N-1th row and N-1th column; represents the real number space; represents complex space;

[0050] Let A=ΓF -1 , and obtain the sparse inversion equation at each moment:

[0051]

[0052] Among them, min means minimization, J represents the objective function, G represents the window function, s i represents the i-th sub-signal, λ represents the regularization operator;

[0053] The sparse inversion equation is solved to obtain the time-frequency spectrum.

[0054] Preferably, the workflow of the window function calculation module includes:

[0055] For each frequency index k, the fixed time window is:

[0056]

[0057] Where w(t) represents the time window, f k represents the frequency corresponding to the index k, f represents the frequency, and t represents the sampling duration of the seismic data;

[0058] The envelope of the time window is:

[0059]

[0060] in, represents the time window envelope;

[0061] The window function is calculated through the envelope according to the sampling duration and the sampling rate characteristics.

[0062] Preferably, the workload flow of the data analysis module includes:

[0063] First, a Toeplitz matrix is ​​constructed based on the window function: the time-frequency spectrum is approximately the convolution of the seismic data amplitude with the envelope of the window function. For each frequency index k of the time-frequency spectrum, there exists:

[0064]

[0065] Among them, S(f k , τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k , τ) takes the absolute value; S(f k , τ) matrix, we get:

[0066]

[0067] Among them, S k represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data,

[0068] Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum.

[0069]

[0070] in, Express minimize;

[0071] Solve the high-resolution time-frequency spectrum The seismic data analysis result is obtained.

[0072] Compared with the prior art, the present invention has the following beneficial effects:

[0073] The present invention is a time-frequency analysis method based on time sparse S transform, which is applied to the time-frequency analysis of seismic data for the first time and has achieved good results in the multi-scale analysis of two-dimensional actual seismic data; through "de-windowing" processing, the time resolution of the time-frequency analysis of seismic data is significantly improved, especially in the high-frequency part, making the thin layer identification of transient signals and seismic data more accurate; the present invention introduces a novel sparse inversion framework for the first time, using non-negative penalty terms and sparse terms to jointly constrain, and solve the optimization problem through the alternating direction multiplier method (ADMM), thereby enhancing the resolvability of the time-frequency spectrum of seismic data; it is not only suitable for the frequency decomposition of seismic data, but also for the analysis of transient signals, and can be combined with other time-frequency analysis methods to provide a variety of high time resolution processing strategies. The present invention can realize modular processing and has a high promotion value. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0075] Figure 1 A schematic diagram of a method flow of an embodiment of the present invention;

[0076] Figure 2 is a flowchart of a method according to an embodiment of the present invention;

[0077] Figure 3: is the time-frequency spectrum obtained by synthesizing seismic data and applying different regularization parameters in the embodiment of the present invention, wherein a is single-channel synthetic test data, b is the time-frequency spectrum obtained based on s-transform, c is the time-frequency spectrum obtained by applying time sparse S-transform when the regularization parameter is 0.05, d is the time-frequency spectrum obtained by applying time sparse S-transform when the regularization parameter is 0.1, e is the time-frequency spectrum obtained by applying time sparse S-transform when the regularization parameter is 0.3, and f is the time-frequency spectrum obtained by applying time sparse S-transform when the regularization parameter is 0.5;

[0078] Figure 4 : The time-frequency spectra obtained by applying different methods to the synthetic seismic data in the embodiment of the present invention, wherein a is single-channel synthetic test data, b is the time-frequency spectra obtained based on STFT transformation, c is the time-frequency spectra obtained based on ST transformation, d is the time-frequency spectra obtained based on SS-STFT transformation, e is the time-frequency spectra obtained based on TSET transformation, and f is the time-frequency spectra obtained based on TSST transformation;

[0079] Figure 5 The single-channel test data and the corresponding time-frequency spectrum synthesized by the random reflection coefficient of the embodiment of the present invention, wherein a is the synthesized single-channel seismic data, b is the time-frequency spectrum obtained based on the STFT transformation, c is the time-frequency spectrum obtained based on the ST transformation, d is the time-frequency spectrum obtained based on the SS-STFT transformation, e is the time-frequency spectrum obtained based on the TSET transformation, and f is the time-frequency spectrum obtained based on the TSST transformation;

[0080] Figure 6 Schematic diagram of reflectivity and synthetic seismic record of the Marmousi model of an embodiment of the present invention, wherein a is reflectivity and b is synthetic seismic record;

[0081] Figure 7 It is a frequency division profile diagram of the data in the Marmousi model using S transformation according to an embodiment of the present invention, wherein a is a frequency division profile diagram under 8 Hz, b is a frequency division profile diagram under 16 Hz, c is a frequency division profile diagram under 24 Hz, d is a frequency division profile diagram under 33 Hz, e is a frequency division profile diagram under 41 Hz, f is a frequency division profile diagram under 49 Hz, g is a frequency division profile diagram under 58 Hz, h is a frequency division profile diagram under 66 Hz, and i is a frequency division profile diagram under 74 Hz;

[0082] Figure 8: This is a frequency division profile diagram of the data in the Marmousi model using TSET transformation according to an embodiment of the present invention, a is a frequency division profile diagram under 8 Hz, b is a frequency division profile diagram under 16 Hz, c is a frequency division profile diagram under 24 Hz, d is a frequency division profile diagram under 33 Hz, e is a frequency division profile diagram under 41 Hz, f is a frequency division profile diagram under 49 Hz, g is a frequency division profile diagram under 58 Hz, h is a frequency division profile diagram under 66 Hz, and i is a frequency division profile diagram under 74 Hz;

[0083] Fig. 9 This is a frequency division profile diagram of the data in the Marmousi model using TSST transformation according to an embodiment of the present invention, a is a frequency division profile diagram under 8 Hz, b is a frequency division profile diagram under 16 Hz, c is a frequency division profile diagram under 24 Hz, d is a frequency division profile diagram under 33 Hz, e is a frequency division profile diagram under 41 Hz, f is a frequency division profile diagram under 49 Hz, g is a frequency division profile diagram under 58 Hz, h is a frequency division profile diagram under 66 Hz, and i is a frequency division profile diagram under 74 Hz;

[0084] Fig.10 is actual seismic data of an embodiment of the present invention;

[0085] Fig.11 The embodiment of the present invention is the 29th data and the corresponding ST and TSST time-frequency spectra, where a is single-channel seismic data, b is the ST method time-frequency analysis result, and c is the TSST time-frequency analysis result;

[0086] Fig.12 It is a comparison diagram of the frequency division profiles of different frequencies obtained by the ST method and the TSST method of the embodiment of the present invention, wherein a is the frequency division profile diagram of 17.5 Hz obtained by the ST method, b is the frequency division profile diagram of 17.5 Hz obtained by the TSST method, c is the frequency division profile diagram of 27.5 Hz obtained by the ST method, d is the frequency division profile diagram of 27.5 Hz obtained by the TSST method, e is the frequency division profile diagram of 37.5 Hz obtained by the ST method, f is the frequency division profile diagram of 37.5 Hz obtained by the TSST method, g is the frequency division profile diagram of 47.5 Hz obtained by the ST method, and h is the frequency division profile diagram of 47.5 Hz obtained by the TSST method. DETAILED DESCRIPTION

[0087] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0088] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0089] Embodiment 1

[0090] In this embodiment, firstly, the time domain and frequency domain S transforms are defined as follows according to the Stockwell phase correction theory:

[0091]

[0092] Among them, t and τ represent time, f and α represent frequency, x(t) represents the seismic data to be transformed, X(f) represents the Fourier transform of x(t), i represents the imaginary unit, and e represents the natural number. At low frequencies, the time resolution is low and the frequency resolution is high; at high frequencies, the time resolution is high and the frequency resolution is low. Because of this characteristic, the S transform is more suitable for the analysis of seismic data signals with high and low frequencies. Compared with the short-time Fourier transform, the S transform has no "low-frequency truncation" phenomenon when decomposing low frequencies, and the time resolution of high frequencies is also higher than that of the short-time Fourier transform.

[0093] like Figure 1 , Figure 2 As shown, a time-frequency analysis method for seismic data based on time sparse S transform includes the following steps:

[0094] S1. Collect seismic data, apply S transform to the collected seismic data for time-frequency analysis, and obtain a time-frequency spectrum containing the time-frequency characteristics of the seismic data.

[0095] According to formula (1), the method of obtaining the time-frequency spectrum includes: Windowing is performed to obtain the windowed sub-signal According to the windowed sub-signal Get the spectrum

[0096] x i ≈ΓF -1 y i (3)

[0097] Γ=[I|O] (4)

[0098]

[0099] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, W in row N-1 and column N-1 N ; represents the real number space; represents the complex space; let A = ΓF -1 , and obtain the sparse inversion equation at each moment:

[0100]

[0101] Among them, min means minimization, J represents the objective function, G represents the window function, s i represents the i-th sub-signal, λ represents the regularization operator; solve the sparse inversion equation to obtain the time-frequency spectrum.

[0102] S2. Based on the sampling duration and sampling rate characteristics of the seismic data, the window function of each frequency of the seismic data is obtained by using S transform time-frequency analysis.

[0103] According to formula (2), it can be found that the size of the window corresponding to each frequency component of the S transform is fixed. Then, consider removing the window effect along the time direction of the time spectrum, that is, adding sparse constraints to improve the resolution of the time spectrum in the time direction. Then, the method of obtaining the window function includes: for each frequency index k, the fixed time window can be written as:

[0104]

[0105] Where w(t) represents the time window, f k represents the frequency corresponding to the index k, f represents the frequency, and in this formula, t represents the sampling time of the seismic data; the envelope of the time window is:

[0106]

[0107] in, represents the time window envelope; as can be seen from formula (8), the window function is calculated by the current frequency and time, and the data needs to provide parameters such as sampling interval and length, but the window is independent of the data itself and is only related to the coordinate parameters. Therefore, in actual processing, the window function is usually obtained by S transform of a delta function of the same length as the seismic data. When the input data is a delta function, the seismic data corresponding to each frequency of the time spectrum is exactly the unit impulse response of the window filter, and the corresponding amplitude spectrum is the envelope of the window function of this frequency, that is, the window function is calculated by the envelope according to the sampling time length and sampling rate characteristics.

[0108] S3. Construct a Toeplitz matrix based on the window function, and use the time-frequency spectrum and the Toeplitz matrix to perform reverse convolution along the time direction. After traversing all frequencies, the seismic data analysis results are obtained.

[0109] The method for obtaining the seismic data analysis result includes: Formula (8) shows that the basis function of the S transform is the cosine / sine signal under Gaussian envelope modulation. First, the Toeplitz matrix is ​​constructed based on the window function: The time-frequency spectrum can be approximated as the convolution of the seismic data amplitude with the window function envelope. For each frequency index k of the time-frequency spectrum, there exists:

[0110]

[0111] Among them, S(f k , τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k , τ) takes the absolute value; S(f k , τ) matrix, we get:

[0112]

[0113] Among them, S k represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data, Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum.

[0114]

[0115] Among them, λ k represents the regularization parameter, when λ k When is larger, the high-resolution time-frequency spectrum of formula (11) The sparser the data, the higher the temporal resolution. Using the optimal convex approximation (l1 norm ||·||1) of the l0 norm (||·||0) to replace the regularization term in equation (11), we get:

[0116]

[0117] Since the time spectrum x * ≥0, it is necessary to add a non-negative penalty term to formula (12). When the weight is less than zero, a penalty is given to make the function value larger, which does not meet the "minimization" principle and forces the optimization to be performed in the feasible domain where the weight is greater than or equal to 0:

[0118]

[0119] in, Express Minimize; It is worth noting that formula (13) only needs to input the amplitude spectrum and the window function to obtain the high time resolution amplitude spectrum, which belongs to the secondary processing of the time-frequency spectrum. In addition to the proposed S transform, this idea is applicable to the secondary processing of most window time-frequency analysis.

[0120] Solving for high-resolution time-frequency spectrogram Get the seismic data analysis results.

[0121] Specifically, the solution method of the optimization formula (13) is as follows: Formula (13) belongs to a constrained optimization problem, and the external penalty function method needs to be considered. By applying a quadratic external penalty constraint, Formula (13) is modified to an unconstrained optimization problem. For each frequency index k:

[0122]

[0123] The role of the external penalty function is to punish solutions that do not satisfy the constraints. Usually, the β value is large. When every element in is less than 0, will be larger when When every element in is greater than or equal to 0, the constraint is 0. Formula (14) can be simplified as:

[0124]

[0125] Because the amplitude spectrum of each frequency in the time spectrum is not necessarily the same, the regularization parameter λ is k Set at k ∝||S k || ∞ , that is, when solving equation (15), the adaptive regularization parameter of the time domain seismic data corresponding to each frequency can ensure the stability of the solution. This embodiment solves equation (15) based on the alternating direction multiplier method (ADMM) to obtain the seismic data analysis result.

[0126] Embodiment 2

[0127] In this embodiment, two tests are performed to demonstrate the effectiveness of the present invention.

[0128] Test 1:

[0129] From the theoretical part, we know that the data resolution after time-frequency analysis is limited by the regularization parameter, so different regularization parameters are given to test the algorithm. This model uses seismic wavelets with different main frequencies to synthesize non-steady-state seismic data. The test parameters are: k =c||S k || ∞, c = 0.05, 0.1, 0.3, 0.5, the algorithm iteration number is 100 times, and the window function is obtained by S transform of the same length function as the seismic data. The model uses zero-phase Ricker wavelet with a main frequency of 25-45Hz and the same positive polarity reflection coefficient to convolve the synthetic seismic records from shallow to deep, such as Figure 3 As shown in (a), Figure 3 (b) is the time-frequency analysis result of S transform. The spectrum data is triangular from low frequency to high frequency, which is consistent with the characteristics of low time resolution at low frequency and high time resolution at high frequency. Figure 3 (c)-3(f) are the processing results of the time sparse S transform with different regularization parameters. As the regularization parameter increases, the sparsity of the time-frequency spectrum gradually increases. When c = 0.5, the resolution requirement can be met.

[0130] This embodiment tests the resolution of different time-frequency analysis methods. Therefore, by synthesizing a simple reflection coefficient, convolution of the Ricker wavelet with a main frequency of 40 Hz, such as Figure 4 (a) is used to test the accuracy of various time-frequency analysis methods. The data length is 400ms and the sampling rate is 1ms. The comparison methods are STFT, ST, SS-STFT, TSET, and the proposed TSST. At this time, there is mutual interference of energy groups in both the time and frequency directions of STFT and ST. The originally separated phase axes also interfere with each other in the time-frequency spectrum, reducing the resolution, such as Figure 4 (b) and Figure 4 (c) shows; SS-STFT is Figure 4 (d) The instantaneous frequency is calculated along the frequency direction and the energy is squeezed. Therefore, its time resolution is affected by the focusing of the window time-frequency analysis, and because the energy is squeezed in the frequency direction, the resolution in the time direction is not improved much; SET and TSST, such as Figure 4 (e) and Figure 4 As shown in (f), both have high time resolution. Compared with the TSET method, TSST has fewer mid- and high-frequency artifacts. This is due to the advantage of the high time resolution of the S transform itself. It is worth noting that the time window of TSET is extremely small at this time, and this approach abandons the frequency resolution; while TSST is more like extracting the S transform ( Figure 4 (c) The ridges in the frequency direction retain the frequency resolution of the S-transform itself and effectively improve the resolution in the time direction.

[0131] Since actual seismic data is usually complex and contains noise, this embodiment simulates the random reflection coefficient, the time length is 1s, the Ricker wavelet with a main frequency of 30Hz is convolved, and random noise with an amplitude of 10% of the seismic data is added. The reflection coefficient and the synthetic seismic record are as follows: Figure 5 (a) shows: compared with Figure 4(a), the complexity of the test data is significantly increased. Similar to the test conclusion of the simple model, due to the average effect of the window, STFT and ST are in the form of energy clusters in the time-frequency spectrum, such as Figure 5 (b) and Figure 5 (c) shows that SS-STFT calculates the instantaneous frequency, assuming that the synthetic seismic record is a multi-component signal. This assumption allows the bandwidth of the seismic data to be compressed and has a higher frequency resolution, such as Figure 5 (d) The temporal resolution of TSET and TSST is relatively high, such as Figure 5 (e) and Figure 5 As shown in (f), TSST has a more ideal frequency resolution. This time-frequency analysis method that takes both time resolution and frequency resolution into consideration is also suitable for frequency decomposition of seismic data in addition to thin layer identification and transient signal analysis.

[0132] Multi-band decomposition of seismic data is helpful to understand the structural information at different scales. In this embodiment, the Marmousi model ( Figure 6 ) to perform frequency decomposition, and the obtained frequency profile is as follows Figure 7 , Figure 8 , Fig. 9 As shown, Figure 7 is the frequency profile of the Marmousi model data using the S transform, Figure 8 This is the frequency profile of the Marmousi model data using the TSET transformation. Fig. 9 The frequency-divided profile using TSST transformation for the data in the Marmousi model has low time resolution at low frequencies and high time resolution at high frequencies. The high and low frequencies of the frequency-divided profile of TSET are similar, and some strata intersect and interfere, resulting in messy reflections and no frequency-related characteristics. Based on ST, TSST performs time windowing in the time direction to obtain a time-frequency spectrum with high time resolution. Compared with ST and TSET, the frequency-divided profile obtained by the recommended TSST is closer to the true reflection coefficient in the medium and high frequencies. At the same time, the strata contact relationships reflected by different frequencies are also different, which also illustrates the importance of frequency decomposition in seismic data analysis.

[0133] Test 2

[0134] We selected SMI's open source seismic data for relevant tests. A two-dimensional line of data is as follows: Fig.10 As shown, there are 120 CDPs in total, with a sampling rate of 2 ms. The data are close to being parallel to the formation above 1.8 s, and the reflections are relatively messy below 1.8 s. Fig.11 The 29th channel data and the corresponding ST and TSST time-frequency spectra are used to test the algorithm parameters. The frequency profile of the data is shown in Fig.12The frequency-division sections of ST and TSST both show that the stratigraphic trend of the data is not layered sediments overlying chaotic reflections. Due to the low resolution of the original data, the internal reflections are covered by thick sediments. These reflection features are displayed on the frequency-division section, and the time resolution of the TSST frequency-division section is much higher than that of the ST. This also shows that the frequency-division properties of TSST seismic data can effectively provide a good reference for stratigraphic interpretation.

[0135] Embodiment 3

[0136] In this embodiment, a seismic data time-frequency analysis system based on time sparse S transform includes: a time-frequency spectrum calculation module, a window function calculation module and a data analysis module.

[0137] The time-frequency spectrum calculation module is used to collect seismic data, apply S transform to the collected seismic data for time-frequency analysis, and obtain a time-frequency spectrum containing the time-frequency characteristics of the seismic data.

[0138] The workflow of the time-frequency spectrum calculation module includes: Windowing is performed to obtain the windowed sub-signal According to the windowed sub-signal Get the spectrum

[0139] x i ≈ΓF -1 yi

[0140] Γ=[I|O]

[0141]

[0142] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, W in row N-1 and column N-1 N ; represents the real number space; represents the complex space; let A = ΓF -1 , and obtain the sparse inversion equation at each moment:

[0143]

[0144] Among them, min means minimization, J represents the objective function, G represents the window function, s i represents the i-th sub-signal, λ represents the regularization operator; solve the sparse inversion equation to obtain the time-frequency spectrum.

[0145] The window function calculation module obtains the window function of each frequency of the seismic data by using S-transform time-frequency analysis based on the sampling duration and sampling rate characteristics of the seismic data.

[0146] The workflow of the window function calculation module includes: for each frequency index k, the fixed time window can be written as:

[0147]

[0148] Where w(t) represents the time window, f k represents the frequency corresponding to index k, f represents the frequency, and t represents the sampling duration of seismic data; the envelope of the time window is:

[0149]

[0150] in, Represents the time window envelope; the window function is calculated through the envelope according to the sampling duration and sampling rate characteristics.

[0151] The data analysis module constructs a Toeplitz matrix based on the window function, and uses the time-frequency spectrum and the Toeplitz matrix to perform reverse convolution along the time direction. After traversing all frequencies, the seismic data analysis results are obtained.

[0152] The workload flow of the data analysis module includes: first, construct the Toeplitz matrix based on the window function: the time-frequency spectrum can be approximated as the convolution of the seismic data amplitude with the envelope of the window function. For each frequency index k of the time-frequency spectrum, there exists:

[0153]

[0154] Among them, S(f k , τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k , τ) takes the absolute value; S(f k , τ) matrix, we get:

[0155]

[0156] Among them, Sk represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data, Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum.

[0157]

[0158] in, Express Minimize; solve for high-resolution time-frequency spectrum Get the seismic data analysis results.

[0159] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A method for time-frequency analysis of seismic data based on time sparse S transform, characterized in that: The following steps are involved: Collecting seismic data, applying S transform to the collected seismic data for time-frequency analysis, and obtaining a time-frequency spectrum containing the time-frequency characteristics of the seismic data; Based on the sampling duration and sampling rate characteristics of the seismic data, using S-transform time-frequency analysis to obtain a window function for each frequency of the seismic data; A Toeplitz matrix is ​​constructed based on the window function, and the time-frequency spectrum and the Toeplitz matrix are used for reverse convolution along the time direction. After traversing all frequencies, the seismic data analysis result is obtained.

2. The method for time-frequency analysis of seismic data based on time sparse S transform according to claim 1, characterized in that: The method for obtaining the time-frequency spectrum includes: The seismic data Windowing is performed to obtain the windowed sub-signal According to the windowed sub-signal Get the spectrum x i ≈ΓF -1 and i Γ=[I|O] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, W in row N-1 and column N-1 N ; represents the real number space; represents complex space; Let A=ΓF -1 , and obtain the sparse inversion equation at each moment: Among them, min means minimization, J represents the objective function, G represents the window function, S i represents the i-th sub-signal, λ represents the regularization operator; The sparse inversion equation is solved to obtain the time-frequency spectrum.

3. The method for time-frequency analysis of seismic data based on time sparse S transform according to claim 2, characterized in that: The method for obtaining the window function includes: For each frequency index k, the fixed time window is: Where w(t) represents the time window, f k represents the frequency corresponding to the index k, f represents the frequency, and t represents the sampling duration of the seismic data; The envelope of the time window is: in, represents the time window envelope; The window function is calculated through the envelope according to the sampling duration and the sampling rate characteristics.

4. The method for time-frequency analysis of seismic data based on time sparse S transform according to claim 3, characterized in that: The method for obtaining the seismic data analysis result includes: First, a Toeplitz matrix is ​​constructed based on the window function: the time-frequency spectrum is approximately the convolution of the seismic data amplitude with the envelope of the window function. For each frequency index k of the time-frequency spectrum, there exists: Among them, S(f k ,τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k ,τ) takes the absolute value; S(f k ,τ) matrix, we get: Among them, S k represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data, Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum. in, Express minimize; Solve the high-resolution time-frequency spectrum The seismic data analysis result is obtained.

5. A time-frequency analysis system for seismic data based on time sparse S transform, the analysis system applying the analysis method according to any one of claims 1 to 4, characterized in that: include: Time-frequency spectrum calculation module, window function calculation module and data analysis module; The time-frequency spectrum calculation module is used to collect seismic data, apply S transformation to the collected seismic data for time-frequency analysis, and obtain a time-frequency spectrum containing the time-frequency characteristics of the seismic data; The window function calculation module obtains the window function of each frequency of the seismic data by using S-transform time-frequency analysis based on the sampling duration and sampling rate characteristics of the seismic data; The data analysis module constructs a Toeplitz matrix based on the window function, and uses the time-frequency spectrum and the Toeplitz matrix to perform reverse convolution along the time direction, and obtains seismic data analysis results after traversing all frequencies.

6. A seismic data time-frequency analysis system based on time sparse S transform according to claim 5, characterized in that: The workflow of the time-frequency spectrum calculation module includes: The seismic data Windowing is performed to obtain the windowed sub-signal According to the windowed sub-signal Get the spectrum x i ≈TF -1 and i Γ=[I|O] Where Γ represents the truncated matrix; I represents the identity matrix, O represents the zero matrix, F represents the Fourier transform matrix; W N represents the Fourier basis, W in row N-1 and column N-1 N ; represents the real number space; represents complex space; Let A=ΓF -1 , and obtain the sparse inversion equation at each moment: Among them, min means minimization, J represents the objective function, G represents the window function, s i represents the i-th sub-signal, λ represents the regularization operator; The sparse inversion equation is solved to obtain the time-frequency spectrum.

7. A seismic data time-frequency analysis system based on time sparse S transform according to claim 6, characterized in that: The workflow of the window function calculation module includes: For each frequency index k, the fixed time window is: Where w(t) represents the time window, f k represents the frequency corresponding to the index k, f represents the frequency, and t represents the sampling duration of the seismic data; The envelope of the time window is: in, represents the time window envelope; The window function is calculated through the envelope according to the sampling duration and the sampling rate characteristics.

8. A seismic data time-frequency analysis system based on time sparse S transform according to claim 7, characterized in that: The workload flow of the data analysis module includes: First, a Toeplitz matrix is ​​constructed based on the window function: the time-frequency spectrum is approximately the convolution of the seismic data amplitude with the envelope of the window function. For each frequency index k of the time-frequency spectrum, there exists: Among them, S(f k ,τ) represents the time f k The time-frequency spectrum under frequency, |·| represents the frequency spectrum of S(f k ,τ) takes the absolute value; S(f k ,τ) matrix, we get: Among them, S k represents the time-frequency spectrum vector; W k represents the Toeplitz matrix, represents the absolute value of the earthquake data, Then, the time-frequency spectrum and the Toeplitz matrix are used to perform a reverse convolution along the time direction to obtain a high-resolution time-frequency spectrum. in, Express minimize; Solve the high-resolution time-frequency spectrum The seismic data analysis result is obtained.

Images