PARAFAC Estimators for Coherent Targets in EMVS-MIMO Radar with Arbitrary Geometry

Abstract

In the past few years, multiple-input multiple-output (MIMO) radar with electromagnetic vector sensor (EMVS) array, or called EMVS-MIMO radar, has attracted extensive attention in target detection. Unlike the traditional scalar sensor-based MIMO radar, an EMVS-MIMO radar can not only provides a two-dimensional (2D) direction finding of the targets but also offers 2D polarization parameter estimation, which may be important for detecting weak targets. In this paper, we investigate into multiple parameter estimations for a bistatic EMVS-MIMO radar in the presence of coherent targets, whose transmitting EMVS and receiving EMVS are placed in an arbitrary topology. Three tensor-aware spatial smoothing estimators are introduced. The core of the proposed estimators is to de-correlate the coherent targets via the spatial smoothing technique and then formulate the covariance matrix into a third-order parallel factor (PARAFAC) tensor. After the PARAFAC decomposition of the tensor, the factor matrices can be obtained. Thereafter, the 2D direction finding can be accomplished via the normalized vector cross-product technique. Finally, the 2D polarization parameter can be estimated via the least squares method. Unlike the state-of-the-art PARAFAC estimator, the proposed estimators are suitable for arbitrary sensor geometries, and they are robust to coherent targets as well as sensor position errors. In addition, they have better estimation performance than the current matrix-based estimators. Moreover, they are computationally efficient than the current subspace methods, especially in the presence of a large-scale sensor array. In addition, the proposed estimators are analyzed in detail. Numerical experiments coincide with our theoretical findings.

1. Introduction

Multiple-input multiple-output (MIMO) radar has attracted widespread attention in remote sensing. MIMO radar is characterized by multiple transmitting (Tx) sensors and multiple receiving (Rx) sensors [1]. A MIMO radar emits mutual orthogonal waveforms and receives the echoes simultaneously [2]. The waveform diversity enables a MIMO radar to achieve better performance than a phased-array radar in noise suppression, interference elimination, and parameter identifiability [3]. Among various kinds of MIMO radars, the bistatic MIMO radar is attractive as the Tx array and the Rx array are separated [4]. Benefitting from such a configuration, a bistatic MIMO radar would have better survival ability than a monostatic one in military applications.

The issue of direction-of-departure (DOD) and direction-of-arrival (DOA) estimation is critical to a bistatic MIMO radar. Many efforts have been devoted to this topic in the past decades. A lot of super-resolution algorithms have been reported, for instance, estimation approach to signal parameters with rotational invariance technique (ESPRIT) [5,6], spatial spectrum peak search [7,8,9], and tensor methods [10,11,12,13]. Generally speaking, tensor-aware methods provide much better performance than the matrix-based one [14,15,16,17,18,19], since the multidimensional structure of the array measurement can be taken into consideration by the latter. However, most of the current algorithms mainly focus on the one-dimensional (1D) direction estimation issue by using the 1D scalar linear array, e.g., uniform linear array (ULA), and non-uniform linear array (nested array, coprime array). In engineering applications, two-dimensional (2D) direction estimation, i.e., 2D-DOD and 2D-DOA estimation, is more attractive than the former. It is well-known that in order to obtain 2D angle estimation using the scalar-array-based MIMO radar, non-linear Tx arrays and Rx arrays must be adopted [20,21,22,23], e.g., rectangular array, circular array, arbitrary array. Unfortunately, these nonlinear scalar geometries are usually too complex to be conformal with the platform. In addition, the calibration of the model errors will become more complicated due to the irregular geometries [24,25].

On the other hand, the electromagnetic vector sensor (EMVS) array has brought new insight into direction estimation [26,27,28,29]. A complete EMVS occupies six components: three of them measure the electric-field and three of them sense the magnetic-field. Unlike the traditional scalar array, an EMVS array not only provides 2D direction estimation, but it also offers the polarization status of the incoming source. The EMVS array enables the MIMO radar to obtain 2D direction estimation as well as 2D polarization estimation of the targets. A bistatic EMVS-MIMO radar was introduced in [30], whose Tx array and Rx array are both ULA. Therein, an ESPRIT algorithm was introduced to estimate the 2D-DOD, the 2D-DOA, the 2D Tx polarization angle (2D-TPA) and the 2D Rx polarization angle (2D-RPA). Another ESPRIT algorithm was proposed in [31], which can avoid the aperture loss in [30]. A parallel factor (PARAFAC) estimator was designed in [32], in which the tensor nature of the array measurement was taken into account, and it offers much better estimation performance than those matrix approaches. In [33], the coprime array geometry-based Tx/Rx configuration was considered for EMVS-MIMO radar. Benefitting from the larger sensor aperture, it offers more accurate estimation performance than the ULA setup. In [34], an arbitrary Tx/Rx geometry was investigated for EMVS-MIMO radar. An ESPRIT-like estimator based on the normalized vector cross-product technique has been given, which is capable of providing closed-form solutions for 2D-DOD, 2D-DOA, 2D-TPA and 2D-RPA estimation. In addition, it is insensitive to the sensor position errors.

The above-mentioned algorithms are only suitable for uncorrelated targets scenario, which is often inappropriate in practice. In the presence of coherent targets, there will be a rank deficient in the array measurement. To alleviate this issue, the polarization smoothing was performed in [35], which resolves the coherent targets via smoothing the array measurement along the polarization domain. A polarization difference smoothing framework was presented in [36], which is suitable for nonstationary noise. Nevertheless, both [35,36] would sacrifice the polarization information of the targets. To avoid such drawback, the generalized spatial smoothing framework was carried out in [37]. Three spatial smoothing approaches, named Tx smoothing (TS), Rx smoothing (RS), and Tx/Rx smoothing (TRS), were derived, which are suitable for arbitrary Tx/Rx geometry. However, the eigendecomposition of TS, RS and TRS are computationally inefficient. In addition, since the tensor nature has been ignored, the estimation accuracy of these approaches can be improved further.

In this paper, tensor-aware approaches are proposed to tackle the coherent targets problem in bistatic EMVS-MIMO radar. More specifically, the contributions of this paper are listed as follows:

(1)

The generalized spatial smoothing-based tensor models are established. The core of the proposed algorithm is to solve the rank-deficiency issue via spatial smoothing. To utilize the multi-dimensional structure of the array measurement, the array data after smoothing is rearranged into a third-order PARAFAC tensor. Unlike the state-of-the-art PARAFAC estimator in [32], the improved PARAFAC approaches are robust to coherent targets owing to the spatial smoothing. In addition, since the third-order PARAFAC decomposition can be quickly accomplished via the existing COMFAC algorithm, the proposed frameworks are more effective than the existing fourth-order PARAFAC algorithms [12].

(2)

The ESPRIT-like strategies are developed for joint 2D-DOD, 2D-DOA, 2D-TPA and 2D-RPA estimation. After PARAFAC decomposition, the factor matrices that form the tensor are achieved. The 2D-DOD and 2D-DOA are estimated via the normalized vector cross-product technique. Thereafter, 2D-TPA and 2D-RPA are achieved via the least squares (LS) with the previously estimated 2D-DOD and 2D-DOA. All the estimated parameters are in closed-form and paired automatically. Since the multi-dimensional structure has been explored, they outperform the matrix-based smoothing methods in [37]. Furthermore, as the 2D-DOD and 2D-DOA estimation rely on the normalized vector cross-product technique instead of the uniformity of the sensor array, the proposed approaches are suitable for arbitrary geometries and sensor position errors, while the PARAFAC estimator in [32] is only effective for the ULA configuration.

(3)

The advantages of the proposed approaches are verified via theoretically analysis and simulations. The proposed approaches are analyzed with respect to computational complexity, identifiability, as well as the Cramer–Rao Bound (CRB). Computer simulations are designed to show its effectiveness.

2. Tensor and Problem Formulation

2.1. Tensor and PARAFAC Decomposition

The following preliminaries concerning tensor and tensor decomposition will be utilized throughout this paper. Readers that are interested in more details on the tensor are referred to [38].

Definition 1.

An $I_1\times I_2\times \cdots \times I_N$ tensor has N degrees-of-freedom or indices. The mode-n unfolding of $\mathcal{X}\in {ℂ}^{I_1\times I_2\times \cdots \times I_N}$ is given by ${\left[\mathcal{X}\right]}_{\left(n\right)}$, which unfolds a tensor into a matrix, where the $(i_1,i_{2,}\cdots ,i_N)$ -entity of $\mathcal{X}$ is the $(i_n,p)$ -th entity of ${\left[\mathcal{X}\right]}_{\left(n\right)}$, where $p=1+{\displaystyle {\sum }_{k=1,k\ne n}^N(i_k-1)J_k}$ with $J_k={\displaystyle {\mathsf{\prod }}_{m=1,m\ne n}^{k-1}{I}_m}$.

Definition 2.

The mode-n product of matrix $\mathit{A}\in {ℂ}^{J_n\times I_n}$ and $\mathcal{X}\in {ℂ}^{I_1\times I_2\times \cdots \times I_N}$ is given by $\mathcal{Y}={\mathcal{X}}_{\times n}\mathit{A}$, which satisfies:

$${\left[\mathcal{Y}\right]}_{\left(n\right)}=\mathit{A}{\left[\mathcal{X}\right]}_{\left(n\right)}$$

(1)

Definition 3.

The PARAFAC decomposition of a rank-K tensor $\mathcal{X}\in {ℂ}^{I_1\times I_2\times \cdots \times I_N}$ is given by

$$\begin{array}{ll}\mathcal{X}& ={\displaystyle \sum _{k=1}^K{\mathit{a}}_k^{(1)}\circ {\mathit{a}}_k^{(2)}\cdots {\mathit{a}}_k^{(N)}}\\ & ={\mathcal{I}}_{\times 1}{\mathit{A}}_{\times 2}^{(1)}{\mathit{A}}_{\times 3}^{(2)}\cdots {\mathit{A}}_{\times N}^{(N-1)}{\mathit{A}}^{(N)}\end{array}$$

(2)

where $\circ$ denotes the outer-product, $\mathcal{I}\in {ℂ}^{I_1\times I_2\times \cdots \times I_N}$ is the identical tensor, ${\mathit{a}}_k^{(n)}\in {ℂ}^{I_n\times 1}(n=1,2,\cdots ,N)$ is a rank-one vector, and ${\mathit{A}}^{(n)}=\left[{\mathit{a}}_1^{(n)},{\mathit{a}}_1^{(n)},\cdots ,{\mathit{a}}_K^{(n)}\right]$ is the so-called factor matrix. In matrix format, (2) is also given by

$${\left[\mathcal{X}\right]}_{\left(n\right)}={\mathit{A}}_n{\left[{\mathit{A}}_{n+1}\odot \cdots \odot {\mathit{A}}_N\odot {\mathit{A}}_1\cdots \odot {\mathit{A}}_{n-1}\right]}^T$$

(3)

where $\odot$ denotes the Khatri–Rao product, and ${(\cdot )}^T$ denotes transpose.

Definition 4.

For the PARAFAC decomposition in (2), define the order set ${\mathbb{O}}_j=\left\{o_{j,1},o_{j,2},\cdots ,o_{j,M_j}\right\}$, $j=1,2,\cdots ,J$, which denotes a partitioning of the dimensions $\mathbb{O}=\left\{1,2,\cdots ,N\right\}$, a new PARAFAC tensor ${\mathcal{X}}_{{\mathbb{O}}_1,{\mathbb{O}}_2,\cdots ,{\mathbb{O}}_J}\in {ℂ}^{T_1\times T_2\times \cdots \times T_J}$ can be obtained via combining the rank-one vectors according to the order set as [39]:

$${\mathcal{X}}_{{\mathbb{O}}_1,{\mathbb{O}}_2,\cdots ,{\mathbb{O}}_J}={\displaystyle \sum _{k=1}^K{\mathit{b}}_k^{(1)}}\circ {\mathit{b}}_k^{(2)}\circ \cdots \circ {\mathit{b}}_k^{(J)}$$

(4)

where $T_j={\mathsf{\Pi }}_{m=1}^{M_j}{I}_{o_{j,m}}$, ${\mathit{b}}_k^{(j)}={\mathit{a}}_k^{\left(o_{j,M_j}\right)}\otimes {\mathit{a}}_k^{\left(o_{j,M_j-1}\right)}\cdots \otimes {\mathit{a}}_k^{\left(o_{j,1}\right)}$, and $\otimes$ denotes the Kronecker product.

2.2. Signal Model

Consider a bistatic EMVS-MIMO radar with arbitrary sensor geometry, as shown in Figure 1. Without loss of generality, we assume that there are M-element colocated transmitting EMVSs and N-element colocated receiving EMVSs. The positions with respect to the m-th ($m=1,2,\cdots ,M$) Tx EMVS and the n-th ($n=1,2,\cdots ,N$) EMVS are denoted by ${\mathit{p}}_{t,m}\triangleq {\left[x_{t,m},y_{t,m},z_{t,m}\right]}^T$ and ${\mathit{p}}_{r,n}\triangleq {\left[x_{r,n},y_{r,n},z_{r,n}\right]}^T$, respectively. Supposing that $K$ far-field targets appearing in same range bin, the 2D-DOD, 2D-DOA, 2D-TPA, and 2D-RPA of the k-th ($k=1,2,\cdots ,K$) target are denoted by $\left({\theta }_{r,k},{\varphi }_{r,k}\right)$, $\left({\theta }_{t,k},{\varphi }_{t,k}\right)$, $\left({\gamma }_{r,k},{\eta }_{r,k}\right)$ and $\left({\gamma }_{t,k},{\eta }_{t,k}\right)$, where ${\theta }_{t,k}$ and ${\theta }_{r,k}$ are the Tx elevation angle and the Rx elevation angle, respectively, ${\varphi }_{t,k}$ and ${\varphi }_{r,k}$ are the Tx azimuth angle and the Rx azimuth angle, respectively, ${\gamma }_{t,k}$ and ${\gamma }_{r,k}$ are the Tx auxiliary polarization angle and the Rx auxiliary polarization angle, respectively, and ${\eta }_{t,k}$ and ${\eta }_{r,k}$ are the Tx polarization phase difference and the Rx polarization phase difference, respectively.

Suppose that the Tx array emit $6M$ mutual orthogonal pulse waveforms $\mathit{w}\left(t,\tau \right)={\left[w_1\left(t,\tau \right),w_2\left(t,\tau \right),\cdots ,w_{6M}\left(t,\tau \right)\right]}^T$, i.e.,

$${\displaystyle \underset{T}{\int }{w}_m\left(t,\tau \right)w_n^{*}\left(t,\tau \right)d\tau =\delta \left(m-n\right)}$$

(5)

where $t$ denotes the pulse index, $\tau$ denotes the fast time index, ${(\cdot )}^{*}$ denotes conjugate, $T$ denotes the pulse duration, and $\delta (\cdot )$ denotes the Kronecker delta. In the presence of orthogonal waveforms assumption, the array outputs from the matched filters are given by [31,32]:

$$\begin{array}{ll}\mathit{x}\left(t\right)& ={\displaystyle \sum _{k=1}^K\left[{\mathit{a}}_{t,k}\otimes {\mathit{b}}_{t,k}\otimes {\mathit{a}}_{r,k}\otimes {\mathit{b}}_{r,k}\right]s_k\left(t\right)+\mathit{n}\left(t\right)}\\ & =\left[{\mathit{A}}_t\odot {\mathit{B}}_t\odot {\mathit{A}}_r\odot {\mathit{B}}_r\right]\mathit{s}\left(t\right)+\mathit{n}\left(t\right)\end{array}$$

(6)

where ${\mathit{A}}_t=\left[{\mathit{a}}_{t,1},{\mathit{a}}_{t,2},\cdots ,{\mathit{a}}_{t,K}\right]$, ${\mathit{A}}_r=\left[{\mathit{a}}_{r,1},{\mathit{a}}_{r,2},\cdots ,{\mathit{a}}_{r,K}\right]$; ${\mathit{B}}_t=\left[{\mathit{b}}_{t,1},{\mathit{b}}_{t,2},\cdots ,{\mathit{b}}_{t,K}\right]$ and ${\mathit{B}}_r=\left[{\mathit{b}}_{r,1},{\mathit{b}}_{r,2},\cdots ,{\mathit{b}}_{r,K}\right]$; $\mathit{s}\left(t\right)={\left[s_1\left(t\right),s_2\left(t\right),\cdots ,s_K\left(t\right)\right]}^T$ denotes the target reflection coefficient vector; $\mathit{n}\left(t\right)$ denotes the additive white Gaussian noise (AWGN) vector with variance is ${\sigma }^2$, namely, $E\left\{\mathit{n}\left(t\right){\mathit{n}}^H\left(t\right)\right\}={\sigma }^2{\mathit{I}}_{36MN}$, where ${(\cdot )}^H$ denotes conjugate transpose, ${\mathit{I}}_{36MN}$ denotes the $36MN\times 36MN$ identity matrix, and $E\{\cdot \}$ is the mathematical expectation; ${\mathit{a}}_{t,k}={\left[1,e^{-j2\pi {\tau }_{t2,k}/\lambda },\cdots ,e^{-j2\pi {\tau }_{tM,k}/\lambda }\right]}^T\in {ℂ}^{M\times 1}$ denotes the k-th Tx spatial steering vector, ${\tau }_{tm,k}={\mathit{p}}_{t,m}^T{\mathit{r}}_{t,k}$ with ${\mathit{r}}_{t,k}={\left[\sin {\theta }_{t,k}\cos {\varphi }_{t,k},\sin {\theta }_{t,k}\sin {\varphi }_{t,k},\cos {\theta }_{t,k}\right]}^T$; ${\mathit{a}}_{r,k}={\left[1,e^{-j2\pi {\tau }_{r2,k}/\lambda },\cdots ,e^{-j2\pi {\tau }_{rN,k}/\lambda }\right]}^T\in {ℂ}^{N\times 1}$ denotes the k-th Rx spatial steering vector, ${\tau }_{rn,k}={\mathit{p}}_{r,n}^T{\mathit{r}}_{r,k}$ with ${\mathit{r}}_{r,k}={\left[\sin {\theta }_{r,k}\cos {\varphi }_{r,k},\sin {\theta }_{r,k}\sin {\varphi }_{r,k},\cos {\theta }_{r,k}\right]}^T$; and ${\mathit{b}}_{t/r,k}\in {ℂ}^{6\times K}$ denote the k-th Tx/Rx polarization response vector, which is given by

$${\mathit{b}}_{t/r,k}=\left[\begin{array}{c}\cos {\varphi }_{t/r,k}\cos {\theta }_{t/r,k}\sin {\gamma }_{t/r,k}{e}^{j{\eta }_{t/r,k}}-\sin {\varphi }_{t/r,k}\cos {\gamma }_{t/r,k}\\ \sin {\varphi }_{t/r,k}\cos {\theta }_{t/r,k}\sin {\gamma }_{t/r,k}{e}^{j{\eta }_{t/r,k}}+\cos {\varphi }_{t/r,k}\cos {\gamma }_{t/r,k}\\ -\sin {\theta }_{t/r,k}\sin {\gamma }_{t/r,k}{e}^{j{\eta }_{t/r,k}}\\ -\sin {\varphi }_{t/r,k}\sin {\gamma }_{t/r,k}{e}^{j{\eta }_{t/r,k}}-\cos {\varphi }_{t/r,k}\cos {\theta }_{t/r,k}\cos {\gamma }_{t/r,k}\\ \cos {\varphi }_{t/r,k}\sin {\gamma }_{t/r,k}{e}^{j{\eta }_{t/r,k}}-\sin {\varphi }_{t/r,k}\cos {\theta }_{t/r,k}\cos {\gamma }_{t/r,k}\\ \sin {\theta }_{t/r,k}\cos {\gamma }_{t/r,k}\end{array}\right]\begin{array}{c}\begin{array}{l}\\ \\ \end{array}\}{\mathit{e}}_{t/r,k}\in {ℂ}^{3\times 1}\\ \begin{array}{l}\\ \\ \end{array}\}{\mathit{m}}_{t/r,k}\in {ℂ}^{3\times 1}\end{array}$$

(7)

where ${\mathit{e}}_{t/r,k}$ denotes the electric-field response vector, and ${\mathit{m}}_{t/r,k}$ denotes the magnetic-field components. Moreover, the vector cross-products with respect to ${\mathit{e}}_{t,k}$, ${\mathit{m}}_{t,k}$ and ${\mathit{e}}_{r,k}$, ${\mathit{m}}_{r,k}$ are given by the normalized Poynting-vector as:

$$\begin{array}{ccc}\frac{{\mathit{e}}_{t,k}}{{\Vert {\mathit{e}}_{t,k}\Vert }_F}⊛\frac{{\mathit{m}}_{t,k}^{*}}{{\Vert {\mathit{m}}_{t,k}\Vert }_F}={\mathit{r}}_{t,k}& ,& \frac{{\mathit{e}}_{r,k}}{{\Vert {\mathit{e}}_{r,k}\Vert }_F}⊛\frac{{\mathit{m}}_{r,k}^{*}}{{\Vert {\mathit{m}}_{r,k}\Vert }_F}={\mathit{r}}_{r,k}\end{array}$$

(8)

where $⊛$ denotes the vector cross-product, and ${\Vert \cdot \Vert }_F$ denotes the Frobenius norm. In addition, the ${\mathit{b}}_{t,k}$ and ${\mathit{b}}_{r,k}$ can be factorized as:

$$\begin{array}{ccc}{\mathit{b}}_{t,k}={\mathit{F}}_{t,k}{\mathit{v}}_{t,k}& ,& {\mathit{b}}_{r,k}={\mathit{F}}_{r,k}{\mathit{v}}_{r,k}\end{array}$$

(9)

where ${\mathit{F}}_{t,k}$ and ${\mathit{F}}_{r,k}$ are direction-dependent-only matrices, and ${\mathit{v}}_{t,k}$ and ${\mathit{v}}_{r,k}$ are polarization dependent-only vectors, which are given below:

$${\mathit{F}}_{t,k}=\left[\begin{array}{cc}\cos {\varphi }_{t,k}\cos {\theta }_{t,k}& -\sin {\varphi }_{t,k}\\ \sin {\varphi }_{t,k}\cos {\theta }_{t,k}& \cos {\varphi }_{t,k}\\ -\sin {\theta }_{t,k}& 0\\ -\sin {\varphi }_{t,k}& -\cos {\varphi }_{t,k}\cos {\theta }_{t,k}\\ \cos {\varphi }_{t,k}& -\sin {\varphi }_{t,k}\cos {\theta }_{t,k}\\ 0& \sin {\theta }_{t,k}\end{array}\right],{\mathit{v}}_{t,k}=\left[\begin{array}{l}\sin {\gamma }_{t,k}{e}^{j{\eta }_{t,k}}\\ \cos {\gamma }_{t,k}\end{array}\right]$$

(10)

Define ${\mathit{C}}_t={\mathit{A}}_t\odot {\mathit{B}}_t$, ${\mathit{C}}_r={\mathit{A}}_r\odot {\mathit{B}}_r$. In addition, suppose that the noise vector $\mathit{n}\left(t\right)$ is uncorrelated with the target reflection coefficient vector $\mathit{s}\left(t\right)$; then, the covariance matrix of $\mathit{x}\left(t\right)$ is given by:

$$\begin{array}{ll}\mathit{R}& =E\left\{\mathit{x}\left(t\right){\mathit{x}}^H\left(t\right)\right\}\\ & =\left[{\mathit{C}}_t\odot {\mathit{C}}_r\right]{\mathit{R}}_s{\left[{\mathit{C}}_t\odot {\mathit{C}}_r\right]}^H+{\sigma }^2{\mathit{I}}_{36MN}\end{array}$$

(11)

where ${\mathit{R}}_s=E\left\{\mathit{s}\left(t\right){\mathit{s}}^H\left(t\right)\right\}$ is the mathematical expectation of $\mathit{s}\left(t\right)$. In the presence of correlated targets, ${\mathit{R}}_s$ is rank deficient, i.e., $rank\left\{{\mathit{R}}_s\right\}<K$, $rank\{\cdot \}$ returns the rank of a matrix. Consequently, the noise $\mathit{R}$ is also rank deficient and yields ill signal subspace. As a result, the subspace methods in [30,31,33,34] will fail to work, and the PARAFAC method in [33] cannot correctly operate either.

3. The Proposed PARAFAC Estimators

3.1. Review of the Generalized Spatial Smoothing Approaches

The core idea of the generalized spatial smoothing approaches is to smooth the array measurement along the spatial domain [37]. According to the mechanism of spatial smoothing, the generalized spatial smoothing approaches can be divided into three categories, namely, the transmitting smoothing-based (TS) approach, the receiving smoothing-based (RS) approach, and the transmitting–receiving-smoothing-based (TRS) approach. The TS approach relies on smoothing the array measurement along the Tx direction. Let ${\mathit{x}}_{t,m}\left(t\right)$ denotes the array measurement corresponding to the m-th Tx EMVS, i.e.,

$${\mathit{x}}_{t,m}\left(t\right)=\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]{\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}\mathit{s}\left(t\right)+{\mathit{n}}_{t,m}\left(t\right)$$

(12)

where ${\mathit{n}}_{t,m}\left(t\right)$ is the corresponding noise vector, and ${\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}$ returns a diagonal matrix whose diagonal elements are the m-th row of ${\mathit{A}}_t$. Accordingly, the covariance matrix of ${\mathit{x}}_{t,m}\left(t\right)$ is given by

$$\begin{array}{ll}{\mathit{R}}_{t,m}& =E\left\{{\mathit{x}}_{t,m}\left(t\right){\mathit{x}}_{t,m}^H\left(t\right)\right\}\\ & =\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]{\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}{\mathit{R}}_s{\mathcal{S}}_m^H\left\{{\mathit{A}}_t\right\}{\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]}^H+{\sigma }^2{\mathit{I}}_{36N}\end{array}$$

(13)

Thereafter, a TS matrix ${\mathit{R}}_t$ can be constructed via averaging all the ${\mathit{R}}_{t,m}$, i.e.,

$$\begin{array}{ll}{\mathit{R}}_t& =1/M{\displaystyle {\sum }_{m=1}^M{\mathit{R}}_{t,m}}\\ & =\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]{\mathit{R}}_{s,t}{\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]}^H+{\sigma }^2{\mathit{I}}_{36N}\end{array}$$

(14)

where ${\mathit{R}}_{s,t}=1/M{\displaystyle {\sum }_{m=1}^M{\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}{\mathit{R}}_s{\mathcal{S}}_m^H\left\{{\mathit{A}}_t\right\}}$. According to [40], ${\mathit{R}}_{s,t}$ is full rank if $M\ge K$ and $\left({\theta }_{t,k},{\varphi }_{t,k}\right)$ are distinct. Accordingly, the rank of the noiseless ${\mathit{R}}_t$ is K, and then the subspace-based ESPRIT is capable for parameter estimation.

Likewise, we can pick up the array measurement associated with the n-th Rx EMVS, which is given by

$${\mathit{x}}_{r,n}\left(t\right)=\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]{\mathcal{S}}_n\left\{{\mathit{A}}_r\right\}\mathit{s}\left(t\right)+{\mathit{n}}_{r,n}\left(t\right)$$

(15)

where ${\mathit{n}}_{r,n}\left(t\right)$ accounts for the associated noise vector. Thereafter, we can construct an RS matrix ${\mathit{R}}_r$ via averaging all the covariance matrices of ${\mathit{x}}_{r,n}\left(t\right)$ (denoted by ${\mathit{R}}_{r,n}$), i.e.,

$$\begin{array}{ll}{\mathit{R}}_r& =1/N{\displaystyle {\sum }_{n=1}^N{\mathit{R}}_{r,n}}\\ & =\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]{\mathit{R}}_{s,r}{\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]}^H+{\sigma }^2{\mathit{I}}_{36M}\end{array}$$

(16)

where ${\mathit{R}}_{s,r}=1/N{\displaystyle {\sum }_{n=1}^N{\mathcal{S}}_n\left\{{\mathit{A}}_r\right\}{\mathit{R}}_s{\mathcal{S}}_n^H\left\{{\mathit{A}}_r\right\}}$. Obviously, the rank of the noiseless ${\mathit{R}}_r$ is K if $N\ge K$ and $\left({\theta }_{r,k},{\varphi }_{r,k}\right)$ are distinct.

In addition, we can choose the array measurement associate with the m-th Tx EMVS and the n-th Rx EMVS, which is given by

$${\mathit{x}}_{tr,m,n}\left(t\right)=\left[{\mathit{B}}_t\odot {\mathit{B}}_r\right]{\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}{\mathcal{S}}_n\left\{{\mathit{A}}_r\right\}\mathit{s}\left(t\right)+{\mathit{n}}_{tr,m,n}\left(t\right)$$

(17)

where ${\mathit{n}}_{tr,m,n}\left(t\right)$ represents the corresponding noise vector. Let the covariance matrix of ${\mathit{x}}_{tr,m,n}\left(t\right)$ is ${\mathit{R}}_{tr,m,n}$. Similarly, the TRS approach is based on the averaging all the ${\mathit{R}}_{tr,m,n}$, i.e.,

$$\begin{array}{ll}{\mathit{R}}_{tr}& =1/MN{\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{m=1}^M{\mathit{R}}_{tr,m,n}}}\\ & =\left[{\mathit{B}}_t\odot {\mathit{B}}_r\right]{\mathit{R}}_{s,tr}{\left[{\mathit{B}}_t\odot {\mathit{B}}_r\right]}^H+{\sigma }^2{\mathit{I}}_{36}\end{array}$$

(18)

where ${\mathit{R}}_{s,tr}=1/NM{\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{m=1}^M{\mathcal{S}}_n\left\{{\mathit{A}}_r\right\}{\mathcal{S}}_m\left\{{\mathit{A}}_t\right\}{\mathit{R}}_s{\mathcal{S}}_n^H\left\{{\mathit{A}}_r\right\}{\mathcal{S}}_m^H\left\{{\mathit{A}}_t\right\}}}$. In addition, the rank of the noiseless ${\mathit{R}}_{tr}$ should be K once $MN\ge K$ and $\left({\theta }_{t,k},{\varphi }_{t,k}\right)$, $\left({\theta }_{r,k},{\varphi }_{r,k}\right)$ are distinct. In practice, ${\mathit{R}}_t$, ${\mathit{R}}_r$, and ${\mathit{R}}_{tr}$ can be estimated via $L$ snapshots as:

$$\begin{array}{c}{\hat{\mathit{R}}}_t=1/ML{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^M{\mathit{x}}_{t,m}\left(l\right){\mathit{x}}_{t,m}^H\left(l\right)}}\\ {\hat{\mathit{R}}}_r=1/NL{\displaystyle \sum _{l=1}^L{\displaystyle \sum _{n=1}^N{\mathit{x}}_{r,n}\left(l\right){\mathit{x}}_{r,n}^H\left(l\right)}}\\ {\hat{\mathit{R}}}_{tr}={\displaystyle \sum _{l=1}^L{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n=1}^N{\mathit{x}}_{tr,m,n}\left(l\right){\mathit{x}}_{tr,m,n}^H\left(l\right)}}}\end{array}$$

(19)

Essentially, the TS approach is only smoothing the data measurements along the Tx spatial direction. The RS approach, on the other hand, is only smoothing them along the Rx spatial direction. Finally, the TRS approach is smoothing the data in both Tx and Rx directions.

The existing subspace-based approaches rely on the eigen decomposition of the estimated covariance matrices, which are computationally inefficient. Since the tensor nature is ignored, the estimation accuracy can be improved further. In what follows, we will introduce the PARAFAC-based estimators, which can avoid the above drawbacks.

3.2. PARAFAC Models and PARAFAC Decomposition

For the covariance data model in (14), it can be formulated as a fourth-order PARAFAC model as:

$$\mathcal{R}={\mathcal{I}}_{\times 1}{\mathit{B}}_{t\times 2}{\mathit{C}}_{r\times 3}{\mathit{B}}_{t\times 4}^{*}\left({\mathit{C}}_r^{*}{\mathit{R}}_{s,t}\right)+\mathcal{N}$$

(20)

where $\mathcal{N}$ is the tensor version of ${\sigma }^2{\mathit{I}}_{36N}$. Actually, ${\mathit{R}}_{s,t}$ can be combined with any of the four factor matrices ${\mathit{B}}_t$, ${\mathit{C}}_r$, ${\mathit{B}}_t^{*}$ and ${\mathit{C}}_r^{*}$. The estimation of the factor matrices can be accomplished via optimizing the following:

$$\begin{array}{ccc}{\mathit{B}}_t,{\mathit{C}}_r& =& \arg \min \end{array}{\Vert \widehat{ℛ}-{\mathcal{I}}_{\times 1}{\mathit{B}}_{t\times 2}{\mathit{C}}_{r\times 3}{\mathit{B}}_{t\times 4}^{*}\left({\mathit{C}}_r^{*}{\mathit{R}}_{s,t}\right)\Vert }_F^2$$

(21)

where $\hat{\mathcal{R}}$ denotes the arranged ${\hat{\mathit{R}}}_t$. Although the above issue can be solved via the quadrilinear alternative least squares (ALS) algorithm, it suffers from the slowness of the convergence speed. Inspired by the fact that there exists a fast algorithm, named COMFAC, for a third-order PARAFAC model, $\mathcal{R}$ can be arranged into a third-order PARAFAC model. For the PARAFAC model in (20), define ${\mathbb{O}}_1=\left\{1\right\}$, ${\mathbb{O}}_2=\left\{2\right\}$, and ${\mathbb{O}}_3=\left\{3,4\right\}$. According to Definition 4, we can obtain another PARAFAC model as:

$${\mathcal{R}}_t={\mathcal{I}}_{\times 1}{\mathit{B}}_{t\times 2}{\mathit{C}}_{r\times 3}{\mathit{D}}_t+{\mathcal{N}}_t$$

(22)

where ${\mathit{D}}_t={\mathit{B}}_t^{*}{\mathit{C}}_r^{*}{\mathit{R}}_{s,t}$, ${\mathcal{N}}_t$ denotes the arranged array noise. Similarly, we can obtain the rearranged array tensor ${\mathcal{R}}_r$ and ${\mathcal{R}}_{tr}$ as:

$${\mathcal{R}}_r={\mathcal{I}}_{\times 1}{\mathit{C}}_{t\times 2}{\mathit{B}}_{r\times 3}{\mathit{D}}_r+{\mathcal{N}}_r$$

(23)

$${\mathcal{R}}_{tr}={\mathcal{I}}_{\times 1}{\mathit{B}}_{t\times 2}{\mathit{B}}_{r\times 3}{\mathit{D}}_{tr}+{\mathcal{N}}_{tr}$$

(24)

where ${\mathit{D}}_r={\mathit{C}}_t^{*}{\mathit{B}}_r^{*}{\mathit{R}}_{s,r}$, ${\mathit{D}}_{tr}={\mathit{B}}_t^{*}{\mathit{B}}_r^{*}{\mathit{R}}_{s,tr}$, ${\mathcal{N}}_r$ and ${\mathcal{N}}_{tr}$ denotes the rearranged array noise.

Taking the model in (22) as an example, we will show how to obtain the factor matrices from the third-order PARAFAC tensor. Obviously, the factor matrices of ${\mathcal{R}}_t$ can be achieved via solving:

$$\begin{array}{ccc}{\mathit{B}}_t,{\mathit{C}}_r,{\mathit{D}}_t& =& \arg \min \end{array}{\Vert {\hat{\mathcal{R}}}_t-{\mathcal{I}}_{\times 1}{\mathit{B}}_{t\times 2}{\mathit{C}}_{r\times 3}{\mathit{D}}_t\Vert }_F^2$$

(25)

On the other hand, the PARAFAC tensor ${\mathcal{R}}_t$ in (22) can be unfolded into matrix format as

$$\begin{array}{l}{\mathit{X}}_1={\left[{\mathcal{R}}_t\right]}_{\left(1\right)}^T=\left[{\mathit{C}}_r\odot {\mathit{D}}_t\right]{\mathit{B}}_t^T+{\mathit{N}}_1\\ {\mathit{X}}_2={\left[{\mathcal{R}}_t^T\right]}_{\left(2\right)}^T=\left[{\mathit{D}}_t\odot {\mathit{B}}_t\right]{\mathit{C}}_r^T+{\mathit{N}}_2\\ {\mathit{X}}_3={\left[{\mathcal{R}}_t\right]}_{\left(3\right)}^T=\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]{\mathit{D}}_t^T+{\mathit{N}}_3\end{array}$$

(26)

where ${\mathit{N}}_1={\left[{\mathcal{N}}_t\right]}_{\left(1\right)}^T$, ${\mathit{N}}_2={\left[{\mathcal{N}}_t\right]}_{\left(2\right)}^T$ and ${\mathit{N}}_3={\left[{\mathcal{N}}_t\right]}_{\left(3\right)}^T$. As a result, the optimization in (25) can be rewritten as joint optimize:

$${\mathit{B}}_t,{\mathit{C}}_r,{\mathit{D}}_t=\{\begin{array}{l}\arg \min \Vert {\hat{\mathit{X}}}_1-[{\mathit{C}}_r\odot {\mathit{D}}_t]{\mathit{B}}_t^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{X}}}_2-\left[{\mathit{D}}_t\odot {\mathit{B}}_t\right]{\mathit{C}}_r^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{X}}}_3-\left[{\mathit{B}}_t\odot {\mathit{C}}_r\right]{\mathit{D}}_t^T{\Vert }_F^2\end{array}$$

(27)

The trilinear alternating least squares (TALS) is a popular solver of the joint optimization issue in (27). TLAS treats each sub-issue as a least squares (LS) problem and assumes two of the factor matrices are known. Then, it obtains the third factor matrix via LS successively. The LS fitting will repeat until algorithm convergence. Mathematically, the LS updates of ${\mathit{B}}_t$, ${\mathit{C}}_r$, and ${\mathit{D}}_t$ are given by:

$$\{\begin{array}{l}{\hat{\mathit{B}}}_t^T=[{\mathit{C}}_r\odot {\mathit{D}}_t{]}^{†}{\hat{\mathit{X}}}_1\\ {\hat{\mathit{C}}}_r^T=[{\mathit{D}}_r\odot {\mathit{B}}_t{]}^{†}{\hat{\mathit{X}}}_2\\ {\hat{\mathit{D}}}_t^T=[{\mathit{B}}_r\odot {\mathit{C}}_t{]}^{†}{\hat{\mathit{X}}}_3\end{array}$$

(28)

Usually, the convergence condition is set to ${\Vert {\hat{\mathit{X}}}_1-\left[{\mathit{C}}_r\odot {\mathit{D}}_t\right]{\mathit{B}}_t^T\Vert }_F^2\le \xi$ or the iteration steps are larger than a given integer; $\xi$ denotes a given threshold. Moreover, the COMFAC algorithm is often adopted [41], which only require a few iteration steps. Denote the Kruskal ranks of ${\mathit{B}}_t$, ${\mathit{C}}_r$, and ${\mathit{D}}_t$ as $KR\left\{{\mathit{B}}_t\right\}$, $KR\left\{{\mathit{C}}_r\right\}$, and $KR\left\{{\mathit{D}}_t\right\}$, respectively. It has been pointed out in [38] that if

$$KR\left\{{\mathit{B}}_t\right\}+KR\left\{{\mathit{C}}_r\right\}+KR\left\{{\mathit{D}}_t\right\}\ge 2K+2$$

(29)

then the estimates of ${\mathit{B}}_t$, ${\mathit{C}}_r$, and ${\mathit{D}}_t$ are unique to permutations and scaling effects of columns, i.e.,

$$\{\begin{array}{c}{\hat{\mathit{B}}}_t={\mathit{B}}_t{\mathsf{\Pi }}_1{\mathit{\Delta }}_1+{\mathit{E}}_1\\ {\hat{\mathit{C}}}_r={\mathit{C}}_r{\mathsf{\Pi }}_1{\mathit{\Delta }}_2+{\mathit{E}}_2\\ {\hat{\mathit{D}}}_t={\mathit{D}}_t{\mathsf{\Pi }}_1{\mathit{\Delta }}_3+{\mathit{E}}_3\end{array}$$

(30)

where ${\mathsf{\Pi }}_1$ denotes a permutation matrix, and ${\mathit{\Delta }}_1$, ${\mathit{\Delta }}_2$, and ${\mathit{\Delta }}_3$ denote the diagonal matrices with ${\mathit{\Delta }}_1{\mathit{\Delta }}_2{\mathit{\Delta }}_3={\mathit{I}}_K$. ${\mathit{E}}_1$, ${\mathit{E}}_2$, and ${\mathit{E}}_3$ account for the associated fitting errors.

For the RS approach, the matrix unfolding of ${\mathcal{R}}_r$ is given by

$$\begin{array}{l}{\mathit{Y}}_1={\left[{\mathcal{R}}_r\right]}_{\left(1\right)}^T=\left[{\mathit{B}}_r\odot {\mathit{D}}_r\right]{\mathit{C}}_t^T+{\mathit{N}}_1^{\prime }\\ {\mathit{Y}}_2={\left[{\mathcal{R}}_r\right]}_{\left(2\right)}^T=\left[{\mathit{D}}_r\odot {\mathit{C}}_t\right]{\mathit{B}}_r^T+{\mathit{N}}_2^{\prime }\\ {\mathit{Y}}_3={\left[{\mathcal{R}}_r\right]}_{\left(3\right)}^T=\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]{\mathit{D}}_r^T+{\mathit{N}}_3^{\prime }\end{array}$$

(31)

where ${\mathit{N}}_1^{\prime }={\left[{\mathcal{N}}_r\right]}_{\left(1\right)}^T$, ${\mathit{N}}_2^{\prime }={\left[{\mathcal{N}}_r\right]}_{\left(2\right)}^T$, and ${\mathit{N}}_3^{\prime }={\left[{\mathcal{N}}_r\right]}_{\left(3\right)}^T$ are the noise matrices. Accordingly, the factor matrices that formed ${\mathcal{R}}_r$ can be formulated as joint optimizers:

$${\mathit{B}}_r,{\mathit{D}}_r,{\mathit{C}}_t=\{\begin{array}{l}\arg \min \Vert {\hat{\mathit{Y}}}_1-[{\mathit{B}}_r\odot {\mathit{D}}_r]{\mathit{C}}_t^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{Y}}}_2-\left[{\mathit{D}}_r\odot {\mathit{C}}_t\right]{\mathit{B}}_r^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{Y}}}_3-\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]{\mathit{D}}_r^T{\Vert }_F^2\end{array}$$

(32)

Similarly, the LS updates of ${\mathit{C}}_t$, ${\mathit{B}}_r$, and ${\mathit{D}}_r$ are given by:

$$\{\begin{array}{c}{\hat{\mathit{C}}}_t^T={\left[{\mathit{B}}_r\odot {\mathit{D}}_r\right]}^{†}{\hat{\mathit{Y}}}_1\\ {\hat{\mathit{B}}}_r^T={\left[{\mathit{D}}_r\odot {\mathit{C}}_t\right]}^{†}{\hat{\mathit{Y}}}_2\\ {\hat{\mathit{D}}}_r^T={\left[{\mathit{C}}_t\odot {\mathit{B}}_r\right]}^{†}{\hat{\mathit{Y}}}_3\end{array}$$

(33)

The estimated factor matrices of ${\mathcal{R}}_r$ can be formulated as:

$$\{\begin{array}{c}{\hat{\mathit{C}}}_t={\mathit{C}}_t{\mathsf{\Pi }}_2{\mathit{\Delta }}_1^{\prime }+{\mathit{E}}_1^{\prime }\\ {\hat{\mathit{B}}}_r={\mathit{B}}_r{\mathsf{\Pi }}_2{\mathit{\Delta }}_2^{\prime }+{\mathit{E}}_2^{\prime }\\ {\hat{\mathit{D}}}_r={\mathit{D}}_r{\mathsf{\Pi }}_2{\mathit{\Delta }}_3^{\prime }+{\mathit{E}}_3^{\prime }\end{array}$$

(34)

where ${\mathit{\Delta }}_1^{\prime }$, ${\mathit{\Delta }}_2^{\prime }$, and ${\mathit{\Delta }}_3^{\prime }$ are diagonal matrices, and ${\mathit{E}}_1^{\prime }$, ${\mathit{E}}_2^{\prime }$, and ${\mathit{E}}_3^{\prime }$ are fitting errors.

For the TRS approach, the matrix unfolding of ${\mathcal{R}}_{tr}$ can be written as:

$$\begin{array}{l}{\mathit{Z}}_1={\left[{\mathcal{R}}_{tr}\right]}_{\left(1\right)}^T=\left[{\mathit{B}}_r\odot {\mathit{D}}_{tr}\right]{\mathit{B}}_t^T+{\mathit{N}}_1^{″}\\ {\mathit{Z}}_2={\left[{\mathcal{R}}_{tr}\right]}_{\left(2\right)}^T=\left[{\mathit{D}}_{tr}\odot {\mathit{B}}_t\right]{\mathit{B}}_r^T+{\mathit{N}}_2^{″}\\ {\mathit{Z}}_3={\left[{\mathcal{R}}_{tr}\right]}_{\left(3\right)}^T=\left[{\mathit{B}}_t\odot {\mathit{D}}_r\right]{\mathit{D}}_{tr}^T+{\mathit{N}}_3^{″}\end{array}$$

(35)

where ${\mathit{N}}_1^{″}={\left[{\mathcal{N}}_{tr}\right]}_{\left(1\right)}^T$, ${\mathit{N}}_2^{″}={\left[{\mathcal{N}}_{tr}\right]}_{\left(2\right)}^T$ and ${\mathit{N}}_3^{″}={\left[{\mathcal{N}}_{tr}\right]}_{\left(3\right)}^T$. Similarly, the factor matrices associated with ${\mathcal{R}}_{tr}$ can be achieved via optimizing the following minimizations simultaneously:

$${\mathit{B}}_t,{\mathit{B}}_r,{\mathit{D}}_{tr}=\{\begin{array}{l}\arg \min \Vert {\hat{\mathit{Z}}}_1-[{\mathit{B}}_r\odot {\mathit{D}}_{tr}]{\mathit{B}}_t^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{Z}}}_2-\left[{\mathit{D}}_{tr}\odot {\mathit{B}}_t\right]{\mathit{B}}_r^T{\Vert }_F^2\\ \arg \min \Vert {\hat{\mathit{Z}}}_3-\left[{\mathit{B}}_t\odot {\mathit{B}}_r\right]{\mathit{D}}_{tr}^T{\Vert }_F^2\end{array}$$

(36)

Therefore, the LS updates of ${\mathit{B}}_t$, ${\mathit{B}}_r$, and ${\mathit{D}}_{tr}$ are given by:

$$\{\begin{array}{c}{\left({\hat{\mathit{B}}}_t^{″}\right)}^T={\left[{\mathit{B}}_r\odot {\mathit{D}}_{tr}\right]}^{†}{\hat{\mathit{Z}}}_1\\ {\left({\hat{\mathit{B}}}_r^{″}\right)}^T={\left[{\mathit{D}}_{tr}\odot {\mathit{B}}_t\right]}^{†}{\hat{\mathit{Z}}}_2\\ {{\hat{\mathit{D}}}_{tr}^{″}}^T={\left[{\mathit{B}}_t\odot {\mathit{B}}_r\right]}^{†}{\hat{\mathit{Z}}}_3\end{array}$$

(37)

Consequently, the estimated factor matrices of and ${\mathcal{R}}_{tr}$ are given by:

$$\{\begin{array}{c}{\hat{\mathit{B}}}_t^{″}={\mathit{B}}_t{\mathsf{\Pi }}_3{\mathit{\Delta }}_1^{″}+{\mathit{E}}_1^{″}\\ {\hat{\mathit{B}}}_r^{″}={\mathit{B}}_r{\mathsf{\Pi }}_3{\mathit{\Delta }}_2^{″}+{\mathit{E}}_2^{″}\\ {\hat{\mathit{D}}}_{tr}={\mathit{D}}_{tr}{\mathsf{\Pi }}_3{\mathit{\Delta }}_3^{″}+{\mathit{E}}_3^{″}\end{array}$$

(38)

where ${\mathit{\Delta }}_1^{″}$, ${\mathit{\Delta }}_2^{″}$, and ${\mathit{\Delta }}_3^{″}$ are diagonal matrices, and ${\mathit{E}}_1^{″}$, ${\mathit{E}}_2^{″}$, and ${\mathit{E}}_3^{″}$ are error matrices.

3.3. 2D-DOD and 2D-DOA Estimation

It has been pointed in [34] that the polarization response vector ${\mathit{b}}_{t/r,k}$ can be represented in the normalized format as:

$${\mathit{b}}_{t/r,k}=b_{t/r,k}\left(1\right){\left[1,{\beta }_{t/r,k}^{\left(1,2\right)},\cdots ,{\beta }_{t/r,k}^{\left(1,6\right)}\right]}^T$$

(39)

where $b_{t/r,k}\left(p\right)$ denotes the p-th entity of ${\mathit{b}}_{t/r,k}$, ${\beta }_{t/r,k}^{\left(1,p\right)}=b_{t/r,k}\left(p\right)/b_{t/r,k}\left(1\right)$. It is easy to conclude that

$$\left(b_{t/r,k}\left(1\right)\left[\begin{array}{c}1\\ {\beta }_{t/r,k}^{\left(1,2\right)}\\ {\beta }_{t/r,k}^{\left(1,3\right)}\end{array}\right]\right)⊛{\left(b_{t/r,k}\left(1\right)\left[\begin{array}{c}{\beta }_{t/r,k}^{\left(1,4\right)}\\ {\beta }_{t/r,k}^{\left(1,5\right)}\\ {\beta }_{t/r,k}^{\left(1,6\right)}\end{array}\right]\right)}^{*}={\Vert b_{t/r,k}\left(1\right)\Vert }^2\left(\left[\begin{array}{c}1\\ {\beta }_{t/r,k}^{\left(1,2\right)}\\ {\beta }_{t/r,k}^{\left(1,3\right)}\end{array}\right]⊛{\left[\begin{array}{c}{\beta }_{t/r,k}^{\left(1,4\right)}\\ {\beta }_{t/r,k}^{\left(1,5\right)}\\ {\beta }_{t/r,k}^{\left(1,6\right)}\end{array}\right]}^{*}\right)$$

(40)

Define ${\tilde{\mathit{e}}}_{t,k}={\left[1,{\beta }_{t,k}^{\left(1,2\right)},{\beta }_{t,k}^{\left(1,3\right)}\right]}^T$, ${\tilde{\mathit{m}}}_{t,k}={\left[{\beta }_{t,k}^{\left(1,4\right)},{\beta }_{t,k}^{\left(1,5\right)},{\beta }_{t,k}^{\left(1,6\right)}\right]}^T$, ${\tilde{\mathit{e}}}_{r,k}={\left[1,{\beta }_{r,k}^{\left(1,2\right)},{\beta }_{r,k}^{\left(1,3\right)}\right]}^T$, and ${\tilde{\mathit{m}}}_{r,k}={\left[{\beta }_{r,k}^{\left(1,4\right)},{\beta }_{r,k}^{\left(1,5\right)},{\beta }_{r,k}^{\left(1,6\right)}\right]}^T$. Consequently, we have:

$$\begin{array}{ll}\frac{{\mathit{e}}_{t,k}}{{\Vert {\mathit{e}}_{t,k}\Vert }_F}⊛\frac{{\mathit{m}}_{t,k}^{*}}{{\Vert {\mathit{m}}_{t,k}\Vert }_F}& =\left(\frac{b_{t,k}\left(1\right)b_{t,k}^{\ast }\left(1\right)}{\Vert b_{t,k}\left(1\right)\Vert \Vert b_{t,k}^{\ast }\left(1\right)\Vert }\right)\left(\frac{{\tilde{\mathit{e}}}_{t,k}}{{\Vert {\tilde{\mathit{e}}}_{t,k}\Vert }_F}⊛\frac{{\tilde{\mathit{m}}}_{t,k}^{*}}{{\Vert {\tilde{\mathit{m}}}_{t,k}\Vert }_F}\right)\\ & =\left(\frac{{\tilde{\mathit{e}}}_{t,k}}{{\Vert {\tilde{\mathit{e}}}_{t,k}\Vert }_F}⊛\frac{{\tilde{\mathit{m}}}_{t,k}^{*}}{{\Vert {\tilde{\mathit{m}}}_{t,k}\Vert }_F}\right)\\ & ={\mathit{r}}_{t,k}\end{array}$$

(41)

$$\begin{array}{ll}\frac{{\mathit{e}}_{r,k}}{{\Vert {\mathit{e}}_{r,k}\Vert }_F}⊛\frac{{\mathit{m}}_{r,k}^{*}}{{\Vert {\mathit{m}}_{r,k}\Vert }_F}& =\left(\frac{b_{r,k}\left(1\right)b_{r,k}^{\ast }\left(1\right)}{\Vert b_{r,k}\left(1\right)\Vert \Vert b_{r,k}^{\ast }\left(1\right)\Vert }\right)\left(\frac{{\tilde{\mathit{e}}}_{r,k}}{{\Vert {\tilde{\mathit{e}}}_{r,k}\Vert }_F}⊛\frac{{\tilde{\mathit{m}}}_{r,k}^{*}}{{\Vert {\tilde{\mathit{m}}}_{r,k}\Vert }_F}\right)\\ & =\left(\frac{{\tilde{\mathit{e}}}_{r,k}}{{\Vert {\tilde{\mathit{e}}}_{r,k}\Vert }_F}⊛\frac{{\tilde{\mathit{m}}}_{r,k}^{*}}{{\Vert {\tilde{\mathit{m}}}_{r,k}\Vert }_F}\right)\\ & ={\mathit{r}}_{r,k}\end{array}$$

(42)

From (41) and (42), we can observe that the direction cosine is related to the normalized polarization response vector. To obtain 2D-DOD and 2D-DOA estimation, we need to estimate the normalized polarization response vector first. For the tensor-based TS smoothing approach, we define ${\mathit{J}}_{r,p}={\mathit{I}}_N\otimes {\mathit{i}}_{6,p}\in {ℂ}^{N\times 6N}$, where ${\mathit{i}}_{6,p}\in {ℂ}^{1\times 6}$ denotes the p-th row of the $6\times 6$ identity matrix. It is easy to find that the following rotational invariance principle is established:

$${\mathit{J}}_{r,p}{\mathit{C}}_r{\mathsf{\Psi }}_p={\mathit{J}}_{r,1}{\mathit{C}}_r$$

(43)

where ${\mathsf{\Psi }}_p=diag\left\{{\beta }_{r,1}^{\left(1,p\right)},{\beta }_{r,2}^{\left(1,p\right)},\cdots ,{\beta }_{r,K}^{\left(1,p\right)}\right\}$, and $diag\{\cdot \}$ returns a diagonal matrix. Usually, ${\beta }_{r,k}^{\left(1,p\right)}$ is called the $\left(k,p\right)$-th rotational invariance factor corresponding to the $k$-th column and the $p$-th polarization component of ${\mathit{C}}_r$, and ${\mathsf{\Psi }}_p$ is called the $p$-th rotational invariance factor corresponding to the $p$-th polarization component of ${\mathit{C}}_r$. Replacing ${\mathit{C}}_r$ with its noiseless estimate ${\hat{\mathit{C}}}_r$ yields

$${\mathit{J}}_{r,p}{\hat{\mathit{C}}}_r{\left({\mathsf{\Pi }}_1{\mathit{\Delta }}_2\right)}^{-1}{\mathsf{\Psi }}_p\left({\mathsf{\Pi }}_1{\mathit{\Delta }}_2\right)={\mathit{J}}_{r,1}{\hat{\mathit{C}}}_r$$

(44)

where ${(\cdot )}^{-1}$ denotes the inverse. Equivalently,

$${\left({\mathsf{\Pi }}_1{\mathit{\Delta }}_2\right)}^{-1}{\mathsf{\Psi }}_p\left({\mathsf{\Pi }}_1{\mathit{\Delta }}_2\right)={\left({\mathit{J}}_{r,p}{\hat{\mathit{C}}}_r\right)}^{†}{\mathit{J}}_{r,1}{\hat{\mathit{C}}}_r$$

(45)

where ${(\cdot )}^{-1}$ denotes the pseudo-inverse. Combining the properties ${\mathsf{\Pi }}_1{\mathit{\Delta }}_2={\mathit{\Delta }}_2{\mathsf{\Pi }}_1$ and ${\mathit{\Delta }}_2^{-1}{\mathsf{\Psi }}_p{\mathit{\Delta }}_2={\mathsf{\Psi }}_p$, (45) can be simplified as:

$${\mathsf{\Pi }}_1^{-1}{\mathsf{\Psi }}_p{\mathsf{\Pi }}_1={\left({\mathit{J}}_{r,p}{\hat{\mathit{C}}}_r\right)}^{†}{\mathit{J}}_{r,1}{\hat{\mathit{C}}}_r$$

(46)

Obviously, performing eigen decomposition on ${\left({\mathit{J}}_{r,p}{\hat{\mathit{C}}}_r\right)}^{†}{\mathit{J}}_{r,1}{\hat{\mathit{C}}}_r$, one can obtain the estimates of ${\mathsf{\Pi }}_1$ and ${\mathsf{\Psi }}_p$ (denoted by ${\hat{\mathsf{\Pi }}}_1$ and ${\hat{\mathsf{\Psi }}}_p$, respectively). From the k-th diagonal of ${\hat{\mathsf{\Psi }}}_p$, we can obtain the estimated ${\beta }_{r,k}^{\left(1,p\right)}$ (denoted as ${\widehat{\beta }}_{r,k}^{\left(1,p\right)}$). Right multiplying ${\hat{\mathit{B}}}_t$ with ${\hat{\mathsf{\Pi }}}_1^{-1}$, followed by normalizing the calculated vector, we can obtain the estimate of ${\beta }_{t,k}^{\left(1,p\right)}$ (denoted as ${\widehat{\beta }}_{t,k}^{\left(1,p\right)}$). Thereafter, we calculate:

$$\begin{array}{ccc}{\hat{\mathit{r}}}_{t,k}=\frac{{\widehat{\mathit{e}}}_{t,k}}{{\Vert {\widehat{\mathit{e}}}_{t,k}\Vert }_F}⊛\frac{{\hat{\mathit{m}}}_{t,k}^{*}}{{\Vert {\hat{\mathit{m}}}_{t,k}\Vert }_F}& ,& {\hat{\mathit{r}}}_{r,k}=\frac{{\widehat{\mathit{e}}}_{r,k}}{{\Vert {\widehat{\mathit{e}}}_{r,k}\Vert }_F}⊛\frac{{\hat{\mathit{m}}}_{r,k}^{*}}{{\Vert {\hat{\mathit{m}}}_{r,k}\Vert }_F}\end{array}$$

(47)

where ${\widehat{\mathit{e}}}_{t,k}$, ${\hat{\mathit{m}}}_{t,k}$, ${\widehat{\mathit{e}}}_{r,k}$, and ${\hat{\mathit{m}}}_{r,k}$ denote the estimates of ${\tilde{\mathit{e}}}_{t,k}$, ${\tilde{\mathit{m}}}_{t,k}$, ${\tilde{\mathit{e}}}_{r,k}$, and ${\tilde{\mathit{m}}}_{r,k}$, respectively. Finally, 2D-DOD and 2D-DOA estimation can be accomplished via

$$\begin{array}{ccc}\{\begin{array}{c}{\widehat{\theta }}_{t,k}=\arcsin \sqrt{{\hat{\mathit{r}}}_{t,k}^2\left(2\right)+{\hat{\mathit{r}}}_{t,k}^2\left(1\right)}\\ {\widehat{\varphi }}_{t,k}=\arctan \left\{{\hat{\mathit{r}}}_{t,k}\left(2\right)/{\hat{\mathit{r}}}_{t,k}\left(1\right)\right\}\end{array}& ,& \{\begin{array}{c}{\widehat{\theta }}_{r,k}=\arcsin \sqrt{{\hat{\mathit{r}}}_{r,k}^2\left(2\right)+{\hat{\mathit{r}}}_{r,k}^2\left(1\right)}\\ {\widehat{\varphi }}_{r,k}=\arctan \left\{{\hat{\mathit{r}}}_{r,k}\left(2\right)/{\hat{\mathit{r}}}_{r,k}\left(1\right)\right\}\end{array}\end{array}$$

(48)

Since the estimated factor matrices share the same permutation matrix, the estimated 2D-DOD and 2D-DOA are automatically paired.

For the tensor-based RS approach, we can construct ${\mathit{J}}_{t,p}={\mathit{I}}_M\otimes {\mathit{i}}_{6,p}\in {ℂ}^{M\times 6M}$; then, we have:

$${\mathsf{\Pi }}_2^{-1}{\mathsf{Ф}}_p{\mathsf{\Pi }}_2={\left({\mathit{J}}_{t,p}{\hat{\mathit{C}}}_t\right)}^{†}{\mathit{J}}_{t,1}{\hat{\mathit{C}}}_t$$

(49)

where ${\mathsf{Ф}}_p=diag\left\{{\beta }_{t,1}^{\left(1,p\right)},{\beta }_{t,2}^{\left(1,p\right)},\cdots ,{\beta }_{t,K}^{\left(1,p\right)}\right\}$, which can be interpreted as the $p$-th rotational invariance factor corresponding to the $p$-th polarization component of ${\mathit{C}}_t$. Likewise, the eigen decomposition of ${\left({\mathit{J}}_{t,p}{\hat{\mathit{C}}}_t\right)}^{†}{\mathit{J}}_{t,1}{\hat{\mathit{C}}}_t$ reveals the estimates of ${\mathsf{\Pi }}_2$ and ${\mathsf{Ф}}_p$ (denoted by ${\hat{\mathsf{\Pi }}}_2$ and ${\widehat{\mathsf{Ф}}}_p$, respectively). From the k-th diagonal of ${\widehat{\mathsf{Ф}}}_p$, we can obtain the estimated ${\beta }_{t,k}^{\left(1,p\right)}$ (denoted as ${\widehat{\beta }}_{r,k}^{\left(1,p\right)}$). Right multiplying ${\hat{\mathit{B}}}_r$ with ${\hat{\mathsf{\Pi }}}_2^{-1}$ yields the estimate of ${\beta }_{r,k}^{\left(1,p\right)}$ (denoted as ${\widehat{\beta }}_{r,k}^{\left(1,p\right)}$). Repeating the steps in (47) and (48), we the achieve the automatically paired 2D-DOD and 2D-DOA.

For the tensor-based TRS approach, since we can obtain the estimates of the polarization vector matrices, ${\beta }_{t,k}^{\left(1,p\right)}$ and ${\beta }_{r,k}^{\left(1,p\right)}$ can be easily obtained via the normalizations of ${\hat{\mathit{B}}}_t^{″}$ and ${\hat{\mathit{B}}}_r^{″}$. In a similar way, automatically paired 2D-DOD and 2D-DOA can be easily obtained via repeating the steps in (47) and (48).

3.4. 2D-TPA and 2D-RPA Estimation

After we have achieved the 2D-DOD and the 2D-DOA, we can construct the matrices ${\hat{\mathit{F}}}_{t,k}$ and ${\hat{\mathit{F}}}_{r,k}$ according to (10). Moreover, the polarization response vectors ${\mathit{b}}_{t,k}$ and ${\mathit{b}}_{r,k}$ can be easily estimated from the factor matrices (denoting the corresponding estimates as ${\hat{\mathit{b}}}_{t,k}$ and ${\hat{\mathit{b}}}_{r,k}$, respectively). Thereafter, we calculate:

$$\{\begin{array}{c}{\hat{\mathit{v}}}_{t,k}={\hat{\mathit{F}}}_{t,k}^{†}{\hat{\mathit{b}}}_{t,k}\\ {\hat{\mathit{v}}}_{r,k}={\hat{\mathit{F}}}_{r,k}^{†}{\hat{\mathit{b}}}_{r,k}\end{array}$$

(50)

Finally, ${\gamma }_{t,k}$, ${\eta }_{t,k}$, ${\gamma }_{r,k}$, and ${\eta }_{r,k}$ can be estimated from

$$\begin{array}{ccc}\{\begin{array}{l}{\widehat{\gamma }}_{t,k}=\arctan \left\{{\hat{\mathit{v}}}_{t,k}\left(2\right)/{\hat{\mathit{v}}}_{t,k}\left(1\right)\right\}\\ {\widehat{\eta }}_{t,k}=\mathrm{angle}\left\{{\hat{\mathit{v}}}_{t,k}\left(2\right)/{\hat{\mathit{v}}}_{t,k}\left(1\right)\right\}\end{array}& ,& \{\begin{array}{l}{\widehat{\gamma }}_{r,k}=\arctan \left\{{\hat{\mathit{v}}}_{r,k}\left(2\right)/{\hat{\mathit{v}}}_{r,k}\left(1\right)\right\}\\ {\widehat{\eta }}_{r,k}=\mathrm{angle}\left\{{\hat{\mathit{v}}}_{r,k}\left(2\right)/{\hat{\mathit{v}}}_{r,k}\left(1\right)\right\}\end{array}\end{array}$$

(51)

where $\mathrm{angle}\{\cdot \}$ returns the phase, in radian. In addition, the estimated polarization parameters are one-to-one pairing with the estimated 2D-DOD and 2D-DOA.

Up to now, we have finished the descriptions of all the proposed approaches. The algorithmic steps of the tensor-based TS approach (marked as T-TS), the RS approach (marked as T-RS), and the TRS approach (marked as T-TRS) are summarized in Algorithm 1, Algorithm 2, and Algorithm 3, respectively.

Algorithm 1: Algorithmic steps of the T-TS approach

Algorithm 2: Algorithmic steps of the T-RS approach

Algorithm 3: Algorithmic steps of the T-TRS approach

4. Algorithm Analysis

4.1. Complexity Analysis

The main complexities of the proposed estimators are the TALS for PARAFAC decomposition. For an $M\times N\times L$ PARAFAC tensor with rank-K, the TALS procedure requires $l\left[MNLK+K^2\left(ML+NL+MN\right)\right]$, where $l$ represents the number of iterations required to achieve convergence. According to (22)–(24), the dimensions of the PARAFAC models associated with the T-TS approach, the T-RS approach, and the T-TRS approach are $6\times 6N\times 36N$, $6M\times 6\times 36M$, and $6\times 6\times 36$, respectively. Therefore, the main complexities of the T-TS approach, the T-RS approach, and the T-TRS approach are $l\left[6^4{N}^2+2K^2\left(6^{3}N+6^3{N}^2\right)\right]$, $l\left[6^4{M}^2+2K^2\left(6^{3}M+6^3{M}^2\right)\right]$, and $l\left[6^4+6^3{K}^2\right]$, respectively. In comparison, the dominant complexity of the TS approach, the RS approach, and the TRS approach in [37] is the eigen decomposition of the covariance matrix. Accordingly, the dominant complexity of the TS approach, the RS approach, and the TRS approach are $O\left\{{36}^3{N}^3\right\}$, $O\left\{{36}^3{M}^3\right\}$, and $O\left\{{36}^3\right\}$, respectively. We list the main complexity comparison in

Table 1. Complexity and identifiability comparison.

Algorithm	Computing Complexity	Identifiability
TS	$O\left\{{36}^3{N}^3\right\}$	$6N$
RS	$O\left\{{36}^3{M}^3\right\}$	$6M$
TRS	$O\left\{{36}^3\right\}$	36
T-TS	$l\left[6^4{N}^2+2K^2\left(6^{3}N+6^3{N}^2\right)\right]$	$21N+2$
T-RS	$l\left[6^4{M}^2+2K^2\left(6^{3}M+6^3{M}^2\right)\right]$	$21M+2$
T-TRS	$l\left[6^4+6^3{K}^2\right]$	23

. Notably, the main computational burden of the matrix-based approaches is on the third-order of the dimension of the covariance matrix, while the complexity may be lower in the tensor-based counterparts (since only the second-order of the Tx/Rx scale is involved), especially in the presence of large scale Tx/Rx EMVS array.

4.2. Identifiability Analysis

The identifiability (maximum K) is another important index to evaluate the estimator. The identifiability of a PARAFAC-model-based approach is usually constrained by the uniqueness condition of the PARAFAC condition, which is presented in (29). As a result, the uniqueness conditions with respect to T-TS, T-RS, and T-TRS are $K\le 21N+2$, $K\le 21M+2$, and $K\le 23$, respectively. In Table 1, we list the identifiabilities of the proposed estimators. In addition, the identifiabilities corresponding to TS, RS, and TRS are added. It is observed from Table 1 that the proposed T-TS and T-RS methods have better identifiability than the matrix-based ones, while the T-TRS offers worse identifiability than the TRS approach.

4.3. CRB

The deterministic CRB on 2D-DOD, 2D-DOA, 2D-TPA, and 2D-TRA are given by [32]

$$CRB=\frac{{\sigma }^2}{2L}{\left[\mathrm{Re}\left\{\left({\tilde{\mathit{D}}}^H{\mathit{\Gamma }}_{\mathit{C}}^{\perp }\tilde{\mathit{D}}\right)\oplus \left({\mathit{R}}_s\otimes 1_{8\times 8}\right)\right\}\right]}^{-1}$$

(52)

where $\tilde{\mathit{D}}=\left[\frac{\partial c_1}{\partial {\theta }_{t,1}},\cdots ,\frac{\partial c_K}{\partial {\theta }_{r,K}},\frac{\partial c_1}{\partial {\phi }_{t,1}},\cdots ,\frac{\partial c_K}{\partial {\eta }_{r,K}}\right]\in {ℂ}^{36MN\times 8K}$ with $\mathit{C}={\mathit{C}}_t\odot {\mathit{C}}_r$, ${\mathit{c}}_k$ denotes the k-th column of $\mathit{C}$; ${\mathit{\Gamma }}_{\mathit{C}}=I_{MN}-\mathit{C}{\mathit{C}}^{†}$; $\oplus$ denotes the Hadamard product; and $1_{8\times 8}$ denotes the $8\times 8$ matrix when all the entities are ones.

5. Simulation Results

In this section, Q = 200 Monte-Carlo trials are carried out to verify the effectiveness of our estimators. Suppose that there exist $K$ far-field targets, and we assume that L snapshots are collected: the targets’ reflected coefficients satisfy a standard complex Gaussian distribution. Herein, we consider two different array configurations scenarios: Case (I) is the arbitrary geometries. The three-dimensional coordinates $(x_{t,m},y_{t,m},z_{t,m})$ of each Tx sensor are randomly chosen from a uniform distribution $unif\left(-{\kappa }_t,{\kappa }_t\right)$, while the three-dimensional coordinates $(x_{r,n},y_{r,n},z_{r,n})$ of each Rx sensor are randomly chosen from another uniform distribution $unif\left(-{\kappa }_r,{\kappa }_r\right)$. Case (II) is a ULA geometry. Both the Tx array and the Rx array are ULAs along the z-axis with half-wavelength spacing. Simulations are performed on a workstation with an Intel(R) Xeon(R) Silver 4214 CPU and 128 GB RAM, and MATLAB 2018b is adopted. The signal-to-noise ratio (SNR) is defined as the power ratio of $\mathit{x}\left(t\right)-\mathit{n}\left(t\right)$ to $\mathit{n}\left(t\right)$, and the two vectors are defined in (1). The estimation accuracy of the parameter $\vartheta$ is evaluated by the root mean square error (RMSE) defined as:

$$RMSE=\sqrt{1/Q{\displaystyle {\sum }_{q=1}^Q{\left({\widehat{\vartheta }}_q-\vartheta \right)}^2}}$$

(53)

where ${\widehat{\vartheta }}_q$ denotes the estimated $\vartheta$ in the q-th trial.

Example 1.

Scatter results of the proposed approaches. In Figure 2, we plot the scatter results of the proposed T-TS approach, the T-RS approach, as well as the T-TRS approach for completely coherent targets in Case (I). Herein, we consider $K=3$ complete coherent targets with ${\theta }_t=({10}^{°},{36}^{°},{78}^{°})$, ${\varphi }_t=({15}^{°},{15}^{°},{15}^{°})$, ${\gamma }_t=({10}^{°},{52}^{°},{75}^{°})$, ${\eta }_t=({18}^{°},-{31}^{°},{67}^{°})$, ${\theta }_r=(5^{°},{45}^{°},{78}^{°})$, ${\varphi }_r=({15}^{°},-23^{°},{30}^{°})$, ${\gamma }_r=({32}^{°},{10}^{°},{65}^{°})$, and ${\eta }_r=({16}^{°},{56}^{°},{38}^{°})$. Moreover, M = 8, N = 12, L = 500, ${\kappa }_t=2\lambda$, ${\kappa }_r=4\lambda$, and SNR = 0 dB are considered. It is obvious that 2D-DOD, 2D-DOA, 2D-TPA, and 2D-RPA can be accurately estimated by the T-TS, T-RS, and T-TRS approaches, and all the results are correctly paired. It is evident that the proposed approaches are effective for EMVS-MIMO radar, even through the two targets have the same Tx azimuth angle ${\varphi }_t$.

Example 2.

RMSE of the proposed approaches versus SNR.Figure 3 gives the average RMSE on direction (2D-DOD and 2D-DOA) estimation and polarization (2D-TPA and 2D-RPA) estimation versus SNR comparison for $K=3$ completely coherent targets in Case (I), where ${\varphi }_t=({15}^{°},{50}^{°},-{53}^{°})$, and other conditions are the same as in Example I. The average RMSE on direction and polarization are marked with suffixes ‘-d’ and ‘-p’, respectively. It is observed that all the approaches provide improved RMSE with the increasing SNR. Obviously, the TS approach and the RS approach offer more accuate direction estimations than the T-TS approach and the T-RS approach, while the proposed T-TRS approach achieves better estimation accuracy than the TRS approach. Interestingly, all the tensor-aware estimators offer much better polarization estimation accuracy than their matrix counterparts. In addition, it should be noticed that there are gaps between all the mentioned approaches and the CRB, implying that there is room to further improve the estimation accuracy.

Example 3.

RMSE comparisons of various algorithms with ULA geometry. Herein, the performances with respect to the PARAFAC estimator in [32] and the ESPRIT-like algorithm in [34] are added. Unless otherwise specified, the simulation conditions are the same as in Example 2. Figure 4 (upper left) displays the average RMSE on direction estimation versus the SNR, from which we can observe that both ESPRIT-like and traditional PARAFAC algorithms fail to work in the presence of coherent targets. In addition, the tensor-based smoothing approaches outperform the matrix-based ones. Figure 4 (upper right) presents the average RMSE versus correlation coefficient between targets, where the SNR is set to 10 dB. As expected, both ESPRIT-like and traditional PARAFAC algorithms have very close (or even better) RMSE performances to the smoothing approaches. However, they would achieve decreased RMSE with the increasing α once α is larger than a given threshold. However, RMSE of the smoothing-based approaches barely change during the entire α regions. Figure 4 (lower left) shows the RMSE results of various estimators versus sensor position error factor δ, where random position error unif ([0, δλ]) is considered in both the Tx and Rx arrays, and the SNR is set to 20 dB. Obvious, the RMSE of all the estimators barely change with δ, except the traditional PARAFAC algorithm. It is evident that the proposed estimators are insensitive to sensor position errors.

Example 4.

RMSE comparisons versus snapshot L in case (I).Figure 5 presents the average RMSE on direction estimation and polarization estimation at different snapshot numbers L in Case (I), where the SNR is fixed at 10 dB and other conditions are the same as in Example 2. The observations in Figure 5 are similar to that in Figure 3. The proposed tensor-based approaches achieve better RMSE than the matrix-based approaches in polarization parameter estimation, but they may perform worse than the latter in direction parameter estimation.

Example 5.

Performance comparisons versus Rx EMVS number N in case (I).Figure 6 illustrates the average RMSE comparison and average running time (ART) comparison for completely coherent targets in Case (I) with different Rx EMVS number N, where M = 12, ${\kappa }_t=2\lambda$, ${\kappa }_r=16\lambda$, L = 500, and SNR = 10 dB are considered. In such a configuration, the MIMO radar system can be interpreted as a massive MIMO radar since there are 72 Tx components and 6N Rx components. It is shown that the tensor approaches offer a much closer direction estimation performance than the matrix approaches, and they outperform the matrix counterparts with respect to polarization parameter estimation. As expected, the T-TS approach is much more efficient than the TS approach, especially in the presence of large N. However, since M is relatively small, the T-RS approach and the T-TRS approach would require more calculation time than the RS approach and the TRS approach, respectively. In sum, there is a trade-off between running time and the estimation accuracy. In practice, the approprite algorithm should be chosen to balance the two indices.

Example 6.

RMSE comparison versus target distance in case (I).Figure 7 illustrates the average RMSE comparison for completely coherent targets in Case (I) with different inter-target distance $\Delta$, where K = 2 complete coherent targets are considered with parameters ${\theta }_t=({10}^{°},{10}^{°}+\Delta )$, ${\varphi }_t=({15}^{°},{15}^{°}+\Delta )$, ${\gamma }_t=({10}^{°},{10}^{°}+\Delta )$, ${\eta }_t=({16}^{°},{16}^{°}+\Delta )$, ${\theta }_r=({12}^{°},{12}^{°}+\Delta )$, ${\varphi }_r=({11}^{°},{11}^{°}+\Delta )$, ${\gamma }_r=(2^{°},2^{°}+\Delta )$ and ${\eta }_r=(7^{°},7^{°}+\Delta )$. Moreover, for L = 200 and SNR = 10 dB, other conditions are the same as in Example 1. It is observed that the RMSE of all the methods will improved quickly with the increasing $\Delta$ when $\Delta$ is smaller than ${10}^{\circ }$; otherwise the RMSE of all the methods would be improved slowly with the increasing $\Delta$. Since the proposed approaches can take advantage of the tensor nature, they outperform the matrix-based estimators, especially in polarization parameters estimation.

6. Conclusions

In this paper, three estimators based on PARAFAC decomposition are introduced for 2D-DOD, 2D-DOA, 2D-TPA, and 2D-RPA estimation in a bistatic EMVS-MIMO radar with coherent targets. The core idea of the proposed estimators is to de-correlate the coherent targets via spatial smoothing and then formulate the covariance matrix into a third-order PARAFAC decomposition model. After performing PARAFAC decomposition, the factor matrices with respect to various PARAFAC models are obtained. Then, the 2D-DOD and 2D-DOA are achieved via the normalized vector cross-product technique. Finally, the 2D-TPA and 2D-RPA are estimated via the LS technique. Simulations coincide with the theorical results. In the near future, attention should be paid to developing more effective estimation algorithms.