基于超拉普拉斯先驗的分數階全變分圖像復原模型

李振

摘要：針對點擴散函數形式未知的模糊圖像復原的問題，提出一種基于超拉普拉斯先驗的分數階全變分正則模型。為了消除的整數階全變分模型中圖像復原中的振鈴效應， 本文引入了分數階差分來作為圖像的先驗信息，在基于最大后驗概率估計（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.

【通聯編輯：光文玲】