約束非線性最優(yōu)化迭代方法解析γ能譜重峰

陳 曄 肖無云 李京倫 張羽中 趙 輝 法 鋒

1（國民核生化災(zāi)害防護國家重點實驗室北京102205）

2（陸軍裝備部北京100072）

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

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

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

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

1 方法原理

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

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

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

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

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

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

1.2 迭代求解方法

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

第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)化方法對重峰數(shù)據(jù)進行非線性最小二乘擬合，得到第n次迭代的結(jié)果，設(shè)為θ(n)：

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

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

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

2 仿真驗證

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

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

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

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

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

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

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

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

圖3 峰參數(shù)解析誤差隨峰間距的變化情況(a)平均峰位誤差，(b)平均FWHM誤差，(c)平均峰面積相對誤差，(d)平均目標函數(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 峰面積比的影響

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

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

3 實驗驗證

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

圖4 峰參數(shù)解析誤差隨峰面積比的變化情況(a)峰1的平均峰位誤差，(b)峰2的平均峰位誤差，(c)峰1的平均峰面積相對誤差，(d)峰2的平均峰面積相對誤差，(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 實測γ能譜重峰解析結(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 實測γ能譜重峰的峰位和FWHM解析結(jié)果Table 1 Peak position and FWHM analysis results of measured spectra

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

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

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

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

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

4 結(jié)語

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

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

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

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