Maximum Likelihood-Based Iterated Divided Difference Filter for Nonlinear Systems from Discrete Noisy Measurements

Abstract

A new filter named the maximum likelihood-based iterated divided difference filter (MLIDDF) is developed to improve the low state estimation accuracy of nonlinear state estimation due to large initial estimation errors and nonlinearity of measurement equations. The MLIDDF algorithm is derivative-free and implemented only by calculating the functional evaluations. The MLIDDF algorithm involves the use of the iteration measurement update and the current measurement, and the iteration termination criterion based on maximum likelihood is introduced in the measurement update step, so the MLIDDF is guaranteed to produce a sequence estimate that moves up the maximum likelihood surface. In a simulation, its performance is compared against that of the unscented Kalman filter (UKF), divided difference filter (DDF), iterated unscented Kalman filter (IUKF) and iterated divided difference filter (IDDF) both using a traditional iteration strategy. Simulation results demonstrate that the accumulated mean-square root error for the MLIDDF algorithm in position is reduced by 63% compared to that of UKF and DDF algorithms, and by 7% compared to that of IUKF and IDDF algorithms. The new algorithm thus has better state estimation accuracy and a fast convergence rate.

1. Introduction

The problem of estimating the state of a nonlinear stochastic system from noisy measurement data has been the subject of considerable research interest during the past few years. Up to now the extended Kalman filter (EKF) has unquestionably been the dominating state estimation technique [1,2]. The EKF linearizes both the nonlinear process and the measurement dynamics with a first-order Taylor series expansion about the current state estimate. However, its accuracy depends heavily on the severity of nonlinearities. The EKF may introduce large errors and even give a divergent estimate when the nonlinearities become severe [3,4]. To improve the estimation accuracy, the second-order EKF proposed retains the Taylor series expansion up to the second term. The second–order EKF generally improves estimation accuracy, but at the expense of an increased computational burden [5]. Another attempt to improve the performance of the EKF involves the use of an iterative measurement update; the resulting algorithm is called the Iterated Extended Kalman filter (IEKF) [6]. The basic idea of IEKF is to linearize the measurement model around the updated state rather than the predicted state. This is achieved iteratively, and it involves the use of the current measurement. The IEKF has been proven to be more accurate on the condition that the state estimate is close enough to the true value, however, this is rarely the case in practice [7]. It was pointed out in [8] that the sequence of iterations generated by the IEKF and that generated by the Gauss-Newton method were identical, thus globally convergence was guaranteed. However, the Gauss-Newton method does not ensure that it goes up the likelihood surface [9,10]. Furthermore, EKF and IEKF require Jacobians, and the second-order KF requires Jacobians and Hessians. Calculation of Jacobians and Hessians is often numerically unstable and computationally intensive. In some system, the Jacobians and Hessians do not exit, which limits the applications of EKF, second-order EKF and IEKF.

Recently, there has been development in derivative-free state estimators. The finite difference has been used in the Kalman filter framework and the resulting filter is referred to as the finite difference filter (FDF) [11]. The FDF uses the first-order difference to approximate the derivative of the nonlinear function; it may introduce large state estimation errors due to a high nonlinearity, similar to the EKF. The unscented Kalman filter (UKF) proposed in [12,13] uses a minimal set of deterministically chosen sample points to capture the mean and covariance of a Gaussian density. When propagated through a nonlinear function, these points capture the true mean and covariance up to a second-order of the nonlinear function. However, the parameters used in the UKF are required to tune finely in order to prevent the propagation of non-positive definite covariance matrix for a state vector's dimension higher than three. Another Gaussian filter, named the divided difference filter (DDF) was introduced in [14] using multidimensional Stirling's interpolation formula. It is shown in [15] that the UKF and DDF algorithms are commonly referred to as sigma point filters due to the properties of deterministic sampling and weighted statistical estimation [16], but the covariance obtained in the DDF is more accurate than that in the UKF. The iterated UKF with the variable step (IUKF-VS) in [10] proposed improved the accuracy of state estimation but its runtime was large due to its computation of the sigma points. Lastly, a relatively new technique called the particle filter (PF) uses a set of randomly chosen samples with associated weights to approximate the posterior density [17] and its variants are presented in [18]. The large number of samples required often makes the use of PF computationally expensive, and the performance of PF is crucially dependent on the selection of the proposal distribution.

Table 1. Pro and cons for various filters.

Various filters	Advantage	Disadvantage
EKF	Less runtime	Calculation of Jacobian
Second-order EKF	high accuracy	Calculation of Jacobian and Hessian
IEKF	high accuracy	Calculation of Jacobian
UKF	Derivative-free	more runtime
IUKF-VS	Derivative-free, high accuracy	more runtime
FDF	Derivative-free, Less runtime	Lower accuracy
DDF	Derivative-free	Low accuracy
PF	Derivative-free	Heavy computational burden

lists the pro and cons of the above filters.

The DDF also shows its weakness in the state estimation due to the large initial error and high nonlinearity in the application for state estimation of maneuvering target in the air-traffic control and ballistic re-entry target. Emboldened by the superiority of DDF, the basic idea of the IEKF and the iteration termination condition based on maximum likelihood, we propose a new filter named the maximum likelihood based iterated divided difference Kalman filter (MLIDDF). The performance of the state estimation for MLIDDF is greatly improved when involving the use of the iteration measurement update in the MLIDDF and the use of the current measurement. The remainder of this paper is organized as follows: in Section 2, we develop the maximum likelihood surface based iterated divided difference Kalman filter (MLIDDF). Section 3 presents the applications of the MLIDDF to state estimation for maneuvering targets in air-traffic control and ballistic target re-entry applications and discuss the simulation results. Finally, Section 4 concludes the paper and presents our outlook on future work.

2. Development of Likelihood Surface Based Iterated Divided Difference Filter

2.1. Divided Difference Filter

Consider the nonlinear function:

$$\mathbf{\text{y}}=\mathbf{\text{f}}(\mathbf{\text{x}})$$

(1)

Assuming that the random variable x ∈ ℝ^n_x has Gaussian density with mean x̄ and covariance P_x. The following linear transformation of x is introduced:

$$\mathbf{\text{z}}={\mathbf{\text{S}}}_x^{-1}\mathbf{\text{x}}$$

(2)

The transformation matrix S_x is selected as a square Cholesky factor of the covariance matrix P_x such that ${\mathbf{\text{P}}}_x={\mathbf{\text{S}}}_x{\mathbf{\text{S}}}_x^T$, so the elements of z become mutually uncorrelated [14]. And the function f̃ is defined by:

$$\mathbf{\text{y}}=\mathbf{\text{f}}({\mathbf{\text{S}}}_x\mathbf{\text{z}})=\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\mathbf{\text{z}})$$

(3)

The multidimensional Stirling interpolation formula of Equation (3) about z̄ up to second-order terms is given by:

$$\mathbf{\text{y}}\approx \stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}})+{\stackrel{\sim }{D}}_{\mathrm{\Delta }z}\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}+\frac{1}{2!}{\stackrel{\sim }{D}}_{\mathrm{\Delta }z}^2\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}$$

(4)

The divided difference operators D̃_Δzf̃, ${\stackrel{\sim }{D}}_{\mathrm{\Delta }z}^2\stackrel{\sim }{f}$ are defined as:

$${\stackrel{\sim }{D}}_{\mathrm{\Delta }z}\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}=\frac{1}{l}(\sum _{j=1}^{n_x}\mathrm{\Delta }{z}_j{\mu }_j{\delta }_j)\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}})$$

(5)

$${\stackrel{\sim }{D}}_{\mathrm{\Delta }z}^2\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}=\frac{1}{l^2}(\sum _{j=1}^{n_x}\mathrm{\Delta }{z}_j^2{\delta }_j^2+\sum _{j=1}^{n_x}\sum _{\begin{array}{l}i=1\\ i\ne j\end{array}}^{n_x}\mathrm{\Delta }{z}_j\mathrm{\Delta }{z}_i({\mu }_j{\delta }_j)({\mu }_i{\delta }_i))\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}})$$

(6)

where Δz ≜ z - Z̄ and Δz_j is the j-th element of Δz. l denotes a selected interval length, the optimal setting of l is selected as $\sqrt{3}$ under the assumption that the estimation error are Gaussian and unbiased.

The partial operators μ and δ are defined as:

$${\mu }_j\mathbf{\text{f}}(\overline{\mathbf{\text{z}}})=\frac{1}{2}\left[\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}}+\frac{l}{2}{\mathbf{\text{e}}}_j)+\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}}-\frac{l}{2}{\mathbf{\text{e}}}_j)\right]$$

(7)

$${\delta }_j\mathbf{\text{f}}(\overline{\mathbf{\text{z}}})=\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}}+\frac{l}{2}{\mathbf{\text{e}}}_j)-\stackrel{\mathbf{\sim }}{\mathbf{\text{f}}}(\overline{\mathbf{\text{z}}}-\frac{l}{2}{\mathbf{\text{e}}}_j)$$

(8)

where e_j is the unit column vector.

We can obtain the approximate mean, covariance and cross-covariance of y using Equation (4):

$$\overline{\mathbf{\text{y}}}=E[\mathbf{\text{y}}]\approx \frac{l^2-n_x}{l^2}\mathbf{\text{f}}(\overline{\mathbf{\text{x}}})+\frac{1}{2l^2}\sum _{j=1}^{n_x}\left[\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})+\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j})\right]$$

(9)

$$\begin{array}{ll}{\mathbf{\text{P}}}_{yy}\hfill & =E[(\mathbf{\text{y}}-\overline{\mathbf{\text{y}}}){(\mathbf{\text{y}}-\overline{\mathbf{\text{y}}})}^T]\hfill \\ \hfill & \approx \frac{1}{4l^2}\sum _{j=1}^{n_x}\left[\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})-\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j})\right]\times {\left[\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})-\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j})\right]}^T\hfill \\ \hfill & +\frac{l^2-1}{4l^4}\sum _{j=1}^{n_x}\left[\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})+\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j})-2\mathbf{\text{f}}(\overline{\mathbf{\text{x}}})\right]\times {\left[\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})+\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j})-2\mathbf{\text{f}}(\overline{\mathbf{\text{x}}})\right]}^T\hfill \end{array}$$

(10)

$${\mathbf{\text{P}}}_{xy}=E[(\mathbf{\text{x}}-\overline{\mathbf{\text{x}}}){(\mathbf{\text{y}}-\overline{\mathbf{\text{y}}})}^T]\approx \frac{1}{2l}\sum _{j=1}^{n_x}{\mathbf{\text{s}}}_{x,j}{(\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}+l{\mathbf{\text{s}}}_{x,j})-\mathbf{\text{f}}(\overline{\mathbf{\text{x}}}-l{\mathbf{\text{s}}}_{x,j}))}^T$$

(11)

where s_x,j is j-th column of S_x.

Consider the state estimation problem of a nonlinear dynamics system with additive noise, the n_x-dimensional state vector x_k of the system evolves according to the nonlinear stochastic difference equation:

$${\mathbf{\text{x}}}_k=\mathbf{\text{f}}({\mathbf{\text{x}}}_{k-1})+{\mathbf{\text{w}}}_{k-1}$$

(12)

and the measurement equation is given as:

$${\mathbf{\text{z}}}_k=\mathbf{\text{h}}({\mathbf{\text{x}}}_k)+{\mathbf{\text{v}}}_k$$

(13)

w_k−1 and v_k are assumed i.i.d. and independent of current and past states, w_k−1 ∼ (0, Q_k−1), v_k ∼ (0, R_k).

Suppose the state distribution at k-1 time instant is x_k−1∼ (x̂_k−1, P_k−1), and a square Cholesky factor of P_k−1 is Ŝ_x,k−1. The divided difference filter (DDF) obtained with Equations (9)–(11) can be described as follows:

• Step 1. Time update

  • Calculate matrices containing the first- and second- divided difference on the estimated state x̂_k−1 at k-1 time:

    $${\mathbf{\text{S}}}_{x\widehat{x},k}^{(1)}=\{\frac{1}{2l}[\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}+l{\widehat{\mathbf{\text{s}}}}_{x,j})-\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}-l{\widehat{\mathbf{\text{s}}}}_{x,j})]\}$$

    (14)
    $${\mathbf{\text{S}}}_{x\widehat{x},k}^{(2)}=\{\frac{\sqrt{l^2-1}}{2l^2}[\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}+l{\widehat{\mathbf{\text{s}}}}_{x,j})+\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}-l{\widehat{\mathbf{\text{s}}}}_{x,j})-2\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1})]\}$$

    (15)

  • Evaluate the predicted state and square root of corresponding covariance:

    $${\overline{\mathbf{\text{x}}}}_k=\frac{l^2-n_x}{l^2}\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1})+\frac{1}{2l^2}\sum _{j=1}^{n_{\mathbf{\text{x}}}}\left[\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}+l{\widehat{\mathbf{\text{s}}}}_{x,j})+\mathbf{\text{f}}({\widehat{\mathbf{\text{x}}}}_{k-1}-l{\widehat{\mathbf{\text{s}}}}_{x,j})\right]$$

    (16)
    $${\overline{\mathbf{\text{S}}}}_{x,k}=\mathit{\text{Tria}}([\begin{array}{ccc}{\mathbf{\text{S}}}_{x\widehat{x},k}^{(1)}& {\mathbf{\text{S}}}_{w,k-1}& {\mathbf{\text{S}}}_{x\widehat{x},k}^{(2)}\end{array}])$$

    (17)
    ŝ_x,j is j-th column of Ŝ_x,k−1. Tria() is denoted as a general triagularization algorithm and S_w,k−1 denotes a square-root factor of Q_k−1 such that ${\mathbf{\text{Q}}}_{k-1}={\mathbf{\text{S}}}_{w,k-1}{\mathbf{\text{S}}}_{w,k-1}^T$.

• Step 2. Measurement update

  • Calculate matrices containing the first- and second-divided difference on the predicted state x̄_k:

    $${\mathbf{\text{S}}}_{z\overline{x},k}^{(1)}=\{\frac{1}{2l}[\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k+l{\overline{\mathbf{\text{s}}}}_{x,j})-\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k-l{\overline{\mathbf{\text{s}}}}_{x,j})]\}$$

    (18)
    $${\mathbf{\text{S}}}_{z\overline{x},k}^{(2)}=\{\frac{\sqrt{l^2-1}}{2l^2}[\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k+l{\overline{\mathbf{\text{s}}}}_{x,j})+\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k-l{\overline{\mathbf{\text{s}}}}_{x,j})-2\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k)]\}$$

    (19)
    where s̄_x,j is the j-th column of S̄_x,k.

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

    $${\overline{\mathbf{\text{z}}}}_k=\frac{l^2-n_x}{l^2}\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k)+\frac{1}{2l^2}\sum _{j=1}^{n_x}\left[\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k+l{\overline{\mathbf{\text{s}}}}_{x,j})+\mathbf{\text{h}}({\overline{\mathbf{\text{x}}}}_k-l{\overline{\mathbf{\text{s}}}}_{x,j})\right]$$

    (20)
    $${\mathbf{\text{S}}}_{zz,k}=\mathit{\text{Tria}}([\begin{array}{ccc}{\mathbf{\text{S}}}_{z\overline{x},k}^{(1)}& {\mathbf{\text{S}}}_{v,k}& {\mathbf{\text{S}}}_{z\overline{x},k}^{(2)}\end{array}])$$

    (21)
    $${\mathbf{\text{P}}}_{xz,k}={\overline{\mathbf{\text{S}}}}_{x,k}{({\mathbf{\text{S}}}_{z\overline{x},k}^{(1)})}^T$$

    (22)
    here S_v,k denotes a square root of R_k such that ${\mathbf{\text{R}}}_k={\mathbf{\text{S}}}_{v,k}{\mathbf{\text{S}}}_{v,k}^T$.

  • Evaluate the gain, state estimation and square root of corresponding covariance at k time:

    $${\mathbf{\text{K}}}_k=({\mathbf{\text{P}}}_{xz,,k}/{\mathbf{\text{S}}}_{zz,k}^T)/{\mathbf{\text{S}}}_{zz,k}$$

    (23)
    $${\widehat{\mathbf{\text{x}}}}_k={\overline{\mathbf{\text{x}}}}_k+{\mathbf{\text{K}}}_k({\mathbf{\text{z}}}_k-{\overline{\mathbf{\text{z}}}}_k)$$

    (24)
    $${\widehat{\mathbf{\text{S}}}}_{x,k}=\mathit{\text{Tria}}([\left[\begin{array}{ccc}({\overline{\mathbf{\text{S}}}}_{x,k}-{\mathbf{\text{K}}}_k{\mathbf{\text{S}}}_{z\overline{x},k}^{(1)})& {\mathbf{\text{K}}}_k{\mathbf{\text{S}}}_{v,k}& {\mathbf{\text{K}}}_k{\mathbf{\text{S}}}_{z\overline{x},k}^{(2)}\end{array}\right]])$$

    (25)
    here, the symbol “/” represents the matrix right divide operator.

2.2. Refining the Measurement Update Based on Divided Difference

Consider x̄_k and current measurement z_k as realization of independent random vectors with multivariate normal distributions, e.g., x̄_k ∼ (x_k, P̄_k) and z_k ∼ (h(x_k), R_k). For convenience, the two vectors are formed to a single augmented one Z = [z_k x̄_k]^T. According to the independent assumption, we have:

$$\mathbf{\text{Z}}\sim \mathcal{N}(\mathbf{\text{g}}({\mathbf{\text{x}}}_k),\stackrel{\mathbf{\sim }}{\mathbf{\text{Q}}})$$

(26)

Here:

$$\mathbf{\text{g}}({\mathbf{\text{x}}}_k)={[\mathbf{\text{h}}(\mathbf{\text{x}}),{\mathbf{\text{x}}}_k]}^T,\stackrel{\mathbf{\sim }}{\mathbf{\text{Q}}}=\left[\begin{array}{cc}{\mathbf{\text{R}}}_k& 0\\ 0& {\overline{\mathbf{\text{P}}}}_k\end{array}\right]$$

(27)

The update measurement problem becomes the one that computing the optimal state estimation and corresponding covariance given Z, g and Q̃.

Defining the objective function:

$$\mathbf{\text{f}}(\mathbf{\text{ξ}})=\frac{1}{2}{\Vert \mathbf{\text{r}}(\mathbf{\text{ξ}})\Vert }^2$$

(28)

where r(ξ) = S(Z − g(ξ)), and r(ξ) is the second-order differentiable. Then the above update measurement becomes clearly a non-linear least squares problem.

Assuming the i-iterate is ${\widehat{\text{x}}}_{\text{k}}^{(i+1)}$, we can obtain the following equation [8]:

$${\widehat{\mathbf{\text{x}}}}_k^{(i+1)}={\overline{\mathbf{\text{x}}}}_k+{\mathbf{\text{K}}}_k^{(i)}\left({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)})+{\mathbf{\text{H}}}_k^{(i)}\left({\overline{\mathbf{\text{x}}}}_k-{\widehat{\mathbf{\text{x}}}}_k^{(i)}\right)\right)$$

(29)

where ${\mathbf{\text{K}}}_k^{(i)}={\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)T}{({\mathbf{\text{H}}}_k^{(i)T}{\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)}+{\mathbf{\text{R}}}_k)}^{-1}$.

We know the sequence of iterates generated by the IEKF and that generated by the Gauss-Newton method were identical, thus globally convergent was guaranteed. The initial state x̄_k is included in the measurement update, the value of x̄_k has a direct and large effect on the final state estimation. When the measurement model fully observes the state, the estimated state ${\widehat{\mathbf{\text{x}}}}_{\text{k}}^{(\text{i})}$ is more approximate to the true state than the predicted state x̄_k [7]. Substituting x̄_k by ${\widehat{\mathbf{\text{x}}}}_{\text{k}}^{(\text{i})}$ into Equation (29), hence, the following iterative formula is obtained:

$${\mathbf{\text{x}}}_k^{(i+1)}={\mathbf{\text{x}}}_k^{(i)}+{\mathbf{\text{K}}}_k^{(i)}({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\mathbf{\text{x}}}_k^{(i)}))$$

(30)

Compared to Equation (29), the Equation (30) is simpler, and the two equations are identical when there is a single iteration.

Now we consider the gain:

$${\mathbf{\text{K}}}_k^{(i)}={\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)T}{\left({\mathbf{\text{H}}}_k^{(i)T}{\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)}+{\mathbf{\text{R}}}_k\right)}^{-1}$$

(31)

where the terms ${\mathbf{\text{H}}}_k^{(i)T}{\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)}+{\mathbf{\text{R}}}_k$ and ${\overline{\mathbf{\text{P}}}}_k{\mathbf{\text{H}}}_k^{(i)T}$ are approximate innovation covariance and cross-covariance obtained by linearizing the measurement equation:

$${\mathbf{\text{P}}}_{zz,k}^{(i)}={\mathbf{\text{H}}}_k^{(i)}{\overline{\mathbf{\text{P}}}}_k{({\mathbf{\text{H}}}_k^{(i)})}^T+{\mathbf{\text{R}}}_k$$

(32)

$${\mathbf{\text{P}}}_{xz,k}^{(i)}={\overline{\mathbf{\text{P}}}}_k{({\mathbf{\text{H}}}_k^{(i)})}^T$$

(33)

The terms of Equations (32) and (33) are achieved by expanding the measurement Equation (13) up to a first-order Taylor term so that the linearized error is introduced due to the high-order truncated terms. As for the highly nonlinear measurement equation, the accuracy for state estimation is decreased if the linearized error is only propagated in the Equation (30). To decrease the propagated error, we can recalculate Equations (18)–(22) to obtain the terms ${\mathbf{\text{P}}}_{zz}^{(i)}$, ${\mathbf{\text{P}}}_{xz}^{(i)}$ in the following way:

$${\mathbf{\text{S}}}_{zx,k}^{(1)(i)}=\{\frac{1}{2l}[\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)}+l{\widehat{\mathbf{\text{s}}}}_{x,j}^{(i)})-\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)}-l{\widehat{\mathbf{\text{s}}}}_{x,j}^{(i)})]\}$$

(34)

$${\mathbf{\text{S}}}_{zx,k}^{(2)(i)}=\{\frac{\sqrt{l^2-1}}{2l^2}[\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)}+l{\widehat{\mathbf{\text{s}}}}_{x,j}^{(i)})+\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)}-l{\widehat{\mathbf{\text{s}}}}_{x,j}^{(i)})-2\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)})]\}$$

(35)

$${\mathbf{\text{S}}}_{zz,k}^{(i)}=\mathit{\text{Tria}}(\left[\begin{array}{ccc}{\mathbf{\text{S}}}_{zx,k}^{(1)(i)}& {\mathbf{\text{S}}}_{v,k}& {\mathbf{\text{S}}}_{zx,k}^{(2)(i)}\end{array}\right])$$

(36)

$${\text{P}}_{zz,k}^{(i)}={\mathbf{\text{S}}}_{zz,k}^{(i)}{\mathbf{\text{S}}}_{zz,k}^{(i)T}$$

(37)

$${\mathbf{\text{P}}}_{xz,k}^{(i)}={\mathbf{\text{S}}}_k^{(i)}{\mathbf{\text{S}}}_{zx,k}^{(1)(i)}$$

(38)

where ${\mathbf{\text{S}}}_k^{(i)}$ is a square Cholesky factor of the covariance ${\mathbf{\text{P}}}_k^{(i)}$.

Hence, we can obtain the following iterative formula:

$${\widehat{\mathbf{\text{x}}}}_k^{(i+1)}={\mathbf{\text{x}}}_k^{(i)}+{\mathbf{\text{P}}}_{xz}^{(i)}{({\mathbf{\text{P}}}_{zz}^{(i)})}^{-1}\left[{\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)})\right]$$

(39)

2.3. Maximum Likelihood Based Iteration Termination Criterion

In the measurement update step of the IEKF algorithm, the inequality $\Vert {\widehat{\mathbf{\text{x}}}}_k^{(i+1)}-{\widehat{\mathbf{\text{x}}}}_k^{(i)}\Vert \le \epsilon$ is used as the criterion to terminate the iteration procedure, where ε is the predetermined threshold. The threshold ε is crucial to successfully using the IEKF algorithm, but selecting a proper value of ε is difficult [10]. The sequence of iterations generated according to the above termination condition has the property of global convergence; however, it is not guaranteed to move up the likelihood surface, so an iteration termination criterion based on maximum likelihood surface is introduced.

Consider ${\widehat{\mathbf{\text{x}}}}_k^{(i)}$ and z_k as the realization of independent random vectors with multivariate normal distributions, i.e., ${\widehat{\mathbf{\text{x}}}}_k^{(i)}\sim \mathcal{N}({\mathbf{\text{x}}}_k,{\mathbf{\text{P}}}_k^{(i)})$, z_k ∼ (h(x_k),R_k). The likelihood function of the two vectors x_k and z_k is defined as:

$$\mathrm{\Lambda }({\mathbf{\text{x}}}_k,{\mathbf{\text{y}}}_k)=\mathit{\text{const}}•exp\{-0.5\left[{({\mathbf{\text{x}}}_k-{\widehat{\mathbf{\text{x}}}}_k^{(i)})}^T{\mathbf{\text{P}}}_k^{(-i)}({\mathbf{\text{x}}}_k-{\widehat{\mathbf{\text{x}}}}_k^{(i)})+{({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\mathbf{\text{x}}}_k))}^T{\mathbf{\text{R}}}_k^{-1}({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\mathbf{\text{x}}}_k))\right]\}$$

(40)

where ${\mathbf{\text{P}}}_{\text{k}}^{(-i)}={({\mathbf{\text{P}}}_{\text{k}}^{(i)})}^{-1}$.

Meanwhile, the likelihood surface is defined as follows:

$$J({\mathbf{\text{x}}}_k)={({\widehat{\mathbf{\text{x}}}}_k^{(i)}-{\mathbf{\text{x}}}_k)}^T{\mathbf{\text{P}}}_k^{(-i)}({\widehat{\mathbf{\text{x}}}}_k^{(i)}-{\mathbf{\text{x}}}_k)+{({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\mathbf{\text{x}}}_k))}^T{\mathbf{\text{R}}}_k^{-1}({\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\mathbf{\text{x}}}_k))$$

(41)

We know that the solution that maximizes the likelihood function is equivalent to minimizing the cost function J(x_k). The optimal value of J(x_k) is difficult to obtain, but the following inequality holds:

$$J({\widehat{\mathbf{\text{x}}}}_k^{(i+1)})<J({\widehat{\mathbf{\text{x}}}}_k^{(i)})$$

(42)

We say that $J({\widehat{\mathbf{\text{x}}}}_k^{(i+1)})$ is close to the maximum likelihood surface than $J({\widehat{\mathbf{\text{x}}}}_k^{(i+1)})$, equivalently, ${\widehat{\mathbf{\text{x}}}}_k^{(i+1)}$ has a more accurate approximation than ${\widehat{\mathbf{\text{x}}}}_k^{(i)}$ to the minimum value of J(x_k) [10]. Extending Equation (42) and using ${\mathbf{\text{P}}}_k^{(i)}={\mathbf{\text{S}}}_k^{(i)}{\mathbf{\text{S}}}_k^{(i)T}$, we immediately obtain the following inequality:

$${({\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)})}^T{({\mathbf{\text{S}}}_k^{(i)}{\mathbf{\text{S}}}_k^{(i)T})}^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)}+{({\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)})}^T{\mathbf{\text{R}}}_k^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)}<{({\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i)})}^T{\mathbf{\text{R}}}_k^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i)}$$

(43)

where ${\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)}$ and ${\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)}$ and are defined as:

$${\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)}={\widehat{\mathbf{\text{x}}}}_k^{(i+1)}-{\widehat{\mathbf{\text{x}}}}_k^{(i)}$$

(44)

$${\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)}={\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i+1)})$$

(45)

The sequence generated is guaranteed to go up the likelihood surface using Equation (43) as the criterion to iteration termination.

2.4. Maximum Likelihood Based Iterated Divided Difference Filter

We have now arrived at the central issue of this paper, namely, the maximum likelihood based iterated divided difference filter. Enlightened by the development of IEKF and the superiority of DDF, we can derive the maximum likelihood based iterated divided difference filter (MLIDDF) which involves the use of the iteration measurement update and the current measurement. But in view of the potential problems exhibited by the IEKF, we shall refine the covariance and cross-covariance based on divided difference and use the termination criterion which guarantees the sequence obtained moves up the maximum likelihood surface. The MLIDDF is described as follows:

• Step 1. Time update

  Evaluate the predicted state x̄_k and square Cholesky factor S̄_k of the corresponding covariance P̄_k using the Equations (14)–(17).

• Step 2. Measurement update

  Let ${\widehat{\mathbf{\text{x}}}}_k^{(0)}={\overline{\mathbf{\text{x}}}}_k$ and ${\mathbf{\text{S}}}_k^{(0)}={\overline{\mathbf{\text{S}}}}_k$. Suppose that the i-th iterates are ${\widehat{\mathbf{\text{x}}}}_k^{(i)}$ and ${\mathbf{\text{S}}}_k^{(i)}$.

  • Evaluate the first- and second-order difference matrices using Equation (34) and (35).

  • Evaluate the square root of innovation covariance and cross-covariance:

    $${\mathbf{\text{S}}}_{zz,k}^{(i)}=\mathit{\text{Tria}}(\left[\begin{array}{ccc}{\mathbf{\text{S}}}_{zx,k}^{(1)(i)}& {\mathbf{\text{S}}}_{v,k}& {\mathbf{\text{S}}}_{zx,k}^{(2)(i)}\end{array}\right])$$

    (46)
    $${\mathbf{\text{P}}}_{xz,k}^{(i)}={\mathbf{\text{S}}}_k^{(i)}{\mathbf{\text{S}}}_{zx,k}^{(1)(i)}$$

    (47)

  • Evaluate the gain:

    $${\mathbf{\text{K}}}_k^{(i)}={\mathbf{\text{P}}}_{xz}^{(i)}/{\mathbf{\text{S}}}_{zz,k}^{(i)T}/{\mathbf{\text{S}}}_{zz,k}^{(i)}$$

    (48)

  • Evaluate the state and the square root of corresponding covariance:

    $${\widehat{\mathbf{\text{x}}}}_k^{(i+1)}={\widehat{\mathbf{\text{x}}}}_k^{(i)}+{\mathbf{\text{K}}}_k^{(i)}[{\mathbf{\text{z}}}_k-\mathbf{\text{h}}({\widehat{\mathbf{\text{x}}}}_k^{(i)})]$$

    (49)
    $${\mathbf{\text{S}}}_k^{(i+1)}=\mathit{\text{Tria}}(\left[\begin{array}{ccc}({\overline{\mathbf{\text{S}}}}_{x,k}-{\mathbf{\text{K}}}_k^{(i)}{\mathbf{\text{S}}}_{zx,k}^{(1)(i)})& {\mathbf{\text{K}}}_k^{(i)}{\mathbf{\text{S}}}_v& {\mathbf{\text{K}}}_k^{(i)}{\mathbf{\text{S}}}_{zx,k}^{(2)(i)}\end{array}\right])$$

    (50)
    where ${\widehat{\mathbf{\text{s}}}}_{x,j}^{(i)}$ is the j-th column of ${\mathbf{\text{S}}}_k^{(i)}$.

• Step 3. If the following inequality holds:

  $${({\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)})}^T{({\mathbf{\text{S}}}_k^{(i)}{\mathbf{\text{S}}}_k^{(i)T})}^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)}+{({\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)})}^T{\mathbf{\text{R}}}_k^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i+1)}<{({\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i)})}^T{\mathbf{\text{R}}}_k^{-1}{\stackrel{\mathbf{\sim }}{\mathbf{\text{z}}}}_k^{(i)}$$

  (51)

  The iteration returns to Step 2; otherwise continue to Step 4. ${\stackrel{\mathbf{\sim }}{\mathbf{\text{x}}}}_k^{(i+1)}$, ${\stackrel{\sim }{\mathbf{\text{z}}}}_k^{(i+1)}$ are defined in the Equations (44) and (45).

• Step 4. If the inequality is not satisfied or if i is too large (i > N_max), and the ultimate state estimation and square root of corresponding covariance at k time instant are:

  $${\widehat{\mathbf{\text{x}}}}_k={\widehat{\mathbf{\text{x}}}}_k^{(N_{max})}$$

  (52)
  $${\mathbf{\text{S}}}_k={\mathbf{\text{S}}}_k^{(N_{max})}$$

  (53)

The MLIDDF algorithm has the virtues of free-derivative and better numerical stability. The measurement update of MLIDDF algorithm is transformed to a nonlinear least-square problem; the optimum state estimation and covariance are solved using Gauss-Newton method, so the MLIDDF algorithm has the same global convergence as the Gauss-Newton method. Moreover, the iteration termination condition that makes the sequence move up the maximum likelihood surface is used in the measurement update process.

3. Simulation and Analysis

In this section, we reported the experimental results obtained by applying the MLIDDF to the nonlinear state estimation of a maneuvering target in an air-traffic control scenario and a ballistic target re-entry scenario. To demonstrate the performance of the MLIDDF algorithm we compared its performance against the UKF, DDF, and the iterated UKF (IUKF) and iterated DDF (IDDF), both using a traditional iteration strategy.

3.1. Maneuvering Target Tracking in the Air-Traffic Control Scenario

We consider a typical air-traffic control scenario, where an aircraft executes a maneuvering turn in a horizontal plane at a constant, but unknown turn rate Ω. The kinematics of the turning motion can be modeled by the following nonlinear process equation [2,19]:

$${\mathbf{\text{x}}}_k=\left(\begin{array}{ccccc}1& \frac{sin\mathrm{\Omega }T}{\mathrm{\Omega }}& 0& -\frac{1-cos\mathrm{\Omega }T}{\mathrm{\Omega }}& 0\\ 0& cos\mathrm{\Omega }T& 0& -sin\mathrm{\Omega }T& 0\\ 0& \frac{1-cos\mathrm{\Omega }T}{\mathrm{\Omega }}& 1& \frac{sin\mathrm{\Omega }T}{\mathrm{\Omega }}& 0\\ 0& sin\mathrm{\Omega }T& 0& cos\mathrm{\Omega }T& 0\\ 0& 0& 0& 0& 1\end{array}\right){\mathbf{\text{x}}}_{k-1}+{\mathbf{\text{w}}}_{k-1}$$

(54)

where the state of the aircraft x = [x ẋ y ẏ Ω]; x and y denote positions, and ẋ and ẏ denote velocities in the x and y directions, respectively; T is the time-interval between two consecutive measurements;

The process nose w_k−1∼ (0, Q) with a nonsingular covariance:

$$\mathbf{\text{Q}}=\left(\begin{array}{ccc}{q}_1\mathbf{\text{π}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& q_1\mathbf{\text{π}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& q_{2}T\end{array}\right),\mathbf{\text{π}}=\left(\begin{array}{cc}{T}^3/3& T^2/2\\ T^2/2& T\end{array}\right)$$

(55)

The parameters q₁ and q₂ are related to process noise intensities. A radar is fixed at the origin of the place and equipped to measure the range, r and bearing, θ. Hence, the measurement equation is written:

$$\left(\begin{array}{c}{r}_k\\ {\theta }_k\end{array}\right)=\left(\begin{array}{c}\sqrt{x_k^2+y_k^2}\\ {tan}^{-1}(y_k/x_k)\end{array}\right)+{\mathbf{\text{v}}}_k$$

(56)

where the measurement noise v_k ∼ (0, R) and $\mathbf{\text{R}}=\text{diag}([{\sigma }_r^2{\sigma }_{\theta }^2])$.

The parameters used in this simulation were the same as those in [19]. To tracking the maneuvering aircraft we use the proposed MLIDDF algorithm and compare its performance against the DDF. We use the root-mean square error (RMSE) of the position, velocity and turn rate to compare the performances of two nonlinear filters. For a fair comparison, we make 250 independent Monte Carlo runs. The RMSE in position at time k is defined as:

$${\mathit{\text{RMSE}}}_p(k)=\sqrt{\frac{1}{N}\sum _{i=1}^N{(x_k^{(i)}-{\widehat{x}}_k^{(i)})}^2+{(y_k^{(i)}-{\widehat{y}}_k^{(i)})}^2}$$

(57)

where $(x_k^{(i)},y_k^{(i)})$ and $({\widehat{x}}_k^{(i)},{\widehat{y}}_k^{(i)})$ are the true and estimated positions at the i-th Monte Carlo run. Similarly to RMSE in position, we may also write formulas of RMSE in velocity and turn rate.

Owing to the fact that the filters is sensitive to initial state estimation, Figures 1–3 show the RMSEs in position, velocity and turn rate, respectively for DDF and MLIDDF in an interval of 50–100 s. As can be seen in Figures 1–3, the MLIDDF significantly outperforms the DDF.

In order to analyze the impact of iteration numbers on performance of the MLIDDF algorithm, the MLIDDFs with various iteration numbers are applied to position estimation of maneuvering target. Figure 4 shows the RMSEs of DDF and MLIDDFs with various numbers in position. From Figure 4, the RMSE in position for the MLIDDF with iteration number 2 begins to decrease largely comparable to that of DDF. The RMSE in position for MLIDDF with iteration number 5 significantly reduce, and the MLIDDF algorithm has very fast convergence rate. The RMSE of the position for MLIDDF decreases slowly when the iteration number is greater than 8; the reason is that the sequence generated has basically reached the maximum likelihood surface when the iteration termination condition is met.

3.2. State Estimation of Reentry Ballistic Target

3.2.1. Simulation Scene

The relative location of the ballistic target re-entry (BTR) and the radar are shown in Figure 5. The inertial coordinate system (ECI-CS) Ox_Iy_Iz_I shown in Figure 5 is a right-handed system with the origin O at the Earth's center, axis Ox_I pointing in the vernal equinox direction, axis Oz_I pointing in the direction of the North Pole N. Its fundamental plane Ox_Iy_I coincides with the Earth's equatorial plane.

Assume that the radar is situated at the surface of the Earth and considering the orthogonal coordinate reference system named East-North-Up coordinates system (ENU-CS) O_sxyz has its origin at the location of the radar. In this system, 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. Assuming that the Earth is spherical and non-rotating and the forces acting on the target are only gravity and drag [20], we can derive the following state equation according to the kinematics of the ballistic re-entry object in the reentry phase in ENU-CS:

$${\mathbf{\text{x}}}_k=\mathbf{\text{Φ}}{\mathbf{\text{x}}}_{k-1}+\mathbf{\text{G}}\mathbf{\psi }({\mathbf{\text{x}}}_{k-1})+{\mathbf{\text{w}}}_{k-1}$$

(58)

where x_k = [x_k ẋ_k y_k ẏ_k z_k ż_k β_k]T is the state of ballistic target re-entry, and:

$$\mathbf{\text{Φ}}=\left[\begin{array}{cccc}\varphi & \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \varphi & \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \varphi & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& 1\end{array}\right],\varphi =\left[\begin{array}{cc}1& T\\ 0& 1\end{array}\right]$$

(59)

$$\mathbf{\text{G}}=\left[\begin{array}{ccc}\tau & \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \tau & \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \tau \\ \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right],\tau =\left[\begin{array}{c}{T}^2/2\\ T\end{array}\right]$$

(60)

$$\mathbf{\psi }({\mathbf{\text{x}}}_{k-1})=\left[\begin{array}{c}-\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}{\dot{x}}_{k-1}-\frac{\mu x_{k-1}}{r_{k-1}^3}\\ -\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}{\dot{y}}_{k-1}-\frac{\mu y_{k-1}}{r_{k-1}^3}\\ -\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}{\dot{z}}_{k-1}-\frac{\mu (z_{k-1}+R_e)}{r_{k-1}^3}\end{array}\right]$$

(61)

where

$$\{\begin{array}{l}{r}_{k-1}=\sqrt{x_{k-1}^2+y_{k-1}^2+{(z_{k-1}+R_e)}^2}\\ V_{k-1}=\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2+{\dot{z}}_{k-1}^2}\\ h_{k-1}=r_{k-1}-R_e\end{array}$$

(62)

T is the time interval between radar measurements, β is the ballistic coefficient (kg/m²), μ and R_e are the Earth's gravitational constant and average Earth radius, respectively. Below 90 km at height, the air density ρ(h) is approximately modeled as an exponentially decaying function of height, ρ(h) = c₁e^−c₂h (c₁, c₂ are constant, specifically, c₁ = 1.227, c₂ = 1.093 × 10⁻⁴ for h < 9,144 m, and c₁ = 1.754, c₂ = 1.49 × 10⁻⁴ for h ≥ 9,144 m) [21].

Process noise w_k is assumed to be white noise with zero mean; its covariance is approximately modeled as [2]:

$${\mathbf{\text{Q}}}_k=\left[\begin{array}{cccc}{q}_1{\mathbf{\text{θ}}}_1& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& q_1{\mathbf{\text{θ}}}_1& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& q_1{\mathbf{\text{θ}}}_1& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& q_{2}T\end{array}\right],\mathbf{\text{θ}}=\left[\begin{array}{cc}\stackrel{T^3}{}/\underset{3}{}& \stackrel{T^2}{}/\underset{2}{}\\ \stackrel{T^2}{}/\underset{2}{}& T\end{array}\right]$$

(63)

where q₁ (in m²/s³) and q₂ (in kg²/m⁴s) are the intensity of noise.

According to relative geometry, the measurement equation in the ENU-CS is described as:

$${\mathbf{\text{z}}}_k=\mathbf{\text{h}}({\mathbf{\text{x}}}_k)+{\mathbf{\text{v}}}_k$$

(64)

here z_k = [R_k E_k A_k]^T, and:

$$R_k=\sqrt{x_k^2+y_k^2+z_k^2}+v_R$$

(65)

$$E_k=ta{n}^{-1}{z}_k/\sqrt{x_k^2+y_k^2}+v_E$$

(66)

$$A_k=ta{n}^{-1}{y}_k/x_k+v_A$$

(67)

The measurement noise v_k is assumed to be white noise with zero mean and covariance:

$${\mathbf{\text{N}}}_k=\mathit{\text{diag}}(\left[\begin{array}{ccc}{\sigma }_R^2& {\sigma }_E^2& {\sigma }_A^2\end{array}\right])$$

(68)

where σ_R, σ_E, σ_A, are the error standard deviations of range, elevation and azimuth, respectively. It is independent of the process noise w_k and initial state x₀.

3.2.2. Numerical Results and Analysis

The parameters used in simulation were: T = 0.1 s, q₁ = 5 m²/s³, q₂ = 5 kg²/m⁴s. The initial position and magnitude of the velocity: x₀ = 232 km, y₀ = 232 km, z₀ = 90 km and v₀ = 3,000 m/s, and the initial elevation and azimuth angle: E₀ = 7π/6 and A₀ = π/4. The ballistic coefficient was selected as β = 4,000 kg/m². The error standard deviations of the measurements were selected as σ_R = 100 m, σ_E = 0.017 rad, and σ_A = 0.017 rad. A threshold ε = 10 was set in the IUKF and IDDF, and a maximum iteration number N_max = 8 was predetermined.

From the above parameters given, we can obtain the initial true state:

$$x_0=[232km-1837m/s232km-1837m/s90km-1500m/s4000kg/m^2]$$

and we select the corresponding covariance as: P₀ = diag[100² 50² 100² 50² 100² 50² 200²].

To compare the performance of the various filter algorithms, we also use RMSE in the position, velocity and ballistic coefficient. The position RMSE at k time of the ballistic target re-entry is defined as:

$${\mathit{\text{RMSE}}}_p(k)=\sqrt{\frac{1}{N}\sum _{i=1}^N({(x_k^{(i)}-{\widehat{x}}_k^{(i)})}^2+{(y_k^{(i)}-{\widehat{y}}_k^{(i)})}^2+{(z_k^{(i)}-{\widehat{z}}_k^{(i)})}^2)}$$

(69)

where $(x_k^{(i)},y_k^{(i)},z_k^{(i)})$ and $({\widehat{x}}_k^{(i)},{\widehat{y}}_k^{(i)},{\widehat{z}}_k^{(i)})$ are the true and estimated position at the i-th Monte Carlo run, N is the Monte Carlo runs. Similarity to the RMSE in position we may also write formulas of the RMSE in velocity and ballistic coefficient.

Figures 6–8 shows the position, velocity and ballistic coefficient RMSEs, respectively, for the various filters in an interval of 15–58 s. The initial state estimate x̂₀ is chosen randomly from x̂₀ ∼ (x₀, P₀) in each run. All the filters are initialized with the same condition in each run. We make 100 independent Monte Carlo runs.

From Figure 6 we can see that the position RMSE of the MLIDDF is much less than the UKF and DDF because of the use of the current measurement in the step of the iteration measurement update of the MLIDDF, and is less than the IUKF and IDDF owing to involving the proposed iteration strategy, so the estimates provided by the MLIDDF are markedly better than those of the UKF and DDF algorithms, and are better than those of IUKF and IDDF algorithms. The MLIDDF also shows a significant improvement over the other filters in the estimation of the velocity, as evidenced by Figure 7.

As to the estimation of the ballistic coefficient, in the Figure 8, we can see that there is no improvement in the RMSE in the initial interval of the observation period (t < 35 s) because there is no effective information about it, while in the remaining period (35 s < t < 58 s) the ballistic coefficient RMSE is decreased because the effective information about the ballistic coefficient from the latest measurement is fully used. Especially, Figure 8 illustrates that toward the end of the trajectory the estimates provided by the MLIDDF are markedly better than those of the UKF, DDF, IUKF and IDDF algorithms.

Meanwhile, we observe from Figures 6–8 that the UKF and DDF have almost the same performance in the problem and the performance of IUKF and IDDF algorithms are almost identical.

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

$${\mathit{\text{AMSRE}}}_p=\frac{1}{N}\sum _{i=1}^N\sqrt{\left[\frac{1}{M}\sum _{k=1}^M[{(x_k^{(i)}-{\widehat{x}}_k^{(i)})}^2+{(y_k^{(i)}-{\widehat{y}}_k^{(i)})}^2+{(z_k^{(i)}-{\widehat{z}}_k^{(i)})}^2]\right]}$$

(70)

where M is the total number of measurement data points. The formulas of the velocity and ballistic coefficient AMSREs can be written similar to the AMSRE_p. The position, velocity and ballistic coefficient AMSREs for the various filters are listed in

Table 2. AMSREs for various filters.

Various algorithms	AMSREp(m)	AMSREv(m/s)	AMSREβ(kg/m2)
UKF	2521.684	329.911	155.735
DDF	2521.573	329.903	155.276
IUKF	1035.340	259.173	149.603
IDDF	1035.273	260.771	149.756
MLIDDF	968.746	255.916	144.953

.

Table 3. Runtimes of various filters.

Various algorithms	Runtime(s)
UKF	1.0840
DDF	0.2888
IUKF	1.9918
IDDF	0.5133
MLIDDF	1.3074

lists the runtimes of the various filters.

From Table 2, it is seen that the position AMSRE for the MLIDDF algorithm is reduced by 62% compared to the UKF and DDF, by 7% compared to IUKF and IDDF. The velocity AMSRE for the MLIDDF algorithm is reduced by 23% compared to the UKF and DDF algorithms. The ballistic coefficient AMSRE for the MLIDDF algorithm is reduced by 7% compared to the UKF and DDF algorithms. Although the ballistic coefficient AMSRE for the MLIDDF algorithm has no significant reduction, we can see that it is reduced significantly compared to the UKF, DDF, IUKF and IDDF algorithms, hence, the MLIDDF is preferred over the other filters in the light of position, velocity and ballistic coefficient AMSREs.

From Table 3, we can see the runtime of MLIDDF is less than that of the IUKF algorithm, and is more than those of UKF, DDF, and IDDF algorithms, so the accuracy of the MLIDDF algorithm is improved at the cost of an increased computational burden. Meanwhile, we can observe the AMSREs for the UKF and DDF in position, velocity and ballistic coefficient are almost identical, and the AMSREs of the IUKF and IDDF are almost the same. Therefore, on the basis of the simulation results presented in Figures 6–8 and Table 2, one can draw the conclusion that the MLIDDF yields a superior performance over the other filters.

4. Conclusions and Future Work

In this study, we provide the maximum likelihood based iterated divided difference filter which inherits the virtues of the divided difference filter and contains the iteration process in the measurement update step. The sequence obtained is guaranteed to move up the likelihood surface using the iteration termination condition based on the maximum likelihood surface. The maximum likelihood based iterated divided difference is implemented easily and is derivative-free. We apply the new filter to state estimation for a ballistic target re-entry scenario and compare its performance against the unscented Kalman filter, divided difference filter, iterated unscented Kalman filter and iterated divided difference filter with the traditional termination criteria. Simulation results demonstrate that the maximum likelihood-based iterated divided difference is much more effective than the other filters. The maximum likelihood-based iterated divided difference greatly improves the performance of state estimation and has a shorter convergence time.

Future work may focus on the applications of the maximum likelihood iteration divided difference filter to remove the outliers which is a serious deviation from the sample and caused by blink and subjective eye movement in video nystagmus signal samples of pilot candidates.