Kalman-Gain Aided Particle
PHD Filter for Multitarget
Tracking
ABDULLAHI DANIYAN , Student Member, IEEE
YU GONG, Member, IEEE
SANGARAPILLAI LAMBOTHARAN, Senior Member, IEEE
Loughborough University, Loughborough, U.K.
PENGMING FENG, Member, IEEE
JONATHON CHAMBERS, Fellow, IEEE
Newcastle University, Newcastle upon Tyne, U.K.
We propose an efficient sequential Monte Carlo probability hy-
pothesis density (PHD) filter which employs the Kalman-gain ap-
proach during weight update to correct predicted particle states by
minimizing the mean square error between the estimated measure-
ment and the actual measurement received at a given time in order
to arrive at a more accurate posterior. This technique identifies and
selects those particles belonging to a particular target from a given
PHD for state correction during weight computation. Besides the im-
proved tracking accuracy, fewer particles are required in the proposed
approach. Simulation results confirm the improved tracking perfor-
mance when evaluated with different measures.
Manuscript received November 11, 2015; revised September 19, 2016 and
December 18, 2016; released for publication March 24, 2017. Date of
publication April 5, 2017; date of current version October 10, 2017.
DOI. No. 10.1109/TAES.2017.2690530
Refereeing of this contribution was handled by B.-N. Vo.
This work was supported in part by the Engineering and Physical Sciences
Research Council under Grant EP/K014307/1, in part by the MOD Uni-
versity Defence Research Collaboration in Signal Processing, U.K., and
in part by the Petroleum Technology Development Fund, Nigeria.
Authors’ addresses: A. Daniyan, Y. Gong, and S. L ambotharan are
with the School of Electronic, Electrical and Systems Engineering,
Loughborough University, Loughborough LE11 3TU, U.K., E-mail:
(a.daniyan@lboro.ac.uk; y.gong@lboro.ac.uk; s.lambotharan@lboro.ac.
uk); P. Feng and J. Chambers are with the School of Electrical and
Electronic Engineering, Newcastle University, Newcastle upon Tyne
NE1 7RU, U.K., E-mail: (p.feng2@newcastle.ac.uk; jonathon.chambers@
newcastle.ac.uk). (Corresponding author: Abdullahi Daniyan.)
0018-9251
C
2017 CCBY
I. INTRODUCTION
Multitarget tracking (MTT) is essential in many appli-
cation areas, such as motion-based recognition, automated
security, navigation and surveillance, medical imaging, traf-
fic control, and human–computer interaction [1]–[3]. MTT
belongs to a class of dynamic state estimation problems
[3]–[5]. In MTT, targets can appear and disappear randomly
in time and this results in a varying and unknown number
of targets and their corresponding states. Furthermore, not
all measurements received by sensors at each time instance
are due to existing targets. The sensor may pick up de-
tections as false alarms due to clutter or may even miss
some detections. As a result, the measurements received at
each time step are corrupted and consist of indistinguish-
able measurements that may be either target originated or
due to clutter. Therefore, the main objective of MTT is to be
able to jointly estimate target states and number of targets
from a set of corrupted observations.
Furthermore, because there is no particular ordering
between measurements received and target states at each
time step in terms of association, both the received mea-
surements and target states can be represented as finite
sets [6]–[9]. The modeling of target states and observa-
tions as a random finite set (RFS) allows for the use of
the Bayesian filtering approach (as an optimal multitar-
get filter) to estimate the multitarget states in the pres-
ence of clutter, missed detections, and association uncer-
tainty [6]–[9]. Tractable alternatives to the optimal multi-
target filters include the RFS-based probability hypothesis
density (PHD) filter, the cardinalized PHD (CPHD) filter
[8]–[10], the multitarget multi-Bernoulli (MeMBer) filter,
and its cardinality-balanced version, the CBMeMBer filter
[8], [9], [11]. Both the CPHD in [10] and the CBMeM-
Ber in [11] have been shown to have better performance
than the MeMBer filter in [8]. The CBMeMBer filter was
proposed specifically to address the pronounced bias in the
cardinality estimate of the MeMBer filter. For more details
on other tractable RFS-based MTT methods, the reader is
referred to [12]–[15].
The PHD filter is a recursion that propagates the
posterior intensity of the RFS of targets in time [6]. The
integral of the PHD is the expected number of targets in a
measurable region, and the peaks of the PHD function pro-
vide the estimates of the target states [6], [8], [9]. The PHD
filter is able to track time-varying multiple targets without
the need to explicitly associate measurements to tracks.
In the literature, the PHD filter has been implemented in
two distinct fashions, that is, as the Gaussian mixture PHD
(GM-PHD) filter [16] and the sequential Monte Carlo
PHD (SMC-PHD) filter [17]. In the GM-PHD filter imple-
mentation, the PHD is assumed to be a GM while in the
SMC-PHD filter implementation, the PHD is approximated
by a set of weighted particles and does not need any further
1
The MeMBer filter is a recursion that propagates (approximately) the
multitarget posterior density and is based on the assumption that every
multitarget posterior is a MeMBer process [8], [11].
IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 53, NO. 5 OCTOBER 2017 2251
assumptions. The SMC-PHD filter is therefore more
suitable for tracking in nonlinear and non-Gaussian
environments.
In SMC filter design, the choice of importance density
function from which samples are drawn to avoid sample de-
generacy and impoverishment is of crucial importance [18].
Furthermore, in MTT which involves multiple modalities,
if particles are in clusters representing the modes of the pos-
terior, the iterative process of randomly drawing samples
from proposal distributions results in random fluctuations
in the total weight attributed to each mode [19]. In addition,
the errors associated with the estimation of the weights of
each mode will increase in magnitude with time [19]. These
errors arise due to the stochastic nature of drawing samples
from the proposal distribution and the stochasticity of
the resampling process [19]. These two processes greatly
influence the performance of SMC filters. SMC filters are
further affected by how well the state space of targets is
populated with samples. In addition, Vo et al. [20] argued
that the mean squared error (MSE) of the SMC-PHD
filter is inversely proportional to the number of samples.
In [18], it is shown that the optimal importance density
function is the posterior. In many cases it is difficult to
sample from the optimal importance density. As an attempt
to solve the importance sampling problem, Maskell and
Julier [19] proposed an optimized proposal distribution for
SMC filters with multiple modes in general. However, this
approach tends to be problem specific. In [21], Yoon et al.
proposed the Gaussian mixture unscented sequential Monte
Carlo probability hypothesis density (GM-USMC-PHD)
filter which uses the GM representation to approximate the
importance sampling f unction and the predictive density
functions via the unscented information filter (UIF). Addi-
tionally, in [22] and [23] the auxiliary SMC-PHD filter and
its improved version, the a uxiliary particle PHD (AP-PHD)
filter are proposed, respectively. Both try to use the
auxiliary particle approach to incorporate the measurement
into the importance sampling function. This, however,
involves double computation on the measurement and more
samples are required to populate the state space in order to
make the importance sampling function more viable.
However, it is also possible to construct suboptimal
approximations to the optimal importance density by using
local linearization techniques [18]. As a realization of this,
the unscented Kalman particle PHD filter was proposed in
[24] for the joint tracking of multitargets. It tries to use the
unscented Kalman filter (UKF) in the prediction step. This
allowed for the inclusion of the latest measurement to draw
particles. Similarly, Ma et al. [25] proposed the Kalman
particle PHD filter for multitarget visual tracking, which
uses the Kalman filter to construct the proposal density
also in the prediction step. Furthermore, Tang et al. [26]
presented an improvement to the SMC-PHD filter, which
incorporates the latest measurements into the resampling
step by using the UKF.
Additionally, there are other works in the literature that
combine the implementation of the GM and particle PHD
filter into one filter such as the GM particle PHD filter in
[27] and [29] and the Gaussian mixture SMC-PHD in [30].
These methods attempt to combine the advantages of both
GM-PHD and SMC-PHD filters. The methods give some
level of performance improvement without easing com-
putational burden or the number of particles. Also, it may
be possible to implement the Markov Chain Monte Carlo
(MC) sampling method in the update stage of the SMC-
PHD filter as a way of asymptotically approximating the
posterior. However, this approach will require even more
particles, as these extra particles will be used to perform
some sort of random walk in order to achieve maximum a
posteriori estimate of target states but no guarantees exist
about it yielding good point estimates [31]. Recently, Zheng
et al. [32] proposed a data-driven SMC-PHD filter for MTT.
The method tries to segment the measurements available at
each time step into measurements due to persistent targets
and measurements due to newborn targets. Again this does
not help to reduce the number of particles but rather, more
particles are required to populate regions of interest.
It is desirable, therefore, to have an efficient filter that
can provide for particle state correction for any proposal
distribution using fewer particles. This gives the motivation
for the Kalman-gain aided sequential Monte Carlo proba-
bility hypothesis density (KG-SMC-PHD) filter. The KG-
SMC-PHD filter provides for the particle state correction of
the predicted multitarget state. This is achieved with the ap-
plication of the Kalman state update technique on selected
particles to minimize the MSE between the estimated mea-
surements and actual measurement.
In this paper, we propose an SMC-PHD filter with a
validation threshold to select promising particles and to
guide them to regions of high likelihood using the Kalman
gain, irrespective of the importance density function. This
method seeks to minimize the MSE between the estimated
measurements due to selected particles and the actual mea-
surements to achieve a more efficient SMC-PHD filter with
less computational complexity. This allows fewer particles
to be used to populate the state space and at the same time
achieve improved tracking performance as opposed to the
standard SMC-PHD filter.
The remainder of this paper is organized as follows. In
Section II, the MTT problem is presented in terms of process
and measurement models. Section III presents the idea of
the importance density function and highlights some com-
mon choices of proposal distributions. In Section IV, the
PHD filter recursion is presented and explained followed by
a description of the standard SMC-PHD filter implemen-
tation. Next, Section V presents our proposed KG-SMC-
PHD filter. Simulation results together with discussions are
presented in Section VI. Finally, conclusions are drawn in
Section VII.
II. MTT PROBLEM FORMULATION
The MTT problem relates to that of modeling a dy-
namical system. Two models are generally used, the state
evolution model and the measurement model.
2252 IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 53, NO. 5 OCTOBER 2017
A. State Model
A nonlinear system governed by the state evolution
model is considered
x
k
= f
k−1
(x
k−1
, v
k
)(1)
where x
k
denotes the tth target state at discrete time k, v
k
is an independent and identically distributed (i.i.d.) process
noise vector, and f
k−1
(·) is the nonlinear system transition
function. Then, the multitarget state at time k can be written
as
X
k
=
x
1,k
,...,x
T,k
∈ E
s
(2)
where T is the number of targets present at each time k,and
E
s
denotes the state space.
B. Measurement M odel
Let the multitarget cumulative measurement sequence
up to time K be Z
1:K
: Z
1
, Z
2
,...,Z
K
⊂ E
o
. Measure-
ments consist of both target-originated measurements and
false alarms due to clutter. Then, the multitarget measure-
ment set at time k in the observation space is
Z
k
=
z
1,k
,...,z
α,k
c
1,k
,...c
β,k
⊂ E
o
(3)
where
z
1,k
,...,z
α,k
denotes the target-originated
measurement set with number of measurements α;
c
1,k
,...c
β,k
denotes the false measurement set with the
number of measurements β,andE
o
denotes the observa-
tion space. The tth target-originated nonlinear measurement
model is given as
z
k
= h
k
(x
k
, n
k
)(4)
where h
k
(·) is a nonlinear function and n
k
is an i.i.d. process
noise vector.
III. IMPORTANCE DENSITY FUNCTION
In this section, we focus on proposal distributions and
their role in S MC methods in general.
A. Importance Sampling
MC methods for numerical integration deal with prob-
lems of the form
g =
n
f (y)π(y)dy (5)
where π(y) is such that π(y) ≥ 0 and integrates to unity
n
π(y)dy = 1(6)
is a pdf.
It is also the assumption that it is possible to generate
N 1 samples distributed according to the probability
density π (y). The MC estimate of the integral (5) is formed
by taking the average over the set of samples
ˆ
g =
1
N
N
i=1
f (y
i
)(7)
where N is assumed to be large. However, π (y) is not
usually a familiar density and so it is difficult to generate
samples directly from it. When the latter is the case, the
integral of (5) can be solved by letting q(y)beaproposal
distribution or importance density which is easy to generate
samples and with the assumption that π(y) > 0 ⇒q(y) > 0
for all y ∈
n
. Under this assumption, (5) becomes
g =
n
f (y)π(y)dy =
n
f (y)
π(y)
q(y)
q(y)dy. (8)
An MC estimate is then computed by generating N 1
samples from q(y) instead of π(y) and forming a weighted
sum
ˆ
g =
1
N
N
i=1
f (y
i
)w(y
i
)(9)
where w(y
i
) ∝
π(y
i
)
q(y
i
)
are the associated weights [33]. To
sum up, given a distribution π(y) that is difficult to sample
from, importance sampling facilitates sampling from π (y)
by sampling from an a lternate distribution q(y) known as
proposal distribution but weighting appropriately.
B. Importance Densities
Some common choices of importance density in SMC
methods are given below.
1) Transitional Prior (TP): This is the most popular
choice of suboptimal proposal distribution for SMC-PHD
filters and particle filters in general because its implemen-
tation is easy and straightforward [34]. This choice requires
sampling from the dynamic prior, i.e.,
q
x
k
|x
l
k−1
, z
k
= p
x
k
|x
l
k−1
. (10)
2) Extended Particle Filter (EPF): Given that the
measurement model of (4) is nonlinear, but Gaussian, it is
possible to use a proposal distribution that exploits a linear
approximation to the posterior [19] in the same way as
the extended Kalman filter (EKF) uses a local linearization
about its estimates. The proposal distribution is then given
as
q
x
k
|x
l
k−1
, z
k
=
N
(
x
k
; u
k
, A
k
)
(11)
where
u
k
= f
k−1
(x
k−1
) + A
k
H
T
k
R
−1
k
(
z
k
− h(f
k−1
(x
k−1
))
)
(12)
H
k
=
∂h
∂x
k
f
k−1
(x
k−1
)
(13)
where A
k
and R
k
denote state and measurement covari-
ances, respectively, and H
k
is the measurement transforma-
tion matrix.
3) Unscented Particle Filter (UPF): As an alter-
native to the EPF, an unscented transform can be used to
calculate the mean h(f
k−1
(x
k−1
)) and covariance H
k
by gen-
erating sigma points and applying a transform such that the
new generated samples have f
k−1
(x
k−1
) as mean and P
k−1
as covariance. h(f
k−1
(x
k−1
)) is then evaluated at each sigma
point and H
k
computed from these samples [19].
DANIYAN ET AL.: KALMAN-GAIN AIDED PARTICLE PHD FILTER FOR MULTITARGET TRACKING 2253
Fig. 1. Schematic representation of the standard SMC-PHD filter showing the 2-D state space of the PHD of two targets populated with particles.
The contours represent the state space of targets. The contour centers and number of centers represent the mode and cardinality of targets, respectively.
Boxes A, B,andC represent various stages of the filter. The square-shaped and diamond-shaped particles are for target 1 and target 2, respectively.
The colours stand for different particle states. The particles marked with “
√
”inB denote particles with higher weight for when the latest observation
arrives.
IV. PROBABILITY HYPOTHESIS DENSITY
A. PHD Filter
The PHD D
of a given RFS is the first-order moment
of and is given by [6], [8], [9] the following equation:
D
(x) = E
{
δ
(x)
}
=
δ
X
(x)P
(dX) (14)
where E
{
·
}
is the statistical expectation operator and
δ
(x) =
y∈
δ
y(x)
is the random density representation
of . P
is the probability measure of the RFS. The PHD
filter is a recursion of the PHD, D
k|k
that is associated with
the multitarget posterior density p(X
k
|Z
k
), and
p
(
X
k
|Z
k
)
∝ p
(
Z
k
|X
k
)
p
(
X
k
|Z
k−1
)
(15)
where p(Z
k
|X
k
)andp(X
k
|Z
k−1
) denote the multitarget
likelihood and prior density, respectively.
The prediction formula of the PHD, D
k|k
is given as [8],
[9] follows:
D
k|k−1
(x
k
|Z
k−1
) = γ
k
(x
k
)
+
φ
k|k−1
(x
k
, x
k−1
)D
k−1|k−1
(x
k−1
|Z
k−1
)dx
k−1
(16)
with the factor
φ
k|k−1
(x
k
, x
k−1
) = p
S
(x
k−1
)f
k|k−1
(x
k
, x
k−1
)
+ b
k|k−1
(x
k
, x
k−1
) (17)
where γ
k
(·) is the PHD of the spontaneous birth, p
S
(·)is
the probability of the target survival, f
k|k−1
(x
k
, x
k−1
)isthe
single target motion model, and b
k|k−1
(x
k
, x
k−1
)isthePHD
of the spawned targets.
The update formula is given as
D
k|k
(x
k
|Z
k
)
=
⎡
⎣
ν(x
k
) +
z∈Z
k
ψ
k,z
(x
k
)
κ
k
(z) +D
k|k−1
,ψ
k,z
⎤
⎦
D
k|k−1
(x
k
|Z
k−1
)
(18)
with ν(x
k
) = 1 − p
D
(x
k
), ψ
k,z
(x
k
) = p
D
(x
k
)g(z|x
k
), and
κ
k
(z) = λ
k
c
k
(z); where p
D
(x
k
)andν(x
k
) denote the proba-
bility of target detection and nondetection for a given (x
k
),
respectively, g(z|x
k
) is the measurement likelihood func-
tion for the single target, κ
k
(z) is the clutter intensity, λ
k
is
the average number of Poisson clutter points per scan, and
c
k
(z) is the probability density over the state-space of the
clutter point; ·, · denotes inner product and is computed
as [8], [9]
D
k|k−1
,ψ
k,z
=
D
k|k−1
(x
k
|Z
k−1
)ψ
k,z
(x
k
)dx
k
. (19)
B. Standard SMC-PHD Filter
The PHD filter can be implemented either as in the SMC
fashion (particle-PHD) or as the GM-PHD. The SMC-PHD
filter approximates the PHD using random samples and
is more specifically an effective scheme in nonlinear and
non-Gaussian scenarios as well as different noise models
[35]. For comparison purposes, the standard SMC-PHD
filter of [17] is briefly presented. The implementation of
the standard SMC-PHD filter usually requires four stages.
These stages are briefly presented in Algorithm 1.
Fig. 1 illustrates how particles are used to represent
and track targets in the standard SMC-PHD filter. The state
2254 IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 53, NO. 5 OCTOBER 2017
Algorithm 1: The Standard SMC-PHD Filter.
1: at k =0,Initialize
x
l
k
,w
l
k
L
k
l=1
.
2: for k = 1:K do
3: Prediction
4: for l = 1:L
k
do
5: Draw samples for existing targets,
˜
x
l
k|k−1
∼
q
k
(·|
˜
x
l
k−1
, Z
k
),
and compute weights,
˜w
l
k|k−1
=
φ
k|k−1
(
˜
x
l
k
,
˜
x
l
k−1
)
q
k
(
˜
x
k|k−1
|
˜
x
l
k−1
,Z
k
)
w
l
k−1
6: end for
7: for l = L
k
+ 1:L
k
do
8: Draw samples for newborn targets,
˜
x
l
k|k−1
∼
p
k
(·|Z
k
),
and compute weights, ˜w
l
k|k−1
=
γ
k
(
˜
x
l
k
)
J
k
p
k
(
˜
x
k|k−1
|Z
k
)
9: end for
10: Update
11: for z ∈ Z
k
do
12: C
k
(z) =
L
k
l=1
p
D
(
˜
x
l
k|k−1
)g(z|
˜
x
l
k|k−1
)˜w
l
k|k−1
13: for l = 1:L
k
do
14: update weight,
˜w
i
k
=
ν +
z∈Z
k
p
D
(
˜
x
l
k|k−1
)g(z|
˜
x
l
k|k−1
)
κ
k
(z)+C
k
(z)
˜w
l
k|k−1
ν = 1 − p
D
(
˜
x
i
k|k−1
)
15: end for
16: end for
17: Resample
18: Compute estimated number of targets,
ˆ
T
k|k
= round
L
k
l=1
˜w
l
k
19: Resample L
k
particles using resampling
techniques such as in [34].
20: return
˜
x
l
k|k−1
,
ˆ
T
k|k
L
k
L
k
l=1
≡
x
l
k
,w
l
k
L
k
l=1
21: end for
space of two targets populated with particles at time k is
shown. In A, during the prediction stage, the PHD is rep-
resented with eight equally weighted particles. In B,as
the latest measurement arrives, the particle weights are up-
dated accordingly. Particles with higher weights are chosen
for resampling. As seen in B, the highly weighted parti-
cles are marked with “
√
,” respectively, five particles for
the first target and six particles for the second target. To
ensure that the number of particles remains eight for each
target, the particles marked with “
√
” are resampled de-
pending on the size of their weights as seen in C. Notice
that the particle positions remain unchanged and the parti-
cles corresponding to high weights are retained and those
with lower weights are discarded. The estimated state of
the targets or the posterior at time k is derived from the re-
sampled particles. It is true that populating the state space
of the targets with many more particles will result in more
particles falling near the modes of the state space. This will
translate to higher weighted particles and a more accurate
posterior. However, doing this will increase computational
complexity.
In the next section, the proposed SMC-PHD filter is
presented.
V. PROPOSED SMC-PHD FILTER
In the standard SMC-PHD filter, the particles appear to
be scattered and it is difficult to guide particles to regions
of interest. The filter’s ability to estimate the posterior at
a given time depends on how densely the state space is
populated with samples and how well the estimated mea-
surements match the actual measurements received in that
time frame. The weights are then updated accordingly. The
SMC-PHD filter does not provide for particle state cor-
rection to achieve particle improvement. In other words, it
does not seek to reduce the error between the actual mea-
surement and the estimated measurements irrespective of
the importance density chosen. The proposed method seeks
to address this problem. The novelty of our approach lies
in the technique behind the Kalman filter. The Kalman fil-
ter is a minimum MSE estimator, which in effect seeks to
recursively minimize the MSE between the estimated mea-
surements and actual measurements using the Kalman gain
[36]. The Kalman gain computes the required correction
from the observation and transforms the correction of the
observation back to the correction of state. The proposed ap-
proach tries to apply particle state correction/improvement
using the Kalman gain to guide validated particles in the
SMC-PHD filter to the region of higher likelihood to better
approximate the posterior at each time step.
A. Measurement Set Partition
Given that T
k
targets exist at time k, the measurements
received at k may consist of target-originated measure-
ments (i.e., measurements due to persistent target or
newborn targets) and clutter. In the standard SMC-PHD
filter, all measurements are used to compute weights to
show the significance of all particles with no attempt to
check for errors. Therefore, a measurement set partition
is needed to separate the measurement set into target-
originated measurements and measurements due to clutter.
We use a statistical distance measure and gating technique
to achieve this. The second step is to identify promising
particles from the predicted target state using a validation
threshold and improve their states using the Kalman gain
while updating weights as measurement arrives.
At time k, measurements assumed to originate from
persistent targets are identified by computing the square
Mahalanobis distance between elements in the measure-
ment set Z
k−1
at time k − 1andZ
k
at time k from (3)
as
d
2
i,j,k
=
z
i
k
− z
j
k−1
T
−1
k
z
i
k
− z
j
k−1
(20)
for i = 1,...,
|
Z
k
|
and j = 1,...,
|
Z
k−1
|
.
k
is the mea-
surement covariance matrix. For target-originated measure-
ments z
i
k
and z
j
k−1
belonging to the same target, the square
Mahalanobis distance d
2
i,j,k
is χ
2
distributed with degree of
freedom equal to the dimension of the measurement vector.
Therefore, a unit-less threshold
˜
d can be computed for a
DANIYAN ET AL.: KALMAN-GAIN AIDED PARTICLE PHD FILTER FOR MULTITARGET TRACKING 2255