A seismic reservoir identification method and system based on sparse Gabor transform

Abstract

The present invention relates to the technical field of seismic reservoir identification, and in particular to a seismic reservoir identification method and system based on sparse Gabor transform, the method comprising: obtaining seismic gather data, superimposing the seismic gather data, obtaining full stacking data and angle stacking data; performing sparse Gabor transform focusing on time resolution on the full stacking data, obtaining a time-frequency spectrum with high time resolution; performing sparse Gabor transform focusing on frequency resolution on the angle stacking data, obtaining a time-frequency spectrum with high frequency resolution; performing two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution, obtaining a time-frequency spectrum after two-dimensional deconvolution; and identifying the seismic reservoir according to the time-frequency spectrum after two-dimensional deconvolution, obtaining the range of the seismic reservoir. The present invention proposes that sparse Gabor transform can effectively improve the time-frequency focusing, thereby improving the analysis accuracy of seismic data.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic reservoir identification, and in particular to a seismic reservoir identification method and system based on sparse Gabor transform. Background Art

[0002] The existing Gabor transform has become one of the important time-frequency analysis tools in the field of signal analysis and processing. The transform coefficient can well reflect the degree of change of the signal frequency over time, and is widely used in speech processing, image processing and communication systems. However, the Gabor transform is limited by the uncertainty principle, and the time and frequency resolution cannot reach the highest at the same time. As a result, the Gabor transform cannot provide high-resolution time-frequency spectrum in time and frequency at the same time in seismic signal analysis, which leads to limited recognition of thin layers or complex geological structures. First, the size and shape of the time-frequency window of the Gabor transform remain unchanged, only the position changes, while in practical applications, it is often hoped that the size and shape of the time-frequency window will change with the change of frequency. Secondly, the window length of the Gabor transform is fixed, and low-frequency or high-frequency data will be damaged during the time-frequency domain processing. In addition, the Gabor transform adds a window function to the signal to be analyzed, which changes the nature of the original signal.

[0003] In summary, the existing Gabor transform is limited by the uncertainty principle, and the time and frequency resolution cannot reach the highest at the same time. In addition, Gabor is non-orthogonal, and there is redundancy between different feature components, so it is not very efficient in analyzing texture images. In seismic reservoir identification, there are problems such as insufficient time and frequency resolution, noise interference, limited transient signal analysis capabilities, and low dispersion attribute extraction accuracy. Summary of the invention

[0004] The purpose of the present invention is to provide a seismic reservoir identification method and system based on sparse Gabor transform, which can accurately locate the reservoir.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A seismic reservoir identification method based on sparse Gabor transform, comprising:

[0007] S1, obtaining seismic gather data, performing full stacking and angle stacking on the seismic gather data, and obtaining full stacking data and angle stacking data;

[0008] S2, performing a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function;

[0009] S3, performing a sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution;

[0010] S4, performing two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution, if all seismic traces are not completed, return to S2, if all seismic traces are completed, enter S5;

[0011] S5. Identify the seismic reservoir according to the time-frequency spectrum after the two-dimensional deconvolution to obtain the range of the seismic reservoir, wherein the time-frequency spectrum after the two-dimensional deconvolution includes a time-frequency spectrum with high time resolution after the two-dimensional deconvolution and a time-frequency spectrum with high frequency resolution after the two-dimensional deconvolution.

[0012] Optionally, performing angle stacking on the seismic gather data to obtain the angle stacking data includes:

[0013] An angle range is preset for the seismic gather data, and the seismic gather data is angle-stacked according to the preset angle range to obtain the angle-stacked data.

[0014] Optionally, improving the Gabor transform by using a point spread function to obtain the sparse Gabor transform includes:

[0015] According to the point spread function and linear convolution, a relationship expression between the high-frequency Gabor time-frequency spectrum and the low-frequency Gabor time-frequency spectrum is obtained;

[0016] Converting the relational expression into an expression in the form of matrix multiplication;

[0017] According to the Fourier transform property, the time domain convolution is equal to the frequency domain product, let represents discrete Fourier transform;

[0018] The Combined with the expression in the form of matrix multiplication, an initial sparse Gabor transform is obtained;

[0019] The initial sparse Gabor transform is constrained by using an L1 norm to obtain the sparse Gabor transform.

[0020] Optionally, the sparse Gabor transform is:

[0021] ;

[0022] in, is a sparse high-resolution time-frequency spectrum, The time-spectrum representing the minimum high-frequency focality, , , , To pass the two-dimensional kernel function The point spread function obtained is, is the absolute value of the time-frequency spectrum through high time-frequency focusing The obtained matrix, is the absolute value of the Gabor spectrum through low-frequency focusing The obtained matrix, and is the regularization parameter, Indicates frequency, Indicates time.

[0023] Optionally, the point spread function is obtained by respectively calculating the Dirac Functions and Perform Gabor transform, obtain the time-frequency spectrum and multiply to obtain, where, represents the instantaneous value of the signal, Indicates the signal frequency;

[0024] The point spread function is:

[0025] ;

[0026] in, is the time window width, is the frequency window width.

[0027] Optionally, performing a sparse Gabor transform focusing on time resolution on the full superimposed data includes: presetting the time window width, and performing a sparse Gabor transform focusing on time resolution on the full superimposed data according to the time window width;

[0028] Performing a sparse Gabor transform focusing on frequency resolution on the angle-stacked data includes: presetting the frequency window width, and performing a sparse Gabor transform focusing on time resolution on the full-stacked data according to the frequency window width.

[0029] Optionally, identifying the seismic reservoir according to the time-frequency spectrum after the two-dimensional deconvolution includes:

[0030] Performing frequency decomposition on the high-time-resolution time-frequency spectrum after the two-dimensional deconvolution to obtain a thin layer identification result;

[0031] Solving the dispersion equation for the high-frequency resolution time-frequency spectrum after the two-dimensional deconvolution to obtain the dispersion AVO property;

[0032] The thin layer identification result and the dispersion AVO attribute are superimposed and smoothed to obtain the seismic reservoir range.

[0033] On the other hand, the present invention also provides a seismic reservoir identification system based on sparse Gabor transform, a data acquisition module, a first data analysis module, a second data analysis module, a deconvolution module and a seismic reservoir identification module;

[0034] The data acquisition module is used to acquire seismic gather data, perform full stacking and angle stacking on the seismic gather data, and acquire full stacking data and angle stacking data;

[0035] The first data analysis module is used to perform a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function;

[0036] The second data analysis module is used to perform a sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution;

[0037] The deconvolution module is used to perform two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution, and if all seismic traces are not completed, return to the first data analysis module and the second data analysis module; if all seismic traces are completed, enter the seismic reservoir identification module;

[0038] The seismic reservoir identification module is used to identify the seismic reservoir according to the time-frequency spectrum after two-dimensional deconvolution and obtain the range of the seismic reservoir, wherein the time-frequency spectrum after two-dimensional deconvolution includes a time-frequency spectrum with high time resolution after two-dimensional deconvolution and a time-frequency spectrum with high frequency resolution after two-dimensional deconvolution.

[0039] The beneficial effects of the present invention are as follows: the present invention aims at the problem of low time-frequency resolution of Gabor transform, proposes sparse Gabor transform, obtains time-frequency spectrum based on Gabor transform, and then obtains high-resolution time-frequency spectrum by sparse inversion, which can effectively improve time-frequency focusing; the present invention aims at the influence of time-frequency window function, adds sparse constraint to time-frequency spectrum by compressed sensing theory, solves two-dimensional deconvolution of Gabor transform time-frequency spectrum by ADMM algorithm, and can remove the influence of time-frequency window function; the present invention aims at the characteristics of fixed time-frequency window size and shape of Gabor transform, proposes an adjustable sparse Gabor transform, when the time window of point spread function is short and the frequency window is long, a time-frequency spectrum with high time resolution can be obtained; the present invention can weaken the limitation of uncertainty principle of Gabor transform, and improve the analysis accuracy of seismic data. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0041] Figure 1 A flowchart of a seismic reservoir identification method based on sparse Gabor transform according to an embodiment of the present invention;

[0042] Figure 2 Schematic diagram of the calculation principle of TSGT and FSGT according to an embodiment of the present invention;

[0043] Figure 3 A schematic diagram of a PSF function used for deconvolution in an embodiment of the present invention;

[0044] Figure 4 Schematic diagrams of time-frequency analysis of a synthetic model of an embodiment of the present invention and different algorithms for the synthetic model, wherein (a) is a synthetic model, (b) is a schematic diagram of time-frequency analysis of the synthetic model by the GT algorithm, (c) is a schematic diagram of time-frequency analysis of the synthetic model by the ST algorithm, (d) is a schematic diagram of time-frequency analysis of the synthetic model by the SST algorithm, (e) is a schematic diagram of time-frequency analysis of the synthetic model by the TSET algorithm, and (f) is a schematic diagram of time-frequency analysis of the synthetic model by the TSGT algorithm;

[0045] Figure 5 Schematic diagram of the Marmousi reflection coefficient model and synthetic seismic record according to an embodiment of the present invention, wherein (a) is the Marmousi reflection coefficient model, and (b) is the synthetic seismic record;

[0046] Figure 6 Schematic diagram of the frequency division cross section of the Marmousi reflection coefficient according to an embodiment of the present invention, wherein (a) is the GT frequency division cross section, (b) is the TSET frequency division cross section, and (c) is the TSET frequency division cross section;

[0047] Figure 7 Schematic diagram of actual data and frequency division profiles of an embodiment of the present invention, wherein (a) is a GT frequency division profile, (b) is a TSET frequency division profile, and (c) is a TSGT frequency division profile;

[0048] Figure 8 A PSF function used for deconvolution in an embodiment of the present invention;

[0049] Figure 9 A schematic diagram of a frequency modulation signal for a synthetic test according to an embodiment of the present invention;

[0050] Figure 10Schematic diagram of time-frequency analysis of a frequency modulation signal according to an embodiment of the present invention, wherein (a) is the time-frequency analysis of the real instantaneous frequency, (b) is the time-frequency analysis of GT, (c) is the time-frequency analysis of ST, (d) is the time-frequency analysis of SST, (e) is the time-frequency analysis of L12-STFT, and (f) is the time-frequency analysis of FSGT;

[0051] Figure 11 The figures are the overlay graphs of the dispersion AVO attributes and thin layer identification results extracted by different time-frequency analysis methods in the embodiments of the present invention, wherein (a) is the overlay graph of the GT dispersion AVO attributes and thin layer identification, (b) is the overlay graph of the SST dispersion AVO attributes and thin layer identification, (c) is the overlay graph of the L12-STFT dispersion AVO attributes and thin layer identification, and (d) is the overlay graph of the FSGT dispersion AVO attributes and thin layer identification. DETAILED DESCRIPTION

[0052] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0053] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0054] Embodiment 1:

[0055] like Figure 1 As shown, this embodiment provides a seismic reservoir identification method based on sparse Gabor transform, including:

[0056] S1. Acquire seismic gather data, perform full stacking and angle stacking on the seismic gather data, and acquire full stacking data and angle stacking data;

[0057] S2, performing a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function;

[0058] S3, performing sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution;

[0059] S4, performing two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution, if all seismic traces are not completed, returning to S2, if all seismic traces are completed, entering S5;

[0060] S5. Identify the seismic reservoir according to the time-frequency spectrum after two-dimensional deconvolution to obtain the range of the seismic reservoir, wherein the time-frequency spectrum after two-dimensional deconvolution includes a time-frequency spectrum with high time resolution after two-dimensional deconvolution and a time-frequency spectrum with high frequency resolution after two-dimensional deconvolution.

[0061] Furthermore, the seismic gather data are angle stacked to obtain angle stacking data, including:

[0062] An angle range is preset for the seismic gather data, and the seismic gather data is angle-stacked according to the preset angle range to obtain the angle-stacked data.

[0063] Specifically, all the traces of the seismic gather are first superimposed to form full stack data, also called stacked profiles. At the same time, the angle range of the seismic gather is manually set, and the traces are superimposed at small, medium, and large angles to form partial stacked profiles. Among them, the angle division is based on the full angle range of the seismic data, which is evenly divided into three angles.

[0064] Furthermore, the Gabor transform is improved by the point spread function to obtain the sparse Gabor transform, including:

[0065] According to the point spread function and linear convolution, the relationship expression between the high-frequency Gabor time-frequency spectrum and the low-frequency Gabor time-frequency spectrum is obtained;

[0066] Convert relational expressions into expressions in matrix multiplication form;

[0067] According to the Fourier transform property, the time domain convolution is equal to the frequency domain product, let represents discrete Fourier transform;

[0068] Will Combined with the expression in matrix multiplication form, the initial sparse Gabor transform is obtained;

[0069] The initial sparse Gabor transform is constrained by using the L1 norm to obtain the sparse Gabor transform.

[0070] Furthermore, the sparse Gabor transform is shown in formula (6).

[0071] Specifically, obtaining the sparse Gabor transform includes:

[0072] Image deblurring and image restoration are important applications of compressed sensing in the field of image processing. For Gabor transform spectrum, we hope to remove the influence of window function through image deblurring. Then we need to use two-dimensional deconvolution to obtain a highly focused time-frequency spectrum when the point spread function is known.

[0073] In this embodiment, the two-dimensional deconvolution of the spectrum during Gabor transformation is solved by using the ADMM algorithm.

[0074] The optimization problem is described as: the absolute value of the Gabor time-frequency spectrum with known low time-frequency focality Under the condition of other prior knowledge, the absolute value of the time-frequency spectrum with high time-frequency focusing is obtained In general, and There is the following linear degradation relationship, that is, the relationship between the high-frequency Gabor time-frequency spectrum and the low-frequency Gabor time-frequency spectrum is expressed as:

[0075] (1);

[0076] Among them, h(f,t) is a two-dimensional kernel function; is the absolute value of the time-frequency spectrum of the original signal noise; Represents a two-dimensional linear convolution. It is further processed by matrix operations (such as convolution, deconvolution, sparse constraints, etc.) , , and ,make , , and ,in, , , and The matrices are , , and At a certain point, To pass the two-dimensional kernel function The point spread function obtained is, is the absolute value of the time-frequency spectrum through high time-frequency focusing The obtained matrix, is the absolute value of the Gabor spectrum through low-frequency focusing The obtained matrix, is the absolute value of the time-frequency spectrum of the original signal noise The obtained matrix, the position of the point is determined by the frequency and time Determine that formula (1) can be written as the following matrix multiplication form:

[0077] (2);

[0078] According to the Fourier transform property, the time domain convolution is equal to the frequency domain product, so that the matrix represents the discrete Fourier transform, that is:

[0079] (3);

[0080] in, It is called the butterfly factor, which has symmetry and periodicity, and N represents the number of matrices.

[0081] (4);

[0082] The high-resolution time-frequency spectrum can be obtained by optimizing equation (5), where: represents circular convolution, Representation Matrix Complex conjugate transpose, when the kernel function is estimated After that, the deconvolution required It is called the point spread function (PSF).

[0083] The high-resolution time-frequency spectrum, i.e. the initial sparse Gabor transform, is:

[0084] (5);

[0085] This embodiment hopes to obtain a sparse high-resolution time-frequency spectrum , with L 1 Norm constraint (5), where The time-frequency spectrum represents the minimum high-frequency focus while adding the positivity constraint on the amplitude spectrum, i.e.,

[0086] (6);

[0087] In the formula, and represents the regularization parameter.

[0088] From formula (6), we know that the high-resolution time-frequency spectrum Mainly determined by the regularization parameter and the point spread function Control, obviously, The larger it is, the sparser the amplitude spectrum is, and the more weak energy is lost; The smaller the value, the lower the resolution of the amplitude spectrum. The PSF required for deconvolution is the kernel function of the Gabor transform. This embodiment requires a reasonable estimation of the shape of the PSF on the amplitude spectrum.

[0089] Furthermore, the point spread function is obtained by respectively calculating the Dirac Functions and Perform Gabor transform, where represents the instantaneous value of the signal, Represents the signal frequency, obtains the time-frequency spectrum and multiplies it; the point spread function is shown in formula (7).

[0090] Specifically, the kernel function of the Gabor transform is calculated by the variance By controlling the time window and frequency window length at the same time, the PSF can be directly obtained based on equation (7), but in fact, the kernel function can be obtained in an easier way , that is, respectively for Dirac Functions and Perform Gabor transform, and the product of the two time spectra is the estimated kernel function ,in ,in, Represents the time sampling interval to ensure that the PSF is symmetric about the center. Although the kernel function required to obtain a high-resolution time-frequency spectrum is obvious, due to the characteristics of the signal itself, the direct application of the two-dimensional kernel function of the Gabor forward transform cannot obtain an accurate high-resolution time-frequency spectrum, such as FM signals. This type of signal has a certain bandwidth and usually requires a larger frequency window to protect these frequency components, such as time-synchronous squeezing transform and time-synchronous extraction transform; when the signal is a completely stable simple harmonic signal, the resolution in the time direction does not need to be improved, that is, the time window width needs to be increased to protect the stability of the signal, such as synchronous squeezing transform, Hilbert-Huang transform (HHT), etc.

[0091] The two-dimensional Gaussian kernel function is:

[0092] (7);

[0093] In the formula, It represents the variance. From the form of formula (8), it can be seen that the time window length and the frequency window length restrict each other.

[0094] Therefore, considering the problem that the kernel function parameters of Gabor transform are not independent, equation (7) is decomposed into a two-dimensional window with independent variance parameters, which is PSF:

[0095] (8)

[0096] in, Indicates the width of the control time window. Indicates controlling the frequency window width.

[0097] Regarding the windowing problem of Gabor transform, its windowing effect is mainly determined by two factors: 1. The time-frequency resolution of the Gabor time-frequency spectrum itself; 2. The size of the PSF used for windowing; Obviously, if high time resolution is required, the time window of the Gabor transform itself needs to be small enough, and then the time window used by the kernel function is also small enough; similarly, if high frequency resolution is required, the frequency window of the kernel function used by the Gabor transform needs to be small enough, and then the frequency window of the PSF is also small enough; in this way, one parameter required for the PSF (time window length or frequency window length) is determined based on the Gabor transform itself, and the other parameter (the corresponding frequency window length or time window length) needs to be determined manually. The selection principle is: if the resolution in another direction is taken into account, the length is smaller than the window length of the Gabor transform itself, otherwise it must be larger than the window length of the Gabor transform itself, but it should not be too small or too large.

[0098] Specifically, according to the time window and frequency window length of the point spread function (PSF) used in formula (8), the resolution of the sparse Gabor transform can be manually controlled. When the time window of the point spread function is short and the frequency window is long, a time-frequency spectrum with high time resolution, namely TSGT, can be obtained; when the time window length and frequency window of the point spread function are short, the sparse Gabor transform has good frequency resolution, namely FSGT. According to formula (5), the control effect of the time window and frequency window length of the PSF on deconvolution is as follows: Figure 2 As shown, when the time window is small (time sparse Gabor transform with high time resolution TSGT), the frequency resolution is infinitely low and the time resolution is improved. Such a time-frequency spectrum no longer has frequency meaning. This approach is similar to the time synchronous squeezing transform and is suitable for transient signal analysis, seismic thin layer identification, etc.; when the frequency window is small (frequency sparse Gabor transform with high frequency resolution FSGT), the data requires frequency resolution and is suitable for extracting frequency-varying AVO attributes (the amplitude spectrum formed is similar to the frequency synchronous squeezing transform).

[0099] Furthermore, performing a sparse Gabor transform focusing on time resolution on the full stack data includes: presetting a time window width, and performing a sparse Gabor transform focusing on time resolution on the full stack data according to the time window width;

[0100] The sparse Gabor transform focusing on frequency resolution of the angle-stacked data includes: presetting a frequency window width, and performing a sparse Gabor transform focusing on time resolution on the full-stacked data according to the frequency window width.

[0101] Specifically, according to the time-frequency window length of the point spread function, the frequency window of the sparse Gabor transform is increased, the time window length is reduced, and the Gabor transform focusing on time resolution is obtained, and the time-frequency spectrum is deconvolved to obtain the seismic profile after frequency decomposition to characterize the reservoir boundary. Similarly, the frequency window length is reduced to obtain the Gabor transform focusing on frequency resolution, and the dispersion equation is solved through the deconvolved time-frequency spectrum to obtain the dispersion AVO attribute that can directly indicate the reservoir.

[0102] Furthermore, based on the time-frequency spectrum after two-dimensional deconvolution, the seismic reservoir is identified by:

[0103] Perform frequency decomposition on the high-resolution time-frequency spectrum after two-dimensional deconvolution to obtain thin layer identification results;

[0104] The dispersion equation is solved for the high-frequency resolution time-frequency spectrum after two-dimensional deconvolution to obtain the dispersion AVO properties;

[0105] The thin layer identification results and the dispersion AVO attributes are superimposed and smoothed to obtain the seismic reservoir range.

[0106] Model testing:

[0107] To facilitate testing, this embodiment provides some time-frequency analysis algorithms and the resolutions emphasized when applying them for analysis, as shown in Table 1.

[0108] Table 1

[0109]

[0110] A. Test of time-frequency analysis method focusing on time resolution:

[0111] This model tests the applicability of several time-frequency analysis methods that focus on temporal resolution to the model. Among them, the sparse Gabor transform is used as follows Figure 3 The PSF shown in the figure has a narrow time window and an infinite frequency window. At this time, the sparse Gabor transform (SGT) changes to TSGT. This time-frequency analysis feature is used to analyze the thin layer of seismic data. The synthetic model is as follows: Figure 4 As shown in (a), this embodiment simulates a three-layer model with reverse polarity, and the layer thickness gradually increases from the shallow layer to the deep layer. The signal-to-noise ratio after adding noise is 11.85dB, as shown in Figure 4(b)-(f) This embodiment uses Gabor transform (GT), S transform (ST), synchronous squeezing transform (SST), time synchronous extraction transform (TSET) and time sparse Gabor transform (TSGT) algorithms to perform time-frequency analysis on this model. It can be found that for GT, since its time window and frequency window lengths are basically the same, the energy cluster is distributed in a circular shape; because ST has multi-resolution characteristics, its time spectrum is a triangular window, and SST squeezes energy to the instantaneous frequency in the frequency direction, so the time resolution is limited; TSET extracts "quasi-instantaneous frequency" in the time direction, which can effectively identify thin layers, but its disadvantage is that it is seriously interfered by noise and has cross-interference items. TSGT can also identify coupled sub-waves, and the strong energy position corresponds to the reflection interface (reflection coefficient), and through sparse inversion, the noise is suppressed by the threshold.

[0112] The Marmousi model is used to test the applicability of the algorithm in complex event conditions, and the stability of the algorithm applied to multi-channel data can also be tested. Figure 5 (a)-(b) are reflection coefficients and synthetic seismic data; Figure 6 (a)-(c) are the results of frequency decomposition using GT, TSET and TSGT, which lose frequency resolution and are equivalent to thin-layer identification. Considering the coupling and strike of events of different frequencies, the data all extract 30Hz spectral decomposition sections. The GT with window function has the lowest resolution but is also the most stable. The section obtained by TSET is close to the reflection coefficient, which improves the identification of geological units such as faults, pinch-outs, and truncation. However, some of its traces have energy crosstalk, resulting in abnormal closure of the event. This is due to the instability of TSET in obtaining "quasi-instantaneous frequency". Although the resolution of TSGT is slightly lower than that of TSET, there is basically no illusion as a whole, which also shows that TSGT has a very high stability in multi-channel data processing.

[0113] Then, for complex actual data, the above method is also applied to thin layer identification. The main frequency of the data is about 20Hz, and the spectrum decomposition section of the dominant frequency of 20Hz is also extracted. Figure 7 (a)-(c) show that both TSET and TSGT can effectively identify the reflection interface (the arrows reflect the interface information). When the matching degree of the well data of several methods is high, TSGT is obviously more stable for the target frame position and can better reflect the formation structure.

[0114] B. Focus on frequency resolution time-frequency analysis test:

[0115] Because the point spread function size of the SGT method is adjustable, when the point spread function is Figure 8As shown in the figure, the SGT time-frequency spectrum will have both high time resolution and frequency resolution, but the frequency resolution is higher. At this time, the frequency resolution of SGT is higher, that is, FSGT, which makes FSGT suitable for time-frequency analysis of continuous multi-component signals and also for frequency-varying AVO analysis of seismic data.

[0116] In this model, the synthesized frequency modulation signal sig is:

[0117] .

[0118] Among them, sig represents the synthesized FM signal, and yes sub-signal.

[0119] The synthetic data is Figure 9 As shown, the real instantaneous frequency is Figure 10 (a) shows the time-frequency analysis of the data. Figure 10 (b)-(f), the test methods are GT, ST, synchronized squeezing transform (SST), inversion time-frequency analysis based on L12 norm constraint (L12-STFT) and frequency sparse Gabor transform (FSGT) with high frequency resolution. It can be seen intuitively that the two windowed time-frequency analysis methods (GT and ST) have the worst time-frequency focusing, while the other three time-frequency analysis methods based on instantaneous frequency have good focusing, but SST has differences in high and low frequency energy; FSGT is basically consistent with SST and L12-STFT, and its time-frequency spectrum reflects the instantaneous frequency of the signal. The adaptability at the intersection is poor, but the overall energy distribution of FSGT is basically consistent with the real one, and no false images appear.

[0120] The absorption and attenuation of underground media leads to dispersion in the process of seismic wave propagation. Dispersion is manifested in that the propagation velocity of seismic waves in the formation and the reflection coefficient of the reflection interface vary with their frequency. This feature can be used to identify gas layers with strong attenuation, which is reflected in the dispersion properties of seismic longitudinal waves, namely the dispersion AVO (FAVO) property. Selecting a time-frequency analysis method with better time-frequency focusing is similar to selecting time-frequency domain deconvolution. If the time-frequency focusing and resolution of the selected time-frequency analysis method are good, it will help to calculate the dispersion parameters more accurately. The well location of this data is at the 294th common depth point gather CDP of the profile. The curve characteristics of 1975-2000ms are low gamma value, low density value and high resistivity value, which shows that it is a sandstone gas reservoir. The overlying 1975ms is a high gamma mudstone cap rock with strong sealing capacity.

[0121] The dispersion AVO attributes were extracted from the angle stacked gathers of the data, and the dispersion AVO attributes (FAVO) and reflection interface superposition diagrams were drawn in combination with the reflection interface identified by TSGT, as shown in the figure. Figure 11As shown in (b)-(d), the dispersion attribute resolutions of the three methods, SST, L12-STFT and FSGT, are all high, with GT having the lowest resolution, as shown in Figure 11 (a) As shown in the box. Combined with the reflection interface obtained by frequency decomposition, only the reservoir identified by FSGT is restricted within the layer by the interface, which also shows the accuracy of the FSGT method. FSGT is more conducive to accurate reservoir identification, while TSGT can clearly depict the reservoir boundary.

[0122] Embodiment 2:

[0123] A seismic reservoir identification system based on sparse Gabor transform, comprising a data acquisition module, a first data analysis module, a second data analysis module, a deconvolution module and a seismic reservoir identification module;

[0124] A data acquisition module is used to acquire seismic gather data, perform full stacking and angle stacking on the seismic gather data, and acquire full stacking data and angle stacking data;

[0125] The first data analysis module is used to perform a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function;

[0126] The second data analysis module is used to perform a sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution;

[0127] A deconvolution module is used to perform two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution. If all seismic traces are not completed, the module returns to the first data analysis module and the second data analysis module. If all seismic traces are completed, the module enters the seismic reservoir identification module.

[0128] The seismic reservoir identification module is used to identify the seismic reservoir and obtain the seismic reservoir range according to the time-frequency spectrum after two-dimensional deconvolution, wherein the time-frequency spectrum after two-dimensional deconvolution includes the time-frequency spectrum with high time resolution after two-dimensional deconvolution and the time-frequency spectrum with high frequency resolution after two-dimensional deconvolution.

[0129] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A seismic reservoir identification method based on sparse Gabor transform, characterized in that: include: S1, obtaining seismic gather data, performing full stacking and angle stacking on the seismic gather data, and obtaining full stacking data and angle stacking data; S2, performing a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function; Improving the Gabor transform by using a point spread function to obtain the sparse Gabor transform includes: According to the point spread function and linear convolution, a relationship expression between the high-frequency Gabor time-frequency spectrum and the low-frequency Gabor time-frequency spectrum is obtained; Converting the relational expression into an expression in the form of matrix multiplication; According to the Fourier transform property, the time domain convolution is equal to the frequency domain product, let represents discrete Fourier transform; The Combined with the expression in the form of matrix multiplication, an initial sparse Gabor transform is obtained; Using L1 norm to constrain the initial sparse Gabor transform to obtain the sparse Gabor transform; The sparse Gabor transform is: ； in, is a sparse high-resolution time-frequency spectrum, The time-spectrum representation of the minimum high-frequency focality, , , , To pass the two-dimensional kernel function The point spread function obtained is, is the absolute value of the time-frequency spectrum through high time-frequency focusing The obtained matrix, is the absolute value of the Gabor spectrum through low-frequency focusing The obtained matrix, and is the regularization parameter, Indicates frequency, Indicates time; The point spread function is obtained by respectively Functions and Perform Gabor transform, obtain the time-frequency spectrum and multiply to obtain, where, represents the instantaneous value of the signal, Indicates the signal frequency; The point spread function is: ； in, is the time window width, is the frequency window width; S3, performing a sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution; S4, performing two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution, if all seismic traces are not completed, return to S2, if all seismic traces are completed, enter S5; S5. Identify the seismic reservoir according to the time-frequency spectrum after the two-dimensional deconvolution to obtain the range of the seismic reservoir, wherein the time-frequency spectrum after the two-dimensional deconvolution includes a time-frequency spectrum with high time resolution after the two-dimensional deconvolution and a time-frequency spectrum with high frequency resolution after the two-dimensional deconvolution.

2. The seismic reservoir identification method based on sparse Gabor transform according to claim 1, characterized in that: Performing angle stacking on the seismic gather data, and obtaining the angle stacking data comprises: An angle range is preset for the seismic gather data, and the seismic gather data is angle-stacked according to the preset angle range to obtain the angle-stacked data.

3. The seismic reservoir identification method based on sparse Gabor transform according to claim 1, characterized in that: The step of performing a sparse Gabor transform focusing on time resolution on the full stacking data includes: presetting the time window width, and performing a sparse Gabor transform focusing on time resolution on the full stacking data according to the time window width; Performing a sparse Gabor transform focusing on frequency resolution on the angle-stacked data includes: presetting the frequency window width, and performing a sparse Gabor transform focusing on time resolution on the full-stacked data according to the frequency window width.

4. The seismic reservoir identification method based on sparse Gabor transform according to claim 1, characterized in that: According to the time-frequency spectrum after the two-dimensional deconvolution, identifying the seismic reservoir includes: Performing frequency decomposition on the high-time-resolution time-frequency spectrum after the two-dimensional deconvolution to obtain a thin layer identification result; Solving the dispersion equation for the high-frequency resolution time-frequency spectrum after the two-dimensional deconvolution to obtain the dispersion AVO property; The thin layer identification result and the dispersion AVO attribute are superimposed and smoothed to obtain the seismic reservoir range.

5. A seismic reservoir identification system based on sparse Gabor transform for implementing the seismic reservoir identification method based on sparse Gabor transform as described in any one of claims 1 to 4, characterized in that: A data acquisition module, a first data analysis module, a second data analysis module, a deconvolution module and a seismic reservoir identification module; The data acquisition module is used to acquire seismic gather data, perform full stacking and angle stacking on the seismic gather data, and acquire full stacking data and angle stacking data; The first data analysis module is used to perform a sparse Gabor transform focusing on time resolution on the full stack data to obtain a time-frequency spectrum with high time resolution, wherein the sparse Gabor transform is obtained by improving the Gabor transform through a point spread function; The second data analysis module is used to perform a sparse Gabor transform focusing on frequency resolution on the angle stacking data to obtain a time-frequency spectrum with high frequency resolution; The deconvolution module is used to perform two-dimensional deconvolution processing on the time-frequency spectrum with high time resolution and the time-frequency spectrum with high frequency resolution to obtain the time-frequency spectrum after two-dimensional deconvolution, and if all seismic traces are not completed, return to the first data analysis module and the second data analysis module; if all seismic traces are completed, enter the seismic reservoir identification module; The seismic reservoir identification module is used to identify the seismic reservoir according to the time-frequency spectrum after two-dimensional deconvolution and obtain the range of the seismic reservoir, wherein the time-frequency spectrum after two-dimensional deconvolution includes a time-frequency spectrum with high time resolution after two-dimensional deconvolution and a time-frequency spectrum with high frequency resolution after two-dimensional deconvolution.

Images