Two-dimensional Underdetermined DOA Estimation of Quasi-stationary Signals via Sparse Bayesian Learning

In order to improve the direction-of-arrival (DOA) estimation performance of quasi-stationary signals (QSS) using a uniform circular array (UCA), this paper addresses novel method in the context of sparse representation framework. Based on the Khatri-Rao transform, UCA can achieve a higher number of degrees of freedom to resolve more signals than the number of sensors. Then, by exploiting the two-dimensional (2-D) joint grid of UCA, the estimations of elevation and azimuth angles can be obtained from the sparse representation perspective. Finally, an expectation-maximization iteration method is developed to estimate DOAs of QSS from a Bayesian perspective. Since SBL makes full use of the sparse structure of QSS, thus the proposed algorithm possesses higher angular resolution and better DOA estimation precision compared with existing methods. Numerical simulations demonstrate the validity of the proposed method.


Introduction
Direction of arrival (DOA) estimation is an important problem in array signal processing, which is widely used in radar, sonar, wireless communication and seismic sensing. Classic subspace-based algorithms, which include multiple signal classification (MUSIC) [1] and estimation of signal parameters via rotational invariance techniques (ESPRIT) [2], have been verified as efficient estimation techniques. However, previous studies mainly focused on the Gaussian signals. In this paper, we address the DOA estimation problem in which the signals are assumed to be quasi-stationary. Quasi-stationary signals (QSS) are a class of nonstationary signals in which the statistics are locally static over a short period of time, but exhibit differences from one local time frame to another. Speech and audio signals are often recognized as QSS [3]. DOA estimation of such signals plays an important role, for example, in microphone array processing of speech signals [4] and birds monitoring systems of the airport [5]. However, it is a big challenge in some scenarios where the number of signals is more than that of sensors, which turns out to be the so-called underdetermined DOA estimation problem.
As is well known, an array with M sensors only can resolve up to M -1 QSS. In order to achieve underdetermined DOA estimation, a Khatri-Rao (KR) subspace method is recently proposed in [6] to tackle the underdetermined DOA estimation of QSS. In particular, by vectorizing the covariance matrix of array output vector, the physically underdetermined DOA estimation problem can be transformed as virtually overdetermined case. The idea behind KR transform lies in achieving the higher number of degrees of freedom (DOFs) by exploiting the difference co-array, whose virtual sensor positions are determined by the lag differences between physical sensors [7], [8].
However, previous works seldom address DOA estimation from the two-dimensional (2-D) perspective. To the best of our knowledge, [9] proposed a 2-D DOA estimator of QSS with the configuration of the L-shape array. Nevertheless, the estimation performance of L-shape array is vulnerable to angle pairing error, which may lead to DOA estimation failure. Since uniform circular array (UCA) is capable of providing 360° azimuthal coverage and identifying both azimuth and elevation angles simultaneously, it is widely employed in the 2-D DOA estimation. Though [8] has proposed DOA estimation of QSS based on the UCA, it assumes that each signal is located at a fixed and known elevation angle. Hence, it is the one-dimensional (1-D) DOA estimation in essence. In addition, the truncated error inherent in the involved manifold separation technique will degrade the estimation performance.
Although subspace-based methods proposed in [6] can be directly utilized to estimate the 2-D DOAs of QSS, its estimation performance may deteriorate significantly in low SNR or small snapshots. In order to circumvent this issue, the emerging sparse representation (SR) methods [10][11][12], which exhibit superiority in estimation precision, robustness to noise and correlation of signals, are tailored to determine the DOAs of QSS. The idea of utilizing SR, which is intrinsically different from subspace-based methods, provides a new perspective for DOA estimation. Sparse Bayesian learning (SBL) is a kind of efficient methods [13][14][15] for the sparse signal recovery in SR, which uses the expectation-maximization (EM) iteration method to estimate DOAs of QSS from a Bayesian perspective [16]. The work in [17] has demonstrated that SBL-based methods can achieve better estimation performance over conventional regularized optimization methods. In addition, the SBL-based methods do not need estimate parameters in performing the algorithms. Therefore, this paper mainly studies the 2-D DOA estimation of QSS via SBL.
The remainder of this paper is organized as follows.
Based on the data model of QSS for UCA, a virtual array with larger aperture is derived by vectorizing the covariance matrix of UCA output vector in Sec. 2. Then, Section 3 proposes novel 2-D underdetermined DOA estimation method based on the SBL. The simulations are carried in Sec. 4. Section 5 concludes the paper.
Throughout this paper, we use boldface lowercase and capital letters to denote vectors and matrices, respectively.
The operators () T , () * , () H ,  and  represent the transpose, conjugate, conjugate transpose, Khatri-Rao product and Kronecker product, respectively. The symbol E() and vec() stand for the mathematical expectation and vectorization operators, respectively. In addition, I M denotes the M  M identity matrix and diag() is a diagonal matrix composed of the elements of a column vector.

UCA Data Model from Sparse
Representation Perspective

DOA Estimation Model of UCA
In this paper, the topological structure of UCA is shown in Fig. 1. We consider a UCA with M sensors and N far-field narrowband QSS impinge on the UCA. The observation vector of the k th frame is modeled as where n(t)   M  1 is the zero-mean white Gaussian vector with covariance  2 k I M . s k (t) = [s 1 (t), s 2 (t), …, s N (t)] T and s n (t) is assumed to be a quasi-stationary process with K non-overlapped frames and the length of each frame is L, i.e., E{s n (t) 2 } = p 2 nk for  , which means the second-order statistics of QSS are static within one frame, but exhibit differences from one local time frame to another. In addition, 2  The corresponding exact local covariance in the k th frame can be written as

Khatri-Rao Transform
In this subsection, by applying KR transform to the covariance matrix k R , we can obtain The vectorized y k behaves like a new signal model.
is the new signal vector.
where the subscripts p, q for p, q =1,2,…,M in (5) are the corresponding element position in a * ( n ,  n ) and a( n ,  n ), respectively. The b pq can be obtained through the p th element in a * ( n ,  n ) multiplying by the q th element in a( n ,  n ). Consequently, we can give a general expression of b pq ( n ,  n ) ( , ) exp (cos( ) cos( )) cos( ) According to (6), we know that virtual sensors located on positions having different radius from the origin, which implies that we synthesize the virtual elements onto a non-uniform concentric circular array by using KR transform with a UCA. In addition, the number of virtual elements is far more than the physical sensors and is less than the M 2 distinct elements, this is because that there exist redundant elements in b pq . It is illustrated in [18] that the virtual elements of UCA with odd number of sensors is M(M -1) + 1 while for the UCA with even number of sensors is M 2 /2 + 1. Therefore, the DOFs of UCA is greatly increased based on the KR transform, and we can apply (4) instead of (1) to achieve underdetermined DOA estimation.
denotes the noise.

2-D Sparse Representation
Following the convention in the context of SR framework, a fixed sampling grid is selected firstly that serves as the set of all candidates of DOA estimates. In this paper, in order to achieve 2-D DOA estimation of QSS, the azimuth and elevation angles are equally sampled into discretized grid sets of  = {θ̅ 1 , θ̅ 2 ,…, θ̅ H θ }(0°, 360°) and  = {̅ 1 , ̅ 2 ,…, ̅ H  }(0°, 90°), respectively, as is shown in Figs. 2(a) and 2(b), where H θ >> K, H  >> K. In addition, it is assumed that the true DOAs are exactly on the sampling grid sets  and . In order to facilitate the writing in subsequent section, we combine θ̅ h θ (1  h θ  H θ ) and Thus, a joint 2-D sampling grid is constructed and the number of total grid points is H = H θ  H  , as is shown in Figs. 2(c), which satisfies the reconstruction condition H >> K. In Fig. 2 and there exists correspondence relationship  1 = (̅ θ 1 ,̅ 1 ), The advantage of representing the 2-D discretized grid as a single vector is that we can handle the 2-D problem in a 1-D angular space. Once  h is solved, we can obtain the corresponding DOA estimate (θ̅ ,̅ ) according to the correspondence relationship in Fig. 2(c).
As a result, equation (7) can be expressed as the following sparse representation model is an over-complete dictionary and it is defined as a collection of the steering vector over the entire grid set Ψ . In addition, it should be noted that 1 2 =[q ,q , q ] has K columns, so equation (9) satisfies the multiple measurement vector (MMV) model.
In order to reduce the computation complexity of the signal reconstruction process and the sensitivity to the measurement noise, we can apply the singular value decomposition (SVD) to (9) and Y can be decomposed as where the columns of U and V are the left-singular and right-singular vectors, respectively, while S   M 2 K is a diagonal matrix and can be expressed as , we can obtain . In addition, Y  contains most signal information of the matrix Y . Based on the SVD, the column dimension of Y  is much smaller than Y , so the computation complexity is reduced.

2-D DOA Estimation Based on SBL
In this section, we use the SBL method to solve equation (12). The reason why we choose the SBL is that the global minimum points of SBL correspond to the most sparse solution and the local minimum points of SBL are very few.

Noise and Sparse Signal Model
For a complex Gaussian distributed random variable ( , ) u CN    with mean  and covariance , the probability density function (PDF) can be expressed as In the SBL, it is usually assumed that the noise satisfies complex Gaussian distribution, thus we can obtain where  denotes the noise precision.
Further, for ease of inference, we assume that  satisfies the Gamma distribution since it is a conjugate prior of the Gaussian distribution.
where a and b are scale parameters.
In addition, the likelihood function of (12) can be written as For the sparse signal matrix Q  of interest, we adopt the two-stage hierarchical sparse prior model to describe Q  , which guarantees that most rows of Q  being zeros.
According to the above analysis, the joint PDF can be expressed as

Bayesian Inference
Since p(Q̃,,Ỹ) cannot be explicitly calculated, we are able to use the EM iteration algorithm to perform the Bayesian inference. By treating Q̃ as a hidden variable, whose posterior distribution is where n μ stands for mean and Σ denotes covariance.
can be rewritten as In (21), p(Ỹ, ) is independent of Q̃ and can be calculated according to  and . Therefore, we have By combining (14)-(23), we can obtain In order to calculate  n and , we need estimate the hyperparameter  and . Based on a maximum a posteriori (MAP) optimal estimate, they can be estimated by denotes the expectation operator with respect to p(Q̃Ỹ, , ) and p(Ỹ, Q̃, , ) has been given in (19). Therefore,  and  satisfy     ( | , , ) , In , , , ( 1)

Simulation and Results
In order to verify the performance of the proposed algorithm, the following simulations are carried. Unless otherwise specified, we set general simulation parameters. The UCA has M = 5 sensors and the radius of the UCA is r = /2. The incident signals are regarded as QSS. Each QSS has K = 30 frames and the length of each frame is L = 1024. In the proposed method, we select scale parameters a = 10 -4 , b = 10 -4 , c = 1, d = 10 -2 , and initialize  

2-D Spatial Spectra Distribution of the Proposed Method
Firstly, we simulate the spatial spectra distribution of both overdetermined and  can correctly estimate DOAs even though the number of signals is larger than the number of sensors, which means the proposed method has the ability to achieve 2-D underdetermined DOA estimation.

DOA Estimation Precision of Different Methods
In this subsection, we compare the DOA estimation precision between the proposed method and state-of-the-art schemes, such as KR-MUSIC [6], KR-CAPON [6] and FO-MUSIC [19] through simulation experiment. The rootmean-square error (RMSE) is introduced as the evaluation standard [20], where P denotes the number of Monte Carlo trials and N is the number of signals. ˆp n  is the estimated value of  n in the p th trial, Meanwhile, the Cramer-Rao Lower Bound (CRB) is also plotted as a benchmark in following simulations [21]. In order to facilitate the precision analysis, we fix elevation angle  at 90°, which means the elevation angle could be ignored in following simulations. The azimuth angle setting is the same as in Sec. 4.1. Unless otherwise specified, we use the discrete grid in the azimuth angle range of   (0°, 360°) with a step size of 1°. By performing 500 times Monte Carlo trials, the overdetermined DOA estimation RMSE versus SNR and snapshots are plotted in Fig. 4, respectively. We record 1024 snapshots and Fig. 4(a) plots the RMSE of various algorithms varying with the SNR. When SNR is 5 dB, the RMSE versus snapshots of each frame is shown in Fig. 4(b). From Fig. 4, we can see that the RMSE of four methods decrease rapidly with the increase of SNR or snapshots. In general, the FO-MUSIC method has the worst estimation accuracy, because there exists error in the estimated four-order cumulants due to the finite sampling snapshots, which deteriorates the estimation performance. By using the Khatri-Rao transform, the virtual array aperture is extended. Therefore, the KR-MUSIC and KR-CAPON have improved accuracy in terms of RMSE compared with FO-MUSIC, and the estimation performance of KR-MUSIC is a little better than KR-CAPON. The proposed method has the highest DOA estimation precision over all the range of SNR or snapshots and the trend of its performance curve is the same as the CRB when the SNR is larger than 5 dB and/or when the snapshots is larger than 1024. This is because the proposed method makes full use of the sparse structure of QSS from a Bayesian perspective and the proposed method does not require the parameter estimation in the calculation process by using the Bayesian inference.
Similarly, we can also obtain the underdetermined DOA estimation RMSE versus SNR and snapshots, as is shown in Fig. 5. It is seen from Fig. 5 that the proposed method still has the best DOA estimation precision compared    with existing three methods and the trend of its performance curve can still reach the CRB at high SNR and/or big snapshots.

Angular Resolution Comparison
In order to verify the angular resolution of the proposed method, we consider two closely spaced QSS from the azimuth angles  1 = 100°-Δ and  2 = 100°+ Δ, respectively, where Δ is varied from 0° to 8° with a 0.5° step.
We define that two signals can be correctly resolved when there are two peaks in the spatial spectrum and the estimated DOAs satisfy 1 1       and 2 2       .
The SNR is 5 dB and the snapshots is 1024. By conducting 200 times Monte Carlo trials, the resolution probability versus Δ is plotted in Fig. 6. It can be seen from Fig. 6 that the proposed method has the best angular resolution and the resolution probability can reach 1 when Δ = 2°. In addition, the angular resolution of KR-MUSIC is a little better than FO-MUSIC.

Conclusion
This paper studies the 2-D DOA estimation of QSS in the context of SR framework. Firstly, the Khatri-Rao transform is applied to the UCA data model, which makes that the virtual array aperture of UCA is extended, so that the proposed method has the ability to estimate more signals than number of sensors. Then, by 2-D joint grid discretization of UCA, the azimuth angle and elevation angle of QSS can be estimated simultaneously from SR perspective. Finally, an expectation-maximization iteration method is developed to estimate DOAs of QSS based on the SBL method. Since this paper makes full use of the sparse structure of QSS from a Bayesian perspective, thus the proposed method has better estimation precision and angular resolution compared with existing methods. In addition, the proposed method does not require the parameter estimation in the calculation process, which facilitates the engineering application. However, the assumption that the true DOAs are located on the predefined spatial grids is not always valid in practical implementations. The DOA estimation performance of the proposed method in this paper is limited by the off-grid effect of signals and mismatch of the SR model. Therefore, future research efforts will aim to solve the problem of off-grid and model mismatch.