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- ( L 1 , L 2 , · ) 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.

1. 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 of parameter estimation in EMVS arrays. In 1981, R.T. Compton put forward an adaptive signal processing method of the tripole antenna array, which can obtain better anti-disturbance performance than traditional array processing methods [5]. In 1992, J. Li et al. introduced the rotational invariance techniques (ESPRIT) algorithm to the signal processing of the EMVS array to present a joint estimation algorithm of the AOA-polarization parameter [6]. In 1994, A. Nehorai and E. Paldi studied the received signal model of the EMVS array and CRB analysis of the parameter estimation comprehensively, in which the mean square angular error (MSAE) and the covariance vector angular error (CVAE) were introduced [7]. In 1996, K.C. Tan et al. analyzed the linear independence of steering vectors of EMVS for different polarization signal types [8]. In 2000, K.T. Wong and M.D. Zoltowski developed an algorithm of AOA and polarization parameter estimation with a closed-form solution for the EMVS array, by which high precision estimation could be achieved [9]. In 2004, D. Rahamim et al. proposed a preprocessing algorithm of coherence signals over the polarization domain, which can be applied to the EMVS array arbitrarily distributed in space [10]. In 2012, X. Yuan constructed a joint estimation algorithm of AOA-polarization parameters for the polynomial-phase signal with full polarization and derived the corresponding analytical expression of CRB [11]. It can be seen that EMVS array signal processing has been established as a relatively complete theoretical system through the past 30 years.

All of the methods mentioned above fall into one-dimensional vector or two-dimensional matrix modeling. However, received signals of the EMVS array have a multidimensional structure naturally, and the structural information associated with parameters exists among each dimension [12,13,14,15]. The traditional methods stretch the multidimensional data into matrices or vectors, which may lead to the loss of the structural information. In recent years, tensor-decomposition-based methods are more and more applied in the signal processing field [16,17,18]. Compared with matrix and vector modeling, the tensor modeling is more suitable to the real signal with a multidimensional structure. Besides, tensor decomposition possesses the blind estimation feature, which can directly estimate the source signal from measurements without the need for any prior information. With the help of such characteristics, even if the noise is non-Gaussian, tensor-based methods can still be applied to achieve parameter estimation. Tensor decomposition was proposed in psychometrics first in the 1970s [19], which also has been obtained by in-depth investigation in phonetics [20] and chemometrics [21]. In the 1990s, the higher order statistics (HOS) attracted attention in the signal processing field [22]. It is found that HOS is essentially interconnected with tensor modeling. In fact, the HOS-based independent component analysis (ICA) [23,24] implied the idea of tensor decomposition. Combining the tensor decomposition model with the array manifold, the canonical polyadic decomposition (CPD) method has been introduced into the parameter estimation of 2D array signal processing [25,26]. In 2005, S. Miron et al. introduced the tensor modeling to the EMVS array and proposed a tensor-based multiple signal classification (MUSIC) algorithm [12]. However, this search algorithm only considered the invariances in the polarization dimension, i.e., only using the uni-mode orthogonality. In 2009, X.F. Gong et al. put forward a two-fold mode-projection MUSIC (TM-MUSIC) algorithm, applying it to the EMVS array [27], which can obtain a better performance than the long-vector MUSIC algorithm by fusing the orthogonality on both dimensions and taking full advantage of the structural information of the tensor data. In 2011, X.J. Guo et al. applied the CPD algorithm to the EMVS array and derived the uniqueness conditions of decomposition, as well as the identifiability of AOAs for coherent signals [28]. In 2014, K. Han and A. Nehorai constructed a nested vector sensor array and introduced the tensor modeling for AOA estimation [29]. In 2015, X.R. Zhang et al. presented an adaptive beamformer for the EMVS array based on tensor modeling with a two-stage technique [30], which has a good robustness even under the condition of steering vector mismatching. In 2016, P. Forster et al. analyzed the theoretical performance of the tensor MUSIC and classic long vector MUSIC in polarized source estimation and gave the corresponding theoretical mean square error (MSE) [31].

The existing tensor-decomposition-based parameter estimation methods for the EMVS array are mostly based on the CPD modeling. Such a model has the advantage of uniqueness, but decomposition factors must be rank-1, which may not be satisfied in practice [32]. Therefore, 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- 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- 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- 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- 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- 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- 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- 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- 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.

2. 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., . The tensor is indicated by calligraphic capital letter, e.g., . is the Frobenius norm of the tensor , 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 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 , is defined as the the unfolded matrix, or matricization of in mode-1, i.e., ; similarly, and denote the matricization in mode-2 and mode-3, respectively. The mode-n rank of tensor is the rank of its unfolded matrix .

2.1. EMVS Array Signal Model

Under the assumption of narrowband and far-field signal [7,8], the measurement of EMVS can be expressed as: where denote the azimuth and elevation angle of arrival signal, respectively, as shown in Figure 1.

Figure 1

The geometry of a single EMVS.

For a fully-polarized signal, , in which is the k-th source signal, and are the variables associated to polarization parameters. where is the polarization orientation angle of the electromagnetic wave and is the ellipticity angle, as shown in Figure 2. Thus, , where is referred to as the polarization steering vector corresponding to the k-th signal.

Figure 2

The electromagnetic polarization ellipse.

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 and considering noise, the received signal model of the EMVS array can be expressed as: where 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 , in which is the wavelength of source signal, we have . is the source signal, vector and 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.

Figure 3

L-shaped array configuration of the EMVS array.

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, cannot be incorporated into the polarization steering vector . For instance, the circularly-polarized signal has , where and 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 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 is the source signal matrix for the N snapshots and are the matrices of AWGN on the subarrays along the x- and z-axis, respectively.

2.2. 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-, rank- and rank- patterns. Different patterns have different forms of decomposition, corresponding to different signal models. This paper mainly pays attention to the rank- BCD, i.e., the Type-2 BCD model.

Type-2 BCD in rank- terms:

The decomposition of a tensor is called the Type-2 block component decomposition, if it can be expressed as the summation of a series rank- factors, where can be termed as the core tensor of the k-th component. The mode-1 rank and mode-2 rank of are and , respectively. Both and are of full column rank. The type-2 BCD model in rank- terms is schematized in

The BCD model can also be expressed as: where , and is defined as: and Equation (10) can be written as (cf. [16]).

3. The Proposed Method

3.1. 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: in which , , , and , where denotes a block diagonal matrix. Note that such a model is the rank- BCD [34], where . Then, Equation (13) can be rewritten as: where is the received signal; is the source signal, of which the mode-3 matricization is , i.e., . is the corresponding AWGN.

Next, we investigate if the rank condition of each factor is satisfied in Equation (14). First, the steering matrix is of full column rank when . The same is true for when . For , the size of mode-1 matricization is and for mode-2. Thus, the mode-1 and mode-2 rank of are and , 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 must be true. Additionally, since , no matter what the snapshots are, multiple or single, Equation (14) will be valid according to Equation (15).

3.1.1. Uniqueness Analysis

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

Essentially unique for rank- BCD:

The rank- BCD described in Equation ( and can be post-multiplied by a nonsingular matrix, i.e., , in which ; (2) corresponding to Condition 1, there exists the transformation of : .

If , and , the two sets of BCD expressions and can be considered as essentially unique.

Uniqueness of BCD modeling for the EMVS array:

The rank- BCD descried in Equation ( and Ψ are of full column rank, and all elements of 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 are surely characterized by jointly continuous probability density functions, i.e., the conditions of factors in Equation (14) are satisfied from Theorem 1. In addition, the array manifold based on the L-shaped array in this paper gives ; the structure of 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 are linearly independent. Hence, . Thus, when , is of full column rank. From the definition of in Section 2.1, we have . Thus, when , 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.

3.2. AOA Estimation Algorithm

Given the rank- BCD model (14) with additive noise, the task is to acquire the estimation of factor matrices , from the received signal and obtain the estimation of spatial frequency . This paper adopts the minimum mean square error (MMSE) criterion [35],

With the mode-n matricizations of the tensor and , i.e., and , , the signal model (14) can be rewritten as: in which , , .

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 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 and fixed, the update of is calculated as follows, then we can obtain: where represents the pseudo-inverse of a matrix. With and fixed, the same procedure gives the update of ,

Then, fixing and leads to the update of ,

During the steps of updating and , the normalization by using QR decomposition can be added to enhance the stability of estimation [35]. It is necessary to guarantee that , and are of full column rank in the iterations.

In addition, the structural information about factor matrices can be considered in the initialization step for accelerating the convergence of ALS. Let represent the n-th frontal slice matrix of the tensor , i.e., ; then given , we have: where is of full rank, known from the previous rank condition of rank- BCD and the characteristics of the source signal. is also full column rank in terms of the uniqueness analysis because Equation (16) is fulfilled. This means that the column space of is the same as that of . A similar conclusion can be deduced that 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 . After finding the steering matrix by means of ALS, the spacial frequency can be obtained. Note that the trivial indeterminacies of exists, i.e., where is a nonsingular matrix. Equation (26) indicates that the BCD estimation result and steering matrix have the same column space. Because 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 can be parted as , where , . In the same way, the estimation of is parted as . Then, combined with Equation (26), we have:

The above equations indicate that 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 , the left singular vector corresponding to the maximum singular value can be obtained. Take the first and last rows of , denoted as and , to construct . Applying the eigenvalue decomposition of , in which: yields the estimation of spacial frequency , 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 . After obtaining the estimation of and , the azimuth angle can be estimated directly,

Since 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 . Such a procedure consists of three parts. The first is to update the factor matrices, which has the complexity of . The second is to solve the linear systems, which has the complexity of . The third is function evaluation, which has the complexity of . In addition, we employ the SVD to calculate the spatial frequencies, which has the the complexity of . For the signal model of the EMVS array in this paper, we have , , , . It can be concluded that the total computational complexity is .

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

4.1. Implementations of 2D AOA Estimation

4.1.1. Different SNRs

The simulation parameters are configured as: , , ; the distance between adjacent sensors is ; two groups of the AOAs are and ; the SNR range is −3 dB∼15 dB; the Monte Carlo trials . The azimuth-elevation angles estimation is shown in Figure 5.

Figure 5

Azimuth-elevation scatter versus SNRs, for two AOAs at , , by 500 Monte Carlo trials (). (a) SNR = –3 dB, (b) SNR = 3 dB, (c) SNR = 9 dB, (d) SNR = 15 dB.

For comparing the accuracy of estimation under different SNRs more directly, we analyze the simulation result for one group of AOA, , 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 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.

Figure 6

The AOA estimation scatter for versus SNRs, by 500 Monte Carlo trials ().

4.1.2. 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 , , . The Monte Carlo simulation trials .

Figure 7

Azimuth-elevation scatter versus angular separations, for two AOAs at , , by 500 Monte Carlo trials (). (a) Angular Separatio = , (b) Angular Separatio = , (c) Angular Separatio = , (d) Angular Separatio = , (e) Angular Separatio = , (f) Angular Separatio = .

From the scatter diagram, the AOA estimation can be fully resolved when the angular separation is large than . Even if the separation degree is very small, 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 , the detection probability of the proposed algorithm is larger than 80%. However, the ESPRIT has only about 40%.

Figure 8

The AOA estimation scatter for versus angular separations, by 500 Monte Carlo trials ().

4.1.3. Different Snapshots

This group of experiments investigates the estimation of AOA under different snapshots, which is set as . 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 , the simulation result is shown in Figure 9.

Figure 9

Azimuth-elevation scatter versus snapshots, for two AOAs at , , by 500 Monte Carlo trials (). (a) snapshots = 5, (b) snapshots = 50, (c) snapshots = 200, (d) snapshots = 500.

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 rises with the snapshot increasing.

Figure 10

The AOA estimation scatter for versus snapshots, by 500 Monte Carlo trials ().

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

4.2.1. RMSE

First, we investigate the RMSE of each method, which is defined as: where 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 and are the diagonal elements of CRB matrix [60,61], corresponding to the parameters and , respectively.

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

Figure 11

RMSE versus SNRs, by 500 Monte Carlo trials ().

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 , 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 for all algorithms with the angular separation larger than . Below this value, the detection probability is diverse for different algorithms. It is noticed that when the angular separation is , 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.

Figure 12

RMSE versus angular separations, by 500 Monte Carlo trials ().

Figure 15

Detection probability versus angular separations, by 500 Monte Carlo trials ().

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 , the uniqueness of BCD can be met as long as snapshots . 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.

Figure 13

RMSE versus snapshots, by 500 Monte Carlo trials ().

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

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

Figure 14

Detection probability versus SNRs, by 500 Monte Carlo trials ().

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 (). For instance, when , 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 , the detection probabilities of all methods have significant improvement. Meanwhile, the proposed algorithm is still optimal. When , all of the methods achieve 100% detection probability.

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.

Figure 16

Detection probability versus snapshots, by 500 Monte Carlo trials ().

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