基于稀疏表示和低秩矩陣逼近的圖像去噪算法的研究*

盧文鋒，佀同光，韓國(guó)勇

（山東建筑大學(xué)管理工程學(xué)院，山東 濟(jì)南 250101）

0 引 言

光學(xué)相干斷層成像（optical coherence tomography，OCT）是一種通過(guò)測(cè)量物體后向反射光強(qiáng)度進(jìn)行成像的技術(shù)，能夠?qū)崟r(shí)快速地對(duì)體內(nèi)組織的表層結(jié)構(gòu)進(jìn)行成像.OCT系統(tǒng)的基本組成光路為邁克爾遜光路［1］，根據(jù)具體使用的光路器件的不同，主要分為時(shí)域OCT（time domain optical coherence tomography，TDOCT）和頻域 OCT（spectral optical coherence tomography，SOCT）.TDOCT 中，寬帶光源發(fā)出的光通過(guò)分束器得到2束相干光，一束照射需要成像的樣品，另一束照射參考鏡片；樣品反射的樣品光與參考鏡片反射的反射光在分束器發(fā)生干涉；光路中的光電探測(cè)器將光信號(hào)轉(zhuǎn)換為電信號(hào)，并發(fā)送到計(jì)算機(jī)進(jìn)行處理.通過(guò)樣品端機(jī)械臂的橫向移動(dòng)，獲得不同位置的信息，參考端機(jī)械臂的移動(dòng)獲得某一位置的深度信息.與TDOCT系統(tǒng)結(jié)構(gòu)類似，SOCT是在TDOCT光路的基礎(chǔ)上用光譜儀替換光電探測(cè)器.光譜儀能夠獲得低相干干涉的光譜信號(hào)，通過(guò)傅里葉變換，能夠一次獲取某一位置縱向的空間結(jié)構(gòu)信息，故SOCT系統(tǒng)只需樣品端機(jī)械臂的橫向掃描即可.

OCT系統(tǒng)中的樣品光，可看作干涉區(qū)內(nèi)被測(cè)樣品所有背向散射的小波之和，樣品中散射點(diǎn)深度分布的隨機(jī)性及樣品中折射率的起伏使這些小波的相位具有隨機(jī)性［2］，因此，OCT圖像不可避免地帶有散斑噪聲.而散斑噪聲與成像光束的波長(zhǎng)和成像物體的細(xì)節(jié)結(jié)構(gòu)有關(guān)［3］，會(huì)降低信噪比、系統(tǒng)分辨率和探測(cè)靈敏度［4］.

圖像的自相似性是指圖像內(nèi)部含有大量的非局部相似塊，相似塊矢量化后的矩陣是低秩的.近年來(lái)基于圖像的自相似性去噪算法發(fā)展迅速.非局部均值（non-local means，NLM）算法首次應(yīng)用圖像自相似先驗(yàn)進(jìn)行圖像去噪［5-6］.在NLM算法中，通常通過(guò)某個(gè)鄰域中所有像素的加權(quán)平均值進(jìn)行濾波.然而，與每個(gè)給定像素相關(guān)聯(lián)的權(quán)重并不取決于其與目標(biāo)像素的幾何距離，而是取決于其與目標(biāo)像素的相似性.本文基于圖像的非局部自相似性，首先尋找圖像中的相似塊［7］，然后根據(jù)稀疏表示方式和低秩矩陣近似進(jìn)行去噪.

1 基于稀疏表示和矩陣逼近的去噪算法

本文提出了基于稀疏表示和矩陣逼近的去噪算法，通過(guò)對(duì)圖像子塊進(jìn)行線性表示，代替對(duì)整幅圖像進(jìn)行線性表示來(lái)減少計(jì)算量，同時(shí)利用奇異值分解（singular value decomposition，SVD）在數(shù)據(jù)表示方面的巨大優(yōu)勢(shì)［8］.（1）通過(guò)對(duì)圖像塊的數(shù)據(jù)矩陣進(jìn)行奇異值分解，進(jìn)行全局子空間分析，確定信號(hào)子空間和噪聲子空間；（2）計(jì)算圖像塊與信號(hào)子空間的距離，將相似的塊分組為訓(xùn)練樣本；（3）對(duì)相似塊矩陣進(jìn)行奇異值分解，并確定表示相似塊的奇異向量；（4）去除屬于噪聲子空間的候選樣本，即去除表示噪聲的基，得到每個(gè)塊的局部基.

1.1 全局子空間分析

設(shè)一幅含有噪聲的灰度圖像中共含有L個(gè)圖像塊，每個(gè)圖像塊的尺寸為其中n為圖像中含有的像素個(gè)數(shù).定義數(shù)據(jù)矩陣即Y中的每一個(gè)列向量為一個(gè)圖像塊的數(shù)據(jù).奇異值分解表示為

鑒于SVD在表示矩陣時(shí)具有的優(yōu)良特性，通過(guò)奇異值分解來(lái)解決低秩矩陣近似問(wèn)題時(shí)，首先計(jì)算矩陣Y的秩r=rank(Y)，僅通過(guò)采用前r個(gè)奇異值所對(duì)應(yīng)的奇異向量，重建矩陣Y［9］.即對(duì)于問(wèn)題

矩陣N經(jīng)過(guò)式（1）的計(jì)算后，該問(wèn)題的解為

為了確定由表示信號(hào)的基構(gòu)成的信號(hào)子空間，需要解決式（4）的最優(yōu)化問(wèn)題，

式中：S為基的數(shù)量，Ω為噪聲標(biāo)準(zhǔn)差，lk為矩陣的第k個(gè)奇異值.即使噪聲塊與去噪塊之間的平均平方殘差接近噪聲方差，來(lái)確定信號(hào)子空間中包含的S，從而確定前S個(gè)奇異值所對(duì)應(yīng)的奇異向量來(lái)構(gòu)成信號(hào)子空間.

在確定由表示噪聲的基構(gòu)成的噪聲子空間時(shí)，運(yùn)用補(bǔ)集的思想先來(lái)確定表示非噪聲信號(hào)的基，即

式中：參數(shù)β用來(lái)控制估計(jì)的誤差，其值越大，表明噪聲子空間包含向量的個(gè)數(shù)越多.確定式（5）的解即可得到表示非噪聲信號(hào)基的數(shù)量C，從而確定奇異值所 對(duì) 應(yīng) 的 奇 異 向 量來(lái)構(gòu)成噪聲子空間 .

式（4和5）中使用的噪聲標(biāo)準(zhǔn)差Ω參考文獻(xiàn)［10］中的魯棒中值估計(jì)器（robust median estimator，RME）進(jìn)行確定，即

式中：Yij∈子帶HH1，HH1為小波變換中的細(xì)節(jié)子圖，代表了輸入圖像Y水平和垂直方向的高頻成分.

1.2 相似塊分組與表示

利用圖像的自相似性，將圖像中包含的相似塊進(jìn)行聚類，找到每個(gè)塊的訓(xùn)練樣本.聚類的方法和實(shí)現(xiàn)算法多種多樣，選取最簡(jiǎn)單的塊匹配（block matching，BM）方法［11］.相比于在全局中尋找相似塊的方式［12］，本文只在信號(hào)子空間中尋找相似塊.圖像塊xi與xj的距離定義為

式中：S為信號(hào)子空間中包含的向量個(gè)數(shù)，已由式（4）確定為圖像塊xi在信號(hào)子空間中的投影向量，定義為

圖像塊的距離越小，表示二者相似程度越高，然后選取閾值Th來(lái)控制樣本數(shù)據(jù)矩陣中含有相似塊的數(shù)量.

相對(duì)于采取奇異值收縮重建矩陣的方法，本文將樣本矩陣看作是無(wú)噪聲數(shù)據(jù)，由零均值加性高斯白噪聲污染后形成的，每一個(gè)樣本是獨(dú)立重復(fù)實(shí)驗(yàn)的觀測(cè)結(jié)果.通過(guò)假設(shè)檢驗(yàn)序列檢測(cè)源個(gè)數(shù)的算法估計(jì)噪聲標(biāo)準(zhǔn)差（σ）［13］，并基于隨機(jī)矩陣?yán)碚撏ㄟ^(guò)一系列假設(shè)檢驗(yàn)估計(jì)秩（r），并取前r個(gè)奇異向量得到局部基的候選向量.

1.3 去除表示噪聲的基

在進(jìn)行秩估計(jì)時(shí)，采取的是基于樣本獨(dú)立同分布的假設(shè)，在實(shí)驗(yàn)中并不成立，使得通過(guò)塊分組得到的樣本可能會(huì)有很大的差異，因此，應(yīng)排除噪聲占主導(dǎo)地位的樣本和噪聲子空間中的候選樣本.提出2種方式來(lái)去除表示噪聲的候選樣本.

考慮到圖像中不包含邊緣結(jié)構(gòu)的圖像塊信噪比往往較低，同時(shí)很難找到合適的樣本來(lái)訓(xùn)練局部基，因此，第1種方式考慮重建這些圖像塊.在完成第1步的全局子空間分析后，對(duì)于任一圖像塊xi，所在組為計(jì)算圖像塊組的能量為

若式（9）成立，即表明該圖像塊所在組的平均能量較低，直接跳過(guò)第2步，用xi的平均值作為對(duì)xi的估計(jì)值；同時(shí)通過(guò)式（10）設(shè)定方差的值為

式中：參數(shù)c＞1且為常數(shù).由于圖像中含有較多的相似塊，如果一個(gè)圖像塊是平滑的，其鄰近的相似塊也可能是平滑的.所以當(dāng)式（9）成立時(shí)，通過(guò)參數(shù)c設(shè)置一個(gè)較大的σ值，從而檢驗(yàn)下一個(gè)圖像塊來(lái)防止過(guò)擬合現(xiàn)象的發(fā)生.如果式（9）不成立，則繼續(xù)執(zhí)行第2步即可；此時(shí)σ的值仍在第2步中進(jìn)行確定.

由于2個(gè)信號(hào)的相似程度可以用向量?jī)?nèi)積來(lái)表示［14］，第2種方式便是通過(guò)計(jì)算向量的內(nèi)積來(lái)去除表示噪聲的基.對(duì)于任一候選向量bi，計(jì)算其與每一個(gè)奇異向量的內(nèi)積，找到使二者內(nèi)積最大的奇異向量αi；如果αi屬于噪聲子空間，則去除該候選向量.

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

選取真實(shí)的OCT圖像進(jìn)行實(shí)驗(yàn)，設(shè)置β=0.000 5，c=1.02，圖像塊的尺寸11×11，搜索窗的尺寸39×39.選取三維塊匹配（block matching 3D，BM3D）算法［15］、NLM 算法［5］、奇異值分解-低秩矩陣逼近（singular value decomposition-low rank approximations，SVD-LRA）算法［16］等作為比較對(duì)象，算法中的參數(shù)均為默認(rèn)值.

通常選取常見(jiàn)的峰值信噪比（peak signal to noise ratio，PSNR）與結(jié)構(gòu)相似度（structural similarity index measure，SSIM）作為評(píng)估去噪效果的客觀評(píng)價(jià)準(zhǔn)則.PSNR通過(guò)計(jì)算對(duì)應(yīng)像素點(diǎn)間的誤差來(lái)評(píng)價(jià)圖像品質(zhì)，去噪效果好的算法往往對(duì)應(yīng)較大的PSNR值，記為RPSN，其計(jì)算公式為

式中：Ems為均方誤差，Gmax為圖像灰度的最大值.SSIM從結(jié)構(gòu)、亮度和對(duì)比度3個(gè)方面來(lái)衡量2幅圖像相似程度.SSIM如下定義，記為ISSM

采用不同方法進(jìn)行降噪處理，并計(jì)算降噪圖像效果，對(duì)應(yīng)圖像和結(jié)果如圖1和表1所示.可知，本文方法對(duì)圖像OCT_1處理后PSNR為71.151 2 dB，對(duì)圖像OCT_2和圖像OCT_3處理后PSNR分別為74.464 9和80.587 7 dB，說(shuō)明本文方法對(duì)圖像OCT_2和圖像OCT_3處理后降噪效果比較好，高于其他3種方法，與目前流行的BM3D方法相比具有很好的競(jìng)爭(zhēng)力.尤其是對(duì)圖像OCT_3處理后的SSIM為0.931 4，最接近于1，說(shuō)明本文方法去噪后能夠與原始圖像保持良好的結(jié)構(gòu)相似性.對(duì)于圖像OCT_1，因?yàn)樯咴肼晣?yán)重地影響了OCT成像品質(zhì)，BM3D算法表現(xiàn)較好，較大程度上減少了散斑噪聲，本文提出的算法和SVD-LRA算法均模糊了OCT圖像，NLM算法仍保留了較多的散斑噪聲；對(duì)于圖像OCT_3，SVD-LRA算法也由于對(duì)圖像過(guò)度平滑，使圖像變得模糊而丟失細(xì)節(jié)部分，本文提出的算法有效地去除了散斑噪聲而且能夠保留圖像中的細(xì)節(jié)信息，從視覺(jué)上優(yōu)于其他3種方法.

圖1 3幅OCT圖像經(jīng)不同方法處理圖像對(duì)比

3 結(jié)束語(yǔ)

相對(duì)于目前比較先進(jìn)的降噪算法，本文提出的基于稀疏表示和矩陣逼近的圖像去噪算法，具有以下優(yōu)勢(shì)：通過(guò)對(duì)圖像子塊進(jìn)行線性表示代替對(duì)整幅圖像進(jìn)行線性表示來(lái)減少計(jì)算量；利用圖像塊與信號(hào)子空間的距離來(lái)尋找相似塊，避免在全局中匹配相似塊；充分利用了奇異值分解在數(shù)據(jù)表示方面的巨大優(yōu)勢(shì)，用來(lái)表示或恢復(fù)矩陣.從對(duì)OCT圖像降噪實(shí)驗(yàn)的結(jié)果來(lái)看，相似塊的選取是關(guān)鍵.當(dāng)圖像中結(jié)構(gòu)相似子塊較多時(shí)，本文提出的方法處理較好.如何準(zhǔn)確而有效地選取圖像相似子塊，提高結(jié)構(gòu)相似性匹配率是下一步要研究的問(wèn)題.