Seismic data high-resolution time-frequency analysis characterization method based on inverse filtering

Abstract

The invention belongs to the technical field of seismic data high-resolution time-frequency analysis and characterization, and particularly relates to a seismic data high-resolution time-frequency analysis and characterization method based on inverse filtering, and the method comprises the steps: calculating the time-varying dominant frequency of a seismic signal, and obtaining a W transformation result of the seismic signal based on the time-varying dominant frequency; parameters of the sliding time window are set, the parameters comprise the time window length and the sliding step length, and the sliding step length is integer power of 2; the sliding time window is applied to the W transformation result, and inverse filtering operation is carried out on the W transformation result in the sliding time window; and moving the sliding time window according to a set sliding step length, and repeating the step 3 until the sliding time window traverses the whole seismic data to obtain a high-resolution time-frequency analysis characterization result.

Description

Technical Field

[0001] This invention belongs to the field of high-resolution time-frequency analysis and characterization technology of seismic data, and particularly relates to a high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering. Background Technology

[0002] In the field of seismic exploration, high-resolution time-frequency analysis is a key technology for characterizing subsurface geological bodies and identifying oil and gas reservoirs and thin layers. Among these methods, linear time-frequency analysis is widely used due to its intuitive principle and high computational efficiency, with the W transform being a typical example. The W transform defines a non-steady-state window function by introducing a time-varying dominant frequency, aiming to balance time resolution across high and low frequency bands. However, due to the piecewise processing characteristics of the window function, its time-frequency results still exhibit energy aliasing, especially in thin-layer identification scenarios, which can easily lead to mutual interference between the time-frequency energies of adjacent strata, making it difficult to accurately distinguish between thin-layer boundaries and thicknesses. Furthermore, traditional linear time-frequency analysis methods are generally constrained by the Gabor uncertainty principle, resulting in a trade-off between time resolution and frequency resolution. Pursuing higher time resolution to identify thin layers often sacrifices frequency resolution, and vice versa, which severely limits the ability to finely characterize complex subsurface geological structures.

[0003] While nonlinear time-frequency analysis methods can improve resolution to some extent, they are prone to cross-term interference and are quite sensitive to noise, resulting in insufficient stability in practical seismic data processing. Current seismic exploration demands increasingly higher data resolution, especially in challenging areas such as nearshore and complex tectonic zones, where existing time-frequency analysis methods struggle to meet the requirements for precise reservoir location and detailed thin-layer characterization. Summary of the Invention

[0004] The purpose of this invention is to address the aforementioned technical problems by providing a high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering.

[0005] In view of this, the present invention provides a high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering, comprising the following steps: Step 1: Calculate the time-varying dominant frequency of the seismic signal, and obtain the W-transform result of the seismic signal based on the time-varying dominant frequency; Step 2: Set the parameters of the sliding window, including the window length and the sliding step size, wherein the sliding step size is an integer power of 2; Step 3: Apply the sliding window to the W-transform result and perform inverse filtering on the W-transform result within the sliding window; Step 4: Move the sliding time window according to the set sliding step size, and repeat Step 3 until the sliding time window has traversed the entire seismic data to obtain high-resolution time-frequency analysis characterization results.

[0006] Preferably, the calculation method for the time-varying dominant frequency of the seismic signal in step one is as follows: Based on the instantaneous frequency and amplitude of the seismic record, the following formula is calculated after filtering with a Gaussian low-pass filter and combining it with a scaling parameter: The specific expression is as follows: ; in, The instantaneous frequency of the earthquake record; The instantaneous amplitude recorded in the earthquake data; It is a Gaussian low-pass filter; r is a scaling parameter.

[0007] Preferably, the expression for calculating the W transformation result in step one is: ; in, For time; This is the window function corresponding to the transformation; t is the time at the center of the window function; f is the frequency; k is a scale parameter that controls the window width; This is the time-varying main frequency calculated in step one.

[0008] Preferably, in step two, the length of the sliding window and the sliding step size are determined according to the characteristics of the seismic data and the actual processing and interpretation requirements, and the sliding step size is 8 or 16.

[0009] Preferably, the inverse filtering operation in step three employs a two-dimensional deconvolution algorithm, specifically including: Search for the maximum value of the W-transform result within the window during the glide, calculate the WVD result of the W-transform window function at the maximum value, and use it as the point spread function (PSF). The square of the W-transform result within the sliding window is used as the initial value for iteration. The PSF is used as the processing operator for iterative operation. The number of iterations does not exceed 30. After the iteration terminates, the inverse filtering result within the sliding window is obtained.

[0010] Preferably, the expression for calculating the WVD result of the W transform window function is: ; Where u represents the time delay and f represents the frequency. It is a conjugate operator.

[0011] Preferably, the scale parameter k for controlling the window width in step one is set to 1.

[0012] Preferably, the number of iterations is set according to the actual resolution requirements, and the value is 5-10 times.

[0013] Preferably, the sliding window is a time window that moves in the time direction, used to perform piecewise inverse filtering on the W-transform result of the non-steady-state features.

[0014] Preferably, the high-resolution time-frequency analysis characterization results are used for oil and gas reservoir distribution characterization and thin-layer identification, wherein the thin-layer identification includes characterizing the morphology and distribution features of gas-bearing thin layers.

[0015] The beneficial effects of this invention are: By deeply integrating inverse filtering technology with W-transform, this method effectively overcomes the key limitations of traditional linear time-frequency analysis: it fully preserves the W-transform's ability to accurately characterize geological features in the low-frequency band, while significantly reducing energy aliasing in the time-frequency results through a sliding time window segmented inverse filtering design, making the boundaries and morphological features of underground strata more clearly discernible. Simultaneously, this method successfully alleviates the resolution constraint imposed by the Gabor uncertainty principle, achieving improved temporal resolution without sacrificing frequency resolution. It completely breaks the traditional dilemma of trade-offs between these two aspects, providing a superior solution for identifying fine geological bodies such as thin layers and thin interbedded layers.

[0016] In practical applications, this method demonstrates strong adaptability and practicality: it not only functions reliably in complex structural zones and near-shore areas where exploration is challenging, accurately characterizing time-frequency anomalies related to oil and gas reservoirs, but also provides more reliable technical support for drilling deployment through clear depiction of subsurface geological morphology. Compared to traditional time-frequency analysis methods, its results more intuitively reflect the true state of subsurface geological structures, effectively reducing exploration misjudgments caused by data analysis biases, and further enhancing the application value of seismic data in oil and gas resource evaluation and geological research. Attached Figure Description

[0017] Figure 1 This is a flowchart of the algorithm of the present invention; Figure 2 This is a single-channel synthesis recording example diagram of the present invention; Figure 3 This is a graph showing the time-frequency analysis results of the wedge model of this invention; Figure 4 This is a graph showing the time-frequency analysis results of actual earthquake data from this invention. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0019] It should be noted that all directional and positional terms used in this invention, such as "up," "down," "left," "right," "front," "back," "vertical," "horizontal," "inner," "outer," "top," "lower," "lateral," "longitudinal," and "center," are only used to explain the relative positional relationships and connections between components in a specific state (as shown in the accompanying drawings). They are merely for the convenience of describing the invention and do not require the invention to be constructed and operated in a specific orientation; therefore, they should not be construed as limitations on the invention. Furthermore, descriptions involving "first," "second," etc., are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated.

[0020] In the description of this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal communication between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0021] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "illustrative embodiment," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0022] A high-resolution time-frequency analysis method based on W-transform is proposed. By performing inverse filtering on the W-transform results, the limitation of the window function on the resolution is reduced, resulting in time-frequency results with less energy aliasing. This improves the characterization ability of oil and gas reservoir distribution and reduces the difficulty of thin-layer identification.

[0023] Because the W transform uses a window function to segment the signal, this type of time-frequency analysis method is called a linear time-frequency analysis method. For a seismic signal... The functional expressions of linear time-frequency analysis methods can generally be summarized into the following form, i.e. ; in For time, Let be the window function of the corresponding transform, t be the time at the center of the window function, and f be the frequency. To achieve ideal time resolution across both high and low frequency bands, the W transform defines the window function as a non-stationary signal form, i.e. ; Where k is a scale parameter that controls the window width, and is usually set to 1. For signal The time-varying dominant frequency is calculated by weighting the instantaneous frequencies of the seismic data. The specific expression is: ; in, The instantaneous frequency of the earthquake record. The instantaneous amplitude of the earthquake record. Here, is a Gaussian low-pass filter, and r is the scaling parameter. The seismic signal can then be derived from formula – The specific expression for the W transform is as follows: ; Besides the linear methods mentioned above, nonlinear time-frequency analysis is also frequently used in seismic exploration. The Wigner-Ville distribution (WVD) is one of the most typical nonlinear methods. (Regarding seismic signals...) Its WVD result is ; Where u represents the time delay and f represents the frequency. This is a conjugate operator. Similarly, we can also calculate the WVD result of the W-transform window function, i.e. ; Combining the formulas, we can represent the two-dimensional convolution of the signal's WVD and the window function's WVD using the W transform, i.e.: ; As can be seen from equation (7), the linear time-frequency analysis method can be characterized as a two-dimensional filter of the seismic data time spectrum by a Gaussian time window function. Therefore, we can use an inverse filtering algorithm to perform inverse filtering on the W-transform time window.

[0024] Furthermore, due to the non-steady-state characteristics of the W-transform window function (the instantaneous frequency in the W-transform changes with time), this invention employs a sliding window during inverse filtering. That is, a sliding window that moves in the time direction is set on the W-transform result, and inverse filtering is performed sequentially on the W-transform within the sliding window as it moves. The flowchart is as follows. Figure 1 As shown, the specific calculation scheme is as follows: Calculating the W-transform result: To calculate the W-transform result, we first use the formula to calculate the time-varying dominant frequency of the signal. Then, the W transformation result is calculated using the formula.

[0025] Setting the parameters for the sliding window: The parameters for the sliding window are the window length and the sliding step size. The length of the sliding window and the step size for each movement are usually set according to the seismic data and the actual processing and interpretation requirements. The step size is generally an integer power of 2, such as 8 or 16.

[0026] Inverse filtering is performed on the W transform results within the sliding window: This invention employs a two-dimensional deconvolution algorithm to effectively separate the seismic data time spectrum from the Gaussian time window PSF in linear time-frequency analysis. Specifically, the operation involves searching within the sliding window... maximum value ,calculate The WVD result of the W transform window function is used as the point spread function (PSF). Then, using the square of the W transform result within the sliding window as the initial value, iterative calculations are performed using the PSF as the processing operator until the set number of iterations is met. The final iterative result is the desired result within the sliding window.

[0027] Moving window: After calculating the result of the inverse filtering within the sliding window, move the sliding window to a position according to the designed step size and repeat step 3 until the sliding window has traversed the entire seismic data.

[0028] To demonstrate the process of this invention and its effectiveness in improving resolution, we performed time-frequency analysis on a simple synthetic seismic record. Figure 2 'a' represents a synthetic seismic signal composed of a series of Ricker waves, with the dominant frequencies of the Ricker waves from top to bottom being [60, 50, 40, 30, 20] Hz. Figure 2 c and Figure 2 d represents the S-transform and W-transform results of the synthesized signal, respectively. In this example, we set the sliding window length to 64ms and the step size to 8ms. Figure 1 The process shown calculates the inverse filtering operator PSF ( Figure 2 b), respectively obtained 5 iterations ( Figure 2 e) and 10 times ( Figure 2 f) The result after processing; from Figure 2 As can be seen, this invention exhibits high time resolution similar to the W transform in the low-frequency band. Simultaneously, this invention improves the time resolution of the W transform results and reduces energy aliasing. With increasing iteration count, the time resolution of the results continuously improves. In particular, the degree of energy cluster aliasing at the dominant frequency of each wavelet is significantly reduced. Comparing the results before and after processing reveals that their frequency resolutions are similar. This indicates that this invention alleviates the resolution limitations imposed by the Gabor uncertainty principle, achieving higher time-frequency results without sacrificing frequency resolution, avoiding trade-offs between resolutions, and thus better enabling tasks such as thin-layer identification. Regarding the iteration count, it is generally set according to actual needs and objectives, typically not exceeding 30 iterations.

[0029] To further verify the effectiveness of this invention, a wedge-shaped model was designed and subjected to time-frequency analysis. Figure 3 'a' represents the velocity model, where the red portion represents the formation velocity of 3500 m / s, and the blue portion represents the formation velocity of 2800 m / s. The time thickness of the thickest part of the red formation is 100 ms, and the time thickness of the thinnest part is 50 ms. The formation density model is derived from an empirical formula. It is concluded that, among them Indicates speed, The density is calculated as follows: There are a total of 100 channels, with a channel spacing of 10 meters, 175 sampling points per channel, and a sampling interval of 1 ms. Figure 3 b represents the simulated seismic record synthesized using a convolution model. The wavelet used is a 25Hz zero-phase Ricker wavelet. In this example, the sliding window width is set to 64, the step size to 8, and the number of iterations to 10. Figure 3 c and Figure 3 d represents the 5Hz profile obtained by W transformation and estimation by the present invention, respectively. Figure 3 e and Figure 3 f represents the 45Hz profile obtained by W transformation and estimation by the present invention, respectively.

[0030] from Figure 3 The comparison of the results shows that the present invention has a reliable and significant effect in reducing energy cluster aliasing and improving temporal resolution. Whether in high-frequency or low-frequency results, the present invention estimates energy clusters with smaller temporal thickness and less aliasing, thus better separating the two energy axes. A comparison of the W-transform and the results of the present invention at low frequencies (…) Figure 3 c and Figure 3d) It can be intuitively seen that this invention inherits the high-resolution characteristics of W-transform in the low-frequency range, while further improving the resolution of W-transform results at low frequencies. This indicates that this invention can not only perform the task of W-transform in characterizing the distribution characteristics of underground gas-bearing reservoirs, but also more clearly display the target geological body (such as the morphology and distribution characteristics of gas-bearing thin layers), guiding the deployment of well drilling operations.

[0031] To verify the effectiveness of this invention in applying actual seismic data, we conducted time-frequency analysis experiments on actual seismic profiles acquired by the OGA Rockall Trough project. This dataset consists of near-shore data from the Atlantic coast of the United Kingdom. The work area comprises 6550... The study consists of 5300 survey lines, each with a time range of 0–10 s and a sampling time of 4 ms. We performed time-frequency analysis on the superimposed profile of survey line 42, whose dominant frequency is around 15 Hz. We used the W transform and our invention to perform time-frequency analysis on this profile, setting the sliding window width to 512, the step size to 256, and the number of iterations to 5. To enhance the clarity of the display, we magnified some local areas in the estimated results; see details below. Figure 4 As shown. Among them Figure 4 a represents the magnified seismic profile. Figure 4 b and Figure 4 c represents the 5Hz profile obtained by W-transformation and estimation based on the present invention. Figure 4 d and Figure 4 e represents the 30Hz profile obtained by W transformation and estimation based on the present invention.

[0032] from Figure 4 As can be seen, this invention can obtain higher resolution time-frequency results in actual seismic records. Especially for... Figure 4 The area enclosed by the red ellipse in the middle is difficult for the W transform to clearly characterize the thin-layer properties. However, this invention reduces energy aliasing in the W transform results, successfully depicting the morphological characteristics of the thin layer. Furthermore, by comparing the high-frequency results of the W transform (… Figure 4 d) and low-frequency results ( Figure 4 (b) It can be observed that high-frequency energy is missing within the red rectangle, while obvious energy clusters appear in the low-frequency results. This phenomenon is often a characteristic of hydrocarbon reservoir distribution. Figure 4 As can be seen, the results of this invention not only show the characteristics of oil and gas reservoir distribution at the red rectangle, but also more clearly show the morphological characteristics of gas-bearing strata.

[0033] In summary, this invention introduces an inverse filtering method into the W-transform, which not only inherits the excellent low-frequency resolution of the W-transform but also achieves results with less energy aliasing and higher resolution. More importantly, this invention alleviates the resolution limitations imposed by the Gabor uncertainty principle, improving temporal resolution without reducing frequency resolution. Due to these advantages, this invention can better characterize the distribution characteristics of oil and gas reservoirs and more clearly characterize thin-layer morphology in actual seismic data interpretation, thereby better guiding well drilling operations.

[0034] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of this application.

Claims

1. A high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering, characterized in that: Includes the following steps: Step 1: Calculate the time-varying dominant frequency of the seismic signal, and obtain the W-transform result of the seismic signal based on the time-varying dominant frequency; Step 2: Set the parameters of the sliding window, including the window length and the sliding step size, wherein the sliding step size is an integer power of 2; Step 3: Apply the sliding window to the W-transform result and perform inverse filtering on the W-transform result within the sliding window; Step 4: Move the sliding time window according to the set sliding step size, and repeat Step 3 until the sliding time window has traversed the entire seismic data to obtain high-resolution time-frequency analysis characterization results.

2. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: The calculation method for the time-varying dominant frequency of the seismic signal in step one is as follows: Based on the instantaneous frequency and amplitude of the seismic record, the following formula is calculated after filtering with a Gaussian low-pass filter and combining it with a scaling parameter: The specific expression is as follows: ； in, The instantaneous frequency of the earthquake record; The instantaneous amplitude recorded in the earthquake data; It is a Gaussian low-pass filter; r is a scaling parameter.

3. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 2, characterized in that: The expression for calculating the W-transform result in step one is: ; in, For time; This is the window function corresponding to the transformation; t is the time at the center of the window function; f is the frequency; k is a scale parameter that controls the window width; This is the time-varying main frequency calculated in step one.

4. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: In step two, the length of the sliding window and the sliding step size are determined based on the characteristics of the seismic data and the actual processing and interpretation requirements. The sliding step size is either 8 or 16.

5. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: Step three, the inverse filtering operation, employs a two-dimensional deconvolution algorithm, specifically including: Search for the maximum value of the W-transform result within the window during the glide, calculate the WVD result of the W-transform window function at the maximum value, and use it as the point spread function (PSF). The square of the W-transform result within the sliding window is used as the initial value for iteration. The PSF is used as the processing operator for iterative operation. The number of iterations does not exceed 30. After the iteration terminates, the inverse filtering result within the sliding window is obtained.

6. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 5, characterized in that: The expression for calculating the WVD result of the W-transform window function is: ； Where u represents the time delay and f represents the frequency. It is a conjugate operator.

7. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: In step one, the scale parameter k for controlling the window width is set to 1.

8. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 5, characterized in that: The number of iterations is set according to the actual resolution requirements, ranging from 5 to 10.

9. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: The sliding window is a time window that moves in the time direction and is used to perform piecewise inverse filtering on the W-transform results of non-steady-state features.

10. The high-resolution time-frequency analysis and characterization method for seismic data based on inverse filtering according to claim 1, characterized in that: The high-resolution time-frequency analysis results are used for oil and gas reservoir distribution characterization and thin-layer identification, the thin-layer identification including the characterization of the morphology and distribution features of gas-bearing thin layers.

Images