Depth domain seismic wavelet extraction method based on generalized seismic wavelet model

Abstract

The invention discloses a depth domain seismic wavelet extraction method based on a generalized seismic wavelet model. The method is characterized in that depth domain logging data information can be fully utilized, depth domain seismic records are directly synthesized in a depth domain, depth domain seismic wavelets can be efficiently and accurately extracted only through 3-5 times of iteration from a depth domain generalized seismic wavelet model, the method has high stability, and the method is suitable for large-scale popularization and application. And the problem of effective information loss caused by domain conversion or data resampling is avoided.

Description

Technical Field

[0001] The invention belongs to the field of oil and gas seismic exploration, and in particular relates to a depth domain seismic wavelet extraction method based on a generalized seismic wavelet model. Background Art

[0002] In the field of oil and gas seismic exploration, with the rapid development and widespread application of pre-stack depth domain migration imaging, it is an urgent practical need to directly realize reservoir parameter inversion, oil and gas reservoir prediction and description in the depth domain. Extracting seismic wavelets from depth domain seismic data is an indispensable and important link and basis for realizing deep domain seismic inversion. As a basic component of seismic data, the characteristics of seismic wavelets are of great significance for revealing underground structures. In the process of seismic data processing, the analysis and application of seismic wavelets run through multiple key links, including but not limited to deconvolution, well seismic calibration, forward modeling and inversion of seismic records, etc. However, due to the non-stationarity of deep domain seismic data, the traditional convolution model method based on the assumption of linear deep invariance system cannot be directly applied to the depth domain, which brings great challenges to the extraction of deep domain seismic wavelets.

[0003] In order to meet this challenge, researchers are constantly exploring new methods and technologies in order to accurately and efficiently extract deep-domain seismic wavelets, so as to better serve the needs of oil and gas exploration and development. At present, there are two main methods for extracting deep-domain seismic wavelets: (1) converting deep-domain seismic data into the time domain or constant-velocity depth domain, extracting seismic wavelets and then converting them into the depth domain. This domain conversion method brings inevitable information loss and cumulative errors to deep-domain seismic data and deep-domain well logging data; (2) using depth-wavenumber spectrum analysis technology to directly extract deep-domain seismic wavelets at each depth position, but this method lacks the constraints of deep-domain well logging data, resulting in multiple solutions for the extraction results, and the seismic data needs to be processed one by one, which is time-consuming and inefficient. Since deep-domain seismic data does not have the "linear depth invariant" characteristic, it is impossible to directly apply the convolution model type method to synthesize deep-domain seismic records. Therefore, the present invention is to better solve this problem, and can directly extract deep-domain seismic wavelets in the depth domain, quickly and efficiently, and provide a basic method for deep-domain reservoir parameter inversion and fluid identification. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a method for extracting seismic wavelets in the depth domain based on a generalized seismic wavelet model, which effectively solves the problems of unreliable seismic wavelets extracted by domain conversion and time-consuming and multi-solution-prone seismic wavelet extraction using depth-wavenumber spectrum. The method of the present invention includes the following main steps:

[0005] (1) The depth domain logging density ρ and the P-wave velocity Vp are used to calculate the depth domain logging reflection coefficient r. The number of sampling points of the logging reflection coefficient is N1, and the sampling interval of the logging depth domain is Δh1;

[0006] (2) Extracting the single-channel seismic record S in the depth domain near the well, the number of sampling points of the single-channel seismic record in the depth domain near the well is N2, and the depth sampling interval is Δh2;

[0007] (3) Using the three-sample interpolation method, the sampling interval in the logging depth domain is Δh1, and 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) Setting the reference wavenumber search range K and order range U of the depth domain generalized seismic wavelet to be extracted, and constructing a parameter grid P(i, j) = [K(i), U(j)] for seismic wavelet extraction;

[0010] (6) The three parameters of the reference wave number k0, order u and amplitude a of the depth domain generalized seismic wavelet are solved by the alternating relaxation iteration 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 under each parameter grid P(i,j)=[K(i),U(j)]. The depth domain generalized seismic wavelet synthesis formula is as follows:

[0013]

[0014] Among them, wavelet Generate is 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 an 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″, the formula is as follows:

[0016]

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

[0018] ④ According to the calculation results of the previous step, calculate the correlation coefficient matrix C between the synthetic depth domain seismic record S″ and the actual wellside seismic record S′ after resampling, select the parameter grid position with the maximum correlation coefficient in the correlation coefficient matrix C, and obtain the optimal depth domain generalized seismic wavelet reference parameter K in the iter-th iteration process iter and order U iter ;

[0019] ⑤ Based on the optimal synthetic seismic record S in the previous step b ″ est , calculate the optimal seismic wavelet amplitude A in each iteration process iter , the formula is as follows:

[0020]

[0021] ⑥ Calculate the residual formula, the formula is as follows:

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

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

[0024]

[0025] in, represents the 2-norm;

[0026] (7) According to the calculation results of the previous step, the seismic wavelets in each iteration process are superimposed to obtain the final seismic wavelet. The formula is as follows:

[0027]

[0028] in, is the final extracted depth domain seismic wavelet. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 is the depth domain logging data used in the embodiment of the present invention. The depth range of the logging data is 3.05km-3.27km, the sampling interval of the logging depth is 0.1m, and the number of sampling points is 2201. Among them, Figure 1 (a) is the logging P-wave velocity, the ordinate is depth in kilometers, and the abscissa is velocity in kilometers per meter (km / s). Figure 1 (b) is the well logging density, the ordinate is the depth in kilometers, and the abscissa is the density in grams per cubic centimeter (gcm 3 ). Figure 1(c) is the logging P-wave impedance, the ordinate is depth in kilometers, and the abscissa is P-wave impedance in grams per cubic centimeter per second (km / s·g / cm 3 ).

[0030] Figure 2 It is the well logging reflection coefficient calculated by the embodiment of the present invention and the resampled wellside seismic data. Figure 1 (a) is the logging reflection coefficient with a depth of 3.05km-3.27km, a sampling interval of 0.1m, and 2201 sampling points. The ordinate is depth in kilometers, and the abscissa is amplitude. Figure 1 (b) is the depth domain wellside seismic record after resampling. The original depth of the 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 ordinate is depth in kilometers, and the abscissa is the track number. The well is located at the 177th track.

[0031] Figure 3 This is the final seismic wavelet extraction result of the embodiment of the present invention. Figure 3 (a) is a comparison between the actual wellbore depth-domain seismic data and the depth-domain seismic data synthesized based on the extracted depth-domain seismic wavelet. The depth-domain seismic data synthesized based on the extracted depth-domain seismic wavelet (gray waveform) is superimposed on the actual wellbore depth-domain seismic data. The ordinate is depth in kilometers, and the abscissa is the track number. The well is located at the 177th track. Figure 3 (b) The black curve is the final extracted depth-domain seismic wavelet, which has 2001 sampling points with a sampling interval of 0.1 m. The ordinate is amplitude and the abscissa is length in meters. Figure 3 (c) is the amplitude of the seismic wavelet in the depth domain extracted during each iteration. The ordinate is the amplitude and the abscissa is the number of iterations. Figure 3 (d) The depth-domain seismic wavelet extracted during each iteration. The ordinate is the length in meters, and the abscissa is the number of iterations. DETAILED DESCRIPTION

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

[0033] (1) The depth domain logging density ρ and the P-wave velocity Vp are used to calculate the depth domain logging reflection coefficient r. The number of sampling points for the logging reflection coefficient is N1 = 2201, and the sampling interval in the depth domain of the logging is Δh1 = 0.1 m. The calculation formula for 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) Extracting the single-channel seismic record S in the depth domain near the well, the number of sampling points of the single-channel 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, the sampling interval in the logging depth domain is Δh1, and the original seismic data is resampled to obtain the resampled seismic data S′;

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

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

[0040] (6) The three parameters of the reference wave number k0, order u and amplitude a of the depth domain generalized seismic wavelet 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 under each parameter grid P(i,j)=[K(i),U(j)]. The depth domain generalized seismic wavelet synthesis formula is as follows:

[0043]

[0044] Among them, wavelet Generate is 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 an 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″, the formula is as follows:

[0046]

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

[0048] ④ According to the calculation results of the previous step, calculate the correlation coefficient matrix C between the synthetic depth domain seismic record S″ and the actual wellside seismic record S′ after resampling, select the parameter grid position with the maximum correlation coefficient in the correlation coefficient matrix C, and obtain the optimal depth domain generalized seismic wavelet reference parameter K in the iter-th iteration process iter and order U iter ;

[0049] ⑤ Based on the optimal synthetic seismic record S in the previous step b ″ est , calculate the optimal seismic wavelet amplitude A in each iteration process iter , the formula is as follows:

[0050]

[0051] ⑥ Calculate the residual formula, the formula is as follows:

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

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

[0054]

[0055] in, represents the 2-norm;

[0056] (7) According to the calculation results of the previous step, the seismic wavelets in each iteration process are superimposed to obtain the final seismic wavelet. The formula is as follows:

[0057]

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

[0059] Figure 1 The present invention uses depth domain logging data, which includes density ( Figure 1 a) Longitudinal wave velocity ( Figure 1 b) and longitudinal wave impedance ( Figure 1 c) Curve, where the logging data depth range is 3.05km-3.27km, the logging depth sampling interval is 0.1m, and the number of sampling points is 2201 sampling points.

[0060] Figure 2The well logging reflection coefficient calculated in the embodiment of the present invention and the resampled wellside seismic data have a depth range of 3.05km-3.27km, a depth sampling interval of 0.1m, and a number of sampling points of 2201. Figure 2 The blue curve in b is the actual wellside seismic record. The main purpose of resampling the seismic data is to ensure the integrity of the logging information. Given that the sampling interval of the original seismic data in the depth domain is 5 m, if the logging data is directly downsampled, the depth domain logging information will be lost. This situation will seriously hinder the accuracy and reliability of high-precision extraction of seismic wavelets in the depth domain.

[0061] Figure 3 is the final seismic wavelet extraction result of the embodiment of the present invention, Figure 3 a shows the high consistency between the actual well-side depth domain seismic data and the depth domain seismic data (gray waveform) synthesized based on the extracted depth domain seismic wavelet, and the correlation coefficient between the two is 0.97. Figure 3 The black curve in b is the final extracted depth domain seismic wavelet, which has 2001 sampling points and a sampling interval of 0.1 m. The method of the present invention uses a sampling interval consistent with the well logging depth to improve the accuracy of depth domain seismic wavelet extraction. Figure 3 c is the amplitude of the depth domain seismic wavelet extracted in each iteration process, Figure 3 d is the depth domain seismic wavelet extracted in each iteration process. It can be seen that the method of the present invention only needs 3 iterations to extract high-precision depth domain seismic wavelets, which shows that the method of the present invention is effective and reliable.

[0062] The advantages of the method of the present invention are: (1) the data at the original logging sampling interval is used without losing the effective information of the 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 application, it has high efficiency and strong stability, and only 3 to 5 iterations are required to extract reliable depth domain seismic wavelets.

[0063] The above embodiments are only used to illustrate the present invention, wherein the various implementation steps of the method can be changed, and any equivalent transformations and improvements based on the technical solution of the present invention should not be excluded from the protection scope of the present invention.

Claims

1. A method for extracting seismic wavelets in depth domain based on a generalized seismic wavelet model, comprising the following main steps: (1) The depth domain logging density ρ and the P-wave velocity Vp are used to calculate the depth domain logging reflection coefficient r. The number of sampling points of the logging reflection coefficient is N1, and the sampling interval of the logging depth domain is Δh1; (2) Extracting the single-channel seismic record S in the depth domain near the well, the number of sampling points of the single-channel seismic record in the depth domain near the well is N2, and the depth sampling interval is Δh2; (3) Using the three-sample interpolation method, the sampling interval in the logging depth domain is Δh1, and the original seismic data is resampled to obtain the resampled seismic data S′; (4) Set the length N3, velocity Vc and sampling interval Δh of the depth domain seismic wavelet to be extracted wave ; (5) Setting the reference wavenumber search range K and order range U of the depth domain generalized seismic wavelet to be extracted, and constructing a parameter grid P(i, j) = [K(i), U(j)] for seismic wavelet extraction; (6) The three parameters of the reference wave number k0, order u and amplitude a of the depth domain generalized seismic wavelet are solved by the alternating relaxation iteration method as follows: ① Set the maximum number of iterations L and the minimum iteration error eps; ② Calculate the depth domain generalized seismic wavelet under each parameter grid P(i,j)=[K(i),U(j)]. The synthesis formula of the depth domain generalized seismic wavelet is as follows: in, wavelet Generate is 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 an imaginary unit; ③ 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″, the formula is as follows: Where S″(i,j,ii) is the ii-th sampling point of the synthesized depth domain seismic record under the parameter grid P(i,j); ④ According to the calculation results of the previous step, calculate the correlation coefficient matrix C between the synthetic depth domain seismic record S″ and the actual wellside seismic record S′ after resampling, select the parameter grid position with the maximum correlation coefficient in the correlation coefficient matrix C, and obtain the optimal depth domain generalized seismic wavelet reference parameter K in the iter-th iteration process iter and order U iter ; ⑤ Based on the optimal synthetic seismic record S in the previous step b ″ est , calculate the optimal seismic wavelet amplitude A in each iteration process iter , the formula is as follows: ⑥ Calculate the residual formula, the formula is as follows: S′=S′-A iter ×S b ″ est ; ⑦ Repeat ② to ⑥ until the number of iterations is reached or the following exit condition is met: in, represents the 2-norm; (7) According to the calculation results of the previous step, the seismic wavelets in each iteration process are superimposed to obtain the final seismic wavelet. The formula is as follows: in, is the final extracted seismic wavelet in the depth domain.

Images