Track Parameter Estimation from Multipath Delay Information

Track Parameter Estimation from Multipath Delay Information (Invited Paper)

Abstract-It is desired to track the location of an underwater acoustic source with range differenee measurements from astationary passive array.Many times, the arrayhas only one or two sensors, andthe multipath and intersensor range differencemeasurements are insufficient to localize and track a source moving along an arbitrary path 111. Here, we propose to track sources with one- or two-sensor stationary passive arrays by making the simplifying assumption that thesource’s path canbe described by a small set of so-called track parameters. Range difference information can then be used to estimate the track parameter set rather than the source location as a function of time. In this paper, we choose the track parameters to specify a straight-line constant-velocity constant-depthpath.Cramer-Rao bounds are presented for estimating these track parameters from the time history of multipath and intersensor range difference measurements. It is shown that this track parameter set cannotbe accurately estimated from the time history of a single multipath range difference without side information (an independent velocity estimate, for instance), although multipath and intersensor range difference measurements from a two-sensor array are generally sufficient to estimate the track parameter set. Computationally efficient techniques are presented which estimate track parameters from range difference measurements taken from one- and two-sensor arrays. Monte-Carlo simulations are presented which show that these techniques have sample mean-square error approximately equal to the Cramer-Raobound when a single multipathrange difference and an independent velocity estimate are available. The sample mean-square error is shown to he in the range of two to ten times the corresponding Cramer-Rao bonnds when these techniques are applied to two-sensor range difference data.

I.

INTRODUCTION

T

HE tracking of a radiating source by a sensor array is a basic problem in underwater acoustics. In this paper, we describe methods for tracking sources using range difference measurements from passive stationary sensor arrays when the number of sensors is small and surface multipath reflections are present. The signals received by the sensors of an array exhibit relative time delays corresponding to a sources’ location. The relative time delay between two sensors, termed time difference of arrival (TDOA) is proportional to the source-sensor range difference (RD), and can be usedto locate a source on a hyperboloid of revolution about a line drawn through both sensors (a constant velocitymediumisassumed). In a p dimensional space, N > p independent RD’s are required to localize a source without additional information [ 11. AccordManuscript received May 12, 1986; revised August 25, 1986. J. S . Abel is with Systems Control Technology,Inc., Palo Alto, CA 94304. K. Lashkari was with Systems Control Technology, Inc., Palo Alto, CA 94304. He is now with Saxpy Computer Corporation, Sunnyvale, CA. IEEE Log Number 8714341.

ingly, a target moving along an arbitrary path can be tracked in p dimensions with N > p RD’s by smoothing source location estimates obtained at various times. Much work has been done on the problem of estimating TDOA (a special issue on the topic is [2] andrecent studies on the performance of time-delay estimators include [3]-[lo]). However, there is relatively little work on the problems of estimating source locations from TDOA measurements [11][15] and tracking moving sources using TDOA measurements from stationary passive arrays [16]-[19]. RD information from one- or two-sensor stationary passive arrays is insufficient to localize and therefore track a source moving along an arbitrary path. However, if it is assumed that the source travels along a path described by a parameter set, the problem of estimating a source’s location at many points in time (i.e., the problem of tracking a source) is reduced to the problem of estimating the parameter set describing the source’s path and may be possible with relatively little RD information. Many underwater sources of interest travel along straightline constant-depthconstant-velocity paths, specified by socalled track parameters-velocity, depth, bearing angle, and the radius and time ofclosest approach to the sensor array. In this paper, computationally efficient techniques for estimating a sources’ track parameters from RDmeasurements taken fromone- or two-sensor stationary passive arrays in the presence of a surface multipath reflection are presented and evaluated. The structure of this paper is as follows. Section II describes the track parameter estimation problem and reviews basic RD relationships. The Cramer-Rao lower bound on the variance of estimating the track parameters given RD and side information is derived in Section ILI. Track parameter estimation methods are developed in Section IV, while Section V reports computer simulation results and Section VI contains concluding remarks.

Iz. THETRACK PARAMETER ESTMATION PROBLEM We consider thefollowing scenario illustrated in Fig. 1. There is a stationary passive one- or two-sensor array listening to a target in the presence of a surface multipath reflection. The source moves by the array at constantvelocityand constant depth along a straight line. Range difference (RD) measurements (between the direct paths and/or multipaths) are available at various points during some interval of time. The

0364-9059/87/0100-0207$01.OO

0 1987 IEEE

.

.

. .

... .

208

-

. . .. .. . . .~

IEEE JOURNAL OF OCEAMC ENGINEEFSNG, VOL. OE-12, NO. 1, JANUARY 1987

ocean 1,11111’ 1111111

111111 111111 11111111

1111’1

3 /*

surface

11‘11111 Ill1 111

& ,

sensor

.e*

...**. +.+

@Jv&* -b

source Fig. 1. RD measurement scenario. A radiating source is shown moving along a constant-depth constant-velocity straight-line path. The sensor is shown receiving simals from the source both directly and from a surface multipath refleition.

track parameter estimation problem is to use the range difference measurements to estimate the parameters describing the target’s path. Below, track parameters are defined and RD relationships are described for three arrays: single sensor, two-sensor vertical, and two-sensor horizontal. Denote LS1and L s 2 as vectors containing the ( x , y , z ) coordinates of the sensor locations for sensors 1 and 2. Denote p as the vector of track parameters, and define

P Auy[ux

xT Y T Z T I T -

Side view

(1)

Here, u, and u, are the x-axis and y-gxis source velocities, is the s o ~ c depth, e and XT and y r determine the time and range of closest approach of the source to the sensor array. At any time ti, the source location may be given as

ZT

-ti

L T = [

8

0 13

-ti

0

]

p = [ x = - u x t i yr- uyti Z ~ I (2)

Top view

where I3 is the 3 x 3 identity matrix. The range difference between the source and the sensors at time tiis given by 4 2 ( 0

=d d l ( 0 - 4

2 U )

(3)

where x ( i ) denotes the ith component of the vector x, the vector d12contains the range differences for times t i , i = I , , N , and ddj(i)is given by

---

= ll M

i ) -L

j

II

(4)

A . Single Sensor In the case of a single sensor in the presence of a multiparh reflection, the track parameter vector, and the source and sensor locations can be given as (see Fig. 2) p = [ u XT

(3)

Y T ZTIT

&(i) = [xT- uti Y T z d L,= [O 0 z,] =.

(5)

Here the x-y plane of the coordinate system is the ocean

Fig. 2. Single-sensor RD measurement and source track. (a) The sensor at Ld located on the z axis and thesource at LT.The direct path and multipath lengths are shown as dd and dm.Note that the multipath le@ dmis the same as the direct path lengthfrom the source to a “virtual” sensor located at L,. (b) The source’s constant-velocitymotion parallel to the x axis and the sensor located at the origin of the x-y plane.

209

ABEL AND LASHKARI: TRACK PARAMETER ESTIMATION

Sideview

surface, positive depth is measured into the ocean, and the coordinate system origin is at the ocean surface above the sensor. Note, due to the radial symmetry of the array, there are only four elements of p . The RD between the direct path and multipathfor the sensor location above is equivalent to the RD between the direct paths to the sensor locations (see Fig. 2) :

Lsd = [0 0 z,] L,= [O 0 -z,lT

ocean surface

(true sensor) (virtual sensor).

The RD vector is given as

d=d,-dd

(6)

where the ith elements of dd and dm,the direct and multipath ranges, are given by

d m =

source

ll~sm-~T(~)ll

=[(XT-Ufi)2+y~+(z~+z,)2]1’2

dd(i)= IILsd-LT(i)ll = E(X7-

Uti)2+Y$.+(Z*-zs)2]1n.

(7)

B. Two-Sensor Vertical Array Here, the track parameter vector and the source and sensor locations can be given as (see Fig. 3)

P=[~xrYrzrlT LT(0 = [x=- uti y r 27-1 L,1= [O 0 zs11= Ls2=[0

zsZ>Zsl*

0

(8)

Again, due to the radial symmetry of the array, there are only four nonzero elements of p . The intersensor RD vector is given as

d12= & I -

dd2

(9)

where the ith elements of ddIand dd2, the direct path ranges, arc

0) Fig. 3. Vertical a m y RD measurement and source track. (a) The sensors at LSIand L,, on thez axis and the source at LT.The direct path lengths from the source to the sensors are shown as ddl and dd2;the multipaths are not shown. (b) The source’s constant-velocity motion parallel to thex axis and the sensors located at the origin of the x-y plane.

locations can be given as (see Fig. 4)

P = [ux

uy X T Y T Z T I T

L,(i) = [x=- uxti yT- urti ZT1 L 1 =

[x, 0 Z S I =

L,a = [ -x, 0 z,] T . Nc measured from a single sensor at (0, 0, 1/2(zS2- zsl))for the same setof target track parameters withp(4) given by zr - 1/ 2(2,2 + zSl)rather than ZT.

The intersensor RD vector is given as

d = ddl - dd2

(12)

where the ith elements of ddl and dd2, the direct path ranges, are given by

C. Two-Sensor Horizontal Array

With a horizontal array there is no longer radial symmetry, and the track Darameter vector andthe source and sensor

(11)

ddl(i)= 1lLs1-LT(i)ll

.

. .

.

.

.

210

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

Sideview ocean surface

estimating the track parameters p = [ u xT yT zT] from the time history ofa single multipath RD is typically large compared to the RD variance unless side information is available. It is also shown that side information may not be needed when a two-sensor array is used and intersensor RD information is available. Lower boundsonthe variance in estimating the track parameters p = [u xT yT zTITare presented for the case of known velocityanda single multipath FtD available at N points intime. Lower bounds on the variance in estimating the track parameters p = [ u xT y~ ZT] are presented for the case of multipath and intersensor RD's available from a twosensor vertical array at N points intime. Finally, lower bounds on the variance in estimating the track parameters p = [u, uy XT yT zT] are presented for the case of the multipath and intersensor RD's available from a two-sensor horizontal array at N points in time. The variance of an unbiased estimator of a parameter 8 is bounded below by the Cramer-Rao lower bound [20]-[231:

7 sensor 2

sensor 1

source (a)

Top view

( 8 ) = ~ ( 8 8 ~ ) - ~ ( 82) I;' ~(8)~

Lc.t,

(14)

where A L B indicates that A - B is positive semidefinite, and Io is the Fisher information matrix for the parameter 8 given by [20]

X

sensor 1

where f ( 0 Ix) is the conditional probability density function of

0 given x, and the estimate 8 is based on data x. If the distribution on 0 is parameterized by t (e.g., 8 = g ( t ) ) , then the Fisher information matrix for t can be written as [20]

Fig. 4. Horizontal array RD measurement and source track. (a)The sensors at L,, and Lr2at the same depth,and the source at LT.The direct paths from the source tothe sensors are shown andhavelengths d d , and dd2; the multipaths are not shown. (b) The source's constant-velocity motionand the Sensors located on the x axis.

Finally, if x and y are drawn from independent distributions, then [20] Ix,=Ix+I,. (17)

A . Single-Sensor Case dd2(i)=

IILS~-~T(~)\~

=[(XT-U,ti-X,)2+(yT-Uyfi)2f(zT-z,)2]112-

(13)

Note that the array is symmetric about the x axis, also the intersensor RD is insensitive to the sign of ( Z T - 2,).

III. BOUNDS ON

THE PERFORMANCE OF

From Section II and (14), the Cramer-Rao lower bound for estimates of p = [u XT Y T ZT] based on N single-sensor RD measurements is var ( f i ) 2 I;'

(18)

where

TRACK PARAMETER

ESTIMATORS

In this section, Cramer-Rao lower bounds on the variance of estimating track parameters from RD measurements are presented. It isshown here that the minimumvariancein

where d is avector composed ofRD observations at times ti, i = 1, -,N. The derivative, ad/ap is easily evaluated from [11, and is given by

--

.

21 1

ABEL AND LASHKARI: TRACK PARAMETER ESTIMATION

Side view

lines of constant RD

ocean surface

source 1

--d, (3)

Fig. 5. Track equivalence. Two sourcesmovingalong different constantvelocity constant-depth straight-linepathsare shown here to produce similar sets of single-sensor multipath RD measurements. Source 1 is shown moving along a constant-velocity constantdepth straight-line path producing RD measurements [d,d&] at times [tlt2t3].Source 2 is shown moving faster at greater depth alongaconstant-velocityconstantdepth straight-line path producing the same RD measuremens at the same times.

where kd and k, are constants depending on the track parameters and the sensor location. Using (23)in (21),we see that the four columns of ad/ap are approximately dependent on the three vectors

where

[

tl

IXT-ut'13 IxT-ut'13 IXT-ut,13

Assumingthe RD estimates are unbiased and Gaussian distributed with correlation matrix & the information in d is given by [20]

I . Ambiguity Issues: As shown below, the Cramer-Rao lower bound in(1 8) is typically very large compared to the RD estimate variances; therefore, low-variance unbiased track parameter estimates cannot be obtainedwithoutadditional information. For many sets of track parameters, during much of the observation time, the term ( X T - ut# will dominate dd(i)and dm(i),and l/dd(i)and l/dm(i)can be approximated by

--1

d,(i)

IxT- uti(

(1+

(xT- ut#

IxT-vfNI IxT-ufN13

; tN

*:]

.

(24)

1;

IxT-UfN13

Therefore, the matrix I~ = (ad/ap)TId(ad/ap)will be nearly singular and the resulting bounds on the minimum variance of track parameter estimators will be large compared to the variance in estimating the RD's. This result can be seen by examining the targedsensor geometry. As illustrated in Fig. 5, there are many sets of parameters which result in roughly the same measurements d . Away from the sensor array, the hyperboloids of constant RD can be approximated as cones. As a result, the differences in d between the case of a target moving slowly close by a sensor and the case of a target moving more quickly further from the sensor can be comparable to the deviation in estimating d . Consequently, there can be many sets of "equivalent" track parameters for a given d and the variance in estimating the track parameter set from measurements of d alone is high. 2. A Priori Estimate of u: Many times, an independent

estimate of the target velocity i s available, through Doppler measurements, for instance [24]. With an a prioriestimate of the target velocity (anunbiased Gaussian-distributed estimate with standard deviation of assumed) is the Cramer-Rao lower bound on the variance of estimating the track parameters is given by var ( p 2 [v^ f T 9 T

ZT]')

2 1;'

(25)

Id2,

and

Id12

are Id]

=Rill

Id2 = R Id12

=R

and the derivatives adl/ap and ad2/ap are given by (20) and ad12/apis

where

where 1 p . =---

1

'

where a d a p = [1 0 0 01 and adlap is given above. When ut is small enough, the columns of Ip are no longer linearly dependent and estimates of the remaining track parameters willhave variance comparable to the RD and velocity estimate variances (depending on the amount of RD data availableandthe track parameters, see Section V, simulations). As a result, the single-sensor track parameter estimation method described herein assumes an independent estimate of the velocity.

B. Two-Sensor Vertical Array When both intersensor RD multipath and RD measurements are available from a two-sensor vertical array, the CramerRao lower bound for estimating the set of track parameters p = [ V XT JJT given Z T ] is by

ddl(i) & 2 ( i ) ' Note here that, by the arguments above, the derivatives a d l / ap, ad2/ap, and ad12/apeach have approximately linearly dependent columns. However, unless the difference zs2 - zsl is smallcompared to the size of the array, the sum Ipl + Ip2 + Ip12will have linearly independent columns, and the CramerRao lower bound variance of the track parameters willbe comparable to the RD estimate variance (depending on the track parameters and the amount of data available, see Section V , simulations). C. Two-Sensor Horizontal Array When both intersensor RD and multipathRD measurements are available from a two-sensor horizontal array, the CramerRao lower bound for estimating the set of track parameters p = [vX vY XT Y T k T ] is given by var (P)2 1;' where,

Ip

is given as

IP'IPl+IP2+IP12

var ( p ) 2 1;' where, Ip is given below using (16) and (17): IP'IPl+IP2+IPlZ

where d l , d2, and d12are the two multipath and intersensor RD functions, Id,, Id2, and Idl2 are

where d l , d2,and d12are the two multipath and intersensor RD vectors (assumed statistically independent), the matrices Idl,

1

213

ABEL A h 9 LASHKARI: TRACK PARAMETERESTIMATION

and the derivatives ddl/dp, dd2/dp, and dd12/dpare given by

(32)

RD d , the sensor location Provided x, is large enough, we expect the sum Idl + Id2 + solution for fi based on the measured IdI2will notcontain linearly dependent terms and the Cramer- L,, and an a priori estimate of the velocity u. Rao lower bound on the track parameter estimate variance will 1. The Equation-Error Method: Recall, be comparable to the RD variance (depending on the track d=dm-dd (34) parameters and the amount of data available, see Section V, simulations). where IV. TRACK PARAMETER ESTIMATION

d,(i)=[(~~--ut~)~+y$+(z~+z,)~]~/~ (354

In this section closed-form equation-error methods are developed for estimating track parameters from a single multipath RD and a priori knowledge of source velocity. Thesemethods are extended to estimate track parameters using intersensor RD measurements from two-sensor vertical and horizontal arrays. Finally, methods are presented for estimating track parameters from intersensor andmultipath RD measurements from two-sensor vertical and horizontal arrays without a priori knowledge of the source velocity. A . Single-Sensor Track Parameter Estimation: The Equation-Error Approach

In the case of a single sensor in the presence of a multipath reflection, the functional form of the multipath RD is knownin terms of the track parameters. Therefore, a fit of a model d based on track parameters @ can be made to the measured d . The parameter estimates would then be chosen as $=min P

J(d, d )

(33)

for some cost function J . The minimization (33) is over a cost functionwhichis typicallynonconvexin p . Therefore, in- general, computationally expensive exhaustive search methods must be used in finding a solution of (33). In this case thecomputational burden can be overcome by choosing the cost function J so that functions of the track parameters p appear as coefficients in a linear least squares minimization, yielding a closed-form

~ ~ ( ~ ) = [ ( X T - U ~ ~ ) ~ + Y ~ T ~ ( Z T - Z , )(35b) ~]”~.

What is desired is an expression relating functions of d ( i ) to dk(i) and d;(i), so that functions of the unknown track parameters appear as linear coefficients whichcanbeestimated using least squares techniques. Manipulating (34) gives the following relationship: [d2(i)-d2,(i)-d2(i)I2-4d%i)d2(i)=0.

(36)

In (36) d,(i) and dd(i)are replaced with d,(i) and dd(i), the values predicted by the estimated track parameters @ = [ D f r j$ P r ] , and as the delays d are not known precisely, an equation error is introduced:

[d2(i)-d2,(i)-d:(i)l2-44d^2,(i)d^:(i)=ei.(37) We assumean a prioriestimate of the velocity (denotedby u ) is available, and assign v^ = u . Using (35) in (37)

d4(i)-d2(i)(z:+ t?2)+ 8d2(i)tiu{2,} -4d2(i)(22,+92,+22,}f 162:{22,} = c i (38) where the remaining track parameters to be estimated are [a, j r d ( i ) is the measured RD at time ti, and ~i is the corresponding equation error 1251. Defining E as the vector with ith element E ~ (38) , becomes

e,]

Sq-r=E

(39a)

... .

2 14

JEEE JOuRh'AL OF OCEANIC ENGINEERING, VOL. OE-12, NO. 1, JANUARY 1987

where

where

w(i) 4 2dm(i)d(i)(d(i)+d(i)). (39b) Note that the cost function J = E ~ isE made small by selecting track parameters which have 11 w 11 as small as possible such that d^ approximates d . The factor w ( i ) is made smaller by selecting track parameters which place the source closer to the (39c)point (0, 0, 0) (recall L, = [0 0 z J T ) .Thus it isexpected thatminimizing the unweighted equation error (minimizing J = E ~ E will ) produce track parameter estimates with a bias in - 8 d ' ( l ) t l4~d 2 ( 1 ) - 162: the direction [a, j T Z"T] = [0 0 01. (39d)Minimizing J = E 'E is equivalent to minimizing a weighted - 8d2(N)tNu4d2(N) - 162: L2 norm of the difference of the measured RD values and the estimated RD values. Specifically, the minimization J = E ~ Defining the cost function J as is equivalent to the minimization J P E~WE (40) J = E*T&* where W is a positive-definite weighting matrix, the set of track parameters minimizing J is given by where E* P d - d and A is a diagonal matrix with ith row/ column entry given as w2(i). If it is desired to minimize J* = E* T ~the , above suggests a weightingmatrixand an iterative procedure for doing so: solve (39) iteratively with weighting matrix

W,,

= diag

{ 1/w2(i)}

E

(45)

where g=(STWS)-'STWr. Clearly, the track parameter estimates givenin (41) are obtained with little computational expense compared to global search methods. The estimate $ given in (41) will be referred to as the equation-error estimate. The following sections discuss properties of the equation-error estimate and present methods for extending theequation-error estimate to the caseof intersensor RD data from two-sensor vertical and horizontal arrays. 2. The Equation Error: Here we study the cost function minimized in obtaining the equation-error estimate. Using equations (36) and (38) theith equation error can be writtenas

E(i)=[d2(i)-d2(i)+2dm(i)dd(i)12-4dK(i)d~(i) (42) where d is the measured RD value, d is the RD estimated by the track parameter estimates, and dd and dmare the direct and multipath ranges given by the track parameter estimates. Rearranging (42) gives

E(i)=[2d,,,(i)dd(i)+(d2(i)-d^2(i))](d2(i)-d2(i)). (43) The factor d2(i)- d2(i)is the difference between the square measured and estimated RD's, and should be much smaller than dm(i)dd(i),the productoftheestimated direct and multipath ranges. Therefore, the ith equation error canbe approximated as

€(i)-2d^~(i)dd(i)(d2(i) -d'(i)) = w ( i ) ( d ( i )- d ( i ) )

(44)

where diag ( x ( i ) } denotes a diagonal matrix with x ( i ) as the ith rowkolumnentry, and w ( i ) is computedusing track parameter estimates from the nth solution of (39) minimizing J = E T W , E . If the process converges, the error minimized is

Note that minimizing J* rather than J places more relative weight on theRD estimates taken when the source is closer to the sensor. The benefit of minimizing J* over J may be evaluated by studying the Fisher information matrix (25) for the relative importance of RD measurements as a function of source range. 3 . Variance of the Equation-ErrorEstimate: In this section we present expressions for the variance of the unweighted equation-error track parameter estimates when RD and a priori velocity information is noisy. a. Velocity known precisely: We first derive an expression for the equation-error estimate variance whenthe velocity isknown precisely andthe RD vector is unbiased with variance &. When the RD variance is small, the variance of the track parameters can be given by [26]

where ap/ad can be evaluated by the chain rule:

215

ABEL AND LASHKAlU: TRACK PARAMETER ESTIMATION

From (39), ap/aq is given as 0 P1

[1 P2

where

0

0 0

1 2P2

--1 2P2

0

1 -

0

(49)

2P3 e. Estimated RD's andvelocity: Finally, if the estimates of d and u are independent, then the variance of the equation error estimate j = [I?iTET can be given (50) approximately as the sum of the variances calculated above:

where pi is the ith element of p . Recall in the noiseless case:

Sq=r. Differentiating (50) by d gives

as -q + s

ad ad

aq ar aq

-=-

B. Two-Sensor Track ParameterEstimation In this sectionmethods for adapting the equation-error

where

".

-=diag ad

{4d3(i)+8d(i)(u2t:+z,2))

as

- q=diag { - 16d(i)tiuq1+8d(i)q2). ad

(53)

estimate to intersensor RD measurements from two-sensor arrays are presented. In addition, methods for estimating track parameters from two-sensor arrays without a priori information are developed. Track ParameterEstimation from Intersensor RD Estimates: The modification of the equation-error estimate (41) to make use of intersensor delay information rather than multipath delay information is straightforward. If the sensors are placed in a vertical array, the intersensor RD is functionally equivalent to the single-sensor multipathRD, and onlythe target and sensor depths need to be adjusted:

b. RD's known precisely: When the RD's are known precisely and the velocity estimate u is unbiased with small variance at, the equation-error estimate variance can be derived in a manner similar to (47)-(53):

1

where where ap/au can be evaluated by the chain rule: q=

I:[

q2 =(STWS)-lSTWr

where S and r areas defined in. (39) with d = d12 the intersensor RD estimate and W is a positive-definiteweighting matrix. In the case of a horizontal array, an equation-error estimate (56) can be derived in a manner similar to (34)-(41). From (12)

The derivative ap/aq is given in (49) and the derivative aq/au can be evaluated as follows. In the noiseless case: S q = r.

[d~,(i)-di,(i)-di2(i)]2=4d~l(i)di2(i)(62)

Differentiating (56) by u gives

as aq ar av q + s -=-

aU aU

where d12is the intersensor R D , ddl is the source range to sensor 1 and dd2is the source range to sensor 2. Substituting (13) into equation (62):

and solving for a q / h

Denoting the track parameters to be estimated by

8

=

.

.

. ..

.

~

. .

.

, ..:~ . ... ,

.

. .. . .

.

'

216

._

.

- -. .

. ._

. .

..

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

[Cx C,, 2T ET

25.1 and dI2as the measured intersensor R D ,

With data available from two sensors placed in a vertical array, sets of estimates, parameterized by u can be constructed using (41) and (61). A line search over u canthen be d ~ 2 ( i ) - 4 d ~ 2 ( i ) [ t ~ u 2 + x , 2 ] + 8 d ~ 2 ( i ) t ; ( C ~ , . 1 1 T + C ~ y 9performed T) for the set of track parameters minimizing some - 4 d : , ( i ) { 2 2 2 , i 9 2 , + ( 2 7 - ~ ~-4xsti{Gx} )~} +4xS{2,}= ~ i cost function. Emperically, it was found that minimizing the following cost function gave relatively unbiased low-variance (64) estimates for a large range ofnoise levels for the track parameter sets used in the Monte-Carlo simulations (see where an equation error E; has been introduced. The set of Section V): track parameters minimizing J = E W are E given by Jv = (2TlZ(U)- 2 T 2 ( u > ) 2 (66)

(63) becomes

where u

+ U;)I/~

= (u:

where 2 T 1 2 ( ~ )isthe depth estimate based on u using the intersensor RD information, and 2 ~ 2 ( u )is the depth estimate based on u using multipath RD information from the deeper sensor. If the data are from two sensors placed in ahorizontal array, sets of estimates, parameterized by u basedonindividual sensor data can be generated by (41). Estimatesof Qu parameterized by u can be made based on intersensor data (65) from (65). Track parameter estimates can then be chosen by is the a priori velocity estimate and finding theset of estimates producing the most consistent j r 2 ,GJu, i.e., the set, minimizing the cost function:

where Cx(u) is the x-axis velocity estimate obtained from the intersensor RD estimate using (65) and TI and gT2are the y axis range estimates from the multipathRD measurements and (41).

where

SA

V. SIMULATION RESULTS

1

d~2(l)-4d:2(l)(t~u2+x,2)

This section reports computer simulation results on the performance of equation-error estimates applied to noisy RD and velocity estimates. Sample bias, standard deviation, and RMS error, defined by

d ; J N ) - 4d;2(N)(r312+x,z)

and W is a positive-deffite weighting matrix. Note thatdue to symmetries in the array, the signs of y r and u,,, and (zT- z,) cannot be determined from RD information alone; however, the sign of Y T U T can be estimated, and here, Cy > 0 and gr takes the sign of jTG,,.

C. Combining Estimates from Two Sensors As shown in Section III, when RD information from more than one sensor is available, side information maynot be needed to obtain low-variance track parameter estimates. Below, methods for estimating track parameters from intersensor and multipathRD information from two-sensor vertical and horizontal arrays are given. These methods are not optimal anddo not produce minimum variance (ML) unbiased track parameter estimates, but provide acomputationally inexpensive alternative to the nonlinear optimization involved in computingthe ML track parameter estimates. If the resources are available to compute the ML estimate via an iterative nonlinear minimization, for example, these methods provide excellent starting points. '

sample bias P

1

h'

-

(&- 0)

Ni=, sample standard deviation A

sample RMS error

A

[N

1 N

-

1/2

(gi- 1 3 ) ~ ]

(68)

i=l

(where 8; is the ith estimate of the parameter 0) were calculated by averaging results of 100-trial Monte-Carlo runs and were compared to theoretically calculated values and Cramer-Rao lower bounds. Simulations were implemented in the Ctrl-C1 language on a VAX 11/785 computer. Results are given for sources moving along straight-line constant-velocity paths by one- and two-sensor arrays. Multipath and intersensor RD estimates were calculated by adding white Gaussian noiseto the true RD vectors; velocity estimates I

Ctrl-C is a registered trademark of SCT, Inc.

217

ABEL AND LASHKARI: TRACK PARAMETER ESTIMATION

TABLE I MONTE-CARLO RUNS Run

NoiseLevels (meters) Ov' Od

Track

Number

Number

1 2 3 4

1 1 1 1

0 0 0.1 1.0

10 100 0 0

5

1 1 1 1

0.1 1.0 0.1 1.0 0.1 1.0 0.1 1.0

10 100 100 10 10 100 100

6 7 8 9 10 11 12

2 2 2

13 14 15 16

1 1 2 2

2

0.1 1.0

0.1 1.0

io

-----

Track parameter set #1: [ v x v y x T yr zT]=[5 2 250 900 1701 Track parameter set #2: [ v x v y x T yr %I=[" 1 140 400 1101 SingleSensorLocation: VerticalArray 'Sensor Locations: HorizontalArray Sensor Locations:

(0,0,300) (O,O,lOO), (0,0,300) (200,0,300), (-200,0,300)

2000

1500

1000

Y (meters) 500

-3000

-2000

-1000

0

1000

I

I

2000

3000

X (meters) Fig. 6. simulation source tracks. The x-y plane source tracks used in the Monte-Carlo simulations are shown as dotted lines, each dot representing the source location atthe time of an RD measurement. The sensorlocations in the x-y plane are shown as squares; the vertical array and single-sensor arrayarerepresented by the centersquare andthehorizontalarrayis represented by the two outer squares.

(used in the one-sensor simulations) were calculated by adding white Gaussian noise to the true velocity. Table I describes the environmental data (sensor locations, track parameters noise levels, etc.) for each of the MonteCarlo runs (the source trajectories and sensor locations are

also shown in Fig. 6). Note that in the case of a single sensor or vertical array positioned at the origin, p = [u, vu XT Y T ZT] = [5 2 250 900 1701 is equivalent to p = [ u XT y~ ~7-= 1 ~[5.38 556.3 747.8 170IT; and [ - 3 1 140 400 llOIT translates to [3.162 8.849 423.7 1101T . Eachrunused100

1BEE JOURNAL 0F.OCEANIC ENGINEERING, VOL. OE-12,

218

NO. I , JANUARY 1987

TABLE II-A TRACK

PARAMETER EsnMATEi BIAS, SINGLE-SENSOR qRRAY

I

(meters) Bias

Rup

.

#

f;T%

b

1

-0.058

-0.090

2

-0.58

-0.93

3

-0.030

-0.084

4

-0.92

-1.8

r

=X

5

1.1

1.7

0.32

-0.009

6

3.5

5.5

0.70

10

-0.15

7

10

14

3.2

100

8

11

17

3.2

100

9

0.20

2.0

0.35

10

10

1.9

4.1

0.52

10

I1

0.36

20

3.5

100

3.5

100 -

=z

=v t

1.0

1.6

0.32

10

2

10

16

3.2

100

3

0.30

0.48

0.60

0

4

3.0

4.8

6.0

0

Run

I TheoreticalDeviation

(meters)

#

4

0,

02

=v t

1

1.0

1.6

0.32

10

2

10

16

3.2

100

3

0.32

0.46

0.60

0

4

3.2

6.0 4.6

0

TABLE III-B v

Run

I

TABLE II-D TRACK PARAMETER

TheoreticalDeviation(meters)

#

4r -

4

5

1.3

6

-0,

=vt

2.0

0.37

10

4.2

6.2

0.91

10

7

11

17

3.2

100

8

14

20

3.8

100

9

0.22

2.3

0.38

10

10

1.9

4.9

0.69

10

11

0.47

20

3.5

100

12

2.2

23

R M S ERROR,SINGLE-SENSOR ARRAY

-3.8

1OQ

TABLE IIIC TRACK PARAMETER CRAMER-R40 BOUND, SINGLE-SENSOR ARRAY

-

RMS Error (meters)

Run

2.0

10

TRACK PARAMETER THEORETICAL DEVIATION, SINGLE-SENSOR ARRAY

TABLE II-C .

21

12

TRACK PARAMETER THEORETICAL DEVIATION,SINGLE-SENSOR ARRAY

I

ovt

0,

1

.

-

CY

4

Y

RMS Error (meters)

Run OX

PARAMETER SAMPLE DEVIATION, SINGLE-SENSOR ARRAY

#

## -

TABLE II-8 TRACK

TABLE III-A

TRACK PARAMETER RMS ERROR, SINGLE-SENSOR W

Run

#

OX

0,

0,

=v t

1

1.0

1.6

0.32

10

CRB Deviation(meters)

# 4 5

-

OY

4

1.1

1.6

0.32

10

o vt

2

10

16

3.2

100

3

0.30

0.48

0.60

0

6

3.0

4.1

0.59

10

4

3.2

5.1

6.1

0

7

10

i6

3.2

99.8

8

11

16

3.2

1OD

9

0.17

2.0

0.35

10

10

1.7

3.3

0.47

10

11

0.36

20

3.4

97.2

1.7

20

3.5

100

points of data: one RD sample every 10 s from time tl = - 490 to tloo= 500 s, and the total observation time is t = tlm - tl = 990 s. So that all biases and deviations are measured in meters, the track parameter velocity is replaced here by the distance traveled during the observation interval u t .

A . Single Sensor Table II shows sample bias, swdard deviation, theoretical deviation (60), and R M S error of the unweighted equation-

12

219

ABEL AND LASHKARI: TRACK PARAMETER ESTIMATION

TABLE N - A TRACK PARAMETER

Run

1

ESIMATE BIAS AND STANDARD DEVIATION, TWO-SENSOR VERTICAL ARRAY

1

Bias (meters)

##

;,-x,

;,-y,

13

0.099

0.10

14

9

OY

0.67

0, t

4.2

0.60

5.5

0.97

34

0.1 0.018 0.22 0.0041 0.017

0.32

0.046

1.4

0.15 4.60.21

3.1

2.2 0.46

14

0.13

16

0.8

0.024 1.4

15

ox

(h)t

;,aT

Standard Deviation (meters)

5.5

6.3 47 4.5

0.12

TABLE N-B TRACK PARAMETER ESTIMATE RMS ERRORAND VERTICAL ARRAY

CRB, TWO-SENSOR

I I Run I #

1

RMS Error (meters)

ox

9

9

o vt

(meters) CRB

4

0,

02

0.31 0.68 0.053 0.34 4.3 0.13 0.61

13

0.53

o vt

2.2 22

14

7.1

8.4

1.7

58

3.4

3.1

15

0.21

0.32

0.046

.14

0.15

0.14

0.96 0.020

0.47 3.1 2.2 15

1.5

1.4

9.6 0.20

16

error estimate (41) applied to multipath RD estimatesand velocity estimates from a single sensor (Monte-Carlo runs 1 thru 4). Monte-Carlo runs 1 and 2 show the case of perfectly known RD estimates and noisy velocity estimates, and runs 3 and 4 showthe case ofnoisy RD estimates andperfectly known velocity. In these cases, the equation-error estimator is seen to be essentiallyunbiasedand to exhibit a standard deviation comparable to the RD and velocity deviations. Note that the theoretical deviations given by (60) agree very well with the sample RMS errors. Accordingly, it was noted that the track parameter estimate standard deviation appears to increase linearly with RD and velocity standard deviation. As expected (by the discussion of the equation error, Section IV-A), the sample bias of the equationerror estimate is negative in runs 3 and 4, indicating that the track parameter estimate pulls the estimated source location closer to the origin than the true source location. Table III shows sample RMS errors, theoretical deviations (60), and Cramer-Rao bounds(25)of the equation-error estimator applied to RD and velocity estimates from a single sensor (Monte-Carlo runs 5-12). The sample bias, not shown here, was noted to be small compared to the sample standard deviation. The R M S error appears to be slightly smaller than theoretically predicted, consistent with a possible small cross correlation between the RD and velocity noises seen in these relatively small sample size Monte-Carlo runs. The sample RMS errors (and theoretical deviations) appear to be very close to the Cramer-Rao bounds indicating that the equationerror estimate is using the RD and velocity information in an efficient manner for these cases.

B. Two-Sensor Vertical Array Table IV shows sample bias, sample standard deviation, sample RMS error, and Cramer-Rao bounds (27) for equationerror estimates applied to RD measurements from a twosensor vertical array; a priori velocity estimates were not given. The track parameter estimates were made from RD measurements taken from the deepest sensor withvelocity estimates determined by (66). The equationerror estimate sample bias is seen to be a strong (perhaps quadratic) function of the RD noise deviation, and is comparable to the sample standard deviation at the larger value of additive RD noise. The sample standard deviation appears to be a linear function of the RD noise deviation, and small compared to the track parameter values. The equationerror estimate sample R M S error is seen to be about twice the Cramer-Rao bound.

C. Two-Sensor Horizontal Array Table V shows sample bias, sample standard deviation, sample R M S error, and Cramer-Rao bounds for the case of the source moving past a horizontal array. The track parameters were estimated by averaging the translated track parameter estimates given by (41) and (35) with a velocity estimated by minimizing the cost function given in (67). Here, the track parameter estimates are essentially unbiasedandhave a sample standard deviation whichissmallcompared to the track parameter value andappears to increase linearly with RD standard deviation. The estimates have an RMSerror of about five to ten times the Cramer-Rao bound, depending on the track parameters. Here, it appears thatthe track parameter estimation method is not using the RD information efficiently.

220

IEEE JOURNAL OF OCEANIC ENGINEERING, VOL. OE-12, NO. 1, JANUARY1987 TABLE V-A TRACK PARAMETERESTIMATE BIAS, TWO-SENSOR ARRAY

Bias (meters)

I

A

XT-X T -0.11

0.49 -0.49 0.26

TABLE V-B TRACK PARAMETER STANDARDDEVIATION,TWO-SENSORHORIZONTAL ARRAY

1

Run

I

Standard Deviation

(meters)

1

TABLE V-C TRACK PARAMETER RMS ERROR, TWO-SENSOR HORIZONTAL ARRAY

Run

RMS Error (meters)

moving along an arbitrary path, and here, the problem of describing the source’s location as afunctionof time was reduced to the problem ofestimating a small set of parameters describing an assumedstraight-line constant-velocity constantdepth source path. Cramer-Rao bounds were presented for estimating the track parameters from the time history of multipath andintersensor range difference measurements. It is shownthat this track parameter set could not be accurately estimated from the time history of a single multipath range difference without side information. However, multipath andintersensor range difference measurements from a two-sensor array were seen to be sufficient to estimate the track parameter set when the sensors are appropriately placed and enough RD data are available. Linear least squares equationerror techniques were presentedwhich estimate track parameters from independent velocity estimates andmultipath range-difference measurementstaken from a one-sensor array. Analytic expressions were developed for the variance of the estimate, and MonteCarlo simulations were presented whichshowthat these estimators are relatively unbiased and have sample RMS error which is given accurately by the analytic expressions and is approximately equal to the Crber-Rao bound. Line-search methods were developed to estimate track parameters from two-sensor RD data whennoindependent velocity estimates were available. The samplemean-square error of the track parameter estimates producedby these methods was shown by Monte-Carlo simulations to be in the rangeoftwo to tentimes the corresponding Cramer-Rao bounds. These methods provide a computationallyinexpensive alternative to the nonlinear optimization required in finding the ML estimate. ACKNOWLEDGMENT The authors would like to thank T. Furukawa for help in preparing the manuscript.

TABLE V-D TRACK PARAMETER CRB,TWO-SENSOR HORIZONTAL ARRAY

CRB Deviation ( m e t e r s )

I

However, depending on the application this estimator might have acceptable performance.

VI. SUMMARY

En this paper, we discussed the problem of tracking sources with multipath andintersensor range difference measurements from one- or two-sensor stationary passive arrays. RD data gathered from such arrays are insufficient to track a source

REFERENCES

J. P. Van Etten, “Navigation systems: Fundamentals of low and very low frequency hyperbolic techniques,” Elec. Commun., vol. 45, no. 3, pp. 192-212,1970.

IEEE Trans. Acoust., Speech, Signal Processing (Special I m e on Time-Delay Estimation), vol. ASSP-29, June 1981. J. P. Ianniello, “The variance of multipath time delay estimation using autocorrelation.” Naval Underwater Systems Center, New London, CT, Tech. Memo 831008, Jan. 1983. J . P. Ianniello, “Large and small error performance limitsfor multipath time delay estimation,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-34, no. 2 , pp. 245-251, Apr. 1986. J. P. Ianniello,“Timedelay estimation via cross-correlation in the presence of large estimation errors,” ZEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-30, pp. 998-1003, 1982. A. Weiss and E. Weinstein, “Fundamental limitations in passive time delay estimation-Part 1: Narrow-band systems,” ZEEE Trans. Acoust., Speech,Signal Processing, vol.ASSP-31, pp. 472-485, 1983. J. P. Ianniello, E. Weinstein, and A. Weiss, “Comparison of the Zivlower bound on time delay estimation with correlator performance,” in Proc. ZCASSP 83 Conf. (Boston, MA), Apr. 1983, pp. 875-878. J. P. Ianniello, E. Weinstein, and A. Weiss, “Lower bounds on worstcase probability of large error for two channel time d e w estimation,” submitted to ZEEE Trans. Acoust., Speech, Signal Processing. B. Friedlander, “OntheCramer-Raobound for time delay and

=ai

22 1

ABEL AND LASHKARI: TRACK PARAMETER ESTIMATION

Doppler estimation,” ZEEE Trans. Znform. Theory, vol. IT-30, no. 3, pp. 575-580, May 1984. [lo] J. C. Hassab and R. Boucher, “A probabilistic analysis of time delay extraction by the cepstrumin stationary Gaussiannoise,” ZEEE Trans. Inform. Theory, vol. IT-22, pp. 444-454, July 1976. [ l l ] R. 0.Schmidt, “A new approach to geometry of range difference location,” ZEEE Trans. Aerosp. Elec. Syst., vol. Am-8, no. 6, pp. 821-835, NOV.1972. [I21 J. M.Delosme, M. Morf, and B. Friedlander, “A linear equation approach to locating sources from timedifference-of-arrival measurements,” in Proc. ZEEEZnt. Conf. Acoust., Speech,andSignal Processing, 1980. [13] H. C.Schau and A. Z . Robinson, “Passive source localization employing intersecting spherical surfaces from time-of-arrival differences,” ZEEE Trans. Acoust., Speech, Signal Processing, accepted for publication. [I41 J. 0.Smith and J . S. Abel, “Closed-form least-squares localization of multiplebroad-band emitters,” SystemsControlTechnology, Inc., Palo Alto, CA, Tech. Memo. 5517-01, Feb. 1986. [I51 J. S. Abeland J. 0. Smith, “Ontheefficiency of thespherical interpolationmethod for closed-form source localization,” Systems Control Technology, Inc., Palo Alto, CA, Tech. Memo. 5517-02, June 1986. [16] R. L. Moose and T. E. Dailey, “Adaptive underwater target tracking usingpassivemultipathtime-delay measurements,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-33, no. 4, pp. 777778, Aug. 1985. [17] D. H.McCabe and R. L. Moose, “Passive source tracking using sonar time delay data,” ZEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 614-617, June 1981. [IS] R. L. Moose, “Passive source tracking using differential timedelays,” in Proc. IEEE 1981 EASCON Conf. (Washington, DC), Nov. 1981. [19] J . C. Hassab, “Passive tracking of a moving source by a single observer in shallow water,” J. Sound and Vibration, vol. 44, no. 1, 1976. [20] E. L.Lehmann, Theory of Point Estimation. New York:Wiley, 1983,pp.123-129. [21] R. A. Fisher, “On the mathematical foundations of theoretical statistics,” Phil. Trans. Roy. SOC.(London), vol. 222, p. 306, 1922. [22]H. Cramer, Mathematical Methods of Statistics. Princeton,NJ: Princeton Univ. Press, 1935. [23] C. R. Rao, “Information and accuracy attainable in the estimation of statistical parameters,” Bull. Calcutta Math. Soc., vol. 37, pp. 8191, 1945. [24] K. Lashkari, B. Friedlander, and J . S. Abel, “Track parameter

estimation using multipath delay and Doppler information,” Systems Control Technology, Inc., Palo Alto, CA, Final Rep. 5517, Feb. 1986. L. Ljungand T. Soderstrom, Theory and Practice of Recursive Identification. Cambridge, MA: MIT Press, 1984. D. R. Coxand D. V. Hinkley, Theoretical Statistics. London, England: Chapman and Hall, 1982.

*

Jonathan S. Abel (”85)

wasbornin Sarasota, FL, in December 1960. He received the B.S. degree in electrical engineering from the Massachusetts ’ Institute ofTechnology, Cambridge, MA, in 1982, and won the 1982 Guilleman Award for his thesis workonconductingpolymers.Hereceivedthe M.S. degree inelectrical engineering from Stanford University, Stanford, CA, in 1984,where he is currently a Ph.D. candidate in theInformation Systems Laboratory. He joined SystemsControlTechnology,Palo Alto,CA, in 1985 as a ResearchEngineerintheAdvancedTechnology Division, where he has been working in the areas of signal processing and information theory. ”

*

Khosrow Lashkari (S’78-M’82) was born in Tehran, Iran, on February 8, 1949. Hereceived the

B.S. degree in electrical engineering from Aryamehr University of Technology in 1972 andthe M.S.and Ph.D. degreesin electrical engineering fromStanfordUniversity, Stanford, CA, in 1977 and 1982, respectively. From 1971 to 1973he worked at Digital Systems Lab., Tehran, Iran, where he was involved in the design and prototype developmentof a voice scrambling system. From1973to1975 hewaswith Electronic SystemsCenter, Tehran, Iran, where he was involved in the design and development of data transmission systems. From 1982 to 1986he was at Systems Control Technology, Inc., Palo Alto, CA. He was responsible for a number of research and development projectsin speech compression, speech recognition, image processingand sonar signal processing. In August 1986 he joined Saxpy ComputerCorporation, Sunnyvale, CA, wherehe is involved in definingsignalprocessingapplications for the SAXPY-IM-a 1000-Mflop supercomputer.