正則化的近似最大公因子的圖像盲復原算法

張萬磊，楊陽，馮濤，李喆

（1.長春理工大學 電子信息工程學院，長春 130022；2.長春理工大學 理學院，長春 130022）

在線性不變的空間中，退化圖像的復原是一個求逆的過程。根據相關的點擴散函數（簡稱PSF）的先驗知識是否已知，將退化圖像的復原分為退化圖像的盲復原和退化圖像的非盲復原[1]。非盲復原就是在退化圖像相關的先驗知識已知的情況下，進行解卷積運算。而圖像盲復原就是在PSF未知的情況下，利用盲復原算法實現退化圖像的清晰化。圖像的退化模型如圖1所示。可以用下面的數學表達式來描述退化過程：

式中，g(x,y)代表退化圖像，f(x,y)代表原始圖像，h(x,y)代表圖像退化過程中的PSF，n(x,y)代表獲取圖像過程中疊加的噪聲，“*”代表卷積運算。

由于在盲復原過程中缺少一些先驗信息，退化圖像中又存在噪聲，這使得圖像的盲復原不再是一個簡單的逆濾波過程。1997年，Liang和Pillai[2]提出了近似最大公因子圖像盲復原算法，該算法把原始圖像看成多幅退化圖像的GCD，將二維GCD問題轉化為一維Sylvester型GCD算法，進而實現圖像盲復原。2010年，支麗紅等人[3]提出了快速近似GCD圖像盲復原算法，降低了盲復原算法的計算量。2015年，Belhaj等人[4]提出了基于Hankel矩陣計算GCD的圖像盲復原算法，將單變量GCD算法運用在連續變換矩陣為上三角Toeplitz矩陣中。2017年，Toca等人[5]提出了基于Bezout矩陣的圖像盲復原算法，該算法使用Barnett方法的廣義Bezout矩陣來求解GCD。

圖1 圖像的退化模型

該盲復原算法在估計出點擴散函數之后，直接進行反卷積，表達式如下：

即直接逆濾波的算法，其中IFFT表示逆快速傅立葉變換，G代表退化圖像的傅立葉變換，代表估計的PSF的傅立葉變換，(x,y)代表復原的圖像。在平面上有些區域上為零或者非常小，F上相應區域的數值會劇烈地變化，這些區域的噪聲會被放大，產生較大的偏差。其次，該算法缺乏約束條件，不能有效抑制因噪聲擾動和數值計算產生的誤差，誤差會不斷累積并向下傳遞。以上這些因素都會影響復原圖像的質量。針對上面的問題，本文提出了相關的改進算法，從抑制噪聲和反卷積運算約束兩個方面去改進近似GCD圖像盲復原算法，提高算法的魯棒性。

1 近似GCD圖像盲復原算法

定義1.二維Z變換將m×n矩陣P的元素映射到雙變量多項式p(x,y)的系數[3]：

其中，x=[1,x,x2,...,xm-1]T，y=[1,y,y2,...,yn-1]T。根據上面的定義可以把式（1）變換成：

式中，G，F，H，N分別代表g(x,y)，f(x,y)，h(x,y)，n(x,y)的二維Z變換。

假設原圖像為f(x,y)，兩個退化圖像為g1(x,y)，g2(x,y)。兩者的PSF分別為h1(x,y)，h2(x,y)，加性噪聲分別為n1(x,y)，n2(x,y)，把g1(x,y)，

通常來說deg(H1)和deg(H2)與deg(F)相比非常小，（其中deg為多項式中各個變量的最高次數），兩者的近似GCD的支撐區域較大，可以通過計算多項式G1，G2的近似GCD實現圖像盲復原。對于退化的RGB圖像，可以對R、G、B三個通道分別描述，如下所示：g2(x,y)進行二維Z變換得到下面的表達式：

1.1 單變量多項式近似GCD算法

假設給出兩個單變量多項式f1,f2∈C[x]{0}，deg(f1)=m和deg(f2)=n，m≥n，

其中，Bezout矩陣定義如下：

式中，J是單位反對角矩陣。

定理1.給定單變量多項式f1,f2∈C[x]，其中deg(f1)=m，deg(f2)=n，m≥n，那么有如下結論：

其中，dim代表矩陣的維度，NullSpace代表矩陣的零空間[6]。

定理2.給定單變量多項式f1(x)，f2(x)，deg(f1)=m，deg(f2)=n，m≥n。令p(x)=gcd(f1,f2)與deg(p)=r，結論如下所示：

（1）對于k≤m-r，rank(B(f1,f2))=m-r和det(B(f1,f2)k)，其中B(f1,f2)k是B(f1,f2)的子矩陣，但對于k＞m-r，det(B(f1,f2)k)=0。

（2）假設y=[1,y,y2,...,ym-1]T滿足Cy=b，其中C=B(f1,f2)m-r，-b是由矩陣B(f1,f2)的第m-r+1個列向量的前m-r列構成。使：

則f1(x)=p(x)u(x)，其中。

根據定理2，可以通過檢查前1×1，2×2，4×4，...，2[log2(m-r+1)]×2[log2(m-r+1)]子矩陣是否是奇異矩陣，從而有效地估計gcd(f1,f2)的次數r。假設r=deg(gcd(f1,f2))，可以通過求解 (m-r)×(m-r)線性系統來計算輔因子u(x)，同時使用基于快速傅立葉變換的近似多項式除法來計算出GCD。顯然，只需要B(f1,f2)的m-r+1階子式，這樣可以在計算GCD的時候節省時間和空間。

1.2 二元多項式近似GCD算法

有兩個二元多項式f1(x,y)和f2(x,y)，假設degx(f1)=degx(f2)=m，degy(f1)=degy(f2)=n，其中degx(f1(x,y))表示f1(x,y)中y為常量時，多項式中變量x的最高次數，degy(f1(x,y))同理。先將和中。每個k對應相應的兩個一元多項式，將前一節的一元近似GCD算法應用到這些一元多項式中，可以得到標量[2]。 將n-1代入到上述的式子中，可得到離散傅立葉變換元素的矩陣，如下所示：

其中，A(k,l)是的GCD估計值。同理可以將代入f1(x,y)和f2(x,y)中，得到單變量多項式的GCD進一步代入，然后得到另一個離散傅立葉變換元素的矩陣：

向量a(k)和b(l)可以通過最小二乘法的方法解得，如下所示：

可以通過用逆離散傅立葉變換來計算p(x,y)=gcd(f1,f2)。

2 近似GCD算法的改進

改進算法主要分為兩個步驟來實現：近似GCD算法求解圖像多項式的輔助因子，全變分正則化迭代解卷積。

2.1 近似GCD算法求解圖像多項式的輔助因子

假設原始圖像g(x,y)產生退化圖像g1(x,y)和g2(x,y)，對于退化圖像多項式來說，一般其PSF多項式各個變量的最高次數較小，可以把g1(x,y)和g2(x,y)的近似GCD多項式的次數看成m和n。用前文引入的二元近似GCD算法，內插輔助因子h1(x,y)或者h2(x,y)，即兩幅退化圖像的PSF，然后通過最小二乘法得到近似PSF。具體的算法流程如下所示：

輸入：degx(g1)=degx(g2)=m，degy(g1)=degy(g2)=n，其中g1(x,y)，g2(x,y)∈?[x,y]，ε∈R＞0，其中ε為給定的容錯度；

輸出：g1(x,y)和g2(x,y)的近似輔助因子h1(x,y)和h2(x,y)，其中h1(x,y),h2(x,y)∈?[x,y]。

步驟1：估計出r=degx(h)和s=degy(h)，即f(x,y)中兩個變量的最高次數；

步驟2：用單變量多項式近似GCD算法去計算矩陣 [h(xk,yl)]∈ ?(m-r+1)×(n-s+1)，其中：

步驟3：h(xk,yl)通過逆快速傅里葉變換計算出h(x,y)，即圖像多項式的輔助因子。

2.2 全變分正則化迭代解卷積

上述的近似GCD圖像盲復原算法中直接使用逆濾波算法去求解復原圖像，會將圖像中有些區域的噪聲放大，嚴重影響最后盲復原的結果[8]。為了解決該問題，在近似GCD算法求解出h(x,y)之后，利用全變分正則化進行非盲解卷積，解出復原圖像。該算法在解卷積的同時可以更好地抑制圖像中的噪聲，保存好圖像的邊緣細節部分，為此引入以下的解卷積約束模型，如下所示：

由于式（16）中點擴散函數h已知，可以根據半二次規整化理論，引入輔助變量z，w將上式轉換成易求解的形式，按照w→z→f的順序來交替最小化[11]。

其中，λ，β為新設置的參數，φ的定義見式（17），滿足λ,β→∞ ，如果‖z- (h?f-g)‖2→0 不成立，這與上式最小化相矛盾，‖z- (h?f-g)‖2→0必成立。同理，‖w- ?f‖2→0必成立。對式（18）進行求偏導，可得：

矩陣范數的導數有下面的性質：

可以通過式（19）和（21）對各個變量進行求解，步驟如下所示：

（1）對w的求解

對w的求解分為兩種情況：當φ=2時，可得到w的表達式如下：

當φ=1時，可得到w的表達式如下：

（2）對z的求解

（3）對f的求解

其中，F-1表示傅里葉逆變換，“*”表示矩陣的共軛，H,G,Z,D,W分別表示h,g,z,?,w的傅里葉變換。

2.3 改進算法的流程

針對原算法對噪聲敏感的問題，改進的算法采用全變分正則化對圖像盲復原的過程進行約束。該算法中使用兩張退化圖像g1、g2，使用近似GCD算法估算出近似PSF后，然后，用全變分正則化迭代解卷積，當滿足停止迭代條件‖fn+1-fn‖＞?時，得到復原圖像。算法流程如下：

輸入：模糊圖像g1、g2；Step1：估計出原圖像多項式兩個變量的最高次數；Step2：計算矩陣h(xk,yl)，即點擴散函數的傅立葉變換形式；Step3：h(xk,yl)進行逆傅立葉變換，得到近似PSF；Step4：利用式(22)或者(23)對w進行更新；Step5：利用式(24)對z進行更新；Step6：利用式(25)對f進行更新；Step7：判斷是否滿足‖ ‖fn+1-fn＞?，如果滿足執行下一步，否則返回Step4；輸出：清晰圖像 fn+1。

3 實驗結果與分析

本文的實驗環境為：處理器為Intel Core i5 3210M，操作系統為Windows7 64位操作系統，運行內存8G，仿真環境為Matlab R2016b。改進算法用全變分正則化算法取代原算法中的逆濾波算法，實驗中閾值Th=1.2δ，δ為圖像中局部的估計的標準差。為了驗證改進算法的性能，下面設置了對比試驗，對原始圖像進行模糊處理，其中原始圖像Lena和Peppers大小為256×256模糊核的大小為7×7，并施加不同水平的噪聲，下面對Lena施加高斯噪聲，對Peppers施加均勻分布噪聲，用PSNR和SSIM來衡量復原圖像的質量。如圖2-圖7給出不同情況下復原的結果，表1-表4給出相應的PSNR和SSIM。下面實驗中GCD算法為文獻[3]中支麗紅等人提出的算法。

圖2 高斯噪聲水平為10-4時的復原結果

圖3 高斯噪聲水平為5×10-4時的復原結果

表1 兩種算法復原結果的PSNR（高斯噪聲）

表2 兩種算法復原結果的SSIM（高斯噪聲）

圖4 高斯噪聲水平為10-3時的復原結果

圖5 均勻分布噪聲水平為10-4時的復原結果

表3 兩種算法復原結果的PSNR（均勻分布噪聲）

表4 兩種算法復原結果的SSIM（均勻分布噪聲）

圖6 均勻分布噪聲水平為5×10-4時的復原結果

圖7 均勻分布噪聲水平為10-3時的復原結果

從實驗結果可以看出文獻[3]中近似GCD算法對噪聲比較敏感，在噪聲水平為10-4時，退化圖像其復原效果不錯，能夠保持圖像中的很多細節信息。但是，隨著噪聲水平的增大，其復原圖像的質量逐漸下降，圖像中的細節變得模糊，同時還出現很多波紋現象，嚴重降低了圖像的視覺效果。這是由于噪聲擾動產生的誤差，在缺乏約束條件和直接逆濾波的作用下被放大。采用改進算法后，可以明顯看出復原圖像的視覺效果得到改善，抗噪性得到提高，可以從表1、表3中的PSNR和表2、表4中的SSIM可以看出改進算法具有較大的優勢，在施加兩種不同噪聲的情況下，改進算法的PSNR提高了1～5dB，SSIM提高了0.09～0.3。

4 結論

本文采用全變分約束解卷積來改進近似GCD圖像盲復原算法。近似GCD圖像盲復原算法魯棒性差的原因在于，缺乏約束條件和后端使用逆濾波算法。在盲復原過程中，缺乏約束條件使數值誤差往下傳遞。其次，算法的后端處理使用了逆濾波算法，實驗表明該步驟會放大圖像局部的噪聲。全變分非盲解卷積不僅可以抑制退化圖像中的噪聲，而且還可以在保存圖像細節的同時完成解卷積。從仿真實驗結果表明，與近似GCD圖像盲復原算法相比，改進算法無論在視覺效果還是在PSNR和SSIM上均有所改進。