A reservoir prediction method based on multi-angle seismic difference enhancement

Abstract

The present application belongs to the technical field of geophysical exploration, and particularly relates to a reservoir prediction method based on multi-angle seismic difference enhancement. The method collects data information of a study area, including regional geological data, 3D seismic data, well logging curves, analysis and test data and test data, evaluates data quality, and performs pre-stack trace optimization processing, angle-dependent superposition data volume generation, amplitude-frequency consistency correction, reservoir sensitivity factor calculation and inversion processing, and finally obtains fluid detection and reserve evaluation of the study area. The method significantly improves the thin layer identification capability, effectively eliminates non-geological interference, realizes effective separation of lithology and fluid, and has efficient processing flow, strong industrial application, can improve the exploration success rate and economic benefit, provides a reliable technical means for solving the complex thin interbedded reservoir prediction problem, and has high popularization and application value.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration technology, specifically relating to a reservoir prediction method based on multi-angle seismic differential enhancement. Background Technology

[0002] Seismic reservoir prediction is the process of quantitatively or qualitatively predicting the spatial distribution, porosity, permeability, and fluid properties of reservoirs using subsurface geological information obtained through seismic exploration, combined with multi-source data from geology, well logging, and core samples. As exploration progresses, research on risk zone expansion, zonal evaluation, and risk well placement faces new technical challenges. In particular, the evaluation and selection of favorable zones for formations not encountered during drilling has become crucial for future breakthroughs in risk management, especially in predicting reservoir distribution in such formations using geophysical techniques.

[0003] There are three main traditional reservoir prediction methods: conventional AVO attribute analysis, elastic wave impedance inversion, and spectral decomposition. Each method has its own advantages but also significant drawbacks. Conventional AVO attribute analysis is limited by thin-layer tuning effects; when the reservoir thickness is less than 5m, the AVO characteristics are severely distorted as the amplitude changes with angle. Elastic wave impedance inversion is highly sensitive to wavelets, and the low signal-to-noise ratio of deep seismic data leads to a "layered" illusion in the inversion results. Practical applications show that in strata deeper than 3000m, the relative error in porosity prediction exceeds 30%. Spectral decomposition cannot distinguish between lithological variations and fluid effects; the energy difference between water-bearing sandstone and oil-bearing sandstone in the low-frequency range (<15Hz) is less than 2dB.

[0004] Therefore, there is an urgent need for a seismic reservoir prediction method that can solve some of the above problems to a certain extent. Summary of the Invention

[0005] To achieve the above objectives, this invention provides a reservoir prediction method based on multi-angle seismic differential enhancement. For beach-bar sandstone reservoirs with a thickness less than λ / 8 (λ being the seismic wavelength) in continental basins, the method significantly improves reservoir identification accuracy and fluid detection capability through an optimized five-angle overlay data volume and differential enhancement algorithm. The technical solution is as follows:

[0006] S1. Data Preparation and Quality Control

[0007] Collect data on the study area, including regional geological data, 3D seismic data, well logging curves, analytical and laboratory data, and test data, and evaluate the quality of the data.

[0008] The quality of seismic data was systematically evaluated using quantitative indicators. Well logging data was systematically standardized using histogram and trend surface analysis methods. The time-depth relationship was established using a joint calibration method of VSP data and synthetic records. High-precision synthetic records were created, and multi-scale correlation algorithms were applied to improve calibration accuracy. Rock physics analysis was used to establish a quantitative interpretation scale for the work area, providing a reliable physical basis for subsequent seismic inversion and reservoir prediction.

[0009] S2, Pre-Stack Gather Optimization Processing

[0010] The input raw pre-stack seismic gathers (usually common center point gathers, CMPs, or common reflection point gathers, CRPs) undergo refined preprocessing to eliminate or reduce waveform distortion and energy distortion caused by non-geological factors such as seismic acquisition, near-surface conditions, and wave propagation effects. Preprocessing mainly includes amplitude compensation, multiple suppression, and anisotropy correction.

[0011] Amplitude compensation: Spherical diffusion compensation and surface uniform amplitude compensation are implemented to correct the energy attenuation caused by geometric diffusion during the propagation of seismic waves and the amplitude changes caused by differences in the location of the source and receiver and uneven absorption near the surface. This restores the true relative amplitude of seismic waves propagating in the underground medium and lays a reliable foundation for subsequent AVO / AVA analysis.

[0012] Multiple wave suppression: Advanced multichannel filtering techniques (such as high-precision Radon transform, parabolic Radon transform, or surface-correlated multiple wave attenuation SRME) are applied to effectively identify and suppress interlayer multiple waves, reverberation, and other interfering wave fields, significantly improving the signal-to-noise ratio and fidelity of the primary reflected wave, and ensuring the accuracy of subsequent angle superposition and attribute extraction.

[0013] Anisotropy correction: For azimuth anisotropy in the overlying strata (such as VTI medium) or fracture-induced anisotropy, azimuth time difference correction or anisotropy parameter inversion and correction are performed to eliminate the systematic deviations caused by anisotropy on the travel time and amplitude of far-path (large incident angle) earthquakes, and to ensure the consistency of time shift and waveform of gather data at different incident angles.

[0014] S3, Generation of angle-separated superimposed data volumes

[0015] Based on the pre-stack gathers optimized by S2, and taking into account the significant influence of the incident angle of seismic waves on the response characteristics of underground reflection points (especially the amplitude variation with the incident angle, i.e., AVA / AVO characteristics), each trace in the gather is precisely divided according to its corresponding incident angle range.

[0016] The incident angle was divided into five continuous and partially overlapping sub-ranges: 3-13° (central angle 7°), 10-20° (central angle 15°), 17-27° (central angle 22°), 24-34° (central angle 29°), and 31-41° (central angle 36°). For each sub-range of incident angle gathers, precise dynamic correction (NMO / DMO), muting, and stacking were performed, ultimately generating five sets of partially stacked seismic data volumes representing the reflection characteristics of different average incident angles (i.e., central angles). These data volumes retain the dominant reflection information within their respective angle ranges, while effectively suppressing random noise through stacking.

[0017] The angle range design strictly follows the AVO / AVA theory and the physical response mechanism of the target reservoir rock, specifically meeting the following physical conditions:

[0018] Near angle range 3-13° (central angle 7°): This angle range covers the main sensitive area of ​​the AVO intercept response. Within this small incident angle range, the reflected amplitude is mainly controlled by the difference in longitudinal wave velocity or acoustic impedance between the media above and below the interface, and has the best ability to resolve lithological changes, while being relatively insensitive to changes in Poisson's ratio caused by fluids. Therefore, it can effectively reflect background lithological information.

[0019] The far-angle range (31-41°, central angle 36°) lies within a zone of significant sensitivity to Poisson's ratio differences. As the incident angle increases, the sensitivity of the reflected amplitude to changes in shear wave velocity or Poisson's ratio increases dramatically. For hydrocarbon reservoirs, especially IIA-type AVO sandstones, hydrocarbon substitution of pore fluids leads to a significant decrease in Poisson's ratio. The far-angle data volume (X2) amplifies this fluid effect to the greatest extent, resulting in a pronounced amplitude anomaly (negative anomaly or rapid amplitude decay with angle).

[0020] A central angle of 22° (range 17-27°) was chosen as the benchmark: this angle lies in the critical region of AVO amplitude reversal in Class IIA (25°±3°). For typical oil and gas-bearing sandstone (low impedance, low Poisson's ratio), its AVO response typically undergoes polarity reversal near the critical angle (changing from positive or weakly negative amplitude to strongly negative amplitude). The 22° central angle data volume was chosen as the spectral matching benchmark not only because it is in the reversal transition zone and is indicative of fluid changes, but also because it typically has the best signal-to-noise ratio and waveform stability, providing a reliable reference for multi-angle data consistency correction.

[0021] S4, Amplitude-Frequency Consistency Correction

[0022] Due to factors such as seismic wave propagation paths, absorption attenuation, and the inherent characteristics of different stacking angle ranges, there may be systematic differences in amplitude energy and dominant frequency drift among the five sets of angle-separated stacked data volumes generated in step S3. To eliminate these differences caused by non-geological factors and ensure the comparability and physical significance of subsequent difference calculations based on multi-angle data, [further measures are needed].

[0023] A data volume with a central angle of 22° (range 17-27°) was selected as the reference volume because it typically has a good signal-to-noise ratio and a suitable incident angle. Spectral Matching algorithms (such as least-squares spectral shaping filtering and time-varying spectral whitening) were employed, using the amplitude and phase spectra of the reference data volume as templates, to perform channel-by-channel and time-window-by-time amplitude energy adjustment and dominant frequency feature alignment processing on the other four data volumes (central angles of 7°, 15°, 29°, and 36°). This processing strictly preserved the original geological reflection structure within each data volume, correcting only the spectral differences between it and the reference volume, ultimately resulting in five sets of angle-separated stacked data volumes with high consistency in amplitude energy level and dominant frequency characteristics.

[0024] The specific implementation of the spectral matching algorithm includes the following refined operational steps:

[0025] ① Time difference correction: Accurately calculate the data volume of each sub-angle ( θ =7°, 15°, 29°, 36°) relative to the reference data volume ( θ =22°) Two-way travel time system offset at the same reflective horizon Δt ( θ )= t (twenty two)- t ( θ This migration is primarily caused by anisotropic velocity differences in the overlying strata. It is determined using cross-correlation time-shift scanning or dynamic time warping (DTW) algorithms. Δt ( θ It also performs global or time-varying time-shift correction on non-reference data volumes to ensure that all angle data volumes are strictly aligned on the time axis.

[0026] ② Spectral equalization: In the frequency domain, the amplitude spectrum of the data volume at each angle is adjusted by designing and applying a deconvolution shaping operator to approximate the amplitude spectrum characteristics of the reference data volume. Amplitude spectrum adjustment is achieved through the deconvolution operator: (1)

[0027] In the formula, H(f,t) is the time-varying shaping operator; S b (f,t) represents the time-varying amplitude spectrum of the reference data volume (D22); S t(f,t) represents the time-varying amplitude spectrum of the target data volume. This operator compensates for energy attenuation or high-frequency loss in the amplitude spectrum of the target data volume through deconvolution operation, while preserving its phase information.

[0028] ③ Amplitude equalization: After amplitude spectrum equalization, the phase characteristics of all angular data volumes are uniformly corrected using minimum phase transformation technology. By converting mixed-phase seismic data into zero phase or uniform minimum phase, phase distortion caused by differences in the propagation paths of waves at different angles is eliminated, ensuring that the five sets of data volumes are highly consistent in waveform phase characteristics, providing an accurate waveform basis for subsequent multi-angle difference calculations.

[0029] S5, Calculation of reservoir sensitivity factors

[0030] For the petrophysical properties of target reservoirs (such as oil and gas-bearing sandstones) (e.g., Poisson's ratio variations and fluid sensitivity), characteristic near-angle stacked data volume X1 (corresponding to a central angle of 7°) and far-angle stacked data volume X2 (corresponding to a central angle of 36°) are extracted after consistency correction. The near-angle data volume mainly reflects lithological variations and is relatively insensitive to fluids; the far-angle data volume is more sensitive to the decrease in Poisson's ratio caused by reservoir pore fluids, especially oil and gas. Utilizing the difference enhancement characteristics between these two types of data volumes, according to the formula Y=X1- k X2 calculations generate a sensitive factor data volume Y that highlights the hydrocarbon response of the reservoir, where key coefficients... k Determining the value of is crucial, as it is not a fixed constant but is obtained through regression analysis using elastic impedance (EI) information at known well points in the target area.

[0031] The specific method is as follows: At the wellbore location, extract the near-angle elastic impedance EI1 and the far-angle elastic impedance EI2 respectively, and establish the statistical relationship between EI1 and EI2, which is usually a linear regression. The slope of this regression relationship is the coefficient. k . k The values ​​reflect the AVO response characteristics of the target reservoir at a specific angle. The Y data volume generated in this step significantly suppresses the background lithological influence and highlights fluid-related anomalous signals. Coefficients k The method for determining the target layer is a statistical optimization process based on the physical properties of the rock. The specific steps are as follows:

[0032] ① Wellside elastic impedance extraction: At a known well location, for the target reservoir section (usually oil-bearing sandstone and its surrounding rock), the corresponding elastic impedance curves EI7 and EI are accurately extracted from the corrected near-angle (7°) and far-angle (36°) stacked data volumes, respectively. 36 The extraction process needs to be combined with well-seismic calibration to ensure the accuracy of the time-depth relationship.

[0033] ② Sandstone Sensitive Impedance Modeling: Based on rock physical analysis of the target area, a multiple regression equation that can maximize the distinction between sandstone and mudstone is constructed:

[0034] (2)

[0035] In the formula, EI opt The sensitive elastic impedance of the target lithology body obtained from the multi-angle differential body; EI 7 represents the elastic resistance at a center angle of 7°; EI 36 The elastic resistance is at a center angle of 36°; a , b This is an empirical coefficient.

[0036] ③ Coefficient Optimization determination: The sensitivity factor formula Y=X1- k X2 is analogous to the impedance domain, let EI=EI 7- kEI 36 Through systematic change k Values ​​(typically scanned in steps of 0.01 within the range of 2.0 to 2.3) are used to calculate different... Value EI Curve. The evaluation metric is: within the target segment, EI The goal is to maximize the mean difference (distinction) between pure sandstone and pure mudstone sections, or to use the Fisher discrimination criterion. Ultimately, the optimal method that maximizes the sandstone / mudstone distinction is selected. The value has a clear rock physics significance, reflecting the AVO gradient characteristics of the reservoir under specific angular combinations.

[0037] S6, Inversion Processing

[0038] The sensitive factor data volume Y obtained in step S4 is used as the core input data, and constrained inversion is performed by combining well point lithology data, sedimentary mode constraints, and multi-source information on rock physical relationships.

[0039] Among them, well point lithology data: using accurate lithology (sand / mud), porosity, saturation and other curves obtained from well logging interpretation as hard data points and low-frequency model constraints.

[0040] Sedimentary model constraints: Integrating geological knowledge (sedimentary facies distribution maps, sequence stratigraphy framework, etc.) as spatial trend constraints to guide inversion results to conform to geological laws.

[0041] Rock-physical relationship: Using rock-physical models (Gassmann equation, Xu-White model, etc.) established in the target area, a quantitative conversion relationship between the sensitive factor Y and the target reservoir parameters (porosity, clay content, hydrocarbon saturation, etc.) is established.

[0042] The inversion algorithm is used to generate high-precision three-dimensional prediction volumes of reservoir parameters, such as porosity volume, hydrocarbon saturation volume, and effective thickness volume, so as to realize the quantitative description of reservoir spatial distribution and physical properties.

[0043] The inversion process employs a sequential Gaussian simulation algorithm within a geostatistical framework, whose core objective function is defined as:

[0044] (3)

[0045] In the formula, m represents model parameters (porosity, clay content, etc.); d represents observed data (Y data volume); G() represents the forward modeling operator; C d Here is the data covariance matrix; m0 is the prior model; C m This is the model covariance matrix.

[0046] Algorithm flow: Under the hard data constraints of well points, sequentially traverse each grid node, and according to the conditional distribution (by C... d C m The reservoir parameters are randomly sampled and realized values ​​are generated with equal probability, while remaining faithful to seismic data, well data and geological laws.

[0047] S7, Fluid Detection and Reserve Assessment

[0048] The fluid spatial distribution prediction and resource estimation are performed using the sensitive factor data volume Y generated in step S4 and the reservoir parameter three-dimensional prediction volume generated in step S6.

[0049] Fluid detection and favorable area delineation: Analyze the planar distribution characteristics of the sensitive factor data volume Y (extracting slices along layers, anomaly amplitude, and anomaly area) and its spatial configuration relationship with the reservoir parameter volume. Based on known hydrocarbon reservoir response patterns, set reasonable threshold values ​​to identify and delineate fluid anomaly regions with significant hydrocarbon-bearing responses, i.e., favorable oil-bearing areas. These anomalies typically manifest as high-value (positive) anomalies in the Y data volume.

[0050] Among them, fluid detection is based on the Class IIA AVO response mechanism, and the gradient is calculated using the Shuey two-point method:

[0051] (4)

[0052] In the formula, R(θ) is the reflection coefficient at angle θ, and θ1 and θ2 are the angle ranges.

[0053] Based on rock physics experiments and calibration using actual well data, an empirical relationship between Y and oil saturation So is established:

[0054] S o=α·e (β·Y) α and β are regional calibration coefficients, determined through nonlinear regression of well logging data with known oil saturation in the target area and the Y-value of the wellbore. This exponential relationship reflects the nonlinear sensitivity of oil saturation to differences in seismic amplitude.

[0055] Geological Reserve Calculation: Within the delineated oil-bearing favorable area, using the three-dimensional reservoir parameter prediction volume obtained from step S5 inversion, combined with oil and gas reservoir engineering parameters (such as formation crude oil volume factor Bo, original dissolved gas-oil ratio Rs, etc.), the basic principle of the volumetric method is adopted to calculate the oil geological reserves unit by unit on the three-dimensional grid. Finally, the total geological reserve assessment results of the target area are obtained by summing them up. The calculation formula of the volumetric method is:

[0056] (5)

[0057] In the formula, OIP This refers to the geological reserves of crude oil. A The oil-bearing area; h The average effective thickness; Φ This represents the average effective porosity. S w The average water saturation level, B O This represents the formation crude oil volume factor. This three-dimensional gridded method for calculating reserves fully considers the spatial variation characteristics of reservoir parameters, greatly improving the accuracy of the assessment.

[0058] The present invention has the following beneficial effects:

[0059] (1) Significantly improves thin layer identification capability. By optimizing the five-angle stacking scheme, the thin layer identification rate is increased by 37% compared with the traditional three-angle method. It effectively identifies ultra-thin reservoirs with a thickness as low as λ / 10 (about 3 meters), breaking through the lower limit of conventional seismic identification.

[0060] (2) Effectively eliminates non-geological interference. By adopting the amplitude-frequency consistency correction algorithm, the system eliminates the difference in amplitude and main frequency between data from different angles caused by wave propagation path, absorption attenuation, etc., and significantly improves the consistency and comparability of multi-angle data.

[0061] (3) Achieve effective separation of lithology and fluid: by constructing a sensitive factor model (Y=X1- k X2), in which the near-angle data volume mainly reflects lithological information, and the far-angle data volume enhances fluid response, achieving high-precision decoupling of lithological background and fluid anomalies, with an oil-bearing prediction accuracy rate of 89.3%;

[0062] (4) The processing flow is efficient and has strong industrial applicability. The method has a clear process and key parameters (such as coefficients) are clearly defined. kBy optimizing well data, the complex and time-consuming process of full waveform inversion is avoided, which greatly saves computing resources and time costs and significantly improves the overall efficiency and operability of reservoir prediction work.

[0063] (5) Improved exploration success rate and economic benefits. In practical application, the drilling success rate reached 100%, the number of invalid wells decreased by 35%, the exploration cycle was shortened by 6 months, and the investment was saved by 150 million yuan, which significantly improved the drilling success rate and overall exploration benefits.

[0064] (6) It has formed an industrial application system, which integrates a series of key technologies such as data processing, inversion modeling and fluid detection. It has good repeatability and adaptability, and provides a reliable technical means to solve the problem of predicting complex thin interbedded reservoirs. It has high application value. Attached Figure Description

[0065] Figure 1 Method flowchart (showing the entire process from S1 to S7);

[0066] Figure 2 Pre-stack gather optimization processing (amplitude compensation, multiple suppression, and anisotropy correction are performed on the input pre-stack seismic gathers to eliminate waveform distortion caused by acquisition factors).

[0067] Figure 3 Data volumes of 3-13° (central angle 7°) are superimposed at various angles;

[0068] Figure 4 Data volumes of 10-20° (central angle 15°) are superimposed at different angles;

[0069] Figure 5 Data volumes of 17-27° (central angle 22°) are superimposed at different angles;

[0070] Figure 6 Data volumes of 24-34° (central angle 29°) are superimposed at various angles;

[0071] Figure 7 Data volumes of 31-41° (central angle 36°) are superimposed at various angles;

[0072] Figure 8 Consistency processing of data volumes superimposed at 31-41° angles (central angle 36°);

[0073] Figure 9 Comparison chart of spectrum analysis of data volumes with 3-13° (central angle 7°) and 31-41° (central angle 36°) superimposed at different angles;

[0074] Figure 10A comparison chart of the spectrum analysis of the 3-13° (central angle 7°) data volume superimposed at different angles and the 31-41° (central angle 36°) data volume after consistency processing;

[0075] Figure 11 AVO characteristic analysis of the target reservoir;

[0076] Figure 12 Analysis of the characteristics of seismic features superimposed from different angles on reservoirs;

[0077] Figure 13 Center impedance curves at different angles and the calculated throughput k The coefficients are used to obtain the EI curve reservoir sensitivity analysis diagram;

[0078] Figure 14 , k Principle diagram for coefficient determination (2.16EI7-EI) 36 Cross-plot analysis of the obtained EI curves).

[0079] Figure 15 Reservoir prediction results (comparison between actual drilled well lithology and predicted reservoir body);

[0080] Figure 16 Fluid detection results (comparison between actual drilling interpretation conclusions and fluid detection data). Detailed Implementation

[0081] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments in this specification without creative effort are within the scope of protection of this application.

[0082] The Weinan Depression, located on the northwestern edge of the continental shelf in the northern South China Sea, is structurally part of the central depression zone of the Beibu Gulf Basin and is a typical Cenozoic rift basin. This region has undergone a complete evolutionary process, experiencing the Paleocene-Eocene rifting stage, the Oligocene depression stage, and the Neogene regional subsidence stage, forming multiple source-reservoir-seal assemblages. Among them, the Weizhou Formation is one of the most important exploration targets in the area, deposited during the lacustrine basin expansion period, and featuring a large-scale delta-shoal-bar sedimentary system. After years of exploration, the Weinan Depression in the South China Sea has abundant usable data, with full coverage of 3D seismic data. Furthermore, in 2023, OBN seismic acquisition was deployed in the Weizhou Oilfield area, directly acquiring P-wave and S-wave seismic data from different azimuths and offsets, improving the imaging and interpretation accuracy of subsurface geological bodies. OBN seismic data also provides full coverage of key fault blocks and structures. The main geological problems faced in reservoir prediction in the work area include: (1) Difficulty in lithology identification: The difference in wave impedance between sandstone and mudstone in the third and fourth sections of Weishan is small. The velocity-density cross plot shows that the overlap of sandstone and mudstone parameters is more than 70%, and conventional wave impedance inversion is difficult to effectively distinguish lithology; (2) Insufficient resolution of thin interbedded layers: The target layer is mostly developed with thin interbedded layers of 1-5m, which is far lower than the seismic tuning thickness (λ / 4≈15-20m). The main frequency of seismic data is 32Hz, and the theoretically distinguishable lower limit of the stratum thickness is 12-15m. The actual thin layer response is submerged in the tuning effect; (3) Uncertainty in fluid detection: The difference in seismic response between oil and gas sandstone and water-bearing sandstone is slight, and the AVO characteristics overlap severely; rock physics analysis shows that the difference in P-wave impedance between oil layer and water layer is less than 5%, and the difference in Poisson's ratio is less than 8%; (4) Poor imaging quality in structurally complex areas: Near the No. 3 fault zone, the continuity of seismic phase axis is poor, and small faults are developed (fault displacement 10-20m), producing a large number of prediction artifacts. Strong anisotropy, difficult time difference correction; (5) rapid sedimentary facies change: rapid lateral change of beach bar sand bodies, thin thickness of single sand bodies, limited distribution, strong ambiguity in seismic attribute prediction. The accuracy needs to be verified, and there is a lack of more effective prediction methods to characterize concentrated sandstone development sections. In response to these geological problems, this invention proposes a complete reservoir prediction method based on multi-angle seismic differential enhancement. This method achieves high-precision identification and fluid detection of thin interbedded reservoirs through innovative data processing and interpretation techniques. The following is a detailed description of the implementation details of each step of this invention in conjunction with the actual application in the southwestern sea area of ​​Weihai. Figure 1 ).

[0083] S1. Data Preparation and Quality Control

[0084] Before the formal data processing flow was started, we systematically implemented comprehensive data preparation and quality control work to ensure the integrity, consistency and reliability of multi-source data. The data foundation of the work area consists of four categories of data: (1) Three-dimensional seismic data, using a 25m×25m grid, with acquisition parameters set to a 2ms sampling rate, a 4s recording length and a maximum offset of 4500m, and coverage times reaching 120 times in shallow layers and 80 times in deep layers, with the original data frequency band ranging from 8 to 80Hz; (2) Complete well logging series, collecting conventional well logging curves (natural gamma ray GR, sonic transit time AC, density DEN, neutron CNL, resistivity RT) and special well logging data (ECS elemental logging, DSI dipole sonic logging, MRI nuclear magnetic resonance) from 4 exploration wells; (3) Regional geological data, mainly including the top and bottom structural map of Weizhou Formation, fault distribution map, beach bar microfacies plane distribution map and core description and experimental analysis data of 120 meters from 4 wells. (4) Laboratory core analysis data and test data. The laboratory core analysis data covers key parameters such as porosity, permeability and saturation. The test data mainly consists of 12 DST test results, 20 PVT fluid sample analysis reports and system formation pressure test data, which provide multi-dimensional support for comprehensive geological interpretation and fluid identification.

[0085] In the data quality assessment and preprocessing stage, quantitative indicators were used to systematically evaluate the quality of seismic data. Analysis showed that the signal-to-noise ratio (SNR) of the raw seismic data exhibited significant longitudinal differences: greater than 3.0 for shallow layers, greater than 2.0 for mid-layers, and greater than 1.5 for deep layers; the dominant frequency distribution ranged from 28 to 35 Hz, with an effective bandwidth of 8–80 Hz; amplitude attribute analysis revealed significant energy differences among shallow, mid-, and deep layers, reaching 35 dB, primarily due to spherical geometric diffusion and formation absorption attenuation effects during seismic wave propagation. The correlation coefficient of gathers at different offsets was greater than 0.6, indicating good waveform consistency across the data. Well logging data underwent systematic standardization using histogram and trend surface analysis methods, achieving sonic transit time errors within 2 μs / ft and density errors less than 0.03 g / cm³. 3 The natural gamma error does not exceed 5 API. After standardization, the consistency of inter-well data is significantly improved, with the correlation coefficient increasing from 0.75 to 0.92. The time-depth relationship is established using a joint calibration method of VSP data and synthetic records. High-precision synthetic records with a main frequency of 35Hz and a frequency band of 10–80Hz are created, and multi-scale correlation algorithms are applied to improve calibration accuracy. Ultimately, the shallow calibration error is less than 2ms, and the deep calibration error is less than 5ms, achieving a time-depth relationship accuracy of R0. 2The value >0.995 lays a solid foundation for high-precision thin-layer interpretation and inversion. Rock physical analysis established a quantitative interpretation scale for the work area, clarifying the longitudinal wave velocity Vp of the mudstone as 2800–3200 m / s, the wave velocity Vs as 1400–1600 m / s, and the density ρ as 2.35–2.50 g / cm³. 3 The water-bearing sandstone has a velocity (Vp) of 3000–3400 m / s, a velocity (Vs) of 1600–1800 m / s, and a density (ρ) of 2.25–2.40 g / cm³. 3 The oil-bearing sandstone has a velocity (Vp) of 2700–3100 m / s, a velocity (Vs) of 1500–1700 m / s, and a density (ρ) of 2.20–2.35 g / cm³. 3 The study also shows that the Lamé coefficient, the product of density λρ and μρ, and Poisson's ratio σ are the most sensitive parameters for fluid and lithology identification, providing a reliable physical basis for subsequent seismic inversion and reservoir prediction.

[0086] S2, Pre-Stack Gather Optimization Processing

[0087] Pre-stack gather optimization is a fundamental and crucial component of our research methodology. Its core objective is to eliminate data distortions caused by non-geological factors such as acquisition, instrumentation, and wave propagation paths, thereby restoring the seismic wavefield response that accurately reflects the elastic properties of the subsurface medium. To achieve this goal, we adopted a step-by-step processing strategy, addressing amplitude, frequency, and time difference issues sequentially to improve data quality.

[0088] Firstly, for amplitude compensation, a series processing approach is adopted, sequentially implementing spherical diffusion compensation and Q-absorption compensation. Spherical diffusion compensation is based on wave theory and employs a time-varying gain method, calculated using the following formula:

[0089]

[0090] In the formula, G ( t (time) t Gain factor at; v 2 rms ( t )for t The root mean square velocity at time t, in m / s, is obtained through velocity analysis; c This is a constant, related to the initial amplitude and reference velocity; in this work area, it is taken as 2.5 × 10⁻⁶. 6 .

[0091] In actual processing, the velocity field was calculated at 50ms intervals with 100m×100m grid points, and the gain value was obtained by sliding within a 100ms time window, completing the compensation while maintaining the relative amplitude relationship. This step significantly reduced the energy difference between shallow, middle, and deep layers from 35dB to 15dB, and increased the effective signal energy of the deep layer by 3 times. Subsequently, based on the time-varying Q model established by VSP dispersion analysis: Q value 80–100 for shallow layer (0–1000m / s), Q value 60–80 for middle layer (1000–2000m / s), and Q value 40–60 for deep layer (>2000m / s), inverse Q filtering compensation was implemented, and its algorithm is expressed as follows:

[0092]

[0093] in: A ( f , t ) represents frequency f In time t Amplitude compensation factor at the location; Q ( t (time) t The quality factor at the location is obtained through dispersion analysis of VSP data; f Frequency, in Hz.

[0094] Within the 5–100 Hz frequency band, three iterations of compensation were performed with a 200 ms time window, which increased the energy of high-frequency components (>40 Hz) by 15 dB, expanded the effective frequency band from 8–80 Hz to 5–100 Hz, and increased the main frequency from 28 Hz to 35 Hz, significantly enhancing the identification potential of deep reservoirs.

[0095] Multiple suppression employs a combination of high-precision Radon transform and surface-correlated multiple attenuation (SRME). In the Radon transform, a parabolic Radon transform is used for velocity filtering in the τ-p domain to eliminate coherent energy with velocities below 1300 m / s, effectively suppressing seafloor rumbling and interlayer multiples commonly found in shallow water working areas. The transform formula is as follows: ,in, m ( τ , p () represents the Radon transform result; d (t, x ) represents the seismic gather; τ represents the intercept time.

[0096] The parameter settings included adaptive filtering with a main frequency of 35Hz and a bandwidth of 60Hz, followed by three iterations of optimization. This resulted in a 23dB energy attenuation for multiples while keeping the effective signal loss below 5%, and improving the signal-to-noise ratio from 2.5 to 4.0. The SRME method, on the other hand, uses a data-driven prediction model for multiples, supplemented by least-squares matched filtering for adaptive subtraction. This ensures a multiple removal rate exceeding 80% while maintaining an effective signal protection rate above 90%.

[0097] To address the time difference distortion problem caused by anisotropy, the Thomsen parameters are iteratively inverted using the layer-stripping method to obtain... ε =0.15、 δ =0.05、 γ =0.10, and based on this, the Tsvankin non-hyperbolic time difference correction formula is applied:

[0098] in: t ( x ) represents the two-way travel time (s) at the offset distance x; t 0 represents a zero-offset two-way trip, measured in seconds. v nmo NMO velocity, in m / s; η For anisotropic parameters, .

[0099] During the calibration process, velocity analysis was performed at every 100m×100m grid point, controlling the time difference accuracy within ±2ms, and the correlation of the stacked long-offset gathers was greater than 0.8. After calibration, the long-offset time difference distortion was reduced from more than 8ms to less than 2ms, the continuity of large-angle gathers was improved by 70%, and the fault relocation accuracy was improved to within 15m.

[0100] Finally, through comparative analysis of the Daoji process before and after ( Figure 2 The processed gather amplitude relative error is less than 5%, the main frequency variation does not exceed ±2Hz, the time difference correction error is less than 2ms, and the overall signal-to-noise ratio is improved by more than 50%, effectively ensuring the reliability of subsequent AVO analysis and inversion processing.

[0101] S3, Generation of angle-separated superimposed data volumes

[0102] Based on the pre-stack gather optimization process, this study innovatively designed and generated a five-angle stacked data volume. This scheme, based on detailed rock physics analysis and AVO forward modeling results, is significantly superior to the traditional three-angle partitioning method and can more accurately capture the differentiated response characteristics of lithology and fluids under different incident angles.

[0103] The design of the angle range is based on a solid foundation in rock physics. Through rock physical analysis of the Weizhou Formation, the elastic parameter ranges of various lithologies were clarified: the longitudinal wave velocity (Vp) of mudstone is 2800–3200 m / s, the transverse wave velocity (Vs) is 1400–1600 m / s, and the density (ρ) is 2.35–2.50 g / cm³. 3 The water-bearing sandstone has a velocity (Vp) of 3000–3400 m / s, a velocity (Vs) of 1600–1800 m / s, and a density (ρ) of 2.25–2.40 g / cm³. 3 The Vp of the oil-bearing sandstone is 2700–3100 m / s, the Vs is 1500–1700 m / s, and the ρ is 2.20–2.35 g / cm³. 3 Based on these parameters, a comprehensive AVO forward model was performed using the exact solution of the Zoeppritz equation. The simulation parameters included an incident angle range of 0–45° and an angle sampling interval of 1°. A forward model was also established that included sandstone-mudstone interfaces and sandstones with different fluid properties.

[0104] Based on forward modeling analysis, five optimal angle ranges were ultimately determined: the near-angle data volume (D07, 3–13°, central angle 7°) covers the AVO intercept sensitive area and has a sensitivity to P-wave impedance exceeding 85%, mainly used to reflect background lithological information; the mid-near-angle data volume (D15, 10–20°, central angle 15°) is located in the AVO response transition zone and is most sensitive to the transition zone response from Type I to Type II AVO characteristics; the mid-angle data volume (D22, 17–27°, central angle 22°) serves as the baseline data volume, targeting Type II... The critical reversal angle of the AVO-like response (25°±3°) has a signal-to-noise ratio as high as 12.8, and the sensitivity of polarity reversal detection for oil and gas sandstone is improved by 2.3 times. The mid-to-far angle data volume (D29, 24–34°, central angle 29°) further enhances the seismic response characteristics of fluids and improves the calculation accuracy of AVO gradients. The far angle data volume (D36, 31–41°, central angle 36°) enters the Poisson's ratio main control zone, with a sensitivity of 78% to Poisson's ratio changes caused by fluids, while effectively suppressing about 60% of lithological background interference.

[0105] During the generation of the stacked data volumes, each subset of gathers within each angular range underwent a rigorous preprocessing procedure. Dynamic correction was performed using an anisotropic NMO correction algorithm, with velocity analysis density reaching 100m×100m grid points at 50ms intervals, maximum stretching controlled within 30%, and maximum offset ratio set to 0.5, ensuring dynamic correction error was less than 2ms. Subsequent partial stacking processing employed differentiated strategies for different angular ranges: near-angle gathers were stacked 120 times, and far-angle gathers 80 times, strictly maintaining the relative amplitude relationship during stacking. This resulted in five sets of angularly stacked data volumes with excellent signal-to-noise ratios (near-angle SNR>8.0, far-angle SNR>5.0), namely: 3-13° (central angle 7°) angularly stacked seismic data volumes (…). Figure 3 ), 10-20° (central angle 15°) sub-angle stacked seismic data volume ( Figure 4 ), 17-27° (central angle 22°) angle-stacked seismic data volume ( Figure 5 ), 24-34° (central angle 29°) angle-stacked seismic data volume ( Figure 6 ), 31-41° (central angle 36°) sub-angle stacked seismic data volume ( Figure 7 ).

[0106] To ensure data volume quality, systematic quality monitoring was implemented. Angle gather quality monitoring showed that the angle range accuracy was controlled within ±2°, the relative error of amplitude fidelity was less than 5%, and the consistency error of in-phase axis time difference was less than 2ms. The quality assessment of the overlaid data volumes showed that each data volume had excellent signal-to-noise ratio (D07>10.0, D22>12.0, D36>8.0), stable dominant frequency characteristics (32±2Hz), and an effective bandwidth of 5–65Hz. Furthermore, the amplitude variation coefficient in the regional mudstone section was less than 0.15, fully demonstrating that this scheme significantly suppressed random noise while preserving effective geological information.

[0107] S4, Amplitude-Frequency Consistency Correction

[0108] After generating the multi-angle stacked data volumes, this study systematically implemented amplitude-frequency consistency correction to eliminate non-geological systematic biases caused by differences in seismic wave propagation paths. Because seismic waves at different incident angles experience different absorption and attenuation paths and anisotropic media during propagation, significant inconsistencies exist between data volumes from different angles, including amplitude energy attenuation, dominant frequency drift, and bandwidth variations. Specifically, the dominant frequency of the far-angle data volume decreases by 5–8 Hz, amplitude energy attenuates by up to 35%, bandwidth narrows by 10–15 Hz, and there is a time difference drift of up to ±4 ms. These factors severely affect the reliability of joint interpretation of multi-angle data, thus requiring consistency correction.

[0109] This study uses the mid-angle data volume D22 (center angle 22°) with the highest signal-to-noise ratio and most stable waveform as a benchmark, and adopts a staged processing strategy to systematically correct time difference, spectrum, and amplitude sequentially. First, time difference correction employs the Dynamic Time Warping (DTW) algorithm. This algorithm calculates the cumulative distance matrix and the Euclidean distance between two points to accurately determine the system time shift along key strata such as the Weishan 2nd, 3rd, and 4th segments. Experimental results show that the near-angle data volume D07 at the top boundary of the Weishan 2nd segment lags by 4 ms, while the far-angle data volume D36 leads by 8 ms. These time shifts are mainly due to the residual influence of anisotropy in the overlying strata. After correction using the DTW algorithm, the in-phase axis alignment error of all angle data volumes is less than 0.5 ms, and the inter-stratum time difference consistency coefficient is greater than 0.95.

[0110] Secondly, in the spectrum matching stage, amplitude spectrum unification is achieved by constructing a time-varying deconvolution operator. Using the time-varying amplitude spectrum of the reference data volume D22 as a template, spectrum shaping is performed on the target data volume, unifying the dominant frequency of all data volumes to 32±0.5Hz and extending the effective frequency band from 8–45Hz to 5–55Hz. This process strictly maintains the geological reflection structure characteristics within each seismic data volume, eliminating only the spectral differences with the reference volume, while simultaneously performing zero-phase processing to ensure the consistency of waveform characteristics.

[0111] Finally, amplitude equalization processing was performed. A time window of 100–200 ms was selected in the regional mudstone section to calculate the root mean square amplitude ratio of each angle data volume to the baseline data volume, obtaining the amplitude scaling factor, which was then applied to the entire data volume. After this processing, the energy attenuation rate of the far-angle data volume was significantly reduced from 35% to less than 8%, the amplitude variation coefficient was less than 0.1, and the overall amplitude fidelity was improved by up to 7 times, significantly improving amplitude fidelity. Figure 8 ).

[0112] To objectively evaluate the effectiveness of consistency correction, this study employed various quantitative indicators for verification. A comparative spectral analysis was conducted using data volumes superimposed at different angles: 3-13° (central angle 7°) and 31-41° (central angle 36°). Figure 9 ), and a comparison chart of the spectrum analysis of the 3-13° (central angle 7°) data volume superimposed at different angles and the 31-41° (central angle 36°) data volume after consistency processing ( Figure 10Analysis shows that the corrected dominant frequency deviation decreased from 7.2Hz to 0.4Hz, the bandwidth was effectively expanded, and the spectral correlation coefficient increased from 0.68 to 0.92. Amplitude consistency analysis shows that the amplitude variation coefficient decreased from 0.25 to 0.08, the cross-angle data correlation increased from 0.65 to 0.90, and the energy balance exceeded 0.95. Geological performance verification results further show that the well-seismic correlation coefficient increased from 0.75 to 0.92, the stratigraphic prediction error decreased from ±15m to ±8m, and the lithology identification accuracy increased from 70% to 88%, confirming the significant effectiveness of this correction method in improving seismic data quality and enhancing the reliability of geological interpretation.

[0113] S5, Calculation of reservoir sensitivity factors

[0114] The calculation of reservoir sensitivity factors is a core innovation in this methodology, aiming to construct a composite seismic attribute body that can effectively enhance the response of hydrocarbon-bearing sandstones while suppressing background lithological interference. Its physical basis lies in the characteristic differences in the variation of seismic reflection amplitude with incident angle (AVA / AVO): near-angle seismic data mainly reflects lithological information and has low sensitivity to fluid changes; far-angle data is extremely sensitive to changes in Poisson's ratio caused by pore fluids. Both rock physical analysis and actual data show that in the Weizhou Formation beach-bar sandstone reservoirs of the study area, high-quality reservoir sections exhibit a trend of significantly decreasing amplitude values ​​with increasing central angle in their seismic response. Figure 12 This characteristic is particularly prominent in oil and gas-bearing areas. Based on this, by constructing weighted difference attributes for near- and far-angle seismic data, effective separation of lithological background and fluid anomalies can be achieved.

[0115] This method is based on detailed rock physical analysis. Through cross-plot analysis of elastic impedance from four exploration wells in the work area, the characteristics of various lithologies at near-angle (EI7, central angle 7°) and far-angle (EI7) were clarified. 36 The spatial distribution of elastic impedance (central angle 36°) is as follows: mudstone exhibits high values ​​and a concentrated distribution; water-bearing sandstone shows a certain decreasing trend; while oil-bearing sandstone exhibits an EI value. 36 A significant decrease in a marked anomaly. The Weizhou Formation beach-bar sandstone exhibits typical Type II AVO characteristics ( Figure 11 Its intercept A is negative ( 0.05~ 0.02), and the gradient B is also negative ( 0.0008~ 0.0003), the product of intercept and gradient A×B is a positive value, which is an important indicator of hydrocarbon content.

[0116] Core coefficient kDetermining the value is a crucial step involving rigorous optimization. Its theoretical range (2.10–2.25) is derived from a linear approximation of the Zoeppritz equation. In the actual determination process, a refined calibration and optimization algorithm based on well data was employed: firstly, using wellside data from four wells (HX, HX1, HX2, HX3), near- and far-angle elastic impedance values ​​were extracted, establishing EI7 and EI... 36 Statistical Relationship Chart ( Figure 13 Furthermore, Fisher's discrimination criterion is used as the objective function to... k A systematic scan was performed on values ​​ranging from 2.0 to 2.3 with a step size of 0.01. Figure 14 The goal is to find the optimal value that maximizes the distinction between sandstone and mudstone. The field calibration results indicate that... k The value remained stable at 2.16±0.07 and had a clear petrophysical significance, particularly in the mudstone area. k The value is approximately 1.0, while in the sandstone area... k The value stabilized at 2.16, clearly reflecting the specific AVO gradient characteristics of the reservoir.

[0117] Based on certainty k The sensitivity factor Y value was calculated using the formula Y=X1-2.16X2 to generate a data volume of the entire work area, where X1 and X2 represent the near-angle (3–13°) and far-angle (31–41°) superimposed data volumes, respectively. Through calibration of the core data, a sensitive quantitative interpretation standard was established: Y>0.10 indicates high-quality sandstone (porosity >15%, thickness >3m); 0.10>Y>0 corresponds to medium-quality sandstone (porosity 10–15%, thickness 1–3m); 0>Y>-0.05 is interpreted as argillaceous sandstone; and Y<-0.05 is predicted as mudstone. In terms of fluid identification, this factor also showed good discriminative ability: the Y value for oil-bearing layers is typically between 0.15–0.25; for oil-water co-layers it is 0.10–0.15; for water-bearing layers it is 0.05–0.10; and for dry layers it is less than 0.05.

[0118] Practical application results show that the generated sensitive factor Y data volume is of excellent quality, with a vertical resolution improved to 3m (approximately λ / 10), a signal-to-noise ratio (SNR) higher than 8.0, and a distinguishing index between sandstone and mudstone reaching 0.86. Significant results have been achieved in geological applications, successfully identifying ultrathin reservoirs with a thickness of only 2.8m. The overall accuracy rate of oil-bearing prediction reached 85.7%, and the reservoir thickness prediction error was less than 1.2m, fully validating the effectiveness and practicality of this method in complex reservoir prediction and fluid detection.

[0119] S6, Inversion Processing

[0120] Inversion processing is a crucial step in achieving three-dimensional quantitative prediction of reservoir parameters. Its core objective is to transform the lithology- and fluid-sensitive factor data volume Y into a three-dimensional reservoir parameter model that is more familiar to geological researchers and can be directly applied, such as porosity and clay content. This study employs a stochastic inversion method based on a geostatistical framework, effectively integrating seismic derived attributes (Y data volume), well core and logging hard data, and prior knowledge of geological sedimentary models, achieving high-precision spatial prediction of reservoir parameters.

[0121] The inversion work is based on reliable rock physical relationships. Through core experimental data calibration, a quantitative conversion relationship between the sensitive factor Y and key reservoir parameters was established: the porosity conversion relationship is as follows: Φ =0.31Y+0.12, its goodness of fit Reaching 0.88; the conversion relationship of clay content is as follows: =0.22-0.28Y, The value is 0.82. These relationships significantly reduce the ambiguity of the inversion problem, providing a solid physical basis for parameter prediction based on seismic attributes.

[0122] Methodologically, the inversion is conducted within a formal Bayesian inference framework, and its objective function is expressed as:

[0123] (3)

[0124] Where m represents model parameters (porosity, clay content, etc.); d represents observation data (Y data volume); G() represents the forward modeling operator; C d Here is the data covariance matrix; m0 is the prior model; C m This is the model covariance matrix. This framework seeks a geologically sound and data-matched optimal solution by balancing the data fit difference with the model's prior information.

[0125] m represents model parameters (porosity, clay content, etc.); d represents observation data (Y data volume); G() represents the forward modeling operator; C d Here is the data covariance matrix; m0 is the prior model; C m The model covariance matrix

[0126] Spatial structure modeling was achieved using a variogram function, and this study employed a spherical model to characterize the spatial correlation of reservoir parameters. Based on the sedimentary pattern of the beach-bar sand bodies in the work area, the horizontal range was set to 3000 meters; based on the formation cycle thickness, the vertical range was set to 12 milliseconds (approximately 20 meters). Furthermore, a nugget value of 0.1 was introduced to characterize the heterogeneity of lithology at the microscale. The prior model weight λ was determined to be 0.35 after optimization through well-seismic calibration tests.

[0127] The inversion was implemented using a step-by-step optimization strategy. First, lithofacies modeling was performed. Using the Y data volume as a secondary variable, a three-dimensional sandstone-mudstone facies model was generated using the Sequential Indicator Simulation (SIS) method. The prior probability P(sand) for sandstone was set to 0.35, and Y=0.05 was used as the threshold value for lithological classification. Under the constraints of the obtained lithofacies model, the second step involved inverting physical property parameters. The Sequential Gaussian Simulation (SGS) algorithm was used to calculate the spatial distribution of porosity. The inversion grid was set to 25m×25m×2ms, with a sampling interval of 4ms, and a total of 100 implementations were performed to fully evaluate the solution space. Statistical analysis of the multiple implementation results assessed the uncertainty of the inversion. The results showed that the standard deviation of porosity prediction was ±2.5%, the thickness prediction error was ±1.2m, and the spatial distribution probability of areas with porosity greater than 12% exceeded 0.8 (…). Figure 15 This provides important quantitative basis for drilling decisions.

[0128] Verification of the inversion results demonstrates their high accuracy and reliability. Blind well testing (HX2 well) results show that the predicted reservoir top depth error is +1 meter, the predicted thickness error is -0.2 meters, the predicted porosity error is +0.5%, and the predicted clay content error is -2%. Cross-validation of four wells in the work area further confirms that the average absolute error of porosity prediction is 1.8%, the average relative error of thickness prediction is 12.5%, and the lithology identification accuracy rate reaches 88.3%. At the full 3D model level, the average absolute error of reservoir thickness prediction is 1.1 meters, the average absolute error of porosity prediction is 1.8%, the boundary error of oil-bearing area is less than 150 meters, and the final reserve calculation error is controlled within 8%, fully meeting the accuracy requirements of exploration and evaluation.

[0129] S7, Fluid Detection and Reserve Assessment

[0130] Building upon the reservoir parameter inversion, this study further conducted fluid detection and reserve assessment. Fluid detection employed an innovative dual-parameter identification mechanism, significantly improving the accuracy and reliability of hydrocarbon identification by comprehensively applying AVO gradient attributes and the sensitivity factor Y data volume. The AVO gradient calculation used the Shuey two-point method, utilizing the reflection coefficients at two angle points, 10° and 35°, to obtain the gradient value G. The oil-bearing criterion was defined as sensitivity factor Y > 0.18 and AVO gradient G < -0.15. This dual-parameter combination effectively overcomes the limitations of a single parameter, particularly eliminating false anomalies caused by the presence of mudstone, increasing the oil-bearing probability to over 85%.

[0131] Oil saturation modeling was based on rock physics experiments and calibration using actual well data. Through systematic rock electrical experiments, the ranges of formation parameters were determined: cementation index m = 1.8–2.2, saturation index n = 1.9–2.3, and rock electrical parameter a = 0.85–1.05. Based on these, an exponential relationship model between oil saturation and the sensitive factor Y was established. S o= ɑ·exp ( β·Y ),in, S o Oil saturation a , β For regional parameters, a =0.38, β =2.15 is a regional calibration coefficient obtained by nonlinear regression of well logging data with known oil saturation in the target area and the Y-value of the wellbore. This model fully reflects the nonlinear response relationship of oil saturation to the characteristics of seismic amplitude differences.

[0132] Reserve assessment employs a modified volumetric method, performing calculations cell-by-cell on a three-dimensional grid. Specific parameter determination methods are as follows: the oil-bearing area is delineated from the anomaly range of the Y data volume (Y>0.18); the effective thickness is extracted from the inversion results (…). Φ >12%, V sh <30%); porosity was obtained directly from the inverted porosity volume; water saturation was calculated from the above saturation model; formation crude oil volume factor B O We take 1.25, a value derived from high-pressure property analysis. The formula for the volumetric method is:

[0133] (5)

[0134] in, OIP This refers to crude oil geological reserves, expressed in tens of thousands of tons. A The oil-bearing area; h The average effective thickness; Φ This represents the average effective porosity. S w The average water saturation level, B O This represents the formation crude oil volume factor. This three-dimensional gridded method for calculating reserves fully considers the spatial variation characteristics of reservoir parameters, greatly improving the accuracy of the assessment.

[0135] Application results show that this method has achieved significant success in actual exploration. Both exploratory wells deployed in the Weixi B-3 block yielded industrial oil flow, achieving a 100% success rate, with an average daily oil production of 1012 tons per well, discovering recoverable reserves equivalent to a medium-sized oilfield. Regarding prediction accuracy, the fluid prediction thickness error is ±2m (…). Figure 16 The porosity prediction error is ±1.8%, the oil-bearing area error is less than 8%, and the reserve calculation error is controlled within 10%. All indicators have reached the leading level in the industry, providing a reliable basis for oilfield development decisions.

[0136] This invention addresses the challenge of predicting thin interbedded reservoirs in the Weizhou Formation in the southwestern Weihai Sea, developing a complete technical solution. Through multi-angle seismic differential enhancement technology, it effectively improves the accuracy of thin interbedded reservoir identification and fluid detection capabilities. Practical application shows that this method significantly improves drilling success rate and exploration efficiency, demonstrating significant value for widespread application.

[0137] The application of this invention in the southwestern sea area of ​​Weihai has achieved remarkable results. It successfully identified an ultrathin reservoir with a thickness of only 2.8m, far below the identification limit of conventional seismic methods. The reservoir prediction accuracy reached 89.3%, and two of the two deployed exploration wells obtained industrial oil flow, achieving a 100% success rate. Ineffective drilling was reduced by 35%, saving 150 million yuan in investment and shortening the exploration cycle by 6 months.

[0138] In summary, this invention provides a reservoir prediction method based on multi-angle seismic differential enhancement. By optimizing the design of five sets of angle-stacked data volumes (3-13°, 10-20°, 17-27°, 24-34°, 31-41°), and after amplitude-frequency consistency correction, an innovative sensitive factor model Y=X1-2.16X2 (X1 and X2 being near / far angle data volumes, respectively) is constructed. This model integrates inversion and fluid detection technologies to form an industrial application system. In the Weizhou Formation beach-bar sandstone strata of the Weixi region, the lower limit of reservoir identification thickness was successfully broken down to 3 meters (λ / 10), the oil-bearing prediction accuracy reached 89.3%, the drilling success rate increased from 62% to 86%, and invalid wells were reduced by 35%, effectively solving the industry problem of insufficient seismic prediction accuracy for thin interbedded marginal oil reservoirs.

[0139] Those skilled in the art should understand that the above embodiments are merely illustrative and are not intended to imply that the scope of the invention is limited to these examples. Within the framework of this invention, technical features of the above embodiments or different embodiments can be combined, steps can be implemented in any order, and many other variations of the different aspects of the invention as described above exist, which are not provided in detail for the sake of brevity. Any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A reservoir prediction method based on multi-angle seismic differential enhancement, characterized in that, Includes the following steps: S1. Data preparation and quality control: Collect data on the study area, including regional geological data, 3D seismic data, well logging curves, analytical and testing data, and evaluate the data quality. S2. Pre-stack gather optimization processing: This process performs fine pre-processing on the input raw seismic pre-stack gathers to reduce waveform distortion and energy distortion caused by non-geological factors. The pre-processing includes amplitude compensation, multiple suppression, and anisotropy correction. S3. Generation of angle-separated stacked data volumes: Based on the pre-stack gathers optimized in S2, and considering the influence of seismic wave incident angle on the response characteristics of underground reflection points, each trace within the gather is divided according to its corresponding incident angle range. The incident angle is divided into five continuous and partially overlapping incident angle range gather subsets, namely 3-13°, 10-20°, 17-27°, 24-34°, and 31-41°, with corresponding central angles of 7°, 15°, 22°, 29°, and 36°, respectively. For each divided incident angle range gather subset, dynamic correction NMO / DMO, muting, and stacking processing are performed respectively, finally generating five sets of partially stacked seismic data volumes representing different average incident angles, i.e., central angle reflection characteristics. S4. Amplitude-frequency consistency correction: The data volume in the range of 17-27° is selected as the reference volume. The Spectral Matching algorithm is used. The amplitude spectrum and phase spectrum of the reference data volume are used as templates to adjust the amplitude energy and align the main frequency characteristics of the other four sets of data volumes channel by channel and time window by time window. Finally, five sets of angle-separated superimposed data volumes with high consistency in amplitude energy level and main frequency characteristics are obtained. S5. Reservoir Sensitivity Factor Calculation: Based on the rock physical characteristics of the target reservoir in the study area, characteristic near-angle stacked data volume X1 and far-angle stacked data volume X2, after consistency correction, are extracted. X1 corresponds to a central angle of 7°, and X2 corresponds to a central angle of 36°. Utilizing the difference enhancement characteristics between these two types of data volumes, the sensitivity factor is calculated according to the formula Y=X1- k X2 calculations generate a sensitive factor data volume Y that highlights the reservoir's hydrocarbon response, including key coefficients. k Its value is not a fixed constant, but is obtained through regression analysis using the elastic impedance (EI) information at known well points in the target area; S6. Inversion Processing: The sensitive factor data volume Y obtained in step S5 is used as the core input data. Constrained inversion is performed by combining well point lithology data, sedimentary mode constraints, and multi-source information on rock physical relationships. The inversion algorithm is used to generate a high-precision three-dimensional prediction volume of reservoir parameters, realizing a quantitative description of reservoir spatial distribution and physical properties. S7. Fluid detection and reserve assessment: Using the sensitive factor data volume Y generated in step S4 and the three-dimensional reservoir parameter prediction volume generated in step S6, fluid spatial distribution prediction and resource estimation are performed.

2. The reservoir prediction method based on multi-angle seismic differential enhancement according to claim 1, characterized in that, In step S3, the incident angle range gather subset design strictly follows the AVO / AVA theory and the target reservoir rock physical response mechanism.

3. The reservoir prediction method based on multi-angle seismic differential enhancement according to claim 1, characterized in that, In step S4, during amplitude energy adjustment and main frequency feature alignment, the original geological reflection structure inside each data volume is strictly maintained, and only the spectral difference between it and the reference volume is corrected.

4. The reservoir prediction method based on multi-angle seismic differential enhancement according to claim 1, characterized in that, In step S5, the regression analysis involves extracting the near-angle elastic impedance EI1 and the far-angle elastic impedance EI2 at the wellbore access location, establishing a linear regression relationship between EI1 and EI2, and the slope of this regression relationship is the coefficient. k ,coefficient k The method for determining the target layer is a statistical optimization process based on the physical properties of the rock. The specific steps are as follows: ① Wellside elastic impedance extraction: At a known well location, for the target reservoir section, extract the corresponding elastic impedance curves EI7 and EI from the corrected near-angle 7° stacked data volume and the far-angle 36° stacked data volume, respectively. 36 This extraction process needs to be combined with well-seismic calibration to ensure accurate time-depth relationship; ② Sandstone sensitive impedance modeling: Based on rock physical analysis of the target area, a multiple regression equation is constructed to maximize the distinction between sandstone and mudstone. （2） In the formula, EI opt To obtain the sensitive elastic impedance of the target lithology body with multi-angle differential body, EI 7 represents the elastic resistance at a center angle of 7°. EI 36 For an elastic impedance with a center angle of 36°, a , b This is an empirical coefficient; ③ Coefficient k Optimization determined that the sensitivity factor formula Y=X1- k X2 is analogous to the impedance domain, let EI=EI 7- k · EI 36 Through systematic change k Values, calculations differ k The EI curves at the specified values ​​were ultimately selected to maximize the distinction between sandstone and mudstone. k The value is the optimal one. k value.

5. The reservoir prediction method based on multi-angle seismic differential enhancement according to claim 1, characterized in that, In step S7, the resource quantity estimation uses the volumetric method, and the calculation formula is as follows: （5） In the formula, OIP This refers to the geological reserves of crude oil. A The oil-bearing area; h The average effective thickness; Φ This represents the average effective porosity. S w The average water saturation level, B o This represents the volume coefficient of crude oil in the formation.

Images