Signal selective doa tracking for multiple moving targets

SIGNAL SELECTIVE DOA TRACKING FOR MULTIPLE MOVING TARGETS Jaume Riba Jason Goldberg * Gregori Vdzquez Miguel Angel Lagunas Department of Signal Theory and Communications, E.T.S.E. Telecomunicaci6, Campus Nord UPC. Edifici D-5, c/. Gran Capiti. s/n, 08034 Barcelona, Spain, Email: { jriba, jason,gregori ,miguel}@tsc .upc. es ABSTRACT A new algorithm for signal selective tracking of the directions-of-arrival (DOA’s) of multiple moving targets with an array of passive sensors is presented. A new method based on the principles of maximum likelihood estimation and cyclostationarity is used to generate initial angle estimates which, in turn, are refined by a Kalman filter. Source angle dynamics are used to achieve correct data association. High performance is obtained with relatively low computational complexity. 1. INTRODUCTION Tracking the DOA’s of moving sources has begun to receive increased attention [l], [2]. The problem consists of estimating the DOA’s of spatially non-stationary sources such that the so-called “measurement/target data association problem” is solved (i.e., each of the DOA estimates is correctly associated with a particular source from one increment of time to the next). In practice, the DOA tracking problem is complicated by the presence of co-channel interfering signals-not-of-interest (SNOI’s). However, a relatively new class of signal selective DOA algorithms has been shown to discriminate against undesired interfering signals and noise by exploiting the known cyclostationary properties of the signals-of-interest (SOI’S) [3]. In this paper, such properties are used to develop an interference tolerant DOA tracking algorithm. Its basis is similar in spirit to that presented in interference free case of [l]. However, unlike [l], the so-called “conditional maximum likelihood” (CML) principle [4] is employed. As will be shown, this likelihood function is highly amenable to the inclusion of signal selective DOA information via the cyclic correlation matrix. 2. PROBLEM FORMULATION Consider N spatially stationarynarrow band desired sources (SOI’s) impinging as plane waves on an array of M passive sensors in the presence of SNOI’s and additive white Gaussian noise. The A4 dimensional snapshot vector at snapshot time index m is modeled as: y(m) = A(e)s(m) + e(m). (1) ~~~ ~ *This work hasbeen supported in part by the Spanish MEC Grant SE92 AC 0890196 and EEC Contract HCM/CHRXCT930405 0-7803-3 192-3196 $5.0001996 IEEE 2559 = [e,, 02,. . . , BN]~ (where (.)T denotes transpose) is the N dimensional vector containing the DOA’s of the desired sources. A(B) is the MzN steering matrix whose columns, a(B,), p E {l,... , N}, consist of the steering vectors of the sources. s(m) is an N dimensional signal vector at time m. e(m) = i(m) +n(m) denotes the sum of the M dimensional interference and noise vectors at time m. With the infinite time averaging operation defined as (.)m = limIC.+m the spatial correlation matrix of the array snapshot vector is given by: I( R E (E [~(m)~~(m)])~ = A@)rAH(8) + Ree (2) (E [s(m)~”(m)])~, Ree z (E [e(~n)e~(m)])~ where E[.] and (.)H respectively denote the expectation operator and the vector/matrix conjugate transpose operation. The “asymmetric” cyclic correlation matrix of the array snapshot vector for cycle frequency a and time lag T is defined as [3]: R*(T) (E [y(m + r)y“(m)] e-Jznam)m (3) ryr) = (E [s(m + r)sH(m)] e-J2?ram)m (4) where I’“(r) is called the “cyclic intersignal coherence matrix,” and it is assumed that the SOI’s are cyclostationary at CY and that the noise and interference are not. In the spatially nonstationary case, it is assumed that the DOA trajectory as a function of time for the pth source is described by the following constant velocity model with update interval T,: xp(k+I) = [BP(k + I),&(k + I)]” = Fx,(k)+w,(k). (5) where i,(.) denotes angular velocity. The angle update time index k corresponds to time t = kT.,, and the state driving noise wP(k) is assumed stationary zero mean Gaussian of covariance Qp. If DOA estimates are required every T, seconds, the estimated (cyclic or ordinary) correlation matrix must also be updated every T, seconds. Provided the DOA’s change slowly with time, it is reasonable to assume that they remain approximately constant over each angle update interval, T, [l]: B(k) s~ fl(t), t E [kT,, (k + l)Ts). Thus, if B snapshots are available every T, seconds, the correlation matrix over the k’th interval can be estimated as: g(k) = (y(kB + m)yH(kB + m))B (6) = A(B)V-)AH@) where the finite time averaging operation is defined as: corresponds to time t = (kB + m)T,/B. In addition, an estimate of the cyclic correlation matrix of the array data for lag T over the time interval t E [ICT,, (IC + l)Ts) can be formed as: Ra(q IC) = (y(t~ + m + T)yH(t~ + m)e--32?ram)B(7) B-1 (.) - - I E,=,(.), and the snapshot time index tB + m 3. CML AND CYCLOSTATIONARITY This section explains how cyclostationary can be used to improve the CML estimates in the presence of interference. To simplify notation, the dependence on the angle vector e and the angle update time index k will be dropped whenever possible. The CML estimator, which models the sources as deterministic unknowns, is given by [4]: where tr(.) denotes the matrix trace operation, and PA and A# are, respectively, the projector and the pseudoinverse of the steering matrix. The CML function can be used in the context of DOA tracking via a first order Taylor series expansion of its derivative about e, a DOA angle vector known to be close to the ML DOA estimate h h e,, M - [L”( E)] -I L’( $) , 8 F5 BML. (9) where the Nz1 gradient and NxN (approximate) Hessian of the likelihood function are given, respectively as [5]: L’(8) = 2 Re {diag (A’GPiD)} (10) L”(E) = -2 Re { (DHPiD) 0 (A#6A#H)T}(11) a a Pi = (I - PA), D = [ afl,a(&), . . . , -a(Oiv)] a0,v with Re{.}, diag(.), and (.) 0 (.) denoting respectively “real part of,” “column vector formed by diagonal elements of,” and the Hadamard (i.e., element-by-element) product oper at ion. 3.1. Effect of SNOI’s on CML Estimates To understand how the presence of SNOI’s affect lem begin by rewriting (8) as follows: L(B) = tr (PAG) = tr (AFA~) h r = A#EA#~. The likelihood function is re-expressed in terms oLA, the angle dependent steering matrix, and the matrix I’ which is an angle and data dependent estimate of the intersignal coherence m2trix I’. When there is no additive noise or interference, I’ is asymptotically unbiased when (13) is IAs seen in Section 4. a can be obtained readily from the DOA estimate of the previous iteration or some function thereof. evaluated at the true DOA’s. Note that the presencz of additive white noise produces an asymptotic bias in l? (even when (13) is evaluated at the true DOA’s). Nevertheless, the CML DOA estimator remains asymptotically unbiased since tr(PAu21) = u2N does not vary with angle. Next, when interfering SNOI’s are present, good performance can be expected only if all of the following conditions are satisfied: (i) The total number of signals is less than the number of sensors. (ii) The DOA’s of all the signals (SNOI’s included) are estimated-resulting in a corresponding increase of the computational complexity of the CML estimator. (iii) Post-processing is performed to determine which DOA estimates correspond to the SOI’s. If the SNOI’s are not taken into account, the performance of the estimator can deteriorate dramatically resulting highly biased DOA estimates. In terms of (12) and (13), when interfering SNOI’s are present, ^r is biased by the additive noise as well as the presence of the interferers. It is this later source of bias which, in turn, results in an asymptotically biased estimator of e. Thus, the influence of the SNOI’s should be able to be reduced if their influence on f; can be reduced. The goal is to modify the CML estimator using the cyclostationary properties of the SOI’s so as to yield signal selective SO1 DOA estimates while avoiding the need to satisfy the conditions listed above. 3.2. An Improved Intersignal Coherence Matrix Estimate based on the Cyclic Correlation Matrix It will now be explained that for a cycle frequency ru of the SOI’s but not of the SNOI’s, the destructive influence of the SNOI’s on f; can be reduced by exploiting information from the estimated cyclic correlation matrix. Begin by denoting the SO1 signals at their respective transmitters as sp(m) = s (m - tP) where the time shifts, {tp}r=l, are due to the combined effect of transmission delay and differences in the phases of the transmitter and receiver clocks. The cyclic cross correlation between sources p and p at their respective transmitters is written as: -P If attention is restricted to the cases of uncorrelated sources and/or frequency non-selective multipath, then I’ can be approximated as [6]: M ---ra(r)TH 1 M -THra(r) 1 (15) Pa (7) Pa (TI where pa(.) = ya (T)/Y and -PP -PPI e--g2?ratl T = [ ... O ] . (16) e-3Zrrat~ 0 Since, as described earlier, the effect of additive noise and interference is nulled out in the cyclic correlation matrix (and, by implication, in the cyclic intersignal coherence matrix), (15) provides the basis for obtaining an estimate of I’ which is less sensitive to (i.e., less biased by) the SNOI’s. In practice, however, any estimate of I’ will require estimates of both T and P(T). For the case of the SOI’s being digital communications signals of the same, linear modulation type, it is further 25 60 shown in [SI that a multiple lag estimate of the intersignal coherence matrix can be implemented efficiently as: 2 Pa =A#(g)6."A#H(g)(17) h Ra=(x(m)xH(m)e-'2nnm)B, x(m) = y(m) *g(m)(18) where i(m) is a sampled version of the transmission pulse, and (.) * (.) denotes convolution. is a good initial estimate of gML (which will be seen to exist the case for the DOA tracking problem). R" is a type of multiple lag estimate of the cyclic correlation matrix formed by temporally pre-filteredsnapshot data. Thus, the multiple lag cyclic correlation matrix estimate of (18) is formed after each sensor output has undergone a simple time-domain matched filtering type operation. Lastly, particularly simple estimates of the time shift matrix of (16) can be obtained by generalizing the technique presented in [7] for the multisensor, multisource case by exploiting the fact that the diagonal elements of FnTH of (17) should be real. Specifically, the diagonal elements of the time shift matrix can be estimated as [6]: h h [*Ipp= ['~]pp/~["]pp~, PE {1,2,.*.,W. (19) 3.3. A Signal Selective CML Function An obvious obvious signal selective modification of (8) would be to use F, of (17) in (8) instead of r, the interference sensitive estimate of (13) normally used. Inserting (17) into (12) yields: h Direct maximization of this cost function would yield interference resistant DOA estimates. Unfortunately, obtaining expressions for the gradient and Hessian of (20) appears - to be difficult. However, a good alternative is to use I?, directly in the expressions for the gradient and Hessian of the original likelihood function of (8) (which assumes absence of interferers). In particular, the expected value of the gradient may be written as: E [L'(g)] = 2 Re{diag (A#(g)E[6.]PA(g)D(g)))) = 2 Re{ diag (A" (E) [A(e)I'AH(&) + Re,] WE) D(u>} M 2 Re{diag(I'AH(e)PA(g) .(E))} (21) +2 Re(diag(A'(8) R,;PA(E)D(g))} where the facts that if E M e then A#(fi) A@) z I and that the additive noise is spatially white have been used. The first and second terms in (21) specify the respective contributions due to the SOI's and the SNOI's. Approximating the desired first term by using (17), the estimate for the cyclic intersignal coherence matrix presented earlier in this section, yields: L'( E) M 2 Re{ diag (F( E) AH(&)Pi( E) D( g) ) )22) M Re{?diag(A#(g) &aHPi(g) D(&)) + ?Hdiag(A#(&)6.aPi(@D($))} where the approximations A#($ A(&) M I and A(&)A#(g) M PA@) both for M & in addition to the second approximation in (15) have been used. Similarly, the Hessian can be approximated as: L" ( 8) M -2 Re { ( DH ( 2) Pi ( 8) D ( g) ) 0 FT ( E) } . The use of (22) and (23) in (9) will be seen yield interference tolerant approximate CML DOA estimates. 4. TRACKING ALGORITHM The details of the Kalman filter based tracking method are now described. Recall the constant velocity state update equation of Eq. 5 used to model the DOA trajectory. The corresponding measurement equation is given as: (23) ~p(k) = hTXp(k) -I- ~p(k) (24) where h = [l 0IT, and the measurement noise vp(k) is zero mean Gaussian of variance rp(k). zp(k), hTxp(k), and +(k) correspond respectively to B^ML, (k), OP(k), and the estimation error associated with the estimate (9) of the DOA of the pth source at time k. If the measurement noises associated with different sources are assumed uncorrelated, the Kalman filter can be decomposed into a bank of computationally efficient scalar Kalman filters [I]. The specific equations are: kp(klk - 1) = Fkp(k - 1(k - l), (25-a) kp(klk) = k,(klk1) (25-b) + kp(k) (zp(k) - hTjt,(klk - 1)) where kp(klk - l), kp(klk), and kpl(k) denote, respectively, the one-step predicted state vector, the filtered state vector estimate, and the Kalman gain all for the pth source at time k. The final angle estimate is given by the first element of At iteration k, the estimator of (9) is evaluated at z(klkl), the predicted DOA source vector such that (25-b), can be rewritten as: jip(kIk) = gp(klk - l)+kp(k) [d~l,,(k) - dp(klk - I)] (26) % (k I = qklk - 1) - kp(k) [ (L"(i(k1k - 1)))-' L'(i(klk - I))] P where [.Ip denotes the p'th elemenl. of a vector. Note that calculation of EM, (k) comprises the main computational burden of the algorithm (e.g., a QR decomposition of the h 2561 A MxN matrix A( e(klk - 1)) in order to compute its pseudoinverse and calculation of the inverse of the (symmetric) The Kalman gain is updated by the following set of stanNXN matrix, L”(~?((LI~ - 1))). dard equations: C,(klk) = Cp(klk - 1) - kp(k)hTCp(klk - 1) qk + ilk) = FC,(~I~)F~ + Q~. It has been found that proper control of rp(k) can be achieved by [6]: where c is a constant directly related to the maximum anticipated total SNOI power. 5. RESULTS Computer simulation results now illustrate the improved performance of the new technique. Consider N = 3 moving SOI’s impinging on a uniform linear array of M = 8 sensors with half-wavelength inter-element spacing. They have initial DOA’s of &(O) = -70°, &(O) = 20°, and f?,(O) = 70’ and are all of unit power. Source three is a multipath reflection of source two. The SOPS consist of four sample per symbol QPSK modulated digital communications signals with unit roll-off square-root Nyquist pulses. Such signals exhibit cyclostationarity at cycle frequency cy = 0.25. The sensor noise power is of n2 = 1 unit variance, and B = 100 snapshots per tracker iteration are used. Fig 1 shows the performance of the ordinary CML tracker (using (9), (lo), and (11)) with the presence of a unit power BPSK SNOI with symbol rate two and a half times smaller than that of the SOI’s and initial DOA d(0) = -30’. CML estimator bias due to the interference results in severe performance degradation with a loss of data association and the tracking of the undesired SNOI instead of SO1 three (from k = 200 onward). Next, as shown in Fig. 2, the newly proposed signal selective algorithm yields greatly improved tracking performance. Lastly, Fig. 3 shows the performance of the algorithm when there are multiple (five) unit power SNOI’s each of the same form as that in Fig. 1. Here ordinary CML based tracking would fail even if the intention were to track all the sources since the total number of sources equals the number of sensors. Although the SINR for each source is nearly -8dB, the new technique is seen to perform well. 6. CONCLUSION A new DOA tracking algorithm for multiple moving targets in the presence of interferers has been presented. The technique applies the output of a modified signal selective CML estimator to a bank of post-processing scalar Kalman filters. High quality performance is obtained for crossing sources and/or fully coherent sources, in low SINR scenarios [I1 for data’records of reasonable length. REFERENCES C.R. Rao, C.R. Sastry, and B. Zhou. “Tracking the direction of arrival of multiple moving targets”. IEEE Transactions on Signal Processing, vol. 42: pp. 11331144, May 1994. C.K. Sword, M. Simaan, and E.W. Kamen. “Multiple target angle tracking using sensor outputs”. IEEE Transactzons on Aerospace and Electronzc Systems, vol. 26: pp. 367-373, March 1990. S.V. Schell and W.A. Gardner. High-resolution direction finding. In N.K. Bose and C.R. Rao, editors, Handbook of Statastzcs 10: Szgnal Processang and its Applzcatrons. North Holland, 1993. P. Stoica and A. Nehorai. “Performance study of conditional and unconditional direction-of-arrival estimation”. IEEE Transactzons on Acoustzcs, Speech and Stgnal Processzng, vol. 38: pp. 1783-1795, October 1990. M. Viberg, B. Ottersten, and T. Kailath. “Detection and estimation in sensor arrays using weighted subspace fitting”. IEEE Transactzons on Szgnal Processzng, vol. 39: pp. 2436-2449, November 1991. J.Riba, J. Goldberg, and G. Vkzquez. “Signal selective direction-of-arrival tracking for multiple moving targets”. submitted IEEE Transactions on Szgnal Processing, June 1995. M. Oerder and H. Meyr. “Digital filter and square timing recovery”. IEEE Transactaons on Communzcations, vol. 36: pp. 605-612, May 1988. r r Figure 1. CML tracker with SNOI (dotted trajectory). 0 SO 100 150 ZM) 2.50 300 350 400 450 500 k Figure 2. Signal selective DOA tracker. ... 0 50 100 150 200 2.50 300 350 400 450 500 k Figure 3. Signal selective tracker with multiple SNOI’s. 2562