Abstract-The problem of associating data with targets in a cluttered multi-target environment is discussed and applied to passive sonar tracking. The probabilistic data association (PDA) method, which is based on computing the posterior probability of each candidate measurement found in a validation gate, assumes that only one real target is present and all other measurements are Poisson-distributed clutter. In this paper, a new theoretical result is presented: the joint probabilistic data association (JPDA) algorithm, in which joint posterior association probabilities are computed for multiple targets (or multiple discrete interfering sources) in Poisson clutter. The algorithm is applied to a passive sonar tracking problem with multiple sensors and targets, in which a target is not fully observable from a single sensor. Targets are modeled with four geographic states, tw'o or more acoustic states, and realistic (Le., low) probabilities of detection at each sample time. A simulation result is presented for two heavily interfering targets illustrating the dramatic tracking improvements obtained by estimating the targets' states using joint association probabilities. HE problems and issues involved in multi-target ocean tracking using a heterogeneous set of passive acoustic measurements were outlined in [ 1 ] ; chief among these are data association and maneuver detection. An approach to solving these problems was also described. T h e resulting experimental algorithm involves an extended Kalman-Bucy fiter (EKF) with both geographic and acoustic states, and handles measurement vectors such as bearing/frequency and delay/Doppler difference. Fundamental to the tracking algorithm is a probabilistic data association (PDA) scheme based on [ 2 ] , in which posterior association probabilities are computed for all current candidate measurements in a validation gate and used to form a weighted sum of innovations for updating the target's state in a suitably modified version of the EKF. A correction term in the Ricatti equation accounts for the uncertainty in measurement origin and the resulting state estimate is the conditional expectation, subject to the assumption that the state estimate is Gaussian prior to each update. This assumption is questionable, but it leads to quite acceptable results in practice. In addition to the basic PDA scheme, a multihypothesis Manuscript received November 2, 1982; revised February 11,1983. This work was supported by the Office of Naval Research under Contracts N00014-80-C-0270 and N00014-78-C-0529. MA 02238. T. E. Fortmann is with Bolt Beranek and Newman, Inc., Cambridge, Computer Science, University of Connecticut, Storrs, CT 06268. M . Scheffe is with the Department of Aeronautical and Mechanical Engineering, Boston University, Boston, MA 02215. structure is used for other decision making: including maneuver detection/correction, track initiation, and particularly difficult data association decisions; the entire tracking system has been run successfully in real-time experiments. However, only the PDA mechanism will be considered here. The basic PDA algorithm assumes that each target is isolated from all other targets: false measurements in a validation gate are modeled as independent clutter points drawn from a Poisson distribution with spatial density C, and detection of a target is an independent event at each sample time with probability Po. In the ocean, the false measurements often originate from other targets and cannot truly be modeled as independent Poisson clutter points. Moreover, with passive measurements, interference is probable at relatively low target densities; for example, geographically distant targets can lie at the same azimuth and be nearly indistinguishable from a single sensor. Nevertheless, our experience with real data indicates that many of these detections are so spurious that the interfering targets are effectively untrackable and, in practice, the Poisson clutter assumption often provides acceptable performance. The Poisson assumption breaks down and performance degrades significantly in the presence of discrete interfering sources, i.e., when interfering targets are detected consistently by one or more sensors. In this case, the interfering sources are usually trackable and the PDA algorithm can be extended to correct for them by computing the posterior probabilities jointly across clusters of targets; only the clutter is then modeled as Poisson. In an initial derivation [ 3 ] , some inappropriate assumptions led to erroneous prior probabilities and the combinatorics of maintaining separate validation gates for each target made it very difficult to implement. The main result of this paper is a new joint probabilistic data association (JPDA) algorithm for multiple targets in Poisson clutter. This is a target-oriented approach, in the sense that a set of established targets is used to form gates in the measurement space and to compute posterior probabilities, in contrast to the measurement-oriented algorithms of Reid [4] and others, where each measurement is considered in turn and hypothesized to have come from some established track, a new target, or clutter. Although these two points of view can be shown (with appropriate assumptions) to yield equivalent expressions for the posterior probabilities, the difference between them is more than just philosophical. Neither position nor velocity is fully observable from a single passive measurement and clutter densities in the ocean environment are typically quite high; these two characteristics make it essentially intractable to hypothesize a new track for each individual measurement. Indeed, in our approach, targets are initiated by a separate operator-interactive process using composite groups of measurements. This is also a nonbackscan (or zero-scan) approach, meaning that all hypotheses are combined after computation of the probabilities, for each target at each time step. While it can be extended to retain n scans as in [5], this would lead in the passive sonar environment to enormous memory requirements because of the relatively high clutter densities (0.2-2.0 false detections per gate) and the large number of separately measured variables (2-12 validation gates per target) that are typically encountered. Moreover, because of the spurious nature of many of the false detections, we believe that the JPDA is more cost-effective than a brute-force n-scan algorithm in this environment. Retaining multiple scans would also cause enormous complications for the maneuver-detection and hypothesis-testing machinery that is built on top of the JPDA tracker. Finally, this is a recursive approach, suited to the efficient maintanance of established tracks, whereas batch-oriented approaches, such as [12], are suitable for initial track formation in an unstructured set of data from a limited time interval and/or from a subspace of the measurement space. The two approaches are complementary and it is appropriate for such a batch algorithm (or the operator-interactive process mentioned above) to initiate tracks and pass them off to the JPDA algorithm for continued data association and tracking. A preliminary version of this paper appeared in [ l o ] . The issues of unresolved measurements and target maneuvers are discussed in [I61 and an algorithm is described that extends the work presented here. Another interesting extension may be found in [ 171, which quantifies the dependence of PDA tracking performance on a detection threshold in the signal processing algorithm that supplies its measurements. The multitarget data association problem is formulated in the next section, followed by a derivation of the new JPDA algorithm in Section 111. A specific application in passive sonar surveillance is described in Section IV, and the performance improvements gained with the JPDA algorithm in this application are illustrated in Section V, using simulated data from two heavily interfering targets in clutter. Consider a dynamic system (target model) of the familiar form where x is the target state vector, y is the measurement vector, w and u are zero-mean mutually dependent white Gaussian noise vectors with covariance matrices Q and R , respectively, and k is a discrete time index. The matrices F , G , H , &, and R are assumed known and their dependence on k is suppressed here for notational convenience. The initial state is assumed Gaussian with mean ;ole and covariance Polo. A specific target model is described below in Section IV.  where the innovation vector has the covariance matrix and the filter gain matrix is The resulting state estimate, under the above assumptions, is the conditional mean where Yk denotes the set of all data vectors y j for i < k . In order to avoid cluttering the discussion to follow, this brief summary ignores a number of complications that arise in practice. If the system is nonlinear, for instance (as it is in the example presented below), then it can usually be linearized and the same basic equations can be applied to deviations from the nominal trajectory [7], [9]. If the target occasionally deviates from the assumed motion model, e.g., by maneuvering, then some decision-making or other machinery must be provided to deal with these instances. A fundamental characteristic of this class of problems is that the size and composition of the measurement vector are unpredictable from one time to the next; in other words,yk comprises a time-varying set of independent subvectors, as discussed below in Section IV. We shall avoid the resulting notational morass by restricting (2.3)-(2.6) to apply to a single measurement subvector y k from a single sensor. In addition, we will suppress the time subscript k from all variables except P and Y , unless it is required for clarity. Without any loss of generality, the data association problem may now be formulated as follows. At each time step, the sensor provides a set of candidate measurements to be associated with targets (or rejected). In most approaches, this is done by forming a "validation gate" around the predicted measurement from each target and retaining only those detections that lie within the gate. There are  and the corresponding innovations are The term measurement will be used interchangeably for yi and j j j , since they contain equivalent information [9]. denotes the event that the jth measurement belongs to that target and x . the event that none of the measurements belongs to it (no detection). The PDA approach builds upon the assumptions that the estimation errors $ and p have Gaussian densities at each time step (this is approximate, since there is an exponentially growing tree of possible measurement sequence hypotheses and the true densities are weighted sums of Gaussians). It is also assumed that the correct measurement is detected with probability Po (independently, at each time) and that all other measurements are Poisson-distributedl with parameter CV, where C is the expected number (or density) of false measurements per unit volume and Vis the volume of the validation gate. The gate is normally a "g-sigma" ellipsoid { y :?IS-'? <2) and PC is the probability that the correct measurement, if de- 1 Equivalently, the number n of false measurements has probability mass function p ( n ) = e-cv(CVp/n! and the location of each false measurement is uniformly distributed in the gate. 2 This is just the Gaussian probability mass in the gate, which is M1/2 + 2. often assumed to be unity in practice, since PC > 0.99 whenever g > accounts for the measurement origin uncertainty. Note that the data-dependent factors Po and Pk transform the original deterministic Riccati equation into a stochastic one. The equations above define the PDA filter for a single target, and additional targets may be handled with multiple copies of the filter. However, with respect to any given target, measurements from interfering targets do not behave at all like the random (Poisson) clutter assumed above. Rather, the probability density of each candidate measurement must be computed based upon the densities of all targets that are close enough (when projected into the measurement space) to interfere. In order to account for this interdependence, consider a cluster3 of targets (established tracks) numbered t = 1, e -, T at a given time k. The set of m candidate measurements associated with this cluster (i.e., found in the validation gates for targets 1. -.. T ) is denoted yi, j = 1 ~ .e., m , as above. Each measurement belongs either to one of the T targets or to the set of false measurements (clutter), which is denoted by target number t = 0. Denoting the predicted measurement for target t b y i f , t h e innovation (2.10) corresponding to measurement j becomes where pit is the posterior probability that measurement j originated from target t and pot is the probability that none of the measurements originated from target t (i.e.: it was not detected). This is used in target t's copy of (2.3) to update the state estimate 5'. In other words, the joint probabilistic data association (JPDA) and PDA approaches utilize the same estimation equations; the difference is in the way the association probabilities are computed. Whereas the PDA algorithm computes pit, j = 0, 1, -0 , m, separately for each t, under the assumption that all measurements not associated with target t are false (i.e., Poisson-distributed clutter), the JPDA algorithm jointly across the set of T targets and clutter. From the point of view of any target, this accounts for false measurements from both discrete interfering sources (other targets) and random clutter. Details are given in the next section. The key to the JPDA algorithm is evaluation of the conditional probabilities of the following joint events: 3 A cluster is a set of targets whose validation gates are "connected" where xiri is the event that the measurement j originated from target ti, 0 < ti < T , and ti is the index of the target to which measurement j is associated. The feasible events are those joint events in which no more than one measurement originates from each target, i.e., It is also convenient to defme the measurement association indicator which indicates whether measurement j is associated with any established target in event x, and the target detection indicator which indicates whether any measurement is associated with target t i n event x (Le. whether target tis detected). For the purpose of deriving expressions for the joint probabilities, no individual validation gates will be assumed for the various targets in a cluster. Instead, each measurement will be assumed validated for each target, i.e., every validation gate coincides with the entire surveillance region. This approach is adopted because the resulting equations are equivalent to but simpler than those derived if different validation gates are assumed for each target. The extra computational burden resulting from the consideration of events with negligible probability can be avoided with suitable logic limiting the probability calculations to events involving validated measurements. This logic, which has a negligible effect on the numerical results, is presented next, followed by a derivation of the joint event probabilities. Finally, the marginal event (measurement-to-target association) probabilities needed by the tracking filter are obtained from the joint event probabilities. The following validation matrix is defmed with binary elements to indicate if measurement j lies in the validation gate for target t. Index t = 0 stands for "no target" and the corresponding column of Sl has all units-each measurement could have originated from clutter or false alarm. A typical set of gates and the corresponding validation matrix is shown in Fig. 1 . Each event x may be represented by a matrix consisting of the units in S2 corresponding to the associations Given considerable patience? one can coqstruct a procedure that scans S2 and enumerates the matrices Q(x) corresponding to feasible events using the following rules: 1) Scan R by rows and pick one unit per row for Sl (i.e., there can be only one origin for a measurement). 2) Only one unit from each column t 2 1 can be taken (i.e., at most one measurement could have originated from a target). The number of units from column t = 0 is not restricted. Using Bayes' rule, the probability of a joint event conditioned on all measurements up to the present time is P ( x I Y~} = P ( X I~~ , . . . , y ' r n , m , y k -l } =~(F~:--,Y,,, I X ,~, yk-lY{xlm, yk-lYc. Using the same assumptions as in the PDA approach, a measurement y j associated with a target t has a Gaussian density5 where Fit is the innovation defined in (2.14) and S, is its covariance. A measurement associated with clutter has a uniform density V-', where V is the volume of the entire surveillance region (note that in Section I1 and the Appendix, V was the volume of a single gate). It follows that the densities of individual measurements in (3.1 1) are given by (3.13) The second factor in (3.10) is the prior (to time k ) probability of a joint event. To obtain it, note that the total number of false measurements in event x is m (3.14) Substituting (3.1 l), (3.13), and (3.16) into (3.10) yields (3.17) Note that @! cancels immediately in (3.17). Furthermore, m ! and e-" appear in (3.1 7) regardless of which event is considered. Since the denominator term c from (3.17) is the sum of all the numerators, it follows that m! and e- '" cancel out as well. The final expression is obtained as where c is the (new) normalization constant. Numerical overflows and underflows are common in the factor Co and the normalization constant c, because the magnitude of C (in units of I/volume in the measurement space) is quite variable and @ can be 10 or more. The problem can be avoided simply by letting 1 /C be the unit volume in calculating (3.18), so that C@is replaced by lo. This change canceis out in the exponential factor, and it causes I Srj I' 1 ' in the denominator to be multiplied by C. Alternatively, one can use logarithms to compress the numerical range. With this, the prior probability of event x is given by P{xIm,Y"-lI where Pot is the probability of detection of target t, the number @ of false measurements is again assumed Poisson distributed with parameter CV, Cis the density of false measurements, and the dependence of @ and 6 t on x has been dropped from the notation. 5 The Gaussian density defmed in (A.7) of the Appendix includes a factor ~/ P G to account for its restriction to the validation gate. Here we assurnePG = 1 (corresponding to a very large gate) for convenience.  Passive sonar tracking is a very challenging application for data association and tracking algorithms; most of the diffculties and idiosyncrasies that arise in the ocean environment have been discussed above and in [ I ] . In this section, we shall briefly describe a particular target/measurement model that has been used successfully in this context. Tracking results using this model in the JPDA filter are presented in the next section. We assume that each target follows a nominal straight-line trajectory with random disturbances (great circle motion can easily be substituted if necessary) and radiates acoustic energy in one or more characteristic frequency bands. Since much of the acoustic energy is produced by the propulsion system (engine, gears, shafts, propellers, etc.), the corresponding frequencies are coupled with target speed. These relationships are summarized for a particular target by the following state equations:6 We consider a set of passive sensors, each of which detects incoming "signals" (acoustic energy) and estimates their bearings (arrival angles) and center frequencies. Thus a single-sensor measurement of one signal from the target described by (4.1) has the form bearing frequency where n is a sensor index, the target's bearing f i n ( -) and aspect angle a, (.) are standard trigonometric functions, Fk is one of associated with a particular target's state consists of subvectors of th'e form (4.2) and (4.3), corresponding to various combinations of sensors, source frequencies, and harmonics. In practice, the received measurements corresponding to these subvectors have little or no a priori association information to link them, and it is appropriate to deal with each one separately. Moreover, each measurement may be absent (undetected) and may be accompanied by false measurements from different targets and/or clutter: this is the essence of the data association problem discussed in Section 11. Equations (4.1) and (4.4) are easily lineariied and sampled' at discrete times, yielding the target model (2.1)-(2.2). Tracking results using this model in the JPDA fdter on simulated data are given in the next section. A simulation program was used to create realistic passive sonar data on which the JPDA algorithm could be tested. One such data set contains measurement files of bearinglfrequency liries and time/Doppler differences for two hypothetical targets, with a common 12-Hz source frequency and courses that result in severe interference. The target-sensor geometry is indicated in Fig. 2 . Targets 1 and 2 travel at 6 kn on courses of 100' and 800, respectively, and cross midway through the 6-h period shown. Measurement data were created by dead reckoning target 6 The symbols C, V , M , f, k, e, m , and p used in this section are un-7 In practice, measurement data are usually integrated over an aprelated to those used previously. propriate period, rather than sampled, in order to improve the SNR. motion (no process noise), adding noise to the computed true measurements, and then adding clutter measurements (false detections). The standard deviation of the measurement noise was 5' for bearings, 80 mHz for frequency, 3.6 s for time difference, and 4 mHz for Doppler difference. The true measurement was detected at any given time with probability Po = 0.7. The number of clutter points was Poisson-distributed and their locations in the measurement space were uniformly distributed over a very broad region about the actual track. The clutter density was C = 0.25/(Hz.deg) for bearing/frequency and C = 0.25/(Hz.s) for time/Doppler difference; with 'varying gate sizes, this ranged from 0.2 to 2.0 false detections per gate. Several tracks were made using this data; the initial 2-sigma confidence ellipse areas were 2% nmi2, with course and speed confidences of +20° and ?3 kn, respectively. In practice, one would use somewhat looser limits to initialize uncertain targets, but the object here was to simulate a situation in which the tracks are already well-established before they intersect. Although the data had no process noise, the filter was given process noise standard deviations of 0.2 kn in speed, 0.2" in course, and 0.01 Hz in source frequency @ order to avoid divergence. Fig. 3 shows tracks made using the well-known "nearest neighbor" data association scheme [ 131 , in which the candidate measurement closest to the center of each validation gate is accepted as correct, all others are ignored, and the covariance equation is not augmented to account for association errors. As one might anticipate, the tracker becomes hopelessly confused and loses one of the targets altogether. The ordinary PDA filter (without multitarget logic) fares only slightly better in this situation. As shown in Fig. 4 , both tracks lock onto a sort of "compromise" and end up lost midway between the two actual targets. This behavior is to be expected, since a basic assumption in the standard PDA algorithm (applied to each target in turn) is that measurements not originating from the target under consideration are random clutter points with a Poisson/uniform distribution. The joint PDA method, described above in Section 111, corrects this situation by allowing the probabilistic weights used in data association to be computed joiyltZy across all known targets. As shown in Fig. 5 , this improves the tracks dramatically, and the 2-sigma confidence ellipses contain the true position in all cases. Note that the ellipses are larger than those of Fig. 4 , particularly near the point of intersection, reflecting the fact that the algorithm "hedges" its decisions and relies more on dead reckoning when the targets are too close to determine with high confidence which measurement belongs to which target. The JPDAcomputed association probabilities of a measurement belonging to the correct target versus the other one typically have ratios of 0.8/0.2 or greater early in the run. As the targets draw nearer each other, the ratios switch rather abruptly to around 0.5/0.5 (and the confidence ellipse areas increase accordingly); then they switch back again when the targets separate. JPDA outperforms the other methods on the basis of how often the correct measurement receives the highest probability [lo], but the difference is not nearly as great at Figs. 3-5 would indicate. These facts suggest that the real value of the JPDA approach lies in the relative probabilities computed at key points in the target trajectories, rather than in the overall percentage of correct binary "decisions." Finally, tracks made using perfect data association are shown in Fig. 6 . Comparing these with Fig. 5 , it is clear that the JPDA tracker performs remarkably well despite the severe interference between targets. Indeed, by the end of the 6-h period, the two sets of tracks are nearly identical. As important as tracking accuracy is the tracker's own assessment of that accuracy, expressed by the covariance matrix. To compare this for the tracks of Figs. 3-6 , the filter consistency measure was computed at each time k , where x k is the true state, ik is the estimate, and Pk is the covariance matrix computed on-line by the tracker. Mean-square and root-mean-square (rms) values of this statistic over the 6-h period (7 1 time steps) are given in Table I , together with tlieir'theoretical values (r2 is a chi-square random variable with 6 degrees of freedom). Table I indicates that the ordinary PDA tracker is overoptimistic, having actual errors much larger than its calculated covariance would indicate. The JPDA tracker is about as consistent as the perfect data association tracker and both are slightly pessimistic in comparison to the theoretical chi-square value. Also shown in Table I are the final position errors (in nautical miles) for the same period. Fortuitously, the fiial position error for Target 2 using the ordinary PDA is slightly better than the others. Nevertheless, the confusion of the ordinary PDA is evident from Fig. 4 . Note that the above results represent only a single "sample" from the space of all possible noise/clutter sequences and target/sensor geometries; a Monte Carlo approach involving many such samples would provide a statistic51 basis for evaluating the algorithm. However, the extensive simulation and realistic clutter densities used here are moderately costly and Monte Carlo simulations were not feasible with the available r'esources. results of this paper demonstrate the importance of correlating (i.e., associating) each report with all tracks simultaneously, rather than separately. The joint probabilistic data association framework is particulary well-suited to this task in the passive sonar environment, where additional measurements arise from both interfering targets and random clutter, detection probability is significantly less than 1, the measurement space is large and heterogeneous, a target's state is not fully observable from any single measurement, and maneuvers are common. Continuing research in this and related areas includes assessment of tracking performance as a function of detection probability and clutter density [17] (this has important implications for the setting of gains and thresholds in the signal processing algorithms that provide measurements to the tracker), extension of the likelihood function and other hypothesistesting machinery (e%., for maneuver detection) to the JPDA situation, the use of more complex models for the detection and clutter processes, and automated schemes for track initiation. In this Appendix we shall derive expressionsfor the posterior association probabilities flj in (2.1 1). Letting rj(m) denote the prior probability of the event xi, conditioned on m, the total probability theorem yields  The first factor is the joint probability density of the m candidate measurements, conditioned on the jth one being correct. According to the PDA assumptions, the correct measurement yj has a Gaussian density with mean 0 and covariance S, where the factor 1/PG accounts for its restriction to the validation gate and each incorrect measurement has a uniform density Y-l. It follows that The second factor in (Ah) is the prior probability of xi, given by (AS). The denominator is the joint prior density of the measurements, conditioned only on m (and the past data) Note that with the above conditioning, the validated measurements are not independent, i.e., (A.9) is not equal to the prod-
More complex approaches, such as Joint Probabilistic Data Association Filter (JPDA) REF combines all of the potential measurements into one weighted average, before associating it to the track, in a single update.
Sonar tracking of multiple targets using joint probabilistic data association
IEEE Journal of Oceanic Engineering
Engineering