（浙江理工大學(xué) 理學(xué)院，浙江 杭州 310018）

0 引言

周期性噪聲和準(zhǔn)周期性噪聲可能源于圖像采集與傳輸過程中的機(jī)電干擾，也可能來(lái)自電子干擾，例如現(xiàn)代數(shù)碼相機(jī)的電子電路和電荷耦合器件等。當(dāng)多成像系統(tǒng)中檢測(cè)器的靈敏度不匹配時(shí)，經(jīng)常會(huì)產(chǎn)生一種特殊的條紋狀擬周期性噪聲，這種噪聲被稱為摩爾紋（moiré）噪聲。摩爾紋噪聲經(jīng)常出現(xiàn)在磁共振成像、顯微鏡圖像及遙感圖像等實(shí)時(shí)圖像的應(yīng)用過程中［1］。由于摩爾紋噪聲在空間上是相關(guān)的，并散布在整個(gè)圖像中，因此使用空間域?yàn)V波算法去除摩爾紋噪聲可能會(huì)在實(shí)體周圍產(chǎn)生偽影，有時(shí)還會(huì)使恢復(fù)的圖像變得模糊。但在傅里葉變換圖像頻譜中，可以很容易觀測(cè)到由摩爾紋噪聲引起的高振幅頻譜分量。因此，通過頻域操作恢復(fù)這些頻譜分量能夠抑制噪聲帶來(lái)的影響，從而恢復(fù)原始圖像信息。

為去除摩爾紋噪聲，可采用的一種方法是使用頻域?yàn)V波器。如Moallem 等［2］提出一種新穎的自適應(yīng)高斯陷波濾波器，可有效降低周期與準(zhǔn)周期噪聲；Varghese［3］使用交換中值濾波器（Switching Median Filter）去除高密度的周期噪聲。在此基礎(chǔ)上，Charkraborty 等［4-5］利用雙閾值法確認(rèn)波峰位置，并使用自適應(yīng)高斯恢復(fù)濾波器（Adaptive Gaussian Restoration Filter）去除周期噪聲，并利用自動(dòng)陷波抑制濾波器對(duì)Gabor 變換后圖像的特定周期信號(hào)進(jìn)行濾波，從而去除對(duì)應(yīng)周期區(qū)間內(nèi)的準(zhǔn)周期噪聲。在醫(yī)學(xué)領(lǐng)域，Ionita 等［6］采用自適應(yīng)閾值，能更好地確定受周期噪聲影響的傅里葉頻域波峰檢測(cè)閾值，然后使用窗口高斯陷波濾波器去除出現(xiàn)在顯微鏡圖像中的周期噪聲。此外，Ionita 等［7］通過分析一階Haar 小波變換頻域，確定了與噪聲對(duì)應(yīng)的頻率區(qū)域，然后通過校正顯微鏡圖像的傅里葉幅度譜分量去除噪聲。另一種方法是使用最小化模型去除摩爾紋噪聲。如Chang等［8］提出一種針對(duì)遙感影像的圖像恢復(fù)算法，該算法首先對(duì)條紋信號(hào)進(jìn)行分類，然后利用全變差（Total variation）模型去除噪聲信號(hào)，以保持影像光譜與空間的一致性。由于該模型需要目標(biāo)圖像的特征梯度小于條紋梯度，因此對(duì)于一般的自然圖像，使用該算法可能會(huì)破壞原始圖像邊緣信息。

上述研究雖然能較為有效地去除摩爾紋噪聲，但恢復(fù)的傅里葉頻域信號(hào)依賴于其領(lǐng)域信號(hào)。因此，若圖像在傅里葉低頻區(qū)域內(nèi)受噪聲影響的系數(shù)范圍較大時(shí)，上述研究無(wú)法準(zhǔn)確地恢復(fù)原始圖像低頻信息。為更好地恢復(fù)原始圖像低頻信息，本文提出一種基于緊小波框架的圖像恢復(fù)算法。算法本質(zhì)是利用緊框架下圖像小波系數(shù)的稀疏性特征，根據(jù)壓縮感知原理建立優(yōu)化模型，通過求解該模型最優(yōu)解得到原始圖像稀疏恢復(fù)結(jié)果。該算法有助于在機(jī)電干擾較為嚴(yán)重的情況下進(jìn)行圖像采集。

1 圖像恢復(fù)原理

1.1 稀疏信號(hào)恢復(fù)

壓縮感知理論表明，利用目標(biāo)（如信號(hào)、圖像等）的稀疏性，可使用遠(yuǎn)遠(yuǎn)少于香農(nóng)采樣定理（Shannon-Nyquist Sampling Theorem）所需的樣本數(shù)重構(gòu)原始信號(hào)。以帶高斯噪聲的信號(hào)恢復(fù)問題為例，若n維列向量x為需要恢復(fù)的原始稀疏信號(hào)，則實(shí)際信號(hào)采集過程可表示為：

其中，y∈?m為采集到的信號(hào)，A∈?m×n為測(cè)量矩陣，w∈?m為獨(dú)立同分布的高斯噪聲項(xiàng)。因此，稀疏信號(hào)恢復(fù)問題可轉(zhuǎn)換為一個(gè)線性方程組的求解問題。在實(shí)際問題中，矩陣A的條件數(shù)通常較大，因此無(wú)法直接采用最小二乘法。一種可行的方法是借助原始信號(hào)x的稀疏性使得問題具有唯一解。

目前普遍采用l0范數(shù)和l1范數(shù)衡量目標(biāo)稀疏度。對(duì)于n維向量x=(x1,x2,…,xn)，l0范數(shù)為x中的非零元素個(gè)數(shù)，用‖x‖0表示。雖然‖x‖0不滿足范數(shù)的定義，但為了統(tǒng)一起見，仍稱其為范數(shù)。l1范數(shù)為x中所有元素絕對(duì)值之和，用‖x‖1表示。如果問題（1）中的噪聲項(xiàng)w滿足‖w‖2＜ε（ε為一個(gè)近似于零的正閾值），且信號(hào)x∈?n×1是稀疏的，即x的l0范數(shù)遠(yuǎn)小于信號(hào)長(zhǎng)度，則問題（1）可轉(zhuǎn)化為帶約束的非凸優(yōu)化模型（模型1）：

求解l0范數(shù)的過程是NP（Non-deterministic Polynomi?al-time）難問題。由于l1范數(shù)可以很好地逼近l0范數(shù)，因此可考慮用l1項(xiàng)替換模型1 中的l0項(xiàng)，從而得到一個(gè)凸的優(yōu)化模型（模型2）：

Donoho 等［8-9］證明了當(dāng)矩陣A滿足一定的零空間性質(zhì)時(shí)，模型1 與模型2 的解是等價(jià)的；Candès 等［10-11］證明了當(dāng)矩陣A滿足一定的約束等距性質(zhì)（Restricted Isometry Property，RIP）時(shí)，模型1 與模型2 的解是穩(wěn)定且等價(jià)的。此外，Candès 等［12］還證明了若矩陣A∈?m×n為隨機(jī)高斯矩陣，且A的行列滿足m=O()kln(n/k)，矩陣A滿足k階約束等距性質(zhì)的概率不低于。在此基礎(chǔ)上，Rul?son 等［13-14］證明了當(dāng)隨機(jī)高斯矩測(cè)量次數(shù)k滿足一定條件時(shí)，通過模型2 能夠精確地恢復(fù)原始信號(hào)。由于l1項(xiàng)是凸的，因此可利用凸優(yōu)化算法快速求解模型2。為便于求解，大部分算法先將模型2 轉(zhuǎn)化為無(wú)約束形式（模型3）：

其中，λ＞0 是一個(gè)正則化參數(shù)。模型2 到模型3 的轉(zhuǎn)換在數(shù)學(xué)上可表述為：對(duì)于任意的ε＞0，存在一個(gè)λ使得模型2 與模型3 的最優(yōu)解相同。

1.2 基于框架的圖像恢復(fù)原理

數(shù)字圖像的退化主要包括像素丟失、圖像變模糊及像素點(diǎn)含有噪聲等。一張灰度圖像在數(shù)學(xué)上通常用向量x=(x1,x2,…,xn) 表示，其中xi∈[0,255] 為對(duì)應(yīng)像素點(diǎn)的值。進(jìn)入21 世紀(jì)之后，基于框架（特別是緊小波框架）的圖像恢復(fù)方法得到了迅速發(fā)展。對(duì)一個(gè)給定的矩陣D∈?n×m(n≤m)，如果D滿足DDT=I，其中I是單位矩陣，DT是D的轉(zhuǎn)置，則D是一個(gè)緊框架。因此，對(duì)任何向量x∈?n，有：

此外，構(gòu)成向量u=DTx的元素是用來(lái)表示x的規(guī)范系數(shù)。如果緊框架D的元素是通過有限多個(gè)小波函數(shù)的擴(kuò)展與轉(zhuǎn)換生成的，則D被稱為一個(gè)緊小波框架。

利用壓縮感知理論的一個(gè)重要前提是測(cè)量目標(biāo)是稀疏的。自然圖像雖無(wú)法滿足稀疏條件，但自然圖像的小波系數(shù)通常是稀疏或者可壓的。小波分析采用細(xì)分函數(shù)構(gòu)造小波函數(shù)，并通過小波函數(shù)的伸縮與平移構(gòu)成函數(shù)空間的一組正交基或框架，空間中任一函數(shù)均可由小波函數(shù)生成的基或框架展開。在數(shù)字圖像處理過程中，大多使用高通和低通濾波器作為圖像分解與重構(gòu)的基本工具，小波系數(shù)可呈現(xiàn)出圖像在不同頻域中的特征?，F(xiàn)實(shí)中的大規(guī)模數(shù)據(jù)往往具有某種意義上的低維結(jié)構(gòu)或特征。對(duì)于一般的自然圖像，小波系數(shù)具有稀疏性特征，例如200 萬(wàn)像素的數(shù)字圖像僅需10 萬(wàn)個(gè)小波系數(shù)便能獲取幾乎所有可見信息，剩余的190 萬(wàn)小波系數(shù)僅貢獻(xiàn)少量、可忽略的細(xì)節(jié)信息。在數(shù)學(xué)上，小波系數(shù)的稀疏性可通過‖DTx‖p0≤p≤1 進(jìn)行刻畫。關(guān)于多尺度分析、細(xì)分函數(shù)、小波函數(shù)、小波濾波器等嚴(yán)格的數(shù)學(xué)理論和應(yīng)用可參考文獻(xiàn)［15］、［16］等。在稀疏性假設(shè)下，由文獻(xiàn)［17］可知，圖像恢復(fù)問題可建模為模型4：

同理，模型4 可轉(zhuǎn)化為無(wú)約束形式（模型5）：

模型5 也被稱為ALASSO 模型［18］。Lin 等［18-19］證明了如果矩陣A滿足一定的D-RIP 條件，則模型5 的最優(yōu)解是穩(wěn)定的。

一般情況下，自然圖像在空間域上具有一定連續(xù)性，因此自然圖像所產(chǎn)生的小波系數(shù)一般較大，而高斯噪聲經(jīng)過小波變換后的系數(shù)具有較強(qiáng)的隨機(jī)性，且噪聲對(duì)應(yīng)的大部分小波系數(shù)都很小。若設(shè)噪聲的小波系數(shù)方差為σ，則根據(jù)切比雪夫不等式，絕大部分噪聲系數(shù)都位于[-3σ,3σ]區(qū)間內(nèi)。因此，信號(hào)經(jīng)過小波變換后，使用閾值法可以較好地去除噪聲信號(hào)。

1.3 摩爾紋噪聲

在數(shù)碼影像中，如果主體中有密紋的紋理，通常會(huì)出現(xiàn)類似水波的條紋，該條紋即為摩爾紋。圖1 表示含有摩爾紋噪聲的小丑圖像去噪結(jié)果。

Fig.1 Denoising result of image“Clown”圖1 小丑圖像去噪結(jié)果

其中，圖1（a）為含有摩爾紋噪聲的小丑圖像，其傅里葉頻譜如圖1（b）所示。圖1（b）中正中心的高亮部分表示直流通分量，圖像中心附近的白色十字狀光譜特征即為摩爾紋噪聲在傅立葉頻譜中的尖峰狀分量。因此，通過校正這些由噪聲改變的頻譜分量，可以有效抑制噪聲的影響。圖1（c）表示使用軟閾值模型的處理結(jié)果，圖1（d）表示提取的摩爾紋噪聲。具體而言，對(duì)于一幅尺寸為M×N的光滑自然圖像f，其傅立葉頻譜為：

其中，H(k,l) 為傅立葉頻譜H的第(k,l)個(gè)系數(shù)。傅立葉頻譜中波峰位置可通過人工標(biāo)注或使用閾值法加以確認(rèn)。對(duì)于傳統(tǒng)閾值法，判定頻譜中第(k,l)點(diǎn)為波峰的條件為：該點(diǎn)系數(shù)Hkl與其周圍I×J個(gè)點(diǎn)系數(shù)的中值MEDI×J(Hkl)滿足不等式：

其中，CL×F為大小為L(zhǎng)×F的低頻置信區(qū)域，T2為給定閾值。在確認(rèn)波峰所在位置之后，選用不同類型濾波器校正由噪聲改變的頻譜分量過程可表示為：

圖2 表示高斯陷波濾波器的頻域處理過程。

Fig.2 Gaussian notch filtering illustration圖2 高斯陷波濾波器頻域處理過程

對(duì)于周期性噪聲，可選用高斯陷波濾波器［20］進(jìn)行濾波，即濾波器GI×J滿足：

2 基于緊小波框架的最小化模型

對(duì)于包含摩爾紋噪聲的圖像f，通過閾值法可確認(rèn)其由噪聲引起頻譜中的波峰位置，然后根據(jù)這些波峰位置可建立一個(gè)掩膜矩陣（Mask matrix）。在數(shù)學(xué)上，掩膜矩陣建立過程可表示為：

其中，Pi,j為掩膜矩陣P的第(i,j)個(gè)元素；s為給定閾值，用于規(guī)范矩陣大小。掩膜矩陣的建立確定了需要恢復(fù)的頻域系數(shù)范圍。由于在緊小波框架下，這部分系數(shù)是稀疏的，根據(jù)模型5，可得到優(yōu)化模型（模型6）：

其中，? 表示傅立葉卷積過程，P為利用公式（13）求得的掩膜矩陣，y=P?f，λ是正則化參數(shù)，D是一個(gè)緊小波框架，x為測(cè)量到的圖像，x為原始圖像。為便于表述，本文稱Θ=1 時(shí)的模型6 為軟閾值模型，Θ=0 時(shí)的模型6 為硬閾值模型。

求解模型6 的難點(diǎn)在于其中的正則項(xiàng)‖Dx‖Θ是不可分離的。DFISTA（Decomposition-basedFast Iterative Shrink?age Thresholding Algorithm，DFISTA）算法［21］提供了一個(gè)可行的分解思路：通過引入輔助變量，從而以簡(jiǎn)單和顯式的方式應(yīng)用迭代收縮閾值算法求解模型。具體而言，對(duì)于最小化問題，引入輔助變量τ，得到損失函數(shù)（lose function），如模型7 所示：

由于l1范數(shù)是不光滑的，在τ為零處不可微，因此需要使用近端梯度法（Proximal Gradient Method，PG）。PG 算法思想是使用近端算子作為近似梯度，從而進(jìn)行梯度下降。近端梯度下降法及其加速算法的收斂性分析可參見文獻(xiàn)［21］。

則函數(shù)?J的利普希茨常數(shù)滿足不等式：

其中，ε是一個(gè)近似于零的正閾值。利用DFISTA 算法求解模型7 的過程如下：

由以上過程可知，影響恢復(fù)結(jié)果的主要因素包括：①檢測(cè)波峰的閾值T2；②規(guī)范矩陣大小的閾值s；③緊小波框架函數(shù)選取；④小波閾值函數(shù)選取與閾值T1。

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

3.1 實(shí)驗(yàn)說明

為了驗(yàn)證本文算法的可行性與有效性，利用Matlab R2019a 軟件進(jìn)行仿真實(shí)驗(yàn)。實(shí)驗(yàn)對(duì)象選取Barbara、Boat、Cameramen 等4 幅灰度圖像以及2 幅使用智能手機(jī)拍攝的彩色自然圖像。對(duì)于純凈的灰度圖像，首先在對(duì)應(yīng)傅里葉頻域中靠近中心的低頻信號(hào)區(qū)域添加干擾噪聲，然后對(duì)使用加速算法與未使用加速算法求解軟閾值模型6 的PSNR（Peak Signal to Noise Ratio）值變化情況，以及分別使用中值濾波器［22］、窗口高斯陷波濾波器、由掩模矩陣直接生成的濾波器與模型7 恢復(fù)結(jié)果的PNSR 值進(jìn)行對(duì)比。由于在傅里葉頻域中，靠近中心的低頻信號(hào)區(qū)域決定了大部分圖像細(xì)節(jié)信息，因此添加的摩爾紋噪聲主要集中在低頻區(qū)域。對(duì)于彩色自然圖像，將其紅色、綠色與藍(lán)色通道圖像分別作為灰度圖進(jìn)行處理，最后再將去噪結(jié)果重新合成彩色圖像［23］。

圖3 表示對(duì)灰度圖像去除摩爾紋噪聲的實(shí)驗(yàn)流程。首先使用快速離散傅里葉逆變換函數(shù)得到與輸入圖像尺寸相同的摩爾紋噪聲圖像，然后向圖像添加摩爾紋噪聲并使用式（8）得到標(biāo)記的傅里葉頻域波峰位置，最后分別利用頻域?yàn)V波器與模型6 對(duì)含有摩爾紋噪聲的圖像進(jìn)行處理，并對(duì)比實(shí)驗(yàn)結(jié)果的峰值信噪比。其中，頻域?yàn)V波器包括高斯陷波濾波器、中值濾波器以及直接使用掩膜矩陣生成的濾波器。

Fig.3 Experimental procedures圖3 實(shí)驗(yàn)流程

為便于比較實(shí)驗(yàn)結(jié)果，設(shè)式（8）中，L=F=8，I=J=11，T2=5；式（10）中，α=1，β=0.1。根據(jù)波峰檢測(cè)結(jié)果，設(shè)s為大于1 且小于6 的整數(shù)。模型7 中，設(shè)當(dāng)Θ=1 時(shí)，μ=0.7；當(dāng)Θ=0 時(shí)，μ=1。使用分段樣條緊小波框架，設(shè)小波函數(shù)分解層次為4 層。

為能客觀地比較算法結(jié)果，采用以下定義的PSNR 值對(duì)圖像恢復(fù)質(zhì)量進(jìn)行評(píng)價(jià)：

其中，x表示恢復(fù)的圖像，表示原始圖像。

3.2 實(shí)驗(yàn)結(jié)果

為驗(yàn)證加速近端梯度算法的有效性與穩(wěn)定性，在使用含有噪聲的Cameramen 圖像作為實(shí)驗(yàn)對(duì)象時(shí)，對(duì)比使用加速算法與未使用加速算法的收斂速度。即在利用DFISTA算法求解模型7 的過程中，設(shè)為加速算法，tk+1=tk為未加速算法。使用加速算法求解軟閾值模型6 的迭代結(jié)果如圖4 所示。

Fig.4 Iterative results of noised image“Cameramen”圖4 含噪聲圖像Cameramen 迭代結(jié)果

由圖4 可以看出，當(dāng)?shù)螖?shù)達(dá)到50 次之后，人眼已難以觀察出實(shí)驗(yàn)結(jié)果與原始圖像的差別。

圖5 給出了隨著迭代次數(shù)增加，含噪聲圖像Camera?men 的PSNR 值變化情況。其中，紅色圈狀軌跡表示未使用加速算法的軟閾值模型6 的處理結(jié)果，藍(lán)色的“+”號(hào)軌跡表示使用加速算法的軟閾值模型6 的處理結(jié)果。從圖中可以看出，圖像Cameramen 在使用加速算法求解軟閾值模型6 時(shí)能更快地收斂到合適的結(jié)果。

圖6 給出了4 種去噪算法的對(duì)比結(jié)果。

Fig.5 Variation of the PSNR for the different algorithms on nosed image“Cameramen”圖5 含噪聲圖像Cameramen PSNR 值變化情況

Fig.6 Comparison of the result for fouralgorithms圖6 4 種去噪算法結(jié)果對(duì)比

由圖6 可以看出，采用不同算法處理Barbara 與Cam?eramen 圖像的結(jié)果各不相同。其中，使用模型6 處理圖像，結(jié)果的邊緣震蕩效果更弱，更符合人的視覺審美。由于高斯陷波濾波器以及掩膜矩陣濾波器對(duì)傅里葉頻域系數(shù)的處理只是進(jìn)行相應(yīng)縮放或截?cái)?，因此結(jié)果的頻域系數(shù)相較于原始頻域系數(shù)可能會(huì)有較大誤差，具體表現(xiàn)為實(shí)驗(yàn)結(jié)果中實(shí)物邊緣產(chǎn)生的震蕩效應(yīng)。實(shí)驗(yàn)圖片處理結(jié)果的PSNR 值如表1 所示。

表1 給出了所有實(shí)驗(yàn)對(duì)象處理結(jié)果的PSNR 值。由表1 可以看出，使用模型6 處理圖像的結(jié)果相比頻域?yàn)V波器有了較大提升，驗(yàn)證了利用稀疏性逼近的傅里葉頻域系數(shù)相較于高斯陷波濾波器的濾波結(jié)果更接近原始頻域系數(shù)。由于l1范數(shù)的凸性質(zhì)，在大部分情況下，使用軟閾值模型6的結(jié)果比使用硬閾值模型6 的結(jié)果更逼近全局最優(yōu)解。

圖7 表示2 張含有摩爾紋噪聲的彩色自然圖像及其處理結(jié)果，圖片來(lái)源于使用智能手機(jī)攝像頭分別拍攝IPS面板與TN 面板計(jì)算機(jī)液晶顯示器得到的彩色圖像。由圖7 可以看出，頻域?yàn)V波器無(wú)法在保護(hù)原始圖像信息的前提下去除摩爾紋噪聲。通過觀察軟閾值模型6 的處理結(jié)果可以發(fā)現(xiàn)，雖然圖像部分細(xì)節(jié)有所丟失，但對(duì)摩爾紋噪聲的總體去除效果較為顯著。

Table 1 Comparison of PSNR of image restoration results表1 圖像恢復(fù)結(jié)果PSNR 值對(duì)比

Fig.7 Denoising results of color images圖7 彩色圖像去噪結(jié)果

4 結(jié)語(yǔ)

本文通過分析摩爾紋噪聲的傅里葉頻域特性，針對(duì)傳統(tǒng)頻域?yàn)V波器無(wú)法準(zhǔn)確恢復(fù)低頻信號(hào)的缺陷，提出一種基于稀疏表示的摩爾紋噪聲去除方法。首先通過閾值法確定頻域中受噪聲影響的區(qū)域，然后根據(jù)這些區(qū)域位置建立掩膜矩陣。由于緊小波框架下自然圖像的小波系數(shù)具有稀疏性，因此根據(jù)壓縮感知原理可建立優(yōu)化模型稀疏逼近原始圖像的傅里葉頻域系數(shù)。通過使用DFISTA 算法的優(yōu)化模型，求解過程能夠收斂到穩(wěn)定的結(jié)果。數(shù)值實(shí)驗(yàn)結(jié)果表明，當(dāng)噪聲影響區(qū)域涉及低頻區(qū)域時(shí)，本文提出的圖像恢復(fù)算法相比中值濾波方法與高斯陷波濾波方法不僅具有更高的峰值信噪比，而且能保留更清晰的細(xì)節(jié)信息。以含有摩爾紋噪聲的彩色自然圖像作為實(shí)驗(yàn)對(duì)象的圖像恢復(fù)數(shù)值實(shí)驗(yàn)也證明了該算法的有效性與可行性。