Poisson multi-Bernoulli mixture filter with general target-generated measurements and arbitrary clutter

A target state is denoted by $x$ and it contains the variables that describe its current dynamics, such as position and velocity, and maybe other attributes, such as orientation and extent. The state $x\in\mathcal{X}$, where $\mathcal{X}$ is a locally compact, Hausdorff and second-countable (LCHS) space [6]. The set of targets at time $k$ is $X_{k}\in\mathcal{F}\left(\mathcal{X}\right)$, where $\mathcal{F}\left(\mathcal{X}\right)$ represents the set of finite subsets of $\mathcal{X}$.

Targets move according to the standard multi-target dynamic model [6]. Given $X_{k}$, each target $x\in X_{k}$ survives with probability $p^{S}\left(x\right)$ and moves with a transition density $g\left(\cdot\left|x\right.\right)$, or disappears with probability $1-p^{S}\left(x\right)$. New targets are born at time step $k$ according to an independent Poisson point process (PPP) with intensity $\lambda_{k}^{B}\left(\cdot\right)$.

A measurement state $z\in\mathcal{Z}$ contains a sensor measurement. The set $Z_{k}\in\mathcal{F}\left(\mathcal{Z}\right)$ of measurements at time $k$ is the union of target-generated measurements and independent clutter measurements, modelled by

• •

  Each target $x\in X_{k}$ generates a set of measurements with density $f\left(\cdot|x\right)$.

• •

  The set of clutter measurements at each time step has density $c\left(\cdot\right)$.

• •

  Given $X_{k}$, the measurements generated by different targets are independent of each other and of the clutter measurements.

It should be noted that $f\left(\cdot|x\right)$ is a general density for target-generated measurements and we can accommodate any probabilistic model for target-generated measurements. The probability that at least one measurement is generated from the target (effective probability of detection [37]) is $1-f\left(\emptyset|x\right)$. For example, the standard point target model in which a target $x$ is detected with probability $p^{D}\left(x\right)$ and generates a measurement with density $l(\cdot|x)$ [6], is obtained by

$\displaystyle f\left(Z|x\right)$

$\displaystyle=\begin{cases}1-p^{D}\left(x\right)&Z=\emptyset\\ p^{D}\left(x\right)l(z|x)&Z=\left\{z\right\}\\ 0&\left|Z\right|>1.\end{cases}$

(1)

In the standard extended target model, we can receive more than one measurement from a target. Specifically, a target $x$ is detected with probability $p^{D}\left(x\right)$ and, if detected, it generates a PPP measurement with intensity $\gamma\left(x\right)l(\cdot|x)$, where $l(\cdot|x)$ is a single-measurement density and $\gamma\left(x\right)$ is the expected number of measurements [20, 38]. It is obtained by

$\displaystyle f\left(Z|x\right)$

$\displaystyle=\begin{cases}1-p^{D}\left(x\right)+p^{D}\left(x\right)e^{-\gamma\left(x\right)}&Z=\emptyset\\ p^{D}\left(x\right)\gamma^{\left|Z\right|}\left(x\right)e^{-\gamma\left(x\right)}\prod_{z\in Z}l(z|x)&\left|Z\right|>0.\end{cases}$

(2)

We can also combine both (1) and (2) to model coexisting point and extended targets, or use any of the probabilistic extended target models in [20]. In addition, the choice of $f\left(\cdot|x\right)$ can also take into account reflectivity models and the propagation conditions in the environment [8, 39, 40, 41].

We will show in this paper that, for the dynamic and measurement models in Section II-A, the predicted and posterior densities are PMBMs. That is, given the sequence of measurements $\left(Z_{1},...,Z_{k^{\prime}}\right)$, the density $f_{k|k^{\prime}}\left(\cdot\right)$ of $X_{k}$ with $k^{\prime}\in\left\{k-1,k\right\}$ is

	$\displaystyle f_{k|k^{\prime}}\left(X_{k}\right)$	$\displaystyle=\sum_{Y\uplus W=X_{k}}f_{k|k^{\prime}}^{\mathrm{p}}\left(Y\right)f_{k|k^{\prime}}^{\mathrm{mbm}}\left(W\right),$		(3)
	$\displaystyle f_{k|k^{\prime}}^{\mathrm{p}}\left(X_{k}\right)$	$\displaystyle=e^{-\int\lambda_{k|k^{\prime}}\left(x\right)dx}\prod_{x\in X_{k}}\lambda_{k|k^{\prime}}\left(x\right),$		(4)
	$\displaystyle f_{k|k^{\prime}}^{\mathrm{mbm}}\left(X_{k}\right)$	$\displaystyle=\sum_{a\in\mathcal{A}_{k|k^{\prime}}}w_{k|k^{\prime}}^{a}\sum_{\uplus_{l=1}^{n_{k|k^{\prime}}}X^{l}=X_{k}}\prod_{i=1}^{n_{k|k^{\prime}}}f_{k|k^{\prime}}^{i,a^{i}}\left(X^{i}\right),$		(5)

where $\lambda_{k|k^{\prime}}\left(\cdot\right)$ is the intensity of the PPP $f_{k|k^{\prime}}^{\mathrm{p}}\left(\cdot\right)$, representing undetected targets, and $f_{k|k^{\prime}}^{\mathrm{mbm}}\left(\cdot\right)$ is a multi-Bernoulli mixture representing targets that have been detected at some point up to time step $k^{\prime}$. The sum in (3) is the convolution sum, which implies that the PPP and MBM are independent [6]. The symbol $\uplus$ denotes the disjoint union and the sum is taken over all mutually disjoint (and possibly empty) sets $Y$ and $W$ whose union is $X_{k}$.

The MBM in (5) has $n_{k|k^{\prime}}$ Bernoulli components (potential targets), each with $h_{k|k^{\prime}}^{i}$ local hypotheses. There is a local hypothesis $a^{i}\in\left\{1,...,h_{k|k^{\prime}}^{i}\right\}$ for each Bernoulli $i>0$. To handle arbitrary clutter, we also introduce a local hypothesis $a^{0}\in\left\{1,...,h_{k|k^{\prime}}^{0}\right\}$ for clutter. Each $a^{i}$ is an index that associates the $i$-th Bernoulli, for $i\geq 0$, or the clutter, for $i=0$, to a specific sequence of subsets of the measurement set, see Section II-C. A global hypothesis is denoted by $a=\left(a^{0},a^{1},...,a^{n_{k|k^{\prime}}}\right)\in\mathcal{A}_{k|k^{\prime}}$, where $\mathcal{A}_{k|k^{\prime}}$ is the set of global hypotheses, see Section II-C. The weight of global hypothesis $a$ is $w_{k|k^{\prime}}^{a}$ and meets

$$w_{k|k^{\prime}}^{a}\propto\prod_{i=0}^{n_{k|k^{\prime}}}w_{k|k^{\prime}}^{i,a^{i}}$$

(6)

where $w_{k|k^{\prime}}^{i,a^{i}}$ is the weight of the $i$-th Bernoulli component, or clutter if $i=0$, with local hypothesis $a^{i}$, and $\sum_{a\in\mathcal{A}_{k|k^{\prime}}}w_{k|k^{\prime}}^{a}=1$. A difference with PMBM filters with PPP clutter [16, 18, 19] is that global hypotheses and weights explicitly consider clutter, with index $i=0$.

The $i$-th Bernoulli component with local hypothesis $a^{i}$ has a density

$\displaystyle f_{k|k^{\prime}}^{i,a^{i}}\left(X\right)$

$\displaystyle=\begin{cases}1-r_{k|k^{\prime}}^{i,a^{i}}&X=\emptyset\\ r_{k|k^{\prime}}^{i,a^{i}}p_{k|k^{\prime}}^{i,a^{i}}\left(x\right)&X=\left\{x\right\}\\ 0&\mathrm{otherwise}\end{cases}$

(7)

where $r_{k|k^{\prime}}^{i,a^{i}}$ is the probability of existence and $p_{k|k^{\prime}}^{i,a^{i}}\left(\cdot\right)$ is the single-target density. It should be noted that in this paper we use the following nomenclature for Bernoulli components, densities and local hypotheses: Each Bernoulli component, which is indexed by $i$ and is initiated by a non-empty subset of measurements at a given time step (see Section III), has $h_{k|k^{\prime}}^{i}$ local hypotheses, each with an associated Bernoulli density, indexed by $i,a^{i}$. Then, the total number of Bernoulli densities in (5) across all global hypotheses is $\sum_{i=1}^{n_{k|k^{\prime}}}h_{k|k^{\prime}}^{i}$, and the number of multi-Bernoulli densities is $|\mathcal{A}_{k|k^{\prime}}|$, which denotes the cardinality of $\mathcal{A}_{k|k^{\prime}}$.

We proceed to describe the set $\mathcal{A}_{k|k^{\prime}}$ of global hypotheses. We denote the measurement set at time step $k$ as $Z_{k}=\left\{z_{k}^{1},...,z_{k}^{m_{k}}\right\}$. We refer to measurement $z_{k}^{j}$ using the pair $\left(k,j\right)$ and the set of all such measurement pairs up to and including time step $k$ is denoted by $\mathcal{M}_{k}$. Then, a local hypothesis $a^{i}$ for $i=\left\{0,1,...,n_{k|k^{\prime}}\right\}$ has an associated set of measurement pairs denoted as $\mathcal{M}_{k}^{i,a^{i}}\subseteq\mathcal{M}_{k}$. The set $\mathcal{A}_{k|k^{\prime}}$ of all global hypotheses meets

	$\displaystyle\mathcal{A}_{k|k^{\prime}}=$	$\displaystyle\left\{\left(a^{0},a^{1},...,a^{n_{k|k^{\prime}}}\right):a^{i}\in\left\{1,...,h_{k|k^{\prime}}^{i}\right\}\,\forall i,\right.$
		$\displaystyle\left.\bigcup_{i=0}^{n_{k|k^{\prime}}}\mathcal{M}_{k^{\prime}}^{i,a^{i}}=\mathcal{M}_{k^{\prime}},\mathcal{M}_{k^{\prime}}^{i,a^{i}}\cap\mathcal{M}_{k^{\prime}}^{j,a^{j}}=\emptyset,\,\forall i\neq j\right\}.$

That is, each measurement must be assigned to a local hypothesis, and there cannot be more than one local hypothesis with the same measurement. In this paper, we construct the (trees of) local hypotheses recursively, and we allow for more than one measurement to be associated to the same local hypothesis at the same time step, i.e., $\mathcal{M}_{k}^{i,a^{i}}$ may contain zero, one or more measurements. At time step zero, the filter is initiated with $n_{0|0}=0$, $w_{0|0}^{0,1}=1$, $h_{0|0}^{0}=1$, and $\mathcal{M}_{0}^{0,1}=\emptyset$.

Given two real-valued functions $a\left(\cdot\right)$ and $b\left(\cdot\right)$ on the target space, we denote their inner product as

$\displaystyle\left\langle a,b\right\rangle$

$\displaystyle=\int a\left(x\right)b\left(x\right)dx.$

(8)

Then, the update of the predicted PMBM $f_{k|k-1}\left(\cdot\right)$ with $Z_{k}$ is given in the following theorem.

Theorem 1 is proved in Appendix A.

We can see that the updated PPP intensity in (9) is the predicted intensity multiplied by the probability of receiving no measurements, as in the PMBM filters for PPP clutter [16, 17, 18, 19]. Contrary to PMBM filters with PPP clutter, the general PMBM filter extends local and global hypotheses to explicitly consider clutter data associations, whose weights depend on the clutter density via (12). As data associations for clutter are separated from the target hypotheses, the probability of existence of new Bernoulli components $r_{k|k}^{i,2}$ is equal to one, as in the PMBM filter for unknown clutter rate in [42]. This is different from the PMBM filters for PPP clutter, where, for each separate measurement, the hypotheses for clutter and the first detection of a new target are merged into one. After merging, the probability of existence of a new Bernoulli component depends on the clutter intensity, and the number of global hypotheses for PPP clutter is lower, see Section II-C.

We also observe that each non-empty subset of $Z_{k}$ generates a new Bernoulli, which has two local hypotheses. It is also possible to instead reformulate the update to only generate $m_{k}$ new Bernoulli components by increasing the number of local hypotheses of the new Bernoulli components [43, Sec. IV]. We proceed to illustrate the differences between the hypotheses generated by this PMBM update, and the update with PPP clutter with an example.

In this subsection, we discuss the spooky action at a distance that may appear when we deal with arbitrary clutter, and an alternative choice of hypotheses when the clutter density has some additional structure.

Theorem 1 also holds for a predicted density of the form MBM or MBM₀₁ by setting $\lambda_{k|k-1}\left(\cdot\right)=0$. Therefore, we obtain the MBM filter by applying the standard MBM prediction step [17, Sec. III.E], which includes multi-Bernoulli birth model, followed by Theorem 1 update. The MBM₀₁ filter is the same as the MBM filter with the difference that the MBM₀₁ filter prediction step expands each Bernoulli density into two hypotheses with deterministic target existence [17, Eq. (36)]. This results in an unnecessary exponential increase in the number of global hypotheses in the prediction step, so it is not recommended.

It is also relevant to notice that the general PMBM update in Theorem 1 also includes Bernoulli densities with deterministic existence (when updating an existing Bernoulli density with a non-empty measurement set or initiating a new Bernoulli component). Nevertheless, these Bernoulli densities with deterministic existence are required to deal with arbitrary clutter. In contrast, the MBM₀₁ filter creates the deterministic Bernoulli densities in the prediction, which is not necessary.

Both MBM and MBM₀₁ can be directly extended to include deterministic, fixed labels assigned to each Bernoulli density of the birth model without changing the filtering recursion [17, Sec. IV.V]. This procedure yields the labelled MBM and labelled MBM₀₁ filters. The $\delta$-GLMB filtering recursion is analogous to the labelled MBM₀₁ filtering recursion, so this procedure provides the $\delta$-GLMB filter for general target-generated measurements and arbitrary clutter.

PMB (and multi-Bernoulli) filters for arbitrary clutter can be obtained by projecting the updated PMBM into a PMB using KLD minimisation [16, 33, 45, 43].

The equations in Theorem 1 are also valid for point targets with arbitrary clutter. In this setting, the target-generated measurement density is given by (1). With the definitions of local and global hypotheses in Section II-C, many hypotheses associate more than one measurement to the same Bernoulli component at the same time. However, according to (1), the weights of these hypotheses will be zero. Therefore, we can directly exclude these hypotheses from the set $\mathcal{A}_{k|k}$ of global hypotheses, by restricting the cardinality of $\mathcal{M}_{k}^{i,a^{i}}$ to zero or one for $i>0$. This approach, combined with the corresponding prediction step, provides the PMBM, MBM, MBM₀₁ ($\delta$-GLMB) filtering recursions for point targets and arbitrary clutter.

In this section, we calculate the number of global hypotheses of PMBM filters to compare: 1) PPP clutter versus arbitrary clutter, and 2) point targets versus the general target-generated measurement model. This analysis helps us understand the differences in the hypothesis structures of the PMBM filters and provides a measure of the complexity of the corresponding data association problems. To simplify the analysis, we first consider the update of a PPP prior in Section III-E1, and then the update of a PMB prior in Section III-E2.

To simplify notation in this section, we drop time indices in the number of measurements $m$ and the measurement set $\left\{z^{1},...,z^{m}\right\}$. In addition, the Bernoulli densities for the predicted global hypothesis, representing a PMB, have parameters $\left(r^{i},p^{i}(\cdot)\right)$ for $i\in\left\{1,...,n_{k|k-1}\right\}$.

Then, we recall that for the point target model, after the update, we have $n_{k|k}=n_{k|k-1}+m$ Bernoulli components, which includes $n_{k|k-1}$ predicted Bernoulli components and the new $m$ Bernoulli components, one generated by each measurement. The data associations can be represented as a vector $\gamma_{1:m}=\left(\gamma_{1},...,\gamma_{m}\right)$ where $\gamma_{j}\in\{0,1,...,n_{k|k}\}$ is the Bernoulli component index associated to the $j$-th measurement and $\gamma_{i}\neq\gamma_{j}$ for all $i\neq j$ such that $\gamma_{i}>0$ and $\gamma_{j}>0$. If $\gamma_{j}=0$, it means that the $j$-th measurement is clutter. The set of vectors $\gamma_{1:m}$ that meet this property is denoted by $\Gamma$.

We also use $1_{\Gamma}\left(\cdot\right)$ to denote the indicator function of set $\Gamma$: $1_{\Gamma}\left(\gamma_{1:m}\right)=1$ if $\gamma_{1:m}\in\Gamma$, and $1_{\Gamma}\left(\gamma_{1:m}\right)=0$ otherwise.

The probability of data association $\gamma_{1:m}$, without accounting for the weight of the predicted global hypothesis, is

$\displaystyle p\left(\gamma_{1:m}\right)$

$\displaystyle\propto 1_{\Gamma}\left(\gamma_{1:m}\right)c\left(Z_{c}\left(\gamma_{1:m}\right)\right)\prod_{j=1:\gamma_{j}>0}^{m}\eta^{j}\left(\gamma_{j}\right)$

(46)

where the measurement set corresponding to clutter is

$\displaystyle Z_{c}\left(\gamma_{1:m}\right)$

$\displaystyle=\left\{z^{j}:\gamma_{j}=0\right\}$

(47)

and

	$\displaystyle\eta^{j}\left(\gamma_{j}\right)$
	$\displaystyle\;=\begin{cases}\frac{r^{\gamma_{j}}\int p^{D}\left(x\right)l\left(z^{j}|x\right)p^{\gamma_{j}}(x)dx}{1-r^{\gamma_{j}}\left\langle p^{\gamma_{j}},p^{D}\right\rangle}&\gamma_{j}\leq n_{k|k-1}\\ \int p^{D}\left(x\right)l\left(z^{j}|x\right)\lambda_{k|k-1}\left(x\right)dx&n_{k|k-1}<\gamma_{j}\leq n_{k|k}.\end{cases}$		(48)

In (48), the second line corresponds to the weight of new Bernoulli components, see (24). The first line corresponds to the weight of the predicted Bernoulli components with a detection $z^{j}$, see (20), divided by the misdetection weight, see (15), both using the point target measurement model (1).

To obtain samples from density (46) using Gibbs sampling, we calculate the conditional density [48]

	$\displaystyle p\left(\gamma_{q}|\gamma_{1:q-1},\gamma_{q+1:m}\right)$	$\displaystyle\propto p\left(\gamma_{1:m}\right)$		(49)
		$\displaystyle\propto 1_{\Gamma}\left(\gamma_{1:m}\right)c\left(Z_{c}\left(\gamma_{1:m}\right)\right)$
		$\displaystyle\quad\times\prod_{j=1:\gamma_{j}>0}^{m}\eta^{j}\left(\gamma_{j}\right)$		(50)

where $q\in\{1,...,m\}$. As $\gamma_{1:q-1},\gamma_{q+1:m}$ are constants in (50), $l^{j}\left(\gamma_{j}\right)$ for $j\neq q$ are constants too, and we can write

	$\displaystyle p\left(\gamma_{q}|\gamma_{1:q-1},\gamma_{q+1:m}\right)$
	$\displaystyle\propto\begin{cases}1_{\Gamma}\left(\gamma_{1:m}\right)c\left(Z_{c}\left(\gamma_{1:m}\right)\right)\eta^{q}\left(\gamma_{q}\right)&\gamma_{q}>0\\ c\left(Z_{c}\left(\gamma_{1:m}\right)\right)&\gamma_{q}=0\end{cases}$		(51)

where we have used that $\left(\gamma_{1:q-1},0,\gamma_{q+1:m}\right)\in\Gamma$. Equation (51) can be evaluated for all possible values $\gamma_{q}\in\{0,1,...,m\}$, and therefore, we can sample from this categorical distribution. To perform Gibbs sampling, we can start with $\gamma_{1:m}$ having all entries equal to zero, and sample $\gamma_{q}$ from $q=1$ to $m$ for a total number of $K$ times. Then, as in [21], we remove the repeated samples of data association hypotheses to yield the newly generated data association hypotheses. In the implementations, we set $K=\left\lceil N_{h}\cdot w_{j}\right\rceil$, where $N_{h}$ is the maximum number of global hypotheses.

In this subsection, we specify (51) for the case where the clutter $c\left(\cdot\right)$ follows an IID cluster process [6] whose the single-measurement density $\breve{c}\left(\cdot\right)$ is uniform in the field-of view $A$. That is,

	$\displaystyle c\left(Z\right)$	$\displaystyle=|Z|!\rho_{c}\left(|Z|\right)\prod_{z\in Z}\breve{c}\left(z\right)$		(52)
	$\displaystyle\breve{c}\left(z\right)$	$\displaystyle=u_{A}\left(z\right)$		(53)

where $\rho_{c}\left(|Z|\right)$ is the cardinality distribution and $u_{A}\left(\cdot\right)$ is a uniform density in $A$.

Let $m_{c}$ denote the number of clutter measurements in $\gamma_{1:q-1},\gamma_{q+1:m}$, i.e.,

$\displaystyle m_{c}$

$\displaystyle=\sum_{j=1:j\neq q}^{m}\delta_{0}\left[\gamma_{j}\right]$

(54)

where $\delta_{j}[\cdot]$ is the Kronecker delta at $j$. Then, considering that all measurements are in the field-of-view, (51) becomes

	$\displaystyle p\left(\gamma_{q}|\gamma_{1:q},\gamma_{q+1:m}\right)$
	$\displaystyle\propto\begin{cases}1_{\Gamma}\left(\gamma_{1:m}\right)\eta^{q}\left(\gamma_{q}\right)\rho_{c}\left(m_{c}\right)m_{c}!\frac{1}{|A|^{m_{c}}}&\gamma_{q}>0\\ \rho_{c}\left(m_{c}+1\right)\left(m_{c}+1\right)!\frac{1}{|A|^{m_{c}+1}}&\gamma_{q}=0\end{cases}$		(55)

where $|A|$ is the area (or the hypervolume) of $A$.

In this section, we consider point targets with clutter being an IID cluster process with negative binomial (NB) cardinality distribution [49]. We compare the following filters: the standard PMBM and PMB filter (with PPP clutter assumption), and the arbitrary clutter PMBM and PMB filters, referred to as A-PMBM and A-PMB. The rest of the parameters are [17]: maximum number of global hypotheses $N_{h}=5,000$, threshold for MBM pruning $10^{-4}$, threshold for pruning the PPP weights $\Gamma_{p}=10^{-5}$, threshold for pruning Bernoulli densities $\Gamma_{b}=10^{-5}$, and ellipsoidal gating with threshold 20. These filters use Estimator 3 as it provides the best performance among the standard PMBM estimators [17, Sec. VI] in this scenario. Estimator 3 obtains the global hypothesis with a deterministic cardinality with highest weight and then reports the mean of the targets in this hypothesis.

Additionally, we have implemented two filters with multi-Bernoulli birth and PPP clutter: the MBM filter [17, 50], and the $\delta$-GLMB filter with joint prediction and update [21]. The $\delta$-GLMB filter has been implemented with 7,000 global hypotheses and pruning threshold $10^{-15}$. All the previous filters have been implemented using Gibbs sampling to select relevant global hypotheses in the update of each predicted global hypothesis. We have also implemented a Gaussian mixture PHD filter for PPP clutter [6, 51] and a Gaussian mixture PHD filter for NB IID clutter in [27], which we refer to as the NB-PHD filter. The PHD filters have parameters: pruning threshold $10^{-5}$, merging threshold 0.1, and maximum number of Gaussian components 30.