Computationally Efficient 2D DOA Estimation with Uniform Rectangular Array in Low-Grazing Angle

In this paper, we propose a computationally efficient spatial differencing matrix set (SDMS) method for two-dimensional direction of arrival (2D DOA) estimation with uniform rectangular arrays (URAs) in a low-grazing angle (LGA) condition. By rearranging the auto-correlation and cross-correlation matrices in turn among different subarrays, the SDMS method can estimate the two parameters independently with one-dimensional (1D) subspace-based estimation techniques, where we only perform difference for auto-correlation matrices and the cross-correlation matrices are kept completely. Then, the pair-matching of two parameters is achieved by extracting the diagonal elements of URA. Thus, the proposed method can decrease the computational complexity, suppress the effect of additive noise and also have little information loss. Simulation results show that, in LGA, compared to other methods, the proposed methods can achieve performance improvement in the white or colored noise conditions.

For spatial smoothing techniques, Yeh et al. [15] developed the spatial smoothing 2D MUSIC algorithm by using the covariance matrices of overlapping rectangular subarrays and Chen et al. [16] presented an analysis of a special smoothing scheme extended in conjunction with the eigenstructure technique. To reduce the computational complexity, a partial spectral search based method (PSS) is proposed for limiting the searching region into a small sector [21]. Then, a tree structure one-dimensional (1D) algorithm [22] was developed based on a repeated use of the 1D MUSIC algorithm. However, it requires a large number of eigenvalue decompositions (EVDs) and does not perform well under a low signal-to-noise ratio (SNR). For spatial differencing techniques, Liu et al. [18] constructed a new spatial differencing matrix to suppress the white noise by using the difference between the first subarray and the spatial backward subarray. By exploiting the difference between the neighboring forward subarrays and backward subarrays, the method in [19] can suppress the colored noise. Additionally, an asymmetric spatial difference smoothing method [20] was used to reduce noise for coherent sources location, especially when the number of targets is odd.
However, due to the 2D peaks search, 2D or 1D EVD operations, all the aforementioned methods suffer from great computational complexity, especially for a large size of subarrays. In other words, these methods can provide better performance at the cost of great computation. Therefore, we concentrate on reducing the computational complexity caused by the 2D operations, while maintaining a high estimation performance. In this paper, we propose a spatial differencing matrix set (SDMS) method for 2D DOA estimation with URA in LGA. Employing the overlapping column or row subarrays along the xor y-direction, we build the SDMSs by rearranging the auto-correlation and cross-correlation matrices in turn among different subarrays. In addition, to suppress the data loss, we only perform difference for the auto-correlations and the cross-correlations are kept completely. Then, the two parameters are estimated independently by using the 1D subspace-based estimation technique, the pair-matching of which is achieved by extracting the diagonal elements of the URA. Simulation results verify the effectiveness of the proposed method.
The advantages of the proposed method can be given as follows: • The methods in [3,4,[15][16][17]21] involve the 2D EVD or 2D peak search, while the proposed method can estimate the parameters with 1D subspace-based estimation techniques.

•
The method in [22] can only use the auto-correlations of different subarrays, while the proposed method can use more information including auto-correlations and cross-correlations.

•
The spatial differencing techniques in [18,19] perform a difference operation on the whole subarrays, while the proposed method is only for the auto-correlations and the cross-correlations are kept completely. Thus, the SDMS method has little data loss.
The rest of this paper is given as follows. In Section 2, the basic received signal model of the URA in LGA is developed. Then, we derive the SDMS method using a 1D subspace-based estimation technique and achieve pair-matching by extracting the diagonal elements of the URA in Section 3. Simulation results are given in Section 4, and we conclude the whole paper in Section 5.
In this paper, operators (·) T , (·) * and (·) H represent transpose, conjugation, and conjugate transpose, respectively. I N denotes an N × N identity matrix and J M denotes an M × M exchange matrix with ones on its anti-diagonal and zeros elsewhere. ⊕ and ⊗ represent Hadamard product and Kronecker product, respectively; diag(·) and blkdiag(·) denote the diagonal matrix or the block diagonal matrix operator. E[·] and vec(·) denote expectation and vectorization, respectively.

System Model
As described in Figure 1, we regard the multipath effect as ideal specular reflection, where both curved earth and atmosphere refraction are not considered. We also consider K narrowband far-field signals s k (t) (k = 1, 2, · · · , K) impinging on a URA with MN well calibrated and identically polarized sensors parallel to the xoy plane. Here, both xand ydirections of the URA are separated by half a wavelength, the height of which is set as h. Since the received signals include two paths, i.e., direct path, reflected path, the output can be given as [11,12] where θ dk and θ rk are the direct and reflected elevation angle for the kth target (θ dk ≈ −θ rk = θ k ), respectively, α k is the azimuth angle and β k is the multipath reflection coefficient from the kth signal to receiver array; Let β k = exp[ j(π − 2π∆R k /λ)], ∆R k ≈ 2hsin θ k for simplicity, ∆R k is the difference value between direct path and reflected path in LGA; a x (α k , θ k ) = a x (u k ) = [1, e −jπu k , . . . , e −jπ(M−1)u k ] T , a y (α k , θ k ) = a y (v k ) = [1, e −jπv k , . . . , e −jπ(N−1)v k ] T , u k = sin θ k cos α k , and v k = sin θ k sin α k . The elements of Z(t) are temporally and spatially complex white Gaussian noises with zero-mean and variance σ 2 . Then, by vectorizing the matrix X(t), a composite data vector can be constructed as where Figure 1. The array geometry model for the URA.
By forming Φ y =diag[e −jπv 1 , e jπv 1 , · · · , e −jπv K , e jπv K ], the model in (2) can be rewritten as From (2) and (3), the incident signals can be divided into K parts and the signals in each part are correlated. With L snapshots (t = 1, 2, · · · , L), the sample covariance matrix can be calculated as where denotes the correlation matrix of coherent signals. Thus, due to rank-deficiency of the covariance matrix R 0 , the classic methods for 2D DOA estimation will lose efficiency [6][7][8].

2D DOA Estimation in LGA
In this section, we proposed a SDMS method that uses the differencing matrices among different spatial smoothing subarrays for 2D DOA estimation in LGA, where the parameters u k and v k are estimated independently.

1D Estimation of Parameter u k
Here, as described in Figure 2a, we divide the URA with M × N sensors into P overlapping forward rectangular subarrays of size Q × N along the x-direction. Each rectangular subarray has N column subarrays with Q sensors, such as the shaded areas for the first rectangular subarray. Then, the nth column of pth rectangular subarray can be given as z pn (t) = X p y n (t) = X p G n y(t), where As a result, using the matrix pencil of auto-correlation and cross-correlation matrices between different column subarrays, we can build the new matrix as  In (5), because of the calculation of covariance matrices in turn, the data information can be used more effectively. In addition, the cross-correlation matrix E[z pd z H pn (t)] (d = n) can also restrain the effect of additive noise. However, the auto-correlation matrix E[z np z H np (t)] has the noise covariance matrix σ 2 I Q , which will certainly decrease the performance. Then, in order to suppress the effect of noise, we first build an initial matrix that has the same noise matrix as R xp , and we have Combining (5) and (6), the forward SDMS for the x-direction (SDMS-x) can be defined as where 1 , e jπu 1 , e −jπu K , · · · , e jπu k ]. Following the forward backward (FB) technique, the new FB SDMS-x can be given as Then, based on the definition in (8), we can prove that the new FB SDMS-x has the following property. Theorem 1. Assume that there are 2K narrowband coherent signals impinging on the URA (M × N sensors). As described in Figure 2a, the URA is divided along the x-direction and the number of sensors in each column subarray is Q. Then, we define the new FB SDMS-x D x as in (8). If Q ≥ 2K, the rank of D x is equal to the number of the signals, namely, rank(D x ) = 2K.

Proof. See the Appendix A.
Under the Theorem 1, we divide A xQ into two submatrices as where A xQ1 and A xQ2 consist of the first 2K rows and the last Q − 2K rows, respectively. Since the matrix A xQ is a Vandermonde matrix with full rank, we can get a 2K × (Q − 2K) linear operator Γ x for A xQ2 = Γ H x A xQ1 . Then, the matrix D x can be divided into two submatrices as where D xQ2 = Γ H x D xQ1 . In addition, the operator Γ x can be calculated as Therefore, by constructing the matrix Ω , which can be used to estimate the parameter u k . By letting a (u k ) = [1, e −jπu k , · · · , e −jπ(Q−1)u k ] T , the parameter u k can be estimated by minimizing the following cost function where

1D Estimation of Parameter v k
Similarly, as described in Figure 2b, we divide the URA into F overlapping forward rectangular subarrays along the y-direction. Each rectangular subarray has M row subarrays with Q sensors. Then, the mth row of f th rectangular subarray can be set as w f m (t) = X f G m T e y(t), where T e is the row permutation matrix, X f = [0 Q×( f −1) I Q 0 Q×(F− f ) ], f = 1, 2, · · · , F, m = 1, 2, · · · , M, F=N-Q+1, and G m = [0 N×(m−1)N I N 0 N×(M−m)N ]. Then, we can build the new FB SDMS for the y-direction (SDMS-y) as where and A yQ is the submatrix of the array response matrix A y consisting of the first Q rows, . Likewise, the new FB SDMS-y in (12) has the following property.
Theorem 2. Assume that there are 2K narrowband coherent signals impinging on the URA (M × N sensors). As described in Figure 2b, the number of sensors in each row subarray is Q, and the new FB SDMS-y is defined in (12). If Q ≥ 2K, then the rank of D y is equal to the number of the signals, namely, rank(D y ) = 2K.
Based on the Theorem 2, we can divide the matrix D y as where D yQ2 = Γ H y D yQ1 and Γ y is the linear operator. By letting a (v k ) = [1, e −jπv k , · · · , e −jπ(Q−1)v k ] T , the parameter v k can be estimated by minimizing the following cost function where

Pair-Matching of Parameters u k and v k
Since the estimated parameters u k and v k are calculated independently, the pair-matching is very important for multiple targets. In the case of M<N, by extracting the diagonal elements of URA, we can write where A (v k ) = A x , A (u k ) is the submatrix of the array response matrix A y consisting of the first M rows, and N 0 is the subarray of N. Since both A (v k ) and A (u k ) are Vandermode matrices, is also a Vandermode matrix. Using the FBSS method, we can divide r diag (t) into P overlapping forward subarrays with Q sensors. Then, the smoothing matrix can be given as where . Hence, we can get the estimated parameters by minimizing the following cost function where , and D d1 and D d2 consist of the first 2K rows and the last Q-2K rows of D diag , respectively. From (20), the parameters u k and v k can be matched by repeating the following minimization for i = 1, · · · , 2K and j = 1, · · · , i − 1 and the constraint condition k i = k j can avoid the different u k paired with the same v k . Then, the azimuth angle and elevation angle can be estimated as

Implementation of the Proposed Method
From (11), (17) and (22), the azimuth and elevation angles can be estimated by using the spatial differencing method, where the computational burdensome EVD is avoided and the effect of additive noise is also suppressed. Then, with the finite array data y(t) L t=1 , the proposed method can be implemented as follows: 1.

2.
Form the FB SDMS-xD x in (8) and the FB SDMS-yD y in (12) aŝ whereD xp andD y f can be calculated by using the covariance matrixR 0 2Q 3 flops.

3.
Estimate the orthogonal projectorsΠ x in Section 3.1 andΠ y in Section 3.2.

4.
Estimate the parameters u k and v k by finding the phases of the p zeros of the polynomial p u (z) and p v (z) using (11) and (17), where p u (z)

5.
Perform the pair-matching of the parameters v k and u k by using (18)- (21) and estimate the azimuth and elevation angles by (22) As shown above, the cost of each step is roughly indicated in terms of the number of MATLAB flops. Furthermore, in the case of Q 2K, the computational complexity of the proposed method mainly includes the calculation of covariance matrix and SDMSs, the estimation of parameters and the pair-matching, which is about Remark 1. The proposed method can estimate the parameters independently, and the null space is calculated by using the linear propagator based on the partition of the array response matrix. The polynomial roots can be obtained by the Linsey-Fox root finding algorithm, which is less than the MATLAB function roots. However, the 2D FBSS-MUSIC method in [15,16], the 2D FBSS based DOA matrix (FBSS-DOAM) method in [23], the conventional spatial differencing (CSD) method in [18,19], and the tree structure one-dimensional (1D) based (TSOD) algorithm in [22] all involve the EVD to obtain the signal subspace or noise subspace. Furthermore, as shown in Table 1, FBSS-MUSIC and CSD both need 2D spectrum peak searching. TSOD needs (2P + 8K + 1) 1D EVD, while FBSS-DOAM needs two times of 1D spectrum peak searching. It is easily seen that the proposed method is computationally more efficient than other recently developed methods.

FBSS-MUSIC FBSS-DOAM CSD TSOD Proposed Method
As described in (7) and (12), the proposed method can form the row and column FB SDMSs by using the auto-correlation and cross-correlation matrices among different subarrays, where we only perform the difference for the auto-correlations, and the cross-correlations are kept completely. In addition, we also employ the FB technique to improve the estimation performance. However, the TSOD method can only use the auto-correlations of different subarrays and CSD performs the difference on the whole subarrays. Thus, SDMS can achieve performance improvement over the methods in [18,19,22].

Cramér-Rao Bounds (CRB)
As described in Section 2, according to [24], the CRB can be obtained as where

Simulation Results
We now evaluate the estimation performance of SDMS by using some numerical experiments. Experiment 1: Effectiveness of proposed method. In this experiment, we mainly examine the effectiveness of SDMS in the presence of white noise and colored noise conditions, where the number of sensors in each subarray is Q = 16 and the total number of snapshots is chosen to be L = 500. Here, the SNR is set as 10 dB. Figure 3 shows the estimation results of the proposed method with 100 Monte Carlo runs. As expected, all the 2D DOAs can be estimated effectively and accurately for the white noise and colored noise conditions. Experiment 2: Performance versus SNR. In Figures 4 and 5, we evaluate the performance in terms of SNR in the white noise and colored noise conditions, where we assume L = 500 and Q = 16. Three methods are performed for comparison, including the CSD method in [18], the TSOD method in [22] and the proposed method. Moreover, the CRB is also provided. It can be observed that, due to the use of more information and the application of the difference operation, the performance of the proposed method is better than those of methods in [18,22]. In the white noise condition, the performance of CSD is better than TSOD for relatively low SNR, whereas it is the opposite with high SNR. It illustrates that the difference operation can reduce the effect of white noise in the low SNR condition, while the data loss will also decrease the performance in a high SNR condition. In the colored noise condition, since the non-diagonal elements of the noise covariance matrix have significant value, CSD is always superior to TSOD. Comparing SDMS and CSD, SDMS can achieve better performance by using more data information. To sum up, SDMS can achieve great performance improvement in the low SNR condition by using more information and performing the difference operation.   Experiment 3: Performance versus the number of snapshots. In this experiment, we evaluate the estimation performance with different methods (described in Experiment 2) in terms of the number of snapshots, where we assume the SNR is 2 dB and Q = 16. Figures 6 and 7 show the RMSE versus the number of snapshots in the white noise and colored noise conditions, respectively. It is shown that, because of the use of difference operation and cross-correlation matrices, even when the number of snapshots is small, the proposed method still outperforms the methods in [18,22]. Then, for the white noise, the performance of CSD is weaker than TSOD for the small number of snapshots due to the data loss caused by the difference operation between forward and backward smoothed matrices. For colored noise, CSD performs better than TSOD. In addition, by using more data information, the performance of SDMS is also less sensitive to the number of snapshots than that of CSD.    Figure 8 describes the performance versus the size of subarrays in the white noise condition. Here, we assume L = 500, SNR = 10 dB. It is easily seen that SDMS still outperforms CSD and TSOD. In addition, due to the increase of array aperture, the performance of those methods becomes better. Then, we can observe that TSOD fails with a small size of subarrays but performs better than CSD with a big size. The reason is that the data loss caused by the difference operation will become larger when the size of subarrays increases. However, the proposed method can use the cross-correlation to reduce the data loss, resulting in a better performance.

Conclusions
A new computationally efficient SDMS method with little information loss is proposed to suppress the effect of white noise or colored noise in LGA. The two parameters are estimated independently by using a 1D subspace-based estimation technique, the pair-matching of which is achieved by extracting the diagonal elements of URA. Simulation results show that, in LGA, the performance of the proposed method is superior to the other recently developed method in low SNR conditions and with a small number of snapshots for white noise and colored noise conditions. In the near future, based on improving the information utilization, the extended spatial differencing method with MIMO radar for joint DOA and DOD estimation will be considered carefully.