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

趙 日，劉立業，曹勤劍

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

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

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

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

1 原理和實驗

1.1 能譜重建技術原理

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

(1)

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

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

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

y=Rx+ε

(2)

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

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

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

1.2 Gold算法原理

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

RTy=RTRx

(3)

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

Mx=y*

(4)

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

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

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

1.3 實驗過程

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

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

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

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

(5)

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

圖4 探測系統完整數字模型Fig.4 Complete digital model of detecting system

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

另外，F8計數卡的能量范圍為20～1 600 keV，能量節點數為1 500個，節點能量值則與實測能譜在20～1 600 keV區間內的各道能量對應。最終計算所得的測量矩陣R(即響應能譜)如圖5所示，其維數為1 500×800。

2 實驗結果

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

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

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

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

另外，實驗發現，當1

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

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

(6)

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

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

(7)

其中，下標1、2分別代表不同射線。

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

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

3 討論

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

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

能譜重建法不僅能準確識別所有主要γ射線，且沒有任何誤識別情況，同時所有核素的活度計算相對誤差均小于10%。總而言之，在全身計數器γ能譜分析中，能譜重建法明顯優于全能峰法。這是易于理解的，因為全能峰法只著眼于能譜的局部形態和計數，且掌握的其他信息極其有限，僅為能量、效率、峰形刻度函數，而能譜重建法不局限于部分能譜，而是將全譜作為分析對象，并通過測量矩陣將探測場景所有已知的物理、幾何信息(如探測器、被測對象的結構、尺寸、材料等及射線與物質相互作用過程等)均包含其中。基于更豐富的全系統信息，其分析結果必然更準確。

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

4 總結

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

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