arrival (DOA) of the signal sources in these subspaces. In establishing the signal and noise model, an assump- tion common to these high resolution methods is ...
Estimation of the directions of arrival of spatially dispersed signals in array processing Y. Meng P. Stoica K.M.Wong
Indexing terms: Array processing, Directions of urrival
Abstract: The problem of estimating the DOA of spatially distributed signals is examined. A mathematical model is first established by making some reasonable assumptions. The correlation matrix of a distributed signal is then derived. The important properties of the correlation matrix are studied, revealing that even if the matrix is of full rank (being equal to the number of sensors), which renders conventional high resolution array processing methods inapplicable, the dimensionality of the signal subspace can be approximated to a number usually much smaller than the number of sensors. From the observation, the quasi-signal and noise subspaces are identified, and, utilising the orthogonality of the signal and noise subspaces, an algorithm (DISPARE) has been developed. Analytic studies and numerical evaluations are carried out to examine the performance of the algorithm under various environments and show that it is indeed effective.
1
Introduction
In recent years, high resolution array processing techniques have been widely applied in areas in which very different wave phenomena occur. These include seismic exploration, radar, passive sonar, radio astronomy, as well as radio communications. A prime objective of array processing is to estimate the location of signal sources. To this end, these high resolution methods [l151 generally start with a signal and noise model, and based on the covariance matrix of the received data, the data space is partitioned into the signal subspace and noise subspace. Various methods are then applied to estimate the parameters concerning the directions of arrival (DOA) of the signal sources in these subspaces. In establishing the signal and noise model, an assumption common to these high resolution methods is that the signals are propagated from distinct point sources. This assumption is, at best, an approximation of the 0IEE, 1996 IEE Proceedings online no. 19960170 Paper first received 20th March and in revised form 20th November 1995 Y. Meng and K.M. Wong are with the Communications Research Laboratory, McMaster University, Hamlton, Ontario, Canada L8S 4Kl P. Stoica is with the Systems and Control Group, Department of Technology, Uppsala University, Uppsala, Sweden IEE Proc.-Radar, Sonur Nuvig., Vol. 143, No. 1, February I996
practical environment. The ‘size’ of a source depends on its distance from the array of sensors. In addition, signal sources, in practice, may often be transmitted by reflection and the reflective medium may often be dispersive and thus violating the ‘point source’ assumption. In the case of communication in the Arctic environment, for instance, the transmission of radio waves 11s often performed by ionospheric scattering. The signal having undergone reflection through a scattering volume reaches the receiver and would appear to be from a distributed source. The receiver, in this application, has to determine the DOA of this distributed signal. ,4second example arises in the case of a radar system performing low angle tracking. In this case the aircraft, if it is sufficiently far away, can be approximately modelled as a point source. However, the diffuse cornponents caused by the reflection from a rough sea surface would appear as a distributed source. Here, the objective may be to estimate the DOA of a point source in the presence of distributed sources. With the exception of very few [16-221, most of the existing methods of high resolution direction estimation rely on the assumption that the signals are from point sources. Thus, the application of the conventional high resolution methods to such problems mentioned above will show grave deterioration in performance. In this paper, we examine the problem of array processing when the signal sources are distributed in space (diffuse) and propose a simple method of estimating their DOA. The performance of the method is analysed and computer simulation results are presented. 2
Signal model
For simplicity, we assume that the distributed signals are on the same plane as the array, although extension to a three-dimensional space is straightforward. We model an entire diffuse narrowband source to be made up of individual point sources closely clustered together. For the kth diffuse source, its lth point source component emits energy at an angle qkl with respect to the array of sensors, and the intensity (standard deviation) of the source is slcl,1 = 1, ..., L,; k = 1 ..., K. Thus, the output of the array at the nth snapshot for K statistically in dependent diffuse sources is given by K
LL
where x(n) is an M-dimensional vector containing the output signals of the A4 sensors in the array. d(&J is 1
the direction manifold, being an M-dimensional complex vector characterised by the parameter Gkkl associated with the Ith component of the kth source, and v(n) is the A4 x 1 complex vector assumed to be ergodic, independent of the signals, with zero mean and covariance q?I, where qv is an unknown constant and I is the identity matrix. The Lk point source components within the kth diffuse source are further assumed to be independent of each other. This assumption is often (and justifiably) used to simplify calculations [23, 241 in many applications and is also supported by recent experimental results in the cases of tropospheric and ionospheric scatterings [25, 261. Under such assumptions, the covariance matrix of the array output is given by K
k=l
where t denotes the Hermitian conjugate of a vector or matrix, and RI, is the covariance matrix of' the output due to the kth diffuse signal, and is given by L&
(3)
2. rk is an M x M Hermitian matrix generally of full rank and is dependent only on the shape of the distribution of the signal but independent of &. 3. R, and rk have the same eigenvalues and their respective eigenvectors are related by m = 1,I . . , M (9) V k m = LkUkm with Ykn, and Ukln being the mth eigenvector of Rk and r,, respectively. The first two properties can be observed easily by using the analogy of the Fourier transform in eqn. 5 , and the third property follows directly from applying the definition of eigendecomposition to eqn. 7. 2.1 Example 7 Suppose that we have a linear array having M sensors uniformly spaced at half the wavelength of the narrowband signal and that the kth signal power density is uniformly distributed within [ F k - o k , + Ok], i.e. I O elsewhere Then, using eqn. 5 , for which d(4) = [1 the mnth element of RI, is given by
... eJ@-*)@IT,
l=l
with (4) which is the average power of the Ith point-source component of the kth difhse signal. If these point-source components of the same diffuse signal are closely clustered together, we replace each of the point sources by a density function spreading over a small angle A$ so that, as A@ -+ 0, eqn. 3 becomes r L~ 1
In this case the mnth element of rkis simply sin(m - n)ok m, n = 1, . . . , A4 (12) ?kmn = ( m- n ) a The mth eigenvector Ukm U,(@ of rk is formed by the elements of the well known discrete prolate spheroidal sequences (DPSS) [27],and the corresponding mth eigenvalue is given by M
2
= 41,
LE@&
d(4)dt(4)Pk(4)d4
(5)
Here, pk($) is a normalised positive density function characterising the distribution of the kth signal such that pk(@)d@
=
(6)
with Qk being the region within which the diffuse signal spreads, and q z is the average power of the kth diffuse signal. A particular distribution function P k ($) gives rise to a particular covariance matrix of the signals at the sensor outputs. In general, a distribution function has two important parameters & and o k , where Fk denotes the mean DOA of the kth signal with respect to the array and ok is a measure of the 'spread' of the signal. Let us examine the covariance matrix Rk in the specific case of a uniform linear array which is one of the most commonly used arrays. In the case of a linear array having A4 sensors uniformly spaced, it can be observed that the covariance matrix RI, has the following properties: 1.
R~= L k r k Q
(7)
where Lk is a diagonal matrix involving only the mean location $, of the signal such that
Akm
=
m = 1,. . . ,Ad (13)
2=1 m
z=-m
where u,,(qJ denotes the ith element in the eigenvector uk,,,. We note that hk, is a function of 0 , and is the ratio of the energy concentrated in the M elements of the DPSS to the energy contained in the entire sequence. Table 1 shows the numerical values of h k m for various values of ok when M = 8. Table 1: Eigenvalues of a uniformly distributed signal ok
0.2682
0.3491
0.5236
hkl h,,
6.5799 1.3654
5.3928 2.3460
hk3
7.1366 8.4596 x IO-' 1.7350 x IO-,
5.4000 x IO-*
2.5322 x IO-'
hk4
1,2209 x IO"
6.8853 x IO4
7.8551 x
hk5 h,, h,,
4.0704 x
4.1023 x
1.0735 x IO4
6.9532 x 10-Io
1.2518 x
7.4453 x
1.1421 x 10-I2
1.6884~
2.5461 x
hks
5.5232 x IO-I3
2.0202 x
1.1466 x
2.2 Example2 For a uniform linear array, if the kth signal power density is Gaussian distributed such that
IEE Proc -Radar, Sonar Navig., Vol. 143, No. 1, February 1996
then, again using eqn. 5 , the mnth element of Rk is given by Tkmn
e-
I
(n--n)22
z
,j(--n)Sb
0.81
(15) Table 2 shows the numerical values of the eigenvalues A,,,., of the matrix Rk for various ‘spread’ of 0, when M = 8. Table 2: Eigenvalues of a Gaussian distributed signal ~~
om
~
0.0583
0.1163
0.1745
7.6969
7.4822
6.9524
2.9846 x IO-’
5.0369 x IO-’
9.8349 x IO-’
4.5987 x I 0-3 4.1575 x 10”
1.3898 x 10-2
6.1805 x IO-*
2.2466 x IO4
2.2810 x IO3
2.4320 x
2.3514 x
5.4649 x
9.2058 x 10-lo
1.5953 x 10-8
8.5213 x
1.7032 x 10-l2
6.4096 x 10-I’
7.9857 x
3.7873 x IO-l3
5.5602 x
3.3071 x IO-’’
Since the covariance matrix Rk of the kth distributed signal is of full-rank, the direct application of conventional high resolution array processing methods such as MUSIC or MLE will produce erroneous results [19, 281. To alleviate the difficulty, we examine the eigenstructure of R,. For a practically reasonable spread 0, of the signal, it is noted that, in general, the A4 eigenvalues of Rk decrease very quickly in magnitude, i.e. the energy of the signal mainly concentrates in the first few eigenvalues. In Example 1 (Table 1) for instance, we note that for a total spread of 0,= 0.5236rad (‘electric angle’) on each side of the mean angle of arrival (corresponding to 9.59” of physical spread on each side for a signal arriving in the direction normal to the array), over 96% of the signal energy is concentrated in the first two eigenvalues. Similarly, for a spread of 0, = 0.2682rad on each side of the mean DOA, the first two eigenvalues total over 99% of the signal energy. This observation is equally evident in the case of a Gaussian distributed signal. In general, the number of eigenvalues in which a dominant proportion of signal energy concentrates depends on the number of sensors M , the extent to which the signal spreads (0,)and, to a lesser degree, the distribution of the signal. The larger is M , and/or the larger is the signal spread, the more eigenvalues have to be accounted for in the concentration of signal energy. Figs. 1 4 show the plot of the eigenvalues against the number of sensors ( M = 1, ..., 20) for the cases of uniformly and Gaussian distributed signals with different spreads. The eigenvalues in each case are normalised to the largest eigenvalue for M = 20. The eigenvalues for other signal distributions have also been evaluated [21]. Similar observations that, for a limited spread, the signal energy concentration is dominated by the first few eigenvalues persist. It is further observed [21] that, with the increase of M and/or o k , the energy concentration spreads more quickly to other eigenvalues in the case of a uniformly distributed signal than signals with other distributions. From the above observations we can infer that, even though the matrix Rk is of full rank, because of the energy concentration in only mok (< A@ of the dominant eigenvalues, the dimension of the signal subspace induced by a distributed signal can be approximated to m o k , This quasi signal subspace is spanned by the eigenvectors corresponding to the largest mOkeigenvalues of Rk. Since the IEE Proc -Radar, Sonar Navig , Vol 143,No I, February 1996
0
2
Fig. 1 Eigenvalues
4 of
6
8
10 M
12
14
16
18
20
uniformly distributed signal with 0, = 0.174Srad
m
-3d 5
.-01
W
6
8
10 M
12
14
16
18
20
Fig.2 Eigenvalues of uniformly distributed signal with 0, = 0.3491rad
- ’- .OI 0.8 m
-2 0.6-
9
-
.-b 0.4-
W
K2
o 0
,
,
.
2
4
6
8
10 M
12
4 14
16
18
20
Fig.3 Eigenvalues of Gaussian distributed signal with 0, = O.OS83rad
M
Fig.4 E&envalues of Gaussian distributed signal with 0, = 0.1164rad 3
uniform distribution represents the widest spread of signal energy among the dominant eigenvalues, i.e. its signal subspace has the largest approximate dimension, we can treat this as 'worst case', so that for a given number of sensors and an estimated spread, we can examine from graphs such as Figs. 1 and 2 the number of dominant eigenvalues and determine the approximate signal subspace dimension accordingly. We further note that, if there are K distributed signals, then the number of sensors must satisfy M > Cf=, mQk.
3 DlSPARE - DOA estimation for distributed signals From the signal model described in the preceding Section we see that there are two essential parameters & and o k which characterise the mean location and the spread of the kth distributed signal, and have to be estimated using the data covariance matrix R, given by eqn. 2. We assume that the number of distributed signals K is known, and that the signals are uncorrelated with each other. We further assume that each of these signals is of a particular distribution and has mOkdominant eigenvalues. Thus, we can express the covariance matrix due to the K distributed signal as
where In practice, we do not know R, exactly; we can only obtain an estimate of it such that
c N
1
R, = /v L
.
x(n)x+(n)
n=l
It is well known [29, 301 that lim,,+- R, R, at a convergence rate of where indicates that the equality holds with probability 1. From R, we can obtain the eigenvectors {w,} corresponding to the eigenvalues {bnl},m = 1, ..., M . Since R, converges asymptotically to R,, and from eqn. 19, we can write lim
A+m
vskWL
= a.s.
o
k = 1,* , . , K
(22)
where
w, = [WJ,f0+l . ' . Wil'M] (23) These asymptotic orthogonality properties are the bases for the development of the algorithm for the estimation of the mean DOA of the distributed signals in this Section. The method is called dispersed signal parametric estimation (DISPARE). Let us exploit the orthogonal relationship of eqn. 22. Due to the asymptotic orthogonality between Vskand M(,, we can write
k=l
Now, we define the correlation matrix formed by the eigenvectors associated with only the dominant eigenvalues such that
but from the approximation of the rank of the signal covariance matrices as indicated by eqn. 17, eqn. 24 can be rewritten as Wo ~ lim R ~ =
N+cc )
k=l
Let MO= CkK,, mOic< M ; since hknzare negligible for MO m M , we can say R, = S. Thus, we can say that the rank of R, is approximately MO.Thus, if we let { w l , ..., wlwo} be the eigenvectors corresponding to the largest MO eigenvalues of R,y,then each w,, rn = I , .., MO,is a linear combination of {vi,,; m = 1, ..., mok;k = 1, ..., K}.In developing the above argument, we have tacitly assumed that none of the Xk, that we have neglected in the approximation from eqn. 16 to 17 is larger than any of the hkn, we have retained. This set of eigenvectors { w l , ..., wMo} spans the quasi-signal subspace S. Now, the covariance matrix R, of the received data is the sum of R, and the covariance matrix of the spatially white noise as indicated by eqn. 2. Therefore, corresponding to the M eigenvalues p,, ..., p M , the eigenvectors of R, are given by {w,, ..., wM},where the first M Oeigenvectors {w,}, m = 1, ..., MO,are identical to those of R,, whereas the last ( M - MO)eigenvectors corresponding to the least significant eigenvalues of R, span the quasi-noise subspace N . If we form the quasi-signal eigenvector matrix such that W,>= [w, ... wMo] and the quasi-noise eigenvector matrix such that W, = [wM0+, ... wM],we have the orthogonal condition that
wiw, = 0
(18) Since the columns of W, are linear combinations of {vkm),we also have the orthogonal condition that
4
(25)
where Rk is the covariance matrix due to the kth distributed signal. Thus, we can conceive an algorithm in which we adjust the parameters Fk and q,in Rk until eqn. 25 is satisfied. Since the square of the Frobenius norm of a matrix is the measure of the sum of the square of each element of the matrix, we can arrive at the estimates such that
{&, G k } = arg min
llRkv27.vll2F
(26)
4Jk> U k
where ll,l\F denotes the Frobenius norm. Eqn. 26 yields the DISPARE estimates of the parameters 6'and o k of the kth signal. If we form the inner products of eigenvectors of Rk and the estimated quasi-noise eigenvectors so that ekirn = vLmWiTZ i = MO + 1,.. . , M ; m = 1,.. . , M (27) then" it can easily be shown by using the relation IIRkW((2= t r { q R k R $ W,) that the DISPARE criterion amounts to the minimisation of E,, hi:, 1 ckim 1 2 , i.e. DISPARE employs a weighting on the inner products of the eigenvectors of Rk and the estimated quasinoise eigenvectors.
4
Application of the DISPARE algorithm
Let us first test the effectiveness of the DISPARE algorithm by carrying out some simulation examples. In each of these examples, an 8-sensor linear array with uniform spacing of a half-wavelength of the signal source is used. All sensors are omnidirectional and have unity gain. Isotropic Gaussian noise with zero mean is added to the received signal. The number of snapshots is taken to be 200, and in each case an IEE Proc.-Radar, Sonur Novip., Vol. 143, No. I , February 1996
exhaustive search for the minimum of the criterion i s located in the & o k plane to obtain the estimates qk and &/'. The root mean square (RMS) of the estimation error (average over 500 trials) are plotted for different signalhoise ratios (SNR).
-
4.1 Example 3 In this example we have the situation in which there is only one uniformly distributed signal arriving at a mean DOA of &( = - 0.5454rad, corresponding to a physical angle of -10" from the normal of the array. We study two cases in which the signal has spreads o k = 0.1745rad and o k = 0.3491rad, corresponding to total physical width of 6.4" and 13.4", respectively. Examining the eigenvalues of the covariance matrix Rk from plots in Figs. 1 and 2, using the threshold that over 95% of the signal energy be included in the quasisignal subspace, we determine that the approximate dimensions of the signal subspace in each case are 1 and 2, respectively. We now apply ;he criterion given by eqn. 26 to obtain the estimates and c f k . Since we are, in general, only interested in the DOA, we only plot the RMS values of the estimation error for the two cases. These are shown in Fig. 5 for the two different signal spreads at various SNRs. As a comparison, we also employ MUSIC [8] to estimate the DOA (given the knowledge that K = 1). It can be observed that for small signal spread, MUSIC and DISPARE perform almost identically (the performances of the two methods coincide in Fig. 5). However, for a large signal spread MUSIC completely breaks down, yielding unacceptably large errors over the whole range of SNR.
$10-1
-
n n C
z 10-2
w
m I (L
m-3
0
5
10
15
20 25 SNR, dB
30
35
LO
Fig.5 Estimation errors of mean DOA of a unlfoimly distributed rignal
Case 1 solid line (ok= 0 1745rad), case 2 dashed line (ok= 0 3491 rad) x DISPARE 0 MUSIC
4.2 Example 4 In this example there is again only one signal arriving at a mean DOA of Fk = -0.5454rad. However, here the signal is Gaussian distributed. The spread o,is chosen to be 0.0583rad and 0.1163rad, corresponding to a physical spread of 1.1 and 2.2", respectively. We note that over 99% of the energy of the signal energy is collected in the range within (& & 30k), making the energy spread of the signals about the same as the corresponding ones in Example 3 . From the plots of the eigenvalues (Figs. 3 and 4), again using the criterion that over 95% be included in the quasisignal subspace, the approximate dimensions of the signal subspace are 1 and 2, respectively. The criterion of the DISPARE algorithmhis applied and the RMS of the estimation error of $k for the two cases is shown in Fig. 6, which O
IEE Proc -Radar, Sonar Navig
,
Vol 143 N o I , Febiuary 1996
also shows the RM[S of the estimation error using MUSIC. Results similar to those in Example 3 are observed here. For small signal spread, MUSIC and DISPARE perform almost identically. For a large signal sprlead, however, the MUSIC method again yields unacceptably large errors. 10 l3
e
0
,o-'l
w
20 25 30 35 40 SNR, dB Fig. 6 Estimation errors ofa mean DOA ofu Gaussian distributc~dsignul Case 1 solid line (ok= 0 0583) case 2 dashed line ( 0 1 ~ = 0 1 Ihlrad) 1op40
5
10
15
x DTSPARE 0
MUSIC
The results in Examples 3 and 4 are not surprising. When the spread of a signal is small, it behaves like a point source (with approximate signal dimension equal to one) Thus, the methods developed here will have similar performance as MUSIC, which is based on the assumplion of a point source. When the spread of a signal iis large, i.e. nzOk 1, there will be erroneous peaks in the MUSIC spectrum, leading to large error in the estiination of the mean DOA, which will not be effective at all. In the application of the DISPARE algorithm and in the two examples above, we have assumed that the distribution of the kth dispersed signal is known. This enables us to establish the functional form of the signal covariance matrix Rk so that the optimisation criterion of eqn. 26 can be correspondingly set up. Furthermore, we also assumed that the spread of the signal is at least approximately known so that the approximate dimension of 1.he signal subspace can be determined. In practice, neither the distribution nor the spread of the signal is usually known for us. In the following example, we test how the DISPARE algorithm performs when these assumptions are not valid.
4.3
Example 5
In this example we study the effect of making a wrong assumption on the distribution of the spread signal. The scenario is identical to that in Example 3, in which a uniformly distributed signal of spread o,arrives at a mean DOA of $k = -0.5454rad to the normal of a uniform linear array. However, we assume that the signal is Gaussian distributed and accordingly we apply the signal covariance matrix RI, as given by eqn. 15 and employ the DISPARE criterion of eqn. 26 to estimate the parameters i&, 9.We assume that we know roughly the signal spread 0,. We tested two cases in which thLe true ok are 0.1745rad and 0.3491 rad, respectively. Now, to encompass 99% of the signal energy within this range of spread under the assumption of a Gaussian distribution, we assume that the Gaussian distribution in the two cases to have 0, equal to 5
0.0583rad and 0.1 163rad, respectively, and therefore, as in Example 4 above, the respective approximate dimensions of the signal subspaces are 1 and 2. The RMS errors of the estimates (&k are plotted in Fig. 7. Compared to the results shown in Fig. 5 , we observe that the performance of the DISPARE algorithm has deteriorated. The wider is the signal spread, the larger are the estimation errors. 10 O
g 10-1 n C 0
E
::lo-* L
v)
. . . . . . . . . . . . . . . . . . . . . . . . .
20 25 30 35 40 S N R , dB Fi .7 Estimation errors of mean DOA of uniformb distributed signal un%er the assumption of Gaussian distribution
0
10
5
15
Case 1: solid line (ok= 0.1745rad); case 2: dashed line (ok= 0.3491rad)
4.4 Example 6 We continue to examine the effect of making a wrong assumption on the distribution of the spread signal in this example. This time the signal is Gaussian distributed. The scenario is identical to that in Example 4, in which the signal with spread (sk arrives at a mean DOA of -0.5454rad to the normal of the array. However, we assume that the signal is uniformly distributed and limited to a spread of &i k 30k. The two cases we study are 0, = 0.0583rad and o k = 0.1163rad. As in Example 3 above, the corresponding approximate dimensions of the signal subspace are decided to be 1 and 2. The RMS of the errors in these two cases is plotted in Fig. 8. Compared to the results in Fig. 6, it is observed that the performance of the DISPARE algorithm shows hardly any deterioration. This is true regardless of the spread of the distribution. Many distributions other than Gaussian have been tested assuming that the signal is uniformly distributed. Similar observations of negligible deterioration in performance persist [21].
The results in Examples 5 and 6 are of great significance. They show us that, no matter what shape the signal distribution may be, if it is assumed to be uniformly distributed, the estimation error will be about the same as in the case where the distribution is known. Other assumptions may lead to poor deterioration of performance if the assumption is incorrect. We now examine the problem of estimating the approximate dimension of the signal subspace. For a given number of sensors, the number of significant eigenvalues of the correlation matrix due to a dispersed signal depends on the spread of the signal. In many applications such as ionospheric scattering we can often roughly estimate how much the spread o k may be [26]. Furthermore, the number of significant eigenvalues chosen usually is not too sensitive to the estimation error of the spread of the signal. This can be illustrated by Fig. 9, in which we show the relative magnitudes of the eigenvalues of Rk for various signal spread values. The signal is uniformly distributed and the number of sensors in the uniform linear array is eight. As can be observed in Fig. 9, for a practical wide range of signal spread (from 0.2rad to 0.55rad on each side of the mean angle of arrival, corresponding to a physical spread of 3.65" to 10.08' on each side), over 95% of the signal energy is contributed by only the first two eigenvalues. Thus, if the spread of signal lies within this wide range, we can choose the approximate dimension of the signal subspace to be 2, and even if there exists some error in the rough estimate of the signal spread, the number of dominant eigenvalues in R, will not be affected.
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 . 9 1 . 0 signal spread on each side, rad fig.9 Eigenvalues of uniformly distributed signal given M = 8 0
0.1
10-1
U
2
,
.....
.
....................................
n
5E l o - * W
m
2
n
. . . . . . . . . . . . 10-
5
20 25 30 35 40 SNR, dB fi .8 Estrmation errors of mean DOA of Guusian dzstributed signal un%r the assumption of unlform distribution Case 1 solid line (ok= 0 0583rad), case 2 dashed line ( G =~ 0 1163rad) 6
0
10
15
Equipped with the above observations, the procedure of applying the DISPARE algorithm to estimate the mean DOA of dispersed signals is clear and is summarised as follows: 1. We assume that the K dispersed signals are all uniformly distributed and are independent of each other. 2. The spread of each signal is roughly estimated so that, for M sensors in the array and a particular array geometry, the number of dominant eigenvalues m0k can be obtained from charts similar to Fig. 9. 3. The dimension of the quasi-noise subspace can now be determined. An eigendecomposition of the estimated covariance matrix of the data is performed, and the orthogonal projector onto the estimated quasi-noise subspace can be formulated using the noise eigenvector matrix W,. IEE Proc-Radar, Sonar Navig., Vol. 143, No. 1, February 1996
4. The DISPARE algorithm can now be applied to estimate the mean DOA of the dispersed signals. Since all the signals are assumed to be uniformly distributed, the functional forms of the covariance matrix R, due to each signal will be identical, and are completely predetermined using eqn. 5 . The criterion for DISPARE can be rewritten as
{&,
~
k
a
and $2 = -0.9708rad for the Gaussian one. We note that there is significant physical overlap between the two signals. The estimation accuracies at various SNRs for the two signals are shown in Fig. 11. We observed that, wihile the results are still acceptable, there is considerable deterioration in the accuracy. 5
Performance analysis of DISPARE
=} arg Inin j ( 4 ,a ) = arg s i n I I R ~ ( ~ , ~ ) W ~ ~ / $ d 'kh',
$ k r'k
= arg p i n
tr{WtRU;(4, a)RL(@, a ) W u } (28)
d'k>uk
where RJq, 0)is the correlation matrix due to a uniformly distributed signal. The mean DOA $ can then be estimated by locating the K minima of eqn. 28. The following example illustrates this procedure.
4.5 Example 7 In this example we apply the procedure as outlined above to estimate the mean DOA of two independent narrowband dispersed signals, one uniformly distributed (0,= 0.3491rad), the other Gaussian distributed (02= 0.0583rad). In the first case, the two signals are well separated with TI = -0.5454rad for the uniform one and g2 = -2.4063 rad for the Gaussian one. The estimation accuracies at various SNRs for the two signals are shown in Fig. 10, which when compared to the corresponding results in Figs. 5 and 8 show little deterioration. In the second case, the signals are much closer to each other with TI = -0.5454rad for the uniform one IO"
In this Section we present an analysis on the performance of the DISPARE algorithm. The statistical properties of the perturbation in the estimated correlation matrix R, are derived and the effects on the DOA estimates are studied. We first present some of these statistical properties. Consider the case when we have K dispersed signals. From eqns. 16 and 17 we can formulate the orthogonal projectors onto the quasi-signal and the quasi-noise subspaces, respectively, such that P,
=vvT,w: = [WI ' . '
where {w,}, rn = 1, ..., M , are the eigenvectors of R,. Since we only have the estimate R, given by eqn. 21, we define the following errors:
AR, = Rx - R, (31) We now present two theorems governing the distribution of pertinent random quantities due to the perturbation ARx.
5.1 Theorem I Given tlhat the observation consists of N independent complex, zero mean, circular Gaussian random vectors, the elements of AR, are asymptotically Gaussian: E{vec(ilR,)} = 0 5
0
15
10
20
25
30
35
40
SNR, dB
Fi . l o Estimation errors of mean DOA for unformly and Gauuian distrifuted signals
(q,,0 , ) = (-0
8484, 0 3491) and
(32)
1
n
E{vec(AR,)vect(AR,)} = -RZ 8 R, = E, N
(33)
(g2,02)= (-2 4063, 0 0883) (34)
10 O
i
1
.. .
......,j...... ..
i
.. . .... ............
I
where I ' is a permutation matrix defined as M
M
(35)
ed
U
l=1 m=l
with El,, being an M x M matrix having unity as its Zmth element and zero elsewhere, vec(.) denoting the vec operation on a matrix, and 0 denoting the Kronecker product of two matrices [30]. The proof of this theorem can be found in [30].
Q
0
5
10
Fig. 11 Estimation errors
15
20 25 S N R dB
30
35
40
of mean DOA for overlapped uniformly and
Gaussian distributed signals (c,, q)= (-0.5454, 0.3491) and
(f2.
IS*)
= (-0.9708, 0.0583)
IEE Proc-Radar, Sonar Navig., Vol. 143, No. 1, February 1996
5.2 Theorem 2 If Y L spanR,, i.e. if the columns of Y are in the subspace spanned by the columns of R,T,then vec(P,Y) is asymptotically Gaussian distributed with mean value and covariance given by E{vec(P,Y)} = vec(P,Y)
N
0
(36) 7
cov{vec(P,Y),vec(P,Y)}
= ( Y T @ I)C,C,Ct,(Y” c 3 I)
cov{vec(i.,Y), vecx (P,Y)
1 = (Y~@I)C,TC,C: (YBI)
(37)
in the Taylor expansion of eqn. 43 become significant. This results in the small discrepancy of the two curves in low SNR. lo-’
(38)
where G , is defined as n
c, = S#T
63 P,
+ PT 8 s#
U
(39)
2
and
Q
0
a
n
E{[u - E(u)][v- E(v)]*} The proof of this theorem neccessitates the properties of perturbed projectors and can be found in [21]. We now introduce the notations that if A is a function of $and 0; then COV(U, V )
L
0 W m
z LL 10- 0 3/
20 25 30 35 LO SNR, dB Fig. 12 Theorerical evalmiion and sirnululion of RMSE of nieun ROA in Example 3
and
5
10
15
x simulation + theory
We want to examine the errors of estimation denoted by A $ h = J k - $k and a o k = 6 k - mk (42) where $k and dk are the co-ordinates of the kth minimum on the q c plane and & and o k are the true parameters of the kth signal. Using eqns. 4 0 4 2 we can expand the partial derivatives off in eqn. 28 evaluated at ( q k , c k ) in a Taylor expansion and equate them to zero. i.e.
-
Assuming AFk and Aok are small, we can ignore the terms involving the second and higher order of A& and A o , in eqn. 43, so that we can write
(44) where
and therefore the covariance matrix of A& and A o ~ can be expressed as
Substituting f from eqn. 28 into eqns 44 and 45, and using Theorems 1 and 2, we obtain 1211 expressions for E[Aa and cov[Afl in terms of C,, Y, C, and R,. The performance can then be easily evaluated numerically. Fig. 12 shows the comparison between the theoretical evaluation of the performance and the simulation results of the DISPARE method. The comparison is based on the same scenario as that in Example 3. The two curves agree very well in the high SNR region. When SNR is low, the second- and higher-order terms 8
6
Conclusion
In this paper we have examined the problem of estimating the mean DOA of spatially distributed signal sources in array processing. From the signal model we arrive at the conclusion that even though the correlation matrix for the kth signal is of full rank, the signal energy usually concentrates in the first mOkeigenvalues, where mOkdepends largely on the spread of the signal. The remaining dimensions of the matrix can be regarded to be in the quasi-noise subspace. Assuming that the number of signals is known, and that each dispersed signal is made up of independent point sources, we then develop the DISPARE method based on the weighted projection of the eigenvectors of the signal covariance matrix onto the estimated quasi-noise subspace. The method is established based on the assumption that the distribution shape is known. However, it has been observed that, in general, if the signal is assumed to be uniformly distributed, the estimation accuracy hardly deteriorates even if the true distribution is otherwise. The problem of the correct partition of the data space into the quasi-signal and noise subspaces, however, still remains if the spread of the signal is not known. Fortunately, it is found that, for a given number of sensors, the signal energy is dominated (taking 95% of the total energy as a threshold) by mOk eigenvalues for a large practical range of signal spread. The procedure of applying the DISPARE algorithm is then established based on these observations. Simulations show that the DISPARE method is highly effective and yields performance much superior to MUSIC if the signal spread is significant. A performance analysis of the method has also been presented. The DISPARE method assumes that each dispersed signal is composed of independent point sources. In some applications this assumption may not be valid. A modified model for correlated constituent point sources can be established [18] and an extension of the DISPARE method to the model is possible [21]. However, in such cases the correlation between the point sources has to be known to obtain the estimates of the DOA. Also, other weighting functions have been used in the projection resulting in various other algorithms. But IEE Proc -Radar, Sonar Navig , Vol 143, No 1, February 1996
the DISPARE algorithm presented in this paper appears to be most robust and most accurate [21]. We can therefore conclude that under the conditions when the assumptions are valid, the DISPARE algorithm is an effective method to estimate the Of dispersed signals. 7
References
1 CAPON, J.: ‘High resolution frequency wavenumber spectrum analysis’, Proc. ZEEE, 1969, 57, pp. 1408-1418 2 PISARENKO, V.F.: ‘The retrieval of harmonics from a covariance function’, Geophys. J. Roy. Astron. Soc., 1973, 33, pp. 347366 3 JOHNSON, D.H., and DEGRAFF, S.R.: ‘Improving the resolution of bearing in passive sonar arrays by eigenvalue analysis’, ZEEE Trans., 1982, ASSP-30, (4), pp. 638-647 4 JOHNSON, D.H.: ‘The application of spectral estimation methods to bearing estimation problems’, Proc. ZEEE, 1982, 70, pp. 1018-1028 S TUFTS, D.W., and KUMARESAN, R.: ‘Estimation of frequencies of multiple sinusoids: Making linear predictions perform like maximum likelihood’, Proc. ZEEE, 1982, 70, (9), pp. 975-989 6 KUMARESAN, R., and TUFTS, D.W.: ‘Estimating the angles of arrival of multiple plane waves’, ZEEE Trans. Acrusp. Electron. Syst., 1983, 19, pp. 134-139 7 BIENVENUE, G., and KOPP, L.: ‘Optimality of high resolution array processing using the eigensystem approach’, ZEEE Trans., 1983, ASSP-31, (3,pp. 1235-1247 8 SCHMIDT, R.O.: ‘Multiple emitter location and signal parameter estimation’, ZEEE Trans., 1986, AP-34, pp. 276-280 9 BRESLER, Y., and MACOVSKI, A.: ‘Exact maximum-likelihood parameter estimation of superimposed exponential signals in,noise’, ZEEE Trans., 1986, ASSP-34, ( S ) , pp. 1081-1089 10 BOHME, J.F.: ‘Estimation of spectral parameters of correlated signals in wavefields’, Signal Process., 1986, 11, pp. 329-337 11 ROY, R., and KAILATH, T.: ‘ESPRIT Estimation of signal parameters via rotational invariance techniques’, ZEEE Trans., 1989, ASSP-37, (7), pp. 984995 12 STOICA, P., and SHARMAN, K.C.: ‘Novel eigenanalysis method for direction estimation’, ZEE Proc. F, 1990, 137, (l), pp. 19-26 13 STOICA, P., and SHARMAN, K.C.: ‘Maximum likelihood methods for direction-of-arrival estimation’, IEEE Trans., 1990, ASSP-38, (7), pp. 1132-1 143 14 VIBERG, M., and OTTERSTEN, B.: ‘Sensor array processing based on subspace fitting’, ZEEE Trans. Signal Process., 1991, 39, (5), pp. 1110-1 121 -
IEE Proc-Radar, Sonar Navig., Vol. 143, No. 1, February 1996
15 VIBERG, M., OTTERSTEN, B., and KAILATH, T.: ‘Detection and estimation in sensor arrays using weighted subspace fitting’, IEEE Trans. Signal Process., 1991, 39, (1 I), pp. 2436-2449 16 KLEIVIM, R.: ‘Use of generalized resolution method to locate sources in random dispersive media’, ZEE Proc. F, 1980, 127, (I), pp. 3 4 4 0 17 KLEpdM, R,: ‘Low-error bearing estimation in shallow water’, IEEE Trans. Aerosp Klectron. Syst., 1982, AES-18, (4), pp. 352357 18 VALAEE, S., KABAL, P., and CHAMPAGNE, B.: ‘Localization of distributed sources’, Proc. Fourteenth GRETSI Symposium on Signal and image processing, Juan Les PIns, France, Sept. 1993 19 MEN’G,Y . , WONG, K.M., and WU, Q.: ‘Estimation of the direction of arrival of spread sources in sensor array processing’, Proc. int. conf. on Sig. Proc., Beijing, China, Oct. 1993 20 WONG, K.M.: ‘Some interesting problems in sensor array processing’, Keynote address, Internat. Conf. on Sig. Proc., Beijing, China, Oct. 1993 21 MENG. Y.: ‘Estimation of directions of arrival of spatially dispersed signals in high resolution array processing’. PhD thesis, McMaster University, Canada, July 1995 22 MENG, Y., WONG, K.M., and WU, Q.: ‘DOA estimation of point and scattered sources Vec-MUSIC’, Proc. IEEE 7th Workshop on Statistical signal and array processing, Quebec City, Canadla, June 1994 23 BECK-MA”, P., and SPIZZICHINO, A.: ‘The scattering of electromagnetic waves from rough surfaces’ (Pergamon Press, 1963) 24 ISHIMARU, A.: ‘Wave propagation and scattering in random media.,vol. 1’ (Academic Press, 1978) 25 JENKINS, R.W.: ‘A field-aligned scattering model for highlatutude propagation’. Tech. Memo. DRL/TM095/92, Comm. Research Center, Ottawa, Dec. 1992 26 JENKINS, R.W.: ‘A simulation study of H F direction finding in the presence of F-region scattering and sporadic-E. Report 93004, Comm. Research Center, Ottawa, 1993 27 SLEPIAN, D.: ‘Prolate spheroidal wave function, Fourier analysis and uncertainty - V: The discrete case’, Bell Syst. Tech. J., 1978, pp. 1371-1430 28 JANTTI, T-P.: ‘The influence of extended sources on the theoretical performance of the MUSIC and ESPRIT methods: Narrowband sources’, Proc. ZCASSPPZ, Mar. 1992, Vol. 2, pp. 429-432 29 WU, Q.,and FUHRMANN, D.R.: ‘A parametric method for determining the number of signals in narrow-band direction-finding’, ZEEE Trans. Signal Process., 1991, 39, (S), pp. 1848-1857 30 WONG, K.M., WU, Q., and STOICA, P.: ‘Generalized correlation decomposition applied to array processing in unknown noise’, in HAYKIN, S.: ‘Advances in spectrum analysis and array processing, vol. 111’ (Prentice Hall, 199S), Chap. 6 -
9
