基于奇異譜分析的γ能譜降噪算法

趙思文，吳 怡，王崇杰

(遼寧師范大學(xué) 物理與電子技術(shù)學(xué)院，遼寧 大連 116029)

γ能譜分析是放射性核素分析及核材料類型識別的重要技術(shù)手段，廣泛應(yīng)用于環(huán)境監(jiān)測、核武器核查、宇宙探測以及核安全保障等領(lǐng)域[1-6].然而，γ能譜中存在的噪聲給γ能譜分析與識別帶來較大困難，甚至導(dǎo)致錯(cuò)誤結(jié)果(特別是對于低水平放射性分析和差異甚微的核材料識別[2-4]).因此，降噪處理是γ能譜分析與識別過程中重要的技術(shù)環(huán)節(jié).

傳統(tǒng)的γ能譜降噪方法主要有道址域的曲線擬合移動平滑法和頻域的濾波法[7-8].移動平滑法容易引起譜線變形、失真，甚至畸變[8].頻域?yàn)V波法是目前較常用的γ能譜降噪方法，其中主要有基于快速傅里葉變換(Fast Fourier transform，F(xiàn)FT)和小波變換(Wavelet transform，WT)的線性濾波法[8-9].γ能譜噪聲的頻譜通常是全域的，與γ能譜成分的頻譜嚴(yán)重重疊，傳統(tǒng)的頻域?yàn)V波法難以對噪聲和譜成分進(jìn)行有效分離.因此，常常存在降噪不充分或降噪過度的現(xiàn)象.另外，除加性噪聲外，γ能譜還含有乘性噪聲，而傳統(tǒng)的降噪方法往往忽視γ能譜的乘性噪聲[9]，導(dǎo)致降噪效果不理想.

奇異譜分析(Singular spectrum analysis，SSA)是在奇異值分解理論基礎(chǔ)上興起的獨(dú)立于信號模型的無參量頻譜估計(jì)技術(shù)[10]，已成功應(yīng)用于信號的降噪處理[11-12].該技術(shù)具有不受噪聲頻譜分布影響和自適應(yīng)降噪的特點(diǎn)[12].本文首先提出了基于SSA的γ能譜降噪方法，并介紹了降噪的基本原理與具體算法.其次，通過分析γ能譜的奇異譜特征，對軌跡矩陣嵌入維數(shù)和重構(gòu)奇異值階數(shù)的最優(yōu)取值方法進(jìn)行了研究和討論.最后，對實(shí)測60CoHPGe γ能譜進(jìn)行了降噪處理，并與傳統(tǒng)的FFT和WT濾波法進(jìn)行比較，驗(yàn)證了降噪算法及其參量選取方法的有效性.

1 基于SSA的γ能譜降噪原理

基于奇異譜分析的γ能譜降噪算法的基本思想是：首先通過構(gòu)建軌跡矩陣將一維γ能譜數(shù)據(jù)映射到高維相空間；然后對軌跡矩陣進(jìn)行奇異值分解，并根據(jù)奇異譜的特征對γ能譜成分和噪聲成分進(jìn)行分離，較大的奇異值對應(yīng)γ能譜成分，較小的奇異值則對應(yīng)噪聲成分；最后利用較大奇異值成分對γ能譜進(jìn)行重構(gòu)，從而達(dá)到降噪的目的.降噪算法包括4個(gè)步驟：軌跡矩陣嵌入、奇異值分解、奇異值分組和對角平均化[10-12].

1)軌跡矩陣嵌入.設(shè)x1,x2,…，xN′為一維γ能譜數(shù)據(jù)，則相應(yīng)γ能譜的L×K軌跡矩陣為

(1)

式中,N′為γ能譜的最大道址數(shù)，L為嵌入維數(shù)，1

2)軌跡矩陣奇異值分解.對軌跡矩陣X進(jìn)行奇異值分解

(2)

式中,d=rank(X)是軌跡矩陣的秩，即X的非零奇異值個(gè)數(shù)，且d≤min(L,K).σ1，σ2，…,σd為降序排列的非零奇異值，向量ui∈RL×1和向量vi∈RK×1(i=1,2,…,d)分別為軌跡矩陣X的左奇異向量和右奇異向量.

3)奇異值分組.按軌跡矩陣X奇異值的大小，將奇異值分成2組，前r(r

(3)

4)對角平均化.對角平均化實(shí)質(zhì)上是將降噪后的軌跡矩陣Y∈RL×K轉(zhuǎn)換成降噪后的γ能譜.設(shè)yij為軌跡矩陣Y的元素，則降噪后γ能譜為

(4)

選擇恰當(dāng)?shù)那度刖S數(shù)L和重構(gòu)階數(shù)p，即能達(dá)到降低γ能譜噪聲的目的.

2 實(shí)測γ能譜降噪處理與結(jié)果

2.1 儀器設(shè)備與測量

選用EG&G ORTEC 918A-HPGe γ能譜儀，多道分析器最大道址為8 192，譜儀系統(tǒng)在60Co 1 332.50 keV處的能量分辨率為2.0 keV.放射源為60Co標(biāo)準(zhǔn)點(diǎn)源，活度為2.75 kBq.測量活時(shí)間設(shè)置為1,2,3,4,5 h，相應(yīng)譜數(shù)據(jù)分別記為Co1h～Co5h.

2.2 γ能譜的奇異譜特征及其降噪算法參量的選取

為了達(dá)到最佳的降噪效果，需選擇最優(yōu)的軌跡矩陣嵌入維數(shù)L和重構(gòu)階數(shù)p.關(guān)于軌跡矩陣嵌入維數(shù)和重構(gòu)階數(shù)的取值問題，目前還沒有通用的準(zhǔn)則和方法，相關(guān)文獻(xiàn)針對不同信號推薦的方法，都具有一定的主觀性[12-15].鑒于此，本文對γ能譜的奇異譜特征進(jìn)行了分析，并研究了γ能譜降噪信噪比隨嵌入維數(shù)的變化規(guī)律.根據(jù)γ能譜降噪剩余噪聲奇異譜的特征，給出了最優(yōu)重構(gòu)階數(shù)的選取方法.

圖1是不同嵌入維數(shù)情況下Co1h γ能譜的奇異譜.由圖1可知，對于不同的嵌入維數(shù)L，奇異譜均存在1個(gè)拐點(diǎn)，拐點(diǎn)左側(cè)的奇異值較大，而且隨著階數(shù)的增加奇異值快速衰減，對應(yīng)γ能譜成分；而拐點(diǎn)右側(cè)的奇異值較小，且隨階數(shù)的變化相對平緩，對應(yīng)噪聲成分.因此，利用奇異譜中的拐點(diǎn)可對譜成分和噪聲進(jìn)行分離，進(jìn)而通過重構(gòu)得到降噪γ能譜.

(a)不同嵌入維數(shù)下的奇異譜

圖2是Co1h γ能譜降噪后信噪比隨軌跡矩陣嵌入維數(shù)L的變化規(guī)律,其中信噪比為

(5)

式中，xi是第i個(gè)道址上的譜計(jì)數(shù)，NCH1和NCH2分別是60Co 1 332.50 keV γ射線康普頓平臺左右邊界能量1 040 keV和1 090 keV所對應(yīng)的多道分析器道址[16-17].

由圖2可知，當(dāng)嵌入維數(shù)較小時(shí)，信噪比相對較低；當(dāng)L>50時(shí)，信噪比變化相對平穩(wěn)；當(dāng)L>110時(shí)，信噪比則開始緩慢下降.研究表明:當(dāng)L在50～110區(qū)間內(nèi)取值時(shí)，均可獲得較大的信噪比和較好的降噪效果，該區(qū)間是嵌入維數(shù)的最佳取值區(qū)間，且與譜數(shù)據(jù)長度基本無關(guān).

圖2 信噪比隨嵌入維數(shù)L的變化規(guī)律

由SSA降噪原理可知，剩余噪聲中不含譜成分時(shí)所對應(yīng)的最小重構(gòu)階數(shù)就是最優(yōu)重構(gòu)階數(shù).上述研究表明，噪聲奇異譜與譜成分奇異譜之間的差異十分明顯.因此，為確定最優(yōu)重構(gòu)階數(shù)，首先在γ能譜奇異譜中分別取拐點(diǎn)r和r-1右側(cè)奇異譜重構(gòu)剩余噪聲譜；然后按下式計(jì)算2個(gè)剩余噪聲譜奇異譜之間的相關(guān)系數(shù)R[18-19]：

(6)

式中，np和np-1分別表示重構(gòu)階數(shù)為p和p-1時(shí)剩余噪聲的奇異譜，D為剩余噪聲軌跡矩陣的秩；最后，逐步減小p值，當(dāng)相關(guān)系數(shù)R明顯減小，且小于一定閾值時(shí)，相應(yīng)的p值就是最優(yōu)重構(gòu)階數(shù).由式(6)計(jì)算可得，當(dāng)p<11時(shí)，p和p-1所對應(yīng)剩余噪聲奇異譜之間的相關(guān)系數(shù)R減小，而當(dāng)p≥11時(shí)，R顯著增大，且均大于0.92.圖3是嵌入維數(shù)L=50，重構(gòu)階數(shù)分別為p=11和p=12時(shí)Co1h的剩余噪聲奇異譜，二者相關(guān)系數(shù)R=0.59.因此，p=11即是最優(yōu)重構(gòu)階數(shù).

圖3 剩余噪聲奇異譜

2.3 γ能譜降噪與結(jié)果

分別采用線性濾波和同態(tài)濾波方式[9]，對實(shí)測60Co γ能譜進(jìn)行了降噪處理，并與傳統(tǒng)的FFT和WT降噪方法進(jìn)行比較.

線性濾波是直接對譜數(shù)據(jù)進(jìn)行變換或分解，然后通過選取適當(dāng)?shù)臑V波參量消除噪聲分量，最后通過逆變換或重構(gòu)得到降噪γ能譜.同態(tài)濾波則用于消除乘性噪聲，首先對譜數(shù)據(jù)進(jìn)行對數(shù)運(yùn)算，將譜成分與噪聲之間的非線性關(guān)系轉(zhuǎn)換成線性關(guān)系，然后通過線性濾波消除噪聲分量，最后通過指數(shù)運(yùn)算得到降噪譜.

在SSA方法中，取軌跡矩陣嵌入維數(shù)L=50，利用上述剩余噪聲奇異譜相關(guān)系數(shù)法確定重構(gòu)階數(shù).在FFT濾波中，根據(jù)γ能譜幅頻特性，將幅值小于一定閾值的傅里葉變換系數(shù)置為0，再進(jìn)行逆變換得到降噪譜[16].在WT濾波中，選擇對稱性和緊支性較好的光滑小波“sym8”為小波基函數(shù)，采用小波分解和低頻系數(shù)進(jìn)行重構(gòu)，實(shí)現(xiàn)降噪[7-9].在確保譜線不變形的情況下，盡量提高降噪幅值閾值和小波分解層數(shù).在線性濾波和同態(tài)濾波中，小波分解層數(shù)分別為2和3.

降噪前后的Co1h譜圖如圖4所示.由圖4可看出，SSA法的降噪效果顯著，且采用同態(tài)濾波降噪效果明顯優(yōu)于線性濾波.

(a)原始譜

Co1h～Co5h譜降噪前后信噪比見表1.除Co1h譜的線性降噪結(jié)果外，SSA法降噪后信噪比均大于傳統(tǒng)降噪方法.采用線性濾波，SSA法降噪后的信噪比提高了20.64%～25.58%，整體高于傳統(tǒng)的FFT法(17.89%～26.50%)和WT法(17.83%～26.49%).采用同態(tài)濾波，SSA法降噪后的信噪比提高了31.57%～39.16%，高于FFT法(21.08%～25.96%)和WT法(28.23%～35.50%)，獲得最佳降噪效果.而采用同態(tài)濾波降噪信噪比均大于相應(yīng)的線性濾波降噪信噪比.研究結(jié)果表明，SSA法能夠降低加性噪聲，通過同態(tài)濾波方式也降低了乘性噪聲，使γ能譜信噪比得到較大幅度提高，較傳統(tǒng)頻域?yàn)V波法具有更強(qiáng)的降噪能力.

表1 降噪前后譜信噪比及其誤差

為了進(jìn)一步探究SSA降噪法對γ能譜的影響，計(jì)算了降噪前后60Co譜中1 332.50 keV γ射線特征峰(全能峰)參量.所用峰形函數(shù)為[17]

(7)

式中，N為道址，H為峰高，Np為峰位，a和b分別為峰本底系數(shù).半高全寬、凈峰面積及峰康比分別為

(8)

(9)

(10)

式中，NCH1和NCH2分別是1 040～1 090 keV能量區(qū)間的左右邊界道址.

降噪前后譜峰高及其誤差如表2所示.由表2可知，在線性濾波方式下，2種傳統(tǒng)頻域?yàn)V波法對峰高無明顯影響，而SSA降噪法對峰高幾乎無影響.在同態(tài)濾波方式下，F(xiàn)FT法降噪后峰高亦無明顯變化，而SSA法降噪后峰高略有增加，WT法降噪后的峰高明顯提高.研究表明:FFT屬于時(shí)-頻全局變換，而γ能譜中噪聲頻譜與譜成分頻譜存在交疊，頻域?yàn)V波過程中容易造成譜成分損失，從而使峰高降低.WT屬于時(shí)-頻局域變換，在同態(tài)濾波方式下能更有效地降低乘性噪聲而使峰高得到較大幅度提高[9]，但降噪效果與小波函數(shù)的對稱性、連續(xù)性等密切相關(guān)，如果小波函數(shù)選擇不當(dāng)，容易引起譜線局部畸變.SSA降噪法則是在相空間通過對譜數(shù)據(jù)的奇異值分解和重構(gòu)實(shí)現(xiàn)降噪，而且噪聲奇異譜與譜成分奇異譜之間的差異明顯，最優(yōu)重構(gòu)階數(shù)的確定不受噪聲頻譜的影響.因此，不會造成過度降噪而使峰高降低和譜線變形等現(xiàn)象，并且在同態(tài)濾波方式下，也可以更有效地分離和降低乘性噪聲，從而使峰高有所增加.

降噪前后譜半高全寬及其誤差如表3所示.由表3可知，采用線性濾波，SSA，F(xiàn)FT和WT 3種降噪法對全能峰半高全寬的影響均很小，尤其是SSA降噪法對半高全寬幾乎沒有影響.采用同態(tài)濾波，F(xiàn)FT法對半高全寬的影響亦很小，而SSA法降噪后半高全寬略有減小，WT法降噪后半高全寬則顯著減小.結(jié)合表2可知，在同態(tài)濾波方式下，SSA法與WT法在使峰高增大的同時(shí)使半高全寬減小，而且峰高與半高全寬的相對變化幅度基本相同.

表2 降噪前后譜峰高及其誤差

表3 降噪前后譜半高全寬及其誤差

降噪前后譜峰康比及其誤差如表4所示.由表4可知，在線性濾波方式下，與傳統(tǒng)頻域?yàn)V波降噪方法相同，SSA法降噪后峰康比基本保持不變.在同態(tài)濾波方式下，F(xiàn)FT法降噪后的峰康比亦保持不變，而SSA法降噪后的峰康比有小幅提高，WT法降噪后的峰康比提高幅度相對較大.結(jié)合表2可知，如果考慮峰本底，峰康比的變化反映了峰高的變化，而康普頓平臺沒有變化.結(jié)果表明：與傳統(tǒng)頻域?yàn)V波方法相同，SSA法降噪不影響譜線整體分布.

表4 降噪前后譜峰康比及其誤差

降噪前后峰凈面積及其誤差如表5所示.由表5可知，與傳統(tǒng)頻域?yàn)V波法相同，無論是線性濾波方式還是同態(tài)濾波方式，SSA法降噪后的凈峰面積保持不變.說明SSA降噪法能夠適用于后續(xù)的γ能譜定量分析，進(jìn)一步表明SSA降噪法是有效的.

表5 降噪前后峰凈面積及其誤差

3 結(jié) 論

SSA降噪法可通過軌跡矩陣將γ能譜映射到相空間，進(jìn)而在相空間通過奇異值分解和重構(gòu)實(shí)現(xiàn)信噪分離，克服了傳統(tǒng)頻域?yàn)V波方法無法解決的噪聲頻譜與譜成分頻譜的重疊問題.同時(shí)，由于軌跡矩陣只與待降噪γ能譜有關(guān)，因此，SSA降噪法不易引起譜線變形或畸變.剩余噪聲奇異譜相關(guān)系數(shù)法可確定降噪譜的最優(yōu)重構(gòu)階數(shù)，有效避免了傳統(tǒng)方法降噪過程中存在的降噪不充分或降噪過度現(xiàn)象.由于受統(tǒng)計(jì)漲落的影響，原始譜的信噪比、半高全寬及峰康比等指標(biāo)的誤差相對較大，而降噪后誤差明顯減小，反映了不同降噪方法的降噪效果.SSA降噪法對全能峰峰高幾乎沒有影響，因此，降噪不會影響γ能譜的能量線性.綜上，SSA降噪法算法簡單，具有較好的降噪性能，在采用線性濾波和同態(tài)濾波2種濾波方式下均可獲得令人滿意的降噪效果，是有效的γ能譜降噪方法.