A Convex Variational Model for Restoring SAR Images Corrupted by Multiplicative Noise

School of Software Engineering, Tongji University, Shanghai 201804, China Shenzhen Key Laboratory of Advanced Machine Learning and Applications, College of Mathematics and Statistics, Shenzhen University, Shenzhen 518060, China Department of Mathematics, Syracuse University, Syracuse, NY 13244, USA Guangdong Key Laboratory of Intelligent Information Processing, Shenzhen University, Shenzhen 518060, China


Introduction
Image restoration problems [1][2][3][4][5][6][7][8][9] have attracted a great amount of attention in real applications. Noise in images is often roughly grouped into additive noise and multiplicative noise. In this paper, we mainly focus on the multiplicative noise removal problem [10], which can be formulated as where f: Ω ⟶ R + is the observed image (Ω ⊂ R 2 is a connected bounded open subset with compact Lipschitz boundary), v ∈ L 2 (Ω) denotes multiplicative noise with mean 1, and K ∈ L(L 2 (Ω)) is a known linear and continuous blurring operator. e goal is to find the unknown true image u from the degraded image f as well as possible.
Multiplicative noise, also known as speckle [11][12][13][14][15], often occurs in active imaging systems, such as synthetic aperture radar (SAR), laser images, and ultrasound imaging. In this paper, the characteristics of speckle which follow a Nakagami distribution in SAR images are considered. In SAR imaging system which are widely used in military reconnaissance, marine monitoring, and other fields, a radar sends a coherent wave. en, it is reflected on the ground and "recorded" by the radar sensor to capture SAR images. e fading of radar echo signal often causes that the SAR images are degraded by speckle, which limits their interpretation and processing.
Over the past decades, there is a great deal of research that has been done on solving the speckled image problem [16][17][18][19]. Early filtering techniques are exploited under the minimum mean square error (MMSE) criterion in the spatial domain [20]. And the bilateral filtering [21], the meshless method [22], the low-rank method [23], and the wavelet-based [24,25] and nonlocal-based [26][27][28] approaches have also been explored for despeckling. Particularly, in some statistics-based methods, the observed image intensity is considered as a random variable following a negative exponential law, and then the intensity of the L-look image which is obtained as the average of L different independent images follows a gamma law with mean equal to one. And thus, the image to be recovered is assumed to be degraded by multiplicative gamma noise. Based on this assumption, a series of variational models, including the ZWM model [29], the AA model, and its variants and convexity improvements (e.g., the exponent-based models [30][31][32][33], the mV model [34], the TwL-mV model [35], and DZ model [36]), are proposed, where the regularization terms generally play a key role to preserve the image structure such as edge detail. On the contrary, in some other statistics-based methods, the image amplitude which is the square root of the reflectivity (intensity) image is taken into consideration. And the pixel amplitudes are modeled as independent and identically distributed according to the Nakagami distribution. Related work includes the probabilistic patch-based (PPB) method [27], the SAR blockmatching 3D (SAR-BM3D) method [28], DD-SRAD method [37], a variation-based model by using fields of experts (FoEs) [38] proposed by Chen et al. [11], and the references cited therein.
Based on the works [11,36,39], we investigate a new convex variational model in this paper, where the image amplitude is considered, and the fidelity term is deduced from the probability density function of the Nakagami distribution and a quadratic penalty term is derived from the statistical property of Nakagami distribution, which guarantees the proposed model to be strictly convex under a mild condition. e convex model guarantees the uniqueness of the solution. Under the framework of continuous functional, we study the existence and uniqueness of the proposed model to the image denoising and deblurring problem. Furthermore, we employ the primal-dual method [36,40] to solve the proposed model. Finally, we arrange some numerical comparisons on simulated images and real SAR images to illustrate the efficiency of the model.

Methods
To recover u for the degraded image f in (1), based on the statistical properties of Nakagami noise, in Section 2.1, we first propose a new variational model (6) for image denoising (i.e., K � I). We analyze the existence and uniqueness of the solution of the new model and discuss the properties of the solution of the model. en, we extend the denoising model (6) to simultaneous deblurring and denoising case and also show the properties of the solution of the model in Section 2.2. And on these bases, we introduce a primal-dual algorithm to solve the proposed model in Section 2.3. By adopting the similar ideas in [36,39], we can get Propositions 1, 2-5, and eorem 1. And the proof of Propositions 1,2, and 4 is given in the Appendices A, C, and D, respectively; and the proof of eorem 1 can be found in Appendix B.

A Convex Variational Denoising Model.
We now introduce a new variational model for image denoising, based on the statistical properties of Nakagami noise. For the multiplicative model f � uv, we assume that v follows a Nakagami distribution, the square root of the reflectivity [11,13]: where L is the number of looks of the image (i.e., number of independent values averaged) and Γ(·) is the classical Gamma function. Based on (2), by using the MAP estimator and combining the resulting data term with the total variation (TV) prior model, we obtain the following AA-like variational model: where λ > 0 and S(Ω): Let z � (1/v) and z denote the instance of random variable Z. It is not hard to get that its PDF is A direct calculation shows that Based on this observation, we therefore introduce a quadratic penalty term Ω ((u/f) − β1) 2 dx into the AA-like model (3), yielding the following optimization problem: where α > 0, λ > 0, and β ≥ 1 vary with the level of the noise. We now show that E(u) in (6) is convex if the parameter α satisfies certain condition.

Proposition 1. If α is chosen satisfying
then the model (6) is strictly convex.
Next, we study the existence and uniqueness of the solution to (6) and the minimum-maximum principle. To this end, we make the following notations. Given a function f ∈ L ∞ (Ω), we denote by sup Ω f (resp., inf Ω f) the supremum (resp., infimum) of f, i.e., sup Ω f � inf C ∈ R; f(x) ≤ C a.e. (resp., inf Ω f � sup C ∈ R; f(x) ≥ C a.e.}).

2
Mathematical Problems in Engineering and inf Ω f > 0, then the problem (6) has a solution u * ∈ BV(Ω) such that Moreover, the solution of (6) is unique when α is chosen to satisfy (7).
We now introduce a comparison principle that characterizes the relationship of solutions of problem (6) with different f values.
in Ω if α and β satisfy

A Simultaneous Deblurring and Denoising Model.
Given a nonnegative, linear, and continuous blurring operator K ∈ L(L 2 (Ω)), the degraded image can be expressed as f � (Ku)v. Here, f is the observed image, u is the image to be recovered, and v is the noise. Similar to (6), the deblurring and denoising model can be formulated as Note that (6) is a special case of (9) with K being the identity operator.
We can show that the above model (9) is convex if α satisfies some certain conditions.

Proposition 3. If α is chosen satisfying
then the model (9) is convex.
In view of Proposition convexdesp and the fact that K is linear, it is not hard to show that the conclusion of Proposition 3holds.
We next show that Ku * is positive; i.e., x ∈ Ω: (Ku * )(x) � 0 { } has a zero Lebesgue measure when ε approaches to zero. And furthermore, we can obtain that Ku * is strictly positive in the discrete situation.
Proposition 4. Suppose f ∈ L ∞ (Ω) and inf Ω f > 0. Let K ∈ L(L 2 (Ω)) be a linear and continuous blurring operator. Assume that u * is the solution of (9). en, there exists a constant C such that, for any 0 < ε < 1, e following Proposition 5 implies that, in general, model (9) cannot automatically preserve the mean of the original image.

Proposition 5.
Let u * be a solution of (7). Suppose that K1 � 1 and inf Ω Ku * > 0. e following statements hold: We can enhance the model (9) to reduce the influence from the bias by keeping the same mean of the recovered image as f by the following model: where E Ω (f) is the mean value of f. Due to the fact that u ∈ S(Ω): E Ω (u) � E Ω (f) is closed and convex, it is easy to see that (13) has exactly one solution.

An Iterative Algorithm.
For the sake of developing numerical algorithms in real-world applications, in this section, we discuss a discrete version of model (6). For simplicity, we use the same notation as the abovementioned continuous version.
and the images u, f ∈ R d , and ∇∈ R 2 d×d be the discrete gradient operator with ∇ x and ∇ y being, respectively, the discrete derivatives in the x-direction and y-direction. With the above notations, the discrete version of model (6) is given as follows: where Mathematical Problems in Engineering e minimization problem (14) can be solved by many numerical algorithms, such as the split-Bregman algorithm [41], the ADMM [42], and proximity algorithm [43]. In this section, we adopt the primal-dual algorithm to solve the problem (14) by formulating it as where div is the divergence operator and P � p ∈ R 2 d : Now, similar to [36], we present the primal-dual algorithm in Algorithm 1. In order to ensure the iterates (u n , p n ) of Algorithm 1 converge to a saddle point of (14), the parameters μ, ρ, and λ should satisfy μρλ 2 < 1/8 [36,44]. And for simplification, we set μ � ρ � 3 in our experiments.

Results and Discussion
In this section, we conduct numerical experiments on simulated images and real SAR images to evaluate the approximation accuracy and computational efficiency of our proposed algorithm. We compare our method in the denoising case with the Lee filters [20], the speckle reducing anisotropic diffusion (SRAD) method [45], the probabilistic patch-based (PPB) method [27], and the model (13) with α � 0 using the algorithm adopted in [39] (here, we call it AA-like method), while in the deblurring and denoising case, our method is compared with the AA-like method. In order to preserve the mean of the observed image, the bias correction technique [36] is utilized for the AA-like method and our method. All the experiments are performed in MATLAB R2017b on a laptop with Intel(R) Core(TM) i7-5500U central processing unit (CPU, 2.40 GHz), 8.0 G memory, and Windows 8.1 operating system.
e quality of recovered images obtained from various despeckling algorithms is evaluated by the structural similarity index (SSIM) [46] and the peak signal-to-noise ratio (PSNR): where u and u denote the original image and the recovered image, respectively; and d is the total number of pixels of u.
All the methods except PPB and Lee filter are terminated by the following stopping criterion: where ε is a fixed threshold. We set ε � 1.5 × 10 − 4 in the denoising experiment and ε � 2 × 10 − 5 in the deblurring experiment.

Denoising Experiment.
e first experiment focuses on the denoising case; i.e., K in model (14) is the identity operator. We first choose the test image of "Boat" and degrade it by multiplicative Nakagami noise with L � 10. With a large value of L, the value of β is 1. We denote Algorithm 2, Algorithm 3, and Algorithm 4 by the same algorithm framework of Algorithm 1 but with the initialization replaced by u 0 � u 0 � u PPB , u Lee , u SRAD , and p 0 � (0, . . . , 0) T ∈ R 2 d . Here, u PPB , u Lee , and u SRAD are the solutions of the PPB method [27], Lee filter [20], and SRAD method [45], respectively. e parameter settings of four algorithms are λ � 0.01, α � 2.7, and β � 1. It can be seen from Figure 2 that although Algorithms 1, 2, 3, and 4 using the same parameters converge to the same stable PSNR value in accordance with the strict convexity and the uniqueness of the solution of the model, the curve of PSNR values associated with Algorithm 2 has a higher peak than that of Algorithms 1, 3, and 4. AA-like method with the same initial guess (named the PPB-AA method) is also added for making the experiments more persuading. Based on this observation, Algorithm 2, PPB-AA, together with Algorithm 1, Lee filter, SRAD method, PPB method, and AA-like method are listed in Table 1 to make a comparison. In Table 1, the images are corrupted by multiplicative noise (Nakagami distribution with L � 10, L � 7, and L � 5). e smaller the L is, the noisier the images are. We tune the parameters to insure a better performance for every image. e parameter settings of proposed algorithms are shown in Table 2. For the Lee filter, the default window size of 5 is used. For the SRAD method, the parameters are set as suggested in [45]. For the PPB method, we use the codes provided by the work in [27]. And for PPB-AA and Algorithm 2, the computing time is expressed as the sum of the CPU cost of getting the initial solution u PPB and true CPU time of Algorithm 2 and AA-like method, respectively. We observe that the PSNR and SSIM values of the restored images by Algorithm 1 are always higher than those of AA-like method, which indicates that the choices of α and β have a significant impact on the quality of denoised images. Moreover, the PSNR and SSIM values of the restored images by Algorithm 2 are generally the highest among all those seven methods tested.
To visually compare the image restoration performance of above seven methods, Figures 3 and 4 report the denoising results of "Lena" and "Remote1" with L � 5, and the certain corresponding partial enlarged denoising images of "Lena" and "Remote1" with L � 10, L � 7, and L � 5 by these methods are shown in Figures 5-10, respectively. Firstly, we find that much noise remains in the results of Lee filter, SRAD method, and AA-like method (e.g., the hat of "Lena" shown in Figures 5-7 and roofs in "Remote1" shown in Figures 8-10). Lee filter leads to a blurred effect on restored images. Since the AA-like method is oversmoothing in the process of despeckling, details of images tend to be reduced. Conversely, more details are preserved by Algorithm 1 and SRAD method. Secondly, the PPB method in [27] can remove the multiplicative noise effectively because it is a patch-based method, which is a different framework from TV-based methods (e.g., AA-like method and Algorithm 1). e images obtained from the PPB method provide smoother regions and better shape preservation (e.g., "Re-mote1" shown in Figure 4 and the backgrounds in "Lena" shown in Figure 3). However, we still notice some artifacts (e.g., the edge of the homogeneous regions in "Lena" and "Remote1"), and the PPB method seems to attenuate the image sharpness (e.g., the hat of "Lena" and trees in "Re-mote1") shown in Figures 5-10. Indeed, those spurious artifacts can be removed by the TV regularization in our method, but Algorithm 1 may suffer from staircasing effects (e.g., the skin of "Lena"). Finally, we observe that Algorithm 2, which combines Algorithm 1 and PPB method, can resolve the problems of PPB method and Algorithm 1 and performs better than the PPB-AA method in attenuating the unpleasant artifacts. Given: a noisy image f ∈ R d + . e positive parameters λ, μ, ρ satisfying ρμλ 2 < 1/8. Initialization: u 0 � u 0 � f and p 0 � (0, . . . , 0) T ∈ R 2 d . repeat Step 1: p n+1 � argmax p∈P λ〈u n , di vp〉 − (1/2μ)‖p − p n ‖ 2 2 .
Step 2: u n+1 � argmin u∈Q F K (u) − λ〈u, di vp n+1 〉 + (1/2ρ)‖u − u n ‖ 2 2 Step 3: u n+1 � 2u n+1 − u n . until u n+1 converges or satisfies a stopping criterion. Output: u * � u n+1 . ALGORITHM 1: Primal-dual iterative algorithm for solving model (14).   blurred, they are corrupted by multiplicative Nakagami noise with L � 10, L � 7, and L � 5, respectively. We list the parameter values in Table 3, and the empirical choice of β is the same as denoising case. e numerical results are presented in Table 4. We can see that the PSNR and SSIM values of the restored images by Algorithm 1 are higher than that of  the AA-like method. Moreover, we show the degraded images and the restored results of both methods for the images "Cameraman" in Figures 11-13. One can see Algorithm 1 preserves more details (e.g., the tripod in "Cameraman"). Due to oversmoothing, the AA-like method tends to lose details and attenuate the image contrast. Note that both methods cost more time than that in the denoising case because of the existence of the blurring operator.

Real SAR Images Restoration Experiment.
Furthermore, we evaluate the performance of a newly proposed model for despeckling based on three real SAR images. e stopping criterion is the same as that of the denoising experiment. Since there are no original noise-free images to be compared, PSNR and SSIM cannot be used to evaluate the quality of recovered images. Assuming L � 5, we tune the parameters until the recovered images show the best results. e parameters settings are reported in Table 5. Figures 14-16   performs well in restoring images, such as smoother regions and better shape preservation. But there are some spurious artifacts in the edge of the images. Comparing the experimental results, we can find that Algorithm 2 reduces the artifacts more effectively than PPB-AA method. At the same time, the staircasing effects also decrease in these two methods.

e Influence of Parameters and Mean.
Here, we study the sensitivities of parameters and the effect of mean of our restoration model. According to the analysis in Section 2.1, the value of parameter β varies with the level of the noise. Experimentally, β is chosen as 1 when L > 5. When L � 5, β is set as 1.1. In the proposed model, parameter α needs to satisfy condition α ≥ (1/12). Parameter λ is set as λ 2 < (1/72) based on the description in Section 2.3, which is changing from 0.01 to 0.11 with an interval 0.01 in experiment. In the denoising case, we take the noisy images "Lena" and "Cameraman" as examples. For each image, we vary the values of α with three levels of the noise, L � 10, 7, 5. e PSNR values of restored images using Algorithm 1 are presented in Figure 17. We observe that there is an upmost point on each curve. With the increasing α, the performance of "Lena" with L � 5 and "Cameraman" with L � 10 decreases rapidly, which means that the denoising results are sensitive to value of α for each L. In the simultaneous deblurring and denoising case, we take the noisy image of "Cameraman" blurred by the motion blur with length 5 and angle 30 and then corrupted by multiplicative Nakagami noise with L � 10 as an example. e parameter settings are λ � 0.01, α � 7, and β � 1 in this case known from Table 3. At first, for the case of λ � 0.01, we study the sensitivity of the parameter α in terms of solution quality and computation time. e PSNR (dB), SSIM, and CPU time (s) used for various values of α are listed in Table 6. We observe that if α is less than 7, the PSNR value of the restored image by Algorithm 1 increases with the increase in α. However, the PSNR and SSIM values decrease as the value of α increases when α is more than 7. e numerical results indicate that there is a certain value of α to make the model get the best quality of recovered image for a fixed value of λ. en, for the case of α � 7, we further analyze the sensitivity of the parameter λ. e PSNR (dB), SSIM, and CPU time (s) used     Table 7. From Table 7, we observe that λ � 0.01 is the best choice for obtaining the highest PSNR value of denoised image. Finally, in order to discuss the influence of preserving the mean of the original image, there are four assumptions, including recovered images obtained from model (9) and model (9) with preserving 1 times, 1.2 times, 0.8 times the mean of the original image, respectively. Table 8 reports the PSNR (dB) and SSIM values on Cameraman image with different mean. We find the fact the PSNR (dB) and SSIM (d) (e) (f )   values obtained from model (9) with preserving mean of the original image are almost higher than the others. As shown in Figure 18, when the mean of recovered image is larger than that of the original image, the restored image by Algorithm 1 looks lighter. In contrast, the restored image is darker. PSNR and SSIM values of both decrease significantly. Although there is not much difference between the recovered images obtained from preserving the mean and not preserving the mean, the PSNR value is higher when the mean is maintained as the original. ese indicate that bias correction technique [36] is helpful to enhance the performance of the model (9).

Conclusion
is paper proposes a new variational model for solving the SAR image restoration problem by investigating the statistical property of multiplicative noise from Nakagami distribution. Under a mild condition, we show that the new model is strictly convex. We further analyze some mathematical properties of the model in detail. We adopt the primal-dual algorithm to solve the convex optimization problem. Two kinds of simulated experiments and real SAR images experiments are conducted. e experimental results demonstrate the excellent performance of the proposed method over some state-of-the-art methods in terms of the noise removal capability and detail preservation capability. Appendix e ideas for proving Propositions 1, 2, and 4 and eorem 1 are similar to those in [36,39]. For readers' convenience, we present them in detail as follows.

Mathematical Problems in Engineering
For each x ∈ Ω, we denote u(x)/f(x) by t, and we can obtain the strict convexity of the first two terms in (6) based on the analysis above and the fact that log f is a constant. In view of the convexity of the TV regularization, we can see that E(u) in (6) is strictly convex when α ≥ (1/12). In addition, the feasible set S(Ω) is convex and we can conclude that the model (6) is strictly convex when α ≥ (1/12). □

B. Proof of Theorem 1
Proof. We first prove that E(u) in (6) is bounded from below. By the monotonic property of the function g(s) � 2 log s + (1/s 2 ) on R + , we have 2 log t + (f(x)/t) 2 ≥ 1 + 2 log f(x) for any t ∈ R + ∪ 0 { } and x ∈ Ω. It then follows that which implies that E(u) is bounded from below. We next show that the solution of model (6) exists and is bounded. Let a � inf Ω f and b � βsup Ω f. We define

Data Availability
All data generated or analysed during this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest.