Compressed sensing spectrum extension method and system based on joint constraints of Noise-DTW and L1 norm

Abstract

The present invention belongs to the field of seismic data perception technology, and discloses a compressed sensing spectrum extension method and system based on the joint constraints of Noise-DTW and L1 norm. This method obtains the time shift of adjacent seismic traces through Noise-DTW, and then constructs a differential operator to ensure that the reflection coefficients of the corresponding time sampling points of the adjacent seismic traces are differentially calculated; the constructed differential operator is used as a regularization operator for lateral constraints, and the longitudinal reflection coefficient is constrained in combination with the L1 norm, thereby improving the longitudinal resolution of the seismic post-stack profile while ensuring good lateral continuity. The present invention has a significantly improved resolution and better lateral continuity of the inverted profile, verifying the correctness and applicability of the method proposed in the present invention.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data perception, and in particular relates to a compressed sensing spectrum spreading method and system based on the joint constraint of Noise-DTW and L1 norm. Background Art

[0002] During seismic data acquisition, the seismic wavelet undergoes spherical diffusion and attenuation due to the earth's filtering, absorbing some of its energy. This results in a certain absorption of high frequencies, with greater absorption occurring at greater depths. This leads to waveform phase distortion and amplitude attenuation, resulting in lower seismic data resolution. This affects the accuracy of characterizing deep, small structures and thin reservoirs. Therefore, widening the bandwidth and improving seismic data resolution are crucial for seismic exploration.

[0003] After years of development, domestic and foreign scholars have proposed many methods to improve the vertical resolution of seismic data, such as deconvolution, spectral whitening, inverse Q filtering, etc. Deconvolution methods are sensitive to high-frequency noise in seismic data, and are not accurate enough when extracting seismic wavelets, which will reduce the signal-to-noise ratio and have limited ability to improve seismic resolution. Spectral whitening methods can only recover spectral information within the effective frequency band of seismic data and cannot widen the seismic frequency band. In actual applications, the inverse Q filtering method does not have enough accuracy in obtaining the Q value, which limits its application in actual data.

[10] The shortcomings of the above three methods limit their application in actual production.

[0004] Based on the fundamental assumption of seismic exploration, the underground reflection coefficient is a random sequence with a full-bandwidth spectrum. After seismic wavelet filtering, the seismic record S becomes a typical limited-bandwidth spectrum, losing much of the valid signal. Because the compressed sensing theory proposed by Donoho, Candes, and others breaks through the traditional sampling theorem to reconstruct sparse signals, and because the reflection coefficient r satisfies sparsity, nonlinear reconstruction using compressed sensing can, to a certain extent, recover the full-bandwidth energy from the limited-bandwidth spectrum. Zhang et al. constructed a reflection coefficient model and added prior information as a constraint term in the target equation, and finally solved the reflection coefficient by solving the sparse inversion equation through the basis pursuit algorithm; Han Liguo and Zhu Xiangyu et al. proposed the sparse inversion problem, constrained it with the minimum error, reconstructed the full bandwidth spectrum through the limited bandwidth of the original seismic data, compensated for the low-frequency and high-frequency information, and realized two-way spectrum expansion; Chen Zuqing et al. proposed a sparse pulse inversion coefficient spectrum method based on compressed sensing sparse reconstruction, and verified the feasibility and practicality of the method through testing; Yin Xingyao et al. added model constraints to the basis pursuit inversion equation, and combined it with the gradient projection sparse reconstruction algorithm to solve the equation, so that the inversion stability was improved; Song Weiqi et al. used the observation matrix under the constraint of a larger reflection signal on the basis of compressed sensing theory, and normalized the seismic data. The method ultimately enhances weak signals. Liu Yang, addressing the problem of thin interbeds, proposed removing strong shielding layers from the formation before inverting to solve for reflection sparsity, and combined it with broadband wavelets to enhance weak signals. Ding Yan et al. combined an improved broadband Yu-type low-pass shaping filter to compensate for low-frequency seismic signals without affecting the original high-frequency signals. Liu et al. expanded the sparse inversion problem from the frequency domain to the Gaussian frequency domain, improving the resolution of seismic data while ensuring the signal-to-noise ratio. Sun et al. added more accurate logging information as a constraint to the frequency domain sparse inversion equation, effectively compensating for low-frequency information. Zhu Yanhui et al. used the reflection coefficient parity decomposition method and Fourier transform to solve reflection sparsity using the sparse inversion equation. This method effectively broadened the effective frequency band of seismic data without reducing the signal-to-noise ratio. In compressed sensing, methods for solving sparse coefficients have been developed, such as orthogonal matching pursuit, basis pursuit, subspace pursuit, and alternating multiplier direction algorithm. The Alternating Direction Method of Multipliers (ADMM) can be considered a combination of the dual ascent method and the Lagrange multiplier method. As a result, the ADMM algorithm boasts excellent robustness, simplifies complex problems, offers low computational cost, and achieves excellent convergence, earning it widespread attention in seismic data processing. However, the aforementioned compressed sensing-based spectrum spreading methods all ignore lateral spatial continuity and are severely affected by high-frequency noise during inversion, resulting in jitter and discontinuity in the profile events. Summary of the Invention

[0005] To overcome the problems existing in the related art, the embodiments disclosed in the present invention provide a compressed sensing spectrum extension method and system based on the joint constraints of Noise-DTW and L1 norm, specifically relating to a compressed sensing spectrum extension technology based on the joint constraints of Noise-DTW and L1 norm.

[0006] The technical solution is as follows: a compressed sensing spectrum spreading method based on the joint constraints of Noise-DTW and L1 norm, including:

[0007] S1, obtain the time shift of adjacent seismic traces through Noise-DTW, construct a difference operator, and perform a difference operation on the reflection coefficients in adjacent seismic traces; if the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, then the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW;

[0008] S2, using the constructed difference operator as a regularization operator for lateral constraints;

[0009] S3, based on lateral constraints, combined with the L1 norm to constrain the longitudinal reflection coefficient, so that the seismic post-stack profile has lateral continuity, while improving the vertical resolution of the seismic post-stack profile and ensuring lateral continuity.

[0010] In step S1, the time shift of adjacent seismic traces is obtained by using Noise-DTW, including:

[0011] The optimal regularization path of two adjacent seismic data is obtained by the Noise-DTW algorithm:

[0012] u={u1,u2…u m}, so that two data f={f1, f2…f m} and g={g1,g2…g m} matches, and the following conditions are met:

[0013] f[i]≈g[i+u[i]]

[0014] Where u is the optimal regularized path, u m is the time shift corresponding to the mth sampling point, there are m sampling points in total; f and g are both single-channel seismic records, f m , g m are all values ​​corresponding to the mth sampling point, f[i] is the value corresponding to the ith sampling point, g[i+u[i]] is the value corresponding to the i+u[i]th sampling point, i is the ith sampling point, and u[i] is the time shift corresponding to the ith sampling point;

[0015] When the two data are matched, the time shift between the two signals is obtained, and the time shift is used to reconstruct the seismic signal.

[0016] Furthermore, when the two data are matched, the time shift between the two signals is obtained, including:

[0017] Step 1: According to the Noise-DTW algorithm, find the alignment error between the two signals f and g:

[0018] e i，j =(f i -g j ) 2

[0019] Where, e i，j is the element in row i and column j of the alignment error matrix, f i is the i-th element in signal f, g j is the jth signal in signal g;

[0020] Step 2: Calculate the error accumulation matrix based on the alignment error matrix. The expression is:

[0021]

[0022] Where c i，j is the element in row i and column j of the error accumulation matrix, e i，j is the element in the i-th row and j-th column of the alignment error matrix, and min() is the minimum value among them;

[0023] Step 3: In the error accumulation matrix, backtrack to find the regular path u with the minimum cumulative error;

[0024] Step 4: To address the effect of noise on the regularized path, a denoising autoencoder is used to add Gaussian noise with a mean of 0 and a variance of n to the f and g time series. Steps 1 to 3 are repeated for the two noised time series. DTW is performed multiple times with noise added. Each DTW execution generates a regularized path. From these multiple regularized paths, the path with the most repetitions of each step is found to obtain the optimal path.

[0025] The fifth step is to obtain a time shift vector U between adjacent traces by regularizing the path u. U aligns the reflection coefficient positions in the two adjacent seismic data traces. The expression is:

[0026] U×r1=r2

[0027] Where U is the time shift vector, r1 is the reflection coefficient corresponding to the first of two adjacent seismic data, and r2 is the reflection coefficient corresponding to the second of two adjacent seismic data.

[0028] Furthermore, backtracking is used to find the regular path u with the minimum cumulative error, which includes satisfying three conditions during the backtracking process:

[0029] (1) Boundary conditions: u1 = (1, 1) and u m =(m, m);

[0030] (2) Monotonicity condition: The vertical and horizontal coordinates of the regular path are monotonically non-decreasing. When two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≥0 and b1-b2≥0;

[0031] (3) Continuity condition: The two adjacent points in the regular path must be adjacent points in C. When the two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≤1 and b1-b2≤1 are satisfied, ensuring that every coordinate of f and g appears in u.

[0032] Furthermore, a differential operator is constructed, including:

[0033] After aligning the two reflection coefficients, add a horizontal first-order difference operator to obtain the lateral constraint operator obtained by Noise-DTW, which is expressed as:

[0034]

[0035] Where D is the first-order horizontal difference operator, is the multiplication operator, is a new reflection coefficient vector composed of the reflection coefficient vectors of adjacent seismic traces.

[0036] In step S2, the constructed differential operator is used as a regularization operator for lateral constraints, including: multi-channel compressed sensing spectrum spreading by jointly inverting the reflection coefficient of multi-channel seismic data and then spreading the spectrum; expanding the single-channel spectrum spreading into multi-channel spectrum spreading as follows:

[0037]

[0038] Where J(R) is the target equation, is the two-norm, λ r is the regularization parameter, is the multi-channel reflection coefficient in series, is the frequency domain seismic data, is a diagonal seismic wavelet matrix composed of n frequency domain seismic wavelet matrices during single-channel inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrices of n single-channel inversion, N is the number of simultaneous seismic traces, and T is the transpose;

[0039] Bringing the lateral constraint operator into the objective equation, we get the compressed sensing spectrum extension objective equation based on lateral constraint, which is expressed as follows:

[0040]

[0041] In the formula, ‖‖1 is a norm, λ h is the regularization parameter of the lateral constraint, U is the time shift vector, and D is the first-order horizontal difference operator;

[0042] For single-channel frequency extension, add a lateral constraint of L2 norm constraint, use alternating direction multiplier method ADMM to solve, and decompose it into two functions to obtain:

[0043]

[0044] J(z)=λ r ‖z‖1

[0045] Where, for Minimum, J(z), J(m) are the corresponding objective equations, and z is the auxiliary variable introduced in the ADMM algorithm;

[0046] And the two variables m and z need to satisfy the following relationship:

[0047]

[0048] The above equation is the ADMM form of compressed sensing with lateral constraints, which can be converted into the augmented Lagrangian form:

[0049]

[0050] Where τ is a penalty parameter and is greater than 0, controlling the iteration step size; b is an auxiliary variable and has the same magnitude as the reflection coefficient after the combination.

[0051] Furthermore, the alternating direction multiplier method ADMM is used to solve the problem, including:

[0052] Step 1: Solve by solving the gradient:

[0053]

[0054] In the formula, arg is the independent variable;

[0055] Step 2: Solve the problem through soft threshold iteration:

[0056]

[0057] Step 3:

[0058]

[0059] Where b k+1is the Lagrange multiplier obtained at the k+1th iteration, is the reflection coefficient obtained in the k+1th iteration, b k is the Lagrange multiplier obtained at the kth iteration, z k+1 is the auxiliary variable obtained in the k+1th iteration;

[0060] The above three steps are summarized as the following solution:

[0061]

[0062] Where, is the diagonal seismic wavelet matrix composed of two frequency domain seismic wavelet matrices during single-channel inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrices of two single-channel inversions, (UD) T is the transpose of the time shift matrix U and the first-order difference matrix after multiplication, T is the transpose, I is the unit matrix, z k is the auxiliary variable obtained in the kth iteration, b k is the Lagrange multiplier obtained at the kth iteration, and soft() is the soft threshold function;

[0063] The soft threshold function soft() is:

[0064]

[0065] In the formula, x is a variable, β is a threshold, and sign(x) is a function that can extract positive and negative values;

[0066] Repeat the above three steps until the stopping condition is met.

[0067] In step S3, based on the lateral constraint, the longitudinal reflection coefficient is constrained in combination with the L1 norm to ensure the lateral continuity of the seismic post-stack profile, thereby improving the longitudinal resolution of the seismic post-stack profile while ensuring the lateral continuity, including:

[0068] The alternating direction multiplier method ADMM iteratively obtains the reflection coefficient of two adjacent channels. Separate the reflection coefficient of the latter channel and get Through Fourier transform F Transformed into the reflection coefficient spectrum in the frequency domain

[0069] Select the inverted full-bandwidth reflection coefficient spectrum The points with the same amplitude as the original data spectrum S at the same frequency in the mid-high frequency part are extracted. The spectrum in and the spectrum in S before this point are combined into the spectrum after high-frequency extension Then, the seismic data after frequency extension is obtained through inverse Fourier transform.

[0070] Furthermore, the reflection coefficient spectrum The expression is:

[0071]

[0072] Where, is the inverted full-bandwidth reflection coefficient spectrum, which is matched with the original seismic data spectrum for high-frequency extension. F is the Fourier transform matrix. is the reflection coefficient;

[0073] earthquake data Expressed as:

[0074]

[0075] Where, is the earthquake data, F -1 is the inverse Fourier transform, L(S) is the original seismic record spectrum retained, is the high-frequency part of the inverted reflection coefficient spectrum, L is the operator for extracting the original data spectrum value, and H is the operator for extracting the high-frequency spectrum value.

[0076] Another object of the present invention is to provide a compressed sensing spectrum spreading system based on the joint constraint of Noise-DTW and L1 norm, the system comprising:

[0077] A differential operator construction module is used to obtain the time shift of adjacent seismic traces through Noise-DTW and construct a differential operator to perform differential operations on the reflection coefficients in adjacent seismic traces. If the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW.

[0078] A lateral constraint module is used to utilize the constructed differential operator as a regularization operator for lateral constraint;

[0079] The longitudinal reflection coefficient constraint module is used to constrain the longitudinal reflection coefficient based on the lateral constraint and the L1 norm, so as to ensure the lateral continuity of the seismic post-stack profile while improving the longitudinal resolution of the seismic post-stack profile.

[0080] Combining all the above technical solutions, the beneficial effects of the present invention are as follows: in view of the limitation that traditional frequency band compensation methods can only widen the high-frequency part or the low-frequency part of seismic data, a compressed sensing extension technology has been developed that can widen both the low-frequency and high-frequency parts at the same time. However, existing compressed sensing extension technologies all perform extension on a single channel and rarely consider lateral continuity. In this paper, a compressed sensing extension technology based on a noise-robust dynamic time warping algorithm (Noise-DTW) and L1 norm joint constraint is proposed. First, the time shift between adjacent seismic channels is obtained by Noise-DTW, and the time shift is combined with the differential operator to form a lateral constraint for multi-channel inversion. Then, the constructed lateral constraint is introduced as a regularization term into the frequency domain compressed sensing extension target equation based on L1 norm constraint. This can solve the problem of poor lateral continuity of the event axis during traditional single-channel extension, while achieving the purpose of improving the resolution of seismic data. The above target equation can be solved by the alternating direction multiplier method (ADMM). The test results of model data and actual seismic data show that the resolution is significantly improved and the lateral continuity of the inverted profile is better, which verifies the correctness and applicability of the method proposed in this invention.

[0081] This invention performs high-frequency and low-frequency compensation on seismic records, improving their resolution and the ability to identify deep, small structures. Single-channel inversion accuracy is easily affected by noise, while multi-channel inversion can enhance the robustness of the inversion. This method can also address the poor lateral continuity of events in traditional single-channel frequency spreading, increasing the vertical resolution of seismic data while strengthening the lateral continuity of the profile. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present disclosure and, together with the description, serve to explain the principles of the present disclosure;

[0083] Figure 1 This is a flow chart of a compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm provided by an embodiment of the present invention;

[0084] Figure 2 The 1-DTW step size constraint provided by the embodiment of the present invention, i.e., the 1-DTW step size constraint and the planning path graph;

[0085] Figure 3 The Noise-DTW optimal path diagram provided by the embodiment of the present invention is selected;

[0086] Figure 4 is an original data diagram of a synthetic seismic record provided by an embodiment of the present invention;

[0087] Figure 5Reflection coefficient diagram of synthetic seismic record obtained by convolution with 15 Hz Ricker wavelet with sampling interval of 2 ms;

[0088] Figure 6 The single-channel seismic data and corresponding reflection coefficient maps corresponding to the first and second channels extracted from the profile;

[0089] Figure 7 Single-channel seismic record diagram after adding 10% noise to the first and second trajectory data extracted from the synthetic seismic data;

[0090] Figure 8 This is the reflection coefficient diagram of the time shift after adding 10% noise;

[0091] Figure 9 Single-channel seismic record diagram after adding 15% noise to the first and second trajectory data extracted from the synthetic seismic data;

[0092] Figure 10 Reflection coefficient diagram of time shift after adding 15% noise;

[0093] Figure 11 Single-channel seismic record diagram after adding 20% ​​noise to the first and second trajectory data extracted from the synthetic seismic data;

[0094] Figure 12 The reflection coefficient diagram after adding 20% ​​noise as a function of time shift. DETAILED DESCRIPTION

[0095] To make the above-mentioned objects, features, and advantages of the present invention more readily apparent, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. The following description sets forth numerous specific details to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0096] Existing compressed sensing spectrum spreading technologies all perform spectrum spreading on a single channel, and seldom consider lateral continuity. In the present invention, a compressed sensing spectrum spreading technology based on a noise-robust dynamic time warping algorithm (Noise-DTW) and L1 norm joint constraint is proposed. The innovation of the present invention is that: first, the time shift between adjacent seismic channels is obtained by Noise-DTW, and the time shift is combined with the differential operator to form a lateral constraint for multi-channel inversion. The constructed lateral constraint is then introduced as a regularization term into the frequency domain compressed sensing spectrum spreading target equation based on L1 norm constraint. This solves the problem of poor lateral continuity of the event axis in traditional single-channel spectrum spreading, while achieving the purpose of improving the resolution of seismic data.

[0097] Example 1, as Figure 1 As shown, the embodiment of the present invention provides a compressed sensing extension method based on the joint constraint of Noise-DTW and L1 norm. Based on the compressed sensing theory, a multi-channel compressed sensing extension method is proposed. It uses the L1 norm to constrain the reflection coefficient while applying a lateral constraint operator. Specifically, it includes:

[0098] S1, obtain the time shift of adjacent seismic traces through Noise-DTW, construct a difference operator, and perform a difference operation on the reflection coefficients in adjacent seismic traces; if the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, then the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW;

[0099] S2, using the constructed difference operator as a regularization operator for lateral constraints;

[0100] S3, based on lateral constraints, combined with the L1 norm to constrain the longitudinal reflection coefficient, so that the seismic post-stack profile has lateral continuity, while improving the vertical resolution of the seismic post-stack profile and ensuring lateral continuity.

[0101] Example 2, as a preferred embodiment, the compressed sensing spectrum extension method based on the joint constraint of Noise-DTW and L1 norm provided by the embodiment of the present invention includes:

[0102] 1.1 Based on the principle of single-channel inversion reflection coefficient in frequency domain compressed sensing.

[0103] According to the seismic convolution model, the original seismic record s is expressed as the convolution of the frequency domain seismic wavelet W and the reflection coefficient r plus the sum of random noise n:

[0104] s=w*r+n (1)

[0105] Where * is the convolution operator;

[0106] By transforming to the frequency domain through Fourier transform, the expression of formula (1) becomes:

[0107] S=W*R=WFr+E (2)

[0108] Where E is the frequency domain noise;

[0109] Where: S, W, R, and F represent the matrix representation of the actual seismic record, seismic wavelet sequence, and reflection coefficient sequence in the frequency domain and the Fourier transform operator, respectively. Formula (2) conforms to the structure of the compressed sensing theoretical expression and the reflection coefficient sequence meets the sparsity condition of compressed sensing. The seismic wavelet sequence W and the Fourier transform operator F are unrelated. The present invention uses the L1 norm to represent the "sparse" characteristics of the signal to be solved, so the objective function of compressed sensing can be constructed as:

[0110]

[0111] Where: is the error function constrained by the L2 norm, and the error is continuously reduced during the iteration process until it reaches the threshold. It is an accurate solution to the L1 norm constraint, and the reflection coefficients r and λ with better sparsity are continuously obtained in the iteration. r is a regularization parameter based on prior information, which is used to adjust the weights of the two norms to ensure a balance between the sparsity of the reflection coefficient and data fidelity.

[0112] 1.2 Noise-DTW algorithm principle

[0113] Dynamic Time Warping (DTW) is mainly used for matching between two time series. It was originally used in speech recognition. The output of dynamic time warping is a series of time shifts that can be used to align two time series. This algorithm is applied to the field of well logging and seismic joint processing technology, proving that DTW accurately estimates the time shift between two time series. The DTW algorithm can be used to obtain the optimal warping path u = {u1, u2…u m} Make two data f={f1,f2…f m} and g={g1,g2…g m} can be best matched, that is, it satisfies:

[0114] f[i]≈g[i+u[i]]

[0115] Where u is the optimal regularized path, u m is the time shift corresponding to the mth sampling point, there are m sampling points in total; f and g are both single-channel seismic records, f m , g m are all values ​​corresponding to the mth sampling point, f[i] is the value corresponding to the ith sampling point, g[i+u[i]] is the value corresponding to the i+u[i]th sampling point, i is the ith sampling point, and u[i] is the time shift corresponding to the ith sampling point;

[0116] When two data are matched, the time shift between the two signals is obtained, and this time shift is used to reconstruct the seismic signal. Specifically, it includes:

[0117] Step 1: First calculate the alignment error between the two signals f and g according to the DTW algorithm:

[0118] e i，j =(f i -g j ) 2 (5)

[0119] Where, e i，j is the element in row i and column j of the alignment error matrix, f i is the i-th element in signal f, g j is the jth signal in signal g;

[0120] Step 2: Then calculate the error accumulation matrix based on the alignment error matrix:

[0121]

[0122] Where c i，j is the element in row i and column j of the error accumulation matrix, e i，j is the element in the i-th row and j-th column of the alignment error matrix, and min() is the minimum value among them;

[0123] Step 3: In the error accumulation matrix, backtrack to find the regular path u with the minimum cumulative error;

[0124] Step 3: Backtrack in the error accumulation matrix obtained in the previous step to find the regular path u with the minimum cumulative error. During the backtracking process, three conditions must be met:

[0125] (1) Boundary conditions: u1 = (1, 1) and u m =(m, m);

[0126] (2) Monotonicity condition: The vertical and horizontal coordinates of the regular path are monotonically non-decreasing. When two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≥0 and b1-b2≥0;

[0127] (3) Continuity condition: The two adjacent points in the regular path must be adjacent points in C. When the two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≤1 and b1-b2≤1 are satisfied, ensuring that every coordinate of f and g appears in u.

[0128] Considering the three properties of regular paths when backtracking, the backtracking path can only come from three directions. If the path passes through point (i, j), then the previous point position may be the left point (i-1, j), the bottom side (i, j-1), or the bottom left side (i-1, j-1). Figure 2 The 1-DTW step size constraint is shown, i.e., the 1-DTW step size constraint and the planned path.

[0129] Step 4: Considering the impact of noise on the regular path, we draw on the idea of ​​the denoising autoencoder algorithm. Add Gaussian noise with a mean of 0 and a variance of n to the f and g time series, and repeat steps 1 to 3 for the two time series after adding noise. Perform DTW multiple times after adding noise. Each time DTW is performed, a regular path is obtained. From the multiple regular paths, the path with the most repetitions of each step is found to obtain the optimal path, thereby increasing the noise resistance when calculating the regular path. Figure 3 The Noise-DTW provided selects the optimal path as shown. After adding noise multiple times, DTW is executed. Each time DTW is executed, a regular path is obtained. From multiple regular paths, the path with the most repetitions per step is found to obtain the optimal path.

[0130] In the fifth step, the present invention can obtain a regularized path u through the first to fourth steps. A time shift vector U between adjacent traces can be obtained through the regularized path u. U can align the reflection coefficient positions in the seismic data of two adjacent traces, that is:

[0131] U×r1=r2 (7)

[0132] Where U is the time shift vector, r1 is the reflection coefficient corresponding to the first of two adjacent seismic data, and r2 is the reflection coefficient corresponding to the second of two adjacent seismic data.

[0133] After aligning the two reflection coefficients through equation (7), adding a horizontal first-order difference operator to the above equation can obtain the lateral constraint operator obtained by Noise-DTW:

[0134]

[0135] Where D is the first-order horizontal difference operator, is the multiplication operator, is a new reflection coefficient vector composed of the reflection coefficient vectors of adjacent seismic traces.

[0136] Figure 4 represents synthetic seismic record data, Figure 5 Represents the reflection coefficient of the synthetic seismic record obtained by convolution with a 15 Hz Ricker wavelet with a sampling interval of 2 ms. Figure 6The single-channel seismic data and corresponding reflection coefficients for the first and second traces extracted from the profile are shown. It can be seen that the time shifts corresponding to the sampling points of the reflection coefficients are identical to those in the single-channel seismic data. Therefore, the time shifts between adjacent seismic traces can be obtained using Noise-DTW. This time shift can then be applied to adjacent reflection coefficient sampling points to align the reflection coefficients at the same level. This can then be combined with the first-order difference matrix to create lateral constraints, enhancing the lateral continuity of the inverted profile.

[0137] Figure 7 、 Figure 9 、 Figure 11 The first and second trajectory data are extracted from the synthetic seismic data, and 10%, 15%, and 20% noise is added to the single-channel seismic records. Noise-DTW is used to obtain the time shift from the noised single-channel seismic records and apply it to the reflection coefficient to obtain the corresponding time shift. Figure 8 、 Figure 10 、 Figure 12 , are respectively the reflection coefficients of the time shift after adding 10% noise, the reflection coefficients of the time shift after adding 15% noise, and the reflection coefficients of the time shift after adding 20% ​​noise;

[0138] from Figure 8 、 Figure 10 、 Figure 12 It can be seen that after adding 10% noise, the time shift obtained by Noise-DTW is relatively accurate, ensuring a one-to-one correspondence between the two reflection coefficients, which is beneficial for combining with the differential operator to create lateral constraints. However, when the noise exceeds 15%, the time shift obtained by Noise-DTW becomes inaccurate due to the noise, affecting the accuracy of the lateral constraints. The greater the noise, the greater the impact on the obtained time shift.

[0139] 1.3 Principle of Noise-DTW-based Time-Shift Picking for Lateral Constraint Inversion of Reflection Coefficients

[0140] Multi-channel compressed sensing frequency spreading is achieved by inverting the reflection coefficients of multiple seismic channel data together and then spreading the frequency. The single-channel frequency spreading in formula (3) is expanded to multi-channel frequency spreading as follows:

[0141]

[0142] Where J(R) is the target equation, is the two-norm, λ r is the regularization parameter, is the multi-channel reflection coefficient in series, is the frequency domain seismic data, is a diagonal seismic wavelet matrix composed of n frequency domain seismic wavelet matrices during single-channel inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrices of n single-channel inversion, N is the number of simultaneous seismic traces, and T is the transpose;

[0143] In the present invention, n is 2, and the reflection coefficient of the next pass is solved using the reflection coefficient constraint in the previous pass.

[0144] Substituting the lateral constraint operator (8) into the target equation (9), the target equation of compressed sensing spectrum extension based on lateral constraint is obtained as follows:

[0145]

[0146] In the formula, ‖‖1 is a norm, λ h is the regularization parameter of the lateral constraint, U is the time shift vector, and D is the first-order horizontal difference operator;

[0147] Compared with single-channel spreading, formula (10) only adds a lateral constraint of the L2 norm constraint. The ADMM algorithm can be used to solve the above formula. The ADMM algorithm first decouples the above formula and decomposes it into two functions to obtain:

[0148]

[0149] J(z)=λ r ‖z‖1

[0150] Where, for Minimum, J(z), J(m) are the corresponding objective equations, and z is the auxiliary variable introduced in the ADMM algorithm;

[0151] And the two variables m and z need to satisfy the following relationship:

[0152]

[0153] The above equations (11) and (12) are the ADMM forms of the compressed sensing spectrum extension (10) based on the lateral constraint, which can be converted into the augmented Lagrangian form:

[0154]

[0155] Where τ is a penalty parameter and is greater than 0, controlling the iteration step size, and b is an auxiliary variable, which is the same as the reflection coefficient after the combination. The solution using the ADMM algorithm is divided into three steps:

[0156] Step 1: Solve by solving the gradient:

[0157]

[0158] In the formula, arg is the independent variable;

[0159] Step 2: Solve the problem through soft threshold iteration:

[0160]

[0161] Step 3:

[0162]

[0163] Where b k+1 is the Lagrange multiplier obtained at the k+1th iteration, is the reflection coefficient obtained in the k+1th iteration, b k is the Lagrange multiplier obtained at the kth iteration, z k+1 is the auxiliary variable obtained in the k+1th iteration;

[0164] The above three steps can be generalized into the following solution:

[0165]

[0166] Where, is the diagonal seismic wavelet matrix composed of two frequency domain seismic wavelet matrices during single-channel inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrices of two single-channel inversions, (UD) T is the transpose of the time shift matrix U and the first-order difference matrix after multiplication, T is the transpose, I is the unit matrix, z k is the auxiliary variable obtained in the kth iteration, b k is the Lagrange multiplier obtained at the kth iteration, and soft() is the soft threshold function;

[0167] The soft threshold function (Soft Tresholding) is:

[0168]

[0169] In the formula, x is a variable, β is a threshold, and sign(x) is a function that can extract positive and negative values;

[0170] Repeat the above three steps until the stopping condition is met.

[0171] 1.4 Bandwidth broadening formula.

[0172] The alternating direction multiplier method ADMM iteratively obtains the reflection coefficient of two adjacent channels. Separate the reflection coefficient of the latter channel and get Through Fourier transform F Transformed into the reflection coefficient spectrum in the frequency domain Right now:

[0173]

[0174] Where, is the inverted full-bandwidth reflection coefficient spectrum, which is matched with the original seismic data spectrum for high-frequency extension. F is the Fourier transform matrix. is the reflection coefficient; the specific method is as follows:

[0175] Select the inverted full-bandwidth reflection coefficient spectrum The points with the same amplitude as the original data spectrum S at the same frequency in the mid-high frequency part are extracted. The spectrum in and the spectrum in S before this point are combined into the spectrum after high-frequency extension Then, the seismic data after frequency extension is obtained through inverse Fourier transform. The formula is:

[0176]

[0177] Where, is the earthquake data, F -1 is the inverse Fourier transform, L(S) is the original seismic record spectrum retained, is the high-frequency part of the inverted reflection coefficient spectrum, L is the operator for extracting the original data spectrum value, and H is the operator for extracting the high-frequency spectrum value. The low-frequency compensation method is similar.

[0178] In Example 3, the compressed sensing spectrum extension system based on the joint constraint of Noise-DTW and L1 norm provided by the embodiment of the present invention is based on the compressed sensing theory and uses L1 norm to constrain the reflection coefficient while applying a lateral constraint operator. Specifically, it includes:

[0179] A differential operator construction module is used to obtain the time shift of adjacent seismic traces through Noise-DTW and construct a differential operator to perform differential operations on the reflection coefficients in adjacent seismic traces. If the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW.

[0180] A lateral constraint module is used to utilize the constructed differential operator as a regularization operator for lateral constraint;

[0181] The longitudinal reflection coefficient constraint module is used to constrain the longitudinal reflection coefficient based on the lateral constraint and the L1 norm, so as to ensure the lateral continuity of the seismic post-stack profile while improving the longitudinal resolution of the seismic post-stack profile.

[0182] In order to verify the accuracy of the method of the present invention, Figure 4 As shown in Figure 5, 10% random noise is added to the original seismic profile, resulting in a noisy seismic profile with 10% random noise. It can be seen that the profile resolution is low, deep interlayers cannot be identified, and the details of the profile are not accurately depicted. In this invention, it is assumed that the seismic wavelet is known. Formula (10) is used to implement lateral constraint-based compressed sensing frequency spectrum spreading (LCSBF-DTW), and the ADMM algorithm is used to solve the target equation.

[0183] The above description is only a preferred specific implementation method of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with this technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered by the scope of protection of the present invention.

Claims

1. A compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm, characterized in that: The method includes: S1, obtain the time shift of adjacent seismic traces through Noise-DTW, construct a difference operator, and perform a difference operation on the reflection coefficients in adjacent seismic traces; if the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, then the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW; S2, using the constructed difference operator as a regularization operator for lateral constraints; S3, based on the lateral constraint, combined with the L1 norm to constrain the longitudinal reflection coefficient, so as to make the seismic post-stack section lateral continuity, while improving the longitudinal resolution of the seismic post-stack section to ensure the lateral continuity; In step S2, the constructed differential operator is used as a regularization operator for lateral constraints, including: multi-channel compressed sensing spectrum spreading by jointly inverting the reflection coefficient of multi-channel seismic data and then spreading the spectrum; expanding the single-channel spectrum spreading into multi-channel spectrum spreading as follows: Where J(R) is the target equation, is the two-norm, λ r is the regularization parameter, is the multi-channel reflection coefficient in series, is the frequency domain seismic data, is the diagonal seismic wavelet matrix composed of the frequency domain seismic wavelet matrix during inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrix during inversion, N is the number of simultaneous seismic traces, and T is the transpose; Bringing the lateral constraint operator into the objective equation, we get the compressed sensing spectrum extension objective equation based on lateral constraint, which is expressed as follows: In the formula, || ||1 is a norm, λ h is the regularization parameter of the lateral constraint, U is the time shift vector, and D is the first-order horizontal difference operator; For single-channel frequency extension, add a lateral constraint of L2 norm constraint, use alternating direction multiplier method ADMM to solve, and decompose it into two functions to obtain: J(z)=λ r ||of||1 Where, for minimum, J(z), J(m) are the corresponding objective equations, and z is the auxiliary variable introduced in the ADMM algorithm; And the two variables m and z need to satisfy the following relationship: The above equation is the ADMM form of compressed sensing with lateral constraints, which can be converted into the augmented Lagrangian form: Where τ is a penalty parameter and is greater than 0, controlling the iteration step size; b is an auxiliary variable and has the same magnitude as the reflection coefficient after the combination.

2. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 1, characterized in that: In step S1, the time shift of adjacent seismic traces is obtained by using Noise-DTW, including: The optimal regularization path of two adjacent seismic data is obtained by the Noise-DTW algorithm: u={u1,u2…u m }, so that two data f={f1, f2…f m } and g={g1,g2,…g m } matches, and the following conditions are met: f[i]≈g[i+u[i]] Where u is the optimal regularized path, u m is the time shift corresponding to the mth sampling point, there are m sampling points in total; f and g are both single-channel seismic records, f m , g m are all values ​​corresponding to the mth sampling point, f[i] is the value corresponding to the ith sampling point, g[i+u[i]] is the value corresponding to the i+u[i]th sampling point, i is the ith sampling point, and u[i] is the time shift corresponding to the ith sampling point; When the two data are matched, the time shift between the two signals is obtained, and the time shift is used to reconstruct the seismic signal.

3. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 2, characterized in that: When the two data are matched, the time shift between the two signals is obtained, including: Step 1: According to the Noise-DTW algorithm, find the alignment error between the two signals f and g: e i,j =(f i -g j ) 2 Where, e i,j is the element in row i and column j of the alignment error matrix, f i is the i-th element in signal f, g j is the jth signal in signal g; Step 2: Calculate the error accumulation matrix based on the alignment error matrix. The expression is: Where c i,j is the element in row i and column j of the error accumulation matrix, e i，j is the element in row i and column j of the alignment error matrix, and min() is the minimum value among them; Step 3: In the error accumulation matrix, backtrack to find the regular path u with the minimum cumulative error; Step 4: To address the effect of noise on the regularized path, a denoising autoencoder is used to add Gaussian noise with a mean of 0 and a variance of n to the f and g time series. Steps 1 to 3 are repeated for the two noised time series. DTW is performed multiple times with noise added. Each DTW execution generates a regularized path. From these multiple regularized paths, the path with the most repetitions of each step is found to obtain the optimal path. The fifth step is to obtain a time shift vector U between adjacent traces by regularizing the path u. U aligns the reflection coefficient positions in the two adjacent seismic data traces. The expression is: U×r1＝r2 Where U is the time shift vector, r1 is the reflection coefficient corresponding to the first of two adjacent seismic data, and r2 is the reflection coefficient corresponding to the second of two adjacent seismic data.

4. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 3, characterized in that: Backtracking to find the regular path u with the minimum cumulative error includes: satisfying three conditions during the backtracking process: (1) Boundary conditions: u1 = (1, 1) and u m =(m, m); (2) Monotonicity condition: The vertical and horizontal coordinates of the regular path are monotonically non-decreasing. When two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≥0 and b1-b2≥0; (3) Continuity condition: The two adjacent points in the regular path must be adjacent points in C. When the two adjacent points are u i =(a1, b1) and u i-1 =(a2, b2), then a1-a2≤1 and b1-b2≤1 are satisfied, ensuring that every coordinate of f and g appears in u.

5. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 1, characterized in that: Construct a difference operator, including: After aligning the two reflection coefficients, add a horizontal first-order difference operator to obtain the lateral constraint operator obtained by Noise-DTW, which is expressed as: Where D is the first-order horizontal difference operator, U is the time shift vector, is the multiplication operator, is the multi-channel reflection coefficient in series.

6. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 1, characterized in that: Solve using the alternating direction multiplier method ADMM, including: Step 1: Solve by solving the gradient: In the formula, arg is the independent variable; Step 2: Solve the problem through soft threshold iteration: Step 3: Where b k+1 is the Lagrange multiplier obtained at the k+1th iteration, is the reflection coefficient obtained in the k+1th iteration, b k is the Lagrange multiplier obtained at the kth iteration, z k+1 is the auxiliary variable obtained in the k+1th iteration; The above three steps are summarized as the following solution: Where, is the diagonal seismic wavelet matrix composed of the frequency domain seismic wavelet matrix during inversion, is the diagonal Fourier transform matrix composed of the Fourier transform matrix during inversion, (UD) T is the transpose of the time shift matrix and the first-order difference matrix after multiplication, T is the transpose, I is the identity matrix, z k is the auxiliary variable obtained in the kth iteration, b k is the Lagrange multiplier obtained at the kth iteration, and soft() is the soft threshold function; The soft threshold function soft( ) is: In the formula, x is a variable, β is a threshold, and sign(x) is a function that can extract positive and negative values; Repeat the above three steps until the stopping condition is met.

7. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 1, characterized in that: In step S3, based on the lateral constraint, the longitudinal reflection coefficient is constrained in combination with the L1 norm to ensure the lateral continuity of the seismic post-stack profile, thereby improving the longitudinal resolution of the seismic post-stack profile while ensuring the lateral continuity, including: The multi-channel reflection coefficients obtained by the alternating direction multiplier method ADMM are the series connection of two adjacent channels. Separate the reflection coefficient of the latter channel and get Through Fourier transform F Transformed into the reflection coefficient spectrum in the frequency domain Select the inverted full-bandwidth reflection coefficient spectrum The points with the same amplitude as the original data spectrum S at the same frequency in the mid-high frequency part are extracted. The spectrum in and the spectrum in S before this point are combined into a high-frequency extended spectrum, and then the extended seismic data is obtained by inverse Fourier transform.

8. The compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to claim 7, characterized in that: Reflection coefficient spectrum The expression is: Where, is the inverted full-bandwidth reflection coefficient spectrum, which is amplitude-matched with the original seismic data spectrum for high-frequency extension, and F is the Fourier transform matrix; earthquake data Expressed as: Where, is the earthquake data, F -1 is the inverse Fourier transform, L(S) is the original seismic record spectrum retained, is the high-frequency part of the inverted reflection coefficient spectrum, L is the operator for extracting the original data spectrum value, and H is the operator for extracting the high-frequency spectrum value.

9. A compressed sensing spectrum spreading system based on the joint constraint of Noise-DTW and L1 norm, characterized in that: The system is implemented by the compressed sensing spectrum spreading method based on the joint constraint of Noise-DTW and L1 norm according to any one of claims 1 to 8, and the system includes: A differential operator construction module is used to obtain the time shift of adjacent seismic traces through Noise-DTW and construct a differential operator to perform differential operations on the reflection coefficients in adjacent seismic traces. If the inverted reflectivity sequence is convolved with the same seismic wavelet used in the inversion process, the adjacent inverted reflection coefficients conform to the spatial relationship picked up from the adjacent original seismic traces using Noise-DTW. A lateral constraint module is used to utilize the constructed differential operator as a regularization operator for lateral constraint; The longitudinal reflection coefficient constraint module is used to constrain the longitudinal reflection coefficient based on the lateral constraint and the L1 norm, so as to ensure the lateral continuity of the seismic post-stack profile while improving the longitudinal resolution of the seismic post-stack profile.

Images