Bias Analysis and Correction for Ill-Posed Inversion Problem with Sparsity Regularization Based on L₁ Norm for Azimuth Super-Resolution of Radar Forward-Looking Imaging

Abstract

The sparsity regularization based on the L₁ norm can significantly stabilize the solution of the ill-posed sparsity inversion problem, e.g., azimuth super-resolution of radar forward-looking imaging, which can effectively suppress the noise and reduce the blurry effect of the convolution kernel. In practice, the total variation (TV) and TV-sparsity (TVS) regularizations based on the L₁ norm are widely adopted in solving the ill-posed problem. Generally, however, the existence of bias is ignored, which is incomplete in theory. This paper places emphasis on analyzing the partially biased property of the L₁ norm. On this basis, we derive the partially bias-corrected solution of TVS and TV, which improves the rigor of the theory. Lastly, two groups of experimental results reflect that the proposed methods with partial bias correction can preserve higher quality than those without bias correction. The proposed methods not only distinguish the adjacent targets, suppress the noise, and preserve the shape and size of targets in visual terms. Its improvement of Peak Signal-to-Noise Ratio, Structure-Similarity, and Sum-Squared-Errors assessment indexes are overall 2.15%, 1.88%, and 4.14%, respectively. As such, we confirm the theoretical rigor and practical feasibility of the partially bias-corrected solution with sparsity regularization based on the L₁ norm.

1. Introduction

Ill-posed models widely exist in fields of geodesy and remote sensing, such as GNSS (Global Navigation Satellite System) fast ambiguity resolution [1,2,3], gravity field determination from the satellite missions of Gravity Recovery And Climate Experiment (GRACE) [4,5,6,7,8,9,10], geometric correction of satellite imagery based on rational function model (RFM) [11,12,13], azimuth super-resolution imaging of Real Aperture Radar (RAR) [13,14,15,16,17,18,19,20,21,22,23,24,25], signal restoration for vision sensing [26,27,28], full waveform Lidar data deconvolution [29,30,31] and optical remote sensing deblurring [32,33,34]. As is known, the solutions of a well-posed model should satisfy three standards: i.e., existence, uniqueness, and stability [35]. For any ill-posed observational equations, the coefficient matrix is full-rank, but its condition number is significantly great. In consequence, very small observation errors can also corrupt the estimated parameters considerably [36,37]. In other words, the solution of an ill-posed model violates the stability of necessary conditions of a well-posed solution, which is existent and unique but unstable. Therefore, regularization methods are proposed to stabilize the ill-posed solutions.

Tikhonov regularization has been successfully proposed to solve the ill-posed problems, which is based on the L₂-norm criterion that the residuals and parameters sum of squares are minimized [38,39]. Tikhonov regularization utilizes the regularization parameter to reduce the influence of the small singular values of the coefficient matrix and further decrease the condition number to get a better regularization solution. But it also introduces biases to the estimated solution. Considering the impact of small singular values, Xu [40] and Hansen [41] directly truncated the small singular values and implemented pseudo inversion to obtain a stable solution, which is the famous Truncated Singular Value Decomposition (TSVD) regularization. Although TSVD can get fine and stable solutions for ill-posed models, it still is biased due to the information loss of the coefficient matrix. To this end, Xu et al. [42], Xu [43], Shen et al. [37], Chen et al. [44], and Ji et al. [36] corrected the biases of the biased parameter estimates with fully-, partially- and adaptively- bias correction strategies. These methods of Tikhonov regularization based on the L₂ norm can depress the noise by smoothing the processed results. Thereby, when latent sparsity of unknown parameters exists, the Tikhonov regularization based on the L₂ norm will be unsuitable. For example, in the imaging region, the interest targets are sparse in ship monitoring or airport surveillance compared with the entire imaging area [16].

Therefore, in order to reconstruct these sparse characteristics, the sparsity regularization terms based on the L₁ norm were applied to estimate sparse parameter solutions. Azadbakht et al. [31] proposed a sparsity-constrained regularization approach for the deconvolution of the returned Lidar waveform and successfully restored the target cross-section. Zhang et al. [23] used the method solved by splitting the Bergman algorithm to get strong-point targets. A sparse denoising-based super-resolution method (SDBSM) was proposed in [18] to avoid losing the shape characteristics treated as strong-point targets. Furthermore, Tuo et al. [45] developed an iterative reweighted least-squares method for the sparsity-constrained regularization method to accelerate the computation speed. Zhang et al. [46] designed a total variation (TV) [47] super-resolution imaging method based on the TSVD strategy. Shortly after, a fast sparse-TSVD version was put forward by Tuo et al. [20] to improve its efficiency. Zhang et al. [22] introduced the TV method to restore the contour information of the target smoother and the L₁ norm to preserve the sparsity to improve the resolution of radar forward-looking imaging, which was called the TV-sparse (TVS) method. Similarly, Quan et al. [19] developed an improved quasi-Newton iteration method based on the Graphics Processing Unit Platform to raise the computation efficiency of sparse reconstruction. Due to the information loss of the coefficient matrix, the restored results of TV are slightly inferior to those of TVS. The Gohberg-Senmencul representation-based fast TV method was proposed by Zhang et al. [17] to improve the computation efficiency. Recently, Huo et al. [48] conceptualized a balanced Tikhonov and TV regularization approach to retain the sparsity of the candidate solution and the smoothness of continuous contours in the regularization solution. Although these methods based on the L₁ norm with different algorithms can get fine results, which can preserve the sparsity of the unknown parameters and its smooth shape of continuous contours, the bias of the regularization solution based on the L₁ norm cannot be analyzed nor mitigated. Obviously, it is theoretically insufficient and incomplete for the analytical solution of the L₁ norm.

To this end, the main purpose of the paper addresses three goals. They are: (1) The bias analysis of the sparsity regularization method based on the L₁ norm is implemented to place emphasis on its theoretical importance and make up for its incompleteness. (2) The comparison with the bias of the sparsity regularization method based on the L₂ norm is realized to explicate its advantage with partial bias. The consistency with the conclusion of partially bias correction regularization methods proposed by Shen et al. [37] is analyzed, and the theoretic reason for its advantage will be mined. (3) A piecewise partially bias-corrected regularization solution based on the L₁ norm is deduced, and the corresponding iterative algorithms based on the two typical TV and TVS models will be designed to achieve a more accurate solution.

The remainder of this paper is organized as follows. Section 2 contains a brief introduction to the azimuth echo convolution model of radar forward-looking imaging. In Section 3, an analysis and a comparison with the regularization solution based on the L₁ norm and L₂ norm are carried out, with their differences being expressed and the nature of partial bias being revealed. In Section 4, taking two typical TV and TVS methods to alleviate the ill-posed problem as examples, we derive their partial bias correction solutions with a piecewise form. The corresponding iterative flowcharts are designed. In Section 5, two experimental examples have been designed to demonstrate the performance of the proposed methods, namely, 1-D point target and 2-D area simulation for azimuth super-resolution of radar forward-looking imaging. A discussion on the superiority of the proposed methods with bias correction is implemented in Section 6. In final, conclusions and remarks are summarized in Section 7.

2. Azimuth Echo Convolution Model of Radar Forward-Looking Imaging

The azimuth echo signal convolution model of radar forward-looking imaging is illustrated in Figure 1, according to which the imaging platform travels along the Y-direction with a velocity of v and at the height of H. R₀ is the initial slant range between the radar and the target P, with the initial azimuth angle being θ₀. φ₀ is the initial pitching angle. As the platform travels with a time interval t, the slant range changes into R, the pitching angle turns to φ, and the azimuth angle becomes θ.

According to the trigonometric relation, the slant range history R(t) satisfies

$$\mathit{R}\left(\mathit{t}\right)=\sqrt{{\mathit{R}}_0^2+{\left(\mathit{vt}\right)}^2-2{\mathit{R}}_0\mathit{vtcos}{\theta }_0\mathit{cos}{\mathit{\phi }}_0}\approx {\mathit{R}}_0-\mathit{vtcos}{\mathit{\theta }}_0\mathit{cos}{\mathit{\phi }}_0$$

(1)

Considering both the range resolution and working distance, the radar transmits linear frequency-modulated signal is expressed as

$$r\left(\tau \right)=\mathit{rect}\left(\frac{\tau }{{\mathit{T}}_{\mathit{p}}}\right)\mathit{exp}\left(\mathit{j}2\pi {\mathit{f}}_{\mathit{c}}\right)\mathit{exp}\left(\mathit{j}\pi \mu {\tau }^2\right)$$

(2)

where τ is the fast times in range direction, the carrier frequency of the transmitted signal is denoted as f_c, T_p is the signal duration, and μ stands for the chirp rate; rect(·) is a rectangle window. After platform traveling and antenna scanning, the received signal of the target P can be evolved as follows,

$$\mathit{r}\left(\tau ,\mathit{t}\right)={\sigma }_{\mathit{p}}\mathit{h}\left(\mathit{t}-{\mathit{t}}_0\right)\mathit{rect}\left(\frac{\tau -{\tau }_{\mathit{d}}}{{\mathit{T}}_{\mathit{p}}}\right)\mathit{exp}\left(\mathit{j}2\pi {\mathit{f}}_{\mathit{c}}\left(\frac{\tau -{\tau }_{\mathit{d}}}{{\mathit{T}}_{\mathit{p}}}\right)\right)\mathit{exp}\left(\mathit{j}\pi \mu {\left(\frac{\tau -{\tau }_{\mathit{d}}}{{\mathit{T}}_{\mathit{p}}}\right)}^2\right)$$

(3)

where h(·) represents the modulation effect of the antenna pattern, σ_p is the target scattering coefficient, ${\tau }_{\mathit{d}}=2\mathit{R}\left(\mathit{t}\right)/\mathit{c}$ is time delay, and t₀ is the moment when the beam center scans P. By matched filtering the widely utilized technology for high range resolution, the received signal becomes

$$\mathit{r}\left(\tau ,\mathit{t}\right)={\sigma }_{\mathit{p}}\mathit{h}\left(\mathit{t}-{\mathit{t}}_0\right)\mathit{sin}\mathit{c}\left({\mathit{B}}_{\mathit{r}}\left(\tau -{\tau }_{\mathit{d}}\right)\right)\times \mathit{exp}\left(-\mathit{j}2\pi {\mathit{f}}_{\mathit{c}}{\tau }_{\mathit{d}}\right)$$

(4)

Radar scanning imaging is accompanied by platform movement, which results in the echoes aliasing in different units. Therefore, its influence has to be removed [16,46,48]. Furthermore, the received echo can be modeled in the following convolution form,

$$\mathit{r}=\mathit{a}\ast \mathit{u}+\mathit{\epsilon }$$

(5)

where r represents the discrete azimuth echo signal vector, a is the convolution kernel standing for antenna pattern signal, “$\ast$” is the convolution symbol, and ε is the random error.

3. Analysis and Comparison between L₂ Norm and L₁ Norm

In order to show the bias of the regularization solution based on the L1 norm intuitively, we derive the formulae of the bias and analyze and compare it with that based on the L₂ norm.

For the linear discrete ill-posed models, its expression can be written by,

$$\mathit{r}=\mathbf{A}\mathit{u}+\mathit{\epsilon }$$

(6)

where r is an m-vector of observations, A is an m × n (m > n) deterministic coefficient matrix, which is assumed as full column rank but with a large condition number, u is the n-vector of unknown parameters, and ε is the random error vector with zero mean and variance-covariance matrix ${\mathit{\sigma }}_0^2{\mathbf{W}}^{-1}$. Here ${\mathit{\sigma }}_0^2$ is the variance of unit weight, and W is a weight matrix, which is set as the identity matrix in this paper. According to the Tikhonov regularization criterion, the objective function is expressed as,

$$\widehat{\mathit{u}}=\underset{\widehat{u}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\alpha }{\Vert \mathit{u}\Vert }_2^2$$

(7)

where α is the positive regularization parameter and ${\Vert ·\Vert }_2^2$ is the L₂ norm. For the solution of the Tikhonov regularization method based on the L₂ norm, namely ${\widehat{\mathit{u}}}_{\mathit{T}}={\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\alpha }\mathbf{I}\right)}^{-1}{\mathbf{A}}^{\mathit{T}}\mathit{r}$, applying the mathematical expectation and matrix inversion formula to it, we get $\mathrm{E}\left({\widehat{\mathit{u}}}_{\mathit{T}}\right)=\overline{\mathit{u}}-\mathit{\alpha }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\alpha }\mathbf{I}\right)}^{-1}\overline{\mathit{u}}$ where E(·) is the operator of mathematical expectation and $\overline{\mathit{u}}$ is the true value of the unknown parameter. Thus, its bias can be obtained according to [37,42],

$$\mathit{bias}\left({\widehat{\mathit{u}}}_{\mathit{T}}\right)=-\mathit{\alpha }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\alpha }I\right)}^{-1}\overline{\mathit{u}}$$

(8)

Since the bias vector can be formally computed with Equation (8), we can formally remove the biases from the regularization solution and obtain the bias-corrected regularized estimate as follows:

$${\mathit{u}}_{\mathit{c}}={\widehat{\mathit{u}}}_{\mathit{T}}-\mathit{bias}\left({\overline{\mathit{u}}}_{\mathit{T}}\right)={\widehat{\mathit{u}}}_{\mathit{T}}+\mathit{\alpha }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\alpha }\mathbf{I}\right)}^{-1}\overline{\mathit{u}}$$

(9)

But the L₂ norm representing the smoothing characters cannot solve the ill-posed problem with sparsity regularization, and in particular, of deconvolution problem. According to Tuo et al. [20], it can be recast as Equation (6), where A is a corresponding circulant Toeplitz matrix from a. In general, the L₀ norm is selected to represent the sparsity of the parameters, but it is a non-convex problem which is hard to solve. Therefore, the L₁ norm is utilized to replace it with the L₀ norm for a balance between sparsity and smoothness. For the ill-posed problem with sparsity regularization, the objective function (7) can be recast as

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\beta }{\Vert \mathit{u}\Vert }_1$$

(10)

where β is a positive regularization parameter. A splitting Bergman algorithm (SBA) can be used to estimate the parameter u [49,50]. An intermediate variable w = u is introduced, and Equation (10) can be rewritten as

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\beta }{\Vert \mathit{w}\Vert }_1,\mathit{s}.\mathit{t}.\hspace{0.17em}\mathit{w}=\mathit{u}$$

According to the SBA, the above-mentioned formula can be further evolved with the following iterative form,

$$\{\begin{array}{l}\left({\widehat{\mathit{u}}}^{\mathit{k}+1},{\mathit{w}}^{\mathit{k}+1}\right)=\underset{\mathit{u},\mathit{d}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\beta }{\Vert \mathit{w}\Vert }_1+\mathit{\gamma } {\Vert \mathit{w}-\mathit{u}-{\mathit{b}}^{\mathit{k}}\Vert }_2^2\\ {\mathit{b}}^{\mathit{k}+1}={\mathit{b}}^{\mathit{k}}+\left({\widehat{\mathit{u}}}^{\mathit{k}+1}-{\mathit{w}}^{\mathit{k}+1}\right)\end{array}$$

(11)

where γ is the penalty parameter with a large number, and b is the gradient of w. The superscript k stands for the kth iteration. The solutions of iterative strategy can be updated with the following formulae,

$$\{\begin{array}{l}{\widehat{\mathit{u}}}^{\mathit{k}+1}={\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\left({\mathbf{A}}^{\mathit{T}}\mathit{r}+\mathit{\gamma }\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}\right)\right)\\ {\mathit{w}}^{\mathit{k}+1}=\mathit{sign}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}\right)\max \left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}-\mathit{\beta }/2\mathit{\gamma },0\right)\\ {\mathit{b}}^{\mathit{k}+1}={\mathit{b}}^{\mathit{k}}+{\widehat{\mathit{u}}}^{\mathit{k}+1}-{\mathit{w}}^{\mathit{k}+1}\end{array}$$

(12)

where sign(·) stands for the signum function. The first formula of Equation (12) is similar to the Tikhonov regularization solution and with the extra term ${\mathit{w}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}$, which is dependent on SBA. Thereby, we derive the equivalent formula according to the matrix inversion formula ${\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}$ as follows,

$$\begin{array}{ll}\widehat{\mathit{u}}& =\left[{\left({\mathbf{A}}^{\mathbf{T}}\mathbf{A}\right)}^{-1}-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}\right)}^{-1}\right]\left({\mathbf{A}}^{\mathit{T}}\mathit{r}+\mathit{\gamma }\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}\right)\right)\\ & ={\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathit{T}}\mathit{r}-\mathit{\gamma }\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right){\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathit{T}}\mathit{r}+\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}\right)\\ & =\overline{\mathit{u}}-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\left(\overline{\mathit{u}}-\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}\right)\right)\end{array}$$

(13)

For convenience to express, the superscript k is omitted. According to the definition of expectation, the bias $\widehat{\mathit{u}}$ can be expressed by,

$$\mathit{bias}\left(\widehat{\mathit{u}}\right)=-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\left(\overline{\mathit{u}}-\left(\mathit{w}-\mathit{b}\right)\right)$$

(14)

When γ is equal to α, through subtracting Equation (8) and Equation (14), one can obtain the difference Δb between the L₁ norm and L₂ norm as,

$$\Delta \mathit{b}=-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\left(\mathit{w}-\mathit{b}\right)$$

(15)

Substituting the intermediate variable solution w in Equation (12) into Δb, we can achieve the following piecewise formulation as

$$\Delta \mathit{b}=\{\begin{array}{l}-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\widehat{\mathit{u}},\widehat{\mathit{u}}-\mathit{b}>\mathit{\beta }/2\mathit{\gamma }\\ \mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\mathit{b} ,\widehat{\mathit{u}}-\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }\end{array}$$

(16)

Referring to the definition of SBA, the intermediate variable b is a small value. Therefore, the above-mentioned Δb can be approximately recast as,

$$\Delta \mathit{b}\approx \{\begin{array}{l}-\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }\mathbf{I}\right)}^{-1}\widehat{\mathit{u}},\widehat{\mathit{u}}-\mathit{b}>\mathit{\beta }/2\mathit{\gamma }\\ 0 ,\widehat{\mathit{u}}-\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }\end{array}$$

(17)

The equation infers that the bias existing in the L₁ norm and L₂ norm are closely equal when $\widehat{\mathit{u}}-\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }$. In other words, the solution of the L₂ norm needs to correct the part of this condition $\widehat{\mathit{u}}-\mathit{b}>\mathit{\beta }/2\mathit{\gamma }$, while that of the L₁ norm need not do. Furthermore, according to the definition of iterative shrinkage threshold algorithms about w, $\overline{\mathit{u}}-\left(\mathit{w}-\mathit{b}\right)$ bias($\widehat{\mathit{u}}$) can be rewritten as

$$\overline{\mathit{u}}-\left(\mathit{w}-\mathit{b}\right)=\{\begin{array}{l}\overline{\mathit{u}}-\widehat{\mathit{u}}+\mathit{\beta }/2\mathit{\gamma },\overline{\mathit{u}}+\mathit{b}>\mathit{\beta }/2\mathit{\gamma }\\ \widehat{\mathit{u}}+\mathit{b}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\overline{\mathit{u}}+\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }\end{array}$$

(18)

With the iterations employed deeply, the true value $\overline{\mathit{u}}$ is gradually approached by the estimated $\widehat{\mathit{u}}$, i.e., Equation (18) can be updated with the form of limit as follows,

$$\underset{\widehat{\mathit{u}}\to \overline{\mathit{u}}}{\lim }\left(\overline{\mathit{u}}-\left(\mathit{w}-\mathit{b}\right)\right)=\{\begin{array}{l}\mathit{\beta }/2\mathit{\gamma },\overline{\mathit{u}}+\mathit{b}>\mathit{\beta }/2\mathit{\gamma }\\ \widehat{\mathit{u}}+\mathit{b}\hspace{0.17em},\overline{\mathit{u}}+\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }\end{array}$$

(19)

The positive penalty parameter γ is a large value and grows with a scale greater than 1 as its iteration implements, then β is a small positive value. Thus, the value of β/γ is close to zero. Thereby, Equation (19) can be approximately recast as,

$$\underset{\widehat{\mathit{u}}\to \overline{\mathit{u}}}{\lim }\left(\overline{\mathit{u}}-\left(\mathit{w}-\mathit{b}\right)\right)\approx \{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\widehat{\mathit{u}}+\mathit{b}>\mathit{\beta }/2\mathit{\gamma }\\ \widehat{\mathit{u}}+\mathit{b}\hspace{0.17em},\widehat{\mathit{u}}+\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }\end{array}$$

(20)

From it, $\widehat{\mathit{u}}$ of Equation (13) based on the L₁ norm is a partially bias-corrected regularization solution, obviously. Namely, the bias existing in the solution $\widehat{\mathit{u}}$ based on the L₁ norm is partially unbiased when $\widehat{\mathit{u}}-\mathit{b}\le \mathit{\beta }/2\mathit{\gamma }$. Shen et al. [37] reported that the partially bias-corrected regularization solution is superior to the fully bias-corrected regularization solution (7) but depends on the inequality with regard to the singular values of coefficient matrix A and the regularization parameter α. It suggests that extra work to fully bias-correct is employed for the Tikhonov regularization solution based on the L₂ norm, while the regularization solution based on the L₁ norm for partial bias correction originates from its own characteristics contained in Equation (20). In other words, from the frequentists’ point of view, the low-frequency part of singular values of the coefficient matrix is not corrected, while the rest part is corrected [36].

4. Bias Correction for TV-Sparse and TV Model

In ill-posed sparsity inversion problems for latent sparsity parameters, like deconvolution, the TVS and TV methods are usually used to relax it for satisfactory results. But their biases are typically ignored. Therefore, we give the following deductions to correct their biases.

4.1. Deduction of TVS Model with Bias Correction

The TVS method proposed by Zhang et al. [22] is based on TV regularization term utilizing the sparsity property of the signal ${\Vert \mathit{u}\Vert }_1$ and its gradient ${\Vert \nabla \mathit{u}\Vert }_1$, whose objective function can be evolved from Equation (10) as follows,

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\alpha }{\Vert \nabla \mathit{u}\Vert }_1+\mathit{\beta }{\Vert \mathit{u}\Vert }_1$$

(21)

According to SBA, the objective function (21) can be recast by introducing the intermediate variables v = ∇u and w = u,

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\arg \min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\alpha }{\Vert \mathit{v}\Vert }_1+\mathit{\beta }{\Vert \mathit{w}\Vert }_1+{\mathit{\gamma }}_1{\Vert \mathit{v}-\nabla \mathit{u}-{\mathit{b}}_1\Vert }_2^2+{\mathit{\gamma }}_2{\Vert \mathit{w}-\mathit{u}-{\mathit{b}}_2\Vert }_2^2$$

(22)

where b₁ and b₂ are the gradient of ∇u and u, γ₁ and γ₂ are the penalty parameters with large values, and α and β are the positive regularization parameters. The corresponding solutions with an iterative strategy can be derived as,

$$\{\begin{array}{l}{\mathit{u}}^{\mathit{k}+1}={\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla +{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left({\mathbf{A}}^{\mathit{T}}\mathit{r}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\left({\mathit{v}}^{\mathit{k}}-{\mathit{b}}_1^{\mathit{k}}\right)+{\mathit{\gamma }}_2\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}_2^{\mathit{k}}\right)\right)\\ {\mathit{v}}^{\mathit{k}+1}=\mathit{sign}\left(\nabla {\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}}\right)\max \left(\nabla {\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}}-\mathit{\alpha }/2{\mathit{\gamma }}_1,0\right)\\ {\mathit{w}}^{\mathit{k}+1}=\mathit{sign}\left({\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}\right)\max \left({\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}}-\mathit{\beta }/2{\mathit{\gamma }}_2,0\right)\\ {\mathit{b}}_1^{\mathit{k}+1}={\mathit{b}}_2^{\mathit{k}}+\nabla {\mathit{u}}^{\mathit{k}+1}-{\mathit{v}}^{\mathit{k}+1}\\ {\mathit{b}}_2^{\mathit{k}+1}={\mathit{b}}_1^{\mathit{k}}+{\mathit{u}}^{\mathit{k}+1}-{\mathit{w}}^{\mathit{k}+1}\end{array}$$

(23)

Referring to the analyses in Section 3 and Equations (13)–(15), the formula of the estimated solution u for TVS can be rewritten as follows,

$$\begin{array}{ll}{\mathit{u}}^{\mathit{k}+1}=& \overline{\mathit{u}}-{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla +{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla \overline{\mathit{u}}-\left({\mathit{v}}^{\mathit{k}}-{\mathit{b}}_1^{\mathit{k}}\right)\right)\\ & -{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla +{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left(\overline{\mathit{u}}-\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}_2^{\mathit{k}}\right)\right)\end{array}$$

(24)

The bias can consist of two parts, including $\mathit{b}\mathit{i}\mathbf{a}{\mathit{s}}_1=-{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla +{\mathit{\gamma }}_2\right)}^{-1}·$${\nabla }^{\mathit{T}}\left(\nabla \overline{\mathit{u}}-\left({\mathit{v}}^{\mathit{k}}-{\mathit{b}}_1^{\mathit{k}}\right)\right)$ and $\mathit{b}i\mathbf{a}\mathit{s}{ }_2=-{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla +{\mathit{\gamma }}_2\right)}^{-1}\left(\overline{\mathit{u}}-\left({\mathit{w}}^{\mathit{k}}-{\mathit{b}}_2^{\mathit{k}}\right)\right)$. For the convenience of representation, $ℝ$ stands for the real number space. We define a set ${\mathrm{\Phi }}_{2,1}\left(\mathit{x}\right)=\left\{\mathit{x}|\mathit{x}>\mathit{\beta }/2{\mathit{\gamma }}_2\right\}$,${\mathrm{\Phi }}_{2,2}\left(\mathit{x}\right)=\left\{\mathit{x}|\mathit{x}\le \mathit{\beta }/2{\mathit{\gamma }}_2\right\}$, and ${\mathrm{\Phi }}_{2.1}\cup {\mathrm{\Phi }}_{2.2}=ℝ$. Similarly, the first term bias₁ can be reconstructed as,

$${\mathit{bias}}_1\approx \{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},{\mathrm{\Phi }}_{1,1}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\\ -{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em},{\mathrm{\Phi }}_{1,2}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\end{array}$$

(25)

Set ${\mathrm{\Phi }}_{1,1}\left(\mathit{x}\right)=\left\{\mathit{x}|\mathit{x}>\mathit{\alpha }/2{\mathit{\gamma }}_1\right\}$, ${\mathrm{\Phi }}_{1,2}\left(\mathit{x}\right)=\left\{\mathit{x}|\mathit{x}\le \mathit{\alpha }/2{\mathit{\gamma }}_1\right\}$, and ${\mathrm{\Phi }}_{1.1}\cup {\mathrm{\Phi }}_{1.2}=ℝ$. Due to the uncertainty of the relationship between sets ${\mathrm{\Phi }}_{1,1}$ and ${\mathrm{\Phi }}_{2,1}$, a necessary discussion can be carried out as follows. By comparison with Equations (20) and (25), the inclusion relation of sets ${\mathrm{\Phi }}_{1,1}$ and ${\mathrm{\Phi }}_{2,1}$ is dependent on the values of α/γ₁ and β/γ₂, namely, α/γ₁ < β/γ₂, α/γ₁ = β/γ₂, and α/γ₁ > β/γ₂. Concretely,

(1)

when $\mathit{\alpha }/{\mathit{\gamma }}_1<\mathit{\beta }/{\mathit{\gamma }}_2$, i.e., ${\mathrm{\Phi }}_{2,1}\subseteq {\mathrm{\Phi }}_{1,1}$ and ${\mathrm{\Phi }}_{1,2}\subseteq {\mathrm{\Phi }}_{2,2}$. The partially bias-corrected solution of Equation (25) can be obtained by,

$${\mathit{u}}_{\mathit{T}\mathit{V}\mathit{S}\mathit{B}\mathit{C}}^{\mathit{k}+1}=\{\begin{array}{l}{\widehat{\mathit{u}}}^{\mathit{k}+1},{\mathrm{\Phi }}_{2,1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right)\\ {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right),{\mathrm{\Phi }}_{2,1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right)\cap {\mathrm{\Phi }}_{1,1}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\\ {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)+{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right),{\mathrm{\Phi }}_{1,2}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\end{array}$$

(26)

(2)

when $\mathit{\alpha }/{\mathit{\gamma }}_1=\mathit{\beta }/{\mathit{\gamma }}_2$, i.e., ${\mathrm{\Phi }}_{2,1}={\mathrm{\Phi }}_{1,1}$ and ${\mathrm{\Phi }}_{1,2}={\mathrm{\Phi }}_{2,2}$. The partially bias-corrected solution of Equation (25) can be simplified as,

$${\mathit{u}}_{\mathit{T}\mathit{V}\mathit{S}\mathit{B}\mathit{C}}^{\mathit{k}+1}=\{\begin{array}{l}{\widehat{\mathit{u}}}^{\mathit{k}+1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},{\mathrm{\Phi }}_{2,1}\left(\widehat{\mathit{u}}+{\mathit{b}}_2^{\mathit{k}+1}\right)\hspace{0.17em}\hspace{0.17em}\mathit{o}\mathit{r}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{\Phi }}_{1,1}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)+{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},{\mathrm{\Phi }}_{1,2}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em} \mathit{o}\mathit{r}\hspace{0.17em}\hspace{0.17em}{\mathrm{\Phi }}_{2,2}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right)\end{array}$$

(27)

(3)

when $\mathit{\alpha }/{\mathit{\gamma }}_1>\mathit{\beta }/{\mathit{\gamma }}_2$, i.e., ${\mathrm{\Phi }}_{1,1}\subseteq {\mathrm{\Phi }}_{2,1}$ and ${\mathrm{\Phi }}_{2,2}\subseteq {\mathrm{\Phi }}_{1,2}$. The partially bias-corrected solution of Equation (25) can be refined as,

$${\mathit{u}}_{\mathit{T}\mathit{V}\mathit{S}\mathit{B}\mathit{C}}^{\mathit{k}+1}=\{\begin{array}{l}{\widehat{\mathit{u}}}^{\mathit{k}+1},{\mathrm{\Phi }}_{1,1}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em}\\ {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right),{\mathrm{\Phi }}_{2,1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right)\hspace{0.17em}\cap {\mathrm{\Phi }}_{1,1}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em}\\ {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{\gamma }}_1{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_1{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}\right)+{\mathit{\gamma }}_2{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+{\mathit{\gamma }}_2\mathbf{I}\right)}^{-1}\left({\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_2^{\mathit{k}+1}\right),{\mathrm{\Phi }}_{1,2}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}_1^{\mathit{k}+1}\right)\hspace{0.17em}\end{array}$$

(28)

Obviously, the solution of the TVS model based on the L₁ norm is partially biased, which is a key point superior to those based on the L₂ norm. The implemented flowchart of the proposed TVS with bias correction (TVSBC) can be found in Figure 2.

4.2. Extension of TV Model with Bias Correction

Considering the TV model proposed by Zhang et al. [17], we find that the TV model is simpler than the TV-sparse model without $\mathit{\beta }{\Vert \mathit{u}\Vert }_1$, the expression of which is displayed as follows,

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\alpha }{\Vert \nabla \mathit{u}\Vert }_1$$

(29)

By employing the SBA to calculate the solution, the intermediate variable vectors v and b can be introduced. The corresponding objective function can be approximately transformed as,

$$\widehat{\mathit{u}}=\underset{\mathit{u}}{\min }{\Vert \mathbf{A}\mathit{u}-\mathit{r}\Vert }_2^2+\mathit{\alpha }{\Vert \mathit{v}\Vert }_1+\mathit{\gamma }{\Vert \nabla \mathit{u}-\left(\mathit{v}-\mathit{b}\right)\Vert }_2^2$$

(30)

And the target parameter and variable vectors can be estimated iteratively with the following formulae,

$$\{\begin{array}{l}{\widehat{\mathit{u}}}^{\mathit{k}+1}={\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }{\nabla }^{\mathit{T}}\nabla \right)}^{-1}\left({\mathbf{A}}^{\mathit{T}}\mathit{r}+\mathit{\gamma }{\nabla }^{\mathit{T}}\left({\mathit{v}}^{\mathit{k}}-{\mathit{b}}^{\mathit{k}}\right)\right)\\ {\mathit{v}}^{\mathit{k}+1}=\mathit{sign}\left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}\right)\max \left(\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}}-\mathit{\alpha }/2\mathit{\gamma },0\right)\\ {\mathit{b}}^{\mathit{k}+1}={\mathit{b}}^{\mathit{k}}+\nabla {\widehat{\mathit{u}}}^{\mathit{k}+1}-{\mathit{v}}^{\mathit{k}+1}\end{array}$$

(31)

As we know, the bias of Equation (31) can be derived easily,

$$\mathit{bias}\left(\widehat{\mathit{u}}\right)\approx \{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\nabla {\mathit{u}}^{\mathit{k}}+{\mathit{b}}^{\mathit{k}}>\mathit{\alpha }/2\mathit{\gamma }\\ -\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla \overline{\mathit{u}}+{\mathit{b}}^{\mathit{k}}\right)\hspace{0.17em},\nabla {\mathit{u}}^{\mathit{k}}+{\mathit{b}}^{\mathit{k}}\le \mathit{\alpha }/2\mathit{\gamma }\end{array}$$

(32)

And the bias-corrected solution of u_TVBC can be calculated according to Equation (32),

$${\mathit{u}}_{\mathit{T}\mathit{V}\mathit{B}\mathit{C}}^{\mathit{k}+1}\approx \{\begin{array}{l}{\mathit{u}}^{\mathit{k}+1}\hspace{0.17em},\nabla {\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}+1}>\mathit{\alpha }/2\mathit{\gamma }\\ {\mathit{u}}^{\mathit{k}+1}+\mathit{\gamma }{\left({\mathbf{A}}^{\mathit{T}}\mathbf{A}+\mathit{\gamma }{\nabla }^{\mathit{T}}\nabla \right)}^{-1}{\nabla }^{\mathit{T}}\left(\nabla {\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}+1}\right)\hspace{0.17em},\nabla {\mathit{u}}^{\mathit{k}+1}+{\mathit{b}}^{\mathit{k}+1}\le \mathit{\alpha }/2\mathit{\gamma }\end{array}$$

(33)

Obviously, by comparison with Equations (26)~(28) and (33) from the form of formulae and judgment criterion, TVBC is simpler than TVSBC. And its iterative process is shorter than TVSBC, and the flowchart can be displayed in Figure 3.

5. Experiments and Results

In addition to theoretical derivation, the performance of the proposed TVSBC and TVBC methods will be verified by two numerical experiments designed in this section. Two aspects of performance, including the ability to distinguish adjacent targets and maintain the shape of the target, are considered. Several traditional deconvolution methods, including Blind deconvolution [51], Richardson–Lucy method [52], Regularized Filter (Tikhonov Regularization) [53], Wiener Filter [53], Truncated Singular Value Decomposition (TSVD) [20], Sparse Denoising-Based Super-Resolution Method (SDBSM) [18], TV [17] and TVS [22], are selected as competitors to the two proposed methods.

5.1. Evaluated Indexes

As we know, in deconvolution problems, the true parameters are impossibly known in practical applications, and only their estimates can be obtained by these methods. However, since the simulation results can be repeated and the true parameters are also known, the sum square errors (SSE) can be numerically calculated as follows,

$$\mathit{M}={\displaystyle \sum _{i=1}^{\mathit{N}}}{\left(\widehat{\mathit{x}}-\overline{\mathit{x}}\right)}^{\mathit{T}}\left(\widehat{\mathit{x}}-\overline{\mathit{x}}\right)$$

(34)

where M denotes the numerical SSE of the estimated parameter, N is the number of repeated experiments, $\widehat{\mathit{x}}$ and $\overline{\mathit{x}}$ stand for the deconvolution estimated and true parameters. The smaller values of M, the closer the parameter estimates are to the true values.

Additionally, the peak signal-to-noise ratio (PSNR) is used to assess the ability of the noise suppression ability, which is defined as

$$\mathrm{PSNR}=20\log \left({\mathbf{A}}_{\mathbf{s}}/{\mathbf{A}}_{\mathbf{n}}\right)$$

(35)

where A_s is the maximum amplitude of the target and A_n is the maximum amplitude of the noise. And then, the Structural Similarity (SSIM) [54] is selected to evaluate the capability to differentiate adjacent targets and the similarity between the restored and original signals, which can be obtained by the following formula,

$$\mathrm{SSIM}\left(\mathit{x},\mathit{y}\right)=\frac{\left(2{\mathit{\mu }}_{\mathit{x}}{\mathit{\mu }}_{\mathit{y}}+{\mathit{C}}_1\right)\left(2{\mathit{\sigma }}_{\mathit{x}\mathit{y}}+{\mathit{C}}_2\right)}{\left({\mathit{\mu }}_{\mathit{x}}^2+{\mathit{\mu }}_{\mathit{y}}^2+{\mathit{C}}_1\right)\left({\mathit{\sigma }}_{\mathit{x}}^2+{\mathit{\sigma }}_{\mathit{y}}^2+{\mathit{C}}_2\right)}$$

(36)

where x and y reference signal and test signal; ${\mathit{\mu }}_{\mathit{x}}$, ${\mathit{\mu }}_{\mathit{y}}$, ${\mathit{\sigma }}_{\mathit{x}}$, and ${\mathit{\sigma }}_{\mathit{y}}$ are average value and standard variance of x and y, respectively; ${\mathit{\sigma }}_{\mathit{x}\mathit{y}}$ is the covariance of x and y; C₁, C₂, and C₃ are the contrast values. The greater PSNR and SSIM, the better the restoration results.

5.2. Experiment1: 1-D Point Target Simulation

In order to assess the performance of the proposed TVBC and TVSBC, referring to Tuo et al. [20], we conduct a 1-D point simulation of azimuth super-resolution for radar forward-looking imaging. The 1-D simulation scene is illustrated in Figure 4a, whose parameters are listed in

Table 1. System parameters of point target simulation.

Parameters	Value	Units
Carrier frequency	10	GHz
Band width	75	MHz
Pulse interval	2 × 10−6	s
Beamwidth	2	°
Antenna scanning velocity	30	°/s
Scanning area	−5~5	°
Pulse repetition frequency	1500	Hz

. There are four-point targets whose centers width are −1.95°, −0.35°, −0.15° and 0.25°, and whose widths are 0.2°, 0.2°, 0.1°, and 0.1°, respectively. The simulated convolution kernel is illustrated in Figure 5, whose parameters can be found in Table 1. For a more suitable actual working environment, Gaussian white noise with different SNRs of the range from 1 dB to 25 dB with an interval of 2 dB. Three typical simulated beam echoes affected by the noise of SNRs = 5 dB, 13 dB, and 21 dB are demonstrated in Figure 4b,c. The obvious features are contaminated by the simulated convolution kernel and random white noise. For example, two peaks in Figure 4b and three peaks in Figure 4c,d are seriously disturbed by high-frequency noise, which makes the signal unrecognizable.

We set up 500 independent simulation experiments with different SNRs. The degraded signals can be restored by these deconvolution methods, whose results can be evaluated by the mean PSNR, SSIM, and SSE listed in

Table 2. Mean PSNR (dB), SSIM and SSE results with 500 simulations.

Methods	Mean PSNR (dB)	Mean SSIM	Mean SSE
Blind Deconvolution	17.211	0.277	0.541
Regularized Filter	10.304	0.156	1.328
Wiener Filter	17.310	0.499	0.531
Richardson–Lucy	15.195	0.698	0.684
TSVD	15.336	0.587	0.666
SDBSM	18.644	0.270	31.292
TV	18.510	0.801	0.469
TVBC	18.972	0.816	0.444
TVS	19.214	0.803	0.4309
TVSBC	19.244	0.809	0.4305

. From it, we can see that Blind Deconvolution, Wiener Filter, and SDBSM can have fine mean PSNR but terrible mean SSIM, and Regularized Filter based on L₂ norm proposed for the smooth signal has the lowest mean PSNR and SSIM due to a loss in the accurate noise level estimation and smooth the significant signal features. With the decreasing SNRs during simulation, the noise level increases, whose effect on these above-mentioned deconvolution methods has been magnified accordingly. It infers that these deconvolution methods are not robust against noise. In addition, the Richardson–Lucy method and TSVD method can keep a balance between mean PSNR and SSIM and with lesser mean SSE, which maintains a middling level compared with TV, TVS, TVBC, and TVSBC. Overall, TV, TVS, TVBC, and TVSBC provide excellent evaluation results that are evidenced by high mean PSNR and SSIM and small mean SSE. By comparing the TV/TVS without and with bias correction, we can find that the proposed bias correction methods can effectively improve the performance of the original TV and TVS. It seems that the proposed methods are not only rigorous in theory but also feasible in practice.

A group of restoration results in 500 simulations when SNR = 21 dB are illustrated in Figure 6. Due to the low noise level, isolated targets can be distinguished by all deconvolution methods. But Blind Deconvolution, Regularized Filter, Wiener Filter, TSVD, and SDBSM cannot distinguish the middle two adjacent targets. From Figure 6d, the Richardson–Lucy method can recover the middle two adjacent targets, but an extra fake peak accompanies the left isolated targets. In addition, as illustrated in Figure 6g–j, TV, TVS, TVBC, and TVSBC have excellent capability of distinguishment to these isolated and adjacent targets, and few significant fake targets are produced.

To better compare TV and TVS methods with and without bias correction, we display three groups of results in Figure 7 with high, middle, and low SNRs. With the increase in the simulated SNRs, these results of TV and TVS methods can restore better-qualified images, and these main targets can be distinguished better. Furthermore, these results after bias correction by TVBC and TVSBC can be higher restoration quality than those of TV and TVS methods. The intensity of these significant targets recovered from TVBC and TVSBC is higher than that of TV and TVS, as demonstrated in Figure 7 with red rectangles. Similarly, the restored contour information is illustrated in Figure 7 with gray rectangles, whose intensities of TVBC and TVSBC are lower than those of TV and TVS. It proves that bias correction can effectively improve the restored quality of the 1-D targets scene in visual.

Furthermore, for better quantitative comparison, the line charts of PSNR, SSIM, and SSE with the increase of simulation SNRs ranging from 1~25 dB with 2 dB interval are shown in Figure 8. As expected, the results of partial bias correction of TVBC and TVSBC are improved. Among them, the PSNRs and SSIMs of TVBC and TVSBC are greater than those of TV and TVS, while the SSEs of TVBC and TVSBC are less than those of TV and TVS. The mean improvement ratio of PSNR, SSIM, and SSE of TV and TVS relative to TVBC and TVSBC are 3.00% and 1.00%, 1.94% and 0.70%, 6.00%, and 2.08%, respectively.

In addition, the statistical results of PSNR, SSIM, and SSE with 500 simulations when the noise of SNR is 15 dB are presented in Figure 9 and

Table 3. Mean PSNR (dB), SSIM and SSE results with 500 simulations of Figure 8.

Methods	PSNR (dB)	SSIM	SSE
TV	19.322	0.808	0.422
TVBC	19.633	0.822	0.407
TVS	19.800	0.814	0.399
TVSBC	20.358	0.823	0.374

, and the results of TVBC and TVSBC are finer than those of TV and TVS. Specifically, the mean improvement ratios of TVBC and TVSBC with regards to TV and TVS are 1.61% and 2.82%, 1.72% and 1.08, 3.57%, and 6.23%, respectively. In summary, the application of bias correction can effectively improve the quality of the restored signals when based on TVS or TV.

5.3. Experiment2: 2-D Area Data Processing

In this section, we simulate the 2-D data to discuss the performance of the proposed bias-corrected methods TVBC and TVSBC. The convolution kernel is the same as the 1-D target simulation. The scanning region is $\pm 5°,$ and the range of breadth is 500 m. Nine random-size targets are designed and distributed in the area scene, as illustrated in Figure 10a. Additionally, the noise under high-, middle- and low-SNR conditions are simulated (5 dB, 15 dB, and 25 dB, respectively). Comparing with Figure 10a and Figure 10b–d, we can see that the adjacent targets at the first and second lines are distinguishable, but those at the third line are blurred.

The results generated by these different methods are demonstrated in Figure 11, Figure 12 and Figure 13. As shown in Figure 11 with SNR = 5 dB, for the Blind Deconvolution method, improvement is poor, whose shapes are blurred, and the noise still exists. Regularized Filter, Wiener Filter, TSVD, and SDBSM can remove noise sufficiently, but the obvious targets are fuzzy. In particular, TSVD and SBDSM produce more fake targets, and even partially recognizable targets are replaced with these fake targets. It shows that these methods are vulnerable to noise. Richardson–Lucy methods have a balance between denoising and deconvolution. In other words, it restores fairly clear targets, which can be helpful in positioning the targets. However, the shape of targets is lost by comparison with the original real echo. As shown in Figure 11g–j, the restored results of TV, TVBC, TVS, and TVSBC can not only distinguish all adjacent targets and suppress the noise but also maintain clear shapes and similar sizes. In particular, for the low intensity of the target, TVBC can get a more obvious rectangle shape than that of TV. TVS and TVSBC have similar results without significant differences in visual, but the size of TVSBC is closer than that of TVBC to the simulated echo. As listed in

Table 4. PSNR, SSIM, and SSE evaluation results of Figure 11, Figure 12 and Figure 13.

Methods	SNR = 5 dB			SNR = 15 dB			SNR = 25 dB
	PSNR (dB)	SSIM	SSE	PSNR (dB)	SSIM	SSE	PSNR (dB)	SSIM	SSE
TV	17.462	0.595	29.978	19.999	0.773	22.387	0.773	0.882	20.152
TVBC	18.377	0.634	26.981	20.431	0.798	21.299	0.798	0.888	19.882
TVS	18.922	0.612	25.341	20.771	0.807	20.481	0.807	0.888	19.796
TVSBC	18.986	0.621	25.153	20.992	0.813	19.967	0.813	0.894	18.979

, all evaluation results of TVBC and TVSBC are superior to those of TV and TVS, respectively. And the improvement ratios of PSNR, SSIM, and SSE with TVBC, TVSBC, and TV, TVS are 5.2% and 0.3%, 6.6% and 1.5%, 10.0%, and 0.7%, respectively.

Similarly, the increases in SNRs in Figure 12 and Figure 13 have fewer effects on these restored results, but the results of Blind Deconvolution, Regularized Filter, Wiener Filter, TSVD, and SDBSM are relatively poor by comparison with those of TV, TVBC, TVS, and TVSBC. The results of the Richardson–Lucy method can get better with the decrease of the noise level, but the shape of it cannot be restored as well as TV, TVBC, TVS, and TVSBC. The line chart of assessment for PSNRs, SSIM, and SSE of these methods is demonstrated in Figure 14. It can be seen that the Richardson–Lucy method is ranked only second to TV, TVBC, TVS, and TVSBC, which is consistent with what is displayed in Figure 11, Figure 12 and Figure 13. For middle SNR, the improvement ratios of PSNR, SSIM, and SSE with TVBC, TVSBC, and TV, TVS are 2.2% and 1.1%, 3.2% and 0.8%, 4.9%, and 2.5%, respectively. For high SNR, the improvement ratios of PSNR, SSIM, and SSE with TVBC, TVSBC, and TV, TVS are 0.6% and 1.7%, 0.6% and 0.7%, 1.3%, and 4.1%, respectively. For middle and high SNRs of noise, TV, TVBC, TVS, and TVSBC can get fine restored results superior to the others. In particular, bias-corrected methods can have better results than those without bias correction, which demonstrates that this work is meaningful and feasible.

6. Discussion

An ill-posed inversion problem is usually solved by the regularization methods with different regularization terms to constrain the parameter space for stable solutions. In typical Tikhonov regularization based on the L₂ norm, the biases of solutions are paid more attention to be analyzed and better corrected. For example, the fully-, partially- and adaptively- bias-corrected methods are raised [36,37,42]. In contrast, the sparsity regularization terms based on the L₁ norm are utilized to address the ill-posed sparsity inversion problems, but the biases of its solution are seldom involved and analyzed. To this end, this paper places emphasis on the analysis of the biases of regularization solutions based on the L₁ norm. From the perspective of the formulae comparison and derivation, the differences between regularization solutions based on the L₁ norm and L₂ norm are obtained. The partially-biased property of the regularization solution based on the L₁ norm is revealed herein, the bias of which is less than that of the L₂ norm by analysis in Section 2 when both are without bias correction. It reflects that the results of the L₁ norm will be better than the L₂ norm, which is consistent with the restored results from TV and TVS displayed in Section 4. The improvements of TVBC and TVSBC for PSNR, SSIM, and SSE overall are 2.15%, 1.8%, and 4.14% to those of TV and TVS, respectively.

Additionally, for the partially biased property, we derive the piecewise form of the partially bias-corrected solution in Section 3. In other words, the proposed partial bias correction can retain original information in low frequency and reduce the negative effect in the high frequency of the coefficient matrix. The corresponding flowchart for the proposed TVSBC and TVBC is designed for better implementation. The idea of bias analysis and correction for the L₁ norm solution is rigorously extended in theory. In practice, for different prior regularization terms based on the L₁ norm, the biases of its solutions can be analyzed and corrected, referring to ours.

Inspired by [36,37,42,43], the quality of the regularization solution is relative to the sparsity regularization parameters selection, which is empirically determined, however, which is empirically selected in this paper. Therefore, we will investigate the relation between parameters selection criterion and bias correction in future works.

7. Conclusions and Outlook

In this paper, an analysis of the bias of the ill-posed sparsity inversion problem based on the L₁ norm is implemented. The partially biased property of the L₁ norm is revealed. Without bias correction, the solution of the L₁ norm is superior to that of Tikhonov regularization because of its own partially biased nature. Furthermore, on the basis of it, we give the partially bias-corrected solutions of the TVS and TV methods, which are dependent on the sparsity property L₁ norm. In a word, the partial bias correction method avoids modifying the part of the unbiased solution, which improves the rigor of the theory. In addition, the bias-corrected method of TVBC and TVSBC are helpful in improving the quality of restored signals by retaining more signals in low-frequency and suppressing noise in high-frequency more effectively. The robustness of TV and TVS has also been enhanced.

Two experimental examples have been presented to demonstrate the performance of the proposed methods, including 1-D point target and 2-D area simulation for super azimuth resolution of radar forward-looking imaging. The bias-corrected methods can preserve higher quality than the corresponding methods without bias correction, whose results not only distinguish the adjacent targets, suppress the noise, and preserve the shape and size of targets in visual, but also have excellent statistical results of PSNR, SSIM, and SSE assessment indexes. It infers that the partially bias correction idea is worth being promoted to improve the theoretic rigor of sparse regularization solution and the quality of the restored signals.