Poisson multi-Bernoulli mixture trackers: continuity through random finite sets of trajectories

TL;DR

This paper introduces Poisson multi-Bernoulli mixture (PMBM) trajectory filters that ensure continuity over time by modeling multi-target trajectories as random finite sets, bridging the gap between Bayesian filtering and hypothesis-based methods.

Contribution

It derives two novel PMBM trajectory filters that maintain hypothesis structure and continuity, extending the PMBM framework to trajectory estimation.

Findings

Efficiently estimates multi-target trajectories.

Establishes continuity similar to multiple hypothesis tracking.

Shares hypothesis structure with existing PMBM filter.

Abstract

The Poisson multi-Bernoulli mixture (PMBM) is an unlabelled multi-target distribution for which the prediction and update are closed. It has a Poisson birth process, and new Bernoulli components are generated on each new measurement as a part of the Bayesian measurement update. The PMBM filter is similar to the multiple hypothesis tracker (MHT), but seemingly does not provide explicit continuity between time steps. This paper considers a recently developed formulation of the multi-target tracking problem as a random finite set (RFS) of trajectories, and derives two trajectory RFS filters, called PMBM trackers. The PMBM trackers efficiently estimate the set of trajectories, and share hypothesis structure with the PMBM filter. By showing that the prediction and update in the PMBM filter can be viewed as an efficient method for calculating the time marginals of the RFS of trajectories, continuity in the same sense as MHT is established for the PMBM filter.

Equations

$$X=\left(\beta ,\epsilon ,x_{\beta :\epsilon }\right)$$

$$x_{\beta },x_{\beta +1},\dots ,x_{\epsilon -1},x_{\epsilon },$$

$${\mathcal{T}}_k={\uplus }_{(\beta ,\epsilon )\in I_k}\{\beta \}\times \{\epsilon \}\times {\mathcal{X}}^{\epsilon -\beta +1},$$

$$p_{k|k^{\prime }}(X)=p_{k|k^{\prime }}(x_{\beta :\epsilon }|\beta ,\epsilon )P_{k|k^{\prime }}(\beta ,\epsilon ),$$

$$\int p(X)\mathrm{d}X=\iiint p(x_{\beta :\epsilon }|\beta ,\epsilon )P(\beta ,\epsilon )\mathrm{d}{x}_{\beta :\epsilon }\mathrm{d}\epsilon \mathrm{d}\beta$$

$$=\sum _{\beta ,\epsilon }\left[\text{idotsint}p(x_{\beta },\dots ,x_{\epsilon }|\beta ,\epsilon )\mathrm{d}{x}_{\beta }\dots \mathrm{d}{x}_{\epsilon }\right]P(\beta ,\epsilon ).$$

$$\begin{gathered} \int g({\mathbf{X}}_k)\delta {\mathbf{X}}_k\triangleq g(\varnothing )+ \\ \sum _{n=1}^{\infty }\frac{1}{n!}\text{idotsint}g(\{X_k^1,\dots ,X_k^n\})\mathrm{d}{X}_k^1\cdots \mathrm{d}{X}_k^n. \end{gathered}$$

$$G_{\mathbf{X}}[h]$$

$$\triangleq f_{\mathbf{X}}(\varnothing )+\sum _{n=1}^{\infty }\frac{1}{n!}\text{idotsint}\left[\prod _{i=1}^{n}h(X^i)\right]$$

$$\times f_{\mathbf{X}}(\{X^1,\dots ,X^n\})\mathrm{d}{X}^1\cdots \mathrm{d}{X}^n,$$

$$f^{\mathrm{ppp}}(\mathbf{X})$$

$$G^{\mathrm{ppp}}[h]$$

$$f^{\mathrm{ber}}(\mathbf{X})$$

$$G^{\mathrm{ber}}[h]$$

$$G_{k|k^{\prime }}[h]=\exp \left\{⟨{\lambda }_{k|k^{\prime }}^{\mathrm{u}};h(X_k)-1⟩\right\}$$

$$\times \sum _{a\in {\mathcal{A}}^{k|k^{\prime }}}{w}_{k|k^{\prime }}^a\prod _{i\in {\mathbb{T}}_{k|k^{\prime }}}\left(1-r_{k|k^{\prime }}^{i,a^i}+r_{k|k^{\prime }}^{i,a^i}⟨f_{k|k^{\prime }}^{i,a^i};h⟩\right).$$

$${\lambda }_{k|k^{\prime }}^{\mathrm{u}},\left\{w_{k|k^{\prime }}^{i,a^i},r_{k|k^{\prime }}^{i,a^i},f_{k|k^{\prime }}^{i,a^i}\right\}$$

$${\lambda }_k^{\mathrm{B}}(X)$$

$${\lambda }_k^{\mathrm{B},x}(x_{k:k}|k,k)$$

$$P_k^{\mathrm{S}}(X)=P^{\mathrm{S}}(x_{\epsilon }){\Delta }_k(\epsilon ).$$

$$f_{k|k-1}^{\mathrm{c}}(\mathbf{X}|{\mathbf{X}}^{\prime })=$$

$$\displaystyle\left\{\begin{array}[]{lcl}1,&&\mathbf{X}^{\prime}=\emptyset,\mathbf{X}=\emptyset\\ 1-P_{k-1}^{\rm S}(X^{\prime}),&\text{if}&\mathbf{X}^{\prime}=\{X^{\prime}\},\mathbf{X}=\emptyset\\ P_{k-1}^{\rm S}(X^{\prime})\pi^{\rm c}(X|X^{\prime}),&&\mathbf{X}^{\prime}=\{X^{\prime}\},\mathbf{X}=\{X\}\\ 0,&&\text{otherwise}\end{array}\right.$$

$${\pi }^{\mathrm{c}}(X|X^{\prime })={\pi }^{\mathrm{c},x}(x_{\beta :\epsilon }|\beta ,\epsilon ,X^{\prime }){\Delta }_{{\epsilon }^{\prime }+1}(\epsilon ){\Delta }_{{\beta }^{\prime }}(\beta ),$$

$${\pi }^{\mathrm{c},x}(x_{\beta :\epsilon }|\beta ,\epsilon ,X^{\prime })={\pi }^x(x_{\epsilon }|x_{{\epsilon }^{\prime }}^{\prime }){\delta }_{x_{{\beta }^{\prime }:{\epsilon }^{\prime }}^{\prime }}(x_{\beta :\epsilon -1}),$$

$$G_{k|k-1}^{\mathrm{c}}[h|\mathbf{X}]=e^{⟨{\lambda }_k^{\mathrm{B}};h-1⟩}{\left(1-P_{k-1}^{\mathrm{S}}+P_{k-1}^{\mathrm{S}}⟨{\pi }^{\mathrm{c}};h⟩\right)}^{\mathbf{X}}$$

$${\lambda }_{k|k-1}^{\mathrm{u}}(X_k)$$

$$n_{k|k-1}$$

$$h_{k|k-1}^i$$

$$w_{k|k-1}^{i,a^i}$$

$$r_{k|k-1}^{i,a^i}$$

$$f_{k|k-1}^{i,a^i}$$

$$f_{k|k-1}^{\mathrm{a}}(\mathbf{X}|{\mathbf{X}}^{\prime })=$$

$$\displaystyle\left\{\begin{array}[]{lcl}1&&\mathbf{X}^{\prime}=\emptyset,\mathbf{X}=\emptyset\\ \pi^{\rm a}(X|X^{\prime})&\text{if}&\mathbf{X}^{\prime}=\{X^{\prime}\},\mathbf{X}=\{X\}\\ 0&&\text{otherwise}\end{array}\right.$$

$${\pi }^{\mathrm{a}}(X|X^{\prime })={\pi }^{\mathrm{a},x}(x_{\beta :\epsilon }|\beta ,\epsilon ,X^{\prime }){\pi }^{\epsilon }(\epsilon |\beta ,X^{\prime }){\Delta }_{{\beta }^{\prime }}(\beta ),$$

$$\displaystyle\pi^{\varepsilon}(\varepsilon|\beta,X^{\prime})=\left\{\begin{array}[]{ll}1,&\varepsilon={\varepsilon}^{\prime}<k-1\\ 1-P_{k-1}^{\rm S}(X^{\prime}),&\varepsilon={\varepsilon}^{\prime}=k-1\\ P_{k-1}^{\rm S}(X^{\prime}),&\varepsilon={\varepsilon}^{\prime}+1=k\\ 0,&\text{otherwise}\end{array}\right.,$$

$${\pi }^{\mathrm{a},x}(x_{\beta :\epsilon }|\beta ,\epsilon ,X^{\prime })=$$

$$\displaystyle\left\{\begin{array}[]{ll}\delta_{x_{\beta^{\prime}:\varepsilon^{\prime}}^{\prime}}(x_{\beta:\varepsilon}),&\varepsilon=\varepsilon^{\prime}\\ \pi^{x}(x_{\varepsilon}|x^{\prime}_{\varepsilon^{\prime}})\delta_{x^{\prime}_{\beta^{\prime}:\varepsilon^{\prime}}}(x_{\beta:\varepsilon-1}),&\varepsilon=\varepsilon^{\prime}+1\end{array}\right..$$

$$G_{k|k-1}^{\mathrm{a}}[h|\mathbf{X}]=e^{⟨{\lambda }_k^{\mathrm{B}};h-1⟩}{⟨{\pi }^{\mathrm{a}};h⟩}^{\mathbf{X}}$$

$${\lambda }_{k|k-1}^{\mathrm{u}}(X_k)$$

$$n_{k|k-1}$$

$$h_{k|k-1}^i$$

Full text

I Introduction

Multi-target tracking (mtt) is a challenging problem due to the unknown correspondence of measurements and targets, referred to as data association. Each new measurement received could be the continuation of some previously detected target, the first detection of a new target, or a false alarm.

The multiple hypothesis tracker (mht) introduced in [1] is a key method for addressing problems of this type. The algorithm maintains a series of global hypotheses for each possible measurement-target correspondence as measurements are received, along with a conditional state distribution for each target under each hypothesis. Each measurement could be hypothesised to represent the first detection of a new target, where the number of newly detected targets was given a Poisson distribution in order to provide a Bayesian prior.

In [2], the model was made rigorous through random finite sequences, under the assumption that the number of targets present is constant but unknown, with a Poisson prior. Target state sequences were formed under each global hypothesis, and the Poisson distribution of targets remaining to be detected provided a Bayes prior for events involving newly detected targets. Hypotheses were constructed as being data-to-data, since no a priori data was assumed on target identity. This formulation was revisited in [3], and compared to random finite sets, and finite point processes.

The framework of random finite sets (rfss) was developed to provide a systematic methodology for dealing with problems involving an uncertain number of targets [4]. In [5], the Bayes filter for Poisson birth models was derived using rfss, obtaining a result somewhat similar to [2], involving a Poisson distribution representing targets which are hypothesised to exist but remain to be detected, and a multi-Bernoulli mixture (mbm) representing targets that have been detected at some stage. Adopting the common modelling assumptions employed in rfs (see, e.g., [6]), target appearance and disappearance was modelled. The resulting, so called, pmbm filter, did not formally establish track continuity, although a hypothesis structure similar to mht has been observed [7, 8].

mht was extended in [9, 10] to address a time-varying number of targets, incorporating birth of targets which are not immediately detected, necessitating equivalence classes of indistinguishable hypotheses. In comparison, in [5] the hypotheses are purely data-to-data assignments, the hypothesis space explicitly includes targets that remain to be detected, and target death subsequent to the final detection.

The labelled rfs formulation of [11] addresses the lack of continuity of unlabelled rfss by incorporating a label (uniqueness of which is ensured through the model) into the target state; continuity is maintained by connecting estimates with the same label at different times. PGFL expressions for the joint target-measurement density over a time interval were given in [12], including modelling of target birth, target death, and target spawning using branching and immigration processes. A solution was proposed using saddle point approximation, however, a practical implementation was not provided in [12]. In [13], a framework for performing inference on problems involving both distinguishable (previously detected) and indistinguishable (never detected) targets, was developed and applied to space situation awareness. Targets were characterized by sequences of measurements, under the non-standard assumption that the targets are “immediately detected upon entering the surveillance scene” [13, A.5].

A formulation of the target tracking problem as an unlabelled rfs of trajectories was provided recently in [14, 15]. Within this set of trajectories framework, the goal of Bayesian mtt is to compute the posterior density over the set of trajectories. One approach to do this is to let the multi-target state be the set of trajectories, and compute the posterior density recursively. Assuming standard point target motion and measurement models, in this paper we formulate two different set of trajectory problems: one in which the set of current (i.e., remaining) trajectories is tracked, and one in which the set of all trajectories up to the current time is tracked. For the two problem formulations, we assume multi-target densities of the pmbm form and derive the pmbm predictions and the pmbm update. This results in two different pmbm set of trajectories filters: we call these tracking algorithms PMBM trackers, to distinguish them from the pmbm filter [5], which is for sets of target states.

The prediction and update of pmbm trackers are the main contribution in the paper. In addition to this, we discuss the pmbm trackers in relation to the pmbm filter [5], the mht [1, 2], and the $\delta$-glmb filter [11]. Interestingly, the relations between the pmbm trackers and the pmbm filter show that the pmbm filter implicitly contains trajectory information. We also show how the pmbm trackers can be implemented for linear Gaussian models, and present results from a simulation study where we compare the results to the $\delta$-glmb filter [11].

Importantly, our paper shows that it is possible to compute the density over the set of trajectories even for complicated problems; in [10, Sec. 2] computing the set of trajectories density was thought to be “extremely challenging…even for small problems”. Further, we compute the density without requiring the assumption that targets are detected with probability one upon entering the scene.

II Problem formulation

To clearly differentiate between a target state at a single time and a sequence of target states, we let target denote the state at some time, and we let trajectory denote a sequence of states. Thus, “set of targets” and “target rfs” refers to formulations involving rfss of target states at a single time (e.g., the unlabelled rfs formulation in [5]), and “set of trajectories” and “trajectory rfs” refers to formulations involving a rfs of trajectories, or sequences of states.

Let $x_{k}$ denote a target state at time $k$, and let $z_{k}$ denote a measurement at time $k$. We utilise the standard multi-target dynamics model, defined in Assumption 1, and the standard point target measurement model, defined in Assumption 2.

Assumption 1**.**

The multiple target state evolves according to the following time dynamics process: targets arrive at each time according to a non-homogeneous Poisson Point Process (PPP) with birth intensity $\lambda^{\mathrm{b}}(x_{k})$, independent of existing targets; targets depart according to iid Markovian processes; the survival probability in state $x_{k}$ is $P^{\rm S}(x_{k})$; target motion follows iid Markovian processes; the single-target transition density is $\pi^{x}(x_{k}|x_{k-1})$.

Assumption 2**.**

The multiple target measurement process is as follows: each target may give rise to at most one measurement; probability of detection in state $x_{k}$ is $P^{\rm D}(x_{k})$; each measurement is the result of at most one target; false alarm measurements arrive according to a non-homogeneous PPP with intensity $\lambda^{\mathrm{FA}}(z_{k})$, independent of targets and target-related measurements; each target-derived measurement is independent of all other targets and measurements conditioned on its corresponding target; the single target measurement likelihood is $f(z_{k}|x_{k})$.

There are many ways in which a Multiple Target Tracking (mtt) problem can be formulated. In this paper, we focus on the following two variants:

The set of current trajectories: the objective is to estimate the trajectories of the targets who are present in the surveillance area at the current time. 2. 2.

The set of all trajectories: the objective is to estimate the trajectories of all targets that have passed through the surveillance area at some point between time [math] and the current time, i.e., both the targets who are present in the surveillance area at the current time, and the targets that have left the surveillance area (but were in the surveillance area at at least one previous time).

Note that both problem formulations are equally valid; which one is relevant depends on the application. Indeed, in some applications, a solution with only current states (no trajectories) is the objective.

III Random finite sets of trajectories

Compact representations of the multiple trajectory density for the standard point target models, cf. Assumptions 1 and 2, are obtained using sets of trajectories, developed in the rfs formalism in [15]. Expressions for integration, Bayes prediction and Bayes update may be found there. In this section, we first review the trajectory state representation, and then present densities and probability generating functionals for sets of trajectories.

III-A Single trajectory

Let $\mathcal{X}$ represent the base state space, e.g., $\mathcal{X}=\mathbb{R}^{4}$, representing position and velocity in two dimensions. We use the trajectory state model proposed in [14, 15], in which the trajectory state is a tuple

[TABLE]

where $\beta$ is the discrete time of the trajectory birth, i.e., the time the trajectory begins; $\varepsilon$ is the discrete time of the trajectory’s most recent state, i.e., the time the trajectory ends. If $k$ is the current time, $\varepsilon<k$ means that the trajectory ended at time $\varepsilon$, and $\varepsilon=k$ means that the trajectory is ongoing; $x_{\beta:\varepsilon}$ is, given $\beta$ and $\varepsilon$, the sequence of states

[TABLE]

where $x_{k}\in\mathcal{X}$ for all $k\in\{\beta,\ldots,\varepsilon\}$. This gives a trajectory of length $\ell=\varepsilon-\beta+1$ time steps. The trajectory state space at time $k$ is [15]

[TABLE]

where $\uplus$ denotes disjoint set union, $I_{k}=\{(\beta,\varepsilon):0\leq\beta\leq\varepsilon\leq k\}$ and $\mathcal{X}^{\ell}$ denotes $\ell$ Cartesian products of $\mathcal{X}$. The trajectory state density factorises as follows111“$p_{k|k^{\prime}}$” is the density of $X_{k}$ given measurements up to and including time $k^{\prime}\leq k$.

[TABLE]

where, if $\varepsilon<\beta$, then $P_{k|k^{\prime}}(\beta,\varepsilon)$ is zero. Integration is performed as follows [15],

[TABLE]

III-B Multiple trajectories

We now move to the case of a set of trajectories. Analogously to a set of targets, a set of trajectories is denoted as $\mathbf{X}_{k}\in\mathcal{F}(\mathcal{T}_{k})$, where $\mathcal{F}(\mathcal{T}_{k})$ is the set of all finite subsets of $\mathcal{T}_{k}$. Let $g(\mathbf{X}_{k})$ be a function on a set of trajectories $\mathbf{X}_{k}$. Integration over sets of trajectories is defined as regular set integration:

[TABLE]

The set of trajectories density function $f_{\mathbf{X}}(\mathbf{X})$ is defined analogously to the set of targets density function. The probability generating functional (pgfl), see, e.g., [4], is a useful tool for manipulating and understanding rfs densities. The pgfl for a trajectory rfs density is defined like the pgfl for a target rfs density: using a test function $h(X):\mathcal{T}\rightarrow\mathbb{R}$, the pgfl is

[TABLE]

where $h^{\mathbf{X}}$ is set power, defined as $h^{\mathbf{X}}=1$ if $\mathbf{X}=\emptyset$ and $h^{\mathbf{X}}=\prod_{X\in\mathbf{X}}h(X)$ if $\mathbf{X}\neq\emptyset$.

A trajectory ppp is analogous to a target ppp, and has set density and pgfl

[TABLE]

where $\left\langle f;g\right\rangle=\int f(X)g(X)\mathrm{d}X$. The trajectory ppp intensity $\lambda(\cdot)$ is defined on the trajectory state space $\mathcal{T}_{k}$, i.e., realisations of the ppp are trajectories with a birth time, a time of the most recent state, and a state sequence.

A trajectory Bernoulli process is analogous to a target Bernoulli process, and has set density and pgfl

[TABLE]

Here, $f(\cdot)$ is a trajectory state density (4), and $r$ is the Bernoulli probability of existence. Together, $r$ and $f(\cdot)$ can be used to find the probability that the target trajectory existed at a specific time, or find the probability that the target state was in a certain area at a certain time. Trajectory mb rfs and trajectory mbm rfs are both defined analogously to target mb rfs and target mbm rfs: a trajectory mb is the disjoint union of a multiple trajectory Bernoulli rfs; trajectory mbm rfs is an rfs whose density is a mixture of trajectory mb densities.

IV PMBM trackers

As in [5], we hypothesise a multi-target conjugate prior222Multi-target conjugate prior was defined in [11] as meaning that ”…if we start with the proposed conjugate initial prior, then all subsequent predicted and posterior distributions have the same form as the initial prior.” of the following Poisson Multi-Bernoulli Mixture (pmbm) form, and observe that the form is preserved through prediction and update. The pmbm pgfl is

[TABLE]

This pgfl form states that the set of trajectories $\mathbf{X}_{k}$ is an independent union of a ppp with intensity $\lambda^{\mathrm{u}}_{k|k^{\prime}}$, and a multi-Bernoulli mixture (mbm) with Bernoulli parameters $r_{k|k^{\prime}}^{i,a^{i}}$ and $f_{k|k^{\prime}}^{i,a^{i}}(\cdot)$.

The structure of the trajectory pmbm (11) is the same as the structure of the target pmbm from [5]. The ppp represents trajectories that are hypothesised to exist, but have never been detected, i.e., no measurement has been associated to them. We present “track oriented” (to) pmbm trackers, where a track is initiated for each measurement. In the mbm in (11), ${\mathbb{T}}_{k|k^{\prime}}$ is a track table with $n_{k|k^{\prime}}$ tracks, $a\in{\cal A}^{k|k^{\prime}}$ is a possible global data association hypothesis, and for each global hypothesis $a$ and each track $i\in{\mathbb{T}}_{k|k^{\prime}}$, $a^{i}$ indicates which track hypothesis is used in the global hypothesis. For each track, there are $n_{k|k^{\prime}}^{i}$ single trajectory hypotheses. The weight of the global hypothesis $a$ is $w_{k|k^{\prime}}^{a}\propto\prod_{i\in{\mathbb{T}}_{k|k^{\prime}}}w_{k|k^{\prime}}^{i,a^{i}}$, where $w_{k|k^{\prime}}^{i,a^{i}}$ is the weight of single trajectory hypothesis $a^{i}$ from track $i$. The set of all measurement indices up to time $k$ is denoted ${\cal M}^{k}$, and ${\cal M}^{k}(i,a^{i})$ is the history of measurements that are hypothesised to belong to hypothesis $a^{i}$ from track $i$. Due to page length constraints, further elaboration is omitted; please refer to [5] for details and discussion.

Depending on the problem formulation, the multi-target state $\mathbf{X}_{k}$ that we are interested in is different. For the set of current trajectories, $\mathbf{X}_{k}$ is the set of trajectories for which $0\leq\beta\leq\varepsilon=k$. For the set of all trajectories, $\mathbf{X}_{k}$ is the set of trajectories for which $0\leq\beta\leq\varepsilon\leq k$. A pmbm process is defined by the parameters of the density/pgfl,

[TABLE]

In the following subsections we will show how the parameters are predicted and updated, in order to track either the set of current trajectories, or the set of all trajectories.

IV-A Prediction step

In this section we describe the time evolution of the set of trajectories. A standard ppp birth model is used, i.e., target birth at time step $k$ is modelled by a ppp, with intensity

[TABLE]

where $\Delta(\cdot)$ denotes Kronecker’s delta function. Note that it is possible to have alternative birth models, such as mb birth, or mbm birth. For those cases, the spatial densities of the Bernoulli birth components would be of the same form as in (13), i.e., for birth at time $k$ we have $\beta=\varepsilon=k$.

For the transition of existing targets, there is one alternative for each of the two problem formulations. For the set of current trajectories, $P^{\rm S}(\cdot)$ is used in a way that is typical for tracking a set of targets, see, e.g., [5, 11]. For the set of all trajectories, $P^{\rm S}(\cdot)$ is used as in [15]. The trajectory state dependent probability of survival at time $k$ is defined as

[TABLE]

IV-A1 Transition model for the set of current trajectories

The Bernoulli rfs transition density without birth is

[TABLE]

where $\delta(\cdot)$ denotes Dirac’s delta function. In this model, $P^{\rm S}(\cdot)$ is used as follows. If the target disappears, or “dies” (case $\mathbf{X}^{\prime}=\{X^{\prime}\},\mathbf{X}=\emptyset$ in (15a)), then the entire trajectory will no longer be a member of the set of current trajectories. If the target survives, then the trajectory is extended by one time step. The pgfl of the multiple object transition density is

[TABLE]

The resulting prediction step is given in the theorem below.

Theorem 1**.**

Assume that the distribution from the previous time step $G_{k-1|k-1}[h]$ is the pgfl of the form given in (11), and that the transition model is of the kind (16). Then the predicted distribution for the next step $G_{k|k-1}[h]$ is of the form (11), with:

[TABLE]

IV-A2 Transition model for the set of all trajectories

The Bernoulli rfs transition density without birth is

[TABLE]

In this model, the interpretation of the probability of survival is that it governs whether or not the trajectory ends, or if it extends by one more time step. However, importantly, regardless of whether or not the trajectory ends, the trajectory remains in the set of all trajectories.

In this case, the pgfl of the multiple object transition density with ppp birth is

[TABLE]

Note the difference to (16), which comes from the way in which the probability of survival is used. In this case, the target always remains in the set of targets, even if its trajectory has ended. The prediction step is presented in the theorem below.

Theorem 2**.**

Assume that the distribution from the previous time step $G_{k-1|k-1}[h]$ is of the form given in (11), and that the transition model is of the kind (19). Then the predicted distribution for the next step $G_{k|k-1}[h]$ is of the form (11), with:

[TABLE]

IV-B Update step

The target measurement model of assumption 2 is extended by defining a Bernoulli measurement density as follows:

[TABLE]

The clutter is modelled as a ppp with intensity $\lambda_{k}^{\rm FA}(z)$, and it follows that the measurement pgfl is

[TABLE]

We see that $G_{k}[g|\mathbf{X}_{k}]$ is a trajectory equivalent to the standard point target pgfl, and thus the trajectory measurement update is analogous to the target measurement update in [5].

Theorem 3**.**

Assume that the predicted distribution $G_{k|k-1}[h]$ is of the form given in (11), and that the measurement model is of the kind (22). Then, the updated distribution $G_{k|k}[h]$ (updated with the measurement set $Z_{k}=\{z_{k}^{1},\dots,z_{k}^{m_{k}}\}$) is of the same form, with $n_{k|k}=n_{k|k-1}+m_{k}$,

[TABLE]

For tracks continuing from previous time steps ($i\in\{1,\dots,n_{k|k-1}\}$), a hypothesis is included for each combination of a hypothesis from a previous time and either a missed detection or an update using one of the $m_{k}$ new measurements, such that the number of hypotheses becomes $h^{i}_{k|k}=h^{i}_{k|k-1}(1+m_{k})$. For missed detection hypotheses ($i\in\{1,\dots,n_{k|k-1}\}$, $a^{i}\in\{1,\dots,h_{k|k-1}\})$:

[TABLE]

For hypotheses updating existing tracks ($i\in\{1,\dots,n_{k|k-1}\}$, $a^{i}=\tilde{a}^{i}+h^{i}_{k|k-1}j$, $\tilde{a}^{i}\in\{1,\dots,h^{i}_{k|k-1}\}$, $j\in\{1,\dots,m_{k}\}$, i.e., the previous hypothesis $\tilde{a}^{i}$, updated with measurement $z_{k}^{j}$):333A hypothesis at the previous time with $r^{i,a^{i}}_{k|k-1}=0$ need not be updated since the posterior weight in (26b) would be zero. For simplicity, the hypothesis numbering does not account for this exclusion.

[TABLE]

Finally, for new tracks, $i\in\{n_{k|k-1}+j\}$, $j\in\{1,\dots,m_{k}\}$ (i.e., the new track commencing on measurement $z_{k}^{j}$),

[TABLE]

IV-C Density/intensity representation and estimation

Note that the birth ppp intensity (13) is an un-normalized mixture density of the form

[TABLE]

where each mixture component is characterised by a weight $w^{j}$, a distinct birth time $b^{j}$, a distinct most recent time $e^{j}$ where $b^{j}\leq e^{j}$ for all $j$, and a state sequence density $p^{j}(\cdot)$. For the weights we have that $\sum_{j}w^{j}=1$ if $p(\cdot)$ is a density, and $\sum_{j}w^{j}\geq 0$ if $p(\cdot)$ is an intensity function, e.g., a ppp intensity. Because of the assumption that the trajectory birth intensity is an un-normalized density of the form (28), it follows that the undetected intensity $\lambda^{\mathrm{u}}_{k|k^{\prime}}(X_{k})$, and all Bernoulli densities $f_{k|k^{\prime}}^{i,a^{i}}(X_{k})$, will be of the form (28). This type of state density facilitates simple representations for the state sequence $x_{\beta:\varepsilon}$, conditioned on $\beta$ and $\varepsilon$.

Trajectory estimation, or trajectory extraction, is the process of obtaining estimates of the set of trajectories (or set of targets) from the multi-target density. A typical approach to mtt estimation is to select a certain global hypothesis and report estimates from it. In this paper, we use a pmbm filter estimator—e.g., any of the three discussed in [16, Sec. VI] is valid—to identify a global hypothesis from which to extract estimates, and then extract the full trajectories directly. Note that alternative estimators are possible, representing different trade-offs between tracking accuracy and computational cost. Comparing estimators will be adressed in future work.

IV-D Structure of PMBM trackers

Two different pmbm trackers result from the theorems: a pmbm tracker for the set of current trajectories is given by the prediction in Theorem 1 and the update in Theorem 3; a pmbm tracker for the set of all trajectories is given by the prediction in Theorem 2 and the update in Theorem 3. Note that both pmbm trackers are to. For each measurement, a potential new track is initiated, see (27). As in the pmbm filter [5], for each track there is a hypothesis tree, where each hypothesis corresponds to different data association sequences for the track. The predictions (17) and (20) preserve the number of tracks and the number of hypotheses444The prediction (20f) results in additional mixture components (28)., meaning that the predictions can be implemented without approximation. In the update, additional hypotheses are created, as indicated in (25) and (26).

The Bernoulli probabilities of existence $r$ have different meanings in the two trackers: for the set of current trajectories problem formulation, $r$ is the probability that the trajectory exists at the current time step $k$ and has not ended yet; in the set of all trajectories problem formulation, $r$ represents the probability that the trajectory existed at any time between [math] and the current time step $k$.

In Theorem 1, the predicted Bernoulli existence probability (17e) is the product of the posterior probability and $\langle P_{k}^{\rm S};f\rangle$. This results in a predicted probability $r\leq 1$, where equality only holds for the atypical choice $P_{k}^{\rm S}(x_{k})=1$. The predicted Bernoulli density (17f) results in a density in which $\varepsilon=k$ with probability $=1$ and the state sequence is augmented by one more time step. This follows directly from the set of current trajectories problem formulation: the filter does not maintain trajectories that ended before time $k$, because those are not current trajectories. Notice the important difference in Theorem 2: here the Bernoulli existence probability is constant, and $P_{k}^{\rm S}(\cdot)$ enters in the integrals with $\pi^{\rm a}$. In the set of all trajectories problem formulation, the existence probability represents existence from time [math] to time $k$, and this is unaffected by whether or not the trajectory ends at this time or not, or if it has already ended. The probabilities of different $\varepsilon$, and the corresponding state sequences, are captured in the Bernoulli density (20f).

In Theorem 3, the updated Bernoulli existence probability decreases following a missed detection (25c), unless the predicted probability is $=1$. For the set of current trajectories, the predicted probability of existence (17e) is typically $<1$, however, for the set of all trajectories the predicted probability of existence can be $=1$. However, remember that trajectory death ($\varepsilon<k$) is modelled in the Bernoulli densities (25d) via the mixture components, cf. (28). For the detection updates (26) and (27), an implication of $P^{\rm D}(\cdot)$ (21g) is that a detected trajectory must necessarily exist at time $k$. It follows that, for both problem formulations, in the updated state densities (26d) and (27f) the probability of $\varepsilon=k$ is $=1$.

V Discussion

V-A Relations to PMBM filter and MHT

Theorem 1 is a direct re-statement of the pmbm filter prediction, see [5, Thm. 1], to the trajectory rfs dynamics model. Further, Theorem 3 is a direct re-statement of the pmbm filter update, see [5, Thm. 2], to the trajectory rfs measurement model. Note that regardless of which problem formulation is considered, current trajectories or all trajectories, the update is the same. An important difference between the work presented here and the work in [5] is that the former is for trajectory sets, whereas the latter is for target sets.

The predictions and update of the pmbm trackers are to, with one track initiated for each measurement. By comparing the track hypotheses in the two to-pmbm trackers with the track hypotheses in the to-pmbm filter [5], we see that the hypothesis structure is identical for the prediction and update steps. This can be seen by comparing the components of the distribution of the form (11), and by considering the marginalisation of the trajectories to leave only the current time step $k$. The result of the marginalisation of a trajectory Bernoulli density with probability of existence $r_{k|k^{\prime}}$ and trajectory density $p_{k|k^{\prime}}(X)$ is a target Bernoulli density, with probability of existence

[TABLE]

If $r^{m}=0$, a state density need not be defined, because the target does not exist at time $k$.

Regarding the general relation between target densities and trajectory densities, a result was presented in [15, Thm. 11] showing that if a trajectory rfs distribution is marginalised to the final time to obtain a target rfs distribution, and prediction and update steps are performed on that target rfs distribution, then the result is equivalent to the marginalisation of the trajectory rfs distribution after prediction and update using the equivalent trajectory models.

Using the marginalisation (29), it can be shown that marginalising each predicted/updated trajectory hypothesis is the same as predicting/updating the marginalised trajectory hypothesis; the details are omitted due to page length constraints. From this, it follows that the marginalisation of the posterior pmbm tracker density is equivalent to the pmbm filter update of the marginalisation of the prior pmbm tracker density. The same type of relation can be shown for the prediction.

Since the system state in the pmbm filter [5] was formally an unlabelled set of target states, track continuity was not provided formally, yet a hypothesis structure similar to mht was observed, which was exploited to provide an implicit form of track continuity. The mht and the pmbm filter were compared in detail in [7, 8], and conditions for equivalence were identified. It follows from the pmbm tracker/filter-relations established here that the pmbm filter [5] can be understood to operate implicitly on the trajectories: it marginalises the state sequence in each prediction, and it only ever explicitly outputs the current time step. This shows the implicit track continuity in the pmbm filter. Note also that the general similarities between the mht and the pmbm filter carry over to the pmbm trackers.

V-B Track initiation

The track initiation in the pmbm trackers is of the same type as in the pmbm filter [5], as well as in the to-mht: a trajectory is defined by the measurement that initiates it. In tracking algorithms based on labelled multi-Bernoulli densities, such as the $\delta$-glmb filter [11], a trajectory is defined by the labelled Bernoulli birth component.555Specifically, each Bernoulli birth component is uniquely labelled at birth, which is modelled through the prediction step; refer to [11] for details.

One disadvantage of labelled birth is that, by its intrinsic modelling properties, it does not permit targets to be indistinguishable up until first detection. This is a significant practical drawback in cases where targets are actually indistinguishable before being detected for the first time, e.g.: when the first detection occurs a long time after target birth; or when target births are iid, which is a common modelling assumption. In fact, in the majority of cases, it is an essential character of the tracking problem that targets are indistinguishable up until first detection, after which continuity is able to be maintained (for as long as sufficient information persists). Thus, we consider it to be an important advantage that in the pmbm tracker, the birth is modelled by a ppp (with intensity (13)).

Note that after the first detection, the targets are distinguishable. From a trajectory and a given data association hypothesis, we can infer at all times the location of the target. Thus, there is no need to label the target upon initialisation.

V-C Estimation

The to-mht provides estimates of a target trajectory by estimating the current state at each time, and connecting estimates from different times. The $\delta$-glmb filter estimates trajectories in a similar fashion: the set of current target states is estimated at each time, and the estimates from different times are connected into trajectories using the labels. Both of these methods can be understood as answering the trajectory question “for the chosen data association hypothesis, what is the current location of the target that was initiated by measurement $z_{\beta}$ at time $\beta$ (in the case of MHT), or was initiated by birth component $l_{\beta}$ at time $\beta$ (in the case of $\delta$-GLMB)?”

The estimation in the pmbm tracker (Section IV-C) can be understood as answering the trajectory question “for the chosen data association hypothesis, what is the state sequence of the target that was initiated by measurement $z_{\beta}$ at time $\beta$?” A benefit of this type of estimation is that instead of appending the sequence of estimates by one more state estimate, it estimates the full state sequence, i.e., smoothed estimates. This is possible, thanks to the fact that each global hypothesis contains full trajectory information.

Thus, the pmbm tracker produces full trajectory estimates upon receipt of each new set of measurements, whereas the pmbm filter, the mht, and the $\delta$-glmb filter produce estimates of the latest time. The beginning of the trajectory will reveal the origin of the target, and thus correspondence between trajectories estimated at different times is established. Practically, these linkages can be made simply through hypothesis metadata, similar to mht, and if all that is required is to estimate the current state at each time, the pmbm filter suffices, since it may be interpreted as implicitly operating on the full trajectory (cf. Section V-A).

VI Linear Gaussian implementation

Let the single target transition density and measurement likelihoods both be linear-Gaussian,

[TABLE]

and let the probabilities of detection and of survival be constant, $P_{k}^{\rm D}(\cdot)=P^{\rm D}$ and $P_{k}^{\rm S}(\cdot)=P^{\rm S}$. Further, let the state sequence, conditioned on $\beta$ and $\varepsilon$, be Gaussian distributed with mean $m$ and covariance $P$. The state sequence density can be expressed in information form,

[TABLE]

with information vector $y_{k|k^{\prime}}=P_{k|k^{\prime}}^{-1}m_{k|k^{\prime}}$ and information matrix $Y_{k|k^{\prime}}=P_{k|k^{\prime}}^{-1}$.

Let $y_{k|k}$ and $Y_{k|k}$ be the posterior information parameters. Given a linear Gaussian motion model (30), the predicted trajectory density has information vector and information matrix [17]

[TABLE]

where $n_{x}$ is the dimension of the single-time-step state $x$, $\mathbb{O}_{m\times n}$ is an $m$ by $n$ all-zero matrix, and $Y^{[a:b,c:d]}$ denotes the part the information matrix with rows for time steps $a$ to $b$ and columns for time steps $c$ to $d$. Given a linear Gaussian measurement model (30), the posterior trajectory density has information vector and information matrix [17]

[TABLE]

Note that in the prediction (32), both the information vector and the information matrix are augmented to account for the additional time step that, following the prediction, is included by the trajectory. A key result is that the bottom left and top right corners of the predicted information matrix (32b) are exactly zero, and that the update (33b) only affects the part of the information matrix that is related to the current state. This means that the information matrix is exactly sparse, a direct consequence of the Markov property associated with $\pi^{x}(\cdot|\cdot)$. The number of non-zero elements in the information matrix increases linearly: a trajectory from $\beta$ to $\varepsilon$ includes $\ell=\varepsilon-\beta+1$ time steps. In the information matrix this leads to $n_{x}^{2}\left(3\ell-2\right)$ non-zero elements. In comparison, with a Gaussian representation (with a mean vector and covariance matrix) the covariance is full and has $n_{x}^{2}\ell^{2}$ non-zero elements. Alternatively, an approximate solution has to be used, in which the trajectory density is marginalised to only contain the most recent time steps.

The pmbm trackers were implemented using Bernoulli state densities and Poisson intensities of the form (28), with information form densities (31) for the state sequences. To limit the number of data associations in each update (cf. Theorem 3), analogously to the pmbm filter, see [16, Sec V.C.3], the $K$ best global hypotheses are found using Murty’s algorithm [18]. Given a Bernoulli state density, an estimate of the trajectory is obtained by selecting the most probable mixture component $j^{\star}=\arg\max_{j}w_{k|k^{\prime}}^{j}$, and computing the expected value of the state sequence $\hat{x}_{\hat{\beta}:\hat{\varepsilon}}=(Y_{k|k^{\prime}}^{j^{\star}})^{-1}y_{k|k^{\prime}}^{j^{\star}}$, where $\hat{\beta}=b_{k|k^{\prime}}^{j^{\star}}$ and $\hat{\varepsilon}=e_{k|k^{\prime}}^{j^{\star}}$. Note that in theory, the weights (26b) and (27d), and the estimate $\hat{x}_{\hat{\beta}:\hat{\varepsilon}}$, involves the inverse of the information matrix $Y_{k|k^{\prime}}^{j^{\star}}$. However, it is not necessary to compute the inverse in practice. Instead, multiplications with the inverse information matrix are solved efficiently as a sparse least squares problem. Further discussion about computationally efficient data association, as well as state recovery (both full and partial), can be found in, e.g., [17].

VII Results

In this section we present results from two simulated scenarios, one with a high number of targets and long trajectories, and one in which there is coalescence, i.e., the targets are all very close at one point. Both scenarios were generated with a 2D constant velocity motion model, see, e.g., [19, Sec. 3], with acceleration standard deviation $\sigma_{v}$, linear position measurements with covariance $R$ The pmbm tracker for all targets is compared to the $\delta$-glmb filter [11], in both the birth density was assumed known.

VII-A Measuring performance

To evaluate tracking performance, the trajectory metric $d(\cdot,\cdot)$ [20] was used with location error cut-off $c=100$, order $p=1$, and switch cost $\gamma=20$; the trajectory metric can be decomposed into a location cost, a missed target cost, a false target cost, and a switch cost. In the simulated scenarios, we apply the metric at each time step, and normalise it by the time step. This allows a comparison of how the metric evolves over time in the scenario, as opposed to only computing the metric at the final time step.

For scenarios with many targets, computing the metric becomes costly; in this case we have used the gospa metric [21] instead, which decomposes into location cost, missed target cost, and false target cost. The gospa metric is applied to the set of target states at each time step, and a trajectory metric is obtained by summing the gospa metric over time.

VII-B Scenario with many targets

A scenario with $200$ time steps and $117$ true trajectories was generated and then processed by the tracker for all trajectories. Average gospa over $100$ Monte Carlo simulations, the location cost per target per time is $0.25$, the number of missed targets per time step is $5\times 10^{-4}$, and the number of false targets per time step is $0.33$. Processing the full scenario took, on averge, $154$ seconds, i.e., $0.77$ seconds per time step.666Matlab implementation on laptop with 3.1GHz processor, 16GB memory. This shows that the pmbm trackers are accurate in terms of tracking performance, with low location error and few cardinality errors (missed and false targets), and that they are computationally feasible even for a high number of long trajectories.

VII-C Scenario with coalescence

As noted in, e.g., [5], a large number of spatially separated targets is not necessarily the most difficult scenario. Therefore, we highlight a challenging scenario with three targets that start separated, come together at the mid point of the scenario, and then separate again. The parameters were $\sigma_{v}=0.5$, $R=\mathrm{diag}\left([100,100]\right)$, $P^{\rm S}=0.99$, $P^{\rm D}=0.98$ and $\lambda^{\mathrm{FA}}=2.5\times 10^{-8}$ in the area $[-10^{3},10^{3}]\times[-10^{3},10^{3}]$. The birth density had a single component, centered in the surveillance area, with a covariance covering the whole surveillance area.

The pmbm tracker for all trajectories is compared to the $\delta$-glmb filter, which provides estimates at each time using measurements up to and including that time. For trajectory estimates extracted at the each time step of the scenario, 100 Monte Carlo simulations gave the results shown in Figure 1. If the results are summed over time, the following is obtained,

[TABLE]

Example results are shown in Figure 2. The pmbm results is three valid trajectories; however, there are switches in the $\delta$-glmb results around $(x,y)=(-350,-50)$ and $(x,y)=(300,-50)$. The switches around $(x,y)=(-350,-50)$ are due to label ambiguity upon birth: all targets are born at the same time and it is not known which birth label belongs to which measurement. The track switches around $(x,y)=(300,-50)$ are a result of ambiguity regarding which target goes where as they separate: following the separation, multiple data association sequences, i.e., global hypotheses, are approximately equally likely. Note that both the pmbm tracker and the $\delta$-glmb filter have target switches around the mid-point of the scenario; however, it is only the $\delta$-glmb filter that sometimes produces these kinds of unrealistic switches. An important difference between the $\delta$-glmb filter and the pmbm trackers is that, at any time, the pmbm trackers always provide a valid trajectory, i.e., not one that is flipping between different hypotheses at different times.

VIII Conclusion

In this paper we have presented two pmbm trackers for the set of target trajectories, and we have established connections to the pmbm filter for the set of targets. The connections show the implicit track continuity in the pmbm filter.