Two-Dimensional Underdetermined DOA Estimation of Quasi-Stationary Signals using Parallel Nested Array

In this paper, a two dimensional underdetermined direction of arrival estimation (DOA) of quasi-stationary signals using a parallel nested array structure is investigated. The quasi-stationary signals have the statistical property that they remain locally static over one frame but exhibit differences from one time frame to others. The special time domain property enables us to perform underdetermined direction-of-arrival estimation in time domain. By exploiting the temporary diversity of the quasi-stationary signals and effective difference coarray virtual array aperture provided inherently in the parallel nested array, more degrees of freedom can be used to resolve DOA estimation. The Khatri-Rao operation for the cross covariance matrix of the subarrays received data is adopted to convert the 2-D DOA estimation problem into two separate onedimensional DOA estimation problems. Then, a subspacebased estimation of signal parameters via rotational invariance technique and a sparsity-based sparse Bayesian learning are proposed to realize the according one-dimensional DOA estimation. And the estimated azimuth and elevation angles can be properly automatically paired. Simulation results are carried out to demonstrate the effectiveness of the proposed algorithms for the 2-D underdetermined DOA estimation.


Introduction
The direction-of-arrival (DOA) estimation problem plays an important role in array signal processing and has wide applications in radar, sonar, mobile communications [1][2][3], etc. The quasi-stationary signals (QSS) are one of most of important signals frequently encountered in the microphone array speech processing and electroencephalogram. And, in an airport, the DOA estimation of the quasistationary characteristic can also be used to prevent colli-sions between the aircrafts and birds. Thus, due to these application backgrounds, it is very meaningful to acquire the accurate direction parameters of the QSS. The onedimensional DOA estimations of the QSS are detailed analyzed in [4][5][6][7][8]. But the researches about the two-dimensional DOA estimation of the QSS are few. In [9], [10], tensor modeling and Khatri-Rao subspace based method are proposed to deal with the 2D DOA estimation problem by using the L-shaped array. However, these proposed methods do not give the detailed analysis after the Khatri-Rao operation about the array steering matrix and the difference coarray domain of the L-shaped array cannot be obtained. A number of two-dimensional (elevation angle and azimuth angle) high-resolution algorithms using uniform array geometry structures (such as the rectangle array, the circular array, L-shaped array and parallel linear array) have been proposed to realize the 2D DOA estimation in [11][12][13][14]. But traditionally uniform array geometry structures suffer from limited degree of freedom and some algorithms need the extra angle pair matching process. For example, the maximum number of the estimated signals using propagator method [15] is limited to half of the sensors.
In order to acquire accurate DOA estimation when the number of sources exceeds the number of sensors, exploiting non-uniform sparse arrays to achieve a larger array aperture without increasing the number of sensors has been a preferred choice in recent years [16][17][18][19][20][21][22][23]. By utilizing the second-order statistics of the received signal, non-uniform arrays can construct virtual difference co-array with larger aperture to improve the DOFs. Among these, the introduction and analysis of the nested array and the coprime array have attracted significant attention due to their ability to estimate the direction-of-arrival (DOA) of more sources than the number of physical sensors. A nested array comprising two uniformly sampling subarrays has many advantages compared to other popular non-uniform arrays, such as minimum redundancy arrays (MRA) [23] , minimum hole arrays (MHA) [24] and coprime arrays [17][18], [21][22]. For a given number of physical sensors, MRAs and MHAs need use exhaustive search to find the optimal sensors position. Besides, finding the suitable covariance matrix corresponding to a large array requires a rather complicated time-consuming iterative process. And the difference coarray domain of coprime array contains some missing elements, thus, cannot exploit all of the DOFs for DOA estimation. Unlike the coprime array, the DOFs and the number of unique and consecutive difference coarray sensors of nested array can be obtained with the exact expression under the condition of fixed sensors.
The recent proposed sparse arrays have a very brilliant application prospect due to their degrees of freedom enhancement capability. However, these array aperture expansion capabilities are mainly limited to the 1-D case. In order to realize the 2D DOA estimation when the number of impinging sources is more than the number of physical sensors, the parallel coprime array structure is investigated in [25][26][27]. But similar to the analysis above, the coprime array cannot eliminate the effect of the "holes" in the difference coarray domain both for 1-D DOA estimation and 2-D DOA estimation. And the above mentioned methods based on parallel coprime array structure cannot deal with the 2-D DOA estimation of the quasi-stationary signals as the cross-covariance matrix of subarray signals behaves like a multiple measurement vector model. Thus, in this paper, an effective 2-D DOA estimation method based on parallel nested array structure is proposed and the proposed method can make full use of the freedoms of the difference coarray in the cross-covariance matrix. Then, a multiple measurement vector model can be constructed by using the time diversity of the quasi-stationary signals. Finally, a subspace-based estimation of signal parameters via rotational invariance technique (ESPRIT) algorithm and a sparsity-based sparse Bayesian learning (SBL) algorithm are derived to realize the 2-D DOAs estimation. The proposed method can convert the 2-D DOA estimation problem into two separate one-dimensional DOA estimation problems, which is quite attractive for decoupled 2-D DOA estimation. Both the proposed ESPRIT algorithm and SBL algorithm can realize the 2-D underdetermined DOA estimation for the quasi-stationary signals. And we also compare the estimated performance for the ESPRIT algorithm and SBL algorithm under different situations. Simulation results verified the effectiveness and superiority of the proposed methods in terms of detectable ability, estimation accuracy, and resolution ability.
The rest of this paper is organized as follows. Section 2 introduces the system model of the parallel nested arrays. In Sec. 3, the subspace-based ESPRIT algorithm and sparsity-based SBL algorithm are proposed to deal with the multiple measurement vector model of the quasi-stationary signals. Simulation results are provided in Sec. 4 and Section 5 concludes this paper.
Notations: Throughout the paper, we use the lowercase (uppercase) boldface symbols to represent vectors (matrices). (·) * ,(·) T , and (·) H denote conjugate, transpose, and conjugate transpose, respectively. Re(·) denotes the real part. The symbols  ,  and  denote Khatri-Rao product, Kronecker product, and Hadamard product, respectively. E(·), vec(·) and angle(·) denote statistical expectation, vectorization, and phase angle operator, respectively. I M and Π M are M  M identity matrix and reverse identity matrix (90° rotation of I M ), respectively. 0 M  M is a M  M zero matrix. The notation diag(a) stands for constructing a diagonal matrix with the elements a on its diagonal.

The Model of Parallel Nested Array
The planform of the parallel nested array is illuminated in Fig. 1, the parallel nested array consisting of two parallel-arranged sparse uniform linear arrays in the xy plane, with a half-wavelength distance between them. The first subarray along the y-axis is a dense uniform linear array with inter-element spacing λ/2. The second subarray with M elements is a sparse uniform linear array with interelement spacing λM/2. Therefore, the positions of the parallel nested array can be denoted as Assume that K narrowband far-field uncorrelated quasi-stationary signals (QSS) impinge on the parallel nested array from ( k ,  k ), k = 1,2,…,K, where  k and  k denote the elevation angle and azimuth angle . To simplify the representation, we redefine ( k ,  k ) to denote the 2D angle, where  k denotes the angle between the incident direction and y-axis, and  k is the angle between the incident direction and x-axis, for all k = 1,2,…,K. Thus, the relationship between ( k ,  k ) and ( k ,  k ) can be denoted as cos( ) sin( ) cos( ).
Then, the received data vectors corresponding to the two subarrays can be expressed as Physical Element in Subarray 1 Physical Element in Subarray 2 -axis x -axis are the steering vectors corresponding to angles ( k ,  k ) for k = 1,2,…,K. A 1 = [a 1 ( 1 ), a 1 ( 2 ),…,a 1 ( K )] and A 2 = [a 2 ( 1 ), a 2 ( 2 ),…,a 2 ( K )] are the manifold steering matrices corresponding to the first subarray and second subarray, respectively.  = () = diag(exp(-jcos 1 ), exp(-jcos 2 ), ,…, exp(-jcos K )), denotes a diagonal matrix corresponding to angle  k , k = 1,2,…,K. The noise n 1 (t) and n 2 (t) are the additive white Gaussian noise with mean zeros and covariance matrix  1 2 I M and  2 2 I M , respectively. And the noises are assumed to be independent to each other and uncorrelated with the signals.
The wide-sense stationary characteristic of QSS within a frame can be expressed as where L denotes the frame length and p denotes the frame index. Equation (8) means that the second-order statistic is static within a frame but time-varying in different frames. Before acquiring the cross-covariance matrix between the first subarray and the second subarray, we firstly use a permutation matrix Π M to (4) By using this permutation operation, the virtual steering matrix corresponding to vectorization of the cross covariance matrix of the two subarrays has no overlapped elements and it can achieve M 2 degree of freedoms. Then, the cross-covariance matrix between x 1z (t) and x 2 (t) within the p-th frame can be expressed as 12 H where is the diagonal signal covariance matrix within the p-th frame. By adopting the cross covariance matrix in (10), the effect of noise can be eliminated due to temporal uncorrelated characteristic.
In practice, the exactly local cross-covariance matrix R x 12 in (10) cannot be obtained due to the finite number of snapshots. Thus, the local sample cross-covariance matrix is estimated using the L available snapshots in each frame, which can be expressed as 12 1 By vectorizing the matrix R x 12 , we can obtain the following M 2  1 measurement vector 12 H H 1 2 The vector y p can be regarded as a new received data vector and b p can be regarded as the virtual impinging sources.
is the new virtual direction matrix whose elements can be denoted as j (  1) Thus, the total number of the degrees of freedom is M 2 , enabling to estimate more impinging sources than the number of sensors. The QSS has the statistical property that it remains locally static over one frame but exhibits difference from one frame to others. By stacking the cross covariance matrices of all P frames and the corresponding vectorization forms can be expressed as 1 2 [ , ,..., ] P   Y y y y AB (15) where B = [b 1 , b 2 ,…,b P ] denotes the new signal matrix, which is formed by the second-order statistics corresponding to different frames.

The Proposed Methods
In this section, we first introduce the subspace based-ESPRIT algorithm and then demonstrate the sparsity-based SBL algorithm for the 2-D underdetermined DOA estimation of the quasi-stationary signals

The Real-valued ESPRIT Algorithm
In this subsection, a real-valued ESPRIT algorithm is proposed. Compared to the complex data multiplication operation, the real data can save about a quarter of the computational complexity. The complex value received data vectors of the array can be transformed to real data by using a unitary matrix. Define the unitary matrix as follows j 1 if is even, Due to the multiple measurement vector model in (15), we can efficiently construct a new virtual covariance matrix as follows H / P  YY R YY (17) where the frame number P behaviors like the sampling snapshot. However, as the number of frames is limited in practice, the obtained covariance matrix in (17) is Hermitian but generally not persymmetric. Thus, the persymmetric Hermitian covariance matrix can be obtained as follows where Q H Π M 2 = Q T , Π M 2 Q= Q * are used from the third row to the fourth row. The eigenvalue decomposition of the real-valued covariance matrix R r can be denoted as s s s n n n  r

R E D E +E D E
where E s is the signal subspace matrix that corresponds to the K largest eigenvalues. Similarly, E n is the noise subspace matrix corresponding to the remaining eigenvectors. D s and D n denote the eigenvalue matrix whose diagonal entries are the eigenvalues in the descending order.
According to (18) and (19), the new direction vector matrix corresponding to R r can be denoted as H  Ω Q A (21) and the relationship between the signal subspace E s and Ω can be denoted as s  E ΩT (22) where T is a nonsingular matrix.
From (15), we can find the virtual direction matrix A is a vandermonde matrix, thus the steering vector a( k ) satisfies are the selection matrices. Then according to [28], the relationship in (23) can be transformed to the real-valued form 1 2 and 1 diag(tan( cos / 2),..., tan( cos / 2)).
From (22), (24), (25), (26) and (27) Thus, the eigenvalues of Ψ will be the diagonal elements of Φ . Denote the eigenvalue of Ψ as u k , k = 1,2,…,K. Then  k can be estimated via If k, k = 1,2,…,K are obtained, the K P matrix B can be estimated by using the least squares fitting method based on the linear model in (15) Finally,  k can be obtained as follows ˆa rccos(angle( ( ,1)) ), 1, 2,..., where B(k,1) denotes the k-th element of the first column in B. In order to improve the performance, all columns of the matrix B can be used to estimate  k and average all the P solutions.

The Sparse-based Sparse Bayesian Learning Algorithm
In this subsection, a sparse-based sparse Bayesian learning algorithm is proposed to improve the 2-D DOA estimation performance for the quasi-stationary signals. According to (15), the multiple measurement vector model can rewritten as where ε reflects the effect of the noise due to the finite number of frames. The signal matrix Ŷ can be sparsely represented over the entire discretized angular grids as where Â is overcomplete dictionary with steering vectors ã( g ) over all possible grids  g , g = 1,2,…,G with G >> K.
It is important to note that true angle parameter  k , k = 1,2,…,K is indicated by the positions of the nonzero entries in matrix X which is jointly row sparse, i.e. all columns of X are sparse and share the same support.
According to [29][30][31], a white complex Gaussian prior is placed on the error matrix E where  0 denotes the noise precision and a Gramma hyperprior which is a conjugate prior of corresponding to the Gaussian distribution is placed on  0 .  [29], [30]. The prior distribution of Ŷ can be expressed as Assume that the entries in jointly sparse matrix X are drawn from the product of the following zero-mean complex Gaussian distributions: A zero-mean complex Gaussian prior and a Gamma prior are placed on the jointly sparse matrix X. The two-stage hierarchical prior can be denoted as: and 0 It can be seen that all columns of X are independent and share the same distribution because of the group sparse property. According to [30], the real and image parts of each column satisfy Laplace distribution and exhibit the same strong peak at the origin. Thus, the two-stage hierarchical prior is suitable to guarantee the sparsity of X.
since p(Ŷ) is independent of the hyperparameters. An expectation-maximization (EM) algorithm is implemented that treats X as a hidden variable and turns to maximizing E{p(X, Ŷ,  0 , )}, where p(X, Ŷ,  0 , ) is given in (42) and E{}denotes an expectation.
. Following a similar procedure as in [31], it is easy to obtain the following updates of  0 and  Once  k are estimated, the corresponding  k can also be obtained according to (32) and (33) for k = 1,2,…,K.

The Complexity Analysis
The computational complexity of the proposed subspace-based ESPRIT algorithm mainly contains the construction of the covariance matrix, the singular value decomposition and least square. Thus, the computational burden for the proposed ESPRIT algorithm is o(M 2 LP + M 6 + 2K 3 + 2K 2 M 2 + KM 2 P).
The proposed sparse Bayesian learning algorithm mainly focuses on the construction of the covariance matrix, the iteration operations and least square. Thus, the total computational complexity is o(M 2 LP + M 2 G 2 + K 3 + 2K 2 M 2 + KM 2 P), where  denotes the number of iterations.

Simulation Results
In the first experiment, we demonstrate the underdetermined spatial spectrum estimation performance of the proposed methods. We consider a parallel nested array with M = 4, the planform of the nested parallel array is demonstrated in Fig. 1. We assume K = 10 narrowband uncorrelated impinging sources that are uniformly distributed between [40°, 45°] and [135°, 140°]. The frame length L is taken to 400 and a total of P = 30 frames are used. The signal-to-noise ratio (SNR) is set to 10 dB and the predefined spatial sampling grids for the sparse Bayesian learning (SBL) algorithm are from 0° to 180° with the sampling interval being 1°. The underdetermined spatial spectrum estimations for the proposed methods are shown in Fig. 2. From Fig. 2, we can find that both the proposed ESPRIT and SBL algorithm can resolve the DOAs when the number of impinging sources is larger than the sensors. But, from Fig. 2(b) and Fig. 2(d), we can find the performance of the ESPRIT algorithm is poorer compared to the SBL algorithm. And from Fig. 2(a) and Fig. 2(c) and Fig. 2(b) and Fig. 2 (d), we can find the performance of estimated  k is better than the estimated  k for k = 1,2,…,K.
In the second experiment, the underdetermined RMSE performance as a function of SNR for the proposed ESPRIT algorithm and the SBL algorithm is illuminated in Fig. 3. The RMSE is defined as where (̂ ik , β^i k ) is the estimated DOA of ( k , β k ) for the i-th Monte Carlo trial, i = 1,2,…,I. And the RMSE is obtained by using I = 200 independent trials in all simulations. We consider K = 10 narrowband uncorrelated impinging sources that are uniformly distributed between [40°, 45°] and [135°, 140°]. In Fig. 3, the SNR increases from -6 to 10 with a step of 2 and the frame length L and frame number P are fixed at 400 and 30, respectively.
The results shown in Fig. 3 indicate that the RMSE of the sparsity based SBL algorithm has the better performance than the ESPRIT algorithm. And the RMSE of  is smaller compared to the RMSE of  for both SBL algorithm and ESPRIT algorithm.
In the third experiment, we consider the RMSE performance versus different frame length and frame number for the SBL algorithm and ESPRTI algorithm. In Fig. 4(a), the frame number increases from 10 to 50 with a step of 4 and the frame length is fixed at 400 for all different frame numbers. The results shown in Fig. 4(a) indicate that the performance of the proposed ESPRIT algorithm is sensitive to the variation of the frame number. In Fig. 4(b), the number of frame length increases from 100 to 1000 with a step of 100 and the frame number is fixed at 30 for all different frame lengths. The results shown in Fig. 4(b) indicate that the performance of the RMSE decreases when the frame length increases. There exists a relatively large error for ESPRIT algorithm when frame number P < 18 due to the error existed in covariance calculated in (16), and the RMSE performance of  using least squares fitting algorithm still has relatively poor performance compared to the RMSE of  using the proposed methods. In the last experiment, we compare the angular superresolution capability of the different algorithms mentioned above, the frame number P is fixed at 30 and the frame length is 400. The SNR is fixed at 10 dB and the 2-D DOA ( 1 ,β 1 ) of the first impinging source is (30°,40°). Two close sources are considered, the second 2-D DOA ( 2 ,β 2 ) is ( 1 + , β 1 + β) where both  and β varies from 1° to 10° with a step of 1°. We define that  and β are correctly resolved when they satisfy  1 -̂ 1   1.5°,  2 -̂ 2   1.5°, and β 1 -β^1  1.5°, β 2 -β^2  1.5°. Two hundred Monte-Carlo trials are conducted for each  and β to acquire the resolution probability versus  illuminated in Fig. 5. We find that when   2°, both proposed methods can achieve a high angular resolution for  while the spacing angle are needed more than 4° for β to achieve a better angular resolution.

Conclusion
In this paper, a parallel nested array is designed to deal with the underdetermined 2D DOA estimation of the quasi-stationary signals. By exploiting the parallel nested array configuration, a large number of degrees of freedom can be obtained from the difference co-array domain with the aid of the cross-covariance matrix. And, two efficient algorithms are proposed to convert the 2D DOA estimation into two 1D problem and the estimated elevation and azimuth angles are automatically paired in the underdetermined case. Both the proposed ESPRIT algorithm and SBL algorithm can exhibit excellent estimation performance in term of achievable DOFs, estimation accuracy and spatial resolution. Besides, the proposed sparsity-based Bayesian learning method has better angle estimation performance and higher spatial resolution compared with the proposed subspace based method.
include inverse synthetic aperture radar imaging and electromagnetic environment effects.
Shunping XIAO was born in Jiangxi, China, in 1964. He received the B.S. and Ph.D. degrees in Electronic Engineer-ing from the National University of Defense Technology (NUDT), Changsha, China, in 1986 and 1995, respectively. He is currently a Professor with the NUDT. His research interests include radar target recognition and radar signal processing. Prof. Xiao is also a senior member of CIE.