## 1.

## Introduction

Radar coincidence imaging (RCI), motivated by classical coincidence imaging in optical systems, is a novel staring imaging technique.^{1}2.^{–}^{3} RCI can obtain focused high-resolution image without the limitation of target relative motion and operate under the observing geometry of forward-looking/staring, with significant potentials for resolution enhancement, interference, and jamming suppression. In RCI, the temporal-spatial stochastic waveforms are transmitted, thus the spatial variety of wavefront is increased, so the super-resolution within a beam emerges.

In RCI, sparse recovery is widely used as the scatterers of target, which are often distributed sparsely in many radar imaging applications. Solving the linear inverse problem with a sparsity constraint depends on the perfect prior knowledge of the system. However, phase error among the transmitter–receiver pairs exists generally since perfect synchronization is impossible for a multitransmitter configuration in RCI,^{4} because of plenty of factors, such as isolated local oscillators. Phase error results in the dictionary mismatch and induces the performance to degrade considerably, as the imaging performance depends on presetting an appropriate sparsifying dictionary based on an accurate prior known model.

Various studies have been presented on phase error. The Cramer–Rao bound for MIMO radar target localization with phase errors has been derived.^{5}^{,}^{6} To compensate the phase error, several eigenstructure-based methods are proposed.^{7}8.9.10.^{–}^{11} These methods are less sensitive to phase error but lack adaptation to demanding scenarios with low signal-to-noise ratio (SNR), limited snapshots, and spatially adjacent sources.^{12} Recently, sparse recovery and compressive sensing^{13} are introduced into signal processing by exploiting the sparsity. A sparsity-driven iterative method for joint synthetic aperture radar (SAR) imaging and phase error correction in a nonquadratic regularization-based framework is proposed in Ref. 14. The method cycles through steps of sparse recovery and phase error compensation, and the idea is widely used in SAR/inverse synthetic aperture radar (ISAR) imaging and MIMO radar imaging with phase error. Inspired by the idea of sparsity-driven iterative method introduced in Ref. 14, a sparse autocalibration imaging (SACI) method is presented to solve the RCI with gain-phase error.^{3} From the Bayesian statistics perspective, sparse imaging via expectation maximization algorithm and sparse self-calibration method via iterative minimization (SSCIM) algorithm are proposed, respectively, to alleviate the influence of phase synchronization mismatch by exploiting the maximum *a posterior* (MAP) criterion.^{15}^{,}^{16} Likewise, MAP estimator is also used into ISAR imaging by exploiting the sparseness prior of ISAR image,^{17}18.^{–}^{19} while the phase error is corrected via modified quasi-Newton algorithm. To realize array calibration and direction-of-arrival estimation, a unified framework based on sparse Bayesian learning (SBL) is formulated and a sparse Bayesian array calibration method is then proposed in Ref. 12. Using variational Bayesian inference (VBI), an array autocalibration SBL algorithm in the full conjugate Bayesian framework is proposed to achieve DOA estimation with gain/phase errors in Ref. 20. In Ref. 21, an autocalibration algorithm via variational Bayesian expectation maximization is presented to solve the problem of multiplicative perturbation.

In this paper, we focus on the sparsity-driven RCI with phase error and propose a self-calibration variational message passing (SC-VMP) algorithm in SBL framework. The merit of SBL is its flexibility in modeling sparse signals that can not only promote the sparsity but also exploit the possible structure of the signal to be recovered.^{22} As exact Bayesian inference is typically intractable, approximations are needed, such as evidence procedure^{23} and VMP.^{24} Using Bayesian hierarchical prior, the scattering coefficient is assigned an appropriate prior, i.e., three-level (3-L) hierarchical Gaussian-Gamma-Gamma (G-Ga-Ga) prior. Then, we propose a self-calibration imaging algorithm for joint imaging and phase error calibration in SBL framework. The algorithm involves an iterative method, which cycles through steps of target reconstruction and phase error estimation, where VMP and Newton’s method are adopted, respectively. Numerical simulations show that the algorithm realizes the imaging robustly and achieves both high resolution and outstanding imaging quality in the presence of phase error and is also simple to implement without changing the algorithm parameters.

The rest of the paper is organized as follows. In Sec. 2, the RCI model with phase error in the range-azimuth space is presented. Section 3 presents the SC-VMP algorithm in detail. In Sec. 4, the performance of the proposed method is verified by numerical examples. Finally, Sec. 5 concludes the paper.

## 2.

## Radar Coincidence Imaging Model with Phase Error

In this paper, RCI can be realized by a multitransmitter configuration to transmit time-independent and group-orthogonal waveforms.^{1} Consider a RCI system with $M$ transmitters and one receiver, each transmitter emits an independent stochastic waveform. The RCI geometry is illustrated in Fig. 1. The imaging plane is a range-azimuth space. In sparsity-driven RCI, the continuous imaging plane is discretized to generate $K$ small grid-cells with uniform size and shape. Denote ${\beta}_{k}$ by the scattering coefficient of the scattering center exactly located at the prediscretized $k$’th grid-cell center, i.e., ${\mathbf{r}}_{k}$ and ${\beta}_{k}=0$ for the grid-cell without scattering center.

As the backscattering of a radar target in the high-frequency region can be approximated as coming from a few dominant scattering centers,^{25} the target is assumed to be composed of a very limited amount of strong scattering centers, the number of which is much smaller than that of grid cells in the image plane.^{19} Thus, the RCI image is spatially sparse, which can be exploited to achieve super-resolution, denoising, and feature extraction.^{26}

The echo is a linear combination of all the scatterers’ reflected waveforms from all the transmitters. Considering the phase error, the echo at the receiver can be expressed as

## (1)

$$y(t)=\sum _{k\in S}\sum _{m=1}^{M}{\beta}_{k}{e}^{j{\varphi}_{m}}S{t}_{m}(t-{\tau}_{m}^{k})+w(t),$$^{1}which is simply structured as

Thus, the echo can be expressed as the superposition of the detecting signals, i.e., $y(t)=\sum _{k\in S}{\beta}_{k}S(t,{\mathbf{r}}_{k})$. After sampling the echo, the imaging equation can be given as follows:

## (3)

$$\mathbf{y}=\mathbf{S}\xb7\mathbf{\beta}+\mathbf{w}\phantom{\rule{0ex}{0ex}}\left[\begin{array}{c}y({t}_{1})\\ y({t}_{2})\\ \vdots \\ y({t}_{N})\end{array}\right]=\left[\begin{array}{cccc}S({t}_{1},{\mathbf{r}}_{1})& S({t}_{1},{\mathbf{r}}_{2})& \cdots & S({t}_{1},{\mathbf{r}}_{K})\\ S({t}_{2},{\mathbf{r}}_{1})& S({t}_{2},{\mathbf{r}}_{2})& \cdots & S({t}_{2},{\mathbf{r}}_{K})\\ \vdots & \vdots & \cdots & \vdots \\ S({t}_{N},{\mathbf{r}}_{1})& S({t}_{N},{\mathbf{r}}_{2})& \cdots & S({t}_{N},{\mathbf{r}}_{K})\end{array}\right]\xb7\left[\begin{array}{c}{\beta}_{1}\\ {\beta}_{2}\\ \vdots \\ {\beta}_{K}\end{array}\right]+\left[\begin{array}{c}w({t}_{1})\\ w({t}_{2})\\ \vdots \\ w({t}_{N})\end{array}\right],$$Under the assumption of sparse prior of the target, the imaging model shown in Eq. (3) is a typical linear model used in most sparse recovery applications. However, the dictionary $\mathbf{S}$ obtained from practical applications is inevitably disturbed by perturbations, which cannot be known accurately in practice, e.g., phase error. Then, $\mathbf{S}$ can be rewritten as $\mathbf{S}(\mathbf{\varphi})$ involving the phase error $\mathbf{\varphi}$, where $\mathbf{\varphi}={[{\varphi}_{1},\cdots ,{\varphi}_{M}]}^{T}$. Then, Eq. (3) can be rewritten as

Involving the generally unknown parameter $\mathbf{\varphi}$, $\mathbf{S}(\mathbf{\varphi})$ is a parameterized unknown dictionary and then $\mathbf{\beta}$ could not be reconstructed directly based on the conventional sparse recovery algorithms. Therefore, a self-calibration imaging algorithm, which is presented in the following section, is proposed to solve the problem.

## 3.

## Self-Calibration Variational Message Passing

Conventional sparsity-driven radar imaging methods assume that the model contains no error and the dictionary is precisely known, which is not realizable in most application scenarios. In the presence of phase error, the structure of the dictionary is destroyed, which leads to the direct use of sparse recovery methods failure. In this section, the SC-VMP algorithm calibrating the phase error is presented. This algorithm is formulated in a Bayesian hierarchical prior model, and works by jointly reconstructing the target and estimating the phase error, where VMP and Newton’s method are adopted, respectively.

## 3.1.

### Hierarchical Prior Model

The hierarchical representation of sparsity-inducing prior allows us to choose simple and analytically tractable probability density functions (PDFs), and results in the construction of efficient yet computationally tractable iterative inference algorithms with analytical derivation of the inference expressions.^{27}

The graphical representation of RCI model priors is shown in Fig. 2. To recover $\mathbf{\beta}$ in Eq. (4), we model it as hidden (unobserved) variable and assign it a 3-L hierarchical G-Ga-Ga prior typically to induce sparsity:

## (5)

$$p(\mathbf{\beta}|\mathbf{\alpha})=\mathcal{CN}(\mathbf{\beta}|\mathbf{0},\mathbf{\Lambda})=\frac{1}{{\pi}^{K}\mathrm{det}(\mathbf{\Lambda})}\text{\hspace{0.17em}}\mathrm{exp}\{-{\mathbf{\beta}}^{H}{\mathbf{\Lambda}}^{-1}\mathbf{\beta}\},$$## (6)

$$p(\mathbf{\alpha}|\eta ;\u03f5)=\prod _{k=1}^{K}\mathrm{\Gamma}({\mathbf{\alpha}}_{k}|\u03f5,\eta ),$$^{27}

Compared with the two-layer (2-L) prior model, the 3-L prior model considers $\eta $ as random variable specified by a distribution and incorporates it into the inference framework. This leads to additional degrees of freedom in controlling the sparsity property of the resulting inference scheme.^{27}

In addition, the measurement noise is assumed to be complex Gaussian with zero-mean and variance ${\alpha}_{0}^{-1}$, to allow conjugate-exponential analysis. Then, the following 2-L hierarchical Gaussian-Gamma (G-Ga) prior is assigned:

## (8)

$$p(\mathbf{w}|{\alpha}_{0})=\mathcal{CN}(\mathbf{w}|\mathbf{0},{\alpha}_{0}^{-1}\mathbf{I})={(\pi {\alpha}_{0}^{-1})}^{-K}\text{\hspace{0.17em}}\mathrm{exp}\{-{\alpha}_{0}{\Vert \mathbf{w}\Vert}_{2}^{2}\},$$Encoded by the graphical model shown in Fig. 2, the joint PDF of the model [Eq. (4)] is of the form:

## 3.2.

### Variational Message Passing

In this part, we present the VMP algorithm to estimate $\mathbf{\beta}$. Denote $\mathbf{\Omega}=\{\mathbf{\beta},\mathbf{\alpha},\eta ,{\alpha}_{0}\}$ as the set of hidden variables. In Bayesian inference framework, the exact posterior distribution $p(\mathbf{\Omega}|\mathbf{y};\mathbf{\varphi})=p(\mathbf{y},\mathbf{\Omega};\mathbf{\varphi})/p(\mathbf{y})$ is intractable since $p(\mathbf{y})$ cannot be expressed explicitly. VBI^{28} could be used to find a tractable distribution $q(\mathbf{\Omega})$ that closely approximates the true posterior distribution $p(\mathbf{\Omega}|\mathbf{y};\mathbf{\varphi})$ by minimizing the Kullback–Leibler divergence (KLD) between them.^{29} A structured mean field approximation over $p(\mathbf{\Omega}|\mathbf{y};\mathbf{\varphi})$ is further assumed as

^{24}in this paper. VMP is an iterative scheme that uses a message passing procedure on a graphical model and attempts to compute the auxiliary PDF. The variational distributions $\{q(\mathbf{\beta}),q(\mathbf{\alpha}),q(\eta ),q({\alpha}_{0})\}$ are iteratively updated in VMP procedure to monotonically decrease the KLD and thus, the convergence is guaranteed.

The updates of $\{q(\mathbf{\beta}),q(\mathbf{\alpha}),q(\eta ),q({\alpha}_{0})\}$ are similar to those in Ref. 27, the derivation detail is shown in Appendix A. It can be concluded that the posterior distributions of $\mathbf{\beta}$ and ${\alpha}_{k}$ are complex Gaussian distribution and generalized inverse Gaussian (GIG) distribution, respectively, while both $\eta $ and ${\alpha}_{0}$ obey a Gamma distribution.

Based on the derivation of $\{q(\mathbf{\beta}),q(\mathbf{\alpha}),q(\eta ),q({\alpha}_{0})\}$, the optimal approximated distribution $q(\mathbf{\Omega})$ can be obtained by iteratively calculating Eqs. (19), (20), (23), (25), and (26) until convergence.

## 3.3.

### Phase Error Estimation

In the presence of phase error, the convergence of SC-VMP is not a direct result of VMP. Consider the phase error $\mathbf{\varphi}$ as an unknown deterministic parameter as it is not varying generally during the entire coherent processing interval.^{6} Then, the resulting algorithm can be interpreted as a variational EM algorithm.^{28}

In VMP framework, the phase error can be estimated from the maximum likelihood:

## (12)

$$\stackrel{\u2322}{\mathbf{\varphi}}=\underset{\mathbf{\varphi}}{\mathrm{argmax}}{\u27e8\mathrm{ln}\text{\hspace{0.17em}}p(\mathbf{y},\mathbf{\Omega};\mathbf{\varphi})\u27e9}_{q(\mathbf{\beta})q({\alpha}_{0})q(\mathbf{\alpha})q(\eta )},$$## (13)

$$\stackrel{\u2322}{\mathbf{\varphi}}=\underset{\mathbf{\varphi}}{\mathrm{argmin}}\{{\Vert \mathbf{y}-\mathbf{S}(\mathbf{\varphi})\mathbf{\mu}\Vert}_{2}^{2}+\mathrm{tr}({[\mathbf{S}(\mathbf{\varphi})]}^{H}\mathbf{S}(\mathbf{\varphi})\mathbf{\Sigma})\}.$$^{30}which reveals the behavior of phase error with no approximation being required. The updated $\stackrel{\u2322}{\mathbf{\varphi}}$ estimate is

## (14)

$$\stackrel{\u2322}{\mathbf{\varphi}}=\mathbf{\varphi}-{[{\nabla}^{2}f(\mathbf{\varphi})]}^{-1}[\nabla f(\mathbf{\varphi})],$$## 3.4.

### Algorithm Description

In SC-VMP procedure, each iteration consists of matrix inversion and matrix-vector multiplication in Eqs. (19) and (20), which leads to expensive computational burden in case of large number of grid cells. Besides, VMP suffers from low rate of convergence. Hence, a fast implementing approach is developed.

In conventional SBL, many of the prior precisions of coefficients, i.e., ${\u27e8{\alpha}_{k}^{-1}\u27e9}_{q(\mathbf{\alpha})}$, typically take on quite large values upon convergence, which implies that the corresponding coefficients are quite small and in turn the contributions of the corresponding bases could be negligible. Thus, the corresponding bases could be removed from the model to realize sparsity when ${\u27e8{\alpha}_{k}^{-1}\u27e9}_{q(\mathbf{\alpha})}$ exceeds a certain large threshold ${\alpha}_{\mathrm{th}}$. Then, the computational complexity of SC-VMP can be drastically reduced. Therefore, to speed up the algorithm, the pruning of the current basis set ${\mathrm{\Theta}}^{i}$ can be achieved via

## (15)

$${\mathrm{\Theta}}^{i+1}=\{k|{\u27e8{\alpha}_{k}^{-1}\u27e9}_{q(\mathbf{\alpha})}<{\alpha}_{\mathrm{th}},k\in {\mathrm{\Theta}}^{i}\},$$In addition, we terminate the algorithm if ${\Vert {\mathbf{\mu}}^{i+1}-{\mathbf{\mu}}^{i}\Vert}_{2}/{\Vert {\mathbf{\mu}}^{i}\Vert}_{2}<\gamma $ or the maximum number of iterations ${I}_{\mathrm{max}}$ is reached, where $\gamma $ is a user-defined tolerance and the superscript $i$ refers to the iteration index.

Based on the above discussions, the procedure of SC-VMP algorithm is given in Algorithm 1.

## Algorithm 1

SC-VMP algorithm.

Input: $\mathbf{y}$, $\mathbf{S}(\mathbf{\varphi}=\mathbf{0})$, ${I}_{\mathrm{max}}$, ${\mathbf{\alpha}}_{\mathrm{th}}$, $\gamma $, $\{a,b,c,d\}$ |

Initialization: $i=0$, $\mathbf{\varphi}=\mathbf{0}$, ${\mathrm{\Theta}}^{0}=\{k|k=1,\cdots ,K\}$ |

while not converged do |

1): VMP procedure: Update $\mathbf{\Sigma}$, $\mathbf{\mu}$, ${\u27e8{\alpha}_{k}^{-1}\u27e9}_{q(\mathbf{\alpha})}$, ${\u27e8{\alpha}_{k}\u27e9}_{q(\mathbf{\alpha})}$, ${\u27e8\eta \u27e9}_{q(\eta )}$, ${\u27e8{\alpha}_{0}\u27e9}_{q({\mathbf{\alpha}}_{0})}$ form in Eqs. (19), (20), (23), (25), and (26), where $k\in {\mathrm{\Theta}}^{i}$ |

2): Newton’s method: Estimate $\mathbf{\varphi}$ from Eq. (14), and update $\mathbf{S}(\mathbf{\varphi})$ |

3): Prune the bases: ${\mathrm{\Theta}}^{i+1}=\{k|{\u27e8{\alpha}_{k}^{-1}\u27e9}_{q(\mathbf{\alpha})}<{\alpha}_{\mathrm{th}},k\in {\mathrm{\Theta}}^{i}\}$ |

4): $i=i+1$ and check for convergence: ${\Vert {\mathbf{\mu}}^{i}-{\mathbf{\mu}}^{i-1}\Vert}_{2}^{2}/{\Vert {\mathbf{\mu}}^{i}\Vert}_{2}^{2}<\gamma $ or $i={I}_{\mathrm{max}}$ |

end while |

Output: reconstructed scattering coefficient vector $\mathbf{\beta \u0361}=\mathbf{\mu}$ |

## 3.5.

### Discussion

After investigating the SC-VMP algorithm in depth, we provide some discussions to gain insight into the algorithm.

1. Compared with related work: Similar with most of the sparsity-based methods for radar imaging with phase error presented in Refs. 3, 12, 1415.16.17.18.19.20.–21, the proposed SC-VMP algorithm also employs the alternating strategy to update the sparse scattering coefficients and phase error iteratively. The methods listed in Refs. 3, 1415.16.17.18.–19 can be summarized as $l1$-based regularization method, where a sparsity-inducing Laplace prior is directly imposed on the signal, then the sparse solution is exploited from the MAP estimation, which corresponds to the point estimate of the sparse scattering coefficient. Comparably, the SC-VMP utilizes the hierarchical modeling procedure to encode signal sparsity and achieve better sparse solutions, which could obtain the approximate posterior distribution and is regarded as a full Bayesian method. During the procedure, the statistical information is used to enhance the estimation performance and avoid converging to a shallow local minimum, due to the utilization of higher order statistical information (such as the estimation covariance matrix) in the full Bayesian framework. Besides, the proposed algorithm can properly utilize uncertainty information during iterations to ameliorate the error propagation problem, which means the estimation error of the sparse signal would degrade the estimation accuracy of phase error during iterations. Inevitably, the error propagation problem exists generally in these methods as the point estimate of the sparse scattering coefficient and phase error are updated alternately. Furthermore, different from the methods in Refs. 14, 20, and 21, we resort to Newton’s method to update the phase error, as it is not tractable to obtain the closed-from expression for updating the parameter when solving the nonlinear least-squares problem expressed in Eq. (13).

2. Convergence: Since the SC-VMP algorithm can be interpreted as a variational EM algorithm, the update of $\mathbf{\beta}$ and $\mathbf{\varphi}$ will monotonically decrease the KLD and the negative expected log-likelihood function, respectively, until convergence. Thus, the (marginal) likelihood monotonically increases during the iterations and the convergence is guaranteed.

^{28}Moreover, the possibility of converging to a local minimum is reduced due to the utilization of higher order statistical information.3. Computational complexity: The main computational burden at each iteration of SC-VMP procedure comes from the update of $\mathbf{\beta}$ in VMP procedure and the update of $\mathbf{\varphi}$ in Newton’s method. Two main time-consuming operations are the matrix inversion and matrix-vector multiplication whose computational complexity are $o({K}^{3})$ and $o({K}^{2})$, respectively. This can incur a high computational cost when $K$ is large. Fortunately, the grid pruning could reduce the computational burden and improve the convergence speed, as the remaining grid number after pruning decreases fast. Additionally, the computational cost can be also controlled by restricting the maximum number of iteration.

## 4.

## Numerical Simulations

In this section, simulations are carried out to verify the SC-VMP algorithm. An X-band RCI radar system with carrier frequency of 10 GHz is considered. The transmitters are configured as a uniform linear array with $M=16$ and the interelement spacing $d=0.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The transmitters emit independent frequency-hopping (FH) waveforms with the bandwidth of 500 MHz. A range-azimuth imaging plane, covering $40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}\times 40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, is discretized to $40\times 40$ grid cells. Further, there are supposed to be 17 point scatterers in the imaging plane.

## 4.1.

### Illustrative Example

In this part, the proposed SC-VMP algorithm is used to reconstruct the target image. For comparison, we also give the results obtained by four sparse-based imaging methods, i.e., VMP-3L, SSCIM,^{16} sparsity-driven autofocus (SDA) method^{14} and SACI^{3} method. VMP-3L is presented in Sec. 3. During the numerical experiments, the phase error randomly varied at $[-45\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg},45\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}]$, and the SNR is fixed at 15 dB. The experiments are performed on a computer with Intel Core CPU i3-4130 at 3.4 GHz and 4 GB of memory.

Figure 3 shows the true image and target images reconstructed by VMP-3L, SSCIM, SDA, SACI, and SC-VMP. It can be concluded that our proposed SC-VMP algorithm achieves superior imaging performance over the other four algorithms tested above. As shown in Fig. 6(f), the target is reconstructed perfectly without any spurious scatterers. Comparably, the other four algorithms suffer from imperfect reconstruction performance, and the reconstructed images are defocused with some spurious scatterers and the strength of scattering centers is not reconstructed exactly, although the target profiles are clear, which makes the images recognizable. Comparing the list algorithms, the difference is mainly focused on the method to reconstruct the scattering coefficients (see Sec. 3.5 for details). For SSCIM, SDA, and SACI, they suffer from the effect of phase error and the signal energy spills over the imaging plane because of the dictionary mismatch caused by phase error, as shown in Figs. 3(c)–3(e), while our proposed SC-VMP utilizes the hierarchical modeling procedure and is considered as a full Bayesian method. More importantly, the method is less sensitive to the phase error, as shown in Fig. 3(b). In addition, the phase error is estimated and compensated perfectly by the SC-VMP algorithm, as shown in Fig. 4. Then, we could conclude that exploiting the phase error improves the performance of SC-VMP. Consequently, benefiting from full utilization of the sparsity prior and phase error calibration, the image quality shown in Fig. 3(f) is improved significantly, compared with the images reconstructed by other algorithms.

To evaluate the numerical complexity of the proposed algorithm, we record the runtime of the five algorithms during the simulations, and the result is shown in Table 1. We can see that SSCIM is time consuming for their heavy computational complexity, because the scattering coefficients are reconstructed by basis pursuit denoising algorithm. As orthogonal matching pursuit is employed in the iterations, SACI shows the best computational efficiency. The other three algorithms, i.e., VMP-3L, SDA, and SC-VMP, are all Bayesian-based algorithms and their computational complexity lies about the same level. Thus, it can be observed that our proposed algorithm does not show superior performance over other algorithms from the perspective of computational complexity. However, the SC-VMP algorithm can be improved to make it more computationally effective, which will be investigated in our future work.

## Table 1

Runtime of the five algorithms.

Algorithm | VMP-3L | SSCIM | SDA | SACI | SC-VMP |
---|---|---|---|---|---|

Runtime (s) | 48.3941 | 411.2841 | 42.6616 | 6.0714 | 62.3943 |

## 4.2.

### Performance of Self-Calibration Variational Message Passing Algorithm

To evaluate the performance of the proposed algorithm, we introduce two criterions, i.e., the relative imaging error (RIE) and probability of successful imaging (PSI). RIE is defined as $20\text{\hspace{0.17em}}{\mathrm{log}}_{10}({\Vert \stackrel{\u2322}{\mathbf{\beta}}-\mathbf{\beta}\Vert}_{2}/{\Vert \mathbf{\beta}\Vert}_{2})$, where $\mathbf{\beta}$ and $\stackrel{\u2322}{\mathbf{\beta}}$ are the true and estimated scattering coefficient vector, respectively. To define PSI, a metric $\mathrm{\Delta}=\mathrm{min}{(\stackrel{\u2322}{\mathbf{\beta}})}_{\mathrm{\Lambda}}/\mathrm{max}{(\stackrel{\u2322}{\mathbf{\beta}})}_{\mathrm{\Lambda}}$^{31} is given first, where ${(\stackrel{\u2322}{\mathbf{\beta}})}_{\mathrm{\Lambda}}$ contains the values that $\stackrel{\u2322}{\mathbf{\beta}}$ carries at the correct basis set $\mathrm{\Lambda}$ and ${(\overline{\stackrel{\u2322}{\mathbf{\beta}}})}_{\mathrm{\Lambda}}$ takes 0 at $\mathrm{\Lambda}$ and takes the same values as $\stackrel{\u2322}{\mathbf{\beta}}$ at every other indices. $\mathrm{\Delta}>1$ can guarantee exact estimation of the scatter’s position, which means a successful imaging trial. Hence, PSI is defined as the percentage of successful imaging trials.

First, we quantitatively analyze the performance of the aforementioned five algorithms versus SNR, and the result is shown in Fig. 5. It can be seen that the RIEs decrease quickly with the increase in SNR, which implies that the image quality is improved as the SNR increases and the algorithms are sensitive to noise, especially in low SNR regimes. Comparing with other algorithms, the RIE obtained by SC-VMP is much lower when $\mathrm{SNR}\ge 0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$. Besides, it is worth noting that in Fig. 5(b), the SC-VMP outperforms other reported algorithms with a higher PSI, especially $\mathrm{SNR}>5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$ can guarantee the target reconstruction for SC-VMP ($\mathrm{PSI}\ge 96\%$ in this case). Thus, SC-VMP exhibits significant performance gains over other algorithms in terms of both metrics, by compensating the phase error, whereas its performance is sensitive to SNR. More precisely, both the scattering coefficients estimation and phase error estimation are seriously affected by noise. The scattering coefficients are reconstructed by VMP-3L whose performance is shown in Fig. 5, while the performance of phase error estimation is given by Fig. 6, which shows the normalized mean square error (NMSE) of phase error, NMSE is defined as $20\text{\hspace{0.17em}}{\mathrm{log}}_{10}({\Vert \stackrel{\u2322}{\mathbf{\varphi}}-\Vert}_{2}/{\Vert \mathbf{\varphi}\Vert}_{2})$. We can see that low SNR results in high estimation error, which means that the phase error cannot be compensated accurately, and the target reconstruction would be degraded accordingly.

Next, we investigate the robustness of SC-VMP against the hopping frequency number and the scope of phase error. In this paper, the sampling interval is setup as the FH duration, which means hopping frequency number is equal to the number of samples $N$. The scope of phase error is $[-{\varphi}_{\mathrm{max}},{\varphi}_{\mathrm{max}}]$, where ${\varphi}_{\mathrm{max}}={[15k]}^{\xb0}$ $(k=1,\cdots ,9)$ is the maximal value of phase error. For each scope of phase error, 100 independent trials are conducted and ${\varphi}_{m}$ is randomly distributed within the scope.

As is clear from Fig. 7, the performance improves as the hopping frequency number increases, for both high imaging quality and high PSI. More independent hopping frequencies could provide more information of resolvability for RCI, which would also lead to the robustness of imaging. Besides, the algorithm is sensitive to the scope of phase error. When the phase error increases, the convergence of SC-VMP degrades, then the phase error could not be estimated exactly and calibrated perfectly. Thus, the imaging performance degrades severely, especially when ${\varphi}_{\mathrm{max}}\ge 90\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$.

## 5.

## Conclusion

RCI is a high-resolution imaging technique, while its performance degrades seriously in the presence of phase error. The SC-VMP algorithm is proposed in this paper to realize the target reconstruction and phase error calibration simultaneously in SBL framework. Bayesian hierarchical prior is used in the proposed algorithm, where the scattering coefficient is modeled statistically and estimated iteratively to achieve sparsity. As an iterative algorithm, SC-VMP cycles through steps of target reconstruction and phase error estimation, where VMP and Newton’s method are adopted, respectively. The proposed algorithm can estimate the phase error accurately and improve the reconstruction performance significantly. Besides, it requires no prior information and provides superior and robust imaging performance for various imaging scenes and system parameters, which shows the potential for the algorithm to be applied in practical RCI or other imaging radar systems.

## Appendices

## Appendix A:

### Implementation of Variational Message Passing Procedure

During the VMP procedure, the hidden variables $\mathbf{\Omega}=\{\mathbf{\beta},\mathbf{\alpha},\eta ,{\alpha}_{0}\}$ can be updated individually.

1. Update of $q(\mathbf{\beta})$. It is obvious to conclude from Eq. (8) that the echo $\mathbf{y}$ obeys a complex Gaussian distribution and the likelihood function of the observation can be formulated as

## (16)

$$p[p(\mathbf{y}|\mathbf{\beta},{\alpha}_{o};\mathbf{\varphi})]=\mathcal{C}\mathcal{N}(\mathbf{y}|\mathbf{S}\mathbf{\beta},{\alpha}_{0}^{-1}\mathbf{I}).$$According to Ref. 27, the approximated posterior of $\mathbf{\beta}$ can be expressed as

Substituting Eqs. (5) and (16) into Eq. (17), we can conclude that $\mathbf{\beta}$ obeys a complex Gaussian distribution: where the mean $\mathbf{\mu}$ and covariance $\mathbf{\Sigma}$ are given by## (17)

$$q(\mathbf{\beta})\propto \mathrm{exp}\{{\u27e8\mathrm{ln}\text{\hspace{0.17em}}p(\mathbf{y}|\mathbf{\beta},{\alpha}_{0};\mathbf{\varphi})\u27e9}_{q({\mathbf{\alpha}}_{0})}{\u27e8p(\mathbf{\beta}|\mathbf{\alpha})\u27e9}_{q(\mathbf{\alpha})}\}.$$2. Update of $q(\mathbf{\alpha})$. The approximated posterior of $\mathbf{\alpha}$ can be given by

Then, we substitute Eqs. (5) and (6) into Eq. (21) and obtain $q(\mathbf{\alpha})$ as## (21)

$$q(\mathbf{\alpha})\propto \mathrm{exp}\{{\u27e8p(\mathbf{\beta}|\mathbf{\alpha})\u27e9}_{q(\mathbf{\beta})}{\u27e8\mathrm{ln}\text{\hspace{0.17em}}p(\mathbf{\alpha}|\eta )\u27e9}_{q(\eta )}\}.$$where ${\u27e8{|{\beta}_{k}|}^{2}\u27e9}_{q(\mathbf{\beta})}={|{\mu}_{k}|}^{2}+{\mathrm{\Sigma}}_{kk}$. Thus, $q(\mathbf{\alpha})$ is the product of GIG PDFs with order $p=\u03f5-1$. The moments of GIG distribution are given in closed form for any $n\in \mathbb{R}$:## (22)

$$q(\mathbf{\alpha})\propto \prod _{k=1}^{K}{\alpha}_{k}^{\u03f5-2}\text{\hspace{0.17em}}\mathrm{exp}(-{\alpha}_{k}^{-1}\u27e8{|{\beta}_{k}|}^{2}{\u27e9}_{q(\mathbf{\beta})}-{\alpha}_{k}{\u27e8\eta \u27e9}_{q(\eta )}),$$where ${\kappa}_{\upsilon}(\xb7)$ is the modified Bessel function of the second kind with order $\upsilon $. The case of $n=-1$ in Eq. (23) gives the evaluation of ${\u27e8{\mathbf{\Lambda}}^{-1}\u27e9}_{q(\mathbf{\alpha})}$ used in Eq. (19), and the case of $n=1$ gives the calculation of ${\u27e8{\alpha}_{k}\u27e9}_{q(\mathbf{\alpha})}$ used in Eq. (25).## (23)

$${\u27e8{\alpha}_{k}^{n}\u27e9}_{q(\mathbf{\alpha})}={\left(\frac{{\u27e8{|{\beta}_{k}|}^{2}\u27e9}_{q(\mathbf{\beta})}}{{\u27e8\eta \u27e9}_{q(\eta )}}\right)}^{\frac{n}{2}}\frac{{\kappa}_{p+n}\left(2\sqrt{{\u27e8\eta \u27e9}_{q(\eta )}{\u27e8{|{\beta}_{k}|}^{2}\u27e9}_{q(\mathbf{\beta})}}\right)}{{\kappa}_{p}\left(2\sqrt{{\u27e8\eta \u27e9}_{q(\eta )}{\u27e8{|{\beta}_{k}|}^{2}\u27e9}_{q(\mathbf{\beta})}}\right)},$$3. Update of $q(\eta )$. The approximated posterior distribution of $\eta $ is obtained by

## (24)

$$q(\mathbf{\alpha})\propto \mathrm{exp}\{{\u27e8p(\mathbf{\alpha}|\eta )\u27e9}_{q(\eta )}p(\eta ;a,b)\}.$$It can be seen that the posterior distribution of $\eta $ is a Gamma distribution: $q(\eta )=\mathrm{\Gamma}(\eta |K\u03f5+a,\sum _{k=1}^{K}{\u27e8{\alpha}_{k}\u27e9}_{q(\mathbf{\alpha})}+b)$, while the mean of $\eta $ is

4. Update of $q({\alpha}_{0})$. It can be shown that $q({\alpha}_{0})=\mathrm{\Gamma}({\alpha}_{0}|N+c,{\u27e8{\Vert \mathbf{y}-\mathbf{S}\xb7\mathbf{\beta}\Vert}_{2}^{2}\u27e9}_{q(\mathbf{\beta})}+d)$. The first moment of ${\alpha}_{0}$ is

## Appendix B:

### Derivation of Estimating Phase Error

In this part, we derive the gradient $\nabla f(\mathbf{\varphi})$ and Hessian ${\nabla}^{2}f(\mathbf{\varphi})$, which are used in Eq. (14) to estimating the phase error $\mathbf{\varphi}$. For convenience of derivation, $\nabla f(\mathbf{\varphi})$ is decomposed by two parts, i.e., $f(\mathbf{\varphi})={f}_{1}(\mathbf{\varphi})+{f}_{2}(\mathbf{\varphi})$, where ${f}_{1}(\mathbf{\varphi})={\Vert \mathbf{y}-\mathbf{S}(\mathbf{\varphi})\mathbf{\mu}\Vert}_{2}^{2}$ and ${f}_{2}(\mathbf{\varphi})=\mathrm{tr}\{{[\mathbf{S}(\mathbf{\varphi})]}^{H}\mathbf{S}(\mathbf{\varphi})\mathbf{\Sigma}\}$. Then, we have

## (27)

$$\nabla f(\mathbf{\varphi})=\nabla {f}_{1}(\mathbf{\varphi})+\nabla {f}_{2}(\mathbf{\varphi}),$$## (28)

$${\nabla}^{2}f(\mathbf{\varphi})={\nabla}^{2}{f}_{1}(\mathbf{\varphi})+{\nabla}^{2}{f}_{2}(\mathbf{\varphi}).$$After derivation and simplification, the gradient and Hessian matrix can be computed as

## (29)

$${[\nabla {f}_{1}(\mathbf{\varphi})]}_{m}=-2\text{\hspace{0.17em}}\mathrm{Im}\{{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})\mathbf{\mu}]}^{H}\stackrel{\u2322}{\mathbf{w}}\},$$## (30)

$${[\nabla {f}_{2}(\mathbf{\varphi})]}_{m}=\mathrm{tr}(j\{{[\mathbf{S}(\mathbf{\varphi})]}^{H}\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})-{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]}^{H}\mathbf{S}(\mathbf{\varphi})\}\mathbf{\Sigma}),$$## (31)

$${[{\nabla}^{2}{f}_{1}(\mathbf{\varphi})]}_{mm}=2\text{\hspace{0.17em}}\mathrm{Re}\{{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})\xb7\mathbf{\mu}]}^{H}\stackrel{\u2322}{\mathbf{w}}\}+2{\Vert \stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})\xb7\mathbf{\mu}\Vert}^{2},$$## (32)

$${[{\nabla}^{2}{f}_{1}(\mathbf{\varphi})]}_{ml}=2\text{\hspace{0.17em}}\mathrm{Re}\{{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})\xb7\mathbf{\mu}]}^{H}[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{l}({\mathbf{\varphi}}_{l})\xb7\mathbf{\mu}]\},\phantom{\rule[-0.0ex]{2em}{0.0ex}}m\ne l,$$## (33)

$${[{\nabla}^{2}{f}_{2}(\mathbf{\varphi})]}_{mm}=\mathrm{tr}(\{2{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]}^{H}[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]-{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]}^{H}\mathbf{S}(\mathbf{\varphi})-{[\mathbf{S}(\mathbf{\varphi})]}^{H}[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]\}\mathbf{\Sigma}),$$## (34)

$${[{\nabla}^{2}{f}_{2}(\mathbf{\varphi})]}_{ml}=\mathrm{tr}(\{{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]}^{H}[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{l}({\mathbf{\varphi}}_{l})]+{[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{l}({\mathbf{\varphi}}_{l})]}^{H}[\stackrel{\u2322}{\mathbf{S}}{\mathbf{t}}_{m}({\mathbf{\varphi}}_{m})]\}\mathbf{\Sigma}),\phantom{\rule[-0.0ex]{2em}{0.0ex}}m\ne l,$$## Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 61302149 and 61302142) and Research Fund for the Doctoral Program of Higher Education of China (20124307110013). The authors would like to thank the editors and reviewers for their insightful comments.

## References

## Biography

**Xiaoli Zhou** received his BS and MS degrees in instrument science and technology from the National University of Defense Technology, Changsha, China, in 2010 and 2012, respectively. He is currently pursuing his PhD degree with the National University of Defense Technology. His current research interests include remote sensing signal processing and coincidence imaging, as well as advanced signal processing with application to radar target imaging and identification.

**Hongqiang Wang** received his BS, MS, and PhD degrees from the National University of Defense Technology, Changsha, China, in 1993, 1999, and 2002, respectively. He is currently a professor with the School of Electronic Science and Engineering, National University of Defense Technology. He has been involved in modern radar signal processing research and development since 1996. His current research interests include automatic target recognition, radar imaging, and target tracking.

**Yongqiang Cheng** received his BS, MS, and PhD degrees in information and communication engineering from the National University of Defense Technology, Changsha, China, in 2005, 2007, and 2012, respectively. He is currently an associate professor with the National University of Defense Technology. From September 2009 to November 2010, he was a visiting research student with Melbourne Systems Laboratory, Department of Electrical and Electronic Engineering, University of Melbourne, Melbourne, Australia. His current research interests include statistical signal processing and information geometry.

**Yuliang Qin** received his BS, MS, and PhD degrees in information and communication engineering from the National University of Defense Technology, Changsha, China, in 2002, 2004, and 2008, respectively. He is currently an associate professor with the School of Electronic Science and Engineering, National University of Defense Technology. His current research interests include SAR imaging and radar signal processing.