Fast Compressed Sensing SAR Imaging based on Approximated Observation
Abstract
In recent years, compressed sensing (CS) has been applied in the field of synthetic aperture radar (SAR) imaging and shows great potential. The existing models are, however, based on application of the sensing matrix acquired by the exact observation functions. As a result, the corresponding reconstruction algorithms are much more time consuming than traditional matched filter (MF) based focusing methods, especially in high resolution and wide swath systems. In this paper, we formulate a new CSSAR imaging model based on the use of the approximated SAR observation deducted from the inverse of focusing procedures. We incorporate CS and MF within an sparse regularization framework that is then solved by a fast iterative thresholding algorithm. The proposed model forms a new CSSAR imaging method that can be applied to highquality and highresolution imaging under subNyquist rate sampling, while saving the computational cost substantially both in time and memory. Simulations and real SAR data applications support that the proposed method can perform SAR imaging effectively and efficiently under Nyquist rate, especially for large scale applications.
I Introduction
Synthetic aperture radar (SAR) is an active microwave radar which can achieve highresolution images in all time of day and weather [1]. In a SAR system, the radar emits a sequence of pulses along its path and receives the echoes (raw data) scattered from the targets. The reconstruction of the scene is traditionally achieved by matched filter (MF) based focusing algorithms which are efficient but need Nyquist rate samples of the echoes. The SAR imaging with increasing resolution and swath requires more and more measurements, storage and downlink bandwidth. The current system hardware, however, frequently hampers such highdimensional application.
The recent development of compressed sensing (CS) brings possibility of reconstructing sparse or compressible signals with fewer measurements than that Nyquist requires[2, 3, 4]. Several applications on radar system appear in recent years, which primarily concern how the data acquisition way can be simplified by using CS[5][6] and what the potential applications will renovate radar imaging with CS technique[7][8]. Further, in the study of CSSAR, much attention has been paid to the effective use of the specific SAR geography and signal form, say, in [9], a SAR raw data compression framework based on CS was suggested by sampling the data in frequency domain. An extension of this work was given in [10] by using the fact that very bright objects are always sparse, resulting in a hybrid sparse model. These works, however, do not apply to the CSSAR system practically where sampling is expected in timedomain. In [11], CS was applied on azimuth after the range compression. By combining range MF, the method was much more efficient, while, the redundant information in range has not been effectively utilized. More general CSSAR model were reported in [12][13] by discretizing the SAR observation function exactly into an observation matrix, while solving by CS straightforwardly.
All those works strongly demonstrated that some exclusive advantages of CSSAR do exist as compared to the traditional SAR imaging methodologies, say, relaxation of required measurements, reduction of side lobe and a further suppression of noise [14]. However, in all applications, a serious drawback has been observed: as compared to the traditional MF based methods, the computational complexity and memory cost of the CSSAR models are much higher, so that it is very inefficient to be applied to highdimensional applications.
In this paper, we formulate a new CSSAR framework within which the computational complexity of the CSSAR imaging can be significantly reduced. Our main idea is to replace the exact observation function in the CSSAR framework with approximated observations derived from the inverse of traditional MF based procedures. Such inversion has ever been applied to yield raw signals (the echoes) in a more economical way [15][16], but requires high accuracy of the adopted method. In this paper, we take a further step by incorporating it into the CS framework, which demands only a well focusing ability to ensure CS reconstruction. We propose to implement the CSSAR imaging through the sparse regularization scheme which is then solved by an iterative thresholding algorithm (ITA). Accordingly, the fast speed and high efficiency of the new method are guaranteed respectively from the use of the approximated observation and from thes CS reconstruction procedure. We show that the new CSSAR imaging method can not only acquire highquality and highresolution images with significantly reduced measurements, but also reduces the memory cost to and computational complexity of onestep iteration to , achieving the same order with the traditional SAR imaging methods.
The reminder of the paper is organized as follows. In Section 2, we introduce the background knowledge on the stripmap mode SAR system and the classical CSSAR model. In section 3, we present the approximated observation by calculating the inverse of MF imaging procedure. In Section 4, we formulate the new CSSAR imaging method through hybridizing the approximated observation and sparse regularization. In Section 5, we show the simulation and application results of the suggested method. Conclusions are then presented in Section 6 with some useful remarks.
Notation: We will use the subsequent notations throughout the paper: Column vectors, matrices and operators will be denoted respectively by bold lower case, , bold upper case, A, and roman upper case, . denotes the transpose, conjugate and Hermitian transpose of , respectively.
Ii CSSAR models based on exact observation
In this section, some preliminary knowledge of CSSAR imaging is summarized. We focus on the general formalization of CSSAR model, with a more detailed introduction of the iterative thresholding procedures for solution of the CSSAR models.
Iia Stripmap Mode SAR Model
In the stripmap mode SAR, the antenna is pointed to a fixed direction and the platform flights with constant velocity . Then, a complex baseband , usually chirp, is modulated to real pulse ( is the carrier frequency, is the range time, is the elevation weight and is the pulse duration) and transmitted at a constant pulse repetition frequency (PRF). The received backscattered energy can then be modeled as a convolution of the pulse waveform with the ground reflectivity function, given by [17]
(1) 
where are respectively the azimuth and range time, denotes the additive noise, is the timevariant convolution kernel which can be composed as:
(2) 
In (2), is the twodimensional azimuth modulation which is responsible for the alongtrack observation while is range convolution kernel that is identical to the transmitted pulses.
Further, we can sample the continuoustime analog echo and discrete the reflectivity map , into twodimensional arrays and . And then we obtain the following observation model for the strip mode SAR:
(3) 
where , , is the observation matrix acquired from the discrete weight of (1) (more detailed information and construction of the observation model can be seen in [17][18]), and is the noise.
IiB Formulation of CSSAR models
In a CSSAR model, the data is sampled and compressed with a proper sampling matrix , resulting in
(4) 
When is a sparse signal, say, most of the entries of are zeros, the theory of CS tells when and how it can be recovered from the above undetermined linear system with fewer measurements than Nyquist criterion requires[2][4]. Generally, considering an illposed linear system () where is sparse enough, if the sensing matrix satisfies some conditions like RIP[19], can be exactly recovered from with the (quasinorm) () optimization:
(5) 
To solve (5), we usually use an equivalent regularization scheme with the following optimization problem
(6) 
where is a regularization parameter. The optimization can be efficiently solved by iterative thresholding algorithm (ITA)[20, 21, 22]. In detail, an ITA generates a sequence of approximates according to:
(7) 
where is a normalized parameter which controlls the convergence of the iteration. In (7), () is a socalled thresholding operator which is componentwisely defined by
(8) 
where can be analytically specified when . For example, the widely used softthresholding, which corresponds to , is
(9) 
The iteration (7) with (8) is the fundamental procedure we suggest to use for the CSSAR imaging.
It can be seen that the main computation load in implementation of (7) comes from the calculation of time domain correlation . From the viewpoint of SAR signal processing, this corresponds to the backscattered projection procedure, which is known to be inefficiency for reconstruction, even implemented by convolution as in (1). On the other hand, we notice that there exists efficient focusing methods using MF in traditional SAR signal processing. This type of processing is in the frequency domain, which is much faster. Moreover, unlike the MF based method which is usually decoupled, the sensing matrix in (7), owns intrinsically twodimensional structure that has to be collected and stored before imaging. Although some compression can be incorporated according to the structure of the matrix, it still consumes a huge memory. All these difficulties then hamper effective applications of the known CSSAR imaging.
The aim of the present research is to suggest a new CSSAR imaging method, which replaces the use of the exact observation model by an welldefined approximation, and then makes it possible to reconstruct the sparse scene via a sequence of 1D operations. Thus, the very high cost of calculation and memory of the existing CSSAR imaging methods can be significantly reduced.
Iii Approximated Observation
In this section, we first explain why an approximated observation is needed and feasible, and then, we provide an example to show how an approximated observation operator can be explicitly constructed by virtue of a concrete example from the inverse of RangeDoppler Algorithm (RDA). A relation between the constructed approximation and the corresponding focus method is analyzed, which then serves as the basis of the development of new method in the next section.
Iiia Why approximation needed
It is known that MF is fast with complexity, which is, among the others, mainly due to frequency domain operations. More precisely, if we denote the imaging procedure by MF, like RDA, the SAR raw data can be well focused in some conditions by
(10) 
where is the traditional MF imaging procedure that can be calculated through decoupling it into a series of 1D operators in the frequency domain. This normally leads to an complexity when fast fourier transformation (FFT) type operations are employed.
Observing these advantages, the purpose of this paper is to accelerate the known CSSAR imaging procedures, so as to achieve a comparable (at the same order) complexity with the traditional MF based methods.
A natural consideration is then to look for the possibility of integrating CS and MF. However, a direct application on the decoupling of is impossible, because intrinsically possesses 2D structure. Nevertheless, it is known from (10) that always approximates , say, , and hence, , whenever exists, approximates the observation . In viewing that is decoupled, it can then naturally be expected that is decoupled, too. Thus, we can expect that under certain conditions, some types of approximations of can be decoupled, so as to bring an complexity.
This is why we would like to approximate the observation, and in the following, we will introduce the details on how to construct and what constraints an appropriate approximation observation.
IiiB How to construct an admissible approximated observation
Fig. 1 draws the main relations between CSSAR observation and MF reconstructions. It can be seen that whenever the imaging procedure is accurate enough, can be viewed as an admissible substitute of . This provides a general principle of how the observation can be remodeled and approximated by any high precision imaging (or reconstruction) procedure. We formalize this principle further as
(11) 
where is any generalized right inverse of , and is any a high precision imaging procedure. We call henceforth an approximated observation.
However, since there are many well known imaging procedures that provided various tradeoffs on imaging accuracy and complexity. We need therefore further to define the extent of accuracy and identify the constraints under the CSSAR framework. To see this, let us compare the exact observation model and the approximated observation model
(12) 
It can be seen that by using approximated observation, other than reconstruct , we actually reconstruct instead, which is assumed to be an approximation of , when it obeys to the following relation
(13) 
where denotes the Hadamard product, denote the phase error while is the error for side lobe, or more severely, the artifacts from unfocusing. Formally, when (13) holds, there exists an acceptable solution with the approximated observation model. However, to find it under CS, we should further emphasis on a better focusing ability of .
As we know, a key parameter in CSSAR, different from traditional SAR, is the sampling rate which measures how a SAR system benefits from CS. The least amount of samples to ensure the reconstruction is incoherently determined by the sparsity of scene , usually irrelevant of the distribution and phases of targets in the scene. Thus, the difference of the sparsity between and determines when and how much the additional measurements does the approximated observation based CSSAR methods need. It can be immediately seen from (13) that the difference is uniquely characterized via . That is to say, whenever the side lobes reconstructed from is low enough, can be ignored, then and can keep the same sparsity. In this situation, the required least sampling rate of the approximated observation model equals to the original model. In turn, to prevent the approximated observation based CSSAR method from demanding more samples, the side lobe should be as low as possible, or equivalently to say, a well focusing capacity should be a criterion to determine whether a specific focusing method can be used to construct the approximated observation.
The above discussions tell a fact that the construction of the approximated observation is quite reflexible, which can be acquired straightforwardly based on well established algorithms with additional requirement on the focusing ability. In the next subsection, we present a concrete example using RangeDoppleralgorithm (RDA)[23] to show how an admissible approximated observation can be simply constructed based on this principle.
IiiC A concrete example
RDA is a very popular procedure for stripmap mode SAR imaging that is simple both in comprehension and in implementation. The procedure (under the low squint case) consists of three main steps (operations): the range compression, RCMC and azimuth compression. In a compact form, the imaging procedure , operated on 2D array, can then be expressed in the following
(14) 
where () is the reconstructed 2D SAR image, respectively are the DFT matrix and inverse DFT matrix (in practice, they are implemented by FFT) to perform, the subscript denotes the direction of azimuth and range where the FFT performs along, and are the frequency domain matched filter operations along azimuth and range, which can be always defined respectively by
(15) 
In (15), are the frequency along Doppler and range, and are the azimuth FM rate and the pulse FM rate. In (14), is the RCMC interpolation operator which is essentially a spacevariant shift, and always approximated by the truncated sinckernel interpolation with as
(16) 
where is the migration (measurement in time) to be corrected, and are the signals before and after RCMC, respectively.
With the so specific operations in RDA procedure, we now can derive the inverse of quite simply by taking the inverse of every subprocedures. The details are as follows:
i) The inverse of Fourier transformations are known as the inverse transformations, which are given by . It is important to keep the throwaway consistent between the pairs.
ii) It is known that phase multiplication is a unitary transformation, so that the inverse is the multiplication of the conjugate phase, , and the Hadamard multiplication can still be applied in order.
iii) The inverse of is difficult to achieve directly. In fact, is approximated from the accurate RCMC defined in continuous range time domain. Because the trajectories of targets with different range gates are disjoint, this shift is a onetoone mapping and the inverse of the origin RCMC exists. We can also approximate through interpolation that
(17) 
Based on the above exposition, the approximated observation deduced from RDA can then be explicitly expressed by:
(18) 
We show that the so constructed approximated observation has an interesting property: It is still a linear operator and its conjugate transposition equals to ^{1}^{1}1It is to say is nearly unitary, since only a minor approximation on calculation of the inverse is included.
Theorem 1. is a linear operator with the property .
Proof:
The linearity of and is obvious because all the suboperations are linear. Let denote the vector form of , namely, . Then, by definition, the linear operators and can be written as matrices, and we then have
(19) 
(20) 
where
(21) 
(22) 
and are real matrices defined by
(23) 
In (23), , is the discretion of . Observing from (23) and is the pulse sampling interval, it is easy to check that . Consequently, comparing (20) and (19), we conclude that . \qed
Theorem 1 shows that we have actually taken the conjugate transposition of as an approximated observation of . Such coincidence plays an important role in the new method to be suggested in the next section.
IiiD Generalization
The approach we have applied to define the approximated observation by the inverse of RDA procedure can be generalized in the following two folds:
1) For high squint cases, we can incorporate secondary compression in RDA or derive an approximated observation from the inverse of other focusing methods like ChirpScaling Algorithm (CSA), algorithms, so as to enhance the focusing ability. With similar suboperations, like FFT, phase multiplication and interpolation, the acquisition of the inverse is just as same as RDA.
This principle can be applied to yield more general algorithms, however, we will not enumerate all possible extensions but remind that the decoupled structure of MF makes the inverse always achievable. This is the reason why MF is fast and why we propose to apply the approximated observation instead of the exact observation in the CSSAR imaging system.
2) In the above derivation, we have assumed that the transmitted signal form is standard chirp and focus both the azimuth and range direction. In fact, the azimuth modulation is the exclusive property and main difficult of SAR signal processing while the range convolution, which possesses a simple 1D structure, can be modeled directly. So, we can apply the approximated observation to nonchirp cases, by replacing with the transmitted pulse, which is always recorded in modern SAR systems. This extension is also of specifical necessity in CSSAR because the design of the pulse form is also a very important issue.
Iv CSSAR imaging based on approximated observation
In this section, we formulate the new CSSAR imaging method based on the use of approximated observation. An regularization model together with the fast iterative thresholding algorithm will be suggested.
Iva The new CSSAR method
By replacing the exact observation using the approximated observation in (18), we can acquire the following CSSAR model:
(24) 
where is the Frobenius norm of a matrix, and are the sampling operators in azimuth and range directions, which corresponds to the general sampling operator in (4). It is well defined because the azimuth signal is of discrete form and the range signal is of continuous form, and thus the sampling procedure of the two types of signals are usually physically separated^{2}^{2}2Although the sampling scheme on range may vary pulsebypulse, we still use this expression which is easily understood.
Then, due to the linearity of , the model (24) can still be very fast solved by ITA, which reads in this case that
(25) 
In this paper, we simply select while parameters will be preset according to the next subsection.
Fig. 2 below shows the flow chart of algorithm (25), which tells that ITA provides an intuitive explanation in terms of SAR signal processing. It is seen that at each iteration, the ITA can be decomposed into mainly three procedures: the compressed data simulation, the matched filter on the residual and the thresholding for new estimation. Physically, this means that in every iteration of the ITA, the useful information in the residue (not the raw data) is first extracted by MF and then added to the current estimate to yield a new update, finally the thresholding procedure enforces the sparsity through regularizing the noise and ambiguity from undersampling.
IvB Parameter Setting
There are two parameters and in (25) that need to be set. First, controls convergence of the ITA that the inverse should obeys.
(26) 
However, it is difficult to calculate directly, where operator is included. As an alternative, we adopt the adaptive step selection strategy in [24] as:
(27) 
where equals to at the support of , and equal to zero elsewhere. It is easy to demonstrate that satisfies (27), and as reported in [24], such choice has an additional advantage of accelerating the algorithm.
Further, the regularization parameter , which functions to compromise the reconstruction precision and the sparsity of the solutions obtained, has a substantial impact on the imaging result. Fortunately, as a part of the regularization theory, the optimally has been resolved in [21], whenever the problem’s sparsity is known. More precisely, assume the considered problem has sparsity (i.e. a sparsity problem), then the optimal setting problem of parameter is shown to satisfy
(28) 
where , is its th largest component in magnitude.
Therefore, we suggest the setting that in the th iteration ( is independent of ). The sparsity , which determines , can be much flexible to be set, say, based on a prior upper estimation on sparsity of the target scene.
IvC Computation Cost
Let us compare the computational complexity and the memory occupation of the suggested CSSAR imaging method (24), as compared with the known CSSAR model (6). The purpose is to see how much reduction of computational cost of the new model has been brought. In the calculation, we have used some standard notations which are: the number of required iteration , the sampling rate , the number of range gates , the number of range lines (), the number of samples of sent pulse , the number of samples of the synthetic aperture time (they equal to the time bandwidth product (TBP) in each direction), and the TBP of radar signal .
With these notations, we can calculate the computational complexity of the approximate observation based CSSAR, , and the computational complexity of the exact observation based CSSAR, as follows. For , it includes calculations of an inverse MF procedure and a MF procedure, which has commonly the computational complexity of , together with a decoupled thresholding operator with complexity in a single step. Thus, the total cost is at the order . For , it includes calculations of a single iteration, two matrix multiplications and the thresholding procedure. Since there are only few nonzero entries in , say, nearly in every column, when coding it using twodimensional convolution, it needs at least complex multiplication. Thus, we find that the total cost is . Then, the ratio between and is given by :
(29) 
It is seen from (29) that the ratio depends linearly on the TBP of radar signal . In SAR applications, the is always designed very large (thousands even millions) to improve the reconstruction signal to noise ratio (SNR), which will bring very high computational cost of the time domain reconstruction method.
The memory loads of the approximate observation based CSSAR, , and the memory loads of the exact observation based CSSAR, can be estimated in the subsequent way. For , it contains only the storage of input, output and several parameter matrices (i.e., azimuth matched filter, range matched filter and the amount of migration in RDA), which is summed up to bytes memory occupation. For , although no filters are stored, it needs additionally to store a sensing matrix, with the number of nonzero entries of . But, because the Doppler history with same range cell share the identical patterns, we only need to store an intact holistic pattern (the convolution kernel) for each range gate to achieve a compression, resulting in an additional memory occupation of bytes (a complex number occupies 16 bytes memory), as compared with . This additional cost can be very large in spaceborne SAR systems. For example, when and , it requires more than 100 GB memory to store the array. But the memory cost of RDA is only a few hundreds MB in the same condition. This will further hamper the application of time domain methods into practice.
Finally, the required number of iteration steps is difficult to compare analytically, but in practice, no obvious difference is observed.
IvD A Summary
From the analysis in the previous subsections, we can see that the suggested new model (24) and method (25) have constituted a more efficient CS based SAR imaging method. While preserving CS features, the new method has the following exclusive advantages:
Lower computational cost: Due to the use of approximated observation, the method only involves 1D operations which makes the imaging process extremely efficient. It has reduced the computational complexity of the existing exact observation based method significantly, as shown in (29). Meanwhile, taking full advantages of the decoupled structure after approximating the observation, the proposed method can save the memory cost with a remarkable amount, which is sometimes of more significances.
Higher Compatibility: As compared with the traditional MF based SAR imaging procedure, the new method uses the same or similar operations, except an additional thresholding operation that yields the sparsity of solution. In particular, the new procedure can be seen as a successive iterative refinement of the well known MF based method, which makes the new method more consistent. As a result, the proposed model requires little modification of the existing SAR imaging algorithms, which makes the combination of MF and CS much simpler.
All these features make the suggested new CSSAR imaging method more useful and efficient, and particularly possible to be applied in high dimensional SAR applications. We will provide simulations and applications in the next section to further support such benefit.
V Simulations and applications
In this section, several simulations and applications are provided to demonstrate the effectiveness and efficiency of the proposed CSSAR imaging method. For abbreviations, we denote by CSRDA the CSSAR imaging method (24) with the approximated observation acquired from the inverse of RDA, and by CSEO the CSSAR method (6) with the exact observation .
We first conduct a series of simulations to compare the performance of the CSRDA method, the CSEO method and the traditional RDA method in terms of reconstruction ability, reconstruction quality and reconstruction cost. Then, we apply the CSRDA to some real SAR imaging tasks from RADARSAT1, which then further demonstrates the outperformance of the suggested method.
The sampling scheme used in the simulations are specified as follows. In the azimuth direction, we employed random downsampling, realized by selecting random rows from the raw data , with sampling rate . In the range direction, we picked up random samples independently on each sampled echoes in azimuth, with sampling rate . And we keep the ratio between and as ^{3}^{3}3The suggested sampling strategy was designed to comprehensively compare the reconstruction algorithms. The proposed model itself is adaptively to more complicated sampling schemes, for example, jitter sampling in the azimuth direction [12] and random demodulation[25] in the range direction. .
Table 2 lists the primary SAR parameters used in both simulations and applications. All the experiments were conducted on a work station of 8core 2.4GHz CPU with 32G memory. The CSRDA was implemented in MATLAB 2012a while the CSEO using optimized convolution was implemented in C++ with parallel codes and careful array operations.
Parameter  Symbol  Simulation  RadarSat1 
Slant range of scene center()  20  1016.7  
Effective radar velocity()  350  7062  
Beam squint angle()  0  0.06  
Radar center frequency()  5000  5300  
Pulse repetition frequency()  175  1256.98  
Range FM rate()  37.5  0.72135  
Pulse duration()  2  41.75  
Sampling rate()  75  32.317 
Va Simulations
In the simulations, the scene was taken as while the scattered coefficients were chosen with unit amplitude and uniform random phases. The raw data were first generated in timedomain by exact slant range and then sampled with different rate, to yield the compressed measurements. The sparsity parameter was kept the same for CSRDA and CSEO and the maximum iteration steps was set to 100 for both methods. The aims of simulations is then to compare the reconstruction ability (RA), reconstruction quality (RQ) and reconstruction cost (RC), of each competitive SAR imaging methods. These are measured respectively by the lowest amount of measurements a method can successfully reconstruct an image, the side lobe and resolution of reconstructed point target, and the computation time cost by a method to recover the image.












1) RA comparison: A set of simulations was made where 9 targets were located at the center with intervals of 6 samples. We varied the sampling rate ranged from 100% to 0.6% and added gaussian noise with level of 20 dB. We applied RDA, CSRDA and CSEO to this experiment with sparsity parameter . Some of the simulation results are shown in Fig. 3 and 4.
It is seen from the top row of Fig. 3 that with full samples (namely, with 100% sampling rate), all the methods RDA, CSEO and CSRDA can successfully recover the scene, say, the amplitude of the target is maintained and no false target is observed. But the reconstruction of RDA is with serious side lobes, which is not observed in CSRDA and CSEO. This shows the exclusive advantage of the sparse regularization based CSSAR imaging methodologies, as reported in [14]. When we reduce the sampling rate, say, 10% samples, as seen in the bottom of Fig. 3 RDA fails to recover the scene, while CSEO and CSRDA both can not only perfectly reconstruct the scene but with significantly reduced side lobes. In this case, no visible difference can be observed for CSEO and CSRDA. Nevertheless, when the sampling rate continues reducing as in Fig. 4, we found that both CSRDA and CSEO can reconstruct the image with only 0.65% of the samples. However, CSRDA fails with 0.55% samples while CSEO fails until the sampling rate takes 0.45%.
All the results consistently show that benefited from sparse regularization, the approximate observation based CSSAR method can reconstruct sparse scenes with far less samples than Nyquist rate requires. However, caused by approximation, it requires a little more samples to reconstruct the scene.
2) RQ comparison: Sparse regularization was demonstrated in SAR and CSSAR imaging the ability of reducing the side lobe and simultaneously improving the resolving ability[14]. Hence, we are interested in whether the enhancement is kept when approximated observation is included, especially when the effect of the accuracy of the observation is excluded. To illustrate it, we compare the reconstruction quality of RDA and CSRDA in terms of side lobe and spatial resolution, when successful reconstruction is achieved. The side lobe is evaluated via the peak side lobe ratio (PSLR), defined as the ratio of the peak intensity of the most prominent side lobe to the peak intensity of the main lobe, say, the smaller the PSLR, the better an algorithm. The spatial resolution is measured by the impulse response width (IRW), defined as the width of the main lobe of the impulse response, measured by 3 dB below the peak value, or the minimum distance an algorithm can separate two targets, which should be also as smaller as better. We have perform a one point simulation with the upsampled factor 16 while the target is analyzed by a chip centered on the peak, to yield a more detailed analysis on both the main lobe and the side lobe. The sampling rate is fixed (10% for CSRDA while 100% for RDA) to ensure a successful reconstruction, and the sparsity parameter is set 600. Some of the simulation results are given in Fig. 5.
Fig. 5 show the contours of the reconstruction results, with contour lines of 3 dB (the boundary of main lobe) and 13 dB (the PSLR of traditional MF output). The comparison intuitively shows that both the area of the main lobe and the PSLR in the side lobe reconstructed from CSRDA is much smaller than those reconstructed from RDA. Further in the table, details on the reconstruction quality are presented in azimuth and range directions. It is seen that the width of main lobe reconstructed from RDA is 15 samples both in range and azimuth, but for CSRDA, the width is only 8 samples in the two directions. For the side lobe, the PSLR reconstructed from RDA is 13.32 dB in range and azimuth direction, which is very obvious. But the CSRDA reduces the PSLR to 21.3 dB and 22.7 dB in azimuth and range, respectively.
We further demonstrate the enhancement of resolution by CSRDA in another way. We add two point targets in the above simulation, with respectively azimuth and range interval of 12 samples from the center. We then test whether the three targets can be separated from CSRDA. It is seen from the simulation result (Fig. 6(a)) that the reconstructed result of traditional RDA exhibits an overlapping of the main lobe, thus the targets are inseparable. By using CSRDA in Fig. 6(b), one can clearly distinguish the locations of the three points, demonstrating an obvious enhancement of the resolution.
All of the simulations in this subsection demonstrate that by combining MF with sparse regularization, we can reduce the side lobe and improve the resolution simultaneously to a great extent. Note that, these two goals have been regarded as tradeoffs traditionally, if only MF is employed.
3) RC comparison: Finally, we compare the CPU time takes by CSRDA and CSEO. According to the analysis in section IV, the computational complexity of CSRDA and CSEO depends on the scene size and TBP of radar signal . So we generated 10 examples for each set of fixed scene size and TBP of radar signal in simulation, with a constant sampling rate 10%. The average computational time in a single iteration of the two methods was then recorded. The comparison results are shown in Fig. 7.
As we can observe from Fig. 7, when is fixed as , the CSRDA scales very well to very high dimensional problems, since even with , it only takes several second to finish an iteration. The CSEO is also insensitive to dimension when is fixed. However, the CPU time of CSEO is consistently higher than that of CSRDA, with a ratio around 100. On the other hand, it is seen from Fig. 7(a) that when is fixed to , the CPU time of CSRDA is constantly 3s, but, the CPU time of CSEO depends linear on , which becomes more and more costly as increases.
This RC comparison shows that the approximate observation based CSRDA is much faster than the time domain method CSEO, benefited from the computational cost of MF. And the computational complexity of CSEO depends linearly on both the TBP of radar pulse and the scene size. So when the is large, i.e. which is commonly in spaceborne cases, CSRDA is expected to accelerate the CSSAR reconstruction more than thousands of times. This acceleration of computational time together with the shown memory saving capacity demonstrate the superiority of the proposed method.
All the above simulations support that the suggested approximated observation based CSSAR method, CSRDA, is capable of high quality imaging under Nyquist rates and with a comparable complexity as the traditional MF based imaging method. Especially, as compared with the exact observation based CSSAR imaging, CSRDA preserves the features of imaging under Nyquist rate and reconstruction with feature enhancement, while reducing the imaging complexity dramatically. Such significant complexity reduction property makes the new method applicable to large scale imagery application (as will be demonstrated in the next subsection). This is, however, at the sacrifice of reconstruction quality or, equivalently, must be compensated with additional measurements. Thus, the new method provides a satisfying tradeoff between the reconstruction complexity and quality.
VB Application
We have applied the new method, CSRDA, along with RDA and CSEO, to some real SAR imaging tasks. RADARSAT1 is a famous satellite SAR launched at 1995, and the data used in this application were collected on June 16, 2002 with Fine Beam 2 about Vancouver region. The related key parameters of SAR system are as in Table. 1.
We first applied the 3 methods to reconstruction of the region of English Bay, in which 6 vessels are sparsely distributed, a very typical sparse scene. The scene was digitalized as image with azimuth resolution 9 m and range resolution 6 m. We then reconstruct the image by the 3 imaging methods with sparsity , and with varied sampling rates from 100% to 5%. Some typical results of reconstructions are shown in Fig. 8. As expected, the application shows a completely similar performance as that in the simulations. For example, Fig. 8(a) exhibits that RDA can only reconstruct the image when samples are fully adopted, but strong side lobe is observed. When sampling rate is 20%, RDA fails with obvious ambiguities, but CSRDA and CSEO both can perfectly recover the image, with much reduced side lobe. Fig. 8(c)(f) then show that when sampling rate is reduced to 5%, the reconstruction of CSEO is with slightly higher precision than that of CSRDA, though both can still recover the targets. On the other hand, as listed in Fig. 7, CSRDA exhibits its dominant advantage in computational cost as compared with CSEO. For example, the computation time of reconstruction, when 20% samples are used, by CSRDA is 1 minutes, while by CSEO is about 9 hours. Also, CSEO needs to store the sensing matrix which occupies about 16Gb memory, while only 100Mb for CSRDA.
We further applied the CSRDA to the large scale imaging problem together with RDA. However, the CSEO is not compared in this experiment because the memory cost is beyond the computational ability of our computer. The scene has a size samples, which is of large scale but not so sparse. CSRDA can be applied in principle because sparse regularization, which is adopted in CSSAR, can be used as a feature enhanced imaging method (with suppressed side lobe and improved resolution), as demonstrated in the simulations. The reconstruction results by RDA and CSRDA (with 100% sampling rate) are shown in Fig. 9. It is seen from Fig. 9 that CSRDA has resulted the reconstruction with improved resolution and reduced side lobes, as demonstrated in the zoomed Fig. 10.
The applications above support that the suggested approximated observation based CSSAR imaging is effective and efficient, especially applicable to high dimensional SAR imaging applications. Such benefit clearly improves on the currently used exact observation based CSSAR imaging methods.
Vi Conclusion
Compressed sensing (CS) has been applied to yield novel SAR imaging methodologies under Nyquist sampling in recent years. The resultant CSSAR models are time domain based and using the exact observation, which then makes it of very high computational cost, and is difficult to be applied in high dimensional applications. In this paper, we have proposed an approximated observation and frequency domain based CSSAR imaging method, with which the computational complexity can be dramatically reduced.
The main contributions of the present work are as follows.
i) Instead of the exact observation matrix, an operator, called the approximated observation is constructed to generate SAR raw data by means of inverse of any traditional MF based imaging procedure (like RDA). Such generic construction makes the SAR imaging approximately decoupled, very fast in processing, while naturally connecting the existing MF based SAR imaging algorithms.
ii) Incorporating the approximated observation into CSSAR framework, an efficient sparse regularization based CSSAR model is formulated. The new model combines naturally CS and MF and is compatible with most existing SAR imaging methods, which needs a little modification of the current SAR imaging technologies.
iii) With the use of approximated observation, an iterative thresholding algorithm is suggested for fast solution of the new CSSAR model, which forms a low complexity, CS featured, new SAR imaging method.
We have tested and applied the new suggested CSSAR method with a series of simulations and applications. The experiments consistently support that the new method is capable of reconstructing sparse scenes with far fewer measurements than Nyquist requires, yielding always a feature enhanced high quality imaging and bringing a speed up of reconstruction hundreds times as compared with the exact observation based CSSAR methods. Due to the fast and feature enhancement features, the new CSSAR method can be accepted as an efficient CS type SAR imaging technique, especially for high dimensional imaging applications.
It is worthwhile, however, to remark that although significant complexity reduction (speedup of reconstruction) can be brought, the use of the approximated observation requires more samples to reconstruct a scene. Thus, how and what extent does the approximation affect the reconstruction deserves a further study. Moreover, since there are many possibilities of concrete realizations of the approximated observation, the criterion on how to select an appropriate one deserved study. All those problems are under our current research.
Acknowledgment
This work was supported by the State Key Development Program for Basic Research of China (973 Program) (Grant No.2010CB731905), Key Program of National Natural Science Foundation of China (Grant No. 11131006), and the National Natural Science Foundations of China (Grants No. 61075054, 60975036, 11171272)
References
 [1] I. G. Cumming, F. H. Wong, U. of British Columbia, and M. D. Dettwiler, Digital Signal Processing of Synthetic Aperture Radar Data: Algorithms and Implementation. Artech House, 2004.
 [2] E. J. Candès, “Compressive sampling,” in Proceedings on the International Congress of Mathematicians, 2006, pp. 1433–1452.
 [3] R. G. Baraniuk, “Compressive sensing,” IEEE Signal Processing Magazine,, vol. 24, no. 4, pp. 118–121, 2007.
 [4] D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
 [5] A. C. Gurbuz, J. H. McClellan, and W. R. Scott Jr, “Compressive sensing for subsurface imaging using ground penetrating radar,” Signal Processing, vol. 89, no. 10, pp. 1959–1972, 2009.
 [6] M. A. Herman and T. Strohmer, “Highresolution radar via compressed sensing,” IEEE Transactions on Signal Processing, vol. 57, no. 6, pp. 2275–2284, 2009.
 [7] J. H. G. Ender, “On compressive sensing applied to radar,” Signal Processing, vol. 90, no. 5, pp. 1402–1414, 2010.
 [8] L. C. Potter, E. Ertin, J. T. Parker, and M. Cetin, “Sparsity and compressed sensing in radar imaging,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1006–1020, 2010.
 [9] S. Bhattacharya, T. Blumensath, B. Mulgrew, and M. Davies, “Fast encoding of synthetic aperture radar raw data using compressed sensing,” in IEEE/SP 14th Workshop on Statistical Signal Processing, 2007, pp. 448–452.
 [10] G. Rilling, M. Davies, B. Mulgrew et al., “Compressed sensing based compression of SAR raw data,” 2009.
 [11] M. Tello Alonso, P. LópezDekker, and J. J. Mallorqui, “A novel strategy for radar imaging based on compressive sensing,” IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 12, pp. 4285–4295, 2010.
 [12] V. M. Patel, G. R. Easley, D. M. Healy, and R. Chellappa, “Compressed synthetic aperture radar,” IEEE Journal of Selected Topics in Signal Processing, vol. 4, no. 2, pp. 244–254, 2010.
 [13] J. S. Zeng, J. Fang, and Z. B. Xu, “Sparse SAR based on regularization,” in SCIENCE CHINA Information Sciences, 2012(In press).
 [14] M. Çetin and W. C. Karl, “Featureenhanced synthetic aperture radar image formation based on nonquadratic regularization,” IEEE Transactions on Image Processing, vol. 10, no. 4, pp. 623–631, 2001.
 [15] A. S. Khwaja, L. FerroFamil, and E. Pottier, “SAR raw data simulation using high precision focusing methods,” in EURDA, 2005, pp. 33–36.
 [16] G. Franceschetti, R. Guida, A. Iodice, D. Riccio, and G. Ruello, “Efficient simulation of hybrid stripmap/spotlight SAR raw signals from extended scenes,” IEEE Transactions on Geoscience and Remote Sensing, vol. 42, no. 11, pp. 2385–2396, 2004.
 [17] C. Wu, K. Y. Liu, and M. Jin, “Modeling and a correlation algorithm for spaceborne SAR signals,” IEEE Transactions on Aerospace and Electronic Systems, no. 5, pp. 563–575, 1982.
 [18] S. J. Wei, X. L. Zhang, J. Shi, and G. Xiang, “Sparse reconstruction for SAR imaging based on compressed sensing,” Progress In Electromagnetics Research, vol. 109, pp. 63–81, 2010.
 [19] E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.
 [20] I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on pure and applied mathematics, vol. 57, no. 11, pp. 1413–1457, 2004.
 [21] Z. B. Xu, X. Y. Chang, F. M. Xu, and H. Zhang, “ regularization: a thresholding representation theory and a fast solver,” IEEE Transaction on Neural Networks and Learning Systems, 2012(In press).
 [22] T. Blumensath and M. E. Davies, “Iterative hard thresholding for compressed sensing,” Applied and Computational Harmonic Analysis, vol. 27, no. 3, pp. 265–274, 2009.
 [23] C. Wu, “A digital system to produce imagery from SAR data,” in Systems Design Driven by Sensors, vol. 1, 1976.
 [24] T. Blumensath and M. E. Davies, “Normalized iterative hard thresholding: Guaranteed stability and performance,” IEEE Journal of Selected Topics in Signal Processing, vol. 4, no. 2, pp. 298–309, 2010.
 [25] J. A. Tropp, J. N. Laska, M. F. Duarte, J. K. Romberg, and R. G. Baraniuk, “Beyond Nyquist: Efficient sampling of sparse bandlimited signals,” IEEE Transactions on Information Theory, vol. 56, no. 1, pp. 520–544, 2010.