基于能譜重建技術(shù)的全身計數(shù)器γ能譜分析

趙 日，劉立業(yè)，曹勤劍

(中國輻射防護研究院，山西 太原 030006)

全身計數(shù)器是專用于測量人體體內(nèi)放射性的儀器。它從體外直接探測體內(nèi)發(fā)出的γ射線，通過分析γ射線能譜實現(xiàn)對放射性核素種類和含量的定量測量。全身計數(shù)器在輻射防護、核安全等領(lǐng)域均有重要應(yīng)用[1]。

目前，全身計數(shù)器的γ能譜分析最常用的方法是全能峰法[2-5]。該方法先用濾波算法尋找能譜中由光電效應(yīng)形成的全能峰，然后對每個峰進行曲線擬合以求得其凈計數(shù)(峰面積)，最后，結(jié)合能量和效率刻度函數(shù)確定測量對象發(fā)出的γ射線的能量和數(shù)量。然而，由于全身計數(shù)器通常采用NaI(Tl)探測器以提高探測效率，而該類型探測器能量分辨率差，且存在溫漂等現(xiàn)象，再加上人體對γ射線的散射以及天然40K干擾等因素，全身計數(shù)器獲取的γ能譜常呈現(xiàn)全能峰辨識度低、重峰現(xiàn)象嚴(yán)重、峰下側(cè)基底高、有假峰等不利性態(tài)[6]。對于這種γ能譜，全能峰法得到的結(jié)果通常準(zhǔn)確性很低：其在尋峰時極易誤識別和漏識別，峰擬合時因重峰和高基底的原因而誤差較大。如何提高全身計數(shù)器γ能譜分析精度一直是長期受關(guān)注的課題。

近期，一種基于能譜重建算法的分析技術(shù)在低分辨率γ能譜中得到越來越多的應(yīng)用。該技術(shù)首先建立核素出射譜到探測器測量譜直接的測量矩陣，然后使用特殊迭代算法從測量譜反解出出射譜。由于出射譜在全譜范圍內(nèi)只在若干能量點處有幅值，這些能量點與γ射線能量一一對應(yīng)，而幅值等于測量時間內(nèi)該γ射線總發(fā)射數(shù)，因此，理論上可直接從反解結(jié)果中讀取γ射線能量和數(shù)量，以此實現(xiàn)對放射性核素的定量分析。該技術(shù)在NaI探測器解譜領(lǐng)域的應(yīng)用始于20世紀(jì)末：1997年Bandzuch等[7]對正電子湮滅譜使用了能譜重建技術(shù)；2000年Meng等[8]在對7.62 cm×7.62 cm NaI(Tl)探測器的γ能譜進行分析時比較了3種不同重建算法的效果。2004年Jandel等[9]對152Eu和56Co的測量譜進行了能譜重建，結(jié)果顯示，全能峰被壓縮入孤立的道，而基底部分被全部移除；2009年Rahman等[10]則將該技術(shù)應(yīng)用于核電廠流出物監(jiān)測系統(tǒng)(使用NaI(Tl)探測器)所獲取的γ能譜，所得到的重建能譜中峰康比大幅提高；2010年Hanka[11]提出了基于動態(tài)指標(biāo)集和梯度法的γ能譜重建算法；Morhac等[12]綜述了目前常用的若干種γ能譜重建算法。本文為提高全身計數(shù)器γ能譜分析精度，開展基于能譜重建技術(shù)的γ能譜分析方法研究。

1 原理和實驗

1.1 能譜重建技術(shù)原理

假設(shè)核素出射的γ能譜為x，探測器測得的γ能譜為y，則理論上有：

(1)

其中，r(E′→E)表示出射譜中能量為E′的單位數(shù)量γ射線在測量譜E處的相應(yīng)計數(shù)。

當(dāng)能譜計數(shù)率較高時，會出現(xiàn)符合相加效應(yīng)和死時間現(xiàn)象，式(1)并不嚴(yán)格成立；不過對于全身計數(shù)器測量，人體體內(nèi)放射性一般非常弱，加上本底輻射是天然環(huán)境輻射，此時γ能譜計數(shù)率較低，符合相加效應(yīng)和死時間現(xiàn)象幾乎可忽略，因此可認(rèn)為式(1)是γ射線出射譜和探測譜之間映射關(guān)系的合理描述。

將式(1)寫為離散形式，并考慮到各道計數(shù)存在統(tǒng)計漲落，則有下式：

y=Rx+ε

(2)

其中：R為測量矩陣，其元素Rij表示出射譜中單位數(shù)量的第i道γ射線在測量譜第j道處的相應(yīng)計數(shù)；ε為由統(tǒng)計漲落帶來的誤差。

圖1 直接矩陣求逆所得的解Fig.1 Solution from direct matrix invert method

由于測量矩陣R高度病態(tài)，且式中存在誤差項，式(2)若用一般的矩陣求逆算法求解，只能得到劇烈振蕩的無意義解，如圖1所示。為從式(2)出發(fā)，穩(wěn)定求解出x且保證x的所有分量均為非負(fù)(能譜幅值不可能是負(fù)值)，提出了若干精細(xì)算法，其中最常用、最有效的主要有兩種，即Gold算法[12-13]和EM算法[14-16]。兩者均屬于迭代算法，且均具有全局收斂性和解非負(fù)性的優(yōu)點，但Gold算法較EM算法收斂更快。本文將采用基于Gold算法的能譜重建技術(shù)。

1.2 Gold算法原理

Gold算法的目標(biāo)是求解式(2)的最小二乘法解，即求解式(3)：

RTy=RTRx

(3)

令y*=RTy，M=RTR，則式(3)可改寫為：

Mx=y*

(4)

基于式(4)，Gold算法的迭代流程可表示如下。

1) 取初值x(0)=[1,1,…,1]T，確定最大迭代次數(shù)K和停機準(zhǔn)則δ。

在上述Gold算法流程基礎(chǔ)上，Jandel等[9]又進一步提出了boosted-Gold算法，其原理是將Gold算法作為內(nèi)循環(huán)，將若干步內(nèi)循環(huán)得到的x進行激勵操作，即令xi=|xi|p(i=1,…,n)，其中p(p>1)為用戶選擇的激勵因子，然后將新的x作為下一次內(nèi)循環(huán)的初始值，如此往復(fù)，直到外循環(huán)達(dá)到預(yù)設(shè)次數(shù)。

1.3 實驗過程

本文選擇了中國輻射防護研究院劑量學(xué)實驗室自主研制的一臺立式NaI型全身計數(shù)器作為實驗對象。首先，利用點源對該全身計數(shù)器進行了能量和全能峰峰形刻度，獲得了該設(shè)備的能量刻度函數(shù)和全能峰展寬函數(shù)。然后，使用該全身計數(shù)器對兩套完全相同但灌裝不同核素的BOMAB人體模型(下文以BOMAB-A和BOMAB-B代表)進行測量。BOMAB模型由中國輻射防護研究院研制，每套由10個不同形狀大小的空心柱形聚乙烯容器組成以代表人體，聚乙烯容器內(nèi)充滿放射性水溶液。BOMAB-A灌裝了4種活度已知的57Co、134Cs、137Cs、60Co放射性核素水溶液，而BOMAB-B僅灌裝了水。測量兩套BOMAB模型是為了更準(zhǔn)確扣減γ能譜中天然輻射本底的貢獻：將BOMAB-A測量譜減去BOMAB-B測量譜(按測量時間的比例折算)則可得到只由BOMAB-A中4種核素貢獻的凈譜，而將天然本底的計數(shù)貢獻消除。本文分別對兩套模型進行2 000 s的測量，經(jīng)過扣減后，最終所得BOMAB-A凈γ能譜如圖2所示。能譜共2 048道，能量范圍為0～2 MeV，圖2還標(biāo)注了能譜中各全能峰對應(yīng)的核素名稱和γ射線能量。

圖2 實測BOMAB-A凈γ能譜Fig.2 Measured net γ spectrum for BOMAB-A

a——探頭接PMT端的透視圖；b——探頭遠(yuǎn)離PMT端的透視圖圖3 NaI探頭的X光透視成像Fig.3 X perspective imaging of NaI detector

采用蒙特卡羅模擬方法計算式(2)中的測量矩陣R。計算方法如下：對該全身計數(shù)器和被測BOMAB人體模型進行完整的數(shù)字建模，其中全身計數(shù)器的支撐結(jié)構(gòu)等對γ射線探測過程影響較小的部分未進行建模；為保證建模精度，對全身計數(shù)器中兩個大體積NaI探頭進行X光透視成像，如圖3所示，通過影像獲得細(xì)節(jié)結(jié)構(gòu)的準(zhǔn)確尺寸。圖4示出了該測量場景的數(shù)字模型，其中，正對體模的為探測模塊，由兩個縱向排列的探測器和屏蔽結(jié)構(gòu)組成?；跀?shù)字模型，使用MCNPX 2.5.1進行測量矩陣的計算。計算時，源粒子在BOMAB模型中均勻分布，且能量在50～1 500 keV等間隔變化，共計算800個能量點，每次模擬的源粒子數(shù)為107個，能譜的展寬函數(shù)使用了實測函數(shù)，形式如下：

(5)

其中：FWHM為全能峰半高寬，keV；E為全能峰峰中心的能量，keV。

圖4 探測系統(tǒng)完整數(shù)字模型Fig.4 Complete digital model of detecting system

圖5 蒙特卡羅模擬得到的響應(yīng)能譜Fig.5 Measure matrix constructed by Monte Carlo simulation

另外，F(xiàn)8計數(shù)卡的能量范圍為20～1 600 keV，能量節(jié)點數(shù)為1 500個，節(jié)點能量值則與實測能譜在20～1 600 keV區(qū)間內(nèi)的各道能量對應(yīng)。最終計算所得的測量矩陣R(即響應(yīng)能譜)如圖5所示，其維數(shù)為1 500×800。

2 實驗結(jié)果

利用所得的BOMAB-A凈γ能譜和測量矩陣R，使用Gold算法求解核素出射譜。圖6比較了實測的BOMAB-A凈γ能譜與重建能譜，其中，實測能譜的幅值放大了500倍?？梢姡亟茏V中高幅值區(qū)域與實測能譜中的全能峰位置吻合，其中，57Co核素的122 keV的全能峰收縮至單獨1道，其余核素的全能峰也均有不同程度的收縮，同時，全能峰區(qū)域外則幾乎不再有計數(shù)。這表明了重建能譜在全能峰收縮和基底消除方面體現(xiàn)了一定優(yōu)勢。然而，這一重建結(jié)果仍不夠理想，因為核素出射譜的形態(tài)應(yīng)為：所有全能峰均只占據(jù)1道而其余道計數(shù)均為0，從而能從中直接讀取γ射線能量和數(shù)量，但基于目前的重建結(jié)果卻無法實現(xiàn)直接讀取定量結(jié)果。

圖6 Gold算法的重建結(jié)果Fig.6 Reconstruction result given by Gold algorithm

為改善重建效果，本文進一步使用boosted-Gold算法進行計算。

當(dāng)參數(shù)p=3時，boosted-Gold算法在第15次外迭代時到達(dá)停機準(zhǔn)則。圖7對比顯示了實測BOMAB-A凈γ能譜與boosted-Gold算法的重建能譜。由圖可見，此時重建能譜已達(dá)理想狀態(tài)，即全譜只有全能峰位置有計數(shù)，且所有全能峰均收縮至1道。

另外，實驗發(fā)現(xiàn)，當(dāng)1

圖7 boosted-Gold算法的重建結(jié)果Fig.7 Reconstruction result given by boosted-Gold algorithm

由boosted-Gold算法的重建能譜，可直接讀取各γ射線的能量和計數(shù)，而核素的活度可通過式(6)得到：

(6)

其中：A為核素活度；N為出射能譜中γ射線的計數(shù)；T為測量時長；η為該γ射線的分支比。

當(dāng)核素有多個γ射線時，其活度可由最小二乘法確定。如134Cs有604.7、797.0 keV兩條射線，60Co有1 173.2、1 332.5 keV兩條射線，則它們的活度可由式(7)求得：

(7)

其中，下標(biāo)1、2分別代表不同射線。

為比較能譜重建法與全能峰法的分析結(jié)果，使用Genie 2000軟件對實測BOMAB-A凈γ能譜進行分析。Genie 2000由Canberra公司開發(fā)并在核探測領(lǐng)域被廣泛使用，其算法采用全能峰法[17]。

最終的計算結(jié)果列于表1、2。

表1 兩種方法的γ射線識別能力比較Table 1 Comparison of γ identification ability between two methods

表2 兩種方法的核素活度計算能力比較Table 2 Comparison of nuclide activity calculation ability between two methods

由表1可見，在γ射線能量確定方面，能譜重建法與全能峰分析法精度相當(dāng)，給出的γ射線能量與實際能量誤差在±3 keV以內(nèi)。但值得注意的是，Genie 2000軟件給出的結(jié)果中還包括一實際并不存在的84.2 keV的峰值。而由表2可見，在核素活度確定方面，能譜重建算法則顯著地優(yōu)于全能峰分析法，具體來說，對于BOMAB-A中的4種核素，能譜重建算法的活度計算值與實際活度的相對誤差小于10%，最小僅為-2.2%，而Genie 2000軟件給出的活度結(jié)果的誤差除60Co外均大于10%，最大達(dá)到17%。另外，兩種方法在57Co核素的活度計算結(jié)果上差異最大。

3 討論

表1顯示Genie 2000軟件誤識別出84.2 keV能量的全能峰，這印證了全能峰法在峰識別環(huán)節(jié)的弊端。全能峰法使用濾波法來識別全能峰，具體對于Genie 2000軟件，其使用高斯函數(shù)二階導(dǎo)函數(shù)作為濾波器[17]，濾波器展寬由全能峰展寬函數(shù)(由峰形刻度實驗給出)確定。當(dāng)能譜局部的形態(tài)與該處高斯函數(shù)相近時，算法判定存在全能峰。參照圖2可見，在84.2 keV能量附近，實測能譜呈現(xiàn)近似全能峰的形態(tài)，這也是Genie 2000誤識別的原因。然而，這一假峰結(jié)構(gòu)事實上是由大體積源(人體模型)的低能散射峰效應(yīng)和57Co核素122 keV能量γ射線的康普頓散射效應(yīng)等綜合因素造成的。由于原理所限，全能峰法無法辨別其真?zhèn)巍?/p>

表2中，Genie 2000軟件計算核素的活度時誤差較大，這也是由全能峰法的原理所致。該方法計算全能峰凈計數(shù)時需扣減峰下側(cè)的基底，當(dāng)扣減的基底與實際情況不符時，就會導(dǎo)致峰凈計數(shù)偏大或偏小，從而造成對核素活度估算的不準(zhǔn)確。對于實測的BOMAB-A凈γ能譜，由圖2可知，57Co核素122 keV全能峰處于康普頓散射計數(shù)堆積區(qū)域，下側(cè)基底幅值較高且形態(tài)復(fù)雜，難以準(zhǔn)確扣減，而一旦扣減不準(zhǔn)確就會導(dǎo)致峰凈計數(shù)誤差很大，這也就是其對57Co核素的活度計算相對誤差達(dá)到17.0%的原因；同時，能譜中134Cs的604.7 keV全能峰與137Cs的661.7 keV全能峰相互重疊，形成重峰，其下側(cè)基底同樣較為復(fù)雜，因此，其給出的這兩個核素活度的相對誤差也都超過了10%；而60Co核素的兩個能量全能峰相互間重疊較少，且處在能譜高能端，下側(cè)基底幅值低、形態(tài)簡單，所以對其活度估算的誤差最小。

能譜重建法不僅能準(zhǔn)確識別所有主要γ射線，且沒有任何誤識別情況，同時所有核素的活度計算相對誤差均小于10%?？偠灾?，在全身計數(shù)器γ能譜分析中，能譜重建法明顯優(yōu)于全能峰法。這是易于理解的，因為全能峰法只著眼于能譜的局部形態(tài)和計數(shù)，且掌握的其他信息極其有限，僅為能量、效率、峰形刻度函數(shù)，而能譜重建法不局限于部分能譜，而是將全譜作為分析對象，并通過測量矩陣將探測場景所有已知的物理、幾何信息(如探測器、被測對象的結(jié)構(gòu)、尺寸、材料等及射線與物質(zhì)相互作用過程等)均包含其中?；诟S富的全系統(tǒng)信息，其分析結(jié)果必然更準(zhǔn)確。

當(dāng)然，能譜重建法也有其明顯的缺點：它的實施依賴于高精度的測量矩陣，而這必須通過細(xì)致的探測器表征及建模工作來實現(xiàn)，但該工作極為耗時耗力，使得該方法的應(yīng)用成本相對傳統(tǒng)方法要大；測量矩陣包含了被測對象的物理、幾何信息以及其與探測器的相對位置信息，這意味著上述條件稍加改變，測量矩陣就會發(fā)生變化，故該方法只適合于被測對象單一且其與探測器相對位置固定的測量場景；復(fù)雜的迭代計算過程導(dǎo)致該方法無法從理論上對其計算結(jié)果給出定量的誤差估計，而僅能通過與真值的比較而給出后驗分析，這無疑也影響了其實施的可靠性。

4 總結(jié)

本文開展了將能譜重建技術(shù)應(yīng)用于全身計數(shù)器γ能譜分析的實驗研究。實驗中，首先使用一臺立式NaI型全身計數(shù)器測量兩套BOMAB人體模型(一套為含源模型，一套為本底模型)，通過扣減本底譜得到了只含模型中放射性核素計數(shù)貢獻的γ能譜；然后對該全身計數(shù)器和人體模型進行精細(xì)數(shù)字建模，并使用蒙特卡羅模擬方法計算了該系統(tǒng)的測量矩陣；基于實測能譜和測量矩陣，使用Gold算法和boosted-Gold算法分別求解出射能譜；最后，將求解結(jié)果與Genie 2000軟件分析結(jié)果進行了比較。實驗結(jié)果顯示，單純的Gold算法效果不佳，而boosted-Gold算法在參數(shù)p取1～10時能較準(zhǔn)確重建出理論的出射能譜。根據(jù)boosted-Gold算法所求結(jié)果，能實現(xiàn)對所有核素的準(zhǔn)確識別，且活度計算相對誤差均小于10%。這一結(jié)果顯著優(yōu)于Genie 2000軟件。

本文還深入分析了能譜重建法在核素定量分析方面優(yōu)于全能峰法的機制。這些結(jié)論顯示了基于能譜重建技術(shù)的γ能譜分析方法在全身計數(shù)器測量以及其他核探測領(lǐng)域進行實際應(yīng)用的潛力。