王 康，趙 琦，李 銘

(1.南京理工大學 電子工程與光電技術(shù)學院，江蘇 南京 210094；2.中國科學院 蘇州生物醫(yī)學工程技術(shù)研究所，江蘇 蘇州 215163)

1 引 言

計算機斷層成像(Computed Tomography,CT)是通過無損方式獲取人體內(nèi)部解剖信息的一種醫(yī)學影像技術(shù)。CT掃描的普及，給臨床診斷帶來了極大方便，但也引發(fā)了受檢者的輻射損傷。由于X射線輻射劑量的累加效應，接受過量的X光照射可能會顯著增加病人誘發(fā)癌癥和白血病等疾病的概率。因此，國際輻射防護委員會建議，在滿足臨床診斷需求的情況下，CT輻射劑量越低越好[1-2]。

低劑量CT成像技術(shù)是減少輻射損傷的有效方法。臨床應用中的CT低劑量掃描成像可大致分為兩類：降低管電流或管電壓和降低采樣數(shù)[3-4]。前者多適用于投影數(shù)據(jù)采集完備的螺旋CT設備，后者則更適用于投影數(shù)據(jù)欠采樣的CBCT(Cone-beam CT,CBCT)設備。目前，國內(nèi)主流的CT設備仍采用濾波反投影法(Filtered Back-projection,FBP)進行重建，采用CT低劑量掃描協(xié)議將導致重建圖像中偽影和噪聲的增加以及對比度的下降，直接影響醫(yī)生對病變組織的診斷分析。因此，如何保證低劑量掃描下的圖像質(zhì)量是CT成像領(lǐng)域的熱點問題之一，而研究迭代重建算法[5-14]則是抑制噪聲和偽影，提高CT重建圖像質(zhì)量的關(guān)鍵。Idris等人[5]首先采用復合泊松分布模型對投影數(shù)據(jù)的噪聲特性進行建模，并提出基于這一模型的統(tǒng)計迭代重建算法(Statistical Iterative Reconstruction，SIR)。Lasio等[6]分析了基于泊松分布模型的統(tǒng)計迭代重建算法在實際CT系統(tǒng)中的性能，并驗證了泊松分布模型對于實際CT設備采集投影數(shù)據(jù)能取得較好的重建效果。在此基礎上，為了進一步提升重建效果，需改進目標函數(shù)中的先驗約束項。Zhang等[7]和Wu等[8]將L1范數(shù)字典學習(DL)正則化項引入統(tǒng)計迭代重建模型(SIR)，減輕了L2最小化帶來的過平滑效果并保留更多圖像細節(jié)。Gou等[9]測試了一種基于Lp字典學習(DL)稀疏約束的新穎重建算法，結(jié)論是較低的p值會呈現(xiàn)出更好的重建性能?；谧值鋵W習的CT低劑量重建對于圖像細節(jié)有所改善，但是其重建時間較長，且由于低劑量投影數(shù)據(jù)會導致迭代前期存在偽影，這些偽影會被字典所表示從而影響最終的重建結(jié)果。Xu等[10]則結(jié)合全變分最小化(Total Variation，TV)，提出了SIR-TV算法。該算法能較好地抑制噪聲和偽影，但對于稀疏性表征差的重建目標，重建圖像中易出現(xiàn)塊狀偽影。此后，Sidky等[11]又提出分數(shù)階全變分算法(TpV)。該算法通過調(diào)整p值可有效改進待重建目標的稀疏性表征，但是TpV算法易導致重建圖像中出現(xiàn)散點噪聲。Lei等[12]對TV結(jié)構(gòu)進行改進，通過引入自適應的加權(quán)權(quán)重，提出了自適應變權(quán)全變分算法(Adaptively Reweighted TV,ARWTV)。ARWTV算法對于噪聲和偽影有很好的抑制效果，但重建結(jié)果中弱對比度組織邊緣模糊。近年來，隨著深度學習的發(fā)展，許多學者將基于深度學習的方法應用于醫(yī)學領(lǐng)域[15-23]。與傳統(tǒng)方法相比，深度學習的圖像質(zhì)量有所提高，但這需要大量的訓練數(shù)據(jù)，且重建質(zhì)量取決于數(shù)據(jù)集的質(zhì)量。目前大多數(shù)深度學習中使用的數(shù)據(jù)是通過噪聲插入來模擬的，這可能無法反映低劑量CT掃描的實際噪聲分布。相反，傳統(tǒng)的基于壓縮感知(Compressed Sensing ,CS)的重建方法不需要大量的訓練數(shù)據(jù)，并且可以很容易地集成到傳統(tǒng)的迭代重建框架中。

針對CT低劑量掃描下出現(xiàn)的噪聲和偽影問題，文中提出一種基于分數(shù)階全變分動態(tài)優(yōu)化(TpV dynamic optimization,TpV-DO)的CT低劑量成像算法。首先通過一族雙曲正切函數(shù)集構(gòu)造TpV的動態(tài)復合函數(shù)模型，以增強待重建目標的稀疏性表征；其次，通過引入動態(tài)優(yōu)化設計以實現(xiàn)對噪聲和偽影的抑制以及對弱對比度組織邊緣的保護；最后，結(jié)合統(tǒng)計迭代重建框架，使用交替優(yōu)化算法實現(xiàn)對目標函數(shù)的優(yōu)化求解。采用數(shù)值模型和動物掃描獲取的投影數(shù)據(jù)對新算法進行了實驗，并與FBP算法、TpV算法、ARWTV算法進行了比較。

2 算法及理論分析

2.1 統(tǒng)計迭代成像模型

基于泊松分布模型來模擬投影數(shù)據(jù)采集過程，將統(tǒng)計迭代模型與正則項結(jié)合的目標函數(shù)可以表示為[10]：

(1)

2.2 基于TpV的動態(tài)優(yōu)化模型描述

在2D圖像空間中，基于TpV模型的稀疏性變換可以用雙索引的方式表示如下[11]：

(2)

其中：

(3)

其中：j=(m-1)×W+n，m=1,…,H，n=1,…,W，H和W分別代表重建圖像x的寬度和高度；p為彈性的范數(shù)模型控制參數(shù)，通過p值調(diào)整可以實現(xiàn)范數(shù)模型在L0范數(shù)與L2范數(shù)之間選擇。

壓縮感知理論指出[24]：重建信號的稀疏性表示越好，越有利于精確地恢復原始信號。為了應用TpV模型增強重建圖像的稀疏表征，并建立動態(tài)優(yōu)化模型。現(xiàn)考慮如下一族雙曲正切函數(shù)集，其具體表示形式如下：

(4)

其對應的范數(shù)求和形式如下：

(5)

其中：σ為動態(tài)調(diào)整參數(shù)，用于控制圖像表征和相鄰迭代過程的近似性。

如圖1所示，當σ參數(shù)取值較大時，對重建圖像的稀疏性表征較差，重建結(jié)果容易導致過平滑，圖像紋理細節(jié)模糊；當σ參數(shù)取值較小時，對重建圖像的稀疏性表征較好，重建結(jié)果容易陷入局部極值，圖像噪聲水平顯著?；谏鲜鲭p曲正切函數(shù)集構(gòu)造TpV動態(tài)復合函數(shù)模型，其具體形式如下：

圖1 雙曲正切函數(shù)曲線動態(tài)演化圖Fig.1 Dynamic evolution of hyperbolic tangent function curve

(6)

其中：p和σ均為動態(tài)優(yōu)化參數(shù)，對式(6)分析，當p趨近于2，σ取較大值時，R(x)接近于L2范數(shù)，對于偽影和噪聲抑制效果較好；當σ趨近于0，p取較小值時，R(x)接近于理想的L0范數(shù)，對圖像紋理細節(jié)的保護效果較好。因此通過對p和σ的合理選擇從而在保持圖像邊緣和細節(jié)的同時能有效克服TV正則項存在的階梯效應問題。

為了使上述的正則化估計同時表現(xiàn)出較好的紋理保護和噪聲抑制性能，本文提出動態(tài)最優(yōu)化過程，即初始迭代優(yōu)化時，采用較大的p和σ值以更好地抑制偽影和噪聲，接下來的優(yōu)化過程中依次降低p和σ值以更好地保存圖像紋理信息。

實驗中，為了保證算法初始時具有較好的偽影和噪聲抑制效果，初始設定：σ0=0.8，p0=1.2；為了保證相鄰兩次迭代的近似性，相鄰迭代的遞減因子Δσ和Δp分別設定為0.016和0.01；同時為了避免算法收斂于局部極值，實驗中設定：σmin≥0.1，pmin≥0.8。

2.2 交替最小化DO-TpV重建算法

如式(1)為典型的統(tǒng)計重建目標函數(shù)，該目標函數(shù)由對數(shù)似然數(shù)據(jù)保真項和正則化項兩部分構(gòu)成，本文采用交替最小化的方法對其進行優(yōu)化求解，具體算法流程如下：

步驟一、最小化對數(shù)似然數(shù)據(jù)保真項：采用典型的可分離拋物面算法[10]進行優(yōu)化，其迭代更新公式為：

(7)

式中：t表示迭代次數(shù)，x0為初始迭代圖像。

步驟二、式(6)的動態(tài)最優(yōu)化求解：采用梯度下降法對其進行優(yōu)化。其中式(6)的梯度公式可以表示為：

(8)

其中：TpV模型的梯度公式可表示為：

(9)

式中：ζ為擾動約束因子，通常取ζ=1×10-8。

基于上述的梯度公式，采用梯度下降法進行迭代更新，其計算公式為：

(10)

上述算法具體實施中，首先進行參數(shù)x、σ、p的初始化，接下來重復執(zhí)行步驟一與步驟二直至滿足預先設定的迭代終止條件。下文實驗中將迭代終止條件設定為t不大于預先設定的迭代總次數(shù)Tmax或相鄰兩次迭代誤差小于設定閾值err，如果t>Tmax，則迭代終止；否則，t=t+1 ，更新σ和p值，并啟動下一輪迭代。本文提出的基于分數(shù)階全變分動態(tài)優(yōu)化算法(DO-TpV)的CT低劑量成像重建流程如下所示：

輸入:A,l^,β,σ,p,err,初始化x0,σ0,p0,Δσ,Δp,t=1;

輸出:重建圖像x.While(t‰Tmaxor‖xt-xt-1‖/‖xt‖≥err){步驟一:優(yōu)化SIR過程xtj=xt-1j-∑Ndi=1ai,jyi[Axt-1]i-^li ∑Ndi=1ai,jyi∑Jk=1ai,k ,j=1,...,J;步驟二:分數(shù)階全變分動態(tài)模型優(yōu)化(TpV-DO)d=‖xt-xt-1‖2;dx=?R(x)?xj;xtj=xtj-β·d·dx;步驟三:參數(shù)更新t=t+1;ifσ≥0.1,p≥0.8σ=σ-Δσ,p=p-Δp.}

3 實驗與結(jié)果

為了驗證本文算法在噪聲抑制和軟組織細節(jié)保護方面的有效性，分別采用數(shù)值仿真實驗和動物掃描實驗對本文提出的CT低劑量成像算法(DO-TpV)進行實驗驗證，并對其與解析重建算法(FBP)、TpV算法、字典學習(DL)、ARWTV算法進行了比較。其中，算法的優(yōu)化參數(shù)設定：本文所提算法采用2.2節(jié)所述方式；TpV算法正則化參數(shù)β=0.2，動態(tài)優(yōu)化參數(shù)p=0.8～1.2；ARWTV算法正則化參數(shù)β=0.2，變?nèi)螖?shù)設定為5；字典學習(DL)算法正則化參數(shù)β=0.2，圖像塊尺寸設定為7×7，過完備字典的稀疏度設定為4。

3.1 數(shù)值模擬實驗

本部分實驗采用數(shù)值模型生成圓軌跡扇形束掃描投影數(shù)據(jù)，投影參數(shù)設定如下：X射線源到旋轉(zhuǎn)中心的距離為54.1 cm，探測器到旋轉(zhuǎn)中心的距離為40.8 cm，CT圖像的像素數(shù)為512×512，尺寸大小為20 cm×20 cm，探測器數(shù)為642，相鄰探測器間距為0.672 mm。掃描范圍為0°～360°，步長為1°和2°，以生成180°和360°兩種采樣角度數(shù)下的投影數(shù)據(jù)。實驗中采用Patrick等[10]提出的Poisson分布模擬光子探測的隨機過程，其中：實驗總光子數(shù)設定為8.0×105。通過對探測器記錄信號進行對數(shù)化處理后，獲得含噪聲的投影數(shù)據(jù)。分別使用FBP算法、TpV算法、字典學習(DL)算法、ARWTV算法和DO-TpV算法進行CT圖像重建。圖2和圖3分別給出了數(shù)值模型在180°和360°兩種采樣數(shù)下的實驗結(jié)果，窗口顯示范圍[0.185,0.205]。

圖3 數(shù)值模型實驗結(jié)果(360采樣數(shù))Fig.3 Numerical model results (360 views)

如圖2和圖3所示，在低劑量掃描條件下，F(xiàn)BP算法重建圖像中出現(xiàn)嚴重的偽影和噪聲，且軟組織細節(jié)完全被噪聲和偽影遮蓋，總體圖像質(zhì)量極差。對比解析算法(FBP)結(jié)果，4種迭代算法重建圖像的質(zhì)量有了明顯提升。對比4種迭代算法結(jié)果，ARWTV算法取得了很好的偽影和噪聲抑制效果，但重建圖像軟組織細節(jié)模糊嚴重，如圖2(e)紅色箭頭所示的放大區(qū)域。相比ARWTV算法，字典學習算法(DL)和TpV算法在軟組織細節(jié)保護方面有了顯著的提升，但字典學習(DL)算法的部分軟組織產(chǎn)生變形模糊，如圖2(d)紅色箭頭所示，且字典學習(DL)從整體看尤其在180°較低采樣數(shù)下會存在條紋狀偽影，而對于TpV算法部分軟組織細節(jié)結(jié)構(gòu)不清晰，如圖2(c)紅色箭頭所示的放大區(qū)域。由圖2和圖3可以看出，DO-TpV算法在軟組織細節(jié)保護方面表現(xiàn)最好，且取得較好的偽影和噪聲抑制效果，總體圖像質(zhì)量的提升顯著。

為了評價不同算法的重建結(jié)果對高對比度邊緣細節(jié)的保護能力，本文提取圖2紅色實線所示位置，并采用剖面密度曲線顯示TpV算法、字典學習(DL)算法、ARWTV算法和DO-TpV算法的重建結(jié)果與金標準圖像的接近程度。圖4實驗結(jié)果表明，與其他算法的重建結(jié)果相比，DO-TpV算法的重建結(jié)果更好地擬合了金標準圖像的密度曲線。

圖4 重建結(jié)果的剖面灰度值比較圖(圖2)Fig.4 Line profiles of reconstructed results of Fig.2

本文采用歸一化平均絕對偏差(Normalized mean absolute deviation,NMAD)和信噪比(Signal noise ratio,SNR)2個標準定量評價重建結(jié)果，NMAD和SNR的定義如下：

(11)

(12)

由表 2 可以看出，DO-TpV算法的NMAD小于FBP算法、TpV算法、字典學習(DL)算法、ARWTV算法的NMAD，表明DO-TpV算法重建的圖像與理想圖像更接近。在180°視角和360°視角的數(shù)值模型實驗中，DO-TpV算法重建圖像的信噪比比FBP算法、TpV算法、字典學習(DL)算法、ARWTV算法要高，表明DO-TpV算法能更好地保護圖像紋理信息和抑制圖像噪聲，重建圖像的質(zhì)量更高。

表2 數(shù)值模型定量分析Tab.2 Quantitative evaluation of the numerical results

3.2 動物數(shù)據(jù)實驗

為了進一步驗證DO-TpV算法的性能，采用真實CT系統(tǒng)掃描小鼠來生成投影數(shù)據(jù)。系統(tǒng)掃描參數(shù)設定如下：X光管電壓50 kVp，X光管電流1 mA，曝光時間0.466 9 s，X射線源到掃描中心距離為22.188 cm，X射線源到探測器中心的距離為65.85 cm，平板探測器陣列數(shù)1 536×880，探測單元尺寸0.15 mm×0.15 mm，掃描間隔為2°，共采集180個等間距投影圖像。分別使用FBP算法、TpV算法、字典學習(DL)算法、ARWTV算法和DO-TpV算法對中心片層投影數(shù)據(jù)進行重建實驗。重建像素大小為512×512，重建像素單元尺寸為0.1 mm×0.1 mm。圖5給出小鼠掃描數(shù)據(jù)重建結(jié)果，窗口顯示范圍為[-0.09，1.19]。

圖5 小鼠掃描數(shù)據(jù)重建結(jié)果Fig.5 Reconstructed results of scanned mouse datasets

從圖5可以看出，低劑量掃描協(xié)議下，解析算法(FBP)重建圖像中出現(xiàn)嚴重的偽影和噪聲，且軟組織邊緣細節(jié)極其模糊，如圖5(b)放大區(qū)域所示。對比解析算法重建結(jié)果，4種迭代算法重建結(jié)果在抑制偽影和噪聲，保護軟組織結(jié)構(gòu)方面取得了明顯改善。對比4種迭代算法結(jié)果，ARWTV算法在抑制偽影和噪聲方面效果顯著，但重建圖像部分細節(jié)模糊不清晰，如圖5(c)放大區(qū)域紅色箭頭所示位置。相比ARWTV算法，TpV算法能夠重建出圖像的邊緣細節(jié)，但重建結(jié)果中存在一些散點噪聲，如圖5(c)放大區(qū)域所示。而對于字典學習(DL)方法，偽影和噪聲得到了一定程度的抑制，但結(jié)果仍然存在由于過完備字典而導致的不規(guī)則偽影和噪聲，如圖5(d)放大區(qū)域所示。對比ARWTV算法和TpV算法的重建圖像，DO-TpV算法重建的圖像細節(jié)更為清晰，且圖像的噪聲水平也顯著降低，總體圖像質(zhì)量最佳。

4 結(jié) 論

本文提出一種基于動態(tài)優(yōu)化的CT低劑量成像算法。該算法通過雙曲正切函數(shù)集構(gòu)造分數(shù)階變分的動態(tài)復合函數(shù)模型，使基于分數(shù)階變分的重建逐步接近理想的L0范數(shù)重建。通過動態(tài)最優(yōu)化過程逐次增強當前重建圖像的稀疏表征，改善了DO-TpV算法在圖像紋理細節(jié)保護和噪聲抑制方面的性能，提升了重建圖像的質(zhì)量。數(shù)值模型實驗結(jié)果表明，在180個采樣角度下，文中算法重建圖像的信噪比分別比FBP算法、TpV算法、字典學習(DL)方法、ARWTV算法高出29.51，8.03，9.15，6.81 dB。動物數(shù)據(jù)實驗結(jié)果表明，DO-TpV算法重建結(jié)果有效抑制了噪聲和偽影，清晰重建出小鼠軟組織細節(jié)，極大地提高了低劑量采集數(shù)據(jù)重建圖像的質(zhì)量。