Sub-Nyquist SAR Imaging and Error Correction Via an Optimization-Based Algorithm

Abstract

Sub-Nyquist synthetic aperture radar (SAR) based on pseudo-random time–space modulation has been proposed to increase the swath width while preserving the azimuthal resolution. Due to the sub-Nyquist sampling, the scene can be recovered by an optimization-based algorithm. However, these methods suffer from some issues, e.g., manually tuning difficulty and the pre-definition of optimization parameters, and a low signal–noise ratio (SNR) resistance. To address these issues, a reweighted optimization algorithm, named pseudo-$ℒ$₀-norm optimization algorithm, is proposed for the sub-Nyquist SAR system in this paper. A modified regularization model is first built by applying the scene prior information to nearly acquire the number of nonzero elements based on Bayesian estimation, and then this model is solved by the Cauchy–Newton method. Additionally, an error correction method combined with our proposed pseudo-$ℒ$₀-norm optimization algorithm is also present to eliminate defocusing in the motion-induced model. Finally, experiments with simulated signals and strip-map TerraSAR-X images are carried out to demonstrate the effectiveness and superiority of our proposed algorithm.

1. Introduction

High-resolution wide-swath (HRWS) synthetic aperture radar (SAR) provides a short repeat cycle so that it has a high-efficiency acquisition capability [1,2]. Although the azimuthal multi-channel SAR [3] and multi-input multi-output (MIMO) SAR [4] can achieve a high azimuthal resolution and wide-range swath, it has a large amount of data and a large antenna. As the compressive sensing (CS) theorem posits [5,6,7], an innovative system concept called sub-Nyquist SAR based on pseudo-random time–space modulation breaks the limitation of the Nyquist sampling theorem with a single channel [8]. Simultaneously, the spatial and temporal phase modulation based on the optimization-based algorithm guarantees azimuthal resolution and wide-swath coverage mosaicked by several range sub-swaths. It adopts sub-Nyquist sampling along the azimuthal dimension, and the observed scene is recovered by sub-Nyquist SAR imaging [8].

In spite of the sub-Nyquist SAR achieving the above merits, sub-Nyquist SAR imaging still has two limitations, as follows:

• Under the assumption of satisfying restricted isometry property (RIP) [9], CS algorithms are vital to sub-Nyquist SAR imaging and include greedy algorithms [10,11], the $ℒ$₁-norm optimization algorithms [12,13], and Bayesian-based methods [14,15,16], where the $ℒ$₁-norm optimization algorithm has a better performance in terms of the recovered error evaluated by the mean square error (MSE) [17,18,19]. Although the $ℒ$₁-norm optimization algorithm has achieved a better-recovered performance, these methods still suffer from some problems, i.e., manually tuning difficulty and the pre-definition of optimization parameters (e.g., regularization parameter and thresholding parameter), and a low signal–noise ratio (SNR) resistance. However, these CS methods do not take full advantage of the scene prior information that we may hold, and sparse property is imposed uniformly and independently on each variable. Some low signal–noise ratio (SNR) targets in the sparse scene cannot be accurately recovered and it often yields false targets by the $ℒ$₁-based method. Although some reweighted optimization-based algorithm has already been proposed [20,21], there is still no knowledge of how and why to select an approximately fair rule in sub-Nyquist SAR imaging to further mitigate the impact of empirical parameter setting on reconstructed performance.
• In addition, the imaging process requires the knowledge of motion parameters, e.g., radar position and radar equivalent velocity [5]. However, the radar platform may deviate from the pre-defined track, and the equivalent velocity is estimated by the curve-fitting method or approximate expression in a practical application [22]; uncertainties and errors may be introduced into the motion-induced model so that the recovered scene may defocus to decrease the image quality [23,24,25,26]. A technology called auto-focusing removes these phase errors [27]. In recent years, many sparsity-driven algorithms [28,29,30,31,32,33,34] have been proposed to solve the defocusing problem and achieve an effective performance. However, the references [28,29,30,31,32,33] do not fully formulate the motion error and adopt an approximate expression so that the error is not removed. The reference [34] integrated the deep SAR imaging algorithm to remove the motion error.

In this paper, we propose a pseudo-$ℒ$₀-norm optimization algorithm based on Bayesian estimation to further improve the sub-Nyquist SAR imaging performance. The proposed algorithm penalizes the regularization item with the scene prior information, i.e., the reciprocal of its previous solution, to nearly acquire the number of nonzero values. Since this method approximates a $ℒ$₀-based algorithm, which needs to be minimized in the sparse recovery, we named it a pseudo-$ℒ$₀-norm optimization algorithm. Sub-Nyquist SAR imaging includes three steps: range compression, range cell migration correction (RCMC), and azimuth compression [5]. The traditional matched filtering (MF) method is adopted for range compression. After RCMC, a pseudo-$ℒ$₀-norm optimization algorithm is adapted to achieve azimuth compression. This method not only takes advantage of the scene prior information but also establishes an approximately fair penalized rule so that it can recover low SNR targets and remove false targets compared to the $ℒ$₁-norm optimization algorithm. In addition, an error correction method integrated with a pseudo-$ℒ$₀-norm optimization algorithm eliminates the influence of phase error and improves the image quality. This method includes two steps: scene recovery based on a pseudo-$ℒ$₀-norm optimization algorithm and error estimation by minimizing the least-squares target function. The two steps are successively iterative and the recovered matrix is updated according to the estimated error. For the exact expression of error is not fully formulated, we use the random phase error regardless of the error expression. The numerical simulation results are detailed in the following sections to make an evident advantage of the proposed algorithm.

The rest of this paper is organized as follows. In Section 2, the observation model for the sub-Nyquist SAR based on the pseudo-random time–space modulation is first described. Then, we build a pseudo-$ℒ$₀-norm regularization model based on Bayesian estimation, and this regularization model can be solved by the Cauchy–Newton method. Additionally, an error correction method integrated with a pseudo-$ℒ$₀-norm optimization algorithm is proposed to remove errors and eliminate defocusing. Simulation experiments and data experiments with real TerraSAR-X images confirm the effectiveness and superiority of our proposed method, presented in Section 3. In Section 4, the discussion shows the performance and advantages of our proposed algorithm. Section 5 concludes this paper.

2. Materials and Methods

In this section, the observation model of the sub-Nyquist SAR based on the pseudo-random time–space modulation is described. To further improve the sub-Nyquist SAR imaging performance, a pseudo-$ℒ$₀-norm regularization model is built and solved by the Cauchy–Newton method to obtain a pseudo-$ℒ$₀-norm optimization algorithm. In addition, an error correction method integrated with our proposed pseudo-$ℒ$₀-norm optimization algorithm is proposed to eliminate the effect of the phase error and achieve autofocusing.

2.1. Sub-Nyquist SAR Imaging and Error Correction Signal Models Based on Pseudo-Random Time–Space Modulation

2.1.1. Sub-Nyquist SAR Imaging Model

For traditional HRWS systems, e.g., the azimuthal multi-channel SAR [3] and MIMO SAR [4], the equivalent sampling still satisfies the Nyquist theorem. To lower the amount of data and break the conflict between high resolution and wide swath, the sub-Nyquist SAR based on pseudo-random time–space modulation has been proposed [8]. The imaging geometry of the sub-Nyquist SAR is shown in Figure 1, where $\eta$ is slow time along the azimuth and $R_i\left(\eta \right)$ represents the range between the radar and the point target located at the coordinate $\left(x_i,y_i,0\right)$ at the azimuth time $\eta$. $x_i$ and $y_i$ denote the azimuth and range coordinates, respectively. Because a two-dimensional image, i.e., azimuth and range, is considered, the coordinate $\left(x_i,y_i,0\right)$ is simplified as $\left(x_i,y_i\right)$.

Similarly to the traditional SAR system, the raw data in the sub-Nyquist SAR for many point targets can be written as:

$$s_c\left(\tau ,\eta \right)={\displaystyle \sum _{x_i,y_i}{\sigma }_i}{W}_i\left(\tau ,\eta \right)\exp \left\{j\pi K_r{\left[\tau -2R_i\left(\eta \right)/c\right]}^2\right\}\cdot \exp \left\{-j4\pi R_i\left(\eta \right)/\lambda \right\}\cdot \exp \left\{j{\phi }_i\left(\eta \right)\right\}+n\left(\tau ,\eta \right)$$

(1)

where $\tau$ is the fast time along the range dimension, i.e., the sampling moment during one pulse width. $\eta$ is the slow time along the azimuth dimension, i.e., the moment of the transmitting pulse. ${\sigma }_i$ and $W_i\left(\tau ,\eta \right)$ are the backscattering coefficient and the weighting pattern corresponding to the i-th target at $\left(x_i,y_i\right)$, respectively. ${\phi }_i\left(\eta \right)$ is the random phase based on the pseudo-random time–space modulation corresponding to the target $\left(x_i,y_i\right)$. $K_r$ denotes the chirp rate of the linear frequency-modulated (LFM) signal, $c$ is the light speed, $\lambda$ is the wavelength, and $n\left(\tau ,\eta \right)$ is the system noise.

To relieve the inherent contradiction in the HRWS system, sub-Nyquist sampling in the sub-Nyquist SAR is implied along the azimuthal dimension, and the sampling method is as shown in Figure 1. The sub-Nyquist SAR system only adopts the traditional MF method to not recover the scene exactly [5]. Accordingly, sub-Nyquist SAR imaging includes three steps: (1) range compression based on MF; (2) range cell migration correction (RCMC): interpolate by zero padding every azimuthal signal after range compression and then sum up all the range-compressed signal of grids on the same range cell to the pre-defined grid; and (3) azimuth reconstruction with sub-Nyquist samples based on the CS algorithm. After implementing range compression and RCMC to Equation (1), the signal at a certain range cell is represented by:

$$\begin{array}{l}{s}_c\left({\tau }_0,\eta \right)\\ ={\displaystyle \sum _{x_i,y_i}{\sigma }_i}{W}_i\left({\tau }_0,\eta \right)T_r\sin c\left\{K_r{T}_r\left({\tau }_0-2\frac{R_i\left(\eta -{\eta }_{ci}\right)}{\mathrm{c}}\right)\right\}\cdot \exp \left\{-j4\pi \frac{R_i\left(\eta -{\eta }_{ci}\right)}{\lambda }\right\}\cdot \exp \left\{j{\phi }_i\left(\eta \right)\right\}+n\left({\tau }_0,\eta \right)\end{array}$$

(2)

where ${\eta }_{ci}$ is the beam center crossing time for the target $\left(x_i,y_i\right)$. $T_r$ is the pulse width and $\sin c\left(\cdot \right)$ denotes the sinc function.

Let $\sigma ={\left[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_M\right]}^T$ be the vectorized backscattering cross-sections of targets on the same range cell, and $s_{N\times 1}=[s_c\left({\tau }_0,{\eta }_1\right),s_c\left({\tau }_0,{\eta }_2\right),\cdots ,s_c\left({\tau }_0,{\eta }_N\right){]}^T$ be the vectored signal after range compression and RCMC; then,

$$s_{N\times 1}=D_{N\times M}{\sigma }_{M\times 1}+n_{N\times 1}$$

(3)

where $N$ is the sampling number on the azimuthal dimension and $M$ is the number of resolution cells at the certain range cell in the observed scene. $D_{N\times M}={\left\{D_i\left({\tau }_0,{\eta }_n\right)\right\}}_{n=1,i=1}^{N,M}$ denotes the mapping relation between the received signal and the scene, and $D_i\left({\tau }_0,{\eta }_n\right)=W_i\left({\tau }_0,{\eta }_n\right)\cdot T_r\cdot \exp \left\{-j4\pi R_i\left({\eta }_n-{\eta }_{ci}\right)/\lambda \right\}\cdot \exp \left\{j{\phi }_i\left(\eta \right)\right\}$ $n_{N\times 1}={\left[n\left({\tau }_0,{\eta }_1\right),n\left({\tau }_0,{\eta }_2\right),\cdots ,n\left({\tau }_0,{\eta }_N\right)\right]}^T$ is the noise.

2.1.2. Motion Error Model

The recovered scene ${\sigma }_{M\times 1}$ can be estimated by Equation (3), while the recovered matrix $D_{N\times M}$ can exactly reflect the relation between the raw data and the recovered scene. However, uncertainties and errors exist in the matrix $D_{N\times M}$, and the inaccuracy of the motion-induced model leads to the phase error [22]. If the preset matrix $D_{N\times M}$ in Equation (3) is still used to recover the scene ${\sigma }_{M\times 1}$ without extra processing steps, it may cause the defocusing of the reconstructed scene. A technology called autofocusing removes these phase errors.

To solve the problem of error correction, it is important to establish the exact mapping model so that we need to discern the uncertain factors. We mainly considered the position error and velocity error leading to the phase error in this paper. The SAR satellite is affected by different perturbations during on-orbit flying; so, it may cause the position error, as demonstrated in the following Figure 2, and the velocity error of the satellite. Firstly, we analyzed the position errors. The position error is mainly the deviation between the realistic position and the hypothetical position on the x-axis and z-axis. In Figure 2a, the solid line denotes the realistic track, and the dashed line denotes the hypothetical track. The hypothetical slant range between the radar and target without position error is:

$$R\left(\eta ,{\zeta }_R\right)=\sqrt{{\left(\frac{H}{\cos {\zeta }_R}\right)}^2+{\left(V_e\eta -V_e{\eta }_0\right)}^2}\approx \frac{H}{\cos {\zeta }_R}+\frac{{V_e}^2\cos {\zeta }_R}{2H}{\left(\eta -{\eta }_0\right)}^2$$

(4)

where $V_e$ is the radar velocity. ${\eta }_0$, the Doppler center moment, is the azimuthal moment when the radar is the nearest to the target.

When a position error exists, the real slant range $R_{\xi }\left(\eta ,{\zeta }_R\right)$ between the radar and target is:

$$\begin{array}{ll}{R}_{\xi }\left(\eta ,{\zeta }_R\right)& =\sqrt{{\left(H\tan {\zeta }_R-\Delta x\left(\eta \right)\right)}^2+{\left(V_e\eta -V_e{\eta }_0\right)}^2+{\left(H+\Delta z\left(\eta \right)\right)}^2}\\ & \approx \frac{H}{\cos {\zeta }_R}+\frac{{V_e}^2\cos {\zeta }_R}{2H}{\left(\eta -{\eta }_0\right)}^2+\Delta z\left(\eta \right)\cdot \cos {\zeta }_R-\Delta x\left(\eta \right)\cdot \sin {\zeta }_R\\ & \approx R\left(\eta ,{\zeta }_R\right)+\Delta R_{\xi }\left(\eta ,{\zeta }_R\right)\end{array}$$

(5)

$\Delta R_{\xi }\left(\eta ,{\zeta }_R\right)=\Delta z\left(\eta \right)\cdot \cos {\zeta }_R-\Delta x\left(\eta \right)\cdot \sin {\zeta }_R$ is the position error. Regardless of the weight $W_i\left({\tau }_0,{\eta }_n\right)\cdot T_r$ in Equation (2), the realistic echo signal $s_{error}\left(\eta ,{\zeta }_R\right)$ of target point P at a certain range cell after range compression and RCMC is denoted as

$$\begin{array}{ll}{s}_{error}\left(\eta ,{\zeta }_R\right)& ={\sigma }_P\cdot \exp \left\{-j\frac{4\pi }{\lambda }{R}_{\xi }\left(\eta ,{\zeta }_R\right)\right\}\\ & ={\sigma }_P\cdot \exp \left\{-j\frac{4\pi }{\lambda }\left(R\left(\eta ,{\zeta }_R\right)+\Delta R_{\xi }\left(\eta ,{\zeta }_R\right)\right)\right\}\\ & =s_0\left(\eta ,{\zeta }_R\right)\cdot H_{\Delta R}\left(\eta ,{\zeta }_R\right)\end{array}$$

(6)

where $s_0\left(\eta ,{\zeta }_R\right)={\sigma }_P\cdot \exp \left\{-j4\pi R\left(\eta ,{\zeta }_R\right)/\lambda \right\}$ is the raw signal without the position error. $H_{\Delta R}\left(\eta ,{\zeta }_R\right)$ is the position error signal.

Secondly, we analyzed the mathematical model for the velocity error. Regardless of the weight $W_i\left({\tau }_0,{\eta }_n\right)\cdot T_r$ in Equation (2), the realistic echo signal $s_{error}\left(\eta ,V_{Ee}\right)$ of target point P at a certain range cell after range compression and RCMC is:

$$s_{error}\left(\eta ,V_{Ee}\right)={\sigma }_P\cdot \exp \left\{-j\frac{4\pi R_{\xi }\left(\eta ,V_{Ee}\right)}{\lambda }\right\}$$

(7)

where $V_{Ee}$ is the equivalent velocity. The raw signal is approximate to the linear frequency modulation signal and is written as:

$$s_{error}\left(\eta ,V_{Ee}\right)=\exp \left\{f_d\eta +f_r\left(V_{Ee}\right){\eta }^2\right\}$$

(8)

where $f_d$ is the Doppler center frequency and the Doppler modulated rate is $f_r\left(V_{Ee}\right)=2{V_{Ee}}^2{\cos }^3{\zeta }_A/\lambda R$, and ${\zeta }_A$ is the squint angle on the azimuthal dimension. When the hypothetical equivalent velocity $V_E$ exists, $f_r\left(V_{Ee}\right)$ is unfolded at $V_{Ee}=V_E$:

$$f_r\left(V_{Ee}\right)=f_r\left(V_E\right)+{f_r}^{\prime }\left(V_E\right)\left(V_{Ee}-V_E\right)+\frac{1}{2!}{f_r}^{″}\left(V_E\right){\left(V_{Ee}-V_E\right)}^2+\cdots +\frac{1}{n!}{f_r}^{\left(n\right)}\left(V_E\right){\left(V_{Ee}-V_E\right)}^n+\cdots$$

(9)

If the error caused by the equivalent velocity is small, the above equation besides the first item $f_r\left(V_E\right)$ is also small and is denoted as $o_r\left(V_{Ee}\right)$:

$$o_r\left(V_{Ee}\right)={f_r}^{\prime }\left(V_E\right)\left(V_{Ee}-V_E\right)+\frac{1}{2!}{f_r}^{″}\left(V_E\right){\left(V_{Ee}-V_E\right)}^2+\cdots +\frac{1}{n!}{f_r}^{\left(n\right)}\left(V_E\right){\left(V_{Ee}-V_E\right)}^n+\cdots$$

(10)

Substituting (24) and (25) into (23),

$$\begin{array}{ll}{s}_{error}\left(\eta ,V_{Ee}\right)& =\exp \left\{f_d\eta +f_r\left(V_{Ee}\right){\eta }^2\right\}=\exp \left\{f_d\eta +f_r\left(V_E\right){\eta }^2+o_r\left(V_{Ee}\right){\eta }^2\right\}\\ & =\exp \left\{f_d\eta +f_r\left(V_E\right){\eta }^2\right\}\cdot \exp \left\{o_r\left(V_{Ee}\right){\eta }^2\right\}\\ & =s_0\left(\eta ,V_E\right)\cdot \exp \left\{o_r\left(V_{Ee}\right){\eta }^2\right\}=s_0\left(\eta ,V_E\right)\cdot H_{V_{Ee}}(\eta )\end{array}$$

(11)

where $s_0\left(\eta ,V_E\right)=\exp \left(f_d\eta +f_r\left(V_E\right){\eta }^2\right)=\exp \left\{-j4\pi R\left(\eta ,{\zeta }_R\right)/\lambda \right\}$ is the azimuthal signal without error and $H_{V_{Ee}}(\eta )=\exp (o_r\left(V_{Ee}\right){\eta }^2)$ is the error signal caused by the equivalent velocity $V_{Ee}$.

According to Equations (21) and (26), the azimuthal signal with the error can be written as the multiplication between the azimuthal signal without the error and the phase error. Considering the above phase error, the received raw data without the weight $W_i\left({\tau }_0,{\eta }_n\right)\cdot T_r$ in Equation (2) after range compression and RCMC at a certain range cell is approximately denoted as:

$$s_{error}\left({\tau }_0,\eta \right)={\displaystyle \sum _{i=1}^M{\sigma }_i\exp \left\{-j4\pi \frac{R_i\left(\eta \right)}{\lambda }\right\}}\cdot \exp \left\{j\vartheta \right\}+n\left({\tau }_0,\eta \right)$$

(12)

where $\vartheta$ is the error phase variation with the azimuthal sampling moment. Similarly to Equation (3), the vector matrix form of the raw signal (27) with an error along the azimuthal sampling moment $\eta$ is:

$$s_{N\times 1}^{error}=E_{N\times N}{D}_{N\times M}{\sigma }_{M\times 1}+n_{N\times 1}$$

(13)

where $E_{N\times N}=diag\left(\exp \left(j{\vartheta }_1\right),\exp \left(j{\vartheta }_2\right),\cdots ,\exp \left(j{\vartheta }_N\right)\right)$. A different error has a different expression form of the error matrix $E_{N\times N}$, and the motion error has been not fully formulated yet due to some approximations in the models. Therefore, we assigned an error matrix $E_{N\times N}$ random phase in the following simulation.

2.2. Sub-Nyquist SAR Imaging and Error Correction Based on the Pseudo-$ℒ$₀-Norm Optimization Algorithm

In this subsection, the CS theorem is first briefly introduced. Then, a pseudo-$ℒ$₀-norm regularization model is presented and solved by the Cauchy–Newton method in detail. Finally, sub-Nyquist SAR imaging and error correction based on pseudo-random space-time modulation are proposed based on our proposed pseudo-$ℒ$₀-norm optimization algorithm.

2.2.1. CS Theorem

The key to the CS theorem is the effective recovered algorithm. While RIP is satisfied, there are three recovered algorithms, i.e., greedy algorithm, $ℒ$₁-norm optimization algorithm, and sparse Bayesian learning method. The $ℒ$₁-norm optimization algorithm is robust in the sense that it can effectively recover nearly sparse signals with/without measurement noise from remarkably few measurements, and its application is so wide that it broadly could be considered the modern least squares [12,13]. The $ℒ$₁-norm optimization algorithm has a better-recovered performance in terms of the recovered error evaluated by MSE [17]. Its optimization equation is to solve the underdetermined problem (3) without a subscript:

$$\widehat{\sigma }=\underset{\sigma }{\arg \min }\left\{{‖s-D\sigma ‖}_2^2+\alpha {‖\sigma ‖}_1\right\}$$

(14)

where $\alpha$ is the regularization parameter. The first item ${‖s-D\sigma ‖}_2^2$ ensures the recovery error, and the second item ${‖\sigma ‖}_1$ guarantees the sparsity of the recovered scene. The parameter $\alpha$ balances between the recovered error and the sparsity, and is empirically chosen by minimizing the recovered error of the whole scene.

Even though the $ℒ$₁-norm optimization algorithm achieved a good recovered performance, some low SNR targets cannot be accurately recovered, and it often yields false targets for an unfair penalization rule. The further optimization of the regularization model contributes to the recovered scene so that false targets can be removed, and then we carried out research about the optimization model. Is there an alternative to the $ℒ$₁-norm regularization model to achieve a better-recovered performance? The problem of further improving the optimization algorithm is thus brought up. To solve this problem, a reweighted $ℒ$₁-based algorithm has been already proposed to improve the recovered performance [20,21]. However, how to select an approximately fair regularization rule still does not have a theoretical analysis.

2.2.2. Sub-Nyquist SAR Imaging Based on the Pseudo-L₀-Norm Optimization Algorithm

For a sparse undetermined equation, the $ℒ$₀-norm optimization algorithm should have the best recovered performance, but this optimization equation is a non-polynomial hard (NP-hard) problem [12]. Based on the idea of adaptive least absolute shrinkage and selection operator (Lasso) technique that is a popular technique for simultaneous estimation and variable selection [21], this paper penalized the regularization item to be close to the $ℒ$₀-norm to achieve a better performance of the $ℒ$₀-based algorithm, and the regularization model can be solved by the Cauchy–Newton method. Since it is closer to the $ℒ$₀-based method in some sense, we named it a pseudo-$ℒ$₀-norm optimization algorithm. In the following, we adopted Bayesian estimation to analyze and deduct a pseudo-$ℒ$₀-norm regularization model by making full use of the scene prior information. The rule was deducted as follows: Usually, the noise $n_{N\times 1}$ is assumed to be a Gaussian distribution with zero mean and variance ${{\sigma }_n}^2$ [14]. For simplicity, all the following matrixes/vectors omit subscripts.

$$p_n\left(n\right)=p_{s/\sigma }\left(s/\sigma \right)={\left(\frac{1}{\sqrt{2\pi }{\sigma }_n}\right)}^N\exp \left\{-\frac{{‖s-D\sigma ‖}_2^2}{2{{\sigma }_n}^2}\right\}$$

(15)

Laplace distribution forces most coefficients to be small so that it can describe the sparse scene [35]. We assume:

$$p\left({\sigma }_i\right)=\left(\frac{{\xi }_i}{2}\right)\exp \left\{-{\xi }_i{‖{\sigma }_i‖}_1\right\}$$

(16)

where ${\xi }_i$ is the scale parameter of Laplace distribution and ${\xi }_i>0$. Then, the probability distribution of the vector ${\sigma }_{M\times 1}$ is:

$$p_{\sigma }\left(\sigma \right)=\prod _{i=1}^M\left\{\left(\frac{{\xi }_i}{2}\right)\exp \left\{-{\xi }_i{‖{\sigma }_i‖}_1\right\}\right\}=\prod _{i=1}^M\left\{\left(\frac{{\xi }_i}{2}\right)\right\}\exp \left\{-{\displaystyle \sum _{i=1}^M{\xi }_i{‖{\sigma }_i‖}_1}\right\}$$

(17)

Based on the Bayesian rule in the information theory, the maximum posterior (MAP) probability of the vector ${\sigma }_{M\times 1}$ is:

$$\widehat{\sigma }=\arg \max \left[p_{\sigma /s}\left(\sigma /s\right)\right]=\arg \max \left[p_{s/\sigma }\left(s/\sigma \right)\cdot p_{\sigma }\left(\sigma \right)\right]$$

(18)

We took the logarithm of the above Formula (18):

$$\begin{array}{ll}\widehat{\sigma }& =\underset{\sigma }{\arg \max }\left[\log \left[p_{s/\sigma }\left(s/\sigma \right)\right]+\log \left[p_{\sigma }\left(\sigma \right)\right]\right]\\ & =\arg \max \left[-\frac{1}{2{{\sigma }_n}^2}{‖s-D\sigma ‖}_2^2-{‖\xi \cdot \sigma ‖}_1\right]\\ & =\arg \min \left[{‖s-D\sigma ‖}_2^2+\beta {‖\xi \cdot \sigma ‖}_1\right]\end{array}$$

(19)

where the reweighting matrix $\xi$ is the diagonal matrix with $\left\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_M\right\}$ on the diagonal and zeros elsewhere. Similarly to $\alpha$ in Equation (4), $\beta$ is the regularization parameter. After the above deduction, the regularization model (19) is the modification to the $ℒ$₁-norm penalization rule.

Additionally, the matrix $\xi$ is first calculated to solve Equation (19). The logarithm function of the vector ${\sigma }_{M\times 1}$ is:

$$\begin{array}{l}\log \left\{p_{\sigma }\left(\sigma \right)\right\}=\log \left\{{\displaystyle \prod _{i=1}^M}\left\{\left(\frac{{\xi }_i}{2}\right)\exp \left\{-{\xi }_i{‖{\sigma }_i‖}_1\right\}\right\}\right\}\\ =\log \left\{{\displaystyle \prod _{i=1}^M}\left(\frac{{\xi }_i}{2}\right)\right\}+\log \left\{\exp \left\{-{\displaystyle \sum _{i=1}^M{\xi }_i{‖{\sigma }_i‖}_1}\right\}\right\}={\displaystyle \sum _{i=1}^M\log }\left(\frac{{\xi }_i}{2}\right)-{\displaystyle \sum _{i=1}^M{\xi }_i{\left|{\sigma }_i\right|}_1}\end{array}$$

(20)

Let the partial derivative function of Formula (10) with respect to ${\xi }_i$ be equal to zero, and the estimation of the scale parameter ${\xi }_i$ is:

$$\frac{\partial \log \left[p_{\sigma }\left(\sigma \right)\right]}{\partial {\xi }_i}=\frac{1}{{\xi }_i}-{‖{\sigma }_i‖}_1=0$$

(21)

$${\xi }_i=\frac{1}{\left|{\sigma }_i\right|}$$

(22)

In the case that ${\sigma }_i=0$, Formula (12) makes no sense. Formula (12) is modified as:

$${\xi }_i=\frac{1}{\left|{\sigma }_i\right|+\iota }$$

(23)

where $\iota >0$ is a very small positive constant.

The above deduction explains the prior distribution, e.g., Laplace distribution in the sub-Nyquist SAR confirms the penalization rule and this method makes full use of the data prior information to achieve a good performance. Certainly, different probability distributions confirm different penalization rules. In the scenario of not knowing the variables themselves, using the iteratively updated method to establish an approximately fair penalization rule allows for the successively better estimation of nonzero variables. According to the optimization Equations (19) and (23), we know that the regularized item is iteratively penalized by itself to nearly acquire the number of nonzero values, and the large coefficients are more heavily penalized to discourage their effects compared to that with a small coefficient and is more likely to be identified as nonzero. Once the nonzero locations are identified, their influence attenuates to allow more sensitivity to identify the remaining small but nonzero elements. It means this algorithm can more accurately recover low SNR targets.

To solve the pseudo-$ℒ$₀-norm regularization equation, the algorithm in Equation (19) is essentially a $ℒ$₁ penalization method so that it is convex and can be solved by the Cauchy–Newton method [36,37]. The solving method is as follows: $ℒ$₁ norm is not differentiable, while ${\sigma }_i=0$. Firstly, the smoothing approximation is introduced:

$${‖\sigma ‖}_1={\displaystyle \sum _{i=1}^N\sqrt{\left({\left|{\sigma }_i\right|}^2+\varsigma \right)}}$$

(24)

where $\varsigma >0$ is a very small positive constant. Equation (19) can be expressed as:

$$\widehat{\sigma }=\underset{\overline{\sigma }}{\arg \min }\left\{{‖s-D\sigma ‖}_2^2+\beta {\displaystyle \sum _{i=1}^M{\xi }_i\sqrt{\left({\left|{\sigma }_i\right|}^2+\varsigma \right)}}\right\}$$

(25)

where $f\left(\widehat{\sigma }\right)={‖s-D\sigma ‖}_2^2+\beta {\displaystyle \sum _{i=1}^M{\xi }_i\sqrt{\left({\left|{\sigma }_i\right|}^2+\varsigma \right)}}$. The conjugate gradient function is written as:

$${\nabla }_{\widehat{\sigma }}f\left(\widehat{\sigma }\right)=2D^{H}D\sigma +\beta U\left(\widehat{\sigma }\right)\cdot \xi \sigma -2D^{H}s=H\left(\widehat{\sigma }\right)\sigma -2D^{H}s$$

(26)

where $H\left(\widehat{\sigma }\right)=2D^{H}D+\beta U\left(\widehat{\sigma }\right)\xi$, $U\left(\widehat{\sigma }\right)=diag\left\{1/\left(\sqrt{{\left|{\sigma }_i\right|}^2+\varsigma }\right)\right\},i=1,2,\cdots ,M$. $D^H$ is the conjugate transpose of the observed matrix $D$.

According to the Newton method [36,37], the iterative equation of the reconstructed result $\sigma$ is:

$$H\left({\widehat{\sigma }}^{\left(g\right)}\right){\widehat{\sigma }}^{\left(g+1\right)}=\left(1-\gamma \right)H\left({\widehat{\sigma }}^{\left(g\right)}\right){\widehat{\sigma }}^{\left(g\right)}+2\gamma D^{H}s$$

(27)

where $\gamma$ is the step size of the iteration. When $\gamma =1$, Equation (27) is simplified as:

$${\widehat{\sigma }}^{\left(g+1\right)}=2H{\left({\widehat{\sigma }}^g\right)}^{-1}{D}^{H}s$$

(28)

Until the iteration terminates, $\widehat{\sigma }$ is the recovered result. Normally, the MSE, the iterative number, etc., can be taken as the iterative criterion [23]. After the above deduction, the flow chart is concluded as follows:

(1)

Initialization: the iterative step $g=1$, ${\widehat{\sigma }}^1=D^{H}s$;

(2)

Updating of the weighting matrix $\xi$ and the matrix $H\left(\widehat{\sigma }\right)$: ${\xi }^{\left(g\right)}=diag\left(\frac{1}{\left|{\sigma }_1\right|+\iota },\frac{1}{\left|{\sigma }_2\right|+\iota },\cdots ,\frac{1}{\left|{\sigma }_M\right|+\iota }\right)$, $H{\left(\widehat{\sigma }\right)}^{\left(g\right)}=2D^{H}D+\beta U\left(\widehat{\sigma }\right){\xi }^{\left(g\right)}$;

(3)

Calculation: ${\widehat{\sigma }}^{\left(g+1\right)}=2{\left(H{\left(\widehat{\sigma }\right)}^{\left(g\right)}\right)}^{-1}{D}^{H}s$;

(4)

$g=g+1$;

(5)

Loop;

(6)

Stopping iteration according to the iterative criterion.

During the deduction of the pseudo-$ℒ$₀-norm optimization algorithm, there are four undetermined parameters, i.e., the regularization factor $\beta$, the parameter $\iota$, the parameter $\varsigma$, and the iterative criterion. Similarly to the regularization $\alpha$, $\beta$ is also an empirical value and is chosen by minimizing the recovered error [12]. In the simulation section, the iterative criterion selects the preset iterative number and this number is chosen empirically. Although the reference [23] demonstrates how to choose the parameter $\iota$ and $\varsigma$, the pseudo-$ℒ$₀-norm optimization algorithm has a good robustness and the reconstructed scene does have not the strictly sparse property; so, the proposed algorithm sets the parameters $\iota$ and $\varsigma$ as ${10}^{-3}$ and ${10}^{-6}$, respectively.

2.2.3. Error Correction Based on the Pseudo-$ℒ$₀-Norm Optimization Algorithm

Based on the observation model of Equation (13), we propose an error correction method integrated with the pseudo-$ℒ$₀-norm optimization algorithm. This method considers the phase error as the model error and removes it during the scene reconstruction, and it enables our method to correct more artifacts due to the robustness of the pseudo-$ℒ$₀-norm optimization algorithm. This error correction method is an update and iteration algorithm, and it includes the scene reconstruction and the error estimation during each iteration. In the first step of every iteration, the cost function is minimized with the scene based on the reweighted $ℒ$₁-norm optimization algorithm, and in the second step, the phase error is estimated given the scene estimate. According to the estimated error, the error matrix $E_{N\times N}$ in (13) is updated and delivered to the next iteration. The procedure is as follows:

(1)

The scene reconstruction

Based on the pseudo-$ℒ$₀-norm optimization algorithm, the optimization equation of (13) without subscripts is:

$$\widehat{\sigma }=\underset{\sigma }{\arg \min }\left[{‖s^{error}-ED\sigma ‖}_2^2+\beta {‖\xi \sigma ‖}_1\right]$$

(29)

The flow chart in Section 2.2.2 can resolve this optimization equation.

(2)

The error estimation

To estimate the phase error, the cost function should minimize the MSE of the reconstructed result since the error exists at each sampling moment, and each element can be handled separately. The cost function at each sampling moment is:

$${\widehat{E}}_n\left(g+1\right)=\underset{E_n}{\arg \min }{{‖s_n^{error}-E_n{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n‖}_2}^2,n\in \left\{1,2,\cdots ,N\right\}$$

(30)

where ${\widehat{E}}_n\left(g+1\right)$ is the n-th diagonal element of the error matrix ${\widehat{E}}^{\left(g+1\right)}$ at the (g + 1)-th iteration. In the following, the cost function is unfolded for analysis:

$$\begin{array}{l}{{‖s_n^{error}-E_n{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n‖}_2}^2={\left(s_n^{error}-E_n{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right)}^H\cdot \left(s_n^{error}-E_n{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right)\\ ={\left(s_n^{error}\right)}^H\cdot s_n^{error}+{{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n-2\cos \left({\vartheta }_n\right)\cdot \mathrm{Re}\left\{{\left(s_n^{error}\right)}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right\}+2\sin \left({\vartheta }_n\right)\cdot \mathrm{Im}\left\{{\left(s_n^{error}\right)}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right\}\end{array}$$

(31)

where $\mathrm{Re}\left\{·\right\}$ and $\mathrm{Im}\left\{·\right\}$ are the real part and the imaginary part of a complex signal, respectively. Assuming that $\rho =\mathrm{Re}\left\{{\left(s_n^{error}\right)}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right\}$ and $\upsilon =\mathrm{Im}\left\{{\left(s_n^{error}\right)}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n\right\}$, the above equation can be written as:

$${{‖s_n^{error}-E_n{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n‖}_2}^2={\left(s_n^{error}\right)}^H\cdot s_n^{error}+{{\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n-2\sqrt{{\rho }^2+{\upsilon }^2}\cdot \cos \left({\vartheta }_n-\arctan \left(\frac{\upsilon }{\rho }\right)\right)$$

(32)

Estimating ${\vartheta }_n$ is equal to minimizing the above cost function (32); so, ${\widehat{\vartheta }}_n$ satisfies:

$${\widehat{\vartheta }}_n=\angle {\left(s_n^{error}\right)}^H\cdot {\left(D^{\left(g\right)}{\sigma }^{\left(g+1\right)}\right)}_n$$

(33)

After the above deduction, the detailed flow chart is as follows:

(1)

Initialization: $g=0$, $E=I$;

(2)

Recovering the scene: ${\widehat{\sigma }}^{\left(g\right)}=\underset{\sigma }{\arg \min }{‖s^{error}-E^{\left(g\right)}D{\sigma }^{\left(g\right)}‖}_2+\beta {‖{\xi }^{\left(g\right)}{\sigma }^{\left(g\right)}‖}_1$;

(3)

Estimating the error matrix: ${\widehat{E}}^{\left(g\right)}=\underset{E}{\arg \min }{‖s^{error}-E^{\left(g\right)}D{\sigma }^{\left(g\right)}‖}_2$;

(4)

Updating $E^{\left(g\right)}$;

(5)

$g=g+1$;

(6)

Loop;

(7)

Stopping iteration according to the iterative criterion.

where ${\widehat{\sigma }}^{\left(g\right)}$ and ${\widehat{E}}^{\left(g\right)}$ denote the recovered result and the estimated error matrix in the g-th iteration, respectively. During each iteration, the cost function focuses part of the defocused scene to generate a relatively accurate error matrix so that a relatively accurate model can lead to a better reconstruction result. The iterative criterion was mentioned in Section 2.2.2 and is not repeated in this section.

2.2.4. Analysis of the Computational Complexity

In order to obtain the total computational complexity of the proposed method, we first analyzed the time complexity per iteration quantitatively. As the notations used above, the total computational cost of the proposed method is in the order of $o\left(G\cdot NM\log \left(NM\right)\right)$, where $G$ is the required number of iterations to recover the result in one range cell. The value of $G$ is difficult to confirm accurately through theoretical analysis, but in practice, the running time of the proposed method is affordable. As for the memory cost of the proposed method, we only need to store the input, the output, and the parameter matrices. In summary, the memory cost is in the order of $o\left(NM\right)$.

3. Experiment

In this section, simulation experiments and data experiments verify the effectiveness and superiority of sub-Nyquist SAR imaging and error correction integrated with a pseudo-$ℒ$₀-norm optimization algorithm.

3.1. Data Description

To verify the validity and effectiveness of sub-Nyquist SAR imaging and error correction based on the pseudo-$ℒ$₀-norm optimization algorithm, we selected real strip-map TerraSAR-X images for the experiments. The reflectivity functions of images were used to simulate raw data according to the simulated parameters, and then the raw data were uniformly received, as shown in Figure 1. During the error correction, the error is added at each sampling moment and the error adopts the random distribution.

Figure 3 features a sea–land interface scene. This scene is more complex and used to verify the performance of our proposed algorithm in the sub-Nyquist SAR imaging compared with the $ℒ$₁-norm optimization algorithm. Figure 4 is a sea containing several boats to verify the availability of error corrections based on the pseudo-$ℒ$₀-norm optimization algorithm.

3.2. The Simulation of the Pseudo-$ℒ$₀-Norm Optimization Algorithm

We used a one-dimensional simulation of point targets to verify the superiority of our proposed algorithm compared with the $ℒ$₁-norm optimization algorithm. For a given under-recovered signal, $\sigma$, the recovered performance can be evaluated by the normalized MSE (NMSE):

$$NMSE=\frac{{‖\widehat{\sigma }-\sigma ‖}_2^2}{{‖\sigma ‖}_2^2}$$

(34)

where $\widehat{\sigma }$ is the recovered result.

The following one-dimensional simulation illustrates that the pseudo-$ℒ$₀-norm algorithm has a better performance compared with the $ℒ$₁-norm optimization algorithm. The simulated parameters are as follows. It selects a sparse signal, $\sigma$, of length $M=256$, with ${‖\sigma ‖}_0=S$. The $S$ nonzero spike positions are selected randomly, and the amplitude of nonzero elements obeys a zero-mean unit-variance Gaussian distribution. This simulation selects the measured number and a $N\times M$ random matrix, $D$, with independent identically distributed (i.i.d.) Gaussian elements. The noise vector, $n$, is drawn from the i.i.d. zero-mean Gaussian function with ${‖n‖}_2=0.5$, so that $SNR=10log10({‖\sigma ‖}_2^2/{‖n‖}_2^2)=21dB$. To recover the signal, $\sigma$, it adopts two algorithms, the $ℒ$₁-norm optimization algorithm and pseudo-$ℒ$₀-norm algorithm, to compare the recovered performance. The NMSE of the two algorithms and the recovered results are compared in Figure 5 and Figure 6.

Figure 5 demonstrates that the pseudo-$ℒ$₀-norm optimization algorithm can achieve the smallest recovered error, although this proposed algorithm has a slower rate of convergence. Our proposed method is sensitive to low SNR targets that are not easy to identify, so the convergence rate is slow. Figure 6 shows the recovered results under different algorithms. In the red circle, ①, the results for the $ℒ$₁-norm optimization algorithm have weak false targets, but two other algorithms do not. In circles ② and ④, the weak target can be recovered under the optimization algorithm, while the $ℒ$₁-norm optimization algorithm cannot recover it. In circle ③, the amplitude of targets can be exactly recovered by the proposed algorithm. It demonstrates that our proposed algorithm is more friendly to the low SNR target and removes false targets to achieve a good recovered performance.

3.3. The Simulation of Sub-Nyquist SAR Imaging

We used a simulation based on the real strip-map TerraSAR-X image in Figure 3 to verify the validity of our proposed algorithm for sub-Nyquist SAR imaging. The SAR image is over the sea–land interface scene. The simulated parameters are shown in

Table 1. Simulated parameters.

Parameter	Data
Average PRF (Hz)	893
Range sampling frequency (MHz)	55
Referred slant range (km)	870
Chirp rate (Hz/s)	1012
Doppler bandwidth (Hz)	2438
Wavelength (mm)	5.55
Velocity (m/s)	7513
Height (km)	693
Squint angle (°)	0

, and the raw data were randomly received, as shown in Figure 1. In the simulation, the interfaced land is recovered without the loss of details under the proposed algorithm. According to the simulated result in Figure 7 and MSE in

Table 2. MSE under different algorithms.

	pseudo-$ℒ$0-norm optimization algorithm	$ℒ$1-norm optimization algorithm
MSE	0.114	0.874

, this simulation also demonstrates that the pseudo-$ℒ$₀-norm optimization algorithm can achieve a better recovered performance compared with the $ℒ$₁-norm optimization algorithm.

3.4. Simulation of the Error Correction

The error correction mainly solves image defocusing for the phase error. As we all know, the better the focusing performance is, the smaller the image entropy is. We adopted image entropy $H\left(\sigma \right)$ to evaluate the reconstruction result.

$$H\left(\sigma \right)=-{\displaystyle \sum _i\frac{{\left|{\sigma }_i\right|}^2}{{\displaystyle \sum _i{\left|{\sigma }_i\right|}^2}}}{\log }_2\frac{{\left|{\sigma }_i\right|}^2}{{\displaystyle \sum _i{\left|{\sigma }_i\right|}^2}}$$

(35)

Firstly, we simulated point targets to illustrate the effectiveness of our method. The simulated parameters are shown in

Table 3. Simulated parameters.

Parameter	Data
Average PRF (Hz) in the sub-Nyquist SAR	155
PRF (Hz) in the traditional SAR	1907
Range sampling frequency (MHz)	120
Referred slant range (km)	888
Pulse width (us)	50
Doppler bandwidth (Hz)	1401
Wavelength (mm)	5.55
Velocity (m/s)	7513
Height (km)	693
Squint angle (°)	0

. When the elements of the error matrix are uniformly distributed $\left[0,17/18\pi \right]$, the simulated results are presented in Figure 8. Figure 8a is the original image, (b) is the reconstructed scene without an error correction in the traditional SAR, (c) is the reconstructed scene with an error correction in the sub-Nyquist SAR based on the $ℒ$₁–norm optimization algorithm, and (d) is the reconstructed scene with the error correction in the sub-Nyquist SAR based on the pseudo-$ℒ$₀-norm optimization algorithm. When the phase error exists, the reconstructed scene without an error correction is defocused (Figure 8b). Although it can also remove errors and focus the scene based on the $ℒ$₁-norm optimization algorithm, there are false targets and lost targets in the red circle (Figure 8c). Our method almost removes errors and recovers the scene based on the pseudo-$ℒ$₀-norm optimization algorithm (Figure 8d). From the quantitative analysis,

Table 4. Image entropy under different algorithms.

	The original image	The reconstructed scene without an error correction	The reconstructed scene with an error correction based on the $ℒ$1-norm optimization algorithm	The reconstructed scene with an error correction based on the pseudo-$ℒ$0-norm optimization algorithm
Image entropy (bit)	4.00	5.90	3.08	3.99

illustrates that the image entropy based on our method achieves the nearly same image entropy as the original image.

The simulation in Figure 9 is based on a strip-map TerraSAR-X image, i.e., sea containing several boats in Figure 4. The reflectivity function of the image was used to simulate raw data, and the simulated parameters are shown in Table 3. The elements of the error matrix are uniformly distributed, $\left[0,\pi /2\right]$. Figure 9a is the recovered result in the traditional SAR, but is defocused for the phase error. The error is partly removed from the Sub-Nyquist SAR with an error correction method based on the $ℒ$₁-norm optimization algorithm in Figure 9b. Based on the pseudo-$ℒ$₀-norm optimization algorithm, it almost removes the error and focuses the scene in Figure 9c. The image entropy in

Table 5. Image entropy under different algorithms.

	The original image	The reconstructed scene without an error correction	The reconstructed scene with an error correction based on the $ℒ$1-norm optimization algorithm	The reconstructed scene with an error correction based on the pseudo-$ℒ$0-norm optimization algorithm
Image entropy (bit)	7.81	12.16	8.35	7.80

illustrates that the image recovered by our method can achieve a similar performance to the original image.

4. Discussion

HRWS SAR imaging has always been the goal of spaceborne SAR systems in remote sensing applications [2]. Since a high resolution and wide swath are inherently conflicting requirements, the HRWS SAR system creates new difficulties and challenges. These requirements are simultaneously satisfied by advanced observing modes, e.g., azimuthal multi-channel SAR [3] and MIMO SAR [4], and these observing modes both have the characteristics of large amounts of data and long antenna. As the CS theorem develops, a novel imaging mode named sub-Nyquist SAR based on the pseudo-random space–time modulation has been proposed without large data and long antenna, and sub-Nyquist SAR imaging based on the CS algorithms can recover the scene [8].

Although the $ℒ$₁-norm optimization algorithm has achieved a quite good performance, some low SNR targets are not accurately recovered and it yields false targets. To further improve the recovered performance of sub-Nyquist SAR, we present a pseudo-$ℒ$₀-norm optimization algorithm based on the Bayesian estimation. This algorithm penalizes the regularization item with the reciprocal of its previous solution to acquire nonzero values so that the rule is fair and makes full use of the data prior information. This proposed algorithm essentially adopts the idea of the adaptive Lasso technique [20,21], and the $ℒ$₁-norm optimization algorithm is the Lasso variable selection method [13]. The Lasso variable selection is consistent with satisfying a necessary condition; sometimes Lasso is not consistent and does not have predictive properties. Adaptive Lasso is the modified version of the Lasso system with important differences. The modification is data-dependent and reasonably selected so that adaptive Lasso has predictive properties [20,21]. Our proposed algorithm is better in terms of consistency. For the recovered error, we deducted the expression of MSE in the following way. Assume that $F_{N\times S}$ is the submatrix constructed by taking $S$ columns from the recovered matrix, $D_{N\times M}$, which are specified by the index vector, $\Lambda$, and each element in $\Lambda$ satisfies the following:

$${\sigma }_{{\Lambda }_l}\ne 0,\left(l=1,2,\cdots ,S\right)$$

(36)

The estimated covariance matrix for the nonzero components, ${\sigma }_l$, is

$$C={\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}{F_{N\times S}}^H{F}_{N\times S}{\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}$$

(37)

where $\Sigma \left({\sigma }_{{\Lambda }_l}\right)=diag\left({\xi }_{{\Lambda }_1}/{\sigma }_{{\Lambda }_1},{\xi }_{{\Lambda }_2}/{\sigma }_{{\Lambda }_2},\cdots ,{\xi }_{{\Lambda }_S}/{\sigma }_{{\Lambda }_S}\right)$ is a diagonal matrix. The estimated error is the trace of the covariance matrix, $C$:

$$\begin{array}{l}{‖{\widehat{\sigma }}_{\Lambda }-{\sigma }_{\Lambda }‖}_2^2={{\sigma }_n}^{2}trace\left(C\right)\\ ={{\sigma }_n}^{2}trace\left({\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}{F_{N\times S}}^H{F}_{N\times S}{\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}\right)\\ ={{\sigma }_n}^{2}trace\left({\left({\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}\right)}^2{F_{N\times S}}^H{F}_{N\times S}\right)\\ \ge {{\sigma }_n}^2\left[2trace\left({\left({F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)\right)}^{-1}\right)-trace\left({\left({F_{N\times S}}^H{F}_{N\times S}\right)}^{-1}\right)\right]\end{array}$$

(38)

Assuming that $A={F_{N\times S}}^H{F}_{N\times S}$, $B={F_{N\times S}}^H{F}_{N\times S}+{\gamma }_2\Sigma \left({\sigma }_{{\Lambda }_l}\right)$, and their eigenvalues are ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_S$ and ${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_S$, respectively. Based on the structure form of the recovered matrix, $D_{N\times M}$, the $ℒ$₂-norm of each column is nearly equal and assumes ${‖D_l‖}_2^2\approx L$, where $L$ is the number of samples in one aperture time [8]. Formula (38) is simplified as:

$${‖{\widehat{\sigma }}_{\Lambda }-{\sigma }_{\Lambda }‖}_2^2\ge {{\sigma }_n}^2\left[2trace\left(B^{-1}\right)-trace\left(A^{-1}\right)\right]$$

(39)

Gail’s circle theorem [38] indicates that

$$\begin{array}{c}\left|{\lambda }_l-L\right|\le \left(S-1\right)\cdot u\cdot L\\ \left(1-\left(S-1\right)\cdot u\right)\cdot L\le {\lambda }_l\le \left(1+\left(S-1\right)\cdot u\right)\cdot L\\ \frac{1}{\left(1+\left(S-1\right)\cdot u\right)\cdot L}\le {{\lambda }_l}^{-1}\le \frac{1}{\left(1-\left(S-1\right)\cdot u\right)\cdot L}\\ \frac{S}{\left(1+\left(S-1\right)\cdot u\right)\cdot L}\le {\displaystyle \sum _{l=1}^S{{\lambda }_l}^{-1}}\le \frac{S}{\left(1-\left(S-1\right)\cdot u\right)\cdot L}\\ \frac{S}{\left(1+\left(S-1\right)\cdot u\right)\cdot L}\le trace\left(A^{-1}\right)\le \frac{S}{\left(1-\left(S-1\right)\cdot u\right)\cdot L}\end{array}$$

(40)

$$\begin{array}{c}\left|{\alpha }_l-\left(L+\frac{{\gamma }_2\cdot {\xi }_{{\Lambda }_l}}{{{\sigma }_{\Lambda }}_l}\right)\right|\le \left(S-1\right)\cdot u\cdot L\\ \frac{{\sigma }_{{\Lambda }_l}\cdot \left[1-\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot {\xi }_{{\Lambda }_l}}{{\sigma }_{{\Lambda }_l}}\le {\alpha }_l\le \frac{{\sigma }_{{\Lambda }_l}\cdot \left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot {\xi }_{{\Lambda }_l}}{{\sigma }_{{\Lambda }_l}}\\ \frac{{\sigma }_{{\Lambda }_l}}{{\sigma }_{{\Lambda }_l}\cdot \left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot {\xi }_{{\Lambda }_l}}\le {{\alpha }_l}^{-1}\le \frac{{\sigma }_{{\Lambda }_l}}{{\sigma }_{{\Lambda }_l}\cdot \left[1-\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot {\xi }_{{\Lambda }_l}}\\ \frac{S}{\left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \max \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}\le {\displaystyle \sum _{l=1}^S{{\alpha }_l}^{-1}}\le \frac{S}{\left[1-\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \min \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}\\ \frac{S}{\left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \max \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}\le trace\left(B^{-1}\right)\le \frac{S}{\left[1-\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \min \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}\end{array}$$

(41)

where $u=\underset{1\le m_1\ne m_2\le M}{\max }\left|〈D_{m_1},D_{m_2}〉\right|/\left({‖D_{m_1}‖}_2\cdot {‖D_{m_2}‖}_2\right)$ is the mutual coherence coefficient that reflects the maximum similarity between any two different columns, $m_1$, $m_2$, in the recovered matrix, $D_{N\times M}$. Combine Formulas (40) and (41):

$$\frac{2S}{\left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \max \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}-\frac{S}{\left(1-\left(S-1\right)\cdot u\right)\cdot L}\le 2trace\left(B^{-1}\right)-trace\left(A^{-1}\right)\le \frac{2S}{\left[1-\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \min \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}-\frac{S}{\left(1+\left(S-1\right)\cdot u\right)\cdot L}$$

(42)

So, we have the MSE of the proposed algorithm:

$${‖{\widehat{\sigma }}_{\Lambda }-{\sigma }_{\Lambda }‖}_2^2\ge {{\sigma }_n}^2\left[2trace\left(B^{-1}\right)-trace\left(A^{-1}\right)\right]\ge {{\sigma }_n}^2\left[\frac{2S}{\left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_2\cdot \max \left(\raisebox{1ex}{${\xi }_{{\Lambda }_l}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}-\frac{S}{\left(1-\left(S-1\right)\cdot u\right)\cdot L}\right]$$

(43)

Similarly, the MSE of the $ℒ$₁-norm optimization algorithm is:

$${‖{\widehat{\sigma }}_{l_1}-{\sigma }_{\Lambda }‖}_2^2\ge {{\sigma }_n}^2\left[2trace\left(B^{-1}\right)-trace\left(A^{-1}\right)\right]\ge {{\sigma }_n}^2\left[\frac{2S}{\left[1+\left(S-1\right)\cdot u\right]\cdot L+{\gamma }_1\cdot \max \left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{{\Lambda }_l}$}\right.\right)}-\frac{S}{\left[1-\left(S-1\right)\cdot u\right]\cdot L}\right]$$

(44)

From the above deduction, we can see the optimal performance is achieved by the two above algorithms. By comparing (43) and (44), it can be observed that the difference is in the dashed block. The pseudo-$ℒ$₀-norm optimization algorithm is more sensitive to small coefficients. For the low SNR targets, i.e., ${\xi }_{{\Lambda }_l}=1/\left({\sigma }_{{\Lambda }_l}+\iota \right)\succ 1$, our proposed algorithm performs better and it can more accurately recover low SNR targets. For example, in the one-dimensional simulation, our algorithm recovers the weak target, and the amplitude of the target is more accurately recovered in Figure 6.

Sub-Nyquist SAR imaging includes three steps: range compression, range cell migration correction (RCMC), and azimuth compression. After range compression and RCMC, a pseudo-$ℒ$₀-norm optimization algorithm is employed to achieve the azimuth compression. Simulation experiments in Figure 5 and Figure 6, and data experiments based on a real TerraSAR-X image in Figure 7, demonstrate that the algorithm can achieve smaller recovered errors and remove false targets compared to that of the $ℒ$₁-norm optimization algorithm. Additionally, considering that CS algorithms themselves have a certain capacity to remove phase errors, an error correction method integrated with a pseudo-$ℒ$₀-norm optimization algorithm eliminates the influence of phase errors and removes defocusing. Different from autofocusing technology as a post-processing method in a traditional SAR system, this method includes two steps: scene reconstruction based on a pseudo-$ℒ$₀-norm optimization algorithm and error estimation by minimizing the least-square target function. These two steps are successively iterative and the recovered matrix is updated to the estimated error. Since the exact expression of the error is not fully formulated, we propose the error matrix with a random phase, whatever the expression for the error. The simulations in the Figure 8 and Figure 9 explain the advantages of our proposed algorithm.

The potential applications of the pseudo-$ℒ$₀-norm optimization algorithm can recover more low SNR targets, e.g., the sea–land interface in the above experiment, so that it is not limited to high SNR targets, e.g., the boats on the sea. This algorithm can also be applied to the traditional HRWS system, e.g., the azimuthal multi-channel SAR and MIMO SAR, to lower the amount of data to relieve the pressure on data storage.

5. Conclusions

In this paper, we propose a signal processing algorithm for an innovative single-channel HRWS system called sub-Nyquist SAR based on a pseudo-random space–time modulation. To further improve sub-Nyquist SAR imaging performance, a pseudo-$ℒ$₀-norm optimization algorithm is proposed. Firstly, the Bayesian estimation explains how to take a more democratic approach with data prior information to acquire nonzero variables, so that it can more accurately recover low SNR targets and remove false targets compared with the prevalent $ℒ$₁-norm optimization algorithm. Then we present an error correction method integrated with a pseudo-$ℒ$₀-norm optimization algorithm to eliminate the effect of phase errors and achieve autofocusing. Finally, the simulated experiments demonstrate the effectiveness of the proposed algorithm.