王青山, 王冬陽(yáng), 張雄杰, 湯 彬, 吳和喜
核技術(shù)應(yīng)用教育部工程研究中心(東華理工大學(xué)), 江西 南昌 330013
γ能譜測(cè)量與分析是放射性信息獲取的重要方式, 在輻射監(jiān)測(cè)、 放射醫(yī)療、 巖性分析等領(lǐng)域[1]都有廣泛的應(yīng)用。 實(shí)際測(cè)量中γ能譜通常會(huì)受到環(huán)境本底與測(cè)量干擾等不可抗力因素的影響, 通常為多種能量和強(qiáng)度不同的γ射線(xiàn)單能譜的疊加, 特別是能量相近的γ射線(xiàn), 往往以重疊峰的形式出現(xiàn), 且強(qiáng)度弱的γ譜線(xiàn)又容易被強(qiáng)峰或本底所掩蓋, 容易對(duì)放射性核素種類(lèi)的判斷及含量的計(jì)算造成不利影響。
重疊峰的數(shù)值解析技術(shù)至關(guān)重要, 近年來(lái), 國(guó)內(nèi)外在重疊峰的分解方面進(jìn)行了較為深入的研究, 如導(dǎo)數(shù)法[2]、 正交投影法[3]、 小波分析法[4]、 遺傳算法[5-6]、 人工神經(jīng)網(wǎng)絡(luò)模型[6]、 純?cè)刈V剝離[7]、 粒子群算法[8]等方法。 都月[9]等改進(jìn)了縱向和平分迭代法, 提出優(yōu)化迭代重疊峰分解方法, 較傳統(tǒng)迭代方法其效率提高了66.7%; 黃洪全等[10]利用高斯混合模型結(jié)合期望最大化迭代算法進(jìn)行γ譜重疊峰分解; 周世融等[11]提出了峰銳化法結(jié)合雙樹(shù)復(fù)小波變換的低分離度重疊峰解析方法, 并實(shí)現(xiàn)了分離度為0.34的X熒光光譜重疊峰分解。
目前, 重疊峰的分解方法眾多且較為成熟, 但是普遍存在卷積過(guò)程復(fù)雜、 計(jì)算資源消耗大、 結(jié)果受算法參數(shù)影響較大的弊端, 不適用于測(cè)量現(xiàn)場(chǎng)能譜快速解析的實(shí)際應(yīng)用。 本文提出一種快速、 低耗的分解重疊峰的新方法, 為放射性能譜測(cè)量中重疊峰的現(xiàn)場(chǎng)分解提供支持。
銳化算法常被應(yīng)用于圖像處理領(lǐng)域中以提高圖像的對(duì)比度或清晰度, 也可在信號(hào)處理領(lǐng)域中用以人為地提高峰的可視分辨率。 在能譜分析中, 對(duì)重疊峰進(jìn)行銳化處理可人為地增大其分離度, 達(dá)到分離重疊峰的目的, 通常采用二階導(dǎo)數(shù)或四階導(dǎo)數(shù)來(lái)確定峰位置[12], 其算法可表示如式(1)
R(x)=s(x)-ks(2)(x)
(1)
式(1)中,R(x)為銳化后的信號(hào),s(x)為原始信號(hào),s(2)(x)為其二階導(dǎo)數(shù),k為銳化因子, 可平衡信號(hào)信噪比和基線(xiàn)平坦度從而調(diào)整其大小達(dá)到最佳銳化效果。
由于核衰變及測(cè)量的統(tǒng)計(jì)性, 能譜數(shù)據(jù)普遍受到統(tǒng)計(jì)漲落的影響, 一般通過(guò)譜光滑處理以增加定性及定量分析的準(zhǔn)確性。 與傅里葉變換法的“頻域?yàn)V波”不同, 褶積滑動(dòng)變換法為“時(shí)域?yàn)V波”, 即為變換函數(shù)與實(shí)驗(yàn)譜線(xiàn)進(jìn)行卷積變換[13]
(2)
式(2)中,F(xiàn)(x)為平滑后的譜線(xiàn),f(x)為原始譜線(xiàn),g(x)為變換函數(shù),m為滑動(dòng)變換的半窗寬。 隨著x的增加,g(x)將逐道向前滑動(dòng),f(x)經(jīng)過(guò)寬度為2m+1的變換從而轉(zhuǎn)換為新譜線(xiàn)F(x), 當(dāng)變換函數(shù)與峰形狀函數(shù)相同且半寬度一致時(shí), 光滑效果最優(yōu)。
銳化的本質(zhì)為對(duì)信號(hào)的微分處理, 當(dāng)一個(gè)噪聲信號(hào)經(jīng)過(guò)微分后, 信號(hào)的噪聲同樣經(jīng)過(guò)了微分。 因此, 不可利用傳統(tǒng)的銳化方法直接對(duì)含噪的信號(hào)進(jìn)行運(yùn)算。 若對(duì)譜線(xiàn)平滑進(jìn)行處理可以提供更高的分辨率, 根據(jù)卷積的性質(zhì)對(duì)式(2)變形, 信號(hào)F(x)的n階導(dǎo)數(shù)可表示為
F(n)(x)=f(n)(x)*g(x)=f(x)*g(n)(x)
(3)
由式(3)可以看出: 對(duì)原信號(hào)而言, 平滑與微分操作二者的先后順序并不影響其計(jì)算結(jié)果。 為了取得平滑后信號(hào)的微分, 我們可先對(duì)變換函數(shù)進(jìn)行微分得到其平滑算子, 然后利用該平滑算子作為濾波器對(duì)原信號(hào)進(jìn)行卷積。 若g(n)(x)選擇高斯導(dǎo)數(shù)小波, 則信號(hào)F(x)的n階導(dǎo)數(shù)便可通過(guò)小波變換實(shí)現(xiàn), 高斯小波變換的本質(zhì)和信號(hào)f(x)經(jīng)過(guò)g(x)平滑的n階導(dǎo)數(shù)運(yùn)算類(lèi)似。
將銳化后的高斯函數(shù)作為高斯銳化算子g(x)可得
(4)
其中,A為高斯銳化算子的峰值,μ為高斯銳化算子的峰位,σ為高斯銳化算子的標(biāo)準(zhǔn)差,k為高斯銳化因子。 為使信號(hào)卷積變換前后峰面積保持不變, 高斯銳化算子g(x)作為平滑窗函數(shù)還需滿(mǎn)足如式(5)關(guān)系
(5)
式(5)中,m為半窗寬。 對(duì)高斯銳化算子按式(5)進(jìn)行歸一化處理后, 分別取不同的k和σ值, 以探究算子參數(shù)對(duì)高斯銳化成形效果的影響, 其結(jié)果如圖1所示。k值越大, 其峰值越大, 峰寬基本保持不變, 且始終過(guò)定點(diǎn)(μ±σ,Ae-1/2), 峰兩側(cè)旁瓣的凹陷也越低, 且當(dāng)k
一葉一滴汗,一果一情懷。高品質(zhì)需要新技術(shù)的支撐,多年來(lái),農(nóng)場(chǎng)不斷推廣新技術(shù),為特色、高效、生態(tài)的芒果產(chǎn)業(yè)夯實(shí)了根基。
圖1 高斯銳化算子(A=1, μ=0, m=200)
結(jié)合式(2)與式(4)可得高斯銳化法(Gauss sharpening method, GSM)遞推公式為
(6)
本方法處理噪聲信號(hào)具有一定的優(yōu)越性, 既能過(guò)濾銳化所引入的噪聲又能分辨峰位置、 突出峰形, 實(shí)際應(yīng)用中可根據(jù)實(shí)際情況選擇進(jìn)行一次或多次高斯銳化成形。
由于核輻射測(cè)量數(shù)據(jù)服從高斯分布, 采用多個(gè)高斯函數(shù)的疊加以仿真能譜重疊峰
(7)
式(7)中,n為重疊的高斯峰個(gè)數(shù),Ai為高斯峰的高度,μi為高斯峰的峰位,σi為高斯峰的標(biāo)準(zhǔn)差, noise為疊加的噪聲信號(hào), 以仿真譜線(xiàn)中的統(tǒng)計(jì)漲落影響。
重疊峰中不同組分的峰高、 峰寬和峰位置參數(shù)將影響峰的重疊程度, 通常采用分離度[11]來(lái)判斷相鄰兩峰的分離程度
(8)
式(8)中,μ1和μ2分別為兩峰的峰位置,σ1和σ2分別為兩峰的峰寬。 即使沒(méi)有噪聲的信號(hào), 當(dāng)分離度過(guò)小時(shí), 重疊的兩個(gè)峰也很難分辨, 本文實(shí)驗(yàn)所用分離度均取0.3以上, 并模擬一組重疊峰(如圖2所示)
(9)
其分離度為0.55, 并利用MATLAB中“awgn”函數(shù)添加60 dB高斯白噪聲。 由于采用的信號(hào)是已知特征參數(shù)的峰信號(hào), 因此可通過(guò)比較處理結(jié)果與原特征參數(shù)來(lái)衡量峰分辨效果。
圖2 重疊峰示意圖(SNR=60 dB)
根據(jù)傳統(tǒng)銳化成形式(1)和高斯銳化成形式(6), 分別對(duì)式(9)所示的重疊峰進(jìn)行處理, 其結(jié)果如圖3所示。 可明顯看出經(jīng)過(guò)銳化的信號(hào)其峰寬變窄、 峰幅值變大、 分離度變大、 峰形狀更加明顯且獨(dú)立。 傳統(tǒng)銳化成形后的信號(hào)分辨率明顯下降, 其信噪比僅為19.5 dB, 而高斯銳化法成形結(jié)果的分辨率有了明顯的提高, 其信噪比約為52 dB。
圖3 重疊峰兩種銳化結(jié)果對(duì)比
高斯銳化法具有較好的噪聲抑制性能, 重疊峰分離效果顯著, 且成形前后峰面積的標(biāo)準(zhǔn)差約為0.008 4, 可認(rèn)為其成形前后峰面積保持不變, 為能譜的定量分析提供實(shí)驗(yàn)結(jié)果支撐。
對(duì)給定數(shù)據(jù), 使得擬合函數(shù)與其誤差平方和最小化原則, 尋求與給定點(diǎn)的距離平方和最小的曲線(xiàn)h(x), 經(jīng)過(guò)GSM成形的信號(hào)服從式(7)分布, 經(jīng)過(guò)化簡(jiǎn)變形可得
(10)
式(10)中,AF,μF,kF和σF均為擬合參數(shù), 本文使用MATLAB中l(wèi)sqcurvefit工具, 以式(10)作為期望函數(shù)進(jìn)行非線(xiàn)性擬合, 并對(duì)其進(jìn)行峰位、 峰面積計(jì)算, 以驗(yàn)證高斯銳化算法的有效及準(zhǔn)確性。
在非線(xiàn)性擬合中初始值的選擇非常重要, 若選擇不當(dāng), 則計(jì)算結(jié)果不一定收斂。 為了避免了曲線(xiàn)擬合的不唯一性, 提高曲線(xiàn)擬合結(jié)果可信度。 本文以各個(gè)子峰峰位為中心, 取7個(gè)點(diǎn)或9個(gè)點(diǎn)進(jìn)行擬合, 得到各個(gè)子峰的初始參數(shù), 再將其作為曲線(xiàn)擬合的初始值, 擬合結(jié)果如圖4所示。 模擬的標(biāo)準(zhǔn)高斯重疊峰均實(shí)現(xiàn)完全分解, 且峰位和峰面積的相對(duì)誤差分別在0.006%和0.03%以?xún)?nèi)。 結(jié)果表明, 峰銳化法結(jié)合褶積滑動(dòng)變換可快速且較準(zhǔn)確地分解標(biāo)準(zhǔn)高斯重疊峰。
圖4 重疊峰GSM非線(xiàn)性擬合結(jié)果
表1 標(biāo)準(zhǔn)高斯重疊峰模型的GSM峰位分解結(jié)果/%
表2 標(biāo)準(zhǔn)高斯重疊峰模型的GSM峰面積分解結(jié)果/%
利用高斯銳化法可對(duì)能譜段進(jìn)行重疊峰分解并提取相關(guān)參數(shù)。 對(duì)于高分辨率(SNR≥50 dB)譜線(xiàn), 即使分離度為0.375時(shí)也可以達(dá)到優(yōu)良的分解效果, 誤差可控制在1%以?xún)?nèi); 對(duì)于低分辨率譜線(xiàn)(SNR≤30 dB), 其誤差均高于5%, 且在分離度為0.375時(shí)已經(jīng)完全失真。 結(jié)合上述分析, 利用該方法可對(duì)SNR≥40dB,RS≥0.375的譜線(xiàn)實(shí)現(xiàn)較好的分離效果。
根據(jù)國(guó)際原子能機(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)混合放射源中包含上述核素時(shí), 特別是在使用NaI(Tl)閃爍探測(cè)器等低分辨率的探測(cè)器時(shí), 能譜譜峰將嚴(yán)重重疊。
為較好地仿真實(shí)驗(yàn)條件, 利用MCNP對(duì)131I,137Cs,214Bi,206Bi和26Al組合放射源進(jìn)行模擬, 根據(jù)實(shí)驗(yàn)時(shí)所用儀器的探頭, 模型中采用φ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殼包裹, 探頭外圍空間用空氣填滿(mǎn), 其結(jié)構(gòu)如圖5所示。
圖5 探測(cè)器模型
利用F8計(jì)數(shù)卡計(jì)算放射源的γ射線(xiàn)在NaI(Tl )晶體中的脈沖高度能譜分布, 所測(cè)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的混合源模擬γ能譜
該混合源γ能譜圖存在兩組重疊峰, 且完全無(wú)法識(shí)別出137Cs核素全能峰的位置。 對(duì)原始譜線(xiàn)直接進(jìn)行GSM處理可快速進(jìn)行全譜分析并對(duì)其中重疊峰進(jìn)行分解, 無(wú)需進(jìn)行能譜平滑以及扣除本底操作。 GSM全譜分解的結(jié)果如圖7所示, 分解誤差如表3所示。
圖7 混合源模擬γ能譜的GSM分解結(jié)果
表3 模擬γ能譜的GSM數(shù)值解析結(jié)果
對(duì)于單峰而言, 處理前后峰位近似不變; 對(duì)于重疊峰而言, 分解前后誤差在1~3道之內(nèi), 二者相對(duì)誤差均在1%以?xún)?nèi)。 實(shí)驗(yàn)結(jié)果證明, 利用高斯銳化法能夠快速、 高效地對(duì)γ能譜進(jìn)行全譜解析且重疊峰分解結(jié)果較為精確。
提出一種高斯銳化法快速全譜分析并分解重疊峰的新方法。 仿真實(shí)驗(yàn)結(jié)果表明, 該方法具備較好的噪聲抑制性能, 且能對(duì)SNR≥40 dB,RS≥0.375的重疊峰有較好的分解效果。 同時(shí), 利用該方法對(duì)MCNP模擬的混合源γ能譜進(jìn)行全譜解析, 實(shí)現(xiàn)了放射性核素全能峰的有效分解。
該方法的特點(diǎn)和優(yōu)勢(shì)為: (1)無(wú)需進(jìn)行能譜平滑及扣除本底操作, 即使在康普頓坪的影響下也可得到較準(zhǔn)確的結(jié)果; (2)計(jì)算資源消耗低, 能譜變換僅需一次滑動(dòng)卷積便可完成, 便于集成化能譜測(cè)量系統(tǒng)中機(jī)器語(yǔ)言的實(shí)現(xiàn), 在工業(yè)應(yīng)用中具有實(shí)用性; (3)受統(tǒng)計(jì)漲落的影響較小且分解結(jié)果較為精確。