# Two-Dimensional DOA Estimation in Compressed Sensing with Compressive-Reduced Dimension-lp-MUSIC

Authors: Weijian Si; Xinggen Qu; Lutao Liu; Zhiyu Qu
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

This paper presents a novel two-dimensional (2D) direction of arrival (DOA) estimation method in compressed sensing (CS) to remove the estimation failure problem and achieve superior performance. The proposed method separates the steering vector into two parts to construct two corresponding noise subspaces by introducing electric angles. Then, electric angles are estimated based on the constructed noise subspaces. In order to estimate the azimuth and elevation angles in terms of estimates of electric angles, arc-tangent operations are exploited. The arc-tangent is a one-to-one function and allows the value of the argument to be larger than unity so that the proposed method never fails. The proposed method can avoid pair matching to reduce the computational complexity and extend the number of snapshots to improve performance. Simulation results show that the proposed method can avoid estimation failure occurrence and has superior performance as compared to existing methods.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation has been a topic of great interest in many fields such as radar, mobile communication systems, and medical imaging [1, 2]. Many conventional DOA estimation methods have been developed in the last years. Among them, multiple signal classification (MUSIC) [3] and estimation of signal parameter via rotational invariance technique (ESPRIT) [4] are regarded as the most popular methods that have high resolution for DOA estimation. However, the performance of these methods can severely degrade when signal to noise ratio (SNR) is low, number of snapshots is small, or sources are coherent.To overcome the aforesaid drawbacks, DOA estimation involving compressed sensing (CS) [5, 6] is investigated due to the development of methods based on CS. The CS-based estimation methods enforce sparsity on the spatial spectrum and perform source localization in an overcomplete dictionary. Malioutov et al. [7] propose the l1-SVD method for DOA estimation, which casts DOA estimation problem as a sparse recovery problem. Stoica et al. [8] present a sparse iterative covariance-based estimation (SPICE) method by exploiting the covariance matching criterion. In [9], an alternative strategy called joint l0 approximation method is proposed to resolve closed spaced and high coherent sources even if the number of sources is unknown.Although one-dimensional (1D) DOA estimation has attracted tremendous interest, two-dimensional (2D) DOA estimation is of more practical importance. Since L-shaped array configuration has higher accuracy than other array configurations [10] such as the parallel uniform linear array (ULA) configuration [11], the rectangular array configuration [12], and the circular array configuration [13], most existing 2D DOA estimation methods are proposed based on the L-shaped array. Tayem and Kwon [14] propose the propagator method (PM) with one or two L-shaped array configurations to remove nonnegligible drawbacks of PM with parallel shape array configuration, but PM with one L-shaped array configuration still has estimation failure problem; that is, estimation method cannot perform estimation correctly in the entire ranges of the azimuth and elevation angles. Based on the L-shaped array configuration, Liang and Liu [15] propose a novel joint azimuth and elevation angles estimation method, which avoids pair matching. It has been proven that conventional 2D-MUSIC [16] method can provide a precise 2D estimation, but the requirement of 2D search needs high computation complexity. To reduce the computational load, a sparse L-shaped array configuration [17] is used, where the azimuth and elevation angles are estimated by the shift-invariance property of ULA and modified total least squares (MTLS) techniques, respectively. In [18], Wang et al. propose 2D-l1-SVD and enhanced 2D-l1-SVD methods, which have several advantages over conventional methods.In this paper, a novel compressive-reduced dimension-lp-MUSIC method called CS-RD-lp-MUSIC, which requires no pair matching, is proposed for 2D DOA estimation in CS. The key idea of the proposed method is that the steering vector is separated into two parts to construct two corresponding noise subspaces by introducing electric angles. Then, based on the constructed noise subspaces, electric angles are estimated by the proposed method, where CS-MUSIC is employed to estimate one electric angle and RD-lp-MUSIC is adopted for other electric angles. What is more, CS-MUSIC can improve performance significantly even if covariance matrix tends to lose rank and RD-lp-MUSIC can reduce computational complexity by reduced dimension (RD). Our objective in this paper is to estimate the azimuth and elevation angles based on estimates of electric angles by arc-tangent operations. The arc-tangent is a one-to-one function and allows the value of its argument to be larger than unity in order that the proposed method never fails. We show that the proposed method can remove the estimation failure problem and achieve good performance due to the application of CS and extension of the number of snapshots. In addition, the proposed array configuration can further remove the estimation failure problem without loss of performance of DOA estimation. Simulation results illustrate the superior performance of the proposed method with comparisons to existing methods.Notations used in this paper are given as follows. Lowercase boldface italic letters are served for vectors and uppercase boldface italic letters are served for matrices.·∗, ·T, and ·H denote the complex conjugate, transpose, and conjugate transpose, respectively. ·-1 denotes the inversion of square matrix or pseudoinversion of nonsquare matrix. ·2 and ·p denote the Euclidean norm for vectors and lp-norm for matrices, respectively. Ai, Ai, and Ai,j denote the ith column, ith row, and (i,j)th entry of matrix A, respectively. diag⁡A is a diagonal matrix with the diagonal elements of matrix A and diag⁡a is a diagonal matrix with a being its diagonal elements. A(K)=E11K,E22K,…,EKKKT is a K2×K matrix, where Eii(K),i=1,2,…,K, is a K×K matrix with one in the (i,i)th entry and zeros elsewhere. eK(i) is a K×1 vector with one in the ith element and zeros elsewhere. IK×K is a K×K identity matrix. 0K×K is a K×K matrix of all zeros. ⊗ and ⊙ denote the Kronecker product and Hadamard product, respectively. unvec⁡· denotes the matrix form of a vector. Re⁡· and Im⁡· denote the real and imaginary parts of a complex variable, respectively.The remainder of the paper is organized as follows. In Section2, we formulate the 2D DOA estimation problem, in which the steering vector is separated into two parts. The proposed method is described in detail in Section 3. Section 4 shows the performance of the proposed method and Section 5 concludes the paper.


## 2. Problem Formulation

The proposed array configuration that consists of three uniform linear arrays (ULAs) ofM sensors with intersensor spacing d is shown in Figure 1. One ULA lies in the x-z plane, another lies in the y-z plane, and the last one lies on the z-axis. The origin of the array configuration is set as the referencing sensor. Consider the array configuration impinged on K narrowband far-field sources, skt,   k=1,2,…,K, where the kth source has the azimuth angle αk and elevation angle βk shown in Figure 1.Figure 1
Proposed array configuration for 2D DOA estimation.The coordinate of the(ix,iy,iz)th sensor is (ixd,iyd,izd) so that the received source on the (ix,iy,iz)th sensor at the tth snapshot can be given by(1)xix,iy,izt=∑k=1Kskte-j2πdcos⁡βkcos⁡αkix+cos⁡βksin⁡αkiy+sin⁡βkiz/λ+nix,iy,izt=∑k=1Ksktejθkixejφkiyejϕkiz+nix,iy,izt,where nix,iy,izt denotes the noise term, λ is the wavelength, j=-1, and three electric angles θk, φk, and ϕk [15], which are the functions of αk and βk, are, respectively, defined as θk=-2πdcos⁡βkcos⁡αk/λ, φk=-2πdcos⁡βksin⁡αk/λ, and ϕk=-2πdsin⁡βk/λ. Since (ix,iy,iz)∈{(1,0,0),(0,0,0),(0,1,0),…,(1,0,M-1),(0,0,M-1),(0,1,M-1)}, the vector form of (1) can be written as the following form:(2)xt=Ast+ntt=1,2,…,T,where x(t)=[x1,0,0(t),x0,0,0(t),x0,1,0(t),…,x1,0,M-1t, x0,0,M-1t,x0,1,M-1(t)]T and n(t)=[n1,0,0(t),n0,0,0(t), n0,1,0t,…,n1,0,M-1t,n0,0,M-1t,n0,1,M-1(t)]T are 3M×1 vectors, s(t)=[s1,s2,…,sK]T, A=[a(θ1,φ1,ϕ1),a(θ2,φ2,ϕ2),…,a(θK,φK,ϕK)] is a 3M×K manifold matrix, and a(θk,φk,ϕk)=[ejθk,1,ejφk,ejθkejϕk,…,ejθkej(M-1)ϕk,ej(M-1)ϕk,ejφkej(M-1)ϕk]T is a 3M×1 steering vector. By the properties of ⊗ and ⊙, the steering vector a(θk,φk,ϕk) can be simplified to(3)aθk,φk,ϕk=a~ϕk⊗a˘θk⊙a⌢φk=a~ϕk⊗a^θk,φk,where a˘(θk)=[ejθk,1,1]T, a⌢(φk)=[1,1,ejφk]T  a^(θk,φk)=a˘(θk)⊙a⌢(φk)=[ejθk,1,ejφk]T, and a~(ϕk)=[1,ejϕk,…,ej(M-1)ϕk]T. Denote A~(ϕ)=[a~(ϕ1),a~(ϕ2),…,a~(ϕK)] and A^(θ,φ)=[a^(θ1,φ1),a^(θ2,φ2),…,a^(θK,φK)] as the manifold matrices that contain information about ϕ and (θ,φ), respectively. The purpose of this simplification is to construct two corresponding noise subspaces to estimate the electric angles. Then, the matrix form of (2) is given by(4)X~t=unvec⁡xt=A~ϕGt+N~tt=1,2,…,T,where G(t)=[s1a^(θ1,φ1),s2a^(θ2,φ2),…,sKa^(θK,φK)]T is a K×3 source matrix and X~(t) and N~(t) are M×3 matrices whose elements in the (m,n)th entry are the (3m+n-3)th elements of vectors x(t) and n(t), respectively. It is easy to see from (4) that ϕ and (θ,φ) span the column and row spaces of x~(t), respectively, and the manifold matrix A~(ϕ) is determined by ϕ while the sources matrix G(t) is decided by (θ,φ). The electric angles estimation can be performed by estimating ϕ and (θ,φ) successively instead of jointly estimating them. Since ϕ is separated from (θ,φ), ϕ can be firstly estimated and then (θ,φ) can be obtained in terms of the estimates of ϕ.Let{ϕ^p}p=1Pϕ be a sampling grid which covers the entire spatial domain so that the true electric angles {ϕk}k=1K are on the sampling grid set where Pϕ  (Pϕ≫K) denotes the number of the sampling grid. This means that if ϕ^p1,ϕ^p2,…,ϕ^pK are one to one corresponding to the true electric angles ϕ1,ϕ2,…,ϕK, we have(5)gp=skp=pkk=1,2,…,K0elsewhere.By denotingX-≜[X~(1),X~(2),…,X~(T)] as M×3T data matrix, (4) can be rewritten as the following sparse form:(6)X-=A-ϕG-+N~1,N~2,…,N~T=A-ϕG-+N-,where A-(ϕ)=[a~(ϕ1),a~(ϕ2),…,a~(ϕPϕ)] is the M×Pϕ manifold matrix corresponding to all potential electric angles which is also defined as an overcomplete dictionary in CS. Since {G-t}t=13T share the same support based on joint sparsity, G-=[G-1,G-2,…,G-Pϕ]T has K nonzero rows, that is, the row K-sparse, where G-pϕ=[gpϕ(1)A^kT(θ,φ),gpϕ(2)A^kT(θ,φ),…,gpϕ(T)A^kT(θ,φ)]T is a 3T×1 vector. To estimate the electric angle, the support needs be determined by the matrix X which is given by(7)X=ΦA-ϕG-+ΦN-=A⌢ϕG-+N,where Φ is the N×M measurement matrix and N is the number of nonadaptive linear projection measurements.


## 3. DOA Estimation

In this section, a 2D DOA estimation problem is solved by two steps in the CS scenario. Based on the above source model, a 1D DOA estimation is firstly performed to estimateϕ and get the information about (θ,φ) which is contained in the source matrix of (7). Secondly, based on the estimates of ϕ, two electric angles θ and φ can be estimated by minimizing the relaxation of the residual fitting error [19](8)Jp,qY=min⁡∑i=1P∑j=1Qyijqp/q,where Y∈CP×Q is assumed to be the residual fitting error matrix. For simplicity, the case q=2 is considered in the rest of this paper; that is,(9)Jp,2Y=min⁡∑i=1P∑j=1Qyij2p/2.The motivation to chooseq=2 is the following two main reasons. As can be seen, the first reason is that the number of nonzero rows can be provided as p approaches zero. A nonzero row can serve as a penalty factor with the reduction of p which promotes a sparse frame among all rows. Secondly, practical issues such as computational complexity are also in favor of this choice. In the case q=2, a low computational complexity is obtained by minimizing (9) instead of (8) [20]. Then, for notational convenience, we denote J(p,2)(Y) by J(p)(Y) which is called lp-norm.


### 3.1. CS-RD-lp-MUSIC

Step 1 (estimateϕ by CS-MUSIC).
CS-MUSIC, which is an extension of MUSIC, can identify the parts of support using CS-based methods, after which the remaining supports are estimated by the generalized MUSIC criterion. The main contribution of CS-MUSIC is to overcome the error caused by losing rank which may cause disastrous consequences in conventional MUSIC.Letsupp⁡G-={1≤g≤Pϕ:G-g≠0} and R(A⌢(ϕ)) denote the support of G- and the range space of A⌢(ϕ), respectively. Due to (4), the number of snapshots is extended to 3T which is one of the advantages of the proposed method. It is obvious that CS-MUSIC can be simplified as MUSIC for K≤3T. However, CS-MUSIC can estimate DOA with success but MUSIC fails in the K>3T case. In CS-MUSIC, if K>3T, K-3T indices of supp⁡G- are determined by CS-based methods such as simultaneous orthogonal matching pursuit (SOPM) [21] so that column vector space R(A⌢ηK-3T(ϕ)) is decided, where ηK-3T is the set of indices and A⌢ηK-3T(ϕ) is the submatrix of A⌢(ϕ) with columns indexed by ηK-3T. The remaining 3T indices of supp⁡G- are obtained by the generalized MUSIC criterion. To make CS-MUSIC applicable for all range of the snapshot, two new orthogonal spaces R(A⌢G-) and R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))) are constructed where PR(Q) is the orthogonal projection onto the noise subspace R(Q). Then, the electric angle ϕ can be estimated by the spectrum search. Now, the major steps of CS-MUSIC are summarized as follows.(1)
FindK-3T indices of supp⁡G- by SOMP and let ηK-3T be the set of indices.(2)
DetermineA⌢ηK-3T(ϕ) in terms of ηK-3T and construct the noise subspace of CS-MUSIC R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))).(3)
Fori∈{1,2,…,Pϕ}, calculate the spatial spectrum, Pspectral=1/A⌢iH(ϕ)[PR(Q)-PR(QQHA⌢ηK-3T(ϕ))]A⌢i(ϕ), and then estimate ϕ by the locations of K highest peaks of the spatial spectrum.Step 2 (estimateθ and φ by RD-lp-MUSIC).
In this step, RD-lp-MUSIC is exploited to estimate electric angles θ and φ in terms of the estimates of ϕ. Let {ϕ-1,ϕ-2,…,ϕ-K} denote the locations of K highest peaks. Then, the row number corresponding to locations of K highest peaks is denoted as {i1,i2,…,iK} and these rows contain information about (θ,φ). It is clear that the ith row vector of source matrix [g⌢1,1(i),g⌢2,1(i),…,g⌢T,1(i),…,g⌢1,3(i),g⌢2,3(i),…,g⌢T,3(i)],  i=i1,i2,…,iK, is a 1×3T vector and thus the matrix form G⌢(i) is a T×3 matrix whose element in the (m,n)th entry is g⌢m,n(i). Denote A~(θ,φ)=A^T(θ,φ) and S⌢=[s(1),s(2),…,s(T)]T so that the lp-norm minimization problem can be expressed as(10)Jpθ,φ=minθ,φ⁡JpG⌢i-S⌢A~θ,φ=minθ,φ⁡∑t=1TG⌢ti-S⌢tA~θ,φ2pi=i1,i2,…,iK.Note thatK  lp-norm minimization problems have the same structure and can be solved in a similar manner. Hence, we only need to find the optimal solution of one minimization problem. The remaining K-1 minimization problems can be handled in the same way.Since the objective function (10) requires an exhaustive 2D search, high computational cost is needed for precise estimation which can result in the reduction of algorithmic efficiency. To avoid heavy computational load, RD-lp-MUSIC is proposed for 2D estimation just through 1D search. A detailed derivation process of exploiting RD-lp-MUSIC for estimating θ and φ is given as follows.DenoteF(θ,φ)=G⌢t(i)-S⌢tA~(θ,φ)22 so that RD can be realized by minimizing F(θ,φ). By the property of ·2, we have(11)Fθ,φ=G⌢ti-S⌢tA~θ,φG⌢ti-S⌢tA~θ,φH=const+S⌢tA⌣θ⊙A⌢φA⌣θ⊙A⌢φHS⌢tH-2Re⁡S⌢tA⌣θ⊙A⌢φG⌢tiH,where A⌣(θ)=[a˘(θ1),a˘(θ2),…,a˘(θK)]T and A⌢(φ)=[a⌢(φ1),a⌢(φ2),…,a⌢(φK)]T. Denote that P(K)=[E11(K),E22(K),…,EKK(K)]T and P(3)=[E11(3),E22(3),E33(3)]T are K2×K and 9×3 selection matrices, respectively, so that the following equation holds based on selection matrices [22]:(12)A⌣θ⊙A⌢φ=PKTA⌣θ⊗A⌢φP3=PKTA⌣θ⊗IK×KI3×3⊗A⌢φP3.During the aforesaid process, selection matrices are of the greatest importance. The purpose of utilizing selection matrices is to separateθ from φ and reduce dimension. Moreover, based on the following equation,(13)Q1I3×3⊗A⌢φQ2=A⌢φ,where Q1=[IK×K0K×2K] and Q2=[I3×303×6]T, (12) can be further rewritten as (14)A⌣θ⊙A⌢φ=PKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3.Subsequently, by substituting (14) into (11), we have(15)Fθ,φ=const-2Re⁡S⌢tPKTA⌣θ⊗IK×K·Q1-1A⌢φQ2-1P3G⌢tiH+S⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3P3HQ2-1HA⌢Hφ·Q1-1HA⌣Hθ⊗IK×KPKTHS⌢tH.To eliminate trivial solutionsA⌣(θ)=0K×3 and A⌢(φ)=0K×3, a constraint condition eKH(1)[A⌣(θ)⊙A⌢(φ)]e3(2)-1=0 is considered. Hence, the cost function L(θ,φ) is given as the following form:(16)Lθ,φ=S⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φ·Q2-1P3P3HQ2-1HA⌢HφQ1-1HA⌣Hθ⊗IK×KPKTHS⌢tH-2ReS⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φ·Q2-1P3G⌢tiH-μeKH1·PKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3e32-1,where μ is a constant. By setting the partial derivation of (16) with respect to A⌢(φ) to zero, we have(17)∂Lθ,φ∂A⌢φ=2WS⌢tHS⌢tWA⌢φVVH-2Re⁡WHS⌢tHG⌢tiHVH-μWHeK1e3H2VH=0,where W=P(K)T[A⌣(θ)⊗IK×K]Q1-1 and V=Q2-1P3. By (17) and the constraint condition eKH(1)[A⌣(θ)⊙A⌢(φ)]e3(2)=1, μ can be expressed as(18)μ=2-2eKH1S⌢t-1Re⁡G⌢tie32eKH1S⌢t-1S⌢t-1HeK1.Therefore, it can be deduced from (17) and (18) that (19)A⌢φ=W-1S⌢t-1ReG⌢ti+1-eKH1G⌢ti-1Re⁡G⌢tie32eK1e3H2eKH1S⌢t-1S⌢t-1HeK1·V-1.As one may note,A~(θ,φ) is transformed into U⌢(θ) by inserting A⌢(φ) into A~(θ,φ) and thus RD is realized so that 2D DOA estimation just requires 1D search. Furthermore, RD can avoid the failure occurrence in separating θ from φ. As a matter of fact, RD also provides an important clue for expounding the relationship between θ and φ. Therefore, RD is the foundation of estimating electric angles θ and φ by RD-lp-MUSIC. Then, we focus on the solver for the last remaining problem to estimate θ; that is,(20)Jpθ=minθ⁡∑t=1TG⌢ti-S⌢tU⌢θ2p.By denotingrt=G⌢t(i)-S⌢tU⌢(θ) as the residual vector, the objective function f(rit) can be expressed as (21)frit=∑t=1T∑i=13rit2p/2.Based on the following equation,(22)∂f∂rit∗=12∂f∂Re⁡rit+j∂f∂Im⁡rit,the partial derivative with respect to (rit)∗ can be given by(23)∂f∂rit∗=p2∑t=1Trit2p/2-1rit.Due to (23), the gradient of f with respect to rt is given by (24)∂f∂rt=∂f∂r1t∗∂f∂r2t∗∂f∂r3t∗=p2rt⊙Λ,where(25)Λ=∑t=1Tr1t2p/2-1∑t=1Tr2t2p/2-1∑t=1Tr3t2p/2-1.Denote F=diag⁡(Λ) so that (24) can be rewritten as(26)∂f∂rt=p2rtF.Note thatF is a positive definite matrix clearly. Then, since the second-order partial derivative of f with respect to rit is given as the following form,(27)∂2f∂ritT∂rjt∗=p2∑t=1Trit2p/2-1p2-1∑i=13rit2-1rit2+1i=j0i≠j,the 3×3 Hessian matrix of f with respect to rt is given by(28)Hrtrt∗=∂2f∂rtT∂rt∗=p2·Fp2-1∑i=13rit2-1diag⁡rtTrt+I3×3.Denote the gradient off with respect to U⌢(θ) as Γ=[Γ1,Γ2,Γ3] so that Γi is computed as(29)Γi=∂f∂U⌢i∗=-S⌢tT⊗e3Tip2rtFT.Then, the Hessian matrix off with respect to U⌢(θ) is (30)HU⌢∗U⌢=∂2f∂U⌢∗∂U⌢T=S⌢tT⊗ΣP3Hrtrt∗P3TS⌢t∗⊗ΣT,where Σ=[e3(1)e3(2)e3(3)]. By applying the Newton method, the following sequence U⌢(z)(θ), z=0,1,2,…, is obtained to find the minimization of f:(31)U⌢z+1θ=U⌢zθ+σzΔzU⌢zθ,where σz is a positive step and Δ(z)(U⌢(z)(θ))=-HU⌢∗U⌢Γ(U⌢(z)(θ)). The initial value of U⌢(z)(θ) is defined as U⌢(0)(θ)=[(S⌢t)HS⌢t]-1(S⌢t)HG⌢t(i). Specifically, the iteration is terminated if the following stopping criterion is satisfied:(32)Jpθz-Jpθz+1Jpθz<εfor some small ε. In the simulations, ε is set to 10-7. After determining U⌢(θ), the projection matrix onto the noise subspace is given by(33)P⌢=I3×3-UθUHθUθ-1UHθ,where U(θ)=U⌢H(θ). Then, by computing the spatial spectrum, Pspetral(θ)=1/a˘H(θ)P⌢a˘(θ), and searching its peaks, we can provide an accurate estimation for θ. Once θ is obtained from the spatial spectrum, φ is easily estimated in terms of (19). Therefore, RD-lp-MUSIC avoids pair matching.After electric angles are estimated, the azimuth angleαk and elevation angle βk of the kth source can be given by αk=arctan⁡φk/θk and βk=arctan⁡ϕksin⁡αk/φk, respectively. The major steps of RD-lp-MUSIC for estimating θ and φ are given as follows.(1)
Thelp-norm minimization problem (10) is given in terms of the estimates of ϕ.(2)
Due to (18) and (19), calculate μ and A⌢(φ) so that RD is realized.(3)
Calculate the gradientΓ and Hessian matrix HU⌢∗U⌢ in terms of (29) and (30). Then, update U⌢(θ) in terms of (31).(4)
If the stopping criterion (32) is satisfied, stop the update and go to the next step. Otherwise, return to step (3).(5)
Estimateθ by searching the spatial spectrum and then φ is estimated in terms of (19).(6)
Estimate the azimuth and elevation angles according to the estimates of electric angles.


### 3.2. Discussion

Note that electric angleϕ is estimated based on R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))) (noise subspace corresponding to ϕ) and the estimates of electric angles θ and φ obtained by R(P⌢) (noise subspace corresponding to θ and φ) rely on the estimate quality of ϕ. Moreover, the elevation angle is estimated based on the estimate of the azimuth angle and both of them use the arc-tangent operator. It is worth pointing out that since the arc-tangent is a one-to-one function and its argument is allowed to be larger than unity, the proposed method never fails. This is the major reason why the azimuth and elevation angles are estimated using three electric angles instead of two electric angles.Regarding the proposed array configuration, we make full use of structural feature for 2D DOA estimation, which shows no failure in the entire ranges ofα and β. Furthermore, when seen from above, the proposed array configuration seems to be an L-array configuration. Referring to [10], the Cramer-Rao bound (CRB) based on the L-shaped array configuration is lower than those of other array configurations. Therefore, the proposed array configuration not only guarantees that the failure probability tends to zero, but it can also have high precision.


## 3.1. CS-RD-lp-MUSIC

Step 1 (estimateϕ by CS-MUSIC).
CS-MUSIC, which is an extension of MUSIC, can identify the parts of support using CS-based methods, after which the remaining supports are estimated by the generalized MUSIC criterion. The main contribution of CS-MUSIC is to overcome the error caused by losing rank which may cause disastrous consequences in conventional MUSIC.Letsupp⁡G-={1≤g≤Pϕ:G-g≠0} and R(A⌢(ϕ)) denote the support of G- and the range space of A⌢(ϕ), respectively. Due to (4), the number of snapshots is extended to 3T which is one of the advantages of the proposed method. It is obvious that CS-MUSIC can be simplified as MUSIC for K≤3T. However, CS-MUSIC can estimate DOA with success but MUSIC fails in the K>3T case. In CS-MUSIC, if K>3T, K-3T indices of supp⁡G- are determined by CS-based methods such as simultaneous orthogonal matching pursuit (SOPM) [21] so that column vector space R(A⌢ηK-3T(ϕ)) is decided, where ηK-3T is the set of indices and A⌢ηK-3T(ϕ) is the submatrix of A⌢(ϕ) with columns indexed by ηK-3T. The remaining 3T indices of supp⁡G- are obtained by the generalized MUSIC criterion. To make CS-MUSIC applicable for all range of the snapshot, two new orthogonal spaces R(A⌢G-) and R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))) are constructed where PR(Q) is the orthogonal projection onto the noise subspace R(Q). Then, the electric angle ϕ can be estimated by the spectrum search. Now, the major steps of CS-MUSIC are summarized as follows.(1)
FindK-3T indices of supp⁡G- by SOMP and let ηK-3T be the set of indices.(2)
DetermineA⌢ηK-3T(ϕ) in terms of ηK-3T and construct the noise subspace of CS-MUSIC R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))).(3)
Fori∈{1,2,…,Pϕ}, calculate the spatial spectrum, Pspectral=1/A⌢iH(ϕ)[PR(Q)-PR(QQHA⌢ηK-3T(ϕ))]A⌢i(ϕ), and then estimate ϕ by the locations of K highest peaks of the spatial spectrum.Step 2 (estimateθ and φ by RD-lp-MUSIC).
In this step, RD-lp-MUSIC is exploited to estimate electric angles θ and φ in terms of the estimates of ϕ. Let {ϕ-1,ϕ-2,…,ϕ-K} denote the locations of K highest peaks. Then, the row number corresponding to locations of K highest peaks is denoted as {i1,i2,…,iK} and these rows contain information about (θ,φ). It is clear that the ith row vector of source matrix [g⌢1,1(i),g⌢2,1(i),…,g⌢T,1(i),…,g⌢1,3(i),g⌢2,3(i),…,g⌢T,3(i)],  i=i1,i2,…,iK, is a 1×3T vector and thus the matrix form G⌢(i) is a T×3 matrix whose element in the (m,n)th entry is g⌢m,n(i). Denote A~(θ,φ)=A^T(θ,φ) and S⌢=[s(1),s(2),…,s(T)]T so that the lp-norm minimization problem can be expressed as(10)Jpθ,φ=minθ,φ⁡JpG⌢i-S⌢A~θ,φ=minθ,φ⁡∑t=1TG⌢ti-S⌢tA~θ,φ2pi=i1,i2,…,iK.Note thatK  lp-norm minimization problems have the same structure and can be solved in a similar manner. Hence, we only need to find the optimal solution of one minimization problem. The remaining K-1 minimization problems can be handled in the same way.Since the objective function (10) requires an exhaustive 2D search, high computational cost is needed for precise estimation which can result in the reduction of algorithmic efficiency. To avoid heavy computational load, RD-lp-MUSIC is proposed for 2D estimation just through 1D search. A detailed derivation process of exploiting RD-lp-MUSIC for estimating θ and φ is given as follows.DenoteF(θ,φ)=G⌢t(i)-S⌢tA~(θ,φ)22 so that RD can be realized by minimizing F(θ,φ). By the property of ·2, we have(11)Fθ,φ=G⌢ti-S⌢tA~θ,φG⌢ti-S⌢tA~θ,φH=const+S⌢tA⌣θ⊙A⌢φA⌣θ⊙A⌢φHS⌢tH-2Re⁡S⌢tA⌣θ⊙A⌢φG⌢tiH,where A⌣(θ)=[a˘(θ1),a˘(θ2),…,a˘(θK)]T and A⌢(φ)=[a⌢(φ1),a⌢(φ2),…,a⌢(φK)]T. Denote that P(K)=[E11(K),E22(K),…,EKK(K)]T and P(3)=[E11(3),E22(3),E33(3)]T are K2×K and 9×3 selection matrices, respectively, so that the following equation holds based on selection matrices [22]:(12)A⌣θ⊙A⌢φ=PKTA⌣θ⊗A⌢φP3=PKTA⌣θ⊗IK×KI3×3⊗A⌢φP3.During the aforesaid process, selection matrices are of the greatest importance. The purpose of utilizing selection matrices is to separateθ from φ and reduce dimension. Moreover, based on the following equation,(13)Q1I3×3⊗A⌢φQ2=A⌢φ,where Q1=[IK×K0K×2K] and Q2=[I3×303×6]T, (12) can be further rewritten as (14)A⌣θ⊙A⌢φ=PKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3.Subsequently, by substituting (14) into (11), we have(15)Fθ,φ=const-2Re⁡S⌢tPKTA⌣θ⊗IK×K·Q1-1A⌢φQ2-1P3G⌢tiH+S⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3P3HQ2-1HA⌢Hφ·Q1-1HA⌣Hθ⊗IK×KPKTHS⌢tH.To eliminate trivial solutionsA⌣(θ)=0K×3 and A⌢(φ)=0K×3, a constraint condition eKH(1)[A⌣(θ)⊙A⌢(φ)]e3(2)-1=0 is considered. Hence, the cost function L(θ,φ) is given as the following form:(16)Lθ,φ=S⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φ·Q2-1P3P3HQ2-1HA⌢HφQ1-1HA⌣Hθ⊗IK×KPKTHS⌢tH-2ReS⌢tPKTA⌣θ⊗IK×KQ1-1A⌢φ·Q2-1P3G⌢tiH-μeKH1·PKTA⌣θ⊗IK×KQ1-1A⌢φQ2-1P3e32-1,where μ is a constant. By setting the partial derivation of (16) with respect to A⌢(φ) to zero, we have(17)∂Lθ,φ∂A⌢φ=2WS⌢tHS⌢tWA⌢φVVH-2Re⁡WHS⌢tHG⌢tiHVH-μWHeK1e3H2VH=0,where W=P(K)T[A⌣(θ)⊗IK×K]Q1-1 and V=Q2-1P3. By (17) and the constraint condition eKH(1)[A⌣(θ)⊙A⌢(φ)]e3(2)=1, μ can be expressed as(18)μ=2-2eKH1S⌢t-1Re⁡G⌢tie32eKH1S⌢t-1S⌢t-1HeK1.Therefore, it can be deduced from (17) and (18) that (19)A⌢φ=W-1S⌢t-1ReG⌢ti+1-eKH1G⌢ti-1Re⁡G⌢tie32eK1e3H2eKH1S⌢t-1S⌢t-1HeK1·V-1.As one may note,A~(θ,φ) is transformed into U⌢(θ) by inserting A⌢(φ) into A~(θ,φ) and thus RD is realized so that 2D DOA estimation just requires 1D search. Furthermore, RD can avoid the failure occurrence in separating θ from φ. As a matter of fact, RD also provides an important clue for expounding the relationship between θ and φ. Therefore, RD is the foundation of estimating electric angles θ and φ by RD-lp-MUSIC. Then, we focus on the solver for the last remaining problem to estimate θ; that is,(20)Jpθ=minθ⁡∑t=1TG⌢ti-S⌢tU⌢θ2p.By denotingrt=G⌢t(i)-S⌢tU⌢(θ) as the residual vector, the objective function f(rit) can be expressed as (21)frit=∑t=1T∑i=13rit2p/2.Based on the following equation,(22)∂f∂rit∗=12∂f∂Re⁡rit+j∂f∂Im⁡rit,the partial derivative with respect to (rit)∗ can be given by(23)∂f∂rit∗=p2∑t=1Trit2p/2-1rit.Due to (23), the gradient of f with respect to rt is given by (24)∂f∂rt=∂f∂r1t∗∂f∂r2t∗∂f∂r3t∗=p2rt⊙Λ,where(25)Λ=∑t=1Tr1t2p/2-1∑t=1Tr2t2p/2-1∑t=1Tr3t2p/2-1.Denote F=diag⁡(Λ) so that (24) can be rewritten as(26)∂f∂rt=p2rtF.Note thatF is a positive definite matrix clearly. Then, since the second-order partial derivative of f with respect to rit is given as the following form,(27)∂2f∂ritT∂rjt∗=p2∑t=1Trit2p/2-1p2-1∑i=13rit2-1rit2+1i=j0i≠j,the 3×3 Hessian matrix of f with respect to rt is given by(28)Hrtrt∗=∂2f∂rtT∂rt∗=p2·Fp2-1∑i=13rit2-1diag⁡rtTrt+I3×3.Denote the gradient off with respect to U⌢(θ) as Γ=[Γ1,Γ2,Γ3] so that Γi is computed as(29)Γi=∂f∂U⌢i∗=-S⌢tT⊗e3Tip2rtFT.Then, the Hessian matrix off with respect to U⌢(θ) is (30)HU⌢∗U⌢=∂2f∂U⌢∗∂U⌢T=S⌢tT⊗ΣP3Hrtrt∗P3TS⌢t∗⊗ΣT,where Σ=[e3(1)e3(2)e3(3)]. By applying the Newton method, the following sequence U⌢(z)(θ), z=0,1,2,…, is obtained to find the minimization of f:(31)U⌢z+1θ=U⌢zθ+σzΔzU⌢zθ,where σz is a positive step and Δ(z)(U⌢(z)(θ))=-HU⌢∗U⌢Γ(U⌢(z)(θ)). The initial value of U⌢(z)(θ) is defined as U⌢(0)(θ)=[(S⌢t)HS⌢t]-1(S⌢t)HG⌢t(i). Specifically, the iteration is terminated if the following stopping criterion is satisfied:(32)Jpθz-Jpθz+1Jpθz<εfor some small ε. In the simulations, ε is set to 10-7. After determining U⌢(θ), the projection matrix onto the noise subspace is given by(33)P⌢=I3×3-UθUHθUθ-1UHθ,where U(θ)=U⌢H(θ). Then, by computing the spatial spectrum, Pspetral(θ)=1/a˘H(θ)P⌢a˘(θ), and searching its peaks, we can provide an accurate estimation for θ. Once θ is obtained from the spatial spectrum, φ is easily estimated in terms of (19). Therefore, RD-lp-MUSIC avoids pair matching.After electric angles are estimated, the azimuth angleαk and elevation angle βk of the kth source can be given by αk=arctan⁡φk/θk and βk=arctan⁡ϕksin⁡αk/φk, respectively. The major steps of RD-lp-MUSIC for estimating θ and φ are given as follows.(1)
Thelp-norm minimization problem (10) is given in terms of the estimates of ϕ.(2)
Due to (18) and (19), calculate μ and A⌢(φ) so that RD is realized.(3)
Calculate the gradientΓ and Hessian matrix HU⌢∗U⌢ in terms of (29) and (30). Then, update U⌢(θ) in terms of (31).(4)
If the stopping criterion (32) is satisfied, stop the update and go to the next step. Otherwise, return to step (3).(5)
Estimateθ by searching the spatial spectrum and then φ is estimated in terms of (19).(6)
Estimate the azimuth and elevation angles according to the estimates of electric angles.


## 3.2. Discussion

Note that electric angleϕ is estimated based on R(PR(Q)-PR(QQHA⌢ηK-3T(ϕ))) (noise subspace corresponding to ϕ) and the estimates of electric angles θ and φ obtained by R(P⌢) (noise subspace corresponding to θ and φ) rely on the estimate quality of ϕ. Moreover, the elevation angle is estimated based on the estimate of the azimuth angle and both of them use the arc-tangent operator. It is worth pointing out that since the arc-tangent is a one-to-one function and its argument is allowed to be larger than unity, the proposed method never fails. This is the major reason why the azimuth and elevation angles are estimated using three electric angles instead of two electric angles.Regarding the proposed array configuration, we make full use of structural feature for 2D DOA estimation, which shows no failure in the entire ranges ofα and β. Furthermore, when seen from above, the proposed array configuration seems to be an L-array configuration. Referring to [10], the Cramer-Rao bound (CRB) based on the L-shaped array configuration is lower than those of other array configurations. Therefore, the proposed array configuration not only guarantees that the failure probability tends to zero, but it can also have high precision.


## 4. Simulation Results

In this section, several simulation results are presented to validate the superior performance of the proposed method as compared to PM of one L-shaped array configuration and 2D-MUSIC. The proposed method and 2D-MUSIC are performed in the proposed array configuration and PM is performed in the L-shaped array configuration of thex-y plane. The total numbers of elements for array configurations and the spacing between the adjacent elements are set to 9 and d=λ/2, respectively. Although the array configurations for three methods are different, the degrees of array configurations are the same and the projections of the array configurations on the x-y plane are L-shaped array configurations. In the typical DOA estimation, since the sources always come from above the sensor, that is, above the x-y plane, the azimuth angle is in the range of 0° and 360° and the elevation angle is in the range of 0° and 90°. Hence, electric angles θ and φ belong to [-180°,180°] and ϕ belongs to [-180°,0°]. All simulation results are obtained from 100 Monte Carlo runs. In the simulations, the root mean squared error (RMSE) of 2D DOA estimation is defined as(34)RMSE=1K∑k=1KRMSEαk+RMSEβk,where RMSE(αk)=∑l=1100(α~k,l-αk)2/1001/2, RMSE(βk)=∑l=1100(β~k,l-βk)2/1001/2, and α~k,l and β~k,l are the estimates of αk and βk in the lth run, respectively. It is well known that the grid interval serves as a tradeoff between precision and computational complexity. For all simulations, we make the coarse grid interval with 1° and perform a local fine grid interval in the locations obtained by utilizing the coarse grid interval.In the first simulation, we show angle estimation results of three methods in the azimuth-elevation plane. Consider three independent sources with DOAs of(α1=3°,β1=2°), (α2=45°,β2=40°), and (α3=80°,β3=85°) impinging on the proposed array configuration. Figure 2 presents angle estimation results of three methods for all three sources with the fixed SNR 5 dB and number of snapshots being 100. It is indicated in Figure 2 that CS-RD-lp-MUSIC and 2D-MUSIC can provide correct estimates for the azimuth and elevation angles of three sources but PM fails. Moreover, although CS-RD-lp-MUSIC slightly outperforms 2D-MUSIC in terms of angle estimation results, it avoids an exhaustive 2D search which is needed in 2D-MUSIC.Figure 2
Angle estimation results of three methods for three sources.TheRMSE of three methods versus SNR and the number of snapshots is investigated in the second simulation. We keep the same source model as in the first simulation. Figure 3 depicts RMSE as a function of SNR of three methods and CRB [23] with the fixed number of snapshots being 100, whereas RMSE versus the number of snapshots with the fixed SNR 5 dB is shown in Figure 4. It can be concluded from Figures 3 and 4 that CS-RD-lp-MUSIC has more precise estimation than 2D-MUSIC and PM with no estimation failure. It is clearly seen that the performance of CS-RD-lp-MUSIC is gradually improving and is close to the CRB with the increase of SNR and the number of snapshots.Figure 3
RMSE versus SNR with the fixed number of snapshots being 100.Figure 4
RMSE versus number of snapshots with the fixed SNR 5 dB.Figure5 illustrates the relation between the bias of DOA estimation and the angle separation of two independent sources. The bias is defined as the difference between the estimated angle and the real angle, which can indicate the degree of deviation from the true angle. The smaller the bias is, the better performance the method has. Therefore, the bias is the significant performance index. Consider two sources impinging from DOAs of (α1=25°,β1=20°) and (α2=α1+δ,β2=β1+δ), where the step of the angle separation δ is 1°. The SNR is 3 dB and the number of snapshots is 50. As can be seen from Figure 5, there is the bias for small angle separation using three methods and the bias of three methods disappears as long as the angle separation is no less than 17°.Figure 5
Bias versus angle separation with the fixed SNR 3 dB and number of snapshots being 50.In the fourth simulation, we compare the performance of three methods for coherent sources by showing theRMSE versus SNR and the number of snapshots. Consider two coherent sources impinging from DOAs of (α1=30°,β1=45∘) and (α2=70°,β2=75°). Since the conventional 2D-MUSIC method is incapable of handing the coherent sources, the forward spatial smoothing method is exploited in the 2D-MUSIC called 2D-FSS-MUSIC to estimate the coherent sources. Figures 6 and 7 plot the RMSE versus SNR and the number of snapshots for coherent sources, respectively. It can be seen from Figures 6 and 7 that the proposed method has the best estimation accuracy among all three methods for coherent sources. Moreover, this performance advantage is gradually improving with SNR or the number of snapshots increasing.Figure 6
RMSE versus SNR with the fixed number of snapshots being 100 for coherent sources.Figure 7
RMSE versus number of snapshots with the fixed SNR 5 dB for coherent sources.Finally, Figure8 presents the RMSE of the proposed method for all possible azimuth and elevation angles with the fixed SNR 3 dB and number of snapshots being 50. One single source is considered in this simulation. The steps of the azimuth and elevation angles are both fixed at 5°. We observe from Figure 8 that no estimation failure occurs for all pair angles with the proposed method.Figure 8
RMSE of the proposed method at different DOAs with the fixed SNR 3 dB and number of snapshots being 50 for one single source.


## 5. Conclusion

In this paper, a novel CS-RD-lp-MUSIC is proposed for 2D DOA estimation in CS. The proposed method introduces electric angles and then separates the steering vector into two parts for constructing two corresponding noise subspaces. The electric angles are estimated by CS-MUSIC and RD-lp-MUSIC based on the constructed noise subspaces, so that the azimuth and elevation angles are obtained by arc-tangent operations in terms of estimates of electric angles. Since the arc-tangent is a one-to-one function and allows the value of its argument to be larger than unity, the proposed method never fails. The proposed method, which requires no pair matching, can reduce computational complexity and extend the number of snapshots to improve performance. Simulation results show that the proposed method never fails for all pair angles and has better estimation performance than PM and 2D-MUSIC in terms of RMSE and bias of DOA estimation.



---

*Source: 792181-2015-06-15.xml*