Method for processing and identifying thin interbeds with high resolution based on compressive sensing with structural constraints

Abstract

The application provides a thin interbed high-resolution processing and identification method based on structural constraint compressed sensing, and improves the resolution of thin layers and thin interbeds by structural constraint; on the basis of traditional compressed sensing for improving seismic resolution, the problems that the spatial connection of seismic data is ignored in the actual processing process and the lateral continuity is difficult to maintain due to the interference of lateral noise are considered; structural information existing in a local structure tensor of seismic profile data is used; a feature vector of the structure tensor is introduced to establish a structural constraint term; and subsequent sparse inversion of seismic reflection coefficients and convolution with a wide-frequency wavelet are completed; frequency width compensation is reasonably performed to improve the main frequency; a high-resolution seismic profile based on structural constraint of the structure tensor is obtained; the applicability of the compressed sensing high-resolution method in a complex geological structure environment is improved to improve the continuity of the processing result; and therefore the identification and prediction ability of the thin layer is effectively improved.

Description

Technical Field

[0001] This invention relates to the field of seismic data target processing and reservoir description, and in particular to a compressed sensing target processing and identification method based on structural tensor tectonic constraints for thin interbedded layers. Background Technology

[0002] With the continuous innovation of oil and gas exploration and development technologies, exploration targets are gradually shifting from conventional to unconventional, from thick to thin, and from continuous to discontinuous reservoirs. The research focus has also shifted to thinner, deeper, and more complex concealed reservoirs. However, the resolution of conventional seismic data is often insufficient to meet the requirements for detailed identification and description of thin reservoirs. At the same time, seismic data is also affected by noise and complex structures, making it difficult to maintain and control its lateral continuity. Conventional frequency upscaling techniques are insufficient to meet the needs of detailed interpretation. Summary of the Invention

[0003] Based on the above-mentioned technical problems, this invention proposes a high-resolution processing and identification method for thin interbedded layers based on structurally constrained compressed sensing. This method is mainly used to improve the seismic resolution of thin interbedded layers while protecting their structural continuity, thereby enhancing the accuracy of subsequent thin reservoir identification and prediction, and providing support for actual seismic production work.

[0004] This application provides a high-resolution processing and recognition method for thin interlayers based on construction-constrained compressed sensing, including:

[0005] S1: Basic information about the seismic data is obtained by performing spectral and waveform feature analysis, and then the wavelet of the original seismic record is extracted using the complex spectrum method;

[0006] S2: Calculate the local structure tensor for the original seismic data, and perform subsequent eigenvalue decomposition on the local structure tensor to obtain the eigenvector at each point of the data.

[0007] S3: Based on the seismic wavelet obtained in step S1, construct the frequency domain diagonal wavelet measurement matrix and further sensing matrix;

[0008] S4: Based on the eigenvectors at the sample points obtained in step S2, combine them with the established horizontal and vertical difference operator matrices to construct the constraint matrix;

[0009] S5: After the multichannel seismic data is rearranged into columns, it is subjected to Fourier transform to obtain frequency domain data. The perception matrix of the entire data is calculated according to step S3 and the construction constraint matrix is ​​obtained according to step S4. Then, the sparse inversion of reflection coefficients is performed using the fast threshold shrinkage iterative algorithm.

[0010] S6: Perform multi-channel cyclic processing on all seismic data to obtain sparse inversion results of reflection coefficients corresponding to all seismic data, and then convolve them with a broadband seismic wavelet to finally obtain seismic data with improved resolution under structural constraints.

[0011] Step S1 further includes the following steps: adjusting the obtained wavelet bandwidth and phase based on the bandwidth and phase of the data and the sliding time window.

[0012] In step S1, the wavelet of the original seismic record is extracted using the complex spectrum method. The specific steps are as follows: First, a time window range is selected for the original seismic data. Then, the seismic traces within this range are subjected to Fourier transform to obtain their amplitude spectrum. Next, the amplitude spectrum is subjected to logarithmic transform and inverse Fourier transform to obtain the complex spectrum result. A low-pass filter is designed by setting the width of the passband to filter the complex spectrum result and perform Fourier and inverse logarithmic transforms. Finally, the extracted wavelet of the original seismic record is obtained.

[0013] In step S2, the local structure tensor and eigenvectors of the original seismic data are calculated. The specific steps are as follows:

[0014] First, the gradient at each sample point within the selected profile data needs to be calculated. As shown in equation (1), then perform smooth outer product processing to obtain the local structure tensor T, as shown in equation (2). After that, perform eigenvalue decomposition on the local structure tensor T at each homogeneous point to obtain the special vectors u and v, as shown in equation (3).

[0015]

[0016]

[0017] In the formula, To calculate the gradient vector at the sample point, <·> is the smoothing operator, and g x For the longitudinal partial derivative, g y Let be the lateral partial derivative, and u and v be the corresponding eigenvectors.

[0018] Step S2 further includes the following step: smoothing the estimated gradient.

[0019] In step S3:

[0020] Based on the seismic wavelet calculated in step S1, construct the frequency domain diagonal wavelet measurement matrix W, as shown in equation (4), and combine it with the partial Fourier basis matrix D, as shown in equation (5), to further obtain the attenuation sensing matrix A, A = GW;

[0021]

[0022] In the formula, W is the measurement matrix, diag is the diagonal function, FT is the forward Fourier transform, w(t) is the estimated wavelet, D is the partial Fourier basis matrix, i represents the imaginary number, f is the frequency, and τ is the time.

[0023] In step S4:

[0024] The feature vector u = (u) at the sample point calculated in step S2 x ,u y ), v = (v x ,v y ), and the established horizontal and vertical difference operator matrices D x D y Combined, we can further construct the constraint matrix s, as shown in equation (6);

[0025]

[0026] In the formula, u and v are eigenvectors, and D x D is the horizontal difference operator matrix. y Let h be the vertical difference operator matrix. u h v These are the constraint parameters for the eigenvectors u and v, respectively, with values ​​between 0 and 1.

[0027] In step S5:

[0028] After rearranging the columns of the multichannel seismic data s(t), we obtain As shown in equation (7), perform Fourier transform on it to obtain frequency domain data y, and calculate the perception matrix A of the data according to step S3 and obtain the construction constraint matrix S according to step S4, and construct the inversion objective function as shown in equation (8).

[0029]

[0030] In the formula, A is the attenuation sensing matrix. λ represents the inversion reflection coefficient, y represents the frequency domain seismic data, λ represents the regularization factor, and S represents the construction constraint matrix.

[0031] In step S5, the sparse inversion of the reflection coefficient is performed using a fast threshold shrinkage iterative algorithm.

[0032] In step S6: the broadband seismic wavelet adopts the higher dominant frequency Reck wavelet.

[0033] The high-resolution processing and recognition method for thin interlayers based on construction-constrained compressed sensing in this application has the following advantages:

[0034] This application addresses the improvement of seismic resolution through structural constraints on thin and interbedded thin layers. Building upon traditional compressed sensing methods for improving seismic resolution, it considers the issues of neglecting spatial relationships in seismic data and the difficulty in maintaining lateral continuity due to lateral noise interference during actual processing. It utilizes structural information (such as dip and azimuth) present in the local structural tensor of seismic profile data, introducing eigenvectors of the structural tensor to establish structural constraint terms. This process then completes the sparse inversion of seismic reflection coefficients and convolution with broadband wavelets, and reasonably performs bandwidth compensation to enhance the dominant frequency. The result is a high-resolution seismic profile based on structural constraints using the structural tensor. This improves the applicability of compressed sensing high-resolution methods in complex geological environments and enhances the continuity of processing results, thereby effectively improving the identification and prediction capabilities of thin layers. Attached Figure Description

[0035] Figure 1 This is a flowchart of the high-resolution processing and recognition method for thin interlayers based on structurally constrained compressed sensing according to the present invention.

[0036] Figure 2 is a two-dimensional theoretical seismic forward model diagram in a specific embodiment of the present invention; wherein, Figure 2(a) is a two-dimensional theoretical wave impedance model diagram, and Figure 2(b) is a two-dimensional theoretical reflection coefficient model diagram;

[0037] Figure 3 is a schematic diagram of a two-dimensional forward modeling seismic record and eigenvector ellipse in a specific embodiment of the present invention; wherein, Figure 3(a) is a two-dimensional forward modeling seismic record with 10% noise, and Figure 3(b) is a schematic diagram of the eigenvector ellipse of a two-dimensional forward modeling seismic record with 10% noise.

[0038] Figure 4 shows the results of constructing constrained sparse inversion and high-resolution processing under 10% noise in a specific embodiment of the present invention; wherein, Figure 4(a) is the result of constructing constrained sparse inversion under 10% noise in two dimensions, and Figure 4(b) is the result of high-resolution processing under 10% noise in two dimensions.

[0039] Figure 5 This is an original seismic profile of a series of wells in a specific embodiment of the present invention;

[0040] Figure 6 This is a schematic diagram of the characteristic vector ellipse of the original seismic profile of the connected wells in a specific embodiment of the present invention;

[0041] Figure 7 This is a high-resolution processed seismic profile of a series of wells in a specific embodiment of the present invention;

[0042] Figure 8 This is a comparison diagram of the spectrum before and after the well profile processing in a specific embodiment of the present invention. Detailed Implementation

[0043] The present application will be further described below with reference to the accompanying drawings and embodiments.

[0044] Traditional compressed sensing-based sparse signal sampling and reconstruction theory provides theoretical support for reconstructing broadband sparse reflection coefficients from band-limited seismic data. However, traditional compressed sensing resolution improvement methods are single-channel inversion algorithms that do not consider the spatial relationships between seismic data, making it difficult to maintain lateral continuity and susceptible to adverse effects from lateral noise. This often results in a "noodle-like" appearance in the inversion profile, which is detrimental to subsequent thin interbedded layer identification and reservoir prediction. Therefore, how to effectively improve resolution while further protecting lateral continuity from band-limited seismic data has become an urgent problem to be solved. Researching a high-resolution processing and identification method for thin interbedded layers based on tectonic-constrained compressed sensing has high theoretical significance and application value. To this end, we have invented a new high-resolution processing and identification method for thin interbedded layers based on tectonic-constrained compressed sensing to solve the above technical problems.

[0045] Example 1

[0046] This application presents a high-resolution processing and identification method for thin interlayered seismic data based on structurally constrained compressed sensing, comprising: S1: obtaining basic information of seismic data through spectral and waveform feature analysis, and then extracting wavelets from the original seismic record using the complex spectrum method; S2: calculating the local structure tensor of the original seismic data, and performing subsequent eigenvalue decomposition on the local structure tensor to obtain the eigenvector at each sample point of the data; S3: constructing a frequency domain diagonal wavelet measurement matrix and a further sensing matrix based on the seismic wavelets obtained in step S1; S4: combining the eigenvectors at the sample points obtained in step S2 with the established horizontal and vertical difference operator matrices to establish a structural constraint matrix;

[0047] S5: After the multichannel seismic data is rearranged columnwise, a Fourier transform is performed to obtain frequency domain data. The sensing matrix of the entire data channel is calculated according to step S3, and the construction constraint matrix is ​​obtained in step S4. Then, a fast threshold shrinkage iterative algorithm is used to perform sparse inversion of the reflection coefficients. S6: Subsequently, multichannel cyclic processing is performed on all seismic data to obtain the sparse inversion results of the reflection coefficients corresponding to all seismic data. These results are then convolved with a broadband seismic wavelet to finally obtain structurally constrained seismic data with improved resolution.

[0048] This application improves the seismic resolution of thin interbedded layers while protecting their structural continuity, enhancing the accuracy of subsequent thin reservoir identification and prediction, and providing support for actual seismic production work.

[0049] Example 2

[0050] like Figure 1 As shown, Figure 1The flowchart below shows the high-resolution processing and recognition method for thin interlayers based on construction-constrained compressed sensing according to this application. The specific steps are as follows:

[0051] Step 1: Extract the wavelet w(t) from the original seismic record using the complex spectrum method.

[0052] First, a suitable time window is selected from the original seismic records, ideally choosing a region with a stable phase axis and gradual changes. Then, a Fourier transform is performed on the seismic traces in this region to obtain their amplitude spectra. Next, a logarithmic transform and inverse Fourier transform are performed on these amplitude spectra to obtain the "complex spectrum." A low-pass filter is designed by setting the passband width to filter the "complex spectrum" result and perform Fourier and inverse logarithmic transforms, ultimately yielding the extracted wavelet from the original seismic records. The filter size should not be too large, aiming to preserve the wavelet information at the origin of the complex spectrum while removing reflection coefficient information from non-origin locations.

[0053] Step 2: Calculate the local structural tensor and eigenvectors of the original seismic data.

[0054] First, the gradient at each sample point within the selected profile data needs to be calculated. As shown in equation (1), the local structure tensor T is obtained by smoothing the outer product, as shown in equation (2). Then, the local structure tensor T at each point is decomposed into eigenvalues ​​to obtain the specific vectors u and v, as shown in equation (3).

[0055]

[0056] In the formula, To calculate the gradient vector at the sample point, <·> is the smoothing operator, and g x For the longitudinal partial derivative, g y Let be the lateral partial derivative, and u and v be the corresponding eigenvectors.

[0057] Step 3: Construct the frequency domain diagonal wavelet measurement matrix and sensing matrix.

[0058] Based on the seismic wavelet calculated in step S1, construct the frequency domain diagonal wavelet measurement matrix W, as shown in equation (4), and combine it with the partial Fourier basis matrix D, as shown in equation (5), to further obtain the attenuation sensing matrix A, A = GW;

[0059]

[0060] In the formula, W is the measurement matrix, diag is the diagonal function, FT is the forward Fourier transform, w(t) is the estimated wavelet, D is the partial Fourier basis matrix, i represents the imaginary number, f is the frequency, and τ is the time.

[0061] Step 4: Construct a constraint matrix based on the eigenvectors.

[0062] The feature vector u = (u) at the sample point calculated in step S2 x ,u y ), v = (v x ,v y ), and the established horizontal and vertical difference operator matrices D x D y Combined, we can further construct the constraint matrix s, as shown in equation (6);

[0063]

[0064] In the formula, u and v are eigenvectors, and D x D is the horizontal difference operator matrix. y Let h be the vertical difference operator matrix. u h v These are the constraint parameters for the eigenvectors u and v, respectively, with values ​​between 0 and 1;

[0065] Step 5: After rearranging the columns of the multichannel seismic data s(t), we obtain... As shown in equation (7), perform Fourier transform on it to obtain frequency domain data y, and calculate the perception matrix A of the data according to step S3 and obtain the construction constraint matrix S according to step S4. Construct the inversion objective function as shown in equation (8). Then, the sparse inversion of the reflection coefficient can be performed using the fast threshold shrinkage iterative algorithm.

[0066]

[0067] In the formula, A is the attenuation sensing matrix. λ represents the inversion reflection coefficient, y represents the frequency domain seismic data, λ represents the regularization factor, and S represents the construction constraint matrix.

[0068] Step 6: Then, perform multi-channel cyclic processing on all seismic data to obtain sparse inversion results of reflection coefficients corresponding to all seismic data, and convolve them with a broadband seismic wavelet to finally obtain seismic data with improved resolution under structural constraints.

[0069] This invention addresses the issue of improving the resolution of thin-layered and interbedded thin-layered seismic data by incorporating structural constraints. Building upon traditional compressed sensing methods for improving seismic resolution, it considers the practical challenges of neglecting spatial relationships in seismic data and the susceptibility to lateral noise interference that makes it difficult to maintain lateral continuity. By utilizing structural information (such as dip and azimuth) present in the local structural tensor of seismic profile data, the invention introduces eigenvectors of the structural tensor to establish structural constraint terms. This process then completes the sparse inversion of seismic reflection coefficients and their convolution with broadband wavelets, and appropriately performs bandwidth compensation to enhance the dominant frequency. The result is a high-resolution seismic profile based on structural constraints using the structural tensor. This improves the applicability of compressed sensing high-resolution methods in complex geological environments and enhances the continuity of processing results, thereby effectively improving the identification and prediction capabilities of thin layers.

[0070] Example 3

[0071] The invention was first applied to a theoretical model to test its effectiveness and applicability. Figure 2 shows the two-dimensional velocity model and reflection coefficient model. This model has 100 channels, 2ms sampling, 151 sampling points, and a 25Hz Ricker wavelet. Figure 3 shows the two-dimensional forward seismic record with 10% noise and the corresponding eigenvector ellipse. Figure 4 shows the structurally constrained sparse inversion and high-resolution processing results under noise conditions. The processing results in Figures 3 and 4 show that the invention can accurately achieve sparse inversion of reflection coefficients while maintaining a certain signal-to-noise ratio and good lateral continuity. It can effectively improve the resolution of seismic data while ensuring the structural continuity of the processing results, providing favorable theoretical model support for practical applications.

[0072] Next, the invention was applied to a thin interbedded section in a certain work area, and the original seismic profile of a main seismic survey line of a certain well-connected well was as follows: Figure 5 As shown, there are 381 sampling points and a sampling interval of 2 milliseconds. (From...) Figure 5 It is evident that the original seismic profile has low resolution. Sandstone reservoirs are sensitive to velocity and gamma curve responses, but their seismic response energy is weak. Thin interbedded layers exhibit waveform superposition, making it difficult to distinguish the thin-layer information. A schematic diagram of the calculated eigenvector ellipse of this profile is shown below. Figure 6 As shown, the calculated feature vectors can provide good structural information such as the dip angle and azimuth of the strata, and are used as structural constraints for subsequent sparse inversion. Then, using the structural constraint-based compressed sensing high-resolution processing and identification method for thin interbedded layers of this invention, an improved resolution profile is obtained, as shown... Figure 7 As shown, high-resolution processing improves the ability to identify thin interbedded layers while maintaining the overall structure of the profile faithfully to the original seismic profile. The characteristics of thin sandstone reservoirs are highlighted, with good lateral continuity and a good well-seismic correspondence, verifying the reliability of this method. Multichannel spectral analysis shows that ( Figure 8The processing results can better expand high-frequency information while retaining low-frequency information. The widening of the frequency band is conducive to the subsequent detailed interpretation of thin reservoirs. This invention can provide a more accurate basis for the exploration and development of thin interbedded oil and gas reservoirs.

[0073] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for processing and identifying thin interbedded layers with high resolution based on compressive sensing with structure constraint, characterized in that, Comprise: S1: by spectrum and waveform feature analysis of seismic data to obtain basic information of data, then using complex cepstrum method to extract wavelet of original seismic record; S2: local structure tensor calculation is carried out to original seismic data, and subsequent eigenvalue decomposition is carried out to local structure tensor, and the eigenvector of each sample point of the data is obtained; S3: according to the seismic wavelet calculated in step S1, the diagonal wavelet measurement matrix in frequency domain and the further perception matrix are constructed; S4: according to the eigenvector of the sample point calculated in step S2, the horizontal and vertical difference operator matrix is combined to establish the constraint matrix; S5: after the multi-channel seismic data is rearranged, the Fourier transform is carried out to obtain the frequency domain data, and the perception matrix of the whole data is calculated according to step S3 and the constraint matrix is obtained according to step S4, then the sparse inversion of reflection coefficient is carried out by using fast threshold shrinkage iteration algorithm; S6: all seismic data are processed in multi-channel cycle, the sparse inversion result of reflection coefficient corresponding to all seismic data is obtained, and it is convolved with a broadband seismic wavelet, and finally the structure constrained high resolution seismic data is obtained.

2. The method of claim 1, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S1, the following steps are further included: according to the frequency band width and phase of the data, the sliding time window is adjusted, and the obtained wavelet frequency band and phase are adjusted.

3. The method of claim 2, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S1, the wavelet of the original seismic record is extracted by using complex cepstrum method, and the specific steps are as follows: first, the time window range of the original seismic data is selected, then the seismic trace in the range is subjected to Fourier transform to obtain its amplitude spectrum, then the amplitude spectrum is subjected to logarithmic transformation and Fourier inverse transformation to obtain the complex cepstrum result, and a low pass filter is designed by setting the width of the pass band to filter the complex cepstrum result, and the Fourier, logarithmic inverse transformation is carried out, and finally the extracted wavelet of the original seismic record is obtained.

4. The method according to any one of claims 1-3, wherein, In step S2, the local structure tensor and the eigenvector are calculated, and the specific steps are as follows: First, the gradient▽g of each sample point in the selected profile data needs to be calculated, as formula (1), then the local structure tensor T is obtained by smoothing the outer product, as formula (2), then the eigenvector u, v is obtained by eigenvalue decomposition of the local structure tensor T of each sample point, as formula (3), wherein is the gradient vector at the sample point, <·> is a smoothing operator, g x is the longitudinal partial derivative, g y is the transverse partial derivative, u, v are the respective eigenvectors.

5. The method according to any one of claims 1-3, wherein, In step S2, the following steps are further included: the estimated gradient is smoothed.

6. The method according to any one of claims 1-3, wherein, In step S3: The seismic wavelet calculated in step S1 is used to construct the diagonal wavelet measurement matrix W in frequency domain, as formula (4), and combined with part of the Fourier basis matrix D, as formula (5), further to obtain the attenuation perception matrix A, A=GW; In the formula, W is the measurement matrix, diag is the diagonal function, FT is the Fourier transform, w(t) is the estimated wavelet, D is the partial Fourier basis matrix, i represents the imaginary number, f is the frequency, and τ is the time.

7. The method according to any one of claims 1-3, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S4: The feature vectors u = (u x ,u y ), v = (v x ,v y ) at the sample points calculated according to step S2 are combined with the established horizontal and vertical difference operator matrices D x , D y to further establish a configuration constraint matrix s, as shown in equation (6); where u, v are eigenvectors, D x is the horizontal difference operator matrix, D y is the vertical difference operator matrix, h u , h v are the eigenvectors u, v constraint parameters, respectively, whose values are between 0 and 1.

8. The method of claim 1-3, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S5: The column rearrangement is completed on the multi-channel seismic data s(t) to obtain As shown in equation (7), Fourier transform is performed on the data y to obtain frequency domain data, and the perceptual matrix A of the data is calculated according to step S3, and the construction constraint matrix S is obtained according to step S4, and the inversion objective function is constructed, as shown in equation (8). where A is the attenuation-aware matrix, is the inverse reflection coefficient, y is the frequency-domain seismic data, λ is the regularization factor, and S is the structural constraint matrix.

9. The method according to any one of claims 1-3, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S5, the sparse inversion of reflection coefficient is carried out by using fast threshold shrinkage iteration algorithm.

10. The method of claim 1-3, wherein the method is based on compressive sensing with structural constraints for thin interbedded high resolution processing and identification. In step S6, the broadband seismic wavelet adopts a higher main frequency Ricker wavelet.

Images