Fast extraction method for mixed phase seismic wavelets

Abstract

The invention discloses a rapid extraction method for mixed phase seismic wavelets. According to the method, main feature information of seismic data is extracted by adopting a high-precision iterative complex field basis tracking method, an instantaneous spectrum of the seismic data is efficiently calculated by utilizing analyzed WVD time-frequency distribution, and frequency and phase information of the seismic data is quickly obtained and used for extracting final seismic wavelets. Under the condition of lack of logging data, the seismic wavelets containing the mixed phase information can be accurately and efficiently extracted from limited seismic data, uncertainty caused by insufficient data is effectively avoided, and a key basis is provided for subsequent seismic data processing.

Description

Technical Field

[0001] The invention belongs to the field of oil and gas seismic exploration and relates to a method for quickly extracting mixed-phase seismic wavelets. Background Art

[0002] In oil and gas seismic exploration, extracting seismic wavelets is the key basis for seismic data processing. It is very important in processing and is used in core links such as deconvolution, well seismic calibration, and forward and inversion of seismic records. Seismic wavelets are the direct response of the source signal propagating in the underground medium, accurately showing the waveform shaping effect of seismic waves from the source to the receiver. The main purpose of extracting seismic wavelets is to extract wavelet characteristic information from complex seismic records, obtain key parameters such as the speed 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 process of seismic wavelet extraction, well logging data plays a vital role and is a key factor in ensuring the accuracy of seismic wavelet extraction. For areas lacking well logging data, the task of accurately extracting seismic wavelets is particularly arduous. We need to use the statistical information of the seismic traces to extract wavelets to ensure the reliability of the results. However, using the statistical seismic wavelet extraction method based on the statistical information of the seismic traces themselves, the acquired wavelets are usually in the form of zero phase, minimum phase or maximum phase, which makes it difficult to directly obtain the true phase information. In addition, seismic signals are often accompanied by significant noise and interference information, which further aggravates the deviation between the extracted wavelet and the actual wavelet, thereby increasing the multi-solution problem of the statistical seismic wavelet extraction method. This situation poses a challenge in practical applications, making it difficult to determine the optimal seismic wavelet extraction scheme. Summary of the invention

[0004] In view of the shortcomings of the prior art, the purpose of the present invention is to provide a method for rapid extraction of mixed-phase seismic wavelets, which is characterized by using a high-precision iterative complex domain basis tracking method and WVD time-frequency analysis technology to directly and accurately extract the average frequency and phase information of the seismic wavelet, which not only quickly responds to the complex situation of the actual seismic signal, but also has high precision and stability. The method of the present invention includes the following main steps:

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

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

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

[0008] ② According to the following formula, calculate the residual SR of the mth seismic record at the kth iteration: k (m):

[0009]

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

[0011] ③ Calculate the residual SR as follows k (m) The maximum position of the amplitude envelope, which is the time shift u of the analytical Morlet wavelet of the kth decomposition k (m):

[0012]

[0013] Among them, Env(·) represents the SR k (m) the amplitude envelope;

[0014] ④Use u k (m), the analytical Morlet wavelet parameter [σ k (m),f k (m)]:

[0015]

[0016] Among them, <·> is the inner product operator symbol, |·| is the absolute value symbol, is the 2-norm symbol;

[0017] ⑤ Use the optimization algorithm to calculate the analytical Morlet wavelet amplitude A of the mth seismic record at the kth iteration k (m):

[0018]

[0019] Among them, ε is the damping factor, I is the unit matrix with dimension N;

[0020] ⑥Use A k (m), the phase of the analytical Morlet wavelet is calculated as follows and amplitude a k (m):

[0021]

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

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

[0024] (3) According to the decomposition result in the previous step, calculate the WVD time-frequency distribution WVD(m,t,f) of the mth channel;

[0025] (4) Calculate the average frequency of the seismic wavelet based on the WVD time-frequency distribution results of the multi-channel seismic data in the previous step Here’s how:

[0026] The WVD time-frequency distribution is normalized by channel, and the formula is as follows:

[0027]

[0028] Among them, max[·] means taking the maximum value;

[0029] use The superimposed amplitude spectrum Amp(f) is calculated as follows:

[0030]

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

[0032]

[0033] (5) Based on the complex domain basis pursuit decomposition results of the high-precision iteration of the seismic records, the phase parameters of the analytical Morlet wavelet extracted from each channel are To perform phase folding:

[0034]

[0035] Among them, φ 1 is the phase after folding;

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

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

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

[0039] Calculate φ 2 The probability density function S of is as follows:

[0040]

[0041] Where f(·) represents φ 2 Weight

[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[·] means taking the average value;

[0045] (7) According to the following generalized seismic wavelet phase formula, the order of the generalized seismic wavelet is obtained:

[0046]

[0047] (8) According to and Calculate the reference frequency of generalized seismic wavelets Parameters, the formula is as follows:

[0048]

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

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

[0051] BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 It is the seismic wavelet result with a phase of 45° extracted from the single-channel theoretical seismic record used in the embodiment 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 is amplitude and the vertical axis is time, with the unit being seconds (s). Figure 1 (b) is the Morlet wavelet atom decomposed by high-precision iterative complex domain basis tracking. The horizontal axis is the number of atoms and the vertical axis is time in seconds (s). Figure 1 (a) Consistency. Figure 1 (c) is the calculated probability density function about the phase, the horizontal axis is the phase, the unit is angle (°), and the vertical axis is the probability. Figure 1(d) is a comparison between the seismic wavelet extracted using this method (black dotted line) and the theoretical seismic wavelet (gray solid line), the ordinate is amplitude, the abscissa is time, and the unit is seconds (s).

[0053] Figure 2 This is a seismic record synthesized by the Marmous model in an embodiment of the present invention. The horizontal axis is the channel number, with a total of 681 channels, and the vertical axis is time, the unit is second (s), the sampling rate is 2 milliseconds (ms), and there are 793 sampling points per channel. The two black boxes in the figure are the positions selected for extracting seismic wavelets. There are 201 channels in area A, 301 sampling points per channel, and the time range is 0.4 seconds to 1 second. There are 21 channels in area B, 171 sampling points per channel, and the time range is 0.46 seconds to 0.8 seconds.

[0054] Figure 3 It is the seismic wavelet result of mixed phase extracted from the two regions A and B in the seismic record synthesized by the Marmous model in the embodiment of the present invention. Figure 3 (a) is the seismic data of area A, the horizontal axis is the track number, the vertical axis is the time, and the unit is second (s). Figure 3 (b) is a partial display of the first 20 Morlet wavelet atoms of the complex domain basis pursuit decomposition result of high-precision iteration in region A. The horizontal axis is the number of atoms and the vertical axis is Figure 3 (a) Consistency. Figure 3 (c) is the calculated probability density function about the phase, the horizontal axis is the phase, the unit is angle (°), and the vertical axis is the probability. Figure 3 (d) is a comparison between the seismic wavelet extracted using this method (black dotted line) and the theoretical seismic wavelet (gray solid line), the ordinate is amplitude, the abscissa is time, and the unit is seconds (s). Figure 3 (e) is the seismic data of area B, with the horizontal axis being the track number and the vertical axis being the time in seconds (s). Figure 3 (f) is a partial display of the first 20 Morlet wavelet atoms of the complex domain basis pursuit decomposition result of high-precision iteration in region B. The horizontal axis is the number of atoms and the vertical axis is Figure 3 (e) consistency. Figure 3 (g) is the calculated probability density function about the phase, the horizontal axis is the phase, the unit is angle (°), and the vertical axis is the probability. Figure 3 (h) is a comparison between the seismic wavelet extracted using this method (black dotted line) and the theoretical seismic wavelet (gray solid line), the ordinate is amplitude, the abscissa is time, and the unit is seconds (s). DETAILED DESCRIPTION

[0055] (1) Importing seismic data S = [S(1), ..., S(i), ..., S(M)] for seismic wavelet extraction, the data has a total of M channels;

[0056] (2) Using a high-precision iterative complex domain basis tracking method, the selected seismic data are decomposed into Morlet wavelets one by one, as follows:

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

[0058] ② According to the following formula, calculate the residual SR of the mth seismic record at the kth iteration: k (m):

[0059]

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

[0061] ③ Calculate the residual SR as follows k (m) The maximum position of the amplitude envelope, which is the time shift u of the analytical Morlet wavelet of the kth decomposition k (m):

[0062]

[0063] Among them, Env(·) represents the SR k (m) the amplitude envelope;

[0064] ④Use u k (m), the analytical Morlet wavelet parameter [σ k (m),f k (m)]:

[0065]

[0066] Among them, <·> is the inner product operator symbol, |·| is the absolute value symbol, is the 2-norm symbol;

[0067] ⑤ Use the optimization algorithm to calculate the analytical Morlet wavelet amplitude A of the mth seismic record at the kth iteration k (m):

[0068]

[0069] Among them, ε is the damping factor, I is the unit matrix with dimension N;

[0070] ⑥Use A k(m), the phase of the analytical Morlet wavelet is calculated as follows and amplitude a k (m):

[0071]

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

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

[0074] (3) According to the decomposition result in the previous step, calculate the WVD time-frequency distribution WVD(m,t,f) of the mth channel;

[0075] (7) Calculate the average frequency of the seismic wavelet based on the WVD time-frequency distribution results of the multi-channel seismic data in the previous step Here’s how:

[0076] The WVD time-frequency distribution is normalized by channel, and the formula is as follows:

[0077]

[0078] Among them, max[·] means taking the maximum value;

[0079] use The superimposed amplitude spectrum Amp(f) is calculated as follows:

[0080]

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

[0082]

[0083] (8) Based on the complex domain basis pursuit decomposition results of high-precision iteration of seismic records, the phase parameters of the analytical Morlet wavelet extracted from each trace are To perform phase folding:

[0084]

[0085] Among them, φ 1 is the phase after folding;

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

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

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

[0089] Calculate φ 2 The probability density function S of is as follows:

[0090]

[0091] Where f(·) represents φ 2 Weight

[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[·] means taking the average value;

[0095] (7) According to the following generalized seismic wavelet phase formula, the order of the generalized seismic wavelet is obtained:

[0096]

[0097] (8) According to and Calculate the reference frequency of generalized seismic wavelets Parameters, the formula is as follows:

[0098]

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

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

[0101]

[0102] Figure 1 is a single-channel seismic recording test used in the embodiment 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). The present invention uses a high-precision iterative complex domain basis tracking method to Figure 1 (a) is decomposed to obtain Figure 1 The average frequency and probability density function of the three Morlet wavelet atoms in (b) are calculated using the method of the present invention, and the average frequency is The average phase is Finally, the parameter order is calculated as and the reference frequency is Figure 1 (d) is a comparison between the seismic wavelet extracted using the present method (black dashed line) and the theoretical seismic wavelet (grey solid line). The two completely coincide with each other, indicating that the present method is effective and reliable.

[0103] Figure 2 This is a seismic record synthesized by the Marmous model in the embodiment of the present invention. There are 681 channels in total, with 793 sampling points in each channel. Two black boxes (A and B) in the figure are selected for extracting seismic wavelets. Among them, there are 201 channels in area A, with 301 sampling points in each channel, which is a complex structural area. There are 21 channels in area B, with 171 sampling points in each channel, which is a relatively flat structural area.

[0104] Figure 3 It is the seismic wavelet result of mixed phase extracted from the two regions A and B in the seismic record synthesized by the Marmous model in the embodiment of the present invention. Figure 3 The seismic signal in area A shown in (a) is complex, so more traces are needed to participate in the seismic wavelet extraction. More iterations are required for each seismic data. The number of iterations per trace is set to 10. Some of the extraction results are shown in Figure 2. Figure 3 As shown in (b), its phase information is Figure 3 As shown in (c), the parameter order is finally calculated as and the reference frequency is Figure 3 (d) shows the comparison between the extracted seismic wavelet (black dashed line) and the theoretical seismic wavelet (gray solid line), and the two completely overlap. Figure 3 The seismic signal in area B shown in (e) is relatively simple, so fewer seismic traces can be used for seismic wavelet extraction, 20 traces in area B. The number of iterations per trace is set to 4, and some of the extraction results are shown in Figure 1. Figure 3 As shown in (f), its phase information is Figure 3 As shown in (g), the parameter order is finally calculated as and the reference frequency is Figure 3(h) shows the comparison between the extracted seismic wavelet (black dotted line) and the theoretical seismic wavelet (gray solid line), and the two completely overlap. The seismic data of two structural areas of different complexity, A and B, are used here to show that the method of the present invention can use fewer seismic traces to extract seismic wavelets, and the high-precision iterative complex domain basis tracking of each trace does not require a complete decomposition of the seismic signal, and the number of iterations can be flexibly set to quickly improve the decomposition of the Morlet wavelet of the seismic data. The seismic wavelet extracted based on the method of the present invention can be further combined with logging data and used in the seismic inversion process to obtain more reliable seismic inversion results.

[0105] The advantages of the present invention are: (1) the present invention adopts a high-precision iterative complex domain basis tracking method for atomic decomposition, as well as WVD time-frequency analysis technology, to deeply process one-dimensional seismic records and extract two-dimensional time-frequency spectrum information, thereby significantly improving the information richness and analysis value of the data; (2) the present invention only requires a very small amount of seismic trace data, as little as one trace, to efficiently and accurately extract seismic wavelets, greatly reducing data acquisition and time costs; (3) the method of the present invention has operational flexibility, allowing the number of high-precision iterative complex domain basis tracking iterations to be set, so as to quickly extract the core information of the seismic records and ensure efficient, stable and reliable extraction results.

[0106] 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 rapidly extracting mixed-phase seismic wavelets, comprising the following main steps: (1) Input seismic data S for seismic wavelet extraction; (2) A high-precision iterative complex domain basis tracking method is used to perform Morlet wavelet decomposition on the selected seismic data S one by one, as follows: ① Calculate the analytical signal AS of the seismic data S; ② According to the following formula, calculate the residual SR of the mth seismic record at the kth iteration: k (m): in, AI is the analytical Morlet wavelet function, HT(·) represents the Hilbert transform, a,t,f,u,σ, are the amplitude, time, average frequency, time shift, scale factor and phase of the analytic Morlet wavelet, respectively, and i is an imaginary unit; ③ Calculate the residual SR as follows k (m) The maximum position of the amplitude envelope, which is the time shift u of the analytical Morlet wavelet of the kth decomposition k (m): Among them, Env(·) represents the SR k (m) Amplitude envelope; ④Use u k (m), the analytical Morlet wavelet parameter [σ k (m),f k (m)]: Among them, <·> is the inner product operator symbol, |·| is the absolute value symbol, is the 2-norm symbol; ⑤ Use the optimization algorithm to calculate the analytical Morlet wavelet amplitude A of the mth seismic record at the kth iteration k (m): Among them, ε is the damping factor, I is the unit matrix with dimension N; ⑥Use A k (m), the phase of the analytical Morlet wavelet is calculated as follows and amplitude a k (m): Among them, Im[·] is the imaginary part of the orientation quantity, Re[·] is the real part of the orientation quantity, and abs[·] is the absolute value of the orientation quantity. ⑦ Repeat the above steps ②-⑥ until the number of iterations is reached; (3) According to the decomposition result in the previous step, calculate the WVD time-frequency distribution WVD(m,t,f) of the mth channel; (4) Calculate the average frequency of the seismic wavelet based on the WVD time-frequency distribution results of the multi-channel seismic data in the previous step Here’s how: The WVD time-frequency distribution is normalized by channel, and the formula is as follows: Among them, max[·] means taking the maximum value; use The superimposed amplitude spectrum Amp(f) is calculated as follows: According to Amp(f), the frequency corresponding to the maximum value of the stacked amplitude spectrum is obtained as the average frequency of the seismic wavelet to be extracted as follows: (5) Based on the complex domain basis pursuit decomposition results of the high-precision iteration of the seismic records, the phase parameters of the analytical Morlet wavelet extracted from each channel are To perform phase folding: Among them, φ1 is the phase after folding; (6) The phase information of the seismic wavelet is obtained by the following method: Get the absolute value of the Morlet wavelet phase φ2: φ2(m,k)=abs[φ1(m,k)]; Calculate the probability density function S of φ2, the formula is as follows: Where f(·) represents the weight of φ2 Calculate the average phase corresponding to the maximum value in the probability density function The formula is as follows: Where j is the index value of the probability density function S, and mean[·] means taking the average value; (7) According to the following generalized seismic wavelet phase formula, the order of the generalized seismic wavelet is obtained: (8) According to and Calculate the reference frequency of generalized seismic wavelets Parameters, the formula is as follows: Where fr is the reference frequency of the generalized seismic wavelet, Γ(·) is the gamma function; (9) According to and The mixed phase seismic wavelet w(t) is reconstructed using the following formula:

Images