2-D DOA Estimation via Matrix Partition and Stacking Technique

A novel approach is proposed for the efficient estimation of the two-dimensional (2-D) direction-of-arrival (DOA) of signals impinging on two orthogonal uniform linear arrays (ULAs). By partitioning the cross-correlation matrix (CCM) between two ULAs data into a great deal of submatrices and making use of the submatrices and the symmetric subarrays, an extended correlation matrix is constructed, and then uses the modified ESPRIT approach to extract out the so-called Kronecker Steering Vectors (KSVs)of which each is the Kronecker product of the elevation and azimuth angle with a one-to-one relationship. Upon that the proposed method yields the estimate of the 2-D DOA efficiently without requiring the additionally computational burden to remove the pair-matching problem. Furthermore, the main idea of the matrix partition and stacking is to much-enhanced subspace estimate. So based on the use of the concept, the proposed method’s performance is better than the existing similar approaches. Meanwhile, unlike the traditional subspace methods, it is shown that the proposed can resolve the same uncorrelated sources as the number of subarray sensor through a delicate partition-and-stacking process. Simulation results demonstrate that the proposed method is superior to the existing techniques in both DOA estimation and the detection capability of sources.


Introduction
The problem of two-dimensional (2-D) DOA (i.e., azimuth and elevation angles) has recently received increasing attention. 2-D DOA problem may be closer to some practical environment than 1-D, for instance, using an airborne or a spaceborne array to observe ground-based sources. Additionally, in the last three decades a number of the highresolution direction finding methods [1] have been studied in the context of 1-D estimation (e.g., the azimuth angle) of multiple plane waves. Among them, MUSIC (Multiple Signal Classification) [2] and ESPRIT [3] are considered the two most popular algorithms. Many 2-D DOA methods are based on the two algorithms. Specially, the latter method has two main advantages over the former. First, the ESPRIT algorithm requires less computational burden and storage space due to that it does not require to search over the whole parameter space. Second, independent of the array response, the ESPRIT algorithm is more robust to array calibration errors. Therefore, the ESPRIT algorithm and its variations [4,5], which are widely devoted to the problem of 2-D DOA estimation with planar arrays, have received considerable attention in the array-processing literature [6][7][8][9][10][11][12][13]. Majority of the planar arrays required to implement() these techniques can be divided into three types: the twoorthogonal uniform linear array (the L-shaped array) [6][7][8][9], the triangular array [10,11], and the rectangular array [12,13]. Although the L-shaped array has a simpler configuration compared to the rectangular and triangular ones, it enjoys higher accuracy among these configurations [6]. Thus the Lshaped array has received increased attention in dealing with 2-D DOA estimation problems recently.
In [9], Tayem and Kwon presented a computationally simple 2-D DOA estimation with the propagator method using one or two L-shaped arrays. They showed that it is possible to decompose the 2-D problem into two independent 1-D problems by using the L-shape array for reducing the computational burden significantly. But the two independent sets of angles would have to be properly paired together using some technique [14]. Different approaches have been put forward in the literature. For example, Kikuchi et al. [15] describe a cross-correlation technique obtaining the correct parameter pairs by constructing a Toeplitz matrix whose first column and row are the diagonal elements of the CCM and its conjugate transpose including the relation between the corresponding combinations of the elevation and azimuth angles. Unfortunately, the Kikuchi's approach still causes the pair-matching problem when the difference of the corresponding combinations of the 2-D angles {cos θ k − cos φ k } K k=1 is small and the signal-to-noise ratio (SNR) is low [8]. Furthermore, it only employs the CCM to deal with the pair-matching problem such as the pairing algorithm suggested in [9], but does not exploit its characteristics to improve the estimation performance. Recently, Gu and Wei [8] described a technique which utilized joint singular value decomposition (JSVD) for two submatrices selected from CCM unaffected by the additive noise on the condition because of the reasonable assumption for noise and signal to achieve automatic pairing and estimate 2-D DOA. In this way, the JSVD technique in [8] takes at least two advantages over the technique suggested in [15]. First, it has been shown that the JSVD technique does not need additional steps to deal with the pair matching problem. Second, the JSVD is superior in estimating the 2-D DOA, especially at low SNR and with a small number of snapshots.
Therefore, the objective of this paper is to develop another technique for 2-D DOA estimation of multiple sources from two orthogonal ULAs (L-shape array) by partition and stacking the CCM between two ULAs, showing new aspects associated to it. It is unlike the algorithm proposed in [8]. The authors divided the CCM into two submatrices. Herein we will partition the CCM into more submatrices and stack them to form an extended correlation matrix. The extended matrix is formed by taking full advantage of the structure in the underlying data models to apply the modified ESPRIT approach to efficiently obtain the KSVs containing the one-to-one information of the azimuth and elevation angles. The 2-D DOAs automatically paired are then found one-by-one by decomposing the KSVs. Actually, not only ESPRIT but also MUSIC or other methods can be modified to process the extended matrix. However, according to the ESPRIT advantages mentioned above, we modify ESPRIT to process here. Finally, through extensive simulations, it can be demonstrated that the proposed technique overcomes the problems encountered in the Kikuchi method and offers improved estimation performance as well.
The organization of this paper is as follows: We formulate the problem in Section 2. In Section 3, the 2-D DOA estimator is proposed. Experimental results are provided in Section 4 to demonstrate the performance of the proposed technique. Finally, some conclusions are given in Section 5.

Problem Formulation
The algorithm presented in this paper applies to sensor arrays with two orthogonally ULAs. For simplicity of presentation, we focus our derivation on an L-shape array of sensors. Consider two ULAs with inter-element spacing d, making up of the L-shape array configuration in the x − z plane as shown in Figure 1. Each linear array consists of M sensor  Figure 1: Array elements configuration for the joint elevation and azimuth DOA estimation [8].
elements and the element placed at the origin belongs to the z-axis. Suppose that there are K narrowband sources with wavelength λ impinging on the array from different directions, where the number of sources K is known or preestimated using some detection techniques (e.g., [16,17] and references therein). The kth source has an elevation angle θ k and an azimuth angle φ k . These sources are assumed to be in the far-field with respect to the sensor location. The observed signals at the ULA along the x-axis and the ULA along the zaxis are given, respectively, by The matrices and vectors in (1) It is also assumed that these Ksources are uncorrelated with each other and the additive noises {n z (t), n x (t)} are independent of the signal samples [8,15,18]. The elements of {n z (t), n x (t)} are white Guassian random processes with zero-mean and variance σ 2 : are the M × K array response matrices in the x subarray and z subarray, respectively. The response vectors in the elevation and azimuth directions are expressed as where α k = 2πd cos θ k /λ and β k = 2πd cos φ k /λ.
Remark A. Although the assumptions above on the additive noises are the spatial white to consistent with those of Kikuchi's, but our method is less strict to the spatial correlation of the noise due to the fact that our method makes use of the CCM between the the different ULA, where the additive noise can be eliminated. Additionally, it is necessary for the Kikuchi's method to assume the same number of sensors in the different ULA, while our method has not the limitation if only both of them are greater than or equal to the number of incident signals due to the fact that our method makes use of the many submatrices selected from the CCM. Therefore, compared to the Kikuchi's method our method is more flexible to design the array configuration.

Proposed Method
Under the assumptions of data model, we easily get a crosscorrelation matrix R xz between the signal vectors x(t) and z(t) described in [8]: where Λ = diag{ρ 1 , . . . , ρ K } and ρ k is the power of the kth source, and E{•} denotes the statistical expectation. The superscript Tdenotes the conjugate transpose. It is easy to see that the cross-correlation matrix R xz is not affected by the additive noise because of the assumptions mentioned above. Now we introduce our proposed method.

Partition and Stacking
Processing. Let us define a selection matrix of dimension p × M as follows: where 0 a×b is a zero matrix of size a × b and I p an identity matrix of dimension p. It is well known that the M×Mmatrix R xz can be divided into (M − p +1) 2 overlapping submatrices of dimension p × p. The corresponding submatrix R i, j are, hence, defined as follows: where A ab is the bth row to the (b + p − 1)th row of A a (a = x, z). According to the ULA assumption, we can get the relationship A ab = A a1 Ω i−1 a , where the matrix Ω i−1 a denotes the (i − 1)th power of a K × K diagonal matrix along the a-axis Ω a = diag(γ a,1 , . . . , γ a,K ) and Therefore, (6) can be written as Next, we can obtain a stacking matrix R of size p(M − p + 1) × p(M − p + 1): where Remark B. In order to avoid the effect of additive noise, the auto-correlation matrices of the submatrices are not used in the computation of the stacking matrix in (10), which contains all information of the elevation and azimuth angles. The methods in [7,9,15] exploit only the auto-correlation matrices to estimate the elevation and azimuth angle in each ULA along the z and x axis, respectively, and the CCM between different ULA is only employed to obtain the pairmatching information described in [7,9]. Then, rewrite (9) in explicit form by take full advantage of array structural information: where the superscript c denotes the complex conjugate and • is the Khatri-Rao matrix product defined by where ⊗ is the Kronecker Product [19], and B x = A x with the last p rows deleted,

EURASIP Journal on Advances in Signal Processing
Then, we can write the stacking matrix R as where b(θ k ) and b(φ k ) are the first (M − p) elements of a(θ k ) and a(φ k ), respectively. d i (θ k , φ k ) (i = 1, 2, k = 1, . . . , K) is defined as the kth Kronecker Steering Vector (KSV) containing the kth one-to-one elevation and azimuth angle information of the ith combination. It can be verified [18] that D 1 (θ, φ) and D 2 (θ, φ) are full column rank if min{p(M − p + 1), M + 1} > K. It is worth noting that our method can detect M signals when p(M − p + 1) ≥ M + 1.
Hence, R is of rank K.
As the array geometry mentioned earlier, A x and A z have the following relationship: where J M is the "reversal identity matrix" obtained by reversing the order of the rows of I M . Now, recall (4) and use (14) to obtain It is easy to find that Π = Ω −M x ΛΩ M z is also a diagonal matrix. Hence, we can construct another stacking matrix following the same steps as (6)-(13) Due to the similarity between R and JR, it can be similarly shown that rank of JR is K if M + 1 > K. (13) and (16), we can see that R and JR are constructed from the same KSVs and the different diagonal matrices (i.e., Λ and Π). The authors of [8] proposed a new scheme to extract the signal subspace by the joint singular value decomposition of the two concatenated matrices with the same left or right singular vectors. Inspired by the idea mentioned above, herein we also exploit the JSVD-based to obtain the KSVs.

Estimating 2-D DOA. From
In the Appendix, we show how the KSVs are estimated from (13) and (16) using the ESPRIT-like algorithm. Note that unlike the standard ESPRIT extracting the signal subspace from the eigendecomposition of the covariance matrix of the concatenated measurements from the main array and its copy, here the KSVs are obtained by the JSVD of (13) and (16). In this way, we can detect the same number of signals as sensors, while the standard ESPRIT will not work in this case.
Having obtained the estimated KSVs D, we can estimate the 2-D DOA of each signal. Obviously, we can implement it by the 2-D "beamforming-like" technique, but searching over 2-D parameter space is a hard work. Therefore, we will suggest the following simple arithmetic to obtain the elevation and azimuth angles, and the kth elevation angle can be got where A denotes D(1 : p(M − p), k) and D(a : b, k) denotes the ath to bth elements of the kth column of D. In order to obtain the corresponding azimuth angle, herein defines a permutation matrix by where e i is the ith column of I p , we then have where B denotes D(M − p + 2 : p(M − p + 1), k).
Remark C. The pairing procedure of [15] is implemented by making good use of the diagonal elements of the cross-correlation matrix between different ULA embodying the corresponding pairs information of the elevation and azimuth angles. However, we find that the parameters with the pair-matching information are ω k = cos θ k − cos φ k k = 1, . . . , K. Therefore, it is easily verify when the difference among the parameters is very small or there is the same elevation (azimuth) angles and the azimuth (elevation) angle are very close, the Toeplitz matrix formed by the diagonal elements of the cross-correlation matrix may be a rankdeficient matrix so as not to obtain the correct pairing information for the elevation and azimuth angles [9]..Like the method in [8], our proposed can also remove these problems completely by the process of estimating every DOAs independently.

The Proposed Algorithm.
In summary, the proposed source estimation method with finite array data can be implemented as follows. (1) Calculate the estimated CCM and its reversal version between the received signal data x(t) and z(t)t = 1, . . . , N: (4) Estimate the KSV matrices D in (A.6) and D in (19), which obtain the K KSVs with one-to-one elevation and azimuth angles information.
(5) Estimate the kth 2-D angle {θ k , φ k } by implementing the kth KSV from (17) and (20). Note that 2-D angles are paired up automatically without extra association needed herein.

Simulation Results
In this section, we compare the estimation performance of the proposed method in estimating the 2-D DOAs and the detection probability of successful pair matching. We consider each ULA with elements spaced λ/2 apart and M = 5 sensors. Monte Carlo simulations based on 1000 independent realizations of the received data were carried out.
In the first test, we examine the DOA performance of two uncorrelated sources with equal power. Submatrices of size where T is the number of the independent trials. Figure 2 shows the RMSE of one of DOA [θ 2 , φ 2 ] = [55 • , 100 • ] estimation versus SNR for different subarray sizes by the proposed method. It can be seen from Figure 2 that RMSE of the elevation (azimuth) angle increases (decreases) with the increasing size of the submatrices. This is mainly due to the fact that the larger the number of dimensions to construct the KSV matrix, the better the alleviation of the noise by SVD to obtain the KSV matrix. Therefore, we will use the proposed method with the submatrix size of 3 in the subsequent experiments because of the largest dimensions. Performance evaluation of the proposed method, the CCM method based on ESPRIT (CCM-ESPRIT) [15], MUSIC [2] and the joint SVD (JSVD) in [8] is conducted along with the Cramer-Rao bound (CRB) [20]. The signals considered herein are the same as those in the first experiment.
To comprehensive comparison of performance of the 2-D DOA estimation, another objective criterion so-called joint  of root mean square error (JRMSE) defined as JRMSEs of the estimated elevation and azimuth angles versus SNR and the number of snapshots are plotted in Figures  3 and 4, respectively. Results of Figure 4 are produced at SNR = 10 dB. We observe that the proposed method delivers the most accurate DOA estimates among four methods, especially at low SNR. Similar results are obtained for smaller EURASIP Journal on Advances in Signal Processing number of snapshots. This may be due to the fact that the CCM-ESPRIT method only uses the CCM to overcome the pair-matching problems, but not for the DOA estimation, and though the JSVD method exploited the CCM to improve the estimation performance significantly, our technique makes good use of the structure of the CCM to make up the extended matrix with higher dimension for eliminating the effect of additive noise further. Performance of these methods at SNR higher than 10 dB is about the same and close to the CRB.
We also assess the estimation performance in terms of the detection probability of successful pair matching, as shown in Figure 5. The results indicate that the proposed method exhibits the similar detection rates to JSVD and higher detection rates than CCM-ESPRIT. It is worth noting that the CCM-ESPRIT method exhibits very low detection rate at low SNR. This phenomenon may be due to the fact that the Toeplitz matrix used in CCM-ESPRIT [15] to get pair-matching while the proposed and JSVD in [8] are paired automatically. Figure 6 shows the scatter plots of the elevation angle versus the azimuth angle for 2-D DOA estimation of five equal power uncorrelated signals from the direction [(5 • 140 • ), (15 •  . It is observed that the proposed method is able to deal with more sources than CCM-ESPRIT and JSVD which can estimate four signals at best, that is, the estimates of five signals approach to their respective true values as shown in Figure 6.

Conclusions
In this paper, we have presented a new automatic pairing method for 2-D DOA Algorithm using two orthogonal ULAs. The proposed method obtains an extended matrix by partition and stacking the CCM, and then estimates the KSVs with pair-matching information to get the one-toone elevation and azimuth angles. In addition, because of the matrix structural rearrangement, more array structural information are utilized to improve the performance of the DOA estimation. Simulation results have demonstrated that the proposed scheme not only provides more accurate estimates than the CCM-ESPRIT and JSVD, especially at low SNR, but also works well when sources equal to sensors of ULA.