See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/3231496

A passive localization algorithm and its accuracy analysis

Article  in  IEEE Journal of Oceanic Engineering · February 1987

DOI: 10.1109/JOE.1987.1145216 · Source: IEEE Xplore

CITATIONS READS

329 1,114

All content following this page was uploaded by Benjamin Friedlander on 28 July 2015.

The user has requested enhancement of the downloaded file.

 . . . . . . .. ~

234 JEEE JOURNAL

OF OCEANIC ENGINEEiUNG, VOL. OE-12, NO. 1, JANUARY 1987

A Passive Localization Algorithm and Its Accuracy

Analysis
BENJAMIN FRIEDLANDER, FELLOW,

(Invited Paper)

Abstract-Theproblem of estimatingsolrrce location from noisy

measurements of range differences W ’ s ) is considered. A localization
technique b v d on solving a set of h e a r equations is presented and its
accuracy properties are analyzed. An optimal weighting matrix for the
leastsquaresestimator is derived.Theanalyticalexpressions for the
varianceandbias of theestimatorarevalidated by Monte-Carlo
simulation. Theproblem of estimating source velocity givenmeasure- -+K source

ments of range differences and range-rate differences is briefly consid-

ered, and a lidear equation technique is derived. sensor 1 2s
I. INTRODUCTION
%sensor 3
A common method for localizing an acoustic or
electromagnetic source is based on measuring the range
differences to several points whose locations are known. In
various navigation systems (e.g., Loran, Decca, Omega) the r13
/ r12

measurement consists of observing differences in the time of Fig. 1. Hyperbolic lines of position.
arrival of signals from sources at known locations to a receiver
at the unknown location. In the passive sonat problem, signals
travel from a source whose location is to beestimated to
sensors whose positions are known.
Sensor
The conventional approach to estimating source location is
based on finding the hyperbolic lines of position (LOP’S);see
[9]-[ll]. Each range (or time) difference determines a
hyperbola, and the point at which these hyperbolas intersect
is the estimated source location (see Fig. 1).
The hyperbolic LOP approach has several serious draw-
backs:thecomputation of the intersection location is quite
cumbersome; these solutions do not appear to be easily
extended to other situations such as calculating source velocity
Fig. 2. Location on the conic axis.
from range differences and their rates of change (cf. Section
V). Furthermore, the complexity of the solutions makes error
analysis very difficult. The method proposed in [11 computed the intersection point
An alternative approach, which circumvents many of these from a set of(:) P (A?/(N - 3)!3!) linear equations, Nbeing
difficulties, was proposed in 113. This approach is based onthe the number of sensors. This approach alleviates some of the
ideathat three sensors with their set of range differences difficulties of the hyperbolic LOP methd. However, it
determine a straight line of position. In fact, this line is the introduces some inconsistencies due to its use of redmdant
major axis of a general conic whichpasses through the information: this problem formulation leads to a set of (;’)
equations while, In fact, there are only N available measure-
sensors. When more than three sensors are available, several
straight LOP’s are generated and their intersection provides an ments (see [2] for a more detailed discussion).
estimate of the source location (see Fig. 2). A different formulation is proposed in [2] and [3] based on
the straight LOP approach, which leads to a set of N
equations. This set of equations contains all the relevant sensor
Manuscript received May 5 , 1986; revised September 5 , 1986. This work
was supported in part by theOffice of NavalResearchunderContract data without redundancy. However, due to a noalinear
N00014-87-C-0147. relationship between certain variables in these equations, the
The author is with Saxpy Computer Corp., Sunnyvale, CA 94086. location can no longer be computed by a linear estimation
IEEE Log Number 8714343.
1 We are considering here the two-dimensional case, i.e., the source and
procedure. An iterative gradient search procedure is proposed
sensors are located in a plane. in [3] to compute the the source location.

0364-9059/87/0100-0234$01.00 0 1987 IEEE

FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM 235

Subtracting one of the equations in [2] and [3] from the (N The various quantities defined above are depicted in Fig. 3.
- 1) remaining equations leads to a set of (N- 1) equations The localization problem is to estimate x, givennoise
containing a single nonlinear term. It was recently noted in[4] measurements of the RD’s { rij,i, j = 1, . ,N)provided by
that this nonlinearity can be eliminated by a two-step least N sensors. Note that there are (f)= (N(N- 1)/2) distinct
squares procedure to yield a closed-form expression for the RD’s. However, all of these RD’s can be completely
source location. This noniterative technique appears to be very determined from only (N- 1) RD’s (e.g., { ril, i = 2, * -, -
attractive both in its modest computational requirements and N))in the noiseless case.
its performance. Using the notation above, we can derive the basic equations
In this paper we present a localization technique based on in a straightforward manner. Note that
the solutionof a set of either (N - 1) or (N - 2) linear
equations. While the derivation and resulting solution appear R ~ = ( r u + R j s ) 2 = ~ + 2 R j s r u + R ~ s (2)
different from the onein [4], it isshownthat the two are and also that
mathematically equivalent. However, themethod presented
here is somewhat simpler and reveals more clearly the (3)X ~ + R
R ~ = I I X ~ - X , ~ ~ ~ = R ~ ~ - ~ X ~
underlying structure of the solution.
The main contribution of this paper isin providing an Subtracting (3) from (2) and rearranging terms, we get
accuracy analysis of the proposed solution method. In Section 2 x , ? x s = R ~ - r ~ . - 2 R j S r u + R & -JSR 2 (4)
III we derive approximate expressions for the bias and
variance of the location estimates. The variance formula is Setting i = j in (4) gives
used to derive an optimalweighting matrix for the least
squares location estimator. In Section VI we verify these +
2xj’X. = R;o Ri0 - R;, (5)
formulas by Monte-Carlo simulations. where we used the fact that rjj = 0. Subtracting (5) from (4)
The methodspresentedin [1]-[4] andinSection II all yields
require that the number of sensors N be greater than (n + 1)
where n is the dimension of the space in which thesensors and 2 ( ~ ~ - x ~ > ~ x ~ = ( R ~ - R ~ , J - r ~ ~(6)
-2R~,
source are located (in other words, N > 3 in the two-
dimensional case and N > 4 in thethree-dimensionalcase). In which can be written in matrix form as
some practical situations no more than N = n + 1 sensors are
(74
available. In Section IV we derive explicit expressions for the
source Iocation for this case. These expressions are somewhat
morecomplicatedthan the one for N > n + 1, and their
accuracy is not analyzed here.
In some situations, it is possible to measure both the range
differences and the range-rate differences (e.g., by measuring , (N-l)x3
Doppler shifts of spectral lines). Using these measurements it
is possible to estimate source velocity as well as source
location. A linear equation technique for estimating source
velocity is presented in Section V. In Section VII we discuss
some additional applications of the results derived in this
paper.

11. LOCALIZATION FROM RANGE DIFFERENCES.

In this section we state and solve the localization problem 1 R;-lo -Rjo -rj-.lj
for the three-dimensional case. The discussion here and in the ’j= 5 R;,lo -R;o - r j i l j , ( N - l ) x 1 (7c)
following sections will be in terms of range difference
measurements. However, the same results apply to localiza-
tion from delay difference measurements. Delaycanbe
converted into range simply by multiplying it by the speed of
propagation.
Let xi = [xi,y i , zi] denote the (x,y , z ) coordinates of the
ith sensor and x, = [xs,y s , z,] denote the source coordinates.
The distance between thesource and the ith sensor is Rk & llxi
- x, 11. The distance from the origin to the ith sensor is Ri, 9
Ilxill, and the distance to the source is R, = Ilx,ll. Let rij be
the range difference (RD)for sensors i and j where
The reference sensor j can be chosen to be any one of the N
sensors. We will comment on this further in Section VI.
-- --.
236 IEEE JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12, NO. 1, JANUARY 1987

z least squares equation solver on the linear set of (1 l), see for

t example [5] and [6]. A closed-form solution which is useful

for some subsequent derivations is given by
Xi sensor i
X, = (STMTMj S,) - STMTMjpj.’ (12)

il It should be noted that Mi is a singular matrix of rank N - 2.

For MjS’ to be nonsingular it is necessary (althoughnot
sufficient) that the number of rows of Sj be larger than the
Y number of columns, i.e., N - 1 > n (where n is the
dimension of x,). Thus, for a unique solution of (1 1) to exist,
wemust have N > n + 1. The case N = n + 1 will be
discussed later.
The approach taken in [4] is based on a two-step procedure.
5.S First it is assumed that Rjs is known, in which case a least
Fig. 3. Sensor-source geometry. squares solution for x, is given by
X,=(STS’)-~S~T(,U~-Rjspj). (13)
Note that given the measurements { rij, i = 1, --
N},the
vectors p j , p j are known; so is the matrix Si,which depends Next Rjs is chosen to minimize IISjxs - ( p j - Rsjpj)ll or
only on the sensor locations (assumed to be hown). The equivalently 11 Pj(pj - Rsjpj)ll where
unknown quantities in (7a)-(7d) are the location vectors x, and pj=r-sj(s;sj)-lsj. (14)
the distance Rj, between the source and the reference sensor,
which is a function of (the unknown) x,. It is the presence of It is straightforward to show that the minimizing Rsj is given
Rjs that complicates the solution of (7). by
Our approach isbased on the ideaof eliminating the
“nuisance parameter” Rjs by premultiplying (7a) by a matrix
M which has pj in its null-space, i.e., Mpj = 0. The matrix
My’= (I- P ) D j (8a) Thus the method presented in [4] computes the location vector
by the least squares solution of
where

rj- Let us denote the solution obtained from ( 11) x i and from (16)
x:. Note that

MjSj~,2=Mjpj=MjSj~B (17)

or
Mjsj(x;- x;) = 0.
Since MjSj is assumed to be nonsingular, it follows that xf =
x:. Establishing this fact by algebraic manipulation of (1 1) and
(16) is quite cumbersome.

clearly has this property for any value of k , since An Alternative Form of (11)
In the following section we w i find it convenient to use a
l
DjPj’l P [l, 1, . e - , 1]T (9) slightly different form of the location estimator. Note that the
and matrix (I - 29 [cf. (8a)l is an (N- 1) X (N- 1) matrix
withrank ( N - 2). Thus the singular valuedecomposition
(I-Zk)l=l-Zkl=1-1=0. (10) (SVD) of this matrix w lihave the form

To simplify the subsequent discussion, we assume that k = 1

=r UkY
and suppress the superscript in My). From the property of Mj (I- Zk) Ukl
discussed above, it follows that TN-2

M j S j X , =Mj pj . (11)

The location vector x, can now be computed by using

any
good = uk diag {T:, - * - y 9i-J V l (18)
FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM 237

where uk, v k are (N- 1) X (N- 2) orthogonal matrices or

-
and { q f, * , are the (N- 2) nonzero singular Values
of ( I - zk).Since ( I - z ~ = I 0, it foUows that VlI = 0.
cov {fs}
=4cFj’ (27a)
where
Therefore, another candidate for the matrix Mj whichwas
used to “annihilate” the second term in (7a) is F j = (STMTM,Sj)- STMTMj (diag { p j } + RjsI). (27b)
A?,= V,‘Dj. (1 9) It should be noted that by performing a first-order Taylor-
series expansion of 2, around x, it follows that Fj is the matrix
Premultiplying (7a) by this matrix leads to a set of (N- 2)
of the first derivatives of 2, with respect to the measurements
linear equations for x,
{ r g } , Le., F j = ( d i S / d p j ) . Thus (27) can be derived by
M.S.X
J J s-
-iff. .
JPJ (20) differentiating (1 1).
whereas (11) consisted of (N - 1) linear equations. In the The derivation of the bias can be done by a second-order
following we will suppress the subscript k in V k . Taylor-series expansion of $. It is straightforward to show
that
111. ACCURACY ANALYSIS OF THE LOCATION ESTIMATE
Inthediscussion so far weassumed the availability of
perfectmeasurements { rij}. In thissectionweanalyzethe
effect of noise on the location estimates. We assume that the and similarly
RD’s are corrupted by additive mutually uncorrelated zero-
mean white Gaussian noise with covariance matrix
C=diag {a’), u 2 = [ a : , a:, -.., u i - , ] ’ . (21)
The second-order derivative matrices of x, and y, are
The measured RD’s are denoted bj and somewhat complicated (cf. [12])andwillnotbe presented
b.-
J-Pj+nj (22) here. Instead we will consider two approximate formulas
which are simpler than the “exact” formula and yet provide
-
where nj X(0, E). The location vector obtained when bj is satisfactory results in many cases.
used instead of pj is denoted is,where 3, = x, + f,, 2, being The first derivation involves approximating the noisy
the error vector. Similarly, we denote Rj, = llxj - $,I1 and ijs
version of Mj (i.e., Mj computed using noisy range differ-
= Rj, + Rjsywhere Rjsjs the error term.
ences) by the true Mi. We then take the expected value of both
Inserting is,bj, and Ris into (7a) we get sides of (24) and premultiply by Mj to get
SjZ, = Sjx, + S,f, o nj + Rj,pj
(
= pj - pi

+-1 nj o nj+Rj,nj) (23)

M,SjE{fs}= -Mj
[I: 3
E - nj
1.
o nj + E { R j , n j } (28)

2 Note that in this case we do not neglect the last two terms of
where denotes element-by-elementmultiplicationof two (24). In fact, it is precisely these t e r n that introduce the bias.
vectors. Alternatively The multiplication. by Mi is necessary to eliminate the term
- 1 E{RjsIpj-
(pi+ Rj,l) o nj+ R j s p j + - nj o nj+ pjsnj] . The first term is the square brackets above is simply u2.To
2 evaluate the second term we insert 6j into (15) to obtain (after
(24) reordering of terms)

To evaluate the covariance matrix of fswemakesome To evaluate R,j we will make asmall error approximation and
approximations. Assuming small error conditions, wewill neglect all terms containing products of noise variables (i.e.,
neglect any terms which contain products of errors. Thus the the last three terms in the numerator and the last term in the

-->
two last terms in (24) are neglected. Next note that denominator). Next we use the standard approximation

MjSj.fs = -Mj(p,+R,,l) o nj (25) 1 --

1 (1 - 2 p?Pjnj (30)
pTPjpj + 2pi’pin; pTPjpj PTPjPj
from which it readily follows that
to get
. . . . . . . . .

238 IEEE JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12, NO. 1, JANUARY 1987

1
Note that A E
K = { MTa.^.} =MTM.
j JPJ J JPJ
MTMj,2
2 J
Pj (L,/, - Rjspj) =Pj( Six,) = 0 (32)
and therefore the second term in (31) reduces to - RjspTnj.It
+MJMjn2+ diag { DjMTMjDJ:I} u2
follows that + diag { DjMTMjDj)u2o p j . (43)
It fbllows that
E{.?,}=H-'STK. (44)
Since 2, = x, + f,, we get
Inserting (33) in (28) we get BIAS=E(fs}=H-'SjK-~,. (45)
Accuracy Optimization
In the discussion so far the location equation (1l ) , or

c
or equivalently (20),was assumedto be solved by astandard least
'
BIAS = E { f,} = - (SiTMjMjSj)- STMTMj - g2+ 6j
1. (35)
squares technique. It is possible, of course, to use instead a
weighted least squares procedure. Let W be a positive definite
(N - 2) X (N - 2) weighting matrix and premultiply (20) by
W'" to yield the weighted version of the location estimator:
Nextwe derive a formula for the biaswhich takes into
account the fact that the Mj used in (11) is corrupted by the Wl/2n;r.X - W1"MjPj
J S- (46a)
measurements noise.
From (22) and (8) it follows that or
x,= (s~MTwMjsj)-ls~M~W~jlLj. (46b)
f i j 9 ZLjj=Z (diag { p j ) + N i ) - ' = Z D j ( I + D j N j ) - ' (36)
The covariance and bias of this estimator are obtained by
where Nj 9 diag { n j } .
replacing Mj with W1I2&%j in (27) and (34) or (35). ,
Using the approximation (I + X ) - ' = 1 - X we get The question naturally arises as to what choice of the
M j G MjCI-DjNj). (37) weighting matrix W wiU lead to the best possible accuracy. To
Next note that answer this question we w li first rewrite the covariance
equation (27) for the weighted case as
aTaj = MTMj -MTMjDjNj - NjDjMTMj
covw{fs} =(ATWA)-1A7WCWA(A*WA)-1 (47a)
+NjDjMTMjDjN,. (38)
where
From (12) it follows that
A iGtjSj (4%)
s;@;.i@jsj]
3, = [ - ls;M;Mjfij. (39)
+
c =Mj (diag { p i } R ~ J ) X (diag { p i } +RjJ)MJ- ( 4 7 ~ )
We will approximate the expected value of the right-hand side
of (39) by Next we show that
E { $ } = [E(S;h?TA?jS'}]-'E{ STi@TA?jfij}. (40) w *= c-' (48)
is the optimal weighting matrix in the sense that
This approximation is a good one when SiTM,TMiSj is a well-
conditioned matrix and the measurement noise is small. From COVw{$}2COVw * {fs}=(ATC-'A)-l (49)
(38) it follows that
for all permissible choices of W. The proof consists of
H 4 E { STMTMjS'} = SjTM,TMjSj verifying that
+ ST diag (DjMTMjDj}C Sj . (41) COV,{f,} -covw * {%}
Next note that = [(A'WA)-'A TW-(A'C-lA)-lATC-l]
,@T$j.^.-MTM.
1
o nj - C[(ATWA)-'ATW-(A*C-'A)-lATC-l]'.
(50)
J P J - J J ~ J 2 MTM.n.
J J J
Since C > 0 the right-hand side of (50) is clearly nonnegative
+MTMjDjNjpj o nj definite and (49) immediately follows.
Notethat the derivation above requires that C and,
+NjDjMTpj o nj+ NjDjMTMjDjNjpj therefore, Mj be nonsingular. This is the reason we used (20)
+terms containing nj rather than (11) in the development above (i.e., Mj is
nonsingular, while Mj is singular).
+terms containing 3rd and 4th powers of nj. The optimal weighting matrix is dependent on the RD's and
(42) on Rjs (which can be computed from the RD's via (15)). In a
Taking the expected value of (42) will eliminate the terms practical situation we will replace the measured RD's into (48)
containing nj, and higher order terms are neglected. Thus which results in an estimate m*
of the optimalweighting
FRIEDLANDER: PASSIVE LOCALJZATION ALGORlTHM 239

matrix. Using the estimate W* rather than W* can be shown where

to have a relatively small effect on the resulting accuracy.
The Variance and Bias in Rotated Coorditlates
The error covariance matrix defines an “uncertainty el- are the elements of p,.
lipse.” As we will see in Section VI, this ellipse is typically Equation (5i) can be rewritten as
highly eccentric (i.e., its major axis is much larger than its
minor axis). Therefore, it is of interest to display the variance
andbiasin the directions of these axes, rather thau the where
prespecified (x, y ) coordinates. Next we derive the necessary
coordinate rotation.
Let
cov {is) = mOT r3lml -r2lm2
8= (52c)
be the singular value decomposition of the error covariance
- r 2 1 ( ~ 3 - ~ 1 ) + ~ 3 1 ( ~ 2’ - ~ 1 )
matrix, where As may be expixted, (11) establishes a singie straight LOP for
the source coordinates. In order to pinpoint the source location
along the LOP we use the first of two equations in (7a):
( ~ Z - - X I ) ~ , + ( Y ~ - Y I ) Y , = ~ I - R ~(53)
,~~I.
is a rotation matrix, and where
lriserting (52) into (53) and rearranging terms we get
1
Rh=- ~ ~ ~ ~ ~ - ~ l ~ + ~ ~ ~ - Y l ~ ~ l ~ , (54)
+ ~ Y z -
r2 1
The angle 8 gives the orientation of the major axis of the
uncertainty ellipse, is the standard deviation of the error in Note also that
the major axis direction, while ii2 is the staadard deviation of R1,=[(xl-x,)2+(Y1-~,)211/2
the error in the minor axis direction.
Finally, the bias in the rotated coordinates is given by = [ ( X ] -xs)2+ ( y ] -ax,- p)2]1’2

= [(I + a2)x5- 2(x1+ (Yl -Ob)&

+(yl-P)2+X;]’/2. (55)
In the case where the source is far from the sensors Squaring (54) md (55) and equating their left-hand sides gives
(compared to the maximum sensor separation), the major axis a quadratic equation for X,
oftheuncertainty ellipse is alongthe line connecting the y2x: + YlX, +y o =0 (56a)
sensors to the source, while the minor axis is perpendicular to
that line. In other words, a small array or cluster of sensors where
provides relatively poor range estimation, while it is capable y2+ ~ a 2 - { ~ ( x 2 - x ~ ) + ( y 2 - y 1 ) a l / r ~(56b)
l)2
of providing good bearing estimation. In the case of a sour&
located near or inside a cluster of sensors, the notions of range y1= -2(x1+(Y,-P)Cr+[(x2-x1)
and bearing are not welldefined, and it is not obvioos how the + ( u 2 - Y 1 ) ~ I [ ( Y 2 - Y I ) P - ~ l l / ~ ~ l ~( 5k)
uncertainty ellipse will be oriented. The formulas derived in
this section provide the answer. Yo=+ {[(Y2-YI)P-mlllr21)2. (564
IV. LOCALIZATION IN THE CASEOF A MINIMAL
NMER OF Thus the explicit formula for the source coordinates is
SENSORS
x,=t-ylf~lyR-4yoy21/(2~2) (574
The localization algorithm presented in Section 11 requires
that the number of sensors exceed the minimal number (n + ys=ax,+p. (5%)
i ) required to solve the problem. A question of considerable An error analysis of these estimates can be camed out by
practical interest is how to localize the source when only the expanding (x5,y,) in afirst-order Taylor expansion around the
minimal number of sensors (N = n + 1) is available. In this
true values of rzl, ~ 3 1 :
section we derive an explicit formula for the source location
for this case.
Consider the two-dimensional case (n = 2) with the first
sensor being the reference sensor ( j = 1). Equation (1 1)
reduces in this case to
[~~I(~Z-~I)--TZI(~~-XI)~~,+[~~I(Y~-YI)
- r 2 1 ( ~ 3 - y l ) l ~ 5 = r 3 ~ m ~ - -(514
2~~2 G
 . . .. ..
. .

240 IEEE JOURNAL OF OCEANIC ENGINEERING,

OE-12, VOL. NO. I , JANUARY 1987

where nl, n2 are the measurement errors corresponding to r21

and r31. Thus
COV {$} = GXG'. (59)
The explicit form of G is given in [121. The bias term can be
evaluated by a second-order Taylor expansion [12].
A very similar procedure can be applied in the three-
dimensional case (n = 3) to derive explicit formulas for the
source coordinates (xs,ys, zs) and to analyze their accuracy
properties. Equation ( 1 1 ) reduces in this case to

[r31(x2-Xl)-r21(x3-xI)Ixs+ [~31(~2-~1)-~21(~3-~1)I~s

AI B1
+[~31(z2-Z1)-~21(z3-Z1)1Zs=~31~1-~21~2
. (59a)
C1 Dl

~~41(~3-~1)-~31(~4-~I~I~s+~~41(~3-~l)-~3jl(~4-~l)l~s -YO
and therefore
-42 B2
+[r41(z3-Z1)-~31(24-Z1)1Zs=-41~2-~3l~3.
Y s = alxs + P I
c
2 D2
2s = a2x, +P 2 .
Multiplying the two equations above by C3 and C2,
respectively, and subtracting, we get An error analysis of the resulting estimates follows the same
steps as in thetwo-dimensional case [cf. the discussion
~ - D ~ C I (58)]:see [12].
( A I C ~ - A ~ C I ) X ~ + ( B I C ~ - B ~ C , ) ~ , = D I Csurrounding
or ESTIMMION
V. VELOCITY FROM RANGE AND RANGE-RATE
YS' a1xs + 01 (59b) DIFFERENCES
where In this section we derive an algorithm for estimating the
~l=(~lCz-A2CI)/(BlC2-~2CI) velocity of a movingsource when measurementsof range-rate
differences are available in addition to the range difference
PI=(~~C~-L)~C,)/(BIC~-BB~C I ) . measurements considered in earlier sections. In other words,
(59c)
Substituting (59b) in (59a), we get
given noisy measurementsof { rij, iij,i, j = 1, . * ,N)where -
f i j 4 8rjj/8t,we wish to estimate the source velocity vector x s

( A I + ~ ~ B ~ ) ~ ~ + C I Z ~ = D I - B I P , P axs/dt. As was mentionedbefore, range-rate measurements

are available in situations where the Doppler shifts of signals
or
emanating from the target can be estimated.
z, = a2xs f P 2 The basic equations for the source velocity are obtained by
CY^ = - ( A1 + Bl)/CI
differentiating (7a) with respect to time:
S.x
J s
=@.-R.
,J
. - R . '.
JSPJ J S ~ J
P2=(01 -BIPI)/CI. (59d)
= [diag { p j } - RjJ] pj- Rjspj (60)
Thus ( 1 1) yields a single straight LOP for the source
location. In order to pinpoint the source location along the
Qj
LOP, we use the first row of (7a)
where pj is the vector of range differences (see (7d)) and pi is
( x 2 - x ~ ) x s + ( ~ 2 - ~ ~ ) ~ s + ( z 2 - ~ 1 ) ~ s = m(59e)
~ - R ~ the
s r 2vector
~ . or range-rate differences, similarly defied.
As in Section II, we eliminate the term pj on the right-hand
Inserting (59c) and (59d) into (59e)and rearranging terms, we
side by premultiplying (60) by Mj (see (8)). Thus the source
get velocity vector is given by the least squares solution of the
1 following set of (N - 1) linear equations:
RIS = - - [(x2-x11 + 4 Y 2 -YJ + az(z2 - z1)lxs
r21
&f.S.i
J J s-
- M .J QJPJ
.'. (611
or
FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM 24 1

X 104 source - sensor

geolnetry
1 I 8 I I 1 I I I

0.a I.............. ............................. ...........................................

!
0. 6 .......................................... ._ ............., ............................ L ........................................................
0. 4 .................................................................................................................................
! .............
i30 i i : x ;
0. 2 - ~ .............” .....^.................................... ”
Ill i lo ;
e o ............................................. ............. ............. .........................................
t
e f 2 0 i f
-
r -0.2 .................................................................... ............. ...........
~ I ................................
S . . . : 40 ;
-0.4 -......................................... ~ ........................... ” ............. ......................................................
~

..
.............. ....................................................................................... * ............
-0.6 _............I........._..... n
. . .
-0.a ........................................... .......... .......-......,..............*.,. .-.................................... .............
I.

-1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
meters X 104
Fig. 4. Thesource-sensor geometry, four-sensor case.

Note that Qj, Mj are determined from the range-difference We note that a problem of considerable practical interest is
measurements. To evaluate Rjs, which is needed for Q,, we to localize a source from measurement of range-rate differ-
solve for x, using (11) or (20) and set RjS = llxj - xsll. An ences only(i.e., f j j is measured, but not rjs).The estimation of
alternative setof (N - 2) linear equations is obtained by the source location is considerably more complicated and it
replacing Mj in (61) by Mj (see (19)). does not seem possible to obtain a linear equation formulation
The analysis of the accuracy of the velocity estimates is in this case. However, some practical solution methods are
straightforward, provided that we assume that the errors in available (see [121).
measuring { rjj>and { f i j ) are mutually uncorrelated.
Let C, and Ci denote the covariance matrices of the VI. NUMERICAL
EXAMPLES
measurement noiseprocesses n, andni corresponding to pj and The performance of the proposedlocationestimation
pi, respectively. It follows from (60) that algorithm wasevaluated by extensive Monte-Carlo simula-
tions. In this section we present a few selected examples to
MjSj&=Mj(Qj+diag {n,.})(bj+ni). (63) illustrate thevalidityof the analytical expressions for the
covariance and bias of the estimates, and the effect of optimal
As before, we ignore the randomness in f i j and replace it by
weighting.
its true value M j . From (63) we get Three different sensor configurations were used: a four-
MjSj.&=Mj(nr o pj+Qjni+nr o ni) sensor cluster depicted in Fig. 4, a six-sensor cluster depicted
in Fig. 5, and an eight-sensor configuration depicted inFig. 6.
and Sensor number 1 was the reference sensor. The target
coordinates were (7000, 3000) in all cases (all location-related
MjSj COV { z ) S T M T = M j (diag {bj}E,diag { b j } parameters are specified in meters). Zero-meanGaussian
+ QjCiQJ+ X,Ci)MT. noise was addedto the range differences with variance chosen
so that all sensors have the same signal-to-noise ratio (SNR).
Thus The variance of the noise added to the ith sensor was c*I rilI,
COV { zs}=HjLjHJ resulting in S N R = l/a2for all sensors.
Table I summarizes the theoretical and sample statistics for
where (Z,,9s)obtained by 500 independent Monte-Carlo sums. The
mean-square error was computed in addition to the mean and
H, = (STMTM,S,) - Si’MjMj ’ standard deviation (mse = [mean2 + std2I1l2).
Lj = diag { b j } C, diag { b j } + QjCiQT+ C,Ci. The location estimates were computed by solving (11) with
M,, = (I - Z ) D (cf. (8)) in the nonweighted case and M,v =
The bias of the velocityestimator can be analyzed usingthe W1’”ll;ri (cf. (19), (46)) in the weighted case. The theoretical
sanie procedures employed in Section HI for evaluating the covariance wascomputedusing(27)with M,,, and M , as
bias of the location vector. We defer the details to [ 121. defined above. The bias was computed from (35) in the four-
 . -. . . . . .

242 I B E JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12, NO, 1, JANUARY 1987

X 104 source - sensor geoluctry

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

meters X 104
Fig. 5. Thesource-sensor geometry, six-sensor case.

X 104 source - s egnesooare t r y

1 1 , I I I 1

........... i..............i.............i...-........i........_..........-.._i.....-......i.........._.i ..-.......... i.._........_

0. 8 _.
..
0,6 _.........".i..............5 ............. .i
I
........4!....-.......i......._....i...-........i.......-._..i........_...i...._......_
i 1
4 :
0,4 -............ i i.........."..i............. $............ 3..........- I ............. L .............?'...
I ^ ........ i............-
io i I 0 ;

0. 2 _ ............ .............. ............. ......-..... ......._.... ........ ............. ............... ..........-........._

: x ;
i i
i 4
j
i j i."
i ".i
2; I
i I i

In
e
t
_
0 ............ i .............. i..... .....i............. 6............. L............ :.............:......0......i..............i ............-
Q

e 1 :i...-.........i:............_
-0. 2 _ ........... t ..............i.............:..............t
; s i i ......... i 4.............ii..............
r
s
-0,4 _...........i.............. i .............ii .......-....
o t
i............. 4:.........0. ; "..._.i...-.........i............_
2 .......-._.
16 i, i 0;
_
-0. 6 ............i..............i...._.......
a :
+.............i.............j............. i.............1..............:..............j .-.-......_
:
-0,8 _............i.............. i.............+-....-...... i......._...
j ......._.....
i .............{..............+.............:..........._
-1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
meters X 104
Fig. 6. Source-sensor geometry for the eight-sensor case.

sensor case and from (45) in the six- apd eight-sensor cases, from (35) are very close to the sample mean. The bias
with M,,, M , as above. In the four-sensor case, the weighted computed from (35) differed significantly from the sample
and unweighted estimates are identical. Therefore, results for mean when more than four sensors were used.
the weighted case are given only for the six- and eight-sensor c) In the six-sensor case the bias from (45) often exceeded
cases. The mean and standard deviations are also presented in the sample mean by a significant amount, but predicted
rotated coordinates with f3 being the angle of rotation. correctly its sign and order of magnitude. In test cases with a
The following observations can be made from Table I, and larger number of sensors the bias computed from (45) was
from other simulations not reported here. much closer to the sample mean.
a) There is excellent agreement between the sample and d) The optimallyweighted estimates were always more
theoretical variances of the estimates. accurate than the nonweighted estimates, but the difference
b) In the four-sensor case the theoretical values of the bias was quite small and not practically significant.
 FRIEDLANDER PASSIVE LOCALIZATION ALGORITHM 243

TABLE I
ERRORSTATISTICSFOR THE LOCATION ESTIMATES (ts,j s )AND THE
ROTATED ESTIMATES (2,y }

TABLE II
LOCALIZATION FOR DIFFERENT CHOICES OF THE REFERENCE SENSOR

1 reference I theoretical std. I theoretical std. I theorerical std. I theoretical std. I

sensor unweighted case weighted case unweighted case weighted case
location ?,I c1 ?,I is e o l z l g e o I z l g
15000.01 40.11 I 13.65 37.93 I 13.02 18.69 I 42.34 I 1.572 18.83 I 40.08 I 1.541
i250014330) ’ 33.04
16.48
11.32
7.038
31.40
13.90
10.84
6.105
18.83
22.95
34.91
17.89 1.019
18.97
23.53
33.20
1.170
15.15
1.139
3955
(-2500,4330)
(-5000.0) 19.24 9.090 18.17 21.26 25.20 8.561
25.13 20.06
.9300 .9100
(-2500:-4330) 17.43 8.085 16.59 7.670 24.77 19.20 ,8668 18.26 24.71 .8902
12500.-43301 17.45 8.173 16.57 7.841 24.97 18.31
,8668 19.25 25.22
9633

e) The bias was always small compared to the standard formulas for the covariance and bias of the resulting estimates.
deviation. Thus localization accuracy was dominated by the These formulas are very usefulfor performance prediction and
variance of the estimator. for system design. Using these formulas, the sensor configura-
f ) The uncertainty ellipse was highly eccentric in all test tion can be optimized for best accuracy for a given surveil-
cases, with the ratio of its axes being on the order of 50-200. lance scenario. Combining these results with formulas for the
g)Asmay be expected, the standard deviation is propo- variance of the range-difference measurements, we obtain a
tional to u, while the bias is proportional to u2. measure of localization accuracy in t e r m of system parame-
Next weconsider the effect of choosingthe reference sensor ters such as signal bandwidth, integration time, and sensor
on localization accuracy. Table II summarizes the theoretical SNR. In the passive sonar problem, simple formulas are
standard deviations of the estimator for the six-sensor case available to the asymptotic variance of the delay (range)
considered earlier, with u = 0.001. These results were estimates (see [ 131).
matched very closely by Monte-Carlo simulation results. The Tracking algorithms for localization from range-difference
first reference sensor was the rightmost sensor in Fig. 6. The measurements usually involves the use of the extended Kalman
other sensors followina counterclockwise direction. From filter (EKF) due to the nonlinearity of the equations relating
examination of Table II wenotethat choosing different the measurements to the state vector. The formulation pre-
reference sensors affects the accuracybyaconsiderablq sented here leads to a linear set ofequations making itpossible
amount. Usingthe expressions developed in this paper, we can to usea linear Kalman fiiter. For example, (11) may be
always choose thebest reference sensor for anygiven interpreted as a linear (although data-dependent) measurement
situation. Note also that the orientation of the uncertainty is equation for the state vector x, where MjSj is the measurement
affected only slightly by the choice of the reference sensor. matrix andMjpj is the “data vector. ” Replacing the EKF by a
linear Kalman filter has several advantages: it is computa-
VII. DISCUSSION tionally simpler andit avoids some of the convergence
In this paper we developed a localization algorithm which problems often encountered in the EKF.
requires solvinga linear set of equations, andpresented In our simulation study we found the localization algorithm
 244 IEEE JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12, NO. 1, JANUARY 1987

X 104
1

0.a

0. 6

0. 4

0.2

e
r -0.2
s
-0. 4

-0. 6
,
-0.8

-1
-1 -0.8 -0.6 - 0 . 4 -0.2 0 0.2 0.4 0.6 0.8 i
meter: X 104
Fig. 7. The effect of bias in the range difference measurements.

to havea relatively small bias (compared to the standard [5] J.J. Dongm, C. B. Moler, J. R. Bunch,and G. W. Stewart,
LZNPACKUser’s Guide. Philadelphia, PA: SIAM, 1979.
deviation of the estimates). This observation is true only if the [61 G . Golub and c. vanLoan, Matrk Complations. Baltimore, MD:
range-difference measurements are essentially unbiased. Any Johns Hopkins Univ. Press, 1984.
bias in the measurement translates into bias in the location [7] J. L. Poirot and G . V. McWilliams, “Applicationof linear statistical
models to radar location techniques,” IEEE Trans. Aerosp. Electron.
estimate. Equations (7) and (12) can be used to study the Syst., vol. AES-IO, no. 6, pp. 830-834, Nov. 1974.
sensitivity of the location estimate to bias in the range [8] W. M. Foy, “Position-location solutions by Taylor series estimation,”
measurements. IEEE Trans. Aerosp. Electron. Syst., vol. AES-12, no. 2, pp. 187-
194, Mar. 1976.
This sensitivity is much greater in the direction of the major [9] J. P.Van&en, “Navigation systems: Fundamentals oflow and very
axis ofthe uncertainty ellipse than inthe direction perpendicu- low frequency hyperbolic techniques,” Electrical Commun., vol. 45,
lar to it. To study this effect we introduce deterministic errors no. 3, pp.
192-212, 1970. ’

[lo] H. B. Lee, “A novel procedure for assessing the accuracy of

of the form

[
hyperbolic multilateration system,” IEEE Trans. Aerosp. Electron.
Syst., vol. AES-11,no. 1, pp. 2-15, Jan. 1975.
[I I] N. Marchand, “Error distributions ofbest estimate ofposition from
multiple time difference hyperbolic networks,” ZEEE Trans. Aerosp.
nj=a Nmig. Electron., vol. ANE-11, no. 2, pp. 96-100, June 1964.
1121B. Friedlander, “Analysis of some passive localization algorithms,”
Saxpy Computer Cop., Sunnyvale, CA, Tech. Rep. 86-010, Mar. 10,
1986.
to the measurements { ru} [cf. (22)]. For each choice of the [13] B. Friedlander, “On the Cramer-Rao bound for Doppler anddelay
error vector we compute a location estimate. A scatter plot of estimation,” ZEEE Trans. Znform. Theory, vol. IT-30, no. 3, pp.
575-580, May1984.
the 2(N-1)location estimates obtained in this way provides a
pictorial indication of the effects of bias. Fig. 7 presents this
plot for the six-sensor case with u = 0.03.
REFERENCES
*
[I] R. 0. Schmidt, “A new approach to geometry of range difference
location,” ZEEE Trans. Aerosp. Electron. Syst., vol. A E S -
,
I 6,
no. $;?,$ Benjamin Friedlander (S’74-M’76-SM’SZ-F’87)
NOV.
821-835, pp. 1972. was born on February 24, 1947. He received the
[2] J. M. Delosme, M. Morf, and B. Friedlander, “Estimating source y;‘; B.Sc.
and the M.Sc. degrees in electrical engineer-
location from timedifference-of-arrival: A linear equation approach,” .* ing from the Technion, Israel Institute of Technol-
Information Systems Lab., Stanford, CA, Tech. Rep. M355-1, Mar. 1 ogy in 1968and 1972, respectively, and the Ph.D.
31, 1979. r degree in electrical engineering and the M.Sc.
[3] J. M. Delosme, M. Morf, andB. Friedlander, “A linear equation .’;. ’ degree in statistics from Stanford University, Stan-
approach to locating sources from time-difference-of-arrival measure- ford, CA,
1976.
in
ments,” in Proc. ZEEE Int. Con$ Acoustics, Speech, and
Signal From 1968 to 1972 he served in the Israel
Processing (Denver, CO), 1980, pp. 818-824. Defence Forces as an Electronic Engineer. From
[4] J. S. Abeland J. 0. Smith, “The spherical interpolation method to 1976 to 1985 he was at Systems Control Technot-
closed form passive source localization using range difference measure- ogy, Inc., Palo Alto, CA. He was responsible for a number of research and
ments,” in Proc. IEEE Int. Conf. Acoustics, Speech, Signal development projects in the areas of signal processing and control, and was
Processing (Dallas, TX), 1987,pp. 471474. the
Managerof the Advanced
Technology
Division.
During this period
was
he
FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM 245

also a Lecturer at Stanford University, teaching graduate courses on linear Dr. Friedlander was an Associate Editor of the IEEE Transactions on
systems,nonlineardetectionandestimation,andsystemidentification.In Automatic Control in 1984, a member of theEditorialBoard of the
Novembcr of 1985 he joined Saxpy Computer Corporation, Sunnyvale, CA, InternafionalJournal of Adoptive Controland Signal Processing,and is a
asthe Director ofAdvancedTechnology. He is responsible for theR&D member of the Administrative Committee of the IEEE Acoustics, Speech, and
activities ofthecompanyand for thedevelopment of signalprocessing Signal Processing (ASSP) Society. He serves as a member of the Technical
applications for the MATRJX I-a 1000-Mflopsupercomputer. He has Committee on Spectrum Estimation of the ASSP, isand the Vice-chairman of
authored more than 170 publications in the areas of signal processing and the Bay Area Chapter of that society. He is a member of the IEEE Admission
estimation.His currentinterests includeparallel and systolicprocessing and Advancement Committee. He is the recipient of the 1983 ASSP Senior
architectures andtheirVLSIimplementations,advancedprocessingtech- Award for the paper “Recursive Lattice Forms for Spectral Estimation,” and
niques for underwater surveillance,high-resolution spectral analysis andarray the 1985 award for best paper of the year from the European Association for
processing, adaptive filtering, radar and communication processing, detection Signal Processingfor the paper“On the Computation of an Asymptotic Bound
tracking and localization of multiple targets, image processing, and knowl- for Estimating Autoregressive Signals in White Noise.” He is a member of
edge-based signal processing. Sigma Xi.

View publication stats