A multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding

Abstract

The application provides a multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding, global frequency division processing is performed on seismic data according to convolution sparse coding; secondly, two-dimensional K-SVD-based seismic data low-frequency component deconvolution processing is performed to learn the horizontal characteristics from two-dimensional low-frequency seismic data; then, convolution sparse coding-based seismic data high-frequency component deconvolution processing is studied; finally, the advantages of single-channel deconvolution and multi-channel deconvolution are complementary, vertical and horizontal constraints are added in the deconvolution objective function, so that the deconvolution result has both resolution improvement and continuity enhancement. The application learns the reservoir reflectivity characteristics and the mapping relationship with the seismic data characteristics from the logging data, has the characteristics of high vertical resolution, simultaneously adds the horizontal distribution law of the seismic data to the deconvolution processing, so that the inversion result presents strong horizontal continuity, and the resolution is significantly enhanced compared with the seismic data.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of intelligent analysis of seismic data, in particular to a multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding. BACKGROUND

[0002] Seismic data high-resolution processing is of great significance to reservoir prediction and geological structure interpretation. Reflection coefficient inversion or deconvolution processing is an important technical means to reflect the position of the stratum structure and restore the parameters of the stratum. Its research has important practical significance. Due to the band-limited characteristics of seismic data and seismic wavelet, the inversion problem faces multiple solutions and uncertainty. Therefore, various constraint terms are added to the deconvolution objective function to narrow the solution range and improve the resolution and stability of the inversion result.

[0003] Currently, the methods for improving the resolution of seismic data mainly include traditional deconvolution methods based on various assumption conditions and model-driven assumption-based methods. Such methods have poor universality, and do not mine features from the data itself, but construct constraint terms based on various assumptions. In addition, the current mainstream deconvolution methods are difficult to balance the vertical resolution and lateral continuity in the process of single-channel and multi-channel deconvolution, often appearing the phenomenon of "trade-off", resulting in "false high frequency", "pulled noodles", and deviation from the actual logging in the deconvolution result. SUMMARY

[0004] In view of the problems in the prior art, the present application provides a multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding, which comprises the following steps:

[0005] S1, performing global frequency division processing on the seismic data, and decomposing the seismic profile into a low-frequency component and a high-frequency component: S=S LFC +S HFC After the global frequency division processing of the seismic data, the corresponding components still satisfy the convolution relationship;

[0006] S2, performing deconvolution processing on the low-frequency component of the seismic data based on two-dimensional K-SVD, learning the lateral characteristics from the two-dimensional low-frequency seismic data, designing a low-pass filter, so that the reflection coefficient obtained can be sparse reconstructed using the K-SVD dictionary after passing through the low-pass filter, and adding L2 norm constraint to the reflection coefficient;

[0007] S3, performing deconvolution processing on the high-frequency component of the seismic data based on convolution sparse coding.

[0008] Preferably, the deconvolution processing on the low-frequency component of the seismic data based on two-dimensional K-SVD comprises the following steps:

[0009] S21, using K-SVD dictionary learning of two-dimensional image space features, by setting a sliding window, along the horizontal and vertical, according to the given step sampling, constructing a two-dimensional dictionary learning sample;

[0010] S22, each two-dimensional sample block y i (i = 1, 2, …, n) is converted to a one-dimensional signal by column merging

[0011] S23, the column vector of each block is normalized from 0 to 1 to obtain And all the are merged to obtain a two-dimensional sample signal Y;

[0012] S24, for the low-frequency component S LFC of the seismic data, construct a two-dimensional sample, learn a two-dimensional K-SVD dictionary D, and apply it to the multi-channel deconvolution regularization constraint term; In the deconvolution process, a plurality of small blocks are processed, and finally the deconvolution of the entire two-dimensional profile is realized by merging the small blocks; The final deconvolution objective function is:

[0013]

[0014] Wherein, and are the low-frequency component small blocks of the seismic data and the corresponding low-frequency component small blocks of the reflection coefficient, A is the sparse coefficient in the sparse representation process, μ vL and μ hL are the regularization parameters corresponding to the vertical constraint term and the horizontal constraint term in the low-frequency component deconvolution objective function, λ is the regularization parameter of the sparse constraint term;

[0015] S25, using ADMM algorithm to solve the deconvolution objective function.

[0016] Preferably, the convolution sparse coding-based deconvolution processing of the high-frequency component of the seismic data comprises the following steps:

[0017] S31, after the frequency processing of the seismic profile, the convolution dictionary learning is performed on the high-frequency component S HFC

[0018]

[0019] Wherein, d i is an n×n two-dimensional convolution dictionary, Z i is a feature mapping diagram, and the original sample S HFC maintains the same dimension, N is the number of dictionaries, and ADMM algorithm is used for solving;

[0020] ​S32, in order to construct the two-dimensional high-resolution deconvolution system of the profile, the desired output reflection coefficient is associated with the high-resolution seismic data through a low-pass filter operator, so that the desired output reflection coefficient and the high-resolution seismic data satisfy the following relationship:

[0021] S H = LR HFC

[0022] Where S H is the high-frequency component of the desired output high-resolution seismic data, L is a low-pass filter operator, and R HFC is the desired output two-dimensional reflection coefficient profile;

[0023] S33, after filtering the desired output, the convolution dictionary is added to the deconvolution objective function in the form of regularization, and the lateral regularization constraint of the two-dimensional profile is:

[0024]

[0025] The constraint is added to the deconvolution objective function in the form of regularization:

[0026]

[0027] By giving a fixed feature mapping matrix The intermediate process variable R HFC is obtained by solving the following objective function:

[0028]

[0029] S34, according to the R HFC obtained in the last step, the feature mapping matrix {Z i} is calculated under the condition of given dictionary d i , that is, the following objective function is optimized:

[0030]

[0031] Where g x and g y are filter operators for calculating gradients along the x-axis and y-axis directions respectively, and τ is a hyperparameter;

[0032] S35, the solution of the objective function is converted into solving by using ADMM algorithm.

[0033] Preferably, the solution by using ADMM algorithm includes the following steps:

[0034] Input: high-frequency component S LFC of seismic data, multi-channel wavelet matrix G, regularization parameters μ vL and μ hL, low-pass filter operator L and two-dimensional K-SVD dictionary D, maximum iteration number K, and maximum allowed reconstruction error ε; output: low-frequency component of reflectivity R LFC ; initialization: and A 0 , Lagrange multipliers U1, U2 and parameter ρ in ADMM algorithm, current iteration number k = 0;

[0035] When and k < K,

[0036] B k = G T G + μ hL I

[0037] x k = G T S LFC + μ hL L T (DA k -U2-U1 / μ vL )

[0038] Solve

[0039]

[0040]

[0041]

[0042] Determine whether the termination condition is met; if met, output R LFC and A; if not met, k = k + 1, continue the loop iteration.

[0043] Preferably, solving using ADMM algorithm includes the following steps:

[0044] Input: high-frequency component of seismic data S HFC , wavelet matrix W, high-frequency component regularization parameter μ vH and μ hH , feature mapping regularization parameter λ, Lagrange parameter ρ in ADMM algorithm, low-pass filter operator L and convolution dictionary d i (i = 1, 2, …, N), gradient regularization operator g x and g y , maximum iteration number K, and maximum allowed reconstruction error ε;

[0045] Output: high-frequency component of reflectivity R HFC , feature mapping map Z i (i = 1, 2, …, N);

[0046] Initialization: Z 0 , ADMM algorithm in Lagrange multiplier B 0 and C 0 , the current iteration number k=0;

[0047] When And k

[0048] According to the formula Update

[0049] J=1,2,…,K

[0050]

[0051]

[0052] C j+1 =C j +Z j+1 -B j+1

[0053] Let

[0054] K=k+1.

[0055] The above technical features can be combined in various suitable ways or replaced by equivalent technical features, as long as the purpose of the application can be achieved.

[0056] The multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding provided by the application has at least the following beneficial effects compared with the prior art:

[0057] The application performs global frequency division processing on seismic data according to convolution sparse coding, studies two-dimensional K-SVD-based seismic data low-frequency component deconvolution processing to learn the horizontal characteristics from two-dimensional low-frequency seismic data, then studies convolution sparse coding-based seismic data high-frequency component deconvolution processing, and finally complements the advantages of single-channel deconvolution and multi-channel deconvolution by adding vertical and horizontal constraints to the deconvolution objective function, so that the deconvolution result has both improved resolution and enhanced continuity, achieving the purpose of high-resolution inversion of vertical axis “separated” and horizontal axis “connected”.

[0058] (1) The method learns the reservoir reflectivity characteristics and the mapping relationship with the seismic data characteristics from the logging data, has the characteristics of high vertical resolution, and adds the horizontal distribution law of the seismic data to the deconvolution processing, so that the inversion result presents strong horizontal continuity and the resolution is significantly enhanced compared with the seismic data.

[0059] (2) The deconvolution result of the application is superior to the model-driven deconvolution method (for example, the sparse pulse deconvolution method).

[0060] (3) The application mines the vertical resolution characteristics from the data itself, and adds the lateral distribution law of the seismic data to the deconvolution processing, so that the inversion result presents strong lateral continuity, and the resolution is significantly enhanced compared with the seismic data, and can be well matched with the actual logging data. BRIEF DESCRIPTION OF DRAWINGS

[0061] In the following, the application will be described in more detail based on the embodiments and with reference to the accompanying drawings. Among them:

[0062] Figure 1 A two-dimensional synthetic seismic record frequency division processing profile low-frequency component result display diagram is shown;

[0063] Figure 2 A two-dimensional synthetic seismic record frequency division processing profile high-frequency component result display diagram is shown;

[0064] Figure 3 A sparse coding principle diagram is shown;

[0065] Figure 4 A signal sparse coding principle diagram is shown;

[0066] Figure 5 A two-dimensional seismic wavelet and reflection coefficient low-frequency component synthetic seismic record schematic diagram is shown;

[0067] Figure 6 A two-dimensional seismic wavelet and reflection coefficient high-frequency component synthetic seismic record schematic diagram is shown;

[0068] Figure 7 A two-dimensional K-SVD dictionary learning sampling schematic diagram is shown;

[0069] Figure 8 An actual seismic data logging seismic wavelet diagram is shown;

[0070] Figure 9 An actual seismic data logging well position distribution diagram is shown; DETAILED DESCRIPTION

[0071] The application will be further described below in combination with the accompanying drawings.

[0072] The application provides a multi-channel seismic deconvolution method based on two-dimensional K-SVD and convolution sparse coding, characterized by comprising the following steps:

[0073] S1, the seismic data is subjected to global frequency division processing, and the seismic profile is decomposed into a low-frequency component and a high-frequency component: S=S LFC +SHFC After the global spectral decomposition, the components still satisfy the convolution relationship;

[0074] S2, based on two-dimensional K-SVD, the low-frequency component of seismic data is deconvolved, and the horizontal characteristics are learned from the two-dimensional low-frequency seismic data. A low-pass filter is designed so that the reflection coefficient obtained can be sparsely reconstructed using the K-SVD dictionary after passing through the low-pass filter, and the L2 norm constraint is added to the reflection coefficient;

[0075] S3, based on convolution sparse coding, the high-frequency component of seismic data is deconvolved.

[0076] In one embodiment, for a seismic profile S, the seismic profile S can be decomposed into a low-frequency component (Low-Frequency Component, LFC) and a high-frequency component (High-Frequency Component, HFC) by solving the following optimization problem

[0077]

[0078] Wherein, Z L is the low-frequency feature map of profile S, γ is the regularization parameter, the larger γ is, the higher the proportion of low-frequency component is, and the high-frequency component is closer to the spectrum of the original seismic data, f L is a low-pass filter with all coefficients of 3*3 being 1 / 9, f dh and f dv represent horizontal and vertical gradient operators [1,-1] and [1,-1] T , * represents a convolution operator. By converting equation (0-1) to the Fourier domain for solving:

[0079]

[0080] Wherein and represent the Fourier and inverse Fourier operators respectively, and represent the Fourier transforms of f L , f dh and f dv , "(·) H " and represent the complex conjugate operator and the corresponding element point multiplication (i.e., corresponding element multiplication) symbol respectively. After obtaining Z L , the original seismic profile S can be decomposed into a low-frequency component S LFC and a high-frequency component S HFC :

[0081] S=S LFC +SHFC (0-3)

[0082] Among them, S LFC =f L *Z L Indicates its low-frequency component, S HFC These are high-frequency components, representing the high-frequency edges and texture structure of seismic data.

[0083] The Marmousi model was used to verify the effect of frequency division processing. A seismic record was synthesized by combining a 30Hz zero-phase Ricker wavelet with the reflection coefficient profile. The parameter γ in the equation was set to 10, and its low-frequency and high-frequency components are shown in the figure. Figure 1 and Figure 2 middle.

[0084] In one embodiment, the signal is described as a linear combination of several basis vectors derived from a redundant matrix, referred to as a dictionary. Simultaneously, it is required to represent the original signal using as few atoms as possible, a technique known as sparse coding.

[0085] K-SVD dictionary learning is essentially a sparse coding process. K-SVD-based dictionary learning methods have advantages such as high computational efficiency and fast convergence speed, and are widely used in signal and image processing. This algorithm aims to find a "supercomplete" dictionary matrix D for the sample set. Any sample in the set can have its corresponding sparse representation obtained from the dictionary. Dictionary learning can be viewed as a form of feature extraction in the form of matrix factorization. For the sample set Y = {y1, y2, ..., y...} M}, where M is the number of samples, and a single sample It is an N-dimensional vector. The dictionary matrix D = {d1, d2, ..., dn} K},in This is called an atomic vector, which is also an N-dimensional vector, where K represents the number of atoms. Based on the relationship between N and K, dictionaries can be divided into the following three types:

[0086]

[0087] Typically, overcomplete dictionaries are learned, meaning the number of atoms is much greater than the dimension of the atoms (K≥N). The goal of dictionary learning is to learn a dictionary D from a sample dataset Y such that it is sparsely represented as Y = DA. This sparsity can be achieved using a sparse coding matrix A = {α1, α2, ..., α...}. M}, α i (i = 1, 2, ..., K) are samples y i The sparse coefficient vector in dictionary D. The schematic diagram of sparse coding is shown below. Figure 3 As shown, the white area represents 0.

[0088] It can be seen that in the sparse coding process, the dimension of the original signal and the dimension of the dictionary atom are consistent, and the intuitive diagram is as shown in the figure: Figure 4 For this process, the objective function of dictionary learning has two forms based on the sparsity formula (0-5) and the error threshold formula (0-6):

[0089]

[0090] Where ||·||0 represents the L0 norm of the vector, that is, the number of non-zero elements in the vector, T0 is the sparsity constraint threshold, and is a constant.

[0091]

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

[0093] The K-SVD algorithm uses a step-by-step optimization strategy to solve D and A. This method takes the size K of the dictionary as the loop variable, and can also use the given decomposition residual ε as the iteration termination condition, and sequentially updates each column d k (k = 1, 2, …, K) in the dictionary D.

[0094] In an embodiment, since the low-frequency component contains the background of the original seismic data, it plays a very important role in deconvolution. The low-frequency component and the high-frequency component of the seismic data are obtained through frequency division processing, and separate processing is used to realize the deconvolution processing of the original seismic data. That is, it is currently assumed that the low-frequency component of the seismic data is obtained by deconvolution processing to obtain the low-frequency component of the reflection coefficient, and the high-frequency component of the seismic data is obtained by deconvolution processing to obtain the high-frequency component of the reflection coefficient.

[0095] For the low-frequency component S LFC of the seismic data, the spatial distribution characteristics of the seismic data are learned by using two-dimensional K-SVD dictionary, and are added to the objective function in a regularized manner to realize multi-channel deconvolution processing of the seismic data.

[0096] The multi-channel convolution model, that is, the multi-channel seismic data and the reflection coefficient are unfolded one by one, can be obtained:

[0097]

[0098] Where s i and r i represent the i-th seismic record and the corresponding reflection coefficient, respectively, and let:

[0099]

[0100] Then the objective function of multi-channel deconvolution is:

[0101]

[0102] To remove the multiple solutions of R in equation (0-9), a regularization constraint term is added in it:

[0103]

[0104] wherein, represents a constraint term added to the multi-channel reflection coefficient R.

[0105] In one embodiment, the entire signal or image is represented as the sum of the convolution of the feature map and the corresponding filter. For the observed signal The mathematical model of the convolution sparse coding is:

[0106]

[0107] wherein * is a convolution operator, is a dictionary filter, and M << L, is a sparse vector convolved with the corresponding filter, referred to as a feature map vector, which has the same size as the original signal y.

[0108] The problems to be solved by the convolution sparse coding mainly include two categories:

[0109] 1. The problem of finding the sparsest solution of the observed signal under a given dictionary, which can be solved by the convolution sparse coding;

[0110] 2. The problem of learning a dictionary to sparsely represent a given sample data set, i.e. convolution dictionary learning, which can be solved by the ADMM algorithm.

[0111] In one embodiment, the deconvolution method based on the K-SVD dictionary learning performs a small block operation on the data and realizes the deconvolution channel by channel, which destroys the consistency of the data. Meanwhile, the deconvolution method based on a single channel does not consider the continuity between adjacent channels, resulting in the "pulled noodles" phenomenon in the deconvolution result. The present application takes improving the lateral continuity of the deconvolution result as the target, and realizes the high-resolution deconvolution processing of the seismic data by taking the two-dimensional sparse coding and the frequency division processing as the basic means.

[0112] Specifically, first, according to the idea of convolution sparse coding, the seismic data is subjected to global frequency division processing, which facilitates subsequent use of K-SVD and convolution dictionary to extract horizontal features. Second, based on two-dimensional K-SVD, the seismic data low-frequency component deconvolution processing is performed to learn horizontal features from two-dimensional low-frequency seismic data. By designing a low-pass filter, the reflection coefficient obtained after passing through the low-pass filter can be sparsely reconstructed using the K-SVD dictionary. At the same time, since it is a low-frequency component, an L2 norm constraint is directly added to the reflection coefficient, so that the deconvolution result is more in line with the physical law. Then, based on the convolution sparse coding of the seismic data high-frequency component deconvolution processing, two CSC-based deconvolution methods are proposed, the main difference being that the gradient regularization constraint is added to the feature mapping graph to further improve the horizontal continuity of the deconvolution result. Finally, taking into account the vertical resolution and horizontal continuity, the advantages of single-channel deconvolution and multi-channel deconvolution are complementary. The vertical and horizontal constraints are added to the deconvolution objective function, so that the deconvolution result has both resolution improvement and continuity enhancement, achieving the purpose of high-resolution inversion of "separating on the vertical axis" and "connecting on the horizontal axis".

[0113] In one embodiment, after the seismic data is subjected to frequency division processing, whether the corresponding components still satisfy the convolution relationship is the key to separate deconvolution processing of each component.

[0114] The frequency division regularization parameter γ is set to 10, and the Marmousi reflection coefficient profile is combined with a 30Hz Ricker wavelet to synthesize a seismic record. According to equation (0-1), the two profiles are subjected to frequency division processing to obtain the corresponding low-frequency components (R LFC and S LFC ) and the corresponding high-frequency components (R HFC and S HFC ). It is assumed that the corresponding components still satisfy the convolution relationship, i.e.:

[0115] S LFC = w * R LFC

[0116] S HFC = w * R HFC (0-12)

[0117] The low-frequency component R LFC and the high-frequency component R HFC of the reflection coefficient are respectively convolved with the 30Hz seismic wavelet, and the obtained corresponding low-frequency and high-frequency synthetic seismic records are shown in Figure 5 , Figure 6 .

[0118] and Figure 1 , Figure 2 and Figure 5 , Figure 6The correlation coefficient and root mean square error between the seismic data components and the corresponding synthesized components are shown in Table 1. It can be seen that the convolution relationship between the corresponding components is still satisfied. This provides a theoretical support for subsequent frequency division deconvolution processing.

[0119] Table 1 Quantitative index analysis between each component of seismic data and the corresponding synthesized component

[0120]

[0121] In one embodiment, for a two-dimensional image, a K-SVD dictionary learning can be used to learn its spatial features. By setting a fixed size sliding window, along the horizontal and vertical directions, sampling according to a given step, a two-dimensional dictionary learning sample is constructed. A sampling diagram of a two-dimensional image is shown in FIG. 1. Figure 7

[0122] In one embodiment, in order to facilitate the use of K-SVD algorithm for dictionary learning, each two-dimensional sample block y i (i = 1, 2, …, n) can be converted into a one-dimensional signal by column merging

[0123]

[0124] wherein represents the kth column data of the ith sample block, and m is the number of all sample blocks.

[0125] In order to extract the features of all sample blocks in the dictionary learning process, the column vector of each block is normalized from 0 to 1 to obtain and all are merged to obtain a two-dimensional sample signal Y:

[0126]

[0127] For the low frequency component S LFC of the seismic data, a two-dimensional sample is constructed, a two-dimensional K-SVD dictionary D is learned, and it is applied to the multi-channel deconvolution regularization constraint term.

[0128] Due to the smoothing characteristics of the low frequency component of the seismic data and the reflection coefficient, a regularization constraint term of L2 norm constraint can be applied in the deconvolution objective function. At the same time, in order to remove the "false high frequency" component introduced in the deconvolution process, a low-pass filtering process is performed on the reflection coefficient to be solved in the dictionary constraint term.

[0129] In a specific implementation, since the two-dimensional dictionary is learned from a two-dimensional block set, in the deconvolution process, multiple blocks are still used for processing, and finally the deconvolution of the entire two-dimensional profile is realized by merging the blocks. The final deconvolution objective function is:​

[0130]

[0131] wherein, and are the low-frequency component patches of the seismic data and the corresponding (to be solved) low-frequency component patches of the reflectivity, respectively, A is the sparse coefficient in the sparse representation process, μ vL and μ hL are the regularization parameters corresponding to the vertical and horizontal constraint terms in the low-frequency component deconvolution objective function, respectively, and λ is the regularization parameter of the sparse constraint term.

[0132] For solving equation (1-4), the following Algorithm 2-1 based on ADMM algorithm can be used to solve it, and this algorithm is called multi-channel seismic deconvolution based on K-SVD.

[0133]

[0134]

[0135] In the above algorithm, two intermediate variables U1 and U2 are introduced, the purpose is to accelerate the convergence speed of ADMM algorithm and improve its stability, by introducing U1 and U2, the solution of the objective function is converted into two independent optimization problems, and these independent optimization problems are easy to solve.

[0136] In one embodiment, after the frequency division processing of the seismic profile, for the high-frequency component S HFC , the convolution dictionary learning of equation (1-10) is used

[0137]

[0138] wherein d i is an n x n two-dimensional convolution dictionary, Z i is a feature map, and the original sample S HFC maintains the same dimension, and N is the number of dictionaries. In this paper, ADMM algorithm is used to solve equation (1-10).

[0139] In order to construct a two-dimensional high-resolution deconvolution system of the profile, the expected output reflectivity is associated with the high-resolution seismic data through a low-pass filter operator, that is, it is made to satisfy the following relationship

[0140] S H = LR HFC (1-11) wherein S H is the expected output high-resolution seismic data (high-frequency component), L is a low-pass filter operator, and R HFC is the expected output two-dimensional reflectivity profile.

[0141] It is important to note the difference between the low-pass filter operator L and the seismic wavelet w. The purpose of equation (1-11) is to limit the spectrum of the obtained full-band reflectivity to a certain range, filter out the "false high frequency" part in R HFC , and the frequency band after convolution by equation (1-11) is expanded compared to the original seismic data S HFC . At the same time, this filter operator is also for the conversion between the "signal domain" (also from the reflectivity domain to the seismic data domain) in order to add a convolution sparsity constraint term based on the seismic data domain in the deconvolution objective function. The convolution of the seismic wavelet w and R HFC will directly obtain S HFC .

[0142] After filtering the desired output, the convolution dictionary is added to the deconvolution objective function in the form of regularization, and equation (1-12) realizes the lateral regularization constraint of the two-dimensional profile

[0143]

[0144] The high-frequency reflectivity itself has a sparse feature, so the L1 norm of the matrix is still used as a constraint term in the vertical direction of the high-frequency deconvolution result to improve its vertical resolution and make the vertical axis "separated". The convolution dictionary regularization constraint term is added in the lateral direction to make the same phase axis "connected". The above two constraints will be added to the deconvolution objective function in the form of regularization to improve the resolution and lateral continuity of the deconvolution result:

[0145]

[0146] Equation (1-13) involves two unknown variables R HFC and the feature mapping matrix Z i (i = 1, 2, …, N), which can be solved in two steps by the ADMM algorithm.

[0147] In the first step, by giving a fixed feature mapping matrix , the intermediate process variable R HFC is obtained by solving the following objective function:

[0148]

[0149] This objective function is a least squares problem with an L1 norm regularization term and a smoothing regularization term, which can be solved using iterative algorithms such as coordinate descent or gradient descent. Let

[0150]

[0151] Equation (1-14) can be expanded as follows:

[0152]

[0153] The gradient of formula (1-16) can be expressed as

[0154]

[0155] Wherein, sign(·) is a sign function. The iterative formula of gradient descent is

[0156]

[0157] Wherein, α is a learning rate, which can be adjusted in a linear or stepwise manner.

[0158] The second step is to obtain R HFC from the previous step, and calculate the feature mapping matrix {Z i} under the given dictionary d i , that is, to optimize the following objective function:

[0159]

[0160] By converting the style of formula (1-19), ADMM algorithm can be used to solve it.

[0161] However, for some dictionary learning-based methods, due to the inaccuracy of the dictionary atoms, it may cause the reconstructed image to lose structure or introduce new features. Similarly, this problem also occurs in convolutional sparse coding. Therefore, scholars have introduced gradient regularization constraints to suppress the emergence of outliers in the dictionary learning process. Adding regularization constraints to the mapping graph is more effective than adding them to the image itself. Therefore, the objective function of the convolutional sparse coding-based seismic data high-frequency component deconvolution is further improved:

[0162]

[0163] Wherein, g x and g y are filter operators for calculating the gradient in the x-axis and y-axis directions, respectively, and τ is a hyperparameter. Similar to the solution method of formula (1-13), by giving the pre-learned dictionary d i , the solution of formula (1-20) can be optimized in two steps.

[0164] The optimization of the first step is similar to the solution of formula (1-14). For the second step, the objective function can be expressed in the following form:

[0165]

[0166] For the newly added gradient constraint term, the convolution operation can be converted into the form of matrix multiplication, that is, let

[0167]

[0168] At this time, the gradient constraint term in formula (1-20) can be expressed as:

[0169]

[0170] By converting the convolution operation into the form of matrix multiplication, that is, introducing the idea of block matrix, formula (1-21) can be expressed as:

[0171]

[0172] Wherein:

[0173]

[0174] In order to solve formula (1-24) using ADMM algorithm, introduce dual variable B, and convert it into the following objective function:

[0175]

[0176] At this time, the solution of the above problem can be converted into solving by using ADMM algorithm. Therefore, the multi-channel seismic deconvolution framework based on convolution sparse coding gradient regularization is shown in algorithm 2-2.

[0177]

[0178]

[0179] In one embodiment, the present application uses three quantitative indicators to evaluate the error between the deconvolution result and the real logging data, and at the same time selects the sparse spike deconvolution (SSD) algorithm as the comparative algorithm of the present application.

[0180] The root mean square error is defined as follows:

[0181]

[0182] Wherein, is the real reflection coefficient at the jth channel, the ith sampling point, and nt and nx represent the number of sampling points and seismic channels respectively. As can be seen from the formula, smaller rmse means better deconvolution result.

[0183] The correlation coefficient represents the similarity between the estimated value and the true value, and the correlation coefficient of the two matrices is defined as follows:

[0184]

[0185] in and This represents the average value of matrices A and B. The correlation coefficient ranges from [0,1]. The larger the value, the more similar the two matrices are.

[0186] Structural similarity is assessed by considering factors such as brightness, contrast, and structure of two images. Its definition formula is as follows:

[0187]

[0188] Where, μ x and μ y Let x and y represent the average values ​​respectively. The constant C1 is introduced to avoid... Instability near 0, σ x and σ y Let x and y represent the standard deviations, respectively, and C2 be a constant.

[0189] In typical image processing based on block operations, the SSIM of the entire image is not calculated. Instead, the image is divided into blocks, and the SSIM of each window is calculated using a sliding window before being averaged. This is called the average SSIM.

[0190]

[0191] Where X and Y represent the reference image (real image) and the processed image, respectively, x k and y k Let M be the image of the k-th local window, and M be the number of local windows in the image.

[0192] In one embodiment, the deconvolution method studied in this chapter is applied to the high-resolution inversion of reflection coefficients from actual seismic data. The data studied comes from a work area in southwestern my country, and its dimensions are 129×110×142 (TimeLine×Crossline×Inline). Figure 8 This demonstrates the seismic wavelet extracted using the well-seismic calibration method. The work area contains 73 wells, and their locations are shown in the diagram. Figure 9 It is worth noting that due to noise and aging hardware during data acquisition, data acquisition was inaccurate. Furthermore, the DDSD-ECJSC method has high requirements for well logging data quality. Therefore, this experiment filtered out 25 wells with a well-seismic calibration correlation coefficient below 75%, and used 47 of these wells for feature extraction. One well was then randomly selected (W). test ) is used for algorithm verification.

[0193] The present application carries out global frequency division processing on seismic data according to convolution sparse coding, and then studies low-frequency component deconvolution processing of seismic data based on two-dimensional K-SVD to learn horizontal features from two-dimensional low-frequency seismic data, and then studies high-frequency component deconvolution processing of seismic data based on convolution sparse coding, and finally, the advantages of single-channel deconvolution and multi-channel deconvolution are complementary, vertical and horizontal constraints are added to the deconvolution objective function, so that the deconvolution result has both resolution improvement and continuity enhancement, realizing the purpose of high-resolution inversion of vertical axis "separated" and horizontal axis "connected".

[0194] (1) The present application learns the reservoir reflectivity characteristics and the mapping relationship with the seismic data characteristics from the logging data, has the characteristics of high vertical resolution, and adds the horizontal distribution law of the seismic data to the deconvolution processing, so that the inversion result presents strong horizontal continuity, and the resolution is significantly enhanced compared with the seismic data.

[0195] (2) The deconvolution result of the present application is better than the model-driven deconvolution method (such as sparse pulse deconvolution method).

[0196] Although the present application is described herein with reference to particular embodiments, it should be understood that these examples are merely illustrative of the principles and applications of the present application. It should therefore be understood that numerous modifications can be made to the illustrative embodiments, and that other arrangements can be devised without departing from the spirit and scope of the present application as defined by the appended claims. It should be understood that the different dependent claims and features described herein can be combined with each other in ways other than those described in the original claims. It should also be understood that features described in connection with a single embodiment can be used in other described embodiments.

Claims

1. A multichannel seismic deconvolution method based on two-dimensional K-SVD and convolutional sparse coding, characterized in that, Includes the following steps: S1. The seismic data underwent global frequency division processing, decomposing the seismic profile into a low-frequency component and a high-frequency component: , Low-frequency components, These are high-frequency components; after global frequency division processing, the corresponding components of the seismic data still satisfy the convolution relationship. S2. Based on the two-dimensional K-SVD, the low-frequency components of seismic data are deconvolved to learn the lateral features from the two-dimensional low-frequency seismic data. By designing a low-pass filter, the reflection coefficients can be sparsely reconstructed using the K-SVD dictionary after passing through the low-pass filter, and the L2 norm constraint is added to the reflection coefficients. S3. Deconvolution processing of high-frequency components of seismic data based on convolutional sparse coding.

2. The multichannel seismic deconvolution method based on two-dimensional K-SVD and convolutional sparse coding according to claim 1, characterized in that, The deconvolution processing of low-frequency components of seismic data based on 2D K-SVD includes the following steps: S21. Use the K-SVD dictionary to learn the spatial features of two-dimensional images. By setting a sliding window, samples are taken along the horizontal and vertical directions according to a predetermined step size to construct samples for two-dimensional dictionary learning. S22, divide each two-dimensional sample into small blocks. Convert to a one-dimensional signal by column merging. ; S23, Column vector for each small block Normalization from 0 to 1 is performed to obtain... and all The two-dimensional sample signal Y is obtained by merging the samples. S24. For the low-frequency components of seismic data Two-dimensional samples are constructed, a two-dimensional K-SVD dictionary D is learned, and it is applied to the multi-channel deconvolution regularization constraint term. In the deconvolution process, multiple small blocks are processed, and the deconvolution of the entire two-dimensional profile is finally achieved by merging the small blocks. The final deconvolution objective function is: ； Where J() is the deconvolution objective function, and D is a two-dimensional K-SVD dictionary. and These are low-frequency component blocks of seismic data and corresponding low-frequency component blocks of reflection coefficients, respectively. A represents the sparsity coefficient in the sparse characterization process. and These are the regularization parameters corresponding to the vertical and horizontal constraint terms in the objective function of the low-frequency component deconvolution, respectively. Let G be the regularization parameter for the sparse constraint term, G be the matrix formed by the seismic wavelets obtained from the seismic records through well-seismic calibration, and L be the low-pass filter operator. The Frobenius norm is such that the matrix With the original data matrix The F-norm of the difference should be as small as possible; S25. Solve the deconvolution objective function using the ADMM algorithm.

3. The multichannel seismic deconvolution method based on two-dimensional K-SVD and convolutional sparse coding according to claim 2, characterized in that, The high-frequency component deconvolution processing of seismic data based on convolutional sparse coding includes the following steps: S31. After frequency division processing of the seismic profile, for the high-frequency components... Learning convolutional dictionaries: ； in, for Two-dimensional convolutional dictionary For the feature map, and the original sample To maintain dimensional consistency, N is the number of words in the dictionary, and the ADMM algorithm is used to solve the problem. S32. To construct a two-dimensional high-resolution deconvolution system for the profile, the desired output reflection coefficient is correlated with the high-resolution seismic data using a low-pass filter operator, such that the desired output reflection coefficient and the high-resolution seismic data satisfy the following relationship: ； in Let L be the high-frequency components of the desired high-resolution seismic data output, and L be the low-pass filter operator. To output a two-dimensional reflection coefficient profile; S33. After filtering the desired output, add the convolution dictionary to the deconvolution objective function in a regularized form to impose lateral regularization constraints on the two-dimensional profile: ； Add constraints to the deconvolution objective function in a regularized manner: ； Where W is the wavelet matrix and the high-frequency component regularization parameter. and Feature mapping regularization parameters ; Given a fixed feature mapping matrix The intermediate process variables are obtained by solving the following objective function. : ； S34. Based on the results obtained in the previous step In a given dictionary In the case of calculating the feature mapping matrix That is, to optimize the following objective function: ； Among them, and These are filter operators that calculate gradients along the x-axis and y-axis, respectively. For hyperparameters; S35. The solution of the objective function is transformed into a solution using the ADMM algorithm.

4. The multichannel seismic deconvolution method based on two-dimensional K-SVD and convolutional sparse coding according to claim 2, characterized in that, Solving the problem using the ADMM algorithm involves the following steps: Input: Low-frequency components of seismic data The multichannel wavelet matrix G, regularization parameters and Low-pass filter operator L and two-dimensional K-SVD dictionary D, maximum number of iterations K, and maximum allowable reconstruction error. Output: Low-frequency component of reflection coefficient ;initialization: Lagrange multipliers in the ADMM algorithm and parameters The current iteration number k=0; ， ， ， ， ， ， ， Determine if the termination condition is met; if so, output... If A is not satisfied, k = k + 1, and continue the loop iteration.

5. The multichannel seismic deconvolution method based on two-dimensional K-SVD and convolutional sparse coding according to claim 3, characterized in that, Solving the problem using the ADMM algorithm involves the following steps: Input: High-frequency components of seismic data Wavelet matrix W, high-frequency component regularization parameters Feature mapping regularization parameters Lagrange parameters in the ADMM algorithm Low-pass filter operator L and convolution dictionary Gradient regularization operator Maximum number of iterations K, maximum allowable reconstruction error ; Output: High-frequency component of reflection coefficient Feature Map ; initialization: Lagrange multipliers in the ADMM algorithm and The current iteration number k=0; ， ， in, For learning rate, , It is a given fixed feature mapping matrix; ， in, These are filter operators that calculate gradients along the x-axis and y-axis, respectively. For hyperparameters, B and C are Lagrange multipliers in the ADMM algorithm; ； ； ； k = k + 1.

Images