POSITIONING USING TIME-DIFFERENCE OF ARRIVAL MEASUREMENTS

Fredrik Gustafsson and Fredrik Gunnarsson

Department of Electrical Engineering

I,['
Linkoping University, SE-581 83 Linkoping, Sweden
Emdl: fredrik@isy.liu.se,fred@isy.liu.se

ABSTRACT Constant TDOA "sin0 WO receivers

The problem of position estimation from Time Difference

Of Arrival (TDOA) measurements occurs in a range of ap-
plications from wireless communication networks to elec-
tronic warfare positioning. Correlation analysis of the trans-
mitted signal to two receivers gives rise to one hyperbolic
-1
function. With more than two receivers, we can compute
more hyperbolic functions, which ideally intersect in one
unique point. With TDOA measurement uncertainty, we
face a non-linear estimation problem. We here suggest and -2
compare both a Monte Carlo based method for position-
ing and a gradient search algorithm using a non-linear least -3
-3 -2 -1 0 1 2 3
squares framework. The former has the feature to be easily
extended to a dynamic framework where a motion model Fig. 1. The hyperbolic function representing constant
of the transmitter is included. A small simulation study is TDOA for three different TDOA's (0.4, 0.6 and 0.9 scale
presented. units, respectively).

1. INTRODUCTION

Figure 1 illustrates how two cooperating receivers can cal-

. Electronic warfare, where the problem is to accurately
locate enemy transmitters to be able to make appro-
culate a path difference from the time-difference of &Val, priate countermeasures. The TDOA approach may
and how this path difference corresponds to a hyperbolic here offer higher accuracy than traditional triangula-
function [I]. We point out two important applications where tion approaches, where the angular measurements po-
this problem occurs: tentially have larger inaccuracy than TDOA measure-
ments.
Mobile terminal positioning. The available measure-
ments are either network-assisted or mobile-assisted; The TDOA measurement is computed as follows:
using up-linkor down-link information. For an overview,
see [2, 31. The current 'yellow page' services are I . The sender transmits a signal s ( t ) which is delayed
T$ to receiver i according to distance to each receiver.
based on Received Signal Strength (RSS) of signals
with known powers and a course Angle measurement The signal can either be a pilot from a mobile (up-
from the sector antenna [4, 51. This gives a quite link), where the mobile's absolute time is unknown,
course estimate. It is well-known [6,4, 11 that a much or it can be unknown, as is the case in electronic war-
higher accuracy rather insensitive to fading is based fare. In either case, ricannot be computed.
on time of arrival (TOA). In future system, only the
2. Correlation analysis provides a time delay ri - r3cor-
time difference of arrival (TDOA) or enhanced ob-
responding to the path difference to receivers i and j .
served time difference (E-OTD) measurements may
be possible to compute. Another source of informa- In [8],a general framework covering all these kind of mea-
tion for tracking moving objects is map information, surements (TOA, TDOA, E-TDOA, RSS, Angle, Map) is
see 171. given, provided that a motion model for the transmitter is
This work was supported by Vinnova's competence center ISlS available. That is, a dynamic framework is assumed, and

0-7803-7663-3/03/$17.00 0 2 0 0 3 IEEE VI - 553 ICASSP 2003

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 30,2025 at 13:52:12 UTC from IEEE Xplore. Restrictions apply.
 the particle filter is suggested. We here study the static case,
where we have only one TDOA measurement for each pair
of receivers. We will compare a particle filter based esti-
mator with a least mean square (LMS) algorithm in a non-
linear least squares framework.

2. TDOA MEASUREMENTS

The received signals are

y i ( t ) = ais(t - 7 1 ) + ei(t), i = 1 , 2 , . . .n (1)

where the receiver i is located at x i , y%and the transmitter

is in x,y. which is unknown.
With known reference s ( t ) and perfect synchronization, Fig. 2. Same as Figure 1, but with measurement uncenainty
we can directly estimate T~ (TOA), and estimate (x,y) using in Ad.
a non-linear least squares framework, similar to GPS.
With unknown reference, the simplest idea is to com- For a general receiver position, we simply translate the
pare the received signals pairwise. Assume a correlation hyperbolic function ( 5 ) in local coordinates (x,y) to global
function that painvise computes an estimate of coordinates ( X ,Y ) using
Adi,j = V(T~ - ~ j ) , 1 5 i < j 5 n. (2)
where zi is the speed of sound, light or water vibrations.
Here, n is the number of receivers and ( i ,j ) is an enumera- + +
where X O = ( X i X j ) / 2 , YO= (K y 3 ) / 2 locates the
tion of all K pair of receivers, where center point of the receiver pair. The hyperbolic function in
global coordinates is thus given by A d i , j = h(z,y, D ) in
(3) (3,with
K=(;j
Each Adi,j corresponds to positions (x, y) along a hyper- D= d(X - y3)2 + (Xi - X .3 ) 2 (84
bola.
Assume first that the receivers are both located at the (8b)
x-axis at x = D / 2 and x.= -012, respectively. The hy-
perbolic function can then be expressed as (8,)

d2 = & + (x + D / 2 ) 2 , We have now a functional form suitable for representing the

di = - J y2 + (Z- D/Z)’, (4b)
measurement uncertainty in TDOA, which implies an un-
certain hyperbolic area rather than a line. This is illustrated
Ad = dz - d i = h(x, y, D ) (4c) in Figure 2. The farther away along the asymptotes, the
= Jy2 + (x + 012)’ - Jy2 + (x - D ) / 2 ) 2 . (4d) larger absolute uncertainty in position.

After some simplifications, this equation can be rewritten in

3. THE NON-LINEAR LEAST SQUARES
a more compact form as
PROBLEM

The general problem is to solve the (possibly over-determined)

non-linear system of K equations
The solution to this equation has asymptotes along the lines
Ad,,, = h ( X , Y ; X i , K , X j , Y , ) , l 5 i < j < n (9)
b
y = *-x =
a
*
.x-/ D2/4 - Ad2/4
Ad2/4
(6) for the sender position ( X ,Y ) ,given the receiver positions
( X i ,Y,).Now, the non-linear least squares estimate of ( X ,Y )
is given by
which defines the angle of amval for far-away transmitters.
Figure 1 illustrates the hyperbolic function in the local co-
ordinate system (x, y).

v1- 554

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 30,2025 at 13:52:12 UTC from IEEE Xplore. Restrictions apply.
 To simplify notation, we use P = (X,Y ) for the position. Algorithm 1 (Stochastic gradient algorithm)
Then, rewrite the minimization problem in vector notation
using a weighted least sqauares criterion p ( m + l ) = p(m)- J m ) h b ( p ( m ) ) ( ~ d- h ( p ( m ) ) )(13)

P = = argmin(Ad - / L ( P ) ) ~ R - ' ( A-~h ( P ) ) (IO) A good step size can be computed by bisection techniques,
P line search 01; as in the simulations, as the normalized LMS
where A d = . .Ad,t-l,n)Tand R = C o v ( A d ) is step size
the covariance matrix for the TDOA measurements. The so-
lution defines the minimum variance estimate. For a Gaus-
sian assumption on the TDOA noise, this coincides with the
maximum likelihood estimate. The particle filter is a static version of the well-known
+
Using the assumption that A d = h(Po) e, where Po SIR algorithm [9,101.
is the true position and the TDOA noise has Cov(e) = R, Algorithm 2 (Static particle filter) .
a first order Taylor expansion around the true value gives
h ( P ) zz h(Po)+hb(Po)(P-Po). Theleastsquarestheory
now gives 1. Randomize N 'particles' (herepossiblepositions)P'.

Cov(P) = (h',(Po))+R((hb(Po))+)T,
(1 1) 2. Choose jittering constants C Rand CQand let the p o -
sition random walk covariance Q = C Q / k 2and jit-
provided that P is sufficiently close to the true position (high +
tering measurement noise R = R G R / ~The ~ .idea
SNR). For Gaussian noise e , this expression also defines the with jittering noise is to explore a smaller andsmaller
Cramer-Rao lower hound. That is, no estimator can per- neighborhoodmore and more accurately.
form better than this bound, given that we have find a small
enough neighborhood of the m e position. From (1 I), we 3. Iterate f o r k = 1 , 2 , .. . until P ( k ) has converged.
can obtain guidelines for what a favourable Po is, or for (a) Compute the particle weights wi using the like-
electronic warfare how to place the receivers in the best pos- lihood
sible way.
touz= cxp ((Ad - h ( P i ) ) T R - ' ( A d - h ( P i ) ) )
4. ALGORITHMS
andnormalize w i= u P / ( C w i ) .
Three different approaches have been compared: ( b ) Compute the estimate @(k) = E, wiPz.
1. Compute the intersection point of each pair of hyper- (cj Resample with replacement the particles, where
bolic functions. There are the probability to pick one particle is p m p o r -
tional to its weight. After the resampling, the
weigths are reset wi= 1 / N .
(dj Spreadout the particles as Pi= Pi + w,where
w E N(0,Q).
pair of hyperbolic functions. The position estimate
can then be the (weighted) average of these points. The resampling step is the key to get a working algorithm.
Since each pair can have no, one or two intersections, In the standard particle filter, k denotes time and there is a
time update step where the particles are moved according to
the logic to find the correct one is non-trivial. Further,
it is a bit complicated to find the correct weighting. a velocity measurement and a movement noise w , otherwise
the algorithms are quite similar. Compare to the particle
2. Applying the stochastic gradient algorithm to the non- filters in [7, 81.
linear least squares problem.
3. Numerical approximation of the non-linear least squares 5. SIMULATIONS
problem using Monte Carlo based techniques. This is
Figure 3(a) illustrates the test scenario, with four receivers
refered to as the particle filter (PF), which is a static
computing in total six TDOA measurements.
version of the general approach described in [8].
Figure 3(b) shows what happens to the hyperbolic func-
The first approach has shown to give inferior results than the tions when Gaussian measurement noise (standard devia-
other two. The latter two ones are described below. tion of 0.1 scale units) is added to the TDOA measurements.
The normalized gradient algorithm can be written as fol- That is, there is no clear cut intersection point. To under-
lows: stand the non-linear least squares criterion, the level curves

VI - 555

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 30,2025 at 13:52:12 UTC from IEEE Xplore. Restrictions apply.
 Fig. 3. (a) Test scenario: four receivers are placed in a square, and six resulting hyperbolic functions from noise-free TD0A:s
intersect at the transmitter position. Also shown is the particle cloud and the resulting position estimate using Algorithm 1 . (b)
Same as (a), but the hyperbolic functions are computed from six different noisy TDOA vectors. This illustrates that in general
there is no unique intersection of all six lines. (c) Contour plot of non-linear least squares criterion &(Adi,) - h ( X ,Y ) ) ’ .
In this sceario, there is no local minima and a gradient algorithm will converge from any initialization. (d) Gradient search
using a normalized least mean square method on Algorithm 2 (compare the path to the contour plat in (c)).

of (10) are plotted in Figure 2(c). This plot explains how [2] C . Drane, M. Macnaughtan, and C. Scott, “Positioning GSM
the weights of the particles are computed, and in which di- telephones,” IEEE Communications Magazine, vol. 36, no.
rection the gradient in LMS points. 4, 1998.
Figure 3(a) shows the particle cloud after a few itera- [31 M. Silventoinen and T. Rantalainen, “Mobile station locating
tions. As might be expected from the level curves, this cloud in GSM:’ in Proc. IEEE Wireless Communication System
is often found beyond the true position, as in this example. Symposium, Nov 1995.
However, in most cases the true position is found. Figure [4] M.A. Spinto, “Further results on GSM mobile station loca-
3(d), finally, shows the leaming path of LMS. For this sce- tion:’ IEE Electronics Lerrm, vol. 35, no. 22, 1999.
nario which lacks local minima, LMS is to prefer. [SI B. Mark and Z. Zaidi, “Robust mobility tracking for cellu-
lar networks,” in Proc. IEEE Intemarional Communicarions
Conference, New York, NY, USA, 2002.
6. CONCLUSIONS
[6] S. Fischer, H. Koorapaty, E. Larsson, and A. Kangas, “Sys-
tem performance evaluation of mobile positioning methods,”
Two algorithms have been suggested for finding the position in Proc. IEEE Vehicular Technology Conference, Houston,
of a transmitter, given TDOA measurements computed from TX, USA, May 1999.
the received signal for at least three receivers. A non-linear
[7] E Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell,
least squares approach was advocated, enabling local anal- J. Jansson, R. Karlsson. and P-J. Nordlund, “Particle filters
ysis yielding a position covariance and a CramBr-Rao lower for positioning, navigation and tracking:’ IEEE Trnnsacrions
bound. A simulation study illustrated the TDOA problem on Signal Pmcessing, vol. 50, no. 2, February 2002.
in general and the performance of the two suggested algo- [X] P-J. Nordlund, E Gunnarsson, and E Gustafsson, “Particle
rithms. On-going work aims at investigating the practical filters for positioning in wireless networks,” in h i r e d ro
performance for acoustic communication, how to use the EUSIPCO, Toulouse, France, September 2002.
CramBr-Rao bound to determine receiver locations and use
[91 N.1. Gordon, D.J. Salmond, and A.F.M. Smith, ” A novel
Monte-Carlo simulations to examine the performance of the
approach to nonlinearlnon-Gaussian Bayesian state estima-
two algorithms. tion:’ in IEE Proceedings on Radar and Signal Processing,
1993, vol. 140, pp. 107-1 13.
7. REFERENCES [IO] A. Doucet, N. de Freitas, and N. Gordon, Eds., Sequential
Monte Carlo Methods in Practice, Springer Verlag, 2001.
[I] M.A. Spinto and A.G. Mattioli, “On the hyperbolic position-
ing of GSM mobile stations,” in Pmc. lnrernational Sympo-
sium on Signals, Systems and Electronics, Sept 1998.

VI - 556

Authorized licensed use limited to: INDIAN INSTITUTE OF TECHNOLOGY DELHI. Downloaded on June 30,2025 at 13:52:12 UTC from IEEE Xplore. Restrictions apply.