Sub-Nyquist SAR via Quadrature Compressive Sampling with Independent Measurements

Abstract

This paper presents an efficient 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, filtering, 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 significantly 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 first-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.

1. 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

$${\mathbf{y}}_l^{cs}={\mathsf{\Phi }}_l{\mathbf{y}}_l+{\mathbf{n}}_l,l=1,\dots ,N_a$$

(1)

for all $N_a$ observations, where l is the index of azimuth position, ${\mathbf{y}}_l$ and ${\mathbf{n}}_l$ are vectors denoting the Nyquist samples of SAR echo and the sub-Nyquist measurements of observation noise, respectively, and ${\mathsf{\Phi }}_l$ is the measurement matrix resulting from the sampling system. Stacking all sequential measurement vectors ${\mathbf{y}}_l^{cs}$ into a large one, we can translate (1) into a structured CS model as

$${\mathbf{y}}^{cs}=\mathrm{diag}\left\{{\mathsf{\Phi }}_1,\cdots ,{\mathsf{\Phi }}_{N_a}\right\}\mathbf{D}\mathbf{x}+\mathbf{n},$$

(2)

where $\mathbf{D}$ is a basis matrix for sparsely representing the Nyquist samples of SAR echo [3], $\mathbf{n}$ is the stacked noise vector, and $\mathbf{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 $\mathbf{x}$ and the restricted isometry property (RIP) of the sensing matrix defined as $\mathbf{A}=\mathrm{diag}\{{\mathsf{\Phi }}_1,\cdots ,{\mathsf{\Phi }}_{N_a}\}\mathbf{D}$. Specifically, the matrix $\mathbf{A}$ satisfies the s-order RIP if there exists a constant ${\delta }_s\in \left(0,1\right)$ such that

$$\left(1-{\delta }_s\right){∥\mathbf{x}∥}_2^2⩽{∥\mathbf{A}\mathbf{x}∥}_2^2⩽\left(1+{\delta }_s\right){∥\mathbf{x}∥}_2^2$$

(3)

holds for all s-sparse vectors $\mathbf{x}$. In addition, the smallest ${\delta }_s$ is called the restricted isometry constant (RIC). If the matrix $\mathbf{A}$ satisfies $2s$-order RIP with RIC ${\delta }_{2s}<{\delta }^{\#}$ for some constant ${\delta }^{\#}<1$, the recovery error between the sparse vector $\mathbf{x}$ and its estimate ${\mathbf{x}}^{\#}$ returned by those algorithms will obey

$${∥\mathbf{x}-{\mathbf{x}}^{\#}∥}_2\le C_1\frac{{\sigma }_s{\left(\mathbf{x}\right)}_1}{\sqrt{s}}+C_2ϵ,$$

(4)

where $ϵ$ is a upper bound of the noise norm $\parallel \mathbf{n}\parallel$. In (4), ${\sigma }_s{\left(\mathbf{x}\right)}_1={inf}_{{∥\mathbf{z}∥}_0\le s}{∥\mathbf{x}-\mathbf{z}∥}_1$ denotes the best s-term approximation error of $\mathbf{x}$ in ${\ell }_1$-norm, and $C_1$ and $C_2$ are positive constants depending on ${\delta }_{2s}$ (their values increase as the increase of ${\delta }_{2s}$). Therefore, it is expected that the RIC ${\delta }_{2s}$ is small for reasonable large s, such that the recovery error ${∥\mathbf{x}-{\mathbf{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 $\mathbf{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 ${\mathsf{\Phi }}_l$ for $l=1,\cdots ,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 ${\mathsf{\Phi }}_l$ vary among different observations so that the sensing matrix $\mathbf{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].

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 ${\mathsf{\Phi }}_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 $\mathbf{A}$ and its adjoint ${\mathbf{A}}^{*}$. However, the large problem size precludes the direct application of these off-the-shelf optimization methods because the matrix $\mathbf{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 $\mathbf{A}$ does not need to be stored explicitly, and the matrix-vector products have only $\mathcal{O}\left(N\log \left(N\right)\right)$ 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 $\mathbb{R}$ and $\mathbb{C}$, respectively. An integer set, $\left\{1,2,\cdots ,N\right\}$, is denoted by $\left[N\right]$. The cardinality of a set S is denoted by $\mathrm{card}\left(S\right)$. Vectors and matrices are denoted by lower case boldface and upper case boldface letters, e.g. $\mathbf{x}$ and $\mathbf{A}$, respectively. The support of a vector $\mathbf{x}$ is denoted by $\mathrm{supp}\left(\mathbf{x}\right)$. The vector output by a reconstruction algorithm is denoted by ${\mathbf{x}}^{\#}$. For a matrix $\mathbf{A}$, $\mathrm{vec}\left(\mathbf{A}\right)$ is the vector obtained by concatenating the columns of $\mathbf{A}$, and $\mathrm{diag}\left(\mathbf{A}\right)$ is the diagonal matrix obtained by taking $\mathrm{vec}\left(\mathbf{A}\right)$ as its diagonal elements. The identity matrix and the zero matrix are denoted by ${\mathbf{I}}_N$ and ${\mathbf{0}}_N\in {\mathbb{R}}^{N\times N}$, respectively.

The real part of a complex-valued signal is denoted by $\mathrm{Re}\left\{·\right\}$. The transpose, conjugate and adjoint of a matrix are denoted by ${\left(·\right)}^T$, $\overline{\left(·\right)}$ and ${\left(·\right)}^{*}$, respectively. The adjoint of an operator is also denoted by ${\left(·\right)}^{*}$. The ${\ell }_1$, Euclidean and Frobenius norm are denoted by ${∥·∥}_1$, ${∥·∥}_2$ and ${∥·∥}_F$. The symbols ⊗ and ∘ denote the Kronecker product and the Hadamard product, respectively.

2. 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.

2.1. 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 $\mathrm{Re}\left\{s\left(t\right)e^{j2\pi f_{0}t}\right\}$ 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

$$y_l\left(t\right)=\sum _{n=1}^K{D}_{l,n}\left(t\right){\sigma }_n.$$

(5)

In (5), ${\sigma }_n$ is the complex reflectivity coefficient of the n-th scattering unit and $D_{l,n}\left(t\right)$ is a waveform related to the scattering unit as

$$\begin{array}{cc}\hfill D_{l,n}\left(t\right)=& w_a\left(lT-{\eta }_n\right)s\left(t-2R_{l,n}/\phantom{2R_{l,n}c}c\right)\hfill \\ & \times \exp \left(-j4\pi f_0{R}_{l,n}/\phantom{-j4\pi f_0{R}_{l,n}c}c\right),\hfill \end{array}$$

(6)

where $w_a\left(·\right)$ is the antenna beam pattern, ${\eta }_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 ${\mathbf{y}}_l$ with length $N_r=BT$. A coherent combination of $N_a$ echo data forms a two-dimensional data matrix of complex samples $\mathbf{Y}=\left[{\mathbf{y}}_1,\cdots ,{\mathbf{y}}_{N_a}\right]$. Let us denote the reflectivity map of the discretized scene by a complex image matrix $\mathbf{X}$ with size $K_r\times K_a$. By stacking the echo data matrix and image matrix into vectors respectively, we obtain the following observation model

$$\mathbf{y}=\mathbf{D}\mathbf{x},$$

(7)

where $\mathbf{y}=vec\left(\mathbf{Y}\right)$ is the vectorized Nyquist rate echo data with length $N=N_r{N}_a$, and $\mathbf{x}=vec\left(\mathbf{X}\right)$ 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, ${\mathbf{x}}_{mf}={\mathbf{D}}^{*}\mathbf{y}$, such as range-Doppler algorithm and chirp scaling algorithm (CSA). A virtue of these algorithms is that they are numerically efficiency with $\mathcal{O}\left(Nlog\left(N\right)\right)$ complexity, owing to range and azimuth decoupled processing and fast Fourier transform (FFT) operations.

2.2. 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

$${\mathbf{y}}_l^{cs}={\mathsf{\Phi }}_l{\mathbf{y}}_l+{\mathbf{n}}_l,l\in \left[N_a\right],$$

(8)

where ${\mathbf{y}}_l^{cs}$ and ${\mathbf{n}}_l$ are the length-$M_r$ sub-Nyquist measurement vector and the measurement noise vector, respectively, and ${\mathsf{\Phi }}_l$ is the $M_r\times N_r$ measurement matrix specified by the sampling system for l-th observation. Let us collect the sample vectors to form a 2D measurement ${\mathbf{Y}}^{cs}=\left[{\mathbf{y}}_1^{cs},\cdots ,{\mathbf{y}}_{N_a}^{cs}\right]$, and denote its vectorization by ${\mathbf{y}}^{cs}=\mathrm{vec}\left({\mathbf{Y}}^{cs}\right)$ with length $M=M_r{N}_a$. We can express the sub-Nyquist measurements of (7) as

$${\mathbf{y}}^{cs}=\mathbf{A}\mathbf{x}+\mathbf{n},$$

(9)

where $\mathbf{A}=diag\left({\mathsf{\Phi }}_1,\cdots ,{\mathsf{\Phi }}_{N_a}\right)\mathbf{D}\in {\mathbb{C}}^{M\times K}$ is the sensing matrix and $\mathbf{n}=\mathrm{vec}\left([{\mathbf{n}}_1,\cdots ,{\mathbf{n}}_{N_a}]\right)\in {\mathbb{C}}^{M\times 1}$ 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 $\mathbf{x}$ is sparse and the sensing matrix $\mathbf{A}$ satisfies the RIP, the image can be estimated by solving the following basis pursuit denoising (BPDN) problem

$$\underset{\mathbf{x}}{min}\frac{1}{2}{∥\mathbf{A}\mathbf{x}-{\mathbf{y}}^{cs}∥}_2^2+\lambda {∥\mathbf{x}∥}_1,$$

(10)

where $\lambda >0$ is a regularization parameter used to trade off between sparsity and data fidelity. For this optimization problem, the very high dimension of $\mathbf{x}$ precludes the direct application of classical optimization methods, such as inner-point methods and greedy algorithms. However, matrix-vector products involving $\mathbf{A}$ and its adjoint ${\mathbf{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.

3. 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 ${\mathsf{\Phi }}_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.

3.1. 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 $\mathrm{Re}\left\{y_l\left(t\right)e^{j2\pi f_{I}t}\right\}$ is modulated by a chipping sequence $p_l\left(t\right)$ with the bandwidth not less than that of $y_l\left(t\right)$. The mixing operation will spread the frequency content of the baseband signal $y_l\left(t\right)$ to full spectrum of $p_l\left(t\right)$. The mixing output is filtered by a bandpass filter $h_{bp}\left(t\right)$ with center frequency $f_0$ and bandwidth $B_{cs}\ll B$. The filter outputs the compressive bandpass signal $r_l\left(t\right)$,

$$\begin{array}{cc}\hfill r_l\left(t\right)=& {\int }_{-\infty }^{\infty }{h}_{bp}\left(t-\tau \right)p_l\left(\tau \right)\mathrm{Re}\left\{y_l\left(\tau \right)e^{j2\pi f_I\tau }\right\}d\tau \hfill \\ & =Re\left\{y_l^{cs}\left(t\right)e^{j2\pi f_{I}t}\right\},\hfill \end{array}$$

(11)

where $y_l^{cs}\left(t\right)$ is the compressive complex envelope given by

$$y_l^{cs}\left(t\right)={\int }_{-\infty }^{\infty }{h}_{bp}\left(t-\tau \right)e^{-j2\pi f_I(t-\tau )}{p}_l\left(\tau \right)y_l\left(\tau \right)d\tau ,$$

(12)

with $I_l^{cs}\left(t\right)=Re\left\{y_l^{cs}\left(t\right)\right\}$ and $Q_l^{cs}\left(t\right)=Im\left\{y_l^{cs}\left(t\right)\right\}$ denoting the compressive I and Q components of $y_l^{cs}\left(t\right)$, respectively. The filter output $y_l^{cs}\left(t\right)$ is then sampled by a low-rate ADC to generate a sub-Nyquist sequence $y_l^{cs}\left[k\right]=r_l\left(k/f_{cs}\right)$, which is the baseband sample of $y_l^{cs}\left(t\right)$. The sampling rate $f_{cs}$ is set according to bandpass sampling theorem as $f_{cs}=\left(4f_L+2B_{cs}\right)/\left(4k+1\right)$, where $f_L=f_I-B_{cs}/2$ and k is a positive integer satisfying $k\le ⌊f_L/2B_{cs}⌋$. The second subsystem is to extract digital I and Q sequences from the sub-Nyquist sequence $y_l^{cs}\left[k\right]$. 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, $I_l^{cs}\left[m\right]=I_l^{cs}\left(m{T}_{cs}\right)$ and $Q_l^{cs}\left[m\right]=Q_l^{cs}\left(m{T}_{cs}\right)$, is half that of $y_l^{cs}\left[k\right]$, $T_{cs}=2/f_{cs}$. In a PRI, we obtain $M_r=⌊T/T_{cs}⌋$ sub-Nyquist complex samples $y_l^{cs}\left[m\right]=I_l^{cs}\left[m\right]+j{Q}_l^{cs}\left[m\right]$ 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_l^{cs}\left[m\right]$ can be characterized as a linear combination of the sparse coefficient vector $\mathbf{x}={\left[{\sigma }_1,\cdots ,{\sigma }_N\right]}^T$. Inserting (5) and (6) into (12), we have

$$\begin{array}{cc}\hfill y_l^{cs}\left(t\right)& =\sum _{n=1}^K{\sigma }_n{\int }_{-\infty }^{\infty }{h}_{bp}\left(t-\tau \right)e^{-j2\pi f_I(t-\tau )}\hfill \\ & \times p_l\left(\tau \right)D_{l,n}\left(\tau \right)d\tau \hfill \\ & =\sum _{n=1}^K{\sigma }_n{A}_{l,n}\left(t\right),\hfill \end{array}$$

(13)

where $A_{l,n}\left(t\right)$ is the compressive signal of $D_{l,n}\left(t\right)$

$$A_{l,n}\left(t\right)={\int }_{-\infty }^{\infty }{h}_{bp}\left(t-\tau \right)e^{-j2\pi f_I(t-\tau )}{p}_l\left(\tau \right)D_{l,n}\left(\tau \right)d\tau .$$

(14)

Then, the sub-Nyquist complex samples can be expressed as

$$y_l^{cs}\left[m\right]=\sum _{n=1}^K{\sigma }_n{A}_{l,n}\left(m{T}_{cs}\right).$$

(15)

To get a standard CS formulation, let us define the sub-Nyquist measurements vector of the l-th observation as

$${\mathbf{y}}_l^{cs}={\left[y_l^{cs}\left[0\right],\cdots ,y_l^{cs}\left[M_r-1\right]\right]}^T,$$

(16)

collect all $N_a$ observations as ${\mathbf{Y}}^{cs}=\left[{\mathbf{y}}_1^{cs},{\mathbf{y}}_2^{cs},\cdots ,{\mathbf{y}}_{N_a}^{cs}\right]$ or its vector form ${\mathbf{y}}^{cs}=\mathrm{vec}\left({\mathbf{Y}}^{cs}\right)$, and define the sensing matrix $\mathbf{A}$ as

$${\left(\mathbf{A}\right)}_{m+(l-1)M_r,n}=A_{l,n}\left(m{T}_{cs}\right),$$

(17)

for $m\in \left[M_r\right]$, $l\in \left[N_a\right]$ and $n\in \left[N\right]$. Then, we have the sub-Nyquist measurements model ${\mathbf{y}}^{cs}=\mathbf{Ax}$ as given by (9) in the noisy case. However, the structure of the sensing matrix $\mathbf{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.

3.2. 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\left(t\right)$ can be expressed in Fourier series as

$$p_l\left(t\right)=\sum _{i=-L_p}^{L_p}{\rho }_{i,l}{e}^{j2\pi f_{p}it},$$

(18)

where $f_p=1/T$, ${\rho }_{i,l}$ is the Fourier coefficient of $p_l\left(t\right)$, and $L_p\ge BT$ is a positive integer. Denoting the discrete Fourier transform (DFT) of ${\mathbf{y}}_l^{cs}$ by ${\widehat{\mathbf{y}}}_l^{cs}$, it can be given by [29]

$${\widehat{\mathbf{y}}}_l^{cs}=\sqrt{\frac{1}{M_r}}{\mathbf{T}}_l{\mathbf{F}}_r{\mathbf{y}}_l,$$

(19)

where ${\mathbf{F}}_r$ is a $N_r\times N_r$ normalized DFT matrix after fftshift (fftshift is to shift zero-frequency component to the center of spectrum), and ${\mathbf{T}}_l$ is a $M_r\times N_r$ Toeplitz matrix defined by the Fourier coefficients as

$${\mathbf{T}}_l=\left[\begin{array}{cccc}{\rho }_{-L_0+N_r-1,l}& {\rho }_{-L_0+N_r-2,l}& \cdots & {\rho }_{-L_0,l}\\ {\rho }_{-L_0+N_r,l}& {\rho }_{-L_0+N_r-1,l}& \dots & {\rho }_{-L_0+1,l}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }_{L_0,l}& {\rho }_{L_0-1,l}& \cdots & {\rho }_{-L_0+M_r-1,l}\end{array}\right],$$

(20)

where $L_0=M_r+N_r-1$. The Toeplitz structure of ${\mathbf{T}}_l$ results from the convolution between the spectrums of the chipping sequence and the input signal. Due to bandpass filtering, ${\mathbf{T}}_l$ does not contain all Fourier coefficients of $p_l\left(t\right)$. Next, according to (19), we can describe the time-domain sub-Nyquist measurements as

$${\mathbf{y}}_l^{cs}={\mathbf{F}}_m^{*}{\widehat{\mathbf{y}}}_l^{cs}=\sqrt{\frac{1}{M_r}}{\mathbf{F}}_m^{*}{\mathbf{T}}_l{\mathbf{F}}_r{\mathbf{y}}_l,$$

(21)

where ${\mathbf{F}}_m$ is a $M_r\times M_r$ normalized DFT matrix after fftshift. Then, the measurement matrix ${\mathsf{\Phi }}_l$ can be described as

$${\mathsf{\Phi }}_l=\sqrt{\frac{1}{M_r}}{\mathbf{F}}_m^{*}{\mathbf{T}}_l{\mathbf{F}}_r,$$

(22)

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, $\mathbf{A}=diag\left({\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2,\cdots ,{\mathsf{\Phi }}_{N_a}\right)\mathbf{D}$. The first factor of $\mathbf{A}$ is a block diagonal matrix with independent diagonal blocks. Therefore, the sensing matrix $\mathbf{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.

3.3. 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 $\pm 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].

4. 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 $\mathbf{A}$ often precludes the direct application of off-the-shelf optimization methods because of two limitations: the matrix $\mathbf{A}$ cannot be stored explicitly, and the matrix-vector products involving $\mathbf{A}$ and its adjoint ${\mathbf{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 $\mathbf{A}$ and decompose the matrix-vector products into several fast operations. Then, the matrix $\mathbf{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.

4.1. Fast Matrix-Vector Products

As we have discussed, the sensing matrix $\mathbf{A}=\mathsf{\Phi }\mathbf{D}$ consists of two parts: measurement matrix $\mathsf{\Phi }$ and basis matrix $\mathbf{D}$. To obtain the desired fast matrix-vector products for $\mathbf{A}$ and ${\mathbf{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 ${\mathbf{T}}_l$ is a Toeplitz matrix, it can be viewed as a submatrix of a $\left(2L_0+1\right)\times \left(2L_0+1\right)$ circulant matrix ${\mathbf{C}}_l$ defined as

$${\mathbf{C}}_l=\left[\begin{array}{cccc}{\rho }_{L_0,l}& {\rho }_{L_0-1,l}& \cdots & {\rho }_{-L_0,l}\\ {\rho }_{-L_0,l}& {\rho }_{L_0,l}& \cdots & {\rho }_{-L_0+1,l}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }_{L_0-1,l}& \cdots & {\rho }_{-L_0,l}& {\rho }_{L_0,l}\end{array}\right].$$

(23)

Specifically, ${\mathbf{T}}_l$ is exactly the right $N_r$ columns and top $M_r$ rows of ${\mathbf{C}}_l$. Let us define row and column truncation matrices as

$${\Delta }_r=\left[\begin{array}{cc}{\mathbf{I}}_{M_r}& {\mathbf{0}}_{N_r-1}\end{array}\right],$$

(24)

and

$${\Delta }_a=\left[\begin{array}{c}{\mathbf{0}}_{M_r-1}\\ {\mathbf{I}}_{N_r}\end{array}\right],$$

(25)

respectively. Then, we embed ${\mathbf{T}}_l$ into ${\mathbf{C}}_l$ as follows:

$${\mathbf{T}}_l={\Delta }_r{\mathbf{C}}_l{\Delta }_a.$$

(26)

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 $\mathcal{O}\left(N\log \right(N\left)\right)$. The goal of (26) is to extend this property to the Toeplitz matrix ${\mathbf{T}}_l$. Specifically, let us define an integer $L=2L_0+1$, and a vector ${\mathbf{\rho }}_l={\left[{\rho }_{L_0,l},{\rho }_{-L_0,l},\cdots ,{\rho }_{L_0-1,l}\right]}^T$ which is the first column vector of ${\mathbf{C}}_l$, and denote the $L\times L$ DFT matrix as ${\mathbf{F}}_L$. Then, the circulant matrix ${\mathbf{C}}_l$ can be diagonalized as

$${\mathbf{C}}_l={\mathbf{F}}_L\mathrm{diag}\left(L{\mathbf{\theta }}_l\right){\mathbf{F}}_L^{*},$$

(27)

where ${\mathbf{\theta }}_l={\mathbf{F}}_L{\mathbf{\rho }}_l$. Inserting (26) and (27) into (21), we have

$${\mathbf{y}}_l^{cs}=\sqrt{\frac{L}{M_r}}{\mathbf{F}}_m^{*}{\Delta }_r{\mathbf{F}}_L\mathrm{diag}\left({\mathbf{\theta }}_l\right){\mathbf{F}}_L^{*}{\Delta }_a{\mathbf{F}}_r{\mathbf{y}}_l.$$

(28)

Defining $\mathsf{\Theta }=\left[{\mathbf{\theta }}_1,\cdots ,{\mathbf{\theta }}_{N_a}\right]$, then the 2D QuadCS measurement ${\mathbf{Y}}^{cs}$ can be written as

$${\mathbf{Y}}^{cs}=\sqrt{\frac{L}{M_r}}{\mathbf{F}}_m^{*}{\Delta }_r{\mathbf{F}}_L\left(\mathsf{\Theta }\circ {\mathbf{F}}_L^{*}{\Delta }_a{\mathbf{F}}_r\mathbf{Y}\right),$$

(29)

which mainly involves FFTs and thus has a low complexity. Note that, from (29), we find that the random matrix $\mathsf{\Theta }$ performs “Hadamard scrambling” to inject randomness into the measurements.

Next, we analyze the inner structure of the basis matrix $\mathbf{D}$. It has been suggested that conventional imaging algorithms, such as CSA, can be exploited to derive its fast matrix-vector product for $\mathbf{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

$${\mathbf{X}}_{CSA}=\left({\mathbf{H}}_3\circ {\mathbf{F}}_r^{*}\left({\mathbf{H}}_2\circ {\mathbf{F}}_r\left({\mathbf{H}}_1\circ \mathbf{Y}{\mathbf{F}}_a\right)\right)\right){\mathbf{F}}_a^{*}.$$

(30)

In this equation, ${\mathbf{X}}_{CSA}$ is the SAR image output by CSA, ${\mathbf{H}}_1$, ${\mathbf{H}}_2$, ${\mathbf{H}}_3$ are phase compensation matrices (i.e., their entries lie in the unit circle in $\mathbb{C}$, for more details, cf. chapter 6 in [35]), and ${\mathbf{F}}_r$, ${\mathbf{F}}_a$ are the range and azimuth DFT matrices, respectively. Because CSA is a high precision algorithm, ${\mathbf{X}}_{CSA}$ is a valid approximation of the reflectivity map $\mathbf{X}$. For ease of discussion, we omit the approximation error, i.e., we take $\mathbf{X}={\mathbf{X}}_{CSA}$, and in this sense, $\mathbf{X}$ and $\mathbf{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]

$$\mathbf{Y}=\left({\overline{\mathbf{H}}}_1\circ {\mathbf{F}}_r^{*}\left({\overline{\mathbf{H}}}_2\circ {\mathbf{F}}_r\left({\overline{\mathbf{H}}}_3\circ \mathbf{X}{\mathbf{F}}_a\right)\right)\right){\mathbf{F}}_a^{*}.$$

(31)

It is easy to see that the linear system (31) is equivalent to $\mathbf{y}=\mathbf{D}\mathbf{x}$ with proper reshape operations. Furthermore, since the involved FFTs and Hadamard products in (30) and (31) are unitary operators, the matrix $\mathbf{D}$ is also unitary [5].

4.2. 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 $\mathbf{A}$ and its adjoint ${\mathbf{A}}^{*}$ followed by a shrinkage/soft-threshold step as follows:

$${\mathbf{x}}_{k+1}={\mathcal{T}}_{\lambda /L_f}\left({\mathbf{z}}_k-{\mathbf{A}}^{*}\left(\mathbf{A}{\mathbf{z}}_k-{\mathbf{y}}_{cs}\right)/L_f\right),$$

(32)

where ${\mathbf{x}}_{k+1}$, ${\mathbf{z}}_k$ are iterative variables, $L_f$ is the Lipschitz constant for the gradient $\Delta f$ of the function $f\left(\mathbf{x}\right)={∥\mathbf{A}\mathbf{x}-{\mathbf{y}}^{cs}∥}_2^2$, and ${\mathcal{T}}_{\lambda /L_f}\left(·\right)$ is the soft shrinkage operator defined as

$${\mathcal{T}}_{\tau }{\left(\mathbf{x}\right)}_i=\left\{\begin{array}{cc}\left(\left|{\mathbf{x}}_i\right|-\tau \right)\frac{\mathbf{x}}{\left|{\mathbf{x}}_i\right|},\hfill & \left|{\mathbf{x}}_i\right|>\tau ;\\ 0,\hfill & \mathrm{otherwise}.\end{array}\right.$$

(33)

From (32), we know that the matrix-vector products involving $\mathbf{A}$ and its adjoint ${\mathbf{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.

Algorithm 1: Fast sparse imaging with QuadCS measurements

Input:  Operators $\mathcal{M}$, $\mathcal{D}$ and their adjoints; QuadCS measurement ${\mathbf{Y}}^{cs}$; number of iterations $Niter$; $\lambda$; Output:  Sparse SAR image ${\mathbf{X}}^{\#}$; 1: initial: take $\lambda$, $L_f\in \mathbb{R}$, ${\mathbf{Z}}_0={\mathbf{X}}_0=\mathbf{0}$, $q_0=1$, $k=1$; 2: repeat 3:  ${\mathbf{Z}}_k={\mathbf{Z}}_{k-1}-\frac{1}{L_f}{\mathcal{D}}^{*}{\mathcal{M}}^{*}\left(\mathcal{MD}{\mathbf{Z}}_{k-1}-{\mathbf{Y}}^{cs}\right)$; 4:  ${\mathbf{X}}_k={\mathcal{T}}_{\lambda /L_f}\left({\mathbf{Z}}_k\right)$; 5:  $q_k=\frac{1}{2}+\frac{1}{2}\sqrt{1+4q_{k-1}^2}$; 6:  ${\mathbf{Z}}_k={\mathbf{X}}_k+\frac{q_{k-1}-1}{q_k}({\mathbf{X}}_k-{\mathbf{X}}_{k-1})$; 7:  $k=k+1$; 8: until$k\ge Niter$; 9: return ${\mathbf{X}}^{\#}={\mathbf{X}}_{k-1}$;

For ease of description, we denote the right sides of (29) and (31) by operators as $\mathcal{M}\mathbf{Y}$ and $\mathcal{D}\mathbf{X}$, respectively, and denote $\mathrm{vec}(·)$ by an operator $\mathcal{V}$. With these notations, we can express the matrix-vector product of $\mathbf{A}$ as

$$\mathbf{Ax}=\mathcal{V}\mathcal{MD}{\mathcal{V}}^{*}\mathbf{x}.$$

(34)

Similarly, the matrix-vector product of the adjoint ${\mathbf{A}}^{*}$ is given by

$${\mathbf{A}}^{*}{\mathbf{y}}^{cs}=\mathcal{V}{\mathcal{M}}^{*}{\mathcal{D}}^{*}{\mathcal{V}}^{*}{\mathbf{y}}^{cs},$$

(35)

where the adjoint operator ${\mathcal{D}}^{*}$ is given by the operations on the right side of (30) and ${\mathcal{M}}^{*}$ is given by

$${\mathcal{M}}^{*}{\mathbf{Y}}^{cs}=\sqrt{\frac{L}{M_r}}{{\mathbf{F}}_r}^{*}{\Delta }_a^T{\mathbf{F}}_L\left(\overline{\mathsf{\Theta }}\circ {\mathbf{F}}_L^{*}{\Delta }_r^T{\mathbf{F}}_m{\mathbf{Y}}^{cs}\right).$$

(36)

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 $\mathbf{A}$ and its adjoint ${\mathbf{A}}^{*}$ both admit fast matrix-vector products.

4.3. Complexity

Without fast matrix-vector products for $\mathbf{A}$ and ${\mathbf{A}}^{*}$, we need to store the $M\times N$ matrix $\mathbf{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 $\mathbf{A}$ with double precision, and roughly $1.3744\times {10}^{11}$ complex multiplications to compute the matrix-vector product $\mathbf{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\times N_a$ matrices ${\mathbf{H}}_1$, ${\mathbf{H}}_2$, ${\mathbf{H}}_3$ and $\mathsf{\Theta }$. It only requires 32 MB to store the four matrices and roughly $8.2\times {10}^4$ complex multiplications to compute $\mathbf{Ax}$ for a scene of $1024\times 1024$ size, which are significantly smaller than the previous requirements. In general, the storage cost is reduced form $\mathcal{O}\left(MN\right)$ to $\mathcal{O}\left(N\right)$, and the computational cost is reduced from $\mathcal{O}\left(MN\right)$ to $\mathcal{O}\left(N\log \right(N\left)\right)$ 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\left(log\left(N_r\right)+log\left(N_a\right)\right)+3N_r{N}_a$ complex products, or in a compact form, $\mathcal{O}\left(Nlog\left(N\right)\right)$ 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.

5. Restricted Isometry Property (RIP) Analysis

In this section, we derive the RIP of the sensing matrix $\mathbf{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 $\mathbf{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 $\mathbf{A}$ because of its special factornizatio. In the following, we first derive a CoM inequality of the matrix $\mathbf{A}$, and then establish a RIP condition on it. The CoM inequality of $\mathbf{A}$ is formulated in Lemma 1.

Lemma 1.

Let $\mathbf{x}$ be a s-sparse vector in ${\mathbb{C}}^N$. For an arbitrary $0<\delta <1$, we have

$$\mathbb{P}\left(\left|{∥\mathbf{A}\mathbf{x}∥}_2^2-{∥\mathbf{x}∥}_2^2\right|\ge \delta {∥\mathbf{x}∥}_2^2\right)\le 4\mathit{exp}\left(-\frac{\delta \xi {\varpi }^{2}M}{16\chi s}\right),$$

(37)

where M is the number of sub-Nyquist measurements, χ denotes two settings of chipping sequences and is defined as

$$\chi =\left\{\begin{array}{cc}1,\hfill & \mathit{independent}\mathit{chipping}\mathit{sequences},\\ N_a,\hfill & \mathit{equal}\mathit{chipping}\mathit{sequences}.\end{array}\right.$$

(38)

ξ is the duty cycle of radar transmitted signals, and ϖ is the maximum value of the normalized antenna beam pattern defined as

$$\varpi =\frac{w_a\left(0\right)T_a}{\sqrt{{\int }_{-T_a/\phantom{T_{a}2}2}^{T_a/\phantom{T_{a}2}2}{w}_a^2\left(\eta \right)d\eta }},$$

(39)

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 $\mathbf{A}$, which is the essence of RIP.

Based on the CoM inequality in Lemma 1, we can derive a uniform tail bound of ${∥{\mathbf{A}}_S^{*}{\mathbf{A}}_S-{\mathbf{I}}_s∥}_2$ for a fixed subset $S\subset \left[N\right]$ with $card\left(S\right)=s$ using the covering number estimate [38]. Then, we take the union bound over all subsets $S\subset \left[N\right]$ of cardinality s, to derive the following RIP condition.

Theorem 1.

The restricted isometry constant ${\delta }_s$ of the sensing matrix $\mathbf{A}$ satisfies ${\delta }_s<\delta$ for some constant $0<\delta <1$ with probability exceeding $1-\eta$ provided that

$$\begin{array}{cc}\hfill M\ge & \frac{32\chi s}{\delta \xi {\varpi }^2}\left(sln\left(\frac{eN}{s}\right)+4sln\left(3\right)+\ln \left(\frac{\eta }{4}\right)\right).\hfill \end{array}$$

(40)

Proof. 

The proof is given in Appendix A.2.□

Theorem 1 gives the RIP condition on the sensing matrix $\mathbf{A}$. The key parameters in (40) are M and $\chi$, 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 $\mathbf{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 $\xi$ and $\varpi$ 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 $\xi$ or $\varpi$, each column of the basis $\mathbf{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 ${\delta }_s$ and ${\gamma }_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 ${\delta }_s<\delta$ and ${\gamma }_s<\delta /N_a$ simultaneously for some constant $0<\delta <1$ with probability exceeding $1-\eta$.

Replacing ${\delta }_s$ by ${\delta }_s/\chi$ 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-\eta$, let $s_1$ and $s_2$ be the maximum integers that satisfy ${\delta }_{2s_1}<{\delta }^{\#}$ and ${\gamma }_{2s_2}<{\delta }^{\#}$ with equal and independent chipping sequences, respectively, where ${\delta }^{\#}=\sqrt{2}-1$. Let ${\mathbf{x}}_1^{\#}$ and ${\mathbf{x}}_2^{\#}$ be the solutions of the constrained BPDN problem

$$\underset{\mathbf{x}}{min}{∥\mathbf{x}∥}_1\mathit{subject}\mathit{to}∥\mathbf{A}\mathbf{x}-{\mathbf{y}}^{cs}∥\le ϵ$$

(41)

in the two settings of chipping sequences (38), respectively, then

$$\begin{array}{cc}\hfill {∥\mathbf{x}-{\mathbf{x}}_2^{\#}∥}_2& \le C_1\frac{{\sigma }_{s_2}{\left(\mathbf{x}\right)}_1}{\sqrt{s_2}}+C_2ϵ\hfill \\ & \approx C_1\frac{{\sigma }_{\sqrt{N_a}{s}_1}{\left(\mathbf{x}\right)}_1}{\sqrt{\sqrt{N_a}{s}_1}}+C_2ϵ\hfill \\ & <C_1\frac{{\sigma }_{s_1}\left(\mathbf{x}\right)}{\sqrt{s_1}}+C_2ϵ,\hfill \end{array}$$

(42)

where the last line is exactly the error bound of ${\parallel \mathbf{x}-{\mathbf{x}}_1^{\#}\parallel }_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 (42), 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\approx \sqrt{N_a}$, which results in much lower best s-term approximation error as shown in (42). 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 ${\mathbf{x}}_2^{\#}$ always has a smaller recovery error than ${\mathbf{x}}_1^{\#}$. Note that we recover images by solving (10) instead of (41). In fact, the two problems are equivalent. If ${\mathbf{x}}^{\#}$ is a minimizer of the (41) with some $ϵ>0$, there necessarily exists $\lambda \ge 0$ such that ${\mathbf{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 $\lambda$ is related to ${\mathbf{x}}^{\#}$ and $ϵ$. In this paper, we assume that the noise level is known for (41), and $\lambda$ 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.

6. 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, ${\parallel \mathbf{x}-{\mathbf{x}}^{\#}\parallel }_2/{\parallel \mathbf{x}\parallel }_2$, and the relative root-mean-square error (RRMSE), $\mathbb{E}\left\{{\parallel \mathbf{x}-{\mathbf{x}}^{\#}\parallel }_2/{\parallel \mathbf{x}\parallel }_2\right\}$, are used as the imaging performance metrics.

The simulation scheme is shown in Figure 3 and main radar parameters are summarized in

Table 1. Simulation parameters.

Parameter	Value
Carrier frequency	5.3 GHz
Fast-time frequency bandwidth	30.11 MHz
Pulse width	41.74 $\mathsf{\mu }$s
Pulse repetition frequency	1256.98 Hz
Effective radar velocity	7062 m/s
Slant range of scene center	150.1 km

. 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

$$x_i^{gray}=\left\{\begin{array}{cc}|x_i|/\tau ,\hfill & |x_i|\le \tau ;\hfill \\ 1,\hfill & \mathrm{else},\hfill \end{array}\right.$$

(43)

where $\tau$ 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 $\alpha$.

6.1. 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 (42); then, the improved performance is gradually diluted.

6.2. 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\times 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\times 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\times 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.

7. Conclusions

In this paper, we presented an efficient sampling system for sub-Nyquist acquisition of SAR data. The proposed system is based on the QuadCS, which is easy to insert randomness into the sensing matrix. With independent measurements in each observation, our system has more randomness than that with fixed measurements in all observations, and thus the corresponding sensing matrix has better image recovery performance. Based on the RIP analysis, we showed that the proposed system has smaller RIC, which can lead to lower sampling/data rate and smaller recover error.

Next, we provided a fast image formation method. Through the frequency-domain representation, we found that the sensing matrix is a block and Toeplitz matrix and its matrix-vector products can be decomposed into several fast operators. In virtue of matrix decomposition, we established a first-order fast SAR imaging algorithm, whose storage and computational complexity are much less than that of the off-the-shelf recovery methods. By using both synthetic data and real data, we made comparisons among QuadCS with independent measurements, QuadCS with fixed measurements, and Xampling. Numerical results showed that the QuadCS with independent measurements achieves much better image qualities than the other simulated schemes.