基于投影域分解的多能譜CT造影劑物質(zhì)識別研究

孫英博， 孔慧華， 張雁霞

(1. 中北大學(xué) 理學(xué)院， 山西 太原 030051； 2. 中北大學(xué) 信息探測與處理山西省重點實驗室， 山西 太原 030051)

0 引 言

近年來， 隨著光子計數(shù)探測器和 CT 制造技術(shù)的發(fā)展[1-4]， 多能譜CT成像技術(shù)日益成熟， 在對小病灶的定性定量診斷、 減少高?；颊咴煊皠┦褂昧康确矫嬗型怀鲐暙I. 雙能 CT作為多能譜CT的一種， 不僅可以確定尿結(jié)石是否含有鈣， 還可以區(qū)分鈣、 骨與碘造影劑， 并提供一系列CT組織表征[5-9]. 但雙能CT僅限于兩種能量， 提供的能譜信息有限導(dǎo)致采集到的物質(zhì)衰減信息有限， 重建圖像中物質(zhì)區(qū)分度不高. 為了更充分地使用材料的多能量衰減特性， 增加能量通道數(shù)量以獲得更豐富的投影信息， 對CT成像技術(shù)有重要意義.

基于光子計數(shù)探測器的多能譜CT可以在一次掃描中得到物體在多個能區(qū)下的成像結(jié)果， 利用材料的能譜特性將材料間的差異最大化進而實現(xiàn)對材料的識別， A.P.Butler[10]等開發(fā)了一種基于Medipix2光子計數(shù)檢測器的多能譜CT(MARS)， 應(yīng)用光譜CT區(qū)分具有不同K-edge的造影劑； Daniel Y. Chong等[11]和J P Schlomka等[12]分別通過實驗對雙能CT中造影劑的可分離性及多能CT光子計數(shù)K-edge成像在臨床實驗中應(yīng)用的可行性進行了評估. 對于物質(zhì)識別重建算法的研究， 投影域分解具有重要的應(yīng)用價值， 現(xiàn)有文獻中已有多種投影分解方法. Alvarez和Macovski等[13]提出了基于雙能CT的基物質(zhì)分解模型和基效應(yīng)分解模型； 吳篤蕃[14]通過實驗針對雙能CT雙物質(zhì)分解及三物質(zhì)分解進行了研究； Bernhard Brendel[15]等在雙物質(zhì)分解模型的基礎(chǔ)上引入第三基礎(chǔ)成分,在將測量的投影數(shù)據(jù)分解成基礎(chǔ)分量投影之后， 執(zhí)行傳統(tǒng)的濾波反投影重建獲得基礎(chǔ)分量圖像， 實現(xiàn)了對造影劑物質(zhì)的識別.

本文對多能譜CT造影劑物質(zhì)識別進行研究， 在雙效應(yīng)分解模型的基礎(chǔ)上增加與待識別造影劑物質(zhì)K-edge相關(guān)的基函數(shù)， 應(yīng)用最小二乘法對不同能量通道下所獲得的待測物體的投影信息進行分解， 再通過傳統(tǒng)的反投影重建算法獲得康普頓效應(yīng)、 光電效應(yīng)及各造影劑物質(zhì)所對應(yīng)的重建圖像， 最后對各造影劑物質(zhì)對應(yīng)的重建圖像衰減值進行分析實現(xiàn)各造影劑物質(zhì)的分離成像、 定量表示及彩色表征.

1 基效應(yīng)分解模型

多能譜CT應(yīng)用光子計數(shù)探測器設(shè)置閾值范圍， 將衰減后的光子能量分能區(qū)進行計算， 得到物體在不同能量通道下的投影， 各能量通道中成像原理與傳統(tǒng)CT相同. 因此， 在多能譜 CT 中， 各能量通道下X 射線穿過物質(zhì)的衰減可以表達為該能量通道中不同能量下衰減的加權(quán)平均， 即

(1)

式中：EBin為當(dāng)前能量通道；I0和I分別為當(dāng)前能量通道下入射X射線的初始光子數(shù)及最終探測器接收到的光子數(shù)；S(E)是當(dāng)前能量通道下光子的等效能譜；μ(E,x)是物質(zhì)的線性衰減系數(shù)隨能量E和位置x的分布， 也就是成像需要求解的量； 積分路徑l表示這條射線在物體中所穿過的路徑.

由式(1)可得， 當(dāng)前能量通道EBin中， 投影函數(shù)為

(2)

直接求解μ(E,x)未知數(shù)太多， 通常需要對它加以一定的假設(shè)， 即假設(shè)物質(zhì)的衰減系數(shù)可以展開成I個基函數(shù)[12]

(3)

式中：fi(E)是給定的基函數(shù)；bi是組合系數(shù)[14].

常用的衰減系數(shù)模型有兩種[10]： 基物質(zhì)分解模型和雙效應(yīng)(光電效應(yīng)、 康普頓散射效應(yīng))分解模型. 常見的 “雙效應(yīng)分解”模型為

μ(E)≈apfph(E)+akfKN(E),

(4)

式中：fph(E)是光電效應(yīng)截面對能量的依賴關(guān)系.

(5)

式中：fKN(E)是 Klein-Nishina 函數(shù)， 表達了康普頓截面對能量的依賴關(guān)系.

(6)

式中：α=E/510.975 keV.

當(dāng)X射線的能量增大到原子K層電子的束縛能之上時， 射線可以與K層電子發(fā)生作用， 導(dǎo)致光電吸收截面突然增大產(chǎn)生K-edge效應(yīng). 原子序數(shù)越高， K 層電子的束縛越大， K-edge對應(yīng)能量也會越高. 當(dāng)待識別物質(zhì) K-edge 對應(yīng)能量足夠高時， 在衰減系數(shù)空間中產(chǎn)生一個新的基函數(shù)[8,14]， 此時雙效應(yīng)分解模型可進一步寫成

μ(E,x)≈apfph(E)+akfKN(E)+

(7)

式中：N是K-edge 物質(zhì)的種類，μn(E,x)是第n種物質(zhì)對應(yīng)的線性衰減系數(shù).

為方便計算， 記：ap=b1;ak=b2;an=bn+2，fph(E)=f1(E);fKN(E)=f2(E);μn(E,x)=fn+2(E)(n=1,2,…,N)， 將式(7)寫為式(3)的形式

(8)

2 基于最小二乘法的投影域分解

本文基于光子計數(shù)探測器， 對能譜CT投影的特性進行研究， 進而尋找投影分解的方法. 設(shè)所選M個能譜通道分別為EBin1,EBin2,EBinM， 各能量通道下的投影分別為P1,P2,…,PM， 具有K-edge的待識別造影劑物質(zhì)為N種， 實驗中通常取N≤3.

由式(2)， 式(8)可得各能量通道下的投影

(9)

式(9)可進一步寫成

(10)

(11)

根據(jù)式(5)， 式(6)及能譜信息， 將f1(EBinm),f2(EBinm)(m=1,2,…,M)表示為

m=1,2,…,M,

(12)

(13)

根據(jù)NIST數(shù)據(jù)庫中各造影劑物質(zhì)的線性衰減系數(shù)及能譜信息， 將fi(EBinm)(m=1,2,…,M;i=3,4,…,N+2)表示為

m=1,2,…,M;n=1,2,…,N;i=n+2,

(14)

式中：B(E)為第m個能量通道EBinm(m=1,2,…,M)中各能量下空掃描時探測器所接收到的光子數(shù)；CoefaN(E)為第m個能量通道EBinm(m=1,2,…,M)中各能量下第n種物質(zhì)所對應(yīng)的衰減系數(shù).

為方便求解， 將式(11)寫為針對投影矩陣單個像素的矩陣形式

P=AP*,

(15)

式中：P=[P1(i,j)P2(i,j) …PM(i,j)]T, (i,j)為投影的像素位置，i=1,2,…,I為采集投影的某一視角，I為采集投影的視角總數(shù)，j=1,2,…,J為采集投影的射線標(biāo)號，J為采集投影的射線總條數(shù).

式(14)有M個方程N+2個未知量個數(shù)， 為獲得更豐富的投影信息， 抑制重建中噪聲放大的影響， 通常M>N+2， 即式(14)為超定方程組.

應(yīng)用最小二乘法解超定方程組有如下定理：x0為超定方程組Ax=b最小二乘解的充分必要條件為x0是ATAx=ATb的解. 由此可得， 式(15)與式(16)同解

ATAP*=ATP.

(16)

對式(16)左右兩邊同乘(ATA)-1可得

P*=(ATA)-1ATP.

(17)

應(yīng)用FBP重建算法對各分解投影進行重建， 得光電效應(yīng)、 康普頓效應(yīng)及各具有K-edge的造影劑物質(zhì)所對應(yīng)的重建圖.

光電效應(yīng)和康普頓效應(yīng)重建圖顯示物體的內(nèi)部結(jié)構(gòu)， 各具有K-edge的造影劑物質(zhì)對應(yīng)的重建圖突出顯示各造影劑所在位置， 達到物質(zhì)識別與區(qū)分的目的.

3 仿真實驗

碘和釓具有足夠高的K-edge， 是醫(yī)療CT檢測中常用的造影劑物質(zhì). 因此， 仿真實驗中選用碘和釓元素作為待識別物質(zhì). 投影重建算法選用計算較為簡單， 運行時間相對較少的FBP重建算法.

3.1 數(shù)值仿真模型

仿真模體大小為10 cm×10 cm， 如圖 1 所示， 含有 8 個圓孔， 5 種材料， 各材料參數(shù)由表 1 列出. 3～8號圓孔外圍為厚度0.5 cm的骨物質(zhì). 所有材料的衰減系數(shù)， 密度等物理參數(shù)由 NIST 數(shù)據(jù)庫獲取.

圖 1 5種不同材料組成的圓孔模體截面圖Fig.1 Cross-sectional view of a circular hole phantom composed of 5 different materials

3.2 數(shù)值仿真環(huán)境及能譜通道范圍

仿真電壓采用120 kVp， 以物質(zhì)K-edge所在能量： 碘 (33 keV)、 釓(50.2 keV)為通道劃分界限的參考值， 將能譜區(qū)域劃分為 6 個能量通道， 在圖 2 中用不同顏色表示， 分別為：

EBin1={25～30 keV}；EBin2={30～40 keV}；

EBin3={40～50 keV}；EBin4={50～60 keV}；

EBin5={60～80 keV}；EBin6={80～120 keV}.

仿真實驗中采用等距扇束掃描， 模體分辨率為 256×256像素. 探測器為由320 個探元組成的光子技術(shù)探測器， 采樣角度為360°， 采樣間隔為1°.

圖 2 用于數(shù)值模擬源光子發(fā)射光譜Fig.2 For numerical simulation of source photon emission spectroscopy

3.3 數(shù)值仿真實驗與結(jié)果

根據(jù)式(17)將6個能量通道下所采集的投影進行分解后所得投影見圖 3， 圖 3(a)～(d)分別對應(yīng)光電效應(yīng)、 康普頓效應(yīng)、 碘圖、 釓圖的分解投影.

應(yīng)用FBP算法對圖3中分解后所得投影進行重建， 重建圖見圖 4， 圖 4(a)～(d)分別對應(yīng)光電效應(yīng)、 康普頓效應(yīng)、 碘圖、 釓圖的重建結(jié)果.

圖 3 6個能量通道下投影分解后所得投影Fig.3 Projection obtained after projection decomposition under 6 energy channels

圖 4 分解后投影重建結(jié)果Fig.4 Projection reconstruction result after decomposition

由圖 4(c) 可見， 3, 4, 5號含釓元素的圓孔較其他位置更亮， 且亮度隨釓濃度的降低逐漸降低， 其中3～8號圓孔位置圖像衰減值分別為： 0.005 3, 0.004 2, 0.003 6, 0.003 0, 0.003 0, 0.003 0. 由此可見， 釓圖中含不同濃度碘元素的圓孔被作為衰減值相同的“背景部分”， 含釓元素的圓孔灰度值隨釓濃度的降低而降低且其與“背景部分”灰度值的差與釓濃度成正比例 ， 見表 2.

表 2 圖4(c)圖像灰度值與釓濃度關(guān)系Tab.2 The relationship between image gray value and gadolinium concentration in Fig.4(c)

由圖 4(d) 可見， 6, 7, 8號含碘元素的圓孔較其他位置更亮， 且亮度隨碘濃度的升高逐漸升高， 其中3～8號圓孔位置圖像灰度值分別為： 0.004 0, 0.004 0, 0.004 0, 0.007 1, 0.005 1, 0.005 0. 由此可見， 碘圖中含不同濃度釓元素的圓孔被作為灰度值相同的“背景部分”， 含碘元素的圓孔灰度值隨碘濃度的升高而升高且其與“背景部分”灰度值的差與碘濃度成正比例 ， 見表 3.

表 3 圖4(d)圖像灰度值與碘濃度關(guān)系Tab.3 The relationship between image gray value and iodine concentration in Fig.4(d)

綜上， 實驗通過對各材料對應(yīng)位置灰度值的對比分析實現(xiàn)對造影劑物質(zhì)濃度的定量分析.

3.4 待識別物質(zhì)的分離成像及彩色表征

為了更方便清晰地識別仿真模體中的造影劑物質(zhì)， 將圖 4(c), 圖4(d)的窗口顯示閾值分別調(diào)整為[0.003 1,0.005 3]和[0.004 1,0.007 1]實現(xiàn)釓和碘的分離成像， 然后應(yīng)用RGB顏色通道對分離出的釓和碘進行彩色表征， 如圖 5(a)， 圖 5(b)， 最后將釓和碘的彩色表征圖像與顯示仿真模體內(nèi)部結(jié)構(gòu)的圖 4(a), 圖 4(b)進行加權(quán)融合. 圖 5(d)既顯示了仿真模體的內(nèi)部結(jié)構(gòu)， 又實現(xiàn)了對釓和碘元素的識別及顏色表征， 且顏色越深表示濃度越高.

圖 5 碘和釓的分離成像與彩色表征Fig.5 Separation imaging and color characterization of iodine and strontium

4 結(jié) 論

醫(yī)療CT成像中造影劑物質(zhì)的識別及定量分析有助于研究物體中造影劑物質(zhì)的分布及其累積程度， 對病灶診斷、 降低輻射劑量等問題的研究有重要意義. 目前基于前處理的投影域分解方法可通過基物質(zhì)分解或基效應(yīng)分解實現(xiàn)對兩種或三種造影劑物質(zhì)的識別. 本文基于光子計數(shù)探測器， 推導(dǎo)了能譜CT投影基效應(yīng)分解模型， 應(yīng)用最小二乘法對多個能譜通道下的投影進行分解和重建， 達到了造影劑物質(zhì)識別與定量分析的目的. 由于光子計數(shù)探測器可以同時獲取多個能量通道的投影， 從而可以實現(xiàn)多個基物質(zhì)的分解， 為能譜CT物質(zhì)識別提供有效支持.