摘要:现有以散射为主的宇宙线缪子成像难以从高原子序数材料中有效区分特殊核材料, 利用缪子在特殊核材料中产生的次级诱发中子标记入射缪子, 可从高原子序数材料中辨别特殊核材料, 但其成像时间长且成像质量较低. 缪子多模态成像利用缪子穿透材料的散射信息以及被材料阻止时产生次级诱发中子的缪子信息, 可有效解决单一成像方法的不足. 基于GEANT4程序设置探测模型, 以Cosmic-ray Shower Library为缪子源, 开发了与诱发中子符合的缪子成像、缪子散射成像、缪子多模态成像模拟程序, 并在模拟数据的基础上分别实现了成像算法, 得到了不同模型的成像结果. 模拟结果表明, 基于散射和次级诱发中子的缪子多模态成像不仅成像快、质量高, 而且能明显区分特殊核材料与其他高原子序数材料, 具有探测特殊核材料的独特优势.
关键词: 宇宙线缪子/
缪子散射成像/
缪子诱发中子/
缪子多模态成像

English Abstract

全文HTML

--> --> -->

入射大气层的初级宇宙射线与大气原子核发生相互作用产生大量次级粒子, 其中大部分为数量相近的不同电性的π介子, 带电π介子又会衰变为缪子^[1]. 宇宙线缪子具有单位电荷, 质量约为质子质量的1/9, 电子质量的207倍, 静止寿命大约为2.2 μs^[2]. 到达海平面的缪子具有广泛的能量分布, 平均能量为3—4 GeV, 通量为1 × 10⁴ m^–2·min^{–1 [3]}. 缪子与物质主要有三种相互作用形式, 分别为能量损失、完全吸收和多次库伦散射^[4]. 自1936年安德森通过云室测量宇宙线粒子能损发现缪子以来^[5], 科研人员相继开发出具有不同探测器结构，适用于不同场景的透射成像和散射成像技术^[6,7], 并将其广泛应用于核不扩散、军控核查、火山监测、空间气象预测、反应堆监控等众多领域^[8,9]. 用作探针的宇宙线缪子作为广泛存在的自然源, 天然具有安全性和极高的穿透性, 能够很容易克服屏蔽问题^[10], 满足军控核查中特殊核材料(special nuclear material, SNM)检测成像的要求.
宇宙线缪子散射成像具有区分高、中、低原子序数材料的能力^[11], 但难以从高原子序数材料中分辨出SNM. 散射成像利用的缪子多重库伦散射信息取决于粒子穿透材料的辐射长度和厚度^[12], 这增加了缪子穿透低原子序数的厚材料与穿透高原子序数的薄材料得到的散射信息相当的可能性. 散射成像的图像重建方法基本还是融合、改进美国洛斯·阿拉莫斯国家实验室(Los Alamos National Laboratory, LANL)提出的最近邻点算法(the point of closest approach algorithm, PoCA)和最大似然散射角及位移算法(the maximum likelihood scattering and displacement algorithm, MLSD)^[13]. 诱发中子符合的缪子成像利用缪子在SNM中产生的次级诱发中子来标记入射缪子, 能够在不泄露核材料的具体数量与同位素组成的前提下, 产生与散射图像质量相匹配的信噪比, 以高穿透性和低剂量实现单侧成像^[14]. 国外对此进行了较为深入的研究: LANL的Guardincerri等^[15]利用宇宙线缪子产生的中子探测SNM的存在, 探索了诱发中子符合的缪子成像在军控核查中的可能应用; Morris等^[16]开发了作为缪子散射成像核材料补充技术的与诱发中子符合的缪子成像技术; Blackwell和Kudryavtsev^[17]利用蒙特卡罗方法模拟测量不同类型靶材料中缪子消失的结果, 证明了利用因产生诱发中子而消失的缪子可以增强缪子散射成像技术检测集装箱中屏蔽高浓铀的能力.
利用缪子与材料间多种相互作用信息进行图像重建及对不同材料进行识别, 称之为缪子多模态成像. 为充分利用缪子的检测及成像能力, 本文结合缪子散射信息以及产生次级诱发中子的缪子信息, 聚焦于不同原子序数材料的检测以及高原子序数材料中SNM的识别, 利用Geant4程序进行模拟研究, 以期通过模拟及数据处理实现基于散射和次级诱发中子的缪子多模态成像, 为进一步研究及实际应用提供可靠依据.

本文分别根据缪子散射成像、诱发中子符合的缪子成像以及综合利用二者信息的多模态成像对不同的探测模型进行成像. 基于Geant4程序, 缪子散射成像设计了两种以²³⁵U材料、²³⁹Pu材料为主的模拟研究场景; 诱发中子符合的缪子成像首先研究了天然宇宙线缪子能谱下HEU (90% ²³⁵U)板、LEU (20% ²³⁵U)板、DU (0.2% ²³⁵U)板的中子增益, 以及三种铀板的中子增益与缪子能量、带电性质的关系, 再利用与诱发中子符合的缪子分别对²³⁵U块、²³⁹Pu 块、铅块、铝块, 以及多层结构的核部件成像; 缪子多模态成像以有无铁屏蔽壳为区别, 设计了两种以检测高原子序数材料中SNM为目的的成像场景, 并将成像结果与单一的缪子散射成像、诱发中子符合的缪子成像结果相比较.
2

-->

2.1.缪子源项、探测模型

模拟中使用独立的物理程序包Cosmic-ray Shower Library (CRY)产生宇宙线缪子, 其功能是产生与宇宙线分布相符的粒子, 包括宇宙线缪子、中子、质子、电子、光子、π介子等^[18]. 利用该软件包获取的缪子能量、位置、方向等数据可以作为源项进行粒子输运和模拟. 探测模型如图1所示, 沿Z轴方向设置了四块边长为1 m的气体探测器平板, 缪子源从一侧探测器板外向另一侧入射, 入射位置为整个探测器平板, 根据研究内容分别在探测空间内放置不同的探测物体. 缪子多模态成像和散射成像需要利用四块探测器平板, 而诱发中子符合的缪子成像仅需入射位置的两块探测器平板即可达到单侧成像的目的.
图 1 探测模型设置
Figure1. Detecting model setting.

2 -->

2.2.缪子散射成像原理

缪子在穿过厚度为L的材料过程中会与材料原子核发生多重库伦散射, 使缪子最终以散射角θ和偏转位移x随机出射材料. 缪子在物质中发生多重库伦散射, 散射角θ近似均值为0, 标准差为σ_θ的高斯分布^[19]:

$ {f}_{\theta }(\theta )\approx \frac{1}{\sqrt{2\text{π}}{\sigma }_{\theta }}{\rm{exp}}\left( {-\frac{{\theta }^{2}}{2{\sigma }_{\theta }{}^{2}}} \right) , $

$ {\sigma _\theta } = \frac{{13.6{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{MeV}}}}{{\beta {\rm{c}}p}}\sqrt {\frac{L}{{{L_0}}}} \left[ {1 + 0.038\log \left(\frac{L}{{{L_0}}}\right)} \right], $

其中, L₀为材料的辐射长度; p为缪子动量; β_C为缪子速度, 对缪子β_C = 1. 缪子在材料中的位移也是近似均值为0, 标准差为σ_x的高斯分布:

$ {\sigma _x} = \frac{L}{{\sqrt 3 }}{\sigma _\theta }, $

辐射长度是电磁相互作用的特征量, 与原子序数Z和质量数A之间的关系如(4) 式所示:

$ {L_0} = \frac{{716.4{\kern 1pt} {\kern 1pt} ({\rm{g}} \cdot {\rm{c}}{{\rm{m}}^{{\rm{ - 2}}}})A}}{{Z(Z + 1)\ln (287/\sqrt Z )}}, $

辐射长度代表粒子能量损失到1/e时所穿过的平均材料长度. 当材料厚度与辐射长度比值不大于10 时, (2)式中对数项对散射角分布宽度的影响不超过4%^[20], 通常忽略对数项的影响. 结合(2)式、(4)式, Z的升高导致L₀下降, 从而使σ_θ增大. 因此测量的散射角分布实际反映了物质原子序数的大小.
为简化数学表达式, 定义与材料厚度、缪子动量和速度无关的散射密度λ, 得到散射角分布与散射密度的关系如(5)式所示, 式中p_r为能量矫正系数, p_r= p₀/p, p₀为名义动量.

$ {\sigma _\theta }^2 = \lambda L{P_r}^2, $

通过测量缪子穿过被测材料前后的径迹, 可以得到缪子偏转角度及位移等信息, 计算其统计分布, 再根据图像重建算法确定被测材料的散射密度, 就能得到材料属性及空间分布.
2 -->

2.3.诱发中子符合的缪子成像原理

缪子诱发裂变最初是惠勒^[21]在1948年提出的, 此后人们对它进行了更深入的探索, 并证明它在研究核能耗散和裂变动力学方面有重要作用^[22]. 遇到裂变物质的宇宙线缪子可能通过三种不同的机制产生诱发中子: 光致中子反应、缪子诱发的裂变和中子增益^[9]. 当高度相对论性的缪子通过高原子序数物质时, 会因发射轫致辐射光子而失去能量, 当重原子核与光子相互作用时, 瞬时光中子可能被释放出来. 失去足够能量的缪子会在物体内部停止, 带正电的缪子会衰变为一个正电子和两个中微子, 带负电的缪子被物质原子中的质子俘获后产生中子和中微子^[23]. 与其他效应相比, 负缪子诱发裂变产生中子的几率大得多(几个数量级), 且主要是低能的负缪子. 在Geant4中, 缪子俘获由物理过程 muMinusCaptureAtRest模拟.

$ {\mu ^ + } \to {{\rm{e}}^ + } + {v_u} + \overline {{v_e}} , $

$ {\mu ^ - } \otimes A \to n + {(A - 1)^*} + \overline {{v_\mu }} , $

在可裂变材料中, 上述过程产生的中子能够触发裂变链, 发射数个中子, 进而增强次级信号, 裂变事件的数量取决于裂变材料的形状、数量以及其他可能反射、吸收或慢化中子的材料. 中子输运方程定义了系统中子增益和损耗的平衡方程^[24]:

$ \frac{1}{v}\frac{{\partial \varPsi }}{{\partial t}} + \boldsymbol\varOmega \cdot \nabla \varPsi + {{{\Sigma }}_t}\varPsi = {Q_s}{\left[ {}\varPsi {} \right]} + {Q_f}{\left[ {}\varPsi {} \right]} + S , $

其中, v是中子速度, $\varPsi (\boldsymbol r, E, \boldsymbol\varOmega, t)$ 为中子角通量密度, ${\varSigma _t} = {\varSigma _s} + {\varSigma _a}$为宏观截面, 等式右侧为源项(散射源、裂变源、外源).
如图2所示, 当缪子在探测物质中因产生诱发中子消失后, 仅能在入射一侧探测器平板记录到两次与诱发中子符合的缪子坐标信息, 结合两个坐标信息得到缪子入射轨迹. 如果入射轨迹被外推到物体所在的平面上, 它们会使物体聚焦. 在物体中产生次级中子的缪子数越多, 缪子轨迹与成像平面的交点在物体位置的密度ρ就越大, 从而实现对探测物质的成像. 研究忽略了缪子入射被测物质前沿其轨迹可能出现的任何散射.
图 2 诱发中子符合的缪子成像图解(蓝色轨迹为被次级中子标记的入射缪子轨迹)
Figure2. Diagram of muon imaging in coincidence with muon induced neutrons (the blue trajectories are the incident muon trajectories tagged by muon induced neutrons).

2 -->

2.4.多模态成像原理

缪子散射成像时间较短、质量较高, 但只利用了穿过被检材料缪子的散射信息. 诱发中子符合的缪子成像具有区分高原子序数材料与SNM的独特优势, 但只利用了被材料俘获的产生了次级中子的缪子, 且产生次级中子的缪子数较少. 若能结合两种效应得到的有效信息, 可在单一信息成像的基础上, 实现基于散射和次级诱发中子的缪子多模态成像. 将成像平面划分为大小一致的网格, 多模态成像信息w_j由同一网格j内、不同权重的散射密度信息λ_j与产生诱发中子的缪子交点密度信息ρ_j加权得到.

$ {w_j} = i{\lambda _j} + k{\rho _j} , $

由于缪子主要在入射SNM时产生次级中子, 从而加强多模态成像信息w_j, 因此基于散射和次级诱发中子的缪子多模态成像在SNM的检测识别方面具有更大优势.

2

-->

3.1.模型设置

设置如图3(a)所示的单物块成像模型对缪子散射成像能力进行了初步检验, 在成像空间中心位置放置了边长为10 cm的²³⁵U块, 通过CRY发射200万个缪子. 进一步设置了如图3(b)所示的距中心等间距放置的边长均为10 cm的²³⁹Pu块、²³⁵U块、铅块、铝块, 以检验缪子散射成像在高原子序数材料中区分SNM的能力, 通过CRY发射500万个缪子.
图 3 散射成像模型　(a) 单物块模型; (b) 四物块模型
Figure3. Scattering imaging model: (a) Single cube model; (b) four cubes model.

2 -->

3.2.成像结果

缪子散射成像利用最近邻点算法对模拟得到的数据进行处理, 得到成像结果如图4和图5所示. 缪子散射成像平面大小为1000 × 1000, 像素大小取为20 × 20. 由图4可以看出, 缪子散射成像对单一²³⁵U块成像能力较好, 重建得到的²³⁵U块位置、大小都较为准确. 由图5可以看出, 缪子散射成像能够明显区分低原子序数铝块和高原子序数²³⁹Pu块、²³⁵U块、铅块, 但同为高原子序数材料的²³⁹Pu块、²³⁵U块、铅块的区分并不明显.
图 4 单物块模型散射成像结果
Figure4. Scattering imaging results of a single cube model.

图 5 四物块模型散射成像结果
Figure5. Scattering imaging results of four cubes model.

2

-->

4.1.三种铀材料的中子增益

长、宽均为50 cm, 厚10 cm铀板放置在探测空间中心位置. 首先研究了天然缪子能谱下不同²³⁵U丰度铀板的中子增益, 对HEU板、LEU板、DU板, 每次模拟都使用CRY发射10万个缪子, 统计出射到铀板外的中子数目, 模拟设置如图6(a)所示. 为进一步研究缪子能量、带电性质与产生中子数目的关系, 通过Geant4内置的通用粒子源(general particle source, GPS)分别设定10万个能量为10 MeV、1 GeV的正缪子、负缪子作为粒子源垂直入射到不同²³⁵U丰度铀板, 并统计出射到铀板外的中子数目, 模拟设置如图6(b)所示. 选择10 MeV能量的缪子是因为该能量处于宇宙线缪子能谱的低能侧, 使缪子因能量损失而完全停止在铀板中; 1 GeV则代表了宇宙线缪子经过30 cm的混凝土屏蔽后宇宙线缪子的平均能量^[25].
图 6 中子增益模拟图　(a) CRY缪子源; (b) GPS缪子源
Figure6. Neutron gain simulation diagram: (a) Simulation with CRY; (b) simulation with GPS.

模拟结果如图7所示, 天然宇宙线缪子能谱下HEU板统计得到的中子数目较LEU板和DU板均有大幅增加, 这是由于²³⁵U有更大的裂变截面, 使缪子诱发裂变产生了更多的次级中子. 在特定能量缪子打靶时, HEU板仍然较LEU板、DU板产生更多的次级中子; 在三种²³⁵U丰度的铀板中, 低能负缪子诱发裂变产生的中子数相比其他设定参数下产生的中子数具有数个量级的优势. 更强的次级信号增加了通过诱发中子标记与其符合的缪子的可能性, 为在高原子序数材料中区分SNM提供了可靠途径.
图 7 中子增益结果
Figure7. Result of neutronic gain.

2 -->

4.2.模型设置

设置了两种成像模型, 模型一与图3(b)所示模型相同, 仍通过CRY发射500万个缪子. 模型二如图8所示, 展示了根据国外典型核武器理论模型设计的简单核部件结构^[26-28], 由内到外依次为空腔、²³⁹Pu层、Be层、²³⁸U层、高能炸药层(黑索金C₃H₆N₆O₆)、Al外壳. 图8(a)为模型具体结构, 图8(b)为Geant4建模结果(为展示部件具体组成, 以半球显示). 核部件结构相对精细, 模拟使用CRY发射1000万个缪子.
图 8 诱发中子符合的缪子成像模型　(a) 核部件具体结构图; (b) 核部件建模示意图
Figure8. Imaging model of muon imaging in coincidence with muon induced neutrons: (a) Detailed structure diagram of nuclear components; (b) nuclear components model.

2 -->

4.3.成像结果

根据成像原理, 分别得到了两种模型在Z = 0平面的二维成像结果以及在不同Z值处断层扫描得到的三维成像结果, 如图9和图10所示. 二维成像平面大小为1000 × 1000, 像素大小为10 × 10, 第j个像素网格内ρ_j越大, 网格颜色越深, 三维成像由序列连续的多张二维图像融合得到. 图9(a)和图9(b)为模型一成像结果, 据二维成像图可知, 利用与诱发中子符合的缪子能够明显区分SNM与其他材料, 成像显示的四个物块所在位置以及²³⁹Pu块、²³⁵U块大小较为准确, 且²³⁹Pu块较²³⁵U块成像更加均匀, 物块边界也更加清晰; 三维成像结果较为立体地重建了²³⁹Pu块和²³⁵U块. 图10(a)和图10(b)为模型二的成像结果, 据二维成像图可知, 核部件整体形状成像较好, 外部的高能炸药和Al壳与内部的反射层分辨清晰, 最内层空腔与裂变材料²³⁹Pu层分辨清晰, ²³⁹Pu层与反射层边界较为模糊; 三维成像图由Z ≥ 0的二维断层重建图像得到, 成像结果显示了²³⁹Pu层与²³⁸U层之间的区域.
图 9 四物块模型的诱发中子符合的缪子成像图　(a) 四物块模型二维成像结果; (b) 四物块模型三维成像结果
Figure9. Muon imaging in coincidence with muon induced neutrons of four cubes model: (a) 2D imaging results of four cubes model; (b) 3D imaging results of four cubes model.

图 10 核部件模型的诱发中子符合的缪子成像图　(a) 核部件模型二维成像结果; (b) 核部件模型三维成像结果
Figure10. Muon imaging in coincidence with muon induced neutrons of nuclear components: (a) 2D imaging results of nuclear components; (b) 3D imaging results of nuclear components.

2

-->

5.1.模型设置

多模态成像设置了两种模型, 模型一如图11(a)所示, 沿X轴由上至下等间距放置了边长均为10 cm的²³⁵U块、铅块、铝块. 进一步设置了如图11(b)所示的屏蔽模型二, 以检验成像方法在有5 cm厚的铁屏蔽外壳存在的情况下, 探测SNM的能力. 模型一、模型二通过链接CRY分别发射200万个缪子和500万个缪子.
图 11 缪子多模态成像模型　(a) 三物块成像模型; (b) 铁外壳屏蔽的三物块成像模型
Figure11. Imaging model of multimodal imaging of muon: (a) Three cubes model; (b) three cubes model with iron shielding shell.

2 -->

5.2.成像结果

两种模型的缪子多模态成像结果均与单独的缪子散射成像以及诱发中子符合的缪子成像结果进行了对比. 据图12可知, 缪子散射成像能够较好地区分高原子序数的²³⁵U块、铅块和低原子序数的铝块, 但难以根据成像结果明显区分²³⁵U块和铅块; 诱发中子符合的缪子成像能够明显区别²³⁵U块和铅块, 并且低原子序数的铝块几乎不会被成像, 但其成像质量较差, 物块重建精度不够高; 缪子多模态成像在有效区分²³⁵U块和铅块的基础上, 实现了较为精确的图像重建, 与图9(a)中500万个缪子打靶物块得到的二维诱发中子符合的缪子成像结果相比, 200万个缪子打靶物块得到的多模态成像结果更加清晰. 由图13可以看出, 当有铁屏蔽壳存在时, 缪子散射成像得到的²³⁵U块与铁屏蔽外壳散射密度相当, 无法凭此检测出SNM; 诱发中子符合的缪子成像能够检测到屏蔽壳中存在SNM, 但成像质量较差; 缪子多模态成像能够明显区分²³⁵U块和铅块、铁屏蔽壳, ²³⁵U块的重建形状较为清晰.
图 12 三物块成像结果　(a) 缪子散射成像; (b) 诱发中子符合的缪子成像; (c) 缪子多模态成像
Figure12. Imaging results of three cubes model: (a) Muon scattering imaging; (b) muon imaging in coincidence with muon induced neutrons; (c) multimodal imaging of muon.

图 13 铁屏蔽下三物块成像结果　(a) 缪子散射成像; (b) 诱发中子符合的缪子成像; (c) 缪子多模态成像
Figure13. Imaging results of three cubes model with iron shielding shell: (a) Muon scattering imaging; (b) muon imaging in coincidence with muon induced neutrons; (c) multimodal imaging of muon

本文通过GEANT4设置了探测模型, 开发了用于缪子散射成像、诱发中子符合的缪子成像、缪子多模态成像的模拟程序, 并在程序中实现了通过链接CRY产生宇宙线缪子的功能. 在模拟程序的基础上, 分别实现了缪子散射成像、诱发中子符合的缪子成像以及缪子多模态成像. 根据研究方法进行了丰富的模拟实验, 由成像结果可知: 缪子散射成像能够较好地区分高原子序数材料和其它材料, 成像时间较短, 重建图像较为清晰, 但散射成像难以从高原子序数材料中区分SNM, 在有屏蔽外壳存在时, 成像效果较差; 诱发中子符合的缪子成像能够准确地从高原子序数材料中识别SNM, 能够在有屏蔽材料覆盖的情况下检测出SNM的存在, 可根据需要重建出二维或三维图像, 但由于次级中子数量较少, 导致成像时间较长且成像质量较低; 缪子多模态成像综合利用了穿透材料的缪子散射信息以及被材料阻止时产生次级中子的缪子信息, 有效融合了缪子散射成像快速、准确的优点和诱发中子符合的缪子成像检测SNM的独特优势, 能够在较短时间内准确重建未经屏蔽的SNM, 即使存在较厚的屏蔽层, 也能有效检测出SNM. 研究表明, 缪子多模态成像满足成像一般物体以及检测SNM的可行性, 能够扩展现有缪子成像应用场景并提高检测效率, 在军控核查领域和维护国土安全方面具有广泛的应用前景.
感谢南华大学王晓冬老师和季选韬、魏鑫同学对本文完成三维诱发中子符合的缪子成像提供的帮助.