Second order multi-object filtering with target interaction using determinantal point processes

TL;DR

This paper introduces a second-order PHD filter based on Determinantal Point Processes to model target interactions, specifically repulsion, improving multi-target tracking accuracy over traditional independent-target assumptions.

Contribution

It develops a novel second-order PHD filter using DPPs that captures target interactions by propagating variance and covariance, extending beyond first-order methods.

Findings

Successfully models target repulsion with DPPs

Provides Monte Carlo simulations demonstrating improved tracking

Estimates correlation between targets in measurement domains

Abstract

The Probability Hypothesis Density (PHD) filter, which is used for multi-target tracking based on sensor measurements, relies on the propagation of the first-order moment, or intensity function, of a point process. This algorithm assumes that targets behave independently, an hypothesis which may not hold in practice due to potential target interactions. In this paper, we construct a second-order PHD filter based on Determinantal Point Processes (DPPs) which are able to model repulsion between targets. Such processes are characterized by their first and second-order moments, which allows the algorithm to propagate variance and covariance information in addition to first-order target count estimates. Our approach relies on posterior moment formulas for the estimation of a general hidden point process after a thinning operation and a superposition with a Poisson Point Process (PPP), and on suitable approximation formulas in the determinantal point process setting. The repulsive properties of determinantal point processes apply to the modeling of negative correlation between distinct measurement domains. Monte Carlo simulations with correlation estimates are provided.

Figures

Equations

$${\mathbf{N}}_{\sigma }:=\{\xi \subseteq {\mathbb{R}}^d:|\xi \cap A|<\infty \text{ for all relatively compact sets }A\subset {\mathbb{R}}^d\}$$

$$\xi \mapsto ⟨f,\xi ⟩:=\sum _{y\in \xi }f(y)$$

$${\mathcal{L}}_{\Phi }(f)={\int }_{{\mathbf{N}}_{\sigma }}{e}^{-⟨f,\xi ⟩}\mathbb{P}(d\xi ).$$

$$\mathrm{IE}[F(\Phi )]:={\int }_{{\mathbf{N}}_{\sigma }}F(\xi )\mathbb{P}(d\xi ).$$

$$\mathrm{IE}\left[F(\Phi )\right]=F(\varnothing )j_{\Lambda }^{(0)}+\sum _{n\ge 1}\frac{1}{n!}{\int }_{{\Lambda }^n}F(\{x_1,\dots ,x_n\})j_{\Lambda }^n\left(x_1,\dots ,x_n\right)\nu (d{x}_1)\cdots \nu (d{x}_n),$$

$$\mathrm{IE}\left[\prod _{i=1}^k\Phi (B_i)\right]={\int }_{B_1\times \cdots \times B_k}{\rho }_{\Phi }^{(k)}(x_1,\dots ,x_k)\nu (d{x}_1)\cdots \nu (d{x}_k),$$

$$\mathrm{IE}\left[\prod _{i=1}^n\frac{\Phi (B_i)!}{(\Phi (B_i)-k_i)!}\right]={\int }_{B_1^{k_1}\times \cdots \times B_n^{k_n}}{\rho }_{\Phi }^{(N)}(x_1,\dots ,x_N)\nu (d{x}_1)\cdots \nu (d{x}_N).$$

$${\int }_{B_1\times \cdots \times B_n}{\rho }_{\Phi }^{(n)}(x_1,\dots ,x_n)\nu (d{x}_1)\cdots \nu (d{x}_n)=\mathrm{IE}\left[\sum _{x_1,\dots ,x_n\in \Phi }{1}_{B_1}(x_1)\cdots 1_{B_n}(x_n)\right],$$

$$j_{\Lambda }^{(n)}(x_1,\dots ,x_n)=\sum _{m\ge 0}\frac{(-1{)}^m}{m!}{\int }_{{\Lambda }^m}{\rho }_{\Phi }^{(n+m)}(x_1,\dots ,x_n,y_1,\dots ,y_m)\nu (d{y}_1)\cdots \nu (d{y}_m),$$

$${\rho }_{\Phi }^{(n)}(x_1,\dots ,x_n)=\sum _{m\ge 0}\frac{1}{m!}{\int }_{{\Lambda }^m}{j}_{\Lambda }^{(m+n)}(x_1,\dots ,x_n,y_1,\dots ,y_m)\nu (d{y}_1)\cdots \nu (d{y}_m),$$

$$h\mapsto {\mathcal{G}}_{\Phi }(h)$$

$$\frac{{\partial }_g}{\partial h}\mathcal{F}(h):=\underset{\epsilon \to 0}{\lim }\frac{\mathcal{F}(h+\epsilon g)-\mathcal{F}(h)}{\epsilon }.$$

$$\frac{{\partial }_{{\delta }_x}}{\partial h}\mathcal{F}(h):=\underset{n\to \infty }{\lim }\frac{{\partial }_{g_n}}{\partial h}\mathcal{F}(h),$$

$$j_{\Phi }^{(n)}(x_1,\dots ,x_n)=\frac{{\partial }_{{\delta }_{x_1}}}{\partial h}\cdots \frac{{\partial }_{{\delta }_{x_n}}}{\partial h}{\mathcal{G}}_{\Phi }(h{)}_{\mid h=0},x_1,\dots ,x_n\in \Lambda ,$$

$${\rho }_{\Phi }^{(n)}(x_1,\dots ,x_n)=\frac{{\partial }_{{\delta }_{x_1}}}{\partial h}\cdots \frac{{\partial }_{{\delta }_{x_n}}}{\partial h}{\mathcal{G}}_{\Phi }(h{)}_{\mid h=1},x_1,\dots ,x_n\in \Lambda ,$$

$$g\mapsto {\mathcal{G}}_{{\Xi }_s}(g\mid x)$$

$$q_d(x):=1-p_d(x)\text{mbox}and{\tilde{l}}_d(z|x):=p_d(x)l_d(z|x),x\in \Lambda .$$

$$x_{1:n}=(x_1,\dots ,x_n)\in {\Lambda }^n,\text{mbox}and\nu (d{x}_{1:n}):=\nu (d{x}_1)\cdots \nu (d{x}_n),$$

$$(h,g)\mapsto{\cal G}_{\Phi,\Xi}(h,g):={\cal G}_{\Xi_{c}}(g){\cal G}_{\Phi}\big{(}h(\cdot){\cal G}_{\Xi_{s}}(g\mid\cdot\ \!)\big{)},$$

$$\displaystyle{\cal G}_{\Phi}\big{(}h(\cdot){\cal G}_{\Xi_{s}}(g\mid\cdot\ \!)\big{)}$$

$$=\sum _{n\ge 0}\frac{1}{n!}{\int }_{{\Lambda }^n}{j}_{\Phi }^{(n)}(x_{1:n})\prod _{j=1}^{n}h(x_j)\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)q_d^{n-k}\prod _{i=1}^k{\int }_{\Lambda }g(z){\tilde{l}}_d(z|x_i)\nu (dz)\nu (d{x}_{1:n})$$

$$=\sum _{k\ge 0}{\int }_{{\Lambda }^k}\prod _{i=1}^k{\int }_{\Lambda }g(z){\tilde{l}}_d(z|x_i)\nu (dz)\sum _{n\ge 0}\frac{q_d^n}{n!k!}{\int }_{{\Lambda }^n}{j}_{\Phi }^{(k+n)}(x_{1:k+n})\prod _{j=1}^{k+n}h(x_j)\nu (d{x}_{1:k+n}).$$

$$g\mapsto {\mathcal{G}}_{\Xi }(g)$$

$$\frac{{\partial }_{{\delta }_x}}{\partial g}{\mathcal{G}}_{\Xi }(g)$$

$${\mu }_{\Xi }^{(1)}(x)=\frac{{\partial }_{{\delta }_x}}{\partial g}{\mathcal{G}}_{\Xi }(g{)}_{\mid g=1}={\mu }_{{\Xi }_c}^{(1)}(x)+{\int }_{\Lambda }{\tilde{l}}_d(x|y){\mu }_{\Phi }^{(1)}(y)\nu (dy).$$

$$\displaystyle{\frac{\partial_{\delta_{x}}}{\partial g}\frac{\partial_{\delta_{y}}}{\partial g}\bigg{(}\mathcal{G}_{\Phi}(\mathcal{G}_{\Xi_{s}}(g\mid\cdot))\mathcal{G}_{\Xi_{c}}(g)\bigg{)}}$$

$$+\frac{{\partial }_{{\delta }_y}}{\partial g}{\mathcal{G}}_{\Phi }({\mathcal{G}}_{{\Xi }_s}(g\mid \cdot ))\frac{{\partial }_{{\delta }_x}}{\partial g}{\mathcal{G}}_{{\Xi }_c}(g)+{\mathcal{G}}_{\Phi }({\mathcal{G}}_{{\Xi }_s}(g\mid \cdot ))\frac{{\partial }_{{\delta }_x}}{\partial g}\frac{{\partial }_{{\delta }_y}}{\partial g}({\mathcal{G}}_{{\Xi }_c}(g)),$$

$${\rho }_{\Xi }^{(2)}(x,y)$$

$$j_{\Phi \mid \Xi =z_{1:m}}^{(n)}(x_1,\dots ,x_n)=\frac{j_{\Phi ,\Xi =z_{1:m}}^{(n,m)}(x_1,\dots ,x_n)}{j_{\Xi }^{(m)}(z_{1:m})},x_1,\dots ,x_n\in \Lambda ,$$

$$j_{\Phi ,\Xi =z_{1:m}}^{(n,m)}(x_1,\dots ,x_n)=j_{\Phi }^{(n)}(x_1,\dots ,x_n)\sum _{\genfrac{}{}{0ex}{}{S\subset \{1,\dots ,m\}}{|S|\le n}}\frac{n!q_d^{n-|S|}}{(n-|S|)!}\prod _{j\notin S}{l}_c(z_j)\sum _{\pi :S\to \{1,\dots ,n\}}\prod _{i\in S}{\tilde{l}}_d(z_i|x_{\pi (i)}),$$

$$j_{\Xi }^{(m)}(z_1,\dots ,z_m)=\sum _{n\ge 0}\sum _{\genfrac{}{}{0ex}{}{S\subset \{1,\dots ,m\}}{|S|\le n}}\prod _{j\notin S}{l}_c(z_j)\frac{q_d^{n-|S|}}{(n-|S|)!}{\int }_{{\Lambda }^n}{j}_{\Phi }^{(n)}(x_{1:n})\prod _{i\in S}{\tilde{l}}_d(z_i|x_i)\nu (d{x}_{1:n}),$$

$$\frac{{\partial }_{{\eta }_1}}{\partial h}\cdots \frac{{\partial }_{{\eta }_n}}{\partial h}\frac{{\partial }_{f_1}}{\partial g}\cdots \frac{{\partial }_{f_m}}{\partial g}{\mathcal{G}}_{\Phi ,\Xi }(h,g)$$

Full text

Keywords: Probability hypothesis density (PHD) filter; higher-order statistics; correlation; second-order moment; determinantal point processes; multi-object filtering; multi-target tracking.

Mathematics Subject Classification (2010): 60G35; 60G55; 62M30; 62L12.

1 Introduction

Probability Hypothesis Density (PHD) filters have been introduced in Mahler (2003) for multi-target tracking in cluttered environments. The construction of the prediction point process $\Phi$ therein uses multiplicative point processes, see e.g. Moyal (1962), Moyal (1964), by thinning and shifting a prior point process $\Psi$, and superposition with a birth point process. The posterior point process $\Phi|\Xi$ is obtained by conditioning $\Phi$ given a measurement point process of targets $\Xi$, also constructed by thinning, shifting and superposition. This step relies on Bayesian estimation with a Poisson point process prior, see e.g. van Lieshout (1995), Mori (1997), Portenko et al. (1997). PHD filters have low complexity, and they allow for explicit update formulas see e.g. Clark et al. (2016) for a review.

While the PHD filter of Mahler (2003) is based on Poisson point processes, several extensions of the PHD filter to non-Poisson prior distributions have been proposed. Cardinalized Probability Hypothesis Density (CPHD) filters have been introduced in Mahler (2007) as a generalization in which the target count is allowed to have an arbitrary distribution. In de Melo and Maskell (2019), discretized Gamma distributions are used to design an efficient approximation of the CPHD filter cardinality distribution. Other generalizations include the Gauss-Poisson point processes that generalize the Poisson point process by allowing for two-point clusters, and have been used in Singh et al. (2009). The PHD filter has been implemented using the Sequential Monte Carlo (SMC) method in Vo et al. (2005), and using Gaussian mixtures in Vo and Ma (2006).

PHD filters approximate the distribution of the number of targets by a Poisson distribution estimated by a single mean (or variance) parameter, which can result into high variance estimates when the estimated mean is high. Second-order PHD filters that can propagate distinct information on mean and variance parameters have been recently proposed in Schlangen et al. (2018), based on the Panjer point process defined therein, where the Panjer cardinality distribution encompasses the binomial, Poisson and negative binomial distributions. Other multi-target filters propagating second-order moment information have also been recently proposed, see e.g. Clark and de Melo (2018) for a filter that propagates second-order point process factorial cumulants.

A common feature of cardinalized filters is to assume that target locations are distributed as $n$ independent random samples according to a reference intensity measure, given that the observation window contains $n$ points. While this hypothesis is natural and facilitates an explicit derivation of prediction formulas, it does not reflect potential interaction between targets. In addition, as observed in the simulations of Section 6, the presence of repulsion between targets can degrade the performance of the Poisson PHD filter.

As a response, we propose to construct a PHD filter based on determinantal point processes introduced in Macchi (1975), which are able to model repulsion among configuration points on a target domain $\Lambda\subset{\mathord{\mathbb{R}}}^{d}$, see also Soshnikov (2000) and Shirai and Takahashi (2003). Taking into account correlation via more general point process-based PHD filters poses several challenges linked to the derivation of closed form filtering formulas. In addition, the distribution of general point processes relies on Janossy densities which may not be characterized by the knowledge of a finite number of moments. Determinantal point processes, on the other hand, are characterized by their kernel functions $(K(x,y))_{(x,y)\in\Lambda^{2}}$, and their Janossy densities can be recovered from their first and second-order moments. In this setting, the knowledge of first and second-order moments can be used to update the Janossy densities that characterize the underlying determinantal point process.

Discrete determinantal point processes have also been recently used for the pruning of Gaussian components in the Gaussian Mixture (GM) PHD filter in Jorquera et al. (2019), see also Jorquera et al. (2018) for other applications to multi-target tracking. Permanental processes have been used in Mahler (2015) to propagate a joint Poisson distribution, however this approach is distinct from the determinantal setting. See Koch (2018) for the use the exclusion principle in multi-target tracking by a fermionic filtering update of anti-symmetric components in the joint probability density functions of states. Note also that determinantal point processes have been originally introduced in Macchi (1975) to represent configurations of fermions.

After recalling general facts and notation on point processes in Section 2, we derive general formulas for the distribution, and for the first and second-order moments, of a posterior point process $\Phi|\Xi$ in Section 3, see also Lund and Rudemo (2000). In Section 4 we review the construction of determinantal point processes, and in Section 5 we present a second-order PHD filtering algorithm based on determinantal point processes, with the computation of the prediction kernel $K_{\Phi}(x,y)$ and of the updated kernel $K_{\Phi\mid\Xi}(x,y)$. An implementation of the Poisson PHD filter that allows for performance evaluation using measurement-estimate associations is presented in Section 6 with numerical illustrations. This simulation is based on the sequential Monte Carlo (or particle filtering) method with a nearly constant turn-rate motion dynamics, see Vo et al. (2009), Li et al. (2017), to which we add a repulsion term. In Section 7 we implement the determinantal PHD filter using the sequential Monte Carlo method. The implementation of the algorithm relies on closed-form filter update expressions obtained from approximation formulas for corrector terms and Janossy densities presented in appendix.

2 Preliminaries on point processes

In this section we review the properties of point processes; see, e.g. Daley and Vere-Jones (2003), Decreusefond et al. (2016), and references therein. For any subset $A\subseteq{\mathord{\mathbb{R}}}^{d}$, let $|A|$ denote the cardinality of $A$, setting $|A|=\infty$ if $A$ is not finite, and let

[TABLE]

denote the set of locally finite point configurations on ${\mathord{\mathbb{R}}}^{d}$, which is identified with the set of all nonnegative integer-valued Radon measures $\xi$ on ${\mathord{\mathbb{R}}}^{d}$ such that $\xi(\{x\})\in\{0,1\}$ for all $x\in{\mathord{\mathbb{R}}}^{d}$. We denote by ${\mathscr{N}_{\sigma}}$ the Borel $\sigma$-field generated by the weakest topology that makes the mappings

[TABLE]

continuous for all continuous and compactly supported functions $f$ on ${\mathord{\mathbb{R}}}^{d}$. Given $\Lambda$ a relatively compact subset of ${\mathord{\mathbb{R}}}^{d}$, we let $\mathbf{N}_{\sigma}(\Lambda)$ be the space of finite configurations on $\Lambda$.

We consider a simple and locally finite point process $\Phi$ on $\Lambda$, defined as a random element on a probability space $(\Omega,{\mathscr{N}_{\sigma}})$ with values in $\mathbf{N}_{\sigma}(\Lambda)$, and denote its distribution by ${\mathbb{P}}$. The point process $\Phi$ is characterized by its Laplace transform $\mathcal{L}_{\Phi}$ which is defined, for any measurable nonnegative function $f$ on ${\Lambda}$, by

[TABLE]

We denote the expectation of an integrable random variable $F$ defined on $(\mathbf{N}_{\sigma},{\mathscr{N}_{\sigma}},{\mathbb{P}})$ by

[TABLE]

Janossy densities

For any relatively compact subset $A\subseteq{\Lambda}$, the Janossy densities of $\Phi$ w.r.t. a reference Radon measure ${\nu}$ on ${\Lambda}$ are symmetric measurable functions $j^{(n)}_{\Lambda}:\Lambda^{n}\rightarrow[0,\infty)$ satisfying

[TABLE]

for all measurable functions $F:\mathbf{N}_{\sigma}(\Lambda)\rightarrow[0,\infty)$; see, e.g., Georgii and Yoo (2005).

For $n\geq 1$, the Janossy density $j^{(n)}_{\Lambda}(x_{1},\ldots,x_{n})$ is proportional, up to a multiplicative constant, to the joint density of the $n$ points of the point process, given that it has exactly $n$ points. For $n=0$, $j^{(0)}_{\Lambda}(\emptyset)$ is the probability that there are no points in $\Lambda$.

Correlation functions

The correlation functions of $\Phi$ w.r.t. the reference measure ${\nu}$ on ${\Lambda}$ are measurable symmetric functions $\rho^{(k)}_{\Phi}:{\Lambda}^{k}\longrightarrow[0,\infty)$ such that

[TABLE]

for any family of mutually disjoint bounded subsets $B_{1},\ldots,B_{k}$ of ${\Lambda}$, $k\geq 1$. More generally, if $B_{1},\ldots,B_{n}$ are disjoint bounded Borel subsets of ${\Lambda}$ and $k_{1},\ldots,k_{n}$ are integers such that $\sum_{i=1}^{n}k_{i}=N$, we have

[TABLE]

In addition, we let $\rho^{(n)}_{\Phi}(x_{1},\ldots,x_{n})=0$ whenever $x_{i}=x_{j}$ for some $1\leq i\neq j\leq n$. In other words, the factorial moment density $\rho^{(n)}_{\Phi}(x_{1},\ldots,x_{n})$ of $\Phi$, $x_{1},\ldots,x_{n}\in\Lambda$, $x_{i}\not=x_{j}$, $1\leq i<j\leq n$, is defined from the relation

[TABLE]

for mutually disjoint measurable subsets $B_{1},\ldots,B_{n}\subset\Lambda$, where ${\bf 1}_{B_{i}}$ denotes the indicator function over $B_{i}$, $i=1,\ldots,n$. Heuristically, $\rho^{(n)}_{\Phi}(x_{1},\ldots,x_{n})\,{\nu}(dx_{1})\cdots{\nu}(dx_{n})$ represents the probability of finding a particle in the vicinity of each $x_{i}$, $i=1,\ldots,n$.

We also recall that the Janossy densities $j^{(n)}_{\Lambda}$ can be recovered from the correlation functions $\rho^{(m)}_{\Phi}$ via the relation

[TABLE]

and vice versa using the equality

[TABLE]

see Theorem 5.4.II of Daley and Vere-Jones (2003).

Probability generating functionals

The Probability Generating Functional (PGFl) of the point process $\Phi$, see Moyal (1962), is defined by

[TABLE]

for $h\in L^{\infty}(\Lambda)$ a bounded measurable function on $\Lambda$. Given ${\cal F}$ a functional on $L^{\infty}(\Lambda)$, we will use the functional derivative $\partial_{g}/\partial h$ of ${\cal F}(h)$ in the direction of $g\in L^{\infty}(\Lambda)$, defined as

[TABLE]

Given $x\in\Lambda$, we also let

[TABLE]

where $(g_{n})_{n\geq 1}$ is a sequence of bounded functions converging weakly to the Dirac distribution $\delta_{x}$ at $x\in\Lambda$.

This construction allows one to recover the Janossy densities $j^{(n)}_{\Phi}(x_{1},\ldots,x_{n})$ and factorial moment densities $\rho^{(n)}_{\Phi}(x_{1},\ldots,x_{n})$ of $\Phi$ from the PGFl ${\cal G}_{\Phi}(h)$ as

[TABLE]

see e.g. § 2.4 of Clark et al. (2016), and as

[TABLE]

with $x_{i}\not=x_{j}$, $1\leq i<j\leq n$; see, e.g., § 3.4 of Clark et al. (2016).

3 Posterior point process distribution

In this section we compute the Janossy densities, and the first and second-order moments, of a posterior point process of targets $\Phi$ given the point process $\Xi$ of sensor measurements. The case of a Poisson prior was treated in Chapter 5 of van Lieshout (1995), see Theorem 29 therein and also Mori (1997) for early sensor fusion applications, or Theorems 6.1-6.2 Portenko et al. (1997) for related derivations based on Laplace transforms.

In Propositions 3.2 and 3.3 below we will use extensions of the corrector terms $l^{(1)}_{z_{1:m}}$, $l^{(2)}_{z_{1:m}}$ introduced in Delande et al. (2014) for the cardinalized PHD filter, see Vo et al. (2007). We start with a review of the thinning and shifting of point processes; see, e.g., Clark et al. (2016) and references therein for details.

Thinning and shifting of point processes

The point process $\Xi$ of sensor measurements is constructed via the following steps.

1. (i)

Thinning and shifting. Every target point $x\in\Phi$ is kept with probability $p_{d}(x)\in(0,1]$ and shifted according to the probability density function $l_{d}(\cdot|x)$, by branching the hidden point process $\Phi$ with a Bernoulli point process $\Xi_{s}$ with PGFl

[TABLE]

where for compactness of notation we take

[TABLE] 2. (ii)

The point process $\Xi$ is obtained by superposing a Poisson point process $\Xi_{c}$ with intensity function $l_{c}(\cdot)$, representing clutter, to the above thinning and shifting of $\Phi$.

In the sequel we use the shorthand notation

[TABLE]

with $x_{1:0}=\emptyset$; see Delande et al. (2014), Schlangen et al. (2018). The joint PGFl of the point process $(\Phi,\Xi)$ is given by

[TABLE]

see, e.g., Theorem 1.1 of Moyal (1964) where, taking $p_{d}(x):=p_{d}$ for simplicity, $x\in\Lambda$, we have

[TABLE]

The marginal PGFl of the point process $\Xi$ is given by

[TABLE]

Marginal moments of $\Xi$

The first derivative of $\mathcal{G}_{\Xi}(g)$ is given by

[TABLE]

from which the first-order moment density of $\Xi$ can be computed after setting $g=1$ as

[TABLE]

Similarly, the second derivative of $\mathcal{G}_{\Xi}(g)$ is given by

[TABLE]

from which the second-order factorial moment density of $\Xi$ can be computed after setting $g=1$ as

[TABLE]

$x,y\in\Lambda$, $x\not=y$.

Posterior distribution

In Lemma 3.1 below, we derive the general expression of the Janossy densities of the posterior point process $\Phi|\Xi$ given the sensor measurements $\Xi$. In the sequel, we let $|S|$ denote the cardinality of subsets $S\subset\{1,\ldots,m\}$, and we use the notation $z_{1:m}=(z_{1},\ldots,z_{n})$, while $z_{1:m}\setminus z$ denotes the sequence $(z_{1},\ldots,z_{n})$ with the omission of $z$ if $z\in\{z_{1},\ldots,z_{n}\}$.

Lemma 3.1

The $n$-th conditional Janossy density of $\Phi$ given that $\Xi=z_{1:m}=(z_{1},\ldots,z_{m})$ satisfies

[TABLE]

$m,n\geq 0$*, where

$(i)$ the $(n,m)$-th joint Janossy density of $(\Phi,\Xi)$ is given by*

[TABLE]

*the above sum being over injective mappings $\pi:S\rightarrow\{1,\ldots,n\}$, and

$(ii)$ the Janossy densities of the measurement point process $\Xi$ are given by*

[TABLE]

$m\geq 0$.

Proof. In order to derive the $(n,m)$-th joint Janossy density of $(\Phi,\Xi)$ as in (2.4), we need to compute

[TABLE]

in the directions of the functions $\eta_{1},\ldots,\eta_{n}$, $f_{1},\ldots,f_{m}\in L^{\infty}(\Lambda)$. For a given set $S\subseteq\{1,\ldots,m\}$ we let $\pi:=\left|S\right|$ and denote the elements of $S$ as $S(1),\ldots,S(\pi)$ in increasing order, where $S$ is identified to the mapping $S:\{1,\ldots,\pi\}\rightarrow\{1,\ldots,m\}$. By the Faà di Bruno’s formula, see, e.g., Clark and Houssineau (2012), and (3.3), we have

[TABLE]

where the index set $S\subseteq\{1,\ldots,m\}$ runs through the collection of $2^{m}$ subsets of $\{1,\ldots,m\}$. Next, the $m^{\text{th}}$ derivative of $\mathcal{G}_{\Phi}(h(\cdot)\mathcal{G}_{\Xi_{s}}(g\ \!|\ \!\cdot))$ can be computed by a standard induction argument as

[TABLE]

Substituting $f_{1},\ldots,f_{m}$ with Dirac delta functions $\delta_{z_{1}},\ldots,\delta_{z_{m}}$ at the distinct configuration points $z_{1:m}\in\Lambda$ as in (2.3) and setting $g=0$ in (3.11), we find

[TABLE]

Next, substituting (3.12) and the relation

[TABLE]

into (3.10), we obtain

[TABLE]

Hence we find

[TABLE]

after setting $h=0$. By substituting $\eta_{1},\ldots,\eta_{n}$ with Dirac delta functions $\delta_{x_{1}},\ldots,\delta_{x_{n}}$ at distinct configuration points $x_{1:n}\in\Lambda$, the $(n,m)$-th joint Janossy density of $(\Phi,\Xi)$ is then given by

[TABLE]

which shows (3.8). By (2.4), (3.3), (3.4) and (3.13), we have

[TABLE]

which shows (3.9). Finally, (3.7) follows from the Bayes formula. $\square$

The combinatorics of Lemma 3.1 is similar to Theorem 1 of Lund and Rudemo (2000), which instead computes the conditional likelihood $j^{(m)}_{\Xi\mid\Phi=x_{1:n}}(z_{1},\ldots,z_{m})$ of the observed point process $\Xi$ given a Poisson point process $\Phi$.

Note that (3.8) and (3.9) admit natural combinatorial interpretations by identifying $S^{c}=\{1,\ldots,m\}\setminus S$ to the set of points created according to the Poisson point process with intensity function $l_{c}(z)$, and by letting $n-|S|$ denote the number of points in $\Phi$ deleted with probability $q_{d}$ by the Bernoulli point process $\Xi_{s}$.

Poisson case

In case $\Phi$ is the Poisson point process with intensity measure $\nu(dx)$ we have $j^{(n)}_{\Phi}=e^{-\nu(\Lambda)}$, $n\geq 0$, hence (3.9) recovers the classical expression

[TABLE]

First-order posterior moment

In the next proposition we express the first-order conditional moment of $\Phi$ given the sensor measurements $\Xi=z_{1:m}=(z_{1},\ldots,z_{m})$, using extensions of the corrector terms $l^{(1)}_{z_{1:m}}$ introduced in Delande et al. (2014) for the cardinalized PHD filter, see Equation (19) in Lemma 1 therein, and also Equation (41) in Theorem IV.7 of Schlangen et al. (2018) for the Panjer-based PHD filter.

Proposition 3.2

The first-order conditional moment of $\Phi$ given that $\Xi=(z_{1},\ldots,z_{m})$ is given by its density

[TABLE]

with respect to $\nu(dx)$, $m\geq 0$, where

[TABLE]

are corrector terms, $j^{(m)}_{\Xi}(z_{1:m})=j^{(m)}_{\Xi}(z_{1},\ldots,z_{m})$ is given by (3.9), and

[TABLE]

$m\geq 0$.

Proof. The first-order joint moment density of $\Phi$ with $\Xi=(z_{1},\ldots,z_{m})$ can be obtained from the PGFl (3.3) as

[TABLE]

or, using the joint Janossy densities (3.8) and denoting by $d\hat{x}_{r}$ the absence of $dx_{r}$, as

[TABLE]

and it remains to divide by $j^{(m)}_{\Xi}(z_{1:m})$. $\square$

Second-order posterior moment

Similarly, the second partial moment of the first-order integral of $\Phi$ when $\Xi=z_{1:m}=(z_{1},\ldots,z_{m})$ is obtained in the next proposition, which uses an extension of the corrector terms $l^{(2)}_{z_{1:m}}$ introduced in Delande et al. (2014) for the cardinalized PHD filter, see Equation (29) in Lemma 2 therein, and also Equation (42) in Theorem IV.8 of Schlangen et al. (2018) for the Panjer-based PHD filter.

Proposition 3.3

The second-order conditional factorial moment of $\Phi$ given that $\Xi=(z_{1},\ldots,z_{m})$ is given by its density

[TABLE]

with respect to $\nu(dx)\nu(dy)$, with the corrector terms

[TABLE]

and

[TABLE]

where $j^{(m)}_{\Xi}(z_{1},\ldots,z_{m})$ is as in (3.9), and

[TABLE]

$x,y\in\Lambda$, $m\geq 0$.

Proof. Factorial moments can be computed using the second derivative of the conditional PGFl (3.3), see (2.5), or equivalently using the joint Janossy densities (3.8) as in the proof of Proposition 3.2. We have

[TABLE]

$x,y\in\Lambda$, $x\not=y$, and it remains to divide by $j^{(m)}_{\Xi}(z_{1:m})$. $\square$

Poisson case

In the case of a Poisson point process with $j^{(n)}_{\Phi}=e^{-\nu(\Lambda)}$, $n\geq 0$, (3.16) reads

[TABLE]

and similarly from (3.21) we find

[TABLE]

Hence, in the Poisson case the corrector terms are given by $l^{(1)}_{z_{1:m}}(x)=l^{(2)}_{z_{1:m}}(x,y)=1$ and

[TABLE]

with

[TABLE]

and the first and second (factorial) moment densities (3.14), (3.3) of $\Phi$ with respect to $\nu(dx)$ given the point process $\Xi$ recover the classical expressions

[TABLE]

of first order moment density, see Relation (2.87) in Clark et al. (2016), and the second-order moment density

[TABLE]

$x,y\in\Lambda$, $x\not=y$, $m\geq 0$. See, e.g., Proposition V.1(a) of Schlangen et al. (2018) and Exercise 4.3.4 in Clark et al. (2016).

Posterior covariance

For $A,B$ measurable subsets of ${\mathord{\mathbb{R}}}^{d}$, let

[TABLE]

denote the posterior covariance, where $\mu^{(1)}_{\Phi\mid\Xi=z_{1:m}}(A)$ is the posterior first order moment

[TABLE]

and $\mu^{(2)}_{\Phi\mid\Xi=z_{1:m}}(A,B)$ is the posterior second-order moment

[TABLE]

Using Relations (3.14) and (3.3) in Propositions 3.2-3.3, we obtain the following representation of the posterior covariance.

Proposition 3.4

The posterior covariance $c^{(2)}_{\Phi\mid\Xi=z_{1:m}}(A,B)$ of $\Phi$ given that $\Xi=(z_{1},\ldots,z_{m})$ is given by

[TABLE]

$m\geq 0$.

When $A=B$, Relation (3.22) becomes the variance identity

[TABLE]

which takes a form similar to the variance update formula obtained for the Panjer-based PHD filter, see Equations (41)-(42) of Theorem IV.8 in Schlangen et al. (2018).

Poisson case

In the case of a Poisson point process with $j^{(n)}_{\Phi}=e^{-\nu(\Lambda)}$, $n\geq 0$, Proposition 3.4 recovers the covariance

[TABLE]

see, e.g., Equation (41) and Proposition V.1(a) of Schlangen et al. (2018), and the variance

[TABLE]

see also Exercise 4.3.4 in Clark et al. (2016).

4 Determinantal point processes

In this section we review the properties of determinantal point processes; see, e.g., Decreusefond et al. (2016) and references therein for additional background.

Kernels and integral operators

For any compact set $\Lambda\subseteq{\mathord{\mathbb{R}}}^{d}$, we denote by $L^{2}(\Lambda,{\nu})$ the Hilbert space of square-integrable functions w.r.t. the restriction of the Radon measure ${\nu}$ on $\Lambda$, equipped with the inner product

[TABLE]

By definition, an integral operator $\mathcal{K}:L^{2}({\Lambda},{\nu})\to L^{2}({\Lambda},{\nu})$ with kernel $K:{\Lambda}^{2}\to{\mathord{\mathbb{R}}}$ is a bounded operator defined by

[TABLE]

It can be shown that $\mathcal{K}$ is a compact operator, which is self-adjoint if its kernel verifies

[TABLE]

Equivalently, this means that the integral operator $\mathcal{K}$ is self-adjoint for any compact set $\Lambda\subseteq{\mathord{\mathbb{R}}}^{d}$. If $\mathcal{K}$ is self-adjoint, by the spectral theorem we have that $L^{2}(\Lambda,{\nu})$ has an orthonormal basis $(\varphi_{n})_{n\geq 1}$ of eigenfunctions of $\mathcal{K}$ with corresponding eigenvalues $({\nu}_{n})_{n\geq 1}$, and the kernel $K$ of $\mathcal{K}$ can be written as

[TABLE]

For $\mathcal{K}$ a self-adjoint integral operator of trace class, i.e.

[TABLE]

we define the trace of $\mathcal{K}$ as $\mathrm{Tr\,}\mathcal{K}=\sum_{n\geq 1}{\nu}_{n}$. Let also $\mathrm{Id}$ denote the identity operator on $L^{2}({\Lambda},{\nu})$ and let $\mathcal{K}$ be a trace class operator on $L^{2}({\Lambda},{\nu})$. We define the Fredholm determinant of $\mathrm{Id}+\mathcal{K}$ as

[TABLE]

with the relation

[TABLE]

where \mathrm{det}\big{(}K(x_{i},x_{j})_{1\leq i,j\leq n}\big{)} is the determinant of the $n\times n$ matrix $(K(x_{i},x_{j}))_{1\leq i,j\leq n}$, see Theorem $2.4$ of Shirai and Takahashi (2003), and also Brezis (1983) for more details on Fredholm determinants.

Determinantal point processes

In the sequel we consider a self-adjoint trace class operator $\mathcal{K}_{\Psi}$ on $L^{2}(\Lambda,\nu)$ with spectrum contained in $[0,1)$, and denote by $K_{\Psi}:\Lambda\times\Lambda\longrightarrow{\mathord{\mathbb{R}}}$ the kernel of ${\mathcal{K}}_{\Psi}$.

By the results in Macchi (1975) and Soshnikov (2000) (see also Lemma 4.2.6 and Theorem 4.5.5 in Hough et al. (2009)) the determinantal point process $\Psi$ on ${\Lambda}$, with integral operator ${\mathcal{K}}_{\Psi}$ is defined as in (2.2) by its correlation functions

[TABLE]

w.r.t. the measure ${\nu}$ on $({\Lambda},\mathcal{B}({\Lambda}))$, $x_{1},\ldots,x_{n}\in\Lambda$, with $x_{i}\not=x_{j}$, $1\leq i<j\leq n$, see also Lemma 3.3 of Shirai and Takahashi (2003). In particular, we have

[TABLE]

and

[TABLE]

$x,y\in\Lambda$, $x\not=y$, i.e.

[TABLE]

with $\rho^{(2)}_{\Psi}(x,x):=0$, $x\in\Lambda$. The covariance of the determinantal point process $\Psi$ is then given by

[TABLE]

which shows that the determinantal point process $\Psi$ is negatively correlated, since when $A\cap B=\emptyset$ we have

[TABLE]

The interaction operator ${\mathcal{J}}_{\Psi}$ on $L^{2}({\Lambda},{\nu})$ is defined as

[TABLE]

and has the kernel

[TABLE]

by (4.1). For $\alpha=\{x_{1},\ldots,x_{n}\}\in\mathbf{N}_{\sigma}(\Lambda)$, we denote by $\mathrm{det}\,J_{\Psi}(\alpha)$ the determinant

[TABLE]

which is ${\nu}^{\otimes n}(x_{1},\ldots,x_{n})$-a.e. nonnegative; see, e.g., the appendix of Georgii and Yoo (2005).

By Lemma 3.3 in Shirai and Takahashi (2003) the determinantal point process $\Psi$ on $\Lambda$ with kernel $K_{\Psi}(x,y)$, $x,y\in\Lambda$, admits the Janossy densities

[TABLE]

In addition, from e.g. Shirai and Takahashi (2003) (see Theorem 3.6 therein) the Laplace transform (2.1) of $\Psi$ is given by

[TABLE]

for each nonnegative $f$ on ${\Lambda}$ with compact support, where $\varphi=1-e^{-f}$ and $\mathcal{K}_{\varphi}$ is the trace class integral operator with kernel

[TABLE]

5 Determinantal PHD filter

In this section we construct a second-order PHD filter based on determinantal point processes. We show that approximate closed-form filter update expressions can be derived using approximation formulas stated in appendix for the corrector terms $l^{(1)}_{z_{1:m}}$, $l^{(2)}_{z_{1:m}}$ and Janossy densities $j^{(n)}_{\Phi}$, when the underlying point process has low cross-correlations.

In the sequel we will restrict the class of determinantal kernels considered to a class of finite range interaction point processes, by enforcing the condition

[TABLE]

as in e.g. Proposition 3.9 in Georgii and Yoo (2005), where $d(\Lambda)$ is the diameter of $\Lambda$ and $\eta\in(0,1)$.

Prediction step

The prediction point process $\Phi$ is constructed by branching the prior point process $\Psi$ with a Bernoulli point process $\Phi_{s}$ with probability of survival $p_{s}(x)$ at the point $x\in\Lambda$, spatial likelihood density $l_{s}(\cdot|x)$ from state $x$, and characterized by the PGFl

[TABLE]

According to (3.3), the PGFl of the prediction point process $\Phi$ is given by

[TABLE]

where $\mathcal{G}_{\Phi_{b}}$ is the PGFl of the Poisson birth point process $\Phi_{b}$ of new targets. In the sequel we use the notation convention (3.2), i.e. $\tilde{l}_{s}(x|u):=p_{s}(u)l_{s}(x|u)$, for compactness of notation.

Proposition 5.1

Assume that the prior point process $\Psi$ is a determinantal point process with kernel $K_{\Psi}(x,y)$. Then, the prediction first and second-order moment densities of $\Phi$ are given by

[TABLE]

and

[TABLE]

$x,y\in\Lambda$, $x\not=y$.

Proof. The expressions (5.3)-(5.4) of the prediction first and second-order (factorial) moment densities are obtained from (3.5) and (3) as

[TABLE]

and

[TABLE]

$\square$

From Proposition 5.1 we can model the prediction point process $\Phi$ as a determinantal process with prediction kernel $K_{\Phi}$, whose diagonal entries are given by

[TABLE]

and whose nondiagonal entries satisfy

[TABLE]

from (4.3), where $\rho^{(2)}_{\Phi}(x,y)$ is given by (5.4) when $x\not=y$, and $\rho^{(2)}_{\Phi}(x,x):=0$, $x\in\Lambda$. The prediction Janossy kernel $J_{\Phi}(x,y)$ of the operator ${\mathcal{J}}_{\Phi}$ is then computed by the formula

[TABLE]

see (4.7).

Update step

First order moment update. Proposition 3.2 and (4.2) show that the diagonal values of the posterior kernel $K_{\Phi\mid\Xi=z_{1:m}}$ are given by

[TABLE]

Based on the approximation of the corrector term $l^{(1)}_{z_{1:m}}$ stated in Proposition A.1, we will estimate the first-order posterior moment density as

[TABLE]

where we choose to approximate $l^{(1)}_{z_{1:m}}(x)\simeq\mu^{(1)}_{\Phi}(x)=K_{\Phi}(x,x)$ for consistency with the standard Poisson PHD filter. The estimate (5.5) allows one to locate the targets by maximizing $K_{\Phi\mid\Xi=z_{1:m}}(x,x)=\mu^{(1)}_{\Phi\mid\Xi=z_{1:m}}(x)$ over $x$, and to estimate the number of targets as

[TABLE] 2. 2.

Cross-diagonal kernel update. As a consequence of (4.4), i.e.

[TABLE]

$x,y\in\Lambda$, the cross-diagonal entries of the posterior kernel $K_{\Phi\mid\Xi=z_{1:m}}(x,y)$ can be estimated from (4.3) as

[TABLE]

$x,y\in\Lambda$. The above relation (5.7) can be rewritten from Propositions 3.2-3.3 as

[TABLE]

$x,y\in\Lambda$, $m\geq 0$. In practice we will estimate the posterior kernel $K_{\Phi\mid\Xi=z_{1:m}}(x,y)$ in (5.7) using (5.5) and the approximation of the corrector term $l^{(2)}_{z_{1:m}}$ in Proposition A.2, to obtain

[TABLE]

with $\rho^{(2)}_{\Phi\mid\Xi=z_{1:m}}(x,x):=0$, $x\in\Lambda$, see (A) which yields the expression of $K_{\Phi\mid\Xi=z_{1:m}}(x,y)^{2}$ in Proposition A.3, where $s_{c}(z)$ is defined in (A.5). After completing the update step, we move to the next prediction step by taking $K_{\Psi}(x,y):=K_{\Phi\mid\Xi=z_{1:m}}(x,y)$, $x,y\in\Lambda$.

6 Implementation

We implement the Determinantal Point Process (DPP) and Poisson Point Process (PPP) PHD filters using the sequential Monte Carlo (or particle filtering) method as in Li et al. (2017), together with the roughening method of Li et al. (2013), which allow us to estimate otherwise intractable integrals using discretized particle summations. Our ground truth dynamics follows the nearly constant turn-rate motion dynamics of Vo et al. (2009), Li et al. (2017), with the addition of a repulsion term. The state of each target at time $t$ is given by $\mathbf{x}_{t}=\left(x_{t},\dot{x}_{t},y_{t},\dot{y}_{t},\theta_{t}\right)^{\top}$, where $x_{t}$, $y_{t}$ are the cartesian coordinates, $\dot{x}_{t}$, $\dot{y}_{t}$ are the respective velocities, and $\theta_{t}$ is the turn rate. At time $t+1$, the location of every target $i$ for $i\in\{1,\ldots,n\}$ is given by

[TABLE]

where

[TABLE]

is a term which models repulsion among targets.

Here, $\mathbf{v}_{t}=(v_{x},v_{y},v_{\theta})^{\top}$ is a zero-mean acceleration noise distributed according to the zero-mean Gaussian noise

[TABLE]

and

[TABLE]

[TABLE]

with $\tau>0$ the time sampling period. When $\mathbf{F}=0$ and $\mathbf{G}=I_{d}$, the repulsive interaction motion dynamics (6.1) has the law the Ginibre DPP for stationary distribution; see, for example, Equation (2.19) in § 2.2 of Osada (2013).

The measurement vector of each target at time $t$ is written as $\mathbf{m}_{t}=(m_{r_{t}},m_{\omega_{t}})^{\top}$ using the range and bearing components $m_{r_{t}}$ and $m_{\omega_{t}}$. The measurement generated by every target at time $t+1$ is then given by

[TABLE]

where \mathbf{p}_{t+1}=\left(\begin{array}[]{c}\sqrt{x^{2}_{t+1}+y^{2}_{t+1}}\\ \arctan(y_{t+1}/x_{t+1})\end{array}\right) and the measurement noise vector $\mathbf{w}_{t+1}=(w_{r_{t+1}},w_{\omega_{t+1}})^{\top}$ is distributed according to the zero-mean Gaussian noise

[TABLE]

The spatial likelihood densities $l_{s}(\mathbf{z}|\mathbf{x})$ and $\tilde{l}_{d}(\mathbf{z}|\mathbf{x})$ from a target state $\mathbf{x}$ to a measurement $\mathbf{z}$ in (3.1) and (5.2) respectively follow the zero-mean multivariate Gaussian distribution (6.8) of $\mathbf{v}_{t}$ and the multivariate Gaussian distribution (6.14) of $\mathbf{w}_{t}$. In addition, the model generates measurement information from every target with a constant probability of detection $p_{d}$, and the spatially distributed clutter measurement points are generated according to a Poisson point process with constant density $l_{c}(z)$.

The implementation of Figures 1 to 3 use $1000$ particles at initialization, $100$ resampling particles per expected target, and $100$ new particles per expected target birth, which follows a time-dependent Poisson birth process. The starting locations are uniformly distributed within a square domain, we take the spatial standard deviations (s.d.) $\sigma_{v_{x}}=\sigma_{v_{y}}=1\ \mathrm{m}/\mathrm{s}^{2}$, the turn-rate noise s.d. $\sigma_{v_{\theta}}=\pi=180\ \mathrm{rad}/\mathrm{s}$, bearing distribution s.d. $\sigma_{w_{\omega}}=\pi=180\ \mathrm{rad}$, and range distribution s.d. $\sigma_{w_{r}}=2\sqrt{2}\ \mathrm{m}$.

For illustration and performance assessment purposes, our code displays the association between target-originated measurements and posterior state estimates. For this, given a measurement we select its associated estimate by minimizing the distance between the measurement and all candidate estimates. In addition, the blue edges show the estimates which improve over the corresponding measurements in terms of Euclidean distances to the ground truth, while the orange edges show the estimates which perform worse than measurements.

Our simulations also display a ratio of good estimate counts against total measurement counts, as well as a gain metric which measures the relative improvement in distance between estimates and measurements. Positive gain correspond to a good estimate ratio above $50\%$, and negative gain is realized when the ratio falls below $50\%$.

Figure 1-(1(b)) shows that when using a nonzero value for the repulsion parameter $\zeta=\zeta_{x}=\zeta_{y}$, the good estimate ratio with repulsive interaction becomes lower as compared with the non repulsive setting of Figure 1-(1(a)), with $50$ time steps. In Figure 1 the Poisson clutter rate is $l_{c}=1$, the probability of detection $p_{d}=0.9$ and the probability of survival is $p_{s}=1$, with four targets.

In Figures 2-(2(a))-(2(b)) we provide further illustrations of three-target interaction and PHD filter output of a single trial at different repulsion values $\zeta=0$ and $30$, with $15$ time steps and $p_{d}=p_{s}=1$.

Figure 2-(2(a)) shows the PHD filter output with $\zeta=0$, where the targets can become closer to each other without repulsion, and with positive gain. For the repulsion value $\zeta=30$ as in Figure 2-(2(b)) the repulsion effect among the three targets become much more evident, the good estimate ratio becomes lower, and the gain becomes negative.

The above results are summarized in Figure 3. Figure 3-(3(a)) presents the SMC-PHD filter first-order moment output with $200$ Monte Carlo runs and $10$ targets across $20$ time steps, with different repulsion parameter values $\zeta=0,4,8$. Figure 3-(3(b)) presents the good estimate ratio for various values of the repulsion coefficient $\zeta$, with $100$ Monte Carlo runs on three targets across $15$ time steps, with $p_{d}=1$ and no clutter.

We note that the PHD filter is correctly estimating the target count when the repulsion coefficient $\zeta$ vanishes; however, the estimation falls short for nonzero values of $\zeta$, showing the performance degradation of the PPP-PHD filter in the presence of target interaction.

7 Determinantal PHD filtering algorithm

Initialization $(t=0)$

The state dynamics of the initial set of $N_{\Phi,0}$ particles is sampled according to a uniform distribution on the state space $\Lambda$. The diagonal entries of the prior discretized determinantal kernel $K_{\Phi,0}$ at time $t=0$ are initialized to $\gamma_{\Phi,0}/N_{\Phi,0}$ with $\gamma_{\Phi,0}$ a prior intensity value. The nondiagonal entries are initialized to $\alpha\gamma_{\Phi,0}/N_{\Phi,0}$ where $\alpha\geq 0$, except for those which are set to zero according to Condition (5.1) with the matrix index threshold $\eta=10\%$.

Using (4.7) we then compute the discretized Janossy kernel $J_{\Phi,0}$ which is needed for the evaluation of the posterior determinantal kernel $\tilde{K}_{\Phi|\Xi,0}$. Letting $P_{p}$ denote the number of resampled particles per target and $\gamma_{0}:=\sum_{i=1}^{N_{\Phi,0}}K_{\Phi|\Xi,0}(\mathbf{x}_{i},\mathbf{x}_{i})$, we resample $\tilde{N}_{0}:=P_{p}\times\lfloor\gamma_{0}\rfloor$ particles $\{\mathbf{x}_{i}\}_{i=1}^{\tilde{N}_{0}}$ that better describe the target locations as in Li et al. (2013) by maximizing the diagonal entries of $\tilde{K}_{\Phi|\Xi,0}$ over $\Lambda$, where $\lfloor\cdot\rfloor$ denotes the integer floor function. Those particles are then used to initialize the post-resampling determinantal kernel $K_{0}$ and to compute the post-resampling Janossy kernel $J_{0}$ in order to estimate the updated kernel $\tilde{K}_{\Phi|\Xi}$ from (5.5), (5.7) and (2). The updated kernel $\tilde{K}_{\Phi|\Xi}$ is then set as the prior kernel $K_{\Psi}$ of the next time step.

Algorithm $(t\geq 1)$

The general algorithm proceeds to compute the prediction state transition dynamics $\{\mathbf{x}_{t+1|t}^{(i)}\}_{i=1}^{\tilde{N}_{t}}$ using (6.1), followed by the computation of the prediction determinantal state transition kernel $K_{\Phi,t+1|t}$ using (5.3) and (5.4). Letting $P_{b}$ denote the number of particles per birth target and $\gamma_{\Phi,t+1}:=\sum_{i=1}^{\tilde{N}_{t}}K_{\Psi,t+1}(\mathbf{x}_{i},\mathbf{x}_{i})\nu(\{x_{i}\})$, we sample the state dynamics of $N_{\Phi,t+1}:=P_{b}\times\lfloor\gamma_{\Phi,t+1}\rfloor$ particles for the target birth process $\Phi$. The discretized prediction determinantal kernel $K_{\Phi,t+1|t}$ is then extended to incorporate the set of additional $N_{\Phi,t+1}$ particles by assigning the diagonal entries corresponding to these new particles to $\gamma_{\Phi,t+1}/N_{\Phi,t+1}$ and the nondiagonal entries to $\alpha\gamma_{\Phi,t+1}/N_{\Phi,t+1}$, and by setting all other new entries to [math] according to Condition (5.1) with the matrix index threshold $\eta:=10\%$. Thereafter, we compute the discretized Janossy kernel $J_{\Phi,t+1|t}$ using (4.7) and then the discretized posterior determinantal kernel $K_{\Phi|\Xi,t+1}$ using (5.5), (5.7) and (2). Next, letting $\gamma_{t+1}:=\sum_{i=1}^{N_{t+1}}K_{\Phi|\Xi,t+1}(\mathbf{x}_{i},\mathbf{x}_{i})$ we resample $\tilde{N}_{t+1}:=P_{p}\times\lfloor\gamma_{t+1}\rfloor$ particles with state dynamics $\{\mathbf{x}_{i}\}_{i=1}^{\tilde{N}_{t+1}}$ as in Li et al. (2013), by maximizing the diagonal entries of $K_{\Phi|\Xi,t+1}$ over $\Lambda$. Those particles are then used to initialize the post-resampling determinantal kernel $K_{t+1}$ by setting diagonal entries to $\gamma_{t+1}/\tilde{N}_{t+1}$ and nondiagonal entries to $\alpha\gamma_{t+1}/\tilde{N}_{t+1}$, except for those which are set to zero according to Condition (5.1) with $\eta=10\%$. Finally, we recompute the post-resampling Janossy kernel $J_{t+1}$ and the posterior determinantal kernel $\tilde{K}_{\Phi|\Xi,t+1}$ using (5.5), (5.7) and (2).

The complexity of this DPP-PHD filter is cubic in the number of discretization steps due to the presence of matrix inversions in the algorithm.

Numerical results

In Figure 4 we assess the spooky effect (see Fränken et al. (2009)) of our DPP-PHD filter, following the approach applied in Schlangen et al. (2018) to second-order PHD filters. Our tracking scenario consists of two disjoint square domains $A$ and $B$ of size $150$ m by $150$ m, which are located 150m diagonally apart. In each domain, $10$ targets are initialized and their state dynamics are centrally distributed at the first time step. The runtime of the experiment is set at $50$ with $100$ Monte Carlo (MC) runs, the targets survive throughout ($p_{s}=1$), and their trajectories remain within the observation domains. In Figure 4 we take the spatial standard deviations (s.d.) $\sigma_{v_{x}}=\sigma_{v_{y}}=1.0$ m/s2, turn-rate noise s.d. $\sigma_{v_{\theta}}=\pi=180$ rad/s, bearing distribution s.d. $\sigma_{w_{\omega}}=\pi=180$ rad, and range distribution s.d. $\sigma_{w_{r}}=\sqrt{2}$ m.

In a similar setting to Schlangen et al. (2018), all targets in domain $B$ are compelled to be misdetected in every cycle of $10$ time steps. We use a constant probability of detection $p_{d}=0.9$ and mean clutter count at $5$ in each measurement space.

At initialization in Figure 4 we set $N_{\Phi,0}=800$, $\gamma_{\Phi,0}=2$ and $\alpha=4$. The DPP-PHD filter implementation uses $P_{p}=30$ resampled particles per target, and $P_{b}=10$ particles per birth target.

Figure 4-(4(a)) shows the estimated intensities in domains $A$ and $B$, where domain $A$ is unaffected by the rapid drop in the intensity of domain $B$. The posterior correlation estimates in Figure 4-(4(b)) are computed by rescaling the covariance expression (4.6) written as

[TABLE]

as in Corollary A.4, where $K_{\Phi\mid\Xi=z_{1:m}}(x,y)$ is estimated as in Proposition A.3 from (5.5) and (5.7). Figure 4-(4(b)) shows negative correlations due to the determinantal point process nature, which leads to a drop in negative correlation during the compelled misdetection at each $10$-steps cycle.

Figure 5 presents miss-distance performance estimates for the experiment of Figure 4, using the $L^{2}$-Optimal Mass Transfer (OMAT, Hoffman and Mahler (2004)) metric, and the $L^{2}$-Optimal Subpattern Assignment (OSPA, Schumacher et al. (2008)) metric with threshould $c=100$, which solves the inconsistencies encountered with the OMAT metric and takes into account differences in cardinalities.

In Figure 6 we compare the robustness of the DPP and PPP-PHD filters when both filters are subjected to sudden death in the number of targets in a single domain of size $100$m by $100$m, beginning with $15$ targets at the first time step.

Figure 6-(6(a)) uses $300$ Monte Carlo runs, while Figure 6-(6(b)) relies on $200$ Monte Carlo runs. The runtime of each Monte Carlo run spans from time $t=0$ to time $t=15$, and the probability of survival is $p_{s}:=1$. The initial $15$ targets are maintained until time $t=9$ when $10$ random targets are forced to die and the remaining $5$ targets survive until the end of the time interval. In Figure 6 we take the spatial standard deviations (s.d.) $\sigma_{v_{x}}=\sigma_{v_{y}}=1.0$ m/s2, turn-rate noise s.d. $\sigma_{v_{\theta}}=\pi=180$ rad/s, with bearing and range distribution s.d. $\sigma_{w_{\omega}}=\pi=180$ rad, $\sigma_{w_{r}}=\sqrt{2}$ m as in Figure 4, with probability of detection $p_{d}=0.95$, mean clutter count at $1$ up to time $t=9$ and then at $0.06$ afterwards for Figure 6-(6(b)), and mean clutter count at $1$ up to time $t=9$ and then at $0.3$ afterwards for Figure 6-(6(a)). At initialization in Figure 6, we set $N_{\Phi,0}=6000$ and $\gamma_{\Phi,0}=0.2$. Both our DPP and PPP-PHD filter implementations use $P_{p}=50$ resampled particles per target in Figure 6, $P_{b}=60$ particles per birth target in Figure 6-(6(b)), and $P_{b}=40$ particles per birth target in Figure 6-(6(a)).

In Figure 7 we compare the robustness and performance recovery of the DPP and PPP-PHD filters when subjected to a rapid birth in the number of targets in a single domain of size $100$ m by $100$ m. The experiment starts with a single target which survives throughout the $45$ time steps, without birth of new targets from time $t=0$ to time $t=9$. At time $t=10$, $9$ new targets are born centrally distributed within the target space and survive through the remaining time steps.

Each Monte Carlo run spans $45$ time steps, with $400$ and $100$ Monte Carlo runs in the experiments of Figures 7-(7(a)) and 7-(7(b)) respectively. In Figure 7 the spatial standard deviations (s.d.) $\sigma_{v_{x}}=\sigma_{v_{y}}=1.0$ m/s2, turn-rate noise s.d. $\sigma_{v_{\theta}}=\pi=180$ rad/s, and bearing and range s.d. $\sigma_{w_{\omega}}=\pi=180$ rad, $\sigma_{w_{r}}=\sqrt{2}$ m are the same as in Figure 6. The model generates measurement information from each target with a constant probability of detection $p_{d}=0.90$, mean clutter count at [math] up to time $t=9$ and then at $5$ afterwards for Figure 7-(7(a)), and mean clutter count at $0.05$ up to time $t=9$ and then at $5$ afterwards for Figure 7-(7(b)). We set $N_{\Phi,0}=300$ and $\gamma_{\Phi,0}=0.2$ at initialization in Figure 7. Both DPP and PPP-PHD filter implementations use $P_{b}=40$ and $P_{b}=50$ resampled particles per target in Figure 7-(7(a)) and Figure 7-(7(b)) respectively. For the target birth process we set $P_{b}=9$ and $P_{b}=15$ particles per birth target in Figure 7-(7(a)) and Figure 7-(7(b)) respectively.

Appendix A Appendix - Janossy density approximation

Since the corrector terms $l^{(1)}_{z_{1:m}}(x)$, $l^{(1)}_{z_{1:m}}(x;z)$, $l^{(2)}_{z_{1:m}}(x,y)$, $l^{(2)}_{z_{1:m}}(x,y;z)$, $l^{(2)}_{z_{1:m}}(x,y;z,z^{\prime})$ in (3.15), (3.19) and the kernel update formula (2) have no closed form expression in the determinantal setting, we propose to use the Janossy density approximations

[TABLE]

$n\geq 1$, which corresponds to a (Poisson) first-order approximation, and

[TABLE]

$n\geq 2$, which corresponds to a second-order (determinant) approximation, obtained from (4.8) by assuming that the off-diagonal entries $J_{\Phi}(x_{i},x_{j})$, $i\not=j$, are small.

This Janossy approximation is specially relevant to $\alpha$-determinantal Ginibre point processes (GPP) which approximate a Poisson point process when $\alpha\in[-1,0)$ tends to [math], see Shirai and Takahashi (2003).

Proposition A.1

Under (A.1) we have the first-order Poisson approximations $l^{(1)}_{z_{1:m}}(x)\simeq J_{\Phi}(x,x)$, $m\geq 0$, and

[TABLE]

$z\in z_{1:m}$, $x\in\Lambda$, $m\geq 1$.

Proof. By (3.16) and (A.1) we have

[TABLE]

by (3.9), which yields the approximation $l^{(1)}_{z_{1:m}}(x)\simeq J_{\Phi}(x,x)$. On the other hand, for $r=1,\ldots,m$, using again (A.1) and (3.9) we have

[TABLE]

We conclude by taking $z_{r}=z$ and noting that by (3.15) and (A.3)-(A.4) we have

[TABLE]

$\square$

Proposition A.2

Under (A.1)-(A) we have the second-order approximations

[TABLE]

$z\in z_{1:m}$, $x,y\in\Lambda$, $m\geq 1$, and

[TABLE]

$z,z^{\prime}\in z_{1:m}$, $z\not=z^{\prime}$, $x,y\in\Lambda$, $m\geq 2$, where

[TABLE]

Proof. By (3.21) and (A) we have

[TABLE]

and for $r,u=1,\ldots,m$, using (A.1)-(A) and (3.9) we find

[TABLE]

We conclude by taking $(z_{r},z_{u})=(z,z^{\prime})$ and noting that by (3.20) and (A.6)-(A.7) we have

[TABLE]

$z,z^{\prime}\in z_{1:m}$, $z\not=z^{\prime}$, $m\geq 2$. Similarly, by (3.19) and (A.4), (A.6) we also have

[TABLE]

$\square$

As a consequence of (3.3) and Proposition A.2, the second-order conditional factorial moment density of $\Phi$ given that $\Xi=z_{1:m}=(z_{1},\ldots,z_{m})$ will be approximated as

[TABLE]

$m\geq 0$, with $\rho^{(2)}_{\Phi\mid\Xi=z_{1:m}}(x,x):=0$, $x\in\Lambda$.

Proposition A.3

The (approximate) kernel update formula is given by

[TABLE]

$m\geq 0$, $x,y\in\Lambda$.

Proof. By (3.14) and Proposition A.1, we have the approximation

[TABLE]

hence by (A) and (A), we find

[TABLE]

$m\geq 0$, and we conclude by (4.4), i.e.

[TABLE]

and (A.4). $\square$

The next result, which provides an approximation formula for the posterior covariance of Proposition 3.4, is a consequence of Proposition A.3 and (A).

Corollary A.4

Under (A.1)-(A) the posterior covariance of $\Phi$ given that $\Xi=z_{1:m}=(z_{1},\ldots,z_{m})$ is approximated as

[TABLE]

Conclusion

Our observations have shown that the performance of the multi-target tracking PPP-based standard PHD filter is degraded in the presence of target interaction such as repulsion. To address this issue, we have constructed a second-order DPP-based PHD filter based on Determinantal Point Processes which are able to model repulsion between targets, and can propagate variance and covariance information in addition to first-order target count estimates. We have derived posterior moment formulas for the estimation of DPPs after thinning and superposition with a Poisson Point Process (PPP), based on suitable approximation formulas. Our numerical experiments include an assessment of the spooky effect on disjoint domains, with negative correlation estimates which are consistent with the nature of DPPs. We have also compared the robustness and performance recovery of the DPP and PPP-PHD filters when subjected to sudden changes in target numbers.