基于動(dòng)態(tài)優(yōu)化的X射線CT低劑量重建

王 康，趙 琦，李 銘

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

1 引 言

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

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

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

2 算法及理論分析

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

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

(1)

2.2 基于TpV的動(dòng)態(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)整可以實(shí)現(xiàn)范數(shù)模型在L0范數(shù)與L2范數(shù)之間選擇。

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

(4)

其對(duì)應(yīng)的范數(shù)求和形式如下：

(5)

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

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

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

(6)

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

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

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

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

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

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

(7)

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

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

(8)

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

(9)

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

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

(10)

上述算法具體實(shí)施中，首先進(jìn)行參數(shù)x、σ、p的初始化，接下來重復(fù)執(zhí)行步驟一與步驟二直至滿足預(yù)先設(shè)定的迭代終止條件。下文實(shí)驗(yàn)中將迭代終止條件設(shè)定為t不大于預(yù)先設(shè)定的迭代總次數(shù)Tmax或相鄰兩次迭代誤差小于設(shè)定閾值err，如果t>Tmax，則迭代終止；否則，t=t+1 ，更新σ和p值，并啟動(dòng)下一輪迭代。本文提出的基于分?jǐn)?shù)階全變分動(dòng)態(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;步驟二:分?jǐn)?shù)階全變分動(dòng)態(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 實(shí)驗(yàn)與結(jié)果

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

3.1 數(shù)值模擬實(shí)驗(yàn)

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

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

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

為了評(píng)價(jià)不同算法的重建結(jié)果對(duì)高對(duì)比度邊緣細(xì)節(jié)的保護(hù)能力，本文提取圖2紅色實(shí)線所示位置，并采用剖面密度曲線顯示TpV算法、字典學(xué)習(xí)(DL)算法、ARWTV算法和DO-TpV算法的重建結(jié)果與金標(biāo)準(zhǔn)圖像的接近程度。圖4實(shí)驗(yàn)結(jié)果表明，與其他算法的重建結(jié)果相比，DO-TpV算法的重建結(jié)果更好地?cái)M合了金標(biāo)準(zhǔn)圖像的密度曲線。

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

本文采用歸一化平均絕對(duì)偏差(Normalized mean absolute deviation,NMAD)和信噪比(Signal noise ratio,SNR)2個(gè)標(biāo)準(zhǔn)定量評(píng)價(jià)重建結(jié)果，NMAD和SNR的定義如下：

(11)

(12)

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

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

3.2 動(dòng)物數(shù)據(jù)實(shí)驗(yàn)

為了進(jìn)一步驗(yàn)證DO-TpV算法的性能，采用真實(shí)CT系統(tǒng)掃描小鼠來生成投影數(shù)據(jù)。系統(tǒng)掃描參數(shù)設(shè)定如下：X光管電壓50 kVp，X光管電流1 mA，曝光時(shí)間0.466 9 s，X射線源到掃描中心距離為22.188 cm，X射線源到探測(cè)器中心的距離為65.85 cm，平板探測(cè)器陣列數(shù)1 536×880，探測(cè)單元尺寸0.15 mm×0.15 mm，掃描間隔為2°，共采集180個(gè)等間距投影圖像。分別使用FBP算法、TpV算法、字典學(xué)習(xí)(DL)算法、ARWTV算法和DO-TpV算法對(duì)中心片層投影數(shù)據(jù)進(jìn)行重建實(shí)驗(yàn)。重建像素大小為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)嚴(yán)重的偽影和噪聲，且軟組織邊緣細(xì)節(jié)極其模糊，如圖5(b)放大區(qū)域所示。對(duì)比解析算法重建結(jié)果，4種迭代算法重建結(jié)果在抑制偽影和噪聲，保護(hù)軟組織結(jié)構(gòu)方面取得了明顯改善。對(duì)比4種迭代算法結(jié)果，ARWTV算法在抑制偽影和噪聲方面效果顯著，但重建圖像部分細(xì)節(jié)模糊不清晰，如圖5(c)放大區(qū)域紅色箭頭所示位置。相比ARWTV算法，TpV算法能夠重建出圖像的邊緣細(xì)節(jié)，但重建結(jié)果中存在一些散點(diǎn)噪聲，如圖5(c)放大區(qū)域所示。而對(duì)于字典學(xué)習(xí)(DL)方法，偽影和噪聲得到了一定程度的抑制，但結(jié)果仍然存在由于過完備字典而導(dǎo)致的不規(guī)則偽影和噪聲，如圖5(d)放大區(qū)域所示。對(duì)比ARWTV算法和TpV算法的重建圖像，DO-TpV算法重建的圖像細(xì)節(jié)更為清晰，且圖像的噪聲水平也顯著降低，總體圖像質(zhì)量最佳。

4 結(jié) 論

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