Off-Grid DOA Estimation Using Sparse Bayesian Learning for MIMO Radar under Impulsive Noise

Abstract

Direction of arrival (DOA) estimation is an essential and fundamental part of array signal processing, which has been widely used in radio monitoring, autonomous driving of vehicles, intelligent navigation, etc. However, it remains a challenge to accurately estimate DOA for multiple-input multiple-output (MIMO) radar in impulsive noise environments. To address this problem, an off-grid DOA estimation method for monostatic MIMO radar is proposed to deal with non-circular signals under impulsive noise. In the proposed method, firstly, based on the property of non-circular signal and array structure, a virtual array output was built and a real-valued sparse representation for the signal model was constructed. Then, an off-grid sparse Bayesian learning (SBL) framework is proposed and further applied to the virtual array to construct novel off-grid sparse model. Finally, off-grid DOA estimation was realized through the solution of the sparse reconstruction with high accuracy even in impulsive noise. Numerous simulations were performed to compare the algorithm with existing methods. Simulation results verify that the proposed off-grid DOA method enables evident performance improvement in terms of accuracy and robustness compared with other works on impulsive noise.

1. Introduction

With the rapid development of wireless communications, such as fifth-generation (5G) communication systems and the upcoming 6G communication systems, various smart technologies have been extended to provide intelligent and effective services to users [1,2,3]. Among these smart technologies, multiple-input multiple-output (MIMO) radar has been extensively applied in many military and civilian fields, such as target detection, intelligent radio monitoring and parameter estimation [4,5,6]. Generally, an MIMO radar system can be divided into two categories: bistatic MIMO radar system and monostatic MIMO radar system. Compared with bistatic MIMO radar, monostatic MIMO radar has larger array aperture and degree of freedom as it is equipped with collocated of phased arrays [4]. Hence, we mainly consider the monostatic MIMO radar system in this paper.

Direction of arrival (DOA) estimation is a fundamental and essential step in target parameter estimation, and it has gained widespread attention in wireless communication and vehicular navigation fields [7]. As MIMO radar can achieve a higher spatial resolution when compared to conventional uniform linear array (ULA), solving DOA estimation with MIMO radar has been a hot topic, especially for monostatic MIMO radar [1,8,9]. Similar to DOA estimation in ULA, subspace based methods are dominated by DOA estimation with MIMO radar, mainly including MUSIC [10], ESPRIT [11], and the parallel factor based method [12]. Nevertheless, these subspace based methods usually suffer performance degradation with less prior knowledge and fewer signal snapshots. Considering these limitations, a promising technique for DOA estimation is further exploited based on sparse representation [13]. Meanwhile, a series of sparse representation-based approaches have been proposed to solve DOA estimation with monostatic MIMO radar. In [14], a reweighted ${\ell }_1$-SVD is utilized for sparse reconstruction to improve the estimation performance. In [15], an improved $l_p$ norm minimization DOA estimation method is proposed to separate two closely spaced sources. Ref. [16] proposes a novel idea of simultaneous source number identification and DOA estimation with high accuracy. Ref. [17] further improved the performance for DOA estimation with a non-circular signal property and nuclear norm minimization. In [18], the iterative reweighted proximal projection is firstly introduced to estimate DOA for MIMO radar with effective performance improvement. Ref. [19] proposed a Bayesian learning based DOA approach which also shows resistance to mutual coupling. Non-circular signals can extend the array aperture with a non-circular signal property [20]. The non-circular signal property is employed to estimate DOA with MIMO radar in [21]. Based on this property, both subspace and sparse representation based methods can improve the angular resolution and estimation accuracy.

Although the aforementioned methods have advantages in MIMO radar DOA estimation, their performance would degrade severely in impulsive noise. This is because the background noise is assumed to be AWGN (additive white Gaussian noise) in these DOA estimation methods. However, it is quite unrealistic to adopt this AWGN assumption to describe the noise in complicated electromagnetic environments. Due to the impact of natural or man-made signal sources, the channel noise usually behaves with non-Gaussian impulsive characteristics [22,23]. Under an impulsive noise environment, the conventional DOA estimation methods, which are based on second-order statistics (SOS) or higher-order statistics (HOS), suffer serious performance degradation as SOS and HOS are not convergent for non-Gaussian impulsive noise. To suppress impulsive noise for DOA estimation, many concepts have been developed such as fractional lower-order statistics [24,25,26], correntropy [27,28], nonlinear functions [29,30] and so on. However, there are few works that focus on solving DOA estimation with MIMO radar based on sparse representation in impulsive noise. Robust and accurate MIMO radar DOA estimation in impulsive noise remains a challenge.

To accurately estimate DOA with MIMO radar in impulsive noise, a novel off-grid DOA estimation method is proposed in this paper. In the proposed method, firstly, a real-valued array model is constructed with the aid of the structure of monostatic MIMO radar and non-circular signals’ property. The constructed real-valued virtual array can effectively reduce the computational complexity for sparse Bayesian learning and extend the array aperture. Then, an off-grid sparse Bayesian learning is proposed and applied to the virtual array to solve DOA estimation. Based on the proposed Bayesian learning framework, off-grid DOA can be estimated and meanwhile the impulsive noise can be effectively suppressed. Finally, comprehensive simulations are performed to verify the superiority of the proposed DOA estimation method. Furthermore, the advantages of the proposed method can be given as follows:

•  The proposed DOA method could accurately estimate DOA for MIMO radar under non-Gaussian impulsive noise without the prior knowledge of its parameters, such as the characteristic exponent of alpha stable distribution. Hence, the proposed method is with wider application scope than conventional DOA methods;
•  The proposed DOA method can estimate DOA more accurately than the conventional DOA methods, including the methods in Refs. [17,18,31], as an off-grid mechanism is adopted to design the proposed DOA method;
•  The proposed DOA method has a more robust performance than the conventional DOA estimation methods under impulsive noise, especially in the low generalized signal to noise ratio condition.

The framework of this paper is organized as follows: Section 2 describes the problem formulation. Section 3 explains the derivation of our method. Numerical experiments are presented in Section 4 and conclusions are drawn in Section 5.

Notation: The operator $\mathrm{diag}\left\{\mathbf{a}\right\}$ denotes a diagonal matrix using the elements in $\mathbf{a}$; ⊙ is the Khatri–Rao product; The operator $\mathrm{Tr}\left(\mathbf{A}\right)$ calculates the trace of $\mathbf{A}$; Matrices ${\mathbf{0}}_M$ and ${\mathbf{I}}_M$ denote the $M\times M$ matrix of zeros and identity, respectively.

2. Problem Formulation

2.1. Signal Model

Considering a monostatic MIMO radar with M transmitted and N received sensors, both of the transmitted and received arrays are assumed to be ULA with a half-wavelength space between the adjacent sensors. K far-field narrow band signals are assumed to be the signal of interest. Since the located arrays are close in monostatic MIMO radar, the directions of targets on the transmitted and received arrays are assumed to be same with ${\theta }_i(i=1,2,\dots ,K)$. Therefore, the output of the array with $MN\times 1$ sensors can be expressed as:

$$\mathbf{x}\left(t\right)=\mathbf{A}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right),$$

(1)

where $\mathbf{x}\left(t\right)={[{\mathbf{x}}_1\left(t\right),{\mathbf{x}}_2\left(t\right),\dots ,{\mathbf{x}}_{MN}\left(t\right)]}^{\mathsf{T}}$, $\mathbf{n}\left(t\right)\in {\mathbb{C}}^{MN\times 1}$ denotes the impulsive noise, $\mathbf{s}\left(t\right)\in {\mathbb{C}}^{K\times 1}$ is the source signal, $\mathbf{A}$ is the steering matrix whose dimension is $MN\times K$, and $\mathbf{A}={\mathbf{A}}_t\odot {\mathbf{A}}_r$ with ${\mathbf{A}}_t=[{\mathbf{a}}_t\left({\theta }_1\right),{\mathbf{a}}_t\left({\theta }_2\right),\dots ,{\mathbf{a}}_t\left({\theta }_K\right)]$, ${\mathbf{A}}_r=[{\mathbf{a}}_r\left({\theta }_1\right),{\mathbf{a}}_r\left({\theta }_2\right),\dots ,{\mathbf{a}}_r\left({\theta }_K\right)]$ and

$${\mathbf{a}}_t\left({\theta }_i\right)={[1,{\mathrm{e}}^{\mathrm{j}\pi sin\left({\theta }_i\right)},\dots ,{\mathrm{e}}^{\mathrm{j}\pi (M-1)sin\left({\theta }_i\right)}]}^{\mathsf{T}},$$

(2)

$${\mathbf{a}}_r\left({\theta }_i\right)={[1,{\mathrm{e}}^{\mathrm{j}\pi sin\left({\theta }_i\right)},\dots ,{\mathrm{e}}^{\mathrm{j}\pi (N-1)sin\left({\theta }_i\right)}]}^{\mathsf{T}},$$

(3)

where the operator ⊙ is the Khatri–Rao product, and for matrices ${\mathbf{F}}_{p\times n}$ and ${\mathbf{G}}_{q\times n}$, their Khatri–Rao product is defined as $\mathbf{F}\odot \mathbf{G}=[f_1\otimes g_1,f_2\otimes g_2,···,f_n\otimes g_n]\in R^{pq\times n}$ where ⊗ denotes the Kronecker product which is a kind of common operation in matrix analysis.

Collecting L snapshots of the received signals, the output can be represented as:

$$\mathbf{X}=\mathbf{A}\mathbf{S}+\mathbf{N},$$

(4)

where

$$\mathbf{X}=\left[\mathbf{x}\right(1),\mathbf{x}(2),\dots ,\mathbf{x}(L\left)\right],$$

(5)

$$\mathbf{S}=\left[\mathbf{s}\right(1),\mathbf{s}(2),\dots ,\mathbf{s}(L\left)\right],$$

(6)

$$\mathbf{N}=\left[\mathbf{n}\right(1),\mathbf{n}(2),\dots ,\mathbf{n}(L\left)\right].$$

(7)

For a non-circular signal, it can be represented by the real-valued signal ${\mathbf{s}}_0$ and the phase $\mathrm{\Phi }$, given as:

$$\mathbf{s}=\mathrm{\Phi }{\mathbf{s}}_0,$$

(8)

$${\mathbf{s}}_0={\left[s_{01},\dots ,s_{0K}\right]}^{\mathsf{T}},$$

(9)

$$\mathrm{\Phi }=\mathrm{diag}\{{\mathrm{e}}^{\mathrm{j}{\phi }_1/2},{\mathrm{e}}^{\mathrm{j}{\phi }_2/2},\dots ,{\mathrm{e}}^{\mathrm{j}{\phi }_K/2}\},$$

(10)

where $\mathrm{\Phi }$ is the non-circular phase of the input signal.

2.2. Impulsive Noise

In practical wireless communication environments, such as the impact of natural or man-made interference sources, Gaussian distribution is usually not suitable to describe the wireless channel noise. Correspondingly, alpha-stable distribution is usually adopted to describe the features of non-Gaussian impulsive noise [3,27]. As there is no closed probability density function of alpha-stable distribution, the characteristic function is introduced to evaluate the impulsiveness of noise as:

$$\phi \left(t\right)=exp\left\{\mathrm{j}\mu t-\gamma {\left|t\right|}^{\alpha }\left[1-\mathrm{j}\beta \mathrm{sgn}\left(t\right)\varpi \left(t,\alpha \right)\right]\right\},$$

(11)

where

$$\varpi \left(t,\alpha \right)=\left\{\begin{array}{c}tan\left(\pi \alpha /\phantom{\pi \alpha 2}2\right)\alpha \ne 1\hfill \\ \left(2/\phantom{2\pi }\pi \right)log\left|t\right|\alpha =1\hfill \end{array}\right.,$$

(12)

$$\mathrm{sgn}\left(t\right)=\left\{\begin{array}{c}1,t>0\hfill \\ 0,t=0\hfill \\ -1,t<0.\hfill \end{array}\right.,$$

(13)

$\alpha \in (0,2]$ denotes the characteristic exponent, and it mainly controls the thickness of distribution tails. The impulsive noise will produce more outliers with smaller $\alpha$. $\gamma$ represents the dispersion parameter with value ranges of $\gamma >0$. $\mu$ represents the location parameter. $\beta$ stands for the symmetry parameter taking values $-1<\beta \le 1$. When $\beta =0$, this distribution would become a symmetric alpha stable ($S\alpha S$) distribution.

In addition, when $\alpha <2$, the alpha stable distribution exhibits heavier tails and this distribution no longer has the finite SOS or HOS. Hence, the performance of conventional DOA estimation methods which are designed on the basis of SOS or HOS degrade seriously in impulsive noise.

2.3. Sparse Representation for MIMO DOA Estimation

Commonly, the dimension reduction is employed in MIMO array output to remove redundant information [11]. The transformation matrix is used to reduce the dimension, given as:

$$\mathbf{U}={\mathbf{F}}^{-1}{\mathbf{H}}^{\mathrm{H}},$$

(14)

where

$$\mathbf{F}=diag([1,2,\dots ,\underset{|M-N|+1}{\underbrace{min(M,N),\dots ,min(M,N)}},\dots ,2,1]),$$

(15)

$$\mathbf{H}={[{\mathbf{h}}_0,{\mathbf{h}}_1,,\dots ,{\mathbf{h}}_M]}^{\mathsf{T}},$$

(16)

and the $m-$th matrix ${\mathbf{h}}_m$ in $\mathbf{H}$ is defined as:

$${\mathbf{h}}_m={[{\mathbf{0}}_{N\times m},{\mathbf{I}}_{N\times N},{\mathbf{0}}_{N\times (M-m-1)}]}^{\mathsf{T}},$$

(17)

where ${\mathbf{I}}_{N\times N}$ is the unit matrix with the dimension of $N\times N$, and ${\mathbf{0}}_{N\times m}$ and ${\mathbf{0}}_{N\times (M-m-1)}$ denote the zero matrix with the dimension of $N\times m$ and $N\times (M-m-1)$, respectively.

Multiplying the output $\mathbf{X}$ with the transformation matrix $\mathbf{U}$, a virtual array output with dimension reduction can be obtained by:

$$\tilde{\mathbf{X}}=\mathbf{U}\mathbf{X}={\mathbf{F}}^{1/2}\mathbf{B}\mathbf{S}+\mathbf{U}\mathbf{N}=\overline{\mathbf{B}}\mathbf{S}+\overline{\mathbf{N}},$$

(18)

where $\tilde{\mathbf{X}}\in {\mathbb{C}}^{(M+N-1)\times L}$, $\overline{\mathbf{N}}=\mathbf{U}\mathbf{N}\in {\mathbb{C}}^{(M+N-1)\times L}$, $\overline{\mathbf{B}}={\mathbf{F}}^{1/2}\mathbf{B}\in {\mathbb{C}}^{(M+N-1)\times K}$, $\mathbf{B}=[\mathbf{b}\left({\theta }_1\right),\mathbf{b}\left({\theta }_2\right),\dots ,\mathbf{b}\left({\theta }_K\right)]$, ${\mathbf{F}}^{1/2}$ denotes the square root of $\mathbf{F}$ and

$$\mathbf{b}\left({\theta }_i\right)={[1,{\mathrm{e}}^{\mathrm{j}\pi sin\left({\theta }_i\right)},\dots ,{\mathrm{e}}^{\mathrm{j}\pi (M+N-2)sin\left({\theta }_i\right)}]}^{\mathsf{T}}.$$

(19)

Then, apply sparse representation method to $\tilde{\mathbf{X}}$, given as:

$$\tilde{\mathbf{X}}=\overline{\mathbf{B}}(\mathrm{\Theta }){\mathbf{S}}_v+\overline{\mathbf{N}},$$

(20)

where $\overline{\mathbf{B}}(\mathrm{\Theta })={\mathbf{F}}^{1/2}\mathbf{B}(\mathrm{\Theta })={\mathbf{F}}^{1/2}[\mathbf{b}\left({\theta }_1\right),\mathbf{b}\left({\theta }_2\right),\dots ,\mathbf{b}\left({\theta }_P\right)]$ and ${\left\{{\theta }_i\right\}}_{i=1}^P,P\gg K$ denote all possible DOAs in the sampling grid, ${\mathbf{S}}_v\in {\mathbb{C}}^{P\times L}$ denotes the sparse representation of $\mathbf{S}$. Only when ${\theta }_i$ equals to the true DOA of received signal, the ith row of ${\mathbf{S}}_v$ has a non-zero value. Therefore, DOA estimation is transformed to a sparse representation problem and DOA can be estimated by finding the non-zero elements in ${\mathbf{S}}_v$ after sparse reconstruction. Moreover, there are also some methods which construct the sparse representation model by exploiting the covariance $\mathbf{R}=\mathrm{E}\left[\tilde{\mathbf{X}}{\tilde{\mathbf{X}}}^{\mathsf{H}}\right]$ based vector. However, the above two models still have some limitations. Specifically, the sparse model of the original output may lead to large computation. The covariance-based model is not useful for coherent or correlated signals and it is also sensitive to the snapshots of signals.

3. Proposed Off-Grid DOA Estimation Method

In this section, we propose a novel off-grid DOA estimation method using sparse Bayesian learning for MIMO radar in impulsive noise. The details are presented as follows.

3.1. Sparse Signal Model

To construct a model with low computational complexity and meanwhile ensure the model is also suitable for coherent signals, a real-valued transformation is firstly utilized as:

$${\mathbf{X}}_R=\frac{1}{2}(\tilde{\mathbf{X}}+{\tilde{\mathbf{X}}}^{*}),$$

(21)

$${\mathbf{X}}_I=\frac{1}{2\mathrm{j}}(\tilde{\mathbf{X}}-{\tilde{\mathbf{X}}}^{*}).$$

(22)

Based on the non-circular signals property, Equation (21) and Equation (22) can be expressed as:

$${\mathbf{X}}_R=\frac{1}{2}(\overline{\mathbf{B}}\mathrm{\Phi }+{\overline{\mathbf{B}}}^{*}{\mathrm{\Phi }}^{*}){\mathbf{S}}_0+Re\left[\overline{\mathbf{N}}\right],$$

(23)

$${\mathbf{X}}_I=\frac{1}{2\mathrm{j}}(\overline{\mathbf{B}}\mathrm{\Phi }-{\overline{\mathbf{B}}}^{*}{\mathrm{\Phi }}^{*}){\mathbf{S}}_0+Im\left[\overline{\mathbf{N}}\right],$$

(24)

where the operators $\mathrm{Re}[·]$ and $\mathrm{Im}[·]$ extract the real and imaginary part of the variable, respectively. After that, we can derive a extended output model as:

$$\overline{\mathbf{X}}=\left[\begin{array}{c}{\mathbf{X}}_R\\ {\mathbf{X}}_I\end{array}\right]={\mathbf{A}}_e{\mathbf{S}}_0+\left[\begin{array}{c}Re\left[\mathbf{N}\right]\\ Im\left[\mathbf{N}\right]\end{array}\right],$$

(25)

where ${\mathbf{A}}_e$ denotes the new steering matrix with

$${\mathbf{A}}_e=\left[{\mathbf{a}}_e\left({\theta }_1\right),\dots ,{\mathbf{a}}_e\left({\theta }_K\right)\right],$$

(26)

$${\mathbf{a}}_e\left({\theta }_i\right)=\left[\begin{array}{c}{\mathbf{F}}_1{\mathbf{a}}_R\left({\theta }_i\right)\\ {\mathbf{F}}_1{\mathbf{a}}_I\left({\theta }_i\right)\end{array}\right],$$

(27)

$${\mathbf{a}}_R\left({\theta }_k\right)={[1,cos\left({\tau }_k\right),\dots ,cos\left((M+N-1){\tau }_k\right)]}^{\mathsf{T}},$$

(28)

$${\mathbf{a}}_I\left({\theta }_k\right)={[0,sin\left({\tau }_k\right),\dots ,sin\left((M+N-1){\tau }_k\right)]}^{\mathsf{T}},$$

(29)

where ${\tau }_k=\pi sin\left({\theta }_k\right)$, ${\mathbf{F}}_1={\mathbf{F}}^{1/2}$. When the received signals are interfered by impulsive noise, the output of signals can be further denoted by:

$$\overline{\mathbf{X}}={\mathbf{A}}_e{\mathbf{S}}_0+\mathbf{Z}+\mathbf{W}=\left[{\mathbf{A}}_e{\mathbf{I}}_{2Q}\right]\left[\begin{array}{c}{\mathbf{S}}_0\\ \mathbf{Z}\end{array}\right]+\mathbf{W},$$

(30)

where $Q=M+N-1$, $\mathbf{Z}$ represents a matrix which consists of outliers and $\mathbf{W}$ is a Gaussian distributed noise. Let $\overline{\mathbf{A}}=\left[{\mathbf{A}}_e{\mathbf{I}}_{2Q}\right]$ and $\overline{\mathbf{S}}={\left[{\mathbf{S}}_0\mathbf{Z}\right]}^{\mathsf{T}}$, the output can be further represented as:

$$\overline{\mathbf{X}}=\overline{\mathbf{A}}\overline{\mathbf{S}}+\mathbf{W}.$$

(31)

To construct a sparse model of Equation (31), we firstly consider the off-grid effect. Given a fixed sampling grid $\tilde{\theta }=[{\tilde{\theta }}_1,\dots ,{\tilde{\theta }}_P]$ from $[-0.5\pi ,0.5\pi ]$, $P\gg K$. Assuming ${\tilde{\theta }}_{nk}$ is the nearest grid point to the true DOA ${\theta }_k$, then the steering vector of ${\mathbf{a}}_e\left({\theta }_k\right)$ can be approximated by:

$${\mathbf{a}}_e\left({\theta }_k\right)\approx {\mathbf{a}}_e\left({\tilde{\theta }}_{nk}\right)+{\mathbf{b}}_e\left({\tilde{\theta }}_{nk}\right)({\theta }_k-{\tilde{\theta }}_{nk}),$$

(32)

where ${\mathbf{b}}_e\left({\tilde{\theta }}_{nk}\right)={\mathbf{a}}_e^{{}^{\prime }}\left({\tilde{\theta }}_{nk}\right)$. Furthermore, the off-grid sparse model of Equation (31) can be expressed as:

$$\overline{\mathbf{X}}={\overline{\mathbf{A}}}_v{\overline{\mathbf{S}}}_v+\mathbf{W},$$

(33)

where ${\overline{\mathbf{S}}}_v=[{\mathbf{S}}_1,\mathbf{Z}]$ and ${\mathbf{S}}_1$ is the sparse representation of ${\mathbf{S}}_0$, ${\overline{\mathbf{A}}}_v=[\overline{\mathbf{A}}(\mathrm{\Theta }),{\mathbf{I}}_{2Q}]+{\mathbf{B}}_1\mathrm{diag}\left({\mathbf{\beta }}_1\right)$, $\overline{\mathbf{A}}(\mathrm{\Theta })=[{\mathbf{a}}_e\left({\theta }_1\right),\dots ,{\mathbf{a}}_e\left({\theta }_P\right)]$, ${\mathbf{B}}_1=[{\mathbf{B}}_v,{\mathbf{0}}_{2Q}]$, ${\mathbf{\beta }}_1=[{\beta }_1,{\beta }_2,\dots ,{\beta }_P,{\mathbf{0}}_{1\times 2Q}]$ and

$${\mathbf{B}}_v=[{\mathbf{b}}_e\left({\tilde{\theta }}_1\right),\dots ,{\mathbf{b}}_e\left({\tilde{\theta }}_P\right)].$$

(34)

Although the array dimension has increased in this model, real-valued transformation can effectively reduce the computational complexity. This model is also suitable for coherent and correlated signals with high efficiency.

3.2. Sparse Bayesian Learning Based Approach

To apply Bayesian inference, we further analyze the distribution of different variables. For the Gaussian distributed noise $\mathbf{W}$, the distribution is assumed as:

$$p\left(\mathbf{W}\right|\epsilon )=\prod _{l=0}^{L-1}\mathcal{CN}\left(\mathbf{W}\left(l\right)\right|{\mathbf{0}}_{2Q\times 1},{\epsilon }^{-1}{\mathbf{I}}_{2Q}),$$

(35)

where $\mathbf{W}\left(l\right)$ is the lth column of $\mathbf{W}$, $\epsilon \simeq {\sigma }^{-2}$ and ${\sigma }^2$ represents the variance of noise. $\epsilon$ is a hyper prior parameter which is assumed with Gamma distribution,

$$p\left(\epsilon \right)=\mathsf{\Gamma }\left(\epsilon \right|1,a).$$

(36)

Then, the distribution of the output can be obtained by:

$$p\left(\overline{\mathbf{X}}\right|{\overline{\mathbf{S}}}_v,\epsilon ,{\mathbf{\beta }}_1)=\prod _{l=0}^{L-1}\mathcal{CN}\left({\overline{\mathbf{X}}}_l\right|{\overline{\mathbf{A}}}_v{\overline{\mathbf{S}}}_v\left(l\right),{\epsilon }^{-1}{\mathbf{I}}_{2Q}).$$

(37)

Since ${\mathbf{S}}_1$ is the sparse representation of ${\mathbf{S}}_0$, it can be assumed with Gaussian distribution as:

$$p\left({\mathbf{S}}_1\right|\mathbf{\delta })=\prod _{l=0}^{L-1}\mathcal{CN}\left({\mathbf{S}}_1\left(l\right)\right|{\mathbf{0}}_{P\times 1},\mathrm{diag}{\left(\mathbf{\delta }\right)}^{-1}),$$

(38)

where $\mathbf{\delta }=[{\delta }_1,{\delta }_2,\dots ,{\delta }_P]$ and ${\delta }_i$ is the variance’s inversion of ith row in ${\mathbf{S}}_1$. Because $\mathbf{\delta }$ is still a hyper parameter, Gamma distribution is also considered for $\mathbf{\delta }$

$$p\left(\mathbf{\delta }\right)=\prod _{n=1}^P\mathsf{\Gamma }\left({\delta }_n\right|1,\xi ),$$

(39)

where $\xi$ is a small value near to zero. Then, an IID (Independent and identically distributed) Gaussian distribution is utilized to model the outlier matrix $\mathbf{Z}$ as:

$$p\left(\mathbf{Z}\right|\mathrm{\Lambda })=\prod _{i=0}^{2Q-1}\prod _{l=0}^{L-1}\mathcal{CN}\left(Z_{i,l}\right|0,{\lambda }_{i,l}^{-1}),$$

(40)

where $Z_{i,l}$ and ${\lambda }_{i,l}$ are the $(i,l)$th element of $\mathbf{Z}$ and $\mathrm{\Lambda }$. $\mathrm{\Lambda }$ is a variance matrix to model the outlier matrix. Gamma distribution is also applied to model the hyper prior of $\mathrm{\Lambda }$ as:

$$p(\mathrm{\Lambda })=\prod _{i=0}^{2Q-1}\prod _{l=0}^{L-1}\mathsf{\Gamma }\left({\lambda }_{i,l}\right|1,\xi ).$$

(41)

As $\mathbf{Z}$ and ${\mathbf{S}}_1$ are considered independent, the distribution of ${\overline{\mathbf{S}}}_v$ could be calculated by:

$$p\left({\overline{\mathbf{S}}}_v\right|\mathrm{\Omega })=\prod _{l=0}^{L-1}\mathcal{CN}\left({\overline{\mathbf{S}}}_v\left(l\right)\right|{\mathbf{0}}_{(N_0+2Q)\times 1},{\mathrm{\Omega }}_l^{-1}),$$

(42)

where ${\mathrm{\Omega }}_l=\mathrm{diag}\left([\mathbf{\delta },{\mathrm{\Lambda }}_l]\right)$. Furthermore, the distribution of ${\mathbf{\beta }}_1$ is modeled by uniform distribution as:

$${\mathbf{\beta }}_1\sim \mathcal{U}\left([-0.5r,0.5r]\right),$$

(43)

where r is the sampling interval of the original coarse grid.

On the basis of these distributions, we can realize sparse reconstruction via maximizing following posterior probability:

$${\mathbf{S}}_{out}=arg\underset{{\overline{\mathbf{S}}}_v}{max}p\left({\overline{\mathbf{S}}}_v\right|\overline{\mathbf{X}},\epsilon ,\mathbf{\delta },\mathrm{\Lambda },{\mathbf{\beta }}_1).$$

(44)

Based on the above distributions, we can get:

$$p\left({\overline{\mathbf{S}}}_v\right|\overline{\mathbf{X}},\epsilon ,\mathbf{\delta },\mathrm{\Lambda },{\mathbf{\beta }}_1)\triangleq \prod _{l=0}^{L-1}\mathcal{CN}\left({\overline{\mathbf{S}}}_v\left(l\right)\right|{\mathbf{\mu }}_l,{\mathrm{\Sigma }}_l),$$

(45)

where

$${\mathbf{\mu }}_l=\epsilon {\mathrm{\Sigma }}_l{\overline{\mathbf{A}}}_v^{\mathsf{T}}\overline{\mathbf{X}}\left(l\right),$$

(46)

$${\mathrm{\Sigma }}_l={(\epsilon {\overline{\mathbf{A}}}_v^{\mathsf{T}}{\overline{\mathbf{A}}}_v+{\mathrm{\Omega }}_l)}^{-1}.$$

(47)

With the aid of Equations (46) and (47), sparse reconstruction can be realized. However, Equations (46) and (47) cannot be directly calculated without estimating the unknown parameters. To solve this problem, expectation maximization (EM) approach is introduced to calculate the unknown parameters. Considering the expectation likelihood function with respect to the posterior of ${\overline{\mathbf{S}}}_v$,

$$\begin{array}{cc}\hfill \mathcal{L}(\epsilon ,\mathbf{\delta },\mathrm{\Lambda },& {\mathbf{\beta }}_1)=\mathbb{E}[lnp(\overline{\mathbf{X}},{\overline{\mathbf{S}}}_v,\epsilon ,\mathbf{\delta },\mathrm{\Lambda },{\mathbf{\beta }}_1)]\hfill \\ \hfill & =\mathbb{E}[lnp\left(\overline{\mathbf{X}}\right|{\overline{\mathbf{S}}}_v,\epsilon )p\left({\overline{\mathbf{S}}}_v\right|\mathbf{\delta },\mathrm{\Lambda })p\left(\epsilon \right)p\left(\mathbf{\delta }\right)p(\mathrm{\Lambda })p\left({\mathbf{\beta }}_1\right)],\hfill \end{array}$$

(48)

where $\mathbb{E}$ denotes an expectation to the posterior of ${\overline{\mathbf{S}}}_v$. On the basis of this likelihood function, the hyper prior parameters can be further calculated as follows.

For $\epsilon$, ignoring the independent terms, the following can be obtained:

$$\begin{array}{cc}\hfill & \mathcal{L}\left(\epsilon \right)\propto \mathbb{E}[lnp\left(\overline{\mathbf{X}}\right|{\overline{\mathbf{S}}}_v,\epsilon ,\mathbf{\delta },\mathrm{\Lambda },{\mathbf{\beta }}_1)p\left(\epsilon \right)]\hfill \\ \hfill & =\mathbb{E}\left[\prod _{l=0}^{L-1}\mathcal{CN}\left(\overline{\mathbf{X}}\left(l\right)\right|{\overline{\mathbf{A}}}_v{\overline{\mathbf{S}}}_v\left(l\right),{\epsilon }^{-1}{\mathbf{I}}_{2Q})\right]+lnp\left(\epsilon \right)\hfill \\ \hfill & =LQln\epsilon -\epsilon \sum _{l=0}^{L-1}\mathbb{E}\left[\parallel \overline{\mathbf{X}}\left(l\right)-{\overline{\mathbf{A}}}_v{\overline{\mathbf{S}}}_v{\left(l\right)\parallel }^2\right]+lnp\left(\epsilon \right)\hfill \\ \hfill & =LQln\epsilon -\epsilon \sum _{l=0}^{L-1}{\parallel \overline{\mathbf{X}}\left(l\right)-{\overline{\mathbf{A}}}_v{\mathbf{\mu }}_l\parallel }^2-\epsilon \sum _{l=0}^{L-1}Tr\left({\overline{\mathbf{A}}}_v^{\mathsf{T}}{\overline{\mathbf{A}}}_v{\mathrm{\Sigma }}_l\right)-a\epsilon .\hfill \end{array}$$

(49)

Let $\frac{\partial \mathcal{L}\left(\epsilon \right)}{\partial \epsilon }=0$; $\epsilon$ will be obtained by

$$\epsilon =\frac{LQ}{a+{\sum }_{l=0}^{L-1}{\parallel \overline{\mathbf{X}}\left(l\right)-{\overline{\mathbf{A}}}_v{\mathbf{\mu }}_l\parallel }^2+Tr\left({\overline{\mathbf{A}}}_v^{\mathsf{T}}{\overline{\mathbf{A}}}_v{\mathrm{\Sigma }}_l\right)}.$$

(50)

For $\mathbf{\delta }$ and $\mathrm{\Lambda }$, the following can be found:

$$\begin{array}{cc}\hfill & \mathcal{L}(\mathbf{\delta },\mathrm{\Lambda })\propto \mathbb{E}[lnp\left({\overline{\mathbf{S}}}_v\right|\mathbf{\delta },\mathrm{\Lambda })p\left(\mathbf{\delta }\right)p(\mathrm{\Lambda })]\hfill \\ \hfill & =\mathbb{E}\left[\prod _{l=0}^{L-1}\mathcal{CN}\left({\overline{\mathbf{S}}}_v\left(l\right)\right|{\mathbf{0}}_{(P+2Q)\times 1},{\mathrm{\Omega }}_l^{-1})\right]+lnp\left(\mathbf{\delta }\right)+lnp(\mathrm{\Lambda })\hfill \\ \hfill & =L\sum _{n=1}^{P}ln{\delta }_n+\sum _{i=0}^{2Q-1}\sum _{l=0}^{L-1}ln{\lambda }_{i,l}-\xi \sum _{n=1}^P{\delta }_n-\xi \sum _{i=0}^{2Q-1}\sum _{l=0}^{L-1}{\lambda }_{i,l}-\sum _{l=0}^{L-1}Tr\left((\mathbf{\mu }{\mathbf{\mu }}_l^{\mathsf{T}}+{\mathrm{\Sigma }}_l){\mathrm{\Omega }}_l\right).\hfill \end{array}$$

(51)

Similarly, differentiate Equation (51) to ${\delta }_n$ and ${\lambda }_{i,l}$, respectively. They will be updated by:

$${\delta }_n=\frac{L}{\xi +{\sum }_{l=0}^{L-1}{\left[{\mathbf{D}}_l\right]}_{n,n}},$$

(52)

$${\lambda }_{i,l}=\frac{1}{\xi +{\left[{\mathbf{D}}_l\right]}_{i+P,i+P}},$$

(53)

where ${\mathbf{D}}_l=Tr(\mathbf{\mu }{\mathbf{\mu }}_l^{\mathsf{T}}+{\mathrm{\Sigma }}_l)$.

For ${\mathbf{\beta }}_1$, it can be updated through ${\mathbf{\beta }}_1={\mathbf{V}}^{-1}\mathbf{\psi }$ if $\mathbf{V}$ is invertible, given as:

$$\mathbf{V}=\left({\overline{\mathbf{B}}}_1^{\mathsf{T}}{\overline{\mathbf{B}}}_1\right)\odot \left(\sum _{l=0}^{L-1}({\mathbf{\mu }}_l{\mathbf{\mu }}_l^{\mathsf{T}}+{\mathrm{\Sigma }}_l)\right),$$

(54)

$$\mathbf{\psi }=\frac{1}{L}\sum _{l=0}^{L-1}\left(\mathrm{diag}\left({\mathbf{\mu }}_l\right){\mathbf{B}}_1^{\mathsf{T}}(\overline{\mathbf{X}}\left(l\right)-{\overline{\mathbf{A}}}_v{\mathbf{\mu }}_l)-\mathrm{diag}\left({\mathbf{B}}_1^{\mathsf{T}}{\overline{\mathbf{A}}}_v{\mathrm{\Sigma }}_l\right)\right).$$

(55)

Additionally, if $\mathbf{V}$ is not invertible, ${\mathbf{\beta }}_1$ can be updated by:

$${\widehat{\mathbf{\beta }}}_1\left(n\right)=\frac{\mathbf{\psi }\left(n\right)-{\mathbf{V}}_{-n}^{\mathsf{T}}{\mathbf{\beta }}_{1,-n}}{{\mathbf{V}}_{n,n}},$$

(56)

$${\mathbf{\beta }}_1{\left(n\right)}^{new}=\left\{\begin{array}{cc}\hfill {\widehat{\mathbf{\beta }}}_1\left(n\right),& if{\widehat{\mathbf{\beta }}}_1\left(n\right)\in [-0.5r,0.5r]\hfill \\ \hfill -0.5r,& if{\widehat{\mathbf{\beta }}}_1\left(n\right)<-0.5r\hfill \\ \hfill 0.5r,& \mathrm{otherwise}\hfill \end{array}\right.,$$

(57)

where ${\mathbf{\beta }}_1\left(n\right),\mathbf{\psi }\left(n\right)$ are the nth elements of ${\mathbf{\beta }}_1$, and $\mathbf{\psi }$, respectively. ${\mathbf{V}}_{-n}$ and ${\mathbf{\beta }}_{1,-n}$ denote the operation which removes the nth entry of $\mathbf{V}$ and ${\mathbf{\beta }}_1$.

Now, sparse reconstruction based SBL is realized by repeating the above procedure until it meets the terminal conditions. The update procedure is given in Algorithm 1. After sparse reconstruction, DOAs will be calculated based on the indexes of non-zero elements in ${\mathbf{S}}_{out}$.

Algorithm 1 Proposed Method

Input:  Input signals $\mathbf{X}$ and sampling grid steering matrix $\overline{\mathbf{A}}(\mathrm{\Theta })$ and its derivation ${\mathbf{B}}_v$. Output:  ${\mathbf{S}}_{out}$ 1: Construct the virtual array signals $\tilde{\mathbf{X}}$ and transform it to a real-valued model $\overline{\mathbf{X}}$. 2: Set the iteration number $t=0$ and initialize the value of different parameters as $a=0.01$, $\xi =0.001$, ${\mathbf{\delta }}^{\left(0\right)}={\mathbf{I}}_{P\times 1}$, ${\mathrm{\Lambda }}^{\left(0\right)}={\mathbf{I}}_{2Q\times L}$, ${\mathbf{\beta }}_1^{\left(0\right)}={\mathbf{0}}_{(P+2Q)\times 1}$. 3: Construct the off-grid steering matrix ${\overline{\mathbf{A}}}_v=\overline{\mathbf{A}}(\mathrm{\Theta })+{\mathbf{B}}_1\mathrm{diag}\left({\mathbf{\beta }}_1^{\left(t\right)}\right)$. 4: Let $t=t+1$ and calculate ${\mathbf{\mu }}_l^{\left(t\right)}$ and ${\mathrm{\Sigma }}_l^{\left(t\right)}$ for $l=1,2,\dots ,L$ by Equations (46) and (47). 5: Calculate ${\mathbf{\delta }}^{\left(t\right)}$ and ${\mathrm{\Lambda }}^{\left(t\right)}$ by Equations (52) and (53). 6: Calculate $\epsilon$ by Equation (50) and update ${\mathbf{\beta }}_1^{\left(t\right)}$ by ${\mathbf{\beta }}_1^{\left(t\right)}={\mathbf{V}}^{-1}\mathbf{\psi }$ or Equation (57). 7: Repeat step 3-6 until $\parallel {\mathbf{\delta }}^{\left(t\right)}-{\mathbf{\delta }}^{(t-1)}{\parallel }_2/{\parallel {\mathbf{\delta }}^{(t-1)}\parallel }_2\le 0.001$ or $t>500$. 8: Output ${\mathbf{S}}_{out}=[{\mathbf{\mu }}_1,{\mathbf{\mu }}_2,\dots ,{\mathbf{\mu }}_L]$.

4. Simulation Results and Analysis

In this section, a series of simulations are provided to evaluate the proposed method. A monostatic MIMO radar with seven transmitted arrays and six received arrays is considered. The simulation dataset has been partly attached in the Supplementary Materials. Moreover, in this dataset, two far-field narrow band BPSK signals are utilized as transmitted non-circular signals. DOAs of signals are considered as ${\theta }_1=45.4^{\circ }$ and ${\theta }_2=31.8^{\circ }$. The interval of the fixed sampling grid is set as $r=1^{\circ }$. Impulsive noise is induced by symmetric alpha-stable distribution. As the variance is not convergent for alpha-stable distribution, a generalized signal to noise ratio (GSNR) is usually utilized to measure the intensity of impulsive noise, given as

$$\mathrm{GSNR}10lg(P_s/\gamma ),$$

(58)

where $P_s$ is the signal power and $\gamma$ is the dispersion parameter of alpha-stable distribution.

The comparison experiments are carried out in different parameters, mainly including characteristic exponent $\alpha$, GSNR and signal snapshots. We compare our method with the three latest methods, mainly including nuclear norm minimization (NUC-NORM) [17], iterative reweighted proximal projection (IRPP) [18], robust variational Bayesian inference (RVSBL) [31] and multiple signal classification (MUSIC). NUC-NORM and IRPP are proposed for MIMO DOA estimation and RVSBL is proposed to deal with impulsive noise via Bayesian learning. Meanwhile, two quantities, named accuracy and root mean square error (RMSE), are calculated with 200 Monte-Carlo experiments. An accurate estimation of DOA is defined as:

$$|{\widehat{\theta }}_k\left(i\right)-{\theta }_k|\le 2^{\circ },k=1,2,$$

(59)

where ${\theta }_k$ and ${\widehat{\theta }}_k\left(i\right)$ denote the true DOA of the kth signal and the estimated DOA in the ith simulation, respectively. Then the accuracy can be obtained utilizing the ratio of the number of accurate estimations to the total number of experiments. RMSE will be achieved by:

$$\mathrm{RMSE}=\sqrt{\frac{1}{N_{s}K}\sum _{i=1}^{N_s}\sum _{k=1}^K{({\widehat{\theta }}_k\left(i\right)-{\theta }_k)}^2},$$

(60)

where $N_s$ is the number of accurate estimations for DOAs in the experiments.

4.1. Simulations versus GSNRs

In this subsection, we conduct the experiments with GSNR in a range from 1 dB to 10 dB. The characteristic exponent of noise is set as $\alpha =1.5$, and the snapshot of signals is set as 10. Simulation results are illustrated in Figure 1 and Figure 2.

It is apparent from Figure 1 and Figure 2 that the proposed method outperforms the contrastive DOA methods under different GSNRs. The accuracy of the proposed method can converge to 1 when GSNR ≥ 3dB, while the contrastive methods require GSNR to be at least 8dB so as to achieve the same performance. Furthermore, RMSE for different methods is illustrated in Figure 2. It shows that our DOA method performs better than other contrastive methods. Among these methods, the proposed method can achieve the best performance in terms of both accuracy and RMSE since it can resist the influence of impulsive noise and off-grid effect. Thus, our method is more suitable for DOA estimation in different GSNRs.

4.2. Simulations versus Characteristic Exponents

In this subsection, we conduct the experiments with characteristic exponents $\alpha$ in a range from 1.1 to 2. GSNR remains unchanged at 5 dB. The snapshot of signals is set as 10. The results are demonstrated in Figure 3 and Figure 4.

As shown in Figure 3 and Figure 4, the performance of all methods become better when $\alpha$ increases. This is because the impulsiveness becomes weaker when $\alpha$ is larger. Clearly, compared with other contrastive methods, our method has higher accuracy in all characteristic exponent conditions. Moreover, we find only our method and RVSBL can keep good accuracies with small $\alpha$ since they have the resistance to impulsive noise. Meanwhile, our method shows a superior performance in terms of RMSE than other methods due to its effective consideration for the off-grid condition.

4.3. Simulations versus Signal Snapshot

With this part, we make the performance comparison with different signal snapshots. Characteristic exponent and GSNR are set as $\alpha$ = 1.5 and 5 dB, respectively. The snapshot ranges from 10 to 100. Simulation results are given in Figure 5 and Figure 6.

As shown in Figure 5 and Figure 6, it can be observed that only our method and RVSBL can keep high accuracies in all conditions because of their suppression of impulsive noise. In addition, our method can even keep a totally successful estimation in all experiments. By contrast, IRPP- and NUC-NORM-based methods show bad and non-stable accuracy in the simulations. Furthermore, our method shows the best performance on RMSE due to its operation for off-grid condition.

4.4. Simulations versus DOAs

In this subsection, we mainly conduct the experiments with different DOAs of signals. Different DOAs and intervals are considered, and ${\theta }_1,{\theta }_2$ are mainly set as $(31.8^{\circ },45.4^{\circ })$, $(0.5^{\circ },8.7^{\circ }),(84.2^{\circ },73.5^{\circ })$, $(0.8^{\circ },87.4^{\circ })$. The characteristic exponents $\alpha$ is set as 1.5, and the snapshot of signals is set as 10. The results are demonstrated in

Table 1. Accuracy with different DOAs of signals.

DOA Method	($31.8^{\circ },45.4^{\circ }$)	($0.5^{\circ },8.7^{\circ }$)	($84.2^{\circ },73.5^{\circ }$)	($0.8^{\circ },87.4^{\circ }$)
Proposed	100%	100%	98%	60%
RVSBL	93%	95%	88%	38%
IRPP	58%	50%	46%	3%
NUC-NORM	12%	10%	8%	2%

.

It is apparent from Table 1 that the proposed method outperforms the contrastive DOA methods under different DOAs. The accuracy of the proposed method can achieve a superior performance under different DOAs and different intervals for a fixed samples grid, which is higher than the other contrastive methods. Among these DOA methods, the proposed method has the highest accuracy in terms of different intervals and DOAs even though DOAs are close to 0 and 90 degrees. Thus, our method has an outstanding performance in different DOAs.

Since the proposed DOA method and RVSBL show a better performance than other methods, the computation complexity of these two methods is compared. According to Ref. [31], the computation complexity of RVSBL is $O\left(P^3\right)$. The computation complexity of the proposed method is $O\left({\left(P+2Q\right)}^3\right)$. Since $P>>Q$, the complexities of two methods are comparable. Otherwise, with the help of real-valued transformation in the proposed method, the complexity can be reduced by a scale of four with the same matrix dimension. Therefore, the proposed method can show a better performance with proper computation complexity.

5. Conclusions

In this paper, an off-grid DOA estimation method for monostatic MIMO radar is proposed under impulsive noise. Based on the property of monostatic MIMO radar and non-circular signals, a real valued virtual array is constructed, which can extend the array aperture and reduce the computational complexity for Bayesian learning. Then, considering the property of impulsive noise, an off-grid sparse model is built for the virtual array. With the help of SBL, DOA estimation is realized by reconstructing the off-grid sparse model. A series of simulations is provided to verify the superiority of the proposed method, and simulation results indicate that the proposed method can lead to high accuracy and robust DOA estimation with monostatic MIMO radar in impulsive noise.