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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

3 壓縮感知理論

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

3.1 稀疏信號(hào)模型

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

3.2 稀疏信號(hào)完整重建

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

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

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

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

4.1 理論模型

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

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

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

4.2 數(shù)值求解算法

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

初始值ρ0=0，對(duì)于第t次迭代，有

5 模擬與驗(yàn)證

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

5.1 數(shù)值模擬

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

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

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

5.2 蒙特卡洛模擬

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

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

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

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

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

6 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

7 結(jié)論

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

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

[3]Wang C J，BaoD M，Cheng S，Zhang A L 2008 Acta Phys.Sin.57 5361(in Chinese)[王崇杰，包東敏，程松，張愛(ài)蓮2008物理學(xué)報(bào)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é)報(bào)60 11701]

[6]BaiX，LiY Q，ZhaoS M 2013 Acta Phys.Sin.62 044209(in Chinese)[白旭，李永強(qiáng)，趙生妹2013物理學(xué)報(bào)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)[侯泉文，曹炳陽(yáng)，過(guò)增元2009物理學(xué)報(bào)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é)報(bào)55 2357]

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