MF-JMoDL-Net: A Deep Network for Azimuth Undersampling Pattern Design and Ambiguity Suppression for Sparse SAR Imaging

We assume that the SAR system operates in single-channel strip mode. The geometric configuration of a typical side-looking SAR is shown in Fig.1 Suppose that the platform carrying the imaging radar moves in a straight line relative to the observation scene at a constant speed $v$. For stripmap SAR, the direction of the antenna beam remains fixed throughout the entire movement. Consequently, the distance between the phase center of the antenna and the target point with coordinates $(x,y)$ in the observation scene is

where $\tau$ and $\eta$ are the fast time and the slow time, respectively, $x$ and $y$ are the spatial coordinates of the target point along the azimuth and range directions, and $R_{0}$ is the closest range between target and platform.

The signal transmitted by the antenna to the observation scene is a linear frequency modulated signal (LFM), therefore, the echo received by the imaging radar of the target point $(x,y)$ in the observation scene is

where $A_{0}$ is the observation gain determined by the antenna; $\sigma\left({x,y}\right)$ is the backscattering coefficient of the target point $(x,y)$; $\omega_{a}(\cdot)$ and $\omega_{r}(\cdot)$ are the azimuth and range signal envelope, respectively; $f_{0}$ is the carrier frequency; $c$ is the speed of light. The SAR echo reflected from the scatterers generated by the antenna beam coverage area $\Omega$ can be expressed as

where $N(\tau,\eta)$ is the additive measurement noise. By discretizing the observation scene $\Omega$ into $M_{r}\times M_{a}$ grids and sampling the echo signal into a matrix, the discrete echo can be expressed as

where $s\left({\tau_{n_{r}},\eta_{n_{a}}}\right)$ is the $n_{r}$-th range sample at the $n_{a}$-th pulse of raw signal $s\left({\tau,\eta}\right)$ , $M_{r}$ and $M_{a}$ are the numbers of the discrete cells in the range and azimuth directions, respectively. The $h\left[\tau_{n_{r}},\eta_{n_{a}},x_{m_{r}},y_{m_{a}}\right]$ is defined by

Write the discretized data $s\left({\tau_{n_{r}},\eta_{n_{a}}}\right)$ and the reflectivity map $\mathbf{\sigma}\left({x_{m_{r}},y_{m_{a}}}\right)$ as two matrices: $\mathbf{S}\in\mathbb{C}^{{M\times N}_{r}}$ and $\mathbf{\sigma}\in\mathbb{C}^{{M_{a}\times M}_{r}}$, respectively, where $M$ is the number of transmitted pulses in the azimuth direction; $N_{r}$ is the number of samples in the range direction. Then, (4) can be rewritten as

where $\mathbf{n}$ is the vector of additive noise, and $\mathbf{H}$ is defined by

Conventionally, the intention of CS-SAR imaging is to reconstruct $\mathbf{\sigma}$ from the measurements $\mathbf{S}$.

In this subsection, we will introduce and discuss the typical MF-Based CS-SAR imaging model to explain the original intention of the proposed method. Then, we discuss the forward model with the sampling pattern, which is the central assumption of the proposed MF-JMoDL-Based SAR imaging model.

1) Typical MF-Based CS-SAR Imaging Model: in the typical CS-SAR imaging model, the measurement vector $\mathbf{H}$ is generally resampled by an appropriate sampling matrix $\Theta$. When the observation scene is sparse and the observation matrix $\mathbf{A}=\Theta\mathbf{H}$ satisfies the RIP condition [15], $\mathbf{\sigma}$ can be exactly recovered from $\mathbf{S}$ with the $L_{q}(0<q<1)$ optimization:

To solve this problem, an equivalent regularization scheme with the following optimization problem can be used:

where $\lambda$ is a regularization parameter. Generally speaking, the solution methods for (8) and (9) can be divided into two categories: greedy methods and iterative methods. However, regardless of the chosen method, the computational and storage costs are enormous due to the size of the measurement matrix $\mathbf{H}$, which significantly reduces computational efficiency. Hence, the corresponding reconstruction algorithms are considerably more time-consuming than traditional matched filter MF-based focusing methods.

To address this issue, the MF-based sparse SAR imaging model [28, 27] offers a general principle on how the observation can be remodeled and approximated by any high-precision imaging procedure:

where $\mathbf{M}$ is the traditional MF imaging procedure that can be calculated through decoupling it into a series of 1-D operators in the frequency domain, like chirp scaling algorithm (CSA) and range doppler algorithm (RDA), $\mathbf{G}$ is any generalized right inverse of $\mathbf{M}$.

2) Forward Model with Sampling Pattern: To optimize the azimuth sampling scheme of MF-based SAR imaging, we employ the model-based imaging algorithms [40][41]. The model-based imaging schemes use a continuous function of $\mathbf{A}$, denoted by the operator $H_{\Theta_{a}}^{*}(\cdot)$, which maps $\mathbf{\sigma}$ to $\mathbf{S}$, and $\Theta_{a}$ represents the azimuth sampling pattern. Then, we construct the forward model as follows:

Due to the recovery of $\mathbf{\sigma}$ from $\mathbf{S}$ is ill-posed, model-based imaging algorithms pose the recovery as a regularized optimization scheme:

here $N\left(\mathbf{\sigma}\right)$ is a regularization penalty. When $\mathbf{\sigma}$ is an artifact-free image, $N\left(\mathbf{\sigma}\right)$ is a small scalar, otherwise, it is a high scalar for noisy images. Commonly-used regularizers in imaging reconstruction include total variation regularization [42], transform domain sparsity [43], and structured low-rank methods [44]. These regularizers aim to enforce specific desirable properties on the reconstructed image, such as smoothness, sparsity in certain domains, or low-rank structures, ultimately leading to improved image quality and reduced artifacts.

In order to address the challenges posed by the proposed model in 2.2, two main problems must be overcome.

Firstly, the large-scale exact measurement matrix $\mathbf{A}$ presents significant computational and storage requirements, which would reduce the efficiency of the proposed SAR imaging model. To tackle this issue, the MF-based SAR imaging method is employed to construct an approximate measurement operator by decoupling azimuth and range processing. Following a similar approach to the work in [28], the chirp scaling algorithm (CSA) is used to build the approximate measurement operator.

Secondly, the sampling pattern proposed in this paper leads to nonuniform samples in the azimuth direction. Consequently, a nonuniform MF operator $H_{\Theta_{a}}(\cdot)$ is needed, which is different from the regular MF operator. In a similar manner to the work in [31], the chirp scaling algorithm (CSA) is adopted to construct the approximate measurement operator. The process for constructing the approximate measurement operator is illustrated in Fig.2(a) and Fig.2(b).

By addressing these two challenges, the proposed SAR imaging model can be solved efficiently while maintaining the benefits of the optimized azimuth sampling pattern.

The formulation of CSA can be expressed as follows:

In the given context, $\hat{\mathbf{\sigma}}$ represents the reconstructed SAR image matrix, while $\mathbf{F}{\tau}$ and $\mathbf{F}{\eta}$ denote the fast Fourier transform (FFT) along the range direction and azimuth direction, respectively. $\mathbf{IF}{\tau}$ is the inverse fast Fourier transform (IFFT) along the range direction. The phase terms for RCMC, range compression, azimuth-range decoupling, and azimuth compression are represented by filters $\theta{1}$, $\theta_{2}$, and $\theta_{3}$, respectively. The expressions for these phase terms can be found in [31], and as such, they are not repeated here.

The equivalent CSA operator is denoted as $H_{\Theta_{a}}(\cdot)$. Therefore, equation (13) can be expressed as $\hat{\mathbf{\sigma}}=H_{\Theta_{a}}\left(\mathbf{S}\right)$. This formulation allows for efficient processing of the SAR imaging problem while taking into account the optimized azimuth sampling pattern.

Now we defined the nonuniform Fourier transform matrix $\mathbf{NF}_{\eta_{\Theta_{a}}}$ in (10), which is used to transform the raw data into the range–doppler domain as

it consists of $M$ column vectors defined by

In equation (15), $\eta_{m}$ represents the $m$-th sample in the azimuth direction, with $m=1,2,\dots,M$, where $M$ is the number of samples in azimuth. The $n$-th sample in Doppler frequency is given by $f_{\eta}^{(n)}=(n-1)PRF/N_{a}$, where $n=1,2,\dots,N_{a}$, and $N_{a}$ is the number of uniform frequency grids in the Doppler domain. $PRF$ denotes the azimuth pulse repetition frequency.

By applying the CSA operator to the raw data $\mathbf{S}$, we can obtain the reconstructed backscattering coefficients of the scatterers in the observed scene. These coefficients represent the true backscattering coefficients we desire, which are approximately equal to $\mathbf{\sigma}$. The expression is as follows:

From (13), it is evident that the CSA operator is invertible, as mentioned in [31]. By applying the inverse of the CSA operator to the true backscattering coefficients of the scene, we can obtain an approximation of the raw data $\mathbf{S}$. The inverse CSA operator can be explicitly expressed as follows:

where $\mathbf{F}_{\eta}$ is the azimuth FFT operator. $\hat{\mathbf{S}}$ is the approximated echo data matrix generated by $H_{\Theta_{a}}(\cdot)$, then the (17) can be expressed as $\hat{\mathbf{S}}=H_{\Theta_{a}}\left(\mathbf{\sigma}\right)$. the nonuniform inverse Fourier transform ${\mathbf{N}\mathbf{I}\mathbf{F}}_{\eta_{\Theta_{a}}}$ composed of $\mathbf{N}_{\mathbf{a}}$ column vectors defined by

which is consist of $\mathbf{N}_{\mathbf{a}}$ column vectors defined by

where $\eta_{1},\ldots,\eta_{M}$ are obtained according to proposed method.

Consider the model proposed in 2.2, we define the regularizer as $N\left(\mathbf{\sigma}\right)=\lambda\left\|{\mathbf{T}\mathbf{\sigma}}\right\|_{\mathcal{L}_{1}}$ in transform domain sparsity, with $\Phi=\left\{\mathbf{T},\lambda\right\}$ denoting the parameters of the transform and regularization. To denote the dependence on the regularization parameters as well as the sampling pattern of (12), we give the optimization problem of the form

To address (13), we treat the image reconstruction as a regularization optimization problem and employ a deep learning approach rather than traditional fixed prior regularization to learn parameters from data samples with optimized azimuth sampling patterns, which are commonly referred to as data-driven methods. Classical data-driven image reconstruction techniques include direct inversion schemes [45][46] and the model-based deep learning framework (MoDL) [38]. Owing to the substantial reduction in the number of network parameters in MoDL, we focus on this deep learning framework, formulating its image recovery as follows:

where $D_{\Phi}\left(\mathbf{\sigma}\right)$ is the “denoised” version of $\mathbf{\sigma}$. To extract the ambiguity and noise in $\mathbf{\sigma}$, we let $D_{\Phi}$ be a residual learning-based CNN. Similar to [38], setting $\mathbf{\sigma}_{n}+\mathrm{\Delta}\mathbf{\sigma}=\mathbf{\sigma}$, the non-linear mapping $D_{\Phi}\left({\mathbf{\sigma}_{n}+\mathrm{\Delta}\mathbf{\sigma}}\right)$ can be approximated using Taylor series around the $n^{th}$ iterate as

where $J_{n}$ is a Jacobian matrix, the penalty term can be approximated as

The above approximation is only valid in the vicinity of $\mathbf{\sigma}_{n}$, hence we obtain the alternating algorithm that approximate (21):

here, (24) is implemented using a conjugate gradient algorithm:

This iterative algorithm is unrolled to obtain a deep recursive network $M_{\Theta_{a},\Phi}$, where the weights of the CNN blocks and data consistency blocks are shared across iterations, as shown in Fig. 3. Specifically, the solution to (21) is given by

Note that once unrolled, the image reconstruction algorithm is essentially a deep network, shown in Fig. 4. Thus, the main distinction between MoDL and direct-inversion scheme is the structure of network $M_{\Theta_{a},\Phi}$.

In this work, we concentrate on optimizing the azimuth sampling pattern denoted by $\Theta_{a}$ in (11) and the parameter $\Phi$ of the reconstruction algorithm (20) to enhance the quality of the reconstructed images. Drawing inspiration from earlier sampling pattern learning approaches employed in MRI based on compressed sensing algorithms [47] [48], we incorporate these models into the SAR imaging field to optimize the azimuth sampling pattern $\Theta_{a}$, such that

is minimized. Here ${\hat{\sigma}}_{i}\left\{{\Theta_{a},\Phi}\right\}$ is the corresponding reconstruction of training image $\sigma_{i};~{}i=1,\ldots,N$ in the optimization.

To improve the reconstruction performance, this paper uses a MF operator model-based deep learning network (MF-JMoDL-Net) framework to jointly optimize both $C_{\Theta_{a}}$ and $D_{\Phi}$ blocks in the MF-JMoDL framework (21), which refers to [37]. Given training data, this framework jointly learns the sampling pattern $\Theta_{a}$ and the CNN parameters $\Phi$ using

where $M_{\Theta_{a},\Phi}$ denotes a general deep learning network architecture that includes forward model as well as unrolled architectures denoted by (27).

Fig.3 illustrates the framework of the proposed MF-JMoDL, which alternates between CNN blocks $D_{\Phi}$ and data consistency blocks $C_{\Theta_{a}}$ solely dependent on the sampling pattern. We unroll the MF-JMoDL framework for five iterations of alternating minimization to solve (21), as depicted in Fig.4. The implementation process of the forward operator $H_{\Theta_{a}}$ is shown in Fig. 2(b), with the trainable parameters $\Theta_{a}=\left[{k_{y_{1}}\ldots k_{y_{i}}\ldots k_{y_{N_{a}}}}\right]^{T}$. Since the operator $\left({H_{\Theta_{a}}H_{\Theta_{a}}+I}\right)$ is not analytically invertible for complex operators such as SAR echo data, we need to identify a specific gradient algorithm to handle complex forward models.

In contrast to conventional proximal gradient (PG) algorithms [49, 50, 51] which alternate between steepest descent and CNN blocks, this paper chooses the CG algorithm to avoid a high total number of iterations. In comparison, CG blocks provide a faster reduction per iteration by enforcing a more accurate data-consistency constraint at each unfolding step. Without trainable parameters, numerous CG steps can be performed at each unfolding with no memory cost during training. Experiments in [38] demonstrate that the CG strategy offers improved performance compared to PG and also facilitates the easy incorporation of other image regularization techniques. In this paper, we set the data consistency block $C_{\Theta_{a}}$ to 10 iterations using the CG algorithm.

Regarding the CNN block $D_{\Phi}$, we implement a U-Net model [39], which involves four pooling and unpooling layers with 3$\times$3 trainable filters. The parameters of the blocks $D_{\Phi}$ and $C_{\Theta_{a}}$ are optimized to minimize (21). We rely on the automatic differentiation capability of TensorFlow to evaluate the gradient of the cost function with respect to $\Theta_{a}$ and $\Phi$.

In this subsection, we discuss the loss function and backpropagation gradient calculation of the proposed method, which ensure that the optimization training strategy is continuous and differentiable.

To demonstrate the performance of image reconstruction and enhancement provided by MF-JMoDL-Net, experiments are conducted under various conditions, such as different sparsity scenes, different sampling rates, and different sampling patterns. Additionally, the azimuth ambiguity suppression capability of the proposed method is compared with the local azimuth ambiguity-to-signal ratio estimation (LAASR) [4] ambiguity suppression proposed by Long. The design of the azimuth sampling pattern using the proposed method is compared with staggered [10, 11, 12, 14, 13], uniform sampling, and Poisson-disk sampling methods [30].