Compressive SAR Imaging with Joint Sparsity and Local Similarity Exploitation

Compressive sensing-based synthetic aperture radar (SAR) imaging has shown its superior capability in high-resolution image formation. However, most of those works focus on the scenes that can be sparsely represented in fixed spaces. When dealing with complicated scenes, these fixed spaces lack adaptivity in characterizing varied image contents. To solve this problem, a new compressive sensing-based radar imaging approach with adaptive sparse representation is proposed. Specifically, an autoregressive model is introduced to adaptively exploit the structural sparsity of an image. In addition, similarity among pixels is integrated into the autoregressive model to further promote the capability and thus an adaptive sparse representation facilitated by a weighted autoregressive model is derived. Since the weighted autoregressive model is inherently determined by the unknown image, we propose a joint optimization scheme by iterative SAR imaging and updating of the weighted autoregressive model to solve this problem. Eventually, experimental results demonstrated the validity and generality of the proposed approach.


Introduction
Due to its day/night and all weather performance capabilities, synthetic aperture radar (SAR) has become one of the most promising remote sensing tools in military and civilian fields, including target recognition, topographic mapping, and environmental monitoring. SAR is capable of producing high-resolution images of stationary surface targets and terrain reflectivity. In general, to achieve high-resolution performance, a wideband transmitted signal as well as a large synthetic aperture size are required. However, such requirements lead to a high sampling rate in both the range and azimuth dimensions, which poses a challenge to the analog-to-digital (AD) converter at the receiver and makes the consequent processing complex.
Based on the assumption that SAR images can be sparsely represented in some spaces, such as wavelet, discrete cosine transform (DCT), or Fourier domain, the recently emerged Compressed Sensing (CS) [1,2] theory demonstrates that the SAR image can be exactly recovered with high probability from very limited measurements by solving a convex optimization. Motivated by this, many works [3][4][5][6] have applied the CS theory into the SAR imaging and developed a new compressive SAR imaging approach (CS-SAR for short), in which a series of advantages, such as high resolution and low sidelobes, are provided. In addition, considering that SAR imaging is inevitably disturbed by uncertain phase errors in practice, sparsity-driven approaches [7,8] are proposed for joint SAR imaging and phase error correction. Though the superior capability of the compressive SAR imaging approach has been extensively studied, only simple scenes with isolated scatterers existing in a low-reflectivity surrounding are considered. Previous studies [9,10] were capable of dealing with more general cases and demonstrated that a satisfactory imaging performance can be achieved when a SAR image can be sparsely represented in an elaborately selected space. However, a fixed space lacks adaptivity in characterizing the varying structure of the scene. Especially when the waveforms vary significantly across the whole image, reconstruction with some fixed space may suffer from an unsatisfactory performance.
In this paper, an adaptive sparse representation scheme derived from the weighted autoregressive (AR) model is introduced in compressive SAR imaging. The work is inspired by [11], where the AR model is integrated into the compressive image recovery for exploiting the structural sparsity of natural images. The premise of the work is that the piecewise statistical stationary assumption holds. However, a more general case, which statistically in a local window is non-stationary, is considered in our paper. Specifically, considering the local non-stationary property within the SAR image, the similarity between any two pixels, evaluated via the structural similarity probability (SSP), is incorporated into the AR model to further promote the capability of the structural sparsity exploration. With this weighted AR model, a new approach of compressive SAR imaging with adaptive sparse representation (ASR-CS-SAR for short) is developed. Since the parameters of the weighted AR model are inherently determined by the unknown SAR image, an alternative optimization scheme is introduced to join SAR imaging and parameters estimation. In the end, experimental results demonstrate that the proposed ASR-CS-SAR approach outperforms general CS-SAR approaches with fixed spaces in terms of both visual quality and other metrics.
The remainder of this paper is organized as follows. Section 2 gives a brief review of the spotlight model SAR. The general CS-SAR approach is presented in Section 3 and the proposed ASR-CS-SAR approach is given in Section 4. In Section 5, simulation results are presented to demonstrate the validity of the proposed approach. Finally, conclusions are made in Section 6.

SAR Imaging Model
In this section, we give a brief review of the spotlight mode SAR. The radar traverses along a straight path with a constant speed, and continuously steers the antenna beam to a fixed ground patch of interest. At each location, radar transmits electromagnetic pulses and records the return signal from the ground patch. Then, the collected SAR echoes from multiple observation locations are processed to image the reflectivity profile. In spotlight mode SAR, the most commonly used pulse is p t e ω = s (1) where p(t) is the transmitted waveform, t represents the signal time, and ω0 is the carrier frequency. Usually, chirp signal is adopted as transmitted waveform where rect(·) denotes a rectangular window, γ is the chirp rate, and T represents the signal duration. Supposing the whole scene consists of a series of point scatterers on a grid, the scattering response can be approximated as a sum of the responses from individual point scatterers. As a consequence, the corresponding echo signal at an aspect angle θ after demodulation can be described as where S denotes the illuminated scene, Q(x,y) represents the reflectivity coefficient of the point scatterer at coordinates (x,y), R(θ,x,y) is the distance between the radar and the point scatterer, and c is the speed of light. Note that the SAR image process is usually operated in the discrete domain, which includes two aspects. On the one hand, the illuminated scene is discretized into a Na × Nr grid. Let q be a lexicographic ordered vector of an unknown sampled reflectivity image Q(x,y) of length N = Na × Nr. On the other hand, temporal samplings in both range and azimuth are also discrete. Let rθ(ts) be the s-th sample of rθ(t), we have where R(θ,i) corresponds to the two-way distance between the radar and the i-th point scatterer, and S is the sampling number in range dimension. Considering that SAR imaging is normally synthesized via multiple aspect angles and the real echoes are contaminated by the additive Gaussian noise, the practical SAR echoes can be expressed as where and L is the total number of aspect angle. C serves as the SAR observation matrix with its element . ε denotes the additive noise.
The procedure to obtain the unknown reflectivity q from the echo r is referred to as SAR imaging [12].

CS-SAR Approach
In this section, a brief review of CS theory is given, and then both an efficient undersampling scheme and the mathematical formation of the CS-SAR approach are presented.

Review of CS Theory
Suppose a signal x of size N × 1 is K-sparse in a basis Ψ, then it can be described as where ϑ is the associated coefficient vector and the number of nonzero in ϑ is K with K << N. The recently emerged CS theory indicates that such a signal can be exactly recovered with high probability from M (M < N) linear measurements, and the related observation procedure can be expressed as where Φ is an M × N measurement matrix and y is the measurement vector of length M. Since M < N, solving Equation (7) is highly underdetermined. CS provides a feasible way to solve this problem by exploiting the sparsity of the signal. According to [1,2], provided that the number of measurements satisfies , CS is capable of recovering the sparse signal x (via its coefficient vector ϑ ) from the measurements y by solving the following constrained optimization problem:

CS-SAR Imaging Approach
To apply the CS framework into the SAR imaging, the random undersampling scheme is implemented in both range and azimuth dimensions firstly. Specifically, a pseudorandom sequence {cl, l = 1, 2, ···, L} is randomly generated for undersampling in azimuth. Note that for the case when cl=1, temporal samplers are available for SAR imaging, while for cl = 0, the samplers are discarded. For each cl=1, a similar pseudorandom sequence {cls, s = 1, 2, ···, S} is produced to reduce the sampler number in range. Obviously for any two indices in {cl}, the corresponding range sampling sequence {cls} varies randomly. To reduce the sampling number, the numbers of the nonzero elements in {cl} and {cls} are much smaller than L and S. By representing the pseudorandom sequences in matrix form, we have With Φ, some of the samplers for SAR imaging are discarded, which reduces the sampling rate. By incorporating the sparsity constraint on the scene q, SAR imaging can be transformed into the following minimization problem 1 arg min , . .
Obviously, SAR imaging quality largely depends on the space Ψ. When Ψ is an identity matrix, the point-like features of the unknown scene can be enhanced. When Ψ is a fixed basis, such as DCT or wavelet, different features of the scene can be enhanced [9]. Choosing a fixed and known basis is appealing due to its easy implementation. However, in most cases, local structure features vary across the scene and usually appear irregularly. For instance, different segments of a SAR image have different waveforms, and such varying second-order statistics exhibit varied sparsity in spatial domain. In such cases, the above fixed bases lack adaptivity to characterize the varied image structures and will suffer from limited representation of the signal.

ASR-CS-SAR Approach
In this section, by exploring the structural similarity among pixels, a weighted AR model is obtained. Based on it, a novel compressive SAR imaging approach with adaptive sparse representation is proposed. Finally, the computational complexity of the proposed approach is discussed.

AR Model and Image Structural Similarity
Based on the fact that an image is a non-stationary random Markov field of a modest order, the AR model [11] has shown to be effective in depicting local image structure. Specifically, the pixel qi can be represented as where Ωi is a local neighborhood around the pixel qi, ζi stands for the associated AR parameter, and ωi is a random perturbation independent of spatial location and the image pixel. The subscript iοk denotes the k-th neighbor of pixel qi in the Ωi stored in a raster scan. Generally speaking, using a large Ω is more efficient to characterize image local waveform. However, considering that the AR model is based on the local stationary assumption, using a large Ω may lead to an over-fitting problem. Thus an empirical 3 × 3 local neighbor is used here to achieve the optimum balance between efficiency and robustness. Normally, the calculation of ζi is implemented in a local window Wi where the statistical stationary assumption holds. Therefore, by adjusting ζi to fit the variant local waveform, the local image structure can be illustrated by spatial contiguous pixels. The AR model performs well in a local stationary window. However, in practice, the waveforms vary significantly across the whole SAR image, and the statistical stationary assumption does not always hold [13]. Especially when the scale of a local structure feature is smaller than the selected local window size, the piecewise stationary assumption is violated. Take Figure 1 for example, where Figure 1b is the close-up view of the selected local area in the green window in Figure 1a. Although the statistical property keeps stationary within each ellipse region, it varies significantly across various ellipses. To be specific, to estimate the pixel qi on the edge, samples within the red region (e.g., pixel qa) can provide more reasonable local structure information than those in the yellow (e.g., pixel qb) and blue (e.g., pixel qc) regions. If pixels in the yellow and blue regions are engaged with the same weights, using the AR model in Equation (11) will lead to a degraded description of the local image structure. Considering that any two pixels share similarity, we introduce the structural similarity probability (SSP) as weight into the AR model to further promote its capability of characterizing the local image structure. Specifically, to achieve an adaptive representation of the pixel qi in Figure 1, we prefer to impose large weight on the pixel that shares a similar local structure (e.g., qa), while small weight or nearly zero on the pixel that has large disparity (e.g., pixels qb and qc). Mathematically, the SSP between any two pixels qi and qj is quantitatively evaluated by the local structure vector Zi and Zj, where Zi denotes a square 8-connected neighborhood centered at pixel qi. To be specific, the SSP between pixels qi and qj is modeled as the Gaussian function of the Euclidean distance between their local structure vectors: (12) where h acts as a degree of filtering that controls the shape of the exponential function. To make the structural similarity independent of the image intensity field, the normalized max ⋅ returns the maximum value. For more details, one can refer to [14]. To facilitate understanding, the SSP between any pixel and the center pixel qi in the local window Figure 1b is calculated based on Equation (12) and the associated result is shown in Figure 1c, where the darker intensity grid means significantly lower weight, and vice versa. It is evident that the weights along the edge are nonzero while most of the others in the local window are zero or nearly zero, which indicates that the pixels along the edge can provide more efficient local structure information compared with other pixels. Thus, the profile characterized by the SSP distribution is consistent with the local image structure.

Adaptive Sparse Representation Scheme
To circumvent the risk of data overfitting, the AR model in Equation (11) is normally split into two lower order models as 1 4 where i m  Remembering that similarity among pixels can be used to further promote the capability of describing the structural property, we incorporate the SSP distribution into the AR model, then a weighted AR model can be obtained. Specifically, for a given local window Wi, the parameters αi and βi corresponding to the center pixel qi can be estimated via where χij is calculated by Equation (12 where di is a vector of length N × 1 and the matrix form of Equation (15) can be expressed as = + q Dq ω (16) where and accounts for the residual. As di is determined by the AR parameters, and only 8 entries are nonzero at most, di is a sparse vector.

The Proposed Approach
By imposing the sparsity constraint on each di, a novel compressive SAR image framework with adaptive sparse representation (ASR-CS-SAR for short) can be obtained as 1 min . .
As the structure described by D is intrinsic to the real physical one, the sparse representation derived from the modified AR model is more adaptive than that with a fixed basis. Though the problem in Equation (17) is convex, such a combinational optimization problem still cannot be easily solved since a better estimation of D returns a more accurate representation of the SAR image, which in turn leads to a smaller residual ω. Motivated by this, we replace the regularized term 1 i i  d with the sparsity constraint on ω, and then a relaxed version of Equation (17) can be obtained as Since D greatly depends on the prior information of the image, SAR imaging as well as the estimation of D should be jointly optimized. Thus, an optimization approach for iterative SAR imaging and updating of D is proposed. Specifically, for a given image q, D can be calculated according to Equations (14) and (15). When given D, Equation (18) where λ is the Lagrangian multiplier, and ν and τ are positive scalar parameters. Since q and ω are unknown, the above minimization problem can be solved by alternatively optimizing the following sub-problems, i.e., Clearly, the minimizing of the sub-problems has a close-form formula and thus can be solved effectively. For a fixed image q, we take the gradient of Equation (20) with respect to ω, and yield the following shrinkage formula , while for a fixed ω, Equation (21) is a quadratic minimization problem and q can be given by Due to the large computation of Equation (23), in our implementation, conjugate gradient (CG) algorithm is adopted to estimate q. During the iterative process, parameters including λ, νand τ should be updated and the detailed updated formulas at l-th iteration can be expressed as [15,16] where ς is a constant. In general, the reconstruction of the image q is implemented iteratively. The initial image is estimated roughly, then the estimated one is used to update ω and D, which are in turn used to improve the reconstruction of q. Such an iterative minimization process terminates until the stopping criterion is satisfied. Specifically, the overall algorithm of the proposed ASR-CS-SAR approach is outlined in Table 1. (a) Compute an initial image q (0) by solving Equation (10)

Computational Complexity
The computational complexity of the proposed approach mainly stems from three parts: the initial reconstruction, the estimation of the weighted AR parameters, and SAR imaging with the adaptive sparse representation. First, a BP algorithm is applied to solve Equation (10) for an initial image, which requires O(NlogN) operations. Next, the computational complexity in the second part includes the calculation of SSP among pixels and the estimation of α and β. The former is solved by Equation (12), which needs O(NU)operations, where U is the order of the AR model, while the latter is obtained by solving a least square problem in Equation (14) MN 2 + MN)operations. In our manuscript, a conjugate gradient (CG) algorithm is adopted to compute Equation (23) to reduce the computer complexity. In fact, matrix-vector multiplication is the dominating operation in the CG algorithm, which can be executed in parallel to speed up the algorithm by taking advantage of the multi-core CPU and GPU structure. In addition, the multiplication A T A can be implemented by fast FFT.As [17] mentioned, if the CG solver converges in k iterations, then the associated computation complexity can be reduced to O(4Nk). Thus the cost of SAR imaging with adaptive sparse representation is O (Nw(4Nk + N 2 )), where Nw denotes the whole iteration number.

Experiments
In this section, experiments over various scenarios are carried out to demonstrate the effectiveness of the proposed approach. In Section 5.1, a visual comparison is conducted between synthetic SAR data. In Section 5.2, quantitative comparisons between our approach and the conventional one using fixed sparse spaces are conducted. In Section 5.3, the application of the proposed approach in real SAR is given.
In our approach, the window size for the estimation of the weighted AR model is empirically set as 11 × 11. During these experiments, the proposed approach is compared with conventional CS-SAR imaging approaches with sparse representation in wavelet [9], DCT [10], and TV [4]. Furthermore, two image quality metrics, including the peak signal to noise ratio (PSNR) and structural similarity (SSIM) [18], are introduced for image quality evaluation, defined as follows: (a) Peak signal-to-noise ratio (PSNR) stands for the ratio between the maximum possible power of a signal and the power of corrupting noise that affects the fidelity of its representation: where q and q denote the original image and the reconstructed image, respectively. MN is the total number of pixels inq. ρmax is the maximum energy of the noisyimage. (b) Structural Similarity (SSIM) is a metric for measuring the similarity between two images, different from PSNR, which is dependent on the average luminance and contrast only. The structural information in an image is defined as those attributes representing an object's structures within the scene; this metric has been proven to be inconsistent with human eye perception. Mathematically, SSIM can be expressed as where μ q and μq represent the mean of the image q and q , respectively. 2 σ q and 2 σ q stand for the corresponding variance. σqq is the covariance between q and q . V1 and V2 are two variables to stabilize the division with weak denominator. The range of values for SSIM lies between 1 and −1. Note that a higher SSIM value corresponds to a better restoration, and vice versa.

Visual Comparison
In this experiment, a visual comparison is presented to show the feasibility of the proposed ASR-CS-SAR imaging approach. The primary radar parameters are enumerated in Table 2. By applying the random undersampling scheme mentioned in Section 2, 33% of Nyquist rate samples in range as well as 25% in azimuth are adopted, which implies that only 1/12 full aperture data is used for reconstruction. Figure 3a,b depict the full echo data and the associated undersampled return, respectively. To facilitate comparison, SAR imaging with full aperture data is displayed in Figure 4.     the image reconstructed by the proposed ASR-CS-SAR approach. From these figures, we can see that most information on the SAR images is preserved. However, for the conventional CS-SAR approaches, they still suffer from a series of artifacts. Specifically, the textures of the roads in the red rectangle window are seriously merged and it is hard distinguish them from the terraced fields. Similarly, the terrace in the pink rectangle window has been blurred. However, those artifacts are suppressed in Figure 5d and a sharper and clearer visual result can be obtained by the proposed approach. The superior performance of the ASR-CS-SAR is also demonstrated by the values of image quality metrics, i.e., PSNR and SSIM, for the maximum gap in PSNR between our approach and the competing methods is up to 7.15 dB.
To further test the stability of the proposed algorithm, we change the sampling rates in azimuth or range with another fixed, and compare the images reconstructed by our proposed algorithm and other algorithms. Figure 6 presents the comparison results. Figure 6a shows that all the PSNR curves rise sharply as the sampling rate in azimuth increases. However, our proposed approach outperforms others for each sampling rate; furthermore, the advantage becomes more and more apparent with the increase of the sampling rate. Similar conclusions can also be drawn from the results in Figure 6b.

Quantitative Comparison and Universality Analysis
To demonstrate the universality of the proposed approach, simulations over a variety of SAR images under different sampling rates (including a low sampling mode and a high sampling mode) were carried out. More specifically, in the low sampling mode, 33% of Nyquist rate samples in range and 25% in azimuth were adopted, while in the high sampling mode, the sampling numbers in range and azimuth increased to 38% and 30%, respectively. All the other radar parameters were the same as in the first experiment. In Figure 7, we list all the SAR images.  Table 3 illustrates the quantitative performance metrics, including PSNR and SSIM, vs. the number of samples for different CS SAR imaging approaches.
It turns out that, for all test images, the proposed approach significantly outperforms other competing approaches in PSNR and SSIM. Especially for the images with rich texture and sharp edges (e.g., mountains and terraced fields), the average gain of proposed approach over other competing approaches can be up to 4 dB in PSNR and 0.02 in SSIM for the low sampling mode. When the sampling number increases, the superiority is enlarged for most complex images. The superiority is mainly due to the capability of the adaptive sparse representation in exploiting the inherent structure of the image. In conclusion, all these quantitative results manifest the fact that the proposed approach is highly competitive and applicable to various types of SAR images.

Application to Real SAR Data
In this experiment, the proposed ASR-CS-SAR imaging approach is applied to real SAR data. The data were acquired by an airborne SAR, and the associated radar parameters are listed in Table 4. In practice, due to the slight deviation of the platform from the nominal trajectory, the inaccuracy of the observation model caused by the additional uncertain factor results in inevitable phase error in the SAR data. Hence, the conventional reconstruction of SAR images suffers from degradation due to the phase error. As shown in Figure 8a, the defocused SAR image has been blurred and is hard to distinguish. To achieve a well-focused SAR image, the unknown phase error should be compensated for first. In our experiment, we adopt the strategy in [19], where the phase error correction and the compressive SAR imaging are iteratively operated. In this strategy, the imaging quality directly influences the estimation accuracy of the phase error, which in turn affects the image reconstruction. To be fair, the same autofocus scheme is carried out in our proposed approach. In this experiment, SAR echoes are undersampled in the azimuth dimension, and only 30% of the original raw data in azimuth is available. For better comparison, the perfectly focused image recovered from full echo data is shown in Figure 8b. The images reconstructed by the CS-SAR imaging approach with motion compensation [19] and our proposed approach are shown in Figure 8c,d, respectively. Compared with the focused image in Figure 8b, most of the image contents have been preserved in Figure 8c, while some regions are still defocused. Specifically, the road marked in the red rectangle has nearly disappeared. In addition, the structure of the target as well as the edge of the lake (in the pink rectangles) is also blurred with the surrounding region. In contrast, all these problems are solved by our proposed approach, and a clearer and sharper image is obtained in Figure 8d. This is due to the fact that a more precise image can be reconstructed by the proposed approach, and thus the unknown phase error can be removed more completely.

Conclusions
By integrating similarity among pixels into the AR model, we propose a weighted AR model to exploit the structural sparsity in a local non-stationary window. A novel compressive SAR imaging approach with adaptive sparse representation was developed. Experimental results over synthetic and real SAR echo data demonstrated the effectiveness and universality of the proposed approach.