苏宁 刘圆圆† 王力 程建平

1) (北京师范大学核科学与技术学院,射线束技术教育部重点实验室,北京 100875)

2) (北京师范大学物理学系,北京 100875)

1 引言

宇宙射线是来自于太空的高能射线.对宇宙射线的观测研究[1,2]不仅推动了天体物理、粒子物理等基础研究的发展,同时研究人员提出可以利用宇宙射线中的缪子(µ子)能量高、穿透能力强的特点,通过探测穿透待测物体的宇宙射线µ子的通量、角分布等信息实现对待测物体的成像,这一技术被称为µ子成像技术.µ子成像技术是一种辐射成像技术,根据所利用的µ子与物质相互作用性质的不同,µ子成像技术可分为两类:散射成像和吸收成像.散射成像利用宇宙射线µ子穿透物质时发生的多重库仑散射的散射角大小与材料的原子序数有关的特性,适用于对小型、高密度、高原子序数的物体的成像[3].吸收成像利用宇宙射线µ子穿透物质时的能量损失率与穿透路径上物质的密度、厚度有关的特性来对物体的密度结构成像.由于宇宙射线µ子能量很高、穿透能力强,µ子吸收成像技术可以对尺度达千米量级的物体成像[4].

具体来说,对µ子吸收成像的研究最早可追溯至1955 年,George[5]使用µ子探测器对澳大利亚一处地下工厂设施上方的岩层厚度进行测量.20 世纪60 年代末期,Alvarez 领导的实验组[6]使用µ子探测器对金字塔中可能存在的墓室进行搜寻,这是此技术在考古领域的首次应用.此后近三十年的时间里,对µ子吸收成像技术的研究进展缓慢,主要是对这一技术应用于火山[7]、地下洞穴[8]、矿体探测[9]的模拟分析和初步的实验应用.进入21 世纪后,得益于探测器技术的进步,µ子吸收成像技术得到了快速发展.日本[4]、意大利[10]等国家的科学家利用µ子吸收成像技术对多座火山的内部结构进行了较好的成像,并开展了对火山内部密度结构动态监测的研究[11],尝试结合重力数据进行反投影成像以提高成像精度[12],成像不再局限于二维,而是发展到三维的密度结构成像[13].在对金字塔[14]、地下洞穴[15]、矿体[13]等不同类型观测目标的成像方面也有了一些重要的研究成果.特别是2017 年,Morishima 等科学家[14]分别将不同的µ子探测器先后放置在胡夫金字塔内部和外侧进行探测,发现并验证了大走廊上方存在一个至少30 m 长的未知墓室.这一研究结果表明µ子吸收成像技术在对大型文物的无损探测方面具备极强的应用潜力.而我国作为一个历史悠久的文明古国,现存大量文物遗迹亟待考古研究.针对帝王陵墓等一些大型文物内部结构的无损探测,目前考古学中传统的地球物理方法存在着一定的局限性,例如电法、磁法易受地面环境的干扰;地震波法成本较高[16];探地雷达法的探测深度范围在几米到十几米之间浮动[17],可探测深度相对较浅.且每一种方法仅对特定的一两种物性敏感,存在多解性的问题[18].将µ子吸收成像技术应用于考古领域,可以成为对传统地球物理方法的重要补充.

目前国内对µ子吸收成像的研究较少,且集中在地质探测方面,在µ子吸收成像应用于大型文物成像方面的研究几乎空白.本文以秦始皇陵地宫为研究对象,利用蒙特卡罗模拟的方法研究µ子吸收成像应用于帝王陵墓无损探测的可行性.

2 µ子吸收成像算法

2.1 µ子与物质相互作用原理

宇宙射线µ子是一种轻子,它主要由原初宇宙射线与大气层中的原子相互作用产生的π 介子和K 介子的衰变产生,能量很高,在相对论效应的作用下能够在发生衰变之前到达海平面.海平面观测到的宇宙射线µ子的平均能量在3—4 GeV,平均通量约为1 cm—2·min—1,是到达海平面的宇宙射线的主要成分[19].到达地面的µ子的天顶角θ可以是0°—90°间的任意值,µ子通量随θ的增大而减小.图1(a)为实验中测得的不同天顶角方向入射的海平面宇宙射线µ子动量分布,图1(b)为天顶角θ、方位角φ示意图.

图1 (a) 实验测量得到的不同方向上的海平面µ子通量[20];(b)探测器探测到的µ子的方位角φ 和天顶角θ,其中xOy为水平面Fig.1.(a) Sea-level muon flux at different zenith angles measured in experiment[20];(b) zenith angle θ and azimuth angle φ of the muon detected by a detector.The xOy plane represents for horizontal plane.

µ子穿过物质会通过电离、轫致辐射、电子对产生、核相互作用等方式损失能量,其平均能量损失率可一般地表达为[4]

其中E代表µ子能量,为µ子穿过物质的密度对穿透长度的积分,a对应电离导致的能量损失率,bE对应其他作用导致的能量损失率.a,b与µ子能量和所穿透物质的种类有关.

2.2 图像重建算法

由(1)式可知,µ子穿透一定量的物质后平均损失的能量 ΔE可表示为

其中X0为µ子穿透物质的密度对穿透路径的积分.X0越大,则µ子穿透物质平均损失的能量越多,能穿透物质的µ子越少,故可根据探测器探测到的穿透待测物体的µ子数量的多少推断µ子穿透物质的X的大小.辅助一些已知信息,即可计算未知量.如已知穿透路径上的平均密度时,可计算出路径长度;已知穿透路径长度时,可计算出穿透路径上的平均密度.

本文基于对待测物体掌握的先验知识(待测物体几何尺寸、部分区域的密度分布),假设待测物体内部存在密度分布不同于先验知识的感兴趣区域(region of interest,ROI),所使用的µ子吸收成像的图像重建算法可概括为以下3 步:

第一步,获取µ子通量.假设Nµ(θ,φ) 表示探测器在 (θ,φ) 方向上探测到的µ子数,则探测器探测到的µ子通量Φ0(θ,φ) 为

其中S为探测器面积,Ω为立体角,t为测量时间.

第二步,与先验知识对比,获取ROI 边界的二维角坐标信息.利用已知的待测物体的几何尺寸和密度信息建立参考模型,模拟得到宇宙射线µ子穿透参考模型后剩余的µ子通量Φs(θ,φ),定义µ子通量差f(θ,φ)=Φ0(θ,φ)-Φs(θ,φ) .当f(θ,φ)/=0 时,代表 (θ,φ) 方向上的密度分布与利用已知信息构建的模型的密度分布不同,(θ,φ) 位于ROI 内.若f(θ,φ)＜0,则代表 (θ,φ) 方向上存在密度偏高的结构(如高密度矿物),若f(θ,φ)＞0 ,则代表(θ,φ)方向上存在密度偏低的结构(如空穴).可根据f(θ,φ) 发生突变处对应的 (θ,φ) 确定ROI 边界的角坐标.

第三步,重构ROI 的三维图像.探测器在单个测量点得到的数据仅能给出ROI 相对该测量点的二维角坐标信息,为了获得ROI 的三维信息,需要有两个或两个以上不同位置的测量点的数据.以两个测量点为例,当ROI 为简单的几何体(如长方体)时,可以很方便地将两个测量点观测到的ROI 边界角坐标对应起来.如图2,在测量点1 观测到ROI 的面abcd的角坐标为 (θ1,φ1),在测量点2 观测到面abcd的角坐标为 (θ1′,φ′1),则可结合两个测量点的位置信息,根据几何关系计算出面abcd的三维坐标,以此类推,最终重建出ROI 的几何大小和三维位置信息.

图2 测量点与ROI 之间的几何关系示意图Fig.2.Geometric relationship between viewpoint and ROI.

3 模拟系统建立

本文基于Geant4 模拟宇宙射线µ子在秦始皇陵中的输运过程.Geant4 是欧洲核子中心使用C++语言开发的一款模拟粒子输运过程的蒙特卡罗软件包,能提供构建几何模型、跟踪粒子径迹、模拟粒子与物质的相互作用等功能[21].本文模拟系统的建立主要包括3 个方面:秦始皇陵地宫模型、µ子探测器、海平面宇宙射线µ子源.

3.1 秦始皇陵地宫模型

本文根据使用综合地球物理方法得到的秦始皇陵考古数据建立了两个秦始皇陵地宫简化模型[22].其中模型1 为待测物体模型,包括了封土堆、土地和地宫,其中地宫由细夯土墙、宫墙、回填夯土和墓室组成;模型2 为参考模型,仅包括封土堆和土地.模型几何及关键结构的尺寸如图3 所示.表1 列出了地宫模型中不同区域的材质和密度,其中空气定义为70%的氮气和30%的氧气,黄土的组分定义参考了文献[23].

表1 秦始皇陵地宫模型材质及密度定义表[22]Table 1.Material and density definition table of the Qinshihuang Mausoleum model[22].

3.2 µ子探测器模型

本文采用了理想µ子探测器模型,参考常见的µ子探测器几何结构,将探测器模型表面积设为1 m × 1 m,其作用为提取到达探测器区域的µ子的速度方向信息.选取探测器位置时考虑了以下4 个因素:1)探测器位置的海拔高度需要低于地宫区域,且避开已知的文物埋藏区域;2)由于µ子通量随天顶角增大而显著减少,为保证µ子的计数率,应避免探测器对待测区域探测方向的天顶角过大;3)目标探测区域要尽可能离探测器近以提高计数率;4)目标探测区域内不同密度的待探测结构在探测方向上尽量没有重叠.为了给出成像目标的三维信息,选定两个测量点,如图3(a4)所示.测量点1 选取在µ子计数率最高的地宫正下方,探测器水平放置,埋深80 m;测量点2 选取在非垂直方向的地宫边缘下方,位于地宫中心正南方80 m 处,探测器埋深89.5 m,探测器表面与水平面成40°夹角,使探测平面法线方向指向墓室中心,从而增大此方向上的µ子计数率.

图3 秦始皇陵模型示意图 (a) 模型1 内部结构示意图;(a1) 模型1 俯视图;(a2) 模型1 正视图;(a3) 模型1 剖面3 示意图;(a4) 模型1 剖面1 示意图;(b) 模型2 示意图(无内部结构);Fig.3.Model of Qinshihuang Mausoleum:(a) Inner structure of Model 1;(a1) top view of Model 1;(a2) front view of Model 1;(a3) profile 3 of Model 1;(a4) profile 1 of Model 1;(b) Model 2 (no inner structure).

3.3 µ子源产生

模拟中µ子源根据海平面µ子能谱经验公式抽样产生.关于海平面µ子能谱分布的模型有很多种,例如Gaisser[24],Reyna[25],以及Smith 和Duller[26]等提出的模型.比较这些模型发现,Reyna 总结的海平面µ子能谱公式能在1 GeV/c ≤p≤2000/cosθGeV/c的范围内拟合全部天顶角范围内的µ子能谱分布,适合于µ子吸收成像中µ子源的抽样[27].其表达式为

I(p,θ) 为海平面µ子微分通量,单位为(GeV/c)-1·cm-2·sr-1·s-1;p为µ子动量,单位为GeV/c;θ为µ子速度方向天顶角大小;c1=0.00253 ,c2=0.2455,c3=1.288,c4=-0.2555 ,c5=0.0209 .

模拟中使用的µ子源通过对(4)式在p=30—1000 GeV/c,θ=0°—70°的范围内抽样产生,以减少对无法穿透土地到达探测器的低能µ子的模拟,所使用的µ子源的动量和天顶角分布如图4.

图4 根据Reyna 公式抽样产生的1000 万个µ子的动量和天顶角分布 (a) µ子数量随µ子动量变化分布;(b) µ子数量随µ子速度方向的天顶角变化分布Fig.4.Momentum spectrum and zenith angle distribution of the 10 million muons sampled by Reyna formula:(a) Momentum spectrum of the sampled muons;(b) zenith angle distribution of the sampled muons.

4 模拟结果

4.1 墓室二维成像结果

对模型1 和模型2 均模拟了等效于实际测量一年的µ子量,对两个测量点在模型1 和模型2 测得的µ子通量作差得到f(θ,φ) 的分布,如图5 所示,图中灰度值表示f(θ,φ)(单位:cm-2·sr-1·d-1) .

图5(a)中间的白色矩形区域f(θ,φ) ＞ 0,说明此方向上存在密度偏低的结构,即墓室,根据白色矩形关于图像中心的对称性可知测量点1 位于墓室正下方;周围的深色矩形区域f(θ,φ) ＜ 0,说明此方向上存在密度偏高的结构,即墓室周围的墙体; t anθy=0 附近深色矩形区域边缘处向内凹陷的部分对应于墙体的不连续部分.图5(b)中的白色扇形区域f(θ,φ) ＞ 0,对应于墓室;周围的深色弧形区域f(θ,φ) ＜ 0,对应于墙体;弧形区域与白色扇形区域之间的区域f(θ,φ) ＜ 0,对应于地宫开挖范围内回填的夯土.对比图5(a)和图5(b)可知,相对于测量点1,测量点2 除可以显示墙体、墓室外,还可以看出地宫开挖范围回填夯土所在区域,这是因为在测量点1 的探测方向上墙体总是与回填夯土区域有重叠,导致两者的位置在图5(a)所示的角度投影图上也是相互重叠的.

为了便于墓室三维重建,选择图3(a1)所示剖面1、剖面2 处的f(θ,φ) 进行分析,从而定位墓室墙壁在这两个剖面处的角坐标.如图5 所示,虚线表示剖面方向,圆点表示剖面处的墙壁的角坐标位置,其中点B(23.75°,270°)和E(38.8°,0°)对应于剖面1 处的南墙,点A(25.64°,90°)和F(58.0°,0°)对应于剖面1 处的北墙,点C(35.75°,180°)对应于剖面2 处的西墙,D(37.23°,0°)对应于剖面2处的东墙.

图5 两个测量点得到的 f (θ,φ) 的二维投影图 (a)测量点1 的 f (θ,φ) 投影图,其中,t anθx=tanθcosφ ,tanθy=tanθsinφ ;(b)测量点2 的 f (θ,φ) 投影图Fig.5.Two-dimensional projection of f (θ,φ) obtained at viewpoint 1 and 2:(a) Distribution of f (θ,φ) obtained at viewpoint 1,where the t anθx=tanθcosφ ,t anθy=tanθsinφ ;(b) distribution of f (θ,φ) obtained at viewpoint 2.

4.2 墓室三维重建结果

墓室的大小和位置是秦始皇陵考古研究最关心的问题之一,利用4.1 节中得到的墓室墙壁的角坐标,结合两个测量点的三维坐标可以重建墓室的三维位置.如图6(a)所示,由几何关系可计算出墓室南北方向长53.35 m,墓室中心埋深22 m.再由图5(a)中墓室东墙、西墙的角坐标结合墓室埋深可计算出墓室东西方向长度为85.82 m,如图6(b)所示.重建得到的墓室墙壁边长相较于理论值的差异约为7%,埋深差异约为6%.该重建结果相对于真值的差异与模拟中探测器位置、统计误差、算法等有关.在未来的实际测量中,还需要考虑探测器的角分辨率、几何接受度、探测效率,以及本底和噪声等对径迹重建和图像重建的影响.

图6 墓室三维重建结果 (a) 剖面1 处重建结果;(b)剖面2 处重建结果Fig.6.Three-dimensional reconstruction results of the chamber:(a) Reconstruction result at Profile 1;(b) reconstruction result at Profile 2.

5 总结

本文使用GEANT4,基于现有的秦始皇陵考古数据以及理想µ子探测器,对秦始皇陵地宫µ子吸收成像进行了模拟研究.成像结果初步验证了秦始皇陵地宫µ子吸收成像的可行性,单视角获得的µ子通量投影数据可以清晰地反映出地宫内部的不同结构,并能给出墓室墙壁的二维角坐标信息,利用两个视角下获得的通量数据可以重建出墓室大小和三维位置,重建得到的墓室边长和墓室中心埋深相对于理论值的差异在7%左右.本文仅是对µ子吸收成像应用于秦始皇陵地宫无损探测的初步研究,后续研究中将进一步细化秦始皇陵模型和µ子探测器模型,调整测量点的位置,增加测量点数量,优化多视角三维图像重建算法,深入分析影响图像重建结果的各种因素.