Gas-containing reservoir characterization time-frequency analysis method based on sparse generalized W transformation

Abstract

The invention discloses a gas-bearing reservoir characterization time-frequency analysis method based on sparse generalized W transformation, and the method comprises the steps: introducing an L1 norm into a mathematical relationship between generalized W transformation and a seismic signal for sparse constraint for the first time, and carrying out the solving through employing a Bregman iterative algorithm, thereby obtaining an analysis result with higher time-frequency resolution. According to the method, the advantages of highlighting seismic signal low-frequency information and avoiding the main frequency splitting problem of generalized W transformation are absorbed, and more sparse time-frequency representation is provided for non-stationary seismic signals. When the method is applied to time-frequency analysis of actual seismic data, a high-precision seismic spectrum decomposition result can be provided for the gas-bearing reservoir, so that the low-frequency abnormal region of the gas-bearing reservoir can be delineated more accurately.

Description

Technical Field

[0001] The invention belongs to the field of signal processing and proposes a sparse generalized W transform time-frequency analysis method for characterizing gas-bearing reservoirs. Background Art

[0002] In seismic exploration, the signals used are usually non-stationary signals. Traditional Fourier transform cannot effectively characterize the local frequency characteristics of such signals, while time-frequency analysis can map the one-dimensional time signal to the two-dimensional time-frequency domain, and reveal the geological structure and reservoir information hidden in the seismic data through spectrum decomposition, thus becoming an important tool for processing non-stationary seismic signals.

[0003] Traditional time-frequency analysis methods include short-time Fourier transform, wavelet transform and S transform. S transform has been widely used in seismic signal processing and analysis because its Gaussian window can change adaptively with frequency. However, the change trend of the window size of S transform with frequency is preset and fixed, which cannot provide a fully flexible and reliable time-frequency representation for seismic signals. In view of this, various generalized S transforms have emerged to adjust this change trend, but unfortunately, they still do not solve the problem of the main frequency shifting to high frequency in S transform.

[0004] In recent years, the W transform has effectively highlighted the low-frequency information of seismic signals by constructing a Gaussian window with a main frequency weight so that the centroid of the spectrum energy coincides with the main frequency of the waveform, and successfully solved the main frequency offset problem in the S transform and its generalized transform. However, the Gaussian window of the W transform contains the absolute value of the main frequency, which leads to the splitting of the time-frequency amplitude peak at the main frequency position. To solve this problem, the generalized W transform parameterized the window function and achieved remarkable results, but its time-frequency resolution is still fundamentally restricted by the Heisenberg uncertainty principle, which makes it difficult to simultaneously obtain high time resolution and high frequency resolution in the time-frequency domain. In particular, for complex seismic signals, the frequency components contained in them span multiple scales. The generalized W transform may be blurred or distorted when capturing these components, which limits the precise analysis of local details.

[0005] With the continuous advancement of seismic exploration technology, the requirements for the accuracy of seismic data interpretation are increasing, and high resolution has become an indispensable key factor. At present, the existing time-frequency analysis technology can hardly meet these more stringent requirements in terms of resolution. Therefore, it has become an urgent and important task to develop more effective time-frequency analysis tools to accurately extract key information from seismic data. Summary of the invention

[0006] In view of the shortcomings of the prior art, the present invention proposes a sparse generalized W transform gas-bearing reservoir characterization time-frequency analysis method that can significantly improve energy concentration and improve the accuracy of oil and gas reservoir characterization. Starting from the theoretical basis of sparse representation, on the basis of generalized W transform, the time-frequency transformation process of the signal is regarded as a linear inversion problem under sparse constraints, and the Bregman iterative algorithm is used to solve it, and finally a time-frequency spectrum with higher time-frequency resolution is obtained. Compared with the existing methods, the sparse generalized W transform can highlight the low-frequency information of the signal, better characterize the time-frequency characteristics of the signal, and significantly improve the time-frequency resolution and energy concentration of the time-frequency spectrum.

[0007] In order to achieve the above object, the present invention adopts the following technical solution to achieve it: a sparse generalized W transform gas reservoir characterization time-frequency analysis method, comprising the following steps:

[0008] S1. Input the original two-dimensional seismic profile x(s,t) to be analyzed, where s is the number of channels and t is the time. Analyze each seismic signal x(t) one by one.

[0009] S2. Perform a generalized W transform on the signal, and obtain the time domain and frequency domain expressions of the generalized W transform respectively:

[0010]

[0011] Where τ and f represent time and frequency respectively, f0(τ) is the time-varying main frequency of the signal and Δf(τ)=f-f0(τ), k represents the scale factor, p is the trend factor, and Z(α) is the Fourier transform of x(t);

[0012] S3. Based on the sparsity theory, a new sparse generalized W transform time-frequency decomposition model is constructed using the L1 norm;

[0013] S4, using the Bregman iterative algorithm to accelerate the solution and obtain the sparse generalized W transform time-frequency result of the signal;

[0014] S5, taking the modulus of the time-frequency transform result to obtain the energy of each time-frequency point, thereby obtaining the time-frequency spectrum of the sparse generalized W transform;

[0015] S6. Use Fourier transform to calculate the average spectrum of the original seismic data. Consider two frequency values ​​as the upper and lower band limits. Extract two high- and low-frequency attribute profiles within the range from the time-frequency spectrum obtained in S5. Characterize the gas-bearing reservoir through the energy attenuation of the seismic signals in these two profiles.

[0016] Preferably, after obtaining the frequency domain expression of the signal through step S2, step S3 constructs the signal sparse generalized W transform spectrum decomposition model process as follows:

[0017] a. Using the frequency domain expression of the generalized W transform, its discrete form is:

[0018]

[0019] Among them, n and l represent discrete time and frequency respectively, and N represents the number of sampling points;

[0020] b. Rewrite formula (2) in matrix form as follows:

[0021] W=AX=AFx (3)

[0022] in, is a column vector of the discrete form of the generalized W transform of the signal GWT[n,l], is the matrix generated by the window function, X∈C N×1 is the column vector generated by X(α), F∈C N×N is the Fourier transform matrix, x is the input signal;

[0023] c. According to formula (3), the new signal sparse generalized W transform time-frequency decomposition model is:

[0024]

[0025] Where G is the window function matrix of the generalized W transform and the Fourier transform matrix F∈C N×N Pseudo-inverse of the product, λ ≥ 0 is the balancing parameter.

[0026] Preferably, the solution process using the Bregman iterative algorithm in step S4 is:

[0027]

[0028] Here, i is the number of iterations.

[0029] The present invention has the following beneficial effects:

[0030] The present invention proposes a sparse generalized W transform gas reservoir characterization time-frequency analysis method, which introduces the L1 norm into the mathematical relationship between the generalized W transform and the seismic signal for sparse constraint for the first time, and uses the Bregman iterative algorithm to solve it, thereby obtaining an analysis result with higher time-frequency resolution. This method absorbs the advantages of the generalized W transform in highlighting the low-frequency information of the seismic signal and avoiding the main frequency splitting problem, and provides a sparser time-frequency representation for non-stationary seismic signals.

[0031] The present invention proposes a sparse generalized W transform method for characterizing gas-bearing reservoirs. The synthetic signal shows that, compared with the generalized S transform, W transform and its generalized transform, the sparse generalized W transform obtains a time-frequency analysis result with higher time-frequency resolution that highlights low-frequency information and avoids main frequency splitting through sparse constraints. After performing sparse generalized W transform on two-dimensional seismic profiles, the analysis of two high- and low-frequency attribute profiles extracted shows that the sparse generalized W transform can effectively characterize gas-bearing reservoirs and can more clearly depict the location and boundaries of the reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 (a) is the synthetic signal;

[0033] Figure 1 (b) is the time-varying main frequency of the synthetic signal;

[0034] Figure 1 (c) is the generalized S transform result of the synthetic signal;

[0035] Figure 1 (d) is the W transformation result of the synthetic signal;

[0036] Figure 1 (e) is the generalized W transform result of the synthetic signal;

[0037] Figure 1 (f) is the sparse generalized W transform result of the synthetic signal;

[0038] Figure 2 A total of 1079 seismic data profiles for an oil and gas field;

[0039] Figure 3 (a) and (b) are the low-frequency profile (28 Hz) and high-frequency profile (42 Hz) extracted by generalized S-transform, respectively;

[0040] Figure 3 (c) and (d) are the low-frequency profile (28 Hz) and high-frequency profile (42 Hz) extracted by W transform, respectively;

[0041] Figure 3 (e) and (f) are the low-frequency profile (28 Hz) and high-frequency profile (42 Hz) extracted by generalized W transform, respectively;

[0042] Figure 3 (g) and (h) are the low-frequency profile (28 Hz) and high-frequency profile (42 Hz) extracted by the sparse generalized W transform, respectively. DETAILED DESCRIPTION

[0043] A sparse generalized W transform time-frequency analysis method for characterizing gas-bearing reservoirs will be described in detail below in conjunction with the accompanying drawings of the present invention.

[0044] Example 1

[0045] The sparse generalized W transform proposed in the present invention comprises the following steps:

[0046] S1. Input the original two-dimensional seismic profile x(s,t) to be analyzed, where s is the number of channels and t is the time. Analyze each seismic signal x(t) one by one.

[0047] S2. Select appropriate parameters k and p to perform generalized W transform on the signal, and the time domain and frequency domain expressions of generalized W transform are:

[0048]

[0049] Where τ and f represent time and frequency respectively, f0(τ) is the time-varying main frequency of the signal and Δf(τ)=f-f0(τ), k represents the scale factor, p is the trend factor, and X(α) is the Fourier transform of x(t);

[0050] S3. According to the frequency domain expression of the generalized W transform, its discrete form is:

[0051]

[0052] Where n and l represent discrete time and frequency respectively, and N represents the total length of the signal;

[0053] S4. Rewrite the discrete form of the generalized W transform into a matrix expression:

[0054] W=AX=AFx (3)

[0055] in, is a column vector of the discrete form of the generalized W transform of the signal GWT[n,l], is the matrix generated by the window function, X∈C N×1 is the column vector generated by X(α), F∈C N×N is the Fourier transform matrix, x is the input signal;

[0056] S5, introduce the L1 norm, select the appropriate parameter λ to construct a new sparse generalized W transform time-frequency decomposition model:

[0057]

[0058] Where G is the window function matrix of the generalized W transform and the Fourier transform matrix F∈C N×N Pseudo-inverse matrix of the product, λ≥0 is the balance parameter;

[0059] S6. Use the Bregman iterative algorithm to accelerate the solution. The solution process is as follows:

[0060]

[0061] Here, i is the number of iterations.

[0062] S7, obtaining the sparse generalized W transform time-frequency result of the signal through step S6, and taking the modulus of the result to obtain the sparse generalized W transform time-frequency spectrum;

[0063] The performance of this method on synthetic signals is as follows Figure 1 shown. Figure 1 (a) is the synthetic signal, which is formed by the convolution of the reflection coefficient and the Ricker wavelet with the main frequencies of 50Hz, 40Hz, 30Hz and 20Hz respectively. Figure 1 (b) shows the time-varying dominant frequency obtained by Hilbert transform, Figure 1 (c)-(f) correspond to the time-frequency analysis results of generalized S transform, W transform, generalized W transform, and sparse generalized W transform, respectively. The dotted line represents the time-varying main frequency. It can be seen that the generalized S transform provides a rough time-frequency representation for the signal. The time resolution in the low-frequency region is low and the low-frequency characteristics of the signal cannot be highlighted. The W transform introduces the time-varying main frequency, which greatly improves the time resolution in the low-frequency region, but there is a splitting phenomenon at the main frequency. The generalized W transform improves the main frequency splitting problem of the W transform by parameterizing the time-varying main frequency function. However, due to the limitations of the uncertainty principle, the above methods all show obvious energy diffusion, and there is serious interference between seismic signals. The sparse generalized W transform has a higher time resolution in the low-frequency region, significantly reduces the interference of seismic signals, and highlights the low-frequency information.

[0064] Example 2

[0065] The sparse generalized W transform proposed in the present invention is applied to the characterization of gas-bearing reservoirs in a gas field in western Sichuan. The specific steps are as follows:

[0066] S1-S7, same as in Example 1;

[0067] S8. Use Fourier transform to calculate the average spectrum of the original seismic data. Consider two frequency values ​​as the upper and lower band limits. Extract two high (42 Hz) and low (28 Hz) frequency attribute sections within the range from the time-frequency spectrum obtained in S7. The energy attenuation of the seismic signals in these two sections is used to characterize the gas-bearing reservoir.

[0068] Figure 2 This is a 2D seismic profile of a gas field in western Sichuan. The profile contains 1079 traces, with 251 sampling points per trace and a sampling interval of 2ms. The two vertical lines A and B are two gas wells, and the box is the gas reservoir development area. The low-frequency and high-frequency profiles corresponding to the generalized S transform, W transform, generalized W transform, and sparse generalized W transform are shown in Figure 2. Figure 3(a) and (b), (c) and (d), (e) and (f), (g) and (h). Note that the generalized S transform cannot well represent the low-frequency information of the location of well A at low frequencies. The energy of the study area shown by the other three methods is strong in the low-frequency band and weak in the high-frequency band, which can well describe the abnormal energy attenuation caused by the gas-bearing reservoir at the location of wells A and B. However, the sparse generalized W transform can more clearly depict the boundary of the reservoir, and further achieve high-precision delineation of the abnormal attenuation area of ​​the gas-bearing reservoir. Compared with the generalized S transform, W transform and generalized W transform, the time-frequency energy concentration of the sparse generalized W transform is significantly improved, and the reservoir characteristics described are more accurate.

[0069] The above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit the same. Although the present invention has been described in detail with reference to the above embodiments, a person skilled in the art can still modify or make equivalent substitutions to the specific implementation schemes of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention are within the scope of protection of the claims of the present invention.

Claims

1. A sparse generalized W transform gas reservoir characterization time-frequency analysis method, characterized in that: The following steps are involved: S1. Input the original two-dimensional seismic profile x(s,t) to be analyzed, where s is the number of channels and t is the time. Analyze each seismic signal x(t) one by one. S2. Perform a generalized W transform on the signal, and obtain the time domain and frequency domain expressions of the generalized W transform respectively: Where τ and f represent time and frequency respectively, f0(v) is the time-varying main frequency of the signal and Δf(τ)=f-f0(τ), k represents the scale factor, p is the trend factor, and X(α) is the Fourier transform of x(t); S3. Based on the sparsity theory, a new sparse generalized W transform time-frequency decomposition model is constructed using the L1 norm; S4, using the Bregman iterative algorithm to accelerate the solution and obtain the sparse generalized W transform time-frequency result of the signal; S5, taking the modulus of the time-frequency transform result to obtain the energy of each time-frequency point, thereby obtaining the time-frequency spectrum of the sparse generalized W transform; S6. Use Fourier transform to calculate the average spectrum of the original seismic data. Consider two frequency values ​​as the upper and lower band limits. Extract two high- and low-frequency attribute profiles within the range from the time-frequency spectrum obtained in S5. Characterize the gas-bearing reservoir through the energy attenuation of the seismic signals in these two profiles.

2. The method according to claim 1, characterized in that Step S3: After obtaining the frequency domain expression of the signal through step S2, the process of constructing a new sparse generalized W transform time-frequency decomposition model is as follows: a. Using the frequency domain expression of the generalized W transform, its discrete form is: Among them, n and l represent discrete time and frequency respectively, and N represents the number of sampling points; b. Rewrite formula (2) in matrix form as follows: W=AX=AFx (3) in, is a column vector of the discrete form of the generalized W transform of the signal GWT[n,l], is the matrix generated by the window function, X∈C N×1 is the column vector generated by X(α), F∈C N×N is the Fourier transform matrix, x is the input signal; c. According to formula (3), the new signal sparse time-frequency spectrum decomposition model is: Where G is the window function matrix of the generalized W transform and the Fourier transform matrix F∈C N×N Pseudo-inverse of the product, λ ≥ 0 is the balancing parameter.

3. The method according to claim 2, characterized in that: In step S4, the solution process using the Bregman iterative algorithm is: Here, i is the number of iterations.

Images