Iterative Regularization Using an Acceleration Technique for Synthetic Aperture Interferometric Radiometer

The inverse problem of synthetic aperture interferometric radiometers (SAIRs) aims at retrieving the brightness temperature map from the visibility function samples. To efficiently obtain an accurate solution, an accelerated iterative regularization technique is proposed to reconstruct the SAIR brightness temperature map. The acceleration technique modifies the conventional least square function using a negative penalty term to speed up the initial iterations. A series of decreasing coefficients must be tuned to prevent noise amplification in subsequent iterations. Several numerical experiments were performed on one-dimensional (1-D) and (2-D) SAIR systems. The experiment results indicate that the accelerated Landweber method reduces the number of iteration steps by more than 50% and effectively improves the computational speed without reducing the reconstruction accuracy compared to the conventional Landweber method.

The measurement outputs of real aperture microwave radiometers are directly the brightness temperature of the target radiation. Unlike real aperture radiometers, SAIRs produce samples of the visibility function in the spatial frequency domain. The relationship between the visibility function and the radiometric brightness temperature distribution can be denoted as the Corbella equation [9].
The aim of SAIR reconstruction is to obtain the radiometric brightness temperature map from the measured visibility function. For an ideal SAIR system, the modified brightness temperature map can be retrieved by an inverse Fourier transform [10]. However, there exist antenna pattern differences and spatial decorrelation effects for a real SAIR instrument. Despite these complexities, the SAIR imaging problem can be solved as a linear inversion problem in general. The Corbella equation that relates the brightness temperature map to the visibility function is discretized as where the vector T ∈ R n denotes the unknown brightness temperature, the vector V ∈ R m indicates the measured visibility function, and the matrix G is the rectangular matrix from R n to R m with m < n. Consequently, the discrete problem to be solved is underdetermined. Since the mapping is not injective, there exists no inverse matrix for the underdetermined system. In addition, the inverse problem is ill-posed since the matrix G is ill-conditioned in practice. For SAIR reconstruction, regularization methods are typically used to alleviate the ill-condition of the inverse problem to produce a stable solution. Currently, the commonly used methods are band-limited regularization [11] and conventional direct regularization methods such as minimum-norm regularization and Tikhonov regularization [12], [13]. Moreover, alternative approaches, such as Bayesian inference and reweighted total variation [14], [15], [16], [17], are also presented to retrieve the brightness temperature distribution in SAIRs.
Because iterative regularization methods can take advantage of the sparsity and structure of the matrix better than the direct regularization methods, the former outperforms the latter when the size of the matrix G is large. However, a disadvantage of iterative methods is that they generally have a slow convergence rate [18]. Some sophisticated techniques, such as preconditioning [19], [20], have been presented to speed up the convergence of iterative algorithms. Besides, a simple acceleration technique This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ based on an "irregularization" concept has recently been proposed [21], [22]. The key idea of the aforementioned acceleration technique is to modify the traditional least square function by using a negative penalty term. In this article, the iterative regularization using this acceleration technique is presented to efficiently reconstruct an accurate brightness temperature map in SAIRs. To our knowledge, this is the first time that iterative regularization combined with this acceleration technique is applied to SAIR reconstruction.
The remainder of this article is organized as follows: after a brief primer on the common minimum-norm regularization in Section II, the iterative regularization using an acceleration technique is described. In Section III, several simulation experiments are performed on one-dimensional (1-D) and 2-D SAIR systems to validate the performance of the proposed method. Section IV presents the relevant discussions, and Section V draws the conclusion.

A. Minimum-Norm Regularization
In this subsection, we briefly introduce the conventional minimum-norm regularization method widely used in the data processing of SAIRs.
The minimum-norm regularization refers to finding the solution with the minimum energy and solving the following optimization problem where the symbol · 2 stands for the Euclidean norm. To solve (2) quickly, singular value decomposition (SVD) of the matrix G is introduced as where = diag(σ 1 , σ 2 , · · ·, σ m ) ∈ R m×n is the diagonal matrix, σ i represents the singular value of the matrix G in descending order; A = (α 1 , α 2 , · · ·, α m ) ∈ R m×m , α i denotes the left singular vector; B = (β 1 , β 2 , · · ·, β n ) ∈ R n×n , β i denotes the right singular vector; B * represents the adjoint of B.
Subsequently, the ill-posed property of matrix G is alleviated by discarding those smaller singular values that have a larger impact on the error propagation. In consequence, with the help of the truncated SVD, the solution of (2) is given by where p < m is the number of singular values kept in the inversion process.
It should be noted that the regularization parameter p indeed plays a crucial role in the performance of the minimum-norm regularization. When the noise cannot be determined, an appropriate approach to estimate the regularization parameters is the generalized cross-validation (GCV) method [23], [24]. The rationale of the GCV approach is that the new equation system generated can predict well the removed component when a component in the original data is removed. To estimate a befitting regularization parameter, the following function is minimized: where tr stands for the trace of the square matrix. Therefore, the GCV method in this article is exploited to choose the regularization parameter of the minimum-norm regularization.

B. Iterative Regularization Using an Acceleration Technique
The primary methods refer to the minimization of the quadratic residual functional Simple iterative regularizations, such as the Landweber and conjugate gradient algorithms, converge very slowly, resulting in high stability in the inversion process [18]. However, for some iterative methods such as the preconditioned conjugate gradient algorithm, early iterations may lead to noise amplification and instability, especially for severely ill-conditioned inverse problems [18], [25]. In this case, to overcome the instability problem, further regularization constraints need to be introduced by adding a new penalty term. For instance, the classical Tikhonov regularization is to minimize the following regularized least square function where η > 0 denotes the regularization parameter. Using simple iterative regularizations like the Landweber method, the kth iteration solution T k can be expressed as [26] where ρ k−1 >0 denotes the step size, and Q k−1 represents an ascent direction at point T k−1 for the least square functional Φ.
For the specific Landweber method, the step length ρ k−1 = ρ ∈ (0, 2/ G 2 2 ) stands for a constant and Consequently, the reconstructed result by the Landweber method can be expressed as The main disadvantage of standard iterative methods is usually a slow convergence rate that limits their practical use, especially when the modeling matrix is on a large-scale. In this article, we present a simple acceleration technique to speed up the convergence of iterative regularizations for SAIR reconstruction. The basic idea of the proposed acceleration approach is described below.
First, by introducing the S-weighted metric T 2 S , the following minimization problem is solved [27] min Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
where the S-weighted metric T 2 S is defined as ST, T with the symbol ·, · being the inner product, S denotes a positive semi-definite self-adjoint operator, and ε < 0 indicates a regularization weight. Here, the parameter ε is negative to increase the convergence rate. As mentioned earlier, the slow iterative method depends on the regularization effect of the iterative process itself. When the iterative method converges very slowly, the regularization effect might be too strong. To improve the convergence speed, we reduce the regularization effect in the first iterations by using a negative penalty term called "irregularization" [21], [22].
It is worth noting that the operator S is now not fixed but depends on the modeling operator G. Here, the operator S is given by [27] where I is the square identity matrix from R n to R n . After substituting (3), (11) can be transformed into where From the perspective of filter analysis, the operator S acts as a model-dependent high-pass filter. As is well known, noise and signal subspaces can be distinguished by the SVD in the inverse problem. When T is largely in the noise subspace concerned with the smallest singular value σ m of G, the metric T 2 S becomes large. Whereas the metric T 2 S becomes small when T is heavily in the signal subspace related to the largest singular value σ 1 of G.
Second, using the Landweber iterative scheme, the solution of (10) can be obtained by where ρ ∈ (0, 2/ G 2 2 ) represents the step size, γ k = ρε k . By substituting (11) into (13), we get The algorithm (14) is called the improved Landweber approach for convenience. By comparing (9) and (14), we notice that the two equations are identical in form, except for the coefficients in front of the matrices or vectors. In all the coefficients, γ k and ρ are real numbers, and the norm G 2 2 in (14) can be calculated in advance. Therefore, the computational complexity of each iteration of the improved Landweber approach can be considered to be the same as that of the traditional one.
Furthermore, we analyze the computational complexity of different algorithms. The computational complexity of the two Landweber methods is O(kmn), where k represents the number of iterations required when the iteration stops. In contrast, the computational complexity of the minimum-norm regularization is O(mn 2 ). In consequence, when k is less than n, the computational complexity of the two Landweber methods is lower than that of the minimum-norm method. One remark to make is that this acceleration technique can be applied to accelerate other gradient-like iterative approaches in addition to the Landweber method described above [28], [29].
Selecting the parameter γ has a certain influence on the effectiveness of the proposed acceleration method. To ensure convergence and prevent noise amplification in subsequent iterations, the parameter γ need to be a function of the iteration number and gradually approaches zero [21]. Consequently, the parameter γ is updated by [21] where γ 0 ≥ 0 is the constant parameter. When γ 0 = 0, the proposed method becomes the traditional Landweber method. Finally, to filter out the Gibbs effect on account of bandlimited measurements, an appropriate apodization function is generally applied to obtain the final solution T w = U * WUT in SAIRs, where U and W are the Fourier transform and the window function operators respectively [12], [30]. Note that the reconstructed maps must be compared to T w = U * WUT, which is the ideal map apodized with the same window function. In addition, to quantitatively evaluate the reconstruction methods' performance, we introduce two commonly used norms. The first one is the root mean square error (RMSE), defined as follows: where Q denotes the number of pixels in the alias-free field of view (AFFOV). The second norm is the peak signal-to-noise ratio (PSNR), which is defined by

III. SIMULATION RESULTS
In this section, the performance of the improved Landweber approach is compared with the conventional Landweber approach and minimum-norm regularization through several numerical experiments that have been conducted.

A. Validation for a 1-D SAIR System
First, simulation experiments are performed on the L-band full polarization interferometric radiometer (FPIR) system [31], a 1-D SAIR system. The main parameters of the simulated 1-D SAIR system are given in Table I. The sixteen antennas in the system are distributed in different positions with an integer multiple of the minimum baseline Δd = 0.589λ 0 , forming full baseline coverage from Δd to 90Δd (please refer to Fig. 1 in [13]). It is worth noting that the voltage patterns of the sixteen antennas are anisotropic, meaning there are amplitude and phase differences. Moreover, to mimic realistic measurements, we also consider the spatial decorrelation effect in the simulation experiments.
In the first experiment, the original distribution concerns the simulated reference brightness temperature profile, shown in Fig. 1. The double rectangle profile in Fig. 1 can be seen as mimicking the abrupt discontinuities related to the sea/land edge. For the FPIR system, the direction cosines ξ shown as the horizontal axis in Fig. 1 denotes sinθ, where θ is the incidence angle in the across-track dimension measured from the nadir. After generating the ideal visibility function samples based on (1), we use the zero-mean Gaussian noise to simulate the noise interference that inevitably appears in actual measurement.
When the standard deviation of the additive Gaussian noise is δ = (T A + T R )/ √ Bτ = 0.078 with the antenna temperature T A = 250K and the receivers' noise temperature T R = 100K, the reconstructed brightness temperature profiles using the improved Landweber method are compared with those using the minimum-norm and Landweber methods. The error profiles between the reconstructed profiles and the initial one in the AFFOV are presented in Fig. 2. The Blackman window has been used for all the reconstructions and the initial profile, which makes all the brightness temperature profiles at the same resolution. It is noted that the above three methods produce satisfactory and similar inverse results with oscillation ripples effectively suppressed. Besides, we calculate the RMSE and PSNR values to quantitatively compare the performance of the three regularization methods, as given in Table II. The results display that the performance (RMSE or PSNR) for the Landweber method is the same order of magnitude as the minimum-norm regularization, and the performance for the improved Landweber method is almost identical to that of the Landweber one.  Furthermore, we analyzed the convergence performance of the Landweber and improved Landweber approaches. The convergent curves of the fitness function with the number of iterations are shown in Fig. 3, where the fitness function is V − GT 2 . When the curve converges, the number of iterations is 217 for the Landweber approach, whereas the number of iterations is reduced to 105 for the improved one. The results reveal that the improved Landweber method reaches convergence in less than half of as many steps as the conventional one.
Note that the convergence performance of the improved Landweber method is related to the parameter γ. For the double rectangle reference profile in Fig. 1, the influence of the parameter γ 0 on the convergence performance is shown in Fig. 4. One can notice that as the number of iteration steps increases, the RMSE drops rapidly until approaching a stable value. This implies that although γ has a specific influence on the convergence speed, if γ 0 is within a reasonable range, such as [2 11 , 2 13 ], the iterations of the improved Landweber algorithm at the convergence will be significantly less than that of the Landweber algorithm, i.e., 217. It's worth noting that when γ 0 ≥    14 , the residual norm of the solution is no longer monotonically decreasing because too much regularity is subtracted within the first iterations.
As already discussed in the previous section, the calculation time of the Landweber and improved Landweber methods is concerned with the number of iterations and the size of the modeling matrix. In the first experiment where the size of G is 241 × 800, the running time of the Landweber and improved Landweber algorithms is separately 0.0162 and 0.0132 s on the MATLAB R2020b with 3.2 GHz AMD R7-5800H processor and 16 GB memory. For improved Landweber method, the running time to calculate the norm G 2 2 is 0.0060s so that the running time of the iteration process is 0.0072 s, which is less than half of 0.0162 s. Moreover, the running time of the minimum-norm method is 0.0481s, most of which comes from calculation process of the SVD and GCV. The results show that the calculation time of the Landweber method is shorter with respect to the minimum-norm method, and the improved Landweber has the shortest time. Consequently, the improved Landweber method performs best in terms of computational speed among the three methods.
After the comparison experiment on the simulated brightness temperature profile, we conducted the second experiment using the real brightness temperature data from the Aquarius satellite, which was collected on August 19, 2013. The reference brightness temperature distributions for the ocean and the sea/land edge representing two typical observed scenes are depicted in Fig. 5(a) and (b), respectively.  In the case of adding the Gaussian noise with the standard deviation δ = 0.078, the error between the reconstructed distributions and the initial distribution in the AFFOV are presented in Fig. 6 via the minimum-norm, Landweber and improved Landweber methods. In Fig. 6, the Blackman window is used for all the distributions. As shown in Fig. 6, the reconstruction error distributions of the three approaches are very close for both the ocean and the sea/land edge. Furthermore, Table III lists the RMSE and PSNR values in connection with reconstruction results in Fig. 6. The results indicate that compared to the minimum-norm regularization, the Landweber approach has slightly better performance with the lower RMSE and the higher PSNR. Besides, the RMSE and PSNR values for the two Landweber approaches are almost the same.
In addition, the convergence performance of the Landweber and improved Landweber methods is depicted in Fig. 7(a) and (b), which correspond to the ocean and the sea/land edge, respectively. In Fig. 7(a), when the convergence is achieved, the number of iterations for the Landweber and improved Landweber approaches are 197 and 84, respectively. In Fig. 7(b), when the convergence is reached, the number of iterations for the Landweber and improved Landweber approaches is separately 207 and 93. The results imply that the improved Landweber method converges distinctly faster than the Landweber method.
Furthermore, we compare the computation time of different regularization approaches for the ocean and the sea/land edge, which is listed in Table IV. It can be noted that for two different observed scenes, the computation time of the improved Landweber approach is the shortest among the three regularization approaches.    Finally, the robustness of regularization methods against noise interference is analyzed quantitatively. The measured visibility functions stained by the Gaussian noise, whose standard deviation δ changes from 0 to 0.2, are reconstructed to obtain the brightness temperature distributions. Figs. 8 and 9 show the performance of regularization methods against noise interference for the ocean and the sea/land edge, respectively. The results reveal that though the level of the added noise varies, the performance (RMSE or PSNR) of the Landweber method is a little better than that of the minimum-norm regularization. Additionally, for both the ocean and the sea/land edge, the two Landweber methods have nearly the same RMSE and PSNR values. This confirms that the improved Landweber method shows similar robustness to noise interference as the Landweber one.

B. Validation for a 2-D SAIR System
Numerical experiments are conducted on a 2-D SAIR system to further prove the proposed method's effectiveness. The antenna array layout is T-shaped with 121 antennas, which generates a rectangular visibility function distribution. Like the aforementioned 1-D SAIR system, the 2-D system also considers the anisotropic antenna patterns with amplitude and phase differences [11] and spatial decorrelation effects to simulate realistic measurements. Table V lists the specific parameters of the 2-D SAIR system. Based on the above 2-D SAIR system, we carry out simulation experiments. The original brightness temperature image is presented in Fig. 10(a), which is obtained by measuring the top of a tower and the corresponding optical image is shown in Fig. 10(b). The noise-corrupted visibility functions, formed by adding the Gaussian noise with zero mean and the standard deviation δ = 0.078 with T A = 250K and T R = 100K, are retrieved using the minimum-norm, Landweber and improved Landweber methods. Fig. 11 shows the error maps between the original and retrieved brightness temperature maps from different reconstruction methods using the Blackman window. By comparing the results in Fig. 11, we can see that the Landweber method produces a slightly better result relative to the minimum-norm method [see Fig. 11(a) and (b)]. Additionally, the error maps are  almost identical for the Landweber and improved Landweber methods [see Fig. 11(b) and (c)]. Moreover, to quantitatively appraise the performance of the reconstructed results, the RMSE and PSNR values are calculated and given in Table VI. The results witness that the Landweber and improved Landweber methods produce the reconstruction results with similar RMSE and PSNR performance, which is a little better than that for the minimum-norm regularization.
Furthermore, the convergence performance of the Landweber and improved Landweber methods is presented in Fig. 12, where the fitness function is V − GT 2 . As can be seen from this figure, when the curve converges, the number of iterations for the Landweber and improved Landweber methods is 2331 and 1124, respectively. The results reveal that the improved Landweber method has less than half the number of iterations as the conventional Landweber method.
We compare the computation time of the three approaches mentioned above. The running time (MATLAB R2020b on a PC with 3.6 GHz Intel i9-9900k processor and 64 GB memory) of the minimum-norm regularization is 2989.3 s, where the SVD takes 2889.8 s. Whereas, the running time of the Landweber and improved Landweber approaches is 1803.1 and 980.2 s, respectively. It can be noted that when the modeling matrix is large, the computation time of the minimum-norm approach is very long, where the SVD occupies the majority. Besides, the computation time of the Landweber method is obviously lower than that of the minimum-norm regularization. Furthermore, compared with the Landweber method, the improved Landweber method can effectively reduce the computation time, demonstrating the effectiveness of the acceleration technique.
Finally, we evaluate the robustness of regularization methods against noise interference quantitatively. The brightness temperature maps are obtained by inverting the noise-corrupted visibility functions, stained by the Gaussian noise with different levels. The performance of regularization methods against noise interference is depicted in Fig. 13. From Fig. 13, one can notice that even though the noise level changes, the performance for the Landweber method slightly outperforms the minimumnorm regularization. Moreover, the Landweber and improved Landweber methods have very close performance. As a result, the improved Landweber method does not reduce the accuracy of the inverse results and the robustness against noise interference with respect to the traditional Landweber one.

IV. DISCUSSION
For SAIRs, the reconstruction from the measured visibility function to the brightness temperature map of the observed scene is an ill-posed inverse problem. When the modeling matrix is large-scale, iterative regularizations are superior to traditional direct regularization approaches such as the minimum-norm and Tikhonov regularizations. An essential issue for iterative methods is the tradeoff between computation speed and smoothing effects. Although conventional iterative methods can achieve excellent smoothing results, they can also lead to slow computation speed due to a slow convergence rate.
To speed up the convergence of iterative regularizations in SAIRs, an acceleration technique is presented in this article. It speeds up the first iterations by modifying the conventional least square function using a negative penalty term. Moreover, the negative penalty term must be tuned by decreasing coefficients to prevent noise amplification in subsequent iterations. Furthermore, we analyze the computational complexity of the accelerated Landweber method relative to the traditional Landweber iterative method and the minimum-norm regularization.
Through several simulation experiments, we compare the performance of the accelerated Landweber method with the conventional minimum-norm and Landweber methods in terms of inverse accuracy and calculation time. The experimental results show that compared with the minimum-norm regularization, the inversion accuracy of the Landweber method is close, sometimes slightly better. The inversion accuracy of the accelerated Landweber and Landweber methods is almost the same. In addition, the calculation time of the Landweber method is lower than that of the minimum-norm regularization. Compared with the Landweber method, the accelerated Landweber method achieves a reduction of over 50% in the number of iterations required. Although the accelerated Landweber method needs to take some extra time to calculate the norm G 2 2 , it takes markedly less time in total than the Landweber method, especially for large-scale SAIR systems. Consequently, the proposed acceleration technique offers significant improvements in computational speed compared to state-of-the-art iterative regularizations.
However, there are still some shortcomings that can be further explored, summarized as follows.
1) An important disadvantage of iterative methods is coming from the fact that the computational time may vary from one set of visibilities to another. This is a serious issue for an operational implementation and will be further studied in future. 2) The Landweber approach is a simple iterative regularization. In this article, we elaborate on how to use the acceleration technique to accelerate the Landweber approach. On this basis, we will further study how to combine the acceleration technique with some complex iterative regularization methods such as conjugate gradient least squares.

V. CONCLUSION
In this article, we present an accelerated iterative regularization technique to reconstruct the SAIR brightness temperature map. The acceleration technique modifies the conventional least square function using a negative penalty term to speed up the initial iterations. A series of decreasing coefficients must be tuned to prevent noise amplification in subsequent iterations. In addition, the computational complexity of regularization methods is analyzed. Several simulation experiments were performed on 1-D and 2-D SAIR systems. The experiment results indicate that compared to the conventional Landweber method, the accelerated Landweber method reduces the number of iterations by more than 50% and effectively improves the computational speed without reducing the reconstruction accuracy.