Tracking Many Objects with Many Sensors

1 Introduction

The problem of tracking multiple objects over time has long been studied in the literature on data asso­ ciation [Bar-Shalom and Fortmann, 1988; Bar-Shalom, 1992]. The problem is defined as that of associating a set of current observations with a set of existing "tracks" or object trajectories, creating new tracks as needed. Radar tracking of multiple aircraft is the canonical application. In AI, the problem of object identification is essentially Statistics Dept.

*Permanent address: Statistics Dept. , Hebrew University of Jerusalem.

University of California, Berkeley { ostland ,ri tov} @stat. ber keley.ed u

the same: deciding if some newly observed object is the same as some previously observed object. Solving this problem is essential for any intelligent agent that rea­ sons about individual objects. Huang and Russell [1997; 1998] provide a fairly general formulation of the prob­ lem and describe an application to traffic surveillance. Other possible applications range from removing "dupli­ cate" entries from databases to re-recognizing locations during exploratory map-building.

The object identification problem is difficult because sensors are noisy, objects look similar, and object be­ haviors are unpredictable. · This leads to a large number of possible assignments specifying identities among ob­ served objects. For example, in the traffic surveillance application studied by Huang and Russell, the sensors are cameras at various locations on a freeway network and the objects are vehicles. Many thousands of ve­ hicles pass each camera, and the system must decide whether each vehicle is a new vehicle or the same as one previously observed at a different location. Over time, these decisions give rise to hypothetical vehicle trajectories. Deriving a set of trajectories is the first step of many traffic surveillance applications, such as the average link travel time between locations or origin­ destination counts along different routes through the sys­ tem. Moreover, sudden changes in these quantities can be indicators of highway incidents, such as accidents or breakdowns.

Adopting the notation used by Huang and Russell, let w E r2 denote an assignment placing pairs (or, more gen­ erally, sets) of observed objects into equivalence classes, where each class represents an existing object, and let 0 denote all observations made to date . Then the poste­ rior probability that two objects a and b are the same is given by

In evaluating terms such as P(wiO), we will make use of the fact that the prior P(w) can be assumed uniform.

This is because the probability of an assignment w, in the absence of observations linking the objects, must be invariant under renaming of the objects. This it the exchangeability assumption of Huang and Russell.

Other quantities can also be calculated by summing over w. For example, in freeway surveillance, w speci­ fies the correspondence between vehicles observed at up­ stream and downstream sensor locations . For a given w, the average link travel time LTT( w) between the two locations can be calculated directly if the observations include the arrival time at each location. Then the pos­ terior expectation of the link travel time is

This paper addresses the two principal difficulties that arise in putting such equations into practice.

Section 2 deals with the computation of the P(wiO) terms-in particular their decomposition into tractable local models that can be estimated from data. We show that the decomposition proposed by Huang and Rus­ sell using appearance probabilities, while adequate for the case of two sensor locations, does not scale up to handle the decomposition of a global model for many sensors. In fact, this appears to require the estimation of intrinsic parameters of the observed objects, which render succes­ sive observations conditionally independent.

Section 3 deals with the intractability of the summa­ tion in Eq. (1), which includes an exponential number of terms. Whereas Huang and Russell describe a heuris­ tic scheme that seems to work well in practice, we ap­ ply the Markov chain Monte Carlo (MCMC) method, a general-purpose approximation algorithm for probabilis­ tic inference that can be shown to converge in polyno­ mial time for the specific inference problem involved in object identification. Furthermore, the algorithm can be adapted easily to incorporate online updating of the probability models required for computing P(wiO) . The overall scheme is in fact an online EM algorithm with MCMC as an approximate E-step. 1

Section 4 describes an applic;ation of the new approach to data extracted from a freeway simulation. We show that the estimation of intrinsic parameters, as described in Section 2, successfully handles some multi-camera sce­ narios for which the appearance probability models of Huang and Russell are not applicable. We also show that the MCMC method allows accurate estimation of long-range origin/destination information even when the individual links in the sensor chain are highly unreliable.

2 Scaling up to multiple sensors

As mentioned in the Section 1, calculation of assignment probabilities P(wiO) is crucial for object identification . The calculation will be done using probability models that, in some way, capture the properties of the sensors and the behavior of the objects being tracked. We hope to find models that allow decomposition of the global assignment probability in much the same way that local causal models allow decomposition of joint probabilities in Bayesian networks.

1 The use of MCMC with EM is well-known in statis­ tics [Wei and Tanner, 1990]; to our knowledge, its use for online learning and for object identification is novel.

We begin by describing the models used by Huang and Russell, explaining how they work only for two sensor lo­ cations. We then describe an alternative approach based on estimation of intrinsic properties of objects.

2.1 Problems with appearance probabilities

We begin with the case of two consecutive cameras A and B as considered by Huang and Russell . Let the observations at each be OA = {of, ... 0~} and 0 8 = {of, ... o~}, and let WAB be an assignment pairing up objects observed at A with objects observed at B. 2 Then

where a is , again, a normalizing constant. The term p ( w AB I oA) is dropped by exchangeability­ conditioning on the initial observations provides no in­ formation about matching with the subsequent objects. The approximate equality in the last line arises from the assumption of approximate independence among vehicle trajectories. 3

In Eq. (2), the terms P(o~lo:, a= b) are called appear­ ance probabilities since they describe "how an object can be expected to appear at subsequent observations given its current appearance" [Huang and Russell , 1997]. The appearance probability models for freeway vehicles are composed of factors such as the arrival time at B given the arrival time at A, the measured colour at B given the measured colour at A, and so on. Huang and Rus­ sell show how these models can be estimated online from matched vehicles in a very st raightforward way, avoiding the need for camera calibration.

2 Since in general m # n and conservation of objects is not assumed, WAB may include unpaired objects from either or both sensors.

3 Trajectory independence is a reasonable assumption only if the individual models are conditioned on some global con­ text variables such as the current link travel time.

Let us now extend this approach to three cameras, using the scenario of (Figure 1) as an example. An as­ signment WABC now specifies sequences of three obser­ vations that belong to a single object, and can be de­ composed into two pairwise assignments WAB and wac. As in Eq . (2), we can apply Bayes' rule and eliminate P(wAsciOA); then we can apply the chain rule:

where WABC is replaced in the last line by WAB because assignments of vehicles at C carry no information about A and B. Now the last term on the RHS of Eq. (3) can be written as the product of appearance probabilities between A and B , as in Eq. (2). However , the first term cannot be simplified to give the appearance probabilities between B and C; that is ,

To see why, consider the ext reme case in which cam­ eras A and C can read the license plate of each vehicle, but camera B is broken. Then the posterior distribu­ tion for assignments should have all its mass on assign­ ments WABC that correct ly match up vehicles at A and C; whereas both the pairwise models will be uninforma­ tive and hence will fail to propaga,te information from A to C.

Two possible fixes are 1) use multicamera. models, e.g., P(Oc lOA , 0 8 , WABC ), and 2) estimate models for all camera pairs , e.g. , P(OciOA ,wAc). The first fix enlarges the model dimension and scales exponentially with the number of cameras. The second fix requires a quadratic number of models; moreover, it is unclear how to combine the predictions of these models. In summary, it seems that , despite their many advantages, appearance probabi li ty models apply only to the two-camera case.

2.2 Spatial decomposition via intrinsic properties

The example of the broken camera. at B raises the fol­ lowing problem: we wish to propagate information be­ tween nonadjacent sensors, yet we do not want to have to employ nonlocal probability models , since such models result in a. combinatorial explosion in the number of pa­ rameters to be estimated. The solution is, essentially, to let the objects themselves carry the necessary informa­ tion. As with Kalman filters and hidden Markov models, hidden state variables can render current observations conditionally independent of previous observations . This provides a decomposition of the global model.

Let H 1 represent the hidden state of the objects ob­ served at location l , and let h 1 range over possible values of H 1 . Notice that. , given the hidden state, the observa­ tions at a camera are independent of all the other ob-. k . . . servations. More formally, ?(0 1 10 ,HJ) = P(0 1 IHJ)

for all j and k. Now, we can introduce hidden state into Eq. (3) by summing over hA, h 8 , and he and simpli­ fying using conditional independence, to yield a nested sum exactly as in the derivation of the forward equations for HMMs. The only significant difference is that the re­ lationship between successive hidden state variables is only meaningful if we condition on the assignment so we know which vehicle is which:

Thus , the introduction of hidden state solves the problem of Section 2.1. Unlike Eq. (3) , Eq. (4) calls only for two-camera models. Moreover , these models need be estimated only for neighbouring camera pairs.

As with appearance probabilities , we can decompose the expressions in Eq. ( 4) into models for individual ve­ hicles by making the appropriate independence assump­ tions (again, assuming the models depend on some global context variables). The models we obtain are the transi­ tion models such as P ( h~ I h: , a = b) and the sensor mod­ els such as P(o:lh:). As with HMMs, these models can be learned online using EM [Dempster et al., 1977), as we show below. The process is complicated by the fact that we must simultaneously estimate both the hidden state variables of the observed objects and the global assignment saying which object is which.

For traffic surveillance and many other applications , some aspects of the hidden state do not change over time.

We call these intrinsic variables; for traffic surveillance, these include colour, length, width, and so on. Intrin­ sic variables have no transition model but often have very noisy sensor models, specific to each sensor loca­ tion. Dynamic variables such as lane, speed, and arrival t ime must be tracked as the vehicle progresses through the freeway network. Often they have relatively noise­ less sensor models. Thus, the hidden variables h can be divided into the intrinsic variables, t, and the dynamic variables, J. Similarly, the observed variables 0 can be divided into i , observations oft, and d, observations of c5. Figure 2 illustrates all the independence assumptions we have made. It is relatively easy to augment this model to include dependencies between , for example, vehicle size and lane.

3 Approximate inference and online model updating

As mentioned in Section 1, the expressions for the proba­ bility of identity (Eq. 1) and other quantities involve ex­ ponentially many terms. It can be shown that the infer­ ence problem is equivalent to computing the permanent of a matrix and hence is #P-complete. It is possible to compute the most probable assignment for n vehicles in O(n 3 ) time using the Hungarian algorithm [Cox and Hin­ gorani, 1994]; as Huang and Russell point out, however, this assignment may be of little interest if the individual matches therein are highly unreliable. They developed a heuristic "leave-one-out" algorithm with runtime 0( n 4 ) that alleviates this problem and works well in practice. The algorithm identifies individual matches that exceed a reliability threshold and then treats those matches as if true. The matches are used to compute link travel time and to update the appearance probability models. A major drawback of this approach is that the if the fraction of reliably matched vehicles on each link is sig­ nificantly below 100%, as often happens, the number of vehicles that can be tracked across a multi-link freeway network is vanishingly small [Huang and Russell, 1998].

In this section, we describe an alternative approach to the inference problem based on Markov chain Monte Carlo (MCMC). Roughly speaking, MCMC approxi­ mates sums of probabilities such as Eq. (1) by sampling a small number of high-probability terms; thus , for ob­ ject identification, it considers a small number of likely assignments and likely values for the hidden state vari­ ables. Given any particular assignment and set of val­ ues for the hidden variables, estimating the transition and sensor models is trivial because the data is com­ plete. This suggests an online EM scheme for learning the models, as shown in Figure 3.

3.1 Introduction to MCMC inference

MCMC inference is a general-purpose method for ap­ proximating E, (! ( x)), the expected value of the func­ tion f when its argument x is drawn from a probability distribution rr. Typically, rr(x) = P(xie), a posterior distribution over x given evidence e.

For each newly detected vehicle

Augment Markov chain state to include the vehi cle Repeat until models converge

E: Run MCMC, sampling from s:1 and H M: Update models from sampled values

Figure 3: The overall inference scheme: online EM using MCMC for the E-step and updating the models direct ly from the states sampled by MCMC.

Ideally, E, (f(x )) can be approximated simply by sam­ pling from rr ; if each sample can be drawn in con­ stant time, then from Chernoff bounds we will have a polynomial-time approximation method. For general rr , we know this cannot exist ; sampling from an arbitrary distribution in constant time is not always possible.

MCMC inference [Metropolis et al., 1953; Gilks et al., 1996] provides a general method to generate sam­ ples from rr by defining a Markov chain whose states are the objects x and whose stationary distribution is rr( x ). Samples are produced by simulating the Markov chain and selecting states from among those visited.

In the Metropolis-Hastings method , t ransitions in the Markov chain are constructed in two steps:

• Given the current state x, a candidate next state is sampled from the proposal distribution q( .r'ix), which may be (more or less) arbi trary.
• The transition to x' is not automat ic, but occurs with an acceptance probability defined by

Notice that to use this rule we need only be able to compute ratios rr( x')/ rr(x) , which conveni ently avoids the need to normalize rr.

Provided q is defined in such a way that the chain is er­ godic, this transition mechanism defines a Markov chain whose stationary distribution is rr( x) , and hence the av­ erage value off over the sampled states will converge to the desired value E,(f(x)).

The complexity of the. original inference problem is reflected in the mixing rate of the Markov chain , which determines the speed at which the sample average con­ verges. Jerrum and Sinclair [1997] have shown that a Markov chain defined on assignments , as described be­ low, yields a fully randomized approximation scheme for estimating expectations over the probability distribution on assignments. This means that, if we want to approxi­ mate some such function f to within ratio 1 + <in a world with n objects, the chain wil l run in time polynomi al in n and 1/< to approximate f with probability > 0.75. The probability may be boosted to > (1 - J) by run­ ning the chain O(log J-l) times and taking the medi an value. Of course, this assumes that proposals and accep­ tance probabilities can be computed in time constant in n , which we show below. Therefore, MCMC provides an efficient approximation scheme for object identiftcation.

3.2 Application to traffic surveillance

As we have seen in Section 1, traffic surveillance, as defined here, involves taking the expected value of a quantity f over all assignments w, given the obser­ vations 0. The calculation of this expected value, E P(wiO) (f(w)), can be approximated by sampling from an ergodic Markov chain with state space nand station­ ary distribution P(wiO). The quantities to be estimated, which include LTTs and origin destination counts, are very simple to derive for any given assignment w. An assignment specifies trajectories through the sensor net­ work with known times at each sensor location, enabling us to read off the desired quantity directly.

The transitions in the Markov chain may be set up in many ways, as long as ergodicity is ensured. In our approach, each transition is simply a swap between two assignment pairs across one pair of sensors. For example, in Figure 4 observed object w at sensor j is originally matched with observed object y at sensor j + 1, while object x at j is matched with z at j + 1. A transition leads to w being matched with z, and x with y. The

observation pairs to be switched are suggested in a very simple manner. The algorithm cycles through all pairs of adjacent sensors, and all matched pairs currently at each sensor. For each such pair, a plausible second pair to swap with is then chosen uniformly at random. Such a chain is provably ergodic.

Once a pair is chosen, it is accepted or rejected based on the acceptance probability as Eq. (3.1). This involves the computation of the ratio

Because of the simple form of the swapping proposals, the q values are trivial to compute. We can also derive a general expression for P(wiO), which will permit us to simplify the P(w'IO)/ P(wiO) ratio and make its calcu­ lation more efficient. Let us demonstrate this using the transition in Figure 4.

Generalizing Eq. ( 4) for the observation sequence 0 1 , 0 2 ... ok, we express P(w1. k IO) as a nested sum

As was noted in Section 2.2, the hidden variables can be divided into dynamic variables c5 , which change over time, and intrinsic variables t, which do not. The observed variables are analogously divided into d and i. Here , we assume that all the dynamics can be observed reliably. Thus, assuming that t and c5 are mutually independent , P(hJ+ 1 Ihj, Wj ,j +1) becomes P(tj+lltj, Wj,j+1)P(JJ+ 1 Ic5j ,wj,j+1), and P(OJ_IhJ) be­ comes P(ijltj)P(dJic5J). Now, we can make use of the assumption that the intrinsic variables remain constant across all sensors to discard the P ( tj +i I tj , w j ,j + 1) terms. Similarly, since we assume that dynamic variables are observed with perfect reliability, we can replace c5s with ds . The P(dJic5J) terms can then be dropped. So, given our assumptions, Eq. (6) simplifies as follows :

Finally, note that the only hidden variables we are sum­ ming over now are the intrinsics t. Moreover, t is the same across all cameras by assumption, so we can com­ bine all the summations. The P(dJidJ- 1 ,wj-1,j ) terms are independent of t, and so can be moved outside the summation

We can now use trajectory independence to factor this equation into separate terms, one for each trajectory. Let the variable t represent a trajectory within a specified w and let it represent the hidden intrinsic variables for the putative object that follows trajectory t. Define it and dt analogously.

Eq. (7) may now be factored into a product of proba­ bilities for individual trajectories:

where

We can now evaluate the ratio P(w'IO)/ P(wiO) in a sim­ plified form. As shown in Figure 4, each Markov chain transition affects only two trajectories, those that in­ clude w, x , y, and z. Clearly, the probabilities along the remaining trajectories remain unchanged. Thus, P(w'IO) and P(wiO) will share many common factors. All these factors, including the normalization constant a , cancel out. If we let xy signify the trajectory which includes x and y, the ratio simplifies to

This ratio can be computed in time proportional to the longest network trajectory, which is constant for any given network. Note that the bottleneck of the com­ putation is the evaluation of the S(t)s, integrals over products of intrinsic noise models. Currently, all of these models are Gaussian, and so the products are Gaussian and integration is not difficult. However, some care is required when introducing other models , to ensure that the computation is still efficient. Similar problems may be introduced if the assumption of noise-free dynamics is abandoned. Without this assumption, we will need to integrate over the P(diJ)s as well as over the P(iilt)s. This can still be efficient, as long as the models are cho­ sen well.

3.3 Model updating

The dynamic variable models at each link and the intrin­ sic variable models at each sensor are all learnt using EM. A model update is performed whenever a new object is observed at a sensor, and the MCMC process has been given time to converge. Each EM iteration proceeds as follows.

The E-step requires computing joint expectations for the hidden variables w and I; with MCMC, this is approximated by sampling w from the Markov chain and extending each sample by calculating E( Jlw, 0). Again, the fact that all our intrinsic models are Gaus­ sian yields a relatively straightforward calculation that requires time proportional to trajectory length. For ex­ ample, consider a sequence of vehicle length measure­ ments with a uniform prior and a sequence of Gaussian noise models with a common variance. In this situa­ tion, the expected true vehicle length is just the mean of all the observations. If other, non-Gaussian , mod­ els are used, the computation may become much more complicated-in which case, sampling over I as well as over w may be the best approach.

The M step is exact. It uses the w and I values of each sample as if they were observed variables to perform conventional parameter learning.

Running EM to its convergence requires a Markov chain run for each EM iteration , and we run EM when­ ever a new observation is made, resulting in many Markov chain runs . In practice, the chain converges very quickly, as the models change only slowly with each in­ coming observation. We are currently investigating the use of online EM, where only a single iteration is per­ formed for each new available data point. Nowlan [1991] has proved that this approach should lead to locally max­ imum likelihood estimates in the limit , and our prelimi­ nary experimental results are encouraging.

4 Experimental results

We have performed two experiments comparing our ap­ proach to that of Huang and Russell. Our data sets were created using a freeway simulator [Forbes et al., 1995] that allowed us to control road configurations, camera characteristics, and complex vehicle behaviour. Each set included the observations generated by approximately one hundred cars, and the models used in the exper­ iments had been learnt by our EM-MCMC algorithm using data generated with the same parameter settings.

4.1 Overcoming faulty sensors

In this experiment , the algorithms were applied to a three-camera network such as that in Figure 1, set up in a manner analogous to the license-plate example of Sec­ tion 2.1. The only intrinsic attribule used was colour, and every car had a unique colour. The outer two cam­ eras , A and C, measured colour exact ly. Gaussian noise was added to the measured colour at B. Thus, the data set contained enough data for a complete matching between cameras A and C, but a pairwise appearance­ based model could be expected to lose accuracy due to its hasty independence assumptions.

The results of the experiment, averaged over three runs, can be seen in Figure 5. The y- axis requires a little explanation. In the case of the Huang and Russell algorithm, which results in discrete matches , it shows simply the percentage of correct A to C matches: in the case of the MCMC, it shows the average percentage of correct matches across all the samples drawn. The ad­ vantage of MCMC is clear: hidden feature estimation helps to maintain an almost constant degree of accu­ racy, whereas the pairwise approach is highly sensitive to sensor noise.

4.2 Estimating origin/ destination counts

Here, we compare the algorithms in the slightly more realistic setting shown in Figure 6. The aim is to esti­ mate the origin-destination counts between the two entry points and the three exit points. This requires tracking each object across the entire network. In this task, the greater number of sensors should be a liability to the pairwise algorithm, as each pair of neighbouring sensors can make its own mistakes independently of the others. Our algorithm, on the other hand, should benefit from the ability to observe each object several times to esti­ mate intrinsics more reliably.

The colour variance was manipulated as before, but at all nine cameras simultaneously. Figure 7 shows the percentage accuracy of the origin-destination counts as a function of the colour noise variance. As expected, the MCMC algorithm substantially outperforms the Huang­ Russell algorithm, holding up well even for levels of noise that essentially wipe out the colour altogether. On the other hand, Figure 8 shows that both methods are un­ able to find exact matches accurately for high levels of noise. The ability of the MCMC algorithm to re­ cover reasonable counts despite the failure of individual matches suggests that its samples contain a reasonable amount of information about the ensemble behaviour of the vehicles .

5 Summary and future work

We have described an improved approach to object iden­ tification based on the estimation of the intrinsic prop­ erties of objects and the use of Markov chain Monte Carlo to approximate the posterior probabilities effi­ ciently. We have shown that this approach works on computer-generated data, and that its computational re­ quirements are not prohibitive. We believe that this ap­ proach should be applicable in many data association applications.

We are currently working towards applying our ap­ proach in the real world . To this end , we have extended the approach to handle realistic problems such as missing observations. In the near future, we will receive data on vehicle observations made by cameras placed above the I-80 freeway as part of the Berkeley Road watch project.

Solving the tracking and data association problems for individual vehicles is only the first part of the solution for large, complex applications such as traffic surveillance. The next step is to connect these low-level computations to high-level, aggregated models, which will permit pre­ diction and control for large networks containing hun­ dreds of thousands of vehicles.