# A 2D Nested Array Based DOA Estimator for Incoherently Distributed Sources via Sparse Representation Utilizing L1-Norm

Authors: Tao Wu; Yiwen Li; Zhengxin Li; Yijie Huang; Jiwei Xu
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

Nested arrays are sparse arrays composed of subarrays with nonuniform sensor spacing. Compared with traditional uniform arrays, nested arrays have more degree of freedoms (DOFs) and larger apertures. In this paper, a nested array has been proposed as well as a direction-of-arrival (DOA) estimation method for two-dimensional (2D) incoherently distributed (ID) sources. A virtual array is firstly obtained through vectorization of the cross-correlation matrix of subarrays. Sensor positions of the virtual array and the optimal configuration of the nested array are derived next. Then rotational invariance relationship for generalized steering matrix of the virtual array with respect to nominal azimuth is deduced. According to the rotational invariance relationship, sparse representation model underl1-norm constraint is established, which is resolved by transferring the objective function to second-order cone constraints and combining a estimation residual error constraint for receive vector of the virtual array. Simulations are conducted to investigate the effectiveness of the proposed method in underdetermined situation and examine different experiment factors includingSNR, snapshots, and angular spreads as well as sensor number of subarrays. Results show that the proposed method has better performance than uniform parallel arrays with the same number of sensors.


---

## Body

## 1. Introduction

In the field of array signal processing, traditional DOA estimation of targets is based on point source models. Point source models suppose signal propagation from a target to receive array is one straight path; geometry of the target can be ignored, which can be considered as a point. In the surrounding of underwater acoustics detection, on one hand, there exist reflection paths from seabed and sea surface between a target and receive array; on the other hand, when the distance from the target to receive array is short, as many parts of the target reflect signal and spatial distribution and the geometry of the target cannot be ignored; signal propagations from the target to receive array cannot be supposed to a single path. Point source models cannot describe spatial properties of sources under such circumstance; estimations based on point source models are no longer effective, which should be replaced by distributed sources [1]. Distributed sources can be regarded as an assembly of point sources which can be called scatterers within a spatial distribution.A distributed source may be coherently distributed (CD) source or incoherently distributed (ID) source according to the coherence of scatterers of a target. A coherently distributed (CD) source means that scatterers of the target are coherent. Spatial distribution of scatterers of a CD source is described by deterministic angular signal distribution function (ASDF). On the contrary, an incoherently distributed (ID) source means that scatterers of the target are uncorrelated. Spatial distribution of scatterers of an ID source is described by angular power density function (APDF). Both ASDF and APDF can be generally modeled as Gaussian or uniform distribution according to the geometry and spatial distribution of a target. In this paper, ID sources are considered.Distributed sources may be one-dimensional (1D) or two-dimensional (2D) which depends on whether impinging signals of scatterers and arrays are in the same plane. As to a 1D source, parameters of ASDF or APDF include nominal angle and angular spread. Nominal angle can be called nominal DOA representing the center of a target while angular spread represents spatial extension of the target. As to a 2D source, parameters contain nominal azimuth, nominal elevation, azimuth spread, and elevation spread. Nominal azimuth and nominal elevation represent center of a target which also can be expressed as DOA; azimuth spread and elevation spread indicating 2D spatial extension of the target are collectively called angular spreads. As impinging signals of scatterers and arrays are not in the same plane practically, 2D sources are more general and reasonable.As to CD sources, utilizing different array configurations representative estimators have been proposed in [2–10], which are mostly based on rotational invariance relationship derived by Taylor approximation of generalized steering vectors. For ID sources, several algorithms have been presented and can be roughly divided into four categories: subspace-based algorithms developed from multiple signal classification (MUSIC) such as distributed signal parameter estimator (DSPE) [1] and dispersed signal parametric estimation (DSPARE) [11], generalizations of Capon’s methods [12–14] involving high-order matrix inversions and spectral searches, the maximum likelihood (ML) approaches [15–17] requiring high computational complexity but with better accuracy, and covariance matching estimation techniques (COMET) [18–21] with lower computational complexity than ML approach but the same large sample behavior.The aforementioned estimations of ID sources suppose the sources and receive arrays are in the same plane, which are special cases modeled as 1D ID sources. Involving more parameters, there have been relatively few studies on estimation of 2D ID sources. An extension of COMET has been proposed in [22] which employs alternating projection principle separating source powers and noise variances from objective functions and then estimates parameters by a four-dimensional nonlinear optimization. Utilizing double parallel uniform linear arrays, a TLS-ESPRIT-like approach has been proposed for DOA estimation in [23], where nominal elevations are firstly estimated via rotational invariance relationship derived by the first-order Taylor series expansions of the generalized steering vectors; then nominal azimuths are obtained through 1D searching. An estimator based on uniform rectangular arrays has been proposed in [24], where DOAs are estimated through an ESPRIT like algorithm by virtue of rotational invariance relationship of receive vectors which are described by generalized signal vector and generalized steering vector composed of nominal steering vector and its first-order partial derivatives.Most distributed sources estimation techniques mentioned above are based on uniform arrays such as uniform linear arrays, uniform rectangular arrays, or circular arrays, which require the distances between sensors less than or equal to half of impinging signal wave. For achieving better accuracy and resolved more sources, estimation should be performed with larger aperture, which needs to increase number of sensors. Consequently, this would lead to an increase in complexity and cost of sonar system. The recent proposed nested arrays [25] possess larger apertures and more degree of freedoms (DOFs) compared with uniform arrays with the same number of sensors. DOA estimations for point sources under different kinds of nested arrays have been proposed in [26–30]. As for distributed sources, based on a linear nested array, a spatial smoothing with spectral search technique has been proposed in [8] for 1D distributed sources supposing a priori knowledge of the angular spreads is known, which is impractical in practice. Using the same array configuration, the authors of [21] have proposed annihilating filter technique combined with structured low rank approximation for 1D distributed sources without a priori knowledge of the angular spreads, where the aperture of virtual array has not been fully utilized as half of the sensors for estimation of DOA and the other half for angular spreads.In this paper, DOA estimation for 2D ID sources is proposed based on a 2D nested array composed of double parallel uniform linear subarrays. According to the concept of difference coarray, the cross-correlation matrix of subarrays is firstly vectorized by Khatri-Rao product to obtain a hole-free virtual array; and the closed form expression for sensor positions of the virtual array is presented. Then rotational invariance relationship with regard to nominal azimuth is derived under the assumption of small angular spreads. Sparse representation withl1-norm constraint for DOA estimation of 2D ID sources is established next. Objection function is then converted to second-order cone constraints to avoid determining the regularization parameters. Based on estimation residual error for receive vector of virtual array, a robust constraint is derived. Followed by the fact that nominal azimuths are resolved by the sparse representation framework, paired nominal elevations are obtained by least-square techniques. The proposed method can detect 2D ID sources more than sensors without angles matching and without information of APDF of sources. Simulations are conducted to investigate the effectiveness and superior of the proposed method lastly.


## 2. The Nested Array and Source Model

Generally, nested arrays are composed of several levels of uniform linear subarrays; the distances between sensors within the first level subarray are less than or equal to half of impinging signal wave; distances between sensors within the rest subarrays are far greater than the first level subarray. As shown in Figure1, the proposed nested array is composed of double parallel uniform subarrays onxoz plane X1 and X2. Subarray X1 is parallel tox axis and consists of the middle sensor placed onz axis andM1 sensors located on both sides. The sensor number of subarray X1 isN1=2M1+1; the distance between sensors of subarray X1 isd. Subarray X2 is onx axis withM2 sensors located on both positive and negative semiaxes. The sensor number of subarray X2 isN2=2M2; the distance between sensors of subarray X2 isD=N1d. The distance between subarrays X1 and X2 is alsod.Suppose that there areq 2D ID narrowband sources with nominal angles (θi, φi) (i=1,2,⋯,q) impinging into the nested array. θi denotes a spatial nominal azimuth betweenx axis and propagation path of theith source; φi is nominal elevation of theith source. θi∈[0,π], φi∈0,π. λ is the wavelength of the impinging signal.Figure 1
(a) Illustration of the nested array configuration in three dimensions; (b) Illustration of the nested array configuration onxoz plane.
(a) (b)Sensor positions set of subarrays X1 and X2 can be expressed, respectively, as(1)L1=lmd,0,d∣lm=-M1,-M1-1,⋯,M1-1,M1,(2)L2=ln+0.5D,0,0∣ln=-M2,-M2-1,⋯,M2-1.N 1×1 dimensional receive vector of subarrays X1 andN2×1 dimensional receive vector of subarray X2 can be expressed as(3)x1t=∑i=1q∬a1θ,φsiθ,φ,tdθdφ+nx1t,(4)x2t=∑i=1q∬a2θ,φsiθ,φ,tdθdφ+nx2t,where nx1(t) and ny1(t) are additive noise vectors which are assumed to be uncorrelated between sensors and Gaussian white with zero mean.a1(θ,φ) anda2(θ,φ) denote the steering vector of subarrays X1 and X2 with respect to point source.si(θ,φ,t) is the complex random angular signal density of theith distributed source representing the scatterer intensity of the source from angle (θ,φ) at the snapshot indext. Unlike a point source, signal of a distributed source exists not only in the single direction (θi, φi) but also in a spatial distribution around (θi,φi). An ID source means different scatterers from one target generate uncorrelated multipath signals. Therefore, one direction fromsi(θ,φ,t) is uncorrelated with other directions fromsi(θ,φ,t); we have(5)Esiθ,φ,tsi∗θ′,φ′,t=σi2fθ,φ;uiδθ-θ′δφ-φ′,where σi2 is power of theith source,δ(•) is the Kronecker delta function, (•)∗ means complex conjugate, andf(θ,φ;u) is the source’s normalized APDF which reflects the distribution of scatterers of a source and is related to the geometry and surface property of the source. APDF can be modeled as 2D Gaussian, uniform, or any other distribution function according to different characteristic of an ID source, which is determined by parameter setui=[θi,ϕi,σθi,σϕi] denoting nominal azimuth, nominal elevation, azimuth spread, elevation spread, respectively, satisfying(6)∬fθ,φ;uidθdφ=1.For a Gaussian ID source, APDF can be expressed as(7)fθ,φ;ui=12πσθiσφiexp⁡-12σθiσφiθ-θiσθi2+φ-φiσφi2.For a uniform ID source, APDF can be expressed as(8)fθ,φ;ui=14σθiσφiθ-θi≤σθiandφ-φi≤σφi0θ-θi≥σθiorφ-φi≥σφi.


## 3. Proposed Method

This section consists of five parts. First of all, a virtual array is obtained by difference of sensor positions in two subarrays. Then, optimal configuration of the nested array is analyzed and the rotational invariance relationship of generalized steering matrix of the virtual array is derived. In the next part, DOA estimation via sparse representation is established and solution process is detailed. Lastly computational procedure is summarized and complexity analysis is analyzed with comparison of two existing methods for 2D ID sources.


### 3.1. Description of Virtual Array

The main advantage of nested array lies in increasing DOFs by forming virtual array through subarrays. According to the assumption of ID sources, different sources are uncorrelated. TheN1×N2 dimensional cross-correlation matrix of subarrays X1 and X2 can be expressed as follows:(9)R=Ex1tx2Ht=∑i=1qσi2∬a1θ,φfθ,φ;uia2Hθ,φdθdφ.VectorizingR by Khatri-Rao product, we obtain(10)r-=recR=∑i=1qσi2∬a2∗θ,φ⊙a1θ,φfθ,φ;uidθdφ,whererec(•) denotes vectorization of a matrix, ⊙ denotes Khatri-Rao product, and (•)H denotes Hermitian transpose. A virtual array can be obtained by vectorization ofR which is shown in Figure 2. Sensor positions of the virtual array are determined by difference of sensor positions in subarrays X1 and X2. The set of sensor positions of the virtual array can be expressed as(11)L=pi-pj∣pi∈L1,pj∈L2.Figure 2
Virtual array formed by the nested array.From (1) and (2) we obtain sensor position set of the virtual array L as(12)L=n-0.5d,0,d∣n=1-2M1M2+M2,2-2M1M2+M2,⋯,2M1M2+M2,which indicates that the virtual array has 2(2M1M2+M2) virtual sensors without repetitive positions.It is noteworthy that elements inr- are not arranged in accordance with sensor order in virtual array. Supposer is receive signal vector from left side of virtual sensor located on [(0.5-2M1M2-M2)d, 0,d] to the right side of sensor located on [(2M1M2+M2-0.5)d, 0,d]. Define a block diagonal matrixJ for sorting order of r-, which is composed ofN2 number of blocks JN1 which areN1 dimensional inverse identity matrix.J can be expressed as follows:(13)J=JN10⋱0JN1.Then,r can be written as(14)r=Jr-=∑i=1qσi2∬bθ,φfθ,φ;uidθdφ,whereb(θ,φ) denotes the steering vector of subarray virtual array with respect to point source, which can be written as(15)bθ,φ=ej2πdcos⁡φ/λe-j2π2M1M2+M2-0.5dcos⁡θ/λ,e-j2π2M1M2+M2-1.5dcos⁡θ/λ,⋯,ej2π2M1M2+M2-0.5dcos⁡θ/λT.


### 3.2. Optimal Configuration of the Nested Array

From (2), DOF of the virtual array is 2(2M1M2+M2); sensor number of the nested array isN=2(M1+M2)+1. It can be concluded that DOF of the proposed nested array can achieve O(N2) which is far more than number of sensors. We investigated the most DOF can be achieved under the given sensor numberN; the optimization problem can be expressed as follows:(16)max22M1M2+M2subjecttoM1+M2=N-12.Equation (16) can be solved through minimization of Lagrange function as follows:(17)LM1,M2=22M1M2+M2-λ2M1+2M2+1.Taking the derivative of (17) with regard toM1,M2 setting the result equal to 0, we obtain (18)M1opt=N4-12M2opt=N4.AsM1,M2, and (N-1)/2 must be integers, optimal solutions are all fractions and need to be modified, which are listed in Table 1.Table 1
The optimal configuration and DOF of the proposed nested arrays.
N(h is a integer) M 1 M 2 DOF 4h +1 h h 2h2+2h 4h +3 h h +1 2h2+4h+2


### 3.3. Rotational Invariance Relationship of Virtual Array

In the point source models case, steering vectors are used to describe response of an array. Nevertheless, as for distribute sources, generalized steering vectors or generalized steering matrices [5, 23] represent response of arrays. Define βi as the generalized steering vector of the virtual array with respect to theith distributed source, which can be expressed as(19)βi=∬bθ,φfθ,φ;uidθdφ,where thekth element of βi means the response ofkth sensor to theith distributed source. Generalized steering matrix of the virtual array can be expressed as(20)B=β1,β2,⋯,βq.r can be expressed as (21)r=Bσ,where σ is the power vector with following expression:(22)σ=σ12,σ22,⋯,σq2T.Suppose the angular spreads are small andd/λ is set at 1/2; we have the following relationship (see Appendix A)(23)βik+1≈ejπcos⁡θiβik,where [•]k represent thekth element of vector. So we can obtain the following rotational invariance relationship:(24)Bk+1≈BkΦ,where Φ=diag[ejπcos⁡θ1,ejπcos⁡θ2,⋯,ejπcos⁡θq] is the rotation matrix,Bk is thekth row ofB representing the response ofkth sensor with respect to distributed sources.


### 3.4. DOA Estimation via Sparse Representation

Representing the response of first sensor to distributed sources,B1 is the first row ofB. Change variables (θi+θ~) for θ and (φi+φ~) for φ, where θ~ and φ~ are the small deviations of θi and φi.B1 can be regarded as Hadamard product of three vectors which are nominal azimuth coefficient vector, nominal elevation coefficient vector α(φ), and real coefficient vectorg.B1 can be expressed as (see Appendix B)(25)B1≈e-jπ2M1M2+M2-0.5cos⁡θ1,e-jπ2M1M2+M2-0.5cos⁡θ2,⋯,e-jπ2M1M2+M2-0.5cos⁡θq⊕αφ⊕g,where ⊕ represent Hadamard product. α(φ) andg can be written as (26)αφ=ejπcos⁡φ1,ejπcos⁡φ2,⋯,ejπcos⁡φq,(27)g=g1,g2,⋯,gq,where gi=∬ejπ(-sin⁡φiφ~)e-jπ(2M1M2+M2-0.5)(-sin⁡θiθ~)f(θi+θ~,φi+φ~;ui)dθ~dφ~. According to the rotational invariance relationship (24), thekth row of B can be obtained as(28)Bk≈e-jπ2M1M2+M2-k+0.5cos⁡θ1,e-jπ2M1M2+M2-k+0.5cos⁡θ2,⋯,e-jπ2M1M2+M2-k+0.5cos⁡θq⊕αφ⊕g.Define a nominal azimuth coefficient matrix as(29)Aθ=αθ1,αθ2,⋯,αθq,where(30)αθi=e-jπ2M1M2+M2-0.5cos⁡θi,e-jπ2M1M2+M2-1.5cos⁡θi,⋯,ejπ2M1M2+M2-0.5cos⁡θiT.Thus, receive vector of virtual array can be expressed as(31)r=AθP,where(32)P=αφT⊕gT⊕σ,Divide the nominal azimuth space intoQ partitions. Suppose set Θ={θ~1,θ~2,⋯,θ~Q} contains all the possible nominal azimuths. Sparse representation of receive vectorr can be expressed as(33)r=Aθ~P~,where P~ is sparse coefficient vector. Thejth element of vector P~ equals theith element ofP if θ~j=θi, (j=1,⋯,Q;i=1,⋯,q). Therefore, θi can be obtained from the positions of nonzero elements of θ~j. Thoughl0 norm is a direct measurement of sparsity of a vector,l0 norm constraint optimization is a NP-hard problem; we usel1 norm constraint to replacel0 norm constraint. Thus, DOA estimation via sparse representation underl1 norm constraint can be expressed as(34)minP~1subjecttor=Aθ~P~.The cross-correlation matrixR can be replaced by sample cross-correlation of subarrays withL snapshots as follows:(35)R^=1L∑t=1Lx1tx2Ht.Then we have the vectorization of receive signal vector of virtual array as follows:(36)r^=JrecR^=r+Δr,where Δr is estimation residual error of virtual array, which has the following probability distribution (see Appendix C) (37)Δr~AsN0,1LJR2T⊗R1JH,whereAsN(μ,C) represents multidimensional normal distribution with mean value μ and covariance matrixC. ⊗ is Kronecker product.R1 andR2 are covariance matrices the subarrays X1 and X2, respectively. According to the nature of normal distribution, we have (38)W-0.5Δr~AsN0,I4M1M2+2M2,(39)W-0.5Δr22~Asχ24M1M2+2M2,whereW=1/LJ(R2T⊗R1)JH,Asχ2(4M1M2+2M2) represents Chi square distribution with DOF as (4M1M2+2M2). Then we introduce a parameter η to restrain Δr as follows: (40)W-0.5Δr22≤η.Let (40) hold with a high probabilityp=0.999. η can be resolved from chi2inv(p, 4M1M2+2M2) function of Matlab. Then (34) can be converted into the following robust form(41)min1TγsubjecttoReP~i,ImP~i2≤γiW-0.5r^-Aθ~P~22≤η,where1 isQ×1 dimensional vector whose elements all are 1, (41) can be resolved by CVX.R1 andR2 can be replaced by sample covariance matrix withL snapshots.(42)R^1=1L∑t=1Lx1tx1HtR^2=1L∑t=1Lx2tx2Ht.Through tracking nonzero position of coefficientγi, θi can be resolved. Next we can obtain the estimated coefficient vectorP as P^=A(θ)†r, where (•)+ denotes pseudo inverse operation. Thus, nominal elevation of theith source can be obtained as follows:(43)φi=arccos⁡arg⁡P^iπ.


### 3.5. Computational Procedure and Complexity Analysis

Now, our algorithm can be summarized as follows.Step 1.
Compute sample cross-correlation matrixR and sample covariance matricesR1 andR2 using (35) and (42).Step 2.
VectorizeR from (10) and reorder elements of the vector to obtain r^ according to (36).Step 3.
Calculate covariance matrix of residual error vector according to (37) and determine upper bound of constraint of η.Step 4.
Discretize nominal azimuth and calculate nominal elevationθi from (41).Step 5.
Calculate nominal elevationφi from (43).We analyze the computational complexity of the proposed method in comparison with COMET [22] which can be applied for uniform 2D arrays and Zhou’s algorithm [23] using double parallel arrays for 2D ID sources. For comparison, the two methods are supposed to be applied to double parallel arrays with the same sensor number as the proposed nested array. COMET [22] is a method of four-dimensional nonlinear optimization under alternating projection algorithm framework, its computational cost is mostly composed of the calculation of the sample covariance matrix which needs O(LN2) and the alternating projection technique with respect to cost functions which is O(3N3q3+ 3N2). Zhou’s algorithm calculates nominal elevation by a TLS-ESPRIT algorithm and finds nominal azimuth by 1D searching. Computational cost of Zhou’s algorithm [23] mainly contains calculation of the sample covariance matrix O(LN2), eigendecomposition, and inversion of the sample covariance matrix O(2N3) and 1D searching O(N3). Computational cost of the proposed algorithm mainly contains calculation of the sample cross-correlation matrix and sample covariance matrix O(LN2), solving second-order core programming O[(4M1M2+2M2)3Q3] and solving the nominal elevation by pseudo inversion O(N3). With more DOF the computational cost of the proposed algorithm is higher than two methods.


## 3.1. Description of Virtual Array

The main advantage of nested array lies in increasing DOFs by forming virtual array through subarrays. According to the assumption of ID sources, different sources are uncorrelated. TheN1×N2 dimensional cross-correlation matrix of subarrays X1 and X2 can be expressed as follows:(9)R=Ex1tx2Ht=∑i=1qσi2∬a1θ,φfθ,φ;uia2Hθ,φdθdφ.VectorizingR by Khatri-Rao product, we obtain(10)r-=recR=∑i=1qσi2∬a2∗θ,φ⊙a1θ,φfθ,φ;uidθdφ,whererec(•) denotes vectorization of a matrix, ⊙ denotes Khatri-Rao product, and (•)H denotes Hermitian transpose. A virtual array can be obtained by vectorization ofR which is shown in Figure 2. Sensor positions of the virtual array are determined by difference of sensor positions in subarrays X1 and X2. The set of sensor positions of the virtual array can be expressed as(11)L=pi-pj∣pi∈L1,pj∈L2.Figure 2
Virtual array formed by the nested array.From (1) and (2) we obtain sensor position set of the virtual array L as(12)L=n-0.5d,0,d∣n=1-2M1M2+M2,2-2M1M2+M2,⋯,2M1M2+M2,which indicates that the virtual array has 2(2M1M2+M2) virtual sensors without repetitive positions.It is noteworthy that elements inr- are not arranged in accordance with sensor order in virtual array. Supposer is receive signal vector from left side of virtual sensor located on [(0.5-2M1M2-M2)d, 0,d] to the right side of sensor located on [(2M1M2+M2-0.5)d, 0,d]. Define a block diagonal matrixJ for sorting order of r-, which is composed ofN2 number of blocks JN1 which areN1 dimensional inverse identity matrix.J can be expressed as follows:(13)J=JN10⋱0JN1.Then,r can be written as(14)r=Jr-=∑i=1qσi2∬bθ,φfθ,φ;uidθdφ,whereb(θ,φ) denotes the steering vector of subarray virtual array with respect to point source, which can be written as(15)bθ,φ=ej2πdcos⁡φ/λe-j2π2M1M2+M2-0.5dcos⁡θ/λ,e-j2π2M1M2+M2-1.5dcos⁡θ/λ,⋯,ej2π2M1M2+M2-0.5dcos⁡θ/λT.


## 3.2. Optimal Configuration of the Nested Array

From (2), DOF of the virtual array is 2(2M1M2+M2); sensor number of the nested array isN=2(M1+M2)+1. It can be concluded that DOF of the proposed nested array can achieve O(N2) which is far more than number of sensors. We investigated the most DOF can be achieved under the given sensor numberN; the optimization problem can be expressed as follows:(16)max22M1M2+M2subjecttoM1+M2=N-12.Equation (16) can be solved through minimization of Lagrange function as follows:(17)LM1,M2=22M1M2+M2-λ2M1+2M2+1.Taking the derivative of (17) with regard toM1,M2 setting the result equal to 0, we obtain (18)M1opt=N4-12M2opt=N4.AsM1,M2, and (N-1)/2 must be integers, optimal solutions are all fractions and need to be modified, which are listed in Table 1.Table 1
The optimal configuration and DOF of the proposed nested arrays.
N(h is a integer) M 1 M 2 DOF 4h +1 h h 2h2+2h 4h +3 h h +1 2h2+4h+2


## 3.3. Rotational Invariance Relationship of Virtual Array

In the point source models case, steering vectors are used to describe response of an array. Nevertheless, as for distribute sources, generalized steering vectors or generalized steering matrices [5, 23] represent response of arrays. Define βi as the generalized steering vector of the virtual array with respect to theith distributed source, which can be expressed as(19)βi=∬bθ,φfθ,φ;uidθdφ,where thekth element of βi means the response ofkth sensor to theith distributed source. Generalized steering matrix of the virtual array can be expressed as(20)B=β1,β2,⋯,βq.r can be expressed as (21)r=Bσ,where σ is the power vector with following expression:(22)σ=σ12,σ22,⋯,σq2T.Suppose the angular spreads are small andd/λ is set at 1/2; we have the following relationship (see Appendix A)(23)βik+1≈ejπcos⁡θiβik,where [•]k represent thekth element of vector. So we can obtain the following rotational invariance relationship:(24)Bk+1≈BkΦ,where Φ=diag[ejπcos⁡θ1,ejπcos⁡θ2,⋯,ejπcos⁡θq] is the rotation matrix,Bk is thekth row ofB representing the response ofkth sensor with respect to distributed sources.


## 3.4. DOA Estimation via Sparse Representation

Representing the response of first sensor to distributed sources,B1 is the first row ofB. Change variables (θi+θ~) for θ and (φi+φ~) for φ, where θ~ and φ~ are the small deviations of θi and φi.B1 can be regarded as Hadamard product of three vectors which are nominal azimuth coefficient vector, nominal elevation coefficient vector α(φ), and real coefficient vectorg.B1 can be expressed as (see Appendix B)(25)B1≈e-jπ2M1M2+M2-0.5cos⁡θ1,e-jπ2M1M2+M2-0.5cos⁡θ2,⋯,e-jπ2M1M2+M2-0.5cos⁡θq⊕αφ⊕g,where ⊕ represent Hadamard product. α(φ) andg can be written as (26)αφ=ejπcos⁡φ1,ejπcos⁡φ2,⋯,ejπcos⁡φq,(27)g=g1,g2,⋯,gq,where gi=∬ejπ(-sin⁡φiφ~)e-jπ(2M1M2+M2-0.5)(-sin⁡θiθ~)f(θi+θ~,φi+φ~;ui)dθ~dφ~. According to the rotational invariance relationship (24), thekth row of B can be obtained as(28)Bk≈e-jπ2M1M2+M2-k+0.5cos⁡θ1,e-jπ2M1M2+M2-k+0.5cos⁡θ2,⋯,e-jπ2M1M2+M2-k+0.5cos⁡θq⊕αφ⊕g.Define a nominal azimuth coefficient matrix as(29)Aθ=αθ1,αθ2,⋯,αθq,where(30)αθi=e-jπ2M1M2+M2-0.5cos⁡θi,e-jπ2M1M2+M2-1.5cos⁡θi,⋯,ejπ2M1M2+M2-0.5cos⁡θiT.Thus, receive vector of virtual array can be expressed as(31)r=AθP,where(32)P=αφT⊕gT⊕σ,Divide the nominal azimuth space intoQ partitions. Suppose set Θ={θ~1,θ~2,⋯,θ~Q} contains all the possible nominal azimuths. Sparse representation of receive vectorr can be expressed as(33)r=Aθ~P~,where P~ is sparse coefficient vector. Thejth element of vector P~ equals theith element ofP if θ~j=θi, (j=1,⋯,Q;i=1,⋯,q). Therefore, θi can be obtained from the positions of nonzero elements of θ~j. Thoughl0 norm is a direct measurement of sparsity of a vector,l0 norm constraint optimization is a NP-hard problem; we usel1 norm constraint to replacel0 norm constraint. Thus, DOA estimation via sparse representation underl1 norm constraint can be expressed as(34)minP~1subjecttor=Aθ~P~.The cross-correlation matrixR can be replaced by sample cross-correlation of subarrays withL snapshots as follows:(35)R^=1L∑t=1Lx1tx2Ht.Then we have the vectorization of receive signal vector of virtual array as follows:(36)r^=JrecR^=r+Δr,where Δr is estimation residual error of virtual array, which has the following probability distribution (see Appendix C) (37)Δr~AsN0,1LJR2T⊗R1JH,whereAsN(μ,C) represents multidimensional normal distribution with mean value μ and covariance matrixC. ⊗ is Kronecker product.R1 andR2 are covariance matrices the subarrays X1 and X2, respectively. According to the nature of normal distribution, we have (38)W-0.5Δr~AsN0,I4M1M2+2M2,(39)W-0.5Δr22~Asχ24M1M2+2M2,whereW=1/LJ(R2T⊗R1)JH,Asχ2(4M1M2+2M2) represents Chi square distribution with DOF as (4M1M2+2M2). Then we introduce a parameter η to restrain Δr as follows: (40)W-0.5Δr22≤η.Let (40) hold with a high probabilityp=0.999. η can be resolved from chi2inv(p, 4M1M2+2M2) function of Matlab. Then (34) can be converted into the following robust form(41)min1TγsubjecttoReP~i,ImP~i2≤γiW-0.5r^-Aθ~P~22≤η,where1 isQ×1 dimensional vector whose elements all are 1, (41) can be resolved by CVX.R1 andR2 can be replaced by sample covariance matrix withL snapshots.(42)R^1=1L∑t=1Lx1tx1HtR^2=1L∑t=1Lx2tx2Ht.Through tracking nonzero position of coefficientγi, θi can be resolved. Next we can obtain the estimated coefficient vectorP as P^=A(θ)†r, where (•)+ denotes pseudo inverse operation. Thus, nominal elevation of theith source can be obtained as follows:(43)φi=arccos⁡arg⁡P^iπ.


## 3.5. Computational Procedure and Complexity Analysis

Now, our algorithm can be summarized as follows.Step 1.
Compute sample cross-correlation matrixR and sample covariance matricesR1 andR2 using (35) and (42).Step 2.
VectorizeR from (10) and reorder elements of the vector to obtain r^ according to (36).Step 3.
Calculate covariance matrix of residual error vector according to (37) and determine upper bound of constraint of η.Step 4.
Discretize nominal azimuth and calculate nominal elevationθi from (41).Step 5.
Calculate nominal elevationφi from (43).We analyze the computational complexity of the proposed method in comparison with COMET [22] which can be applied for uniform 2D arrays and Zhou’s algorithm [23] using double parallel arrays for 2D ID sources. For comparison, the two methods are supposed to be applied to double parallel arrays with the same sensor number as the proposed nested array. COMET [22] is a method of four-dimensional nonlinear optimization under alternating projection algorithm framework, its computational cost is mostly composed of the calculation of the sample covariance matrix which needs O(LN2) and the alternating projection technique with respect to cost functions which is O(3N3q3+ 3N2). Zhou’s algorithm calculates nominal elevation by a TLS-ESPRIT algorithm and finds nominal azimuth by 1D searching. Computational cost of Zhou’s algorithm [23] mainly contains calculation of the sample covariance matrix O(LN2), eigendecomposition, and inversion of the sample covariance matrix O(2N3) and 1D searching O(N3). Computational cost of the proposed algorithm mainly contains calculation of the sample cross-correlation matrix and sample covariance matrix O(LN2), solving second-order core programming O[(4M1M2+2M2)3Q3] and solving the nominal elevation by pseudo inversion O(N3). With more DOF the computational cost of the proposed algorithm is higher than two methods.


## 4. Results and Discussion

In this section, the effectiveness of the proposed algorithm is investigated through four numerical simulation experiments. We set the nested array withM1=2 andM2=3, the total sensor number isN=11,d=λ/2. The nominal azimuth space is divided from 0° to 180° in step of 1°. Root mean squared error (RMSE) is applied for the evaluation of the performance of estimation. TheRMSE of theith source with respect to DOA is defined as(44)RMSEθi=1Mc∑ςMcθ^iς-θi2,(45)RMSEφi=1Mc∑ςMcφ^iς-φi2,where θ^iς and φ^iς are the estimated nominal azimuth and estimated nominal elevation of theith ID source, respectively. The superscript ς denotes the estimated parameters in ςth Monte Carlo run.Mc is the number of Monte Carlo simulations. θi and φi are the true nominal azimuth and nominal elevation, respectively.In the first experiment, we examine the performance considering underdetermined scenario. 2D ID sources to be estimated are distributed from (15°, 15°) to (165°, 165°) uniformly. The center ofkth sources is located on [15°+(k-1)10°, 15°+(k-1)10°]. APDF of sources are Gaussian whenk is even while sources are uniform whenk is odd. Angular spreads are all set at 2°,SNR is set at 5db, and number of snapshots is 200. 100 independent Monte Carlo simulations are run to get the result. Figure 3 shows the positions of given 16 sources as well as average estimated positions of the 16 sources, which indicate that the proposed method can estimate sources with different APDF effectively in underdetermined scenario.Figure 3
Estimation for 2D ID sources with different APDF.In the second experiment, we examine the performance of the proposed method versus angular spreads of ID sources. As the rotational invariance relationship of generalize steering matrix is obtained under the small angular spreads assumption, it is necessary to investigate the influence of angular spreads on the estimation of the proposed algorithm. We consider both Gaussian and uniform sources, respectively.SNR is set at 5dB, number of snapshots is 200, andMc=100. Sources are centered at (30°, 45°) and (150°, 70°) with both angular spreads ranging from 0° to 10°.RMSE takes averages of two sources. Figure 4 shows thatRMSE reaches 0.6 as angular spreads are 5° andRMSE reaches 0.8 when angular spreads are 10° for Gaussian sources, whileRMSE reaches 0.75 and 1.1 when angular spreads reach 5° and 10° for uniform sources, respectively. It can be concluded that estimation accuracy deteriorates for both Gaussian and uniform sources as the angular spreads increase. The estimated results are still satisfactory within 5° and the proposed method shows robustness under small angular spreads.Figure 4
RMSE estimated for 2D ID sources with different angular spreads.In the third experiment, we examine the performance of the proposed method versusSNRand number of snapshots. Two Gaussian ID sources with parameters [30°, 45°, 2°, 2°] and [150°, 70°, 2°, 2°] are set to be estimated. Mc=100. Figure 5(a) showsRMSE of estimation withSNR ranging from -5dB to 30dB while number of snapshots is set at 200.RMSE takes averages of RMSEθi and RMSEφi of two sources. Figure 5(b) shows the influence of number of snapshots withSNRset at 5 dB. Sources are also estimated by COMET [22] and Zhou’s algorithm [23] based on two kinds of double parallel uniform linear arrays for comparison. One is composed of 12 sensors with 6 sensors separated by λ/2 meters in each subarray, which has one more sensor than nested array we proposed. The other has 30 sensors with each subarray containing 15 sensors separated by λ/2 meters. Physical aperture of the parallel uniform linear array with 30 sensors is similar to the proposed nested array. As can be seen from Figure 5, estimation accuracy increases asSNR or number of snapshots increases. As the aperture of nested array extended by the virtual array, the proposed method shows better performance than the two methods with the same number of sensors. Considering the parallel uniform linear array with 30 sensors, as shown in Figure 5(b),RMSE of COMET method is smaller than the proposed method when number of snapshots exceeds 200 and bigger than the proposed method when number of snapshots is less than 100. COMET separates unknown variables of each source based on alternating projection technique and formulates equations set of unknown variables from covariance matrix of receive vectors. Using more sensors, COMET shows better performance if the samples are sufficient. In the case of number of snapshots less than 100, the proposed method performs better as applying a robust sparse representation. In summary, with the same aperture, the proposed method has a better performance than uniform arrays in less snapshots circumstance; with same number of sensors, the proposed method has a better performance than uniform arrays under the same experimental conditions.Figure 5
(a)RMSE estimated by three algorithms for 2D ID sources versusSNR; (b)RMSE estimated by three algorithms for 2D ID sources versus number of snapshots.
(a) (b)In the fourth experiment, we investigate the impact of sensors number of subarrays on estimation. Two Gaussian ID sources with parameters [30°, 45°, 2°, 2°] and [150°, 70°, 2°, 2°] are set to be estimated.Mc=100,SNR=5dB, number of snapshots is 200. Figure 6 shows four curves representing 4 sensor number settings. The distances between adjacent sensorsd are all set at λ/2. Estimation accuracy will be improved as sensor number in subarrays increases; but this improvement is not significant when number of sensors in subarrays increases to a certain extent. Rotational invariance relationship is an approximate resolution under small angular spreads assumption; when sensors number is bigger, though aperture is larger, expression of receive signal with respect to later sensors in virtual array is a result of more transitive relationships, which means existing more errors.Figure 6
(a)RMSE estimated by different sensor number of subarrays versusSNR; (b)RMSE estimated by different sensor number of subarrays versus number of snapshots.
(a) (b)


## 5. Conclusions

In this paper a 2D nested array is proposed aiming at DOA estimation for 2D ID sources. Through vectorization of cross-correlation matrix by Khatri-Rao product, a virtual array is obtained. The close form of sensor position of virtual array and the optimal configuration of the nested array are derived. Rotational invariance relationship of generalized steering matrix with regard to nominal azimuth is deduced under the small angular spreads assumption. Sparse representation of DOA estimation for 2D ID sources based on the nested array is established and the solution process is detailed. Simulations investigate the influence of experiment conditions, angular spread, and sensor number on estimation. Results indicate that the proposed method can estimate more sources than sensors without information of APDF of sources and has better performance than uniform parallel arrays with the same sensor number.



---

*Source: 6941963-2019-07-03.xml*