Sparse Bayesian Learning Based Three-Dimensional Imaging Algorithm for Off-Grid Air Targets in MIMO Radar Array

Abstract

In recent years, the development of compressed sensing (CS) and array signal processing provides us with a broader perspective of 3D imaging. The CS-based imaging algorithms have a better performance than traditional methods. In addition, the sparse array can overcome the limitation of aperture size and number of antennas. Since the signal to be reconstructed is sparse for air targets, many CS-based imaging algorithms using a sparse array are proposed. However, most of those algorithms assume that the scatterers are exactly located at the pre-discretized grids, which will not hold in real scene. Aiming at finding an accurate solution to off-grid target imaging, we propose an off-grid 3D imaging method based on improved sparse Bayesian learning (SBL). Besides, the Bayesian Cramér-Rao Bound (BCRB) for off-grid bias estimator is provided. Different from previous algorithms, the proposed algorithm adopts a three-stage hierarchical sparse prior to introduce more degrees of freedom. Then variational expectation maximization method is applied to solve the sparse recovery problem through iteration, during each iteration joint sparsity is used to improve efficiency. Experimental results not only validate that the proposed method outperforms the existing off-grid imaging methods in terms of accuracy and resolution, but have compared the root mean square error with corresponding BCRB, proving effectiveness of the proposed method.

1. Introduction

In the last few years, three-dimensional (3D) radar imaging systems and algorithms have received significant attention among researchers worldwide. Different from the 2D radar image, which is a projection of the 3D model, 3D radar image can provide shape and energy distribution of dominated scatterers of the target [1,2]. This information is of great significance to target identification and diagnostic analysis. Multiple-input–multiple-output (MIMO) radar system has drawn much attention and shows its potential in 3D imaging in the last decade [3]. For one thing, it can overcome the problem of limited number of antennas and limited aperture size. For another, the coherent processing interval (CPI) is reduced in MIMO radar by using space sampling instead of time sampling [4]. In traditional interferometric inverse synthetic aperture radar (InISAR), the imaging performance is severely degraded for high maneuvering target because the rotation axis of the target relative to the radar is time varying during the CPI [5,6]. Compared with the InISAR technique, the 3D image can be reconstructed in one snapshot [7] with a wideband monostatic MIMO radar system, thus free from the target’s high mobility. In the last few years, compressed sensing (CS) has been widely used in sparse signal recovery with fewer samples [8]. It has shown great advantage in radar imaging due to its super-resolution ability. It is proven that CS-based imaging algorithms provide a better resolution enhancement effect than the RELAX algorithm [9]. Besides, it can also exploit the sparsity of signal and reconstruct the signal from limited samples with high probability. Thus, the combination of MIMO radar and CS technique becomes a hot topic, and much existing research focus on CS recovery techniques for MIMO radar imaging [4,10,11,12]. In Ding et al. [4], a CS-based 2D imaging method based on sparse array is proposed. In Hu et al. [12], a narrowband 3D imaging method based on Kronecker CS is discussed. As a matter of fact, the sparse recovery algorithms severely affect the imaging performances. Hence, sparse recovery algorithm with high accuracy is of fundamental importance in MIMO 3D imaging.

The CS-based imaging methods, such as orthogonal matching pursuit (OMP) or basis pursuit (BP), can achieve a better imaging performance in the cases of low signal-to-noise-ratio (SNR) and limited snapshots. However, most of the existing CS-based imaging methods require that the scatterers of target to be exactly on the discretized sampling grid [13]. Otherwise, the off-grid problem will introduce basis mismatch and leakage of energy over all grids, finally lead to the deterioration of imaging performance. In practical MIMO radar imaging, the targets are distributed in the continuous space, and the off-grid problem always exists. There are mainly two ways to solve the problem. One is to compensate for the mismatch, the other is to use the gridless sparse recovery [14]. The latter is more accurate but the guaranteed theoretical resolution of this kind of method is a few Rayleigh limits [15]. So in practical CS imaging, we still use the former method. A denser grid may mitigate the effect of mismatch to some extent, but it will increase the mutual coherence of the sensing matrix [16]. This may cause violation of the restricted isometry property (RIP) property, which is the guarantee for reliable recovery. The computation cost will also increase sharply with a denser grid. Hence, CS recovery with a perturbed sensing matrix has been a hot topic. In Chi et al. [17], the performance of CS recovery is analyzed when the mismatch problem exists, nevertheless, it does not provide any algorithm for off-grid CS recovery. Sparsity-cognizant total least-squares (S-TLS) is the first recovery method for perturbed compressed sensing problem [18]. The S-TLS algorithm supposes that the off-grid bias obeys the Gaussian distribution. However, without any prior information, the uniform distribution is more suitable for the off-grid bias (the distance between the true scatterer and its nearest grid). A perturbation approach is established for compressive radar imaging based on OMP [13], but it does not use the statistical information and the performance is not satisfying. An off-grid direction of arrival (DOA) estimator based on sparse Bayesian inference (OGSBI) is proposed in Yang et al. [19], which adopts the uniform distribution. However, research on MIMO radar CS imaging are mainly restricted to one dimensional(1D) resolution enhancement [12]. Study on problems with high dimensionality has just begun. Hence, this study sets out to find an algorithm for high dimensionality MIMO imaging for off-grid air target.

Starting from the purpose of finding a 3D imaging method with high accuracy for maneuvering off-grid air target, we propose a novel algorithm based on improved sparse Bayesian learning to overcome the aforementioned problems. Under this framework, we first derive the on-grid and off-grid 3D imaging model using a sparse array. Then the three-stage hierarchical sparse prior is adopted, which introduces more degrees of freedom. Finally, this algorithm can get the approximated analytic expressions of the unknown parameters by using the variational expectation maximization (EM) method. Besides, both the imaging result and estimations of unknown parameters such as off-grid biases can be obtained at the same time. Simulation results prove the superiority of the proposed algorithm over other existing off-grid sparse recovery algorithms such as S-TLS and OGSBI. In addition, the BCRB and the Mean Square Error (MSE) of off-grid bias estimator are compared to verify effectiveness of the proposed method.

Notations used in this paper are as follows. Bold-case letters are reserved for vectors and matrices. $\widehat{a}$ is the expectation of a. $diag(\mathit{x})$ is a matrix with its main diagonal being $\mathit{x}$. ⊙ denotes the Hardmard product and $\mathsf{\Re }\{\mathit{x}\}$ means the real part of complex $\mathit{x}$. ${\mathcal{K}}_p(·)$ is the modified Bessel function of the 3rd kind. ${(·)}^{\mathit{T}}$, ${(·)}^{\mathit{H}}$ are the transpose and conjugate transpose operator, respectively. $\overline{\mathit{A}}$ is the conjugate of $\mathit{A}$. $E{[·]}_p$ is the expectation operator which is taken with respect to the probability density function p. $x^{(i)}$ is the update of x in the ith iteration. ${\mathit{A}}_i,{\mathit{A}}^j,{\mathit{A}}_{i,j}$ are the ith column, the jth row and the $(i,j)$th element of $\mathit{A}$, respectively.

The remaining parts of this paper proceed as follows: the off-grid target 3D imaging problem based on sparse array is formulated in Section 2. In Section 3, an algorithm for 3D imaging based on improved sparse Bayesian learning is proposed and the BCRB for off-grid bias is also derived. In Section 4, experimental results and analyses of the imaging performance are presented. Finally, conclusions are drawn in Section 5.

2. Problem Formulation of 3D Imaging

2.1. Ideal Imaging Model Based on Sparse Antenna Array

The geometry of the 3D imaging based on sparse antenna array is depicted in Figure 1. In this paper, we consider a wideband MIMO radar system with M transmitters and N receivers. Considering the need for a good orthogonal characteristic of transmitted signals, we adopt a group of M orthogonal phase-code modulation (PCM) signals with the same carrier frequency and bandwidth in this system.

Transmitted signal from the mth transmitter can be expressed as follows:

$$s_m(t)=rect(\frac{t-T_r}{T_r})\exp (\mathrm{j}{\phi }_m(t))\exp (\mathrm{j}2\pi f_{c}t)(m=1,2,\dots ,M)$$

(1)

In which $T_r$ is the pulse width, $f_c$ is the carrier frequency and ${\phi }_m(t)$ is the phase-code function. $rect(\frac{t}{T})$ is defined as follows:

$$rect(\frac{t}{T})=\{\begin{array}{cc}0\hfill & \left|t\right|>\frac{T}{2}\hfill \\ 1\hfill & \left|t\right|\le \frac{T}{2}\hfill \end{array}$$

(2)

This group of transmitted signals are orthogonal to each other.

$${\int }_0^{T_r}{s}_p(t)s_q^{*}(t)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}p\ne q\hfill \\ Constant\hfill & \mathrm{if}p=q\hfill \end{array}\right.$$

(3)

Assuming that the transmitters and receivers are in the XOY plane, $(x_m,y_m)$ and $(x_n,y_n)$ are the coordinates of mth transmitter and nth receiver, respectively. Suppose that there are K dominant scatterers in the imaging scene. The coordinate of the kth scatterer is ($r_k,{\theta }_k,{\phi }_k$), and its RCS is ${\sigma }_k$. According to the sparse radar array configuration in Figure 1, the transmit and receive steering vector can be described:

$$\mathit{t}(k)=[1\dots \exp (\mathrm{j}\frac{2\pi }{\lambda }{\gamma }_t(m,k))\dots \exp (\mathrm{j}\frac{2\pi }{\lambda }{\gamma }_t(M,k))]$$

(4)

$$\mathit{r}(k)=[1\dots \exp (\mathrm{j}\frac{2\pi }{\lambda }{\gamma }_r(n,k))\dots \exp (\mathrm{j}\frac{2\pi }{\lambda }{\gamma }_r(N,k))]$$

(5)

$${\gamma }_t(m,k)=(x_m\sin {\theta }_k\cos {\phi }_k+y_m\sin {\theta }_k\sin {\phi }_k)$$

(6)

$${\gamma }_r(n,k)=(x_n\sin {\theta }_k\cos {\phi }_k+y_n\sin {\theta }_k\sin {\phi }_k)$$

(7)

In which $\lambda =\frac{c}{f_c}$ is the carrier wavelength, and c is the speed of light. Define the transmitting vector as $\mathit{s}(t)={[s_1(t),s_2(t),\dots ,s_M(t)]}^T$, then the received signal can be expressed as:

$$\mathit{g}(t)=\sum _{k=1}^K{\sigma }_k{\mathit{r}}^T(k)\mathit{t}(k)\mathbf{s}(t-\frac{2r_k}{c})$$

(8)

After the matched filters are adopted, we can get the received signal from the mth transmitter and the nth receiver

$$g_{mn}(t)=\sum _{k=1}^K{\sigma }_k{r}_n(k)t_m(k){\psi }_m(t-\frac{2r_k}{c})$$

(9)

In Equation (9), $r_n(k)$ is the nth element of $\mathit{r}(k)$, and $t_m(k)$ is similar. ${\psi }_m(t)$ is the autocorrelation function of $s_m(t)$, which is the 1D range image. Since we adopt a wideband sparse antenna array, the range resolution can be achieved through matched filtering. After above operations, we can get $MN$ one-dimensional range images from the N receivers. Because the range resolution provided by the bandwidth is usually sufficient, we only use CS technique to enhance resolutions of the other two dimensions. After range alignment, we transfer the 3D imaging problem into a corresponding lower dimensional one. For the ith range cell, assuming that there are I out of K dominant scatterers falling in this range cell, echo return of the ith range cell can be formulated as follows:

$$y_{mn}=\sum _{i=1}^I{\sigma }_i{r}_n(i)t_m(i)$$

(10)

For the application of air targets imaging, the true scatterers are always sparse compared with the background. Thus, the sparsity of those targets are exploited and CS-based methods show its advantages. According to traditional CS-based imaging method, we then discretize the imaging scene into uniform grids. Assuming there are P grids for $\theta$ and Q grids for $\phi$, i.e., $Z(=PQ)$ grids in total. Then the observed formula with noise can be described as:

$$y_{mn}=\sum _{p=1}^P\sum _{q=1}^Q{\sigma }_{p,q}\exp \{\mathrm{j}\frac{2\pi }{\lambda }[(x_m+x_n)\sin {\theta }_p\cos {\phi }_q+[y_m+y_n]\sin {\theta }_p\sin {\phi }_q]\}+e_{mn}$$

(11)

After data rearrangement and vectorization, i.e., ${\sigma }_i={\sigma }_{p,q},i=(q-1)\times P+p$, the signal model for sparse recovery is as follows:

$$y_{mn}=\sum _i^Z{\sigma }_i\exp \{\mathrm{j}\frac{2\pi }{\lambda }[(x_m+x_n)\sin {\theta }_i\cos {\phi }_i+[y_m+y_n]\sin {\theta }_i\sin {\phi }_i]\}+e_{mn}$$

(12)

$${\mathit{y}=\mathbf{\Phi }}_0\mathbf{\sigma }+\mathit{e}$$

(13)

In Equation (13) $\mathit{y}$ is a $MN\times 1$ vector standing for the echo return after range alignment. ${\mathbf{\Phi }}_0$ is a $MN\times Z$ matrix with its $(i,j)$th element being $\exp \{\mathrm{j}\frac{2\pi }{\lambda }[(x_m+x_n)\sin {\theta }_j\cos {\phi }_j+[y_m+y_n]\sin {\theta }_j\sin {\phi }_j]\}$, in which $i=(n-1)\times M+m$. $\mathbf{\sigma }$ for different range cells can be extracted from Equation (13) via sparse recovery, which finally construct the 3D image.

2.2. Off-Grid Imaging Model Using Taylor Expansion

The imaging model in Equation (12) is under the assumption that all the dominant scatterers are exactly located on the pre-discretized grid. In practical radar imaging this assumption does not hold since the scatterers of real target are distributed in a continuous space. This means that the off-grid problem always exists. There are two off-grid models to solve the problem, one is based on the Taylor series expansion and the other based on linear interpolation. The performances of these two models are compared in Das [20]. The result is that the former has a better accuracy, so we choose the Taylor series expansion to formulate the off-grid imaging model. Considering the mismatch between scatterer’s true position and the grids, the true position of a scatterer can be expressed as:

$$\{\begin{array}{c}{\theta }_{i0}={\theta }_i+{\delta }_{{\theta }_i}\hfill \\ {\phi }_{i0}={\phi }_i+{\delta }_{{\phi }_i}\hfill \end{array}$$

(14)

In which ${\theta }_{i0}$ and ${\phi }_{i0}$ are the true coordinates. ${\theta }_i$ and ${\phi }_i$ stand for the nearest grid point while ${\delta }_{{\theta }_i}$ and ${\delta }_{{\phi }_i}$ stand for the off-grid bias. Let $\mathsf{\Psi }=\{({\theta }_1,{\phi }_1),({\theta }_2,{\phi }_2),\dots ,({\theta }_Z,{\phi }_Z)\}$ be the uniform discretization of the 2D scene after vectorization and let the matrix ${\mathbf{\Phi }}_0$ in Equation (13) expressed as ${\mathbf{\Phi }}_0=[\mathit{a}({\theta }_1,{\phi }_1),\mathit{a}({\theta }_2,{\phi }_2),\dots ,\mathit{a}({\theta }_Z,{\phi }_Z)]$. The true position set of the scatterers is $\{({\theta }_{10},{\phi }_{10}),\dots ,({\theta }_{k0},{\phi }_{k0}),\dots ,({\theta }_{K0},{\phi }_{K0})\}$. Supposing for a certain $k\in \{1,2,\dots ,K\}$, its corresponding coordinate $({\theta }_{k0},{\phi }_{k0})\notin \mathsf{\Psi }$. The nearest grid point is $({\theta }_k,{\phi }_k),k\in \{1,2,\dots ,Z\}$. We adopt the first order Taylor series expansion approximation of $\mathit{a}({\theta }_{k0},{\phi }_{k0})$ with respect to $({\delta }_{{\theta }_k},{\delta }_{{\phi }_k})$ as:

$$\mathit{a}({\theta }_{k0},{\phi }_{k0})\approx \mathit{a}({\theta }_k,{\phi }_k)+\mathit{b}({\theta }_k,{\phi }_k)({\theta }_{k0}-{\theta }_k)+\mathit{c}({\theta }_k,{\phi }_k)({\phi }_{k0}-{\phi }_k)$$

(15)

In which $\mathit{b}({\theta }_k,{\phi }_k)$ and $\mathit{c}({\theta }_k,{\phi }_k)$ are the partial derivatives of $\mathit{a}({\theta }_k,{\phi }_k)$ with respect to ${\theta }_k$ and ${\phi }_k$.

Denote ${\mathbf{\Phi }}_{\theta }=[\mathit{b}({\theta }_1,{\phi }_1),\mathit{b}({\theta }_2,{\phi }_2),\dots ,\mathit{b}({\theta }_Z,{\phi }_Z)]$, ${\mathbf{\Phi }}_{\phi }=[\mathit{c}({\theta }_1,{\phi }_1),\mathit{c}({\theta }_2,{\phi }_2),\dots ,\mathit{c}({\theta }_Z,{\phi }_Z)]$, ${\mathbf{\delta }}_{\theta }=[{\delta }_{{\theta }_1},{\delta }_{{\theta }_2},\dots ,{\delta }_{{\theta }_Z}]$ and ${\mathbf{\delta }}_{\phi }=[{\delta }_{{\phi }_1},{\delta }_{{\phi }_2},\dots ,{\delta }_{{\phi }_Z}]$. The off-grid bias vector ${\mathbf{\delta }}_{\theta }$ satisfies

$${\delta }_{{\theta }_k}=\left\{\begin{array}{cc}{\theta }_{k0}-{\theta }_k\hfill & \mathrm{if}{\sigma }_k\ne 0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(16)

where $k=1,2,\dots ,Z$. It is the same for ${\mathbf{\delta }}_{\phi }$. Hence, the off-grid imaging model corresponding to Equation (13) can be expressed as:

$$\mathit{y}=({\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta }+{\mathbf{\Phi }}_{\phi }{\Delta }_{\phi })\mathbf{\sigma }+\mathit{e}$$

(17)

with ${\Delta }_{\theta }=diag({\mathbf{\delta }}_{\theta })$ and ${\Delta }_{\phi }=diag({\mathbf{\delta }}_{\phi })$. With this off-grid imaging model, our goal is to find the optimal $\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$ simultaneously, after which the accurate 3D image can be reconstructed.

3. The Proposed Off-Grid Imaging Algorithm

According to Equation (17), the off-grid 3D imaging problem is presented. The off-grid bias can be seen as a kind of multiplicative noise, thus the traditional sparse recovery algorithms, such as OMP and BP, fail to solve this problem. In Yang et al. [19], OGSBI is proposed and proves its potential to the off-grid DOA estimation problem. Inspired by the idea in Bishop [21], we propose an algorithm based on sparse Bayesian learning using a three-stage sparse prior to solve the off-grid imaging problem. In this section, we first describe the SBL-based off-grid imaging algorithm. Then the Bayesian Cramér-Rao bounds for off-grid bias estimator are presented.

3.1. Algorithm Description

3.1.1. The Three-Stage Sparse Prior Model

The overall graphical model of the three-stage sparse prior is shown in Figure 2.

Starting from a Bayesian perspective, all the knowns in Equation (17) are assigned with probability distributions. Then we can estimate $\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$ based on the maximum a posteriori (MAP) criterion. Note that the probability density function of a complex Gaussian distributed vector $\mathit{x}\sim \mathcal{CN}(\mathbf{\mu },\Sigma )$ with its mean $\mathbf{\mu }$ and covariance matrix $\Sigma$ is:

$$\mathcal{CN}(\mathit{x}|\mathbf{\mu },\Sigma )=\frac{1}{{\pi }^Z|\Sigma |}\exp (-{(\mathit{x}-\mathbf{\mu })}^H{\Sigma }^{-1}(\mathit{x}-\mathbf{\mu }))$$

(18)

First, we start from the sparse prior model for the radar cross section (RCS) vector $\mathbf{\sigma }$. Here we adopt a three-stage sparse prior for $\mathbf{\sigma }$. Supposing that the RCS vector $\mathbf{\sigma }$ has a complex Gaussian distribution, $\mathbf{\sigma }\sim \mathcal{CN}(0,\Sigma )$, which corresponds to the well known Swerling-1 model. In addition, the covariance matrix $\Sigma =diag(\mathbf{\alpha })$, in which $\mathbf{\alpha }={[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_Z]}^T$. The reason why we choose this probability distribution function is that a sparse prior is needed under the Bayesian framework. Since $\mathbf{\sigma }$ is sparse, that is to say, most of the elements in $\mathbf{\sigma }$ are equal or close to zero, we choose the Gaussian zero-mean model as the sparse prior. Along with the three-stage sparse prior model in the following, $\mathbf{\sigma }$ is strongly peaked at the origin. As a result, this Gaussian zero-mean prior favors that most elements of $\mathbf{\sigma }$ being zero. It is shown in Bishop [21] that both Gaussian distribution and Gamma distribution belong to the exponential distribution family and they are a pair of conjugate priors. This property leads to the posterior function having the same functional form as the prior and to a simplified Bayesian analysis. So we assume a Gamma distribution for $\mathbf{\alpha }$:

$$p(\mathbf{\alpha }|\mathbf{\beta })=\prod _{i=1}^Z\mathsf{\Gamma }({\alpha }_i|1,{\beta }_i)$$

(19)

with $\mathsf{\Gamma }(x|a,b)={\left[\mathsf{\Gamma }(a)\right]}^{-1}{b}^a{x}^{a-1}{e}^{-bx}$ and $\mathsf{\Gamma }(·)$ being the Gamma function. The first stage of the sparse prior model is that $\mathbf{\beta }$ also follows a Gamma distribution with $a,b>0$:

$$p(\mathbf{\beta }|a,b)=\prod _{i=1}^Z\mathsf{\Gamma }({\beta }_i|a,b)$$

(20)

As is shown in Figure 2, this three-stage prior is a sparse prior of $\mathbf{\sigma }$. Compared with the method in Yang et al. [19], this model introduces Z more hyperparameters into the sparse recovery model to offer more degrees of freedom.

Secondly, assuming that the noise is white complex Gaussian, we have:

$$p(\mathit{e}|\eta )=\mathcal{CN}(\mathit{e}|0,\eta \mathit{I})$$

(21)

where $\eta$ being the noise variance. Similarly, a conjugate prior is assigned that $p(\eta |c,d)=\mathsf{\Gamma }(\eta |c,d)$. It is worth mentioning here why the sensing noise $\mathit{e}$ is assumed zero-mean. In Equation (15), a first order Taylor series expansion approximation is adopted and higher order items are neglected. So sensing noise in Equation (17) only contains the measurement noise which is assumed white Gaussian distributed with zero mean.

Finally, for the off-grid bias ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\varphi }$, the uniform distribution is more suitable than the Gaussian distribution since we do not have any information about them. Supposing that the pre-discretized grid intervals are ${\rho }_{\theta }$ and ${\rho }_{\phi }$, we get:

$${\mathbf{\delta }}_{\theta }\sim U\left({[-\frac{{\rho }_{\theta }}{2},\frac{{\rho }_{\theta }}{2}]}^Z\right){\mathbf{\delta }}_{\phi }\sim U\left({[-\frac{{\rho }_{\phi }}{2},\frac{{\rho }_{\phi }}{2}]}^Z\right)$$

(22)

3.1.2. Variational Inference EM Based Sparse Recovery Algorithm

Under the framework of sparse Bayesian learning, we always seek to find the optimal estimations from the MAP criterion. The posterior distribution can be expressed by the following equations:

$$p(\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta |\mathit{y})=p(\mathit{y},\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta )/p(\mathit{y})$$

(23)

$$p(\mathit{y},\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta )=p(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )p(\mathbf{\sigma }|\mathbf{\alpha })p(\mathbf{\alpha }|\mathbf{\beta })p(\mathbf{\beta })p({\mathbf{\delta }}_{\theta })p({\mathbf{\delta }}_{\phi })p(\eta )$$

(24)

$$p(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )=\mathcal{CN}(\mathit{y}|({\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta }+{\mathbf{\Phi }}_{\phi }{\Delta }_{\phi })\mathbf{\sigma },\eta \mathit{I})$$

(25)

However, according to Equation (24) the marginalized likelihood function

$$p(\mathit{y})=\oint p(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )p(\mathbf{\sigma }|\mathbf{\alpha })p(\mathbf{\alpha }|\mathbf{\beta })p(\mathbf{\beta })p({\mathbf{\delta }}_{\theta })p({\mathbf{\delta }}_{\phi })p(\eta )d\mathbf{\sigma }d{\mathbf{\delta }}_{\theta }d{\mathbf{\delta }}_{\phi }d\mathbf{\alpha }d\mathbf{\beta }d\eta$$

(26)

which cannot be calculated analytically, thus the normalized constant in Equation (23) cannot be computed. In Blei et al. [22], the author gives a solution to this kind of problem using variational inference EM algorithm. Hence, we adopt the variational EM algorithm in this sparse recovery problem. According to the jargon in the EM algorithm, we treat $\mathbf{\sigma },\mathbf{\alpha },\mathbf{\beta },\eta$ as the hidden variables while ${\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$ being the parameters. Since the variational inference EM algorithm is a two-stage iterative optimization algorithm, assuming that ${\mathbf{\sigma }}^{(i)},{\mathbf{\alpha }}^{(i)},{\mathbf{\beta }}^{(i)},{\eta }^{(i)},{\mathbf{\delta }}_{\theta }^{(i)},{\mathbf{\delta }}_{\phi }^{(i)}$ are obtained in the ith iteration, we now seek for their updates in the $(i+1)$th iteration.

(1) The E-step

Since the posterior distribution $p(\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta |\mathit{y})$ cannot be calculated analytically, we use the variational inference technique in Tzikas et al. [23] to find an approximation for the posterior distribution:

$$p(\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta |\mathit{y})\approx q(\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\mathbf{\alpha },\mathbf{\beta },\eta |\mathit{y}):=q(\mathbf{\sigma })q({\mathbf{\delta }}_{\theta })q({\mathbf{\delta }}_{\phi })q(\mathbf{\alpha })q(\mathbf{\beta })q(\eta )$$

(27)

which minimizes the Kullback-Leibler divergence. This approximation comes from the variational methods used in Bayesian inference. Among these applications, a particular form that has been used with great success is the factorized one which is used in this paper. In a broad sense, this approximation is always valid, but the performance is effected by the form of the prior distributions. In our method, we use the conjugate priors in order to get an analytical expression, and the performance of the approximation in Equation (27) is satisfactory. Based on the above equation, we can get the estimation of the hidden variables. The detailed procedures are as follows.

For $\mathbf{\sigma }$, we have:

$$\begin{array}{cc}\hfill lnq(\mathbf{\sigma })& =E{[lnp(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )+lnp(\mathbf{\sigma }|\mathbf{\alpha })]}_{q({\mathbf{\delta }}_{\theta })q({\mathbf{\delta }}_{\phi })q(\mathbf{\alpha })q(\eta )}+C\hfill \\ & \propto -{\mathbf{\sigma }}^H[\widehat{\eta }{\mathbf{\Phi }}^H\mathbf{\Phi }+diag(1/\widehat{\mathbf{\alpha }})]\mathbf{\sigma }-2\widehat{\eta }{\mathit{y}}^H\mathbf{\Phi }\mathbf{\sigma }\hfill \end{array}$$

(28)

where $\mathbf{\Phi }={\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }diag({\widehat{\mathbf{\delta }}}_{\theta })+{\mathbf{\Phi }}_{\phi }diag({\widehat{\mathbf{\delta }}}_{\phi })$. According to the above equation, $q(\mathbf{\sigma })$ is a complex Gaussian distribution, i.e., $\mathbf{\sigma }\sim \mathcal{CN}(\mathbf{\mu },\Sigma )$ with

$$\begin{array}{cc}\hfill \mathbf{\mu }& =\widehat{\eta }\Sigma {\mathbf{\Phi }}^H\mathit{y}\hfill \\ \hfill \Sigma & ={[\widehat{\eta }{\mathbf{\Phi }}^H\mathbf{\Phi }+diag(1/\widehat{\mathbf{\alpha }})]}^{-1}\hfill \end{array}$$

(29)

Since we have already obtained both hidden variables and parameters in the ith iteration, the expectations in Equation (29) are all replaced by the ith updates, i.e.,

$$\begin{array}{cc}\hfill {\mathbf{\mu }}^{(i+1)}({\mathbf{\sigma }}^{(i+1)})& ={\eta }^{(i)}{\Sigma }^{(i+1)}{\mathbf{\Phi }}^H\mathit{y}\hfill \\ \hfill {\Sigma }^{(i+1)}& ={[{\eta }^{(i)}{\mathbf{\Phi }}^H\mathbf{\Phi }+diag(1/{\mathbf{\alpha }}^{(i)})]}^{-1}\hfill \\ \hfill \mathbf{\Phi }& ={\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }diag({\mathbf{\delta }}_{\theta }^{(i)})+{\mathbf{\Phi }}_{\phi }diag({\mathbf{\delta }}_{\phi }^{(i)})\hfill \end{array}$$

(30)

In the following part, the expectations of these unknowns are replaced by their corresponding latest updates.

For $\mathbf{\alpha }$, we have:

$$\begin{array}{cc}\hfill lnq(\mathbf{\alpha })& =E{[lnp(\mathbf{\sigma }|\mathbf{\alpha })+lnp(\mathbf{\alpha }|\mathbf{\beta })]}_{q(\mathbf{\sigma })q(\mathbf{\beta }})+C\hfill \\ \hfill lnq({\mathbf{\alpha }}_j)& \propto -\widehat{{\sigma }_j^2}{\alpha }_j^{-1}-ln{\alpha }_j-{\widehat{\beta }}_j{\alpha }_j\hfill \end{array}$$

(31)

Thus, the elements in $\mathbf{\alpha }$ are independent and their distribution is a generalized inverse Gaussian distribution [24]. Here we use the expectation of ${\mathbf{\alpha }}_j$ to get the $(i+1)$th update for ${\mathbf{\alpha }}^{(i+1)}$. For the jth element in $\mathbf{\alpha }$:

$${\alpha }_j^{(i+1)}={(\frac{{\xi }_j^{(i+1)}}{{\beta }_j^{(i)}})}^{1/2}\frac{{\kappa }_1(2\sqrt{{\beta }_j^{(i)}{\xi }_j^{(i+1)})}}{{\kappa }_0(2\sqrt{{\beta }_j^{(i)}{\xi }_j^{(i+1)})}}$$

(32)

In which ${\xi }_j^{(i+1)}={({\sigma }_j^2)}^{(i+1)}={|{\sigma }_j^{(i+1)}|}^2+{\Sigma }_{j,j}^{(i+1)}$.

For $\mathbf{\beta }$, we have:

$$\begin{array}{cc}\hfill lnq(\mathbf{\beta })& =E{[lnp(\mathbf{\alpha }|\mathbf{\beta })+lnp(\mathbf{\beta })]}_{q(\mathbf{\alpha })}+C\hfill \\ \hfill lnq({\mathbf{\beta }}_j)& \propto aln{\beta }_j-(b+\widehat{{\alpha }_j}){\beta }_j\hfill \end{array}$$

(33)

Similarly, elements in $\mathbf{\beta }$ are independent and $p({\beta }_j)=\mathsf{\Gamma }({\beta }_j|a+1,b+{\alpha }_j)$. The $(i+1)$th update of $\mathbf{\beta }$ are as follows:

$${\beta }_j^{(i+1)}=\frac{a+1}{b+{\alpha }_j^{(i+1)}}$$

(34)

For $\eta$, we have:

$$\begin{array}{cc}\hfill lnq(\eta )& =E{[lnp(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )+lnp(\eta )]}_{q(\mathbf{\alpha })q({\mathbf{\delta }}_{\theta })q({\mathbf{\delta }}_{\phi })}+C\hfill \\ & \propto (Z+c-1)ln\eta -\left({\parallel \mathit{y}-\mathbf{\Phi }\mathbf{\mu }\parallel }_F^2+tr({\mathbf{\Phi }}^H\mathbf{\Phi }\Sigma )+d\right)\eta \hfill \end{array}$$

(35)

where $\mathbf{\mu }$ and $\Sigma$ are defined in Equation (29). Similar to $\mathbf{\beta }$, $q(\eta )$ is also a Gamma distribution. Combining the expression for $\Sigma$ in Equation (29) and the expectation for a Gamma distribution, we can get:

$${\eta }^{(i+1)}=\frac{Z+c}{\parallel \mathit{y}-{\mathbf{\Phi }\mathbf{\mu }}^{(i+1)}{\parallel }_F^2+d+\frac{1}{{\eta }^{(i)}}{\sum }_{n=1}^Z(1-\frac{1}{{\alpha }_n^{(i+1)}}{\Sigma }_{n,n}^{(i+1)})}$$

(36)

(2) The M-step

For ${\mathbf{\delta }}_{\theta }$, its estimate maximize $E[lnp(\mathit{y}|\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi },\eta )+lnp({\mathbf{\delta }}_{\theta })]$, which equivalent to minimizing

$$\begin{array}{c}E\left\{\parallel \mathit{y}-\left[{\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }diag({\mathbf{\delta }}_{\theta })+{\mathbf{\Phi }}_{\phi }diag({\mathbf{\delta }}_{\phi })\right]{\mathbf{\sigma }\parallel }_2^2\right\}\hfill \\ ={\mathbf{\delta }}_{\theta }^T{\mathit{A}}_1{\mathbf{\delta }}_{\theta }+2{\mathit{a}}_1^T{\mathbf{\delta }}_{\theta }+C\hfill \end{array}$$

(37)

in which

$$\begin{array}{cc}\hfill {\mathit{A}}_1& =\mathsf{\Re }\left\{\overline{{\mathbf{\Phi }}_1^H{\mathbf{\Phi }}_1}\odot ({\Sigma }^{(i+1)}+{\mathbf{\mu }}^{(i+1)}{({\mathbf{\mu }}^{(i+1)})}^H)\right\}\hfill \\ \hfill {\mathit{a}}_1& =\mathsf{\Re }\left\{diag({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_1{\Sigma }^{(i+1)})-diag({\mathbf{\mu }}^{(i+1)}){\mathbf{\Phi }}_{\theta }^H(\mathit{y}-{\mathbf{\Phi }}_1{\mathbf{\mu }}^{(i+1)})\right\}\hfill \\ \hfill {\mathbf{\Phi }}_1& ={\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\phi }diag({\mathbf{\delta }}_{\phi }^{(i)})\hfill \end{array}$$

(38)

Hence, the update of ${\mathbf{\delta }}_{\theta }$ in the $(i+1)$th iteration can be calculated by minimizing the expression in Equation (37).

For ${\mathbf{\delta }}_{\phi }$, it is the same as ${\mathbf{\delta }}_{\theta }$, so we directly give the result:

$${\mathbf{\delta }}_{\phi }^{(i+1)}=\underset{{\mathbf{\delta }}_{\phi }}{argmin}\left\{{\mathbf{\delta }}_{\phi }^T{\mathit{A}}_2{\mathbf{\delta }}_{\phi }+2{\mathit{a}}_2^T{\mathbf{\delta }}_{\phi }\right\}$$

(39)

where ${\mathit{A}}_2$, ${\mathit{a}}_2$ and ${\mathbf{\Phi }}_2$ are similar to those for ${\mathbf{\delta }}_{\theta }$. Notice that in the maximization step, both ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\phi }$ are jointly sparse with $\mathbf{\sigma }$, leading to the dimensionality reduction and significant decrease of the computation load. Based on the above analysis, the variational EM-based off-grid imaging algorithm can be described in

Table 1. Main steps of the proposed imaging method.

Processing Steps of the Off-Grid Imaging Method.
Input: $\mathit{y},{\mathbf{\Phi }}_0,{\mathbf{\Phi }}_{\theta },{\mathbf{\Phi }}_{\phi }$
Initialization:
${\mathbf{\sigma }}_i^{(0)}={\mathbf{\Phi }}_i^H\mathit{y}/{\mathbf{\Phi }}_i^H{\mathbf{\Phi }}_i$, ${\eta }^0=\left|\right|\mathit{y}-\mathbf{\Phi }{\mathbf{\sigma }}^{(0)}{\left|\right|}^2/MN$, ${\mathbf{\delta }}_{\theta }={\mathbf{\rho }}_{\theta }/2$, ${\mathbf{\delta }}_{\varphi }={\mathbf{\rho }}_{\phi }/2$, $a=b=c=d={10}^{-5}$
Iteration: let i denotes the iteration counter
(1) update ${\mathbf{\sigma }}^{(i+1)}$ by Equation (30)
(2) update ${\mathbf{\alpha }}^{(i+1)}$ by Equation (32)
(3) update ${\mathbf{\beta }}^{(i+1)}$ by Equation (34)
(4) update ${\eta }^{(i+1)}$ by Equation (36)
(5) update ${\mathbf{\delta }}_{\theta }^{(i+1)}$ and ${\mathbf{\delta }}_{\phi }^{(i+1)}$ by Equations (37) and (39)
Termination condition:
The iteration ends when $\left|\right|{\mathbf{\sigma }}^{(i+1)}-{\mathbf{\sigma }}^{(i)}{\left|\right|}^2/\left|\right|{\mathbf{\sigma }}^{(i)}{\left|\right|}^2<{10}^{-5}$.
Output:
The imaging result ${\mathbf{\sigma }}^{(i+1)}$, the off-gird bias estimation ${\mathbf{\delta }}_{\theta }^{(i+1)}$ and ${\mathbf{\delta }}_{\phi }^{(i+1)}$

. It is worth mentioning that this algorithm does not need the sparsity K as prior information. This property makes it different from traditional super-resolution methods such as MUSIC or ESPRIT. Since the true sparsity of the imaging scene is often not available, this algorithm has vast application prospect.

The overall 3D imaging method is shown in Figure 3. First, the channel separation and range compression are conducted. In practical applications, the range resolution provided by the wide bandwidth is sufficient. After we get the echo signal, the range compression technique is used to get the one dimensional range profile first. For range cells where scatterers exist, the algorithm will further solve the two dimensional imaging problem in the elevation and azimuth direction. Thus, the 3D imaging problem is reduced to a series of off-grid imaging problems in the elevation and azimuth direction. By using the proposed method, these problems can be solved and the 3D image is finally reconstructed. The proposed 3D imaging algorithm adopts the uniform prior, which is more suitable than the Gaussian distribution in S-TLS. This algorithm seeks to find the optimal estimation according to the MAP criterion and it can be seen as an extension to the imaging problem with higher dimensionality. Moreover, more degrees of freedom are introduced into the model via a three-stage sparse prior, which will improve the imaging performance.

3.2. Bayesian Cramér-Rao Bounds For Off-Grid Biases

It is well known that the Cramér-Rao Lower Bound(CRLB) is an effective indicator for the MSE performance of unbiased estimators. In Prasad and Murthy [25], the author provides an analogous bound called Bayesian Cramér-Rao Bound to provide a lower bounds for estimation problems in sparse Bayesian learning. Different from the traditional CRLB, the prior distribution of the unknowns are considered. The MSEs of ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\phi }$ are compared with their corresponding BCRBs. This estimation algorithm is seen as effective if their MSEs approach the BCRB with the increase of signal-to-noise ratio (SNR).

In this off-grid imaging algorithm, $\mathbf{\sigma }$, ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\phi }$ are the unknowns to be estimated. We denote a new vector $\Theta =[\mathbf{\sigma };{\mathbf{\delta }}_{\theta };{\mathbf{\delta }}_{\phi }]$ containing the all the unknowns to be estimated and the MSE matrix is defined as ${\mathit{E}}^{\Theta }\triangleq E_{\mathit{y},\Theta }\left[(\Theta -\widehat{\Theta }(\mathit{y})){(\Theta -\widehat{\Theta }(\mathit{y}))}^T\right]$. The first step to calculate the Cramér-Rao Lower Bound is to derive the Fisher Information Matrix(FIM) ${\mathit{I}}^{\Theta }$. Usually it is convenient to express ${\mathit{I}}^{\Theta }$ in terms of submatrices, in which the $(i,j)$th block is as follows [25]:

$${\mathit{I}}_{ij}^{\Theta }\triangleq -E_{\mathit{y},\Theta }\left[{\nabla }_{{\Theta }_{\mathit{i}}}{\nabla }_{{\Theta }_{\mathit{i}}}^T\log p(\mathit{y},\Theta ;\eta ,\mathbf{\alpha })\right]$$

(40)

Thus, the Fisher information matrix can be written as:

$${\mathit{I}}^{\Theta }=\left(\begin{array}{ccc}{\mathit{I}}_{11}& {\mathit{I}}_{12}& {\mathit{I}}_{13}\\ {\mathit{I}}_{21}& {\mathit{I}}_{22}& {\mathit{I}}_{23}\\ {\mathit{I}}_{31}& {\mathit{I}}_{32}& {\mathit{I}}_{33}\end{array}\right)$$

(41)

Only consider the items that are relevant to $\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$, we have:

$$-\log p(\mathit{y},\Theta ;\eta ,\mathbf{\alpha })=\frac{1}{\eta }{\parallel \mathit{y}-({\mathbf{\Phi }}_0+{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta }+{\mathbf{\Phi }}_{\phi }{\Delta }_{\phi })\mathbf{\sigma }\parallel }_2^2+{\mathbf{\sigma }}^H{\Sigma }^{-1}\mathbf{\sigma }$$

(42)

Substituting Equation (42) into Equation (40), we can get the following results. For ${\mathit{I}}_{11}$, noting that only quadratic items of $\mathbf{\sigma }$ are useful and the expectations of ${\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$ are $\mathbf{0}$, we can get the equation:

$$\begin{array}{c}{\mathit{I}}_{11}=-E_{\mathit{y},\Theta }\left[{\nabla }_{\mathbf{\sigma }}{\nabla }_{\mathbf{\sigma }}^T\log p(\mathit{y},\Theta ;\eta ,\mathbf{\alpha })\right]\hfill \\ =E_{\mathit{y},\Theta }\left[\frac{2}{\eta }({\mathbf{\Phi }}_0^H{\mathbf{\Phi }}_0+{\Delta }_{\theta }{\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta }+{\Delta }_{\phi }{\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi }{\Delta }_{\phi }+\eta {\Sigma }^{-1})\right]\hfill \\ =\frac{2}{\eta }[{\mathbf{\Phi }}_0^H{\mathbf{\Phi }}_0+\frac{{\rho }_{\theta }^2}{12}(\mathit{I}\odot {\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })+\frac{{\rho }_{\phi }^2}{12}(\mathit{I}\odot {\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })+\eta {\Sigma }^{-1}]\hfill \end{array}$$

(43)

As for ${\mathit{I}}_{12},{\mathit{I}}_{13},{\mathit{I}}_{21},{\mathit{I}}_{31}$, since $\mathbf{\sigma },{\mathbf{\delta }}_{\theta },{\mathbf{\delta }}_{\phi }$ are independent, it is easy to find that they are all zero matrices. For ${\mathit{I}}_{22}$ and ${\mathit{I}}_{33}$, the derivations are similar. Here we take ${\mathit{I}}_{22}$ as an example. Similar to that of $\mathbf{\sigma }$, only quadratic items of ${\mathbf{\delta }}_{\theta }$ contribute to the final result.

$$\begin{array}{c}{\mathit{I}}_{22}=E_{\mathit{y},\Theta }\left[{\nabla }_{{\mathbf{\delta }}_{\theta }}{\nabla }_{{\mathbf{\delta }}_{\theta }}^T\frac{1}{\eta }{\mathbf{\sigma }}^H({\Delta }_{\theta }{\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta })\mathbf{\sigma }\right]\hfill \\ =\frac{1}{\eta }{\nabla }_{{\mathbf{\delta }}_{\theta }}{\nabla }_{{\mathbf{\delta }}_{\theta }}^T{E}_{\mathit{y},\Theta }\left[tr({\mathbf{\sigma }}^H({\Delta }_{\theta }{\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta })\mathbf{\sigma })\right]\hfill \end{array}$$

(44)

By exchanging the expectation and the trace operator, the above equation reduces to $\frac{1}{\eta }{\nabla }_{{\mathbf{\delta }}_{\theta }}{\nabla }_{{\mathbf{\delta }}_{\theta }}^{T}tr(\Sigma {\Delta }_{\theta }{\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta })$. Notice that $tr(\Sigma diag({\mathbf{\delta }}_{\theta }){\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta }diag({\mathbf{\delta }}_{\theta }))={\mathbf{\delta }}_{\theta }^T(\Sigma \odot {({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}^T){\mathbf{\delta }}_{\theta }$. Besides, by observing that $\Sigma \odot {({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}^T$ is positive semi-definite and that ${\mathbf{\delta }}_{\theta }$ is real-valued, we can get ${\mathbf{\delta }}_{\theta }^T(\Sigma \odot {({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}^T){\mathbf{\delta }}_{\theta }={\mathbf{\delta }}_{\theta }^T\mathsf{\Re }(\Sigma \odot {({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}^T){\mathbf{\delta }}_{\theta }$. This leads to the result that ${\mathit{I}}_{22}=\frac{2}{\eta }\mathsf{\Re }(\Sigma \odot {\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })$. The derivation process of ${\mathit{I}}_{33}$ is the same.

Similarly, we only present the derivation process of ${\mathit{I}}_{23}$ since that of ${\mathit{I}}_{32}$ is the same. For ${\mathit{I}}_{23}$:

$$\begin{array}{c}{\mathit{I}}_{23}=\frac{1}{\eta }{E}_{\mathit{y},\Theta }\left[{\nabla }_{{\mathbf{\delta }}_{\theta }}{\nabla }_{{\mathbf{\delta }}_{\phi }}^T{\mathbf{\sigma }}^H({\Delta }_{\theta }{\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi }{\Delta }_{\phi }+{\Delta }_{\phi }{\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\theta }{\Delta }_{\theta })\mathbf{\sigma }\right]\hfill \\ =\frac{1}{\eta }({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi }+{\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\theta })\odot \Sigma \hfill \\ =\frac{2}{\eta }\Sigma \odot \mathsf{\Re }({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi })\hfill \end{array}$$

(45)

The submatrices of the FIM are listed here:

$$\begin{array}{c}{\mathit{I}}_{11}=\frac{2}{\eta }[{\mathbf{\Phi }}_0^H{\mathbf{\Phi }}_0+\frac{{\rho }_{\theta }^2}{12}(\mathit{I}\odot {\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })+\frac{{\rho }_{\phi }^2}{12}(\mathit{I}\odot {\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })+\eta {\Sigma }^{-1}]\hfill \\ {\mathit{I}}_{12},{\mathit{I}}_{13},{\mathit{I}}_{21},{\mathit{I}}_{31}=\mathbf{0}\hfill \\ {\mathit{I}}_{22}=\frac{2}{\eta }\mathsf{\Re }(\Sigma \odot {\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })\hfill \\ {\mathit{I}}_{23}={\mathit{I}}_{32}=\frac{2}{\eta }\Sigma \odot \mathsf{\Re }({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi })\hfill \\ {\mathit{I}}_{33}=\frac{2}{\eta }\mathsf{\Re }(\Sigma \odot {\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })\hfill \end{array}$$

(46)

Then a lower bound on the MSE matrix ${\mathit{E}}^{\Theta }$ is presented by the inversion of the FIM ${\mathit{I}}^{\Theta }$:

$${\mathit{E}}^{\Theta }-{({\mathit{I}}^{\Theta })}^{-1}⪰\mathbf{0}$$

(47)

where $⪰\mathbf{0}$ is interpreted as meaning that the matrix is positive semidefinite. Noting that ${\mathit{I}}_{22}$, ${\mathit{I}}_{23}$, ${\mathit{I}}_{32}$ and ${\mathit{I}}_{33}$ are diagonal matrices, ${({\mathit{I}}^{\Theta })}^{-1}$ can be expressed as follows:

$${\mathit{J}}^{\Theta }={({\mathit{I}}^{\Theta })}^{-1}=\left(\begin{array}{ccc}{\mathit{J}}_{11}& {\mathit{J}}_{12}& {\mathit{J}}_{13}\\ {\mathit{J}}_{21}& {\mathit{J}}_{22}& {\mathit{J}}_{23}\\ {\mathit{J}}_{31}& {\mathit{J}}_{32}& {\mathit{J}}_{33}\end{array}\right)$$

(48)

where ${\mathit{J}}_{11}={\mathit{I}}_{11}^{-1}$ and ${\mathit{J}}_{12},{\mathit{J}}_{13},{\mathit{J}}_{21},{\mathit{J}}_{31}=\mathbf{0}$ according to Equation (46). Since ${\mathit{I}}_{22}$, ${\mathit{I}}_{23}$, ${\mathit{I}}_{32}$ and ${\mathit{I}}_{33}$ are all diagonal matrices, we find that ${\mathit{J}}_{22}$, ${\mathit{J}}_{33}$ and ${\mathit{J}}_{23}$ are also diagonal matrices and their elements can be easily calculated by solving a set of linear equations. For a positive semidefinite matrix, its diagonal elements are nonnegative. Thus, the BCRBs of the off-grid bias ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\phi }$ are obtained. For the ith element of ${\mathbf{\delta }}_{\theta }$ and ${\mathbf{\delta }}_{\phi }$, the expressions are:

$$\begin{array}{c}\hfill E(|{\widehat{\delta }}_{{\theta }_i}-{\delta }_{{\theta }_i}{|}^2)\ge \frac{\mathsf{\Re }[{({\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })}_{i,i}]}{2\left\{\mathsf{\Re }[{({\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })}_{i,i}]\mathsf{\Re }[{({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}_{i,i}]-{\mathsf{\Re }}^2[{({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi })}_{i,i}]\right\}\mathrm{SNR}}\\ \hfill E(|{\widehat{\delta }}_{{\phi }_i}-{\delta }_{{\phi }_i}{|}^2)\ge \frac{\mathsf{\Re }[{({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}_{i,i}]}{2\left\{\mathsf{\Re }[{({\mathbf{\Phi }}_{\phi }^H{\mathbf{\Phi }}_{\phi })}_{i,i}]\mathsf{\Re }[{({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\theta })}_{i,i}]-{\mathsf{\Re }}^2[{({\mathbf{\Phi }}_{\theta }^H{\mathbf{\Phi }}_{\phi })}_{i,i}]\right\}\mathrm{SNR}}\end{array}$$

(49)

where the SNR is defined as $\frac{{\alpha }_i}{\eta }$.

4. Experimental Results

In this section, we present some simulation results to analyze performance of the proposed algorithm. In Section 4.1, a few sparse recovery algorithms are applied to the 3D imaging problem and the results are compared. The 3D imaging performance, i.e., the normalized mean square error(NMSE) of target image recovery and off-grid bias estimation are presented. Besides, the imaging results of a complex model is provided to further validate the feasibility of proposed method. Then super-resolution ability is analyzed with respect to SNR in Section 4.2. Finally the off-grid bias estimations are compared with corresponding BCRBs in Section 4.3. The simulation parameters in the sparse array 3D imaging system are listed in

Table 2. Simulation Conditions.

Parameter	Symbol	Value
Bandwidth	B	500 MHz
Carrier frequency	$f_c$	10 GHz
Baseline in the X-direction	$L_x$	6 m
Baseline in the Y-direction	$L_y$	6 m
Target range	R	2500 m
Pulse width	$T_r$	10 $\mathsf{\mu }$s
Number of transmitters	M	4
Number of receivers	N	225

.

4.1. Validation of The Proposed Algorithm

A simulated target with 24 scatterers is considered in this part. Since the range resolution is guaranteed through range compression, for simplicity, the scatterers of the target are located at three different range cells. The original 3D model is shown in Figure 4a, and both on-grid and off-grid scatterers are included. The sparse antenna array is distributed on a circular area with its diameter being 6 m and center being the origin. Four transceivers are located at the intersections of the boundary and the axes. The other receivers are first uniformly distributed on this area, then random disturbances are added to their positions in order to alleviate the ambiguity problem which emerges in sparse array. The center of target is located at $(2500$ m, $\frac{\pi }{4},\frac{\pi }{4})$.

3D imaging results using different sparse recovery algorithms with SNR equals 20 dB are presented in Figure 4, with blue circle and red circle representing the true scatterer and the reconstructed scatterer, respectively. The size of the circle is decided by the scatterer’s RCS. As discussed before, the classic sparse recovery algorithms including OMP and BP fail to reconstruct the 3D imaging because of the mismatch. Compared with the classic algorithms, S-TLS takes the off-grid bias into consideration and gets a better imaging performance. However, there exist some false scatterers and the off-grid bias estimator does not perform well. This can be explained by the Gaussian prior used for the off-grid biases while their true distribution being uniform. Different from S-TLS, OGSBI adopts a uniform prior and its imaging performance is better. Yet there still exists mismatch between the scatterers’ true position and reconstructed position. In contrast, the imaging performance of the proposed algorithm is the best among all methods because of the application of the three-stage sparse prior.

The NMSEs using these different algorithms are summarized in

Table 3. NMSEs of 3D image and off-grid bias estimation.

Algorithm	OMP	BP	S-TLS	OGSBI	Proposed Method
NMSE of $\mathbf{\sigma }$	1.3371	1.1129	0.6921	0.5397	0.3016
NMSE of ${\mathbf{\delta }}_{\theta }$	*	*	0.6048	0.2452	0.1333
NMSE of ${\mathbf{\delta }}_{\phi }$	*	*	0.6321	0.2283	0.1266

. NMSE is defined as follows and I is the number of Monte Carlo trails:

$$NMSE(\mathbf{\sigma })=\frac{1}{I}\sum \frac{\left|\right|\widehat{\mathbf{\sigma }}-\mathbf{\sigma }{\left|\right|}^2}{{\left|\right|\mathbf{\sigma }\left|\right|}^2}NMSE({\mathbf{\delta }}_{\theta ,\phi })=\frac{1}{I}\sum \frac{\left|\right|{\widehat{\mathbf{\delta }}}_{\theta ,\phi }-{\mathbf{\delta }}_{\theta ,\phi }{\left|\right|}^2}{\left|\right|{\mathbf{\delta }}_{\theta ,\phi }{\left|\right|}^2}$$

(50)

In Table 3, the asterisks mean that these algorithms cannot estimate the off-grid biases. From this table, it is clear that the proposed method obtains the best performance among all the algorithms. However, we also find that the NMSEs of 3D image is larger than those of off-grid bias estimation, which implies that current algorithms fail to estimate the RCS as accurate as position. It is also shown that the differences between $\theta$ and $\phi$ are small. So in the following part we only present the NMSE of ${\mathbf{\delta }}_{\theta }$ for brevity.

NMSEs of sparse recovery with different SNRs and grid intervals are presented. Figure 5a,b shows the NMSEs under different SNR levels. The grid interval is set to half the Rayleigh resolution. The target model is the same and the SNR increases from 0dB to 45 dB with an interval of 5 dB. Number of the Monte-Carlo trails is set to 100. Since the imaging procedure can be seen as a joint estimation problem, the correlation between $\theta$ and $\phi$ in Equation (17) makes the performances more sensitive to noise. As we can see, the NMSEs stay at a high level for all algorithms when the SNR is low. However, the image recovery error of both OMP and BP are always high even if the SNR increases. The reason is that these two algorithms fail when the mismatch problem emerges. The performance of S-TLS is better than the former two algorithms. However, it is not satisfying because of the mismatch between Gaussian distribution and the true model. OGSBI can achieve a relatively good performance. The proposed algorithm can get the best performance among all the algorithms when the SNR is larger than 20 dB because the three-stage hierarchical model introduces more degrees of freedom and is a better approximation of the $l_0$-norm optimization.

In our study, we find that the imaging performances vary with different grid sizes. So we analyze the performances versus different grid intervals in Figure 5c,d. The simulation parameters are all the same but the grid interval varies. $\chi =\frac{\rho }{R_{Rayleigh}}$ is the ratio of grid interval to the Rayleigh limit. It is shown in Figure 5 that the NMSEs of all algorithms decrease first and then increase. The reason for this phenomenon is that denser grids can alleviate the mismatch problem to some extent, however, much too dense grids will lead to a false recovery result due to violation of RIP. It is interesting that the NMSEs of OMP and BP get the minimum when $\chi$ equals 0.5 while for the rest 0.4. This is caused by the fact that these algorithms incorporate the mismatch factor into sparse recovery model. So the off-grid biases can be estimated and compensated, leading to a better performance when small grid interval is used. However, the NMSEs increase with $\chi$ increasing because the approximation error caused by the first order Taylor expansion also increases.

To further validate the practicality of proposed method for complex target 3D reconstruction, the imaging results based on a complex plane model is presented. This plane model is based on the RCS reconstruction result of a real aircraft at the airport. We use the model to generate echo signal and then use the proposed method to reconstruct the 3D plane model. In the future, we would further validate the proposed method using real-world radar data. In this simulation, the target center is located at $(2400,\frac{\pi }{4},\frac{\pi }{4})$ and the system configuration is the same as before. The reconstructed 3D image of the plane is presented in Figure 6. Figure 6a is an overall view of the 3D image while Figure 6b–d are the three views of the reconstructed target model. The size of the target can be calculated using these figures. Furthermore, features of the plane, such as the nose, engines, vertical and horizontal stabilizer, can be recognized from the reconstructed 3D image as is shown in Figure 6a,c. These results show its potential for target identification.

4.2. Super-Resolution Performance Versus SNR

When it comes to imaging algorithm, one key problem is its resolution. Since the range resolution is obtained through pulse compression, we consider the resolutions of the other two dimensions. In this part, the grid interval is set to the best in Section 4.1 and only two scatterers are considered. One experiment is defined as successful if the two scatterers are separated. Then the success rate is calculated through 100 Monte Carlo simulations.

Since the resolutions in elevation and azimuth directions are similar, only the success rate for $\theta$ are presented for brevity. In Figure 7, $\omega =\frac{d}{R_{Rayleigh}}$ is the normalized distance between the two scatterers. Simulations are conducted with different SNRs and scatterer distances. This plot gives us some perspectives about the problem. First, the super-resolution ability of proposed method is better than that of OGSBI. The reason is that the off-grid bias estimation is more accurate, making the separation of two closely spaced off-grid scatterers possible. Second, with the increasing of SNR, super-resolution ability of the proposed algorithm is increasing. Here we define the resolution as the distance between two scatterers when the success rate is more than $50\%$. The proposed algorithm can realize super-resolution with its resolutions of elevation being a quarter of the Rayleigh limit when the SNR is 20 dB. Considering the shape of scatter plot, we adopt the Boltzmann fit. The red line and black line are fitted curves for proposed method and OGSBI. The model of Boltzmann fit is:

$$y=A_2+\frac{A_1-A_2}{1+exp(\frac{x-x_0}{dx})}$$

(51)

The fitting parameters for the two algorithms under different SNRs are listed in

Table 4. Fitting Parameters.

	${\mathit{A}}_1$	${\mathit{A}}_2$	${\mathit{x}}_0$	$\mathit{dx}$
OGSBI, SNR = 10 dB	0.0064	0.9932	0.5682	0.0490
OGSBI, SNR = 20 dB	−0.0064	0.9949	0.2857	0.0412
Proposed method, SNR = 10 dB	0.0014	0.9948	0.4726	0.0418
Proposed method, SNR = 20 dB	−0.0080	0.9975	0.2553	0.0435

. According to Equation (51), both $A_1$ and $A_2$ are the normalization factors and $A_2$ can be regarded as the max value of the curve. $dx$ is a constant which influences the sharpness of the curve. It is insensitive to the SNR level according to Table 3. The success rate is around 0.5 when x equals to $x_0$, which can be seen as the reciprocal of the super-resolution factor. From Table 3 we can see it is inversely proportion to the SNR.

4.3. BCRB for Off-Grid Biases

In this subsection, the Bayesian Cramér-Rao Bound for off-grid bias is presented, then the root mean square error(RMSE) of the proposed algorithm is compared with it. The expression of RMSE is:

$$RMSE({\delta }_{\theta })=\sqrt{\frac{1}{I}\sum _{i=1}^I{({\widehat{\delta }}_{\theta }-{\delta }_{\theta })}^2}$$

(52)

where I is the number of Monte Carlo trails. The BCRB for either ${\mathbf{\delta }}_{\theta }$ or ${\mathbf{\delta }}_{\phi }$ is a function of two variables, i.e., $\theta$ and $\phi$. Since $\theta \in [0,\frac{\pi }{2}]$ and $\phi \in [0,\pi ]$, the overall BCRB for them is shown in Figure 8 with SNR = 20 dB. From these two figures we can find that the two BCRBs are insensitive to $\phi$. However, for $\theta$, the corresponding BCRB increases with $\theta$ increasing. When $\theta =\frac{\pi }{2}$, the BCRB is approaching infinity, which means that the imaging performance is bad when $\theta$ tends to $\frac{\pi }{2}$. As for BCRB of ${\mathbf{\delta }}_{\phi }$, it approaches infinity when $\theta$ tends to zero. The reason is that the information of $\phi$ is lost since sin$\theta \to 0$ according to Equation (12).

To verify the theory in Section 3.2, the RMSE of proposed algorithm and square root of the BCRB are compared in this part. Only one point scatterer located at $(2500$ m, $\frac{\pi }{4},\frac{\pi }{4})$ is considered. The SNR increases from 0 dB to 40 dB with a stepsize of 5dB. The BCRB is calculated according to Equation (49). The RMSE is calculated using 100 Monte Carlo trails. The results are presented in Figure 9. The green stars in the plots mean that the sparse recovery failed because of the noise. With the increase of SNR, the estimation accuracy increases and the RMSE of proposed algorithm is approaching the BCRB. Thus, the proposed off-grid bias estimator is proven to be effective.

5. Conclusions

In order to solve the problems of off-grid target 3D imaging, a novel sparse Bayesian learning-based algorithm using sparse antenna array is proposed in this paper. The main impact of the off-grid problem is that it will lead to energy leakage, which will finally spoil the reconstruction results. In consideration of the characteristics of the off-grid target, a sparse Bayesian learning-based imaging algorithm is proposed to estimate not only the RCS but also the off-grid biases simultaneously. A three-stage hierarchical sparse prior is introduced and the BCRB for off-grid bias is presented, providing a lower bound for methods based on the MAP criterion. Quantitative analyses are provided to compare the 3D imaging performances of the proposed method with other state of art methods. The results show that the imaging performance of proposed method is better than other popular methods, i.e., a higher reconstruction precision and a better resolution. The effectiveness of the algorithm is verified by comparing the RMSE with its corresponding BCRB. These results all directly or indirectly validate the feasibility of proposed method. The proposed method can reconstruct the 3D model accurately and different components of the target can be recognized by the result. These findings contribute in several ways to our understanding of off-grid air target 3D reconstruction and identification, and also provide a reference for further applications of aircraft 3D imaging based on sparse antenna array.