A method for correcting the amplitude of seismic waves reflected from inclined interfaces based on forward modeling of wave equation

Abstract

The invention relates to a method for correcting the amplitude of seismic waves reflected from an inclined interface based on wave equation forward modeling, and belongs to the technical field of reservoir prediction in oil and gas field exploration. The method firstly establishes three geological models with different dip characteristics based on three-dimensional post-stack seismic data of a study area, logging data of a completed well, and seismic layer data of a target layer and its upper portion, then adopts a wave equation forward modeling method to perform seismic forward modeling, extracts the amplitude and dip data of seismic waves reflected from the target layer on the forward modeling record, adopts the least squares method to fit to obtain a quadratic fitting function of amplitude and dip, then adopts the quadratic fitting function of amplitude and dip obtained by forward modeling analysis to calculate a simulated reflected seismic wave amplitude array of the target layer, corrects a real reflected wave amplitude array of the target layer, and obtains an amplitude array of the target layer that can be effectively used for reservoir prediction.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas field exploration, and in particular to a method for correcting amplitude of seismic event axes of inclined interfaces based on wave equation forward modeling. Background Art

[0002] In seismic exploration in the field of oil and gas field exploration, reservoirs are the main target of search. For reservoirs developed near the top, their response characteristics usually weaken the seismic event amplitude in the reservoir development area due to the existence of the strong reflection interface at the top. In the case of horizontal or small formation inclination, the seismic event amplitude of the top interface can qualitatively indicate the development of the top reservoir. Conventional amplitude correction methods (wavelet decomposition and reconstruction, etc.) highlight the seismic response characteristics of the reservoir developed near the top by removing the influence of the strong seismic reflection of the top interface without considering the formation inclination. This technology has achieved good application results.

[0003] However, in inclined strata, in addition to considering the impact of the difference in wave impedance above and below the interface, it is also necessary to consider the impact of the inclination of the strata on the amplitude of the reflection event. The reflection event amplitude caused by the inclined interface becomes weaker, and its amplitude is often much greater than the impact of the reservoir on the reflection event. At this time, the accuracy of predicting the reservoir using the conventional top interface reflection amplitude attribute is low. At present, the research on the method of removing the top interface seismic strong reflection from the horizontal or small-angle reflection interface for reservoir prediction at home and abroad is mainly focused on the matching pursuit-based wavelet decomposition and reconstruction method, which is the mainstream method with good results. Recently, domestic scholars have proposed a wavelet decomposition and reconstruction technology based on a strong reflection model, which has achieved good application results when the distance between the reservoir and the strong reflection interface is less than 1 / 4 wavelength of the wavelet.

[0004] 1. Wavelet decomposition and reconstruction technology based on matching pursuit

[0005] This method is based on a multi-wavelet seismic trace model. Through the decomposition of seismic wavelets, the seismic trace is decomposed into a set of wavelets with different main frequencies and different amplitudes, and the corresponding relationship between the reflection coefficients of different strata and the seismic wavelets is obtained. Through seismic wavelet reconstruction, all or part of the wavelets with different main frequencies and different amplitudes obtained by wavelet decomposition are selected to be re-superimposed to form a new seismic trace, and the effective signal with weaker energy in the seismic reflection is strengthened. The wavelet decomposition and reconstruction based on the matching pursuit algorithm projects the original signal onto a series of orthogonal atoms, and the seismic signal is represented by a linear combination of different atoms. By analyzing and screening different atoms and then reconstructing them, the shielding effect of the strong reflection interface on the underlying strata is corrected.

[0006] 2. Wavelet decomposition and reconstruction technology based on strong reflection model

[0007] Based on the strong reflection model, domestic scholars cross-correlated the seismic records and superimposed them along the source direction. After convolving the superposition results with the original data, they superimposed them along the detector direction to obtain the reconstructed first arrival wave. Compared with the original signal, the strong reflection first arrival signal and the weak reflection first arrival signal can be distinguished, and the amplitude spectrum of the weak reflection signal can be enhanced. The reconstructed first arrival signal is sparsely represented and the basis matrix is ​​constructed for solution. For the seismic signals around the strong reflection, after the sparse decomposition of the signal, it is only necessary to cut off the components with larger amplitudes in the dynamic dictionary around the strong reflection to achieve the purpose of removing the strong reflection.

[0008] 3. Deficiencies of existing methods

[0009] 1. The wavelet decomposition and reconstruction based on the matching pursuit algorithm is mainly based on the multi-wavelet seismic trace model. When the reflection coefficient and seismic wavelet are unknown, the matching pursuit algorithm takes longer and requires more calculations. The matching pursuit algorithm assumes that atoms are orthogonal to each other. The reservoir developed on the top of the inclined interface will be affected by the interference of seismic reflection waves because the range of the reservoir from the top interface changes vertically and horizontally. The orthogonality of the atoms cannot be guaranteed, and the atoms directly decomposed by the matching pursuit method on the original data body will produce errors. The amplitude changes vertically and horizontally due to the tilt of the interface, and the matching pursuit algorithm will reduce the correlation of the amplitude components of atoms near the same interface, and reduce the accuracy of the calculation.

[0010] 2. When the inclination angle of the formation interface changes greatly laterally, the reliability of the result in the calculation time window of the wavelet decomposition and reconstruction based on the strong reflection model is reduced, and the error of the reconstructed first arrival wave increases. Since the amplitude of the phase axis of the inclined interface in the calculation time window is not uniformly distributed, the large amplitude component removed from the signal after sparse decomposition will lead to greater errors.

[0011] Therefore, in order to address the phenomenon of seismic wave amplitude variation at inclined interfaces, that is, the phenomenon that the amplitude of seismic reflection event axes varies unevenly along the same wave impedance interface and is greatly affected by the degree of inclination, the present invention has invented a method for correcting the amplitude of seismic waves reflected at inclined interfaces based on forward modeling of wave equations, and a correction function of the target layer interface inclination array to the reflection seismic wave amplitude array is established in a certain area. Summary of the invention

[0012] The present invention relates to a method for correcting the amplitude of seismic waves reflected from an inclined interface based on wave equation forward modeling, and belongs to the technical field of reservoir prediction in oil and gas field exploration. The method first establishes three geological models with different dip characteristics based on three-dimensional post-stack seismic data of the study area, logging data of completed wells, and seismic horizon data of the target layer and above, and then uses the elastic wave equation post-stack forward modeling method to perform seismic forward modeling, and extracts the amplitude and dip data of the seismic waves reflected from the target layer on the forward modeling record, and uses the least squares method to fit to obtain a quadratic fitting function of the amplitude and dip; then extracts the amplitude array and formation dip array of the target layer interface of the three-dimensional post-stack seismic data along the seismic horizon of the target layer, and uses the quadratic fitting function of the amplitude and dip obtained by forward simulation analysis to calculate the simulated reflected seismic wave amplitude array of the target layer, and corrects the real reflected wave amplitude array of the target layer to obtain an amplitude array of the target layer that can be effectively used for reservoir prediction.

[0013] The specific steps of the present invention include:

[0014] (1) Input the 3D post-stack seismic data of the study area and the seismic horizons of the target layer and its upper strata, calculate the stratigraphic dip of the target layer interface, and then select the seismic profiles of the target layer with three different dip angles. The dip angle characteristics of the target layer are seismic profiles with large dip angle, medium dip angle and small dip angle. Based on the seismic horizons of the target layer and its upper strata, establish three corresponding stratigraphic framework models.

[0015] (2) Input the logging data of the wells drilled through the target layer in the study area, statistically analyze the interval velocity and density of the target layer and the upper strata, and fill them into the three stratigraphic grid models in step 1 to obtain three geological models for forward simulation;

[0016] (3) using the wave equation numerical simulation method, forward modeling the three geological models in step 2 is performed to obtain three corresponding simulated seismic records, and the amplitude array of the reflected seismic waves of the target layer of the three simulated records is extracted along the seismic horizon of the target layer;

[0017] (4) Based on the inclination of the target layer in step 1 and the amplitude array in step 3, a quadratic fitting function of amplitude and inclination is obtained by using the least squares method;

[0018] (5) extracting the amplitude array and formation dip array of the target layer interface of the three-dimensional post-stack seismic data in step 1 along the seismic horizon of the target layer, and normalizing the amplitude array to obtain the real reflected seismic wave amplitude array of the target layer;

[0019] (6) Based on the dip array in step 5 and the amplitude-dip fitting function in step 4, the amplitude array corresponding to the dip array in step 5 is calculated and normalized, i.e., the simulated reflected seismic wave amplitude array of the target layer;

[0020] (7) using the simulated reflected seismic wave amplitude array of the target layer in step 6 to correct the real reflected seismic wave amplitude array of the target layer in step 5, the reflected seismic wave amplitude array of the target layer after correction = the real reflected seismic wave amplitude array of the target layer + 1.0 - the simulated reflected seismic wave amplitude array of the target layer;

[0021] (8) Outputting the reflected seismic wave amplitude array of the target layer after correction in step 7 for reservoir prediction based on amplitude variation characteristics.

[0022] A method for correcting the amplitude of tilted interface seismic event axis based on wave equation forward modeling has the following characteristics, which are mainly manifested as follows:

[0023] (1) This invention, for the first time in the field of post-stack reservoir prediction, highlights the true response characteristics of the lithology combination by correcting the influence of the formation dip angle on the amplitude of the reflection wave event axis. Through the forward modeling results of the wave equation of the geological model, the relationship between the formation dip angle array and the target layer interface reflection wave amplitude array is established. The original interface reflection wave amplitude array is corrected by calculating the forward modeling reflection wave amplitude array of the actual target layer interface dip angle, thereby achieving the amplitude correction of the seismic reflection interface, which is highly innovative.

[0024] (2) The core of the present invention lies in the necessity of formation dip correction, that is, whether the correction result can accurately indicate the response characteristics (such as amplitude, etc.) of the reservoir developed near the top interface. After investigating a large number of literatures, it is found that the superposition and offset processing of the pre-stack CMP gather cannot well superimpose the reflected wave energy of the inclined interface. This leads to the impedance interface with the same reflection coefficient in the post-stack seismic data. Due to different degrees of inclination, the amplitude strength of the event axis on the post-stack seismic profile is different. The amplitude attribute of the top interface of the reservoir developed near the top is often used to indicate the development of the top reservoir. Due to the superposition and offset of the pre-stack CMP gather processing, the response characteristics of the top reservoir of the inclined interface are weakened without amplitude correction. At this time, the conventional amplitude attribute extraction cannot accurately reflect the real impedance difference above and below the interface. The present invention constructs a forward geological model based on actual three-dimensional post-stack seismic data, seismic layer data and logging data, and analyzes the influence of the dip angle on the interface reflection amplitude by controlling the change of the single factor of the formation dip angle, thereby establishing a suitable correction function, which has a good methodological theoretical basis and feasibility. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 It is a technical flow chart of the present invention;

[0026] Figure 2 Seismic profiles selected to assist in explaining the construction of the geological model;

[0027] Figure 3To assist in explaining the forward geological model constructed based on seismic sections, seismic horizons and well logging data;

[0028] Figure 4 To assist in explaining the post-stack seismic profile obtained by forward modeling of the forward geological model using the wave equation;

[0029] Figure 5 This is the post-stack reservoir prediction plane map before the top interface reflection amplitude array correction;

[0030] Figure 6 This is the post-stack reservoir prediction plan after the dip amplitude correction of the top interface reflection amplitude array. DETAILED DESCRIPTION

[0031] Example 1

[0032] A method for correcting the amplitude of seismic waves reflected from an inclined interface based on forward modeling of a wave equation, comprising the following steps:

[0033] Step 1: Input the 3D post-stack seismic data volume SA(x, y, t) of the study area, the seismic horizon HT(x, y) of the target layer in the study area, and the seismic horizon HT(x, y) of the overlying stratum of the target layer. 1 (x,y),HT 2 (x,y),…,HT nc (x, y), nc is the total number of strata above the target layer in the study area, which is determined by the stratigraphic structure of the study area;

[0034] Step 2: Use the instantaneous frequency and instantaneous wavenumber in the x and y directions in the complex channel analysis to calculate the apparent dip array Px(x,y,t) in the x direction and the apparent dip array Py(x,y,t) in the y direction of the three-dimensional data body SA(x,y,t). The calculation method adopts the method disclosed in the article "Application of three-dimensional multi-scale volume curvature in fracture identification of deep granite buried hill reservoirs" by Zeng Xianghao et al. (2023), and calculate the dip array P(x,y,t).

[0035]

[0036] Step 3, using the target layer seismic layer array HT(x,y) to extract the dip array HP(x,y) of the target layer in the dip array P(x,y,t), and according to the HP(x,y) value distribution range, the dip of the target layer is divided into three scales: large, medium, and small. In the three-dimensional post-stack seismic data volume, three y-direction sections Sx1(x,y,t), Sx2(x,y,t) and Sx3(x,y,t) with obvious dip characteristics of the target layer are selected. Based on the stratigraphic layer array of the target layer in the study area and the stratigraphic layers above it, three stratigraphic framework models M1(x,t), M2(x,t) and M3(x,t) corresponding to the sections Sx1(x,y,t), Sx2(x,y,t) and Sx3(x,y,t) are established.

[0037] Step 4: Input the logging data of the completed wells that penetrate the target layer in the study area, use the average value analysis method to statistically analyze the layer velocity and density of the target layer and the upper strata, and fill them into the stratigraphic framework model in step 3 to obtain three velocity models M1_vel(x, t), M2_vel(x, t), M3_vel(x, t) and density models M1_den(x, t), M2_den(x, t), M3_den(x, t) for forward simulation;

[0038] Step 5, based on the velocity model and density model, obtain the simulation records of the three models, M1_s(x, t), M2_s(x, t) and M3_s(x, t). The calculation method adopts the method disclosed in the article "Earthquake Structure Detailed Interpretation Technology Based on Wave Equation Numerical Simulation" by Xiong Xiaojun et al. in 2011, and uses the target layer seismic layer array HT(x, y) and the dip array P(x, y, t) of step 2 to extract the amplitude array M1_amp(x, t), M2_amp(x, t), M3_amp(x, t) and dip array M1_dip(x), M2_dip(x), M3_dip(x) of the target layer in the simulation records of the three models;

[0039] Step 6, based on the amplitude array and dip array of the target layer in step 5, M1_amp(x,t), M2_amp(x,t), M3_amp(x,t) and M1_dip(x), M2_dip(x), M3_dip(x), the least squares method is used to fit the amplitude and dip, and the fitting coefficients a, b and c in formula 2 are obtained by using the method disclosed by Gao Qiuying et al. in the article "Research on Least Squares Curve Fitting and Optimization Algorithm" in 2021.

[0040] amp(θ)=a×θ 2 +b×θ+c (2)

[0041] Step 7, using the target layer seismic layer array HT(x,y), the dip array P(x,y,t) of step 2 and the three-dimensional post-stack seismic data volume SA(x,y,t) of the study area in step 1, extract the amplitude array HAmp(x,y) and dip array HDmp(x,y) of the target layer, and use the fitting relationship in formula 2 to obtain the amplitude array HAmp_S(x,y),

[0042] HAmp_S(x,y)=a×HDmp 2 (x,y)+b×HDmp(x,y)+c (3)

[0043] Step 8: Normalize the arrays HAmp(x,y) and HAmp_S(x,y).

[0044] HAmp ′ (x,y)=HAmp(x,y)÷max(HAmp(x,y)) (4)

[0045] HAmp_S ′ (x,y)=HAmp_S(x,y)÷max(HAmp_S(x,y)) (5)

[0046] Among them, HAmp ′ (x,y) and HAmp_S ′ (x,y) is the normalized array, and max() is the function for calculating the maximum value;

[0047] Step 9, use HAmp_S from step 8 ′ (x,y) array pair HAmp ′ (x,y) is corrected,

[0048] HAmp_E(x,y)=HAmp ′ (x,y)+1.0-HAmp_S ′ (x,y) (6)

[0049] Among them, HAmp_E(x,y) is the reflected seismic wave amplitude array of the target layer after formation dip correction,

[0050] Step 10, outputting the reflected seismic wave amplitude array of the target layer corrected in step 9, for reservoir prediction based on amplitude variation characteristics.

[0051] Example 2

[0052] Figure 5 and Figure 6 They are examples of post-stack reservoir prediction for a target layer in a certain study area. Figure 5It is a plan view of reservoir prediction based on the reflection amplitude array extracted from the three-dimensional seismic data of the study area along the interface of the target layer. The result is only 67% consistent with the reservoir development of the completed wells. Figure 6 This is a reservoir prediction plane diagram of the target layer reflection amplitude array obtained by using the present invention to perform dip correction. The result is consistent with the reservoir development situation after drilling at a rate of 92%. It can be seen that the reservoir prediction effect of the present invention is better.

[0053] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for correcting the amplitude of seismic waves reflected from an inclined interface based on forward modeling of a wave equation comprises the following specific steps: Step 1: Input the 3D post-stack seismic data volume SA(x,y,t) of the study area, input the seismic horizon HT(x,y) of the target layer in the study area and the seismic horizons HT1(x,y), HT2(x,y), ..., HT nc (x, y), nc is the total number of strata above the target layer in the study area, which is determined by the stratigraphic structure of the study area; Step 2, using the instantaneous frequency and instantaneous wave number in the x and y directions in the complex channel analysis, calculate the x-direction apparent dip array Px(x,y,t) and the y-direction apparent dip array Py(x,y,t) of the three-dimensional data body SA(x,y,t), and calculate the dip array P(x,y,t). Step 3, using the target layer seismic layer array HT(x,y) to extract the dip array HP(x,y) of the target layer in the dip array P(x,y,t), and according to the HP(x,y) value distribution range, the dip of the target layer is divided into three scales: large, medium, and small. In the three-dimensional post-stack seismic data volume, three y-direction sections Sx1(x,y,t), Sx2(x,y,t) and Sx3(x,y,t) with obvious dip characteristics of the target layer are selected. Based on the stratigraphic layer array of the target layer in the study area and the stratigraphic layers above it, three stratigraphic framework models M1(x,t), M2(x,t) and M3(x,t) corresponding to the sections Sx1(x,y,t), Sx2(x,y,t) and Sx3(x,y,t) are established. Step 4: Input the logging data of the completed wells that penetrate the target layer in the study area, use the average value analysis method to statistically analyze the layer velocity and density of the target layer and the upper strata, and fill them into the stratigraphic framework model in step 3 to obtain three velocity models M1_vel(x, t), M2_vel(x, t), M3_vel(x, t) and density models M1_den(x, t), M2_den(x, t), M3_den(x, t) for forward simulation; Step 5, based on the velocity model and density model, obtain the simulation records of the three models, M1_s(x, t), M2_s(x, t) and M3_s(x, t), and use the target layer seismic layer array HT(x, y) and the dip array P(x, y, t) of step 2 to extract the amplitude array M1_amp(x, t), M2_amp(x, t), M3_amp(x, t) and dip array M1_dip(x), M2_dip(x), M3_dip(x) of the target layer in the simulation records of the three models; Step 6: Based on the amplitude array and dip array of the target layer in step 5, M1_amp(x, t), M2_amp(x, t), M3_amp(x, t) and M1_dip(x), M2_dip(x), M3_dip(x), the least squares method is used to fit the amplitude and dip to obtain the fitting coefficients a, b and c in formula 2. amp(θ)=a×θ 2 +b×θ+c (2) Step 7, using the target layer seismic layer array HT(x,y), the dip array P(x,y,t) of step 2 and the three-dimensional post-stack seismic data volume SA(x,y,t) of the study area in step 1, extract the amplitude array HAmp(x,y) and dip array HDmp(x,y) of the target layer, and use the amplitude array HAmp_S(x,y) of the fitting coefficient in formula 2, HAmp_S(x,y)=a×HDmp 2 (x,y)+b×HDmp(x,y)+c (3) Step 8: Normalize the arrays HAmp(x,y) and HAmp_S(x,y). HAmp′(x,y)=HAmp(x,y)÷max(HAmp(x,y)) (4) HAmp_S′(x,y)=HAmp_S(x,y)÷max(HAmp_S(x,y)) (5) Among them, HAmp′(x,y) and HAmp_S′(x,y) are normalized arrays, and max() is the maximum value calculation function; Step 9: Use the HAmp_S′(x,y) array from step 8 to calibrate HAmp′(x,y). HAmp_E(x,y)=HAmp′(x,y)+1.0-HAmp_S′(x,y) (6) Among them, HAmp_E(x,y) is the reflected seismic wave amplitude array of the target layer after formation dip correction; Step 10, outputting the reflected seismic wave amplitude array of the target layer corrected in step 9, for reservoir prediction based on amplitude variation characteristics.

Images