Random Stepped Frequency ISAR 2D Joint Imaging and Autofocusing by Using 2D-AFCIFSBL

Abstract

With the increasingly complex electromagnetic environment faced by radar, random stepped frequency (RSF) has garnered widespread attention owing to its remarkable Electronic Counter-Countermeasure (ECCM) characteristic, and it has been universally applied in inverse synthetic aperture radar (ISAR) in recent years. However, if the phase error induced by the translational motion of the target in RSF ISAR is not precisely compensated, the imaging result will be defocused. To address this challenge, a novel 2D method based on sparse Bayesian learning, denoted as 2D-autofocusing complex-value inverse-free SBL (2D-AFCIFSBL), is proposed to accomplish joint ISAR imaging and autofocusing for RSF ISAR. First of all, to integrate autofocusing into the ISAR imaging process, phase error estimation is incorporated into the imaging model. Then, we increase the speed of Bayesian inference by relaxing the evidence lower bound (ELBO) to avoid matrix inversion, and we further convert the iterative process into a matrix form to improve the computational efficiency. Finally, the 2D phase error is estimated through maximum likelihood estimation (MLE) in the image reconstruction iteration. Experimental results on both simulated and measured datasets have substantiated the effectiveness and computational efficiency of the proposed 2D joint imaging and autofocusing method.

1. Introduction

Inverse synthetic aperture radar (ISAR) [1] can achieve 2D high-resolution images of noncooperative moving targets in all-day and all-weather scenarios, which plays a significant role in both civil and military fields [2]. For traditional ISAR imaging, it is usually necessary to transmit a large bandwidth signal to achieve high resolution in the range direction. However, transmitting a large bandwidth signal is not only susceptible to interference but also puts tremendous pressure on the radar system hardware. To overcome these difficulties, random stepped frequency (RSF) radar was introduced into the field of ISAR imaging [3]. Because the carrier frequency of RSF radar jumps rapidly and randomly [4], making it unpredictable, RSF radar exhibits excellent Electronic Counter-Countermeasure (ECCM) characteristics in modern complex electromagnetic environments [5]. RSF ISAR can synthesize wideband signals by transmitting a large number of narrowband pulses in a burst, thereby achieving high range resolution. RSF radar significantly reduces the instantaneous bandwidth requirements and overall cost of the radar system. Therefore, continuing to study RSF ISAR imaging is valuable. However, RSF ISAR requires a long coherent processing interval (CPI), which invalidates the assumption of a constant rotation rate of the target and further affects the quality of imaging. To mitigate this issue, sparse RSF (SRSF) technology [6] has been introduced, reducing CPI by transmitting fewer subpulses and bursts. At the same time, random stepped frequency radar can track and image multiple targets, which can result in the loss of some pulses. Additionally, in complex electromagnetic environments, some subpulses and bursts may be unavailable, making the study of SRSF ISAR imaging highly significant.

Considering the inherent sparsity of radar observation targets in space, ISAR imaging can be seen as a problem of sparse signal reconstruction. Compressive sensing (CS) [7], capable of reconstructing sparse signals from insufficient measurements, has been widely applied in the ISAR imaging field in recent years. It has been demonstrated to be more effective than many other methods [8]. Similar to traditional range-Doppler (RD) methods, the most common methods based on CS synthesize high-resolution images separately for the range and cross-range directions [6,9]. Consequently, a number of CS-based reconstruction algorithms have been employed in RSF ISAR imaging, including orthogonal matching pursuit (OMP) [10], smoothed ${\ell }_0$ norm (SL0) [11], ${\ell }_1$ norm optimization [9] and alternating direction method of multipliers (ADMM) [12]. Nevertheless, the imaging quality of these algorithms relies on the choice of model parameters and is influenced by the SNR of the echo and sparse aperture. Unlike other methods that require manual parameter adjustment, the sparse Bayesian learning (SBL) [13] algorithm produces superior reconstruction outcomes without the need for manual parameter adjustment. Many ISAR imaging methods based on the SBL algorithm have been introduced [14,15,16,17]. However, these SBL algorithms typically involve matrix inversion with high computational complexity, resulting in the computational complexity of SBL-based ISAR imaging algorithms being impractical. Zhang et al. [18] proposed a sparse Bayesian algorithm based on ADMM, treating the problem of solving the posterior expectation of sparse vectors as a convex optimization problem. ADMM is utilized to solve this convex optimization problem, thus avoiding matrix inversion. While the complexity of the algorithm is reduced, it introduces manually adjusted regularization parameters.

However, the sequential imaging method does not fully account for the two-dimensional sparsity of ISAR images, resulting in a decrease in imaging quality. To enhance imaging sparsity, a joint imaging framework is proposed, which can acquire high-resolution profiles in both range and cross-range directions simultaneously. It is common practice to vectorize a 2D ISAR image [19,20]. However, vectorized optimization encounters challenges in terms of high computational complexity and substantial memory requirements due to the large dictionary dimension. To directly process sparse matrix data, various 2D compressive sensing (2D-CS) algorithms, including 2D smoothed ${\ell }_0$ norm (2D-SL0) [21], 2D-FISTA [22], 2D-ADMM [12], have been proposed. Unlike the above iterative optimization algorithm, which requires an extensive manual adjustment of model parameters to achieve optimal results in experiments, the SBL-based method does not. However, many Bayesian algorithms struggle with directly handling two-dimensional sparse matrices and face challenges such as matrix inversion and high computational complexity. Zhang et al. [23] proposed a 2D-Fast SBL for ISAR joint imaging. However, this method requires converting the complex matrix into a real matrix, leading to increased computational complexity.

The relative motion of a target to the radar can generally be decomposed into two parts: translation and rotation [24]. Achieving high-resolution ISAR imaging requires precise compensation for the target’s translational motion. However, accurately estimating the motion parameters of the target can be challenging, particularly in SRSF ISAR scenarios with sparse subpulses and gapped apertures, which further affect algorithm performance. Huang et al. [3] introduced an RSF ISAR motion compensation and imaging method, approximating the target’s translational motion as uniform. However, this simplistic approximation leads to residual phase errors, making focusing challenging. An approach to mitigate this is integrating autofocusing into the ISAR image reconstruction process. Zhao et al. [25] used maximum likelihood estimation (MLE) criterion to estimate phase errors during ISAR imaging, reducing the phase error impact. Zhang et al. [26] used the minimum entropy criterion as the optimization goal for the autofocusing process, aiming to improve the quality of sparse aperture ISAR imaging. Li et al. [27] proposed a joint imaging and autofocuing method using MLE criterion to estimate the phase error and a 2D-ADMM method for ISAR imaging. The aforementioned algorithm can mitigate the influence of phase errors on the cross-range direction and produce a well-focused ISAR image. RSF radar synthesizes high-resolution profiles in the distance dimension through multi-pulse synthesis, so there is not only a phase error in the cross-range direction but also a phase error in the distance dimension, posing challenges for existing algorithms regarding effectively addressing range-dimensional defocus. Thus, Shao et al. [28] proposed a 2D joint ISAR imaging and motion compensation algorithm for sparse stepped frequency radar, addressing translation-induced phase errors. However, it cannot handle phase errors from model inaccuracies or system defects. To tackle this issue, Wang et al. [29] introduced imaging autofocusing and motion compensation algorithms. However, vectorizing ISAR images leads to an increased computational burden. Using the 2D-ADMM algorithm, Lv et al. [30] proposed 2D-IADIA to alternately obtain ISAR imaging and phase estimation. However, the quality of the imaging results is highly dependent on the selection of model regularization parameters. To avoid manual parameter tuning, Lv et al. [31] further unfolded 2D-IADIA into a deep neural network, which has better robustness than traditional algorithms.

In this paper, we propose a 2D joint imaging and autofocusing method named 2D-autofocusing and complex-value inverse-free SBL (2D-AFCIFSBL) for SRSF ISAR, aiming to achieve high-resolution imaging from 2D sparse received data. To begin, to leverage the full potential of 2D coupling information of echoes, we develop a vectorized representation model that integrates high-resolution ISAR imaging and 2D phase error estimation. After that, the probabilistic modeling is performed on the ISAR image, which has a gamma-complex Gaussian hierarchical imposed prior to ensure the sparsity of the model, and then variational Bayesian inference is used for image reconstruction in the complex domain. Due to the enormous computational complexity caused by matrix inversion at each iteration in Bayesian inference, matrix inversion is avoided by relaxing the evidence lower bound function (ELBO), which transforms the covariance matrix to be inverted into a diagonal matrix. So as to accomplish phase autofocusing during the ISAR image reconstruction process, the 2D phase error is estimated through MLE in the image reconstruction iteration. To further alleviate the storage pressure and improve the imaging efficiency, the iterative process is reformulated into a matrix form by leveraging the properties of the Kronecker product. Finally, both simulated and measured dataset experiments demonstrate that the proposed 2D joint method has excellent performance and efficiency compared to other joint imaging and autofocusing methods.

The remainder of this article is organized as follows. Section 2 formulates a 2D imaging model incorporating a 2D phase error for SRSF ISAR. Section 3 proposes the joint imaging and autofocusing probabilistic model and provides a detailed derivation and discussion of the 2D-AFCIFSBL method for solving this model. Section 4 evaluates the performance of the proposed method by simulated and measured datasets compared to the 2D-Fast SBL, 2D-IADIA and 2D-UADN methods. Section 5 concludes this article and discusses future work. Unless specifically defined, the following mathematical symbols will be used in this article. The bold letters denote matrices or vectors. Given a matrix $\mathit{A}$, ${\mathit{A}}_{i·}$, ${\mathit{A}}_{·j}$, ${\mathit{A}}_{ij}$, ${\mathit{A}}^{-1}$, ${\mathit{A}}^{*}$, ${\mathit{A}}^{\mathrm{T}}$ and ${\mathit{A}}^{\mathrm{H}}$ represent the ith row, jth column, $(i,j)$th entry, inverse, conjugate, transpose, and conjugate transpose of $\mathit{A}$, respectively. ⊗, ⊙ and ⊘ represent the Kronecker, Hadamard products and element-wise division, respectively. $\mathbb{R}$ and $\mathbb{C}$ denotes the real and complex domains, respectively.

2. Imaging Model for RSF ISAR

The RSF ISAR system transmits N subpulses in a burst to acquire high resolution in the range direction, and it transmits $N_a$ bursts in a CPI to acquire high resolution in the cross-range direction. Assuming $n=[0,1,\cdots ,N-1]$ and $n_a=[0,1,\cdots ,N_a-1]$, we can express the nth subpulse in the $n_a$th burst of the transmitted signal as follows

$$s_t(t,n,n_a)=\mathrm{rect}\left(\frac{t-n_{a}N{T}_r-n{T}_r}{T_p}\right)exp\left[j2\pi f_n(t-N_{a}N{T}_r-n{T}_r)\right]$$

(1)

where $T_r$, $T_p$, and $f_n$ denote the subpulses repetition interval, the subpulses bandwidth and the carrier frequency of the nth subpulse, respectively. Random step frequency signals show stronger anti-interference characteristics compared to traditional step frequency signals, which is primarily due to the random frequency hopping step size. Here, the carrier frequency of the nth subpulse is given by $f_n=f_c+\Gamma \left(n\right)\Delta f$, where $f_c$ is the initial carrier frequency, $\Gamma \left(n\right)$ is the nth randomized modulation code chosen from the set of $[0,1,2,\cdots ,N-1]$, and $\Delta f$ is the stepped frequency. Moreover, $\Delta f$ is defined as $\Delta f=1/T_p$, resulting in the synthetic bandwidth $B=N\Delta f$. $\mathrm{rect}(·)$ is the rectangular function denoted as

$$\mathrm{rect}\left(t\right)=\left\{\begin{array}{cc}\hfill 1\hfill & \hfill 0\le t\le 1\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$$

(2)

Figure 1 shows the ISAR imaging turntable model. We assume that the electromagnetic scattering characteristics of the target can be modeled as the superposition of K strong scattering centers. Then, the received echo signal can be denoted as

$$s_r(\widehat{t},n,n_a)=\sum _{k=1}^K{\sigma }_k\mathrm{rect}\left(\frac{\widehat{t}-2R_k\left(t\right)/c}{T_p}\right)exp\left[j2\pi f_n(\widehat{t}-2R_k\left(t\right)/c)\right]$$

(3)

where $\widehat{t}=t-n_{a}N{T}_r-n{T}_r$ is fast time and ${\sigma }_k$ denotes the scattering coefficient of the kth strong scattering center. $R_k\left(t\right)$ denotes the instantaneous distance between the kth strong scattering center and the radar at time t.

Based on the “stop-and-go” assumption, it holds that $R_k\left(t\right)\approx R_k\left(t_{n,n_a}\right)$, where $t_{n,n_a}=n_{a}N{T}_r+n{T}_r$. And $R_k\left(t_{n,n_a}\right)$ can be deduced as follows

$$R_k\left(t_{n,n_a}\right)=R_O\left(t_{n,n_a}\right)+x_{k}sin\left(\omega t_{n,n_a}\right)-y_{k}cos\left(\omega t_{n,n_a}\right)$$

(4)

where $R_O\left(t_{n,n_a}\right)\approx r_0+v{t}_{t_{n,n_a}}+a/2t_{t_{n,n_a}}^2$ is the instantaneous radar distance of the specific point O at time $t_{n,n_a}$, which represents the translational motion of the target. And $r_0$, v, a is the initial range between radar and the point O, the speed of the target and the acceleration of the target, respectively. $x_{k}sin\left(\omega t_{n,n_a}\right)+y_{k}cos\left(\omega t_{n,n_a}\right)$ denotes the rotational motion of the kth scatter, while $(x_k,y_k)$ is the coordinate of the kth scattering point on the imaging plane. $\omega$ represents the angular speed. Due to the assumption of a short CPI, the equivalent rotation angle during the CPI is not very large. Therefore, Equation (4) can be expanded using a Taylor series and neglecting terms of the second order and higher, leading to

$$R_k\left(t_{n,n_a}\right)\approx R_O\left(t_{n,n_a}\right)+x_k\omega t_{n,n_a}-y_k$$

(5)

Before imaging, compensating for translational motion is necessary. However, traditional motion compensation methods suffer from low accuracy due to the downsampling of the echo signal [32]. The random selection of subpulse carrier frequencies poses a challenge in accurately estimating the parameters of translational motion using parametric methods. Since the translational motion parameters of the target cannot be accurately estimated, there will still be residual phase error after translation compensation. Therefore, the sampled echo neglecting the constant term can be denoted as

$$s_r(n,n_a)\approx exp\left(j\varphi (n,n_a)\right)\sum _{k=1}^K{\sigma }_{k}exp[j4\pi \Gamma \left(n\right)\Delta f\frac{y_k}{c}]exp[-j4\pi f_c\frac{x_k\omega n_{a}N{T}_r}{c}]$$

(6)

where $\varphi (n,n_a)$ consists of two parts: the phase residual after translation motion compensation and the model error [29]. It can be expressed as

$$\varphi (n,n_a)=\frac{4\pi f_n((v-\widehat{v})t_{n,n_a}+(a-\widehat{a})t_{n,n_a}^2/2)}{c}+\Delta \varphi (n,n_a)$$

(7)

where $\widehat{v}$, $\widehat{a}$ denote the estimated speed of acceleration of the target, and $\Delta \varphi (n,n_a)$ represents the model error.

The SRSF echo signal matrix is obtained by random sampling from the RSF echo signal matrix. Assume the SRSF echo signal contains H$(H\le N_a)$ bursts and M$(M\le N)$ subpulses in each burst, and the echo signal data matrix $\mathbf{S}\in {\mathbb{C}}^{M\times H}$ can be written as

$$\mathbf{S}={\mathsf{\Phi }}_r{\mathbf{E}}^{\prime }{{\mathsf{\Phi }}_a}^T\odot \left[{\mathsf{\Phi }}_r{{\mathit{F}}^{\prime }}_r\mathbf{X}{{{\mathit{F}}^{\prime }}_a}^T{{\mathsf{\Phi }}_a}^T\right]+\mathit{G}=\mathbf{E}\odot \left[{\mathit{F}}_r\mathbf{X}{{\mathit{F}}_a}^T\right]+\mathit{G}$$

(8)

where ${\mathbf{E}}^{\prime }\in {\mathbb{C}}^{N\times N_a}$, $\mathbf{X}\in {\mathbb{C}}^{P\times Q}$, $\mathit{G}\in {\mathbb{C}}^{M\times H}$ represent the phase error matrix, the 2D ISAR image and the noise matrix, respectively. Among them, $\mathbf{S}$ and ${\mathbf{E}}^{\prime }$ can be denoted as

$$\mathbf{S}={\left[\begin{array}{cccc}{s}_r(1,1)& s_r(1,2)& \cdots & s_r(1,H)\\ s_r(2,1)& s_r(2,2)& \cdots & s_r(2,H)\\ ⋮& ⋮& \ddots & ⋮\\ s_r(M,1)& s_r(M,2)& \cdots & s_r(M,H)\end{array}\right]}_{M\times H}$$

(9)

$${\mathbf{E}}^{\prime }={\left[\begin{array}{cccc}e(1,1)& e(1,2)& \cdots & e(1,N_a)\\ e(2,1)& e(2,2)& \cdots & e(2,N_a)\\ ⋮& ⋮& \ddots & ⋮\\ e(N,1)& e(N,2)& \cdots & e(N,N_a)\end{array}\right]}_{N\times N_a}$$

(10)

where $e(2,2)=exp\left(j\varphi (n,n_a)\right)$. ${{\mathit{F}}^{\prime }}_r\in {\mathbb{C}}^{N\times P}(N\le P)$ and ${{\mathit{F}}^{\prime }}_a\in {\mathbb{C}}^{N_a\times Q}(N_a\le Q)$ denote the range and the cross-range sensing matrix, respectively, which can be written as

$${{\mathit{F}}^{\prime }}_r={\left[\begin{array}{cccc}{f}_r(1,1)& f_r(1,2)& \cdots & f_r(1,P)\\ f_r(2,1)& f_r(2,2)& \cdots & f_r(2,P)\\ ⋮& ⋮& \ddots & ⋮\\ f_r(N,1)& f_r(N,2)& \cdots & f_r(N,P)\end{array}\right]}_{N\times P}$$

(11)

$${{\mathit{F}}^{\prime }}_a={\left[\begin{array}{cccc}{f}_a(1,1)& f_a(1,2)& \cdots & f_a(1,Q)\\ f_a(2,1)& f_a(2,2)& \cdots & f_a(2,Q)\\ ⋮& ⋮& \ddots & ⋮\\ f_a(N_a,1)& f_a(N_a,2)& \cdots & f_a(N_a,Q)\end{array}\right]}_{N_a\times Q}$$

(12)

where $f_r(n,p)=exp(j2\pi \Gamma \left(n\right)\frac{p-1}{P})$ and $f_a(n_a,q)=exp(-j2\pi n_a\frac{q-1}{Q})$. ${\mathsf{\Phi }}_r\in {\mathbb{R}}^{M\times N}$ and ${\mathsf{\Phi }}_a\in {\mathbb{R}}^{H\times N_a}$ represent the randomly downsampling matrices, which can be acquired by randomly selecting M and H rows of two identity matrices with sizes of $N\times N$ and $N_a\times N_a$, respectively. ${\mathit{F}}_a={{\mathit{F}}^{\prime }}_a{\mathsf{\Phi }}_a$ and ${\mathit{F}}_r={\mathsf{\Phi }}_r{{\mathit{F}}^{\prime }}_r$.$\mathbf{E}={\mathsf{\Phi }}_r{\mathbf{E}}^{\prime }{{\mathsf{\Phi }}_a}^T$ represent the two-dimensional phase error matrix after downsampling.

Due to the impact of the two-dimensional phase error $\mathbf{E}$, several 2D reconstruction algorithms like 2D-SL0 and 2D-FISTA, as well as algorithms that only consider one-dimensional phase error such as 2D-UADN [27], struggle to effectively reconstruct ISAR images. Although the 2D-IADIA [30] algorithm, which accounts for two-dimensional phase errors, can reconstruct ISAR images, it is difficult to adjust the regularization coefficient. While the SBL algorithm does not require manual parameter adjustment, most current SBL algorithms are not directly applicable to matrix solving and have high computational complexity. A straightforward approach involves using the Kronecker product [33] to transform the 2D ISAR model from Equation (8) into a one-dimensional model, representing it in an equivalent vector form, and subsequently applying Bayesian algorithms for reconstruction. Then, to improve computational efficiency and alleviate storage pressure, the vectorized iteration is transformed into matrix iteration. Taking into account the effect of the two-dimensional phase error $\mathbf{E}$, the vectorized equivalent model of Equation (8) can be formulated as follows

$$\mathbf{s}=\mathbf{e}({\mathit{F}}_a\otimes {\mathit{F}}_r)\mathbf{x}+\mathit{g}=\mathbf{e}\mathit{f}\mathbf{x}+\mathit{g}$$

(13)

where $\mathbf{s}=\mathrm{vec}\left(\mathbf{S}\right)\in {\mathbb{C}}^{MH\times 1}$, $\mathbf{x}=\mathrm{vec}\left(\mathbf{X}\right)\in {\mathbb{C}}^{PQ\times 1}$, $\mathit{g}=\mathrm{vec}\left(\mathit{G}\right)\in {\mathbb{C}}^{MH\times 1}$. The dictionary matrix $\mathit{f}={\mathit{F}}_a\otimes {\mathit{F}}_r\in {\mathbb{C}}^{MH\times PQ}$, $\mathrm{vec}(·)$ denotes the vectorization operation. $\mathbf{e}\in {\mathbb{C}}^{MH\times MH}$ is a diagonal matrix, which can be denoted as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{e}& \hfill =\mathrm{diag}\left\{e^{j{\varphi }_{1,1}},\cdots ,e^{j{\varphi }_{M,1}},e^{j{\varphi }_{1,2}},\cdots ,e^{j{\varphi }_{M,H}}\right\}\hfill \\ \hfill & \hfill ={\left[\begin{array}{cccc}{e}^{j{\varphi }_{1,1}}& 0& \cdots & 0\\ 0& e^{j{\varphi }_{2,1}}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & e^{j{\varphi }_{M,H}}\end{array}\right]}_{MH\times MH}\hfill \end{array}\end{array}$$

(14)

where $\mathrm{diag}(·)$ is an operation creating a diagonal matrix using the elements inside.

3. Proposed 2D-AFCIFSBL Method

In this section, we first derive the 1D complex inverse-free SBL algorithm based on the model (13). To improve computational efficiency, the properties of the Kronecker product are then used to transform the iterative process into a matrix form. Finally, two-dimensional self-focusing is integrated into the imaging process to develop the 2D-AFCIFSBL algorithm.

3.1. ISAR Imaging Based on Variational Bayesian Inference

In this section, a probabilistic model is established for the ISAR imaging shown in (13), and then variational Bayesian inference is used for image reconstruction.

In the SBL framework, we usually assume that the noise follows zero-mean complex Gaussian distribution $\mathcal{CN}\left(g\right|0,{\delta }^{-1}\mathit{I})$, where $\mathcal{CN}(·)$ denotes complex Gaussian distribution and $\delta$ is the reciprocal of the noise variance. Therefore, the likelihood function of the received echo signal $\mathbf{s}$ can be written as

$$p\left(\mathbf{s}\right)=\mathcal{CN}\left(\mathbf{s}\right|\mathbf{e}\mathit{f}\mathbf{x},{\delta }^{-1}{\mathit{I}}_{MH})$$

(15)

where ${\mathit{I}}_{MH}$ is the identity matrix with size $MH\times MH$. $\delta$ has a gamma distribution applied

$$p\left(\delta \right)=Gamma\left(\delta \right|a,b)$$

(16)

where $Gamma(·)$ denotes gamma distribution, and a, b are the position and scale parameter, respectively. To promote the sparsity of the ISAR image, we impose a hierarchical Bayesian prior on $\mathbf{x}$. In the first layer, a distinctively Gaussian distribution is placed on each pixel of the sparse vector $\mathbf{x}$

$$p\left(\mathbf{x}\right)=\mathcal{CN}\left(\mathbf{x}\right|0,\mathsf{\alpha })$$

(17)

where $\mathsf{\alpha }=\mathrm{diag}({\mathsf{\alpha }}_1^{-1},{\mathsf{\alpha }}_2^{-1},\dots ,{\mathsf{\alpha }}_{PQ}^{-1})$. ${\mathsf{\alpha }}_i^{-1}$ denotes the variance of ${\mathbf{x}}_i$. ${\mathsf{\alpha }}_i$ and ${\mathbf{x}}_i$ are the ith element of $\mathsf{\alpha }$ and $\mathbf{x}$. In the second layer, we impose Gamma distribution with position parameter c and scale parameter d on the reciprocal of the variance ${\mathsf{\alpha }}_i$

$$p\left(\mathsf{\alpha }\right)=\prod _i^{PQ}Gamma\left({\mathsf{\alpha }}_i\right|c,d)$$

(18)

Finally, the corresponding probabilistic graphical model is illustrated in Figure 2. Our aim is to estimate the posterior of all latent variables $\mathsf{\theta }=\{\mathbf{x},\mathsf{\alpha },\delta \}$, whose posterior can be computed through the variational Bayesian inference algorithm [34]. In the traditional VB method, the evidence lower bound (ELBO) is expressed by

$$\mathrm{ELBO}\left(q\right)={\int }_{\mathsf{\theta }}q\left(\mathsf{\theta }\right)log\frac{p\left(\mathbf{s},\mathbf{x},\mathsf{\alpha },\delta \right)}{q\left(\mathsf{\theta }\right)}$$

(19)

where $p\left(\mathbf{s},\mathbf{x},\mathsf{\alpha },\delta \right)=p\left(\mathbf{s}\right|\mathbf{x},\delta )p\left(\mathbf{x}\right|\mathsf{\alpha })p\left(\mathsf{\alpha }\right)p\left(\delta \right)$, $q\left(\mathsf{\theta }\right)=q_{\mathbf{x}}\left(\mathbf{x}\right)q_{\mathsf{\alpha }}\left(\mathsf{\alpha }\right)q_{\sigma }\left(\sigma \right)$.

To speed up computation, we employed relaxed-ELBO as a replacement for the ELBO, which eliminates the need to invert a matrix [35]. Assume f is a Lipschitz continuous function, so there exists a number $T\ge \frac{\nabla f\left(\mathbf{x}\right)}{2},\forall \mathbf{x}$ such that for any $\mathit{u},\mathit{v}$, we have:

$$\begin{array}{c}\hfill f\left(\mathit{u}\right)⩽f\left(\mathit{v}\right)+{(\mathit{u}-\mathit{v})}^H\nabla f\left(\mathit{v}\right)+{(\mathit{u}-\mathit{v})}^{H}T(\mathit{u}-\mathit{v})\end{array}$$

(20)

Then, the lower bound on $p\left(\mathbf{s}\right|\mathbf{x},\delta )$ can be obtained as

$$p\left(\mathbf{s}\right|\mathbf{x},\delta )=\frac{{\delta }^{PQ}}{{\pi }^{PQ}}exp\left\{-\frac{\delta }{2}{\parallel \mathbf{s}-\mathbf{e}\mathit{f}\mathbf{x}\parallel }_2^2\right\}\ge G(\mathbf{s},\mathbf{x},\delta ;\mathbf{z})$$

(21)

where $G(\mathbf{s},\mathbf{x},\delta ;\mathbf{z})$ is defined as

$$\begin{array}{c}\hfill G(\mathbf{s},\mathbf{x},\delta ;\mathbf{z})=\frac{{\delta }^{PQ}}{{\pi }^{PQ}}exp\left(\begin{array}{c}-2\delta \Re \left\{{(\mathbf{x}-\mathbf{z})}^{\mathrm{H}}\left({\mathit{f}}^{\mathrm{H}}\mathit{f}\mathbf{z}-{\mathit{f}}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}\right)\right\}\\ {-\delta \parallel \mathbf{s}-\mathbf{e}\mathit{f}\mathbf{z}\parallel }_2^2-\frac{T}{2}\delta {\parallel \mathbf{x}-\mathbf{z}\parallel }_2^2\end{array}\right)\end{array}$$

(22)

where $\Re (·)$ denotes taking the real parts and T is slightly larger than the eigenvalue of matrix $2{\mathit{f}}^{\mathrm{H}}\mathit{f}$. When $\mathbf{x}=\mathbf{z}$, the inequality becomes an equality. Then, the relaxed-ELBO can be written as

$$\mathrm{relaxed}-\mathrm{ELBO}\left(q,\mathbf{z}\right)={\int }_{\mathsf{\theta }}q\left(\mathsf{\theta }\right)log\frac{F\left(\mathbf{s},\mathsf{\theta },\mathbf{z}\right)}{q\left(\mathsf{\theta }\right)}$$

(23)

where $F\left(\mathbf{s},\mathsf{\theta },\mathbf{z}\right)=G(\mathbf{s},\mathbf{x},\delta ;\mathbf{z})p\left(\mathbf{x}\right|\mathsf{\alpha })p\left(\mathsf{\alpha }\right)p\left(\delta \right)$ and $\mathsf{\theta }=\{\mathbf{x},\mathsf{\alpha },\delta \}$ is the latent variable to be estimated. In the VBE-step, $\mathbf{z}$ is held fixed, and Equation (24) is used iteratively to approximately compute the posterior probability density function for each latent variable.

$$ln{q}_{{\mathsf{\theta }}_i}\left({\mathsf{\theta }}_i\right)={\langle lnF(\mathbf{s},\mathsf{\theta },\mathbf{z})\rangle }_{q_{j\ne i}}+\mathrm{const}$$

(24)

Then, the posterior of $\mathbf{x}$ can be regarded a complex Gaussian distribution

$$\begin{array}{cc}\hfill q_{\mathbf{x}}\left(\mathbf{x}\right)& \hfill =\mathcal{CN}\left(\mathbf{x}\right|\mathsf{\mu },\mathsf{\Sigma })\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \mathsf{\mu }& \hfill =\langle \delta \rangle \mathsf{\Sigma }\left(\frac{T\mathbf{z}}{2}+{\mathit{f}}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}-{\mathit{f}}^{\mathrm{H}}\mathit{f}\mathbf{z}\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \mathsf{\Sigma }& \hfill ={\left(\frac{T\langle \delta \rangle }{2}\mathit{I}+\mathit{D}\right)}^{-1}\hfill \end{array}$$

(27)

where $\mathit{D}=\mathrm{diag}(\langle {\mathsf{\alpha }}_1\rangle ,\langle {\mathsf{\alpha }}_2\rangle ,\dots ,\langle {\mathsf{\alpha }}_{PQ}\rangle )$ and $\langle ·\rangle$ denotes the exception operation. It can be seen that the covariance matrix $\mathsf{\Sigma }$ is a diagonal matrix. The calculation of Equation (27) does not require the inversion operation of the matrix. Therefore, the computational efficiency is significantly improved compared to the traditional Bayesian algorithm.

For noise precision $\delta$, the approximate posterior $q_{\delta }\left(\delta \right)$ is the Gamma distribution, which can be written as

$$\begin{array}{cc}\hfill q_{\delta }\left(\delta \right)& \hfill =Gamma\left(\delta \right|\tilde{a},\tilde{b})\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \tilde{a}& \hfill =a+MH\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \tilde{b}& \hfill =b+\langle g(\mathbf{x},\mathbf{z})\rangle \hfill \end{array}$$

(30)

with

$$\begin{array}{cc}\hfill \langle g(\mathbf{x},\mathbf{z})\rangle & ={\parallel \mathbf{s}-\mathbf{e}\mathit{f}\mathbf{z}\parallel }_2^2+2\Re \left\{{(\mathbf{x}-\mathbf{z})}^{\mathrm{H}}\left({\mathit{f}}^{\mathrm{H}}\mathit{f}\mathbf{z}-{\mathit{f}}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}\right)\right\}\hfill \\ \hfill & +T(\parallel \mathsf{\mu }-\mathbf{z}{\parallel }_2^2+tr\left(\mathsf{\Sigma }\right))/2\hfill \end{array}$$

(31)

Thus, $\langle \delta \rangle$ can be expressed as the expectation of $\delta$

$$\langle \delta \rangle =\frac{\tilde{a}}{\tilde{b}}$$

(32)

The approximate posterior of $\mathsf{\alpha }$ is the distinctively Gamma distribution, which can be written as

$$\begin{array}{cc}\hfill q_{\mathsf{\alpha }}\left(\mathsf{\alpha }\right)& \hfill =\prod _{i=1}^{PQ}Gamma\left({\mathsf{\alpha }}_i\right|\tilde{c},\widetilde{d_i})\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill \tilde{c}& \hfill =c+1\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\tilde{\mathit{d}}}_i& \hfill =d+\langle {\mathbf{x}}_i^{*}{\mathbf{x}}_i\rangle \hfill \end{array}$$

(35)

where $\langle {\mathbf{x}}_i^{*}{\mathbf{x}}_i\rangle ={\mathbf{x}}_i^{*}{\mathbf{x}}_i+\left|{\mathsf{\Sigma }}_{i,i}\right|$. ${\mathsf{\Sigma }}_{i,i}$ denotes the ith element on the diagonal of matrix $\mathsf{\Sigma }$. Thus, $\langle \mathsf{\alpha }\rangle$ can be expressed as the expectation of $\mathsf{\alpha }$

$$\langle {\mathsf{\alpha }}_i\rangle =\frac{\tilde{c}}{\widetilde{{\mathit{d}}_i}}$$

(36)

In the VBM step, the approximate posterior $q\left(\mathsf{\theta }\right)$ remains fixed.

$${\mathbf{z}}^{new}=arg\underset{\mathbf{z}}{max}{\langle lnF(\mathbf{s},\mathsf{\theta },\mathbf{z})\rangle }_{q(\mathsf{\theta };{\mathbf{z}}^{old})}=Q\left(\mathbf{z}\right|{\mathbf{z}}^{old})$$

(37)

By differentiating Equation (37) and setting it equal to zero, we can find the value of ${\mathbf{z}}^{new}$.

$$\frac{\partial Q\left(\mathbf{z}\right|{\mathbf{z}}^{old})}{\partial \mathbf{z}}=(T{\mathit{I}}_{PQ}-2{\mathit{f}}^{\mathrm{H}}\mathit{f})(\mathbf{z}-\mathsf{\mu })=0$$

(38)

Since T is slightly larger than the maximum eigenvalues of matrix $2{\mathit{f}}^{\mathrm{H}}\mathit{f}$, $T{\mathit{I}}_{PQ}-2{\mathit{f}}^{\mathrm{H}}\mathit{f}$ is a positive-definite matrix. Therefore, we can obtain the value of ${\mathbf{z}}^{new}$ as

$${\mathbf{z}}^{new}=\mathsf{\mu }$$

(39)

By alternately and iteratively executing the VBE step and VBM step, the vectorized ISAR image $\mathbf{x}=\mathsf{\mu }$ is finally obtained. Although the above process of reconstructing the sparse vector $\mathbf{x}$ avoids calculating the inverse of the matrix $\mathsf{\Sigma }$, converting the 2D echo and ISAR images into vectors imposes significant computational and storage burdens. As the dimensions N, $N_a$, P, and Q increase, the observation matrix $\mathit{f}$ also rapidly expands, leading to computational challenges. In contrast, matrix operations offer higher computational efficiency and require less storage space than vectorized calculations. Therefore, it is necessary to convert the above calculation process into matrix operations to efficiently address the matrix sparse problem. Assuming $\mathsf{\mu }=\mathrm{vec}\left(\mathit{U}\right),\mathbf{z}=\mathrm{vec}\left(\mathbf{Z}\right)$, we can rewrite the third part of Equation (26) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathit{f}}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}-{\mathit{f}}^{\mathrm{H}}\mathit{f}\mathbf{z}& \hfill ={({\mathit{F}}_a\otimes {\mathit{F}}_r)}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}-{({\mathit{F}}_a\otimes {\mathit{F}}_r)}^{\mathrm{H}}({\mathit{F}}_a\otimes {\mathit{F}}_r)\mathbf{z}\hfill \\ \hfill & \hfill =({\mathit{F}}_a^{\mathrm{H}}\otimes {\mathit{F}}_r^{\mathrm{H}}){\mathbf{e}}^{\mathrm{H}}\mathbf{s}-\left({\mathit{F}}_a^{\mathrm{H}}{\mathit{F}}_a\right)\otimes \left({\mathit{F}}_r^{\mathrm{H}}{\mathit{F}}_r\right)\mathbf{z}\hfill \\ \hfill & \hfill =\mathrm{vec}\left({\mathit{F}}_r^{\mathrm{H}}\left({\mathbf{E}}^{*}\odot \mathbf{S}-{\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}\right){\mathit{F}}_a^{*}\right)\hfill \end{array}\end{array}$$

(40)

Then, Equations (26) and (27) can be expressed as the matrix form

$$\begin{array}{cc}\hfill & \hfill \mathit{U}=\langle \delta \rangle \left(\frac{T\mathbf{Z}}{2}+{\mathit{F}}_r^{\mathrm{H}}\left({\mathbf{E}}^{*}\odot \mathbf{S}-{\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}\right){\mathit{F}}_a^{*}\right)\odot {\mathsf{\Sigma }}_{P\times Q}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & \hfill {\mathsf{\Sigma }}_{P\times Q}={\mathbf{1}}_{P\times Q}\oslash (\frac{T\langle \delta \rangle }{2}{\mathbf{1}}_{P\times Q}+{\mathit{D}}_{P\times Q})\hfill \end{array}$$

(42)

where ${\mathsf{\Sigma }}_{P\times Q}=\aleph \left(\mathsf{\Sigma }\right)$, ${\mathit{D}}_{P\times Q}=\aleph \left(\mathit{D}\right)$. ℵ represent the extraction of the diagonal elements of matrix $\mathsf{\Sigma }$ to form a matrix of size $P\times Q$. ${\mathbf{1}}_{P\times Q}$ is a matrix with all elements equal to 1, and its size is $P\times Q$. $T=2{\lambda }_{max}\left({\mathit{F}}_r^{\mathrm{H}}{\mathit{F}}_r\right)·{\lambda }_{max}\left({\mathit{F}}_a^{\mathrm{H}}{\mathit{F}}_a\right)+{10}^{-5}$ is slightly larger than the Lipschitz constant. ${\lambda }_{max}(·)$ denotes the extraction of the maximum eigenvalue of a matrix.

We can rewrite Equations (32) and (36) in matrix form

$$\langle \delta \rangle =\frac{a+MH}{b+\langle g(\mathbf{X},\mathbf{Z})\rangle }$$

(43)

$${\mathit{D}}_{P\times Q}=(c+1)\oslash (d+\mathit{U}\odot {\mathit{U}}^{*}+|{\mathsf{\Sigma }}_{P\times Q}\left|\right)$$

(44)

where $\langle g(\mathbf{X},\mathbf{Z})\rangle$ can be rewritten in matrix form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \langle g(\mathbf{X},\mathbf{Z})\rangle & ={∥\mathbf{S}-\mathbf{E}\odot \left({\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}\right)∥}_F^2+2·\mathrm{sum}\left(\Re \left\{{(\mathit{U}-\mathbf{Z})}^{*}\odot \left[{\mathit{F}}_r^{\mathrm{H}}({\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}-{\mathbf{E}}^{*}\odot \mathbf{S}){\mathit{F}}_a^{*}\right]\right\}\right)\hfill \\ \hfill & +\frac{{\parallel \mathit{U}-\mathbf{Z}\parallel }_F^2+\mathrm{sum}\left({\mathsf{\Sigma }}_{P\times Q}\right)}{2}T\hfill \end{array}\end{array}$$

(45)

The second part of (45) can be obtained from the second part of (31)

$$\begin{array}{c}\begin{array}{cc}\hfill & 2\Re \left\{{(\mathbf{x}-\mathbf{z})}^{\mathrm{H}}\left({\mathit{f}}^{\mathrm{H}}\mathit{f}\mathbf{z}-{\mathit{f}}^{\mathrm{H}}{\mathbf{e}}^{\mathrm{H}}\mathbf{s}\right)\right\}\hfill \\ \hfill & =2\Re \left\{\mathrm{vec}{\left(\mathit{U}-\mathbf{Z}\right)}^{\mathrm{H}}\left[{({\mathit{F}}_a\otimes {\mathit{F}}_r)}^{\mathrm{H}}\left(({\mathit{F}}_a\otimes {\mathit{F}}_r)\mathbf{z}-{\mathbf{e}}^{\mathrm{H}}\mathbf{s}\right)\right]\right\}\hfill \\ \hfill & =2\Re \left\{\mathrm{vec}{\left(\mathit{U}-\mathbf{Z}\right)}^{\mathrm{H}}\mathrm{vec}\left[{\mathit{F}}_r^{\mathrm{H}}({\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}-{\mathbf{E}}^{*}\odot \mathbf{S}){\mathit{F}}_a^{*}\right]\right\}\hfill \\ \hfill & =2·\mathrm{sum}\left(\Re \left\{{(\mathit{U}-\mathbf{Z})}^{*}\odot \left[{\mathit{F}}_r^{\mathrm{H}}({\mathit{F}}_r\mathbf{Z}{\mathit{F}}_a^{\mathrm{T}}-{\mathbf{E}}^{*}\odot \mathbf{S}){\mathit{F}}_a^{*}\right]\right\}\right)\hfill \end{array}\hfill \end{array}$$

(46)

where ${\parallel ·\parallel }_F$ represents the Frobenius norm and $\mathrm{sum}(·)$ denotes the summation of matrix elements. And Equation (39) can rewritten as

$$Z^{new}=U$$

(47)

In the above derivation, while we assume the phase error matrix $\mathbf{E}$ is known, it is, in fact, unknown. Therefore, it becomes necessary to estimate the phase error during ISAR imaging.

3.2. Autofocusing Based on MLE

In this subsection, we estimate the two-dimension phase error matrix $\mathbf{E}$ through MLE in the image reconstruction iteration. Assuming that the ISAR image is known during the autofocusing process, the maximum likelihood estimation (MLE) problem can be converted into the following unconstrained optimization problem

$$\widehat{\mathbf{e}}=arg\underset{\mathbf{e}}{min}{∥\mathbf{s}-\mathbf{e}\mathit{f}\mathbf{x}∥}_2^2$$

(48)

For the lth iteration, the aforementioned Equation (48) can be transformed into an optimized form as follows

$$arg\underset{\mathbf{E}}{min}\sum _{m=1}^M\sum _{h=1}^H{∥{\mathbf{S}}_{m,h}-e^{j{\varphi }_{m,h}}{\left({\mathit{F}}_r{\left({\mathbf{X}}^{\left(l\right)}{\mathit{F}}_a^{\mathrm{T}}\right)}_{·,h}\right)}_{m,·}∥}_2^2$$

(49)

where ${\mathbf{X}}^{\left(l\right)}$ denotes the imaging result of the lth iteration. We initially decompose Equation (49) into the following minimization problem form

$${\widehat{e}}_{l+1}^{m,h}=arg\underset{\mathbf{E}}{min}{∥{\mathbf{S}}_{m,h}-e^{j{\varphi }_{m,h}}{\left({\mathit{F}}_r{\left({\mathbf{X}}^{\left(l\right)}{\mathit{F}}_a^{\mathrm{T}}\right)}_{·,h}\right)}_{m,·}∥}_2^2$$

(50)

where ${\widehat{e}}_{l+1}^{m,h}$ represents the estimated phase error of the mth subpulse in the hth burst. Therefore, the optimal solution can be obtained by setting the first derivative of Equation (50) to zero

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\partial }{\partial {\widehat{e}}_{l+1}^{m,h}}{∥{\mathbf{S}}_{m,h}-e^{j{\varphi }_{m,h}}{\left({\mathit{F}}_r{\left({\mathbf{X}}^{\left(l\right)}{\mathit{F}}_a^{\mathrm{T}}\right)}_{·,h}\right)}_{m,·}∥}_2^2=0\hfill \\ \hfill & ⟹{\widehat{e}}_{l+1}^{m,h}=exp\left\{j·\left[\mathrm{angle}\left[{\mathbf{S}}_{m,h}·{\left({\mathit{F}}_r{\mathbf{X}}^{\left(l\right)}{\mathit{F}}_a^{\mathrm{T}}\right)}_{m,h}^{*}\right]\right]\right\}\hfill \end{array}\end{array}$$

(51)

where $\mathrm{angle}(·)$ is the operator of acquiring the phase information. To reduce memory requirements and improve computational efficiency, the estimation of the 2D phase error matrix can be further denoted in the matrix form

$${\mathbf{E}}^{\left(l\right)}=exp\left\{j·\left[\mathrm{angle}\left[\mathbf{S}\odot {\left({\mathit{F}}_r{\mathbf{X}}^{\left(l\right)}{\mathit{F}}_a^{\mathrm{T}}\right)}^{*}\right]\right]\right\}$$

(52)

In summary, detailed information on the proposed 2D fast ISAR imaging and autofocusing method has been presented in Algorithm 1.

Algorithm 1 2D-AFCIFSBL method

Input: $\mathbf{S}\in {\mathbb{C}}^{M\times H}$ Output: $\mathbf{E},\mathbf{X}$   1: Initialization: $\langle {\delta }^{\left(0\right)}\rangle \leftarrow 1,{\mathit{D}}_{P\times Q}^{\left(0\right)}\leftarrow {\mathbf{1}}_{P\times Q},{\mathbf{E}}^{\left(0\right)}\leftarrow {\mathbf{1}}_{P\times Q},{\mathbf{X}}^{\left(0\right)}\leftarrow {\mathit{F}}_r^{\mathrm{H}}\mathbf{S}{\mathit{F}}_a^{*}$   2: for $l=1,2,3,\cdots ,L_{max}$ do   3:    Update ${\mathbf{X}}^{\left(l\right)}$ by ${\mathbf{X}}^{\left(l\right)}\leftarrow \mathit{U}$ using (41), (42)   4:    Update ${\mathit{D}}_{P\times Q}^{\left(l\right)}$ using (44)   5:    Update $\langle {\delta }^{\left(l\right)}\rangle$ using (43)   6:    Update ${\mathbf{Z}}^{\left(l\right)}$ using (47)   7:    Update ${\mathbf{E}}^{\left(l\right)}$ using (52)   8:    if ${∥{\mathbf{X}}^{\left(l\right)}-{\mathbf{X}}^{(l-1)}∥}_F^2/{∥{\mathbf{X}}^{(l-1)}∥}_F^2\le {10}^{-4}$ then   9:       break 10:    end if 11: end for $\mathbf{X}={\mathbf{X}}^{\left(l\right)},\mathbf{E}={\mathbf{E}}^{\left(l\right)}$

Here, $L_{max}$ is the maximum number of iterations of the algorithm. Finally, we provide a brief analysis of the complexity of the proposed 2D-AFCIFSBL algorithm. The complexity of each iteration of the 2D-AFCIFSBL algorithm mainly comes from the update of latent variables. For updating $\mathbf{X}$ and $\langle \delta \rangle$ by (41), (42), (32), and (45), the computational complexity, which are $O\left(PQ\right(M+H+1)+MH(P+Q+1\left)\right)$ and $O\left(PQ\right(M+2H+2)+MH(2P+Q+3\left)\right)$, respectively, mainly come from four matrix multiplications, element-wise multiplication and element-wise division. The update process for the 2D phase error matrix $\mathbf{E}$ by Equation (52), which has a computational complexity of $O(PQH+MH(P+1\left)\right)$, primarily involves three matrix multiplications. Then, assuming $P=Q$ and $M,H<P$, the worst case complexity of 2D-AFCIFSBL is $O(L_{max}\left(P^2(2M+4H+2)\right)$.

4. Experiments and Analysis

In this section, we verified the effectiveness of the proposed RSF-ISAR imaging algorithm through simulated and measured datasets. Sparse data are obtained by randomly sampling in the range direction and the cross-range direction, respectively. The 2D SPR [30] is defined as $\gamma =(M/N,H/Na)$. For comparison, the 2D-Fast SBL method [23], 2D-UADN method [27] and 2D-IADIA method [30] are also adopted. For $P=Q$ and $M,N<P$, the worst case complexity of 2D-Fast SBL, 2D-UADN and 2D-IADIA is $O\left(L_{max}\left(4P^2(4M+4H+1)\right)\right)$, $O(L_{max}\left(P^2(2M+2H+1)\right)$ and $O\left(L_{max}\left(P^2(M+2H)\right)\right)$, respectively, where $L_{max}$ is the maximum number of iterations. In addition, we utilized the image entropy (IE) [36] and target-to-background ratio (TBR) [37] to evaluate the imaging performance of different algorithms. For a well-focused image, the lower IE and the higher TBR values are expected. All experiments were conducted on a computer with an Inter Core i9-10900K CPU and 32GB RAM.

4.1. Experimental Results on Simulated Dataset

We preliminarily validated the imaging performance of the proposed 2D-AFCIFSBL algorithm by using the simulated dataset. The simulated dataset is generated by an X-band RSF radar with 10 GHz initial carrier frequency and 640 MHz bandwidth, and its main parameters are presented in

Table 1. The main parameters of the X-band RSF radar.

Parameters	Value
Initial carrier frequency $f_c$	10 GHz
Synthetic bandwidth B	640 MHz
Number of sub-pules N	256
Number of bursts $N_a$	256
Stepped frequency $\Delta f$	2.5 MHz
Angular speed $\omega$	0.01 rad/s

. The target model in the simulated dataset consists of 74 independent scattering points presented in Figure 3a and its ISAR imaging result by using the 2D-AFCIFSBL method with complete echo data and no phase error added, as illustrated in Figure 3b. It is evident from the ISAR imaging result that the range between the target and the radar is 5060 m.

In the following experiments, we will add three types of 2D phase errors, namely random, linear and mixed phase errors to the raw simulation data separately as illustrated in Figure 4, to compare the phase error estimation capabilities of four methods. Among the three types of 2D phase errors, the amplitude variance of the random phase error is $\pi /2$, the linear phase error varies linearly between $-\pi /2$ and $\pi /2$ and the mixed phase error is the mean of the random phase error and linear phase error. Additionally, we set the SNR of the experiment to 20 dB by adding complex-valued Gaussian noise into the raw simulated data.

The first experiment is conducted under the $(0.8,0.8)$ SPR and the imaging results acquired by aforementioned four methods, as illustrated in Figure 5. From the first column, it is evident that the imaging results produced by 2D-Fast SBL are contaminated by a large number of pseudo-scatterers under linear phase error. This limitation stems from the 2D-Fast SBL method’s incapacity to estimate and compensate for linear phase error. However, the imaging effect under random phase error is better than that under linear phase error because the random phase error can be regarded as random noise, and the SBL-based method can utilize prior sparsity to denoise. For the 2D-UADN method, it can estimate the phase error in the cross-range direction but does not account for the phase error in the range direction. Consequently, the resulting images acquired by the 2D-UADN method still have some pseudo-scatterers and noise floor under the linear phase error. In contrast, the proposed 2D-AFCIFSBL method consistently produces neat and well-focused images under all given types of phase errors, which verifies its robustness against phase errors. Particularly under the linear phase error condition, the imaging results acquired by using the 2D-Fast SBL method and 2D-UADN method exhibit significant noise floor, while the proposed 2D-AFCIFSBL method still achieves excellent results, implying its superior performance relative to other methods. In addition, through comparison of the third and fourth columns of images, it is evident that the 2D-IADIA and 2D-AFCIFSBL methods yield similar imaging results. But the 2D-AFCIFSBL method has lower image entropy than the 2D-IADIA method, as demonstrated in Figure 6. It is evident that due to the lack of 2D phase errors estimation, the image entropy of the 2D-Fast SBL method and 2D-UADN method is extremely high, especially in the condition of linear phase errors.

Figure 7 and Figure 8 show the imaging results under $(0.6,0.6)$ and $(0.4,0.4)$ SPR, and the corresponding image entropy of each SPR is presented in Figure 9 and Figure 10, respectively. It is evident from the image entropy curves that the four methods exhibit rapid convergence to a stable image entropy value. In particular, the 2D-UADN and 2D-Fast SBL methods acquire the highest image entropy under the random phase error and linear phase error, respectively, while the 2D-Fast SBL method has similar levels of image entropy to the 2D-UADN method under the mixed phase error. The imaging results in the first and second rows of Figure 7 and Figure 8 also demonstrate that under the linear phase error, the 2D-UADN method has slightly better imaging results than the 2D-Fast SBL method due to its 1D phase error estimation capability, while the 2D-Fast SBL method lacks any phase estimation. But the imaging results of the 2D-Fast SBL method and 2D-UADN method are both poor under all SPRs compared to the 2D-IADIA method and 2D-AFCIFSBL method. Upon comparing the imaging results between the 2D-IADIA and 2D-AFCIFSBL methods, it is evident that the image quality achieved by the 2D-IADIA method is comparable to that of the 2D-AFCIFSBL method. Only when the SPR is down to $(0.4,0.4)$ do both methods begin to have few pseudo-scatterers and noise floor, but the phase error estimation ability of the 2D-IADIA method decreases more obviously. Additionally, the convergence speed of the 2D-IADIA method notably decreases at the SPR of $(0.4,0.4)$, while the 2D-AFCIFSBL method remains relatively stable, as shown in Figure 10. Therefore, the aforementioned experiments further substantiate the excellent imaging performance of the proposed 2D-AFCIFSBL method under three types of phase error with different SPR.

Subsequently, we assessed the imaging performance of the aforementioned four methods under different SNR conditions and only compared the imaging results under the random phase error, where the variance of the amplitude is $\pi /2$. Additionally, the SPR was set to $(0.5,0.5)$, which means that 128 subpulses and 128 bursts were randomly sampled from the complete echo data. The imaging results obtained by the four methods are demonstrated in Figure 11, in which the three rows from top to bottom represent the imaging results at 10 dB, 5 dB, and 0 dB, respectively. It is evident from the results that the proposed 2D-AFCIFSBL method obtains the lowest image entropy under all the given SNRs shown in Figure 12. Although the 2D-UADN method possesses the capability to estimate the 1D phase error in the cross-range direction, its imaging performance is adversely impacted by the 2D random phase error [27]. This results in an increase in pseudo-scatterers and noise floor in the imaging results, leading to the highest image entropy and worst image quality. Compared to the 2D-UADN method, the 2D-Fast SBL method achieve better results by leveraging the prior sparsity to resist random noise as mentioned in the previous experimental analysis. Particularly at the SNR of 5 dB, both the 2D-Fast SBL and 2D-UADN methods exhibit significant degradation in imaging performance, while the proposed 2D-AFCIFSBL method consistently produces neat and well-focused images, denoting its resistance to the noise. We can also observe that when the SNR is down to 0 dB, the noise floor begins to appear in the imaging result of the 2D-IADIA method, while only few pseudo-scatterers appear in the imaging result of the 2D-AFCIFSBL method. In addition, as the SNR gradually decreases, the convergence speed of the 2D-IADIA method also gradually decreases compared with the 2D-AFCIFSBL method shown in Figure 12.

Figure 13 illustrates the quantitative performance comparisons of aforementioned four methods under two evaluation metrics of image entropy and computation time over 100 Monte Carlo trials with the SNR ranging from 0 to 20 dB. The curves in Figure 13a reveal that compared to the 2D-Fast SBL and 2D-UADN methods, the proposed 2D-AFCIFSBL algorithm achieves imaging results with the lowest image entropy, further demonstrating its excellent imaging performance. Additionally, under high SNR conditions, the proposed 2D-AFCIFSBL method exhibits comparable image entropy to the 2D-IADIA method, whereas the 2D-IADIA method has higher image entropy at low SNR levels, suggesting that the 2D-AFCIFSBL method has a certain resistance to noise. In fact, SBL-based methods like 2D-Fast SBL and 2D-AFCIFSBL exhibit consistent performance regardless of SNR variations. In contrast, ADMM-based iterative optimization algorithms, such as 2D-UADN and 2D-IADIA, demonstrate increased image entropy under low SNR conditions. In terms of processing time, as shown in Figure 13b, both the 2D-AFCIFSBL and 2D-Fast SBL methods require more computational time compared to the 2D-IADIA and 2D-UADN methods. This can be attributed to the fact that the 2D-UADN and 2D-IADIA methods are iterative optimization algorithms based on ADMM, which typically exhibit faster computational speeds than sparse Bayesian algorithms like 2D-Fast SBL and 2D-AFCIFSBL. However, compared to the 2D-Fast SBL method, the 2D-AFCIFSBL method has higher computational efficiency in accomplishing the reconstruction task. It demonstrates that the 2D-AFCIFSBL method can improve computational efficiency without compromising imaging performance.

Finally, we analyzed the performance of different algorithms in estimating phase errors, as shown in Figure 14. The range of variance of the phase error is set to $0-\pi$, and the 2D SPR is set to $(0.8,0.8)$. Figure 14a shows the convergence curve of the estimated phase error with the variance of the phase error set to $\pi /2$. It can be seen that 2D-UADN has a larger NMSE due to its inability to estimate the phase error between subpulses. Both 2D-IADIA and 2D-AFCIFSBL have the ability to estimate 2D phase errors; therefore, the NMSE is relatively small. Figure 14b shows the impact of the variance of phase error on the entropy of ISAR images. Figure 14c shows the impact of the variance of phase error on the phase error estimation. As the phase error increases, 2D-UADN significantly deteriorates due to its inability to estimate the phase error between subpulses. Figure 14d shows the impact of the variance of phase error on runing time.

4.2. Experimental Results on Measured Dataset

The measured dataset of the Yak-42 aircraft is collected by an X-band radar which transmits LFM signals with a central frequency of 10 GHz and a bandwidth of 400 MHz [30]. The complete radar echoes consist of 256 pulses with each pulse comprising 256 samples. After traditional range alignment and phase compensation, the received echo dimensions are $256\times 256$. However, these echo data are not a stepped frequency signal, so we obtain the equivalent stepped frequency signal by conducting the inverse Fourier transform on the range dimension of the original data. Typically, the linear stepped frequency signal can be considered a specific instance of a random stepped frequency signal. Therefore, we utilize the Yak-42 aircraft measured dataset to validate the imaging performance of the proposed 2D-AFCIFSBL method in this subsection. In this experiment, we only compared the imaging results under the random phase error, where the variance of the amplitude is $\pi /2$. In addition, we set the SNR of the experiment to 20 dB by adding complex-valued Gaussian noise into measured data.

First, the real image of the Yak-42 aircraft and its ISAR imaging result acquired by the 2D range-Doppler method with a complete measured dataset are demonstrated in Figure 15a and Figure 15b, respectively. It is evident from Figure 15b that the phase error that is not fully compensated can lead to severe blurring in the imaging results. In order to acquire sparse data with the SPR of $(0.8,0.8)$, $(0.6,0.6)$ and $(0.4,0.4)$, we randomly selected few subpulses and bursts from the complete echo data. And the imaging results acquired by the aforementioned four methods are demonstrated in Figure 16, in which three rows from top to bottom represent the imaging results under the SRP of $(0.8,0.8)$, $(0.6,0.6)$ and $(0.4,0.4)$, respectively, and the imaging results of each column are acquired by the corresponding aforementioned methods.

It is evident from Figure 16 that as the SPR decreases, the completeness of the aircraft’s imaging gradually diminishes, and the image of the aircraft’s nose part gradually disappears, which also causes the TBR of the algorithms to gradually decrease, as shown in

Table 2. Numerical evaluation results of measured dataset of Yak-42 aircraft under different SPRs.

SPR	Methods	IE	TRB	Time (s)
$(0.8,0.8)$	2D-Fast SBL	6.01	5.55	1.54
	2D-UADN	9.14	−2.96	0.93
	2D-IADIA	5.05	22.29	0.94
	2D-AFCIFSBL	4.86	28.27	1.17
$(0.6,0.6)$	2D-Fast SBL	5.22	7.85	1.48
	2D-UADN	8.62	−3.05	0.91
	2D-IADIA	4.81	20.38	0.9
	2D-AFCIFSBL	4.68	29.29	1.1
$(0.4,0.4)$	2D-Fast SBL	5.35	4.86	1.6
	2D-UADN	7.93	−2.23	0.94
	2D-IADIA	4.55	17.42	0.92
	2D-AFCIFSBL	4.51	19.18	1.07

. In addition, the 2D-UADN method exhibits the poorest imaging performance with numerous pseudo-scatterers and noise floor across all SPR cases compared to other methods. The 2D-Fast SBL method has slightly better imaging performance than the 2D-UADN method due to its sparse prior resistance to noise. Further, the 2D-IADIA method becomes almost ineffective in case the SPR is $(0.4,0.4)$ due to the appearance of noise floor, while the 2D-AFCIFSBL method continues to produce neat and well-focused images. Table 2 summarizes the experimental results of three evaluation indicators under three types of SPR. Through the analysis of the evaluation indicators in the table, it is evident that the 2D-AFCIFSBL method surpasses the 2D-Fast SBL and 2D-UADN methods in terms of image entropy and TBR by a large margin. The 2D-IADIA method has a similar image entropy to the 2D-AFCIFSBL method, but there is still a gap in TRB. Regarding computational time, the processing time of the four methods on the measured dataset is similar to that on the simulated dataset, as the dimensions of both datasets are the same. Additionally, the imaging speed of the 2D-AFCIFSBL method is observed to be slower than that of the 2D-IADIA method, because the iterative optimization algorithms based on ADMM are faster than the sparse Bayesian learning algorithms. In SBL-based methods, the 2D-AFCIFSBL method is faster than the 2D-Fast SBL method. Nonetheless, the running times of the 2D-AFCIFSBL method remain within the same order of magnitude as the iterative optimization algorithms.

The final experiment was conducted with the SPR of $(0.7,0.7)$ under different SNR conditions, as illustrated in Figure 17, in which three rows from top to bottom represent the imaging results under 20 dB, 10 dB and 0 dB SNR, respectively. It is evident that as the SNR decreases, the imaging results of the 2D-UADN method deteriorate sharply, while the imaging effect of the 2D-Fast SBL method is not as bad as that of the 2D-UADN method due to the sparse prior resistance to noise. For the 2D-IADIA method, as the SNR decreases, many image details are lost, and instead some noise floor begins to appear in the imaging result, which leads to a decrease in the TBR, as shown in

Table 3. Numerical evaluation results of measured dataset of Yak-42 aircraft under different SNRs.

SNR	Method	IE	TRB	Time (s)
20 dB	2D-Fast SBL	5.16	10.32	1.38
	2D-UADN	8.70	−2.33	0.93
	2D-IADIA	5.02	19.99	0.91
	2D-AFCIFSBL	5.01	23.94	1.08
10 dB	2D-Fast SBL	5.11	10.6	1.46
	2D-UADN	8.75	−2.72	0.90
	2D-IADIA	5.01	16.51	0.98
	2D-AFCIFSBL	4.97	24.47	1.34
0 dB	2D-Fast SBL	5.10	6.14	1.54
	2D-UADN	9.21	−5.61	0.92
	2D-IADIA	4.91	8.85	0.94
	2D-AFCIFSBL	4.63	10.38	1.22

. It can also be observed that even with a decrease in the SNR, the 2D-AFCIFSBL method maintains lower image entropy and a higher TBR compared to other methods. Regarding the running time, each method exhibits similar performance under different SNRs.

5. Conclusions

In this article, we propose a 2D joint imaging and autofocusing method named 2D-AFCIFSBL for SRSF ISAR which can process the 2D echoes directly through matrix operations. First, we model SRSF ISAR imaging using sparse Bayesian learning by imposing a gamma-complex Gaussian hierarchical prior, thereby ensuring the sparsity of the probabilistic imaging model. Then, we increase the speed of Bayesian inference by relaxing the evidence lower bound (ELBO) to avoid matrix inversion, and we further convert the iterative process into a matrix form to improve computational efficiency. After that, phase error estimation is integrated into the ISAR imaging stage, which can be calculated through maximum likelihood estimation. Experimental results on both the simulated and measured datasets have demonstrated that even under low SNR conditions, well-focused images can be achieved by the proposed 2D-AFCIFSBL method in high computational efficiency. In the future, we will concentrate on devising probabilistic models that better align with real-world scenarios, leverage the advantages of deep learning methodologies to find more precise solutions, and explore the potential combination between sparse Bayesian inference and deep learning techniques. Additionally, this paper does not consider the structural information of weak scattering points. We will utilize block sparse Bayesian learning to extract structured information from weak scattering regions in subsequent work.