Passive Multi-Target Tracking Using the Adaptive Birth Intensity PHD   Filter

Christopher Berry; Donald J. Bucci; Samuel W. Schmidt

arXiv:1906.02396·eess.SY·August 11, 2021

Passive Multi-Target Tracking Using the Adaptive Birth Intensity PHD Filter

Christopher Berry, Donald J. Bucci, Samuel W. Schmidt

PDF

TL;DR

This paper introduces an adaptive multi-sensor PHD filter for passive multi-target tracking, improving particle sampling efficiency by leveraging online measurements, demonstrated with TDOA and FDOA data.

Contribution

It extends the adaptive birth intensity PHD filter to multi-sensor scenarios, enabling more efficient target birth particle sampling based on real-time sensor data.

Findings

01

Enhanced particle density efficiency over uniform sampling

02

Effective integration of TDOA and FDOA measurements

03

Improved tracking performance in computationally limited platforms

Abstract

Passive multi-target tracking applications require the integration of multiple spatially distributed sensor measurements to distinguish true tracks from ghost tracks. A popular multi-target tracking approach for these applications is the particle filter implementation of Mahler's probability hypothesis density (PHD) filter, which jointly updates the union of all target state space estimates without requiring computationally complex measurement-to-track data association. Although this technique is attractive for implementation in computationally limited platforms, the performance benefits can be significantly overshadowed by inefficient sampling of the target birth particles over the region of interest. We propose a multi-sensor extension of the adaptive birth intensity PHD filter described in (Ristic, 2012) to achieve efficient birth particle sampling driven by online sensor…

Figures4

Click any figure to enlarge with its caption.

Equations82

X_{k} = {x_{k, 1}, \dots, x_{k, n_{k}}} \in F (X)

X_{k} = {x_{k, 1}, \dots, x_{k, n_{k}}} \in F (X)

Z_{k} = {z_{k, 1}, \dots, z_{k, m_{k}}} \in F (Z)

Z_{k} = {z_{k, 1}, \dots, z_{k, m_{k}}} \in F (Z)

D_{k ∣ k - 1} (x) = γ_{k ∣ k - 1} (x) + p_{s} (x) \int π_{k ∣ k - 1} (x ∣ x^{'}) D_{k - 1∣ k - 1} (x) d x^{'}

D_{k ∣ k - 1} (x) = γ_{k ∣ k - 1} (x) + p_{s} (x) \int π_{k ∣ k - 1} (x ∣ x^{'}) D_{k - 1∣ k - 1} (x) d x^{'}

D_{k ∣ k} (x) = (1 - p_{d} (x)) D_{k ∣ k - 1} (x) + z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) p _{d} ( x )}{κ _{k} ( z _{k} ) + C ( z _{k} )} D_{k ∣ k - 1} (x)

D_{k ∣ k} (x) = (1 - p_{d} (x)) D_{k ∣ k - 1} (x) + z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) p _{d} ( x )}{κ _{k} ( z _{k} ) + C ( z _{k} )} D_{k ∣ k - 1} (x)

C (z_{k}) = \int g_{k} (z_{k} ∣ x^{'}) p_{d} (x^{'}) D_{k ∣ k - 1} (x^{'}) d x^{'} .

C (z_{k}) = \int g_{k} (z_{k} ∣ x^{'}) p_{d} (x^{'}) D_{k ∣ k - 1} (x^{'}) d x^{'} .

D_{k ∣ k - 1} (x, 0) = p_{s} (x) \int π_{k ∣ k - 1} (x ∣ x^{'}) [D_{k - 1∣ k - 1} (x, 1) + D_{k - 1∣ k - 1} (x, 0)] d x^{'}

D_{k ∣ k - 1} (x, 0) = p_{s} (x) \int π_{k ∣ k - 1} (x ∣ x^{'}) [D_{k - 1∣ k - 1} (x, 1) + D_{k - 1∣ k - 1} (x, 0)] d x^{'}

D_{k ∣ k - 1} (x, 1) = γ_{k ∣ k - 1} (x) .

D_{k ∣ k - 1} (x, 1) = γ_{k ∣ k - 1} (x) .

D_{k ∣ k} (x, 0) = (1 - p_{d} (x)) D_{k ∣ k - 1} (x, 0) + z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) p _{d} ( x ) D _{k ∣ k - 1} ( x , 0 )}{κ _{k} ( z _{k} ) + C ^{'} ( z _{k} )}

D_{k ∣ k} (x, 0) = (1 - p_{d} (x)) D_{k ∣ k - 1} (x, 0) + z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) p _{d} ( x ) D _{k ∣ k - 1} ( x , 0 )}{κ _{k} ( z _{k} ) + C ^{'} ( z _{k} )}

D_{k ∣ k} (x, 1) = z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) γ _{k ∣ k - 1} ( x )}{κ _{k} ( z _{k} ) + C ^{'} ( z _{k} )}

D_{k ∣ k} (x, 1) = z_{k} \in Z_{k} \sum \frac{g _{k} ( z _{k} ∣ x ) γ _{k ∣ k - 1} ( x )}{κ _{k} ( z _{k} ) + C ^{'} ( z _{k} )}

C^{'} (z_{k}) = \int g_{k} (z_{k} ∣ x^{'}) γ_{k ∣ k - 1} (x^{'}) d x^{'} + \int g_{k} (z_{k} ∣ x^{'}) p_{d} (x^{'}) D_{k ∣ k - 1} (x^{'}, 0) d x^{'} .

C^{'} (z_{k}) = \int g_{k} (z_{k} ∣ x^{'}) γ_{k ∣ k - 1} (x^{'}) d x^{'} + \int g_{k} (z_{k} ∣ x^{'}) p_{d} (x^{'}) D_{k ∣ k - 1} (x^{'}, 0) d x^{'} .

D_{k ∣ k} (x, 0) \approx i = 1 \sum N_{k}^{p} w_{k, p}^{(i)} δ_{x_{k, p}^{(i)}} (x)

D_{k ∣ k} (x, 0) \approx i = 1 \sum N_{k}^{p} w_{k, p}^{(i)} δ_{x_{k, p}^{(i)}} (x)

D_{k ∣ k} (x, 1) \approx i = 1 \sum N_{k}^{b} w_{k, b}^{(i)} δ_{x_{k, b}^{(i)}} (x)

D_{k ∣ k} (x, 1) \approx i = 1 \sum N_{k}^{b} w_{k, b}^{(i)} δ_{x_{k, b}^{(i)}} (x)

D_{k ∣ k} (x, 0) + D_{k ∣ k} (x, 1) \approx i = 1 \sum N_{k} w_{k}^{(i)} δ_{x_{k}^{(i)}} (x)

D_{k ∣ k} (x, 0) + D_{k ∣ k} (x, 1) \approx i = 1 \sum N_{k} w_{k}^{(i)} δ_{x_{k}^{(i)}} (x)

\cup_{i = 1}^{N_{k}} (x_{k}^{(i)}, w_{k}^{(i)}) = [\cup_{i = 1}^{N_{k}^{p}} (x_{k, p}^{(i)}, w_{k, p}^{(i)})] ⋃ [\cup_{i = 1}^{N_{k}^{b}} (x_{k, b}^{(i)}, w_{k, b}^{(i)})]

\cup_{i = 1}^{N_{k}} (x_{k}^{(i)}, w_{k}^{(i)}) = [\cup_{i = 1}^{N_{k}^{p}} (x_{k, p}^{(i)}, w_{k, p}^{(i)})] ⋃ [\cup_{i = 1}^{N_{k}^{b}} (x_{k, b}^{(i)}, w_{k, b}^{(i)})]

w_{k ∣ k - 1, p}^{(i)} = p_{s} (x_{k, p}^{(i)}) \frac{π _{k ∣ k - 1} ( x _{k, p}^{(i)} ∣ x _{k - 1}^{(i)} )}{q _{k} ( x _{k, p}^{(i)} ∣ x _{k - 1}^{(i)} , Z _{k} )} w_{k - 1}^{(i)} .

w_{k ∣ k - 1, p}^{(i)} = p_{s} (x_{k, p}^{(i)}) \frac{π _{k ∣ k - 1} ( x _{k, p}^{(i)} ∣ x _{k - 1}^{(i)} )}{q _{k} ( x _{k, p}^{(i)} ∣ x _{k - 1}^{(i)} , Z _{k} )} w_{k - 1}^{(i)} .

w_{k ∣ k - 1, b}^{(i)} = \frac{ν _{k ∣ k - 1}^{b}}{N _{k}^{b}}

w_{k ∣ k - 1, b}^{(i)} = \frac{ν _{k ∣ k - 1}^{b}}{N _{k}^{b}}

w_{k ∣ k, p}^{(i)} = (1 - p_{d} (x_{k ∣ k - 1, p}^{(i)})) w_{k ∣ k - 1, p}^{(i)} + z \in Z_{k} \sum \frac{p _{d} ( x _{k ∣ k - 1, p}^{(i)} ) g _{k} ( z ∣ x _{k ∣ k - 1, p}^{(i)} ) w _{k ∣ k - 1, p}^{(i)}}{L ( z )}

w_{k ∣ k, p}^{(i)} = (1 - p_{d} (x_{k ∣ k - 1, p}^{(i)})) w_{k ∣ k - 1, p}^{(i)} + z \in Z_{k} \sum \frac{p _{d} ( x _{k ∣ k - 1, p}^{(i)} ) g _{k} ( z ∣ x _{k ∣ k - 1, p}^{(i)} ) w _{k ∣ k - 1, p}^{(i)}}{L ( z )}

L (z) = κ_{k} (z) + i = 1 \sum N_{k}^{b} w_{k ∣ k - 1, b}^{(i)} + i = 1 \sum N_{k - 1} p_{d} (x_{k ∣ k - 1, p}^{(i)}) g_{k} (z ∣ x_{k ∣ k - 1, p}^{(i)}) w_{k ∣ k - 1, p}^{(i)} .

L (z) = κ_{k} (z) + i = 1 \sum N_{k}^{b} w_{k ∣ k - 1, b}^{(i)} + i = 1 \sum N_{k - 1} p_{d} (x_{k ∣ k - 1, p}^{(i)}) g_{k} (z ∣ x_{k ∣ k - 1, p}^{(i)}) w_{k ∣ k - 1, p}^{(i)} .

w_{k ∣ k, b}^{(i)} = z \in Z_{k} \sum \frac{w _{k ∣ k - 1, b}^{(i)}}{L ( z )}

w_{k ∣ k, b}^{(i)} = z \in Z_{k} \sum \frac{w _{k ∣ k - 1, b}^{(i)}}{L ( z )}

Δt_{k}^{(l_{1} - l_{2})} = {Δ t_{k, 1}^{(l_{1} - l_{2})}, \dots, Δ t_{k, m_{k}^{(l_{1} - l_{2})}}^{(l_{1} - l_{2})}}

Δt_{k}^{(l_{1} - l_{2})} = {Δ t_{k, 1}^{(l_{1} - l_{2})}, \dots, Δ t_{k, m_{k}^{(l_{1} - l_{2})}}^{(l_{1} - l_{2})}}

Δf_{k}^{(l_{1} - l_{2})} = {Δ f_{k, 1}^{(l_{1} - l_{2})}, \dots, Δ f_{k, m_{k}^{(l_{1} - l_{2})}}^{(l_{1} - l_{2})}}

Δf_{k}^{(l_{1} - l_{2})} = {Δ f_{k, 1}^{(l_{1} - l_{2})}, \dots, Δ f_{k, m_{k}^{(l_{1} - l_{2})}}^{(l_{1} - l_{2})}}

Δ t_{k, i}^{(l_{1} - l_{2})} \sim N (\frac{1}{c} (∣ r_{l_{1}} ∣ - ∣ r_{l_{2}} ∣), σ_{Δ t}^{2})

Δ t_{k, i}^{(l_{1} - l_{2})} \sim N (\frac{1}{c} (∣ r_{l_{1}} ∣ - ∣ r_{l_{2}} ∣), σ_{Δ t}^{2})

Δ f_{k, i}^{(l_{1} - l_{2})} \sim N (\frac{f _{c}}{c} (∣ \overset{r}{˙}_{l_{1}} ∣ - ∣ \overset{r}{˙}_{l_{2}} ∣), σ_{Δ f}^{2})

Δ f_{k, i}^{(l_{1} - l_{2})} \sim N (\frac{f _{c}}{c} (∣ \overset{r}{˙}_{l_{1}} ∣ - ∣ \overset{r}{˙}_{l_{2}} ∣), σ_{Δ f}^{2})

x_{k ∣ k - 1, b}^{(n)} \sim b (x_{k ∣ k - 1, b}^{(n)} ∣ (Δt_{k}^{(l_{1} - l_{2})}, Δf_{k}^{(l_{1} - l_{2})})) .

x_{k ∣ k - 1, b}^{(n)} \sim b (x_{k ∣ k - 1, b}^{(n)} ∣ (Δt_{k}^{(l_{1} - l_{2})}, Δf_{k}^{(l_{1} - l_{2})})) .

Δ r_{(l_{1}, l_{2})} = ∣ r_{l_{1}} ∣ - ∣ r_{l_{2}} ∣

Δ r_{(l_{1}, l_{2})} = ∣ r_{l_{1}} ∣ - ∣ r_{l_{2}} ∣

Δ \overset{r}{˙}_{(l_{1}, l_{2})} = ∣ \overset{r}{˙}_{l_{1}} ∣ - ∣ \overset{r}{˙}_{l_{2}} ∣

Δ \overset{r}{˙}_{(l_{1}, l_{2})} = ∣ \overset{r}{˙}_{l_{1}} ∣ - ∣ \overset{r}{˙}_{l_{2}} ∣

x = s_{l_{1}} - ∣ r_{l_{1}} ∣ [cos (α - α_{0}) sin (α - α_{0})]

x = s_{l_{1}} - ∣ r_{l_{1}} ∣ [cos (α - α_{0}) sin (α - α_{0})]

∣ r_{l_{2}} ∣^{2}

∣ r_{l_{2}} ∣^{2}

= ∣ r_{l_{1}} ∣^{2} + B^{2} + 2 B ∣ r_{l_{1}} ∣ cos (α) .

∣ r_{l_{1}} ∣ = \frac{Δ r _{(l_{1}, l_{2})}^{2} - B ^{2}}{2 ( Δ r _{(l_{1}, l_{2})} + B cos ( α ))} .

∣ r_{l_{1}} ∣ = \frac{Δ r _{(l_{1}, l_{2})}^{2} - B ^{2}}{2 ( Δ r _{(l_{1}, l_{2})} + B cos ( α ))} .

α = cos^{- 1} (\frac{Δ r _{(l_{1}, l_{2})}^{2} - B ^{2} - 2Δ r _{(l_{1}, l_{2})} ∣ r _{l_{1}} ∣}{2 B ∣ r _{l_{1}} ∣})

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Passive Multi-Target Tracking

Using the Adaptive Birth Intensity PHD Filter

Christopher Berry1, Donald J. Bucci1, Samuel Watt Schmidt1

1Lockheed Martin Advanced Technology Laboratories, Cherry Hill, NJ, USA

Email: {christopher.m.berry, donald.j.bucci.jr, samuel.w.schmidt}@lmco.com

Abstract

Passive multi-target tracking applications require the integration of multiple spatially distributed sensor measurements to distinguish true tracks from ghost tracks. A popular multi-target tracking approach for these applications is the particle filter implementation of Mahler’s probability hypothesis density (PHD) filter, which jointly updates the union of all target state space estimates without requiring computationally complex measurement-to-track data association. Although this technique is attractive for implementation in computationally limited platforms, the performance benefits can be significantly overshadowed by inefficient sampling of the target birth particles over the region of interest. We propose a multi-sensor extension of the adaptive birth intensity PHD filter described in (Ristic, 2012) to achieve efficient birth particle sampling driven by online sensor measurements from multiple sensors. The proposed approach is demonstrated using distributed time-difference-of-arrival (TDOA) and frequency-difference-of-arrival (FDOA) measurements, in which we describe exact techniques for sampling from the target state space conditioned on the observations. Numerical results are presented that demonstrate the increased particle density efficiency of the proposed approach over a uniform birth particle sampler.

I Introduction

Passive localization and tracking of multiple co-channel targets is an important research topic in speaker localization [1], sonar [2], passive radar [3], and intrusion detection in wireless sensor networks [4]. By passive, we mean that the signal shape emitted by the target is unknown, but observed similarly across multiple spatially distributed sensors. As a result, the localization and tracking algorithms for these applications rely on difference measurements between pairs of sensors. Time-difference-of-arrival (TDOA) and frequency-difference-of-arrival (FDOA) measurements, for example, use the generalized cross correlation of raw signal measurements between sensor pairs at different time and frequency offsets to estimate range and range-rate differences [5, 6]. Power-difference-of-arrival (PDOA) measurements are generated by observing the difference in average receive power measurements between sensors, which is closely related to the difference in the cross correlation peak magnitudes. In applications such as underwater acoustics and passive radar a direction finding solution based on angle-of-arrival (AOA) is also possible [2].

Although single-target localization and tracking algorithms using passive measurements are well understood [7], the problem becomes significantly more complex for multiple co-channel targets. For passive localization and tracking algorithms in particular, the presence of multiple targets leads to state space estimation ambiguities known as ghost targets. These ghost targets manifest from the intersection of plausible state space estimates suggested by each passive sensor measurement [8]. Resolution of these ghost targets requires specification of measurement-to-track association algorithms that in general do not scale well with the number of targets, especially in the presence of high clutter and high false alarm rate environments.

A number of multi-target tracking approaches have been proposed for the passive multi-target tracking problem based on the random finite sets (RFS) framework [9, 10]. These techniques extend the calculus and statistics used to derive the standard recursive Bayes filter to the case where the cardinality of the underlying state space is itself a random variable. The resulting multi-target recursive Bayes filter provides a mechanism for estimating the number of targets and their state space estimates without having to specify a measurement-to-track association technique. Although this filter is practically difficult to implement, a number of moment-based approximations have been proposed that are widely used in computationally limited settings. The probability hypothesis density (PHD) filter results from a first-order moment approximation where the cardinality distribution is assumed Poisson. The cardinalized probability hypothesis density (CPHD) filter results from a first-order moment approximation where the cardinality distribution is estimated in parallel. Both of these approximation techniques are defined for single sensor tracking problems, as the original multi-sensor update equation requires the combinatorial evaluation of all false-alarm and missed detection hypotheses per sensor. In practice, the iterated-corrector PHD/CPHD method is used which involves applying the update step sequentially to the set of measurements from each sensor [11, 3]. Although this approach is suboptimal [12, 13], it remains practical in many applications due to the marginal increase in complexity.

In this paper, we focus on the particle filter implementation of the iterated-corrector PHD filter [14, 11, 3] applied to the passive multi-target tracking problem using a stationary wireless sensor network. A set of spatially distributed sensors observe co-channel, non-overlapping (in time) signals from a time varying number of targets. Each sensor observes pairs of TDOA and FDOA measurements that are generated from at most one of the targets, or the result of environmental clutter. The goal is to use the measurements at each time step to track the position and velocity of all targets. In this work, we accomplish this task using the adaptive birth intensity extension to the PHD filter detailed in [15]. This approach constructs an augmented state space consisting of the original state-space variables concatenated with a label that distinguishes a newborn target from a persistent target. The resulting predict and update steps of the PHD filter after applying this augmentation removes the need for prior specification of the target birth intensity, and instead replaces it with a particle sampling procedure based on the measurements themselves. More importantly, this data-adaptive sampling technique enables very accurate multi-target state estimation without requiring an excessive number of particles. This consideration is important, especially when considering implementation on computationally limited software defined radio (SDR) platforms.

In order to implement this technique for the TDOA/FDOA passive localization problem, we provide the following contributions. First, we describe an extension of the adaptive birth intensity technique provided in [15] to multiple sensor measurements using the iterated-corrector implementation of the particle PHD filter. The approach of [15] provides a more rigorous justification for an adaptive birth intensity, in contrast to similar heuristics previously suggested in [16, 17]. As opposed to the linear Gaussian approximation of the birth intensity sampler suggested in [15], we also provide a mechanism for directly sampling the plausible state-space particles using the TDOA and FDOA measurement tuples. Our technique is based on the geometric models presented previously in [18]. The specification of these sampling functions is very important, as the relationship between the state space estimates and the TDOA and FDOA measurements are non-linear and non-invertible.

The remainder of the paper is organized as follows. In Section II, a brief overview of RFS for multi-target tracking is provided. The adaptive birth intensity particle PHD filter from [15] alongside our extension to multi-sensor measurements using the iterated-corrector update rule is summarized in Section III. We then detail the TDOA and FDOA measurement models alongside the relevant position and velocity sampling techniques in Section IV. Finally, Section V provides a simulation comparing the performance of the proposed technique with a standard uniform placement of birth particles.

II Random Finite Set Preliminaries

A full treatment of the multi-target tracking problem with RFS is outside the scope of this paper. Instead, we provide a pragmatic description using the notation provided in [15]. A detailed description of RFS can be found in [19], and its application to the particle PHD filter in [14].

II-A Random Finite Sets and Multi-target Tracking

Consider a multi-target tracking problem where the goal is to estimate the state vector, $\mathbf{x}_{k}\in\mathscr{X}$ , for $N_{k}$ targets at the discrete time step $k=0,1,\ldots$ . At each time step, the number of targets is random, and the states for each target evolve according to a Markov process. The multi-target state at time $k$ is defined as

[TABLE]

where $\mathscr{F}(\mathscr{X})$ denotes a finite collection of subsets defined on $\mathscr{X}$ . A single sensor generates a set of measurements, $\mathbf{z}_{k}\in\mathscr{Z}$ , for each target in addition to false alarm measurements (i.e., clutter) that are uncorrelated with the any of the target states. In addition, missed detections may occur in which no measurements are observed from a given target. The total set of $m_{k}$ measurements at time step $k$ , consisting of true target and clutter returns, is defined as

[TABLE]

Both $\mathbf{X}_{k}$ and $\mathbf{Z}_{k}$ are described as random finite sets (RFS), where the values of the elements and the cardinality of the sets are random variables with no specific element permutation (i.e., ordering of the targets). In a manner similar to single-target filtering, a RFS multi-target filter seeks to estimate the posterior multi-target state distribution, $f_{k|k}(\mathbf{X}_{k}|\mathbf{Z}_{1:k})$ , which relates the multi-target state estimates $\mathbf{X}_{k}$ given the measurements up until time $k$ , denoted $\mathbf{Z}_{1:k}=\{\mathbf{Z}_{1},\ldots,\mathbf{Z}_{k}\}$ . This RFS posterior is fully specified using the multi-target recursive Bayes filter derived in [19].

II-B Probability Hypothesis Density Filtering

Although particle filter implementations of the multi-target recursive Bayes filter exist [1], they are not often used in practice because their complexity becomes intractable for large numbers of targets. Instead, a number of RFS moment approximations are applied to arrive at a suboptimal but more tractable realization. The moment approximation of interest for this work is the probability hypothesis density (PHD). The PHD, denoted $D(\mathbf{x}):\mathscr{X}\to\Re$ , defines a multi-target analogue to an expected value in the sense that it provides a point-wise estimate of the density of targets located at $\mathbf{x}_{k}\in\mathscr{X}$ [15, 19]. As a consequence, the standard integral of the PHD function across the state space $\mathscr{X}$ evaluates to the expected number of targets in the state space.

In a multi-target tracking system the PHD function can be loosely interpreted as an estimate of the union (i.e., weighted summation) of single-target tracker posteriors and a prior density function consisting of a birth probability distribution over the state space multiplied by the expected number of newborn targets. Denote the most recent estimate of the PHD intensity function at time $k-1$ as $D_{k-1|k-1}(\mathbf{x})$ . The PHD filter prediction step is given by

[TABLE]

where $p_{s}(\mathbf{x})$ is the probability of survival for a target, $\pi_{k|k-1}(\mathbf{x}|\mathbf{x}^{\prime})$ is the standard state-transition probability density function that describes how target kinematics propagate from time step $k-1$ to $k$ , and $\gamma_{k|k-1}(\mathbf{x})$ is a PHD function representing the birth process (i.e., the expected number of targets that appear at each state space location). The PHD filter update step is given by

[TABLE]

where $p_{d}(\mathbf{x})$ is the probability of detection of a measurement being generated from a target, $g_{k}(\mathbf{z}_{k}|\mathbf{x})$ is the standard single-target measurement likelihood, $\kappa_{k}(\mathbf{z}_{k})$ is a PHD function representing the measurement clutter process (i.e., the expected number of measurement false alarms that appear at each measurement space location). The normalization term $C(\mathbf{z}_{k})$ is given as

[TABLE]

Since the PHD filter involves multiple standard integrals over the state space, it is typically implemented using a sequential Monte Carlo (i.e., particle filter) realization [14], or using Gaussian mixture distributions [20]. When multiple sensors are present, the iterated-corrector approach applies the update step of Equation (2) sequentially for the set of measurements obtained at each sensor.

III Multi-sensor Adaptive Birth Intensity Filter

III-A Adaptive Birth Intensity PHD Filter

The adaptive birth intensity technique of [15] augments the original state space by appending a label, $\beta\in\{0,1\}$ , which distinguishes a newborn target ( $\beta=1$ ) from a persisting target ( $\beta=0$ ). Augmenting the original state vector of the PHD predict and update equations with this label decouples the birth PHD intensity from the predicted PHD intensity of persisting targets. The PHD filter predict step for this augmented state space, denoted $D_{k|k-1}(\mathbf{x},\beta)$ is given in [15] as

[TABLE]

and

[TABLE]

The PHD update step for this augmented state space is

[TABLE]

and

[TABLE]

where the new normalization term $C^{\prime}(\mathbf{z}_{k})$ is given as

[TABLE]

It is important to note that the addition of these target labels do not change the mechanics of the PHD filter, but instead provide insight as to how the birth PHD function affects the predict and update steps. Newborn targets from the birth process are brought into the updated PHD function based on their product with the measurement likelihood function. As discussed in [15], this insight allows for construction of a particle filtering solution where the birth intensity is set such that the particle density is high at areas of high measurement likelihood (i.e., $\int g_{k}(\mathbf{z}_{k}|\mathbf{x}^{\prime})\gamma_{k|k-1}(\mathbf{x}^{\prime})d\mathbf{x}^{\prime}$ is large).

The particle filter implementation of the adaptive birth intensity PHD filter for a single sensor is specified by constructing two particle systems: one for persisting targets and one for newborn targets. These approximations are expressed as

[TABLE]

where the sets $\cup_{i=1}^{N^{p}_{k}}(\mathbf{x}_{k,p}^{(i)},w_{k,p}^{(i)})$ and $\cup_{i=1}^{N^{b}_{k}}(\mathbf{x}_{k,b}^{(i)},w_{k,b}^{(i)})$ represent the particle systems for the persistent and newborn targets respectively. The Dirac delta function $\delta_{\mathbf{x}^{\prime}}(\mathbf{x})=1$ when $\mathbf{x}=\mathbf{x}^{\prime}$ and 0 otherwise. An additional particle system is defined to represent the sum of the persistent and newborn PHD intensity functions used in the predict step of Equation (4),

[TABLE]

which is defined as the union of the persistent and newborn particle systems,

[TABLE]

where $N_{k}=N^{b}_{k}+N^{p}_{k}$ .

With these three particle system approximations, the particle filter variation of the single target adaptive birth intensity PHD filter is given as follows [15]. In the predict step, the persisting particle weights are determined by sampling all $\mathbf{x}_{k,p}^{(i)}$ from an importance density function $q_{k}(\mathbf{x}^{(i)}_{k,p}|\mathbf{x}^{(i)}_{k-1},\mathbf{Z}_{k})$ and setting the weights according to

[TABLE]

The newborn particles and their weights are determined by sampling $\mathbf{x}_{k+1,b}^{(i)}$ from some measurement conditioned probability distribution function $b(\mathbf{x}_{k+1,b}|\mathbf{z}_{k})$ for each $\mathbf{z}_{k}\in\mathbf{Z}_{k}$ . The probability density function $b(\mathbf{x}_{k+1,b}|\mathbf{z}_{k})$ is selected such that all particles generated are in regions of space where the likelihood $g_{k+1}(\mathbf{z}_{k}|\mathbf{x}_{k+1,b})$ is high. The corresponding weights for the newborn particles are uniform, and given as

[TABLE]

where $\nu^{b}_{k|k-1}$ is the expected number of target births at each predict step. In the update step, the predicted persistent particle weights are refined as

[TABLE]

where the normalization constant is given as

[TABLE]

Similarly, the birth particle weights are refined from the predicted birth particle set according to

[TABLE]

Once the update step is finished, a resampling step is applied to the persistent111As [15] suggests, the resampling step can also be applied for the updated newborn particles, if necessary. In this work, we apply resampling only to the updated persistent particle systems. PHD particle systems in order to remove samples with low weights. The usual approach is to compute the weight sum of each system, normalize the particle weights by each system’s weight sum, resample the points using each particle system’s normalized weights, and finally set the new particle weights uniformly such that they sum to the original weight sums of each system. A more detailed description is provided in [14].

III-B Multi-sensor Adaptive Birth Intensity PHD Filter

Given the particle filter implementation of the single-sensor adaptive birth intensity PHD filter [15], our extension to multi-sensor based on the iterated-corrector rule is straight forward. The pseudocode for the resulting multi-sensor adaptive birth intensity PHD filter is shown in Figure 1. Each step of the filter operates on the combined persistent (after resampling) and newborn particle system from the previous step, using the measurements generated by each sensor $l\in{1,\ldots,L}$ . The predict step (lines 1-5) is similar to what is described in [15], and based on target survival and kinematic propagation. Applying the iterated-corrector to the update step, the proposed approach samples a particle system using a birth intensity based on the measurements observed for a given sensor (lines 6-23). The predicted persistent and newborn particle systems are updated with that sensor’s observed measurements, denoted $\mathbf{Z}_{k}^{(l)}$ , using the adaptive birth intensity PHD update steps defined by Equations (11) and (13) (lines 17-18). Then, an intermediate “predict” step occurs that converts all of the birth particles generated by that sensor (after performing the update of Equation (13)) to persisting particles (line 19-22). The next sensor’s birth particle system is sampled using its measurement set, and the update equations are repeated using the persistent particle set generated after the previous sensor’s update.

It should be noted that the iterated-corrector multi-sensor combination rule is not commutative. As a result, sensors that are combined earlier in the update step will have a larger effect on the resulting persistent particle system than those combined later. For passive tracking applications, this could have the effect of introducing persistent ghost targets in the filter output if sensors with higher false alarm rates are combined first. A common strategy used in passive radar applications is to randomize the combination order, favoring earlier selection of sensors with lower expected false alarm rates and higher detection rates. It should also be noted that the complexity of the proposed iterator-corrector update step grows with the number of measurements detected at each of the sensors, which can be very large for high false alarm rate sensors. However, this complexity growth is largely offset by the increase in particle placement efficiency provided by this technique over a uniform placement of particles.

IV Particle Sampling from Passive Measurements

We now describe the application of the multi-sensor adaptive birth intensity filter to a passive multi-target target tracking problem using TDOA and FDOA measurements. The TDOA and FDOA geometry for a single sensor pair and a single target is shown in Figure 2. Let

[TABLE]

be the measurement RFS for TDOA measurements generated at time-step $k$ for sensor pair $(l_{1},l_{2})$ . Similarly, let

[TABLE]

be the measurement RFS for FDOA measurements generated at time-step $k$ for sensor pair $(l_{1},l_{2})$ . The notation $(l_{1}-l_{2})$ represents that the measurement difference is taken with respect to sensor $l_{2}$ . We assume that the signal-level detection mechanism provides joint TDOA and FDOA detections that are independent and normally distributed for each true target such that

[TABLE]

and

[TABLE]

where $\sigma^{2}_{\Delta t}$ and $\sigma^{2}_{\Delta f}$ represent the expected TDOA and FDOA error variances, $c$ is the speed of light, and $f_{c}$ is the carrier frequency of the observation. The quantities $|r_{l}|$ and $|\dot{r}_{l}|$ represent the range and range-rate from sensor $l$ to the target, where we have dropped the dependence on the location of the target for convenience. Because the measurements are obtained jointly, their detections are characterized by a single TDOA/FDOA correct detection and false alarm rate probability per sensor. In the following sections, we describe a measurement driven birth particle sampling procedure using TDOA/FDOA measurements

[TABLE]

Without loss of generality, we consider the sampling technique for a single target at a fixed time step and for a single pair of sensors in the following subsections, denoting

[TABLE]

and

[TABLE]

for convenience.

IV-A Birth Particle Sampling from TDOA Measurements

The $x,y$ coordinates of the target assuming the coordinate system given in Figure 2 is given as

[TABLE]

where the quantities $|r_{l_{1}}|$ and $\alpha$ must be sampled in order to be consistent with the observed TDOA value (i.e., the range difference $\Delta r_{({l_{1}},{l_{2}})}$ ). Applying law of cosines to Figure 2, we have that

[TABLE]

Noting that $|r_{l_{2}}|=|r_{l_{1}}|-\Delta r_{({l_{1}},{l_{2}})}$ , the relationship between $|r_{l_{1}}|$ and $\alpha$ is given as [18]

[TABLE]

where $|r_{l_{1}}|\geq(\Delta r_{({l_{1}},{l_{2}})}+B)/2$ . This lower bound is achieved when the target is on the line that intersects sensors $l_{1}$ and $l_{2}$ (i.e., $\alpha=\pi$ ). In order to sample points uniformly along the hyperbola defined by Equations (18) and (19), we first use the observed TDOA measurement to sample an estimate of the true range difference measurement, $\Delta r_{({l_{1}},{l_{2}})}$ , according to the expected measurement variance and Gaussian assumption of Equation (14). This sampled TDOA value is used to sample $|r_{1}|$ uniformly between $(\Delta r_{({l_{1}},{l_{2}})}+B)/2$ and some maximum value that is set based on the detection range of the sensor or the bounding box that defines the region of interest. Equation (19) is then solved for $\alpha$ to yield

[TABLE]

In this construction, the corresponding range on $\alpha$ covers the lower half of the hyperbola (i.e., positions beneath the line connecting the two sensors in Figure 2). To additionally consider the upper half of the hyperbola, we invert the sign of the inverse cosine argument in Equation (20) with probability $0.5$ . The sampled values of $|r_{1}|$ and $\alpha$ are substituted into Equation (18) to generate the corresponding particle positions in the desired $x,y$ coordinate frame.

IV-B Birth Particle Sampling from FDOA Measurements

Given a particle location generated from a TDOA measurement, the angles $\alpha$ and $\beta$ can be determined uniquely by applying law of cosines. Using Figure 2, the observed range-rate difference of Equation (17) for that particle is expressed as

[TABLE]

where the final step results by noting that the linear combination of two cosines with arbitrary phases is a sinusoid with a scaled magnitude and phase shift. The resulting FDOA relationship between the velocity magnitude $|v|$ and the platform heading $\theta$ is then given as

[TABLE]

where

[TABLE]

with equality when

[TABLE]

In order to sample a velocity magnitude and bearing for each position particle generated by the TDOA measurements, we use the observed FDOA measurement to sample an estimate of the true range-rate difference, $\Delta\dot{r}_{(l_{1},l_{2})}$ , according to the expected measurement variance and Gaussian assumption of Equation (15). A velocity magnitude is uniformly sampled between the lower bound defined in Equation (23) and some maximum possible velocity estimate. The corresponding heading of the target is then found by solving Equation (22) for $\theta$ to yield

[TABLE]

Similar to the TDOA case, we invert the sign of the inverse sine argument in Equation (24) with probability $0.5$ in order to evenly samples headings with values greater than and less than $\frac{\alpha-\beta}{2}$ . The final velocity vector expressed in the $x,y$ coordinate system of Figure 2 can then be fully determined as

[TABLE]

V Tracking Simulation

V-A Simulation Setup

In this section we numerically compare the tracking performance of the proposed TDOA/FDOA adaptive birth intensity PHD filter (denoted PHDF-M) versus a uniform birth intensity PHD filter of comparable particle density (denoted PHDF-U). Figure 3 shows the transmitter trajectories and receiver locations in a 2km $\times$ 2km area of interest. The velocity magnitude of all targets was fixed at 15 m/s. The state vector of each target consisted of position and velocity in $x$ and $y$ coordinates (i.e., $\mathbf{x}=[x,\dot{x},y,\dot{y}]^{\text{T}}$ ). A constant velocity transitional density, $\pi_{k|k-1}(\mathbf{x}|\mathbf{x^{\prime}})=\mathcal{N}(\mathbf{x};\mathbf{Fx^{\prime}},\mathbf{Q})$ , was assumed for persisting targets where

[TABLE]

$T$ is the sampling interval, $\otimes$ is the Kronecker product, $\mathbf{I}_{m}$ is the $m\times m$ identity matrix, and $q$ is the process noise intensity. The simulation used $T=1$ second and $q=0.3$ . The probability of survival was set to $p_{s}(\mathbf{x})=.98$ .

TDOA and FDOA measurements were obtained from all sensor pair combinations $\{(0,1),(0,2),(0,3),(1,2),(1,3),(2,3)\}$ . As described in Equations (14) and (15), the measurement likelihood function for each sensor pair was set to $g_{k}(\mathbf{z}^{(l_{1}-l_{2})}|\mathbf{x})=\mathcal{N}(\mathbf{z}^{(l_{1}-l_{2})};\mathbf{h}^{(l_{1}-l_{2})}(\mathbf{x}),\mathbf{R})$ , where

[TABLE]

and $\mathbf{R}=\text{diag}(\sigma^{2}_{\Delta t},\sigma^{2}_{\Delta f})$ . The notations $|r_{l}(\mathbf{x})|$ and $|\dot{r}_{l}(\mathbf{x})|$ denote the range and range-rate from sensor $l$ to the target located in the state space at $\mathbf{x}$ . The timing and spectral measurements standard deviations were set as $\sigma_{\Delta t}=20\text{ns}$ and $\sigma_{\Delta f}=2.5\text{Hz}$ , respectively. The clutter was uniformly distributed in time from $-cB^{(l_{1}-l_{2})}$ to $cB^{(l_{1}-l_{2})}$ and frequency from $-2|v|f_{c}/c$ to $2|v|f_{c}/c$ , where $|v|=25\text{m/s}$ and $f_{c}=2.4\text{GHz}$ . Recall from the previous section that the quantity $B^{(l_{1}-l_{2})}$ represents the straight-line distance between sensors $l_{1}$ and $l_{2}$ . The number of clutter points per receive window was Poisson distributed with a mean value of $\lambda=2$ . For reference, a single realization of the clutter and true target measurements under this value of $\lambda$ is shown in Figure 4. The probability of detection was set as $p_{d}(\mathbf{x})=0.99$ . For every measurement $\mathbf{z}^{(l_{1}-l_{2})}\in\mathbf{Z}^{(l_{1}-l_{2})}_{k}$ , newborn particles were birthed using the method described in Section IV, with $|r_{l_{1}}|\leq 2000\text{m}$ , $|v|\leq 25\text{m/s}$ , and $M_{b}=500\text{ particles/measurement}$ . To ensure a fair comparison between filters, the particle density of the PHDF-U was fixed to be the same as the PHDF-M density when all targets were present. Finally, the parameter $\nu_{k|k-1}^{b}=0.0001$ was selected so that the average number of newborn targets per scan was $\hat{\nu}_{k}^{b}\approx 0.25$ .

V-B Results

To capture median performance, $100$ Monte Carlo (MC) realizations of the tracks shown in Figure 3 were simulated. State extraction was performed using the method described in [21], where the persistent particle system before resampling was used to estimate per particle measurement association weights and target existence probabilities. The optimal subpattern assignment (OSPA) metric [22] was used to compare performance of the two filters. Figure 5 shows the median OSPA metric alongside the OSPA distance and cardinality components at each time step for the PHDF-M and PHDF-U filters. The OSPA cutoff parameter was set to $c=20\text{m}$ . The PHDF-M cardinality performance significantly outperforms the PHDF-U. The PHDF-M was able to track all targets while the PHDF-U was only able to consistently track one of the three targets over the duration of the entire experiment. For all tracked target states, the magnitude of the OSPA distance component of the PHDF-M and PHDF-U was similar. To achieve equivalent OSPA cardinality and distance performance the number of particles of the PHDF-U had to be increased substantially.

VI Conclusions and Future Work

We have presented a technique for passive TDOA/FDOA tracking of multiple targets using the adaptive birth intensity PHD filter described in [15]. A multi-sensor extension of adaptive birth intensity PHD filter using the iterated-corrector technique was described alongside an exact technique for performing the state-space sampling of newborn particles directly from TDOA and FDOA measurements. This technique was shown to adequately track the number and state of multiple targets with significantly less particles as compared to the uniform birth sampling technique typically used with the PHD filter. Future work will involve the evaluation of the proposed tracking algorithm on passive measurements collected from a wireless sensor network of software defined radios. Additionally, we will also be investigating sampling techniques for PDOA and extensions to the product multi-sensor fusion rule.

Bibliography22

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] W.-K. Ma, B.-N. Vo, S. Singh, and A. Baddeley, “Tracking an unknown time-varying number of speakers using TDOA measurements: A random finite set approach,” IEEE Transactions on Signal Processing , vol. 54, no. 9, pp. 3291–3304, 2006.
2[2] Q. Zhang, Y. Xie, and T. L. Song, “Distributed multi-target tracking in clutter for passive linear array sonar systems,” in 20th International Conference on Information Fusion (Fusion) . IEEE, Jul 2017, pp. 1–8.
3[3] S. Tian, Y. He, and G. Wang, “PHD Filter of Multi-target Tracking With Passive Radar Observations,” in 8th international Conference on Signal Processing . IEEE, 2006.
4[4] A. Vempaty, Y. S. Han, and P. K. Varshney, “Target localization in wireless sensor networks using error correcting codes,” IEEE Transactions on Information Theory , vol. 60, no. 1, pp. 697–712, Jan 2014.
5[5] N. Vankayalapati and S. Kay, “Asymptotically optimal detection of low probability of intercept signals using distributed sensors,” IEEE Transactions on Aerospace and Electronic Systems , vol. 48, no. 1, pp. 737–748, Jan 2012.
6[6] C. Holt, “Two-channel likelihood detectors for arbitrary linear channel distortion,” IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 35, no. 3, pp. 267–273, Mar 1987.
7[7] N. Vankayalapati, S. Kay, and Quan Ding, “TDOA based direct positioning maximum likelihood estimator and the Cramer-Rao bound,” IEEE Transactions on Aerospace and Electronic Systems , vol. 50, no. 3, pp. 1616–1635, Jul 2014.
8[8] A. El Gonnouni and F. Lehmann, “Deghosting method for multiple target tracking in a single frequency network,” in International Radar Conference . IEEE, Oct 2014, pp. 1–6.