Correntropy Based Divided Difference Filtering for the Positioning of Ships

Abstract

In this paper, robust first and second-order divided difference filtering algorithms based on correntropy are proposed, which not only retain the advantages of divided difference filters, but also exhibit robustness in the presence of non-Gaussian noises, especially when the measurements are contaminated by heavy-tailed noises. The proposed filters are then applied to the problem of ship positioning. In order to improve the accuracy and reliability of ship positioning, the positioning method combines the Dead Reckoning (DR) algorithm and the Global Positioning System (GPS). Experimental results of an illustrative example show the superior performance of the new algorithms when applied to ship positioning.

1. Introduction

With the development of microcomputers and electronic integration technology, the requirements for accuracy and reliability of ship positioning are becoming more and more complex. As a key node in the air-space-ground of integrated information networks, ships can obtain the positioning information from many systems. The Dead Reckoning (DR) is a common technology used in ship positioning, which calculates the position of ship in real time based on the speed of ship and direction information, it also considers the influence of wind and flow. However, errors in the speed and direction information, and incomplete compensations of wind and flow cause positioning errors of DR which accumulate with time [1,2]. Global Positioning System (GPS) is a satellite based system, which provides the accurate velocity and position of a ship by making use of the GPS receiver. But GPS has limitations such as a low sampling rate as well as being susceptible to interference [3,4,5]. The integration of DR and GPS takes advantage of two techniques, this integration system has a better performance than the single techniques in terms of accuracy and reliability [6,7,8].

Filtering is a key problem in the integrated positioning system. The Kalman filter (KF) [9,10] is a well known method to estimate the state of linear systems. However, the model of DR/GPS integrated positioning system is nonlinear. To solve the nonlinear filtering problem, extensions of the KF using some approximations have been proposed. The extended Kalman filter (EKF) [11,12,13] approximates the nonlinear equation by its first-order linearization. The unscented Kalman filter (UKF) [13,14,15,16] approximates the probability distribution of the state by a set of deterministically chosen sigma points and propagates the distribution through the nonlinear equation, which provides higher-order approximation than the EKF. Nevertheless, the parameters used in the UKF have a great effect on performance of the algorithm. If they are not tuned finely, it is easy to face the problem of numerical instability in practical applications due to the propagation of the non-positive definite covariance matrix. Another effective, alternative way is the divided difference filter (DDF). The DDF [17,18] is derived from polynomial approximations of the nonlinear transformations using multidimensional Stirling’s interpolation formula and can be classified into the first-order divided difference (DD1) filter and the second-order divided difference (DD2) filter. The DDF can guarantee a positive semi-definiteness of the covariance matrix. The DD2 filter not only has fewer parameters, but also has nearly the same performance as the UKF. Therefore, we use the DDF type filters for ship positioning. The traditional DDF is suitable for the Gaussian noise. However, in this real application, the measuring instruments can be affected by extreme sea environments. They sometimes break down, or suffer from operator error, which cause the measurement noise to be of the heavy-tailed non-Gaussian form. In these cases, the traditional DDF, which is based on the minimum mean square error (MMSE) criterion, cannot play well because the MMSE criterion is very sensitive to the heavy-tailed non-Gaussian noise [19].

To solve the non-Gaussian filtering problems, some robust algorithms exist. The Huber methodology is a method that combines minimum ${\ell }_1$ and ${\ell }_2$-norm [20,21,22]. The Student t technique assumes that the process and measurement noises obey Student t distribution [23,24]. Another effective approach is the information theoretic learning (ITL). In particular, the correntropy, which is one of the ITL measures, can capture high-order statistics of the data rather than the common second-order statistics and has been widely used in many fields [25,26,27,28,29,30,31,32]. In recent years, some correntropy-based Kalman filterings have been proposed [33,34,35], which are mainly applied to linear models. The correntropy-based unscented Kalman filters can solve some nonlinear problems [36,37], but it is easy to have problems with numerical instability for integrated positioning.

In this paper, two novel nonlinear filters, the correntropy-based first-order divided difference (CDD1) filter and the correntropy-based second-order divided difference (CDD2) filter, are proposed to solve the problem of ship positioning. The proposed algorithms not only retain the advantages of DDF algorithms, but also exhibit the robustness in the presence of heavy-tailed non-Gaussian noise. Different from the works [38] which adopts the linear regression model, the proposed algorithms adopt the nonlinear regression model.

The rest of the paper is organized as follows. Section 2 provides a short review of the correntropy, the DD1, and DD2 filters. In Section 3, the the CDD1 and CDD2 filters are derived. In Section 4, the example of ship positioning shows the performance of the proposed algorithms. The conclusion is given in Section 5.

2. Preliminaries

2.1. Correntropy

The correntropy is an important measure in ITL, which measures the similarity between two random variables $X\in \mathbb{R}$ and $Y\in \mathbb{R}$. Given the joint distribution function of X and Y, $F_{XY}(x,y)$, the correntropy is defined by

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

(1)

where $E[·]$ represents the expectation operator, and $\kappa (·,·)$ is a shift-invariant Mercer kernel. In this paper, the Gaussian kernel is chosen as the kernel function:

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

(2)

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

In most practical applications, only a limited number of data samples are available and the joint distribution $F_{XY}$ is usually unknown. In this case, we often use the sample mean estimator to estimate the correntropy:

$$\widehat{V}(X,Y)=\frac{1}{N}\sum _{i=1}^N{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 $F_{XY}$.

Taking the Taylor series expansion for the Gaussian kernel yields

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

(4)

Note that the correntropy is the weighted sum of all even order moments of the error variable $X-Y$ and the kernel bandwidth appears as a parameter to weight the second-order and higher-order moments. In particular, with a very large kernel bandwidth, the second-order moment will play a major role.

2.2. DD1 Filter

In this paper, the nonlinear state and measurement equations are described by

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

(5)

$${\mathbf{y}}_k=h\left({\mathbf{x}}_k,{\mathbf{r}}_k\right)$$

(6)

where $\mathbf{x}\left(k\right)\in \mathbb{R}{}^n$ and $\mathbf{y}\left(k\right)\in \mathbb{R}{}^m$ are the n-dimensional state vector and m-dimensional measurement vector at time step k. $f\left(·\right)$ and $h\left(·\right)$ are the continuous and differentiable state functions and the measurement function. ${\mathbf{q}}_{k-1}$ and ${\mathbf{r}}_k$ are the process and measurement noises, which are assumed i.i.d and independent of states, and with means denoted by ${\overline{\mathbf{q}}}_{k-1}$ and ${\overline{\mathbf{r}}}_k$ and covariance matrices denoted by ${\mathbf{Q}}_{k-1}$ and ${\mathbf{R}}_k$.

The square-root decompositions of the predicted state error covariance ${\overline{\mathbf{P}}}_k$, update state error covariance ${\widehat{\mathbf{P}}}_k$, process noise covariance ${\mathbf{Q}}_k$, and measurement noise covariance ${\mathbf{R}}_k$ are introduced as

$${\overline{\mathbf{P}}}_k={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}^T$$

(7)

$${\widehat{\mathbf{P}}}_k={\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}{\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}^T$$

(8)

$${\mathbf{Q}}_k={\mathbf{S}}_{{\mathbf{q}}_k}{\mathbf{S}}_{{\mathbf{q}}_k}^T$$

(9)

$${\mathbf{R}}_k={\mathbf{S}}_{{\mathbf{r}}_k}{\mathbf{S}}_{{\mathbf{r}}_k}^T$$

(10)

The DD1 filter uses the first-order divided differences to approximate. Let the element of i-th row, j-th column of ${{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_k}$ be denoted as ${{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k}i,j}$, i.e., ${{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_k}=\left\{{{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k}i,j}\right\}$, and similarly for other matrices. The four matrices containing first-order divided differences are defined as

$${{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k-1}}=\left\{\left[f_i\left({\widehat{\mathbf{x}}}_{k-1}+c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)-f_i\left({\widehat{\mathbf{x}}}_{k-1}-c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)\right]/\phantom{\left[f_i\left({\widehat{\mathbf{x}}}_{k-1}+c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)-f_i\left({\widehat{\mathbf{x}}}_{k-1}-c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)\right]2c}2c\right\}$$

(11)

$${{\mathbf{S}}^{\prime }}_{\mathbf{x}{\mathbf{q}}_{k-1}}=\left\{\left[f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}+c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)-f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}-c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)\right]/\phantom{\left[f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}+c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)-f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}-c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)\right]2c}2c\right\}$$

(12)

$${{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}=\left\{\left[h_i\left({\overline{\mathbf{x}}}_k+c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)-h_i\left({\overline{\mathbf{x}}}_k-c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)\right]/\phantom{\left[h_i\left({\overline{\mathbf{x}}}_k+c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)-h_i\left({\overline{\mathbf{x}}}_k-c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)\right]2c}2c\right\}$$

(13)

$${{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}=\left\{\left[h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k+c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)-h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k-c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)\right]/\phantom{\left[h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k+c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)-h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k-c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)\right]2c}2c\right\}$$

(14)

where ${\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j}$ is the j-th column of ${\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1}}$, c is the interval length, and is generally set as $c^2=3$.

The predicted state and the Cholesky factor of corresponding covariance are given by

$${\overline{\mathbf{x}}}_k=f\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}\right)$$

(15)

$${\overline{\mathbf{S}}}_{{\mathbf{x}}_k}=\mathrm{H}\left(\left[\begin{array}{cc}{{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k-1}}& {{\mathbf{S}}^{\prime }}_{\mathbf{x}{\mathbf{q}}_{k-1}}\end{array}\right]\right)$$

(16)

where $\mathrm{H}\left(·\right)$ denotes a Householder transformation of the augment matrix [17].

The predicted measurement and the Cholesky factor of corresponding covariance are given by

$${\overline{\mathbf{y}}}_k=h\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k\right)$$

(17)

$${\mathbf{S}}_{{\mathbf{y}}_k}=\mathrm{H}\left(\left[\begin{array}{cc}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}\end{array}\right]\right)$$

(18)

The Kalman gain is computed as

$${\mathbf{K}}_k={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}^T{\left({\mathbf{S}}_{{\mathbf{y}}_k}{\mathbf{S}}_{{\mathbf{y}}_k}^T\right)}^{-1}$$

(19)

The updated state and the Cholesky factor of corresponding covariance are given by

$${\widehat{\mathbf{x}}}_k={\overline{\mathbf{x}}}_k+{\mathbf{K}}_k\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)$$

(20)

$${\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}=\mathrm{H}\left(\left[\begin{array}{cc}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}-{\mathbf{K}}_k{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {\mathbf{K}}_k{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}\end{array}\right]\right)$$

(21)

2.3. DD2 Filter

The DD2 filter uses the second-order divided differences to approximate. Four matrices containing second-order divided differences are defined as

$$\begin{array}{cc}\hfill {{\mathbf{S}}^{″}}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k-1}}=& \left\{\frac{\sqrt{c^2-1}}{2c^2}\left[f_i\left({\widehat{\mathbf{x}}}_{k-1}+c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)+f_i\left({\widehat{\mathbf{x}}}_{k-1}-c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)\right.\right.\hfill \\ & \left.\left.-2f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}\right)\right]\right\}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {{\mathbf{S}}^{″}}_{\mathbf{x}{\mathbf{q}}_{k-1}}=& \left\{\frac{\sqrt{c^2-1}}{2c^2}\right.\left[f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}+c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)\right.+f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}-c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)\hfill \\ & \left.\left.-2f_i\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}\right)\right]\right\}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {{\mathbf{S}}^{″}}_{\mathbf{y}{\overline{\mathbf{x}}}_k}=& \left\{\frac{\sqrt{c^2-1}}{2c^2}\left[h_i\left({\overline{\mathbf{x}}}_k+c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)\right.\right.+h_i\left({\overline{\mathbf{x}}}_k-c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)\hfill \\ & \left.\left.-2h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k\right)\right]\right\}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {{\mathbf{S}}^{″}}_{\mathbf{y}{\mathbf{r}}_k}=& \left\{\frac{\sqrt{c^2-1}}{2c^2}\left[h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k+c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)\right.\right.+h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k-c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)\hfill \\ & \left.\left.-2h_i\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k\right)\right]\right\}\hfill \end{array}$$

(25)

The predicted state and the Cholesky factor of corresponding covariance are given by

$$\begin{array}{cc}\hfill {\overline{\mathbf{x}}}_k=& \left(\frac{c^2-n_x-n_q}{c^2}\right)f\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}\right)\hfill \\ & +\frac{1}{2c^2}\sum _{j=1}^{n_x}\left[f\left({\widehat{\mathbf{x}}}_{k-1}+c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)+f\left({\widehat{\mathbf{x}}}_{k-1}-c{\widehat{\mathbf{S}}}_{{\mathbf{x}}_{k-1},j},{\overline{\mathbf{q}}}_{k-1}\right)\right]\hfill \\ & +\frac{1}{2c^2}\sum _{j=1}^{n_q}\left[f\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}+c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)+f\left({\widehat{\mathbf{x}}}_{k-1},{\overline{\mathbf{q}}}_{k-1}-c{\mathbf{S}}_{{\mathbf{q}}_{k-1},j}\right)\right]\hfill \end{array}$$

(26)

$${\overline{\mathbf{S}}}_{{\mathbf{x}}_k}=\mathrm{H}\left(\left[\begin{array}{cccc}{{\mathbf{S}}^{\prime }}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k-1}}& {{\mathbf{S}}^{\prime }}_{\mathbf{x}{\mathbf{q}}_{k-1}}& {{\mathbf{S}}^{″}}_{\mathbf{x}{\widehat{\mathbf{x}}}_{k-1}}& {{\mathbf{S}}^{″}}_{\mathbf{x}{\mathbf{q}}_{k-1}}\end{array}\right]\right)$$

(27)

The predicted measurement and the Cholesky factor of corresponding covariance are given by

$$\begin{array}{cc}\hfill {\overline{\mathbf{y}}}_k=& \left(\frac{c^2-n_x-n_r}{c^2}\right)h\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k\right)\hfill \\ & +\frac{1}{2c^2}\sum _{j=1}^{n_x}\left[h\left({\overline{\mathbf{x}}}_k+c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)+h\left({\overline{\mathbf{x}}}_k-c{\overline{\mathbf{S}}}_{{\mathbf{x}}_k,j},{\overline{\mathbf{r}}}_k\right)\right]\hfill \\ & +\frac{1}{2c^2}\sum _{j=1}^{n_r}\left[h\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k+c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)+h\left({\overline{\mathbf{x}}}_k,{\overline{\mathbf{r}}}_k-c{\mathbf{S}}_{{\mathbf{r}}_k,j}\right)\right]\hfill \end{array}$$

(28)

$${\mathbf{S}}_{{\mathbf{y}}_k}=\mathrm{H}\left(\left[\begin{array}{cccc}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}& {{\mathbf{S}}^{″}}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {{\mathbf{S}}^{″}}_{\mathbf{y}{\mathbf{r}}_k}\end{array}\right]\right)$$

(29)

The Kalman gain is computed as

$${\mathbf{K}}_k={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}^T{\left({\mathbf{S}}_{{\mathbf{y}}_k}{\mathbf{S}}_{{\mathbf{y}}_k}^T\right)}^{-1}$$

(30)

The updated state and the Cholesky factor of corresponding covariance are given by

$${\widehat{\mathbf{x}}}_k={\overline{\mathbf{x}}}_k+{\mathbf{K}}_k\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)$$

(31)

$${\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}=\mathrm{H}\left(\left[\begin{array}{cccc}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}-{\mathbf{K}}_k{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {\mathbf{K}}_k{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}& {\mathbf{K}}_k{{\mathbf{S}}^{″}}_{\mathbf{y}{\overline{\mathbf{x}}}_k}& {\mathbf{K}}_k{{\mathbf{S}}^{″}}_{\mathbf{y}{\mathbf{r}}_k}\end{array}\right]\right)$$

(32)

3. Correntropy-Based DDF

This section derives the DD1 and DD2 filters under maximum correntropy criterion (MCC). First, the measurement Equation (6), whose noise is additive, would be written as ${\mathbf{y}}_k=h\left({\mathbf{x}}_k\right)+{\mathbf{r}}_k$, and the predicted state and the Choleskey factor of corresponding covariance are obtained by Equations (15) and (16). Assuming the true state is denoted as ${\mathbf{x}}_k$, the predicted state error is written as ${\xi }_k={\mathbf{x}}_k-{\overline{\mathbf{x}}}_k$. Then, a nonlinear regression model is reconstructed as follows:

$$\left[\begin{array}{c}{\overline{\mathbf{x}}}_k\\ {\mathbf{y}}_k\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}_k\\ h\left({\mathbf{x}}_k\right)\end{array}\right]+{\delta }_k$$

(33)

where ${\delta }_k=\left[\begin{array}{c}-{\xi }_k\\ {\mathbf{r}}_k\end{array}\right]$, with the Cholesky factor of covariance of ${\delta }_k$ is given by

$${\mathbf{D}}_k=\left[\begin{array}{cc}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}& \mathbf{0}\\ \mathbf{0}& {\mathbf{S}}_{{\mathbf{r}}_k}\end{array}\right]$$

(34)

Left multiplying both sides of Equation (33) by ${\mathbf{D}}_k^{-1}$, we have the following regression model:

$${\mathbf{z}}_k=m\left({\mathbf{x}}_k\right)+{\mathbf{e}}_k$$

(35)

where ${\mathbf{z}}_k={\mathbf{D}}_k^{-1}\left[\begin{array}{c}{\overline{\mathbf{x}}}_k\\ {\mathbf{y}}_k\end{array}\right]$, $m\left({\mathbf{x}}_k\right)={\mathbf{D}}_k^{-1}\left[\begin{array}{c}{\mathbf{x}}_k\\ h\left({\mathbf{x}}_k\right)\end{array}\right]$ and ${\mathbf{e}}_k={\mathbf{D}}_k^{-1}{\delta }_k$.

The problem associated with Equation (35) can be solved by making use of the MCC, and the corresponding cost function is given by

$$J\left({\mathbf{x}}_k\right)=\sum _{i=1}^L{G}_{\sigma }\left(e_{k,i}\right)$$

(36)

where $e_{k,i}$ is the i-th component of the vector ${\mathbf{e}}_k={\mathbf{z}}_k-m\left({\mathbf{x}}_k\right)$, and L is the dimension of ${\mathbf{e}}_k$.

Then, the solution of the problem mentioned previously can be found by setting the first derivation of Equation (36) to be zero:

$$\sum _{i=1}^L\varphi \left(e_{k,i}\right)\frac{\partial e_{k,i}}{\partial {\mathbf{x}}_k}=\mathbf{0}$$

(37)

where $\varphi \left(e_{k,i}\right)=G_{\sigma }\left(e_{k,i}\right)·e_{k,i}$. By defining $C\left(e_{k,i}\right)=\varphi \left(e_{k,i}\right)/\phantom{\varphi \left(e_{k,i}\right)e_{k,i}=G_{\sigma }}{e}_{k,i}=G_{\sigma }\left(e_{k,i}\right)$ and $\mathbf{C}=diag\left[C\left(e_{k,i}\right)\right]=\left[\begin{array}{cc}{\mathbf{C}}_{\mathbf{x}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_{\mathbf{y}}\end{array}\right]$, Equation (37) can be written in matrix form as

$${\left(\frac{\partial m\left({\mathbf{x}}_k\right)}{\partial {\mathbf{x}}_k}\right)}^T\mathbf{C}\left(m\left({\mathbf{x}}_k\right)-{\mathbf{z}}_k\right)=\mathbf{0}$$

(38)

In fact, the MCC uses $\mathbf{C}$ to re-weight the covariance matrix of the residual error ${\mathbf{e}}_k$ and reconstruct the measurement information. Thus, the updated covariance can be written as

$$\tilde{\mathsf{\Psi }}={\mathbf{D}}_k{\mathbf{C}}^{-1}{\mathbf{D}}_k^T$$

(39)

and decomposed into two portion so that

$$\tilde{\mathsf{\Psi }}=\left[\begin{array}{cc}{\tilde{\mathsf{\Psi }}}_{\mathbf{x}}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathsf{\Psi }}}_{\mathbf{y}}\end{array}\right]$$

(40)

The initial value can be set as ${\widehat{\mathbf{x}}}_k^{\left(0\right)}={\overline{\mathbf{x}}}_k$ or be equal to the updated state computed from the corresponding standard DD1 or DD2. Then, we can obtain the following equations

$${\tilde{\mathsf{\Psi }}}_{\mathbf{x}}={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}·{\mathbf{C}}_{\mathbf{x}}^{-1}·{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}^T$$

(41)

$${\tilde{\mathsf{\Psi }}}_{\mathbf{y}}={\mathbf{S}}_{{\mathbf{r}}_k}·{\mathbf{C}}_{\mathbf{y}}^{-1}·{\mathbf{S}}_{{\mathbf{r}}_k}^T$$

(42)

It is noted that the aforementioned derivation is the extension of Equation (25) in Reference [33]. To reduce the computation time, only one iteration would work to obtain the solution.

Then, a one-step correntropy update for the DD1 filter can be written as

$${\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)}=\mathrm{H}\left(\left[\begin{array}{cc}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}\end{array}\right]\right)$$

(43)

$${\mathbf{K}}_k^{\left(1\right)}={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}^T{\left({\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)}{\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)T}\right)}^{-1}$$

(44)

$${\widehat{\mathbf{x}}}_k={\overline{\mathbf{x}}}_k+{\mathbf{K}}_k^{\left(1\right)}\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)$$

(45)

$${\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}=\mathrm{H}\left(\left[\begin{array}{cc}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}-{\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}\end{array}\right]\right)$$

(46)

Similarly, a one-step correntropy update for the DD2 filter can be written as

$${\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)}=\mathrm{H}\left(\left[\begin{array}{cccc}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}& {{\mathbf{S}}^{″}}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {{\mathbf{S}}^{″}}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}\end{array}\right]\right)$$

(47)

$${\mathbf{K}}_k^{\left(1\right)}={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}^T{\left({\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)}{\mathbf{S}}_{{\mathbf{y}}_k}^{\left(1\right)T}\right)}^{-1}$$

(48)

$${\widehat{\mathbf{x}}}_k={\overline{\mathbf{x}}}_k+{\mathbf{K}}_k^{\left(1\right)}\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)$$

(49)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}=& \mathrm{H}\left(\left[\begin{array}{cc}{\overline{\mathbf{S}}}_{{\mathbf{x}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}-{\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{\prime }}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}\end{array}\right.\right.\hfill \\ & \left.\left.\begin{array}{cc}{\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{″}}_{\mathbf{y}{\overline{\mathbf{x}}}_k}{\mathbf{C}}_{\mathbf{x}}^{-1/2}& {\mathbf{K}}_k^{\left(1\right)}{{\mathbf{S}}^{″}}_{\mathbf{y}{\mathbf{r}}_k}{\mathbf{C}}_{\mathbf{y}}^{-1/2}\end{array}\right]\right)\hfill \end{array}$$

(50)

Remark 1.

The proposed correntropy-based divided difference filters utilize the correntropy theory to improve the performance in the presence of heavy-tailed non-Gaussian noises, in which a nonlinear regression problem is introduced to update the measurement information. It is noted that the kernel bandwidth plays a key role in the proposed algorithm. In general, a smaller kernel bandwidth exhibits more robust properties of the correntropy. However, when the kernel bandwidth is too small, it will lead to an accretion of estimation error or even filtering divergence. A sufficient condition for guaranteeing the convergence of filter was introduced in Reference [39]. Moreover, when kernel bandwidth $\sigma \to \infty$, the matrix $\mathbf{C}\to \mathbf{I}$, and the CDD1 and CDD2 filters would reduce to the original DD1 and DD2 filters. The choice of the kernel bandwidth in practical applications is discussed in next section.

In practical applications, the measurement system may sometimes obtain extremely large measurements. In this case, the CDD1 and CDD2 filters may face numerical problems since ${\mathbf{C}}_{\mathbf{y}}$ will be nearly singular. In view of this problem, a method is introduced as

$$a_k={\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)}^T{\left({\mathbf{S}}_{{\mathbf{y}}_k}{\mathbf{S}}_{{\mathbf{y}}_k}^T\right)}^{-1}\left({\mathbf{y}}_k-{\overline{\mathbf{y}}}_k\right)$$

(51)

If $\left|a_k\right|>\theta$, with $\theta$ being a positive threshold, only the predicted step is worked, that is ${\widehat{\mathbf{x}}}_k={\overline{\mathbf{x}}}_k$, ${\widehat{\mathbf{S}}}_{{\mathbf{x}}_k}={\overline{\mathbf{S}}}_{{\mathbf{x}}_k}$. If $\left|a_k\right|⩽\theta$, the entire steps are worked.

The flow of the CDD1 filter algorithm is shown in Figure 1. Since the flow of the CDD2 filter is similar, we omit it.

4. Positioning of Ships

In this section, to demonstrate the performance of the proposed algorithms, we apply them to solve the problem of ship positioning. The model combines DR and GPS technology to improve the accuracy of positioning.

4.1. The State and Measurement Models

The motion of a ship can be denoted as a nonlinear function in terms of many factors, such as speed, course, shape of earth, sea current, wind, and so on. There are two kinds of maneuvering motions of a ship, which are speed maneuver in a straight line and direction maneuver. Since the acceleration of the ship is generally small and the sampling period of the integrated positioning system is relatively short, the speed maneuver can be ignored or regarded as the speed noise. The direction maneuver can be approximately described by uniform circular motion with a constant change rate of the ship’s course. Therefore, there are two kinds of states related to the motion of a ship, which are uniform linear motion and uniform circular motion. In view of the above, the state vector is chosen as

$$\mathbf{x}={\left[\begin{array}{ccccccc}\phi & \lambda & v_N& v_E& s& K& \mathsf{\Omega }\end{array}\right]}^T$$

where $\phi$ and $\lambda$ denote the arc lengths of latitude and longitude, $v_N$ and $v_E$ are the northward velocity component and the eastward velocity component of the ocean current, s denotes the velocity of that ship relative to the water, K is the ship’s course, and $\mathsf{\Omega }$ represents the change rate of the ship’s course. Correspondingly, the state equations of the ship are written as [40]

$$\dot{\phi }=v_N+scosK+q_1$$

(52)

$$\dot{\lambda }=v_E+ssinK+q_2$$

(53)

$${\dot{v}}_N=-\beta v_N+q_3$$

(54)

$${\dot{v}}_E=-\beta v_E+q_4$$

(55)

$$\dot{s}=q_5$$

(56)

$$\dot{K}=\mathsf{\Omega }+q_6$$

(57)

$$\dot{\mathsf{\Omega }}=q_7$$

(58)

where $\beta$ denotes the inverse correlation time of ocean current, $q_i$ are independent Gaussian white noises. By discretizing these equations, we obtain the following equations

$$\begin{array}{cc}\hfill {\phi }_k=& {\phi }_{k-1}+{\beta }^{-1}\left(1-exp(-\beta T)\right)v_{N,k-1}+\hfill \\ & s_{k-1}cos\left(K_{k-1}+0.5T{\mathsf{\Omega }}_{k-1}\right)T+q_{1,k-1}\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill {\lambda }_k=& {\lambda }_{k-1}+{\beta }^{-1}\left(1-exp(-\beta T)\right)v_{E,k-1}+\hfill \\ & s_{k-1}sin\left(K_{k-1}+0.5T{\mathsf{\Omega }}_{k-1}\right)T+q_{2,k-1}\hfill \end{array}$$

(60)

$$v_{N,k}=exp(-\beta T)v_{N,k-1}+q_{3,k-1}$$

(61)

$$v_{E,k}=exp(-\beta T)v_{E,k-1}+q_{4,k-1}$$

(62)

$$s_k=s_{k-1}+q_5$$

(63)

$$K_k=K_{k-1}+{\mathsf{\Omega }}_{k-1}T+q_{6,k-1}$$

(64)

$${\mathsf{\Omega }}_k={\mathsf{\Omega }}_{k-1}+q_{7,k-1}$$

(65)

where T is sampling period.

Then, the measurement vector is chosen as

$${\mathbf{y}}_k={\left[\begin{array}{cccc}{\phi }_{G,k}& {\lambda }_{G,k}& s_{L,k}& K_{G,k}\end{array}\right]}^T$$

where ${\phi }_G$ and ${\lambda }_G$ denote the arc lengths of latitude and longitude provided by GPS, $s_L$ is the velocity of that ship relative to the water provided by electromagnetic log, $K_G$ represents the ship’s course provided by electric gyrocompass. Accordingly, the measurement equations are given as

$${\phi }_{G,k}={\phi }_k+{\eta }_{1,k}$$

(66)

$${\lambda }_{G,k}={\lambda }_k+{\eta }_{2,k}$$

(67)

$$s_{L,k}=s_k+{\eta }_{3,k}$$

(68)

$$K_{G,k}=K_k+{\eta }_{4,k}$$

(69)

where ${\eta }_{i,k}$ are independent white noises.

Some filters are used for comparison, including the EKF, DD1, UKF, DD2, the Huber-based first-order divided difference (HDD1) filter and the Huber-based second-order divided difference (HDD2) filter, and the HDD1 and HDD2 adopt the Huber methodology. The following benchmarks are used to show the estimation performance:

$${\mathrm{MSE}}_k=\frac{1}{M}\sum _{k=1}^M{\left(x_k-{\widehat{x}}_k\right)}^2$$

(70)

$$\mathrm{TMSE}=\frac{1}{k_2-k_1+1}\sum _{k=k_1}^{k_2}{\mathrm{MSE}}_k$$

(71)

In this simulation, 100 independent Monte Carlo experiments were conducted, and the lasting time of each experiment was 1200 s. The parameters of this example and initial conditions are summarized in

Table 1. Parameters of the example.

Parameter	Value
${\beta }^{-1},\mathrm{s}$	27,780/s
$T,\mathrm{s}$	12

and

Table 2. Initial conditions.

Parameter	Value
${\phi }_0,\mathrm{m}$	$2.2239\times {10}^6$
${\lambda }_0,\mathrm{m}$	$1.2565\times {10}^7$
$v_{N,0},\mathrm{m}/\mathrm{s}$	1
$v_{E,0},\mathrm{m}/\mathrm{s}$	1
$s_0,\mathrm{m}/\mathrm{s}$	10.289
$K_0,\mathrm{rad}$	$\pi /4$
${\mathsf{\Omega }}_0,\mathrm{rad}/\mathrm{s}$	0

, and the process noises satisfy the Gaussian distributions:

$$\begin{array}{c}{q}_1\sim \mathrm{N}(0,0.684{\mathrm{m}}^2),q_2\sim \mathrm{N}(0,0.684{\mathrm{m}}^2)\hfill \\ q_3\sim \mathrm{N}(0,0.000158{(\mathrm{m}/\mathrm{s})}^2),q_4\sim \mathrm{N}(0,0.000158{(\mathrm{m}/\mathrm{s})}^2)\hfill \\ q_5\sim \mathrm{N}(0,0.00158{(\mathrm{m}/\mathrm{s})}^2),q_6\sim \mathrm{N}(0,0.0026\mathrm{ra}{\mathrm{d}}^2)\hfill \\ q_7\sim 0\hfill \end{array}$$

4.2. Simulation Results

First, we assume the measurement noises satisfy the Gaussian distributions:

$$\begin{array}{c}{r}_1\sim \mathrm{N}(0,10000{\mathrm{m}}^2),r_2\sim \mathrm{N}(0,10000{\mathrm{m}}^2)\hfill \\ r_3\sim \mathrm{N}(0,0.0423{(\mathrm{m}/\mathrm{s})}^2),r_4\sim \mathrm{N}(0,0.0000395\mathrm{ra}{\mathrm{d}}^2)\hfill \end{array}$$

$\mathrm{TMSE}$s of position in Gaussian noises are revealed in

Table 3. $\mathrm{TMSE}$s of position in Gaussian noises.

Filter	$\mathbf{TMSE}$sof$\mathit{\phi }$	$\mathbf{TMSE}$sof$\mathit{\lambda }$
EKF	10.3197	11.1708
DD1	10.3196	11.1707
HDD1	12.1968	13.1074
CDD1 $\left(\sigma =2\right)$	23.8353	22.6151
CDD1 $\left(\sigma =5\right)$	10.4517	11.2693
CDD1 $\left(\sigma =20\right)$	10.3215	11.1690
UKF	10.3093	11.1557
DD2	10.3086	11.1549
HDD2	12.1613	13.0604
CDD2 $\left(\sigma =2\right)$	23.8503	22.3792
CDD2 $\left(\sigma =5\right)$	10.4366	11.2484
CDD2 $\left(\sigma =20\right)$	10.3102	11.1529

. It can be seen that the DD2 filter and UKF have a similar performance, likewise for the DD1 filter and EKF. The DD2 filter and UKF are superior to the DD1 filter and EKF, exhibiting smaller errors. Meanwhile, the robust filters do not perform as well as their non-robust counterparts under Gaussian noise conditions. Moreover, the CDD1 and CDD2 filters achieve almost the same performance as the DD1 and DD2 filters when the kernel bandwidth is large enough. It is noted that the UKF sometimes stops executing because the parameters of the UKF may not be finely tuned enough to bring about the problem of propagation of the non-positive definite covariance matrix.

Second, we assume the measurement noises are heavy-tailed non-Gaussian, satisfying the following distributions:

$$\begin{array}{c}{r}_1\sim \mathrm{N}(0,10000{\mathrm{m}}^2)+\mathrm{N}(0,1000000{\mathrm{m}}^2)\hfill \\ r_2\sim \mathrm{N}(0,10000{\mathrm{m}}^2)+\mathrm{N}(0,1000000{\mathrm{m}}^2)\hfill \\ r_3\sim \mathrm{N}(0,0.0423{(\mathrm{m}/\mathrm{s})}^2)+\mathrm{N}(0,4.23{(\mathrm{m}/\mathrm{s})}^2)\hfill \\ r_4\sim \mathrm{N}(0,0.0000395\mathrm{ra}{\mathrm{d}}^2)+\mathrm{N}(0,0.00395\mathrm{ra}{\mathrm{d}}^2)\hfill \end{array}$$

Figure 2, Figure 3, Figure 4 and Figure 5 show the $\mathrm{MSE}$s of different filters for the non-Gaussian case, and

Table 4. $\mathrm{TMSE}$s of position in non-Gaussian noises.

Filter	$\mathbf{TMSE}$sof$\mathit{\phi }$	$\mathbf{TMSE}$sof$\mathit{\lambda }$
EKF	91.7839	84.7678
DD1	91.7776	84.7557
HDD1	44.2747	37.6608
CDD1 $\left(\sigma =2\right)$	39.4284	33.0989
CDD1 $\left(\sigma =5\right)$	64.5026	56.9180
CDD1 $\left(\sigma =20\right)$	89.3531	82.3220
UKF	91.4480	84.5617
DD2	91.4269	84.5294
HDD2	43.8700	37.4722
CDD2 $\left(\sigma =2\right)$	39.0126	32.8825
CDD2 $\left(\sigma =5\right)$	64.1485	56.7222
CDD2 $\left(\sigma =20\right)$	89.0010	82.0980

summarizes the corresponding $\mathrm{TMSE}$s. The performances of the DD2, UKF, DD1, and EKF follow the behavior from the Gaussian case. In the non-Gaussian case, the robust filters outperform the corresponding non-robust filters. With a very large kernel bandwidth, the CDD1 and CDD2 filters achieve a similar estimation to the DD1 and DD2 filters; with a proper kernel bandwidth, the CDD1 and CDD2 give smaller errors than the non-robust filters; in particular, when the kernel bandwidth is set to 2, the CDD2 exhibits the smallest errors.

5. Conclusions

This paper proposes two correntropy-based divided difference filtering methods, namely CDD1 and CDD2, which show strong robustness against heavy-tailed non-Gaussian noises. The proposed algorithms recast the nonlinear regression models and use the maximum correntropy criterion to obtain the solution. The two robust filters are then applied to the DR/GPS integrated positioning system of ships. The filters used for comparison include the EKF, DD1, HDD1, UKF, DD2, and HDD2. The results show that with a very large kernel bandwidth, the performances of the CDD1 and CDD2 filters are similar to those of standard DD1 and DD2 filters; with a proper kernel bandwidth, the CDD2 filter exhibits the best performance for the non-Gaussian noise case. Moreover, extensions of this research might include combining it with adaptive filtering methods, considering the problem of continuous systems, and applying it to other practical examples.