IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 50, NO. 2, FEBRUARY 2002 309
Sequential Monte Carlo Methods for Multiple Target
Tracking and Data Fusion
Carine Hue, Jean-Pierre Le Cadre, Member, IEEE, and Patrick Pérez
Abstract—The classical particle filter deals with the estimation
of one state process conditioned on a realization of one observation
process. We extend it here to the estimation of multiple state pro-
cesses given realizations of several kinds of observation processes.
The new algorithm is used to track with success multiple targets in
a bearings-only context, whereas a JPDAF diverges. Making use
of the ability of the particle filter to mix different types of observa-
tions, we then investigate how to join passive and active measure-
ments for improved tracking.
Index Terms—Bayesian estimation, bearings-only tracking,
Gibbs sampler, multiple receivers, multiple targets tracking,
particle filter.
I. INTRODUCTION
M
ULTITARGET tracking (MTT) deals with the state esti-
mation of an unknown number of moving targets. Avail-
able measurements may both arise from the targets, if they are
detected, and from clutter. Clutter is generally considered to be
a model describing false alarms. Its (spatio–temporal) statistical
properties are quite different from those of the target, which
makes the extraction of target tracks from clutter possible. To
perform multitarget tracking, the observer has at his disposal a
huge amount of data, possibly collected on multiple receivers.
Elementary measurements are receiver outputs, e.g., bearings,
ranges, time-delays, Dopplers, etc.
The main difficulty, however, comes from the assignment of
a given measurement to a target model. These assignments are
generally unknown, as are the true target models. This is a neat
departure from classical estimation problems. Thus, two distinct
problems have to be solved jointly: the data association and the
estimation.
As long as the association is considered in a deterministic
way, the possible associations must be exhaustivelyenumerated.
This leads to an NP-hard problem because the number of pos-
sible associations increases exponentially with time, as in the
multiple hypothesis tracker (MHT) algorithm [28]. In the joint
probabilistic data association filter (JPDAF) [11], the associa-
tion variables are considered to be stochastic variables, and one
needs only to evaluate the association probabilities at each time
step. However, the dependence assumption on the associations
implies the exhaustive enumeration of all possible associations
Manuscript received January 31, 2001; revised October 11, 2001. The asso-
ciate editorcoordinating the review of this paper and approving it for publication
was Dr. Petar M. Djuric
´
.
C. Hue is with Irisa/Université de Rennes 1, Rennes, France (e-mail:
chue@irisa.fr).
J.-P. Le Cadre is with Irisa/CNRS, Rennes, France (e-mail: lecadre@irisa.fr).
P. Pérez is with Microsoft Research, Cambridge, U.K. (e-mail: pperez@mi-
crosoft.com).
Publisher Item Identifier S 1053-587X(02)00571-8.
at the current time step. When the association variables are in-
stead supposed to be statistically independent like in the prob-
abilistic MHT (PMHT [12], [32]), the complexity is reduced.
Unfortunately, the above algorithms do not cope with nonlinear
models and non-Gaussian noises.
Under such assumptions (stochastic state equation and non-
linear state or measurement equation non-Gaussian noises), par-
ticle filters are particularly appropriate. They mainly consist of
propagating a weighted set of particles that approximates the
probability density of the state conditioned on the observations.
Particle filtering can be applied under very general hypotheses,
is able to cope with heavy clutter, and is very easy to implement.
Such filters have been used in very different areas for Bayesian
filtering under different names: The bootstrap filter for target
tracking in [15] and the Condensation algorithm in computer vi-
sion [20] are twoexamples,among others. In theearliest studies,
the algorithm was only composed of two periods: The particles
were predicted according to the state equation during the pre-
diction step; then, their weights were calculated with the likeli-
hood of the new observation combined with the former weights.
A resampling step has rapidly been added to dismiss the parti-
cles with lower weights and avoid the degeneracy of the particle
set into a unique particle of high weight [15]. Many ways have
been developed to accomplish this resampling, whose final goal
is to enforce particles in areas of high likelihood. The frequency
of this resampling has also been studied. In addition, the use of
kernel filters [19] has been introduced to regularize the sum of
Dirac densities associated with the particles when the dynamic
noise of the state equation was too low [26]. Despite this long
history of studies, in which the ability of particle filter to track
multiple posterior modes is claimed, the extension of the par-
ticle filter to multiple target tracking has progressively received
attention only in the five last years. Such extensions were first
claimed to be theoretically feasible in [2] and [14], but the ex-
amples chosen only dealt with one single target. In computer
vision, a probabilistic exclusion principle has been developed
in [24] to track multiple objects, but the algorithm is very de-
pendent of the observation model and is only applied for two
objects. In the same context, a Bayesian multiple-blob tracker
(BraMBLe) [6] has just been proposed. It deals with a varying
number of objects that are depth-ordered thanks to a 3-D state
space. Lately, in mobile robotic [29], a set of particle filters for
each target connected by a statistical data association has been
proposed. We propose here a general algorithm for multitarget
tracking in the passive sonar context and take advantage of its
versatility to extend it to multiple receivers.
This work is organized as follows. In Section II, we recall
the principles of the basic particle filter with adaptive resam-
pling for a single target.We beginSection III with a presentation
of the multitarget tracking problem and its classical solutions
1053–587X/02$17.00 © 2002 IEEE
310 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 50, NO. 2, FEBRUARY 2002
and of the related work on multitarget tracking by particle fil-
tering methods. Then, we present our multitarget particle filter
(MTPF). The new algorithm combines the two major steps (pre-
diction and weighting) of the classical particle filterwithaGibbs
sampler-based estimation of the assignment probabilities. Then,
we propose to add two statistical tests to decide if a target has
appeared or disappeared from the surveillance area. Carrying
on with the approach of the MTPF, we finally present an exten-
sion to multireceiver data in the context of multiple targets (the
MRMTPF). Section IV is devoted to experiments. Simulation
results in a bearings-only context with a variable clutter den-
sity validate the MTPF algorithm. A comparison with a JPDAF
based on the extended Kalman filter establishes the superiority
of the MTPF for the same scenario. Then, a simulation with oc-
clusion shows how the disappearance of a target can be handled.
The suitable quantity and distribution of active measurements
are then studied in a particular scenario to improve the perfor-
mance obtained with passive measurements only.
As far as the notational conventions are concerned, we always
use the index
to refer to one among the tracked targets. The
index
designates one of the observationsobtainedat instant
. The index is devoted to the particles. The index is used
for indexing the iterations in the Gibbs sampler, and
is used
for the different receivers. Finally, the probability densities are
denoted by
if they are continuous and by if they are discrete.
II. B
ASIC SINGLE-TARGET PARTICLE FILTER
For the sake of completeness, the basic particle filter is now
briefly reviewed. We consider a dynamic system represented by
the stochastic process
, whose temporal evolution
is given by the state equation:
(1)
It is observed at discrete times via realizations of the stochastic
process
governed by the measurement model
(2)
The two processes
and in (1) and (2)
are only supposed to be independent white noises. Note that the
functions
and are not assumed linear. We will denote by
the sequence of the random variables and by
one realization of this sequence. Note that throughout the
paper, the first subscript of any vector will always refer to the
time.
Our problem consists of computing at each time
the condi-
tional density
of the state , given all the observations ac-
cumulated up to
, i.e., and
of estimating any functional
of the state by the expecta-
tion
as well. The recursive Bayesian filter, which
is also called the optimal filter, resolves exactly this problem in
two steps at each time
.
Suppose we know
. The prediction step is done ac-
cording to the following equation:
(3)
The observation
enables us to correct this prediction using
Bayes’s rule:
(4)
Under the specific assumptions of Gaussian noises
and
and linear functions and , these equations lead to the
Kalman filter’s equations. Unfortunately, this modeling is not
appropriate in many problems in signal and image processing,
which makes the calculation of the integrals in (3) and (4)
infeasible (no closed-form).
The original particle filter, which is called the bootstrap
filter [15], proposes to approximate the densities
by a finite
weighted sum of
Dirac densities centered on elements of
, which are called particles.
The application of the bootstrap filter requires that one knows
how to do the following:
• sample from initial prior marginal
;
• sample from
for all ;
• compute
for all through a known
function
such that , where
missing normalization must not depend on
.
The algorithm then evolves the particle set
, where is the particle and
its weight, such that the density can be approximated by
the density
. In the bootstrap filter, the
particles are “moved” by sampling from the dynamics (1), and
importance sampling theory shows that the weighting is only
based on likelihood evaluations. In the most general setting
[9], the displacement of particles is obtained by sampling
from an appropriate density
, which might depend on the
data as well. The complete procedure is summarized in Fig. 1.
Some convergence results of the empirical distributions to the
posterior distribution on the path space have been proved when
the number
of particles tends toward infinity [25], [10].
In the path space
, each particle at time can be
considered to be a discrete path of length
. Compared
with the particle filter presented in Fig. 1, particle filtering
in the space of paths consists of incrementing the particle
state space at each time step and representing each particle
by the concatenation of the new position at time
and the set
of previous positions between times 0 and
. In [25], the
fluctuations on path space of the so-called interacting particle
systems are studied. In the context of sequential Monte Carlo
methods [10] that cover most of the particle filtering methods
proposed in the last few years, the convergence and the rate of
convergence of order
of the average mean square error
is proved. Under more restrictive hypotheses, the almost-sure
convergence is proved as well [10].
To evaluate the degeneracy of the particle set, the effective
sample size has been defined [21], [23]. As advocated in [9], a
resampling step is performed inthe algorithm presented in Fig. 1
in an adaptive way when the effective sample size, estimated by
, is under a given threshold. It avoids to obtain a degenerate
HUE et al.: SEQUENTIAL MONTE CARLO METHODS FOR MTT AND DATA FUSION 311
Fig. 1. Basic particle filter with adaptive resampling.
particle set composed of only few particles with high weights
and all the others with very small ones.
Beside the discretization of the filtering integrals, the use of
such particles enables the maintenance of multiple hypotheses
on the position of the target and to keep in the long term only
the particles whose position is likely given the whole sequence
of observations.
We find more details on the algorithm in [9] or [15] and on
adaptive resampling in [9] and [21]. After these recalls, let us
present briefly the multitarget tracking problem and its clas-
sical solutions, as well as the existing works on particle filtering
methods for MTT. Then, we will propose the MTPF.
III. M
ULTITARGET PARTICLE FILTER
A. MTT Problem and Its Classical Treatment
Let
be the number of targets to track that are assumed to
be known and fixed for the moment (the case of a varying un-
known number will be addressed in Section III-C). The index
designates one among the targets and is always used as first
superscript. Multitarget tracking consists of estimating the state
vector made by concatenating the state vectors of all targets. It is
generally assumed that the targets aremoving according to inde-
pendent Markovian dynamics. At time
,
follows the state equation (1) decomposed in partial equa-
tions
(5)
The noises
and are supposed only to be white both
temporally and spatially and independent for
.
The observation vector collected at time
is denoted by
. The index is used as first superscript to refer
to one of the
measurements. The vector is composed of
detection measurements and clutter measurements. The false
alarms are assumed to be uniformly distributed in the obser-
vation area. Their number is assumed to arise from a Poisson
density of parameter
, where is the volume of the obser-
vation area, and
is the number of false alarms per unit volume.
As we do not know the origin of each measurement, one has to
introduce the vector
to describe the associations between the
measurements and the targets. Each component
is a random
variable that takes its values among
. Thus,
indicates that is associated with the th target. In this case,
is a realization of the stochastic process
if (6)
Again, the noises
and are supposed only to be
white noises, independent for
. We assume that the func-
tions
are such that they can be associated with functional
forms
such that
We dedicate the model 0 to false alarms. Thus, if ,
the
th measurement is associated with the clutter, but we do
not associate any kinematic model to false alarms.
As the indexing of the measurements is arbitrary, all the mea-
surements have the same a priori probability to be associated
with a given model
. At time , these association probabilities
define the vector
. Thus,
for
, for all is
the discrete probability that any measurement is associated with
the
th target.
To solve the data association, some assumptions are com-
monly made [3].
A1) One measurement can originate from one target or from
the clutter.
A2) One target can produce zero or one measurement at one
time.
The assumption A1) expresses that the association is exclusive
and exhaustive. Consequently,
.
Assumption A2) implies that
may differ from and,
above all, that the association variables
for
are dependent.
Under these assumptions, the MHT algorithm [28] builds re-
cursively the association hypotheses. One advantage of this al-
gorithm is that the appearance of a new target is hypothesized
at each time step. However, the complexity of the algorithm in-
creases exponentially with time. Some pruning solutions must
be found to eliminate some of the associations.
The JPDAF begins with a gating of the measurements. Only
the measurements that are inside an ellipsoid around the pre-
dicted state are kept. The gating assumes that the measurements
are distributed according to a Gaussian law centered on the pre-
dicted state. Then, the probabilities of each association
are estimated. As the variables are assumed dependent by
312 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 50, NO. 2, FEBRUARY 2002
A2), this computation implies the exhaustive enumeration of all
the possible associations
for .
The novelty in the PMHT algorithm [12], [32], [33] consists
of replacing the assumption A2) by A3):
A3) One target can produce zero or several measurements at
one time.
This assumption is often criticized because it does not match the
physical reality. However, from a mathematical point of view, it
ensures the stochastic independence of the variables
and it
drastically reduces the complexity of the
vector estimation.
The assumptions A1) and A3) will be kept in the MTPF pre-
sented later. Let us present now the existing works solving MTT
with particle filtering methods.
B. Related Work: MTT With Particle Filtering Methods
In the context of multitarget tracking, particle filtering
methods are appealing: As the association needs only to be
considered at a given time iteration, the complexity of data
association is reduced. First, two extensions of the bootstrap
filter have been considered. In [2], a bootstrap-type algorithm
is proposed in which the sample state space is a “(multitarget)
state space.” However, nothing is said about the association
problem that needs to be solved to evaluate the sample weights.
It is, in fact, the ability of the particle filtering to deal with
multimodality due to (high) clutter that is pointed out compared
with deterministic algorithms like the nearest neighbor filter
or the probabilistic data association (PDA) filter. No examples
with multiple targets are presented. The simulations only deal
with a single target in clutter with a linear observation model.
In [14], a hybrid bootstrap filter is presented where the particles
evolve in a single-object state space. Each particle gives a
hypothesis on the state of one object. Thus, the a posteriori
law of the targets, given the measurements, is represented by
a Gaussian mixture. Each mode of this law then corresponds
to one of the objects. However, as pointed out in [14], the
likelihood evaluation is possible only under the availability of
the “prior probabilities of all possible associations between”
the measurements and the targets. It may be why the simulation
example only deals with one single target in clutter. Even if
the likelihood could be evaluated, the way to represent the a
posteriori law by a mixture can lead to the loss of one of the
targets during occlusions. The particles tracking an occluded
target get very small weights and are therefore discarded during
the resampling step. This fact has been pointed out in [29].
In image analysis, the Condensation algorithm has been ex-
tended to the case of multiple objects as well. In [24], the case
of two objects is considered. The hidden state is the concate-
nation of the two single-object states and of a binary variable
indicating which object is closer to the camera. This latter vari-
able solves the association during occlusion because the mea-
surements are affected to the foreground object. Moreover, a
probabilistic exclusion principle is integrated to the likelihood
measurement to penalize the hypotheses with the two objects
overlapping. In [6], the state is composed of an integer equal to
the number of objects and of a concatenation of the individual
states. A three-dimensional (3-D) representation of the objects
gives access to their depth ordering, thus solving the association
issue during occlusions. Finally, in mobile robotics [29], a par-
ticle filter is used for each object tracked. The likelihood of the
measurements is written like in a JPDAF. Thus, the assignment
probabilities are evaluated according to the probabilities of each
possible association. Given these assignment probabilities, the
particle weights can be evaluated. The particle filters are then
dependent through the evaluation of the assignment probabili-
ties. Independently of the two latter works [6] and [29], we have
developed the MTPF, where the data association is approached
in the same probabilistic spirit as the basic PMHT [12], [32].
First, to estimate the density
,
with particle filtering methods, we must
choose the state space for the particles. As mentioned before,
a unique particle filter with a single-target state space seemed
to us a poor choice as the particles tracking an occluded object
would be quickly discarded. We have considered using one
particle filter per object but without finding a consistent way
to make them dependent. The stochastic association vector
introduced in Section III-A could also be considered to be an
additional particle component. However, as the ordering of the
measurements is arbitrary, it would not be possible to devise a
dynamic prior on it. Moreover, the state space would increase,
further making the particle filter less effective. Finally, we have
chosen to use particles whose dimension is the sum of those
of the individual state spaces corresponding to each target, as
in [6] and [24]. Each of these concatenated vectors then gives
jointly a representation of all targets. Let us describe the MTPF.
Further details on the motivations for the different ingredients
of the MTPF can be found in [18].
C. MTPF Algorithm
Before describing the algorithm itself, let us first notice that
the association probability
that a measurement is associated
with the clutter is a constant that can be computed
(7)
(8)
(9)
where
is the number of measurements arising from the
clutter at time
. Assuming that there are clutter originated
measurements among the
measurements collected at
time
, the a priori probability that any measurement comes
from the clutter is equal to
; hence, we get the equality
used to derive (9) from (8).
The initial set
of particles is
such that each component
for is sampled
from
independently from the others. Assume we have
obtained
with .
Each particle is a vector of dimension
, where we de-
note by
the th component of and where designates
the dimension of target
.
HUE et al.: SEQUENTIAL MONTE CARLO METHODS FOR MTT AND DATA FUSION 313
The prediction is performed by sampling from some proposal
density
. In the bootstrap filter case, coincides with the dy-
namics (5)
.
.
.
(10)
Examine now the computation of the likelihood of the obser-
vations conditioned on the
th particle. We can write for all
(11)
Note that the first equality in (11) is true only under the
assumption of conditional independence of the measure-
ments, which we will make. To derive the second equality
in (11), we have used the total probability theorem with
the events
and under the supplementary
assumption that the normalization factors between
and
is the same for all .
We still need to estimate at each time step the association
probabilities
, which can be seen as the stochastic
coefficients of the
-component mixture. Two main ways have
been found in the literature to estimate the parameters of this
model: the expectation maximization (EM) method (and its sto-
chastic version the SEM algorithm [5]) and the data augmenta-
tion method. The second one amounts in fact to a Gibbs sampler.
In [12], [32], and [33], the EM algorithm is extended and ap-
plied to multitarget tracking. This method implies that the vec-
tors
and are considered as parameters to estimate. The
maximization step can be easily conducted in the case of de-
terministic trajectories, and the additional use of a maximum
a posteriori (MAP) estimate enables the achievement of it for
nondeterministic trajectories. Yet, the nonlinearity of the state
and observation functions makes this step very difficult. Finally,
the estimation is done iteratively in a batch approach we would
like to avoid. For these reasons, we have not chosen an EM al-
gorithm to estimate the association probabilities.
The data augmentation algorithm is quite different in its prin-
ciple. The vectors
, , and are considered to be random
variables with prior densities. Samples are then obtained iter-
atively from their joint posterior using a proper Markov chain
Monte Carlo (MCMC) technique, namely, the Gibbs sampler.
This method has been studied in [4], [8], [13], [30], and [31],
for instance. It can be run sequentially at each time period. The
Gibbs sampler is a special case of the Metropolis–Hasting algo-
rithm with the proposal densities being the conditional distribu-
tions and the acceptance probability being consequently always
equal to one. See [7] for an introduction to MCMC simulation
methods and also for a presentation of the EM algorithm.
Let
be the stochastic variable associated with the proba-
bility
.For , the method consists of gener-
ating a Markov chain that converges to the stationary distribu-
tion
, which cannot be sampled directly. For that, we
must get a partition
of and to sample alternatively
from the conditional posterior distribution of each component of
the partition. Let the index
denote the iterations in the Gibbs
sampler. The second subscript of the vectors refers to the itera-
tion counter. Assume that the
first elements of the Markov
chain
have been drawn. We sample the com-
ponents of
as follows:
.
.
.
.
.
.
In our case, at a given instant , we follow this algorithm with
for
for
for
(12)
The initialization of the Gibbs sampler consists of assigning uni-
form association probabilities, i.e.,
for all
, and taking , i.e., the cen-
troid of the predicted particle set. The
variables do not need
initializing because at the first time step of the Gibbs sampler,
they will be sampled conditioned on
and
. Then, suppose that at instant , we have already simulated
. The th iteration is handled as follows.
• As the
are assumed to be independent,
their individual discrete conditional density reads
(13)
Assignment variables
are discrete, and we can write
1
if
if
(14)
The realizations
of the vector are then sam-
pled according to the weights
for
• Mixture proportion vector is drawn from the con-
ditional density
Dirichlet
(15)
where we denote by
the number of equal to
and where
Dirichlet
denotes the Dirichlet
distribution on the simplex
, ,
1
Using Bayes’s rule,
p
(
a
j
b; c
)=(
p
(
b
j
a; c
)
p
(
a
j
c
))
=p
(
b
j
c
)
.
