# DOA Estimation of Two-Dimensional Coherently Distributed Sources Based on Spatial Smoothing of Uniform Rectangular Arrays

Authors: Tao Wu; Yiwen Li; Xiaofeng Zhang; Yijie Huang; Qingyue Gu; Chaoqi Fu
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

Aiming at the direction-of-arrival (DOA) estimation of two-dimensional (2D) coherently distributed (CD) sources which are coherent with each other, we explore the propagator method based on spatial smoothing of a uniform rectangular array (URA). The rotational invariance relationships with respect to the nominal azimuth and nominal elevation are obtained under the small angular spreads assumption. A propagator operator is constructed through spatial smoothing of sample covariance matrices firstly. Then, combination of propagator and identical matrix is divided according to rotational operators, and the nominal angles can be obtained through eigendecomposition lastly. Realizing angle matching automatically, the proposed method can estimate multiple DOAs of 2D coherent CD sources without spectral peak searching and prior knowledge of deterministic angular signal distribution function. Simulations are conducted to verify the effectiveness of the proposed method.


---

## Body

## 1. Introduction

In the field of array signal processing, traditional DOA estimation is based on point source models. In the real surroundings of radar and sonar systems, because of multipath propagation between receive arrays and targets, especially when the distances of targets and receive arrays are short, the spatial scattering of targets cannot be ignored, and the assumed condition of point source models is no longer valid. In such a condition, DOA estimation based on distributed source models has presented better accuracy [1]. A distributed source can be regarded as an assembly of point sources within a spatial distribution where point sources can be called scatterers.According to coherence of scatterers, distributed sources can be classified as coherently and incoherently distributed sources [1]. Incoherently distributed (ID) sources are defined as scatterers within a source which are uncorrelated. On the contrary, coherently distributed (CD) sources are defined as scatterers within a source which are coherent. According to spatial distribution of sources, distributed sources can be classified as one-dimensional (1D) distributed sources and two-dimensional (2D) distributed sources. Under the assumption that scatterers of a target and receive array are not in the same plane, the 2D distributed sources should be more generally.No matter CD or ID sources, most estimators suppose that signal from different sources is uncorrelated. In the field of underwater detection, the general process of detection is emitting narrowband pulse sound signal firstly, followed by receiving and analyzing the target’s backscatter signal. Thus, different CD sources are supposed to be coherent, which is more reasonable on account of the coherent multipath characteristics of underwater acoustic channel.As to ID sources, utilizing different array configurations representative estimators have been proposed in [2–11]. Supposing that signals from different CD sources are coherent, sample covariance matrices are rank deficient. Thus, subspace-based algorithms [1, 2, 12–19] and ESPRIT class algorithms [20–23] which are based on eigendecomposition of sample covariance matrices cannot be applied for DOA estimation. PM class algorithms [24–27] based on linear operation of full rank sample covariance matrices are no longer applicable too. Considering signals from different CD sources are coherent, authors of [28] have proposed a DOA estimator by means of Toeplitz operation of sample covariance matrices. Authors of [29] have proposed a generalized MUSIC algorithm based on spatial smoothing. Both of the two algorithms above deal with 1D CD sources. Authors of [30] have proposed a general model of 2D CD sources coherent with each other which are defined as 2D coherent CD sources and proposed two spatial smoothing methods with double parallel linear arrays (DPLA).In [1, 6–23], CD sources have been regarded as uncorrelated with each other; estimators cannot be applied for CD sources which are correlated. In [24], CD sources correlated with each other are discussed while sources are supposed to be 1D case. Literatures [10, 31] have considered 2D DOA estimators under uniform rectangular arrays, whereas they deal with ID sources and point sources, respectively. In [30], authors have considered 2D coherent CD sources and DOA estimators based on DPLA. In this paper, we are concerned on DOA estimation for 2D coherent CD sources under uniform rectangular array (URA). 2D coherent CD source and URA are introduced first. Then, the relationships of rotational invariance within and between subarrays have been derived under the assumption of small angular spreads. Afterwards, for decoherence of 2D coherent CD sources, spatial smoothing approach of URA is proposed, and a modified propagator is detailed for DOA estimation of 2D coherent CD sources. The method proposed need not angles matching and prior knowledge of deterministic angular signal distribution function. The simulation outcomes indicate that the proposed method is effective and advanced as for DOA estimation of 2D coherent CD sources.


## 2. Array Configuration and Signal Model

Figure1 shows the URA configuration placed in the xoz plane. Arrays consist of M × K sensors separated by d meters along the direction of the x and z axis. Suppose that there are q 2D CD sources with nominal angles (θi, ϕi) (i=1,2,…,q) impinging the arrays, where θi is nominal azimuth angle and ϕi is nominal elevation angle of ith source, θi ∈ (0, π), ϕi ∈ (0, π). The noise is assumed to be additive Gaussian white with zero mean and uncorrelated between sensors.Figure 1
The uniform rectangular arrays configuration.Denoteymk(t) as the signal received by (m, k)th sensor, which can be expressed as follows:(1)ymkt=∑i=1qvit∑il=1iLαilej2πdm−1cosθilsinϕil/λej2πdk−1cosϕil/λ+nmkt,where vit is the impinging signal to the ith CD source, iL is the number of scatterers of the ith source, and (θil, ϕil) is the azimuth and elevation of the ilth scatterer of the ith source. Suppose different scatterers of the same source differ by one phase delay and a random amplitude gain. αil is complex gain of the ilth scatterer reflecting the reflection coefficient. E(αil) = αi. nmk(t) are noises received.g
i(θ, ϕ; ui) denotes the deterministic angular signal distribution function (ASDF) of the ith source which reflects the distribution of scatterers of a source and is related to the geometry and surface property of the source. A 2D ASDF is generally characterized by the parameter set ui = (θi, ϕi, σθi, σϕi) with four elements denoting the nominal azimuth, nominal elevation, azimuth angular spread, and elevation angular spread, respectively. Nominal azimuth and nominal elevation represent the center of the target which also can be expressed as DOA; azimuth spread and elevation spread indicating 2D spatial extension of a target are collectively called angular spreads.For a Gaussian CD source, which means scatterers of the source obey Gaussian distribution, deterministic ASDF of the source can be expressed as(2)giθ,ϕ;ui=12πσθiσϕiexp−0.5θ−θiσθi2+ϕ−ϕiσϕi2.For a uniform CD source, which means scatterers of the source obey uniform distribution, deterministic ASDF of the source can be expressed as(3)giθ,ϕ;ui=14σθiσϕi,θ−θi≤σθiandϕ−ϕi≤σϕi,0,θ−θi≥σθiorϕ−ϕi≥σϕi.Then, the received signal of (m, k)th sensor can be expressed as (see Appendix A)(4)ymkt=∑i=1qsit∫∫ejπm−12dcosθsinϕ/λejπk−12dcosϕ/λgiθ,ϕ;uidθdϕ+nmkt,where si(t) = vitiLαi is the reflected signal of the ith CD source.Apparently, whenσθi = σϕi = 0, giθ,ϕ;ui can expressed as a 2D Dirac function δθ−θi,ϕ−ϕi. Then, equation (4) can written as(5)ymkt=∑i=1qsitejπm−12dcosθisinϕi/λejπk−12dcosϕi/λ+nmkt,which is equivalent to the received signal model with regard to the point source. Most CD source models assume that scatterers within a source are coherent but scatterers between different sources are uncorrelated. According to characteristics of underwater acoustic channel, the assumption that scatterers between different sources are also coherent is more appropriate. So the reflected signal of the ith CD source can be written as(6)sit=s1tρi,where s1(t) reflects of time behavior of 1st source and ρi is the complex correlation coefficient between the ith source and 1st source. Define the generalized steering coefficient ηmk(θi, ϕi) of (m, k)th sensor as(7)ηmkθi,ϕi=∫∫ej2πdm−1cosθsinϕ/λej2πdk−1cosϕ/λgiθ,ϕ;uidθdϕ.Reflecting the responses of (m, k)th sensor to all CD sources, the generalized steering vector of (m, k)th sensor can be written as(8)amkθ,ϕ=ηmkθ1,ϕ1,ηmkθ2,ϕ2,…,ηmkθq,ϕq.Thus, the received signal of (m, k)th sensor can be written as(9)ymkt=amkθ,ϕs1tϱ+nmkt,where ρ = ρ1,ρ2,…,ρkT is the correlation coefficient vector. When the angular spreads of 2D CD sources are small and d/λ = 1/2, the following rotational invariance relationships can be obtained (see Appendix B):(10)ηmkθi,ϕi≈ejπcosθisinϕiηm−1kθi,ϕi,(11)ηmkθi,ϕi≈ejπcosϕiηmk−1θi,ϕi,(12)ηmkθi,ϕi≈ejπcosθisinϕiejπcosϕiηm−1k−1θi,ϕi.


## 3. Proposed Method

### 3.1. Subarray Signal Model

As shown in Figure2, each subarray has M0 and K0 sensors along the direction of the x and z axes, respectively. Thus, subarrays smooth Ms = M − M0 + 1 times along the direction of the x axis and Ks = K − K0 + 1 times along the direction of the z axis.Figure 2
The subarray grouping of URA on thexoz plane.Considering the (1, 1)th subarray, the receive vector ofM0 sensors along the x axis can be expressed as(13)x111t=Ax11s1tϱ+nx111t,where superscript 11 denotes the variable or vector with regard to (1, 1)th subarray. Define x111 array as the closest and parallel to x axis of the (1, 1)th subarray. The array x1mk can be defined similarly as the (m, k)th subarray. nx111t is the noise vector of array x111.Reflecting the responses of thex111 array to all CD sources, Ax11 is the generalized steering vector of x111, which can be defined as(14)Ax11=a11θ,ϕa21θ,ϕ⋯aM01θ,ϕ=η11θ1,ϕ1η11θ2,ϕ2⋯η11θq,ϕqη21θ1,ϕ1η21θ2,ϕ2⋯η21θ2,ϕ2⋯⋯⋯⋯ηM01θ1,ϕ1ηM01θ1,ϕ1⋯ηM01θq,ϕq.According to equations (10)–(12), the received vector of closest kth and parallel to x axis array xk11 in the (1, 1)th subarray can be expressed as(15)xk11t≈Ax11Φzk−1s1tϱ+nxk11t,where Φz is the rotation operator, which can be written as(16)Φz=diagej2πdcosϕ1/λ,ej2πdcosϕ2/λ,…,ej2πdcosϕq/λ.The received vector of the (1, 1)th subarray along thex axis can be expressed as(17)X11t=x111tx211t⋯xK011t≈Ax11Ax11Φz⋯Ax11ΦzK0−1s1tϱ+nX11t,where nX11t is the noise vector of the (1, 1)th subarray along x axis, which can be expressed as(18)nX11t=nx111tT,nx211tT,…,nxK011tTT.From equations (10)–(12) we can conclude that the received vector of (m, k)th subarray along x axis can be obtained as(19)Xmkt=x1mktx2mkt⋯xK0mkt≈Ax11Ax11Φz⋯Ax11ΦzK0−1Φxm−1Φzk−1s1tϱ+nXmkt.Considering the (1, 1)th subarray and definingz111 array as the closest and parallel to z axis of the (1, 1)th subarray, the receive vector of K0 sensors along the z axis can be expressed as(20)z111t=Az11s1tϱ+nz111t,where Az11 is the generalized steering vector of z111, which can be defined as(21)Az11=a11θ,ϕa12θ,ϕ⋯a1Koθ,ϕ=η11θ1,ϕ1η11θ2,ϕ2⋯η11θq,ϕqη12θ1,ϕ1η12θ2,ϕ2⋯η12θ2,ϕ2⋯⋯⋯⋯η1K0θ1,ϕ1η1K0θ1,ϕ1⋯η1K0θq,ϕq.According to equations (10)–(12), the received vector of closest mth and parallel to z axis array zm11 can be expressed as(22)zm11t≈Az11Φxm−1s1tϱ+nzm11t,where Φx is the rotation operator, which can be written as(23)Φx=diagej2πdcosθ1sinϕ1/λ,ej2πdcosθ2sinϕ2/λ,…,ej2πdcosθqsinϕq/λ.The received vector of the (1, 1)th subarray along thez axis can be expressed as(24)Z11t=z111tz211t⋯zM011t≈Az11Az11Φx⋯Az11ΦxM0−1s1tϱ+nZ11t,where nZ11t is the noise vector of the (1, 1)th subarray along z axis, which can be expressed as(25)nZ11t=nz111tT,nz211tT,…,nzK011tTT.From equations (10)–(12), the received vector of (m, k)th subarray along x axis can be obtained as follows:(26)Zmkt=z1mktz2mkt⋯zM0mkt≈Az11Az11Φx⋯Az11ΦxM0−1Φxm−1Φzk−1s1tϱ+nZmkt.DefineX111t as a vector selecting from 1 to M0K0 − M0 row of X11t. X211t is a vector selecting from M0 + 1 to M0K0 row of X11t. X111t and X211t can be expressed as(27)X111t=X11t1:,MoKo−Mo:=Bx11s1tϱ+nX111t,X211t=X11tM0+1:,M0K0:≈Bx11Φzs1tϱ+nX211t,where X11t (1:, M0K0 − M0:) denotes vectors from 1 to M0K0 − M0 row of X11t. nX111t is noise vector selecting from 1 to M0K0 − M0 row of nX11t. nX211t is noise vector selecting from M0 + 1 to M0K0 row of nX11t. The generalized steering vector of X111t can be expressed as(28)Bx11=Ax11Ax11Φz⋯Ax11ΦzK0−2.Considering the (m, k)th subarray, X1mkt and X2mkt are the sub vectors of Xmkt, which can be expressed as follows:(29)X1mkt=Xmkt1:,M0K0−M0:≈Bx11Φxm−1Φzk−1s1tϱ+nX1mkt,X2mkt=XmktM0+1:,M0K0:≈Bx11Φxm−1Φzk−1Φzs1tϱ+nX2mkt.Similarly, sub vectorsZ111t and Z211t of Z11t can be obtained as(30)Z111t=Z11t1:,M0K0−K0:=Bz11s1tϱ+nZ111t,Z211t=Z11tK0+1:,M0K0:≈Bz11Φxs1tϱ+nZ211t.The generalized steering vector ofZ111t can be expressed as(31)Bz11=Az11Az11Φx⋯Az11ΦxM0−2.Considering the (m, k)th subarray, Z1mkt and Z2mkt are the sub vectors of Zmkt, which can be expressed as(32)Z1mkt=Zmkt1:,M0K0−K0:≈Bz11Φxm−1Φzk−1s1tϱ+nZ1mkt,Z2mkt=ZmktK0+1:,M0K0:≈Bz11Φxm−1Φzk−1Φxs1tϱ+nZ2mkt.


### 3.2. Spatial Smoothing and the Propagator Operator

Define the covariance matrices of the received signal of the (1, 1)th subarray as(33)R111=EX111tZ111tH=σ12Bx11ϱϱHBz11H+N1,R211=EX211tZ111tH=σ12Bx11ΦzϱϱHBz11H+N2,R311=EX111tZ211tH=σ12Bx11ΦxHϱϱHBz11H+N3,R411=EX211tZ211tH=σ12Bx11ΦzΦxHϱϱHBz11H+N4,where σ12=Es1t2, (·)H denotes the Hermitian transpose and N1, N2, N3, and N4 are covariance matrices of the noise vector. The spatial smoothing covariance matrices of the received signal can be obtained as follows:(34)R1S=1MsKs∑m=1Ms∑k=1KsR1mk=σ12Bx11RBz11H+N1S,R2S=1MsKs∑m=1Ms∑k=1KsR2mk=σ12Bx11ΦzRBz11H+N2S,R3S=1MsKs∑m=1Ms∑k=1KsR3mk=σ12Bx11ΦxHRBz11H+N3S,R4S=1MsKs∑m=1Ms∑k=1KsR4mk=σ12Bx11ΦzΦxHRBz11H+N4S,where(35)R=1MsKs∑m=1Ms∑k=1KsΦxm−1Φzk−1ϱϱHΦzk−1HΦxm−1H.N
1
S,N2S, N3S, and N4S are spatial smoothing covariance matrices of the noise vector. The author of [31] has proved the spatial smoothing covariance matrices are full rank matrices when M0 > q, K0 > q and MS ≥ q, KS ≥ q.Construct matrixB as(36)B=Bx11Bx11ΦzBx11ΦxHBx11ΦzΦxH=B1B2,where B1 is a q × q dimensional nonsingular matrix and B2 is a (4M0K0 − 4M0 − q) × q matrix. There exists a q × (4M0K0 − 4M0 − q) dimensional propagator operator P satisfying PHB1=B2.Let(37)E=Iq×qPH,where Iq×q is q × q identity matrix. B can be expressed as follows:(38)B=EB1=Iq×qPHB1.Sample covariance matrices withN snapshots of the (m, k)th subarray can be defined as(39)R^1mk=1N∑t=1NX1mktZ1mktH,(40)R^2mk=1N∑t=1NX2mktZ1mktH,(41)R^3mk=1N∑t=1NX1mktZ2mktH,(42)R^4mk=1N∑t=1NX2mktZ2mktH.Thus, spatial smoothing of sample covariance matrices can be obtained as(43)R^1S,R^2S,R^3S,R^4S=1MsKs∑m=1Ms∑k=1KsR^1mk,∑m=1Ms∑k=1KsR^2mk,∑m=1Ms∑k=1KsR^3mk,∑m=1Ms∑k=1KsR^4mk.CombineR^1S,R^2S,R^3S, and R^4S into(44)R^S=R^1SR^2SR^3SR^4S.Divide (4M0K0 − 4M0) × (M0K0 − K0) dimensional matrix R^s into q × (M0K0 − K0) dimensional matrix G and (4M0K0 − 4M0 − q) × (M0K0 − K0) matrix H:(45)R^S=GH.The propagator operatorP can be obtained as(46)P=GH+HH,where (·)+ denotes the pseudoinverse. Divide E into four (M0K0 − M0) × q dimensional matrices E1, E2, E3, and E4.(47)E=IPH=E1E2E3E4=Bx11Bx11ΦzBx11ΦxHBx11ΦzΦxHB1−1.From this equation, the following relationship can be obtained:(48)E1,E2,E3,E4=Bx11B1−1,Bx11ΦzB1−1,Bx11ΦxHB1−1,Bx11ΦzΦxHB1−1.So, we have(49)Ω1=E1E2+E3E4=B1ΦxHB1−1,(50)Ω2=E1E3+E2E4=B1ΦzB1−1.We can obtain eigenvalueμi (i=1,2,…,q) and it’s corresponding eigenvectors ξi of Ω1 by means of eigendecomposition. From equations (49) and (50), we can conclude that the corresponding eigenvalue of Ω1 and Ω2 have the same eigenvector; so the eigenvector of Ω2 can be obtained as follows:(51)νi=1q11×q⋅Ω2ξi./ξi,where (./) denotes element-wise division operation. When d/λ = 1/2, nominal elevation angle ϕi and nominal azimuth angle θi of ith source can be obtained as follows:(52)ϕi=arccosangleνiπ,θi=arccos−angleμiπsinϕi,i=1,2,…,q,where angle⋅ denotes argument of complex variable.


### 3.3. Computational Procedure Complexity Analysis

Now, our algorithm can be summarized as follows:Step 1: compute spatial smoothing of sample covariance matrices using equations (39)–(43).Step 2: estimate the propagator operatorP using equation (46). Construct E and divide E into E1, E2, E3, E4 from equation (47).Step 3: find eigenvalueμi (i=1,2,…,q) and it’s corresponding eigenvectors ξi through eigendecomposition of Ω1.Step 4: obtain eigenvalueνi of Ω2 from equation (50).Step 5: calculate the nominal azimuthθi and nominal elevation ϕi from equation (51).The computation cost of the proposed method mainly consists of three parts: the calculation of the sample covariance matrixONMsKsM02K02 which is the number of snapshots times the number of spatial smoothing multiply multiplication number with regard to calculation of covariance matrices of subarray, calculation of Ω1 and Ω2 where computation cost mainly lies in pseudoinverse operation O[(2M0K0 − 2M0)3], and eigendecomposition of Ω1 and Ω2O(q3). Suppose 2D CD sources are incoherent, which means decoherence by spatial smoothing is unnecessary. Then, computation cost mainly consists of calculation of the sample covariance matrix ONM2K2, pseudoinverse operation O[(2MK − 2M)3], and eigendecomposition O(q3). It can be concluded that using spatial smoothing to realize decoherence of 2D coherent CD sources does not change the magnitude of computational complexity in essence.


## 3.1. Subarray Signal Model

As shown in Figure2, each subarray has M0 and K0 sensors along the direction of the x and z axes, respectively. Thus, subarrays smooth Ms = M − M0 + 1 times along the direction of the x axis and Ks = K − K0 + 1 times along the direction of the z axis.Figure 2
The subarray grouping of URA on thexoz plane.Considering the (1, 1)th subarray, the receive vector ofM0 sensors along the x axis can be expressed as(13)x111t=Ax11s1tϱ+nx111t,where superscript 11 denotes the variable or vector with regard to (1, 1)th subarray. Define x111 array as the closest and parallel to x axis of the (1, 1)th subarray. The array x1mk can be defined similarly as the (m, k)th subarray. nx111t is the noise vector of array x111.Reflecting the responses of thex111 array to all CD sources, Ax11 is the generalized steering vector of x111, which can be defined as(14)Ax11=a11θ,ϕa21θ,ϕ⋯aM01θ,ϕ=η11θ1,ϕ1η11θ2,ϕ2⋯η11θq,ϕqη21θ1,ϕ1η21θ2,ϕ2⋯η21θ2,ϕ2⋯⋯⋯⋯ηM01θ1,ϕ1ηM01θ1,ϕ1⋯ηM01θq,ϕq.According to equations (10)–(12), the received vector of closest kth and parallel to x axis array xk11 in the (1, 1)th subarray can be expressed as(15)xk11t≈Ax11Φzk−1s1tϱ+nxk11t,where Φz is the rotation operator, which can be written as(16)Φz=diagej2πdcosϕ1/λ,ej2πdcosϕ2/λ,…,ej2πdcosϕq/λ.The received vector of the (1, 1)th subarray along thex axis can be expressed as(17)X11t=x111tx211t⋯xK011t≈Ax11Ax11Φz⋯Ax11ΦzK0−1s1tϱ+nX11t,where nX11t is the noise vector of the (1, 1)th subarray along x axis, which can be expressed as(18)nX11t=nx111tT,nx211tT,…,nxK011tTT.From equations (10)–(12) we can conclude that the received vector of (m, k)th subarray along x axis can be obtained as(19)Xmkt=x1mktx2mkt⋯xK0mkt≈Ax11Ax11Φz⋯Ax11ΦzK0−1Φxm−1Φzk−1s1tϱ+nXmkt.Considering the (1, 1)th subarray and definingz111 array as the closest and parallel to z axis of the (1, 1)th subarray, the receive vector of K0 sensors along the z axis can be expressed as(20)z111t=Az11s1tϱ+nz111t,where Az11 is the generalized steering vector of z111, which can be defined as(21)Az11=a11θ,ϕa12θ,ϕ⋯a1Koθ,ϕ=η11θ1,ϕ1η11θ2,ϕ2⋯η11θq,ϕqη12θ1,ϕ1η12θ2,ϕ2⋯η12θ2,ϕ2⋯⋯⋯⋯η1K0θ1,ϕ1η1K0θ1,ϕ1⋯η1K0θq,ϕq.According to equations (10)–(12), the received vector of closest mth and parallel to z axis array zm11 can be expressed as(22)zm11t≈Az11Φxm−1s1tϱ+nzm11t,where Φx is the rotation operator, which can be written as(23)Φx=diagej2πdcosθ1sinϕ1/λ,ej2πdcosθ2sinϕ2/λ,…,ej2πdcosθqsinϕq/λ.The received vector of the (1, 1)th subarray along thez axis can be expressed as(24)Z11t=z111tz211t⋯zM011t≈Az11Az11Φx⋯Az11ΦxM0−1s1tϱ+nZ11t,where nZ11t is the noise vector of the (1, 1)th subarray along z axis, which can be expressed as(25)nZ11t=nz111tT,nz211tT,…,nzK011tTT.From equations (10)–(12), the received vector of (m, k)th subarray along x axis can be obtained as follows:(26)Zmkt=z1mktz2mkt⋯zM0mkt≈Az11Az11Φx⋯Az11ΦxM0−1Φxm−1Φzk−1s1tϱ+nZmkt.DefineX111t as a vector selecting from 1 to M0K0 − M0 row of X11t. X211t is a vector selecting from M0 + 1 to M0K0 row of X11t. X111t and X211t can be expressed as(27)X111t=X11t1:,MoKo−Mo:=Bx11s1tϱ+nX111t,X211t=X11tM0+1:,M0K0:≈Bx11Φzs1tϱ+nX211t,where X11t (1:, M0K0 − M0:) denotes vectors from 1 to M0K0 − M0 row of X11t. nX111t is noise vector selecting from 1 to M0K0 − M0 row of nX11t. nX211t is noise vector selecting from M0 + 1 to M0K0 row of nX11t. The generalized steering vector of X111t can be expressed as(28)Bx11=Ax11Ax11Φz⋯Ax11ΦzK0−2.Considering the (m, k)th subarray, X1mkt and X2mkt are the sub vectors of Xmkt, which can be expressed as follows:(29)X1mkt=Xmkt1:,M0K0−M0:≈Bx11Φxm−1Φzk−1s1tϱ+nX1mkt,X2mkt=XmktM0+1:,M0K0:≈Bx11Φxm−1Φzk−1Φzs1tϱ+nX2mkt.Similarly, sub vectorsZ111t and Z211t of Z11t can be obtained as(30)Z111t=Z11t1:,M0K0−K0:=Bz11s1tϱ+nZ111t,Z211t=Z11tK0+1:,M0K0:≈Bz11Φxs1tϱ+nZ211t.The generalized steering vector ofZ111t can be expressed as(31)Bz11=Az11Az11Φx⋯Az11ΦxM0−2.Considering the (m, k)th subarray, Z1mkt and Z2mkt are the sub vectors of Zmkt, which can be expressed as(32)Z1mkt=Zmkt1:,M0K0−K0:≈Bz11Φxm−1Φzk−1s1tϱ+nZ1mkt,Z2mkt=ZmktK0+1:,M0K0:≈Bz11Φxm−1Φzk−1Φxs1tϱ+nZ2mkt.


## 3.2. Spatial Smoothing and the Propagator Operator

Define the covariance matrices of the received signal of the (1, 1)th subarray as(33)R111=EX111tZ111tH=σ12Bx11ϱϱHBz11H+N1,R211=EX211tZ111tH=σ12Bx11ΦzϱϱHBz11H+N2,R311=EX111tZ211tH=σ12Bx11ΦxHϱϱHBz11H+N3,R411=EX211tZ211tH=σ12Bx11ΦzΦxHϱϱHBz11H+N4,where σ12=Es1t2, (·)H denotes the Hermitian transpose and N1, N2, N3, and N4 are covariance matrices of the noise vector. The spatial smoothing covariance matrices of the received signal can be obtained as follows:(34)R1S=1MsKs∑m=1Ms∑k=1KsR1mk=σ12Bx11RBz11H+N1S,R2S=1MsKs∑m=1Ms∑k=1KsR2mk=σ12Bx11ΦzRBz11H+N2S,R3S=1MsKs∑m=1Ms∑k=1KsR3mk=σ12Bx11ΦxHRBz11H+N3S,R4S=1MsKs∑m=1Ms∑k=1KsR4mk=σ12Bx11ΦzΦxHRBz11H+N4S,where(35)R=1MsKs∑m=1Ms∑k=1KsΦxm−1Φzk−1ϱϱHΦzk−1HΦxm−1H.N
1
S,N2S, N3S, and N4S are spatial smoothing covariance matrices of the noise vector. The author of [31] has proved the spatial smoothing covariance matrices are full rank matrices when M0 > q, K0 > q and MS ≥ q, KS ≥ q.Construct matrixB as(36)B=Bx11Bx11ΦzBx11ΦxHBx11ΦzΦxH=B1B2,where B1 is a q × q dimensional nonsingular matrix and B2 is a (4M0K0 − 4M0 − q) × q matrix. There exists a q × (4M0K0 − 4M0 − q) dimensional propagator operator P satisfying PHB1=B2.Let(37)E=Iq×qPH,where Iq×q is q × q identity matrix. B can be expressed as follows:(38)B=EB1=Iq×qPHB1.Sample covariance matrices withN snapshots of the (m, k)th subarray can be defined as(39)R^1mk=1N∑t=1NX1mktZ1mktH,(40)R^2mk=1N∑t=1NX2mktZ1mktH,(41)R^3mk=1N∑t=1NX1mktZ2mktH,(42)R^4mk=1N∑t=1NX2mktZ2mktH.Thus, spatial smoothing of sample covariance matrices can be obtained as(43)R^1S,R^2S,R^3S,R^4S=1MsKs∑m=1Ms∑k=1KsR^1mk,∑m=1Ms∑k=1KsR^2mk,∑m=1Ms∑k=1KsR^3mk,∑m=1Ms∑k=1KsR^4mk.CombineR^1S,R^2S,R^3S, and R^4S into(44)R^S=R^1SR^2SR^3SR^4S.Divide (4M0K0 − 4M0) × (M0K0 − K0) dimensional matrix R^s into q × (M0K0 − K0) dimensional matrix G and (4M0K0 − 4M0 − q) × (M0K0 − K0) matrix H:(45)R^S=GH.The propagator operatorP can be obtained as(46)P=GH+HH,where (·)+ denotes the pseudoinverse. Divide E into four (M0K0 − M0) × q dimensional matrices E1, E2, E3, and E4.(47)E=IPH=E1E2E3E4=Bx11Bx11ΦzBx11ΦxHBx11ΦzΦxHB1−1.From this equation, the following relationship can be obtained:(48)E1,E2,E3,E4=Bx11B1−1,Bx11ΦzB1−1,Bx11ΦxHB1−1,Bx11ΦzΦxHB1−1.So, we have(49)Ω1=E1E2+E3E4=B1ΦxHB1−1,(50)Ω2=E1E3+E2E4=B1ΦzB1−1.We can obtain eigenvalueμi (i=1,2,…,q) and it’s corresponding eigenvectors ξi of Ω1 by means of eigendecomposition. From equations (49) and (50), we can conclude that the corresponding eigenvalue of Ω1 and Ω2 have the same eigenvector; so the eigenvector of Ω2 can be obtained as follows:(51)νi=1q11×q⋅Ω2ξi./ξi,where (./) denotes element-wise division operation. When d/λ = 1/2, nominal elevation angle ϕi and nominal azimuth angle θi of ith source can be obtained as follows:(52)ϕi=arccosangleνiπ,θi=arccos−angleμiπsinϕi,i=1,2,…,q,where angle⋅ denotes argument of complex variable.


## 3.3. Computational Procedure Complexity Analysis

Now, our algorithm can be summarized as follows:Step 1: compute spatial smoothing of sample covariance matrices using equations (39)–(43).Step 2: estimate the propagator operatorP using equation (46). Construct E and divide E into E1, E2, E3, E4 from equation (47).Step 3: find eigenvalueμi (i=1,2,…,q) and it’s corresponding eigenvectors ξi through eigendecomposition of Ω1.Step 4: obtain eigenvalueνi of Ω2 from equation (50).Step 5: calculate the nominal azimuthθi and nominal elevation ϕi from equation (51).The computation cost of the proposed method mainly consists of three parts: the calculation of the sample covariance matrixONMsKsM02K02 which is the number of snapshots times the number of spatial smoothing multiply multiplication number with regard to calculation of covariance matrices of subarray, calculation of Ω1 and Ω2 where computation cost mainly lies in pseudoinverse operation O[(2M0K0 − 2M0)3], and eigendecomposition of Ω1 and Ω2O(q3). Suppose 2D CD sources are incoherent, which means decoherence by spatial smoothing is unnecessary. Then, computation cost mainly consists of calculation of the sample covariance matrix ONM2K2, pseudoinverse operation O[(2MK − 2M)3], and eigendecomposition O(q3). It can be concluded that using spatial smoothing to realize decoherence of 2D coherent CD sources does not change the magnitude of computational complexity in essence.


## 4. Simulation Results

In this section, five simulation experiments are conducted to verify the effectiveness of the algorithm we proposed. All simulation experiments are based on array configuration as shown in Figure1. The distance of adjacent sensor d is set at λ/2.RMSE
i denotes root mean squared error (RMSE) of the ith source, which can be expressed as(53)RMSEi=1Mc∑ςMcθ^iς−θi2+1Mc∑ςMcϕ^iς−ϕi2,where Mc is the Monte Carlo simulations number which is set at 100. θ^iς and ϕ^iς are the estimated nominal azimuth and nominal elevation of ith source in ςth Monte Carlo simulation.In the first example, we investigate the performance of the proposed method versus three 2D CD sources which are uncorrelated with each other. The first source is Gaussian CD source with parameter set (30°, 45°, 2°, 2°). The second and third sources are uniform with parameter sets (50°, 45°, 2°, 2°) and (50°, 60°, 2°, 2°). The number of snapshots is set at 200. The subarray number parameter isMs = Ks = 4, subarray elements parameter is M0 = K0 = 4. RMSE is the mean RMSEi of three sources. As shown in Figure 3, the proposed method presents better estimation as the SNR increases. Figure 3 also shows estimation of the DSPE [18] which uses L-shaped arrays containing sensors in x and z axes, and the ESPRIT [23] method uses DPLA. It can be concluded that DOA estimation of all the three methods is effective with respect to 2D CD sources which are uncorrelated with each other.Figure 3
RMSE estimated by three methods for 2D CD sources uncorrelated with each other.In the second example, we investigate the estimation of 2D coherent CD sources versus SNR and number of snapshots. We consider three full coherence sources with equal power. The first source is Gaussian with parameter set [30°, 45°, 2°, 2°]. The second and third sources are uniform with parameter sets [50°, 45°, 2°, 2°] and [50°, 60°, 2°, 2°].Ms = Ks = 4 and M0 = K0 = 4. RMSE takes mean values of the three sources. Figure 4(a) shows RMSE with SNR varying from 0 dB to 30 dB while the number of snapshots is set at 200. Figure 4(b) shows RMSE with the number of snapshots ranging from 100 to 1000 while SNR is fixed at 15 dB. It can be observed that the proposed algorithm performs better as SNR or number of snapshots increases. Nevertheless, the DSPE [18] and ESPRIT [23] present big errors as shown in Figures 4(a) and 4(b). It can be concluded that utilizing traditional methods for 2D coherent CD sources is invalid.Figure 4
(a) RMSE estimated by three algorithms versus SNR; (b) RMSE estimated by three algorithms versus number of snapshots.
(a)
(b)In the third example, we investigate the performance of the proposed method versus angular spreads. We consider two scenarios. One scenario has three full coherence Gaussian sources with equal power, and the nominal angles are (30°, 45°), (50°, 45°), and (50°, 60°). The other has three uniform sources with the same nominal angles as the first scenario. To simplify the analysis, we assume that azimuth angular spread is equal to elevation angular spread within each source, and angular spreads of three sources are the same;σ is used to replace σθi and σϕi for convenience. The number of snapshots is set at 200 and SNR is 15 dB. Ms = Ks = 4 and M0 = K0 = 4. RMSE of the two trails is defined as the mean values of three sources. From Figure 5, it can be observed that estimation errors increase as angular spreads increase with respect to both Gaussian and uniform sources. RMSE reaches 0.027° as angular spreads increase to 5° and reaches 0.4° as angular spreads increase to 10°, which still is a satisfactory result. It can be seen that the proposed method has a good performance under the small angular spreads condition. It is noteworthy that when the angular spreads is 0°, 2D CD sources degenerate into point sources; this experiment also manifests that the method proposed can also be applied to point sources.Figure 5
RMSE estimated by the proposed versus the angular spread for 2D coherent CD sources.In the fourth example, we investigate the performance of the proposed method versus subarray parameters. We consider three full coherence sources with equal power. The source parameter sets are the same as the second example. The number of snapshots is set at 200, and SNR is 15 dB. For the convenience of analysis,Ms is supposed equal to Ks, and M0 is supposed equal to K0. Figure 6 displays the performance of the proposed method versus the subarray number parameter Ms with M0 = 4, while Figure 7 shows the performance versus subarray elements parameter M0 with Ms = 4. From Figures 6 and 7, it can be concluded that the proposed method can estimate DOA of the 2D coherent CD sources effectively when M0 > q, K0 > q, MS ≥ q, and KS ≥ q and provides better performance as subarray number or subarray elements increases, but RMSE changes slightly when M0 or Ms increases to a certain extent.Figure 6
RMSE estimated by the proposed versusMs with M0 = 4.Figure 7
RMSE estimated by the proposed versusM0 with Ms = 4.In the fifth example, we investigate the performance of the proposed method near the boundary region in comparison with the method [30] utilizing DPLA. We consider double full coherence Gaussian sources with equal power. The number of snapshots is set at 200, and SNR is 15 dB. Angular spreads are set at 2°. In the first trail, nominal elevations of two sources are set at 30° and 50°, and nominal azimuths of two sources are varying from 0° to 180° simultaneously. Similarly, in the second trail, nominal azimuths of two sources are set at 30° and 50°, and nominal elevations of two sources are varying from 0° to 180° simultaneously. Ms = Ks = 3 and M0 = K0 = 3, and the total sensor number of URA is 25. DPLA is also set with total sensor number as 25. RMSE takes the mean values of the two sources. Figures 8(a) and 8(b) show the first and second trails, respectively. As shown in Figures 8(a) and 8(b), the method proposed has a better performance near the boundary regions and more significant near the elevation boundary region. Aperture sizes of URA are balanced in azimuth and elevation dimensions. Compared with URA, DPLA has smaller aperture in the elevation dimension. Therefore, errors of elevation estimation near the elevation boundary region are greater; consequently, estimation deterioration will be more significant in these areas as solving elevation and azimuth of both methods are coupled.Figure 8
(a) RMSE estimated versus azimuth boundary region; (b) RMSE estimated versus elevation boundary region.
(a)
(b)


## 5. Conclusions

In this study, we have considered the problem of estimation of 2D coherent CD sources utilizing URA. The rotation invariance relationships within and between subarrays have been deduced under the small angular spreads assumption. Decoherence can be realized by virtue of 2D spatial smoothing of URA. Propagator method of sample covariance matrices based on spatial smoothing has been introduced in detail. Simulation investigates experiment conditions containing SNR, number of snapshots, angular spreads, and subarray parameters. Simulation outcomes indicated that the proposed method is effective for DOA estimation of 2D coherent CD sources and has a better performance near the boundary regions. The method proposed is also effective for the 2D CD sources which are incoherent as well as point sources.



---

*Source: 4816183-2019-10-07.xml*