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

(中國原子能科學(xué)研究院,北京 102413)

宇宙射線繆子具有穿透力強、對重核材料敏感的特點,近年來被廣泛應(yīng)用于核材料檢查等領(lǐng)域.繆子與不同原子序數(shù)材料發(fā)生的多重庫侖散射效果不同,利用這點可以對被測物體進行成像以及材料鑒別,而在該過程中引入繆子的能量信息可以提高材料鑒別的準確度.本文基于原子能院研制的繆子成像裝置開展了5 種樣品的材料鑒別實驗,使用離散能量擬合近似連續(xù)能量的繆子散射角分布,進而測量出各材料的輻射長度并以此作為特征量進行材料鑒別.實驗結(jié)果表明,在約1400 個有效繆子事例下,Al-Fe 和Fe-Pb 可以在95%的置信水平下有效區(qū)分,該方法相比于不引入繆子能量信息對Pb-W 鑒別的準確率提高了18.5%.

1 引言

繆子是一種由宇宙射線與大氣層相互作用產(chǎn)生的高能帶電粒子,質(zhì)量是電子的207 倍,平均能量為3—4 GeV[1].繆子的穿透能力強,可以有效地穿過屏蔽層對物體內(nèi)部結(jié)構(gòu)進行成像探測,并且無放射性危害.目前繆子成像技術(shù)已被應(yīng)用于核材料檢查[2,3]、金字塔探測[4,5]、陵墓探測[6]、火山監(jiān)測[7,8],以及對反應(yīng)堆[9,10]和核廢料[11]的監(jiān)測等多個領(lǐng)域.

2003 年美國洛斯·阿拉莫斯國家實驗室(LANL)首次提出利用繆子散射成像方法進行核材料檢查[12],該技術(shù)利用繆子穿過物體時發(fā)生多重庫侖散射的角度偏轉(zhuǎn)與材料原子序數(shù)、繆子能量的規(guī)律對被測物進行成像[13]和識別[14].在實際應(yīng)用中,天然繆子的能量是連續(xù)的且實時測量非常困難,因此大多數(shù)研究使用繆子的平均能量代替未知繆子能量,這種近似導(dǎo)致了成像圖像質(zhì)量和材料鑒別準確度的降低.Vanini 等[15]開展了不同能量精度下的高爐成像模擬研究,Bae 和Chatzidakis [16]開展了多種特殊核材料的材料鑒別模擬研究,Morris 等[17]提出了繆子能量的多群模型等.這些模擬結(jié)果均表明繆子能量信息對于成像和材料鑒別有著重要影響.此外,部分學(xué)者也開展了繆子能量測量的實驗研究.例如,加拿大團隊[18]研制出了帶能量測量結(jié)構(gòu)的CRIPT 探測器,通過測量已知材料阻擋繆子產(chǎn)生的散射角還原繆子動量;清華大學(xué)羅志飛[19]利用多氣隙阻性板室(multi-gap resistive plate chamber,MRPC)時間分辨率高的特性,使用飛行時間法計算出繆子的分段能量信息.

本文基于原子能院研制的繆子成像裝置開展了5 種不同原子序數(shù)樣品的材料鑒別實驗研究,采用離散能量繆子散射角分布擬合近似連續(xù)能量的繆子散射角分布,以輻射長度作為特征量進行材料鑒別,并分析實驗測量誤差對材料鑒別準確度的影響.

2 基于繆子離散能量信息的材料鑒別方法

當繆子穿過靶材料時,會受到原子核庫侖力的作用發(fā)生多重庫侖散射,其散射角在同一平面內(nèi)的投影近似服從均值為0、標準差為σθ的高斯分布[20].根據(jù)Molière 理論[21,22],散射角分布寬度與繆子能量、材料輻射長度有如下關(guān)系:

其中,輻射長度Lrad是材料的特征量,通常原子序數(shù)越大的材料,輻射長度越小.βc是繆子的速度,約等于光速,p是繆子的動量,L為材料的厚度.(1)式表明繆子散射角與繆子能量有關(guān).天然繆子的能譜是連續(xù)的,實驗測量到的繆子散射角是不同能量的繆子產(chǎn)生的散射角集合.單能的繆子散射角服從高斯分布,對天然繆子散射角分布的貢獻為其能量對應(yīng)的高斯概率密度函數(shù),

由此推廣到整個繆子能量區(qū)間,天然繆子散射角滿足的分布應(yīng)為全能量區(qū)間上的繆子產(chǎn)生的散射角貢獻加權(quán)求和,權(quán)重為各微分能量點上出現(xiàn)的繆子概率.該分布理論上為高斯分布與繆子能譜的卷積.由于繆子能譜較為復(fù)雜,本研究采用繆子的9 個特征能量點代替繆子能譜.在天然繆子能量區(qū)間上以對數(shù)線性選取0.25,0.50,1.00,2.00,4.00,8.00,16.00,32.00,64.00 GeV 這9 個能量點,它們的跨度0.25—64.00 GeV 區(qū)間上包含了94.3%的天然繆子,覆蓋了繆子的主要相互作用能量區(qū)間,具有理論意義.用9 個特征能量點上的繆子散射角分布加權(quán)求和近似連續(xù)能量繆子的散射角分布,對實驗測量到的天然繆子散射角統(tǒng)計點進行擬合,可得

擬合函數(shù)(3)式為9 個能量點對應(yīng)的高斯概率密度函數(shù)加權(quán)求和形式,常數(shù)項已歸納到待定系數(shù)中.通過測量標準材料的繆子散射角標定權(quán)重系數(shù)A1—A9,在已知材料種類時各能量點對應(yīng)的σi為已知量,可以由(1)式計算.標定得到的權(quán)重A1—A9代表各自對應(yīng)能量點上的繆子數(shù)量占比.將A1—A9代入(3)式得到一個近似天然繆子散射角分布的表達式.此后測量未知材料的繆子散射角并與標定后的(3)式耦合進而推算出未知材料的輻射長度.該方法從統(tǒng)計角度引入了繆子的離散能量信息.

綜上所述,基于離散能量信息進行材料鑒別分為兩個步驟: 1)利用已知厚度的鉛作為標準材料標定各能量點的權(quán)重系數(shù);2) 測量待鑒別樣品的散射角,通過標定后的繆子散射角分布計算出各材料的輻射長度,并以此對材料進行鑒別.

3 實驗配置與蒙特卡羅模擬

基于繆子離散能量的材料鑒別實驗在中國原子能科學(xué)研究院研制的宇宙射線繆子成像裝置上進行.該裝置的主體部分由六層探測器陣列構(gòu)成,每層由各兩排垂直相交的氣體漂移管緊密排列組成.上、下部分各三層分別用于測量入射和出射繆子的徑跡坐標,中間約1 m 間距的探測區(qū)域擺放實驗樣品.本文采用Geant4 對實驗過程進行模擬,Geant4 是歐洲核子中心開發(fā)的蒙特卡羅軟件[23],用于模擬微觀粒子與宏觀物質(zhì)相互作用的全過程.

3.1 離散能量點權(quán)重的實驗標定

權(quán)重標定實驗選用已知厚度的鉛標準材料測量散射角并計算各離散能量所占權(quán)重.由于鉛的厚度過大會阻擋低能繆子,厚度過小又使高能繆子產(chǎn)生的散射角過小,因此本實驗選用5,10,15 cm 三種厚度的鉛塊,取標定結(jié)果的平均值作為最終權(quán)重.將三種厚度的鉛塊同時放入宇宙射線繆子成像裝置中進行測量,實驗設(shè)置和成像結(jié)果如圖1 所示.

圖1 三種不同厚度鉛塊及散射成像結(jié)果Fig.1.Lead cube with different thickness and the image of scattering tomography.

根據(jù)材料在探測區(qū)域內(nèi)放置位置的先驗知識篩選出穿過完整厚度鉛塊的繆子.選取其中穿透每種厚度的繆子各90000 個,計算上述繆子在XZ和YZ兩正交平面上形成的平面散射角.兩平面內(nèi)的散射角獨立、同分布,可以混合為一個樣本集.在–200—200 mrad 的范圍內(nèi),以1 mrad 為間隔劃分區(qū)間,統(tǒng)計測量散射角的歸一化頻數(shù)后做出散點圖并用(3)式進行擬合.在Geant4 中構(gòu)建與上述實驗相同的鉛塊模型,探測器陣列面積和位置均與實際裝置一致.繆子從最頂層探測器表面1 m ×1 m 的正方形區(qū)域內(nèi)均勻抽樣發(fā)射.入射方向和繆子能量均由CRY 庫抽樣產(chǎn)生[24].模擬計算繆子穿透鉛塊產(chǎn)生的散射角分布,圖2 為本實驗的Geant4 模型示意圖.

圖2 權(quán)重標定實驗的Geant4 模型Fig.2.Geant4 model for weight calibration experiment.

三種不同厚度鉛塊的散射角分布如圖3(a)所示,其中10 cm 組的實驗結(jié)果與模擬結(jié)果對比如圖3(b)所示,實驗結(jié)果已扣除本底影響.每個厚度的實驗組可以擬合出一組權(quán)重A1—A9,將三組權(quán)重取平均值,得到的實驗結(jié)果及模擬結(jié)果見表1,權(quán)重系數(shù)A1—A9代表以對應(yīng)的特征能量點為中心的能量分段中包含的繆子數(shù)量占比.由表1 可知,A3—A5的數(shù)值在實驗和模擬結(jié)果中均較高,表明在天然繆子中1—4 GeV 的繆子數(shù)量較多,對總體散射角分布貢獻較大.由于實驗測量存在一定誤差,部分權(quán)重系數(shù)的實驗與模擬結(jié)果存在差異.

表1 離散能量權(quán)重模擬結(jié)果與實驗結(jié)果Table 1.Discrete energy’s weights of experiment and simulation.

圖3 (a) 繆子穿過不同厚度鉛塊的散射角分布(點)與擬合曲線(線);(b) 10 cm 鉛塊的實驗結(jié)果與模擬結(jié)果對比Fig.3.(a) Measured scattering angle distributions of lead cubes under different thicknesses (dot) and fitting curve (line);(b) the comparison between experiments and simulations of 10 cm Pb.

3.2 引入離散能量的材料鑒別實驗與分析

材料鑒別實驗選用C,Al,Fe,Pb,W 5 種樣品作為待測材料,利用標定后的繆子散射角分布和各材料散射角測量值計算材料輻射長度.5 種材料的厚度均為10 cm,其中Al,Fe,Pb,W 為底面積10 cm × 10 cm 的立方體,C 為底面直徑10 cm 的圓柱體.實驗設(shè)置和成像結(jié)果如圖4 所示.將5 種材料樣品放入繆子成像裝置中,經(jīng)一段時間測量后,穿透每種材料完整厚度的粒子數(shù)各40000 個.在–200—200 mrad 內(nèi),以1 mrad 為間隔劃分區(qū)間,統(tǒng)計各材料散射角的歸一化頻數(shù)并繪制散點圖,如圖5 所示.

圖4 測量五種材料的散射角及散射成像結(jié)果Fig.4.Measurement of scattering angles for the different materials and the image of scattering tomography.

根據(jù)(3)式可知,離散能量點近似的天然繆子散射角分布在θ=0 時取得最大值,將該值與測量到的天然繆子散射角頻率峰值N0(即圖5 中的峰值)對應(yīng),

圖5 繆子穿過不同材料的散射角分布(點)與擬合曲線(線)Fig.5.Measured scattering angle distribution of different materials (dot) and fitting curve (line).

從而計算出材料的輻射長度,

將穿過各材料的繆子事例等分為8 組做平行實驗,按上述方法分別計算5 種材料在平行實驗組的輻射長度值,對平行實驗的結(jié)果取均值以減小隨機誤差.在Geant4 中構(gòu)建C,Al,Fe,Pb,W 5 種材料模型,按照與實驗相同步驟進行模擬計算.表2為5 種材料輻射長度的實驗結(jié)果、模擬結(jié)果以及與勞倫斯伯克利實驗室提供的輻射長度標準值[25]比較得到的相對誤差.由表2 可知,重建5 種材料輻射長度的模擬結(jié)果與輻射長度標準值也存在一些偏差.由于模擬中不存在繆子徑跡的測量誤差,模擬結(jié)果的誤差為重建材料輻射長度算法帶來的方法誤差.該方法誤差的來源為兩部分,一是用離散能量繆子散射角分布近似連續(xù)能量的繆子散射角分布產(chǎn)生的誤差,二是對散射角測量值擬合的誤差.方法誤差通常在15%以內(nèi)浮動.在實驗結(jié)果中,Pb 和W 兩種高原子序數(shù)材料的重建輻射長度誤差較小,分別為4.7%和9.7%;而低原子序數(shù)材料受到實驗測量誤差的影響,重建輻射長度的誤差較大,第4 節(jié)將進一步分析測量誤差對計算不同材料輻射長度值的影響.

表2 不同材料輻射長度的估計值與相對誤差Table 2.Estimated radiation lengths of different materials with relative error.

本文引入受試者工作特征曲線(ROC)評價兩種材料之間鑒別的準確度,該曲線可反映各數(shù)據(jù)集整體的區(qū)分能力.其中AUC 值為ROC 曲線下方面積,AUC 值越接近1,代表各材料之間的鑒別準確度越高.本實驗以材料輻射長度作為特征量,用二分類法對5 種材料中密度相鄰的材料進行兩兩鑒別.將5 種材料分為C-Al,Al-Fe,Fe-Pb,Pb-W四個實驗組,計算每種材料的輻射長度樣本值30 個.待鑒別的兩種材料的輻射長度樣本混合后由小到大進行排序,從最小值到最大值之間等間隔選取閾值對混合樣本逐個判斷鑒別材料種類.對于單個閾值,材料分類準確率為準確分類的樣本個數(shù)與混合樣本總數(shù)之比.隨著閾值的改變,分類準確率也不斷變化.材料鑒別準確率定義為在所有閾值下分類準確率的最大值,對應(yīng)的閾值為最優(yōu)閾值.

繪制各實驗組材料鑒別的ROC 曲線并計算曲線下面積AUC 值.四個實驗組的ROC 曲線及材料鑒別準確率見圖6 和表3.結(jié)果表明,在各樣本值對應(yīng)1400 個繆子事例下,四個實驗組材料鑒別的AUC 值均在0.9 以上.其中Al-Fe 和Fe-Pb 在95%的置信水平下可以區(qū)分;C-Al,Pb-W 在85%的置信水平下可以區(qū)分.

圖6 材料鑒別實驗ROC 曲線Fig.6.ROC curve of material distinguishment experiment.

表3 各材料鑒別準確率Table 3.Distinguishment accuracy of the different materials.

材料鑒別的準確度也與兩種材料間輻射長度的差值大小(材料密度差異)有關(guān).本研究以Pb 材料為參照,模擬了Ag,Cd,Sn,Cu 四種輻射長度遞增的材料與Pb 進行區(qū)分,各材料的鑒別準確率與輻射長度的關(guān)系如表4 所列.使用繆子離散能量的材料鑒別方法可在95%置信水平下區(qū)分輻射長度差值大于0.7 cm 的兩種材料.

表4 材料鑒別準確率與輻射長度的關(guān)系Table 4.Materials distinguishment accuracy versus the radiation length.

對比分析引入繆子離散能量信息與不引入繆子能量信息對材料鑒別準確度的提升效果.在不引入繆子能量信息時,需要假定入射繆子能量為常數(shù),材料鑒別方法只能利用繆子散射角大小作為判斷標準.計算本實驗中繆子穿過Pb 和W 散射角的平方平均值作為特征量,在不引入繆子能量信息下對Pb-W 進行材料鑒別,并與上文中引入繆子離散能量信息的Pb-W 鑒別進行對比,結(jié)果如圖7 所示.引入離散繆子能量信息后,Pb-W 鑒別的AUC值與準確率分別相比前者提升了16.7%和18.5%,由此表明引入離散繆子能量信息可以提高重核材料鑒別的準確度.

圖7 兩種方法鑒別鉛-鎢的ROC 曲線對比Fig.7.ROC curves of the lead-tungsten distinguishment performed by different method.

4 徑跡測量誤差對鑒別準確度的影響

在實驗測量過程中,由于探測器存在一定的位置分辨率,電子學(xué)讀出以及徑跡擬合等過程也存在一定誤差,多種因素影響到測量的繆子坐標與繆子真實位置產(chǎn)生了偏差.為分析探測器徑跡測量誤差對材料鑒別準確度的影響,本節(jié)模擬計算徑跡測量誤差分別為0.5,1.0,1.5,2.0 mm 時5 種材料的輻射長度值.各材料輻射長度相對誤差與徑跡測量誤差之間的關(guān)系如圖8 所示.

圖8 徑跡測量誤差與各材料輻射長度誤差的關(guān)系Fig.8.Relationship between measurement error and the radiation length error of different materials.

圖8 表明各材料的輻射長度誤差隨徑跡測量誤差的增加而增加.在相同水平的測量誤差下,代入繆子離散能量計算低原子序數(shù)材料輻射長度的誤差大于高原子序數(shù)材料輻射長度的誤差.探測器的徑跡測量誤差影響材料輻射長度的計算精度,進而影響材料鑒別的準確度.隨著測量誤差的增大,C-Al 實驗組和Pb-W 實驗組的ROC 曲線變化如圖9 所示.圖9 表明,隨著徑跡測量誤差的增加,兩實驗組的鑒別準確率均不斷降低.與精準測量相比,當繆子徑跡測量存在1 mm 的誤差時,C-Al 鑒別的AUC 值降低了11.1%;Pb-W 鑒別的AUC值降低了9.8%.當徑跡測量存在2 mm 的誤差時,C-Al 鑒別的AUC 值降低了21.6%;Pb-W 鑒別的AUC 值降低了16.6%.由此表明徑跡測量誤差對高原子序數(shù)材料間鑒別的影響較小.

圖9 不同徑跡測量誤差下鑒別ROC 曲線 (a) C-Al;(b) Pb-WFig.9.ROC curve of distinguishment under the different measurement error: (a) C-Al;(b) Pb-W.

5 總結(jié)

本文介紹了基于繆子離散能量進行材料鑒別的實驗方法.基于繆子散射現(xiàn)象,通過已知厚度的鉛塊標定出離散能量的權(quán)重系數(shù),測量樣品的散射角計算得到各材料的輻射長度值,并以此為特征量進行不同材料的鑒別.Pb 和W 兩種材料輻射長度的實驗值與標準值相比,分別相差4.7%和9.7%.在約1400 個有效繆子事例下,分組對比實驗表明,Al-Fe 和Fe-Pb 在95%的置信水平下可以區(qū)分.進一步模擬表明,引入繆子離散能量可在95%置信水平下區(qū)分輻射長度差值大于0.7 cm 的兩種材料,并且相比于不引入繆子能量信息時,對Pb-W鑒別的準確率提高了18.5%.此外,高原子序數(shù)材料的鑒別準確度受測量誤差的影響較小.