Stochastic Feedback Based Continuous-Discrete Cubature Kalman Filtering for Bearings-Only Tracking

Abstract

Bearings-only tracking only adopts measurements from angle sensors to realize target tracking, thus, the accuracy of the state prediction has a significant influence on the final results of filtering. There exist unpredictable approximation errors in the process of filtering due to state propagation, discretization, linearization or other adverse effects. The idea of online covariance adaption is proposed in this work, where the post covariance information is proved to be effective for the covariance adaption. With theoretical deduction, the relationship between the posterior covariance and the priori covariance is investigated; the priori covariance is modified online based on the feedback rule of covariance updating. The general framework integrates the continuous-discrete cubature Kalman filtering and the feedback rule of covariance updating. Numerical results illustrated that the proposed method has advantages over decreasing unpredictable errors and improving the computational accuracy and efficiency.

1. Introduction

Bearings-only tracking (BOT) [1] has universal applications in the field of navigation, especially in passive target tracking. In the BOT process, only angle parameters measured by the angle sensors (such as passive angle sensor, sonar, antenna array, sensor network, infrared sensor) are adopted to achieve the localization and tracking of the target object. The accuracy of state prediction and measurement information are two main factors that influence results of the target tracking. In recent years increasing accuracy of bearings with the development of the passive angle sensors has been witnessed. Generally, filtering methods describe the state process in a discretized form, such as extended Kalman filtering (EKF) [2], unscented Kalman filtering (UKF) [3], particle filter (PF) [4], cubature Kalman filtering (CKF) [5] and so on. The form is propitious for calculation and is convenient to implement mathematical deductions. However, the state process of a BOT system is continuous, yet the measurement process is discrete, and discretization of state process would inevitably lead to the increase of estimation errors. Thus, the continuous-discrete (CD) filtering methods are more suitable to solve a BOT problem than those traditional filtering methods. The CD filtering methods can also be applied for target tracking [6], finance [7], stochastic control [8], etc.

While state models for most target tracking algorithms are described as stochastic differential equations (SDEs), CD filtering methods are quite different from traditional methods. Thus, mathematical models for CD filtering methods are more complicated than traditional ones, while the CD methods are potentially more accurate [9]. The common form of the stochastic system can be expressed in the form of stochastic differential equation (SDE):

$$d\mathit{x}(t)=\mathit{F}(\mathit{x}(t),t)\mathit{dt}+\mathit{G}(t)d\mathit{w}(t),$$

(1)

where $\mathit{x}(t)\in {ℝ}^n$ is the n-dimensional target state vector, $\mathit{F}:{ℝ}^n\times {ℝ}^p\to {ℝ}^n$ is known as the drift function, $\mathit{G}(t)$ is the diffusion matrix, and $\mathit{w}(t)$ is a Brownian motion, which is known as a Wiener process. The alternative form is shown below:

$$\frac{d\mathit{x}(t)}{\mathit{dt}}=\mathit{F}(\mathit{x}(t),t)+\mathit{G}(t)\frac{d\mathit{w}(t)}{\mathit{dt}}$$

(2)

This form known as Langevin form in Physics is widely used, while the derivative of the Wiener term $\frac{d\mathit{w}(t)}{\mathit{dt}}$ actually does not exist formally. It should be written as Equation (1). Moreover, Equation (1) cannot be solved simply as the following integral function if the diffusion matrix is a random process related to noises.

$$\mathit{x}(t)=\mathit{x}(t_0)+{\displaystyle {\int }_{t_0}^t\mathit{F}(\mathit{x}(s),s)\mathit{ds}}+{\displaystyle {\int }_{t_0}^t\mathit{G}(\mathit{x}(s),s)d\mathit{w}(t)}$$

(3)

where the covariance (suppose it is 1-dimension) of the Wiener process is related to t, the value of the covariance will trend to be infinitely large; therefore, it makes no sense for the second item to be calculated based on the rules of Riemann--Stieltjes calculus. The SDEs problem mentioned above can be solved by Itô calculus and Stratonovich calculus [10], the difference between them can be referred to in Reference [11]. These characteristics distinguish the CD filtering methods from the discrete time domain filtering methods.

Many numerical methods are applied to solve approximately the SDEs problem in continuous time domain filtering methods. Taylor approximation and Runge-Kutta approximation are two main numerical methods which are widely used to design continuous time domain filtering methods. In Reference [12], the order 0.5 Euler-Maruyama method was proposed, the Euler scheme of SDEs was applied in the method to achieve the 0.5 rate of convergence. Similarly, the order 1.5 strong Taylor approximations and CKF are integrated in Reference [13] to propose continuous-discrete CKF (CD-CKF); however, the computation efficiency is compromised by too many derivative operations in the process of state propagating. The deterministic Runge-Kutta methods and a moment matching technique are used for continuous-time cubature Kalman filters in Reference [14], where the non-additive process noise is considered. The accurate and effective implementation of the filtering method is a key to most researches. In Reference [15], an efficient embedded Runge-Kutta pair and automatic global error control are proposed to improve performance of the complex computational procedure. Similarly, the global error control is proposed to enhance the accuracy and the robustness of the continuous-discrete extended Kalman filter (CD-EKF) in Reference [6]. The Explicit Singly Diagonally Implicit Runge-Kutta (ESDIRK) integrator with sensitivity analysis capabilities in Reference [16] was proved to be efficient to solve nonlinear continuous-discrete filtering problems for stochastic systems, in which the internal integration step is chosen by the step-size controller for Runge-Kutta method. The step-size control technology was also adopted in Reference [17], and the variable-step-size Gauss- and Lobatto-type nested implicit Runge-Kutta formulas of orders four and six are built within the EKF framework. Automatic local and global error regulation mechanisms were implemented in Reference [16,17], which are similar to the notion of accurate continuous-discrete extended Kalman filtering (ACD-EKF) proposed by Kulikov and Kulikova [6,15]. There exist the unpredictable errors due to discretization, linearization, several state predictions in measurement intervals or other causes, which are not considered in the above methods. In Reference [18], Kalman fitler (KF) and CKF are integrated to reduce the problem of the amount of calculations, low accuracy and poor convergence in discrete-time domain, and the adaptive fading factor to adjust the error covariance. In Reference [19], the covariance adaption scheme within the EKF framework was presented to solve the problem. However, the details of the state propagation in continuous-discrete framework are not fully considered, and undoubtedly, larger errors occur when EKF is applied to deal with nonlinear problems.

To alleviate the effects of unpredictable errors, the stochastic feedback of continuous discrete cubature Kalman filter (SFCD-CKF) is proposed in this paper. Compared with the CD-CKF and other CD filtering methods, SFCD-CKF adopts posterior information in the stochastic feedback framework to adapt the priori error covariance online, and online adaption decreases the costly computation during the approximation process and improves the accuracy of state estimation. The paper is structured as follows: Section 2 presents the BOT model, Section 3 briefly introduces general types of continuous-discrete filters, and Section 4 deduces covariance updating and presents stochastic feedback framework of CD-CKF. Simulations results are shown in Section 5, and Section 6 concludes the main work.

2. The Model of BOT in Continuous-Discrete Form

Suppose that the positions are set in a 2-dimensional Cartesian coordinate, the measurement information is provided by the observer with the passive angel sensor, the state of target is: ${\mathit{x}}_{}^t(t)={[x_{}^t(t),y_{}^t(t),{\dot{x}}^t(t),{\dot{y}}^t(t)]}^T$, where the position and velocity are represented by vector: ${\mathit{p}}_{}^t(t)={[x_{}^t(t),y_{}^t(t)]}^T$ and ${\mathit{v}}_{}^t(t)={[{\dot{x}}^t(t),{\dot{y}}^t(t)]}^T$ respectively. The state of the observer is: ${\mathit{x}}_{}^{\mathrm{s}}(t)={[x_{}^{\mathrm{s}}(t),y_{}^{\mathrm{s}}(t),{\dot{x}}^{\mathrm{s}}(t),{\dot{y}}^{\mathrm{s}}(t)]}^T$, and the related state is: $\mathit{x}(t)={\mathit{x}}_{}^t(t)-{\mathit{x}}_{}^{\mathrm{s}}(t)={[x(t),y(t),\dot{x}(t),\dot{y}(t)]}^T$.

The state equation of the bearings-only tracking system is given by a stochastic differential equation:

$$d\mathit{x}(t)=\mathit{f}(\mathit{x}(t),t)\mathit{dt}+\sqrt{\mathit{Q}}d\mathit{w}(t),$$

(4)

where $\mathit{Q}$ is the diffusion matrix, the other parameters are the same as Equation (1). The measurements can be described in a discrete time form:

$$z_k=h({\mathit{x}}_k,k)+{\mathsf{\omega }}_k,k=1,2,\cdots$$

(5)

$$h({\mathit{x}}_k,k)=\arctan \frac{x_k}{y_k}$$

(6)

where the measurement noise ${\mathsf{\omega }}_k$ is assumed to follow independent Gaussian distribution with zero mean and known covariance matrix ${\mathsf{\sigma }}^2$.

3. General Types of Continuous-Discrete Filters

3.1. The Taylor Approximation

The Taylor approximation and Runge-Kutta approximation are two main approaches applied to CD filters. The forms of the approximation equations used in these approaches are quite different, and details of Taylor approximation are illustrated below.

The CD-EKF [20] shows the general form of the Taylor approximation. For the time interval $(t,t+\delta )$, the state is:

$$\mathit{x}(t+\delta )=\mathit{x}(t)+\delta \mathit{f}(\mathit{x}(t),t)+\sqrt{\mathit{Q}}\mathit{\beta },$$

(7)

where $\mathit{\beta }$ is the Gaussian random variable, which is independent of the state. The expectation can be described as:

$$E[\mathit{x}(t+\delta )]=E[\mathit{x}(t)]+\delta E[\mathit{f}(\mathit{x}(t),t)].$$

(8)

Let $t=t_k$ and $t+{\delta =t}_{k+1}$, the covariance matrix is:

$$\mathrm{var}[\mathit{x}(t+\delta )]=\mathrm{var}[\mathit{x}(t)+\delta \mathit{f}(\mathit{x}(t),t)]+\delta \mathit{Q},$$

(9)

where $\delta =T$ is the measurement sampling interval. $\mathit{f}(\mathit{x}(t),t)$ is extended around the known estimate ${\widehat{\mathit{x}}}_{k|k}=E[{\mathit{x}}_k|z_{1:k}]$, which is:

$$\mathit{f}(\mathit{x}(t),t)=\mathit{f}({\widehat{\mathit{x}}}_{k|k},t)+{\mathit{f}}_{\mathrm{x}}(\mathit{x}(t)-{\widehat{\mathit{x}}}_{k|k})+R_n(\mathit{x}(t)),$$

(10)

where ${\mathit{f}}_x(\mathit{x}(t)-{\widehat{\mathit{x}}}_{k|k})$ is the Jacobian of $\mathit{f}$, and $R_n(\mathit{x}(t))$ are the high-order terms.

The predicted state estimate is:

$${\widehat{\mathit{x}}}_{k+1|k}=E[{\mathit{x}}_{k+1}|z_{1:k}]\approx {\widehat{\mathit{x}}}_{k|k}+T\mathit{f}({\widehat{\mathit{x}}}_{k+1|k},k).$$

(11)

Then the predicted covariance is:

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathrm{k}+1|\mathrm{k}}& =\mathrm{var}[{\mathit{x}}_{k+1}|z_{1:k}]\hfill \\ & \approx (\mathit{I}+T{\mathit{f}}_{\mathrm{x}}(k))P_{\mathrm{k}|\mathrm{k}}{(\mathit{I}+T{\mathit{f}}_{\mathrm{x}}(k))}^{\mathrm{T}}+T\mathit{Q}\hfill \end{array}.$$

(12)

The 1.5 order approximation is more accurate and the equation is more complex than the method shown above, but the basic form and the perception for deduction are similar. The state is:

$$\begin{array}{cc}\hfill \mathit{x}(t+\delta )=& \mathit{x}(t)+\delta \mathit{f}(\mathit{x}(t),t)\hfill \\ & +\frac{1}{2}{\delta }^2({\mathsf{\Gamma }}_0\mathit{f}(\mathit{x}(t),t)+\sqrt{\mathit{Q}}\mathit{\beta }+(\mathsf{\Gamma }\mathit{f}(\mathit{x}(t),t)\mathit{\gamma }\hfill \end{array},$$

(13)

where $(\mathit{\beta },\mathit{\gamma })$ are the pair of correlated Gaussian random variables. The noise-free process function is defined as:

$${\mathit{f}}_{\mathrm{d}}(\mathit{x}(t),t)={\mathit{x}}_k+\delta \mathit{f}(\mathit{x}(t),t)+\frac{1}{2}{\delta }^2({\mathsf{\Gamma }}_0\mathit{f}(\mathit{x}(t),t),$$

(14)

where

$${\mathsf{\Gamma }}_0=\frac{\partial }{\partial t}+{\displaystyle \sum _{i=1}^n{\mathit{f}}_i}\frac{\partial }{\partial {\mathit{x}}_i}+\frac{1}{2}{\displaystyle \sum _{j,p,q=1}^n{\sqrt{\mathit{Q}}}_{p,j}{\sqrt{\mathit{Q}}}_{q,j}}\frac{{\partial }^2}{\partial {\mathit{x}}_p\partial {\mathit{x}}_q}$$

$${\mathsf{\Gamma }}_j={\displaystyle \sum _{i=1}^n{\sqrt{\mathit{Q}}}_{i,j}}\frac{\partial }{\partial {\mathit{x}}_i}.$$

3.2. CD-CKF

The cubature rule used in CKF has the advantage in solving nonlinear filtering problems, in which the linearization is not needed in filtering. The 1.5 order Itô-type approximation CD-CKF is easily understood based on the details shown in the last section. For the interval $j\in (k,k+1)$, the predicted state estimate is:

$$\begin{array}{cc}\hfill {\widehat{\mathit{x}}}_{k|k}^1& =E[{\mathit{x}}_k^1|z_{1:k}]\hfill \\ & \approx E[{\mathit{f}}_d({\mathit{x}}_k,\mathit{kT})+\sqrt{\mathit{Q}}\mathit{\beta }+\mathsf{\Gamma }\mathit{f}({\mathit{x}}_k,\mathit{kT})\mathit{\gamma }|z_{1:k}]\hfill \end{array}.$$

(15)

For interval $j$, the error covariance matrix will also be propagated after the state propagating.

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathrm{k}|\mathrm{k}}^1\approx & {\displaystyle {\int }_{{ℝ}^n}{\mathit{f}}_{\mathrm{d}}({\mathit{x}}_k,\mathit{kT}){\mathit{f}}_{\mathrm{d}}^{\mathrm{T}}({\mathit{x}}_k,\mathit{kT})}N({\mathit{x}}_k;{\widehat{\mathit{x}}}_{k|k},{\mathit{P}}_{\mathrm{k}|\mathrm{k}})d{\mathit{x}}_k\hfill \\ & +\frac{{\delta }^3}{3}(\mathsf{\Gamma }\mathit{f}({\widehat{\mathit{x}}}_{k|k},\mathit{kT})(\mathsf{\Gamma }\mathit{f}{({\widehat{\mathit{x}}}_{k|k},\mathit{kT})}^{\mathrm{T}}\hfill \\ & +\frac{{\delta }^2}{2}[\sqrt{\mathit{Q}}(\mathsf{\Gamma }\mathit{f}{({\widehat{\mathit{x}}}_{k|k},\mathit{kT})}^{\mathrm{T}}+(\mathsf{\Gamma }\mathit{f}({\widehat{\mathit{x}}}_{k|k},\mathit{kT}){\sqrt{\mathit{Q}}}^{\mathrm{T}}]\hfill \\ & -({\widehat{\mathit{x}}}^1{}_{k|k}){({\widehat{\mathit{x}}}^1{}_{k|k})}^{\mathrm{T}}+\delta \mathit{Q}\hfill \end{array}$$

(16)

Based on the cubature rule, the predicted state estimate can be expressed as:

$${\widehat{\mathit{x}}}_{k|k}^1\approx \frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}{\mathit{X}}_{i,k|k}^{\ast (1)}},$$

(17)

$${\mathit{X}}_{i,k|k}^{\ast (1)}{=f}_{\mathrm{d}}({\widehat{\mathit{x}}}^{}{}_{k|k}+\sqrt{{\mathit{P}}_{k|k}}{\xi }_i,\mathit{kT}).$$

(18)

Then, the predicted state error covariance can be expressed as:

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathrm{k}|\mathrm{k}}^1\approx & (\Delta {\mathit{X}}_{k|k}^{\ast (1)}){(\Delta {\mathit{X}}_{k|k}^{\ast (1)})}^{\mathrm{T}}\hfill \\ & +\frac{{\delta }^3}{3}(\mathsf{\Gamma }\mathit{f}({\widehat{\mathit{x}}}_{k|k},\mathit{kT})(\mathsf{\Gamma }\mathit{f}{({\widehat{\mathit{x}}}_{k|k},\mathit{kT})}^{\mathrm{T}}\hfill \\ & +\frac{{\delta }^2}{2}[(\mathsf{\Gamma }\mathit{f}({\widehat{\mathit{x}}}_{k|k},\mathit{kT}){\mathit{Q}}^{\mathrm{T}}+\mathit{Q}(\mathsf{\Gamma }\mathit{f}{({\widehat{\mathit{x}}}_{k|k},\mathit{kT})}^{\mathrm{T}}]+\delta \mathit{Q}\hfill \end{array},$$

(19)

where

$$\Delta {\mathit{X}}_{k|k}^{\ast (1)}=\frac{1}{\sqrt{2n}}\left[\Delta {\mathit{X}}_{1,k|k}^{\ast (1)}-{\widehat{\mathit{x}}}_{k|k}^1,\cdots ,\Delta {\mathit{X}}_{2n,k|k}^{\ast (1)}-{\widehat{\mathit{x}}}_{k|k}^1\right].$$

(20)

The details of CD-CKF algorithm can be referred to Reference [13].

4. Stochastic Feedback Framework of CD-CKF

Although the accuracy of the CD-CKF algorithm is superior over CD-EKF, the unpredictable error is inevitable. In this work, the stochastic feedback framework of CD-CKF is proposed, and the key is that the posteriori information is utilized to adapt the priori error covariance online.

4.1. Covariance Adaption

First, the states, the gain, and the covariance at each measurement interval are considered.

Estimate the updated state is:

$${\widehat{\mathit{x}}}^{}{}_{k+1|k+1}={\widehat{\mathit{x}}}^{}{}_{k+1|k}+{\mathit{W}}_{k+1}{\mathit{e}}_{k+1}.$$

(21)

The continuous-discrete cubature gain is:

$${\mathit{W}}_{k+1}={\mathit{P}}_{\mathit{xz},k+1|k}{\mathit{P}}_{\mathit{zz},k+1|k}^{-1},$$

(22)

where ${\mathit{P}}_{\mathit{xz},k+1|k}$ and ${\mathit{P}}_{\mathit{zz},k+1|k}^{}$ are cross-covariance matrix and innovations covariance matrix respectively.

The error covariance matrix is:

$${\mathit{P}}_{k+1|k+1}={\mathit{P}}_{k+1|k}-{\mathit{W}}_{k+1}{\mathit{P}}_{\mathit{zz},k+1|k}{\mathit{W}}_{k+1}^{\mathrm{T}}.$$

(23)

The innovation is:

$${\mathit{e}}_{k+1}=z_{k+1}-{\widehat{z}}_{k+1|k}.$$

(24)

The predicted measurement is:

$${\widehat{z}}_{k+1|k}=\frac{1}{2n}{\displaystyle \sum _{i=1}^{2n}h({\mathit{X}}_{i,k+1|k}^{},k+1)}\approx \mathit{H}{\widehat{\mathit{x}}}^{}{}_{k+1|k},$$

(25)

where ${\mathit{X}}_{i,k+1|k}^{}$ is the cubature points, and $\mathit{H}$ is the Jacobian matrix of partial derivatives of h. The innovations covariance matrix is:

$${\mathit{P}}_{\mathit{zz},k+1|k}^{}={\mathit{Z}}_{k+1|k}{\mathit{Z}}_{k+1|k}^{\mathrm{T}}+{\mathit{R}}_{k+1},$$

(26)

$${\mathit{Z}}_{k+1|k}=\frac{1}{\sqrt{2n}}\left[Z_{1,k+1|k}^{}-{\widehat{z}}_{k+1|k}^{},\cdots ,Z_{2n,k+1|k}^{}-{\widehat{z}}_{k+1|k}^{}\right].$$

(27)

Substituting Equations (24) and (25) into Equation (26), the innovations covariance matrix can be modified as:

$${\mathit{P}}_{\mathit{zz},k+1|k}^{}\approx \mathit{H}{\mathit{P}}_{k+1|k}^{}{\mathit{H}}^{\mathrm{T}}+{\mathit{R}}_{k+1},$$

(28)

where ${\mathit{P}}_{k+1|k}^{}$ is ${\mathit{P}}_{k|k}^j$ when $j\to k+1$.

The cross-covariance matrix is:

$$\begin{array}{cc}\hfill {\mathit{P}}_{\mathit{xz},k+1|k}^{}& ={\mathit{X}}_{k+1|k}{\mathit{Z}}_{k+1|k}^{\mathrm{T}}\hfill \\ & \frac{1}{\sqrt{2n}}\left[\Delta {\mathit{X}}_{1,k|k}^{\ast (k+1)}-{\widehat{\mathit{x}}}_{k|k}^{(k+1)},\cdots ,\Delta {\mathit{X}}_{2n,k|k}^{\ast (k+1)}-{\widehat{\mathit{x}}}_{k|k}^{(k+1)}\right]·\hfill \\ & \frac{1}{\sqrt{2n}}{\left[Z_{1,k+1|k}^{}-{\widehat{z}}_{k+1|k}^{},\cdots ,Z_{2n,k+1|k}^{}-{\widehat{z}}_{k+1|k}^{}\right]}^{\mathrm{T}}\hfill \end{array}.$$

(29)

Substituting Equation (25) into Equation (29) yields:

$${\mathit{P}}_{\mathit{xz},k+1|k}^{}\approx {\mathit{P}}_{k|k}^j{\mathit{H}}^{\mathrm{T}},j=k+1.$$

(30)

Namely, ${\mathit{P}}_{\mathit{xz},k+1|k}^{}\approx {\mathit{P}}_{k+1|k}^{}{\mathit{H}}^{\mathrm{T}}$. Then, the continuous-discrete cubature gain is:

$${\mathit{W}}_{k+1}\approx {\mathit{P}}_{k+1|k}^{}{\mathit{H}}^{\mathrm{T}}{\mathit{P}}_{\mathit{zz},k+1|k}^{-1}.$$

(31)

Secondly, the covariance updated online is considered.

Here, an assumption is presented at first: The covariance ${\mathit{P}}_{j+1|j}^{}$ is nearly a constant. Actually, the assumption is not stringently valid within the whole filtering process, but it is a perfect simplification method that ensures the explicit and real time online estimation of the covariance. The method was also applied in References [21,22].

Let’s consider a maximum likelihood estimation problem.

$$\begin{array}{cc}\hfill L({\mathit{P}}_{k+1|k}^{})& =\ln p({\mathit{e}}_{k-N+1},\cdots ,{\mathit{e}}_k|{\mathit{P}}_{k+1|k}^{})\hfill \\ & ={\displaystyle \sum _{j=k-N+1}^k\ln p({\mathit{e}}_j|{\mathit{P}}_{k+1|k}^{})}\hfill \end{array},$$

(32)

where $N$ is the time window, ${\mathit{e}}_j$ can be regarded as a normally distributed random variable, so $p(\cdot )$ can be written as:

$$p({\mathit{e}}_j|{\mathit{P}}_{k+1|k}^{})=\frac{1}{\sqrt{{(2\pi )}^m\left|{\mathit{P}}_{\mathit{zz},j+1|j}^{}\right|}}\exp \left(-\frac{1}{2}{\mathit{e}}_j^T{\mathit{P}}_{\mathit{zz},j+1|j}^{}{\mathit{e}}_j\right).$$

(33)

$L({\mathit{P}}_{k+1|k}^{})$ is a scalar and ${\mathit{P}}_{k+1|k}^{}$ is a symmetric with size $n\times n$. For derivative ${\mathsf{\Omega }}_{k+1}=\partial L({\mathit{P}}_{k+1|k}^{})/\partial {\mathit{P}}_{k+1|k}^{}$, the s-th row and t-th column element of the derivative is:

$${\mathsf{\Omega }}_{k+1}^{s,t}=-\frac{1}{2}\mathrm{tr}\left\{{\displaystyle \sum _{j=k-N+1}^k\left[\mathsf{\Theta }·\mathsf{\Psi }\right]}\right\},$$

(34)

$$\mathsf{\Theta }={\mathit{P}}_{\mathit{zz},j+1|j}^{-1}-{\mathit{P}}_{\mathit{zz},j+1|j}^{-1}{\mathit{e}}_j{\mathit{e}}_j^T{\mathit{P}}_{\mathit{zz},j+1|j}^{-1},$$

(35)

$$\mathsf{\Psi }=\frac{\partial {\mathit{P}}_{\mathit{zz},j+1|j}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}},$$

(36)

where ${\mathit{P}}_{k+1|k}^{s,t}$ is the s-th row and t-th column element of ${\mathit{P}}_{k+1|k}^{}$. Then ${\widehat{\mathit{P}}}_{k+1|k}^{}$ can be yielded by setting ${\mathsf{\Omega }}_{k+1}$ to zero, i.e., ${\mathsf{\Omega }}_{k+1}^{s,t}=0$.

Because ${\mathit{R}}_{k+1}$ and H are independent of ${\mathit{P}}_{k+1|k}^{}$, with Equations (28) and (36),

$$\mathsf{\Psi }=\mathit{H}\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}{\mathit{H}}^{\mathrm{T}},$$

(37)

$$\mathrm{tr}\left\{{\displaystyle \sum _{j=k-N+1}^k\left[\mathsf{\Theta }·\mathit{H}\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}{\mathit{H}}^{\mathrm{T}}\right]}\right\}=0.$$

(38)

Then pre- and post-multiply the matrix inside $\mathrm{tr}\{·\}$ by ${\mathit{H}}^{\mathrm{T}}$ and its inverse (or the generalized inverse) can be expressed as:

$$\mathrm{tr}\left\{{\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }·\mathit{H}\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}\right]}\right\}=0,$$

(39)

because $\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}$ is a constant matrix and its s-th row and t-th column element is one while other elements are zero. Based on the multiplication rule, the t-th column of the matrix ${\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }\mathit{H}·\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}$ cannot be zero while the other column must be zero. So, the value of $\mathrm{tr}\{·\}$ is equal to the t-th diagonal element of $\{·\}$:

$${\left\{{\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }·\mathit{H}\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}\right]}\right\}}^{t,t}=0.$$

(40)

Furthermore, because the s-th row and t-th column element of $\frac{\partial {\mathit{P}}_{k+1|k}^{}}{\partial {\mathit{P}}_{k+1|k}^{s,t}}$ is one and other elements are zero to ensure elements of $\{·\}$ to be zero, we can obtain:

$${\displaystyle \sum _{j=k-N+1}^k{\left\{\left[{\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }·\mathit{H}\right]\right\}}^{t,s}}=0,$$

(41)

where t, s can be any value within (0,n), so:

$${\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }·\mathit{H}\right]}=0.$$

(42)

Let us multiply both sides of Equation (42) by ${\mathit{P}}_{j+1|j}^{}$, we get:

$${\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{P}}_{j+1|j}^{}{\mathit{H}}^{\mathrm{T}}\mathsf{\Theta }·\mathit{H}{\mathit{P}}_{j+1|j}^{}\right]}=0.$$

(43)

With Equations (31) and (35), we have:

$${\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{W}}_{j+1}\mathit{H}{\mathit{P}}_{j+1|j}^{}-{\mathit{W}}_{j+1}{\mathit{e}}_j{\mathit{e}}_j^{\mathrm{T}}{\mathit{W}}_{j+1}^{\mathrm{T}}\right]}=0.$$

(44)

Based on Equations (21) and (23), we can get:

$${\mathit{W}}_{j+1}{\mathit{e}}_j=\Delta {\widehat{\mathit{x}}}_{j+1}={\widehat{\mathit{x}}}^{}{}_{j+1|j+1}-{\widehat{\mathit{x}}}^{}{}_{j+1|j},$$

(45)

$${\displaystyle \sum _{j=k-N+1}^k\left[{\mathit{P}}_{j+1|j}^{}-{\mathit{P}}_{j+1|j+1}^{}-\Delta {\widehat{\mathit{x}}}_{j+1}\Delta {\widehat{\mathit{x}}}_{j+1}^{\mathrm{T}}\right]}=0.$$

(46)

Then:

$${\displaystyle \sum _{j=k-N+1}^k{\mathit{P}}_{j+1|j}^{}={\displaystyle \sum _{j=k-N+1}^k\left({\mathit{P}}_{j+1|j+1}^{}+\Delta {\widehat{\mathit{x}}}_{j+1}\Delta {\widehat{\mathit{x}}}_{j+1}^{\mathrm{T}}\right)}}.$$

(47)

Considering the assumption, ${\mathit{P}}_{k+1|k}^{}$ can be approximated by:

$${\widehat{\mathit{P}}}_{k+1|k}^{}=\frac{1}{N}{\displaystyle \sum _{j=k-N+1}^k{\mathit{P}}_{j+1|j}^{}}={\displaystyle \sum _{j=k-N+1}^k{\mathit{P}}_{j+1}^{\ast }},$$

(48)

where ${\mathit{P}}_{j+1}^{\ast }$ is defined as an intermediate matrix. Hence:

$${\widehat{\mathit{P}}}_{k+1|k}^{}={\widehat{\mathit{P}}}_{k|k-1}^{}+\left({\mathit{P}}_k^{\ast }-{\mathit{P}}_{k-N}^{\ast }\right).$$

(49)

4.2. Stochastic Feedback Framework

The functions, Equations (47–49), constitute a major part of the covariance adaption framework. As shown in Figure 1, the state prediction and measurement updating are completed within CD-CKF, the difference is the covariance prediction process. The post-covariance information ${\mathit{P}}_{k-N}^{\ast }$ in the memorizer is used to update the covariance ${\widehat{\mathit{P}}}_{k+1|k}^{}$ online, and it is also used to update the cubature gain ${\mathit{W}}_{k+1}$. Then the $\Delta {\widehat{\mathit{x}}}_{k+1}$ is calculated to generate the state estimation ${\widehat{\mathit{x}}}_{k+1|k+1}^{}$, and $\Delta {\widehat{\mathit{x}}}_{k+1}$ can be used for covariance generation from feedback channel.

Compared with CD-EKF and CD-CKF, the advantages of the stochastic feedback CD-CKF (SFCD-CKF) framework can be concluded as:

(1)

Post covariance information is used to decrease the influence of the unpredictable error within continuous-time domain state prediction process, and as a result, it enhances the accuracy of the filtering.

(2)

Online covariance updating can decrease the computational complexity in derivative or matrix operation compared with traditional methods.

(3)

The cubature rule improves the performance of the filtering when dealing with non-linear problems, and it also provides a more accurate innovation covariance matrix than CD-EKF and other methods.

5. Numerical Simulation

To illustrate the performance of the proposed method, the linear and non-linear continuous time domain models with discrete measurements are considered in this section. Specifically, “linear” and “non-linear” described here refer to the type of state models, because the measurement models in BOT are all non-linear in this study. Itô-1.5 order CD-CKF [13], CD-EKF [20] and continuous-discrete adaptive Kalman filter (CD-AKF, which is also called SFCD-CKF) [19] methods are compared in the simulations. The Monte Carlo simulation is set to be 200 times for each scenario.

5.1. Linear Model Tracking

The constant velocity (CV) model is used as the linear state model of BOT. The model can be rewritten as:

$${\mathit{x}}_k=\mathit{F}{x}_{k-1}+{\mathit{v}}_{k-1},$$

(50)

where $v_{k-1}$ is a zero-mean Gaussian noise with covariance matrix $\mathit{Q}$, $\mathit{F}$ is the state transition matrix. The Gaussian noise distribution can approximately describe the state error in real physical space, and Gaussian noise is convenient to be expressed and deduced mathematically. The measurements function is Equation (5). $q={10}^{-5} {\mathrm{km}}^2{/\mathrm{s}}^3$ is the intensity of process noise, and the measurement noise is $e_k~N(0,{10}^{-2})$. The time window N is 5. The observer’s initial state is ${\left[0 \mathrm{m},0 \mathrm{m},100 \mathrm{m}/\mathrm{s},50 \mathrm{m}/\mathrm{s}\right]}^T$. Target’s initial state is ${\left[1\times {10}^5 \mathrm{m},1\times {10}^5 \mathrm{m},40 \mathrm{m}/\mathrm{s},-190 \mathrm{m}/\mathrm{s}\right]}^T$.

$$\mathit{Q}=\left[\begin{array}{cccc}\raisebox{1ex}{$\Delta t^3$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& 0& \raisebox{1ex}{$\Delta t^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.& 0\\ 0& \raisebox{1ex}{$\Delta t^3$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& 0& \raisebox{1ex}{$\Delta t^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \raisebox{1ex}{$\Delta t^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.& 0& \Delta t& 0\\ 0& \raisebox{1ex}{$\Delta t^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.& 0& \Delta t\end{array}\right]·q,$$

$$\mathit{F}=\left[\begin{array}{cccc}1& 0& \Delta t& 0\\ 0& 1& 0& \Delta t\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right].$$

The root mean square error (RMSE) of position is applied to evaluate the accuracy of BOT.

$$RMSE(k)={\left(\frac{1}{N_{\mathrm{MC}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{MC}}}{\Vert {\mathit{x}}_k-{\widehat{\mathit{x}}}_k^i\Vert }_2^2}\right)}^{1/2},$$

(51)

where $N_{\mathrm{MC}}$ is the number of Monte Carlo simulations, ${\widehat{\mathit{x}}}_k^i$ is the estimation at time k for Monte Carlo simulation i, ${\mathit{x}}_k$ is the real state of the target.

As is shown in Figure 2a, the RMSE of CD-CKF decreases faster than the other three methods in the first 10 min. Because of the superior performance of CD-CKF in the nonlinear problem (actually, it is a nonlinear measurement problem in this CV BOT model), CD-CKF is more accurate for filtering in the whole BOT process. The performance of SFCD-CKF is not as good as CD-CKF at the initial moment, while RMSE of SFCD-CKF remains in decline to the end of the process, and the accuracy is higher than the other methods. CD-AKF is a CD-EKF method modified by the idea of covariance online updating, and its accuracy is higher than CD-EKF; however, when compared with CD-CKF and SFCD-CKF, the framework of EKF used in CD filtering is suboptimal, and the accuracy of CD-AKF and CD-EKF are lower comparatively. Figure 3 shows the RMSE with different sampling intervals. In general, the accuracy can be improved with smaller sampling interval, while the computational efficiency would decrease. Figure 4 shows the state estimate with these different filtering methods, where the real state is the track (or path) of the target, and the path length is 48.53 km.

Figure 2b shows the relative ratio of the computational time, and the CD-EKF with sampling interval 0.01 is the datum point. In general, the computational time will decrease with the increase of sampling interval. Because of the complex matrix operation, CD-CKF is the most time-consuming method computationally. CD-AKF and SFCD-CKF exhibit better computational efficiency comparatively, which shows that the framework of covariance online updating has an advantage in improving the computational efficiency.

5.2. Nonlinear Model Tracking

In this section, the classical nonlinear model Van der Pol oscillator is used as the state model [15].

$$d\left[\begin{array}{c}x(t)\\ y(t)\end{array}\right]=\left[\begin{array}{c}y(t)\\ \lambda (1-x^2(t))y(t)-x(t)\end{array}\right]\mathit{dt}+\left[\begin{array}{c}0\\ 1\end{array}\right]\mathit{dw}(t),$$

(52)

where $y(t)=\dot{x}(t)$ in the Van der Pol oscillator model, and the initial state is ${[1,1]}^T$. The state $\mathit{x}(t)=\mathit{x}(0)+{[t,x(t)]}^T$, where the initial position is $\mathit{x}(0)={\left[1\times {10}^5 \mathrm{m},1\times {10}^5 \mathrm{m}\right]}^T$. The parameters are λ = 0.3 and $w(t)~N(0,{10}^{-5})$. The measurement is defined by Equation (5).

Figure 5a shows the different RMSEs of four methods when dealing with a nonlinear filtering problem. CD-CKF and SFCD-CKF show the advantage for nonlinear filtering. In the last stage of the target tracking process, the accuracy of SFCD-CKF is slightly higher than CD-CKF. As we can see from Figure 2a and Figure 5a, the performance of SFCD-CKF is superior over CD-CKF only when the error approaches zero. This indicates that the online covariance updating has an advantage in state estimation with micro variation of the covariance because the approach is sensitive to unpredictable errors in this moment than in the moments when the covariance changes rapidly. The computational efficiency of SFCD-CKF and CD-AKF are demonstrated to be better than the other two methods in Figure 5b. Figure 6 shows the RMSE with different sampling intervals. The accuracy of CD-EKF and CD-AKF are improved rapidly along with the decrease of the sampling intervals, because the error of EKF framework used by CD-EKF and CD-AKF is large when dealing with the nonlinear filtering problem, and the impact of sampling interval on estimation accuracy is significant. Figure 7 shows the state estimate with these different filtering methods, where the real state is the track (or path) of the target. The path length is 133.01 km. With superior computational efficiency, SFCD-CKF generally has an advantage in solving unpredictable error problems.

6. Conclusions

Considering the unpredictable error within the process of CD filtering, the framework of online covariance updating is proposed. The key is to use the post information to modify the state estimation. The post covariance information is proved to be valuable for covariance updating by theoretical deduction. The proposed framework integrates CD-CKF with the covariance updating approach, which enhance the performance for the method when implementing the nonlinear filtering, the efficient computation and the correction of unpredictable error. While the approach does not work well when the covariance changes significantly, the adaptive approaches which combine CD-CKF with SFCD-CKF will be the focus of our future research.