Adaptive and robust fractional gain based interpolatory cubature Kalman filter

Abstract

In this study, we put forward the robust fractional gain based interpolatory cubature Kalman filter (FGBICKF) and the adaptive FGBICKF (AFGBICKF) for the development of the state estimators for stochastic nonlinear dynamics system. FGBICKF introduces a fractional gain to interpolatory cubature Kalman filter to increase the robustness of state estimation. AFGBICKF is developed to enhance the state estimation adaptive to stochastic nonlinear dynamics system with unknown process noise covariance through recursive estimation. The simulations on re-entry target tracking system have shown that the performance of FGBICKF is superior to that of cubature Kalman filter and interpolatory cubature Kalman filter, and standard deviation of FGBICKF is closer to posterior Cramér-Rao lower bound. Moreover, our simulations have also demonstrated that AFGBICKF remains stable even when the initial process noise covariance increase, proving its adaptiveness, robustness, and effectiveness on state estimation.

Introduction

Bayesian filtering (BF) has been intensively researched in various applications such as communication, state estimation, and signal processing. ¹ Extended Kalman filter (EKF) is a common state estimation algorithms used in the target tracking problem. ² However, the implementation of EKF requires computing Jacobian derivative, so EKF may lead to large errors and even divergence in some systems. Then free-derivative Kalman filters were proposed, including unscented Kalman filter (UKF), ³ cubature Kalman Filters (CKF), ⁴ improved CKF, ⁵ robust and adaptive CKF, ⁶ interpolatory CKF (ICKF), ⁷ the improved central difference Kalman filter, ⁸ etc. In the approaches above, the filter gain is achieved by calculating the state covariance and cross covariance between state and measurement. In the course of target tracking, the filter gain lags behind the true state of the target when the target maneuvers, which causes severe tracking errors. So researchers have studied the gain based robust Bayesian filter, such as an optimal control approach to designing constant gain filters, ⁹ modified gain EKF,^10,11 adaptive-gain tracking filters based on minimization of the innovation variance, ¹² UKF using modified filter gain, ¹³ analysis of the characteristic of the filter gain in cubature Kalman filter for 1D chaotic maps, ¹⁴ and distributed adaptive high-gain EKF. ¹⁵

Generally, the prior knowledge of process noise covariance increases the performance of estimation. In practice, however, the prior knowledge of process noise covariance is usually unknown, resulting in model mismatching. Adaptive filtering is an effective method to solve the model mismatching problem. For instance, adaptive Bayesian filters were proposed for unknown process noise covariance, including adaptive Kalman Filtering for dynamic system with outliers,^16–19 variational Bayesian based Kalman filtering, ²⁰ and prior probability statistics based robust estimation algorithm. ²¹ These algorithms were developed for the linear and time-invariant systems, assuming prior probability statistics (such as sampling distributions). Some Kalman or extended Kalman based filters were developed such as ELM based adaptive Kalman filter, ²² Kalman filter for dynamic state estimation based on adaptive adjustment of noise covariance, ²³ adaptively estimate Q and R based on innovation and residual extended Kalman filter, and sample-based adaptive Kalman filtering for accurate camera pose tracking. ²⁴ For specific applications, the adaptive square-root sigma-point Kalman ²⁵ and adaptive embedded CKF, ²⁶ improved CKF for spacecraft attitude estimation ⁵ were proposed. Then the robust Huber-Based Cubature Kalman Filter for GPS Navigation Processing ²⁷ was put forward for the non-Gaussian noise. Despite the adaptive filters mentioned above can well address the model mismatching to some extent, they still have limitations such as non-positive covariance matrices and heavy computational load.

Since fractional Kalman filter (FKF) was developed for tracking vehicles, ²⁸ quite a few fractional KFs have been proposed recently as fractional calculus gives more accurate results for system analysis. These filters include FCKF for fractional-order nonlinear stochastic systems, ²⁹ fractional central difference Kalman filter, ³⁰ fractional ICKF,^31,32 innovation-based fractional adaptive KF, ³³ fractional feedback KF for vehicle tracking in video ³⁴ and CKF for continuous-time nonlinear fractional-order systems. ³⁵

It can be seen from the filters mentioned above that ICKF is practical for state estimation. However, the filter gain in the ICKF is obtained by calculating the covariance of state and cross-covariance between state and measurement, which causes the filter gain to drop behind the target state in rapid change situations. Especially, when the target highly maneuvers, which causes severe tracking error, the tracking performance of the ICKF becomes worse and may diverge. To avoid such divergence, in this paper we put forward the fractional gain based ICKF (FGBICKF), in which the filter gain uses fractional derivative and the gain of the present state depends upon the previous ones. The gain will never be too large, and the FGBICKF performs better even when the target highly maneuvers. Thus, the proposed filter improves the state estimation performance by modifying the filter gain using factional calculus. Moreover, we take a further step to propose adaptive FGBICKF (AFGBICKF) with recursive estimation of unknown process noise covariance. The simulations on state estimation for re-entry ballistic target (RBT) tracking system have demonstrated the effectiveness and robustness of our proposed filters. The main contributions are summarized in the following:

• We propose FGBICKF, which incorporates fractional derivative of previous filter gains to estimate state of nonlinear systems under Gaussian noise. Meanwhile, we analyze the errors of FGBICKF with posterior Cramér-Rao lower bound (PCRLB). ³⁶

• We propose AFGBICKF to enhance FGBICKF for estimation under unknown process noise covariance.

• We conduct the simulation on RBT tracking with FGBICKF and AFGBICKF, and we make errors analysis of FGBICKF using PCLRB. The results show that FGBICKF outperforms CKF and ICKF.

• We analyze the influence of process noise covariance on the performance of AFGBICKF. The simulation results demonstrate AFGBICKF’s adaptiveness, robustness, and effectiveness on RBT tracking.

The remainder of this paper is organized as follows. We review some preliminaries on the Grünwald-Letnikov (G-L) fractional difference, Bayesian filtering, and interpolatory cubature rule (ICR) in Section 2. The main algorithms are derived in Section 3, where we first developed FGBICKF by introducing a fractional gain to ICKF, and then put forward AFGBICKF with adaption to unknown process noise covariance using recursive method. In Section 4, we apply the proposed filters to RBT tracking and show the performance of FGBICKF with analysis of its PCRLB. Meanwhile, we present the performance of AFGBICKF and analyze the influence of various process noise covariance on AFGBICKF. Concluding remarks on the results are drawn in Section 5.

Preliminaries

G-L fractional difference

The G-L fractional difference concept can be defined as follows:

$${\mathrm{\Delta }}^{\alpha }{x}_k=\frac{1}{h^{\alpha }}\sum _{j=0}^k{(-1)}^i\left(\begin{array}{c}\alpha \hfill \\ \multicolumn{1}{c}{j}\end{array}\right)x_{k-j}$$

(1)

where $\mathrm{\Delta }$ is the operator of fractional order system, and $\alpha \in \mathbb{R}$ ( $\mathbb{R}$ is the set of real numbers) is fractional difference order. And $h$ (it is considered to be unity in the paper) and $k$ are the sampling interval and the sampling number, respectively. The coefficient $\left(\begin{array}{c}\alpha \hfill \\ \multicolumn{1}{c}{j}\end{array}\right)$ can be calculated as:

$$\left(\begin{array}{c}\alpha \\ j\end{array}\right)=\{\begin{array}{cc}1& ifj=0\\ \frac{\alpha (\alpha -1)\cdots (\alpha -j+1)}{j!}& ifj>0\end{array}$$

(2)

Equation (2) is the discrete equivalent of derivative when $\alpha$ is greater than zero.

Bayesian filtering

We consider the stochastic nonlinear dynamics system (SNDS) with additive Gaussian noise, which is modelled using the following state and the measurement equation:

$$x_k=\mathrm{f}(x_{k-1})+w_{k-1}$$

(3)

$$l_k=\mathrm{h}(x_k)+v_k$$

(4)

where $x_k\in {\mathbb{R}}^{n_{\mu }}$ and $l_k\in {\mathbb{R}}^{m_l}$ are the state and measurement vector, respectively, $\mathrm{f}$ and $\mathrm{h}$ are known nonlinear functions; $w_{k-1}~N(0,Q_{k-1})$ is Gaussian process noise, and $v_k$ is measurement noise with zero mean and $R_k$ , respectively, $w_{k-1}$ and $v_k$ are mutually uncorrelated noises. $L_{1:k}=\{l_1,l_2,\cdots ,l_k\}$ is measurement data set obtained by one or more sensors from time step 1 to k.

In the Bayesian filter under Gaussian noise, the complete statistical description of the state ( $x_k$ ) can be obtained by the posterior density of the state. ${\hat{\mu }}_{k-1}$ and $P_{k-1}$ are denoted as the state estimation and covariance at time $k-1$ , ${\hat{\mu }}_k$ , and $P_k$ as the state estimation and covariance at time $k$ . When a new measurement ( $l_k$ ) at time $k$ is received, the posterior density $p(x_k|L_k)=N(x_k;{\hat{\mu }}_k,P_k)$ from the posterior density $p(x_{k-1}|L_{k-1})=N(x_{k-1};{\hat{\mu }}_{k-1},P_{k-1})$ of the state at time $k-1$ is obtained in two steps:

• (1) Time update, which involves computing the predictive state ${\overline{\mu }}_k$ and covariance ${\overline{P}}_k$ :

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{\mu }}_k=E[x_k|L_{k-1}]}\\ =E[\mathrm{f}(x_{k-1})+w_{k-1}|L_{k-1}]=E[\mathrm{f}(x_{k-1})|L_{k-1}]\hfill \\ ={\int }_{{\mathbb{R}}^{n_x}}\mathrm{f}(x_{k-1})p(x_{k-1}|L_{k-1})d{x}_{k-1}\hfill \\ ={\int }_{{\mathbb{R}}^{n_x}}\mathrm{f}(x_{k-1})\times N(x_{k-1};{\hat{\mu }}_{k-1},P_{k-1})d{x}_{k-1}\hfill \end{array}$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{P}}_k=E[(x_k-{\overline{\mu }}_k){(x_k-{\overline{\mu }}_k)}^T|L_{1:k}]}\\ ={\int }_{{\mathbb{R}}^{n_x}}\mathrm{f}(x_{k-1})\mathrm{f}{(x_{k-1})}^T\times N(x_{k-1};{\hat{\mu }}_{k-1},P_{k-1})d{x}_{k-1}\hfill \\ -{\overline{\mu }}_k{\overline{\mu }}_k^T+Q_{k-1}\hfill \end{array}$$

(6)

• (2) Measurement update, which involves computing the measurement prediction ${\overline{l}}_k$ , innovation covariance ${\overline{P}}_{\mathrm{ll},k}$ , cross covariance $P_{\mathrm{xl},k}$ based on predictive posterior density $x_k~N(x_k;{\overline{\mu }}_k,{\overline{P}}_k)$ obtained in the time update.

$${\overline{l}}_k={\int }_{{\mathbb{R}}^{n_x}}\mathrm{h}(x_{k-1})\times N(x_k;{\overline{\mu }}_k,{\overline{P}}_k)d{x}_k$$

(7)

$${\overline{P}}_{\mathrm{ll},k}={\int }_{{\mathbb{R}}^{n_x}}\mathrm{h}(x_{k-1})\mathrm{h}{(x_{k-1})}^T\times N(x_k;{\overline{\mu }}_k,{\overline{P}}_k)d{x}_k-{\overline{l}}_k{\overline{l}}_k^T+R_k$$

(8)

$$P_{\mathrm{xl},k}={\int }_{{\mathbb{R}}^{n_x}}{x}_k{\mathrm{h}}^T(x_k)N(x_k;{\overline{\mu }}_k,{\overline{P}}_k)-{\overline{\mu }}_k{\overline{l}}_k$$

(9)

The Kalman gain $G_k$ , estimated state ${\hat{\mu }}_k$ , and covariance $P_k$ at $k$ time instant are calculated as:

$$G_k=P_{\mathrm{xl},k}{\overline{P}}_{\mathrm{ll},k}^{-1}$$

(10)

$${\hat{\mu }}_k={\overline{\mu }}_k+G_k(l_k-{\overline{l}}_k)$$

(11)

$$P_k={\overline{P}}_k-G_k{\overline{P}}_{\mathrm{ll},k}{G}_k^T$$

(12)

Interpolatory Cubature rule

As can be seen in equations (5)–(9), Gaussian filter can be represented as weighted Gaussian integral under Bayesian framework. The product of a nonlinear function $\tau (x)$ and a Gaussian probability density function (PDF) $N(x;0,\mathbf{I})$ ( $\mathbf{I}$ is identity covariance) is described as:

$$\mathrm{Integral}[\tau ]=\int \tau (x)N(x;0,\mathbf{I})\mathrm{dx}$$

(13)

where $\mathrm{Integral}[\tau ]$ is an integration and $\tau (x)$ is an arbitrary non-linear function. $\mathrm{Integral}[\tau ]$ can be approximated by ${\varphi }^{(m,n)}(\tau )$ , which is a $2m+1$ th-degree fully symmetric interpolatory cubature rule (ICR) for a n-dimensional Gaussian weighted integral ³⁷ :

$$\mathrm{Integral}[\tau ]\approx {\varphi }^{(m,n)}[\tau ]=\sum _{p\in P^{(m,n)}}{\phi }_p^{(m,n)}\tau [\lambda ]$$

(14)

Here, $P^{(m,n)}=\{(p_1,\cdots ,p_n)|m\ge p_1\ge \cdots \ge p_n\ge 0,\left|p\right|\le m\}$ , which denotes a set of all distinct $n$ -partitions of the integers $\{0,1,\cdots ,m\}$ , $p\in \{0,1,\cdots ,m\}$ , and $|p|=\sum _{i=1}^n{p}_i$ . $\lambda$ is defined as a generator composed by $[{\lambda }_{p_1},{\lambda }_{p_2},\cdots ,{\lambda }_{p_n}]$ , ${\lambda }_0=0$ , and ${\lambda }_{p_i}\ge 0$ . The fully symmetric sum $\tau [\lambda ]$ is defined as

$$\tau [\lambda ]=\sum _{q\in {\Pi }_p}\sum _s\tau [s_1{\lambda }_{q_1},s_2{\lambda }_{q_2},\cdots ,s_n{\lambda }_{q_n}]$$

(15)

where ${\Pi }_p$ denotes all distinct permutations of $p$ and the inner sum is taken over all of the sign combinations that occur when $s_i=\pm 1$ for those values of $i$ where ${\lambda }_{q_i}\ne 0$ . The weight ${\phi }_p^{(m,n)}$ of generator $[\lambda ]$ is given by

$${\phi }_p^{(m,n)}=2^{-K}\sum _{|k|\le m-|p|}\underset{i=1}{\overset{n}{\Pi }}\frac{a_{k_i+p_i}}{\underset{j=0,\ne p_i}{\overset{k_i+p_i}{\Pi }}({\lambda }_{p_i}^2-{\lambda }_j^2)}$$

(16)

where $K$ is the number of non-zero entries in $p$ and $a_0=1$ , and $a_i$ is derived as following:

$$a_i=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }{e}^{-x^2/2}\underset{j=0}{\overset{i-1}{\Pi }}(x^2-{\lambda }_j^2)\mathrm{dx}(i>0)$$

(17)

The arbitrary degree ICR in equation (14) can be used to numerically compute the Gaussian weighted integrals in Gaussian filters. In this paper, the third-degree ICR ( $m=1$ ) is used to develop the proposed filters to balance computational load and accuracy. The third-degree ICR corresponds to $p\in \{0,1\}$ and $|p|=\sum _{i=1}^n{p}_i$ , $|p|\le 1$ (i.e. $|p|=0$ or $|p|=1$ ). When $|p|=0$ and $\left|p\right|=1$ , the basic interpolatory cubature points (ICPs) ${\xi }_j$ and the weights ${\varpi }_j$ can be calculated from (16) as

$$\begin{array}{l}{\xi }_j=\{\begin{array}{l}[0]j=1\\ {\lambda }_1{e}_{j}j=2,\cdots ,n+1\\ -{\lambda }_1{e}_{j}j=n+2,\cdots ,r=2n+1\end{array}\\ {\varpi }_j=\{\begin{array}{c}1-\frac{n}{{\lambda }_1^2}j=1\\ \frac{1}{2{\lambda }_1^2}j=2,\cdots ,r=2n+1\end{array}\end{array}$$

(18)

here, $e_i$ denotes the ith column of a unit matrix.

Using the equation (18), $\mathrm{Integral}[\tau ]$ can be calculated as:

$$\mathrm{Integral}[\tau ]\approx \sum _{j=1}^r{\varpi }_j\tau ({\xi }_j)$$

(19)

Further, we expressed the product of a nonlinear function $\tau (x)$ and a Gaussian PDF $N(x;\hat{\mu },P)$ as follows ³⁸ :

$$\mathrm{Integra}{l}_N[\tau ]=\int \tau (x)N(x;\hat{\mu },P)\mathrm{dx}$$

(20)

Using equation (19) and $P=S{S}^T$ , equation (20) can be approximated as ⁷ :

$$\mathrm{Integra}{l}_N[\tau ]\approx \sum _{j=1}^r{\varpi }_j\tau (S{\xi }_j+\hat{\mu })$$

(21)

Proposed methods

FGBICKF

We assume that the state estimate ${\hat{\mu }}_{k-1}$ and its corresponding covariance $P_{k-1}$ have been obtained at the time step $k-1$ . First, we factorize the covariance $P_{k-1}=S_{k-1}{S}_{k-1}^T$ , evaluate the ICPs using ${\xi }_j$ defined in equation (18) and propagate the ICPs through nonlinear state equation:

$$X_{j,k-1}=S_{k-1}{\xi }_j+{\hat{\mu }}_{k-1}$$

(22)

$$X_{j,k}^{*}=\mathrm{f}(X_{j,k-1})$$

(23)

We obtained the state prediction ${\overline{\mu }}_k$ and the prediction error covariance ${\overline{P}}_k$ conditioned on measurements $L_{1:k-1}$ by using the equation (21), $E[w_k]=0$ and $E(w_k{w}_k^T)=Q_k$ ,

$${\overline{\mu }}_k=E\left[x_k|L_{1:k-1}\right]\approx {\sum }_{j=1}^r{\varpi }_j{X}_{j,k}^{*}$$

(24)

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{P}}_k=E\left[({\overline{\mu }}_k-x_k){(({\overline{\mu }}_k-x_k))}^T|L_{1:k-1}\right]}\\ =E\left[(\mathrm{f}(x_{k-1})-\mathrm{f}({\hat{\mu }}_{k-1})){(\mathrm{f}(x_{k-1})-\mathrm{f}({\hat{\mu }}_{k-1}))}^T|L_{1:k-1}\right]\hfill \\ +E(w_k{w}_k^T)\hfill \\ \approx \sum _{j=1}^r{\varpi }_i[(X_{j,k}^{*}-{\overline{\mu }}_k){(X_{j,k}^{*}-{\overline{\mu }}_k)}^T]+Q_{k-1}\hfill \end{array}$$

(25)

Then, factorize the predicted covariance ${\overline{P}}_k={\overline{S}}_k{\overline{S}}_k^T$ , calculate the predicted cubature points and propagate ICPs as

$${\mathrm{\Upsilon }}_{i,k}={\overline{S}}_k{\xi }_j+{\overline{\mu }}_k$$

(26)

$${\mathrm{\Upsilon }}_{i,k}^{*}=\mathrm{h}({\mathrm{\Upsilon }}_{i,k})$$

(27)

Evaluate the predicted measurement, the cross-covariance and innovation covariance as the following:

$${\overline{l}}_k={\sum }_{i=1}^r{\omega }_i{\mathrm{\Upsilon }}_{i,k}^{*}$$

(28)

$$\begin{array}{c}\multicolumn{1}{c}{P_{\mathrm{xl},k}=E[(x_k-{\overline{x}}_k)(l_k-(h({\overline{x}}_k)))]}\\ \approx \sum _{i=1}^r{\varpi }_i{\mathrm{\Upsilon }}_{i,k}{\mathrm{\Upsilon }}_{i,k}^{*}-{\overline{x}}_k{\overline{l}}_k^T\hfill \end{array}$$

(29)

$$\begin{array}{c}\multicolumn{1}{c}{P_{\mathrm{ll},k}=E[(l_k-h({\overline{x}}_k))(l_k-{(h({\overline{x}}_k))}^T]}\\ \approx {\sum }_{i=1}^r{\varpi }_i{\mathrm{\Upsilon }}_{i,k}^{*}{\mathrm{\Upsilon }}_{i,k}^{*T}-{\overline{l}}_k{\overline{l}}_k^T+R_k\hfill \end{array}$$

(30)

In the Bayesian filtering, the Kalman gain ( $G_k$ ) is calculated as:

$$G_k=P_{\mathrm{xl},k}{P}_{\mathrm{ll},k}^{-1}$$

(31)

The filter gain is obtained by calculating the covariance of state and cross-covariance between state and measurement. This may worsen the performance of filters when the target dramatically changes its motion. With a high gain, the filter gives more weight to the measurements and thus follows the target more closely. With a low gain, the filter depends on the model predictions more closely and the state estimation accuracy decreases. To improve the performance of filters, we have proposed the fractional gain which the filter gain of the present state depends upon the previous ones. The filter gain value will never be too large so that the proposed filter is more robust to the change of the target’s motion.

Now using the fractional derivative, we define the new filter gain named fractional gain ( $G_{\mathrm{new}}$ ) as:

$$G_{\mathrm{new}}=G_k-{\sum }_{j=1}^k{(-1)}^j{\mathrm{\Upsilon }}_j{G}_{k-j}$$

(32)

where ${\mathrm{\Upsilon }}_j=\mathrm{diag}[\left(\begin{array}{c}{\alpha }_1\hfill \\ \multicolumn{1}{c}{j}\end{array}\right),\cdots ,\left(\begin{array}{c}{\alpha }_{m_l}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)]$ is fractional order. $G_{\mathrm{new}}$ includes previous filter gains using fractional derivative, that is, $G_{\mathrm{new}}$ is related to all of the previous gains ( $1,2,...,k$ time instant) and it means that $G_{\mathrm{new}}$ is with k-length memory size.

So, the state estimation and corresponding covariance are

$${\hat{x}}_k={\overline{x}}_k+G_{\mathrm{new}}(l_k-{\overline{l}}_k)$$

(33)

$$P_k={\overline{P}}_k-G_{\mathrm{new}}{P}_{\mathrm{ll},k}{G}_{\mathrm{new}}^T$$

(34)

Proof of the fractional gain ( $G_{\mathrm{new}}$ ) in equation (32).

The predicted state estimation ${\overline{\mu }}_k$ is evaluated as

$${\overline{\mu }}_k=\mathrm{f}({\hat{\mu }}_{k-1})$$

(35)

where ${\hat{\mu }}_{k-1}$ is posteriori estimated state at k  − 1 instant time.

The estimated state ${\hat{\mu }}_k$ with fractional gain is given as

$${\hat{\mu }}_k={\overline{\mu }}_k+(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-h({\overline{\mu }}_k))$$

(36)

Here, $G_g$ is gain variable and ${\mathrm{\Delta }}^{\alpha }{G}_k$ is the fractional derivative of the previous filter gain, and it is denoted as:

$${\mathrm{\Delta }}^{\alpha }{G}_k=\sum _{j=1}^k{(-1)}^j{\mathrm{\Upsilon }}_j{G}_{k-j}$$

(37)

The covariance ( $P_k$ ) of state is defined as:

$$\begin{array}{l}{P}_k=E\left\{(x-{\widehat{\mu }}_k){(x-{\widehat{\mu }}_k)}^T\right\}\\ =E\{(x-{\overline{\mu }}_k-(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-\mathrm{h}({\overline{\mu }}_k)))\\ {(x-{\overline{\mu }}_k-(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-\mathrm{h}({\overline{\mu }}_k)))}^T\}\end{array}$$

(38)

To prove the equation (32), the minimum of a posteriori error covariance has to be obtained in the following:

$$\begin{array}{c}\multicolumn{1}{c}{G_k=\underset{G_g}{\arg min}E\{(x-{\overline{\mu }}_k-(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-\mathrm{h}({\overline{\mu }}_k)))}\\ {(x-{\overline{\mu }}_k-(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-\mathrm{h}({\overline{\mu }}_k)))}^T\}\hfill \end{array}$$

(39)

To obtain $G_k$ , we solve the following equation in which the left-hand side is the first derivative of function in equation (39):

$$E[(x-{\overline{\mu }}_k-(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)(l_k-\mathrm{h}({\overline{\mu }}_k)){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]=0$$

(40)

Reformulating equation (40), we have:

$$\begin{array}{c}E[(x-{\overline{\mu }}_k){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]\hfill \\ -(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)E[(l_k-\mathrm{h}({\overline{\mu }}_k)){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]=0\hfill \end{array}$$

(41)

This yields:

$$\begin{array}{c}(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)E[(l_k-\mathrm{h}({\overline{\mu }}_k)){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]\hfill \\ =E[(x-{\overline{\mu }}_k){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]\hfill \end{array}$$

(42)

Since $P_{\mathrm{xl},k}$ and $P_{\mathrm{ll},k}$ are defined as

$$P_{\mathrm{xl},k}=E[(x_k-{\overline{\mu }}_k){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]$$

(43)

$$P_{\mathrm{ll},k}=E[(l_k-\mathrm{h}({\overline{\mu }}_k)){(l_k-\mathrm{h}({\overline{\mu }}_k))}^T]$$

(44)

Equation (41) is transformed into

$$(G_g+{\mathrm{\Delta }}^{\alpha }{G}_k)P_{\mathrm{ll},k}=P_{\mathrm{xl},k}$$

(45)

With equation (35), we have

$$G_{\mathrm{new}}=P_{\mathrm{xl},k}{P}_{\mathrm{ll},k}^{-1}-{\mathrm{\Delta }}^{\alpha }{G}_k=G_k-{\mathrm{\Delta }}^{\alpha }{G}_k$$

(46)

So $G_{\mathrm{new}}$ in equation (32) is then derived.

We summarize FGBICKF in Algorithm 1 below.

Algorithm 1.

FGBICKF algorithm.

Given the state estimates ${\hat{\mu }}_0$ and its associated error covariance $P_0$ at time $k=0$ , the state estimation procedure can be recursively implemented as follows. Initialize parameters: ${\hat{\mu }}_0$ , $P_0$ For $k=1,2,\cdots$ Step 1. Time update Calculate the state prediction ${\overline{\mu }}_k$ and covariance ${\overline{P}}_k$ using equations (24) and (25). Step 2. Measurement update Step 2.1 Calculate the predicted measurement, the cross-covariance, and innovation covariance using equations (28)–(30) Step 2.2 Compute the fractional gain $G_{\mathrm{new}}=G_k-\sum _{j=1}^k{(-1)}^j{\mathrm{\Upsilon }}_j{G}_{k-j}$ Step 2.3 Evaluate the state estimate ${\hat{\mu }}_k$ and $P_k$ covariance using equations (33) and (34). End.

Computational complexity

To analyze the numerical complexity of FGBICKF using floating-point operations (flops), we define basic arithmetic operations such as addition, subtraction, multiplication, division, comparison, or square root as one flop.

The number of flops for vector-vector operations, matrix-vector product, and matrix-matrix product are explained in Ref. ³⁰ Table 1 lists the specific flops of each step of the FGBICKF. In applications, the memory size (k) is too large as the filters is evaluated recursively, so we select the fixed constant value as memory size. In the paper, we select $M_L$ as the memory size (as in equation (32)) under the guarantee of filters’ accuracy. $F(n_{\mu })$ and $H(n_{\mu },m_l)$ are assumed to be the required flops of two nonlinear functions $\mathrm{f}(x_{k-1})$ and $\mathrm{h}(x_k)$ , respectively. Their exact computational complexity is significant but difficult to evaluate. For FGBICKF we have the total complexity:

$$\begin{array}{c}\multicolumn{1}{c}{T_{\mathrm{FGBICKF}}=6n_{\mu }^3+16n_{\mu }^2+6n_{\mu }^2{m}_l+(2F(n_{\mu })}\\ +2H(n_{\mu },m_l)+7)n_{\mu }\hfill \\ +(5+M_L+M_L{m}_l)n_{\mu }+10m_l^2{n}_{\mu }\hfill \\ +F(n_{\mu })+H(n_{\mu },m_l)+5m_l^2+2m_l^3\hfill \end{array}$$

(47)

Table 1.

Computational complexity of each step.

Step	Flops
${\overline{\mu }}_k$	$2n_{\mu }^3+7n_{\mu }^2+2(F(n_{\mu })+1)n_{\mu }+F(n_{\mu })$
${\overline{P}}_k$	$2n_{\mu }^3+2n_{\mu }^2+2n_{\mu }$
${\overline{l}}_k$	$2n_{\mu }^3+7n_{\mu }^2+2(H(n_{\mu },m_l)+1)n_{\mu }+H(n_{\mu },m_l)$
$P_{\mu l,k}$	$4n_{\mu }^2{m}_l+4n_{\mu }{m}_l$
$P_{\mathrm{ll},k}$	$4n_{\mu }{m}_l^2+5m_l^2$
$G_{\mathrm{new}}$	$M_L{n}_{\mu }{m}_l+m_l^3+2n_{\mu }{m}_l^2$
${\hat{\mu }}_k$	$(2+M_L)n_{\mu }{m}_l+2n_{\mu }{m}_l^2+n_{\mu }+m_l^3$
$P_k$	$2m_l{n}_{\mu }^2+(2m_l^2-m_l)n_{\mu }$

The numerical complexity of FGBICKF is $max\{O(n_{\mu }^3),O(m_l^3),O(n_{\mu }F(n_{\mu })),O(n_{\mu }H(m_l,n_{\mu }))\}$ .

Adaptive FGBICKF (AFGBICKF)

In pactice, the covariance of the process noise is not usually available as described in equation (3), so we develop a recursive estimation strategy based on covariance matching principle. ³⁹

Rewrite equation (25):

$${\overline{P}}_k\approx P_{\mathrm{xx},k}+Q_{k-1}$$

(48)

Here, $P_{\mathrm{xx},k}=\sum _{j=1}^r{\varpi }_i[(X_{i,k}^{*}-{\overline{\mu }}_k){(X_{i,k}^{*}-{\overline{\mu }}_k)}^T]$ .

Given one-step prediction of state ${\overline{\mu }}_k$ at time $k-1$ and the estimate state ${\hat{\mu }}_k$ at time $k$ , the residual between them can be represented by:

$${\zeta }_k={\hat{\mu }}_k-{\overline{\mu }}_k$$

(49)

Given the residual data from time $k-M_Q+1$ ( $M_Q$ is an adjustable parameter) to time $k$ , the mean and the covariance of ${\zeta }_k$ can be estimated by

$${\overline{\zeta }}_k=\frac{1}{M_Q}\sum _{j=k-M_Q+1}^k{\zeta }_j$$

(50)

$$\sum _{{\zeta }_k}=\frac{1}{M_Q-1}\sum _{j=k-M_Q+1}^k({\zeta }_j-{\overline{\zeta }}_k){({\zeta }_j-{\overline{\zeta }}_k)}^T$$

(51)

Provided that the covariance of process noise stays constant, the expectation of equation (51) is given by:

$$E\left({\sum }_{{\overline{\zeta }}_k}\right)=\frac{1}{M_Q}\sum _{j=k-M_Q+1}^k(P_{xx,j}-P_j)+Q_k$$

(52)

By combining equations (51) and (52), the covariance of process noise can be approximately calculated by

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{Q}}_k=\frac{1}{M_Q-1}\sum _{j=k-M_Q+1}^k({\zeta }_i-{\overline{\zeta }}_k){({\zeta }_i-{\overline{\zeta }}_k)}^T}\\ -\frac{1}{M_Q}\sum _{j=k-N+1}^k(P_{\mathrm{xx},j}-P_j)\hfill \end{array}$$

(53)

Similarly, $Q_{k-1}$ can be computed as

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{Q}}_{k-1}=\frac{1}{M_Q-2}\sum _{j=k-M_Q+1}^{k-1}({\zeta }_i-{\overline{\zeta }}_k){({\zeta }_i-{\overline{\zeta }}_k)}^T}\\ -\frac{1}{M_Q-1}\sum _{j=k-M_Q+1}^k(P_{\mathrm{xx},j}-P_j)\hfill \end{array}$$

(54)

Through simple mathematical manipulations, equation (53) can be rewritten as the following form

$${\hat{Q}}_k=\frac{M_Q-1}{M_Q}{W}_k+\mathrm{\Delta }{Q}_k$$

(55)

where

$$\begin{array}{c}{W}_k=\frac{M_Q}{{(M_Q-1)}^2}\sum _{j=k-M_Q+1}^{k-1}\hfill \\ ({\zeta }_i-{\overline{\zeta }}_k){({\zeta }_i-{\overline{\zeta }}_k)}^T-\frac{1}{M_Q-1}\sum _{j=k-M_Q+1}^k(P_{\mathrm{xx},j}-P_j)\hfill \end{array}$$

(56)

If $M_Q$ is large enough, the difference between $M_Q/{(M_Q-1)}^2$ and $1/(M_Q-2)$ is negligible. We set M_Q = 10 according to the Remark 4 in the reference Jiang and Cai. ³⁹ So, $W_k$ can be approximated as

$$\begin{array}{c}\multicolumn{1}{c}{W_k\approx \frac{1}{M_Q-2}\sum _{i=k-M_Q+1}^{k-1}({\zeta }_j-{\overline{\zeta }}_k){({\zeta }_j-{\overline{\zeta }}_k)}^T}\\ -\frac{1}{M_Q-1}\sum _{j=k-M_Q+1}^k(P_{\mathrm{xx},j}-P_j)\hfill \end{array}$$

(57)

and

$$\mathrm{\Delta }{Q}_k=\frac{1}{M_Q-1}({\zeta }_k-{\overline{\zeta }}_k){({\zeta }_k-{\overline{\zeta }}_k)}^T-\frac{1}{M_Q}(P_{\mathrm{xx},k}-P_k)$$

(58)

Obviously, the above equation has the equivalent form with $Q_{k-1}$ ; therefore, equation (55) can be duduced. Similarly, the recursive estimation equation of ${\zeta }_k$ is given by:

$${\overline{\zeta }}_k=\frac{M_Q-1}{M_Q}{\overline{\zeta }}_{k-1}+\frac{1}{M_Q}{\zeta }_k$$

(59)

So, we can use equations (49), (55), (56), and (59) to update the covariance of process noise. Provided that the covariance of process noise stays constant, the recursive relation between $Q_{k-1}$ and $Q_k$ can then be described as:

$${\hat{Q}}_k=\frac{M_{Q_k}-1}{M_{Q_k}}{\hat{Q}}_{k-1}+\mathrm{\Delta }{Q}_k$$

(60)

AFGBICKF algorithm is formulated in the Algorithm 2.

Algorithm 2.

AFGBICKF algorithm.

Initialize ${\hat{\mu }}_0$ , $P_0$ , ${\hat{Q}}_0$ , $M_Q$ . For $k=1,2,\cdots$ Step 1. Time Update Calculate the predicted state ${\overline{\mu }}_k$ and covariance ${\overline{P}}_k$ , using equations (24) and (25). Step 2. Measurement update. Step 2.1 Evaluate the predicted measurement estimate, square root of the cross-covariance, and innovation covariance, using equations (28)–(30). Step 2.2 Evaluate the modified gain $G_{\mathrm{new}}$ using equation (32), state estimate ${\hat{\mu }}_k$ and covariance $P_k$ using equations (33) and (34). Step 3. Estimate the process noise covariance recursively. ${\hat{Q}}_k=(M_Q-1){\hat{Q}}_{k-1}/M_Q+\mathrm{\Delta }{Q}_k$ where $\mathrm{\Delta }{Q}_k$ is defined as equation (58). End.

Case study: Re-entry ballistic target tracking

In the following, we apply the proposed FGBICKF and AFGBICKF algorithm to the RBT tracking. ⁴⁰ Firstly, the effect of fractional order on FGBICKF is analyzed. Then we compare FGBICKF with CKF and ICKF. Furthermore, we compare the standard deviation (STD) of FGBICKF with PCRLB and analyze errors from FGBICKF. Afterwards, we compare AFGBICKF with FGBICKF, while analyzing the influence of various initial process noise covariance on the performance of AFGBICKF algorithm.

RBT dynamics model

Figure 1 shows the geometry of the RBT. P is re-entry ballistic target and the radar is situated at the surface of Earth. The two coordinate systems: the Earth-centered inertial coordinate system (ECI-CS, $O{x}_I{y}_I{z}_I$ ) and East-North-Up coordinates system (ENU-CS, $O_s\mathrm{xyz}$ ) are presented in the Figure 1. The ECI-CS is a right-handed system with the origin O at Earth’s center, $O{x}_I$ pointing to the vernal equinox direction, and the axis $O{z}_I$ pointing to the direction of the North Pole N. Its fundamental plane $O{x}_I{y}_I$ coincides with Earth’s equatorial plane. The ENU-CS has its origin at the location of the radar, and $z$ is directed along the local vertical. $x$ and $y$ lie in a local horizontal plane, $x$ pointing to the east, and $y$ pointing to the north, respectively.

Figure 1.

Geometry of radar and RBT (the ECI-CS (OxIyIzI) and ENU-CS (Osxyz).

In order to obtain the state equation, we make two hypotheses: one is that Earth is spherical and non-rotating, and another is that only the gravity and drag are the forces acting on the RBT. ⁴⁰ According to the transformation relationship of ECI-CS and ENU-CS, we derive the kinematics of the RBT with unknown ballistic coefficient in the ENU-CS. And we model the discrete-time stochastic nonlinear state equation of RBT:

$${\eta }_k=\mathrm{\Phi }{\eta }_{k-1}+\mathrm{\Gamma }\mathrm{\Psi }({\eta }_{k-1})+w_{k-1}$$

(61)

here, ${\eta }_k=[x_{k}v{x}_k{y}_{k}v{y}_k{z}_{k}v{z}_k{\beta }_k{]}^T$ is the RBT’s state, ${\beta }_k$ (kg/m²) is ballistic coefficient, $\mathrm{\Phi }$ and $\mathrm{\Gamma }$ are described in the following:

$$\begin{array}{c}\mathrm{\Phi }=\left[\begin{array}{cccc}\delta \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \delta \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \delta \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],\delta =\left[\begin{array}{cc}1\hfill & t\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right],\hfill \\ \mathrm{\Gamma }=\left[\begin{array}{ccc}\kappa \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \kappa \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \kappa \hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right],\kappa =\left[\begin{array}{c}{t}^2/2\hfill \\ \multicolumn{1}{c}{t}\end{array}\right],\hfill \end{array}$$

(62)

$$\mathrm{\Psi }({\eta }_{k-1})=\left[\begin{array}{c}-\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}v{x}_{k-1}-\frac{\mu x_{k-1}}{r_{k-1}^3}\hfill \\ \multicolumn{1}{c}{-\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}v{y}_{k-1}-\frac{\mu y_{k-1}}{r_{k-1}^3}}\\ \multicolumn{1}{c}{-\frac{\rho (h_{k-1})}{2{\beta }_{k-1}}{V}_{k-1}v{z}_{k-1}-\frac{\mu (z_{k-1}+R_e)}{r_{k-1}^3}}\end{array}\right]$$

(63)

Here t (in second, i.e. s) is the time interval between radar measurements, $\mu =3.986005\times {10}^{14}{m}^3/s^2$ and $R_e=6371004m$ are Earth’s gravitational constant and Earth radius, respectively. And $r_{k-1}=\sqrt{x_{k-1}^2+y_{k-1}^2+{(z_{k-1}+R_e)}^2}$ , $V_{k-1}=\sqrt{{\mathrm{vx}}_{k-1}^2+{\mathrm{vy}}_{k-1}^2+{\mathrm{vz}}_{k-1}^2}$ , and $h_{k-1}=r_{k-1}-R_e$ · $\rho (h)$ (kg/m³) is the air density, which can be approximately modeled as $\rho =c_1{e}^{-c_{2}h}$ at height below 90 km (c₁, c₂ are dimensionless and constant). ⁴¹ $w_k$ is Gaussian noise with zero mean and covariance matrix ( $Q_k$ ) ⁴² :

$$Q_k=\left[\begin{array}{cccc}{q}_1\vartheta \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& q_1\vartheta \hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & q_1\vartheta \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & q_{2}t\hfill \end{array}\right],\vartheta =\left[\begin{array}{cc}{t}^3/3\hfill & t^2/2\hfill \\ \multicolumn{1}{c}{t^2/2}& t\hfill \end{array}\right]$$

(64)

Here $q_1$ (in m²/ s²) and $q_2$ (in kg²/(m⁴ s)) control the amount of process noise and RBT’s ballistic coefficient, respectively.

The radar collects the measurements: the range (R), elevation(E), and azimuth (A), which are presented in Figure 2.

Figure 2.

Range (R), elevation(E), and azimuth(A) collected by radar.

The measurement equation in ENU-CS is depicted as follows:

$$l_k=h(x_k)+v_k$$

(65)

where $l_k=[R_k{E}_k{A}_k{]}^T$ , $h(x_k)=[\begin{array}{cc}\sqrt{x_k^2+y_k^2+z_k^2}\arctan z_k/\sqrt{x_k^2+y_k^2}\hfill & \arctan y_k/x_k\hfill \end{array}{]}^T$ . The range $R_k$ , elevation $E_k$ , and azimuth $A_k$ at time instant k are defined: $R_k=\sqrt{x_k^2+y_k^2+z_k^2}+v_R$ , $E_k=\arctan z_k/\sqrt{x_k^2+y_k^2}+v_E$ , $A_k=\arctan y_k/x_k+v_A$ . The measurement noise $v_k=[\begin{array}{ccc}{v}_R\hfill & v_E\hfill & v_A\hfill \end{array}{]}^T$ is the white Gaussian noise, which is with the zero-mean and the covariance matrix $N_k=\mathrm{diag}([\begin{array}{ccc}{\sigma }_R^2\hfill & {\sigma }_E^2\hfill & {\sigma }_A^2\hfill \end{array}])$ , ${\sigma }_R$ , ${\sigma }_E$ , and ${\sigma }_A$ are the standard deviations of range, elevation, and azimuth, respectively.

To evaluate the performance of the proposed filters, we use the performance metrics: root mean-square error (RMSE) and average accumulated mean-square root error (AMSRE). ⁴³ The RMSE and AMSRE in position at $k$ time instant are defined as follows:

$$\begin{array}{c}\mathrm{RMS}{E}_p(k)=\frac{1}{M_t}\sum _{i=1}^{M_t}\hfill \\ \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}\hfill \end{array}$$

(66)

$$\begin{array}{c}\mathrm{AMSR}{E}_p=\frac{1}{N_t}\sum _{i=1}^{N_t}\hfill \\ \sqrt{\frac{1}{M_t}\sum _{k=1}^{M_t}{(x_k^{(i)}-{\hat{x}}_k^{(i)})}^2+{(y_k^{(i)}-{\hat{y}}_k^{(i)})}^2+{(z_k^{(i)}-{\hat{z}}_k^{(i)})}^2}\hfill \end{array}$$

(67)

Here $(x^{(i)},y^{(i)},z^{(i)})$ and $({\hat{x}}^{(i)},{\hat{y}}^{(i)},{\hat{z}}^{(i)})$ are the true and estimated position at the $i$ th Monte Carlo run, M_t is the total Monte Carlo runs, and N_t is the total number of measurement data points. Similarly, we can obtain RMSE and AMSRE in velocity and ballistic coefficient. The obtained RMSEs and AMSREs in the position, velocity and ballistic coefficient in this paper were averaged over 100 independent Monte Carlo runs.

Simulations and analysis

Influence of fractional difference order (α) on FGBICKF

In the simulation, we set the parameters λ₁, T, q₁, and q₂ as 2, 0.1 s, 5 m²/s³, and 5 kg²/(m⁴s). We initialize the position x₀, y₀, and z₀ as 232, 232, and 90 km and module of the velocity v₀ as 3000 m/s. We select the initial elevation E₀, azimuth angle A₀ and initial ballistic coefficient β₀ as 210°, 45°, and 4000 kg/m². Then the true initial state x₀ = [232 km–1837 m/s, 232 km–1837 m/s, 90 km–1500 m/s, 4000 kg/m²]^T is obtained, and the covariance P₀ = diag([100²  50²  100²  50² 100²  50²  200²]) is set. The initial state estimate is randomly chosen from ${\hat{x}}_0~N(x_0,P_0)$ . we set the standard deviations of the range errors (σ_R) as 100 m, the elevation angle (σ_E), azimuth angle (σ_A) as 0.017 rad.

In the FGBICKF, the Kalman gain is modified with fractional gain. We have done experiments and found that FGBICKF converged when α was less than 1(α < 1). The simulation is carried out by selecting the fractional order α as 0.05, 0.1, 0.2, 0.5, and 0.8, respectively. Figures 3 to 5 presents the RMSEs of FGBICKF with various fractional orders for the position, the velocity, and the ballistic coefficient.

Figure 3.

RMSE of FGBICKF for various values of α in position.

Figure 4.

RMSE of FGBICKF for various values of α in velocity.

Figure 5.

RMSE of FGBICKF for various values of α in ballistic coefficient.

As can been seen from Figures 3 to 5, FGBICKF clearly converges. It is observed that there is a tradeoff between convergence rate and estimation accuracy from RMSEs in position, velocity, and ballistic coefficient. For α = 0.1, the RMSE of FGBICKF in the position is lower than that of other values of α except α = 0.05, and the RMSE of FGBICKF in the velocity is lower than that of other values of α except α = 0.05 before the estimation time 50 s. The RMSE of FGBICKF in the ballistic coefficient is higher than that of other values of α except α = 0.05 after the estimation time 50 s. It has been shown that FGBICKF has better performance when the fractional order is set as α = 0.1.

Comparison of FGBICKF and CKF, ICKF

In the subsection, we compare the performance of FGBICKF with fractional order $\alpha =0.1$ with that of CKF and ICKF when they are applied to state estimation in the RBT tracking problem. Here, the values of the parameters ( ${\lambda }_1$ , ${\hat{x}}_0$ , $P_0$ , $Q_k$ , and $R_k$ ) are set the same as those in the first subsection. Figures 6 to 8 present RMSEs in position, velocity, and the ballistic coefficient for FGBICKF, CKF, and ICKF.

Figure 6.

RMSEs in position for various filters.

Figure 7.

RMSEs in velocity for various filters.

Figure 8.

RMSEs in ballistic coefficient for various filters.

Obviously, in terms of effectiveness, Figure 6 to 8 show the high accuracy of FGBICKF, compared with CKF and ICKF. Particularly, Figure 8 illustrates that the estimates of FGBICKF in ballistic coefficient are more accurate than those of CKF and ICKF after 40 s. In FGBICKF, the Kalman gain is modified with fractional characteristics. The information of variations in the previous gain is used to evaluate the next state with the help of fractional derivative. Modified gain gives the better performance for state estimation. From Figures 6 to 8, we can also see that CKF and ICKF achieves mostly the same level of estimation performance in position, velocity, and ballistic coefficient.

Moreover, the runtimes of UKF, CKF, ICKF, and FGBICKF are about 0.76066, 0.24731, 0.26729, and 0.36222 s, respectively. We see the runtime of FGBICKF is less than that of UKF due to its requirement for computing sigma points, and FGBICKF has a slightly higher computational complexity, compared to CKF and ICKF. FGBICKF is slightly slower due to its requirements on evaluating more past gain.

Error analysis on FGBICKF

Cramér-Rao lower bound (CRLB) provides the best error analysis for filter and is used to evaluate the performance of filters. ⁴⁴ $\mathrm{CRL}{B}_k$ (i.e. CRLB at $k$ time instant) is defined as:

$$\mathrm{CRL}{B}_k=J_k^{-1}$$

(68)

where $J_k$ is the fisher information matrix and is defined as

$$J_k=E[-\frac{{\partial }^2\ln p(x_{1:k},l_{1:k})}{\partial x_k^2}]$$

(69)

Here, the posterior distribution function of $x_{1:k},l_{1:k}$ is $p(x_{1:k},l_{1:k})$ . In practice, $J_k$ in equation (69) is hard to be calculated, so we use posterior Fisher information matrix $J_k$ , which is calculated using the recursive form in equation (70) to obtain the posterior CRLB (PCRLB). ³⁶

$$J_{k+1}=D_k^{22}-D_k^{21}{[D_k^{11}+J_k]}^{-1}{D}_k^{12}$$

(70)

Here, $D_k^{11}$ , $D_k^{12}$ , $D_k^{21}$ , and $D_k^{22}$ are defined as:

$$D_k^{11}=E\{-\frac{{\partial }^2\ln p(x_k|x_{k-1})}{{\partial }^2{x}_{k-1}^2}\}$$

(71)

$$D_k^{12}=(D_k^{21}{)}^T=E\{-\frac{{\partial }^2\ln p(x_k|x_{k-1})}{\partial x_{k-1}\partial x_k}\}$$

(72)

$$D_k^{22}=E\{-\frac{{\partial }^2\ln p(x_k|x_{k-1})}{{\partial }^2{x}_k^2}\}$$

(73)

For nonlinear system with addictive Gaussian noise in equations (3) and (4), $D_k^{11}$ , $D_k^{12}$ , $D_k^{21}$ , and $D_k^{22}$ are calculated as

$$D_k^{11}=E\{[{\nabla }_{x_k}{\mathrm{f}}^T(x_k)]Q_k^{-1}[{\nabla }_{x_k}{\mathrm{f}}^T(x_k){]}^T\}$$

(74)

$$D_k^{12}=-E\{{\nabla }_{x_k}{\mathrm{f}}^T(x_k)\}{Q}_k^{-1}$$

(75)

$$D_k^{21}=(D_k^{12}{)}^T$$

(76)

$$D_k^{22}=Q_k^{-1}+E\{[{\nabla }_{x_{k+1}}{\mathrm{h}}^T(x_{k+1})]R_{k+1}^{-1}[{\nabla }_{x_{k+1}}{\mathrm{h}}^T(x_{k+1}){]}^T\}$$

(77)

The initial information matrix ( $J_0$ ) is calculated as:

$$J_0=E\{-{\mathrm{\Delta }}_{x_0}^{x_0}\log p(x_0)\}$$

(78)

Here, ${\nabla }_x=[\partial /\partial x_1,\partial /\partial x_2,\cdots ,\partial /\partial x_n]$ and ${\mathrm{\Delta }}_x^x={\nabla }_x{\nabla }_x^T$ .

Next, we compare the STD of FGBICKF with PCRLB, as shown in Figure 9. Here, the values of the parameters ( ${\lambda }_1$ , $\alpha$ , ${\hat{x}}_0$ , $P_0$ , $Q_k$ , and $R_k$ ) are selected the same as those in the first subsection.

Figure 9.

Comparison of STD of various filters and PCRLB: (a) x-axis position, (b) x-axis velocity, (c) y-axis position, (d) y-axis velocity, (e) z-axis position, (f) z-axis velocity, and (g) ballistic coefficient.

From the Figure 9(a) to (g), we see that the standard deviations of the position, velocity, and ballistic coefficient of FGBICKF are smaller than that of CKF and ICKF, and they are very close to PCRLB. These results further prove that FGBICKF is a state estimation algorithm with better performance.

Comparison of FGBICKF and AFGBICKF

In the simulation, the values of parameters ( ${\lambda }_1$ , $\alpha$ , ${\hat{x}}_0$ , $P_0$ , $Q_k$ , and $R_k$ ) are set the same as those in the first subsection. Figures 10 to 12 present RMSEs in position, velocity and ballistic coefficient for FGBICKF and AFGBICKF when the initial process noise covariance is $Q_0=1000*Q_k$ .

Figure 10.

RMSEs in position for FGBICKF and AFGBICKF.

Figure 11.

RMSEs in velocity for FGBICKF and AFGBICKF.

Figure 12.

RMSEs in ballistic coefficient for FGBICKF and AFGBICKF.

We see the obvious difference on the performance between the FGBICKF and AFGBICKF. Figures 10 to 12 shows that AFGBICKF keeps converging and its performance is better than that of FGBICKF. The simulation results prove that AFGBICKF can suppress the influence of the initial process noise covariance on state estimation.

Moreover, we compute the AMSREp, AMSREv (AMSRE in velocity), and AMSREb (AMSRE in ballistic coefficient) for FGBICKF and AFGBICKF, respectively, as shown in Figure 13.

Figure 13.

AMSREs for FBCICKF and AFGBICKF.

From Figure 13, we can see that AFGBICKF has achieved comparable estimation accuracy to FGBICKF. The simulation results have demonstrated the prominent improvement over FGBICKF because the AFGBICKF incorporates the recursive procedures of estimating the process noise covariance.

AFGBICKF’s adaptiveness to initial process noise covariance

In the simulation, the values of the parameters ( ${\lambda }_1$ , $\alpha$ , ${\hat{x}}_0$ , $P_0$ , and $R_k$ ) are set the same as those in the first subsection and initial process covariance $Q_0=Q_k$ . Figures 14 to 16 present RMSEs in position, velocity and the ballistic coefficient using the AFGBICKF algorithm when the initial process noise covariances are set to $Q_0/1000$ , $Q_0$ , and $1000*Q_0$ , respectively.

Figure 14.

Position RMSE for AFGBICKF with various initial process noise covariances.

Figure 15.

Velocity RMSE for AFGBICKF with various initial process noise covariances.

Figure 16.

Ballistic coefficient RMSE for AFGBICKF with various initial process noise covariances.

From Figures 14 to 16, we can see that when the initial process noise covariance changes in a big range, from being small (1/1000 of $Q_0$ ) to big (1000 times of $Q_0$ ), AFGBICKF always converges and achieves almost same estimation performance. This validates the superiority of AFGBICKF and the correctness of our theoretical analysis in section 3.2. Moreover, a bigger initial estimate of process noise covariance gives better estimation performance. Thus, when accurate statistic of the process noise is unknown, a larger initial value is recommended to ensure the stability and convergence of filtering algorithms.

Conclusion

This paper proposes a novel filter FGBICKF to estimate the states for SNDS with Gaussian noise. Utilizing the fractional order gain, we further develop the FGBICKF with recursive process noise covariance estimation. FGBICKF uses fractional derivative of previous filter gains as feedback to current filter gain for estimating next states. The simulations on RBT tracking have proven FGBICKF’s superior estimation performance. Meanwhile, we have computed the PCRLB and compared it with relevant standard deviations of different filters, proving the effectiveness of FGBICKF. Moreover, the application of AFGBICKF to RBT tracking problem with unknown process noise covariance has shown AFGBICKF’s adaptiveness even when the initial process noise covariance changes dramatically. However, FGBICKF has some limitations. Firstly, the accuracy of the proposed filters depends on the accuracy of the system dynamics and measurement process. If the model does not accurately represent the real system, the filter’s estimates may deviate significantly from true states. Secondly, as the dimensionality of state increases, the matrix operations involved in the calculations become computationally intensive and require significantly additional computational resources. In the future, we will continue our research on addressing these limitations to estimate the states for SNDS more accurately and robustly.