Generalized s-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition

Abstract

The application provides a generalized S transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition, comprising the following steps: step 1, inputting an original seismic signal set, and constructing a wavelet library based on K-SVD dictionary learning; step 2, combining the fast complex domain matching pursuit algorithm with the wavelet library constructed based on K-SVD dictionary learning to improve the speed and accuracy of seismic signal decomposition; step 3, obtaining the time-frequency spectrum of independent wavelets based on the generalized S transform, and superimposing the time-frequency spectrum of all independent wavelets to perform time-frequency joint analysis of the high-resolution time-frequency spectrum; step 4, studying the change of reservoir characteristics by using single-frequency attributes based on the high-resolution time-frequency spectrum; and step 5, eliminating the high-resolution time-frequency spectrum representing the change of reservoir characteristics after strong reflection. The generalized S transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition forms a high-resolution time-frequency analysis technology according to the development of actual application requirements and the characteristics of non-stationary signals.

Description

Technical Field

[0001] This invention relates to the field of oilfield development technology, and in particular to a generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition. Background Technology

[0002] Seismic signals are typically non-stationary signals. Time-frequency analysis (TFA) techniques can convert seismic signals from the time domain to the frequency domain, revealing the frequency variation over time. Observing these TFA changes in the frequency domain allows for the identification of anomalous variations at certain times, which often reflect underlying geological structures and oil and gas reservoirs. Currently, TFA techniques are widely used in seismic exploration. However, in practical data processing, a crucial problem to solve is how to efficiently and accurately obtain high-resolution TFA results from the massive amounts of high-density seismic data in western China's Zhunzhong region. This data is essential for analyzing changes in reservoir structure before and after the elimination of strong reflection effects from coal seams in the Zhunzhong area.

[0003] Currently, methods for constructing wavelet libraries based on dictionary learning have been widely applied in image denoising and image compression. With the development of machine learning, research on seismic data denoising based on sparse representation and dictionary learning has received widespread attention in recent years. In 2008, Zhao Tianzi et al. proposed a method for reconstructing seismic signals using a matching pursuit algorithm based on sparse representation in the time-frequency domain, achieving noise removal without compromising the effective wave.

[0004] Conventional matching and tracing algorithms require scanning the entire atom library in each match, but the atom library is overcomplete, and the scanning process is computationally intensive and slow.

[0005] Time-varying and non-stationary characteristics are universal laws governing real-world signals. Traditional Fourier transform time-frequency analysis methods are insufficient to describe the detailed local information of signals. Therefore, it is essential to use relevant time-frequency joint distribution functions and transform methods to analyze and process the original signals. Moreover, the development of practical techniques for analyzing original signals in both the time and frequency domains better meets the needs of real-world applications. Traditional time-frequency analysis methods oversimplify real-world signals, leading many to assume that most signals are stationary. However, many natural signals in real life, especially in the field of geophysical exploration, are non-stationary in their original form.

[0006] Chinese patent application CN201810494025.X discloses a high-precision synchronous extraction method for generalized S-transform time-frequency analysis. First, a generalized S-transform is performed on the signal to obtain the generalized S-transform value, which in turn yields the time-frequency spectrum of the generalized S-transform. Based on the generalized S-transform value, the instantaneous frequency is calculated. Then, a synchronous extraction operator SEO is derived from the instantaneous frequency. The SEO is then used to extract the time-frequency spectrum of the generalized S-transform, retaining only the energy near the time-frequency ridge in the time-frequency plane while discarding all other divergent energy, resulting in the synchronously extracted generalized S-transform value. Finally, this value is converted into the time-frequency spectrum of the synchronously extracted generalized S-transform. This invention allows for flexible adjustment of the window function size according to actual needs, achieving higher time-frequency decomposition accuracy and greater flexibility, significantly improving the time and frequency resolution of the signal.

[0007] Chinese patent application CN201910672552.X discloses a high-precision multi-synchronous compressed generalized S-transform time-frequency analysis method, comprising: S1, inputting the original one-dimensional signal to be analyzed; S2, performing a three-parameter generalized S-transform on the input signal to obtain the generalized S-transform value, and calculating the modulus of the generalized S-transform value to obtain the time spectrum of the generalized S-transform; S3, obtaining preliminary estimates of the instantaneous frequency at each point on the time-frequency surface based on the generalized S-transform value; S4, calculating more accurate N-fold instantaneous frequency estimates at each point on the time-frequency surface through N iterations based on the above instantaneous frequency estimates; S5, using the frequency set on the time spectrum as the center frequency set, squeezing the time-frequency point values ​​corresponding to the instantaneous frequencies in the intervals near each center frequency to the center frequency point to obtain the N-fold synchronous compressed generalized S-transform value; S6, calculating the modulus of the N-fold synchronous compressed generalized S-transform value to obtain the time spectrum of the N-fold synchronous compressed generalized S-transform.

[0008] Chinese patent application CN201710501413.1 discloses a method for seismic wavelet extraction. This method includes: obtaining reflection coefficient sequences and corresponding well-side seismic traces for multiple wells based on well data and seismic data, wherein the reflection coefficient sequences are represented by a reflection coefficient matrix, and the corresponding well-side seismic traces are represented by a well-side seismic trace matrix; establishing a dictionary learning equation for the seismic wavelet matrix based on the reflection coefficient matrix and the well-side seismic trace matrix; establishing an objective function based on the reflection coefficient matrix, the well-side seismic trace matrix, and the seismic wavelet matrix; calculating the seismic wavelet matrix that minimizes the objective function value, which is then used as the final seismic wavelet matrix; and extracting and processing the final seismic wavelet matrix to obtain the average wavelet. This invention employs a dictionary learning method, enabling accurate and efficient seismic wavelet extraction from multiple wells simultaneously.

[0009] The existing technologies described above are significantly different from our invention and have failed to solve the technical problem we want to address. Therefore, we have invented a new generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition. Summary of the Invention

[0010] The purpose of this invention is to provide a generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition, which is based on K-SVD dictionary learning for fast matching and tracing molecular waves, performs time-frequency analysis based on generalized S-transform, and uses single-frequency properties to study the changes in reservoir characteristics before and after eliminating strong reflections.

[0011] The objective of this invention can be achieved through the following technical measures: a generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition, which includes:

[0012] Step 1: Input the original seismic signal set and construct the wavelet library based on K-SVD dictionary learning;

[0013] Step 2: Based on the wavelet library constructed by K-SVD dictionary learning, combined with the fast complex domain matching pursuit algorithm, the speed and accuracy of seismic signal decomposition are improved;

[0014] Step 3: Obtain the time spectrum of the independent wavelets based on the generalized S-transform, and superimpose the time spectra of all independent wavelets to perform time-frequency joint analysis with high resolution time spectrum;

[0015] Step 4: Based on the high-resolution time spectrum, study the changes in reservoir characteristics using single-frequency attributes;

[0016] Step 5: Eliminate strong reflections and characterize the reservoir features at high resolution using spectral analysis.

[0017] The objective of this invention can also be achieved through the following technical measures:

[0018] In step 1, based on the seismic waveform characteristics of the target layer, appropriate feature wavelets are extracted using K-SVD dictionary learning technology to form a wavelet library.

[0019] Step 1 includes:

[0020] Step 1a: Initialize the sample Y. The original sample Y belongs to the sample space R, which has dimensions of m rows and n columns. k represents any k column vectors drawn from space R as the initial dictionary. m×n Randomly select k column vectors, or select the first k column vectors {d1, d2, ..., d} of its left singular matrix. k} as the atoms of the initial dictionary, resulting in dictionary D. (0) ∈R m×kLet j = 0, and repeat steps 1b-1c below until the specified number of iterations is reached, or the error converges to the specified error.

[0021] Step 1b: Perform sparse coding. Using the dictionary D(j) obtained in step 1a and the sparse coding, obtain X(j)∈R. k×N ;

[0022] Step 1c, perform dictionary update: update dictionary D(j) column by column, column d of the dictionary. k ∈{d1,d2,L d k}

[0023] Step 1c includes:

[0024] Step c1, when updating d k ∈{d1,d2,L d k When calculating the error matrix...

[0025] Step c2: Extract the k-th row vector of the sparse matrix. A set of non-zero indices

[0026] Step c3, from E k Take out the corresponding ω k For columns that are not zero, we get E′. k ;

[0027] Step c4, for E′ k Singular value decomposition Take the first column of U to update the kth column of the dictionary, i.e., d k ==U(·,1), let get Then, update it accordingly to the original

[0028] In step c5, j = j + 1, the process returns to step 1a.

[0029] In step 2, the complex domain fast matching pursuit decomposition and reconstruction of the wavelet library is based on the K-SVD dictionary learning. The seismic signal is decomposed using the complex domain fast matching pursuit technique based on the K-SVD wavelet library.

[0030] In step 2, the Fast Complex Domain Matching Pursuit (CFMP) method is used to improve computational efficiency. The CFMP method requires obtaining complex seismic records based on real seismic records. Prior information, including the time, instantaneous frequency, and instantaneous phase at the maximum amplitude envelope, is obtained from the complex signal. Then, under the constraint of the prior information, a small portion of the atomic library is scanned to obtain the best match. After the complex signal is reconstructed, the result is returned to the real domain, thereby reducing the amount of computation in the matching process and improving the algorithm efficiency to a certain extent.

[0031] In step 2, the formulas for calculating the amplitude envelope, instantaneous frequency, and instantaneous phase of the complex signal are as follows:

[0032]

[0033]

[0034]

[0035] Where f(t) is the real part of the signal, f * (t) represents the imaginary part of the signal; based on this prior information, the matching atom of the response is quickly found in the dictionary.

[0036] In step 2, let H be the Hilbert space, and G = {g r |||g r ||=1,r∈M} are basic functions in H, and each g r All signals have been standardized, and M represents the total number of time-frequency atoms in the time-frequency atom library. Assuming a complex signal F ∈ H, the goal of the CFMP algorithm is to represent F as a linear combination of complex time-frequency atoms selected from G. The complex signal is obtained through a Hilbert transform.

[0037] F(t)=f(t)+i·H(f(t)) (4)

[0038] In the formula: f(t) is the real signal; the imaginary part H(f(t)) is the Hilbert transform of the real signal, which is essentially a 90° phase shift of the real part.

[0039] In step 2, assuming g(t) is a time-frequency atom and H(g(t)) represents the Hilbert transform of the time-frequency atom, the expression for the complex time-frequency atom is:

[0040] G(t)=g(t)+i·H(g(t)) (5).

[0041] In step 2, the matching trace starts from the initial residual R. 0 Starting with F = F, first calculate the instantaneous amplitude envelope, instantaneous frequency, and instantaneous phase of the residual. Then, assign the time, instantaneous frequency, and instantaneous phase at the maximum value of the amplitude envelope to t. j f j and Through n matching calculations, we obtain:

[0042]

[0043] R n F is a complex signal; where... It is a complex time-frequency atom, and R n F is the residual quantity. Let represent the inner product of the signal and the residual; therefore, we have:

[0044]

[0045] In order to make ||R n+1 If F|| is minimized, then the selected complex wavelet is required. Satisfaction makes The optimal selection criterion for the largest, or best, atom is:

[0046]

[0047] The time-frequency atom selected by equation (8) can minimize the residual obtained in each iteration; the best matching wavelet atom can be calculated by equation (8), that is, the optimal amplitude, time shift, main frequency and phase parameters can be calculated.

[0048] In step 3, a series of seismic wavelets characterizing the local features of the target layer are obtained by using a complex domain fast matching pursuit decomposition algorithm based on the K-SVD dictionary learning wavelet library. Then, the time spectrum of each individual wavelet is calculated, and the time spectra of the individual wavelets are superimposed to form a high-resolution time spectrum, which can effectively eliminate the influence of non-stationary signals on the accuracy of video analysis.

[0049] In step 3, the mathematical formula for the generalized S-transform is:

[0050]

[0051] In the formula, gst(τ,f) is the time spectrum obtained by the generalized S-transform, where τ is the time parameter, f is the frequency parameter, x(t) is the time-domain seismic signal, and the parameters λ and p are used to adjust the rate at which the width of the time window function changes with the frequency.

[0052] When the window width is compressed in the time domain, the corresponding bandwidth will naturally be stretched in the frequency domain, and vice versa. According to the Heisenberg uncertainty criterion, time and frequency resolution cannot be optimal at the same time. In order to obtain excellent time resolution, the window width of the window function is reduced, but at this time, excellent frequency resolution cannot be obtained, and vice versa. Here, λ is the time resolution adjustment factor and p is the frequency resolution adjustment factor.

[0053] In step 4, the time spectra of all independent wavelets calculated based on formula (32) are superimposed to form the time spectrum of a certain seismic signal. The three-dimensional data volume formed by extracting the energy of a single frequency from all channels in the three-dimensional data volume is one of the single frequency volumes. Through the above time-frequency analysis, a series of single frequency volumes can be formed. Based on the single frequency volume, the changes in reservoir characteristics can be analyzed by single frequency profile or slice analysis, thereby making reservoir prediction.

[0054] In step 5, during the time-frequency analysis of seismic signals, strong reflection seismic signals, including those from special lithological bodies such as coal seams and igneous rocks, reflect strong energy and affect the time-frequency characteristics of the reservoirs above and below them, reducing the time-frequency resolution. Therefore, when processing such signals, when obtaining the time spectrum of a single wavelet, based on the strong reflection characteristics of special lithological bodies such as coal seams and igneous rocks, strong reflection wavelets can be removed or weakened. The remaining wavelets can then be used to obtain the time spectrum and superimposed to form a high-resolution time spectrum with de-reflection response, thereby improving the accuracy of reservoir prediction.

[0055] The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition in this invention reduces the scanning workload in matching by first constructing a wavelet library using K-SVD dictionary learning, and then using the Complex Domain Fast Matching Pursuit Decomposition (CFMP) method to improve computational efficiency. The CFMP method requires obtaining complex seismic records from real seismic records, extracting prior information (such as the time, instantaneous frequency, and instantaneous phase at the maximum amplitude envelope) from the complex signal, and then scanning a small portion of the wavelet library under the constraints of the prior information to obtain the best match. After reconstructing the complex signal, the result is returned to the real domain, thereby reducing the computational workload in matching to a certain extent and improving algorithm efficiency. Based on the development of practical application needs and targeting the characteristics of non-stationary signals, this invention forms a high-resolution time-frequency analysis technique.

[0056] The key point of this invention is based on the K-SVD dictionary learning algorithm. By learning from actual seismic data, a wavelet library is established. Then, based on the dictionary learning of the wavelet library, the seismic signal is decomposed using the complex domain fast matching pursuit decomposition algorithm. The seismic signal is decomposed into independent wavelets that can express the local characteristics of the seismic signal. The time spectrum of the independent wavelets is then obtained using the generalized S-transform. Finally, the time spectra of all independent wavelets are superimposed to form a time-frequency joint analysis technique with high resolution time spectrum.

[0057] Compared with conventional methods, the present invention utilizes the above technical solutions and has the following advantages: 1. The present technical solution provides an efficient and reliable K-SVD dictionary learning algorithm, which learns local features of seismic signals based on the target layer and constructs a wavelet library. Compared with the redundant wavelet library constructed using analytical expressions in the traditional method, the wavelets extracted using K-SVD dictionary learning have higher accuracy and clearer geological meaning; 2. The present technical solution provides a wavelet library constructed based on K-SVD dictionary learning combined with a fast complex domain matching pursuit algorithm, laying the foundation for efficient and rapid seismic wavelet decomposition and accurate feature extraction; 3. The generalized S-transform time-frequency analysis method provided by the present technical solution is based on seismic signal decomposition and can effectively eliminate the influence of non-stationary signals, improve the accuracy of time-frequency analysis, and obtain high-resolution time-frequency spectra; 4. When strong reflection signals (such as coal seam reflection) affect the reservoir time-frequency characteristics, the present technical solution can combine the coal seam waveform characteristics to selectively remove the reflection waveforms of strong reflection layers (such as coal seams) and use other reflection waveforms to carry out time-frequency attribute analysis, thereby eliminating the influence of strong reflection layers (such as coal seams). Attached Figure Description

[0058] Figure 1 This is a flowchart of a specific embodiment of the generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition of the present invention;

[0059] Figure 2 This is a schematic diagram of single-frequency attribute analysis in a specific embodiment of the present invention;

[0060] Figure 3 This is a schematic diagram of the original seismic profile in a specific embodiment of the present invention;

[0061] Figure 4 This is a schematic diagram of a high-resolution single-frequency attribute profile (30Hz) corresponding to the original seismic profile in a specific embodiment of the present invention;

[0062] Figure 5 This is a schematic diagram of a high-resolution single-frequency attribute profile (40Hz) corresponding to the original seismic profile in a specific embodiment of the present invention;

[0063] Figure 6 This is a schematic diagram of a high-resolution single-frequency attribute profile (50Hz) corresponding to the original seismic profile in a specific embodiment of the present invention;

[0064] Figure 7 This is a schematic diagram of a high-resolution single-frequency attribute profile (30Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile in a specific embodiment of the present invention.

[0065] Figure 8This is a schematic diagram of a high-resolution single-frequency attribute profile (40Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile in a specific embodiment of the present invention.

[0066] Figure 9 This is a schematic diagram of a high-resolution single-frequency attribute profile (50Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile in a specific embodiment of the present invention.

[0067] Figure 10 This is a single-frequency attribute planar plot showing the change in reservoir characteristics before and after eliminating strong reflections at high resolution in a specific embodiment of the present invention. Detailed Implementation

[0068] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0069] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0070] Time-varying and non-stationary characteristics are universal laws governing real-world signals. Traditional Fourier transform time-frequency analysis methods are insufficient to describe the detailed local information of signals. Therefore, it is essential to analyze and process the original signal using relevant time-frequency joint distribution functions and transform methods. Furthermore, the development of practical techniques for analyzing the original signal in both the time and frequency domains better meets the needs of real-world applications. Traditional time-frequency analysis methods oversimplify real-world signals, leading many to assume most signals are stationary. However, many natural signals in real life, especially in geophysical exploration, are non-stationary. The purpose of this invention is to develop a high-resolution time-frequency analysis technique tailored to the characteristics of non-stationary signals, based on the evolving needs of practical applications.

[0071] The generalized S-transform time-frequency analysis technique based on K-SVD dictionary learning of this invention mainly targets the seismic waveform characteristics of the target layer. It utilizes K-SVD dictionary learning to extract wavelets that reflect the local waveform characteristics of the target layer, forming a wavelet library. A complex-domain fast matching pursuit decomposition algorithm based on the K-SVD dictionary-learned wavelet library can obtain a series of seismic wavelets characterizing the local features of the target layer. The time spectrum of each individual wavelet is then calculated using the generalized S-transform, and the time spectra of the individual wavelets are superimposed to form a high-resolution time spectrum, effectively eliminating the influence of non-stationary signals on the accuracy of time-frequency analysis. To address the impact of strong reflections (such as coal seams) on the accuracy of time-frequency analysis, when calculating the time spectrum of individual wavelets, coal seam wavelets can be removed based on their reflection characteristics. The remaining wavelets are then used to calculate their time spectra, which are then superimposed to form a high-resolution time spectrum without the coal seam response. Actual data processing has fully verified the effectiveness and feasibility of the technical solution described in this invention.

[0072] The key point of this invention is based on the K-SVD dictionary learning algorithm. By learning from actual seismic data, a wavelet library is established. Then, based on the dictionary learning of the wavelet library, the seismic signal is decomposed using the complex domain fast matching pursuit decomposition algorithm. The seismic signal is decomposed into independent wavelets that can express the local characteristics of the seismic signal. The time spectrum of the independent wavelets is then obtained using the generalized S-transform. Finally, the time spectra of all independent wavelets are superimposed to form a time-frequency joint analysis technique with high resolution time spectrum.

[0073] Example 1

[0074] In a specific embodiment 1 of this invention, the present invention provides a wavelet library construction method based on K-SVD dictionary learning. By learning wavelet characteristics from actual seismic data, the redundancy of the wavelet library is effectively reduced, and the constructed wavelet library is closer to the actual data, truly reflecting the local characteristics of the actual data. Based on this, fast complex domain matching pursuit is used to effectively improve the accuracy and speed of wavelet decomposition. Then, the time spectrum of independent wavelets is obtained using the generalized S-transform, and the time spectra of all independent wavelets are superimposed to form a time-frequency joint analysis technique with high-resolution time spectra. The method includes:

[0075] (1) Wavelet library construction technology based on K-SVD dictionary learning;

[0076] Primarily targeting the seismic waveform characteristics of the target layer, K-SVD dictionary learning technology is used to extract appropriate feature wavelets and construct a wavelet library; specifically including:

[0077] Step 1a: Initialize the sample Y. The original sample Y belongs to the sample space R, which has dimensions of m rows and n columns. k represents any k column vectors drawn from space R as the initial dictionary. m×nRandomly select k column vectors, or select the first k column vectors {d1, d2, ..., d} of its left singular matrix. k} as the atoms of the initial dictionary, resulting in dictionary D. (0) ∈R m×k Let j = 0, and repeat steps 1b-1c below until the specified number of iterations is reached, or the error converges to the specified error.

[0078] Step 1b: Perform sparse coding. Using the dictionary D(j) obtained in step 1a and the sparse coding, obtain X(j)∈R. k×N ;

[0079] Step 1c, perform dictionary update: update dictionary D(j) column by column, column d of the dictionary. k ∈{d1,d2,L d k}

[0080] Step 1c includes:

[0081] Step c1, when updating d k ∈{d1,d2,L d k When calculating the error matrix...

[0082] Step c2: Extract the k-th row vector of the sparse matrix. The set of non-zero indices

[0083] Step c3, from E k Take out the corresponding ω k For columns that are not zero, we get E′. k ;

[0084] Step c4, for E′ k Perform singular value decomposition to obtain E k =U∑V T Take the first column of U and update the kth column of the dictionary, i.e., d. k =U(·,1), let get Then, update it accordingly to the original

[0085] In step c5, j = j + 1, the process returns to step 1a.

[0086] (2) The wavelet library constructed based on K-SVD dictionary learning, combined with the fast complex domain matching pursuit algorithm, improves the speed and accuracy of seismic signal decomposition;

[0087] Rapid matching pursuit decomposition and reconstruction in the complex domain based on K-SVD dictionary learning wavelet library: Based on the K-SVD wavelet library, the seismic signal is decomposed using rapid matching pursuit technology in the complex domain.

[0088] The formulas for calculating the amplitude envelope, instantaneous frequency, and instantaneous phase of a complex signal are as follows:

[0089]

[0090]

[0091]

[0092] Where f(t) is the real part of the signal, f * (t) represents the imaginary part of the signal; based on this prior information, the matching atom of the response is quickly found in the dictionary.

[0093] Let H be a Hilbert space, and G = {g r |||g r ||=1,r∈M} are basic functions in H, and each g r All signals have been standardized, and M represents the total number of time-frequency atoms in the time-frequency atom library. Assuming a complex signal F ∈ H, the goal of the CFMP algorithm is to represent F as a linear combination of complex time-frequency atoms selected from G. The complex signal is obtained through a Hilbert transform.

[0094] F(t)=f(t)+i·H(f(t)) (4)

[0095] In the formula: f(t) is the real signal; the imaginary part H(f(t)) is the Hilbert transform of the real signal, which is essentially a 90° phase shift of the real part.

[0096] Assuming g(t) is a time-frequency atom, and H(g(t)) represents the Hilbert transform of the time-frequency atom, then the expression for the complex time-frequency atom is:

[0097] G(t)=g(t)+i·H(g(t)) (5).

[0098] Matching and tracking from the initial residual R 0 Starting with F = F, first calculate the instantaneous amplitude envelope, instantaneous frequency, and instantaneous phase of the residual. Then, assign the time, instantaneous frequency, and instantaneous phase at the maximum value of the amplitude envelope to t. j f j and Through n matching calculations, we obtain:

[0099]

[0100] R n F is a complex signal; where... It is a complex time-frequency atom, and R n F is the residual quantity. Let represent the inner product of the signal and the residual; therefore, we have:

[0101]

[0102] In order to make ||R n+1 If F|| is minimized, then the selected complex wavelet is required. Satisfaction makes The optimal selection criterion for the largest, or best, atom is:

[0103]

[0104] The time-frequency atom selected by equation (8) can minimize the residual obtained in each iteration; the best matching wavelet atom can be calculated by equation (8), that is, the optimal amplitude, time shift, main frequency and phase parameters can be calculated.

[0105] (3) The time spectrum of independent wavelets is obtained based on the generalized S-transform, and the time spectra of all independent wavelets are superimposed to form a time-frequency joint analysis technique with high resolution time spectrum;

[0106] A fast matching pursuit decomposition algorithm in the complex domain based on K-SVD dictionary learning wavelet library is used to obtain a series of seismic wavelets that characterize the local features of the target layer. The time spectrum of each individual wavelet is then calculated, and the time spectra of the individual wavelets are superimposed to form a high-resolution time spectrum, which can effectively eliminate the influence of non-stationary signals on the accuracy of video analysis.

[0107] (4) Based on high-resolution time spectrum, the reservoir characteristic changes are studied by using single-frequency attributes; the time spectrum of all independent wavelets calculated based on generalized S-transform is superimposed to form the time spectrum of a certain seismic signal. The three-dimensional data volume formed by extracting the single-frequency energy of all traces in the three-dimensional data volume is one of the single-frequency volumes. Through the above time-frequency analysis, a series of single-frequency volumes can be formed. Based on the single-frequency volume, the changes in reservoir characteristics can be analyzed by single-frequency profile or slice, thereby making reservoir prediction.

[0108] (5) High-resolution time spectrum characterization of reservoir characteristics after eliminating strong reflections: When performing time-frequency analysis of seismic signals, strong reflection seismic signals, including special lithological bodies such as coal seams and igneous rocks, reflect strong energy, which will affect the time-frequency characteristics of the reservoirs above and below them and reduce the time-frequency resolution. Therefore, when processing such signals, when obtaining the time spectrum of a single wavelet, based on the strong reflection characteristics including special lithological bodies such as coal seams and igneous rocks, the strong reflection wavelet can be eliminated or weakened, and the time spectrum can be obtained with the remaining wavelets and superimposed to form a high-resolution time spectrum with the strong reflection response eliminated, thereby improving the accuracy of reservoir prediction.

[0109] Example 2

[0110] In a specific embodiment 2 of the present invention, such as Figure 1 As shown, Figure 1 This is a flowchart of the generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition of the present invention. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition includes the following steps:

[0111] Step 1, K-SVD dictionary learning

[0112] In sparse representation-based methods, the quality of the dictionary directly impacts the sparsity of sparse coding and plays a crucial role in selective signal reconstruction. Generally, overcomplete dictionaries increase the number of atomic bases, exceeding the signal dimension. Therefore, the signal representation in the spatial domain of an overcomplete dictionary is not unique. However, the advantage is that we can adaptively search for the optimal combination of dictionary atoms through learning and training, thus representing the target signal more sparsely. Dictionary learning is currently the most commonly used method for constructing overcomplete dictionaries. It can adaptively learn redundant dictionaries that match the signal characteristics using training data for different types of signals. Thus, during sparse signal representation, the best-reflecting atoms in the trained overcomplete dictionary can be found, and the optimal sparse approximation of the signal can be achieved through linear combination.

[0113] The dictionary learning process can be viewed from the perspective of matrix factorization: given a sample dataset Y, each column represents a sample; the goal of dictionary learning is to decompose the Y matrix into matrices D and X.

[0114] Y≈DX(9)

[0115] The requirements are as follows: X should be as sparse as possible, and each column of D should be a normalized vector. D is called a dictionary, and each column of D is called an atom; X is called the encoding vector, feature, or coefficient matrix; dictionary learning can have three forms of objective function.

[0116] The first form:

[0117]

[0118] Because L0 is difficult to solve, L1 regularization terms are often used as approximations.

[0119] The second form:

[0120]

[0121] ε is the maximum value allowed for the reconstruction error.

[0122] The third form:

[0123] D,X=argmin||Y-DX||2 st.||X||0≤L (12)

[0124] L is a constant, a sparsity constraint parameter, and the three forms above are equivalent to each other.

[0125] Let Y = [y1, y2, ..., y] N ]∈R n×N If N signals are used for dictionary learning training, then the mathematical model for dictionary learning can be expressed as:

[0126]

[0127] In the formula, D is an overcomplete dictionary, and X = [x1, x2, ..., x... N ]∈R k×N Let be the sparse representation coefficients of Y in dictionary D, where y i =Dx i , To express a sparsity constraint on X, you can use ||X||0 or ||X||1, or other l-values. p The norm is used as a constraint, and λ is a constant. For the dual-criteria problem, a trade-off needs to be made between the sparse representation error and the sparsity of the representation coefficients. Since there are two unknown variables in the objective function, the objective function can be solved by iteratively fixing the dictionary D and the sparse coefficients X in sequence, thus reducing the difficulty of the solution. The specific solution process can be divided into two main steps:

[0128] 1) Solving for sparsity coefficients: Assuming dictionary D is known, the above problem becomes a sparse decomposition problem, which can be solved using sparse decomposition algorithms such as Lasso (Least Absolute Shrinkage and Selection Operator) and OMP (Orthogonal Matching Purpose).

[0129]

[0130] 2) Dictionary Training: Assuming the sparse representation coefficients X are known, train an adaptive dictionary D:

[0131]

[0132] In this way, the two steps are continuously updated and iterated, and after multiple cycles, an overcomplete dictionary D that matches the features of signal Y and the corresponding sparse coefficients X can be obtained.

[0133] Step 2, Fast Matching and Pursuit Decomposition in Complex Fields

[0134] Conventional matching pursuit algorithms require scanning the entire atomic library in each match. However, the atomic library is overcomplete, resulting in a large computational burden and slow computation speed. To reduce the scanning workload in matching, the Complex Domain Fast Matching Pursuit Decomposition (CFMP) method can be used to improve computational efficiency. The CFMP method requires obtaining complex seismic records from real seismic records, extracting prior information (such as the time, instantaneous frequency, and instantaneous phase at the maximum amplitude envelope) from the complex signal, and then scanning a small portion of the atomic library under the constraints of the prior information to obtain the best match. After reconstructing the complex signal, the result is returned to the real domain, thereby reducing the computational burden in matching to a certain extent and improving the algorithm efficiency. Let H be a Hilbert space, satisfying... These are the basic functions in H, each g r All have been standardized, and Γ is the index. The primitive function G is a time-frequency atom library, where each g... r Each time-frequency atom is a single element. Through iterative matching calculations on the complex domain signal, each match yields a complex time-frequency atom. Assuming a complex signal F ∈ H, the goal of the CFMP algorithm is to represent F as a linear combination of complex time-frequency atoms selected from G. The complex signal can be obtained through the Hilbert transform:

[0135] F(t)=f(t)+i·H(f(t)) (16)

[0136] In the formula: f(t) is the real signal; the imaginary part H(f(t)) is the Hilbert transform of the real signal, which is essentially a 90° phase shift of the real part.

[0137] Similarly, assuming g(t) is a time-frequency atom, and H(g(t)) represents the Hilbert transform of the time-frequency atom, then the expression for the complex time-frequency atom is:

[0138] G(t)=g(t)+i·H(g(t)) (17)

[0139] Matching and tracking from the initial residual R 0 Starting with F = F, first calculate the instantaneous amplitude envelope, instantaneous frequency, and instantaneous phase of the residual. Then, assign the time, instantaneous frequency, and instantaneous phase at the maximum value of the amplitude envelope to t. j f j and Through n matching calculations, we obtain:

[0140]

[0141] R n F is a complex signal. It is a complex time-frequency atom, and R n F is the residual quantity. This represents the inner product of the signal and the residual. It is easy to derive... Orthogonal to R n+1 F, therefore we have:

[0142]

[0143] In order to make ||R n+1 If F|| is minimized, then the selected complex wavelet... Make The maximum, or optimal, atom selection criterion is:

[0144]

[0145] The time-frequency atom energy selected by equation (20) minimizes the residual obtained in each iteration. When the local properties of the signal are exactly the same as or completely match those of the matching atom, the projection value of the complex signal onto the atom is... If it is a real number, then we have Therefore, the matching condition is improved as follows:

[0146]

[0147] The optimal matching wavelet atom can be calculated using equation (21), which yields the optimal time shift, dominant frequency, phase, and other parameters. To achieve the best matching effect, this paper constructs a time-frequency atom dictionary W within a certain neighborhood of the previously calculated parameters (delay, dominant frequency, and phase) when building the time-frequency atom library. W = {W γ |||W γ ||=1}, where: the index set γ is an element in the set Γ, and the complex time-frequency atoms are all normalized. After applying the matching condition (21), the result of the first matching is obtained as:

[0148]

[0149] Among them, R 1 F is the complex residual from the first iteration, which can be directly used in the next iteration, saving the computational burden of the Hilbert transform. The calculation process is the same as the previous one. By setting the number of iterations or an error threshold, the final result of complex matching pursuit is expressed as:

[0150]

[0151] Due to various noises and other interference factors in the signal, the matching parameters obtained using instantaneous information contain some errors. This paper uses the least squares method in the complex domain to correct the obtained matching parameters, so that the complex Ricker wavelet determined by the matching parameters can more closely approximate the local information of the seismic signal. Let R...n F = A j W γ +R n+1 F, where W γ To determine the obtained complex time-frequency atom using the above instantaneous parameters, A j =|A j |exp(i·θ j ) represents the correction term for the wavelet atom, which includes amplitude correction |A j |With phase correction θ j To ensure that the corrected atomic wavelet best matches the seismic signal, ||R n+1 Find the minimum value of F|| using the least squares method over the complex field:

[0152]

[0153] In the formula: W γ The conjugate transpose of ; ε is the damping factor; I is the identity matrix.

[0154] A is obtained by equation (24). j After obtaining the amplitude and phase correction terms, the matched wavelet atom becomes G. γ =|A j |W γ (t), where the matching parameter is φ j =α j +θ j The phase after correction. The optimal matching wavelet atom can be calculated from equation (24), that is, the optimal time shift, dominant frequency, phase and other parameters are calculated. The prior information of the complex signal is obtained, including: amplitude envelope, instantaneous frequency and instantaneous phase. The calculation formulas for the amplitude envelope, instantaneous frequency and instantaneous phase of the complex signal are as follows:

[0155]

[0156]

[0157]

[0158] Where f(t) is the real part of the signal, f * (t) represents the imaginary part of the signal. The amplitude coefficient after n iterations is:

[0159]

[0160] Step 3, Generalized S-transform

[0161] The STFT transform of signal x(t) is:

[0162]

[0163] The window function used in signal analysis and its width are closely related to the time-frequency resolution of the STFT. According to the Heisenberg criterion, the product of the time width and frequency bandwidth of the Gaussian window function can reach a minimum value of 1 / 2. Therefore, the window function expression is defined as follows:

[0164]

[0165] make

[0166]

[0167] Where λ>0, p>0.

[0168] Therefore, the scaling factor σ is a function of frequency. Formula (32) not only adaptively adjusts the time window function according to the frequency distribution characteristics in the time-frequency domain, but also embodies the idea of ​​multi-resolution analysis. Combining the above derivation, the mathematical formula of the generalized S-transform is obtained:

[0169]

[0170] In the formula, gst(τ,f) is the time spectrum obtained by the generalized S-transform, where τ is the time parameter, f is the frequency parameter, x(t) is the time-domain seismic signal, and the parameters λ and p are used to adjust the rate at which the width of the time window function changes with the frequency.

[0171] When the time window width is compressed in the time domain, the corresponding bandwidth will naturally be stretched in the frequency domain, and vice versa. According to the Heisenberg uncertainty criterion, time and frequency resolution cannot be optimal simultaneously. To achieve excellent time resolution, we should reduce the time window width of the window function, but this will not yield excellent frequency resolution, and vice versa. Therefore, by adding two adjustable variables, λ and p, we obtain another form of GST representation. When λ = 1 and p = 1, the GST transform becomes an S-transform. Like the standard S-transform, the GST transform also possesses lossless and invertible properties. Based on the inverse Fourier transform, the inverse transform of the GST time spectrum can be achieved, as expressed below:

[0172]

[0173] Step 4: Based on the high-resolution time spectrum, study the changes in reservoir characteristics using single-frequency attributes;

[0174] The time spectra of all independent wavelets calculated by (32) are superimposed to form the time spectrum of a certain seismic signal. The three-dimensional data volume formed by extracting the energy of a single frequency of all channels in the three-dimensional data volume is one of the single frequency volumes. Through the above time-frequency analysis, a series of single frequency volumes can be formed. Based on the single frequency volume, the changes in reservoir characteristics can be analyzed by single frequency profile or slice analysis, thereby making reservoir prediction.

[0175] Step 5: Eliminate strong reflections and characterize the reservoir features at high resolution using spectral analysis.

[0176] When performing time-frequency analysis of seismic signals, the strong reflection seismic signals (such as those reflected by special lithological bodies like coal seams and igneous rocks) have high energy, which can affect the time-frequency characteristics of the reservoirs above and below them, reducing the time-frequency resolution. Therefore, when processing such signals, we can eliminate or weaken the strong reflection wavelets based on the strong reflection characteristics (such as those reflected by special lithological bodies like coal seams and igneous rocks) when obtaining the time spectrum of a single wavelet, and use the remaining wavelets to obtain the time spectrum and superimpose them to form a high-resolution time spectrum with the strong reflection response removed, thereby improving the accuracy of reservoir prediction.

[0177] Example 3

[0178] In a specific embodiment 3 of the present invention, the technology of the present invention was applied to actual data from a certain oil field. Figure 2 This is a schematic diagram of single-frequency attribute analysis. Figure 3 This is the original seismic profile. Figure 4 This is the high-resolution single-frequency attribute profile (20Hz) corresponding to the original seismic profile. Figure 5 This is the high-resolution single-frequency attribute profile (40Hz) corresponding to the original seismic profile. Figure 6 This is the high-resolution single-frequency attribute profile (50Hz) corresponding to the original seismic profile. Figure 7 High-resolution single-frequency attribute profile (20Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile. Figure 8 High-resolution single-frequency attribute profile (40Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile. Figure 9 The high-resolution single-frequency attribute profile (50Hz) corresponding to the elimination of strong coal seam reflections from the original seismic profile is shown in the figure. As can be seen from the figure, the single-frequency profile provides a clearer and more reliable characterization of the reservoir after the strong coal seam reflections are eliminated. Figure 10 (a) is the root mean square amplitude plane attribute map extracted along the J1s21 layer from the high-resolution 30Hz single-frequency attribute volume corresponding to the strong coal seam reflection without elimination. Figure 10(b) is the root mean square amplitude plane attribute map extracted along the J1s21 layer from the high-resolution 30Hz single-frequency attribute volume corresponding to the elimination of strong coal seam reflection. Among them, Z3, Z110, Z102, and Z109 are all actual producing wells. Compared with the actual drilling verification, the time-frequency attributes of the second section of the Sangonghe Formation (J1s21) after coal removal can better characterize the spatial distribution of the reservoir.

[0179] The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition in this invention constructs a wavelet library based on K-SVD dictionary learning in seismic data processing, obtains the time spectrum of independent wavelets based on generalized S-transform, and superimposes the time spectra of all independent wavelets to form a time-frequency joint analysis technique with high resolution time spectrum. It also uses single-frequency attributes to study the changes in reservoir characteristics before and after the elimination of strong reflection in coal seams.

[0180] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0181] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition, characterized in that, This generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition includes: Step 1: Input the original seismic signal set and construct the wavelet library based on K-SVD dictionary learning; Step 2: Based on the wavelet library constructed by K-SVD dictionary learning, combined with the fast complex domain matching pursuit algorithm, the speed and accuracy of seismic signal decomposition are improved; Step 3: Obtain the time spectrum of the independent wavelets based on the generalized S-transform, and superimpose the time spectra of all independent wavelets to perform time-frequency joint analysis with high resolution time spectrum. Step 4: Based on the high-resolution time spectrum, study the changes in reservoir characteristics using single-frequency attributes; Step 5: Eliminate strong reflections and characterize the reservoir features at high resolution using spectral analysis. In step 3, a series of seismic wavelets characterizing the local features of the target layer are obtained by using a complex domain fast matching pursuit decomposition algorithm based on the K-SVD dictionary learning wavelet library. Then, the time spectrum of each individual wavelet is calculated, and the time spectra of the individual wavelets are superimposed to form a high-resolution time spectrum, which effectively eliminates the influence of non-stationary signals on the accuracy of time-frequency analysis. In step 3, the mathematical formula for the generalized S-transform is: （32） In the formula, It is the time spectrum obtained by the generalized S-transform, where It is a time parameter. It is a frequency parameter. It is a time-domain seismic signal, parameters and Used to adjust the rate at which the width of the time window function changes with frequency; where It is the time resolution adjustment factor. It is the frequency resolution adjustment factor; In step 4, the time spectra of all independent wavelets calculated based on formula (32) are superimposed to form the time spectrum of a certain seismic signal. The three-dimensional data volume formed by extracting the energy of a single frequency from all channels in the three-dimensional data volume is one of the single frequency volumes. Through the above time-frequency analysis, a series of single frequency volumes are formed. Based on the single frequency volumes, the changes in reservoir characteristics are analyzed by single frequency profiles or slices to make reservoir prediction.

2. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 1, characterized in that, In step 1, based on the seismic waveform characteristics of the target layer, appropriate feature wavelets are extracted using K-SVD dictionary learning technology to form a wavelet library.

3. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 2, characterized in that, Step 1 includes: Step 1a: Initialize the sample Y. The original sample Y belongs to the sample space R, which has dimensions of m rows and n columns. k represents any k column vectors drawn from space R as the initial dictionary. Random selection The column vectors or the first part of its left singular matrix. column vectors As atoms in the initial dictionary, the dictionary is obtained. ;make Repeat steps 1b-1c below until the specified number of iterations is reached, or the error converges to the specified error. Step 1b: Perform sparse coding using the dictionary obtained in step 1a. And sparse coding, to obtain ; Step 1c, update the dictionary: Update the dictionary column by column. , dictionary columns .

4. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 3, characterized in that, Step 1c includes: Step c1, when updating At that time, calculate the error matrix. ; Step c2, extract the sparse matrix of the first... row vectors A set of non-zero indices ; Step c3, from Take out the corresponding For columns that are not zero, we get ; Step c4, for Singular value decomposition ,Pick The first column updates the dictionary's first column. column, i.e. ,make ,get Then, update it accordingly to the original ; Step c5, The process returns to step 1a.

5. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 1, characterized in that, In step 2, the complex domain fast matching pursuit decomposition and reconstruction of the wavelet library is based on the K-SVD dictionary learning. The seismic signal is decomposed using the complex domain fast matching pursuit technique based on the K-SVD wavelet library.

6. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 5, characterized in that, In step 2, the Fast Complex Domain Matching Pursuit (CFMP) method is used to improve computational efficiency. The CFMP method requires obtaining complex seismic records based on real seismic records. Prior information, including the time, instantaneous frequency, and instantaneous phase at the maximum amplitude envelope, is obtained from the complex signal. Then, under the constraint of the prior information, a small portion of the atomic library is scanned to obtain the best match. After the complex signal is reconstructed, the result is returned to the real domain, thereby reducing the amount of computation in the matching process and improving the algorithm efficiency to a certain extent.

7. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 6, characterized in that, In step 2, the formulas for calculating the amplitude envelope, instantaneous frequency, and instantaneous phase of the complex signal are as follows: （1） （2） （3） in: Let be the real part of the signal. The imaginary part of the signal is denoted by ; based on this prior information, the matching atom of the response is quickly found in the dictionary.

8. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 7, characterized in that, In step 2, let It is a Hilbert space, G = {gr ||gr|| = 1, r ∈ M}. The basic functions in the text, each All have been standardized. It is the total number of time-frequency atoms in the time-frequency atom library; assuming a complex signal The purpose of the CFMP algorithm is to... Indicates from The complex signal of the selected complex time-frequency atom linear combination is obtained through Hilbert transform: （4） In the formula: imaginary part The Hilbert transform of a real signal is essentially a 90° phase shift of the real part.

9. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 8, characterized in that, In step 2, it is assumed that For time-frequency atoms, The Hilbert transform of a time-frequency atom represents the expression for a complex time-frequency atom: （5）。 10. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 9, characterized in that, In step 2, the matching trace starts from the initial residual. First, the instantaneous amplitude envelope, instantaneous frequency, and instantaneous phase of the residual are calculated. The time, instantaneous frequency, and instantaneous phase at the maximum value of the amplitude envelope are then set as follows: , and Through n matching calculations, we obtain: （6） Let be the complex residual signal after the nth iteration; where It is a complex time-frequency atom, and In the nth iteration, the best matching complex time-frequency atom is selected from the complex atom library G; This represents the new complex residual signal in the (n+1)th iteration; Let represent the inner product of the signal and the residual; therefore, we have: （7） In order to To minimize this, the selected complex wavelet is required. Satisfaction makes The maximum, or optimal, atom selection criterion is: （8） The time-frequency atom selected by equation (8) minimizes the residual obtained in each iteration; the best matching wavelet atom is calculated by equation (8), that is, the optimal amplitude, time shift, main frequency and phase parameters are calculated.

11. The generalized S-transform time-frequency analysis method based on dictionary learning and matching pursuit decomposition according to claim 1, characterized in that, In step 5, during the time-frequency analysis of seismic signals, strong reflection seismic signals, including those from special lithological bodies such as coal seams and igneous rocks, reflect strong energy and can affect the time-frequency characteristics of reservoirs above and below them, reducing the time-frequency resolution. Therefore, when processing such signals, when calculating the time spectrum of a single wavelet, based on the strong reflection characteristics of special lithological bodies such as coal seams and igneous rocks, strong reflection wavelets are removed or weakened. The remaining wavelets are used to calculate the time spectrum and superimposed to form a high-resolution time spectrum with de-reflection response, thereby improving the accuracy of reservoir prediction.

Images