Three-dimensional earthquake high-resolution strengthening method based on adaptive reflection coefficient

Abstract

The invention discloses a three-dimensional seismic high-resolution strengthening method based on an adaptive reflection coefficient. The method comprises the following steps: S10, acquiring a three-dimensional seismic data volume; s20, constructing a reflection coefficient prior probability model through multi-scale dictionary construction and statistical feature extraction; s30, adaptive matched filtering: performing time-varying filter and directional compensation on the original three-dimensional seismic data volume and the Q-field model to obtain a data volume subjected to Q compensation and reflection feature enhancement; s40, three-dimensional anisotropy enhancement: carrying out structure tensor field calculation and tensor driving diffusion on the data volume in the step S30 to obtain enhanced data; and S50, joint inversion optimization: performing joint inversion by using the reflection coefficient prior probability model and the enhanced data to obtain a final high-resolution seismic body. According to the three-dimensional seismic data processing method, a physical mechanism and data driving are fused, a three-dimensional processing technology considering vertical and transverse resolution improvement is adopted, and the contradiction between resolution improvement and geological authenticity maintenance in three-dimensional seismic data processing is solved.

Description

Technical Field

[0001] The present invention belongs to the field of geophysical exploration technology, and in particular relates to a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient, which is suitable for high-precision seismic data processing in the fields of oil and gas exploration, mineral resource exploration, etc. Background Art

[0002] As oil and gas exploration gradually expands into deeper layers, complex structures, and unconventional reservoirs, the requirements for seismic data resolution are increasing. Traditional seismic data is limited by factors such as geodetic filtering, formation absorption attenuation, and acquisition noise. The effective frequency band is typically narrow, making it difficult to identify thin interbeds, detect small faults, and provide detailed reservoir descriptions. Current technologies for improving seismic resolution have the following main limitations:

[0003] Currently, the most commonly used frequency band expansion methods include inverse Q filtering and spectrum blueing technologies, which broaden the frequency band by compensating for high-frequency energy. However, they have obvious defects: Insufficient physical constraints: Most methods use a stable Q value model and fail to consider the spatiotemporal variation characteristics of formation absorption, resulting in excessive or insufficient compensation of deep high frequencies; noise amplification problems, which will produce an artificial ringing effect when the signal-to-noise ratio is lower than 2:1; isotropic processing defects: Traditional methods use a channel-by-channel processing mode in three-dimensional applications, ignoring the directional characteristics of wavefield propagation.

[0004] There are also inversion methods. Although the inversion method based on sparse constraints can improve the vertical resolution, it has the following problems: simplified reflection coefficient assumption: 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: lack of geological structure constraints in two-dimensional / three-dimensional processing, resulting in phase axis fracture or false anomalies.

[0005] The deep learning super-resolution technology that has emerged in recent years faces the following challenges: dependence on training samples: there is a domain offset problem between synthetic data and field data, and the phenomenon of "pseudo-details" occurs in actual applications; poor physical interpretability: the black box characteristics of the network make the processing results difficult to verify geologically.

[0006] Therefore, the resolution enhancement modules currently used in this technical field generally have the following problems: the contradiction between frequency band expansion and amplitude preservation, insufficient three-dimensional anisotropy processing, and limited thin layer recognition capabilities. Summary of the Invention

[0007] To address the above issues, the present invention proposes a 3D seismic high-resolution enhancement method based on adaptive reflection coefficients. This method integrates physical mechanisms with data-driven 3D processing technology, taking into account both vertical and lateral resolution enhancement. It aims to resolve the contradiction between resolution enhancement and geological authenticity preservation in 3D seismic data processing.

[0008] To achieve the above object, the present invention adopts a technical solution: a 3D seismic high-resolution enhancement method based on adaptive reflection coefficient, comprising the steps of:

[0009] S10, acquiring a three-dimensional seismic data volume;

[0010] S20, reflection coefficient feature modeling: using 3D seismic data and well logging to calibrate reflection coefficients, a priori probability model of reflection coefficients is constructed through multi-scale dictionary construction and statistical feature extraction;

[0011] S30, adaptive matched filtering: The original 3D seismic data volume and Q field model are subjected to a time-varying filter and directionality compensation to obtain a Q-compensated and reflection-enhanced data volume;

[0012] S40, 3D anisotropic enhancement: The data volume obtained in S30 is subjected to structural tensor field calculation and tensor-driven diffusion to obtain enhanced data with geological structure preservation;

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

[0014] Furthermore, the multi-scale dictionary is constructed by decomposing the logging reflection coefficients using complex-valued wavelet packets to establish a scale-lithology correlation dictionary.

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

[0016] Furthermore, the calculation of the kurtosis of the reflection coefficient of the work area includes the following steps:

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

[0018] Sliding window statistics: open a time window along the target layer, with the window length ≥ λ / 2, where λ is the main wavelength; establish a 3D work area grid;

[0019] Lithologic correlation analysis: In the three-dimensional work area, the kurtosis K of different lithologies is calculated separately;

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

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

[0022] Directional calculation: Calculate the variogram along the formation dip α, strike β, and vertical γ respectively;

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

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

[0025] S301: Data input includes 3D velocity field, lithologic interpretation data, and vertical seismic profiles (VSPs). Q-value estimation is performed: energy attenuation rates are calculated for each frequency band of the seismic data to estimate the Q value of each layer. The Q value is calibrated using well logging acoustic wave data to ensure consistency with rock physical properties. 3D modeling is performed: discrete Q values are interpolated into a continuous 3D field using geostatistical methods to obtain a time-varying Q field model.

[0026] S302, reflection gradient weight calculation:

[0027] Reflection coefficient gradient extraction: Seismic traces are calibrated using well logging reflection coefficients to generate high-precision reflection coefficient volumes. A noise-resistant Sobel operator is used to calculate three-dimensional spatial gradients to identify areas of lithologic abrupt changes. Adaptive weight generation: Calculates the gradient amplitude across the entire work area and automatically determines the gradient threshold, reducing the filter intensity in areas with large gradients and enhancing the filter in areas with small gradients.

[0028] S303, perform directional matched filtering: convert the seismic data to the frequency-wavenumber domain; rotate the coordinate axis according to the local formation dip to make the filtering direction consistent with the formation; apply frequency / Q / gradient triple weighted filtering; convert back to the time-space domain output;

[0029] S304, Feedback Optimization: Residual Analysis: Compare the differences between the filtered data and the synthetic logging records; Parameter Adjustment: If the residual is too large, automatically reduce the Q compensation intensity or adjust the gradient weight threshold. Update the global parameters after processing every 5 profiles.

[0030] Furthermore, the filter design in the directional matched filtering is as follows: frequency compensation: dynamically adjusting the high-frequency enhancement amplitude according to the main frequency and reflection coefficient kurtosis of the target layer; Q compensation: combining the time-varying Q field to adopt stronger energy recovery for deep formations; direction control: rotating the filter operator along the formation dip to avoid cross-fault or cross-stratification filtering.

[0031] Furthermore, the structure tensor field calculation includes the steps of:

[0032] S411, gradient field calculation: perform point-by-point gradient calculation on the 3D seismic data volume:

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

[0034] Vertical direction Z: Use central difference method to retain thin layer information;

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

[0036] S412, construct component integration: Generate a gradient outer product matrix for each calculation point: the matrix elements reflect the correlation of gradients in different directions; average within the local window to suppress the influence of random noise;

[0037] S413, eigendecomposition: Perform eigenvalue decomposition on the structure tensor matrix of each point to obtain:

[0038] Eigenvalue: characterizes structural strength;

[0039] Eigenvector: indicates the main direction of the structure;

[0040] Geologically significant mapping includes: homogeneous areas; layered structures; faults / fractures.

[0041] Furthermore, the tensor driven diffusion comprises the steps of:

[0042] S421, Diffusion Tensor Construction: Generates a diffusion tensor based on eigenvalues and eigenvectors: Primary Diffusion Direction: Strong diffusion is set along the stratigraphic direction to smooth noise; Secondary Diffusion Direction: Moderate diffusion is set to maintain lateral continuity; Weak Diffusion Direction: Suppresses diffusion perpendicular to the bedding direction to protect faults / thin layers; Diffusion strength is inversely proportional to the eigenvalue; the more pronounced the structure, the weaker the diffusion.

[0043] S422, iterative diffusion execution: 3D diffusion is performed using the explicit Euler method: the diffusion flux at each point is calculated by combining the diffusion tensor and the local gradient; the data volume is updated according to the time step; after every 3 iterations, the eigenvector field changes are checked and the diffusion intensity is dynamically adjusted;

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

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

[0046] Furthermore, the joint inversion optimization comprises the steps of:

[0047] S501, input data preparation: seismic data volume: 3D data processed by the first three steps; reflection coefficient prior model: statistical features from well logging calibration and 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: Parameter to be solved: high-resolution reflection coefficient model; Observation data: pre-processed seismic data volume; Forward operator: matrix including Q compensation and wavefield propagation effects;

[0050] Multiple constraint design: Data fidelity: forces the inversion results to match the seismic data; Sparse constraint: uses the L1 norm to enhance the impulse of the reflection coefficient; 3D total variation constraint: maintains smoothness along the main direction of the structural tensor, and allows mutations in the vertical direction;

[0051] S503, alternating direction optimization: variable splitting: decompose the original problem into two subproblems: the data fitting subproblem: solved using the conjugate gradient method, and the constraint term subproblem: dealing with sparsity and spatial continuity respectively; iterative steps: (a) fix the constraint term and update the reflection coefficient model to match the seismic data, (b) fix the reflection coefficient and impose the sparsity constraint, (c) apply 3D anisotropic smoothing; (d) calculate the residual and adjust the Lagrange multiplier;

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

[0053] S504, dynamic parameter adjustment: weight adaptation: sparse weight α: dynamically adjusted according to the kurtosis of the reflection coefficient; smoothing weight β: according to the autocorrelation length; prior model guidance: in the early stage of inversion, the well logging calibration section is preferentially fitted and gradually expanded to the global section.

[0054] The beneficial effects of adopting this technical solution are:

[0055] The reflection coefficient-based 3D seismic high-resolution enhancement method proposed in this paper significantly improves the vertical resolution, lateral continuity, and geological interpretation reliability of seismic data through three core technologies: adaptive matched filtering, 3D anisotropic diffusion, and joint inversion optimization. The following details its technical advantages from the perspectives of quantitative indicators, geological application effects, and computational efficiency.

[0056] Resolution improvement effect: (1) Vertical resolution breaks through the thin layer identification limit: Traditional methods: Limited by the wavelength λ of seismic waves, conventional processing is difficult to distinguish strata with a thickness of less than λ / 4 (such as when the main frequency is 30Hz, λ / 4≈25m). The present invention: effectively identifies thin layers of λ / 8 to λ / 10 (such as when the main frequency is increased to 50Hz, λ / 10≈6m). Taking the detection results of a certain oil field as an example: the identification thickness of thin interbedded sandstone and mudstone layers is increased from 12m to 5m, and the drilling verification rate is increased from 52% to 88%. (2) Band widening and main frequency improvement.

[0057] Enhanced reliability of geological interpretation: (1) Fault and fracture system detection: Fault identification accuracy: Traditional method: 65% to 75% (limited by noise and resolution), present invention: 88% to 93% (based on anisotropic diffusion to protect fault boundaries); Fracture prediction consistency rate is improved; (2) Lithology identification and reservoir description: Sand body connectivity analysis: Sand body distribution prediction consistency rate is increased from 68% to 91%; Carbonate rock fracture and cave detection: The fracture and cave body detection diameter is reduced from >30m to 15m, and the drilling hit rate is increased by 40%.

[0058] Computational efficiency and engineering applicability: (1) Processing speed is greatly improved; (2) Automation and stability: Parameter adaptation: Dynamic adjustment of parameters such as reflection coefficient kurtosis, Q value, and structure tensor to reduce manual intervention; Noise resistance: In deep data with a signal-to-noise ratio of <2:1, the resolution can still be improved (traditional methods tend to amplify noise). BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 This is a flow chart of a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient according to the present invention;

[0060] Figure 2 Flowchart of the adaptive matched filtering method in an embodiment of the present invention. DETAILED DESCRIPTION

[0061] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention is further described below with reference to the accompanying drawings.

[0062] In this embodiment, see Figure 1 As shown, the present invention proposes a three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient, comprising the steps of:

[0063] S10, acquiring a three-dimensional seismic data volume;

[0064] S20, reflection coefficient feature modeling: using 3D seismic data and well logging to calibrate reflection coefficients, a priori probability model of reflection coefficients is constructed through multi-scale dictionary construction and statistical feature extraction;

[0065] S30, adaptive matched filtering: The original 3D seismic data volume and Q field model are subjected to a time-varying filter and directionality compensation to obtain a Q-compensated and reflection-enhanced data volume;

[0066] S40, 3D anisotropic enhancement: The data volume obtained in S30 is subjected to structural tensor field calculation and tensor-driven diffusion to obtain enhanced data with geological structure preservation;

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

[0068] As an optimization solution of the above embodiment, the multi-scale dictionary is constructed by decomposing the logging reflection coefficient using complex-valued wavelet packets to establish a scale-lithology association dictionary.

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

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

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

[0072] Sliding window statistics: open a time window along the target layer, with the window length ≥ λ / 2, where λ is the main wavelength; establish a 3D work area grid;

[0073] Lithologic correlation analysis: In the three-dimensional work area, the kurtosis K of different lithologies is calculated separately;

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

[0075] Among them, the calculation process of the work area autocorrelation length is:

[0076] Directional calculation: Calculate the variogram along the formation dip α, strike β, and vertical γ respectively;

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

[0078] This paper proposes: directional statistics: directly embedding the three-dimensional anisotropic autocorrelation length into the inversion regularization term; and dynamic thresholding that automatically switches between L1 and L2 constraints based on the kurtosis value. This method, through simulation validation and analysis, has reduced the thickness error of thin sand bodies from ±12m with traditional methods to ±4m.

[0079] As an optimization solution of the above embodiment, Figure 2 As shown, the adaptive matched filtering includes the steps of:

[0080] S301: Data input includes a 3D velocity field, lithologic interpretation data, and vertical seismic profiles (VSPs). Q-value estimation is performed: the energy attenuation rate is calculated for each frequency band of the seismic data, estimating the Q value of each layer (higher frequencies have faster attenuation, resulting in lower Q values). The Q value is calibrated using well logging acoustic wave data to ensure consistency with rock physical properties. 3D modeling is performed: geostatistical methods are used to interpolate discrete Q values into a continuous 3D field, resulting in a time-varying Q field model.

[0081] S302, reflection gradient weight calculation:

[0082] Reflection coefficient gradient extraction: Seismic traces are calibrated using well logging reflection coefficients to generate high-precision reflection coefficient volumes. A noise-resistant Sobel operator is used to calculate three-dimensional spatial gradients to identify areas of lithologic abrupt changes. Adaptive weight generation: Calculates the gradient amplitude across the entire work area and automatically determines the gradient threshold (e.g., 1.5 times the median value). This reduces the filtering intensity in areas with large gradients (e.g., faults) and enhances the filtering in areas with small gradients (e.g., homogeneous layers).

[0083] S303, perform directional matched filtering: convert the seismic data to the frequency-wavenumber domain; rotate the coordinate axis according to the local formation dip to make the filtering direction consistent with the formation; apply frequency / Q / gradient triple weighted filtering; convert back to the time-space domain output;

[0084] Preferably, the filter design in the directional matching filter is as follows: frequency compensation: dynamically adjust the high-frequency enhancement amplitude according to the main frequency and reflection coefficient kurtosis of the target layer (the compensation is reduced if the kurtosis is high to avoid noise amplification); Q compensation: combine the time-varying Q field to adopt stronger energy recovery for deep formations; direction control: rotate the filter operator along the formation dip to avoid cross-fault or cross-stratification filtering.

[0085] S304, Feedback Optimization: Residual Analysis: Compare the differences between the filtered data and the synthetic logging records; Parameter Adjustment: If the residual is too large, automatically reduce the Q compensation intensity or adjust the gradient weight threshold. Update the global parameters after processing every 5 profiles.

[0086] The adaptive matched filter proposed in this paper is a technique for dynamically adjusting filter parameters to compensate for energy attenuation (Q absorption effect) during seismic wave propagation while protecting geological structural features (such as faults and lithologic boundaries). This paper proposes reflection coefficient constraints: using the statistical characteristics of the reflection coefficient calibrated by well logging (such as the strength of impulses) to guide the filter intensity; three-dimensional directional sensing: dynamically adjusting the filter direction based on the dip / strike of the formation; and closed-loop feedback optimization: automatically correcting parameters based on the difference between the processing results and the well logging data.

[0087] This invention proposes: a dual-driven approach using physics and data: dynamically combining Q compensation theory with the statistical characteristics of reflection coefficients; intelligent directional sensing: real-time adjustment of filter direction through structural tensors; closed-loop optimization: introducing a residual feedback mechanism in matched filtering; and engineered acceleration: optimizing GPU-specific memory access for 3D seismic data. Simulation tests show that this solution improved thin-bed identification by an average of 41% in seven shale gas blocks. Resolution was improved: the effective frequency band was widened from 8-40Hz to 6-60Hz; geological preservation: fault identification accuracy increased from 68% to 91%. Compared to conventional methods, processing time was reduced from 18 hours to 2.5 hours.

[0088] The essential differences from the traditional method are shown in Table 1.

[0089] Table 1

[0090]

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

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

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

[0094] Vertical direction Z: Use central difference method to retain thin layer information;

[0095] Output the three-dimensional gradient vector of each sampling point

[0096] S412, construct component integration: generate the gradient outer product matrix (3×3 symmetric matrix) for each calculation point: the matrix elements reflect the correlation of gradients in different directions (such as Indicates XY structural correlation); average within a local window (e.g., 5×5×5 sampling points) to suppress the influence of random noise;

[0097] S413, eigendecomposition: Perform eigenvalue decomposition on the structure tensor matrix of each point to obtain:

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

[0099] Eigenvectors: Indicates the main direction of the structure (eigenvectors v1, v2, v3, v1 = stratum dip, v3 = vertical bedding direction)

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

[0101] In order to smooth the noise along the formation strike, enhance the effective signal perpendicular to the bedding direction, and protect the fault boundary, the tensor-driven diffusion includes the following steps:

[0102] S421, Diffusion Tensor Construction: Generates a diffusion tensor based on eigenvalues and eigenvectors: Primary Diffusion Direction: Strong diffusion is set along the stratigraphic direction to smooth noise; Secondary Diffusion Direction: Moderate diffusion is set to maintain lateral continuity; Weak Diffusion Direction: Suppresses diffusion perpendicular to the bedding direction to protect faults / thin layers; Diffusion strength is inversely proportional to the eigenvalue; the more pronounced the structure, the weaker the diffusion.

[0103] S422, iterative diffusion execution: 3D diffusion is performed using the explicit Euler method: the diffusion flux at each point is calculated by combining the diffusion tensor and the local gradient; the data volume is updated according to the time step; after every 3 iterations, the eigenvector field changes are checked and the diffusion intensity is dynamically adjusted;

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

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

[0106] The method proposed in the present invention achieves the following: GPU acceleration: dividing the three-dimensional data into blocks (such as 64×64×64), with each CUDA thread block processing one sub-volume, and using shared memory to cache the gradient field to reduce global memory access latency; parameter adaptation; and geological consistency checking: comparing with the well logging synthetic record after each iteration to prevent over-smoothing.

[0107] The essential differences from the traditional method are shown in Table 2.

[0108] Table 2

[0109] Comparison Dimension Traditional isotropic diffusion This patented method Direction perception No direction distinction 3D geological structure drive Fault protection Fuzzy fault boundaries Automatically detect and freeze diffusion Thin layer treatment Vertical information loss Targeted enhancement of vertical resolution

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

[0111] Table 3

[0112] index Before treatment After processing Vertical resolution 12m 7m Fault recognition rate 65％ 88％ Improved signal-to-noise ratio 8dB 12dB

[0113] As an optimization scheme for the above embodiment, by integrating seismic data, prior knowledge of reflection coefficients, and geological structural constraints, the resolution is improved while ensuring geological rationality. The main problems to be solved are: data fitting: matching the inversion results with the original seismic data; sparsity constraints: reflecting the pulse characteristics of the reflection coefficient; spatial continuity: maintaining the lateral ductility of the formation and avoiding false anomalies. The joint inversion optimization includes the following steps:

[0114] S501, input data preparation: seismic data volume: 3D data processed by the first three steps; reflection coefficient prior model: statistical features from well logging calibration and step S20; structure tensor field: output information from step S30; Q compensation operator: attenuation compensation parameters generated in step S20;

[0115] S502, inversion framework construction:

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

[0117] Multiple constraint design: Data fidelity: forces the inversion results to match the seismic data; Sparse constraint: uses the L1 norm to enhance the impulse of the reflection coefficient; 3D total variation constraint: maintains smoothness along the main direction of the structural tensor, and allows mutations in the vertical direction;

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

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

[0120] S504, dynamic parameter adjustment: weight adaptation: sparse weight α: dynamically adjusted according to the kurtosis of the reflection coefficient (high kurtosis increases α); smoothing weight β: based on the autocorrelation length (horizontally β is larger, vertically β is smaller); prior model guidance: in the early stage of inversion, the well logging calibration segment is preferentially fitted and gradually expanded to the global segment.

[0121] The actual application effect is shown in Table 4.

[0122] Table 4

[0123] index Traditional inversion This patented method Vertical resolution 10m 6m Fault recognition rate 72％ 91％ Thin sand body detection rate 53％ 82％

[0124] The present invention proposes: geologically driven adaptive constraints: automatic association of reflection coefficient sparsity with stratum lithology (e.g., α = 0.8 for sandstone layers and α = 0.5 for mudstone layers), and real-time guidance of anisotropic smoothing intensity by the structural tensor; a closed-loop feedback system: feedback of inversion residuals to step S20Q compensation and S30 diffusion parameters; and engineered acceleration: the use of mixed-precision calculations.

[0125] In simulation experiments, this technology increased the prediction consistency rate of 2-5 meter thin sand layers from 38% to 76%, and the drilling consistency rate reached 89%.

[0126] The basic principles, main features, and advantages of the present invention are shown and described above. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely illustrative of the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and modifications are intended to fall within the scope of the present invention. The scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. A 3D seismic high-resolution enhancement method based on adaptive reflection coefficient, characterized in that: Including steps: S10, acquiring a three-dimensional seismic data volume; S20, reflection coefficient feature modeling: using 3D seismic data and well logging to calibrate reflection coefficients, a priori probability model of reflection coefficients is constructed through multi-scale dictionary construction and statistical feature extraction; S30, adaptive matched filtering: The original 3D seismic data volume and Q field model are subjected to a time-varying filter and directionality compensation to obtain a Q-compensated and reflection-enhanced data volume; S40, 3D anisotropic enhancement: The data volume obtained in S30 is subjected to structural tensor field calculation and tensor-driven diffusion to obtain enhanced data with geological structure preservation; S50, joint inversion optimization: Utilize the reflection coefficient prior probability model and the enhanced data obtained in S40 to perform joint inversion, achieve global consistency, and obtain the final high-resolution seismic volume.

2. The three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient according to claim 1, characterized in that: The multi-scale dictionary is constructed by decomposing the logging reflection coefficients using complex-valued wavelet packets to establish a scale-lithology correlation dictionary.

3. The three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient according to claim 1, characterized in that: During the statistical feature extraction process, the kurtosis of the work area reflection coefficient and the autocorrelation length of the work area are calculated.

4. The method for 3D seismic high-resolution enhancement based on adaptive reflection coefficient according to claim 3, characterized in that: The calculation of the kurtosis of the reflection coefficient of the work area includes the following steps: Data normalization: eliminate the influence of amplitude dimension; Sliding window statistics: Calculate the time window along the target layer, the window length ≥ λ / 2, λ is the main wavelength; To establish a three-dimensional work area, a block grid is required; Lithologic correlation analysis: In the three-dimensional work area, the kurtosis K of different lithologies is calculated separately; K>3: activate strong sparsity constraint and use L1 regularization constraint; K<2: use L1 regularization constraint.

5. The method for 3D seismic high-resolution enhancement based on adaptive reflection coefficient according to claim 3, characterized in that: The calculation process of the autocorrelation length of the work area: Directional calculation: Calculate the variogram along the formation dip α, strike β, and vertical γ respectively; The autocorrelation length field of the work area is constructed using Kriging interpolation.

6. The 3D seismic high-resolution enhancement method based on adaptive reflection coefficient according to claim 1, characterized in that: The adaptive matched filtering comprises the steps of: S301, data input 3D velocity field, lithologic interpretation data, vertical seismic profile VSP; Q-value estimation: Calculate the energy attenuation rate of seismic data in different frequency bands to estimate the Q value of each layer; calibrate the Q value using well logging acoustic wave data to ensure consistency with rock physical properties; 3D modeling: Use geostatistical methods to interpolate discrete Q values into a continuous 3D field to obtain a time-varying Q field model; S302, reflection gradient weight calculation: Reflection coefficient gradient extraction: Seismic traces are calibrated using well logging reflection coefficients to generate high-precision reflection coefficient volumes. A noise-resistant Sobel operator is used to calculate three-dimensional spatial gradients to identify areas of lithologic abrupt changes. Adaptive weight generation: Calculates the gradient amplitude across the entire work area and automatically determines the gradient threshold, reducing the filter intensity in areas with large gradients and enhancing the filter in areas with small gradients. S303, perform directional matched filtering: convert the seismic data to the frequency-wavenumber domain; rotate the coordinate axis according to the local formation dip to make the filtering direction consistent with the formation; apply frequency / Q / gradient triple weighted filtering; convert back to the time-space domain output; S304, Feedback Optimization: Residual Analysis: Compare the differences between the filtered data and the synthetic logging records; Parameter Adjustment: If the residual is too large, automatically reduce the Q compensation intensity or adjust the gradient weight threshold. Update the global parameters after processing every 5 profiles.

7. The method for 3D seismic high-resolution enhancement based on adaptive reflection coefficient according to claim 6, characterized in that: The filter design in the directional matched filtering is as follows: frequency compensation: dynamically adjusting the high-frequency enhancement amplitude according to the main frequency and reflection coefficient kurtosis of the target layer; Q compensation: combining the time-varying Q field to adopt stronger energy recovery for deep formations; direction control: rotating the filter operator along the formation dip to avoid cross-fault or cross-bedding filtering.

8. The three-dimensional seismic high-resolution enhancement method based on adaptive reflection coefficient according to claim 1, characterized in that: The structure tensor field calculation includes the steps of: S411, gradient field calculation: perform point-by-point gradient calculation on the 3D seismic data volume: Horizontal X / Y direction: Use Sobel operator or Scharr operator to enhance noise resistance; Vertical direction Z: Use central difference method to retain thin layer information; Output the three-dimensional gradient vector of each sampling point S412, construct component integration: generate the gradient outer product matrix for each calculation point: the matrix elements reflect the correlation of gradients in different directions; Averaging within a local window to suppress the influence of random noise; S413, eigendecomposition: Perform eigenvalue decomposition on the structure tensor matrix of each point to obtain: Eigenvalue: characterizes structural strength; Eigenvector: indicates the main direction of the structure; Geologically significant mapping includes: homogeneous areas; layered structures; faults / fractures.

9. The method for 3D seismic high-resolution enhancement based on adaptive reflection coefficient according to claim 1, characterized in that: The tensor driven diffusion comprises the steps of: S421, Diffusion Tensor Construction: Generates a diffusion tensor based on eigenvalues and eigenvectors: Primary Diffusion Direction: Strong diffusion is set along the stratigraphic direction to smooth noise; Secondary Diffusion Direction: Moderate diffusion is set to maintain lateral continuity; Weak Diffusion Direction: Suppresses diffusion perpendicular to the bedding direction to protect faults / thin layers; Diffusion strength is inversely proportional to the eigenvalue; the more pronounced the structure, the weaker the diffusion. S422, iterative diffusion execution: 3D diffusion is performed using the explicit Euler method: the diffusion flux at each point is calculated by combining the diffusion tensor and the local gradient; the data volume is updated according to the time step; after every 3 iterations, the eigenvector field changes are checked and the diffusion intensity is dynamically adjusted; S423, boundary protection mechanism: fault detection, thin layer enhancement; S424, termination condition: maximum number of iterations; Energy change rate threshold.

10. The method for 3D seismic high-resolution enhancement based on adaptive reflection coefficient according to claim 1, characterized in that: The joint inversion optimization comprises the steps of: S501, input data preparation: seismic data volume: 3D data processed by the first three steps; reflection coefficient prior model: statistical features from well logging calibration and step S20; structure tensor field: output information from step S30; Q compensation operator: attenuation compensation parameters generated in step S20; S502, inversion framework construction: Variable definition: Parameter to be solved: high-resolution reflection coefficient model; Observation data: pre-processed seismic data volume; Forward operator: matrix including Q compensation and wavefield propagation effects; Multi-constraint design: Data fidelity: forcing the inversion results to match the seismic data; Sparse constraint: using the L1 norm to enhance the impulse nature of the reflection coefficient; Three-dimensional total variation constraint: maintain smoothness along the main direction of the structure tensor, and allow mutations in the vertical direction; S503, alternating direction optimization: variable splitting: decompose the original problem into two subproblems: the data fitting subproblem: solved using the conjugate gradient method, and the constraint term subproblem: dealing with sparsity and spatial continuity respectively; iterative steps: (a) fix the constraint term and update the reflection coefficient model to match the seismic data, (b) fix the reflection coefficient and impose the sparsity constraint, (c) apply 3D anisotropic smoothing; (d) calculate the residual and adjust the Lagrange multiplier; Termination conditions: relative error change < 0.5% or reaching the maximum number of iterations; S504, dynamic parameter adjustment: weight adaptation: sparse weight α: dynamically adjusted according to the kurtosis of the reflection coefficient; smoothing weight β: according to the autocorrelation length; prior model guidance: in the early stage of inversion, the well logging calibration section is preferentially fitted and gradually expanded to the global section.

Images