李自維，白立新，張一云

(四川大學物理學院，成都 610065)

1 引 言

使用Monte Carlo方法模擬高純鍺探測器對γ光子的探測效率相比傳統(tǒng)的物理實驗方法能夠大大地節(jié)約時間和費用[1]，在高純鍺探測器參數(shù)標定和修正中有著廣泛的應用. 探測器生產(chǎn)廠家給出探測器參數(shù)為標稱值，而冷指尺寸、真空層厚度等亦未給出，且探測器老化后探測器參數(shù)也會變化. 因此根據(jù)廠家給出的探測器參數(shù)，進行探測效率模擬，得到的結果往往與實驗值存在較大的偏差[2]. 為了得到準確模擬結果需要對計算模型的輸入?yún)?shù)進行調(diào)整. 探測效率由探測器的多個參數(shù)共同影響，目前文獻中對探測器參數(shù)的調(diào)整方法還處于人工試探、窮舉法等，沒有一種理論性強、簡便、通用的方法. 本文以計算模型的統(tǒng)計性不確定度評價方法為基礎，通過對探測效率的蒙特卡洛計算模型進行敏感性分析和不確定度分析，實現(xiàn)對高純鍺探測器重要參數(shù)的修正，在方法上具有通用性和一致性.

2 原理及實驗

2.1 高純鍺探測器原理

高純鍺探測器是一種具有高分辨率的半導體核輻射探測器，其用于探測射線的靈敏區(qū)域是PN結(又稱耗盡層). 當γ射線進入結區(qū)時，因γ光子與結區(qū)物質(zhì)發(fā)生光電效應，康普頓散射或電子對效應從而損失能量，在結區(qū)產(chǎn)生電子-空穴對，在外加反偏電場的作用下，在輸出回路中形成能代表入射γ射線能量的輸出電信號. 探測器結構如圖1所示，高純鍺探測器需要工作在85～100 K的低溫環(huán)境，來保證鍺晶體的禁帶寬度.

圖1 探測器結構示意圖Fig.1 Sketch map of HPGe detector

2.2 與探測效率有關的探測器參數(shù)

探測效率是探測器的重要指標之一，由多參數(shù)共同影響，其準確性直接影響到放射源活度等物理測量結果的準確度[3].

對于特定的放射源和確定的高純鍺探測器，探測效率主要與以下幾個探測器參數(shù)相關：(1)探測器晶體尺寸，包括晶體半徑和晶體高度，探測器靈敏區(qū)域的大小主要由晶體的尺寸決定；(2)死層厚度，死層的厚度直接影響探測器晶體的有效體積；(3)冷指尺寸，包括冷指長和冷指半徑，入射粒子進入到冷指井中不產(chǎn)生空穴-電子對，冷指的大小也直接影響到探測器晶體的有效體積；(4)空隙高度，入射窗到探測器托架之間的距離稱為空隙高度，空隙高度的偏差影響放射源到探測器的距離.

2.3 計算模型的不確定度評價方法

對計算模型的不確定度評價，一般指對計算模型的輸出結果與實際值之間的差異進行評價. 包括分析輸入?yún)?shù)對輸出結果的不確定度貢獻和輸出參數(shù)對輸入?yún)?shù)的敏感性.

2.3.1 不確定度評價方法分類和選擇 不確定度評價分析方法根據(jù)輸入不確定度到輸出不確定度的傳遞方向，分為基于輸入?yún)?shù)的不確定度傳遞和基于輸出的不確定度傳遞的不確定度評價分析.

根據(jù)對計算模型的輸出結果的不確定度計算方法不同，不確定度評價方法分為統(tǒng)計性方法和確定性方法. 目前，基于輸入?yún)?shù)不確定度傳遞的統(tǒng)計性方法的研究和應用相對更為成熟和廣泛.

統(tǒng)計性不確定度評價方法首先確定輸入?yún)?shù)的分布類型和分布參數(shù)，然后對各輸入?yún)?shù)抽樣生成N個輸入樣本，用此樣本集進行模擬計算得到N個輸出集，最后分析輸出與輸入的關系. 統(tǒng)計性不確定度評價方法是不確定度評價的一種通用方法，示意圖如圖2所示.

圖2 統(tǒng)計性不確定度評價方法Fig.2 Statistical uncertainty evaluation method

探測器效率模擬的模型是一個多輸入?yún)?shù)的復雜模型，輸出的結果與各輸入?yún)?shù)并非簡單的線性關系，因此需采用基于輸入?yún)?shù)不確定度傳遞的統(tǒng)計性不確定度評價方法，來分析評價探測器真實效率與模擬效率的關系.

2.3.2 敏感性分析方法 敏感性分析是指分析計算模型的每個輸入變量和參數(shù)在模型輸出變量中的相對重要性[4]. 敏感性分析可用于輸入變量和輸入?yún)?shù)篩選，輸出變量的不確定度分析，減少輸出變量不確定度的方法等，其廣泛應用于機械、交道運輸網(wǎng)絡建模、熱工水力、電力系統(tǒng)分析和控制以及市場營銷、運籌學等領域[5].

敏感性分析分為局部敏感性分析和全局敏感性分析[6]. 局部敏感性分析法關注單個因子對模型輸出結果的影響，優(yōu)點在于簡單快捷，可操作性強，但是由于未考慮整個參數(shù)空間內(nèi)的響應，當模型是非線性的或者影響輸入變量的不確定度處于不同數(shù)量級時，局部技術可能不能夠提供有效的分析結果.

全局敏感性分析法分析各因素的概率密度函數(shù)的分布和形狀的影響，在計算分析時，所有因素都可同時變化[7-8]. 所以，本文采用簡相關系數(shù)、全局敏感性分析方法Sobol′一階和總敏感性分析方法對模型輸入?yún)?shù)進行敏感性分析.

2.4 不確定評價方法修正參數(shù)的原理

在高純鍺探測器對放射源的探測效率模擬中，不確定度的傳遞示意圖如圖3所示.探測器的參數(shù)是模擬模型的輸入?yún)?shù)，探測器參數(shù)的不確定度直接反映到模擬的結果中，即模擬效率的不確定度來源于Monte Carlo N-Particle Transport Code，Version 5(MCNP5)程序的模型輸入?yún)?shù)(探測器參數(shù))的不確定度.

圖3 探測效率模擬中的不確定度傳遞Fig.3 Uncertainty transfer in detection efficiency simulation process

用模擬程序模擬高純鍺探測器的探測效率，模擬結果與實驗效率存在一定偏差，其主要原因是建模的輸入?yún)?shù)值與探測器真實參數(shù)值的存在偏差. 而輸入?yún)?shù)值與真實參數(shù)值的偏差程度決定了模擬探測效率與實驗效率的偏差程度. 所以，我們認為當模擬探測效率與實驗效率偏差最小時，此模型的輸入?yún)?shù)與真實參數(shù)偏差最小.

本文提出一種不確定度分析法，即對有代表性的多個空間點(對輸入?yún)?shù)敏感性高的點)的點源探測效率，采用簡單蒙卡抽樣方法對重要輸入?yún)?shù)的每個參數(shù)同時進行N次均勻抽樣，使抽樣值覆蓋探測器真實參數(shù)值，由抽樣數(shù)據(jù)生成模擬程序的N個輸入卡. 再由探測器效率模擬程序進行N次模擬計算，求出模擬計算結果與真實探測效率的偏差，找出偏差最小值所對應的探測器效率模擬程序輸入卡，即可知此最小偏差對應的模擬模型的輸入?yún)?shù)，完成探測器參數(shù)修正.

2.5 方法驗證實驗

2.5.1 實驗假設和實驗目的 本實驗合理地假設一組高純鍺探測器參數(shù)，作為模擬實驗的高純鍺探測器的真實參數(shù). 實驗中由這些假設的真實探測器參數(shù)作為模型參數(shù)模擬計算出的探測效率即認為是真實效率或稱為實驗效率. 假設另一組高純鍺探測器參數(shù)，作為模擬實驗的高純鍺探測器的標稱值. 由標稱值作為模型參數(shù)模擬計算出的探測效率稱為標稱值探測效率.

實驗選用137Cs，241Am，60Co做標準γ源，選取9個空間點做為放射源所在位置，這9個點與探測器位置關系如同4所示，第1～9個點位置分別記為位置1，位置2 ,直到位置9.

圖4 實驗選取的空間點Fig.4 Space points in the experiment

實驗中，用MCNP5程序模擬高純鍺探測器探測效率，應用前文提出的不確定度評價方法來修正高純鍺探測器參數(shù)，最終比較修正的探測器參數(shù)與實驗假設的探測器真實參數(shù)的誤差，并驗證此不確定度評價方法的正確性.

2.5.2 參數(shù)修正 根據(jù)前文所提出的不確定度評價方法，對探測器參數(shù)進行修正. 第一步，分析高純鍺探測器參數(shù)，選取目標參數(shù)，分析探測效率對各目標參數(shù)的敏感性，確定重要參數(shù).

第二步，對各重要輸入?yún)?shù)同時進行N次簡單蒙卡抽樣，得到N組模型輸入?yún)?shù)，根據(jù)敏感性分析結果，選取n個標定點，對不同放射源在這n個標定點，用MCNP5程序進行同時的探測器效率模擬，得到模擬計算結果.

第三步，計算探測器對各空間點對應放射源的真實探測效率與模擬效率的偏差的平方之和，表示為ΔEff.

其中，εsi表示模擬效率，εti表示真實效率.找到最小ΔEff，這個ΔEff對應的那組模型參數(shù)值即可認為是最佳參數(shù)值，也就是探測器參數(shù)的修正值.

本實驗采用針對多個有代表性的空間點，對不同的放射源進行同時的探測器效率模擬，是因為探測器探測效率是由多參數(shù)共同影響的，在進行探測器對確定位置的單一放射源進行探測效率模擬時，容易出現(xiàn)多個模型輸入?yún)?shù)均與真實值偏移較大，而模擬效率卻接近真實效率的情況. 而采用多空間點多放射源進行同時模擬，能很好地避免這種情況.

2.5.3 參數(shù)驗證 根據(jù)修正的高純鍺探測器參數(shù)，令137Cs，241Am，60Co這3個γ放射位于除標定點外的其他m個點，用MCNP5程序模擬探測器對其探測效率，并與用實驗假設的探測器真實參數(shù)模擬計算的探測效率進行對比，以此檢驗探測器參數(shù)的修正結果.

2.6 探測效率的MCNP程序模擬[9]

本文采用Monte Carlo光子和電子耦合輸運程序MCNP5對高純鍺探測器探測效率進行模擬.MCNP程序是一款應用廣泛[10]，由Los Alamos國家實驗室研發(fā)的基于蒙特卡羅方法的大型的多功能Monte Carlo粒子輸運程序. 該程序目前可用于計算與中子、光子和電子或耦合中子、 光子、 電子相關的物理輸運問題[11].

在光子和電子的耦合輸運過程中，模擬計算的光子能量為1 keV～100 MeV，適合于γ射線的模擬計算[12]. 在探測效率模擬中，用脈沖幅度卡F8計算γ射線在HPGe晶體中的脈沖高度能譜分布，聯(lián)合使用E8計數(shù)卡，即可求得探測器的模擬效率.

2.7 敏感性分析

2.7.1 確定輸入?yún)?shù)及參數(shù)抽樣 分析影響探測效率的模型參數(shù)，選取探測器晶體半徑，空隙高度，死層厚度，冷指半徑，冷指長度，晶體高度共6個輸入?yún)?shù)作為進行敏感性分析的目標參數(shù).

各輸入?yún)?shù)的相對偏差并無嚴格的統(tǒng)計數(shù)據(jù)指導. 一般由標稱值模擬的探測效率與真實效率相對誤差在20%以上，因此可由MCNP5程序模擬計算出探測器參數(shù)的最大偏移量. 結合實驗室經(jīng)驗數(shù)據(jù)，假設各輸入?yún)?shù)的分布均服從正態(tài)分布，得到各參數(shù)的標準差，如表1所示.

采用拉丁超立方抽樣方法對輸入?yún)?shù)進行N次同時抽樣產(chǎn)生N個輸入樣本，并保存抽樣結果.

表1 輸入?yún)?shù)分布Tab.1 The distribution of input parameters

2.7.2 模型計算及敏感性計算 由抽樣數(shù)據(jù)編寫N個MCNP5程序的輸入卡，由自編程序調(diào)用MCNP5程序進行模擬計算，保存目標參數(shù). 這里目標參數(shù)只有一個，即模型計算的結果模擬探測效率.

根據(jù)抽樣時保存的輸入?yún)?shù)抽樣結果和模擬計算得到的探測器效率，先后計算各輸入?yún)?shù)的簡相關系數(shù)、Sobol′一階敏感系數(shù)和Sobol′總敏感系數(shù).

3 敏感性分析及參數(shù)修正

3.1 敏感性分析結果

選取位置2，位置4，位置5，位置8，位置9放置不同放射源做敏感性分析. 首先，對6個輸入?yún)?shù)分別抽樣500、1 000、2 000、3 000次，令137Cs放射源處于位置4，用MCNP5程序進行探測效率模擬，計算出簡相關系數(shù)，Sobol′一階敏感系數(shù)和Sobol′總敏感系數(shù). 當抽樣數(shù)是2 000次時，各個參數(shù)的敏感性計算結果與1 000次抽樣的計算結果偏差在5%之內(nèi)，3 000次抽樣計算結果與2 000次抽樣的計算結果偏差在2%之內(nèi)，故此認為當抽樣數(shù)為3 000時，敏感性分析結果收斂. 下文敏感性分析結果均為抽樣3 000次時的數(shù)據(jù).

這里列出137Cs源位于位置3、241Am源位于位置5，60Co源位于位置6的敏感性分析結果，如表2所示.

表2 敏感性分析結果Tab.2 The results of sensitivity analysis

由結果可見相關系數(shù)和Sobol′敏感系數(shù)具有一致性，這證明了在MCNP5程序模擬探測器效率中應用Sobol′敏感性分析的正確性. 結果表明，模擬探測效率對探測器參數(shù)的敏感度與放射源位置和入射γ光子能量均有關系：當入射粒子能量低時，探測效率對死層厚度較為敏感，而入射粒子能量高，且放射源位置不在探測器晶體軸向時，探測效率對晶體半徑、高度敏感度比較大.

綜合分析不同放射源在其他位置的敏感性研究結果，高純鍺探測器探測效率對死層厚度、晶體半徑、晶體高度這3個探測器參數(shù)的敏感性始終大于對冷指半徑、冷指長和空隙高度這3個參數(shù)的敏感性，且當放射源偏離軸心時，探測效率對晶體高度更敏感. 因此，晶體的半徑和死層厚度應該作為探測器重要參數(shù). 高純鍺探測器探測效率對冷指長度和冷指半徑敏感性很小，所以冷指長度和冷指半徑不作為重要參數(shù).

3.2 模型參數(shù)修正

3.2.1 尋找最佳模型參數(shù) 根據(jù)敏感性分析的結果，選取晶體高度、晶體半徑，死層厚度這三個重要參數(shù)進行修正，用于參數(shù)修正的標定點應該選偏離軸心的點.

表3 重要模型參數(shù)值Tab.3 Important model parameter values

選取位置3，位置5，位置6作為標定點，修正參數(shù)所用γ放射源仍為241Am，60Co，137Cs.241Am源在位置3，137Cs源在位置5，60Co源在位置6. 對這三個位置的點源，用MCNP5程序模擬高純鍺探測器對他們的探測效率. 實驗所假設的探測器重要參數(shù)真實值與標稱值如表3所示.

經(jīng)MCNP5程序模擬計算的探測器對各放射源的真實探測效率和使用標稱值計算的模擬效率結果如表4所示.

表4 不同γ源在不同位置的探測器探測效率Tab.4 The detection efficiencies for different γ sources located at different spatial points

采用簡單蒙卡抽樣方法對晶體高度、晶體半徑死層厚度共3個輸入?yún)?shù)，同時進行均勻抽樣3 000次，生成3 000張MCNP5程序輸入卡并進行模擬計算.

表5 探測器參數(shù)修正結果Tab.5 The results of corrected important detector parameter values

求出3個γ源在3個不同位置的模擬探測效率與真實探測效率的偏差的平方和ΔEff，并找出最小值，得出對應輸入?yún)?shù)的值. 這些參數(shù)值即視為探測器參數(shù)的修正值，修正值如表5所示. 結果表明，修正的探測器參數(shù)非常接近真實值，各個參數(shù)與真實值的相對誤差均在0.5%內(nèi). 采用修正后的探測器參數(shù)建模計算高純鍺探測器對不同γ點源在3個標定位置的探測效率，結果如表6所示，模擬探測效率與真實的探測效率相對誤差均在0.5%之內(nèi).

表6 參數(shù)修正后的探測器探測效率Tab.6 The detection efficiencies after parameters corrected

重復上述參數(shù)修正步驟12次，共有11次修正得到重要參數(shù)的修正值與真實值相對誤差均未超過2%. 有1次探測器參數(shù)修正值與真實值相對誤差在10%左右，而模擬效率與真實效率相對誤差在1%內(nèi). 這是前文所提到的修正的偶然性. 對于這種修正結果直接舍棄.

3.2.2 修正參數(shù)驗證 選取位置2，4，8作為驗證點，采用修正后的探測器參數(shù)，使用MCNP5程序建模分別計算137Cs，60Co，241Am放射源處在驗證點的探測器探測效率，以此驗證探測器重要參數(shù)修正是否正確. 對各源，高純鍺探測器的探測效率真實值與模擬的探測效率結果如表7所示.

結果表明，使用修正后的探測器參數(shù)對不同點源在驗證點進行探測效率模擬，結果與真實效率相對誤差均在0.5%之內(nèi)，這遠好于傳統(tǒng)修正方法修正的結果(傳統(tǒng)修正方法的修正結果一般相對誤差在2%～5%之內(nèi)).這很好地說明了應用此方法修正探測器重要參數(shù)的正確性和準確性.

表7 探測器對驗證源的探測效率Tab.7 The detection efficiencies for verification sources

4 結 論

使用MCNP5程序能快速地模擬高純鍺探測器的探測效率，節(jié)約實驗成本，縮短實驗時間. 應用本文提出的不確定度分析方法來輔助標定或修正探測器參數(shù)，修正后的探測器參數(shù)與真實探測器參數(shù)相對誤差在2%以內(nèi).此方法具有良好的重復性，參數(shù)修正效果好于傳統(tǒng)的高純鍺探測器參數(shù)修正方法. 此方法同樣也適用于其他類型的核探測的探測器參數(shù)標定與修正.

不確定評價方法在具有多輸入?yún)?shù)且復雜的模型計算的熱工水力、機械結構、電力輸運等工程領域已經(jīng)有廣泛的應用和研究[13].在實驗室領域，對于較為簡單的計算模型，應用不確定度分析的思想，應用統(tǒng)計性不確定度評價方法來分析輸入?yún)?shù)或模型參數(shù)與輸出目標參數(shù)的關系也是實用且便捷的.