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

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

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

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

1 本底來(lái)源分析

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

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

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ù)的時(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算法之前對(duì)每一道址的能譜計(jì)數(shù)進(jìn)行變換。在2009年，Morhá?等［22］將原始SNIP算法中的函數(shù)用高階濾波函數(shù)替代，來(lái)緩解因?qū)⒃糞NIP算法應(yīng)用于不同能量區(qū)時(shí)，固定窗寬相同而產(chǎn)生的本底扣除誤差，從而得到一個(gè)更精確的本底扣除效果。2012年，尹旺明等［23］通過(guò)借助半峰全寬（Full Width at Half Maxima，F(xiàn)WHM）刻度對(duì)SNIP濾波窗寬度進(jìn)行自適應(yīng)調(diào)整，進(jìn)一步解決固定窗寬帶來(lái)的影響。2019年，Li等［24］針對(duì)SNIP算法的演變給出了準(zhǔn)確評(píng)價(jià)，提出新的算法演變思路，并指出基于SNIP的改進(jìn)算法，都在避免因變換寬度而出現(xiàn)扣除本底過(guò)程中在峰的高能端與低能端出現(xiàn)的較大誤差。而在SNIP光滑的本身，可以采取迭代的思想優(yōu)化。

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

通過(guò)比較三種方法得出的峰面積與實(shí)際特征峰面積的平均相對(duì)誤差，由表1可知，采用原始SNIP法、改進(jìn)SNIP法、剝峰法得出的峰面積的平均相對(duì)誤差分別為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%。通過(guò)對(duì)比，證明了在本文中改進(jìn)的SNIP算法對(duì)單峰的凈計(jì)數(shù)保留和本底扣除方面相對(duì)原始SNIP算法有細(xì)小提升，對(duì)特征峰的本底扣除效果最好，但兩種方法的評(píng)價(jià)系數(shù)r和t均超過(guò)95%，具有足夠高的精度，而剝峰法在該方面要稍次于這兩種方法。

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

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

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

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

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

4 結(jié)語(yǔ)

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

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

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

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

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