張鵬程，謝世朋

(南京郵電大學(xué) 通信與信息工程學(xué)院，江蘇 南京 210003)

0 引 言

輻射劑量問題已成為CT領(lǐng)域亟待解決的熱點問題，當(dāng)前的研究主要集中在降低電流和減少投影數(shù)方面。降低電流的方法，雖然降低了X射線的輻射量，但同時引入了較多的統(tǒng)計噪聲。為了去除低劑量CT中的統(tǒng)計噪聲，L.L.Chen等提出了使用改進的BM3D算法來提高圖像的質(zhì)量[1]，而SayedMasoud Hashemi等采用了NCRE的迭代方法減小低劑量CT中的噪聲[2]。減少投影數(shù)的方法可以減少CT中的大量輻射，但同時給圖像帶來了嚴重的偽影。為了去除重建中的偽影，Shaoting Zhang等利用字典學(xué)習(xí)的方法稀疏重建[3]，張瀚銘等利用非局部的TV正則化稀疏重建[4]，張俊峰等利用伽瑪正則化稀疏重建[5]，在圖像去偽影方面取得了良好的效果。由于數(shù)據(jù)噪聲的存在和不理想的先驗懲罰項的選取，使得基于TV正則化的重建結(jié)果存在梯形偽影和片狀偽影。在臨床診斷過程中，片狀偽影會影響到診斷的誤判和正確性。其它稀疏重建的方法都需要大量的正則化約束，且恢復(fù)出來的圖像有較低的PSNR(峰值信噪比)和SSIM(結(jié)構(gòu)相似度)值，仍然存在較嚴重的偽影。卷積神經(jīng)網(wǎng)絡(luò)(CNN)已經(jīng)在圖像去噪、圖像識別和圖像分類領(lǐng)域取得了很好的效果且誤差率都很小，在醫(yī)學(xué)圖像處理中得到了廣泛地應(yīng)用。針對這一現(xiàn)象，本文提出了基于殘差學(xué)習(xí)的CT稀疏重建偽影校正的方法，能更好估計這些偽影的特征達到去偽影的目的。

1 方 法

本文提出了基于殘差學(xué)習(xí)的CT稀疏重建偽影校正方法，首先通過FBP(濾波反投影)方法對CT產(chǎn)生的稀疏投影數(shù)據(jù)進行重建，重建后的圖像帶有嚴重的偽影；然后通過建立殘差神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)對偽影的特征進行學(xué)習(xí)，得到偽影圖；最后通過稀疏重建的圖與偽影的圖做殘差，恢復(fù)出清晰的CT圖像。

1.1 FBP重建

濾波反投影重建(filtered back-projection，F(xiàn)BP)算法具有重建速度快、易于硬件實現(xiàn)等優(yōu)點,當(dāng)投影數(shù)據(jù)比較完備時,能快速重建出高質(zhì)量圖像,因此它是目前應(yīng)用最為廣泛的一種重建算法。對于平行束CT系統(tǒng)，遵循bear定律。這個投影的過程實質(zhì)上是一個拉東變化的過程。探測器的距離為t，投影角度為θ，公式可以表示為

(1)

f(x,y)代表當(dāng)前的圖像，當(dāng)t=xcosθ+ysinθ，采樣投影數(shù)據(jù)充足那么FBP則有

(2)

其中，w表示斜坡濾波器，Pθ(w)表示沿著探測器方向的一維傅里葉變換，wPθ(w)就是濾波的投影。第二重積分本質(zhì)上是‘濾波反投影’，即對濾波投影的數(shù)據(jù)進行傅里葉反變換，然后累積得到反投影的圖像。

使用CT對患者進行不同投影數(shù)的掃描，產(chǎn)生60個(3600圓軌道掃描，掃描間隔為60)、90個(3600圓軌道掃描，掃描間隔為40)、120個投影(3600圓軌道掃描，掃描間隔為30)和全投影(3600圓軌道掃描，掃描間隔為10)的數(shù)據(jù)，而后對這些數(shù)據(jù)進行FBP反投影重建。由CT掃描患者肺部切片全投影成像和稀疏重建成像如圖1所示。

圖1 全投影和稀疏投影重建

FBP稀疏重建的CT圖像存在明顯的條狀偽影，投影數(shù)減少能夠降低X射線對人體造成的輻射危害，但同時帶來了偽影的增加，為了去除偽影引入了基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)。

1.2 基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)

由于深度學(xué)習(xí)在學(xué)習(xí)圖像噪聲特點方面取得了很好的效果，張凱等[6]利用殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)去高斯噪聲、解決SISR和JPEG圖像塊效應(yīng)等方面較傳統(tǒng)的方法效果更好。另外，Yo Seob Han等[7]根據(jù)張凱的思路提出了基于殘差學(xué)習(xí)的壓縮感知CT重建方法，該方法在去稀疏重建偽影方面同樣取得了很好的效果?；诙叩乃枷?，本文利用基于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)對稀疏重建的圖像進行學(xué)習(xí)，然后通過稀疏重建的圖像減去學(xué)習(xí)得到的偽影圖像，最后恢復(fù)出清晰的校正圖像。殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)的輸入為稀疏重建圖像f=x+n。在深度學(xué)習(xí)卷積網(wǎng)絡(luò)當(dāng)中，典型的去噪模型如MLP[8]和CSF[9]旨在學(xué)習(xí)一個映射函數(shù)F(f)=x預(yù)測一個清晰的圖像。對于殘差學(xué)習(xí)的卷積神經(jīng)網(wǎng)絡(luò)，訓(xùn)練的是一個映射μ(f)≈n函數(shù)。本文提出的基于殘差的卷積神經(jīng)網(wǎng)絡(luò)旨在解決以下問題[8]

(3)

(4)

(5)

該殘差學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)即為本文提出的去偽影校正方法，不同于何愷明等[10]提出的用多個殘差單元預(yù)測殘差圖像，本文只使用單個殘差單元預(yù)測殘差圖像。由殘差學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)得到的殘差圖即為偽影特征圖，然后根據(jù)x=f-n，得到最終清晰的圖像，實現(xiàn)過程如圖3所示。

圖2 殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

圖3 稀疏重建圖中去偽影恢復(fù)CT的過程

1.2.1 殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的具體結(jié)構(gòu)

該網(wǎng)絡(luò)的深度為D，它有3種不同的網(wǎng)絡(luò)層，分別為：卷積網(wǎng)絡(luò)+ReLU[11](修正線性單元)、卷積網(wǎng)絡(luò)+批規(guī)范化(batch normalization，BN)和卷積網(wǎng)絡(luò)層。第一層為卷積網(wǎng)絡(luò)和ReLU構(gòu)成的網(wǎng)絡(luò)層，它由64個3*3*1的濾波器組成，用來產(chǎn)生64個特征映射。這一層加入ReLU，目的是在卷積網(wǎng)絡(luò)中引入非線性，因為神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)的實際數(shù)據(jù)是非線性的，卷積是一個線性操作——元素級別的矩陣相乘和相加，因此需要通過使用非線性函數(shù)ReLU來引入非線性。然而經(jīng)實踐證明，ReLU具備引導(dǎo)適度稀疏的能力，訓(xùn)練后的網(wǎng)絡(luò)完全具備適度的稀疏性，而且訓(xùn)練后的可視化效果和傳統(tǒng)預(yù)處理的效果很相似。正因為它有這樣的作用使得ReLU的速度非常快，而且精確度很高。第二層到D-1層，在卷積網(wǎng)絡(luò)和修正線性單元的基礎(chǔ)上引入了BN(批規(guī)范)[12]，使用64個3*3*64的濾波器，將BN加在卷積和ReLU之間。深層神經(jīng)網(wǎng)絡(luò)在做非線性變換前激活輸入值隨著網(wǎng)絡(luò)深度的加深，分布逐漸發(fā)生偏移或者變動，逐漸往非線性函數(shù)取值區(qū)間的上下限兩端靠近，這導(dǎo)致后向傳播時低層神經(jīng)網(wǎng)絡(luò)梯度的消失，進而使得深層神經(jīng)網(wǎng)絡(luò)的訓(xùn)練速度收斂越來越慢。BN通過一定的規(guī)范化手段，把每層神經(jīng)網(wǎng)絡(luò)中任意神經(jīng)元輸入值的分布強行拉回到均值為0、方差為1的標(biāo)準正態(tài)分布中(即把越來越偏的分布強制拉回比較標(biāo)準的分布中)，使得激活輸入值落在非線性函數(shù)對輸入比較敏感的區(qū)域。這樣輸入的小變化引起損失函數(shù)較大的變化，梯度隨之變大，于是避免了梯度消失的問題，同時學(xué)習(xí)收斂速度加快，大大提高了訓(xùn)練速度。第三層為一個卷積層，使用3*3*64的濾波器用來重建輸出結(jié)果。

1.2.2 優(yōu)化加速

在訓(xùn)練中優(yōu)化過程一般采用梯度下降法，梯度下降法有隨機梯度下降法(stochastic gradient descent，SGD)和Adam算法等。其中Adam算法根據(jù)損失函數(shù)對每個參數(shù)梯度的一階矩估計和二階矩估計進行動態(tài)調(diào)整，使其調(diào)整到一個合適的學(xué)習(xí)速率。Adam算法在每次迭代中參數(shù)的學(xué)習(xí)步長都有一個確定的范圍，不會因為很大的梯度導(dǎo)致學(xué)習(xí)步長的增大，參數(shù)的值也比較穩(wěn)定。批量梯度下降算法的缺點是每次都會使用全部訓(xùn)練樣本，而且每次都使用完全相同的樣本集，這會使得計算相對的冗雜。相反隨機梯度下降算法的實現(xiàn)是可以在線更新的，每次只隨機選擇一個樣本來更新模型參數(shù)，學(xué)習(xí)速度是非?？斓?。在神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程當(dāng)中，SGD和Adam優(yōu)化方法的引入使網(wǎng)絡(luò)結(jié)構(gòu)訓(xùn)練得到了更好的結(jié)果，其中Adam算法在圖像高斯去噪過程中效果要強與SGD算法[7]。

2 結(jié) 果

2.1 實驗數(shù)據(jù)

實驗的數(shù)據(jù)由南京某醫(yī)院提供，來自10個患者肺部的CT掃描。實驗設(shè)備為臨床使用的西門子SOMATOM Sensation 16排CT，曝光量為100 mAs，管電壓為120 kVp,射線源距離探測器陣列為1040 mm，射線源到旋轉(zhuǎn)中心的距離為570 mm。重建的CT圖像大小為512像素×512像素，每個像素的大小為1.2 mm×1.2 mm。

設(shè)定不同的間隔掃描角度，從CT掃描中分別得獲得全投影(360個投影)、120個、90個和60個投影數(shù)的原始數(shù)據(jù)，利用FBP方法對這些數(shù)據(jù)重建得到對應(yīng)的重建圖像。全投影重建的圖像作為完整的圖像，用于和稀疏投影重建的圖像作對比。針對某一個投影數(shù)的重建圖像，從重建的圖像中隨機選取400張512*512的圖像作為訓(xùn)練數(shù)據(jù)，再從剩余中選取20張512*512的圖像作為測試數(shù)據(jù)。對于殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)的輸入，使用全投影圖像和稀疏重建后圖像之間的殘差作為訓(xùn)練初始化值，同時對不同投影數(shù)的重建圖像分別進行訓(xùn)練，經(jīng)過20次訓(xùn)練后的殘差圖作為輸出。

2.2 網(wǎng)絡(luò)訓(xùn)練和參數(shù)的設(shè)定

本文為了訓(xùn)練的準確性而得到較好的去偽影效果，不同投影數(shù)重建的訓(xùn)練樣本分別為400張，測試圖像分別為20張，同時把網(wǎng)絡(luò)的深度D設(shè)置為20。采用損失函數(shù)學(xué)習(xí)殘差映射μ(f)估計預(yù)測殘差n。優(yōu)化方法選擇Adam下降法，梯度的一階矩估計α為0.01，二階矩估計β1為0.9，β2為0.999，更新幅度的校正參數(shù)epsilon為1e-8。該殘差網(wǎng)絡(luò)的訓(xùn)練速率為10-3～10-4，每次訓(xùn)練逐漸遞減。為了產(chǎn)生較好的去偽影效果，提取出條狀偽影的特征圖，訓(xùn)練次數(shù)設(shè)置為60。

本文在訓(xùn)練的過程中使用MatConvNet工具包，運行環(huán)境為MATLAB 2015a，處理器為Intel(R) Core(TM) i7_2600 CPU@ 3.4 GHz，內(nèi)存為8 GB，GPU顯卡為GeForce GTX 1070。在此配置環(huán)境中，通過該網(wǎng)絡(luò)訓(xùn)練樣本大約需要24個小時。對不同的投影數(shù)重建，通過測試程序恢復(fù)出重建的圖像，所需的測試時間見表1。

表1 不同投影數(shù)經(jīng)過稀疏重建恢復(fù)過程的測試時間

2.3 實驗結(jié)果

2.3.1 定性分析

分別利用120個投影、90個投影和60個投影重建的圖像作為樣本進行訓(xùn)練，而后通過殘差圖像恢復(fù)清晰的CT圖。圖4顯示了不同投影重建圖像訓(xùn)練后恢復(fù)的結(jié)果，應(yīng)用本文提出的方法達到了很好的去偽影效果。

圖4 不同投影數(shù)重建圖像經(jīng)過殘差學(xué)習(xí)恢復(fù)的結(jié)果

隨著稀疏重建投影數(shù)的減少，仍能很好地恢復(fù)出重建的圖像，效果接近全投影的原圖像。使用殘差網(wǎng)絡(luò)訓(xùn)練方法達到了預(yù)期的處理目的了。

2.3.2 定量分析

為了驗證本文提出的方法，觀察圖像像素強度尤為重要。本實驗采用120個投影的重建圖像、90個投影的重建圖像和60個投影的重建圖像與全投影、恢復(fù)出的圖像作對比。選取全投影重建的圖像與對應(yīng)的不同投影個數(shù)重建的圖像和恢復(fù)出的圖像，測試其在同一幅CT圖像中200行，120列到250列像素的強度。

圖5中120個、90個和60個投影重建的圖像強度偏離了全投影的像素強度，尤其是60個投影重建的圖像偏離的比較嚴重，通過本文提出的方法恢復(fù)的圖像強度已經(jīng)很接近全投影的圖像強度了。

圖5 全投影、不同投影個數(shù)的重建和恢復(fù)的圖像(圖像矩陣的第200行)在不同像素位置的強度

2.3.3 圖像評價

在評價重建圖像的好壞主要有兩個指標(biāo)：峰值信噪比(peak signal to noise ratio，PSNR)和結(jié)構(gòu)相似性(structu-ral similarity，SSIM)[13]。

2.3.3.1 PSNR

PSNR的計算如下

(6)

(7)

MSE表示稀疏重建圖像X和通過殘差學(xué)習(xí)恢復(fù)的圖像Y的均方誤差(mean square error)，H、W分別為圖像的高度和寬度；n為像素的位深，在CT圖像中位深為12。

2.3.3.2 SSIM

SSIM是一種全參考的圖像質(zhì)量評價指標(biāo)，它分別從亮度、對比度、結(jié)構(gòu)3方面度量圖像相似性

SSIM(X,Y)=[l(X,Y)]α·[c(X,Y)]β·[s(X,Y)]γ

(8)

l(X,Y)為亮度對比函數(shù)，c(X,Y)為對比度對比函數(shù)，s(X,Y)為結(jié)構(gòu)對比度函數(shù)，X和Y分別表示稀疏重建的圖像和殘差學(xué)習(xí)恢復(fù)的圖像，α、β、γ是3個對比函數(shù)所占的比例因子。

本實驗選取了20個原始數(shù)據(jù)作為測試數(shù)據(jù)，針對不同投影個數(shù)計算使用本文方法與全投影圖像比較的平均PSNR和SSIM值，結(jié)果見表2。

表2 不同投影數(shù)稀疏重建后使用殘差學(xué)習(xí)方法與全投影重建圖像比較的PSNR和SSIM平均值

表2反映的是訓(xùn)練恢復(fù)出的圖像與全投影重建后的圖像峰值信噪比和結(jié)構(gòu)相似度值。由上表的數(shù)據(jù)對比，隨著投影數(shù)的減少重建效果在一定程度上下降了。PSNR介于20 dB～40 dB之間，效果是最佳的。當(dāng)投影數(shù)太少，PSNR較小，效果要變得稍微差一些。SSIM的值也在逐漸下降。其中SSIM的值越接近于1，效果越好。

2.3.4 實驗比較

為了驗證本文提出方法的效果，實驗部分采用了與傳統(tǒng)方法的比較。使用張俊峰等的伽瑪正則化[14]的CT稀疏重建方法分別對120個、90個和60個投影的數(shù)據(jù)進行重建，觀察重建圖的效果和相應(yīng)的PSNR、SSIM值。圖6從左到右依次是稀疏投影重建的圖像、使用伽瑪正則化方法重建的圖像和本文提出的方法重建的圖像。

從測試數(shù)據(jù)的平均PSNR和SSIM分析，可以明顯的得出本文提出的方法要好于傳統(tǒng)的伽瑪正則化方法，見表3、表4。在一定程度上，利用殘差學(xué)習(xí)的方法去偽影較伽瑪正則化方法好。

表3 不同投影數(shù)的平均PSNR/dB比較

表4 不同投影數(shù)的平均SSIM比較

圖6 伽瑪正則化方法和本文方法的比較

3 結(jié)束語

與傳統(tǒng)方法不同的是本文研究的方法在稀疏重建的基礎(chǔ)上引入了殘差學(xué)習(xí)的結(jié)構(gòu)框架，在一定的程度上大大降低了X射線的輻射劑量。同時采用GPU加速，縮短了訓(xùn)練的時間，加速了實驗的訓(xùn)練過程。

總之，本文提出的校正方法不需要大面積掃描人體，在現(xiàn)有的條件下不增加設(shè)備，不影響圖像的成像質(zhì)量；降低了輻射對人體不必要的傷害，達到了很好的效果。另外，該殘差學(xué)習(xí)卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)比較單一，需要進一步優(yōu)化和改進，這將是本文下一步的重點工作。