牛善洲，梁禮境，李 碩，張夢(mèng)真，邱 洋，劉漢明，李 楠

1.贛南師范大學(xué)數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，江西贛州341000

2.贛南師范大學(xué)贛州市計(jì)算成像重點(diǎn)實(shí)驗(yàn)室，江西贛州341000

3.贛南師范大學(xué)經(jīng)濟(jì)管理學(xué)院，江西贛州341000

X 射線計(jì)算機(jī)斷層（computed tomography，CT）成像在疾病防治與診斷方面取得了巨大成就，其作為一種無(wú)痛、無(wú)創(chuàng)傷性的檢查方法，可用于全身多臟器檢查，是現(xiàn)代影像學(xué)的杰出代表[1]。CT 圖像的質(zhì)量與X 射線的輻射劑量緊密相關(guān)，劑量越高圖像質(zhì)量越好，但是過(guò)高劑量的X 射線照射會(huì)造成患者脫發(fā)、皮膚灼傷，甚至可能誘發(fā)癌癥或者遺傳性疾病[2-4]。針對(duì)CT檢查過(guò)程中的X 射線照射劑量問(wèn)題，世界衛(wèi)生組織、國(guó)際放射委員會(huì)以及國(guó)際醫(yī)學(xué)物理組織制定了X 射線照射劑量保證和劑量控制標(biāo)準(zhǔn)，極力主張X 線CT 檢查應(yīng)遵循實(shí)踐正當(dāng)性、防護(hù)最優(yōu)化的原則，希望以最小的代價(jià)和劑量獲取最好的CT 影像診斷效果。因此，在保證圖像質(zhì)量的前提下最大限度地減少X射線的輻射劑量已成為CT成像領(lǐng)域亟待解決的關(guān)鍵問(wèn)題。

降低管電流或者管電壓是實(shí)現(xiàn)低劑量CT成像最簡(jiǎn)單有效的方法，但是投影數(shù)據(jù)會(huì)受到量子噪聲和電子噪聲的污染，從而使得FBP 算法重建圖像的質(zhì)量嚴(yán)重退化，無(wú)法滿足臨床診斷的需求[5-6]。低劑量CT 圖像重建方法主要分為基于圖像域重建和基于投影域重建兩大類?；趫D像域的重建方法是在圖像域中對(duì)待重建的圖像進(jìn)行系統(tǒng)建模，然后通過(guò)迭代算法求解目標(biāo)函數(shù)進(jìn)行圖像重建[7-8]。但是，由于迭代重建需反復(fù)進(jìn)行投影與反投影運(yùn)算，導(dǎo)致運(yùn)算量大，重建速度慢，無(wú)法滿足快速診斷的需要。基于投影域的方法利用投影數(shù)據(jù)的噪聲統(tǒng)計(jì)特性建立投影數(shù)據(jù)恢復(fù)模型，得到高質(zhì)量的投影數(shù)據(jù)后利用FBP 算法重建最終CT 圖像[9-15]。該類方法與圖像域迭代重建方法相比計(jì)算時(shí)間大幅減少，在去除噪聲的同時(shí)可以很好地保持圖像的空間分辨率。Lu 等[16]通過(guò)大量的臨床低劑量CT 實(shí)驗(yàn)以及相應(yīng)的理論分析，提出經(jīng)對(duì)數(shù)變換后的投影數(shù)據(jù)噪聲近似滿足均值和方差呈非線性關(guān)系的高斯分布?；谕队皵?shù)據(jù)這一重要統(tǒng)計(jì)特性，Wang等[12-13，17]提出多種基于吉波斯隨機(jī)場(chǎng)先驗(yàn)的懲罰加權(quán)最小二乘（penalized weighted least-squares，PWLS）低劑量CT 投影數(shù)據(jù)恢復(fù)方法，如基于KL 變換的懲罰加權(quán)最小二乘方法[13]、多尺度加權(quán)懲罰加權(quán)最小二乘方法[12]等。但Wang等提出的基于吉波斯隨機(jī)場(chǎng)先驗(yàn)的懲罰加權(quán)最小二乘方法的目標(biāo)函數(shù)使用的是二次懲罰函數(shù)，而在投影數(shù)據(jù)恢復(fù)過(guò)程中，二次懲罰函數(shù)對(duì)投影數(shù)據(jù)的均勻區(qū)域和邊緣區(qū)域使用相同的懲罰權(quán)重會(huì)造成圖像的邊緣和結(jié)構(gòu)細(xì)節(jié)信息的丟失，從而會(huì)導(dǎo)致重建圖像的空間分辨率下降。

圖像塊流形的一個(gè)重要特性是其具有低維結(jié)構(gòu)，并且包含豐富的結(jié)構(gòu)信息。圖像塊可以視為嵌在高維空間中低維流形的點(diǎn)云，極小化圖像塊流形的維數(shù)可以作為先驗(yàn)信息進(jìn)行圖像恢復(fù)。低維流形先驗(yàn)將線性對(duì)象的低秩先驗(yàn)擴(kuò)展到具有更復(fù)雜結(jié)構(gòu)的數(shù)據(jù)對(duì)象上，其可以更好地恢復(fù)圖像的結(jié)構(gòu)細(xì)節(jié)信息。因此，將低維流形的維數(shù)作為先驗(yàn)信息引入到低劑量CT 投影數(shù)據(jù)恢復(fù)過(guò)程中，提出了一個(gè)基于低維流形懲罰加權(quán)最小二乘（penalized weighted least-squares based on low-dimensional manifold，PWLS-LDMM）的低劑量CT 投影數(shù)據(jù)恢復(fù)模型，通過(guò)極小化相應(yīng)的目標(biāo)函數(shù)可以得到理想的投影數(shù)據(jù)，再使用FBP 算法對(duì)恢復(fù)后的投影數(shù)據(jù)進(jìn)行CT 重建。實(shí)驗(yàn)結(jié)果表明該方法在有效抑制低劑量CT圖像噪聲和偽影的同時(shí)可以很好地保持圖像的邊緣和結(jié)構(gòu)細(xì)節(jié)特征。

1 模型與方法

1.1 PWLS重建模型

先前的研究工作表明[18]，在低劑量CT 掃描協(xié)議下經(jīng)過(guò)系統(tǒng)校準(zhǔn)和對(duì)數(shù)轉(zhuǎn)換的投影數(shù)據(jù)近似遵循高斯分布，數(shù)據(jù)樣本均值和方差之間滿足如下的公式：

其中，I0是入射X線強(qiáng)度，是第i個(gè)探測(cè)器單元上投影數(shù)據(jù)（對(duì)數(shù)變換后）的均值，是電子噪聲的方差。設(shè)y=(y1,y2,…,yM)T為系統(tǒng)校正和對(duì)數(shù)變換后的投影數(shù)據(jù)測(cè)量值，q=(q1,q2,…,qM)T為待估計(jì)投影數(shù)據(jù)的理想值，則基于PWLS的低劑量CT重建模型為：

1.2 低維流形先驗(yàn)

設(shè)投影數(shù)據(jù)q為m×n的矩陣，選取q中的任意一個(gè)像素點(diǎn)(i,j)，在像素點(diǎn)(i,j)處選取大小為s1×s2的圖像塊Ψx(q)，其中，s1和s2分別為圖像塊Ψx(q)的行數(shù)和列數(shù)，為s1×s2矩陣左上角的像素。對(duì)q中所有像素點(diǎn)的圖像塊Ψx(q)組合成q的點(diǎn)云Ψ(q)：

其中，d=s1×s2,Rd為d維的歐幾里德空間。

基于微分流形理論，點(diǎn)云Ψ(q)可以模擬為一個(gè)嵌在Rd空間中的低維光滑流形M=∪l Ml，其中Ml是對(duì)應(yīng)q中不同區(qū)域的流形，M為q的圖像塊流形。投影數(shù)據(jù)q的低維流形先驗(yàn)可以寫為[19]：

其中，dim(M)為流形M的維數(shù)，αi(w)=wi為指標(biāo)函數(shù)，且w=(w1,w2,…,wd)∈M,?Mαi(w)為指標(biāo)函數(shù)αi(w)在流形M上的梯度且為歐幾里德范數(shù)。

1.3 PWLS-LDMM重建模型

PWLS-LDMM的投影數(shù)據(jù)恢復(fù)模型可以寫為：

其中，p和q都為待恢復(fù)投影數(shù)據(jù)，T 為轉(zhuǎn)置運(yùn)算，Σ為對(duì)角矩陣，對(duì)角矩陣對(duì)角線上的元素為方差β2>0 為正則化參數(shù)為q的低維流形先驗(yàn)。

使用如下的交替優(yōu)化算法求解PWLS-LDMM投影數(shù)據(jù)恢復(fù)模型：

其中，k為迭代次數(shù)。

由于式（6）的目標(biāo)函數(shù)為光滑的二次凸函數(shù)，可以對(duì)其求導(dǎo)并令導(dǎo)數(shù)為零來(lái)進(jìn)行求解。得到其解為：

其中，I為單位矩陣，Σ-1為對(duì)角矩陣Σ的逆矩陣，β1為正則化參數(shù)，(Σ-1+β1I)為對(duì)角矩陣。

使用交替優(yōu)化法求解式（7），先固定流形更新投影數(shù)據(jù)，再固定投影數(shù)據(jù)更新流形。具體求解步驟如下：

（1）先固定流形初始值M(k)，通過(guò)求解下式得到指標(biāo)函數(shù)α(k+1)i(i=1,2,…,d)和投影數(shù)據(jù)q(k+1)：

其中，Ψ*為Ψ的伴隨算子，Ψ*Ψ為對(duì)角矩陣。

（2）通過(guò)指標(biāo)函數(shù)對(duì)流形進(jìn)行更新，公式如下：交替求解公式（6）、（7）得到恢復(fù)后的投影數(shù)據(jù)后，再使用FBP算法，即可重建出最終CT圖像。

綜上所述，基于低維流形先驗(yàn)加權(quán)最小二乘的低劑量CT重建方法的計(jì)算步驟可以寫成如下所示：

步驟1 初始化令k=0；

步驟2 將qk代入式（6），通過(guò)式（8）得到pk+1；

步驟3 將pk+1代入式（7），通過(guò)式（9）~（13）得到αk+1和qk+1，并通過(guò)式（14）對(duì)流形進(jìn)行更新；

步驟4 重復(fù)步驟2~步驟3，直至滿足迭代終止條件；

步驟5 利用FBP算法將步驟3得到的投影數(shù)據(jù)qk+1重建出CT圖像。

1.4 對(duì)比方法

為了驗(yàn)證PWLS-LDMM 方法的有效性，將其與傳統(tǒng)的FBP方法，基于二次模先驗(yàn)的懲罰加權(quán)最小二乘方法[14]（penalized weighted least-squares via quadratic membrane，PWLS-QM）和基于字典學(xué)習(xí)先驗(yàn)的懲罰加權(quán)最小二乘（penalized weighted least-squares via dictionary learning，PWLS-DL）方法[20]進(jìn)行比較。

PWLS-QM方法的目標(biāo)函數(shù)為：

其中，Ni是投影數(shù)據(jù)中第i個(gè)像素的一階鄰域，wim在水平方向取值為1，在垂直方向取值為0.25。

PWLS-DL方法的目標(biāo)函數(shù)為：

其中，D∈RL×K為訓(xùn)練字典，αs∈RL×1是含有少量非零元的向量，Es∈RL×K是從圖像f提取的塊矩陣，vs是拉格朗日乘子，γ是超參數(shù)。

1.5 實(shí)驗(yàn)設(shè)置

分別使用Shepp-Logan 體膜（如圖1（a）所示）和XCAT 體膜（如圖1（b）所示）以及臨床數(shù)據(jù)實(shí)驗(yàn)來(lái)對(duì)所提出的PWLS-LDMM方法進(jìn)行驗(yàn)證，并使用文獻(xiàn)[21]中的方法仿真生成低劑量CT 投影數(shù)據(jù)。數(shù)值體膜實(shí)驗(yàn)中，CT成像幾何采用扇形束和弧形探測(cè)器，X射線源到旋轉(zhuǎn)中心和探測(cè)器的距離分別為570 mm 和1 040 mm，旋轉(zhuǎn)角在[]0,2π 間采樣值為1 160，每個(gè)采樣角對(duì)應(yīng)672個(gè)探測(cè)器單元，探測(cè)器單元的大小為1.407 mm。圖像維數(shù)為512×512，像素大小為1 mm×1 mm。Shepp-Logan和XCAT 體膜室驗(yàn)的入射光子強(qiáng)度I0分別為5.0×104、3.0×104，電子噪聲的方差σ2e都為10.0。在臨床數(shù)據(jù)實(shí)驗(yàn)中的腹部投影數(shù)據(jù)，由商業(yè)單排CT 設(shè)備采集。X 射線源到旋轉(zhuǎn)中心和探測(cè)器的距離分別為570 mm 和1 040 mm，旋轉(zhuǎn)角在間采樣值為1 160，每個(gè)采樣角對(duì)應(yīng)672個(gè)探測(cè)器單元，探測(cè)器單元的大小為1.407 mm。在球管電壓為120 kV 下，分別獲取400 mAs（高劑量）和50 mAs（低劑量）的CT 投影數(shù)據(jù)。重建圖像的維數(shù)為512×512，像素大小為0.6 mm×0.6 mm。分別使用相對(duì)均方根誤差（relative root mean square error，RRMSE），結(jié)構(gòu)相似性指標(biāo)[22]（structural similarity，SSIM）和特征相似性指標(biāo)[23]（feature similarity，F(xiàn)SIM）對(duì)重建圖像進(jìn)行定量分析。使用MatlabR2020b（Ubuntu）編程，所有程序都在Intel?Reon?Gold 6234 2.33 GHz，16 核處理器，64 GB內(nèi)存PC機(jī)上實(shí)現(xiàn)。

圖1 數(shù)值體膜圖像Fig.1 Digital phantom

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

2.1 Shepp-Logan體膜實(shí)驗(yàn)

圖2為不同方法重建的Shepp-Logan 體膜圖像。圖2（a）為FBP算法重建的圖像；圖2（b）為PWLS-QM方法重建的圖像；圖2（c）為PWLS-DL 方法重建的圖像；圖2（d）為PWLS-LDMM方法重建的圖像。從圖2可以看出，F(xiàn)BP 方法重建的圖像質(zhì)量嚴(yán)重退化，存在大量的噪聲和條形偽影。PWLS-GS重建的圖像相比于FBP重建圖像、噪聲和偽影得到抑制，但仍含有一些噪聲和條形偽影。PWLS-DL方法重建的圖像噪聲得到了一定程度上的抑制，但仍含有條形偽影。PWLS-LDMM 方法重建的圖像中噪聲和條形偽影都被有效抑制，同時(shí)可以保持邊緣和細(xì)節(jié)結(jié)構(gòu)特征。因此，本文提出的PWLSLDMM方法在抑制噪聲偽影和保持圖像分辨率方面都有良好表現(xiàn)，可以保持良好的圖像一致性。

圖2 不同方法重建的Shepp-Logan體膜圖像Fig.2 Shepp-Logan images reconstructed by different methods

表1為Shepp-Logan 體膜重建圖像的結(jié)構(gòu)相似性、特征相似性和相對(duì)均方根誤差。與FPB、PWLS和PWLSDL 方法相比，PWLS-LDMM 方法的相對(duì)均方根誤差分別降低了64.87%、54.81%和7.02%。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了16.78%、1.88%和1.91%；FSIM 值分別提高了12.68%、2.73%和1.65%。PWLS-LDMM 方法重建的圖像具有最高的SSIM 值和FSIM 值，即最接近于真實(shí)的體膜圖像。

表1 Shepp-Logan體膜圖像的RRMSE、SSIM和FSIMTable 1 RRMSE，SSIM and FSIM of Shepp-Logan

為了更清晰地比較各種方法對(duì)噪聲和偽影的抑制效果，圖3給出了圖1中感興趣區(qū)域（如圖1（a）中紅色方框所示）的局部放大圖。與FBP、PWLS-QM 和PWLSDL 三種方法相比，PWLS-LDMM 方法能在有效去除噪聲和偽影的同時(shí)可以保持圖像的邊緣信息和結(jié)構(gòu)細(xì)節(jié)特征。進(jìn)一步，使用SSIM 和FSIM 值對(duì)局部放大圖進(jìn)行定量分析。圖4 給出了圖3 中各重建結(jié)果的SSIM 值和FSIM 值。與FBP、PWLS-QM、PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了329.71%、11.74%和23.41%；FSIM 值分別提高了111.84%、8.31%和13.48%。PWLS-LDMM方法的SSIM和FSIM值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實(shí)體膜。

圖3 Shepp-Logan體膜的局部放大圖Fig.3 Zoomed-in views of Shepp-Logan phantom

圖4 圖3中圖像的SSIM值和FSIM值Fig.4 SSIM and FSIM values of images in Figure 3

2.2 XCAT體膜實(shí)驗(yàn)

圖5為不同方法重建的XCAT 體膜圖像。圖5（a）為FBP方法重建的圖像；圖5（b）為PWLS-QM方法重建的圖像；圖5（c）為PWLS-DL 方法重建的圖像；圖5（d）為PWLS-LDMM 方法重建的圖像。從圖5 可以看出，F(xiàn)BP 方法重建的圖像存在較多的噪聲和條形偽影；PWLS-GS 重建的圖像中仍含有一些噪聲和條形偽影；PWLS-DL方法重建的圖像噪聲得到了一定程度上的抑制，但仍含有條形偽影；PWLS-LDMM方法重建的圖像中噪聲和條形偽影都被有效抑制。因此，本文提出的PWLS-LDMM方法在抑制噪聲偽影和保持圖像分辨率方面都有良好表現(xiàn)，對(duì)比其他三種方法得到了更優(yōu)質(zhì)的重建圖像。

圖5 不同方法重建的XCAT體膜圖像Fig.5 XCAT images reconstructed by different methods

表2為XCAT 體膜重建圖像的結(jié)構(gòu)相似性、特征相似性和相對(duì)均方根誤差。如表2所示，PWLS-LDMM方法與其他方法相比，在降低相對(duì)均方根誤差方面有很好的表現(xiàn)。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLS-LDMM方法的相對(duì)均方根誤差分別降低了37.46%、22.17%和11.48%。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLS-LDMM方法的SSIM值分別提高了3.33%、0.73%和1.22%；FSIM 值分別提高了0.83%、0.68%和0.71%。PWLS-LDMM 方法的SSIM 和FSIM 值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實(shí)體膜。

表2 XCAT體膜圖像的RRMSE、SSIM和FSIM Table 2 RRMSE，SSIM and FSIM of XCAT images

同樣，為了更清晰地比較各種方法對(duì)噪聲和偽影的抑制效果，圖6給出了圖1中感興趣區(qū)域（如圖1（b）中大小為60×100 的紅色方框所示）的局部放大圖，并使用SSIM 和FSIM 值來(lái)對(duì)結(jié)果的感興趣區(qū)域進(jìn)行定量分析。圖7 給出了圖6 中各重建結(jié)果的SSIM 值和FSIM值。與FBP、PWLS-QM 和PWLS-DL 方法相比，PWLSLDMM 方法的SSIM 值分別提高了80.17%、8.97%和20.46%；FSIM 值分別提高了26.17%、4.98%和8.30%。PWLS-LDMM方法的SSIM和FSIM值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM 方法重建的圖像最接近真實(shí)體膜。

圖6 XCAT體膜圖像的局部放大Fig.6 Zoomed-in views of XCAT phantom

圖7 圖6中圖像的SSIM值和FSIM值Fig.7 SSIM and FSIM values of images in Figure 6

2.3 臨床數(shù)據(jù)實(shí)驗(yàn)

圖8為臨床腹部CT投影數(shù)據(jù)重建圖像。圖8（a）為400 mAs條件下FBP算法重建得到的高劑量圖像，其作為金標(biāo)準(zhǔn)圖像來(lái)對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行定量分析；圖8（b）為50 mAs 條件下FBP 方法重建的圖像；圖8（c）為50 mAs條件下PWLS-QM 方法重建的圖像；圖8（d）為50 mAs條件下PWLS-DL方法重建的圖像；圖8（e）為50 mAs條件下PWLS-LDMM 方法重建的圖像。從圖中可以看出，本文提出方法在噪聲濾除和邊緣細(xì)節(jié)以及圖像分辨率保持方面表現(xiàn)更為突出，優(yōu)于對(duì)比方法。

圖8 臨床數(shù)據(jù)重建圖像Fig.8 Images reconstructed by clinical data

表3為低劑量臨床數(shù)據(jù)重建的CT圖像的結(jié)構(gòu)相似性、特征相似性和相對(duì)均方根誤差。如表3所示，PWLSLDMM 方法與其他方法相比，在降低相對(duì)均方根誤差方面有很好的表現(xiàn)。與FBP、PWLS-QM 和PWLS-DL方法相比，PWLS-LDMM 方法的相對(duì)均方根誤差分別降低了45.96%、25.61%和15.87%。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了19.12%、7.46%和8.63%；FSIM 值分別提高了4.37%、1.83%和1.14%。PWLS-LDMM 方法的SSIM和FSIM 值均高于其他方法的SSIM 值和FSIM 值，即PWLS-LDMM方法重建的圖像最接近高劑量圖像。

表3 低劑量臨床數(shù)據(jù)重建的CT圖像的RRMSE、SSIM和FSIMTable 3 RRMSE，SSIM and FSIM of low-dose CT images reconstructed by clinical data

進(jìn)一步，為了更清晰地比較各種方法對(duì)噪聲和偽影的抑制效果，圖9給出了重建圖像感興趣區(qū)域（如圖8（a）中紅色方框所示）的局部放大圖。圖10給出了圖9中各重建結(jié)果的SSIM值和FSIM值。與FBP、PWLS-QM和PWLS-DL 方法相比，PWLS-LDMM 方法的SSIM 值分別提高了46.88%、1.02%和26.71%；FSIM值分別提高了8.34%、2.59%和3.78%。

圖9 臨床數(shù)據(jù)重建CT圖像的局部放大圖Fig.9 Zoomed-in views of CT images reconstructed by clinical data

圖10 圖9中圖像的SSIM值和FSIM值Fig.10 SSIM and FSIM values of images in Figure 9

3 結(jié)論

為了提高低劑量CT 重建圖像的質(zhì)量，將低維流形作為先驗(yàn)信息引入投影數(shù)據(jù)的恢復(fù)過(guò)程中，提出了一種基于低維流形先驗(yàn)的低劑量CT重建方法。數(shù)值仿真和臨床數(shù)據(jù)實(shí)驗(yàn)結(jié)果表明，本文提出的低劑量CT 投影數(shù)據(jù)恢復(fù)模型在噪聲和條形偽影抑制方面都取得了很好效果。為了驗(yàn)證PWLS-LDMM 方法的有效性，使用多種定量指標(biāo)對(duì)重建結(jié)果進(jìn)行分析與比較。與FBP、PWLS-QM 和PWLS-DL 對(duì)比方法相比，PWLS-LDMM方法在相對(duì)均方根誤差、結(jié)構(gòu)相似性和特征相似性等方面均有上佳的表現(xiàn)。PWLS-LDMM方法在噪聲去除，條形偽影抑制已經(jīng)結(jié)構(gòu)細(xì)節(jié)保持方面明顯優(yōu)于對(duì)比方法。