阅读说明：本技术 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 (Sparse time-frequency spectrum decomposition method, device and equipment based on mixed norm and wavelet transform and storage medium ) 是由 王治国 杨阳 高静怀 刘乃豪 张兵 于 2021-08-12 设计创作，主要内容包括：本发明公开了一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质,首先提出了一种可以更好匹配地震子波的母小波,然后将求解稀疏时频谱分解方法表示为一个带有非凸稀疏约束和L2范数共同约束的反问题。通过迭代分割的算法来获得稀疏的时频谱分解方法,最后基于这个谱分解方法计算高低频之间的差异值,从而识别海洋水合物储层的位置。通过合成数据和实际数据对比,本发明提出的稀疏谱分解方法具有较高的时频分辨率,能够更加准确的识别海洋水合物储层的位置。(The invention discloses a sparse time spectrum decomposition method, a device, equipment and a storage medium based on mixed norm and wavelet transformation, which firstly provides a mother wavelet capable of better matching seismic wavelets, and then expresses a method for solving sparse time spectrum decomposition as an inverse problem with non-convex sparse constraint and L2 norm common constraint. And finally, calculating difference values between high and low frequencies based on the spectral decomposition method so as to identify the position of the marine hydrate reservoir. By comparing the synthetic data with the actual data, the sparse spectrum decomposition method provided by the invention has higher time-frequency resolution, and can more accurately identify the position of the marine hydrate reservoir.)

1. A sparse time-frequency spectrum decomposition method based on mixed norm and wavelet transformation is characterized in that a mother wavelet which is completely analyzed and can be matched with seismic wavelets is provided, then a time-frequency spectrum decomposition method with localization is provided according to the ideas of sparse representation and mixed norm, and finally the time-frequency spectrum decomposition method is used for hydrate detection and predicting the reservoir position of free gas.

2. The sparse-time spectral decomposition method based on mixed norm and wavelet transform according to claim 1, comprising the steps of:

step 1) obtaining post-stack observation data

Collecting original seismic data, preprocessing the original seismic data to obtain post-stack observation data, and recording the post-stack observation data asWhere N is the number of time sampling points, N is the [1, N ]]Indicating that it is currently the nth sample point;

step 2), obtaining a localized time-frequency spectrum decomposition method, and constructing wavelet transformation and inverse transformation of a frame operator F and post-stack observation data y (n) by adopting the post-stack observation data obtained in the step 1);

step 3), constructing a sparse time-frequency spectrum decomposition model based on the mixed norm according to the pricing operator F, the wavelet transformation and the inverse transformation obtained in the step 2);

and 4) solving the sparse time spectrum decomposition model in the step 3) by using a split Bregman iterative algorithm to obtain time-frequency spectrum coefficients with time-frequency localization.

3. The sparse-time spectral decomposition method based on mixed norm and wavelet transform of claim 2, wherein in step 2), first, a mother wavelet of wavelet transform is constructed, the mother wavelet being defined in frequency domain as:

ψ(f)＝U(f)F(f)^α(1-F(f))^β (1)

wherein U (f) is a unit step signal;is a constructed mother wavelet basis function; f is within 0, fc]Is the frequency; f. of_cIs unit orderA cutoff frequency of the jump signal; alpha and beta are adjusting parameters for adjusting the mother wavelet form;

secondly, after determining the alpha and beta adjusting parameters of the mother wavelet, constructing the wavelet transform into a representation form of a close frame, and knowing that the expressions of the wavelet transform and the inverse transform of the post-stack observation data y (n) are as follows:

x＝F^*y (2)

y＝Fx (3)

wherein x is [ x ]_j,k]J ═ 1,2,. ·, J; k1, 2,. K; x is the coefficient of wavelet transform, J represents the scale sampling length of wavelet transform, and K represents the time sampling length of wavelet transform; y is a column vector generated by y (n) and represented as the acquired seismic signals; f is a frame operator generated from the mother wavelet psi (F), F^*Its companion operator; x is the number of_j,kExpressed as the (j, k) -th coefficient of the wavelet coefficients.

4. The sparse temporal spectral decomposition method based on mixed norm and wavelet transform of claim 2, wherein in step 3), first, according to the frame theory, when the frame operator F and the seismic signal y are known, the coefficient x for solving the wavelet transform is expressed as an inverse problem solution with constraint, and the sparse model is as shown in formula (4):

wherein, the lambda is a regularization parameter,is a penalty function;

in order to obtain the time-frequency spectrum coefficient with localization, a mixed norm is introduced into the sparse model, and the mixed norm consists of a non-convex sparse constraint and an L2 norm, so that the sparse time-frequency spectrum decomposition model is shown as formula (5):

wherein y represents the post-stack observed data collectedA generated column vector; x ═ x_j,k]J ═ 1,2,. ·, J; k1, 2,. K; x is the coefficient of wavelet transform, J represents the scale sampling length of wavelet transform, and K represents the time sampling length of wavelet transform; lambda [ alpha ]_jAnd λ₂Respectively, a regularization parameter, J ═ 1, 2.., where J denotes a scale, and J denotes a scale sampling length of wavelet transform; k represents time; | x | non-conducting phosphor₂An L2 norm regularization term is expressed to prevent the time-frequency spectrum coefficient from being too sparse; phi (x)_j,k,a_j) Is formed by a known variable a_jA penalty function of the determined sparse constraint;

compared with the traditional L1 norm sparse constraint, the non-convex penalty function is non-convex, so that the problem caused by the L1 norm not being the most sparse constraint can be avoided; the non-convex penalty function definition used is shown in equation (6):

5. the sparse time-frequency spectrum decomposition method based on mixed norm and wavelet transform of claim 4, wherein in step 3), different types of time-frequency spectrum coefficients can be obtained according to different types of penalty functions.

6. The sparse-time spectral decomposition method based on mixed norm and wavelet transform of claim 4, wherein in step 3), when the variable a in formula (6) is_jSatisfy the requirement ofThe sparse temporal spectral decomposition model (5) is convex, utilizingThe theory of convex optimization enables the optimization problem to be solved.

7. The sparse-time spectral decomposition method based on mixed norm and wavelet transform of claim 2, wherein in step 4), first, a sparse regularization parameter λ is determined₁And λ₂And an initial value x₀If the intermediate variable u is introduced, the above equation (5) becomes as shown in equation (7):

wherein u is [ u ]_j,k]J ═ 1,2,. J; k is an intermediate variable introduced having dimensions similar to the wavelet coefficients x, u_j,kThe (j, k) th element representing the intermediate variable u;

then, according to the convex optimization theory, the constrained optimization problem is changed into an unconstrained optimization problem:

mu is expressed as a regularization parameter, and the above unconstrained optimization problem (8) is changed into two sub optimization problems according to the variable segmentation principle, as shown in formula (9):

in the formula, k represents the number of iterations; u ═ u_j,k]J ═ 1,2,. J; k is an intermediate variable introduced; b ═ b_j,k]J ═ 1,2,. J; k1, 2, K also represents an intermediate variable introduced; j represents the scale sampling length of wavelet transformation, and K represents the time sampling length of wavelet transformation;

and (4) aiming at the two sub-optimization problems in the formula (9), respectively solving the sub-optimization problems, and solving the final optimal solution through alternate iteration between the two sub-optimization problems.

8. The apparatus of sparse temporal spectral decomposition method based on mixed norm and wavelet transform according to any one of claims 1 to 7, comprising:

the seismic data acquisition unit is used for preprocessing the seismic data and acquiring post-stack observation data;

the time spectrum acquisition unit is used for constructing a frame operator;

the sparse time spectrum decomposition model obtaining unit is used for introducing a mixed norm into the sparse time spectrum decomposition model and avoiding the optimization problem caused by the fact that the mixed norm is not the sparsest constraint;

and the time-frequency spectrum coefficient acquisition unit is used for acquiring the localized time-frequency spectrum coefficient.

9. A computer device comprising a memory storing a computer program and a processor, wherein the processor when executing the computer program implements the steps of the method of sparse time spectral decomposition based on mixed norm and wavelet transform of any one of claims 1 to 7.

10. A computer readable storage medium storing a computer program, wherein the computer program when executed by a processor implements the steps of the method of sparse temporal spectral decomposition based on mixed norm and wavelet transform of any one of claims 1 to 7.

Technical Field

The invention belongs to the technical field of seismic exploration, and relates to a sparse time-frequency spectrum decomposition method based on sparse representation, in particular to a sparse time-frequency spectrum decomposition method, a sparse time-frequency spectrum decomposition device, sparse time-frequency spectrum decomposition equipment and a sparse time-frequency spectrum decomposition storage medium based on mixed norm and wavelet transformation.

Background

After the seismic waves pass through the oil and gas-containing reservoir, the attenuation of high-frequency components is faster, the attenuation of low-frequency components is slower, so that the local main frequency of the seismic waves in the area is reduced, and the low-frequency shadow under the amplitude abnormality is often used for indicating the position of the oil and gas reservoir. However, the anomaly is not obvious on the original seismic data, but can be obviously found by a frequency slice obtained by a time-frequency analysis method. Therefore, time-frequency analysis methods are often used to detect where these amplitudes are abnormal, thereby indicating the location of the hydrocarbon reservoir. The reservoir of the marine hydrate is thin, and free gas in the hydrate has the characteristics of low-frequency shadow and the like, so that time-frequency transformation can be used for predicting the reservoir position in the marine hydrate. Time-frequency analysis methods have been widely applied to seismic data processing and interpretation, such as short-time fourier transform, wavelet transform, S-transform and its modified generalized S-transform, and so on. However, the time-frequency analysis method is limited by the uncertainty principle of Heisenberg, so that the time-frequency resolution is low, and the positioning of the oil and gas reservoir cannot be accurately positioned.

In order to improve the resolution of the time-frequency analysis method, a sparse inversion method is proposed by many researchers, and the theory is that the time-frequency spectrum decomposition method is expressed as an inverse problem with constraints according to the principle of sparse representation. The method has the advantage that corresponding constraint conditions can be added according to different application scenes, so that a more appropriate time-frequency spectrum decomposition method is obtained. Gholoami (2013) proposes a method based on l₁-l₂A norm-constrained sparse time-frequency spectrum decomposition method introduces l on the basis of the traditional short-time Fourier₁-l₂Norm, so as to obtain a sparse time-frequency spectrum decomposition method. Based on Gholoami' S work, Sattari (2017) proposes a method based on S transformation and l₁-l₂And (3) a sparse spectral decomposition method of norm. Chen et al (2019) propose a method based on l_pA sparse temporal spectral decomposition method of norm. The above-mentioned sparse spectrum decomposition methodHowever, the time-frequency methods have the following disadvantages, and a more accurate sparse time-frequency result cannot be obtained.

The above techniques have the following disadvantages:

(1) the traditional linear time-frequency analysis method is limited by a Heisenberg uncertainty principle, so that the time-frequency resolution is low, and the accurate identification of the marine hydrate reservoir is influenced.

(2) Although the time-frequency analysis method based on sparse representation can improve the resolution of the time-frequency method, l₁Norm is only for l₀A relaxation of the norm, hence based on l₁The sparse time-frequency spectrum decomposition method of norm cannot obtain the best result; at the same time, based on l_pThe norm optimization problem is non-convex and is easy to fall into a local optimal value, so that the frequency spectrum result in sparse is influenced.

Disclosure of Invention

In order to overcome the disadvantages of the prior art, the present invention provides a sparse temporal spectral decomposition method, device, equipment and storage medium based on mixed norm and wavelet transform. The invention provides a sparse time-frequency spectrum decomposition method, a sparse time-frequency spectrum decomposition device, a sparse time-frequency spectrum decomposition equipment and a sparse time-frequency spectrum decomposition storage medium based on mixed norm and wavelet transformation, and aims to solve the problems that in the prior art, a linear time-frequency analysis method needs to meet the uncertainty principle, so that the time-frequency resolution is low and the time-frequency resolution is based on l₁The sparse time-spectrum decomposition method of norm does not achieve the best results and is based on l_pThe norm optimization problem is non-convex and is easy to fall into a local optimal value, so that the defect problem of a frequency spectrum result during sparsity is influenced.

In order to achieve the purpose, the invention adopts the following technical scheme to realize the purpose:

the invention provides a sparse time-frequency spectrum decomposition method based on mixed norm and wavelet transformation, which comprises the steps of firstly providing a mother wavelet which is completely analyzed and can be matched with seismic wavelets, then providing a time-frequency spectrum decomposition method with localization according to the ideas of sparse representation and mixed norm, and finally applying the time-frequency spectrum decomposition method to hydrate detection to predict the reservoir position of free gas.

Preferably, the method comprises the following steps:

step 1) obtaining post-stack observation data

Collecting original seismic data, preprocessing the original seismic data to obtain post-stack observation data, and recording the post-stack observation data asWhere N is the number of time sampling points, N is the [1, N ]]Indicating that it is currently the nth sample point;

step 2), obtaining a localized time-frequency spectrum decomposition method, and constructing wavelet transformation and inverse transformation of a frame operator F and post-stack observation data y (n) by adopting the post-stack observation data obtained in the step 1);

step 3), constructing a sparse time-frequency spectrum decomposition model based on the mixed norm according to the pricing operator F, the wavelet transformation and the inverse transformation obtained in the step 2);

and 4) solving the sparse time spectrum decomposition model in the step 3) by using a split Bregman iterative algorithm to obtain time-frequency spectrum coefficients with time-frequency localization.

Preferably, in step 2), first, a mother wavelet of the wavelet transform is constructed, the mother wavelet being defined in the frequency domain as:

ψ(f)＝U(f)F(f)^α(1-F(f))^β (1)

wherein U (f) is a unit step signal;is a constructed mother wavelet basis function; f is an element of 0, f_c]Is the frequency; f. of_cIs the cut-off frequency of the unit step signal; alpha and beta are adjusting parameters for adjusting the mother wavelet form;

secondly, after determining the alpha and beta adjusting parameters of the mother wavelet, constructing the wavelet transform into a representation form of a close frame, and knowing that the expressions of the wavelet transform and the inverse transform of the post-stack observation data y (n) are as follows:

x＝F^*y (2)

y＝Fx (3)

wherein x is [ x ]_j,k]J ═ 1,2,. ·, J; k1, 2,. K; x is the coefficient of wavelet transform, J represents the scale sampling length of wavelet transform, and K represents the time sampling length of wavelet transform; y is a column vector generated by y (n) and represented as the acquired seismic signals; f is a frame operator generated from the mother wavelet psi (F), F^*Its companion operator; x is the number of_j,kExpressed as the (j, k) -th coefficient of the wavelet coefficients.

Preferably, in step 3), first, according to the frame theory, when the frame operator F and the seismic signal y are known, the coefficient x for solving the wavelet transform is expressed as an inverse problem solution with constraints, and the sparse model is shown as formula (4):

wherein, the lambda is a regularization parameter,is a penalty function;

in order to obtain the time-frequency spectrum coefficient with localization, a mixed norm is introduced into the sparse model, and the mixed norm consists of a non-convex sparse constraint and an L2 norm, so that the sparse time-frequency spectrum decomposition model is shown as formula (5):

wherein y represents the post-stack observed data collectedA generated column vector; x ═ x_j,k]J ═ 1,2,. ·, J; k1, 2,. K; x is the coefficient of wavelet transform, J represents the scale sampling length of wavelet transform, and K represents the time sampling length of wavelet transform; lambda [ alpha ]_jAnd λ₂Respectively, a regularization parameter, J ═ 1, 2.., where J denotes a scale, and J denotes a scale sampling length of wavelet transform; k represents time; | x | non-conducting phosphor₂An L2 norm regularization term is expressed to prevent the time-frequency spectrum coefficient from being too sparse; phi (x)_j,k,a_j) Is formed by a known variable a_jA penalty function of the determined sparse constraint;

compared with the traditional L1 norm sparse constraint, the non-convex penalty function is non-convex, so that the problem caused by the L1 norm not being the most sparse constraint can be avoided; the non-convex penalty function definition used is shown in equation (6):

preferably, in step 3), different types of time-frequency spectral coefficients can be obtained according to different types of penalty functions.

Preferably, in step 3), when the variable a in the formula (6) is_jSatisfy the requirement ofThe sparse time spectrum decomposition model (5) is convex, and the optimization problem can be solved by using the convex optimization theory.

Preferably, in step 4), first, a sparse regularization parameter λ is determined₁And λ₂And an initial value x₀If the intermediate variable u is introduced, the above equation (5) becomes as shown in equation (7):

wherein u is [ u ]_j,k]J ═ 1,2,. J; k is an intermediate variable introduced having dimensions similar to the wavelet coefficients x, u_j,kThe (j, k) th element representing the intermediate variable u;

then, according to the convex optimization theory, the constrained optimization problem is changed into an unconstrained optimization problem:

mu is expressed as a regularization parameter, and the above unconstrained optimization problem (8) is changed into two sub optimization problems according to the variable segmentation principle, as shown in formula (9):

in the formula, k represents the number of iterations; u ═ u_j,k]J ═ 1,2,. J; k is an intermediate variable introduced; b ═ b_j,k]J ═ 1,2,. J; k1, 2, K also represents an intermediate variable introduced; j represents the scale sampling length of wavelet transformation, and K represents the time sampling length of wavelet transformation;

and (4) aiming at the two sub-optimization problems in the formula (9), respectively solving the sub-optimization problems, and solving the final optimal solution through alternate iteration between the two sub-optimization problems.

The invention also provides a device of the sparse time spectrum decomposition method based on the mixed norm and wavelet transform, which comprises the following steps:

the seismic data acquisition unit is used for preprocessing the seismic data and acquiring post-stack observation data;

the time spectrum acquisition unit is used for constructing a frame operator;

the sparse time spectrum decomposition model obtaining unit is used for introducing a mixed norm into the sparse time spectrum decomposition model and avoiding the optimization problem caused by the fact that the mixed norm is not the sparsest constraint;

and the time-frequency spectrum coefficient acquisition unit is used for acquiring the localized time-frequency spectrum coefficient.

The invention proposes a computer device comprising a memory storing a computer program and a processor implementing the steps of a sparse time-spectrum decomposition method based on mixed norms and wavelet transforms when executing the computer program.

The present invention proposes a computer readable storage medium, storing a computer program which, when executed by a processor, implements the steps of a sparse temporal spectral decomposition method based on mixed norms and wavelet transforms.

The invention has the following beneficial effects:

the invention discloses a sparse time-frequency spectrum decomposition method, a device, equipment and a storage medium based on mixed norm and wavelet transformation, belonging to the technical field of seismic exploration. The invention firstly provides a mother wavelet which is completely analyzed and can be better matched with seismic wavelets, then provides a time-frequency spectrum decomposition method with high time-frequency localization according to the ideas of sparse representation and mixed norm, and finally uses the time-frequency spectrum decomposition method for hydrate detection to predict the reservoir position of free gas.

Furthermore, the shape of the mother wavelet can be adjusted by adjusting parameters alpha and beta of the mother wavelet in the frequency domain, so that the mother wavelet is more matched with the seismic wavelet, and high-time-frequency localized time-frequency spectrum decomposition is obtained.

Furthermore, the penalty functions in the sparse time spectrum decomposition model are of various types, and the different types of the penalty functions can obtain different types of time spectrum coefficients; in order to obtain a time-frequency spectrum coefficient with higher time-frequency localization, a mixed norm is introduced into a sparse time-frequency spectrum decomposition model, the L1 norm can realize sparsity, and the L1 is widely applied due to the optimized solving characteristic; from the perspective of a learning theory, the L2 norm can prevent overfitting and improve the generalization capability of the model; the penalty function can convert a constrained nonlinear problem into an unconstrained nonlinear programming, the unconstrained linear programming can be solved by a gradient method and the like, and the computer algorithm can be made more conveniently by using the penalty function.

Furthermore, the sparse time-frequency spectrum decomposition method represents the coefficient for solving the wavelet transformation as an inverse problem solution on the basis of the wavelet transformation, and introduces a mixed norm penalty function into the inverse problem model, thereby obtaining a time-frequency analysis method with higher time-frequency localization.

Drawings

FIG. 1 illustrates a comparison of different time-frequency methods for noise-free synthesis of seismic signals (a) synthetic seismic records, (b) wavelet transforms, (c) short-time Fourier transforms, (d) squeeze wavelet transforms, (e) time-frequency transforms based on L1 constraints, (f) time-frequency transforms proposed by the present invention);

FIG. 2 illustrates a comparison of different time-frequency transforms for a noise-containing synthetic seismic signal (a) noise-containing synthetic seismic records, (b) wavelet transforms, (c) time-frequency transforms based on L1 constraints, (d) time-frequency transforms proposed by the present invention);

in the three-dimensional seismic data of FIG. 3, the number of tracks in the Xline direction and the Inline direction is 1306 and 95 respectively, and the time sampling interval is 2 ms;

FIG. 4 Inline11 section with the black line being the BSR and free gas below the BSR;

fig. 5 shows that the black region is a region with relatively severe attenuation along the BSR attenuation profile calculated by the time-frequency spectrum decomposition method provided by the present invention.

Detailed Description

In order to make the technical solutions of the present invention better understood, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

It should be noted that the terms "first," "second," and the like in the description and claims of the present invention and in the drawings described above are used for distinguishing between similar elements and not necessarily for describing a particular sequential or chronological order. It is to be understood that the data so used is interchangeable under appropriate circumstances such that the embodiments of the invention described herein are capable of operation in sequences other than those illustrated or described herein. Furthermore, the terms "comprises," "comprising," and "having," and any variations thereof, are intended to cover a non-exclusive inclusion, such that a process, method, system, article, or apparatus that comprises a list of steps or elements is not necessarily limited to those steps or elements expressly listed, but may include other steps or elements not expressly listed or inherent to such process, method, article, or apparatus.

The invention is described in further detail below with reference to the accompanying drawings:

the invention provides a sparse time-frequency spectrum decomposition method based on mixed norm and wavelet transformation, which comprises the steps of firstly providing a mother wavelet which is completely analyzed and can be matched with seismic wavelets, then providing a time-frequency spectrum decomposition method with localization according to the ideas of sparse representation and mixed norm, and finally applying the time-frequency spectrum decomposition method to hydrate detection to predict the reservoir position of free gas.

The invention provides a sparse time-frequency spectrum decomposition method based on data driving, which comprises the following steps:

step 1) obtaining post-stack observation data

Collecting original seismic data, preprocessing the original seismic data to obtain post-stack observation data, and recording the post-stack observation data asWhere N is the number of time sampling points, N is the [1, N ]]Indicating that it is currently the nth sample point.

Step 2), obtaining a localized time-frequency spectrum decomposition method, and adopting the post-stack observation data obtained in the step 1) to construct wavelet transformation and inverse transformation of a frame operator F and post-stack observation data y (n):

first, a mother wavelet of the wavelet transform is constructed, and the mother wavelet is defined in a frequency domain as shown in formula (1):

ψ(f)＝U(f)F(f)^α(1-F(f))^β (1)

wherein U (f) is a unit step signal;is a constructed mother wavelet basis function; f is an element of 0, f_c]Is the frequency; f. of_cIs the cut-off frequency of the unit step signal; alpha and beta are adjusting parameters for adjusting the mother wavelet form, so that the mother wavelet is more matched with the seismic wavelet, and high-time frequency localized time-frequency spectrum decomposition is obtained.

Secondly, after determining the alpha and beta adjusting parameters of the mother wavelet, constructing the wavelet transform into a representation form of a close frame, and knowing that the expressions of the wavelet transform and the inverse transform of the post-stack observation data y (n) are shown as formula (2) and formula (3):

x＝F^*y, (2)

y＝Fx, (3)

wherein x is [ x ]_j,k]J ═ 1,2,. ·, J; k is 1, 2., K is a coefficient of the wavelet transform, J denotes a scale sampling length of the wavelet transform, and K denotes a wavelet transform time sampling length; y is the column vector generated by y (n); f is a frame operator generated from the mother wavelet psi (F), F^*Is the adjoint operator of F, FF is the unit matrix, and in the case of known F, F can be calculated^*，x_j,kExpressed as the (j, k) -th coefficient of the wavelet coefficients.

Step 3), constructing a sparse time-frequency spectrum decomposition model based on the mixed norm according to the pricing operator F, the wavelet transformation and the inverse transformation obtained in the step 2):

firstly, according to the frame theory, when a frame operator F and a seismic signal y are known, a coefficient x for solving the wavelet transformation is expressed as an inverse problem solution with constraint, and a sparse model is shown as a formula (4):

wherein λ is a regularization parameter.Is a penalty function. The types of the penalty functions are different, and different types of time-frequency spectrum coefficients can be obtained.

In order to obtain the time-frequency spectrum coefficient with higher time-frequency localization, a mixed norm is introduced into the sparse model, and the mixed norm consists of a non-convex sparse constraint and an L2 norm, so that the new sparse time-frequency spectrum decomposition model is shown as a formula (5):

wherein y represents the post-stack observed data collectedA generated column vector; x ═ x_j,k]J ═ 1,2,. ·, J; k1, 2,. K; x is the coefficient of wavelet transform, J represents the scale sampling length of wavelet transform, and K represents the time sampling length of wavelet transform; lambda [ alpha ]_jAnd λ₂Respectively representing regularization parameters, j represents a scale, and k represents time; | x | non-conducting phosphor₂An L2 norm regularization term is represented to prevent the time-frequency spectral coefficients from being too sparse. Phi (x)_j,k,a_j) Is formed by a known variable a_jA penalty function of the determined sparsity constraint. Compared with the traditional L1 norm sparse constraint, the penalty function is a non-convex penalty function, can avoid the problem caused by the L1 norm not being the most sparse constraint, and can also be used in the variable a_jAnd when a certain condition is met, ensuring that the optimization problem is convex optimization, and solving the optimization problem by utilizing the convex optimization. There are many kinds of non-convex penalty functions, and this patent uses a non-convex penalty function similar to the arctan form, which is defined as shown in equation (6):

when the variable a in equation (6)_jSatisfy the requirement ofThe sparse time spectrum decomposition model (5) is convex, and the sparse time spectrum decomposition model problem can be easily solved by utilizing a convex optimization theory.

And 4) solving the sparse time spectrum decomposition model in the step 3) by using a split Bregman iterative algorithm to obtain time-frequency spectrum coefficients with time-frequency localization.

First, a sparse regularization parameter λ is determined₁And λ₂And an initial value x₀. Introducing the intermediate variable u, the above equation (5) becomes the result shown in equation (7):

wherein u is [ u ]_j,k]J ═ 1,2,. J; k is an intermediate variable introduced having dimensions similar to the wavelet coefficients x, u_j,kThe (j, k) th element representing the intermediate variable u.

Then, according to the convex optimization theory, the constrained optimization problem is changed into an unconstrained optimization problem:

according to the variable segmentation principle, the above unconstrained optimization problem (8) can be changed into two sub-optimization problems as shown in formula (9):

in the formula, k represents the number of iterations; u ═ u_j,k]J ═ 1,2,. J; k is an intermediate variable introduced; b ═ b_j,k]J ═ 1,2,. J; k1, 2, K also represents an intermediate variable introduced; j represents the scale sampling length of wavelet transformation, and K represents the time sampling length of wavelet transformation;

aiming at the two sub-optimization problems in the formula (9), the sub-optimization problems are solved respectively, and the final optimal solution is solved through alternate iteration between the two optimization problems.

There are generally two ways to stop the iteration: the first is to reach the maximum number of iterations; and secondly, when the error between the current iteration result and the last iteration result is smaller than a certain threshold value, the maximum iteration times are reached.

Numerical simulation results-synthetic seismic record data:

first, the effectiveness of the present invention is verified using a noise-free synthetic seismic signal, as shown in FIG. 1 (a). The time sampling interval is 1ms, and the time sampling times are 512. The first wavelet is generated by adding a Ricker wavelet with 60Hz dominant frequency at 0.05 s. The second wavelet is generated by convolving a negative reflection coefficient with a Ricker wavelet with a dominant frequency of 50hz at 0.15 s. The third wavelet consists of a positive 40hz Ricker wavelet at 0..25s and a negative 40hz Ricker wavelet at 0.275 s. The last wavelet consists of 3 Ricker wavelets with the same dominant frequency, mainly comprising two 30Hz positive Ricker wavelets (0.35 and 0.41s) and one 30Hz negative Ricker wavelet (0.38 s). The experiment respectively compares wavelet transformation, short-time Fourier transformation, extrusion wavelet transformation, a time-frequency spectrum decomposition method based on L1 constraint and the time-frequency spectrum decomposition method provided by the invention. The maximum iteration times of the time-frequency spectrum decomposition method based on the L1 constraint and the time-frequency spectrum decomposition method provided by the invention are set to be 50. Fig. 1(b) -fig. (f) show the results of the time-frequency transformation methods.

As shown in fig. 1(b) and 1(c), the wavelet transform has better time-frequency localization compared to the short-time fourier transform. However, both of these transformations are limited by the uncertainty principle and time-frequency localization is limited. The squeeze wavelet transform has better frequency resolution as shown in fig. 1(d), but does not consider time resolution. Compared with the time-frequency spectrum decomposition result based on the L1 constraint in FIG. 1(e), the time-frequency spectrum decomposition method provided by the invention has a sparser time-frequency spectrum decomposition result as shown in FIG. 1 (f).

Next, fig. 2 shows the result of the spectral decomposition in the case of a noisy synthetic earthquake. FIG. 2(a) is generated by adding white Gaussian noise to FIG. 1(a), and the SNR is 10 dB. Since the traditional time-frequency spectrum decomposition method has poor noise resistance, the experiment only verifies the wavelet transformation, the sparse time-frequency spectrum decomposition method based on the L1 constraint and the time-frequency spectrum decomposition method provided by the invention, and the results are shown in figures 2(b) -2 (d). Obviously, the noise immunity of wavelet transformation is poor, and two inversion-based methods have better noise immunity. Compared with a sparse time-frequency spectrum decomposition method based on L1 constraint, the time-frequency spectrum decomposition method provided by the invention has better noise resistance and sparsity.

Actual seismic data profile: three-dimensional hydrate seismic data are utilized to further verify the effectiveness of the invention. FIG. 3 shows three-dimensional seismic data with Inline numbered 95 and Xline numbered 1306. The sampling time interval for each trace is 2 ms. As shown in fig. 4, in the seismic profile Inline11, the BSR, i.e., red line, can be easily traced due to high data quality. Because the free gas has the attenuation characteristic, the position of the free gas reservoir can be detected by utilizing the time-frequency spectrum decomposition method provided by the invention. The result of the attenuation profile along the BSR calculated by the time-frequency spectrum decomposition method provided by the invention is shown in FIG. 5, and the black part is the place with serious attenuation, namely the position of the free gas reservoir. Obviously, the time-frequency spectrum decomposition method provided by the invention can be used for effectively detecting the position of free gas, and an effective detection method is provided for hydrate reservoir detection.

The above-mentioned contents are only for illustrating the technical idea of the present invention, and the protection scope of the present invention is not limited thereby, and any modification made on the basis of the technical idea of the present invention falls within the protection scope of the claims of the present invention.