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

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

在線性不變的空間中，退化圖像的復(fù)原是一個(gè)求逆的過程。根據(jù)相關(guān)的點(diǎn)擴(kuò)散函數(shù)（簡稱PSF）的先驗(yàn)知識(shí)是否已知，將退化圖像的復(fù)原分為退化圖像的盲復(fù)原和退化圖像的非盲復(fù)原[1]。非盲復(fù)原就是在退化圖像相關(guān)的先驗(yàn)知識(shí)已知的情況下，進(jìn)行解卷積運(yùn)算。而圖像盲復(fù)原就是在PSF未知的情況下，利用盲復(fù)原算法實(shí)現(xiàn)退化圖像的清晰化。圖像的退化模型如圖1所示。可以用下面的數(shù)學(xué)表達(dá)式來描述退化過程：

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

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

圖1 圖像的退化模型

該盲復(fù)原算法在估計(jì)出點(diǎn)擴(kuò)散函數(shù)之后，直接進(jìn)行反卷積，表達(dá)式如下：

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

1 近似GCD圖像盲復(fù)原算法

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

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

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

假設(shè)原圖像為f(x,y)，兩個(gè)退化圖像為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為多項(xiàng)式中各個(gè)變量的最高次數(shù)），兩者的近似GCD的支撐區(qū)域較大，可以通過計(jì)算多項(xiàng)式G1，G2的近似GCD實(shí)現(xiàn)圖像盲復(fù)原。對于退化的RGB圖像，可以對R、G、B三個(gè)通道分別描述，如下所示：g2(x,y)進(jìn)行二維Z變換得到下面的表達(dá)式：

1.1 單變量多項(xiàng)式近似GCD算法

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

其中，Bezout矩陣定義如下：

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

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

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

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

（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）假設(shè)y=[1,y,y2,...,ym-1]T滿足Cy=b，其中C=B(f1,f2)m-r，-b是由矩陣B(f1,f2)的第m-r+1個(gè)列向量的前m-r列構(gòu)成。使：

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

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

1.2 二元多項(xiàng)式近似GCD算法

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

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

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

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

2 近似GCD算法的改進(jìn)

改進(jìn)算法主要分為兩個(gè)步驟來實(shí)現(xiàn)：近似GCD算法求解圖像多項(xiàng)式的輔助因子，全變分正則化迭代解卷積。

2.1 近似GCD算法求解圖像多項(xiàng)式的輔助因子

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

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

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

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

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

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

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

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

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

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

矩陣范數(shù)的導(dǎo)數(shù)有下面的性質(zhì)：

可以通過式（19）和（21）對各個(gè)變量進(jìn)行求解，步驟如下所示：

（1）對w的求解

對w的求解分為兩種情況：當(dāng)φ=2時(shí)，可得到w的表達(dá)式如下：

當(dāng)φ=1時(shí)，可得到w的表達(dá)式如下：

（2）對z的求解

（3）對f的求解

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

2.3 改進(jìn)算法的流程

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

圖2 高斯噪聲水平為10-4時(shí)的復(fù)原結(jié)果

圖3 高斯噪聲水平為5×10-4時(shí)的復(fù)原結(jié)果

表1 兩種算法復(fù)原結(jié)果的PSNR（高斯噪聲）

表2 兩種算法復(fù)原結(jié)果的SSIM（高斯噪聲）

圖4 高斯噪聲水平為10-3時(shí)的復(fù)原結(jié)果

圖5 均勻分布噪聲水平為10-4時(shí)的復(fù)原結(jié)果

表3 兩種算法復(fù)原結(jié)果的PSNR（均勻分布噪聲）

表4 兩種算法復(fù)原結(jié)果的SSIM（均勻分布噪聲）

圖6 均勻分布噪聲水平為5×10-4時(shí)的復(fù)原結(jié)果

圖7 均勻分布噪聲水平為10-3時(shí)的復(fù)原結(jié)果

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

4 結(jié)論

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