Maximum Correntropy Unscented Kalman Filter for Spacecraft Relative State Estimation

Abstract

A new algorithm called maximum correntropy unscented Kalman filter (MCUKF) is proposed and applied to relative state estimation in space communication networks. As is well known, the unscented Kalman filter (UKF) provides an efficient tool to solve the non-linear state estimate problem. However, the UKF usually plays well in Gaussian noises. Its performance may deteriorate substantially in the presence of non-Gaussian noises, especially when the measurements are disturbed by some heavy-tailed impulsive noises. By making use of the maximum correntropy criterion (MCC), the proposed algorithm can enhance the robustness of UKF against impulsive noises. In the MCUKF, the unscented transformation (UT) is applied to obtain a predicted state estimation and covariance matrix, and a nonlinear regression method with the MCC cost is then used to reformulate the measurement information. Finally, the UT is adopted to the measurement equation to obtain the filter state and covariance matrix. Illustrative examples demonstrate the superior performance of the new algorithm.

1. Introduction

The relative state estimation problem is rather important in space communication networks, including inter-satellite link (ISL) establishment, communication support, multiple spacecraft formation flying, and spacecraft rendezvous and docking [1,2,3,4]. Many relative motion scenes make use of the radar or/and ladar to obtain the relative measurement information among spacecrafts, which is necessary to determine the relative motions.

As for the state estimation problem, the Kalman filter (KF) is often used for a linear dynamic system, which is, in essence, a recursive minimum $l_2$-norm linear filter [5,6,7]. However, from [8], we know that the relative state estimation in a circular reference orbit involves a linear state Clohessy-Wiltshire (CW) equation and a nonlinear measurement equation. Especially, when the reference orbit is elliptic, the dynamic state equation is also nonlinear [9]. To handle the nonlinear problems, the widely used methods are the extended Kalman filter (EKF) [10,11,12] and unscented Kalman filter (UKF) [11,12,13]. The EKF is a popular nonlinear extension of KF, which makes use of the first order Taylor series expansions to approximate the nonlinear system and applies the KF to this approximation. Nevertheless, the approximation is crude and it may lead to filter divergence if the function exhibits highly non-linear characteristics. Furthermore, the derivation of the Jacobian matrices is cumbersome and often results in the implementation difficulties. In UKF, the probability distribution of the state is approximated by a set of deterministically selected sigma points and propagated though the non-linear process and measurement equations. UKF does not need to calculate the cumbersome Jacobian matrices and provides derivative-free higher-order approximations to achieve more accurate performances than the EKF. Unfortunately, the EKF and UKF may fail to perform well when the measurements are disturbed by some heavy-tailed non-Gaussian noises, which arises naturally in many real applications of engineering. The main reason for this is that both methods are based on the minimum $l_2$-norm technique and thus exhibit sensitivity to heavy-tailed noises [14]. To address this problem, some robust methods have been proposed. The Huber’s generalized maximum likelihood methodology is very common [15,16], which is a combined minimum $l_1$ and $l_2$-norm technique.

Besides Huber’s robust statistics, the optimization criteria in information theoretic learning (ITL) [17,18] provide an alternative effective way, using information theoretic quantities estimated directly from the data as the optimization costs instead of using the usual second-order statistical measures (such as variance and covariance). The ITL costs can capture higher-order statistics of data, and can be used as a robust adaptation cost to achieve excellent performance in a number of applications [19,20,21,22,23,24,25,26,27,28,29,30,31]. Particularly, an optimization criterion based on correntropy, called the maximum correntropy criterion (MCC) [17,18], has recently been successfully applied in robust adaptive filtering in the presence of heavy-tailed non-Gaussian noises [17,28,29,30,31,32].

In this paper, a non-linear Kalman type filter based on MCC has been derived, namely the maximum correntropy unscented Kalman filter (MCUKF). In MCUKF, the UT is applied to obtain a predicted state estimation and covariance matrix, and a non-linear regression model under MCC is used to reformulate the measurement information. Then, the UT is adopted to the measurement equation to obtain the filter state and covariance matrix. Similar to the original UKF, the new filter also has a recursive structure and is suitable for online implementation. It is worth noting that the proposed MCUKF is different from the algorithm in [33].

The organization of the paper is as follows. In Section 2, we provide a short review of the MCC and the traditional UKF. In Section 3, the MCUKF algorithm is derived. Two illustrative examples are then provided in Section 4 to show the excellent performance of the MCUKF. Finally, Section 5 gives the conclusion.

2. Preliminaries

2.1. Maximum Correntropy Criterion

Given two random variables $X,Y\in \mathbb{R}$ with joint distribution function ${\mathrm{F}}_{XY}(x,y)$, the correntropy is a generalized similarity measure between X and Y defined by

$$V(X,Y)=\mathrm{E}\left[\kappa (X,Y)\right]=\int \kappa (x,y)d{\mathrm{F}}_{XY}(x,y),$$

(1)

where E is the expectation operator, and $\kappa (·,·)$ denotes a shift-invariant kernel following the Mercer condition. In this paper, the kernel function is chosen as the Gaussian kernel given by

$$\kappa (x,y)={\mathrm{G}}_{\sigma }\left(e\right)=\exp \left(-\frac{e^2}{2{\sigma }^2}\right),$$

(2)

where $e=x-y$, and $\sigma >0$ is the kernel bandwidth.

In many practical applications, we have only a limited amount of data and the joint distribution ${\mathrm{F}}_{XY}$ is unknown. In these cases, the correntropy can be estimated by the sample mean estimator:

$$\widehat{V}(X,Y)=\frac{1}{N}\sum _{i=1}^N{\mathrm{G}}_{\sigma }\left(e\left(i\right)\right),$$

(3)

where $e\left(i\right)=x\left(i\right)-y\left(i\right)$, with ${\left\{x\left(i\right),y\left(i\right)\right\}}_{i=1}^N$ being N samples drawn from ${\mathrm{F}}_{XY}$.

Taking the Taylor series expansion of the Gaussian kernel yields

$$V(X,Y)=\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{2^n{\sigma }^{2n}n!}\mathrm{E}\left[{(X-Y)}^{2n}\right].$$

(4)

It can be seen that the correntropy is a weighted sum of all even order moments of the error variable $X-Y$. In addition, the kernel bandwidth appears as a parameter weighting the second order and higher order moments. It is noted that the correntropy will be dominated by the second order moment when the kernel bandwidth is very large compared to the dynamic range of the data.

Given a sequence of error data ${\left\{e\left(i\right)\right\}}_{i=1}^N$, the cost function of MCC is given by

$$J_{MCC}=\frac{1}{N}\sum _{i=1}^N{\mathrm{G}}_{\sigma }\left(e\left(i\right)\right).$$

(5)

Assume $\mathbf{W}$ is a parameter vector of an adaptive system to learn, and let $x\left(i\right)$ and $y\left(i\right)$ denote the model output and the desired response, respectively. The MCC based learning can be formulated as the following optimization problem:

$$\widehat{\mathbf{W}}=\underset{\mathbf{W}\in \mathbf{\Omega }}{\arg max}\frac{1}{N}\sum _{i=1}^N{\mathrm{G}}_{\sigma }\left(e\left(i\right)\right),$$

(6)

where $\widehat{\mathbf{W}}$ denotes the optimal solution, and Ω denotes a feasible set of the parameter.

2.2. Unscented Kalman Filter

The UKF provides a suitable alternative for EKF to deal with nonlinear systems. In this paper, we discuss a discrete-time nonlinear system described as

$$\mathbf{x}\left(k\right)=\mathrm{f}\left(k-1,\mathbf{x}(k-1)\right)+\mathbf{q}(k-1),$$

(7)

$$\mathbf{y}\left(k\right)=\mathrm{h}\left(k,\mathbf{x}\left(k\right)\right)+\mathbf{r}\left(k\right),$$

(8)

where $\mathbf{x}\left(k\right)\in {\mathbb{R}}^n$ denotes the n-dimensional state vector of the system at time step k, $\mathbf{y}\left(k\right)\in {\mathbb{R}}^m$ denotes an m-dimensional measurement vector. f represents a nonlinear system function, and h is a nonlinear measurement function and they are assumed to be continuously differentiable. The process noise $\mathbf{q}(k-1)$ and measurement noise $\mathbf{r}\left(k\right)$ are mutually uncorrelated, with zero mean values and covariance matrices

$$\mathrm{E}\left[\mathbf{q}(k-1){\mathbf{q}}^T(k-1)\right]=\mathbf{Q}(k-1),\mathrm{E}\left[\mathbf{r}\left(k\right){\mathbf{r}}^T\left(k\right)\right]=\mathbf{R}\left(k\right).$$

(9)

In general, similar to KF, the UKF includes two steps, namely predict and update:

2.2.1. Predict

First, a set of $2n+1$ sigma points are computed by the estimated state and covariance matrix of last time step:

$$\begin{array}{cc}\hfill {\chi }^0(k-1|k-1)=& \widehat{\mathbf{x}}(k-1|k-1),\hfill \\ \hfill {\chi }^i(k-1|k-1)=& \widehat{\mathbf{x}}(k-1|k-1)\hfill \\ \hfill & +{\left(\sqrt{(n+\lambda )\mathbf{P}(k-1|k-1)}\right)}_i,i=1,\dots ,n,\hfill \\ \hfill {\chi }^i(k-1|k-1)=& \widehat{\mathbf{x}}(k-1|k-1)\hfill \\ \hfill & -{\left(\sqrt{(n+\lambda )\mathbf{P}(k-1|k-1)}\right)}_{i-n},i=n+1,\dots ,2n,\hfill \end{array}$$

(10)

where ${\left(\sqrt{(n+\lambda )\mathbf{P}(k-1|k-1)}\right)}_i$ is the i-th column of the matrix square root of $(n+\lambda )\mathbf{P}(k-1|k-1)$, n denotes the dimension of the state vector, and λ represents a compound scaling factor defined by

$$\lambda ={\alpha }^2(n+\epsilon )-n,$$

(11)

where α determines the extent of the spread of the sigma points, which is usually chosen as $0<\alpha \le 1$, ε is a scaling factor and is usually set to $3-n$.

The transformed set is given through the process model

$${{\chi }^i}^{*}\left(k\right|k-1)=\mathrm{f}\left(k-1,{\chi }^i(k-1|k-1)\right),i=0,\dots ,2n.$$

(12)

Then, the predicted state mean and associated covariance matrix are calculated by

$$\widehat{\mathbf{x}}\left(k\right|k-1)=\sum _{i=0}^{2n}{w}_m^i{{\chi }^i}^{*}\left(k\right|k-1),$$

(13)

$$\begin{array}{cc}\hfill \mathbf{P}\left(k\right|k-1)=& \sum _{i=0}^{2n}{w}_c^i\left[{{\chi }^i}^{*}\left(k\right|k-1)-\widehat{\mathbf{x}}\left(k\right|k-1)\right]\hfill \\ & \times {\left[{{\chi }^i}^{*}\left(k\right|k-1)-\widehat{\mathbf{x}}\left(k\right|k-1)\right]}^T+\mathbf{Q}(k-1),\hfill \end{array}$$

(14)

where the corresponding weights of the state and covariance matrix are defined as

$$\begin{array}{c}{w}_m^0={\displaystyle \frac{\lambda }{(n+\lambda )}},\hfill \\ w_c^0=\lambda /(n+\lambda )+(1-{\alpha }^2+\beta ),\hfill \\ w_m^i=w_c^i={\displaystyle \frac{1}{2(n+\lambda )}},i=1,\dots ,2n.\hfill \end{array}$$

(15)

with β being related to prior knowledge of the distribution of $\mathbf{x}\left(k\right)$ (for Gaussian noises, $\beta =2$ is optimal).

2.2.2. Update

Similarly, a set of $2n+1$ sigma points are computed by the predicted state mean and covariance matrix

$$\begin{array}{cc}\hfill {\chi }^0\left(k\right|k-1)=& \widehat{\mathbf{x}}\left(k\right|k-1),\hfill \\ \hfill {\chi }^i\left(k\right|k-1)=& \widehat{\mathbf{x}}\left(k\right|k-1)\hfill \\ \hfill +& {\left(\sqrt{(n+\lambda )\mathbf{P}\left(k\right|k-1)}\right)}_i,i=1,\dots ,n,\hfill \\ \hfill {\chi }^i\left(k\right|k-1)=& \widehat{\mathbf{x}}\left(k\right|k-1)\hfill \\ \hfill -& {\left(\sqrt{(n+\lambda )\mathbf{P}\left(k\right|k-1)}\right)}_{i-n},i=n+1,\dots ,2n.\hfill \end{array}$$

(16)

The transformed set is given through the measurement model

$${\gamma }^i\left(k\right)=\mathrm{h}(k,{\chi }^i\left(k\right|k-1)),i=0,\dots ,2n.$$

(17)

Then, the predicted measurement mean and covariance matrix are calculated by

$$\widehat{\mathbf{y}}\left(k\right)=\sum _{i=0}^{2n}{w}_m^i{\gamma }^i\left(k\right),$$

(18)

$${\mathbf{P}}_{\mathbf{yy}}=\sum _{i=0}^{2n}{w}_c^i\left[{\gamma }^i\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right]{\left[{\gamma }^i\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right]}^T+\mathbf{R}\left(k\right).$$

(19)

In addition, the state-measurement cross-covariance matrix is computed by

$${\mathbf{P}}_{\mathbf{xy}}\left(k\right)=\sum _{i=0}^{2n}{w}_c^i\left[{\chi }^i\left(k\right|k-1)-\widehat{\mathbf{x}}\left(k\right|k-1)\right]{\left[{\gamma }^i\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right]}^T.$$

(20)

Next, the Kalman gain is given as

$$\mathbf{K}\left(k\right)={\mathbf{P}}_{\mathbf{xy}}{\mathbf{P}}_{\mathbf{yy}}^{-1}.$$

(21)

Finally, the filter state in this case is obtained by

$$\widehat{\mathbf{x}}\left(k\right|k)=\widehat{\mathbf{x}}\left(k\right|k-1)+\mathbf{K}\left(k\right)\left(\mathbf{y}\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right).$$

(22)

Additionally, the filter covariance is recursively updated as

$$\mathbf{P}\left(k\right|k)=\mathbf{P}\left(k\right|k-1)-\mathbf{K}\left(k\right){\mathbf{P}}_{\mathbf{yy}}{\mathbf{K}}^T\left(k\right).$$

(23)

3. Unscented Kalman Filter under MCC

To improve the robustness of the UKF, we present the idea to combine the MCC cost function with the UKF framework to derive a novel UKF, which may perform much better in non-Gaussian noise environments, since correntropy contains second and higher order moments of the error.

First, we consider a nonlinear model given by Equation (7) and Equation (8), combine expression Equations (13) and (14) and recast the nonlinear regression model as

$$\left[\begin{array}{c}\widehat{\mathbf{x}}\left(k\right|k-1)\\ \mathbf{y}\left(k\right)\end{array}\right]=\left[\begin{array}{c}\mathbf{x}\left(k\right)\\ \mathrm{h}(k,\mathbf{x}(k\left)\right)\end{array}\right]+\psi \left(k\right),$$

(24)

where $\psi \left(k\right)$ denotes by

$$\psi \left(k\right)=\left[\begin{array}{c}\eta \left(\mathbf{x}\right(k\left)\right)\\ \mathbf{r}\left(k\right)\end{array}\right],$$

(25)

with

$$\begin{array}{cc}& \mathrm{\Psi }\left(k\right)\hfill \\ \hfill =& \mathrm{E}\left[\psi \left(k\right){\psi }^T\left(k\right)\right]\hfill \\ \hfill =& \left[\begin{array}{cc}\mathbf{P}\left(k\right|k-1)& 0\\ 0& \mathbf{R}\left(k\right)\end{array}\right]\hfill \\ \hfill =& \left[\begin{array}{cc}{\mathbf{T}}_p\left(k\right|k-1){\mathbf{T}}_p^T\left(k\right|k-1)& 0\\ 0& {\mathbf{T}}_r\left(k\right){\mathbf{T}}_r^T\left(k\right)\end{array}\right]\hfill \\ \hfill =& \mathbf{T}\left(k\right){\mathbf{T}}^T\left(k\right),\hfill \end{array}$$

(26)

where $\mathbf{T}\left(k\right)$ can be obtained by Cholesky decomposition of $\mathrm{\Psi }\left(k\right)$. Left multiplying both sides of Equation (24) by ${\mathbf{T}}^{-1}\left(k\right)$, we have the following equation

$$\mathbf{D}\left(k\right)=\mathrm{g}(k,\mathbf{x}(k\left)\right)+\mathbf{e}\left(k\right),$$

(27)

where $\mathbf{D}\left(k\right)={\mathbf{T}}^{-1}\left(k\right)\left[\begin{array}{c}\widehat{\mathbf{x}}\left(k\right|k-1)\\ \mathbf{y}\left(k\right)\end{array}\right]$, $\mathrm{g}(k,\mathbf{x}\left(k\right))={\mathbf{T}}^{-1}\left(k\right)\left[\begin{array}{c}\mathbf{x}\left(k\right)\\ \mathrm{h}(k,\mathbf{x}(k\left)\right)\end{array}\right]$, $\mathbf{e}\left(k\right)={\mathbf{T}}^{-1}\left(k\right)\psi \left(k\right)$ and the covariance of $\mathbf{e}\left(k\right)$ is the identity matrix.

Then, we define the following cost function $J_L\left(\mathbf{x}\left(k\right)\right)$ based MCC:

$$J_L\left(\mathbf{x}\left(k\right)\right)=\sum _{i=1}^L{\mathrm{G}}_{\sigma }\left(d_i\left(k\right)-{\mathrm{g}}_i(k,\mathbf{x}\left(k\right))\right),$$

(28)

where $d_i\left(k\right)$ is the i-th element of $\mathbf{D}\left(k\right)$, ${\mathrm{g}}_i(k,\mathbf{x}\left(k\right))$ is the i-th row of $\mathrm{g}(k,\mathbf{x}(k\left)\right)$, and $L=n+m$ is the dimension of $\mathbf{D}\left(k\right)$.

Under the MCC, the optimal estimate of $\mathbf{x}\left(k\right)$ can be found by maximizing the cost function

$$\widehat{\mathbf{x}}\left(k\right)=\arg \underset{\mathbf{x}\left(k\right)}{max}{J}_L\left(\mathbf{x}\left(k\right)\right)=\arg \underset{\mathbf{x}\left(k\right)}{max}\sum _{i=1}^L{G}_{\sigma }\left(e_i\left(k\right)\right),$$

(29)

where $e_i\left(k\right)$ is the i-th element of $\mathbf{e}\left(k\right)$ and

$$e_i\left(k\right)=d_i\left(k\right)-{\mathrm{g}}_i\left(k,\mathbf{x}\left(k\right)\right).$$

(30)

Hence, the optimal solution for $\mathbf{x}\left(k\right)$ can be found from the following equation

$$\frac{\partial J_L\left(\mathbf{x}\left(k\right)\right)}{\partial \mathbf{x}\left(k\right)}=0,$$

(31)

which can also be expressed as

$$\sum _{i=1}^L\phi \left(e_i\left(k\right)\right)\frac{\partial e_i\left(k\right)}{\partial x_i\left(k\right)}=0,$$

(32)

when $x_i\left(k\right)$ is the i-th element of $\mathbf{x}\left(k\right)$, and $\phi \left(e_i\left(k\right)\right)={\mathrm{G}}_{\sigma }\left(e_i\left(k\right)\right)·e_i\left(k\right)$. Then, we define $C(k,i)=\phi \left(e_i\left(k\right)\right)/e_i\left(k\right)$ and have

$$C(k,i)={\mathrm{G}}_{\sigma }\left(e_i\left(k\right)\right).$$

(33)

Equation (32) further writes in matrix form as

$${\left(\frac{\partial \mathrm{g}(k,\mathbf{x}(k\left)\right)}{\partial \mathbf{x}\left(k\right)}\right)}^T\mathbf{C}\left(k\right)\left[\mathbf{D}\left(k\right)-\mathrm{g}(k,\mathbf{x}(k\left)\right)\right]=0,$$

(34)

where

$$\mathbf{C}\left(k\right)=\left[\begin{array}{cc}{\mathbf{C}}_x\left(k\right)& \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_y\left(k\right)\end{array}\right],$$

with

$${\mathbf{C}}_x\left(k\right)=diag\left({\mathrm{G}}_{\sigma }\left(e_1\left(k\right)\right),...,{\mathrm{G}}_{\sigma }\left(e_n\left(k\right)\right)\right),{\mathbf{C}}_y\left(k\right)=diag\left({\mathrm{G}}_{\sigma }\left(e_{n+1}\left(k\right)\right),...,{\mathrm{G}}_{n+m}\left(e_{n+m}\left(k\right)\right)\right).$$

Based on Equation (27) with only one iteration in [32], we could obtain similar results in the UKF framework by using $\mathbf{C}\left(k\right)$ to reformulate the measurement information. Similar to [16], there are two ways to do it. One is to re-weight the residual error covariance depending on the value of $e_i\left(k\right)=d_i\left(k\right)-{\mathrm{g}}_i\left(k,\mathbf{x}\left(k\right)\right)$, the other is to reconstruct the ’pseudo-observation’. In general, the derived results based on two ways are the same. In this paper, we just describe the algorithm based on the first way for simplicity.

Defining $\tilde{\mathrm{\Psi }}\left(k\right)$ is the modified covariance, and could be given as

$$\tilde{\mathrm{\Psi }}\left(k\right)=\mathbf{T}\left(k\right){\mathbf{C}}^{-1}\left(k\right){\mathbf{T}}^T\left(k\right).$$

(35)

For our next discussion, we write $\tilde{\mathrm{\Psi }}\left(k\right)$ into two portions so that

$$\tilde{\mathrm{\Psi }}\left(k\right)=\left[\begin{array}{cc}{\tilde{\mathrm{\Psi }}}_x\left(k\right)& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathrm{\Psi }}}_y\left(k\right)\end{array}\right].$$

(36)

In fact, since the true state $\mathbf{x}\left(k\right)$ is unknown, we set $\eta \left(\mathbf{x}\left(k\right)\right)=\widehat{\mathbf{x}}\left(k\right|k-1)-\mathbf{x}\left(k\right)=\mathbf{0}$. In this case, we can see easily that

$$\begin{array}{cc}\hfill {\tilde{\mathrm{\Psi }}}_x\left(k\right)=& {\mathbf{T}}_p\left(k\right|k-1)\ast \mathbf{I}\ast {\mathbf{T}}_p^T\left(k\right|k-1)\hfill \\ \hfill =& \mathbf{P}\left(k\right|k-1).\hfill \end{array}$$

(37)

The modified measurement covariance can be obtained

$$\tilde{\mathbf{R}}\left(k\right)={\tilde{\mathrm{\Psi }}}_y\left(k\right).$$

(38)

In the above derivation, the MCUKF algorithm can be summarized as follows:

• Choose a proper kernel bandwidth σ; set an initial estimate $\widehat{\mathbf{x}}\left(0\right|0)$ and corresponding covariance matrix $\mathbf{P}\left(0\right|0)$; and let $k=1$;
• Use Equations (10) and (12) to calculate the sigma points and the propagated sigma points with respect to function f;
• Compute predicted estimate $\widehat{\mathbf{x}}\left(k\right|k-1)$ and covariance matrix $\mathbf{P}\left(k\right|k-1)$ by Equations (13) and (14), and adopt Cholesky decomposition to get ${\mathbf{T}}_p\left(k\right|k-1)$;
• Use Equations (16) and (17) to calculate the sigma points and the propagated sigma points with respect to function h;
• Use Equations (24)–(38) to derive the modified $\tilde{\mathbf{R}}\left(k\right)$; compute the predicted measurement mean by Equation (18); and replace $\mathbf{R}\left(k\right)$ with $\tilde{\mathbf{R}}\left(k\right)$ in Equation (19) and have the following covariance matrix:

  $${\mathbf{P}}_{\mathbf{yy}}=\sum _{i=0}^{2n}{w}_c^i\left[{\gamma }^i\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right]{\left[{\gamma }^i\left(k\right)-\widehat{\mathbf{y}}\left(k\right)\right]}^T+\tilde{\mathbf{R}}\left(k\right).$$

  (39)
• Use Equations (20)–(23) to compute the filter state mean and covariance matrix; go back to 2.

Remark 1. 

As one can see, different from the UKF algorithm, the MCUKF uses a non-linear regression model combined MCC to update the measurement information. Note that the kernel bandwidth σ is a key parameter in MCUKF. In general, a smaller kernel bandwidth makes the algorithm more robust (with respect to outliers). However, when the kernel bandwidth is too small, it may result in the filter divergence or accuracy deterioration. The reason for this can be seen in [34]. On the other hand, when σ becomes larger and larger, the MCUKF will behave more and more like the UKF algorithm. In particular, the following theorem holds.

Theorem 1. 

If the kernel bandwidth σ approaches infinity, the MCUKF will reduce to the traditional UKF.

Proof of Theorem 1. 

The proof is simple, we give a short explanation. Note that as $\sigma \to \infty$, we have $\mathbf{C}\left(k\right)\to \mathbf{I}$. In this case, the MCUKF is equal to the traditional UKF. ☐

Thus, it is seen that the kernel bandwidth has significant influences on the behavior of MCUKF. In practical applications, the kernel bandwidth can be set manually or optimized by trial and error methods.

4. Illustrative Examples

This section demonstrates the performances of the proposed algorithm by two illustrative examples. In the paper, we mainly compare the estimation performance based on the following benchmarks:

$${\mathrm{MSE}}_1\left(m\right)=\frac{1}{K}\sum _{k=1}^K{(x\left(k\right)-\widehat{x}\left(k\right|k))}^2,m=1,\dots ,M,$$

(40)

$$\mathrm{MSE}=\frac{1}{M}\sum _{m=1}^M{\mathrm{MSE}}_1\left(m\right),$$

(41)

$${\mathrm{MSD}}_{\mathrm{p}}\left(k\right)=\frac{1}{M}\sum _{m=1}^M\left({(x\left(k\right)-\widehat{x}\left(k\right|k))}^2+{(y\left(k\right)-\widehat{y}\left(k\right|k))}^2+{(z\left(k\right)-\widehat{z}\left(k\right|k))}^2\right),k=1,\dots ,K,$$

(42)

$${\mathrm{MSD}}_{\mathrm{v}}\left(k\right)=\frac{1}{M}\sum _{m=1}^M\left({(\dot{x}\left(k\right)-\widehat{\dot{x}}\left(k\right|k))}^2+{(\dot{y}\left(k\right)-\widehat{\dot{y}}\left(k\right|k))}^2+{(\dot{z}\left(k\right)-\widehat{\dot{z}}\left(k\right|k))}^2\right),k=1,\dots ,K,$$

(43)

$${\mathrm{TAMSD}}_{\mathrm{p}}=\frac{1}{k_2-k_1}\sum _{k=k_1+1}^{k_2}{\mathrm{MSD}}_{\mathrm{p}}\left(k\right),$$

(44)

$${\mathrm{TAMSD}}_{\mathrm{v}}=\frac{1}{k_2-k_1}\sum _{k=k_1+1}^{k_2}{\mathrm{MSD}}_{\mathrm{v}}\left(k\right),$$

(45)

where K is the entire time steps in every Monte Carlo run and M represents the total Monte Carlo runs.

4.1. Example 1

Now, a univariate nonstationary growth model (UNGM) that is often used as a benchmark example for nonlinear filtering is considered, whose state and measurement equations are given by

$$x\left(k\right)={\alpha }_{1}x(k-1)+{\alpha }_2\frac{x(k-1)}{1+x{(k-1)}^2}+{\alpha }_3\cos \left(1.2(k-1)\right)+q(k-1),$$

(46)

$$y\left(k\right)=\frac{x{\left(k\right)}^2}{20}+r\left(k\right).$$

(47)

The parameters are set to ${\alpha }_1=0.5,{\alpha }_2=25,{\alpha }_3=8$.

First, we consider the case in which all the noises are Gaussian

$$\begin{array}{c}q(k-1)\sim N(0,1),\hfill \\ r\left(k\right)\sim N(0,1).\hfill \end{array}$$

In this example, the parameters are chosen as $K=500,M=100$.

Table 1. $\mathrm{MSE}$s of x in Gaussian noises.

Filter	MSE of x
UKF	67.6974
MCUKF $(\sigma =2.0)$	87.6836
MCUKF $(\sigma =3.0)$	80.8406
MCUKF $(\sigma =5.0)$	74.0286
MCUKF $(\sigma =10)$	72.3362
MCUKF $(\sigma =20)$	68.6795

illustrates the $\mathrm{MSE}$s of x, defined in Equation (41), in Gaussian noise. Since all the noises are Gaussian, the UKF gives the smallest $\mathrm{MSE}$ in all filters. Moreover, it is noted that when the kernel bandwidth is small, the MCUKF may result in a worse estimation; in contrast, when the kernel bandwidth becomes larger, its performance will approach that of the UKF. Actually, it has been proved easily that when $\sigma \to \infty$, the MCUKF will reduce to the traditional UKF (see Theorem 1). Therefore, one should choose a larger kernel bandwidth in Gaussian noises.

Second, we change the observation noise into a heavy-tailed non-Gaussian noise, with a mixed-Gaussian distribution

$$\begin{array}{c}q(k-1)\sim N(0,1),\hfill \\ r\left(k\right)\sim 0.8N(0,1)+0.2N(0,500).\hfill \end{array}$$

Table 2. $\mathrm{MSE}$s of x in non-Gaussian measurement noise.

Filter	MSE of x
UKF	85.8439
MCUKF $(\sigma =1.0)$	84.1944
MCUKF $(\sigma =2.0)$	82.6933
MCUKF $(\sigma =3.0)$	83.1098
MCUKF $(\sigma =5.0)$	84.7173
MCUKF $(\sigma =10)$	85.4411

shows the $\mathrm{MSE}$s of x in non-Gaussian measurement noise. As one can see, in this case, when kernel bandwidth is too small or too large, the performance of MCUKF will be not good. However, with a proper kernel bandwidth (say $\sigma =2.0$), the MCUKF can outperform the UKF, achieving the smallest $\mathrm{MSE}$. Again, when σ is very large, MCUKF achieves almost the same performance as the UKF.

4.2. Example 2

Finally, we consider a practical example with respect to the relative motion of two spacecrafts, which is illustrated in Figure 1. One of the spacecrafts is called the chief spacecraft, which is moving on the reference orbit, and the other is the deputy spacecraft. They all revolve around the earth and thus the inertial orbital equations of two spacecraft are given as

$${\ddot{\mathbf{r}}}_c=-\frac{\mu }{r_c^3}{\mathbf{r}}_c,$$

(48)

$${\ddot{\mathbf{r}}}_d=-\frac{\mu }{r_d^3}{\mathbf{r}}_d,$$

(49)

where ${\mathbf{r}}_c$ and ${\mathbf{r}}_d$ are the position vectors of the chief spacecraft and the deputy spacecraft in ECI coordinate frame, $r_c$ and $r_d$ are the norms of ${\mathbf{r}}_c$ and ${\mathbf{r}}_d$, respectively, and μ is the gravitational parameter of the earth.

The position vector of the deputy spacecraft relative to the chief spacecraft is

$${\mathbf{\ae }}_{dc}={\mathbf{r}}_d-{\mathbf{r}}_c.$$

(50)

Here, we use the Hill coordinate frame, which is a rectangular, Cartesian, dextral rotating frame centered on the chief spacecraft and refer to it as the local-vertical-local-horizontal (LVLH) frame. Then, the motion of the deputy spacecraft relative to the chief spacecraft in the Hill coordinate frame can be described in [35] as

$$\ddot{x}=2\omega \dot{y}+\dot{\omega }y+{\omega }^{2}x+\frac{\mu }{r_c^2}-\frac{\mu (r_c+x)}{{\left[{(r_c+x)}^2+y^2+z^2\right]}^{3/2}}+q_1,$$

(51)

$$\ddot{y}=-2\omega \dot{x}-\dot{\omega }x+{\omega }^{2}y-\frac{\mu y}{{\left[{(r_c+x)}^2+y^2+z^2\right]}^{3/2}}+q_2,$$

(52)

$$\ddot{z}=-\frac{\mu z}{{\left[{(r_c+x)}^2+y^2+z^2\right]}^{3/2}}+q_3,$$

(53)

with

$${\ddot{r}}_c=r_c{\omega }^2-\frac{\mu }{r_c^2},$$

(54)

$$\dot{\omega }=-\frac{2{\dot{r}}_c\omega }{r_c},$$

(55)

where x, y and z are the radial, in-track and cross-track coordinates of the Hill frame, and ω denotes the orbital angular velocity of the chief spacecraft.

The radar is set on the chief spacecraft to obtain the measurements, and the measurement coordinate system is showed in Figure 2. The measurement equations are given as

$$\rho =\sqrt{(x^2+y^2+z^2)}+r_1,$$

(56)

$$\theta ={\tan }^{-1}\frac{y}{x}+r_2,$$

(57)

$$\varphi ={\tan }^{-1}\frac{z}{\sqrt{x^2+y^2}}+r_3,$$

(58)

where ρ is the relative range between the chief spacecraft and deputy spacecraft, θ is the azimuth angle, and $\varphi$ is the elevation angle.

The two spacecrafts are on the elliptic low Earth orbits. The initial six orbital elements of the chief spacecraft are showed in

Table 3. Initial orbital elements of chief spacecraft.

Orbital Elements	Chief Spacecraft
Semi-major axis	8000 km
Eccentricity	0.150
Orbit inclination	$\pi /6$ rad
Argument of perigee	$\pi /6$ rad
Right ascension of the ascending node	$\pi /18$ rad
True anomaly	0 rad

, and the trajectory of which in ECI coordinate frame is propagated by the fixed-step fourth-order Runge–Kutta algorithm. The prediction of filters is performed every 0.1 s and the measurements record from the radar every 1 s. The state vector $\mathbf{x}\left(k\right)={\left[\begin{array}{cccccc}x\left(k\right)& y\left(k\right)& z\left(k\right)& \dot{x}\left(k\right)& \dot{y}\left(k\right)& \dot{z}\left(k\right)\end{array}\right]}^T$ contains the relative position and velocity components in Hill frame, respectively. The initial true state is

$$\begin{array}{c}x\left(0\right)=31.9262\mathrm{km},\dot{x}\left(0\right)=-0.005583\mathrm{km}/\mathrm{s},\hfill \\ y\left(0\right)=-7.1384\mathrm{km},\dot{y}\left(0\right)=-0.071774\mathrm{km}/\mathrm{s},\hfill \\ z\left(0\right)=33.4729\mathrm{km},\dot{z}\left(0\right)=0.026249\mathrm{km}/\mathrm{s},\hfill \end{array}$$

and the initial estimates and covariance matrix of the states are chosen as

$$\begin{array}{c}\widehat{x}\left(0\right|0)=31.9262\mathrm{km},\widehat{\dot{x}}\left(0\right|0)=-0.004416\mathrm{km}/\mathrm{s},\hfill \\ \widehat{y}\left(0\right|0)=-8.1384\mathrm{km},\widehat{\dot{y}}\left(0\right|0)=-0.061774\mathrm{km}/\mathrm{s},\hfill \\ \widehat{z}\left(0\right|0)=32.4729\mathrm{km},\widehat{\dot{z}}\left(0\right|0)=0.036249\mathrm{km}/\mathrm{s},\hfill \\ \mathbf{P}\left(0\right|0)=diag\left([1,1,1,{0.01}^2,{0.01}^2,{0.01}^2]\right).\hfill \end{array}$$

In this example, 100 independent Monte Carlo runs have been conducted, and, in each case, an elapsed time of 7200 s is considered, that is $K=7200,M=100$. Since correntropy is a local similarity measure, the MCUKF may converge to the optimal solution slowly, especially when the initial values deviate greatly from the true values. Thus, we use the UKF during the first 100 s to make the process settle down and then switch to the MCUKF to continue. The performance of EKF, Huber-EKF (HEKF) [36], UKF and novel robust UKF (NRUKF) [16] are shown for comparison with that of the proposed algorithm.

First, we consider the case in which all the noises are Gaussian, that is,

$$\begin{array}{c}{q}_1\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),r_1\sim N(0,{\left({10}^{-3}\mathrm{km}\right)}^2),\hfill \\ q_2\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),r_2\sim N(0,{(0.05\pi /180)}^2),\hfill \\ q_3\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),r_3\sim N(0,{(0.05\pi /180)}^2).\hfill \end{array}$$

Figure 3 describes the relative motion of the deputy spacecraft in the Hill coordinate frame, Figure 4 and Figure 5 demonstrate the ${\mathrm{MSD}}_{\mathrm{p}}$ and ${\mathrm{MSD}}_{\mathrm{v}}$, defined in Equations (42) and (43), with different filters in Gaussian noises, and

Table 4. ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$ in Gaussian noises.

Filter	${\mathbf{TAMSD}}_{\mathbf{p}}$	${\mathbf{TAMSD}}_{\mathbf{v}}$
EKF	$9.6397\times {10}^{-6}$	$3.7704\times {10}^{-11}$
HEKF	$1.0742\times {10}^{-5}$	$3.9391\times {10}^{-11}$
UKF	$8.4386\times {10}^{-6}$	$2.6552\times {10}^{-11}$
NRUKF	$9.0277\times {10}^{-6}$	$2.7508\times {10}^{-11}$
MCUKF $\left(\sigma =2.0\right)$	$9.3516\times {10}^{-6}$	$2.8038\times {10}^{-11}$
MCUKF $\left(\sigma =5.0\right)$	$8.4909\times {10}^{-6}$	$2.6625\times {10}^{-11}$
MCUKF $\left(\sigma =20\right)$	$8.4403\times {10}^{-6}$	$2.6554\times {10}^{-11}$

summarizes the corresponding ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$ , defined in Equations (44) and (45) (the parameter is set to $k_1=1000,k_2=7200$). Those results illustrate that the UKF play the best performance in all filters in this case and the UKF type filters have the better performance than the EKF type counterparts. One can also observe that the robust filters do not perform as well as their non-robust counterparts in the Gaussian noises. Moreover, it is worth noting that when the kernel bandwidth is small, the MCUKF may achieve a worse performance; whereas, when the kernel bandwidth becomes larger, it has similar results to the UKF.

Second, we consider the case in which the measurement noises are heavy-tailed, with mixed-Gaussian distributions

$$\begin{array}{c}{q}_1\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),\hfill \\ q_2\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),\hfill \\ q_3\sim N(0,{({10}^{-7}\mathrm{km}/{\mathrm{s}}^2)}^2),\hfill \\ r_1\sim 0.9N(0,{({10}^{-3}\mathrm{km})}^2)+0.1N(0,{({10}^{-2}\mathrm{km})}^2),\hfill \\ r_2\sim 0.9N(0,{(0.05\pi /180)}^2)+0.1N(0,{(0.5\pi /180)}^2),\hfill \\ r_3\sim 0.9N(0,{(0.05\pi /180)}^2)+0.1N(0,{(0.5\pi /180)}^2).\hfill \end{array}$$

Figure 6 and Figure 7 reveal the ${\mathrm{MSD}}_{\mathrm{p}}$ and ${\mathrm{MSD}}_{\mathrm{v}}$ with different filters in non-Gaussian noises, and

Table 5. ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$ in non-Gaussian noises.

Filter	${\mathbf{TAMSD}}_{\mathbf{p}}$	${\mathbf{TAMSD}}_{\mathbf{v}}$
EKF	$1.0680\times {10}^{-4}$	$2.2751\times {10}^{-10}$
HEKF	$3.2302\times {10}^{-5}$	$7.5898\times {10}^{-11}$
UKF	$4.2884\times {10}^{-5}$	$8.4705\times {10}^{-11}$
NRUKF	$2.5510\times {10}^{-5}$	$5.7120\times {10}^{-11}$
MCUKF $\left(\sigma =2.0\right)$	$2.4842\times {10}^{-5}$	$5.6557\times {10}^{-11}$
MCUKF $\left(\sigma =5.0\right)$	$3.1545\times {10}^{-5}$	$6.6535\times {10}^{-11}$
MCUKF $\left(\sigma =20\right)$	$4.1800\times {10}^{-5}$	$8.3054\times {10}^{-11}$

lists the corresponding ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$. As one can observe again, the UKF type filters give a smaller ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$ than the EKF type counterparts. All the robust filters are superior to their non-robust counterparts in impulsive noises. It is noted that when the kernel bandwidth is very large, MCUKF achieves almost the same performance as the UKF. However, with a proper kernel bandwidth, the MCUKF can outperform the UKF significantly. Particularly, when $\sigma =2.0$, the MCUKF exhibits the smallest ${\mathrm{TAMSD}}_{\mathrm{p}}$ and ${\mathrm{TAMSD}}_{\mathrm{v}}$.

Moreover, to compare the computational cost, the computation times of different filters in this example are shown in

Table 6. Computation time comparison.

Filter	Computation Ratio
UKF	1
HEKF	$6.81$
NRUKF	$4.62$
MCUKF	$3.91$

. We can see that the computation time of the MCUKF is moderate compared with the UKF, and is superior to that of the HEKF and NRUKF.

5. Conclusions

In this paper, we propose a new unscented Kalman filter (UKF), namely the maximum correntropy unscented Kalman filter (MCUKF), which shows strong robustness against heavy-tailed non-Gaussian noises. The proposed algorithm is recast in the form of a non-linear regression model and makes use of the MCC to obtain the filter estimates. The MCUKF is then applied to the spacecraft relative state estimation, compared with EKF, HEKF, UKF and NRUKF. Simulation results confirm that, with a large kernel bandwidth, the MCUKF will behave like the original UKF. With a proper kernel bandwidth, the new filter can achieve better performance than other filters, especially when the measurement system is disturbed by some impulsive non-Gaussian noises.