錢(qián)祎劍，張立軍，陳靈新，王冠鷹,*

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

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

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

1 原理

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

(1)

(2)

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

(3)

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

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

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

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

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

(4)

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

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

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

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

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

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

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

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

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

(5)

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

(6)

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

圖2 探測(cè)空間和待測(cè)物體大小關(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為繆子探測(cè)數(shù)據(jù)；Qi為隱變量，計(jì)算式為：

Q(n)=

(8)

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

(9)

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

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

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

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

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

(10)

對(duì)于ρ不同的繆子散射探測(cè)數(shù)據(jù)，待測(cè)材料的繆子散射探測(cè)數(shù)據(jù)分布峰存在不同。而峭度常用于與正態(tài)分布比較數(shù)據(jù)分布峰的大小，在信號(hào)處理領(lǐng)域中已有廣泛應(yīng)用[17]，對(duì)于GMM這類非高斯分布，嘗試將高階統(tǒng)計(jì)量作為分布特征進(jìn)行分類。為簡(jiǎn)化計(jì)算，一般用歸一化的四階矩代替式(10)計(jì)算峭度。式(5)描述了繆子散射探測(cè)數(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ù)此推測(cè)峭度在ρ≥0.01時(shí)其值仍主要受σobj的影響，可較好描述繆子散射探測(cè)數(shù)據(jù)分布。根據(jù)式(11)可得出ρ≥0.01時(shí)與繆子散射探測(cè)數(shù)據(jù)理論峭度KGMM的關(guān)系(取表1中U的σobj=40.98 mrad，N=10 000)，并將其與同等條件下的繆子仿真探測(cè)數(shù)據(jù)用式(10)計(jì)算的峭度值Kdata對(duì)比，如圖5所示。Kdata與KGMM基本符合，但Kdata離散程度較大。由于ρ=0時(shí)數(shù)據(jù)分布為高斯分布，此時(shí)KGMM=3。但相比于σEM和ρEM，ρ較小時(shí)KGMM和Kdata均有出現(xiàn)明顯減小，據(jù)此進(jìn)一步測(cè)試ρ≥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的繆子仿真探測(cè)數(shù)據(jù)并計(jì)算(σEM，ρEM，K)。如圖6所示，不同材料的K在(σEM，ρEM)重疊的部分具有明顯區(qū)分。同樣使用SVM對(duì)U、Pb、Fe的(σEM，ρEM，K)計(jì)算分類模型并測(cè)試，結(jié)果如圖7所示。在ρ≥0.01時(shí)，U、Pb、Fe的分類準(zhǔn)確率均在99%以上，在ρ=0.001時(shí)，由于Fe作為一般物與U和Pb的繆子庫(kù)倫散射能力差異較大，其分類準(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)確率所需的繆子散射探測(cè)數(shù)據(jù)中來(lái)自待測(cè)材料的數(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

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

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

4 總結(jié)

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