基于壓縮感知理論的非線性γ譜分析方法*

馮丙辰 方 晟 張立國 李 紅 童節(jié)娟 李文茜

(清華大學(xué)核能與新能源技術(shù)研究院，北京 100084)

(2013年1月23日收到;2013年2月26日收到修改稿)

1 引言

γ譜分析方法通過測量和分析不同能量的γ射線，獲得被測對象中放射性核素的種類及含量信息，是一種重要的放射性核素定量分析方法，在核科學(xué)與技術(shù)、環(huán)境監(jiān)測、天體物理和地質(zhì)學(xué)等領(lǐng)域有著廣泛的應(yīng)用[1-3].然而，對于放射性活度微弱的放射性核素，其特征峰有可能被γ譜的本底統(tǒng)計漲落所淹沒.同時，由于探測器分辨率的限制以及探測器電荷收集不完全導(dǎo)致的拖尾效應(yīng)，射線能量接近的γ能峰容易疊加在一起，形成重疊峰[4].

這兩個問題對γ能譜分析的準(zhǔn)確性有著明顯的影響.傳統(tǒng)γ譜分析方法將γ譜分析過程分割成線性平滑、尋峰、峰擬合等獨立步驟分別處理上述問題，沒有納入有關(guān)物理機理，缺乏嚴(yán)密的數(shù)學(xué)框架，易造成低計數(shù)峰丟失、重疊峰分辨困難等問題.而近年來誕生的壓縮感知理論指出：對于稀疏信號，可以在降低采樣率的同時，通過特定的非線性信號恢復(fù)算法，恢復(fù)原始信號.壓縮感知理論在遙感成像、地震勘探、表面測量、多媒體數(shù)據(jù)編碼等領(lǐng)域得到了廣泛的應(yīng)用[5-11].

為了解決傳統(tǒng)γ譜分析中弱峰和重疊峰探測能力差的問題，基于壓縮感知理論框架，提出了一種新的γ譜分析方法，并通過數(shù)值模擬和蒙特卡洛模擬進行了驗證.該方法從譜儀對真實γ譜信號調(diào)制的物理機理出發(fā)，通過數(shù)學(xué)建模，將γ譜分析轉(zhuǎn)化為一個求逆問題;基于壓縮感知理論將γ譜的稀疏信號模型與譜儀調(diào)制模型耦合起來，通過非線性優(yōu)化方法求解，在降低統(tǒng)計漲落的同時，減小了能譜展寬，有效提高了弱峰和重疊峰的分辨能力.該方法為γ譜分析方法的發(fā)展提供了新的思路.

2 γ譜儀調(diào)制的物理建模

γ射線與物質(zhì)的相互作用會導(dǎo)致探測器內(nèi)的原子電離，γ譜儀通過收集電離產(chǎn)生的電荷測量樣品的γ能譜.然而由于電離產(chǎn)生電荷的統(tǒng)計漲落、探測器的邊緣效應(yīng)以及彈道效應(yīng)[12]都會造成帶電粒子能量的波動，從而導(dǎo)致γ譜產(chǎn)生展寬.實驗觀察到的γ譜展寬通常符合高斯分布，可以用如下的高斯函數(shù)進行描述：

其中a，μ，σ均為大于0的實常數(shù).

根據(jù)γ譜儀的展寬分布，可以建立下述譜儀調(diào)制模型：

其中，ρ是以不同道址計數(shù)為元素的列向量，代表實際的γ信號.d是經(jīng)過γ譜儀調(diào)制后的列向量，代表觀測到的γ信號.Ψ為描述γ譜儀對于真實γ譜調(diào)制作用的系統(tǒng)矩陣：

可以看出Ψ的列向量以對角線元素為中心呈歸一化的高斯分布，hkl的值代表第l道計數(shù)因高斯展寬對第k道計數(shù)的貢獻.

根據(jù)譜儀調(diào)制模型，在最小二乘意義下，真實γ譜數(shù)據(jù)的求解可以表示為

3 壓縮感知理論

壓縮感知理論提供了利用非線性優(yōu)化方法完整重建稀疏信號的框架，為許多逆問題的求解提供了有力的途徑[14,15].壓縮感知理論包括稀疏信號模型、數(shù)據(jù)采集和稀疏信號完整重建算法三個方面[16,17].對于γ譜分析而言，上述三個模型中獲得直接應(yīng)用的是稀疏信號模型和稀疏信號完整重建算法.

3.1 稀疏信號模型

其中‖·‖0是0范數(shù).則信號向量ρ是一個k稀疏的信號向量，Φ為對應(yīng)的稀疏變換矩陣.對于某些本身即具備稀疏性的信號，不需要進行稀疏變換.此時，Φ退化成單位矩陣，(8)式可以改寫成

3.2 稀疏信號完整重建

對于稀疏信號ρ，如果采樣矩陣Θ滿足約束等距(Restricted Isometry Property)條件，則可以在只采集少量數(shù)據(jù)d的情況下，通過求解如下的優(yōu)化問題，完整重建信號[18].

然而(10)式是NP難度的問題，直接求解較為困難.故通常用1范數(shù)代替(10)式中的0范數(shù)，從而將(10)式的非凸優(yōu)化問題轉(zhuǎn)化為更易于求解的凸優(yōu)化問題[19]：

其中，λ是正則化參數(shù).

4 基于壓縮感知的非線性γ譜分析

4.1 理論模型

由于γ譜中的特征峰僅出現(xiàn)在幾個少數(shù)的能量段上，所以γ譜是典型的稀疏信號.對照(2)式與(11)式，不難發(fā)現(xiàn)，譜儀調(diào)制模型可以作為約束條件，直接納入到壓縮感知框架中，系統(tǒng)矩陣Ψ即對應(yīng)采樣矩陣Θ.由此，即可以實現(xiàn)稀疏信號模型和譜儀調(diào)制模型的耦合，將γ譜分析轉(zhuǎn)換成以真實γ譜為求解目標(biāo)的數(shù)值最優(yōu)化問題，并通過下式求解：

由于受譜儀性能的影響，測量所得γ譜存在展寬，拖尾等現(xiàn)象，降低了γ譜的稀疏性.因此，本文選用一階差分算子對γ譜進行變換，以獲得更為稀疏的變換系數(shù)，通過這些變換系數(shù)來衡量γ譜的稀疏程度.相應(yīng)地，(13)式可以進一步寫成

其中，?是一階差分變換算子矩陣.上式等號右邊第一項衡量了γ譜與譜儀調(diào)制模型的偏差;第二項衡量了γ譜與稀疏信號模型的偏差.正則化參數(shù)λ則對兩種誤差進行均衡，以使所得到的求解結(jié)果與譜儀調(diào)制模型和稀疏信號模型都有較好的符合.

4.2 數(shù)值求解算法

其中Δt是人工演化的時間步長，滿足CFL條件，以保證人工演化過程穩(wěn)定收斂.因此，(14)式的迭代求解過程可以表示如下：

初始值ρ0=0，對于第t次迭代，有

5 模擬與驗證

為了驗證所建立的γ譜分析方法，分別進行了數(shù)值模擬和γ譜測量的蒙特卡洛模擬.

5.1 數(shù)值模擬

通過計算機模擬產(chǎn)生了長度為1024道的γ譜信號，并運用隨機抽樣方法產(chǎn)生了7個γ特征峰.如圖1(a)中短箭頭所示，所模擬的理想γ譜的第301道存在計數(shù)較低的弱特征峰.第536和549道則存在兩個十分鄰近的特征峰(圖1(a)中長箭頭所指).理想γ譜的各個特征峰的平均半高寬為6個道址.

圖1 所模擬的γ譜信號 (a)理想γ譜信號;(b)考慮了譜儀展寬和統(tǒng)計漲落的γ譜信號

為了模擬實際譜儀造成的展寬效應(yīng)，采用半高寬為15個道址的高斯卷積核對理想γ譜進行了卷積.在此基礎(chǔ)上，進一步疊加了高斯白噪聲，模擬實際γ譜中的統(tǒng)計漲落現(xiàn)象.如圖1(b)所示，理想γ譜信號經(jīng)過高斯卷積后，各個特征峰均出現(xiàn)了不同程度的展寬.其中，位于第536和549道的特征峰由于卷積的展寬作用而重疊在一起(如圖1(b)長箭頭所示);位于第301道的弱特征峰則在展寬和統(tǒng)計漲落的共同作用下，幾乎淹沒在統(tǒng)計漲落之中.

5.2 蒙特卡洛模擬

蒙特卡洛方法可以逼真地模擬光子與探測系統(tǒng)各部件物質(zhì)中的原子發(fā)生的物理作用過程(包括光電效應(yīng)、康普頓效應(yīng)和電子對效應(yīng))，并記錄光子沉積在探頭上的能量，生成與實驗符合較好的脈沖高度譜.MCNP即是一種通過蒙特卡洛方法進行粒子(包括光子，中子和電子)輸運模擬的大型軟件包，是由美國Los Alamos國家實驗室研制開發(fā)的，在射線測定領(lǐng)域有著廣泛的應(yīng)用[20,21].

為了驗證所提出γ譜分析方法，運用MCNP程序模擬了ORTEC公司的γ射線探測系統(tǒng).該系統(tǒng)包括GEM30P70型HPGe探測器、DSPEC-Plus數(shù)字化光譜儀和X-Cooler-Ⅱ電致冷器.通過對放射性燃料球的分析計算，根據(jù)探測器探頭可能探測到的γ射線強度建立了采樣過程，并運用蒙特卡洛方法對探測過程進行了模擬.

改進的模型預(yù)測直接轉(zhuǎn)矩控制方案的結(jié)構(gòu)框圖如圖2所示,與傳統(tǒng)模型預(yù)測直接轉(zhuǎn)矩控制相比,增加了延時補償環(huán)節(jié)和優(yōu)化矢量選擇器來提高控制系統(tǒng)的性能.

蒙特卡洛模擬中的幾何模型如圖2所示.燃料球發(fā)出的γ射線經(jīng)過密封部件和準(zhǔn)直系統(tǒng)，被探測器接收，探測器根據(jù)γ射線的強度和能量將接收到的信號轉(zhuǎn)化為能譜圖.所采用的蒙特卡洛方法根據(jù)放射性燃料球的γ光子發(fā)射率和準(zhǔn)直系統(tǒng)的角度因子進行源抽樣，記錄能夠通過準(zhǔn)直器系統(tǒng)到達探測器探頭的γ光子數(shù)，模擬γ射線能量沉積的幅度分布，采樣計算的γ光子總數(shù)約為104.

圖2 γ射線探測系統(tǒng)示意圖(部件1，2，3，4分別代表裝有放射性燃料球的定位器、密封部件、準(zhǔn)直儀和GEM30P70型HPGe探測器)

6 結(jié)果分析

圖3顯示了數(shù)值模擬結(jié)果.如圖3(a)所示，非線性γ譜估計方法在去除統(tǒng)計漲落的同時，有效減小了特征峰的展寬.圖3(b)放大了γ譜中含有重疊峰和低計數(shù)特征峰的能段，更清晰地顯示了非線性γ譜估計方法對特征峰信息的恢復(fù)效果.由于展寬而重疊在一起的兩個特征峰被成功分離開(圖3(b)中長箭頭所指)，且恢復(fù)后的兩個特征峰的峰形與理想γ譜中的對應(yīng)特征峰峰形基本一致，準(zhǔn)確度較高.同時，整個γ譜中的統(tǒng)計漲落被有效降低，低計數(shù)特征峰信息則保留完好，避免了傳統(tǒng)線性平滑方法中低計數(shù)峰易丟失的問題.

表1比較了數(shù)值模擬中理想γ譜、存在展寬和漲落的γ譜和非線性γ譜估計結(jié)果的峰面積和半高寬.由表可見，展寬和統(tǒng)計漲落會明顯改變γ譜的峰面積和半高寬的信息，給γ譜定量帶來明顯的誤差.非線性γ譜估計方法有效減少了這種誤差，其估計結(jié)果的峰面積和半高寬與理想γ譜符合較好.

圖4顯示了γ譜測量的蒙特卡洛模擬結(jié)果.與數(shù)值模擬結(jié)果類似，整體上，非線性γ譜估計方法有效降低了所采集γ譜中的統(tǒng)計漲落;同時，各個γ譜的特征峰信息也得到了很好地保留(圖4(a)).圖4(b)放大了蒙特卡洛模擬譜中含有重疊峰和低計數(shù)特征峰的能段，以進一步評估非線性γ譜估計方法的效果.如圖4(b)中所示，由于譜儀的展寬作用，662 keV的Cs-137特征峰、665 keV的Ce-143特征峰、670 keV和672 keV的I-132特征峰均被淹沒在了668 keV的I-132特征峰中.經(jīng)過非線性γ譜估計方法處理后，很好地將模擬測量中的五重重疊峰分解，恢復(fù)了強特征峰周圍的弱特征峰，消除了譜儀展寬帶來的γ特征峰丟失問題.同時，計數(shù)較低的弱特征峰的峰形狀也與實際采集γ譜符合地很好，這說明非線性γ譜估計方法對于弱特征峰信息破壞較少.

圖3 數(shù)值模擬結(jié)果 (a)對1024道模擬γ譜進行非線性估計的結(jié)果;(b)估計所得γ譜的局部放大圖

表1 數(shù)值模擬中的γ譜特征峰面積和半高寬計算結(jié)果

圖4 γ譜測量蒙特卡洛模擬結(jié)果 (a)非線性γ譜估計的結(jié)果;(b)非線性γ譜估計的結(jié)果的局部放大圖

對于HPGe探測器，只考慮統(tǒng)計漲落影響時，γ譜的理論最小展寬可以用下式表示：

其中F為Ge的法諾因子，取值0.13;E為特征峰中心能量;W為Ge的平均電離能，取值2.98×10-3keV.蒙特卡洛模擬中以此為標(biāo)準(zhǔn)，對非線性γ譜分析結(jié)果進行了定量評估.

表2比較了由(19)式計算出的HPGe探測器理論最小能譜展寬、蒙特卡洛模擬譜的能譜展寬和非線性γ譜估計結(jié)果的能譜展寬.可以看出，蒙特卡洛模擬譜主要的特征峰展寬約為理論最小展寬的2倍，經(jīng)過非線性γ估計方法處理后，特征峰的展寬得到了有效的降低，比模擬測量結(jié)果更接近理論最小展寬;與蒙特卡洛模擬譜相比，非線性γ譜估計結(jié)果最大減小了16.7%的特征峰展寬，最小減小了10.9%的特征峰展寬.以上結(jié)果表明，非線性γ譜估計方法能夠較為穩(wěn)定的減小譜儀調(diào)制形成的能譜展寬，更真實地反映出實際γ譜的特征.

表2 蒙特卡洛模擬γ譜特征峰展寬對比結(jié)果

7 結(jié)論

本文基于γ譜儀對于γ譜信號的物理調(diào)制過程和壓縮感知理論框架，提出了一種新的γ譜分析方法.該方法從測量譜儀中γ譜信號形成的物理過程出發(fā)，通過數(shù)學(xué)建模，將γ譜分析轉(zhuǎn)化為一個求逆問題，并建立以信號稀疏性作為約束條件的非線性優(yōu)化求解方法.數(shù)值模擬和蒙特卡洛模擬結(jié)果表明：該方法能在降低本底漲落的同時，有效減小γ譜中特征峰的展寬，從而提高弱峰和重疊峰的分辨能力，有望在環(huán)境監(jiān)測、核燃料燃耗監(jiān)測等γ譜應(yīng)用領(lǐng)域獲得進一步應(yīng)用.

[1]HaoF H，Hu G C，Liu S P，Gong J，Xiang Y C，Huang R L 2005 Acta Phys.Sin.54 3523(in Chinese)[郝樊華，胡廣春，劉素萍，龔建，向永春，黃瑞良2005物理學(xué)報54 3523]

[2]ZhaoH F，Du L，He L，BaoJ L 2011 Acta Phys.Sin.60 028501(in Chinese)[趙鴻飛，杜磊，何亮，包軍林2011物理學(xué)報60 028501]

[3]Wang C J，BaoD M，Cheng S，Zhang A L 2008 Acta Phys.Sin.57 5361(in Chinese)[王崇杰，包東敏，程松，張愛蓮2008物理學(xué)報57 5361]

[4]Grenier G，Poussier C 1970 Nuclear Instruments and Methods 89 199

[5]Sun B，Jiang J J 2011 Acta Phys.Sin.60 110701(in Chinese)[孫彪，江建軍2011物理學(xué)報60 11701]

[6]BaiX，LiY Q，ZhaoS M 2013 Acta Phys.Sin.62 044209(in Chinese)[白旭，李永強，趙生妹2013物理學(xué)報62 044209]

[7]LustigM，DonohoD，PaulyJM2007MagneticResonanceinMedicine 58 1182

[8]Yu S W，Khwaja A S，Ma J W 2012 SignalProcess.92 357

[9]Ma J W，HussainiM Y 2011 IEEE.T.Instrum.Meas.60 3128

[10]Ma J W 2010 IEEE.T.Instrum.Meas.59 1600

[11]Ma J W，Plonka G，HussainiM Y 2012 IEEE.T.Circ.Syst.Vid.22 1354

[12]Hou Q W，CaoB Y，GuoZ Y 2009 Acta Phys.Sin.58 7809(in Chinese)[侯泉文，曹炳陽，過增元2009物理學(xué)報58 7809]

[13]Briesmeister J F 1997 MCNP-A GeneralMonte CarloN-Particle Transport Code Version 4B，LA-12625-M.Los Alamos NationalLaboratory

[14]Candes E J，TaoT 2005 IEEE.T.Inform.Theory.24 118

[15]DonohoD L，Elad M，Temlyakov V N 2006 IEEE.T.Inform.Theory.52 6

[16]Candes E J，Romberg J，TaoT 2006 IEEE.T.Inform.Theory.52 489

[17]DonohoD L 2006 IEEE.T.Inform.Theory.52 1289

[18]Candes E J，Romberg J，TaoT 2006 Commun.Pur.Appl.Math.59 1207

[19]Candes E J，Wakin M B 2008 IEEE.SignalProc.Mag.25 21

[20]Xu M H，Liang T J，Zhang J 2006 Acta Phys.Sin.55 2357(in Chinese)[徐妙華，梁天驕，張杰2006物理學(xué)報55 2357]

[21]Xu H B，Peng X K，Chen C B 2010 Chin.Phys.B 19 062901