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

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

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

1 引言

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)配置與蒙特卡羅模擬

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

3.1 離散能量點(diǎn)權(quán)重的實(shí)驗(yàn)標(biāo)定

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

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

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

圖2 權(quán)重標(biāo)定實(shí)驗(yàn)的Geant4 模型Fig.2.Geant4 model for weight calibration experiment.

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

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

圖3 (a) 繆子穿過不同厚度鉛塊的散射角分布(點(diǎn))與擬合曲線(線);(b) 10 cm 鉛塊的實(shí)驗(yàn)結(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 引入離散能量的材料鑒別實(shí)驗(yàn)與分析

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

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

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

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

從而計(jì)算出材料的輻射長度,

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

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

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

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

圖6 材料鑒別實(shí)驗(yàn)ROC 曲線Fig.6.ROC curve of material distinguishment experiment.

表3 各材料鑒別準(zhǔn)確率Table 3.Distinguishment accuracy of the different materials.

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

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

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

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

4 徑跡測量誤差對鑒別準(zhǔn)確度的影響

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

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

圖8 表明各材料的輻射長度誤差隨徑跡測量誤差的增加而增加.在相同水平的測量誤差下,代入繆子離散能量計(jì)算低原子序數(shù)材料輻射長度的誤差大于高原子序數(shù)材料輻射長度的誤差.探測器的徑跡測量誤差影響材料輻射長度的計(jì)算精度,進(jìn)而影響材料鑒別的準(zhǔn)確度.隨著測量誤差的增大,C-Al 實(shí)驗(yàn)組和Pb-W 實(shí)驗(yàn)組的ROC 曲線變化如圖9 所示.圖9 表明,隨著徑跡測量誤差的增加,兩實(shí)驗(yàn)組的鑒別準(zhǔn)確率均不斷降低.與精準(zhǔn)測量相比,當(dāng)繆子徑跡測量存在1 mm 的誤差時(shí),C-Al 鑒別的AUC 值降低了11.1%;Pb-W 鑒別的AUC值降低了9.8%.當(dāng)徑跡測量存在2 mm 的誤差時(shí),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é)

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