Reduction of the Reconstruction Error With Lower and Upper Bounds in Synthetic Aperture Imaging Radiometers

Synthetic aperture imaging radiometers (SAIRs) are powerful instruments for high-resolution Earth observation by use of small-aperture antennas sparsely arranged to achieve a large-aperture antenna. High-precision reconstruction algorithm is one of the key contents of SAIRs. Owing to the ill-posed problem and band-limited physical characteristic, there is a still large residual error for traditional regularization methods. It should be noted that the prior information like the lower and upper bounds of the brightness temperature distributions has not been utilized in the reconstruction procedure, especially for the open ocean with relatively small brightness temperature difference. In order to reduce the reconstruction error in SAIRs, a reconstruction method based on active set algorithm is presented by solving the least squares problems with lower and upper bounds. The simulation experiment results show that the proposed method can more effectively reduce the reconstruction error and better improve the accuracy of retrieved brightness temperature distributions in SAIRs than the band-limited regularization method, demonstrating the effectiveness of the proposed method.


I. INTRODUCTION
Synthetic aperture imaging radiometers (SAIRs) are passive microwave sensors by use of synthetic aperture technique. For real aperture radiometers, increasing the antenna aperture is the only way to improve the spatial resolution. However, large aperture antenna brings the difficulties of mechanical scanning, the weight as well as the bulky volume. In order to overcome the above shortcomings of real aperture radiometers, SAIRs improve the spatial resolution by sparsely arranging small antennas. With the deepening of theoretical and systematic research, SAIRs have entered practical applications [1] and researchers have developed some instruments such as Soil Moisture and Ocean Salinity (SMOS) [2], the Geostationary Synthetic Thinned Array The associate editor coordinating the review of this manuscript and approving it for publication was Stefania Bonafoni .
It has been demonstrated that the inverse problem of SAIRs is mathematically ill-posed so that the solution is neither unique nor stable [6]. Therefore, the inverse problem needs to be regularized in order to provide a stable and unique solution. Currently, regularization methods have become the dominant methods to solve the inverse problem of SAIRs and get excellent results [7]- [9]. The direct regularizations, such as Tikhonov regularization and truncated singular value decomposition (TSVD), overcome the ill-conditioned property of the inverse problem by use of numerical methods, but they need to choose the suitable regularization parameter. Band-limited regularization proposed by [6] cures the ill-posed property of the inverse problem by considering the physical characteristic of limited bandwidth for SAIRs. However, owing to the ill-posed property and band-limited physical characteristic, there is still a large residual error after reconstructing [10], [11].
An important application of SAIRs is to detect Sea Surface Salinity (SSS), which requires that the measurement accuracy of the ocean brightness temperature is within 0.1K. Residual error for traditional reconstruction methods needs to be effectively reduced to meet the detection needs of geophysical parameters such as SSS. It should be noted that the three regularization methods mentioned above do not make use of any prior information about brightness temperature images of the observed scene in the spatial domain. Actually, it is easy to obtain some prior information such as the upper and lower bounds, especially for the observed scene like the open ocean with relatively small dynamic range of the brightness temperatures.
The principle of SSS remote sensing is based on the salinity sensitivity of sea surface brightness temperatures at microwave frequencies [12]. The dielectric constant models have been widely used to study the sensitivity of microwave radiation to the water salinity. It has been shown that the Ellison model proposed by [13] can be well consistent with satellite microwave observations of sea surfaces [14]. The salinity for open oceans is in the range of 32-37 practical salinity unit (psu). According to the Ellison model, for a smooth water surface of the open oceans at 40 • incidence angle and 5 degrees Celsius sea surface temperature, the average of the vertically and horizontally polarized brightness temperatures (T h + T v )/2 is within the scope of 92-96 K ( Figure 2 in [15]).
After considering the prior information, the inverse problem in SAIRs can be transformed into the least squares problem with lower and upper bounds. A direct method called active set algorithm has been suggested for solving linear least squares problem with lower and upper bounds [16], [17]. The prominent characteristic of this method is that the iteration point can follow the constraint of the upper and lower bounds until it reaches the optimal solution of the inverse problem. In this paper, active set algorithm is applied to reduce the reconstruction error in SAIRs.

II. THE IMAGING PRINCIPLE
By measuring the complex correlation between the signals collected by two spatially separated antennas, which have an overlapping field of view, SAIRs yield samples of the complex visibilities of the brightness temperature of the observed scene. The relationship between the measured visibility function V (u) and the brightness temperature images T B ( ξ ) is given by [18] V (u kl ) where u kl stands for the spatial frequency associated with the two antennas, and the Cartesian coordinates ξ = (sinθ cosφ, sinθ sinφ) are direction cosines (θ and φ are the traditional spherical coordinates). F k (ξ ) and F l ( ξ ) stand for the normalized voltage patterns of the two antennas with equivalent solid angles k and l , and the overbar is the operator of the complex conjugate.r kl (t) is called as the fringe-washing function, t = -u kl ξ/f 0 stands for the spatial delay, and f 0 is the observed central frequency. T r is the physical temperature of the receivers and assumed to be the same for all the receivers. After discretization, (1) can be written in the matrix equation where G stands for the modeling operator, T and V represent the discrete brightness temperature distribution and the visibility function samples, respectively. Since the number of pixels in T is larger than the number of samples of V, the linear system is underconstrained and has multiple solutions. Hence, the minimum of the least squares criterion is not unique owing to the singularity. Moreover, small disturbance in the observed visibility samples would lead to a very large perturbation in reconstructed brightness temperature distribution. Therefore, the inverse problem is ill-conditioned and has to be regularized so as to provide a unique and stable solution.
By taking into account the physical constraint of the limited resolution for SAIRs, the band-limited regularization method is to solve the constrained optimization problem [6] where I is the identity matrix, U H = F * ZZ * F is the projection operator (F * is the adjoint of F), F is the Fourier transform operator, and Z is the zero-padding operator beyond the experimental frequency coverage H . The unique solution of (4) is given by where J + = (J * J) −1 J * is the Moore-Penrose pseudoinverse of the matrix J = GF * Z.

III. A RECONSTRUCTION METHOD WITH LOWER AND UPPER BOUNDS
Error sources for microwave remote sensing of SSS include the sea surface temperature, sea surface roughness, ionospheric Faraday rotation, atmospheric gases and solar system and galactic radiation. It has been pointed out by [15] that most of the error sources are expected to produce the brightness temperature effects of a few kelvin. In this paper, the Ellison model, which is combined with the effects of the error sources such as sea surface roughness and Faraday rotation described in [15], is used to predict the range of the sea surface brightness temperatures at 1.4 GHz, namely the upper and lower bounds.
The cumulated brightness temperature in ocean from May 6-9, 2010 for SMOS satellite [19] is shown in Figure 1, which shows stable range of the brightness temperatures compatible with global salinity maps. From Figure 1, we can find that the upper and lower bounds of the brightness temperature distributions in the open oceans are 85K and 95K, respectively. Therefore, the brightness temperature distribution in the open oceans is relatively stable and has upper and lower limits.
By considering the bounded property of the brightness temperature distribution, a reconstruction method based on active set algorithm is proposed. The optimization problem can be expressed as where β 1 is the constant vector with the lower bound β 1 and β 2 is the constant vector with the upper bound β 2 . The least-squares problem can be transformed into a quadratic programming problem The bounded constraints can be converted into where I is the identity matrix. If we define the constraints are written as follows: Thus, (6) can be given by standard form of the quadratic programming problem [16] min where D = 2G T G is the symmetric matrix, C = -2G T V, A= (a 1 , a 2 ,. . . , a m ), When D is a positive semidefinite matrix, the objective function in (11) is convex. If the feasible region of the quadratic programming problem is not empty, it must be a convex set. In consequence, (11) is the convex quadratic programming problem. Under the circumstance, the local optimal solutions are the global optimal solutions [17], [20]. In addition, when D is a positive definite matrix, the global optimal solution is unique.
Under general circumstances, the quadratic programming problem with inequality constraints needs to be transformed into the equality constraint problem. In order to solve the quadratic programming problem with the equality constraints, Lagrange multipliers γ = (γ 1 , γ 2 , . . ., γ m ) T are introduced w represents the set including all the effective constraints.
The equality constraint problem can be solved by finding the stationary point of the Lagrange function. Hence, the T and γ satisfy: We can get T and γ by solving (13). If T satisfies the Kuhn-Tucker rules, T is also the optimal point. In this situation, γ should satisfy For the active set algorithm, it is assumed that after k iterations, we can get the feasible point T k and the working set w k . And then check whether the current iteration point is optimal under the current working set. If not, we define a step q [17] By substituting (15) into (11), we get The solution of (16) is given by q k . Consequently, the iteration point is updated by where s k ∈ [0,1] is the step-length parameter. In order to ensure that the new iteration point meets the original constraints, s k needs to satisfy [17] If s k < 1, it means the step q k is blocked by the constraint that is not in the working set. So we construct a new working set w k+1 by adding this constraint that satisfies When q k = 0 and γ i ≥ 0, the current iteration point T k is the optimal solution of (11). Moreover, it has been pointed out by [17] that if the direction q k is non-zero, the objective function will decrease along this direction, which indicates the convergence of the active set method.
The algorithm procedure is summarized as follows in Table 1. In this paper, for purpose of speeding up the convergence, the constant vector with the mean value of the upper and lower limits is selected as the initial solution T 0 .
In order to quantitatively analyze the performance of the retrieved brightness temperature distribution, the root mean square error (RMSE) and peak signal to noise ratio (PSNR) are calculated as follows: (19) where N is the number of the brightness temperature pixels in the alias-free field of view (AFFOV), m T is the maximum of the brightness temperature distribution.

IV. RESULTS AND DISCUSSION
Numerical simulation experiments have been performed in order to verify the effectiveness and performance of the above inverse method. The experiments are based on the L-band FPIR system [21] and the specific parameters of the system are listed in Table 2. Antenna array configuration of the system is shown in Figure 2, where the antennas elements n1, n 2, . . . , n16 are arranged at different positions with the shortest baseline set to d = 0.589λ 0 . After applying the hermiticity property, there are 241 visibility function samples in total. The antenna has been simulated with nonisotropic antenna voltage patterns, shown in Figure 3. In Figure 3, the horizontal axes ξ is the direction cosines in the Cartesian coordinates, and n1, n2, . . . , n16 denote the element antennas in sequence. In addition, the fringe-washing function is set to sinc(Bt) for modeling the operator G.    As shown in Figure 4, the initial distribution used in the simulations derives from the observed ocean brightness temperature of H polarization in the 1.4 GHz channel for the SMAP satellite. In Figure 4, the horizontal axes ξ is the direction cosines in the Cartesian coordinates. When the matrix size of the brightness temperature distribution is 500 × 1, the size of the matrices G is 241×500. The complex visibilities are obtained by simulating the test distribution based on the above FPIR system. In actual measurement, the measurement error or noise interference is generally unavoidable. Gaussian white noises with zero mean and the variance σ 2 are added to the complex visibilities to mimic measurement error and noise interference.
When the variance of the added noise is σ 2 = 0.1 max(V i ), the retrieved distributions in the alias-free field of view by active set method and band-limited regularization method are shown in Figure 5. For active set method, the upper and lower bounds of the brightness temperature distribution are separately set to 105K and 85K according to the Ellison model, and the initial solution T 0 is the constant vector with the value 95. As can be seen from Figure 5, reconstruction result of the band-limited regularization has obvious oscillation ripples especially at the edge of the brightness temperature distribution, which is the Gibbs effect owing to limited bandwidth coverage in the frequency domain for SAIR system. In addition, compared with the band-limited regularization, the active set method produces better result, which has weaker ripples and exhibits better suppression of oscillation ripples, demonstrating the effectiveness of the active set method. The RMSE and PSNR are calculated for quantitative analysis on the reconstruction error. In Figure 5, the RMSE and PSNR for band-limited method are respectively 5.51K and 33.30dB, and the RMSE and PSNR for active set method are respectively 4.19K and 35.68dB.
In order to analyze quantitatively the impact of noises on the retrieved distributions, the Gaussian white noises of different levels are added to the complex visibilities and then reconstructed to obtain the brightness temperature distribution. The performance comparison of the two inverse methods in different noise level is shown in Figure 6, where the upper and lower bounds and the initial distribution for active set method are the same as those in Figure 5. Figure 6(a) indicates that the RMSE for the active set method is lower obviously than the RMSE for the band-limited regularization in the case of high-intensity noises such as σ 2 = 0.1 max(V i ) and the two RMSE values gradually approach when the variance σ 2 decreases from 0.1 max(V i ) to 0.01 max(V i ). It should be noted that in the case of low variance, the RMSE absolute values for the two methods are small so that the relative difference between the two RMSE becomes very small. For example, when the variance σ 2 = 0.01 max(V i ), the RMSE for the band-limit regularization and active set method are 1.86K and 1.63K, respectively. From Figure 6(b), we can find that whether the noise level is high or low, the PSNR for the active set method is 1.8dB higher than that for the bandlimited regularization. Therefore, the results show that the active set method is more robust to the noise interference than the band-limit regularization, proving that it can better improve the accuracy of the retrieved brightness temperature distributions.
In addition, the convergence performance of active set method in different level noises is presented in Figure 7. As shown in Figure 7, when the iteration steps increase, the fitness function decreases rapidly until a stable value, demonstrating the convergence of the active set method. In Figure 7, when the number of iterations is about 10, the fitness function converges to a stable value.  For the active set method, the upper and lower bounds of the brightness temperature distribution in the open oceans are predicted according to the Ellison model combined with the influence of other error sources. However, there may be the errors of a few kelvin or even more than a dozen kelvin between the estimated results and the actual results. Simulations have been performed to quantitatively analyze the impact of different upper and lower bounds on the reconstruction results. The original distribution, where the minimum and maximum brightness temperatures are respectively T min = 89.05K and T max = 101.73K, is shown in Figure 4. In an ideal situation, the lower and upper limits for the active set method are T min and T max . The boundary error of the brightness temperature distribution e is defined as the difference of the minimum T min and the lower bound β 1 . For convenience, it is assumed to e = T min -β 1 = β 2 -T max . When the measured visibility function samples are corrupted by Gaussian white noise with the variance σ 2 = 0.1 max(V i ), the retrieved distributions in the alias-free field of view for different boundary errors are shown in Figure 8. The results show that the boundary error has an impact on the retrieved distribution. The greater the boundary error is, the poorer quality of the retrieved distribution is, particularly on the edge of the distribution.
In order to determine the impact of the boundary error on the retrieved distributions for the active set method, the RMSE and PSNR are calculated. The relation between the boundary error and the performance (RMSE or PSNR) of the retrieved distribution is presented in Figure 9. From Figure 9, we can find that when the boundary error increases, the RMSE gradually increases and the PSNR gradually decreases. Moreover, compared with the band-limit regularization method, the active set method has a lower RMSE and a higher PSNR even when the boundary error is 20K, that is, the upper and lower limits of the brightness temperatures are set to 121.73K and 79.05K, respectively. Therefore, although the boundary error has an impact on the performance of the reconstruction results, the active set method can in general produce better reconstruction results compared to the bandlimit regularization.
Compared with the band-limited regularization method, the computation time of the active set method as an iterative method is mainly related to the iteration count. When the iteration count is 10, the MATLAB runtime (MATLAB2019b on a PC equipped with 2.4GHz Intel i5-4210U processors, 12 GB of memory) of the active set method is 3.65 s. In addition, the MATLAB runtime of the band-limited regularization is 0.78 s. The results show that the computational cost of the active set method is higher compared to the band-limited regularization. In order to speed up the performance of the active set method, the reconstruction code can be run on the GPU hardware platform in the future.

V. CONCLUSION
High-precision reconstruction algorithm is one of the key contents of SAIRs. Different from real aperture radiometers, the output of SAIRs is the sampled data of the visibility function. The purpose of the reconstruction algorithm is to transform the visibility function samples in the frequency domain into the brightness temperature distribution in the spatial domain. Due to the ill-posed property and band-limited physical characteristic, there is still a large residual error for the reconstruction result of the traditional reconstruction algorithm. However, the prior information has not been utilized in the restructuring procedure. A novel reconstruction method, which makes use of the lower and upper bounds of the brightness temperature distributions, is proposed to reduce the reconstruction error. The proposed method is based on active set algorithm, which solves sparse least squares problems with lower and upper bounds. The simulation experiment results show that the proposed method can more effectively reduce the reconstruction error, especially in the case of high-intensity noise, compared to the band-limited regularization. In addition, the upper and lower bounds should be set as close as possible to the minimum and maximum values of the original brightness temperature distribution. He is currently an Associate Professor with the School of Information Science and Technology, Zhejiang Sci-Tech University. His research interests include remote sensing image processing, sparse representation, and machine learning. VOLUME 8, 2020