多重库透仑射散成射像对精小度尺的度影物响体*缪子

张建鸣 李志伟 刘芳‡ 李景太 冒鑫 程雅苹 庞捷 冯鑫茁 倪四道 欧阳晓平 韩然

1) (华北电力大学核科学与工程学院,北京 102206)

2) (北京卫星环境工程研究所,可靠性与环境工程重点实验室,北京 100094)

3) (中国科学院精密测量科学与技术创新研究院,大地测量与地球动力学国家重点实验室,武汉 430077)

4) (西北核技术研究所,西安 710024)

缪子透射成像方法是一种基于宇宙线缪子穿过目标物体前后的通量变化,进而获得其内部密度结构的无损探测成像方法.缪子透射成像方法假设缪子在低原子序数物质中沿直线运动,但实际上多重库仑散射作用会使缪子路径一定程度上偏离直线,有可能对成像精度造成影响.为此,本文使用Geant4 软件包开展了缪子透射成像蒙特卡罗模拟,针对数米尺度多种密度结构的模型,定量分析了打开和关闭多重库仑散射物理过程时对成像精度的影响.结果表明: 对于数米尺度标准岩石物质,缪子透射成像方法能够很好地恢复内部密度异常几何特征；但多重库仑散射作用对目标物体内部区域近垂向缪子通量造成的偏差可达5%,而在目标物体边界区域的偏差可达13%.因此,需要宇宙线缪子透射成像中考虑多重库仑散射作用的影响,以获得数米尺度目标物体更准确的绝对密度值成像结果.

1 引言

宇宙线缪子来自于初级宇宙射线与大气相互作用产生的次级粒子的衰变,具有衰变时间长、能量高、穿透性强(穿透数百米厚地层)等特点[1].缪子与物质主要发生库仑散射和电离能损两种相互作用,由此发展出两类成像方法: 散射成像和透射成像[2−5].散射成像技术利用缪子多重库仑散射角大小与物质的原子序数相关的性质,适用于区分几厘米至几十厘米、密度相似、原子序数有较大差异的不同材料.2003 年,美国洛斯·阿拉莫斯国家实验室[6]首次实现了缪子多重库仑散射成像,可区分圆柱体钨块和支撑用的两条钢轨支架.随后该实验室先后发展了多种缪子散射的图像重建算法,如最邻近点 (point of closest approach,PCA)、最大似然估计 (maximum likehood scattering and displacement,MLSD)、最大后验估计 (maximum a posteriori,MAP)等算法[7−11],可有效区分钨W、铁Fe、铝Al 等材料.缪子透射成像方法是基于缪子穿过物体时产生的能量损失与穿过路径的密度长度相关的原理,通过分析穿过物体前后的缪子通量分布,即可实现对数米级至数百米尺度物体内部密度结构成像.1955 年,George[12]使用探测器测量了地表和矿井隧道内的缪子通量,根据缪子通量的衰减测定了隧道上方的岩层厚度.20 世纪60 年代末,Alvarez[13]首次将缪子应用到考古领域,测量了金字塔的内部结构.2017 年,Morishima 等[14]基于缪子成像技术发现并验证了金字塔内存在一个至少30 m 长的巨大密室.缪子透射成像技术还被广泛应用于火山监测,隧道探测等领域.中国空间技术研究院于2019 年完成了隧道上覆层结构的缪子成像,并分析了观测时长及数据误差对缪子透射成像的影响[15,16],继而完成了黑龙江五大连池火山的缪子成像观测,获得了火山内部高分辨的密度结构图像[17],显示出缪子成像在地球科学研究中的良好潜力.北京师范大学苏宁等[18]使用宇宙线缪子对秦始皇陵地宫透射成像开展了模拟研究,利用两个视角的投影重建了模型的大小和三维位置,成像结果显示缪子成像可获得皇陵内部的主要结构特征.

鉴于对物体内部密度结构的敏感性,缪子透射成像技术也被应用于小尺度物体内部密度结构成像.单一缪子散射成像在垂直方向易造成图像边界模糊和伪影,而利用缪子吸收信息能限制目标物体的厚度,从而对图像重建加以修正.中国科学技术大学何伟波[19]使用缪子多模态成像方法,同时联合利用缪子的散射和透射信息,可在较短时间内得到质量更好的图像重建结果.而缪子透射模式和散射模式均可完成不同工况下的反应堆堆芯成像[20−22].

针对缪子在介质中直线传播假设的合理性及多重库仑散射作用对于缪子透射成像精度的影响,本文开展了数米级小尺度物体的透射成像模拟,定量分析了多重库仑散射对小尺度物体透射成像精度的影响,可以为小尺度物体缪子透射成像工作提供理论参考.

2 缪子成像方法

带电缪子穿过物体时,缪子受到物体内原子核的库仑电场作用偏转,偏转的方向与入射方向的夹角称为散射角.多次库仑散射平面角分布近似于均值为零的高斯分布,其平面角的均方根 (root mean square,RMS)可表示为

其中,θ为散射角,p,β分别为缪子的动量和相对速度,c为光速,Q为电荷数,X0为物质的辐射长度,

其中,ρ为物体密度,Z为物质的原子序数,A为相对原子质量.在10–3

缪子穿过物体时,同时伴随着电离作用能量损失过程,平均能量损失率可以描述为

a项与电离损失相关,b项为韧致辐射、µ核反应以及电子对产生等辐射能量损失项,−dE为缪子的能量损失,x为质量厚度,是密度分布ρ(ε)在缪子路径L上的线积分:

从(3)式可以看出,缪子的能量损失与穿过物体的密度和路径有关.缪子穿过一定质量厚度的物体后,其剩余能量存在一个最小值Emin,即最小穿透能量,当缪子入射能量Eμ

通过分析缪子穿过物体前后通量的变化,可以计算缪子穿透物质的质量厚度x,再根据先验信息,得到物体的长度或密度分布.

对小尺度物体进行透射成像时,将待测物体放在探测器中间,待测物体上方探测器用于记录入射缪子位置信息,下方探测器用于记录缪子出射位置信息.当探测器间不放置任何模型时,通过上层探测器挑选近似垂直入射的缪子,在下层探测器探测到这些近似竖直入射的缪子的通量分布,记为N(o;i,j);探测器间放置模型时,用同样的方法得到下层探测器的缪子分布,记为N(a;i,j).定义各像素缪子计数的比值 Ratio(i,j)为

Ratio 越小,说明该像素对应竖直方向上的缪子被吸收的越多,该像素对应方向上的密度与厚度之积也越大.因此,当待测物体厚度一定时,Ratio 的分布可反映物体内部的密度结构分布差异.当利用缪子进行小尺度物体的透射成像时,缪子的传播路径示意图如图1 所示: 一部分缪子被吸收,一部分缪子穿透物体,穿透缪子由于散射作用,出射方向与位置可能会发生一定程度的偏移.

图1 缪子穿透物体的路径Fig.1.Path of the muon through the object.

本文首先通过模拟不同能量的缪子穿过不同厚度岩石的出射分布,验证多重库仑散射模拟过程的正确性,并对不同厚度的岩石-空洞模型进行透射成像,对比开关散射物理过程的成像结果,从而定量分析多重库仑散射作用对透射成像精度的影响,完成小尺度物体缪子透射成像和多重库仑散射作用对透射成像精度影响的研究.

3 模拟环境搭建和可靠性验证

本文基于Geant4 (Geometry and Tracking)软件包,针对以标准岩石为主体的多种结构模型进行了宇宙线缪子成像模拟,包括Geant4 散射模块的可靠性验证、单能缪子穿过不同厚度岩石的散射影响分析、开关多重库仑散射过程对透射成像精度的影响分析等.模拟计算使用的宇宙线缪子源由宇宙射线生成程序(Cosmic-ray shower generator,CRY)产生,通过设置不同时间、不同纬度、高度缪子的能量、位置、方向等信息,产生相应分布的宇宙线缪子源[25],模拟产生的107个缪子能量分布见图2.

图2 CRY 随机抽样产生缪子的能量分布Fig.2.Energy distribution of muons generated by random sampling of CRY.

3.1 多重库仑散射模块验证

图3 为模拟计算模型,仿真空间和探测器材料设置为真空,待测物体尺寸为3.0 m×3.0 m×10 cm.在Geant4 程序,设置缪子为单能点源,入射能量为3 GeV,垂直入射到待测物质,缪子数量为106个,物理过程包括缪子散射物理过程.计算统计缪子穿过物体后的散射角分布,并将模拟结果与(1)式计算的散射角RMS 值进行对比.

图3 散射模块验证模型(缪子为单能点源,垂直穿过待测物体)Fig.3.Verification model of scattering module (A single energy muon at a point source,passing vertically through the object to be measured).

对比结果见表1,模拟使用的缪子数目为106个,统计误差为取散射角中心98%分布计算RMS 值,得到材料U 的散射角分布RMS 值为31.7 mrad,与(1)式计算值29.30 mrad 相对偏差约8.20%.同时本文多次改变平板材料并计算得到缪子穿过Cu,Al,U,Fe 等物质后的散射角分布.鉴于(1)式计算值的误差约为11%,模拟值均在经验公式计算值的误差范围内,验证了散射模拟所建立的物理模型及物理过程的可靠性.

表1 3 GeV 的缪子穿过10 cm 厚不同材料的散射角Table 1.Multiple scattering for 3 GeV muons passing through 10 cm of various materials.

3.2 缪子最小穿透能量

对于一定密度长度的物质,穿透该物质需要的缪子最小能量Emin,可通过缪子的连续阻止本领(continuous slowing down ability,CSDA)[26]表格给出,CSDA 结合缪子能损公式和实际观测数据给出了缪子能量和不同材料的平均穿透厚度的关系.本文使用Geant4 模拟不同能量的缪子能穿透的标准岩石的厚度,重新研究了缪子能量和标准岩石的最小穿透厚度的关系.模拟设置单能点源缪子在岩石上表面中心位置垂直入射,缪子能量范围设置为0.4—20.0 GeV,对于每个单能粒子模型入射106次,并记录缪子能量沉积为零时竖直方向的位置坐标.结果显示: 500 MeV 的缪子就可穿透1.5—1.8 m厚的标准岩石(图4(a)),10 GeV 的缪子穿透的标准岩石厚度在15—22 m 间(图4(c)),同能量的缪子所能穿透的厚度是一个区间范围.标准岩石的厚度与最小穿透能量Emin基本呈线性关系(图5).

图4 单能缪子竖直穿透的标准岩石厚度分布 (a) 0.5 GeV;(b) 1.0 GeV;(c) 10.0 GeVFig.4.Standard rock thickness distribution for vertical penetration of muons with different energies: (a) 0.5 GeV;(b) 1.0 GeV;(c) 10.0 GeV.

图5 不同能量缪子能竖直穿透的标准岩石的平均厚度(误差棒给出1σ 误差)Fig.5.Average thickness of standard rocks penetrated vertically by muons of different energies (Error bar is 1σ).

3.3 库仑散射特征分析和缪子透射成像模型

为了分析低能缪子穿透物体后的多重库仑散射影响,实验模拟了不同能量的缪子穿透不同厚度的岩石的出射点分布.计算模型见图6,探测器中间放置100 m×100 m×h的标准岩石,其厚度h为1,3,5,10,30 m.待测物体上方放置探测器1 和2,待测物体下方紧贴着放置探测器3,用于记录缪子的出射位置,并结合探测器4 记录到的缪子位置信息计算缪子的出射角度,探测器1 和2,探测器3 和4 间距均为1 m.缪子在探测器1 中心位置上方2 cm 位置点垂直入射,缪子入射能量根据待测物体的厚度h,选取入射能量大于Emin的缪子.

图6 多重散射特征分析模型 (100 个1.5 GeV 的缪子,垂直穿过3 m 厚度的标准岩石(红色线条),蓝色为缪子径迹)Fig.6.Multiple scattering feature analysis model (100 muons of 1.5 GeV vertically passes through a standard rock with a thickness of 3 m (red line),the blue is the muon track).

为了定量研究多重库仑散射对透射成像的影响,实验设置了使用缪子对不同厚度模型透射成像模拟.模型示意图如图7 所示,其中待测物体尺寸为1.8 m×1.8 m×0.8h的岩石体,内部包括一个尺寸为0.9 m×0.9 m×0.6h的空气腔体,岩石和空气腔体中心重合,h=1,3,5 m.岩石模型放置在探测空间坐标原点,模型上下方各是两层间距为1 m,尺寸为3 m×3 m×2 cm 的探测器,缪子在探测器1上表面入射.为了提高模拟效率,入射缪子的天顶角限定在3°以内.模拟粒子物理过程分别包括开启和关闭库仑散射过程的模拟.

图7 岩石-空腔模型透射成像模拟Fig.7.Rock-cavity model transmission imaging simulation.

4 模拟结果与讨论

4.1 单能缪子穿透岩石的多重库仑散射影响

能量分别为0.7,1.0,3.0 GeV 的缪子垂直穿透1 m 厚的标准岩石的出射位置分布,如图8(a)—(c)所示.从图8(a)可明显看出,700 MeV 的缪子穿透1 m 厚度的标准岩石,多重库仑散射使得部分缪子偏离入射方向,粒子出射位置呈现中间密集,越往边缘越稀疏的分布,多重库仑散射引起的部分缪子路径偏移距离最大可达十几厘米.而3.0 GeV 的缪子穿透1 m 厚度的标准岩石,多重库仑散射引起的部分缪子路径偏移距离在1 cm 内,散射引起的径迹偏移很小.通过不同能量点源缪子垂直穿过不同厚度岩石的路径偏移量(图9(a))和散射角图(图9(b))可以看出,能量较低的缪子散射角和累计路径偏移距离较大,累计路径散射偏移距离与累计散射角随入射缪子能量的增加而迅速减小.如图9(a)所示,3 GeV 的缪子穿透5 m 厚度的岩石引起的路径平均偏移距离为18 cm,占穿过路径长度的3%左右,平均立体散射角接近0.13 rad.总体来说,对于一定厚度的物体,能量接近Emin的缪子(即剩余能量比较小的缪子)的散射作用最为明显.

图8 点源缪子竖直穿过1 m 厚标准岩石的出射位置分布 (a) 0.7 GeV;(b) 1.0 GeV;(c) 3.0 GeVFig.8.Distribution of emission positions of point source muons vertically passing through 1 m thick standard rock: (a) 0.7 GeV;(b) 1.0 GeV;(c) 3.0 GeV.

图9 单能点源缪子穿过不同厚度的标准岩石的多重库仑散射影响 (a) 路径平面平均偏移长度(误差棒给出0.5σ 误差);(b) 累计散射角Fig.9.Multiple Coulomb scattering effects of single-energy point source muons passing through standard rocks with different thicknesses: (a) Average offset length of the path plane (the error bar is 0.5σ);(b) cumulative scattering angle.

4.2 散射作用对小尺度物体透射成像的影响

以不同厚度模型为例,定量分析了多重库仑散射作用对缪子透射成像精度影响.图10 自上而下为岩石厚度为0.8,2.4,4.0 m 的模型,使用的缪子为经过筛选的天顶角小于3°的宇宙线缪子,开启与关闭多重库仑散射物理过程的透射成像结果以及各像素上两种结果的差值.从不同厚度的模型成像结果(图10(a)—(c)和图10(d)—(f))可以看出,不论有无散射过程,岩石与内部空腔的边界能够清晰地看到,缪子透射成像方法可以很好地恢复待测物体密度异常的空间分布和几何形态.但多重库仑散射作用使得穿过高密度区域边缘的缪子更容易散射到低密度区域,增大高密度区域与低密度区域的通量衰减差异(图10(g)—(i)),高密度区域缪子减少得更多,低密度区域缪子减少得更少,同时造成密度异常边界附近出现一些通量增加的假象(图10(a)—(c)).

图10 缪子穿过不同厚度模型前后各像素上计数的比值(俯视图) (a),(d),(g) 0.8 m;(b),(e),(h) 2.4 m;(c),(f),(i) 4.0 m (黑色实线为岩石模型的外边界,黑色虚线为岩石与空腔的分界)Fig.10.The ratio of the counts on each pixel before and after the muons pass through the models with different thicknesses(top view): (a),(d),(g) 0.8 m;(b),(e),(h) 2.4 m;(c),(f),(i) 4.0 m (Solid black line is the outer boundary of the rock model,the dashed black line is the boundary between the rock and the cavity).

图11(a)—(c)和图11(d)—(f)给出不同厚度下A-A′(图10(a)—(c)黄色虚线)和B-B′(图10(a)—(c)红色虚线)切线经过像素上缪子穿过模型前后计数的比值分布,A-A′切线在岩石和外边界附近,B-B′切线横穿了岩石-空腔模型.从图11(a)—(c)可以看出,不开启库仑散射过程,A-A′切线上各像素位置的缪子穿过物体前后的通量比值几乎为1.0,而开启库仑散射过程,A-A′切线上的各像素位置的缪子穿过物体前后的通量比值普遍大于1.0,多重库仑散射作用使得缪子散射到高密度边界区域外.由图11(d)—(f)可以看出,开启散射作用,岩石区域的缪子计数减少比值比不开散射过程的多,而空腔区域的缪子计数减少比值比不开散射过程的整体偏多,散射作用增大了不同密度物质的缪子的通量衰减差异,同时使缪子的出射位置分布更加离散.例如,对于穿过2.4 m 厚度的标准岩石,近垂直方向的缪子因多重库仑散射作用造成的通量衰减占总通量衰减的5%左右,边界附近的占比达到了13%,引起的路径平均偏移量约为6 cm,占穿透路径长度的2%.

图11 不同厚度模型切线经过像素的计数比值分布 (a),(d) 0.8 m;(b),(e) 2.4 m;(c),(f) 4.0 m (A-A′ 线在岩石和外边界附近,B-B′ 切线横穿了岩石-空腔模型 (图10 (a)—(c)))Fig.11.Distribution of count ratios of the tangent lines passing through pixels for different thickness models: (a),(d) 0.8 m;(b),(e) 2.4 m;(c),(f) 4.0 m (A-A′ line is near the rock and the outer boundary,and B-B′ tangent traverses the rock-cavity model(Fig.10 (a)–(c))).

4.3 缪子散射角度挑选后的透射成像

增加探测器可以准确地测出缪子穿透物体前后的方向变化,方向变化较小的缪子直线性更好,有利于提高成像的位置精度.本文使用探测器1 和2 测出缪子的入射天顶角,探测器3 和4 测出缪子的出射天顶角,只使用两者天顶角相差小于1°的缪子用于成像,筛选掉散射角较大的缪子.对于不同厚度模型的岩石-空气模型,成像结果如图12 所示,可以看出经过方向挑选后的缪子成像结果并没有出现岩石外边缘周边通量增加的假象,空气与岩石的边界更加清晰,模型的异常的位置结构更接近于不开启散射作用的图像(图10(d)—(f)),但岩石区域的缪子通量变化更加明显,与不做散射角度筛选的成像结果对比,缪子的通量减少得更多.这是因为在没有放置模型的情况下,缪子穿过的是空气,缪子散射角普遍小于1°,能统计到的缪子几乎没有变化,但在探测区域放置岩石-空腔模型后,一部分缪子穿过模型后的散射角超过了1°,没有被统计到,导致该像素上的N(a;i,j)减小,通量的比值相较于未经过散射角筛选的成像结果减小.

图12 经过缪子散射角筛选后不同厚度模型缪子透射成像结果(俯视图) (a) 0.8 m;(b) 2.4 m;(c) 4.0 mFig.12.Transmission imaging results of the model of different thicknesses after the screening of the muon scattering angle(top view): (a) 0.8 m;(b) 2.4 m;(c) 4.0 m.

5 结论

本文基于Geant4 软件包,开展了多重库仑散射作用对小尺度物体透射成像精度影响的研究.研究结果显示缪子透射成像方法能够很好地恢复待测物体密度异常的空间分布和几何形态,初步验证了使用近垂直缪子对小尺度物体密度结构进行透射成像的可行性,表明缪子透射成像方法在小尺度物体内部密度结构探测中也有很好的成像效果.但是,多重库仑散射作用一定程度上改变了缪子的运动方向和径迹,影响了高密度区域与低密度区域的通量衰减差异,造成在密度差异较大物质边界附近出现通量增加的假象,增大了界面附近的通量衰减差异.通过挑选散射角小于1°的缪子透射成像的方法,可以改善密度差异较大物质边界附近通量增加的假象,但缪子的通量变化更加明显,对于由通量衰减比值恢复密度绝对值的准确性也会产生一定影响.因此,在基于宇宙线缪子进行小尺度物体透射成像时,使用通量衰减比值来恢复物体绝对密度值要注意多重库仑散射引起的误差.