1. Introduction

Synthetic Aperture Radar (SAR) is widely used to obtain images under all weather conditions at any time of a day^[1]. With the growing demand for SAR images and the development of SAR imaging technology, the image resolution has increased significantly. Traditionally, the scene is reconstructed by Matched Filter (MF)-based focusing algorithms. However, these methods are efficient above Nyquist sample rate, which would bring a lot of burden to current SAR system hardware with too much data.

The recent theory of Compressed Sensing (CS) appears to be a powerful method to solve the problem above. CS theory proves that a sparse signal can be reconstructed with much fewer measurements than Nyquist sample rate needs^[2-5]. In the past several years, CS has been used in many radar applications and shown very good results. The theory was applied to eliminate the need of the pulse compression matched filter and reduce the required receiver's analog-to-digital conversion bandwidth^{[6, 7]}. In Ref. ^[8], after range compression, CS was then used on azimuth, but the signal was operated in One-Dimensional (1-D) form. For Two-Dimensional (2-D) operations, some work has been done in Refs. [9-12]. Many different algorithms have been used to solve SAR imaging problem, including matrix completion^[13], 1-bit compressive sensing^[14], truncated Singular Value Decomposition (SVD)^[15], and Bayesian compressive sensing^[16].

Traditional MF-based methods have relatively lower computational complexity and memory cost compared to CS-SAR methods, but they suffer from high sidelobes. Range Doppler Algorithm (RDA), a traditional MF-based method of SAR imaging, has some drawback in correcting the coupling between range and azimuth. Another traditional MF-based method, the Chirp Scaling Algorithm (CSA), uses an approximated equation of signal when operating range IFFT. The omega-K algorithm adopts a special operation to correct the coupling of range and azimuth to overcome the drawback of RDA. Besides, omega-K algorithm uses a precise expression of signal compared to CSA^[1].

In this study, we propose a new CS-SAR method for sparse scene based on omega-K algorithm, which is named as Compressed Sensing omega-K Algorithm (CS-OKA). This new method could reduce computational complexity and memory cost while achieving high resolution and eliminating sidelobes with reduced sampling rate. First, we introduce the procedure of omega-K algorithm and represent it in the form of multiplications of matrices. Then we try to derive its inverse in order to get raw signal from SAR image. Similar practice has been explored before finding the inverse of MF-based method^[17]. Finally, we combine the inverse of omega-K algorithm and CS framework to formulate a method from which undersampled data can recover the original SAR image fast and precisely. Different from previous CS-SAR algorithm, we get the echoes from the inverse of omega-K algorithm, not from the original observation. This results in the reduction of computational complexity. The whole procedure is performed two-dimensionally. To get the expected sparse scene, we use Iterative Thresholding Algorithm^{[18, 19]} (ITA) to solve the sparse problem. The CS-OKA is much faster than previous CS-SAR algorithms since we substitute the large matrix-vector multiplication with omega-K algorithm and its inverse which can be realized efficiently. Also, memory cost is significantly reduced. Simulations show that CS-OKA can reconstruct the sparse SAR images effectively below Nyquist rate. Refs. [10] and [11] use similar process but with RDA and CSA, respectively. As mentioned above, omega-K has some advantages compared to RDA and CSA, and that is why we chose omega-K in our algorithm.

This paper is organized as follows: In Section 2, the CS theory is introduced. In Section 3, the inverse of omega-K algorithm is described. In Section 4, our new CS-SAR method CS-OKA is presented. In Section 5, simulation results are shown to evaluate the performance of CS-OKA. Finally, Section 6 concludes the paper.

2. CS Theory

Suppose x$\in {{\mathbb{C}}^{N}}$ is the signal to be estimated. We call it K-sparse if its nonzero elements are no more than K^{[2, 5]}. Let $x\in {{\mathbb{C}}^{M}}$ be the collected data, and matrix $\mathit{\Theta }\in {{\mathbb{R}}^{M\times N}}(M < N)$ be the observation matrix. According to CS theory, it is possible to recover x if x is sparse:

where n $\in {{\mathbb{C}}^{M}}$ represents the additive noise. For those signals which are not sparse themselves, we sometimes can find a transform or basis so that their coefficients under the basis are sparse:

f is the non-sparse signal to be estimated; Ψ is the basis; x and is the sparse coefficient vector under the basis. In this case, the expression would be

where A is sensing matrix.

To recover the sparse signal x, CS theory puts forward some requirements on the sensing matrix. One sufficient condition is the Restricted Isometry Property (RIP). If there exists a constant ${{\delta }_{K}}\in \left( 0, 1 \right)$ such that

holds for all K-sparse signal x, the matrix A is said to satisfy the RIP of order K. When the RIP constant satisfies some conditions, the recovery can be achieved by solving a convex optimization problem:

where ‖·‖₁ denotes the l₁ norm and e bounds the energy of noise (${{\left\| \mathit{\boldsymbol{n}} \right\|}_{2}}\le \varepsilon$). An equivalent regularization scheme has the form of

where λ is a regularization parameter. ITA is an general algorithm to solve the optimization problem. It uses the iteration

to implement the procedure. In Eq. (7), μ is the step size that controls the convergence. Let σ=λμ, then E_{q, σ} is the thresholding operator which is defined in terms of components by

where e_{q, σ} can be analytically specified when q=0, 1/2, 2/3, 1. If q=1, then it corresponds to soft-thresholding, and E_{q, σ} would have the following expression:

Solving Eq. (5) or Eq. (6) with ITA is very time consuming, especially for large scale SAR imaging. Because it involves large matrix-vector multiplications. In Section 3, we will derive a new way to form the sensing matrix, which would accelerate the process significantly.

3. Inverse Omega-K Algorithm

As mentioned in Section 2, when using Eq. (7) to solve the optimization problem, we have to conduct multiplications of very large scale matrices. However, if we can substitute A and A^H with much simpler operators, we will avoid such time consuming multiplications. In this study, we derive such operators based on omega-K algorithm. Detailed information is further discussed later.

The omega-K algorithm^[1] has a sufficient accuracy with a moderate computational cost. In this section, after a brief review of omega-K algorithm, we will derive the inverse of omega-K algorithm.

3.1 Omega-K algorithm

Let s₀(τ, η) be the received data, where τ and η represent fast range time and slow azimuth time, respectively:

where A₀ is a complex constant; w_r(τ)=rect (τ/T_r) and T_r is the transmitted signal duration of each pulse; w_a(η)=rect (η/Tη) and T_η is the integration time; R(η) is the slant range; f₀ is the carrier frequency; K_r is chirp rate. To perform omega-K algorithm, we need to go through four steps.

The first step is 2-D Fourier transform. Let ${{\cal F}_\tau }\left\{ \cdot \right\}$ and ${{\cal F}_\eta }\left\{ \cdot \right\}$ be the Fourier transform operator for range and azimuth, respectively. The transform would be:

The second step is multiplying Eq. (11) with a reference function ${H_{{\rm{ref}}}}\left( {{f_\tau }, {f_\eta }} \right)$. This step is to compensate the phase for the target in the scene center:

where R_ref is the slant range of reference target, and ${V_{{{\rm{r}}_{{\rm{ref}}}}}}$ is radar velocity at reference target position. We usually assume that radar velocity remains the same (V_r) during the procedure. After reference function multiplication, Eq. (11) becomes

The phase of ${S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right)$ is

The third step is Stolt interpolation. If the square root term in Eq. (14) is substituted by

then the coupling effect between range and azimuth is thus removed. After the substitution, the frequency spectrum of ${S_{{\rm{RFM}}}}\left( {{f_\tau }, {f_\eta }} \right)$ is changed. For the convenience of the last step, a truncated sinc-kernel interpolation is carried out. We use ${\rm{C}}\left\{ \cdot \right\}$ to represent the Stolt interpolation operator, then Eq. (13) becomes

The last step is 2-D inverse Fourier transform to get the final image $g\left( {\tau, \eta } \right)$:

We can sample the continuous-time analog signals, ${s_0}\left( {\tau, \eta } \right)$, $g\left( {\tau, \eta } \right)$, into two-dimensional arrays, S $\in {{\mathbb{C}}^{{{N}_{\tau }}\times {{N}_{\eta }}}}$, G $\in {{\mathbb{C}}^{{{N}_{\tau }}\times {{N}_{\eta }}}}$, then the omega-K algorithm can be written as

where F_τ and F_η represent Fourier transform matrices in range and azimuth, respectively; ${\left\{ \cdot \right\}^{\rm{H}}}$ denotes the conjugate transposition.

3.2 Inverse omega-K algorithm

At the beginning of this section, we have introduced that if there exists a substitution for A, then the efficiency of the algorithm can be improved. Inspired by this idea, we derive the inverse of omega-K algorithm. Since omega-K algorithm is used to get SAR image from echoes, then it is intuitive that we can get echoes from SAR image using the inverse of omega-K algorithm. From Eq. (3), we are able to say that the inverse of omega-K algorithm has the same effect as A, thus can substitute A in Eq. (7).

We can easily conclude that omega-K is a linear algorithm that only contains Fourier transform, complex multiplication, and Stolt interpolation. The inverse of these three operations are also linear and the general procedures are quite obvious:

where D represents the inverse of Stolt interpolation. It takes a similar form as Eq. (15):

We substitute the square root with ${f_0} + {f_\tau }$. The inverse of truncated sinc-kernel interpolation is hard to achieve directly. In our method, we also use a truncated sinc-kernel interpolation to approximate its inverse.

Note that matrices in Eq. (19) could be very large, but we can perform Fourier transform and inverse Fourier transform directly instead of con-ducting matrices multiplications. It can be a signi-ficant save of time consumption. Meanwhile, we do not have to restore the large matrices either, which is also a save of memory cost.

The inverse of omega-K algorithm can be concluded accordingly: Fourier transform in range and azimuth first; then the inverse of Stolt interpolation; followed by is multiplication with the conjugate transposition of the reference function; at last IFFT in range and azimuth.

We use T to represent the procedure of inverse omega-K algorithm above, then Eq. (19) becomes

Our goal is to get G from S using T(·)·.

4. Compressed Sensing Inverse Omega-K Algorithm

We can conclude from Eq. (21) that it is possible to get echoes from the scene directly and precisely if no noise exists. If there is noise in the signal, we can draw the following objective function for sparse scene with the help of CS theory:

We know that there is a lot of redundancy in SAR signals and that CS theory has proved to be capable of reconstructing a sparse signal with much fewer measurements than Nyquist theory needs. So we sample the signals to reduce redundancy:

where ${{\mathit{\Theta }}_{\tau }}\in {{\mathbb{C}}^{N_{\tau }^{'}\times {{N}_{\tau }}}}, \left( N_{\tau }^{'} < {{N}_{\tau }} \right)$ and ${{\mathit{\Theta }}_{\eta }}\in {{\mathbb{C}}^{{{N}_{\eta }}\times N_{\eta }^{'}}}, \left( N_{\eta }^{'}\times {{N}_{\eta }} \right)$ represent random sampling operators in range and azimuth, respectively. Then the objective function becomes:

Since T is linear, Eq. (24) is a convex problem that can be solved with many convex optimization algorithms. In this study, we choose ITA to get G. In Section 2, we have introduced that Eq. (7) can be a solution to Eq. (6). Then in Section 3, we derived the inverse of omega-K, T, to substitute A. In order to solve the objective function Eq. (24) effectively without multiplications of large matrices, we still need to find a simpler operator to substitute A^H. Since A^H can be viewed as the inverse of A, it is clear that omega-K algorithm can be a substitution for A^H. We use M to represent omega-K algorithm, then from Eq. (7) we can get the following solution to Eq. (24):

M has a similar form as T, so it can save time and memory cost too. From Eq. (25), we can give the following algorithm procedure:

Algorithm 1: Iterative thresholding algorithm for proposed CS-SAR imaging
Input: SAR raw echoes SS,
omega-K algorithm M and inverse omega-K algorithm T,
sampling operator Θτ and Θη
Initial: G(0), λ, μ, and max iteration Imax
1: for i = 1 to Imax do
2: residue: ${{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}={{\mathit{\boldsymbol{S}}}_{S}}-{\mathit{\Theta }_\tau }T\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right){{\mathit{\Theta }}_{\eta }}$
3: omega-K on residue: $\Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}=M\left( \mathit{\Theta } _{\tau }^{\rm{T}}{{\mathit{\boldsymbol{R}}}^{\left( i-1 \right)}}\mathit{\Theta } _{\eta }^{\rm{T}} \right)$
4: Thresholding: ${{\mathit{\boldsymbol{G}}}^{\left( i \right)}}={{E}_{1, \lambda \mu }}\left( {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}}+\mu \Delta {{\mathit{\boldsymbol{G}}}^{\left( i-1 \right)}} \right)$
5: end for
Output: The recovery image G*=G(i)

 | Show Table

DownLoad: CSV

Let's generally compare the computational complexity of CS-OKA and other CS-SAR methods which also use ITA to solve the problem. First, we denote $N = {N_\tau } \times {N_\eta }$. In each iteration, the main computational complexity of CS-OKA results from omega-K, inverse omega-K, which have commonly the computational complexity of ${\mathcal O}\left( {N{\rm{log}}N} \right)$, and the thresholding operator with complexity of ${\mathcal O}\left( N \right)$. The total complexity is ${\mathcal O}\left( {{I_{\max }}N{\rm{log}}N} \right)$. For other CS-SAR methods that use Eq. (7) to solve the problem, multiplication of two large matrices brings the complexity of ${\mathcal O}\left( {NsU} \right)$ in each iteration, where represents the sampling rate, and U is the product of samples of sent pulses and samples of azimuth time. Then the total complexity of other CS-SAR methods is ${\mathcal O}\left( {{I_{\max }}NsU} \right)$. In practice, U is usually very large to acquire high Signal-to-Noise Ratio (SNR). As a result, ${\mathcal O}\left( {{I_{\max }}NsU} \right)$ is much higher than ${\mathcal O}\left( {{I_{\max }}N{\rm{log}}N} \right)$. CS-OKA obviously has lower computational complexity.

As for memory cost, since other CS-SAR methods that use ITA have to restore the big matrices A and A^H, the memory cost would be ${\mathcal O}\left( {{N^2}} \right)$. For CS-OKA, we only need to restore the reference function H_ref. That leads to the memory cost of ${\mathcal O}\left( N \right)$. Apparently, CS-OKA has much lower memory cost than other CS-SAR methods.

CS-OKA uses the redundancy of SAR signal to reduce the data rate compared with traditional MF-based methods. Furthermore, sparse optimization method is adopted to explore the sparse property of the signal, which could significantly suppress the sideslobes and improve the resolution. Also, it applies the inverse of omega-K algorithm to substitute the large matrix operations in other CS-SAR algorithms, which significantly contributes to the improvement of time complexity.

5. Simulation Results

In this section, we will carry out several simulations to demonstrate the effectiveness of our method. The simulations were taken under different sampling rates and SNRs. In our simulations, signal in SNR refers to the echo. All simulations are performed under side-looking SAR. Tab. 1 lists the parameters of our simulations.

Table  1.  Parameters used in the simulation

Parameter	Value
Pulse duration (μs)	1.33
Bandwidth in range (MHz)	150
Carrier frequency (MHz)	600
Sampling rate (MHz)	225
Slant range of scene center (m)	1200
Length of synthetic aperture (m)	300
Pulse repeat frequency (Hz)	150 000 000
Radar velocity in azimuth (m/s)	150

 | Show Table

DownLoad: CSV

In our simulations, we put 4 targets in the scene, all with the same unit amplitude. We first recover the scene using omega-K algorithm and CS-OKA with full sampled data. As shown in Fig. 1, both algorithms can recover the scene under different SNRs from 20 dB to 5 dB. Fig. 2 shows the contours of magnitude of the lower left target in Fig. 1. It is clear that results of CS-OKA has much lower sidelobes than omega-K algorithm. Then we undersample the signal (echo) and try to recover the scene with CS-OKA. In Fig. 3, we randomly sample 50% data in range and azimuth, respectively (i.e., for ${\mathit{\Theta }_{\tau }}\in {{\mathbb{C}}^{N_{\tau }^{'}\times {{N}_{\tau }}}}, N_{\tau }^{'}=0.5{{N}_{\tau }};\text{for}\ {\mathit{\Theta }_{\eta }}\in {{\mathbb{C}}^{{{N}_{\eta }}\times N_{\eta }^{'}}}, N_{\eta }^{'}=0.5{{N}_{\eta }}$). The results proved the validity of our method. We further undersample the signal to 10% data in both range and azimuth. Results are shown in Fig. 4. The scene can still be recovered perfectly even with 1% sampled data. In Fig. 5, we put two targets that are very close to each other. omega-K algorithm failed to separate the two targets because of high sidelobes [Figs. 5(a) and 5(b)], while CS-OKA can clearly distinguish the two targets [Figs. 5(c) and 5(d)].

Figure  1.  Reconstruction results with full sampled data but different SNRs. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

DownLoad: Full-Size Img PowerPoint

Figure  2.  Contours of magnitude of lower left target in Fig. 1. The left column is recovered by omega-K algorithm, and the right column is recovered by CS-OKA. From top to bottom, SNRs are 20, 10, and 5 dB, respectively

DownLoad: Full-Size Img PowerPoint

Figure  3.  Reconstruction results with 50% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

DownLoad: Full-Size Img PowerPoint

Figure  4.  Reconstruction results with 10% samples in both range and azimuth. From left to right, SNRs are 20, 10, and 5 dB, respectively

DownLoad: Full-Size Img PowerPoint

Figure  5.  Reconstruction results of two close targets

DownLoad: Full-Size Img PowerPoint

It is obvious that CS-OKA works very well even when the sample rate is low. Besides, CS-OKA maintains a high quality of the recovered images under different SNRs. It proves that our method is robust against noise because of the intrinsic property of Eq. (6). Also, CS-OKA has a higher resolution compared to omega-K algorithm.

6. Conclusion

In this study, we have proposed a new method for SAR imaging called CS-OKA. Different from previous CS-SAR methods, we combine the omega-K algorithm and CS to recover the scene.

We first derived the inverse of omega-K algorithm to get echoes from the scene. It can be a substitution for a large matrix in the original ITA. The inverse of omega-K algorithm can be realized efficiently with the help of several linear operations, which is also a save for memory cost. Then we proved that omega-K algorithm can substitute the other large matrix for convenience. It has a similar form as its inverse. These two substitutions together have significantly improved the computational complexity since the multiplications of large matrices are avoided. We have analyzed the detailed computational complexity and proved the efficiency of CS-OKA.

Simulations have shown that our method per-formed well under low sampling rate and is robust under different SNRs when the scene is sparse.