CN104022757B - A kind of linear expansion method of the multilamellar Unscented kalman filtering device of High Order Moment coupling - Google Patents

Abstract

The invention discloses the linear expansion method of the Unscented kalman filtering device of a kind of High Order Moment coupling, belong to nonlinear filtering technique field.This method comprises the following steps: 1) sets up the state equation of nonlinear system and measures equation；2) initial state value of system is determined；3) state estimation based on previous step and state equation, use linear expansion Unscented transform calculates the distribution characteristics of the stochastic variable of a step status predication；4) linear expansion Unscented transform is used to calculate the stochastic variable of status predication distribution characteristics after measuring equation transform；5) prediction of Kalman gain Fusion Strain and actual measurement data is used to calculate the distribution characteristics of optimum state；6) judge whether iteration terminates.The present invention fills the thought of use ratio correction and the structure of Symmetric Orthogonal sample and multiple sampling is combined, and by mating more High Order Moment so that close approximation precision significantly improves, reduces computation complexity, greatly improves computational efficiency.

Description

A kind of linear expansion method of the multilamellar Unscented kalman filtering device of High Order Moment coupling

Technical field

The present invention relates to the information fusion technology fields such as nonlinear filtering, Digital Signal Processing, target locating, specifically Relating to is the linear expansion method (Linear-Extension of multilamellar Unscented kalman filtering device of a kind of High Order Moment coupling Unscented Kalman Fliter, LUKF).

Background technology

Almost all of reality system is all nonlinear, especially in aircraft navigation, target following and Industry Control etc. Field.Such as: during target locating, when utilizing radar to be observed aerial target, radar is obtained in that in the air Target is relative to the azimuth of self, but this observation contains noise, and in observational equation, the azimuthal observation amount of radar is target to be estimated The nonlinear function of location parameter, it is impossible to directly utilizing linear filter method and obtain the kinestate of target, it is in the nature non-thread Property filtering problem, is the common difficulty of the research field such as target following, Digital Signal Processing.

For Nonlinear Filtering Problem, frequently with two class filtering methods: a class is that nonlinear function is carried out linearisation is near Seemingly, higher order term is used the measure ignored or approach, is wherein most widely used that extended Kalman filter (Extended Kalman Filter, EKF), its basic ideas are that the Taylor expansion to nonlinear function carries out first-order linear and blocks, Thus nonlinear problem is converted into linearly；Another kind of is to use the distribution of method of sampling approximate non-linear, and the conventional particle that has is filtered Ripple device (Particle Filer, PF) and Unscented kalman filtering device (Unscented Kalman Filter, UKF), it is basic Principle is the distribution that use sample point combines the stochastic variable of its weight Nonlinear Function Approximation.

Comparing with Unscented kalman filtering device, it is not enough that EKF has following three points: (1) is when nonlinear function Taylor exhibition When the higher order term of open type cannot be ignored, it is stable that linearisation can make the error that system generation is bigger, even wave filter be difficult to；(2) In many practical problems, hardly result in the Jacobian matrix of nonlinear function, even do not exist；(3) EKF needs derivation, so It must be appreciated that understand the concrete form of nonlinear function, cannot accomplish that black box encapsulates, thus be difficult to modular applications.At present, though So EKF being had numerous improved methods, as high-order blocks EKF, iteration EKF etc., but these defects are still difficult to overcome.Research table The estimated result that bright UKF is given with accurately, can reach the computational accuracy of higher order than EKF, and amount of calculation and the same order of EFK.With Unscented kalman filtering device compares, the random sample point that particle filter uses, and the quantity of needs is very big, and its sample Point quantity along with the dimension of problem be geometrical progression increase, its calculation cost is sufficiently expensive.The determining that property that UKF method is taked Be distributed closely-related typical sample point with it, its quantity greatly reduces relatively, and general dimension is that the random distribution of n has only to 2n+1 sample point, the most only needs n+1 sample point i.e. to can reach the precision of coupling second moment.

Common, the main flow sampling policy of Unscented kalman filtering device is broadly divided into four kinds: (1) symmetric sampling, its feature Be sample point be symmetrical about average point；(2) simple form sampling, is characterized in that sampling sigma point distribution is not about all Being worth point-symmetric, sample point is few, more more than dimension；(3) ratio correction, is characterized in that comparing the point of sample Example scales；(4) high-order sampling, is characterized in that the High Order Moment using various strategy matching to be distributed.

At present, the UKF that people apply in target tracking domain is second order UKF method, for Gauss nonlinear system, The estimated accuracy of second order UKF, can only achieve three Taylor expansions of nonlinear function, limited precision, and in actual application, urgently Need the sampling policy taking into account precision with computational efficiency.

Summary of the invention

In view of this, it is an object of the invention to provide the line of the multilamellar Unscented kalman filtering device of a kind of High Order Moment coupling Property extended method, the method for solving nonlinear filter precision in actual application and computational efficiency problem, knot Close existing sampling policy, use High Order Moment and the sampling of multiple single file to improve precision, meanwhile, use symmetric sampling and simplification Ratio modification method reduces computation complexity, thus realizes the Synchronous lifting of precision and computational efficiency.

For reaching above-mentioned purpose, the present invention provides following technical scheme:

The linear expansion method of the multilamellar Unscented kalman filtering device of a kind of High Order Moment coupling, specifically comprises the steps of

1) according to practical engineering application, set up the state equation of nonlinear system and measure equation；

2) determine the random distribution characteristic of the initial state value of system, i.e. original state, including its average, covariance and High Order Moment, the distribution characteristics of noise, and initial measurement；

3) a step status predication: state estimation based on previous step and state equation, uses linear expansion Unscented transform meter Calculate the distribution characteristics of the stochastic variable of a step status predication；

4) one step surveying prediction: based on step 3) status predication and measure equation, use linear expansion Unscented transform meter The stochastic variable of calculation status predication distribution characteristics after measuring equation transform；

5) prediction of Kalman gain Fusion Strain and actual measurement data is used to calculate the distribution characteristics of optimum state, complete Become nonlinear system unified step estimation task；

6) judge whether iteration terminates, if do not terminated, then bring the feature of the stochastic variable currently walked into step 3) make For the state estimation of previous step, carry out the next step calculating.

Further, step 1) described in state equation and measure equation be:

Wherein: x_k+1N for kth+1 step ties up state vector, z_k+1M for kth+1 step ties up measurement vector, f () and h () For nonlinear function, w_kFor n tie up stochastic system noise, and system noise obey average be zero, covariance be Q_kGauss distribution； v_k+1For m dimension random measurement noise, wherein measure noise obey average be zero, covariance be R_kGauss distribution, and system Noise and measurement noise are orthogonal；Function f (x_k) be system mode conversion mathematical model, function h (x_k+1) it is system mode The mathematical model measured.

Further, step 3) described state estimation based on previous step and state equation, use linear expansion to become without mark Change the distribution characteristics of stochastic variable calculating a step status predication, specifically include following four step:

3-1) according to status predication stochastic variable and the associating (x of system noise stochastic variable of previous step_k,w_k) point Cloth feature, determines the High Order Moment of required coupling, and then determines layering number of plies l, according to one positive number sequence of particular problem feature selection Row 0 ＜ r₁＜ r₂＜ ... ＜ r_lDetermine the layering of sample point, based on this layering, determine sample point, sigma point choose formula such as Under:

$\left\{\begin{array}{ccc}{\mathrm{\&chi;}}_{0,k}=\left[\begin{array}{c}{\left(x_{0,k}\right)}_{n\&times;1}\\ {\left(W_{0,k}\right)}_{n\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k|k}\right)}_{n\&times;1}\\ 0_{n\&times;1}\end{array}\right]& & \\ {\mathrm{\&chi;}}_{i,k}^j=\left[\begin{array}{c}{\left(x_{i,k}^j\right)}_{n\&times;1}\\ {\left(W_{i,k}^j\right)}_{n\&times;1}\end{array}\right]=\left[\begin{array}{c}{({\tilde{x}}_{k|k}+{\left(\sqrt{r_j{P}_{x,k|k}}\right)}_i)}_{n\&times;1}\\ 0_{n\&times;1}\end{array}\right]& 1\&le;i\&le;n,1\&le;j\&le;l& \\ {\mathrm{\&chi;}}_{i,k}^j=\left[\begin{array}{c}{\left(x_{i,k}^j\right)}_{n\&times;1}\\ {\left(W_{i,k}^j\right)}_{n\&times;1}\end{array}\right]\text{=}\left[\begin{array}{c}{({\tilde{x}}_{k|k}-{\left(\sqrt{r_j{P}_{x,k|k}}\right)}_{i-n})}_{n\&times;1}\\ 0_{n\&times;1}\end{array}\right]& n<i\&le;2n,1\&le;j\&le;l& \\ {\mathrm{\&chi;}}_{i,k}^j=\left[\begin{array}{c}{\left(x_{i,k}^j\right)}_{n\&times;1}\\ {\left(W_{i,k}^j\right)}_{n\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k|k}\right)}_{n\&times;1}\\ {\left({\left(\sqrt{r_j{Q}_k}\right)}_{i-2n}\right)}_{n\&times;1}\end{array}\right]& 2n<i\&le;3n,1\&le;j\&le;l& \\ {\mathrm{\&chi;}}_{i,k}^j=\left[\begin{array}{c}{\left(x_{i,k}^j\right)}_{n\&times;1}\\ {\left(W_{i,k}^j\right)}_{n\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k|k}\right)}_{n\&times;1}\\ {(-{\left(\sqrt{r_j{Q}_k}\right)}_{i-3n})}_{n\&times;1}\end{array}\right]& 3n<i\&le;4n,1\&le;j\&le;l& \end{array}\right.---\left(2\right)$

Wherein:Corresponding stochastic variable (x in corresponding kth step_k,w_k) jth layer i-th sigma point vector, by shape State sample point vectorWith system noise sample point vectorComposition；χ_0,kIt is the sigma point vector of innermost layer,Represent Kth step optimal estimation mean vector, P_x,k|kOptimal estimation covariance is walked for kth；

3-2) give step 3-1) the sigma point that obtained composes weight, and it is regular as follows: χ in formula (2)_0,kCorresponding power It is heavily W_0,k,Corresponding weight isThen require that the weighted value of the sample point of each layer is the biggest, i.e.And 1≤j≤l, 1≤η≤2n, 1≤λ≤2n；

3-3) coupling High Order Moment: solve system of linear equations according to formula (3), try to achieve weight:

$\left\{\begin{array}{c}\left[\begin{array}{c}{\tilde{x}}_{k|k}\\ 0_{n\&times;1}\end{array}\right]=W_{0,k}{\mathrm{\&chi;}}_0+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{\mathrm{\&chi;}}_{i,k}^j\\ P_{\mathrm{xw},k|k}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j({\mathrm{\&chi;}}_{i,k}^j-{\mathrm{\&chi;}}_{0,k}){({\mathrm{\&chi;}}_{i,k}^j-{\mathrm{\&chi;}}_{0,k})}^T\\ m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k|k}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{({\mathrm{\&chi;}}_{i,k}^j-{\mathrm{\&chi;}}_{0,k})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\\ \underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{4n}{\mathrm{\&Sigma;}}}{W}_{i,k}^j=1\end{array}\right.---\left(3\right)$

Wherein, $P_{\mathrm{xw},k|k}=\left[\begin{array}{cc}{P}_{x,k|k}& 0_{n\&times;n}\\ 0_{n\&times;n}& Q_k\end{array}\right],$ For random vector x_k|kThe α rank square of β stochastic variable.

3-4) distribution characteristics of stochastic variable state equation conversion calculates:

According to transforming function transformation function, calculate sigma point conversion sigma point after state equation convertsIts computational methods For:

Y_0,k=f (x_0,k)+W_0,k, $Y_{i,k}^j=f\left(x_{i,k}^j\right)+W_{i,k}^j,i=1,...,4n,1\&le;j\&le;l---\left(4\right)$

Conversion stochastic variable x is calculated according to formula (5)_k+1|k=f (x_k)+w_kMean vector:

${\tilde{x}}_{k+1|k}=W_{0,k}{Y}_{0,k}+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{Y}_{i,k}^j---\left(5\right)$

Conversion stochastic variable x is calculated according to formula (6)_k+1|kCovariance:

$P_{x,k+1|k}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j(Y_{i,k}^j-{\tilde{x}}_{k+1|k}){(Y_{i,k}^j-{\tilde{x}}_{k+1|k})}^T---\left(6\right)$

Conversion stochastic variable x is calculated according to formula (7)_k+1|kHigh Order Moment:

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k+1|k}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{(Y_{i,k}^j-{\tilde{x}}_{k+1|k})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(7\right)$

Further, step 4) described based on step 3) status predication and measure equation, use linear expansion to become without mark Change the stochastic variable calculating status predication through measuring the distribution characteristics after equation transform, specifically include following four step:

4-1) according to status predication stochastic variable and the associating (x of measurement noise stochastic variable_k+1|k,v_k+1) distribution special Levy, determine the High Order Moment of required coupling, and then determine layering number of plies l, according to one sequence of positive numbers 0 ＜ of particular problem feature selection r₁＜ r₂＜ ... ＜ r_lDetermine the layering of sample point, based on this layering, determine sample point, sigma point to choose formula as follows:

$\left\{\begin{array}{ccc}{X}_{0,k+1}=\left[\begin{array}{c}{\left(X_{0,k+1}\right)}_{n\&times;1}\\ {\left(V_{0,k+1}\right)}_{m\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k+1|k}\right)}_{n\&times;1}\\ 0_{m\&times;1}\end{array}\right]& & \\ X_{i,k+1}^j=\left[\begin{array}{c}{\left(X_{i,k+1}^j\right)}_{n\&times;1}\\ {\left(V_{i,k+1}^j\right)}_{m\&times;1}\end{array}\right]=\left[\begin{array}{c}{({\tilde{x}}_{k+1|k}+{\left(\sqrt{r_j{P}_{x,k+1|k}}\right)}_i)}_{n\&times;1}\\ 0_{m\&times;1}\end{array}\right]& 1\&le;i\&le;n,1\&le;j\&le;l& \\ X_{i,k+1}^j=\left[\begin{array}{c}{\left(X_{i,k+1}^j\right)}_{n\&times;1}\\ {\left(V_{i,k+1}^j\right)}_{m\&times;1}\end{array}\right]\text{=}\left[\begin{array}{c}{({\tilde{x}}_{k+1|k}-{\left(\sqrt{r_j{P}_{x,k+1|k}}\right)}_{i-n})}_{n\&times;1}\\ 0_{m\&times;1}\end{array}\right]& n<i\&le;2n,1\&le;j\&le;l& \\ X_{i,k+1}^j=\left[\begin{array}{c}{\left(X_{i,k+1}^j\right)}_{n\&times;1}\\ {\left(V_{i,k+1}^j\right)}_{m\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k+1|k}\right)}_{n\&times;1}\\ {\left({\left(\sqrt{r_j{R}_k}\right)}_{i-2n}\right)}_{m\&times;1}\end{array}\right]& 2n<i\&le;2n+m,1\&le;j\&le;l& \\ X_{i,k+1}^j=\left[\begin{array}{c}{\left(X_{i,k+1}^j\right)}_{n\&times;1}\\ {\left(V_{i,k+1}^j\right)}_{m\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k+1|k}\right)}_{n\&times;1}\\ {(-{\left(\sqrt{r_j{R}_k}\right)}_{i-3n})}_{m\&times;1}\end{array}\right]& 2n+m<i\&le;2(n+m),1\&le;j\&le;l& \end{array}\right.---\left(8\right)$

Wherein,Represent (x in kth step_k+1|k,v_k+1) the i-th sigma point vector of corresponding jth layer, by measurement sample This some vectorVectorial with measuring noise sample pointComposition, X_0,k+1It is the sigma point of corresponding innermost layer,Table Show that the status predication that kth walks estimates mean vector, P_x,k+1|kStatus predication estimate covariance for kth step.

4-2) give step 3-1) the sigma point that obtained composes weight, and it is regular as follows: in formula (8), X_0,k+1Corresponding Weight is W_0,k+1,Corresponding weight isThen require that the weighted value of the sample point of each layer is the biggest, i.e.And 1≤j≤l, 1≤η≤2n, 1≤λ≤2n；

4-3) coupling High Order Moment: solve system of linear equations according to formula (9), try to achieve weight:

$\left\{\begin{array}{c}\left[\begin{array}{c}{\tilde{x}}_{k+1|k}\\ 0_{m\&times;1}\end{array}\right]=W_{0,k+1}{X}_{0,k+1}+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j{X}_{i,k+1}^j\\ P_{\mathrm{xv},k+1|k}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j(X_{i,k+1}^j-X_{0,k+1}){(X_{i,k+1}^j-X_{0,k+1})}^T\\ m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k+1|k}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j{(X_{i,k+1}^j-X_{0,k+1})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\\ \underset{j=1}{\overset{l}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{2(m+n)}{\mathrm{\&Sigma;}}}{W}_{i,k+1}^j=1\end{array}\right.---\left(9\right)$

Wherein, $P_{\mathrm{xv},k+1|k}=\left[\begin{array}{cc}{P}_{x,k+1|k}& 0_{n\&times;m}\\ 0_{m\&times;n}& R_k\end{array}\right],$ For random vector x_k+1|kThe α rank of β stochastic variable Square；

4-4) stochastic variable distribution characteristics after measuring equation transform calculates:

Calculate the conversion sigma point after the sigma measured equation transform of pointIts computational methods are:

Z_0,k+1=h (X_0,k+1)+V_0,k+1, $Z_{i,k+1}^j=h\left(X_{i,k+1}^j\right)+V_{i,k+1}^j,i=1,...,2(m+n),1\&le;j\&le;l---\left(10\right)$

Conversion stochastic variable z is calculated according to formula (11)_k+1=h (x_k+1)+v_k+1Mean vector:

${\tilde{z}}_{k+1|k}=W_{0,k+1}{Z}_{0,k+1}+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j{Z}_{i,k+1}^j---\left(11\right)$

Conversion stochastic variable z is calculated according to formula (12)_k+1|kCovariance:

$P_{z,k+1|k}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j(Z_{i,k+1}^j-{\tilde{z}}_{k+1|k}){(Z_{i,k+1}^j-{\tilde{z}}_{k+1|k})}^T---\left(12\right)$

High Order Moment according to formula (13) calculating conversion stochastic variable:

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(z_{k+1|k}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j{(Z_{i,k+1}^j-{\tilde{z}}_{k+1|k})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(13\right)$

Based on sigma pointWithIntercept about x_kAnd x_k+1|kSamplingWithCalculate according to formula (14) x_k+1|kAnd z_k+1|kCross covariance:

$P_{\mathrm{xz},k+1|k}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j(X_{i,k+1}^j-{\tilde{x}}_{k+1|k}){(Z_{i,k+1}^j-{\tilde{z}}_{k+1|k})}^T---\left(14\right)$

Further, step 5) described in the prediction of use Kalman gain Fusion Strain and measurement data calculate optimum shape The distribution characteristics of state, the detailed process completing nonlinear system unified step estimation task is:

X is calculated according to formula (15)_k+1Average:

${\tilde{x}}_{k+1}={\tilde{x}}_{k+1|k}+K_{k+1}(z_{k+1}-{\tilde{z}}_{k+1|k})---\left(15\right)$

In formula $K_{k+1}=P_{\mathrm{xz},k+1|k}{P}_{z,k+1|k}^{-1},$ For Kalman gain；

X is calculated according to formula (16) and (17)_k+1The covariance P of optimal estimation_x,k+1|k+1And High Order Moment

$P_{x,k+1|k+1}=P_{x,k+1|k}-K_{k+1}{P}_{z,k+1|k}{K}_{k+1}^T---\left(16\right)$

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k+1}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k+1}^j{(X_{i,k+1}^j+K_{k+1}(y_c-(h\left(X_{i,k+1}^j\right)+V_{i,k+1}^j))-{\tilde{z}}_{k+1|k})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(17\right)$

Wherein: y_cIt it is actual measurement data value.

It is an advantage of the current invention that: the present invention makes full use of the geometric properties of random distribution, use the thought of ratio correction The structure of Symmetric Orthogonal sample and multiple sampling is combined, by mating more High Order Moment so that close approximation essence Degree significantly improves, and meanwhile, the method uses symmetric sampling and simplifying rate modification method to reduce computation complexity, greatly Improve computational efficiency.

The further advantage of the present invention, target and feature will be illustrated to a certain extent in the following description, and And to a certain extent, will be apparent to those skilled in the art based on to investigating hereafter, or can To be instructed from the practice of the present invention.

Accompanying drawing explanation

In order to make the object, technical solutions and advantages of the present invention clearer, below in conjunction with accompanying drawing the present invention made into The detailed description of one step, wherein:

Fig. 1 is the method detailed flow chart of the present invention；

The model of growth application that Fig. 2 is the LUKF method provided based on the present invention and one-dimensional non-stationary based on UKF method In absolute error curve chart, wherein, chain-dotted line is the curve of LUKF, and solid line is the curve of UKF method；

The model of growth application that Fig. 3 is the LUKF method provided based on the present invention and one-dimensional non-stationary based on UKF method In mean square error curve chart, wherein, chain-dotted line is the curve of LUKF, and solid line is the curve of UKF method.

Detailed description of the invention

Below with reference to accompanying drawing, the preferred embodiments of the present invention are described in detail；Should be appreciated that preferred embodiment Only for the explanation present invention rather than in order to limit the scope of the invention.

During target tracking, observation station is obtained in that the target bearing information comprising noise, and azimuth information with treat Between the target position information estimated, relation is non-linear, commonly used EKF, second order UKF, second order CKF or high-order CKF nonlinear filtering Method obtains the kinestate of target.But requirement can not be met at the above-mentioned filtering method of the occasion that target positioning requirements is higher. What the present invention provided method possesses higher estimated accuracy than existing method, it is possible to increase the precision of target following.

When carrying out kth stepping and calculating, illustrate with example below, specifically comprise the following steps that

Step 1): according to the model of growth (UNGM) of one-dimensional non-stationary in economics, set up target following system described below The state equation (18) of system and observational equation (19):

$x_k=\frac{x_{k-1}}{2}+\frac{25\&CenterDot;x_{k-1}}{1+x_{k-1}^2}+8\&CenterDot;\cos \left(\frac{6(k-1)}{5}\right)+w_k---\left(18\right)$

$z_k=\frac{x_{k-1}^2}{20}+v_k---\left(19\right)$

Wherein, k=1,2 ..., N, for emulation step number, in this example, N takes 1000 steps；Stochastic system noise vector w_kObey Average is zero, variance Q_k-1Gauss distribution for unit matrix；Random measurement noise vector v_kObeying average is 0, variance R_kFor list The Gauss distribution of bit matrix.

Step 2): determine the random distribution characteristic of the initial state value of system, i.e. original state, including its average, association side Difference and High Order Moment:

Initial state value: x₀=0.1；Initial covariance matrix P₀=1, high-order away from: $m_1^8\left(x_0\right)=105.$

Step 3): state estimation based on previous step and state equation, use linear expansion Unscented transform to calculate a step shape The distribution characteristics of the stochastic variable of state prediction.

Step 3 specifically includes following four steps:

3-1) according to state estimation stochastic variable and the associating (x of noise stochastic variable of kth-1 step_k-1,w_k-1) point Cloth feature, i.e. averageCovariance P_x,k-1|k-1And variance Q of noise_k-1, special, as k=1, P_x,0|0=P₀；According to random distribution characteristic, owing to there are 8 rank High Order Moment, and then determine that the layering number of plies is 4 layers, according to specifically asking Topic one sequence of positive numbers 0 ＜ r of feature selection₁=1 ＜ r₂=2 ＜ r₃=3 ＜ r₄=4 are layered；

Choosing of sigma point is carried out by formula (20):

$\left\{\begin{array}{ccc}{\mathrm{\&chi;}}_{0,k-1}=\left[\begin{array}{c}{\left(x_{0,k-1}\right)}_{1\&times;1}\\ {\left(W_{0,k-1}\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k-1|k-1}\right)}_{1\&times;1}\\ 0_{1\&times;1}\end{array}\right]& & \\ {\mathrm{\&chi;}}_{i,k-1}^j=\left[\begin{array}{c}{\left(x_{i,k-1}^j\right)}_{1\&times;1}\\ {\left(W_{i,k-1}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{({\tilde{x}}_{k-1|k-1}+{\left(\sqrt{r_j{P}_{x,k-1|k-1}}\right)}_i)}_{1\&times;1}\\ 0_{1\&times;1}\end{array}\right]& i=1,1\&le;j\&le;4& \\ {\mathrm{\&chi;}}_{i,k-1}^j=\left[\begin{array}{c}{\left(x_{i,k-1}^j\right)}_{1\&times;1}\\ {\left(W_{i,k-1}^j\right)}_{1\&times;1}\end{array}\right]\text{=}\left[\begin{array}{c}{({\tilde{x}}_{k-1|k-1}-{\left(\sqrt{r_j{P}_{x,k-1|k-1}}\right)}_{i-1})}_{1\&times;1}\\ 0_{1\&times;1}\end{array}\right]& i=2,1\&le;j\&le;4& \\ {\mathrm{\&chi;}}_{i,k-1}^j=\left[\begin{array}{c}{\left(x_{i,k-1}^j\right)}_{1\&times;1}\\ {\left(W_{i,k-1}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k-1|k-1}\right)}_{1\&times;1}\\ {\left({\left(\sqrt{r_j}\right)}_{i-2}\right)}_{1\&times;1}\end{array}\right]& i=3,1\&le;j\&le;4& \\ {\mathrm{\&chi;}}_{i,k-1}^j=\left[\begin{array}{c}{\left(x_{i,k-1}^j\right)}_{1\&times;1}\\ {\left(W_{i,k-1}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\left({\tilde{x}}_{k-1|k-1}\right)}_{1\&times;1}\\ {(-{\left(\sqrt{r_j}\right)}_{i-3})}_{1\&times;1}\end{array}\right]& i=4,1\&le;j\&le;4& \end{array}\right.---\left(20\right)$

Wherein,Represent corresponding stochastic variable (x in kth-1 step_k-1,w_k-1) jth layer i-th sigma point vector, By state sample point vectorWith system noise sample point vectorComposition；χ_0,k-1It it is the sigma point of corresponding innermost layer Vector,Represent the optimal estimation vector of kth-1 step, P_x,k-1|k-1Optimal estimation covariance matrix for kth-1 step；

3-2) according to following rule, to step 3-1) gained sigma point tax weight: χ in formula (20)_0,k-1Corresponding weight For W_0,k-1,Corresponding weight isWherein i=1 ..., 4, j=1 ..., 4, generally we require the sample of each layer The weighted value of point is the biggest, i.e. to any 1≤j≤l, and 1≤η≤2n, 1≤λ≤2n,Amount to 5 non-rights to know Weight.

3-3) coupling High Order Moment: solve system of linear equations according to formula (3), try to achieve weight；Owing to using balanced sample point, public Formula (3) is rewritten as formula (21), altogether 5 systems of linear equations, existence and unique solution:

$\left\{\begin{array}{c}{P}_{\mathrm{xw},k-1|k-1}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k-1}^j({\mathrm{\&chi;}}_{i,k-1}^j-{\mathrm{\&chi;}}_{0,k-1}){({\mathrm{\&chi;}}_{i,k-1}^j-{\mathrm{\&chi;}}_{0,k-1})}^T\\ m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k-1|k-1}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k-1}^j{({\mathrm{\&chi;}}_{i,k-1}^j-{\mathrm{\&chi;}}_{0,k-1})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\\ \underset{j=1}{\overset{4}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{W}_{i,k-1}^j=1\end{array}\right.---\left(21\right)$

Wherein, $P_{\mathrm{xw},k-1|k-1}=\left[\begin{array}{cc}{P}_{x,k-1|k-1}& 0\\ 0& 1\end{array}\right],$ For random vector x_k-1|k-1The α of β stochastic variable Rank square, α=4 in this example, 6,8.

3-4) calculate the distribution characteristics that stochastic variable state equation converts:

According to transforming function transformation function, calculate sigma point conversion sigma point after state equation convertsAccording to formula (22), its computational methods are:

Y_0,k-1=f (x_0,k-1)+W_0,k-1, $Y_{i,k-1}^j=f\left(x_{i,k-1}^j\right)+W_{i,k-1}^j,i=1,...,4,1\&le;j\&le;4---\left(22\right)$

Then, conversion stochastic variable x is calculated according to formula (23)_k|k-1=f (x_k-1)+w_k-1Mean vector:

${\tilde{x}}_{k|k-1}=W_{0,k-1}{Y}_{0,k-1}+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k-1}^j{Y}_{i,k-1}^j---\left(23\right)$

Conversion stochastic variable x is calculated according to formula (24)_k|k-1Covariance:

$P_{x,k|k-1}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k-1}^j(Y_{i,k-1}^j-{\tilde{x}}_{k|k-1}){(Y_{i,k-1}^j-{\tilde{x}}_{k|k-1})}^T---\left(24\right)$

The x of conversion stochastic variable is calculated according to formula (25)_k|k-1High Order Moment:

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k|k-1}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k-1}^j{(Y_{i,k-1}^j-{\tilde{x}}_{k|k-1})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(25\right)$

Step 4): based on step 3) status predication and measure equation, use linear expansion Unscented transform calculate state pre- The stochastic variable surveyed distribution characteristics after measuring equation transform；

Step 4) specifically include following four step:

4-1) according to step 3) status predication stochastic variable and the associating (x of noise stochastic variable_k|k-1,v_k) distribution Feature, i.e. averageCovariance P_x,k|k-1And measure variance R of noise_k.According to random distribution characteristic, owing to there are 8 rank High Order Moment, and then determine that the layering number of plies is 4 layers, according to one sequence of positive numbers 0 ＜ r of particular problem feature selection₁=1 ＜ r₂=2 ＜ r₃=3 ＜ r₄=4 are layered；Based on this layering, carry out choosing of sigma point by formula (26).

$\left\{\begin{array}{ccc}{X}_{0,k}=\left[\begin{array}{c}{\left(X_{0,k}\right)}_{1\&times;1}\\ {\left(V_{0,k}\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\tilde{x}}_{k|k-1}\\ 0\end{array}\right]& & \\ X_{i,k}^j=\left[\begin{array}{c}{\left(X_{i,k}^j\right)}_{1\&times;1}\\ {\left(V_{i,k}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\tilde{x}}_{k|k-1}+{\left(\sqrt{r_j{P}_{x,k|k-1}}\right)}_i\\ 0\end{array}\right]& i=1,1\&le;j\&le;4& \\ X_{i,k}^j=\left[\begin{array}{c}{\left(X_{i,k}^j\right)}_{1\&times;1}\\ {\left(V_{i,k}^j\right)}_{1\&times;1}\end{array}\right]\text{=}\left[\begin{array}{c}{\tilde{x}}_{k|k-1}-{\left(\sqrt{r_j{P}_{x,k|k-1}}\right)}_{i-1}\\ 0\end{array}\right]& i=2,1\&le;j\&le;4& \\ X_{i,k}^j=\left[\begin{array}{c}{\left(X_{i,k}^j\right)}_{1\&times;1}\\ {\left(V_{i,k}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\tilde{x}}_{k|k-1}\\ {\left(\sqrt{r_j}\right)}_{i-2}\end{array}\right]& i=3,1\&le;j\&le;4& \\ X_{i,k}^j=\left[\begin{array}{c}{\left(X_{i,k}^j\right)}_{1\&times;1}\\ {\left(V_{i,k}^j\right)}_{1\&times;1}\end{array}\right]=\left[\begin{array}{c}{\tilde{x}}_{k|k-1}\\ -{\left(\sqrt{r_j}\right)}_{i-3}\end{array}\right]& i\&le;4,1\&le;j\&le;4& \end{array}\right.---\left(26\right)$

Wherein,Represent stochastic variable (x in kth-1 step_k|k-1,v_k) the i-th sigma point of corresponding jth layer is vectorial, By measuring sample point vectorVectorial with measuring noise sample pointComposition；X_0,kIt it is the sigma point of corresponding innermost layer.

4-2) according to following rule, to step one gained sigma point tax weight: X in formula (26)_0,kCorresponding weight For W_0,k,Corresponding weight isWherein i=1 ..., 4, j=1 ..., 4, generally we require the sample point of each layer Weighted value the biggest, i.e. to any 1≤j≤l, 1≤η≤2n, 1≤λ≤2n,Amount to 5 unknown weights.

4-3) coupling High Order Moment: solve system of linear equations according to formula (9), try to achieve weight.Owing to using balanced sample point, public Formula (9) is rewritten as formula (27), altogether 5 systems of linear equations, existence and unique solution:

$\left\{\begin{array}{c}{P}_{\mathrm{xv},k|k-1}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j(X_{i,k}^j-X_{0,k}){(X_{i,k}^j-X_{0,k})}^T\\ m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_{k|k-1}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{(X_{i,k}^j-X_{0,k})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\\ \underset{j=1}{\overset{4}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{W}_{i,k}^j=1\end{array}\right.---\left(27\right)$

Wherein, $P_{\mathrm{xv},k|k-1}=\left[\begin{array}{cc}{P}_{x,k|k-1}& 0\\ 0& 1\end{array}\right],$ For random vector x_k|k-1The α rank square of β stochastic variable.

4-4) calculate the random distribution characteristics becoming measurement amount equation transform:

According to transforming function transformation function, calculate the conversion sigma point after the sigma measured equation transform of pointAccording to formula (28), its computational methods are:

Z_0,k=h (X_0,k)+V_0,k, $Z_{i,k}^j=h\left(X_{i,k}^j\right)+V_{i,k}^j,i=1,...,4,1\&le;j\&le;4---\left(28\right)$

Then, conversion stochastic variable z is calculated according to formula (29)_k=h (x_k)+v_kMean vector:

${\tilde{z}}_{k|k-1}=W_{0,k}{Z}_{0,k}+\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{Z}_{i,k}^j---\left(29\right)$

Conversion stochastic variable z is calculated according to formula (30)_k|k-1Covariance:

$P_{z,k|k-1}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j(Z_{i,k}^j-{\tilde{z}}_{k|k-1}){(Z_{i,k}^j-{\tilde{z}}_{k|k-1})}^T---\left(30\right)$

Conversion stochastic variable z is calculated according to formula (31)_k|k-1High Order Moment:

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(z_{k|k-1}\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{(Z_{i,k}^j-{\tilde{z}}_{k|k-1})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(31\right)$

Based on sigma pointWithIntercept about x_k-1And x_k|k-1SamplingWithWe are according to formula (32) x is calculated_k|k-1And z_k|k-1Cross covariance:

$P_{\mathrm{xz},k|k-1}=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j(X_{i,k}^j-{\tilde{x}}_{k|k-1}){(Z_{i,k}^j-{\tilde{z}}_{k|k-1})}^T---\left(32\right)$

Step 5): use the prediction of Kalman gain Fusion Strain and actual measurement data to calculate the distribution spy of optimum state Levy, complete nonlinear system unified step estimation task, and iteration returns to step 3), carry out subsequent time and estimate task.

Particularly as follows: calculate x according to formula (33)_kAverage,

${\tilde{x}}_k={\tilde{x}}_{k|k-1}+K_k(z_k-{\tilde{z}}_{k|k-1})---\left(33\right)$

In formula, $K_k=P_{\mathrm{xz},k|k-1}{P}_{z,k|k-1}^{-1}$ For Kalman gain；

X is calculated according to formula (16)_kThe covariance P of optimal estimation_x,k|k:

$P_{x,k|k}=P_{x,k|k-1}-K_k{P}_{z,k|k-1}{K}_k^T---\left(34\right)$

X is calculated according to formula (16)_kThe High Order Moment of optimal estimation

$m_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}\left(x_k\right)=\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{W}_{i,k}^j{(X_{i,k}^j+K_k(y_c-(h\left(X_{i,k}^j\right)+V_{i,k}^j))-{\tilde{z}}_{k|k-1})}_{\mathrm{\&beta;}}^{\mathrm{\&alpha;}}---\left(35\right)$

Then judge whether iteration terminates, if do not terminated, then by stochastic variable x_kFeature bring step 3 into as upper The state estimation of one step, carries out the calculating of kth+1 step.

Implementation process: use the absolute error performance indications being defined below in simulation process, is shown in formula (36), mean square error Performance indications, are shown in formula (37), relatively various filtering methods:

$\mathrm{AE}\left(k\right)=|{\tilde{x}}_k-x_k|---\left(36\right)$

$\mathrm{MSE}\left(k\right)=\frac{1}{k}\sqrt{\underset{n=1}{\overset{k}{\mathrm{\&Sigma;}}}{({\tilde{x}}_k-x_k)}^2}---\left(37\right)$

The mean square error estimating target location is the least, shows that target tracking accuracy is the highest, and effect is the best.Wherein, emulation Time is 1000 steps.Under identical simulated conditions, being respectively adopted UKF method and LUKF method is tested, every kind of method is all entered The Monte Carlo emulation that row is 500 times, to compare.

Implementation result: Fig. 2, Fig. 3 sets forth UNGM application in, based on the present invention provide LUKF method with based on Absolute error curve in the UNGM application of UKF method, the LUKF method provided based on the present invention and UNGM based on UKF method Mean square error curve in application.It is clear that LUKF method is relative to the fluctuation maximum of the absolute error curve of UKF method Much smaller, and it is more steady to fluctuate；LUKF method is less relative to the mean square error of UKF method, and all levels off to one stably Error amount.

Finally illustrate, preferred embodiment above only in order to technical scheme to be described and unrestricted, although logical Cross above preferred embodiment the present invention to be described in detail, it is to be understood by those skilled in the art that can be In form and it is made various change, without departing from claims of the present invention limited range in details.

Claims (3)

1. the linear expansion method of the multilamellar Unscented kalman filtering device of a High Order Moment coupling, it is characterised in that specifically comprise Following steps:

1) according to practical engineering application, set up the state equation of nonlinear system and measure equation；

2) random distribution characteristic of the initial state value of system, i.e. original state is determined, including its average, covariance and high-order Square, the distribution characteristics of noise, and initial measurement；

3) a step status predication: state estimation based on previous step and state equation, uses linear expansion Unscented transform to calculate one The distribution characteristics of the stochastic variable of step status predication；

4) one step surveying prediction: based on step 3) status predication and measure equation, use linear expansion Unscented transform calculate shape The stochastic variable of state prediction distribution characteristics after measuring equation transform；

5) use the prediction of Kalman gain Fusion Strain and actual measurement data to calculate the distribution characteristics of optimum state, complete non- Linear system one step estimates task；

6) judge whether iteration terminates, if do not terminated, then bring the feature of the stochastic variable currently walked into step 3) as upper The state estimation of one step, carries out the next step calculating；

Step 1) described in state equation and measure equation be:

Wherein: x_k+1N for kth+1 step ties up state vector, z_k+1M for kth+1 step ties up measurement vector, and f () and h () is non- Linear function, w_kFor n tie up stochastic system noise, and system noise obey average be zero, covariance be Q_kGauss distribution；v_k+1For M dimension random measurement noise, wherein measure noise obey average be zero, covariance be R_kGauss distribution, and system noise Noise is orthogonal with measuring；Function f (x_k) be system mode conversion mathematical model, function h (x_k+1) it is that system mode is measured Mathematical model；

Step 3) described state estimation based on previous step and state equation, use linear expansion Unscented transform to calculate a step state The distribution characteristics of stochastic variable of prediction, specifically includes following four step:

3-1) according to status predication stochastic variable and the associating (x of system noise stochastic variable of previous step_k,w_k) distribution special Levy, determine the High Order Moment of required coupling, and then determine layering number of plies l, according to one sequence of positive numbers 0 ＜ of particular problem feature selection r₁＜ r₂＜ ... ＜ r_lDetermine the layering of sample point, based on this layering, determine sample point, sigma point to choose formula as follows:

Wherein:Corresponding stochastic variable (x in corresponding kth step_k,w_k) jth layer i-th sigma point vector, by state sample Point vectorWith system noise sample point vectorComposition；χ_0,kIt is the sigma point vector of innermost layer,Represent that kth walks Excellent estimation mean vector, P_x,k|kOptimal estimation covariance is walked for kth；

3-2) give step 3-1) the sigma point that obtained composes weight, and it is regular as follows: χ in formula (2)_0,kCorresponding weight is W_0,k,Corresponding weight isThen require that the weighted value of the sample point of each layer is the biggest, i.e.And 1≤j≤ L, 1≤η≤2n, 1≤λ≤2n；

3-3) coupling High Order Moment: solve system of linear equations according to formula (3), try to achieve weight:

Wherein, For random vector x_k|kThe α rank square of β stochastic variable；

3-4) distribution characteristics of stochastic variable state equation conversion calculates:

According to transforming function transformation function, calculate sigma point conversion sigma point after state equation convertsIts computational methods are:

Y_0,k=f (x_0,k)+W_0,k,

Conversion stochastic variable x is calculated according to formula (5)_k+1|k=f (x_k)+w_kMean vector:

Conversion stochastic variable x is calculated according to formula (6)_k+1|kCovariance:

Conversion stochastic variable x is calculated according to formula (7)_k+1|kHigh Order Moment:

The linear expansion method of the multilamellar Unscented kalman filtering device of a kind of High Order Moment the most according to claim 1 coupling, It is characterized in that, step 4) described based on step 3) status predication and measure equation, use linear expansion Unscented transform to calculate The stochastic variable of status predication, through measuring the distribution characteristics after equation transform, specifically includes following four step:

4-1) according to status predication stochastic variable and the associating (x of measurement noise stochastic variable_k+1|k,v_k+1) distribution characteristics, really The High Order Moment of fixed required coupling, and then determine layering number of plies l, according to one sequence of positive numbers 0 ＜ r of particular problem feature selection₁＜ r₂ ＜ ... ＜ r_lDetermine the layering of sample point, based on this layering, determine sample point, sigma point to choose formula as follows:

Wherein,Represent (x in kth step_k+1|k,v_k+1) the i-th sigma point vector of corresponding jth layer, by measurement sample point VectorVectorial with measuring noise sample pointComposition, X_0,k+1It is the sigma point of corresponding innermost layer,Represent kth The status predication of step estimates mean vector, P_x,k+1|kStatus predication estimate covariance for kth step；

4-2) give step 3-1) the sigma point that obtained composes weight, and it is regular as follows: in formula (8), X_0,k+1Corresponding weight For W_0,k+1,Corresponding weight isThen require that the weighted value of the sample point of each layer is the biggest, i.e.And 1≤j≤l, 1≤η≤2n, 1≤λ≤2n；

4-3) coupling High Order Moment: solve system of linear equations according to formula (9), try to achieve weight:

Wherein, For random vector x_k+1|kThe α rank square of β stochastic variable；

4-4) stochastic variable distribution characteristics after measuring equation transform calculates:

Calculate the conversion sigma point after the sigma measured equation transform of pointIts computational methods are:

Z_{0, k+1}=h (X_0,k+1)+V_0,k+1,

Conversion stochastic variable z is calculated according to formula (11)_k+1=h (x_k+1)+v_k+1Mean vector:

Conversion stochastic variable z is calculated according to formula (12)_k+1|kCovariance:

High Order Moment according to formula (13) calculating conversion stochastic variable:

Based on sigma pointWithIntercept about x_kAnd x_k+1|kSamplingWithX is calculated according to formula (14)_k+1|k And z_k+1|kCross covariance:

。

The linear expansion method of the multilamellar Unscented kalman filtering device of a kind of High Order Moment the most according to claim 1 coupling, It is characterized in that, step 5) described in use Kalman gain Fusion Strain prediction and measurement data calculate optimum state point Cloth feature, the detailed process completing nonlinear system unified step estimation task is:

X is calculated according to formula (15)_k+1Average:

In formulaFor Kalman gain；

X is calculated according to formula (16) and (17)_k+1The covariance P of optimal estimation_x,k+1|k+1And High Order Moment

Wherein: y_cIt it is actual measurement data value.