基于高斯銳化法的重疊峰分解方法研究

王青山， 王冬陽， 張雄杰， 湯 彬， 吳和喜

核技術(shù)應(yīng)用教育部工程研究中心(東華理工大學(xué))， 江西 南昌 330013

引 言

γ能譜測量與分析是放射性信息獲取的重要方式， 在輻射監(jiān)測、 放射醫(yī)療、 巖性分析等領(lǐng)域[1]都有廣泛的應(yīng)用。 實際測量中γ能譜通常會受到環(huán)境本底與測量干擾等不可抗力因素的影響， 通常為多種能量和強度不同的γ射線單能譜的疊加， 特別是能量相近的γ射線， 往往以重疊峰的形式出現(xiàn)， 且強度弱的γ譜線又容易被強峰或本底所掩蓋， 容易對放射性核素種類的判斷及含量的計算造成不利影響。

重疊峰的數(shù)值解析技術(shù)至關(guān)重要， 近年來， 國內(nèi)外在重疊峰的分解方面進(jìn)行了較為深入的研究， 如導(dǎo)數(shù)法[2]、 正交投影法[3]、 小波分析法[4]、 遺傳算法[5-6]、 人工神經(jīng)網(wǎng)絡(luò)模型[6]、 純元素譜剝離[7]、 粒子群算法[8]等方法。 都月[9]等改進(jìn)了縱向和平分迭代法， 提出優(yōu)化迭代重疊峰分解方法， 較傳統(tǒng)迭代方法其效率提高了66.7%； 黃洪全等[10]利用高斯混合模型結(jié)合期望最大化迭代算法進(jìn)行γ譜重疊峰分解； 周世融等[11]提出了峰銳化法結(jié)合雙樹復(fù)小波變換的低分離度重疊峰解析方法， 并實現(xiàn)了分離度為0.34的X熒光光譜重疊峰分解。

目前， 重疊峰的分解方法眾多且較為成熟， 但是普遍存在卷積過程復(fù)雜、 計算資源消耗大、 結(jié)果受算法參數(shù)影響較大的弊端， 不適用于測量現(xiàn)場能譜快速解析的實際應(yīng)用。 本文提出一種快速、 低耗的分解重疊峰的新方法， 為放射性能譜測量中重疊峰的現(xiàn)場分解提供支持。

1 高斯銳化方法(GSM)

1.1 峰銳化方法

銳化算法常被應(yīng)用于圖像處理領(lǐng)域中以提高圖像的對比度或清晰度， 也可在信號處理領(lǐng)域中用以人為地提高峰的可視分辨率。 在能譜分析中， 對重疊峰進(jìn)行銳化處理可人為地增大其分離度， 達(dá)到分離重疊峰的目的， 通常采用二階導(dǎo)數(shù)或四階導(dǎo)數(shù)來確定峰位置[12]， 其算法可表示如式(1)

R(x)=s(x)-ks(2)(x)

(1)

式(1)中，R(x)為銳化后的信號，s(x)為原始信號，s(2)(x)為其二階導(dǎo)數(shù)，k為銳化因子， 可平衡信號信噪比和基線平坦度從而調(diào)整其大小達(dá)到最佳銳化效果。

1.2 褶積滑動變換方法

由于核衰變及測量的統(tǒng)計性， 能譜數(shù)據(jù)普遍受到統(tǒng)計漲落的影響， 一般通過譜光滑處理以增加定性及定量分析的準(zhǔn)確性。 與傅里葉變換法的“頻域濾波”不同， 褶積滑動變換法為“時域濾波”， 即為變換函數(shù)與實驗譜線進(jìn)行卷積變換[13]

(2)

式(2)中，F(xiàn)(x)為平滑后的譜線，f(x)為原始譜線，g(x)為變換函數(shù)，m為滑動變換的半窗寬。 隨著x的增加，g(x)將逐道向前滑動，f(x)經(jīng)過寬度為2m+1的變換從而轉(zhuǎn)換為新譜線F(x)， 當(dāng)變換函數(shù)與峰形狀函數(shù)相同且半寬度一致時， 光滑效果最優(yōu)。

1.3 高斯銳化方法

銳化的本質(zhì)為對信號的微分處理， 當(dāng)一個噪聲信號經(jīng)過微分后， 信號的噪聲同樣經(jīng)過了微分。 因此， 不可利用傳統(tǒng)的銳化方法直接對含噪的信號進(jìn)行運算。 若對譜線平滑進(jìn)行處理可以提供更高的分辨率， 根據(jù)卷積的性質(zhì)對式(2)變形， 信號F(x)的n階導(dǎo)數(shù)可表示為

F(n)(x)=f(n)(x)*g(x)=f(x)*g(n)(x)

(3)

由式(3)可以看出： 對原信號而言， 平滑與微分操作二者的先后順序并不影響其計算結(jié)果。 為了取得平滑后信號的微分， 我們可先對變換函數(shù)進(jìn)行微分得到其平滑算子， 然后利用該平滑算子作為濾波器對原信號進(jìn)行卷積。 若g(n)(x)選擇高斯導(dǎo)數(shù)小波， 則信號F(x)的n階導(dǎo)數(shù)便可通過小波變換實現(xiàn)， 高斯小波變換的本質(zhì)和信號f(x)經(jīng)過g(x)平滑的n階導(dǎo)數(shù)運算類似。

將銳化后的高斯函數(shù)作為高斯銳化算子g(x)可得

(4)

其中，A為高斯銳化算子的峰值，μ為高斯銳化算子的峰位，σ為高斯銳化算子的標(biāo)準(zhǔn)差，k為高斯銳化因子。 為使信號卷積變換前后峰面積保持不變， 高斯銳化算子g(x)作為平滑窗函數(shù)還需滿足如式(5)關(guān)系

(5)

式(5)中，m為半窗寬。 對高斯銳化算子按式(5)進(jìn)行歸一化處理后， 分別取不同的k和σ值， 以探究算子參數(shù)對高斯銳化成形效果的影響， 其結(jié)果如圖1所示。k值越大， 其峰值越大， 峰寬基本保持不變， 且始終過定點(μ±σ，Ae-1/2)， 峰兩側(cè)旁瓣的凹陷也越低， 且當(dāng)k

一葉一滴汗，一果一情懷。高品質(zhì)需要新技術(shù)的支撐，多年來，農(nóng)場不斷推廣新技術(shù)，為特色、高效、生態(tài)的芒果產(chǎn)業(yè)夯實了根基。

圖1 高斯銳化算子(A=1, μ=0, m=200)

結(jié)合式(2)與式(4)可得高斯銳化法(Gauss sharpening method, GSM)遞推公式為

(6)

本方法處理噪聲信號具有一定的優(yōu)越性， 既能過濾銳化所引入的噪聲又能分辨峰位置、 突出峰形， 實際應(yīng)用中可根據(jù)實際情況選擇進(jìn)行一次或多次高斯銳化成形。

2 仿真實驗

2.1 重疊峰模型構(gòu)建及銳化效果驗證

由于核輻射測量數(shù)據(jù)服從高斯分布， 采用多個高斯函數(shù)的疊加以仿真能譜重疊峰

(7)

式(7)中，n為重疊的高斯峰個數(shù)，Ai為高斯峰的高度，μi為高斯峰的峰位，σi為高斯峰的標(biāo)準(zhǔn)差， noise為疊加的噪聲信號， 以仿真譜線中的統(tǒng)計漲落影響。

重疊峰中不同組分的峰高、 峰寬和峰位置參數(shù)將影響峰的重疊程度， 通常采用分離度[11]來判斷相鄰兩峰的分離程度

(8)

式(8)中，μ1和μ2分別為兩峰的峰位置，σ1和σ2分別為兩峰的峰寬。 即使沒有噪聲的信號， 當(dāng)分離度過小時， 重疊的兩個峰也很難分辨， 本文實驗所用分離度均取0.3以上， 并模擬一組重疊峰(如圖2所示)

(9)

其分離度為0.55， 并利用MATLAB中“awgn”函數(shù)添加60 dB高斯白噪聲。 由于采用的信號是已知特征參數(shù)的峰信號， 因此可通過比較處理結(jié)果與原特征參數(shù)來衡量峰分辨效果。

圖2 重疊峰示意圖(SNR=60 dB)

根據(jù)傳統(tǒng)銳化成形式(1)和高斯銳化成形式(6), 分別對式(9)所示的重疊峰進(jìn)行處理， 其結(jié)果如圖3所示。 可明顯看出經(jīng)過銳化的信號其峰寬變窄、 峰幅值變大、 分離度變大、 峰形狀更加明顯且獨立。 傳統(tǒng)銳化成形后的信號分辨率明顯下降， 其信噪比僅為19.5 dB， 而高斯銳化法成形結(jié)果的分辨率有了明顯的提高， 其信噪比約為52 dB。

圖3 重疊峰兩種銳化結(jié)果對比

高斯銳化法具有較好的噪聲抑制性能， 重疊峰分離效果顯著， 且成形前后峰面積的標(biāo)準(zhǔn)差約為0.008 4， 可認(rèn)為其成形前后峰面積保持不變， 為能譜的定量分析提供實驗結(jié)果支撐。

2.2 基于非線性曲線擬合的重疊峰參數(shù)提取

對給定數(shù)據(jù)， 使得擬合函數(shù)與其誤差平方和最小化原則， 尋求與給定點的距離平方和最小的曲線h(x)， 經(jīng)過GSM成形的信號服從式(7)分布， 經(jīng)過化簡變形可得

(10)

式(10)中，AF，μF，kF和σF均為擬合參數(shù)， 本文使用MATLAB中l(wèi)sqcurvefit工具， 以式(10)作為期望函數(shù)進(jìn)行非線性擬合， 并對其進(jìn)行峰位、 峰面積計算， 以驗證高斯銳化算法的有效及準(zhǔn)確性。

在非線性擬合中初始值的選擇非常重要， 若選擇不當(dāng)， 則計算結(jié)果不一定收斂。 為了避免了曲線擬合的不唯一性， 提高曲線擬合結(jié)果可信度。 本文以各個子峰峰位為中心， 取7個點或9個點進(jìn)行擬合， 得到各個子峰的初始參數(shù)， 再將其作為曲線擬合的初始值， 擬合結(jié)果如圖4所示。 模擬的標(biāo)準(zhǔn)高斯重疊峰均實現(xiàn)完全分解， 且峰位和峰面積的相對誤差分別在0.006%和0.03%以內(nèi)。 結(jié)果表明， 峰銳化法結(jié)合褶積滑動變換可快速且較準(zhǔn)確地分解標(biāo)準(zhǔn)高斯重疊峰。

圖4 重疊峰GSM非線性擬合結(jié)果

表1 標(biāo)準(zhǔn)高斯重疊峰模型的GSM峰位分解結(jié)果/%

表2 標(biāo)準(zhǔn)高斯重疊峰模型的GSM峰面積分解結(jié)果/%

利用高斯銳化法可對能譜段進(jìn)行重疊峰分解并提取相關(guān)參數(shù)。 對于高分辨率(SNR≥50 dB)譜線， 即使分離度為0.375時也可以達(dá)到優(yōu)良的分解效果， 誤差可控制在1%以內(nèi)； 對于低分辨率譜線(SNR≤30 dB)， 其誤差均高于5%， 且在分離度為0.375時已經(jīng)完全失真。 結(jié)合上述分析， 利用該方法可對SNR≥40dB，RS≥0.375的譜線實現(xiàn)較好的分離效果。

3 基于蒙特卡洛方法的模擬能譜解析

根據(jù)國際原子能機構(gòu)(IAEA)提供的核數(shù)據(jù)資料，214Bi核素的γ衰變能量為609.32 keV，137Cs核素的γ衰變能量為661.675 keV， 二者相差52.355 keV；206Bi核素的γ衰變能量為1 718.7 keV，26Al核素的γ衰變能量為1 808.65 keV， 二者相差89.95 keV。 當(dāng)混合放射源中包含上述核素時， 特別是在使用NaI(Tl)閃爍探測器等低分辨率的探測器時， 能譜譜峰將嚴(yán)重重疊。

為較好地仿真實驗條件， 利用MCNP對131I，137Cs，214Bi，206Bi和26Al組合放射源進(jìn)行模擬， 根據(jù)實驗時所用儀器的探頭， 模型中采用φ30 mm×50 mm的圓柱型NaI(Tl)晶體作為探頭材料， 密度為3.67 g·cm-3。 晶體外側(cè)包裹一層0.05 cm厚的MgO反射層， 底部的SiO2光學(xué)玻璃厚度為0.2 cm， 側(cè)面和前面用0.2 cm的Al殼包裹， 探頭外圍空間用空氣填滿， 其結(jié)構(gòu)如圖5所示。

圖5 探測器模型

利用F8計數(shù)卡計算放射源的γ射線在NaI(Tl )晶體中的脈沖高度能譜分布， 所測5種放射源的γ能譜如圖6所示， 求得131I譜峰(364.489 keV)、137Cs譜峰(661.675 keV)、214Bi譜峰(609.32和1 120.294 keV)、206Bi譜峰(1 718.7 keV)和26Al譜峰(1 808.65 keV)峰位分為451， 757， 823， 2 145和2 258道。

圖6 基于MCNP的混合源模擬γ能譜

該混合源γ能譜圖存在兩組重疊峰， 且完全無法識別出137Cs核素全能峰的位置。 對原始譜線直接進(jìn)行GSM處理可快速進(jìn)行全譜分析并對其中重疊峰進(jìn)行分解， 無需進(jìn)行能譜平滑以及扣除本底操作。 GSM全譜分解的結(jié)果如圖7所示， 分解誤差如表3所示。

圖7 混合源模擬γ能譜的GSM分解結(jié)果

表3 模擬γ能譜的GSM數(shù)值解析結(jié)果

對于單峰而言， 處理前后峰位近似不變； 對于重疊峰而言， 分解前后誤差在1～3道之內(nèi)， 二者相對誤差均在1%以內(nèi)。 實驗結(jié)果證明， 利用高斯銳化法能夠快速、 高效地對γ能譜進(jìn)行全譜解析且重疊峰分解結(jié)果較為精確。

4 結(jié) 論

提出一種高斯銳化法快速全譜分析并分解重疊峰的新方法。 仿真實驗結(jié)果表明， 該方法具備較好的噪聲抑制性能， 且能對SNR≥40 dB，RS≥0.375的重疊峰有較好的分解效果。 同時， 利用該方法對MCNP模擬的混合源γ能譜進(jìn)行全譜解析， 實現(xiàn)了放射性核素全能峰的有效分解。

該方法的特點和優(yōu)勢為： (1)無需進(jìn)行能譜平滑及扣除本底操作， 即使在康普頓坪的影響下也可得到較準(zhǔn)確的結(jié)果； (2)計算資源消耗低， 能譜變換僅需一次滑動卷積便可完成， 便于集成化能譜測量系統(tǒng)中機器語言的實現(xiàn)， 在工業(yè)應(yīng)用中具有實用性； (3)受統(tǒng)計漲落的影響較小且分解結(jié)果較為精確。