Direction-Of-Arrival Estimation and Tracking Based on a Sequential Implementation of C-SPICE with an Off-Grid Model

This paper focuses on the problem of estimating and tracking time-varying direction-of-arrivals (DoAs) with an antenna array. A sequential DoA estimation method is proposed by extending the capon and sparse iterative covariance-based estimation (C-SPICE) method, which is an iterative off-grid method for estimating constant DoAs. Then, a moving average initialization technique is introduced such that the spatial spectrum information estimated in this snapshot can be utilized in the next one. In uniform linear arrays (ULAs), we replace the uniform grid in direction domain with that in a “frequency” domain, to improve estimation accuracy without additional complexity in practical applications. The validity and efficiency of the proposed methods are demonstrated through numerical experiments.


Introduction
Tracking directions of multiple moving targets by antenna arrays has widely applications in both military and civilian fields, such as surveillance, air traffic control, and wireless communication [1,2]. In recent two decades, many direction-of-arrival (DoA) tracking algorithms have been proposed. For example, the subspace-based methods (see [3] and references therein), methods for wideband signals tracking (see [4] and references therein), and methods for two-dimensional direction tracking (see [5] and references therein). Comparing with DoA estimation, tracking algorithms need not only to estimate DoAs repeatedly and sequentially over time, but also to associate DoA estimates to tracks according to their sources, which is referred to as data association. Hence, many data association methods for multiple target tracking are also proposed [6][7][8][9][10].
Recently, the research on DoA estimation has been evolved due to the developments in sparse signal reconstruction (SSR) methods and many algorithms are proposed accordingly (see [11][12][13] and references therein). However, the SSR method designed for DoA tracking is rare. Some methods [14,15] can be used for DoA tracking, but they have a relatively high computational complexity. In [16], an efficient algorithm named by tracking via low-rank plus sparse recovery (TvLSR) is proposed. Different from traditional tracking algorithms, where DoAs are estimated and tracked repeatedly and sequentially over snapshots (or time windows), TvLSR estimates tracks of sources over all snapshots once for all. Hence, it is actually a block trajectory estimation method. A Bayesian sparse-plus-low-rank matrix decomposition (BSPLR) method is also proposed to improve the estimation accuracy of TvSLR [17]. These methods are based on on-grid models, whose estimation accuracy is limited by size of the grid [11]. Although an off-grid model can be incorporated, BSPLR is computationally expensive.
A capon and sparse iterative covariance-based estimation (C-SPICE) method, which estimates DoAs of stable sources based on an off-grid model, is proposed in [11] for a better estimation accuracy and computational complexity tradeoff. In this paper, we extend C-SPICE to track DoAs of moving targets by exploiting the iterative structures in it. Specifically, an efficient sequential implementation of C-SPICE is introduced to estimate DoAs in each snapshot. To save number of iterations, a moving average technique is proposed such that the spatial spectrum information estimated in this snapshot can be utilized to construct a good initial point for iterations in the next one. In uniform linear arrays (ULAs), we replace the uniform grid in direction domain by the grid in a "frequency" domain, whose density adapts automatically to the resolution ability of a ULA and is greater when DoA is small. Since the effective scanning range of a ULA may be no more than (−60 • , 60 • ) in some practical applications [1], the new grid can improve estimation accuracy without additional complexity.
Notations: The operators (·) T and (·) H denote transpose and conjugate transpose, respectively, and diag(x) returns a diagonal matrix with main diagonal x. a and A denote the l 2 -norm and Frobenius norm of a and A, respectively. C N×1 represents the space of N-dimensional complex column vectors. For an integer N, [N] is defined as {1, 2, · · · , N}. I represents an identity matrix and 1 N an N × 1 vector with entries 1. Re(·) denotes real part of a complex and tr(·) calculates the trace of a matrix.

Signal Model
Consider the problem of estimating DoAs of M far-field narrow-band sources based on the outputs of an N-element array. A nonparametric model of the array output can be expressed by [14]: where y(t) ∈ C N×1 and n(t) denote the snapshot and the noise at time t, respectively, {θ k } K k=1 is a uniform and discrete grid covering the possible range of DoAs, a(θ k ) ∈ C N×1 is the steering vector ofθ k , x k (t) ∈ C is the unknown signal waveform of a possible source atθ k , and M < N K. Obviously, the vector x(t) = [x 1 (t), · · · , x K (t)] T is M-sparse where a nonzero entry corresponds to a source.
In a tracking scenario, suppose that DoAs of sources vary continuously and slowly over neighbouring snapshots. Then the actual DoA of m-th source at time t, i.e., θ m (t), may not locate on but close to some grid pointθ m k , and thus θ m (t) =θ m k + m k (t). With this modeling mismatch, the signal model (1) can be modified into [18,19] where e(θ k ) is the first order derivative of a(θ) at θ =θ k , k (t) ∈ [−l , l ], and 2l denotes the distance between neighbouring grid points. Equation (2) can be expressed in a matrix form where

The C-SPICE Algorithm
C-SPICE is a grid-based off-grid sparse estimation method for constant DoA estimation [11]. It enjoys an effective tradeoff between modeling accuracy and computational complexity. Here, we introduce it briefly with the time index "(t)" omitted for simplicity. First, the following assumptions are required in deriving the method: 1.
The noise is spatially and temporally white Gaussian with covariance σI; 2.
the noise is uncorrelated with the signals; 3.
the signals are uncorrelated with each other.
Note that the third assumption is not essential, but only required during the derivation [11,14]. With these assumptions, there is where DefineR = ∑ T t=1 y(t)y H (t)/T. When T > N such thatR −1 exists, the following optimization problem can be constructed under covariance fitting criterion [11] min p,σ, Problem (10) is in general nonlinear and nonconvex (see Proposition 1 in [11]) and C-SPICE is a method to approach its optimal solution. In order to make the paper self-contained, below we provide a different and new derivation of C-SPICE.
First, it is obvious that the grid size K is usually large and thus | k | ≤ l 1, ∀k ∈ [K]. Hence, the objective function f (p, σ, ) can be simplified by letting = 0 in the first term, i.e., tr(RR −1 ). Then substituting (7) into (10) and omitting the terms independent of , we obtain min p, If p k is given, problem (11) is convex and k can be estimated by the closed-form solution In (12), it is seen that can be determined without the information of the tuple (p, σ).
The second step estimates (p, σ). Let k =ˆ k and define g k = a(θ k ) +ˆ k e(θ k ), ∀k. Then problem (10) is simplified to where Rˆ = ∑ K k=1 p k g k g H k + σI. Problem (13) can be solved by SPICE. For self-containedness, a derivation of SPICE, which is different from that of [14], is provided in the Appendix.
Finally, with the estimates of and p, DoA estimates can be obtain bŷ where k m denotes the index of the m-th highest peak in p.

A Sequential Implementation of C-SPICE
In this section, we extend C-SPICE to the scenario of estimating moving targets in the hope of deriving a fast and accurate tracking algorithm.

Sequential Implementation of C-SPICE
First, we introduce a sequential implementation of C-SPICE. Suppose that DoA estimation at time t is based on the most recent Since ln |R| is nonconvex in (p, σ, ), we approximate it by the first-order Taylor expansion at some point R 0 (which is a positive definite matrix close to the true R) and then f L (p, σ, ) by If W > N, then R 0 =R(t) is a good approximation of R and f R 0 (p, σ, ) becomes f (p, σ, ) in (10b). So DoAs can be estimated directly by using C-SPICE. Otherwise, when W < N or source signals are highly correlated,R(t) does not coincide with the definition of R in (7). In this case, we can minimize f R 0 (p, σ, ) to estimate (p, σ, ) with a given R 0 and then update R 0 with the estimated (p, σ, ), iteratively. This iteration is referred to as outer iteration, where R 0 can be initialized by I. In each outer iteration, f R 0 (p, σ, ) can be minimized by a modification of C-SPICE, in whichã (12) and the iterations updating (p, σ), i.e., Equation (A5), are modified to We call the iterations updating (p, σ) by (17) inner iteration. When outer iteration converges, DoAs can be obtained by (14). We summarize the sequential implementation of C-SPICE (SIC-SPICE) in Algorithm 1, where I 2 denotes number of outer iterations in each snapshot and I 1 the number of inner iterations. If W N, Step 12 can be simplified by replacingR 0.
Equation (17), which can be readily derived by modifyingR 0.5 (t) to Y H (t)/ √ W in the constraint of (A1).
Since SIC-SPICE is an iterative method, its convergence behavior needs further discussion. In each snapshot, there are two types of iterations. In inner iterations, the tuple (p, σ) is estimated with =ˆ . According to the derivation in Appendix, inner iterations can converge to a global optimal point (which depends onˆ ) from any initial value in the interior of feasible set. On the other hand, if function (16) can be minimized in each outer iteration, then the obtained optimal solutions reduce the value of f L (p, σ, ) monotonically (see (28)-(31) in [11]). However, C-SPICE is not guaranteed to provide a point minimizing (16). In this case, we use the solution of outer iterations that minimizes f L (p, σ, ) for DoA estimation.
The complexity of SIC-SPICE in each snapshot is O(N 3 KI 1 I 2 ) [11]. When N is given, we may pursuit an accuracy and complexity trade-off by choosing K, I 1 , and I 2 intelligently.

A Moving-Average Initialization Technique
To obtain an accurate estimation of (p, ), the required number of inner iterations in SIC-SPICE, i.e., I 1 , will be large, which means that a large processing time is required for each snapshot. This is undesirable, especially when the processing time is larger than the snapshot interval.
On the other hand, in the tracking scenario, one can assume that the variations of actual DoAs between neighbouring snapshots are small such that |θ m (t) − θ m (t − 1)| ≤ 2Bl , ∀m, t, where B is an integer. This assumption is reasonable since time interval between neighbouring snapshots is usually tiny. Therefore, the optimal (p, σ) of this snapshot is close to that of the next one, which means that outputs of this outer iteration can be a good initial point for the next. Unfortunately, these outputs cannot be used directly to initialize (p, σ), which can be readily observed by setting p k = 0 for some k. Then, according to (17), p k still equals zero even if a target has moved toθ k . To overcome this problem, we propose a moving-average (MA) technique for initializing p and σ in each outer iteration, which is, where c is a small number preventing p k and σ being zeros and B controls the size of the block dispersing amplitude of p k . Then the information obtained in previous iteration can be used and the algorithm will not converge to a trivial solution.

SIC-SPICE for ULAs
In a ULA, a steering vector can be expressed by a(θ) = 1, e j 2π λ d sin(θ) , . . . , e j 2π λ (N−1)d sin(θ) T , where λ and d denote wavelength and element spacing, respectively. Then, instead of using the uniform grid in θ-domain, i.e., {θ k } K k=1 , we use a uniform grid in the "frequency" domain, i.e., { f k } K k=1 , where f = 2d λ sin(θ). Without loss of generality, let us assume d = λ/2. Thus the domains of θ and f are (− π 2 , π 2 ) and (−1, 1), respectively. By mapping { f k } K k=1 into the θ-domain, the resulting grid is more fine around θ = 0 and more coarse as θ approaching ± π 2 . Since the resolution ability of a ULA decreases monotonically as the DoA θ approaching ± π 2 [1], the density of grid in "frequency" domain adapts to resolution of the ULA automatically. Moreover, in practical applications, the effective working scope of a ULA is usually less than (− π 2 , π 2 ) [1]. Hence, for a given grid size K, { f k } K k=1 enjoys a more fine grid and thus higher estimation accuracy for DoAs close to 0.
In addition, since the correlation matrix R can be constructed based on the result of outer iterations in SIC-SPICE, subspace methods, such as multiple signal classification algorithm (MUSIC) [21], can be used to estimate DoAs in Step 17 of Algorithm 1. The advantages of subspace methods are two folds: first, estimation accuracy of the algorithm is no longer limited by the capacity of grid; second, it may avoid missing true DoAs at nearby grids, which produce only one peak.

Simulations
In this section, we evaluate the performance of the proposed algorithms by comparing them with existing SSR-based methods. The compared algorithms include SPICE [14], Algorithm 1 (SIC-SPICE), Algorithm 1 with moving average initialization (SIC-SPICE+), Algorithm 1 with moving average initialization based on uniform grid in "frequency" domain (SIC-SPICE++), tracking via low-rank plus sparse matrix decomposition (TvLSR) [16], and Bayesian sparse-plus-low-rank matrix decomposition method (BSPLR) [17]. For all SPICE kind of methods, we set W = 1 such that DoAs are estimated based on the latest snapshot, apply MUSIC algorithm to estimate DoAs, and use B = 1 in the moving average technique (18). In SPICE, the maximum number of iterations is 100 for the first snapshot and 6 for the others. In SIC-SPICE methods, we set I 1 = 100 for the first outer iteration in the first snapshot, I 1 = 3 for the others, and I 2 = 2. All the parameters in TvLSR and BSPLR are chosen according to those in [16,17], respectively, and the maximum number of iterations for both is 150.
In all simulations, ULAs with half-wavelength spaced sensors is used, and the number of the signal sources in space is M = 3. Two kinds of trajectories are considered: • Case-1: The first two sources move from 35 • to 5 • and from 15 • to −5 • , respectively, over 100 snapshots. The third one moves following the function θ 3 (t) = −25 • -8 • sin(0.05t).

•
Case-2: The trajectories of the first two sources are the same as those in Case-1. The third one moves according to θ 3 (t) = −5 + 0.25t + 5 sin(0.05t).
In both cases, the signal of the m-th source is generated by x m (t) = ς m exp(−jυ m (t)), where υ m (t) is uniformly distributed in [0, 2π) and amplitudes have the relationship ς 2 m = 3ς 2 m−1 . The signal to noise ratio (SNR) is defined by 10 log 10 (ς 2 1 /σ). A prerequisite for multiple-target tracking is associating DoA estimates between neighbouring snapshots. This problem is nontrivial and has been discussed in many literatures (see [6][7][8][9][10] and literatures there in). Hence, we assume that the estimates are well associated in simulations. For example, we may assume that powers of different sources are different such that DoA estimates can be associated according to signal power estimation. Alternatively, in ULAs, DoA estimates can be updated by using Newton's method to the MUSIC spectrum with estimates in the previous snapshot as initial values. Then, the estimates of neighbouring snapshots are associated automatically.
In simulations, we assume that sources have different powers, i.e., ς 2 m = 3ς 2 m−1 , such that data association can be easily achieved based on power estimation. This assumption is not necessary when a more powerful data association method, such as those proposed in [6][7][8][9][10] are used.
In first experiment, the trajectories estimated by proposed methods are illustrated. The parameter settings are N = 20, K = 10 N, and SNR = 10 dB. Figure 1 shows estimated trajectories of SIC-SPICE+ in both cases. The estimated trajectories of SIC-SPICE and SIC-SPICE++ are similar and omitted here. In Figure 1a, we can see that SIC-SPICE+ can estimate DoAs correctly. In Figure 1b, fails may happen at cross points. However, by careful observation, it can be found that these fails are caused by the naive data association technique, which does not work at the crossing point. Since data association is out of the scope of this paper, we will focus on Case-1 in the following simulations. The second experiment compares DoA estimation performance of the proposed methods and SPICE in Case-1. The root mean square error (RMSE) versus snapshot index obtained through 50 trials are shown in Figure 2. The RMSE performance of the second source is similar to that of the first one and omitted. We can see that SIC-SPICE outperforms SPICE slightly. The reason is that DoA estimations of both method are not limited by grid points when MUSIC is applied, and off-grid model of SIC-SPICE is more accurate. Moreover, SIC-SPICE+ and SIC-SPICE++ can reduce RMSE significantly. The reason is that SPICE and SIC-SPICE cannot converge to the optimal solution with a small number of iterations, while methods based on moving average have a better initial point and thus converge faster to the optimal point.
The third experiment evaluates performance of the proposed methods with different grid sizes. ULAs with N = 20 and N = 30 are considered, respectively. The RMSE curves of different algorithms are shown in Figure 3, where each point is based on results of 50 trials. Again, RMSE curve of the second source is omitted. It is seen that TvSLR performs worst. The failure of TvSLR is caused by two reasons: when the SNR is not high enough, e.g., for the first source, TvSLR may suffer from false DoA estimates; when the SNR is large, e.g., for the third source, its RMSE is lower bounded by a number on the same order of 2l , for the reason that its DoA estimates are restricted by grid points. BSPLR can overcome false estimation, but its estimates are still restricted by the grid. It is seen in Figure 3a that the performance improvement caused by moving average initialization becomes larger as the increase of K. The reason is that moving average can provide a more precise information on spatial spectrum based on estimates of the previous outer iteration. Furthermore, SIC-SPICE++ is attractive when K is small. This is because the grid of SIC-SPICE++ is more fine around zero and more coarse around ±90 • , and the range of true DoAs is around zero. Although performance loss may be induced if DoAs are large, the effective working area of a ULA is usually in (−60 • , 60 • ) in many practical applications [1]. In Figure 3b, when K is too large such that the assumption |θ m (t) − θ m (t − 1)| ≤ 2Bl with B = 1 is not satisfied, moving average leads to performance loss in SIC-SPICE+ and SIC-SPICE++. At this point, B > 1 is required, such that this assumption can be satisfied. In addition, the average running times of different algorithms are plotted in Figure 4. It is seen that BSPLR is time consuming and the average running time of the others are similar. It should be mentioned that when N and K are both large, e.g., N = 90 and K = 20 N (with 2l = 0.1 • ) as in [16], TvSLR is faster than SPICE. Luckily, according to Figure 3, the latter does not require such a large K to achieve an accuracy on the order of 0.1 • in a ULA. For other kinds of arrays, such as uniform circular arrays, SIC-SPICE methods, which are based on off-grid models, can achieve a high estimation accuracy with a relatively small K.   The last simulation assesses RMSE of algorithms versus SNR in Case-1 when N = 20 and K = 10 N. BSPLR is not considered for the heavy workload. RMSE curves of the first and third sources are illustrated in Figure 5. It is seen that TvSLR is lower bounded, for the reason of error DoA estimation as well as on grid constraint. SPICE is better than SIC-SPICE at the low SNR region, the reason is that SIC-SPICE is not converged. Note that the number of inner iterations is I 1 = 3 in SIC-SPICE and 6 in SPICE. As SNR increases, SIC-SPICE outperforms SPICE, since the estimation error is dominated by modeling errors at the moment. SIC-SPICE+ and SIC-SPICE++ outperform the others, especially at the high SNR region.

Conclusions
In this paper, DoA estimation and tracking algorithms are discussed in the framework of SSR. The proposed methods are extensions and improvements of C-SPICE, which is an efficient off-grid method for estimating constant DoAs. Numerical experiments demonstrate the validity and efficiency of the proposed methods. An accurate data association algorithm under the framework of SSR can be investigated in the future.
Since problem (A3) is convex, its optimal solution can be attained from any initial point in the interior of feasible set, one of which is provided by the periodogram (PER) method [14], i.e., p k = R 0.5 g k