宇宙线缪子散射成像模拟与算法研究

智 宇，周 静，陈 雷，李沛玉，赵明锐，刘雯迪，贾世海，张昀昱，胡守扬，于伟翔，李笑梅

(中国原子能科学研究院 核数据重点实验室，北京 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。本文研究结果将为今后宇宙线缪子成像实验研究提供重要的参数依据。