Sparse constrained least square earthquake time-frequency analysis method and device

Abstract

The embodiment of the invention relates to a sparse constraint least square seismic time-frequency analysis method and device, and the method comprises the steps: carrying out the complex processing of a real seismic channel, and constructing a complex seismic channel; constructing a diagonal matrix with a data weight as a window function as a diagonal line, and initializing a model weight into a unit matrix; constructing an Euler formula kernel matrix formed by different cut-off frequencies and time ranges in a time domain; applying the data weight constraint to the complex seismic trace, and applying the data weight constraint and the model weight constraint to the Euler formula kernel matrix; and constructing an error function among a constrained Euler formula kernel matrix, a model parameter and a constrained complex seismic trace, applying L1 regularization sparse constraint to the model parameter, and performing inversion by applying a greedy iterative shrinkage threshold method to obtain a sparse constrained least square seismic time-frequency analysis result. Therefore, the seismic signal can be efficiently and accurately decomposed to obtain the high-precision time-frequency spectrum.

Description

Technical Field

[0001] The embodiments of the present invention relate to the fields of seismic exploration and signal processing technology, and in particular to a sparse-constrained least-squares seismic time-frequency analysis method and apparatus. Background Technology

[0002] As exploration targets increasingly shift towards thinner, smaller, and unconventional complex oil and gas reservoirs, it becomes necessary to extract attributes reflecting subtle structural and lithological variations from complex seismic wavegroup characteristics. However, existing time-domain conventional seismic data has relatively low resolution, which is insufficient for resolving more detailed and complex reservoir descriptions. This difficulty in achieving precise descriptions of target layers is one of the fundamental problems restricting the current level of fine-grained oil and gas reservoir description.

[0003] To address the problem of low resolution in time-domain seismic data making it difficult to interpret complex and fine structures, spectral decomposition techniques can be used to analyze the frequency-time-space characteristics of seismic signals and extract amplitude anomaly attributes to indicate thin layers or attenuating geological bodies. However, conventional methods are limited by the length of the time window, suffer from window tailing effects, or are affected by energy interference from surrounding seismic events, resulting in either low resolution or high computational costs, which cannot meet the current requirements for high-precision reservoir identification. Therefore, how to more accurately invert the amplitude spectrum of seismic reflection waves while maintaining high computational efficiency and achieving stable decomposition of seismic data is one of the core issues in improving the interpretability of seismic data.

[0004] Existing similar time-frequency decomposition methods directly calculate the Fourier series coefficients as functions of time by inverting the basis of the truncated Euler formula kernel within a moving time window. This method can obtain time-spectrum data with high time and frequency resolution, improving upon the limitations of traditional short-time Fourier transforms due to time window length constraints, reducing window tailing effects and interference from surrounding seismic events, and also addressing the time resolution problem of continuous wavelet transforms at low frequencies. However, this method has weak noise resistance, and its computational efficiency for large-scale seismic data requires improvement.

[0005] An existing energy-separation-based Wegener-Ville time-frequency decomposition method mentions a similar time-frequency decomposition approach: it expands the original one-dimensional signal x(t) into a Fourier-Bessel series, obtaining the zeroth-order Fourier-Bessel series expansion and the Fourier-Bessel coefficient sequence. An energy separation algorithm is used to calculate the amplitude envelope of the Fourier-Bessel coefficient sequence, obtaining M sets of Fourier-Bessel coefficient subsequences. An inverse Fourier-Bessel transform is performed on these subsequences to obtain the time-domain sub-signal xi(t) of each subsequence. The Wegener-Ville distribution of each sub-signal xi(t) is then obtained. Finally, the Wegener-Ville distributions of all sub-signals are superimposed to obtain the Wegener-Ville time-frequency spectrum of the one-dimensional signal x(t) without cross terms. This invention can suppress the generation of cross terms, making the obtained time-frequency spectrum more accurate. However, this method uses amplitude envelopes in the energy separation algorithm to calculate subsequences of the Fourier-Bessel coefficient sequence, which introduces certain errors, especially for low-energy or noisy signals, posing a risk of information loss. Furthermore, it cannot completely and accurately suppress the generation of cross terms for non-stationary signals or signals with abrupt changes.

[0006] An existing time-frequency decomposition method for multi-impact vibration signals involves performing a short-time Fourier transform (SFT) on the multi-impact vibration signals to convert them into time-frequency domain signals. A quadratic convex optimization model is then established with the objective of minimizing the square of the 2-norm of the time-frequency two-dimensional amplitude moments of the time-frequency signals. This transforms the multi-impact signal time-frequency decomposition problem into a variational problem involving the impact time-frequency dimension center and the decomposed time-frequency sub-signals. The alternating direction multiplier method is then used to further transform the solution problem into iterative solutions for the time-frequency center and the decomposed time-frequency sub-signals. Finally, an inverse short-time Fourier transform is performed on the decomposed time-frequency sub-signals to obtain their time-domain waveforms, thus completing the multi-impact vibration signal time-frequency decomposition process. While this method uses a quadratic convex optimization model to improve the accuracy of time-frequency decomposition, it increases computational complexity and workload, especially for large-scale impact signal decomposition problems.

[0007] An existing fractional-order synchronous extraction method for generalized S-transform time-frequency decomposition and reconstruction includes the following steps: Inputting x(t), performing fractional-order Fourier transform time-frequency analysis on x(t) using different rotation angles α, selecting the optimal rotation angle αopt, performing a fractional-order generalized S-transform on the signal, obtaining the fractional-order generalized S-transform value, taking its modulus, and then calculating a fractional-order instantaneous frequency estimate. Based on the instantaneous frequency estimate, a synchronous extraction operator is obtained. This operator is used to extract the fractional-order generalized S-transform time-frequency spectrum, retaining effective energy, and finally reconstructing the signal. This method can extend the time-frequency representation of a signal to a fractional-order time-frequency representation. By adaptively selecting the optimal rotation angle, it can identify the details of the time-frequency information. Furthermore, by adjusting the window function parameters and the "extraction" operation, the energy of the fractional-order synchronously extracted generalized S-transform time-frequency spectrum is more concentrated, thereby improving the time-frequency resolution of the signal. However, the operation of retaining effective energy during the extraction process may introduce certain errors, especially for low-energy signals or noisy signals, posing a risk of information loss. For complex signals or high-frequency time-frequency information, it cannot completely and accurately identify and reconstruct the signal.

[0008] An existing seismic inversion spectral decomposition method obtains a zero-phase seismic wavelet from the amplitude spectrum of well-side seismic data, then processes this wavelet to obtain a single-frequency seismic wavelet and establishes a wavelet library matrix, and finally performs time-spectral decomposition by combining it with a priori information matrix. This method is computationally complex, its performance depends heavily on the accuracy of the wavelet estimation, and it struggles to achieve the desired results with few or no wells. Furthermore, its high computational cost makes it unsuitable for large-scale data processing.

[0009] An existing ideal seismic spectral decomposition method based on variable-phase Ricker wavelet matched pursuit (RPP) constructs a variable-phase Ricker wavelet library, then uses the PRP algorithm to obtain the matched wavelet and its corresponding amplitude coefficient, and multiplies the amplitude coefficient with the wavelet atom to obtain the decomposed signal. Next, a Fourier transform is performed on the sub-signal corresponding to the single-reflection interface to calculate the amplitude and phase spectra of the seismic wave reflected from that interface, ultimately obtaining the ideal spectral decomposition result. However, the PRP algorithm used in this method is computationally intensive and requires the seismic signal morphology to meet Ricker wavelet characteristics, which differs from the wavelet morphology of real seismic data, thus lacking universal applicability. Summary of the Invention

[0010] In view of this, in order to solve the above-mentioned technical problems or some of the technical problems, the present invention provides a sparse-constrained least-squares seismic time-frequency analysis method and apparatus.

[0011] In a first aspect, embodiments of the present invention provide a sparse-constrained least-squares seismic time-frequency analysis method, comprising:

[0012] Real seismic traces are converted to complex numbers to construct complex seismic traces;

[0013] Construct a diagonal matrix with data weights and window functions as its diagonal, and initialize the model weights as an identity matrix;

[0014] Construct an Euler formula kernel matrix in the time domain consisting of different cutoff frequencies and time ranges;

[0015] The data weight constraints are applied to the complex seismic traces, and the data weight constraints and model weight constraints are applied to the Euler formula kernel matrix;

[0016] The error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces is constructed. L1 regularized sparsity constraints are applied to the model parameters. The least squares seismic time-frequency analysis results with sparse constraints are obtained by applying a greedy iterative shrinkage threshold method.

[0017] In one possible implementation, the method further includes:

[0018] Representing real seismic traces in terms of real and imaginary parts, and adding a Hilbert transform as an additional constraint, complex seismic traces are constructed using the first formula: d = d r +id i ;

[0019] Where d is the windowed segment of the complex seismic trace, d r It is a windowed segment of a real seismic trace, d i It is a window segment that performs a Hilbert transform on the seismic trace.

[0020] In one possible implementation, the method further includes:

[0021] Data weights are constructed using a second formula, which is:

[0022] Among them, W d Here, nΔt is the time relative to the center of the window, abs(d0) is the instantaneous amplitude at the center of the window function, and l is the window length.

[0023] The model weights are initialized to an identity matrix using the third formula, which is: W m =I;

[0024] Among them, W m Here, I represents the model weights, and I is the identity matrix.

[0025] In one possible implementation, the method further includes:

[0026] The Euler formula kernel matrix with different cutoff frequencies and time ranges is constructed in the time domain using the fourth formula, which is: F(t,f)=cos(2πkΔfmΔt)+isin(2πkΔfmΔt);

[0027] In this context, the Euler formula kernel matrix F represents a complex sine signal with a window length, the number of columns k in F is the number of frequencies, and the number of rows m in F is the number of samples in the time window.

[0028] In one possible implementation, the method further includes:

[0029] Data weight constraints W are applied to the complex seismic traces. d d;

[0030] Applying data weight constraints and model weight constraints to the Euler formula kernel matrix using the fifth formula, the fifth formula is: F w =W d FW m ;

[0031] Among them, F w W is the Euler's formula kernel matrix after applying data weight constraints and model weight constraints. d For data weights, W m Here, F represents the model weights, and F is the Euler formula kernel matrix.

[0032] The model parameter vector is weighted using a fifth formula, which is:

[0033] Where m is the unknown model parameter vector, i.e., the amplitude spectral coefficients that need to be obtained. w It is the weighted model parameter vector.

[0034] In one possible implementation, the method further includes:

[0035] The error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces is constructed, and L1 regularized sparsity constraints are applied to the model parameters, as shown in Equation 7:

[0036] In one possible implementation, the method further includes:

[0037] The solution x of the iteration k y k Initialize by setting step size γ0 = 1, S > 1, ξ < 1, and iteratively calculate x using the eighth and ninth formulas. k+1 The eighth formula is: y k =x k+(x k -x k-1 The ninth formula is:

[0038] If (y k -x k+1 ) T (x k+1 -x k If y ≥ 0, then y k =x k ;

[0039] If ||x k+1 -x k ||≥S||x1-x0||, then γ=max(ξ,γ1);

[0040] If ||x k+1 -x k If ||<ε, then output the iterative solution x. k+1 Otherwise, repeat the iteration steps.

[0041] Secondly, embodiments of the present invention provide a sparse-constrained least-squares seismic time-frequency analysis apparatus, comprising:

[0042] The module is used to convert real seismic traces into complex seismic traces.

[0043] The building module is used to construct a diagonal matrix with data weights and window functions as the diagonal, and to initialize the model weights as an identity matrix;

[0044] The module is used to construct the Euler formula kernel matrix in the time domain, which consists of different cutoff frequencies and time ranges.

[0045] The constraint module is used to apply the data weight constraints to the complex seismic traces and to apply data weight constraints and model weight constraints to the Euler formula kernel matrix.

[0046] The analysis module is used to construct the error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces, and to apply L1 regularized sparse constraints to the model parameters. The greedy iterative shrinkage threshold method is used to invert and obtain the least squares seismic time-frequency analysis results with sparse constraints.

[0047] Thirdly, embodiments of the present invention provide an electronic device, including: a processor and a memory, wherein the processor is configured to execute a sparse-constrained least-squares seismic time-frequency analysis program stored in the memory, so as to implement the sparse-constrained least-squares seismic time-frequency analysis method described in the first aspect above.

[0048] Fourthly, embodiments of the present invention provide a storage medium, comprising: the storage medium storing one or more programs, the one or more programs being executable by one or more processors to implement the sparse-constrained least-squares seismic time-frequency analysis method described in the first aspect above.

[0049] The sparse-constrained least-squares seismic time-frequency analysis scheme provided in this invention constructs complex seismic traces by complexifying real seismic traces; constructs a diagonal matrix with data weights and window functions as its diagonal, and initializes the model weights as an identity matrix; constructs an Euler's formula kernel matrix in the time domain composed of different cutoff frequencies and time ranges; applies the data weight constraints to the complex seismic traces, and applies both data weight constraints and model weight constraints to the Euler's formula kernel matrix; constructs an error function between the constrained Euler's formula kernel matrix, model parameters, and constrained complex seismic traces; applies L1 regularized sparse constraints to the model parameters; and uses a greedy iterative shrinking threshold method to invert and obtain the sparse-constrained least-squares seismic time-frequency analysis results. Compared with conventional methods, which are limited by the time window length, suffer from window tailing effects, or are affected by the energy interference of surrounding seismic events, resulting in lower resolution or higher computational costs, this approach offers a significant advantage. This scheme is based on the standard least squares inversion principle and adopts rigorous theoretical derivation, resulting in high accuracy. The Hilbert transform additional constraints, model weight constraints, data weight constraints, and sparsity constraints involved in the theory all have a rigorous theoretical basis. By applying the greedy fast iterative shrinking threshold method, the inversion results under short window conditions conform to the characteristics of seismic signals and have strong noise resistance. It can be used to efficiently decompose seismic signals to obtain high-precision time-frequency spectra. Attached Figure Description

[0050] Figure 1 A flowchart illustrating a sparse-constrained least-squares seismic time-frequency analysis method provided in an embodiment of the present invention;

[0051] Figure 2 This is a schematic diagram of a synthetic seismic signal and time-frequency analysis results provided in an embodiment of the present invention;

[0052] Figure 3 This invention provides a synthetic seismic signal and time-frequency analysis result as an embodiment of the invention.

[0053] Figure 4 The results are from the traditional least squares time-frequency analysis method;

[0054] Figure 5 The results are from the least-squares fast time-frequency analysis with sparse constraints.

[0055] Figure 6 This is an image of the original seismic data;

[0056] Figure 7 This is a 10Hz profile image obtained using the traditional least squares time-frequency analysis method.

[0057] Figure 8 A 10Hz profile image of the least-squares fast time-frequency analysis method with sparse constraints;

[0058] Figure 9 This is a 20Hz profile image obtained using the traditional least squares time-frequency analysis method.

[0059] Figure 10 A 20Hz profile image of the least-squares fast time-frequency analysis method with sparse constraints;

[0060] Figure 11 This is a 60Hz profile image obtained using the traditional least squares time-frequency analysis method.

[0061] Figure 12 A 60Hz profile image of the least-squares fast time-frequency analysis method with sparse constraints;

[0062] Figure 13 A schematic diagram of a sparse-constrained least-squares seismic time-frequency analysis device provided in an embodiment of the present invention;

[0063] Figure 14 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0064] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0065] To facilitate understanding of the embodiments of the present invention, further explanations and descriptions will be provided below with reference to the accompanying drawings and specific embodiments. These embodiments do not constitute a limitation on the embodiments of the present invention.

[0066] Example 1

[0067] Figure 1 This is a flowchart illustrating a sparse-constrained least-squares seismic time-frequency analysis method provided in an embodiment of the present invention, as shown below. Figure 1 As shown, the method specifically includes:

[0068] S11. Perform complexification on real seismic traces to construct complex seismic traces.

[0069] In this embodiment of the invention, firstly, actual seismic data is input, and the real seismic traces are processed into complex numbers, representing them as real and imaginary parts to construct complex seismic traces. Complex seismic traces possess excellent mathematical properties in subsequent signal processing and can be accurately and conveniently combined with subsequent processing methods.

[0070] The steps for constructing complex seismic traces with additional constraints are as follows:

[0071] The expression for constructing the complex seismic trace is shown in equation (1):

[0072] d = d r +id i (1)

[0073] Where d is the windowed segment of the complex seismic trace; d r It is a windowed segment of a real seismic trace; d i It is a window segment that performs a Hilbert transform on the seismic trace.

[0074] S12. Construct a diagonal matrix with data weights and window functions as the diagonal, and initialize the model weights as an identity matrix.

[0075] Data weights and model weights are constructed. Data weight constraints are applied by weighting the data using a diagonal matrix with the Hanning window function as its diagonal line. This emphasizes the spectral information of the data at the center of the window during the inversion process, thereby improving the accuracy and stability of the inversion results. Model weight constraints, on the other hand, are applied by weighting the inversion results using a model weight matrix initialized as an identity matrix. This makes the inversion results smoother, avoiding excessive detail or noise, and effectively preventing overfitting. Therefore, in this sparse-constrained least-squares fast time-frequency analysis method, both data weight constraints and model weight constraints aim to improve the accuracy and stability of the inversion results, leading to more accurate and reliable time-frequency analysis results.

[0076] The steps for constructing data weights and model weights are as follows:

[0077]

[0078] Among them, W d Here, is the data weight; d0 is the tracking sample at the center of the window; nΔt is the time relative to the center of the window; abs(d0) is the data size at the center of the window, which scales the basis function to the data size at the center of the window; l is the window length. The initialized model weights are:

[0079] W m =I (3)

[0080] Among them, W m is the model weight; I is the identity matrix.

[0081] S13. Construct the Euler formula kernel matrix in the time domain, consisting of different cutoff frequencies and time ranges.

[0082] The Euler kernel matrix is ​​constructed in the time domain by consisting of sinusoidal signals with different cutoff frequencies and time ranges. The construction steps are as follows:

[0083] F(t,f)=cos(2πkΔfmΔt)+isin(2πkΔfmΔt) (4)

[0084] Wherein, the kernel matrix F is a complex sinusoidal signal of window length; the number of columns k in F is the number of frequencies; and the number of rows m in F is the number of samples in the time window.

[0085] S14. Apply the data weight constraints to the complex seismic traces, and apply data weight constraints and model weight constraints to the Euler formula kernel matrix.

[0086] Applying data weight constraints W to complex seismic traces d d. Apply data weight constraints and model weight constraints to the Euler formula kernel matrix, following these steps:

[0087] F w =W d FW m (5)

[0088] Among them, F w This is the kernel matrix after applying data weight constraints and model weight constraints.

[0089]

[0090] Where m is the unknown model parameter vector, i.e., the frequency coefficients that need to be calculated; m w It is the weighted model parameter vector.

[0091] S15. Construct the error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces, and apply L1 regularized sparse constraints to the model parameters. Use the greedy iterative shrinkage threshold method to invert and obtain the least squares seismic time-frequency analysis results with sparse constraints.

[0092] Construct the constrained Euler formula kernel matrix F w Model parameters m w and constrained complex seismic trace W d The error function between d is given by applying L1 regularization sparsity constraints to the model parameters as shown in equation (7):

[0093]

[0094] By applying a greedy fast iterative threshold shrinkage method, the least squares fast time-frequency analysis results of sparse constraints are obtained through inversion. The inversion first involves converting the iterative solution x... k y k Initialize by setting step size γ0 = 1, S > 1, ξ < 1, and iteratively calculate x. k+1 :

[0095] y k =x k +(x k -x k-1 (8)

[0096]

[0097] like

[0098] (y k -x k+1 ) T (x k+1 -x k If y ≥ 0, then y k =x k (10)

[0099] like

[0100] ||x k+1 -x k ||≥S||x1-x0||, then γ=max(ξ,γ1) (11)

[0101] If ||x k+1 -x k If ||<ε, then output the iterative solution x. k+1 Otherwise, repeat iteration steps (8)-(11).

[0102] The sparse-constrained least-squares seismic time-frequency analysis method provided in this embodiment is based on the standard least-squares inversion principle and employs rigorous theoretical derivation, resulting in high accuracy. The Hilbert transform additional constraints, model weight constraints, data weight constraints, and sparse constraints involved in the theory all have a sound theoretical foundation, and the inversion results under short window conditions conform to the characteristics of seismic signals. It can be used to efficiently decompose seismic signals to obtain high-precision time-frequency spectra, and can also extract amplitude anomaly attributes to indicate thin-layered or attenuated geological bodies.

[0103] The sparse-constrained least-squares seismic time-frequency analysis method provided in this invention constructs complex seismic traces by complexifying real seismic traces; constructs a diagonal matrix with data weights and window functions as its diagonal, and initializes the model weights as an identity matrix; constructs an Euler's formula kernel matrix in the time domain composed of different cutoff frequencies and time ranges; applies the data weight constraints to the complex seismic traces, and applies both data weight constraints and model weight constraints to the Euler's formula kernel matrix; constructs an error function between the constrained Euler's formula kernel matrix, model parameters, and constrained complex seismic traces; applies L1 regularized sparse constraints to the model parameters; and uses a greedy iterative shrinking threshold method to invert and obtain the sparse-constrained least-squares seismic time-frequency analysis results. Compared with conventional methods, which are limited by the time window length, suffer from window tailing effects, or are affected by the energy interference of surrounding seismic events, resulting in lower resolution or higher computational complexity, this method offers a significant advantage. This method is based on the standard least squares inversion principle and employs rigorous theoretical derivation, resulting in high accuracy. The Hilbert transform additional constraints, model weight constraints, data weight constraints, and sparsity constraints involved in the theory all have a rigorous theoretical foundation. By applying a greedy fast iterative threshold shrinkage method, the inversion results under short window conditions conform to the characteristics of seismic signals and have strong noise resistance. It can be used to efficiently decompose seismic signals to obtain high-precision time-spectrum data.

[0104] Example 2

[0105] like Figure 1 As shown in S11, actual seismic data is input, and the real seismic traces are processed into complex numbers, representing them as real and imaginary parts to construct complex seismic traces. Complex seismic traces possess excellent mathematical properties in subsequent signal processing and can be accurately and conveniently combined with subsequent processing methods.

[0106] In S12, data weights and model weights are constructed. Data weight constraints use a diagonal matrix with the Hanning window function as its diagonal to weight the data, emphasizing the spectral information of the data at the center of the window during the inversion process, thereby improving the accuracy and stability of the inversion results. Model weight constraints use a model weight matrix initialized as the identity matrix to weight the inversion results, making them smoother and avoiding excessive detail or noise, while also effectively preventing overfitting. Therefore, in this sparse-constrained least-squares fast time-frequency analysis method, both data weight constraints and model weight constraints aim to improve the accuracy and stability of the inversion results, leading to more accurate and reliable time-frequency analysis results.

[0107] In S13, a kernel matrix is ​​constructed in the time domain, consisting of a cutoff frequency and a sinusoidal signal within the time range.

[0108] In S14, data weight constraints are applied to the complex seismic traces, and data weight constraints and model weight constraints are applied to the kernel matrix.

[0109] In S15, the error function of the error function between the constraint kernel matrix, model parameters, and the complex seismic traces of the constraints is used to apply L1 regularized sparse constraints to the model parameters. The sparse-constrained least-squares fast time-frequency analysis results are obtained through inversion. This invention's sparse-constrained least-squares fast time-frequency analysis method is based on the standard least-squares inversion principle, employs rigorous theoretical derivation, and possesses high accuracy. The Hilbert transform additional constraints, model weight constraints, data weight constraints, and sparse constraints involved in the theory all have a rigorous theoretical foundation, and the inversion results under short window conditions conform to the characteristics of seismic signals. It can be used to efficiently decompose seismic signals to obtain high-precision time-frequency spectra, and can also extract amplitude anomaly attributes to indicate thin-layer or attenuated anomalous geological bodies.

[0110] Example 3

[0111] In a specific embodiment of this method, synthetic seismic records and actual seismic data are used. The synthetic seismic records are a convolutional model of Ricker wavelet superposition, such as... Figure 2 Left image and Figure 3 As shown in the left figure, and with the addition of a certain amount of Gaussian white noise, this embodiment includes the following steps:

[0112] As an optimization, real seismic traces are complexified, representing them as real and imaginary parts to construct complex seismic traces. Complex seismic traces possess excellent mathematical properties in subsequent signal processing, allowing for accurate and convenient integration with subsequent processing methods.

[0113] Add a Hilbert transform as an additional constraint to construct a complex seismic trace. The steps for constructing a complex seismic trace with additional constraints are as follows:

[0114] The expression for constructing the complex seismic trace is shown in equation (1):

[0115] d = d r +id i (1)

[0116] Where d is the windowed segment of the complex seismic trace; d r It is a windowed segment of a real seismic trace; d i It is a window segment that performs a Hilbert transform on the seismic trace.

[0117] As an optimization, data weight constraints use a diagonal matrix with the Hanning window function as its diagonal to weight the data, emphasizing the spectral information of the data at the center of the window during the inversion process, thereby improving the accuracy and stability of the inversion results. Model weight constraints, on the other hand, use a model weight matrix initialized as an identity matrix to weight the inversion results, making them smoother and avoiding excessive detail or noise, while also effectively preventing overfitting. Therefore, in this sparse-constrained least-squares fast time-frequency analysis method, both data weight constraints and model weight constraints aim to improve the accuracy and stability of the inversion results, leading to more accurate and reliable time-frequency analysis results.

[0118] The data weights are constructed as a diagonal matrix with the Hanning window function as the diagonal, and the model weights are initialized as an identity matrix, which can be updated by the spectrum of the data at the center of the window.

[0119] The steps for constructing data weights and model weights are as follows:

[0120]

[0121] Among them, W d Here, is the data weight; d0 is the tracking sample at the center of the window; nΔt is the time relative to the center of the window; abs(d0) is the data size at the center of the window, which scales the basis function to the data size at the center of the window; l is the window length. The initialized model weights are:

[0122] W m =I (3)

[0123] Among them, W m is the model weight; I is the identity matrix.

[0124] Construct a kernel matrix in the time domain consisting of sinusoidal signals with a cutoff frequency and a specified time range. The steps for constructing the kernel matrix in the time domain are as follows:

[0125] F(t,f)=cos(2πkΔfmΔt)+isin(2πkΔfmΔt) (4)

[0126] Wherein, the kernel matrix F is a complex sinusoidal signal of window length; the number of columns k in F is the number of frequencies; and the number of rows m in F is the number of samples in the time window.

[0127] As an optimization, data weight constraints W are applied to the complex seismic traces. d d. Apply data weight constraints and model weight constraints to the kernel matrix, following these steps:

[0128] F w =W d FWm (5)

[0129] Among them, F w This is the kernel matrix after applying data weight constraints and model weight constraints.

[0130]

[0131] Where m is the unknown model parameter vector, i.e., the frequency coefficients that need to be calculated; m w It is the weighted model parameter vector.

[0132] As an optimization, an error function is constructed between the constraint kernel matrix, model parameters, and complex seismic traces of the constraints. L1 regularized sparse constraints are applied to the model parameters, and the least squares fast time-frequency analysis results of sparse constraints are obtained through rapid inversion under a short window.

[0133] Constraint kernel matrix F w Model parameters m w and constrained complex seismic traces W d The error function between d is given by applying L1 regularization sparse constraints to the model parameters as shown in equation (7):

[0134]

[0135] By applying a greedy fast iterative threshold shrinkage method, the least squares fast time-frequency analysis results of sparse constraints are obtained through inversion. The inversion first involves converting the iterative solution x... k y k Initialize by setting step size γ0 = 1, S > 1, ξ < 1, and iteratively calculate x. k+1 :

[0136] y k =x k +(x k -x k-1 (8)

[0137]

[0138] like

[0139] (y k -x k+1 ) T (x k+1 -x k If y ≥ 0, then y k =x k (10)

[0140] like

[0141] ||x k+1 -xk ||≥S||x1-x0||, then γ=max(ξ,γ1) (11)

[0142] If ||x k+1 -x k If ||<ε, then output the iterative solution x. k+1 Otherwise, repeat iteration steps (8)-(11).

[0143] Figure 2 The figures show the synthetic seismic signal and time-frequency analysis results provided in the embodiments of the present invention; the left figure shows the synthetic seismic signal of the Ricker wavelet convolution of the reflection coefficient sequence; the middle figure shows the results of the traditional least squares time-frequency analysis method; and the right figure shows the results of the least squares fast time-frequency analysis with sparse constraints. Figure 3 The figures show the results of time-frequency analysis of the synthesized seismic signal provided in this embodiment of the invention; the left figure shows the synthesized seismic signal of Ricker wavelet convolution of reflection coefficient sequence, with 10dB of Gaussian white noise added; the middle figure shows the results of the traditional least squares time-frequency analysis method; and the right figure shows the results of the sparse-constrained least squares fast time-frequency analysis. Compared with the traditional method, the traditional method takes about 10 times longer to run, resulting in a significant improvement in computational efficiency and noise resistance.

[0144] Figure 4 The results are for the traditional least squares time-frequency analysis method; the left figure shows the result of the traditional least squares time-frequency analysis method after 2 iterations; the middle figure shows the result of the traditional least squares time-frequency analysis method after 5 iterations; and the right figure shows the result of the traditional least squares time-frequency analysis method after 10 iterations. Figure 5 The results of sparse-constrained least-squares fast time-frequency analysis are shown in the left and right figures. The left figure shows the result with a weight of 0.5 for the L1 constraint term; the middle figure shows the result with a weight of 5 for the L1 constraint term; and the right figure shows the result with a weight of 15 for the L1 constraint term. Traditional methods increase resolution by controlling the number of iterations, but this significantly increases the computational cost. This method improves time-frequency resolution by changing the weight of the L1 constraint term, accurately reflecting the time position and dominant frequency information of the reflected wave, without significantly increasing the computational cost.

[0145] Figure 6 The image shows the original seismic data. The red box indicates the selected target analysis area, and the upper right corner shows a magnified image of the area within the red box. Figure 7 This is a 10Hz profile of the traditional least squares time-frequency analysis method. The red box represents the selected target analysis area, and the upper right corner shows a magnified image of the part within the red box. Figure 8 This is a 10Hz profile of the least squares fast time-frequency analysis method with sparse constraints. The red box represents the selected target analysis region, and the upper right corner shows a magnified image of the part within the red box. Figure 9This is a 20Hz profile of the traditional least squares time-frequency analysis method. The red box represents the selected target analysis area, and the upper right corner shows a magnified image of the part within the red box. Figure 10 This is a 20Hz profile of the least squares fast time-frequency analysis method with sparse constraints. The red box represents the selected target analysis region, and the upper right corner shows a magnified image of the part within the red box. Figure 11 This is the 60Hz profile of the traditional least squares time-frequency analysis method. The red box represents the selected target analysis area, and the upper right corner shows a magnified image of the part within the red box. Figure 12 This image shows the 60Hz profile of the sparse-constrained least-squares fast time-frequency analysis method. The red box indicates the selected target analysis region, and the upper right corner shows a magnified image of the area within the red box. Compared to traditional least-squares spectral decomposition methods, this method is more flexible, more sensitive to anomalies, and can highlight local structural information. Real-world data validates the correctness of the method. Spectral profiles at different frequencies demonstrate that, compared to traditional methods, this method is more sensitive to frequency response, which is beneficial for characterizing small structures and shows good application potential.

[0146] This method addresses the problems of low resolution in seismic data, window tailing effect in time-frequency analysis, and low computational efficiency. Under the constraints of Hilbert transform, model weight constraints, data weight constraints, and sparsity constraints, it solves the least squares problem. In practical applications, the time consumption is only one-tenth that of traditional constrained least squares spectral decomposition, and the high-precision time spectrum of seismic signals is obtained through efficient decomposition.

[0147] Figure 13 This is a schematic diagram of a sparse-constrained least-squares seismic time-frequency analysis device provided in an embodiment of the present invention. Figure 13 As shown, the device includes:

[0148] Module 1301 is used to perform complexification processing on real seismic traces to construct complex seismic traces. For detailed explanations, please refer to the relevant descriptions in the above method embodiments; they will not be repeated here.

[0149] Module 1301 is used to construct a diagonal matrix with data weights and window functions as its diagonals, and to initialize the model weights as an identity matrix. For detailed explanations, please refer to the relevant descriptions in the above method embodiments; they will not be repeated here.

[0150] Module 1301 is used to construct an Euler's formula kernel matrix in the time domain, consisting of different cutoff frequencies and time ranges. For detailed explanations, please refer to the relevant descriptions in the above method embodiments; they will not be repeated here.

[0151] The constraint module 1302 is used to apply the data weight constraints to the complex seismic traces and to apply both data weight constraints and model weight constraints to the Euler formula kernel matrix. For detailed explanations, please refer to the relevant descriptions in the above method embodiments; they will not be repeated here.

[0152] Analysis module 1303 is used to construct the error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces, and to apply L1 regularized sparse constraints to the model parameters. A greedy iterative shrinkage threshold method is then used to invert and obtain the least-squares seismic time-frequency analysis results with sparse constraints. For detailed explanations, please refer to the relevant descriptions in the above method embodiments; they will not be repeated here.

[0153] The sparse-constrained least-squares seismic time-frequency analysis device provided in this embodiment of the invention is used to execute the sparse-constrained least-squares seismic time-frequency analysis method provided in the above embodiment. Its implementation method and principle are the same. For details, please refer to the relevant description of the above method embodiment, which will not be repeated here.

[0154] Figure 14 An electronic device according to an embodiment of the present invention is shown, such as... Figure 14 As shown, the electronic device may include a processor 1401 and a memory 1402, wherein the processor 1401 and the memory 1402 may be connected via a bus or other means. Figure 14 Taking the example of a connection between China and Israel via a bus.

[0155] Processor 1401 may be a central processing unit (CPU). Processor 1401 may also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, or combinations thereof.

[0156] The memory 1402, as a non-transitory computer-readable storage medium, can be used to store non-transitory software programs, non-transitory computer-executable programs, and modules, such as the program instructions / modules corresponding to the methods provided in the embodiments of the present invention. The processor 1401 executes various functional applications and data processing of the processor by running the non-transitory software programs, instructions, and modules stored in the memory 1402, thereby implementing the methods in the above-described method embodiments.

[0157] The memory 1402 may include a program storage area and a data storage area. The program storage area may store the operating system and applications required for at least one function; the data storage area may store data created by the processor 1401, etc. Furthermore, the memory 1402 may include high-speed random access memory and may also include non-transitory memory, such as at least one disk storage device, flash memory device, or other non-transitory solid-state storage device. In some embodiments, the memory 1402 may optionally include memory remotely located relative to the processor 1401, and these remote memories may be connected to the processor 1401 via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0158] One or more modules are stored in memory 1402 and, when executed by processor 1401, perform the methods described in the above method embodiments.

[0159] The specific details of the aforementioned electronic device can be understood by referring to the relevant descriptions and effects in the above method embodiments, and will not be repeated here.

[0160] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by hardware related to computer program instructions. The program can be stored in a computer-readable storage medium, and when executed, it can include the processes of the embodiments of the methods described above. The storage medium can be a magnetic disk, optical disk, read-only memory (ROM), random access memory (RAM), flash memory, hard disk drive (HDD), or solid-state drive (SSD), etc.; the storage medium can also include combinations of the above types of memory.

[0161] Although embodiments of the invention have been described in conjunction with the accompanying drawings, those skilled in the art can make various modifications and variations without departing from the spirit and scope of the invention, and such modifications and variations all fall within the scope defined by the appended claims.

Claims

1. A sparse-constrained least-squares seismic time-frequency analysis method, characterized in that, include: Real seismic traces are converted to complex numbers to construct complex seismic traces; Construct a diagonal matrix with data weights and window functions as its diagonal, and initialize the model weights as an identity matrix; Construct an Euler formula kernel matrix in the time domain consisting of different cutoff frequencies and time ranges; The data weight constraints are applied to the complex seismic traces, and the data weight constraints and model weight constraints are applied to the Euler formula kernel matrix; The error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces is constructed. L1 regularized sparsity constraints are applied to the model parameters. The least squares seismic time-frequency analysis results with sparse constraints are obtained by applying a greedy iterative shrinkage threshold method.

2. The method according to claim 1, characterized in that, The process of converting real seismic traces to complex numbers to construct complex seismic traces includes: Representing real seismic traces in terms of real and imaginary parts, and adding a Hilbert transform as an additional constraint, complex seismic traces are constructed using the first formula: d = d r +id i ; Where d is the windowed segment of the complex seismic trace, d r It is a windowed segment of a real seismic trace, d i It is a window segment that performs a Hilbert transform on the seismic trace.

3. The method according to claim 2, characterized in that, The construction of data weights is a diagonal matrix with window functions as its diagonal, and the model weights are initialized to an identity matrix, including: Data weights are constructed using a second formula, which is: Among them, W d Here, nΔt is the time relative to the center of the window, abs(d0) is the instantaneous amplitude at the center of the window function, and l is the window length. The model weights are initialized to an identity matrix using the third formula, which is: W m =I; Among them, W m Here, I represents the model weights, and I is the identity matrix.

4. The method according to claim 3, characterized in that, The construction of the Euler formula kernel matrix in the time domain, consisting of different cutoff frequencies and time ranges, includes: The Euler formula kernel matrix with different cutoff frequencies and time ranges is constructed in the time domain using the fourth formula, which is: F(t,f)=cos(2πkΔfmΔt)+isin(2πkΔfmΔt); In this context, the Euler formula kernel matrix F represents a complex sine signal with a window length, the number of columns k in F is the number of frequencies, and the number of rows m in F is the number of samples in the time window.

5. The method according to claim 4, characterized in that, The application of the data weight constraints to the complex seismic traces, and the application of data weight constraints and model weight constraints to the Euler formula kernel matrix, include: Data weight constraints W are applied to the complex seismic traces. d d; Applying data weight constraints and model weight constraints to the Euler formula kernel matrix using the fifth formula, the fifth formula is: F w =W d FW m ; Among them, F w W is the Euler's formula kernel matrix after applying data weight constraints and model weight constraints. d For data weights, W m Here, F represents the model weights, and F is the Euler formula kernel matrix. The model parameter vector is weighted using a fifth formula, which is: Where m is the unknown model parameter vector, i.e., the amplitude spectral coefficients that need to be obtained. w It is the weighted model parameter vector.

6. The method according to claim 5, characterized in that, The process involves constructing the constrained Euler formula kernel matrix, the error function between the model parameters and the constrained complex seismic traces, and applying L1 regularized sparsity constraints to the model parameters, including: The error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces is constructed, and L1 regularized sparsity constraints are applied to the model parameters, as shown in Equation 7:

7. The method according to claim 6, characterized in that, The application of the greedy iterative shrinkage threshold method to invert and obtain sparsely constrained least-squares seismic time-frequency analysis results includes: The solution x of the iteration k y k Initialize by setting step size γ0 = 1, S > 1, ξ < 1, and iteratively calculate x using the eighth and ninth formulas. k+1 The eighth formula is: y k =x k +(x k -x k-1 The ninth formula is: If (y k -x k+1 ) T (x k+1 -x k If y ≥ 0, then y k =x k ; If ||x k+1 -x k || ≥ S||x1 - x0||, then γ = max(ξ, γ1); If ||x k+1 -x k If ||<ε, then output the iterative solution x. k+1 Otherwise, repeat the iteration steps.

8. A sparse-constrained least-squares seismic time-frequency analysis device, characterized in that, include: The module is used to convert real seismic traces into complex seismic traces. The building module is used to construct a diagonal matrix with data weights and window functions as the diagonal, and to initialize the model weights as an identity matrix; The module is used to construct the Euler formula kernel matrix in the time domain, which consists of different cutoff frequencies and time ranges. The constraint module is used to apply the data weight constraints to the complex seismic traces and to apply data weight constraints and model weight constraints to the Euler formula kernel matrix. The analysis module is used to construct the error function between the constrained Euler formula kernel matrix, model parameters, and constrained complex seismic traces, and to apply L1 regularized sparse constraints to the model parameters. The greedy iterative shrinkage threshold method is used to invert and obtain the least squares seismic time-frequency analysis results with sparse constraints.

9. An electronic device, characterized in that, include: A processor and a memory, the processor being configured to execute a sparse-constrained least-squares seismic time-frequency analysis program stored in the memory to implement the sparse-constrained least-squares seismic time-frequency analysis method according to any one of claims 1 to 7.

10. A storage medium, characterized in that, The storage medium stores one or more programs, which can be executed by one or more processors to implement the sparse-constrained least-squares seismic time-frequency analysis method according to any one of claims 1 to 7.

Images