Source Localization Using Distributed Electromagnetic Vector Sensors

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.


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][2][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][8][9], sparsity-aware strategy [10,11], and tensor approach [12][13][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 tech-nique 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 highdimensional 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][22][23][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][33][34][35][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 higherorder 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; I N is a N × N identity matrix; ⊗ , ⊙ , and * denote the Kronecker product, the Khatri-Rao product, and the vector cross product, respectively; the superscript ðXÞ T , ðXÞ H , ðXÞ −1 , and ðXÞ † indicate the operations of transpose, Hermitian transpose, inverse, and pseudoinverse, respectively; D n fAg denotes a diagonal matrix with the diagonal entities which are the nth row of A; aðmÞ accounts for the mth entity of the vector a; k⋅k F accounts for the Frobenius norm; phaseð•Þ is to get the phase; and Ef⋅g returns the mathematical expectation.

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.
The ði 1 ,⋯,i n Þth entity of X is the ði n , jÞth element of ½X n , j Definition 2 (mode-n tensor-matrix product). The mode-n product of X ∈ C I 1 ×I 2 ×⋯I N by a matrix A ∈ C J n ×I n , is denoted by Definition 3 (important principles of the mode-n tensormatrix product). The following principles are important for mode-n tensor-matrix product: Definition 4 (HOSVD). The HOSVD of an N-order tensor X ∈ ℂ I 1 ×⋯I N is denoted by where G ∈ ℂ J 1 ×⋯J N is the so-called core tensor and U ðnÞ ∈ ℂ I n ×⋯J N is the unitary matrix. In mode-n tensor unfolding format, (2) is rewritten as 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 parameterdepend only vector. An important characteristic is that e and h are mutually orthogonal to each other. Moreover, the normalized Poynting-vector fulfills [20] u v w 2 6 6 4 2.3. Problem Formulation. In this paper, we consider an M -element EVS array with arbitrary geometry, as illustrated in Figure 1. Suppose that the mðm = 1, 2, ⋯, MÞth EVS is located at coordinate ðx m , y m , z m Þ. 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] y where a k ≜ ½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 = x m sin θ k cos ϕ k + y m sin θ k sin ϕ k + z m cos θ k . b k is the associated polarization steering vector, s k ðtÞ is the kth source signal, and sðtÞ = ½s 1 ðtÞ, s 2 ðtÞ,⋯,s K ðtÞ T is the source signal vector. nðtÞ is the array noise, which is assumed to be spatially correlated, temporally independent Gaussian process, i.e., where Q is the covariance matrix of the noise, which is a Her-mitian matrix. In the presence of Gaussian white noise, Q = σ 2 I.

The Proposed Algorithm
3.1. Spatially Colored Noise Suppression. Suppose the source signals are uncorrelated and sðtÞ is uncorrelated with nðtÞ; then, the covariance of yðtÞ is given by where R s is the covariance matrix of sðtÞ. 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: where Λ = Efsðt 1 Þs H ðt 2 Þg 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 viâ Obviously, the noise is eliminated in R. In practice, R can be approximated by its truncated eigenvalue decomposition as where Σ d is the dominate eigenvalue matrix, U s consist of the K dominate eigenvectors.

HOSVD.
Actually, R can be arranged as a fourth-order tensor as [35] R m, n, p, q The HOSVD of R is given by where G ∈ ℂ M×6×M×6 is the core tensor, U 1 ∈ ℂ M×M , U 2 ∈ ℂ 6×6 , U 3 ∈ ℂ M×M , and U 4 ∈ ℂ 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 where U 1s , U 2s , U 3s , and U 4s consist of the eigenvectors associated to the K dominate eigenvalues. G s is called the signal component of G, which is given by Inserting G s into (14) yields By Hermitian unfolding R a , we can construct an enhanced covariance matrix as Combined with the result of (11), we can obtain Since R is Hermitian, we have U 1s = U * 3s and U 2s = U * 4s . Performing eigendecomposition on R a , we can obtain a signal subspace as Since the virtual response matrix A ⊙ B spans the same subspace with E s , we have where T ∈ ℂ K×K is a nonsingular matrix. . We can get the following relationship: with where β 1,q k . Define the following selection matrices: where i 6,q is the qth row of I 6 . Therefore, the relationships in (21) becomes Inserting (24) into (20) establishes Equivalently, According to the first line of (26), we can perform eigendecomposition of ðJ 2 E s Þ † J 1 E s , 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.

Wireless Communications and Mobile Computing
In fact, b k can be written in the following normalized format: k ; then, we can get the estimation of u k , v k , and w k , which are denoted byû k ,v k , and w k , respectively. Then, 2D-DOA estimation can be achieved via Obviously, the estimated 2D-DOA are automatically paired.  Table 1: Algorithmic steps of the proposed algorithm.

Polarization Parameter Estimation. It is easy to find
Step No. Operation Step 1 Estimate the noiseless covariance matrix R through (10) Step 2 Rearrange the covariance measurement into a fourth-order tensor R via (12), and perform HOSVD on R Step 3 Construct R a from (17), and then, get R a via the Hermitian unfolding of R a Step 4 Perform eigendecomposition on R a to get the signal subspace E s Step 5 Calculate J 2 E s ð Þ † J 1 E s , and perform eigendecomposition on it to get Φ 1,2 ð Þ and T Step 6 To get Φ 1,3 ð Þ , Φ 1,4 ð Þ , Φ 1,5 ð Þ , and Φ 1,6 ð Þ via (26) Step 7 Estimate u k v k w k ½ T via (30), and then, get the 2D-DOA through (31) Step 8 EstimateB m via (37), and computeĝ k through (38). Finally, get the polarization parameters via (39)

Wireless Communications and Mobile Computing
According to (4), we should havẽ where D k and g k are the direction-only and polarization-only variables. Define Then,B m can also be expressed as Therefore,B m can be estimated viâ the kth column of which is denoted byb k,m . After the 2D-DOA have been estimated, we can construct the matrixD k according to (4); then, we can compute the following least squares:ĝ Obviously,ĝ k is the estimation of e −jπτ m,2 /λ g k . Finally, the Also, the polarization parameters are paired automatically.

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 Table 1.

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 Of6 3 M 3 g. The main computational burden of the PARAFAC estimator in [28] is the iteration, which requires OfMK 2 + 6K 2 + LK 2 g 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

4.2.
Identifiability. The identifiability of the proposed approach is equal to the maximum K, 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.

Simulation Results
Numerical simulations are carried out to show the superiority of the proposed algorithm. In the following simulations, we assumed an EVS array with M = 8 elements, the coordinates of the EVSs ðx m , y m , z m Þ fulfill uniform distribution in ½−2λ, 2λ. Assume that there are K = 3 uncorrelated farfield 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 ω = ½w 1 , w 2 ,⋯, ω K T is defined as whereŵ t,k is the estimation of w k 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 versus L. 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.

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 stateof-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.

Data Availability
There is no available data for this paper.

Conflicts of Interest
The authors declare no conflict of interest.