【作 者】 孫智鵬,李明,馬健,馬瑾瑾,梁國棟
東軟醫(yī)療系統(tǒng)股份有限公司,沈陽市,110167
PET數(shù)據(jù)采集過程中通過時(shí)間窗和能量窗對(duì)探測(cè)到的γ光子對(duì)進(jìn)行篩選,滿足篩選條件的稱為符合事件。符合事件分為隨機(jī)符合、散射符合和真符合。散射符合和隨機(jī)符合都是噪聲數(shù)據(jù),需要估計(jì)其數(shù)量和分布,從而對(duì)數(shù)據(jù)進(jìn)行校正。
常見的散射校正方法分為四大類:能窗法[1]、卷積法[2]、物理模型法和AI模型法。前兩種方法出現(xiàn)的時(shí)間比較早,且都存在比較多的近似。
能窗法通過掃描一個(gè)特定放射源,估算不同空間位置上不同能窗內(nèi)散射符合和真符合的比例,從而估算出每一條LOR上真符合的計(jì)數(shù)。但是實(shí)際上不同位置上的散射與真符合的比例,與放射源和物質(zhì)的空間分布有關(guān),理論上只有在被掃描物體的放射源分布、物質(zhì)空間分布、擺位等屬性與前面提到的特定放射源的分布等完全一致時(shí),能窗法才是準(zhǔn)確的。臨床條件下被掃描物體更加復(fù)雜,難以確保散射估計(jì)的準(zhǔn)確性。
卷積法假設(shè)空間中每一個(gè)點(diǎn)的散射分布都呈高斯分布,所以真符合與高斯核函數(shù)的卷積被認(rèn)為是散射分布。但實(shí)際上散射的空間分布比較復(fù)雜,和物體的密度分布、形狀關(guān)系較大,所以卷積法也不夠準(zhǔn)確。
基于上面提到的問題,以單散射模擬[3]為代表的模型法被提出,通過建立物理模型和統(tǒng)計(jì)模型,單散射模擬可以較準(zhǔn)確地計(jì)算出不同源分布、不同物質(zhì)分布下的單散射分布。但實(shí)際上臨床條件下約有15%的散射是二次散射或多散射,所以單散射模擬被推廣到了雙散射模擬。雙散射模擬雖然可以估算出雙散射,但是解決不了視野外散射、多散射問題。
此外,還有一些基于人工神經(jīng)網(wǎng)絡(luò)的方法。這類方法[4]一般基于卷積神經(jīng)網(wǎng)絡(luò)(或者擴(kuò)展網(wǎng)絡(luò),如U-Net等)。使用核素和物質(zhì)的空間分布信息作為神經(jīng)網(wǎng)絡(luò)的輸入,前者可以是PET符合數(shù)據(jù),也可以是PET未經(jīng)散射校正的圖像;后者指的是衰減分布圖或者是該圖的投影。網(wǎng)絡(luò)的輸出是物理模型法計(jì)算的結(jié)果。通過迭代調(diào)整神經(jīng)網(wǎng)絡(luò)的參數(shù),使得網(wǎng)絡(luò)的輸出結(jié)果與預(yù)設(shè)的結(jié)果越來越接近,直到滿足收斂條件。但是因?yàn)樯窠?jīng)網(wǎng)絡(luò)本身難以解釋,且存在訓(xùn)練樣本不能保證覆蓋所有可能情形,導(dǎo)致某些情況下結(jié)果不可預(yù)知的問題,所以目前此類方法尚未用于臨床,故暫不介紹。
下面按照技術(shù)發(fā)展的時(shí)間順序,分別介紹單散射模擬、視野外散射處理方法以及完全蒙特卡羅計(jì)算法,并進(jìn)行效果對(duì)比。
單散射模擬(single scatter simulation,SSS)是一種結(jié)合了物理模型和統(tǒng)計(jì)模型的方法[5-6]。如圖1所示,一個(gè)正電子湮滅事件在橢圓柱中發(fā)生,生成了一對(duì)γ射線。其中向下的一條在S點(diǎn)處發(fā)生了康普頓散射,損失一部分能量,行進(jìn)方向也發(fā)生了改變,形成了一個(gè)散射事件。
圖1 散射事件示意Fig.1 Schematic diagram of scattering events
晶體對(duì)AB的散射可以用式(1)~式(3)來表示:
式中:VS是某個(gè)散射點(diǎn)S對(duì)應(yīng)的體素;是晶體A和晶體B的立體角,表示空間幾何上的探測(cè)概率;μ是線性衰減系數(shù);σC是康普頓截面,兩者之商是單位體積內(nèi)的原子數(shù)(其物理單位是atom/cm3),與發(fā)生康普頓作用的概率成正比;是康普頓差分截面,表示發(fā)生當(dāng)前散射角的概率;是探測(cè)器效率,是發(fā)生散射,光子能量降低后的效率;s為光子傳播路徑。
SSS的總體思路是,先確定所有可能的散射點(diǎn),計(jì)算每個(gè)散射點(diǎn)與任意兩個(gè)晶體形成折線上的散射值,遍歷所有散射點(diǎn)和晶體對(duì)的組合后,計(jì)算出每個(gè)晶體對(duì)上的散射值。因?yàn)镾SS算法計(jì)算出的是與真實(shí)散射空間等比例的散射分布,所以通常會(huì)通過尾部擬合的方式對(duì)計(jì)算結(jié)果進(jìn)行縮放,得到與實(shí)際散射匹配的結(jié)果。
在實(shí)際使用中,尾部擬合有時(shí)會(huì)不穩(wěn)定,特別是當(dāng)被掃描目標(biāo)體積較大時(shí),留給擬合用的數(shù)據(jù)量變少、噪聲增加,容易出現(xiàn)欠估計(jì)或過估計(jì)的情況。圖2中展示了數(shù)據(jù)中的過估計(jì)和欠估計(jì)情況。其中,實(shí)線是真符合和散射的分布圖,點(diǎn)狀波動(dòng)較大的虛線是真實(shí)的散射分布。三角形、點(diǎn)劃線和星號(hào)分別是散射過估計(jì)、準(zhǔn)確估計(jì)和欠估計(jì)這3種情況。散射估計(jì)誤差對(duì)圖像的影響,如圖3所示。其中散射過估計(jì)會(huì)使得一部分真符合被丟棄,造成圖像中間部分偏低(見圖3(c)):較低虛線是背景部分的理論值,較高虛線是熱區(qū)理論值。圖像中心部分的背景部分偏低,圖像對(duì)比度過高;而散射欠估計(jì)(見圖3(a))則會(huì)導(dǎo)致圖像中心的背景部分偏高,圖像對(duì)比度降低。
圖2 散射過估計(jì)和欠估計(jì)Fig.2 Scatter over-fitting and under-fitting
圖3 散射估計(jì)誤差對(duì)圖像的影響Fig.3 Effects in PET images caused by scatter estimation error
為了防止這樣的情況發(fā)生,一種基于蒙特卡羅的散射分?jǐn)?shù)計(jì)算方法被引入SSS擬合過程中[7]。該方法通過建立探測(cè)器和核素、物質(zhì)分布的模型,通過光子傳播模型,只需要比較小的計(jì)算量,就可以準(zhǔn)確地計(jì)算出數(shù)據(jù)中散射事件的比例。這樣擬合的結(jié)果就可以被這個(gè)比例系數(shù)約束,輸出的散射分布估計(jì)結(jié)果更準(zhǔn)確、穩(wěn)定。
蒙特卡羅方法通過對(duì)探測(cè)器、放射性核素和衰減系數(shù)圖進(jìn)行建?;蛳袼鼗?,結(jié)合采樣技術(shù)和物理模型,模擬出正電子湮滅和γ光子傳播過程。甚至可以適當(dāng)?shù)叵蜉S向視野外擴(kuò)展,將一部分來自視野外的散射符合事件也考慮進(jìn)來。
蒙特卡羅方法是一類把概率現(xiàn)象作為研究對(duì)象的數(shù)值模擬方法,適用于對(duì)離線系統(tǒng)進(jìn)行計(jì)算仿真。對(duì)于某個(gè)過程,只要我們理清了其背后的概率特性,數(shù)值仿真邏輯并不會(huì)非常復(fù)雜,但往往伴隨著非常大的計(jì)算量。對(duì)于模擬正電子湮滅現(xiàn)象,主要分為以下幾個(gè)步驟,其基本流程如圖4所示。
可以看到圖4中有一個(gè)環(huán)路,其每個(gè)循環(huán)代表一個(gè)湮滅事件的模擬,直到模擬的湮滅事件達(dá)到預(yù)設(shè)值。
圖4 蒙特卡羅方法的基本流程Fig.4 Flow-charts of Monte Carlo method
早在2000年前后,就有學(xué)者意識(shí)到了基于蒙特卡羅方法是一種準(zhǔn)確的散射估計(jì)方法[8],但是受限于當(dāng)時(shí)的硬件水平,其計(jì)算耗時(shí)無法被臨床接受。但隨著近20年來硬件算力的快速提升,特別是GPU計(jì)算性能的提升,使得基于蒙特卡羅方法越來越實(shí)用化[9]。以目前比較新的Nvidia RTX 3080為例,其計(jì)算能力可以達(dá)到兩百萬計(jì)數(shù)每秒的仿真速度或更高。假設(shè)PET系統(tǒng)靈敏度為10 cps/kBq,那么單個(gè)GPU完成蒙特卡羅計(jì)算需要約2200 s。如果使用專業(yè)的顯卡,這個(gè)計(jì)算耗時(shí)可以縮短到幾分鐘甚至更短,基本上可以滿足臨床的要求。
早在20世紀(jì)90年代,3D PET剛剛開始成為主流方向的時(shí)候,就有學(xué)者注意到了視野外部的放射性核素對(duì)于散射的影響[10]。
如1.1節(jié)、1.2節(jié)中所述,目前主流的散射估計(jì)方法都需要核素空間分布信息和衰減分布信息作為輸入。受限于PET有限的軸向視野,往往無法獲知視野外的信息。即便是多個(gè)床位的連續(xù)掃描,也不能得到下一個(gè)床位所在位置的核素空間分布(衰減分布可知,因?yàn)镃T掃描可以先于PET)。所以,如何估計(jì)視野外部分的核素大致分布是估計(jì)視野外散射的關(guān)鍵環(huán)節(jié)。
ANDREYEV等[11]提出了一種以散射估計(jì)結(jié)果為目標(biāo)的方法,其主要算法如圖5所示。該方法使用擴(kuò)展了軸向視野的核素分布圖像和衰減分布圖像作為單散射模擬的輸入,其目標(biāo)是通過不斷更新視野外部分的核素分布圖像,使基于該圖像計(jì)算出的散射分布與實(shí)采數(shù)據(jù)接近。但是該方法存在一個(gè)比較明顯的問題:估計(jì)結(jié)果的精度不可控。因?yàn)檎麄€(gè)圖像更新過程隨機(jī),無法保證結(jié)果收斂;同時(shí)受限于該方法的目標(biāo)ü ü 只關(guān)心散射估計(jì)的結(jié)果,而不關(guān)心視野外圖像的估計(jì)是否準(zhǔn)確,所以其精度很難保證。
圖5 視野外散射估計(jì)方法Fig.5 Estimation of out-of-FOV (field of view) scatter
因?yàn)楣烙?jì)散射的時(shí)刻,PET視野內(nèi)、外的衰減信息已知,視野內(nèi)的核素分布可以估計(jì),所以在實(shí)驗(yàn)中可以通過擬合的方式(如多項(xiàng)式擬合、神經(jīng)網(wǎng)絡(luò)擬合等)計(jì)算出衰減系數(shù)與核素之間的大致映射關(guān)系,粗略估計(jì)出視野外一定距離內(nèi)的核素分布,如式(4)所示。其中θ是擬合參數(shù),μ和I分別是衰減分布圖像和核素分布圖像,下標(biāo)in和out表示軸向視野內(nèi)或外,F(xiàn)(·)是映射函數(shù),L(·)是似然函數(shù)。
在計(jì)算出映射函數(shù)F的參數(shù)之后,我們可以用視野外的衰減分布圖像μout作為函數(shù)輸入,估算出視野外的核素分布圖像Iout。這樣就可以得到軸向視野內(nèi)外的衰減分布圖像和核素分布圖像,然后使用單散射模擬或蒙特卡羅方法進(jìn)行散射估計(jì)。這種方法可能不是最優(yōu),但穩(wěn)定性有保障、計(jì)算量不高,比較適合在臨床場(chǎng)景中使用。
我們做了一組仿真實(shí)驗(yàn),探測(cè)器模型是東軟醫(yī)療的NeuWise Pro PET/CT產(chǎn)品,其軸向視野約為224 mm。放射源的軸向長度分別為200 mm和400 mm,半徑都是100 mm,位于探測(cè)器視野中心。從GATE符合數(shù)據(jù)中取出散射符合,并進(jìn)行單層重組,按層求和,得到如圖6所示的視野外核素對(duì)視野內(nèi)散射分布的影響曲線。
圖6 視野外核素對(duì)視野內(nèi)散射分布的影響Fig.6 Effect on scatter distribution from out-of-FOV nuclide
圖6中曲線橫坐標(biāo)代表不同的軸向?qū)游恢?,縱坐標(biāo)代表每一層上的總散射符合計(jì)數(shù)??梢钥闯霾煌L度的放射源,雖然視野內(nèi)部分的核素分布完全一致,視野外的核素對(duì)視野內(nèi)的散射分布存在明顯的影響。
所以,為了給出一個(gè)準(zhǔn)確且全面的散射校正效果評(píng)估,我們對(duì)單散射模擬和蒙特卡羅方法進(jìn)行了實(shí)驗(yàn),并且這兩種方法又分為考慮視野外散射和不考慮視野外散射兩種情形。
圖7(a)中的實(shí)線是真實(shí)的散射分布;長虛線和斷虛線分別對(duì)應(yīng)單散射模擬和蒙特卡羅方法;星號(hào)或三角標(biāo)記的,表示考慮了視野外的散射。在散射的軸向分布上,如果不考慮視野外的核素衰變產(chǎn)生的散射,那么散射估計(jì)的精度會(huì)呈現(xiàn)從探測(cè)器軸向視野中心到兩邊遞減的趨勢(shì)。這一現(xiàn)象是由康普頓散射角的概率分布導(dǎo)致的:隨著康普頓散射角度的增大,其發(fā)生康普頓散射的概率會(huì)降低。所以視野外部分的散射事件分布在探測(cè)器兩端比較多,中心比較少。在圖7(b)散射的徑向分布對(duì)比中,可以看出無論是單散射模擬,還是蒙特卡羅方法,視野外部分的影響都非常小。在相同的視野條件下(即都考慮或都不考慮視野外散射),因?yàn)镾SS模型是針對(duì)單次散射設(shè)計(jì)的簡化模型,而蒙特卡羅方法則沒有散射次數(shù)的限制,模型與真實(shí)情況更接近。故蒙特卡羅方法比SSS方法更準(zhǔn)確一些。
圖7 不同散射估計(jì)方法在徑向/軸向上散射分布的估計(jì)結(jié)果Fig.7 Axial &radial scatter distribution estimation of different methods
在散射估計(jì)方法的實(shí)際使用中,通常會(huì)在散射擬合過程中使用散射分?jǐn)?shù)對(duì)散射做歸一化處理,所以這里使用KL距離來度量不同散射估計(jì)方法的結(jié)果SE與真實(shí)散射分布ST的差異。
表1是4種散射估計(jì)方法與真實(shí)散射分布的KL距離,越小越好。
表1 不同散射估計(jì)方法與真實(shí)散射分布之間的偏差Tab.1 Deviation between different scatter estimation methods and real scatter distribution
從實(shí)驗(yàn)結(jié)果來看,視野外部分的核素產(chǎn)生的散射事件,如果沒有被散射估計(jì)算法納入模型,則會(huì)引起散射估計(jì)出現(xiàn)偏差。因?yàn)閷?shí)驗(yàn)中的仿真模體在軸向上的形狀、活性濃度變化不大,所以視野外的核素對(duì)散射的軸向分布影響較大,徑向影響較小。但是在實(shí)際臨床掃描條件下,如果視野外部分剛好是活性濃度較高的區(qū)域,如腦部、肝部、腎臟和膀胱等,則很可能會(huì)導(dǎo)致徑向的散射也出現(xiàn)偏差。從理論上來講,這種偏差會(huì)引起圖像均勻性、圖像定量精度等指標(biāo)的下降。
此外,更準(zhǔn)確的模型,如基于蒙特卡羅的散射估計(jì)算法,有助于提高散射分布的估計(jì)的精度。蒙特卡羅(含視野外)方法取得了最接近真實(shí)散射的結(jié)果。
筆者簡要介紹了PET散射校正的發(fā)展歷程。通過對(duì)目前主流方法的原理進(jìn)行剖析,引出了主流方法中存在的問題,如SSS尾部擬合穩(wěn)定性問題、視野外核素散射的影響等。
實(shí)際上,目前主流的散射校正還是局限在單散射的處理上,雙散射甚至多散射的處理技術(shù)尚未鋪開。WATSON等[6]提出了將單散射模擬擴(kuò)展到雙散射模擬的方法,該方法可以比較準(zhǔn)確地計(jì)算出散射次數(shù)不大于2的散射事件分布;相比單(雙)散射模擬法,蒙特卡羅方法減少了空間抽樣、數(shù)據(jù)近似,同時(shí)可以對(duì)探測(cè)器、物質(zhì)、核素甚至電子線路進(jìn)行建模,模型精度更高,并且可以支持大于2次的散射。所以,蒙特卡羅方法是目前業(yè)內(nèi)公認(rèn)的最接近真實(shí)情形的散射計(jì)算方法。隨著算法的進(jìn)步和硬件性價(jià)比的提升,相信這類方法會(huì)出現(xiàn)在商業(yè)化的PET產(chǎn)品中。
對(duì)于PET/CT或PET/MR來講,只要下一個(gè)床位中存在核素分布,就需要考慮來自視野外的散射。在掃描完當(dāng)前床位時(shí),通常我們不清楚下一個(gè)床位位置的核素分布。如何準(zhǔn)確地估計(jì)下一個(gè)床位位置的核素分布,是實(shí)現(xiàn)視野外散射估計(jì)的一個(gè)重要前提。隨著人工智能的快速發(fā)展,相關(guān)技術(shù)已經(jīng)被引入醫(yī)療領(lǐng)域,相信今后會(huì)有更優(yōu)的方法來估算未知的核素分布。此外,超長軸向視野PET的逐步商業(yè)化也使得視野外散射的估計(jì)更加準(zhǔn)確、可信。