A reservoir gas-bearing property identification method based on frequency-domain basis stretching Chirplet transform

Abstract

The present invention discloses a method for identifying gas-bearing properties of reservoirs based on frequency-domain basis-stretched Chirplet transform, comprising the following steps: S1, sequentially extracting single-channel seismic data; S2, performing Fourier transform on each signal extracted in S1 to obtain a corresponding frequency-domain seismic signal; S3, constructing a kernel phase function for matching complex signals; S4, windowing the frequency-domain seismic signal in S2 based on the kernel phase function in S3 to obtain a frequency-domain basis-stretched Chirplet transform result corresponding to each data channel; S5, calculating the dominant frequency range of the wellbore signal of the seismic profile to be measured using Fourier transform, and selecting high and low frequencies based on this; extracting low-frequency and high-frequency profiles from the result in S4, and identifying gas-bearing properties of the reservoir by comparing the attenuation characteristics of the high-frequency and low-frequency profiles. The present invention uses a frequency-domain basis-stretched Chirplet transform method to process seismic data, which can significantly improve time-frequency focusing and provide a reliable basis for identifying gas-bearing properties of reservoirs.

Description

Technical Field

[0001] The present invention belongs to the field of oil and gas seismic signal processing and analysis, and in particular relates to a method for identifying reservoir gas content based on frequency domain basis stretching Chirplet transform. Background Art

[0002] With the rapid development of society and the economy, energy issues have become a primary constraint on global development. As oil and gas exploration and development continue to deepen, readily accessible reservoirs have been depleted, prompting researchers to focus on harder-to-find complex reservoirs. However, the complex and variable physical properties of rocks, such as significant heterogeneity and lateral discontinuities, pose significant challenges to the exploration and identification of complex reservoirs. Time-frequency analysis is a powerful tool for processing seismic data, providing a two-dimensional approach to interpreting the time-frequency characteristics of seismic signals. Related research has demonstrated that oil and gas responses exhibit a "low-frequency resonance, high-frequency attenuation" pattern (Yin et al. 2015; Castagna et al. 2020). Time-frequency analysis can effectively identify the response characteristics of geological bodies in different frequency bands (Liu et al. 2017) and is now widely used in the identification of complex reservoirs.

[0003] Traditional time-frequency analysis methods mainly fall into two categories: classical and parametric. Classical methods, such as the short-time Fourier transform (STFT), continuous wavelet transform (CWT), and S-transform (ST), use windowed Fourier transforms to analyze the local time-frequency characteristics of nonstationary signals and have been widely applied in many fields. To improve the adaptability of strong frequency signals, several parametric time-frequency analysis methods have been proposed. These methods use parameterized phase functions to estimate parameters based on signal characteristics. When these parameters match the signal's high-order phases, the signal's time-frequency focusing performance is optimal. These methods include the Chirplet transform (CT) and polynomial Chirplet transform (PCT) for single-component signals, and the general linear Chirplet transform (GLCT) and synchronous compensated Chirplet transform (SCCT) for multi-component signals. Although these methods offer high time-frequency resolution, most cannot process complex multi-component signals without prior knowledge, and they also struggle to precisely characterize complex seismic signals. Summary of the Invention

[0004] To address the above-mentioned shortcomings of the existing technology, the present invention provides a method for identifying reservoir gas content based on frequency-domain basis scaling Chirplet transform. This method significantly improves the time-frequency energy focusing of complex multi-component signals, enhancing its adaptability to seismic signals and the accuracy of reservoir gas content identification.

[0005] To achieve the above object, the specific implementation steps adopted by the present invention include:

[0006] S1. Extract single-channel earthquake signals from seismic data in sequence, and extract each earthquake signal as x. i (t), where i = (1, 2, ..., N), N is the total number of seismic profile data, and t is time;

[0007] S2, perform Fourier transform on each signal extracted in S1 to obtain each frequency domain seismic signal X i (f), is calculated as follows:

[0008]

[0009] Where t and f represent time and frequency respectively;

[0010] S3. Construct a kernel phase function T Γ (t,f,f0) is used to match complex signals, where t and f0 are the time and frequency centers, and f is the frequency;

[0011] S4, according to the kernel phase function in S3, the frequency domain seismic signal X in S2 i (f) Perform windowing calculation to obtain the frequency domain basis stretching Chirplet transform result corresponding to each data

[0012] S5. Use Fourier transform to calculate the main frequency range of the wellbore signal of the seismic profile to be measured, and find the frequency value F at the peak amplitude in the Fourier spectrum. max As the low frequency, select 1.5*F max As the high frequency; extract the low-frequency and high-frequency profiles respectively from the results obtained in S4, and identify the gas content of the reservoir by comparing the attenuation characteristics of the high-frequency and low-frequency profiles.

[0013] Preferably, in step S3, the kernel phase function T Γ The expression is as follows:

[0014] T Γ (t,f,f0)=t×(f+β1(f-f0) 2 +β2(f-f0) 3 +…+β n (f-f0) n+1 )

[0015] Where t and f0 are the time and frequency centers, f is the frequency, Γ=(β1,β2,β3,…,β n ) are the polynomial coefficients of the kernel phase function.

[0016] Preferably, in step S3, since the group delay GD(f) of an arbitrary signal is defined as the inverse of the frequency derivative of the phase function, we have:

[0017]

[0018] Among them, T Γ The specific expression of ′(t,f,f0) is:

[0019] T Γ ′(t,f,f0)=t×(1+2β1(f-f0)+3β2(f-f0) 2 +…+(n+1)β n (f-f0) n )

[0020] Its second-order derivative T Γ ″ is T Γ ′With respect to the derivative of frequency, there is the following relationship:

[0021]

[0022] Where α is the rotation angle of the TF basis at the frequency center, and α changes with time t. Γ The specific expression of " is:

[0023] T Γ ″(t,f,f0)=t×(2β1+6β2(f-f0)+…+n(m+1)β n (f-f0) n-1 )When f=f0, T Γ "for:

[0024] T Γ =t×2β1=-tan(α)

[0025] Where α is the rotation angle of the TF basis at the frequency center f0, then:

[0026] α=-arctan(g(f0)×2β1)

[0027] where g(f0) is the function value of the GD trajectory at f0. Let the inclination angle of the GD trajectory at f0 be θ. When β1 is accurately estimated using a certain method, the TF basis matches the GD trajectory, and α = θ.

[0028] As an example, in step S3, the kernel phase function T Γ , when n is 1, T Γ ′ is expressed as:

[0029] T Γ ′(t,f,f0)=t×(1+2β1(f-f0))

[0030] At this point, the linear group delay trajectory can be characterized. In order to better match the complex signal containing the nonlinear group delay trajectory and reduce the amount of calculation as much as possible, n is generally taken as 2. At this time, T Γ ′ is expressed as:

[0031] T Γ ′(t,f,f0)=t×(1+2β1(f-f0)+3β2(f-f0) 2 )

[0032] As an example, in step S4, the frequency domain seismic signal X i (f) Perform window calculation to obtain the corresponding result for each data Finally, the polynomial coefficient Γ of the kernel phase function is determined according to the kurtosis theory, and finally the result with the highest time-frequency focusing is obtained. The specific expression of the frequency domain basis stretching Chirplet transform is:

[0033]

[0034] where X i (f) is the frequency domain seismic signal obtained by S3, T Γ (t,f,f0) is the kernel phase function, H σ (f) is the Gaussian window function h σ (t) is represented in the frequency domain, H σ (f) and h σ The specific expression of (t) is:

[0035]

[0036]

[0037] Preferably, in step S4, the method for determining the polynomial coefficient Γ of the kernel phase function according to the kurtosis theory is as follows:

[0038]

[0039] When G max When the maximum value is reached, the estimated parameters Γ=(β1,β2,β3,…,β n ) takes the optimal solution, where N is the total number of seismic profile data, m4 is the fourth-order sample moment about the mean, m2 is the second-order sample moment about the mean, is the mean, is the processing result of the i-th seismic profile data.

[0040] The idea of the present invention is:

[0041] First, extract the single-channel earthquake signal from the seismic data in sequence, and extract each seismic signal as x i (t), where i = (1, 2, ..., N), N is the total number of seismic profile data, and t is time;

[0042] Second, perform Fourier transform on each signal extracted from S1 to obtain each frequency domain seismic signal X i (f), is calculated as follows:

[0043]

[0044] Where t and f represent time and frequency respectively;

[0045] Third, construct a kernel phase function T Γ (t,f,f0) is used to match complex signals, where t and f0 are the time and frequency centers, and f is the frequency;

[0046] Fourth, according to the kernel phase function in the third step, the frequency domain seismic signal X in the second step is i (f) Perform windowing calculation to obtain the frequency domain basis stretching Chirplet transform result corresponding to each data

[0047] Fifth, use Fourier transform to calculate the main frequency range of the wellbore signal of the seismic profile to be measured, and find the frequency value F at the peak amplitude in the Fourier spectrum. max As the low frequency, select 1.5*F max As high frequency; extract low-frequency and high-frequency profiles respectively from the results obtained in the fourth step, and identify the gas content of the reservoir by comparing the attenuation characteristics of the high-frequency and low-frequency profiles.

[0048] The working principle of the present invention is to extract single-channel earthquake signals from seismic data in sequence, and extract each seismic signal as x i (t); Perform Fourier transform on each signal extracted in the first step to obtain each frequency domain seismic signal X i (f); Construct a kernel phase function T Γ (t,f,f0) is used to match complex signals; according to the kernel phase function in the third step, the frequency domain seismic signal X in the second step i (f) Perform windowing calculation to obtain the frequency domain basis stretching Chirplet transform result corresponding to each data Use Fourier transform to calculate the main frequency range of the wellbore signal of the seismic profile to be measured, and find the frequency value F at the peak amplitude in the Fourier spectrum. max As the low frequency, select 1.5*F maxAs the high frequency, the results obtained in step 4 are used to extract low-frequency and high-frequency profiles, respectively. Gas-bearing properties of the reservoir are identified by comparing the attenuation characteristics of the high- and low-frequency profiles. This invention approaches seismic data from a frequency-domain perspective. Based on the abnormal energy and frequency attenuation experienced when seismic waves pass through gas-bearing reservoirs, a frequency-domain basis-stretched Chirplet transform method is used to process seismic data. This method significantly improves time-frequency focusing and provides a reliable basis for identifying gas-bearing properties in reservoirs.

[0049] This paper addresses the multi-component, non-stationary, and complex nonlinear signals found in actual seismic data from a frequency-domain perspective, proposing a method for identifying reservoir gas content based on the frequency-domain basis-stretching Chirplet transform. This method effectively processes signals whose time-frequency ridges are biased toward the frequency direction, effectively matching the signal's group delay (GD) trajectory. Because seismic waves experience abnormal energy and frequency attenuation when passing through gas-bearing reservoirs, this method first processes the seismic data using the frequency-domain basis-stretching Chirplet transform to produce a high-frequency energy focus. This method then uses these characteristics to better identify reservoir gas content, providing a reliable basis for practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 Flowchart of the present invention;

[0051] Figure 2 is a multi-component frequency domain synthetic signal and its time-frequency spectrum comparison diagram, where Figure 2 (a) is a multi-component frequency domain synthetic signal, Figure 2 (b)- Figure 2 (f) Several time-frequency analysis methods are used to Figure 2 The results obtained after processing the signal in (a) are: Figure 2 (b) Short frequency Fourier transform (SFFT), Figure 2 (c) Continuous Wavelet Transform (CWT), Figure 2 (d) Synchro-Squeeze Transform (SST), Figure 2 (e) Frequency domain polynomial Chirplet transform (FPCT), Figure 2 (f) Frequency domain basis stretching Chirplet transform (FSBCT);

[0052] Figure 3 This is the original seismic profile of the actual data;

[0053] Figure 4 For Figure 3Comparison of the low-frequency (28 Hz) and high-frequency (40 Hz) profiles generated after the data are processed. The left column shows the low-frequency case, and the right column shows the high-frequency case. The methods used are: (a)-(b) short-frequency Fourier transform (SFFT), (c)-(d) continuous wavelet transform (CWT), (e)-(f) synchronized squeezing transform (SST), (g)-(h) frequency-domain polynomial chirplet transform (FPCT), and (i)-(j) frequency-domain basis stretching chirplet transform (FSBCT). DETAILED DESCRIPTION

[0054] The present invention will be further described below with reference to the accompanying drawings.

[0055] Example 1: See Figure 1 A method for identifying gas content in reservoirs based on frequency domain basis stretching Chirplet transform includes the following steps:

[0056] S1. Extract single-channel earthquake signals from seismic data in sequence, and extract each earthquake signal as x. i (t), where i = (1, 2, ..., N), N is the total number of seismic profile data, and t is time;

[0057] S2, perform Fourier transform on each signal extracted in S1 to obtain each frequency domain seismic signal X i (f), is calculated as follows:

[0058]

[0059] Where t and f represent time and frequency respectively;

[0060] S3. Construct a kernel phase function T Γ (t,f,f0) is used to match complex signals, where t and f0 are the time and frequency centers, and f is the frequency;

[0061] S4, according to the kernel phase function in S3, the frequency domain seismic signal X in S2 i (f) Perform windowing calculation to obtain the frequency domain basis stretching Chirplet transform result corresponding to each data

[0062] S5. Use Fourier transform to calculate the main frequency range of the wellbore signal of the seismic profile to be measured, and find the frequency value F at the peak amplitude in the Fourier spectrum. max As the low frequency, select 1.5*F max As the high frequency; extract the low-frequency and high-frequency profiles respectively from the results obtained in S4, and identify the gas content of the reservoir by comparing the attenuation characteristics of the high-frequency and low-frequency profiles.

[0063] Example 2: See Figure 2 This example shows a synthetic multi-component frequency domain signal and uses this model to test the performance of the proposed method for time-frequency focusing. Figure 3 Several existing methods were selected for comparison with the method of the present invention. Figure 2 The methods selected are SFFT, CWT, SST, FPCT and FSBCT in order. Figure 2 (b) with Figure 2 As can be seen from the results in (c), the processing results of SFFT and CWT only show slightly better focusing at 2Hz-4Hz, and the effects in other parts are poor. This is because these two methods are constrained by the Heisenberg uncertainty principle and cannot achieve both time resolution and frequency resolution. Figure 2 (d) shows the processing result of SST. Although the energy in the entire frequency band is relatively focused, the existence of its tuning effect will have a significant impact on the time-frequency analysis results. Figure 2 (e) with Figure 2 (f) shows the processing results of FPCT and FSBCT respectively. It can be seen that both methods show good results, but the energy of the GD trajectory at both ends of FPCT at 0-2Hz and 5-8Hz diverges slightly. In comparison, FSBCT shows better time-frequency focusing. This model proves that the method of the present invention has better performance.

[0064] Example 3, see Figure 3-Figure 4 , taking the tight sandstone gas reservoir of a gas field as the seismic data to be analyzed, Figure 3 For the original seismic profile, each seismic data contains 551 sampling points, with a sampling interval of 2ms, and the entire data contains 286 seismic records. It is known that well A is a well-developed gas well and well B is a dry well. For each seismic data, Figure 1 The steps shown in the figure are processed. According to the spectrum of the wellbore reflection layer signal, the frequency range corresponding to the abnormal frequency area is found (28 Hz-42 Hz in this case). The frequency value of the high-frequency and low-frequency profile comparison is determined. Then, the high and low frequencies are selected for single-frequency profile extraction. Then, the gas content of the reservoir is directly identified based on the characteristics of "low-frequency energy enhancement and high-frequency energy attenuation".

[0065] Figure 4 The high and low frequency profiles after processing by several time-frequency analysis methods are shown. It can be seen that the common feature of the five processing methods is the low frequency profile (see Figure 4 The energy in the rectangular box intersecting with well A on (a)(c)(e)(g)(i)) is higher, while the high-frequency profile (see Figure 4There is obvious energy attenuation at the corresponding positions on (b)(d)(f)(h)(j)). Based on the characteristics of "low-frequency energy enhancement and high-frequency energy attenuation", the above methods can all verify that well A is a gas-bearing well and well B is not a gas-bearing well, indicating that all five methods can effectively identify the gas-bearing properties of tight sandstones. Comparing SFFT, CWT, and FPCT, it is not difficult to find that the method of the present invention has a more refined characterization, its time-frequency resolution is higher, and the energy attenuation of FSBCT is more obvious. Comparing the results with the classic post-processing method SST, the advantages of FSBCT can also be found. Therefore, through the above analysis, it can be seen that the frequency domain basis stretching Chirplet transform can more accurately extract the abnormal response characteristics of seismic signals and can more accurately describe tight sandstone gas layers.

[0066] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for identifying reservoir gas content based on frequency domain basis scaling Chirplet transform, characterized in that: The following steps are involved: S1. Extract single-channel earthquake signals from seismic data in sequence, and extract each earthquake signal as x. i (t), where i = (1, 2, ..., N), N is the total number of seismic profile data, and t is time; S2, perform Fourier transform on each signal extracted in S1 to obtain each frequency domain seismic signal X i (f), is calculated as follows: Where t and f represent time and frequency respectively; S3. Construct a kernel phase function T Γ (t,f,f0) is used to match complex signals, where t and f0 are the time and frequency centers, and f is the frequency; S4, according to the kernel phase function in S3, the frequency domain seismic signal Xi(f) in S2 is windowed and the frequency domain basis stretching Chirplet transform result corresponding to each data is obtained. The calculation formula is as follows: where X i (f) is the frequency domain seismic signal obtained by S2, T Γ (t,f,f0) is the kernel phase function obtained by S3, H σ (f) is the Gaussian window function h σ (t) is represented in the frequency domain, σ is the standard deviation, H σ (f) and h σ The specific expression of (t) is: S5. Use Fourier transform to calculate the main frequency range of the wellbore signal of the seismic profile to be measured, and find the frequency value F at the peak amplitude in the Fourier spectrum. max As the low frequency, select 1.5*F max As the high frequency; extract the low-frequency and high-frequency profiles respectively from the results obtained in S4, and identify the gas content of the reservoir by comparing the attenuation characteristics of the high-frequency and low-frequency profiles.

2. The method for identifying reservoir gas content based on frequency domain basis scaling Chirplet transform according to claim 1, characterized in that: In step S3, the kernel phase function T Γ The expression is as follows: T Γ (t,f,f0)=t×(f+β1(f-f0) 2 +β2(f-f0) 3 +…+β n (f-f0) n+1 ) Where t and f0 are the time and frequency centers, f is the frequency, Γ=(β1,β2,β3,…,β n ) are the polynomial coefficients of the kernel phase function.

3. The method according to claim 2, characterized in that Since the group delay GD(f) of any signal is defined as the inverse of the frequency derivative of the phase function, we have: Among them, T Γ The specific expression of ′(t,f,f0) is: T Γ ′(t,f,f0)=t×(1+2β1(f-f0)+3β2(f-f0) 2 +…+(n+1)β n (f-f0) n ) Its second-order derivative T Γ ″ is T Γ ′With respect to the derivative of frequency, there is the following relationship: Where α is the rotation angle of the TF basis at the frequency center, and α changes with time t, T Γ The specific expression of " is: T Γ ″(t,f,f0)=t×(2β1+6β2(f-f0)+…+n(n+1)β n (f-f0) n-1 ) Among them, when f=f0, T Γ "for: T Γ ″=t×2β1=-tan(α) Where α is the rotation angle of the TF basis at the frequency center f0, then: α=-arctan(g(f0)×2β1) Where g(f0) is the function value of the GD trajectory at f0; let the inclination angle of the GD trajectory at f0 be θ. When β1 is accurately estimated using a certain method, the TF basis matches the GD trajectory, and at this time α = θ.

4. The method according to claim 3, characterized in that When n is 1, T Γ ′ is expressed as: T Γ ′(t,f,f0)=t×(1+2β1(f-f0)) At this point, the linear group delay trajectory can be characterized. In order to better match the complex signal containing the nonlinear group delay trajectory and reduce the amount of calculation as much as possible, n is generally taken as 2. At this time, T Γ ′ is expressed as: T Γ ′(t,f,f0)=t×(1+2β1(f-f0)+3β2(f-f0) 2 ) 5. The method according to claim 1, wherein The method for determining the polynomial coefficient Γ of the kernel phase function is as follows: When G max When the maximum value is reached, the estimated parameters Γ=(β1,β2,β3,…,β n ) takes the optimal solution, where N is the total number of seismic profile data, m4 is the fourth-order sample moment about the mean, m2 is the second-order sample moment about the mean, is the mean, is the processing result of the i-th seismic profile data.

Images