宇宙射線繆子核材料快速檢測算法研究

錢祎劍，張立軍，陳靈新，王冠鷹,*

(1.中國科學院 微電子研究所，北京 100029；2.中國科學院大學 微電子學院，北京 100190)

核材料走私威脅著世界和平發展，有效監管核材料、打擊核走私一直是世界各國的關注重點。宇宙射線繆子具有穿透性強、無額外輻照等特點，相比γ射線等傳統射線探測技術，能穿透Pb等屏蔽材料，不會破壞待測物質和危害人員安全，是核材料的理想檢測技術。2003年，美國LANL根據繆子在物質中的庫倫散射特性，提出用于核材料檢測的PoCA[1](point of closest approach)算法，成功對兩條鋼軌和鎢塊進行成像。隨后LANL基于極大似然估計開發了成像精度更高的MLSD(maximum likelihood scattering and displacement)[2-4]算法，主要應用于對核反應堆和核燃料存儲狀態的監測[5-6]。對這兩類算法的改進研究[7-9]主要以高精度成像為目的，為保證成像區域有足夠的繆子探測數據進行觀測，需較長的探測時間，在港口、鐵路等貨運量大，需快速檢出待測物是否為核材料的場合并不適用。近幾年，為提高核材料檢測速度，不通過成像來判別是否存在核材料的宇宙射線繆子快速檢測算法被提出[10-13]。2019年，中國科學技術大學的何偉波用灰色關聯分析[13]實現了僅用繆子散射角數據即可對U、Pb、Fe進行分類，2 min的探測時間對U和Pb的區分可達到95%以上的準確率。

為進一步提高用宇宙射線繆子檢測核材料的速度，研究采用較少的繆子探測數據實現核材料檢測，本文利用Geant4[14]軟件仿真繆子庫倫散射數據進行分析驗證，研究繆子探測數據的概率分布函數和高階統計量，使用支持向量機測試這些參數對核材料和檢測干擾物的分類性能，并提出一種新的宇宙射線繆子核材料快速檢測算法。

1 原理

使用宇宙射線繆子對核材料進行檢測成像是根據繆子在材料中的庫倫散射特性實現的。宇宙射線繆子產生于宇宙射線與地球大氣的粒子反應，在地球海平面處平均通量約為10 000 min-1·m-2。繆子穿過物體時會發生能量損失和庫倫散射，其中繆子在材料中多次庫倫散射后偏轉角度近似于零均值的高斯分布，其平面角和立體角時的分布均方根RMS(mrad)分別為式(1)[1]和式(2)[2]。

(1)

(2)

其中：βcp為繆子入射能量，MeV；x為繆子穿過物體時的路徑長度；X0為物體材料的輻射長度(g·cm-2)，其定義[2]為：

(3)

其中，Z、A分別為材料的原子序數和相對原子質量。式(2)在10-3>x/X0>102時能較好描述繆子穿過材料的庫倫散射角分布，其誤差好于11%[15]。

式(1)～(3)表明，繆子在穿過物體時庫倫散射角分布與物體材料原子序數有關。通過統計繆子探測數據，計算出穿過待測物體的繆子庫倫散射角分布均方根值，即可據此判斷空間內是否存在核材料。檢測算法所需的繆子探測數據量越少，對應所需的探測時間就越少，核材料檢測速度就越快。

2 仿真環境的建立和數據獲取

為快速獲取大量繆子穿過不同材料的散射數據，進一步開展宇宙射線繆子核材料快速檢測算法研究，使用Geant4軟件包獲取繆子穿過待測材料時的庫倫散射角仿真數據。Geant4仿真環境設置如下：在空氣環境中以坐標原點為中心放置一塊厚度10 cm、底面為100 cm×100 cm的待測材料，在待測材料中心上、下50 cm處放置3層間隔5 cm、厚度為1cm、面積為100 cm×100 cm的塑料閃爍體(材質為聚苯乙烯，密度為1.05 g/cm3，H/C比為1.12)作為繆子探測器。仿真環境如圖1所示，處于上、下兩組繆子探測器中間部分的1 m3立方體為探測空間。通過記錄繆子穿過上、下兩組探測器的位置，分別擬合入射和出射直線，計算其夾角即可得到繆子穿過探測空間時的散射角。

式(2)描述的繆子庫倫散射角分布均方根包含繆子入射能量，對于不同入射能量、入射角度的繆子散射數據，忽略式(2)中較小的對數項，可將其歸一化得到相同入射能量、垂直入射繆子的散射角：

(4)

其中：pin為繆子入射動量；θin為繆子入射角度；θs為繆子穿過探測空間時的累積庫倫散射角。式(4)是歸一化入射天然源繆子的有效方式，參照已有研究[4,13]，假設探測器的接收效率為100%，在Geant4仿真環境中，直接采用能量為3 GeV、垂直入射的繆子點源代替天然繆子源項，位于最上層探測器中心上方5 cm處，點源與探測器之間的空間為真空環境。該設置簡化了天然源繆子的數據處理流程，利于快速獲取不同材料中的繆子庫倫散射數據，開展對繆子散射探測數據分布的研究。

圖1 仿真實驗環境示圖Fig.1 Illustration of simulation experiment environment

選擇元素U(Z=92)作為被檢測的核材料，用Pb(Z=82)作為檢測干擾物，另外選擇生活中常見的Fe(Z=26)作為一般物。將U、Pb、Fe分別設置為上述環境中的待測材料，分別仿真獲取500 000個繆子散射數據，作為不同材料的繆子庫倫散射角數據集；同時在不放置待測材料時，獲取相同數量的繆子散射數據作為空氣的庫倫散射角數據集。最終仿真得到各材料繆子庫倫散射角數據集，其統計均方根列于表1，由于仿真環境的探測空間中為空氣背景，實際在Geant4中測得的繆子庫倫散射角為空氣和待測材料共同作用的結果，不同材料的庫倫散射角數據集均方根略大于式(2)計算的理論值。

3 繆子庫倫散射探測數據特征分析

已知探測空間內的材料情況，則可理論推導繆子散射探測數據的統計量和概率分布函數；反之，通過繆子散射探測數據估計概率分布函數參數或計算其統計量，可反演出探測空間內的材料情況。結合機器學習分類算法，可根據這些數值對不同材料進行分類，將這類數值稱為繆子散射探測數據的分布特征。分別利用繆子散射探測數據估計概率分布函數參數和高階統計量作為分布特征進行分析，使用支持向量機(SVM)測試分布特征對不同材料分類的性能。

表1 Geant4仿真的10 cm厚度材料 繆子庫倫散射角數據均方根Table 1 RMS of Geant4 simulation data for Muon Coulomb scattering angle of 10 cm thick material

3.1 基于概率分布函數參數的分布特征分析

穿過探測空間的繆子可分為兩類：穿過待測物體發生了庫倫散射的繆子和未穿過待測物體在空氣中散射的繆子。因此探測得到的繆子散射數據的概率分布函數為高斯混合模型(GMM)，由待測物體的繆子庫倫散射分布和空氣中繆子庫倫散射分布組合而成。探測空間內僅1種材料時，繆子散射探測數據服從2分量混合的GMM：

P(θ)=(1-ρ)G0,σair(θ)+ρG0,σobj(θ)

(5)

其中：σair為空氣背景中繆子庫倫散射分布均方根；σobj為穿過待測材料的繆子庫倫散射分布均方根；G0，σ(θ)為高斯分布函數；ρ為混合高斯分量的占比，表示了繆子探測數據中穿過待測物體的數據占比，根據GMM的定義，有：

(6)

其中：N為探測到的繆子探測數據總量；Nobj為N個探測數據中來自待測物體的繆子數據量。待測物體體積越小，ρ越低，分辨待測物的難度就越大，如圖2所示，ρ可近似用待測物體在探測器上的投影面積與探測器面積之比表示。在1 m3立方體的探測空間中，探測10 cm邊長的立方體時，ρ=0.01。因此ρ可描述探測空間和待測目標之間的大小關系，以及在此條件下探測獲取的繆子散射探測數據構成，同時也可作為區分不同材料時分辨率的指標。

圖2 探測空間和待測物體大小關系Fig.2 Size relationship of detection space and object to be measured

a——U、Pb、Fe的σEM和ρEM ；b——使用σEM和ρEM分類不同材料的準確率圖3 使用σEM和ρEM分類不同材料Fig.3 Classification for different materials by σEM and ρEM

(7)

其中：θi為繆子探測數據；Qi為隱變量，計算式為：

Q(n)=

(8)

空氣中的繆子散射分布均方根σair已知，最終σobj收斂結果即為EM算法對待測物體的繆子庫倫散射分布均方根估計σEM。ρ也可由EM算法估計的σEM進一步計算隱變量得到：

(9)

由式(7)～(9)，使用EM算法可直接估計出式(5)的所有未知參數得到σEM、ρEM。根據式(6)，代入ρ和N計算出繆子散射探測數據中分別來自空氣和材料的數據量，從Geant4仿真獲取的不同材料庫倫散射數據集抽樣得到繆子仿真探測數據。取0.01≤ρ≤0.4，N=10 000的U、Pb、Fe的繆子仿真探測數據并用EM算法計算(σEM，ρEM)，如圖3a所示。對于標簽已知的分類問題，使用SVM對圖3a中不同材料的分布特征計算分類模型并測試，結果如圖3b所示。(σEM，ρEM)對Fe的分類準確率為100%；在ρ>0.01時，對U和Pb的分類準確率大于95%；而在ρ=0.01時，由圖3a可見，U和Pb的(σEM，ρEM)存在較多重疊，Pb的分類準確率為98.3%時，U的準確率僅60.1%，兩種材料分布特征無法區分。

為進一步分析U和Pb的(σEM，ρEM)出現重疊的原因，用Geant4仿真數據測試EM的估計效果。取0.01≤ρ≤0.4，N=10 000的U的繆子仿真探測數據計算(σEM，ρEM)并與實際值對比。如圖4所示，由于(σEM，ρEM)隨著ρ趨近于0收斂于(σair，0.5)，相當于無材料的情況，使得ρ<0.05時EM估計結果主要受空氣中的繆子庫倫散射數據影響，這與圖3a中U、Pb、Fe的(σEM，ρEM)隨σEM減小趨向同一方向并出現重疊一致。綜上所述，使用(σEM，ρEM)作為分布特征對不同材料進行分類，在繆子散射探測數據中來自待測材料的數據占比大于1.2%(ρ≥0.012)時，對U、Pb、Fe的分類準確率可達到95%以上。但在實際檢測中ρ通常會更小，對10 cm邊長的立方體待測材料(ρ=0.01)，(σEM，ρEM)并不能有效區分U和Pb。

圖4 U繆子探測數據σobj和ρ的EM估計Fig.4 EM estimation of U Muon detection data of σobj and ρ

3.2 基于高階統計量的分布特征分析

上述討論主要使用繆子散射探測數據二階統計量對不同材料進行分類，然而二階統計量局限于對高斯假設的分析，對非高斯數據的估計存在偏差。更高階的統計量相比均值和方差包含有更豐富的非高斯信息[16]，其中峭度是描述數據分布峰的四階累計量，對于零均值分布的數據，峭度計算式[17]為：

(10)

對于ρ不同的繆子散射探測數據，待測材料的繆子散射探測數據分布峰存在不同。而峭度常用于與正態分布比較數據分布峰的大小，在信號處理領域中已有廣泛應用[17]，對于GMM這類非高斯分布，嘗試將高階統計量作為分布特征進行分類。為簡化計算，一般用歸一化的四階矩代替式(10)計算峭度。式(5)描述了繆子散射探測數據分布的GMM，其峭度理論值KGMM可由ρ和σobj表達：

(11)

由式(11)，代入U、Pb、Fe的σobj計算KGMM取最大值時ρ為0.001 4、0.002 7和0.009 1，據此推測峭度在ρ≥0.01時其值仍主要受σobj的影響，可較好描述繆子散射探測數據分布。根據式(11)可得出ρ≥0.01時與繆子散射探測數據理論峭度KGMM的關系(取表1中U的σobj=40.98 mrad，N=10 000)，并將其與同等條件下的繆子仿真探測數據用式(10)計算的峭度值Kdata對比，如圖5所示。Kdata與KGMM基本符合，但Kdata離散程度較大。由于ρ=0時數據分布為高斯分布，此時KGMM=3。但相比于σEM和ρEM，ρ較小時KGMM和Kdata均有出現明顯減小，據此進一步測試ρ≥0.001時使用σEM、ρEM和Kdata作為分布特征分類不同材料的準確率。

圖5 ρ和2分量高斯混合模型峭度的關系Fig.5 Relationship between ρ and kurtosis of 2-component Gaussian mixture model

取0.001≤ρ≤0.4，N=10 000的U、Pb、Fe的繆子仿真探測數據并計算(σEM，ρEM，K)。如圖6所示，不同材料的K在(σEM，ρEM)重疊的部分具有明顯區分。同樣使用SVM對U、Pb、Fe的(σEM，ρEM，K)計算分類模型并測試，結果如圖7所示。在ρ≥0.01時，U、Pb、Fe的分類準確率均在99%以上，在ρ=0.001時，由于Fe作為一般物與U和Pb的繆子庫倫散射能力差異較大，其分類準確率仍有96.2%，U和Pb兩種相近的物質隨著K收斂分類準確率下降到71.5%和71.1%，在ρ≥0.006時，3種材料的分類準確率均在95%以上。相比使用(σEM，ρEM)，增加K作為分布特征，達到相同分類準確率所需的繆子散射探測數據中來自待測材料的數據占比減小了50%。

圖6 不同材料的σEM、ρEM和K分布Fig.6 σEM, ρEM and K distributions of different materials

圖7 使用σEM、ρEM和K分類不同材料的準確率Fig.7 Accuracy of classifying different materials using σEM, ρEM and K

綜上所述，將繆子散射探測數據的概率分布函數參數和統計量作為分布特征可以對U、Pb、Fe 3種材料實現準確分類，選取的分布特征應在待測材料的繆子散射數據較少時仍能較好描述繆子散射探測數據分布的參數或統計量。前述仿真環境中設置了上下各3層的1 m×1 m的繆子探測器、探測空間為1 m3的立方體，以此構成的一個探測器單元，對于更大體積的探測空間，可用陣列的方式擴展繆子探測器的覆蓋范圍，數據處理時可將每個探測器單元的數據單獨處理，如圖8所示(圖中1 GP=0.304 8 m)。由此提出基于分布特征的宇宙射線繆子核材料快速檢測算法，通過Geant4仿真構建當前環境下一個探測器單元內不同材料的繆子庫倫散射數據集，計算各材料繆子散射探測數據的分布特征并使用SVM計算獲得材料分類模型，實現對核材料的快速檢測。該算法可用N=10 000的繆子散射探測數據對100 cm2×10 cm的U、Pb、Fe進行區分，分類準確率在98.9%以上。

圖8 用于20 GP集裝箱的 宇宙射線繆子核材料探測系統Fig.8 Cosmic ray Muon nuclear material detection system used in 20 GP container

4 總結

本文分析了繆子散射探測數據的分布特征，探究了使用繆子散射探測數據的高斯混合分布模型參數和峭度對不同材料的分辨能力。為大量獲取繆子仿真數據，利用Geant4仿真獲取不同材料的仿真繆子庫倫散射數據集，并以此提出一種快速獲取繆子仿真探測數據的方法。基于上述分析和仿真數據，提出一種基于分布特征的宇宙射線繆子核材料快速檢測算法，并對U、Pb、Fe的繆子仿真散射探測數據進行測試。最終用數據量為10 000的繆子散射探測數據實現對核材料(U)和檢測干擾物(Pb)的區分，證明了使用繆子散射探測數據的分布特征作為材料分類依據的有效性。下一步將引入繆子入射能量和能量損失信息，進一步探究分布特征的選取方法，對多材料情況下不同厚度的材料分類進行研究，并逐步展開在現實環境下的測試實驗。