# Source Localization Using Distributed Electromagnetic Vector Sensors

Authors: Tingping Zhang; Di Wan; Xinhai Wang; Fangqing Wen
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Electromagnetic vector sensor (EVS) array has drawn extensive attention in the past decades, since it offers two-dimensional direction-of-arrival (2D-DOA) estimation and additional polarization information of the incoming source. Most of the existing works concerning EVS array are focused on parameter estimation with special array architecture, e.g., uniform manifold and sparse arrays. In this paper, we consider a more general scenario that EVS array is distributed in an arbitrary geometry, and a novel estimator is proposed. Firstly, the covariance tensor model is established, which can make full use of the multidimensional structure of the array measurement. Then, the higher-order singular value decomposition (HOSVD) is adopted to obtain a more accurate signal subspace. Thereafter, a novel rotation invariant relation is exploited to construct a normalized Poynting vector, and the vector cross-product technique is utilized to estimate the 2D-DOA. Based on the previous obtained 2D-DOA, the polarization parameter can be easily achieved via the least squares method. The proposed method is suitable for EVS array with arbitrary geometry, and it is insensitive to the spatially colored noise. Therefore, it is more flexible than the state-of-the-art algorithms. Finally, numerical simulations are carried out to verify the effectiveness of the proposed estimator.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation using sensor array is one of the most popular methods for source localization, and it has struck a series of technical elevations in wireless communications [1–3], radars, sonars, etc. As a nonlinear issue, many efforts have been devoted and many superresolution methodologies have been developed, e.g., the estimation method of signal parameters via rotational invariance technique (ESPRIT) [4, 5], maximum-likelihood [6], multiple signal classification (MUSIC) [7–9], sparsity-aware strategy [10, 11], and tensor approach [12–14]. Nevertheless, most of the current works are concentrating on one-dimensional (1D) DOA estimation. A small number of works have paid attention to the two-dimensional (2D) DOA estimation problem, by cooperating with nonlinear scaler arrays [15, 16], e.g., L-shaped array, circular array, and rectangular array. A common characteristic of these arrays is that sensors must be strictly nonlinear. Besides, the sensors’ locations must be known priors, otherwise algorithms will fail to work. As an alternative of nonlinear scaler array, electromagnetic vector sensor (EVS) array has drawn more and more interest in the past decades. Unlike scaler arrays, a signal EVS consists of six collocated components, i.e., three orthogonally oriented dipoles and three orthogonally oriented loops, which measure the electric field and magnetic field, respectively. EVS array can provide not only 2D-DOA estimation but also the polarization status of the incoming sources, which maybe particularly important in detecting weak signal [17].Many attempts have been done with respect to 2D-DOA estimation using EVS array, and various estimators have been proposed. The most popular method is the vector cross product [18], which is based on the Poynting vector of the polarization steering vector. The vector cross-product technique is suitable for a single EVS. Another important branch is combining the traditional subspace method with vector cross-product. In [19], MUSIC-like algorithm was presented. However, it is computationally inefficient owing to high-dimensional spectrum grid search. To avoid such drawback, the ESPRIT-like approach was investigated in [20]. It firstly obtains the signal subspace via eigendecomposition, and then, a normalized Poynting vector is constructed; finally, the vector cross product is performed to achieve 2D-DOA estimation. ESPRIT-like approach is suitable for arbitrary array geometry, and it offers closed-form solution for parameter estimation. Similar works concerning subspace methods and vector cross-product technique have been done in [21–24]. To further lower the computational burden, the propagator method was carried out in [25], which replaces the eigendecomposition with least squares. To exploit the tensor nature, the tensor-aware estimators were derived in [26, 27]. Since the tensor algebra can suppress the noise much better than the matrix-based tools, tensor methods provide more accurate estimation performance than the matrix algorithms. Some efforts have been devoted to further improve the estimation accuracy or reduce the mutual coupling effort. In [28], the coprime array (see [29] and the references therein) geometry was introduced into EVS, which occupies much larger aperture than the traditional uniform manifold. In [30, 31], the EVS architecture has been redesigned. The component number in the improved EVS is much smaller than six as to release the mutual coupling. In addition, the EVS array has been extended to MIMO radar in [32–36], which bring new insights into parameter estimation to MIMO systems. It should be emphasized that most of the existing algorithms are only effective to specific array geometry, e.g., uniform linear array and coprime array. In the presence of limited space, these array geometries may be not practical. Besides, none of the existing works has taken the spatially colored noise into consideration.In this paper, we generalize the issue of source localization using the EVS array. We consider that the EVSs are displaced in arbitrary geometry. Moreover, spatially colored noise appears in the array. An improved higher-order singular value decomposition (HOSVD) estimator is proposed. Firstly, the temporal cross-correlation scheme is utilized to eliminate the spatially colored noise. Thereafter, the crosscovariance measurement is arranged into a fourth-order tensor. The HOSVD is adopted to obtain an enhanced signal subspace. Then, a normalized Poynting vector is constructed, and 2D-DOAs are obtained via the shift invariance property of the previously constructed Poynting vector. Finally, the polarization parameters are achieved via the least squares technique. Our algorithm is insensitive to the sensors’ position and spatially colored noise. Numerical examples are given to show the superiority of the proposed algorithm.


### 1.1. Notations

Throughout the paper, lowercase letters represent vectors and uppercase letters denote matrices, respectively;IN is a N×N identity matrix; ⊗, ⊙, and ∗ denote the Kronecker product, the Khatri-Rao product, and the vector cross product, respectively; the superscript XT, XH, X−1, and X† indicate the operations of transpose, Hermitian transpose, inverse, and pseudoinverse, respectively; DnA denotes a diagonal matrix with the diagonal entities which are the nth row of A; am accounts for the mth entity of the vector a; ⋅F accounts for the Frobenius norm; phase• is to get the phase; and E⋅ returns the mathematical expectation.


## 1.1. Notations

Throughout the paper, lowercase letters represent vectors and uppercase letters denote matrices, respectively;IN is a N×N identity matrix; ⊗, ⊙, and ∗ denote the Kronecker product, the Khatri-Rao product, and the vector cross product, respectively; the superscript XT, XH, X−1, and X† indicate the operations of transpose, Hermitian transpose, inverse, and pseudoinverse, respectively; DnA denotes a diagonal matrix with the diagonal entities which are the nth row of A; am accounts for the mth entity of the vector a; ⋅F accounts for the Frobenius norm; phase• is to get the phase; and E⋅ returns the mathematical expectation.


## 2. Preliminaries and Problem Formulation

### 2.1. Tensor Preliminaries

Let us first introduce some preliminaries concerning tensors. Interested readers are recommend to refer to [37] for more details.Definition 1 (tensor unfolding).
LetX∈ℂI1×⋯IN be an N-order tensor. The n-mode n=1,2,⋯,N tensor unfolding of X is denoted by Xn, which is a In×I1⋯IN matrix. The i1,⋯,inth entity of X is the in,jth element of Xn, j=1+∑k=1,k≠nNik−1Jk, andJk=∏m=1,m≠nk−1Im.Definition 2 (mode-n tensor-matrix product).
The mode-n product of X∈CI1×I2×⋯IN by a matrix A∈CJn×In, is denoted by Y=X×nA, where Y∈CI1×I2×⋯In−1×Jn×In+1⋯×IN−1×IN and Yi1×⋯in−1×in×in+1⋯iN=∑in=1InXi1×⋯in−1×in×in+1⋯iN⋅Ajn,in.Definition 3 (important principles of the mode-n tensor-matrix product).
The following principles are important for mode-n tensor-matrix product:
(1)X×n⋅A×m⋅B=X×m⋅B×n⋅A,m≠n,X×n⋅A×m⋅B=X×n⋅B⋅A,m=n,X×1⋅A1×2⋯ANn=An⋅Xn⋅AN⊗⋯An+1⊗An−1⋯⊗A1T.Definition 4 (HOSVD).
The HOSVD of anN-order tensor X∈ℂI1×⋯IN is denoted by
(2)X=G×1U×21U×32⋯×NUN,where G∈ℂJ1×⋯JN is the so-called core tensor and Un∈ℂIn×⋯JN is the unitary matrix. In mode-n tensor unfolding format, (2) is rewritten as
(3)Xn=UnGnUn+1⊗⋯UN⊗U1⋯⊗Un−1T.


### 2.2. Preliminaries of EVS

A single EVS consists of six-component electromagnetic antennas. For an incoming source that impinges on an isolated EVS, the electromagnetic response of the EVS is given by [20]
(4)b≜eh≜cosϕcosθ−sinϕsinϕcosθ−sinθ−sinϕcosϕ0cosϕ0−cosϕcosθ−sinϕcosθsinθ⏟Dsinγejηcosγ⏟g,where θ, ϕ, γ, and η are the elevation angle, azimuth angle, auxiliary polarization angle, and polarization phase difference of the source, respectively; e∈ℂ3×1 and h∈ℂ3×1 are measurements that sense the electric field and magnetic field, respectively. D∈ℂ6×2 is the direction parameter-depend only matrix, and g∈ℂ2×1 is the polarization parameter-depend only vector.An important characteristic is thate and h are mutually orthogonal to each other. Moreover, the normalized Poynting-vector fulfills [20]
(5)uvw≜ee∗h∗h≜sinθcosϕsinθsinϕcosθ.


### 2.3. Problem Formulation

In this paper, we consider anM-element EVS array with arbitrary geometry, as illustrated in Figure 1. Suppose that the mm=1,2,⋯,Mth EVS is located at coordinate xm,ym,zm. Assume that there are K far-field source signals that impinge on the EVS array. Let θk,ϕk,γk,ηk be an angle pair that parameterized the kth (k=1,2,⋯,K) source signal. The received array signal from the EVS array can be written as [20]
(6)yt=∑k=1Kak⊗bkskt+nt=A⊙Bst+nt,where ak≜1,e−jπτ2,k/λ,⋯,e−jπτM,k/λT is the spatial steering vector of the kth source signal, λ is the carrier wavelength, and τm,k=xmsinθkcosϕk+ymsinθksinϕk+zmcosθk. bk is the associated polarization steering vector, skt is the kth source signal, and st=s1t,s2t,⋯,sKtT is the source signal vector. nt is the array noise, which is assumed to be spatially correlated, temporally independent Gaussian process, i.e.,
(7)Ent1nHt2=0,t1≠t2,Q,t1=t2,where Q is the covariance matrix of the noise, which is a Hermitian matrix. In the presence of Gaussian white noise, Q=σ2I.Figure 1
Illustration of the EVS array with arbitrary geometry.


## 2.1. Tensor Preliminaries

Let us first introduce some preliminaries concerning tensors. Interested readers are recommend to refer to [37] for more details.Definition 1 (tensor unfolding).
LetX∈ℂI1×⋯IN be an N-order tensor. The n-mode n=1,2,⋯,N tensor unfolding of X is denoted by Xn, which is a In×I1⋯IN matrix. The i1,⋯,inth entity of X is the in,jth element of Xn, j=1+∑k=1,k≠nNik−1Jk, andJk=∏m=1,m≠nk−1Im.Definition 2 (mode-n tensor-matrix product).
The mode-n product of X∈CI1×I2×⋯IN by a matrix A∈CJn×In, is denoted by Y=X×nA, where Y∈CI1×I2×⋯In−1×Jn×In+1⋯×IN−1×IN and Yi1×⋯in−1×in×in+1⋯iN=∑in=1InXi1×⋯in−1×in×in+1⋯iN⋅Ajn,in.Definition 3 (important principles of the mode-n tensor-matrix product).
The following principles are important for mode-n tensor-matrix product:
(1)X×n⋅A×m⋅B=X×m⋅B×n⋅A,m≠n,X×n⋅A×m⋅B=X×n⋅B⋅A,m=n,X×1⋅A1×2⋯ANn=An⋅Xn⋅AN⊗⋯An+1⊗An−1⋯⊗A1T.Definition 4 (HOSVD).
The HOSVD of anN-order tensor X∈ℂI1×⋯IN is denoted by
(2)X=G×1U×21U×32⋯×NUN,where G∈ℂJ1×⋯JN is the so-called core tensor and Un∈ℂIn×⋯JN is the unitary matrix. In mode-n tensor unfolding format, (2) is rewritten as
(3)Xn=UnGnUn+1⊗⋯UN⊗U1⋯⊗Un−1T.


## 2.2. Preliminaries of EVS

A single EVS consists of six-component electromagnetic antennas. For an incoming source that impinges on an isolated EVS, the electromagnetic response of the EVS is given by [20]
(4)b≜eh≜cosϕcosθ−sinϕsinϕcosθ−sinθ−sinϕcosϕ0cosϕ0−cosϕcosθ−sinϕcosθsinθ⏟Dsinγejηcosγ⏟g,where θ, ϕ, γ, and η are the elevation angle, azimuth angle, auxiliary polarization angle, and polarization phase difference of the source, respectively; e∈ℂ3×1 and h∈ℂ3×1 are measurements that sense the electric field and magnetic field, respectively. D∈ℂ6×2 is the direction parameter-depend only matrix, and g∈ℂ2×1 is the polarization parameter-depend only vector.An important characteristic is thate and h are mutually orthogonal to each other. Moreover, the normalized Poynting-vector fulfills [20]
(5)uvw≜ee∗h∗h≜sinθcosϕsinθsinϕcosθ.


## 2.3. Problem Formulation

In this paper, we consider anM-element EVS array with arbitrary geometry, as illustrated in Figure 1. Suppose that the mm=1,2,⋯,Mth EVS is located at coordinate xm,ym,zm. Assume that there are K far-field source signals that impinge on the EVS array. Let θk,ϕk,γk,ηk be an angle pair that parameterized the kth (k=1,2,⋯,K) source signal. The received array signal from the EVS array can be written as [20]
(6)yt=∑k=1Kak⊗bkskt+nt=A⊙Bst+nt,where ak≜1,e−jπτ2,k/λ,⋯,e−jπτM,k/λT is the spatial steering vector of the kth source signal, λ is the carrier wavelength, and τm,k=xmsinθkcosϕk+ymsinθksinϕk+zmcosθk. bk is the associated polarization steering vector, skt is the kth source signal, and st=s1t,s2t,⋯,sKtT is the source signal vector. nt is the array noise, which is assumed to be spatially correlated, temporally independent Gaussian process, i.e.,
(7)Ent1nHt2=0,t1≠t2,Q,t1=t2,where Q is the covariance matrix of the noise, which is a Hermitian matrix. In the presence of Gaussian white noise, Q=σ2I.Figure 1
Illustration of the EVS array with arbitrary geometry.


## 3. The Proposed Algorithm

### 3.1. Spatially Colored Noise Suppression

Suppose the source signals are uncorrelated andst is uncorrelated with nt; then, the covariance of yt is given by
(8)Ry=A⊙BRsA⊙BH+Q,where Rs is the covariance matrix of st. In the presence of spatially colored noise, the signal subspace will be corrupted by the noise subspace; thus, the performance of the traditional subspace methods would be decreased or even fail to work. Therefore, it is necessary for us to eliminate the effect of the spatially colored noise. As illustrated in (7), the noise between different snapshots is uncorrelated; thus, we can construct the following cross-correlation matrix:
(9)R≜Eyt1yHt2≜A⊙BΛA⊙BH,where Λ=Est1sHt2 is the cross-correlation matrix of the source signal. In the presence of temporally correlated source signal, e.g., sinusoidal signal, Λ is a diagonal matrix. When L snapshots are available, R can be estimated via
(10)R^≜1L‐1∑l=1L−1ylyHl+1.Obviously, the noise is eliminated inR. In practice, R can be approximated by its truncated eigenvalue decomposition as
(11)R≈UdΣdUdH,where Σd is the dominate eigenvalue matrix, Us consist of the K dominate eigenvectors.


### 3.2. HOSVD

Actually,R can be arranged as a fourth-order tensor as [35]
(12)Rm,n,p,q=∑k=1KΛk,kAm,kBn,kA∗p,kB∗q,km,p∈1,⋯,M;n,q∈1,⋯,6.The HOSVD ofR is given by
(13)R=G×1⋅U1×1⋅U2×2⋅U3×3⋅U4,where G∈ℂM×6×M×6 is the core tensor, U1∈ℂM×M, U2∈ℂ6×6, U3∈ℂM×M, and U4∈ℂ6×6 are the unitary matrices, which consist of the left singular vectors of the mode-1, mode-2, mode-3, and mode-4 unfolding of R, respectively. Since the rank of R is K, we can approximate R by the truncated HOSVD as [35]
(14)Ra=Gs×1⋅U1s×1⋅U2s×2⋅U3s×3⋅U4s,where U1s, U2s, U3s, and U4s consist of the eigenvectors associated to the K dominate eigenvalues. Gs is called the signal component of G, which is given by
(15)Gs=R×1U1s×2HU2s×3HU3s×4HU4sH.InsertingGs into (14) yields
(16)Ra=R×1⋅U1sU1sH×2⋅U2sU2sH×3⋅U3sU3sH×4⋅U4sU4sH.By Hermitian unfoldingRa, we can construct an enhanced covariance matrix as
(17)Ra=U1sU1sH⊗U2sU2sHRU3sU3sH⊗U4sU4sH∗.Combined with the result of (11), we can obtain
(18)Ra=U1sU1sH⊗U2sU2sHUdΣdUdHU3sU3sH⊗U4sU4sH∗.SinceR is Hermitian, we have U1s=U3s∗ and U2s=U4s∗. Performing eigendecomposition on Ra, we can obtain a signal subspace as
(19)Es=U1sU1sH⊗U2sU2sHUd.Since the virtual response matrixA⊙B spans the same subspace with Es, we have
(20)Es=A⊙BT,where T∈ℂK×K is a nonsingular matrix.


### 3.3. 2D-DOA Estimation

Letbk=bk1,bk2,⋯,bk6T. We can get the following relationship:
(21)AD1B=AD2BΦ1,2,AD1B=AD3BΦ1,3,⋮AD1B=AD6BΦ1,6,with
(22)Φ1,q=β11,qβ21,q⋱βK1,q,q=2,3,⋯,6,where βk1,q=bkq/bk1. Define the following selection matrices:
(23)Jq≜IM⊗i6,q,where i6,q is the qth row of I6. Therefore, the relationships in (21) becomes
(24)J1A⊙B=J2A⊙BΦ1,2,J1A⊙B=J3A⊙BΦ1,3,⋮J1A⊙B=J6A⊙BΦ1,6.Inserting (24) into (20) establishes
(25)J1Es=J2EsT−1Φ1,2T,J1Es=J3EsT−1Φ1,3T,⋮J1Es=J6EsT−1Φ1,6T,Equivalently,(26)J2Es†J1Es=T−1Φ1,2T,TJ3Es†J1EsT−1=Φ1,3,⋮TJ6Es†J1EsT−1=Φ1,6.According to the first line of (26), we can perform eigendecomposition of J2Es†J1Es, the eigenvalues of which reveal the estimation of Φ1,2, and the associated eigenvectors form the estimation of T. Compute the left parts of the reminder equations, and then, we can get the estimation of Φ1,3, Φ1,4, Φ1,5, and Φ1,6, respectively.In fact,bk can be written in the following normalized format:
(27)bk=bk11βk1,2⋮βk1,6.Then, we have(28)bk11βk1,2βk1,3∗bk1βk1,4βk1,5βk1,6∗=bk121βk1,2βk1,3∗βk1,4βk1,5βk1,6∗.Sincebk12 is a constant, it can be easily removed from the normalized operation. Let ek=bk1,bk2,bk3T and hk=bk4,bk5,bk6T. According to the property in (5), we have
(29)ekek∗hk∗hk=sinθkcosϕksinθksinϕkcosθk.Accordingly, we can get(30)ukvkwk≜1βk1,2βk1,31βk1,2βk1,3∗βk1,4βk1,5βk1,6∗βk1,4βk1,5βk1,6=sinθkcosϕksinθksinϕkcosθk.Denote the estimates ofβk1,q as β^k1,q; then, we can get the estimation of uk, vk, and wk, which are denoted by u^k, v^k, and w^k, respectively. Then, 2D-DOA estimation can be achieved via
(31)θ^k=arccosw^k,ϕ^k=arctanv^ku^k.Obviously, the estimated 2D-DOA are automatically paired.


### 3.4. Polarization Parameter Estimation

It is easy to find(32)B~m≜BDmA=Be−jπτm,1/λe−jπτm,2/λ⋱e−jπτm,K/λ.So thekth column of B~m is
(33)b~k,m=e−jπτm,2/λbk.According to (4), we should have
(34)b~k,m=e−jπτm,2/λDkgk,where Dk and gk are the direction-only and polarization-only variables. Define
(35)Cm≜iM,mT⊗IM.Then,B~m can also be expressed as
(36)B~m=CmA⊙B.Therefore,B~m can be estimated via
(37)B^m=CmEsT−1,the kth column of which is denoted by b^k,m. After the 2D-DOA have been estimated, we can construct the matrix D^k according to (4); then, we can compute the following least squares:
(38)g^k=D^k†b^k,m.Obviously,g^k is the estimation of e−jπτm,2/λgk. Finally, the polarization parameters can be achieved via
(39)γ^k=arctang^k1g^k2,η^k=phaseg^k1g^k2.Also, the polarization parameters are paired automatically.


## 3.1. Spatially Colored Noise Suppression

Suppose the source signals are uncorrelated andst is uncorrelated with nt; then, the covariance of yt is given by
(8)Ry=A⊙BRsA⊙BH+Q,where Rs is the covariance matrix of st. In the presence of spatially colored noise, the signal subspace will be corrupted by the noise subspace; thus, the performance of the traditional subspace methods would be decreased or even fail to work. Therefore, it is necessary for us to eliminate the effect of the spatially colored noise. As illustrated in (7), the noise between different snapshots is uncorrelated; thus, we can construct the following cross-correlation matrix:
(9)R≜Eyt1yHt2≜A⊙BΛA⊙BH,where Λ=Est1sHt2 is the cross-correlation matrix of the source signal. In the presence of temporally correlated source signal, e.g., sinusoidal signal, Λ is a diagonal matrix. When L snapshots are available, R can be estimated via
(10)R^≜1L‐1∑l=1L−1ylyHl+1.Obviously, the noise is eliminated inR. In practice, R can be approximated by its truncated eigenvalue decomposition as
(11)R≈UdΣdUdH,where Σd is the dominate eigenvalue matrix, Us consist of the K dominate eigenvectors.


## 3.2. HOSVD

Actually,R can be arranged as a fourth-order tensor as [35]
(12)Rm,n,p,q=∑k=1KΛk,kAm,kBn,kA∗p,kB∗q,km,p∈1,⋯,M;n,q∈1,⋯,6.The HOSVD ofR is given by
(13)R=G×1⋅U1×1⋅U2×2⋅U3×3⋅U4,where G∈ℂM×6×M×6 is the core tensor, U1∈ℂM×M, U2∈ℂ6×6, U3∈ℂM×M, and U4∈ℂ6×6 are the unitary matrices, which consist of the left singular vectors of the mode-1, mode-2, mode-3, and mode-4 unfolding of R, respectively. Since the rank of R is K, we can approximate R by the truncated HOSVD as [35]
(14)Ra=Gs×1⋅U1s×1⋅U2s×2⋅U3s×3⋅U4s,where U1s, U2s, U3s, and U4s consist of the eigenvectors associated to the K dominate eigenvalues. Gs is called the signal component of G, which is given by
(15)Gs=R×1U1s×2HU2s×3HU3s×4HU4sH.InsertingGs into (14) yields
(16)Ra=R×1⋅U1sU1sH×2⋅U2sU2sH×3⋅U3sU3sH×4⋅U4sU4sH.By Hermitian unfoldingRa, we can construct an enhanced covariance matrix as
(17)Ra=U1sU1sH⊗U2sU2sHRU3sU3sH⊗U4sU4sH∗.Combined with the result of (11), we can obtain
(18)Ra=U1sU1sH⊗U2sU2sHUdΣdUdHU3sU3sH⊗U4sU4sH∗.SinceR is Hermitian, we have U1s=U3s∗ and U2s=U4s∗. Performing eigendecomposition on Ra, we can obtain a signal subspace as
(19)Es=U1sU1sH⊗U2sU2sHUd.Since the virtual response matrixA⊙B spans the same subspace with Es, we have
(20)Es=A⊙BT,where T∈ℂK×K is a nonsingular matrix.


## 3.3. 2D-DOA Estimation

Letbk=bk1,bk2,⋯,bk6T. We can get the following relationship:
(21)AD1B=AD2BΦ1,2,AD1B=AD3BΦ1,3,⋮AD1B=AD6BΦ1,6,with
(22)Φ1,q=β11,qβ21,q⋱βK1,q,q=2,3,⋯,6,where βk1,q=bkq/bk1. Define the following selection matrices:
(23)Jq≜IM⊗i6,q,where i6,q is the qth row of I6. Therefore, the relationships in (21) becomes
(24)J1A⊙B=J2A⊙BΦ1,2,J1A⊙B=J3A⊙BΦ1,3,⋮J1A⊙B=J6A⊙BΦ1,6.Inserting (24) into (20) establishes
(25)J1Es=J2EsT−1Φ1,2T,J1Es=J3EsT−1Φ1,3T,⋮J1Es=J6EsT−1Φ1,6T,Equivalently,(26)J2Es†J1Es=T−1Φ1,2T,TJ3Es†J1EsT−1=Φ1,3,⋮TJ6Es†J1EsT−1=Φ1,6.According to the first line of (26), we can perform eigendecomposition of J2Es†J1Es, the eigenvalues of which reveal the estimation of Φ1,2, and the associated eigenvectors form the estimation of T. Compute the left parts of the reminder equations, and then, we can get the estimation of Φ1,3, Φ1,4, Φ1,5, and Φ1,6, respectively.In fact,bk can be written in the following normalized format:
(27)bk=bk11βk1,2⋮βk1,6.Then, we have(28)bk11βk1,2βk1,3∗bk1βk1,4βk1,5βk1,6∗=bk121βk1,2βk1,3∗βk1,4βk1,5βk1,6∗.Sincebk12 is a constant, it can be easily removed from the normalized operation. Let ek=bk1,bk2,bk3T and hk=bk4,bk5,bk6T. According to the property in (5), we have
(29)ekek∗hk∗hk=sinθkcosϕksinθksinϕkcosθk.Accordingly, we can get(30)ukvkwk≜1βk1,2βk1,31βk1,2βk1,3∗βk1,4βk1,5βk1,6∗βk1,4βk1,5βk1,6=sinθkcosϕksinθksinϕkcosθk.Denote the estimates ofβk1,q as β^k1,q; then, we can get the estimation of uk, vk, and wk, which are denoted by u^k, v^k, and w^k, respectively. Then, 2D-DOA estimation can be achieved via
(31)θ^k=arccosw^k,ϕ^k=arctanv^ku^k.Obviously, the estimated 2D-DOA are automatically paired.


## 3.4. Polarization Parameter Estimation

It is easy to find(32)B~m≜BDmA=Be−jπτm,1/λe−jπτm,2/λ⋱e−jπτm,K/λ.So thekth column of B~m is
(33)b~k,m=e−jπτm,2/λbk.According to (4), we should have
(34)b~k,m=e−jπτm,2/λDkgk,where Dk and gk are the direction-only and polarization-only variables. Define
(35)Cm≜iM,mT⊗IM.Then,B~m can also be expressed as
(36)B~m=CmA⊙B.Therefore,B~m can be estimated via
(37)B^m=CmEsT−1,the kth column of which is denoted by b^k,m. After the 2D-DOA have been estimated, we can construct the matrix D^k according to (4); then, we can compute the following least squares:
(38)g^k=D^k†b^k,m.Obviously,g^k is the estimation of e−jπτm,2/λgk. Finally, the polarization parameters can be achieved via
(39)γ^k=arctang^k1g^k2,η^k=phaseg^k1g^k2.Also, the polarization parameters are paired automatically.


## 4. Algorithm Analysis

Now, we have accomplished the estimation algorithm for distributed EVS array. To help the reader better understand the proposed algorithm, we summarize its algorithmic steps in Table1.Table 1
Algorithmic steps of the proposed algorithm.
Step No.OperationStep 1Estimate the noiseless covariance matrixR through (10)Step 2Rearrange the covariance measurement into a fourth-order tensorR via (12), and perform HOSVD on RStep 3ConstructRa from (17), and then, get Ra via the Hermitian unfolding of RaStep 4Perform eigendecomposition onRa to get the signal subspace EsStep 5CalculateJ2Es†J1Es, and perform eigendecomposition on it to get Φ1,2 and TStep 6To getΦ1,3, Φ1,4, Φ1,5, and Φ1,6 via (26)Step 7EstimateukvkwkT via (30), and then, get the 2D-DOA through (31)Step 8EstimateB~m via (37), and compute g^k through (38). Finally, get the polarization parameters via (39)


### 4.1. Complexity

Next, we summarize the main complexity (number of complex multiplication that is required) of typical algorithms. The main complexity of ESPRIT in [20] is eigendecomposition, which is on the order O63M3. The main computational burden of the PARAFAC estimator in [28] is the iteration, which requires OMK2+6K2+LK2 complex multiplications. In the proposed algorithm, the main complexity is the HOSVD, which is the same as that of the eigendecomposition. We listed the main complexity of the above algorithms in Table 2. From this perspective, the proposed approach is less efficient than PARAFAC. However, it should be noticed, which will be shown in the simulation section, that the proposed approach is suitable for colored noise, while the others are sensitive to the colored noise.Table 2
Comparison of the complexity.
MethodComplexityESPRITO63M3PARAFACOMK2+6K2+LK2ProposedO63M3


### 4.2. Identifiability

The identifiability of the proposed approach is equal to the maximumK, which is also the same as that of ESPRIT. According to (26), the proposed approach is effective if K≤6. Therefore, the proposed algorithm can identify K=6 sources at most. In comparison, PARAFAC can identify M+6+L−3/2 sources. Usually, L is much larger than 6M; thus, PARAFAC provides much better identifiability than ESPRIT and the proposed approach.


## 4.1. Complexity

Next, we summarize the main complexity (number of complex multiplication that is required) of typical algorithms. The main complexity of ESPRIT in [20] is eigendecomposition, which is on the order O63M3. The main computational burden of the PARAFAC estimator in [28] is the iteration, which requires OMK2+6K2+LK2 complex multiplications. In the proposed algorithm, the main complexity is the HOSVD, which is the same as that of the eigendecomposition. We listed the main complexity of the above algorithms in Table 2. From this perspective, the proposed approach is less efficient than PARAFAC. However, it should be noticed, which will be shown in the simulation section, that the proposed approach is suitable for colored noise, while the others are sensitive to the colored noise.Table 2
Comparison of the complexity.
MethodComplexityESPRITO63M3PARAFACOMK2+6K2+LK2ProposedO63M3


## 4.2. Identifiability

The identifiability of the proposed approach is equal to the maximumK, which is also the same as that of ESPRIT. According to (26), the proposed approach is effective if K≤6. Therefore, the proposed algorithm can identify K=6 sources at most. In comparison, PARAFAC can identify M+6+L−3/2 sources. Usually, L is much larger than 6M; thus, PARAFAC provides much better identifiability than ESPRIT and the proposed approach.


## 5. Simulation Results

Numerical simulations are carried out to show the superiority of the proposed algorithm. In the following simulations, we assumed an EVS array withM=8 elements, the coordinates of the EVSs xm,ym,zm fulfill uniform distribution in −2λ,2λ. Assume that there are K=3 uncorrelated far-field sinusoidal signals with frequency which is 0.2 MHz, 0.4 MHz, and 0.85 MHz, respectively. The associated parameters are θ=10°,30°,50°T, ϕ=15°,‐25°,35°, γ=10°,30°,60°, and η=20°,−40°,55°, respectively. L samples are collected. The signal-to-noise ratio (SNR) in the simulation is defined as the ratio of signal power to noise power in (6). The spatially correlated noise is modeled as an autoregressive (AR) process with the coefficients 1,β,αT. All the results rely on 200 independent trials.Example 1.
Scatter results of the proposed algorithm. In such simulation, we considerβ=−1, α=0.8, and L =1000 and SNR is set to 10 dB. The scatter results of 2D-DOA estimation and polarization parameter estimation are depicted in Figure 2. It is obvious to us that all the parameters can be accurately estimated and correctly paired. It is evident that the proposed algorithm is suitable for distributed EVSs.Example 2.
The root mean square error (RMSE) performance versus SNR. Herein, RMSE for a parameterω=w1,w2,⋯,ωKT is defined as
(40)RMSE=1K∑k=1K1200∑t=1200w∧t,k−wk2,where w^t,k is the estimation of wk in the tth trial. In this simulation, we consider β=−1, α=0.8, and L =1000. For comparison purposes, we add the performance of the traditional ESPRIT algorithm [20] and the PARAFAC algorithm [28]. Figure 3 presents the average RMSE on 2D-DOA estimation (marked with suffix “-d”) and average RMSE on polarization parameter estimation (marked with suffix “-p”) in the presence of white Gaussian noise. It is seen that the proposed algorithm provides better RMSE than ESPRIT, but a slight worse RMSE performance than PARAFAC. Figure 4 displays the RMSE results with colored noise. Obviously, all the algorithms achieve better RMSE with higher SNR. As expected, both ESPRIT and PARAFAC perform worse at low SNR regions since they are affected by the spatially colored noise. On the other hand, the proposed method offers much better performance than ESPRIT and PARAFAC at low SNR regions. However, all the algorithms provide very similar RMSE performance at high SNR regions (SNR is larger than 10 dB).Example 3.
The RMSE performance versusL. In this simulation, we consider β=−1 and α=0.8, and SNR is fixed at 0 dB. The RMSE curves are illustrated in Figure 5. It seems that RMSE can be slightly improved with L increasing. The ESPRIT method provides the worst RMSE performance. Since the tensor nature has been exploited, the PARAFAC estimator offers better RMSE than ESPRIT. However, the proposed method outperforms all the compared methods since it is free from the spatially colored noise.Example 4.
The RMSE performance versus colored parameterβ. In this simulation, we consider α=0 and L=1000 and SNR is set to 0 dB. Figure 6 gives the RMSE curves. It should be pointed that when β≤0.2, the noise can be regarded as spatially uncorrelated (or with very weak relevance). Under such condition, all the algorithms offer very close RMSE performance. However, RMSE associated with ESPRIT and PARAFAC become worse when β>0.2, since they are sensitive to the spatially colored noise. By comparison, the proposed algorithm offers much better RMSE than the compared methods. It is evident that the proposed method is insensitive to the spatially colored noise.Figure 2
Scatter results of the proposed algorithm: (a) 2D-DOA estimation, (b) polarization parameter estimation.
(a)(b)Figure 3
Average RMSE performance comparison versus SNR with Gaussian white noise: (a) average RMSE on 2D-DOA estimation; (b) average RMSE on polarization parameter estimation.
(a)(b)Figure 4
Average RMSE performance comparison versus SNR with colored noise: (a) average RMSE on 2D-DOA estimation; (b) average RMSE on polarization parameter estimation.
(a)(b)Figure 5
Average RMSE performance comparison versusL with colored noise: (a) average RMSE on 2D-DOA estimation; (b) average RMSE on polarization parameter estimation.
(a)(b)Figure 6
Average RMSE performance comparison versusβ: (a) average RMSE on 2D-DOA estimation; (b) average RMSE on polarization parameter estimation.
(a)(b)


## 6. Conclusion

In this paper, a novel HOSVD algorithm is proposed for distributed EVS array, the core of which is to estimate 2D-DOA by the vector cross product of the normalized Poynting vector, which is achieved from cross-correlation tensor signal subspace. The proposed algorithm is better than the state-of-the-art algorithms since it is insensitive to the spatially colored noise and sensor position. Numerical simulation results verify the effectiveness of the proposed algorithm. It should have a bright future in radar, sonar, and wireless communication as well as Internet of Things.



---

*Source: 9973253-2021-05-11.xml*