A method for extracting a seismic wavelet in a depth domain based on a generalized seismic wavelet model

Abstract

The application discloses a depth domain seismic wavelet extraction method based on a generalized seismic wavelet model. The method is characterized in that, the depth domain logging data information can be fully utilized, the depth domain seismic record is directly synthesized in the depth domain, the depth domain seismic wavelet can be efficiently and accurately extracted from the depth domain generalized seismic wavelet model only through 3 to 5 iterations, the method has high stability, and the effective information loss problem caused by domain conversion or data resampling is avoided.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas seismic exploration, specifically relating to a depth domain seismic wavelet extraction method based on a generalized seismic wavelet model. Background Technology

[0002] In the field of oil and gas seismic exploration, with the rapid development and widespread application of pre-stack depth-domain migration imaging, directly realizing reservoir parameter inversion, oil and gas reservoir prediction and description in the depth domain is an urgent practical need. Extracting seismic wavelets from depth-domain seismic data is an indispensable and fundamental step in achieving depth-domain seismic inversion. As a basic component of seismic data, the characteristics of seismic wavelets are of great significance for revealing subsurface structures. In the process of seismic data processing, the analysis and application of seismic wavelets run through many key stages, including but not limited to deconvolution, well-seismic calibration, and forward and inverse modeling of seismic records. However, due to the non-stationarity of depth-domain seismic data, traditional convolution model methods based on the assumption of linear depth-invariant systems cannot be directly applied to the depth domain, posing a significant challenge to depth-domain seismic wavelet extraction.

[0003] To address this challenge, researchers are constantly exploring new methods and technologies to accurately and efficiently extract depth-domain seismic wavelets, thereby better serving the needs of oil and gas exploration and development. Currently, depth-domain seismic wavelet extraction mainly employs two methods: (1) converting depth-domain seismic data to the time domain or constant velocity depth domain, extracting the seismic wavelets, and then converting them back to the depth domain. This domain conversion method inevitably introduces information loss and accumulated errors into both depth-domain seismic data and depth-domain logging data; (2) using depth-wavenumber spectrum analysis technology to directly extract depth-domain seismic wavelets at each depth location. However, this method lacks the constraints of depth-domain logging data, resulting in multiple solutions in the extraction results, and requires processing seismic data channel by channel, which is time-consuming and inefficient. Since depth-domain seismic data does not possess the "linear depth-invariant" characteristic, it is impossible to directly apply convolution model-like methods to synthesize depth-domain seismic records. Therefore, this invention aims to better solve this problem by enabling direct, rapid, and efficient extraction of depth-domain seismic wavelets in the depth domain, providing a fundamental method for depth-domain reservoir parameter inversion and fluid identification. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention proposes a depth-domain seismic wavelet extraction method based on a generalized seismic wavelet model. This method effectively solves the problems of unreliability of seismic wavelets extracted through domain transformation and the time-consuming and multi-solution issues associated with using depth-wavenumber spectrum extraction. The method of this invention includes the following main steps:

[0005] (1) The depth domain logging reflection coefficient r is calculated using the depth domain logging density ρ and the longitudinal wave velocity Vp. The number of sampling points for the logging reflection coefficient is N1, and the sampling interval in the depth domain is Δh1.

[0006] (2) Extract the single-track seismic record S in the depth domain next to the well. The number of sampling points for the single-track seismic record in the depth domain next to the well is N2, and the depth sampling interval is Δh2.

[0007] (3) Using the three-sample interpolation method with a well depth domain sampling interval of Δh1, the original seismic data is resampled to obtain the resampled seismic data S′.

[0008] (4) Set the length N3, velocity Vc, and sampling interval Δh of the depth domain seismic wavelet to be extracted. wave ;

[0009] (5) Set the reference wavenumber search range K and order range U for the depth domain generalized seismic wavelet to be extracted, and construct the parameter grid P(i,j)=[K(i),U(j)] for seismic wavelet extraction;

[0010] (6) The three parameters of the depth-domain generalized seismic wavelet, namely the reference wave number k0, the order u, and the amplitude a, are solved by the alternating relaxation iterative method as follows:

[0011] ① Set the maximum number of iterations L and the minimum iteration error eps;

[0012] ② Calculate the depth-domain generalized seismic wavelet for each parameter grid P(i,j)=[K(i),U(j)], and the synthesis formula for the depth-domain generalized seismic wavelet is as follows:

[0013]

[0014] Among them, wavelet Generate For the depth domain generalized seismic wavelet model, k and h are the wave number and depth of the depth domain generalized seismic wavelet, respectively, and i is the imaginary unit;

[0015] ③ Use the depth-domain generalized seismic wavelet under each parameter grid P(i,j) constructed in the previous step to synthesize the depth-domain seismic record S″, as shown in the following formula:

[0016]

[0017] Where S″(i,j,ii) is the ii-th sampling point of the depth domain seismic record synthesized under the parameter grid P(i,j);

[0018] ④ Based on the calculation results of the previous step, calculate the correlation coefficient matrix C between the synthesized depth-domain seismic record S″ and the resampled actual well-side seismic record S′. Select the parameter grid position with the maximum correlation coefficient in the correlation coefficient matrix C to obtain the optimal depth-domain generalized seismic wavelet reference parameter K in the iter-th iteration process. iter Sum of order U iter ;

[0019] ⑤ Based on the optimal synthetic seismic record S from the previous step b " est Calculate the optimal seismic wavelet amplitude A in each iteration. iter The formula is as follows:

[0020]

[0021] ⑥ The formula for calculating residuals is as follows:

[0022] S′=S′-A iter ×S b " est ;

[0023] ⑦ Repeat steps ② through ⑥ until the required number of iterations is reached or the following exit condition is met:

[0024]

[0025] in, Represents the 2-norm;

[0026] (7) Based on the calculation results of the previous step, the seismic wavelets in each iteration are superimposed to obtain the final seismic wavelet, as shown in the following formula:

[0027]

[0028] in, This is the final extracted depth-domain seismic wavelet. Attached Figure Description

[0029] Figure 1 This refers to the depth-domain logging data used in this embodiment of the invention. The logging data has a depth range of 3.05 km to 3.27 km, a sampling interval of 0.1 m, and a total of 2201 sampling points. Figure 1 (a) represents the P-wave velocity in well logging. The vertical axis represents depth in kilometers, and the horizontal axis represents velocity in kilometers per second (km / s). Figure 1 (b) represents the logging density, with the vertical axis representing depth in kilometers and the horizontal axis representing density in grams per cubic centimeter (g / cm³). 3 ). Figure 1(c) P-wave impedance in well logging, with the vertical axis representing depth in kilometers and the horizontal axis representing P-wave impedance in grams per cubic centimeter per second (km / s·g / cm). 3 ).

[0030] Figure 2 These are the well logging reflection coefficients calculated in this embodiment of the invention and the resampled wellside seismic data. Figure 1 (a) The well logging reflection coefficient depth is 3.05km-3.27km, the sampling interval is 0.1m, the number of sampling points is 2201, the vertical axis is the depth in kilometers, and the horizontal axis is the amplitude. Figure 1 (b) is the resampled depth domain well-side seismic record. The original depth of this seismic data is 3.05km-3.27km, the original sampling interval is 5m, and the number of sampling points is 45. After resampling, the depth is 3.05km-3.27km, the sampling interval is 0.1m, and the number of sampling points is 2201. The vertical axis is the depth in kilometers, and the horizontal axis is the trace number. This well is located at trace number 177.

[0031] Figure 3 This is the final seismic wavelet extraction result of the embodiment of the present invention. Figure 3 (a) Comparison between actual well-side depth domain seismic data and depth domain seismic data synthesized based on extracted depth domain seismic wavelets. The depth domain seismic data synthesized based on extracted depth domain seismic wavelets (grey waveform) is overlaid on the actual well-side depth domain seismic data. The vertical axis represents depth in kilometers, and the horizontal axis represents the well number, which is located at well number 177. Figure 3 (b) The black curve represents the final extracted depth domain seismic wavelet, which has 2001 sampling points with a sampling interval of 0.1m. The vertical axis represents amplitude and the horizontal axis represents length, with the unit being meters. Figure 3 (c) represents the amplitude of the depth-domain seismic wavelet extracted in each iteration, with the vertical axis representing the amplitude and the horizontal axis representing the iteration number. Figure 3 (d) The depth domain seismic wavelet extracted in each iteration, with the vertical axis representing the length in meters and the horizontal axis representing the iteration number. Detailed Implementation

[0032] Figure 1-3 The method for extracting seismic wavelets in the depth domain presented in this paper is demonstrated, and its specific implementation includes the following steps:

[0033] (1) The depth domain logging reflection coefficient r is calculated using the depth domain logging density ρ and the P-wave velocity Vp. The number of sampling points for the logging reflection coefficient is N1 = 2201, and the sampling interval in the logging depth domain is Δh1 = 0.1m. The formula for calculating the reflection coefficient is:

[0034]

[0035] Where r(i) is the depth domain reflection coefficient when the index number is i, ρ(i) is the depth domain logging density when the index number is i, and Vp(i) is the depth domain logging velocity when the index number is i.

[0036] (2) Extract the single-track seismic record S in the depth domain near the well. The number of sampling points for the single-track seismic record in the depth domain near the well is N2 = 45, and the depth sampling interval is Δh2 = 5m.

[0037] (3) Using the three-sample interpolation method with a well depth domain sampling interval of Δh1, the original seismic data is resampled to obtain the resampled seismic data S′.

[0038] (4) Set the length N3 = 2001, velocity Vc = 3000m / s, and sampling interval Δh of the depth domain seismic wavelet to be extracted. wave =0.1m;

[0039] (5) Set the reference wavenumber search range K and order range U for the depth domain generalized seismic wavelet to be extracted, and construct the parameter grid P(i,j)=[K(i),U(j)] for seismic wavelet extraction;

[0040] (6) The three parameters of the depth-domain generalized seismic wavelet, namely the reference wave number k0, order u, and amplitude a, are solved by the alternating relaxation iteration method as follows:

[0041] ① Set the maximum number of iterations L = 3, and the minimum iteration error eps = 1 × 10 -6 ;

[0042] ② Calculate the depth-domain generalized seismic wavelet for each parameter grid P(i,j)=[K(i),U(j)], and the synthesis formula for the depth-domain generalized seismic wavelet is as follows:

[0043]

[0044] Among them, wavelet Generate For the depth domain generalized seismic wavelet model, k and h are the wave number and depth of the depth domain generalized seismic wavelet, respectively, and i is the imaginary unit;

[0045] ③ Use the depth-domain generalized seismic wavelet under each parameter grid P(i,j) constructed in the previous step to synthesize the depth-domain seismic record S″, as shown in the following formula:

[0046]

[0047] Where S″(i,j,ii) is the ii-th sampling point of the depth domain seismic record synthesized under the parameter grid P(i,j);

[0048] ④ Based on the calculation results of the previous step, calculate the correlation coefficient matrix C between the synthesized depth-domain seismic record S″ and the resampled actual well-side seismic record S′. Select the parameter grid position with the maximum correlation coefficient in the correlation coefficient matrix C to obtain the optimal depth-domain generalized seismic wavelet reference parameter K in the iter-th iteration process. iter Sum of order U iter ;

[0049] ⑤ Based on the optimal synthetic seismic record S from the previous step b " est Calculate the optimal seismic wavelet amplitude A in each iteration. iter The formula is as follows:

[0050]

[0051] ⑥ The formula for calculating residuals is as follows:

[0052] S′=S′-A iter ×S b " est ;

[0053] ⑦ Repeat steps ② through ⑥ until the required number of iterations is reached or the following exit condition is met:

[0054]

[0055] in, Represents the 2-norm;

[0056] (7) Based on the calculation results of the previous step, the seismic wavelets in each iteration are superimposed to obtain the final seismic wavelet, as shown in the following formula:

[0057]

[0058] in, This is the final extracted depth-domain seismic wavelet.

[0059] Figure 1 For the purposes of this invention, depth-domain logging data is used, which includes density ( Figure 1 a) Longitudinal wave velocity ( Figure 1 b) and longitudinal wave impedance ( Figure 1 c) Curve, where the well logging data depth ranges from 3.05km to 3.27km, the well logging depth sampling interval is 0.1m, and the number of sampling points is 2201.

[0060] Figure 2The well logging reflection coefficients calculated for this embodiment of the invention and the resampled well-side seismic data both have a depth range of 3.05km-3.27km, a depth sampling interval of 0.1m, and 2201 sampling points. Figure 2 The blue curve in b represents the actual well-side seismic record. The main purpose of resampling the seismic data is to ensure the integrity of the well logging information without omission. Given that the original seismic data has a sampling interval of 5m in the depth domain, directly downsampling the well logging data will result in the loss of depth domain well logging information. This situation will seriously hinder the accuracy and reliability of high-precision extraction of seismic wavelets in the depth domain.

[0061] Figure 3 This is the final seismic wavelet extraction result of the embodiment of the present invention. Figure 3 This demonstrates a high degree of agreement between actual well-side depth-domain seismic data and depth-domain seismic data synthesized from extracted depth-domain seismic wavelets (grey waveform), with a correlation coefficient of 0.97. Figure 3 The black curve in figure b represents the final extracted depth-domain seismic wavelet, which consists of 2001 sampling points with a sampling interval of 0.1m. The method of this invention uses a consistent sampling interval throughout the depth of the domain logging to improve the accuracy of depth-domain seismic wavelet extraction. Figure 3 c represents the amplitude of the depth-domain seismic wavelet extracted in each iteration. Figure 3 d represents the depth domain seismic wavelet extracted in each iteration. It can be seen that the method of the present invention only requires 3 iterations to extract the high-precision depth domain seismic wavelet, indicating that the method of the present invention is effective and reliable.

[0062] The advantages of the method of the present invention are: (1) it uses data under the original well logging sampling interval without losing the effective information of the well logging data; (2) in the process of extracting the depth domain seismic wavelet, only three model parameters of the depth domain generalized seismic wavelet need to be determined, and their upper and lower limits are clear and easy to define; (3) in practical applications, it has high efficiency and strong stability, and only 3 to 5 iterations are needed to extract reliable depth domain seismic wavelets.

[0063] The above embodiments are only used to illustrate the present invention. The implementation steps of the method can be varied. Any equivalent transformations and improvements made on the basis of the technical solution of the present invention should not be excluded from the protection scope of the present invention.

Claims

1. A depth-domain seismic wavelet extraction method based on a generalized seismic wavelet model, comprising the following steps: (1) Using depth domain logging density and longitudinal wave velocity Calculate depth domain logging reflection coefficient The number of sampling points for well logging reflection coefficient is The sampling interval in the well logging depth domain is ; (2) Extraction of single-channel seismic records in the depth domain near the well. The number of sampling points for a single-channel seismic record in the depth domain near the well is The depth sampling interval is ; (3) Using a three-sample interpolation method, the sampling interval in the logging depth domain is utilized. The original seismic data is resampled to obtain resampled seismic data. ; (4) Set the length of the depth domain seismic wavelet to be extracted. ,speed and sampling interval ; (5) Set the search range of the reference wavenumber for the depth-domain generalized seismic wavelet to be extracted. and order range And construct a parametric grid for seismic wavelet extraction. ; (6) The reference wavenumber of the depth-domain generalized seismic wavelet is solved by the alternating relaxation iterative method. order and amplitude These three parameters are obtained as follows: ① Set the maximum number of iterations Minimum iteration error ; Calculate each parameter grid The depth-domain generalized seismic wavelet synthesis formula is as follows: ， in, For the depth domain generalized seismic wavelet model, These are the wavenumber and depth of the generalized seismic wavelet in the depth domain, respectively, where I is the imaginary unit; ③ Use each parameter mesh built in the previous step Depth Domain Generalized Seismic Wavelet Synthetic depth domain seismic records The formula is as follows: ， in, It is a parametric mesh The first layer of the depth domain seismic record synthesized below One sampling point; ④ Based on the calculation results from the previous step, calculate the synthetic depth-domain seismic record. Compared with the actual well-side seismic records after resampling Correlation coefficient matrix Select the correlation coefficient matrix The parameter grid position with the maximum correlation coefficient is obtained to obtain the first... The optimal depth-domain generalized seismic wavelet reference wavenumber in the next iteration process sum order ; ⑤ Based on the optimal synthetic seismic record from the previous step Calculate the optimal seismic wavelet amplitude in each iteration. The formula is as follows: ； ⑥ Calculate the residual using the following formula: ； ⑦ Repeat Continue until step 6 is reached, or the following exit condition is met: ， in, Represents the 2-norm; (7) Based on the calculation results of the previous step, the seismic wavelets in each iteration are superimposed to obtain the final seismic wavelet, as shown in the following formula: ， in, This is the final extracted depth-domain seismic wavelet.

Images