CN101876546B - MEMS (Micro Electronic Mechanical System) gyro data processing method based on wavelet threshold de-noising and FAR (Finite Automaton Recognizable) model - Google Patents

Abstract

The invention discloses an MEMS (Micro Electronic Mechanical System) gyro data processing method based on wavelet threshold de-noising and an FAR (Finite Automaton Recognizable) model. For solving the problem that large amount of noise exists in the MEMS gyro output signal, the wavelet threshold de-noising method is used to process the gyro output signal to filter the noise and improve the signal-noise ratio; aiming at a static signal formed by de-noising gyro wavelet, a polynomial fitting method is used for the compensation to the gyro deterministic drift, wherein a residual error after the compensation is a random drift of the gyro, that is to say, a sample sequence required by modeling the random drift of the gyro is obtained; and aiming at modeling the random drift of the MEMS gyro at higher precision, the FAR model is used for modeling the random drift. The invention solves the problem of high output noise of the MEMS gyro, improves the signal-noise ratio and can accurately model the random drift of the gyro, thereby improving the output precision of the MEMS gyro.

Description

MEMS gyro data disposal route based on wavelet threshold denoising and FAR model

Technical field

The present invention relates to MEMS (Micro Electro Mechanical System; MEMS) gyro output data process field is specifically related to a kind of MEMS gyro data splitting disposal route based on wavelet threshold denoising and FAR (autoregression of Functional coefficient Autoregressive. function coefficients) model.

Background technology

Along with the development of MEMS technology, the MEMS gyro has obtained application more and more widely with characteristics such as its size are little, in light weight, cost is low, power consumption is little in low-cost inertial navigation system.But receive the restriction of manufacturing process and technical merit, the traditional relatively gyro of MEMS Gyro Precision is low at present, and its error becomes the main cause of INS errors.The MEMS gyro error generally can be divided into calibration factor error, deterministic drift, random drift and random noise.The compensation of calibration factor sum of errors deterministic drift can realize through calibration technique.So emphasis that is treated as the processing of MEMS gyro data of noise and random drift.

Contain a large amount of noises in the MEMS gyro output signal, need at first carry out denoising.Mean filter is the simplest signal antinoise method, and its algorithm is realized simple, and denoising effect is better, but generally only is applicable to static state or low current intelligence; Finite Impulse Response filter has been inherited the characteristics of analog filter, and available fast Fourier transform (FFT) algorithm realizes trap signal, has improved operation efficiency greatly.The FIR wave filter carries out filtering and noise reduction to signal to be accomplished in frequency domain, relies on the different spectral characteristic of signal and noise to realize noise filtering, all is applicable to quiet, Dynamic Signal denoising.But denoising effect is general, and it is good to be not so good as mean filter.

Wavelet transformation is the signal analysis and processing instrument that develops rapidly in recent years; It has overcome the frequecy characteristic that Fourier transform can only expression signal but can not reflect the defective of local message on the time domain; And have the partial analysis characteristic and the multiresolution analysis characteristic of time and frequency simultaneously, obtained at aspects such as compression of images, signal filtering and feature extractions to use widely at present.The wavelet threshold denoising algorithm is the basis with the wavelet transformation, based on the corresponding different qualities that wavelet coefficient had after signal and the noise process wavelet decomposition, through selected appropriate threshold criterion, wavelet coefficient is handled, and can realize the noise separation well.The wavelet threshold denoising effect is superior to mean filter and FIR wave filter greatly; And all be suitable for for static, Dynamic Signal; Therefore use it for the MEMS gyro data and handle, not only can be used for the denoising of gyro static state and Dynamic Signal, and can remove noise effectively; Signal to noise ratio (S/N ratio) is improved greatly, for follow-up random drift modeling lays the foundation.

The random drift of MEMS gyro generally shows as than the complex random process, has uncertainty, can't compensate through demarcating.And generally be used for constituting integrated navigation system with GPS by the miniature inertial navigation system that MEMS gyro and mems accelerometer constitute, utilize the long-term high precision of GPS to remedy the problem that the miniature inertial navigation system error accumulates in time; Utilize miniature inertial navigation system short-term high precision to remedy the shortcoming that gps signal is subject to disturb, be prone to losing lock.Carry out the two information fusion through junction filter, obtain reliable, accurate localization output, realize having complementary advantages.Therefore, in integrated navigation, can the random drift of MEMS gyro be modeled as certain probabilistic model, with its state equation as junction filter, the estimation and the compensation of random drift are carried out in the observed quantity that relies on GPS to provide; When gps signal was lost, the probabilistic model that can pass through to be set up carried out the prediction and the real-Time Compensation of the drift of MEMS Gyro Random, thereby reduces the gyro output error, finally can improve the bearing accuracy of integrated navigation system.

First-order Markov process is generally adopted in early stage MEMS Gyro Random drift modeling, because of it can describe many physical accidental processes more exactly, and has relative simple mathematical equation.And its modeling process is simple, the modeling time is short.But the parameter of first-order Markov process confirms to need to rely on the autocorrelation function of sample; And there is very big difference in form in the autocorrelation function curve of the autocorrelation function of first-order Markov process and the drift of MEMS Gyro Random; Therefore the drift of MEMS Gyro Random is modeled as first-order Markov process and can has the model error problem, modeling accuracy is restricted.

Be used for the model error problem that MEMS Gyro Random drift modeling exists for solving first-order Markov process, the S.Nassar of Canadian CALGARY in 2003 is used for AR (Autoregressive, autoregression) model the random drift modeling of gyro first.Experimental result shows: compare with the single order markov, the bearing accuracy when adopting the AR model modeling that gps signal is lost improves 15% ~ 35%.After this, the AR model is widely used in the Gyro Random drift modeling gradually.The AR model belongs to a kind of linear model in the time series analysis method, is to utilize the dependence between the adjacent time series to carry out analysis and modeling, and the exponent number of model is variable.Compare with first-order Markov process, it can reflect time series better by changing exponent number, therefore has bigger modeling flexibility and the modeling accuracy of Geng Gao.But Gyro Random drift often exists non-linear and stationarity not, adopt linear modelling can introduce certain error, and the desired stationarity of AR model generally is difficult to real satisfied.Function coefficients autoregression (FAR) model is a kind of based on nonparametric, nonlinear statistical model, and its coefficient is the function that relies on variable.Compare with AR, it seeks proper model in a bigger model family, can effectively reduce the error of introducing because of Model Selection is improper.At present, the FAR model has been succeedd in fields such as computing machine, electric power, finance, biomedicines and has been used, but also is not used to as yet in the random drift modeling of MEMS gyro.

Therefore, wavelet threshold denoising and FAR model are combined, constitute a kind of novel MEMS gyro data splitting disposal route; Can effectively solve the big problem of MEMES gyro output noise; Improve signal to noise ratio (S/N ratio), and can carry out accurate modeling, thereby improve MEMS gyro output accuracy the Gyro Random drift.

Summary of the invention

The objective of the invention is to: the deficiency that overcomes prior art; A kind of MEMS gyro data disposal route based on wavelet threshold denoising and FAR model is provided; Solve the big problem of MEMES gyro output noise; Improved signal to noise ratio (S/N ratio), and can carry out accurate modeling, thereby improved MEMS gyro output accuracy the Gyro Random drift.

The objective of the invention is to realize through following technical scheme: based on the MEMS gyro data disposal route of wavelet threshold denoising and FAR model, performing step is following:

(1) at first, for solving the problem that contains much noise in the MEMS gyro output signal, adopt the wavelet threshold denoising method that gyro output signal is handled, filtering noise improves signal to noise ratio (S/N ratio), and the step of said wavelet threshold denoising method is following:

(1.1) decomposable process: selected a kind of wavelet function, and utilize the level of confirming decomposition based on the adaptive algorithm of singular spectrum analysis, the original output signal of gyro is carried out N layer wavelet decomposition;

(1.2) threshold value process: to decomposing each layer wavelet coefficient that obtains, select threshold value according to certain criterion, mould keeps greater than the wavelet coefficient of threshold value, and mould is made as zero less than the wavelet coefficient of threshold value;

(1.3) process of reconstruction: recover original signal according to each layer coefficients after the denoising through wavelet reconstruction, obtain the gyro output after the denoising;

(2) to the stationary singnal after the denoising of gyro small echo, adopt the method for fitting of a polynomial to carry out the compensation of gyro deterministic drift, the residual error after the compensation is the random drift of gyro, promptly obtains the required sample sequence of Gyro Random drift modeling;

(3) purpose for realizing MEMS Gyro Random drift motion degree of precision is carried out modeling adopts the FAR model to the random drift modeling, and said modeling process is following:

(3.1) with the MEMS Gyro Random drift data that obtains in the step (2) as modeling sample, with sample data carry out zero-mean, normalization is handled, and the Gyro Random wander sequences { X after will handling ₁, X ₂..., X _nMaximal value be designated as: $u_{\mathrm{Max}}=\underset{t}{\mathrm{Max}}\left(X_t\right),$ Minimum value is designated as $u_{\mathrm{Min}}=\underset{t}{\mathrm{Min}}\left(X_t\right),$ If { u _j| j=0,1 ... N} is Along ent set, wherein u ₀=u _Min, u _N=u _Max, rely on Along ent can the interval uniformly-spaced be divided into N+1 segment:

[u _N-b ₀，u _N)，[u ₀，u ₀+b ₀)，[u ₁-b ₀，u ₁+b ₀)，[u ₂-b ₀，u ₂+b ₀)，…，[u _N-1-b ₀，u _N-1+b ₀)，

Wherein: b ₀=(u _Max-u _Min)/2N,

The higher limit of N can be not less than 2p with reference to the sample point number in arbitrary segment and come givenly, and p is a model order;

(3.2) at model order p, under the situation that the parameter d of appointment dependence variable is confirmed, behind selected window function width b, minimize error of fitting, try to achieve each interval upward local linear approximation coefficient of coefficient function by weighted least-squares;

(3.3) model parameter N, p, choosing of d and b can realize through defining average mean square prediction error APE, lets N=1 successively ..., N _MaxP=1 ..., p _MaxD=1 ..., p; B=b ₀..., 1, N _MaxBe the maximal value of interval division N, p _MaxMaximal value (p for exponent number p _MaxGenerally be taken as 5 ~ 8), b ₀=(u _N-u ₀)/2N calculates APE respectively, makes the minimum N of APE value, and p, d, b calculate the local linear approximation coefficient α of the FAR coefficient function of this moment as the optimal parameter of model _{I, j}, β _{I, j}(i=1 ... P; J=0,1 ... N), set up the FAR model of final Gyro Random drift.

The present invention's beneficial effect compared with prior art is mainly reflected in: compare with existing MEMS gyro data disposal route; The present invention proposes a kind of novel MEMS gyro data combination treatment method; This method had both efficiently solved the big problem of MEMS gyro output noise, had solved the problem of its random drift modeling again.And overcome existing method, thereby can improve MEMS gyro output accuracy greatly, the positioning performance of final enhanced navigation system to the big problem of MEMS Gyro Random drift modeling error.

Description of drawings

Fig. 1 is a data splitting disposal route synoptic diagram of the present invention;

Fig. 2 is a wavelet threshold denoising method flow diagram of the present invention;

Fig. 3 is that the decomposition number of plies self-adaptation based on singular spectrum analysis of the present invention is chosen algorithm flow chart;

Fig. 4 is a FAR modeling method process flow diagram of the present invention;

Fig. 5 is the definition mode synoptic diagram of average mean square prediction error.

Embodiment

Data splitting disposal route among the present invention is that the noise that exists to MEMS gyro output signal is big, and the problem that random drift is serious proposes.Wherein, wavelet threshold denoising can be used for the noise remove of gyro static state, Dynamic Signal; And the random drift modeling of FAR model needs the static state output based on gyro, i.e. modeling under static state, but can under current intelligence, use after the modelling.

Figure is as shown in Figure 1 for this data splitting process flow, and its practical implementation step is following:

1, at first gather the original output signal of MEMS gyro of a period of time, beginning a period of time is required to be stationary singnal.

2, the original output signal of MEMS gyro is carried out wavelet threshold denoising, Fig. 2 is the process flow diagram of wavelet threshold denoising, and the practical implementation process is following:

2.1 according to the actual conditions of denoising, select to have the suitable support length (5 ~ 9) and the wavelet basis function of higher vanishing moment, choose the db5 small echo among the present invention;

Confirm the optimal Decomposition number of plies 2.2 adopt based on the decomposition number of plies adaptive algorithm of singular spectrum analysis, its algorithm flow chart is as shown in Figure 3, and concrete steps are following:

1) the initial number of plies j=1 that decomposes of order;

2) gyro output signal is carried out j layer wavelet decomposition, obtain each layer scale coefficient c _{J, k}With wavelet coefficient d _{J, k}, k=1 ... [n/ (2j)], n is for needing the sample number of denoising, and [n/ (2j)] representes n/ (2j) round numbers;

3) with wavelet coefficient d _{J, k}It is m that sequence is embedded into dimension by certain delay τ ₀Phase space, by Takens embedding theorems restructural attractor track matrix:

M=[n/ (2j)] wherein, the selection of delay time T and wavelet coefficient sequence { d _{J, k}The auto-correlation degree relevant, the auto-correlation degree of sequence is more little, then the value of τ is more little; And m ₀Size determined the obvious degree that singular spectrum changes, m ₀Big more, the singular spectrum intensity of variation is obvious more.τ and m here ₀All be taken as 10.

To attractor track matrix

Carry out singular spectrum analysis,, calculate singular value e promptly through svd _iAnd singular spectrum $S_i=e_i/\left(\frac{1}{{\mathrm{<figure-callout\; id="0"\; label="m"\; filenames="H2009102412382E00021.png,H2009102412382E00041.png"\; state="\{\{state\}\}">m</figure-callout>}}_0}\underset{i=1}{\overset{m_0}{\mathrm{\&Sigma;}}}{e}_i\right),i=\mathrm{1,2},<figure-callout\; id="0"\; label="...\; m"\; filenames="H2009102412382E00021.png,H2009102412382E00041.png"\; state="\{\{state\}\}">.</figure-callout>..m_{}{}_0,$ And the slope of calculating singular spectrum $K={\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="H2009102412382E00041.png"\; state="\{\{state\}\}">e</figure-callout>}}_1/e_{m_0}$ (ratio of maximum, minimum singular value), e ₁Be maximum singular value,

Be minimum singular value;

4) compare singular spectrum slope K and empirical value $\mathrm{Th}=\frac{2}{1+\mathrm{Exp}(-{\mathrm{<figure-callout\; id="1"\; label="e"\; filenames="H2009102412382E00041.png"\; state="\{\{state\}\}">e</figure-callout>}}_1*4)}+0.05$ Size, if K＜th, make j add 1, continue execution in step 2); Otherwise execution in step 5);

5) the best number of plies N that decomposes is taken as j-1.

2.3, select the appropriate threshold criterion according to signal and noise behavior.Big according to the gyro output noise among the present invention, characteristics such as signal to noise ratio (S/N ratio) is lower are selected soft-threshold function and rigrsure threshold value criterion, have confirmed corresponding threshold value th _j

Soft-threshold is a kind of processing policy to wavelet coefficient, and the soft-threshold function can use following expression formula to describe:

$\mathrm{th}\left(d\right)=\left\{\begin{array}{cc}\mathrm{sign}\left(x\right)\left(\right|x\left|\right)-\mathrm{th}& \left|x\right|>\mathrm{<figure-callout\; id="0"\; label="th"\; filenames="H2009102412382E00021.png,H2009102412382E00041.png"\; state="\{\{state\}\}">th</figure-callout>}\\ 0& \left|x\right|\&le;\mathrm{th}\end{array}\right.$

Wherein th is the size of threshold value, and sign (x) represents the symbol of wavelet coefficient x.

The rigrsure threshold value criterion confirms that specifically algorithm is:

(a) take absolute value the element in the pending sequence, ascending ordering with each element square, obtains new sequence Y then, and its length is the length n of former sequence;

(b) corresponding each element subscript (being the sequence number of element) k is if get k the element Y that thresholding is a sequence _kSquare root, then the risk algorithm is:

$\mathrm{Risk}\left(k\right)=\frac{n-2k+\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{Y}_j+(n-k)*Y_{n-k}}{n}$

Y wherein _k, Y _j, Y _N-kBe respectively k element, a j element, the n-k element of sequence Y.

According to following formula, draw corresponding to the risk curve of different value of K, find out minimum risk point and corresponding with it k value, then its threshold value th does $\mathrm{Th}=\sqrt{Y_k}$

2.4 utilize selected wavelet basis function that gyro output signal is carried out N layer wavelet decomposition, obtain each layer scale coefficient c _{J, k}With wavelet coefficient d _{J, k}, its computation process is following:

If the scaling function that selected wavelet basis function is corresponding is established gyro output signal for be f (t),

calculate initial scale coefficient

earlier for needing the number of samples of denoising.By iterative algorithm, can calculate the scale coefficient of each layer: $c_{j,k}=\underset{m=1}{\overset{[n/2(j-1)]}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{h_{m-2k}}{c}_{j-1,m},$ K=1,2 ... [n/2j], c _J-1, m is the scale coefficient on the j-1 decomposition scale, each layer wavelet coefficient: $d_{j,k}=\underset{m=1}{\overset{[n/2(j-1)]}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{g_{m-2k}}{c}_{j-1,m},$ K=1,2 ... [n/2j].Wherein: [n/2j] representes the n/2j round numbers; h _M-2k, g _M-2kBe respectively wavelet basis function corresponding low-pass filter coefficients and Hi-pass filter coefficient,

g _m-2k＝(-1) ^m-2kh _1-(m-2k)

Wherein:

2.5 with each layer wavelet coefficient d _{J, k}With corresponding threshold value th _jCompare, if | d _{J, k}|＜th _j, d then _{J, k}Be made as zero; If | d _{J, k}|>=th _j, d then _{J, k}Handle by the soft-threshold function;

2.6 according to c _{J, k}With the d after the processing _{J, k}Utilize wavelet reconstruction algorithm reconstruction signal, promptly obtain the gyro output after the denoising;

Restructing algorithm is following: $c_{j-1,m}=\underset{k=1}{\overset{[n/2j]}{\mathrm{\&Sigma;}}}{h}_{m-2k}{c}_{j,k}+\underset{k=1}{\overset{[n/2j]}{\mathrm{\&Sigma;}}}{g}_{m-2k}{d}_{j,k},m=1,...[n/2(j-1)],$

Finally can get initial gauges coefficient c by above iterative algorithm _{0, k}, k=1,2 ... N then can get the signal after the denoising by wavelet inverse transformation:

Wherein n is the number of samples that carries out denoising.

3, to the gyro static output signal sequence { x behind the process wavelet threshold denoising ₁, x ₂..., x _n, adopt the method for fitting of a polynomial to obtain the determinacy drift of gyro, and compensate.Residual error after the deterministic drift compensation is the random drift { X of gyro ₁, X ₂..., X _n;

If the gyro deterministic drift shows as polynomial form: f (t)=a _kt ^k+ a _K-1t ^K-1+ ... + a ₀, wherein f (t) is the function of time t, t=1 ... N; K is the degree of polynomial, generally gets 2 ~ 4, a _k, a _K-1A ₀Be multinomial coefficient; The random drift X of gyro then _i=x _i-f (i), i=1 ... N.

4, utilize the Gyro Random drift { X of last step acquisition ₁, X ₂..., X _nAs sample, set up the FAR model of MEMS Gyro Random drift.

For zero-mean, normalization time series { X _t, t ∈ N}, the equation of FAR model is:

X _t＝f ₁(X _t-d)X _t-1+…+f _p(X _t-d)X _t-p+ε _t

Wherein, f (.) is a coefficient function, does not stipulate its concrete form (need satisfy flatness); P is the autoregression order; D ∈ 1 ... P}, designated model relies on variable X _T-d{ ε _t, t ∈ N} is independent identically distributed random series, and supposition ε _t～N (0, σ _ε ²); ε _tWith X _T-iSeparate.Thus, with model brief note for FAR (p, d).

Concrete modeling process is following:

4.1 at first to Gyro Random drift sample sequence { X ₁, X ₂..., X _nCarry out normalization, wherein X _t∈ { X ₁, X ₂..., X _n}: $X_t=\frac{X_t-\underset{t}{\mathrm{Min}}\left(X_t\right)}{\underset{t}{\mathrm{Max}}\left(X_t\right)-\underset{t}{\mathrm{Min}}\left(X_t\right)}$

Again it carry out zero-meanization:

$X_t=X_t-\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\&Sigma;}}}{X}_t$

Obtain the random drift sequence { X of normalizing, zero-mean ₁, X ₁..., X _n.

4.2 set the maximal value N of interval division N _Max, the maximal value p of exponent number p _Max, the increase step delta b of window function width b.N, p, d add 1 since 1 successively successively, and b begins to increase (b0 is relevant with the value of N, and it is defined in the following step and provides) successively with step delta b from b0, and circulation carried out for 4.3 and 4.4 steps;

Estimate the coefficient function of FAR model 4.3 adopt local linear regression technique.

1) with the Gyro Random wander sequences { X of normalization, zero-meanization ₁, X ₂..., X _nMaximal value be designated as: $u_{\mathrm{Max}}=\underset{t}{\mathrm{Max}}\left(X_t\right),$ Minimum value is designated as $u_{\mathrm{Min}}=\underset{t}{\mathrm{Min}}\left(X_t\right),$ And note u=X _T-d, u ∈ [u then _Min, u _Max].

2) establish { u _j| j=0,1 ... N} is Along ent set, wherein u ₀=u _Min, u _N=u _MaxRelying on Along ent can be N+1 segment with interval division:

[u _N-b ₀，u _N)，[u ₀，u ₀+b ₀)，[u ₁-b ₀，u ₁+b ₀)，[u ₂-b ₀，u ₂+b ₀)，…，[u _N-1-b ₀，u _N-1+b ₀)，

Wherein: b ₀=(u _Max-u _Min)/2N.

3) with branch { u _j| j=0,1 ... N} is the center, b (b ₀≤b≤1) is the interior { u of regional area of the radius of neighbourhood _Sel-b, u _Sel+ b}, coefficient function f _i(u) can be through the single order Taylor expansion

f _i(u)≈f _i(u _sel)+f′ _i(u _sel)(u-u _sel)＝β _i+α _i.(u-u _sel)，u∈{u _sel-b，u _sel+b}

Wherein: u _Sel∈ { u _j| j=0,1 ... N}, β _i=f _i(u _Sel) and α=f ' _i(u _Sel) be the local linear approximation coefficient of the corresponding coefficient function of FAR model.

4) through minimizing (weighted least-squares)

$\underset{t=p+1}{\overset{n}{\mathrm{\&Sigma;}}}\{{\{X_t-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i+{\mathrm{\&alpha;}}_i.(X_{t-d}-u_{\mathrm{sel}}).X_{t-i}\}}^2.K_b(X_{t-d}-u_{\mathrm{sel}})\}$

Can obtain α _iAnd β _iEstimated value, i=1,2 ..., p.K in the formula _bCan adopt the Epanechnikov window function:

K _b(u)＝(0.75/b).[1-(u/b) ²] ₊

Its effect is only to allow sample data (X _t, X _T-1... X _T-p) | X _T-d∈ (u _Sel-b, u _Sel+ b) computing is estimated in participation, and rely on variable-value the closer to the center, and respective sample is to estimating α _iAnd β _iContribution maximum.

5) then coefficient function at segment { u _Sel-b ₀, u _Sel+ b ₀(work as u _Sel=u ₀Or u _Sel=u _NThe time, segment does

[u ₀, u ₀+ b ₀) or [u _N-b ₀, u _N)) on linear-apporximation do

f _i(u)＝β _i+α _i·(u-u _sel)，u∈(u _sel-b ₀，u _sel+b ₀)，i＝1，2，…，p

u _Sel∈ { u _j| j=0,1 ... N}, { u _j| j=0,1 ... N} is the Along ent set

Can obtain the local linear-apporximation of coefficient function on all segments according to said method.

4.4 calculate the corresponding average mean square prediction error (APE) of model parameter this moment, Fig. 5 is the definition mode synoptic diagram of average mean square prediction error (APE), its computing method are:

1) be the Gyro Random wander sequences of n to length, selected positive integer L and m require n＞mL, and n-mL random drift sample carries out the FAR modeling before at first adopting, and a following L sample is predicted; Use L sample of preceding n-(m-1) to carry out the FAR modeling again, and L sample below the prediction, by that analogy, carry out m modeling and prediction altogether;

2) calculate every section mean square prediction error respectively:

${\mathrm{MSPE}}_j(p,d,b,N)=\frac{1}{L}\underset{t=n-\mathrm{jL}+1}{\overset{n-\mathrm{jL}+L}{\mathrm{\&Sigma;}}}{[X_t-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{f}_i^{j,b,N}\left(X_{t-d}\right)X_{t-i}]}^2$

Wherein: p is a model parameter, d ∈ 1 ... P} is used for designated model and relies on variable X _T-d, b is the window function width, N is the interval division number, f _i ^{J, b, N}(X _T-d) (i=1 ... P) be the coefficient function of FAR model, X _t, X _T-d, X _T-iBe respectively the element that is designated as t, t-d, t-i in the Gyro Random drift sample sequence down, m is the Gyro Random wander sequences prediction number of times described in Fig. 5, and L is the number of samples of each prediction, j=1,2 ..., m.

3) m mean square prediction error averaged be average mean square prediction error (APE):

$\mathrm{APE}(p,d,b,N)=\frac{1}{m}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{MSPE}}_j(p,d,b,N)$

4.5 work as N=N _Max, p=p _Max, d=p, during b=1, might situation APE (p, d, b, N) value calculating finishes;

4.6 the APE of value minimum (p, d, b, N) pairing parameter p, d, b, N is as the optimal parameter of FAR model;

4.5 storage parameter p at this moment, d, b, N calculates the alpha of the local linear-apporximation of FAR model coefficient function at this moment _{I, j}, β _{I, j}(i=1 ... P; J=0,1 ... N).Zero-meanization before considering need be done reduction and handle, so the FAR model of final Gyro Random drift of setting up is:

X _t+1＝f ₁(u)(X _t-μ _X)+f ₂(u)(X _t-1-μ _X)+…+f _p(u)(X _t-p+1-μ _X)+μ _X

$u=\frac{X_{t-d+1}-{\mathrm{\&mu;}}_X}{\underset{t}{\max }\left(X_t\right)-\underset{t}{\min }\left(X_t\right)}$

μ wherein _XBe former Gyro Random wander sequences { X ₁, X ₂..., X _nAverage, f ₁(u), f ₂(u) ..., f _p(u) be the coefficient function of FAR model, X _T+1, X _t, X _T-1, X _T-p+1, X _T-d+1The Gyro Random that is respectively the t+1 moment in the Gyro Random wander sequences, the t moment, the t-1 moment, the t-p+1 moment, the t-d+1 moment is drifted about,

Be former Gyro Random wander sequences { X ₁, X ₂..., X _nMaximal value,

Be former Gyro Random wander sequences { X ₁, X ₂..., X _nMinimum value.

In sum, the present invention at first adopts the wavelet threshold denoising method to come the much noise in effective filtering MEMS gyro output signal, improves its output signal-to-noise ratio; Utilize fitting of a polynomial to carry out the compensation of gyro deterministic drift after the denoising; Utilize the FAR model to carry out the modeling of MEMES Gyro Random drift at last, obtain the accurate model of its random drift, for the compensation of random drift lays the foundation.Finally can improve the bearing accuracy of navigational system, be a kind of novel, effective MEMS gyro data splitting disposal route.

Below only be concrete exemplary applications of the present invention, protection scope of the present invention is not constituted any limitation.The technical scheme that forms but the application that its expanded application is handled in all gyro datas, all employing equivalents or equivalence are replaced all drops within the rights protection scope of the present invention.The part that the present invention does not set forth in detail belongs to techniques well known.

Claims (2)

1. based on the MEMS gyro data disposal route of wavelet threshold denoising and FAR model, it is characterized in that performing step is following:

(1) adopt the wavelet threshold denoising method that gyro output signal is handled, filtering noise improves signal to noise ratio (S/N ratio), and the step of said wavelet threshold denoising method is following:

(1.1) decomposable process: selected a kind of wavelet function, and utilize the level of confirming decomposition based on the adaptive algorithm of singular spectrum analysis, the original output signal of gyro is carried out N layer wavelet decomposition;

(1.2) threshold value process: to decomposing each layer wavelet coefficient that obtains, select threshold value according to certain criterion, mould is handled by the soft-threshold function greater than the wavelet coefficient of threshold value, and mould is made as zero less than the wavelet coefficient of threshold value;

(1.3) process of reconstruction: recover original signal according to each layer coefficients after the denoising through wavelet reconstruction, obtain the gyro output after the denoising;

(2) to the stationary singnal after the denoising of gyro small echo, adopt the method for fitting of a polynomial to carry out the compensation of gyro deterministic drift, the residual error after the compensation is the random drift of gyro, promptly obtains the required sample sequence of Gyro Random drift modeling;

(3) adopt the FAR model to the random drift modeling, said modeling process is following:

(3.1) with the MEMS Gyro Random drift data that obtains in the step (2) as modeling sample, with sample data carry out zero-mean, normalization is handled, and the Gyro Random wander sequences { X after will handling ₁, X2 ..., X _nMaximal value be designated as:

Minimum value is designated as

X _t∈ { X ₁, X ₂..., X _n, establish { u _j| j=0,1 ... N} is Along ent set, wherein u ₀=u _Min, u _N=u _Max, rely on Along ent can the interval uniformly-spaced be divided into N+1 segment:

[u _N-b ₀, u _N), [u ₀, u ₀Ten b ₀), [u ₁-b ₀, u ₁+ b ₀), [u ₂-b ₀, u ₂+ b ₀) ..., [u _N-1-b ₀, u _N-1+ b ₀),

Wherein: b ₀=(u _Max-u _Min)/2N, the higher limit of N can be not less than 2p with reference to the sample point number in arbitrary segment and come givenly, and p is a model order;

(3.2) at model order p, under the situation that the parameter d of appointment dependence variable is confirmed, behind selected window function width b, minimize error of fitting, try to achieve each interval upward local linear approximation coefficient of coefficient function by weighted least-squares;

(3.3) model parameter N, p, choosing of d and b can realize through defining average mean square prediction error APE, lets N=1 successively ..., N _MaxP=1 ..., p _MaxD=1 ..., p; B=b ₀..., 1, N wherein _MaxBe the maximal value of interval division N, p _MaxBe the maximal value of exponent number p, b ₀=(u _N-u ₀)/2N calculates APE respectively, makes the minimum N of APE value, and p, d, b calculate the local linear approximation coefficient α of the FAR coefficient function of this moment as the optimal parameter of model _{I, j}, β _{I, j}(i=1 ... P; J=0,1 ... N), set up the FAR model of final Gyro Random drift;

Adopt the N that makes average mean square prediction error APE value minimum in the said step (3.3), p, d, b are as the model optimal parameter, and its step is following:

The maximal value N of A, setting Gyro Random drift sample sequence interval division N _Max, the maximal value p of exponent number p _Max, the increase step delta b of window function width b, d ∈ 1 ... P}, N, p, d add 1 since 1 successively successively, and b begins to increase successively with step delta b from b0, circulation execution in step B, C, D;

B, be the Gyro Random wander sequences of n to length, selected positive integer L and m require n＞mL, and n-mL random drift sample carries out the FAR modeling before at first adopting, and a following L sample is predicted; Use L sample of preceding n-(m-1) to carry out the FAR modeling again, and L sample below the prediction, by that analogy, carry out m modeling and prediction altogether;

C, calculate every section mean square prediction error respectively:

${\mathrm{MSPE}}_j(p,d,b,N)=\frac{1}{L}\underset{t=n-\mathrm{jL}+1}{\overset{n-\mathrm{jL}+L}{\mathrm{\&Sigma;}}}{[X_t-\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{f}_i^{j,b,N}\left(X_{t-d}\right)X_{t-i}]}^2$

Wherein: p is a model parameter, d ∈ 1 ... P} is used for designated model and relies on variable X _T-d, b is the window function width, N is the interval division number,

Be the coefficient function of FAR model, X _t, X _T-d, X _T-iBe respectively the element that is designated as t, t-d, t-i in the Gyro Random drift sample sequence down, m is the prediction number of times of Gyro Random wander sequences among the B, and L is the number of samples of each prediction, j=1,2 ..., m;

D, m mean square prediction error averaged is average mean square prediction error APE:

$\mathrm{APE}(p,d,b,N)=\frac{1}{n}\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{MSPE}}_j(p,d,b,N);$

E, work as N=N _Max, p=p _Max, d=p, during b=1, might situation APE (p, d, b, N) value calculating finishes;

The APE of F, value minimum (p, d, b, N) pairing parameter p, d, b, N is as the optimal parameter of FAR model.

2. the MEMS gyro data disposal route based on wavelet threshold denoising and FAR model according to claim 1 is characterized in that: said step (1.1) adopts confirms that based on the decomposition number of plies adaptive algorithm of singular spectrum analysis the step of the optimal Decomposition number of plies is following:

A. number of plies j=1 is initially decomposed in order;

B. gyro output signal is carried out j layer wavelet decomposition, obtain each layer scale coefficient c _{J, k}With wavelet coefficient d _{J, k}, k=1 ... [n/2j)], wherein n is for needing the number of samples of denoising, and [n/ (2j)] representes n/ (2j) round numbers;

C. with wavelet coefficient d _{J, k}Sequence is arranged in m by postponing τ ₀Matrix, carry out singular spectrum analysis, calculate singular value e _iAnd singular spectrum

And the slope of calculating singular spectrum

E wherein ₁Be maximum singular value,

Be minimum singular value;

D. the size that compares singular spectrum slope K and empirical value

; If K＜th; Make j add 1, continue execution in step B; Otherwise execution in step E;

E. best decomposition number of plies N is taken as j-1.

Images