.

. . . .

..

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-The problem of estimating solrrce 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 variance bias and of the estimator validated are by Monte-Carlo simulation. Theproblem of estimating source velocity givenmeasurements of range differences and range-rate differences is briefly considered, and a lidear equation technique is derived.

- K +
source sensor 1

2 s
%sensor

I. INTRODUCTION
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 measurement consists of observing differences in the time of 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. The conventional approach to estimating source location is based on finding the hyperbolic lines of position (LOPS);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 drawbacks: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 from range differences and their rates of change (cf. Section V). Furthermore, the complexity of the solutions makes error analysis very difficult. An alternative approach, which circumvents many of these difficulties, was proposed in 113. This approach is based onthe ideathat three sensors with their set of range differences determine a straight line of position. In fact, this line is the major axis o a general conic which f passes through the sensors. When more than three sensors are available, several straight LOPs are generated and their intersection provides an estimate of the source location (see Fig. 2).
Manuscript received May 5 , 1986; revised September 5 , 1986. This work was supported in part by theOffice of NavalResearchunderContract N00014-87-C-0147. The author is with Saxpy Computer Corp., Sunnyvale, CA 94086. IEEE Log Number 8714343. 1 We are considering here the two-dimensional case, i.e., the source and sensors are located in a plane.

r13

r12

Fig. 1. Hyperbolic lines of position.

Sensor

Fig. 2.

Location on the conic axis.

The method proposed in [11 computed the intersection point from a set o() P (A?/(N - 3)!3!) linear equations, Nbeing f: the number of sensors. This approach alleviates some of the difficulties of the hyperbolic LOP methd. However, it introduces some inconsistencies due to its use of redmdant information: this problem formulation leads to a set of ( ; ) equations while, In fact, there are only N available measurements (see [2] for a more detailed discussion). A different formulation is proposed i [2] and [3] based on n the straight LOP approach, which leads to a set of N equations. This set of equations contains all the relevant sensor data without redundancy. However, due to a noalinear relationship between certain variables in these equations, the location can no longer be computed by a linear estimation procedure. An iterative gradient search procedure is proposed in [3] to compute the the source location.

0364-9059/87/0100-0234$01.00

0 1987 IEEE

FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM

235
The various quantities defined above are depicted in Fig. 3. The localization problem is to estimate x, given noise measurements of the RDs { rij,i, j = 1, . ,N)provided by N sensors. Note that there are (f)= (N(N- 1)/2) distinct RDs. However, all of these RDs can be completely determined from only (N- 1) RDs (e.g., { ril, i = 2, * -, N)) the noiseless case. in Using the notation above, we can derive the basic equations in a straightforward manner. Note that

Subtracting one of the equations in [2] and [3] from the (N 1) remaining equations leads to a set of (N- 1) equations containing a single nonlinear term. It was recently noted in [4] that this nonlinearity can be eliminated by a two-step least squares procedure to yield a closed-form expression for the source location. This noniterative technique appears to be very attractive both in its modest computational requirements and its performance. In this paper we present a localization technique based on the solutionof a set of either (N - 1) or (N - 2) linear equations. While the derivation and resulting solution appear different from the onein [4], it isshownthat the two are mathematically equivalent. However, themethod presented here is somewhat simpler and reveals more clearly the underlying structure of the solution. The main contribution of this paper isin providing an accuracy analysis of the proposed solution method. In Section III we derive approximate expressions for the bias and variance of the location estimates. The variance formula is used to derive an optimalweighting matrix for the least squares location estimator. In Section VI we verify these formulas by Monte-Carlo simulations. The methodspresentedin [1]-[4] and Section in II all require that the number of sensors N be greater than (n + 1) where n is the dimension of the space in which thesensors and source are located (in other words, N > 3 in the twodimensional case and N > 4 in thethree-dimensionalcase). In some practical situations no more than N = n + 1 sensors are available. In Section IV we derive explicit expressions for the source Iocation for this case. These expressions are somewhat 1, and their morecomplicatedthan the one for N > n 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 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.
-

R ~ = ( r u + R j s ) 2 = ~ + 2 R j s r u + R ~ s (2)
and also that

R ~ = I I X ~ - X , ~ ~ ~ = R ~ ~ - ~ X ~ X ~ + R (3)

Subtracting (3) from (2) and rearranging terms, we get

2 x , ? x s = R ~ - r ~ . - 2 R j S r u + R & -JS 2 R
Setting i
=

(4)

j in (4) gives
2xjX. = R;o Ri0 - R;,

(5)

where we used the fact that rjj = 0. Subtracting (5) from (4) yields

2 ( ~ ~ - x ~ > ~ x ~ = ( R ~ - R ~ , J - r ~ ~(6)2 R ~ , which can be written in matrix form as (74

, (N-l)x3

11. LOCALIZATION

FROM RANGE DIFFERENCES.

In this section we state and solve the localization problem for the three-dimensional case. The discussion here and in the following sections will be in terms of range difference measurements. However, the same results apply to localization from delay difference measurements. Delay be can converted into range simply by multiplying it by the speed of propagation. Let xi = [xi,y i , zi] denote the (x, , z ) coordinates of the y ith sensor and x, = [xs, s , z,] denote the source coordinates. y 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) sensors i and j where for

1
j=

R;-lo R;,lo

-Rjo -R;o

-rj-.lj - r j i l j , ( N - l ) x 1 (7c)

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
z

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

least squares equation solver on the linear set of (1 l), see for example [5] and [6]. A closed-form solution which is useful for some subsequent derivations is given by
Xi sensor i
X ,

i l
Y
5S .

= (STMTMj S,) - STMTMjpj.

(12)

Fig. 3. Sensor-source geometry.

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 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. First it is assumed that Rjs is known, in which case a least squares solution for x, is given by
X,=(STS)-~S~T(,U~Rjspj).

Note that given the measurements { rij, i = 1, N}, the vectors p j , p j are known; so is the matrix Si, which depends only on the sensor locations (assumed to be hown). The unknown quantities in (7a)-(7d) are the location vectors x, and the distance Rj, between the source and the reference sensor, which is a function of (the unknown) x,. It is the presence of Rjs that complicates the solution of (7). Our approach is based on the idea of 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

--

(13)

Next Rjs is chosen to minimize IISjxs - ( p j - Rsjpj)ll or equivalently 11 Pj(pj - Rsjpj)ll where

pj=r-sj(s;sj)-lsj.

(14)

It is straightforward to show that the minimizing Rsj is given by

My= (I- P ) D j
where

(8a)

Thus the method presented in [4] computes the location vector by the least squares solution of

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 DjPjl and
(I-Zk)l=l-Zkl=1-1=0.

P [l, 1,

. e - ,

1]T

An Alternative Form of (11) In the following section we w find it convenient to use a l i (9) slightly different form of the location estimator. Note that the matrix (I - 29 [cf. (8a)l is an (N- 1) X (N- 1) matrix withrank ( N - 2). Thus the singular valuedecomposition (10) (SVD) of this matrix whave the form l i

To simplify the subsequent discussion, we assume that k = 1 and suppress the superscript in My). From the property of Mj discussed above, it follows that
M j S j X , =Mj pj

(I- Zk) r U k Y =

Ukl
TN-2

(11)
=

The location vector x, can now be computed by using good any

u diag {T:, k

- y

9i-J V l

(18)

FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM

237 or

where uk, v k are (N- 1) X (N- 2) orthogonal matrices are the (N- 2) nonzero singular Values and { q f, * , of ( I - z ) Since ( I - z ~ = 0, it foUows that VlI = 0. k. I Therefore, another candidate for the matrix Mj whichwas used to annihilate the second term in (7a) is

cov {fs} =4cFj

(27a)

where F j = (STMTM,Sj)- STMTMj (diag { p j } + RjsI). (27b) It should be noted that by performing a first-order Taylorseries expansion of 2, around x, it follows that Fj is the matrix of the first derivatives of 2 with respect to the measurements , { r g } , Le., F j = ( d i S / d p j ) . Thus (27) can be derived by differentiating (1 1). The derivation of the bias can be done by a second-order Taylor-series expansion of $. It is straightforward to show that

A , V,Dj. ?=

(1 9)

Premultiplying (7a) by this matrix leads to a set of (N- 2) linear equations for x, M.S.X -iff. JPJ . J J s(20) whereas (11) consisted of (N - 1) linear equations. In the following we will suppress the subscript k in V k . 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 RDs are corrupted by additive mutually uncorrelated zeromean white Gaussian noise with covariance matrix C=diag { , a)
u2=[a:, a : ,

and similarly

The second-order derivative matrices of x, and y, are presented somewhat complicated (cf. [12])andwillnotbe b.- P j + n j J(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 bj, Inserting is, and Ris into (7a) we get sides of (24) and premultiply by Mj to get The measured RDs are denoted bj and

-.., u i - , ] .

(21)

SjZ, = Sjx, + S,f,

= pj - pi

o nj + Rj,pj

+-1 nj o nj+Rj,nj)
2

M,SjE{fs}= -Mj
(23)
two

[I: 3
E
- nj

o nj + E { R j , n j }

1.

(28)

where denotes element-by-element multiplication of vectors. Alternatively

(pi+

1 Rj,l) o nj+ R j s p j + - nj o nj+ pjsnj] . 2 (24)

Note that in this case we do not neglect the last two terms of (24). In fact, it is precisely these t e r n that introduce the bias. The multiplication. by Mi is necessary to eliminate the term E{RjsIpjThe first term is the square brackets above is simply u2.To evaluate the second term we insert 6j into (15) to obtain (after reordering of terms)

To evaluate the covariance matrix of fswe make some approximations. Assuming small error conditions, wewill neglect any terms which contain products of errors. Thus the two last terms in (24) are neglected. Next note that

To evaluate R,j we will make asmall error approximation and neglect al terms containing products of noise variables (i.e., l the last three terms in the numerator and the last term in the denominator). Next we use the standard approximation

MjSj.fs = -Mj(p,+R,,l) o nj
from which it readily follows that

(25) to get

1 pTPjpj + 2pipin;

- -(1 1
pTPjpj

- 2 p?Pjnj PTPjPj

-->

(30)

. .

. .

. .

238
Note that Pj (L - Rjspj) =Pj( Six,) = 0 ,, /

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

A K = E{

MTa. =MTM. ^.}

j
JPJ
J

JPJ

(32)

1 MTMj,2 2 J

and therefore the second term in (31) reduces to - RjspTnj.It follows that

+MJMjn2+ diag { DjMTMjDJ:I} u2

+ diag { DjMTMjDj)u2o p j .
It fbllows that

(43)

E . , =H-'STK. {?}
Since 2, = x, Inserting (33) in (28) we get

(4 4)
(45)

+ f,, we get
BIAS=E(fs}=H-'SjK-~,.

Accuracy Optimization In the discussion so far the location equation (1l ) , or or equivalently (20),was assumedto be solved by astandard least squares technique. It is possible, of course, to use instead a BIAS = E { f,} = - (SiTMjMjSj)- STMTMj - g2+ 6j weighted least squares procedure. Let W be a positive definite (N - 2) X (N - 2) weighting matrix and premultiply (20) by (35) 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. or From (22) and (8) it follows that 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 (27) and (34) or (35). , in 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 answer this question we w first rewrite the covariance l i Next note that equation (27) for the weighted case as = MTMj -MTMjDjNj - NjDjMTMj covw{fs} =(ATWA)-1A7WCWA(A*WA)-1 (47a) +NjDjMTMjDjN,. (38) where From (12) it follows that A iGtjSj (4%) 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 Next we show that of (39) by w *= (48) E { $ } = [E(S;h?TA?jS'}]-'E{ STi@TA?jfij}. (40)

'

1.

aa Tj

s;@;.i@jsj]

c-'

This approximation is a good one when SiTM,TMiSj is a wellconditioned matrix and the measurement noise is small. From (38) it follows that

is the optimal weighting matrix in the sense that

COVw{$}2COVw

* {fs}=(ATC-'A)-l

(49)

H 4 E { STMTMjS'} = SjTM,TMjSj

for all permissible choices of W. The proof consists of verifying that

+ ST diag (DjMTMjDj}C Sj .
Next note that 1 ,@T$j.^.-MTM. J 2 MTM.n. o nj JPJ- J J ~ J J J

(41)

COV,{f,}

-covw

* {%}

= [(A 'WA)-'A

TW-(A'C-lA)-lATC-l] C[(ATWA)-'ATW-(A*C-'A)-lATC-l]'. (50)

Since C > 0 the right-hand side of (50) is clearly nonnegative definite and (49) immediately follows. +MTMjDjNjpj o nj Note that the derivation above requires that C and, +NjDjMTpj o nj+ NjDjMTMjDjNjpj therefore, Mj be nonsingular. This is the reason we used (20) rather than (11) in the development above (i.e., Mj is +terms containing nj 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) which results in an estimate of the optimalweighting containing nj, and higher order terms are neglected. Thus

m*

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,. Equation (5i) can be rewritten as lipse. As we will see in Section VI, this ellipse is typically 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 r3lml -r2lm2 cov {is) = mOT 8= (52c) -r21(~3-~1)+~31(~2-~1) be the singular value decomposition of the error covariance As may be expixted, (11) establishes a singie straight LOP for matrix, where the source coordinates. In order to pinpoint the source location along the LOP we use the first of two equations in (7a): is a rotation matrix, and where

(~Z--XI)~,+(Y~-YI)Y,=~I-R~,~~I. (53)
lriserting (52) into (53) and rearranging terms we get

The angle 8 gives the orientation of the major axis of the uncertainty ellipse, is the standard deviation of the error in the major axis direction, while ii2 is the staadard deviation of the error in the minor axis direction. Finally, the bias in the rotated coordinates is given by

Rh=r2 1

~ ~ ~ ~ ~ - ~ l ~ + ~ ~ ~ - Y l ~ ~ l ~ , (54) Y z + ~

Note also that

R1,=[(xl-x,)2+(Y1-~,)211/2
= [ ( X ] -xs)2+ ( y ] -ax,= [(I
p)2]12

+ a2)x5- 2(x1+ (Yl -Ob)&

(55)

+(yl-P)2+X;]/2.
In the case where the source is far from the sensors (compared to the maximum sensor separation), the major axis oftheuncertainty ellipse is alongthe line connecting the sensors to the source, while the minor axis is perpendicular to that line. In other words, a small array or cluster of sensors provides relatively poor range estimation, while it is capable of providing good bearing estimation. In the case of a sour& located near or inside a cluster of sensors, the notions of range and bearing are not welldefined, and it is not obvioos how the uncertainty ellipse will be oriented. The formulas derived in this section provide the answer.

Squaring (54) md (55) and equating their left-hand sides gives a quadratic equation for X,

y2x: + YlX,
where

+y o =0

(56a)

y2+ ~ a 2 - { ~ ( x 2 - x ~ ) + ( y 2 - y 1 ) a l / r ~(56b) l)2

y1=

-2(x1+(Y,-P)Cr+[(x2-x1)
+ ( u 2 - Y 1 ) ~ I [ ( Y 2 - Y I ) P - ~ l l / ~ ~ l ~(

5k)

Yo=+

{[(Y2-YI)P-mlllr21)2.

(564

x,=t-ylf~lyR-4yoy21/(2~2) (574 The localization algorithm presented in Section 1 requires 1 ys=ax,+p. (5%) that the number of sensors exceed the minimal number (n + 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

IV. LOCALIZATION IN THE CASE A MINIMAL OF NMER SENSORS

OF

Thus the explicit formula for the source coordinates is

[~~I(~Z-~I)--TZI(~~-XI)~~,+[~~I(Y~-YI)
- r 2 1 ( ~ 3 - y l ) l ~ 5 = r 3 ~ m ~ - -(514 2 2~~

. .

..

..

. .

240

IEEE JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12,

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 threedimensional 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

B 1

+[~31(z2-Z1)-~21(z3-Z1)1Zs=~31~1-~21~2
C 1

(59a)
-YO

Dl
and therefore

~~41(~3-~1)-~31(~4-~I~I~s+~~41(~3-~l)-~3jl(~4-~l)l~s

-42

B 2
Y s = alxs + P I
2s = a2x,

+[r41(z3-Z1)-~31(24-Z1)1Zs=-41~2-~3l~3.

c 2

D2

+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 the two-dimensional case [cf. the discussion surrounding ( A I C ~ - A ~ C I ) X ~ + ( B I C ~ - B ~ C , ) ~ , = D I C ~ - D ~ C I (58)]:see [12].
a1xs

or
YS'

+ 01

(59b)

V. VELOCITY ESTIMMION RANGE AND RANGE-RATE FROM DIFFERENCES

where

In this section we derive an algorithm for estimating the velocity of a movingsource when measurementsof range-rate ~l=(~lCz-A2CI)/(BlC2-~2CI) differences are available in addition to the range difference P I = ( ~ ~ C ~ - L ) ~ C , ) / ( B I C ~ - B B ~ C I ) . measurements considered in earlier sections. In other words, (59c) given noisy measurements { rij, iij, j = 1, . * ,N)where of i, Substituting (59b) in (59a), we get 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 differentiating (7a) with respect to time: CY^ = - ( A1 + Bl)/CI

P2=(01

-BIPI)/CI.

(59d)

S.xs [email protected]. . - R J. S '~. J J ,J JSPJ

= [diag { p j } - RjJ]
Qj pj-

Thus ( 1 1) yields a single straight LOP for the source

location. In order to pinpoint the source location along the LOP, we use the first row of (7a)

Rjspj

(60)

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) ~ the 2 ~ . or range-rate differences, similarly defied. ~ - R s r vector As in Section II, we eliminate the term pj on the right-hand Inserting (59c) and (59d) into (59e) rearranging terms, we and 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 r1 2 &f.S.i - M .J QJPJ .'. J J s(611 or

FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM

X

24 1
source
I
8

104

- sensor geolnetry
I

0.a
0. 6
..................... .....................

I . . . . . . . . . . . . . .

.............................
L

...........................................
.............
x ;

._

............. ............................ ,
i30

! ........................................................

0. 4
0. 2
Ill

! .................................................................................................................................

............. .....^....................................
i lo

e
t
S

o ............................................. .............

;
i
f
~

............. .........................................
I

r -0.2 ....................................................................
~

. . . -0.4 -.........................................
-0.6
_............I........._.....

f 2 0

: 40 ;

............. ...........

................................

...........................

............. ......................................................
*

.. .............. ....................................................................................... . . .
n
I.

............

-0.a

........................................... .......... .. .....-......,.............*,. .-.................................... ............. . .

-0.8
-0.6

-1
-1

-0.4

-0.2

0 meters

0.2

0.4

0.6

0.8
X

1
104

Fig. 4. Thesource-sensor geometry, four-sensor case.

Note that Qj, Mj are determined from the range-difference measurements. To evaluate Rjs, which is needed for Q,, we solve for x, using (11) or (20) and set RjS = llxj - xsll. An alternative setof (N - 2) linear equations is obtained by replacing Mj in (61) by Mj (see (19)). The analysis of the accuracy of the velocity estimates is straightforward, provided that we assume that the errors in measuring { rjj>and { f i j ) are mutually uncorrelated. Let C, and Ci denote the covariance matrices of the measurement noise processes n, andni corresponding to pj and pi, respectively. It follows from (60) that MjSj&=Mj(Qj+diag {n,.})(bj+ni).

We note that a problem of considerable practical interest is to localize a source from measurement of range-rate differences only(i.e., f j j is measured, but not rjs).The estimation of the source location is considerably more complicated and it does not seem possible to obtain a linear equation formulation in this case. However, some practical solution methods are available (see [121).

VI. NUMERICAL EXAMPLES

The performance of the proposed location estimation algorithm wasevaluated by extensive Monte-Carlo simulations. In this section we present a few selected examples to illustrate thevalidityof the analytical expressions for the covariance and bias of the estimates, and the effect of optimal weighting. Three different sensor configurations were used: a foursensor cluster depicted in Fig. 4, a six-sensor cluster depicted in Fig. 5, and an eight-sensor configuration depicted inFig. 6. Sensor number 1 was the reference sensor. The target coordinates were (7000, 3000) in all cases (all location-related parameters are specified in meters). Zero-mean Gaussian noise was addedto the range differences with variance chosen so that all sensors have the same signal-to-noise ratio (SNR). The variance of the noise added to the ith sensor was c*I ril I, resulting in S N R = l/a2for all sensors. Table I summarizes the theoretical and sample statistics for ( , 9s)obtained by 500 independent Monte-Carlo sums. The Z, mean-square error was computed in addition to the mean and standard deviation (mse = [mean2 + std2I1l2). The location estimates were computed by solving (11) with M,, = (I - Z ) D (cf. (8)) in the nonweighted case and M,v = W1ll;ri (cf. (19), (46)) in the weighted case. The theoretical covariance wascomputedusing(27)with M,,, and M , as defined above. The bias was computed from (35) in the four-

(63)

As before, we ignore the randomness in f i j and replace it by its true value M j . From (63) we get

MjSj.&=Mj(nr o pj+Qjni+nr o ni) and MjSj COV { z ) S T M T = M j (diag {bj}E,diag { b j }

+ QjCiQJ+
Thus
COV {

X,Ci)MT.

zs} HjLjHJ =

where

H, = (STMTM,S,) - SiMjMj
Lj = diag { b j } C, diag { b j } + QjCiQT+ C,Ci.

The bias of the velocityestimator can be analyzed usingthe sanie procedures employed in Section H for evaluating the I bias of the location vector. We defer the details to [ 121.

-.

. . .

242
X 104

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

source

- sensor

geoluctry

-1

-0.8

-0.6 -0.4

-0.2

0.2

0.4

0.6

0.8

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

X 104

104
1

source

- s egnesooare t r y
,
I

........... i..............i.............i ...-........ ........_..........-.._.....-......i.........._.i..-.......... i.._ i i ........_ 0. 8 _.

i 0,6 _........."...............5 ............. .i
0,4 -............ i
i

..

........4!....-.............._.......-........i.......-._.. i i i........_... i...._......_

i
1
^

4 : I i.........."..i............. $............ 3..........- I ............. L ............. ........ i............?'...

io
i
I
0 ;

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

i
j

i 4

i."

" i .

2;

In
e t

r
s

_ ; s i..............t i -0. 2 _ ........... t ..............i.............:

Q

................... ......i..............i ............: 0 0 ............ i .............. i..... .....i............. 6............. L............ :

......... i

1 i...4.............ii..............: .........i............_ :

io t : 0 .......... i i............. 4.......... ; "..._. i...-....................._ i -0,4 _. i.............. i ............. .......-....2 .......-._. 16 i i 0; , a : i + i .....j .............1..............:..............j .-.-......_ -0. 6 ..........................i...._.................... ........ ............. i

......................... i............. -....-...... i......._... i + j ......._..... .......................... .............: i { .+ ..........._ -0,8 _.

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

0.4

0.6

0.8
X

1
104

meters

Fig. 6 Source-sensor geometry for the eight-sensor case. .

sensor case and from (45) in the six- apd eight-sensor cases, with M,,, M , as above. In the four-sensor case, the weighted and unweighted estimates are identical. Therefore, results for the weighted case are given only for the six- and eight-sensor cases. The mean and standard deviations are also presented in rotated coordinates with f3 being the angle of rotation. The following observations can be made from Table I, and from other simulations not reported here. a) There is excellent agreement between the sample and theoretical variances of the estimates. b) In the four-sensor case the theoretical values of the bias

from (35) are very close to the sample mean. The bias computed from (35) differed significantly from the sample mean when more than four sensors were used. c) In the six-sensor case the bias from (45) often exceeded the sample mean by a significant amount, but predicted correctly its sign and order of magnitude. In test cases with a larger number of sensors the bias computed from (45) was much closer to the sample mean. d) The optimallyweighted estimates were always more accurate than the nonweighted estimates, but the difference was quite small and not practically significant.

FRIEDLANDER PASSIVE LOCALIZATION ALGORITHM

243
TABLE I

ERROR STATISTICS FOR THE LOCATION ESTIMATES ROTATED ESTIMATES (2, } y

(ts, s )AND j

THE

TABLE II
LOCALIZATION FOR DIFFERENT CHOICES OF THE REFERENCE SENSOR

1 reference
sensor location
15000.01 i250014330) (-2500,4330) (-5000.0) (-2500:-4330) 12500.-43301

I theoretical std. I theoretical std. I

unweighted case

?,I

weighted case

c1
I
13.65 11.32 7.038 9.090 8.085 8.173

?,I

theorerical std. unweighted case

theoretical std. weighted case

is

40.11 33.04 16.48 19.24 17.43 17.45

37.93 31.40 13.90 18.17 16.59 16.57

e o l z l g
18.69 18.83 22.95

e o I z l g
18.83 18.97 23.53 8.561 25.13

13.02 10.84 6.105 7.670 7.841

I 42.34 I 1.572

I 40.08 I 1.541
33.20 1.170 15.15 20.06 .9300 1.139 3955 .9100 ,8668 18.26 24.71 .8902

34.91 17.89 1.019 21.26 25.20 24.77 19.20 24.97 19.25 9 6 3 ,8668 18.31 25.22 3

e) The bias was always small compared to the standard deviation. Thus localization accuracy was dominated by the variance of the estimator. f ) The uncertainty ellipse was highly eccentric in a l test l cases, with the ratio of its axes being on the order of 50-200. g)Asmay be expected, the standard deviation is propotional to u, while the bias is proportional to u2. Next weconsider the effect of choosingthe reference sensor on localization accuracy. Table II summarizes the theoretical standard deviations of the estimator for the six-sensor case considered earlier, with u = 0.001. These results were matched very closely by Monte-Carlo simulation results. The first reference sensor was the rightmost sensor in Fig. 6. The other sensors followina counterclockwise direction. From examination of Table II we that note choosing different reference sensors affects the accuracy a by considerablq amount. Usingthe expressions developed in this paper, we can always choose the best reference sensor for any given situation. Note also that the orientation of the uncertainty is affected only slightly by the choice of the reference sensor.
VII. DISCUSSION

In this paper we developed a localization algorithm which requires solving linear set of equations, and a presented

formulas for the covariance and bias of the resulting estimates. These formulas are very usefulfor performance prediction and for system design. Using these formulas, the sensor configuration can be optimized for best accuracy for a given surveillance scenario. Combining these results with formulas for the variance of the range-difference measurements, we obtain a measure of localization accuracy in t e r m of system parameters such as signal bandwidth, integration time, and sensor SNR. In the passive sonar problem, simple formulas are available to the asymptotic variance of the delay (range) estimates (see [ 131). Tracking algorithms for localization from range-difference measurements usually involves use of the extended Kalman the filter (EKF) due to the nonlinearity of the equations relating the measurements to the state vector. The formulation presented here leads to a linear set ofequations making itpossible to usea linear Kalman fiiter. For example, (11) may be interpreted as a linear (although data-dependent) measurement equation for the state vector x, where MjSj is the measurement matrix andMjpj is the data vector. Replacing the EKF by a linear Kalman filter has several advantages: it is computationally simpler and it avoids some of the convergence problems often encountered in the EKF. In our simulation study we found the localization algorithm

244
X

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

104

1
0a .

06 .
0. 4

02 .

r -0.2 s
-0. 4

-0. 6

,
-0.8
-1

-1

-0.8

-0.6 - 0 . 4

-0.2

0.2

0.4

0.6

0.8
X

i
104

meter:
Fig. 7. The effect of bias in the range difference measurements.

to havea relatively small bias (compared to the standard deviation of the estimates). This observation is true only if the range-difference measurements are essentially unbiased. Any bias in the measurement translates into bias in the location estimate. Equations (7) and (12) can be used to study the sensitivity of the location estimate to bias in the range measurements. This sensitivity is much greater in the direction of the major axis ofthe uncertainty ellipse than inthe direction perpendicular to it. To study this effect we introduce deterministic errors of the form

nj=a

to the measurements { ru} [cf. (22)]. For each choice of the error vector we compute a location estimate. A scatter plot of 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

[5] J. J. Dongm, C. B. Moler, J. R. Bunch, and G. W. Stewart, LZNPACKUsers Guide. Philadelphia, PA: SIAM, 1979. [61 G . Golub and c. vanLoan, Matrk Complations. Baltimore, MD: Johns Hopkins Univ. Press, 1984. [7] J. L. Poirot and G . V. McWilliams, Applicationof linear statistical models to radar location techniques, IEEE Trans. Aerosp. Electron. Syst., vol. AES-IO, no. 6, pp. 830-834, Nov. 1974. [8] W. M. Foy, Position-location solutions by Taylor series estimation, IEEE Trans. Aerosp. Electron. Syst., vol. AES-12, no. 2, pp. 187194, Mar. 1976. [9] J. P.Van&en, Navigation systems: Fundamentals oflow and very low frequency hyperbolic techniques, Electrical Commun., vol. 45, no. 3, pp. 192-212, 1970. [lo] H. B. Lee, A novel procedure for assessing the accuracy of 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. Nmig. Electron., vol. ANE-11, no. 2, pp. 96-100, June 1964. 1121 B. Friedlander, Analysis of some passive localization algorithms, Saxpy Computer Cop., Sunnyvale, CA, Tech. Rep. 86-010, Mar. 10, 1986. [13] B. Friedlander, On the Cramer-Rao bound for Doppler anddelay estimation, ZEEE Trans. Znform. Theory, vol. IT-30, no. 3, pp. 575-580, May1984.

*
Benjamin Friedlander (S74-M76-SMSZ-F87) was born on February 24, 1947. He received the B.Sc. and the M.Sc. degrees in electrical engineer.* ing from the Technion, Israel Institute of Technol1 ogy in 1968 and 1972, respectively, and the Ph.D. r degree in electrical engineering and the M.Sc. .; . degree in statistics from Stanford University, Stanford, CA, in 1976. From 1968 to 1972 he served in the Israel Defence Forces as an Electronic Engineer. From 1976 to 1985 he was at Systems Control Technotogy, Inc., Palo Alto, CA. He was responsible for a number of research and development projects in the areas of signal processing and control, and was the of Manager the Advanced Technology Division. During this period was he

R. 0. Schmidt, A new approach to geometry of range difference location, ZEEE Trans. Aerosp. Electron. Syst., vol. A S 6, E, no. I NOV. 821-835, pp. 1972. B. Friedlander, Estimating source [2] J. M. Delosme, M. Morf, and location from timedifference-of-arrival: A linear equation approach, Information Systems Lab., Stanford, CA, Tech. Rep. M355-1, Mar. 31, 1979. [3] J. M. Delosme, M. Morf, and B. Friedlander, A linear equation approach to locating sources from time-difference-of-arrival measurements, i Proc. ZEEE Int. Con$ Acoustics, Speech, n and Signal Processing (Denver, CO), 1980, pp. 818-824. [4] J. S. Abeland J. 0. Smith, The spherical interpolation method to closed form passive source localization using range difference measurements, in Proc. IEEE Int. Conf. Acoustics, Speech, Signal Processing (Dallas, TX), 1987, pp. 471474.

[I]

$;?,$ y;;

FRIEDLANDER: PASSIVE LOCALIZATION ALGORITHM

245
Dr. Friedlander w s an Associate Editor of the IEEE Transactions on a Automatic Control in 1984, a member of the Editorial Board of the InternafionalJournal of Adoptive Control and Signal Processing,and is a member of the Administrative Committee of the IEEE Acoustics, Speech, and Signal Processing (ASSP) Society. He serves as a member of the Technical Committee on Spectrum Estimation of the ASSP, is the Vice-chairman of and the Bay Area Chapter of that society. He is a member of the IEEE Admission and Advancement Committee. He is the recipient of the 1983 ASSP Senior Award for the paper Recursive L t i e Forms for Spectral Estimation, and atc the 1985 award for best paper of the year from the European Association for Signal Processingfor the paperOn the Computation of an Asymptotic Bound for Estimating Autoregressive Signals in White Noise. He is a member of Sigma Xi.

also a Lecturer at Stanford University, teaching graduate courses on linear systems,nonlineardetectionandestimation,andsystemidentification.In Novembcr of 1985 he joined Saxpy Computer Corporation, Sunnyvale, CA, asthe Director ofAdvancedTechnology. He is responsible for theR&D activities ofthecompanyand for thedevelopment of signalprocessing applications for the MATRJX I-a 1000-Mflop supercomputer. He has authored more than 170 publications in the areas of signal processing and estimation. His currentinterests include parallel and systolic processing architectures andtheirVLSIimplementations,advancedprocessingtechniques for underwater surveillance,high-resolution spectral analysis andarray processing, adaptive filtering, radar and communication processing, detection tracking and localization of multiple targets, image processing, and knowledge-based signal processing.