陳 曄 肖無(wú)云 李京倫 張羽中 趙 輝 法 鋒

1（國(guó)民核生化災(zāi)害防護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室北京102205）

2（陸軍裝備部北京100072）

γ能譜測(cè)量分析可以實(shí)現(xiàn)快速、可靠的放射性核素定性識(shí)別和定量分析，在核輻射監(jiān)測(cè)、輻射防護(hù)、核物理研究等方面應(yīng)用廣泛［1］。由于γ譜儀系統(tǒng)受探測(cè)器本征能量分辨率的限制，當(dāng)存在能量接近的射線時(shí)，能譜中會(huì)產(chǎn)生重峰，不利于射線能量和強(qiáng)度計(jì)算。因此復(fù)雜譜處理通常需要進(jìn)行重峰解析。

當(dāng)前，雖然已經(jīng)有基于直接解調(diào)、壓縮感知、凸優(yōu)化等理論的全譜反演方法用于γ能譜解析［2-6］，但是該類方法需要對(duì)系統(tǒng)響應(yīng)函數(shù)進(jìn)行精細(xì)刻度，刻度過(guò)程中的誤差會(huì)影響測(cè)量分析結(jié)果；此外，全譜反演方法往往需要更大的工作量和運(yùn)算量，而在很多應(yīng)用場(chǎng)景中，可能只關(guān)注能譜的局部區(qū)域，全譜分析是不必要的［7］。

重峰解析的目的在于從能譜數(shù)據(jù)中確定各個(gè)峰的峰位和峰面積（分別對(duì)應(yīng)射線的能量和強(qiáng)度）等參數(shù)。非線性擬合是一種應(yīng)用廣泛的數(shù)據(jù)分析和參數(shù)估計(jì)方法，它利用優(yōu)化算法最小化擬合函數(shù)與能譜之間的均方誤差（或加權(quán)均方誤差），實(shí)現(xiàn)對(duì)峰參數(shù)的估計(jì)，是一種經(jīng)典的γ能譜重峰解析手段［8］。擬合函數(shù)通常由峰形函數(shù)和本底函數(shù)兩部分組成。例如，對(duì)NaI（Tl）和LaBr3（Ce）等常用閃爍譜儀測(cè)量的γ能譜，其全能峰可用高斯函數(shù)描述，而本底可用線性或者二次函數(shù)的形式［7，9-10］。通常本底函數(shù)可以利用峰區(qū)兩端的數(shù)據(jù)確定［7，11］，或采用統(tǒng)計(jì)敏感的非線性迭代剝峰（Statistics-sensitive Nonlinear Iterative Peak-clipping，SNIP）［12-13］等本底估計(jì)方法獲得［14-15］。因此可提前將本底扣除，后續(xù)的非線性擬合只關(guān)注高斯峰的解析［16］。在優(yōu)化算法上，常用方法包括Levenberg-Marquardt（LM）算 法［14，16-18］、信 賴 域（Trust Region）算 法［11，19］和Nelder-Mead單 純 形法［15］等。

通常情況下，非線性最小二乘可以得到較準(zhǔn)確的擬合結(jié)果。然而，重峰解析問題的復(fù)雜度高，受統(tǒng)計(jì)漲落、重峰個(gè)數(shù)、峰間距等多種因素影響，在統(tǒng)計(jì)漲落水平較高（如測(cè)量時(shí)間短、計(jì)數(shù)少等）或待測(cè)輻射場(chǎng)較復(fù)雜的應(yīng)用場(chǎng)景中，傳統(tǒng)擬合方法的解析精度相對(duì)較低。研究表明：充分利用已知的先驗(yàn)知識(shí)，為擬合問題增加約束或減少待定參數(shù)的個(gè)數(shù)，有利于提高擬合的準(zhǔn)確度［20-21］。為解決傳統(tǒng)方法存在的上述問題，本文提出一種改進(jìn)的重峰解析方法，利用峰寬刻度信息對(duì)非線性最小二乘問題增加約束，并采用迭代求解法，實(shí)時(shí)更新峰寬和權(quán)重等信息，保證解析結(jié)果收斂到更準(zhǔn)確的位置，提高重峰解析精度。

1 方法原理

1.1 約束非線性最優(yōu)化解析重峰

設(shè)扣除本底后的待解析能譜中第i道的計(jì)數(shù)為Si，i1≤i≤iN，i1和iN為峰區(qū)的起始和結(jié)束道址，N為峰區(qū)內(nèi)數(shù)據(jù)點(diǎn)數(shù)。采用的峰區(qū)擬合函數(shù)為高斯峰疊加的形式，如式（1）所示：

式中：Gaj，uj，σj(i)為高斯函數(shù)；aj、uj和σj分別為高斯函數(shù)的高度、中心位置和標(biāo)準(zhǔn)差，如式（2）所示；M為全能峰的個(gè)數(shù)；j為全能峰的序號(hào)，1≤j≤M；θ為待定參數(shù)向量，如式（3）所示，共有3M個(gè)待定參數(shù)。

擬合的目的是求得待定參數(shù)θ的值，從而獲得各個(gè)峰的參數(shù)。在傳統(tǒng)非線性最小二乘方法中，擬合問題轉(zhuǎn)化為最小化加權(quán)均方誤差的形式［8，10］，如式（4）所示。

其中：wi為權(quán)重系數(shù)，根據(jù)能譜計(jì)數(shù)泊松分布的統(tǒng)計(jì)特性，權(quán)重系數(shù)設(shè)置為：

在解析高復(fù)雜度、高統(tǒng)計(jì)漲落的重峰時(shí)，上述傳統(tǒng)方法的誤差較大。針對(duì)該問題，本文使用γ譜儀的峰寬刻度信息對(duì)擬合問題進(jìn)行約束。γ能譜全能峰的半高寬（Full Width at Half Maximum，F(xiàn)WHM）可以表示為能量或者道址的函數(shù)，該函數(shù)通過(guò)實(shí)測(cè)譜刻度得到［22］。設(shè)FWHM與道址i的關(guān)系為函數(shù)F(i)，據(jù)此可將式（4）中待求解的問題變?yōu)橐粋€(gè)帶約束的非線性最小二乘優(yōu)化問題，如式（6）所示：

1.2 迭代求解方法

為了求解式（6）中帶約束的非線性優(yōu)化問題，提出了一種簡(jiǎn)單的迭代方法，具體步驟如下。

第1步：設(shè)置待定參數(shù)的初值θ(0)：

第2步：設(shè)置待定參數(shù)取值范圍的上下界θL、θU：

其 中：aL，j=1，aU，j=2max(Si)，uL，j=i1，uU，j=iN，的上下界都與初值相同，，即每次迭代中，將σj作為常數(shù)）；

第3步：開始迭代擬合，設(shè)n為迭代的次數(shù)；

第4步：在第n次迭代中，利用權(quán)重系數(shù)wi、參數(shù)初值θ(0)和參數(shù)取值范圍θL、θU，使用LM或信賴域等最優(yōu)化方法對(duì)重峰數(shù)據(jù)進(jìn)行非線性最小二乘擬合，得到第n次迭代的結(jié)果，設(shè)為θ(n)：

第5步：判斷是否滿足迭代終止條件：如果連續(xù)兩次迭代擬合結(jié)果之間的歐式距離大于閾值ε，即|θ(n)-θ(n-1)|＞ε，則不滿足終止條件，執(zhí)行第6步；如果|θ(n)-θ(n-1)|≤ε，則滿足迭代終止條件，結(jié)束迭代，最終確定的參數(shù)為θ(n)，執(zhí)行第7步；

第6步：更新參數(shù)并進(jìn)入下一次迭代：首先根據(jù)第n次迭代的擬合結(jié)果θ(n)以及峰半高寬刻度F(i)，更新待定參數(shù)初值、權(quán)重系數(shù)和σj的取值范圍，然后令n=n+1，返回第4步，繼續(xù)進(jìn)行下一次迭代；

第7步：根據(jù)擬合結(jié)果計(jì)算各個(gè)全能峰的峰位、峰半高寬和峰面積信息，完成重峰的解析：設(shè)pj、FWHMj和Aj分別為第j個(gè)全能峰的峰位、峰半高寬和峰面積，則有

2 仿真驗(yàn)證

首先使用仿真數(shù)據(jù)驗(yàn)證方法的性能。使用MATLAB軟件，在400～600道范圍內(nèi)，根據(jù)峰位、峰面積和FWHM參數(shù)仿真產(chǎn)生兩個(gè)疊加的高斯峰，然后添加泊松噪聲，得到仿真γ能譜重峰數(shù)據(jù)。

能譜計(jì)數(shù)的統(tǒng)計(jì)漲落水平、峰之間的距離以及峰之間的面積比等參數(shù)不同，重峰解析的難度也不同。因此在生成仿真數(shù)據(jù)時(shí)，調(diào)整上述參數(shù)，研究其對(duì)方法性能的影響。由于能譜計(jì)數(shù)服從泊松分布，總計(jì)數(shù)越大，統(tǒng)計(jì)性越好。因此可通過(guò)調(diào)整全能峰面積模擬統(tǒng)計(jì)漲落的影響。

每種條件下生成2 000組重峰數(shù)據(jù)，解析得到峰位、FWHM以及峰面積信息，統(tǒng)計(jì)其平均誤差，并與傳統(tǒng)方法的結(jié)果對(duì)比。為求解式（4）和式（6）中的最優(yōu)化問題，采用了MATLAB中的fit函數(shù)，其底層使用的最優(yōu)化方法為信賴域算法。

圖1 展示了三組典型的仿真重峰數(shù)據(jù)及其解析結(jié)果。雖然兩種方法得到的擬合曲線都與模擬譜符合，但是重峰解析結(jié)果卻存在較大差別?？梢?，擬合曲線與譜數(shù)據(jù)的符合程度并不能準(zhǔn)確反映重峰解析的質(zhì)量，應(yīng)該使用峰位、峰面積、FWHM等峰參數(shù)的誤差判斷解析結(jié)果的優(yōu)劣。

2.1 統(tǒng)計(jì)漲落的影響

固定兩個(gè)高斯峰的峰位分別為480道和520道，兩者的峰面積相同，并由100逐步增加到10 000。解析得到的峰位、FWHM和峰面積的平均誤差隨峰面積的變化情況分別如圖2（a）、圖2（b）和圖2（c）所示?？梢钥闯?，隨著峰面積的增大，由于能譜計(jì)數(shù)的統(tǒng)計(jì)漲落減少，兩種方法的誤差都呈下降趨勢(shì)，該結(jié)果與理論預(yù)期一致。同時(shí)，本方法解析結(jié)果的誤差一直低于傳統(tǒng)方法。特別是在峰面積較小，即數(shù)據(jù)的統(tǒng)計(jì)漲落較高時(shí)，本方法可以有效提高解析精度，在峰面積為100時(shí)，峰位誤差降低了27%，峰面積相對(duì)誤差降低了71%，說(shuō)明其抗統(tǒng)計(jì)漲落能力強(qiáng)。

圖2 （d）為兩種方法解析結(jié)果對(duì)應(yīng)的目標(biāo)函數(shù)取值情況。可見，傳統(tǒng)方法的擬合曲線與模擬譜的殘差比本方法低；然而，圖2（a）至圖2（c）的結(jié)果已經(jīng)證明本方法解析的峰參數(shù)精度更高，這說(shuō)明擬合殘差小并不代表實(shí)際解析效果好，傳統(tǒng)方法存在過(guò)擬合問題。本方法利用了峰寬刻度信息為最優(yōu)化問題增加了約束，雖然擬合殘差比傳統(tǒng)方法高，但是能夠得到更符合實(shí)際物理規(guī)律的結(jié)果，從而提升解析精度。在峰間距影響的研究結(jié)果圖3（d）中也得到了類似的結(jié)果，這些結(jié)果符合理論預(yù)期，驗(yàn)證了為優(yōu)化問題增加約束的效果。

圖1 不同條件下的仿真重峰數(shù)據(jù)及其解析結(jié)果示例(a)峰面積為100，峰間距40道，(b)左右兩峰面積分別為1 000和100，峰間距40道，(c)峰面積為10 000，峰間距20道Fig.1 Examples of simulated data and analysis results under different conditions(a)The peak area is 100,the peak distance is 40 channels,(b)The peak area of the left and right peaks are 1 000 and 100 respectively,the peak distance is 40 channels,(c)The peak area is 10 000,the peak distance is 20 channels

圖2 峰參數(shù)解析誤差隨峰面積的變化情況(a)平均峰位誤差，(b)平均FWHM誤差，(c)平均峰面積相對(duì)誤差，(d)平均目標(biāo)函數(shù)值Fig.2 Variation of peak parameter resolving errors with peak area(a)Average peak position error,(b)Average FWHM error,(c)Average peak area relative error,(d)Average objective function value

2.2 峰間距的影響

固定兩個(gè)高斯峰的峰面積為10 000，研究峰間距（10～30道）對(duì)于解析結(jié)果的影響。由于第500道的FWHM值為32.1道，以此為基準(zhǔn)，峰間距變化范圍為0.31～0.93倍FWHM。圖3（a）、圖3（b）和圖3（c）分別為峰位、FWHM和峰面積的平均誤差隨峰間距的變化情況。峰間距越小，解析復(fù)雜度和難度就越大。因此隨著峰間距的減小，兩種方法的誤差逐漸增大。與圖2中的結(jié)果類似，本方法的誤差低于傳統(tǒng)方法，在峰間距為10道，即0.31倍FWHM時(shí)，峰位誤差降低了77%，峰面積相對(duì)誤差降低了72%。

圖3 峰參數(shù)解析誤差隨峰間距的變化情況(a)平均峰位誤差，(b)平均FWHM誤差，(c)平均峰面積相對(duì)誤差，(d)平均目標(biāo)函數(shù)值Fig.3 Variation of peak parameter resolving errors with peak distance(a)Average peak position error,(b)Average FWHM error,(c)Average peak area relative error,(d)Average objective function value

2.3 峰面積比的影響

固定兩個(gè)高斯峰的峰位分別為480道和520道，分別稱之為峰1和峰2。固定峰1的面積為1 000，改變峰2的面積，范圍由100～10 000，即峰2與峰1的峰面積比的范圍為0.1～10。

圖4 （a～f）分別為解析得到的兩個(gè)峰的峰位、FWHM和峰面積的平均誤差隨峰面積比的變化情況?？梢钥闯觯瑢?duì)于峰2，由于其峰面積一直在增大，峰參數(shù)解析誤差呈現(xiàn)出隨峰面積增大而減小的趨勢(shì)，與圖2中結(jié)果一致；峰1的情況與峰2相反，其峰參數(shù)誤差隨著峰面積比的增大比呈現(xiàn)出增大的趨勢(shì)?？梢?，雖然峰1的面積一直固定，但是附近峰的峰面積增大同樣會(huì)導(dǎo)致峰1的解析誤差增大，這說(shuō)明在強(qiáng)峰背景下，弱峰的解析難度較大。當(dāng)峰1和峰2的面積分別為1 000和100時(shí)，相比于傳統(tǒng)方法，本方法解析的峰1和峰2的峰位誤差分別降低了26%和49%，峰面積相對(duì)誤差分別降低了60%和59%。

3 實(shí)驗(yàn)驗(yàn)證

為了驗(yàn)證該方法對(duì)實(shí)測(cè)譜的解析能力，分別采用SAINT-GOBAIN公司的BrilLanCeTM380 76S76/3.5型LaBr3（Ce）探 測(cè) 器 和CANBERRA公 司 的BICRONTM802-3X3型NaI（Tl）探測(cè)器，配合ORTEC digiBASE-E數(shù)字多道，獲取實(shí)測(cè)能譜。兩個(gè)探測(cè)器在661.7 keV的能量分辨率分別為3.1%和7.3%。實(shí)驗(yàn)用放射源為152Eu和133Ba點(diǎn)源，活度分別約為7 000 Bq和25 000 Bq，源與探測(cè)器的距離為25 cm，能譜測(cè)量時(shí)間為30 min。其中，使用LaBr3（Ce）探測(cè)器 測(cè) 量 了152Eu的γ能 譜，并 對(duì) 其1 085.8 keV、1 112.1 keV射線［23］的2重峰區(qū)進(jìn)行解析，得到的擬合曲線和重峰解析結(jié)果分別如圖5（a）和圖5（b）所示；使用NaI（Tl）探測(cè)器測(cè)量了133Ba的γ能譜，并對(duì)其276.4 keV、302.9 keV、356.0 keV和383.8 keV射線［23］的4重峰區(qū)進(jìn)行解析，得到的擬合曲線和重峰解析結(jié)果分別如圖5（c）和圖5（d）所示。與圖1中的結(jié)果類似，在圖5（a）和圖5（c）中，傳統(tǒng)方法和本方法得到的擬合曲線非常接近，且都與實(shí)測(cè)譜保持一致；但是由圖5（b）和圖5（d）可以看出，兩種方法解析出的全能峰，特別是133Ba重峰區(qū)內(nèi)較弱的三個(gè)峰，存在較為明顯的差別。

圖4 峰參數(shù)解析誤差隨峰面積比的變化情況(a)峰1的平均峰位誤差，(b)峰2的平均峰位誤差，(c)峰1的平均峰面積相對(duì)誤差，(d)峰2的平均峰面積相對(duì)誤差，(e)峰1的平均FWHM誤差，(f)峰2的平均FWHM誤差Fig.4 Variation of peak parameter resolving errors with peak area ratio(a)Average peak position error of peak 1,(b)Average peak position error of peak 2,(c)Average peak area relative error of peak 1,(d)Average peak area relative error of peak 2,(e)Average FWHM error of peak 1,(f)Average FWHM error of peak 2

圖5 實(shí)測(cè)γ能譜重峰解析結(jié)果(a)152Eu能譜2重峰擬合曲線，(b)152Eu能譜2重峰解析結(jié)果，(c)133Ba能譜4重峰擬合曲線，(d)133Ba能譜4重峰解析結(jié)果Fig.5 Overlapping peak analysis results of measuredγspectrum(a)Fitting curve of 2 overlapping peaks in 152Eu spectrum,(b)Resolved 2 overlapping peaks in 152Eu spectrum,(c)Fitting curve of 4 overlapping peaks in 133Ba spectrum,(d)Resolved 4 overlapping peaks in 133Ba spectrum

表1 實(shí)測(cè)γ能譜重峰的峰位和FWHM解析結(jié)果Table 1 Peak position and FWHM analysis results of measured spectra

表2 實(shí)測(cè)γ能譜重峰的峰面積解析結(jié)果Table 2 Peak area analysis results of measured spectra

表1 列出了解析得到的峰位和FWHM，表2列出了解析得到的峰面積。在表1和表2中，為了評(píng)價(jià)解析質(zhì)量，利用國(guó)際原子能機(jī)構(gòu)（International Atomic Energy Agency，IAEA）提供的放射性核素?cái)?shù)據(jù)庫(kù)［23］，獲得了152Eu和133Ba的γ射線能量和分支比信息；再利用譜儀系統(tǒng)的能量刻度和FWHM刻度信息，確定了各個(gè)全能峰的峰位、FWHM和峰面積的參考值。其中，峰面積與γ射線的分支比和譜儀的探測(cè)效率有關(guān)，由于重峰區(qū)內(nèi)射線能量接近，探測(cè)效率變化不大，因此主要考慮峰面積與射線分支比的對(duì)應(yīng)關(guān)系。因?yàn)闇?zhǔn)確的峰面積值無(wú)法獲得，所以在每個(gè)峰區(qū)內(nèi)，將分支比最大的峰作為基準(zhǔn)，計(jì)算其它峰與該峰的峰面積比，用于評(píng)價(jià)峰面積的準(zhǔn)確度。此外，除了1 085.8 keV和1 112.1 keV射線，152Eu核素還有一個(gè)1 089.7 keV的γ射線，分支比為0.017 3，在LaBr3（Ce）能譜中無(wú)法與1 085.8 keV射線區(qū)分，因此在計(jì)算峰面積比時(shí)將1 089.7 keV射線也考慮在內(nèi)，如表2所示。

由表1整體看來(lái)，傳統(tǒng)方法解析實(shí)測(cè)譜的峰位和FWHM與參考值誤差較大，特別是存在低能量峰的FWHM值大于高能量峰的FWHM值的情況，不符合實(shí)際規(guī)律；而本方法解析的峰位和FWHM的準(zhǔn)確度在多數(shù)情況下優(yōu)于傳統(tǒng)方法。

表2 中兩種方法解析的峰面積比的相對(duì)誤差如圖6所示，本方法的峰面積比相對(duì)誤差相比于傳統(tǒng)方法平均降低了80%。NaI（Tl）探測(cè)器的能量分辨性能相對(duì)較差，而133Ba的4重峰區(qū)由于峰個(gè)數(shù)多，計(jì)數(shù)相對(duì)較少，且峰面積差別大，解析的難度較大?？梢钥闯鰝鹘y(tǒng)方法對(duì)其解析的峰面積誤差明顯大于152Eu的2重峰區(qū)；而本方法對(duì)兩個(gè)峰區(qū)解析的峰面積誤差沒有明顯區(qū)別，都保持在較低水平。該結(jié)果驗(yàn)證了本方法比傳統(tǒng)方法解析復(fù)雜能譜的能力更強(qiáng)。

圖6 峰面積比的相對(duì)誤差對(duì)比Fig.6 Comparison of the relative errors of peak area ratio

4 結(jié)語(yǔ)

γ能譜重峰的解析受到統(tǒng)計(jì)漲落、峰間距、峰面積比等多種因素的影響。實(shí)驗(yàn)結(jié)果表明：較高的統(tǒng)計(jì)漲落水平、較小的峰間距、較大的峰面積差距以及較多的重峰個(gè)數(shù)都會(huì)導(dǎo)致重峰解析問題復(fù)雜度和難度的增加。

得益于對(duì)峰寬刻度等已知信息的充分利用，在處理高復(fù)雜度、高統(tǒng)計(jì)漲落的能譜數(shù)據(jù)時(shí)，通過(guò)將非線性擬合的結(jié)果約束在更準(zhǔn)確的范圍內(nèi)，可有效提高重峰解析的精度。這說(shuō)明，盡可能多地利用已知的物理?xiàng)l件或先驗(yàn)知識(shí)，為擬合問題增加約束，有利于增強(qiáng)抗統(tǒng)計(jì)漲落能力，使結(jié)果收斂到更準(zhǔn)確的位置。

在傳統(tǒng)重峰擬合解析方法中，往往把擬合殘差等作為判斷擬合結(jié)果好壞的依據(jù)。然而，殘差實(shí)際上反映的是由擬合結(jié)果重建的曲線與測(cè)量數(shù)據(jù)在整體上的吻合程度；本文研究結(jié)果表明，殘差值小并不能充分說(shuō)明擬合結(jié)果與真實(shí)值的誤差一定小。因此，對(duì)解析結(jié)果的評(píng)價(jià)除了要考慮擬合殘差外，還要考慮其是否與已知條件一致，是否符合實(shí)際的物理規(guī)律。單純追求殘差最小化，容易陷入過(guò)擬合。

應(yīng)該注意到，本方法適用于NaI（Tl）、LaBr3（Ce）等具有對(duì)稱的高斯峰形的探測(cè)器。對(duì)于HPGe、CdZnTe等半導(dǎo)體探測(cè)器，其譜峰一般具有明顯的低能拖尾，在使用本方法時(shí)應(yīng)根據(jù)峰形特點(diǎn)調(diào)整約束條件。