Angle Estimation for MIMO Radar in the Presence of Gain-Phase Errors with One Instrumental Tx/Rx Sensor: A Theoretical and Numerical Study

Abstract

Ideal transmitting and receiving (Tx/Rx) array response is always desirable in multiple-input multiple-output (MIMO) radar. In practice, nevertheless, Tx/Rx arrays may be susceptible to unknown gain-phase errors (GPE) and yield seriously decreased positioning accuracy. This paper focuses on the direction-of-departure (DOD) and direction-of-arrival (DOA) problem in bistatic MIMO radar with unknown gain-phase errors (GPE). A novel parallel factor (PARAFAC) estimator is proposed. The factor matrices containing DOD and DOA are firstly obtained via PARAFAC decomposition. One DOD-DOA pair estimation is then accomplished from the spectrum searching. Thereafter, the remainder DOD and DOA are achieved by the least squares technique with the previous estimated angle pair. The proposed estimator is analyzed in detail. It only requires one instrumental Tx/Rx sensor, and it outperforms the state-of-the-art algorithms. Numerical simulations verify the theoretical advantages.

1. Introduction

Array radar is one of the most important techniques in remote sensing [1,2,3,4,5]. As a branch of active array radar, multiple-input multiple-output (MIMO) radar has gained extensive interest owing to the capability of waveform diversity (it emits mutual orthogonal waveforms using sensor array and receives the echoes using sensor array). To accurately estimate the direction-of-departure (DOD) and direction-of-arrival (DOA) is crucial to MIMO radar. In the past decades, many algorithms have emerged. Typical estimators include multiple signal classification (MUSIC) [6], Capon [7], estimation of signal parameters via rotational invariance techniques (ESPRIT) [8,9], parallel factor (PARAFAC) [10], etc. As reported in the literature, their excellent performances rely on ideal sensor response, i.e., transmitting and receiving (Tx/Rx) sensors are well-calibrated. In practice, unfortunately, unknown gain-phase errors (GPE) may exist in Tx/Rx arrays owing to hardware nonideality, which would bring perturbation to direction matrix, leading to a decline in estimation performance.

The GPE self-calibration issue has been extensively discussed in sensor arrays, but such a topic has always been overlooked in MIMO radar. Usually, the calibration issue in MIMO radar is much more complex than that in sensor array, because it often involves high-dimension array data processing. Several self-calibration methodologies have been proposed for MIMO radar with unknown GPE. In [11], a MUSIC-like framework was developed, which only needs one well-calibrated Tx/Rx sensor. The angles and GPE vectors are achieved via an alternative iteration method. Similar to the traditional counterparts, MUSIC-like is computationally unacceptable owing to exhaustive parameter search and slow convergence performance. Besides, it may occasionally fail since it usually gets stuck in the local minima. In [12], an ESPRIT-like method was introduced. The problem of DOD, DOA, and GPE vectors estimation is linked to a constrained optimization issue, and it is solved by Lagrange multiplier approach. It offers closed-form solutions to DOA and DOD estimation, so it is much more efficient than MUSIC-like. Another reduced-dimension MUSIC (RD-MUSIC) approach was derived in [13]. Therein, the estimation problem was transformed into a quadratic optimization issue, and DOD and DOA were separately obtained via two one-dimensional spatial spectrum searches. In [14], an improved ESPRIT method was presented, where rotational invariance property was only applicable to the data collected from the well-calibrated sensors. Therefore, the virtual aperture of MIMO radar was sacrificed. Recently, tensor algebra has brought new insight into this area. Benefiting from the tensor gain, the tensor methods usually provide more accurate decomposition results than the matrix approaches. By exploiting the tensor nature, the PARAFAC-like estimators were introduced in [15,16]. In both works, the factor matrices containing DOD and DOA were firstly achieved via PARAFAC decomposition. To estimate the angles and GPE vectors from the factor matrices, the Lagrange multiplier method was followed in [15], while the element-wise division operation was chosen in [16]. Another PARAFAC estimator was introduced in [17], which first estimates the gain error via element-wise division, and then it obtains the angle estimation via the least squares fitting. Finally, it achieves the phase errors with the estimated angles. More recently, an improved PARAFAC model was established in [18], which formulates the errors into an isolated factor matrix, and the iterative method is adopted to estimate the factor matrix. Unfortunately, the iteration is very sensitive to the initial values. Besides, the polarimetric signal model was investigated in [19], which calibrates the errors using an additional polarization resource. It should be pointed out that the methods in [12,13,14,16,17] require at least two auxiliary Tx/Rx sensors. Although the estimator in [15] works well with one well-calibrated Tx/Rx sensor, the array manifold must be strictly nonlinear, which limits its implementation.

To avoid the drawbacks pointed out previously, an improved PARAFAC-based estimator is proposed in this paper; the novelties of the proposed estimator are listed as follows:

• Unlike the existing frameworks, the proposed estimator is suitable for MIMO radar with only one instrumental Tx/Rx sensor. This improvement benefits from the fact that the stochastic feature of the phase error is taken into account in the proposed estimator. Moreover, it is adaptive to arbitrary Tx/Rx sensor geometries;

• The proposed estimator is computationally friendly. The DOD/DOA estimation in the proposed estimator can be accomplished via the combination of one-dimensional grid searching and least squares (LS) fitting. It does not involve eigen decomposition or high-dimension spectrum searching.

Owing to the above reasons, the proposed estimator is attractive from the perspectives of robustness, complexity, and accuracy. The identification ability as well as Cramer-Rao bound (CRB) of the proposed algorithm is analyzed. Numerical simulations are designed to verify the improvement of the proposed algorithm.

2. Problem Formulation

Throughout this paper, all the notations are the same to that in [8], unless otherwise specified.

We consider a bistatic MIMO radar configuration, which is equipped with M-element transmitting sensors and N-element receiving sensors. Both Tx/Rx arrays are uniform linear arrays (ULA) with the inter-element spacing d. We suppose the transmitting sensors simultaneously emit M orthogonal pulse waveforms. Assume there are K uncorrelated far field slow-moving targets appearing in the same range bin, let ${\varphi }_k$ and ${\theta }_k$ be the DOD and DOA of the k-th ($k=1,2,\cdots ,K$) target, respectively. The echoes of the targets are collected by the receiving array and passed through the matched filters. The noisy array measurement after matched filtering is given by [10].

(1)$x\left(\tau \right)={\displaystyle \sum _{k=1}^K\left[{\mathit{a}}_r\left({\theta }_k\right)\otimes {\mathit{a}}_t\left({\varphi }_k\right)\right]b_k\left(\tau \right)+\mathit{n}\left(\tau \right)}$

where $\tau$ denotes the snapshot index, and $b_k\left(\tau \right)$ denotes the complex envelope of the k-th target. $\mathit{n}\left(\tau \right)$ denotes the array noise, which is assumed to be Gaussian white. ${\mathit{a}}_t\left({\varphi }_k\right)\in {ℂ}^{M\times 1}$ and ${\mathit{a}}_r\left({\theta }_k\right)\in {ℂ}^{N\times 1}$ are the k-th transmitting response vector (the phase shifts with respect to the reference sensor) and the receiving response vector, which are given by

(2)$\{\begin{array}{c}{\mathit{a}}_t\left({\varphi }_k\right)={\left[1,\exp \left(-\mathrm{j}2\pi d\sin \left({\varphi }_k\right)/\lambda \right),\cdots ,\exp \left(-\mathrm{j}2\pi d\left(M-1\right)\sin \left({\varphi }_k\right)/\lambda \right)\right]}^T\\ {\mathit{a}}_r\left({\theta }_k\right)={\left[1,\exp \left(-\mathrm{j}2\pi d\sin \left({\theta }_k\right)/\lambda \right),\cdots ,\exp \left(-\mathrm{j}2\pi d\left(N-1\right)\sin \left({\theta }_k\right)/\lambda \right)\right]}^T\end{array}$

where $\lambda$ is the wavelength of the carrier frequency. Specifically, suppose that GPE exists in both the Tx/Rx arrays. Then the model in (1) should be modified into [12,13,14,15,16,17]

(3)$x\left(\tau \right)={\displaystyle \sum _{k=1}^K\left[\left({\mathit{C}}_r{\mathit{a}}_r\left({\theta }_k\right)\otimes {\mathit{C}}_t{\mathit{a}}_t\left({\varphi }_k\right)\right)\right]b_k\left(\tau \right)+\mathit{n}\left(\tau \right)}$

where ${\mathit{C}}_t={\mathit{G}}_t{\mathit{P}}_t\in {ℂ}^{M\times M}$ and ${\mathit{C}}_r={\mathit{G}}_r{\mathit{P}}_r\in {ℂ}^{N\times N}$ are two diagonal matrices, whose diagonal elements denote the GPE of the transmitting array and receiving array, respectively. ${\mathit{G}}_t=\mathrm{diag}\left\{{\mathit{g}}_t\right\}$,$P_t=\mathrm{diag}\left\{\exp \left(j{\mathit{p}}_t\right)\right\}$, ${\mathit{G}}_r=\mathrm{diag}\left\{{\mathit{g}}_r\right\}$,${\mathit{P}}_r=\mathrm{diag}\left\{\exp \left(j{\mathit{p}}_r\right)\right\}$,${\mathit{g}}_t={\left[g_{t1},g_{t2},\cdots ,g_{tM}\right]}^T$,${\mathit{p}}_t={\left[p_{t1},p_{t2},\cdots ,p_{tM}\right]}^T$, ${\mathit{g}}_r={\left[g_{r1},g_{r2},\cdots ,g_{rN}\right]}^T$,${\mathit{p}}_r={\left[p_{r1},p_{r2},\cdots ,p_{rN}\right]}^T$, where $g_{tm}\left(m=1,2,\cdots ,M\right)$ and $g_{tn}\left(n=1,2,\cdots ,N\right)$ denote the gain errors associated with the m-th transmitting antenna and the n-th receiving antenna, and $p_{tm}$ and $p_{rn}$ denote the associated phase errors, respectively. Herein, the reference Tx/Rx sensor is assumed to be well-calibrated, e.g., $g_{t1}=g_{r1}=1$, $p_{t1}=p_{r1}=0$. Once L $\left(\tau =1,2,\cdots ,L\right)$ snapshots are collected, (3) can be formulated in matrix format as

(4)$\begin{array}{c}\mathit{X}=\left[\left({\mathit{C}}_r{\mathit{A}}_r\right)\odot \left({\mathit{C}}_t{\mathit{A}}_t\right)\right]{\mathit{B}}^T+{\mathit{N}}_x\\ =\left[{\tilde{\mathit{A}}}_r\odot {\tilde{\mathit{A}}}_t\right]{\mathit{B}}^T+{\mathit{N}}_x\end{array}$

where ${\mathit{A}}_t=\left[{\mathit{a}}_t\left({\varphi }_1\right),{\mathit{a}}_t\left({\varphi }_2\right)\cdots ,{\mathit{a}}_t\left({\varphi }_K\right)\right]\in {ℂ}^{M\times K}$ denotes the transmitting direction matrix; ${\mathit{A}}_r=\left[{\mathit{a}}_r\left({\theta }_1\right),{\mathit{a}}_r\left({\theta }_2\right)\cdots ,{\mathit{a}}_r\left({\theta }_K\right)\right]\in {ℂ}^{N\times K}$ denotes the receiving direction matrix; $\mathit{B}=\left[{\mathit{b}}_1,{\mathit{b}}_2,\cdots ,{\mathit{b}}_K\right]\in {ℂ}^{L\times K}$ denotes the target reflection coefficient matrix; and ${\mathit{b}}_k={\left[b_k\left(1\right),b_k\left(2\right),\cdots ,b_k\left(L\right)\right]}^T$, ${\mathit{N}}_x=\left[\mathit{n}\left(1\right),\mathit{n}\left(2\right),\cdots ,\mathit{n}\left(L\right)\right]$, ${\tilde{\mathit{A}}}_t={\mathit{C}}_t{\mathit{A}}_t$ and ${\tilde{\mathit{A}}}_r={\mathit{C}}_r{\mathit{A}}_r$ denote the associated perturbed direction matrices. Since the original direction matrices ${\mathit{A}}_t$ and ${\mathit{A}}_t$ are corrupted by the GPE matrices, thus the traditional subspace algorithms will fail to work.

It was pointed out that $\mathit{X}$ can be formulated into a third-order tensor $\mathcal{X}\in {ℂ}^{N\times M\times L}$ as [10]

(5)$\mathcal{X}={\displaystyle \sum }_{k=1}^K{\tilde{\mathrm{a}}}_{r,k}◦{\tilde{\mathrm{a}}}_{t,k}◦{\mathit{b}}_k+\mathcal{N}$

where $\mathcal{N}$ is the tensor format of the noise. According to [20] $\mathit{X}$ can be interpreted as the Mode-3 unfolding of $\mathcal{X}$ Similarly, we can achieve another two matrices $\mathit{Y}\in {ℂ}^{ML\times N}$ and $\mathit{Z}\in {ℂ}^{NL\times L}$, respectively, via Mode-1 and Mode-2 unfolding of $\mathcal{X}$ , which are given by

(6)$\mathit{Y}=\left[{\tilde{\mathit{A}}}_t\odot \mathit{B}\right]{\tilde{\mathit{A}}}_r^T+{\mathit{N}}_y$

and

(7)$\mathit{Z}=\left[\mathit{B}\odot {\tilde{\mathit{A}}}_r\right]{\tilde{\mathit{A}}}_t^T+{\mathit{N}}_z$

where ${\mathit{N}}_y$ and ${\mathit{N}}_z$ denote the rearranged array noise.

3. The Proposed Framework

The main steps of the proposed estimator are shown in Figure 1. It performs PARAFAC decomposition to acquire the corrupted direction matrices. Thereafter, it acquires the DOD/DOA via the combination of grid searching and the LS method. The details of the proposed approach will be shown next.

3.1. Estimation of the Corrupted Direction Matrices

Since the GPE matrices ${\mathit{C}}_t\in {ℂ}^{M\times M}$ and ${\mathit{C}}_r\in {ℂ}^{N\times N}$ are coupled with the direction matrices ${\mathit{A}}_t$ and ${\mathit{A}}_r$, we need to achieve the corrupted response matrices ${\tilde{\mathit{A}}}_t$ and ${\tilde{\mathit{A}}}_r$ first. From the model in (5), they can be estimated via optimizing

(8) ${\underset{{\tilde{\mathit{A}}}_t,{\tilde{\mathit{A}}}_r,\mathit{B}}{\min }\Vert X-{\mathcal{I}}_{\times 1}{\tilde{\mathit{A}}}_{r\times 2}{\tilde{\mathit{A}}}_{t\times 3}\mathit{B}\Vert }_F$

To tackle the above optimization issue, the trilinear alternating least squares (TALS) technique is commonly utilized. Since the tensor can be unfolded into matrices, the optimization in (8) is equivalently expressed as

(9)$\{\begin{array}{c}\underset{{\tilde{\mathit{A}}}_t,{\tilde{\mathit{A}}}_r,\mathit{B}}{\min }{\Vert \mathit{X}-\left[{\tilde{\mathit{A}}}_r\odot {\tilde{\mathit{A}}}_t\right]{\mathit{B}}^T\Vert }_F\\ \underset{{\tilde{\mathit{A}}}_t,{\tilde{\mathit{A}}}_r,\mathit{B}}{\min }{\Vert \mathit{Y}-\left[{\tilde{\mathit{A}}}_t\odot \mathit{B}\right]{\tilde{\mathit{A}}}_r^T\Vert }_F\\ \underset{{\stackrel{˜˜}{\mathit{A}}}_t,{\stackrel{˜˜}{\mathit{A}}}_r,\mathit{B}}{\min }{\Vert \mathit{Z}-\left[\mathit{B}\odot {\tilde{\mathit{A}}}_r\right]{\tilde{\mathit{A}}}_t^T\Vert }_F\end{array}$

The main principle of TALS is to successively update one of ${\tilde{\mathit{A}}}_r$,${\tilde{\mathit{A}}}_t$, and $\mathit{B}$ via the LS technique and treat the remainder as known prior. For instance, if ${\tilde{\mathit{A}}}_r$ and ${\tilde{\mathit{A}}}_t$ are known, then $\mathit{B}$ can be estimated via fitting the first item of (9), the LS solution of which is then given by

(10)${\hat{\mathit{B}}}^T={\left[{\tilde{\mathit{A}}}_r\odot {\tilde{\mathit{A}}}_t\right]}^{†}\mathit{X}$

Similarly, once ${\tilde{\mathit{A}}}_t$ and $\mathit{B}$ are estimated, we can further estimate ${\tilde{\mathit{A}}}_r$ via LS fitting the second item of (9), and the LS calculation of ${\tilde{\mathit{A}}}_r$ is

(11)${\hat{\mathit{A}}}_r^T={\left[{\tilde{\mathit{A}}}_t\odot \mathit{B}\right]}^{†}\mathit{Y}$

In addition, with the previously estimated ${\tilde{\mathit{A}}}_r$ and $\mathit{B}$, we can estimate ${\tilde{\mathit{A}}}_t$ via the LS fitting of the last item in (9), from which one can obtain the LS update of ${\tilde{\mathit{A}}}_t$ as

(12)${\hat{\mathit{A}}}_t^T={\left[\mathit{B}\odot {\tilde{\mathit{A}}}_r\right]}^{†}\mathit{Z}$

Usually, the LS updates in (10)–(12) will repeat before the convergence conditions, e.g., the fitting error is smaller than a fixed threshold. Although the TALS is guaranteed to converge, it suffers from the slow convergence speed. Herein, the well-known COMFAC algorithm in [21,22] is adopted to speed the iteration, which undergoes some pre-processing and can be quickly converged with only a few iterations.

Unlike the traditional matrix decomposition, PARAFAC decomposition is usually unique. For the model in (5), it was pointed out that if the Kruskal ranks of ${\tilde{\mathit{A}}}_r$, ${\tilde{\mathit{A}}}_t$ and $\mathit{B}$ (denoted by $kr\left({\tilde{\mathit{A}}}_r\right)$, $kr\left({\tilde{\mathit{A}}}_t\right)$ and $kr\left(\mathit{B}\right)$, respectively) fulfill [10]

(13)$kr\left({\tilde{\mathit{A}}}_r\right)+kr\left({\tilde{\mathit{A}}}_t\right)+kr\left(\mathit{B}\right)\ge 2K+2$

then the estimations of ${\tilde{\mathit{A}}}_r$, ${\tilde{\mathit{A}}}_t$, and $\mathit{B}$ are unique with permutations and scaling effect, i.e.,

(14)$\{\begin{array}{c}{\hat{\mathit{A}}}_t={\tilde{\mathit{A}}}_t\mathit{\Pi }{\mathsf{\Delta }}_1+{\mathit{N}}_1\\ {\hat{\mathit{A}}}_r={\tilde{\mathit{A}}}_r\mathit{\Pi }{\mathsf{\Delta }}_2+{\mathit{N}}_2\\ \hat{\mathit{B}}=\mathit{B}\mathit{\Pi }{\mathsf{\Delta }}_3+{\mathit{N}}_3\end{array}$

where $\mathit{\Pi }$ is a permutation matrix, and ${\mathsf{\Delta }}_1$,${\mathsf{\Delta }}_2$, and ${\mathsf{\Delta }}_3$ denote the corresponding scaler effects, which are diagonal scaling matrices and satisfy ${\mathsf{\Delta }}_1{\mathsf{\Delta }}_2{\mathsf{\Delta }}_3={\mathrm{I}}_k$.${\mathit{N}}_1$, ${\mathit{N}}_2$, and ${\mathit{N}}_3$ are the fitting error matrices.

3.2. DOD and DOA Estimation

In what follows, we will show how to estimate the DOD from ${\hat{\mathit{A}}}_t$, and the DOA can be estimated in a similar way from ${\hat{\mathit{A}}}_r$. Before the detailed derivation, we assume that the phase errors are zero-mean random variables, i.e.,

(15)$\{\begin{array}{c}E\left\{{\rho }_t\right\}=0\\ E\left\{{\rho }_r\right\}=0\end{array}$

The above hypothesis is always held for large scale arrays. It is easy to verify that

(16)$0_1=0_k\mathit{\Pi }$

Besides, let ${\hat{\mathrm{a}}}_{t,k}$ and ${\tilde{\mathrm{a}}}_{t,k}$ stand for, respectively, the k-th column of ${\tilde{\mathit{A}}}_t$ and ${\tilde{\mathit{A}}}_t$. Ignoring the noise item in ${\widehat{\mathit{a}}}_{t,k}$,we can get ${\widehat{\mathit{a}}}_{t,k}=\delta {\widehat{\mathit{a}}}_{t,k}$, where $\delta$ is a scalar. Define $p_m\left(\gamma \right)$ as

(17)$p_m\left(\gamma \right)={\left\{\mathrm{angle}\left({\left[{\widehat{\mathit{a}}}_{t,k}\right]}_m{\left[{\mathit{a}}_t^H\left(\gamma \right)\right]}_m\right)\right\}}^2$

where ${\left[\mathit{a}\right]}_m$ denotes the m-th entity of $\mathit{a}$. Consequently, if $\gamma ={\varphi }_k$, then $p_m\left(\gamma \right)=p_{tm}^2$. Similar to [23], we construct the following spatial spectrum function

(18)$f\left(\gamma \right)={\displaystyle \sum _{m=2}^M{p}_m\left(\gamma \right)}$

Notably, once $\gamma ={\varphi }_k$, we have $f\left(\gamma \right)=p_{t2}^2+p_{t3}^2+\cdots +p_{tM}^2$. Replace $\sin \left(\gamma \right)$ with x, and then $f\left(\gamma \right)$ can be rewritten as

(19)$g\left(x\right)={{\displaystyle \sum _{m=2}^M\left[p_{tm}+2\pi d_m\left(x-\sin \left({\varphi }_k\right)\right)/\lambda \right]}}^2$

where $d_m=d\left(m-1\right)$. To find out the minimum value of $g\left(x\right)$, one can calculate the first derivative of $g\left(x\right)$ and enforce it to zero, and then we can get the minimum point

(20)$x_0=-\frac{\lambda {\displaystyle \sum _{m=2}^M{d}_m{p}_{tm}}}{2\pi {\displaystyle \sum _{m=2}^M{d}_m^2}}+\sin {\varphi }_k$

In addition, let $\tilde{m}$ denote the arithmetic mean of $d_m{p}_{tm}$, i.e.,

(21)$\tilde{m}=\frac{1}{M}{\displaystyle \sum _{m=2}^M{d}_m{p}_{tm}}$

The mathematical expectation of $\tilde{m}$ is

(22)$E\left\{\tilde{m}\right\}=\frac{1}{M}{\displaystyle \sum _{m=2}^{M}E\left\{d_m{p}_{tm}\right\}}=E\left\{d{p}_t\right\}$

It is easy to prove that $\tilde{m}$ is an unbiased estimation. Then the mathematical expectation of $x_0$ is given by

(23)$E\left\{x_0\right\}=-\frac{\lambda \cdot M\cdot E\left\{d{p}_t\right\}}{2\pi {\displaystyle \sum _{m=2}^M{d}_m^2}}+\sin {\varphi }_k$

Since d and $p_t$ are independent, then we have $E\left\{d{p}_t\right\}=E\left\{d\right\}E\left\{p_t\right\}$. In combination with (15), it is easy to have

(24)$E\left\{x_0\right\}=\sin {\varphi }_k$

Therefore, by setting a series of search grids ${\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_Q$, the first DOD can be estimated via finding the minimum value of $f\left(\gamma \right)$, e.g., ${\widehat{\varphi }}_1=\arg \underset{\gamma }{\max }f\left(\gamma \right)$.

Although the rest DOD can be achieved in a similar method, it is computationally inefficient. Fortunately, the columns of ${\tilde{\mathit{A}}}_t$ suffer from the same GPE. Let $\mathit{h}=-\mathrm{angle}\left({\tilde{\mathit{a}}}_{t,p}./{\tilde{\mathit{a}}}_{t,q}\right)$, where $./$ represents the element-wise division [16]. Define

(25)${\mathit{P}}_t={\left[\begin{array}{l}1 1 \cdots 1\\ 0 \pi \cdots \left(M-1\right)\pi \end{array}\right]}^T$

It is easy to find ${\mathit{P}}_t\mathit{v}=\mathit{h}$, where $\mathit{v}\in {ℂ}^{2\times 1}$ is a vector with the second entity ${\left[\mathit{v}\right]}_2$$\sin \left({\varphi }_p\right)-\sin \left({\varphi }_q\right)$. Accordingly, let ${\widehat{\mathit{h}}}_{t,k}=-\mathrm{angle}\left({\widehat{\mathit{a}}}_{t,p}./{\widehat{\mathit{a}}}_{t,q}\right)$, and compute

(26)${\widehat{\mathit{v}}}_{t,k}={\mathit{P}}_t^{†}{\hat{\mathit{h}}}_{t,k}$

After that, we can obtain the remainder (k = 2, 3…, K) DOD via

(27)${\hat{\varphi }}_k=\arcsin \left({\left[\hat{v}\right]}_2+\sin \left({\hat{\varphi }}_1\right)\right)$

Similarly, we can estimate all the DOA via one spatial spectrum search and K-1 LS fittings.

4. Algorithmic Analysis

4.1. Related Remarks

Remark 1.

From (14), one can observe that both ${\hat{\mathit{A}}}_t$ and ${\hat{\mathit{A}}}_t$ suffer from the same permutation $\mathit{\Pi }$ Therefore, the estimated k-th DOD-DOA pair is paired automatically.

Remark 2.

Obviously, the proposed algorithm does not involve array geometry, but it can be easily extended to configuration with arbitrary array geometry, by adjusting the fitting matrix ${\mathit{P}}_t$

.

Remark 3.

To simplify the derivation, we ignore the noise item in ${\hat{\mathit{a}}}_{t,k}$ , so the proposed algorithm works well at high signal-to-noise ratio (SNR), but may fail at low SNR.

Remark 4.

Once the DOD and DOA are obtained, the direction matrices ${\mathit{A}}_t$and ${\mathit{A}}_r$can be constructed. With the element-wise division technique in [16], one can easily achieve the estimated GPE matrices, and consequently the GPE self-calibration can be carried out by left multiplying X with ${\hat{\mathit{C}}}_r^{-1}\otimes {\hat{\mathit{C}}}_t^{-1}$ , where ${\hat{\mathit{C}}}_r$ and ${\hat{\mathit{C}}}_t$ denote the estimations of ${\mathit{C}}_r$ and ${\mathit{C}}_t$, respectively.

Remark 5.

The target reflection coefficient matrix B contains the Doppler information of target. Although the topic of estimate target Doppler frequency is interesting, it is beyond the scope of this paper.

4.2. Identifiability

The identifiability is equal to the maximum value K. (13) reveals the identifiability of the proposed estimator. Since we have $\max \left(kr\left({\tilde{\mathit{A}}}_t\right)\right)=M$, $\max \left(kr\left({\tilde{\mathit{A}}}_r\right)\right)=N$, when B is full Kruskal rank, (13) becomes

(28)$M+N+\min \left(L,K\right)\ge 2\left(K+1\right)$

Generally, we have L ≥ K, so (28) is simplified as $M+N\ge K+2$. It means that the proposed estimator can identify at most $M+N-2$ targets.

4.3. Deterministic CRB

The whole unknown parameter vector can be written as $\mathsf{\xi }={\left[{\mathit{\gamma }}^T,{\mathit{\eta }}^T\right]}^T$, where $\mathit{\gamma }={\left[{\mathit{g}}_t^T,{\mathit{p}}_t^T,{\mathit{g}}_r^T,{\mathit{p}}_r^T\right]}^T\in {ℝ}^{2\left(M+N\right)\times 1}$, $\mathit{\eta }={\left[{\theta }_{t,1},\cdots {\theta }_{t,K},{\varphi }_{t,1},\cdots {\theta }_{r,k},\cdots {\varphi }_{r,K}\right]}^T\in {ℝ}^{4K\times 1}$. According to [12], the CRB on angle estimation is given by

(29)$\mathrm{CRB}\left(\mathit{\eta }\right)=\frac{{\sigma }^2}{2}{\left[{\mathit{J}}_{\eta \eta }-{\mathit{J}}_{\gamma \eta }{\mathit{J}}_{\gamma \gamma }^{-1}{\mathit{J}}_{\gamma \eta }^T\right]}^{-1}$

where ${\mathit{\Pi }}_{\tilde{\mathit{A}}}^{\perp }=\mathit{I}-\tilde{\mathit{A}}{\tilde{\mathit{A}}}^{†}$, $\tilde{\mathit{A}}={\tilde{\mathit{A}}}_r\odot {\tilde{\mathit{A}}}_t$, ${\mathit{J}}_{\eta \eta }={\displaystyle \sum _{\tau =1}^L{\mathit{D}}_{\tau }^H{\mathit{\Pi }}_{\tilde{\mathit{A}}}^{\perp }{\mathit{D}}_{\tau }}$, ${\mathit{J}}_{\gamma \eta }={\displaystyle \sum _{\tau =1}^L{\mathit{D}}_{\tau }^H{\mathit{\Pi }}_{\tilde{\mathit{A}}}^{\perp }{\mathit{H}}_{\tau }}$, ${\mathit{J}}_{\gamma \gamma }={\displaystyle \sum _{\tau =1}^L{\mathit{H}}_{\tau }^H{\mathit{\Pi }}_{\tilde{A}}^{\perp }{\mathit{H}}_{\tau }}$ with

(30)${\mathit{D}}_{\tau }=\left[\frac{\partial {\tilde{\mathit{a}}}_1{b}_1\left(\tau \right)}{\partial {\theta }_{t,1}},\cdots ,\frac{\partial \tilde{\mathit{a}}{b}_k\left(\tau \right)}{\partial {\varphi }_{t,k}}\cdots ,\frac{\partial {\tilde{\mathit{a}}}_k{b}_k\left(\tau \right)}{\partial {\theta }_{r,k}},\cdots ,\frac{\partial {\tilde{\mathit{a}}}_K{b}_K\left(\tau \right)}{\partial {\varphi }_{t,K}}\right]$

(31)${\mathit{H}}_{\tau }=\left[\frac{\partial \mathit{C}}{\partial g_{t1}},\cdots ,\frac{\partial \mathit{C}}{\partial p_{tm}},\cdots ,\frac{\partial \mathit{C}}{\partial g_{rn}},\cdots ,\frac{\partial \mathit{C}}{\partial p_{rn}}\right]blkdi\mathit{a}g\left\{\underset{2\left(M+N\right)}{\underbrace{\mathit{A}b\left(\tau \right)\cdots \mathit{A}b\left(\tau \right)}}\right\}$

where $blkdi\mathit{a}g\{·\}$ denotes a block diagonal matrix.

5. Simulation Results

Herein, Q = 200 Monte-Carlo simulation experiments are carried out to show the improvement of our estimator. Suppose that there are M = 11 transmitting sensors and N = 11 receiving sensors, both of which are half-wavelength spacing ULA. Consider that K = 3 targets are located at $\left({\theta }_1,{\varphi }_1\right)=\left({10}^{\circ },{15}^{\circ }\right)$, $\left({\theta }_2,{\varphi }_2\right)=\left({20}^{\circ },{25}^{\circ }\right)$ and $\left({\theta }_3,{\varphi }_3\right)=\left({30}^{\circ },{35}^{\circ }\right)$, respectively. The following two cases are considered. Case (I): one instrumental Tx/Rx sensor with ${\mathit{g}}_t$=[1, 1.26, 0.62, 1.09, 0.84, 0.56, 0.64, 1.03, 0.93, 1.06, 0.52]^T and ${\mathit{g}}_r$=[1, 0.6, 1.27, 1.34, 0.63, 1.42, 0.57, 1.15, 0.89, 0.75, 1.34]^T; Case (II): three instrumental transmitting sensors and two instrumental receiving sensors with ${\mathit{g}}_t$=[1, 1, 1, 1.07, 0.69, 1.17, 1.15, 1.43, 1.11, 1.3, 0.53]^T and ${\mathit{g}}_r$=[1, 1, 0.6, 1.27, 0.67, 1.17, 0.68, 1.15, 1.41, 0.72, 0.72, 1.31, 0.52, 0.78]^T. In both cases, phase errors are randomly chosen from ±0.1, and L = 100. SNR in the simulation is defined as $\mathrm{SNR}=10\lg \left({\Vert X-N_x\Vert }_F^2/{\Vert N_x\Vert }_F^2\right)$ [dB]. The root mean square error (RMSE) is utilized to assess the estimation accuracy, the definition of which is

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

where ${\vartheta }_q$ is the estimated parameter in the q-th trial, and ${\vartheta }_0$ is the true value. For comparison purpose, the performances of ESPRIT-like [12], RD-MUSIC [13], PARAFAC-like [15], and CRB are added (Matlab demo of the proposed algorithm can be found here).

Firstly, the average RMSE performances in Case I with various SNR are depicted in Figure 2. It is seen that the RMSE performances corresponding to ESPRIT-like, RD-MUSIC, and PARAFAC-like are at least an order of magnitude higher than the proposed estimator, and they barely vary with the increasing SNR. The observation indicates that ESPRIT-like, RD-MUSIC, and PARAFAC-like fail to work in such cases, since they need at least two well-calibrated Tx/Rx sensors. In contrast, average RMSE of the proposed estimator improves when SNR increases, which proves that the proposed estimator is suitable for the scenario with one instrumental Tx/Rx sensor.

Secondly, we plot the average RMSE of the proposed estimator vs. GPE ranges of our estimator in Case I. Herein, we consider the scenario that additional gain error δ is added to each sensor or additional phase error η is added to the mean of the phase errors. Figure 3 indicates the results. Obviously, the proposed estimator is insensitive to the gain error while it is sensitive to the mean of phase error. The mean of the phase errors in the proposed estimator must be strictly zero, otherwise it may fail to work.

Thirdly, we repeat the above simulation in Case II, and the result is shown in Figure 4. It is clearly observed that all the RMSE improve once SNR increases. Compared with ESPRIT-like, RD-MUSIC, and PARAFAC-like, the RMSE of the proposed estimator is about 10 dB lower. The observation may be caused by two reasons. On the one hand, the estimation accuracy of MUSIC is related to its search grid while ESPRIT-like suffers from the inaccurate GPE error matrices estimation. On the other hand, the eigen decomposition in both MUSIC and ESPRIT cannot exploit the tensor nature, leading to degraded subspace estimation. Notably, the proposed estimator always provides the lowest RMSE, which proves that the proposed estimator is superior to all comparison algorithms. However, it should be pointed out, as illustrated in Figure 2 and Figure 4, the distance between the proposed approach and CRB is big, so there is much room for further improvement.

Finally, the average run time of different algorithms versus N is given. In this experiment, M = 21, SNR = 15 dB, snapshot L = 500, gain errors, and phase errors are randomly chosen from [0.5,1.5] and ±0.1, respectively. It can be found in Figure 5 that the run time that RD-MUSIC required is at least 10 dB higher than the proposed estimator. Besides, the run time of our estimator is similar to that of PARAFAC-like (the gap is smaller than 5 dB), and the proposed estimator is more efficient than ESPRIT-like for large scale array.

6. Conclusions and Future Research

We have stressed the issue of joint angle and GPE self-calibration in bistatic MIMO radar system in this paper, where only one well-calibrated Tx/Rx sensor is considered. An improved PARAFAC estimator was proposed, the RMSE of which is an order of magnitude lower than that of the state-of-the-art PARAFAC-like method [15], while they offer comparable complexity. It should be a good choice for engineering application, especially in massive MIMO configuration. In the near future, more efforts should be devoted to verifying the proposed approach in the presence of real-life data.

Author Contributions

F.W.—designed the simulations and provided the draft; J.S. and X.W. and L.W.—helped to improve the writing. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China under grant number 62071476.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Illustration of the proposed estimator.

Figure 2. Average RMSE vs. SNR in Case I.

Figure 3. (a) Average RMSE vs. SNR with various gain error range δ, and (b) Average RMSE vs. SNR with various phase error range η.

Figure 4. Average RMSE vs. SNR in Case II.

Figure 5. CPU run time versus N in Case II.