# Direction-of-Arrival Estimation for 2D Coherently Distributed Sources with Nested Array Based on Matrix Reconstruction

Authors: Yu Xiao; Tao Wu; Yiwen Li; Xinping Ma; Yijie Huang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper has made proposition of a nested array and an estimation algorithm for direction-of-arrival (DOA) of two-dimensional (2D) coherently distributed (CD) sources. According to the difference coarray concept, double parallel hole-free virtual uniform linear arrays are generated by virtue of vectorization operation on cross-correlation matrices of subarrays. Sensor coordinates of virtual arrays are derived. Rational invariance relationships of virtual arrays are derived. According to the rotational invariance relationships, matrices satisfying rotation invariance are constructed by extracting and regrouping the receive vectors of the virtual arrays, and then an estimation of signal parameters via rotational invariance techniques- (ESPRIT-) like framework on matrix reconstruction is deduced. Optimal configuration of the nested array as well as computational complexity are analyzed. Without pair matching, the proposed method can resolve more sources than the sensor number. Simulation outcomes indicate that the proposed method tends to have a better performance as compared to the traditional uniform arrays that have similar number of sensors.


---

## Body

## 1. Introduction

Traditional direction-of-arrival (DOA) estimators consider a target as a point source. In underwater acoustic detection, with the reduction of the distance between the receive array and a target, as many parts of the target reflect signals; a point source model cannot describe the target effectively and the supposed conditions of a point source model are invalid. In this case, a distributed source model is proposed [1]. A distributed source can be considered as an aggregation of point sources within a spatial distribution, where the point sources can be called scatterers.Generally, distributed sources can be possibly categorized as incoherent distributed (ID) and coherent distributed (CD) sources according to the coherence of scatterers. The scatterers of an ID source can be assumed to be incoherent, while scatterers of a CD source are coherent. Distributed sources can be categorized into one-dimensional (1D) and two-dimensional (2D) according to the spatial distribution dimension. The 2D distributed sources with the assumption that scatterers of a distributed source and the receive array are not in a same plane is more general and accordant with practical circumstances. In this paper, we are concerned with 2D CD sources.The distributions of scatterers of CD sources can be described by deterministic angular signal distribution function (ASDF). ASDF of a source can be modeled as Gaussian, uniform, or any other form of distribution based on the geometric and surface characteristics of the source. ASDF of a 1D CD source involves parameters which are nominal angle and angular spread. ASDF of a 2D source is described by nominal azimuth, nominal elevation, azimuth spread, and elevation spread. Nominal azimuth as well as nominal elevation are collectively called nominal angles representing the target centers; these can also be possibly presented as DOA. Azimuth spread and elevation spread are called angular spreads denoting the spatial extension of a target.Considering point sources, some important achievements of DOA estimation have been presented in a mixed signal field and multiangle estimation recently. The 2D DOA estimation problem for a mixture of circular and noncircular signals is considered based on two parallel uniform linear arrays (ULAs) [2] and the uniform rectangular array (URA) [3]. In [4], authors have achieved estimating 2D DOA, 2D direction-of-departure, 2D receive polarization angle, and 2D transmit polarization angle simultaneously under an electromagnetic vector sensor (EMVS) MIMO system. A dimensional reduction noncircular MUSIC algorithm with a low computational complexity is proposed in [5] for DOA and polarization estimation of circular and noncircular mixture signals in a virtual MIMO system. As for ID sources, 2D DOA estimation of ID sources in massive multiple-input multiple-output (MIMO) systems using URAs is proposed in [6].As for CD sources, based on spectral search, several estimators have been provided from the classical point source technique multiple signal classification (MUSIC) such as distributed signal parameter estimator (DSPE) [1], dispersed signal parametric estimation (DISPARE) [7], and vec-multiple signal classification (vec-MUSIC) [8]. Shahbazpanahi et al. [9] have extended another classical point sources technique estimation of signal parameters via rotational invariance techniques (ESPRIT) to make an estimation of the nominal angles of sources in the first place and thereafter make calculations of the angular spreads by means of spectral search. Zoubir et al. have explored a generalized beamforming estimator in [10]. Xiong et al. [11] have analyzed the performance of DSPE algorithm considering that estimators for CD sources mostly approximate conclusions made on assumptions of small angular spreads. All these estimators are concerned on 1D CD sources where scatterers and receive arrays are in the same plane. However, scatterers and arrays are in three-dimensional space generally, which should be modeled as 2D CD sources.Involving more parameters, estimators for 2D CD sources have higher computational complexity. Extended form DSPE for 1D CD sources, Zhang et al. [12] have proposed a spectral search approach for 2D CD sources with L-shape arrays. Utilizing double parallel linear arrays (DPLA), treble parallel linear arrays (TPLA), and conformal arrays, several low-complexity algorithms have been opined in [13–18], where estimations of DOAs are performed under ESPRIT or propagator framework based on approximate rotational invariance relationships of steering vectors. Lee et al. [19] have proposed a sequential one-dimensional searching algorithm based on double uniform circular arrays (DUCA). In [20], DOAs can be estimated from proposed symmetric properties of a centrosymmetric crossed array. Considering CD sources are correlated with each other, Wu et al. [21] have made a proposition of a 2D estimator utilizing DPLA. Based on L-shaped arrays, Wu et al,[22] have explored the estimation of a 2D nonsymmetric CD source, where nonsymmetric ASDF is established by a Gaussian mixture model.Aforementioned methods mostly are proposed based on uniform arrays where the separations between sensors are limited to the value no more than half proportion of the wavelength of impinging signal. In the context of point source estimation, in order to achieve the greatest levels of accuracy as well as to have more sources resolved, larger apertures should be used; consequently, a larger number of sensors are required. Although point sources and distributed sources are different targets, the traditional estimators for distributed source also requires the separations between sensors not greater than half wavelength, and the estimation accuracy is closely related to the aperture size. The nested arrays which have been proposed in [23] tend to have greater degree of freedoms (DOFs) as well as huge apertures. Aiming at point sources, DOA estimators with diverse types of nested arrays have been explored in [24–29]. In [29], Zheng et al. have presented an estimator resolving near-field (NF) and far-field (FF) point sources simultaneously using a symmetric double nested array (SDNA). Zheng and Mu [30] have proposed a 2D DOA estimator using two parallel nested arrays where an augmented covariance matrix is constructed to reduce computational complexity. The coprime arrays composed of ULAs with sensor spacing related by a cprime integer have been proposed in [31]. Shi et al. [32] have presented 2D DOA estimation with coprime planar arrays (CPPAs) which achieve a great increase of aperture and reduce computational complexity by virtue of iterative scheme for transforming 2D grids searching into 1D searching during sparse representation. In [33], Zheng et al. have proposed a new sparse array called configuration maximum interelement spacing constraint (MISC) array where the array structure is designed in terms of the interelement spacing set, which is given in a closed-form expression as a function of the number of sensors.Compared with estimators for a point source model, the researches of estimators for distributed sources under sparse arrays are relatively few and most of these studies deal with 1D distributed source. Estimators for 1D exponential distributed sources based on the linear nested array have been presented in [34] with a prior knowledge of angular spreads and without the prior knowledge in [35]. Considering 2D ID sources, Wu et al. [36] have developed a sparse representation method under nested array, where computational complexity is extremely high.In this paper, based on a proposed nested array containing double parallel linear subarrays, a DOA estimator for 2D CD sources is presented. By means of vectorization operation on the cross-correlation matrices, double hole-free virtual uniform linear arrays are realized and the closed expressions of virtual sensor coordinates are deduced based on the difference coarray concept. Next, rotational invariance relationships of virtual arrays are deduced. Based on rotational invariance relationships, matrices satisfying rotation invariance are constructed by extracting and regrouping the receive vectors of the virtual arrays. Then, an ESPRIT-like framework based on matrix reconstruction is proposed. The proposed method can detect 2D CD sources more than sensors without angles matching and without prior knowledge of ASDF. Lastly, the effectiveness of the proposed method is investigated through simulations. To show the contributions of this paper clearly, the main differences between the state-of-the-art methods and our work are listed as follows:(i)
Though DOA estimations for 2D CD sources are presented in [12–22], the arrays are all uniform. In [27, 28, 30, 32], DOA estimators for 2D sources using sparse arrays are considered, which aim at point sources. In [34, 35], estimators for 1D exponential distributed sources based on nested arrays are established, while 2D distributed sources are more general and accordant with practical circumstances. In this paper, estimation for 2D CD sources with the nested array is considered.(ii)
A nested array used for 2D sources is proposed. Configuration of virtual arrays containing close form of virtual sensor positions, maximum DOF, and the optimal parameters get deduced. Compared with uniform linear arrays for 2D sources, the nested array proposed has larger aperture as well as more DOFs.(iii)
A matrix reconstruction algorithm is detailed based on an ESPRIT-like framework on matrices which is constructed by extracting receive vectors of the virtual arrays. Without angles matching and without prior knowledge of ASDF, the algorithm proposed has less computation complexity than existing methods using uniform arrays.Notations. Scalar variables are denoted by italic letters, and vectors and matrices are denoted by bold letters. (•)−1, (•)∗, (•)T, and (•)H mean inverse, complex conjugate, transpose, and Hermitian transpose of a matrix. E[•], (./), vec(•), ⊙ and (•)+ denote expectation, element-wise division, vectorization, Khatri–Rao product and pseudoinverse operation. [•]k is the kth element of a vector. diag[•] means diagonal matrix. arg(•) is the phase of a complex number.


## 2. The Nested Array and Source Model

As illustrated in Figure1, double parallel uniform subarrays A1 and A2 on the xoz plane constitute a nested array. Parallel to x-axis, A1 is made up of a middle sensor on z-axis and M1 sensors are separated by d meters on both sides. Subarray A2 is on x-axis containing M2 sensors located on both semiaxes with sensor spacing D = (2M1 + 1)d. The interval of the two subarrays equals d/2. q 2D CD sources with wavelength λ and nominal angles (θi, φi) (i = 1, 2, …, q) in the far field are set to be estimated, θi∈0,π, φi∈0,π. θi is a denotation of the spatial nominal azimuth that exists between x positive semiaxes and the propagation path of the ith source. φi is the nominal elevation of the source. The consideration of noise is the Gaussian white zero mean which is an additive as well as uncorrelated with sensors.Figure 1
(a) Configuration of the nested array proposed; (b) illustration of the array proposed on thexoz plane.
(a)
(b)The sensor coordinate set of subarraysA1 and A2 can be expressed, respectively, as(1)C1=lmd,0,0.5dlm=−M1,−M1−1,…,M1−1,M1,(2)C2=ln+0.5D,0,0ln=−M2,−M2−1,…,M2−1.Receive vectors of subarraysA1 and A2 have the following expressions:(3)x1t=∑i=1q∬η1θ,φsiθ,φ,tdθdφ+nx1t,(4)x2t=∑i=1q∬η2θ,φsiθ,φ,tdθdφ+nx2t,where nx1(t)and nx2(t) are noise vectors. η1(θ,φ) and η2(θ,φ) denote the steering vectors of subarrays A1 and A2, which have the following expressions:(5)η1θ,ϕ=ejπdcosϕ/λe−j2πM1dcosθ/λ,e−j2πM1−1dcosθ/λ,…,1,…,ej2πM1−1dcosθ/λ,ej2πM1dcosθ/λT.(6)η2θ,ϕ=e−j2πM2−0.5Dcosθ/λ,e−j2πM2−1.5Dcosθ/λ,…,ej2πM2−1.5Dcosθ/λ,ej2πM2−0.5Dcosθ/λT.s
i(θ, φ, t) in equations (3) and (4) represents angular signal density function of the ith 2D CD source. Based on the assumption of the CD source, si(θ, φ, t) can be written as(7)siθ,φ,t=sitgθ,φ;ui,where si(t) reflects time behavior of the ith source. giθ,φ;ui represents deterministic ASDF of the ith source reflecting spatial distribution of scatterers with respect to the ith source. A 2D deterministic ASDF is generally characterized by the parameter vector ui = [θi, φi, σθi, σφi]T with four elements which are nominal azimuth θi, nominal elevation φi, azimuth angular spread σθi, and elevation angular spread σφi.Deterministic ASDF of a Gaussian CD source has the following expression:(8)gθ,φ;ui=12πσθiσϕiexp−0.5θ−θiσθi2+ϕ−ϕiσϕi2.Deterministic ASDF of a uniform CD source have the following expression:(9)gθ,φ;ui=14σθiσϕi,θ−θi≤σθiandϕ−ϕi≤σϕi,0,θ−θi≥σθiorϕ−ϕi≥σϕi.


## 3. Proposed Method

This part is composed of four portions. Firstly, double parallel virtual uniform linear arrays are realized through difference operation. Then, the rotational invariance relationship of virtual arrays gets deduced. Afterwards, matrices satisfying rotation invariance are constructed using the receive vectors of the virtual arrays and a matrix reconstruction algorithm based on an ESPRIT-like framework is proposed. The procedure of computation gets summarized as well as the analysis of the complexity is performed by comparing other methods with traditional uniform arrays.


### 3.1. Description of Virtual Arrays

The nested array increases DOFs by forming virtual arrays through difference operation between sensor coordinates of subarrays. (2M1 + 1)×2M2 dimensional matrix R1 is denoted as the cross-correlation matrix of subarrays A1 and A2 and the 2M2 × (2M1 + 1) dimensional matrix R2 as the cross-correlation matrix of subarrays A2 and A1. As different CD sources are uncorrelated, R1 and R2 can be expressed as follows:(10)R1=Ex1tx2Ht=∑i=1qσi2∫∫η1θ,φgθ,φ;uigθ,φ;uiη2Hθ,φdθdφ,(11)R2=Ex2tx1Ht=∑i=1qσi2∫∫η2θ,φgθ,φ;uigθ,φ;uiη1Hθ,φdθdφ,where σi2=Esit2 is the power of the ith source. Vectorizing R1 and R2 by the Khatri–Rao product, we have(12)r¯1=vecR1=∑i=1qσi2∫∫η2∗θ,φ⊙η1θ,φgθ,φ;uigθ,φ;uidθdφ,(13)r¯2=vecR2=∑i=1qσi2∫∫η1∗θ,φ⊙η2θ,φgθ,φ;uigθ,φ;uidθdφ.Double parallel virtual uniform linear arrays can be obtained by vectorization of cross-correlation matricesR1 and R2 which is illustrated in Figure 2. Through difference operation between sensor coordinates of subarrays A1 and A2, virtual arrays VA and VB can be realized with sensor coordinate sets ℂA and ℂB expressed as(14)ℂA=pi−pjpi∈ℂ1,pj∈ℂ2,(15)ℂB=pj−pipi∈ℂ1,pj∈ℂ2.Figure 2
Virtual arrays formed by nested array.From equations (1) and (2), we have sensor coordinate sets of the virtual arrays ℂA and ℂB as(16)ℂA=n−0.5d,0,0.5dn=1−2M1M2+M2,2−2M1M2+M2,…,2M1M2+M2,(17)ℂB=n−0.5d,0,−0.5dn=1−2M1M2+M2,2−2M1M2+M2,…,2M1M2+M2.Equations (16) and (17) show that each virtual array contains L = 2(2M1M2 + M2) virtual sensors without repetitive positions, which indicates both virtual arrays formed by nested array are hole-free. It is worth noting that the element sequence in r¯1 does not correspond to sensor orders of virtual array VA; meanwhile, the element sequence in r¯2 does not correspond to sensor orders of VB. rA is denoted as receive signal vector with elements reflecting signal received from the virtual sensor located on [(0.5 − 2M1M2 − M2)d, 0, 0.5d] to the sensor located on [(2M1M2 + M2 − 0.5)d, 0, 0.5d] and rB as receive signal vector with elements reflecting signal received from the virtual sensor located on [(0.5 − 2M1M2 − M2)d, 0, −0.5d] to the sensor located on [(2M1M2 + M2 − 0.5)d, 0, −0.5d]. After reordering r¯1 and r¯2, rA and rB can be expressed as(18)rA=∑i=1qσi2∫∫aθ,φgθ,φ;uigθ,φ;uidθdφ,(19)rB=∑i=1qσi2∫∫bθ,φgθ,φ;uigθ,φ;uidθdφ,where a(θ, φ) and b(θ, φ) denote the L × 1 dimensional steering vector of virtual arrays VA and VB, which have the following expressions:(20)aθ,φ=ejπdcosφ/λe−j2π2M1M2+M2−0.5dcosθ/λ,e−j2π2M1M2+M2−1.5dcosθ/λ,…,ej2π2M1M2+M2−0.5dcosθ/λT,(21)bθ,φ=e−jπdcosφ/λe−j2π2M1M2+M2−0.5dcosθ/λ,e−j2π2M1M2+M2−1.5dcosθ/λ,…,ej2π2M1M2+M2−0.5dcosθ/λT.In the process of estimation, the sample of covariance matrix withN snapshots can be substituted for R1 and R2 as follows:(22)R^1=1N∑t=1Nx1tx2Ht,R^2=1N∑t=1Nx2tx1Ht.


### 3.2. Rotational Invariance Relationships of Virtual Arrays

Representing response of arrays to theith distributed source, generalized steering vectors are defined in [10]. The generalized steering vector αi reflects response of array VA to the ith CD source while the generalized steering vector βi reflects response of array VB to the source. αi and βi can be written as(23)αi=∬aθ,φgθ,φ;uigθ,φ;uidθdφ,(24)βi=∬bθ,φgθ,φ;uigθ,φ;uidθdφ.Each element ofαi or βi represents the response of corresponding virtual sensor of VA and VB to the ith CD source. Representing response of arrays to all distributed sources, A and B are defined as L × q dimensional generalized steering matrices of the virtual arrays VA and VB, which can be written as(25)A=α1,α2,…,αq,(26)B=β1,β2,…,βq.σ is denoted as a power vector of CD sources, which can be expressed as(27)σ=σ12,σ22,…,σq2T.Thus,rA and rB can be expressed as(28)rA=Aσ,(29)rB=Bσ.Ifd/λ is set at 1/2, the following relationships can be obtained under the condition that angular spreads of CD sources are small (the proof can be seen in Appendix):(30)αik+1≈ejπcosθiαik,(31)βik+1≈ejπcosθiβik,(32)αik≈ejπcosφiβik.Therefore, rotational invariance relationships can be drawn as follows:(33)Ak+1≈AkΦ,(34)Bk+1≈BkΦ,(35)Ak≈BkΨ,where the rotational invariance operators Φ and Ψ can be written as(36)Φ=diagejπcosθ1,ejπcosθ2,…,ejπcosθq,(37)Ψ=diagejπcosφ1,ejπcosφ2,…,ejπcosφq.


### 3.3. DOA Estimation via Matrix Reconstruction

Receive vectors described by equations (28) and (29) can be considered q coherent power signals impinging into virtual arrays VA and VB. In this paper, according to the concept of ESPRIT, matrices satisfying rotation invariance are constructed by extracting and regrouping the receive vectors of the virtual arrays.DefiningrA1 as the receive vector of the first q elements of rA and rB1 as the receive vector of the first q elements of rB, we have(38)rA1=A1σ,(39)rB1=B1σ,where A1 and B1 are generalized steering matrices of rA1 and rB1. A1 and B1 can be expressed as(40)A1=Iq×q0q×L−qA,(41)B1=Iq×q0q×L−qB,where Iq×q is the q × q dimensional identity matrix and 0q×(L-−q) is the q × (L − q) dimensional null matrix. The operation Iq×q0q×L−q left multiply A or B can be considered as extracting the first q rows of a matrix.According to equations (33)–(35), a generalized steering matrix of the receive vector of ith q elements of rA and rB can be expressed as(42)Ai≈A1Φi−1,i=1,2,…,q+1,(43)Bi≈B1Φi−1,i=1,2,…,q+1,(44)B1≈A1Ψ.Then, the receive vector ofith q elements of rA and rB can be expressed as(45)rAi=Aiσ≈A1Φi−1σ,(46)rBi=Biσ≈B1Φi−1σ.Aq × q dimensional matrix RA′ is constructed in the way that rAi is converted into the ith column of RA′. Another q × q dimensional matrix RA″ is reconstructed in the way that rAi+1 is converted into the ith columns of RA″. Thus, RA′ and RA″ can be written as follows:(47)RA′=∑i=1qrAieiT=A1Λ∑i=1qQi−1eiT=A1ΛA1H,(48)RA″=∑i=2q+1rAieiT=A1Λ∑i=2q+1Qi−1eiT=A1ΦΛA1H,where Λ=diagσ12,σ22,…,σq2 and Q=ejπcosθ1,ejπcosθ2,…,ejπcosθqT. ei=0,…,1,0,… is the q dimensional vector with 1 as the ith element and 0 as others, which can be seen as an operator transforming a vector into a column in a matrix.Similarly,RB′ is constructed in the way that rBi is converted into the ith column of the matrix.(49)RB′=∑i=1qrBieiT=B1Λ∑i=1qQi−1eiT=B1ΛA1H=A1ΨΛA1H,where A1H is a Vandermonde full rank matrix and Λ is a diagonal matrix. From equations (48) and (49), we have(50)RA′+RA″=A1ΦA1−1,(51)RA′+RB′=A1ΨA1−1.R
A
′
+
R
A
″ and RA′+RB′ are denoted as Ω1 and Ω2, respectively. The rotational invariance operators Φ and Ψ are diagonal matrices. According to equations (50) and (51), the eigenvectors corresponding to the elements in the same position in the two diagonal matrices Φ and Ψ are the same, which means eigenvalue ejπcosθi of Ω1 and eigenvalue ejπcosφi of Ω2 have a same eigenvector. Therefore, once we obtain ejπcosθi through eigendecomposition of Ω1, eigenvectors corresponding to ejπcosθi can be substituted in (51) to calculate the ejπcosφi paired with ejπcosθi.We can obtain eigenvalueμi (i = 1, 2, …, q) of Ω1 and its corresponding eigenvectors ξi by means of eigendecomposition. Then, we can obtain the nominal azimuth as follows:(52)θi=arccosargμiπ,i=1,2,…,q.As corresponding eigenvalues ofΩ1 and Ω2 have the same eigenvector, without angles matching, eigenvalue of Ω2 can be obtained as follows:(53)νi=1q11×q⋅Ω2ξi⋅/ξi,where 11×q is the 1 × q dimensional vector with all elements as 1 and (./) denotes element-wise division. Ω2ξi equals viξi as the ith eigenvalue μi of Ω1, and the ith eigenvalue vi of Ω2 has the same eigenvector. The operation of element-wise division between Ω2ξi and ξi obtains a q × 1dimensional vector with all elements equal to vi theoretically.Then, we can obtain the nominal elevation as follows:(54)ϕi=arccosargviπ,i=1,2,…,q.


### 3.4. Analysis of the Nested Array

Based on the analysis of Section 3.3, maximum number of unknown CD sourcesq can achieve half of the virtual sensors of each virtual arrays; thus, DOF of the arrays proposed is (2M1M2 + M2); total sensor number of the array is M = 2(M1 + M2) + 1. We can draw the conclusion that DOF of the nested array proposed can reach O(M2) and the DOF is far more than the total sensor number. We examine the maximum DOF that can be achieved, and the maximum goal under the given total sensor number can be written as(55)max2M1M2+M2subjecttoM=2M1+M2+1.We obtain(56)M1opt=M4−12,M2opt=M4.AsM1opt and M2opt could be fractions, M1, M2, and (M − 1)/2 are all integers, and optimal solutions should be modified. Refer to the definition in [29], coarray aperture of the proposed nested array is defined as physical length in x dimension multiply physical length in z dimension of virtual arrays. Table 1 shows maximum DOF and the coarray aperture of the optimal configuration.Table 1
Maximum DOF and the coarray aperture of the optimal configuration.
N (n is an integer)
M
1
M
2
DOF
Coarray aperture
4n + 1
N
n
2n2 + n
(4n2 + 2n − 1)d × d
4n + 3
N
n +1
2n2 + 3n + 2
(4n2 + 6n + 1)d × dWe analyze DOFs of the method proposed in comparison with representative existing methods for 2D CD sources with the same sensor numberM which are DSPE using a L-shape array [12], modified propagator (MP) using TPLA, [15] and SOS using DUCA [19]. According to the optimal configuration described in Table 1, DOFs of the method proposed are shown in Figure 3. It can be seen that the nested arrays proposed have a larger number of DOFs than uniform arrays with the same sensor number.Figure 3
DOFs of different methods.The algorithm proposed can be performed according to Table2.Table 2
Summary of the proposed method.
Input:N snapshots of receive vectors of the subarrays A1 and A2 array x1tt=1N and x2tt=1N
Output: nominal angles ofq 2D CD sources (θi, φi) (i = 1, 2, …, q).
(1) Computing the sample cross-correlation matrix of subarraysR1 and R2 from equations (22).
(2) VectorizingR1 and R2 from equations (12) and (13) and reordering sequences of elements in vectors r¯1 and r¯2 to obtain rA and rB.
(3) Extracting the elements of vectorsrA and rB to get rAi and rBi and constructing RA′, RA″, and RB′ form equations (47)–(49).
(4) Calculation and eigendecomposition ofΩ1 and Ω2.
(5) Computing nominal elevationθi from equation (52) and nominal elevation φi from equation (54).We analyze the computational complexity of the method proposed in comparison with representative existing methods for 2D CD sources based on uniform arrays with the same sensor numberM which are DSPE using a L-shape array [12], Modified Propagator (MP) using TPLA [15], and SOS using DUCA [19]. Computational cost of DSPE using a L-shape array [12] mainly lies in computing the sample covariance matrix O(NM2), pseudoinverse operation of the covariance matrix O(M3), and four-dimensional searching O(M3K4), where K denotes search points of unknown parameters. Computational cost of MP using TPLA [15] mainly lies in computing the modified propagator O(NM2) and inversion as well as eigendecomposition of the sample covariance matrix O(M3), obtaining rotational invariance matrices O(q3). SOS using DUCA [19] is an algorithm based on 1D searching; computational cost mainly lies in computing the sample covariance matrix O(NM2) and inversion as well as eigendecomposition of the sample covariance matrix O(M3), one-dimensional searching O(M3K). The computational cost of the method proposed mostly lie in computing the sample covariance matrix O(NM2), calculation of pseudoinversion of RA′ which is O(q3), and eigendecomposition of Ω1 and Ω2 which is O(q3). As M >q, the computation complexity of the proposed method is far less than existing methods for 2D CD sources using uniform arrays.Remark 1.
The advantages of the proposed method lies in that the proposed arrays have larger aperture as well as more DOFs compared with traditional estimators for CD sources based on uniform arrays with similar number of sensors. Compared with existing estimators for CD sources based on nested arrays, the proposed method can estimate 2D CD sources and the estimation can be performed without prior knowledge of ASDF of sources as well as spectral searching.Remark 2.
The method proposed is not suitable for 2D CD sources which are coherent with each other. A CD source is defined as the scatterers of a source which are coherent; meanwhile, the sources are incoherent. In other words signals emitted from scatterers within a source are coherent while signals from scatterers from different sources are incoherent. Under the condition that CD sources are incoherent, we can obtainR1 and R2 in the forms of (10) and (11) which are the basis of vectorization to obtain the virtual receive vector of virtual arrays.Remark 3.
Although the array structure proposed achieves aperture expansion in both the elevation and azimuth dimensions, the aperture expansion in the elevation direction is limited. In addition, rotational invariance relationship is an approximate resolution under small angular spread assumption; when the sensor number is bigger, later sensors in the virtual array are a result of more transitive relationships. The abovementioned two factors indicate that the proposed method can improve the estimation accuracy of DOA of 2D CD sources to a certain extent.


## 3.1. Description of Virtual Arrays

The nested array increases DOFs by forming virtual arrays through difference operation between sensor coordinates of subarrays. (2M1 + 1)×2M2 dimensional matrix R1 is denoted as the cross-correlation matrix of subarrays A1 and A2 and the 2M2 × (2M1 + 1) dimensional matrix R2 as the cross-correlation matrix of subarrays A2 and A1. As different CD sources are uncorrelated, R1 and R2 can be expressed as follows:(10)R1=Ex1tx2Ht=∑i=1qσi2∫∫η1θ,φgθ,φ;uigθ,φ;uiη2Hθ,φdθdφ,(11)R2=Ex2tx1Ht=∑i=1qσi2∫∫η2θ,φgθ,φ;uigθ,φ;uiη1Hθ,φdθdφ,where σi2=Esit2 is the power of the ith source. Vectorizing R1 and R2 by the Khatri–Rao product, we have(12)r¯1=vecR1=∑i=1qσi2∫∫η2∗θ,φ⊙η1θ,φgθ,φ;uigθ,φ;uidθdφ,(13)r¯2=vecR2=∑i=1qσi2∫∫η1∗θ,φ⊙η2θ,φgθ,φ;uigθ,φ;uidθdφ.Double parallel virtual uniform linear arrays can be obtained by vectorization of cross-correlation matricesR1 and R2 which is illustrated in Figure 2. Through difference operation between sensor coordinates of subarrays A1 and A2, virtual arrays VA and VB can be realized with sensor coordinate sets ℂA and ℂB expressed as(14)ℂA=pi−pjpi∈ℂ1,pj∈ℂ2,(15)ℂB=pj−pipi∈ℂ1,pj∈ℂ2.Figure 2
Virtual arrays formed by nested array.From equations (1) and (2), we have sensor coordinate sets of the virtual arrays ℂA and ℂB as(16)ℂA=n−0.5d,0,0.5dn=1−2M1M2+M2,2−2M1M2+M2,…,2M1M2+M2,(17)ℂB=n−0.5d,0,−0.5dn=1−2M1M2+M2,2−2M1M2+M2,…,2M1M2+M2.Equations (16) and (17) show that each virtual array contains L = 2(2M1M2 + M2) virtual sensors without repetitive positions, which indicates both virtual arrays formed by nested array are hole-free. It is worth noting that the element sequence in r¯1 does not correspond to sensor orders of virtual array VA; meanwhile, the element sequence in r¯2 does not correspond to sensor orders of VB. rA is denoted as receive signal vector with elements reflecting signal received from the virtual sensor located on [(0.5 − 2M1M2 − M2)d, 0, 0.5d] to the sensor located on [(2M1M2 + M2 − 0.5)d, 0, 0.5d] and rB as receive signal vector with elements reflecting signal received from the virtual sensor located on [(0.5 − 2M1M2 − M2)d, 0, −0.5d] to the sensor located on [(2M1M2 + M2 − 0.5)d, 0, −0.5d]. After reordering r¯1 and r¯2, rA and rB can be expressed as(18)rA=∑i=1qσi2∫∫aθ,φgθ,φ;uigθ,φ;uidθdφ,(19)rB=∑i=1qσi2∫∫bθ,φgθ,φ;uigθ,φ;uidθdφ,where a(θ, φ) and b(θ, φ) denote the L × 1 dimensional steering vector of virtual arrays VA and VB, which have the following expressions:(20)aθ,φ=ejπdcosφ/λe−j2π2M1M2+M2−0.5dcosθ/λ,e−j2π2M1M2+M2−1.5dcosθ/λ,…,ej2π2M1M2+M2−0.5dcosθ/λT,(21)bθ,φ=e−jπdcosφ/λe−j2π2M1M2+M2−0.5dcosθ/λ,e−j2π2M1M2+M2−1.5dcosθ/λ,…,ej2π2M1M2+M2−0.5dcosθ/λT.In the process of estimation, the sample of covariance matrix withN snapshots can be substituted for R1 and R2 as follows:(22)R^1=1N∑t=1Nx1tx2Ht,R^2=1N∑t=1Nx2tx1Ht.


## 3.2. Rotational Invariance Relationships of Virtual Arrays

Representing response of arrays to theith distributed source, generalized steering vectors are defined in [10]. The generalized steering vector αi reflects response of array VA to the ith CD source while the generalized steering vector βi reflects response of array VB to the source. αi and βi can be written as(23)αi=∬aθ,φgθ,φ;uigθ,φ;uidθdφ,(24)βi=∬bθ,φgθ,φ;uigθ,φ;uidθdφ.Each element ofαi or βi represents the response of corresponding virtual sensor of VA and VB to the ith CD source. Representing response of arrays to all distributed sources, A and B are defined as L × q dimensional generalized steering matrices of the virtual arrays VA and VB, which can be written as(25)A=α1,α2,…,αq,(26)B=β1,β2,…,βq.σ is denoted as a power vector of CD sources, which can be expressed as(27)σ=σ12,σ22,…,σq2T.Thus,rA and rB can be expressed as(28)rA=Aσ,(29)rB=Bσ.Ifd/λ is set at 1/2, the following relationships can be obtained under the condition that angular spreads of CD sources are small (the proof can be seen in Appendix):(30)αik+1≈ejπcosθiαik,(31)βik+1≈ejπcosθiβik,(32)αik≈ejπcosφiβik.Therefore, rotational invariance relationships can be drawn as follows:(33)Ak+1≈AkΦ,(34)Bk+1≈BkΦ,(35)Ak≈BkΨ,where the rotational invariance operators Φ and Ψ can be written as(36)Φ=diagejπcosθ1,ejπcosθ2,…,ejπcosθq,(37)Ψ=diagejπcosφ1,ejπcosφ2,…,ejπcosφq.


## 3.3. DOA Estimation via Matrix Reconstruction

Receive vectors described by equations (28) and (29) can be considered q coherent power signals impinging into virtual arrays VA and VB. In this paper, according to the concept of ESPRIT, matrices satisfying rotation invariance are constructed by extracting and regrouping the receive vectors of the virtual arrays.DefiningrA1 as the receive vector of the first q elements of rA and rB1 as the receive vector of the first q elements of rB, we have(38)rA1=A1σ,(39)rB1=B1σ,where A1 and B1 are generalized steering matrices of rA1 and rB1. A1 and B1 can be expressed as(40)A1=Iq×q0q×L−qA,(41)B1=Iq×q0q×L−qB,where Iq×q is the q × q dimensional identity matrix and 0q×(L-−q) is the q × (L − q) dimensional null matrix. The operation Iq×q0q×L−q left multiply A or B can be considered as extracting the first q rows of a matrix.According to equations (33)–(35), a generalized steering matrix of the receive vector of ith q elements of rA and rB can be expressed as(42)Ai≈A1Φi−1,i=1,2,…,q+1,(43)Bi≈B1Φi−1,i=1,2,…,q+1,(44)B1≈A1Ψ.Then, the receive vector ofith q elements of rA and rB can be expressed as(45)rAi=Aiσ≈A1Φi−1σ,(46)rBi=Biσ≈B1Φi−1σ.Aq × q dimensional matrix RA′ is constructed in the way that rAi is converted into the ith column of RA′. Another q × q dimensional matrix RA″ is reconstructed in the way that rAi+1 is converted into the ith columns of RA″. Thus, RA′ and RA″ can be written as follows:(47)RA′=∑i=1qrAieiT=A1Λ∑i=1qQi−1eiT=A1ΛA1H,(48)RA″=∑i=2q+1rAieiT=A1Λ∑i=2q+1Qi−1eiT=A1ΦΛA1H,where Λ=diagσ12,σ22,…,σq2 and Q=ejπcosθ1,ejπcosθ2,…,ejπcosθqT. ei=0,…,1,0,… is the q dimensional vector with 1 as the ith element and 0 as others, which can be seen as an operator transforming a vector into a column in a matrix.Similarly,RB′ is constructed in the way that rBi is converted into the ith column of the matrix.(49)RB′=∑i=1qrBieiT=B1Λ∑i=1qQi−1eiT=B1ΛA1H=A1ΨΛA1H,where A1H is a Vandermonde full rank matrix and Λ is a diagonal matrix. From equations (48) and (49), we have(50)RA′+RA″=A1ΦA1−1,(51)RA′+RB′=A1ΨA1−1.R
A
′
+
R
A
″ and RA′+RB′ are denoted as Ω1 and Ω2, respectively. The rotational invariance operators Φ and Ψ are diagonal matrices. According to equations (50) and (51), the eigenvectors corresponding to the elements in the same position in the two diagonal matrices Φ and Ψ are the same, which means eigenvalue ejπcosθi of Ω1 and eigenvalue ejπcosφi of Ω2 have a same eigenvector. Therefore, once we obtain ejπcosθi through eigendecomposition of Ω1, eigenvectors corresponding to ejπcosθi can be substituted in (51) to calculate the ejπcosφi paired with ejπcosθi.We can obtain eigenvalueμi (i = 1, 2, …, q) of Ω1 and its corresponding eigenvectors ξi by means of eigendecomposition. Then, we can obtain the nominal azimuth as follows:(52)θi=arccosargμiπ,i=1,2,…,q.As corresponding eigenvalues ofΩ1 and Ω2 have the same eigenvector, without angles matching, eigenvalue of Ω2 can be obtained as follows:(53)νi=1q11×q⋅Ω2ξi⋅/ξi,where 11×q is the 1 × q dimensional vector with all elements as 1 and (./) denotes element-wise division. Ω2ξi equals viξi as the ith eigenvalue μi of Ω1, and the ith eigenvalue vi of Ω2 has the same eigenvector. The operation of element-wise division between Ω2ξi and ξi obtains a q × 1dimensional vector with all elements equal to vi theoretically.Then, we can obtain the nominal elevation as follows:(54)ϕi=arccosargviπ,i=1,2,…,q.


## 3.4. Analysis of the Nested Array

Based on the analysis of Section 3.3, maximum number of unknown CD sourcesq can achieve half of the virtual sensors of each virtual arrays; thus, DOF of the arrays proposed is (2M1M2 + M2); total sensor number of the array is M = 2(M1 + M2) + 1. We can draw the conclusion that DOF of the nested array proposed can reach O(M2) and the DOF is far more than the total sensor number. We examine the maximum DOF that can be achieved, and the maximum goal under the given total sensor number can be written as(55)max2M1M2+M2subjecttoM=2M1+M2+1.We obtain(56)M1opt=M4−12,M2opt=M4.AsM1opt and M2opt could be fractions, M1, M2, and (M − 1)/2 are all integers, and optimal solutions should be modified. Refer to the definition in [29], coarray aperture of the proposed nested array is defined as physical length in x dimension multiply physical length in z dimension of virtual arrays. Table 1 shows maximum DOF and the coarray aperture of the optimal configuration.Table 1
Maximum DOF and the coarray aperture of the optimal configuration.
N (n is an integer)
M
1
M
2
DOF
Coarray aperture
4n + 1
N
n
2n2 + n
(4n2 + 2n − 1)d × d
4n + 3
N
n +1
2n2 + 3n + 2
(4n2 + 6n + 1)d × dWe analyze DOFs of the method proposed in comparison with representative existing methods for 2D CD sources with the same sensor numberM which are DSPE using a L-shape array [12], modified propagator (MP) using TPLA, [15] and SOS using DUCA [19]. According to the optimal configuration described in Table 1, DOFs of the method proposed are shown in Figure 3. It can be seen that the nested arrays proposed have a larger number of DOFs than uniform arrays with the same sensor number.Figure 3
DOFs of different methods.The algorithm proposed can be performed according to Table2.Table 2
Summary of the proposed method.
Input:N snapshots of receive vectors of the subarrays A1 and A2 array x1tt=1N and x2tt=1N
Output: nominal angles ofq 2D CD sources (θi, φi) (i = 1, 2, …, q).
(1) Computing the sample cross-correlation matrix of subarraysR1 and R2 from equations (22).
(2) VectorizingR1 and R2 from equations (12) and (13) and reordering sequences of elements in vectors r¯1 and r¯2 to obtain rA and rB.
(3) Extracting the elements of vectorsrA and rB to get rAi and rBi and constructing RA′, RA″, and RB′ form equations (47)–(49).
(4) Calculation and eigendecomposition ofΩ1 and Ω2.
(5) Computing nominal elevationθi from equation (52) and nominal elevation φi from equation (54).We analyze the computational complexity of the method proposed in comparison with representative existing methods for 2D CD sources based on uniform arrays with the same sensor numberM which are DSPE using a L-shape array [12], Modified Propagator (MP) using TPLA [15], and SOS using DUCA [19]. Computational cost of DSPE using a L-shape array [12] mainly lies in computing the sample covariance matrix O(NM2), pseudoinverse operation of the covariance matrix O(M3), and four-dimensional searching O(M3K4), where K denotes search points of unknown parameters. Computational cost of MP using TPLA [15] mainly lies in computing the modified propagator O(NM2) and inversion as well as eigendecomposition of the sample covariance matrix O(M3), obtaining rotational invariance matrices O(q3). SOS using DUCA [19] is an algorithm based on 1D searching; computational cost mainly lies in computing the sample covariance matrix O(NM2) and inversion as well as eigendecomposition of the sample covariance matrix O(M3), one-dimensional searching O(M3K). The computational cost of the method proposed mostly lie in computing the sample covariance matrix O(NM2), calculation of pseudoinversion of RA′ which is O(q3), and eigendecomposition of Ω1 and Ω2 which is O(q3). As M >q, the computation complexity of the proposed method is far less than existing methods for 2D CD sources using uniform arrays.Remark 1.
The advantages of the proposed method lies in that the proposed arrays have larger aperture as well as more DOFs compared with traditional estimators for CD sources based on uniform arrays with similar number of sensors. Compared with existing estimators for CD sources based on nested arrays, the proposed method can estimate 2D CD sources and the estimation can be performed without prior knowledge of ASDF of sources as well as spectral searching.Remark 2.
The method proposed is not suitable for 2D CD sources which are coherent with each other. A CD source is defined as the scatterers of a source which are coherent; meanwhile, the sources are incoherent. In other words signals emitted from scatterers within a source are coherent while signals from scatterers from different sources are incoherent. Under the condition that CD sources are incoherent, we can obtainR1 and R2 in the forms of (10) and (11) which are the basis of vectorization to obtain the virtual receive vector of virtual arrays.Remark 3.
Although the array structure proposed achieves aperture expansion in both the elevation and azimuth dimensions, the aperture expansion in the elevation direction is limited. In addition, rotational invariance relationship is an approximate resolution under small angular spread assumption; when the sensor number is bigger, later sensors in the virtual array are a result of more transitive relationships. The abovementioned two factors indicate that the proposed method can improve the estimation accuracy of DOA of 2D CD sources to a certain extent.


## 4. Results and Discussion

Three numerical simulation experiments are carried out to investigate the effectiveness of the proposed algorithm in this section. The root mean squared error (RMSE) of theith source with regard to DOA can be written as(57)RMSEθi=1Mc∑ζMcθ^iζ−θi2,(58)RMSEφi=1Mc∑ζMcφ^iζ−φi2,where θ^ςi and φ^ςi are the estimated nominal azimuth and elevation of the ith source where ς means the estimated value in the ςth Monte Carlo run. Mc is the total number of Monte Carlo simulations. The signal-to-noise ratio (SNR) is defined as 10log(1/σ2) where the noise is assumed to be the Gaussian white zero mean with variance σ2.In the first example, we investigate the performance of the method proposed, considering the underdetermined scenario. Figure4(a) demonstrates that 16 2D CD sources distributed uniformly from diamond area from (30°, 81°), (72°, 40°) to (114°, 81°), (72°, 122°) are set to be estimated. The nested array is set with M1 = 2, M2 = 3, M = 11, and d = λ/2. Angular spreads are 2°, SNR = 5 dB, the number of snapshots equal to 300, and Mc = 100. An indicator reflecting the number of sources detected is needed to measure the performance of different estimators. In the first trial, all 16 sources are Gaussian while in the second trial are uniform. Detection probability is defined as Nd/16 where Nd is the number of sources estimated effectively. Estimation is regarded as effective when the estimated angles satisfy θ^i−θi2+ϕ^i−ϕi20.5≤5∘. Figure 4(b) exhibits the detection probability with change of SNR. As can be seen from Figure 4(b), all the 2D CD sources can be detected when SNR exceeds 5 dB whether sources are Gaussian or uniform. Detection probability decreases as SNR decreases. The difference of detection probability between the two types of CD sources is not significant. This example demonstrated that the nested array and algorithm proposed can resolve 2D CD sources in undermining circumstances.Figure 4
(a) Location of 2D CD sources to be estimated; (b) detection probability versus SNR.
(a)
(b)In the second example, we investigate the performance of the method proposed with regard to angular spreads of sources. Based on the assumption that angular spreads of sources are small rotational invariance relationships that get deduced, the influence of angular spreads on the estimation should be investigated. Both Gaussian and uniform sources are considered, respectively. Experiments are set withM1 = 1, M2 = 1, M = 5, d = λ/2, Mc = 100, SNR = 5 dB, and the number of snapshots equal to 300. Nominal angles of sources are (35°, 35°) and (145°, 145°). Angular spreads of sources are set from 0° to 10°. RMSE is replaced by the mean values of RMSEθi and RMSEφi of two sources. As shown in Figure 5, RMSE reaches 0.58° when angular spreads are 10° for Gaussian sources, while RMSE reaches 0.73° when angular spreads reach 10° for uniform sources. It is obvious that estimated results deteriorate for both kinds of sources as the angular spreads increase. The estimation performance is still satisfactory within 10°. The method proposed exhibits robustness under the condition of small angular spreads. Estimation of the method proposed can be conducted without prior knowledge of ASDF and can also be demonstrated.Figure 5
RMSE estimated for 2D CD sources versus angular spreads.In the third example, we investigate the performance of the method proposed versus SNR and the number of snapshots. Gaussian CD sources with centers (35°, 35°) and (145°, 145°) are set to be estimated. Angular spreads are 2°. Nested array configuration is the same as example 2, Mc = 100. RMSE takes mean values of sources. Experiments are conducted in comparison with existing methods for 2D CD sources based on the L-shape array [12] with 5 sensors, MP using TPLA [15] with 7 sensors, and SOS using DUCA [19] with 6 sensors, respectively. Figure 6(a) shows simulation results as well as the Cramer–Rao bound (CRB) [37] while SNR is ranging from −5 dB to −30 dB with the number of snapshots set at 300. Figure 6(b) shows simulation results as well as CRB with the number of snapshots ranging from 200 to 1000 while SNR is fixed at 5 dB. Figure 6 shows that RMSE of all methods decrease with either SNR increase or number of snapshots increase, and the method proposed has lower RMSE than others. The method proposed illustrates greater performance than the traditional methods with similar sensor number under the same experimental conditions. The reason probably lies in that the nested array possesses larger aperture than uniform arrays with similar number of sensors.Figure 6
(a) Estimation of the method proposed versus SNR; (b) estimation of the method proposed versus number of snapshots.
(a)
(b)In the fourth example, we investigate the performance of the proposed method versus sensor numberM. The sensor number of subarrays A1 and A2 is set according to the optimal configuration described in Table 1. Both Gaussian and uniform sources are considered, respectively. The sources to be estimated in the first trial are Gaussian with nominal angles (35°, 35°) and (145°, 145°) while in the second trial are uniform with the same nominal angles as the first trail. RMSE takes the mean values of sources. Figure 7 shows simulation results with the number of total sensors ranging from 5 to 19 while SNR is fixed at 5 dB and the number of snapshots fixed at 300. Angular spreads are 2°, d = λ/2. It can be concluded that the proposed method provides better performance for both kinds of 2D CD sources as the sensor number increases, but RMSE changes slightly when the total sensor number increases to a certain extent. The reason mainly lies in that the rotational invariance relationship is an approximate resolution under small angular spreads assumption. When the sensor number is bigger, on one hand, the coarray aperture increases which is advantageous to the estimation performance; on the other hand, later sensors in the virtual array is a result of more transitive relationships, which means the proposed method can improve the estimation accuracy to a certain extent.Figure 7
Estimation of the method proposed versus the total sensor numberM.


## 5. Conclusion

Aiming at DOA of 2D CD sources, a 2D nested array as well as an estimation algorithm is proposed in this paper. According to the difference coarray concept, double parallel hole-free virtual uniform linear arrays are obtained by vectorization operation on cross-correlation matrices of subarrays. Configuration of virtual arrays containing close form of sensor coordinates, maximum DOF, and the optimal parameters get deduced. Rotational invariance relationships of virtual arrays are derived under the condition that angular spreads of sources are small. A matrix reconstruction algorithm is detailed based on an ESPRIT-like framework on matrices constructed by extracting and regrouping the receive vectors of the virtual arrays. Numerical simulations considering underdetermined scenario, angular spread, the influence of experiment, and total sensor number are conducted. Outcomes illustrated that the method proposed can deal with more sources than sensors without prior knowledge of ASDF as well as exhibit higher accuracy than uniform arrays under the same experimental conditions.



---

*Source: 6494967-2020-05-12.xml*