基于數(shù)字信號(hào)處理的γ儀器譜合成方法研究

殷經(jīng)鵬，宋衛(wèi)杰，賴萬昌，馮天成

（1.西北核技術(shù)研究所，陜西 西安 710024；2.成都理工大學(xué)，四川 成都 610059）

單個(gè)探測(cè)器通常由于工藝或價(jià)格等原因不能做得太大，使單個(gè)探測(cè)器的探測(cè)效率有限。將多個(gè)探測(cè)器組成一個(gè)陣列可有效提高整個(gè)譜儀系統(tǒng)的探測(cè)效率［1－3］。例如，在車載放射性測(cè)量、航空放射性測(cè)量和肺部計(jì)數(shù)器測(cè)量系統(tǒng)中，單個(gè)探測(cè)器的探測(cè)效率往往達(dá)不到所需的要求，一般將多個(gè)NaI（Tl）探測(cè)器或多個(gè) HPGe探測(cè)器組成一個(gè)探測(cè)器陣列，即探測(cè)器組，可成倍提高探測(cè)器的有效體積，從而提高整個(gè)譜儀系統(tǒng)的探測(cè)效率，節(jié)省測(cè)量時(shí)間，提高工作效率。由于核輻射測(cè)量的相對(duì)性和隨機(jī)性，兩套γ譜儀在相同測(cè)量條件下，測(cè)量的儀器譜會(huì)有差異［4－6］，其原因是：1）γ射線與探測(cè)器相互作用的過程具有隨機(jī)性；2）由于探測(cè)器和電子線路單元性能指標(biāo)的差異，造成每路γ譜儀測(cè)量的儀器譜存在差異［7－8］。

在相同的測(cè)量條件下，由不同探測(cè)器測(cè)得的儀器譜主要差別［4－5，9－11］為：1）有效能量起始位置不同；2）相同能量的射線在不同儀器譜中產(chǎn)生全能峰對(duì)應(yīng)的道址不一致，相同的脈沖幅度經(jīng)前放及主放后，由于不同放大器的參數(shù)很難達(dá)到一致，造成相同能量的射線產(chǎn)生的脈沖波形及幅度有一定差異。以上差異使探測(cè)器陣列的儀器譜合成變得困難，因此，一般采用參數(shù)盡可能一致的探測(cè)器及其電子線路系統(tǒng)進(jìn)行基于硬件方法的儀器譜合成，或采用混合放大器進(jìn)行多路探測(cè)器 信號(hào)匹配合成［1－2，11］，但實(shí)現(xiàn)較為困難，不同譜儀的下閾值和主放大器的參數(shù)很難調(diào)節(jié)一致，同時(shí)，在高計(jì)數(shù)率情況下采用混合放大器時(shí)，造成儀器的死時(shí)間變大，不利于測(cè)量［11］。

本文提出一種將多個(gè)探測(cè)器獨(dú)立測(cè)得的、有較大差異的儀器譜通過數(shù)字信號(hào)處理的方法使其在能量坐標(biāo)上重合，實(shí)現(xiàn)儀器譜的合成。

1 基于數(shù)字信號(hào)處理的儀器譜合成法

將探測(cè)器組成一陣列，每個(gè)探測(cè)器擁有獨(dú)立的前放、主放及多道脈沖幅度分析器，測(cè)量時(shí)每路探 測(cè) 器 獨(dú) 立 獲 取 儀 器 譜［3，10－11］，將 其 通 過數(shù)字信號(hào)處理算法進(jìn)行儀器譜合成，得到合成譜，其基本原理如圖1所示。

圖1 陣列γ譜儀系統(tǒng)及儀器譜合成原理Fig.1 Principle of arrayγray spectrum analyzer system and instrument spectrum synthesis

1.1 譜合成的數(shù)學(xué)模型

對(duì)給定隨時(shí)間變化的離散信號(hào)x［n］，改變其抽樣率：一是L倍增抽樣，即把數(shù)字信號(hào)的抽樣率增加L倍，如原抽樣率為fs1，則L倍增抽樣后的抽樣率為fs2＝Lfs1，增抽樣變換原理如圖2所示；二是M倍減抽樣，即離散時(shí)間變量n變?yōu)镸n，且限定M為正整數(shù)，即把數(shù)字信號(hào)的抽樣率減少M(fèi)倍，為增抽樣的逆過程［10－12］。

對(duì)于離散信號(hào)x［n］，若希望將抽樣率轉(zhuǎn)變L／M倍，由于抽取使x［n］的數(shù)據(jù)點(diǎn)減少，會(huì)產(chǎn)生信息的丟失，故先作L倍的增抽樣，再作M倍的抽取，過程如圖3所示。

假設(shè)原抽樣率和抽樣間隔為fs1和T1，抽樣率轉(zhuǎn)換后的抽樣間隔為fs2和T2，抽樣率的比值為fs1／fs2＝T1／T2＝L／M稱為有理比L／M倍抽樣率轉(zhuǎn)換，可用L倍增抽樣級(jí)聯(lián)M倍減抽樣實(shí)現(xiàn)。即先對(duì)離散序列x［n］進(jìn)行L倍增抽樣得到序列z［n］，然后再經(jīng)過M倍減抽樣，實(shí)現(xiàn)x［n］的有理比L／M倍抽樣率轉(zhuǎn)換［10－12］。

圖3中內(nèi)插濾波器的單位激沖響應(yīng)hL（n）和抽取濾波器的帶限低通濾波器hd（n）均在Lfs1工作的級(jí)聯(lián)系統(tǒng)中，因此，可將它們等效為一組合濾波器h（n），其離散時(shí)間傅里葉變換分別為L（Ω）、d（Ω）和（Ω）。

可得到主值期間（－π，π］內(nèi)的頻率響應(yīng)為：

此組合低通濾波器理論上應(yīng)具有理想低通特性，其低通截止頻率等于內(nèi)插濾波器和抽取濾波器低通截止頻率的最小值。

一般抽樣率轉(zhuǎn)換系統(tǒng)的輸入輸出頻率的關(guān)系為：

圖2 信號(hào)增抽樣及頻譜的變換Fig.2 Signal increase sampling and frequency spectrum transform

圖3 插值（a）和抽取（b）的級(jí)聯(lián)實(shí)現(xiàn)Fig.3 Cascade realization of interpolation（a）and extraction（b）

有理比抽樣率轉(zhuǎn)換系統(tǒng)的輸入輸出關(guān)系為：

需注意的是，欲使有理比抽樣率轉(zhuǎn)換時(shí)不產(chǎn)生混疊，L倍內(nèi)插必須先于M倍抽取。否則可能產(chǎn)生混疊。

1.2 儀器譜合成原理

利用數(shù)字信號(hào)處理的方法對(duì)γ射線儀器譜按能量特征進(jìn)行一致化變換，變換后的數(shù)據(jù)相加得到合成譜數(shù)據(jù)。儀器譜數(shù)據(jù)可看作是在能量軸上的離散能量信號(hào)，每一離散能量對(duì)應(yīng)的數(shù)據(jù)是離散能量間距內(nèi)的脈沖計(jì)數(shù)，由此建立數(shù)字信號(hào)與譜數(shù)據(jù)間的關(guān)系模型［10］。

由不同的譜儀測(cè)得的儀器譜差異為：能量起始位置不同；相同能量的全能峰對(duì)應(yīng)道址不一致。采用抽樣率轉(zhuǎn)換，使差異較大的儀器譜按能量特征進(jìn)行平移、尺度變換，使之在相同的道址上有相同的能量，最后將相同能量道的計(jì)數(shù)相加得到合成儀器譜。以NaI（Tl）γ譜儀測(cè)得的天然γ射線儀器譜為例，其合成的步驟如下。

1）根據(jù)所關(guān)注的能量段，選用相近能量的標(biāo)準(zhǔn)源或天然核素的全能峰，對(duì)每條儀器譜進(jìn)行能量刻度。以天然鉀峰和釷峰為例，尋找第m個(gè)儀器譜鉀峰峰位和釷峰峰位。其中鉀峰峰位對(duì)應(yīng)的能量為1 461keV，釷峰峰位對(duì)應(yīng)的能量為2 615keV，進(jìn)行能量刻度，得到儀器譜的能量刻度系數(shù)。

設(shè)m條晶體測(cè)得儀器譜的能量刻度曲線為：

其中：E為能量；C為道址；am、bm為第m＃儀器譜的能量刻度系數(shù)。

令E＝30keV，根據(jù)式（4）可得1組道址Cm。對(duì)各組儀器譜數(shù)據(jù)進(jìn)行移位運(yùn)算，設(shè)第m個(gè)譜數(shù)據(jù)的序列為xm［n］，則對(duì)該序列左移Cm位，得到一新的序列ym［n］，該新序列的第零道對(duì)應(yīng)的能量即為30keV：

2）對(duì)新的序列ym［n］進(jìn)行抽樣率轉(zhuǎn)換。由于新鉀峰峰位是在原來峰位基礎(chǔ)上左移Cm，則可得到1組新的鉀峰峰位：

設(shè)將鉀峰移至目標(biāo)道A，則有比值式：

其中，（L／M）m為第m個(gè)譜線的擴(kuò)展系數(shù)。

對(duì)ym［n］進(jìn)行（L／M）m倍抽樣率轉(zhuǎn)換，先將ym［n］作插值運(yùn)算，再作抽取運(yùn)算得到新的譜數(shù)據(jù)wm［n］，此時(shí)鉀峰峰位轉(zhuǎn)換至A道。將m個(gè)儀器譜按目標(biāo)鉀峰峰位進(jìn)行抽樣率轉(zhuǎn)換，儀器譜的相同道址對(duì)應(yīng)的能量也相同。

3）對(duì)轉(zhuǎn)換后的序列wm［n］進(jìn)行和運(yùn)算，即：

w［n］即為最后合成的儀器譜數(shù)據(jù)序列。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)及結(jié)果

利用NaI（Tl）γ譜儀和HPGeγ譜儀測(cè)得的儀器譜進(jìn)行譜合成實(shí)驗(yàn)。實(shí)驗(yàn)采用4個(gè)10.16cm×10.16cm×40.64cm大體積 NaI（Tl）γ譜儀在實(shí)驗(yàn)室本底下測(cè)得的儀器譜及4個(gè)HPGeγ譜儀在屏蔽室里對(duì)弱標(biāo)準(zhǔn)源進(jìn)行測(cè)量所得的儀器譜作為實(shí)驗(yàn)對(duì)象。將測(cè)得的儀器譜利用抽樣率轉(zhuǎn)換的方法進(jìn)行合成。NaI（Tl）γ儀器譜及HPGeγ儀器譜合成前后效果示于圖4。

圖4 NaI（Tl）（a）和 HPGe（b）探測(cè)器測(cè)得的儀器譜及合成譜Fig.4 Instrument spectra and synthetic spectra measured by NaI（Tl）（a）and HPGe（b）detector

圖4a中計(jì)數(shù)最大的譜線是以1＃晶體測(cè)得的儀器譜為標(biāo)準(zhǔn)的合成譜，其它儀器譜經(jīng)抽樣率轉(zhuǎn)換后，在能量坐標(biāo)上與1＃晶體測(cè)得的儀器譜重合。圖4b中計(jì)數(shù)最大的譜線是以探測(cè)器A測(cè)得的儀器譜為標(biāo)準(zhǔn)的合成譜，其它儀器譜通過抽樣率轉(zhuǎn)換后，在能量坐標(biāo)上與探測(cè)器A測(cè)得的儀器譜重合。轉(zhuǎn)換后的數(shù)據(jù)在相同的道址上有相同的能量，將相同能量的計(jì)數(shù)相加得到合成譜。

通過對(duì)合成前后儀器譜特征能量的峰位、峰面積和半高寬的分析可知，原來差別較大的4條儀器譜經(jīng)一致化后，各條譜線已較好地在能量坐標(biāo)上重合。

表1、2所列數(shù)據(jù)（由Genie2000軟件計(jì)算）分別為NaI（Tl）及HPGeγ射線儀器譜一致化前后特征峰參數(shù)。其中，原譜為γ譜儀直接測(cè)得的儀器譜，轉(zhuǎn)換譜為經(jīng)抽樣率轉(zhuǎn)換后的儀器譜。通過儀器譜和測(cè)試數(shù)據(jù)可看出，對(duì)于陣列NaI（Tl）和HPGeγ譜儀，合成前各儀器譜差異均較大，但經(jīng)抽樣率轉(zhuǎn)換后各儀器譜能較好地在能量坐標(biāo)上重合。

2.2 結(jié)果分析

對(duì)于陣列NaI（Tl）探測(cè)器，鉀、釷峰合成的目標(biāo)峰位與合成的實(shí)際峰位相差小于1道（約3keV）；反映各譜儀能量分辨率優(yōu)良的鉀峰和釷峰的半高寬相對(duì)變化小于1%；合成前各鉀峰的凈峰面積之和與合成后鉀峰凈峰面積相差小于0.7%。

表1 NaI（Tl）譜儀合成前后各儀器譜40K、Th峰參數(shù)Table 1 Peak parameters of 40K and Th in different instrument spectra before and after synthesis by NaI（Tl）spectrometer

表2 HPGe譜儀合成前后各儀器譜241Am、137Cs峰參數(shù)Table 2 Peak parameters of 241Am and 137Cs in different instrument spectra before and after synthesis by HPGe spectrometer

對(duì)于陣列 HPGe探測(cè)器，241Am、137Cs全能峰合成的目標(biāo)峰位與合成的實(shí)際峰位相差小于1道（約為0.1keV）。由于HPGeγ譜儀的計(jì)數(shù)較低，統(tǒng)計(jì)漲落較大，計(jì)算峰的半高寬時(shí)誤差較大，但合成譜的半高寬與單個(gè)HPGeγ譜儀經(jīng)約4倍測(cè)量時(shí)間降低漲落后的相應(yīng)全能峰的半高寬一致。合成前各59.54keV全能峰的凈峰面積之和與合成后59.54keV全能峰凈峰面積相差小于0.1%；合成前各661.66keV全能峰的凈峰面積之和與合成后661.66keV全能峰凈峰面積相差小于4.5%。

3 結(jié)論

采用數(shù)字信號(hào)處理的方法分別對(duì)陣列NaI（Tl）及HPGeγ譜儀進(jìn)行了儀器譜合成方法研究及合成實(shí)驗(yàn)。測(cè)試結(jié)果表明，在陣列NaI（Tl）、HPGeγ譜儀中4個(gè)未合成譜的峰面積之和與合成譜的峰面積相差均小于1%；合成前后NaI（Tl）譜儀半高寬小于1%。HPGe譜儀由于全能峰的計(jì)數(shù)低，造成在計(jì)算半高寬時(shí)誤差較大，但與單個(gè)譜儀測(cè)量約4倍時(shí)間后的半高寬一致。

該儀器譜合成方法解決了硬件方法對(duì)陣列探測(cè)器的儀器譜合成時(shí)，不同探測(cè)器及硬件參數(shù)不一致造成的相同能量的γ射線產(chǎn)生的脈沖幅度不一致的問題。該方法在車載掃描搜索系統(tǒng)及肺部計(jì)數(shù)系統(tǒng)中進(jìn)行了相關(guān)實(shí)驗(yàn)，具有較好的應(yīng)用效果。

［1］王紅艷，劉森林，潘大金.車載式NaI（Tl）大晶體組的多道譜儀［J］.原子能科學(xué)技術(shù)，2004，38（7）：252－255.WANG Hongyan，LIU Senlin，PAN Dajin.A vehicle－based NaI（Tl） multiple channel spectrometer［J］.Atomic Energy Science and Technology，2004，38（7）：252－255（in Chinese）.

［2］陸士立，翟玉貴，馬艷芳，等.CZD－6型車載多道γ能譜測(cè)量系統(tǒng)［J］.鈾礦地質(zhì)，1998，14（2）：112－114.LU Shili，ZHAI Yugui，MA Yanfang，et al.Car－borne multichannel gamma－ray spectrometric system model CZD－6［J］.Uranium Geology，1998，14（2）：112－114（in Chinese）.

［3］Explorainum.GR－820user manual［M］.［S.l.］：［s.n.］，1996.

［4］王經(jīng)瑾.核電子學(xué)［M］.北京：原子能出版社，1983.

［5］復(fù)旦大學(xué)，清華大學(xué)，北京大學(xué)三校合編.原子核物理實(shí)驗(yàn)方法［M］.北京：原子能出版社，1997.

［6］汲長松.核輻射探測(cè)器及其實(shí)驗(yàn)技術(shù)手冊(cè)［M］.北京：原子能出版社，1990.

［7］HASHEM M H.Nonlinear response function of a 3×3in.NaI scintillation detector［J］.Asian J Exp Sci，2007（5）：127－130.

［8］γ－ray spectroscopy using a sodium iodide scintillation detector［G］.UK：University of Sheffield，2008：259－263.

［9］殷經(jīng)鵬，賴萬昌，王廣西，等.陣列 NaI（Tl）探測(cè)器γ能譜合成方法探討［J］.計(jì)測(cè)技術(shù)，2009，29（3）：4－6.YIN Jingpeng，LAI Wanchang，WANG Guangxi，et al.The discussion on synthesis method of array NaI（Tl）detectorγ－ray spectrum［J］.Metrology ＆ Measurement Technology，2009，29（3）：4－6（in Chinese）.

［10］殷經(jīng)鵬.陣列NaI（Tl）探測(cè)器γ能譜合成方法研究［D］.成都：成都理工大學(xué)，2009.

［11］喻名德，楊春才.核試驗(yàn)場(chǎng)及其治理［M］.北京：國防工業(yè)出版社，2007.

［12］徐守時(shí).信號(hào)與系統(tǒng)［M］.北京：清華大學(xué)出版社，2008.

［13］薛年喜.Matlab在數(shù)字信號(hào)處理中的應(yīng)用［M］.北京：清華大學(xué)出版社，2003.