Abstract
A novel direction of arrival (DOA) estimation method in compressed sensing (CS) is presented, in which DOA estimation is considered as the joint sparse recovery from multiple measurement vectors (MMV). The proposed method is obtained by minimizing the modified-based covariance matching criterion, which is acquired by adding penalties according to the regularization method. This minimization problem is shown to be a semidefinite program (SDP) and transformed into a constrained quadratic programming problem for reducing computational complexity which can be solved by the augmented Lagrange method. The proposed method can significantly improve the performance especially in the scenarios with low signal to noise ratio (SNR), small number of snapshots, and closely spaced correlated sources. In addition, the Cramér-Rao bound (CRB) of the proposed method is developed and the performance guarantee is given according to a version of the restricted isometry property (RIP). The effectiveness and satisfactory performance of the proposed method are illustrated by simulation results.
1. Introduction
Direction of arrival (DOA) estimation of multiple narrow-band sources has been an interesting research topic in array signal processing. Its applications span many fields including radar, communication systems, and medical imaging [1, 2]. Many effective algorithms are proposed for DOA estimation, which are mainly classified into three categories. The beamforming algorithms such as MVDR [3] and MEM [4] obtain a nonparametric spatial spectrum by optimizing the filter weights. The subspace algorithms such as MUSIC [5] and ESPRIT [6] and their derivatives [7, 8] exploit the orthogonality of signal subspace and noise subspace for DOA estimation. The subspace fitting algorithms including Maximum Likelihood (ML) [9] and Weighted Subspace Fitting (WSF) [10] solve a multidimensional nonlinear optimization problem to obtain a precise estimation, but a good initialization is required to ensure global convergence and computational complexity is very high. All these algorithms focus on two important issues: resolution (i.e., the ability to correctly resolve two closely spaced sources) and precision (i.e., the deviation from true DOA) which are considered to be the theoretical bases of evaluating certain algorithms. Many high resolution methods suffer from serious performance degradation in the scenarios with low SNR, small number of snapshots, and closely spaced corrected sources. More recently, many research applications involving compressed sensing (CS) [11], especially DOA estimation [12, 13], have become more and more popular in the signal processing. Moreover, the restricted isometry property (RIP) based on the perfect theoretical framework given with modern probability theory by Donoho [11] and Candés et al. [14] provides the performance guarantee in CS. Hence, in this paper, DOA estimation is posed as a joint sparse recovery problem where we recover jointly sparse sources from multiple measurement vectors (MMV) under CS framework.
CS is an emerging area which can break through the limit of Nyquist sampling theorem. On one hand, CS can simultaneously capture and store compressed or sparse source at a rate much lower than that of Nyquist sampling. On the other hand, it can recover the original source using nonadaptive linear projection measurements onto a suitable measurement matrix which satisfies the RIP. In CS, the joint sparse recovery aims to find a common support shared by unknown sparse vectors from MMV, which is obtained by the measurement matrix. Support denotes the indices of the nonzero elements in the unknown sparse vectors. A sparse solution can be obtained as long as the support is determined.
CS theory has been widely applied to DOA estimation according to the source sparsity, which results from the fact that there are much fewer source directions than all potential directions in the spatial domain. The DOA estimation methods in CS are attractive since they have much better estimation performance than conventional estimation methods. In [15], Gürbüz et al. firstly formulate the DOA estimation problem under CS framework in the time domain. Wang et al. [16] propose a new architecture to estimate DOA by exploiting compressive sampling in the spatial domain. Stoica et al. [17] make full use of covariance matching criterion and present a semiparametric/sparse estimation method and its derivative called LIKES [18] for the separable model. Malioutov et al. [19] present the -SVD algorithm for DOA estimation which combines the SVD of the data matrix with a sparse recovery method based on -norm minimization. A new class of subspace-base algorithms such as compressive MUSIC (CS-MUSIC) [20] and subspace-augmented MUSIC (SA-MUSIC) [21] is proposed in recent years.
The RIP and various modified versions of it have been used as a foundation of performance guarantee [21–24] for the joint sparse recovery. The performance guarantee of MUSIC based on joint sparse recovery is given for identifying the unique support in a favorable case [25]. However, in the unfavorable case where the number of measurement vectors is smaller than the sparsity or the covariance matrix tends to lose rank due to correlated sources, performance guarantee fails. Lee et al. [21] propose a new performance guarantee in terms of a version of the RIP under such unfavorable conditions. The performance guarantees of other methods such as greedy algorithms and convex relaxation have been developed to find the sparse solution in [23, 24]. However, the guarantees of these methods cannot be simply extended to the MMV case to obtain a better bound for the sparse solution.
In this paper, we propose a novel augmented Lagrange based on modified covariance matching criterion method called AULMC for DOA estimation in CS. This method utilizes the minimization of the modified covariance matching criterion which is acquired by using regularization method to add penalties in order to obtain the stable sparse parameter estimation, especially in the low sparsity case. This minimization problem is shown to be a semidefinite program (SDP) and transformed into the constrained quadratic programming problem for the sake of reducing computational complexity. The augmented Lagrange function is formed to solve this problem by the use of augmented Lagrange method. AULMC has a number of advantages over the other methods. For example, it provides more precise estimation and higher resolution in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources, and it does not need priori knowledge about the number of sources or to choose the regularization parameter of the -optimization framework which is very difficult to select in the DOA estimation. In addition, we give a detailed derivation process of the closed-form expression for the Cramér-Rao bound (CRB) of the new method and discuss an explicit condition that guarantees performance of the new method. This performance guarantee is given in terms of weaker version of the RIP which is referred to as weak-1 RIP. Simulation results illustrate the performance of the proposed method.
It is worth noting that covariance matching criterion has been used for DOA estimation [26]. In [26], a sparse iterative covariance-base estimation method, abbreviated as SPICE, is proposed. Our approach is different from SPICE because it utilizes modified covariance matching criterion which can guarantee the stability of solution even if the source sparsity is rather low. In the future work, we will focus on the application of our approach to data-driven design [27–30]. Now we briefly summarize the contributions of this work as follows.(i)A modified covariance matching criterion is proposed by adding penalties according to the regularization method.(ii)The original SDP problem is transformed into the constrained quadratic programming problem. The motivation to transform the problem is that it can reduce computational complexity.(iii)Augmented Lagrange based on modified covariance matching criterion method is devised to solve the resulting programming problem.(iv)The CRB and performance guarantee of the new method are given in detail.
The rest of this paper is organized as follows. In Section 2, we formulate the DOA estimation problem. Section 3 presents a modified covariance matching criterion. A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is described in detail in Section 4; the performance of which is analyzed in Section 5. The performance of the proposed method is evaluated by simulations in Section 6. Conclusions are provided in Section 7.
2. Problem Formulation
Consider narrow-band sources impinging on a sensor array that consisted of omnidirectional sensors from far field. At time instant , the array received source can be given by where denotes a noise term, is the unknown direction, of the th source where denotes the entire spatial range and is the steering vector. Although the DOA estimation of the single snapshot, which is a typical single measurement vector (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 MMV model.
Let denote a fine grid which covers where there exist potential directions of the sources so that the true directions are aligned or are close to the grids. This means that if are equal to , respectively, we have
Hence, the multiple snapshots model can be written as the following sparse form: where is the manifold matrix corresponding to all the potential directions which is also referred to as an overcomplete dictionary in CS. is a -sparse vector since it has at most nonzero elements in elements, and is defined as sparsity where the operator denotes transport. are jointly -sparse if they share a common support. Hence, the matrix has no more than nonzero rows in order to be called row -sparse. The MMV problem is that of identifying the row support of the unknown matrix from the matrix that consisted of MMV which is given by with a common measurement matrix of the size with where is the number of nonadaptive linear projection measurements, such as random Gaussian matrix or random partial Fourier matrix, and noise matrix .
3. Modified Covariance Matching Criterion
In this section, the modified covariance matching criterion is developed in the CS scenario. A conventional covariance matrix of the compressed measurement source with the size is given by where is a covariance matrix of the sparse source whose off-diagonal elements denote the source correlation and diagonal elements denote the source powers. Since the powers are one-to-one corresponding to all the potential directions and our focus is on the source angle parameter estimation, can be reduced to a diagonal matrix . According to (5), the measurement matrix can change the covariance matrix of the noise to render the noise colored even if the noise is white. Therefore, this adverse factor must be considered before recovering jointly sparse sources. In addition, the operators and denote conjugate transpose and expectation, respectively.
Since is a positive definite Hermitian matrix, a prewhitened method is given by the Cholesky decomposition. Let be the Cholesky factor that satisfies where is an upper triangular matrix of positive line. Hence, a prewhitened process is implemented by multiplying by in order to obtain a pure source whose covariance matrix is written as where and is an identity matrix of the size . It is important to note that the unknown covariance matrix of the noise is transformed into the known identity matrix after the prewhitened process. Therefore, the prewhitened method improves the robustness to the noise. Then, the covariance matrix of the compressed measurement source denoising is realized by
The parameter can be estimated by a class of the covariance matching estimation techniques (COMET) based on covariance matching criterion [31]. This parameter estimation method is well known to be a large-sample approximation of ML method and provide a more attractive solution than ML estimator.
The principle of COMET is that of using the right data to minimize its data model by the weighted least-squares (WLS) method. However, the lower the source sparsity is, the more likely it is to be ill-posed for the covariance matrix estimation error meaning that the optimal solution obtained directly by minimizing the conventional covariance matching criterion is instable. Hence, we employ regularization method to add penalties in order to sufficiently exploit prior knowledge to reduce the solution space for determining the stable optimal solution. The modified covariance matching criterion is proposed as the following form: where parameter controls the solution smoothness (guarantee precision), parameter controls the solution scale (guarantee sparsity), the inverses of the matrix , the sample covariance matrix , and are assumed to exist, and the matrix is the inverse of the covariance matrix of the residuals, . Since is equal to as , it follows from [32] that satisfies the asymptotic normal distribution with mean zero and variance . In addition, the operators , , and denote matrix trace, column stacking operation, and Kronecker product, respectively. Then, the matrix can be given by
We consider the normalized data model in (9) and choose the regularization parameters based on perceptual criterion [33]. By substituting (10) into (9), we have
By the properties of , , and , the data model can be further simplified to where denotes the Frobenius norm for matrices or the Euclidean norm for vectors. The data model in (12) is considered to be the modified covariance matching criterion. It can be seen from (12) that a positive definite Hermitian matrix is added to the covariance matrix estimation error by applying penalties according to regularization method in order to guarantee the stability of solution.
4. DOA Estimation
In this section, we will utilize the minimization of the modified covariance matching criterion to estimate parameter . Let be the optimal solution of in the structure model of (8). Our goal is to utilize the modified covariance matching criterion for an estimate that is asymptotically equal to . By the properties of the trace and the Hermitian matrix, a derivation process is shown as follows, where we omit the dependence on for notational convenience:
Due to we can deduce that the minimization of is equal to the minimization of :
Then, we will demonstrate that the minimization of in (16) with respect to is an SDP. To prove this fact, One assumes
The in (16) can be rewritten as
By the Schur complement, let be auxiliary variables satisfying
The equivalent form of this minimization problem is expressed as
It is easy to see that (20) is an SDP [34]. Many software packages can solve an SDP, but solving (16) as an SDP is not a good choice because SDP solvers have generally rather high computational complexity for the values of , , and in the DOA estimation. To solve it effectively, we transform it into the constrained quadratic programming problem for reducing computational complexity, as described next.
Since a consistent estimation of is given by (15), we can reformulate the minimization of in (16) as the following constrained minimization by the Schur inequality of the trace: where is a diagonal matrix of and is a -dimensional vector with . By substituting (8) into (21), the objective function in (21) can be rewritten as
Based on the following equation: where denotes the Schur-Hadamard product, and are both matrices, , and , the objective function (22) can be written as where and . Hence, the minimization problem in (21) is transformed into the following form: which is a typical constrained quadratic programming problem.
It is well known that an important class of methods for solving the constrained quadratic programming problem is to form the auxiliary function. To solve (25), we form the following augmented Lagrange function with respect to (25): where is the asymptotical solution of the Lagrange multiplier of (25) and is a penalty factor.
By setting the gradient and the Hessian matrix of (26) with respect to to zero, we have where and . One assumes that
Applying the Newton method, we obtain
For notational convenience, we assume that , , , , and we get the following equation by (29):
Assuming that the inverse of exists, by left-multiplying (30) by , we have
It follows from (31) that where satisfies
By substituting (32) into (30), we have Both the multiplier factor and the penalty factor are determined with difficulty for utilizing the augmented Lagrange function to solve the constrained quadratic programming problem. Hence, an updated sequence for the multiplier factor is given in terms of Proposition 1.
Proposition 1. One assumes that is the optimal solution of the problem . Then, the following equation holds: where .
Proof. See Appendix A.
It can be deduced from Proposition 1 that is referred to as the next iteration of . We apply a heuristic update sequence for the penalty factor to achieve a stable solution. If the th iterative solution is closer to the feasible region than the previous solution , the penalty factor is decreased. Inversely, we increase the penalty factor when is not closer to the feasible region.
The specific steps of solving by the augmented Lagrange method are given as follows.
Initialization: set , , , , and .
(1)Calculate and in terms of (33) and (34). If , is the KKT point of the problem (25); then stop iteration.(2)We use Armijo search method to find the maximum of satisfying (3)Set and update the multiplier factor and penalty factor, respectively: (4)Set and return to step (1).
5. Performance Analysis
5.1. Cramér-Rao Bound
In this subsection, the closed-form expression for the CRB of the proposed method with complex white Gaussian noise after the prewhitened process is illustrated. The bound of the noise variance estimation can be computed separately as (see [35]). The remaining parameters consist of an unknown vector . It is not an easy task to get the CRB of the unknown parameters. However, fortunately, [36, 37] have provided a critical inspiration for the derivation in this paper.
First, the likelihood function is given by
Thus, the log-likelihood function is
Then, the partial derivatives of (39) with respect to , and are given by where and . Following [35, 37], we can obtain the Fisher information matrix (FIM) as follows: where
It is well known that
It can be deduced from (41) and (43) that where is a matrix and denotes pseudoinverse. Note that the CRB in CS is affected not only by the conventional factors, for example, SNR, array structure and signal relative location, but also by the measurement matrix.
5.2. Performance Guarantee
In CS, the RIP has been deeply studied for the joint sparse recovery by minimizing the norm. We say that the matrix obeys the RIP of the order if there exists a constant satisfying
Therefore, all submatrices of with columns are uniformly well conditioned. The restricted isometry constant (RIC) of the order , described as , is the smallest that satisfies (45) and satisfies where is a -dimensional subsupport of and denotes the submatrix of with columns indexed by . Note that the condition satisfying the RIP is so demanding that its applications are limited. Therefore, we should make use of a new version of the RIP, which is called the weak-1 RIP [38], to control the size of the recovery error. The weak-1 RIP is given in the following form: for all supported on , where the cardinality of the set is . If the matrix satisfies the weak-1 RIP, it can be deduced that , where denotes the th largest eigenvalue of . The corresponding weak-1 RIC is given by
In this paper, when the estimation quality is imperfect, especially in the unfavorable case, the support is no longer identified by the algebraic property of . Hence, a new performance guarantee is given in terms of the weak-1 RIP in the following proposition.
Proposition 2. One assumes that and where . is an estimation that is asymptotically equal to such that for . Let and be the L-dimensional supports that consisted of the indexes of largest elements in and , respectively. If the matrix satisfies for the support can be identified.
Proof. See Appendix B.
It follows from Proposition 2 that the performance guarantee of the proposed method requires a mild condition in the unfavorable case. However, in the favorable case, the performance guarantee only requires a much milder condition, , which is an algebraic condition.
6. Simulation Results
In this section, the performance of the proposed method is illustrated by several simulation results and compared with that of CS-MUISC, SPICE, and CRB under the condition that the number of sources is unknown. We consider the spatial signal impinging on the uniform linear array (ULA) with interspacing where denotes the wavelength of source. In the ULA case, the steering vector corresponding to the DOA equal to is given by where the number of the array elements is set to be .
In the simulation, the average root mean square error (RMSE) of the DOA estimation with 50 Monte Carlo runs is defined as the significant performance index: where 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 can lead to poor precision, but a too fine grid increases computational complexity. Therefore, an adaptive grid refinement method is used to balance the tradeoff between precision and computational complexity. In the simulation, we make a coarse grid with step in the range of to and perform a local fine grid in the vicinity of locations obtained by using the coarse grid.
In the first simulation, we display the superimposed spatial spectra of three algorithms in 10 Monte Carlo runs in the scenario with low SNR, small number of snapshots, and five sources impinging from , respectively where two most closely spaced sources are correlated and the remaining sources are uncorrelated. The spatial spectra are shown in Figure 1 with 3 dB SNR and 50 snapshots. The following facts can be acquired from Figure 1 as follows: the spatial spectrum obtained by CS-MUSIC suffers from severe interference at the true directions, especially at two most closely spaced correlated sources whose bias is clearly seen from the insert. Two most closely spaced correlated sources can be resolved by SPICE (note that peaks in the spatial spectrum are much larger than 2, but they are cut off at 2 to use the same scale as the other figures in Figure 1), but SPICE can yield false peaks and slight bias in the vicinity of the correlated sources and uncorrelated sources, respectively. The proposed method AULMC yields a nearly ideal spatial spectrum and provides a precise estimation for all the sources. In summary, AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum in the scenario with low SNR, small number of snapshots, and closely spaced correlated sources.
(a)
(b)
(c)
We analyze the RMSE of three algorithms under different conditions in the second simulation. The source model is the same as the first simulation. Figure 2 shows the RMSE as a function of SNR of all the algorithms and CRB in 50 Monte Carlo runs for the fixed number of snapshots 50, whereas the RMSE versus number of snapshots is shown in Figure 3 for the fixed SNR 5 dB in 50 Monte Carlo runs. Based on Figures 2 and 3, we can draw the conclusions that the RMSE of AULMC is smaller than those of other two algorithms and AULMC has the more significant performance advantages than the other two algorithms, especially in the scenarios with low SNR or small number of snapshots. One possible explanation is that AULMC can give the stable estimation in every Monte Carlo run. It can be also seen that the RMSE is close to the CRB with the increase of SNR and the number of snapshots.
In Figure 4, we display the relation between the RMSE and angel separation of correlated sources which can illustrate the resolving ability. Let two correlated sources at angles and , where the step of is , be impinged on the ULA. The SNR is 13 dB and the number of snapshots is 100. It can be seen from Figure 4 that when angle separation is , AULMC fails; however, AULMC can still provide a precise estimation as long as the angle separation is no less than and has higher resolution than the other two algorithms.
Finally, the computation time of different algorithms versus number of snapshots is shown in Table 1 for comparing the efficiency of these algorithms and SDP solver. Two correlated sources impinge on the ULA at and . The SNR is fixed at 13 dB. The computation time is obtained by the MATLAB 7.8 (R2009a) on a 2.8 GHz 4 GB PC. For AULMC, the computation time is mainly spent on the iterations of augmented Lagrange.
It can be seen from Table 1 that the computation time of SDP solver is the longest, and although the computation time of AULMC is longer than that of other two algorithms, it is comparable. Moreover, it is worth noting that the performance of AULMC is much better than that of CS-MUSIC or SPICE.
7. Conclusion
A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is proposed in CS. It is proved that the problem of minimizing the modified covariance matching criterion is an SDP, which can be transformed into the constrained quadratic programming problem solved by the augmented Lagrange method. A detailed derivation for the CRB and a theoretical performance guarantee for identifying the support are provided. Simulation results show that AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum and has more precise estimation as well as higher resolution, especially in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources.
Appendices
A. Proof of Proposition 1
Due to , we have With (34) and (A.1), we get
Therefore, it is implied by (A.2) that (35) holds and the proof of Proposition 1 is completed.
B. Proof of Proposition 2
The support can be identified if
By the Cauchy-Schwarz inequality of the trace, we have
Then, for all we have
Thus, an upper bound on the left-band side of (B.1) is given by
Similarly, we can obtain for all and hence describe a lower bound on the right-hand side of (B.1) as
Combining (B.4) with (B.6), the support can be identified if satisfies for
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported in part by the National Science Foundation of China under Grant 61201410 and in part by Fundamental Research Funds for the Central Universities under Grand no. HEUCF130804.