# Enhanced 2D-DOA Estimation for Large Spacing Three-Parallel Uniform Linear Arrays

Authors: Dong Zhang; Yongshun Zhang; Cunqian Feng
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

An enhanced two-dimensional direction of arrival (2D-DOA) estimation algorithm for large spacing three-parallel uniform linear arrays (ULAs) is proposed in this paper. Firstly, we use the propagator method (PM) to get the highly accurate but ambiguous estimation of directional cosine. Then, we use the relationship between the directional cosine to eliminate the ambiguity. This algorithm not only can make use of the elements of the three-parallel ULAs but also can utilize the connection between directional cosine to improve the estimation accuracy. Besides, it has satisfied estimation performance when the elevation angle is between 70° and 90° and it can automatically pair the estimated azimuth and elevation angles. Furthermore, it has low complexity without using any eigen value decomposition (EVD) or singular value decompostion (SVD) to the covariance matrix. Simulation results demonstrate the effectiveness of our proposed algorithm.


---

## Body

## 1. Introduction

Two-dimensional DOA estimation has attracted extensive attention due to its wide range of applications [1–4], and there are many algorithms been proposed for DOA estimation. Among these algorithms, the subspace-based algorithms, such as MUSIC and ESPRIT, have received a lot of attention for its accurate estimation performance [5, 6]. However, the complexity of subspace-based algorithms is often too large because of the EVD or SVD. In [7], Marcos et al. put forward the well-known propagator method (PM) algorithm for 1D-DOA estimation. It uses linear partitioning instead of any EVD or SVD to reduce complexity. Then, Chen et al. [8] extend it to the DOA estimation of noncircular signal. References [9–11] extend the PM algorithm to 2D direction estimation. But these algorithms have some drawbacks. Reference [9] requires angle search operations. Reference [10] may fail in practical situation when elevation angle is between 70° and 90°. Reference [11] has worse estimation accuracy, because some element information is missing when calculating the propagation matrix (PMA). Based on this, Chen et al. [12] put forward an improved 2D angle estimation algorithm for three-parallel ULAs. It can solve all the problems mentioned in the above literature. But, it still has some shortcomings, that is, the array aperture is reduced because of the adoption of three parallel arrays with half wavelength spacing. And the estimation accuracy is reduced. References [13, 14] establish the “array of subarrays” idea and “cyclic ambiguity” idea to improve the estimation accuracy of 2D-DOA. According to this, we expanded the spacing between the ULA and proposed an enhanced 2D-DOA estimation algorithm. What needs to be stressed here is that our method is different from the method in [13, 14], although our algorithm and the method in [13, 14] both use the large array spacing to obtain high accuracy but ambiguous estimation. However, the principle of resolving ambiguity is different. References [13, 14] use coarse estimation without ambiguity or based on eigen space to resolve ambiguity. However, our algorithm uses the triangular relation between the three directional cosine to resolve the ambiguity. By doing this, we do not have to limit ourselves to using subarrays or eigen spaces to solve ambiguity problems. That is the main difference between our method and the method in [13, 14]. And it is also the main innovation of this paper. The flow of our algorithm is as follows: Firstly, we use the method in [12] to obtain the estimation of three directional cosines. Because of the large spacing between array elements, the estimated directional cosine is high precision and ambiguous. Then, we use the triangle relationship between the three directional cosines to eliminate the ambiguity. Then, we can get the true 2D-DOA of targets. Simulation results show that it cannot only avoid the problem in [8–11] simultaneously but also has better estimation accuracy than algorithm of [12] because of adopting large aperture. And the complexity of the algorithm is comparable to that algorithm of [12].Notations.
Superscripts⋅∗, ⋅T, ⋅†, ⋅−1, and ⋅H denote complex conjugation, transpose, pseudo-inverse, inverse, and conjugate transpose, respectively. IK is K×K identity matrix. diag⋅ denotes the diagonalization of the entity inside. ⋅ denotes take the absolute value of the element. ⋅ denotes round the element to the nearest integers less than or equal to the element. ⋅ denotes round the element to the nearest integers greater than or equal to the element. E⋅ denotes the expectation operation. arg⋅ and Re⋅ denote the phase and the real part of a complex number separately.


## 2. Problem Formulation

As shown in Figure1, assume that there are three-parallel ULAs, namely, X,Y, and Z. Array X contains N+1 sensors. Array Y and array Z have N sensors, respectively. The spacing between adjacent elements is dx. The distance between array X and array Y is dy, and the distance between array X and array Z is dz.λ is the wavelength. It is assumed that P far-field narrowband uncorrelated signals are incident onto the array. The elevation angle and azimuth angle of ith target are θi and φi, respectively. Here, we assume that the range of 2D angle is the same with that of [12]. That is to say that the range value of θi is 0,1/2π and the range value of φi is −1/2π,1/2π. Then, the output of the three ULAs can be expressed as follows:
(1)Xt=AxSt+Wxt,Yt=AyΩySt+Wyt,Zt=AyΩzSt+Wzt,where Ax=axθ1,φ1,axθ2,φ2,…,axθP,φP and axθi,φi=1,e‐j2πdxsinθisinφi/λ,…,e−j2πNdxsinθisinφi/λ.Ay contains the first N row of Ax.Ωy=diage−j2πdycosθ1/λ,e−j2πdycosθ2/λ,…,e−j2πdycosθP/λ and Ωz=diage−j2πdzsinθ1cosφ1/λ,e−j2πdzsinθ2cosφ2/λ,…,e−j2πdzsinθPcosφP/λ.Wxt,Wyt, and Wzt are assumed to be Gaussian white noise vectors whose mean value is zero and variance is σ2. Then a new vector Wt=XtTYtTZtTT, with L snapshots, and W=W1,W2,…,WL can be represented as
(2)W=XYZ=AxAyΩyAyΩzS+WxWyWz=AS+N,where A=AxTAyΩyTAyΩzTT,S=S1,S2,…,SL∈CP×L, and N=WxTWyTWzTT∈C3N+1×L.Figure 1
Sketch map of array structure.


## 3. Enhanced 2D-DOA Estimation

### 3.1. Highly Accurate but Ambiguous Estimation of Directional Cosine

Assuming thatdx=dy=d=λ/2,dz=Kd, where K is a positive integer and K>1. The difference between our model and the model in [12] is that the spacing between the ULA is larger than half wavelength. The large spacing will lead to ambiguity of the directional cosine, but it would not affect the estimation method for them. So the first several steps of our algorithm is same with the algorithm of [12]. We abbreviate the first several procedures as follows:Step 1.
Compute the covariance matrix ofWt by RW=EWWH.Step 2.
DivideRW into two parts, that is, RW=RW1RW2, where RW1∈C3N+1×P and RW2∈C3N+1×3N+1−P. Compute the propagator matrix (PMA) P by P^=RW1HRW1−1RW1HRW2 and get the extended PMA by Pe=IPHP^=PxT⏟P×N+1PyT⏟P×NPzT⏟P×NT.Step 3.
Implement EVD onΨz=Px1+Pz=A1ΩzA1−1 to get eigenvalues β^i and eigenvectors A1′. They correspond to the diagonal elements of Ωz and the estimation of A1, separately. And Px1+ is the first N rows of Px.Step 4.
LetB1=Px1A1′,B2=PyA1′. Attain α^i by extracting the ith diagonal elements of B1+B2. Similarly, let C1=Px1TPy1TPz1TTA1′ and C2=Px2TPy2TPz2TTA1′. Attain γ^i by extracting the ith diagonal elements of C1+C2. And Py1 and Pz1 are the first N−1 rows of Py and Pz.Py2 and Pz2 are the last N−1 rows of Py and Pz.Px2 represents the last N rows of Px.Until now, we get the paired estimation ofα^i,β^i, and γ^i which represent cosθi,sinθicosφi, and sinθisinφi, respectively. Because dz is larger than half wavelength, so the β^i is high accuracy but ambiguous. And all the ambiguous values are listed as follows:
(3)β^ik=β^i+kλdz,−1−β^idzλ≤k≤1−β^idzλ,where i=1,…,P.


### 3.2. Highly Accurate Estimation of 2D-DOA

Before estimating the 2D-DOA, we need to eliminate the ambiguity first. Note that there exists a triangle relationship between these directional cosines:(4)cosθi2+sinθicosφi2+sinθisinφi2=1.And the true directional cosine of targets satisfies this equation. So, we can use this equation to select the true directional cosine from the ambiguous directional cosine.But because there are square terms in the constraint, the mirror image of the true value also satisfies this equation. If the mirror image is exactly equal to an ambiguous value of the true value, we may make mistakes and treat the mirror image as the true value. Here, we set the range value ofθ and φ as 0,1/2π and −1/2π,1/2π, separately. So the true value of sinθicosφi is limited in the range of 0,1. And the mirror image can be avoided.And the procedure of eliminating ambiguity can be got as follows:Step 1.
Get all the ambiguous value ofsinθicosφi by
(5)β^ik=β^i+kλdz,0−β^idzλ≤k≤1−β^idzλ.Step 2.
Find the true value by(6)β^ik^=β^i+k^λdz,k^=argminkαi2+β^ik2+γi2−1.Lastly, we can get the 2D-DOA estimation for each source by(7)φ^i=tan−1argγ^iargβ^ik^,θ^i=tan−1argβ^ik^argα^icosφ^i.


### 3.3. Algorithm Analysis

Here, the analysis of the complexity, estimation performance, and some notices of our algorithm are as follows:(1)
Our algorithm has an additional ambiguity elimination operation compared with the algorithm of [12]. And the complexity of ambiguity elimination is far less than the complexity of constructing the covariance. Thus, we can say that the algorithm has comparable computational complexity with the algorithm of [12].
(2)
Because our array aperture is larger than that in [12], so our proposed algorithm can obtain a better estimation accuracy.
(3)
Note that here we expand the spacing between theX-axis and Z-axis. Similarly, we can also expand the spacing between the X-axis and Y-axis or simultaneously expand the spacing in the X-axis, Y-axis, and Z-axis. Then, we can get high accuracy but ambiguous estimation of sinθisinφi and cosθi, respectively. And we can still use (4) to eliminate the ambiguity. But considering the mirror image problem, the range of angle needs to be adjusted accordingly.
(4)
Because the first several steps of our algorithm are the same with those of [12], so our algorithm inherits the merits of algorithm of [12]. For example, (1) it can take full advantage of the elements of the three-parallel ULAs to estimate the PMA; (2) it has satisfied estimation accuracy in actual mobile elevation angles; (3) it can automatically pair the estimated 2D-DOA; (4) it has low complexity without using any EVD or SVD.


## 3.1. Highly Accurate but Ambiguous Estimation of Directional Cosine

Assuming thatdx=dy=d=λ/2,dz=Kd, where K is a positive integer and K>1. The difference between our model and the model in [12] is that the spacing between the ULA is larger than half wavelength. The large spacing will lead to ambiguity of the directional cosine, but it would not affect the estimation method for them. So the first several steps of our algorithm is same with the algorithm of [12]. We abbreviate the first several procedures as follows:Step 1.
Compute the covariance matrix ofWt by RW=EWWH.Step 2.
DivideRW into two parts, that is, RW=RW1RW2, where RW1∈C3N+1×P and RW2∈C3N+1×3N+1−P. Compute the propagator matrix (PMA) P by P^=RW1HRW1−1RW1HRW2 and get the extended PMA by Pe=IPHP^=PxT⏟P×N+1PyT⏟P×NPzT⏟P×NT.Step 3.
Implement EVD onΨz=Px1+Pz=A1ΩzA1−1 to get eigenvalues β^i and eigenvectors A1′. They correspond to the diagonal elements of Ωz and the estimation of A1, separately. And Px1+ is the first N rows of Px.Step 4.
LetB1=Px1A1′,B2=PyA1′. Attain α^i by extracting the ith diagonal elements of B1+B2. Similarly, let C1=Px1TPy1TPz1TTA1′ and C2=Px2TPy2TPz2TTA1′. Attain γ^i by extracting the ith diagonal elements of C1+C2. And Py1 and Pz1 are the first N−1 rows of Py and Pz.Py2 and Pz2 are the last N−1 rows of Py and Pz.Px2 represents the last N rows of Px.Until now, we get the paired estimation ofα^i,β^i, and γ^i which represent cosθi,sinθicosφi, and sinθisinφi, respectively. Because dz is larger than half wavelength, so the β^i is high accuracy but ambiguous. And all the ambiguous values are listed as follows:
(3)β^ik=β^i+kλdz,−1−β^idzλ≤k≤1−β^idzλ,where i=1,…,P.


## 3.2. Highly Accurate Estimation of 2D-DOA

Before estimating the 2D-DOA, we need to eliminate the ambiguity first. Note that there exists a triangle relationship between these directional cosines:(4)cosθi2+sinθicosφi2+sinθisinφi2=1.And the true directional cosine of targets satisfies this equation. So, we can use this equation to select the true directional cosine from the ambiguous directional cosine.But because there are square terms in the constraint, the mirror image of the true value also satisfies this equation. If the mirror image is exactly equal to an ambiguous value of the true value, we may make mistakes and treat the mirror image as the true value. Here, we set the range value ofθ and φ as 0,1/2π and −1/2π,1/2π, separately. So the true value of sinθicosφi is limited in the range of 0,1. And the mirror image can be avoided.And the procedure of eliminating ambiguity can be got as follows:Step 1.
Get all the ambiguous value ofsinθicosφi by
(5)β^ik=β^i+kλdz,0−β^idzλ≤k≤1−β^idzλ.Step 2.
Find the true value by(6)β^ik^=β^i+k^λdz,k^=argminkαi2+β^ik2+γi2−1.Lastly, we can get the 2D-DOA estimation for each source by(7)φ^i=tan−1argγ^iargβ^ik^,θ^i=tan−1argβ^ik^argα^icosφ^i.


## 3.3. Algorithm Analysis

Here, the analysis of the complexity, estimation performance, and some notices of our algorithm are as follows:(1)
Our algorithm has an additional ambiguity elimination operation compared with the algorithm of [12]. And the complexity of ambiguity elimination is far less than the complexity of constructing the covariance. Thus, we can say that the algorithm has comparable computational complexity with the algorithm of [12].
(2)
Because our array aperture is larger than that in [12], so our proposed algorithm can obtain a better estimation accuracy.
(3)
Note that here we expand the spacing between theX-axis and Z-axis. Similarly, we can also expand the spacing between the X-axis and Y-axis or simultaneously expand the spacing in the X-axis, Y-axis, and Z-axis. Then, we can get high accuracy but ambiguous estimation of sinθisinφi and cosθi, respectively. And we can still use (4) to eliminate the ambiguity. But considering the mirror image problem, the range of angle needs to be adjusted accordingly.
(4)
Because the first several steps of our algorithm are the same with those of [12], so our algorithm inherits the merits of algorithm of [12]. For example, (1) it can take full advantage of the elements of the three-parallel ULAs to estimate the PMA; (2) it has satisfied estimation accuracy in actual mobile elevation angles; (3) it can automatically pair the estimated 2D-DOA; (4) it has low complexity without using any EVD or SVD.


## 4. Computer Simulations

Here, we will conduct several simulation experiments to test the performance of our proposed method. The 2D-DOA of two uncorrelated targets areθ1,φ1=35°,40° and θ2,φ2=80°,75°, separately. N is set to 6. So the total number of array elements is 19. The spacing between X-axis and Z-axis is set to dz=3d. The other spacing is set to dx=dy=d.In the first simulation, Figure2 shows the scattergram of 300 times 2D-DOA estimate results by our algorithm and the algorithm of [12]. The snapshot is set to 50. The SNR = 10 dB. As we can see, both the two algorithms can clearly observe the two targets. But our algorithm has better statistical properties. That is because we make better use of the information between ULAs.Figure 2
2D-DOA scattergram at SNR = 10 dB.In the second simulation, we compare the estimation performance of the two algorithms under different SNR which vary from 0 dB to 30 dB. The snapshot is set to 300. 200 times Monte Carlo simulations are conducted under each SNR. As illustrated in Figure3, the RMSE curve of azimuth estimated by our algorithm is better than that estimated by the algorithm of [12]. The elevation estimation performance of our algorithm is comparable with that estimated by the algorithm of [12]. The reason is that when estimating azimuth angle, the equivalent array aperture of our algorithm is larger than that of [12]. While when estimating elevation angle, the equivalent array aperture of our algorithm is equal with that of [12]. Note that if we want to improve the estimation accuracy of elevation angle, we can expand dy instead of expandingdz.Figure 3
RMSE contrast of two algorithms with different SNR.
(a)
RMSE of azimuth estimation
(b)
RMSE of elevation estimationIn the third simulation, we compare the estimation performance of the two algorithms under different snapshots which vary from 100 to 2100. The SNR is set to 20 dB. 1000 times Monte Carlo simulations are conducted under each snapshot. As illustrated in Figure4, the RMSE curve of azimuth estimated by our algorithm is better than that estimated by the algorithm of [12]. The elevation estimation performance of the two algorithms is similar. The reason is the same with the previous one. From the second and third simulations, we can also see that our algorithm has satisfied estimation performance when elevation is between 70° and 90°.Figure 4
RMSE contrast of two algorithms with different snapshots.
(a)
RMSE of azimuth estimation
(b)
RMSE of elevation estimationIn the fourth simulation, we compare the complexity of our algorithm and that of [12]. We choose the runtime as the evaluation criteria. Two kinds of comparison were conducted. In the first comparison, the number of elements in subarray, that is, N, is set to vary from 10 to 100 and the snapshot is set to 300. The result is shown in Figure 5(a). In the second comparison, the snapshot is set to vary from 100 to 2100. The number of elements in subarray is set to 20. And the result is shown in Figure 5(b). From the two figures, we can see that the runtime of our algorithm is similar with that of [12] which has shown the efficiency of our proposed algorithm.Figure 5
Runtime comparison of two algorithms.
(a)
Runtime versus element number
(b)
Runtime versus number of snapshotsIn the fifth simulation, we compare the azimuth estimation performance of the two algorithms with different separation of angle. The 2D-DOA of two targets areθ1,φ1=30°,60°+Δφ1 and θ2,φ2=50°,70°+Δφ2, respectively. Both Δφ1 and Δφ2 vary from 0° to 18°. The SNR is equal to 20 dB. The snapshot is set to 500. 200 times Monte Carlo simulations are conducted under each separation of angle. Here, the “angle separation” denotes Δφ1 and Δφ2 which vary from 0° to 18°. It is different from the common definition of “angle separation.” In fact, the separation angle of the two targets is fixed, what changes is the azimuth angle of the two targets. The purpose of doing this is to make the azimuth angle traverse all angles between 70° and 88° to prove the effectiveness of the algorithm in practical application. The result is shown in Figure 6. We can see that the performance of our proposed algorithm is better than that of the algorithm of [12].Figure 6
RMSE of azimuth estimation versus angular separation.In the sixth simulation, we compare the elevation estimation performance of the two algorithms with different separation of angle. Here, we setdx=dz=d and dy=3d. The 2D-DOA of the two targets are θ1,φ1=60°+Δθ1,30° and θ2,φ2=70°+Δθ2,50°, respectively. Both Δθ1 and Δθ2 vary from 0° to 18°. The other parameters are the same with the fifth simulation. The result is shown in Figure 7. We can see that the performance of our proposed algorithm is better than that of the algorithm of [12].Figure 7
RMSE of elevation estimation versus angular separation.


## 5. Conclusions

An enhanced 2D-DOA estimation algorithm with three-parallel large spacing ULAs is proposed in this paper. It cannot only make use of the elements of the three-parallel ULAs to estimate the PMA but it can also utilize the information between three-parallel ULAs to improve the estimation accuracy. Besides, it has satisfied estimation performance when the elevation angle is between 70° and 90° and it can automatically pair the estimated 2D-DOA. Furthermore, it has low complexity without using any EVD or SVD to the covariance matrix. Simulations proved the effectiveness of the proposed algorithm.



---

*Source: 3679494-2018-03-04.xml*