μ子断层成像技术凭借其无需额外放射源、安全系数高、穿透能力强、对重核材料灵敏等众多优势成为了辐射成像领域一个新兴的研究热点[1-3]。在利用辐射成像手段进行目标对象结构诊断相关领域中，常需要研究重核材料（例如铀）及其元件经过一段时间放置或特定工艺处理后的内部结构变化问题，并且目标对象的初始结构组成通常是可以准确预知的。对于此类问题，传统的μ子散射成像反演法的基本思路是：分别对变化前、后两个目标对象进行μ子散射特征测量并利用特定算法反演出其内部三维结构图像，通过对比二者图像的差异即可获知目标对象之间的结构差异。这是一种直接、自然的方法，对于某些类型的目标对象也是行之有效的，但该方法往往涉及不确定性大、复杂性高的图像反演过程。如果目标对象组成结构较为复杂、尤其是具有大量孔洞、狭缝或低密度内填充物等各种类空心结构时，传统的成像反演方法极有可能受到严重的伪影干扰而难以给出清晰、准确的图像，因而使其应用范围大为受限。基于上述原因，亟需发展一种适用范围更广、实现难度更低的基于μ子散射特征来识别目标对象之间结构差异的新方法。
在这种背景下，我们提出了一种基于目标对象发生结构变化前后对天然宇宙射线μ子的散射特征变化来识别目标对象所发生的结构变化类型以及评估该特定结构变化的尺度的新方法。由于该方法不涉及逆向的图像反演（重建）过程，而是采用了与预先建立的参比数据库进行比对的方法。并且，由于该方法规避了繁琐且易被干扰的图像重建过程，与目标对象的初始结构复杂程度无关，因而对于重核材料对象结构变化诊断相关领域而言，其适用性要比传统成像反演方法更为广泛。
1 方法原理为了解决传统的μ子散射成像反演法用于特殊目标对象结构差异诊断上的难题，我们尝试了一种新的方法，如图 1所示。其基本思想是：首先，建立目标对象标准状态下的理论模型 (P1)，然后充分考虑目标对象所有可能发生的一系列结构变化（包括形变类型、尺度、位置等多个参量），建立一组或多组与之分别对应的微扰模型 (P2)。在这一步骤中，微扰模型体系必须具有很好的完备性。其次，通过蒙特卡罗模拟获得上述各个微扰模型对应的μ子散射特征量 (P3)，然后将这些散射特征量与相应的结构变化形成一一对应的映射关系，从而组成一套完备的参比数据库待用 (P4)。P1-P4这几个步骤（图 1中虚线框内）构成了该方法的先决条件。而在日常应用（步骤1-3）中，对于实际待测目标对象，采用μ子成像平台实验测量方法获得其μ子散射特征量，然后将该散射特征与数据库中的数据进行比对，寻找出相似度最高的条目，即可在给定的显著性水平α下，判定待测目标对象发生了相应的结构变化。由于该方法摒弃了较为繁琐、易引入误差的图像反演过程，试图通过直接比较两种状态下的目标对象对μ子的散射特征的差异来获取目标对象的结构变化信息，因而仅涉及正向的散射测量过程，故称之为“正向参比法”。
|
图 1 正向参比法的实现及其所涉及的关键技术 Figure 1 Implemention of forward comparison method and major techniques involved |
以水平放置的铀板的GEANT4模拟结果为例，对正向参比法用于铀板内部水平狭缝宽度识别的实现以及其中的关键技术作详细说明。
2.1 μ子散射特征量参比数据库的建立 2.1.1 微扰模型的建立本文采用的铀板狭缝模型如图 2所示，左、右两边分别为理想铀板模型和基于该模型扩展得到的一组分别带有宽度为di（1≤i≤N，N为指定的有限正整数）的狭缝的铀板模型。为简单起见，假定铀板只可能产生水平方向的狭缝，并且所有可能的狭缝宽度的集合D={di}为有限集，至少是可以近似看作有限集。由于微扰模型必须具有完备性，因此需要针对集合D中的每一个元素di分别设计相应的狭缝模型，然后进行后续的μ子散射蒙特卡罗模拟。
|
图 2 铀板狭缝模型 Figure 2 Schematic illustration of U slabs with a series of slit widths |
为了获得足够多的μ子散射径迹，我们编制了一套基于GEANT4的μ子成像系统来模拟μ子在目标对象中的多重库仑散射 (Multiple Coulomb Scattering, MCS) 行为。天然宇宙射线μ子源项采用CRY软件包[4]来生成。模拟计算中考虑了μ子与材料的库仑散射、电离以及轫致辐射等主要物理作用。
2.1.3 数据库的构建根据Moliere经验公式，μ子散射极角概率分布函数f(θ) 对不同目标对象具有较好的特异性，故本文选取f(θ) 作为μ子散射特征量。然而，由于天然μ子能谱范围极广、入射方向各异，f(θ) 往往难以采用精确的解析函数进行描述。并且由于μ子与物质多重库仑散射的随机性，即使对于同一目标，两次实验得到的散射角分布也可能存在差异。考虑这两方面的因素，实际应用中，需要进行同一条件下的多次重复实验来构建总体f(θ) 的多个独立样本。
2.2 目标对象结构变化的判定——基于假设检验的二元分类器根据上述分析，利用正向参比法评估目标物体结构变化的核心问题可归结为待测目标对象μ子散射特征量与预建数据库中μ子散射特征量之间的比对问题。考虑到实际采用的是多次重复实验获得的多个相互独立的散射角分布样本，因此问题的实质是：在总体分布形式f(Θ) 未知的条件下，检验两个独立样本是否取自同一总体的非参数检验问题。这里采用比较累积频数分布函数的Kolmogorov-Smirnov (KS) 检验方法[5]，其具体过程如下：
令Θ为单个μ子的散射角度量：
| $\mathit{\Theta } = {\left( {\Delta {\theta _x}} \right)^2} + {\left( {\Delta {\theta _y}} \right)^2}$ | (1) |
式中：Δθx、Δθy分别为XZ和YZ平面上的累积偏转角。假定两次独立实验的目标物体对应的天然μ子散射角总体分别服从f1(Θ) 和f2(Θ) 分布，则两次实验测得的散射角可以看作分别来自这两个总体的样本。考虑假设检验问题：
| $\begin{array}{*{20}{l}} {{H_0}:{f_1}\left( \mathit{\Theta } \right) \equiv {f_2}\left( \mathit{\Theta } \right),\quad {\rm{Same}}\;\;{\rm{object}}}\\ {{H_1}:\;{f_1}\left( \mathit{\Theta } \right) \ne {f_2}\left( \mathit{\Theta } \right),\quad {\rm{Different}}\;\;{\rm{object}}} \end{array}$ | (2) |
在给定的显著性水平α下，计算KS检验的p值（零假设H0可被拒绝的最小显著性水平）。如果p≤α，则拒绝H0，认为两次实验目标对象存在差异；反之，如果p > α，则接受H0，认为两次实验目标对象完全相同。
基于上述KS检验对H0拒绝/接受的结论，可以构建一个判定目标对象是否发生变化的二元分类器。令显著性水平α在一定范围内变化，即可在接受者操作特征 (Receiver Operating Characteristic, ROC) 空间内绘制出相应的ROC曲线。ROC曲线可以对该分类器的性能进行可视化，给出分类器的真阳性率 (True Positive Rate, TPR)（收益）与假阳性率 (False Positive Rate, FPR)（开销）之间的相对平衡关系[6-7]。下文中，我们会通过ROC曲线来评估该基于KS检验的二元分类器性能。
3 结果与讨论 3.1 不同μ子通量下可识别的最小狭缝宽度由于所用μ子通量不同会导致散射样本容量不同，进而可能影响假设检验的判定结果，因此有必要对上述方法在各种不同μ子通量条件下对目标对象是否存在差异的识别能力进行研究。
这里以厚度1 cm的铀质均匀平板作为研究对象，在平板中心水平面上分别设置宽度为d1= 0.1mm、d2=0.2 mm、d3=0.5 mm和d4=1.0 mm的4种水平狭缝。以40组理想铀板（狭缝宽度d0=0）的μ子散射模拟结果获取的40个散射角分布样本，以及40组带有宽度di（i为[1, 4]范围内某一指定正整数）狭缝的铀板的μ子散射模拟结果获取的40个散射角分布样本共同组成参比数据库，以带有宽度di狭缝的铀板作为待测目标对象，依次选取μ子通量为50 cm-2、125 cm-2、250 cm-2和1000 cm-2，分别采用上述正向参比法考察在一定的显著性水平α下，各种μ子通量条件下该方法对狭缝的识别能力。对于给定的α，一对 (TPR, FPR) 就对应了ROC空间的一个点。在 (0, 1) 范围内改变α，就可以得到一系列ROC空间点组成的一条ROC曲线，如图 3-5所示。根据ROC曲线的含义可知，曲线越靠近左上角，即曲线下面积 (Area Under Curve, AUC) 越接近1，相应的分类器性能也越好。由图 3可知，50cm-2的μ子通量下该分类器已可实现对宽度低至1.0 mm的水平狭缝的完美分辨（AUC约为1），相应的阈值p ≈ 0。
|
图 3 μ子通量为50 cm-2时，分类器对于带有宽度di的水平狭缝的铀板的识别性能ROC曲线 (a) d1=0.1 mm，AUC 0.60，(b) d2=0.2 mm，AUC 0.75，(c) d3=0.5 mm，AUC 0.96，(d) d4=1.0 mm，AUC 1.00 Figure 3 ROC curves for the FCM classifier to identify U slabs with a slit width of di at the muon flux of 50 cm-2 (a) d1=0.1 mm, AUC 0.60, (b) d2=0.2 mm, AUC 0.75, (c) d3=0.5 mm, AUC 0.96, (d) d4=1.0 mm, AUC 1.00 |
随着进一步增大入射μ子通量，每个散射样本的容量也会增大，从而对判定结果也可能产生影响。图 4中结果显示，μ子通量增大到125cm-2时，选择适当的阈值p，该分类器可近似完美识别最小宽度为0.5mm的狭缝（AUC约为1）。
|
图 4 μ子通量为125 cm-2时，分类器对于带有宽度di的水平狭缝的铀板的识别性能ROC曲线 (a) d1=0.1 mm，AUC 0.55，(b) d2=0.2 mm，AUC 0.86，(c) d3=0.5 mm，AUC 1.00 Figure 4 ROC curves for the FCM classifier to identify U slabs with a slit width of di at the muon flux of 125 cm-2 (a) d1=0.1 mm, AUC 0.55, (b) d2=0.2 mm, AUC 0.86, (c) d3=0.5 mm, AUC 1.00 |
图 5表明，当μ子通量增大到1000 cm-2时，选择适当的阈值p，该分类器可较完美识别最小宽度为0.2 mm的狭缝（AUC约为0.99）。
|
图 5 μ子通量为1000 cm-2时，分类器对于带有宽度d的水平狭缝的铀板的识别性能ROC曲线 (a) d1=0.1 mm，AUC 0.83，(b) d2=0.2 mm，AUC 0.99 Figure 5 ROC curves for the FCM classifier to identify U slabs with a slit width of di at the muon flux of 1000 cm-2 (a) d1=0.1 mm, AUC 0.83, (b) d2=0.2 mm, AUC 0.99 |
根据§3.1的结果，在μ子通量1000cm-2的条件下，采用正向参比法可完美分辨出1cm铀板内宽度最小可达0.2mm的狭缝。然而，在实际应用中，除了需要对目标对象是否产生狭缝进行简单定性的Y/N判断以外，目标对象的结构变化可能性比较多，即使对于同一类型和位置的变化（如中心层的水平狭缝），往往同时存在多种尺度（如狭缝宽度）的可能性。为了评估这种情况下正向参比法对目标对象结构变化的识别性能，构建如下场景：设计狭缝宽度分别为d2、d3、d4的三个铀板作为待测目标对象，以狭缝宽度分别为d2、d3、d4的铀板各40组（共计120组）的模拟结果获取的120个散射角分布样本组成参比数据库，选取μ子通量为1000 cm-2，考察该情形下正向参比法对于不同狭缝宽度的识别能力。其结果如图 6所示。
|
图 6 μ子通量1000 cm-2条件下，铀板目标对象狭缝宽度同时存在图中三种可能时，正向参比法对狭缝宽度di的目标对象的识别性能ROC曲线 (a) d2=0.2 mm，AUC 0.99，(b) d3= 0.5 mm，AUC 1.00，(c) d4=1.0 mm，AUC 1.00 Figure 6 ROC curves for the FCM classifier to discriminate three coexistent U slabs, each with a slit width of di at the muon flux of 1000 cm-2 (a) d2=0.2 mm, AUC 0.99, (b) d3= 0.5 mm, AUC 1.00, (c) d4=1.0 mm, AUC 1.00 |
由图 6可知，在μ子通量达到1000 cm-2的条件下，对于任意一种狭缝宽度为di(i=2, 3, 4) 的铀板对象来说，其ROC曲线的AUC值（分别为0.99、1和1）都接近1的理想值。由此可见，在该铀板对象结构变化的可能性被限定为上述三种宽度的水平狭缝的前提下，正向参比法可以较为理想地实现其中任意一种宽度的狭缝与另外两种狭缝的区分与识别。这表明，在适当的限定条件下，正向参比法可以实现目标对象特定结构变化尺度的准确评估。
4 结语相对于经典μ子成像方法用于特殊对象结构探测的技术难点，本文提出的正向参比法更适用于初始结构已知的特殊对象内部结构变化探测，有望实现更好的空间分辨率和时间效率。从带有0.1-1.0mm范围内不同宽度水平狭缝的1 cm厚度铀板模型的GEANT4模拟结果来看，在合理的天然μ子通量条件下，该方法对于目标对象的微小结构变化尺度进行准确判定和评估，未来有望在复杂大型核对象的无损探测等相关领域取得显著成效。
| [1] | Borozdin K, Hogan G, Morris C, et al. Surveillance:radiographic imaging with cosmic-ray muons[J]. Nature, 2003, 422(6929): 277. DOI: 10.1038/422277a |
| [2] | Schultz L, Borozdin K, Gomez J, et al. Image reconstruction and material Z discrimination via cosmic ray muon radiography[J]. Nuclear Instruments and Methods in Physics ResearchA, 2004, 519(3): 687–694. DOI: 10.1016/j.nima.2003.11.035 |
| [3] | Morris C, Alexander C, Bacon J, et al. Tomographic imaging with cosmic ray muons[J]. Science & Global Security, 2008, 16(1-2): 37–53. DOI: 10.1080/08929880802335758 |
| [4] | Hagmann C, Lange D, Wright D. Monte Carlo simulation of proton-induced cosmic-ray cascades in the atmosphere (UCRL-TM-229452)[R]. USA:Lawrence Livermore National Laboratory, 2007. |
| [5] |
以恒冠. 面向特殊对象的宇宙线缪子散射成像方法研究[D]. 北京: 清华大学, 2015.
YI Hengguan. Research on cosmic ray muon scattering tomography methods for special objects inspection[D]. Beijing:Tsinghua University, 2015. |
| [6] | Fawcett T. An introduction to ROC analysis[J]. Pattern Recognition Letters, 2006, 27: 861–874. DOI: 10.1016/j.patrec.2005.10.010 |
| [7] |
邹洪侠, 秦锋, 程泽凯, 等. 二类分类器的ROC曲线生成算法[J].
计算机技术与发展, 2009, 19(6): 109–112.
ZOU Hongxia, QIN Feng, CHENG Zekai, et al. Algorithm for generating ROC curve of two-classifier[J]. Computer Technology & Development, 2009, 19(6): 109–112. |