Likelihood-based iteration square-root cubature Kalman filter with applications to state estimation of re-entry ballistic target

Abstract

A new algorithm named the likelihood-based iteration square-root cubature Kalman filter (LISRCKF) is provided in this study. The LISRCKF inherits the virtues of the square-root cubature Kalman filter (SRCKF), which uses the cubature rule-based numerical integration method to calculate the mean and square root of covariance for the non-linear random function. The LISRCKF involves the use of the iterative measurement update and the use of the latest measurement, and the iteration termination criterion based on maximum likelihood is introduced in the measurement update. The LISRCKF algorithm is applied to the state estimation for re-entry ballistic target with unknown ballistic coefficient. Its performance is compared against that of the unscented Kalman filter and SRCKF. Moreover, the suitable choice of iteration number is studied; iteration number 5 is the most appropriate for the LISRCKF algorithm. Simulation results indicate that the LISRCKF algorithm has the features of short run time and fast convergence rate; the advantage in robustness is also demonstrated through the numerical simulation, and it is an effective state estimation method.

Keywords

Cubature Kalman filter maximum likelihood surface non-linear state estimation re-entry ballistic target

Introduction

A series of estimators have been developed to estimate the state of a non-linear system for the last few decades. Up to now, the commonly used state estimator is the extended Kalman filter (EKF) (Bar-Shalom et al., 2001; Grewal and Andrews, 2008). The EKF is based on first-order Taylor approximations of the state transition and observation equation about the estimated state trajectory under Gaussian assumption, so its accuracy depends on the ‘severity’ of the non-linearities. The EKF may introduce significant bias or even convergence problems due to the overly crude approximation (Julier and Uhlmann, 2004). To improve the estimation accuracy, the sophisticated second-order Kalman filter (second-order KF) is available, but it happens at the cost of an increased computational burden (Athans et al., 1968). Another attempt to improve the performance of the EKF, the iterated extended Kalman filter (IEKF), was introduced in Gelb (1974). The IEKF algorithm involves the use of an iterative measurement update. The IEKF is more accurate on the condition that the state estimate is close enough to the true value. However, this is rarely the case in practice (Lefebvre et al., 2004). In reality, the sequence of iterates generated by the IEKF and that generated by the Gauss–Newton method are identical; thus global convergence is guaranteed (Bell and Cathey, 1993). However, the Gauss–Newton method does not ensure that it goes up the likelihood surface (Johnston and Krishnamurthy, 2001; Zhan and Wan, 2007). Moreover, the second-order KF requires Jacobians and Hessians, whereas the EKF requires Jacobians only. Calculation of Jacobians and Hessians often introduces numerical instability and limits their applications.

The derivative-free state estimators have been developed in the past few years. The unscented transformation has been used in the KF framework and the resulting filter is referred to as the unscented Kalman filter (UKF) (Julier et al., 1995, 2000). The basic idea in the UKF is to choose a minimal set of deterministic sample points that capture the mean and covariance of Gaussian density. When propagated through a non-linear function, these points capture the true mean and covariance accurately to the second-order for any non-linearity. However, the UKF may require some fine-tuning in order to prevent the propagation of a non-positive definite covariance matrix for a state vector dimension higher than three. The divided difference filters (DDFs) were introduced in Nørgaard et al. (2000a, b) using Stirling’s interpolation formula. The UKF and DDF algorithms are commonly referred to as sigma point filters due to the properties of deterministic sampling and weighted statistical estimation (Merwe, 2004). Then a new technique called the particle filter (PF) and its variants were introduced in Gordon et al. (1993) and Arulampalam et al. (2002). The basic idea of the PF uses a set of randomly chosen samples with associated weights to approximate the posterior density. In practice, the large number of samples required often makes the use of PF computationally expensive, and the performance of the PF is crucially dependent on the selection of the proposal distribution.

Recently, another suboptimal non-linear filter called the quadrature Kalman filter (QKF) was introduced in Ito and Xiong (2000) and Arasaratnam et al. (2007). The QKF is based on the Gauss–Hermite quadrature rule and calculates the mean and covariance of the states under Gaussian assumption. The QKF introduces the curse of dimensionality to some extent because the number of quadrature points increases exponentially with the state dimension. Lastly, Arasaratnam and Haykin (2009) introduced the cubature Kalman filter (CKF) to ease the curse of dimensionality. The CKF uses the cubature rules for those multi-dimensional integrals to approximate the recursive Bayesian estimation integrals under the Gaussian assumption. The CKF can be applied to solve high-dimensional non-linear filtering problems with minimal computational effort. The computational complexity scales linearity with the state dimension using the third-degree cubature rule. However, the error covariance evaluated in the CKF is required to be decomposed for the Cholesky factors to calculate the cubature points. To mitigate ill effects that may eventually lead to an unstable or even divergent behaviour, the square-root cubature Kalman filter (SRCKF) was developed. The SRCKF propagates the square-root factors of the predictive and posterior error covariance to avoid matrix square-rooting operations and to preserve the positive semi-definiteness of the error covariance and doubled-order precision (Kaminski et al., 1971).

The SRCKF method has proven to be far superior and has more numerical stability. However, because of the large initial error and measurement error in the application for the state estimation of re-entry ballistic targets (RBTs) with unknown ballistic coefficient, the SRCKF also shows its weakness in the robustness and state estimation accuracy. Emboldened by the basic idea of the iterate method in the IEKF as well as the superiority of the SRCKF algorithm, we develop a new filter named the likelihood-based iteration square-root cubature Kalman filter (LISRCKF). The accuracy of the state estimation in the new filter is greatly improved when involving the use of the iterative measurement update and the use of the iterative termination condition based on the maximum likelihood.

The rest of the paper is organized as follows. We begin with a description of the SRCKF. Then, the LISRCKF is developed, after which we apply the LISRCKF algorithm to state estimation for RBT with unknown ballistic coefficient and discuss the simulation results. Finally, the conclusions are given.

Square-root cubature Kalman filter

Consider the state estimation problem of a non-linear dynamics system with additive noise, the $n_{x}$ -dimensional state vector $x_{k}$ of the system evolves according to the non-linear stochastic difference equation and the measurement equation is given by

x_{k} = f (x_{k - 1}) + w_{k - 1}

(1)

z_{k} = h (x_{k}) + v_{k}

(2)

where $f$ : $ℝ^{n_{x}} \to ℝ^{n_{x}}$ and $h$ : $ℝ^{n_{x}} \to ℝ^{n_{z}}$ are some known non-linear functions; $z_{k} \in ℝ^{n_{z}}$ is the measurement; $w_{k - 1}$ and $v_{k}$ are process and measurement Gaussian noise sequences with zero means and covariance $Q_{k - 1}$ and $R_{k}$ , respectively, and ${w_{k - 1}}$ and ${v_{k}}$ are mutually uncorrelated.

The heart of the SRCKF is a three-order spherical-radial cubature rule, which makes it possible to numerically compute multivariate moment integrals encountered in the non-linear Bayesian filter. The SRCKF based on the Equations (1) and (2) can be described as follows (Arasaratnam and Haykin, 2009).

Step 1. Evaluate the basic cubature points and weights. The basic cubature points and corresponding weights are obtained using the three-order cubature rule

ξ_{i} = \sqrt{m / 2} {[1]}_{i}, ω_{i} = 1 / m, i = 1, \dots m

(3)

where $m = 2 n_{x}$ , the ${[1]}_{i}$ is a $n_{x}$ -dimensional vector and is generated in the following way. We denote the $n_{x}$ dimensional unit vector as $e = {[1, 0, \dots 0]}^{T}$ and call it a generator. We suppress the $(n_{x} - 1)$ zero co-ordinates for brevity and use the notation $[1]$ to represent a complete fully symmetric set of points that can be obtained by permutating and changing the sign of the generator $e$ in all possible ways. The ${[1]}_{i}$ is denoted the ith point from the set $[1]$ .

Step 2. Time update. Suppose that the state distribution is $x_{k - 1} ~ N ({\hat{x}}_{k - 1}, P_{k - 1})$ , and a square root $S_{k - 1}$ of $P_{k - 1}$ such that $P_{k - 1} = S_{k - 1} S_{k - 1}^{T}$ is obtained. The cubature points at $k - 1$ instant time are generated and propagated through the non-linear state equation; the predicted state and the square root of the predicted covariance are evaluated as follows.

Calculate the cubature points

X_{i, k - 1} = S_{k - 1} ξ_{i} + {\hat{x}}_{k - 1}

(4)

Propagate the cubature points through the state equation

X_{i, k}^{*} = f (X_{i, k - 1})

(5)

Evaluate the predicted state and square root of the predicted covariance

{\bar{x}}_{k} = \sum_{i = 1}^{m} ω_{i} X_{i, k}^{*}

(6)

{\bar{S}}_{k} = Tria ([χ_{k}^{*} S_{Q, k - 1}])

(7)

where $S_{Q, k - 1}$ denotes a square-root factor of $Q_{k - 1}$ such that $Q_{k - 1} = S_{Q, k - 1} S_{Q, k - 1}^{T}$ and Tria() is denoted a general triagularization algorithm. The weighted, centred (prior mean is subtracted off) matrix $χ_{k}^{*}$ is defined as

χ_{k}^{*} = 1 / 1 \sqrt{m} [X_{1, k}^{*} - {\bar{x}}_{k} X_{2, k}^{*} - {\bar{x}}_{k}, \dots, X_{m, k}^{*} - {\bar{x}}_{k}]

(8)

Step 3. Measurement update.

Calculate the predicted cubature points

X_{i, k} = {\bar{S}}_{k} ξ_{i} + {\bar{x}}_{k}

(9)

Evaluate the propagated cubature points

Z_{i, k} = h (X_{i, k})

(10)

Evaluate the predicted measurement, square root of the innovation covariance and cross-covariance

{\bar{z}}_{k} = \sum_{i = 1}^{m} ω_{i} Z_{i, k}

(11)

S_{zz, k} = Tria ([ϒ_{k} S_{R, k}])

(12)

P_{xz, k | k - 1} = χ_{k} ϒ_{k}^{T}

(13)

Here, $S_{R, k}$ denotes a square-root factor of $R_{k}$ such that $R_{k} = S_{R, k} S_{R, k}^{T}$ and the weighted, centred matrices $ϒ_{k}$ and $χ_{k}$ are defined as

ϒ_{k} = 1 / 1 \sqrt{m} [Z_{1, k} - {\bar{z}}_{k} Z_{2, k} - {\bar{z}}_{k}, \dots, Z_{m, k} - {\bar{z}}_{k}]

(14)

χ_{k} = 1 / 1 \sqrt{m} [X_{1, k} - {\bar{x}}_{k} X_{2, k} - {\bar{x}}_{k}, \dots, X_{m, k} - {\bar{x}}_{k}]

(15)

Evaluate the state estimation and the square root of the covariance at $k$ instant time.

W_{k} = (P_{xz, k} / S_{zz, k}^{T}) / S_{zz, k}

(16)

{\hat{x}}_{k} = {\bar{x}}_{k} + W_{k} (z_{k} - {\bar{z}}_{k})

(17)

S_{k} = Tria ([χ_{k} - W_{k} ϒ_{k} W_{k} S_{R, k}])

(18)

where the symbol ‘/’ represents the matrix right divide operator.

The above Steps 2 and 3 provide a general summary of the SRCKF algorithm. Given the initial condition ${\hat{x}}_{0}$ and $S_{0}$ , the state estimation procedure can be recursively implemented.

Development of likelihood-based iteration square-root cubature Kalman filter

Cubature rule based iterative measurement update

Assuming the predicted state ${\bar{x}}_{k}$ and the current measurement $z_{k}$ at time instant $k$ as realization of independent random vectors with multivariate normal distributions, e.g. ${\bar{x}}_{k} ~ N (x_{k}, {\bar{P}}_{k})$ and $z_{k} ~ N (h (x_{k}), R_{k})$ , the two vectors are formed to a single augmented one $Z = {[z_{k} {\bar{x}}_{k}]}^{T}$ . According to the independent assumption, so we have

Z ~ N (g (x_{k}), \tilde{Q})

(19)

where $g (x_{k}) = {[h (x_{k}), x_{k}]}^{T}$ , $\tilde{Q} = diag ([R_{k} {\bar{P}}_{k}])$ . The update measurement process becomes the one that computing the optimal state estimation and corresponding covariance given $Z$ , $g$ and $\tilde{Q}$ .

Defining the objective function as:

f (ξ) = 0.5 ‖ r (ξ) ‖^{2}

(20)

where $r (ξ) = S (Z - g (ξ))$ , which is the second-order differentiable, and $S^{T} S = {\tilde{Q}}^{- 1}$ . Then Equation (20) is a typical non-linear least squares problem with the initial ${\bar{x}}_{k}$ . Newton’s solution to the non-linear least squares problem is the iterative sequence

{\hat{x}}_{k}^{(i + 1)} = {\hat{x}}_{k}^{(i)} - {(\nabla^{2} r ({\hat{x}}_{k}^{(i)}))}^{- 1} \nabla r ({\hat{x}}_{k}^{(i)})

(21)

For small residual problems, a useful approximation to Equation (21) is the Gauss–Newton method, which approximates the Hessian $\nabla^{2} r ({\hat{x}}_{k}^{(i)})$ by $r' {({\hat{x}}_{k}^{(i)})}^{T} r' ({\hat{x}}_{k}^{(i)})$ , and the gradient of $r (ξ)$ is $\nabla r ({\hat{x}}_{k}^{(i)}) = r' {({\hat{x}}_{k}^{(i)})}^{T} r ({\hat{x}}_{k}^{(i)})$ . Hence, the Gauss–Newton method defines the sequence of iterates

{\hat{x}}_{k}^{(i + 1)} = {\hat{x}}_{k}^{(i)} - {(r' {({\hat{x}}_{k}^{(i)})}^{T} r' ({\hat{x}}_{k}^{(i)}))}^{- 1} r' {({\hat{x}}_{k}^{(i)})}^{T} r ({\hat{x}}_{k}^{(i)})

(22)

where $r' ({\hat{x}}_{k}^{(i)}) = - {SG}_{k}^{(i)}$ is the Jacobian of $r (ξ)$ on ${\hat{x}}_{k}^{(i)}$ , and $G_{k}^{(i)} = \partial g ({\hat{x}}_{k}^{(i)}) / \partial {\hat{x}}_{k}^{(i)} = {[H_{k}^{(i)}, I]}^{T}$ , $H_{k}^{(i)}$ is the Jacobian of $h (x_{k})$ about ${\hat{x}}_{k}^{(i)}$ . Hence, Equation (22) becomes

{\hat{x}}_{k}^{(i + 1)} = {(G_{k}^{(i) T} {\tilde{Q}}^{- 1} G_{k}^{(i)})}^{- 1} G_{k}^{(i) T} {\tilde{Q}}^{- 1} (Z - g (x_{k}^{(i)}) + G_{k}^{(i)} x_{k}^{(i)})

(23)

Noting that ${(G_{k}^{(i) T} {\tilde{Q}}^{- 1} G_{k}^{(i)})}^{- 1} = {(H_{k}^{(i) T} R_{k}^{- 1} H_{k}^{(i)} + {(P_{k}^{(i)})}^{- 1})}^{- 1}$ , we need two forms of matrix inversion lemma 1):

{(H^{T} R^{- 1} H + P^{- 1})}^{- 1} H^{T} R^{- 1} = P H^{T} {(HP H^{T} + R)}^{- 1}

(24)

and Sherman–Morrison–Woodbury matrix inversion lemma 2)

{(H^{T} R^{- 1} H + P^{- 1})}^{- 1} = P - P H^{T} {(HP H^{T} + R)}^{- 1} HP

(25)

where P is an $n \times n$ matrix, R is an $m \times m$ matrix, and $H$ is an $m \times n$ matrix.

Expanding the Equation (23) and making use of above matrix inversion lemma, we obtain

{\hat{x}}_{k}^{(i + 1)} = {\bar{x}}_{k} + {\bar{P}}_{k}^{(i)} H_{k}^{(i) T} {(H_{i}^{T} {({\bar{P}}_{k}^{(i)})}^{- 1} H_{k}^{(i)} + R_{k})}^{- 1} (z_{k} - h ({\hat{x}}_{k}^{(i)}) + H_{k}^{(i)} ({\bar{x}}_{k} - {\hat{x}}_{k}^{(i)})) = {\bar{x}}_{k} + K_{k}^{(i)} (z_{k} - h ({\hat{x}}_{k}^{(i)}) + H_{k}^{(i)} ({\bar{x}}_{k} - {\hat{x}}_{k}^{(i)}))

(25)

where $K_{k}^{(i)}$ is the gain, which is defined as $K_{k}^{(i)} = {\bar{P}}_{k} H_{k}^{(i) T} {(H_{k}^{(i) T} {\bar{P}}_{k}^{(i)} H_{k}^{(i)} + R_{k})}^{- 1}$ .

Equation (26) includes the predicted state ${\bar{x}}_{k}$ , which has a direct and large impact on the final state estimation.

When the measurement model fully observes the state, i.e. if the measurement sensitivity matrix is invertible at any time, then the system state can be uniquely determined. Considering the noiseless measurement equation, $z_{k} = h (x_{k})$ , the first-order Taylor expansion can be given as $z_{k} = h ({\hat{x}}_{k}) + H_{k} (x_{k} - {\hat{x}}_{k})$ . When measurement sensitivity matrix $H_{k}$ is invertible at any time, i.e. $H_{k}^{- 1}$ exists, then the system state $x_{k}$ can be uniquely determined, linearization errors will be small in the uncertainty region around ${\hat{x}}_{k}^{(i)}$ , and the estimated state ${\hat{x}}_{k}^{(i)}$ is more approximate to the true state than the predicted state ${\bar{x}}_{k}$ (Lefebvre et al., 2004).

Substituting ${\bar{x}}_{k}$ by ${\hat{x}}_{k}^{(i)}$ into Equation (26) in order to decrease the influence of ${\bar{x}}_{k}$ on the final state estimation, the following iterative formula is obtained

{\hat{x}}_{k}^{(i + 1)} = {\hat{x}}_{k}^{(i)} + K_{k}^{(i)} (z_{k} - h ({\hat{x}}_{k}^{(i)}))

(27)

Compared with Equation (26), Equation (27) is simpler, noting that Equations (26) and (27) are identical when one iteration is run.

Now we consider the iterative gain $K_{k}^{(i)} = {\bar{P}}_{k} H_{k}^{(i) T} {(H_{k}^{(i) T} {\bar{P}}_{k}^{(i)} H_{k}^{(i)} + R_{k})}^{- 1}$ ; the terms $H_{k}^{(i)} {\bar{P}}_{k} {(H_{k}^{(i)})}^{T} + R_{k}$ and ${\bar{P}}_{k} {(H_{k}^{(i)})}^{T}$ are approximate innovation covariance and cross-covariance, and are expressed as

P_{z z, k}^{(i)} = H_{k}^{(i)} {\bar{P}}_{k} {(H_{k}^{(i)})}^{T} + R_{k}

(28)

P_{xz, k}^{(i)} = {\bar{P}}_{k} {(H_{k}^{(i)})}^{T}

(29)

Equations (28) and (29) are achieved by linearizing the measurement equation (2) with a first-order Taylor expansion so that the linearized error is introduced due to high-order terms being truncated. As for the higher non-linear measurement equation, the accuracy for state estimation is degraded if the linearized error is only propagated for calculating Equations (28) and (29). To decrease the propagated error, we can recalculate equations (9)–(13) to obtain the terms $P_{zz}^{(i)}$ and $P_{xz}^{(i)}$ based on the cubature rule in the following way

X_{j, k}^{(i)} = S_{k}^{(i)} ξ_{j} + {\hat{x}}_{k}^{(i)}

(30)

Z_{j, k}^{(i)} = h (X_{j, k}^{(i)})

(31)

{\hat{z}}_{k}^{(i)} = \sum_{j = 1}^{m} ω_{j} Z_{j, k}^{(i)}

(32)

S_{z z, k}^{(i)} = T r i a ([Υ_{k}^{(i)} S_{R, k}])

(33)

P_{xz, k}^{(i)} = χ_{k}^{(i)} {(ϒ_{k}^{(i)})}^{T}

(34)

P_{zz, k}^{(i)} = S_{zz, k}^{(i)} {(S_{zz, k}^{(i)})}^{T}

(35)

where $ϒ_{k}^{(i)}$ and $χ_{k}^{(i)}$ are defined as

ϒ_{k}^{(i)} = 1 / 1 \sqrt{m} [Z_{1, k}^{(i)} - {\hat{z}}_{k}^{(i)} Z_{2, k}^{(i)} - {\hat{z}}_{k}^{(i)}, \dots, Z_{m, k}^{(i)} - {\hat{z}}_{k}^{(i)}]

(36)

χ_{k}^{(i)} = 1 / 1 \sqrt{m} [X_{1, k}^{(i)} - {\hat{x}}_{k}^{(i)} X_{2, k}^{(i)} - {\hat{x}}_{k}^{(i)}, \dots, X_{m, k}^{(i)} - {\hat{x}}_{k}^{(i)}]

(37)

So $K_{k}^{(i)}$ is expressed as

K_{k}^{(i)} = P_{xz}^{(i)} {(P_{zz}^{(i)})}^{- 1}

(38)

Hence, we can obtain the following iterative formula

{\hat{x}}_{k}^{(i + 1)} = x_{k}^{(i)} + P_{xz}^{(i)} {(P_{zz}^{(i)})}^{- 1} [z_{k} - h ({\hat{x}}_{k}^{(i)})]

(39)

Maximum likelihood-based iteration termination criteria

The commonly used termination criterion in the measurement update of the IEKF algorithm is that the inequality $‖ x_{k}^{(i + 1)} - x_{k}^{(i)} ‖ \leq ε$ is met, where $ε$ is the predetermined threshold. The sequence of iterates generated according to the above termination condition has the property of the global convergence, but it is not guaranteed to go up the likelihood surface and the threshold $ε$ is crucial to using the iterative filtering algorithm successfully; however, selecting a proper choice of $ε$ is difficult. So the iteration termination criterion based on the maximum likelihood surface is introduced.

Under the Gaussian and independent assumption, we consider $x_{k}^{(i)}$ and $z_{k}$ as the realization of independent random vectors with multivariate normal distributions, i.e. $x_{k} ~ N (x_{k}^{(i)}, P_{k}^{(i)})$ and $z_{k} ~ N (h (x_{k}), R_{k})$ . The likelihood function of the two vectors $x_{k}$ and $z_{k}$ is defined as

Λ (x_{k}, z_{k}) = const • \exp {- 0.5 [{(x_{k} - x_{k}^{(i)})}^{T} {(P_{k}^{(i)})}^{- 1} (x_{k} - x_{k}^{(i)}) + {(z_{k} - h (x_{k}))}^{T} R_{k}^{- 1} (z_{k} - h (x_{k}))]}

(40)

Defining the likelihood surface function (LSF) as

\begin{matrix} J (x_{k}) = {(x_{k} - x_{k}^{(i)})}^{T} {(P_{k}^{(i)})}^{- 1} (x_{k} - x_{k}^{(i)}) \\ + {(z_{k} - h (x_{k}))}^{T} R_{k}^{- 1} (z_{k} - h (x_{k})) \end{matrix}

(41)

As we know, obtaining the maximum likelihood estimate is equivalent to solving the minimum value of $J (x_{k})$ . The optimal value is difficult to obtain, but we know that the following inequality holds

J ({\hat{x}}_{k}^{(i + 1)}) < J ({\hat{x}}_{k}^{(i)})

(42)

We say that $J ({\hat{x}}_{k}^{(i + 1)})$ is closer to the maximum likelihood surface than $J ({\hat{x}}_{k}^{(i)})$ , and equivalently, ${\hat{x}}_{k}^{(i + 1)}$ has a more accurate approximation than ${\hat{x}}_{k}^{(i)}$ to the optimal solution. According to Equation (42), we have the following inequality

\begin{matrix} {[{\hat{x}}_{k}^{(i + 1)} - {\hat{x}}_{k}^{(i)}]}^{T} {(P_{k}^{(i)})}^{- 1} [{\hat{x}}_{k}^{(i + 1)} - {\hat{x}}_{k}^{(i)}] \\ + {[z_{k} - h ({\hat{x}}_{k}^{(i + 1)})]}^{T} R_{k}^{- 1} [z_{k} - h ({\hat{x}}_{k}^{(i + 1)})] \\ < {[{\hat{x}}_{k}^{(i)} - {\hat{x}}_{k}^{(i)}]}^{T} {(P_{k}^{(i)})}^{- 1} [{\hat{x}}_{k}^{(i)} - {\hat{x}}_{k}^{(i)}] \\ + {[z_{k} - h ({\hat{x}}_{k}^{(i)})]}^{T} R_{k}^{- 1} [z_{k} - h ({\hat{x}}_{k}^{(i)})] \end{matrix}

(43)

Simplifying Equation (43) and using $P_{k}^{(i)} = S_{k}^{(i)} S_{k}^{(i) T}$ , we immediately generate the following inequality

\begin{matrix} {({\tilde{x}}_{k}^{(i + 1)})}^{T} {(S_{k}^{(i)} S_{k}^{(i) T})}^{- 1} {\tilde{x}}_{k}^{(i + 1)} + {({\tilde{z}}_{k}^{(i + 1)})}^{T} R_{k}^{- 1} {\tilde{z}}_{k}^{(i + 1)} \\ < {({\tilde{z}}_{k}^{(i)})}^{T} R_{k}^{- 1} {\tilde{z}}_{k}^{(i)} \end{matrix}

(44)

here

{\tilde{x}}_{k}^{(i + 1)} = {\hat{x}}_{k}^{(i + 1)} - {\hat{x}}_{k}^{(i)}

(45)

{\tilde{z}}_{k}^{(i + 1)} = z_{k} - h ({\hat{x}}_{k}^{(i + 1)})

(46)

Using the inequality (44) as the criterion to iteration termination, the iteration sequence generated is guaranteed to go up the likelihood surface function.

Likelihood-based iteration square-root cubature Kalman filter

Now we have arrived at the central issue of this paper, namely the LISRCKF. Enlightened by the superiority of the SRCKF algorithm and using the Gauss–Newton iterative method, we derive the LISRCKF. The LISRCKF algorithm involves the use of the iterative measurement update and the use of the current measurement. Moreover, the iteration termination criterion, which guarantees the sequence obtained moves up the maximum likelihood surface, is used in the LISRCKF algorithm.

The procedure of the LISRCKF algorithm is described as follows.

Step 1. Calculate the basic cubature points and corresponding weights according to Equation (3).

Step 2. Time update: evaluate the predicted state ${\bar{x}}_{k}$ and square root ${\bar{S}}_{k}$ of the predicted covariance using (6) and (7).

Step 3. Measurement update:

Let ${\hat{x}}_{k}^{(0)} = {\bar{x}}_{k}$ and $S_{k}^{(0)} = {\bar{S}}_{k}$ . Suppose that the ith iterations are ${\hat{x}}_{k}^{(i)}$ and $S_{k}^{(i)}$ , respectively.

Use Equations (33)–(34) to calculate square root of innovation covariance $S_{zz, k}^{(i)}$ and the cross-covariance $P_{xz, k}^{(i)}$ . Evaluate the estimated state and square root of the corresponding covariance at the k time instant using the following formulas

K_{k}^{(i)} = (P_{xz, k}^{(i)} / {(S_{zz, k}^{(i)})}^{T}) / S_{zz, k}^{(i)}

(47)

{\hat{x}}_{k}^{(i + 1)} = {\hat{x}}_{k}^{(i)} + K_{k}^{(i)} [z_{k} - h ({\hat{x}}_{k}^{(i)})]

(48)

S_{k}^{(i + 1)} = Tria ([χ_{k}^{(i)} - K_{k}^{(i)} ϒ_{k}^{(i)} K_{k}^{(i)} S_{R, k}])

(49)

where $ϒ_{k}^{(i)}$ and $χ_{k}^{(i)}$ are defined as (36) and (37).

If the following inequality holds

\begin{matrix} {({\tilde{x}}_{k}^{(i + 1)})}^{T} {(S_{k}^{(i)} S_{k}^{(i) T})}^{- 1} {\tilde{x}}_{k}^{(i + 1)} + {({\tilde{z}}_{k}^{(i + 1)})}^{T} R_{k}^{- 1} {\tilde{z}}_{k}^{(i + 1)} \\ < {({\tilde{z}}_{k}^{(i)})}^{T} R_{k}^{- 1} {\tilde{z}}_{k}^{(i)} \end{matrix}

(50)

the iterate return to 2); otherwise continue to 4). ${\tilde{x}}_{k}^{(i + 1)}$ and ${\tilde{z}}_{k}^{(i + 1)}$ are defined as (45) and (46).

If the inequality is not satisfied or if i is too large ( $i > N_{max}$ , $N_{max}$ is maximum iteration number predetermined). The iteration number is assumed to be N when the iteration is terminated; the ultimate state estimation and square root of corresponding covariance at the k time instant are

{\hat{x}}_{k} = {\hat{x}}_{k}^{(N)}

(51)

S_{k} = S_{k}^{(N)}

(52)

The LISRCKF algorithm inherits the virtues of the SRCKF algorithm, which is simply implemented and free derivative, and has better numerical stability. In the LISRCKF algorithm, the covariance and cross-covariance based on cubature rule are refined to make the propagated error decrease; the iteration termination condition that makes the sequence move up the maximum likelihood surface is used in the process of measurement update of the LISRCKF algorithm. Hence, the performance of LISRCKF is superior to SRCKF algorithms. Moreover, the measurement update of the LISRCKF algorithm is transformed into the non-linear least-square problem; the optimum state estimation and covariance are solved using Gauss–Newton method, and the sequences obtained have global convergence (Bell and Cathey, 1993). So the LISRCKF algorithm has the same global convergence as the Gauss–Newton method.

Applications to state estimation for re-entry ballistic target

To demonstrate the performance of the LISRCKF algorithm, we apply the LISRCKF to the state estimation for the RBT with unknown ballistic coefficient and compare its performance against the other algorithms.

Model of target motion and of radar measurements

The relative locations of the RBT and the radar are depicted in Figure 1. Two co-ordinate systems are shown in Figure 1.

One orthogonal co-ordinate reference system is the Earth-centred inertial co-ordinate system (ECI-CS) $O x_{I} y_{I} z_{I}$ , which is a right-handed system with the origin O at the Earth centre, axis $O x_{I}$ pointing in the vernal equinox direction and axis $O z_{I}$ pointing in the direction of the North Pole N. Its fundamental plane $O x_{I} y_{I}$ coincides with the Earth’s equatorial plane.

Assume that the radar is situated at the surface of the Earth. Another orthogonal co-ordinate reference system named the East-North-Up co-ordinates system (ENU-CS) $O_{s} xyz$ has its origin at the location of the radar, z is directed along the local vertical and x and y lie in a local horizontal plane, with x pointing east and y pointing north.

The kinematics of the RBT with unknown ballistic coefficient is derived under the two hypotheses. One assumption is that the Earth is spherical and non-rotating. Another is that the forces acting on the target are gravity and drag. The effects of centrifugal acceleration, Coriolis acceleration etc. are ignored, due to their small effect on the trajectory (Li and Jilkov, 2001).

Under the above two assumptions, the kinematics of the RBT with unknown ballistic coefficient in the ENU-CS is described by the following discrete-time non-linear dynamic state equation according to the transformation relationship of ECI-CS and ENU-CS.

x_{k} = Φ x_{k - 1} + G ψ (x_{k - 1}) + w_{k - 1}

(53)

where $x_{k} = {[x_{k} {\overset{\cdot}{x}}_{k} y_{k} {\overset{\cdot}{y}}_{k} z_{k} {\overset{\cdot}{z}}_{k} β_{k}]}^{T}$ is the state of the RBT, $β_{k}$ ( $kg / kg m^{2} m^{2}$ ) is the unknown ballistic coefficient as an estimated state, and

\begin{matrix} Φ = [\begin{matrix} ϕ & 0 & 0 & 0 \\ 0 & ϕ & 0 & 0 \\ 0 & 0 & ϕ & 0 \\ 0 & 0 & 0 & 1 \end{matrix}], ϕ = [\begin{matrix} 1 & T \\ 0 & 1 \end{matrix}], \\ G = [\begin{matrix} τ & 0 & 0 \\ 0 & τ & 0 \\ 0 & 0 & τ \\ 0 & 0 & 0 \end{matrix}], τ = [\begin{matrix} T^{2} / 2 \\ T \end{matrix}] \end{matrix}

(54)

ψ (x_{k - 1}) = [\begin{matrix} - \frac{ρ (h_{k - 1})}{2 β_{k - 1}} V_{k - 1} {\overset{\cdot}{x}}_{k - 1} - \frac{μ x_{k - 1}}{r_{k - 1}^{3}} \\ - \frac{ρ (h_{k - 1})}{2 β_{k - 1}} V_{k - 1} {\overset{\cdot}{y}}_{k - 1} - \frac{μ y_{k - 1}}{r_{k - 1}^{3}} \\ - \frac{ρ (h_{k - 1})}{2 β_{k - 1}} V_{k - 1} {\overset{\cdot}{z}}_{k - 1} - \frac{μ (z_{k - 1} + R_{e})}{r_{k - 1}^{3}} \end{matrix}]

(55)

where

{\begin{matrix} r_{k - 1} = \sqrt{x_{k - 1}^{2} + y_{k - 1}^{2} + {(z_{k - 1} + R_{e})}^{2}} \\ V_{k - 1} = \sqrt{{\overset{\cdot}{x}}_{k - 1}^{2} + {\overset{\cdot}{y}}_{k - 1}^{2} + {\overset{\cdot}{z}}_{k - 1}^{2}} \\ h_{k - 1} = r_{k - 1} - R_{e} \end{matrix}

(56)

T (in s) is the time interval between the radar measurements, µ (=3.986005×10¹⁴ m³/s²) and R_e (=6371004 m) are Earth’s gravitational constant and Earth radius, respectively. $ρ (h)$ (kg/m³) is the air density. Below 90 km at height, it can be approximately modelled as an exponentially decaying function of height, i.e. $ρ = c_{1} e^{- c_{2} h}$ ( $c_{1}, c_{2}$ are constants (dimensionless), $c_{1} = 1.227$ , $c_{2} = 1.093 \times 10^{- 4}$ for $h < 9144 m$ , and $c_{1} = 1.754$ , $c_{2} = 1.49 \times 10^{- 4}$ for $h \geq 9144 m$ ) (Farina et al., 2002).

Process noise $w_{k}$ is modelled as a zero-mean white Gaussian process with non-singular covariance matrix (Bar-Shalom et al., 2001).

Q_{k} = [\begin{matrix} q_{1} θ & 0 & 0 & 0 \\ 0 & q_{1} θ & 0 & 0 \\ 0 & 0 & q_{1} θ & 0 \\ 0 & 0 & 0 & q_{2} T \end{matrix}], θ = [\begin{matrix} T^{3} / 3 & T^{2} / 2 \\ T^{2} / 2 & T \end{matrix}]

(57)

where parameters $q_{1}$ (in m²/s³) and $q_{2}$ (in kg²/(m⁴s)) control the amount of process noise in target dynamics and the ballistic coefficient, respectively.

The measurements, collected by the radar for state estimation, are the range $R_{k}$ , elevation $E_{k}$ and azimuth $A_{k}$ . According to the relative geometry, the measurement equation in the ENU-CS is given by

z_{k} = h (x_{k}) + v_{k}

(58)

where $z_{k} = {[\begin{matrix} \begin{matrix} R_{k} & E_{k} \end{matrix} & A_{k} \end{matrix}]}^{T}$ and $h (x_{k}) = {[\begin{matrix} \sqrt{x_{k}^{2} + y_{k}^{2} + z_{k}^{2}} & \arctan z_{k} / \sqrt{x_{k}^{2} + y_{k}^{2}} & \arctan y_{k} / x_{k} \end{matrix}]}^{T}$ , so

R_{k} = \sqrt{x_{k}^{2} + y_{k}^{2} + z_{k}^{2}} + v_{R}

(59)

E_{k} = \arctan z_{k} / \sqrt{x_{k}^{2} + y_{k}^{2}} \sqrt{x_{k}^{2} + y_{k}^{2}} + v_{E}

(60)

A_{k} = \arctan y_{k} / x_{k} x_{k} + v_{A}

(61)

The measurement process $v_{k} = {[\begin{matrix} v_{R} & v_{E} & v_{A} \end{matrix}]}^{T}$ is modelled as the zero-mean white Gaussian noise with the covariance matrix

N_{k} = diag ([\begin{matrix} σ_{R}^{2} & σ_{E}^{2} & σ_{A}^{2} \end{matrix}])

(62)

where $σ_{R}$ , $σ_{E}$ and $σ_{A}$ are the error standard deviations of range, elevation and azimuth. It is independent of the process noise $w_{k}$ and initial state $x_{0}$ .

Simulation and discussion

In the simulation, the parameter used is selected as $T = 0.1 s$ , $q_{1} = 5 m^{2} / s^{3}$ , $q_{2} = 5 k g^{2} / (m^{4} s)$ , the initial position and module of the velocity: $x_{0} = 232 km$ , $y_{0} = 232 km$ , $z_{0} = 90 km$ and $v_{0} = 3000 s$ ; the initial elevation and azimuth angle: $E_{0} = 210^{\circ}$ and $A_{0} = 45^{\circ}$ . The initial ballistic coefficient: $β_{0} = 4000 kg / m^{2}$ .The ballistic trajectory is generated using the fourth-order Runge–Kutta method. Then the true initial state can be obtained as x₀=[232 km −1837 m/s 232 km −1837 m/s 90 km −1500 m/s 4000 kg/m² ]^T, and the corresponding covariance is selected as $P_{0} = diag [100^{2} 5 0^{2} 1 00^{2} 5 0^{2} 1 00^{2} 5 0^{2} 20 0^{2}]$ .

The initial state estimate ${\hat{x}}_{0}$ is chosen randomly from ${\hat{x}}_{0} ~ N (x_{0}, P_{0})$ in each run. The ballistic coefficient $β$ is modelled as an independent-increment Wiener process with the same variance. The error standard deviations of the measurements are selected as $σ_{R} = 100 m$ , $σ_{E} = 0.017 rad$ and $σ_{A} = 0.017 rad$ .

To compare the performance of the various non-linear algorithms, we use the absolute error (AE) in the position, velocity and ballistic coefficient. The performance metrics of AE in the position at k time is defined as

A E_{p} (k) = \frac{1}{N} \sum_{i = 1}^{N} \sqrt{{(x_{k}^{(i)} - {\hat{x}}_{k}^{(i)})}^{2} + {(y_{k}^{(i)} - {\hat{y}}_{k}^{(i)})}^{2} + {(z_{k}^{(i)} - {\hat{z}}_{k}^{(i)})}^{2}}

(63)

where $(x^{(i)}, y^{(i)}, z^{(i)})$ and $({\hat{x}}^{(i)}, {\hat{y}}^{(i)}, {\hat{z}}^{(i)})$ are the true and estimated position at the ith Monte Carlo run, N is the Monte Carlo runs. Similarly to the AE in position, we may also write formulas of the AE in velocity and ballistic coefficient. Figures 2, 3 and 4 show the AEs for the UKF, SRCKF, iterative square-root cubature Kalman filter (ISRCKF) and LISRCKF algorithms in position, velocity and ballistic coefficient in the interval of 15–58 s. The parameters in the UKF algorithm are selected as $κ = 0$ , $α = 1$ , $β_UKF = 2$ (Julier, 2002). The ISRCKF algorithm uses the traditional iteration termination condition; and the iteration number selected in the LISRCKF and ISRCKF algorithms is 4. All performance curves were obtained by averaging over 100 independent Monte Carlo runs. All the filters were initialized with the same condition in each run. All the simulations were done in MATLAB on a ThinkPad PC with an Intel (R) CORE i5 M480 processor with a 2.67 GHz clock speed and 3 GB physical memory.

Observing from Figures 2 and 3, we can see that the errors of LISRCKF in position and velocity are far lower than those of SRCKF and UKF, and are lower than those of ISRCKF. Moreover, the LISRCKF needs 17.7 s to reduce the errors in position to below 500 m, the ISRCKF algorithm needs 34.6 s, and UKF and SRCKF algorithms need about 47.6 s, so the LISRCKF algorithm has faster convergence rate than the UKF, SRCKF and ISRCKF algorithms. The reason is that the latest measurement is involved in the step of measurement update of the LISRCKF and the iterative termination condition based on maximum likelihood surface is used. So the estimates provided by the LISRCKF in the position and velocity are markedly better than those of UKF, SRCKF and ISRCKF algorithms.

As to the estimation of the ballistic coefficient, in Figure 4, the AEs of the four filters have no improvement in the interval time (0<t<35 s) because of having less effective information about it from the noisy measurement. In the remained period (35 s<t<58 s), their errors for the ballistic coefficient are decreased. We can see that the AE in the LISRCKF algorithm decreases most rapidly of the four algorithms from Figure 4. So the ballistic coefficient estimate in the LISRCKF algorithm has the great improvement.

For further comparison of performance for the various filters, the average accumulated mean-square root error (AMSRE) in position is defined as follows:

\begin{matrix} AMSR E_{p} = \\ \frac{1}{N} \sum_{i = 1}^{N} \sqrt{[\frac{1}{M} \sum_{k = 1}^{M} [{(x_{k}^{(i)} - {\hat{x}}_{k}^{(i)})}^{2} + {(y_{k}^{(i)} - {\hat{y}}_{k}^{(i)})}^{2} + {(z_{k}^{(i)} - {\hat{z}}_{k}^{(i)})}^{2}]]} \end{matrix}

(64)

where M is the total number of measurement data points. The formulas of the AMSRE in velocity and ballistic coefficient can be written similarly to the AMSRE_p. The AMSREs in position, velocity and ballistic coefficient for the four filters are listed in Table 1, and the runtime of four algorithms is listed in Table 2. These figures were averaged over 100 runs.

Table 1.

Accumulated mean-square root error (AMSREs) of various algorithms.

Various filters	AMSRE_p, m	AMSRE_v, m/s	AMSRE_β, kg/m²
UKF	2647.837	300.525	152.060
SRCKF	2647.805	300.523	151.983
ISRCKF	1432.837	246.272	149.538
LISRCKF	1070.875	228.891	146.074

UKF, unscented Kalman filter; SRCKF, square root cubature Kalman filter; ISRCKF, iterative square root cubature Kalman filter; LISRCKF, likelihood-based iteration square root cubature Kalman filter.

Table 2.

Mean runtime for various algorithms.

Various filters	UKF	SRCKF	ISRCKF	LISRCKF
Time, s	1.0778	0.3415	0.5556	0.7468

UKF, unscented Kalman filter; SRCKF, square root cubature Kalman filter; ISRCKF, iterative square root cubature Kalman filter; LISRCKF, likelihood-based iteration square root cubature Kalman filter.

From Table 1, it is seen that the LISRCKF’s AMSRE in position reduces by about 60%, and its AMSRE in velocity reduces by about 25% compared with the UKF and SRCKF. Furthermore, compared with the ISRCKF, the AMSRE of the LISRCKF algorithm in position decreases by about 25% and its AMSRE in velocity decreases by about 7%. Table 1 shows LISRCKF’s AMSRE in the ballistic coefficient reduces marginally, but Figure 4 shows the LISRCKF’s error is lower than the other three filters in the interval of 40–58 s. Hence, the LISRCKF is to be preferred over the other filters in light of the AMSREs in the position, velocity and ballistic coefficient. The LISRCKF has better performance because it involves use of the latest measurement and use of iteration termination criteria based on the maximum likelihood surface.

Table 2 shows the runtime of the four algorithms, the runtime of the LISRCKF algorithm is far less than that of UKF, and is slightly larger than that of ISRCKF; however, the estimation accuracy of the LISRCKF is superior to that of ISRCKF. Though the runtime of SRCKF is shortest, its AMRSE is large to make the state estimation accuracy lower.

Meanwhile, from Figures 2, 3, 4 and Table 1, we can see the AEs and AMSREs of the SRCKF algorithm in the position, velocity and ballistic coefficient is marginally lower than those of the UKF.

Therefore, on the basis of the simulation results presented in Figures 2, 3, 4 and Table 1, one can draw a conclusion that the LISRCKF algorithm yields a superior performance over the UKF, SRCKF and ISRCKF algorithms, due to the use of the iterative measurement update and iterative termination criterion based on the maximum likelihood, and the deviation of UKF and SRCKF is very small.

Analysis of iteration number’s impact on the LISRCKF algorithm

In order to analyse the suitable iteration numbers selected in the LISRCKF algorithm, the LISRCKFs with various iteration numbers are applied to state estimation of RBT. Table 3 lists the mark of LISRCKF with various iteration numbers. Figures 5–7 show the errors of SRCKF and LISRCKFs with various numbers in position, velocity and ballistic coefficient. All performance curves and figures were also obtained by averaging over 100 independent Monte Carlo runs. Table 4 lists the AMSREs in position, velocity and ballistic coefficient for SRCKF and LISRCKFs with various iteration numbers.

Table 3.

The likelihood-based iteration square root cubature Kalman filter’s (LISRCKF’s) mark with various iteration numbers.

Iteration numbers	LISRCKF’s mark	Iteration number	LISRCKF’s mark
it_no=2	LISRCKF₁	it_no=8	LISRCKF₄
it_no=5	LISRCKF₂	it_no=9	LISRCKF₅
it_no=6	LISRCKF₃

Table 4.

Accumulated mean-square root error (AMSREs) of the likelihood-based iteration square root cubature Kalman filter (LISRCKF) with various iteration numbers.

Various algorithms	AMSRE_p, m	AMSRE_v, m/s	AMSRE_β, kg/m²
SRCKF	2674.638	306.322	160.452
LISRCKF₁	1944.108	275.680	156.286
LISRCKF₂	1262.888	239.957	153.202
LISRCKF₃	1157.747	233.486	152.808
LISRCKF₄	1008.717	223.695	152.295
LISRCKF₅	953.181	219.833	152.120

SRCKF, square root cubature Kalman filter.

From Figures 5–7, the errors in position and velocity for the LISRCKF with iteration number 2 begin to decrease largely compared with that of SRCKF. The errors in position, velocity and ballistic coefficient for LISRCKF with iteration number 5 are significantly reduced, and the LISRCKF algorithm has a very fast convergence rate. The errors of position, velocity and ballistic coefficient for LISRCKF decrease slowly when the iteration number is greater than 6; the reason is that the sequence generated has basically reached the maximum likelihood surface when the iteration termination condition is met. The more iteration number, the longer that LISRCKF algorithm runs. From the above analysis, iteration number 5 selected is suitable for LISRCKF, the more accurate state estimation is obtained and its runtime is shorter.

Observe from Table 4, the analysed results obtained are identical to those of Figures 5–7.

From the above analysis, iteration number 5 is the most appropriate for LISRCKF algorithm. The computational complexity of LISRCKF algorithm has the same level as the UKF algorithm when additive noise is considered in the system, but its speed is faster than that of the UKF.

Conclusion

In this study, we developed a more accurate and robust non-linear state estimator named the LISRCKF. The proposed algorithm inherits the virtues of the SRCKF and involves the iterative process in the step of the measurement update. The sequence obtained in the iterative measurement update is guaranteed to move up the likelihood surface using the iteration termination condition based on the maximum likelihood surface. So the LISRCKF has the advantage of free-derivative, easy computer and numerical stability. The LISRCKF is applied to state estimation of the RBT. Simulation results demonstrate that the performance of LISRCKF is superior to the UKF, SRCKF and ISRCKF, and the impact of iteration number on the state estimation was analysed. Iteration number 5 selected is the most suitable for the LISRCKF. So the LISRCKF is much more effective than the other filters and improves the performance of state estimation to a marked degree and has the fast convergence.

Footnotes

Acknowledgements

The authors would like to thank the anonymous reviewers and the editors for their many helpful suggestions.

Funding

The authors would like to thank the support of the National Natural Science Foundation of China (NSFC) under contract no. 60972146.

References

Arasaratnam

Haykin

(2009) Cubature Kalman filters. IEEE Transactions on Automatic Control 54: 1254–1269.

Arasaratnam

Haykin

Elliott

(2007) Discrete-time nonlinear filtering algorithms using Gauss–Hermite quadrature. Proceedings of the IEEE 95: 953–977.

Arulampalam

Maskell

Gordon

. (2002) A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing 50: 174–188.

Athans

Wishner

Bertolini

(1968) Suboptimal state estimation for continuous-time nonlinear system from discrete noisy measurements. IEEE Transactions on Automatic Control 13: 504–514.

Bar-Shalom

Kirubarajan

(2001) Estimation with Applications to Tracking and Navigation. New York: John Wiley & Sons, Inc.

Bell

Cathey

(1993) The iterated Kalman filter update as a Gauss–Newton method. IEEE Transactions on Automatic Control 38: 294–297.

Farina

Ristic

Benvenuti

(2002) Tracking a ballistic target: comparison of several nonlinear filters. IEEE Transactions on Aerospace and Electronic Systems 38: 854–867.

Gelb

(1974) Applied Optimal Estimation. Cambridge, MA: The MIT Press.

Gordon

Salmond

Smith

AFM

(1993) Novel approach to nonlinear non-Gaussian Bayesian state estimation. IEE Proceedings – F Radar and Signal Processing 140: 107–113.

10.

Grewal

Andrews

(2008) Kalman Filtering: Theory and Practice using Matlab. New York: John Wiley & Sons, Inc.

11.

Ito

Xiong

(2000) Gaussian filters for nonlinear filtering problems. IEEE Transactions on Automatic Control 45: 910–927.

12.

Johnston

Krishnamurthy

(2001) Derivation of a sawtooth iterated extended Kalman smoother via the AECM algorithm. IEEE Transactions on Signal Processing 49: 1899–1909.

13.

Julier

Uhlmann

Durrant-Whyte

(2000) A new method for the nonlinear transformation of means and covariances in filters and estimations. IEEE Transactions on Automatic Control 45: 477–482.

14.

Julier

(2002) The scaled unscented transformation. In: Proceedings of the American Control Conference, Anchorage, AK, pp. 4555–4559.

15.

Julier

Uhlmann

(2004) Unscented filtering and nonlinear estimation. Proceedings of the IEEE 92: 1958–1958.

16.

Julier

Uhlmann

Durrant-Whyte

(1995) A new approach for filtering nonlinear systems. Proceedings of the American Control Conference, Seattle, WA; American Automatic Control Council, pp. 1628–1632.

17.

Kaminski

Bryson

Schmidt

(1971) Discrete square root filtering: survey of current techniques. IEEE Transactions on Automatic Control AC16: 727–736.

18.

Lefebvre

Bruyninckx

De Schutter

(2004) Kalman filters for non-linear systems: a comparison of performance. International Journal of Control 77: 639–653.

19.

Jilkov

(2001) A survey of maneuvering target tracking – part II: Ballistic target models. Signal and Data Processing of Small Targets 4473: 559–581.

20.

Merwe

RVD

(2004) Sigma-point filters for probabilistic inference in dynamic state-space models. Ph.D., OGI School of Science and Engineering.

21.

Nørgaard

Poulsen

Ravan

(2000a) Advances in derivative-free state estimation for nonlinear systems. Technical University Denmark, Lyngby.

22.

Nørgaard

Poulsen

Ravn

(2000b) New developments in state estimation for nonlinear systems. Automatica 36: 1627–1638.

23.

Zhan

Wan

(2007) Iterated unscented Kalman filter for passive target tracking. IEEE Transactions on Aerospace and Electronic Systems 43: 1155–1163.