CN102323615B - Method for reservoir prediction and fluid identification with earthquake data and device - Google Patents

Abstract

The invention provides a method for reservoir prediction and fluid identification with earthquake data and a device. The method comprises the steps that: the earthquake data is read in; fractional Gabor conversion is carried out to the earthquake data, so that the time frequency distribution of signals is obtained; fluid factors are extracted from the obtained time frequency distribution; and reservoir prediction and fluid identification are carried out through the fluid factors. The processing and the application of actual data show that: the introduction of fractional Gabor conversion into the field of earthquake signal processing has the advantages of good local time frequency energy collection, multi-scale decomposition capability under different orders and the like, and can provide a spectrum profile with high resolution and fine features of the difference of distributed fluid. The fluid attribute factors which are extracted based on the invention can objectively reflect the reservoir information, so that the prediction precision is improved.

Description

Utilize geological data to carry out the method and apparatus of reservoir prediction and fluid identification

Technical field

The present invention relates to the method and apparatus of a kind of reservoir prediction and fluid identification, specifically, relate to and a kind ofly extract fluid factor by the fractional order Gabor conversion that geological data is carried out under optimal order, and utilize described fluid factor to carry out the method and apparatus of reservoir prediction and fluid identification.

Background technology

Nowadays, utilize seismic data directly to carry out reservoir prediction and fluid identification, become one of the mainstream development direction of Earthquake Reservoir and the study hotspot in exploration geophysics field.It is directly the acoustic feature (amplitude) of rock due to the seismological observation data, be not the direct reflection of reservoir, although can pass through various inversion methods, obtain various Reservoir Parameters, but this is the Nonlinear Mapping relation of a height, and inverting has multi-solution.In addition, oil in reservoir, gas, water or other fluid nature difference are extremely faint.Therefore, utilizing seismic data directly to carry out reservoir prediction and fluid identification is a difficult point.

At present, relevant scholar has carried out a large amount of research work to the problem of utilizing seismic data directly to carry out reservoir prediction and fluid identification both at home and abroad.Relatively effective method mainly comprises the poststack attributive analysis, and the prestack inversion technology such as conventional AVO analysis, elastic impedance, elastic parameter are based on the parametric inversion of attenuation by absorption theory and various spectral factorization methods based on time frequency analysis etc.

Time-Frequency Analysis Method is the important method that non-stationary signal is processed, and it utilizes the Copula of time and frequency to represent non-stationary signal, and it is analyzed and processes.The conjoint analysis method of time-frequency domain has become a powerful tool of Analysis of Complex signal.Seismic signal is typical non-stationary signal, and after seismic data decomposed by time-frequency spectrum, the result that time-frequency spectrum is decomposed was subjected to character in reservoir thickness, deposition characteristics, factor of porosity, hole and the factor controllings such as natural attenuation of seismic signal.Therefore, can pass through the time-frequency spectrum decomposition based on the seismic data of time frequency analysis, utilize the statistical property of different frequency composition in seismic signal to detect the important information relevant with reservoir structure and fluid properties, suitable Time-Frequency Analysis Method is carry out seismic signal Wigner distribution, detection reservoir structure important information basic and crucial.

The below analyzes several Time-Frequency Analysis Method commonly used simply.

1, linear time-frequency distribution

(1) Fourier (Fourier) conversion in short-term

Short time discrete Fourier transform first adds time window with time signal, obtains like this signal frequency in window, then time window is slided and does Fourier transform, just obtains the time varying spectrum of signal.Therefore, short time discrete Fourier transform is that a segment signal with time window represents its characteristic sometime.But due to the constraint that is subjected to uncertain theorem, its time resolution and frequency resolution can not improve simultaneously: window is wider, and temporal resolution is poorer, for improving temporal resolution, window is narrowed down, but frequency resolution decreases.

(2) wavelet transformation

Wavelet transformation is to characterize this signal with signal at the space projection that cluster basis function (small echo) forms, and this cluster function is to consist of with translation by stretching of the different scale of basic mother and sons' wave function.The method is yardstick form changeable Time-Frequency Localization analytical approach all of the time window of fixing but its window function of a kind of window size (being the window total area) and frequency window.Compare with short time discrete Fourier transform, it preferably resolves the contradiction of temporal resolution and frequency resolution, and utilized cleverly the resolution of non-uniform Distribution, with high frequency resolution and low temporal resolution, adopt low frequency resolution and high temporal resolution at high band in low-frequency range.In case but wavelet basis function selected after, its characteristic is just fixing, the wavelet function on each yardstick obtains by yardstick and translation transformation, due to the every decomposition of signal once, the length of approximation signal and details reduces half.There are differences between the approximation signal feature that obtains on different scale, adopt the wavelet function of deriving with a plurality of basis functions to be difficult to approach exactly the local signal feature on different scale during wavelet transformation, so can lose original temporal signatures during reconstruction signal.

2, Bilinear TFD

The Wigner-Vile of signal x (t) is distributed as:

${\mathrm{WVD}}_x(t,f)=\&Integral;x(t+\frac{\mathrm{\&tau;}}{2})x^{*}(t-\frac{\mathrm{\&tau;}}{2})e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000509703400012.png,HSA00000509703400032.png"\; state="\{\{state\}\}">j</figure-callout>}2\mathrm{\&pi;f\&tau;}}\mathrm{d\&tau;}$

It is basis and the core of numerous non-linear time-frequency analysis technologies that Wigner-Vile distributes.The Wigner-Vile of signal do not contain any window function in distributing, and avoided the problem that in the linear time-frequency distribution, temporal resolution and frequency resolution restrict mutually.Industry generally admits, the time frequency resolution that distributes without any a kind of time-frequency combination can surpass the Wigner-Vile distribution, because the time-frequency bandwidth product of this distribution has reached the lower bound that uncertainty principle provides.Yet there is serious intersection distracter in this analytical approach when many component signals are analyzed, limited to a great extent the application of the method.

3, parametrization time-frequency distributions

(1) Adaptive matching projection turriform decomposition algorithm

In linear time-frequency distribution (basis function decomposition) method, if basis function is similar to the principal ingredient of signal, only need the linear combination of minority basis function just can represent that more accurately signal, the result of decomposition will be sparse.Otherwise, if the textural difference of the shape of basis function and signal is very large, so just needing a large amount of or even infinite many basis functions, ability is accurate expression original signal enough, and the information of signal will be distributed on too many basis function, be unfavorable for effectively representing signal.The core concept of Parametric Time-frequency Analysis method be adopt that a warp is flexible, " atom " collection that time shift and warbled Gauss function form, seek the linear combination of best base function according to maximum matching pursuit principle on this collection, to reach the purpose of coming decomposed signal with the least possible basis function.This method to the time constant frequency component effect very good, but when signal and chirp Signal Matching, this coupling is equivalent to zeroth order approaches, cause exist in decomposable process many block and component between the mixing distortion.

(2) Chirplet conversion

In order to overcome the shortcoming of Adaptive matching projection turriform decomposition algorithm, the Chirplet conversion adopts that a warp is flexible, the Gauss function of time shift, frequency displacement and frequency ramps is as basis function, to replace the indeclinable Gauss basis function of frequency.This conversion shows remarkable performance when analyzing the chirp signal.Yet the chirplet conversion based on Adaptive matching projection turriform decomposition algorithm is in fact that any energy trace on time-frequency plane is approached with one group of line segment that tilts arbitrarily.Obviously, single order approach than zeroth order approach can be compacter expression chirp class signal.And the signal of occurring in nature that such was the case with is simple, in time during nonlinearities change, certainly will cause the increase of basis function number, thereby impact is to understanding and the annotation of decomposition result when the radio-frequency component of signal.

Therefore, a kind of more suitably Time-Frequency Analysis Method of needs is carried out the time-frequency spectrum decomposition to the seismic signal of complexity, thereby carries out more accurately reservoir prediction and fluid identification.

Summary of the invention

For overcoming the shortcoming of above-mentioned Time-Frequency Analysis Method, the present invention is incorporated into the seismic data processing field with score field Time-Frequency Analysis Method (FrGT).The invention provides and a kind ofly extract fluid factor by geological data being carried out optimum fractional order Gabor conversion, and utilize described fluid factor to carry out the method and apparatus of reservoir prediction and fluid identification.Optimum fractional order Gabor conversion can provide for us the fine feature of high-resolution frequency spectrum section and spread fluid difference.The fluid that extracts on its basis belongs to the factor, can objectively react reservoir information, improves precision of prediction.

According to an aspect of the present invention, provide a kind of method of utilizing geological data to carry out reservoir prediction and fluid identification, can comprise: read in geological data; Geological data is carried out fractional order Gabor conversion with the time-frequency distributions of picked up signal; Extract fluid factor on the time-frequency distributions that obtains; Utilize fluid factor to carry out reservoir prediction and fluid identification.

The step of wherein, carrying out fractional order Gabor conversion can comprise: set order p; The window function h (n) of the fractional order Gabor conversion of setting signal; Ask for the quadrature analysis window function γ (n) corresponding with described window function h (n); Each road earthquake data x (n) is carried out windowing process carrying out fractional order Gabor conversion, thus the time-frequency distributions of picked up signal.

Wherein, the step of carrying out fractional order Gabor conversion also can comprise: after setting order p, determine the optimal order of fractional order Gabor conversion.

Wherein, determine that the step of the optimal order of fractional order Gabor conversion can comprise: (1) with order p as the initialization order, and definite iteration step length dp, smallest degree p _minWith maximum order p _max(2) the broad sense time-frequency bandwidth product under the corresponding order p of calculating signal; (3) order p is added the resulting result of iteration step length dp is as current order p; (4) as order p less than maximum order p _maxThe time, repeating step (2) and (3); Otherwise, carry out step (5); (5) the order p corresponding with minimum broad sense time-frequency bandwidth product is defined as the optimal order of geological data together.

Wherein, can adopt the method for averaging to calculate the optimal order of whole seismic section.

Wherein, the step of extracting fluid factor on the time-frequency distributions that obtains can comprise: according to the priori of seismic signal, extract fluid factor on the time-frequency distributions that obtains.

Wherein, fluid factor can comprise single-frequency section and energy attenuation gradient.

According to a further aspect in the invention, provide a kind of device that utilizes geological data to carry out reservoir prediction and fluid identification, can comprise: read in the unit, be used for reading in geological data; Converter unit is used for geological data is carried out fractional order Gabor conversion with the time-frequency distributions of picked up signal; Extraction unit is used for extracting fluid factor on the time-frequency distributions that obtains; Recognition unit utilizes fluid factor to carry out reservoir prediction and fluid identification.

Wherein, converter unit can be carried out: set order p; The window function h (n) of the fractional order Gabor conversion of setting signal; Ask for the quadrature analysis window function γ (n) corresponding with described window function h (n); Each road earthquake data x (n) is carried out windowing process carrying out fractional order Gabor conversion, thus the time-frequency distributions of picked up signal.

Wherein, converter unit also can comprise: the optimal order determining unit is used for determining the optimal order of fractional order Gabor conversion after setting order p.

Wherein, but optimal order determining unit execution in step is as follows: (1) with order p as the initialization order, and definite iteration step length dp, smallest degree p _minWith maximum order p _max(2) the broad sense time-frequency bandwidth product under the corresponding order p of calculating signal; (3) order p is added the resulting result of iteration step length dp is as current order p; (4) as order p less than maximum order p _maxThe time, repeating step (2) and (3); Otherwise, carry out step (5); (5) the order p corresponding with minimum broad sense time-frequency bandwidth product is defined as the optimal order of geological data together.

Wherein, can adopt the method for averaging to calculate the optimal order of whole seismic section.

Wherein, but extraction unit extracts fluid factor according to the priori of seismic signal on the time-frequency distributions that obtains.

Wherein, fluid factor can comprise single-frequency section and energy attenuation gradient.

Description of drawings

By the description of exemplary embodiment of the present invention being carried out below in conjunction with accompanying drawing, above and other purpose of the present invention and characteristics will become apparent, wherein:

Fig. 1 illustrates the non-rectangle coordinate grid corresponding with fractional order Gabor conversion (FrGT) according to exemplary embodiment of the present invention;

Fig. 2 illustrates the two dimensional surface coordinate according to the FrGT characterization signal of exemplary embodiment of the present invention;

Fig. 3 is that the geological data that utilizes that illustrates according to exemplary embodiment of the present invention carries out the process flow diagram of the method for reservoir prediction and fluid identification;

Fig. 4 is the process flow diagram that illustrates according to the FrGT of exemplary embodiment of the present invention;

Fig. 5 is the process flow diagram according to the optimal order of the FrGT of exemplary embodiment of the present invention;

Fig. 6 is the diagram that illustrates according to the FrGT Numerical Simulation Results of exemplary embodiment of the present invention;

Fig. 7 is the diagram that illustrates according to certain survey line CDP stacked section of carbonate roe beach, the Northeast, the river reservoir of exemplary embodiment of the present invention;

Fig. 8-Figure 11 is the diagram that the single-frequency section of certain survey line CDP stacked section under different score field frequencies in Fig. 7 is shown respectively;

Figure 12 is the diagram that the energy attenuation gradient profile of certain the survey line CDP stacked section in Fig. 7 is shown;

Figure 13 is that the geological data that utilizes that illustrates according to exemplary embodiment of the present invention carries out the block diagram of the device of reservoir prediction and fluid identification.

Embodiment

Below, describing embodiments of the invention in detail with reference to accompanying drawing, its example represents in the accompanying drawings, wherein, identical label represents identical parts all the time.

Fig. 1 illustrates the non-rectangle coordinate grid corresponding with fractional order Gabor conversion (FrGT) according to exemplary embodiment of the present invention.

The FrGT result that to be fractional order thought obtain in the popularization of traditional Gabor conversion.Improve the time-frequency distributions performance of traditional Gabor conversion by the kernel function of introducing FrFT (fraction Fourier conversion).With reference to Fig. 1, with respect to the rectangle time frequency grid characteristic of traditional Gabor conversion, FrGT obtains is non-rectangle time frequency grid more flexibly, and wherein, the t axle is time shaft, and the ω axle is frequency axis, and T and Ω are respectively time and frequency period.In this sense FrGT is the classical popularization of Gabor conversion in non-time frequency grid.When the frequency time-varying characteristics of basis function and signal did not mate, traditional Gabor conversion is the time-frequency local characteristics of reaction signal well, and fractional order Gaobr conversion can overcome this shortcoming.

Fig. 2 illustrates the two dimensional surface coordinate according to the FrGT characterization signal of exemplary embodiment of the present invention.

Discrete Fourier transform (DFT) (DFT) conversion is a kind of function f (t) to be rotated

By the representation of t principal axis transformation to the ω axle, namely function with the time shaft angle is being The ω axle on the expression be exactly the Fourier transform of this function.With reference to Fig. 2, what traditional time-frequency analysis technology characterized is the feature of signal in (t, f) two dimensional surface.According to exemplary embodiment of the present invention, FrGT passes through with time frequency plane anglec of rotation α (wherein,

P is the exponent number of fractional order conversion), with characterization in (t, u) two dimensional surface.According to the rotation relationship between score field sampling thheorem and score field and Fourier domain, the mapping relations between the frequency u of score field and Fourier domain frequency f are:

u＝f×sin(pπ/2)+t×cos(pπ/2)

As long as select the suitable anglec of rotation, FrGT just can make the time-frequency bandwidth product of signal reach minimum, thereby it has very high time-frequency aggregation, can be good at detecting the local feature information of signal.In addition, FrGT is a kind of linear time-frequency distribution, it can not produce the staggered problem of disturbing in similar Bilinear TFD, for traditional linear time-frequency distribution method, it has Duoed an adjustable factors (order p) than traditional linear time-frequency distribution method, when order arranges when reasonable, the time frequency analysis result than the better time-frequency aggregation of traditional time-frequency distributions method can be provided, better portray the local time-frequency characteristics (seeing Fig. 4) of signal.

The below carries out the method for reservoir prediction and fluid identification with reference to the geological data that utilizes that Fig. 3 describes in detail according to exemplary embodiment of the present invention.

Fig. 3 is that the geological data that utilizes that illustrates according to exemplary embodiment of the present invention carries out the process flow diagram of the method for reservoir prediction and fluid identification.

With reference to Fig. 3, in step 300, read in geological data, that is, and input seismic signals x (n).

In step 302, the seismic signals x (n) that reads in is carried out FrGT with the time-frequency distributions of picked up signal x (n).

In step 304, extract fluid factor on the time-frequency distributions that obtains.According to the priori of seismic signal x (n), extract fluid factor at time-frequency domain.Specifically, according to dominant frequency distribution range and the attenuation characteristics of seismic event when passing through hydrocarbon zone of seismic signal, extract single-frequency section and energy attenuation gradient at time-frequency domain.

In step 306, utilize the fluid factor that extracts to carry out reservoir prediction and fluid identification.

The below is with reference to the FrGT of Fig. 4 specific descriptions according to exemplary embodiment of the present invention.Fig. 4 is the process flow diagram that illustrates according to the FrGT of exemplary embodiment of the present invention.

With reference to Fig. 4, in step 400, set order p.

In step 402, the window function h (n) of the FrGT of setting signal x (n).

In step 404, ask for quadrature analysis window function γ (n) corresponding to window function h (n).

In exemplary embodiment of the present invention, with fractional order kernel function K _p(n, k) replaces traditional Gabor kernel function Obtain the FrGT expansion of signal x (n), as follows:

$x\left(n\right)=\underset{m=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{K-1}{\mathrm{\&Sigma;}}}{a}_{m,k,p}{\tilde{h}}_{m,k,p}\left(n\right)---\left(1\right)$

Wherein:

${\tilde{h}}_{m,k,p}\left(n\right)=\tilde{h}(n-\mathrm{mMb})K_p(n,k)---\left(2\right)$

$K_p(n,k)=\sqrt{\frac{\sin \mathrm{\&alpha;}-j\cos \mathrm{\&alpha;}}{N}}{e}^{\frac{j}{2}[{\mathrm{<figure-callout\; id="2"\; label="n"\; filenames="HSA00000509703400012.png,HSA00000509703400032.png"\; state="\{\{state\}\}">n</figure-callout>}}^2{\mathrm{\&Delta;<figure-callout\; id="2"\; label="t"\; filenames="HSA00000509703400012.png,HSA00000509703400032.png"\; state="\{\{state\}\}">t</figure-callout>}}^2+{\left(\mathrm{kKb}\right)}^2{\mathrm{\&Delta;u}}^2]\cot \mathrm{\&alpha;}}\&times;e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000509703400012.png,HSA00000509703400032.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;kKb}}{N}n}---\left(3\right)$

Wherein, n is the sampled point sequence number of signal, and Δ t is the discrete time interval, and Δ u is the discrete interval of score field,

For the cycle of window function h (n) expands, K _pBe the p rank kernel function under score field, α is rotation angle, order P be called the fractional order conversion exponent number (due to

Only appear on the parameter position of trigonometric function, so the definition take p as parameter is take 4 as the cycle, therefore only needs to investigate signal and get final product in one-period inner conversion result).K, M are respectively the sampling number of score field and time domain, and Kb, Mb are respectively the sampling interval of score field and time domain, and N is signal length, and satisfied relational expression is as follows:

M×Mb＝K×Kb＝N (4)

Corresponding FrGT is:

$a_{m,k,p}=\underset{n=-(N-1)/2}{\overset{(N-1)/2}{\mathrm{\&Sigma;}}}x\left(n\right){{\widetilde{\mathrm{\&gamma;}}}^{*}}_{m,k,p}\left(n\right)---\left(5\right)$

Wherein,

For

Corresponding quadrature analysis window function,

For Conjugation.

In formula (5) substitution formula (1), obtain:

$x\left(n\right)=\underset{i=-(N-1)/2}{\overset{(N-1)/2}{\mathrm{\&Sigma;}}}x\left(i\right)\underset{m=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{K-1}{\mathrm{\&Sigma;}}}{{\widetilde{\mathrm{\&gamma;}}}^{*}}_{m,k,p}\left(n\right){\tilde{h}}_{m,k,p}\left(i\right)---\left(6\right)$

By formula (6), have according to the completeness condition of fractional order basis function:

$\underset{m=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{k=0}{\overset{K-1}{\mathrm{\&Sigma;}}}{{\widetilde{\mathrm{\&gamma;}}}^{*}}_{m,k,p}\left(n\right){\tilde{h}}_{m,k,p}\left(i\right)=\mathrm{\&delta;}\left(i\right)---\left(7\right)$

Will

With

K _p(n, k) substitution formula (7) has according to Poisson summation formula:

$\frac{1}{\mathrm{Mb}\&times;\mathrm{Kb}}\underset{n=-(N-1)/2}{\overset{(N-1)/2}{\mathrm{\&Sigma;}}}\tilde{h}\left(n\right){\widetilde{\mathrm{\&gamma;}}}^{*}(n+\mathrm{mK})e^{-\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="HSA00000509703400012.png,HSA00000509703400032.png"\; state="\{\{state\}\}">j</figure-callout>}\frac{2\mathrm{\&pi;k}}{N}n}={\mathrm{\&delta;}}_m{\mathrm{\&delta;}}_k---\left(8\right)$

Wherein, i and n are respectively subscript corresponding in original window function h (n) and analysis window function gamma (n), δ (i) the expression i shock response of correspondence constantly, δ _mδ _kThe product of two shock responses of expression, only having the product when m is identical with k is not just 0.

Can try to achieve according to formula (8)

Corresponding quadrature analysis window function

In step 406, each road earthquake data x (n) is carried out windowing process carrying out FrGT, thus the time-frequency distributions of picked up signal.

Fig. 5 is the process flow diagram according to the optimal order of the FrGT of exemplary embodiment of the present invention.

According to exemplary embodiment of the present invention, after the step 400 of Fig. 4 is set order p, also can determine the optimal order p of FrGT _opt

With reference to Fig. 5, in step 500, the order p that will set in the step 400 of Fig. 4 is as the initialization order, and definite iteration step length dp, maximum order p _maxWith smallest degree p _min, wherein, make smallest degree p _minEqual initialization order p.

In step 502, the broad sense time-frequency bandwidth product GTBP{x under the corresponding order p of calculating signal x (n) _p(t) }, its formula is as follows:

$\mathrm{GTBP}\left\{x_p\left(t\right)\right\}=T_{x_p}\&times;B_{x_p}---\left(9\right)$

Formula (9) is at traditional TBP{x (t) }=T _x* B _xPromote a formula that obtains on formula, wherein,

With

Be respectively finite time width and the finite frequency width of signal x (n), its formula is as follows:

$T_{x_p}=\frac{{[\&Integral;{(t-{\mathrm{\&eta;}}_t)}^2{\left|x\left(t\right)\right|}^2\mathrm{dt}]}^{1/2}}{\left|\right|x\left|\right|},$ ${\mathrm{\&eta;}}_t=\frac{\&Integral;t{\left|x\left(t\right)\right|}^2\mathrm{dt}}{{\left|\right|x\left|\right|}^2},$

$B_{x_p}=\frac{{[\&Integral;{(u-{\mathrm{\&eta;}}_u)}^2{\left|X\left(u\right)\right|}^2\mathrm{du}]}^{1/2}}{\left|\right|x\left|\right|},$ ${\mathrm{\&eta;}}_u=\frac{\&Integral;f{\left|X\left(u\right)\right|}^2\mathrm{du}}{{\left|\right|x\left|\right|}^2},$

Wherein, ‖ x ‖ is the mould of signal, and u is the score field frequency, and u and f essence are same parameters, so η u can be written as again

X (u) is the result of (for example) p fraction Fourier conversion of signal x (t).

In step 504, order p is added the resulting result of iteration step length dp as current order p, that is, and p=p+dp.

In step 506, determine that whether current order p is more than or equal to maximum order p _maxIf determine that current order p is less than maximum order p _max, turn back to step 502, and execution in step 502.If determine that current order p is more than or equal to maximum order p _max, execution in step 508.

In step 508, determine the optimal order p of FrGT _optGTBP{x _p(t) } less, illustrate that time frequency resolution is higher; Otherwise, illustrate that time frequency resolution is lower.Therefore, according to the Compact support principle as can be known, signal the time meta-bandwidth product the size quality that can portray signal time-frequency aggregation, thereby can obtain by following formula the optimal order of one geological data:

$p_{\mathrm{opt}}=\arg \underset{0\&le;a<4}{\mathrm{<figure-callout\; id="0"\; label="min"\; filenames="HSA00000509703400032.png,HSA00000509703400043.png"\; state="\{\{state\}\}">min</figure-callout>}}\mathrm{GTBP}\left\{x_p\left(t\right)\right\}---\left(10\right)$

Wherein, the definition take p as parameter is take 4 as the cycle.

Because whole seismic section comprises multichannel seismic data, therefore, adopt the method for averaging to calculate the optimal order of whole seismic section.

Fig. 6 is the diagram that illustrates according to the FrGT Numerical Simulation Results of exemplary embodiment of the present invention.

To test signal x (t)=exp (j π (t ²+ 0.1 (t+0.5) ³)) exp (π t ²/ 9) carry out respectively traditional Gabor conversion and FrGT, and compare its performance, wherein, time t gets-8s～8s, and sampling rate is 64.With reference to Fig. 4, corresponding from left to right is respectively that result, the p of signal amplitude, traditional Gabor conversion is 0.7,1.0,1.4,1.175 FrGT result.

Calculate GTBP{x corresponding under different rank according to formula (9) _p(t) }, its result is as shown in table 1.Shown the time-frequency bandwidth product under different orders in expression 1.

Table 1

As can be seen from Table 1, when p=1.175, GTBP{x _p(t) } be 20.55, be minimal support set.Therefore, the optimal order p of the FrGT of this signal _optBe 1.175.As can be seen from Figure 4, when exponent number p=1, the result of FrGT is consistent with traditional Gabor transformation results; At optimal order p _opt=1.175 times, corresponding FrGT time-frequency aggregation obviously improves (seeing rightmost figure in Fig. 6), thereby has verified the remarkable result of the FrGT under the optimal order.

Below, be described on the time-frequency distributions that adopts FrGT to obtain with reference to Fig. 7-Figure 12 and extract fluid factor to carry out the method for reservoir prediction and fluid identification.Fig. 7 carries out the diagram of FrGT example before to geological data, Fig. 8-Figure 12 carries out the diagram of FrGT example afterwards to geological data.

Fig. 7 is the diagram that illustrates according to certain survey line CDP stacked section of carbonate roe beach, the Northeast, the river reservoir of exemplary embodiment of the present invention.With reference to Fig. 7, reservoir is positioned at the zone of demarcating with oval solid line.

Fig. 8-Figure 11 is the diagram that the single-frequency section of certain survey line CDP stacked section under different score field frequencies in Fig. 7 is shown respectively.Figure 12 is the diagram that the energy attenuation gradient profile of certain the survey line CDP stacked section in Fig. 7 is shown.

According to the priori of seismic signal, extract fluid factor to carry out reservoir prediction and fluid identification at time-frequency domain.Specifically, dominant frequency distribution range and the attenuation characteristics of seismic event when passing through hydrocarbon zone according to seismic signal, extract single-frequency section and energy attenuation gradient at time-frequency domain, thereby single-frequency section and energy attenuation gradient are applied to reservoir prediction and fluid identification.

With reference to Fig. 8-Figure 11, the dominant frequency cut-away section self-energy aggregation of seismic signal is very high.The single-frequency section of contrast different frequency, at HFS, the reservoir energy decay is obvious, and energy accumulating weakens gradually.

With reference to Figure 12, find out that obviously the energy attenuation gradient in reservoir zone is local peaking.This has effectively verified the front end attenuation characteristic of seismic event, for oil-gas recognition provides foundation.

Figure 13 is that the geological data that utilizes that illustrates according to exemplary embodiment of the present invention carries out the block diagram of the device 1300 of reservoir prediction and fluid identification.

With reference to Figure 13, according to exemplary embodiment of the present invention, a kind ofly utilize the device that geological data carries out reservoir prediction and fluid identification to comprise: to read in unit 1302, be used for reading in geological data; Converter unit 1304 is used for geological data is carried out FrGT with the time-frequency distributions of picked up signal; Extraction unit 1306 is used for extracting fluid factor on the time-frequency distributions that obtains; Recognition unit 1308 is used for utilizing fluid factor to carry out reservoir prediction and fluid identification.

In exemplary embodiment of the present invention, converter unit 1304 also can comprise: optimal order determining unit 1310 is used for determining the optimal order of fractional order Gabor conversion after setting order p.

In exemplary embodiment of the present invention, optimal order determining unit 1310 execution in step are as follows: (1) with order p as the initialization order, and definite iteration step length dp, smallest degree p _minWith maximum order p _max(2) the broad sense time-frequency bandwidth product under the corresponding order p of calculating signal; (3) order p is added the resulting result of iteration step length dp is as current order p; (4) as order p less than maximum order p _maxThe time, repeating step (2) and (3); Otherwise, carry out step (5); (5) the order p corresponding with minimum broad sense time-frequency bandwidth product is defined as the optimal order of geological data together.

In exemplary embodiment of the present invention, adopt the method for averaging to calculate the optimal order of whole seismic section.

In exemplary embodiment of the present invention, extraction unit 1306 extracts fluid factor according to the priori of seismic signal on the time-frequency distributions that obtains.Specifically, extraction unit 1306 extracts single-frequency section and energy attenuation gradient according to dominant frequency distribution range and the attenuation characteristics of seismic event when passing through hydrocarbon zone of seismic signal at time-frequency domain.

Although described the present invention with reference to exemplary embodiment, but it should be appreciated by those skilled in the art, in the situation that do not break away from the spirit or scope of the present invention, can make various variants and modifications to the present invention, scope of the present invention is limited by claim and equivalent thereof.

Claims (8)

1. a method of utilizing geological data to carry out reservoir prediction and fluid identification, comprise the steps:

Read in geological data;

Geological data is carried out fractional order Gabor conversion with the time-frequency distributions of picked up signal;

Extract fluid factor on the time-frequency distributions that obtains;

Utilize fluid factor to carry out reservoir prediction and fluid identification,

The step of wherein, carrying out fractional order Gabor conversion comprises:

Set order p;

Determine the optimal order of fractional order Gabor conversion;

The window function h (n) of the fractional order Gabor conversion of setting signal;

Ask for the quadrature analysis window function γ (n) corresponding with described window function h (n);

Each road earthquake data x (n) is carried out windowing process carrying out fractional order Gabor conversion, thus the time-frequency distributions of picked up signal,

Wherein, the step of determining the optimal order of fractional order Gabor conversion comprises:

(1) with order p as the initialization order, and definite iteration step length dp, smallest degree p _minWith maximum order p _max

(2) the broad sense time-frequency bandwidth product under the corresponding order p of calculating signal;

(3) order p is added the resulting result of iteration step length dp is as current order p;

(4) as order p less than maximum order p _maxThe time, repeating step (2) and (3); Otherwise, carry out step (5);

(5) the order p corresponding with minimum broad sense time-frequency bandwidth product is defined as the optimal order of geological data together.

2. the method for claim 1, wherein adopt the method for averaging to calculate the optimal order of whole seismic section.

3. the step of the method for claim 1, wherein extracting fluid factor on the time-frequency distributions that obtains comprises:

According to the priori of seismic signal, extract fluid factor on the time-frequency distributions that obtains.

4. the method for claim 1, wherein fluid factor comprises single-frequency section and energy attenuation gradient.

5. device that utilizes geological data to carry out reservoir prediction and fluid identification comprises:

Read in the unit, be used for reading in geological data;

Converter unit is used for geological data is carried out fractional order Gabor conversion with the time-frequency distributions of picked up signal;

Extraction unit is used for extracting fluid factor on the time-frequency distributions that obtains;

Recognition unit utilizes fluid factor to carry out reservoir prediction and fluid identification,

Wherein, converter unit is carried out: set order p; The window function h (n) of the fractional order Gabor conversion of setting signal; Ask for the quadrature analysis window function γ (n) corresponding with described window function h (n); Each road earthquake data x (n) is carried out windowing process carrying out fractional order Gabor conversion, thus the time-frequency distributions of picked up signal,

Wherein, converter unit also comprises: the optimal order determining unit, after setting order p, determine the optimal order of fractional order Gabor conversion,

Wherein, optimal order determining unit execution in step is as follows:

(1) with order p as the initialization order, and definite iteration step length dp, smallest degree p _minWith maximum order p _max

(2) the broad sense time-frequency bandwidth product under the corresponding order p of calculating signal;

(3) order p is added the resulting result of iteration step length dp is as current order p;

(4) as order p less than maximum order p _maxThe time, repeating step (2) and (3); Otherwise, carry out step (5);

(5) the order p corresponding with minimum broad sense time-frequency bandwidth product is defined as the optimal order of geological data together.

6. device as claimed in claim 5, wherein, adopt the method for averaging to calculate the optimal order of whole seismic section.

7. device as claimed in claim 5, wherein, extraction unit extracts fluid factor according to the priori of seismic signal on the time-frequency distributions that obtains.

8. device as claimed in claim 5, wherein, fluid factor comprises single-frequency section and energy attenuation gradient.

Images