# An Improved Gaussian Mixture CKF Algorithm under Non-Gaussian Observation Noise

Authors: Hongjian Wang; Cun Li
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

In order to solve the problems that the weight of Gaussian components of Gaussian mixture filter remains constant during the time update stage, an improved Gaussian Mixture Cubature Kalman Filter (IGMCKF) algorithm is designed by combining a Gaussian mixture density model with a CKF for target tracking. The algorithm adopts Gaussian mixture density function to approximately estimate the observation noise. The observation models based on Mini RadaScan for target tracking on offing are introduced, and the observation noise is modelled as glint noise. The Gaussian components are predicted and updated using CKF. A cost function is designed by integral square difference to update the weight of Gaussian components on the time update stage. Based on comparison experiments of constant angular velocity model and maneuver model with different algorithms, the proposed algorithm has the advantages of fast tracking response and high estimation precision, and the computation time should satisfy real-time target tracking requirements.


---

## Body

## 1. Introduction

With the universality of nonlinear problems, the principle and method of nonlinear filtering are being widely used for the nonlinear systems. One of the applications is target states estimation.For the currently nonlinear systems, the possible solution approaches are some approximate methods, such as approximating the probability density of system state as the Gaussian density, which is called Gaussian filtering. The Gaussian filtering could sum up three classes according to the different approximate methods: the first is function approximation, which approximates the nonlinear system function using the low-order expansion such as the extended Kalman filtering (EKF) [1] and the improved algorithms (adaptive fading EKF [2], strong tracking EKF [3], and central difference Kalman filtering (CDKF) [4, 5]). The EKF is a suboptimal filter that demands not only that the system has accurate states and accurate observation models but also that the observation noise is Gaussian in nature. However, the Gaussian Hypothesis is poor when the system model is highly nonlinear, and the estimation results are divergent. The CDKF based on Stirling’s polynomial interpolation has better performance including accuracy, efficiency, and stability than EKF for the nonlinear problems but would cause the greater computation burden. The second is deterministic sampling approximation method, which approximates system state and the probability density using deterministic sampling such as the Unscented Kalman Filter (UKF) [6, 7] and the improved algorithms [8, 9]. In principle, the UKF is simple and easy to implement, as it does not require the calculation of Jacobians at each time step. However, the UKF requires updating the matrix with sigma points with weighted negative values, which makes it difficult to preserve the positive definiteness during the iterative filter process. The third is approximation using quadrature, which approximates the multidimensional integrals of Bayesian recursive equation using some numerical technologies, such as the Gaussian-Hermite Kalman Filter (GHKF) [10] and the Cubature Kalman Filter (CKF) [11–13]. GHKF acquires the statistic characteristics after nonlinear transformation by Gaussian-Hermite numerical integral, which has higher accuracy. The shortcoming of GHKF is that it may not be suitable for addressing the high-dimensional nonlinear filtering issue. The CKF is proposed based on the spherical-radial cubature rule and can solve the Bayesian filter integral problem using cubature points with the same weight. The high-degree CKFs have high accuracy and stability but at a high computational cost.In actual application on target tracking, the noise of system processes or observations do not have ideal Gaussian density. In addition, the various Gaussian approximate filtering algorithms based on Gaussian noise do not display ideal performance because of the mismatched model. Besides Gaussian filter, other solution approaches for the nonlinear estimation problem involve the particle filter (PF) [14] based on the Monte Carlo numerical integral theory and sequence importance sampling and Gaussian sum filter [15, 16] based on the mixture of several Gaussian components. Although the PF do not require any assumption about the probability density function, they inevitably face enormous undesirable computational demands owing to the numerous stochastic sampling particles for fulfilling the estimation accuracy. The PF algorithm is unable to avoid the disadvantages of particle degradation and sample dilution; hence, the research on PF is focused on the revolution of particle degradation and sample dilution [17, 18]. A Gaussian sum CKF has been proposed for bearings-only tracking problems in the literature [19], and this CKF displays comparable performance to the particle filter. An improved Gaussian mixture filter algorithm has been proposed for highly nonlinear passive tracking systems in the literature [20], and the limited Gaussian mixture model has been used to approximate the posterior density of the state, process noise, and measurement noise. However, in all of these methods, the weights of the Gaussian components are kept constant while propagating the uncertainty through a nonlinear system and are updated only in the stage of measurement update. This assumption is valid if the system is marginally nonlinear or measurements are precise and available very frequently. The same is not true for the general nonlinear case. Terejanu et al. [21, 22] proposed a new Gaussian sum filter by adapting the weights of Gaussian mixture model in case of both time update and measurement update, which could obtain a better approximation of the posterior probability density function.In this paper, we focus on the target tracking problem on offing and design an improved Gaussian mixture CKF based on the GMCKF algorithm. Firstly, the formulation of target tracking is described, and the observation model with glint noise based on sensor is introduced. Secondly, the derivate steps of IGMCKF based integral square difference and Gaussian mixture density are designed. At last, the comparison experiments of constant angular velocity model and maneuver model with different algorithms are present, respectively.


## 2. Problem Formulation

Consider the target tracking problem on offing, giving the motion model and observation model:(1)
x
k
=
F
x
k
-
1
+
w
k
-
1
,
z
k
=
h
x
k
+
v
k
,where x
k
=
x
,
x
˙
,
y
,
y
˙
T
k denotes the target states, including the position and velocity of the target, and w
k
-
1 denotes the state process noise, which is assumed as Gaussian density with covariance Q. v
k denotes the observation noise which is assumed as glint noise because of the complicated environment on offing and the characteristics of sensors. The observation noise could be approximated as glint noise which is composed of Gaussian noise and Laplace noise. We modelled the observation noise as 2 Gaussian components with different covariance:(2)
p
v
=
1
-
ε
N
v
;
u
1
,
R
1
+
ε
N
v
;
u
2
,
R
2
,where N
v
;
u
,
R denotes Gaussian noise with mean u and covariance R. ε
∈
0,1 is the glint frequency factor.Consider that the position of sensor isx
0
,
y
0, and the sensor in this paper is Mini RadaScan; the observation information of the target including distance and bearing is as follows:(3)
h
x
k
=
x
k
-
x
o
2
+
y
k
-
y
o
2
arctan
⁡
y
k
-
y
o
x
k
-
x
o
.


## 3. Design of Improved Gaussian Mixture Cubature Kalman Filter Algorithm

Lemma 1.
The probability density function of then-dimension vector x can be approximated by the following equation:(4)
p
x
=
∑
i
=
1
I
ω
i
N
x
;
μ
i
,
P
i
.The approximate error can be arbitrarily small as soon as the number of components I is large enough. ω
i denotes the weighted value of the ith component; N
x
;
μ
i
,
P
i denotes the Gaussian distribution with mean μ
i and covariance P
i. Consider(5)
∑
i
=
1
I
ω
i
=
1
,
N
x
;
μ
i
,
P
i
=
1
2
π
P
i
exp
⁡
-
1
2
x
-
μ
i
T
P
i
-
1
x
-
μ
i
.


### 3.1. Time Update

Consider that the discrete nonlinear system with Gaussian mixture added noise and the prior and posterior density can be indicated by a Gaussian sum. The process noise and observation noise could be approximated as follows:(6)
p
w
k
≈
∑
j
=
1
J
ω
k
p
j
N
w
k
;
w
-
k
,
Q
k
j
,
p
v
k
≈
∑
l
=
1
L
ω
k
m
l
N
v
k
;
v
-
k
,
R
k
l
,where α
k
j denotes the weighted value of the jth component of the process noise and β
k
l is the weighted value of the lth component of the observation noise. Consider(7)
∑
j
=
1
J
ω
k
p
j
=
1
,
∑
l
=
1
L
ω
k
m
l
=
1
.Then, consider that the posterior density can be approximated with a Gaussian mixture model at k
-
1:(8)
p
x
k
-
1
∣
z
1
:
k
-
1
=
∑
i
=
1
I
ω
k
-
1
s
i
N
x
k
-
1
;
x
^
k
-
1
i
,
P
k
-
1
i
.We can obtain the prior density of state transition equation by the following equation:(9)
p
x
k
∣
x
k
-
1
=
∑
j
=
1
J
ω
k
-
1
p
j
N
x
k
;
f
x
k
-
1
+
w
-
k
j
,
Q
k
j
.According to the Bayesian formula, the predicted state density function can be approximated using a Gaussian mixture model:(10)
p
x
k
∣
z
1
:
k
-
1
=
∑
j
=
1
J
∑
i
=
1
I
ω
k
-
1
s
i
ω
k
p
j
∫
R
n
x
N
x
k
-
1
;
x
^
k
-
1
i
,
P
k
-
1
i
·
N
x
k
;
f
x
k
-
1
+
w
-
j
,
Q
k
j
d
x
k
-
1
≈
∑
j
=
1
J
∑
i
=
1
I
ω
k
∣
k
-
1
i
,
j
N
x
k
;
x
^
k
∣
k
-
1
i
,
j
,
P
k
∣
k
-
1
i
,
j
≈
∑
r
=
1
I
·
J
ω
k
∣
k
-
1
t
r
N
x
k
;
x
^
k
∣
k
-
1
r
,
P
k
∣
k
-
1
r
,where ω
k
∣
k
-
1
t
r
=
ω
k
-
1
s
i
ω
k
p
j. R is the set of real numbers, and n
x is the dimension of the state vector. I and J are the number of Gaussian components of system state and process noise, respectively. x
^
k
∣
k
-
1
r and P
k
∣
k
-
1
r could be calculated by the time update steps of CKF:(11)
x
^
k
∣
k
-
1
r
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
i
+
w
-
k
j
,
P
k
∣
k
-
1
r
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
i
ξ
c
,
k
∣
k
-
1
T
i
-
x
^
k
∣
k
-
1
r
-
w
-
k
j
x
^
k
∣
k
-
1
r
-
w
-
k
j
T
+
Q
k
j
,where m is the number of cubature points.


### 3.2. Adaptive Weight Update

Consider the following nonlinear system (1) with the probability density function of the initial conditions p
x
0. According to formula (4), the Gaussian mixture approximation of the probability density function is given:(12)
p
^
x
k
=
∑
i
=
1
I
ω
k
s
i
N
x
k
i
;
x
^
k
i
,
P
k
i
.The true probability density function of system state by the Chapman-Kolmogorov equation is given:(13)
p
x
k
=
∫
R
n
p
x
k
∣
x
k
-
1
p
x
k
-
1
d
x
k
-
1
.Mean-square optimal new weights can be obtained by minimizing the following integral square difference between the true probabilityp
x
k and the approximation one p
^
x
k in the least square algorithm:(14)
w
k
∣
k
-
1
=
arg
⁡
min
⁡
1
2
∫
R
n
p
x
k
-
p
^
x
k
2
d
x
k
,where w
k
∣
k
-
1
=
ω
k
∣
k
-
1
1
ω
k
∣
k
-
1
2
⋯
ω
k
∣
k
-
1
I
T denotes the vector of the weights of every Gaussian component at time k. In order to resolve formula (17), the cost function is given as follows [23]:(15)
J
w
k
∣
k
-
1
=
J
p
p
w
k
∣
k
-
1
-
2
J
p
p
^
w
k
∣
k
-
1
+
J
p
^
p
^
w
k
∣
k
-
1
,where the first term represents the self-likeness of the true probability density function of the system state and J
p
p
w
k
∣
k
-
1
=
∫
R
n
p
x
k
p
x
k
d
x
k. The first term is not needed in the optimization process and it is used only to provide an overall magnitude of the uncertainty propagation error [24]. The second represents the cross-likeness of the true probability and the approximation one, and J
p
p
^
w
k
∣
k
-
1
=
∫
R
n
p
x
k
p
^
x
k
d
x
k. The last term is the self-likeness of the approximation probability of the system state, and J
p
^
p
^
w
k
∣
k
-
1
=
∫
R
n
p
^
x
k
p
^
x
k
d
x
k.Formula (16) is based on the assumption that the Gaussian mixture approximation is equal to the true probability density function at time k
-
1; namely, p
^
x
k
-
1
=
p
x
k
-
1. Consider(16)
p
x
k
=
∫
R
n
p
x
k
∣
x
k
-
1
p
^
x
k
-
1
d
x
k
-
1
.Now the derivation of the terms of cost function is given as follows:(17)
J
p
p
w
k
∣
k
-
1
=
∫
R
n
p
x
k
2
d
x
k
=
∫
R
n
∫
R
n
p
x
k
∣
x
k
-
1
p
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
p
x
k
∣
x
k
-
1
p
^
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
p
^
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
=
∑
i
=
1
I
ω
k
-
1
s
i
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
d
x
k
-
1
2
d
x
k
=
w
k
-
1
T
J
p
p
w
k
-
1
.Similarly,(18)
J
p
p
^
w
k
∣
k
-
1
=
w
k
∣
k
-
1
T
J
p
p
^
w
k
-
1
,
J
p
^
p
^
w
k
∣
k
-
1
=
w
k
∣
k
-
1
T
J
p
^
p
^
w
k
∣
k
-
1
.The elements of the matrix are given as follows:(19)
J
i
j
p
p
=
∫
∫
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
d
x
k
-
1
2
d
x
k
=
∫
E
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
2
d
x
k
,
J
i
j
p
p
^
=
∫
N
f
x
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
+
Q
k
-
1
N
x
k
-
1
;
x
k
-
1
j
,
P
k
-
1
j
d
x
k
-
1
=
E
N
x
k
-
1
;
x
k
-
1
j
,
P
k
-
1
j
N
f
x
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
+
Q
k
-
1
=
∑
c
=
1
m
ω
c
N
ξ
c
,
k
∣
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
,where ω
c
,
ξ
c is the cubature points. Consider(20)
J
i
j
p
^
p
^
=
N
x
k
∣
k
-
1
i
;
x
k
∣
k
-
1
j
,
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
j
=
2
π
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
i
-
1
/
2
exp
⁡
-
1
2
x
k
∣
k
-
1
i
-
x
k
∣
k
-
1
j
T
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
i
-
1
x
k
∣
k
-
1
i
-
x
k
∣
k
-
1
j
.According to the above relations, the final formulation of formula (14) could be obtained as follows:(21)
w
k
∣
k
-
1
=
arg
⁡
min
⁡
1
2
w
k
∣
k
-
1
T
J
p
^
p
^
w
k
∣
k
-
1
-
w
k
∣
k
-
1
T
J
p
p
^
w
k
-
1
;the result of (21) is the weights of components after update.


### 3.3. Measurement Update

Consider the observation equation; the likelihood density of observation equation can be approximated as follows:(22)
p
z
k
∣
x
k
=
∑
l
=
1
L
ω
k
m
l
N
z
k
;
h
x
k
+
v
-
k
l
,
R
k
l
.Then, we can obtain the following equation:(23)
p
x
k
∣
z
1
:
k
=
∑
l
=
1
L
∑
r
=
1
I
·
J
ω
k
m
l
ω
k
∣
k
-
1
t
r
∫
R
n
x
N
x
k
;
x
^
k
∣
k
-
1
r
,
P
k
∣
k
-
1
r
N
z
k
;
h
x
k
+
v
-
k
l
,
R
k
l
d
x
k
=
∑
p
=
1
L
∑
r
=
1
I
·
J
ω
k
m
l
ω
k
∣
k
-
1
t
r
N
x
k
;
x
^
k
r
,
l
,
P
k
r
,
l
=
∑
n
=
1
I
·
J
·
L
ω
k
n
N
x
k
;
x
^
k
n
,
P
k
n
,where (24)
ω
k
n
=
ω
k
∣
k
-
1
t
r
ω
k
m
l
p
z
k
∣
x
k
,
n
∑
r
=
1
I
·
J
∑
l
=
1
L
ω
k
∣
k
-
1
t
r
ω
k
m
l
p
z
k
∣
x
k
,
n
.The term ω
k
-
1
s
i in ω
k
∣
k
-
1
t
r
=
ω
k
-
1
s
i
ω
k
p
j is calculated by the adaptive weight update (21). Consider(25)
p
z
k
∣
x
k
,
n
=
1
2
π
σ
n
2
exp
⁡
-
1
2
z
k
-
z
^
k
∣
k
-
1
r
,
l
σ
n
2
.The estimation of state of Gaussian components is as follows:(26)
x
^
k
n
=
x
^
k
∣
k
-
1
r
+
K
k
r
,
l
z
^
k
-
z
^
k
∣
k
-
1
r
,
l
,
P
k
n
=
P
k
∣
k
-
1
r
-
K
k
r
,
l
P
z
z
,
k
∣
k
-
1
r
,
l
K
k
T
r
,
l
.The Kalman gain K is as follows:(27)
K
k
r
,
l
=
P
x
z
,
k
∣
k
-
1
r
,
l
P
z
z
,
k
∣
k
-
1
-
1
r
,
l
.The covariance P
z
z and cross covariance P
x
z are, respectively, as follows:(28)
P
z
z
,
k
∣
k
-
1
r
,
l
=
∑
c
=
1
m
ω
c
ε
c
,
k
∣
k
-
1
r
ε
c
,
k
∣
k
-
1
T
r
-
z
^
k
∣
k
-
1
r
-
v
-
k
l
z
^
k
∣
k
-
1
r
-
v
-
k
l
T
+
R
k
l
,
(29)
P
x
z
,
k
∣
k
-
1
r
,
l
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
r
ε
c
,
k
∣
k
-
1
T
r
-
x
^
k
∣
k
-
1
r
z
^
k
∣
k
-
1
r
-
v
-
k
l
T
.Hence, the output of the filter is as follows:(30)
x
^
k
=
∑
n
=
1
I
·
J
·
L
ω
k
n
x
^
k
n
,
(31)
P
k
=
∑
n
=
1
I
·
J
·
L
ω
-
k
n
P
k
n
+
x
^
k
n
-
x
^
k
x
^
k
n
-
x
^
k
T
,where n
=
r
-
1
L
+
l.


### 3.4. Merging the Gaussian Components

As we know from (31), the number of Gaussian components is I
·
J
·
L. If I
·
J
·
L
>
I, some of the Gaussian components must be merged after the measurement update stage so that the total number of Gaussian components calculated in the next time index can be reduced to I.There is no harm in considering the index of Gaussian components to bet
=
1
,
…
,
I
·
J
·
L, according to their weight values in descending order, such that the Gaussian component with the largest weight value has an index t
=
1. The Gaussian components with indices t
=
1
,
…
,
I
-
1 are selected first. The other Gaussian components with indices t
=
I
,
…
,
I
·
J
·
L will be merged. The weight value of the Ith component is as follows:(32)
ω
k
I
=
∑
t
=
I
J
·
L
ω
k
t
,and the mean and covariance are, respectively, as follows:(33)
x
^
k
I
=
∑
t
=
I
J
·
L
ω
~
k
t
x
^
k
t
,
(34)
P
k
I
=
∑
t
=
I
J
·
L
ω
~
k
t
·
P
k
t
+
x
^
k
t
-
x
^
k
I
x
^
k
t
-
x
^
k
I
T
,where ω
~
k
t
=
ω
k
t
/
ω
k
I are the normalised weight values of the Gaussian components to be merged.


## 3.1. Time Update

Consider that the discrete nonlinear system with Gaussian mixture added noise and the prior and posterior density can be indicated by a Gaussian sum. The process noise and observation noise could be approximated as follows:(6)
p
w
k
≈
∑
j
=
1
J
ω
k
p
j
N
w
k
;
w
-
k
,
Q
k
j
,
p
v
k
≈
∑
l
=
1
L
ω
k
m
l
N
v
k
;
v
-
k
,
R
k
l
,where α
k
j denotes the weighted value of the jth component of the process noise and β
k
l is the weighted value of the lth component of the observation noise. Consider(7)
∑
j
=
1
J
ω
k
p
j
=
1
,
∑
l
=
1
L
ω
k
m
l
=
1
.Then, consider that the posterior density can be approximated with a Gaussian mixture model at k
-
1:(8)
p
x
k
-
1
∣
z
1
:
k
-
1
=
∑
i
=
1
I
ω
k
-
1
s
i
N
x
k
-
1
;
x
^
k
-
1
i
,
P
k
-
1
i
.We can obtain the prior density of state transition equation by the following equation:(9)
p
x
k
∣
x
k
-
1
=
∑
j
=
1
J
ω
k
-
1
p
j
N
x
k
;
f
x
k
-
1
+
w
-
k
j
,
Q
k
j
.According to the Bayesian formula, the predicted state density function can be approximated using a Gaussian mixture model:(10)
p
x
k
∣
z
1
:
k
-
1
=
∑
j
=
1
J
∑
i
=
1
I
ω
k
-
1
s
i
ω
k
p
j
∫
R
n
x
N
x
k
-
1
;
x
^
k
-
1
i
,
P
k
-
1
i
·
N
x
k
;
f
x
k
-
1
+
w
-
j
,
Q
k
j
d
x
k
-
1
≈
∑
j
=
1
J
∑
i
=
1
I
ω
k
∣
k
-
1
i
,
j
N
x
k
;
x
^
k
∣
k
-
1
i
,
j
,
P
k
∣
k
-
1
i
,
j
≈
∑
r
=
1
I
·
J
ω
k
∣
k
-
1
t
r
N
x
k
;
x
^
k
∣
k
-
1
r
,
P
k
∣
k
-
1
r
,where ω
k
∣
k
-
1
t
r
=
ω
k
-
1
s
i
ω
k
p
j. R is the set of real numbers, and n
x is the dimension of the state vector. I and J are the number of Gaussian components of system state and process noise, respectively. x
^
k
∣
k
-
1
r and P
k
∣
k
-
1
r could be calculated by the time update steps of CKF:(11)
x
^
k
∣
k
-
1
r
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
i
+
w
-
k
j
,
P
k
∣
k
-
1
r
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
i
ξ
c
,
k
∣
k
-
1
T
i
-
x
^
k
∣
k
-
1
r
-
w
-
k
j
x
^
k
∣
k
-
1
r
-
w
-
k
j
T
+
Q
k
j
,where m is the number of cubature points.


## 3.2. Adaptive Weight Update

Consider the following nonlinear system (1) with the probability density function of the initial conditions p
x
0. According to formula (4), the Gaussian mixture approximation of the probability density function is given:(12)
p
^
x
k
=
∑
i
=
1
I
ω
k
s
i
N
x
k
i
;
x
^
k
i
,
P
k
i
.The true probability density function of system state by the Chapman-Kolmogorov equation is given:(13)
p
x
k
=
∫
R
n
p
x
k
∣
x
k
-
1
p
x
k
-
1
d
x
k
-
1
.Mean-square optimal new weights can be obtained by minimizing the following integral square difference between the true probabilityp
x
k and the approximation one p
^
x
k in the least square algorithm:(14)
w
k
∣
k
-
1
=
arg
⁡
min
⁡
1
2
∫
R
n
p
x
k
-
p
^
x
k
2
d
x
k
,where w
k
∣
k
-
1
=
ω
k
∣
k
-
1
1
ω
k
∣
k
-
1
2
⋯
ω
k
∣
k
-
1
I
T denotes the vector of the weights of every Gaussian component at time k. In order to resolve formula (17), the cost function is given as follows [23]:(15)
J
w
k
∣
k
-
1
=
J
p
p
w
k
∣
k
-
1
-
2
J
p
p
^
w
k
∣
k
-
1
+
J
p
^
p
^
w
k
∣
k
-
1
,where the first term represents the self-likeness of the true probability density function of the system state and J
p
p
w
k
∣
k
-
1
=
∫
R
n
p
x
k
p
x
k
d
x
k. The first term is not needed in the optimization process and it is used only to provide an overall magnitude of the uncertainty propagation error [24]. The second represents the cross-likeness of the true probability and the approximation one, and J
p
p
^
w
k
∣
k
-
1
=
∫
R
n
p
x
k
p
^
x
k
d
x
k. The last term is the self-likeness of the approximation probability of the system state, and J
p
^
p
^
w
k
∣
k
-
1
=
∫
R
n
p
^
x
k
p
^
x
k
d
x
k.Formula (16) is based on the assumption that the Gaussian mixture approximation is equal to the true probability density function at time k
-
1; namely, p
^
x
k
-
1
=
p
x
k
-
1. Consider(16)
p
x
k
=
∫
R
n
p
x
k
∣
x
k
-
1
p
^
x
k
-
1
d
x
k
-
1
.Now the derivation of the terms of cost function is given as follows:(17)
J
p
p
w
k
∣
k
-
1
=
∫
R
n
p
x
k
2
d
x
k
=
∫
R
n
∫
R
n
p
x
k
∣
x
k
-
1
p
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
p
x
k
∣
x
k
-
1
p
^
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
p
^
x
k
-
1
d
x
k
-
1
2
d
x
k
=
∫
R
n
∫
R
n
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
=
∑
i
=
1
I
ω
k
-
1
s
i
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
d
x
k
-
1
2
d
x
k
=
w
k
-
1
T
J
p
p
w
k
-
1
.Similarly,(18)
J
p
p
^
w
k
∣
k
-
1
=
w
k
∣
k
-
1
T
J
p
p
^
w
k
-
1
,
J
p
^
p
^
w
k
∣
k
-
1
=
w
k
∣
k
-
1
T
J
p
^
p
^
w
k
∣
k
-
1
.The elements of the matrix are given as follows:(19)
J
i
j
p
p
=
∫
∫
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
d
x
k
-
1
2
d
x
k
=
∫
E
N
x
k
-
1
i
;
x
^
k
-
1
i
,
P
k
-
1
i
N
x
k
;
f
x
k
-
1
+
w
-
k
,
Q
k
2
d
x
k
,
J
i
j
p
p
^
=
∫
N
f
x
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
+
Q
k
-
1
N
x
k
-
1
;
x
k
-
1
j
,
P
k
-
1
j
d
x
k
-
1
=
E
N
x
k
-
1
;
x
k
-
1
j
,
P
k
-
1
j
N
f
x
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
+
Q
k
-
1
=
∑
c
=
1
m
ω
c
N
ξ
c
,
k
∣
k
-
1
;
x
k
∣
k
-
1
i
,
P
k
∣
k
-
1
i
,where ω
c
,
ξ
c is the cubature points. Consider(20)
J
i
j
p
^
p
^
=
N
x
k
∣
k
-
1
i
;
x
k
∣
k
-
1
j
,
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
j
=
2
π
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
i
-
1
/
2
exp
⁡
-
1
2
x
k
∣
k
-
1
i
-
x
k
∣
k
-
1
j
T
P
k
∣
k
-
1
i
+
P
k
∣
k
-
1
i
-
1
x
k
∣
k
-
1
i
-
x
k
∣
k
-
1
j
.According to the above relations, the final formulation of formula (14) could be obtained as follows:(21)
w
k
∣
k
-
1
=
arg
⁡
min
⁡
1
2
w
k
∣
k
-
1
T
J
p
^
p
^
w
k
∣
k
-
1
-
w
k
∣
k
-
1
T
J
p
p
^
w
k
-
1
;the result of (21) is the weights of components after update.


## 3.3. Measurement Update

Consider the observation equation; the likelihood density of observation equation can be approximated as follows:(22)
p
z
k
∣
x
k
=
∑
l
=
1
L
ω
k
m
l
N
z
k
;
h
x
k
+
v
-
k
l
,
R
k
l
.Then, we can obtain the following equation:(23)
p
x
k
∣
z
1
:
k
=
∑
l
=
1
L
∑
r
=
1
I
·
J
ω
k
m
l
ω
k
∣
k
-
1
t
r
∫
R
n
x
N
x
k
;
x
^
k
∣
k
-
1
r
,
P
k
∣
k
-
1
r
N
z
k
;
h
x
k
+
v
-
k
l
,
R
k
l
d
x
k
=
∑
p
=
1
L
∑
r
=
1
I
·
J
ω
k
m
l
ω
k
∣
k
-
1
t
r
N
x
k
;
x
^
k
r
,
l
,
P
k
r
,
l
=
∑
n
=
1
I
·
J
·
L
ω
k
n
N
x
k
;
x
^
k
n
,
P
k
n
,where (24)
ω
k
n
=
ω
k
∣
k
-
1
t
r
ω
k
m
l
p
z
k
∣
x
k
,
n
∑
r
=
1
I
·
J
∑
l
=
1
L
ω
k
∣
k
-
1
t
r
ω
k
m
l
p
z
k
∣
x
k
,
n
.The term ω
k
-
1
s
i in ω
k
∣
k
-
1
t
r
=
ω
k
-
1
s
i
ω
k
p
j is calculated by the adaptive weight update (21). Consider(25)
p
z
k
∣
x
k
,
n
=
1
2
π
σ
n
2
exp
⁡
-
1
2
z
k
-
z
^
k
∣
k
-
1
r
,
l
σ
n
2
.The estimation of state of Gaussian components is as follows:(26)
x
^
k
n
=
x
^
k
∣
k
-
1
r
+
K
k
r
,
l
z
^
k
-
z
^
k
∣
k
-
1
r
,
l
,
P
k
n
=
P
k
∣
k
-
1
r
-
K
k
r
,
l
P
z
z
,
k
∣
k
-
1
r
,
l
K
k
T
r
,
l
.The Kalman gain K is as follows:(27)
K
k
r
,
l
=
P
x
z
,
k
∣
k
-
1
r
,
l
P
z
z
,
k
∣
k
-
1
-
1
r
,
l
.The covariance P
z
z and cross covariance P
x
z are, respectively, as follows:(28)
P
z
z
,
k
∣
k
-
1
r
,
l
=
∑
c
=
1
m
ω
c
ε
c
,
k
∣
k
-
1
r
ε
c
,
k
∣
k
-
1
T
r
-
z
^
k
∣
k
-
1
r
-
v
-
k
l
z
^
k
∣
k
-
1
r
-
v
-
k
l
T
+
R
k
l
,
(29)
P
x
z
,
k
∣
k
-
1
r
,
l
=
∑
c
=
1
m
ω
c
ξ
c
,
k
∣
k
-
1
r
ε
c
,
k
∣
k
-
1
T
r
-
x
^
k
∣
k
-
1
r
z
^
k
∣
k
-
1
r
-
v
-
k
l
T
.Hence, the output of the filter is as follows:(30)
x
^
k
=
∑
n
=
1
I
·
J
·
L
ω
k
n
x
^
k
n
,
(31)
P
k
=
∑
n
=
1
I
·
J
·
L
ω
-
k
n
P
k
n
+
x
^
k
n
-
x
^
k
x
^
k
n
-
x
^
k
T
,where n
=
r
-
1
L
+
l.


## 3.4. Merging the Gaussian Components

As we know from (31), the number of Gaussian components is I
·
J
·
L. If I
·
J
·
L
>
I, some of the Gaussian components must be merged after the measurement update stage so that the total number of Gaussian components calculated in the next time index can be reduced to I.There is no harm in considering the index of Gaussian components to bet
=
1
,
…
,
I
·
J
·
L, according to their weight values in descending order, such that the Gaussian component with the largest weight value has an index t
=
1. The Gaussian components with indices t
=
1
,
…
,
I
-
1 are selected first. The other Gaussian components with indices t
=
I
,
…
,
I
·
J
·
L will be merged. The weight value of the Ith component is as follows:(32)
ω
k
I
=
∑
t
=
I
J
·
L
ω
k
t
,and the mean and covariance are, respectively, as follows:(33)
x
^
k
I
=
∑
t
=
I
J
·
L
ω
~
k
t
x
^
k
t
,
(34)
P
k
I
=
∑
t
=
I
J
·
L
ω
~
k
t
·
P
k
t
+
x
^
k
t
-
x
^
k
I
x
^
k
t
-
x
^
k
I
T
,where ω
~
k
t
=
ω
k
t
/
ω
k
I are the normalised weight values of the Gaussian components to be merged.


## 4. Simulation

### 4.1. Simulation 1

#### 4.1.1. Simulation Parameter

The proposed algorithm is simulated by Monte Carlo withMC
=
50, the simulation time is Time
=
300
s, simulation step size is T
=
1
s, the covariance of observation noise in (2) is u
1
=
u
2
=
0, R
1
=
R, R
2
=
100
R
1, and the glint frequency factor ε
=
0.1; the other filter parameters are set as follows:(35)
F
=
1
sin
⁡
ω
T
ω
0
cos
⁡
ω
T
-
1
ω
0
cos
⁡
ω
T
0
-
sin
⁡
ω
T
0
1
-
cos
⁡
ω
T
ω
1
sin
⁡
ω
T
ω
0
sin
⁡
ω
T
0
cos
⁡
ω
T
,
P
0
=
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
,
Q
=
0.05
0
0.1
0
0
0.05
0
0.1
,
R
=
5
0
0
5
.The other parameters are in Table 1.Table 1
Target and observation station state.
Target
Observation
Initialize position (m)
1
00
,
100
50
,
150
Initialize velocity (m/s)
2.8
0.3
Initialize heading (°)
45
160
Motion mode
CAV
CV


#### 4.1.2. Simulation Results and Analysis

We compare the performances of GMEKF, GMUKF, GMCKF, and IGMCKF algorithms with the conditions mentioned above. Figure1 shows the results of the trajectory estimations of the target from the different algorithms versus the truth trajectory. Figure 2 shows the position error of different algorithms, and Figure 3 compares the position errors in X direction and Y direction among the different algorithms.Figure 1
Comparison of the trajectory estimations.Figure 2
The curve of position error.Figure 3
The curve of position error inX and Y.In Figure1, the red dotted-dashed line denotes the trajectory of observation station, and the arrow denotes the movement heading.From Figure2, it can be observed that the GMEKF and GMUKF have bad stability, and the error can be up to 25 m. The performance improvement of GMCKF is much better, and the error could be up to 13 m. The IGMCKF has the best tracking performance with an error up to 10 m, and the error could be limited to 3 m during stabilization.Figure4 shows the tracking performance of target velocity of different algorithms. Figure 5 shows the tracking performance of velocity components in X direction and Y direction, respectively.Figure 4
The tracking curve of target velocity.Figure 5
The tracking curve of velocity components.From Figure4, it can be observed that the GMCKF algorithm could track the target velocity after about 30 s, while the other algorithms need about 80 s.Figure6 shows the comparison of velocity errors of different algorithms.Figure 6
The curve of velocity error.Figure7 shows the comparison of target headings of different algorithms.Figure 7
The estimation curve of target heading.The accuracy achieved by the filtering methods was analyzed using the following metrics.The Root-Mean-Square Error (RMSE) is defined as(36)
RMSE
=
1
t
step
∑
k
=
1
t
step
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
,where t
step is the simulation steps and x
k
,
y
k and x
^
k
,
y
^
k are the true position and estimated position, respectively, of target at time k.The average Noncredibility Index (NCI) [25, 26] is defined as(37)
NCI
=
1
t
step
∑
k
=
1
t
step
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
P
p
,
k
-
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
RMSE
,where P
p
,
k is the uncertainty information of the position estimation and P
p
,
k
=
P
x
,
k
2
+
P
y
,
k
2.Table2 gives the NCI, RMSE, and the computing time cost of different tested algorithms.Table 2
RMSE and time of different algorithms.
NCI
RMSE (m)
Time (ms)
GMEKF
1.1534
6.3300
36.105
GMUKF
0.7591
5.4350
130.631
GMCKF
0.6522
3.9051
163.330
IGMCKF
0.6280
2.0042
257.302From Table2, we can see that IGMCKF has the smallest RMSE which illuminates that the algorithm has the highest precision. The running time of IGMCKF is the longest because IGMCKF needs to update the weight during the time update stage and split the observation noise into several Gaussian components and merge them based on CKF estimating the cubature points. But the running time could satisfy the requirement of real-time target tracking. The NCI of the algorithms except the GMEKF are below 1, which illuminate that the credibility of the various filter algorithms is high.


### 4.2. Simulation 2

#### 4.2.1. Simulation Parameter

We will compare the performance of proposed algorithm (interacting multiple model CKF, IMMCKF) and interacting multiple model Gaussian mixture CKF (IMMGMCKF) algorithms. The filter parameters are set as follows. The initial position of target is20,150, the initial velocity is 0
,
-
1.5
m
/
s, and the acceleration is 0
m
/
s
2. The target moves with constant initial velocity in a straight line from 0 s to 100 s. Then it maneuvers and moves with constant acceleration 0.0075,0.0075
m
/
s
2 from 101 s to 200 s. From 201 s to 300 s, it moves with constant velocity 0.7575
,
-
0.7425
m
/
s. From 301 s to 380 s, it maneuvers with constant acceleration -
0.025
,
-
0.015
m
/
s
2, and it moves with constant velocity -
1.2175
,
-
1.9275
m
/
s for the last 20 s.The observation noise is set the same as simulation 1. Singer model is employed for the IGMCKF algorithm, the state vector isx
=
x
,
x
˙
,
x
¨
,
y
,
y
˙
,
y
¨, and other filter parameters are set as follows:(38)
F
=
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
0
0
0
0
0
0
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
,
P
0
=
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
0
0
0
0
0
0
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
,
Q
=
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
0
0
0
0
0
0
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
,
R
=
10
0
0
10
.For the IMMCKF and IMMGMCKF algorithms, three models are employed: two constant acceleration models with different acceleration and constant velocity model. The model initial distribution isu
=
0.8,0.1,0.1, and Markov transition matrix is as follows:(39)
Π
=
0.95
0.025
0.025
0.025
0.95
0.025
0.025
0.025
0.95
.


#### 4.2.2. Simulation Results and Analysis

The simulation results are obtained from 50 Monte Carlo runs. Figure8 shows the position trajectory of various algorithms. Figures 9 and 10 denote the RMSE of position tracking in X direction and Y direction, respectively. Figures 11 and 12 denote the RMSE of velocity tracking in X direction and Y direction, respectively.Figure 8
Comparison of the trajectory estimations.Figure 9
X position RMSE.Figure 10
Y position RMSE.Figure 11
X velocity RMSE.Figure 12
Y velocity RMSE.From the above simulation results and Table3, it can be seen that the IMMCKF has the worst accuracy of estimation, and because of that the observation noise is modelled as glint noise, and the IMMCKF algorithm is suitable for the Gaussian noise. The IMMGMCKF and IGMCKF have the similar accuracy of estimation and because of that the Gaussian mixture algorithm could deal with the non-Gaussian noise well. The running time of the IMMGMCKF is the longest and because of that the algorithm needs to switch between different motion models. The proposed algorithm has the best estimation accuracy and best computation complexity and could satisfy the real-time target tracking requirement.Table 3
Comparisons of average RMSE and running time of different algorithms.
X position (m)
Y position (m)
X velocity (m/s)
Y velocity (m/s)
Running time (ms)
IMMCKF
0.6713
0.6909
0.4741
0.4930
222
IMMGMCKF
0.1074
0.1238
0.0909
0.0920
315
IGMCKF
0.0983
0.0895
0.0463
0.0462
182


## 4.1. Simulation 1

### 4.1.1. Simulation Parameter

The proposed algorithm is simulated by Monte Carlo withMC
=
50, the simulation time is Time
=
300
s, simulation step size is T
=
1
s, the covariance of observation noise in (2) is u
1
=
u
2
=
0, R
1
=
R, R
2
=
100
R
1, and the glint frequency factor ε
=
0.1; the other filter parameters are set as follows:(35)
F
=
1
sin
⁡
ω
T
ω
0
cos
⁡
ω
T
-
1
ω
0
cos
⁡
ω
T
0
-
sin
⁡
ω
T
0
1
-
cos
⁡
ω
T
ω
1
sin
⁡
ω
T
ω
0
sin
⁡
ω
T
0
cos
⁡
ω
T
,
P
0
=
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
,
Q
=
0.05
0
0.1
0
0
0.05
0
0.1
,
R
=
5
0
0
5
.The other parameters are in Table 1.Table 1
Target and observation station state.
Target
Observation
Initialize position (m)
1
00
,
100
50
,
150
Initialize velocity (m/s)
2.8
0.3
Initialize heading (°)
45
160
Motion mode
CAV
CV


### 4.1.2. Simulation Results and Analysis

We compare the performances of GMEKF, GMUKF, GMCKF, and IGMCKF algorithms with the conditions mentioned above. Figure1 shows the results of the trajectory estimations of the target from the different algorithms versus the truth trajectory. Figure 2 shows the position error of different algorithms, and Figure 3 compares the position errors in X direction and Y direction among the different algorithms.Figure 1
Comparison of the trajectory estimations.Figure 2
The curve of position error.Figure 3
The curve of position error inX and Y.In Figure1, the red dotted-dashed line denotes the trajectory of observation station, and the arrow denotes the movement heading.From Figure2, it can be observed that the GMEKF and GMUKF have bad stability, and the error can be up to 25 m. The performance improvement of GMCKF is much better, and the error could be up to 13 m. The IGMCKF has the best tracking performance with an error up to 10 m, and the error could be limited to 3 m during stabilization.Figure4 shows the tracking performance of target velocity of different algorithms. Figure 5 shows the tracking performance of velocity components in X direction and Y direction, respectively.Figure 4
The tracking curve of target velocity.Figure 5
The tracking curve of velocity components.From Figure4, it can be observed that the GMCKF algorithm could track the target velocity after about 30 s, while the other algorithms need about 80 s.Figure6 shows the comparison of velocity errors of different algorithms.Figure 6
The curve of velocity error.Figure7 shows the comparison of target headings of different algorithms.Figure 7
The estimation curve of target heading.The accuracy achieved by the filtering methods was analyzed using the following metrics.The Root-Mean-Square Error (RMSE) is defined as(36)
RMSE
=
1
t
step
∑
k
=
1
t
step
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
,where t
step is the simulation steps and x
k
,
y
k and x
^
k
,
y
^
k are the true position and estimated position, respectively, of target at time k.The average Noncredibility Index (NCI) [25, 26] is defined as(37)
NCI
=
1
t
step
∑
k
=
1
t
step
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
P
p
,
k
-
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
RMSE
,where P
p
,
k is the uncertainty information of the position estimation and P
p
,
k
=
P
x
,
k
2
+
P
y
,
k
2.Table2 gives the NCI, RMSE, and the computing time cost of different tested algorithms.Table 2
RMSE and time of different algorithms.
NCI
RMSE (m)
Time (ms)
GMEKF
1.1534
6.3300
36.105
GMUKF
0.7591
5.4350
130.631
GMCKF
0.6522
3.9051
163.330
IGMCKF
0.6280
2.0042
257.302From Table2, we can see that IGMCKF has the smallest RMSE which illuminates that the algorithm has the highest precision. The running time of IGMCKF is the longest because IGMCKF needs to update the weight during the time update stage and split the observation noise into several Gaussian components and merge them based on CKF estimating the cubature points. But the running time could satisfy the requirement of real-time target tracking. The NCI of the algorithms except the GMEKF are below 1, which illuminate that the credibility of the various filter algorithms is high.


## 4.1.1. Simulation Parameter

The proposed algorithm is simulated by Monte Carlo withMC
=
50, the simulation time is Time
=
300
s, simulation step size is T
=
1
s, the covariance of observation noise in (2) is u
1
=
u
2
=
0, R
1
=
R, R
2
=
100
R
1, and the glint frequency factor ε
=
0.1; the other filter parameters are set as follows:(35)
F
=
1
sin
⁡
ω
T
ω
0
cos
⁡
ω
T
-
1
ω
0
cos
⁡
ω
T
0
-
sin
⁡
ω
T
0
1
-
cos
⁡
ω
T
ω
1
sin
⁡
ω
T
ω
0
sin
⁡
ω
T
0
cos
⁡
ω
T
,
P
0
=
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
,
Q
=
0.05
0
0.1
0
0
0.05
0
0.1
,
R
=
5
0
0
5
.The other parameters are in Table 1.Table 1
Target and observation station state.
Target
Observation
Initialize position (m)
1
00
,
100
50
,
150
Initialize velocity (m/s)
2.8
0.3
Initialize heading (°)
45
160
Motion mode
CAV
CV


## 4.1.2. Simulation Results and Analysis

We compare the performances of GMEKF, GMUKF, GMCKF, and IGMCKF algorithms with the conditions mentioned above. Figure1 shows the results of the trajectory estimations of the target from the different algorithms versus the truth trajectory. Figure 2 shows the position error of different algorithms, and Figure 3 compares the position errors in X direction and Y direction among the different algorithms.Figure 1
Comparison of the trajectory estimations.Figure 2
The curve of position error.Figure 3
The curve of position error inX and Y.In Figure1, the red dotted-dashed line denotes the trajectory of observation station, and the arrow denotes the movement heading.From Figure2, it can be observed that the GMEKF and GMUKF have bad stability, and the error can be up to 25 m. The performance improvement of GMCKF is much better, and the error could be up to 13 m. The IGMCKF has the best tracking performance with an error up to 10 m, and the error could be limited to 3 m during stabilization.Figure4 shows the tracking performance of target velocity of different algorithms. Figure 5 shows the tracking performance of velocity components in X direction and Y direction, respectively.Figure 4
The tracking curve of target velocity.Figure 5
The tracking curve of velocity components.From Figure4, it can be observed that the GMCKF algorithm could track the target velocity after about 30 s, while the other algorithms need about 80 s.Figure6 shows the comparison of velocity errors of different algorithms.Figure 6
The curve of velocity error.Figure7 shows the comparison of target headings of different algorithms.Figure 7
The estimation curve of target heading.The accuracy achieved by the filtering methods was analyzed using the following metrics.The Root-Mean-Square Error (RMSE) is defined as(36)
RMSE
=
1
t
step
∑
k
=
1
t
step
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
,where t
step is the simulation steps and x
k
,
y
k and x
^
k
,
y
^
k are the true position and estimated position, respectively, of target at time k.The average Noncredibility Index (NCI) [25, 26] is defined as(37)
NCI
=
1
t
step
∑
k
=
1
t
step
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
P
p
,
k
-
10
ln
⁡
x
k
-
x
^
k
2
+
y
k
-
y
^
k
2
RMSE
,where P
p
,
k is the uncertainty information of the position estimation and P
p
,
k
=
P
x
,
k
2
+
P
y
,
k
2.Table2 gives the NCI, RMSE, and the computing time cost of different tested algorithms.Table 2
RMSE and time of different algorithms.
NCI
RMSE (m)
Time (ms)
GMEKF
1.1534
6.3300
36.105
GMUKF
0.7591
5.4350
130.631
GMCKF
0.6522
3.9051
163.330
IGMCKF
0.6280
2.0042
257.302From Table2, we can see that IGMCKF has the smallest RMSE which illuminates that the algorithm has the highest precision. The running time of IGMCKF is the longest because IGMCKF needs to update the weight during the time update stage and split the observation noise into several Gaussian components and merge them based on CKF estimating the cubature points. But the running time could satisfy the requirement of real-time target tracking. The NCI of the algorithms except the GMEKF are below 1, which illuminate that the credibility of the various filter algorithms is high.


## 4.2. Simulation 2

### 4.2.1. Simulation Parameter

We will compare the performance of proposed algorithm (interacting multiple model CKF, IMMCKF) and interacting multiple model Gaussian mixture CKF (IMMGMCKF) algorithms. The filter parameters are set as follows. The initial position of target is20,150, the initial velocity is 0
,
-
1.5
m
/
s, and the acceleration is 0
m
/
s
2. The target moves with constant initial velocity in a straight line from 0 s to 100 s. Then it maneuvers and moves with constant acceleration 0.0075,0.0075
m
/
s
2 from 101 s to 200 s. From 201 s to 300 s, it moves with constant velocity 0.7575
,
-
0.7425
m
/
s. From 301 s to 380 s, it maneuvers with constant acceleration -
0.025
,
-
0.015
m
/
s
2, and it moves with constant velocity -
1.2175
,
-
1.9275
m
/
s for the last 20 s.The observation noise is set the same as simulation 1. Singer model is employed for the IGMCKF algorithm, the state vector isx
=
x
,
x
˙
,
x
¨
,
y
,
y
˙
,
y
¨, and other filter parameters are set as follows:(38)
F
=
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
0
0
0
0
0
0
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
,
P
0
=
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
0
0
0
0
0
0
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
,
Q
=
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
0
0
0
0
0
0
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
,
R
=
10
0
0
10
.For the IMMCKF and IMMGMCKF algorithms, three models are employed: two constant acceleration models with different acceleration and constant velocity model. The model initial distribution isu
=
0.8,0.1,0.1, and Markov transition matrix is as follows:(39)
Π
=
0.95
0.025
0.025
0.025
0.95
0.025
0.025
0.025
0.95
.


### 4.2.2. Simulation Results and Analysis

The simulation results are obtained from 50 Monte Carlo runs. Figure8 shows the position trajectory of various algorithms. Figures 9 and 10 denote the RMSE of position tracking in X direction and Y direction, respectively. Figures 11 and 12 denote the RMSE of velocity tracking in X direction and Y direction, respectively.Figure 8
Comparison of the trajectory estimations.Figure 9
X position RMSE.Figure 10
Y position RMSE.Figure 11
X velocity RMSE.Figure 12
Y velocity RMSE.From the above simulation results and Table3, it can be seen that the IMMCKF has the worst accuracy of estimation, and because of that the observation noise is modelled as glint noise, and the IMMCKF algorithm is suitable for the Gaussian noise. The IMMGMCKF and IGMCKF have the similar accuracy of estimation and because of that the Gaussian mixture algorithm could deal with the non-Gaussian noise well. The running time of the IMMGMCKF is the longest and because of that the algorithm needs to switch between different motion models. The proposed algorithm has the best estimation accuracy and best computation complexity and could satisfy the real-time target tracking requirement.Table 3
Comparisons of average RMSE and running time of different algorithms.
X position (m)
Y position (m)
X velocity (m/s)
Y velocity (m/s)
Running time (ms)
IMMCKF
0.6713
0.6909
0.4741
0.4930
222
IMMGMCKF
0.1074
0.1238
0.0909
0.0920
315
IGMCKF
0.0983
0.0895
0.0463
0.0462
182


## 4.2.1. Simulation Parameter

We will compare the performance of proposed algorithm (interacting multiple model CKF, IMMCKF) and interacting multiple model Gaussian mixture CKF (IMMGMCKF) algorithms. The filter parameters are set as follows. The initial position of target is20,150, the initial velocity is 0
,
-
1.5
m
/
s, and the acceleration is 0
m
/
s
2. The target moves with constant initial velocity in a straight line from 0 s to 100 s. Then it maneuvers and moves with constant acceleration 0.0075,0.0075
m
/
s
2 from 101 s to 200 s. From 201 s to 300 s, it moves with constant velocity 0.7575
,
-
0.7425
m
/
s. From 301 s to 380 s, it maneuvers with constant acceleration -
0.025
,
-
0.015
m
/
s
2, and it moves with constant velocity -
1.2175
,
-
1.9275
m
/
s for the last 20 s.The observation noise is set the same as simulation 1. Singer model is employed for the IGMCKF algorithm, the state vector isx
=
x
,
x
˙
,
x
¨
,
y
,
y
˙
,
y
¨, and other filter parameters are set as follows:(38)
F
=
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
0
0
0
0
0
0
1
T
T
2
2
0
0
0
0
1
T
0
0
0
0
0
1
,
P
0
=
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
0
0
0
0
0
0
10
10
10
0
0
0
10
20
30
0
0
0
10
30
60
,
Q
=
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
0
0
0
0
0
0
1
20
1
8
1
6
0
0
0
1
8
1
3
1
2
0
0
0
1
6
1
2
1
,
R
=
10
0
0
10
.For the IMMCKF and IMMGMCKF algorithms, three models are employed: two constant acceleration models with different acceleration and constant velocity model. The model initial distribution isu
=
0.8,0.1,0.1, and Markov transition matrix is as follows:(39)
Π
=
0.95
0.025
0.025
0.025
0.95
0.025
0.025
0.025
0.95
.


## 4.2.2. Simulation Results and Analysis

The simulation results are obtained from 50 Monte Carlo runs. Figure8 shows the position trajectory of various algorithms. Figures 9 and 10 denote the RMSE of position tracking in X direction and Y direction, respectively. Figures 11 and 12 denote the RMSE of velocity tracking in X direction and Y direction, respectively.Figure 8
Comparison of the trajectory estimations.Figure 9
X position RMSE.Figure 10
Y position RMSE.Figure 11
X velocity RMSE.Figure 12
Y velocity RMSE.From the above simulation results and Table3, it can be seen that the IMMCKF has the worst accuracy of estimation, and because of that the observation noise is modelled as glint noise, and the IMMCKF algorithm is suitable for the Gaussian noise. The IMMGMCKF and IGMCKF have the similar accuracy of estimation and because of that the Gaussian mixture algorithm could deal with the non-Gaussian noise well. The running time of the IMMGMCKF is the longest and because of that the algorithm needs to switch between different motion models. The proposed algorithm has the best estimation accuracy and best computation complexity and could satisfy the real-time target tracking requirement.Table 3
Comparisons of average RMSE and running time of different algorithms.
X position (m)
Y position (m)
X velocity (m/s)
Y velocity (m/s)
Running time (ms)
IMMCKF
0.6713
0.6909
0.4741
0.4930
222
IMMGMCKF
0.1074
0.1238
0.0909
0.0920
315
IGMCKF
0.0983
0.0895
0.0463
0.0462
182


## 5. Conclusion

An alternative efficient nonlinear filtering algorithm is designed in this paper, which is an improved Gaussian mixture filter within the Cubature Kalman Filter frame. The algorithm could deal with the nonlinear filtering problem with non-Gaussian noise using Gaussian mixture density approximating the observation noise and update the weight adaptively with integral square difference during the time update stage. The feasibility of the algorithm was proven using a simulation based on the target tracking on offing, and the simulation shows that the proposed algorithm has higher estimation precision and faster tracking response.



---

*Source: 1082837-2016-07-17.xml*