CN102768365A - High resolution seismic wavelet extracting method based on high-order statistics and ARMA (autoregressive moving average) model - Google Patents

Abstract

The invention relates to a high resolution seismic wavelet extracting method based on high-order statistics and an ARMA (autoregressive moving average) model, which belongs to the field of seismic signal processing. The high resolution seismic wavelet extracting method provided by the invention is characterized in that: under the precondition of performing ARMA parameter simplified modeling on seismic wavelets, SVD (singular value decomposition) based on autocorrelation function is adopted to determine an order of AR part, an MA order determining method is provided for integrating information content criterion function in a high-order cumulant MA order determining method, and the MA order determining accuracy rate in a seismic wavelet ARMA model is improved; an SV-TLS (singular value decomposition - total least squares estimation) and a cumulant method are respectively adopted to estimate wavelet parameters; and under the precondition of ensuring wavelet precision, the model order is decreased as far as possible to improve the operation efficiency and to finally realize seismic wavelet extraction in high efficiency and high precision. Through the data simulation verification and the practical seismic data processing demonstration, the method provided by the invention is proved to effectively improve the estimated precision and extracting efficiency for the seismic wavelets and to have obvious effect even under short-time seismic data and strong noise pollution.

Description

High-resolution seismic exploration wavelet extraction method based on high-order statistic and arma modeling

Technical field

The invention belongs to the seismic data processing field.

Background technology

Nowadays exploration of oil and gas field develops with the direction of exploitation forward small-scale, thin reservoir, and is increasingly high to the accuracy requirement of seismic prospecting.Be the performance prediction and searching complex structure and lithologic deposit that adapts to hydrocarbon-bearing pool, require the section after earthquake is handled to have high s/n ratio, high resolving power and Hi-Fi " three height " characteristics.

Seismic prospecting comprises field acquisition, indoor processing and seismic data interpretation three big links.The resolution that improves seismologic record just must be excavated the ability of each link from gathering, handle and explaining that each link works hard, and makes all to reach high-quality on each link the seismic section that could obtain to have " three height " characteristics at last.And if known seismic wavelet is done deconvolution with this wavelet and seismologic record, just can try to achieve the reflection coefficient on stratum, thereby can get a high-resolution seismic section.Seismic wavelet estimation (wavelet extraction) is as the basis of seismic data deconvolution processing, question of seismic wave impedance inversion and earthquake forward model, and it extracts reliability and accuracy that precision directly influences latter earthquake Data Processing and seismic data interpretation.

In recent years; High-order statistic is demonstrating unique advantage aspect analysis that solves non-linear, non-gaussian sum non-minimum phase system and the processing gradually; For high-precision seismic wavelet extraction provides possibility, emerging nonlinear optimization method has further strengthened statistical seismic wavelet extraction The Application of Technology potentiality; The System Discrimination technology is also increasing to be incorporated in the seismic data processing.The present invention proposes to adopt higher order statistical theoretical; Seismic wavelet is carried out ARMA (Autoregressive Moving Average autoregressive moving average) parameter simplify modeling; Research and develop a kind of seismic wavelet arma modeling of practicality and decide the rank method, under the prerequisite that guarantees estimated accuracy, improve operation efficiency as far as possible.

Statistical wavelet extraction method is to utilize signal processing technology that limited seismic data recording is carried out statistical treatment; Make full use of the information such as amplitude, phase place, frequency of seismic data recording; Extract the estimation of the seismic wavelet that more is suitable for, but shortcoming is to carry out certain hypothesis (minimum phase, zero phase, maximum phase) to the distribution of used seismic data and underground reflection coefficient sequence.In fact a kind of often mixed phase wavelet of seismic wavelet is once this only simply based on not obtaining result accurately on the autocorrelative statistical wavelet extraction theoretical method.

Along with the higher order statistical theory further develops and perfect, the higher order statistical signal Processing has been penetrated into each application of signal Processing, and has emerged in large numbers a large amount of theories and applied research achievement.High-order statistic (higher-order spectrum) than second-order statistic (power spectrum) except the abundanter quantity of information (like complete phase information) that comprises signal; Also has characteristic to any Gaussian noise signal " blind "; Can in coloured Gaussian noise, extract non-Gaussian signal; This has solved the problem that a lot of other methods cann't be solved, like reconstruct non-minimum phase system, effective identification non-minimum phase system and NLS.To many characteristics of high-order statistic, geophysics worker recognizes gradually and uses the higher order statistical metering method can extract a non-minimum phase from seismic data actual seismic wavelet.

Summary of the invention

The objective of the invention is to propose a kind of high-resolution seismic exploration wavelet extraction method, simplify the foundation of arma modeling comprising the seismic wavelet parameter based on Higher Order Cumulants and arma modeling; The rank method of accurately deciding of seismic wavelet arma modeling.

One of characteristic of the present invention is that the arma modeling that adopts parameter to simplify first carries out modeling to seismic wavelet.

Seismic trace y (n) is assumed to be the stationary stochastic process of zero-mean usually, be formulated as a convolution and Gaussian noise v (n) and form, often write the seismic trace convolution model that we knew:

$y\left(n\right)=\underset{i=0}{\overset{q}{\mathrm{\Σ}}}w\left(i\right)r(n-i)+v\left(n\right)=w\left(n\right)*r\left(n\right)+v\left(n\right)---\left(1\right)$

Wherein w (n) is the mixed-phase seismic wavelet, and q is expressed as wavelet cause and effect partial-length; Reflection coefficient sequence r (n) be stably, zero-mean, independent identically distributed nongausian process; V (n) is the additivity coloured noise, and its gauss component is much larger than non-gauss component and separate with r (n).

Earthquake convolution model shown in the formula (1) is a typical FIR system (MA model), and digital q represents true wavelet length in the model, and model parameter value is represented the amplitude of wavelet.Wear semi-invariant matrix equation method that people such as Yongshou adopts the MA model and estimate that wavelet finds, when the length of wavelet greater than 7 the time, the maximum square error of model parameter increases sharply.But real seismic wavelet length is many more than 70ms, and under the SI of 1ms, MA model parameter to be estimated will be more than 70.With counting yield, MA model and impracticable in high precision seismic is handled.

According to the WOLD decomposition theorem, wear people such as Yongshou and at first arma modeling is applied among the seismic wavelet.Compare with the MA model, the exponent number p of arma modeling and the true length of q and seismic wavelet directly do not concern, can use accurate wavelet of less parametric description.We claim that these characteristics of arma modeling are " parameter is stingy ".

Suppose significant wave x (n) be stably, ergodic linear process, its difference equation is expressed as:

$\underset{i=0}{\overset{p}{\mathrm{\Σ}}}a\left(i\right)x(n-i)=\underset{j=0}{\overset{q}{\mathrm{\Σ}}}b\left(j\right)r(n-j)---\left(2\right)$

Significant wave is observed in additive noise, that is:

y(n)＝x(n)+v(n) (3)

V (n) is the coloured Gauss's observation noise that satisfies following formula:

$\underset{i=0}{\overset{p_v}{\mathrm{\Σ}}}c\left(i\right)v(n-i)=\underset{j=0}{\overset{q_v}{\mathrm{\Σ}}}d\left(j\right)e(n-j)---\left(4\right)$

Above-mentioned signal satisfies following hypothesis:

AS1) stratum reflection coefficient sequence r (n) is zero-mean, independent identically distributed nongausian process, and exist at least semi-invariant 0＜| γ _Kw|＜∞, k＞2;

AS2) neighbourhood noise v (n) is additivity Gauss's coloured noise of zero-mean, cause and effect, and independent with r (n) statistics, therefore also separate with x (n);

AS3) the seismic wavelet system is the cause and effect mixed-phase, and its transport function does not have pole zero cancellation, and definition is as follows:

$H\left(z\right)=\underset{i=0}{\overset{\∞}{\mathrm{\Σ}}}h\left(i\right)z^{-i}=\left[\underset{j=0}{\overset{q}{\mathrm{\Σ}}}b\left(j\right)z^{-j}\right]/\left[\underset{i=0}{\overset{p}{\mathrm{\Σ}}}a\left(i\right)z^{-i}\right]---\left(5\right)$

Wherein: a (p) ≠ 0, b (q) ≠ 0 supposes a (0)=b (0)=1.

Can find out that from above-mentioned model seismologic record can be considered the noisy output of stratum reflection coefficient sequence excitation seismic wavelet ARMA system.If will obtain seismic wavelet, only demand gets the arma modeling parameter and gets final product.The present invention will adopt based on the auto-correlation of sample function and fourth order cumulant method and confirm model order, simultaneously the estimation model parameter.

Two of characteristic of the present invention is in order under the prerequisite that guarantees the wavelet precision, to improve operation efficiency to greatest extent, has proposed the rank method of accurately deciding of arma modeling.

In the System Discrimination field,, estimate that accurately the exponent number of arma modeling and parameter are two key issues if obtain the output of given sequence through the ARMA process.A lot of ripe methods can the estimation model parameter, but all is under the known prerequisite of hypothesized model exponent number, so definite model order is the top priority of estimation model parameter.Owe to decide the distortion that wavelet can be caused in rank, the surge of operand can be caused in the overdetermination rank.Existing Causal model is decided the rank method has relevant function method, final predicated error method (FPE), red pond quantity of information criterion method (AIC), singular value decomposition method (SVD) etc., said method to decide the rank result still not ideal enough, remain further to be improved.The present invention will analyze the SVD method, and incorporate quantity of information criterion method and be applied to seismic wavelet estimation.

The present invention adopts and to confirm autoregression (AR) model order based on the autocorrelative SVD method of sample function, compares research with SVD method based on Higher Order Cumulants simultaneously; When the Higher Order Cumulants method is confirmed running mean (MA) model order, introduce quantity of information criterion method, effectively improved the MA model and decided the rank accuracy.

In the present invention, will write down the accuracy and the practicality of testing institute's application process through synthetic different earthquake.Suppose that the synthetic required reflection coefficient sequence of seismic trace is independent identically distributed stochastic process, its probability density is obeyed Bernoulli Jacob and is distributed.Seismic wavelet is the cause and effect mixed-phase, and its ARMA difference model is expressed as:

x(t)-1.6x(t-1)+0.8x(t-2)＝r(t)-0.8r(t-1)-1.2r(t-2) (6)

And the noise in the seismologic record is assumed to be additivity Gauss coloured noise, and the original wavelet of theogram is as shown in Figure 1.

The legal rank of svd (SVD) are through the structure Special matrix; Model is decided rank convert effective order of utilizing singular value decomposition method to find the solution Special matrix into; Theoretical proof SVD method shows very strong robust property in confirming the AR model order, in suitable noise circumstance, also have effect preferably.And the SVD method is not only applicable to pure AR model, can be used for confirming the AR part exponent number of arma modeling yet.Existing employing SVD confirms that the method for model order mainly contains the related function method that is applied to Gaussian process that Cadzow proposes, and the Higher Order Cumulants method that is applied to nongausian process of Giannakis and Mendel proposition.Can suppress Gaussian noise fully on the Higher Order Cumulants law theory, section has susceptibility but Special matrix is to semi-invariant, and operation time is long; The related function law technology is ripe, simple to operate, and operation time is short.

(exponent number of its MA part is q, has: b for p, q) process for the limited ARMA of exponent number _i=0 i＞q, must the Yule-Walker equation do,

$R_x\left(l\right)+\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{R}_x(l-i)=0,\∀l>q---\left(7\right)$

Following formula abbreviates the MYW equation as.

Can know by AR parameter identifiability theorem, if (p, q) the system function H (z) of model does not have zero limit cancellation to ARMA, and a _p≠ 0, then, the AR parameter a of this model ₁, a ₂... A _pCan be by p MYW equation

$\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i{R}_x(l-i)=-R_x\left(l\right)---\left(8\right)$

Only definite or identification.

Above-mentioned theorem is told us, at known AR model order p and autocorrelation function R _xUnder the prerequisite (τ), the AR model parameter can be revised the Yule-Walker equation and try to achieve through finding the solution p, and real data handle in AR model order p and autocorrelation function R _xAll be unknown (τ).We often adopt the ARMA (p with expansion exponent number _e, q _e) (p, q) process is simultaneously with autocorrelation function R to replace former ARMA _x(τ) use its estimated value

Replace, revise the Yule-Walker equation and still set up this moment, is expressed as with matrix form:

R _ea _e＝0 (9)

Wherein:

$R_e=\left[\begin{array}{cccc}{R}_x(q_e+1)& R_x\left(q_e\right)& ...& R_x(q_e+1-p_e)\\ R_x(q_e+2)& R_x(q_e+1)& ...& R_x(q_e+2-p_e)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ R_x(q_e+M)& R_x(q_e+M-1)& ...& R_x(q_e+M-p_e)\end{array}\right]---\left(10\right)$

$a_e={[1,a_1,...a_p,a_{p+1},...a_{p_e}]}^T---\left(11\right)$

Theoretical analysis proves, as M>=p _e>=p, q _e>=q, and q _e-p _eDuring>=q-p, rank (R _e)=p.

Constitute R _eThe autocorrelation function R of element samples function _x(τ) be second-order statistic, and second-order statistic has noise susceptibility in theory, possibly affect deciding the rank precision.Compare with second-order statistic, high-order statistic just has the characteristic of " blind " to Gaussian noise, can suppress Gaussian noise fully in theory, and Higher Order Cumulants can be accomplished all work that second-order statistic is done.Three rank, fourth order cumulant and frequency spectrum thereof are used morely in signal Processing in the Higher Order Cumulants, and fourth order cumulant is compared with three rank semi-invariants, can be used to analyze the signal of symmetrical distribution, and this paper will adopt the fourth order cumulant of sample function to make up the wavelet model.In view of the SVD based on autocorrelation function decides the rank method, we decide rank with the AR model and convert the order of confirming by the special section structural matrix of the fourth order cumulant of sample function into.

Structure is based on the expansion rank M of fourth order cumulant ₂(N ₂-N ₁+ 1) * M ₂Matrix C _e:

$C_e=\left[\begin{array}{ccc}{c}_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1,N_1)& ...& c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+M_2-1,N_1)\\ .& .& .\\ .& .& .\\ .& .& .\\ c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1,N_2)& ...& c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+M_2-1,N_2)\\ .& .& .\\ .& .& .\\ .& .& .\\ c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+{\mathrm{<figure-callout\; id="2"\; label="M"\; filenames="HDA0000058994520000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_2,N_1)& ...& c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+2M_2-1,N_1)\\ .& .& .\\ .& .& .\\ .& .& .\\ c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+{\mathrm{<figure-callout\; id="2"\; label="M"\; filenames="HDA0000058994520000031.png"\; state="\{\{state\}\}">M</figure-callout>}}_2,N_2)& ...& c_{4x}({\mathrm{<figure-callout\; id="1"\; label="M"\; filenames="HDA0000058994520000022.png,HDA0000058994520000041.png"\; state="\{\{state\}\}">M</figure-callout>}}_1+2M_2-1,N_2)\end{array}\right]---\left(12\right)$

Wherein, get M ₁>=q+1-p, M ₂>=p, N ₁≤q-p and N ₂>=q.As ARMA (p, q) the no pole zero cancellation of transfer function H (z), and a (p) ≠ 0, b (q) ≠ 0 o'clock, theoretical analysis proof rank (C _e)=p.

Real data has only limited sample, in formula (10), uses sampled correlation function

Replace ebsemble correlation function R _x(τ), in formula (12), use the sample fourth order cumulant

Replace overall fourth order cumulant c _4x(m, n).

Make up Special matrix R _eOr C _eAfter, try to achieve effective order of Special matrix through singular value decomposition method, can obtain the exponent number p of AR model.Matrix is carried out svd, is meant, then have a m * m unitary matrix U and m * n unitary matrix V, matrix A is decomposed into if matrix A is the complex matrix of m * n:

A＝H∑V ^H (13)

Wherein, subscript H is the conjugate transpose of complex matrix, and ∑ is that a m * n ties up diagonal matrix, and the element of its principal diagonal is non-negative, i.e. singular value, put in order into:

${\mathrm{\&sigma;}}_{11}^2\&GreaterEqual;{\mathrm{\&sigma;}}_{22}^2\&GreaterEqual;...\&GreaterEqual;{\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="hh"\; filenames="HDA0000058994520000031.png"\; state="\{\{state\}\}">hh</figure-callout>}}^2\&GreaterEqual;0---\left(14\right)$

H=min in the formula (m, n).

When adopting singular value decomposition method to confirm the effective order of matrix, singular value zero constantly corresponding with effective order of matrix.Because R during real data is handled _x(τ) and C _4x(m, n) quilt

With

Replace.Because real data length is shorter, at matrix R _eAnd C _eSingular value often non-vanishing, generally adopt normalization F norm method or normalization singular value method to obtain the effective order of matrix.

It is ARMA (2,2) that the present invention adopts realistic model, for considering that the semi-invariant section is to Special matrix C _eIt is the nothing of the 10000ms composite traces of making an uproar that the effectively influence confirmed of order, experiment produce length, respectively to Special matrix C _eAnd R _eCarry out svd, carry out Monte Carlo random experiments 40 times.It is as shown in table 1 to get five singular value results:

Table 1 is decided the singular value in the rank based on fourth order cumulant and the autocorrelative SVD of sample

Can know that through the data in the comparison sheet 1 unstable based on the singular value that the SVD of fourth order cumulant decides in the rank, it is bigger fluctuate, shows the susceptibility to cutting into slices.Adopt the method for normalization F norm (threshold value gets 0.995) and normalization singular value (threshold value gets 0.005) that singular value is handled, to obtain matrix R _eOr C _eEffective order, as shown in table 2.

Singular value after table 2 normalization is handled

Can know that by table 2 the SVD method based on fourth order cumulant under nothing is made an uproar environment adopts normalization F norm and normalization singular value method all can not obtain the accurate exponent number of AR.In other words, because Special matrix is relatively harsher to the requirement of semi-invariant section, Higher Order Cumulants is to the characteristic of Gaussian noise " blind " crucial point of having no way of that becomes.Compare with Higher Order Cumulants, in the environment that nothing is made an uproar, can confirm the AR exponent number exactly based on the SVD method of autocorrelation function.

Table 3 svd operation time

The svd kind	Based on fourth order cumulant	Based on autocorrelation function
			Operation time/s	6.4688	0.0313

Table 3 is that two kinds of methods are carried out the required time of svd one time, and data can be found out from table, based on the SVD method arithmetic speed of autocorrelation function than based on fast three one magnitude of the SVD method of Higher Order Cumulants.Responsive in view of SVD method to the semi-invariant section based on fourth order cumulant, decide the rank weak effect, and operation time is slow, the present invention will adopt based on the SVD method of autocorrelation function and confirm the AR model order.

The AR model is decided the influence of rank method for considering additive noise and data length; It is the composite traces of 2000ms, 10000ms that experiment generates length, this record is added jamtosignal be respectively behind Gauss's coloured noise of 30% and 50% application and carry out the AR model based on the method for autocorrelation function and decide rank.Adopt Matlab 7.1 to carry out Monte Carlo random experiments 40 times, decide the rank number of success and be 40 times.Hence one can see that, and it is less influenced by data length based on autocorrelative SVD method, has better anti noise ability.

Open SVD method that the prominent personage decides rank with the AR model and be incorporated into MA and decide in the rank, MA is decided rank convert into and find the solution matrix F _eOrder.Matrix F _eWith expansion exponent number q _e＞q is relevant, is defined as:

$F_e=\left[\begin{array}{cccc}{f}_k\left(0\right)& f_k\left(1\right)& ...& f_k\left(q_e\right)\\ f_k\left(1\right)& ...& f_k\left(q_e\right)& \\ .& & & \\ .& & & \\ .& & & \\ f_k\left(q_e\right)& & & 0\end{array}\right]---\left(15\right)$

$f_k\left(\mathrm{\&tau;}\right)=\underset{i=0}{\overset{p}{\mathrm{\&Sigma;}}}a\left(i\right)C_{\mathrm{ky}}(\mathrm{\&tau;}-i,0,...,0)---\left(16\right)$

Theoretical proof rank (F _e)=q+1.

Consider F from another point of view _eBe a upper triangular matrix structure, the confirming of its effective order is equivalent to test:

${\hat{f}}_k^{m+1}\left(m\right)\&NotEqual;0---\left(17\right)$

The maximum integer m that following formula is set up is MA exponent number q, and this test is called diagonal element product (Product Of Diagonal Elements is called for short PODE) test.

F in actual process _eIn element by the sample estimated value

Replace, be to avoid the situation that when confirming the MA exponent number of arma modeling, the overdetermination rank occurs or owe to decide rank, open prominent personage's suggestion SVD and PODE are tested to integrate application, abbreviate " Zhang algorithm " as.

For considering the influence of data length and neighbourhood noise to the Zhang algorithm; It is the composite traces of 2000ms, 10000ms that experiment generates data length; To this record add jamtosignal be respectively Gauss's coloured noise of 30% and 50% after application Zhang algorithm on Matlab 7.1 platforms, carry out Monte Carlo random experiments 40 times, it is as shown in table 4 to decide the rank result.

Table 4Zhang algorithm model is decided the rank result

Can know that by table 4 because Higher Order Cumulants estimated value error under the short time data situation is bigger, it is unstable that the Zhang algorithm is decided rank, near slight fluctuations accurate rank q=2; Under than the long data situation, the rank error also can occur deciding, possibly produce the result on overdetermination rank.The present invention uses quantity of information criterion method correction Zhang algorithm to decide the rank result, and expectation further improves the MA model and decides the rank precision.

The quantity of information that quantity of information criterion method utilizes time series to contain defines different quantity of information criterion functions, gets the parameter p that makes quantity of information minimum, and q is as the exponent number of arma modeling.The BIC criterion function definition is:

$\mathrm{BIC}(p,q)\stackrel{\mathrm{\&Delta;}}{=}\ln {\widehat{\mathrm{\&sigma;}}}_{\mathrm{\&omega;p}}^2+(p+q)\frac{\ln N}{N}---\left(18\right)$

Wherein: N is a data length;

is the variance of linear prediction error, and computing formula is:

${\widehat{\mathrm{\&sigma;}}}_{\mathrm{\&omega;p}}^2\stackrel{\mathrm{\&Delta;}}{=}\underset{i=0}{\overset{p}{\mathrm{\&Sigma;}}}{\hat{a}}_i{\hat{R}}_x(q-i)---\left(19\right)$

Quantity of information criterion method need be attempted various exponent numbers possibilities, if the order information of considering is imperfect, is absorbed in locally optimal solution easily.Can know by table 4,, decide the rank result all near accurate exponent number though that the Zhang algorithm is decided the rank result is unstable, the present invention the Zhang algorithm decide incorporate the quantity of information criterion on the basis, rank and come it is revised, title the method is " MA decides the rank new method ".Concrete performing step is: obtain initial MA exponent number q=N through the Zhang algorithm, calculate q=N ± 1 then respectively, q=N ± 2; Q=N (guarantees q>=1; If q≤0, then abandon) quantity of information criterion function value, get and make the minimum q value of quantity of information criterion function as the MA exponent number.

Get data length 10000ms in the table 4 at present, under jamtosignal 30% situation, decide the rank result decides the rank new method for the case verification MA of q=3 validity.Suppose that AR exponent number p=2 is definite, find the solution q=1,2,3,4,5 o'clock BIC criterion function result is as shown in table 5.

Table 5BIC criterion function result

q	1	2	3	4	5
						BIC(2，q)	0.1332	-0.1096	0.0134	1.2372	0.0038

Can be known that by table 5 during q=2, the BIC criterion function reaches minimum value, be the MA exponent number so we get q=2.Through the correction of quantity of information criterion function, can obtain accurate MA exponent number.

Consider that " MA decides the rank new method " introducing quantity of information criterion function possibly cause the low problem of operation efficiency; Experiment generated data length is the composite traces of 10000ms; It is as shown in table 6 that Zhang algorithm and MA decide the once required time of rank new method (supposing to carry out 5 quantity of information criterion functions judges) operation

Two kinds of algorithms of table 6 compare operation time

Algorithm	The Zhang algorithm	MA decides the rank new method
			t/s	0.0513	0.0514

Can know that by table 6 introducing of quantity of information criterion method is very little to the influence of operation time, does not influence operation efficiency, has verified the high efficiency of new method.

MA is decided the influence of rank new method for considering data length and neighbourhood noise; It is the composite traces of 2000ms, 10000ms that experiment generates length; Add jamtosignal be respectively Gauss's coloured noise of 30% and 50% after New Application Method MA is carried out deciding rank; Adopt Matlab 7.1 to carry out Monte Carlo random experiments 40 times, it is as shown in table 7 to decide the rank result:

Table 7 new method is decided the rank result

Can know by table 7 that invention adopts quantity of information criterion method correction Zhang algorithm to decide the rank result, improve the accuracy that the MA model is decided rank effectively, verify the validity of new method.

Present parameter identification theory comparative maturity, wherein (q, n) algorithm and semi-invariant algorithm application are more based on total least square (SVD-TLS) algorithm, R-C algorithm, the C of svd.The present invention adopts based on the SVD-TLS method of Higher Order Cumulants and estimates the AR model parameter, adopts semi-invariant matrix equation method to estimate the MA model parameter simultaneously.

Try to achieve before matrix R _eRight singular matrix V and effective order p of AR model, following surface construction (p+1) * (p+1) matrix S ^(p):

$S^{\left(p\right)}=\underset{j=1}{\overset{p}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{n+1-p}{\mathrm{\&Sigma;}}}{\mathrm{\&sigma;}}_{\mathrm{jj}}^2{v}_j^i{\left(v_j^i\right)}^H---\left(20\right)$

Wherein: (i j) is the capable j column element of i among the V to

v.

The AR model parameter can be got by following formula:

a(i)＝S ^(-p)(i+1，1)/S ^(-p)(1，1) (21)

Wherein: S ^(-p)Be S ^(p)Inverse matrix.

The MA model parameter estimation method that exists at present mainly contains: R-C algorithm, w microtomy and C (q, n) algorithm etc.The present invention adopts the Higher Order Cumulants algorithm, can suppress Gaussian noise fully in theory, and this method has two kinds of different expression-forms as follows:

$\underset{i=0}{\overset{q}{\mathrm{\&Sigma;}}}b\left(i\right)c_{\mathrm{kx}}^3(q,i-\mathrm{\&tau;})=c_{\mathrm{kx}}(\mathrm{\&tau;},0)c_{\mathrm{kx}}(q,0)c_{\mathrm{kx}}(q,q)---\left(22\right)$

$\underset{i=0}{\overset{q}{\mathrm{\&Sigma;}}}{b}^3\left(i\right)c_{\mathrm{kx}}(q,i-\mathrm{\&tau;})=c_{\mathrm{kx}}(\mathrm{\&tau;},0)c_{\mathrm{kx}}(q,q)/c_{\mathrm{kx}}(q,0)---\left(23\right)$

In real data is handled, must two kinds of forms of algorithm be moved simultaneously.B (i) estimated value in the IF expression (22) is much littler than 1, and method form one is inadvisable; If the b in the formula (23) ³(i) estimated value is much smaller than 1, and method form two is inadvisable.The MA parameter is through b (i), b in estimator (22) and formula (23) formula ³(i) obtain.

Under the prerequisite that obtains the accurate model exponent number, we will estimate model parameter.Be the effect of certificate parameter method of estimation to the different length data; Each experiment all making an uproar the letter energy than being under 40% the situation, generates length at random and is respectively 20000ms, 10000ms; 5000ms; The theogram of 2000ms is extracted wavelet by the model parameter estimation method that the present invention adopts, and the result is as shown in table 9.

Can know that by table 9 the AR parameter all can obtain estimated result preferably under long data and short data situation; (5000ms, estimated result and actual value have than large deviation in the time of 2000ms) and the MA parameter is in short time data.The wavelet waveform that but under the different data lengths situation, extracts has continuity preferably, and is as shown in Figure 2.

Seismic wavelet parameter estimation result under table 9 different data lengths

For investigating the influence of neighbourhood noise to the model parameter estimation method; Each experiment all generates the theogram that length is 5000ms; It is 10% than (NSR) that this record is added the letter energy of making an uproar respectively; Utilize the inventive method antithetical phrase wave pattern to carry out parameter estimation behind Gauss's coloured noise of 30%, 50%, 100%.

Seismic wavelet parameter estimation result under the table 10 varying strength Gauss coloured noise environment

Can know by table 10; Along with the increase AR parameter estimating error of noise intensity slightly increases, even the energy ratio of noise and signal reaches 1: 1, the AR parameter estimation still can obtain good result; Verified that the SVD-TLS method has good anti-noise and holds the ability of making an uproar, can in suitable noise circumstance, use.Because Higher Order Cumulants is responsive to data length; Under the 5000ms data length; MA parameter and actual value have certain deviation, and still along with the increase of noise intensity, the MA estimates of parameters changes little; This is because Higher Order Cumulants has the characteristic of " blind " to Gaussian noise, can suppress the cause of Gaussian noise effectively.

As shown in Figure 3, the seismic wavelet that varying strength Gauss coloured noise environment extracts down has good continuity, has verified that method for parameter estimation of the present invention has certain practicality.

Accurately deciding rank is the keys that guarantee wavelet precision and operation efficiency.The wavelet model that emulation experiment adopts is ARMA (2,2), so, get ARMA (1,2) as owing to decide the rank research object, get ARMA (4,5) as overdetermination rank research object.It is that the nothing of 5000ms is made an uproar behind the seismologic record that ARMA (2,2) wavelet model generated data length is adopted in experiment, utilizes respectively to owe to decide rank and overdetermination rank model carries out wavelet extraction.

Can know by Fig. 4, owe to decide the serious distortion of wavelet of extracting under the situation of rank; Wavelet and the original wavelet extracted under the situation of overdetermination rank are very approaching, but can increase rapidly along with the increase of exponent number operation time, make efficiency of algorithm reduce.

Be the validity of checking the inventive method under known mixed-phase seismic wavelet and non-Gauss's reflection coefficient condition; Get the mixed phase wavelet sequence that contains 28 sample datas: { 0.0060 0.0058 0.0031 0.0023-0.0016 0.0009 0.0008-0.0003-0.0004 0.0001 0.0002-0.0000-0.0001-0.4092 0.3186 1.0000-0.0523-0.4435-0.0158 0.1933 0.0239-0.0827-0.0177 0.0348 0.0108-0.0143 0.0108-0.0143} is shown in Fig. 5-a; The reflection coefficient sequence of statistics independent same distribution (Bernoulli Jacob's distribution) adopts matrix of coefficients to be expressed as: { (0 4 0.3888), (0 20-0.24159), (0 25 0.044189); (0 30-0.01322), (0 53-0.1 1698), (0 58 0.13312); (0 62 0.88258); (0 94 0.14891), (0 96-1.5124) }, shown in Fig. 5-(b); Adopt Robinson seismic trace convolution model to generate composite traces, shown in Fig. 5-(c).

Adopt the inventive method from Fig. 5-(c) extract seismic wavelet the theogram, the wavelet model be ARMA (2,2), model parameter is: a=[1-0.0524 0.4604]; B=[1-0.7752-1.6249]; Relatively extract wavelet and original wavelet is as shown in Figure 6, extract wavelet and original wavelet and have good consistance, thereby verified the validity of the inventive method.

To a data length is 201 roads, and sampling rate is 1ms, and record length is the seismic section of 2s, uses the inventive method and carries out wavelet extraction, and the result is as shown in Figure 7.

In the real data processing procedure, because reflection coefficient sequence departs from assumed condition, the waveform wild effect can appear in the wavelet that some seismic trace extracts.For this reason, adopt multiple tracks to unite and find the solution the equal method of making even to obtain than stable result.Fig. 8 is some adjacent seismic trace application the inventive method of actual seismic section is carried out wavelet extraction, and parameter is averaged the back gained.

Can find out that by above analysis the present invention adopts arma modeling to describe wavelet, model is carried out accurately deciding rank, adopt based on the SVD-TLS method of Higher Order Cumulants and estimate the AR model parameter, adopt semi-invariant matrix equation method to estimate the MA model parameter simultaneously.Handle through data simulation and real data, the method that proves among the present invention to be proposed can precise and high efficiency extracts seismic wavelet.

Description of drawings

Fig. 1, the original wavelet of theogram

The wavelet waveform that extracts under Fig. 2, the different data lengths

Fig. 3, varying strength Gauss coloured noise environment extract the wavelet waveform down

Fig. 4, non-ly accurately decide to extract wavelet under the situation of rank

Fig. 5, emulation seismologic record building-up process

Fig. 6, original wavelet and extraction wavelet comparison diagram

Fig. 7, actual seismic section

The wavelet waveform that extracts in Fig. 8, the actual seismic data

Fig. 9, based on the high-resolution seismic exploration wavelet extraction method flow of high-order statistic and arma modeling

Embodiment

The present invention proposes high-resolution seismic exploration wavelet extraction method based on high-order statistic and arma modeling.The characteristics of this method are: break through the seismic wavelet extraction technology of MA model hypothesis, the arma modeling of first parameter being simplified is applied to the seismic wavelet parametric modeling.The wavelet arma modeling is decided the rank method to be furtherd investigate; At first adopt based on autocorrelative singular value decomposition method and confirm AR part exponent number; Then the quantity of information criterion function being incorporated Higher Order Cumulants MA decides in the method for rank; Proposed a kind of new MA and decide the rank method, that has realized wavelet arma modeling MA part accurately decides rank.Employing is estimated the AR model parameter based on the SVD-TLS method of Higher Order Cumulants, adopts semi-invariant matrix equation method to estimate the MA model parameter, has finally realized the high-precision seismic wavelet extraction of high-level efficiency.Implementing procedure of the present invention is as shown in Figure 9.

The present invention implements according to following steps:

Step 1: set up the arma modeling of seismic wavelet, satisfy the cause and effect mixed-phase hypothesis of seismic wavelet;

Step 2: adopt and find the solution AR exponent number p based on the autocorrelative singular value decomposition method of sample; Make up autocorrelation matrix

M>=p wherein _e>=p, q _e>=q, and q _e-p _e>=q-p finds the solution rank (R _e)=p.

Step 3: adopt " Zhang algorithm " MA exponent number according to a preliminary estimate; Make up matrix

Q wherein _e＞q, and

Find the solution rank (F _e)=q+1.And F _eBe a upper triangular matrix structure, the confirming of its effective order is equivalent to the PODE test:

Make the maximum integer m of its establishment be MA exponent number q.The two is combined the preliminary MA exponent number q of confirming.

Step 4: that uses quantity of information criterion check " Zhang algorithm " decides the rank result; Obtain initial MA exponent number q=N through the Zhang algorithm, calculate q=N ± 1 then respectively, q=N ± 2, the quantity of information criterion function value of q=N (guaranteeing q>=1, if then abandon q≤0) is got and is made the minimum q value of quantity of information criterion function as the MA exponent number.

Step 5: adopt based on the SVD-TLS method of Higher Order Cumulants and estimate the AR model parameter, adopt semi-invariant matrix equation method to estimate the MA model parameter simultaneously.

Claims (2)

1. based on the high-resolution seismic exploration wavelet extraction method of high-order statistic and arma modeling, it is characterized in that the arma modeling of first parameter being simplified is applied in the middle of the seismic wavelet modeling; Employing is confirmed AR part exponent number based on autocorrelative singular value decomposition method, then the quantity of information criterion function is incorporated Higher Order Cumulants MA and decides in the method for rank, proposes a kind of new MA and decides the rank method, improves wavelet arma modeling MA and decides the rank accuracy rate; Adopt SVD-TLS algorithm and The Cumulant Method Using estimator wave parameter respectively; Under the prerequisite that guarantees the wavelet precision, reduce model order as much as possible, improve operation efficiency, finally realize the high-precision seismic wavelet extraction of high-level efficiency.

2. based on the high-resolution seismic exploration wavelet extraction method of high-order statistic and arma modeling, it is characterized in that this method contains following steps successively:

Step (1) is set up the stingy arma modeling of seismic wavelet parameter, satisfies the cause and effect mixed-phase hypothesis of seismic wavelet;

Step (2) adopts and finds the solution AR exponent number p based on the autocorrelative singular value decomposition method of sample:

Make up autocorrelation matrix

M>=p wherein _e>=p, q _e>=q, and q _e-p _e>=q-p finds the solution rank (R _e)=p;

Step (3) adopts " Zhang algorithm " MA exponent number according to a preliminary estimate:

Make up matrix

Q wherein _e＞q, and

Find the solution rank (F _e)=q+1, and F _eBe a upper triangular matrix structure, the confirming of its effective order is equivalent to the PODE test

Make the maximum integer m of its establishment be MA exponent number q, the two is combined the preliminary MA exponent number q of confirming;

What step (4) was used quantity of information criterion check " Zhang algorithm " decides the rank result:

Obtain initial MA exponent number q=N through the Zhang algorithm, calculate q=N ± 1 then respectively, q=N ± 2, the quantity of information criterion function value of q=N (guaranteeing q>=1, if then abandon q≤0) is got and is made the minimum q value of quantity of information criterion function as the MA exponent number;

Step (5) adopts based on the SVD-TLS method of Higher Order Cumulants estimates the AR model parameter, adopts semi-invariant matrix equation method to estimate the MA model parameter simultaneously.

Images