I. INTRODUCTION Comparison of Recursive Algorithms for Emitter Localisation using TDOA Measurements from a Pair of UAVs NICKENS OKELLO, Member, IEEE NICTA Australia FIONA FLETCHER DSTO Australia ² DARKO MUSICKI, Member, IEEE Hanyang University BRANKO RISTIC DSTO Australia This paper presents a comparative analysis of three nonlinear filters for estimation of the location and velocity of a moving emitter, using time difference of arrival (TDOA) measurements received by two unmanned aerial vehicles (UAVs) as they traverse the surveillance region. The TDOA measurements are generated over time by comparing and subtracting leading edge time of arrivals (TOAs) of signals. The Cramér Rao lower bound (CRLB) of estimation errors is derived and used as the benchmark in performance analysis. The three nonlinear filters considered in the comparison are: a Gaussian mixture measurement integrated track splitting filter (GMM-ITSF), a multiple model filter with unscented Kalman filters (UKFs) and a multiple-model filter with extended Kalman filters (EKFs). Manuscript received March 8, 2009; revised April 21, 2010; released for publication June 1, 2010. IEEE Log No. T-AES/47/3/941759. Refereeing of this contribution was handled by W. Koch. This work was completed under the Research Contract “Advanced Methods of Tracking and Fusion,” between ISR Division of DSTO and The University of Melbourne. D. Mu²sicki was partially supported by FGAN-FKIE under Contract 4500035757. Authors’ current addresses: N. Okello, National ICT Australia (NICTA), Victoria Research Laboratory, Dept. of Electrical & Electronic Engineering, The University of Melbourne, Parkville, VIC 3010, Australia, E-mail: (nickens.okello@nicta.com.au); F. Fletcher, ISR Division, DSTO, PO Box 1500, Edinburgh, SA 5111, Australia; D. Mu²sicki, Dept. of Electronic Systems Engineering, Sa 3 Dong 1271, Hanyang University, Room 407, 3rd Eng. Building, Ansan, Kyunggido 426-791, South Korea; B. Ristic, ISR Division, DSTO, 506 Lorimer St., Melbourne, VIC 3207, Australia. c 2011 IEEE 0018-9251/11/$26.00 ° Emitter geolocation using unmanned aerial vehicles (UAVs) is an important application [1]. The UAVs have passive sensors that are only able to measure the time of arrival (TOA) of signals at the receiver. A single sensor of this type is unable to infer any emitter location information from this measurement without knowing the time of emission of the signal. However, when two such sensors are available, and subjected to the same signal, the time difference of arrival (TDOA) may be calculated. In the absence of measurement noise, this TDOA defines a hyperbola in two-dimensional space1 on which the emitter must be located. When measurements from n such sensors are available from a given emission, the emitter position is found as the intersection of n ¡ 1 independent hyperbolae defined by the TDOA measurements formed by selecting a reference sensor (referred to as sensor 1 in the remainder of the paper) and then pairing each of the remaining sensors with the reference sensor.2 For a two-dimensional geolocation problem, only three sensors are required to find the emitter position, provided that the sensors are not located collinearly with the emitter [2, 3]. When measurements are noise corrupted, the hyperbolae will no longer intersect exactly at the emitter location. An estimate of emitter location can be found using least squares [4, 5]. The presence of further additional sensors can help improve the estimation of emitter location. Examples of applications include mobile positioning of wireless communication systems [6, 7] and geolocation of stationary or moving emitters using multiple UAVs [1, 4, 5, 8, 9]. These emitters may be tracked over time from position estimates calculated using these localisation techniques at each emission as measurements in a linear Kalman filter or other smoothing technique [6, 7, 9]. These localization techniques together with many others [10—14] can be classified as traditional methods in that they employ at least three noncollocated sensors that collect TDOAs of a single pulse from which emitter location estimate is obtained. This paper considers a scenario with two noncollocated mobile sensors where sensors measure TDOA over a number of emissions from a possibly moving emitter. The time history of received measurements, coupled with the assumption that emitter motion model is known (in the paper we assume constant velocity emitter motion), is sufficient to estimate target location over time. 1 In three dimensional space, the TDOA defines a hyperboloid on which the emitter must be located. 2 In practice, there are n(n ¡ 1)=2 pairings of sensors, but only n ¡ 1 are independent. IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011 1723 The batch localization techniques, such as maximum likelihood and least squares estimation [4, 5] are restricted to deterministic emitter motion models. These techniques are also typically computationally intensive, and if an updated estimate of target location is required each time a signal is received, the processing must be repeated for each measurement. Hence a recursive estimation procedure is preferred. Since the TDOA is a nonlinear measurement, nonlinear filtering techniques are considered in this paper. Sathyan et al. [15] investigate the problem of tracking an unknown number of emitters using TDOA measurements. They use an interacting multiple model tracker with unscented Kalman filters (UKFs) for estimation, but focus on solving the data association problem. They do not specifically address the problem of tracking with only two UAVs as sensors, and do not analyse the tracking performance separate from the data association procedure. The problem of tracking stationary or moving emitters using multiple TDOA measurements collected over time by two instead of three UAVs was tackled by Okello and Mu²sicki [8]. In their approach, the hyperbolic measurement error region is approximated by a sum of weighted Gaussian distributions. Furthermore, an integrated track splitting (ITS) framework is used for state estimation, and some form of pruning or merging is required to avoid a combinatorial increase in the number of components comprising the estimate at any time due to the number of measurement components. They also describe a bank of UKFs initiated with Gaussian mixture measurement (GMM) components but using the nonlinear measurement process after initialisation. However they do not provide detailed analysis of the performance of either algorithm. Later this algorithm was extended to consider ambiguous association of measurements between UAVs [16]. Mu²sicki et al. [17] have extended this algorithm to include frequency difference of arrival in addition to TDOA. Fletcher et al. [18] have investigated this problem using a simple extended Kalman filter (EKF) and UKF, without considering the track initiation issues. In this paper we derive the Cramér Rao lower bound (CRLB) for this problem, and compare it with the performance of three algorithms for tracking using TDOA measurements obtained by two UAVs. The first filter uses the GMM approximation to linearize each measurement before applying the result to a dynamic Kalman filter bank. The second and third tracking algorithms are the static bank of EKFs and UKFs, respectively. The paper is organized as follows. The problem statement is presented in Section II. Section III describes the method of approximating the likelihood of a TDOA measurement by a Gaussian mixture. Section IV describes three candidate tracking 1724 algorithms that will be used in the comparison. The calculation of a CRLB is also discussed in Section IV. Finally, simulation studies are presented in Section V and conclusions in Section VI. II. PROBLEM DESCRIPTION The emitter localisation problem considered here is based on two networked UAVs which measure TDOA of received signals. For simplicity, we assume that the UAVs and emitters are coplanar and use two-dimensional models, but the techniques presented may be extended to higher dimensional spaces. The UAVi, i = 1, 2 is located at Cartesian position (i) (i) s(i) k = (´k , ³k ), where k is a discrete-time index which refers to the instances when TDOA measurements are formed. UAV locations are assumed known. The emitter is a stationary or moving radar located at (xk , yk ) with one main beam and negligible side lobes. Times of arrival at each UAV are measured with additive errors, which form a zero-mean white Gaussian sequence with a standard deviation of ¾t . The measured arrival time at the ith UAV is tk(i) = tk + rk(i) =c + n(i) k , i = 1, 2 (1) where tk is the transmission time, rk(i) = q (xk ¡ ´k(i) )2 + (yk ¡ ³k(i) )2 is the Euclidean distance between the emitter and UAVi, and n(i) k is a zero-mean Gaussian measurement noise with standard deviation ¾t . The emitter position and the time of emission of the signal tk are unknown. The times of arrival measurements to two UAVs are subtracted and the TDOA measurement ¿k is ¿k = tk(2) ¡ tk(1) = (rk(2) ¡ rk(1) )=c + nk (2) (1) 2 where nk = n(2) k ¡ nk » N (n; 0, 2¾t ) and N (n; m, P) denotes the Gaussian probability density function (pdf) of variable n with mean m and variance P. A single noiseless TDOA measurement defines a hyperbola on which the emitter must lie and the additive measurement noise defines an uncertainty area around this hyperbola. It is convenient to use the range difference of arrival (RDOA), calculated by multiplying the TDOA by the speed of signal propagation c (the speed of light in the case of electromagnetic emissions). The measurement equation becomes zk = c ¢ ¿k = (tk(2) ¡ tk(1) )c = rk(2) ¡ rk(1) + wk = hk (xk ) + wk (3) (1) hk (xk ) = k¡ xk ¡ s(2) k k ¡ k¡ xk ¡ sk k (4) where wk is white zero-mean Gaussian, wk » N (w; 0, ¾w2 ), ¾w2 = 2c2 ¾t2 , xk is the emitter state vector defined by its position and velocity in the Cartesian coordinates, IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011 xk = [xk x_ k yk y_ k ]T , while ¡ xk = [xk , yk ]T is the Cartesian position of the emitter (¡ projects the emitter state into the Cartesian position space). The emitter motion model is xk = Fk¡1 xk¡1 + vk¡1 where 2 1 ±k 60 6 Fk = 6 40 0 1 0 0 0 0 3 07 7 7 1 ±k 5 0 0 (5) (6) 1 is the state transition matrix, vk¡1 is white zero-mean Gaussian process noise, vk¡1 » N (v; 0, Qk¡1 ), independent of wk , and ±tk is the time between consecutive emissions that were used to form TDOA measurements, which we approximate by the time between consecutive TOA measurement stages at (1) UAV1, ±tk ¼ tk(1) ¡ tk¡1 , as the TDOAs of a given pulse is usually very small compared with the time interval between measurement stages. The emitter periodically transmits a signal and if both UAVs are within the beamwidth of the emitter and are able to receive the same signal transmission, then a TDOA measurement can be formed. We assume that the TDOA is small compared with the time interval between consecutive emitter transmissions so that the signals arriving at each sensor can be unambiguously associated. The sensors may use any motion model, provided that their positions are known at the arrival time of the wavefront. Given a sequence of TDOA measurements Z k = hz1 , z2 , : : : , zk i, our goal is to recursively estimate the posterior pdf p(xk j Z k ) from which we can then form estimates of the state vector x̂k . From (3), the likelihood function of a RDOA measurement zk is (7) Likelihood (7) is a function of ¡ xk and defines a hyperbolic uncertainty region in the Cartesian emitter position space ¡ xk . We approximate likelihood (7) by (up to the multiplicative constant) [17] p(yk j xk ) = Gk X °k (g)N (ŷk (g); ¡ xk , Rk (g)) g=1 where each element g is termed a measurement component, Gk denotes the number of measurement components, and relative measurement component probabilities satisfy °k (g) ¸ 0. In other words, we replace measurement zk by a set of triplets °k (g), ŷk (g), Rk (g), with g = 1, : : : , Gk . IV. ALGORITHMS FOR TRACKING WITH TDOA MEASUREMENTS All three candidate algorithms adopt a Gaussian mixture representation of the posterior pdf [19]. Each element of the emitter state pdf is termed a “track component.” At time k ¡ 1 (after processing measurement zk¡1 ) the posterior pdf is given by p(xk¡1 j Z k¡1 ) = Jk¡1 X »k¡1 (j)N (xk¡1 ; x̂k¡1jk¡1 (j), Pk¡1jk¡1 (j)) j=1 (9) where j is the index of track components, Jk¡1 denotes the number of track components, x̂k¡1jk¡1 (j) and Pk¡1jk¡1 (j) are the mean and covariance of track component j, and »k¡1 (j) is the weight (probability) of track component j at time k ¡ 1 and satisfies Jk¡1 X »k¡1 (j) ¸ 0: »k¡1 (j) = 1, j=1 III. TDOA MEASUREMENT GAUSSIAN MIXTURE APPROXIMATION p(zk j xk ) = N (z; h(xk ), ¾w2 ): For each RDOA measurement zk , zk ¡ ¾w and zk + ¾w define two hyperbolae which bound the §¾ emitter location region. By segmenting the region bounded by these two hyperbolae, we generate a set of two-dimensional measurement components. Each measurement component is represented by a Gaussian pdf whose mean is the centre of the segment and whose covariance is defined by a 1-¾ ellipse that is bounded by the segment [17]. Relative measurement probability is proportional to the ellipse area, °k / p jRk (g)j. The described procedure is referred to here as the GMM approximation. An application of GMM approximation is illustrated in Fig. 3. (8) (10) The algorithms are initialized identically, using the first TDOA/RDOA measurement z1 and the GMM technique described in Section III. The GMM technique forms G1 triples f°1 (j), ŷ1 (j), R1 (j)g, j = 1, : : : , G1 , corresponding to z1 . Each track component j = 1, : : : , J1 = G1 is then initialized by x̂1j1 (j) = [ŷ1 (j)[1] ºx ŷ1 (j)[2] ºy ]T 2 R1 (j)[1, 1] 0 R1 (j)[1, 2] 6 0 ¾v2 0 6 P1j1 (j) = 6 4 R1 (j)[2, 1] 0 R1 (j)[2, 2] 0 0 0 (11) 0 3 07 7 7 05 ¾v2 (12) where y[i] denotes the ith element of vector y, and similarly Y[i, j] the (i, j)th element of matrix Y. Initial velocities ºx and ºy in (11) are typically zeros, with variance ¾v2 in (12) chosen based on the expected maximum emitter velocity. The initial OKELLO ET AL.: COMPARISON OF RECURSIVE ALGORITHMS FOR EMITTER LOCALISATION 1725 track component probabilities are »1 (i) = °1 (i) for i = 1, : : : , J1 = G1 . Due to high uncertainty of the initial location of the emitter, as depicted in Fig. 3, standard non-Gaussian-mixture-based nonlinear estimators are likely to suffer from track divergence. This includes the EKFs and UKFs. A particle filter implementation is also an option, however as the algorithms presented here converge to the Cramer-Rao bound, it is deemed to be an overkill. A. Gaussian Mixture Measurement–Integrated Track Splitting Filter In GMM-ITSF both the state estimate and the observation space measurement pdfs are modelled by Gaussian mixtures. Let the posterior pdf at time k ¡ 1 be given by (9). The target state transition is assumed linear, and thus the Chapman-Kolmogorov equation [20] is Jk¡1 X »k¡1 (j)N (xk ; x̂kjk¡1 (j), Pkjk¡1 (j)) p(xk j Z k¡1 ) = j=1 subject to Jk¡1 Gk X X »k (j + ) = 1: (22) j=1 g=1 The output of the GMM-ITSF can be either the mean and the covariance of the highest probability component j¤+ = arg maxj + »k (j + ), or the combined mean and covariance [21] of the Gaussian mixture X x̂kjk = »k (j + )x̂kjk (j + ) (23) j+ Pkjk = X j+ »k (j + )(Pkjk (j + ) + x̂kjk (j + )x̂Tkjk (j + )) ¡ x̂kjk x̂Tkjk : (24) If left unchecked the number of track components increases exponentially which quickly exhaust computational resources. Thus track component management in the form of pruning and/or merging is a prerequisite for any practical implementation of the GMM-ITSF. (13) B. A Bank of Extended Kalman Filters where each component j propagates in a linear manner x̂kjk¡1 (j) = Fk¡1 xk¡1jk¡1 (j) (14) Pkjk¡1 (j) = Fk¡1 Pk¡1jk¡1 (j)FTk¡1 + Qk¡1 : (15) The TDOA or RDOA measurement is converted to the Cartesian emitter location space using the GMM approximation of Section III. Each measurement component g in (8) uses the Kalman filter update of each track component j to create a new component of the a posteriori track pdf. The number of a posteriori track components at time k is Jk¡1 £ Gk and the a posteriori track state pdf is Jk¡1 Gk k p(xk j Z ) = X + + + »k (j )N (xk ; x̂kjk (j ), Pkjk (j )) j + =1 (16) + where j ´ fg, jg is the index into a posteriori track components, and x̂kjk (j + ) = x̂kjk¡1 (j) + Kk (j + )(ŷk (g) ¡ ¡ x̂kjk¡1 (j)) (17) + + + + T Pkjk (j ) = Pkjk¡1 (j) ¡ Kk (j )Sk (j )Kk (j ) (18) The proposed algorithm consists of a bank of EKFs (EKFB), each contributing a track component in (9). The prediction step of each EKF follows the Kalman filter equations (14) and (15). The jth track component is updated at time k by [21] x̂kjk (j) = x̂kjk¡1 (j) + Kk (j)(zk ¡ hk (x̂kjk¡1 (j))) (25) Pkjk (j) = Pkjk¡1 (j) ¡ Kk (j)Sk (j)Kk (j) T where Sk (j) is the innovation covariance, and Kk (j) is the filter gain defined by Sk (j) = Hk (j)Pkjk¡1 (j)Hk (j)T + ¾w2 (27) Kk (j) = Pkjk¡1 (j)Hk (j)T Sk (j)¡1 : (28) Matrix Hk (j) is the Jacobian of the measurement function (4), evaluated at the predicted emitter state, ¯ @hk (x) ¯¯ Hk (j) = @x ¯x=x̂kjk¡1 (j) à ! T T (¡ x̂kjk¡1 (j) ¡ s(1) (¡ x̂kjk¡1 (j) ¡ s(2) k ) k ) ¡ ¡: = k¡ x̂kjk¡1 (j) ¡ s(1) k¡ x̂kjk¡1 (j) ¡ s(2) k k k k (29) where Kk (j + ) = Pkjk¡1 (j)¡ T Sk (j + )¡1 (19) Sk (j + ) = ¡ Pkjk¡1 (j)¡ T + Rk (g): (20) + Relative probabilities »k (j ) ´ »k (fg, jg) satisfy »k (j + ) / »k¡1 (j)°k (g)N (ŷk (g); ¡ x̂kjk¡1 (j), Sk (j + )) (21) 1726 (26) The likelihood of the jth track component is evaluated as [21] ¤k (j) = N fzk ; hk (x̂kjk¡1 (j)), Sk (j)g: (30) Track component probabilities are updated according to » (j)¤k (j) : (31) »k (j) = PJk k¡1 j=1 »k¡1 (j)¤k (j) IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011 As in the GMM-ITSF, the output of the EKFB can be either the mean and the covariance of the highest probability track component j ¤ = arg maxj »k (j), or the combined mean and covariance of the Gaussian mixture based on (23) and (24). The number of track components in the EKFB Jk = Jk¡1 . C. A Bank of Unscented Kalman Filters The UKF bank (UKFB) propagation is identical to the EKFB propagation. During the update of filter j, a set of sigma points, Â(i) kjk¡1 (j), are generated from the predicted state by 8 x̂kjk¡1 (j) i=0 > > > > ´T ³ q > > > > x̂kjk¡1 (j) + (nx + ·)Pkjk¡1 (j) > > i < (i) i = 1, : : : , nx Âkjk¡1 (j) = > > > ´T ³q > > > x̂ (n + ·)P (j) (j) ¡ > x kjk¡1 kjk¡1 > > i > : i = nx + 1, : : : , 2nx (32) p where ( A)i denotes the ith row of the matrix square root L, such that A = LT L and nx is the dimension of the state vector. The corresponding weights are 8 · i=0 > < nx + · , (i) (33) Wkjk¡1 (j) = > : 1 , i = 1, : : : , 2nx nx + · as described in [22]. The predicted measurement is zkjk¡1 (j) = 2nx X (i) Wkjk¡1 (j)hk (Â(i) kjk¡1 (j)): (34) i=0 The update of track component j is x̂kjk (j) = x̂kjk¡1 (j) + Kk (j)(zk ¡ zkjk¡1 (j)) (35) Pkjk (j) = Pkjk¡1 (j) ¡ Kk (j)Sk (j)Kk (j)T (36) where Kk (j) = Pxz (j)=Sk (j) Sk (j) = Pxz (j) = Pzz (j) + ¾w2 2nx X i=0 2nx X i=0 (38) D. Cramér-Rao Lower Bound The CRLB [23], defined as the inverse of the Fisher information matrix (FIM) Jk defines theoretical minimum of the mean squared error (MSE) of state estimate x̂kjk : Ef(xk ¡ x̂kjk )(xk ¡ x̂kjk )T g ¸ J¡1 k : (42) In the TDOA scenario, the CRLB is independent of the actual measurement set Z k and depends only on the emitter/UAV geometry, sampling interval ±tk , measurement noise variance ¾w2 , and process noise covariance Qk . In the absence of process noise, the CRLB may be found using the covariance recursions of the EKF, where the Jacobian (29) is evaluated at the true state, xk , rather than the predicted state [24]. The CRLB computation in this case corresponds to recursions of the FIM given by T ¡1 Jk = [F¡1 k¡1 ] Jk¡1 Fk¡1 + 1 T H̃ H̃ ¾w2 k k (43) with Jacobian ! à ¯ T T (¡ xk ¡ s(2) (¡ xk ¡ s(1) @hk (x) ¯¯ k ) k ) ¡: ¡ = H̃k = @x ¯x=xk k¡ xk ¡ s(1) k¡ xk ¡ s(2) k k k k (44) The recursion (43) is initialized by J1 = §1¡1 , where two types of initial covariances §1 are considered. One uses the initial combined covariance formed by the GMM technique described in Section III. (Note that this does not require measurement z1 ; it is based on true h(x1 ) and ¾w ). The second uses the initial covariance corresponding to the initial component that is closest to the true emitter location. V. COMPARISON RESULTS A. Scenario Description (i) (j)(Â(i) Wkjk¡1 kjk¡1 (j) ¡ x̂kjk¡1 (j)) £ (hk (Â(i) kjk¡1 (j)) ¡ zkjk¡1 (j)) Pzz (j) = (37) and the probabilities are updated according to (31). The number of components in the UKFB never increases and there is no need for component management. The output of the UKFB may be the highest probability component or the combined mean and covariance of the Gaussian mixture. (39) (i) 2 Wkjk¡1 (j)(hk (Â(i) kjk¡1 (j)) ¡ zkjk¡1 (j)) : (40) The likelihood at each UKF is ¤k (j) = N fzk ; zjkjk¡1 (j), Sk (j)g (41) A target is initially located at (0, ¡1700) m and moves in a northerly direction at 15 m/s. Its radar initially transmits in a northerly direction and scans clockwise with a constant scan time ±tk = 2:5 s, transmitting 666 pulses per second. The radar beamwidth is 15 deg including the antenna sidelobes. Fig. 1 shows a scenario in which UAVs travel in a southwesterly direction from their initial position with speed 30 m/s. After 75 s both UAVs manoeuvre maintaining the same speed turning to the southeasterly direction. OKELLO ET AL.: COMPARISON OF RECURSIVE ALGORITHMS FOR EMITTER LOCALISATION 1727 Fig. 1. Simulation scenario. Initial positions shown with triangles. Fig. 3. A GMM approximation of first RDOA measurement. Asterisk at (0, ¡1700)m is true radar position. example of the GMM approximation for the test scenario is shown in Fig. 3. The GMM-ITSF uses converted measurements with Gk = 40 components, and the number of components is limited to 100. All three algorithms prune track components with probabilities below 0.001. Algorithm divergence was defined using the normalized estimation error squared (NEES), defined by NEES(k) = (xk ¡ x̂kjk )T P¡1 kjk (xk ¡ x̂kjk ): Fig. 2. True RDOA and a realization of RDOA measurements. The first signal at each scan that is received by both UAVs is used to calculate the TDOA. Since the separation between the UAVs is known to be 300 m, then the maximum possible magnitude of time difference between associated TOAs at each UAV is 300=c s (c ¼ 3 £ 108 m/s). At each radar scan, multiple transmissions may be received by both UAVs. However, only the first transmission at which both UAVs are within the beamwidth of the signal is used. A total of 60 TDOA measurements are received over a period of approximately 150 s. Fig. 2 shows theptrue RDOA and the measured RDOAs (with ¾w = 50 m). UAV2 is initially closer to the target, and from around 75 s UAV1 is closer to the emitter and the RDOA changes sign. The search region is limited to [¡4000, 8000] m along the x-axis and [¡7000, 5000] m along the y-axis. B. Numerical Results The first measurement is approximated by a Gaussian mixture of G1 = 39 measurement components (Section III), used to initialize all three algorithms, GMM-ITSF, EKFB, and UKFB. An 1728 If the filter is consistent, the NEES should be Â2 -distributed with Nx degrees of freedom, where Nx is the dimension of the emitter state. A track was deemed to have diverged if the NEES was outside the 99% confidence limits, given by [0:21, 14:9], for 10 consecutive iterations updates beyond the 30th. The 30th update was chosen as the CRLB, shown in Figs. 4 and 5, had approximately converged by this iteration. Tracks were also deemed to have diverged if the relative probabilities of all track components were zero. Table I shows the number of tracks that diverged during 1000 Monte Carlo runs. The GMM-ITSF suffered the smallest number of divergent tracks. Less than 2% of the 1000 Monte Carlo runs were divergent for any tracker. The CPU time required to complete the 1000 Monte Carlo simulations for each algorithm is summarized in Table II. The GMM-ITSF required significantly more computational effort than the other two algorithms, however even GMM-ITSF simulations were finalized in only 0.8% of the total simulated time. All presented algorithms may be comfortably used in a real time application. The rms estimation error has been calculated and p compared with CRLB for position and velocity. The diverged tracks were not included in these statistics. Two different types of outputs are compared against the CRLBs. Fig. 4 shows the rms error for IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011 Fig. 4. RMS errors of highest probability components. Fig. 5. RMS errors of combined estimate. TABLE I Number of Diverged Tracks (out of 1000) EKFB UKFB GMM-ITSF Highest Probability Track Components Combined Tracks 18 16 4 10 12 3 TABLE II CPU Time (s) Required for Each Tracker from 1000 Monte Carlo Runs EKFB UKFB GMM-ITSF 140.22 349.19 1177.25 the highest probability track component of each algorithm, while Fig. 5 displays the rms errors for the combined estimates via (23) and (24). Two types of CRLBs (closest component and combined initial covariance) are also presented. Overall all three algorithms demonstrate similar error performance. In the early stages of tracking, the GMM-ITSF performs slightly better than the other two algorithms, due to the aggressive mixture reduction technique. Towards the end of the simulation, after the CRLB for position has converged (around the 30th scan), the UKFB displays marginally better than EKFB and GMM-ITSF followed closely by the EKFB and GMM-ITSF. VI. CONCLUSION This paper considers tracking of an emitter using a time sequence of TDOA measurements, obtained by a single pair of UAVs. The CRLB is calculated, and three algorithms of various complexity are compared with the CRLB. The algorithms are GMM-ITSF, EKFB, and UKFB. The rms estimation performance of all algorithms were close to the CRLB, with the UKFB having a slight edge. On the other hand, GMM-ITSF had a smaller number of diverged tracks. As all algorithms presented are compatible with the real time requirements, the priorities of individual applications should determine the choice. OKELLO ET AL.: COMPARISON OF RECURSIVE ALGORITHMS FOR EMITTER LOCALISATION 1729 [13] REFERENCES [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] 1730 Finn, A., Brown, K., and Lindsay, T. Miniature UAVs and future electronic warfare. In Proceedings of the Land Warfare Conference 2002, Brisbane, Australia, Oct. 2002, 93—106. Fang, B. T. Simple solutions for hyperbolic and related position fixes. IEEE Transactions on Aerospace and Electronic Systems, 26, 5 (Sept. 1990), 748—753. Torrieri, D. J. Statistical theory of passive location systems. IEEE Transactions on Aerospace and Electronic Systems, 20, 2 (Mar. 1984), 183—198. Drake, S. and Dogancay, K. Geolocation by time difference of arrival using hyperbolic asymptotes. In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, May 2004, 361—364. Hashemi-Sakhtsari, A. and Dogancay, K. Recursive least squares solution to source tracking using time difference of arrival. In Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing, May 2004, 385—388. Wann, C. and Chen, Y. Mobile location tracking with velocity estimation. In Proceedings of the IEEE 5th International Conference on Intelligent Transportation Systems, Sept. 2002, 566—571. Yu, K., et al. UWB location and tracking for wireless embedded networks. Signal Processing, 86 (2006), 2153—2171. Okello, N. and Mu²sicki, D. Emitter geolocation with two UAVs. In Proceedings of Information, Decision and Control Conference, Adelaide, Australia, Feb. 2007. Dogancay, K. and Hashemi-Sakhtsari, A. Target tracking by time difference of arrival using recursive smoothing. Signal Processing, 85 (2005), 667—679. Ho, K. and Chan, Y. Geolocation of an unknown altitude object from TDOA and FDOA measurements. IEEE Transactions on Aerospace and Electronic Systems, 23, 3 (July 1997), 770—783. Ho, K. C. and Chan, Y. C. Solution and performance analysis of geolocation by TDOA. IEEE Transactions on Aerospace and Electronic Systems, 29, 4 (Oct. 1993), 1311—1322. Savage, C., Cramer, R., and Schmitt, H. TDOA geolocation with the unscented Kalman filter. In Proceedings of 2006 IEEE International Conference on Networking, Sensing and Control, 2006, 602—606. [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] Mikhalev, A. and Ormondroyd, R. Passive emitter geolocation using agent-based data fusion of AOA, TDOA and FDOA measurements. In Proceedings of 10th International Conference on Information Fusion, Quebec, Canada, Oct. 2007. Mikhalev, A. Comparison of Hough transform and particle filter methods of emitter geolocation using fusion of TDOA data. In Proceedings of 4th Workshop on Positioning, Navigation and Communication, Hannover, Germany, 2007, 121—127. Sathyan, T., Sinha, A., and Kirubarajan, T. Passive geolocation and tracking of an unknown number of emitters. IEEE Transactions on Aerospace and Electronic Systems, 42, 2 (Apr. 2006), 740—750. Okello, N. and Mu²sicki, D. Measurement association for emitter geolocation with two UAVs. In Proceedings of 10th International Conference on Information Fusion, Quebec, Canada, July 2007. Mu²sicki, D., Kaune, R., and Koch, W. Mobile emitter geolocation and tracking using TDOA and FDOA measurements. IEEE Transactions on Signal Processing, 58, 3 (Mar. 2010), 1863—1874. Fletcher, F., Ristic, B., and Mu²sicki, D. Recursive estimation of emitter location using TDOA measurements from two UAVs. In Proceedings of 10th International Conference on Information Fusion, Quebec, Canada, July 2007. Alspach, D. and Sorenson, H. Nonlinear Bayesian estimation using Gaussian sum approximation. IEEE Transactions on Automatic Control, AC-17, 4 (Aug. 1972), 439—448. Jazwinski, A. H. Stochastic Processes and Filtering Theory. New York: Academic Press, 1970. Bar-Shalom, Y., Li, X. R., and Kirubarajan, T. Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software. Hoboken, NJ: Wiley, 2001. Julier, S. and Uhlman, J. Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 92, 3 (Mar. 2004), 401—422. Trees, H. L. V. Detection, Estimation and Modulation Theory, Part I. Hoboken, NJ: Wiley, 1968. Taylor, J. H. The Cramer-Rao estimation error lower bound computation for deterministic nonlinear systems. IEEE Transactions on Automatic Control, AC-24, 2 (Apr. 1979), 343—344. IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011 Nickens Okello (S’87–M’94) received the B.S., M.S., and Ph.D. degrees all in electrical engineering from Makerere University, Kampala, in 1978, Southern Illinois University at Carbondale, in 1988, and The Pennsylvania State University, University Park, in 1994, respectively. From 1995 to 1999 he was a Research Fellow at the Cooperative Research Centre for Sensor Signal and Information Processing (CSSIP) at Technology Park, Adelaide, Australia. He was a founding member of CSSIP’s Data and Information Fusion Program. During this period he developed advanced track and sensor registration algorithms, and delivered shorts courses to graduate students and industry R&D participants. From 2000 to 2004 he was with CSSIP’s Melbourne node and from 2005 to 2007 he was with Melbourne Systems Laboratory both within the Department of Electrical and Electronic Engineering, University of Melbourne. During these periods he worked on industry research contracts on registration, air picture compilation, situation assessment, geolocation, and network-centric warfare problems. He is now a senior researcher with National ICT Australia (NICTA) where he works on the application of data and information fusion on water information networks. His research interests include, multisensor data and information fusion, track and sensor registration, emitter localisation using UAVS, assessment of threat using Bayesian networks, time series analysis, and multispectral and hyperspectral signal processing. He has published numerous journal and conference papers on these topics. Fiona Fletcher received the B.Sc. with honours in applied mathematics from the University of Adelaide, Australia, in 1999. She received a Ph.D. in electrical and electronic engineering from the University of Melbourne, Australia in 2005. She worked in the Maritime Operations Division of the Defence Science and Technology Organisation, Australia from 1999 until 2005 and is currently employed in Intelligence, Surveillance and Reconnaissance Division of the Defence Science and Technology Organisation. Her research interests include passive and active tracking, multi-target tracking and sensor fusion. Darko Mu²sicki (M’90) was born in Belgrade, Serbia, in 1957. He completed the B.E. and M.Sc. degrees in electrical engineering in 1979 and 1984, respectively, at the University of Belgrade where he was a member of the academic staff between 1985 and 1988. He moved to Australia in 1988 where he completed a Ph.D. in 1994 at the University of Newcastle, NSW. Between 2002 and 2007 he was with the University of Melbourne, Australia. He is now with the Department of Electronic Systems Engineering, Hanyang University, Kyunggido, Korea. Since 1979 he worked in the area of radar system design and simulations, as well as fast implementation of digital signal processing algorithms, and from 1984 Darko Mu²sicki was involved in target tracking research and design, sensor information fusion and nonlinear estimation, as well as resource allocation and wireless sensor networks. Dr. Mu²sicki is member of Board of Directors of ISIF (International Society for Information Fusion), and has served as ISIF president during 2008. OKELLO ET AL.: COMPARISON OF RECURSIVE ALGORITHMS FOR EMITTER LOCALISATION 1731 Branko Ristic received all his degrees in electrical engineering: Ph.D. from Queensland University of Technology, Australia, in 1995, M.Sc. from Belgrade University, Serbia, in 1991, and B.Eng. from The University of Novi Sad, Serbia, in 1984. He has been involved in R&D related to signal processing, estimation, tracking and data fusion for over 25 years. He began his career as a research engineer in 1984 at Vinca Institute, Belgrade, until he emigrated to Australia in 1989. From 1989 to 1994 he was with The University of Queensland and the Queensland University of Technology. During 1995 he was a DSP engineer in GEC Marconi Systems in Sydney. Since 1996 he has been with Defence Science and Technology Organisation (DSTO). His role in DSTO has been to carry out research, develop new capabilities, and provide technical advice to the Australian Defence Organisation on topics of target tracking and data fusion. He is currently a DSTO Fellow. During 2003/2004 he spent a year in IRIDIA (Université libre de Bruxelles, Belgium) on a study leave. Since 2006 he is an Honorary Fellow of The University of Melbourne. Dr. Ristic coauthored two books, over 150 technical papers and presented several invited talks, short courses, and tutorials. He served on technical committees of several international conferences, and was the Chair of the Fourth Australian Data Fusion Symposium in 2007. He is best know as a coauthor of the book Beyond the Kalman Filter: Particle Filters for Tracking Applications (Artech House, 2004). Three times he has received “the best-paper award” at international conferences (Fusion 05 and DICTA 05 and 09). 1732 IEEE TRANSACTIONS ON AEROSPACE AND ELECTRONIC SYSTEMS VOL. 47, NO. 3 JULY 2011