Sub-Nyquist SAR via Quadrature Compressive Sampling with Independent Measurements

: This paper presents an efﬁcient sampling system for the acquisition of synthetic aperture radar (SAR) data at sub-Nyquist rate. The system adopts a quadrature compressive sampling architecture, which uses modulation, ﬁltering, sampling and digital quadrature demodulation to produce sub-Nyquist or compressive measurements. In the sequential transmit-receive procedure of SAR, the analog echoes are modulated by random binary chipping sequences to inject randomness into the measurement projection, and the chipping sequences are independent from one observation to another. As a result, the system generates a sequence of independent structured measurement matrices, and then the resulting sensing matrix has better restricted isometry property, as proved by theoretical analysis. As a standard recovery problem in compressive sensing, image formation from the sub-Nyquist measurements has signiﬁcantly improved performance, which in turn promotes low sampling/data rate. Moreover, the resulting sensing matrix has structures suitable for fast matrix-vector products, based on which we provide a ﬁrst-order fast image formation algorithm. The performance of the proposed sampling system is assessed by synthetic and real data sets. Simulation results suggest that the proposed system is a valid candidate for sub-Nyquist SAR.


Introduction
Synthetic aperture radar (SAR) is an imaging radar mounted on a moving platform.As the platform moves, electromagnetic waves are sequentially transmitted to illuminate an area and the backscattered echoes are collected by the radar receiver to form raw data for subsequent digital processing.With increasing demands for SAR products, modern SAR systems are facing several problems.On the one hand, the raw data has resulted in a large burden on board and subsequent storage as well as downlink transmission, especially for spaceborne repeat-pass observations.Nowadays, numerous on-orbit earth observation missions are collecting SAR data continuously and have accumulated huge historical archives.On the other hand, it is often required to use large bandwidth signals for high-resolution applications, which thus need high-speed analog-to-digital converters (ADCs).These requirements present a big challenge for SAR development.
Compressive sensing (CS) theory [1], which emerged in the past decade, brings us new concepts for sub-Nyquist acquisition of SAR data.For a sparse or compressible radar scene, it is sufficient to recover a high-resolution SAR image with a small number of samples.Past investigations find that CS offers various benefits for SAR, such as reduced storage cost, increased resolvability of point targets [2], and reduced speckle [3].In sub-Nyquist SARs, CS can be applied both in azimuth and range dimensions.Azimuth-dimensional CS is easily realized through observation at a random subset of azimuth positions [4,5], while range-dimensional CS falls into the scope of analog sub-Nyquist sampling [6,7].Since CS in range dimension is an important bottleneck in conventional high-resolution SARs, building analog samplers at sub-Nyquist rates is crucial for the sub-Nyquist acquisition of SAR data.In this paper, we therefore focus on the range-dimensional CS.Now, several sub-Nyquist sampling systems have been suggested for sampling radar echoes.The original sub-Nyquist sampler is random sampling [8] or non-uniform sampling [9], which digitizes analog signals at random or irregular time points.However, random/non-uniform sampling is not conducive to hardware implementation and is sensitive to jitter noise.To circumvent these disadvantages, Kirolos et al. established the random demodulation system, which modulates large bandwidth lowpass signals with a pseudo-random binary sequence, followed by lowpass filtering and low-rate ADC sampling [6].Thereafter, Xi et al. merged the random demodulation into classical quadrature sampling and thus presented a quadrature compressive sampling (QuadCS) system to gather intermediate frequency radar echoes [10].On the other hand, Eldar et al. provided a Xampling system [11,12] by accomplishing the random sampling in frequency domain.Among these systems, Xampling has been successfully applied to SAR systems [4].Therefore, these sampling systems are promising in SAR data acquisition at sub-Nyquist rate.
Different from classical SARs via Nyquist sampling, sub-Nyquist SARs collect sub-Nyquist or compressive measurements, and the imaging performance greatly depends on the sensing matrix resulting from the sampling system.Since SAR observations are sequentially conducted, the sub-Nyquist measurements can be modeled as for all N a observations, where l is the index of azimuth position, y l and n l are vectors denoting the Nyquist samples of SAR echo and the sub-Nyquist measurements of observation noise, respectively, and Φ l is the measurement matrix resulting from the sampling system.Stacking all sequential measurement vectors y cs l into a large one, we can translate (1) into a structured CS model as where D is a basis matrix for sparsely representing the Nyquist samples of SAR echo [3], n is the stacked noise vector, and x is the sparse reflecting coefficient vector.An example of the sparse reflecting coefficient is that, in a marine environment, scenes consisting of ships or islands usually lead to sparse SAR images [13,14].In this paper, we focus on this kind of sparse scenes.
In the context of CS, the linear system (2) is typically ill-posed, and image recovery from the sub-Nyquist measurements is equivalent to a sparse recovery problem, which can be handled by many methods, such as greedy algorithms [15] and convex optimization-based algorithms [16].The recovery errors between the ground truth and the estimates returned by these algorithms greatly depend on the sparsity of the vector x and the restricted isometry property (RIP) of the sensing matrix defined as holds for all s-sparse vectors x.In addition, the smallest δ s is called the restricted isometry constant (RIC).If the matrix A satisfies 2s-order RIP with RIC δ 2s < δ # for some constant δ # < 1, the recovery error between the sparse vector x and its estimate x # returned by those algorithms will obey where is a upper bound of the noise norm n .In (4), σ s (x) 1 = inf z 0 ≤s x − z 1 denotes the best s-term approximation error of x in 1 -norm, and C 1 and C 2 are positive constants depending on δ 2s (their values increase as the increase of δ 2s ).Therefore, it is expected that the RIC δ 2s is small for reasonable large s, such that the recovery error x − x # 2 is well bounded, or, in other words, the recovery is stable and robust in the sense that (4) is bounded in the presence of noise and mismodeling (for example, the scene is not reasonable sparse).Unfortunately, previous sub-Nyquist SAR works lack in RIP-based design and analysis.
The original motivation of this paper is to design a RIP-based efficient sampling system for sub-Nyqusit SAR.Here, the RIP-based efficiency means that the sensing matrix A has the smaller RIC, which thus provides better imaging performance with theoretical guarantees.As mentioned above, Xampling has been exploited for sub-Nyquist SAR [4], where a fixed measurement matrix is taken to sample radar echoes for all observations, i.e., the measurement matrices Φ l for l = 1, • • • , N a are the same.However, in the context of CS, the sensing matrix with more randomness usually results in a smaller RIC, so Xampling is not in favor of an efficient sensing matrix.Along this analysis, a natural question arises: may the measurement matrix Φ l vary among different observations so that the sensing matrix A will have smaller RIC?In this paper, we provide an affirmative answer.By designing time-varying measurement matrices, the sampling systems that inject randomness by random modulation are competitive candidates, and they have been well studied, such as random demodulation [6].A variant of random demodulation is QuadCS [10], as shown in Figure 1, which can effectively acquire sub-Nyquist in-phase and quadrature components of sparse bandpass signals.The QuadCS uses random modulation by chipping sequence, bandpass filtering, intermediate-frequency sampling and digital quadrature demodulation to collect sub-Nyquist measurements.It has been successfully applied to sample pulse-Doppler radar echoes with a fixed chipping sequence between different pulse repetition intervals [17][18][19].The architecture of QuadCS.The first subsystem is to perform low-rate random measurement; the second subsystem is a classical digital quadrature demodulation system [20] used for extract compressive in-phase and quadrature measurement.
Based on the QuadCS architecture, this paper proposes an efficient sampling system for sub-Nyquist SAR, which was partially presented in [21].In the sequential observations, the initial seeds that control the chipping sequences are assigned to different values to generate independent chipping sequences.As a result, the measurement matrices Φ l are different from one observation to another, which is significantly different from Xampling [4].By this strategy, more randomness is injected into the sampling system, which better meets the principle of random measurement in CS, and consequently the resulting sensing matrix has a smaller RIC.Theoretical analysis will demonstrate that the RIP attains a diversity gain N a by this strategy, and thus the recovery performance will be improved.In addition, it is important to note that the proposed QuadCS system does not put further burden on hardware implementation because the independence of chipping sequences can be easily set beforehand by maximal-length linear feedback shift registers (MLFSR) with a set of unique initial seeds.
Another contribution of this paper is the development of a fast sparse imaging algorithm for the proposed sub-Nyquist SAR system.As an image inverse problem, the recovery of sparse images from sub-Nyquist measurements have a large problem size.Optimization methods designed for large-scale problems, such as truncated Newton interior-point method [22] and the more popular first-order methods [23][24][25][26], numerically rely on computing the matrix-vector products involving A and its adjoint A * .However, the large problem size precludes the direct application of these off-the-shelf optimization methods because the matrix A cannot be stored explicitly and it is costly, even impractical to calculate the matrix-vector products directly.In the proposed sampling system, we will decompose the matrix-vector products into several fast Fourier transform (FFT) based and Hadamard product-based fast operations.Then, the matrix A does not need to be stored explicitly, and the matrix-vector products have only O (Nlog(N)) complexity and thus much less computational overhead, where N is the Nyquist sampling size of SAR echo.By exploiting the decomposition into typical first-order methods, we provide a first-order fast imaging algorithm.
The rest of this paper is organized as follows.Section 2 reviews some basic knowledge of conventional SAR and sub-Nyquist SAR.Section 3 presents our QuadCS-based sampling system for SAR.Section 4 gives a fast sparse imaging algorithm for sub-Nyquist measurements.Section 5 provides the RIP analysis of proposed system.Section 6 gives several numerical simulations to illustrate the improved performance, and finally we conclude the paper in Section 7.

Notation
The sets of real and complex numbers are denoted by R and C, respectively.An integer set, {1, 2, • • • , N}, is denoted by [N].The cardinality of a set S is denoted by card(S).Vectors and matrices are denoted by lower case boldface and upper case boldface letters, e.g.x and A, respectively.The support of a vector x is denoted by supp(x).The vector output by a reconstruction algorithm is denoted by x # .For a matrix A, vec (A) is the vector obtained by concatenating the columns of A, and diag(A) is the diagonal matrix obtained by taking vec (A) as its diagonal elements.The identity matrix and the zero matrix are denoted by I N and 0 N ∈ R N×N , respectively.
The real part of a complex-valued signal is denoted by Re {•}.The transpose, conjugate and adjoint of a matrix are denoted by (•) T , (•) and (•) * , respectively.The adjoint of an operator is also denoted by (•) * .The 1 , Euclidean and Frobenius norm are denoted by • 1 , • 2 and • F .The symbols ⊗ and • denote the Kronecker product and the Hadamard product, respectively.

SAR Imaging and Sub-Nyquist SAR Model
In this section, we give a basic introduction to SAR imaging, and then discuss the sub-Nyquist SAR and associated image formation method.

SAR Imaging
For simplicity, we take the stripmap SAR as an example.As shown in Figure 2, the radar platform travels along azimuth direction at a constant velocity while the radar illuminates a scene.As the radar moves, the SAR transmits a pulse Re s (t) e j2π f 0 t with carrier frequency f 0 and bandwidth B to illuminate the scene and collect its echo data, and this procedure is sequentially conducted every pulse repetition interval (PRI) T to collect multiple coherent measurements that contains the reflectivity information about the scene.Supposing that the illuminated scene is discretized into K scattering units with predesigned resolution, then the received backscattered signal can be modeled as superpositions of the echo signals due to these scattering units.For the l-th observation, after demodulation, the complex baseband envelope can be expressed as In ( 5), σ n is the complex reflectivity coefficient of the n-th scattering unit and D l,n (t) is a waveform related to the scattering unit as where w a (•) is the antenna beam pattern, η n is the beam center crossing time, and R l,n is the slant range of the n-th scattering unit.For subsequent digital processing, the baseband echo ( 6) is sampled according to the Nyquist theorem to collect its complex sample vector y l with length N r = BT.
A coherent combination of N a echo data forms a two-dimensional data matrix of complex samples Let us denote the reflectivity map of the discretized scene by a complex image matrix X with size K r × K a .By stacking the echo data matrix and image matrix into vectors respectively, we obtain the following observation model where y = vec (Y) is the vectorized Nyquist rate echo data with length N = N r N a , and x = vec (X) is the image vector with length K = K r K a .The SAR image can be well focused by traditional matched-filter-based focusing methods, x m f = D * y, such as range-Doppler algorithm and chirp scaling algorithm (CSA).A virtue of these algorithms is that they are numerically efficiency with O (N log (N)) complexity, owing to range and azimuth decoupled processing and fast Fourier transform (FFT) operations.

Sub-Nyquist SAR
Since its invention, CS theory has been applied to SAR systems to reduce the data rate and subsequent demands on storage [27,28].In classical SAR systems, radar echoes are sampled at Nyqusit rate in azimuth and range dimensions, respectively.Therefore, CS can be conducted both in azimuth and range dimensions.In the azimuth dimension, CS is often realized through randomly transmitting a reduced number of pulses at irregular azimuth positions.A shortcoming of this strategy is that the underlying mathematical model falls into the paradigm of block sampling, which limits the performance of its sensing matrix.In the range dimension, CS of the analog radar signals involves structured sampling systems, such as the Xampling system reported in [4,12].This paper focuses on the sub-Nyquist sampling in range dimension, so the azimuth-dimension CS is not considered.Our work can be naturally extended to the case of azimuth-dimension CS.
In a sub-Nyquist SAR, the sub-Nyquist measurements of the l-th observation can be expressed as where y cs l and n l are the length-M r sub-Nyquist measurement vector and the measurement noise vector, respectively, and Φ l is the M r × N r measurement matrix specified by the sampling system for l-th observation.Let us collect the sample vectors to form a 2D measurement and denote its vectorization by y cs = vec (Y cs ) with length M = M r N a .We can express the sub-Nyquist measurements of (7) as where is the contaminating measurement noise.Typically, the measurement dimension M is much smaller than the image size K, and then image formation from the linear system ( 9) is an ill-posed problem.However, assuming that the image vector x is sparse and the sensing matrix A satisfies the RIP, the image can be estimated by solving the following basis pursuit denoising (BPDN) problem min where λ > 0 is a regularization parameter used to trade off between sparsity and data fidelity.
For this optimization problem, the very high dimension of x precludes the direct application of classical optimization methods, such as inner-point methods and greedy algorithms.However, matrix-vector products involving A and its adjoint A * , if they exist, can be used to speed up many popular algorithms for solving (10), such as fast iterative shrinkage-thresholding algorithm (FISTA) [25].
In this sense, the sampling system is expected to have special structures for fast matrix-vector products.
Then, the recovery algorithm can be adapted for fast image formation.

Quadrature Compressive Sampling for SAR with Independent Measurements
In this section, we present the QuadCS-based sampling system for sub-Nyquist SAR.For sequential SAR observations, we use independent chipping sequences to modulate the echoes; thus, the measurement matrix Φ l is different from one to another.We call this strategy as independent measurements.Then, we provide the frequency-domain representation of this strategy for the convenience of subsequent RIP analysis.As will be shown in Section 5, this strategy will lead to smaller RICs and therefore better recovery performance.

Sampling SAR Echoes via QuadCS
The fundamentals of the QuadCS system [10] are shown in Figure 1, which consists of two subsystems: a sub-Nyquist sampling subsystem and a digital quadrature demodulation subsystem.Before inputting them into the system, the radar signal is down converted to intermediate frequency (IF) f I .With an IF input signal, the first subsystem generates a compressive bandpass signal and its bandpass samples, and then the second subsystem translates the bandpass samples into the complex baseband samples by quadrature demodulation.
Let us take the l-th observation to illustrate its operating principle.In the first subsystem, the IF radar echo Re y l (t)e j2π f I t is modulated by a chipping sequence p l (t) with the bandwidth not less than that of y l (t).The mixing operation will spread the frequency content of the baseband signal y l (t) to full spectrum of p l (t).The mixing output is filtered by a bandpass filter h bp (t) with center frequency f 0 and bandwidth B cs B. The filter outputs the compressive bandpass signal r l (t), where y cs l (t) is the compressive complex envelope given by with I cs l (t) = Re y cs l (t) and Q cs l (t) = Im y cs l (t) denoting the compressive I and Q components of y cs l (t), respectively.The filter output y cs l (t) is then sampled by a low-rate ADC to generate a sub-Nyquist sequence y cs l [k] = r l (k/ f cs ), which is the baseband sample of y cs l (t).The sampling rate f cs is set according to bandpass sampling theorem as f cs = (4 f L + 2B cs )/(4k + 1), where f L = f I − B cs /2 and k is a positive integer satisfying k ≤ f L /2B cs .The second subsystem is to extract digital I and Q sequences from the sub-Nyquist sequence y cs l [k].Its operation is the same as classic quadrature sampling [20].Because of the down-sampling operation, the rate of the digital I and Q sequences, In a PRI, we obtain M r = T/T cs sub-Nyquist complex samples y cs l [m] = I cs l [m] + jQ cs l [m] or 2M r samples of I and Q components, which is much less than 2BT by the classical Nyquist sampling.
Although the QuadCS system works on analog bandpass signals, its output y cs l [m] can be characterized as a linear combination of the sparse coefficient vector (5)  and ( 6) into ( 12), we have where A l,n (t) is the compressive signal of D l,n (t) Then, the sub-Nyquist complex samples can be expressed as To get a standard CS formulation, let us define the sub-Nyquist measurements vector of the l-th observation as collect all N a observations as Y cs = y cs 1 , y cs 2 , • • • , y cs N a or its vector form y cs = vec (Y cs ), and define the sensing matrix A as . Then, we have the sub-Nyquist measurements model y cs = Ax as given by ( 9) in the noisy case.However, the structure of the sensing matrix A here is not as clear as that in (2).To decompose the sensing matrix into an explicit structured form, we next analyze the QuadCS in frequency domain, which is also important for subsequent discussions on the fast imaging algorithm and RIP.

Frequency-Domain Representation
From (11), we can see that the QuadCS involves a time-domain product and a bandpass filtering operation.The time-domain product is equivalent to frequency-domain convolution, which can be described by the product of a Toeplitz matrix by a vector, and the bandpass filtering can be modelled by a truncation operation in frequency domain under the assumption of idea frequency response.Based on the above-mentioned view, we have the frequency-domain measurement model of QuadCS as studied in [10].
Since each observation interval has a duration T, the chipping sequences p l (t) can be expressed in Fourier series as where f p = 1/T, ρ i,l is the Fourier coefficient of p l (t), and L p ≥ BT is a positive integer.Denoting the discrete Fourier transform (DFT) of y cs l by ŷcs l , it can be given by [29] where F r is a N r × N r normalized DFT matrix after fftshift (fftshift is to shift zero-frequency component to the center of spectrum), and T l is a M r × N r Toeplitz matrix defined by the Fourier coefficients as where The Toeplitz structure of T l results from the convolution between the spectrums of the chipping sequence and the input signal.Due to bandpass filtering, T l does not contain all Fourier coefficients of p l (t).Next, according to (19), we can describe the time-domain sub-Nyquist measurements as where F m is a M r × M r normalized DFT matrix after fftshift.Then, the measurement matrix Φ l can be described as which is independent from one observation to another due to the adoption of different chipping sequences.Finally, according to (8) and ( 9), we can rewrite the sensing matrix as a structured form, The first factor of A is a block diagonal matrix with independent diagonal blocks.Therefore, the sensing matrix A has more randomized structures, which better meets the principle of random measurement in CS.In addition, the block diagonal and Toeplitz structures are also helpful for deriving a fast sparse recovery algorithm, as will be shown later.

Remarks
The QuadCS generates independent measurement matrices while does not complicate hardware implementation.In the system, chipping sequences are realized by a pseudo-random number generator, which can be implemented using a MLFSR.The MLFSR has the benefit of providing a random sequence of ±1 with zero average, and it offers the ability of regenerating the same sequence by a given initial seed [6].This technique has been widely used in modern radio communication, such as CDMA.During sequential observations of SAR, the initial seeds are easily assigned to different values, which can be set in advance or in real-time.On the contrary, in the Xampling system [4,12], it is not allowed to generate varied measurement matrices among different observations, because four bandpass filters with different passbands are utilized to perform sub-Nyquist measurement and the measurement matrix is inherently fixed when the hardware is built.Moreover, bandpass filters cannot sample SAR echoes randomly enough in frequency domain and thus often generate bad sensing matrices, as observed in the following simulations.As for non-uniform sampling, it often imposes strict timing requirements on the digital circuits to preserve timing information [8,30], which makes it difficult to interface with synchronous digital circuits.In fact, it is easier to build a high-rate modulator via mixer/chipping sequence combination than a specially designed non-uniform sampler [31].

Fast Sparse Imaging with QuadCS Measurements
As a standard sparse recovery problem in CS, we model our imaging problem as the BPDN formulation in (10).For the large-scale problem, the very high dimensions of A often precludes the direct application of off-the-shelf optimization methods because of two limitations: the matrix A cannot be stored explicitly, and the matrix-vector products involving A and its adjoint A * cannot be quickly calculated.If these two limitations are well resolved, several popular solvers can be used to solve (10), such as FISTA [32], SpaRSA [33], TFOCS [26] and L1LS [22].In the proposed sampling system, we exploit the block and Toeplitz structure of the sensing matrix A and decompose the matrix-vector products into several fast operations.Then, the matrix A does not need to be stored explicitly, and the matrix-vector products can be fast calculated.By inserting the decomposition into above-mentioned solvers, we derive a FISTA-based fast SAR imaging algorithm for the proposed sub-Nyquist SAR system, and analyze the algorithm complexity.

Fast Matrix-Vector Products
As we have discussed, the sensing matrix A = ΦD consists of two parts: measurement matrix Φ and basis matrix D. To obtain the desired fast matrix-vector products for A and A * , we are going to exploit the inner structures of the two matrices.The following content is to expand the derivations in detail.
Let us first discuss the Toeplitz structure as shown in the frequency-domain model (21).Since T l is a Toeplitz matrix, it can be viewed as a submatrix of a (2L 0 + 1) × (2L 0 + 1) circulant matrix C l defined as Specifically, T l is exactly the right N r columns and top M r rows of C l .Let us define row and column truncation matrices as and respectively.Then, we embed T l into C l as follows: Since circulant matrices can be diagonalized by DFT matrices, there exist fast matrix-vector products for circulant matrices via FFT and the computational complexity is O(Nlog(N)).The goal of ( 26) is to extend this property to the Toeplitz matrix T l .Specifically, let us define an integer L = 2L 0 + 1, and a vector T which is the first column vector of C l , and denote the L × L DFT matrix as F L .Then, the circulant matrix C l can be diagonalized as where θ l = F L ρ l .Inserting ( 26) and ( 27) into ( 21), we have then the 2D QuadCS measurement Y cs can be written as which mainly involves FFTs and thus has a low complexity.Note that, from (29), we find that the random matrix Θ performs "Hadamard scrambling" to inject randomness into the measurements.Next, we analyze the inner structure of the basis matrix D. It has been suggested that conventional imaging algorithms, such as CSA, can be exploited to derive its fast matrix-vector product for D [5,34], which demonstrate that the matrix is also highly structured.We assume that the SAR operates in stripmap mode and the transmitted pulses are linear frequency modulated (LFM); then, CSA can be used to obtain the SAR image, which can be described in a compact form as In this equation, X CSA is the SAR image output by CSA, H 1 , H 2 , H 3 are phase compensation matrices (i.e., their entries lie in the unit circle in C, for more details, cf.chapter 6 in [35]), and F r , F a are the range and azimuth DFT matrices, respectively.Because CSA is a high precision algorithm, X CSA is a valid approximation of the reflectivity map X.For ease of discussion, we omit the approximation error, i.e., we take X = X CSA , and in this sense, X and Y have the same size, i.e., K r = N r and K a = N a .After simple derivation, the inversion of ( 30) is given by [5] It is easy to see that the linear system (31) is equivalent to y = Dx with proper reshape operations.Furthermore, since the involved FFTs and Hadamard products in (30) and (31) are unitary operators, the matrix D is also unitary [5].

FISTA-Based Fast Imaging Algorithm
FISTA is an accelerated variant of the iterative shrinkage-thresholding algorithms (ISTA), and it has a quadratic rate of convergence thanks to the generalized Nesterov's method [23,25].Each iteration of FISTA involves matrix-vector product involving A and its adjoint A * followed by a shrinkage/soft-threshold step as follows: where x k+1 , z k are iterative variables, L f is the Lipschitz constant for the gradient ∆ f of the function , and T λ/L f (•) is the soft shrinkage operator defined as From (32), we know that the matrix-vector products involving A and its adjoint A * dominate the computational cost in each iteration.The previous subsection shows that the matrix-vector products can be fast calculated through matrix decomposition.Now, we merge this calculation strategy into FISTA and thus provide a FISTA-based fast imaging algorithm, as shown in Algorithm 1.

Input:
Operators M, D and their adjoints; QuadCS measurement Y cs ; number of iterations Niter; λ; Output: Sparse SAR image X # ; 1: initial: take λ, 5: For ease of description, we denote the right sides of ( 29) and ( 31) by operators as MY and DX, respectively, and denote vec(•) by an operator V.With these notations, we can express the matrix-vector product of A as Ax = V MDV * x. (34) Similarly, the matrix-vector product of the adjoint A * is given by where the adjoint operator D * is given by the operations on the right side of ( 30) and M * is given by It is clear that the involved operations in the right sides of ( 34) and ( 35) mainly consist of FFTs and Hadamard products.Therefore, we can conclude that A and its adjoint A * both admit fast matrix-vector products.

Complexity
Without fast matrix-vector products for A and A * , we need to store the M × N matrix A explicitly and compute the involved matrix-vector products directly.For a small scene of 1024 × 1024 pixels and compression ratio M/N = 1/8, it requires 1024 GB to store the matrix A with double precision, and roughly 1.3744 × 10 11 complex multiplications to compute the matrix-vector product Ax.Obviously, the image recovery will be numerically intractable.With the fast matrix-vector multiplications (34) and (35), we only need to store N r × N a matrices H 1 , H 2 , H 3 and Θ.It only requires 32 MB to store the four matrices and roughly 8.2 × 10 4 complex multiplications to compute Ax for a scene of 1024 × 1024 size, which are significantly smaller than the previous requirements.In general, the storage cost is reduced form O(MN) to O(N), and the computational cost is reduced from O(MN) to O(Nlog(N)) per iteration, respectively.
Looking back to the traditional imaging algorithm, CSA, it only needs to evaluate (30) once to form a SAR image.Specifically, it involves 2N r N a (log (N r ) + log (N a )) + 3N r N a complex products, or in a compact form, O (N log (N)) complexity.The proposed iterative algorithm exhausts more overhead compared to CSA, depending on the number of iterations.Indeed, the essence of CS is to sacrifice the computational cost in return for the reduction of sampling/data rate.

Restricted Isometry Property (RIP) Analysis
In this section, we derive the RIP of the sensing matrix A resulting from the proposed sampling system.The established RIP condition shows that the strategy of independent measurements provides a "diversity gain", which is up to N a .Therefore, the sensing matrix A has smaller RICs with independent measurements compared to that with fixed measurements (i.e., with equal chipping sequences), and thus leads to lower sampling/data rate and smaller recovery error.
Our derivation is based on the concentration of measure (CoM) inequality [36], which is one of the leading techniques for the RIP analysis of random sensing matrices in CS.The CoM inequalities have been established for the random matrix with i.i.d.entries and other structured matrices.However, these results are not applicable to the proposed sensing matrix A because of its special factornizatio.
In the following, we first derive a CoM inequality of the matrix A, and then establish a RIP condition on it.The CoM inequality of A is formulated in Lemma 1.

Lemma 1.
Let x be a s-sparse vector in C N .For an arbitrary 0 < δ < 1, we have where M is the number of sub-Nyquist measurements, χ denotes two settings of chipping sequences and is defined as ξ is the duty cycle of radar transmitted signals, and is the maximum value of the normalized antenna beam pattern defined as where T a is the antenna footprint time.
Proof.The proof is given in Appendix A.1.
Lemma 1 points out that, with independent measurements, the probabilistic tail bounds of the CoM decays exponentially N a times faster than that with fixed measurements.Since the CoM inequalities quantify how well a random matrix preserve the norm of a high-dimensional signal in the mapped low-dimensional space [37], the independent measurements result in a better distance-preserving property for the sensing matrix A, which is the essence of RIP.
Based on the CoM inequality in Lemma 1, we can derive a uniform tail bound of A * S A S − I s 2 for a fixed subset S ⊂ [N] with card (S) = s using the covering number estimate [38].Then, we take the union bound over all subsets S ⊂ [N] of cardinality s, to derive the following RIP condition.
Theorem 1.The restricted isometry constant δ s of the sensing matrix A satisfies δ s < δ for some constant 0 < δ < 1 with probability exceeding 1 − η provided that M ≥ 32χs Proof.The proof is given in Appendix A.2.
Theorem 1 gives the RIP condition on the sensing matrix A. The key parameters in (41) are M and χ, which describe the measurement amount and the setting of chipping sequences, respectively.If M (or equivalently, compression ratio) is properly chosen according to the sparsity s and the imaged scene size N, the sensing matrix A will satisfy RIP with high probability in both settings of chipping sequences.However, the measurement amounts required by the two settings vary somewhat surprisingly, up to a multiplicative factor N a .This variance suggests that the strategy of independent measurements provides a diversity gain, which is similar to the waveform diversity gain in MIMO radar.For a given imaged scene, the diversity gain results in lower sampling/data rates.On the other hand, the parameters ξ and are often fixed in a given SAR operation mode, and they suggest that the duty cycle and the antenna beam pattern have a direct influence on the RIP.Intuitively, for smaller ξ or , each column of the basis D will be localized in smaller regions, which is bad for the sensing matrix.Another benefit of the diversity gain is that, with the same number of measurements, the RIC is expected to be much smaller, as stated in the following corollary.
Corollary 1.Let δ s and γ s be the s-order RICs of the sensing matrices for the setting of equal and independent chipping sequences, respectively.For fixed M and η, the RICs satisfy δ s < δ and γ s < δ/N a simultaneously for some constant 0 < δ < 1 with probability exceeding 1 − η.
Replacing δ s by δ s /χ in Theorem 1 will prove the corollary.According to CS theory, smaller RIC means smaller recovery error.Thus, independent measurements have better performance than fixed measurements, as revealed by Theorem 2.
Theorem 2. For a fixed number of measurements M and a fixed probabilistic guarantee 1 − η, let s 1 and s 2 be the maximum integers that satisfy δ 2s 1 < δ # and γ 2s 2 < δ # with equal and independent chipping sequences, respectively, where δ # = √ 2 − 1.Let x # 1 and x # 2 be the solutions of the constrained BPDN problem in the two settings of chipping sequences (39), respectively, then where the last line is exactly the error bound of x − x # 1 2 .
Proof.The proof is shown in Appendix A.3.
Theorem 2 makes a comparison between the recovery error bounds about the two settings of chipping sequences.In (43), the error bound consists of two parts: one is the best s-term approximation error and the other comes from the measurement noise.As shown in Appendix A.3, we have a gain factor s 2 /s 1 ≈ √ N a , which results in much lower best s-term approximation error as shown in (43).Therefore, we can recover more sparse targets in the case of high signal-to-noise ratio (SNR).For example, when N a = 4096, the gain factor will be 64.Thus, the stability of image recovery is significantly improved provided that the scene is not reasonable sparse.When the SNR decreases, the noisy contribution C 2 gradually dominates the error bounds, and consequently, the relative errors get closer.Namely, the advantage of QuadCS with independent measurement is weakened in this case.However, there is no doubt that x # 2 always has a smaller recovery error than x # 1 .Note that we recover images by solving (10) instead of (42).In fact, the two problems are equivalent.If x # is a minimizer of the (42) with some > 0, there necessarily exists λ ≥ 0 such that x # is a minimizer of (10) [39,40].It should be pointed out that is related to the noise level and the optimal choice of λ is related to x # and .In this paper, we assume that the noise level is known for (42), and λ is empirically chosen in (10) according to the noise level.In real SAR applications, the noise level can be measured when no SAR echo arrives at the radar receiver.In future works, we will consider how to deal with unknown noise level adaptively.
In summary, the above theoretical analysis shows that, using independent measurements, the proposed sampling system has smaller recovery error and enables lower sampling/data rate.This conclusion is of great importance since our goal is to achieve low-rate sampling and meanwhile provide robust recovery performance for sub-Nyquist SAR.

Simulations
In this section, we conduct several simulations to evaluate the proposed QuadCS-based SAR imaging performance.Results on both synthetic data and real SAR data collected by TerraSAR-X and RADARSAT-1 missions are presented.Comparisons are made between Xampling and QuadCS with equal or independent chipping sequences.The relative error, x − x # 2 / x 2 , and the relative root-mean-square error (RRMSE), E x − x # 2 / x 2 , are used as the imaging performance metrics.The simulation scheme is shown in Figure 3 and main radar parameters are summarized in Table 1.Given a synthetic image, real SAR single-look complex (SLC) image (as shown in Figure 4) or a raw data matrix, we use the operations (29) or (34) to generate sub-Nyquist measurements.Similar methods for SAR raw data generation can be found in [41,42].To model the measurement noise, an additive white Gaussian noise (AWGN) is added.The regularization parameter and the maximum number of iterations are empirically set to 10 −3 and 200, respectively.For visualization, the recovered images are mapped to grayscale images by the following map where τ is a threshold for preventing the image from being almost black due to the presence of strong scattering points.A default value provided by a SAR image processing toolbox (code "plotimage_sar.m" in PPB toolbox [43]) is adopted in our simulations, which is given by the sum of mean and triple standard deviation of the amplitude image.For ease of description, we refer to QuadCS-IndSeq and QuadCS-EquSeq as QuadCS with independent and equal chipping sequences, respectively, and denote the compression ratio M/N by α.

Simulations with Synthetic Data
Firstly, we use some synthetic geometries as a SLC image, as shown in Figure 4a, to evaluate the imaging performance.The compression ratio is set to 1/8 and 1/16, respectively, and the SNR equals 10 dB.As shown in the first column of Figure 5, QuadCS-IndSeq recovers these geometries very well, while QuadCS-EquSeq and Xampling cannot reconstruct horizontal lines clearly.Namely, QuadCS-EquSeq and Xampling cannot recover images in range dimension because of low sampling rate and the usage of fixed measurement matrix in each observation.At the smaller compression ratio 1/16, QuadCS-IndSeq can still recover all geometries with some performance loss, but QuadCS-EquSeq and Xampling cannot reconstruct the geometries correctly.Obviously, the proposed sampling system has the best performance.
Next, we investigate the relationship between RRMSE and the image sparsity under different compressive ratios and SNRs.We choose random subsets of the test image and set its entries to be i.i.d.random variables obeying uniform distribution on [0, 1].The image recovery conducts 100 times and the resulting RRMSE curves are shown in Figure 6.All the figures show that, as the sparsity increases, the recovery errors of all sampling schemes increases.This can be explained by ( 4) that the best s-term approximation error increases as the sparsity increases.On the other hand, the figures also show that QuadCS-IndSeq has the lowest recovery error in all simulation cases.This is because QuadCS-IndSeq has better RIP, as disclosed in Theorem 2. Furthermore, QuadCS is superior to Xampling in all cases, no matter what the independent measurement strategy is utilized.The reason is that Xampling conducts frequency-domain random sampling by several bandpass filters, so that the randomness is not enough and thus the RIP of corresponding sensing matrix is not very well.Hence, QuadCS is superior to Xampling for sub-Nyqusit SAR.According to each row subfigures, we find that QuadCS-IndSeq has more performance advantages in the high SNR case (SNR = 20 dB), since the RMMSE relies on the best s-term approximation error in this case.Specifically, as seen from the middle left subfigure in Figure 6, where SNR is 20 dB and compression ratio is 1/16, when the sparsity is 13‰, the RRMSE of QuadCS-IndSeq is about −14.2 dB, while those of QuadCS-EquSeq and Xampling are about −8.8 dB and −4.4 dB, respectively.On the other hand, according to each column subfigure, we know that the image recovery performance degrades faster with the lower compressive ratio, which is the inherent characteristic of sub-Nyquist sampling systems.However, QuadCS-IndSeq still has the best image recovery performance.In summary, our sampling system considerably improves the imaging performance under relatively large sparsity and low sampling rate.However, with the decrease of SNR, the recovery errors mainly come from the noise, as shown in (43); then, the improved performance is gradually diluted.

Simulations with Real SAR Data
Now, we use a real SLC image, shown in Figure 4b, to evaluate the proposed sampling system.The image with size 2048 × 2048 is collected by TerraSAR-X on 4 July 2008, which contains several ships near the coast of Barcelona.For a clear display, we show 200 × 200 area of the recovered images in Figure 7.We find that QuadCS-IndSeq has a smaller recovery error than QuadCS-EquSeq and Xampling.
We also use RADARSAT-1 raw data collected on 16 June 2002 to evaluate the proposed system.Imaging results using the raw data is shown in Figure 4c.Because the scene is not reasonably sparse, we consider larger settings of compression ratios, namely, 1/4 and 1/8 , respectively.Figure 8 shows a 500 × 500 area of the recovered images.As we noted in Section 1, the recovery error greatly depends on the sparsity of the image.Due to lack of sparsity, the recovery errors are large for all simulated schemes as shown in Figure 8.However, QuadCS-IndSeq still has the lowest recovery errors.B S = x ∈ C N , supp (x) ⊂ S, x 2 ≤ 1 , which satisfies card (Z) ≤ (1 + 1/ρ) 2s and min z∈Z x − z 2 ≤ ρ, for all x ∈ B S .Let us choose one element from each subset Z i , i.e., z i ∈ Z i for i = 1, 2, • • • , card (Z).According to Lemma 1, the concentration inequality gives for t ∈ (0,

(A35)
Taking the maximum over all x ∈ B S gives  [38], taking the union bound over all ( N s ) subsets S gives

Figure 1 .
Figure 1.The architecture of QuadCS.The first subsystem is to perform low-rate random measurement; the second subsystem is a classical digital quadrature demodulation system[20] used for extract compressive in-phase and quadrature measurement.

Figure 3 .
Figure 3. Block diagram of simulations.We take synthetic, real images and real raw data as input data, respectively.The simulated data output by QuadCS and Xampling are processed by a sparse recovery algorithm.

Figure 4 .
Figure 4. Scenes of interest in simulations.(a) a synthetic image; (b) a real SAR SLC image, collected by TerraSAR-X on 1 May 2008, and the imaged area is near the coast of Barcelona; (c) an SLC image focused from RADARSAT-1 raw data and the red-framed area is the Tsawwassen Ferry Port.

Figure 5 .
Figure 5. Recovery of the synthetic image in 10 dB AWGN.From top to bottom: QuadCS with independent chipping sequences (Quad-IndSeq), QuadCS with equal chipping sequences (Quad-EquSeq), Xampling.From left to right: recovered images and corresponding residuals with compression ratios 1/8 and 1/16.Range and azimuth are along the horizontal and vertical axes, respectively.The relative errors ( measured in dB) of these images are shown at the lower left corners of their residual images, respectively.

Figure 7 .Figure 8 .
Figure 7. Recovery of the TerraSAR-X SLC image.From top to bottom: QuadCS-IndSeq, QuadCS-EquSeq, Xampling.From left to right: recovered images and corresponding residuals with compression ratios 1/8 and 1/16.Range and azimuth are along the horizontal and vertical axes, respectively.The relative errors for these images are shown at the lower left corners of their residual images, respectively.For the sparse scene, QuadCS-IndSeq achieves −5.7 dB even at the compression ratio 1/16.
1), For a proper choice of ρ, t and δ, we prove that (A34) implies Ax2 2 − x 2 2 ≤ δ for all x ∈ B S , i.e., A * S A S − I s 2 ≤ δ.For simplicity, denote Q = A * S A S − I s .Now, for any x ∈ B S , there exists a vector z ∈ {z i |i = 1, • • • , card (Z)} satisfying x − z 2 ≤ ρ ≤ 1/2.It then follows that Let us choose t = (1 − 2ρ)δ, such that Q 2 < δ.According to (A34), we haveP A * S A S − I s 2 ≥ δNext, we derive the RIP condition.Since δ s can be equivalently described asδ s = sup S⊂[N], card(S)=s A * S A S − I s 2