CN103888100B - A kind of non-gaussian linear stochaastic system filtering method based on negentropy - Google Patents

Abstract

A kind of non-gaussian linear stochaastic system filtering method based on negentropy.First, for general linear stochaastic system, design wave filter；Secondly, it is contemplated that error equation is a non-gaussian stochastic process, and target function elects the linear combination of estimation error covariance battle array and negentropy as so that it is the probabilistic statistical characteristics of the most complete sign estimation difference；Again, utilize the character of characteristic function, ask for the probability density function of estimation difference, thus obtain the negentropy expression formula of estimation difference；Finally, filter gain, minimization target function are solved；This method can estimate the linear stochaastic system state affected by non-Gaussian noise, can be used for the fields such as inertial navigation, guidance system, target following, signal processing.

Description

A kind of non-gaussian linear stochaastic system filtering method based on negentropy

Technical field

The present invention relates to a kind of non-gaussian linear stochaastic system filtering method based on negentropy, particularly a kind of at random The probability density density function of input has filtering method during the strong non-Gaussian features such as asymmetric, multi-peak, and the method can be used In fields such as inertial navigation, guidance system, target following, signal processing.

Background technology

In recent years, along with the development of aeronautical and space technology, the requirement to aircraft autonomy and quick-reaction capability more comes The highest, it is meant that a lot of complicated estimation tasks must be solved fast and effectively.And filtering theory, it is to solve these problems Key technology.Use which kind of filtering method, from the angle estimation system mode that probability statistics are optimum, in existing hardware condition On the basis of, the control of raising aircraft and the precision of navigation are significant, are the important guarantees of all kinds of Mission Success.20 generation Recording the sixties, Kalman filtering successfully solves the multi-sensor integrated navigation difficult problem in american apollo moonfall, causes The extensive attention of engineering circles.Kalman filtering is a kind of Recursive Filtering algorithm, as long as knowing the estimated value in a upper moment and working as The observation of front state just can calculate the estimated value of current state, it is not necessary to hourly observation or the historical information of estimation, Meet requirement of real-time, be therefore widely used in navigation system.Kalman filter method takes full advantage of for information about And the physical characteristic of the amount of being estimated itself, it is better than least squares estimate, but is only applicable to the line of noise Gaussian distributed Sexual system.Owing to navigation system working environment is more complicated, and there is uncertainty highly, be therefore generally operational in non-gaussian and make an uproar Under acoustic environment.If the most still using kalman filter method can lower navigation accuracy.

For the problems referred to above, on the basis of kalman filtering theory, many advanced filtering methods occur in succession.1993 Year, the Monte Carlo filtering algorithm that Grodon et al. proposes is rooted in Bayesian filter framework as Kalman filtering algorithm, But the most do not require that Posterior distrbutionp meets the form of Gauss distribution, therefore can process non-gaussian filtering problem. Condensation, particle filter, it is all Monte Carlo filtering algorithm title in different field that bootstrap filtering waits, Its core concept is to represent to be asked by the random sample obtaining the sampling of a certain probability distribution and corresponding probability distribution Posterior distrbutionp, when sample number tends to infinity, the unlimited approaching to reality of algorithm is distributed, and the result obtained is by infinite approach optimal solution. But owing to hardly resulting in optional sampling distribution in practice, therefore the precision of Monte Carlo filtering algorithm is the highest and algorithm is multiple Miscellaneous, computationally intensive, poor real.In sum, existing filtering method or suppose random noise Gaussian distributed or Need substantial amounts of sample, such processing method mainly to have following inferior position: for the system affected by non-Gaussian noise, permitted Multi information is included in inside higher order statistical characteristic, and variance can not characterize its probabilistic statistical characteristics completely.Kalman filter method With the minimum criterion of variance of estimaion error, precision will certainly be reduced；Monte Carlo filtering algorithm needs substantial amounts of sample Just can obtain reasonable result, but owing to real system is by various restrictions, hardly result in optional sampling sample, and do not have There is the statistical information making full use of noise.

Summary of the invention

The technology of the present invention solves problem: for non-gaussian linear stochaastic system, overcomes the deficiencies in the prior art, and fills Divide the statistical information utilizing stochastic inputs, it is provided that a kind of filtering method based on negentropy, solve non-gaussian linear stochaastic system Filtering problem, improves the filtering accuracy of system.

The technical solution of the present invention is: a kind of non-gaussian linear stochaastic system filtering method based on negentropy, in fact Existing step is as follows:

The first step, design wave filter:

Discrete-time linear non-gaussian stochastic system for general:

x_k+1=A_kx_k+G_kω_k

y_k=C_kx_k+υ_k

It is constructed as follows the wave filter of form:

${\hat{x}}_{k+1}=A_k{\hat{x}}_k+L_k(y_{k+1}-C_{k+1}{A}_k{\hat{x}}_k)$

Wherein, x_k∈RⁿFor system mode, y_k∈R^mFor measuring output；A_k∈R^n×nFor known state-transition matrix, G_k∈R^{n ×s}For known interference transfer matrix, C_k∈R^m×nFor known output matrix；For system mode x_kEstimated value, L_k∈R^n×mFor Filter gain to be determined；ω_k∈R^s, υ_k∈R^mIt is respectively known to zero-mean, bounded, separate and probability density function Nongausian process noise and measurement noise, its probability density function is respectively Initial vector x₀With ω_k, υ_kSeparate, probability density function isα₁, β₁, α₂, β₂For known real number； R represents real number field；

Definition estimation differenceThen error equation meets:

e_k+1=(I-L_kC_k+1)(A_ke_k-G_kω_k)+L_kυ_k+1

Wherein I is that n × n ties up unit matrix；

Second step, target function elects as:

$V_{k+1}\left(L_k\right)=R_{1k}\∫{\mathrm{\γ}}_{e_{k+1}}(x,L_k)\log ({\mathrm{\γ}}_{e_{k+1}}(x,L_k)/{\mathrm{\γ}}_{k+1}(x,L_k))+R_{2k}{c}_{k+1}$

Wherein, c_k+1Represent estimation difference e_k+1Covariance function, γ e_k+1(x) and γ_k+1X () represents estimation difference respectively e_k+1Probability density function and and e_k+1There is the probability density function of the Gaussian random variable of identical covariance matrix；R_1k, R_2kFor known weight matrix；

Asking for of 3rd step, estimation difference covariance matrix and probability density function:

The covariance matrix of error equation and characteristic function meet following two formulas respectively:

$c(k+1)=(I-L_k{C}_{k+1})A_{k}c\left(k\right)A_k^T{(I-L_k{C}_{k+1})}^T+(I-L_k{C}_{k+1})G_k{c}_{{\mathrm{\ω}}_k}{G}_k^T{(I-L_k{C}_{k+1})}^T+L_k{c}_{{\mathrm{\υ}}_k}{L}_k^T$

Wherein,Represent random noise ω respectively_k, υ_kCovariance matrix,Represent respectively with Machine noise ω_k, υ_kCharacteristic function andRepresent The characteristic function of estimation difference, t is time variable, and span is (-∞, ∞), and ∞ represents infinite.

Owing to characteristic function the most uniquely determines with probability density function, estimation difference e_k+1Probability density function can lead to Cross following inverse Fourier transform to try to achieve:

With estimation difference e_k+1The probability density function of the Gaussian random variable with identical covariance matrix is: ${\mathrm{\γ}}_{k+1}(x,L_k)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^n\left|c_{k+1}\right|}}{e}^{-\frac{1}{2}{x}^T{c}_{k+1}^{-1}x}$

4th step, solves filter gain:

Owing to target function is nonlinear function, generally can not be by solvingObtain L_k Analytical expression, filter gain is tried to achieve by following gradient algorithm the most here:

L_k=L_k-1+λ_kd_k

Wherein d_kIt is from L_k-1The direction of search set out, i.e.λ_kIt is from L_k-1The edge side set out To d_kThe step-length of the linear search carried out.

Present invention advantage compared with prior art is: present invention utilizes the statistical information of stochastic inputs, selected finger What scalar functions was the most complete characterizes the probabilistic statistical characteristics of estimation difference.By traditional filtering method based on minimum variance It is generalized to, based on minimum negentropy and the filtering method of variance, extend the range of application of Kalman filtering and improve precision.And Asked for the probability density function of estimation difference by the character of characteristic function, reduce amount of calculation.

Accompanying drawing explanation

Fig. 1 is the design flow diagram of a kind of non-gaussian filtering method based on negentropy of the present invention.

Detailed description of the invention

As it is shown in figure 1, the present invention implements following (as a example by a certain movable body the moved along a straight line side of explanation of step Implementing of method):

1, design wave filter

When movable body is for linear motion, following kinetics equation can be obtained:

$\left\{\begin{array}{c}{s}_k=s_{k-1}+{\mathrm{Tv}}_{k-1}+\frac{T^2}{2}{a}_{k-1}\\ v_k=v_{k-1}+{\mathrm{Ta}}_{k-1}\\ a_k=a_{k-1}+{\mathrm{T\ω}}_{k-1}\end{array}\right.$

Position detection value is:

y_k=s_k+υ_k

Wherein, s_k, v_k, a_kBeing respectively movable body at the displacement in k moment, speed, acceleration, T is the sampling period；y_kFor seeing Survey output；ω_kFor acceleration, for the spoorer of movable body, acceleration ω_kIt is random quantity, it is assumed herein that it is in interval [-0.5,0.5] upper obedience is uniformly distributed；υ_kFor observation noise, its probability density function is:

${\mathrm{\γ}}_{{\mathrm{\υ}}_k}\left(x\right)=\left\{\begin{array}{cc}-{6x}^2+3/2& x\∈[-\mathrm{0.5,0.5}]\\ 0& x\∈(-\∞,0.5)\∪(0.5,\∞)\end{array}\right.$

Take state vector $x_k=\left[\begin{array}{c}{s}_k\\ v_k\\ a_k\end{array}\right],$ The kinetics equation that then movable body is for linear motion can be converted into following form:

x_k+1=Ax_k+Gω_k

y_k=Cx_k+υ_k

Wherein:

$A=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right],G\left[\begin{array}{c}0\\ 0\\ T\end{array}\right],C=\left[\begin{array}{ccc}1& 0& 0\end{array}\right],$

ω_k, υ_kFor separate non-gaussian random variables, probability density function is it is known that its characteristic function can simply be asked , it is respectivelyState initial value is taken as x₀ =a₀, then its characteristic function is

For above-mentioned discrete-time linear non-gaussian stochastic system, it is constructed as follows the wave filter of form:

Wherein,For the estimated value of xk, Lk is filter gain to be determined, is tried to achieve by subsequent step 4.Definition is estimated Meter errorThen error equation meets:

$e_{k+1}={\hat{x}}_{k+1}-x_{k+1}=(I-L_{k}C)({\mathrm{Ae}}_k-{\mathrm{G\ω}}_k)+L_k{\mathrm{\υ}}_{k+1}$

Wherein I is 3 × 3-dimensional unit matrix.

2, target function is elected as:

$V_{k+1}\left(L_k\right)=R_{1k}\∫{\mathrm{\γ}}_{e_{k+1}}(x,L_k)\log ({\mathrm{\γ}}_{e_{k+1}}(x,L_k)/{\mathrm{\γ}}_{k+1}(x,L_k))+R_{2k}{c}_{k+1}$

Wherein, c_k+1Represent estimation difference e_k+1Covariance matrix,And γ_k+1X () represents estimation difference respectively e_k+1Probability density function and and e_k+1There is the probability density function of the Gauss distribution of identical covariance matrix；R_1k, R_2kFor Known weight matrix.

3, the asking for of estimation difference covariance matrix and probability density function:

Due to the function that target function is estimation difference probability density function and covariance, it is therefore desirable to first ask for estimating The covariance of error and the expression formula of probability density function.

Due to

$e_x={\hat{x}}_0-x_0$

e₁=(I-L₀C)Ae₀-(I-L₀C)Gω₀+L₀υ₁

$\begin{array}{c}{e}_2=(I-L_{1}C){\mathrm{Ae}}_1-(I-L_{1}C)G{\mathrm{\ω}}_1+L_1{\mathrm{\υ}}_2\\ =(I-L_{1}C)A[(I-L_{0}C){\mathrm{Ae}}_0-(I-L_{0}C){\mathrm{G\ω}}_0+L_0{\mathrm{\υ}}_1]-(I-L_{1}C){\mathrm{G\ω}}_1+L_1{\mathrm{\υ}}_2\\ \·\\ \·\\ \·\end{array}$

e_k+1=(I-L_kC)(Ae_k-Gω_k)+L_kυ_k+1

Can be seen that e_kIt is by e₀；ω₀, ω₁..., ω_k-1；υ₁..., υ_kForm Deng linear combination；Meanwhile, because e₀, ω_i, υ_i+1(i=0,1 ..., k) separate, so e_k, ω_k, υ_k+1Separate.Therefore estimation error covariance function can be by Following formula obtains:

$c(k+1)=(I-L_k{C}_{k+1})A_{k}c\left(k\right)A_k^T{(I-L_k{C}_{k+1})}^T+(I-L_k{C}_{k+1})G_k{c}_{{\mathrm{\ω}}_k}{G}_k^T{(I-L_k{C}_{k+1})}^T+L_k{c}_{{\mathrm{\υ}}_k}{L}_k^T$

Wherein,Represent random noise ω respectively_k, υ_kCovariance matrix.

Owing to characteristic function the most uniquely determines with probability density function, the characteristic function of estimation difference therefore can be utilized Ask for its probability density function.For mutually independent random variables, the characteristic function expression formula of its sum is:

WhereinRepresent random noise ω respectively_k, υ_kCharacteristic function.It is assumed here thatEstimate Initial value

The characteristic function expression formula of known estimation difference, its probability density function can be by following inverse Fourier transform Try to achieve:

With estimation difference e_k+1The probability density function of the Gauss distribution with identical covariance is:

${\mathrm{\γ}}_{k+1}(x,L_k)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^n\left|c_{k+1}\right|}}{e}^{-\frac{1}{2}{x}^T{c}_{k+1}^{-1}x}$

4, filter gain is solved:

Owing to target function is nonlinear function, generally can not be by solvingObtain L_k Analytical expression, filter gain is tried to achieve by following gradient algorithm the most here:

L_k=L_k-1+λ_kd_k

Wherein d_kIt is from L_k-1The direction of search set out, i.e.λ_kIt is from L_k-1The edge side set out To d_kThe step-length of the linear search carried out.

The content not being described in detail in description of the invention belongs to prior art known to professional and technical personnel in the field.

Claims (1)

1. a non-gaussian linear stochaastic system filtering method based on negentropy, it is characterised in that comprise the following steps: first, pin To general discrete-time linear stochastic system, design wave filter；Secondly, target function elects estimation error covariance function as Linear combination with negentropy；Again, utilize the character of characteristic function, ask for the probability density function of estimation difference；Finally, solve Filter gain, minimization target function；Specifically comprise the following steps that

The first step, design wave filter:

Discrete-time linear non-gaussian stochastic system for general:

x_k+1=A_kx_k+G_kω_k

y_k=C_kx_k+υ_k

It is constructed as follows the wave filter of form:

${\hat{x}}_{k+1}=A_k{\hat{x}}_k+L_k(y_{k+1}-C_{k+1}{A}_k{\hat{x}}_k)$

Wherein, x_k∈RⁿFor system mode, y_k∈R^mFor measuring output；A_k∈R^n×nFor known state-transition matrix, G_k∈R^n×sFor Known interference transfer matrix, C_k∈R^m×nFor known output matrix；For system mode x_kEstimated value, L_k∈R^n×mFor wanting The filter gain determined；ω_k∈R^sAnd υ_k∈R^mIt is respectively known to zero-mean, bounded, separate and probability density function Nongausian process noise and measurement noise, its probability density function is respectively Initial vector x₀With ω_k, υ_kSeparate, probability density function isα₁, β₁, α₂, β₂For known real number； R represents real number field；

Definition estimation differenceThen error equation meets:

e_k+1=(I-L_kC_k+1)(A_ke_k-G_kω_k)+L_kυ_k+1

Wherein I is that n × n ties up unit matrix；

Second step, target function elects as:

$V_{k+1}\left(L_k\right)=R_{1k}\∫{\mathrm{\γ}}_{e_{k+1}}(x,L_k)\log ({\mathrm{\γ}}_{e_{k+1}}(x,L_k)/{\mathrm{\γ}}_{k+1}(x,L_k))+R_{2k}{c}_{k+1}$

Wherein, c_k+1Represent estimation difference e_k+1Covariance function,And γ_k+1X () represents estimation difference e respectively_k+1General Rate density function and and e_k+1There is the probability density function of the Gaussian random variable of identical covariance matrix；R_1k, R_2kFor The weight matrix known；

Asking for of 3rd step, estimation difference covariance matrix and probability density function:

The covariance matrix of error equation and characteristic function meet following two formulas respectively:

$c(k+1)=(I-L_k{C}_{k+1})A_{k}c\left(k\right)A_k^T{(I-L_k{C}_{k+1})}^T+(I-L_k{C}_{k+1})G_k{c}_{{\mathrm{\ω}}_k}{G}_k^T{(I-L_k{C}_{k+1})}^T+L_k{c}_{{\mathrm{\υ}}_k}{L}_k^T$

Wherein,Represent random noise ω respectively_k, υ_kCovariance matrix,Represent respectively and make an uproar at random Sound ω_k, υ_kCharacteristic function andRepresent and estimate The characteristic function of error, t is time variable, and span is (-∞, ∞), and ∞ represents infinite；

Owing to characteristic function the most uniquely determines with probability density function, estimation difference e_k+1Probability density function can be by such as Under inverse Fourier transform try to achieve:

With estimation difference e_k+1The probability density function of the Gaussian random variable with identical covariance matrix is:

${\mathrm{\γ}}_{k+1}(x,L_k)=\frac{1}{\sqrt{{\left(2\mathrm{\π}\right)}^n\left|c_{k+1}\right|}}{e}^{-\frac{1}{2}{x}^T{c}_{k+1}^{-1}x}$

4th step, solves filter gain:

Owing to target function is nonlinear function, generally can not be by solvingObtain L_kSolution Analysis expression formula, filter gain is tried to achieve by following gradient algorithm the most here:

L_k=L_k-1+λ_kd_k

Wherein d_kIt is from L_k-1The direction of search set out, i.e.λ_kIt is from L_k-1Set out along direction d_kEnter The step-length of the linear search of row.