A fast extraction method of mixed phase seismic wavelet

Abstract

This invention discloses a rapid method for extracting mixed-phase seismic wavelets. The method employs a high-precision iterative complex domain basis pursuit method to extract the main feature information of seismic data. It utilizes the analytical WVD time-frequency distribution to efficiently calculate the instantaneous spectrum of the seismic data, rapidly obtaining the frequency and phase information for the final extraction of the seismic wavelet. In the absence of well logging data, this invention can accurately and efficiently extract seismic wavelets containing mixed-phase information from limited seismic data, effectively avoiding uncertainties caused by insufficient data and providing a crucial foundation for subsequent seismic data processing.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas seismic exploration and relates to a rapid extraction method for mixed-phase seismic wavelets. Background Technology

[0002] In seismic exploration for oil and gas, seismic wavelet extraction is a crucial foundation for seismic data processing. It is vital in processing, applied to core steps such as deconvolution, well-seismic calibration, and forward and inverse modeling of seismic records. The seismic wavelet is the direct response of the source signal propagating through the subsurface medium, accurately revealing the waveform shaping effect of seismic waves from the source to the receiver. The main purpose of seismic wavelet extraction is to extract wavelet characteristic information from complex seismic records, obtain key parameters such as the velocity and attenuation characteristics of seismic waves propagating underground, and provide data and theoretical support for seismic data interpretation and oil and gas exploration.

[0003] In the seismic wavelet extraction process, well logging data plays a crucial role and is a key element in ensuring the accuracy of seismic wavelet extraction. For areas lacking well logging data, the task of accurately extracting seismic wavelets is particularly challenging. We need to utilize the statistical information of seismic traces for wavelet extraction to ensure the reliability of the results. However, statistical seismic wavelet extraction methods based on the statistical information of the seismic traces themselves often yield wavelets that exhibit zero phase, minimum phase, or maximum phase, making it difficult to directly obtain the true phase information. Furthermore, seismic signals are often accompanied by significant noise and interference, which further exacerbates the deviation between the extracted wavelet and the actual wavelet, thus increasing the ambiguity problem of statistical seismic wavelet extraction methods. This situation poses a challenge in practical applications, making it difficult to determine the optimal seismic wavelet extraction scheme. Summary of the Invention

[0004] To address the shortcomings of existing technologies, the present invention aims to provide a rapid extraction method for mixed-phase seismic wavelets. The method is characterized by utilizing a high-precision iterative complex domain basis pursuit method and WVD time-frequency analysis to directly and accurately extract the average frequency and phase information of the seismic wavelet. This not only rapidly addresses the complexities of actual seismic signals but also possesses high accuracy and stability. The method of the present invention includes the following main steps:

[0005] (1) Input the seismic data S used for seismic wavelet extraction;

[0006] (2) The selected seismic data S is decomposed into Morlet wavelet one channel at a time using a high-precision iterative complex domain basis tracing method, as follows:

[0007] ① Calculate the analytical signal AS of the seismic data S;

[0008] ② Calculate the residual SR of the m-th seismic record at the k-th iteration using the following formula. k (m):

[0009]

[0010] Where AI is the analytic Morlet wavelet function, HT(·) denotes the Hilbert transform, and a,t,f,u,σ, These are the amplitude, time, average frequency, time shift, scale factor, and phase of the analytical Morlet wavelet, respectively, where i is the imaginary unit;

[0011] ③ Calculate the residual SR using the following formula k (m) The maximum position of the amplitude envelope, which serves as the time shift u of the analytic Morlet wavelet in the k-th decomposition. k (m):

[0012]

[0013] Where Env(·) represents calculating SR k The amplitude envelope of (m);

[0014] ④ Using u k (m), the analytical Morlet wavelet parameters [σ] of the k-th decomposition of the m-th seismic record are calculated using the following formula. k (m),f k (m)]:

[0015]

[0016] Where <·> is the dot product operator, and |·| is the absolute value symbol. The symbol for the 2-norm;

[0017] ⑤ Calculate the analytical Morlet wavelet amplitude A of the m-th seismic record in the k-th iteration using an optimization algorithm. k (m):

[0018]

[0019] Where ε is the damping factor and I is the identity matrix of dimension N;

[0020] ⑥Use A k (m), calculate the phase of the analytical Morlet wavelet using the following formula. and amplitude a k (m):

[0021]

[0022] Among them, Im[·] is the imaginary part operation of the orientation quantity, Re[·] is the real part operation of the orientation quantity, and abs[·] is the absolute value operation of the orientation quantity;

[0023] ⑦ Repeat steps ②-⑥ above until the required number of iterations is reached;

[0024] (3) Based on the decomposition results in the previous step, calculate the WVD time-frequency distribution WVD(m,t,f) of the m-th channel;

[0025] (4) Calculate the average frequency of the seismic wavelet based on the WVD time-frequency distribution results of the multichannel seismic data from the previous step. The method is as follows:

[0026] The WVD time-frequency distribution is normalized by channel, as shown in the following formula:

[0027]

[0028] Where max[·] represents taking the maximum value;

[0029] use The superimposed amplitude spectrum Amp(f) is calculated using the following formula:

[0030]

[0031] Based on Amp(f), the frequency corresponding to the maximum amplitude of the superimposed amplitude spectrum is obtained as the average frequency of the seismic wavelet to be extracted, as follows:

[0032]

[0033] (5) Based on the high-precision iterative complex domain basis tracing decomposition results of the seismic records, the phase parameters of the analytical Morlet wavelet extracted channel by channel are analyzed. Perform phase folding:

[0034]

[0035] Where φ1 is the folded phase;

[0036] (6) The phase information of the seismic wavelet is obtained using the following method:

[0037] Obtain the absolute value of the Morlet wavelet phase φ2:

[0038] φ2(m,k)=abs[φ1(m,k)];

[0039] The probability density function S of φ2 is calculated using the following formula:

[0040]

[0041] Where f(·) represents the weight of φ2.

[0042] Calculate the average phase corresponding to the maximum value in the probability density function. The formula is as follows:

[0043]

[0044] Where j is the index value of the probability density function S, and mean[·] represents the average value;

[0045] (7) Determine the order of the generalized seismic wavelet based on the phase formula of the generalized seismic wavelet.

[0046]

[0047] (8) According to and Calculating the reference frequency of the generalized seismic wavelet The parameters and formula are as follows:

[0048]

[0049] Where fr is the reference frequency of the generalized seismic wavelet, and Γ(·) is the gamma function;

[0050] (9) According to and The mixed-phase seismic wavelet w(t) is obtained by reconstruction using the following formula:

[0051] Attached Figure Description

[0052] Figure 1 This is the seismic wavelet result with a 45° phase extracted from the single-channel theoretical seismic record used in the embodiments of the present invention. Figure 1 (a) is a single-channel synthetic seismic record with 201 sampling points and a sampling rate of 2 milliseconds (ms). The horizontal axis represents amplitude and the vertical axis represents time, with the unit being seconds (s). Figure 1 (b) shows the Morlet wavelet atoms decomposed using high-precision iterative complex-domain basis tracing. The x-axis represents the number of atoms, and the y-axis represents time in seconds (s). Figure 1 (a) Consistent. Figure 1 (c) is the calculated probability density function with respect to the phase, where the horizontal axis represents the phase in degrees (°) and the vertical axis represents the probability. Figure 1 (d) compares the seismic wavelet extracted using this method (black dashed line) with the theoretical seismic wavelet (gray solid line). The vertical axis represents amplitude, and the horizontal axis represents time, with the unit being seconds (s).

[0053] Figure 2 This is a seismic record synthesized using the Marmous model in an embodiment of the present invention. The horizontal axis represents the trace number, with a total of 681 traces, and the vertical axis represents time in seconds (s). The sampling rate is 2 milliseconds (ms), with 793 sampling points per trace. The two black boxes in the figure indicate the selected locations for extracting seismic wavelets. Region A contains 201 traces, with 301 sampling points per trace and a time range of 0.4 seconds to 1 second. Region B contains 21 traces, with 171 sampling points per trace and a time range of 0.46 seconds to 0.8 seconds.

[0054] Figure 3 The results are seismic wavelet results of mixed phase extracted from regions A and B in the seismic record synthesized by the Marmous model in this embodiment of the invention. Figure 3 (a) shows the seismic data for region A, with the horizontal axis representing the trace number and the vertical axis representing time, in seconds (s). Figure 3 (b) is a partial display of the first 20 Morlet wavelet atoms from the high-precision iterative complex domain basis tracing decomposition results in region A. The horizontal axis represents the number of atoms, and the vertical axis represents the number of atoms. Figure 3 (a) Consistent. Figure 3 (c) is the calculated probability density function with respect to the phase, where the horizontal axis represents the phase in degrees (°) and the vertical axis represents the probability. Figure 3 (d) compares the seismic wavelet extracted using this method (black dashed line) with the theoretical seismic wavelet (gray solid line). The vertical axis represents amplitude, and the horizontal axis represents time, with the unit being seconds (s). Figure 3 (e) shows the seismic data for region B, with the horizontal axis representing the trace number and the vertical axis representing time, in seconds (s). Figure 3 (f) is a partial display of the first 20 Morlet wavelet atoms from the high-precision iterative complex domain basis tracing decomposition results in region B. The horizontal axis represents the number of atoms, and the vertical axis represents the number of atoms. Figure 3 (e) is consistent. Figure 3 (g) is the calculated probability density function with respect to the phase, where the horizontal axis represents the phase in degrees (°) and the vertical axis represents the probability. Figure 3 (h) is a comparison between the seismic wavelet extracted using this method (black dashed line) and the theoretical seismic wavelet (gray solid line). The vertical axis represents amplitude, and the horizontal axis represents time, with the unit being seconds (s). Detailed Implementation

[0055] (1) Import the seismic data S=[S(1),…,S(i),…,S(M)] for seismic wavelet extraction. This data has M channels.

[0056] (2) A high-precision iterative complex domain basis tracing method is used to perform Morlet wavelet decomposition on the selected seismic data channel by channel, as follows:

[0057] ① Calculate the analytical signal AS of the seismic data S;

[0058] ② Calculate the residual SR of the m-th seismic record at the k-th iteration using the following formula. k (m):

[0059]

[0060] Where AI is the analytic Morlet wavelet function, HT(·) denotes the Hilbert transform, and a,t,f,u,σ, These are the amplitude, time, average frequency, time shift, scale factor, and phase of the analytical Morlet wavelet, respectively, where i is the imaginary unit;

[0061] ③ Calculate the residual SR using the following formula k (m) The maximum position of the amplitude envelope, which serves as the time shift u of the analytic Morlet wavelet in the k-th decomposition. k (m):

[0062]

[0063] Where Env(·) represents calculating SR k The amplitude envelope of (m);

[0064] ④ Using u k (m), the analytical Morlet wavelet parameters [σ] of the k-th decomposition of the m-th seismic record are calculated using the following formula. k (m),f k (m)]:

[0065]

[0066] Where <·> is the dot product operator, and |·| is the absolute value symbol. The symbol for the 2-norm;

[0067] ⑤ Calculate the analytical Morlet wavelet amplitude A of the m-th seismic record in the k-th iteration using an optimization algorithm. k (m):

[0068]

[0069] Where ε is the damping factor and I is the identity matrix of dimension N;

[0070] ⑥Use A k (m), calculate the phase of the analytical Morlet wavelet using the following formula. and amplitude a k (m):

[0071]

[0072] Among them, Im[·] is the imaginary part operation of the orientation quantity, Re[·] is the real part operation of the orientation quantity, and abs[·] is the absolute value operation of the orientation quantity;

[0073] ⑦ Repeat steps ②-⑥ above until the required number of iterations is reached;

[0074] (3) Based on the decomposition results in the previous step, calculate the WVD time-frequency distribution WVD(m,t,f) of the m-th channel;

[0075] (7) Calculate the average frequency of the seismic wavelet based on the WVD time-frequency distribution results of the multichannel seismic data from the previous step. The method is as follows:

[0076] The WVD time-frequency distribution is normalized by channel, as shown in the following formula:

[0077]

[0078] Where max[·] represents taking the maximum value;

[0079] use The superimposed amplitude spectrum Amp(f) is calculated using the following formula:

[0080]

[0081] Based on Amp(f), the frequency corresponding to the maximum amplitude of the superimposed amplitude spectrum is obtained as the average frequency of the seismic wavelet to be extracted, as follows:

[0082]

[0083] (8) Based on the high-precision iterative complex domain basis tracing decomposition results of the seismic records, the phase parameters of the analytical Morlet wavelet extracted channel by channel are analyzed. Perform phase folding:

[0084]

[0085] Where φ1 is the folded phase;

[0086] (9) The phase information of the seismic wavelet is obtained using the following method:

[0087] Obtain the absolute value of the Morlet wavelet phase φ2:

[0088] φ2(m,k)=abs[φ1(m,k)];

[0089] The probability density function S of φ2 is calculated using the following formula:

[0090]

[0091] Where f(·) represents the weight of φ2.

[0092] Calculate the average phase corresponding to the maximum value in the probability density function. The formula is as follows:

[0093]

[0094] Where j is the index value of the probability density function S, and mean[·] represents the average value;

[0095] (7) Determine the order of the generalized seismic wavelet based on the phase formula of the generalized seismic wavelet.

[0096]

[0097] (8) According to and Calculating the reference frequency of the generalized seismic wavelet The parameters and formula are as follows:

[0098]

[0099] Where fr is the reference frequency of the generalized seismic wavelet, and Γ(·) is the gamma function;

[0100] (9) According to and The mixed-phase seismic wavelet w(t) is obtained by reconstruction using the following formula:

[0101]

[0102] Figure 1 This is a single-channel seismic record test used in embodiments of the present invention. Figure 1 (a) is a single-channel seismic record synthesized using a 45° phase seismic wavelet, with 201 sampling points and a sampling rate of 2 milliseconds (ms). This invention uses a high-precision iterative complex domain basis tracing method for... Figure 1 (a) Decompose to obtain Figure 1 (b) The average frequency and probability density function of the three Morlet wavelet atoms in this invention are calculated using the method of this invention, and the average frequency is obtained as follows: The average phase is obtained as Finally, the parameter order is calculated to be... and reference frequency Figure 1 (d) shows a comparison between the seismic wavelet extracted using this method (black dashed line) and the theoretical seismic wavelet (gray solid line). The two completely overlap, indicating that the method of this invention is effective and reliable.

[0103] Figure 2 This is a seismic record synthesized using the Marmous model in an embodiment of the present invention. It contains 681 traces, each with 793 sampling points. Two black boxes (A and B) in the figure were selected for seismic wavelet extraction. Region A contains 201 traces, each with 301 sampling points, and represents a complex tectonic region. Region B contains 21 traces, each with 171 sampling points, and represents a relatively gentle tectonic region.

[0104] Figure 3 The results are seismic wavelet results of mixed phase extracted from regions A and B in the seismic record synthesized by the Marmous model in this embodiment of the invention. Figure 3 (a) shows that the seismic signal in region A is complex, therefore requiring more traces for seismic wavelet extraction. This necessitates more iterations for each trace; with 10 iterations per trace, some extraction results are shown below. Figure 3 As shown in (b), its phase information is as follows: Figure 3 As shown in (c), the final calculated parameter order is: and reference frequency Figure 3 (d) shows a comparison between the extracted seismic wavelet (black dashed line) and the theoretical seismic wavelet (gray solid line), which completely overlap. Figure 3 (e) shows that the seismic signal in region B is relatively simple, therefore fewer seismic traces can be used for seismic wavelet extraction. Region B has 20 traces. The iteration count for each trace is set to 4. Partial extraction results are shown below. Figure 3 As shown in (f), its phase information is as follows: Figure 3 As shown in (g), the final calculated parameter order is: and reference frequency Figure 3 (h) shows a comparison between the extracted seismic wavelet (black dashed line) and the theoretical seismic wavelet (gray solid line), which completely overlap. Using seismic data from two tectonic regions of different complexity, A and B, this demonstrates that the method of this invention can extract seismic wavelets using fewer seismic traces. Furthermore, the high-precision iterative complex domain basis tracing for each trace does not require complete decomposition of the seismic signal, allowing for flexible setting of the iteration count to rapidly improve the decomposition of the Morlet wavelet in seismic data. The seismic wavelet extracted based on this method can be further combined with well logging data for use in seismic inversion processes, yielding more reliable seismic inversion results.

[0105] The advantages of this invention are: (1) This invention uses a high-precision iterative complex domain basis tracing method for atomic decomposition and WVD time-frequency analysis technology to deeply process one-dimensional seismic records and extract two-dimensional time-frequency information, which significantly improves the information richness and analytical value of the data; (2) This invention requires only a very small amount of seismic trace data, and can use as little as one trace to efficiently and accurately extract seismic wavelets, greatly reducing data acquisition and time costs; (3) The method of this invention has operational flexibility, allowing the setting of the number of iterations of high-precision iterative complex domain basis tracing, quickly extracting the core information of seismic records, and ensuring efficient, stable and reliable extraction results.

[0106] 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 rapid extraction method for mixed-phase seismic wavelets, comprising the following main steps: (1) input seismic data for seismic wavelet extraction ; (2) The selected seismic data S is decomposed into Morlet wavelet one channel at a time using a high-precision iterative complex domain basis tracing method, as follows: ① Calculate the analytic signal of seismic data S ; The residual of the m-th seismic record at the k-th iteration is calculated using the following formula. : ， in, To analyze the Morlet subwavefunction, Represents the Hilbert transform. These are the amplitude, time, average frequency, time shift, scale factor, and phase of the analytical Morlet wavelet, respectively, where i is the imaginary unit; ③ Calculate the residuals using the following formula. The position of the maximum amplitude envelope is used as the time shift of the analytic Morlet wavelet in the k-th decomposition. : ， in, Expressing the request The amplitude envelope; ④Use The analytical Morlet wavelet parameters of the k-th decomposition of the m-th seismic record are calculated using the following formula. : ， in, This is the dot product operator. It is the absolute value symbol. The symbol for the 2-norm; (5) using an optimization algorithm to calculate the analytical Morlet wavelet analytical amplitude of the mth seismic record at the kth iteration : ， wherein is a damping factor, is an identity matrix of dimension N; ⑥Use The phase of the analytical Morlet wavelet is calculated using the following formula. and amplitude : ， in, Operations on the imaginary part of the orientation quantity. Operations on the real part of the orientation quantity. Absolute value calculation of orientation quantities; ⑦ Repeat the above steps -⑥, until the required number of iterations is reached; (3) According to the decomposition result of the last step, the WVD time-frequency distribution of the mth channel is calculated ; (4) According to the WVD time-frequency distribution result of the multi-channel seismic data in the last step, the average frequency of the seismic wavelet is calculated The method is as follows: The WVD time-frequency distribution is normalized by channel, as shown in the following formula: wherein denotes taking the maximum value; Utilizing , the superimposed amplitude spectrum is calculated as : ， According to , the frequency corresponding to the maximum value of the superimposed amplitude spectrum amplitude is taken as the average frequency of the seismic wavelet to be extracted as follows: ， (5) Based on the high-precision iterative complex domain basis tracing decomposition results of the seismic records, the phase parameters of the analytical Morlet wavelet extracted channel by channel are... Perform phase folding: ， wherein is the folded phase; (6) The phase information of the seismic wavelet is obtained using the following method: Obtain the absolute value of the Morlet wavelet phase : ； calculate probability density function The formula is as follows: ， wherein represents the weight of Computing the average phase corresponding to the maximum of the probability density function The formula is as follows: ， where j is a probability density function index value, denotes taking the average. (7) According to the phase formula of the generalized seismic wavelet below, the order of the generalized seismic wavelet is obtained : ， (8) According to and , the reference frequency of the generalized seismic wavelet is calculated parameters, as follows: ， wherein, is a reference frequency of the generalized seismic wavelet, is the gamma function; (9) According to and The mixed-phase seismic wavelet is obtained by reconstructing using the following formula. : 。

Images