Method and system for identifying gas reservoir based on array acoustic logging time-frequency analysis

By combining frequency domain bandpass filtering, time domain median filtering, and two-dimensional variational mode decomposition with a smooth pseudo-Wigner-Ville distribution algorithm, the problems of weak signal and noise interference in gas layer identification in array acoustic logging were solved, and high-precision identification of gas layers in complex reservoirs was achieved.

CN121069495APending Publication Date: 2025-12-05ZHANJIANG BRANCH OF CHINA NATIONAL OFFSHORE OIL CORP
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202511263019.X
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2025-09-05
Publication Date
2025-12-05

AI Technical Summary

Technical Problem

Existing array acoustic logging technology has difficulty identifying gas layers in complex reservoirs. The signals are weak and noise interference is significant, failing to effectively separate gas-sensitive parts, resulting in inaccurate identification results.

Method used

A combination of frequency domain bandpass filtering and time domain median filtering is used for noise reduction. Two-dimensional variational mode decomposition and smooth pseudo-Wigner-Ville distribution algorithm are combined to extract high-frequency mode components that are sensitive to the gas content of the formation. Gas layers are identified by time-spectrum energy attenuation characteristics.

Benefits of technology

It effectively filters out high-frequency reflection signals in a strong noise background, improving the extraction accuracy of gas layer signals and the accuracy of formation gas content assessment.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN121069495A_ABST
    Figure CN121069495A_ABST
Patent Text Reader

Abstract

The invention discloses a method and a system for identifying a gas reservoir based on array acoustic logging time-frequency analysis. The method comprises the following steps of: acquiring array acoustic logging data; multi-stage noise combined suppression is carried out; variational mode decomposition and high-frequency signal component extraction are carried out; identifying a gas layer based on time-frequency characteristics and the like. The system comprises an acquisition module, a denoising module, a decomposition module, an extraction module and a time frequency module. According to the method, array acoustic logging information is utilized, enhanced extraction of weak gas reservoir signals under the strong noise background is achieved through staged processing of frequency domain-time domain multi-stage filtering, two-dimensional variational mode decomposition and time-frequency analysis, the problems that under the complex geological condition, the gas-bearing characteristic response of a reservoir is weak, and the gas reservoir recognition precision of a conventional method is low are effectively solved, and the method is suitable for being used for gas reservoir recognition. The method has good application value in gas bearing evaluation of low-permeability reservoirs, and provides technical support for exploration reserve evaluation and well location implementation.
Need to check novelty before this filing date? Find Prior Art

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of geophysical well logging, and particularly relates to a method and system for identifying gas layers based on array acoustic logging time-frequency analysis. BACKGROUND

[0002] With the global oil and gas exploration and development expanding to deep layers and unconventional reservoirs, complex reservoirs such as low-permeability sandstone and shale gradually become important fields for resource replacement. Such reservoirs generally have strong heterogeneity, complex pore structure, and diverse fluid occurrence states, which lead to challenges in gas layer identification, such as weak signals, significant noise interference, and difficulty in extracting gas-bearing characteristics. Array acoustic logging technology can obtain elastic parameters and fluid response information of the surrounding formation by collecting full-wave train data (P-wave, S-wave, and Stoneley wave), which provides an important basis for reservoir gas-bearing evaluation.

[0003] In recent years, the rapid development of multi-scale signal processing technology provides a new technical path for the analysis of array acoustic data. Frequency domain filtering technology can effectively suppress specific noise interference through band separation, while variational mode decomposition (VMD) as an adaptive signal decomposition method can separate mode components of different frequency bands according to the intrinsic characteristics of the signal. Meanwhile, time-frequency analysis technology (such as Wigner-Ville distribution) can accurately depict the time-frequency energy spectrum of the signal due to its high resolution characteristics. However, under complex reservoir conditions, the response signal of the gas layer is weak and the noise is diverse. Most current array acoustic logging time-frequency analysis methods directly perform time-frequency analysis on the full waveform without considering the separate separation of the part sensitive to gas-bearing. In addition, in the noise processing part, the received formation reflection signals are not separated, which may interfere with the final identification results.

[0004] In order to better realize the identification of gas layers in complex reservoirs, a more accurate, objective, and effective method for comprehensive evaluation of the gas-bearing property of the formation is needed. SUMMARY

[0005] The present application is proposed to solve the problems in the prior art, and the purpose is to provide a method and system for identifying gas layers based on array acoustic logging time-frequency analysis.

[0006] A method for identifying gas layers based on array acoustic logging time-frequency analysis, comprising the following steps:

[0007] S1, acquiring array acoustic full-wave train data of a target well;

[0008] S2, for the array acoustic full-wave train data of the target well, first performing frequency domain band-pass filtering, and then performing time domain median filtering on the data after frequency domain band-pass filtering;

[0009] S3, performing two-dimensional variational modal decomposition on the data processed in step S2 to obtain K modal signal components of different center frequencies through adaptive decomposition;

[0010] S4, extracting a high-frequency modal component sensitive to gas-bearing property of a formation from the K modal signal components based on an energy-frequency joint screening criterion, and reconstructing the high-frequency modal component meeting the criterion;

[0011] S5, calculating a time-frequency spectrum of the high-frequency modal component reconstructed in step S4 by using a smooth pseudo Wigner-Ville distribution algorithm to obtain a time-frequency spectrum of the array sonic logging at different depths, and identifying a gas layer according to a spectral energy attenuation characteristic of the time-frequency spectrum.

[0012] In the above technical solution, the frequency range of the frequency domain band-pass filtering in step S2 is 1 kHz-20 kHz; the window length of the time domain median filtering is 5-15 sampling points; and K is a positive integer not less than 3 in step S3.

[0013] In the above technical solution, the energy-frequency joint screening criterion in step S4 is that the center frequency of the separated modal signal component is greater than 6 kHz and the energy proportion is greater than 2%.

[0014] In the above technical solution, the calculation formula of the center frequency of the modal signal component is:

[0015]

[0016] In the formula, f k is the center frequency of the modal component, unit: Hz, ω k =(ω k,x ,ω k,y ) is the angular frequency corresponding to the energy focus center of the modal component in the two-dimensional frequency domain, unit: rad / s, ω k,x and ω k,y are the center angular frequencies of the modal component in the x direction and the y direction, unit: rad / s, ||ω k ||is the modulus of ω k , unit: rad / s, k=1-K, representing the kth modal signal component;

[0017] The calculation formula of the energy proportion of the modal signal component is:

[0018]

[0019] In the formula, E k % is the energy proportion of the modal signal component, unit: %; E k is the energy of the kth modal component; and E total is the sum of the energies of the K modal components decomposed.

[0020] E k The calculation formula of u is:

[0021]

[0022] In the formula, u k (x,y) is the time domain expression of the kth modal signal component;

[0023] E total The calculation formula of u is:

[0024]

[0025] In the formula, E k is the energy of the kth modal component.

[0026] In the above technical solution, the calculation formula of the high-frequency modal component reconstruction is:

[0027]

[0028] In the formula, u n is the time domain representation of the nth high-frequency modal component, and N is the total number of screened high-frequency modal components.

[0029] In the above technical solution, the calculation formula of the time-frequency spectrum calculated by the smoothing pseudo Wigner-Ville distribution algorithm in step S5 is:

[0030]

[0031] In the formula, z(t) is the analytic signal of the real signal x(t) (obtained by Hilbert transform), t is the time variable, f is the frequency variable, g(u) is the time domain smoothing window, which controls the smoothing degree in the time direction, h(τ) is the frequency domain smoothing window, which controls the smoothing degree in the frequency direction, u is the time integral variable, τ is the time delay variable (the offset of signal autocorrelation), and * is the complex conjugate operator.

[0032] In the above technical solution, the calculation formula of the spectral energy corresponding to the time-frequency spectrum P(t,f) at different depths d of the array acoustic logging in step S5 is:

[0033]

[0034] In the formula, d is the calculation depth, with the unit of m; t is the time variable, with the unit of s; and f is the frequency variable, with the unit of Hz.

[0035] In the above technical solution, the principle of identifying gas layers according to the attenuation characteristics of the time-frequency spectrum energy in step S5 is that the spectral energy attenuation rate of the high-frequency modal component of the gas layer section is greater than 30% compared with the reference layer spectral energy.

[0036] In the technical solution, the reference layer is the target layer section or an oil layer or a water layer with similar physical properties, the spectral energy of the corresponding high-frequency modal component of the reference layer does not attenuate or the attenuation rate is less than 10%; the layer with an attenuation interval in the range of 10% to 30% is considered to be a layer with poor physical properties or poor gas-bearing properties, such as a poor gas layer, which does not reach the standard of a gas layer.

[0037] In the technical solution, the calculation formula of the spectral energy attenuation rate of the high-frequency modal component of the gas layer section compared with the spectral energy of the reference layer is as follows:

[0038]

[0039] In the formula, E_spec% is the spectral energy attenuation rate at the calculation depth, the unit is %, E_spec is the spectral energy at the calculation depth, E_spec is the spectral energy of the selected reference layer. 基准层

[0040] A system for identifying a gas layer based on array acoustic logging time-frequency analysis for implementing the foregoing method, the system comprises:

[0041] An acquisition module is configured to acquire array acoustic full-wave train data of a target well.

[0042] A denoising module is configured to perform frequency domain band-pass filtering and time domain median filtering on the array acoustic full-wave train data of the target well.

[0043] A decomposition module is configured to decompose the denoised array acoustic full-wave train data based on a two-dimensional variational modal decomposition method to obtain K modal signal components with different center frequencies.

[0044] An extraction module is configured to extract a high-frequency modal component sensitive to the gas-bearing property of the formation from the K modal signal components based on an energy-frequency joint screening criterion.

[0045] A time-frequency module is configured to calculate the time-frequency spectrum of the high-frequency modal component based on a smoothed pseudo Wigner-Ville distribution algorithm and obtain the corresponding spectral energy.

[0046] The present application has the following advantages:

[0047] ​The application provides a gas layer identification method and system based on array acoustic logging time-frequency analysis, random noise and reflected noise are filtered through joint frequency domain band pass filtering and time domain median filtering; K mode signal components of different center frequencies are obtained by 2D VMD decomposition of the denoised full wave train data, the high frequency mode component sensitive to gas bearing is selected through energy-frequency joint screening, the time-frequency spectrum of the high frequency mode component is calculated through the method of smooth pseudo Wigner-Ville distribution, and the corresponding spectral energy is obtained. Compared with the traditional array acoustic time-frequency analysis method which directly performs time-frequency analysis on the full waveform and does not consider removing the high frequency reflected signal, the application extracts the high frequency component sensitive to the formation gas bearing in the full wave train, filters out the influence of high frequency reflected noise, realizes the enhancement extraction and identification of weak gas layer signal in strong noise background, and improves the precision of formation gas bearing evaluation. BRIEF DESCRIPTION OF DRAWINGS

[0048] Figure 1 is a method flowchart of the application;

[0049] Figure 2 is a comparison schematic diagram of array acoustic logging full wave train original data before and after band pass filtering to remove random noise and median filtering to remove reflected noise in embodiment 1 of the application;

[0050] Figure 3 is a schematic diagram of K mode signal components decomposed by two-dimensional variational mode decomposition of the denoised array acoustic full wave train data in embodiment 1 of the application;

[0051] Figure 4 is a schematic diagram of the high frequency signal reconstructed based on the energy-frequency joint screening criterion in embodiment 1 of the application;

[0052] Figure 5 is a schematic diagram of the time-frequency spectrum calculated by using the smooth pseudo Wigner-Ville distribution on the reconstructed high frequency signal in embodiment 1 of the application.

[0053] For those skilled in the art, other related drawings can be obtained from the above drawings without creative labor. DETAILED DESCRIPTION

[0054] In order for those skilled in the art to better understand the technical solutions of the application, the technical solutions of the application will be further described below by combining with the drawings of the specification and through specific embodiments.

[0055] Embodiment 1

[0056] As shown in Figure 1 , a gas layer identification method based on array acoustic logging time-frequency analysis includes the following steps:

[0057] S1, array acoustic logging data acquisition:

[0058] Obtaining array acoustic full wave train data of the target well;

[0059] S2, multi-stage noise joint suppression

[0060] The array acoustic full wave train data of the target well is first filtered by frequency domain band-pass filtering to eliminate environmental noise and random interference, and then the data after band-pass filtering is filtered by time domain median filtering to filter out the received formation reflection wave signal, so that the waveform data after filtering out the formation reflection wave signal can be obtained.

[0061] The frequency range of the frequency domain band-pass filtering is 1kHz-20kHz;

[0062] The window length of the time domain median filtering is generally 5-15 sampling points;

[0063] Figure 2 A comparison diagram before and after band-pass filtering and reflection removal of array acoustic logging full wave train original data provided by the embodiment of the present application is provided. Among them, Figure 2 The left graph shown in the figure is the array acoustic logging full wave train original data (selecting a certain measured well 3460-3480m sandstone reservoir depth section as an example, wherein 3465-3472m is a gas layer, and 3472m below is a water layer), Figure 2 The right graph shown in the figure is the waveform data after band-pass filtering to remove random noise and median filtering (window length is 7 sampling points) to remove reflection noise.

[0064] S3, two-dimensional variational mode decomposition

[0065] The data processed in step S2 is subjected to two-dimensional variational mode decomposition, and K mode signal components of different center frequencies are adaptively decomposed; wherein K is a positive integer not less than 3, generally 6-10;

[0066] K is the number of modes decomposed, and too small K value will cause low frequency and high frequency components to alias, reducing the extraction accuracy of high frequency components; too large K value will cause high frequency information to be scattered to multiple modes, reducing the calculation efficiency. Therefore, K value needs to be selected according to the complexity of the signal, and in the embodiment of the present application, K is an integer not less than 3. Preferably, the value range of K is [6, 10].

[0067] S4, high frequency signal component extraction

[0068] Based on the energy-frequency joint screening criterion, the high frequency mode component sensitive to the gas bearing property of the formation is extracted from the K mode signal components, and the high frequency mode component meeting the criterion is reconstructed;

[0069] Because the gas-bearing formation absorbs high-frequency signals more strongly, the high-frequency energy in the full wave train of the acoustic wave is attenuated significantly, while the low-frequency energy is relatively stable, therefore, based on the energy-frequency joint screening principle, the high-frequency modal component sensitive to the gas-bearing formation is extracted;

[0070] The energy-frequency joint screening principle is that the center frequency of the separated modal signal component is greater than 6 kHz and the energy proportion is greater than 2%;

[0071] The center frequency of the modal signal component represents the energy concentration point of the modal component in the frequency domain, and represents the dominant frequency, the calculation formula of the center frequency of the modal signal component is:

[0072]

[0073] In the formula, f k is the center frequency of the modal component, the unit is Hz, ω k =(ω k,x ,ω k,y ) is the angle frequency corresponding to the energy focusing center of the modal component in the two-dimensional frequency domain, the unit is rad / s, ω k,x and ω k,y are the center angle frequencies of the modal component in the x direction and the y direction, the unit is rad / s, ||ω k || is the modulus of ω k , the unit is rad / s, k=1~K, representing the kth modal signal component.

[0074] The energy proportion of the modal signal component represents the proportion of the energy of the modal component in the total energy of the signal, the calculation formula of the energy proportion of the modal signal component is:

[0075]

[0076] In the formula, E k % is the energy proportion of the modal signal component, the unit is %; E k is the energy of the kth modal component; E total is the sum of the energies of the K modal components decomposed.

[0077] The calculation formula of E k is:

[0078]

[0079] In the formula, u k (x,y) is the time domain expression of the kth modal signal component;

[0080] The calculation formula of E total is:

[0081]

[0082] E = ∑k=1Kak k is the energy of the kth modal component.

[0083] The calculation formula of the high-frequency modal component reconstruction is:

[0084]

[0085] u = x * z n is the time-domain representation of the nth high-frequency modal component, and N is the total number of screened high-frequency modal components.

[0086] Figure 3 is a schematic diagram of K modal signal components decomposed from the denoised array acoustic full-wave data by two-dimensional variational modal decomposition according to an embodiment of the present application (in this embodiment, K is 7), as shown in Figure 3 The first image in the upper left corner is the denoised array acoustic full-wave data, and the subsequent ones are the kth modal signal components (k = 1, 2,..., 7) decomposed in turn; Figure 4 is a schematic diagram of the reconstructed high-frequency signal based on the energy-frequency joint screening criterion, wherein the screened high-frequency modal components k are 4, 5, 6, and 7.

[0087] S5, gas layer identification based on time-frequency characteristics

[0088] The time-frequency spectrum is calculated by using the smoothed pseudo Wigner-Ville distribution algorithm on the reconstructed high-frequency modal components of step S4, to obtain the time-frequency spectrum P(t, f) of the array acoustic logging at different depths d (indicating the time-frequency energy distribution at depth d), and the gas layer is identified according to the spectral energy attenuation characteristics of the time-frequency spectrum.

[0089] The calculation formula of the smoothed pseudo Wigner-Ville distribution algorithm for calculating the time-frequency spectrum is:

[0090]

[0091] wherein z(t) is the analytic signal of the real signal x(t) (obtained by Hilbert transform), t is the time variable, f is the frequency variable, g(u) is the time-domain smoothing window, which controls the smoothing degree in the time direction, h(τ) is the frequency-domain smoothing window, which controls the smoothing degree in the frequency direction, u is the time integral variable, τ is the time delay variable (the offset of the signal autocorrelation), and * is the complex conjugate operator.

[0092] After the smoothed pseudo Wigner-Ville distribution is calculated, the time-frequency spectrum P(t, f) of the array acoustic logging at different depths d (indicating the time-frequency energy distribution at depth d) is obtained, and the spectral energy corresponding to the time-frequency spectrum of the array acoustic logging at different depths is:

[0093]

[0094] wherein d is the calculated depth, unit is m; t is the time variable, unit is s; f is the frequency variable, unit is Hz.

[0095] The principle of identifying the gas layer according to the attenuation characteristics of the time-frequency spectrum energy is:

[0096] The spectrum energy of the high-frequency mode component of the gas layer section is compared with the spectrum energy of the reference layer, and the spectrum energy attenuation rate E_spec% is greater than 30%;

[0097] The reference layer is the target layer section or an oil layer or a water layer similar in physical property to the target layer section, and the spectrum energy of the corresponding high-frequency mode component of the reference layer is basically not attenuated or the attenuation rate is less than 10%;

[0098] The calculation formula of the spectrum energy attenuation rate E_spec% of the high-frequency mode component of the gas layer section compared with the spectrum energy of the reference layer is:

[0099]

[0100] wherein E_spec% is the spectrum energy attenuation rate at the calculated depth, unit is %, E_spec is the spectrum energy at the calculated depth, E_spec 基准层 is the spectrum energy of the selected reference layer.

[0101] Figure 5 According to the embodiment of the present application, a time-frequency spectrum diagram of the extracted high-frequency signal is obtained by using a smooth pseudo Wigner-Ville distribution, as shown in Figure 5 The left graph is the time-frequency spectrum of the high-frequency signal at the depth of 3468m (gas layer), and the right graph is the time-frequency spectrum of the high-frequency signal at the depth of 3476m (water layer). It can be seen that, taking the spectrum energy of the water layer as the reference, the time-frequency spectrum energy of the high-frequency signal of the gas layer is obviously weaker, and the calculated spectrum energy attenuation rate of the high-frequency signal of the gas layer is 73%.

[0102] The array acoustic logging received array acoustic full wave train data includes direct wave, reflected wave and random noise, the formation gas bearing identification method provided in the application is based on the high frequency signal in the direct wave for evaluation, therefore the high frequency signal in the direct wave is the target signal, the reflected wave and the random noise are interference signals and need to be filtered out; the time-frequency analysis gas bearing evaluation method provided in the application can extract the high frequency signal sensitive to the formation gas bearing in the acoustic full wave train for subsequent research; the application realizes the extraction and identification of the weak gas layer signal in the strong noise background by the method of the frequency domain-time domain multi-stage filtering and the 2D VMD decomposition, and the corresponding spectral energy is obtained by calculating the time-frequency spectrum of the high frequency signal through the smoothing pseudo Wigner-Ville distribution method. Compared with the traditional array acoustic time-frequency analysis method, the application extracts the high frequency component sensitive to the formation gas bearing in the full wave train, and filters out the influence of the reflected noise, realizes the extraction and identification of the weak gas layer signal in the strong noise background, and the precision of the reservoir gas bearing evaluation is improved.

[0103] Embodiment 2

[0104] A system for identifying gas layers based on array acoustic logging time-frequency analysis of the method in embodiment 1, the system comprises:

[0105] An acquisition module is configured to acquire array acoustic full wave train data of a target well.

[0106] A denoising module is configured to perform frequency domain band-pass filtering and time domain median filtering on the array acoustic full wave train data of the target well.

[0107] A decomposition module is configured to decompose the denoised array acoustic full wave train data based on a two-dimensional variational mode decomposition method to obtain K modal signal components with different center frequencies.

[0108] An extraction module is configured to extract high frequency modal components sensitive to formation gas bearing from the K modal signal components based on an energy-frequency joint screening criterion.

[0109] A time-frequency module is configured to calculate the time-frequency spectrum of the high frequency modal components based on a smoothing pseudo Wigner-Ville distribution algorithm and obtain the corresponding spectral energy.

[0110] The applicant declares that the above description is only a specific embodiment of the application, but the protection scope of the application is not limited thereto, and those skilled in the art should understand that any changes or replacements within the technical scope disclosed in the application can be easily thought of by those skilled in the art, and all fall within the protection scope and disclosure scope of the application.

Claims

1. A method for identifying gas zones based on array sonic time-frequency analysis, characterized in that: The method comprises the following steps: S1, acquiring array acoustic full wave train data of a target well; S2, performing frequency domain band pass filtering on the array acoustic full wave train data of the target well, and then performing time domain median filtering on the data filtered in the frequency domain; S3, performing two-dimensional variational mode decomposition on the data processed in step S2, and adaptively decomposing to obtain K mode signal components with different center frequencies; S4, extracting high-frequency mode components sensitive to formation gas content from the K mode signal components based on an energy-frequency joint screening criterion, and reconstructing the high-frequency mode components meeting the criterion; S5, calculating a time-frequency spectrum of the high-frequency mode components reconstructed in step S4 by using a smooth pseudo Wigner-Ville distribution algorithm, obtaining a time-frequency spectrum of array acoustic logging at different depths, and identifying a gas layer according to the spectral energy attenuation characteristics of the time-frequency spectrum.

2. The method of claim 1, wherein: In the step S2, the frequency range of the frequency domain band pass filtering is 1 kHz-20 kHz; the window length of the time domain median filtering is 5-15 sampling points; and K is a positive integer not less than 3 in the step S3.

3. The method of claim 1, wherein: In the step S4, the energy-frequency joint screening criterion is that the center frequency of the separated mode signal component is greater than 6 kHz and the energy proportion is greater than 2%.

4. The method of claim 3, wherein: The calculation formula of the center frequency of the mode signal component is: wherein: f k is the center frequency of the modal component, in Hz, ω k = (ω k,x , ω k,y ) is the energy focus center of the modal component in the two-dimensional frequency domain, in rad / s, ω k,x and ω k,y are the center angular frequencies of the modal component in the x direction and the y direction, respectively, in rad / s, ||ω k || is the modulus of ω k , in rad / s, k = 1 ~ K, representing the kth modal signal component; The calculation formula of the energy proportion of the mode signal component is: In the formula, E k is the energy percentage of the kth modal signal component, in %; E k is the energy of the kth modal component; E total is the sum of the energies of the K decomposed modal components; E The calculation formula is: k The calculation formula is: wherein: u k (x,y) is a time domain representation of the kth modal signal component; E The calculation formula is: total The calculation formula is: wherein: E k is the energy of the kth modal signal component.

5. The method of claim 3, wherein: The calculation formula of the high-frequency mode component reconstruction is: wherein: u n is the time domain representation of the nth high frequency modal component, and N is the total number of the screened high frequency modal components.

6. The method of claim 1, wherein: The calculation formula of the time-frequency spectrum calculated by the smooth pseudo Wigner-Ville distribution algorithm in the step S5 is: In the formula, z(t) is an analytic signal of a real signal x(t), t is a time variable, f is a frequency variable, g(u) is a time domain smoothing window, h(τ) is a frequency domain smoothing window, u is a time integral variable, τ is a time delay variable, and * is a complex conjugate operator.

7. The method of claim 1, wherein: The calculation formula of the spectral energy corresponding to the time-frequency spectrum P(t,f) of array acoustic logging at different depths d in the step S5 is: In the formula, d is the calculation depth, with the unit of m; t is a time variable, with the unit of s; and f is a frequency variable, with the unit of Hz. In the step S5, the principle of identifying a gas layer according to the attenuation characteristics of the time-frequency spectrum energy is that the spectral energy attenuation rate of the high-frequency mode component of a gas layer section is greater than 30% compared with the spectral energy of a reference layer; the reference layer is a target layer section or an oil layer or a water layer similar in physical property thereto, and the spectral energy of the high-frequency mode component corresponding to the reference layer does not attenuate or has an attenuation rate less than 10%.

8. The method of claim 1, wherein: The calculation formula of the spectral energy attenuation rate of the high-frequency mode component of the gas layer section compared with the spectral energy of the reference layer is:

9. The method of claim 8, wherein: The system comprises: wherein: E_spec% is the spectral energy attenuation rate at the calculated depth, in %, E_spec is the spectral energy at the calculated depth, E_spec 基准层 is the spectral energy of the selected reference layer.

10. A system for identifying gas zones based on array sonic time-frequency analysis for carrying out the method of any one of claims 1 to 9, characterized in that it comprises: an acquisition module configured to acquire array acoustic full wave train data of a target well; a denoising module configured to perform frequency domain band pass filtering and time domain median filtering on the array acoustic full wave train data of the target well; a decomposition module configured to decompose the denoised array acoustic full wave train data based on a two-dimensional variational mode decomposition method to obtain K mode signal components with different center frequencies; an extraction module configured to extract high-frequency mode components sensitive to formation gas content from the K mode signal components based on an energy-frequency joint screening criterion; and a reconstruction module configured to reconstruct the high-frequency mode components. A time-frequency module is configured to calculate a time-frequency spectrum of the high-frequency modal component based on a smoothed pseudo Wigner-Ville distribution algorithm, and obtain corresponding spectral energy.