Radar Imaging by Sparse Optimization Incorporating MRF Clustering Prior

TL;DR

This paper introduces a novel radar imaging method that combines sparse optimization with MRF clustering priors, improving image reconstruction by exploiting clustered sparse structures.

Contribution

It proposes integrating an MRF prior with FISTA for enhanced sparse signal reconstruction in radar imaging, which is a new approach in this context.

Findings

Outperforms standard FISTA in image quality.

Demonstrates advantages over existing MRF-based methods.

Validated through simulations and experiments.

Abstract

Recent progress in compressive sensing states the importance of exploiting intrinsic structures in sparse signal reconstruction. In this letter, we propose a Markov random field (MRF) prior in conjunction with fast iterative shrinkagethresholding algorithm (FISTA) for image reconstruction. The MRF prior is used to represent the support of sparse signals with clustered nonzero coefficients. The proposed approach is applied to the inverse synthetic aperture radar (ISAR) imaging problem. Simulations and experimental results are provided to demonstrate the performance advantages of this approach in comparison with the standard FISTA and existing MRF-based methods.

Figures

a that the MRF-FISTA improves clearly over FISTA and LaSAL2 more than at least 3 dB and 4 dB, respectively, when the data ratios are above 20%. Also, the entropy of the results obtained by the MRF-FISTA is much smaller than those of FISTA and LaSAL2, indicating that the imaging results of MRF-FISTA are best focused among these three algorithms, as demonstrated in Fig. [7](#S4.F7)LABEL:sub@fig7b.

Tables

Table 1. TABLE I: Simulation Parameters for Full Data Set

Parameters	Values
Distance $R_0$	5 km
Start frequency	35 GHz
Stop frequency	36 GHz
Number of frequency steps	64
Total rotation angle	${1.7}^{\circ }$
Number of angle steps	64

Equations

$$\mathbf{y}=\mathbf{\Phi x}+ϵ$$

$$\hat{\mathbf{x}}=\arg \underset{\mathbf{x}}{\min }\left\{\frac{1}{2}\Vert \mathbf{y}-\mathbf{\Phi x}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1\right\}$$

$$\hat{\mathbf{s}}$$

$$=\arg \underset{\mathbf{s}}{\max }p(\mathbf{x}|\mathbf{s})p(\mathbf{s}).$$

$$p(\mathbf{s};\alpha,\beta)\!=\!\frac{1}{Z}e^{-\frac{1}{T}\big{[}\sum_{i\in{V}}V_{1}(s_{i})\!+\!\sum_{(i,j)\in{E}}{V_{2}(s_{i},s_{j})}\big{]}}$$

$$\displaystyle\mathbf{\hat{s}}\!=\!\arg\max\limits_{\mathbf{s}}\bigg{\{}\!\!-\!\!\sum_{(i,j)\in{E}}V_{2}(s_{i},s_{j})\!\!-\!\!\sum_{i\in{V}}\big{[}V_{1}(s_{i})\!\!+\!\!\log(p(x_{i}|s_{i}))\big{]}\bigg{\}}.$$

$$\hat{\mathbf{x}}$$

$$=\arg \underset{\mathbf{x}}{\max }p(\mathbf{y}|\mathbf{x})p(\mathbf{s}|\mathbf{y},\mathbf{x})p(\mathbf{x})=\arg \underset{\mathbf{x}}{\max }p(\mathbf{y}|\mathbf{x})p(\mathbf{s}|\mathbf{x})p(\mathbf{x})$$

$$=\arg \underset{\mathbf{x}}{\max }p(\mathbf{y}|\mathbf{x})p(\mathbf{x}|\mathbf{s})p(\mathbf{s})=\arg \underset{\mathbf{x}}{\max }p(\mathbf{y}|\mathbf{x})p(\mathbf{x}|\mathbf{s}).$$

$${\hat{\mathbf{x}}}^{\prime }$$

$$\hat{\mathbf{x}}$$

$${\hat{\mathbf{x}}}^{\prime }=\arg \underset{{\mathbf{x}}^{\prime }}{\min }\left\{\frac{1}{2}\Vert \mathbf{y}-\mathbf{\Phi }{\mathbf{x}}^{\prime }{\Vert }_2^2+{\lambda }_1\Vert {\mathbf{x}}^{\prime }{\Vert }_1+{\lambda }_{\text{TV}}\Vert {\mathbf{x}}^{\prime }{\Vert }_{\mathrm{TV}}\right\}$$

$${\hat{\mathbf{x}}}^{\prime }=\arg \underset{{\mathbf{x}}^{\prime }}{\min }\left\{\frac{L}{2}\Vert {\mathbf{x}}^{\prime }-{\mathbf{v}}_k{\Vert }_2^2+{\lambda }_1\Vert {\mathbf{x}}^{\prime }{\Vert }_1+{\lambda }_{\text{TV}}\Vert {\mathbf{x}}^{\prime }{\Vert }_{\mathrm{TV}}\right\}$$

$${\hat{\mathbf{x}}}_1^{\prime }=\arg \underset{{\mathbf{x}}^{\prime }}{\min }\left\{\frac{w_{1}L}{2}\Vert {\mathbf{x}}^{\prime }-{\mathbf{v}}_k{\Vert }_2^2+{\lambda }_1\Vert {\mathbf{x}}^{\prime }{\Vert }_1\right\},$$

$${\hat{\mathbf{x}}}_2^{\prime }=\arg \underset{{\mathbf{x}}^{\prime }}{\min }\left\{\frac{w_{2}L}{2}\Vert {\mathbf{x}}^{\prime }-{\mathbf{v}}_k{\Vert }_2^2+{\lambda }_{\text{TV}}\Vert {\mathbf{x}}^{\prime }{\Vert }_{\mathrm{TV}}\right\},$$

$$Y(f_n,t_m)=\sum _{i=1}^I{\sigma }_i\exp \left[-j2\pi f_n\frac{2R_i(t_m)}{c}\right],$$

$$R_i(t_m)=\sqrt{R_0^2+x_i^2+y_i^2+2R_0[-x_i\sin \omega t_m+y_i\cos \omega t_m]},$$

Full text

I Introduction

The theory of compressive sensing (CS) [1, 2] has generated great interest in radar imaging [3, 4, 5, 6, 7, 8, 9, 10]. Compressive sensing-based radar imaging methods can achieve better image quality than traditional Fourier transform-based techniques under image sparsity condition and with incomplete observations, as viewed from Nyquist sampling requirements.

However, in many real world scenarios, there exists other important signal structure information, beside sparsity, which can be exploited for enhanced imaging [11, 12, 13, 14]. For instance, most of the wavelet coefficients of a natural image are both sparse and exhibit a tree structure [11]. Clustering is a structure property in which the image assumes high values in neighboring pixels. It may be exhibited in the image canonical basis or under sparsifying basis. From sparse reconstruction perspective, this property is known as block sparse, or clustered sparse[12]. In essence, the support of the image spatial pixels assumes high correlations among closely separated pixels which conforms to the block sparse model.

Image reconstruction approaches incorporating block sparse models were developed in [15, 16]. These approaches combine the $\ell_{1}$ and $\ell_{2}$ norms for feature selections, and assume knowledge of group partition information. This assumption is not mandated in the block sparse Bayesian frame work proposed in [17], though performance could still be affected by the block size.

Using graph theory, the radar image scene can be cast as an undirected graphical model, also referred to as Markov random field (MRF). MRF provides an effective way for modeling spatial context dependent entities, such as image pixels and correlated features. It has been widely employed in a variety of areas including computer vision, image segmentation[18], and sparse signal recovery[12, 13, 14].

The use of MRF prior in CS based magnetic resonance imaging was recently introduced in [19]. The authors extended a constrained split augmented Lagrangian shrinkage algorithm (C-SALSA) with an MRF prior [20], and solved the problem using the alternating direction method of multipliers (ADMM) [21]. The analysis of the ADMM reveals the sensitivity of this type of algorithms to the penalty parameter underpinning the augmented Lagrangian where the primal and dual residuals exhibit conflicting requirements.

In this letter, we incorporate the MRF prior with fast iterative shrinkage/thresholding algorithm (FISTA)[22]. Unlike ADMM, the penalty parameter of FISTA is used to balance the twin objectives of minimizing both the $l_{2}$ error and signal sparsity, without significantly affecting the algorithm convergence. FISTA also preserves the computational simplicity of ISTA family, and achieves a quadratic convergence rate. To demonstrate performance, we apply the proposed MRF-based FISTA algorithms to inverse synthetic aperture radar (ISAR) imaging, in which scatterers are typically located contiguously, satisfying the clustering property.

The remainder of this letter is as follows. Section II presents the formulation of the MRF-based FISTA. In Section III, we apply the MRF-FISTA to ISAR imaging. Numerical simulations and experimental results are shown in Section IV. Finally, concluding remarks are presented in Section V.

II MRF-based FISTA

In compressive sensing, we aim to recover a signal $\mathbf{x}\in\mathbb{C}^{\mathnormal{N}}$ from $M\leq N$ noisy linear measurements $\mathbf{y}\in\mathbb{C}^{\mathnormal{M}}$,

[TABLE]

where $\mathbf{\Phi}\in\mathbb{C}^{\mathnormal{M}\times\mathnormal{N}}$ denotes a known sensing matrix, and $\mathbf{\epsilon}\in\mathbb{C}^{\mathnormal{M}}$ is white Gaussian noise. In the following analysis, we assume $\mathbf{x}$ is sparse in its canonical basis for convenience. If not, we can apply a sparsifying transform $\mathbf{\Psi}$, such as $\mathbf{x}\!=\!\mathbf{\Psi}\bm{\alpha}$, where $\bm{\alpha}$ is a sparse vector.

The estimation of $\mathbf{x}$ from $\mathbf{y}$ in (1) is an ill-posed linear inverse problem, since the sensing matrix $\mathbf{\Phi}$ is singular or ill-conditioned. However, if $\mathbf{x}$ is sparse or has sparse coefficients in some basis $\mathbf{\Psi}$, then exact solution of $\mathbf{x}$ can be obtained by solving the following relaxed convex optimization problem,

[TABLE]

where $\|\mathbf{x}\|_{1}=\sum_{i}|{x}_{i}|$ denotes the $\ell_{1}$ norm of $\mathbf{x}$, and $\|\mathbf{y}\|_{2}=(\sum_{i}|{y}_{i}|^{2})^{1/2}$ represents the $\ell_{2}$ norm of $\mathbf{y}$.

To invoke the clustering property, we assume the probability distribution of the prior structural information as $p(\mathbf{s})$, with $\mathbf{s}$ representing the signal support. The relationship among the signal $\mathbf{x}$, the signal support $\mathbf{s}$, and the measurement $\mathbf{y}$ can be illustrated by an undirected graphical model, shown in Fig. 1. Independence between two variables is displayed as a lack of connection between their corresponding vertices in the graph; conversely, dependent variables should correspond to connected vertices. Further, the absence of a direct link between two variables reflects an interaction that is conditional. In this regard, and for the underlying problem, given signal $\mathbf{x}$, the signal support $\mathbf{s}$ and the measurement $\mathbf{y}$ are independent. Accordingly, the maximum a posterior (MAP) estimation of the signal support $\mathbf{s}$ becomes

[TABLE]

The aforementioned undirected graph, referred to as Markov random field, can be described by an Ising model [19],

[TABLE]

where each vertex has two possible states $\mathbf{s}\in\{0,1\}^{N}$, $Z$ and $T$ are constants. The content in the “[]” is the energy function, which is a sum of potentials over the single-site cliques $V$ and the pair-site cliques $E$, defined as $V_{1}(s)\!=\!\pm\alpha$ corresponding to $s\!=\!0$ and $1$, respectively, and $V_{2}(s_{i},s_{j})\!=\!\pm\beta$ corresponding to $s_{i}\!\neq\!s_{j}$ and $s_{i}\!=\!s_{j}$, respectively. It is noted that a higher value of $\alpha$ enforces a sparser signal activity, and a higher value of $\beta$ implies a stronger spatial correlation.

Substituting (4) into (II) yields,

[TABLE]

In compressive sensing, $p(\mathbf{x})$ is typically selected as Laplacian $p(\mathbf{x})=\frac{\lambda}{2}\exp\left(-\frac{\lambda}{2}\|{\mathbf{x}}\|_{1}\right)$ to preserve sparsity.

Exact inference on (5) is hard because the graph in Fig. 1 is not a tree, so one cannot use any of the sum-product algorithms. However, the variables of the Ising model can be estimated using approximate inference algorithms [23]. In this letter, we use the Markov Chain Monte Carlo (MCMC) sampler introduced in [18].

Upon obtaining the intermediate value of $\mathbf{s}$, the MAP estimation of $\mathbf{x}$ can proceed as,

[TABLE]

For white Gaussian noise with variance $\sigma^{2}$, and $p(\mathbf{x})$ as Laplacian, (II) can be represented by the following two steps,

[TABLE]

where $\mathbf{\hat{x}}^{\prime}$ denotes the recovered signal without MRF prior, “$\circ$” represents an element-to-element multiplication. Therefore, we can obtain an estimate of signal $\mathbf{x}$ from the solutions of (5) and (II). Towards this end, we employ convex fast iterative shrinkage/thresholding algorithm (FISTA) [22] combined with the MRF prior. FISTA has been shown to have a quadratic convergence rate, whereas the traditional ISTA family only achieves a first-order one.

In this letter, we apply two norms, namely, the $\ell_{1}$ norm prior and the TV norm prior, to improve the recovered signal quality. Whereas the $\ell_{1}$ norm is typically used to enhance the sparsity of the signal, the TV norm is utilized as an edge preserving regularization of the signal. We consider the signal reconstruction problem,

[TABLE]

where the parameter $\sigma^{2}$ in (II) is absorbed into $\lambda_{1}$ and $\lambda_{\textrm{TV}}$.

Based on approximate Taylor expansion [22], (8) can be transformed to,

[TABLE]

where $\mathbf{v}_{k}=\mathbf{x^{\prime}}_{k-1}-\frac{1}{L}\nabla f(\mathbf{x^{\prime}}_{k-1})$, $\mathbf{x}^{\prime}_{k}$ represents the estimated value of $\mathbf{x}^{\prime}$ at the $k\textrm{th}$ iteration, $\nabla f(\mathbf{x^{\prime}})=\mathbf{\Phi}^{H}(\mathbf{\Phi x^{\prime}}-\mathbf{y})$ representing the gradient of $f(\mathbf{x^{\prime}})$, $L$ is the Lipschitz constant of $\nabla f$, and $f(\mathbf{x^{\prime}})=\frac{1}{2}\|{\mathbf{y}-\mathbf{\Phi x^{\prime}}}\|_{2}^{2}$.

To solve for $\mathbf{x}^{\prime}$, we apply FISTA. In particular, we employ the composite splitting approach in [24] to decompose the above problem into two sub problems. The first problem is

[TABLE]

and the second is

[TABLE]

where $w_{1}$ and $w_{2}$ are the decomposition weights that separate $\frac{L}{2}\|{\mathbf{x^{\prime}}-\mathbf{v}_{k}}\|_{2}^{2}$ into two groups associated with the $\ell_{1}$ norm and the TV norm, respectively, and satisfying $w_{1}+w_{2}=1$.

Equations (10) and (11) can be readily solved by the algorithms in [22]. The final solution is $\mathbf{\hat{x}}^{\prime}=w_{1}\mathbf{\hat{x}}_{1}^{\prime}+w_{2}\mathbf{\hat{x}}_{2}^{\prime}$. The different steps in the proposed MRF based FISTA algorithm are given in Algorithm 1. Here, we assume $w_{1}\!=\!w_{2}\!=\!1/2$.

III MRF-FISTA Applied to ISAR Imaging

The ISAR geometry, after compensating for the translational motion, can be shown in Fig. 2 in which the target rotates around its centroid with a constant angular velocity $\omega$. The scattered electromagnetic (EM) wave after demodulation is,

[TABLE]

where $\sigma_{i}$ denotes the backscattering coefficient of the $i\mathrm{th}$ point scatterer, $I$ is the number of total pixels of the target scene, $f_{n}$ represents the operating frequency: $f_{n}\!=\!f_{0}\!+\!n\Delta f$, $f_{0}$ and $\Delta f$ denote the start frequency and the frequency step of the transmitted signal, respectively, and $n$ is an integer from the set $[0,N-1]$. The variable $t_{m}$ represents the discrete slow time in the cross-range dimension, with $m$ is an integer from the set $[0,M-1]$, $c$ is the speed of EM wave propagating in free space, and $R_{i}(t_{m})$ is given by,

[TABLE]

representing the instantaneous range between the $i\mathrm{th}$ scatterer and the radar, where $R_{0}$ denotes the distance from the target centroid to the radar.

According to (12), we can formulate the element of the sensing matrix $\mathbf{\Phi}$ as, \phi_{l,i}=\exp\big{[}\mathrm{-j2\pi}f_{n}\frac{2R_{i}(t_{m})}{c}\big{]}, where $n=(l\mod N)$, and $m=\lfloor\frac{l}{M}\rfloor$. Here, “mod” means the modulo operation, and $``\lfloor\rfloor"$ denotes rounding toward negative infinity, and $l=0,1,\cdots,MN-1$, $i=0,1,\cdots,I-1$.

The target scene can be viewed as a matrix $\mathbf{X}$ with a size of $N\times M$ (assuming $I=N\times M$). Thus, the forward model of radar imaging can be expressed as, $\textrm{vec}(\mathbf{Y})=\mathbf{\Phi}\textrm{vec}(\mathbf{X})$, where $\mathbf{Y}$ denotes the matrix form of the received data (with zeros padding at the positions where there are no samples), and $\textrm{vec}(\cdot)$ represents vectorization of a matrix. Then, the image $\mathbf{X}$ can be reconstructed by the proposed MRF-FISTA.

If the size of the sensing matrix is too large, we can use the method in [25] to change $\mathbf{\Phi}$ into a 2-D Fourier transform with a compensation for the spherical wavefront. In addition, if the distance $R_{0}$ is sufficiently large, (13) can be approximated as $R_{i}(t_{m})\!\approx\!R_{0}\!-\!x_{i}\sin{\omega t_{m}}\!+\!y_{i}\cos{\omega t_{m}}$. In this case, the sensing matrix can also be replaced by a Fourier transform.

IV Simulations and Experimental Results

This section shows the performance of the MRF-FISTA based ISAR imaging method using both simulation and real data.

IV-A Simulation results

The simulation parameters for the full data set, according to the Nyquist-Shannon sampling theorem, are shown in Table I. The target model is shown in Fig. 3.

First, we present the imaging results of four techniques, two of which use MRF prior. We consider the proposed MRF-FISTA, FISTA, LaSAL2, which is a recently introduced ADMM based algorithm with the MRF prior[19], and the traditional back-projection (BP) algorithm, respectively, in Fig. 4. The sampling rate (i.e., the ratio of the utilized data size over the full data set size) is 30%. The signal-to-noise ratio (SNR) is 20 dB. Note that the imaging result of MRF-FISTA is vividly of higher quality than those of the other three algorithms.

In order to show the performance sensitivity of the proposed method and LaSAL2 with respect to their corresponding penalty parameters, we compare the convergence of the algorithms with different penalties. Here, we employ the root mean squared error (RMSE) to assess convergence, which is given by, \textrm{RMSE}\!\!=\!\!\sqrt{\frac{1}{NM}\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}\big{[}\mathbf{X}(n,m)\!\!-\!\!\mathbf{\hat{X}}(n,m)\big{]}^{2}}, where $\mathbf{X}$ and $\mathbf{\hat{X}}$ denote the target model and the recovered image, respectively.

Fig. 5 shows the RMSEs of LaSAL2 with different starting value of the penalty parameters, $\mu$, which underlies the augmented Lagrangian, for different data sets. The SNR is 20 dB for this simulation. Note that smaller penalty parameters can result in divergence for larger data sets. A larger penalty has a better solution, however, with a slower convergence rate. Thus, for such methods, it becomes critical to choose a proper starting penalty parameter.

Unlike the ADMM based algorithms, the penalty parameter $\lambda$ of FISTA is only employed to tradeoff the $l_{2}$ error and the signal sparsity.

Fig. 6 shows the RMSEs of MRF-FISTA with different $\lambda$s. Clearly, the convergence rate almost remains invariant with different penalty parameters.

Next, we include another performance metric – entropy, for the quantitative comparisons. The Shannon entropy of an image $\mathbf{X}$ is defined as, $H(\mathbf{X})=-\sum_{k=1}^{K}p_{k}(\mathbf{X}_{k})\log_{2}p_{k}(\mathbf{X}_{k})$, where $(p_{1},p_{2},\cdots,p_{K})$ is a finite discrete probability distribution, and $p_{k}$ contains the normalized histogram counts within a fixed bin $\mathbf{X}_{k}$, i.e., it is a function of pixel intensity in an image. The distribution of pixel intensity is commensurate with the degree of image focus.

Fig. 7 shows the RMSEs and entropy of the imaging results obtained by the MRF-FISTA, FISTA, and the LaSAL2, respectively, for different data ratios. The SNR is fixed to 10 dB. We perform 50 independent runs for each data ratio to obtain the mean values of the RMSEs and entropy. It is noted from Fig. 7LABEL:sub@fig7a that the MRF-FISTA improves clearly over FISTA and LaSAL2 more than at least 3 dB and 4 dB, respectively, when the data ratios are above 20%. Also, the entropy of the results obtained by the MRF-FISTA is much smaller than those of FISTA and LaSAL2, indicating that the imaging results of MRF-FISTA are best focused among these three algorithms, as demonstrated in Fig. 7LABEL:sub@fig7b.

The parameters of the Ising model were chosen as: $\alpha=0.01$ and $\beta=0.3$ for the above simulations.

IV-B Experimental results

Here, we use measurements from an unmanned aerial vehicle. The radar frequencies in the experiments vary from 34.5 GHz to 35.3 GHz with 128 steps. The total observation angles are 2 degrees, also with 128 steps.

Fig. 8 depicts the imaging results of MRF-FISTA, FISTA, LaSAL2 and BP, using 30% of the full sampled data. It is seen that the result of MRF-FISTA assumes the best quality with much smaller noisy like clutter among all the algorithms.

Since we do not have a bench mark for the real data imaging, we present the entropy results of the images to assess the quantitative performance of the algorithms when processing the real data. This is shown in Fig. 9. Clearly, the precision of MRF-FISTA is still higher than those of FISTA and LaSAL2, with much smaller entropy.

V Conclusion

This letter proposed a MRF based FISTA as an effective image reconstruction technique for cluster-type images, including those of ISAR. The MRF is a powerful tool for modeling spatial context dependent entities. It was used in this letter as a prior to describe the high spatial correlations in an image. The MRF prior was incorporated into sparsity based optimization algorithm, and shown to yield higher image reconstruction precision compared to existing techniques. Simulations and experimental results demonstrated that the proposed algorithm improves over the original FISTA and another recently introduced MRF based ADMM algorithm when considering imaging errors as well as convergence rate.