李佳慶，雷 蕾

（1.山西經(jīng)濟(jì)管理干部學(xué)院，山西 太原 030024；2.山西城市職業(yè)技術(shù)學(xué)院，山西 太原 030027）

1 綜述

X-CT 低劑量薄層掃描，即“薄掃”，能更精準(zhǔn)地檢出肺癌早期關(guān)鍵病灶。通過X 射線斷層掃描和低劑量計(jì)算機(jī)斷層掃描（Computed Tomography，CT）技術(shù)，對(duì)病灶進(jìn)行篩查，發(fā)現(xiàn)病灶后再進(jìn)行局部高分辨率掃描的檢查。分析時(shí)，把電子探測(cè)器接收信號(hào)轉(zhuǎn)化為數(shù)字信號(hào)輸入計(jì)算機(jī)，再由計(jì)算機(jī)轉(zhuǎn)化為圖像。這種掃描已成為肺癌早期診斷最有效、最直接的影像學(xué)方式。

本文研究低輻射劑量、薄層掃描技術(shù)的CT 肺部序列結(jié)構(gòu)影像，對(duì)連續(xù)高密掃描影像進(jìn)行檢測(cè)分析。由于成像設(shè)備和獲取條件等因素的影響，需要對(duì)部分出現(xiàn)質(zhì)量問題的序列圖像進(jìn)行快速去噪，利用圖像去噪技術(shù)改善CT 圖像的質(zhì)量，在最大可能地保留和增強(qiáng)醫(yī)學(xué)細(xì)微特征和信息的基礎(chǔ)上，解決不完備投射帶來(lái)的影像數(shù)據(jù)質(zhì)量問題；去除低劑量方式帶來(lái)的圖像噪聲，在去噪過程中避免將微小結(jié)節(jié)誤清除。

由于醫(yī)學(xué)圖像去噪的特殊性，本文擬在研究圖像去噪算法的基礎(chǔ)上，研究保留醫(yī)學(xué)細(xì)微特征信息，盡可能地去除掉圖像的干擾成分，同時(shí)保證原始圖像有用信息的完整度。本文研究算法實(shí)現(xiàn)快速有效的圖像去噪處理方法來(lái)支持肺部病灶的檢出，實(shí)現(xiàn)CT 薄層掃描下序列醫(yī)學(xué)圖像快速去噪。

目前，有許多學(xué)者在圖像去噪方面做出了研究。MUKHOPADHYAY 等人[1]提出了一種利用遺傳算法的隨機(jī)化技術(shù)對(duì)閾值進(jìn)行設(shè)定與優(yōu)化的圖像去噪方法，利用該方法將含噪聲的圖像分成固定大小的小塊，通過小波變換對(duì)每個(gè)小塊圖像進(jìn)行去噪。YANG H 等人[2]針對(duì)圖像去噪時(shí)對(duì)邊緣和紋理保護(hù)的問題，提出了一種基于雙支持向量機(jī)的圖像去噪方法，在對(duì)圖像進(jìn)行去噪的過程中充分利用Shearlet 變換的多分辨率分析和最優(yōu)逼近的特點(diǎn)，對(duì)圖像進(jìn)行去噪的同時(shí)很好地保留了圖像的邊緣特征。PROTTER 等人[3]提出了一個(gè)通過稀疏和冗余表示對(duì)序列圖像去噪的方法，比循環(huán)地處理單張圖像有更好的峰值信噪比（Peak Signal to Noise Ratio，PSNR）。LI H 等人[4]根據(jù)小波變換、稀疏表示和冗余表示，提出了一個(gè)單尺度小波K-SVD 算法，對(duì)于強(qiáng)噪聲圖像處理有很好的峰值信噪比。

本文考慮薄層掃描切片及切片上肺結(jié)節(jié)的分析，在分析肺部CT 序列上下層切片相關(guān)性的基礎(chǔ)上，結(jié)合專家提供的經(jīng)驗(yàn)知識(shí)，設(shè)計(jì)非局部信息稀疏和冗余表示的圖像去噪方法，將圖像序列轉(zhuǎn)化成一個(gè)既包含時(shí)間相關(guān)度又包含空間相關(guān)度的體素，從而將圖像的去噪問題轉(zhuǎn)換成一個(gè)基于冗余字典的最優(yōu)稀疏表示的二次尋優(yōu)問題。通過序列圖像信息與稀疏表示方法對(duì)圖像進(jìn)行去噪，從而將圖像的去噪問題轉(zhuǎn)換成基于冗余字典的最優(yōu)稀疏表示的尋優(yōu)問題[5]。

2 相關(guān)工作

2.1 信號(hào)處理

時(shí)域和頻域是信號(hào)的基本性質(zhì)，為信號(hào)分析提供了不同的角度。對(duì)于需要處理的信號(hào)來(lái)說，在時(shí)域上本身就具有稀疏性的信號(hào)很少。為了尋找更好的稀疏，總能找到某種變換，使得在某種變換域中稀疏，即能夠用較少的系數(shù)來(lái)表達(dá)。

常見的變換有正交變換、余弦基變換（Discrete Cosine Transformation，DCT）、小 波 變 換（Wavelet Transformation，WT）、多分辨率分析（小波變換）、多 尺 度 幾 何 分 析[6]（包 括Ridgelet、Curvelet、Bandelet、Contourlet 等一系列多尺度分析工具）。這些工具符合人類視覺皮層對(duì)圖像有效表示的要求，通過變換的數(shù)學(xué)分析方法把圖像分解在不同的尺度（不同的精細(xì)程度）上來(lái)處理。

2.2 稀疏表示去噪原理

含噪的圖像是有一定的結(jié)構(gòu)的，可以從信號(hào)中提取出能表達(dá)圖像性質(zhì)的原子，即將信號(hào)表示為字典和稀疏的乘積，將信號(hào)提取出來(lái)，即為去噪后的圖像。而噪聲是隨機(jī)、不相關(guān)的，沒有結(jié)構(gòu)特性，因此提取不出來(lái)。通過這樣的原理可以將圖像提取出來(lái)，以達(dá)到圖像去噪的目的。

2.3 稀疏表示

稀疏表示[7]是獲得對(duì)信號(hào)的良好近似，優(yōu)化模型是從信號(hào)重建的角度來(lái)設(shè)計(jì)的，將信號(hào)表示為字典的線性表示。即，y=Ax，其中y是待處理的信號(hào)，A是字典，x是稀疏系數(shù)。稀疏（性）是指x的大部分元素為0 或接近0，只有少部分元素不為0。其中，x滿足稀疏性。稀疏表示的目的是選擇最少的系數(shù)，從而重構(gòu)信號(hào)y。已知量是輸入信號(hào)y，未知量是字典A與稀疏系數(shù)x。對(duì)信號(hào)稀疏表示的目的是尋找一個(gè)字典使得信號(hào)在該字典下的系數(shù)（表達(dá)）最稀疏。

稀疏表示中稀疏的原因是：字典是在一個(gè)足夠大的訓(xùn)練樣本空間內(nèi)，對(duì)于一個(gè)類別的物體，可以大致由訓(xùn)練樣本中的樣本子空間線性表示，因此當(dāng)該物體有整個(gè)樣本空間表示時(shí)，其表示的系數(shù)是稀疏的[8]。

稀疏表示的兩個(gè)主要的問題是字典的生成和稀疏分解（編碼）。

2.3.1 字典的生成

字典是一個(gè)矩陣，一個(gè)原子（基）是字典A矩陣中的一列。由這樣一列一列的原子排列成的矩陣就是一個(gè)字典。字典表達(dá)信號(hào)的能力取決于信號(hào)的特征是否與字典中基的特征相吻合。字典是基的組合，有表達(dá)更多、更復(fù)雜的信號(hào)的能力。

字典有固定字典和自適應(yīng)字典兩種。常用的固定字典有超完備字典DCT、小波字典Contourlet 和曲波字典Wavelet 等。固定字典不能隨著信號(hào)的變化做出相應(yīng)的變化，一經(jīng)選定，對(duì)所有的信號(hào)一視同仁，信號(hào)的表達(dá)形式單一；而學(xué)習(xí)字典是通過大量圖像數(shù)據(jù)學(xué)習(xí)得到的，通過訓(xùn)練學(xué)習(xí)得到的字典可以根據(jù)輸入信號(hào)的不同做出自適應(yīng)的變化，能夠更好地適應(yīng)不同的圖像數(shù)據(jù)。

2.3.2 稀疏分解（編碼）

求信號(hào)在字典上的稀疏表示，即求出信號(hào)y在字典A上的稀疏表示x，常用的方法是正交匹配追蹤 算 法（Orthogonal Matching Pursuit，OMP）。主要目標(biāo)是找出x中最主要的k個(gè)分量和值。假設(shè)x=(x0,0,0,x1,x2,…)。為尋求與y最為接近的x，首先需計(jì)算信號(hào)與原子（字典中的每一列）之間的內(nèi)積。內(nèi)積可理解為計(jì)算矢量在某一指定方向上的投影長(zhǎng)度。其次，確定最大值對(duì)應(yīng)的原子，并利用最小二乘法來(lái)確立其稀疏系數(shù)。在接下來(lái)的步驟中，將余下的原子與字典中相應(yīng)的列進(jìn)行內(nèi)積計(jì)算，以識(shí)別出與最大值匹配的字典列，不斷迭代，直至達(dá)到算法的收斂狀態(tài)。

稀疏分解算法的概念最早由MALLAT 提出，名為匹配追蹤算法（Matching Pursuit，MP），以其簡(jiǎn)潔和易于實(shí)現(xiàn)的特性得到了廣泛的應(yīng)用[9]。此后，PATI 等學(xué)者在MP 算法的基礎(chǔ)上提煉出OMP，其特點(diǎn)是具備更大的收斂速度。隨后的研究進(jìn)一步優(yōu)化了OMP 算法，孕育出如壓縮采樣匹配追蹤算法（Compressive Sampling Matching Pursuit，CoSaMP）、正則化正交匹配追蹤算法（Regularized Orthogonal Matching Pursuit，ROMP）、分段式正交匹配追蹤算法（Stagewise OMP，StOMP）以及子空間追蹤算法（Subspace Pursuit，SP）等多種變種。

2.3.3 優(yōu)化問題

2.3.3.1 凸優(yōu)化問題

字典A是一個(gè)N×M的矩陣。由于M遠(yuǎn)小于N，用M個(gè)方程求解N個(gè)未知數(shù)，即用少的方程求解多的未知數(shù)，因此y=A*x是個(gè)欠定方程（有無(wú)窮多解），說明使用字典的線性組合表達(dá)信號(hào)將是不唯一的，所以求解此問題需要找到最優(yōu)解。

解優(yōu)化問題，如果加上適當(dāng)?shù)南薅l件即有唯一解。那么需要解決的問題就是尋找是否存在一種最好的表達(dá)方式，使得在字典的表示下系數(shù)最稀疏。||x||0為x的稀疏度，表示x中非零系數(shù)的個(gè)數(shù)。要求解x盡可能地稀疏，即||x||0盡可能地小。0 范數(shù)的優(yōu)化問題很難求解，而1 范數(shù)優(yōu)化問題容易求解（存在且唯一）。在滿足一定條件下將問題轉(zhuǎn)換為1范數(shù)的優(yōu)化問題（最小化問題），此時(shí)優(yōu)化問題變成一個(gè)凸優(yōu)化問題[10]，應(yīng)用已有的理論就可以解決。

2.3.3.2 最大期望算法

優(yōu)化的模型是兩個(gè)參數(shù)相乘，需要同時(shí)估計(jì)字典A和稀疏表示系數(shù)x，基本思想是：第一步，將字典A固定，求出x的值，即求信號(hào)y在A上稀疏表示；第二步，使用上一步得到的x來(lái)更新字典A，即字典更新，使得y和Ax最接近。如此反復(fù)迭代幾次，即可得到優(yōu)化的A和x。

在稀疏表示里，已知y和A求x的方法是OMP。已知y，求一個(gè)比較好的A的方法是Module（MOD）和K-SVD 等。字典的更新方法主要有最 大 似 然（Maximum Likelihood，ML）、最 大 后驗(yàn) 概 率（Maximum A posteriori Probability,MAP）、MOD、FOCUSS 字典學(xué)習(xí)算法、主成分分析（Principal Component Analysis，PCA）、K-SVD 以 及Online 等。本文以應(yīng)用最廣泛的K-SVD[11]算法為例做對(duì)比實(shí)驗(yàn)。

輸入含噪序列圖像后，K-SVD[12]訓(xùn)練字典用于圖像去噪的算法步驟如下。

（1）選擇字典作為過完備字典，訓(xùn)練字典；選取字典集，并更新字典。

（2）將含噪圖像矩陣以m×m為一塊，將一個(gè)大的矩陣劃分為多個(gè)子矩陣。對(duì)于每個(gè)子矩陣，將元素按列連接為一個(gè)列向量，一個(gè)m×m的塊轉(zhuǎn)換為m2×1 的列向量，生成一個(gè)新的矩陣B。最后生成的矩陣B就是由A個(gè)塊矩陣生成的A列矩陣。矩陣B就作為字典學(xué)習(xí)的輸入。

（3）先取B中的M列，更新其中的每一個(gè)元素，B[i][j]B[i][j]-means{B[:,j]}B[i][j]=B[i][j]-means{B[:,j]}，式中B[i][j]表示B[i][j]矩陣中第i行第j列元素，means{B[:,j]}表示元素(i,j)所在列（第j列）的平均值。

（4）固定字典，利用OMP 算法，求解該B在字典下面的稀疏系數(shù)cofes。

（5）更新B的每一列：B[:][j]=B[:][j]+means{B[:,j]}，在此上下文中，B[:][j]第j列可被視為原子的代表。K-SVD 策略涉及逐一地更新原子（即字典的某一特定列）及其相應(yīng)的稀疏系數(shù)。只有當(dāng)所有的原子都經(jīng)歷了更新過程后，該算法才進(jìn)入下一迭代。經(jīng)過多次的精細(xì)迭代，最終能夠獲得經(jīng)過優(yōu)化的字典及其稀疏系數(shù)。

（6）取下一個(gè)M列，重復(fù)步驟（3）～（5），直到B矩陣的所有列都處理完成。

（7）將B的每一列(m2×1)轉(zhuǎn)換為m×m的小塊，逆運(yùn)算得到去噪后的圖像。

最終，輸出去噪后的序列圖像。

3 方法描述

本文采用模擬退火算法[13]得到最優(yōu)值，算法是一種隨機(jī)性組合優(yōu)化方法，來(lái)源于對(duì)固體退火降溫的過程。即將固體加溫至充分高，再讓其溫度慢慢退卻。固體加熱時(shí)，固體內(nèi)粒子隨溫度的升高而變?yōu)闊o(wú)序狀，而冷卻時(shí)，粒子逐漸趨于有序。這一過程中，每個(gè)溫度下都達(dá)到平衡狀態(tài)，不但能夠解優(yōu)化問題，而且能夠優(yōu)化字典中原子的組合，使字典具備更好的表達(dá)信息的能力，從而更好地去噪。

算法是從整體上考慮的，概率與系統(tǒng)能量E、溫度T密切相關(guān)。在高溫下，接受與當(dāng)前狀態(tài)能量差較大的新狀態(tài)；在低溫下，接受與當(dāng)前狀態(tài)能量差較小的新狀態(tài)。這樣避免陷入局部最優(yōu)，有的放矢地高效找出最優(yōu)解。所提方法描述如圖1 所示。模擬退火算法[14]稀疏表示去噪算法具體步驟如下。

圖1 方法描述

首先輸入y=A*x中A與x的無(wú)窮多個(gè)解，限定條件||x||0和||x||1盡可能?。ㄏ∈瑁?。

（1）初始化狀態(tài)信息，隨機(jī)產(chǎn)生初始解x0（即在溫度T條件下的解為x0）。

（2）設(shè)目標(biāo)（能量）函數(shù)E=y-y*，其中y*=Ax，計(jì)算差值ΔE=Ej-Ei，i為當(dāng)前狀態(tài)，j為下一狀態(tài)。若ΔE<0，則接受j為當(dāng)前狀態(tài)，即更接近最優(yōu)解（內(nèi)能變小更趨于有序狀態(tài)），否則以一定的概率來(lái)接 受 狀 態(tài)。min{1,exp(-ΔE/KT)}>random[0,1]，K是Metropolis 接受準(zhǔn)則，是以概率接受新狀態(tài)。溫度T為控制參數(shù)。對(duì)于exp(A)函數(shù)，當(dāng)A<0，其值小于1；當(dāng)A>0，其值大于1。ΔE<0，-ΔE>0，-Δf/TK>0，exp(Δf/TK)>1，min{1,exp(-Δf/TK)}取1，接受全部新值。當(dāng)解為差解時(shí)ΔE>0，-ΔE<0，-ΔE/KT<0，exp(-ΔE/KT)<1，min{1,exp(-ΔE/KT)}中取exp(-ΔE/KT)，最后通過概率p=exp(-ΔE/KT)>random[0,1]來(lái)接受新值，即概率p=exp[-(Ej-Ei)/KT]若大于[0,1]區(qū)間的隨機(jī)數(shù)，接受狀態(tài)j為當(dāng)前狀態(tài)；若不成立，則保留狀態(tài)i為當(dāng)前狀態(tài)。

（3）不斷產(chǎn)生新解，計(jì)算目標(biāo)函數(shù)差，判斷是否接受新解。經(jīng)過大量的解變化后，可以求得給定控制參數(shù)的優(yōu)化問題的相對(duì)最優(yōu)解。

（4）設(shè)定閾值n是最優(yōu)的n個(gè)解。保留這些解，尋找與圖像最接近的若干解，并對(duì)這些解進(jìn)行標(biāo)記。利用模擬退火算法來(lái)優(yōu)化自適應(yīng)字典，以找到每個(gè)解的最佳自適應(yīng)字典。最終，在所有解中選擇最優(yōu)解。

輸出最優(yōu)解，得到最優(yōu)的字典A和稀疏系數(shù)x，逆運(yùn)算Ax=y*得到去噪后的圖像，從而完成圖像去噪。

與廣泛使用的K-SVD 算法相比，K-SVD 算法需要先更新每一個(gè)元素再更新每一列，更新速度慢，迭代次數(shù)多。而本文算法是從宏觀的角度動(dòng)態(tài)自適應(yīng)地更新協(xié)調(diào)全部的值，更加整體和智能，在對(duì)比實(shí)驗(yàn)中得到的去噪效果較好。

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

算法的實(shí)驗(yàn)環(huán)境是Matlab R2015b 軟件，計(jì)算機(jī)處理器為Intel Core i7-3770，主頻3.40 GHz，內(nèi)存8 GB。用數(shù)學(xué)軟件進(jìn)行算法研究和數(shù)據(jù)分析，從評(píng)價(jià)圖像質(zhì)量的指標(biāo)峰值信噪比PSNR 和算法性能的時(shí)間復(fù)雜度兩方面進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)中，將提出的基于模擬退火算法的稀疏表示圖像去噪算法與K-SVD 算法進(jìn)行去噪效果對(duì)比。實(shí)驗(yàn)數(shù)據(jù)為采集的山西省人民醫(yī)院患者數(shù)據(jù)，采集設(shè)備是醫(yī)院現(xiàn)有X-Ray 薄掃設(shè)備。依次標(biāo)記20 組序列圖像，將每組序列圖像數(shù)據(jù)集作為一組測(cè)試集。現(xiàn)列出一組序列圖像的實(shí)驗(yàn)數(shù)據(jù)來(lái)說明算法[15]。

訓(xùn)練的字典如圖2 所示。圖2（a）是K-SVD算法的自適應(yīng)字典，圖2（b）是本文算法所得出的字典。

圖2 訓(xùn)練字典

序列圖像中其中一幅圖像去噪后的實(shí)驗(yàn)結(jié)果對(duì)比如圖3 所示。從圖3 可以看出，本文算法的實(shí)驗(yàn)結(jié)果有較高的峰值信噪比，從直觀感受方面看得到了很好的去噪效果。

圖3 噪聲圖像與去噪后的圖片對(duì)比

從整體看，一組序列圖像的實(shí)驗(yàn)數(shù)據(jù)如圖4、圖5 所示。

圖4 本文方法與K-SVD 方法的信噪比對(duì)比

圖5 本文方法和K-SVD 方法時(shí)間復(fù)雜度的對(duì)比

從實(shí)驗(yàn)數(shù)據(jù)可以看出，對(duì)于同一序列圖像，傳統(tǒng)的K-SVD 方法在序列圖像的信噪比上有間歇性的不穩(wěn)定，對(duì)于某些圖像的去噪效果不好，但本文方法有很好的信噪比與穩(wěn)定性。傳統(tǒng)的K-SVD 方法時(shí)間復(fù)雜度較高，而本文方法在時(shí)間復(fù)雜度上有所降低[16]。

接下來(lái)，擴(kuò)展至更多的數(shù)據(jù)集，使實(shí)驗(yàn)更具廣泛性和穩(wěn)定性[17-18]。以下是一組序列圖像的實(shí)驗(yàn)結(jié)果，運(yùn)行程序共運(yùn)行了20 組實(shí)驗(yàn)數(shù)據(jù)，現(xiàn)列出5組數(shù)據(jù)平均值，以說明實(shí)驗(yàn)結(jié)果及去噪效果。將本文算法與現(xiàn)在廣泛使用的K-SVD 的相關(guān)數(shù)據(jù)對(duì)比，結(jié)果如表1、表2 所示。表1 給出了其中5 組序列圖像的去噪峰值信噪比平均值?？梢钥闯?，在對(duì)薄掃圖像進(jìn)行去噪時(shí)，本文提出的方法具有較高的信噪比。

表1 峰值信噪比平均值對(duì)比 單位：dB

時(shí)間復(fù)雜度是評(píng)估算法性能的重要參數(shù)。實(shí)驗(yàn)中使用的算法對(duì)樣本的時(shí)間復(fù)雜度如表2 所示。在實(shí)驗(yàn)比較的算法中，本文算法的速度比K-SVD 算法快。在測(cè)試階段，本文算法識(shí)別一個(gè)手勢(shì)圖像序列消耗的平均時(shí)間為99.39 s。

在將快速去噪算法應(yīng)用于整個(gè)肺部結(jié)節(jié)的檢出實(shí)驗(yàn)中，根據(jù)先驗(yàn)知識(shí)與診斷結(jié)果，對(duì)是否檢出病灶的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比對(duì)，進(jìn)而算出結(jié)節(jié)識(shí)別率進(jìn)行對(duì)比，對(duì)于病灶不明顯的圖像具有一定的作用，實(shí)驗(yàn)數(shù)據(jù)存在偶然性，就目前所具有的數(shù)據(jù)，識(shí)別的結(jié)果有存在不一致的情況，這樣降低漏診率，在預(yù)處理方面對(duì)此課題有一定的輔助作用[19]。相對(duì)于傳統(tǒng)的K-SVD 方法的性能優(yōu)勢(shì)的描述，包括信噪比的提高、穩(wěn)定性的增強(qiáng)以及時(shí)間復(fù)雜度的降低。這些描述強(qiáng)調(diào)了本文方法在圖像去噪任務(wù)中的優(yōu)越性，并提供了實(shí)驗(yàn)數(shù)據(jù)支持，進(jìn)一步強(qiáng)調(diào)了本文方法的創(chuàng)新性和實(shí)用性，使讀者更容易理解本文方法是一個(gè)有價(jià)值的改進(jìn)，為圖像去噪帶來(lái)了重要的進(jìn)展。因此，本文提出的算法更適用于稀疏表示去噪，是一個(gè)行之有效的方法，在效果和性能方面都有一定的提高。

5 結(jié)語(yǔ)

將稀疏表示用于去噪，分為多尺度幾何分析、稀疏表示和優(yōu)化問題。其中，稀疏表示的兩大任務(wù)是字典的生成和信號(hào)的稀疏分解。運(yùn)用多尺度幾何分析將圖像轉(zhuǎn)換到某個(gè)域中，使稀疏盡可能地稀疏，能夠達(dá)到好的效果。本文運(yùn)用優(yōu)化方法進(jìn)行解優(yōu)化問題。在字典生成中，下一步的工作是用online 子空間聚類去批量去噪。對(duì)稀疏表示去噪還會(huì)繼續(xù)進(jìn)行研究，以為醫(yī)學(xué)診斷做好堅(jiān)實(shí)的前期工作。肺癌死亡率高的主要原因是其早期的細(xì)小病變特征很難被發(fā)現(xiàn)，臨床漏診率高。肺癌早期的檢出對(duì)提高治愈率起著至關(guān)重要的作用，能夠使達(dá)到醫(yī)學(xué)圖像診斷誤差率盡可能小。