A Type-2 Block-Component-Decomposition Based 2D AOA Estimation Algorithm for an Electromagnetic Vector Sensor Array

This paper investigates a two-dimensional angle of arrival (2D AOA) estimation algorithm for the electromagnetic vector sensor (EMVS) array based on Type-2 block component decomposition (BCD) tensor modeling. Such a tensor decomposition method can take full advantage of the multidimensional structural information of electromagnetic signals to accomplish blind estimation for array parameters with higher resolution. However, existing tensor decomposition methods encounter many restrictions in applications of the EMVS array, such as the strict requirement for uniqueness conditions of decomposition, the inability to handle partially-polarized signals, etc. To solve these problems, this paper investigates tensor modeling for partially-polarized signals of an L-shaped EMVS array. The 2D AOA estimation algorithm based on rank-(L1,L2,·) BCD is developed, and the uniqueness condition of decomposition is analyzed. By means of the estimated steering matrix, the proposed algorithm can automatically achieve angle pair-matching. Numerical experiments demonstrate that the present algorithm has the advantages of both accuracy and robustness of parameter estimation. Even under the conditions of lower SNR, small angular separation and limited snapshots, the proposed algorithm still possesses better performance than subspace methods and the canonical polyadic decomposition (CPD) method.


Introduction
Array signal processing is an important branch of the information acquisition and detection field, which has been studied and applied extensively in the academic and industrial communities for nearly half a century [1][2][3]. The traditional array usually consists of scalar sensors, which are able to receive a one-dimensional signal only. However, electromagnetic waves have diverse magnetic and electric field components. How to make the best of the complete information of the electromagnetic wave has been an important subject that has been paid close attention by researchers. The EMVS is a kind of polarization-sensitive antenna, which can receive magnetic and electric field information of the signal simultaneously [4]. Compared with the traditional scalar sensor, EMVS can sense electromagnetic waves with different polarization forms and fully exploit information carried by the wave itself. Since the 1980s, many investigations have been presented to study the problem the signal processing community hopes to find a new framework that not only can maintain the uniqueness of decomposition, but also that does not need the severe requirement for rank constraints. In 2008, L. De Lathauwer proposed a tensor decomposition model combined with the features of CPD and Tucker decomposition (TKD), called the block component decomposition (BCD), or the block term decompositions (BTD) [33][34][35]. This model has various forms of decomposition. For instance, a rank-(L, L, 1) BCD expresses the factors as the outer product of one rank-L matrix with a vector. Such a model has been effectively applied in blind source separation [36,37], in which a rank-(2, 2, 1) BCD model is utilized to match data with limited samples. When the reflections result in the inter-symbol-interference (ISI) occurring simultaneously on both the near and far field, a more flexible block-CPD model can be applied [38]. This shows that the BCD model has more flexible applicability than the CPD. Another one is rank-(L 1 , L 2 , ·) BCD, which expresses the factors as the product of a core tensor with two factor matrices. In 2008, D. Nion and L. De Lathauwer introduced such a model into the DS-CDMA (Direct Sequence-Code Division Multiple Access) system to construct a blind receiver used for convolutive mixtures. It can achieve a good performance on both channel identification and symbol estimation [39]. In 2013, combining the partial least squares (PLS) with BCD modeling, Q.B. Zhao et al. put forward a higher order PLS method that can achieve good results in 3D movement trajectories decoding for the electrocorticogram (ECoG) signal [40]. In 2014, B. Hunyadi et al. introduced the BCD into electroencephalogram signal processing [41], in which a better performance is gained in the ictal pattern identification than the CPD algorithm. In addition, these types of BCD models also include the PARALIND (PARAllel factors with LINear Dependency) [42,43], which has mainly been utilized in blind beamforming and multi-antenna coding [44]. In 2012, M. Sorensen and L. De Lathauwer proposed a rank-(L 1 , L 2 , ·) BCD algorithm with Vandermonde factors, which has good estimation performance in the uniform linear array (ULA) [45]. In 2013, L. Sorber et al. generalized the CPD and rank-(L r , L r , 1) BCD method to propose a optimization-based tensor decomposition algorithm that can diminish computation source requirement by means of Jacobian's Gramian structure [46]. In 2015, P. Tichavsky et al. proposed a non-orthogonal tensor diagonalization BCD algorithm, which has a nearly quadratic convergence performance [47]. In 2016, X.F. Gong et al. put forward two kinds of coupled block simultaneous generalized Schur decomposition (SGSD)-based BCD algorithms, which have advantages of low computational resource requirement and high computation accuracy [48]. So far, there has been little public research on the application of the BCD in the vector sensor array. As discussed above, the received signal model of the EMVS array has a natural multidimensional structure. Additionally, for 2D AOA, the angle pair-matching is a troublesome problem. The earlier research works often need complicated searching methods to achieve pair-matching [49][50][51][52]. Some researchers put forward the automatic pairing method by using the eigenvalue of the transfer matrix, but the issue of eigenvalue ordering has not been well solved with such an algorithm [53]. Since the BCD method possesses the blind estimation feature, it can solve such a problem properly. This paper takes this as the research motivation to investigate the parameter estimation algorithm based on the BCD modeling for the EMVS array.
The key contribution of this paper is the rank-(L 1 , L 2 , ·) BCD modeling of the partial polarization signal incident on an EMVS array. For achieving the AOA parameters estimation, an iterative algorithm is proposed to compute each factor of such a BCD model. This algorithm adopts a block diagonal constraint criterion of the steering matrix and constructs a subspace-like estimator to achieve the estimation of spatial frequencies. The BCD-based methods have a common feature that the performance of the algorithm can be guaranteed only if the rank conditions are satisfied. Thus, we analyze whether or not the rank conditions can be met. Additionally, for the partial polarization signal, the decomposition uniqueness of rank-(L 1 , L 2 , ·) BCD is investigated in this paper. The proposed method can accomplish the pair-matching of azimuth and elevation angles automatically by means of the decomposition uniqueness of the BCD model.
The other sections of this paper are organized as follows: Section 2 describes the received signal model of the EMVS array and introduces the BCD method; Section 3 develops the algorithm based on rank-(L 1 , L 2 , ·) BCD for 2D AOA estimation; Section 4 presents numerical simulations for verifying the proposed algorithm and demonstrates the comparison with the existing methods; finally, the last section concludes this paper.

Data Model
For convenience of statement hereafter, several notions are explained as follows. The scalar is indicated by italic letters, e.g., x. The vector and matrix are indicated by bold italic lower-and upper-case letters, e.g., x, X. The tensor is indicated by calligraphic capital letter, e.g., X . X F is the Frobenius norm of the tensor X , expressed as the square root of the sum of the absolute squares of its elements. The operator ⊗ denotes the Kronecker product between matrices. The operator indicates the column-wise Kronecker product, or the Khatri-Rao product [16]. The operator × n denotes the mode-n product between the tensor and matrix, where the mode is considered as the order of the tensor data. For a third-order tensor X ∈ C M×N×P , X (1) is defined as the the unfolded matrix, or matricization of X in mode-1, i.e., X (1) = [X 1 , X 2 , ..., X P ] ∈ C M×(PN) ; similarly, X (2) and X (3) denote the matricization in mode-2 and mode-3, respectively. The mode-n rank of tensor X is the rank of its unfolded matrix X (n) .

EMVS Array Signal Model
Under the assumption of narrowband and far-field signal [7,8], the measurement of EMVS can be expressed as: where φ k , θ k denote the azimuth and elevation angle of arrival signal, respectively, as shown in Figure 1. For a fully-polarized signal, ξ k (t) = Q k w k s k (t) = g k s k (t), in which s k (t) is the k-th source signal, Q k and w k are the variables associated to polarization parameters.
where µ k ∈ [−π/2, π/2] is the polarization orientation angle of the electromagnetic wave and ν k ∈ [−π/4, π/4] is the ellipticity angle, as shown in Figure 2. Thus, y k (t) = Θ k g k s k (t) = b k s k (t), where b k ∈ C 6×1 is referred to as the polarization steering vector corresponding to the k-th signal. The above signal model is only for single EMVS at the receiver. For an EMVS array, the manifold configuration of sensors must be dealt with. This paper considers the L-shaped array [54,55] shown in Figure 3, which consists of a couple of orthogonal uniform linear arrays, and each of them has M sensors. Given the polarization steering matrix as B = [b 1 , . . . , b K ] ∈ C 6×K and considering noise, the received signal model of the EMVS array can be expressed as: where A x = [a x,1 , a x,2 , ..., a x,K ] ∈ C M×K is the steering matrix of the subarray along the x-axis. Without loss of generality, assuming the array element spatial interval as d and letting α k = − 2πd λ cos φ k , in which λ is the wavelength of source signal, we have a x,k = 1, e jα k , · · · , e j(M−1)α k T . s (t) ∈ C 1×K is the source signal, vector and n x (t) ∈ C M×1 is the additive white Gaussian noise (AWGN) vector. The signal model of the subarray along the z-axis is consistent with the form of the x-axis.  In practical applications, the fully-polarized signal appears rarely. Most signals are partially polarized. The trajectory of the endpoint of the electric field vector at the spatial frequency is no longer a certain ellipse, but a curve of shape and direction that changes over time. In this case, the polarization parameters are not constant. Thus, g k cannot be incorporated into the polarization steering vector b k . For instance, the circularly-polarized signal has ξ k (t) = Q k ws (k) 2 (t) are the k-th signal complex envelope [8]. For such a partially-polarized signal, the received signal model of the EMVS array can be rewritten as: where b denotes the Khatri-Rao product with block form [33]; is the vector of partially-polarized signals. The signal model of the subarray along the z-axis has the same form as that along the x-axis. Without ambiguity, the argument t is omitted in the following discussions for simplicity.
For N snapshots, the received signal models (5) and (6) become: where S ∈ C N×2K is the source signal matrix for the N snapshots and N x , N z ∈ C M×N are the matrices of AWGN on the subarrays along the xand z-axis, respectively.

Block Component Decomposition
As the combination of CPD and TKD, the BCD framework maintains the advantages of both, the former with decomposition uniqueness and the latter with the relaxing rank requirement for the factors [33][34][35]. In addition, BCD has the feature of structural diversity, such as rank-(L, L, 1), rank-(L 1 , L 2 , L 3 ) and rank-(L 1 , L 2 , ·) patterns. Different patterns have different forms of decomposition, corresponding to different signal models. This paper mainly pays attention to the rank-(L 1 , L 2 , ·) BCD, i.e., the Type-2 BCD model.
The decomposition of a tensor X ∈ C M×P×N is called the Type-2 block component decomposition, if it can be expressed as the summation of a series rank-(L 1 , L 2 , ·) factors, where C k ∈ C L 1 ×L 2 ×N can be termed as the core tensor of the k-th component. The mode-1 rank and mode-2 rank of C k are L 1 and L 2 , respectively. Both A k ∈ C M×L 1 and B k ∈ C P×L 2 are of full column rank. The type-2 BCD model in rank-(L 1 , L 2 , ·) terms is schematized in The BCD model can also be expressed as: where and Equation (10) can be written as X = A, B, C (cf. [16]).

The BCD Modeling of the EMVS Array
Consider the observation data on both subarrays [8]; we have the received signal model, The corresponding matrix form is: Then, Equation (13) can be rewritten as: where Y ∈ C (ML 1 )×6×N is the received signal; S k ∈ C L 1 ×L 2 ×N is the source signal, of which the mode-3 matricization is S T k , i.e., S T k = [S k ] (3) . N is the corresponding AWGN. Next, we investigate if the rank condition of each factor is satisfied in Equation (14). First, the steering matrix {A k } K k=1 ∈ C ML 1 ×L 1 is of full column rank when M ≥ 1. The same is true , the size of mode-1 matricization is (L 1 × NL 2 ) and (L 2 × NL 1 ) for mode-2. Thus, the mode-1 and mode-2 rank of {S k } K k=1 are L 1 and L 2 , respectively, when snapshot N is sufficiently large. Consequently, the rank conditions of each decomposition factor in Equation (14) can be satisfied if: For the L-shaped array in this paper, on each subarray, the number of sensors M > 1 must be true. Additionally, since L 1 = L 2 = 2, no matter what the snapshots are, multiple or single, Equation (14) will be valid according to Equation (15).

Uniqueness Analysis
Based on the EMVS array signal model, the decomposition uniqueness of Equation (14) is investigated in this section.

Theorem 1. Uniqueness of BCD modeling for the EMVS array:
The rank-(L 1 , L 2 , ·) BCD descried in Equation (14) is essentially unique, if A and Ψ are of full column rank, N ≥ 3 and all elements of {S k } K k=1 are characterized by jointly continuous probability density functions.
For the proof of Theorem 1, reference can be made to [34]. Since the EMVS model in this paper is based on the assumption that source signals obey the zero mean Gaussian random process [7], this implies that the elements of {S k } K k=1 are surely characterized by jointly continuous probability density functions, i.e., the conditions of factors {S k } K k=1 in Equation (14) are satisfied from Theorem 1. In addition, the array manifold based on the L-shaped array in this paper gives L 1 = L 2 = 2; the structure of {A k } K k=1 is block diagonal; and each block factor is the steering vector of the subarray, which has the Vandermonde structure. As a consequence, the arbitrary two columns of A are linearly independent. Hence, rank (A) = min (2M, 2K). Thus, when M > K, A is of full column rank. From the definition of {Θ k } K k=1 in Section 2.1, we have rank (Ψ) = min (6, 2K). Thus, when K < 3, Ψ is of full column rank. The above analysis reveals that the decomposition from the tensor model (14) of the EMVS array is unique, if the following conditions are met, It should be noted that Equation (16) is only the sufficient condition. In practical applications, if the constraints are imposed on factor matrices, the decomposition (14) can still be unique, even if the conditions limited by Equation (16) may not be satisfied.

AOA Estimation Algorithm
Given the rank-(L 1 , L 2 , ·) BCD model (14) with additive noise, the task is to acquire the estimation of factor matrices {A k } K k=1 , {Θ k } K k=1 from the received signal Y ∈ C (2M)×6×N and obtain the estimation of spatial frequency {α k } K k=1 . This paper adopts the minimum mean square error (MMSE) criterion [35], With the mode-n matricizations of the tensor Y and N , i.e., Y (n) and N (n) , n = 1, 2, 3, the signal model (14) can be rewritten as: Several algorithms can be employed to solve the above problems [35,46,56,57]. The typical one is alternating least squares (ALS) [19,46]. With two fixed components of A, Ψ, S in each iteration, the remaining one is updated by the least squares approach. Then, the same scheme is repeated in turn until the convergence condition is satisfied.
For instance, with {S k } K k=1 and {Θ k } K k=1 fixed, the update of A is calculated as follows, then we can obtain:Â where (·) † represents the pseudo-inverse of a matrix. With {S k } K k=1 and {A k } K k=1 fixed, the same procedure gives the update of Ψ,Ψ Then, fixing {A k } K k=1 and {Θ k } K k=1 leads to the update of S, During the steps of updating A and Ψ, the normalization by using QR decomposition can be added to enhance the stability of estimation [35]. It is necessary to guarantee that F T 1 ,F whereS k = [S k ] :,:,2 · [S k ] † :,:,1 is of full rank, known from the previous rank condition of rank-(L 1 , L 2 , ·) BCD and the characteristics of the source signal. A is also full column rank in terms of the uniqueness analysis because Equation (16) is fulfilled. This means that the column space of Y 2 Y † 1 is the same as that of A. A similar conclusion can be deduced that Y T 2 Y T 1 † and Ψ have the same column space. Therefore, we can employ singular value decomposition (SVD) to initiate the factor matrices for the ALS algorithm.
The analysis of Section 3.1 shows that solving Equation (14) can result in the estimation of A k , Θ k , k = 1, . . . K. After finding the steering matrixÂ k by means of ALS, the spacial frequency can be obtained. Note that the trivial indeterminacies ofÂ k exists, i.e., where Γ k ∈ C L 1 ×L 1 is a nonsingular matrix. Equation (26) indicates that the BCD estimation result and steering matrix have the same column space. Because A k has the block diagonal structure, we can extract the correct AOA information by exploiting this feature.
For the array manifold configuration in this paper, the A k can be parted as In the same way, the estimation ofÂ k is parted aŝ . Then, combined with Equation (26), we have: The above equations indicate thatÂ k,m is a rank-1 matrix and has the same column space as the steering vector. Hence, the spacial frequency can be estimated by subspace methods. Giving the SVD ofÂ k,m ,Â the left singular vector u k,1 corresponding to the maximum singular value can be obtained. Take the first and last (M − 1) rows of u k,1 , denoted as u a and u b , to construct U ab = [u a , u b ]. Applying the eigenvalue decomposition of U H ab U ab , in which: yields the estimation of spacial frequencyα k = ∠ − w 1,2 w 2,2 , where ∠ (·) denotes the phase angle of a complex variable. Finally, the estimation of the arrival angle on the array along the x-axis can be obtained,φ The same method can be used to get the estimation of the elevation angleθ k . After obtaining the estimation of φ k and θ k , the azimuth angle φ k can be estimated directly, SinceÂ k includes the steering vector information on both subarrays and the Type-2 BCD is unique, the estimation of azimuth and elevation angles accomplish the matching pairing automatically. It should be noted that even if the array manifold is arbitrarily configured and no longer has the Vandermonde structure, the spatial frequency can still be estimated by utilizing the one-dimensional spectral peak search method.
The algorithm flow is shown in Algorithm 1. The ALS procedure is the key step to estimate the steering matrix, which has the complexity of O 2K (P + 1) ∏ P p=1 I p . Such a procedure consists of three parts. The first is to update the factor matrices, which has the complexity of O 2KP ∏ P p=1 I p . The second is to solve the linear systems, which has the complexity of O 1 3 K 3 P + 2K 2 ∑ P p=1 I p . The third is function evaluation, which has the complexity of O 2K ∏ P p=1 I p . In addition, we employ the SVD to calculate the spatial frequencies, which has the the complexity of O 2L 3 1 . For the signal model of the EMVS array in this paper, we have P = 3, I 1 = 2M, I 2 = 6, I 3 = N. It can be concluded that the total computational complexity is O 2K (3 + 1) × (2M × 6 × N) + 2L 3 1 = O 96KMN + 2L 3 1 .

Algorithm 1 2D AOA estimation algorithm for the EMVS array based on Type-2 BCD.
Input: the received signal Y ∈ C 2M×6×N ; the number of source signals, K; BCD model parameters, L 1 = L 2 = 2; the threshold of error, .

Numerical Simulations and Discussion
In this section, several simulation experiments are demonstrated for proving the validation of the proposed method. The experimental scenario is consistent with the received signal model (14) described in Section 3.1. It is assumed that the L-shaped array consists of two ULAs, and each one has M sensors. The number of source signals is K; the number of snapshots is N; the Monte Carlo independent trials is L.
Since the signal studied in this paper is partially-polarized waves and the polarization state is time-varying, so the simulations in this section only estimate the AOA of the signal. In the numerical experiments, the polarization orientation angle is set as a random variable obeying the zero mean Gaussian distribution; the simulation results are given in the following subsections. If the distribution function changes to obeying the uniform distribution, the corresponding simulation results would have considerable consistency with the Gaussian distribution.

Different SNRs
The simulation parameters are configured as: M = 5, K = 2, N = 200; the distance between adjacent sensors is λ/2; two groups of the AOAs are (φ 1 , θ 1 ) = (100 • , 65 • ) and (φ 2 , θ 2 ) = (120 • , 45 • ); the SNR range is −3 dB∼15 dB; the Monte Carlo trials L = 500. The azimuth-elevation angles estimation is shown in Figure 5.  For comparing the accuracy of estimation under different SNRs more directly, we analyze the simulation result for one group of AOA, (φ 2 , θ 2 ), as shown in Figure 6. The result indicates that the proposed algorithm has better estimation accuracy under different noise conditions. Even if the SNR is low, different AOAs can be distinguished, which can be verified by the result of detection probability in Section 4.2.2. When the SNR is −3 dB, the detection probability of the proposed algorithm is larger than 85%. However, the value of the existing method, e.g., ESPRIT, is less than 50% only.

Different Angular Separations
This group of experiments investigates the estimation of AOA under different angular separations for fixed SNR, which is shown in Figure 7. Set SNR as 5 dB, and the remaining parameters agree with those of Section 4.1.1. Two groups of AOA are selected as (φ 1 , From the scatter diagram, the AOA estimation can be fully resolved when the angular separation is large than 3 • . Even if the separation degree is very small, (φ 2 , θ 2 ) can still be estimated correctly with high probability, as shown in Figure 8. This conclusion can also be proven by the result of detection probability analysis in Section 4.2.2. When ∆ = 1 • , the detection probability of the proposed algorithm is larger than 80%. However, the ESPRIT has only about 40%.

Different Snapshots
This group of experiments investigates the estimation of AOA under different snapshots, which is set as N ∈ [5, 500]. The SNR is fixed at 5 dB, and the remaining simulation parameters are the same as given in Section 4.1.1. With the given Monte Carlo trials L = 500, the simulation result is shown in Figure 9. Even if the snapshots are finite, the simulation result shows that the proposed algorithm can achieve better angular resolution. For the different snapshots, the estimation of the two AOA groups appears to be hardly aliasing. When the snapshot is five, the detection probability with the proposed algorithm is more than 95% from the result in Section 4.2.2, while the detection probability of ESPRIT is less than 75%. One group of angle estimation results is shown in Figure 10. The simulation result reveals that the estimation accuracy of (φ 2 , θ 2 ) rises with the snapshot increasing.

Performance Comparison
The comparison of the proposed algorithm with three existing algorithms has been investigated in this section, which shows the performance advantage of it in both root mean square error (RMSE) and detection probability. The reference methods include subspace-like algorithms, i.e., TM-MUSIC [27], ESPRIT [58,59], and the CPD-based 2D AOA estimation algorithm [26]. The ESPRIT algorithm is based on matrix modeling, the TM-MUSIC and CPD algorithms are based on tensor modeling.

RMSE
First, we investigate the RMSE of each method, which is defined as: where φ k (l) ,θ k (l) are the estimation values of the k-th AOA for the l-th Monte Carlo trial.
The Cramer-Rao lower bound is defined as where CRB (φ k ) and CRB (θ k ) are the diagonal elements of CRB matrix [60,61], corresponding to the parameters φ k and θ k , respectively.
(1) RMSE versus different SNRs: The parameter settings are the same as those in the Section 4.1.1. The results are shown in Figure 11. The RMSE curves interpret the advantage of the algorithm developed in this paper over the reference methods in the whole SNR range. The three tensor-based algorithms are better than the matrix-based algorithm, i.e., ESPRIT, since tensor modeling can make the best of the multidimensional structure information in the received signal. It is noticed that when the SNR is high, e.g., 15 dB, the performance of TM-MUSIC is close to the CPD algorithm. Compared with the CPD algorithm, as analyzed in Section 3.2, the BCD algorithm can solve the steering matrix A k directly and then automatically achieve the matching pair of the azimuth and elevation angles. Thus, the detection probability of BCD is higher than CPD, which utilizes the pair-matching technique based on the cross-correlation matrix (CCM) [26]. The results in Section 4.2.2 provide the same conclusion. Therefore, the performance of the proposed algorithm is better than that of the CPD algorithm. As a classic matrix subspace algorithm, the RMSE of the ESPRIT algorithm is the highest. This may be because ESPRIT has to lose a part of array aperture for constructing the rotation-invariant matrix. (2) RMSE versus different angular separations: The parameter settings are taken as those in the Section 4.1.2. Figure 12 shows the simulation results under different angular separations. The RMSE curves vary with the steeper slope as the separation degree is less than 5 • , which reveals the sensitivity to angular separation within such an interval. The results of detection probability analysis in Section 4.2.2 can also verify this conclusion. The curves in Figure 15 indicate that the detection probability is larger than 95% for all algorithms with the angular separation larger than 5 • . Below this value, the detection probability is diverse for different algorithms. It is noticed that when the angular separation is 1 • , of the three reference methods, the performance of TM-MUSIC algorithm approaches the CPD. When the angular separation is increasing, the RMSE of TM-MUSIC is becoming larger than the tensor-decomposition-based algorithms and so is the ESPRIT algorithm. Essentially, the TM-MUSIC is a kind of smoothing approach along a certain mode. Therefore, it can achieve a better performance than the classic subspace method, which may be influenced by outliers when the parameters to be estimated are close. Of all of the methods, the proposed algorithm has the lowest RMSE, which shows its high-resolution feature under the condition of small angular separation. (3) RMSE versus different snapshots: The experimental parameters are set as those in the Section 4.1.3. The simulation results are depicted in Figure 13. The RMSE curves indicate that the tensor-based methods outperform the classic matrix-based subspace method ESPRIT during the entire snapshots' range. The RMSE of TM-MUSIC approximates that of CPD algorithm when the snapshots are larger than 300. The analysis in Section 3.1.1 indicates that when K = 2, the uniqueness of BCD can be met as long as snapshots N > 3. The similar conclusion is also suitable for the CPD algorithm. The results of detection probability analysis in Section 4.2.2 also provide the same conclusion that even for limited snapshots, BCD modeling can keep the detection probability on a higher level. Consequently, the performance of the proposed algorithm is better than the reference algorithms.

Detection Probability
The detection probability measures the success rate of pair-matching for two group of azimuth-elevation angles. In this section, detection probabilities of the four methods are compared in accordance with three simulation experiments, grouped by SNR, angular separation and snapshot.
(1) Detection probability versus different SNRs: The experiment parameters as set in the Section 4.1.1 yield the simulation results shown in Figure 14. The detection probability curves show the obvious advantage of the proposed algorithm in the success rate of angle pair-matching at low SNR. For example, when SNR = −3 dB, the detection probability of the proposed algorithm reaches 88% and increased by 69.2%, 49.1% and 35.4% compared with ESPRIT (52%), TM-MUSIC (59%) and CPD (65%), respectively. This benefit by automatic angle pairing resulted from estimating steering vectors and the decomposition uniqueness of BCD, which makes the proposed algorithm obtain more powerful detection ability and possess the best estimation performance under the worse noise environment. (2) Detection probability versus different angular separations: The experimental parameters are the same as those in Section 4.1.2. Figure 15 shows the simulation results. We can find that for different angular separations, the algorithm developed in this paper always has the highest success rate of angle pair-matching, especially at a small angle separation degree (∆ ≤ 3 • ). For instance, when ∆ = 1 • , the success rate of the proposed algorithm is 84% and increased by about 82.6%, 52.7% and 25.4% compared with ESPRIT (46%), TM-MUSIC (55%) and CPD (67%), respectively. When ∆ ∈ (3 • , 5 • ), the detection probabilities of all methods have significant improvement. Meanwhile, the proposed algorithm is still optimal. When ∆ > 5 • , all of the methods achieve 100% detection probability. (3) Detection probability versus different snapshots: The experimental parameters are selected as set in Section 4.1.3. The simulation results can be seen in Figure 16. The detection probability curves show that the proposed algorithm has a higher success rate of angle pair-matching versus different snapshots, especially at small snapshots. As an example, when the snapshot is five, the detection probability of the proposed algorithm reaches to 97% and increased by about 31.1%, 15.4% and 7.8% compared with ESPRIT (74%), TM-MUSIC (84%) and CPD (90%), respectively. What is more, one more interesting characteristic different from the reference methods is that the detection probability of the proposed algorithm is slightly influenced by snapshots, which is embodied by the evenly-varied curve. This characteristic makes the proposed algorithm possess the highest detection probability even under the limited snapshots. Hence, the algorithm is more suitable for the scenario with the real-time requirement.

Conclusions
The application of BCD tensor modeling in parameter estimation of the EMVS array has been studied in this paper. In allusion to the issue of lacking effective methods for solving parameter estimation of partially-polarized waves, we put forward an algorithm for 2D AOA estimation of the EMVS array based on the rank-(L 1 , L 2 , ·) BCD model. This algorithm can make full use of the multiway structural information of the received signal and achieve a high-resolution estimation. Benefiting from the decomposition uniqueness of such a Type-2 BCD model, the pair-matching of azimuth and elevation angles can be automatically achieved. According to the possible situations in practical applications, we carry out several groups of numerical experiments to verify the effectiveness of the proposed algorithm. The experimental results show that both the estimation accuracy and detection success rate of this algorithm are superior to the existing subspace method based on matrix decomposition and the CPD method based on tensor decomposition. The proposed algorithm manifests robust and good performance under severe conditions, such as low SNR, small angular separation and limited snapshots. This makes the algorithm possess potential application value in practical engineering applications.
Author Contributions: Yu-Fei Gao conceived of and verified the feasibility of tensor modeling for the EMVS array and proposed the 2D AOA estimation algorithm based on rank-(L 1 , L 2 , ·) BCD. Yu-Fei Gao, Wei Xie and Yan-Bin Zou performed the experiments. Guan Gui and Qun Wan contributed the design of numerical simulations. Yu-Fei Gao wrote this paper. Wei Xie, Yan-Bin Zou, Yue Yang and Guan Gui checked the manuscript and contributed to the rearrangement of the materials.

Conflicts of Interest:
The authors declare no conflict of interest.