基于微結構氣體探測器對單能和連續譜快中子的模擬解譜*

程凱 魏鑫 曾德凱 季選韜 朱坤 王曉冬

(南華大學核科學技術學院， 衡陽 421000)

本工作探索了基于Triple GEM探測器對快中子能譜的測量， 利用MCNPX和Geant4軟件分別模擬了兩種在Triple GEM陰極耦合由多層聚乙烯組成的堆棧式中子轉化質子的轉化模型， 研究對象包含了5種單能中子源和Am-Be連續譜中子源.模擬得到了探測系統對160條單能中子的響應函數和上述源的反沖質子譜分布， 使用GRAVEL和MLEM算法對模擬得到的6種快中子源的反沖質子譜進行了解譜研究， 并把解譜結果與標準輸入譜進行了對比， 結果顯示與標準輸入譜均符合較好， 解譜的相對不確定度為10%—15%; 并研究了氣體探測器能量分辨率對解譜精度影響的關系， 結果表明微結構氣體探測器的能量分辨率好于30%時， 快中子的解譜精度就可以滿足實際應用需求.本研究在以前的實驗基礎上提出了一種新的轉化結構，并通過模擬結果得出微結構氣體探測器可以應用于快中子探測， 并能夠利用得到的反沖質子譜結合合適的反演算法實現入射中子源的能譜重建.本文積累的建模和解譜算法為將來微結構氣體探測器組成的快中子探測系統應用于未知快中子源探測， 能譜重建， 實現源項識別提供了新的辦法.

1 引 言

中子能譜[1，2]是每個中子源的獨有特征， 可用于源項的識別[3]， 精準獲得中子能譜信息在空間探測、輻射防護、核不擴散和國土安全、核材料管制等領域具有重要實際意義[4，5].近十年， 國內外對中子能譜的測量進行了大量的研究.Avdic等[6]基于液體閃爍體探測器測量了反沖質子的脈沖高度譜， 使用了非線性人工智能算法進行了解譜， 通過模擬和實驗驗證， 達到了識別未知中子源的目的.Hosseini等[7]使用人工神經網絡(artificial neural network， ANN)算法和改進的最小二乘法(modified least square， MLSQR)對NE-213 閃爍體探測器收集的反沖質子能譜進行了中子解譜， 論證了ANN 算法對中子解譜的適用性.Pehlivanovic等[8]使用NE-213液體閃爍體進行了中子能譜探測研究， 使用了最大似然估計(maximum-likelihood expectation maximization， MLEM)算法和一步延時法(one step later， OSL)對反沖質子的脈沖高度譜進行了解譜， 結果表明對于連續中子源， OSL算法解譜效果優于MLEM算法.國內在快中子解譜方面也做了一系列研究工作， 并取得了非常好的成果， 如王冠鷹等[9]利用液體閃爍體探測器和金硅面壘探測器獲得反沖質子， 并使用勢能還原內點算法對中子能譜進行解譜， 該算法能夠快速準確地求解中子能譜.言杰等[10]使用簡單的不含隱藏層的線性神經網絡結構， 對BC501A液體閃爍體探測器測的反沖質子譜進行了解譜， 但由于受響應函數的限制， 最后解得的能譜能點較少， 能譜數據較為稀疏， 能譜不光滑.王冬等[11]基于遺傳算法對多球中子譜儀開展了中子解譜的仿真測試， 提出了一種由距離項和懲罰項組成的新型適應度函數， 驗證了基于該新型適應度函數的遺傳算法的可行性.林存寶[12]使用了信賴域算法對 BC501A液體閃爍體探測器測的反沖質子譜進行解譜， 得到了Am-Be， D-T和252Cf 等中子源的中子能譜.燕奕宏等[13]使用Bonner球譜儀測量中子源能譜， 并使用GRAVEL算法進行解譜， 中能區解譜結果與輸入譜符合較好， 低能區和高能區符合較差.

綜上， 當前用于中子解譜的探測器多為閃爍體探測器和Bonner球， 氣體探測器應用于此方面較少， 主要是因為氣體介質的能量阻止本領較低， 中子轉化為質子的效率低， 導致探測效率低， 不利于中子探測.

本工作使用了一種微孔狀的氣體探測器-氣體電子倍增器(gas electron multiplier， GEM)， 1997年由歐洲核子中心CERN的Sauli[14]研制， 具有成本低廉、耐輻照、增益高(102—106)、空間分辨率好(70 μm)[15，16]等諸多優點， 且可以通過多層GEM膜級聯的方式提高信噪比， 將GEM探測器改裝成一種中子敏感的探測器， 操作比較簡單.近年來國際上雖然有多個科研機構針對基于GEM結構的中子測量進行了研究， 也取得了很好的成果[17-20]，但大多數集中在使用含硼的轉化材料提高熱中子探測效率方面的研究.本文的研究工作是利用Triple GEM探測器耦合一種新的快中子轉化結構， 其將快中子轉化為可直接探測的反沖質子， 而氣體探測器對質子的探測效率近乎于100%， 因此可以用反沖質子數的產額表征入射中子的探測效率.利用MCNPX和Geant4分別模擬了基于Triple GEM設計的中子探測器對不同中子源的響應， 獲得了由160個單能中子響應函數組成的探測器響應矩陣， 以及探測器對Am-Be中子源的反沖質子能譜.由于響應矩陣通常情況下為病態矩陣，因此通過反沖質子譜求解入射中子譜(解譜)的過程實際上是求解不適定問題， 不能用直接求解的方式， 需要增加一些條件獲得其近似解.國內外研究團隊發展了GRAVEL算法[21]、最小二乘法[22]、極大似然算法[23](MLEM)、一步延時算法、信賴域算法和遺傳算法等多種解譜算法.本文使用GRAVEL算法和MLEM算法分別對以上模擬數據進行了聯合求解， 得到了入射中子源的能譜， 并和標準能譜進行了比較.

2 探測系統

為了實現快中子高效率的探測， 本文提出了一個堆棧式高密度聚乙烯(high density polyethylene， HDPE)快中子轉化結構， 原理圖如圖1所示.快中子垂直于聚乙烯薄膜進入探測器， 與其發生彈性碰撞得到反沖質子， 能量大的反沖質子逃出聚乙烯表面， 進入到氣體間隔內和氣體原子或分子通過碰撞電離產生電離電子， 這些電離電子通過漂移區電場的引導進入到微孔中發生雪崩效應， 最后受收集區電場的作用， 在向陽極移動過程中使陽極板產生了感應信號， 即反沖質子能量信息.為了準確求解入射中子能譜， 探測器的響應矩陣是一個關鍵因素， 能否準確地獲得響應矩陣將直接決定中子解譜的正確性和準確性， 由于單能中子源非常稀缺， 因此需要模擬大量的單能中子與轉化結構相互作用得到響應函數， 不同能量的單能中子通過探測器后得到的反沖質子數是一一對應的， 分別對應于探測器響應矩陣中的一個元素.

圖1 中子進入Triple GEM探測器原理圖(藍色為HDPE，綠色部分為CO2和Ar混合氣體)Fig.1.Schematic diagram of the Triple GEM-based neutron detector (The blue color is HDPE， and the green part is the mixture of CO2 and Ar).

2.1 堆棧式多層快中子轉化結構模擬

使用MCNPX和Geant4對聚乙烯轉化結構進行模擬， 主要參數為: 入射中子能量范圍2—17 MeV; 入射中子數為108個; 中子源距離聚乙烯薄膜距離為4 cm; 聚乙烯薄膜的密度為0.95 g/cm3.

由于快中子與聚乙烯轉化結構中的氫元素發生相互作用轉化為質子時會損失部分能量， 與之前入射中子相比能量減少， 若采用同一厚度聚乙烯進行轉化， 可能會導致這些中子產生的質子因能量較小而不能逃出聚乙烯.圖2所示為單能中子探測效率與最優厚度之間的關系.

圖2 不同能量的單能中子探測效率與聚乙烯厚度之間的關系及最優厚度Fig.2.The function of the detector efficiency with different polyethylene thicknesses.

圖3 所示為MCNPX軟件對中子與聚乙烯發生彈性碰撞過程建模圖， 一部分反沖質子如圖中藍色徑跡所示， 其穿過聚乙烯后進入氣隙內， 能量沉積在氣隙間隔; 圖4為MCNPX和Geant4軟件模擬DT， Am-Be， DD源的探測效率的對比.

由圖4可知， 當聚乙烯厚度較薄時， 探測效率隨著厚度的增加而增加; 當厚度達到一定時， 探測效率達到最大值， 此時的厚度為最優厚度.當厚度超過最優厚度， 探測效率趨于飽和.實際上不管多厚的聚乙烯膜， 反沖質子僅能從聚乙烯膜表面幾十個μm出射.Geant4和MCNPX對三種不同中子源與聚乙烯相互作用的最佳厚度和探測效率的模擬結果對比如表1所列.

由于Geant4和MCNPX兩者物理過程中使用的截面庫和中子通量的計算方式不同， 因此探測效率會有一些差異[24.25].尤其是在能量較低的區域， 探測效率相差較大， Geant4比MCNPX的探測效率在DD源處高133%， Am-Be源處高33%，DT源處高12.5%.

圖3 中子與單層聚乙烯彈性散射圖探(紅色點代表中子，藍色為質子)Fig.3.Screenshot of the elastic collision between neutrons and polyethylene in MCNPX (the red dots represent neutrons and the blue are protons).

圖4 MCNPX與Geant4對不同中子源的測效率的模擬對比Fig.4.Simulation comparison of detection efficiency between MCNPX and Geant4 for the same thickness of polyethylene.

表1 Geant4和MCNPX對不同中子源與聚乙烯相互作用的最佳厚度和探測效率的模擬結果Table 1.Simulation results of Geant4 and MCNPX for the optimal thickness and detection efficiency of different neutron sources interacting with polyethylene.

由于單層聚乙烯與中子發生彈性散射幾率很低， 絕大多數中子會穿過聚乙烯薄膜， 為了提高探測效率， 采用多層聚乙烯結構， 以提高中子發生彈性散射的概率.將最優厚度聚乙烯和一定厚度氣隙作為一個轉化單元， 根據文獻[18]可知， 氣體間隔采用500 μm較為合理.中子與多層聚乙烯轉化結構發生相互作用如圖5所示， 反沖質子出現了小角度散射， 反沖質子徑跡與中子源在同一水平線.單元數與轉化效率的關系如圖6所示， 模擬結果表明DT， Am-Be和DD源的探測效率最大值分別為1.50%， 0.44%和0.14%， 此時聚乙烯轉化單元均為200層左右， 比單層的探測效率高近5倍， 由此可見， 采用了多層聚乙烯的方法可以有效提高探測效率.

圖5 MCNPX模擬的中子與多層聚乙烯相互作用(紅色的點代表中子， 藍色為質子)Fig.5.MCNPX simulated neutron interactions with multilayered polyethylene (the red dots for neutrons， the blue for protons).

圖6 多層聚乙烯結構轉化效率Fig.6.Conversion efficiency of multilayer polyethylene structures.

使用MCNPX對多層聚乙烯轉化結構進行了模擬， 由于初始入射中子與聚乙烯發生彈性散射產生反沖質子后會損失了一部分能量， 以致后續的聚乙烯厚度和層數不是最優， 為了獲得更多反沖質子， 采用聚乙烯厚度遞減的堆棧式的優化結構.第一種模型有17個堆棧， 每個堆棧內有20層相同厚度的聚乙烯， 考慮到3D打印技術加工精度為100 μm，因此將相鄰堆棧內聚乙烯厚度按100 μm的厚度遞減， 與Geant4軟件模擬的結構相同[26]， 分別用M1和G1表示(下同).為了提高反沖質子的產額，對此提出了新的聚乙烯轉化結構， 第二種模型采用了能量截斷的方法， 將2—17 MeV的入射中子劃分為15個間隔， 每個能量間隔為1 MeV， 求出每個能量截斷值所對應的最優聚乙烯厚度及穿過的最大層數， 將每個能量間隔穿過的最大層數作為一個堆棧， 這15個堆棧組成了新的聚乙烯轉化結構，用M2表示(下同).如圖7所示， 為M2模擬的部分聚乙烯堆棧結構示意圖.

圖7 M2模擬的堆棧結構示意圖(藍色、綠色、橘黃色和紫色分別代表不同厚度聚乙烯， 紅色代表氣隙厚度)Fig.7.Schematic diagram of the stack structure simulated by M2 (blue， green， orange， and purple represent different thicknesses of polyethylene respectively， red represents air gap thickness).

2.2 響應函數的計算

探測器響應函數描述了入射中子能量與其引起的脈沖幅值之間的隨機關系， 需要根據具體的實驗條件來確定.探測器的響應函數取決于許多可變的因素， 主要包括中子源特性以及探測器的材料組成、幾何特性、工作狀態和計數率等.

利用MCNPX軟件對GEM探測器的響應矩陣進行模擬計算.響應矩陣由入射的單能中子和聚乙烯轉化結構相互作用后得到的反沖質子能譜組成， 將每一個入射的單能中子與其相對應的反沖質子作為一個行向量將這些向量存入一個矩陣中一行， 該矩陣即為n條行向量組成的響應函數矩陣，共160條.通過對響應矩陣和反沖質子譜聯合求解即可得出入射中子源能譜信息， 其關系如下所示.

響應矩陣:

反沖質子能譜:

反沖質子譜g與響應矩陣A之間的關系可以寫為下式的形式:

其中f即為所求入射中子譜的能譜信息.

圖8僅展示了16條整數單能中子的響應函數.使用的入射中子能量為2—17 MeV， 入射中子能量間隔為0.1 MeV， 因此得到了 1 60×170 的響應矩陣.由于模擬單能中子的能量間隔小， 在解譜時便于區間劃分， 得到的響應函數越精細， 求解得出的中子譜的點越密集.

圖8 不同中子源對應的響應函數Fig.8.Response functions corresponding to different neutron energies.

3 解譜算法

3.1 GRAVEL解譜算法

GRAVEL解譜算法由SAND-Ⅱ算法演化而來， 是最早由德國 PTB 實驗室提出的一種交互式迭代算法， SAND-Ⅱ常用于反沖質子解譜， 其迭代公式為

式中，fj是入射中子在j能量區間的輻射強度;fk為第k次迭代得到的第j個能量間隔的輻射強度，當k為0時， 初始中子譜常取常數譜;Aij是第i個脈沖幅度間隔與第j個能量區間的耦合的響應函數;是權重因子;gi是第i個脈沖幅度處的強度，σi為gi的統計誤差， 一般為gi的平方根， 若gi的數值為0， 則σi取1.GRAVEL具有對初始值不敏感， 解為非負值且噪聲小等優點， 缺點是解不完全收斂， 迭代次數越多， 越容易發散.

3.2 MLEM算法

MLEM算法是利用樣本結果信息， 反推最大概率導致這些樣本結果出現的模型參數值.理論上， 質子和中子譜均為連續函數， 測量裝置允許在有限的多個點測量質子能譜的值.因此將實軸劃為有限的多個能量區間， 區間越小， 落在區間的粒子數具備泊松分布， 且相互獨立.

MLEM算法具有以下特點: 算法收斂， 對于每個k， 總有l(fk)＜l(fk+1).若初始值f0非負， 那么整個響應矩陣A和質子譜g均是非負， 所有的解fk也都為非負.對于每個近似的fk總有

由此可見， 脈沖總和保持不變， 只在不同能量組中重新分配.

4 結果分析

4.1 單能中子源解譜

通過對單能入射中子得到的反沖質子譜的電離過程進行模擬， 得到反沖質子在氣隙(氣體成分為70% Ar和30% CO2)中的沉積能量， 如圖9所示.隨著入射中子能量增加， 反沖質子最大概率沉積的能量(峰位所對應的能量)逐漸減少， 這主要是由于氣隙間隔固定(500 μm)， 遠小于質子在氣體中完全沉積時所穿行的距離， 因此低能質子更容易沉積到氣隙里.圖10顯示了入射中子能量與反沖質子沉積能量峰位的關系.

圖9 不同入射中子對應的反沖質子在氣隙中的沉積能量Fig.9.Deposition energy of recoil protons in the air gap corresponding to different incident neutrons.

圖11 為使用兩種算法分別對Geant4軟件和MCNPX軟件的模擬數據進行解譜的結果.將沉積能量與響應矩陣聯合求解得到單能入射中子源的能譜信息.從圖11中可知， 兩種算法都正確地預測了能量分辨率約為0.1 MeV的中子的解譜能量，但其峰值強度高低不等， 以峰值強度最大值為標準進行歸一化處理.在4 MeV處的峰值與參考數據相匹配， 但在2.4， 8.0， 12.0和14.0 MeV處的峰值分別低估了32%， 15%， 5%和20%， 這因為在解譜算法中， 在峰值處出現了不同程度的能譜展寬.

圖10 入射中子能量與反沖質子沉積能量關系Fig.10.Relationship between neutron incidence energy and recoil proton deposition energy.

4.2 連續中子源解譜

當入射源為Am-Be源時， 通過GRAVEL算法和MLEM算法對模擬的探測器的響應矩陣和脈沖幅度譜進行解譜， 解譜的結果如下圖12所示.

從圖12可以看出， GRAVEL算法和MLEM算法均可以用于Am-Be源的解譜， 兩種算法解譜的峰位均與輸入譜符合較好.

圖11 兩種算法對Geant4軟件和MCNPX軟件模擬數據的解譜結果 (a) GRAVEL算法對單能中子解譜圖; (b) MLEM算法對單能中子解譜圖Fig.11.Results of two algorithms for solving the spectrum of simulated data from Geant4 software and MCNPX software:(a) GRAVEL algorithm on single-energy neutron unfolding spectra; (b) MLEM algorithm for single-energy neutron unfolding spectra.

圖12 兩種解譜算法對Geant4軟件和MCNPX軟件模擬Am-Be源的解譜結果 (a) MCNPX與Geant4基于GRAVEL算法的解譜結果; (b) MCNPX與Geant4基于MLEM算法的解譜結果Fig.12.Results of the unfolding of the simulated Am-Be sources by two solution algorithms for Geant4 software and MCNPX software: (a) Unfolding results of MCNPX and Geant4 based on GRAVEL algorithm; (b) unfolding results of MCNPX and Geant4 based on MLEM algorithm.

G1結構解譜結果表明， GRAVEL算法在低能區出現較大幅度振蕩， 在中高能區振蕩幅度較小， 這是由于GRAVEL算法的解不完全收斂， 因此會出現不同程度的能量展寬.MLEM算法則比較平滑， 可見MLEM算法抗干擾能力強于GRAVEL算法.MLEM算法與輸入譜相比， 整體符合較好.

對M1結構和M2結構解譜， 兩種算法在2.7—11.0 MeV與標準譜符合較好， 但在1.0—2.6 MeV處的值低于輸入中子譜在該范圍內的值， 由此可見， MCNPX更適合模擬中高能區中子.

使用氣體探測器對55Fe低能的5.9 keV的X射線測量的能量分辨率的好于15%.現將探測器應用于反沖質子的測量， 通過模擬沉積能量峰值為13.8 keV的反沖質子在探測器內部電離作用以及電離電子在GEM孔內雪崩效應， 得到了探測器的能量分辨率為13.4%， 如圖13所示， 相比于對X射線的測量， 能量分辨率得到改善和提高， 因此Triple GEM結構的中子探測器在中子探測及測量方面應用是可行的.

為了使基于氣體探測器的中子探測系統能給將來的實驗提供更可靠的參考， 本工作在模擬仿真時考慮了氣體探測器的能量分辨率對解譜精度的影響， 模擬時采用的氣體探測器的能量分辨率為10%—30%， 如圖14所示， 結果表明探測器的能量分辨率越小， 求解的中子譜的相對不確定度就越小， 結果越準確.同為MCNPX軟件模擬的兩種不同模型之間， 相對不確定度隨著能量分辨率變差而逐漸增大， 且趨勢大致相同.Geant4模擬結果解譜的相對不確定度在能量分辨率為10%—18%時為均勻增加， 在能量分辨率大于18%時， 相對不確定度擴大了1%.兩種不同的變化趨勢原因在于兩者模擬結果數據的不同， Geant4模擬結果為單個粒子數分別耦合能量分辨率的高斯分布， 然后對其進行不同能量間隔的劃分， 因此隨著能量分辨率不斷變差， 能譜的展寬程度愈加劇烈， 導致相對不確定度不是均勻變化.而MCNPX模擬結果則先進行能量區間劃分， 然后再耦合關于能量分辨率的高斯分布， 因此， MCNPX的相對不確定度為均勻變化.對于同一軟件或模型， GRAVEL算法和MLEM算法解譜結果的相對不確定度相差不大， 約為0.2%.

圖13 反沖質子在氣體探測器中的能量分辨率Fig.13.Energy resolution of recoil protons in gas detectors.

圖14 相對不確定度與能量分辨率的關系Fig.14.Plot of relative uncertainty versus energy resolution.

表2 MCNPX兩種結構和Geant4模擬結果經MLEM算法和GRAVEL算法解譜后誤差水平Table 2.Error levels of the two MCNPX structures and Geant4 simulation results after unfolding by the MLEM and GRAVEL algorithms.

為了準確地評價兩種算法對MCNPX的兩種結構和Geant4的結構模擬結果的求解效果， 分別用均方誤差(mean square error， MSE)表征預測的中子譜和真實中子譜之間的差異程度， 理想狀態下接近于0， 平均相對偏差(average relative deviation， ARD)表示預測的中子光譜與真實中子能譜之間的偏差程度， 中子能譜質量(quality of neutron spectrum， Qs)表示預測的中子譜與真實中子譜的接近程度.作為評價解譜結果誤差水平的標準， 公式如下所示:

其中fj為所求的入射中子譜解譜結果;為標準中子源.對MCNPX兩種結構和Geant4模擬的結果進行誤差分析， 并列于表2中.

由表2可得G1比M1的解譜結果在MSE，ARD和Qs三個誤差方面分別好1.7%， 3%和16%.造成這種誤差的原因主要在于MCNPX在低能區1.0—2.6 MeV模擬結果較差; M2比M1的解譜結果在MSE， ARD和Qs三個誤差方面分別好0.7%，1%和9%.由于兩種結構響應函數不同， M2在低能區可以獲得更多的反沖質子， 因而對于低能區解譜結果略好于結構M1.比較GRAVEL算法與MLEM算法的解譜結果， GRAVEL算法比MLEM算法在MSE， ARD和Qs分別高約0.3%， 1.0%和2.0%， 誤差是由于GRAVEL算法解譜結果在低能區存在較大振蕩， MLEM算法解譜結果相對平滑所致.

5 結 論

本工作利用MCNPX和Geant4軟件對基于Triple GEM結構快中子探測器陰極轉化結構進行了研究模擬， 獲得了該轉化結構對應的響應函數矩陣， 由于響應矩陣為病態矩陣， 因此響應矩陣和脈沖幅度譜之間微小的變化都會導致解譜結果和輸入譜有很大偏差.通過利用GRAVEL算法、MLEM算法對M1， G1和M2的模擬數據進行了解譜.解譜結果表明， 兩種算法在Triple GEM結構快中子探測器的解譜過程中均有很好的效果.GRAVEL算法在求解低能部分時存在解譜精度低， 結果不穩定， 有較大幅度的振蕩.MLEM算法解譜精度相較于GRAVEL算法高， 且解譜曲線更加平滑， 與輸入譜符合較好， 是一個實用的中子解譜方法.兩種蒙特卡羅軟件， 對于單能中子源模擬均比較準確，可以用于對未知中子源的模擬.對于連續中子源，Geant4的低能部分精確度相較于MCNPX要高.MCNPX軟件更適合用于中高能中子的模擬.本工作完成了MCNPX對于探測器結構的模擬與中子能譜解譜的計算， 后續將展開相關實驗， 進一步探究探測器轉化結構對于解譜的精度的影響， 并考慮低能部分用G4模擬， 中高能部分用MCNPX模擬.