基于稀疏表示的摩爾紋噪聲去除方法

（浙江理工大學 理學院，浙江 杭州 310018）

0 引言

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

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

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

1 圖像恢復原理

1.1 稀疏信號恢復

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

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

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

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

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

其中，λ＞0 是一個正則化參數。模型2 到模型3 的轉換在數學上可表述為：對于任意的ε＞0，存在一個λ使得模型2 與模型3 的最優解相同。

1.2 基于框架的圖像恢復原理

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

此外，構成向量u=DTx的元素是用來表示x的規范系數。如果緊框架D的元素是通過有限多個小波函數的擴展與轉換生成的，則D被稱為一個緊小波框架。

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

同理，模型4 可轉化為無約束形式（模型5）：

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

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

1.3 摩爾紋噪聲

在數碼影像中，如果主體中有密紋的紋理，通常會出現類似水波的條紋，該條紋即為摩爾紋。圖1 表示含有摩爾紋噪聲的小丑圖像去噪結果。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實驗結果與分析

3.1 實驗說明

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

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

Fig.3 Experimental procedures圖3 實驗流程

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

為能客觀地比較算法結果，采用以下定義的PSNR 值對圖像恢復質量進行評價：

其中，x表示恢復的圖像，表示原始圖像。

3.2 實驗結果

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

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

由圖4 可以看出，當迭代次數達到50 次之后，人眼已難以觀察出實驗結果與原始圖像的差別。

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

圖6 給出了4 種去噪算法的對比結果。

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 種去噪算法結果對比

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

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

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

Table 1 Comparison of PSNR of image restoration results表1 圖像恢復結果PSNR 值對比

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

4 結語

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