1

The Bin-Occupancy Filter and its Connection to the PHD Filters Ozgur Erdinc, Peter Willett and Yaakov Bar-Shalom ECE Department, University of Connecticut {ozgur,willett,ybs}@engr.uconn.edu

Abstract—An algorithm that is capable not only of tracking multiple targets but also of “track management” – meaning that it does not need to know the number of targets as a user input – is of considerable interest. In this paper we devise a recursive trackmanaged filter via a quantized state-space (“bin”) model. In the limit, as the discretization implied by the bins becomes as refined as possible (infinitesimal bins) we find that the filter equations are identical to Mahler’s probability hypothesis density (PHD) filter, a novel track-managed filtering scheme that is attracting increasing attention. Thus one contribution of this paper is an interpretation of, if not the PHD itself, at least what the PHD is doing. This does offer some intuitive appeal, but has some practical use as well: with this model it is possible to identify the PHD’s “target-death” problem, and also the statistical inference structures of the PHD filters. To obviate the target death problem, PHD originator Mahler developed a new “cardinalized” version of PHD (CPHD). The second contribution of this paper is to extend the “bin-occupancy” model such that the resulting recursive filter is identical to the cardinalized PHD filter. Index Terms—Probability hypothesis density filter, PHD, cardinalized, CPHD, multi-target tracking, MTT. EDICS: SSP-TRAC, SAM-SONR, SAM-RADR.

- p(i|xj ): probability that a target located at xj moves to bin i. - f (xi |xj ): probability density that target moves from xj to xi . - Z1k−1 : all measurements from time 1 up to and including time k − 1. - zs : sth (single) measurement at time k (we suppress the notation k). - Zk : the set of measurements received at time k. Zk = {z1 , z2 , . . . , zm }. - λ: false alarm density. - f (Zk ): joint probability density function of measurements at time k, i.e. f (z1 , z2 , . . . , zm ). - f (z): probability density2 function (pdf) of single measurement z at time k. - pk−1|k−1 (Ui ): P{Uk−1 (i) = 1|Z1k−1 } – Dk−1|k−1 (xi ) is the limiting density, as |νi | → 0, of pk−1|k−1 (Ui ), through (??). - pk|k−1 (Ui ): P{Uk (i) = 1|Z1k−1 }) – Dk|k−1 (xi ) is the limiting density, as |νi | → 0, of pk|k−1 (Ui ), through (??).

A. Notation Our notation is the following: -

d: the dimension of the state vector. S: surveillance region. |V |: the volume of the surveillance region. |V | < ∞. xi : the middle point1 (location) of bin i in the state space. xi ∈ S. νi : the set of all points inside bin i. |νi |: volume of bin i, identical for all i. Uk (i): indicator that bin i contains a target at time k, Uk (i) ∈ {0, 1}. Vk (i): indicator that an extant target in bin i is detected, Vk (i) ∈ {0, 1}. Pd (xi ): probability of detection of a target located at xi . Ps (xi ): probability of survival of a target located at xi . P{A}: probability that event A occurs. b(i) : “birth” probability, that bin i, which did not contain a target on the prior scan, now does – b(xi ) is the limiting density, as |νi | → 0, of b(i), through (??).

This manuscript, an expanded and improved version of [?], was supported by the Office of Naval Research under contracts N00014-07-1-0429 and N00014-07-1-0055, and by ARO. 1 By definition. Any point inside bin i could be chosen as x . i

I. INTRODUCTION The multitarget tracking (MTT) problem [?], [?], [?] aims to estimate the states of an unknown number of targets by processing of measurements. Many of the well-established MTT algorithms attack the problem in two phases, as by answering the question: “Assuming there are N targets, where are they?” followed by “What is N ?”. In this paper, we propose a multitarget tracking algorithm that approaches the problem by asking a different question: “Is there a target at a given point?”. Our approach is via a model of the surveillance region using small “bins” wherein a target can (or may not) lie. Then, by considering the limiting case in which the volume of these bins goes to zero, we reach the continuous form of this “bin-occupancy” filter. The resulting prediction and update equations for the binoccupancy filter will be immediately recognizable as that of the probability hypothesis density (PHD) filter, which was originally proposed by Mahler in a sequence of papers (e.g. [?], [?], [?], [?]) and is attracting increasing attention3 [?], R 2 We use f (.) to denote the probability density function of its argument, f (γ)dγ = 1. a comprehensive literature survey, the reader may refer to [?].

3 For

2

[?], [?], [?], [?], [?], [?], [?]. The PHD is a surface whose integration over any region yields the expected number of targets in that region [?]. In its full generality the PHD requires a full particle-based update [?], [?], [?], [?]. However, in the important special case of a linear/Gaussian model there is an explicit update [?], [?]; we take from this that the PHD, already elegant, has become practical as well. Our derivation of the bin-occupancy filter is based on a physical bin model: it does not use tools from finite point processes (FPP) as does [?], nor does it require set-derivatives or higher-mathematical concepts from, for example, [?], and it needs no approximations. We also must stress that because two different sets of assumptions or criteria suggest the same procedure, it does not mean that these assumptions or criteria are equivalent. An example of this is MAP estimation in linear/Gaussian problems, which has the same Kalman result as LMMSE estimation in more general situations; but it is clear that there is insight to be gained by looking at one from the other’s perspective. The development of the PHD by Mahler is elegant, correct and really needs no comment by us, but we hope that a contribution here is that the proposed (algorithmically-equivalent) bin-occupancy filter would increase understanding and accessibility of the PHD filter amongst engineers who might prefer a simpler physical model from which the development is easier to grasp. As pointed out in a recent conference paper [?], the PHD filter tends to exaggerate the impact of missed observations, and declares target-death prematurely. This leads to an unattractive behavior of the estimate of the number of targets, and, as a solution to this problem, Mahler proposed a modified version, the cardinalized PHD filter [?], [?]. The mathematics behind the CPHD is of even greater depth than that supporting the original PHD. In the second part of this paper, we derive a modified version of the “bin-occupancy” filter that results in the CPHD filter equations. The organization of the paper is as follows. Section ?? describes the bin model and gives the derivation of the binoccupancy filter: the filter itself is the same as Mahler’s PHD filter, just derived directly from a physical model4 . When one has a (physical) model, one can discuss, criticize and perhaps improve it. For example, in section ?? we discuss an approximation that underlies the bin occupancy filter, an approximation that may be relevant to the PHD as well. Further, our ability to describe the PHD in physical terms led us to identify its target-death problem: that is discussed in section ?? together with the bin model derivation of the CPHD filter. We note that this manuscript shares many of the conclusions that can be found in [?]; however, in that conference paper the mathematical development was both expurgated and in some cases a little flaccid; this manuscript is, we hope, both more complete and more correct.

4 A part of this work, based on [?], can be found in Mahler’s recent book [?]. This manuscript has significantly improved since the original publication, most significantly that the physical bin-model is more exact in the current version.

II. T HE BIN - OCCUPANCY FILTER We assume the surveillance region is partitioned into bins such that each bin (i.e., partition) has the same volume. Moreover, these bins are sufficiently small that each bin is potentially occupied by at most one target. By “sufficiently small” we mean the following: regardless of the dimension of the targets, if two targets happen to be very closely located such that one bin is occupied by both of them, one can always reduce the bin volume accordingly so that the assumption of one target per bin remains valid. Our approach is to devise a

Fig. 1. Illustration of the model of “infinitesimal” bins: the target in bin j moves to bin i with probability p(i|xj ). The other possibilities are that the target in bin i dies with probability 1 − Ps (xi ) and a new-born target appears at bin i.

recursive filter (prediction and update) for the bin-occupancy probabilities and, by considering the limiting case where the bin-widths go to zero (i.e. the partition gets finer and finer – see Appendix ??), to devise the continuous version of the filter. We make the following assumptions: 1) The bins are sufficiently small so that each bin is occupied by at most one target5 [?]. 2) One target gives rise to only one measurement. 3) Each target generates measurements independently. 4) False alarms are independent of target originated measurements. 5) False alarms are Poisson distributed. We do not require, of course, that the dimension of the observation zs be d. Also, b(x) may be taken as intensity of a Poisson point process, one of target births (in S): the definition of b(x) will be given in (??).

A. Prediction Given knowledge of all the occupied bins at time k − 1, the event that bin i contains a target at current time k is possible in two6 ways: 1) a new-born target appears in bin i, 2) one of the targets in another bin moves to bin i (this includes the event that target in bin i stays in bin i). That is, we have the 5 In [?] page 9, the same assumption is evident in a context that describes inevitability of Poisson distribution. 6 We ignore the target spawning for simplicity. However, it is straightforward to include the spawning model into this derivation. The equations that ensue when spawn is considered are exactly those of the PHD and CPHD, and the reader is referred to Mahler’s work (e.g., in [?]) for those.

3

and

event “bin i contains a target” = {Uk (i) = 1} [  = “target from bin j went to bin i”

(1)

j

[

“birth at bin i”

Our first assumption above, a bin is occupied by at most one target, indicates that two or more of the events above cannot occur simultaneously. In other words, as the bin widths are chosen “sufficiently small”, the probability of more than one target occupying a bin (e.g., two existing targets move to the same bin) evaporates, and as such we can consider the events in (??) to be mutually exclusive. Hence, the bin-occupancy probabilities can be calculated as P{Uk (i) = 1|Z1k−1 } = (2) X k−1 p(i|xj )Ps (xj )P{Uk−1 (j) = 1|Z1 } b(i) + j

where the second term in the right hand side of (??) is via the total probability theorem. The target in bin j at time (k−1) exists at time k with probability Ps (xj ), and it moves to bin i with probability p(i|xj ). The bin-occupancy density, Dk−1|k−1 (xi ), is defined such that ∆

Dk−1|k−1 (xi ) =

P{Uk−1 (i) = 1|Z1k−1 } |νi | |νi |→0 lim

(3)

Similarly, we define the following: P{Uk (i) = 1|Z1k−1 } |νi | |νi |→0 b(i) ∆ b(xi ) = lim |νi |→0 |νi | p(i|xj ) ∆ f (xi |xj ) = lim |νi |→0 |νi | ∆

Dk|k−1 (xi ) =

lim

(4) (5) (6)

Dk|k−1 (xi ) = b(xi ) +

f (xi |γ)Ps (γ)Dk−1|k−1 (γ)dγ (9)

The second term in (??) constitutes a Riemann sum [?] and since the function inside the sum is Riemann-integrable, convergence to integral is ensured as bin-volume, |νj |, approaches to zero. Equation (??) is the prediction step.

B. Update We use the Bayes equation to update the bin-occupancy probability of bin i. First, we define a “visibility” indicator Vk (i), a Bernoulli random variable with probability Pd (xi ), which is unity if the target in bin i is indeed detected, i.e. P{Vk (i) = 1} = Pd (xi ). We have |Zk | = number of detections at time k X = Uk (j)Vk (j) + number of false alarms (10) j

where the number of false alarms is assumed Poisson with parameter µF A . The bin-occupancy events, Uk (i), are Bernoulli random variables with non-identical probabilities (from bin to bin). As the volume of the bin goes to zero and the number of bins P goes to infinity, a Poisson approximation with parameter j pk|k−1 (Uj )Pd (xj ) for the first term in the sum is adequate [?], if one assumes that the bin-occupancy events are independent. This assumption has important consequences as we discuss them in the subsequent sections, most importantly section ?? reveals our intuition. Then, the number of measurements is Poisson distributed with parameter X ∆ µ= pk|k−1 (Uj )Pd (xj ) + µF A (11) j

and we assume that the surfaces defined by these quantities do not have point masses7 , and they are continuous8 . Also note that the surveillance region is bounded. We rewrite (??), by dividing both sides with |νi |, P{Uk (i) = 1|Z1k−1 } b(i) = + |νi | |νi | X p(i|xj ) P{Uk−1 (j) = 1|Z1k−1 } Ps (xj ) |νj | |νi | |νj | j

Z

i.e. |Zk | ∼ P(µ). Notice that the number of measurements is random, as well as the locations of individual measurements. Suppose that at the current time m measurements are received. Then, using Bayes equation, we can update the bin-occupancy probability of bin i by using not only the locations of the measurements but also the received number of measurements at time k,

(7)

and taking the limit where all the bin volumes go to zero P{Uk (i) = 1|Z1k−1 } b(i) = lim (8) |νi | |νi |→0,∀i |νi |→0,∀i |νi | X p(i|xj ) P{Uk−1 (j) = 1|Z1k−1 } + lim Ps (xj ) |νj | |νi | |νj | |νi |→0,∀i j

P{Uk (i)=1|Z k } 1

= =

lim

7 In point process theory, this condition is well-recognized as it defines an “orderly” process [?] page 33 or, equivalently, a “non-atomic” measure [?] p. 13. 8 This property ensures that each of these functions are Riemann-integrable since they are continuous, real-valued and bounded [?].

=

k−1 } 1 k−1 k−1 k−1 )P{Uk−1 (i)=1|Z }f (Z ) 1 1 1 k−1 k−1 f (Zk ,|Zk |=m|Z )f (Z ) 1 1 k−1 f (Zk ,|Zk |=m|Uk (i)=1,Z ) k−1 1 P{Uk−1 (i)=1|Z } k−1 1 f (Zk ,|Zk |=m|Z ) 1

P{Uk (i)=1|Zk ,|Zk |=m,Z

f (Zk ,|Zk |=m|Uk (i)=1,Z

(12) (13) (14)

Next, we derive explicit expressions for the terms in equation (??). But first, we look at the special case that there are no measurements at time k, |Zk | = m = 0: f (Zk , |Zk | = m|Uk (i) = 1, Z1k−1 )

=

(1 − Pd (xi ))e−µ (15)

f (Zk , |Zk | = m|Z1k−1 )

=

e−µ

(16)

4

since the probability of having |Zk | = 0 is e−µ and, given there is a target9 in bin i, its detection must have been missed. Plugging (??), (??) into (??), −µ

pk|k (Ui ) =

(1 − Pd (xi ))e pk|k−1 (Ui ) e−µ

(17)

By substituting (??) and (??) into (??), the update equation becomes pk|k (Ui ) 

pk|k (Ui )

=

(1 − Pd (xi ))pk|k−1 (Ui )



+

we arrive at the update equation for the case |Zk | = 0:

µm −µ k−1 )(1−Pd (xi )) e f (Zk |Z m! 1 µm−1 −µ k−1 )Pd (xi ) e f (zs |xi )f (Zk (sc )|Z m! 1 s=1

Pm



=

(23)

(18)

This equation defines the update of the bin-occupancy probabilities when there are no measurements. In the next section, we show that this result indicates what we call the “targetdeath” problem. In the general case where there are |Zk | = m observations, we can write the numerator of (??) combining the cases of that the target in bin i is respectively detected and missed. Employing the visibility Bernoulli random variable Vk (i), we have

=

 Pm pk|k−1 (Ui )

(1−Pd (xi ))+ k−1 c ) Pd (xi ) f (zs |xi )f (Zk (s )|Z1 µ k−1 s=1 f (Zk |Z ) 1

h =

f (Zk , |Zk | = m|Uk (i) = 1, Z1k−1 )

∆

f (zs |xi )f (Zk (sc )|Z1k−1 )

s=1

µm−1 e−µ m!

where (??) follows from (??) due to the vanishing bin volume: the exclusion of bin i has no effect on the likelihood of the other measurements when as the bin reduces to a singularity. We also have the same above given the target in bin i is a missed detection9 , f (Zk , |Zk | = m|Uk (i) = 1, Vk (i) = µm −µ = f (Zk |Z1k−1 ) e m!

0, Z1k−1 ) (22)

be exact, in the likelihood equations conditioned on Uk (i) = 1, the µ0 ,

∆ µ0 =

P

µ should be replaced by such that p (Uj )Pd (xj ) + j,j6=i k|k−1 µF A . However, with the “sufficiently small” bin assumption µ ≈ µ0 , and as the volume of the bins go to zero µ0 → µ. In other words, in the limiting case the approximation becomes exact.

(27)

= f (zs |Z1k−1 )

(28)

i

µ P FA c(zs )µF A + j pk|k−1 (Uj )Pd (xj ) X pk|k−1 (Ui )Pd (xi ) P + f (zs |xi ) (30) µ + FA j pk|k−1 (Uj )Pd (xj ) i  X 1 c(zs )µF A + f (zs |xi )pk|k−1 (Ui )Pd (xi ) (31) µ i

=

c(zs )

using the conditional probability of events (see Appendix ?? for the details) where the total P number of measurements are distributed with parameter µ = j pk|k−1 (Uj )Pd (xj ) + µF A . The update equation for the bin occupancy probabilities hence simplifies to pk|k (Ui ) =

 pk|k−1 (Ui )

Pm

(1−Pd (xi ))+

s=1 c(zs )µ F A+

(32) 

PPd (xi )f (zs |xi ) j

Pd (xj )f (zs |xj )pk|k−1 (Uj )

Now, by dividing both sides by |νi |, and taking the limit for all bin volumes → 0, lim|ν |→0,∀i i

pk|k (Ui ) |νi |

=lim|ν |→0,∀i i

pk|k−1 (Ui )

(33) #

|νi |

" (1−Pd (xi ))+

9 To

f (zs |xi )P(Hi (s)|Z1k−1 )

c(zs )P{zs is false alarm|Z1k−1 } X + f (zs |xi )P{zs from target in bin i|Z1k−1 }(29)

=

=

m X

X

=

f (zs |Z1k−1 )

or

=

(25)

i

f (Zk , |Zk | = m|Uk (i) = 1, Vk (i) = 1, Z1k−1 ) (20) m m−1 −µ X µ e 1 f (zs |xi )f (Zk (sc )|Z1k−1 , Hi (s)) = m (m − 1)! s=1

lim f (Zk , |Zk | = m|Uk (i) = 1, Vk (i) = 1, Z1k−1(21) )

i

to (??). pdf f (zs |Z1k−1 ) can be calculated using the total probability theorem where the spatial densities are weighted by the corresponding probabilities of the possible sources for the measurement zs :

Assuming all of these are a-priori equally likely over s, we have the joint probability density of measurements given the target in bin i is detected9 ,

|νi |→0

Pd (xi ) f (zs |xi ) µ k−1 s=1 f (zs |Zk (sc ),Z ) 1

i

lim f (zs |Zk (sc ), Z1k−1 )

(19)

Hi (s): zs originates from the target in bin i, and Zk (sc ) = {z1 , . . . , zs−1 , zs+1 , . . . , zm } originate from other sources (either from other targets or clutter) for s = 1, 2, . . . , m.

(24)



in which (??) follows from (??) since X f (zs |Zk (sc ), Z1k−1 ) = f (zs |xi )P(Hi (s)|Zk (sc ), Z1k−1 (26))

f (Zk , |Zk | = m|Uk (i) = 1, Vk (i) = 1, Z1k−1 )Pd (xi ) as in the transition of (??) +f (Zk , |Zk | = m|Uk (i) = 1, Vk (i) = 0, Z1k−1 )(1 − Pd (xi )) The single measurement

Now, we consider the following events (hypotheses):



Pm

pk|k−1 (Ui ) (1−Pd (xi ))+

|νi |→0

=

pk|k−1 (Ui )

k−1 µm −µ ) f (Zk |Z e m! 1

Pm s=1

λ+

P j

Pd (xi )f (zs |xi ) pk|k−1 (Uj )

Pd (xj )f (zs |xj )

|νj |

|νj |

The summation term in the denominator of (??) is a Riemann sum, and the function inside is Riemann-integrable since it is continuous, real-valued, and bounded [?]. Hence, as the

5

volumes of the bins go to zero, the summation term converges to an integral (see Appendix ??), so we get lim

X

|νj |→0

Pd (xj )f (zs |xj )

j

pk|k−1 (Uj ) |νj | |νj |

Z =

Pd (γ)f (zs |γ)Dk|k−1 (γ)dγ

(34)

Fig. 2. Evolution of PHD surface (represented by particles, as from [?], [?], [?], [?], [?]) in 3 consecutive scans: there are two targets heading east, and there are, respectively, 2, 2 and 1 false alarms. In this scenario, the PHD’s  estimate for the number of targets is close to 2. The integration of PHD in each of the circled regions (see leftmost plot) would be close to unity, indicating Pd (γ)f (zs |γ)Dk|k−1 (γ)dγ the tracker’s opinion that one target is present inside the corresponding region.

Hence, the continuous form of the bin-occupancy filter update equation is given by  Pm RPd (xi )f (zs |xi ) Dk|k (xi )=Dk|k−1 (xi ) (1−Pd (xi ))+ s=1 c(zs )µF A +

(35) C. The Approximation The bin-occupancy filter has a hidden approximation, and it involves re-inserting the result of (??) to (??) for the next frame of data. The fact is that there is no reason that Dk|k (x) from (??) should be the intensity of a Poisson process: that is, reverting to the bin-occupancy probabilities, although {Ui } may have been independent a-priori (before the frame of data), given Zk they no longer are. In particular, although (??) correctly updates the occupancy probability of each (even infinitesimal) bin pk|k (Ui ), it is not true that pk|k (Ui , Uj ) = pk|k (Ui )pk|k (Uj ). As such, in this section we stress the necessary approximation that after each prediction step the bins be independent, conditioned on the available observations: even the most basic prediction equation (??) would not work without this. We note that all practical target tracking systems that we aware of make some approximations or simplifying assumptions: The probabilistic data association filter (PDAF) [?] makes the posterior Gaussian; multi-hypothesis trackers (MHTs) [?] assume their kept hypothesis set is sufficiently rich to represent all hypotheses; the probabilistic multi-hypothesis tracker (PMHT) [?], while it does not approximate, makes a simplifying association-model assumption. Even Mahler’s PHD [?] is the first-moment of the joint multi-target density that would, ideally, be of greatest interest. D. Bin Occupancy Filter vs. Probability Hypothesis Filter It is immediately recognizable that the continuous form of the bin-occupancy filter prediction and update equations, (??) and (??), are the same as the probability hypothesis (PHD) filter [?] equations. The PHD filter is the first factorial moment density m[1] (x) ([?] p.147, or [?] p.1163) and the physical interpretation of the first factorial moment density can be found in [?] p.133 (see also [?]),

In figure ??, the PHD surface is depicted for a hypothetical scenario in which two targets (one at (400,200) and the other at (400,1500)) are traveling east with constant speed. The PHD surface peaks around the estimates of each target’s location. Moreover, the integral of the PHD in the proximity of each peak (i.e., the circled regions in the leftmost plot) would yield a number around unity. This would indicate that there exists in the corresponding region a target whose location – and the uncertainty of which – can be inferred from the shape of the PHD surface. In the leftmost plot, there are two false detections around coordinates (1500, 1500) and (1800, 1500). As seen in the plot, the PHD surface has non-zero regions with peak values considerably lower than those of the targets. These regions imply the possibility that those detections are due to new-born targets, and presumably the masses near there would either increase if these were corroborated later by further detections or die away if not. III. T HE C ARDINALIZED B IN O CCUPANCY F ILTER A. The Missed Detection Problem: Why Cardinalized? Since the cardinalized bin occupancy filter is directly motivated by – and will come to the same algorithm as – Mahler’s development of the CPHD from the PHD, we shall begin with some discussion of why Mahler created the CPHD: it is considerably more complicated than the original PHD. The concern is that the PHD handles missed detections in a way that is unfamiliar and somewhat unsatisfying to most target tracking practitioners. Recall that if there is no measurement at the current scan (|Zk | = 0), then Dk|k (x) = [1 − Pd (x)] Dk|k−1 (x)

(37)

infinitesimal region (x, x + dx), divided by dx

is the PHD update equation (??). Now, consider a Markov chain model for a generic “target existence” process, Yk . It can easily be shown [?], [?] that when there is no detection at time k, one obtains (1 − Pd )Q (38) P{Yk|k = 1 |Zk | = 0} = 1 − Pd Q

Since the bin-model underlies the surface defined by the first factorial moment of a finite point process, it is perhaps natural that we arrive the same filtering equations as the PHD filter. The approximation discussed in the previous section presumably applies to the PHD as well.

where Q = P{Yk = 1|Yk−1 = 1} and the probability of detection is Pd (naturally, given existence). Consider the simplistic case where there is at most one target and the target existence probability, the integrated PHD over the surveillance region, is approximately unity at time k − 1. If a missed

m[1] (x)

= Probability that one target is located in the(36)

∆

6

detection occurs at k, the PHD might according to intuition be expected to have equation (??) as its update. But instead it follows equation (??). This was pointed out in [?], and to address it Mahler proposed a modified version of the PHD filter, the cardinalized PHD (CPHD) [?], [?]. The CPHD filter relaxes the Poisson assumption on the number of targets and, instead, propagates the full probability mass function of the number of targets. The CPHD (which we are motivating here but have not begun to explain) does indeed solve the PHD’s missed detection problem, as we next illustrate by an example. We must note that the CPHD is not perfect, and there is a tendency [?] for a missed detection in one location to provoke apparent increase in the number of targets elsewhere; but we do not treat that concern here. Scan no 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Detection Yes Yes Yes Yes Yes Yes Yes Yes Yes No No No No No Yes No No No No No Yes

ˆP HD N 1.095 1.1026 1.1032 1.1033 1.1033 1.1033 1.1033 1.1033 0.0883 0.00706 0.000565 4.5e-005 3.6e-006 1.015 0.0812 0.00650 0.000520 4.2e-005 3.3e-006 1.015

ˆCP HD N 1.255 1.2509 1.247 1.2468 1.2468 1.2468 1.2468 1.2468 0.215 0.0193 0.00174 0.000159 1.5e-005 1.017 0.213 0.0196 0.00177 0.000161 1.5e-005 1.017

TABLE I S INGLE TARGET - E STIMATED NUMBER OF TARGETS FOR EACH SCAN .

Consider a simulation in which a target is detected at scans 1-9, 15 and 20-30. There are no measurements observed during scans 10-14, and 16-19. The parameters are Pd = 0.9, Q = 0.99 and birth density = 0.15 [/surveillance region] (birth density is a Poisson process with parameter 0.15, appropriately normalized). As seen in Table ??, the estimated number of ˆP HD , is confidently close to targets for the PHD filter, N 1 (which is the truth) until scan 10. With the first missed detection at scan 11, it drops down to 0.08, i.e. (1 − Pd )Q. As long as there is no detection, it decreases by the factor ˆP HD jumps 0.08 at each scan. At scan 15, with a detection, N back to 1 and at scan 16, it drops back to 0.08; and so on. ˆCP HD , The estimated number of targets for the CPHD filter, N appears to obey equation (??). The Poisson assumption [?] on the number of targetoriginated measurements implies that the bins are occupied independently. The Poisson assumption is apparent in the original derivation of PHD as well ([?] page 1166). This is the reason behind the target-death problem. Most trackers estimate a target’s state as an aggregated region in the state space. That is, a missed detection in a conventional “lumped” target tracker

Fig. 3. The underlying hidden Markov model of the inference process of PHD (left) and CPHD (right) filters. Observe that for the PHD the target number is derived from (and completely dependent on) the PHD itself, whereas for the CPHD it is exogenous and affects the PHD surface.

is not problematic, since the probability that a target remains is given by (??) with relatively large Q. On the contrary, the PHD’s bins each have infinitesimal Q, and since the bins are not linked in such a way that, speaking loosely, would mean that the lack of a target in a given bin does not significantly alter the fact that there is a target in at least one of the surrounding bins, all PHD bins are reduced abruptly according to (??) with Q ≈ 0 – this amounts to (??). To reiterate: the culprit is the implied assumption of independence between bins. The CPHD filter overcomes the target-death problem. The difference between the two inference processes is depicted in figure ??. In the CPHD filter, the target number process is above that of the individual bin occupancies in the hierarchy; in the PHD filter the number of targets is inferred from the PHD surface, i.e. from the bin probabilities. The target number process imposes the kind of dependence on the bin probabilities that the PHD lacks. B. The Cardinalized Bin Occupancy Filter: Statistical Model and the Notation Consider the right panel of figure ??. We want to update the probability mass function (pmf) → − → − p( U k , Nk |Z1k ) = p( U k , Nk |Zk , Z1k−1 ) (39) → − = p( U k |Nk , Zk , Z1k−1 )p(Nk |Zk , Z1k−1 ) S → − where U k = i Uk (i) represents the vector of bin-occupancy → − events at time k. Hence, p( U k |Z1k ) is the (discrete) surface defined by the PHD. As is common, we assume the previous → − estimates of p( U k−1 ) and p(Nk−1 ) are sufficient, → − → − → − p( U k |Nk , Z1k ) ≈ p( U k |Nk , Zk , p( U k−1 |Z1k−1 )) (40) p(Nk |Z1k ) ≈

p(Nk |Zk , p(Nk−1 |Z1k−1 ))

(41)

We need both, and clearly (??) must precede (??). C. The Cardinalized Bin Occupancy Filter: Prediction Equations We adopt the additional notations: c : number of clutter detections (false alarms) at time k c(zj ) : clutter spatial density ψ : number of new-born targets at time k n : number of targets at time k, i.e., Nk = n

7

the received observation set Zk = {z1 , z2 , . . . , zm } via Bayes’ formula, namely,

Throughout, we use, as before, pk|k−1 (Ui ) pk|k−1 (n) pk|k−1 (c)

= P{Uk (i) = 1|Z1k−1 } = P{there exists a target in bin i}

(42)

= P{Nk = n|Z1k−1 } = P{there are n targets at time k}

pk|k (Ui )

=

(43)

pk|k (n)

=

= P{there are c false alarms at time k}(44)

The prediction equation for the bin probabilities is the same as in the PHD case (see (??)). P{Uk (i) = 1|Z1k−1 } = X b(i) + Ps (xj )p(i|xj )P{Uk−1 (j) = 1|Z1k−1 } (45) j

k−1 P{Uk (i)=1|Z } b(i) 1 = + |νi | |νi |

p(i|xj )

P j

|νi |

Ps (xj )

P{Uk−1 (j)=1|Z |νj |

k−1 } 1 |νj |

f (Zk ,Ui |Z

= (46)

=

P P

k−1 )= 1

P{Uk (i) = 1|Z1k−1 } b(i) lim = lim (47) |νi | |νi |→0,∀i |νi |→0,∀i |νi | X p(i|xj ) P{Uk−1 (j) = 1|Z1k−1 } + lim Ps (xj ) |νj | |νi | |νj | |νi |→0,∀i j Z Dk|k−1 (xi ) = b(xi ) +

f (xi |γ)Ps (γ)Dk−1|k−1 (γ)dγ (48)

The second term in (??) constitutes a Riemann sum [?] and since the term inside the sum is Riemann-integrable, convergence to integral is ensured as bin-volume, |νj |, approaches to zero. Equation (??) is the prediction step. We define the average probability that a target survives as10 P p (Uj ) Pss = j Ps (xj ) P k−1|k−1 . Then, the prediction step p (U ) l

k−1|k−1

l

for the pmf of the number of targets, pk|k−1 (n), is given by Pn pk|k−1 (n)= p(ψ=n−i) (49) i=0     0 P∞ n 0 i pk−1|k−1 (n0 )(1−Pss )n −i Pss n0 =i n0 − i Assuming that the pmf of the birth cardinality, p(ψ), is known, the above equation is derived by applying of the total probability theorem to each of the following hypotheses: given that n − i new targets are born at time k, and given that there were n0 targets at time k − 1, each hypothesis identifies the n0 − i dead-targets, and all of these hypotheses are assumed to be equally likely to occur. D. The Cardinalized Bin Occupancy Filter: Update equations The modified bin-occupancy filter propagates the bin probabilities, pk−1|k−1 (Ui ), and the probability mass function (pmf) of the number of targets, pk−1|k−1 (n), and updates them by 10 limiting case where all the bins’ volumes go to zero, Pss = R In the D (α) dα Ps (α) R k−1|k−1 Dk−1|k−1 (γ)dγ

(51)

c

n

f (Zk ,Ui ,n,c|Z

(52)

k−1 ) 1

P P

k−1 k−1 f (Zk |Ui ,n,c,Z )p(Ui ,n,c|Z ) 1 1 n

Pc P

Pc Pn c

n

k−1 k−1 f (Zk |Ui ,n,c,Z )p(Ui |n,c,Z )pk|k−1 (n)pk|k−1 (c) 1 1 k−1 f (Zk |Ui ,n,c,Z ) 1

P j

Using the definitions (??)-(??) and taking the limit where all the bin volumes go to zero,

(50)

Next, we express the quantities in the right hand side(s) of the above equations (??) and (??) in more explicit forms. We assume that the pmf of number of clutter (false alarm) observations p(c) is given.

=

We rewrite (??), by dividing both sides with |νi |,

f (Zk , Ui |Z1k−1 ) f (Zk |Z1k−1 ) f (Zk |n, Z1k−1 ) pk|k−1 (n) f (Zk |Z1k−1 )

n pk|k−1 (Uj )

pk|k−1 (Ui )pk|k−1 (n)p(c)

The last equation follows from the following assumptions: 1) The number of false alarms, c, and the number of targets, n, are independent. 2) The number of false alarms, c, is independent of the measurement sequence, Z1k−1 . 3) The location of the targets are independent given the number of targets. 4) A target can generate at most one measurement. 5) Each target evolves and generates measurements independently. 6) False alarms are independent from target-originated measurements. The third assumption above amounts to pk|k−1 (Ui |n, c)

= =

pk|k−1 (Ui |n) n P pk|k−1 (Ui ) (53) j pk|k−1 (Uj )

and follows from the fact that given the PHD surface (the vector of the bin probabilities), and given that there is only one target, the probability that bin i has the target is pk|k−1 (Ui ) P j pk|k−1 (Uj )

(54)

If there are n such targets, probability that bin i has one of these targets is given by 1−P{bin i has none of these targets}. Then, for sufficiently small bins11 , we have  n pk|k−1 (Ui ) pk|k−1 (Ui | n) = 1 − 1 − P (55) j pk|k−1 (Uj ) n pk|k−1 (Ui ) (56) ≈ P p j k|k−1 (Uj ) assuming the bin volume goes to zero, |νi | → 0. Note that this assumption imposes a constraint on the integration of the 11 Note that the approximation in (??) becomes exact as the bins’ volumes go to zero, i.e. pk|k−1 (Ui ) → 0

8

bin-occupancy probabilities, pk|k−1 (Ui )

∞ X

pk|k−1 (Ui |n)pk|k−1 (n)

(57)

=

pk|k−1 (Ui ) npk|k−1 (n) P j pk|k−1 (Uj ) n=0

(58)

=

P

E[n] pk|k−1 (Ui ) p j k|k−1 (Uj )

(59)

=

n=0 ∞ X

Hence, the expected number of targets, E[n] = P p (U ), is can be found by integrating the binj k|k−1 j occupancy probabilities. The modified bin-occupancy filter uses likelihood ratios. In other words, each hypothesis12 , f (Z|Ui , n, c) or f (Z|n, c), is compared to the null hypothesis that all measurements are false alarms: f (Z|all clutter). The likelihood ratio is formed as L(Z|n, c)

f (Z|n, c) f (Z|all clutter)

=

(60)

Fig. 4. Evolution of CPHD surface and the number of targets pdf in 3 consecutive scans: there are two targets heading east, and there are, respectively, 2, 2 and 1 false alarms. As seen in the lower plots, the numberof-targets posterior pdf has a significant mass at n = 2, indicating there are two targets (p(n = 2) = 0.987 in the lower left-most plot). The integration of the CPHD surface over whole surveillance region is close to 2 as well, which is equal to the mean of the target number posterior pdf, pk|k (n). Similar to the PHD filter, the integration of CPHD surfaces in each of the circled regions (see upper-leftmost plot) is close to unity.

Then the bin probabilities p(Ui ) and the cardinality pmf p(n) are updated by p(Ui |Z)

=

L(Z|Ui ) p(Ui ) L(Z)

=

(1−Pd )·

L(Z|Ui ,Vi =1) L(Z|Ui ,Vi =0) p(Ui )+Pd · p(Ui ) L(Z) L(Z)

P|Z| h P∞ c=0

=

n=0

L(Z|Ui ,n,c,Vi =0)

P j

(1−Pd )·

P|Z| h P∞ c=0

P|Z| h P∞ c=0

n=0

n=0

c=0

n

p(n) p(c)

p(Uj ) p(Ui )

i

L(Z|Ui ,n,c,Vi =1)

P|Z| h P∞

i

P

i

n

p(n) p(c)

p(Uj )

(62)

p(Ui )

i

n=0

L(Z|n,c)p(n) p(c)

and p(n|Z)

=

L(Z|n) p(n) L(Z)

P|Z| =

c=max{|Z|−n,0} P|Z| h P∞ c=0

n=0

L(Z|n,c)p(c)

i

p(n)

(63)

L(Z|n,c)p(n) p(c)

In the appendix ??, we explicitly derive the terms in the numerators and denominators of the above equations. The resulting forms of (??) and (??) are given in (??) and (??). It can be noticed that these modified continuous bin-occupancy filter equations are identical to the original CPHD filter [?] (predictor and corrector) equations. In figure ??, the exact same scenario used in figure ?? is considered for the CPHD filter. The surfaces of binprobabilities for 3 consecutive scans are given in the upper plots, and the corresponding target number probability mass functions (pmf) are given in the lower plots. As seen in the figure, the target number pmfs have only one significant mass, and that is at n = 2. This is the correct number of targets as the scenario described earlier in section ??. Also, note that this pmf is truncated at n = 10, assuming there will be at most 10 targets. This is what needs to be done in implementation to avoid the infinite summation from, for example, equation 12 In

IV. S UMMARY

L(Z|n,c)p(n) p(c)

j

+Pd ·

(61)

(??). In other words, in practical situations it is reasonable to assume that there will be at maximum Nmax targets in the surveillance region. This number can be as large as one needs to have, we chose it very low for illustration purposes.

the sequel, we drop the time subscripts for simplicity.

Mahler’s PHD is an attempt to simplify the calculation of family of joint probability densities by computing only the first moment, and has the useful characteristic that its integral over a volume of space is the (posterior) expected number of targets in that volume. The PHD filter is a large step forward, but we believe its development, being via setderivative and mathematical concepts, is perhaps unfamiliar and uncomfortable to many engineers. Thus the first aim of the paper: we have provided a “bin” model of probability (motivated by the physical interpretation of first factorial moment density) and showed that it gives rise (in the limit of infinitesimal bins) to exactly the PHD prediction and update equations. In other words, our development is not a stepby-step interpretation of Mahler’s original derivation, but an alternative route that shares the start- and the end-points of the original. A weakness of the (original) PHD is that a missed detection gives rise to an unreasonably-speedy target demise, as measured in PHD terms. The reason behind this is likewise exposed via the physical-space (bin) model we have contributed: the problem is an assumed independence between occupancy events of the bins. Mahler recognized this problem, and his “fix” – the CPHD – is both effective and clever. Despite its excellent performance the CPHD is, unfortunately, rather abstruse. Thus our second aim has been to extend the physical “bin-occupancy” model such that the CPHD equations arise: it turns out that we need an additional inference level by which the number of targets is estimated separately from their locations.

9

The underlying Hidden Markov models are presented (figure ??), and they reveal important intuitive aspects of the PHD filters: the PHD filter’s estimate of the number of targets relies only on the bin-occupancy probabilities, hence the independence (i.e. Poisson) assumption results in noisy estimates on number of targets. On the other hand, the CPHD filter imposes a constraint on the bin-occupancy probabilities such that the average number of targets, which is estimated in parallel by the target number process, must be equal to the integral of the bin-occupancy surface. Some consequences of this constraint are discussed in [?] as the “spooky” behavior of the CPHD filter. The PHD and CPHD offer a tractable and elegant means to estimate both target number and location. There has been considerable recent interest in applying the PHD to a variety of tracking systems. We must credit the work in [?] and [?] for considerably simplifying implementation: a particle approach is no longer necessary in the case of linear Gaussian systems that can be represented by Gaussian mixtures. We hope that we have been able to contribute to (C)PHD understanding, and acceptance as well. A PPENDIX

A more mathematical treatment of the bin model is given here, and is particularly relevant to (??) and (??). As defined earlier, the set of points bin i contains are defined as νi . The (bounded) state space S ∈ Rd , where d is the dimension of the state vector, is partitioned such that νi = S

(64)

i=1

T

and νi νj = ∅, ∀i 6= j for i, j = 1, 2, . . . , n. We assume n is large enough to ensure “sufficiently small” partitions. The volume of the surveillance region is |V |, and since we assume |V | . Also, for this each bin has the same volume, |νi |, n = |ν i| partition, the generalized function below constitutes a Riemann sum, hn

= =

n X j=1 n X

P{

An event belongs to T ype − I |m events occurred} Pm =

n1 =0

P{The event is T ype − I|N1 =n1 ,N2 =m−n1 }

(67)

×P{N1 =n1 ,N2 =m−n2 } =

Pm n1 =0

P{The event is T ype − I|N1 =n1 ,N2 =m−n1 } ×

P{N1 =n1 }P{N2 =m−n2 } P{N1 +N2 =m}

(68)

n m−n1 γ γ 1 e−γ1 1 e−γ2 2 m n1 ! n1 (m−n1 )! = · (γ1 +γ2 )m n1 =0 m e−(γ1 +γ2 ) m! m n m−n1 γ γ2 1· m! 1 = 1 n1 · m γ1 +γ2 γ1 +γ2 n1 !(m−n1 )! n1 =0 γ 1 = 1 ·m· m γ1 +γ2

P

P

(

)

(

)

γ1 = γ1 +γ2

A. The Bin Model

n [

B. Conditional Probability of Random Events Here we give details on the conditional probability that is used to go from (??) to (??). We are interested in two types of events, {T ype−I, T ype− II}. N1 and N2 are independent random variables such that realization N1 = n1 implies there are n1 occurrences of T ype−I events. Assume that total number of events are known for a given realization i.e., n1 +n2 = m. Moreover N1 and N2 are Poisson distributed with parameters γ1 and γ2 , respectively. We want to calculate the probability that a randomly picked event (among these m events) belongs to T ype − I. Then,

ξ(xj )|νj |

(65)

f (xi |xj )Ps (xj )D(xj )|νj |

(66)

j=1

where13 xi ∈ νi . The function ξ(·) is continuous, real-valued and bounded by [0, 1] since each term that forms ξ(·) defines a probability. Therefore, ξ(·) is Riemann-integrable [?]. A finer partitioning of S is such that given n0 > n, a new set of S n0 n0 {xi , νi }i=1 are chosen where i=1 νi = S. Hence as n → ∞ the Riemann sum converges to an integral – for example, (??) converges to (??). 13 x is an arbitrary point inside the bin i. It is not necessary, but for intuition i we choose it to be the middle point.

(69) (70) (71) (72)

Equation (??), evidently and unsurprisingly, is the expectation 1 divided of a binomial random variable with probability γ1γ+γ 2 by m. Hence (??) holds. The transition from equation (??) to (??), with γ1 = µF A , P and γ2 = j Pd (xj )pk|k−1 (Uj ), uses P{zs is false alarm|Z1k−1 } µF A P = µF A + j pk|k−1 (Uj )Pd (xj )

(73)

Similarly, the probability that zs is originated from the target in bin i is given by P{zs from target in bin i|Z1k−1 } pk|k−1 (Ui )Pd (xi ) P = µF A + j pk|k−1 (Uj )Pd (xj )

(74)

C. Details of CPHD derivation We derive the joint measurement likelihoods under any given hypothesis, f (Z|Ui , n, c) and f (Z|n, c) and then form the corresponding likelihood ratios to achieve the explicit forms of update equations (??) and (??). For the sake of clarity we make a simplifying (and really quite unnecessary) assumption that Pd is constant over the surveillance region. 1) Likelihood ratios for general case: m − j false alarms, n targets: To arrive the general form of the modified binoccupancy filter update equations, we need to derive the terms in the numerator and the denominator of (??). We assume there are m − j false alarms, where |Z| = m is the number of measurements received at current time. For notational simplicity, we define the likelihood of a target originated measurement as X p(Uj ) ∆ (75) f (z|xj ) P Lz = l p(Ul ) j

10

First we express the term f (Z|n, c) f (Z|n,c=m−j)

=

  n

1



c(z1 )...c(zm−j )Lz



m

m−j+1

...Lzm +···

j

m−j ···+c(zm−j+1 )...c(zm )Lz1 ...Lz



(76)

j P (1−Pd )n−j d

j

Since it is known  that there are m − j false alarms we m consider all combinations to be equally likely, and m−j

the remaining j measurements from some of the n  originate  n targets. This can happen in different ways, and means j

that the remaining (n − j) targets are not detected. In the case that bin i has a target but it is not detected (i.e., Vi = 0), the joint likelihood of the measurements is the same as above with the knowledge that (n − 1) more targets are there and all detections are missed. (77)

f (Z|Ui ,n,c=m−j,Vi =0)=f (Z|n−1,c=m−j)



=

n−1

1



m

h

P

∞

L(Z|Ui ,Vi =0)=p(c=m)

c(z1 )...c(zm−j )Lz



Fig. 5. An illustration for the update of the bin probabilities p(Ui ) with a measurement

m−j+1

n=1

...Lzm +

(1−Pd )n−1

l

j



Pn

p(n)

p(Ul )

P

∞

1 ·(1−P )n−2 (n−1)p(n) d n=2 m

+p(c=m−1)

m−j

l

i

h

j P (1−Pd )n−1−j d

···+c(zm−j+1 )...c(zm )Lz ...Lz 1 j

Pn

×Pd

Lz2 Lz3 Lz1 Lzm + + +···+ c(z1 ) c(z2 ) c(z3 ) c(zm )



p(Ul )

i

P

Given bin i has a target and that it is detected, then m possibilities exist: any one of z1 , . . . zm is from the target in bin i. In each case, there are m − j false alarms and the remaining j − 1 measurements are from some of the remaining (n − 1) targets. Hence the joint measurement likelihood function is given by

∞

1 (1−Pd )n−3 (n−1)(n−2)p(n) n=3 m(m−1)

+p(c=m−2)

×P 2 d

h

Pn p(Ul ) l i

Lz L Lz Lz Lz Lz Lz Lz m−1 zm 1 2 3 + 3 +···+ 1 2 + c(z1 )c(z2 ) c(z1 )c(z3 ) c(z2 )c(z3 ) c(zm−1 )c(zm )

.. . P

∞

1 p(n)(n−1)···(n−m)(1−P )n−m−1 d n=m+1 m!

+p(c=0)

Pn l

×P m d

f (Z|Ui ,n,c=mj,Vi =1)

=

 1 f (z |x ) 1 i m

1



m−1



n−1

Lz Lz Lz ···Lzm 1 2 3 c(z1 )c(z2 )c(z3 )···c(zm )

Pm

L(Z|Ui ,Vi =0)=

h i j−1

∞

×

n=j+1

1 m−1



n−1

p(c=m−j)



n(n−1)···(n−j)p(n)(1−Pd )(n−1)−j

m−j

P

j−1 (1−Pd )(n−1)−(j−1) d

(78)



p(Ul )

(80)

0