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

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

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

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

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

1 原理和實(shí)驗(yàn)

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

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

(1)

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

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

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

y=Rx+ε

(2)

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

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

由于測(cè)量矩陣R高度病態(tài)，且式中存在誤差項(xiàng)，式(2)若用一般的矩陣求逆算法求解，只能得到劇烈振蕩的無意義解，如圖1所示。為從式(2)出發(fā)，穩(wěn)定求解出x且保證x的所有分量均為非負(fù)(能譜幅值不可能是負(fù)值)，提出了若干精細(xì)算法，其中最常用、最有效的主要有兩種，即Gold算法[12-13]和EM算法[14-16]。兩者均屬于迭代算法，且均具有全局收斂性和解非負(fù)性的優(yōu)點(diǎn)，但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和停機(jī)準(zhǔn)則δ。

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

1.3 實(shí)驗(yàn)過程

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

圖2 實(shí)測(cè)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

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

(5)

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

圖4 探測(cè)系統(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計(jì)數(shù)卡的能量范圍為20～1 600 keV，能量節(jié)點(diǎn)數(shù)為1 500個(gè)，節(jié)點(diǎn)能量值則與實(shí)測(cè)能譜在20～1 600 keV區(qū)間內(nèi)的各道能量對(duì)應(yīng)。最終計(jì)算所得的測(cè)量矩陣R(即響應(yīng)能譜)如圖5所示，其維數(shù)為1 500×800。

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

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

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

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

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

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

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

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

(6)

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

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

(7)

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

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

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

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

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

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

3 討論

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

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

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

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

4 總結(jié)

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

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