# Track Coalescence and Repulsion:
MHT, JPDA, and BP

Thomas Kropfreiter1, Florian Meyer1, Stefano Coraluppi2, Craig Carthel2, Rico
Mendrzik3, and Peter Willett4
1Scripps Institution of Oceanography and Department of Electrical and Computer
Engineering  University of California San Diego, La Jolla, CA, USA ({flmeyer,
<EMAIL_ADDRESS> 2Systems & Technology Research, Woburn, MA, USA
({stefano.coraluppi<EMAIL_ADDRESS> 3Ibeo Automotive Systems
GmbH, Hamburg, Germany<EMAIL_ADDRESS> 4Department of ECE,
University of Connecticut, Storrs, CT, USA<EMAIL_ADDRESS>
###### Abstract

Joint probabilistic data association (JPDA) and multiple hypothesis tracking
(MHT) introduced in the 70s, are still widely used methods for multitarget
tracking (MTT). Extensive studies over the last few decades have revealed
undesirable behavior of JPDA and MHT type methods in tracking scenarios with
targets in close proximity. In particular, JPDA suffers from the track
coalescence effect, i.e., estimated tracks of targets in close proximity tend
to merge and can become indistinguishable. Interestingly, in MHT, an opposite
effect to track coalescence called track repulsion can be observed. In this
paper, we review the estimation strategies of the MHT, JPDA, and the recently
introduced belief propagation (BP) framework for MTT and we investigate if BP
also suffers from these two detrimental effects. Our numerical results
indicate that BP-based MTT can mostly avoid both track repulsion and
coalescence.

###### Index Terms:

Multitarget tracking, joint probabilistic data association, multiple
hypothesis tracking, belief propagation

##  I Introduction

Multitarget tracking (MTT) [1, 2, 3, 4, 5] is a key signal processing task
that enables situational awareness in a variety of applications. Historically,
the most important MTT scenarios are aerospace surveillance [6, 7, 8] and
applied ocean sciences [9, 10, 11]. However, in the context of autonomous
navigation, robotic applications are recently drawing increased attention [12,
13]. MTT methods are typically derived in the Bayesian estimation framework
and they rely on different modeling paradigms, which can be broadly classified
as vector-type [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24] and set-type [25,
26, 27, 28, 29].

Vector-type MTT methods describe both multitarget states and measurements by
random vectors. The well-known joint probabilistic data association (JPDA)
filter [7] introduced in the 70s established the probabilistic data
association (DA) paradigm for MTT. The JPDA aims to calculate the minimum mean
square error (MMSE) state estimates for a known number of targets. Data
association hypotheses are employed in a soft manner by performing weighted
sums over all possible associations at each time step. Performing these
weighted sums is equivalent to calculate the marginal posterior probability
density functions needed for state estimation by “marginalizing out” the
unknown data association vector. In contrast, the very popular multiple
hypothesis tracking (MHT) methods [6] seek to calculate the maximum a
posteriori (MAP) estimate of an entire sequence of data association vectors.
In particular, MHT methods perform hard data association based on the single
most likely measurement-target combination over multiple consecutive instances
of time [6]. After an approximate MAP estimate of the data association vectors
has been fixed, MAP or MMSE estimates of the target states can be calculated
individually by means of sequential Bayesian estimation, e.g., extended Kalman
filtering.

Extensive studies over the last few decades have revealed undesirable behavior
of JPDA and MHT type methods in tracking scenarios with targets in close
proximity. In particular, JPDA type methods suffer from the track coalescence
effect [30], i.e., estimated tracks of closely spaced targets that move in
parallel will tend to come together, merge, and can become indistinguishable.
Interestingly, in MHT, an opposite effect to track coalescence called track
repulsion [31] can be observed. In particular, estimated tracks of close
targets that move in parallel repel each other in the sense that their
separation is larger compared to the distance between the actual targets.
Methodological adaptations tailored to alleviate these detrimental effects
have been introduced in [32, 33, 34, 35]. Track coalescence can be avoided by
using a special hypothesis pruning strategy [32, 33] and track repulsion by
selecting a global hypothesis that is not the single MAP solution, but an
element of a MAP equivalence class [34]. While these adaptations can reduce
the effects of track coalescence and repulsion, they are unsuitable for large-
scale tracking scenarios with an unknown and time-varying number of targets.

Set-type methods describe both multitarget states and measurements by random
finite sets. Some approaches in this class including probability hypothesis
density (PHD) methods [25, 26, 27, 29] avoid data association at the cost of a
reduced target detection and tracking performance as well as a suboptimal
postprocessing step for track formation. Some other RFS methods rely on
association strategies similar to the ones performed by the JPDA filter and
thus also suffer from track coalescence [36, 28].

Recently developed graph signal processing methods for MTT [37, 38, 5, 13] fit
in the vector-type family. Here, filtering and data association for randomly
appearing and disappearing targets are described by a joint factor graph. This
factor graph provides the blueprint for sum-product message passing algorithms
also known as belief propagation (BP). By following the MMSE estimation
approach of the JPDA, BP performs efficient marginalization operations and can
provide scalable solutions to two pertinent problems: the data association
problem that arises due to measurement-origin uncertainty and the presence of
clutter as well as the nonlinear filtering problem with appearing and
disappearing target states.

In this paper, we review the track coalescence and repulsion effects of the
JPDA and MHT methods and investigate them in the context of BP-based MTT. The
contributions of this paper are as follows.

  * •

We review the estimation strategies of the MHT, JPDA, and BP approaches to
MTT.

  * •

We discuss the track coalescence and track repulsion effects of MHT, JPDA, and
BP.

  * •

We demonstrate significantly reduced track coalescence and repulsion of BP in
three different case studies.

_Notation:_ Random variables are displayed in sans serif, upright fonts; their
realizations in serif, italic fonts. Vectors and matrices are denoted by bold
lowercase and uppercase letters, respectively. For example, a random variable
and its realization are denoted by $\mathsf{x}$ and $x$ and a random vector
and its realization by $\bm{\mathsf{x}}$ and $\bm{x}$. Furthermore,
${\bm{x}}^{\text{T}}$ denotes the transpose of vector $\bm{x}$; $\propto$
indicates equality up to a normalization factor; $f(\bm{x})$ denotes the
probability density function (pdf) of random vector $\bm{\mathsf{x}}$.
Finally, $1(a)$ denotes the unit sample function (i.e., $1(a)=1$ if $a=0$ and
$0$ otherwise).

##  II General System Model

In this section, we introduce a general MTT system model.

###  II-A Potential Target States and State-Transition PDF

In most practical applications, the number of targets is time-varying and
unknown. Here, we introduce potential target (PT) states [38, 5]. At discrete
time $k\geq 0$, the number of PTs $j_{k}$ is the maximum possible number of
targets that have generated a measurement. The state of PT
$j\in\\{1,\dots,j_{k}\\}$ consists of a binary existence variable
$\mathsf{r}^{(j)}_{k}\in\\{0,1\\}$ and a kinematic state
$\bm{\mathsf{x}}^{(j)}_{k}\in\mathbb{R}^{n}$. The existence variable
$\mathsf{r}^{(j)}_{k}$ models the existence/nonexistence of PT $j$ in the
sense that PT $j$ exists at time $k$ if and only if
$\mathsf{r}^{(j)}_{k}\\!=\\!1$. The state $\bm{\mathsf{x}}^{(j)}_{k}$ of PT
$j$ consists of the PT’s position and possibly further motion-related
parameters.

The kinematic state $\bm{\mathsf{x}}^{(j)}_{k}$ of PTs with
$\mathsf{r}^{(j)}_{k}\\!=\\!0$ is not defined and can be represented by an
arbitrary “dummy” pdf $f_{\text{D}}\big{(}\bm{x}^{(j)}_{k}\big{)}$ with $\int
f_{\text{D}}\big{(}\bm{x}^{(j)}_{k}\big{)}\mathrm{d}\bm{x}^{(j)}_{k}=1$. All
pdfs of PT states, $f\big{(}\bm{x}^{(j)}_{k},r^{(j)}_{k}\big{)}$, have the
property that
$f\big{(}\bm{x}^{(j)}_{k},0\big{)}=f^{(j)}_{k}f_{\text{D}}\big{(}\bm{x}^{(j)}_{k}\big{)}$,
where $f^{(j)}_{k}\\!\in\\![0,1]$ is the probability of nonexistence. For each
PT state $\big{(}\bm{\mathsf{x}}^{(j)}_{k-1},\;\mathsf{r}^{(j)}_{k-1}\big{)}$,
$j\in\\{1,\dots,j_{k-1}\\}$ at time $k-1$, there is one “legacy” PT state
$\big{(}\underline{\bm{\mathsf{x}}}^{(j)}_{k},\;\underline{\mathsf{r}}^{(j)}_{k}\big{)}$
at time $k$. The motion and disappearance of each PT is modeled by the single-
target state-transition pdf
$f\big{(}\underline{\bm{x}}_{k}^{(j)},\underline{r}_{k}^{(j)}\big{|}\bm{x}_{k-1}^{(j)},r_{k-1}^{(j)}\big{)}$
that involves the kinematic state-transition pdf
$f\big{(}\underline{\bm{x}}_{k}^{(j)}\big{|}\bm{x}_{k-1}^{(j)}\big{)}$ and the
probability of target survival $p_{\text{s}}$ (e.g., see [5, Section VIII-C]).
It is assumed that each target evolves independently. At time $k=0$, the
target states are assumed statistically independent across PTs $j$. However,
often no prior information for PT states is available, i.e., $j_{0}=0$.

###  II-B New PTs, Data Association, and Measurement Likelihood Function

At each time (or “scan”) $k\geq 1$, a sensor produces measurements
$\bm{\mathsf{z}}^{(m)}_{k}$, $m\in\big{\\{}1,\dots,\mathsf{m}_{k}\big{\\}}$.
(Note that up to time $k$, the number of measurements $\mathsf{m}_{k}$ is
random.) Each measurement can have one of three possible sources: (i) a new PT
representing a target that generates a measurement for the first time; (ii) a
legacy PT representing a target that has generated at least one measurement
before; and (iii) clutter. At time $k$, $\mathsf{m}_{k}$ new PT states are
introduced, i.e.,
$\big{(}\overline{\bm{\mathsf{x}}}^{(m)}_{k}\;\overline{\mathsf{r}}^{(m)}_{k}\big{)}$,
$m\in\\{1,\dots,\mathsf{m}_{k}\\}$. Here, $\overline{\mathsf{r}}^{(m)}_{k}=1$
means that measurement $\bm{\mathsf{z}}^{(m)}_{k}$ has source (i) and
$\overline{\mathsf{r}}^{(m)}_{k}=0$ means that measurement
$\bm{\mathsf{z}}^{(m)}_{k}$ has either source (ii) or (iii). The birth of new
targets is modeled by a Poisson point process with mean $\mu_{\text{b}}$ and
pdf $f_{\text{b}}\big{(}\overline{\bm{x}}_{k}\big{)}$. The joint vectors of
all existence variables and kinematic states at time $k$ are denoted as
$\bm{\mathsf{x}}_{k}=\big{[}\underline{\bm{\mathsf{x}}}^{(1)}_{k}\cdots\underline{\bm{\mathsf{x}}}^{(j_{k-1})}_{k}\hskip
0.85358pt\hskip
0.85358pt\overline{\bm{\mathsf{x}}}^{(1)}_{k}\cdots\overline{\bm{\mathsf{x}}}^{(m_{k})}_{k}\big{]}^{\mathrm{T}}$
and
$\mathsf{r}_{k}=\big{[}\underline{\mathsf{r}}^{(1)}_{k}\cdots\underline{\mathsf{r}}^{(j_{k-1})}_{k}\hskip
0.85358pt\hskip
0.85358pt\overline{\mathsf{r}}^{(1)}_{k}\cdots\overline{\mathsf{r}}^{(m_{k})}_{k}\big{]}^{\mathrm{T}}$,
respectively.

The target represented by a PT $j$ is detected (in the sense that it generates
a measurement $\bm{\mathsf{z}}_{k}^{(m)}$) with probability $p_{\text{d}}$.
The statistical relationship of a measurement $\bm{\mathsf{z}}_{k}^{(m)}$ and
the kinematic state $\bm{\mathsf{x}}_{k}^{(j)}$ of the corresponding detected
PT, is described by the conditional pdf
$f\big{(}\bm{z}^{(m)}_{k}|\bm{x}^{(j)}_{k}\big{)}$. This conditional pdf can
be derived from the measurement model of the sensor. Clutter measurements are
modeled by a Poisson point process with mean $\mu_{\text{c}}$ and pdf
$f_{\text{c}}\big{(}\bm{z}^{(m)}_{k}\big{)}$. At time $k$, the joint
measurement vector is fixed and denoted as
$\bm{z}_{k}\\!\triangleq\big{[}\bm{z}_{k}^{(1)\hskip
0.85358pt\text{T}}\cdots\hskip 0.85358pt\bm{z}_{k}^{(m_{k})\hskip
0.85358pt\text{T}}\big{]}^{\text{T}}.$ Thus, the total number of (legacy and
new) PTs is $j_{k}\\!=j_{k-1}\\!+m_{k}$.

In MTT, the source of measurements is uncertain. In particular, it is unknown
which measurement is generated by which PT, and a measurement may also be
clutter, i.e., not generated by any PT. The considered MTT problem follows the
conventional DA assumption, i.e., a target can originate at most one
measurement and a measurement can be originated from at most one target [1, 3,
5]. The association between $m_{k}$ measurements and $j_{k-1}$ legacy PTs at
time $k$ is typically modeled by a “target-oriented” DA vector
$\bm{\mathsf{a}}_{k}=\big{[}\mathsf{a}_{k}^{(1)}\cdots\hskip
0.85358pt\mathsf{a}_{k}^{(j_{k-1})}\big{]}^{\text{T}}$. The target-oriented
association variable $\mathsf{a}_{k}^{(j)}$ is $m\in\\{1,\dots,m_{k}\\}$ if PT
$j$ originated measurement $m$ and zero if PT $j$ did not originate any
measurement [1, 5]. A DA vector $\bm{a}_{k}$ is valid if it does not violate
the conventional DA association assumption.

We also introduce the “measurement-oriented” DA vector
$\bm{\mathsf{b}}_{k}=\big{[}\mathsf{b}_{k}^{(1)}\cdots\hskip
0.85358pt\mathsf{b}_{k}^{(m_{k})}\big{]}^{\text{T}}$. The measurement-oriented
association variable $\mathsf{b}_{k}^{(m)}$ is $j\in\\{1,\dots,j_{k-1}\\}$ if
measurement $m$ is generated by legacy PT $j$ and zero if it is generated by
clutter or by a newly detected PT. The DA representation by means of both
target-oriented and measurement-oriented DA vectors is redundant in the sense
that for any valid $\bm{a}_{k}$, there is exactly one valid $\bm{b}_{k}$ and
for any valid $\bm{b}_{k}$, there is exactly one valid $\bm{a}_{k}$. However,
this representation based on both $\bm{a}_{k}$ and $\bm{b}_{k}$ makes it
possible to develop scalable MTT methods that exploit the structure of the DA
problem to reduce computational complexity.

###  II-C Joint Posterior PDF and Factor Graph

Following common assumptions [1, 6, 3, 28, 37, 39, 38, 5], the joint posterior
pdf of $\bm{\mathsf{x}}_{0:k}$, $\bm{\mathsf{r}}_{0:k}$,
$\bm{\mathsf{a}}_{1:k}$, and $\bm{\mathsf{b}}_{1:k}$ conditioned on fixed
measurements $\bm{z}_{1:k}$ can be obtained as

$\displaystyle
f\big{(}\bm{x}_{0:k},\bm{r}_{0:k},\bm{a}_{1:k},\bm{b}_{1:k}\big{|}\bm{z}_{1:k}\big{)}$
$\displaystyle\propto\Bigg{(}\prod^{j_{0}}_{j^{\prime\prime}=1}f\Big{(}\bm{x}^{(j^{\prime\prime})}_{0},r^{(j^{\prime\prime})}_{0}\Big{)}\Bigg{)}\prod^{k}_{k^{\prime}=1}$
$\displaystyle\hskip 12.80373pt\times\hskip
2.84526pt\Bigg{(}\prod^{j_{k^{\prime}-1}}_{j^{\prime}=1}f\big{(}\underline{\bm{x}}^{(j^{\prime})}_{k^{\prime}},\underline{r}^{(j^{\prime})}_{k^{\prime}}\big{|}\bm{x}^{(j^{\prime})}_{k^{\prime}-1},r^{(j^{\prime})}_{k^{\prime}-1}\big{)}\Bigg{)}$
$\displaystyle\hskip
12.80373pt\times\Bigg{(}\prod^{j_{k^{\prime}-1}}_{j=1}q\big{(}\underline{\bm{x}}^{(j)}_{k^{\prime}}\\!,\underline{r}^{(j)}_{k^{\prime}}\\!,a^{(j)}_{k^{\prime}};\bm{z}_{k^{\prime}}\big{)}\prod^{m_{k^{\prime}}}_{m^{\prime}=1}\Psi_{j,m^{\prime}}\big{(}a_{k^{\prime}}^{(j)},b_{k^{\prime}}^{(m^{\prime})}\big{)}\Bigg{)}$
$\displaystyle\hskip 12.80373pt\times\hskip
2.84526pt\prod^{m_{k^{\prime}}}_{m=1}v\big{(}\overline{\bm{x}}^{(m)}_{k^{\prime}}\\!,\overline{r}^{(m)}_{k^{\prime}}\\!,b^{(m)}_{k^{\prime}};\bm{z}_{k^{\prime}}^{(m)}\big{)}.$
(1)

Here, the new PT factor
$v\big{(}\overline{\bm{x}}^{(m)}_{k}\\!,\overline{r}^{(m)}_{k},$
$b^{(m)}_{k};\bm{z}_{k}^{(m)}\big{)}$ is given by

$\displaystyle
v\big{(}\overline{\bm{x}}^{(m)}_{k}\\!,1,b^{(m)}_{k};\bm{z}_{k}^{(m)}\big{)}$
$\displaystyle\triangleq\begin{cases}0\hskip
0.85358pt,&b^{(m)}_{k}\in\\{1,\dots,j_{k-1}\\}\\\\[5.69054pt] {\hskip
0.85358pt\frac{p_{\mathrm{d}}\hskip 0.85358pt\mu_{\text{b}}\hskip
0.85358ptf_{\text{b}}\big{(}\overline{\bm{x}}^{(m)}_{k}\big{)}}{\mu_{\text{c}}f_{\text{c}}\big{(}\bm{z}_{k}^{(m)}\big{)}}}f\big{(}\bm{z}_{k}^{(m)}\big{|}\overline{\bm{x}}^{(m)}_{k}\big{)}\hskip
0.85358pt,&b^{(m)}_{k}\\!=0\end{cases}$

and
$v\big{(}\overline{\bm{x}}^{(m)}_{k}\\!,0,b^{(m)}_{k};\bm{z}_{k}^{(m)}\big{)}\triangleq
f_{\text{D}}\big{(}\overline{\bm{x}}^{(m)}_{k}\big{)}$. In addition, the
legacy PT factor
$q\big{(}\underline{\bm{x}}^{(j)}_{k}\\!,\underline{r}^{(j)}_{k},$
$a^{(j)}_{k};\bm{z}_{k}\big{)}$ reads

$\displaystyle
q\big{(}\underline{\bm{x}}^{(j)}_{k}\\!,1,a^{(j)}_{k};\bm{z}_{k}\big{)}$
$\displaystyle\triangleq\begin{cases}\frac{p_{\text{d}}}{\mu_{\text{c}}f_{\text{c}}\big{(}\bm{z}_{k}^{(m)}\big{)}}f\big{(}\bm{z}_{k}^{(m)}\big{|}\hskip
0.85358pt\underline{\bm{x}}_{k}^{(j)}\big{)},&a^{(j)}_{k}=m\in\\{1,\dots,m_{k}\\}\\\\[5.69054pt]
1-p_{\text{d}}\hskip 0.85358pt,&a^{(j)}_{k}=0\end{cases}$

and
$q\big{(}\underline{\bm{x}}^{(j)}_{k}\\!,0,a^{(j)}_{k};\bm{z}_{k}\big{)}\triangleq
1(a^{(j)}_{k})$.

Finally, consistency of a pair of target-oriented and measurement-oriented
variables $\big{(}a_{k}^{(j)},b_{k}^{(m)}\big{)}$ is checked by the binary
indicator function $\Psi_{j,m}\big{(}a_{k}^{(j)},b_{k}^{(m)}\big{)}$. In
particular, $\Psi_{j,m}\big{(}a_{k}^{(j)},b_{k}^{(m)}\big{)}$ is zero if
$a_{k}^{(j)}=m,b^{(m)}_{k}\neq j$ or $b^{(m)}_{k}=j,a_{k}^{(j)}\neq m$ and one
otherwise (see [37, 5] for details). A detailed derivation of the joint
posterior in (1) has been introduced in [5, Section VIII-G].

After marginalizing out $\bm{\mathsf{b}}_{1:k}$ in (1), the resulting
posterior distribution

$\displaystyle
f\big{(}\bm{x}_{0:k},\bm{r}_{0:k},\bm{a}_{1:k},\big{|}\bm{z}_{1:k}\big{)}=\sum_{\bm{b}_{1:k}}f\big{(}\bm{x}_{0:k},\bm{r}_{0:k},\bm{a}_{1:k},\bm{b}_{1:k}\big{|}\bm{z}_{1:k}\big{)}$

represents a system model that is identical to the one used by the MHT111In
the MHT literature, the two vectors $\bm{\mathsf{r}}_{0:k}$ and
$\bm{\mathsf{a}}_{1:k}$ are typically represented by an equivalent single
vector $\bm{\mathsf{q}}_{0:k}$ referred to as global hypothesis. and a more
general version of the one used by JPDA methods. Note that this
marginalization is trivial since for $\bm{a}_{1:k}$ there is exactly one
$\bm{b}_{1:k}$ for which (1) is nonzero.

##  III Multitarget Tracking Strategies

In this section, we will review the different estimation strategies of the
MHT, the JPDA, and the BP methods and discuss how they can cause track
repulsion and track coalescence.

###  III-A MHT for Sequence MAP Estimation

MHT tracking aims to compute the sequence MAP estimate of
$\bm{\mathsf{x}}_{0:k},\bm{\mathsf{r}}_{0:k},\bm{\mathsf{a}}_{1:k}$ given
$\bm{z}_{1:k}$. In particular, first the sequence MAP of
$\bm{\mathsf{r}}_{0:k},\bm{\mathsf{a}}_{1:k}$ given $\bm{z}_{1:k}$ is computed
as

$\displaystyle\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{)}=\operatorname*{arg\,max}_{\bm{r}_{0:k},\bm{a}_{1:k}}\int
f\big{(}\bm{x}_{0:k},\bm{r}_{0:k},\bm{a}_{1:k}\big{|}\bm{z}_{1:k}\big{)}\hskip
0.85358pt\text{d}\bm{x}_{0:k}$

This is equivalent to searching for the most likely target-measurement
association or global hypothesis [6]. The MHT exploits the fact that the
marginalization operation $\int\text{d}\bm{x}_{0:k}$ can be performed
sequentially and, assuming a linear-Gaussian state-transition and measurement
model, also in closed form. The optimal MAP solution however is a “multi-scan
solution” in the sense that the optimal estimate
$\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{)}$
at time $k$ can potentially be completely different than the optimal estimate
$\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k-1},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k-1}\big{)}$
at time $k-1$. Thus, for optimal MAP estimation it is necessary to keep track
of an exponentially growing number of global hypotheses
$\big{(}\bm{r}_{0:k},\bm{a}_{1:k}\big{)}$.

Given the most likely target-measurement association
$\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{)}$,
the MHT then computes a MAP estimate of $\bm{\mathsf{x}}_{0:k}$ as

$\displaystyle\hat{\bm{x}}_{0:k}^{\text{MAP}}$ $\displaystyle=\hskip
0.85358pt\operatorname*{arg\,max}_{\bm{x}_{0:k}}f\big{(}\bm{x}_{0:k}\big{|}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k},\bm{z}_{1:k}\big{)}$
$\displaystyle=\hskip
0.85358pt\operatorname*{arg\,max}_{\bm{x}_{0:k}}f\big{(}\bm{x}_{0:k},\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{|}\bm{z}_{1:k}\big{)}\hskip
0.85358pt$ (4)

Note that for
$\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{)}$
fixed, the joint posterior (LABEL:eq:jointPosteriorComplete1) can be split
into independent posterior distributions (one for each detected target
according to $\hat{\bm{r}}^{\mathrm{MAP}}_{0:k}$, cf. (1)). The joint target
state estimate $\hat{\bm{x}}_{0:k}^{\text{MAP}}$ is then obtained by
performing Bayesian smoothing for each detected target in parallel and
independent of all other targets, e.g., by means of the Kalman smoother.

In a practical implementation, MHT computes the MAP estimates of the target-
measurement association in (LABEL:eq:MHT0) and of the target state in (4) only
over a sliding window of $N$ consecutive times steps in order to reduce the
computational burden. A hard decision is then made only on the association at
the oldest time step. MHT can be formulated in form of hypothesis-oriented MHT
[6] and track-oriented MHT [21, 22, 23]. The computational complexity of the
original hypothesis-oriented MHT formulation is still problematic due to the
high number of hypotheses. However, it can be reduced by discarding unlikely
hypotheses using, e.g., an efficient $m$-best assignment algorithm [19, 20].

The more efficient track-oriented MHT methods [21, 22, 23] represent
association hypotheses via a series of tree structures. Here, each tree
represents the possible association histories of a single target. The most
likely hypothesis is then found by choosing a leaf node for each single-target
tree in such a way that no measurement is used by more than one target. Fast
hypothesis-search is enabled by means of combinatorial optimization methods
[40, 24, 41, 42]. Track-oriented MHT offers a more compact representation of
the data association problem than hypothesis-oriented MHT. The two
formulations are equivalent under the assumption that target births and
clutter follow a Poisson point process [43].

A limitation of MHT tracking is known as track repulsion, which can arise when
targets come in close proximity. It is a direct consequence of performing hard
data association by considering only a single global hypothesis
$\big{(}\hat{\bm{r}}^{\mathrm{MAP}}_{0:k},\hat{\bm{a}}^{\mathrm{MAP}}_{1:k}\big{)}$
for target state estimation according to (4). Estimated tracks are displaced
at a greater distance than the targets themselves and thus result in an
increased estimation error. This displacement can also lead to track swaps (or
identity switches) [44, 45]. In a crossing target scenario the corresponding
estimated tracks have the tendency to “bounce” more often than they cross. A
method for mitigating the detrimental effects of track repulsion has been
presented in [34].

Figure 1: Factor graph for an unknown, time-varying number of objects,
corresponding to the joint posterior pdf
$f\big{(}\bm{x}_{0:k},\bm{r}_{0:k}\bm{a}_{1:k},\bm{b}_{1:k}\big{|}\bm{z}_{1:k}\big{)}$
in (1). The factor graph depicts one time step $k$. The time index $k$ is
omitted, and the following short notations are used: $n_{\mathrm{m}}\triangleq
m_{k}$, $n_{\mathrm{p}}\triangleq j_{k-1}$,
$\underline{\bm{y}}^{j}\triangleq[\underline{\bm{x}}^{(j)}_{k},\underline{r}^{(j)}_{k}]^{\text{T}}$,
$\overline{\bm{y}}^{m}\triangleq[\overline{\bm{x}}^{(m)}_{k},\overline{r}^{(m)}_{k}]^{\text{T}}$,
$\bm{a}^{j}\triangleq\bm{a}^{(j)}_{k}$,
$\bm{b}^{m}\triangleq\bm{b}^{(m)}_{k}$, $f^{j}\triangleq
f\big{(}\underline{\bm{x}}^{(j)}_{k},\underline{r}^{(j)}_{k}\big{|}\bm{x}^{(j)}_{k-1},r^{(j)}_{k-1}\big{)}$,
$q^{j}\triangleq
q\big{(}\underline{\bm{x}}^{(j)}_{k},\underline{r}^{(j)}_{k},a_{k}^{(j)};\bm{z}_{k}\big{)}$,
$v^{m}\triangleq
v\big{(}\overline{\bm{x}}^{(m)}_{k},\overline{r}^{(m)}_{k},b_{k}^{(m)};\bm{z}^{(m)}_{k}\big{)}$,
$\Psi^{j,m}\triangleq\Psi_{j,m}\big{(}a^{(j)}_{k},b^{(m)}_{k}\big{)}$.

###  III-B JPDA for Marginal MMSE-Estimation

JPDA methods are based on the MMSE estimation criterion. (An exception is the
set JPDA (SJPDA) filter [35], which aims to minimize the mean optimal sub-
pattern assignment metric.) In particular, the JPDA aims to calculate MMSE
estimates for each individual target state $\bm{\mathsf{x}}_{k}^{(j)}$, i.e.,

$\displaystyle\hat{\bm{x}}_{k}^{(j)\text{MMSE}}=\hskip
0.85358pt\int\bm{x}^{(j)}_{k}\hskip
0.85358ptf(\bm{x}^{(j)}_{k}|\bm{z}_{1:k})\hskip 0.85358pt\hskip
0.85358pt\text{d}\bm{x}_{k},\quad j=1,\dots,l_{k}\hskip 0.85358pt.$ (5)

Contrary to MHT, in conventional JPDA it is assumed that the true sequence of
existence variables $\hat{\bm{r}}_{0:k}$ (and thus the number of targets
$l_{k}$ at time $k$) is known, i.e.,

$f(\bm{x}_{0:k},\bm{a}_{1:k}|\bm{z}_{1:k})\propto
f(\bm{x}_{0:k},\hat{\bm{r}}_{0:k},\bm{a}_{1:k}|\bm{z}_{1:k}).\vspace{0mm}$ (6)

The marginals needed for optimum state estimation in (5) are then calculated
from (6) according to

$\displaystyle f(\bm{x}^{(j)}_{k}|\bm{z}_{1:k})=\hskip
0.85358pt\int\sum_{\bm{a}_{1:k}}\hskip
0.85358ptf(\bm{x}_{0:k},\bm{a}_{1:k}|\bm{z}_{1:k})\hskip
0.85358pt\text{d}\bm{x}^{\sim(j)}_{0:k}\hskip 0.85358pt,$ (7)

where $\bm{x}^{\sim(j)}_{0:k}$ denotes the vector that consists of all entries
in $\bm{x}_{0:k}$ except $\bm{x}^{(j)}_{k}$. In practical implementations a
track initialization and termination heuristic provides an approximate
sequence of existence variables $\tilde{\bm{r}}_{0:k}$. Under the assumption
of linear and Gaussian measurement and motion models as used in the original
formulation of the JPDA, approximate closed-form expressions for the
marginalization (7) can be found. However, the summation over all joint
associations hypotheses performed in (7) remains computationally infeasible in
tracking scenarios with more than 7 targets. Further complexity reduction can
be obtained by applying heuristic approximations [46, 16], by using efficient
hypothesis management [17], or by limiting the number of target-measurement
associations [18].

A deficiency of JPDA tracking is known as track coalescence, which can arise
when targets come close to each other. In such situations, JPDA is unable to
maintain awareness of target positions and the estimated tracks tend to merge,
and can become indistinguishable. Track coalescence is a direct consequence of
calculating marginal distributions $f(\bm{x}^{(j)}_{k}|\bm{z}_{1:k})$ for the
individual target states according to (7). In particular, performing soft data
association by means of the sum over all possible associations
$\sum_{\bm{a}_{k}}$ for each time step $k$, has the effect that targets with
identical prior distribution will also have identical posterior distribution.
Consequently, as targets are in close proximity and their prior distributions
become similar, their posterior distributions tend to merge. Mitigation
methods include the pruning of association hypotheses [32, 33], approximate
calculation of marginal pdfs via information-theoretic based optimization
techniques [35, 47] or mean-field approximations [48].

###  III-C BP for Marginal Detection and MMSE-Estimation

In BP-based multitarget tracking [38, 5, 49, 13], we aim to compute the
marginal posterior pdfs
$f\big{(}\bm{x}^{(j)}_{k},\bm{r}^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}$ for each
PT $j\in\\{1,\dots,j_{k}\\}$. From these marginal posterior pdfs, detection
and MMSE estimation can be performed individually for each PT
$j\in\\{1,\dots,j_{k}\\}$.

In particular, for each PT, marginal posterior probability mass functions
$p\big{(}r^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}$ of the existence indicators
$\mathsf{r}^{(j)}_{k}$ are obtained from marginal posterior pdfs
$f\big{(}\bm{x}^{(j)}_{k},\bm{r}^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}$ as

$p\big{(}r^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}=\int
f\big{(}\bm{x}^{(j)}_{k},r^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}\hskip
0.85358pt\mathrm{d}\bm{x}^{(j)}_{k}.$ (8)

Targets are then declared to exist if
$p\big{(}r^{(j)}_{k}\\!=\\!1\big{|}\bm{z}_{1:k})$ is greater than a predefined
threshold $P_{\text{th}}$. For all targets declared to exist, state estimation
is performed by the MMSE estimator according to

$\hat{\bm{x}}^{(j)\text{MMSE}}_{k}\hskip 0.85358pt\triangleq\hskip
0.85358pt\int\bm{x}_{k}\hskip
0.85358pt\frac{f\big{(}\bm{x}^{(j)}_{k},r^{(j)}_{k}\\!=\\!1\big{|}\bm{z}_{1:k}\big{)}}{p\big{(}r^{(j)}_{k}\\!=\\!1\big{|}\bm{z}_{1:k}\big{)}}\hskip
0.85358pt\mathrm{d}\bm{x}_{k}.\vspace{0mm}$ (9)

Let us first investigate the calculation of the joint marginal pdfs
$f\big{(}\bm{x}_{k},r_{k}$ $\big{|}\bm{z}_{1:k}\big{)}$ that involves the
joint state and existence variable of all targets at time $k$. A naive
marginalization of (1) can be performed as

$\displaystyle f\big{(}\bm{x}_{k},\bm{r}_{k}\big{|}\bm{z}_{1:k}\big{)}$
$\displaystyle=\int\sum_{\bm{r}_{0:k-1}}\sum_{\bm{a}_{1:k}}\hskip
0.85358pt\hskip 0.85358pt\sum_{\bm{b}_{1:k}}\hskip 0.85358pt\hskip
0.85358ptf\big{(}\bm{x}_{0:k},\bm{r}_{0:k},\bm{a}_{1:k},\bm{b}_{1:k}\big{|}\bm{z}_{1:k}\big{)}\hskip
0.85358pt\mathrm{d}\bm{x}_{0:k-1}.$

By using (1) in (LABEL:eq:marginalization1) and marginalizing out the vectors
$\bm{x}_{1:k-2}$ $\bm{r}_{1:k-2}$, $\bm{a}_{1:k-1}$, and $\bm{b}_{1:k-1}$, we
obtain

$\displaystyle f\big{(}\bm{x}_{k},r_{k}\big{|}\bm{z}_{1:k}\big{)}$
$\displaystyle\hskip 17.07164pt=\int\sum_{\bm{r}_{k-1}}\sum_{\bm{a}_{k}}\hskip
0.85358pt\hskip 0.85358pt\sum_{\bm{b}_{k}}\hskip
0.85358ptf\big{(}\bm{x}_{k-1},\bm{r}_{k-1}\big{|}\bm{z}_{1:k-1}\big{)}$
$\displaystyle\hskip 34.1433pt\times
g(\bm{x}_{k},\bm{r}_{k},\bm{x}_{k-1},\bm{r}_{k-1},\bm{a}_{k},\bm{b}_{k})\hskip
0.85358pt\hskip 0.85358pt\mathrm{d}\bm{x}_{k-1}$ (11)

where
$g(\bm{x}_{k},\bm{r}_{k},\bm{x}_{k-1},\bm{r}_{k-1},\bm{a}_{k},\bm{b}_{k})$ is
equal to the last three lines in (1) (with $k^{\prime}$ replaced by $k$). Now
marginals for individual targets used in (8) and (9) can be obtained according
to
$f\big{(}\bm{x}^{(j)}_{k},r^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}=\int\sum_{\bm{r}^{\sim(j)}_{k}}f\big{(}\bm{x}_{k},r_{k}\big{|}\bm{z}_{1:k}\big{)}\mathrm{d}\bm{x}^{\sim(j)}_{k}$,
where $\bm{r}^{\sim(j)}_{k}$ denotes the vector that consists of all entries
in $\bm{r}_{k}$ except $\bm{r}^{(j)}_{k}$ and $\bm{x}^{\sim(j)}_{k}$ denotes
the vector that consists of all entries in $\bm{x}_{k}$ except
$\bm{x}^{(j)}_{k}$.

There are two interesting observations that can be made. (i) As in
conventional filtering methods, by exploiting the temporal factorization
structure of (1), the computational complexity can be strongly reduced. (Note
that, e.g., the summation over $\sum_{\bm{r}_{0:k-1}}$ in
(LABEL:eq:marginalization1) consists of $2^{\hskip
0.28453pt\prod^{k-1}_{k^{\prime}=0}j_{k^{\prime}}}$ terms, while the $k$
individual summations performed sequentially according to (11) only consist of
$\sum^{k-1}_{k^{\prime}=0}\hskip 0.85358pt2^{\hskip 0.28453ptj_{k^{\prime}}}$
terms.) (ii) The optimal MMSE solution is a “single scan solution” in the
sense that the optimal marginal pdf
$f\big{(}\bm{x}_{k},r_{k}\big{|}\bm{z}_{1:k}\big{)}$ at time $k$ can be
directly obtained from the optimal marginal pdf
$f\big{(}\bm{x}_{k-1},r_{k-1}\big{|}\bm{z}_{1:k-1}\big{)}$ at time $k-1$,
i.e., $f\big{(}\bm{x}_{k-1},r_{k-1}\big{|}\bm{z}_{1:k-1}\big{)}$ provides all
relevant information from past time steps.

BP makes it possible to also exploit the spatial factorization structure of
(1) to further reduce computational complexity [38, 5]. It operates on a
factor graph that represents the statistical model of the considered
estimation problem [50]. In particular, marginal pdfs are calculated
efficiently by performing local operations for the graph nodes and exchanging
their results called “messages” along the graph edges. Since the factor graph
in Fig. 1 representing this structure has cycles, BP can only provide
approximate marginal pdfs
$f\big{(}\bm{x}^{(j)}_{k},\bm{r}^{(j)}_{k}\big{|}\bm{z}_{1:k}\big{)}$,
$j\in\\{1,\dots,j_{k}\\}$ [50]. These approximate marginal pdfs are typically
referred to as “beliefs”. Furthermore, the cycles in the factor graph lead to
overconfident messages and beliefs [51], i.e., their spread “downplays” their
true uncertainty.

Since BP performs marginalization similarly to JPDA, it is expected to suffer
from track coalescence. However, as will be seen in Section IV, track
coalescence is reduced significantly compared to JPDA. This can be explained
by the particle-based implementation used by BP-based MTT that makes it
possible to more accurately represent multimodal marginal pdfs. Furthermore,
the overconfident nature of BP-based data association “upsells” the most
likely measurement-to-target associations hypothesis and thus leads to a
hybrid form of soft and hard data association and further reduces coalescence.

(a)

(b)

(c)

Figure 2: True target tracks used in scenario 1 (a), scenario 2 (b), and
scenario 3 (c). Initial target locations are indicated by a black circle. Note
that the two coordinate axis are scaled differently.

(a)

(b)

Figure 3: Scenario 1: D-Center (a) and D-Tracks (b) for the three considered
multitarget tracking methods.

##  IV Results

We simulate three different scenarios with two targets whose states consist of
two-dimensional (2D) position and velocity, i.e.,
$\bm{\mathsf{x}}_{k}^{(j)}=[\mathsf{x}_{1,k}^{(j)}\;\hskip
0.85358pt\mathsf{x}_{2,k}^{(j)}\;\hskip
0.85358pt\dot{\mathsf{x}}_{1,k}^{(j)}\;\hskip
0.85358pt\dot{\mathsf{x}}_{2,k}^{(j)}]^{\text{T}}$, $j\in\\{1,2\\}$. The
targets are assumed to move according to the near constant-velocity motion
model, i.e., $\bm{\mathsf{x}}_{k}^{(j)}=\bm{A}\hskip
0.85358pt\bm{\mathsf{x}}_{k-1}^{(j)}+\bm{\mathsf{u}}_{k}^{(j)}$ where
$\bm{\mathsf{u}}_{k}^{(j)}\sim\mathcal{N}(\bm{0},\bm{\Sigma}_{\bm{u}})$ is a
sequence of independent and identically distributed (iid) 2D Gaussian random
vectors. The matrixes $\bm{A}$ and $\bm{\Sigma}_{\bm{u}}$ represent
discretized continuous-time kinematics [52], i.e.,

$\small\bm{A}={\begin{pmatrix}1&0&T&0\\\ 0&1&0&T\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix}}\hskip
0.85358pt,\quad\,\bm{\Sigma}_{\bm{u}}={\begin{pmatrix}\frac{T^{3}}{3}&0&\frac{T^{2}}{2}&0\\\
0&\frac{T^{3}}{3}&0&\frac{T^{2}}{2}\\\ \frac{T^{2}}{2}&0&T&0\\\
0&\frac{T^{2}}{2}&0&T\end{pmatrix}}\sigma^{2}_{\bm{u}}\vspace{-1mm}$

with $\sigma^{2}_{\bm{u}}\\!=0.1\hskip 0.85358pt\hskip
0.85358pt\text{m}^{2}/\text{s}^{4}$ and $T=1\hskip 0.85358pt\text{s}$. The
true target tracks for the three different scenarios as well as the initial
target positions are shown in Fig. 2. Note that in scenarios 1 and 2, the
minimum distance between the two targets is equal to $\sigma_{\bm{v}}$
(defined in what follows). The region-of-interest (ROI) is given by
$[-750\hskip 0.85358pt\text{m},\hskip 0.85358pt750\hskip
0.85358pt\text{m}]\times[-750\hskip 0.85358pt\text{m},\hskip
0.85358pt750\hskip 0.85358pt\text{m}]$ for scenarios 1 and 3 and by $[0\hskip
0.85358pt\text{m},\hskip 0.85358pt1500\hskip
0.85358pt\text{m}]\times[-750\hskip 0.85358pt\text{m},\hskip
0.85358pt750\hskip 0.85358pt\text{m}]$ for scenario 2.

A single sensor performs noisy measurements of the positions of the targets.
More specifically, the target-generated measurements are given by

$\bm{\mathsf{z}}^{(m)}_{k}=\hskip 0.85358pt[\mathsf{x}_{1,k}^{(m)}\;\hskip
0.85358pt\mathsf{x}_{2,k}^{(m)}]^{\text{T}}\\!+\hskip
0.85358pt\bm{\mathsf{v}}_{k}^{(m)}\vspace{-.8mm}$

where
$\bm{\mathsf{v}}_{k}^{(m)}\sim\mathcal{N}(\bm{0},\sigma_{\bm{v}}^{2}\bm{I}_{2})$
with $\sigma_{\bm{v}}=10$m is an iid sequence of 2D Gaussian random vectors.
The clutter pdf $f_{\mathrm{c}}\big{(}\bm{z}_{k}^{(m)}\big{)}$ is uniform on
the ROI. The mean number of clutter measurements is
$\mu_{\mathrm{c}}\\!=\\!10$. A birth distribution
$f_{\mathrm{b}}(\bm{x}_{n,k})$ that is uniform on the ROI is considered.
Furthermore, the mean number of newborn targets per time step and the survival
probability are set to $\mu_{\mathrm{b}}=0.01$ and $p_{\mathrm{s}}=0.995$,
respectively.

We use the particle-based implementation of BP for MTT presented in [38] with
$J=5000$ particles for each PT state. The threshold for target confirmation is
$P_{\text{th}}=0.5$. Potential target states are pruned if their existence
probability is below threshold $P_{\text{pr}}=10^{-5}$. We simulated 300 time
steps $n$ and 1000 Monte Carlo realizations for each of the three scenarios.
We used the track-oriented variant of the MHT method. The JPDA and the MHT
methods used a gate validation threshold of 13.82, which corresponds to a
probability of 0.999 that target-originated measurements (when the
corresponding target exists) are within the gate. A gate validation threshold
of 9.21, which corresponds to a probability of 0.99 that target-originated
measurements are within the gate, led to similar results (not shown). The
track confirmation logics of JPDA and MHT are set to 10/16 and 12/24,
respectively. Furthermore, the JPDA and the MHT allow at most 13 missed
detections prior to track termination. These values are tuned to best
represent the algorithms across all three scenarios. The hypothesis depth of
the MHT filter is set to five scans.

To visualize the coalescence and repulsion effects related to the estimated
tracks provided by the JPDA, MHT, and BP methods, we compute the mean distance
of the estimated y-coordinates to the center between the two tracks’
(“D-Center”) as well as the distance between the estimated y-coordinates of
the two tracks (“D-Tracks”). In particular the D-Center is the mean of
$(|\hat{x}_{2,k}^{(1)}|+|\hat{x}_{2,k}^{(2)}|)/2$, averaged over all time
steps and realizations. If more than two tracks are estimated at a certain
time step, then only the two estimates which are closest (in a global nearest
neighbor sense) to the true tracks are considered. Similarly, if one or no
tracks are estimated at certain time steps, the missing state estimates are
excluded from the mean. Furthermore, D-Tracks is the mean of
$|\hat{x}_{2,k}^{(1)}-\hat{x}_{2,k}^{(2)}|$, averaged over all time steps and
iterations. Additional and missing tracks are handled similarly as for
D-Center. D-Tracks can alert to coalescence, while (assuming tracks symmetric
about the origin) D-Center examines the quality of the tracks’ centroid. For
example, if both are zero, the tracks may have coalesced but at least their
joint track is good; but if D-Center is large then even the joint track is
poor.

Fig. 3 shows D-Center and D-Tracks for scenario 1. In Fig. 3(a) it can be seen
that when the two targets move in parallel and in close proximity, the tracks
provided by the MHT yield an increased D-Center compared to the true tracks
due to track repulsion. On the other hand, for tracks provided by the JPDA,
D-Center is reduced; for time steps $200\leq n\leq 250$, the reduced D-Center
due to track coalescence is particularly pronounced after the true tracks have
split again. Furthermore, it can also be observed that the tracks provided by
BP only suffer from a very weak form of repulsion (when true tracks are in
close proximity) and coalescence (after the true tracks have split). For the
D-Tracks in 3(b), it can be seen that the track coalescence effect leads to
target tracks that move towards each other until they have almost identical
y-coordinates. In scenarios 2 and 3, we have observed a similar behavior of
D-Center and D-Tracks (not shown).

As an additional performance metric, we use the general optimal subpattern
assignment (GOSPA) [53] based on the 2-norm as inner metric as well as
parameters $p=1$ and $c=50$. The mean GOSPA—averaged over all simulation
runs—for the three tracking methods in scenario 1 is shown in Fig. 4. We show
the mean total GOSPA in Fig. 4(a) as well as the individual GOSPA
contributions, namely the mean total localization error, the mean missed
target error, and the mean false target error in Figs. 4(b)–(d), respectively.
The estimated MHT tracks yield a slightly increased total GOSPA for times
$n<100$, due to an increase number of missed tracks compared to the estimated
JPDA and BP tracks. This slightly increased number of missed tracks is the
result of MHT parameters that were tuned to yield satisfying performance in
all three scenarios. (Note that for the JPDA method it was not possible to
find parameters that lead to a satisfying performance in scenario 2.) For the
times where the targets are in close proximity, i.e., $100\leq n\leq 200$,
coalescence of estimated JPDA tracks results in a lower localization error
compared to the other time steps and compared to the other methods. However,
right after the true target tracks split, i.e., at times $200<n<240$, the JPDA
yields a significant localization error peak. Due to track repulsion, the
localization error of estimated MHT tracks is increased when targets are in
close proximity. Notably, the localization error related to BP is
significantly reduced compared to the other two methods.

Figure 4: Scenario 1: Mean total GOSPA (a) mean localization error (b), mean
missed target error (c), and mean false target error (d) versus time for the
three considered multitarget tracking methods.

Fig. 5 shows the mean GOSPA error of the JPDA, MHT, and BP methods in scenario
2. For time steps where targets are in close proximity, i.e., $0\leq n\leq
150$, due to track repulsion, the localization and total GOSPA errors of the
MHT are increased compared to those of the BP method. In addition, it can be
seen that for times $n<150$, the missed target and total GOSPA errors of the
JPDA result are significantly increased compared to the other two methods.
This is because the confirmation logic of the JPDA is unable to initialize two
tracks as long the true target states remain in close proximity [1]. The fact
that the JPDA only creates a single track for these times steps, results in a
reduced localization error compared to MHT and BP. As expected this
localization error is 5m on average, which is equal to half the distance
between the two tracks. In scenario 2, BP yields a slightly increased false
tracks error compared to JPDA and MHT. This increased error is mostly due the
generation of short tracks that only consists of few time steps and could be
removed in post-processing. (Note that the tracking confirmation logics used
for JPDA and MHT automatically suppress such short tracks.)

Figure 5: Scenario 2: Mean total GOSPA (a) mean localization error (b), mean
missed target error (c), and mean false target error (d) versus time for the
three considered multitarget tracking methods. Figure 6: Scenario 3: Mean
total GOSPA (a) mean localization error (b), mean missed target error (c), and
mean false target error (d) versus time for the three considered multitarget
tracking methods.

Finally, Fig. 6 shows the mean GOSPA error of the JPDA, MHT, and BP methods in
scenario 3. Here we can again see that for time steps where targets are in
close proximity, i.e., $125\leq n\leq 175$, due to track repulsion, the
localization and total GOSPA error of MHT is increased compared to the other
methods. On the other hand, track coalescence of the JPDA decreases the
localization and total GOSPA error when targets are in close proximity
($125\leq n\leq 175$) and increases them right after the two targets have
split ($175\leq n\leq 225$). An important aspect in a scenario with crossing
tracks are track swaps [45]. A quantification of track swaps based on track-
level metrics [44] will be reported in future work.

##  V Conclusion

In this paper, we reviewed the estimation strategies of multiple hypothesis
tracking (MHT), joint probabilistic data association (JPDA), and the recently
introduced belief propagation (BP) approach for MTT. Since BP performs
marginalization similarly to JPDA, it is expected to suffer from track
coalescence. Our results however indicated that BP-based MTT can mostly avoid
coalescence and repulsion, problems that respectively afflict JPDA and MHT.
This can be explained by the particle-based implementation used by BP-based
MTT that makes it possible to accurately represent multimodal marginal pdfs.
Furthermore, the overconfident nature of BP-based data association “upsells”
the most likely measurement-to-target associations hypothesis and thus leads
to a hybrid form of soft and hard data association that reduces track
coalescence.

## References

  * [1] Y. Bar-Shalom, P. K. Willett, and X. Tian, _Tracking and Data Fusion: A Handbook of Algorithms_. Storrs, CT: Yaakov Bar-Shalom, 2011.
  * [2] S. Challa, M. R. Morelande, D. Mušicki, and R. J. Evans, _Fundamentals of Object Tracking_. Cambridge, UK: Cambridge University Press, 2011.
  * [3] R. Mahler, _Statistical Multisource-Multitarget Information Fusion_. Norwood, MA: Artech House, 2007\.
  * [4] W. Koch, _Tracking and Sensor Data Fusion: Methodological Framework and Selected Applications_. Berlin, Germany: Springer, 2014.
  * [5] F. Meyer, T. Kropfreiter, J. L. Williams, R. Lau, F. Hlawatsch, P. Braca, and M. Z. Win, “Message passing algorithms for scalable multitarget tracking,” _Proc. IEEE_ , vol. 106, no. 2, pp. 221–259, 2018.
  * [6] D. B. Reid, “An algorithm for tracking multiple targets,” _IEEE Trans. Autom. Control_ , vol. 24, no. 6, pp. 843–854, Dec. 1979.
  * [7] Y. Bar-Shalom, “Extension of the probabilistic data association filter in multi-target tracking,” in _Proc. SNETA-74_ , San Diego, CA, Sep. 1974.
  * [8] Y. Bar-Shalom, F. Daum, and J. Huang, “The probabilistic data association filter,” _IEEE Control Syst. Mag._ , vol. 29, pp. 82–100, Dec. 2009.
  * [9] T. Fortmann, Y. Bar-Shalom, and M. Scheffe, “Sonar tracking of multiple targets using joint probabilistic data association,” _IEEE J. Ocean. Eng._ , vol. 8, no. 3, pp. 173–184, Jul. 1983.
  * [10] G. Ferri, A. Munafò, A. Tesei, P. Braca, F. Meyer, K. Pelekanakis, R. Petroccia, J. Alves, C. Strode, and K. LePage, “Cooperative robotic networks for underwater surveillance: an overview,” _IET Radar Sonar Navig._ , vol. 11, no. 12, pp. 1740–1761, 2017.
  * [11] D. Y. Kim, B. Ristic, R. Guan, and L. Rosenberg, “A Bernoulli track-before-detect filter for interacting targets in maritime radar,” _IEEE Trans. Aerosp. Electron. Syst._ , pp. 1–10, Jan. 2021.
  * [12] J. Levinson and al., “Towards fully autonomous driving: Systems and algorithms,” in _Proc. IEEE IV 2011_ , Baden-Baden, Germany, Jun. 2011, pp. 163–168.
  * [13] F. Meyer and J. L. Williams, “Scalable detection and tracking of geometric extended objects,” 2021, arXiv:2103.11279.
  * [14] D. Musicki and R. Evans, “Joint integrated probabilistic data association: JIPDA,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 40, no. 3, pp. 1093–1099, Jul. 2004.
  * [15] J. Vermaak, S. Maskell, and M. Briers, “A unifying framework for multi-target tracking and existence,” in _Proc. FUSION-05_ , Philadelphia, PA, Jul. 2005, pp. 250–258.
  * [16] Musicki, D. and La Scala, B., “Multi-target tracking in clutter without measurement assignment,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 44, no. 3, pp. 877–896, Jul. 2008.
  * [17] P. Horridge and S. Maskell, “Real-time tracking of hundreds of targets with efficient exact JPDAF implementation,” in _Proc. FUSION-06_ , Florence, Italy, Jul. 2006.
  * [18] K. Pattipati, R. L. Popp, and T. Kirubarajan, “Survey of assignment techniques for multitarget tracking,” in _Multitarget-Multisensor Tracking: Applications and Advances_ , Y. Bar-Shalom and W. D. Blair, Eds. Norwood, MA: Artech-House, 2000, vol. 3, ch. 2, pp. 77–159.
  * [19] I. J. Cox and S. L. Hingorani, “An efficient implementation of Reid’s multiple hypothesis tracking algorithm and its evaluation for the purpose of visual tracking,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 18, no. 2, pp. 138–150, Feb. 1996.
  * [20] R. Danchick and G. E. Newnam, “Reformulating Reid’s MHT method with generalised Murty K-best ranked linear assignment algorithm,” _IEE Radar Sonar Navig._ , vol. 153, no. 1, pp. 13–22, Feb. 2006.
  * [21] T. Kurien, “Issues in the design of practical multitarget tracking algorithms,” in _Multitarget-Multisensor Tracking: Advanced Applications_ , Y. Bar-Shalom, Ed. Norwood, MA: Artech-House, 1990, pp. 43–83.
  * [22] S. S. Blackman, “Multiple hypothesis tracking for multiple target tracking,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 19, pp. 5–18, Jan. 2004.
  * [23] C. Morefield, “Application of 0-1 integer programming to multitarget tracking problems,” _IEEE Trans. Autom. Control_ , vol. 22, no. 3, pp. 302–312, Jun 1977.
  * [24] S. Deb, M. Yeddanapudi, K. Pattipati, and Y. Bar-Shalom, “A generalized S-D assignment algorithm for multisensor-multitarget state estimation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 33, no. 2, pp. 523–538, Apr. 1997.
  * [25] R. P. S. Mahler, “Multitarget Bayes filtering via first-order multitarget moments,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 39, no. 4, pp. 1152–1178, Oct. 2003.
  * [26] B.-N. Vo, S. Singh, and A. Doucet, “Sequential Monte Carlo methods for multitarget filtering with random finite sets,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 41, no. 4, pp. 1224–1245, Oct. 2005.
  * [27] B.-T. Vo, B.-N. Vo, and A. Cantoni, “Analytic implementations of the cardinalized probability hypothesis density filter,” _IEEE Trans. Signal Process._ , vol. 55, no. 7, pp. 3553–3567, Jul. 2007.
  * [28] J. L. Williams, “Marginal multi-Bernoulli filters: RFS derivation of MHT, JIPDA and association-based MeMBer,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 51, no. 3, pp. 1664–1687, Jul. 2015.
  * [29] S. Nannuru, S. Blouin, M. Coates, and M. Rabbat, “Multisensor CPHD filter,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 52, no. 4, pp. 1834–1854, Aug. 2016.
  * [30] R. J. Fitzgerald, “Track biases and coalescence with probabilistic data association,” _IEEE Trans. Signal Process._ , vol. 21, no. 6, pp. 822–825, Nov. 1985.
  * [31] S. Coraluppi, C. Carthel, P. Willett, M. Dingboe, O. O’Neill, and T. Luginbuhl, “The track repulsion effect in automatic tracking,” in _Proc. FUSION-09_ , Seattle, WA, USA, Jul. 2009, pp. 2225–2230.
  * [32] H. Blom and E. Bloem, “Probabilistic data association avoiding track coalescence,” _IEEE Trans. Autom. Control_ , vol. 45, no. 2, pp. 247–259, Feb. 2000.
  * [33] H. A. P. Blom, E. A. Bloem, and D. Musicki, “JIPDA: Automatic target tracking avoiding track coalescence,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 51, no. 2, pp. 962–974, Apr. 2015.
  * [34] S. Coraluppi and C. Carthel, “An equivalence-class approach to multiple-hypothesis tracking,” in _Proc. AeroConf-12_ , Big Sky, MT, USA, Mar. 2012, pp. 1–8.
  * [35] L. Svensson, D. Svensson, M. Guerriero, and P. Willett, “Set JPDA filter for multitarget tracking,” _IEEE Trans. Signal Process._ , vol. 59, no. 10, pp. 4677–4691, Oct. 2011.
  * [36] S. Reuter, B.-T. Vo, B.-N. Vo, and K. Dietmayer, “The labeled multi-Bernoulli filter,” _IEEE Trans. Signal Process._ , vol. 62, no. 12, pp. 3246–3260, Jun. 2014.
  * [37] J. L. Williams and R. Lau, “Approximate evaluation of marginal association probabilities with belief propagation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 50, no. 4, pp. 2942–2959, Oct. 2014.
  * [38] F. Meyer, P. Braca, P. Willett, and F. Hlawatsch, “A scalable algorithm for tracking an unknown number of targets using multiple sensors,” _IEEE Trans. Signal Process._ , vol. 65, no. 13, pp. 3478–3493, Mar. 2017.
  * [39] F. Meyer, O. Hlinka, H. Wymeersch, E. Riegler, and F. Hlawatsch, “Distributed localization and tracking of mobile networks including noncooperative objects,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 2, no. 1, pp. 57–71, Mar. 2016.
  * [40] K. R. Pattipati, S. Deb, Y. Bar-Shalom, and R. B. J. Washburn, “A new relaxation algorithm and passive sensor data association,” _IEEE Trans. Autom. Control_ , vol. 37, no. 2, pp. 198–213, Feb. 1992.
  * [41] A. P. Poore and N. Rijavec, “A Lagrangian relaxation algorithm for multidimensional assignment problems arising from multitarget tracking,” _SIAM J. Optimiz._ , vol. 3, no. 3, pp. 544–563, 1993.
  * [42] A. B. Poore and S. Gadaleta, “Some assignment problems arising from multiple target tracking,” _Math. Comput. Model._ , vol. 43, no. 9–10, pp. 1074–1091, 2006.
  * [43] S. Mori, C.-Y. Chong, and K. C. Chang, “Evaluation of data association hypotheses: Non-Poisson i.i.d. cases,” in _Proc. FUSION-04_ , Stockholm, Sweden, Jun. 2004.
  * [44] S. Coraluppi, D. Grimmett, and P. de Theije, “Benchmark evaluation of multistatic trackers,” in _Proc. FUSION-06_ , Florence, Italy, Jul. 2006\.
  * [45] P. Willett, T. Luginbuhl, and E. Giannopoulos, “MHT tracking for crossing sonar targets,” in _Proc. SPIE_ , San Diego, CA, USA, Sep. 2007, pp. 469–480.
  * [46] R. Fitzgerald, “Development of practical PDA logic for multitarget tracking by microprocessor,” in _Multitarget-Multisensor Tracking: Advanced Applications_ , Y. Bar-Shalom, Ed. Norwood, MA: Artech-House, 1990.
  * [47] J. L. Williams, “An efficient, variational approximation of the best fitting multi-Bernoulli filter,” _IEEE Trans. Signal Process._ , vol. 63, no. 1, pp. 258–273, Jan. 2015.
  * [48] R. A. Lau and J. L. Williams, “A structured mean field approach for existence-based multiple target tracking,” in _Proc. FUSION-16_ , Heidelberg, Germany, Jul. 2016, pp. 1111–1118.
  * [49] G. Soldi, F. Meyer, P. Braca, and F. Hlawatsch, “Self-tuning algorithms for multisensor-multitarget tracking using belief propagation,” _IEEE Trans. Signal Process._ , vol. 67, no. 15, pp. 3922–3937, Aug. 2019.
  * [50] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 2, pp. 498–519, Feb. 2001.
  * [51] Y. Weiss and W. T. Freeman, “Correctness of belief propagation in Gaussian graphical models of arbitrary topology,” _Neural Comput._ , vol. 13, no. 10, pp. 2173–2200, 2001.
  * [52] Y. Bar-Shalom, T. Kirubarajan, and X.-R. Li, _Estimation with Applications to Tracking and Navigation_. New York City, NY: Wiley, 2002.
  * [53] A. S. Rahmathullah, Á. F. García-Fernández, and L. Svensson, “Generalized optimal sub-pattern assignment metric,” in _Proc. FUSION-17_ , Xi’an, China, Jul. 2017.

  *[DA]: data association
  *[pdf]: probability density function
  *[PT]: potential target
  *[PTs]: potential target
  *[pdfs]: probability density function
  *[GOSPA]: general optimal subpattern assignment