A gas content prediction method and system based on hilbert-huang transform time-frequency spectrum

Abstract

The application discloses a gas content prediction method and system based on Hilbert-Huang transform time-frequency spectrum, and belongs to the technical field of oil and gas field exploration, and comprises the following steps: performing empirical mode decomposition on a to-be-analyzed section of each seismic trace to obtain a linear steady signal IMF; performing Hilbert transform on each IMF component to obtain an instantaneous frequency and an envelope; synthesizing the frequency spectrum of all the IMF components to obtain a two-dimensional Hilbert-Huang transform time-frequency spectrum; extracting an instantaneous frequency spectrum at a sample point of a certain seismic trace, recording the highest and the second highest frequencies and the corresponding amplitudes; calculating the spectral area of a trapezoid formed by the two spectrum lines, dividing the highest amplitude frequency, and constructing a gas content factor; and calculating the gas content factors of all the sample points of the to-be-analyzed section to obtain evaluation data of the gas content of a reservoir. The method improves the reliability of gas content prediction through Hilbert transform and construction of the gas content factor.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic interpretation technology for oil and gas field exploration and development, and particularly relates to a gas-bearing property prediction method and system based on Hilbert-Huang transform time-frequency spectrum. BACKGROUND

[0002] In the seismic exploration interpretation work, for the continental clastic rock sedimentary strata, in addition to the conventional structure, reservoir interpretation and prediction, the gas-bearing property prediction of the reservoir is also an important work. In recent years, the natural gas exploration of a certain basin has strengthened the research and review of the shallow tight sandstone gas reservoir, and further characterization of the gas-bearing reservoir in these strata is carried out to implement the exploration potential of each block. Among them, the gas-bearing property evaluation of the reservoir is a difficulty and a key point. The existing seismic technology can better depict the planar morphology, thickness and scale of the channel, but the reliability of the prediction of the gas-bearing property is not satisfactory, which affects the formulation of the well deployment decision. Therefore, it is important to know the gas-bearing property of the reservoir in the area without well control, which is an important decision-making prerequisite for the deployment of subsequent well targets.

[0003] The post-stack data of the earthquake is the most widely used, relatively small in size, and the fastest implemented data type in production application, and a large number of reservoir seismic response feature analysis, channel depiction, and reservoir quantitative evaluation are based on the post-stack seismic data. The existing post-stack seismic amplitude and waveform attribute technology can better depict the planar morphology, thickness and scale of the channel, but the evaluation of the gas-bearing condition of the channel reservoir still needs to rely on the frequency attribute. According to a large amount of theoretical and experimental research, the gas-bearing reservoir has a strong attenuation effect on the high-frequency component of the seismic wave, and on the received seismic signal, the reservoir with good gas-bearing property shows the characteristics of relatively strong low frequency and weak high frequency. Therefore, the existing gas-bearing evaluation mainly relies on the signal analysis technology of the earthquake and the analysis and transformation of the frequency component, such as the instantaneous frequency attribute, the instantaneous Q value, the frequency division attribute, the two-dimensional time-frequency transformation and other technologies. The basis of the instantaneous frequency technology is the Hilbert transform, and then the instantaneous frequency is obtained by taking the derivative of the instantaneous phase with respect to time, but the negative frequency without physical meaning is obtained by directly calculating according to the analytical signal formula, which is not conducive to frequency evaluation; the instantaneous Q value can evaluate the attenuation characteristics of the stratum, but the same is based on the calculation of the instantaneous frequency, and cannot analyze the proportion of different frequency components at the same time; the frequency division attribute is to filter out the signals of different frequency bands of the seismic signal through filtering technology, and analyze the differences between different frequency bands, so as to evaluate the tuning frequency and gas-bearing property of the stratum, but since it is still based on the analysis of the amplitude and waveform characteristics, no factor directly related to the gas-bearing property is found; the wavelet transform, S transform and other two-dimensional time-frequency analysis technologies can decompose the non-stationary time domain signal and analyze the signal composition from the time domain and the frequency domain, and are the best tools for non-stationary signal analysis and the commonly used basic technologies for predicting the gas-bearing property of the earthquake, but since the energy of the two-dimensional time-frequency spectrum formed by the time-frequency decomposition is scattered, the energy group is large, the time and frequency resolution cannot be considered, and it is difficult to accurately estimate the frequency composition at a certain time point, and usually the gas-bearing evaluation effect is not good. Therefore, under the existing technical means, the gas-bearing prediction of the conventional post-stack seismic data is often not good, and can only be used as a reference. SUMMARY

[0004] The present application is aimed at the estimation problem of the post-stack seismic signal in the gas-bearing prediction, and proposes a gas-bearing prediction method and system based on the Hilbert-Huang transform time-frequency spectrum, which clearly splits the frequency component by using the Hilbert-Huang transform, and constructs a gas-bearing factor representing the lithological difference and the gas-bearing difference of the reservoir, thereby improving the prediction persuasiveness of the reservoir gas-bearing property.

[0005] In order to achieve the above-mentioned application purposes, the technical scheme of the present application is as follows:

[0006] A gas-bearing prediction method based on the Hilbert-Huang transform time-frequency spectrum, comprising the following steps:

[0007] Step one, for the seismic data to be analyzed section, for each seismic trace to be analyzed section of the empirical mode decomposition (EMD), to obtain linear stable signal (IMF);

[0008] Step two, for each IMF component to do Hilbert transform, and to obtain the instantaneous frequency and envelope;

[0009] Step three, all IMF component of the instantaneous frequency and envelope, in the need to analyze the frequency range synthesis, get Hilbert-Huang transform two-dimensional time-frequency spectrum P(t, f);

[0010] Step four, extract the instantaneous spectrum P(t i , f) of a seismic trace sample point, record the highest and the second highest frequency and the corresponding amplitude; the trapezoidal spectrum area formed by the two spectrum lines is calculated, and then divided by the highest amplitude frequency to construct the gas factor T;

[0011] Step five, calculate the gas factor of all sample points in the section to be analyzed, and get the evaluation data volume of the reservoir gas content.

[0012] Further, in step four, the frequency range to be analyzed is the main frequency range of the seismic data, which is 1-100Hz.

[0013] Further, the method for obtaining the Hilbert-Huang transform two-dimensional time-frequency spectrum comprises: constructing a time-frequency spectrum matrix with frequency as the horizontal axis, time as the vertical axis, and amplitude as the Z-axis; positioning each IMF component according to its time and instantaneous frequency in the time-frequency spectrum matrix, and then summing the envelope amplitudes of all IMF components at the corresponding positions.

[0014] Further, for each seismic trace sample point, the extracted instantaneous spectrum P(t i , f) is the spectrum of the Hilbert-Huang transform spectrum P(t, f) at the t i moment corresponding to the seismic sample point.

[0015] Further, in step four, according to the amplitude, the highest amplitude F1 and the corresponding frequency f1, the second highest amplitude E2 and the corresponding frequency f2 are found; the two frequencies form a trapezoidal on the spectrum, and the area is calculated by S=(E1+E2)*(f1-f2) / 2.

[0016] Further, the seismic data in step one is the pre-stack time migration seismic waveform data of the working area.

[0017] Further, the empirical mode decomposition in step one is an important part of the Hilbert-Huang transform proposed by Dr. Norden E. Huang. The specific steps of empirical mode decomposition for each seismic trace to be analyzed include:

[0018] Find all extreme points of signal x(t);

[0019] Fit the envelope lines s1(t) and s2(t) of the upper and lower extreme points with cubic spline curve, and calculate the average value m(t) = [s1(t)-s2(t)] / 2, and subtract it from x(t): h(t) = x(t)-m(t);

[0020] Determine whether h(t) is IMF; if the number of extreme points of h(t) is not equal to the number of zero-crossing points, then it is not IMF;

[0021] If h(t) is not IMF, replace x(t) with h(t), and repeat the above steps until h(t) meets the criteria, then h(t) is the IMF to be extracted;

[0022] Each time an IMF is obtained, subtract it from the original signal, repeat the above steps until the last remaining part r n is a monotonic sequence or a constant sequence;

[0023] In this way, the original signal x(t) is decomposed into a series of IMF and the linear superposition of the remaining part, where C i represents the i-th IMF component, and r n represents the remaining signal after EMD decomposition:

[0024]

[0025] Further, the steps of Hilbert transform for each IMF component include:

[0026] For signal x(t), its Hilbert transform can be represented in time domain as:

[0027]

[0028] where represents the Hilbert transform of signal x(t), i.e. the convolution of x(t) and signal .

[0029] According to the formula, the envelope, instantaneous phase, and instantaneous frequency can be calculated:

[0030]

[0031]

[0032]

[0033] Wherein, a(t) is signal envelope, namely instantaneous amplitude;φ(t) is instantaneous phase, and derivation of φ(t) can obtain instantaneous frequency f(t) of signal x(t).

[0034] The application further provides a gas-bearing property prediction system based on Hilbert-Huang transform time-frequency spectrum, comprising:

[0035] The data processing module one is used for carrying out empirical mode decomposition on the to-be-analyzed segment of each seismic trace, and obtaining linear steady signal IMF;

[0036] The data processing module two is used for carrying out Hilbert transform on each IMF component, and obtaining instantaneous frequency and envelope;

[0037] The data processing module three is used for synthesizing the instantaneous frequency and envelope of all IMF components in the frequency range to be analyzed, and obtaining Hilbert-Huang transform two-dimensional time-frequency spectrum P(t, f);

[0038] The data processing module four is used for extracting the instantaneous spectrum P(t i , f) of a sample point of each seismic trace, recording the frequency with the highest amplitude and the second highest amplitude and the corresponding amplitude;the spectrum area of the trapezoid formed by the two spectrum lines is obtained, and the spectrum area is divided by the frequency with the highest amplitude, and a gas-bearing factor T is constructed;

[0039] The data processing module five is used for calculating the gas-bearing factor of all sample points of the to-be-analyzed segment, and obtaining the evaluation of the gas-bearing property of the reservoir.

[0040] In conclusion, the application has the following advantages:

[0041] 1. In the application, the IMF component is processed by using Hilbert-Huang transform, so that the main frequency component of the signal is decomposed, the frequency composition of the signal at a certain moment is clearly indicated, and a good basis is provided for frequency analysis related to gas-bearing property, and the frequency characteristic analysis related to gas-bearing property is beneficial to be analyzed;

[0042] 2. The application proposes and constructs a gas-bearing factor for representing the lithology difference and gas-bearing property difference of the reservoir, the gas-bearing factor has strong correlation with the actual well test production, and therefore the prediction of the gas-bearing property of the reservoir is more convincing;The gas-bearing factor obtained by the gas-bearing factor is composed of the spectrum area of the molecule and the frequency with the highest amplitude of the denominator, can represent the lithology difference intensity of the reservoir segment and the high-frequency attenuation condition at the same time, has high conformity with the actual gas production, has good reliability in gas-bearing property measurement, and is a reliable means for evaluating the gas-bearing property in low-speed sandstone reservoirs. BRIEF DESCRIPTION OF DRAWINGS

[0043] Figure 1 The Hilbert-Huang transform of the signal t i The spectral area diagram of the time instant;

[0044] Figure 2 The flow chart of the gas content prediction method of the present application;

[0045] Figure 3 The relationship between the Hilbert transform time-frequency spectrum of the seismic trace near the well and the test production;

[0046] Figure 4 The comparison of the gas content factor profile of the high-yield well and the low-yield well;

[0047] Figure 5 The gas content factor profile of the high-yield well in the same channel. DETAILED DESCRIPTION

[0048] In order to more clearly illustrate the present application, the present application will be further described below in conjunction with the preferred embodiments and the accompanying drawings. Those skilled in the art should understand that the specific description below is illustrative rather than limiting, and should not limit the protection scope of the present application.

[0049] The technical solution of the instantaneous frequency and the instantaneous Q value for evaluating the reservoir gas content calculates the Hilbert transform of the signal, estimates the instantaneous frequency, analyzes the attenuation of the formation, and thus represents the characteristics of the fluid in the formation in the seismic signal. However, this kind of method only uses the Hilbert transform, and for the actual seismic signal which is complex and non-stationary, the frequency component in the signal cannot be analyzed, and the single instantaneous frequency cannot accurately represent the real composition of the signal at a time instant. In order to solve this problem, the present application provides a gas content prediction method based on the Hilbert-Huang transform time-frequency spectrum, which comprises:

[0050] Step one, input the seismic data body, determine the data range and the frequency range of the segment to be analyzed;

[0051] Step two, perform empirical mode decomposition on the segment to be analyzed of each seismic trace, and obtain the linear and stationary signal IMF;

[0052] Step three, perform Hilbert transform on each IMF component, and obtain the instantaneous frequency and the envelope;

[0053] Step four, collect the instantaneous frequency and the envelope of all IMF components, synthesize in the frequency range to be analyzed, and obtain the two-dimensional Hilbert-Huang transform time-frequency spectrum (t, f);

[0054] Step five, for each seismic trace sample point, extract the instantaneous spectrum at the sample point, record the highest and second highest frequency f1, f2 and the corresponding amplitude E1, E2, the area S of the trapezoid formed by the two spectrum lines is called the spectral area, the area is calculated and divided by the highest amplitude frequency, and the gas content factor T = S / f1 is constructed;

[0055] Step six, calculate the gas content factor of all sample points in the analyzed section to obtain the evaluation data volume of the reservoir gas content.

[0056] The overall idea of the method is based on two-dimensional time-frequency decomposition method, and the Hilbert-Huang transform is selected on the time-frequency decomposition method, which has the characteristics of clear instantaneous frequency component splitting of non-stationary signal, which is convenient for calculating high frequency attenuation, analyzing instantaneous frequency component, evaluating reservoir lithology difference and high frequency attenuation, and evaluating reservoir gas content, which avoids the problem that the energy dispersion phenomenon of wavelet transform and S transform leads to unclear frequency composition analysis, and is beneficial to the analysis of frequency characteristics related to gas content.

[0057] At the same time, by implementing the technical scheme, the evaluation data volume of the reservoir gas content is finally obtained, which can directly show the gas-bearing favorable positions in the reservoir, and according to these favorable positions, data basis can be provided for geological modeling or well deployment argument. The method can be used for time domain and depth domain seismic data, and is suitable for low-velocity and low-impedance channel sand gas reservoirs.

[0058] Embodiment 1

[0059] The following provides a specific application example to illustrate the gas content prediction method based on Hilbert-Huang transform time-frequency spectrum. The embodiment is the application of the present application in the Shaximiao Formation gas reservoir of a structure in a basin. The gas reservoir is a structural-lithologic gas reservoir, and the channel deposition develops several main channels in the work area, and the channel sand body is the main reservoir with low velocity and low wave impedance characteristics. The morphology of the river can be described by conventional seismic amplitude attributes, but the conventional seismic frequency attribute related to the fluid cannot establish a good relationship with the actual gas production of the well, and the channel gas content prediction effect is poor. Therefore, under the condition of only post-stack data, the gas content prediction method based on the seismic time-frequency spectrum area method is carried out in the work area, and the gas content in the main channel is described.

[0060] Referring to the description of the accompanying drawings Figure 2 , the implementation process is as follows:

[0061] Step one, input the seismic data volume, determine the data range and frequency range of the analyzed section according to the purpose layer tracking; the input seismic data volume in this step is the pre-stack time migration seismic waveform data volume of the work area, and the seismic waveform data volume is referred to as seismic data hereinafter.

[0062] Step two, the experience model decomposition is carried out to the segment to be analyzed of each seismic trace, and linear steady signal IMF is obtained, and the specific steps are as follows:

[0063] 1) find all extreme points of signal x(t);

[0064] 2) the envelope lines s1(t) and s2(t) of upper and lower extreme points are fitted by using 3-order spline curve, and the average value m(t) of the upper and lower envelope lines is calculated: m(t) = [s1(t)-s2(t)] / 2, and it is subtracted from x(t): h(t) = x(t)-m(t);

[0065] 3) whether h(t) is IMF is judged. If the number of extreme points of h(t) is not equal to the number of zero-crossing points, h(t) is not IMF;

[0066] 4) if h(t) is not IMF, h(t) is replaced by x(t), and the above steps are repeated until h(t) meets the criterion, then h(t) is the IMF to be extracted;

[0067] 5) each obtained IMF is subtracted from the original signal, and the above steps are repeated until the last remaining part r n of the signal is a monotonic sequence or a constant sequence.

[0068] 6) in this way, the original signal x(t) is decomposed into a series of IMFs and the linear superposition of the remaining part, wherein C i represents the i-th IMF component, and r n represents the remaining signal after EMD decomposition:

[0069]

[0070] Step three, Hilbert transform is carried out to each IMF component, and the instantaneous frequency and envelope are calculated;

[0071] Hilbert transform is a common means in signal processing. For a signal x(t), the Hilbert transform in time domain can be expressed as:

[0072]

[0073] Wherein represents the Hilbert transform of the signal x(t), that is, the convolution of x(t) and the signal .

[0074] According to the formula, the envelope, instantaneous phase and instantaneous frequency can be calculated:

[0075]

[0076]

[0077]

[0078] Wherein, a(t) is signal envelope, namely instantaneous amplitude; φ(t) is instantaneous phase, and derivation of φ(t) can obtain instantaneous frequency f(t) of the signal x(t).

[0079] Step four, the instantaneous frequency and envelope of all IMF components are summarized, and are synthesized in the frequency range to be analyzed to obtain the Hilbert-Huang transform two-dimensional time-frequency spectrum (t, f); the spectrum represents the frequency component composition at each time point; in this step, the frequency range to be analyzed is the main frequency range of the seismic data, and is generally taken as 1-100 Hz;

[0080] The method for calculating the Hilbert-Huang transform two-dimensional time-frequency spectrum P(t, f) is that: a time-frequency spectrum matrix with the frequency as the horizontal axis, the time as the vertical axis and the amplitude as the Z axis is constructed; each IMF component obtained in steps two to three is positioned in the time-frequency spectrum matrix according to the time and the instantaneous frequency thereof, and then the envelope amplitudes of all the IMF components are summed at the corresponding positions.

[0081] Step five, for each seismic trace sample, the instantaneous spectrum at the sample is extracted, the highest and the second highest frequencies f1, f2 and the corresponding amplitudes E1, E2 are recorded, the area S of the trapezoid formed by the two spectrum lines is called the spectrum area, the area is calculated and divided by the frequency of the highest amplitude to construct the gas-bearing factor: T=S / f1.

[0082] As shown in Figure 1 , for each seismic sample, the instantaneous spectrum P(t i , f) extracted is the spectrum of the Hilbert-Huang transform spectrum P(t, f) at the t i moment corresponding to the seismic sample. Since the number of IMF is small, the peak value of the spectrum P(t i , f) is sparse and easy to determine. According to the amplitude, the highest amplitude E1 and the corresponding frequency f1, the second highest amplitude E2 and the corresponding frequency f2 are found. The two frequencies can form a trapezoid on the spectrum, and the area S=(E1+E2)*(f1-f2) / 2 is calculated, which is called the spectrum area.

[0083] As shown in Figure 3 , the construction basis of the gas-bearing factor T is that: the law is obtained from the comparison and analysis of the actual well production and the Hilbert-Huang transform two-dimensional time-frequency spectrum of the seismic trace beside the well, the spectrum area S of the reservoir part of the high-yield gas well is larger, and the f1 value is lower, on the contrary, the spectrum area of the low-yield well, dry well or non-reservoir is smaller, and the f1 value is higher. According to this law, the spectrum area S is taken as the numerator, and the f1 is taken as the denominator to construct the gas-bearing factor.

[0084] Step six, calculate the gas factor of all sample points in the analyzed section, and obtain the evaluation of the reservoir gas content.

[0085] The Figure 4 The time-domain gas factor profiles of two over-high and over-low production wells are shown. The line marked J2s2 in the profile represents the bottom boundary of the upper Shaximiao Formation of the Jurassic System, and the channel reservoir develops in the range of 15ms-25ms below the interface. From the calculated gas content data body, firstly, the high value area of the data body can better reflect the position of the reservoir and will not deviate from the main channel; secondly, the favorable gas-bearing area reflected by the gas factor corresponds to the yield of the well: profile 1 is the gas factor profile of over-high production well A, which has a higher gas factor (greater than 0.2) in the target layer position and a larger lateral extension; profile 2 is the gas factor profile of over-low production well C, which has a lower gas factor value (less than 0.1).

[0086] The Figure 5 The gas content difference of three wells in the same channel is shown. Well D has lower production than wells A and B, and its trajectory has a long section in the low gas factor section on the gas factor profile; well B has higher production than well A, one reason is that its horizontal section is longer, and the other reason is that the gas factor of the stratum it passes through is higher than that of well A, and the gas factor of well A is about 0.22, and the gas factor of well B is greater than 0.24 in some areas. It can be seen that the gas factor calculated by the technical solution is highly consistent with the actual gas production, and the reliability of the gas content is good.

[0087] Embodiment 2

[0088] The embodiment also provides a gas content prediction system based on the Hilbert-Huang transform time-frequency spectrum, comprising:

[0089] A data processing module one is configured to perform empirical mode decomposition on the to-be-analyzed section of each seismic trace for the to-be-analyzed section of the seismic data, and obtain a linear and stable signal IMF;

[0090] A data processing module two is configured to perform Hilbert transform on each IMF component, and obtain an instantaneous frequency and an envelope;

[0091] A data processing module three is configured to synthesize the instantaneous frequency and the envelope of all IMF components in a frequency range to be analyzed, and obtain a two-dimensional Hilbert-Huang transform time-frequency spectrum P(t, f);

[0092] A data processing module four is configured to extract an instantaneous spectrum P(t i , f) at a sample point of a seismic trace, record the highest and second-highest frequencies and the corresponding amplitudes, calculate the spectral area of the trapezoid formed by the two spectrum lines, and divide the highest amplitude frequency to construct a gas factor T;

[0093] The data processing module five is used for calculating the gas-bearing factors of all sample points in the analyzed section, and obtaining the evaluation of the gas-bearing property of the reservoir.

[0094] The function implementation of each data processing module of the prediction system in the embodiment corresponds to each method step in the above embodiment, and will not be described here.

[0095] The above is only the preferred embodiment of the present application, and does not limit the present application in any form. Any simple modification and equivalent change of the above embodiment according to the technical essence of the present application falls within the protection scope of the present application.

Claims

1. A method for predicting gas content based on the time-spectrum of Hilbert-Huang transform, characterized in that, Includes the following steps: Step 1: For the seismic data segment to be analyzed, perform empirical mode decomposition on the segment to be analyzed for each seismic trace and obtain the linear steady-state signal IMF; Step 2: Perform a Hilbert transform on each IMF component and calculate the instantaneous frequency and envelope; Step 3: Combine the instantaneous frequencies and envelopes of all IMF components within the frequency range to be analyzed to obtain the Hilbert-Huang transform two-dimensional time spectrum P(t,f); Step 4: Extract the instantaneous spectrum P(t) at a certain seismic trace sampling point. i Record the frequencies of the highest and second highest amplitudes and their corresponding amplitudes; calculate the spectral area of ​​the trapezoid formed by the two spectral lines and divide it by the frequency of the highest amplitude to construct the gas-containing factor T. Step 5: After calculating the gas-bearing factor of all samples in the section to be analyzed, the gas-bearing capacity of the reservoir is evaluated. The method for obtaining the two-dimensional time spectrum of the Hilbert-Huang transform includes: constructing a time spectrum matrix with frequency on the horizontal axis, time on the vertical axis, and amplitude on the Z-axis; locating each IMF component in the time spectrum matrix according to its time and instantaneous frequency; and then summing the envelope amplitudes of all IMF components at the corresponding positions. The specific steps for empirical model decomposition of each seismic trace segment to be analyzed include: Find all the extreme points of the signal x(t); Use cubic spline curves to fit the envelopes s1(t) and s2(t) of the upper and lower extreme points, and calculate the average value of the upper and lower envelopes m(t) = [s1(t) - s2(t)] / 2. Subtract it from x(t): h(t) = x(t) - m(t); Determine if h(t) is an IMF; if the number of extreme points of h(t) is not equal to the number of zero crossings, then it is not an IMF. If h(t) is not an IMF, then replace x(t) with h(t) and repeat the above steps until h(t) satisfies the criterion. Then h(t) is the IMF that needs to be extracted. For each order IMF obtained, it is subtracted from the original signal. This process is repeated until the last remaining part r of the signal is obtained. n It is a monotonic sequence or a constant sequence; Thus, the original signal x(t) is decomposed into a series of IMFs and a linear superposition of the remaining parts, where C i Let r represent the i-th IMF component. n This represents the remaining signal after EMD decomposition: ； The steps for performing a Hilbert transform on each IMF component include: For a signal x(t), its Hilbert transform in the time domain can be expressed as: ； in Indicates signal The Hilbert transform, i.e. With signal Convolution; The envelope, instantaneous phase, and instantaneous frequency can be calculated using the formula: ； ； ； in, This refers to the signal envelope, i.e., the instantaneous amplitude. For instantaneous phase, The signal can be obtained by taking the derivative. instantaneous frequency .

2. The method according to claim 1, characterized in that, In step three, the frequency range to be analyzed is 1-100Hz.

3. The method according to claim 1, characterized in that, For each seismic trace point, the extracted instantaneous spectrum P(t) i f) is the Hilbert-Huang transform spectrum P(t,f) at the seismic sample point corresponding to t. i The spectrum at any given moment.

4. The method according to claim 1, characterized in that, In step four, based on the amplitude, find the highest amplitude E1 and its corresponding frequency f1, and the second highest amplitude E2 and its corresponding frequency f2; these two frequencies form a trapezoid on the spectrum, and the area is calculated as S=(E1+E2)*(f1-f2) / 2.

5. The method according to claim 1, characterized in that, The seismic data in step one is the pre-stack time-migrated seismic waveform data of the work area.

6. A gas-bearing prediction system based on the Hilbert-Huang transform time-spectrum, characterized in that, include: Data processing module one is used to perform empirical mode decomposition on the seismic data segment to be analyzed for each seismic trace and obtain the linear steady-state signal IMF. Data processing module two is used to perform Hilbert transform on each IMF component and obtain the instantaneous frequency and envelope; Data processing module three is used to synthesize the instantaneous frequencies and envelopes of all IMF components within the frequency range to be analyzed, and obtain the Hilbert-Huang transform two-dimensional time spectrum P(t,f); Data processing module four is used to extract the instantaneous spectrum P(t) at a certain seismic trace sampling point. i Record the frequencies of the highest and second highest amplitudes and their corresponding amplitudes; calculate the spectral area of ​​the trapezoid formed by the two spectral lines and divide it by the frequency of the highest amplitude to construct the gas-containing factor T. Data processing module five is used to calculate the gas-bearing factor of all samples in the section to be analyzed, and to obtain an assessment of the gas-bearing properties of the reservoir. The method for obtaining the two-dimensional time spectrum of the Hilbert-Huang transform by the data processing module three includes: constructing a time spectrum matrix with frequency on the horizontal axis, time on the vertical axis, and amplitude on the Z-axis; locating each IMF component in the time spectrum matrix according to its time and instantaneous frequency; and summing the envelope amplitudes of all IMF components at the corresponding positions. The specific steps of the data processing module in performing empirical mode decomposition on each seismic trace segment to be analyzed include: Find all the extreme points of the signal x(t); Use cubic spline curves to fit the envelopes s1(t) and s2(t) of the upper and lower extreme points, and calculate the average value of the upper and lower envelopes m(t) = [s1(t) - s2(t)] / 2. Subtract it from x(t): h(t) = x(t) - m(t); Determine if h(t) is an IMF; if the number of extreme points of h(t) is not equal to the number of zero crossings, then it is not an IMF. If h(t) is not an IMF, then replace x(t) with h(t) and repeat the above steps until h(t) satisfies the criterion. Then h(t) is the IMF that needs to be extracted. For each order IMF obtained, it is subtracted from the original signal. This process is repeated until the last remaining part r of the signal is obtained. n It is a monotonic sequence or a constant sequence; Thus, the original signal x(t) is decomposed into a series of IMFs and a linear superposition of the remaining parts, where C i Let r represent the i-th IMF component. n This represents the remaining signal after EMD decomposition: ； The steps of the data processing module two in performing a Hilbert transform on each IMF component include: For a signal x(t), its Hilbert transform in the time domain can be expressed as: ； in Indicates signal The Hilbert transform, i.e. With signal Convolution; The envelope, instantaneous phase, and instantaneous frequency can be calculated using the formula: ； ； ； in, This refers to the signal envelope, i.e., the instantaneous amplitude. For instantaneous phase, The signal can be obtained by taking the derivative. instantaneous frequency .

Images