Multiple Sparse Measurement Gradient Reconstruction Algorithm for DOA Estimation in Compressed Sensing

A novel direction of arrival (DOA) estimation method in compressed sensing (CS) is proposed, in which the DOA estimation problem is cast as the joint sparse reconstruction from multiple measurement vectors (MMV). The proposed method is derived through transforming quadratically constrained linear programming (QCLP) into unconstrained convex optimization which overcomes the drawback that l 1 -norm is nondifferentiable when sparse sources are reconstructed by minimizing l 1 -norm. The convergence rate and estimation performance of the proposedmethod can be significantly improved, since the steepest descent step and Barzilai-Borwein step are alternately used as the search step in the unconstrained convex optimization. The proposed method can obtain satisfactory performance especially in these scenarios with low signal to noise ratio (SNR), small number of snapshots, or coherent sources. Simulation results show the superior performance of the proposed method as compared with existing methods.


Introduction
Direction of arrival (DOA) estimation of multiple narrowband sources is an important research topic in array signal processing.It has been extensively studied in acoustic source localization, radar, and medical imaging [1][2][3].Many effective DOA estimation algorithms have been proposed and developed, which mainly include beamforming algorithms such as MVDR [4] and subspace-based algorithms such as MUSIC [5].To obtain preferable estimation performance, the Nyquist sampling theorem must be used in these conventional methods of data acquisition.However, high-speed sampling rate can impose so enormous pressure on capturing and storing data that requirements on both hardware and software are increased.Moreover, these methods suffer from serious performance degradation in these scenarios with low signal to noise ratio (SNR), small number of snapshots, or coherent sources.
Recently, many applications involving compressed sensing (CS) [6][7][8], especially DOA estimation, have been attracting tremendous research interest in the signal processing.CS is an emerging area, and it can not only capture and store compressed or sparse sources simultaneously at a rate much lower than the Nyquist sampling rate but also reconstruct original sources using nonadaptive linear projection measurements onto a suitable measurement matrix, which satisfies the restricted isometry property (RIP) [9][10][11].The sparse reconstruction aims to find the support which is shared by the unknown sparse vectors from multiple measurement vectors (MMV).The support denotes the indices of the nonzero elements in the unknown sparse vectors.
CS has been widely applied to DOA estimation since sources are sparse in the spatial domain which results from the fact that there are much fewer true sources directions than all potential directions.Stoica et al. [12] proposed a sparse iterative covariance-based estimation method (SPICE) for array processing which is semiparametric estimation method and can avoid parameter selection.Hyder and Mahata [13] proposed an alternative strategy called joint  0 approximation (JLZA-DOA) algorithm based on spatial sparsity, which can resolve closely spaced and high correlated sources even if the number of sources is unknown.Figueiredo et al. [14] proposed gradient projection algorithm to solve the boundconstrained quadratic programming formulation.Although it can be simply implemented, the regularization parameter is selected with difficulty and it is only suitable for the single measurement vector (SMV) model, which limits its practical engineering application.
In this paper, we propose a novel multiple sparse measurement gradient reconstruction method called MSMGR for DOA estimation in CS.The method is derived through transforming quadratically constrained linear programming (QCLP) into unconstrained convex optimization to overcome the drawback that  1 -norm is nondifferentiable when minimizing the  1 -norm for the sparse reconstruction.The search steepest descent step [15] and Barzilai-Borwein step [16] are alternately used as the search step to improve the convergence rate and estimation performance significantly.Furthermore, the singular value decomposition (SVD) is incorporated into the proposed method to reduce the computational complexity and the sensitivity to the noise.The proposed method is suitable for both SMV and MMV, and it has higher estimation accuracy and resolution than existing methods especially in these scenarios with low SNR, small number of snapshots, or coherent sources.Simulation results show the superior performance of the proposed method as compared with existing methods.

Problem Formulation
Consider  narrowband far-field sources, s (),  = 1, 2, . . ., , impinging on the sensor array consisting of  omnidirectional sensors from different directions, θ ,  = 1, 2, . . ., .The array observation model at th snapshot can be formulated as where n() ∈ C ×1 is a complex Gaussian white noise vector with zero mean and covariance matrix  2 I and a(  ) is the  × 1 steering vector of the source from the direction   .
Although the DOA estimation based on the single snapshot, which is a typical SMV model, has its value, the number of snapshots is larger than one in the most practical applications.Correspondingly, the multiple snapshots model is a typical MMV model.In order to cast the DOA estimation as a sparse reconstruction, let {  }  =1 denote a fine enough grid which covers the entire spatial domain where there are  ( ≫ ) potential directions of the sources so that the true directions { θ }  =1 are aligned or are close to the grids.It means that there exists   1 ,   2 , . . .,    being equal to θ1 , θ2 , . . ., θ , respectively.Thus, we have The multiple snapshots model can be written as the following sparse form: where  is the number of snapshots and Α = [a( 1 ) a( 2 ) ⋅ ⋅ ⋅ a(  )] is the  ×  array manifold matrix corresponding to all the potential directions which is also defined as an overcomplete dictionary in CS.
is the  × 1 sparse vector with  nonzero elements at positions corresponding to the true directions and zero elements at the remaining  −  positions, where [⋅]  denotes the transpose operation.Hence, the matrix S = [s(1) s(2) ⋅ ⋅ ⋅ s()] ∈ C × has  nonzero rows, that is, row -sparse, since {s()}  =1 share the common support.Obviously, the DOA estimation problem of multiple snapshots is that of identifying the row support of the unknown matrix S from the matrix Y ∈ C × which is given by with sensing matrix Θ, noise matrix N, and a common measurement matrix Φ of the size  ×  with  <  where  is the number of nonadaptive linear projection measurements.
It is well known that sparse sources can be reconstructed by solving the  0 -norm minimization problem.However, the optimization problem is nonconvex and the optimization method is both numerically unstable and computationally unaccepted [17].Then, this problem is transformed into the  1 -norm minimization problem [18] so that we can accurately reconstruct the matrix S by solving the following QCLP problem: min S      s (2) where s (2) is the  × 1 unknown sparse vector and ‖ ⋅ ‖ 2 represents the Frobenius norm of matrices or the Euclidean norm of vectors.The th entry of s (2) is equal to Euclidean norm of the th row of S; that is, s (2)   = ‖S(, :)‖ 2 .

DOA Estimation
The SVD is employed on the matrix Y to reduce the computational complexity and the sensitivity to the noise.Hence, we have where V is the orthogonal matrix and [⋅]  denotes the conjugate transpose operation.U  and U  denote the signal subspace and noise subspace, respectively.The  eigenvalues of the matrix U are arranged from the largest to the smallest; that is, As one may note, the dimension of the matrix Y is reduced from  ×  to  ×  which can significantly reduce the computational complexity and the sensitivity to the noise especially in the scenario with large number of snapshots.The essence of the dimension reduction is to keep the signal subspace and discard the noise subspace.By using the SVD decomposition, (5) can be rewritten as the following form: where s (2)    is also a sparse vector and shares the same support with s (2) .To overcome the drawback that  1 -norm is nondifferentiable by minimizing the  1 -norm for solving the sparse sources, ( 8) is transformed into the unconstrained convex optimization by using Lagrange multiplication [19,20]: where  is nonnegative and it is called regularized factor that can serve as a tradeoff between the ability of suppressing noise and source sparsity.Note that the search path, which is obtained by projecting the negative gradient of the objective function in (9) onto the feasible set, cannot perform a backtracking line search well.Therefore, we adopt the  1norm of matrices to change the search direction which results in where ‖ ⋅ ‖  1 denotes the  1 -norm of matrices.A detailed derivation process of using MSMGR for DOA estimation is shown as follows.Assume that r Γ −1 denotes the residual of the th iteration, Γ  denotes the support of the th iteration, Θ Γ  denotes the submatrix of Θ with columns indexed by Γ  , and ŜΓ  denotes the reconstructed source after the th iteration.Therefore, the objective function of the th iteration can be written as The purpose of the current iteration is to find the sparse reconstructed source ŜΓ  which can minimize the objective function ( ŜΓ  ); that is, the residual is minimum after the current iteration.The expansion of (11) can be expressed as where (, ) denotes the element of the th row and th column of the matrix ŜΓ  .With further derivation, the minimum of the objective function ( 12) is equal to the maximum of ( 13) which can be given as the following form: Based on the properties of matrix trace, ( 13) can be further simplified to Then, the negative gradient is obtained by the partial derivative of ( 14) with respect to ŜΓ  which is given by where D is referred to as polarity matrix that can judge the polarity of nonzero elements: where (, ) denotes the element of the th row and th column of the matrix D. Since the conventional search step is too small based on the orthogonality, the steepest descent step and Barzilai-Borwein step are alternately exploited as the search step   in order to improve the convergence rate and estimation performance.Then, we have where  SD and  BB are the steepest descent step and Barzilai-Borwein step, respectively.The specific steps of MSMGR are given as follows.
Step 3. Calculate the negative gradient d  and the search step   in terms of ( 15) and ( 17), respectively.
Step 4. Update the constructed source ŜΓ  = ŜΓ −1 + d  ⋅   in terms of the negative gradient and the search step.
Step 5. Update the polarity matrix where O 1× is the zero matrix of the size 1 × .
Step 6. Update the residual r  = r −1 − Θ Γ  ⋅ d  ⋅   .If the residual satisfies the stopping criterion, stop the iteration; otherwise, set  =  + 1 and return to Step 1.
The core of the new method is that of updating the polarity matrix by zero-padding process in Step 5 since the dimensions of the support and corresponding submatrix are both expanded in Step 2.Moreover, zero-padding process can guarantee the precision of the DOA estimation.The spectrum of the proposed method is obtained by estimating constructed source power from all potential directions.Like other spectral-based methods, the true directions are estimated by the locations of the highest peaks of the spectrum.

Simulation Results
In this section, the superior performance of the proposed method is shown as compared with existing JLZA-DOA and MUSIC methods by several numerical simulations.Consider the spatial sources impinging on the uniform linear array (ULA) with interspacing /2, where  denotes the wavelength of source.In the ULA case, the steering vector corresponding to the source from the direction   is given by where the number of array elements is set to  = 16.In the simulation, the regularized factor can be chosen as suggested in [21] Following [22], it is easy to see that the unique minimum of (10) is the zero matrix for  > ‖Θ    ‖  ∞ .In the simulation, the average root mean square error (RMSE) of the DOA estimation is defined as the significant performance index: where  is the number of independent Monte Carlo runs and θ is the estimate of   in the th run.The resolution of the grid is closely related to the precision of the DOA estimation.
A coarse grid leads to poor precision, but a too fine grid increases computational complexity.Therefore, an adaptive grid refinement method is used for the tradeoff between precision and computational complexity.In the simulation, we set a coarse grid with 1 ∘ step in the range of −90 ∘ to 90 ∘ and make a local fine grid in the vicinity of the estimated angle.In the first simulation, we show the spatial spectra of three methods in the scenario with low SNR, small number of snapshots, and two uncorrelated sources impinging from  1 = 0 ∘ and  2 = 2 ∘ , respectively.The spatial spectra are shown in Figure 1 with SNR = 3 dB and 50 snapshots.The following conclusion can be acquired from Figure 1 that since the spatial spectra obtained by MUSIC and JLZA-DOA have only one peak, MUSIC and JLZA-DOA cannot identify the closely spaced sources accurately.However, the proposed method MSMGR has a nearly ideal spectrum and a precise estimation for the closely spaced sources.Therefore, MSMGR outperforms JLZA-DOA and MUSIC in terms of the spatial spectrum.
The RMSEs of three methods are analyzed under different conditions in the second simulation.Consider three sources impinging from  1 = 0 ∘  2 = 49.5 ∘ and  3 = 57.3∘ , respectively, where the latter two closely spaced sources are correlated and the first source is uncorrelated to them.The forward spatial smoothing method is exploited on the MUSIC called FSS-MUSIC to resolve the correlated sources.Figure 2 shows the RMSE as a function of SNR of three methods in 100 Monte Carlo runs for the fixed number of 50 snapshots whereas the RMSE versus the number of snapshots is shown in Figure 3 for the fixed SNR −6 dB in 100 Monte Carlo runs.The conclusions can be drawn from Figures 2 and  3 that the RMSE of MSMGR is smaller than those of other two methods and MSMGR has better estimation performance than the other two methods, especially in the scenarios with low SNR or small number of snapshots.The reason is that MSMGR can overcome the nondifferentiable drawback and exploit the alternate search step to improve the convergence rate and estimation performance.Moreover, the RMSE of MSMGR is close to those of the other two methods with the increase of SNR and the number of snapshots.
In Figure 4, the relation between the RMSE and angle separation of correlated sources is shown, which can illustrate the resolving capability.Let two correlated sources at directions 20 ∘ and 20 ∘ + Δ, where the step of Δ is 1 ∘ , be impinged on the ULA.It can be seen from Figure 4 that MSMGR suffers from serious performance degradation if the angle separation is 4 ∘ .However, MSMGR can still provide the most precise estimation as long as the angle separation is no less than 5 ∘ .Simulation results show that MSMGR has higher resolution than the other two methods.Finally, we compare the computation time of different methods versus number of snapshots in Table 1.Two correlated sources located at 20 ∘ and 27 ∘ impinge on the ULA.MSMGR-nosvd denotes that the SVD is not adopted in the process of employing MSMGR for DOA estimation.
It is easy to see from Table 1 that the computation time of MSMGR-nosvd is the longest and the computation time of MSMGR is longer than that of JLZA-DOA and FSS-MUSIC, but it is important to note that the performance of MSMGR is much better than that of these two methods.Moreover, it is proved that the SVD can significantly reduce the computation time.

Conclusion
In this paper, a novel MSMGR method for DOA estimation is proposed in CS.The proposed method is obtained by transforming QCLP into unconstrained convex optimization to overcome the drawback that  1 -norm is nondifferentiable when sparse sources are reconstructed by minimizing the  1 -norm.An alternate search step is used to improve the convergence rate and estimation performance.The SVD is used to reduce the computational complexity and the sensitivity to the noise.Simulation results show that MSMGR outperforms JLZA-DOA and MUSIC in terms of the spatial spectrum and has more precise estimation as well as higher resolution; in particular when SNR is low, the number of snapshots is small and sources are coherent.
where  large eigenvalues are dominant.U  and U  , respectively, consist of left singular eigenvectors corresponding to  big eigenvalues and  −  small eigenvalues.Denote Y  = U  = UΣΒ = YVΒ ∈ C × , S  = SΣΒ, and N  = NΣΒ where Β = [I  O ×(−) ] , I  is an identity matrix of the size  × , and O ×(−) is a zero matrix of the size  × ( − ) so that we have

Figure 3 :Figure 4 :
Figure 3: RMSE of the DOA estimation versus number of snapshots for −6 dB SNR.

Table 1 :
Computation time (sec) of methods.