缪子对小尺寸中低原子序数物质三维成像技术的模拟研究

季选韬，罗思远，朱 坤，彭肖宇，祝 锦，肖 敏，罗凤娇，王晓冬

(南华大学 核科学技术学院，湖南 衡阳 421001)

宇宙射线缪子是初级宇宙射线与大气层通过级联簇射作用产生的高能带电粒子[1]，是大自然的本底射线源，质量约为电子的207倍，到达海平面的平均能量为3～4 GeV(最高能量可达TeV)，服从与cos2θ呈正比的分布[2]，通量约1 cm-2·min-1。因为缪子穿透性强，所以缪子成像技术被广泛应用于各类交叉学科研究中。当前宇宙射线缪子成像技术主要可分为两大类。第1类是应用在考古、地质勘察领域的缪子透射成像技术，1970年，Alvarez利用缪子穿过物质前后通量的改变，发现了金字塔内的暗室[3]。随后该技术被广泛应用于活火山监测[4-6]、煤矿探索[7]、冰川测量[8]、月球阴影调查[9]、胡夫金字塔内未知的密室成像[10]。2020年，国内研究团队通过测量常熟地下隧道内外的缪子通量分布，首次成功重建了隧道的图像[11]。第2类是应用于核材料检测领域的缪子散射成像技术，2003年美国洛斯阿拉莫斯国家实验室(LANL)首次提出了基于缪子穿过物质前后的散射角信息和POCA(point of closest approach)算法，在短时间内实现了对高原子序数Z物质的成像[12]。随后该技术应用在缪子对铀、钚等重核材料的监测[13-14]和集装箱检测[15]等核安保方面的相关研究。为实现好的成像效果，更多的先进探测器被研制[16-19]。福岛核电站事故后，美国研究团队利用散射成像的原理重建了反应堆内散落的核物质图像[20]，成为散射成像技术的一次重要应用。然而，当前的缪子透射成像和散射成像技术均仅使用了缪子本身的数量和其径迹信息，且主要应用于高Z或大尺寸的物质。实际上，缪致次级粒子也携带了丰富的待测物体信息。Mrdja等提出了基于缪子与缪致次级粒子符合探测技术的缪子成像技术[21-23]，在实验上对牛股骨重建了三维图像，首次实现了缪子对中低Z物质的成像，拓宽了应用领域。

为进一步研究基于缪子与缪致次级粒子符合探测技术的缪子对中低Z物质的成像技术，本文模拟缪致次级粒子的产生和输运行为，研究基于缪子与缪致次级粒子符合探测技术的符合缪子径迹筛选，并结合有限角度成像算法ASD-POCS(adaptive steepest descent projection on convex sets)，实现缪子对中低Z物质的三维成像。

1 缪致次级粒子的研究

1.1 Geant4模拟参数

1) 探测器系统的几何建模

本文利用Geant4软件建立的成像系统的几何模型如图1所示。整体由两部分组成，第1部分是缪子径迹探测器模块，其作用是获取入射缪子的径迹。该模块由系统上方3层位置灵敏探测器构成，如图1中灰色部分所示，相邻两个位置灵敏探测器距离5 cm，单个探测器是尺寸为40 cm×40 cm×1 cm的理想探测器(探测效率100%)，均设置了100 μm的位置分辨率。第2部分是次级粒子探测器模块，用于获取次级粒子的各类信息。该模块由4个闪烁体探测器按照四面环绕的方式排布。单个探测器的尺寸是50 cm×50 cm×5 cm，并依据EJ200型号将其定义为C9H10的塑料闪烁体探测器，如图1黄色部分所示。白色立方体是待测物体，被放置在中央，其尺寸为3 cm×3 cm×3 cm。在模拟中，当缪子进入系统并与待测物体发生作用时，依据作用截面产生缪致次级粒子并从待测物体中射出，入射缪子的位置信息和缪致次级粒子信息分别被两个模块在1个符合时间窗内被记录。

图1 探测器成像系统的几何模型示意图Fig.1 Schematic diagram of geometric model of imaging system

2) 入射粒子及物理过程

Geant4模拟中，粒子源是CRY天然缪子源[24]。模拟时物理过程被限定为G4Decay-Physics、G4RadioactiveDecayPhysics、G4Em-standardPhysics、G4OpticalPhysics 4种模块。

1.2 缪致次级粒子产生过程

缪子作为高能的带电粒子，主要通过4种相互作用损失能量：缪子的电离、缪子的轫致辐射、缪子的电子对效应和缪子的光核反应。当缪子穿过铁时，4种相互作用的能量损失随缪子能量的变化曲线如图2所示。从图2可知，对于平均能量为3～4 GeV的天然缪子源，主要是通过电离作用损失能量。由于缪子本身的能量非常高，通过电离作用产生的电子大部分是能量较高的δ电子。因此缪致次级粒子主要包括两种粒子：δ电子经过电磁级联簇射作用产生的次级电子和次级γ。

图2 不同电磁相互作用下缪子的能量损失随缪子能量的变化曲线[25]Fig.2 Variation of energy loss with muon energy for different electromagnetic interactions[25]

1.3 缪致次级粒子数量随缪子穿过的待测物体厚度的变化趋势

图3示出了不同材料中原初电离电子、次级γ和次级电子的数量与缪子穿过待测物体厚度的关系。其中，原初电离电子是通过Geant4软件限定了物理过程后获取的，是一种理想情况下缪致原初电离电子数据，次级电子指缪子电离产生的电子，包括从待测物体出射的原初电离电子和原初电离电子继续电离产生的电子。为获得全部的次级粒子，模拟中利用离开待测物体的粒子数据取代了闪烁体探测器获取的，但在图像重建时利用的仍是闪烁体探测器获取的数据。图3a、b、c分别展示了在水、铁、铅 3种材料中原初电离电子数量、次级粒子总数、次级γ数量、次级电子数量随缪子穿过的待测物体厚度的变化。从图中可知，在缪子穿过的厚度不超过10 cm时，原初电离电子、次级电子、次级γ的数量与缪子穿过的厚度近似呈线性关系。这是因为，原初电离电子的产额正比于缪子的沉积能量。缪子能量与物质相互作用主要由电离损失能量，因此缪子的能量损失率用-dE/dx表示，其中dE为缪子的微分能量损失，dx为缪子的单位路径。因为在dx路径上缪子的能量变化非常小，即dE几乎保持不变，缪子的能量损失仅与穿过的总厚度相关，所以缪致的原初电离电子随缪子穿过的厚度的增加而线性增加。

a——水中原初电离电子、次级γ和次级电子的数量与缪子穿过待测物体厚度的关系；b——铁中原初电离电子、次级γ和次级电子的数量与缪子穿过待测物体厚度的关系；c——铅中原初电离电子、次级γ和次级电子的数量与缪子穿过待测物体厚度的关系；d——不同材料中原初电离电子数量对比；e——不同材料中次级γ数量对比；f——不同材料中次级电子数量对比图3 不同材料中原初电离电子、次级γ和次级电子的数量与缪子穿过待测物体厚度的关系 Fig.3 Variation of muonic primary ionized electrons, secondary gamma and secondary electrons counts with thickness of muon passes through in different materials

次级电子数量远小于原初电离电子的数量，这是因为大部分的低能原初电离电子因为自吸收效应而被待测物体吸收。从待测物体出射的次级电子可看作由两部分构成，第1部分是产生于待测物体表层(小于0.4 mm)的次级电子，对于水，这部分占比小于25%；第2部分是产生于待测物体内部的，来源于高能的原初电离电子。由于产生于表层的次级电子数量与待测物体的形状、大小无关，是大致不变的值。而高能的原初电离电子同样与待测物体厚度近似呈线性关系，所以次级电子也呈现类似线性关系。对于次级γ，其主要来源于原初电离电子的电磁级联簇射中的轫致辐射。轫致辐射的能量损失率正比于原初电离电子能量，即正比于缪子的沉积能量，因为缪子的沉积能量仅与穿过的总厚度相关，从而次级γ数量也与缪子穿过的厚度近似存在线性关系。次级γ也有可能来源于缪子直接轫致辐射产生，由于轫致辐射的截面与粒子的质量m2呈反比，而缪子的质量约为电子的207倍，所以，缪子直接轫致辐射产生的次级γ的数量相较于原初电离电子的可忽略不计。考虑到实验上难仅通过闪烁体探测器区分次级电子与次级γ，计算两种次级粒子的总和与实际的实验条件更相符。

图3d、e、f对比了不同材料产生的原初电离电子、次级γ和次级电子数量。对于原初电离电子，由于缪子的电离能损与待测物体的Z呈正比，所以水产生原初电离电子数量少于铁和铅；对于次级γ，因为原初电离电子的电磁级联中的轫致辐射截面与Z2呈正比，从而水产生的数量远小于铁和铅；对于次级电子，水的线性程度较好，这是因为次级电子主要产生于水的内部，来源于高能的原初电离电子。铁和铅的线性程度较差，其次级电子主要产生于待测物体表层。

1.4 缪致次级粒子能谱

图4示出了水、铁、铅的次级粒子能谱。图4a为次级γ能谱，从图可知，随Z的增加，能谱的峰值逐渐右移，这是因为缪子的电离损失截面和轫致辐射的截面均与待测物体的Z呈正比。图4b是次级电子能谱，从图可知，铁和铅产生的缪致次级电子能谱趋势变化不大，仅是计数的变化，这是因为铁和铅产生的次级电子主要来源于待测物体表层。水的次级电子能谱的峰值低一些，这是因为水的次级电子主要来源于内部产生的高能原初电离电子，经过自吸收后其能量减小。3种材料的缪致次级电子能谱的峰值均在0.1 MeV以上，在实验上可与本底环境噪声区分开。水产生的次级γ的能量较低，在实验上难以从环境噪声中区分，但由于水的缪致次级粒子是以次级电子为主，所以水仍可通过次级粒子的信息成像。

2 成像原理

2.1 符合缪子径迹

图5a示出了基于符合缪子径迹对物体进行成像的原理图。在1个很短的时间窗内，若缪子径迹探测模块与次级粒子探测模块同时响应，则该入射的缪子径迹被定义为符合缪子径迹。这个很短的时间窗取决于缪致次级粒子从产生到打到闪烁体探测器上的时间。为更好探测次级粒子，本研究使用了具有快时间响应的EJ200塑料闪烁体探测器，其上升时间0.9 ns，衰减时间2.1 ns，脉冲宽度FWHM=2.5 ns。但实验上会存在两种情况导致无法获取全部的次级粒子，第1种是如果1个入射缪子产生了多个次级粒子，且打在了同一片闪烁体探测器上，那么该片闪烁体探测器的时间分辨难以区分多个次级粒子；第2种是若多个次级粒子打在不同的探测器上，且都产生了可观测的信号，由于是在1个时间窗内，所以在实验上也仅被认为是1个信号。因此在筛选符合缪子径迹的过程中，如果1个入射缪子产生了多个次级粒子，符合缪子径迹数量也仅增加1个。在整个过程中，假定入射缪子径迹不发生改变，即符合缪子径迹是一直线。

a——缪致次级γ能谱；b——缪致次级电子能谱图4 缪致次级粒子能谱Fig.4 Energy spectrum of muonic secondary particles

a——符合缪子径迹定义的示意图；b——沿某一方向下符合缪子径迹的数量随沿该方向下待测物体厚度的变化图5 符合缪子径迹的定义与特性Fig.5 Definition and characteristic of coincide muon trajectory

结合图3的结论，在缪子穿过的厚度不超过10 cm的范围内，次级γ和次级电子的数量之和随缪子穿过的厚度呈线性增长，可推出，当沿某一方向下的待测物体厚度增加时，沿该方向入射的缪子则更易于产生缪致次级粒子，从而沿该方向的符合缪子径迹数量更多。利用Geant4软件模拟了沿某一方向下符合缪子径迹的数量与沿该方向下待测物体厚度的关系，结果如图5b所示。水、铁和铅3种材料的某一方向下符合缪子径迹数量均随沿该方向下待测物体厚度的增加而近似线性增加。但符合缪子径迹的线性关系相较于次级粒子数量的差一些，特别是铅的线性关系。这是因为在产生多个次级粒子的符合事件中，符合缪子径迹数量仅增加1个，这使得次级粒子数量的线性关系没有被充分利用，这种情况在次级粒子数量最多的铅中最明显。所以当获取了沿某一方向下符合缪子径迹的数量，则可反推出沿该方向下待测物体的厚度，如果获取的径迹足够多，则可对待测物体重建三维图像，此原理与医学成像的断层成像原理类似。但由于缪子的天顶角服从cos2θ呈正比的分布，即沿水平或接近水平入射方向的天然缪子较少，从而符合缪子径迹被限定在一有限角度的天顶角范围内，因此需结合有限角度的医学成像ASD-POCS算法重建图像。

2.2 ASD-POCS算法

ASD-POCS算法[26]是一种将联合代数重建算法[27](SART)与总变分方法[28](TV)相结合的算法，其主要应用于医学成像中稀疏角度投影问题或有限角度投影问题，与天然缪子角度不完善的情况符合。该算法通过受约束的TV最小化得到离散图像f，如下式所示：

f*=argmin‖f‖TVsubject to

Af=p,f≥0

其中：f为参与迭代的离散待重建图像；A为系统矩阵；p为投影数据。ASD-POCS算法的成像过程可分为以下3个步骤。

第1步，创建投影数据p，和迭代停止条件ε，如下式所示：

投影数据p相当于密度函数μ(x,y,z)沿某射线做线积分，应用到缪子成像领域时，可用沿某方向的符合缪子径迹密度作为投影数据p，而迭代停止条件ε的取值取决于投影数据质量以及配置参数。

第2步，通过SART将所有投影数据进行迭代，不断减少重建图像的估计投影数据与真实投影数据的残差，如下式所示：

当获得180°的投影数据时，可直接通过迭代将图像完整重建出来。但由于探测器的面积有限，而符合缪子径迹相当于一同时穿过缪子径迹探测器和待测物体的直线，所以能探测到的符合缪子径迹会被限定在小于180°的有限角度内，如图6所示，所以仅通过第2步无法获得完整的图像。

图6 在有限角度下收集符合缪子径迹的示意图Fig.6 Schematic of collection of coincide muon trajectories at limited angle range

第3步，采用最速下降算法(SD)使TV最小化，如下式所示：

min TV(f) subject toAf=p,f≥0

完成第3步后，再回到第2步，这样交替进行第2步和第3步，直到满足迭代停止条件，最终获得一离散图像f的最优解。

3 缪致次级粒子的研究

3.1 不同材料的成像效果对比

利用ASD-POCS算法和符合缪子径迹信息重建的三维图像如图7所示，从左到右3个立方体填充的材料为铅、铁、水，体积均为3 cm×3 cm×3 cm。入射的缪子数量为5×107个，换算为天然宇宙射线缪子的通量约为20 d的实际成像时间。由于铁和铅产生的符合缪子径迹数量远大于水，对3种物质同时成像可能导致重建图像中的水接近透明。为更好比较不同材料的成像效果，Geant4模拟时是在相同条件下分别进行模拟再汇于同一三维矩阵中。当待测物体为水时，次级电子对成像的贡献达到90%。当待测物体为铁和铅时，次级γ对成像的贡献占主要，分别为72.5%和82.1%。从图7可知以下结论。

1) 铁的重建图像内部无空洞，噪声点最少，成像效果最好。

2) 水的内部存在较少中空，噪声点较多，这是由两个原因导致，一是因为符合缪子径迹数量较少，从而SART迭代的驱动矩阵的赋值略小，再经过ASD-POCS算法迭代后仍存在大量离散的噪声点，不过该问题可通过增加成像时间来解决；二是因为部分次级电子集中产生于待测物体表层(小于0.4 mm)，导致边缘一层像素的数据量大，从而导致内部略微中空。模拟结果可发现，水的成像效果仍好于铅。

3) 铅的内部中空的部分相较另两种材料较大，噪声点也较多，这是因为在铅中心产生的次级粒子更难逃逸出铅块，所以穿过铅中心部分的符合缪子径迹较少，导致成像的中空和大量的噪声点。因此，该技术适合于对低、中序数物质成像，对中Z物质的成像效果最好，高Z物质的重建图像中间可能会有中空。

3.2 复杂待测物体的成像效果

本文构建了“USC”的几何模型，其高度均为2 cm，字母“U”和“C”的凹槽宽度为2 cm，字母“S”的凹槽宽度为0.8 cm，3个字母全部填充材料为铁，入射缪子数量同样是5×107，重建的图像如图8所示。从图中可知，重建的图像可清晰区分“USC”，尺寸与位置与几何建模匹配，且可区分“S”中0.8 cm的凹槽。但成像的不足也很明显，3个字母均有微小的倾斜，且字母周围的噪声点较多，这是物理模型与算法的数学模型不匹配导致的，算法中假定沿某一方向的符合缪子径迹的数量与沿该方向下的待测物体厚度呈线性关系，但实际的过程中仅是近似呈线性关系，这导致ASD-POCS算法最后的收敛结果噪声点较多，甚至导致了重建图像的轻微倾斜。可通过优化算法的数学模型改善噪声点，实现更好成像精度。

a——铅、铁、水的三维成像的俯视图；b——铅、铁、水的三维成像的三维图图7 不同材料的三维成像结果对比图Fig.7 Comparison chart about 3D imaging result of different materials

a——“USC”重建图像的俯视图；b——“USC”重建图像的三维图图8 几何模型“USC”的三维成像结果Fig.8 3D imaging result of "USC" geometric model

4 小结

本文研究了基于缪子与缪致次级粒子符合探测技术与有限角度成像算法的缪子对小尺寸中低Z的三维成像技术，包括模拟次级粒子的特性、研究基于符合缪子径迹的成像原理、编写ASD-POCS有限角度的三维成像算法。结果表明，次级粒子由次级γ和次级电子组成；其能谱均大部分大于0.1 MeV，可与本底区分；两种次级粒子数量之和与缪子穿过的待测物体厚度近似呈线性关系，所以基于缪子及缪致次级粒子的符合探测技术建立的某一方向下的符合缪子径迹数量也与该方向下的待测物体厚度近似呈线性关系，以此为成像原理再结合ASD-POCS成像算法进行三维成像。重建的三维图像表明，该技术更适合于对小尺寸中、低Z物质成像，对小尺寸中Z物质的成像效果最好，高Z的重建图像内部会有中空。图像中可区分字母“S”里0.8 cm的凹槽，通过进一步的算法优化可实现更好成像精度。