A three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient

Abstract

The application discloses a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient, comprising the following steps: S10, obtaining a three-dimensional seismic data body; S20, constructing a reflection coefficient prior probability model through multi-scale dictionary construction and statistical feature extraction; S30, adaptive matched filtering: obtaining a data body with Q compensation and reflection feature enhancement through time-varying filter and directionality compensation on the original three-dimensional seismic data body and a Q field model; S40, three-dimensional anisotropy enhancement: obtaining enhanced data through structure tensor field calculation and tensor-driven diffusion on the data body in S30; S50, joint inversion optimization: using the reflection coefficient prior probability model and the enhanced data to perform joint inversion to obtain a final high-resolution seismic body. The three-dimensional processing technology fuses physical mechanism and data driving and considers vertical and horizontal resolution improvement, and aims to solve the contradiction between resolution improvement and geological authenticity maintenance in three-dimensional seismic data processing.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of geophysical exploration, in particular to a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficients, which is suitable for high-precision seismic data processing in the fields of oil and gas exploration, mineral resources exploration and the like. BACKGROUND

[0002] With the development of oil and gas exploration to deep layers, complex structures and unconventional reservoirs, the requirement for seismic data resolution is increasingly improved. The effective frequency band range of the traditional seismic data is usually narrow due to the factors such as the earth filter effect, stratum absorption attenuation and acquisition noise, which is difficult to meet the needs of thin interbed identification, small fault detection and reservoir fine description. The current techniques for improving seismic resolution mainly have the following limitations:

[0003] At present, the methods of band expansion such as inverse Q filtering and spectral bluing are usually adopted to widen the frequency band by compensating high-frequency energy, but there are obvious defects: lack of physical constraints: most methods adopt a stable Q value model, which fails to consider the space-time variation characteristics of stratum absorption, resulting in excessive or insufficient compensation of high frequencies in deep layers; noise amplification problem: artificial ringing effect will occur when the signal-to-noise ratio is lower than 2:1; isotropic processing defects: the traditional method adopts a trace-by-trace processing mode in three-dimensional application, ignoring the directional characteristics of wave field propagation.

[0004] There are also inversion methods, and although the inversion method based on sparse constraints can improve the vertical resolution, it has the following problems: reflection coefficient assumption simplification: it is usually assumed that the reflection coefficient obeys Gaussian or Laplace distribution, which is inconsistent with the statistical characteristics of the measured logging data; lack of spatial continuity: the lack of geological structure constraint in two-dimensional / three-dimensional processing leads to the fracture of the same phase axis or false anomalies.

[0005] The deep learning super-resolution technology appeared in recent years faces the following problems: training sample dependence: there is a domain offset problem between synthetic data and field data, and the phenomenon of "pseudo detail" occurs in actual application; poor physical interpretability: the network black box characteristics make it difficult to geologically verify the processing results.

[0006] Therefore, the resolution enhancement modules currently used in this technical field generally have the following problems: contradiction between band expansion and amplitude preservation, lack of three-dimensional anisotropic processing, and limited thin layer identification capability. SUMMARY

[0007] In order to solve the above problems, the present application proposes a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficients, which combines physical mechanism and data-driven, and takes into account the vertical and horizontal resolution improvement of three-dimensional processing technology, aiming to solve the contradiction between resolution improvement and geological authenticity in three-dimensional seismic data processing.

[0008] To achieve the above object, the technical scheme adopted by the present application is: a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient, comprising the steps of:

[0009] S10, obtaining a three-dimensional seismic data body;

[0010] S20, reflection coefficient feature modeling: using three-dimensional seismic data body and logging calibrated reflection coefficient, constructing a reflection coefficient prior probability model through multi-scale dictionary construction and statistical feature extraction;

[0011] S30, adaptive matched filtering: obtaining a data body enhanced in reflection characteristics and compensated by Q through time-varying filter and directionality compensation, by using the original three-dimensional seismic data body and the Q field model;

[0012] S40, three-dimensional anisotropic enhancement: obtaining an enhanced data with geological structure preservation through structure tensor field calculation and tensor-driven diffusion by using the data body obtained in S30;

[0013] S50, joint inversion optimization: using the reflection coefficient prior probability model and the enhanced data obtained in S40 to perform joint inversion to achieve global consistency and obtain a final high-resolution seismic body.

[0014] Further, the multi-scale dictionary construction uses complex wavelet packet decomposition logging reflection coefficient to establish a scale-lithology correlation dictionary.

[0015] Further, in the statistical feature extraction process, the kurtosis of the work area reflection coefficient and the work area autocorrelation length are calculated.

[0016] Further, the kurtosis calculation of the work area reflection coefficient comprises the steps of:

[0017] Data normalization: eliminating the influence of amplitude dimension;

[0018] Sliding window statistics: calculating along the target horizon opening time window, the window length is greater than or equal to λ / 2, and λ is the main wavelength; the three-dimensional work area needs to be divided into block grids;

[0019] Lithology correlation analysis: in the three-dimensional work area block grid, the kurtosis K of different lithologies is respectively calculated;

[0020] K>3: activate strong sparse constraint and use L1 regularization constraint; K<2: use L1 regularization constraint.

[0021] Further, the calculation process of the work area autocorrelation length is:

[0022] Directionality calculation: calculating the variogram along the stratigraphic dip α, strike β and vertical γ respectively;

[0023] Using Kriging interpolation to construct the work area autocorrelation length field.

[0024] Further, the adaptive matched filtering comprises the steps of:

[0025] S301, inputting three-dimensional velocity field, lithology interpretation data, and vertical seismic profile (VSP); estimating Q value: calculating energy attenuation rate for each frequency band of seismic data to estimate Q value of each layer; calibrating Q value in combination with logging acoustic data to ensure consistency with petrophysical properties; three-dimensional modeling: interpolating discrete Q value to continuous three-dimensional field by using geostatistics to obtain time-varying Q field model;

[0026] S302, calculating reflection gradient weight:

[0027] Reflection coefficient gradient extraction: calibrating seismic traces with logging reflection coefficients to generate high-precision reflection coefficient volume; calculating three-dimensional spatial gradient by using noise-resistant Sobel operator to identify lithology mutation area; adaptive weight generation: statistically analyzing gradient amplitude of the entire work area to automatically determine gradient threshold, reducing filtering strength in areas with large gradient and enhancing filtering in areas with small gradient;

[0028] S303, performing directional matched filtering: converting seismic data to frequency-wavenumber domain; rotating coordinate axes according to local stratigraphic dip angle to make filtering direction consistent with strata; applying frequency / Q / gradient triple-weighted filtering; converting back to time-space domain for output;

[0029] S304, feedback optimization: residual analysis: comparing differences between filtered data and logging synthetic records; parameter adjustment: if residual is too large, automatically adjusting Q compensation strength or adjusting gradient weight threshold; updating global parameters once every 5 sections processed.

[0030] Further, in the directional matched filtering, the filter design comprises the following steps:

[0031] Further, the structure tensor field calculation comprises the steps of:

[0032] S411, gradient field calculation: calculating point-by-point gradient of three-dimensional seismic data volume:

[0033] Horizontal direction X / Y: using Sobel operator or Scharr operator to enhance noise resistance;

[0034] Vertical direction Z: using central difference method to preserve thin layer information;

[0035] Outputting three-dimensional gradient vector of each sampling point

[0036] S412, component integration: generate gradient cross product matrix for each point: matrix elements reflect correlation of different directional gradients; average within local window to suppress random noise effects;

[0037] S413, eigen decomposition: perform eigenvalue decomposition on each point's structure tensor matrix to obtain:

[0038] eigenvalue: represents structural strength;

[0039] eigenvector: indicates structural principal direction;

[0040] geological significance mapping includes: homogeneous region; stratigraphic structure; fault / fracture.

[0041] Further, the tensor-driven diffusion includes the steps of:

[0042] S421, diffusion tensor construction: generate diffusion tensor from eigenvalue and eigenvector: primary diffusion direction: set strong diffusion along stratigraphic strike, smooth noise; secondary diffusion direction: moderate diffusion, maintain lateral continuity; weak diffusion direction: suppress diffusion perpendicular to bedding, preserve faults / thin layers; diffusion strength is inversely proportional to eigenvalue, the more obvious the structure, the weaker the diffusion;

[0043] S422, iterative diffusion execution: use explicit Euler method for three-dimensional diffusion: calculate diffusion flux for each point by combining diffusion tensor and local gradient; update data volume according to time step; check eigenvector field change every 3 iterations, dynamically adjust diffusion strength;

[0044] S423, boundary protection mechanism: fault detection, thin layer enhancement;

[0045] S424, termination condition: maximum number of iterations; energy change rate threshold.

[0046] Further, the joint inversion optimization includes the steps of:

[0047] S501, input data preparation: seismic data volume: three-dimensional data after the first three steps; reflectance prior model: from well calibration and statistical characteristics of step S20; structure tensor field: output information from step S30; Q compensation operator: attenuation compensation parameters generated in step S20;

[0048] S502, inversion framework construction:

[0049] variable definition: to-be-solved parameter: high-resolution reflectance model; observation data: preprocessed seismic data volume; forward operator: matrix containing Q compensation and wavefield propagation effects;

[0050] Multi-constraint term design: data fidelity term: force the inversion result to match the seismic data; sparsity constraint term: enhance the pulse of the reflectivity using L1 norm; 3D total variation constraint: keep smooth along the principal direction of the structure tensor, allow sudden change in the perpendicular direction;

[0051] S503, Alternating Direction Optimization: variable splitting: decompose the original problem into two sub-problems, data fitting sub-problem: solve using conjugate gradient method, and constraint term sub-problem: handle sparsity and spatial continuity respectively; iterative steps: (a) fix the constraint term, update the reflectivity model to match the seismic data, (b) fix the reflectivity, apply the sparsity constraint, (c) apply 3D anisotropic smoothing; (d) calculate the residual and adjust the Lagrange multiplier;

[0052] Termination condition: relative error change <0.5% or reach the maximum number of iterations;

[0053] S504, Dynamic parameter adjustment: weight self-adaptation: sparse weight α: dynamically adjust according to the kurtosis of reflectivity; smoothing weight β: according to the autocorrelation length; prior model guidance: preferentially fit the logging calibration section in the early stage of inversion, gradually expand to the whole.

[0054] The beneficial effects of the technical solution are:

[0055] The three-dimensional seismic high-resolution enhancement method based on reflectivity proposed by the application significantly improves the vertical resolution, lateral continuity and geological interpretation reliability of seismic data through three core technologies of adaptive matched filtering, 3D anisotropic diffusion and joint inversion optimization. The technical advantages are explained in detail from three dimensions of quantitative indicators, geological application effect and calculation efficiency.

[0056] Resolution improvement effect: (1) Breakthrough in vertical resolution limit for thin layer identification: traditional method: limited by the wavelength λ of seismic wave, conventional processing is difficult to identify the formation with thickness < λ / 4 (such as λ / 4 ≈ 25m when the main frequency is 30Hz). The application: effectively identifies thin layers of λ / 8-λ / 10 (such as λ / 10 ≈ 6m when the main frequency is increased to 50Hz). Taking the detection results of an oilfield as an example: the identification thickness of sand-shale thin interbedding is increased from 12m to 5m, and the drilling verification rate is increased from 52% to 88%. (2) Bandwidth widening and main frequency improvement.

[0057] Geological interpretation reliability is enhanced: (1) fault and fracture system detection: fault identification accuracy: traditional method: 65%~75% (limited by noise and resolution), the present application: 88%~93% (based on anisotropic diffusion preserving fault boundary); fracture prediction coincidence rate is improved; (2) lithology identification and reservoir description: sand body connectivity analysis: sand body distribution prediction coincidence rate is improved from 68% to 91%; carbonate fracture-cave detection: fracture-cave detection diameter is reduced from > 30m to 15m, and drilling hit rate is increased by 40%.

[0058] Computing efficiency and engineering applicability: (1) processing speed is greatly improved; (2) automation and stability: parameter self-adaptation: dynamic adjustment of parameters such as reflectance kurtosis, Q value and structure tensor, reduces manual intervention; noise resistance: in deep data with signal-to-noise ratio < 2:1, resolution is still improved (traditional method is prone to amplify noise). BRIEF DESCRIPTION OF DRAWINGS

[0059] Figure 1 It is a three-dimensional seismic high-resolution strengthening method based on adaptive reflectance coefficient for the present application.

[0060] Figure 2 It is a flow chart of the adaptive matching filter method in the embodiment of the present application. DETAILED DESCRIPTION

[0061] In order to make the purpose, technical scheme and advantages of the present application clearer, the present application will be further described below with reference to the drawings.

[0062] In this embodiment, referring to the drawing, Figure 1 The present application proposes a three-dimensional seismic high-resolution strengthening method based on adaptive reflectance coefficient, which comprises the following steps:

[0063] S10, obtaining a three-dimensional seismic data body;

[0064] S20, reflectance coefficient feature modeling: using three-dimensional seismic data body and logging calibrated reflectance coefficient, constructing a reflectance coefficient prior probability model through multi-scale dictionary construction and statistical feature extraction;

[0065] S30, adaptive matching filtering: obtaining a data body with Q compensation and reflectance feature enhancement by time-varying filter and directionality compensation on the original three-dimensional seismic data body and Q field model;

[0066] S40, three-dimensional anisotropy enhancement: obtaining enhanced data with geological structure preservation by structure tensor field calculation and tensor driven diffusion on the data body obtained in S30;

[0067] S50, joint inversion optimization: using the reflection coefficient prior probability model and the enhanced data obtained in S40, joint inversion is carried out to achieve global consistency, and finally high-resolution seismic volume is obtained.

[0068] As an optimization scheme of the above embodiment, a complex-valued wavelet packet decomposition logging reflection coefficient is used for multi-scale dictionary construction to establish a scale-lithology correlation dictionary.

[0069] In the statistical feature extraction process, the kurtosis of the work area reflection coefficient and the work area autocorrelation length are calculated to complete feature extraction.

[0070] The kurtosis calculation of the work area reflection coefficient includes the following steps:

[0071] Data normalization: eliminate the influence of amplitude dimension;

[0072] Sliding window statistics: calculate along the target horizon opening time window, the window length is greater than or equal to λ / 2, and λ is the main wavelength; a block grid is established in the three-dimensional work area;

[0073] Lithology correlation analysis: in the three-dimensional work area block grid, the kurtosis K of different lithologies is calculated respectively;

[0074] K>3: activate strong sparse constraint, use L1 regularization constraint; K<2: use L1 regularization constraint.

[0075] The calculation process of the work area autocorrelation length is as follows:

[0076] Directionality calculation: calculate the variogram along the stratigraphic dip α, strike β and vertical γ respectively;

[0077] The autocorrelation length field of the work area is constructed by using Kriging interpolation.

[0078] The present application proposes: directional statistics: propose to directly embed the three-dimensional anisotropic autocorrelation length into the inversion regularization term; the dynamic threshold automatically switches the L1 / L2 constraint according to the kurtosis value. This method reduces the thickness error of thin sand body prediction from ±12m of the traditional method to ±4m through simulation verification and analysis.

[0079] As an optimization scheme of the above embodiment, as shown in Figure 2 The adaptive matched filtering includes the following steps:

[0080] S301, data input three-dimensional velocity field, lithology interpretation data, and vertical seismic profile VSP; Q value estimation: calculate the energy attenuation rate of the seismic data in the frequency band to estimate the Q value of each layer (high frequency attenuation is fast, and the Q value is low); combine the logging acoustic data to calibrate the Q value, and ensure the consistency with the petrophysical properties; three-dimensional modeling: using the method of geological statistics, the discrete Q value is interpolated into a continuous three-dimensional field to obtain a time-varying Q field model;

[0081] S302, reflection gradient weight calculation:

[0082] Reflection coefficient gradient extraction: use logging reflection coefficient to calibrate seismic trace, generate high-precision reflection coefficient volume; adopt anti-noise Sobel operator to calculate three-dimensional space gradient to identify lithology mutation area; adaptive weight generation: statistics of gradient amplitude of the whole work area, automatically determine gradient threshold (such as 1.5 times of median), reduce filter strength in large gradient area (such as fault), enhance filter in small gradient area (homogeneous layer);

[0083] S303, execute directional matched filter: convert seismic data to frequency-wave number domain; rotate coordinate axis according to local stratigraphic dip angle, make filter direction consistent with stratum; apply frequency / Q / gradient triple weighted filter; convert back to time-space domain output;

[0084] Preferably, in the execution of the directional matched filter, the filter design: frequency compensation: dynamically adjust high-frequency enhancement amplitude according to the main frequency and reflection coefficient kurtosis of the target layer (high kurtosis reduces compensation to avoid noise amplification); Q compensation: combined with time-varying Q field, stronger energy recovery is adopted for deep stratum; direction control: rotate filter operator along stratigraphic trend to avoid cross-fault or cross-bedding filtering.

[0085] S304, feedback optimization: residual analysis: compare the difference between filtered data and logging synthetic record; parameter adjustment: if the residual is too large, automatically reduce Q compensation strength or adjust gradient weight threshold, update global parameters once every 5 sections processed.

[0086] In the present application, adaptive matched filtering is proposed as a technology for dynamically adjusting filter parameters, aiming to compensate for energy attenuation (Q absorption effect) during seismic wave propagation, while protecting geological structure characteristics (such as faults and lithological boundaries). The present application proposes reflection coefficient constraint: using logging calibrated reflection coefficient statistical characteristics (such as pulse strength) to guide filter strength; three-dimensional direction perception: dynamically adjust filter direction according to stratigraphic trend / heading; closed-loop feedback optimization: automatically correct parameters through the difference between processing results and logging data.

[0087] The present application proposes: using physical-data dual driving: dynamically combining Q compensation theory and reflection coefficient statistical characteristics; directionally intelligent perception: real-time adjustment of filter direction through structure tensor; closed-loop optimization: introducing residual feedback mechanism in matched filtering; engineering acceleration: GPU-specific memory access optimization for three-dimensional seismic data. According to simulation tests in 7 shale gas blocks, the average thin layer identification capability is improved by 41%. Resolution improvement: effective frequency band is widened from 8-40Hz to 6-60Hz; geological preservation: fault identification accuracy is improved from 68% to 91%; efficiency comparison: processing time is shortened from 18 hours of traditional method to 2.5 hours.

[0088] The essential difference from the conventional method is shown in Table 1.

[0089] Table 1

[0090]

[0091] As an optimization of the above embodiment, in order to quantify the continuity of seismic data in different directions, identify the stratigraphic dip, strike and fault system, the structural tensor field calculation includes the steps of:

[0092] S411, Gradient field calculation: point-by-point gradient calculation on 3D seismic data volume:

[0093] Horizontal direction X / Y: Sobel operator or Scharr operator is used to enhance noise resistance;

[0094] Vertical direction Z: central difference method is used to preserve thin layer information;

[0095] Output the 3D gradient vector of each sampling point

[0096] S412, Component integration: generate a gradient outer product matrix (3x3 symmetric matrix) for each calculation point: matrix elements reflect the correlation of gradients in different directions (such as representing the structural correlation in XY direction); average in a local window (such as 5x5x5 sampling points) to suppress the influence of random noise;

[0097] S413, Eigenvalue decomposition: eigenvalue decomposition is performed on the structural tensor matrix of each point to obtain:

[0098] Eigenvalue: represents the structural strength (eigenvalue λ1≥ λ2≥ λ3, the larger λ1 is, the stronger the continuity is);

[0099] Eigenvector: indicates the main direction of the structure (eigenvector v1, v2, v3, v1 = stratigraphic dip, v3 = vertical bedding direction)

[0100] Geological significance mapping includes: homogeneous area: λ1≈ λ2≈ λ3≈ 0 (no significant directionality); stratified structure: λ1>> λ2≈ λ3 (v1 indicates stratigraphic dip); fault / fracture: λ1≈ λ2>> λ3 (v3 indicates the normal direction of the fracture surface).

[0101] In order to smooth noise along the stratigraphic strike, enhance effective signal in the vertical bedding direction, and preserve fault boundaries, the tensor-driven diffusion includes the steps of:

[0102] S421, diffusion tensor construction: generate diffusion tensor according to eigenvalues and eigenvectors: main diffusion direction: set strong diffusion along the formation trend, smooth noise; secondary diffusion direction: moderate diffusion, maintain lateral continuity; weak diffusion direction: suppress diffusion perpendicular to the bedding direction, preserve the discontinuity / thin layer; diffusion intensity is inversely proportional to eigenvalue, the more obvious the structure, the weaker the diffusion;

[0103] S422, iterative diffusion execution: use explicit Euler method for three-dimensional diffusion: calculate the diffusion flux of each point by combining the diffusion tensor and the local gradient; update the data volume according to the time step; check the change of the eigenvector field every 3 iterations, and dynamically adjust the diffusion intensity;

[0104] S423, boundary protection mechanism: fault detection: when λ1 / λ3>10, it is determined as a fault area, and the diffusion is frozen; thin layer enhancement: local anti-diffusion in the vertical λ3 direction (equivalent to sharpening);

[0105] S424, termination condition: maximum number of iterations (usually 10-15 times); energy change rate threshold (stop when the energy change of the entire data volume is less than 0.1%).

[0106] The above method realizes GPU acceleration: the three-dimensional data is divided into blocks (such as 64x64x64), and each CUDA thread block processes a sub-volume, and the shared memory is used to cache the gradient field, reducing the global memory access delay; parameter self-adaptation; geological consistency check: compare with the logging synthetic record after each iteration, to prevent over-smoothing.

[0107] The essential difference from the traditional method is shown in Table 2.

[0108] Table 2

[0109] Comparison dimension Traditional isotropic diffusion The method of the present patent Directional perception No directional distinction Three-dimensional geological structure driving Fault protection Fuzzy fault boundary Automatic detection and freezing diffusion Thin layer processing Vertical information loss Targeted enhancement of vertical resolution

[0110] The actual application effect is shown in Table 3.

[0111] Table 3

[0112] Index Before processing After processing Vertical resolution 12m 7m Fault recognition rate 65％ 88％ Signal-to-noise ratio improvement 8dB 12dB

[0113] As an optimization scheme of the above embodiment, by integrating seismic data, reflectivity prior knowledge and geological structure constraints, the resolution is improved while ensuring geological rationality, and the following problems are mainly solved: data fitting: make the inversion result match the original seismic data; sparsity constraint: reflect the pulse characteristics of reflectivity; spatial continuity: maintain the lateral extension of the formation, and avoid false anomalies. The joint inversion optimization includes the following steps:

[0114] S501, Input data preparation: Seismic data volume: 3D data after previous three steps; Reflection coefficient prior model: statistical features from well calibration and step S20; Structural tensor field: output information from step S30; Q compensation operator: attenuation compensation parameters generated from step S20;

[0115] S502, Inversion framework construction:

[0116] Variable definition: Parameters to be solved: high-resolution reflection coefficient model; Observation data: preprocessed seismic data volume; Forward operator: matrix containing Q compensation and wavefield propagation effects;

[0117] Multi-constraint term design: Data fidelity term: force the inversion result to match the seismic data; Sparse constraint term: enhance the pulse nature of reflection coefficients using L1 norm; 3D total variation constraint: maintain smoothness along the principal direction of the structural tensor and allow abrupt changes in the perpendicular direction;

[0118] S503, Alternating Direction Optimization (ADMM algorithm): Variable splitting: decompose the original problem into two sub-problems, data fitting sub-problem: solve using the conjugate gradient method, and constraint term sub-problem: handle sparsity and spatial continuity respectively; Iterative steps: (a) fix the constraint term, update the reflection coefficient model to match the seismic data, (b) fix the reflection coefficient, apply the sparse constraint, (c) apply 3D anisotropic smoothing; (d) calculate the residual and adjust the Lagrange multiplier;

[0119] Termination condition: relative error change <0.5% or reach the maximum number of iterations (usually 20-30 times);

[0120] S504, Dynamic parameter adjustment: Weight adaptation: sparse weight α: dynamically adjust according to the reflection coefficient kurtosis (high kurtosis increases α); Smoothing weight β: according to the autocorrelation length (horizontal direction β is larger, vertical direction β is smaller); Prior model guidance: preferentially fit the well calibration section in the early stage of inversion, and gradually expand to the global.

[0121] Actual application effect, as shown in Table 4.

[0122] Table 4

[0123] Index Traditional inversion The method of the present patent Vertical resolution 10m 6m Fault recognition rate 72％ 91％ Thin sand body detection rate 53％ 82％

[0124] The present application proposes: geological-driven adaptive constraints: automatic association of reflection coefficient sparsity and formation lithology (such as sandstone layer α=0.8, mudstone layer α=0.5), real-time guidance of anisotropic smoothing strength by structural tensor; Closed-loop feedback system: feedback the inversion residual to steps S20 Q compensation and S30 diffusion parameters; Engineering acceleration: use mixed precision calculation.

[0125] The technology makes the prediction coincidence rate of 2-5 meter thin sand layer increase from 38% to 76% in simulation experiment test, and the drilling coincidence rate reaches 89%.

[0126] The above shows and describes the basic principles and main features of the present application and the advantages of the present application. Those skilled in the art should understand that the present application is not limited to the above-mentioned embodiments, and the above-mentioned embodiments and descriptions in the specification are only to illustrate the principles of the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficients, characterized in that, The method comprises the steps of: S10, obtaining a three-dimensional seismic data volume; S20, reflection coefficient feature modeling: using the three-dimensional seismic data volume and the logging calibrated reflection coefficient, a reflection coefficient prior probability model is constructed through multi-scale dictionary construction and statistical feature extraction; Multi-scale dictionary construction adopts complex wavelet packet decomposition of logging reflection coefficients to establish a scale-lithology correlation dictionary; In the statistical feature extraction process, the kurtosis of the work area reflection coefficient and the work area autocorrelation length are calculated; The kurtosis calculation of the work area reflection coefficient comprises the steps of: Data normalization: eliminating the influence of amplitude dimension; Sliding window statistics: calculating along the target horizon opening time window, the window length is greater than or equal to λ / 2, and λ is the main wavelength; The three-dimensional work area needs to be divided into block grids; Lithology correlation analysis: in the three-dimensional work area block grid, the kurtosis K of different lithologies is respectively calculated; K>3: activate strong sparse constraint, use L1 regularization constraint; K<2: use L2 regularization constraint; S30, adaptive matched filtering: the original three-dimensional seismic data volume and the Q field model are subjected to time-varying filter and directionality compensation to obtain a data volume subjected to Q compensation and reflection feature enhancement; S40, three-dimensional anisotropy enhancement: the data volume obtained in S30 is subjected to structure tensor field calculation and tensor driven diffusion to obtain enhanced data with geological structure preservation; S50, joint inversion optimization: using the reflection coefficient prior probability model and the enhanced data obtained in S40, joint inversion is carried out to achieve global consistency and obtain a final high-resolution seismic volume.

2. The method of claim 1, wherein, The calculation process of the work area autocorrelation length comprises the steps of: Directionality calculation: calculating the variogram along the stratigraphic dip α, strike β and vertical γ respectively; The work area autocorrelation length field is constructed by using Kriging interpolation.

3. The method of claim 1, wherein the method is based on adaptive reflection coefficients. The adaptive matched filtering comprises the steps of: S301, data input three-dimensional velocity field, lithology interpretation data and vertical seismic profile VSP; Q value estimation: calculating the energy attenuation rate of the seismic data in different frequency bands to estimate the Q value of each layer; combining the logging acoustic data to calibrate the Q value to ensure consistency with the petrophysical properties; three-dimensional modeling: using the method of geological statistics to interpolate discrete Q values into a continuous three-dimensional field to obtain a time-varying Q field model; S302, reflection gradient weight calculation: Reflection coefficient gradient extraction: using logging reflection coefficients to calibrate seismic traces to generate a high-precision reflection coefficient volume; using the anti-noise Sobel operator to calculate the three-dimensional spatial gradient to identify the lithology mutation area; adaptive weight generation: statistically calculating the gradient amplitude of the entire work area to automatically determine the gradient threshold, reducing the filtering strength in the area with large gradient and enhancing the filtering in the area with small gradient; S303, performing directional matched filtering: converting the seismic data into the frequency-wavenumber domain; rotating the coordinate axes according to the local stratigraphic dip to make the filtering direction consistent with the strata; applying frequency / Q / gradient triple weighting filtering; converting back to the time-space domain for output; S304, feedback optimization: residual analysis: comparing the difference between the filtered data and the logging synthetic record; parameter adjustment: if the residual is too large, automatically adjusting the Q compensation strength or adjusting the gradient weight threshold; updating the global parameters once every 5 profiles processed.

4. The method of claim 3, wherein, The filter design in the directional matched filtering: frequency compensation: dynamically adjust the high frequency enhancement amplitude according to the main frequency and the kurtosis of the reflection coefficient of the target layer; Q compensation: combined with the time-varying Q field, stronger energy recovery is adopted for deep strata; Direction control: rotate the filter operator along the strata inclination to avoid cross-fault or cross-bedding filtering.

5. The method of claim 1, wherein, The structural tensor field calculation includes the steps of: S411, gradient field calculation: point-by-point gradient calculation is performed on the three-dimensional seismic data volume: Horizontal direction X / Y: Sobel operator or Scharr operator is used to enhance noise resistance; Vertical direction Z: central difference method is used to preserve thin layer information; Output the three-dimensional gradient vector of each sampling point (∂I / ∂x, ∂I / ∂y, ∂I / ∂z); S412, component integration: generate a gradient outer product matrix for each calculation point: the matrix elements reflect the correlation of gradients in different directions; Average in the local window to suppress the influence of random noise; S413, eigenvalue decomposition: eigenvalue decomposition is performed on the structural tensor matrix of each point to obtain: Eigenvalue: characterizes the structural strength; Eigenvector: indicates the main direction of the structure; Geological significance mapping includes: homogeneous area; stratified structure; fault / crack.

6. The method of claim 1, wherein, The tensor-driven diffusion includes the steps of: S421, diffusion tensor construction: generate a diffusion tensor according to the eigenvalue and eigenvector: main diffusion direction: set strong diffusion along the strata trend to smooth noise; secondary diffusion direction: moderate diffusion to maintain lateral continuity; weak diffusion direction: suppress diffusion perpendicular to the bedding direction to preserve faults / thin layers; the diffusion intensity is inversely proportional to the eigenvalue, and the more obvious the structure, the weaker the diffusion; S422, iterative diffusion execution: three-dimensional diffusion is performed using the explicit Euler method: calculate the diffusion flux of each point by combining the diffusion tensor and the local gradient; update the data volume according to the time step; check the change of the eigenvector field every 3 iterations to dynamically adjust the diffusion intensity; S423, boundary protection mechanism: fault detection, thin layer enhancement; S424, termination condition: maximum number of iterations; Energy change rate threshold.

7. The method of claim 1, wherein the method is a method of adaptive reflectivity based 3D seismic high resolution enhancement. The joint inversion optimization includes the steps of: S501, input data preparation: seismic data volume: three-dimensional data processed by the previous three steps; reflection coefficient prior model: from logging calibration and statistical characteristics of step S20; structural tensor field: output information of step S30; Q compensation operator: attenuation compensation parameters generated in step S20; S502, inversion framework construction: Variable definition: parameters to be solved: high-resolution reflection coefficient model; observation data: preprocessed seismic data volume; forward operator: matrix containing Q compensation and wave field propagation effect; Multi-constraint term design: data fidelity term: force the inversion result to match the seismic data; sparse constraint term: use L1 norm to enhance the pulse nature of the reflection coefficient; Three-dimensional total variation constraint: keep smooth along the main direction of the structural tensor, and allow abrupt change in the vertical direction; S503, Alternating Direction Optimization: Variable Splitting: decompose the original problem into two sub-problems, data fitting sub-problem: solve with conjugate gradient method, and constraint term sub-problem: handle sparsity and spatial continuity separately; Iterative steps: (a) fix the constraint term, update the reflectivity model to match the seismic data, (b) fix the reflectivity, apply sparsity constraint, (c) apply 3D anisotropic smoothing; (d) calculate the residual and adjust the Lagrange multiplier; Termination condition: relative error change <0.5% or reach the maximum number of iterations; S504, Dynamic parameter adjustment: Weight adaptation: sparse weight α: dynamically adjust according to the reflectivity kurtosis; smoothing weight β: according to the autocorrelation length; Prior model guidance: prioritize fitting the logging calibration section at the initial stage of inversion, gradually expand to the global.

Images