Correntropy Based Divided Difference Filtering for the Positioning of Ships

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.


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 1 and 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.

Correntropy
The correntropy is an important measure in ITL, which measures the similarity between two random variables X ∈ R and Y ∈ R. Given the joint distribution function of X and Y, F XY (x, y), the correntropy is defined by where E[·] represents the expectation operator, and κ(·, ·) is a shift-invariant Mercer kernel. In this paper, the Gaussian kernel is chosen as the kernel function: where e = x − y, and σ > 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: being N samples drawn from F XY . Taking the Taylor series expansion for the Gaussian kernel yields 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.

DD1 Filter
In this paper, the nonlinear state and measurement equations are described by where x(k) ∈ R n and y(k) ∈ R m are the n-dimensional state vector and m-dimensional measurement vector at time step k. f (·) and h (·) are the continuous and differentiable state functions and the measurement function. q k−1 and r k are the process and measurement noises, which are assumed i.i.d and independent of states, and with means denoted by q k−1 and r k and covariance matrices denoted by Q k−1 and R k . The square-root decompositions of the predicted state error covariance P k , update state error covariance P k , process noise covariance Q k , and measurement noise covariance R k are introduced as The DD1 filter uses the first-order divided differences to approximate. Let the element of i-th row, j-th column of S x x k be denoted as S x x k i,j , i.e., S x x k = S x x k i,j , and similarly for other matrices. The four matrices containing first-order divided differences are defined as where S x k−1 ,j is the j-th column of S 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 where H (·) denotes a Householder transformation of the augment matrix [17]. The predicted measurement and the Cholesky factor of corresponding covariance are given by The Kalman gain is computed as The updated state and the Cholesky factor of corresponding covariance are given by

DD2 Filter
The DD2 filter uses the second-order divided differences to approximate. Four matrices containing second-order divided differences are defined as The predicted state and the Cholesky factor of corresponding covariance are given by The predicted measurement and the Cholesky factor of corresponding covariance are given by The Kalman gain is computed as The updated state and the Cholesky factor of corresponding covariance are given by

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 y k = h (x k ) + 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 x k , the predicted state error is written as ξ k = x k − x k . Then, a nonlinear regression model is reconstructed as follows: where δ k = −ξ k r k , with the Cholesky factor of covariance of δ k is given by Left multiplying both sides of Equation (33) by D −1 k , we have the following regression model: where The problem associated with Equation (35) can be solved by making use of the MCC, and the corresponding cost function is given by where e k,i is the i-th component of the vector e k = z k − m (x k ), and L is the dimension of e k . Then, the solution of the problem mentioned previously can be found by setting the first derivation of Equation (36) to be zero: where φ (e k,i ) = G σ (e k,i ) · e k,i . By defining C (e k,i ) = φ (e k,i ) e k,i = G σ (e k,i ) and C = diag [C (e k,i )] = C x 0 0 C y , Equation (37) can be written in matrix form as In fact, the MCC uses C to re-weight the covariance matrix of the residual error e k and reconstruct the measurement information. Thus, the updated covariance can be written as and decomposed into two portion so that The initial value can be set as x (0) k = x k or be equal to the updated state computed from the corresponding standard DD1 or DD2. Then, we can obtain the following equations 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 S (1) Similarly, a one-step correntropy update for the DD2 filter can be written as 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 σ → ∞, the matrix C → 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 C y will be nearly singular. In view of this problem, a method is introduced as If |a k | > θ, with θ being a positive threshold, only the predicted step is worked, that is x k = x k , S x k = S x k . If |a k | θ, 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.

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.

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 where ϕ and λ 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 Ω represents the change rate of the ship's course. Correspondingly, the state equations of the ship are written as [40] where β denotes the inverse correlation time of ocean current, q i are independent Gaussian white noises. By discretizing these equations, we obtain the following equations where T is sampling period. Then, the measurement vector is chosen as where ϕ G and λ 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 where η 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: 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 Tables 1 and 2, and the process noises satisfy the Gaussian distributions:

Simulation Results
First, we assume the measurement noises satisfy the Gaussian distributions: TMSEs of position in Gaussian noises are revealed in Table 3. 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:   Table 4 summarizes the corresponding TMSEs. 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.

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.