# Two-Dimensional DOA Estimation for Coprime Planar Arrays Based on Self-Correlation Tensor

Authors: Hao Li; Weijia Cui; Chunxiao Jian; Haiyun Xu; Fengtong Mei
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In the coprime planar array (CPA), the existing tensor DOA estimation has the problem that the statistics are not fully utilized. We propose a two-dimensional DOA estimation method based on tensor self-correlation, which realizes the high-resolution and high-precision joint estimation of elevation angle and azimuth angle. Firstly, we represent the received signals of two subarrays with tensors and then obtain the self-correlation covariance tensor of the subarrays themselves and the cross-correlation covariance tensor of the two subarrays. Then, we extract the covariance tensor corresponding to the maximum continuous virtual array and prove the expression of the maximum continuous virtual array aperture of the proposed method. Compared with the existing methods, the proposed method effectively improves the maximum aperture of the continuous virtual array. Finally, the signal subspace is solved by tensor expansion and tensor decomposition. Simulation results show that under the same conditions, the proposed method has higher estimation accuracy and degree of freedom than the cross-correlation tensor method, and the resolution is also improved significantly.


---

## Body

## 1. Introduction

Direction of arrival (DOA) is a technology that processes and analyzes the data received from sensor array to obtain target location information, which is an important research topic in array signal processing [1]. It has been widely applied in radar, sonar, wireless communication [2–6], and other areas. Most subspace DOA methods were initially proposed based on non-sparse uniform linear arrays [7–10], which can complete DOA estimation in the case of overdetermination, that is, M array sensors can complete M−1 signal estimation at most [11–13]. However, in the actual situation, the number of located target sources is greater than the number of array sensors, that is, underdetermined estimate.To improve the degree of freedom of the array and achieve a practical estimation of multiple signal sources under unknown conditions, experts and scholars have proposed sparse array [14–17], which can effectively expand the aperture of the array and improve the degrees of freedom and estimation accuracy, and typical sparse arrays include nested arrays [18, 19], coprime arrays [20], and minimum redundant arrays (MRAs) [21]. Also, many algorithms are proposed based on these sparse arrays, such as spatial smoothing MUSIC (SSM) [22], compressive sensing (CS) [23], and discrete Fourier transform (DFT) [24]. However, the above arrays are all one-dimensional arrays. Usually, one-dimensional arrays can only estimate one-dimensional DOA. To achieve two-dimensional DOA estimation, many experts and scholars have generalized traditional subspace methods based on one-dimensional arrays to two-dimensional arrays. Aiming to realize DOA estimation with high degrees of freedom as much as possible, two-dimensional MUSIC methods (2-D MUSIC), two-dimensional ESPRIT methods (2-D ESPRIT) [25], and PM methods [26] based on coprime planar arrays and coprime L-arrays have also been proposed, which combine the advantages of sparse formations and realize two-dimensional DOA estimation. However, these methods are based on matrix operation processing, which will blur the data structure and bring loss to the accuracy of the DOA estimation. This operation will blur the data structure and bring loss to the accuracy of the DOA estimation.In this regard, experts and scholars have introduced tensors into DOA estimation [27–32]. Tensors [33–35] are high-dimensional data, which can save information of different physical meanings in other dimensions. By comparing the performance of tensor-based DOA estimation method with matrix-based DOA estimation method, the authors in [36–38] showed that the former has better estimation accuracy. However, these methods are mainly based on the cross-correlation tensor (CCT) for processing the subarray received signal, which fails to make full use of the position information of the array elements, which will lead to the narrowing of the aperture of the virtual array, resulting in loss of degrees of freedom and estimation accuracy.In this paper, we propose a self-correlation tensor-based method for estimating 2-D DOA in the coprime planar array. Compared with the cross-correlation tensor-based method, the proposed method can achieve a higher degree of freedom and higher accuracy of joint estimation of elevation angle and azimuth angle. The proposed method uses tensors to represent the received signal. Through the self-correlation processing of the signal, we can obtain the covariance tensor of the signal and then extract the covariance tensor to get the covariance tensor of the corresponding maximum continuous virtual uniform rectangular array. Finally, we can estimate the azimuth angle and elevation angle of the signal by CP decomposition (CANDECOMP/PARAFAC decomposition (CPD)) [33, 39]. Compared with CCT and SSM, simulation results show that the proposed method improves estimation accuracy and degree of freedom.The main contribution of this paper is summarized as follows:(1)
We propose a new method for estimating the two-dimensional (2-D) DOA in the coprime planar array based on SCT, which achieves the high-resolution and high-precision joint estimation of elevation angle and azimuth angle. It uses tensors to model signals, taking into account the dimensional structure of the data. Compared with the CCT and SSM methods, the estimation accuracy of the proposed method is better than that of CCT and SSM. Compared with the PSS methods, when the number of sources is less than the number of minimum subarray elements, the performance of the proposed method is close to that of the PSS method, and when the number of sources is more than the number of small subarray elements, the proposed method is much better than the PSS method.(2)
Aiming at the problem that the virtual array is not constructed sufficiently in the existing tensor DOA estimation method, the self-covariance information and the cross-covariance information are combined to increase the maximum degree of freedom and estimation accuracy of the DOA estimate. However, in order to apply this method to coprime planar arrays of different size, the maximum continuous virtual planar array aperture expression must be obtained, so we derive the expression of the largest continuous virtual planar array in the appendix, which makes this method more practical. By proof derivation, the aperture expression of the maximum continuous virtual URA ofM2+N2−1 physical array is W=W1∪W2, where W1=x,y|1−M≤x≤M+N−1,1−M−N≤y≤M−1, W1=−W2.The rest of the paper is organized as follows. We give the signal model in Section2. In Section 3, we introduce the two-dimensional DOA estimation based on the proposed method and analyze the performance of the method. The simulation results are provided in Section 4, and we conclude this paper in Section 5.Notations. Scalars are denoted by italic letters (M, m), matrices and column vectors by bold letters C,c, tensors by ℛ,N, and collections by hollow letters ℕ,W. Let •T, •∗, •H, and E• correspond to transpose, conjugate, Hermitian transpose, and statistical expectation. The Kronecker product of two matrices A,B is denoted by A⊗B, and the outer product of two tensors A,ℬ is denoted by A°ℬ.


## 2. Preliminaries

### 2.1. Signal Model

The model of the coprime planar array shown in Figure1 comprises two sparse uniform planar arrays, named as D1 and D2, and the spacing between adjacent elements is Nλ/2, Mλ/2, where N and M are coprime integers and λ is the wavelength of impinging signals. D1 and D2 only coincide at the origin, the total number of physical elements is M2+N2−1, and the reference position of the ith element to the origin is xi,yi. The mathematical expression of the position of the whole array element is Ω=mNd0,mNd0∪nMd0,nMd0|m∈0,M−1,n∈0,N−1. Suppose there are k far-field narrowband independent signals impinging on the array with the signals from θ1,φ1,θ2,φ2,…,θK,φK, where θi∈−π,π, φi∈0,π/2 are the azimuth angle and elevation angle of the ith signal.Figure 1
Geometry of CPA whenM  = 3 and N  = 2.According to the geometric relation of the element and the tensor theory [27], the received signals of D1 can be represented as a 3-D tensor:(1)X1=∑k=1Ka1xuk°a1yvk°Sk+N1,where a1xuk=1,ej2π/λNd0uk,…,ej2π/λM−1Nd0uk, a1yvk=1,ej2π/λNd0vk,…,ej2π/λM−1Nd0vk are steering vectors of D1 and uk=cosθkcosφk, vk=sinθkcosφk are directional factors in the x and y directions. The signal data vector(2)Sk=Sk1,Sk2,…,SkL,where L denotes number of snapshots and N1 is additive white Gaussian noise of D1. Similarly, the expression of D2 received signal can be represent as(3)X2=∑k=1Ka2xuk°a2yvk°Sk+N2,where a2xuk=1,ej2π/λMd0uk,…,ej2π/λN−1Md0uk, a2yvk=1,ej2π/λMd0vk,…,ej2π/λN−1Md0vk are steering vectors of D2 and N2 is additive white Gaussian noise.Then, we can slice the tensorsXi along the snap dimension [40] and get Xi1,Xi2,…,Xil,…,XiL, i=1,2. The 4-D covariance tensor RXiXj can be expressed as(4)ℛXiXj=EXil°Xj∗l=∑k=1Kσk2aix=uk°aiyvk°ajx∗uk°ajy∗vk+Nij=PXiXj+Nij,where i=1,2, j=1,2, σk2=Eskl°sk∗l is the power of kth signal, Ni,j=EN1°N2∗.However, the number of snapshots is limited. Here, we use statistical average instead of time average, and 4-D tensor of the received signal covariance can be represented as(5)ℛ^XiXj=1L∑l=1LXil°Xj∗l,and obtain P^XiXj.


## 2.1. Signal Model

The model of the coprime planar array shown in Figure1 comprises two sparse uniform planar arrays, named as D1 and D2, and the spacing between adjacent elements is Nλ/2, Mλ/2, where N and M are coprime integers and λ is the wavelength of impinging signals. D1 and D2 only coincide at the origin, the total number of physical elements is M2+N2−1, and the reference position of the ith element to the origin is xi,yi. The mathematical expression of the position of the whole array element is Ω=mNd0,mNd0∪nMd0,nMd0|m∈0,M−1,n∈0,N−1. Suppose there are k far-field narrowband independent signals impinging on the array with the signals from θ1,φ1,θ2,φ2,…,θK,φK, where θi∈−π,π, φi∈0,π/2 are the azimuth angle and elevation angle of the ith signal.Figure 1
Geometry of CPA whenM  = 3 and N  = 2.According to the geometric relation of the element and the tensor theory [27], the received signals of D1 can be represented as a 3-D tensor:(1)X1=∑k=1Ka1xuk°a1yvk°Sk+N1,where a1xuk=1,ej2π/λNd0uk,…,ej2π/λM−1Nd0uk, a1yvk=1,ej2π/λNd0vk,…,ej2π/λM−1Nd0vk are steering vectors of D1 and uk=cosθkcosφk, vk=sinθkcosφk are directional factors in the x and y directions. The signal data vector(2)Sk=Sk1,Sk2,…,SkL,where L denotes number of snapshots and N1 is additive white Gaussian noise of D1. Similarly, the expression of D2 received signal can be represent as(3)X2=∑k=1Ka2xuk°a2yvk°Sk+N2,where a2xuk=1,ej2π/λMd0uk,…,ej2π/λN−1Md0uk, a2yvk=1,ej2π/λMd0vk,…,ej2π/λN−1Md0vk are steering vectors of D2 and N2 is additive white Gaussian noise.Then, we can slice the tensorsXi along the snap dimension [40] and get Xi1,Xi2,…,Xil,…,XiL, i=1,2. The 4-D covariance tensor RXiXj can be expressed as(4)ℛXiXj=EXil°Xj∗l=∑k=1Kσk2aix=uk°aiyvk°ajx∗uk°ajy∗vk+Nij=PXiXj+Nij,where i=1,2, j=1,2, σk2=Eskl°sk∗l is the power of kth signal, Ni,j=EN1°N2∗.However, the number of snapshots is limited. Here, we use statistical average instead of time average, and 4-D tensor of the received signal covariance can be represented as(5)ℛ^XiXj=1L∑l=1LXil°Xj∗l,and obtain P^XiXj.


## 3. Tensor DOA Estimation

In this section, we propose a new tensor-based two-dimensional DOA estimation method that effectively incorporates auto-related information of signals. Compared with the methods in [38], the proposed method in this paper effectively increases the aperture of the virtual array and improves the estimation accuracy.


### 3.1. Construction of Virtual Array

We define the dimensional setsℚ1=1,3, ℚ2=2,4 and then use tensor expansion [33, 41] to merge dimensions with the same directional factor:(6)Ui,j=ℛXiXj1,32,4=∑k=1Kσk2ajx∗uk⊗aixuk°ajy∗vk⊗aiyvk+Nij,where i=1,2, j=1,2, the terms ajx∗uk⊗aixuk and ajy∗vk⊗aiyvk deduce a difference coarray set named as ℕ, and ℕ=x,y|x∈ℤi,j,y∈ℤi,j; ℤi,j can be represented as(7)ℤ1,1=m1N−m2Nd0|0≤mi≤M−1,i=1,2ℤ1,2=m1N−n1Md0|0≤n1≤N−1,0≤m1≤M−1ℤ2,1=n1M−m1Nd0|0≤n1≤N−1,0≤m1≤M−1ℤ2,2=n1M−n2Md0|0≤ni≤N−1,i=1,2..Figure2 shows the virtual array position of a coprime planar array with M  = 3 and N  = 2. Consider that coarray tensor Ui is the single snapshot, which cannot meet the uniqueness principle of tensor decomposition, and the number of sources that can be estimated is limited. To solve this problem, we generate new structural coarray tensors ℋi by structural recombination of coarray tensors Ui. However, only continuous virtual planar arrays can be used for structural recombination of coarray tensors Ui. Thus, we divide the continuous sensors in the virtual array into four parts, which are represented as(8)U˜1=x1,y1x1=k1d0,y1=k2d0,1−M−N≤k1≤M−1,1−M≤k2≤M+N−1,U˜2=x2,y2x2=k3d0,y2=k4d0,1−N≤k3≤M−1,1−N≤k4≤M−1,U˜3=−U˜1,U˜4=−U˜2.Figure 2
The location of virtual sensors and real sensors whenM  = 3 and N  = 2.In Figure2, we show four continuous virtual area arrays wrapped in rectangular boxes, respectively. The elements in Ui,j are reorganized to map the augmented virtual URA U˜1, U˜2, U˜3, and U˜4, where the equivalent 2-D coarray tensor U1,U2,U3,U4 can be represented as(9)Ui=∑k=1Kσk2bixuk°biyvk+Nij,where b1x=ejπ1−M−Nμk,ejπ2−M−Nμk,…,ejπM−1μk, b1y=ejπ1−Mμk,ejπ2−Mμk,…,ejπN+M−1μk, b2x=ejπ1−Nμk,ejπ2−Nμk,…,ejπM−1μk,b2y=ejπ1−Nμk,ejπ2−Nμk,…,ejπM−1μk.


### 3.2. Tensor Decomposition of Structure Coarray

Coarray tensorUi is the single snapshot, which cannot meet the uniqueness principle of tensor decomposition, the identification of tensor decomposition is limited, and only a single source can be detected. In order to estimate more sources, the idea of spatial smoothing algorithm can be introduced.We divide the virtualUi∈ℂni×ni into Ji×Ji overlapping subarrays with P×P sensors, which is denoted as Gx,yi, where 1≤x≤P,1≤y≤P, n1=n3=2M+N−1, n2=n4=M+N−1, and Ji+P=ni. The position of the first sensors of each subarray can be represented as −M−N+x,M+N−y, so the coarray tensor of the first subarray can be represented as(10)ℋ1i=Ei∑k=1Kσk2cxuk°cyvk+Nij,where E=1,ejπMuk−Nvk,ejπNuk−Mvk,ejπNuk−Nvk is used to compensate for the phase difference caused by the different positions of Gx,yi in the four virtual URA, cxuk=e−jπM+N−1uk,e−jπM+N−2uk,…,e−jπM+N−Puk, and cyvk=ejπM+N−1vk,ejπM+N−2vk,…,ejπM+N−Pvk are the steering vectors of ℋ1i. Therefore, by space smoothing for virtual URA u˜i, smooth coarray tensor ℋi can be represented as(11)ℋi=Ei∑k=1K∑a=1ℒi∑b=1ℒiσk2cxuk×dixauk°cyvk×diybvk+Nij,where dixauk=ejaπuk,a∈1,2,…,Ji, and diybuk=e−jbπvk,b∈1,2,…,Ji serve as the shifting factor vectors along the direction of x-axis and y-axis, respectively. We can concatenate ℋ1, ℋ2, ℋ3, and ℋ4 along the third dimension as 3-D tensor T∈ℂP×P×4, which is represented as(12)T=ℋ1°w+Nij,and then we do a tensor decomposition of T; we define the ranka as Kruskal’s rank of a. According to the tensor decomposition uniqueness criterion, we have(13)rankCx+rankCy+rankw≥2K+2,where Cx=cxu1,cxu2,…,cxuK, Cy=cyv1,cyv2,…,cyvK, w=w1,w2,w3,w4. It is noted that rankCx=minK,fP,M,N, rankCy=minK,fP,M,N, and fP,M,N=minP2,2×M+N−P2+2×M+N−P2. In order to make full use of the DOFs of the structured coarray tensor T, we should maximize fP,M,N. Take M  = 2, N  = 3 for example; when P=4, degree of freedom (DOF) can reach the maximum 15, and the signal subspace Us can be represented as(14)Us=orthcxu1⊗cyv1,cxu2⊗cyv2,⋯,cxuk⊗cyvk,Then, the noise subspace can be represented as(15)UnUnH=I−UsUsH,where I is the identity matrix. With the steering function of the virtual ℕ for the scanning bin u,v defined as Au,v=cxu⊗cyv, the coarray tensor MUSIC spectrum of MUSIC spectrum based on SCT can be calculated as(16)Pu,v=1Au,vHUnUnHAu,v.As such, the 2-D DOA estimatesθk,φk can be obtained via spectrum searching on Pu,v. The detailed steps of the method are shown in Table 1.Table 1
The steps of the proposed method.
Step 1Construct 3-D tensor of the received signalXiStep 2Construct 4-D self-correlation coarray tensorℛ^XiXj from the received signal Xi of the subarrayStep 3Transform the dimension order of the 4-D tensorℛ^XiXj and then merge the dimensions with the same directional factor to get Ui,jStep 4Take the coarray tensorU1,U2,U3,U4 corresponding to the largest continuous virtual array from Ui,jStep 5Transform the dimension of the coarray tensorUi by tensor expansion to obtain structured coarray tensor TStep 6Perform CPD on the structured coarray tensorT to obtain matrix factors cxuk, cyvk containing direction information and calculate θk,φk by spectral peak search


### 3.3. Complexity Analysis

We analyze the computational complexity of the proposed method in the coprime planar array and compare it with the CCT method, SSM method, and PSS method. Let us first analyze the computational complexity of the proposed method. For a coprime planar array, it can be divided into two uniform planar subarrays of sizesM×M and N×N, the proposed method constructs covariance tensors by received signals, and two subarrays can construct four covariance tensors, so the computational complexity of this step is OLM4+N4+2M2N2. Then, by tensoring the covariance tensor, the structure covariance matrix T is obtained, and the computational complexity of this step is O2KP22M+N−P2. Then, we do a tensor decomposition of T, the computational complexity of this step is 16KP22M+N−P2+3K3, and the complexity required for the spectral peak search to calculate the elevation angle and azimuth angle is OngnlP4, where ng=360°/0.01°, nl=90°/0.01°, and ng×nl represents the number of meshes searched by the spectral peak. Therefore, the total complexity of the SCT method is OLM4+N4+2M2N2+2KP22M+N−P2+16KP22M+N−P2+3K3+ngnlP4. The complexity of the PSS method is OngnlN4/M2+M4/N2. The complexity of the SSM method is OM2+N22L+2KP22M+N−P2+P6+ngnlP4. The complexity of the CCT method is OLKM2N2+2KP2M+N−P2+16KP2M+N−P2+3K3+ngnlP4.For a more intuitive representation, Table2 summarizes the complexity of the different methods. As can be seen from Table 2, the computational complexity of the four methods is mainly derived from the spectral peak search, and the computational complexities of the SCT, CCT, and SSM methods are roughly equal and higher than the PSS method.Table 2
Comparison of computational complexity of different methods.
SCTOLM4+N4+2M2N2+2KP22M+N−P2+16KP22M+N−P2+3K3+ngnlP4CCTOLKM2N2+2KP2M+N−P2+16KP2M+N−P2+3K3+ngnlP4SSMOM2+N22L+2KP22M+N−P2+P6+ngnlP4PSSOngnlN4/M2+M4/N2


## 3.1. Construction of Virtual Array

We define the dimensional setsℚ1=1,3, ℚ2=2,4 and then use tensor expansion [33, 41] to merge dimensions with the same directional factor:(6)Ui,j=ℛXiXj1,32,4=∑k=1Kσk2ajx∗uk⊗aixuk°ajy∗vk⊗aiyvk+Nij,where i=1,2, j=1,2, the terms ajx∗uk⊗aixuk and ajy∗vk⊗aiyvk deduce a difference coarray set named as ℕ, and ℕ=x,y|x∈ℤi,j,y∈ℤi,j; ℤi,j can be represented as(7)ℤ1,1=m1N−m2Nd0|0≤mi≤M−1,i=1,2ℤ1,2=m1N−n1Md0|0≤n1≤N−1,0≤m1≤M−1ℤ2,1=n1M−m1Nd0|0≤n1≤N−1,0≤m1≤M−1ℤ2,2=n1M−n2Md0|0≤ni≤N−1,i=1,2..Figure2 shows the virtual array position of a coprime planar array with M  = 3 and N  = 2. Consider that coarray tensor Ui is the single snapshot, which cannot meet the uniqueness principle of tensor decomposition, and the number of sources that can be estimated is limited. To solve this problem, we generate new structural coarray tensors ℋi by structural recombination of coarray tensors Ui. However, only continuous virtual planar arrays can be used for structural recombination of coarray tensors Ui. Thus, we divide the continuous sensors in the virtual array into four parts, which are represented as(8)U˜1=x1,y1x1=k1d0,y1=k2d0,1−M−N≤k1≤M−1,1−M≤k2≤M+N−1,U˜2=x2,y2x2=k3d0,y2=k4d0,1−N≤k3≤M−1,1−N≤k4≤M−1,U˜3=−U˜1,U˜4=−U˜2.Figure 2
The location of virtual sensors and real sensors whenM  = 3 and N  = 2.In Figure2, we show four continuous virtual area arrays wrapped in rectangular boxes, respectively. The elements in Ui,j are reorganized to map the augmented virtual URA U˜1, U˜2, U˜3, and U˜4, where the equivalent 2-D coarray tensor U1,U2,U3,U4 can be represented as(9)Ui=∑k=1Kσk2bixuk°biyvk+Nij,where b1x=ejπ1−M−Nμk,ejπ2−M−Nμk,…,ejπM−1μk, b1y=ejπ1−Mμk,ejπ2−Mμk,…,ejπN+M−1μk, b2x=ejπ1−Nμk,ejπ2−Nμk,…,ejπM−1μk,b2y=ejπ1−Nμk,ejπ2−Nμk,…,ejπM−1μk.


## 3.2. Tensor Decomposition of Structure Coarray

Coarray tensorUi is the single snapshot, which cannot meet the uniqueness principle of tensor decomposition, the identification of tensor decomposition is limited, and only a single source can be detected. In order to estimate more sources, the idea of spatial smoothing algorithm can be introduced.We divide the virtualUi∈ℂni×ni into Ji×Ji overlapping subarrays with P×P sensors, which is denoted as Gx,yi, where 1≤x≤P,1≤y≤P, n1=n3=2M+N−1, n2=n4=M+N−1, and Ji+P=ni. The position of the first sensors of each subarray can be represented as −M−N+x,M+N−y, so the coarray tensor of the first subarray can be represented as(10)ℋ1i=Ei∑k=1Kσk2cxuk°cyvk+Nij,where E=1,ejπMuk−Nvk,ejπNuk−Mvk,ejπNuk−Nvk is used to compensate for the phase difference caused by the different positions of Gx,yi in the four virtual URA, cxuk=e−jπM+N−1uk,e−jπM+N−2uk,…,e−jπM+N−Puk, and cyvk=ejπM+N−1vk,ejπM+N−2vk,…,ejπM+N−Pvk are the steering vectors of ℋ1i. Therefore, by space smoothing for virtual URA u˜i, smooth coarray tensor ℋi can be represented as(11)ℋi=Ei∑k=1K∑a=1ℒi∑b=1ℒiσk2cxuk×dixauk°cyvk×diybvk+Nij,where dixauk=ejaπuk,a∈1,2,…,Ji, and diybuk=e−jbπvk,b∈1,2,…,Ji serve as the shifting factor vectors along the direction of x-axis and y-axis, respectively. We can concatenate ℋ1, ℋ2, ℋ3, and ℋ4 along the third dimension as 3-D tensor T∈ℂP×P×4, which is represented as(12)T=ℋ1°w+Nij,and then we do a tensor decomposition of T; we define the ranka as Kruskal’s rank of a. According to the tensor decomposition uniqueness criterion, we have(13)rankCx+rankCy+rankw≥2K+2,where Cx=cxu1,cxu2,…,cxuK, Cy=cyv1,cyv2,…,cyvK, w=w1,w2,w3,w4. It is noted that rankCx=minK,fP,M,N, rankCy=minK,fP,M,N, and fP,M,N=minP2,2×M+N−P2+2×M+N−P2. In order to make full use of the DOFs of the structured coarray tensor T, we should maximize fP,M,N. Take M  = 2, N  = 3 for example; when P=4, degree of freedom (DOF) can reach the maximum 15, and the signal subspace Us can be represented as(14)Us=orthcxu1⊗cyv1,cxu2⊗cyv2,⋯,cxuk⊗cyvk,Then, the noise subspace can be represented as(15)UnUnH=I−UsUsH,where I is the identity matrix. With the steering function of the virtual ℕ for the scanning bin u,v defined as Au,v=cxu⊗cyv, the coarray tensor MUSIC spectrum of MUSIC spectrum based on SCT can be calculated as(16)Pu,v=1Au,vHUnUnHAu,v.As such, the 2-D DOA estimatesθk,φk can be obtained via spectrum searching on Pu,v. The detailed steps of the method are shown in Table 1.Table 1
The steps of the proposed method.
Step 1Construct 3-D tensor of the received signalXiStep 2Construct 4-D self-correlation coarray tensorℛ^XiXj from the received signal Xi of the subarrayStep 3Transform the dimension order of the 4-D tensorℛ^XiXj and then merge the dimensions with the same directional factor to get Ui,jStep 4Take the coarray tensorU1,U2,U3,U4 corresponding to the largest continuous virtual array from Ui,jStep 5Transform the dimension of the coarray tensorUi by tensor expansion to obtain structured coarray tensor TStep 6Perform CPD on the structured coarray tensorT to obtain matrix factors cxuk, cyvk containing direction information and calculate θk,φk by spectral peak search


## 3.3. Complexity Analysis

We analyze the computational complexity of the proposed method in the coprime planar array and compare it with the CCT method, SSM method, and PSS method. Let us first analyze the computational complexity of the proposed method. For a coprime planar array, it can be divided into two uniform planar subarrays of sizesM×M and N×N, the proposed method constructs covariance tensors by received signals, and two subarrays can construct four covariance tensors, so the computational complexity of this step is OLM4+N4+2M2N2. Then, by tensoring the covariance tensor, the structure covariance matrix T is obtained, and the computational complexity of this step is O2KP22M+N−P2. Then, we do a tensor decomposition of T, the computational complexity of this step is 16KP22M+N−P2+3K3, and the complexity required for the spectral peak search to calculate the elevation angle and azimuth angle is OngnlP4, where ng=360°/0.01°, nl=90°/0.01°, and ng×nl represents the number of meshes searched by the spectral peak. Therefore, the total complexity of the SCT method is OLM4+N4+2M2N2+2KP22M+N−P2+16KP22M+N−P2+3K3+ngnlP4. The complexity of the PSS method is OngnlN4/M2+M4/N2. The complexity of the SSM method is OM2+N22L+2KP22M+N−P2+P6+ngnlP4. The complexity of the CCT method is OLKM2N2+2KP2M+N−P2+16KP2M+N−P2+3K3+ngnlP4.For a more intuitive representation, Table2 summarizes the complexity of the different methods. As can be seen from Table 2, the computational complexity of the four methods is mainly derived from the spectral peak search, and the computational complexities of the SCT, CCT, and SSM methods are roughly equal and higher than the PSS method.Table 2
Comparison of computational complexity of different methods.
SCTOLM4+N4+2M2N2+2KP22M+N−P2+16KP22M+N−P2+3K3+ngnlP4CCTOLKM2N2+2KP2M+N−P2+16KP2M+N−P2+3K3+ngnlP4SSMOM2+N22L+2KP22M+N−P2+P6+ngnlP4PSSOngnlN4/M2+M4/N2


## 4. Simulation Results

### 4.1. Underdetermined Estimation

First, we simulated the performance of the proposed method in the case of imperfection, whenN  = 2, M  = 3. There are 12 physical sensors; K = 13 sources are from φ=18,28,38,48,22,32,42,52,62,26,36,46,56 for elevation angle and θ=−150,−125,−100,−75,−50,−25,0,25,50,75,100,125,150 for azimuth angle, respectively. The estimation performance of the proposed method is shown in Figure 3, where the number of snapshots is 5000, SNR = 15 dB, and the search step is 0.1°.Figure 3
The proposed method of estimation performance under underdetermined conditions whenM  = 3 and N  = 2.It can be seen that the DOA estimated by the proposed method matches the true ones. Therefore, the proposed method also has good estimation performance under underdetermined conditions.


### 4.2. Comparison of Angle Resolution Capabilities

We compare the proposed method with the method based on CCT and SSM, respectively, and the number of sources is 2. First, we compare the case where the direction of the signals is close to each other, and the azimuth and elevation angles of the two sources are20°,20° and 21°,21°, respectively. To avoid the loss of accuracy of the algorithm estimate due to the small aperture of the array, we set M=7,N=8, the signal-to-noise ratio (SNR) is 10 dB, the number of snapshots is 2000, and the search step is 0.01°. It can be seen from Figure 4 that the proposed method in this paper can distinguish the signal direction very well, while the other two methods cannot distinguish two sources. Compared with SSM, the proposed method has the same virtual array aperture, but the proposed method uses tensors to model the signal and estimate the DOA from multiple dimensions. Compared with the CCT, both the proposed method and CCT are based on tensor modeling of signals, but the proposed method has a larger virtual array aperture in DOA estimation, thereby improving the resolution and accuracy, so the proposed method has higher estimation accuracy and angle resolution.Figure 4
Spectral peak search comparison. (a) SCT. (b) CCT. (c) SSM.
(a)(b)(c)


### 4.3. Comparison of RMSE Performance under the Condition of Different Angle Intervals

We compare the RMSE performance of the proposed method with the other three methods under the condition of different angle intervals. ConsiderM=3,N=2, the search step is 0.01°, the number of snapshots = 2000, and SNR = 10 dB. Firstly, we set the DOAs of the two sources as 20°,20°,23°,23°, and each time, the angle interval change is 1°, so the next DOAs of the sources are 20°,20°,24°,24°. By analogy, the last DOAs are 20°,20°,35°,35°, respectively. The simulation results are shown in Figure 5. When the angle interval is close, the advantage of high tensor dimension makes the RMSE performance of CCT method better than that of SSM, and the proposed method not only has the advantage of high tensor dimension but also has a larger array aperture, making the RMSE performance better than CCT and SSM. The PSS method exchanges a smaller degree of freedom for a larger array aperture, and the RMSE of proposed method is close to that of the PSS method when the angle interval is not less than 5°.Figure 5
The RMSE performance varies with angle interval. (a) Azimuth RMSE varying with the angle interval. (b) Elevation RMSE varying with the angle interval.
(a)(b)


### 4.4. Comparison of RMSE Performance under Two Sources

We compare the estimation accuracy of the proposed method with the subspace-based SSM, PSS [42], and CCT and define the root mean square error as the performance metric, and the RMSE of azimuth can be calculated as(17)RMSE=1RK∑r=1R∑k=1Kφk−φ^r,k2,where R stands for the number of Monte Carlo trails, R = 100, and φ^r,k is the estimated azimuth angle of the kth source for the rth trial, respectively. Consider M=3,N=2, the search step is 0.01°, and the azimuth and elevation angles of the two sources are 20°,20° and 30°,30°, respectively. Figure 6 shows the RMSE performance varying with the number of snapshots, where the simulation conditions of the number of snapshots are 100, 200, 500, 1000, 1500, and 2000 and the SNR = 15 dB. Figure 7 shows the RMSE performance varying with SNR, where the simulation conditions of SNR are -10 dB, -5 dB, 0 dB, 5 dB, 10 dB, and 15 dB, and the number of snapshots = 2000.Figure 6
The RMSE performance varies with the number of snapshots whenK is 2. (a) Azimuth RMSE varying with the number of snapshots. (b) Elevation RMSE varying with the number of snapshots.
(a)(b)Figure 7
The RMSE performance varies with SNR whenK is 2. (a) Azimuth RMSE varying with the SNR. (b) Elevation RMSE varying with the SNR.
(a)(b)As shown in Figures6 and 7, the proposed method has better estimation performance than CCT and SSM. Under the condition of low SNR and low number of snapshots, the proposed method has more advantages and better robustness. PSS is slightly better than the proposed method in precision. The specific reasons are as follows:(1)
Compared with CCT, the proposed method increases the self-correlation of the subarrays, which can increase the aperture of the maximum continuous virtual array, the degrees of freedom, and the estimation accuracy.(2)
Compared with SSM, the proposed method effectively incorporates the multidimensional information, which is not utilized in the matrix-based methods.(3)
Under the condition of the same physical sensor, although PSS has higher estimation accuracy and its estimation accuracy is related to the array aperture of the large subarray, its degree of freedom is limited by the number of array elements of the small subarray, that is,DOF≤minM2,N2−1. The DOA estimation method of PSS is to estimate the signal direction with two sparse uniform subarrays, and the interval between the elements of the subarrays is greater than half a wavelength, so ambiguous peaks will appear to interfere with the DOA estimation, but the corresponding MUSIC spectrums of the two subarrays uniquely generate a common peak at the direction of true DOA, so the estimated accuracy of the PSS depends on the array aperture of the large subarray. When estimating two sources, the main reason why the PSS estimation accuracy is higher than the proposed algorithm when estimating the two sources is that the array aperture of the PSS in the DOA estimate is M−1Nλ2. In the proposed method, according to (11), the virtual array aperture is P−1λ2. In the RMSE performance simulation, P=4,M=3,N=2, and the array aperture of the proposed method is smaller than that of the PSS method, but the theoretical degree of freedom of the PSS algorithm is only 3, which is much lower than the theoretical degrees of freedom of the proposed method.


### 4.5. Comparison of RMSE Performance under Four Sources

We compare the RMSE performance of the proposed method with the other three methods under the condition of estimating four sources. ConsiderM=3,N=2, the search step is 0.01°, and the azimuth and elevation angles of the two sources are 15°,15°, 30°,30°, 45°,45°, and 60°,60°, respectively. Figure 8 shows the RMSE performance varying with the number of snapshots where the simulation conditions of the number of snapshots are 100, 500, 1000, 2000, and 5000 and the SNR = 15 dB. Figure 9 shows the RMSE performance varying with SNR, where the simulation conditions of SNR are -5 dB, 0 dB, 5 dB, 10 dB, and 15 dB and the number of snapshots = 5000.Figure 8
The RMSE performance varies with the number of snapshots whenK is 4. (a) Azimuth RMSE varying with the number of snapshots. (b) Elevation RMSE varying with the number of snapshots.
(a)(b)Figure 9
The RMSE performance varies with SNR whenK is 4. (a) Azimuth RMSE varying with the SNR. (b) Elevation RMSE varying with the SNR.
(a)(b)The PSS method reduces the degrees of freedom in exchange for a higher estimation accuracy, and the maximum number of predictable sources isDOF≤minM2,N2−1. As shown in Figures 8 and 9, we can see that the PSS method fails in the DOA estimation when the number of sources is not less than the number of array elements of the small array. The proposed method has high degrees of freedom after adding the self-covariance information, and the estimation accuracy is better than the other three methods.


## 4.1. Underdetermined Estimation

First, we simulated the performance of the proposed method in the case of imperfection, whenN  = 2, M  = 3. There are 12 physical sensors; K = 13 sources are from φ=18,28,38,48,22,32,42,52,62,26,36,46,56 for elevation angle and θ=−150,−125,−100,−75,−50,−25,0,25,50,75,100,125,150 for azimuth angle, respectively. The estimation performance of the proposed method is shown in Figure 3, where the number of snapshots is 5000, SNR = 15 dB, and the search step is 0.1°.Figure 3
The proposed method of estimation performance under underdetermined conditions whenM  = 3 and N  = 2.It can be seen that the DOA estimated by the proposed method matches the true ones. Therefore, the proposed method also has good estimation performance under underdetermined conditions.


## 4.2. Comparison of Angle Resolution Capabilities

We compare the proposed method with the method based on CCT and SSM, respectively, and the number of sources is 2. First, we compare the case where the direction of the signals is close to each other, and the azimuth and elevation angles of the two sources are20°,20° and 21°,21°, respectively. To avoid the loss of accuracy of the algorithm estimate due to the small aperture of the array, we set M=7,N=8, the signal-to-noise ratio (SNR) is 10 dB, the number of snapshots is 2000, and the search step is 0.01°. It can be seen from Figure 4 that the proposed method in this paper can distinguish the signal direction very well, while the other two methods cannot distinguish two sources. Compared with SSM, the proposed method has the same virtual array aperture, but the proposed method uses tensors to model the signal and estimate the DOA from multiple dimensions. Compared with the CCT, both the proposed method and CCT are based on tensor modeling of signals, but the proposed method has a larger virtual array aperture in DOA estimation, thereby improving the resolution and accuracy, so the proposed method has higher estimation accuracy and angle resolution.Figure 4
Spectral peak search comparison. (a) SCT. (b) CCT. (c) SSM.
(a)(b)(c)


## 4.3. Comparison of RMSE Performance under the Condition of Different Angle Intervals

We compare the RMSE performance of the proposed method with the other three methods under the condition of different angle intervals. ConsiderM=3,N=2, the search step is 0.01°, the number of snapshots = 2000, and SNR = 10 dB. Firstly, we set the DOAs of the two sources as 20°,20°,23°,23°, and each time, the angle interval change is 1°, so the next DOAs of the sources are 20°,20°,24°,24°. By analogy, the last DOAs are 20°,20°,35°,35°, respectively. The simulation results are shown in Figure 5. When the angle interval is close, the advantage of high tensor dimension makes the RMSE performance of CCT method better than that of SSM, and the proposed method not only has the advantage of high tensor dimension but also has a larger array aperture, making the RMSE performance better than CCT and SSM. The PSS method exchanges a smaller degree of freedom for a larger array aperture, and the RMSE of proposed method is close to that of the PSS method when the angle interval is not less than 5°.Figure 5
The RMSE performance varies with angle interval. (a) Azimuth RMSE varying with the angle interval. (b) Elevation RMSE varying with the angle interval.
(a)(b)


## 4.4. Comparison of RMSE Performance under Two Sources

We compare the estimation accuracy of the proposed method with the subspace-based SSM, PSS [42], and CCT and define the root mean square error as the performance metric, and the RMSE of azimuth can be calculated as(17)RMSE=1RK∑r=1R∑k=1Kφk−φ^r,k2,where R stands for the number of Monte Carlo trails, R = 100, and φ^r,k is the estimated azimuth angle of the kth source for the rth trial, respectively. Consider M=3,N=2, the search step is 0.01°, and the azimuth and elevation angles of the two sources are 20°,20° and 30°,30°, respectively. Figure 6 shows the RMSE performance varying with the number of snapshots, where the simulation conditions of the number of snapshots are 100, 200, 500, 1000, 1500, and 2000 and the SNR = 15 dB. Figure 7 shows the RMSE performance varying with SNR, where the simulation conditions of SNR are -10 dB, -5 dB, 0 dB, 5 dB, 10 dB, and 15 dB, and the number of snapshots = 2000.Figure 6
The RMSE performance varies with the number of snapshots whenK is 2. (a) Azimuth RMSE varying with the number of snapshots. (b) Elevation RMSE varying with the number of snapshots.
(a)(b)Figure 7
The RMSE performance varies with SNR whenK is 2. (a) Azimuth RMSE varying with the SNR. (b) Elevation RMSE varying with the SNR.
(a)(b)As shown in Figures6 and 7, the proposed method has better estimation performance than CCT and SSM. Under the condition of low SNR and low number of snapshots, the proposed method has more advantages and better robustness. PSS is slightly better than the proposed method in precision. The specific reasons are as follows:(1)
Compared with CCT, the proposed method increases the self-correlation of the subarrays, which can increase the aperture of the maximum continuous virtual array, the degrees of freedom, and the estimation accuracy.(2)
Compared with SSM, the proposed method effectively incorporates the multidimensional information, which is not utilized in the matrix-based methods.(3)
Under the condition of the same physical sensor, although PSS has higher estimation accuracy and its estimation accuracy is related to the array aperture of the large subarray, its degree of freedom is limited by the number of array elements of the small subarray, that is,DOF≤minM2,N2−1. The DOA estimation method of PSS is to estimate the signal direction with two sparse uniform subarrays, and the interval between the elements of the subarrays is greater than half a wavelength, so ambiguous peaks will appear to interfere with the DOA estimation, but the corresponding MUSIC spectrums of the two subarrays uniquely generate a common peak at the direction of true DOA, so the estimated accuracy of the PSS depends on the array aperture of the large subarray. When estimating two sources, the main reason why the PSS estimation accuracy is higher than the proposed algorithm when estimating the two sources is that the array aperture of the PSS in the DOA estimate is M−1Nλ2. In the proposed method, according to (11), the virtual array aperture is P−1λ2. In the RMSE performance simulation, P=4,M=3,N=2, and the array aperture of the proposed method is smaller than that of the PSS method, but the theoretical degree of freedom of the PSS algorithm is only 3, which is much lower than the theoretical degrees of freedom of the proposed method.


## 4.5. Comparison of RMSE Performance under Four Sources

We compare the RMSE performance of the proposed method with the other three methods under the condition of estimating four sources. ConsiderM=3,N=2, the search step is 0.01°, and the azimuth and elevation angles of the two sources are 15°,15°, 30°,30°, 45°,45°, and 60°,60°, respectively. Figure 8 shows the RMSE performance varying with the number of snapshots where the simulation conditions of the number of snapshots are 100, 500, 1000, 2000, and 5000 and the SNR = 15 dB. Figure 9 shows the RMSE performance varying with SNR, where the simulation conditions of SNR are -5 dB, 0 dB, 5 dB, 10 dB, and 15 dB and the number of snapshots = 5000.Figure 8
The RMSE performance varies with the number of snapshots whenK is 4. (a) Azimuth RMSE varying with the number of snapshots. (b) Elevation RMSE varying with the number of snapshots.
(a)(b)Figure 9
The RMSE performance varies with SNR whenK is 4. (a) Azimuth RMSE varying with the SNR. (b) Elevation RMSE varying with the SNR.
(a)(b)The PSS method reduces the degrees of freedom in exchange for a higher estimation accuracy, and the maximum number of predictable sources isDOF≤minM2,N2−1. As shown in Figures 8 and 9, we can see that the PSS method fails in the DOA estimation when the number of sources is not less than the number of array elements of the small array. The proposed method has high degrees of freedom after adding the self-covariance information, and the estimation accuracy is better than the other three methods.


## 5. Conclusion

In this paper, we propose a two-dimensional DOA estimation method based on self-correlation tensor for the problem of statistics of signals underutilized in tensor DOA estimation and solve the problem that tensors are insufficient in virtual array aperture construction in DOA estimation. First, we use tensor to represent the received signal, then construct the self-correlation coarray tensor of the received signal, then merge the dimensions with the same direction factor and extract the coarray tensor corresponding to the largest continuous virtual array, obtain the structural coarray tensor with multirank by tensor expansion, and finally approach the CPD to the structural coarray tensor to obtain DOA estimation with super-resolution. Through simulation and proof, the proposed method effectively expands the maximum virtual array aperture. The proposed method can effectively enlarge the aperture of the maximum continuous virtual planar array and increase the degree of freedom. Compared with CCT and SSM, simulation results show that under the same conditions, the proposed method has better performance.



---

*Source: 7999641-2022-08-08.xml*