顏筱宇 韓紀(jì)鋒 黃宇 宋瑞強(qiáng) 徐子虛 劉星泉 曲國峰 張藝蓉 冷強(qiáng)鐘 任飛旭 張鑫 陳蕾 劉吉珍

本文提出了一種基于碳化硼慢化體和涂硼微結(jié)構(gòu)中子探測器的通用中子能譜儀設(shè)計方法，獲得了計算最優(yōu)慢化體厚度和探測器響應(yīng)函數(shù)的通用公式，并用蒙特卡羅方法進(jìn)行了驗證和修正，可實現(xiàn)對中子能譜儀的快速設(shè)計. 該類能譜儀可實現(xiàn)各中子響應(yīng)函數(shù)之間的解耦，使每個探測器對各個分立能量區(qū)間的中子最為敏感，具有較強(qiáng)適用性和靈活性. 基于該方法，我們設(shè)計了一個用于硼中子俘獲治療（BNCT）超熱中子能譜測量的能譜儀，通過Gravel算法實現(xiàn)了中子能譜解析，并提出了一種具有良好普適性的基于響應(yīng)函數(shù)的預(yù)置譜設(shè)置方法. 結(jié)果表明，該能譜儀對單能中子的峰位解析精度約為1%，對BNCT連續(xù)譜的解析均方差約為0.76%，具有較大技術(shù)優(yōu)勢和可行性.

中子譜儀；蒙特卡羅方法；少道解譜算法；硼中子俘獲治療

TL816+.3A2023.024003

收稿日期： 2022-10-17

基金項目： 四川省科技計劃（2022NSFSC0831，2020YJ0313）; 中央高?；究蒲袠I(yè)務(wù)費

作者簡介： 顏筱宇（1996-）， 男， 山東曲阜人， 碩士研究生， 研究領(lǐng)域為中子探測.

通訊作者： 韓紀(jì)鋒. E-mail： hanjf@scu.edu.cn

A general neutron spectrometer design method based on boron moderator and boron-coated micro-structure neutron detectors

YAN Xiao-Yu1， HAN Ji-Feng1， HUANG Yu1， SONG Rui-Qiang1， XU Zi-Xu2， LIU Xing-Quan1， QU Guo-Feng1， 3， ZHANG Yi-Rong1， LENG Qiang-Zhong1， REN Fei-Xu1， ZHANG Xin1， CHEN Lei1， LIU Ji-Zhen4

（1. Key Laboratory of Radiation Physics and Technology of the Ministry of Education， Institute of Nuclear Science and Technology， Sichuan University， Chengdu 610064， China;

2. Graduate School of Engineering， Osaka University， Osaka 565-0871， Japan;

3. Helmholtz Institute， Johannes Gutenberg-Universitt Mainz， Mainz 55099， Germany;

4. Nuclear Power Institute of China， Chengdu 610041， China）

In this paper， a general neutron spectrometer design method based on boron carbide moderators and boron-coated micro-structure neutron detectors is proposed. The formula for calculating the optimal moderator thickness and detector response function according to the neutron energy is obtained， which is verified and modified by the Monte Carlo simulation. The formula is able to realize the rapid design of the neutron spectrometer effectively. This kind of energy spectrometer can realize the decoupling of neutron response functions， so that each detector is most sensitive to neutrons in each discrete energy range， and has strong applicability and flexibility. An epithermal neutron spectrometer for boron neutron capture therapy （BNCT） is designed based on the proposed method， the response function of the spectrometer is calculated and the neutron spectrum unfolding is achieved by the Gravel algorithm. The results show that the resolution accuracy of the energy spectrometer to the peak position of monoenergetic neutrons is about 1%， and the analytical mean square error of the BNCT continuous spectrum is about 0.76%， which has great technical advantage and feasibility.

Neutron spectrometer; Monte Carlo method; Less channel spectrum unfolding algorithm; Boron neutron capture therapy

1 引 言

中子能譜的準(zhǔn)確測量，在核醫(yī)學(xué)、核物理基礎(chǔ)研究、核工程、輻射防護(hù)和空間探測等領(lǐng)域有重要意義［1， 2］，也是反應(yīng)堆、加速器等中子輻射場的關(guān)鍵參數(shù)［3］. 因此針對中子能譜儀的研發(fā)與中子能譜的測量等工作一直是該領(lǐng)域的研究重點和難點. 基于核反應(yīng)法的Bonner多球譜儀具有靈敏度較高以及γ排除能力較好等優(yōu)點，在中子能譜測量領(lǐng)域應(yīng)用較為廣泛［4-8］. 而中子響應(yīng)函數(shù)是實現(xiàn)能譜測量的關(guān)鍵. 現(xiàn)有多球譜儀內(nèi)各個探測器均對很寬能量區(qū)間內(nèi)的中子有響應(yīng)，響應(yīng)曲線的能量區(qū)間重合度較高，僅在響應(yīng)靈敏度上有所差異. 即， 各響應(yīng)函數(shù)存在耦合，難以解耦，導(dǎo)致能譜解析困難、能量分辨率不高等問題. 南華大學(xué)的鄒益晟等人［9］基于多球譜儀設(shè)計出一款解譜性能較好的基于閃爍體探測器的單球譜儀. 此儀器具有結(jié)構(gòu)緊湊和操作簡便的優(yōu)點，但由于其響應(yīng)函數(shù)能量區(qū)間重合度較高等問題導(dǎo)致能譜解析困難，其性能的進(jìn)一步提升也受限于響應(yīng)函數(shù)的固有特征. 目前，針對多球/單球中子譜儀的改進(jìn)主要集中在提高熱中子的探測效率和提高測量中子能量探測上限等方面［10-12］. 能量分辨率的進(jìn)一步提高則存在較大困難，而且解譜的準(zhǔn)確度受預(yù)制譜的影響較大. 尤其在硼中子俘獲治療（BNCT）領(lǐng)域——主要使用超熱中子和含硼靶向藥物用于腫瘤的定向治療，超熱中子能譜的準(zhǔn)確探測對能譜儀提出了新的要求. 日本國家防衛(wèi)大學(xué)的Takada等人［13］提出使用碳化硼（B4C）對中子進(jìn)行慢化與吸收、利用LiF探測器進(jìn)行中子計數(shù)測量的方法，能夠較好地使各個探測器響應(yīng)曲線的峰分開，使得探測器的分辨率有所提高，并顯著提升了熱中子和超熱中子區(qū)域的探測效率，驗證了該方法具有較大技術(shù)優(yōu)勢.

在中子能譜解析算法方面，國內(nèi)外相關(guān)團(tuán)隊發(fā)展了類似SAND-Ⅱ算法、Gravel算法、遺傳算法（GA）和極大似然算法（MLEM），以及基于神經(jīng)網(wǎng)絡(luò)的多種解譜算法，用于對中子能譜進(jìn)行解析［14-17］，已較為成熟. 其中， Gravel算法由SAND-Ⅱ算法改進(jìn)而來，是1996年由德國 PTB 實驗室提出的一種交互式迭代算法. 它具有低噪聲和解為非負(fù)等優(yōu)點，自誕生以來在少道算法方面應(yīng)用廣泛，不僅可以對于多球譜儀的響應(yīng)進(jìn)行中子能譜解析［8，12，17］，而且對反沖質(zhì)子法和活化法測量中子能譜［18，19］的解譜精度也很高. 該算法同樣適用于γ能譜解析［20］. 中國輻射防護(hù)研究院的李建偉等人［21］利用MCNP模擬獲得了幾種輻射防護(hù)領(lǐng)域的常見中子能譜，而后使用迭代算法進(jìn)行解譜研究. 其結(jié)果表明， 基于迭代法的解譜算法精度較高，但同時對預(yù)置譜精度的要求也較高. 而各種場景下的中子能譜可能非常復(fù)雜且具有顯著差異，通常難以獲得高精度的預(yù)置譜.這制約了解譜精度的進(jìn)一步提升.

本文針對硼、鋰元素的中子反應(yīng)截面規(guī)律，提出了一種使用硼、鋰作為慢化體和靈敏探測器的通用中子能譜儀設(shè)計方法，并使用碳化硼作為慢化體、涂硼微結(jié)構(gòu)中子探測器作為靈敏體進(jìn)行了設(shè)計驗證. 該方法可實現(xiàn)各個靈敏中子探測器響應(yīng)函數(shù)之間的解耦，在通用中子能譜設(shè)計上具有較大技術(shù)優(yōu)勢. 基于該方法，我們設(shè)計了一個用于BNCT超熱中子能譜探測的中子能譜儀，用碳化硼作為慢化體，包括15個靈敏體探測器. 利用MCNP模擬我們獲得了該能譜儀在各種中子輸入下的探測器響應(yīng)，驗證了設(shè)計方法的準(zhǔn)確性. 最后，本文提出了基于響應(yīng)函數(shù)的預(yù)置譜設(shè)置方法，并利用Gravel算法實現(xiàn)了能譜解析，可反解得到相應(yīng)的中子能譜. 我們用標(biāo)準(zhǔn)BNCT能譜進(jìn)行了測試和驗證［22］ ，結(jié)果表明該設(shè)計方法具有較好的適用性和技術(shù)優(yōu)勢.

10B（n，α）7Li的熱中子反應(yīng)截面為3837b，反應(yīng)截面隨能量增長近似符合“1/v規(guī)律”，同10B類似，6Li也具有類似的中子反應(yīng)截面趨勢規(guī)律. 其反應(yīng)截面可以用式（2）近似表示，其中a0為系數(shù)，E為中子能量（單位eV）. 如圖1所示，截面數(shù)據(jù)來自于IAEA ENDF/B-VI. 式（2）在低能區(qū)同實際截面數(shù)據(jù)符合較好，在高能區(qū)則存在一定差異. 對于10B和6Li，系數(shù)a0分別為610.3122和149.5162.本文研究結(jié)果表明，式（2）對截面的近似對低能區(qū)（E＜75 eV）的影響可忽略，對中等能區(qū)（75 eV

σ≈a0 E（2）

當(dāng)采用B/Li元素同時作為慢化體和探測器靈敏體材料時，對于垂直入射的中子束，在慢化體深度D處的通量如式（3）所示［23］.

I=I0e-σNvD（3）

其中， I0為入射中子注量率; σ為反應(yīng)截面; D為慢化體厚度; Nv為慢化體單位體積內(nèi)的原子核數(shù). 若在慢化體深度D處放置一個靈敏體厚度為Dd的中子探測器，Nvd為探測器靈敏體中單位體積的原子核數(shù)，則該探測器可收集到的中子計數(shù)N如式（4）所示［23］.

N=σNv dDdI=I0σNv dDde-σNvD≈a EI0e-bD E（4）

在慢化體材料以及探測器靈敏體涂層材料和厚度（Dd 取0.0137 cm）確定后，式（4）中系數(shù)a、b均為常數(shù). 若探測器靈敏體和慢化體都為碳化硼材料，則a、b分別為0.182和13.292. 需注意式（4）在計算探測器靈敏體厚度較薄時計算準(zhǔn)確，若采用的探測器靈敏體厚度較厚，則需要考慮中子能量和通量在探測器靈敏體中的連續(xù)降低（通量計算方法如式（3）），并對式（4）進(jìn)行修正. 此外，當(dāng)存在多個探測器時，需考慮前面探測器由于吸收或散射中子而對后面探測器產(chǎn)生影響，以及慢化體和探測器的B含量不同的修正等. 本文采用的探測器靈敏體厚度為137 μm，首先用式（4）對探測器響應(yīng)進(jìn)行簡化計算，然后通過蒙卡模擬方法精確計算探測器的響應(yīng)，并根據(jù)蒙卡模擬結(jié)果對式（4）進(jìn)行修正，來補(bǔ)償這些因素的影響，進(jìn)而給出更為準(zhǔn)確的修正后的公式. 通過簡化公式（4）可給出各深度D處靈敏探測器的計數(shù)率N/I0與入射中子能量E之間的關(guān)系，即各靈敏探測器的中子響應(yīng)函數(shù). 為達(dá)到預(yù)測慢化體厚度的目的，需要建立慢化體厚度D與探測器響應(yīng)曲線峰值能量EPeak的關(guān)系. 通過設(shè)置式（4）對能量E的偏導(dǎo)為零，可獲得計數(shù)率N最大時EPeak和D的關(guān)系，結(jié)果如式（5）.

EPeak=b2×D2（5）

根據(jù)式（5）可由中子能量E計算最優(yōu)慢化體厚度D，并通過式（4）獲得在各個慢化體深度D處探測器的中子響應(yīng)函數(shù). 當(dāng)設(shè)置探測中子能量為4.42E-1、4.42E1、4.42E3和4.42E5 eV時，對應(yīng)的最優(yōu)慢化體厚度D分別為0.05、0.5、5和50 cm. 對應(yīng)探測器的計數(shù)率同入射中子能量的關(guān)系曲線如圖2所示. 圖中已對4條曲線根據(jù)峰值進(jìn)行了歸一化處理. 從圖2可見，每條曲線都具有單峰結(jié)構(gòu)，不同深度處的靈敏探測器對不同能量的中子最為敏感，且對某一能量的中子具有最高響應(yīng)，即不同深度處的靈敏探測器對不同能量的中子最敏感.

隨后我們進(jìn)一步研究了各響應(yīng)曲線的峰位能量（EPeak）、慢化體厚度、峰值計數(shù)率、半高寬（FWHM）和相對寬度（FWHM/EPeak）隨慢化體厚度的變化規(guī)律，結(jié)果如表1所示. 結(jié)果表明，隨著慢化體厚度的增加，響應(yīng)曲線的峰值計數(shù)率逐漸降低，峰位能量和半高寬都逐漸增加，相對寬度保持不變，約為18.45.

3 用于BNCT中子能譜測量的能譜儀

3.1 能譜儀設(shè)計

基于上述設(shè)計方法，針對BNCT超熱中子能譜的測量需求，使用天然豐度的碳化硼作為慢化體，使用涂硼硅基微結(jié)構(gòu)中子探測器設(shè)計了一個超熱中子能譜儀，共包括15個靈敏探測器. 由于BNCT主要關(guān)注在0.1 eV～10 keV的超熱能區(qū)，因此將15個探測器可測量的中子峰位能量設(shè)計為近似均勻分布在0.1 eV～10 keV范圍內(nèi)，并通過公式（5）獲得了每個探測器對應(yīng)的慢化體厚度，得到能譜儀的結(jié)構(gòu)示意圖如圖3所示. 圖中每個探測器對應(yīng)的可測量中子峰值能量依次為：0.40、0.80、1.60、3.30、6.20、13.00、24.00、46.00、90.00、175.00、340.00、1300.00、2500.00、5000.00和10 000.00 eV. 在計算過程中，已考慮了前端探測器和慢化體對后端探測器響應(yīng)的影響.

15個探測器的響應(yīng)函數(shù)波形可通過式（4）計算，結(jié)果如圖4a所示. 圖中將所有區(qū)間的峰值進(jìn)行了歸一化以便于對比，每條響應(yīng)曲線的峰值計數(shù)率和半高全寬（FWHM）值如圖4b所示. 圖中可見，15條響應(yīng)曲線具有類似的形狀，其峰值隨慢化體厚度增大而逐漸增加，均勻分布在0.1 eV～10 keV范圍內(nèi)，達(dá)到了設(shè)計預(yù)期要求. 由于碳化硼具有較大的中子吸收截面，因此峰值計數(shù)率隨慢化體厚度的增加而快速衰減. 同截面曲線類似，峰值計數(shù)率隨能量按1/sqrt（E）規(guī)律變化，第1個探測器（峰位0.40 eV）的峰值計數(shù)率為0.121，而第15個探測器（峰位10 keV）的峰值計數(shù)率為6.91E-4，降低了約175倍. 響應(yīng)曲線的FWHM隨探測峰位的增大而增大，但相對寬度固定不變. 相對寬度為R=FWHM/E，其中FWHM為中心能量為E的探測器所對應(yīng)響應(yīng)函數(shù)的半高寬. 此處的相對寬度為中子譜儀對單能中子響應(yīng)函數(shù)的半高寬除以峰值能量，是15個計數(shù)型探測器綜合獲得的結(jié)果，同模擬中的能譜展寬等因素?zé)o關(guān). 同時相對寬度同中子譜儀的中子探測能量分辨率不同，雖然響應(yīng)函數(shù)的相對寬度較大，但通過15個計數(shù)型探測器聯(lián)合解譜，可顯著提高最終的中子能量分辨率，具體見第4節(jié). 3.2 蒙卡模擬驗證

針對設(shè)計的能譜儀，通過MCNP模擬對響應(yīng)函數(shù)進(jìn)行了驗證，檢查設(shè)計方法的可行性. 此外，在能譜儀外圍用厚度1 cm的聚乙烯進(jìn)行包裹，以減弱周圍環(huán)境中散射中子的干擾，結(jié)構(gòu)如圖3所示. 此部分包裹結(jié)構(gòu)以及探測器材料與空氣的影響難以通過公式直接計算，需要通過蒙卡模擬方法進(jìn)行修正. 公式（4）中對于探測器計數(shù)的簡化計算也可以通過蒙卡模擬的結(jié)果進(jìn)行修正.

MNCP5采用美國國家核數(shù)據(jù)中心（NNDC）維護(hù)的ENDF/B-Ⅵ子截面數(shù)據(jù)庫，針對B/Li截面，該數(shù)據(jù)庫與IAEA推薦的截面數(shù)據(jù)相吻合. 使用MCNP程序的F4卡對15個靈敏探測器柵元的中子平均體通量進(jìn)行統(tǒng)計，與FM4乘子卡結(jié)合來獲得相應(yīng)柵元內(nèi)發(fā)生核反應(yīng)的計數(shù)［24］. 為盡可能準(zhǔn)確地模擬得到15個探測器的中子響應(yīng)函數(shù)，采用能量為0.1 eV～20 MeV范圍內(nèi)共計177個單能中子束垂直入射能譜儀.每個能量點模擬1E9個事例，相對誤差大多低于0.1%，僅在低能中子入射時會導(dǎo)致個別位于中子射程截止處的探測器誤差較大，約為1%～10%. 利用FM4卡進(jìn)行計數(shù)，得到15個中子響應(yīng)函數(shù)結(jié)果如圖5a所示，圖中已對峰值計數(shù)率進(jìn)行了歸一化處理以便于對比.在低能區(qū)時，響應(yīng)函數(shù)同式（4）給出的曲線形狀較為一致，但高能區(qū)的差異則較大，這是由于反應(yīng)截面在高能區(qū)存在多個峰結(jié)構(gòu)（見圖1），導(dǎo)致高能區(qū)探測器的響應(yīng)曲線存在多個峰結(jié)構(gòu). 各區(qū)間的峰值計數(shù)率和FWHM如圖5b所示，峰值計數(shù)率與能量近似按1/sqrt（E）規(guī)律變化. 對比圖4b可以得到，在低能區(qū)（小于100 eV）峰值計數(shù)率和FWHM與計算值符合得較好；在中高能區(qū)（100 eV～0.1 MeV）峰值計數(shù)率逐漸偏低，F(xiàn)WHM逐漸變寬；在0.1 MeV左右時，由于截面吸收峰的影響導(dǎo)致FWHM變窄. MCNP模擬結(jié)果與式（4）和式（5）計算得到的峰值計數(shù)率和FWHM之間的平均偏差約為26.7%和57.8%，結(jié)果與預(yù)期偏差較大，詳細(xì)分析如下.

圖6中展示了公式（4）計算得到的相對寬度與MCNP模擬得到的相對寬度曲線.相對寬度R=FWHM/E，其中FWHM為中心能量為E的探測器所對應(yīng)響應(yīng)函數(shù)的半高寬. 二者在低能區(qū)基本一致;在中能區(qū)模擬結(jié)果則遠(yuǎn)大于計算結(jié)果，表明式（4）對反應(yīng)截面的近似較大低估了相對寬度；而在高能區(qū)相對寬度又迅速降低，此為模擬的半高寬在高能端趨于穩(wěn)定所致（見圖5b）. 說明截面公式對于高能區(qū)的近似以及能譜儀外增加的1 cm包裹層的貢獻(xiàn)，即公式（4）更適用于75 eV以下能區(qū)的計算. 在高能時響應(yīng)函數(shù)展寬較大.如對于峰值能量為90 eV的探測器，其響應(yīng)曲線中峰值高度一半對應(yīng)的能量分別為8.25 eV和5200 eV. 因此，雖然公式（2）在100 keV以下都可較好地描述截面，但響應(yīng)曲線對應(yīng)半高寬的誤差依然較大. 由于誤差隨慢化體厚度的增加而增大，導(dǎo)致二者偏差隨慢化體厚度增加而增大. 圖6中在高能時二者趨于接近，是由于高能區(qū)下的響應(yīng)曲線從單峰結(jié)果演化為多峰結(jié)果，導(dǎo)致FWHM計算出現(xiàn)偏差所致. 為了更準(zhǔn)確地描述響應(yīng)曲線，在中高能區(qū)需用蒙卡模擬結(jié)果進(jìn)行修正，修正方法和結(jié)果見下一節(jié).

3.3 探測器響應(yīng)曲線修正與擬合

對比圖4和圖5的響應(yīng)曲線，可以發(fā)現(xiàn)二者在低能區(qū)可較好符合，在中高能區(qū)則存在顯著差異.因此需要用蒙卡模擬結(jié)果對式（4）和式（5）進(jìn)行修正，提高計算得到響應(yīng)曲線的準(zhǔn)確性.

計算設(shè)定的15個探測器的實際峰位能量和MCNP模擬計算得到的實際峰位能量結(jié)果如圖7所示. 圖中可見，二者結(jié)果在低能區(qū)較為一致，在高能區(qū)時式（5）高估了所需的慢化體厚度.這是由于式（2）的截面公式高估了高能區(qū)的中子截面. 通過擬合MCNP結(jié)果可獲得修正后的中子能量和慢化體厚度的計算公式，使之更適合用于中子能譜儀的設(shè)計. 為更好擬合MCNP模擬結(jié)果，對圖7數(shù)據(jù)進(jìn)行分段擬合.擬合公式如式（6）所示. 擬合相關(guān)系數(shù)R2值為0.9999和0.9998，表明擬合效果較好. 基于式（6），可根據(jù)設(shè)定的探測器峰位能量來實現(xiàn)對慢化體厚度的快速設(shè)計.

EPeak=207.670D2.649+0.100，D<0.714 cm

265.852D3.529+4.206，D≥0.714 cm（6）

由于模擬和計算的響應(yīng)曲線間存在一定差異，需要對公式（4）進(jìn)行修正.同樣的，僅需要對中高能區(qū)結(jié)果進(jìn)行修正. 在能量E低于75 eV時，不需要進(jìn)行修正；在能量E高于75 eV時，需要進(jìn)行修正. 修正方法為：在公式（4）的基礎(chǔ)上增加一個高斯函數(shù)，高斯函數(shù)的系數(shù)為p， 平均值位于0.01 MeV， 標(biāo)準(zhǔn)差為1 MeV. 針對不同慢化體厚度D，修正參數(shù)有2個，分別為系數(shù)m，p. 修正后的響應(yīng)函數(shù)計算方法如式（7）所示，修正系數(shù)m， p以及4個典型響應(yīng)函數(shù)如圖8所示.

N=0.182 Ee-13.292D E，E≤75 eVm*0.182 Ee-13.292D E+Gauss（p，1E4，1E6），

E>75 eV（7）

圖8c～8f展示了其中4個探測器進(jìn)行修正后的響應(yīng)函數(shù)和MCNP模擬結(jié)果的對比圖.對于能量低于75 eV的部分未進(jìn)行修正；能量高于75 eV的部分在修正后，蒙卡模擬數(shù)據(jù)同式（7）擬合結(jié)果更為一致. 對于峰位能量在75 eV以下的探測器，擬合R2值均在0.95以上；峰位能量在75～5000 eV之間的探測器，擬合R2值從不修正時的低于0.9，提高到修正后的0.95以上；而峰位能量高于5 keV的探測器在修正后擬合效果依然較差，但擬合R2值仍較修正前有所提高.

圖8 （a， b）擬合參數(shù)m、p擬合曲線;（c～f）探測器響應(yīng)曲線擬合舉例

Fig.8 （a， b） Fitting curves of parameters m and p; （c～f） examples of detector response curve fitting

式（7）中設(shè)置75 eV的截斷是為了消除疊加的高斯函數(shù)修正項在低能端的貢獻(xiàn). 在不加截斷時，高斯修正項將導(dǎo)致低能端曲線異常偏高.通過設(shè)置75 eV的截斷可降低其對低能端曲線的影響，有助于改善擬合效果.但該截斷將導(dǎo)致修正項在75 eV處不連續(xù). 對于低能區(qū)探測器，式（7）中針對第一項的修正系數(shù)m可補(bǔ)償該截斷的影響，使截斷可忽略不計. 但對于高能區(qū)探測器，式（7）的第一項修正無法補(bǔ)償該截斷的影響，導(dǎo)致在75 eV存在明顯的截斷，使修正效果較差. 即式（7）更適合用于計算低能區(qū)（5 keV以下）探測器的響應(yīng)函數(shù)，對于高能區(qū)探測器則存在誤差，需要借助蒙卡模擬程序來提高精度.

圖8a和圖8b展示了15個探測器修正參數(shù)m、p的數(shù)值，并對修正系數(shù)數(shù)據(jù)進(jìn)行了擬合，擬合公式見式8. 兩個修正參數(shù)均可用冪函數(shù)進(jìn)行擬合，擬合的R2值分別為0.9100和0.8977.

m=11.32E0.45-2.8p=0.000202E-0.123（8）

綜上，可根據(jù)公式（6）（7）（8）進(jìn)行中子能譜儀的快速設(shè)計. 通過設(shè)定探測器可探測的能量中心值，利用公式（6）獲得其對應(yīng)的慢化體厚度，利用公式（7）（8）獲得其響應(yīng)函數(shù). 基于該方法可實現(xiàn)對中子能譜儀的快速初步設(shè)計，可快速獲得一個接近最優(yōu)方案的探測器結(jié)構(gòu)，有效避免需要大量蒙卡模擬計算才可進(jìn)行中子能譜儀設(shè)計的缺點. 但該方法給出的響應(yīng)函數(shù)存在一定的誤差，適合用于低能區(qū)探測器的響應(yīng)函數(shù)計算，尤其是中心能量在5 keV以下的探測器響應(yīng)函數(shù)計算.在計算高能區(qū)探測器的響應(yīng)函數(shù)時則存在一定的誤差，需通過蒙卡方法進(jìn)行修正. 由于本文設(shè)計的中子能譜儀針對BNCT超熱中子，不關(guān)注快中子能區(qū)，因此是否對5 keV以上能區(qū)探測器響應(yīng)矩陣進(jìn)行修正并不影響超熱中子能區(qū)的解譜結(jié)果. 研究結(jié)果表明，在直接使用公式（7）（8）計算的響應(yīng)函數(shù)且不用蒙卡結(jié)果進(jìn)行修正時，在超熱中子能區(qū)，二者解譜結(jié)果沒有差異.

4 中子能譜解析

在獲得探測器陣列的響應(yīng)矩陣R后，針對探測器陣列的輸出計數(shù)率向量N，可解譜得到入射中子能譜φ，三者關(guān)系可以用式（9）表示［25］.

N1N2Nn=R11R12…R1mR21R12…R2mRn1Rn2…Rnm×φ1φ2φm （9）

式（9）中N的每一行為單個探測器的輸出計數(shù)率，計算方法見式（10）.

Ni=∑mk=1Rikφk i=1，2，…，n.（10）

其中，n為能譜儀中探測器的數(shù)量; m為待求解能譜的分區(qū)數(shù)量. 解譜過程是針對每一組探測器計數(shù)率輸出向量N，通過響應(yīng)矩陣R求解入射中子能譜φ的過程.

對于本文設(shè)計的BNCT能譜儀，其響應(yīng)矩陣可通過修正后的公式（7）獲得，并對5 keV以上的高能區(qū)探測器使用蒙卡模擬得到的響應(yīng)函數(shù)進(jìn)行代替. 其中n=15，m=177，表示通過15個中子探測器進(jìn)行能譜測量和解析，可獲得0.1 eV～20 MeV能量范圍內(nèi)共計177個能量區(qū)間內(nèi)的中子強(qiáng)度. 因此，響應(yīng)矩陣R為15×177的矩陣.

由于各探測器響應(yīng)函數(shù)的峰值計數(shù)率存在顯著差異，尤其是高能端探測器的峰值計數(shù)率較低，導(dǎo)致其響應(yīng)靈敏度較差，從而可能較大延長能譜探測所需要的時間. 提高高能端探測器的響應(yīng)靈敏度，縮短能譜探測所需要的時間，可通過提高高能端探測器的權(quán)重來解決. 本文將探測器權(quán)重設(shè)置為各探測器響應(yīng)函數(shù)峰值的倒數(shù)，將各個探測器的響應(yīng)靈敏度調(diào)整為基本一致. 將修正后的響應(yīng)函數(shù)同時用于能譜探測和能譜解析，則對解譜結(jié)果無影響.

響應(yīng)矩陣R通常存在n<

φ（k+1）j=φ（k）jexp ∑ni=1Wkijln Ni∑ml=1Rilφkl∑ni=1Wkij ， i=1， 2，…， n；j=1， 2，…， m（11）

其中，φj為第j個能區(qū)的迭代能譜注量率；Ni為第i個探測器的實測計數(shù)率；W稱作權(quán)重因子，其計算公式見式（12）［26］.

W（k）ij=NiRijφkj∑ml=1Rilφkl（12）

Gravel的收斂條件如式（13）～（15）所示.

J（k）=∑ni=1（Ni-qki）2∑ni=1qki（13）

J（k+1）=J（k）-J（k+1）（14）

ΔJ（k+1）=J（k）-J（k+1）（15）

當(dāng)ΔJ（k+1）小于某一自定義的值時，即判斷為收斂. 其中，式（13）中q（k）i的計算公式如式（16）所示.

q（k）i=∑mj=1Rijφ（k）j（16）

對于類似Gravel之類的迭代算法，迭代精度受預(yù)置譜影響較大，而在實際應(yīng)用中的大多數(shù)情況下很難得到較為準(zhǔn)確的預(yù)置譜. 因此，本文將探測器輸出計數(shù)率向量N（1×n）與探測器陣列響應(yīng)矩陣的轉(zhuǎn)置RT（n×m）點乘得到對應(yīng)預(yù)置譜φ0（1×m） ，如式（17）所示. 結(jié)果表明，由于能譜儀各探測器的響應(yīng)函數(shù)實現(xiàn)了解耦，通過該方法設(shè)置的預(yù)置譜具有良好的普適性，可用于各種類型中子能譜的解析. 主要原因是各探測器的響應(yīng)函數(shù)之間實現(xiàn)了基本解耦.

φ0=N×RT （17）

4.1 單能中子能量解譜

通過蒙卡模擬獲得了能譜儀對單能中子輸入時的響應(yīng)情況，輸入單能中子能量分別為5E-1、5E0、5E1、5E2、5E3、5E4和5E5 eV，基本覆蓋了所有熱中子、超熱中子和快中子能區(qū). 用MCNP計算得到能譜儀的15個探測器在各種單能中子輸入時的輸出計數(shù)率向量，輸入Gravel算法進(jìn)行能譜解析，結(jié)果如圖9所示. 本文所有解譜輸出均為Gravel算法的直接輸出結(jié)果，未進(jìn)行進(jìn)一步處理，設(shè)置Gravel算法收斂條件為兩次迭代間的相對偏差小于1E-10.

解譜得到的單能中子峰與輸入譜的峰位一致程度較好，峰位偏差小于1%. 隨著輸入單能中子能量的增加，解譜帶來的展寬逐漸增加，即相對分辨率逐漸變差，從低能端的約30%逐漸增加到高能端的約90%. 在0.5 MeV單能中子輸入時，解譜結(jié)果中出現(xiàn)了兩個峰.在0.5 MeV的真實中子峰左側(cè)，出現(xiàn)了峰位在0.23 MeV的假峰. 其原因是此處響應(yīng)函數(shù)存在多個峰結(jié)構(gòu)，如圖5所示，尤其是對此能區(qū)中子最為靈敏的幾個探測器的響應(yīng)函數(shù)都存在假峰，導(dǎo)致能譜解析時出現(xiàn)假峰. 將最敏感的最后兩個探測器的響應(yīng)函數(shù)進(jìn)行平滑處理，再利用平滑后的響應(yīng)函數(shù)進(jìn)行解譜，可以避免產(chǎn)生此假峰. 確認(rèn)假峰的出現(xiàn)是響應(yīng)函數(shù)中的多峰結(jié)構(gòu)所致，解譜結(jié)果如圖9所示. 解譜中使用探測器響應(yīng)靈敏度一致的響應(yīng)矩陣，可有助于提高解譜效率和準(zhǔn)確度. 解析得到能譜的峰位、半高寬和能量分辨率如表2所示.

4.2 連續(xù)中子能譜解譜

設(shè)置了一個擁有復(fù)雜雙峰結(jié)構(gòu)的中子能譜作為輸入， 通過響應(yīng)矩陣快速獲得中子能譜儀的輸出向量，用Gravel算法進(jìn)行能譜解析，來驗證解譜算法精度，輸入能譜和解譜能譜結(jié)果如圖10所示.

由圖10可得，反解能譜與輸入雙峰能譜符合程度較好，兩個能峰的峰位基本重合，但受探測器分辨率的限制，導(dǎo)致雙峰之間波谷部分的符合程度較差，同時在高能端細(xì)節(jié)結(jié)構(gòu)上存在一定誤差. 為了準(zhǔn)確評價解譜的精確程度，分別利用輸入能譜與解析能譜之間的均方差值（MSE）表征解析譜和真實中子譜之間的差異程度， 用平均相對偏差值（ARD）表示解析譜與真實中子能譜之間的偏差程度，用中子能譜質(zhì)量值（Qs）表示解析譜與真實中子譜的接近程度.理想狀態(tài)下上述三個指標(biāo)應(yīng)接近于0［8］. 圖10中的MSE為0.896%、ARD為7.54%、Qs為7.35%，表明解譜符合程度較好，從而驗證算法的可行性.

為反解得到BNCT能譜，采用中國科學(xué)院上海應(yīng)用物理研究所的朱益楠等人［22］基于氘氚中子源經(jīng)BSA慢化整形后的BNCT中子能譜作為MCNP的輸入. 該中子能譜符合IAEA推薦的BNCT治療用中子束指標(biāo)［27］. 模擬獲得了能譜儀探測器陣列的輸出計數(shù)率向量，通過GRAVEL算法進(jìn)行解譜，結(jié)果如圖11所示.

圖中可見，GRAVEL算法給出的能譜和輸入能譜符合較好，可較好的重現(xiàn)BNCT能譜結(jié)構(gòu). 僅在0.1～20 MeV的高能范圍內(nèi)存在差異，原因是該能譜儀的15個中子探測器均是針對超熱中子能區(qū)的，對快中子響應(yīng)靈敏度較差，限制了在該區(qū)域的靈敏度. 此外，10B（n，α）7Li反應(yīng)截面在該范圍內(nèi)存在多個峰結(jié)構(gòu)，也導(dǎo)致解析難度的增加. 增加一個對快中子響應(yīng)敏感的探測器可有效提高在快中子能區(qū)的測量精度. 由于BNCT能譜主要位于超熱中子能區(qū)，快中子影響相對較小，因此本文設(shè)計的能譜儀主要采用對超熱中子能區(qū)敏感的探測器，其體積更小，具有便攜性和較好的靈敏度.

圖11中的MSE為0.76%、ARD為8.11%、Qs為6.67%，說明解譜的精度較高. 結(jié)果表明，GRAVEL算法反解得到的中子能譜與輸入譜符合度較高，說明該能譜儀適用于BNCT的能譜測量和解析. 對比中國原子能科學(xué)研究院的陳軍等人［28］利用多球譜儀求解的BNCT能譜，本文的中子能譜儀反解得到的能譜精度更高，峰的展寬更小，從而驗證了本結(jié)構(gòu)的優(yōu)越性.

5 結(jié) 論

本文提出了一種使用B/Li作為慢化體和靈敏探測器的通用中子能譜儀快速設(shè)計方法.此方法可通過公式快速實現(xiàn)能譜儀最優(yōu)慢化體厚度和中子探測器響應(yīng)函數(shù)的計算， 結(jié)果同蒙特卡羅模擬結(jié)果較為一致. 該類能譜儀各個探測器的響應(yīng)函數(shù)可實現(xiàn)解耦，在預(yù)置譜設(shè)置、解譜精度等方面具有較大的優(yōu)勢. 利用B4C作為中子探測器的慢化吸收體，設(shè)計了一個針對BNCT超熱中子探測的能譜儀，利用少道解譜算法GRAVEL算法實現(xiàn)了能譜解析. 解譜結(jié)果表明， 該能譜儀具有較高的中子能譜測量精度，可以用于BNCT超熱中子的能譜探測，驗證了該方法的可行性.

參考文獻(xiàn)：

［1］ 蔣世倫， 祁建敏， 王立宗， 等. 反沖質(zhì)子磁分析技術(shù)用于氘氚中子能譜測量研究［J］. 物理學(xué)報，2011， 61： 072902.

［2］ 祁建敏， 周林， 蔣世倫， 等.聚變中子能譜測量系統(tǒng)脈沖中子靈敏度的實驗研究［J］. 物理學(xué)報， 2013，62： 245203.

［3］ 陳永浩. 反沖質(zhì)子法快中子能譜測量及解譜技術(shù)研究［D］. 蘭州： 蘭州大學(xué)， 2013.

［4］ Esposito A， Bedogni R， Domingo C， et al. Measurements of leakage neutron spectra from a high-energy accumulation ring using extended range bonner sphere spectrometers ［J］. Radiat Meas， 2010， 45： 1522.

［5］ Goldhagen P， Reginatto M， Kniss T， et al. Measurement of the energy spectrum of cosmic-ray induced neutrons aboard an ER-2 high-altitude airplane ［J］. Nucl Instrum Meth A， 2002， 476： 42.

［6］ Vega-Carrillo H R， Manzanares-Acua E. Background neutron spectrum at 2420 m above sea level ［J］. Nucl Instrum Meth A， 2004， 524： 146.

［7］ Wiegel B， Agosteo S， Bedogni R， et al. Intercomparison of radiation protection devices in a high-energy stray neutron field， Part II： bonner sphere spectrometry ［J］. Radiat Meas， 2009， 44： 660.

［8］ 程凱， 魏鑫， 曾德凱， 等. 基于微結(jié)構(gòu)氣體探測器對單能和連續(xù)譜快中子的模擬解譜［J］.物理學(xué)報， 2021， 70： 112901.

［9］ 鄒益晟， 肖德濤， 張偉華， 等. 基于單球中子譜儀的能量響應(yīng)模擬及解譜方法研究［J］.輻射防護(hù)， 2019， 39：? 280.

［10］ Birattari C， Esposito A， Ferrari A， et al. Calibration of the neutron rem counter LINUS in the energy range from thermal to 19 MeV ［J］. Nucl Instrum Meth A， 1993， 324： 232.

［11］ Birattari C， Ferrari A， Nuccetelli C， et al. An extended range neutron rem counter ［J］. Nucl Instrum Meth A， 1990， 29： 250.

［12］ Tolymkhan Y， 王忠海， 陳秀蓮， 等. 以水作為慢化體的多球中子譜儀模擬與解譜研究［J］. 四川大學(xué)學(xué)報： 自然科學(xué)版， 2020， 57： 735.

［13］ Takada M， Abe Y， Nakamura S， et al. Spectrometer design of low energy neutrons for boron neutron capture therapy ［J］. Nucl Instrum Meth A， 2021， 1020： 165848.

［14］ Hosseini S A， Alvarez-Galvan M C. Study of physical-chemical properties and catalytic activities of ZnCr2O4 spinel nano oxides obtained from different methods—modeling the synthesis process by response surface methodology and optimization by genetic algorithm ［J］. J Taiwan Inst Chem E， 2016， 61： 261.

［15］ 莫雙榮， 劉鈺， 幸浩洋， 等. Elman 神經(jīng)網(wǎng)絡(luò)在中子解譜中的應(yīng)用［J］. 四川大學(xué)學(xué)報： 自然科學(xué)版，2020， 57： 531.

［16］ Pehlivanovic B， Avdic S， Marinkovic P， et al. Comparison of unfolding approaches for monoenergetic and continuous fast-neutron energy spectra ［J］. Radiat Meas， 2013， 49： 109.

［17］ 王松林， 王琦， 徐小三， 等. 閾探測器法測量 Am-Be 中子源屏蔽輻照腔內(nèi)的中子能譜 ［J］. 原子能科學(xué)技術(shù)， 2009， 43： 16.

［18］ Suman V， Tripathy S P， Sunil C， et al. Measurement of neutron energy distributions from p+Be reaction at 20 MeV using threshold activation foils ［J］. IEEE T Nucl Sci， 2016， 63： 2283.

［19］ Wang G Y， Han R， Liu J L. Comparison and research on the GRAVEL and PRIP algorithms of neutron energy spectrum unfolding ［J］. Radiat Detect Techno， 2017， 1： 5.

［20］ 張馳， 王玉東， 周榮， 等. G（E）法與Gravel法處理能譜劑量轉(zhuǎn)換效果研究［J］. 核電子學(xué)與探測技術(shù)，2017， 37： 268.

［21］ 李建偉， 李德源， 劉建忠， 等. 三種解譜算法求解中子能譜的解譜效果比較［J］. 核電子學(xué)與探測技術(shù)， 2017， 37： 147.

［22］ 朱益楠，林作康，盧林遠(yuǎn)，等. 基于氘氚中子源硼中子俘獲治療的中子慢化整形研究［J］.核技術(shù)， 2022， 45： 010202.

［23］ 盧希庭. 原子核物理（修訂版）［M］. 北京：原子能出版社， 2000： 242.

［24］ 燕奕宏， 譚新建， 翁秀峰， 等. 三氟化硼B(yǎng)onner多球譜儀的中子源能譜測量方法［J］. 西安交通大學(xué)學(xué)報， 2020， 54： 136.

［25］ Banerjee K， Ghosh T K， Kundu S， et al. Variation of neutron detection characteristics with a dimension of BC501A neutron detector ［J］. Nucl Instrum Meth A， 2009， 608： 440.

［26］ Liu B， Lv H， Li L， et al. Study on iterative regularization method and application to neutron spectrum unfolding of multi-sphere spectrometer measurement ［J］. Nucl Instrum Meth A， 2021， 992： 165027.

［27］ Levin V A P， Dodd B. Current status of neutron capture therapy ［R］. Vienna， Austria： IAEA， 2001.

［28］ 陳軍， 李春娟， 李瑋， 等. 多球譜儀測量BNCT醫(yī)院中子照射器中子束能譜［J］. 原子能科學(xué)技術(shù)，2014， 48： 2127.

引用本文格式：