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

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

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

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

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

1 原理

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

(1)

(2)

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

(3)

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

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

2 仿真環(huán)境的建立和數(shù)據(jù)獲取

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

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

(4)

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

圖1 仿真實(shí)驗(yàn)環(huán)境示圖Fig.1 Illustration of simulation experiment environment

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

3 繆子庫倫散射探測數(shù)據(jù)特征分析

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

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

3.1 基于概率分布函數(shù)參數(shù)的分布特征分析

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

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

(5)

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

(6)

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

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

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

(7)

其中：θi為繆子探測數(shù)據(jù)；Qi為隱變量，計(jì)算式為：

Q(n)=

(8)

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

(9)

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

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

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

3.2 基于高階統(tǒng)計(jì)量的分布特征分析

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

(10)

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

(11)

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

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

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

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

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

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

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

4 總結(jié)

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