高春宇，湯秀章，陳欣南，范澄軍，陳雁南，李雨芃，呂建友

(中國原子能科學研究院，北京 102413)

近年來，核材料的非法轉移、核擴散等問題時刻威脅著國土安全，宇宙射線μ子散射成像作為一種新興的核材料檢測技術逐漸引起世界各國的重視。μ子是一種天然輻射源，沒有引入額外輻照危害，同時對高原子序數(shù)物質(zhì)敏感，穿透能力強，在核材料檢測技術應用中具備天然優(yōu)勢。在集裝箱、貨物材料走私檢測等要求時效性的應用場景，需要在較短時間內(nèi)(幾分鐘)快速準確地檢查出是否藏匿特殊核材料，從而發(fā)展了一種不需要重建圖像而只依據(jù)少量μ子數(shù)據(jù)特征快速實現(xiàn)材料識別的方法。英國核武器研究所(AWE)與布里斯托爾大學合作開發(fā)了基于分箱聚類理論的散射成像材料識別算法，1 min可以區(qū)分10 cm×10 cm×10 cm鈾塊與鉛塊[1]。美國洛斯阿拉莫斯國家實驗室(LANL)與決策公司開發(fā)了多模態(tài)被動式成像算法，可以在1 min內(nèi)檢測出40英尺集裝箱內(nèi)藏匿的核材料[2]。印度薩哈核物理研究所采用了模式識別方法，識別具有相似模式的材料，拒絕不同模式的材料[3]。目前國內(nèi)的研究主要基于聚類算法。清華大學研究了基于支持向量機的聚類分析和分類器，根據(jù)散射密度的平均值和標準偏差進行物質(zhì)識別[4]。中國科學技術大學通過計算待識別目標物與參考物之間散射角的灰色關聯(lián)度快速識別材料[5]。

本文提出一種將卷積神經(jīng)網(wǎng)絡(CNN)應用于μ子成像材料識別的方法，以滿足快速準確檢測出藏匿特殊核材料的實際應用需求。利用宇宙射線μ子穿過待測材料后的散射角統(tǒng)計信息，根據(jù)其與材料原子序數(shù)的正相關性實現(xiàn)材料的分類識別。材料識別模型是主要由卷積層、池化層和全連接層構成的卷積神經(jīng)網(wǎng)絡結構，實驗測量的散射角數(shù)據(jù)或蒙特卡羅模擬Geant4數(shù)據(jù)通過卷積層進行數(shù)據(jù)特征提取，池化層降低數(shù)據(jù)維度，最后在全連接層給出模型分類結果。利用Geant4搭建μ子散射成像實驗模型，并在中國原子能科學研究院的CIAE-MNT宇宙射線μ子成像檢測實驗裝置上進行驗證實驗，分別測試鋁、鐵、鎢3種不同原子序數(shù)的10 cm×10 cm×10 cm樣品，實驗材料識別模型結構上增加卷積層數(shù)以提取更深層次的數(shù)據(jù)特征，進一步引入殘差結構、特征矩陣改進網(wǎng)絡結構來提高模型的準確度。并分別測試模擬模型與實驗模型在1、3、5和10 min不同測量時間下的識別準確度。

1 μ子散射成像原理

宇宙射線μ子穿過待測樣品時，會發(fā)生多次庫侖散射現(xiàn)象，依據(jù)其散射角分布與材料原子序數(shù)呈正比的規(guī)律，實現(xiàn)對材料的檢測識別。根據(jù)莫里埃理論[6-7]，投影到平面上的μ子的散射角θ近似服從期望為0的高斯分布，高斯分布的寬度為：

(1)

輻射長度Lrad由經(jīng)驗公式[8]計算：

(2)

其中：A為相對原子質(zhì)量；ρ為材料密度；Z為材料原子序數(shù)。

定義散射密度λLrad為標準μ子(βc=1，動量為p0)穿過單位厚度物質(zhì)后散射角分布的均方根，則由式(3)可知[9]：

(3)

則有：

(4)

即散射密度λLrad與材料原子序數(shù)呈正相關，可將散射密度作為材料識別的特征量[10]。PoCA算法對μ子的多次庫侖散射做出了簡化，假設為一次散射的結果，對經(jīng)過探測區(qū)域的有效μ子重建入射、出射徑跡，并計算散射角[11]。將探測區(qū)域均勻劃分為若干網(wǎng)格，獲取每個網(wǎng)格內(nèi)所有μ子的散射角。根據(jù)散射密度的定義，得到第i個網(wǎng)格散射密度的估計值λi為：

(5)

2 實驗裝置與仿真模擬

2.1 實驗裝置

樣品測試實驗在中國原子能科學研究院研制的一套基于漂移管探測器的CIAE-MNT宇宙射線μ子成像檢測實驗裝置上進行。該裝置主要由漂移管μ子探測器、信號電子學、數(shù)據(jù)采集系統(tǒng)、成像處理軟件等構成。漂移管探測器陣列共6層，由上千支漂移管互成90°垂直交叉排列構成，上部分3層探測器組用于測量入射μ子的位置信息，下部分3層用于測量μ子經(jīng)過待測樣品后的散射角度和位置，中間100 cm高擺放待測樣品。漂移管是長250 cm、直徑3 cm的鋁管，中間陽極絲引出測量信號，管內(nèi)充Ar/CO2混合氣體，探測面積最大為220 cm×220 cm。本實驗樣品較小，利用100 cm×100 cm的探測面積，將10 cm×10 cm×10 cm的樣品放入探測區(qū)域中心。

μ子穿過待測樣品前后的位置和角度信息通過上下兩組漂移管探測器組分別測量得到，并通過測量μ子穿過漂移管電離產(chǎn)生電子的漂移時間得到其精確位置信息。μ子徑跡重建包括徑跡尋找和徑跡擬合兩步，徑跡尋找是搜索一定時間窗口內(nèi)漂移管陣列記錄的擊中信號的集合，通過一定的規(guī)則進行判斷，并挑選出μ子徑跡信息；徑跡擬合采用基于最小距離平方和擬合二維投影面上的徑跡，使徑跡到各層信號位置的距離最小。根據(jù)重建的μ子入射、出射徑跡，采用PoCA算法計算得到μ子的散射角分布，依據(jù)散射角均方根與材料原子序數(shù)的對應關系得到待測樣品的空間分布。

2.2 仿真模擬

采用Geant4程序進行宇宙射線μ子散射成像的蒙特卡羅模擬[12]，根據(jù)實驗裝置的實際尺寸構建物理模型，模擬高能粒子與探測器和被探測材料的相互作用，跟蹤粒子的輸運過程，顯示粒子徑跡，如圖1所示。模型由上下兩部分共6層探測器陣列構成，上部分和下部分各3組探測器分別用來記錄入射和出射μ子徑跡，每層探測器高度間隔30 cm，探測面積100 cm×100 cm，上下兩部分探測器之間為探測區(qū)域，間距高度100 cm。選擇鋁、鐵、鎢作為低中高原子序數(shù)的代表材料，將10 cm×10 cm×10 cm的鋁、鐵、鎢塊樣品放置在探測區(qū)域中心，與實驗時的擺放位置一致。

圖1 Geant4模擬實驗裝置結構Fig.1 Experimental facility structure simulated by Geant4

CRY軟件庫由美國勞倫斯利弗莫爾國家實驗室(LLNL)研發(fā)，可鏈接到Geant4程序，根據(jù)模擬的需要設置生成不同時間、不同緯度與不同海拔高度上的宇宙射線粒子簇射分布。使用CRY軟件庫引入海平面宇宙線μ子能譜及角分布源項，以更真實地模擬天然μ子源[13]?？紤]到地球海平面處μ子的平均注量率為10 000 m-2·min-1[14]，模擬顯示在100 cm×100 cm的探測區(qū)域上，1 min約有1 600個μ子能穿過6層探測器，用于重建有效徑跡。因此，在本模型中，整個探測區(qū)域上1 600條μ子散射角數(shù)據(jù)等效于1 min的測量時間，依次模擬不同的測量時間，輸出各層探測器記錄的μ子位置信息。

3 神經(jīng)網(wǎng)絡材料識別方法

3.1 用于模擬的卷積神經(jīng)網(wǎng)絡模型

模擬模型定義一個包含兩層卷積層、兩層池化層和1層全連接層的卷積神經(jīng)網(wǎng)絡結構，如圖2所示。網(wǎng)絡的輸入數(shù)據(jù)由PoCA算法對Geant4模擬數(shù)據(jù)進行計算得到，為50×50的散射密度矩陣。

圖2 模擬條件下基于卷積神經(jīng)網(wǎng)絡的材料識別模型Fig.2 Material discrimination model based on convolutional neural network under simulation condition

第1層卷積層定義64個卷積核，第2層卷積層定義32個卷積核，卷積核大小均為3×3，卷積層通過卷積核與輸入數(shù)據(jù)進行卷積計算，以掃描的方式對輸入數(shù)據(jù)進行特征提取。使用非線性激活函數(shù)ReLU函數(shù)，對卷積層輸出結果做非線性映射，定義如下：

f(x)=max(0,x)

(6)

ReLU函數(shù)保留大于等于0的值，其余所有小于0的值直接改寫為0，以修正線性單元。池化層使用3×3的池化窗口掃描特征圖，采用最大值池化函數(shù)選取池化窗口中所有值的最大值，進一步降低數(shù)據(jù)維度，輸入數(shù)據(jù)中微小的變化、冗余將不會改變池化的輸出，使網(wǎng)絡更加具有健壯性。經(jīng)過多層卷積層和池化層的處理后，50×50的散射密度矩陣被抽象成高階特征，由全連接層對提取的特征進行非線性組合以得到預測輸出，給出分類結果。由于線性模型的表達能力不夠，本文加入了非線性因素將卷積層輸出結果做非線性映射，計算得到概率矩陣[P鋁,P鐵,P鎢]，并進一步根據(jù)預測值與真實值計算模型誤差，使用交叉熵損失函數(shù)進行反向傳播，更新網(wǎng)絡權重，優(yōu)化模型[15]。

3.2 實驗條件下的卷積神經(jīng)網(wǎng)絡模型

散射密度是材料識別的判據(jù)，根據(jù)宇宙射線μ子成像檢測實驗裝置測得的μ子位置及角度信息擬合μ子入射、出射徑跡，進一步計算統(tǒng)計信息得到散射密度矩陣。為增強散射數(shù)據(jù)本身的特征，更多地保留原始數(shù)據(jù)的信息，擴大數(shù)據(jù)維度，將100 cm×100 cm探測區(qū)域劃分為100×100的網(wǎng)格，每個網(wǎng)格大小為1 cm×1 cm。對整個探測區(qū)域計算散射密度得到100×100的散射密度矩陣。為了與模擬模型進行定量對比，對于探測區(qū)域，同樣以1 min的測量時間等效于1 600條有效μ子位置信息。分別基于不同測量時間下實驗樣品鐵和鎢立方體材料區(qū)域的散射密度矩陣，計算散射密度大于0的像素數(shù)量的概率密度，并繪制正態(tài)分布曲線，如圖3所示。可看出，測量時間為1 min時，非零數(shù)據(jù)點較少，矩陣提供的有效數(shù)據(jù)有限，這給分類模型的特征提取工作帶來了難度。測量時間為10 min時，μ子的散射角數(shù)據(jù)較多，鐵塊、鎢塊的散射密度矩陣內(nèi)大于0的像素數(shù)量增加，呈現(xiàn)出較好的分布。同時，相對于中等原子序數(shù)材料鐵，重核材料鎢的非零數(shù)據(jù)點更多。

a——鐵；b——鎢圖3 不同測量時間下材料區(qū)域內(nèi)散射密度大于0的像素數(shù)量的概率密度Fig.3 Probability density of number of pixel whose scattering density is greater than zero in material area at different measurement time

神經(jīng)網(wǎng)絡的訓練過程往往需要大量樣本數(shù)據(jù)，而實驗數(shù)據(jù)的總樣本數(shù)量有限，為此進一步改進網(wǎng)絡結構彌補數(shù)據(jù)量少的不足，本文從多角度提取了μ子的散射數(shù)據(jù)特征，包括最大值、最小值、均值、標準差，在傳統(tǒng)結構上加入了殘差結構及特征矩陣，充分利用散射角數(shù)據(jù)。同時設置測試集盡可能大，樣本量是訓練集的25%，以測試模型泛化能力，使得測試結果可信。

實驗模型經(jīng)過優(yōu)化后為包含4層卷積層、4層池化層、1層全連接層和殘差結構的卷積神經(jīng)網(wǎng)絡，其中卷積層從輸入的散射密度矩陣中提取數(shù)據(jù)特征，獲得局部感知信息，池化層減少網(wǎng)絡參數(shù)，控制過擬合現(xiàn)象，相對于模擬模型結構更加復雜。網(wǎng)絡結構如圖4所示，散射密度矩陣輸入到卷積神經(jīng)網(wǎng)絡中，要依次經(jīng)過4層卷積-池化結構，其中第1層卷積層有8個卷積核，卷積核大小為3×3，分別與散射密度矩陣進行卷積運算以提取不同特征，輸出100×100×8的特征圖，并對輸出結果做非線性映射。在池化層使用最大值池化函數(shù)，對特征圖降維到50×50×8。同樣地，在第2、3、4層進行卷積池化運算，最終輸出維度為25×6×6的特征圖。

圖4 應用于實驗環(huán)境的卷積神經(jīng)網(wǎng)絡材料識別模型Fig.4 Material discrimination model based on convolutional neural network applied to experiment

其中，經(jīng)過第2層的卷積-池化操作后，考慮到可能會丟失原始散射密度矩陣X的部分數(shù)據(jù)特征，引入第2條線，加入殘差結構，第2層最終輸出的結果H(X)為：

H(X)=F(X)+X

(7)

其中，F(xiàn)(X)為第2層卷積-池化結構的運算結果。與傳統(tǒng)卷積神經(jīng)網(wǎng)絡的輸出H(X)=F(X)相比，殘差結構將F(X)與原始數(shù)據(jù)X相加，相當于對輸入X加入一個微小變化，使得殘差在網(wǎng)絡的反向傳播過程中梯度信息更易傳播。

為進一步加強數(shù)據(jù)特征，提高網(wǎng)絡的分類準確度，引入第3條線，將提取的原始數(shù)據(jù)的特征矩陣C=[最小值，最大值，均值，標準差]與輸出的特征圖相加后進入全連接層。全連接層對所有得到的特征進行加權組合，計算得到概率矩陣[P鋁,P鐵,P鎢]，其中概率最大的類別為模型分類結果。對測試集數(shù)據(jù)輸入到神經(jīng)網(wǎng)絡后的識別分類過程進行追蹤，以更加清晰地了解模型中每步過程的效果。圖5為測試集中的1個散射密度矩陣從輸入模型到給出識別結果所經(jīng)過的所有操作。可看出，淺層的卷積池化結構更多的是對圖像邊緣提取特征，提取到的信息全面而概括，同時也會有部分關鍵信息提取，隨著層數(shù)的加深，特征圖越來越抽象，相當于對特征圖的提純過程。經(jīng)過4層的卷積-池化后，100×100的散射密度矩陣被抽象成6×6的高階特征，在全連接層輸出概率矩陣[P鋁,P鐵,P鎢]=[0,0.01,0.99]，模型預測該樣本是鋁的概率為0，鐵的概率為0.01，鎢的概率為0.99，取其中的最大值0.99，即該樣本識別為鎢。

圖5 散射密度矩陣輸入到卷積神經(jīng)網(wǎng)絡后的識別分類過程Fig.5 Classification process when scattering density matrix is input into convolutional neural network

當成像區(qū)域存在多個、置于不同位置的待測物塊時，對整個探測區(qū)域的散射密度矩陣提取異常數(shù)據(jù)位置，進而輸入材料識別模型進行判別。

4 結果與分析

定義材料xi被正確識別的準確度taccuracy為：

(8)

其中：xi∈X，X={鋁，鐵，鎢}；n為待測材料的數(shù)量，n=3。定義總體準確度ttotal為：

(9)

利用Geant4對10 cm×10 cm×10 cm的鋁、鐵和鎢3種材料分別進行1、3、5和10 min測量時間的模擬。將每種材料識別的準確度以及總體準確度作為評價模型的指標，得到的結果列于表1。

表1 模擬條件下卷積神經(jīng)網(wǎng)絡模型的材料識別準確度Table 1 Discrimination accuracy by convolutional neural network model under simulation condition

測量時間為1 min時，μ子的散射事例數(shù)據(jù)量有限，散射密度矩陣中非零數(shù)據(jù)比例較大，導致準確度不高。由于鐵的原子序數(shù)介于鋁、鎢之間，相對于二者來說較易混淆，在3種材料中識別準確度最低，僅為82.1%。3 min時，3種材料的識別準確度均可達到90%以上。5 min時，μ子穿過待測物體的散射角數(shù)據(jù)較多，用來訓練網(wǎng)絡的散射密度矩陣中非零網(wǎng)格的比例有所提高，3種材料的準確度均在98%之上，總體準確度達到98.5%。測量時間為10 min時，鋁與鎢塊的識別準確度均已達到100.0%，鐵塊的誤報率也僅為0.1%，總體材料識別準確度達到99.9%。統(tǒng)計10 min時模擬與實驗條件下μ子穿過鎢塊的投影角分布，實驗結果與模擬結果吻合良好(圖6a)。進一步計算不同測量時間下模擬與實驗模型對鎢塊識別的準確度如圖6b所示，由于在實驗中探測器具有一定的位置分辨率，同時網(wǎng)絡的訓練數(shù)據(jù)集較少，導致鎢塊的識別準確度略低于模擬結果。但是隨著測量時間的增加，實驗模型的準確度逐漸接近于模擬模型，差距逐漸縮小，在10 min時兩個模型對鎢塊的準確度都達到了100.0%，由此說明了模擬與實驗平臺是相對可靠的。

a——投影角分布；b——材料識別準確度圖6 模擬與實驗結果的對比Fig.6 Comparison between simulation and experiment results

在實驗裝置上，同樣對10 cm×10 cm×10 cm的鋁塊、鐵塊和鎢塊進行1、3、5和10 min的測量實驗，結果列于表2。當測量時間為3 min時，實驗模型對重核材料鎢的準確度達到了98.1%，由于低、中原子序數(shù)材料的散射密度較小，探測器的位置分辨率對其散射角的測量影響更大，導致鋁、鐵樣品的識別準確度相對較低。測量時間為5 min時，3種材料的準確度均達到90%以上，總體準確度為94.6%。測量時間為10 min時，該模型對鎢的準確度達到100.0%。結果表明，基于卷積神經(jīng)網(wǎng)絡的材料識別方法能在較短的測量時間內(nèi)實現(xiàn)對低、中、高原子序數(shù)材料的檢測，并有效識別重核材料鎢。

表2 實驗條件下卷積神經(jīng)網(wǎng)絡模型的材料識別準確度Table 2 Discrimination accuracy by convolutional neural network model under experiment condition

為研究不同層數(shù)的卷積-池化結構對神經(jīng)網(wǎng)絡準確度的影響，以實驗條件下測量時間10 min獲得的散射密度矩陣作為訓練集，分別建立包含2層、3層、4層以及5層卷積-池化結構的網(wǎng)絡模型，測試結果如圖7所示，可看出，4層、5層卷積-池化結構能更加有效地提取到散射密度矩陣的特征，提升材料識別的準確度。此外，由于5層結構相比于4層結構準確度提升效果并不顯著，且復雜的網(wǎng)絡結構可能會過度學習訓練集數(shù)據(jù)，易出現(xiàn)過擬合現(xiàn)象。因此實驗模型選擇使用4層卷積-池化結構。

圖7 不同數(shù)量卷積結構的卷積神經(jīng)網(wǎng)絡模型的材料識別準確度Fig.7 Material discrimination accuracy of convolutional neural network models with different numbers of convolutional structures

使用機器學習中的隨機森林、經(jīng)典BP網(wǎng)絡分類算法，與本文提出的卷積神經(jīng)網(wǎng)絡算法進行材料識別效能的比較(表3)，結果表明，卷積神經(jīng)網(wǎng)絡材料識別模型在μ子散射的材料識別上具有一定的優(yōu)勢。

表3 不同算法模型材料識別準確度的對比Table 3 Comparison of material discrimination accuracy of different algorithm models

5 結論

根據(jù)宇宙射線μ子穿過待測物體時散射密度與材料原子序數(shù)的相關性，將卷積神經(jīng)網(wǎng)絡的方法應用于μ子散射成像中的材料識別。通過卷積進行特征提取，通過池化降低數(shù)據(jù)維度，不斷迭代調(diào)整權重來優(yōu)化神經(jīng)網(wǎng)絡模型，并分別構建了不同結構的卷積神經(jīng)網(wǎng)絡，分析網(wǎng)絡的測試結果，選擇其中最優(yōu)的網(wǎng)絡結構。本文構建了不同的卷積神經(jīng)網(wǎng)絡材料識別模型以分別應用于模擬數(shù)據(jù)及實驗測量數(shù)據(jù)，當測量時間為1、3、5和10 min時，模擬模型的總體識別準確度分別為86.1%、95.6%、98.5%和99.9%；實驗模型的總體識別準確度分別為68.4%、81.3%、94.6%和97.1%。本文采用的卷積神經(jīng)網(wǎng)絡不失為一種有效的μ子成像材料識別方法，在實驗裝置1 m3的探測區(qū)域上對10 cm×10 cm×10 cm的鎢塊5 min即可實現(xiàn)99.1%的準確度，但由于實際應用需要大量學習訓練，還要考慮到更復雜的探測情況，要克服這些限制，還需要更深入地進一步研究，優(yōu)化模型結構。