2D and 3D Angles-Only Target Tracking Based on Maximum Correntropy Kalman Filters

In this paper, angles-only target tracking (AoT) problem is investigated in the non Gaussian environment. Since the conventional minimum mean square error criterion based estimators tend to give poor accuracy in the presence of large outliers or impulsive noises in measurement, a maximum correntropy criterion (MCC) based framework is presented. Accordingly, three new estimation algorithms are developed for AoT problems based on the conventional sigma point filters, termed as MC-UKF-CK, MC-NSKF-GK and MC-NSKF-CK. Here MC-NSKF-GK represents the maximum correntropy new sigma point Kalman filter realized using Gaussian kernel and MC-NSKF-CK represents realization using Cauchy kernel. Similarly, based on the unscented Kalman filter, MC-UKF-CK has been developed. The performance of all these estimators is evaluated in terms of root-mean-square error (RMSE) in position and % track loss. The simulations were carried out for 2D as well as 3D AoT scenarios and it was inferred that, the developed algorithms performed with improved estimation accuracy than the conventional ones, in the presence of non Gaussian measurement noise.


Introduction
In state estimation, Kalman filter (KF) is a recursive solution used in various applications, such as information fusion, system control, integrated navigation, target tracking, and GPS solutions [1][2][3][4]. Kalman filter gives optimal estimates provided the dynamical system is linear and the noises assumed are Gaussian. However, it is extended to nonlinear systems through suitable approximation of the nonlinear functions. Using the Taylor series to linearize the nonlinear functions, the popular extended Kalman filter (EKF) [5] was derived. Also various sigma point filters have been proposed in the literature such as unscented Kalman filter (UKF) [6], cubature Kalman filter (CKF) [7], new sigma point Kalman filter (NSKF) [8], to obtain improved estimation accuracy than the EKF.
Since these filters are based on minimum mean square error criterion, their performance is likely to get deteriorated in the presence of non Gaussian noises such as heavy tailed and impulsive noises [9]. This makes state estimation a very challenging problem in the presence of nonlinear models and non Gaussian noise. Other possible solutions that can provide robust state estimates are Gaussian sum filter (GSF) [10,11], particle filter (PF) [12], Huber's KF (HKF, also known as M-estimation) [13], H ∞ filter [14] etc.
In order to improve the robustness of nonlinear state estimators in the presence of non Gaussian noise, a local similarity measure called correntropy [15,16], based filter called correntropy filter (C-Filter) was first proposed in [17]. Since it was developed by replacing the minimum mean square error (MMSE) criterion with maximum correntropy criterion (MCC), it proved beneficial for non Gaussian systems, but only for linear systems [18]. This algorithm made use of least squares and fixed point iteration, but failed to incorporate covariance estimation. In order to avoid this, a maximum correntropy Kalman filter (MCKF) involving fixed point iteration and covariance propagation was proposed in [19]. Similar issue was also addressed in [20], which used a cost function consisting of weighted least square (WLS) and Gaussian kernel function, and hence was named as maximum correntropy criterion-Kalman filter (MCC-KF).
To deal with nonlinear systems, extensions to the existing conventional algorithms based on MCC criterion were also developed and were named as maximum correntropy EKF (MC-EKF) [21], maximum correntropy UKF (MC-UKF) [22] and maximum correntropy sparse grid Gauss-Hermite quadrature filter (MC-SGHQF) [23]. But in the presence of large outliers in measurements, these filters incurred analytical problems in calculating inverse of matrices. Thus, new algorithms involving new cost function, WLS and statistical linearisation were proposed in [24], which were called as MC-UKF-constant and MC-UKFadaptive [25]. In developing the above mentioned estimators, Gaussian kernel played an important role in suppressing the non Gaussian measurement noise. In target tracking applications, we may receive measurements which have larger outliers. This could prove to be a challenging task in successful estimation of states using Gaussian kernel as it may be difficult to select a proper kernel bandwidth. Hence, Gaussian kernel may not always prove to be the best choice for a kernel function. To overcome this drawback, a Cauchy kernel function is constructed which gives reasonable estimation accuracy for a wide range of kernel bandwidth [26,27].
This paper deals with angles only target tracking problem in 2D and 3D. The literature presents with many variations of this tracking problem such as when the target is a curvilinear manoeuvring target [28,29]. However, as is common in passive sonar target tracking applications, the objective here is to estimate the states of a moving constant velocity target with the help of angles-only measurements, but corrupted with non Gaussian noise. The observer continuously monitors for the signals, that are generated due to the sound radiated by the target. The AoT can also be performed with other measurement sources like IRST sensor [4], radar [30] and also through video tracking [31]. Any irregularities in these signals received by the observer can be considered as glint noise. A mixture model of two zero-mean Gaussian for glint noise has been proposed in [32]. This consists of one Gaussian density with high probability and small variance while the other has small probability of occurrence and large variance. Alternatively, it is also modelled in [33] as a mixture of zero mean with small variance. In this work, the non Gaussian noise in angular measurements is modelled as a mixture of Gaussian densities plus shot noise.
The main contribution of this paper is the development of three new nonlinear filters for AoT problem, MC-UKF-CK, MC-NSKF-GK and MC-NSKF-CK, and their performance evaluation in the context of angles-only tracking. Accordingly, conventional filters UKF and NSKF have been reformulated based on maximum correntropy criterion. MC-UKF and MC-NSKF based on Gaussian kernel (MC-UKF-GK, MC-NSKF-GK) and Cauchy kernel (MC-UKF-CK, MC-NSKF-CK) functions have been derived. The performance evaluation of these estimators are conducted considering RMSE in position and track loss as the two performance metrics and a comparative discussion is presented. The simulation results highlight that the existing solutions behave poorly in comparison to the proposed algorithms.
The rest of the paper is organised as follows. Section 2 describes the problem formulation for AoT in 2D as well as 3D. Section 3 illustrates the correntropy, its properties for two random variables and MCC. In Section 4, the already existing Gaussian kernel based MC state estimation framework is revisited. In Section 5, the Cauchy kernel based MC state estimation framework for nonlinear systems is derived. Section 6 briefly discuss about the state estimators on which the developed MCC framework is incorporated. Section 7 describes the realization of non Gaussian noise, followed by simulation study in Section 8. Finally, the concluding remarks are given in Section 9.

Problem Formulation
The aim of the angles only tracking problem is to track the target trajectory using the noise corrupted angular measurements. The dynamics of the target is assumed to be a constant velocity motion. The observer motion is deterministic, implying that the position and velocity of the observer is known to us. The 2D and 3D target observer dynamics is illustrated below.

Process Model
The target and observer state vector with position and velocity as its states is given as The discrete time linear process model representing the target motion is given as Now, the relative state vector dynamics is where X k , the relative vector is defined as F is the state transition matrix and w k−1 is zero mean Gaussian process noise with Q as the covariance matrix. For problem formulation in the two dimensional space (let n = 2), F, Q matrices are defined as, The target observer dynamics in 2D for a moderately nonlinear scenario, is shown in Figure 1. Similarly for n = 3, the state and the associated matrices are where T is the sampling time interval and q x , q y , q z are the power spectral densities of the process noise along the X, Y, and Z axes respectively.

2D AoT problem:
The only available measurements are the bearing angles through which the states of the relative state vector can be estimated. The measurement model and the true angle measurements for the problem can be represented as where v k shall be modelled as the non Gaussian noise. 3D AoT problem: Figure 3 represents the target observer dynamics in Cartesian coordinate.
The range vector r is defined as Here, the nonlinear noise corrupted measurements are bearing (β) and elevation ( ) angles respectively, where β ∈ [0, 2π] and ∈ [− π 2 , π 2 ]. The measurement model involving the bearing and elevation angle is where, Here, v k is to be modelled as the non Gaussian noise.

Correntropy Measure
Correntropy is directly related to the probability of how similar two random variables are in the joint space controlled by the kernel bandwidth. The kernel bandwidth controls the window in which the similarity has to be assessed, and hence provides a way to eliminate the detrimental effect of outliers [16]. If X and Y are random variables, correntropy is defined as where k σ denotes a positive definite kernel function, p XY (.) denotes the joint PDF of X and Y and E is the expectation operator. Since the joint density is not accessible and if only a finite number of data points N are available, a sample estimator can be defined as Here G σ (·) is the Gaussian kernel, defined as which is bounded, positive and reaches its maximum only when X = Y, leading to the maximum correntropy criterion (MCC). By taking the Taylor series expansion of the Gaussian kernel, correntropy can also be expressed as a weighted sum of all even order moments of (x i − y i ), i.e., On the other hand, Cauchy kernel based non-linear state estimators can be developed using Cauchy kernel instead of Gaussian kernel function. It is defined as [34] Here δ is a positive scalar, representing the Cauchy kernel bandwidth. Similar to the Gaussian kernel, it can be shown that the Cauchy kernel also incorporates the higher order moments, given as A detailed derivation of the above equation is given in Appendix A.

Gaussian Kernel Based Maximum Correntropy Estimation Framework
Let us consider the process model described in Equations (1) and (5). To accommodate for the large outliers in measurements, the noise v k is considered non-Gaussian. Hence for MC based estimation framework [24], the Gaussian assumption of v k is relaxed.
In order to deal with the non Gaussian noises in the measurement update step, a statistical linearisation approach is employed. Consider that the nonlinear function h(·), operating on vector random variables X k is evaluated at N-points χ k , k = 1, · · · , N, with Then the prior and cross covariance P k|k−1 and P Xz are given as The nonlinear measurement function is represented in terms of measurement slope matrix H k , and a constant term c k as h(X k ) ≈ H k X k + c k . Here H k and c k are computed by minimizing the weighted least squares, arg min Then the solutions are H k = P −1 k|k−1 P Xz and c k = z k|k−1 − H k X k|k−1 . As the Thus, the linearised measurement equation is given as Accordingly, a cost function is formulated with the help of weighted least squares (WLS) to handle Gaussian process noise. To handle non-Gaussian measurement noise, statistical linearisation approach was used to define WLS function which in turn is used in MCC. Hence the cost function can be defined as , ℘ and are adjustable weights. In order to find the optimal estimate of X k , the cost function has to be minimized i.e., and the solution can be obtained as ∂J where In order to guarantee the convergence of the algorithm to a corresponding state estimator that follows a complete Gaussian assumption (when the kernel bandwidth σ becomes infinity), the values for weights in J are taken as ℘ = 1 and = −2σ 2 . Then Equation (9) becomes Since L G k is related to X k , Equation (11) represents a fixed point equation that can be solved using the fixed point iteration algorithm considering X k equal to X k|k−1 in Equation (10). But as mentioned in [19,22,24], for a satisfactory estimation performance, a single iteration is sufficient. Hence, adopting the same approach leads to the modification of Equation (11) as A more appropriate form for K G k , in terms of reduced computational complexity, can be derived using the matrix inversion lemma (detailed derivation is given in Appendix B) as Now, the corresponding posterior error covariance matrix is given as

Cauchy Kernel Based Maximum Correntropy Estimation Framework
In this section, we derive a maximum correntropy estimation framework using Cauchy kernel for potential improvement in estimation accuracy, in the presence of large multi dimensional non Gaussian noise. Hence the cost function becomes where ℘ C and C are adjustable weights, and with ℵ being the same as that mentioned in Section 4. To obtain the optimal estimate of X k , We set ℘ C = 1 and C = δ so as to guarantee the convergence of the estimator when kernel bandwidth δ tends to ∞. Rearranging, Here also, L C k is related to X k and hence Equation (13) is a fixed point equation that is to be solved using fixed point iteration algorithm, considering X k equal to X k|k−1 . Using the same justification that was adopted in Gaussian kernel case that only a single iteration is required, the expression for posterior mean is obtained as As per the proof given in Appendix B, Kalman gain can be modified as Then the posterior error covariance matrix shall be calculated as Theorem 1. As the kernel bandwidth δ → ∞, the Cauchy kernel based MC estimator reduces to the standard nonlinear state estimation algorithm.
Proof. As the time update is the same for the developed algorithms with respect to the standard nonlinear state estimators, the prior mean and covariance is unchanged. Hence the focus shall be on the posterior mean and covariance. This implies that the Kalman gain equation has to be revisited. When δ → ∞, Substituting the Equations (7) and (16) and H k in K C k , we have Since the expression of K C k is similar to the Kalman gain of standard nonlinear state estimator, posterior mean is also the same. Now, for the posterior covariance P k|k , consider Equation (15), Post multiplying Equation (14) by (R k + H k P k|k−1 L C k H k ) on both sides give Using Equations (16) and (18) Substituting H k , we get P k|k = P k|k−1 − K C k P Xz . For the given condition, K C k = P Xz P zz −1 , then P Xz = P zz K C k . Thus P k|k will become P k|k = P k|k−1 − K C k P zz K C k , which matches with the posterior error covariance of standard nonlinear estimator.

Remark 1.
For systems with non-Gaussian noise with large probability of abnormal noise, small value of δ is likely to provide more robustness. If the occurrence of abnormal noise is less, large value of δ could be considered.

Remark 2.
Cauchy kernel based nonlinear estimator with different δ performs with more estimation accuracy than Gaussian kernel based nonlinear estimator with different σ. Hence it is easier to select a value for δ that can provide accurate and robust estimates in the presence of abnormal noise.

Nonlinear State Estimators
This section deals with the nonlinear state estimators UKF and NSKF with a generalized algorithm for kernel based MC estimator.

Unscented Kalman Filter (UKF)
In the Bayesian framework, when the functions are nonlinear, the integrals encountered are intractable in nature and has to be evaluated using suitable numerical approximation methods. The UKF, through its unscented transformation, provides a way to numerically evaluate these integrals. Assuming that the integral to be approximated is I(X) = h(X)p X (X)dX, and X ∼ N (X, P), the unscented transformation defines a set of sigma points (X i ) and weights (W i ) such that [35] where n is the dimension of the state space and N = 2n + 1.
The sigma points and weights are defined aŝ with κ being the tuning parameter andX is the mean.

New Sigma Point Kalman Filter (NSKF)
From Equation (19), it can observed that in the unscented transformation, the maximum weight is assigned to the mean value. All the other sigma points are assigned equal weights, i.e., same probability of occurrence. In NSKF, a new approach was considered such that the sigma points closer to the mean will have more probability of occurrence. To realize this, a new method was formulated for defining the sigma points and weights, stated as [8] , i = n + 1, · · · , 2n Now the total number of sigma points N = 4n + 1, P i and S i denote the i th column of P and S respectively, and SS = P. The scalar variables are defined as b > { 1 4 max(mα i ) − 1 2 ∑ n i=1 α i }, m ∈ (0.5, 1) and α i = |<X,P i >| X 2 P i 2 . The Algorithm 1 for the developed estimators, both Gaussian kernel and Cauchy kernel based is given below. In this algorithm, K k and L k can be defined as per the chosen kernel function. R k is the noise covariance matrix which is assumed to be known in case there are no measurement outliers.
Posterior covariance: P k|k = (I − K k H k )P k|k−1 (I − K k H k ) + K k R k K k .

Modelling of Non Gaussian Noise in Angular Measurements
As mentioned in [36], a suitable way of modelling glint noise is to assume a Gaussian mixture. It is observed that the glint is more like Gaussian around the mean but has a non-Gaussian nature towards the tail region. The tail region represents the outliers, termed as glint spikes [32]. But shot noise, on the other hand, is modelled as an impulse with fixed amplitude at specific time steps. The mixture density of glint noise is modelled as The non Gaussian noises for angular measurements have been modelled by taking appropriate values for µ, σ 1 and σ 2 .

Simulation Results
The scenario for angles-only tracking problem in 2D as well as 3D Cartesian coordinate frame is considered in this section. The parameters required for generating the targetobserver dynamics, and simulation results are discussed. For simulations, moderately nonlinear tracking scenario for 2D as well as for 3D is considered. The simulation is carried for 1000 Monte Carlo runs with sampling time interval denoted as T. The entire tracking scenario is implemented and simulated in MATLAB software. Figure 4 shows the tracking performance of MC-NSKF-CK, where the estimated target path is plotted along with the truth target path, and the observer path for a single Monte Carlo run. It should be noted that for each run, observer path remains the same where as the target path varies due to the process noise. Further, the filter initialisation is also changing because of the randomness introduced in each run, as mentioned in Equation (21).

2D Scenario and Filter Initialisation
The filter is initialised as given in [10]. It is to be noted that for filter initialisation, we need an initial guess for speed, initial course and range of the target. Considering the problem at hand, these estimates have to be obtained from the initial angle measurement received. From these initial guess for parameters, the initial estimate for the states are obtained which are the positions and velocities. Accordingly, the initial range, target course and speed values are considered and mentioned in the Table 1. They are defined as s = N (s, σ 2 s ), c r = N (c r , σ 2 c ) and r = N (r, σ 2 r ) where c r can be defined as c r = z 0 + π with z 0 as the first bearing measurement. Finally the initial state vector X 0|0 and the initial covariance P 0|0 is calculated as where    Figure 5 shows the estimated target path obtained from MC-NSKF-CK and truth target path with observer trajectory. The initial parameter values required for generating the 3D scenario is given in the Table 2. Assuming that there are no outliers in the measurement, R k is defined as R k = diag(σ β , σ ). The bearing angle β is calculated with reference to the true North. For each Monte Carlo run, according to the new measurement received, initial range r and speed of the target s is assumed as mentioned in Table 2. According to these values, the relative state is initialised using the range estimate r ∼ N (r, σ r 2 ), initial bearing and elevation estimateβ 1 and 1 with headings α 1 = β 1 + π rad/s and γ 1 = 0 rad/s, and the initial speed estimate s ∼ N (s, σ s 2 ) with s as 0.258 km/s. The σ α 1 = π/ √ 12 and σ γ 1 = π/60 respectively [37]. The initial relative state vector X 0|0 is given as [38]

Performance Metrics
Performance analysis of the estimators formulated is evaluated by considering the below mentioned error statistics.
1. RMSE: Root-mean-square error in resultant target position is computed as follows where k denotes the time steps and M the total number of Monte Carlo runs.

Track Divergence:
In order to identify if a track is divergent or not, a certain threshold value (T b ) is set according to the position error computed at the final time instant of observation (k max ) as So, if the difference between estimated and truth target position is more than the threshold value (T b ), then we can say that the estimated path is moving away from the truth path. Thus, the track is considered to be divergent, and the number of such tracks are counted over M Monte Carlo runs.

Performance Analysis
The performance analysis of the developed filters is evaluated in the presence of glint plus shot noise in angle measurements. The accuracy of the estimators are evaluated by computing root mean square error (RMSE) in position at the end of the simulation period by imposing a track loss condition of 1 km.
The measurement noise v k for both the scenarios are given in Equations (22) and (23), respectively. Here, σ θ 1 = 0.5 • , σ θ 2 = 5 • , (σ β 1 , σ 1 ) as 0.0001 rad and (σ β 2 , σ 2 ) as 0.01 rad. The noise corrupted angle measurement for 2D is plotted as Figure 6. For 3D, the noise corrupted bearing and elevation angle are as shown in Figures 7 and 8 respectively. For illustration, in the figures, we have also plotted the angle measurements with Gaussian noise. With track loss condition of less than 1 km, the observed RMSE in position at the last time instant and percentage track loss for 2D as well as 3D is given in the Tables 3 and 4 respectively. Also, RMSE in resultant position (after excluding the diverged tracks) is evaluated and plotted in Figures 9 and 10. From these figures it can be inferred that in the presence of non Gaussian noise the estimation accuracy of UKF deteriorates, whereas filters based on MC framework performed with superior estimation accuracy. From the tabulation results it is evident that the Cauchy kernel based MC-UKF and MC-NSKF gives 108.9 m, 108.8 m RMSE and 1.1 and 0.5% track loss which is much less than that of the conventional UKF and NSKF which gives 152.8 m and 151.1 m RMSE with 4.4 and 2.8% track loss in 2D scenario. Similar observations can be made with respect to Gaussian kernel based MC-UKF and MC-NSKF giving much better accuracy but slightly less than Cauchy kernel MC framework. However, in 3D scenario, the MC based filters gave even superior estimation efficiency than that of the 2D scenario. UKF and NSKF in 3D with non Gaussian noise resulted in 100% track loss. Hence it can be inferred that for the given problem set up and noise statistics, UKF and NSKF failed to give estimates that met the track loss condition set, where as the Gaussian and Cauchy kernel based maximum correntropy filters gave more robust and accurate estimates. This can be inferred from the simulation results where the developed filters incurred only 13 to 14% track loss, with a final error in range of around 500 m. All these simulations are carried out by assuming the bandwidth (σ, δ) for 2D as (9,70) and for 3D as (11,75) such that the estimators can achieve maximum estimation accuracy. Also, the tuning parameter value of NSKF, m = 0.6 is assumed for simulation.

Conclusions
Since, measurements obtained in target tracking scenarios are corrupted with non Gaussian noise, this paper presents a maximum correntropy framework for 2D as well as 3D angles-only target tracking problem. The reformulation of UKF and NSKF in terms of Gaussian and Cauchy kernel based MC framework was realized. The non Gaussian noise is modelled as a Gaussian mixture (glint noise) plus shot noise. Finally, the performance of the estimators were evaluated and a comparative analysis is presented on the basis of RMSE in position and % track loss. From the comparative analysis, it can be concluded that the Gaussian and Cauchy kernel based MC framework provides improved estimation accuracy than UKF and NSKF in non Gaussian noise environments. Thus, it can be inferred that MC based estimators have the potential to give accurate and robust state estimates in the presence of non Gaussian noises in angle measurements. As a future work, the proposed estimation framework can be extended to track a manoeuvring target in the presence of angles-only measurements corrupted with non Gaussian noise.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Power Series Expansion of Cauchy Kernel Function
The binomial expansion of (1 + x) −N for negative integer −N is given as follows: Now the correntropy measure, by taking the binomial series expansion of Cauchy

Appendix B. Derivation of Kalman Gain
The Kalman gain for Gaussian kernel based nonlinear estimator is K G k with L G k as a scalar term. Similarly, for Cauchy kernel, it is K C k and L C k respectively. A general expression for Kalman gain is given as K k = (P −1 k|k−1 + L k H R −1 k H) −1 L k H R −1 k . Applying matrix inversion lemma After certain algebraic manipulations, we get K k = P k|k−1 L k H (I + R −1 k HP k|k−1 L k H ) −1 R −1 k = P k|k−1 L k H k (R k + H k P k|k−1 L k H k ) −1 .
From the above equation, K G k and K C k can be defined by making necessary substitution for L G k and L C k .