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

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

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

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

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

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

1 宇宙線繆子與物質相互作用

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

(1)

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

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

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

(2)

(3)

(4)

(5)

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

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

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

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

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

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

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算法中，成像區域被劃分為許多空間體素，如圖4所示。每一個空間體素認為是單一材料，具有唯一確定的散射密度。圖4中，fN為物質的散射密度。

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

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

3.3 成像模擬結果

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

圖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 米級宇宙線繆子成像系統對被屏蔽核材料的成像結果Fig.9 Imaging result of shielded nuclear material by meter-scaled cosmic ray muon imaging system

4 極大似然算法

4.1 成像原理

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

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

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

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

(6)

出射位置偏移Xi服從均值為0、標準差為Di的高斯分布。其中Di滿足下式：

(7)

根據貝葉斯公式：

(8)

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

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

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

(9)

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

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

(10)

則：

(11)

對每次測量結果都可列出上面的待優化表達式。為求解這樣復雜表達式的最優解(極值點)，可利用ROOT提供的數值優化工具類TminuitMinimizer[10]來實現。

4.2 成像結果

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

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

5 宇宙線繆子成像實驗室原型系統模擬

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

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

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

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

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

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

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

6 小結

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

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