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

趙思文，吳 怡，王崇杰

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

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

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

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

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

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

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

(1)

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

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

(2)

式中,d=rank(X)是軌跡矩陣的秩，即X的非零奇異值個數，且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)對角平均化.對角平均化實質上是將降噪后的軌跡矩陣Y∈RL×K轉換成降噪后的γ能譜.設yij為軌跡矩陣Y的元素，則降噪后γ能譜為

(4)

選擇恰當的嵌入維數L和重構階數p，即能達到降低γ能譜噪聲的目的.

2 實測γ能譜降噪處理與結果

2.1 儀器設備與測量

選用EG&G ORTEC 918A-HPGe γ能譜儀，多道分析器最大道址為8 192，譜儀系統在60Co 1 332.50 keV處的能量分辨率為2.0 keV.放射源為60Co標準點源，活度為2.75 kBq.測量活時間設置為1,2,3,4,5 h，相應譜數據分別記為Co1h～Co5h.

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

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

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

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

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

(5)

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

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

圖2 信噪比隨嵌入維數L的變化規律

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

(6)

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

圖3 剩余噪聲奇異譜

2.3 γ能譜降噪與結果

分別采用線性濾波和同態濾波方式[9]，對實測60Co γ能譜進行了降噪處理，并與傳統的FFT和WT降噪方法進行比較.

線性濾波是直接對譜數據進行變換或分解，然后通過選取適當的濾波參量消除噪聲分量，最后通過逆變換或重構得到降噪γ能譜.同態濾波則用于消除乘性噪聲，首先對譜數據進行對數運算，將譜成分與噪聲之間的非線性關系轉換成線性關系，然后通過線性濾波消除噪聲分量，最后通過指數運算得到降噪譜.

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

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

(a)原始譜

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

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

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

(7)

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

(8)

(9)

(10)

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

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

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

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

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

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

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

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

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

3 結 論

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