基于超拉普拉斯先验的分数阶全变分图像复原模型

李振

摘要：针对点扩散函数形式未知的模糊图像复原的问题，提出一种基于超拉普拉斯先验的分数阶全变分正则模型。为了消除的整数阶全变分模型中图像复原中的振铃效应， 本文引入了分数阶差分来作为图像的先验信息，在基于最大后验概率估计（MAP）的图像盲去模糊模型的基础上，提出了一种[Lp]范数分数阶正则项的全变分复原模型，利用分裂Bregman算法交替求解得到对点扩散函数的估计，再利用非盲复原的方法求得原始清晰图像。实验结果表明，此模型能够有效地处理不同形态的模糊核，得到了较为理想的复原结果，复原图像中的振铃效应有效地降低。

关键词：超拉普拉斯先验;全变分正则模型;分数阶微积分;图像盲去卷积

中图分类号：TP311      文献标识码：A

文章编号：1009-3044（2020）25-0043-0c

20世纪末，数字图像技术被广泛应用于包括医疗诊断、航空遥感、工业检测、显微成像、军事用途和公共交通中。数字图像在产生、采集、传输、处理的过程中，会受到各种因素的影响，导致得到的图像与原始图像不符，这一过程我们称之为图像的退化，需要通过各种形式上的数学方法从退化图像的先验信息和退化方式恢复出原始图像。在这其中本文的研究重点是基于稀疏先验正则化的模糊图像盲复原。

图像复原是一种线性病态问题，无论点扩散函数是否已知，复原图像的计算在数学上都是不适定的。在某些特定的情形中，图像的退化形式可以被先验地假设为参数化的点扩散函数并加以求得。但是在更广泛的应用场景中，图像的点扩散函数很难被获取到。所以图像的盲复原有很大的研究价值。

在减轻盲反卷积病态性的过程中，模糊核先验和图像先验常常作为正则项加入图像盲去模糊模型中。吉洪诺夫（Tikhonov）正则化方法[1]基于平滑准则，采用了[H1Ω]中的半范数平方作为正则项，也就是图像梯度的[L2]范数，将病态的求解过程转换为稳定的良态过程，所以吉洪诺夫正则化成为图像去噪中常用的工具：

该方法在有效消除噪声的同时也很好地保护了图像边缘，因此在此基础上人们又做了大量的研究。You[3]提出了扩散率固定点迭代数值算法，在迭代求解全变分模型时近似线性化;Chan[4]具有很强的鲁棒性。基于L1范数提出了全变分正则化的盲复原模型，联合求解了清晰图像和点扩散函数。Fergus[5]等提出的复原算法利用零均值混合高斯分布（Mixture of Gaussians）來模拟自然场景图像梯度域中的重尾分布，用混合指数作为模糊核先验知识，对相机抖动模糊问题取得了良好的效果。Shan[6]连续分段函数模拟自然场景梯度图的重尾分布先验，同时求解模糊核和清晰图像，抑制了振铃效应对模糊核估计的影响。Levin[7]将图像盲复原分成清晰图像估计和模糊核估计两个独立的子过程，证明了算法的有效性。Krishnan和Fergus[8]将图像的超拉普拉斯先验（ Hyper-Laplacian） 作为正则项，得到的结果保留了相当的高频信息，同时一定程度上抑制了振铃效应。Kotera[9]同样采用超拉普拉斯先验对图像进行建模，为了模糊核的稀疏性，采用L1 范数作为模糊核的先验信息. Krishnan等人[10]在稀疏表示压缩感知领域提出一种L1 /L2 范数正则化模型，并用迭代收缩阈值算法（Iterative Shrinkage-Thresholding Algorithm，ISTA）进行求解。Goldstein和Li[11]等人将图像复原问题转化为鞍点问题，提出自适应PDHG算法加以求解。同时期Chen[12]结合全变分正则先验和分数阶微分提出了一种新的分数阶总变分正则化函数：

而图像的复原过程就是通过观测到的模糊图像[g]，对形式未知或者可以参数化的模糊核[h]进行估计，最终得到清晰图像[u]。由式（5）可以看到，图像复原过程可以视作清晰图像退化过程的逆问题，所有的反问题都有不同程度的病态性（Ill-posed problem），而正则化模型则是在病态性问题解决上的比较好用的框架。

从Bayesian统计框架的角度来看，清晰图像[u]和模糊核[h]的同时求解相当于求解标准最大后验（Maximum A Posteriori，MAP）估计：

惯常使用的图像灰度值梯度是在微小邻域内定义的（通常是二方向或四方向的一阶差分），容易造成复原图像整体呈现“阶梯效应[17]”，而在图像复原的过程中我们非常需要一种具有长时记忆的灰度值梯度的度量方法来描述图像的中频信息来解决阶梯效应。为此考察微积分理论中高阶微分[18]与分数阶微分[19]两个数学工具。高阶微分可以减少图像复原中的阶梯效应，同时会造成图像地过平滑，破坏图像的边界。分数阶微分，可以视作整数阶微分在复数域上的自然延拓，在数值计算上具有记忆性和全局性。

前人给出了不同角度给出了分数阶微积分的定义，包括基于广义函数的分数阶微积分定义、Caputo分数阶微积分定义，R-L（Riemann-Liouville）分数阶微积分定义以及G-L（Gru?wald-Letnikov）分数阶微积分定义这几类。本文选择Gru?wald-Letnikov定义[α]阶微分定义：

式（12）中第一项是数据保真项，用于保证求解时[u*h]的正确性和图像[u]解的范围，[Q（u）]和[R（h）]分别为空域中对图像[u]的分数阶[Lp]范数约束以及对模糊核[h]的非负[Lp]范数约束，[（i，j）]为图像像素坐标二元组，[γ]为正则化参数，用于平衡数据保真项和正则项。

2模型求解

一般来说，超拉普拉斯先验可以为正则化模型带来更好的反问题解效果，但是令模型本身变成高度非凸的。并且如果不考虑先求解点扩散函数再通过非盲去卷积求得清晰图像的策略，则模型求解目标为清晰图像和点扩散函数两个子问题。Goldstein和Osher[20]在2009年提出的分裂Bregman算法是近年来图像处理领域里运用广泛的迭代方法，可以看作是缩放版的交替乘子算法， 收敛速度快，在求解过程中的正则参数保持不变。为了从数值上找到[u]，[h]的解，借鉴文献[20]利用分裂Bregman（Split-Bregman）迭代的思路，在求解过程中交替地固定一个变量，最小化另外一个变量，将模型（12）分解为[u]步骤和[h]步骤两个子问题。在每个最小化子问题中，使用增广拉格朗日乘子法（ALM）进行求解。主迭代完成后，利用求解得到的模糊核通过非盲去卷积得到最终的复原清晰图像。

2.1U步骤

在更新图像[u]变量的过程中，我们用上一次迭代中求得的模糊核结果[H]固定另外一个待求变量[h]。在最小化問题的求解中忽略常量，此时U步骤子问题整理为：

3 实验结果和结论

我们在文献[7]提供的数据集上测试了我们模型的可计算型，该数据集由四个[255 × 255]灰度清晰图像和八个尺寸为[13 × 13 ～ 27 × 27]真实运动模糊的PSF组成，生成了32个测试降质模糊图像。实验给出了本文和其他几种方法的图像恢复结果及估计PSF，以及恢复图像对应原始图像的峰值信噪比（ PSNR） 。此外，所有实验均在MacBook Pro （15-inch， 2018）笔记本上完成，运行macOSMojave10.14.6系统以及MATLAB 2019a计算平台。

在式（15）和式（21）中各参数选取为[γ= 100，α= 100，β=1， p = 0.5]，分数阶[α=1.25]，主迭代和u、h步骤的迭代次数分别为10次，得出结果并和近年来盲去模糊较为经典和效果比较好的算法进行比较，图3是以Painting+5号模糊核为例各模型的计算结果。

在基于MAP 估计的图像盲复原框架下，本文提出了一种基于拉普拉斯先验分数阶变分的图像盲复原算法。针对整数阶微分在图像处理中信息利用过少的问题，引入了分数阶微积分的处理方式，同时对待估模糊核L1稀疏正则项进行约束，优化求解过程中通过多尺度处理框架下的分裂Bregman法进行处理。在实验中本文提出的算法和其他算法相比，不仅复原效果更好，计算速度也得到了提高。

参考文献：

[1] Tikhonov， Andrej Nikolaevich， Arsenin， VasiliyYakovlevich. Solutions of ill-posed problems[J]. mathematics of computation， 1977， 32（144）：491-491.

[2] Rudin L I，OsherS，FatemiE.Nonlinear total variation based noise removal algorithms[J].Physica D：NonlinearPhenomena，1992，60（1/2/3/4）：259-268.

[3] You Y L，Kaveh M.A regularization approach to joint blur identification and image restoration[J].IEEE Transactions on Image Processing，1996，5（3）：416-428.

[4] Chan T F，Wong C K.Total variation blind deconvolution[J].IEEE Transactions on Image Processing，1998，7（3）：370-375.

[5] Fergus R，SinghB，HertzmannA，etal.Removing camera shake from a single photograph[J].ACM Transactions on Graphics，2006，25（3）：787-794.

[6] Shan Q，Jia J Y，AgarwalaA.High-quality motion deblurringfrom asingle image[J].ACM Transactions on Graphics，2008，27（3）：557-566.

[7] Levin A ， Weiss Y ， Durand F ， et al. Understanding and evaluating blind deconvolution algorithms[C]// 2009 IEEE Conference on Computer Vision and Pattern Recognition （CVPR）. IEEE， 2009.

[8] Krishnan D ， Fergus R . Fast Image Deconvolution using Hyper-Laplacian Priors[C]// Advances in Neural Information Processing Systems 22： Conference on Neural Information Processing Systems A Meeting Held December. Curran Associates Inc. 2009.

[9] KoteraJ，?roubekF，MilanfarP.Blinddeconvolutionusing alternating maximum a posteriori estimation with heavy-tailed priors[J].Computer Analysis of Images and Patterns，2013.

[10] Krishnan D，BrunaJ，FergusR.Blinddeconvolution with non-local sparsityreweighting[J].2013.

[11] Goldstein T，Li M，Yuan X ，et al.Adaptive Primal-Dual Hybrid Gradient Methods for Saddle-Point Problems[J].2013.

[12] Chen D L，Sun S S，ZhangCR，etal.Fractional-order TV-L2 model for image denoising[J].Central European Journal of Physics，2013，11（10）：1414-1422.

[13] XuL ，Zheng S ， Jia J . Unnatural l0 sparse representation for natural image deblurring[C]// 2013 IEEE Conference on Computer Vision and Pattern Recognition. IEEE， 2013.

[14] XuL，Lu C W，XuY，etal.Imagesmoothing via L0 gradient minimization[J].ACM Transactions on Graphics，2011，30（6）：1-12.

[15] ChartrandR ， Yin W . Iteratively reweighted algorithms for compressive sensing[C]// Acoustics， Speech and Signal Processing， 2008. ICASSP 2008.IEEE International Conference on.IEEE， 2008.

[16] MazumderR，Friedman J H，HastieT.SparseNet：coordinatedescent with nonconvexpenalties[J].Journal of the American Statistical Association，2011，106（495）：1125-1138.

[17] GuichardF，MoisanL，Morel J M.A review of P.D.E.models in image processing and image analysis[J].Journal De Physique IV （Proceedings），2002，12（1）：137-154.

[18] Chan T F，Shen J .Image Processing and Analysis （Variational， PDE， Wavelet， and Stochastic Methods） || 2. Some Modern Image Analysis Tools[J]. 2005， 10.1137/1.9780898717877：31-89.

[19] ZhangJ，WeiZH.A class of fractional-order multi-scale variational models and alternating projection algorithm for image denoising[J].Applied Mathematical Modelling，2011，35（5）：2516-2528.

[20] Goldstein T，OsherS.Thesplit bregman method for L1-regularized problems[J].SIAM Journal on Imaging Sciences，2009，2（2）： 1-21.

[21] ChartrandR.Nonconvexsplitting for regularized low-rank + sparse decomposition[J].IEEE Transactions on Signal Processing，2012，60（11）：5810-5819.

[22] She Y Y.An iterative algorithm for fitting nonconvex penalized generalized linear models with grouped predictors[J].Computational Statistics & Data Analysis，2012，56（10）：2976-2990.

[23]Emmanuel J.Candès，Wakin M B，Boyd S P.Enhancing Sparsity by Reweighted L1Minimization[J].Journal of Fourier Analysis & Applications， 2007，14（5）：877-905.

[24] Levin A ， Weiss Y ， Frédo Durand， et al. Efficient Marginal Likelihood Optimization in Blind Deconvolution[C]// The 24th IEEE Conference on Computer Vision and Pattern Recognition， CVPR 2011， Colorado Springs， CO， USA， 20-25 June 2011. IEEE， 2011.

【通聯编辑：光文玲】