小波閾Wiener濾波器在航空伽瑪能譜降噪中的應(yīng)用

羅耀耀 葛良全 熊 超 徐立鵬 花永濤

(成都理工大學(xué) 成都 610059)

航空γ能譜測(cè)量是一種快速、經(jīng)濟(jì)、有效的核地球物理方法[1,2]，航空γ射線能譜儀，記錄0–3 MeV的多道(256 道或更多道)全譜數(shù)據(jù)。但這些單次測(cè)量的原始能譜的統(tǒng)計(jì)漲落明顯，對(duì)這些能譜數(shù)據(jù)進(jìn)行降噪處理，能降低數(shù)據(jù)中統(tǒng)計(jì)漲落帶來的影響，是提高放射性核素含量精度的關(guān)鍵。

小波分析是近20多年發(fā)展起來的交叉學(xué)科，在信噪分離、提取弱信號(hào)等方面應(yīng)用較廣。小波去噪有硬閾值去噪和軟閾值去噪[3–6]，但這兩個(gè)函數(shù)存在一定的缺陷，如硬閾值函數(shù)的不連續(xù)性,軟閾值函數(shù)估計(jì)的小波系數(shù)與被處理信號(hào)的小波系數(shù)之間存在恒定偏差等，限制了該方法的應(yīng)用效果[7,8]。

本文分析了最小二乘法及小波的軟硬閾值去噪方法，并采用雙小波基與Wiener濾波器相結(jié)合的方法，實(shí)現(xiàn)了航空γ能譜的降噪新方法。

1 小波多分辨率分析閾值降噪方法

小波變換能同時(shí)在時(shí)域和頻域中對(duì)信號(hào)進(jìn)行分析，具有自動(dòng)變焦功能，故能有效區(qū)分信號(hào)中的突變部分和噪聲，從而實(shí)現(xiàn)信號(hào)的降噪。

航空γ能譜數(shù)據(jù)f(n)可用式(1)表示

式中，f(i)為能譜信號(hào)，s(i)為純信號(hào)，n(i)為高斯白噪聲信號(hào)，n(i)=σv(i)為加性隨機(jī)噪聲，其強(qiáng)度是σ。噪聲信號(hào)n(i)與采樣信號(hào)s(i)不相關(guān)。由于小波變換是線性的，所以經(jīng)過小波變換后的小波系數(shù)模型可寫成：

其中，y=wf，θ=ws，z=wn分別為信號(hào)、純凈信號(hào)以及噪聲的小波系數(shù)，w為小波變換。

噪聲和信號(hào)在小波變換下具有如下表現(xiàn)形式:從能量觀點(diǎn)看，噪聲能量分布在所有小波系數(shù)上，而信號(hào)能量僅分布在一小部分系數(shù)上；噪聲經(jīng)小波變換后幅值較小，數(shù)目較多，而信號(hào)經(jīng)小波變換后幅值較大，數(shù)目較少。通過分析，可以構(gòu)造閾值去噪方法，基于小波的消除噪聲的方法可分為三個(gè)步驟[8]：

(1) 選擇小波和小波分解的層次，計(jì)算信號(hào)s第1層到第N層的小波分解系數(shù)。

(2) 高頻系數(shù)的閾值選取，對(duì)于從第 1層到第N層的每一層，選取1個(gè)閾值(也可用同1個(gè)閾值)，對(duì)高頻系數(shù)用軟、硬閾值進(jìn)行處理(用硬閾值處理信號(hào)較粗糙，一般采用軟閾值)。

(3) 根據(jù)第N層的低頻系數(shù)和從第1層到第N層經(jīng)修改后的高頻系數(shù)，重構(gòu)信號(hào)。

2 改進(jìn)的經(jīng)驗(yàn)Wiener濾波器的設(shè)計(jì)

Wiener濾波器在最小均方誤差意義上是最優(yōu)的，其前提是已知信號(hào)和噪聲的統(tǒng)計(jì)特性。因此，實(shí)際應(yīng)用中需要先對(duì)信號(hào)和噪聲作相應(yīng)估計(jì)，為此，本文根據(jù)Wiener濾波器和小波變換的特點(diǎn)，采用小波域Wiener濾波器對(duì)航空γ能譜譜線數(shù)據(jù)進(jìn)行處理。

由于小波變換對(duì)信號(hào)有去相關(guān)的特性，對(duì)小波系數(shù)可采用Wiener濾波器，小波的Wiener濾波器設(shè)計(jì)值為：其中，σ2是噪聲z(i)的方差，因?yàn)棣?i)是未知的，所以，使用其估計(jì)值(i) ，于是得到經(jīng)驗(yàn)的Wiener濾波器的值為：

在設(shè)計(jì)經(jīng)驗(yàn)Wiener濾波器時(shí)，須從信號(hào)中估計(jì)出θ2(i)和σ2，其中噪聲σ2可由小波系數(shù)估計(jì)，而θ2(i)則較難估計(jì)。故采用兩個(gè)小波基W1和W2，W1對(duì)信號(hào)進(jìn)行降噪，從而得到原始信號(hào)較為準(zhǔn)確的估計(jì)(i)，再用經(jīng)驗(yàn)Wiener濾波器在W2小波閾中進(jìn)行處理，由此獲得較好的降噪效果。采用兩個(gè)小波基可保證在 W1閾內(nèi)由于閾值函數(shù)被置為零的系數(shù)在W2小波閾中得以恢復(fù)，達(dá)到降噪效果。根據(jù)要求設(shè)計(jì)的小波域經(jīng)驗(yàn)Wiener濾波器如圖1所示。

圖1 小波閾Wiener濾波器結(jié)構(gòu)圖Fig.1 Structure of the wavelet domain Wiener filter.

3 對(duì)航空γ能譜數(shù)據(jù)處理與分析

為研究該方法在航空能譜數(shù)據(jù)處理中的有效性和準(zhǔn)確性，將內(nèi)蒙某測(cè)區(qū)航空γ能譜原始譜線數(shù)據(jù)按工作流程分為早晚校正數(shù)據(jù)集和測(cè)區(qū)數(shù)據(jù)集。早晚校正數(shù)據(jù)集是指航空γ能譜儀在飛行前和飛行后對(duì)儀器穩(wěn)定性、譜線峰位等參數(shù)檢查過程中獲得的數(shù)據(jù)；測(cè)區(qū)數(shù)據(jù)是指航空γ能譜儀在測(cè)區(qū)工作時(shí)獲得的原始數(shù)據(jù)。

3.1 早晚校正數(shù)據(jù)集的處理與分析

取航空γ能譜儀位于內(nèi)蒙某機(jī)場(chǎng)地面靜止測(cè)量時(shí)的數(shù)據(jù)，即早晚校正數(shù)據(jù)進(jìn)行分析，比較了不同降噪方法后的鉀窗、鈾窗、釷窗內(nèi)的計(jì)數(shù)率均值和相應(yīng)的標(biāo)準(zhǔn)差(表1)；與原始譜線相比，采用最小二乘法和小波軟閾值濾波法及小波域 Wiener濾波器方法降噪后，各個(gè)窗口的計(jì)數(shù)率平均值基本不變，但標(biāo)準(zhǔn)差發(fā)生了變化，經(jīng)小波域Wiener濾波器方法處理后誤差小于其他兩種方法?？梢?，采用小波方法應(yīng)用于單次測(cè)量取得了比較理想的效果。

3.2 測(cè)區(qū)數(shù)據(jù)處理與分析

圖2顯示飛行高度位于142 m航空γ能譜儀獲得的原始能譜數(shù)據(jù)與經(jīng)小波域 Wiener濾波器去噪后的譜線對(duì)比圖。可以看出，能譜經(jīng)小波域Wiener濾波器后，降低了噪聲能譜數(shù)據(jù)統(tǒng)計(jì)漲落帶來的影響，特征峰均能清晰顯露，如208Tl的2.62 MeV、214Bi的1.76 MeV和1.12 MeV、40K的1.46 MeV、228Ac的0.911 MeV的特征γ射線峰以及由214Bi的0.609 MeV和208Tl的0.583 MeV構(gòu)成的重疊峰。

圖3為采用小波域Wiener濾波器和5點(diǎn)3次多項(xiàng)式擬合法去噪后得到的譜線。觀察208Tl的0.609 MeV峰和40K的1.46 MeV以及208Tl的2.62 MeV的峰可以看出，小波域Wiener濾波方法的降噪效果明顯優(yōu)于傳統(tǒng)的能譜降噪方法(多項(xiàng)式光滑方法)，糾正了傳統(tǒng)方法導(dǎo)致的譜線畸變。

表1 采用不同降噪方法處理全能譜數(shù)據(jù)前后鉀窗、鈾窗、釷窗內(nèi)的計(jì)數(shù)率均值和相應(yīng)的標(biāo)準(zhǔn)差Table 1 Counts of K, U and Th windows before and after denoising by different methods.

4 結(jié)語

本文介紹了采用小波域Wiener濾波器對(duì)航空γ能譜數(shù)據(jù)進(jìn)行降噪的方法，并應(yīng)用于某測(cè)區(qū)的數(shù)據(jù)中，與常規(guī)方法進(jìn)行了對(duì)比分析。通過對(duì)航空γ能譜儀在測(cè)區(qū)工作時(shí)早晚校正數(shù)據(jù)三個(gè)窗口數(shù)據(jù)的平均值和標(biāo)準(zhǔn)差分析，對(duì)窗口數(shù)據(jù)的標(biāo)準(zhǔn)差有所改善；航空γ能譜儀在高空中(142 m)的譜線數(shù)據(jù)表明，該方法能有效降低航空γ能譜數(shù)據(jù)的統(tǒng)計(jì)漲落，凸顯能譜的特征峰。

圖2 采用小波域Wiener濾波器降噪前(*)后(●)譜線比較Fig.2 Energy spectra before(*) and after (●)denoising.

圖3 采用多項(xiàng)式擬合法(*)和小波域Wiener濾波器(▲)去噪后的譜線比較Fig.3 Denoising result comparison between conventional and wavelet domain Wiener filter.* Smoothed by Polynomial fitting, ▲ Smoothed by wavelet domain Wiener filter

1 Bruce L D. Recent advances in aerial gamma-ray surveying[J]. Journal of Environmental Radioactivity,2004, 76(1–2): 225–236

2 International atomic energy agency[S]. Guidelines for radioelement mapping using gamma-ray spectrometry data, 2003: 52–71

3 于國梁, 張一云, 賀梁國. 小波分析法檢測(cè)譜特征峰[J].原子能科學(xué)技術(shù), 2008, 42(11): 984–987 YU Guoliang, ZHANG Yiyun, HE Liangguo. Characteristic peak detection by using wavelet analysis in γ-ray spectra[J]. Atomic Energy Science and Technology, 2008,42(11): 984–987

4 鄒永祥, 吳建平. 基于多分辨率分析的 X 熒光能譜的除噪與弱峰檢測(cè)[J]. 核技術(shù), 2008, 31(4): 260–264 ZOU Yongxiang,WU Jianping. Multi-resolution analysis for noise reduction and weak peak detection in X-ray fluorescence spectra[J]. Nucl Tech, 2008, 31(4): 260–264

5 肖 剛, 周春林, 李天柁, 等. 一種自適應(yīng)的射線能譜數(shù)據(jù)光滑方法[J]. 原子能科學(xué)技術(shù), 2001, 35(1): 69–72 XIAO Gang, ZHOU Chunlin, LI Tiantuo,et al. An adaptive method for γ-spectra smoothing[J]. Atomic Energy Science and Technology, 2001, 35(1): 69–72

6 史東生, 弟宇鳴, 周春林. 小波變換與傅里葉變換在能譜降噪處理中的比較研究[J]. 核電子學(xué)與探測(cè)技術(shù),2006, 26(6): 134–137, 107 SHI Dongsheng, DI Yuming, ZHOU Chunlin. Compara-tive study on Cenergy spectrum denoise by fourier and wavelet transforms[J]. Nuclear Electr onics &Detection Technology, 2006, 26(6): 134–137, 107

7 Mallat S. A wavelet tour of signal processing[M]. London:Academic Press, 1998

8 Donoho D L. De-noising by soft-thresholding[J]. IEEE Trans on Information Theory, 1995, 41(3): 613–627