宇宙線繆子散射成像模擬與算法研究

智 宇，周 靜，陳 雷，李沛玉，趙明銳，劉雯迪，賈世海，張昀昱，胡守?fù)P，于偉翔，李笑梅

(中國原子能科學(xué)研究院 核數(shù)據(jù)重點實驗室，北京 102413)

宇宙線繆子是由高能宇宙線(絕大部分為高能質(zhì)子、α粒子)與大氣中的分子發(fā)生核反應(yīng)及次級反應(yīng)產(chǎn)生的，其平均能量高達(dá)4 GeV，可穿透2 m左右的鋼鐵及1 m左右的鉛和鈾[1-2]。由于宇宙線繆子是隨機散布在三維空間中，可實現(xiàn)物體的三維成像。此外，根據(jù)宇宙線繆子入射前后傳播方向的變化(即散射角)可反映出繆子所穿過空間體素的材料的原子序數(shù)，由此可實現(xiàn)輕重元素的鑒別。宇宙線繆子散射成像技術(shù)具有高穿透、三維成像及鑒別原子序數(shù)的能力，使得它成為解決核擴散問題的潛在的有力工具。

宇宙線繆子散射成像技術(shù)最早由美國洛斯阿拉莫斯國家實驗室(LANL)提出[3]，利用漂移管探測器測量宇宙線繆子徑跡，其位置分辨率可達(dá)400 μm，首次證明了宇宙線繆子散射成像的物質(zhì)鑒別能力。2009年，美國FIT(Florida Institute of Technology)實驗室首次利用微結(jié)構(gòu)氣體探測器、GEM探測器通過宇宙線繆子散射成像技術(shù)進行了物質(zhì)鑒別和成像實驗，其位置分辨率可達(dá)130 μm[4]。

本文對兩種主要的宇宙線繆子成像算法(PoCA算法和極大似然算法)進行改進，根據(jù)Geant4模擬數(shù)據(jù)驗證米級宇宙線繆子成像系統(tǒng)快速對被屏蔽的重元素鑒別與成像的能力，根據(jù)對10 cm×10 cm的宇宙線繆子成像實驗室原型系統(tǒng)的模擬結(jié)果，對有效事例觸發(fā)率、探測器位置分辨、增益及成像時間等重要實驗參數(shù)進行估算。

1 宇宙線繆子與物質(zhì)相互作用

繆子與物質(zhì)的相互作用主要包括電離、韌致輻射和庫侖散射。其中電離和韌致輻射作用引起宇宙線繆子的能量損失，庫侖散射相互作用引起宇宙線繆子的軌跡偏移。庫侖散射微分截面滿足下式：

(1)

其中：σ為反應(yīng)截面；θ為散射角；Ω為立體角；z為繆子電荷大小，z=1；Z為物質(zhì)原子序數(shù)；e為元電荷量；p和v分別為繆子的動量和速度。

圖1示出多次庫侖散射示意圖。物體沿繆子入射方向的厚度為x，出射位置偏離入射點的距離為yplane，出射方向偏離入射方向的夾角為θplane，入射位置與出射位置連線的夾角為ψplane。對于給定能量的繆子，yplane、θplane和ψplane服從均值為0的高斯分布，它們的標(biāo)準(zhǔn)差分別為：

圖1 多次庫侖散射示意圖Fig.1 Schematic diagram of multiple Coulomb’s scattering

(2)

(3)

(4)

(5)

圖2示出不同材料對應(yīng)的動能為3 GeV的宇宙線繆子的散射密度。由圖2可見，用散射密度可很好地區(qū)分低(Z≤20)、中(20

圖2 散射密度隨材料原子序數(shù)的變化Fig.2 Scattering density vs. atom number

2 宇宙線繆子散射成像模擬

本文基于Geant4[5]開發(fā)了宇宙線繆子散射系統(tǒng)模擬程序。其中宇宙線繆子產(chǎn)生程序用的CRY工具包是由LANL開發(fā)，其功能是產(chǎn)生與宇宙線分布相符的粒子，包括宇宙線繆子、中子、質(zhì)子、電子、光子、π介子等[6]。由CRY生成的宇宙線繆子的能譜和角分布如圖3所示。

為較方便地改變探測器與被測物體的幾何結(jié)構(gòu)、空間位置及材料，采用GDML文件[7-8]來描述宇宙線繆子成像系統(tǒng)的幾何與材料。

宇宙線繆子在穿過探測器和物體的過程中會發(fā)生電離和多次散射，這些物理規(guī)律可用Geant4中定義的物理模型計算。而繆子在探測器靈敏體積內(nèi)引起電離，電子離子對漂移、電子雪崩放大及產(chǎn)生感應(yīng)信號的過程可用Garfield[9]提供的相關(guān)接口進行計算。通過定義G4VFastSimulationModel的子類GarfieldG4-FastSimulationModel可實現(xiàn)Geant4和Garfield物理模型的轉(zhuǎn)換。

a——能譜；b——天頂角分布；c——繆子角度分布圖3 CRY生成的宇宙線繆子能譜和角分布Fig.3 Energy and angular spectra of cosmic ray muon generated by CRY

3 PoCA算法

3.1 成像原理

PoCA算法是由LANL提出的宇宙線繆子散射成像算法[3]。在PoCA算法中，成像區(qū)域被劃分為許多空間體素，如圖4所示。每一個空間體素認(rèn)為是單一材料，具有唯一確定的散射密度。圖4中，fN為物質(zhì)的散射密度。

如圖5所示，對于任意1條宇宙線繆子，可通過上下設(shè)置的探測器得到它的入射和出射軌跡。假定其在成像區(qū)域只發(fā)生1次較大的散射，散射發(fā)生的位置是入射軌跡和出射軌跡的垂直平分線中間點即最近點所在的體素。則該假定折線軌跡所穿過的空間體素中，最近點所在體素貢獻(xiàn)了全部散射密度，將該射線的散射角平方除以射線在該體素內(nèi)穿透距離的值賦予該體素，其余體素沒有貢獻(xiàn)散射密度，將這些體素賦予0。將測量到的每一個宇宙線繆子對應(yīng)的散射密度分布圖求平均即得到整個成像區(qū)域的散射密度分布情況。

圖4 成像區(qū)域體素劃分示意圖Fig.4 Schematic diagram of dividing imaging area into volxel

圖5 PoCA算法示意圖Fig.5 Schematic diagram of PoCA algorithm

圖6 對散射角篩選前(a)、后(b)的最近點分布情況Fig.6 Distribution of the nearest point before (a) and after (b) choosing for scattering angle

3.2 最近點位置分布與成像算法改進

PoCA算法一個重要的假定是用宇宙線繆子的入射軌跡與出射軌跡的最近點所在體素作為散射發(fā)生的位置。為驗證該假定成立的條件，在Geant4模擬程序中對成像區(qū)域為50 cm×50 cm×50 cm，成像物體為分開放置的3個邊長為10 cm的正方體，其材料分別為鐵、鉛和鈾的宇宙線繆子成像系統(tǒng)進行了模擬。在不對宇宙線繆子軌跡進行篩選的條件下，計算得到每一個宇宙線繆子的入射軌跡和出射軌跡所對應(yīng)的最近點的分布情況，如圖6a所示。由圖6a可見，原本應(yīng)集中在成像物體內(nèi)部的最近點，其位置在豎直方向上存在較大噪聲，甚至部分最近點超出了成像區(qū)域范圍。分析模擬數(shù)據(jù)發(fā)現(xiàn)這些在豎直方向上偏差較大的最近點對應(yīng)的宇宙線繆子散射角均較小(基本小于1 mrad)，因此在最終的PoCA算法中需將散射角小于0.8 mrad的宇宙線繆子事例排除。篩選掉小散射角事例后的最近點分布如圖6b所示。由圖6b可見，排除小散射角事例可有效地解決最近點在豎直方向上偏差較大的問題，從而降低由此引起的圖像噪聲。圖7為宇宙線繆子穿過10 cm不同物體后的散射角分布，對鋁這樣的輕元素，散射角小于0.8 mrad的事例占總數(shù)的13.5%，對于鉛這樣的重元素其占比為3.2%。

3.3 成像模擬結(jié)果

為研究米級宇宙線繆子成像系統(tǒng)對被屏蔽物體的成像能力，將成像物體改為被鉛圓柱體包圍的鈾圓柱體，圓柱體底面直徑為10 cm、高為10 cm，探測器尺寸為1 m×1 m，如圖8所示。在模擬程序中共產(chǎn)生了100 000個宇宙線繆子事例，對應(yīng)10 min的照射量，且取探測器的效率為100%。成像區(qū)域被劃分為40×40×40的空間體素(邊長為2.5 cm)。成像結(jié)果如圖9所示，可見當(dāng)成像物體在豎直方向存在互相遮擋的情況時，PoCA算法的成像質(zhì)量顯著下降，具體表現(xiàn)為在豎直方向上兩個物體間的區(qū)域存在較大噪聲。盡管成像質(zhì)量有所下降，但用PoCA算法依舊可很明顯發(fā)現(xiàn)被鉛圓柱體包圍的鈾圓柱體。針對PoCA算法在豎直方向存在屏蔽的條件下成像質(zhì)量不佳的問題，發(fā)展出了基于非線性優(yōu)化的極大似然算法。

圖7 宇宙線繆子穿過10 cm不同材料后的散射角分布 Fig.7 Distribution of cosmic ray muonscattering angle after traversing 10 cm in different materials

圖8 被屏蔽核材料檢測示意圖Fig.8 Testing schematic diagram of shielded nuclear material

圖9 米級宇宙線繆子成像系統(tǒng)對被屏蔽核材料的成像結(jié)果Fig.9 Imaging result of shielded nuclear material by meter-scaled cosmic ray muon imaging system

4 極大似然算法

4.1 成像原理

對于任意一條宇宙線繆子事例，利用探測器可測量得到其入射軌跡和出射軌跡。根據(jù)入射軌跡和出射軌跡可對宇宙線繆子在成像區(qū)域內(nèi)的軌跡進行估算，即以入射軌跡和出射軌跡的最近點所在空間體素為分界，其上方的繆子軌跡與入射軌跡重合，其下方的繆子軌跡與出射軌跡重合，而該空間體素內(nèi)的繆子軌跡為入射軌跡在該體素的入射點與出射軌跡在該體素的出射點的連線。則入射軌跡與出射軌跡間的變化，包括傳播方向、出射位置等，完全由繆子徑跡所經(jīng)過的空間體素的性質(zhì)及徑跡在相應(yīng)體素內(nèi)的長度決定。

圖10 宇宙線繆子軌跡和在空間體素內(nèi)的穿透距離Fig.10 Cosmic ray muon track and transmit length in voxel

對于第i條宇宙線繆子事例，需知道其穿過哪些體素及在相應(yīng)體素內(nèi)的穿透距離。假設(shè)空間體素在豎直方向上將成像區(qū)域劃分N層，則對于任意第j層成像區(qū)域，可首先計算出宇宙線繆子軌跡在該層成像區(qū)域上表面及下表面的交點。由此可得到第i條宇宙線繆子軌跡在第j層成像區(qū)域的水平投影，如圖10所示。然后依次計算繆子軌跡的水平投影與空間體素邊界的每一交點。每產(chǎn)生1個交點則意味著該繆子軌跡在第j層成像區(qū)域內(nèi)多穿透了1個空間體素。記錄每一被穿透空間體素的編號n及相應(yīng)空間體素區(qū)域內(nèi)繆子軌跡的水平投影的長度Pijn。則實際的繆子軌跡在相應(yīng)的空間體素內(nèi)的穿透距離Lijn=Pijn/cosφ，φ為宇宙線軌跡與豎直方向的傾斜角。

對于第i條宇宙線繆子事例，其散射角θi服從均值為0、標(biāo)準(zhǔn)差為Si的高斯分布。其中Si滿足下式：

(6)

出射位置偏移Xi服從均值為0、標(biāo)準(zhǔn)差為Di的高斯分布。其中Di滿足下式：

(7)

根據(jù)貝葉斯公式：

(8)

其中：λijn為第i條繆子射線在第j層成像區(qū)域中穿過的第n個體素的散射密度；P(λ|θ，X)為給定測量結(jié)果散射角θ及位置偏移X，散射密度為λ的概率；P(θ，X|λ)為給定散射密度λ，測量結(jié)果為散射角θ及位置偏移X的概率；P(λ)為散射密度λ的概率；P(θ，X)為測量結(jié)果為θ、X的概率。

成像的目的是找到給定測量結(jié)果對應(yīng)的最可能的散射密度分布，即最大化P(λ|θ，X)。當(dāng)默認(rèn)對被成像區(qū)域內(nèi)部一無所知時，可認(rèn)為P(λ)為均勻分布，即任何散射密度出現(xiàn)的概率相等，則最大化P(λ|θ，X)等價于最大化P(θ，X|λ)，即：

λopt=argmaxP(λ|θ,X)=argmaxP(θ,X|λ)

(9)

其中，λopt為成像區(qū)域散射密度的最可幾分布。

對P(θ,X|λ)取對數(shù)：

(10)

則：

(11)

對每次測量結(jié)果都可列出上面的待優(yōu)化表達(dá)式。為求解這樣復(fù)雜表達(dá)式的最優(yōu)解(極值點)，可利用ROOT提供的數(shù)值優(yōu)化工具類TminuitMinimizer[10]來實現(xiàn)。

4.2 成像結(jié)果

在模擬程序中將成像物體設(shè)置為3個在豎直方向上重合的立方體(邊長為10 cm)，材料分別為銅(上)、鉛(中)和鈾(下)。成像區(qū)域為1 m×1 m×1 m。宇宙線繆子事例總數(shù)為100 000，對應(yīng)的照射時間約為10 min。PoCA算法和極大似然算法(迭代30次)得到的圖像結(jié)果如圖11所示。由圖11可見，極大似然算法可有效消除PoCA算法中存在的在豎直方向上兩物體間產(chǎn)生的噪聲。但在成像時間上PoCA算法只需不到1 min，而極大似然算法約需30 min。

圖11 PoCA算法(a)和極大似然算法(b)成像結(jié)果對比Fig.11 Comparison of imaging result between PoCA algorithm (a) and maximum likelihood algorithm (b)

5 宇宙線繆子成像實驗室原型系統(tǒng)模擬

圖12示出宇宙線繆子成像實驗室原型系統(tǒng)示意圖。在宇宙線繆子成像實驗室原型系統(tǒng)中，Micromegas探測器的尺寸為10 cm×10 cm，成像區(qū)域為10 cm×10 cm×5 cm。成像物體為一長4 cm、寬2 cm的M形鉛塊。產(chǎn)生的宇宙線繆子事例總數(shù)為30 000，對應(yīng)約5 h的測量時間。在4個探測器內(nèi)全部沉積能量的事例有2 934個，占總數(shù)的9.78%，對應(yīng)的平均事例觸發(fā)率為0.163 s-1。

圖12 宇宙線繆子成像實驗室原型系統(tǒng)示意圖Fig.12 Schematic diagram of cosmic ray muon imaging system prototype in laboratory

根據(jù)實驗測得的探測器各通道電容分布和噪聲電容變化曲線[11]，128讀出通道、放大區(qū)間距為100 μm的Micromegas探測器的每個讀出通道的電容平均為46.84 pF，總電容為5 995.52 pF。其對應(yīng)的數(shù)字獲取系統(tǒng)等效噪聲電子數(shù)為33 734.27。根據(jù)Garfield模擬結(jié)果，70%氬氣與30%二氧化碳混合氣體的平均電離能為27.04 eV，取放大區(qū)電子增益為1 000，則為了產(chǎn)生超過3倍噪聲標(biāo)準(zhǔn)差的信號需至少沉積2.73 keV的能量。圖13示出宇宙線繆子事例在各探測器內(nèi)沉積能量分布。由圖13可見，只有24%的觸發(fā)事例滿足這一沉積能量要求。這意味著在實際實驗中需至少20 h的照射時間才能達(dá)到模擬程序中5 h的照射時間對應(yīng)的宇宙線繆子事例總數(shù)。

在不考慮探測器位置分辨率的條件下得到的成像結(jié)果如圖14所示。由圖14可見，雖然圖像存在較大的散粒噪聲，但依舊可辨認(rèn)出由重元素體素組成的M形狀。圖15示出位置分辨率分別為300 μm和1 mm時的成像結(jié)果。由圖15可見，當(dāng)位置分辨率達(dá)到300 μm時，在xy平面上依舊可分辨出M形狀，而在xz和yz平面上噪聲情況較為嚴(yán)重，已不能分辨出M形物體側(cè)視圖形狀。而當(dāng)位置分辨率達(dá)到1 mm時，則在任何方向均無法分辨出成像物體。

a——探測器1；b——探測器2；c——探測器3；d——探測器4圖13 宇宙線繆子事例在各探測器內(nèi)沉積能量分布Fig.13 Energy deposition distribution of cosmic ray muon event in each detector

圖14 不考慮探測器位置分辨率的成像結(jié)果Fig.14 Imaging result without considering position resolution

位置分辨率：a——300 μm；b——1 mm圖15 考慮位置分辨率后的成像結(jié)果Fig.15 Imaging result considering position resolution

6 小結(jié)

本工作開發(fā)了用于宇宙線繆子散射成像的模擬程序。該程序?qū)崿F(xiàn)了基于CRY的宇宙線繆子產(chǎn)生，可較方便地對探測器及成像物體的形狀、位置及材料進行修改，實現(xiàn)了Geant4-Garfield的程序接口，允許在輸出文件中保存Garfield的模擬數(shù)據(jù)，定制了較完整的模擬數(shù)據(jù)結(jié)構(gòu)。在模擬程序的基礎(chǔ)上實現(xiàn)了兩種宇宙線繆子散射成像算法，其中PoCA算法耗時短(1～2 min)，但當(dāng)被成像物體在豎直方向上互相遮擋時存在較大噪聲；極大似然算法耗時較長(30 min)，但對于豎直方向上互相遮擋的物體成像質(zhì)量更好。

根據(jù)模擬和成像研究結(jié)果，對于米級宇宙線繆子成像系統(tǒng)，預(yù)計10 min的照射時間即可獲得較為清晰的成像結(jié)果，且能清楚地發(fā)現(xiàn)被鉛屏蔽的鈾材料。對于10 cm×10 cm成像系統(tǒng)，繆子事例觸發(fā)率為0.16 s-1，要想獲得較為清晰的成像結(jié)果，要求探測器位置分辨率達(dá)到300 μm，探測器增益為1 000時實際測量時間至少需要20 h。本文研究結(jié)果將為今后宇宙線繆子成像實驗研究提供重要的參數(shù)依據(jù)。