孫英博， 孔慧華， 張雁霞

(1. 中北大學(xué) 理學(xué)院， 山西 太原 030051； 2. 中北大學(xué) 信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室， 山西 太原 030051)

0 引 言

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

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

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

(4)

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

(5)

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

(6)

式中：α=E/510.975 keV.

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

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

(7)

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

為方便計(jì)算， 記：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)寫(xiě)為式(3)的形式

(8)

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

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

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

(9)

式(9)可進(jìn)一步寫(xiě)成

(10)

(11)

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

m=1,2,…,M,

(12)

(13)

根據(jù)NIST數(shù)據(jù)庫(kù)中各造影劑物質(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個(gè)能量通道EBinm(m=1,2,…,M)中各能量下空掃描時(shí)探測(cè)器所接收到的光子數(shù)；CoefaN(E)為第m個(gè)能量通道EBinm(m=1,2,…,M)中各能量下第n種物質(zhì)所對(duì)應(yīng)的衰減系數(shù).

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

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)號(hào)，J為采集投影的射線總條數(shù).

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

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

ATAP*=ATP.

(16)

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

P*=(ATA)-1ATP.

(17)

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

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

3 仿真實(shí)驗(yàn)

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

3.1 數(shù)值仿真模型

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

圖 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 個(gè)能量通道， 在圖 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}.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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