改進(jìn)辛普森-SNIP算法在航空γ儀器譜本底扣除中的應(yīng)用

楊 津 李 飛 葛良全 張慶賢 熊茂淋 孫 坤 谷 懿

（成都理工大學(xué) 核技術(shù)與自動化工程學(xué)院 成都 610059）

航空γ能譜測量（Airborne Gamma-ray Spectrometry Survey，AGS）是指將航空γ能譜儀系統(tǒng)安裝在飛機(jī)、飛艇或其他飛行器上，測量地空界面的γ射線能譜，從而獲得地表介質(zhì)中放射性鈾、釷、鉀等核素的含量（比活度）。該方法也可以應(yīng)用在貴金屬礦、稀土礦、多金屬礦礦床等與放射性核素相關(guān)的其他礦產(chǎn)的普查上［1］；近年來被廣泛應(yīng)用于巖性地質(zhì)填圖、地下水資源搜尋、放射性環(huán)境污染監(jiān)測等諸多方面［2-4］。通過航空γ儀器譜的道址-計數(shù)，可以計算當(dāng)前譜線反映的放射性核素的種類和含量。因此，航空γ儀器譜的解譜工作一直是該領(lǐng)域研究的重點(diǎn)和難點(diǎn)。其中對全譜的本底扣除是解譜過程中最重要的步驟之一［5］。本底扣除的可信度和精確度將直接影響到后續(xù)的尋峰、峰面積計算、核素含量反演等。

理想的扣除本底方法應(yīng)滿足以下特征：1）原理簡單、運(yùn)算快速且本底扣除的設(shè)定參數(shù)少；2）計算機(jī)能進(jìn)行快速自動處理且運(yùn)行成本低；3）能有效扣除因康普頓效應(yīng)對譜峰計數(shù)造成的影響。針對散射本底扣除這一步驟，國內(nèi)外提出了大量方法，其中包括統(tǒng)計敏感的非線性迭代剝峰算法（Statisticssensitive Nonlinear Iterative Peak-clipping，SNIP）及其相關(guān)衍生算法［6］、傅里葉變換法［7-8］、剝峰法［9］、迭代法［10］等。已知最簡便的方法是線性法，即將能量特征峰的兩波谷用一條直線連接，將其視為基線，基線以下為本底部分，但該方法會將有效計數(shù)當(dāng)作本底扣除，相對誤差大，只適用于粗略估計。SNIP方法被認(rèn)為是其中的一種較為理想的選擇，因為它具有原理簡單、運(yùn)行速度快、成本低等優(yōu)點(diǎn)，國外核數(shù)據(jù)計算系統(tǒng)中也已將它運(yùn)用于全譜本底評估和扣除。張慶賢等［11］結(jié)合SNIP本底扣除方法與傅里葉變換低通濾波方法提出了傅里葉變換本底逐步逼近法。然而，SNIP算法容易受到峰寬參數(shù)的影響，采用固定峰區(qū)寬度進(jìn)行本底扣除會帶來較大的誤差；其次，算法迭代過程中的收斂速度緩慢，導(dǎo)致運(yùn)行時間延長。因此，本文基于具有較高代數(shù)精度的辛普森公式的思想，針對原始算法中的不足，提出了新的改進(jìn)方法，并將其應(yīng)用于航空γ儀器譜本底扣除中。改進(jìn)后的SNIP算法在保留原始算法簡潔、高效的優(yōu)點(diǎn)的同時，減少了迭代次數(shù)，提高了運(yùn)行效率，使得譜峰的本底扣除效果和準(zhǔn)確度進(jìn)一步改善。

1 本底來源分析

通常而言，本底是指除被測量的放射源之外的其他能夠在能譜區(qū)間造成計數(shù)的干擾因素，如放射性污染、宇宙射線及儀器的電子學(xué)噪聲等。本底扣除是任何輻射測量中不可避免的問題，如果要提高放射性核素含量和放射性活度分析的準(zhǔn)確性，則必須有效地識別和扣除與全能峰面積無關(guān)的計數(shù)［12］。

γ能譜中的本底主要來自康普頓散射效應(yīng)和γ射線在探測器靈敏體積內(nèi)產(chǎn)生的小角度散射以及其他干擾因素在能譜的能量區(qū)間造成的無用計數(shù)［13］。由于γ射線同物質(zhì)發(fā)生光電效應(yīng)、康普頓散射和電子對效應(yīng)等一系列相互作用，當(dāng)γ射線同周邊散射和探測器自身發(fā)生散射，產(chǎn)生能量更低的散射光子和反沖電子［14］。當(dāng)散射光子被探測器探測或者散射光子逃離探測器時，輸出的γ譜中會出現(xiàn)康普頓平臺［15］。此外一些外在的環(huán)境因素，例如環(huán)境中的放射性核素、電子儀器引起的噪聲、宇宙射線、電磁干擾等會在能譜中產(chǎn)生環(huán)境本底［16-17］。而放射源的內(nèi)在因素會導(dǎo)致放射性本底的出現(xiàn)，因為在放射源中除了待測核素外，還存在大量其他的無關(guān)核素，如果這些核素也有放射性，就會產(chǎn)生計數(shù)干擾［18］。在對γ能譜進(jìn)行分析的過程中，散射本底加大了γ能譜的復(fù)雜程度，使得一些弱全能峰難以被識別，從而降低能譜分析數(shù)據(jù)的精度［19］。在γ能譜的分析過程中，核心問題是要提高能譜儀的測量精度，其關(guān)鍵在于選擇合適的本底扣除方法，才能有效抑制康普頓散射對譜線的影響。

2 改進(jìn)型辛普森-SNIP方法

2.1 SNIP算法及其改進(jìn)

目前，SNIP算法被認(rèn)為是最優(yōu)的本底扣除算法之一。在1988年，Ryan等［20］首次提出SNIP算法扣除γ能譜本底的原理，但SNIP算法在細(xì)節(jié)描述方面，很大程度受到峰區(qū)寬度參數(shù)的影響，當(dāng)選用不同的參數(shù)的時候，扣除本底的效果存在差異。隨后Morhá?等［21］在 1997年運(yùn)用LLS運(yùn)算符（taking the Square root and using the Ln operator twice）改進(jìn)了在弱峰下的SNIP算法，在運(yùn)用SNIP算法之前對每一道址的能譜計數(shù)進(jìn)行變換。在2009年，Morhá?等［22］將原始SNIP算法中的函數(shù)用高階濾波函數(shù)替代，來緩解因?qū)⒃糞NIP算法應(yīng)用于不同能量區(qū)時，固定窗寬相同而產(chǎn)生的本底扣除誤差，從而得到一個更精確的本底扣除效果。2012年，尹旺明等［23］通過借助半峰全寬（Full Width at Half Maxima，F(xiàn)WHM）刻度對SNIP濾波窗寬度進(jìn)行自適應(yīng)調(diào)整，進(jìn)一步解決固定窗寬帶來的影響。2019年，Li等［24］針對SNIP算法的演變給出了準(zhǔn)確評價，提出新的算法演變思路，并指出基于SNIP的改進(jìn)算法，都在避免因變換寬度而出現(xiàn)扣除本底過程中在峰的高能端與低能端出現(xiàn)的較大誤差。而在SNIP光滑的本身，可以采取迭代的思想優(yōu)化。

利用SNIP算法進(jìn)行本底扣除，求解凈峰面積主要分為4個過程。首先需要對能譜的每一道址計數(shù)變換，采用LLS算符（式（1）），通過自然對數(shù)運(yùn)算壓縮計數(shù)范圍，計算高計數(shù)譜線，同時平方根運(yùn)算可以提取弱峰信息。

式中：i表示道址；y（i）表示i道址的計數(shù)；V（i）表示經(jīng)LLS運(yùn)算后的變換值。

其次，利用迭代的方式，依次計算出V1（i）、V2（i）、…、Vm（i）。假設(shè)迭代到p次，則通過比較當(dāng)前Vp-1（i）與Vp-1（i+p）和Vp-1（i-p）的平均值的大小，取兩者中較小值作為下一次的迭代值（式（2））。

式中：p表示第p次迭代，即從1開始每次迭代運(yùn)算后加1直到等于給定的m值為止（變換窗寬逐增法）；Vp（i）表示第p次迭代的LLS變換值。

迭代完畢后，再通過逆LLS變換（式（3）），求得本底譜數(shù)據(jù)。即作式（1）的逆運(yùn)算：對變換值y（i）使用兩次以e為底的指數(shù)運(yùn)算，再進(jìn)行求平方。

最后，用原始譜數(shù)據(jù)減去本底譜數(shù)據(jù)，得到γ能譜的凈峰計數(shù)譜（式（4））。

本底扣除方法示意圖如圖1所示。

圖1 SNIP本底扣除方法示意圖Fig.1 The schematic diagram of SNIP background subtraction algorithm

在第二步SNIP算法迭代的過程中，實(shí)質(zhì)利用的是矩形公式的思想。雖然計算簡單，方法可靠，但是收斂速度太慢。本文采用辛普森公式的思想，即利用區(qū)間二等分的三個點(diǎn)來進(jìn)行積分插值，其科特斯系數(shù)分別為1/6、4/6、1/6，由此來替代第二步的SNIP迭代計算。由于辛普森公式具有更高的代數(shù)精度，計算結(jié)果的準(zhǔn)確度較之SNIP迭代更優(yōu)秀。在迭代過程中，設(shè)定步長為單位距離，經(jīng)過LLS轉(zhuǎn)換的數(shù)值為對應(yīng)的函數(shù)值（式（5））。

2.2 測試數(shù)據(jù)

本文采用XPG-3000A航空γ能譜勘查系統(tǒng)的實(shí)測數(shù)據(jù)。系統(tǒng)采用兩箱探測器，每箱由兩條40 cm×10 cm×5 cm陣列式NaI晶體探測器組成。單晶體能量分辨率在137Cs的0.661 MeV＜8.0%，能譜峰漂（208Tl的2.62 MeV）＜0.5道。能量范圍最大1～10 MeV，最大線性計數(shù)率可達(dá)200 ks-1（kilo counts per second），每道最高計數(shù)264。

應(yīng)用XPG-3000A系統(tǒng)對137Cs人工核素源進(jìn)行測試響應(yīng)。137Cs源活度為9 795 Bq，為了更好體現(xiàn)特征γ信息，采用探測器近距離測量的方式，即雷達(dá)高度為0 m。由于閃爍探測器的陣列式結(jié)構(gòu)，γ射線與物質(zhì)反應(yīng)在探測器內(nèi)部的電子沉積情況不盡相同，因此采用不同的探測點(diǎn)位分析，并記錄譜線數(shù)據(jù)。每次實(shí)驗測量時間為200 s，總測量次數(shù)為50次。

3 結(jié)果與討論

Morhá?等［21］在1997年提出單峰凈峰計數(shù)保留率r和本底扣除率t兩個系數(shù)，用于評價算法的本底扣除效果。

式中：PNC（Peak Net Counts）為扣除本底后的凈峰計數(shù)；PNC0為未扣除本底的凈峰計數(shù)；PTC（Peak Total Counts）為扣除本底后的總峰計數(shù)。

本文采用的是137Cs放射性核素的特征航空γ譜數(shù)據(jù)，其特征峰對應(yīng)能量為0.661 MeV。在確定m值的時候，一般采用公式：

式中：W表示濾波窗的寬度，即變換寬度；m為SNIP迭代的次數(shù)。對于不同寬度W的單能峰來說，m值取在（W-1）/2附近，能獲得最佳的譜峰本底扣除效果，且不會造成譜峰形狀的改變。如圖2所示，本文運(yùn)用SNIP算法對譜線進(jìn)行自適應(yīng)扣除，先確定特征峰峰位道址為393道和譜峰區(qū)間為355～433道，再根據(jù)特征峰的峰區(qū)確定SNIP的迭代次數(shù)。

由圖2，迭代次數(shù)會影響譜線的本底的剝離結(jié)果，當(dāng)?shù)螖?shù)m=7時，會導(dǎo)致原始譜出現(xiàn)“過度剝離”現(xiàn)象，SNIP法能削弱本底對特征峰的貢獻(xiàn)，并減少本底干擾因素（如本底峰、屏蔽材料產(chǎn)生的特征X射線等）的影響，但同時過度削減了全能峰面積；當(dāng)?shù)螖?shù)m=12時，譜線的本底扣除會出現(xiàn)“欠剝離”現(xiàn)象，SNIP方法不僅不能有效抑制本底計數(shù)的影響，還會因本底扣除不徹底影響到凈峰計數(shù)。由此，當(dāng)?shù)螖?shù)m=9時最為合理，能夠有效抑制康普頓散射對譜線的影響。

圖2 不同迭代次數(shù)下SNIP本底扣除效果Fig.2 Background subtraction effect shown on spectrum for SNIP algorithm with different iteration times

為了驗證改進(jìn)算法和原始的SNIP算法的本底扣除效果，本文引用了近年來常用的剝峰法［25-26］，同時作為本底扣除效果的參照。為減少剝峰法的剝離次數(shù)，提高運(yùn)行效率，在剝離前對譜數(shù)據(jù)進(jìn)行取對數(shù)，剝峰后，反變換可得到本底的形狀。由于核素活度和元素含量的定量分析與相應(yīng)特征峰的凈峰面積有很直接的關(guān)系，故本文采用這三種方法對譜線進(jìn)行本底扣除后得到特征峰譜線（圖3），并通過對比扣除本底后得到的特征峰面積與實(shí)際凈峰面積的平均相對誤差來評價本底扣除效果（表1）。圖3中曲線由上至下分別表示原始譜線用高斯函數(shù)擬合后的特征峰、經(jīng)過原始SNIP計算得到的凈計數(shù)峰、改進(jìn)SNIP計算后的凈計數(shù)峰以及經(jīng)剝峰法處理得到的凈計數(shù)峰。最后任意選取15次不同時間下測得的譜線，通過單峰凈峰計數(shù)保留率r和本底扣除率t兩個系數(shù)，用于進(jìn)一步評價原始SNIP算法、改進(jìn)SNIP算法以及剝峰法這三種算法的本底扣除效果（表2）。

表1 經(jīng)本底扣除后得到的特征峰面積以及平均相對誤差Table 1 The characteristic peak area and relative error table obtained after deducting from the background

通過比較三種方法得出的峰面積與實(shí)際特征峰面積的平均相對誤差，由表1可知，采用原始SNIP法、改進(jìn)SNIP法、剝峰法得出的峰面積的平均相對誤差分別為7.60%、4.42%、-9.72%，改進(jìn)方法的準(zhǔn)確度較原始SNIP提高了41.8%。由表2可知，在本底扣除效果上，改進(jìn)SNIP算法的r范圍為97.79%～99.83%，原始SNIP算法的r范圍為96.98%～99.45%，剝峰法的r范圍為96.51%～99.18%，平均值分別為98.71%、98.29%、97.93%。改進(jìn)SNIP算法的t范圍為96.59%～99.47%，原始的SNIP算法的t范圍為95.69%～98.78%，剝峰法的t范圍為94.73%～98.76%，平均值分別為97.87%、97.15%、97.06%。通過對比，證明了在本文中改進(jìn)的SNIP算法對單峰的凈計數(shù)保留和本底扣除方面相對原始SNIP算法有細(xì)小提升，對特征峰的本底扣除效果最好，但兩種方法的評價系數(shù)r和t均超過95%，具有足夠高的精度，而剝峰法在該方面要稍次于這兩種方法。

SNIP算法最大的優(yōu)勢是運(yùn)算簡單、迅速，因此需要對改進(jìn)的SNIP算法的工作效率做評價。通過記錄改進(jìn)前后兩種算法的迭代次數(shù)，其中原始的SNIP算法的最優(yōu)迭代次數(shù)為12次，而基于辛普森公式迭代的改進(jìn)型SNIP算法的最優(yōu)迭代次數(shù)僅為9次，相對于原始的SNIP算法減少了三次，有一定的提升。進(jìn)一步通過tic-toc的計時方式對兩種算法的運(yùn)行時間進(jìn)行評價。

圖3 扣除本底后的凈計數(shù)峰Fig.3 The comparison of the two background subtraction algorithms

圖4 比較SNIP算法改進(jìn)前后的運(yùn)行時間Fig.4 Running time comparison before and after the SNIP algorithm improvement

如圖4所示，得到改進(jìn)SNIP算法運(yùn)行時間為53.99 ms，而原始SNIP算法運(yùn)行時間為55.39 ms，從對比中可以發(fā)現(xiàn)，改進(jìn)型SNIP算法比原始SNIP算法的運(yùn)行效率提高了2.53%，但仍然有較高的運(yùn)算速度。

表2 三種算法本底扣除評價系數(shù)Table 2 The evaluation coefficient of background subtraction of the three algorithms

4 結(jié)語

本文利用辛普森公式替換SNIP算法中的迭代公式，對SNIP算法進(jìn)行改進(jìn)并應(yīng)用在航空γ探測的儀器譜本底扣除中，并得到以下結(jié)論：

1）通過比較原始SNIP算法、改進(jìn)SNIP算法以及剝峰法得出的峰面積與實(shí)際特征峰面積的平均相對誤差，發(fā)現(xiàn)改進(jìn)SNIP法的相對誤差為4.42%，其準(zhǔn)確度相比于原始算法提高了41.8%。

2）對比三種方法的評價系數(shù)，證明改進(jìn)的SNIP算法對單峰的凈計數(shù)保留和本底扣除方面相對原始SNIP算法有細(xì)小提升，本底扣除效果更準(zhǔn)確，且改進(jìn)前后的兩種算法在該方面均優(yōu)于剝峰法。

3）在實(shí)際應(yīng)用中，改進(jìn)SNIP算法完整保存了原始算法簡潔、高效的優(yōu)勢，同時選取了合適的變換寬度來達(dá)到縮減誤差的效果，且最優(yōu)迭代次數(shù)減少三次，運(yùn)行效率提高了2.53%。

4）本文通過對算法的改進(jìn)，在完成γ譜預(yù)處理工作的同時，對后續(xù)核素的定量分析提供了良好的基礎(chǔ)。