Method for extracting depth domain seismic wavelets through iterative optimization

Abstract

The invention discloses a method for extracting depth domain seismic wavelets through iterative optimization. Depth domain data and original logging data are directly utilized to quickly and efficiently extract the depth domain seismic wavelets. According to the method, the depth domain seismic wavelet basis is established, the objective function for extracting the depth domain seismic wavelet is constructed, and the control parameters of the depth domain seismic wavelet basis are obtained through optimization iteration, so that the final depth domain seismic wavelet is obtained, the parameter adjustment process is simplified, and the operation difficulty and subjectivity are reduced. In practical application, depth domain seismic wavelets can be accurately and reliably extracted through iteration of limited times, and high robustness and accuracy can still be kept under complex geological conditions.

Description

Technical Field

[0001] The invention belongs to the field of oil and gas seismic exploration, and in particular relates to a method for iteratively optimizing and extracting seismic wavelets in a depth domain. Background Art

[0002] Extracting seismic wavelets from deep-domain seismic data is an indispensable and important cornerstone for achieving deep-domain seismic inversion. With the widespread application of deep-domain seismic migration data in earth exploration and oil and gas resource exploration, inversion using deep-domain seismic data is an important basic means to achieve fine imaging of underground stratigraphic structure and important structural information, characterize reservoir distribution, and identify oil and gas. However, since deep-domain seismic data does not meet the "linear depth invariance" condition of the traditional convolution model, the existing time-domain seismic wavelet extraction method cannot be directly applied to the depth domain, and the work of directly extracting seismic wavelets in the depth domain faces great challenges.

[0003] At present, the extraction of seismic wavelets in the depth domain mainly relies on two technical implementation paths: one method is to convert the depth domain seismic data to the time domain or constant velocity depth domain, use the existing time domain method to extract the seismic wavelet, and then convert the extracted seismic wavelet to the depth domain (Zhang et al., 2018; 2019; 2020; 2022). However, this domain conversion process inevitably leads to information loss and error accumulation of depth domain seismic data and depth domain well logging data, which affects the accuracy of seismic wavelet extraction and reduces the computational efficiency. Another method is to use depth-wavenumber spectrum analysis technology to extract a separate seismic wavelet for each seismic trace (Zhang and Deng, 2018; 2023; Andreas et al., 2019; Sengupta et al., 2021; Cai et al., 2022; 2024; Sun et al., 2025). al., 2024), due to the lack of constraints of logging data in the depth domain, the seismic wavelet extraction results have multiple solutions, and because the seismic data needs to be processed one by one, the calculation cost is large and the efficiency is low. The present invention is to solve the limitations and problems of these two methods and directly realize the seismic wavelet extraction with higher efficiency and accuracy in the depth domain. Summary of the invention

[0004] The purpose of the method of the present invention is to provide a method for quickly and efficiently extracting deep-domain seismic wavelets by directly using deep-domain data and original well logging data. The method comprises the following main steps:

[0005] (1) Input the depth domain logging data of the well, and use the depth domain logging P-wave velocity and V density ρ to calculate the depth domain logging reflection coefficient r; input the actual depth domain seismic record S near the well, the length of the seismic record is N, and the depth sampling interval of the seismic record is Δh; according to the P-wave velocity V, set the depth domain seismic wavelet velocity V to be extracted c ;

[0006] (2) Using the characteristics of the actual depth domain seismic record S near the well and the depth logging reflection coefficient r of the well, the depth domain seismic wavelet basis ψ(h|k m ,σ,φ,a), construct the objective function for deep domain seismic wavelet extraction:

[0007]

[0008] Where h represents the depth, represents the L2 norm, Rs j represents the residual of the jth iteration and sets Rs 1 =S,(k m ,σ,φ,a) are the control parameter sets of the seismic wavelet basis in the depth domain, which represent the average wave number, scale factor, phase and amplitude respectively;

[0009] In the above formula, the wavenumber domain seismic wavelet calculated by the depth domain seismic wavelet is used to synthesize the depth domain seismic record S′, and the formula is as follows:

[0010] S′(k m ,σ,φ,a)=real{B[A⊙Φ(k|k m ,σ,φ,a)]}·r,

[0011] Where ⊙ represents the Hadamard product operation of the vector, real{·} represents the real part operation, k is the wave number, Φ is the seismic wavelet basis in the wave number domain, A and B are the established phase extension matrix and its corresponding conjugate matrix respectively;

[0012] (3) Preset the number of iterations and the minimum residual E, and use the objective function of the previous step to iteratively solve the control parameter set of the j-th depth-domain seismic wavelet basis. The specific process is as follows:

[0013] ① Calculate the average wave number and scale factor of the seismic wavelet basis in the depth domain:

[0014]

[0015] Where H(·) represents the Hilbert transform operation, and <·,·> represents the vector inner product operation;

[0016] ② Calculate the amplitude coefficient P in the complex domain j :

[0017]

[0018] in, The superscript T indicates the matrix transpose operation;

[0019] ③ Calculate the phase and amplitude of the seismic wavelet basis in the depth domain:

[0020]

[0021] Among them, Angle(·) is the argument operation of complex numbers, |·| is the modular operation of complex numbers;

[0022] ④ Calculate the residual according to the following formula:

[0023]

[0024] ⑤ Repeat ① to ④ until the preset number of iterations is reached or the following formula is satisfied, then exit the iteration process:

[0025]

[0026] (4) According to the calculation in the previous step, the depth domain seismic wavelet basis extracted in each iteration is superimposed to obtain the final depth domain seismic wavelet:

[0027]

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

[0029] Figure 1 The depth domain well logging data and wellside seismic data used in the embodiment of the present invention have a depth range of 2914m-3133m, a sampling interval of 0.15m, and 1461 sampling points. Figure 1 (a) is the logging P-wave velocity, the ordinate is depth in meters, and the abscissa is velocity in meters per second (m / s). Figure 1 (b) is the well logging density, the ordinate is the depth in meters, and the abscissa is the density in grams per cubic centimeter (g / cm 3 ). Figure 1 (c) is the logging reflection coefficient, the ordinate is depth, and the unit is meter. Figure 1 (d) is the wellside seismic record in the depth domain. The ordinate is the depth in meters, and the abscissa is the track number. The well is located in the 6th track (marked with an inverted triangle).

[0030] Figure 2 The depth domain seismic wavelet extraction process and results of the embodiment of the present invention. Figure 2(a) is a comparison between the actual wellbore depth-domain seismic data (black waveform) and the depth-domain seismic data (gray waveform) reconstructed using the depth-domain seismic wavelet extracted by the present invention. Note: To facilitate direct observation and comparison, the depth-domain seismic data (gray waveform) reconstructed using the depth-domain seismic wavelet extracted by the present invention is superimposed on the actual wellbore depth-domain seismic data (black waveform) profile. The ordinate is depth in meters, and the abscissa is the track number. Figure 2 (b) is the final extracted seismic wavelet in the depth domain, where the ordinate is amplitude and the abscissa is length in meters. Figure 2 (c) is the amplitude of the seismic wavelet basis in the depth domain extracted during each iteration. The ordinate is the amplitude and the abscissa is the number of iterations. Figure 2 (d) is the seismic wavelet basis in the depth domain extracted during each iteration. The ordinate is the length in meters, and the abscissa is the number of iterations. DETAILED DESCRIPTION

[0031] Figure 1 and Figure 2 The depth domain seismic wavelet extraction method of the present invention is demonstrated, and its corresponding specific implementation method includes the following steps:

[0032] (1) Input the depth domain logging data of the well (such as Figure 1 (a) and Figure 1 (b)) using the depth domain logging P-wave velocity ( Figure 1 (a)) and V density ρ( Figure 1 (b)) to calculate the depth domain logging reflection coefficient r( Figure 1 (c)); Input the actual depth domain seismic record S( Figure 1 (d)), the well is located in the 6th track (inverted triangle mark), the length of the seismic record is 1461, and the depth sampling interval of the seismic record is 0.15m; according to the longitudinal wave velocity V, the depth domain seismic wavelet velocity V to be extracted is set c =3000m / s;

[0033] (2) Using the characteristics of the actual depth domain seismic record S near the well and the depth logging reflection coefficient r of the well, the depth domain seismic wavelet basis ψ(h|k m ,σ,φ,a), construct the objective function for deep domain seismic wavelet extraction:

[0034]

[0035] Where h represents the depth, represents the L2 norm, Rs j represents the residual of the jth iteration and sets Rs 1 =S,(k m,σ,φ,a) are the control parameter sets of the seismic wavelet basis in the depth domain, which represent the average wave number, scale factor, phase and amplitude respectively;

[0036] In the above formula, the wavenumber domain seismic wavelet calculated by the depth domain seismic wavelet is used to synthesize the depth domain seismic record S′, and the formula is as follows:

[0037] S′(k m ,σ,φ,a)=real{B[A⊙Φ(k|k m ,σ,φ,a)]}·r,

[0038] Where ⊙ represents the Hadamard product operation of the vector, real{·} represents the real part operation, k is the wave number, Φ is the seismic wavelet basis in the wave number domain, A and B are the established phase extension matrix and its corresponding conjugate matrix respectively;

[0039] (3) The number of iterations is set to 5 and the minimum residual E = 0.001. The control parameter set of the j-th depth domain seismic wavelet basis is iteratively solved using the objective function of the previous step. The specific process is as follows:

[0040] ① Calculate the average wave number and scale factor of the seismic wavelet basis in the depth domain:

[0041]

[0042] Where H(·) represents the Hilbert transform operation, and <·,·> represents the vector inner product operation;

[0043] ② Calculate the amplitude coefficient P in the complex domain j :

[0044]

[0045] in, The superscript T indicates the matrix transpose operation;

[0046] ③ Calculate the phase and amplitude of the seismic wavelet basis in the depth domain:

[0047]

[0048] Among them, Angle(·) is the argument operation of complex numbers, |·| is the modular operation of complex numbers;

[0049] ④ Calculate the residual according to the following formula:

[0050]

[0051] ⑤ Repeat ① to ④ until the preset number of iterations is reached or the following formula is satisfied, then exit the iteration process:

[0052]

[0053] (6) According to the calculation in the previous step, the actual iteration is completed after 2 times ( Figure 2 (c), i.e., j = 2), the depth domain seismic wavelet basis extracted by 2 iterations ( Figure 2 (d)) are stacked to obtain the final depth domain seismic wavelet ( Figure 2 (b)):

[0054]

[0055] exist Figure 2 In (a), the depth domain seismic data (gray waveform) reconstructed using the final depth domain seismic wavelet extracted in an embodiment of the present invention is highly consistent with the actual depth domain seismic data (black waveform) near the well, and the correlation coefficient between the two reaches 0.96, indicating that the method of the present invention is efficient and reliable.

[0056] The advantages of the method of the present invention are: (1) the method of the present invention directly uses the original logging data to extract the deep domain seismic wavelet without resampling the data, thereby avoiding the problem of missing detailed information caused by resampling, and ensuring that the effective information of the logging data is fully utilized; (2) in the process of extracting the deep domain seismic wavelet, the method of the present invention only needs to take four control parameters, namely, the average wave number, the scale factor, the phase and the amplitude. Compared with the traditional seismic wavelet extraction method, the parameter adjustment process is greatly simplified, and the operation difficulty and subjectivity are reduced; (3) in actual application, the method of the present invention only needs 2-3 iterative calculations to accurately and reliably extract the deep domain seismic wavelet, showing strong robustness and high computational efficiency.

[0057] 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 iteratively optimizing and extracting seismic wavelets in depth domain, comprising the following main steps: (1) Input the depth domain logging data of the well, and use the depth domain logging P-wave velocity V and density ρ to calculate the depth domain logging reflection coefficient r; input the actual depth domain seismic record S near the well, the length of the seismic record is N, and the depth sampling interval of the seismic record is Δh; according to the P-wave velocity V, set the depth domain seismic wavelet velocity V to be extracted c ; (2) Using the characteristics of the actual depth domain seismic record S near the well and the depth logging reflection coefficient r of the well, the depth domain seismic wavelet basis ψ(h|k m ,σ,φ,a), construct the objective function for deep domain seismic wavelet extraction: in, h represents the depth, represents the L2 norm, Rs j represents the residual of the jth iteration and sets Rs 1 =S,(k m ,σ,φ,a) are the control parameter sets of the seismic wavelet basis in the depth domain, which represent the average wave number, scale factor, phase and amplitude respectively; In the above formula, the wavenumber domain seismic wavelet calculated by the depth domain seismic wavelet is used to synthesize the depth domain seismic record S′, and the formula is as follows: S′(k m ,σ,φ,a)=real{B[A⊙Φ(k|k m ,σ,φ,a)]}·r, Where ⊙ represents the Hadamard product operation of the vector, real{·} represents the real part operation, k is the wave number, Φ is the seismic wavelet basis in the wave number domain, A and B are the established phase extension matrix and its corresponding conjugate matrix respectively; (3) Preset the number of iterations and the minimum residual E, and use the objective function of the previous step to iteratively solve the control parameter set of the j-th depth-domain seismic wavelet basis. The specific process is as follows: ① Calculate the average wave number and scale factor of the seismic wavelet basis in the depth domain: Where H(·) represents the Hilbert transform operation, and <·,·> represents the vector inner product operation; ② Calculate the amplitude coefficient P in the complex domain j : in, The superscript T indicates the matrix transpose operation; ③ Calculate the phase and amplitude of the seismic wavelet basis in the depth domain: Among them, Angle(·) is the argument operation of complex numbers, |·| is the modular operation of complex numbers; ④ Calculate the residual according to the following formula: ⑤ Repeat ① to ④ until the preset number of iterations is reached or the following formula is satisfied, then exit the iteration process: (4) According to the calculation in the previous step, the depth domain seismic wavelet basis extracted in each iteration is superimposed to obtain the final depth domain seismic wavelet: in, is the final depth-domain seismic wavelet.

Images