A decoupled algorithm for geolocation of multiple emitters

ARTICLE IN PRESS

Signal Processing 87 (2007) 2348–2359 www.elsevier.com/locate/sigpro

A decoupled algorithm for geolocation of multiple emitters$ Alon Amar, Anthony J. Weiss School of Electrical Engineering—Systems Department, Tel Aviv University, Tel Aviv 69978, Israel Received 1 August 2006; received in revised form 22 March 2007; accepted 22 March 2007 Available online 30 March 2007

Abstract The problem of determining the positions of multiple emitters by widely separated arrays is considered. We propose an iterative algorithm for estimating the positions based on the maximum likelihood (ML) criterion. The position of each emitter is decoupled from the other emitters, and in each iteration is determined by a two dimensional grid search. This approach substantially reduces the computation load of a straightforward implementation of the ML estimator, which requires a multidimensional search. Numerical examples are provided to illustrate the performance. We demonstrate that only a few iterations are required to achieve convergence. The results are compared with the Crame´r–Rao lower bound. r 2007 Elsevier B.V. All rights reserved. Keywords: Array-processing; Emitter-localization; Maximum-likelihood estimation

1. Introduction The position estimation of multiple emitters, intercepted by several receiving stations, is usually performed by a sub-optimal two step procedure. Recently, several papers discussed optimal estimation methods based on exploiting the raw data intercepted by all stations in a fusion center. The recent increase of cellular communications generated a renewed interest in this field. Cellular providers may use the location of subscribers for a

$

This research was supported by the Israel Science Foundation (Grant No. 1232/04) and by the Institute for Future Technologies Research Named for the Medvadi, Shwartzman and Gensler. Corresponding author. Tel.: +972 36408605; fax: +972 36407095. E-mail addresses: [email protected] (A. Amar), [email protected] (A.J. Weiss).

variety of applications, such as commercial advertisement, self navigation, rescue requests, etc. Most publications on decentralized localization methods (two-steps methods) concentrate on the first step, i.e., Angle-of-Arrival (AoA) estimation, or Time-of-Arrival (ToA) estimation, performed by each station independently [1]. Several mathematical algorithms were proposed for performing the second step: position determination based on triangulation using the results of the first step. Details can be found in [2], and the reference therein. Wax and Kailath [3] proposed an intermediate approach between decentralized and centralized (single-step methods) geolocation. Instead of transmitting the raw data to the fusion center, as is done in the centralized approach, each station transmits only the sufficient statistic (the sample correlation matrix). The geolocation is determined using a cost function based on eigendecomposition of the complete set of sufficient statistics.

0165-1684/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2007.03.008

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

Surprisingly, the number of publications on centralized geolocation is rather small. Chen et al. [4] discussed Maximum-Likelihood (ML) geolocation estimation for near-field emitters. Kozick and Sadler [5] examined the geolocation of a single aeroacoustical emitter, focusing on the correlation between the signals intercepted by different stations. Recently, the authors presented several geolocation algorithms for unknown and known signals [6]. The analysis and simulations showed that the proposed method (Direct Position Determination (DPD)) outperforms conventional, two-steps methods, based on AoA or ToA. The objective of the current correspondence is to present a Decoupled Geolocation Maximum Likelihood (DGML) algorithm for centralized geolocation. Up to date, a centralized ML position determination of multiple emitters with unknown waveforms has not been presented. A straightforward solution of the ML requires a multidimensional search with high complexity. We resort to an iterative least squares algorithm to overcome this difficulty. The proposed algorithm estimates the positions of Q emitters in a plane by a two-dimensional search in each iteration, instead of a single 2Q multidimensional search. It is shown that convergence to the Crame´r–Rao Bound (CRB) is achieved with a small number of iteration steps. It is worth mention that the algorithm, as any other iterative algorithms, depends on the initial estimate of the emitters’ positions. However, the algorithm can be applied to tracking of multiple emitters, where there is an initial estimate of the positions. Such initial estimate can be provided using the centralized algorithms presented in [6]. 2. Problem formulation Consider L stations located at geographically separated sites within a plane. Each station is equipped with an antenna array. Each array consists of M sensors. One of the sensors in each array is used as a reference sensor and its location is denoted by d‘ 9ðX ‘ ; Y ‘ Þ. Assume that the Q sources are located in the far-field of the arrays. This implies that the wave-front can be considered planar w.r.t. the array aperture. However, since the stations are widely separated, the wave-front curvature should be taken into account when considering all the arrays together. Denote by pq 9ðxq ; yq Þ, q 2 f1; . . . ; Qg, the position of the qth source. The propagation time from the qth source to the reference sensor of the ‘th array is given by

2349

tq;‘ 9ð1=cÞkpq  d‘ k, where c is the wave propagation speed. Let the M  1 column vector, aq;‘ denote the ‘th array response to the qth source. The complex envelopes of the signals observed by the ‘th array are given by r‘ ðtÞ ¼

Q X

bq;‘ aq;‘ sq ðt  tq;‘ Þ þ n‘ ðtÞ;

0ptpT,

q¼1

(1) where T is the observation time, bq;‘ is an unknown non-zero complex scalar representing the attenuation of the qth signal as observed by the ‘th station, sq ðÞ is the qth signal waveform, and n‘ ðÞ represents the noise and interference observed by the array. We assume that: (1) The signals are realizations of a continuoustime, ergodic, zero-mean, wide-sense stationary, complex Gaussian random processes. (2) The noise at the output of each of the array elements is wide-sense stationary, zero mean, complex Gaussian process, uncorrelated with the noise at the other antenna elements and uncorrelated with the signals. (3) The signals and noise processes have the same bandwidth, W, and the same center frequency o0 , (both in rad/s). (4) The length of the observation interval T is long compared with the correlation time of the signals and the correlation time of the noise processes, and also long compared to the propagation time between the stations. (5) Without loss of generality the attenuations to the reference (first) station assumed real, i.e., Pare L 2 Imfbq;1 g ¼ 0, 8q and jb ‘¼1 q;‘ j ¼ 1. These constraints eliminate the trivial solution. The first stage of the processing consists of partitioning the observation interval to K sections, each of length T 0 9T=K. Since the observed signal is a stationary process, it can be represented by a set of Fourier coefficients Z kT 0 1 r‘ ðk; jÞ9 pffiffiffiffi r‘ ðtÞeioj t dt, (2) T ðk1ÞT 0 where k ¼ 1; . . . ; K, and oj 9ð2p=T 0 Þj, j ¼ 0; 1; . . . ; J. The Fourier coefficients of (1) satisfy r‘ ðk; jÞ ¼

Q X q¼1

bq;‘ aq;‘ sq ðk; jÞeioj tq;‘ þ n‘ ðk; jÞ,

(3)

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

2350

where r‘ ðk; jÞ, sq ðk; jÞ, and n‘ ðk; jÞ are the jth Fourier coefficients of the kth section of r‘ ðtÞ, sq ðtÞ, and n‘ ðtÞ. The noise vectors fn‘ ðk; jÞg are i.i.d., zero-mean complex Gaussian vectors with covariance s2 I. Arranging the observations from all arrays in a column vector yields

3. Decoupled geolocation ML estimation

rðk; jÞ ¼ AðjÞBsðk; jÞ þ nðk; jÞ,

where TðlÞ is the cost function

(4)

where rðk; jÞ9½rT1 ðk; jÞ; . . . ; rTL ðk; jÞT , nðk; jÞ9½nT1 ðk; jÞ; . . . ; nTL ðk; jÞT , sðk; jÞ9½s1 ðk; jÞ; . . . ; sQ ðk; jÞT ,

¼ ð5Þ

B9Diagfb1 ; . . . ; bQ g, ð6Þ

DiagfX1 ; . . . ; XR g, denotes a HR  TR block-diagonal matrix where Xn is a H  T matrix. Collecting all frequency coefficients yields (7)

where IJ is the J  J identity matrix, and rðkÞ9½rT ðk; 1Þ; . . . ; rT ðk; JÞT , s9½sT ðk; 1Þ; . . . ; sT ðk; JÞT , n9½nT ðk; 1Þ; . . . ; nT ðk; JÞT , A9DiagfAð1Þ; . . . ; AðJÞg.

ð8Þ

This concludes the basic mathematical formulation of the problem. The problem considered in this paper may be stated briefly: Given a collection of samples rð1Þ; . . . ; rðKÞ and the model in (7), estimate the parameters l9½sT ; bT ; pT T ,

(9)

where s9½¯sT ð1Þ; . . . s¯T ðKÞ; s~T ð1Þ; . . . s~T ðKÞT , T

T

T

T

b9½b¯ 1 ; . . . ; b¯ Q ; . . . ; b~ 1 ; . . . ; b~ Q T , p9½pT1 ; . . . ; pTQ T ,

(11)

l

K X

krðkÞ  AðIJ  BÞsðkÞk2

k¼1

AðjÞ9½W1 ðjÞ; . . . ; WQ ðjÞ, Wq ðjÞ9Diagfaq;1 ðjÞ; . . . ; aq;L ðjÞg,

rðkÞ ¼ AðIJ  BÞsðkÞ þ nðkÞ,

l^ ¼ arg minfTðlÞg,

TðlÞ ¼

and

bq 9½bq;1 ; bq;2 ; . . . ; bq;L T .

The ML estimator of the parameter vector given the observation model in (7) is

ð10Þ

and x; ¯ x~ represent the real part and the imaginary part of x, respectively. We are mainly interested in estimating the position vector, p, while the vectors, s, and b, are nuisance parameters.

K X J X

krðk; jÞ  AðjÞBsðk; jÞk2 .

ð12Þ

k¼1 j¼1

The major computational challenge is that minimizing (11) requires a multidimensional search over the parameter space. To reduce the search dimension we adopt, and modify the technique presented in [7], to the geolocation model with wideband signals. Assume we have an initial rough estimate of the vectors fpq ; bq gQ q¼1 . As in many iterative algorithms, the convergence to the global minimum of the cost function depends on the initial estimate (such an initial estimate can be provided by the DPD method based on the MUSIC approach presented in [6]). Also, note that the initial estimate of the positions can be used to relax Assumption 4 in Section 2. The original requirement on the DFT window (T=K) requires that T=Kb maxftmax ; Dtmax g, where tmax is the maximal correlation time of the signals and noise and Dtmax is the propagation time between the two most spatially separated stations. However, while the restriction T=Kbtmax is required to obtain uncorrelated frequency samples, the restriction T=KbDtmax ensures that the DFT edge effect is negligible. The edge effect can be eliminated by adjusting the observation interval in all stations according to the assumed propagation time between the emitter and the stations using the available (initially rough) emitter location. As the algorithm progresses and the emitter location becomes precise the edge effect can be fully eliminated. We next compute the matrices AðjÞ, and B in (6). Now, we can minimize (11) w.r.t. sðkÞ. The result is s^ ðkÞ ¼ ½^sT ðk; 1Þ; . . . ; s^T ðk; JÞT , ^ BÞ ^ þ rðk; jÞ, s^ ðk; jÞ ¼ ðAðjÞ

ð13Þ

where Xþ 9ðXH XÞ1 XH is the pseudo-inverse ma^ ^ are the matrices AðjÞ, B trix of X , and AðjÞ, B computed with the initial estimate of fpq ; bq gQ q¼1 .

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

Define the Q  Q diagonal matrix, Eq , where the entries on the main diagonal are all ones, except for the ðq; qÞth entry, which is zero. Define the vector: s^ðqÞ ðk; jÞ9Eq s^ðk; jÞ.

(14)

Also define a vector, associated with the qth emitter, ^ B^ ^ sðqÞ ðk; jÞ. wq ðk; jÞ9rðk; jÞ  AðjÞ

(15)

Assume that the positions and attenuations of all the emitters, except those of the qth emitter, are equal to their true values. In such a case, using the definitions in (6), Eq. (15) reduces to wq ðk; jÞ ¼ Wq ðjÞbq sq ðk; jÞ þ nðk; jÞ, that is, the qth column vector AðjÞ with the additive noise. This motivates a decoupled search. Instead of minimizing (11), we minimize, F¼

K X J X

kwq ðk; jÞ  Wq ðjÞbq s^q ðk; jÞk2 ,

(16)

k¼1 j¼1

where s^q ðk; jÞ is the qth element of s^ðk; jÞ in (13). The minimization of (16) w.r.t. s^q ðk; jÞ is given by s^q ðk; jÞ ¼ ðWq ðjÞbq Þþ wq ðk; jÞ.

(17)

Substituting (17) in (16), it is straightforward to show that the estimate of pq , and bq is given by ( ) K X J X 2 H H 2 ^ ð^pq ; bq Þ ¼ arg max kbq k jbq Wq ðjÞwq ðk; jÞj . pq ;bq

k¼1 j¼1

(18) According to the assumptions kbq k ¼ 1 in (18). Define the sample correlation matrix associated with the qth signal and the jth Fourier coefficient, K X b q ðjÞ9 1 wq ðk; jÞwH W q ðk; jÞ, K k¼1

Dq 9

that is: b^ q ¼ umax ðDq Þ.

(22)

Substituting this choice of bq in (21) we get p^ q ¼ arg max lmax fDq g,

(23)

pq

where lmax fDq g is the largest eigenvalue of Dq . Therefore, the maximum of the cost function is found by performing a two dimensional grid search over the area of interest. For each grid point the largest eigenvalue of Dq is the local value of the cost function. The grid point associated with the largest value of the cost function is the temporary estimate of the qth source position. The DGML procedure is therefore given by: Initialization stage: ð0Þ Q (1) Select an initial set of the vectors fpð0Þ q ; bq gq¼1 . (2) Define a convergence threshold 40. Reset the iteration counter to m ¼ 1. (3) Compute the initial matrices Að0Þ ðjÞ; Bð0Þ , as given in (6). Calculation stage: For q ¼ 1 to Q, repeat steps 4–9: (4) Compute s^ðqÞ ðk; jÞ as expressed in (14). (5) Compute the vector wq ðk; jÞ as expressed in b q ðjÞ in (19). (15), and the matrix W (6) For any position ðx; yÞ in the search grid: (6.1) Compute the matrix Wq ðjÞ expressed in (6). (6.2) Evaluate the matrix Dq expressed in (20). (6.3) Calculate the maximal eigenvalue, and the corresponding eigenvector of Dq . ðmÞ (7) Determine the vectors p^ ðmÞ , b^ associated with

(19)

the grid point having the largest cost function value. (8) Update the matrices AðmÞ ðjÞ, and BðmÞ , ðmÞ expressed in (6), using p^ ðmÞ , b^ . k

b WH q ðjÞWq ðjÞWq ðjÞ.

(20)

j¼1

Substituting (20) in (18) we get ð^pq ; b^ q Þ ¼ arg maxfbH q Dq bq g.

q

k

and the L  L Hermitian matrix J X

2351

(21)

pq ;bq

The right side of (21) should be maximized w.r.t. the attenuation vector and w.r.t. the emitter position. For any given emitter position, the vector bq , that maximizes the expression is equal to the eigenvector, umax , corresponding to the largest eigenvalue of Dq ,

q

Check convergence: P ^ qðm1Þ kp. If no, repeat steps (9) Is Q pðmÞ q p q¼1 k^ 4–8. If yes, stop. This summarizes the proposed algorithm for the geolocation of multiple emitters. Notice that at each updating steps (that is step 4, and step 7) the ML cost function, TðlÞ, in (12) is decreased. Since this cost function is always greater or equal to zero, the proposed algorithm will converge at least to a local minimum of TðlÞ.

ARTICLE IN PRESS 2352

A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

Obviously, good initial estimates improve the probability of converging to a global minimum. Since the algorithm is based on the ML criterion, a good predictor of the estimation performance for large number of snapshots (i.e., large K) is the CRB. The CRB provides a lower bound on the covariance matrix of any unbiased estimator. In fact, under mild regularity conditions, the ML estimator achieves the CRB asymptotically, as the number of snapshots tends to infinity. In Appendix A we present the derivation of the CRB for the model in hand. 4. Numerical examples In this section we demonstrate the performance of the DGML via Monte Carlo simulations. In all simulations (except otherwise stated) we used the following parameters:









 

Four base-stations (L ¼ 4), located at the corners of a square at positions: ðX 1 ; Y 1 Þ ¼ ½0; 0, ðX 2 ; Y 2 Þ ¼ ½1; 0, ðX 3 ; Y 3 Þ ¼ ½1; 1, and ðX 4 ; Y 4 Þ ¼ ½0; 1, all units are kilometers. A uniform linear array with three sensors (M ¼ 3) is located at each station. The spacing between adjacent sensors is D ¼ l=2, where l is the wavelength of the transmitted signals. The rotation angles of the arrays, measured between the array baseline and the west–east axis, denoted by f‘ , are: f1 ¼ 135 , f2 ¼ 45 , f3 ¼ 135 , and f4 ¼ 45 . The signals, and the noise, are random, circular, complex Gaussian vectors with covariance matrices s2s I, and s2 I, respectively. We assume the signals are narrowband, hence, J ¼ 1. We generated 50 independent Fourier coefficients ðK ¼ 50Þ of the signals and noise processes. Our definition of signal to noise ratio (SNR) is SNR½dB910 log10 ðs2s =s2 Þ. The initial emitters positions were set randomly according to: pð0Þ q ¼ pq þ 100  nq , q ¼ 1; . . . ; Q where nq is a 2  1, Gaussian vector with zeromean and covariance equal to the identity matrix. The initial values of the attenuations were chosen randomly subject to the restrictions described in the assumptions. The grid search was performed in two steps: a coarse search with 10 (m) resolution over the full area of interest and a fine search with 1 (m) resolution near the peaks found in the course search.



The convergence threshold of the algorithm was set to  ¼ 1 (m).

In the following examples we consider several values of SNR. For each SNR we perform 50 independent trials. In each trial the initial positions, and initial attenuations are randomly chosen. We then estimate the positions and attenuations of the emitters in each iteration step (note that the number of iteration steps in each trial is not necessarily equal). The Root Mean Square Error (RMSE) for each iteration step is defined by N qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1X RMSEq ðiÞ9 ðx^ n;q ðiÞ  xq Þ2 þ ðy^ n;q ðiÞ  yq Þ2 , N n¼1 i ¼ 1; . . . ; I, where ðx^ n;q ðiÞ; y^ n;q ðiÞÞ is the estimated position of the qth emitter at the nth experiment and at the ith iteration step, I is the maximum number of iterations in all trials, and N is the number of experiments with at least i iterations. Denote by N exp the number of experiments. As indicated previously N exp ¼ 50. The value of N=N exp in percent is shown in the plots. 4.1. Geolocation of a single emitter A single emitter ðQ ¼ 1Þ is positioned at p1 ¼ ½304:52; 699:57 (m). In Fig. 1 we plot the RMSE versus the number of iterations, for two values of SNR: 5 (dB) (left pane), and 10 (dB) (right pane). Also shown is the CRB. As can be seen, the algorithm converges in one iteration. For the high SNR the estimation error approaches the CRB. 4.2. Geolocation of two emitters: linear arrays Assume two emitters ðQ ¼ 2Þ are positioned at p1 ¼ ½202:79; 596:29 (m), and p2 ¼ ½801:26; 499:30 (m). In Fig. 2 we show the RMSE of the first emitter for two values of SNR: 5 (dB) (left pane), and 20 (dB) (right pane). The CRB is also shown. The numbers indicate N=N exp in percent. As can be seen, in both cases, after only two iterations the CRB is reached. 4.3. Geolocation of two emitters: circular arrays Assume two emitters ðQ ¼ 2Þ are positioned at p1 ¼ ½202:78; 603:27 (m), and p2 ¼ ½803:50; 496:25

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

SNR = 5[dB]

SNR = 10[dB]

160

120

120

100

100

80

Initial

DGML CRB

140

100%

RMSE [m]

RMSE [m]

160

DGML CRB

Initial 140

2353

80

60

60

40

40

20

20

0

100%

0 0

1

0

Iteration Step

1

Iteration Step

Fig. 1. Algorithm convergence for single emitter.

SNR = 5[dB]

SNR = 20[dB]

140

140 DGML CRB

120

120

Initial

100 RMSE [m]

100

RMSE [m]

DGML CRB

Initial

80 60

100%

40

80 60 40

100%

100% 96% 30% 20

20 100%

0

82% 28%

0 0

1

2

3

4

Iteration Step

0

1

2

3

4

Iteration Step

Fig. 2. Algorithm convergence for two emitters.

(m). In Fig. 3 we show the RMSE of the first emitter using circular arrays in all stations. Each array consists of five elements, four elements equally

distributed on a circle with radius of l=2, and the fifth element (the reference) is located at the circle center. The SNR is 20 (dB). The numbers indicate

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

2354

SNR = 20[dB] 140 DGML CRB

Initial 120

RMSE [m]

100 80 100% 60 40 100% 20

90%

72%

6%

16%

40%

0 0

1

2

3

4

5

6

7

8

Iteration Step Fig. 3. Algorithm convergence for two emitters using circular arrays.

120 DGML CRB 100

RMSE [m]

80

60

40

20

0

-6

-5

-4

-3

-2

-1

0

1

2

3

4

5

6

SNR [dB] Fig. 4. RMSE versus SNR.

N=N exp in percent. As can be seen after only four iterations the CRB is reached. 4.4. Geolocation of two emitters: RMSE versus SNR Assume two emitters ðQ ¼ 2Þ are positioned at p1 ¼ ½201:99; 595:06 (m), and p2 ¼ ½800:14; 502:42

(m). In Fig. 4 we show the RMSE of the first emitter versus the SNR. The SNR varies between 5 and 5 (dB) with 1 (dB) step. The number of sections is K ¼ 50. As expected, the RMSE converges to the CRB as the SNR increases. SNR of 2 (dB) is sufficient for approaching the CRB.

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

2355

80 DGML CRB

70 60

RMSE [m]

50 40 30 20 10 0 2

4

6

8

10

12

14

16

18

20

22

K, Number of sections Fig. 5. RMSE versus number of sections, K.

4.5. Geolocation of two emitters: RMSE versus number of sections ðKÞ Assume two emitters ðQ ¼ 2Þ are positioned at p1 ¼ ½198:41; 600:34 (m), and p2 ¼ ½802:27; 498:09 (m). In Fig. 5 we show the RMSE of the first emitter versus the number of sections, K. The SNR is 10 (dB). As can be seen the RMSE converges to the CRB, even for a small number of sections. 4.6. Geolocation of three emitters Assume that three emitters ðQ ¼ 3Þ are positioned at p1 ¼ ½204:56; 95:88 (m), p2 ¼ ½397:45; 803:04 (m), and p3 ¼ ½702:03; 304:19 (m). In Fig. 6 we plotted the RMSE of the three emitters for SNR ¼ 20 (dB). As can be seen, after five iterations the CRB is reached in 92% of the trials. Only in 2% of the cases 10 iterations are required for convergence. Fig. 7 illustrates the convergence of the algorithm to the true position for a specific single experiment. The labels indicate the iteration number. As can be seen after seven iterations the algorithm converges. The residual estimation error is slightly larger than the error predicted by the CRB. 5. Conclusions An algorithm is proposed for estimating the positions of multiple emitters. The emitters transmit

unknown waveforms that are intercepted by widely separated arrays. The positions are estimated in the fusion center using a direct position determination approach. The algorithm implements the ML criterion in an iterative approach. In each iteration the position of one of the emitters is updated by a planar grid search. This substantially reduces the complexity of a straightforward implementation of the ML estimator. It is shown that convergence to the CRB is achieved with only a small number of iterations. Numerical examples demonstrate the effectiveness of the algorithm.

Appendix A A.1. Derivation of the CRB Given K independent samples of a circular Gaussian process with mean m and covariance matrix K, both depend on a parameter vector l, the ðm; nÞth elements of the Fisher information matrix (FIM), J, are given by [8]   qK 1 qK Jmn ¼ K Tr K1 K qlm qln  H  K X qm ðkÞ 1 qmðkÞ Re K þ2 . ð24Þ qlm qln k¼1

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

2356

First Emitter

RMSE [m]

200 150

Initial

DGML CRB

100%

100

100% 100%

50

100%

92%

42%

66%

14%

0 0

1

2

3

4

5

6

7

8

6% 9

2% 10

11

Iteration Step Second Emitter

RMSE [m]

200 150

DGML CRB

Initial

100 100%

50

100%

100%

100%

0 0

1

2

3

4

92% 5

66% 6

42% 7

14% 8

6% 9

2% 10

11

Iteration Step Third Emitter

RMSE [m]

200 150

DGML CRB

Initial 100%

100

100%

100%

50

100%

92%

42%

66%

14%

0 0

1

2

3

4

5

6

7

8

6% 9

2% 10

11

Iteration Step Fig. 6. The convergence of the algorithm for three emitters.

The CRB equals the inverse of the FIM. The FIM can be partitioned into block matrices, where each block is associated with one or two of the parameter vectors included in l, thus, we can write the FIM as 2

Jss 6 JT J94 sb JTsp

Jsb Jbb JTbp

3

Jsp Jbp 7 5. Jpp

(25)

The block matrices of the FIM are summarized in the following.

qmðkÞ=q~sðkÞ ¼ jAðIJ  BÞ, and obtain that Js¯ðkÞ;¯sðk0 Þ ¼ Js~ðkÞ;~sðk0 Þ 2 ¼ 2 RefðAðIJ  BÞÞH AðIJ  BÞgdk;k0 , s

Js¯ðkÞ;~sðk0 Þ ¼  Js~ðkÞ;¯sðk0 Þ 2 ¼  2 ImfðAðIJ  BÞÞH AðIJ  BÞgdk;k0 .ð27Þ s Define: H9

A.2. The blocks Js;s , Js;b , and Js;p (1) The block Js;s : We use the derivatives: qmðkÞ=q¯sðkÞ ¼ AðIJ  BÞ,

ð26Þ

2 ðAðIJ  BÞÞH AðIJ  BÞ, s2

"

¯ H G9 ~ H

# ~ H . ¯ H

(28)

(29)

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

320 76 300

CRB Radius

2357

True Position Iteration Position

5 4

Y Axis [m]

280 3

260 Initial

240

2

220

1

200 180 660

670

680

690

700

710

720

730

740

X axis [m] Fig. 7. Illustration of the algorithm iterations.

Combine the matrices into a single matrix "

Js¯ðkÞ;¯sðkÞ JsðkÞ;sðkÞ 9 Js~ðkÞ;¯sðkÞ

Js¯ðkÞ;~sðkÞ Js~ðkÞ;~sðkÞ

#

" ¼

¯ H ~ H

~ H ¯ H

Define:

# ¼ G. (30)

The matrix Js;s is a block diagonal matrix with JsðkÞ;sðkÞ along the main diagonal. (2) The block Js;b : We use the derivatives: qmðkÞ=qb¯ q ¼ Wq ðsq ðkÞ  IL Þ, qmðkÞ=qb~ q ¼ jWq ðsq ðkÞ  IL Þ, where sq ðkÞ9½sq ðk; 1Þ; . . . ; sq ðk; JÞT ,

(31)

Wq 9DiagfWq ð1Þ; . . . ; Wq ðJÞg,

(32)

then we obtain that

Dk 9

2 ðAðIJ  BÞÞH WðSðkÞ  IL Þ, s2

W9½W1    WQ ,

(36)

SðkÞ9Diagfs1 ðkÞ; . . . ; sQ ðkÞg.

(37)

Combining all matrices yields: " # " ¯k Js¯ðkÞ;b¯ Js¯ðkÞ;b~ D JsðkÞ;b 9 ¼ ~ Js~ðkÞ;b¯ Js~ðkÞ;b~ Dk

  2 Re ðAðIJ  BÞÞH Wq ðsq ðkÞ  IL Þ , 2 s

Dx;q 9 ð33Þ

qWq ; qxq

(38)

Dy;q 9

qWq . qyq

(39)

Then,   2 Re ðAðIJ  BÞÞH Dx;q ðIJ  bq Þsq ðkÞ , s2   2 ¼ 2 Re ðAðIJ  BÞÞH Dy;q ðIJ  bq Þsq ðkÞ , s   2 ¼ 2 Im ðAðIJ  BÞÞH Dx;q ðIJ  bq Þsq ðkÞ , s

Js¯ðkÞ;xq ¼ Js¯ðkÞ;b~ q ¼  Js~ðkÞ;b¯ q   2 ¼  2 Im ðAðIJ  BÞÞH Wq ðsq ðkÞ  IL Þ . s ð34Þ

# ~k D . D¯ k

(3) The block Js;p : We now use the following derivatives: qmðkÞ=qxq ¼ Dx;q ðIJ  bq Þsq ðkÞ, qmðkÞ= qyq ¼ Dy;q ðIJ  bq Þsq ðkÞ, where

Js¯ðkÞ;b¯ q ¼ Js~ðkÞ;b~ q ¼

(35)

Js¯ðkÞ;yq Js~ðkÞ;xq

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

2358

Js~ðkÞ;yq ¼

  2 Im ðAðIJ  BÞÞH Dy;q ðIJ  bq Þsq ðkÞ . 2 s ð40Þ

where K 2X ðSðkÞ  IL ÞH WH WðSðkÞ  IL Þ. s2 k¼1

R9

(51)

Define: Dq 9½Dx;q ; Dy;q ,

(41)

2 Xk;q 9 2 ðAðIJ  BÞÞH Dq ðI2  ðIJ  bq ÞÞðsq ðkÞ  I2 Þ. s (42) Combining the terms we get " # " # ¯ k;q Js¯ðkÞ;xq Js¯ðkÞ;yq X JsðkÞ;pq 9 ¼ ~ . Js~ðkÞ;xq Js~ðkÞ;yq Xk;q Next define: h D9 D1   

(2) The block Jb;p : In this case we get that Jb¯ q ;xq0 ¼

Jb¯ q ;yq0 ¼

(43)

K 2X H ImfðsH q ðkÞ  IL ÞWq Dx;q0 ðIJ  bq0 Þsq0 ðkÞg, 2 s k¼1

Jb~ q ;yq0 ¼ 

K 2X H ImfðsH q ðkÞ  IL ÞWq Dy;q0 ðIJ  bq0 Þsq0 ðkÞg. s2 k¼1

ð52Þ

i DQ ,

Define:

2 Xk 9 2 ðAðIJ  BÞÞH DWðSðkÞ  I2 Þ, s

(45)

P9

W9DiagfI2  ðIJ  b1 Þ; . . . ; I2  ðIJ  bQ Þg.

(46)

Combining all matrices in a single matrix: " # h i ¯k X JsðkÞ;p 9 JsðkÞ;p1    JsðkÞ;pQ ¼ ~ . Xk

K X

2 H RefðsH q ðkÞ  IL ÞWq Wq0 ðsq0 ðkÞ  IL Þg, s2 k¼1 ð48Þ

Jb~ q ;b¯ q0 ¼  Jb¯ q ;b~ q0 K 2X H ImfðsH q ðkÞ  IL ÞWq Wq0 ðsq0 ðkÞ  IL Þg. s2 k¼1 ð49Þ

Collecting all matrices: Jb;¯ b~ Jb;~ b~

#

" ¼

¯ R ~ R

# R~ , ¯ R

(54)

In this case

Jb¯ q ;b¯ q0 ¼ Jb~ q ;b~ q0

Jb;¯ b¯ Jb;b 9 Jb;~ b¯

Arranging in a matrix form: " # " # Jb;p ¯ ¯ P Jb;p 9 ¼ . Jb;p ~ ~ P

(53)

A.4. The block Jp;p

(1) The block Jb;b : In this case we obtain that

"

K 2X ðSðkÞ  IL ÞH WH DWðSðkÞ  I2 Þ. s2 k¼1

(47)

A.3. The blocks Jb;b and Jb;p

¼

K 2X H RefðsH q ðkÞ  IL ÞWq Dy;q0 ðIJ  bq0 Þsq0 ðkÞg, s2 k¼1

Jb~ q ;xq0 ¼ 

(44)

¼

K 2X H RefðsH q ðkÞ  IL ÞWq Dx;q0 ðIJ  bq0 Þsq0 ðkÞg, 2 s k¼1

(50)

Jxq ;xq0 ¼

K 2X H H RefsH q ðkÞðIJ  bq Þ Dx;q Dx;q0 ðIJ  bq0 Þsq0 ðkÞg, 2 s k¼1

Jxq ;yq0 ¼

K 2X H H RefsH q ðkÞðIJ  bq Þ Dx;q Dy;q0 ðIJ  bq0 Þsq0 ðkÞg, 2 s k¼1

Jyq ;xq0 ¼

K 2X H H RefsH q ðkÞðIJ  bq Þ Dy;q Dx;q0 ðIJ  bq0 Þsq0 ðkÞg, 2 s k¼1

Jyq ;yq0 ¼

K 2X H H RefsH q ðkÞðIJ  bq Þ Dy;q Dy;q0 ðIJ  bq0 Þsq0 ðkÞg. s2 k¼1 ð55Þ

Arranging the terms in a matrix form yields: 2 3 Jxq ;xq0 Jxq ;yq0 5 Jpq ;pq0 94 Jyq ;xq0 Jyq ;yq0 ¼

2 H H H RefðsH q ðkÞÞI2J  I2 bq ÞDq Dq0 ðI2 bq  I2J sq ðkÞÞg. s2 ð56Þ

ARTICLE IN PRESS A. Amar, A.J. Weiss / Signal Processing 87 (2007) 2348–2359

Define: X9

K   2X Re ðSðkÞ  I2 ÞH WH DH DWðSðkÞ  I2 Þ . s2 k¼1

(57) 2

Arranging all Q matrices in a matrix single form: 2 3 Jp1 ;p1    Jp1 ;pQ 6 . .. 7 .. 6 7 ¼ X. (58) Jp;p 96 .. . 7 . 4 5 JpQ ;p1

JpQ ;pQ



Arranging all results in a single matrix yields: 2

¯ H 6 ~ 6 H 6 6 T 6 0 6 6 . 6 .. 6 J96 T 6 0 6 6 ¯T 6 D1 6 6 ~T 6 D1 4 ¯T X 1

~ H ¯ H

0 ..



0

.

.. .

D~ 1 D¯ 1 .. .

¯1 D ~1 D

~ H ¯ H ~T D

~K D ¯DK ~ R ¯ R

~T D 1



¯ H ~ H ¯T D

K

D¯ K D~ K ¯ R

¯T D 1 ~T X 1



~ D K

¯T D K

~ R



¯T X K

~T X K

K T

¯ P

T

~T P

3 ¯1 X ~1 7 X 7 7 .. 7 . 7 7 7 ¯K7 X 7 . ~K7 7 X 7 ¯ 7 P 7 7 ~ 7 7 P 5 X

(59) This concludes the derivation. References [1] D.J. Torrieri, Statistical theory of passive location systems, IEEE Trans. Aerospace Electronic Systems AES-20 (2) (March 1984) 183–198.

2359

[2] W. Foy, Position–location solutions by Taylor-series estimation, IEEE Trans. Aerospace Electronic Systems AES-12 (March 1976) 183–193. [3] M. Wax, T. Kailath, Decentralized processing in sensor arrays, IEEE Trans. Acoust. Speech Signal Process. ASSP-33 (4) (October 1985) 1123–1129. [4] J.C. Chen, R.E. Hudson, K. Yao, Maximum likelihood source localization and unknown sensor location estimation for wideband signals in the near-field, IEEE Trans. Signal Process. 50 (8) (August 2002) 1843–1854. [5] R. Kozick, B.M. Sadler, Source localization with distributed sensor arrays and partial spatial coherence, IEEE Trans. Signal Process. 52 (3) (March 2004) 601–616. [6] A. Amar, A.J. Weiss, Direct position determination of multiple radio signals, EURASIP J. Appl. Signal Process. 2005 (1) (2005) 37–49 (Special issue: Advances in Sensor Array Signal Processing Technology). [7] A.J. Weiss, A.S. Willsky, B.C. Levy, Maximum likelihood array processing for the estimation of superimposed signals, Proc. IEEE 76 (2) (February 1988) 203–205. [8] H.L. Van Trees, Optimum Array Processing, Detection Estimation, and Modulation Theory—Part IV, Wiley, New York, 2002.