缪子散射的成像方法及成像系统

技术领域

本发明涉及宇宙射线成像的技术领域，尤其是涉及一种缪子散射的成像方法及成像系统。

背景技术

入射大气层的初级宇宙射线与大气原子核发生相互作用产生大量次级粒子，其中大部分为数量相近的不同电性的π介子，带电的π介子又会衰变为缪子。

作为一种穿透性极强的高能带电粒子，缪子可以很轻易地穿透厚的屏蔽层，其静止质量为105.7MeV，约为电子质量的207倍，既可以带正电，也可以带负电。缪子为带电粒子，质量介于电子和质子之间，所以缪子与物质发生相互作用时，μ

自1936年安德森通过云室测量宇宙线粒子能损发现缪子以来，缪子成像可以应用于一些独特的领域：(1)可以对超大型结构，如火山、金字塔等进行成像；(2)对集装箱是否含高Z物质进行检测；(3)对核反应堆等难穿透设施进行无损成像探伤。

关于缪子成像的反演方法，有透射方法、最邻近点(Point of Closest Approach，PoCA)方法及简化改版、统计类方法和快速检测方法等。目前运用缪子进行成像时，缪子散射事件与探测板大小和面积、天顶的夹角以及总探测时间等因素相关，成像数据往往很有限，在此条件下，这些成像方法往往面临成像质量低，成像结果不稳定等问题。

发明内容

为了解决现有技术存在的上述问题，本发明提供了一种缪子散射的成像方法，在基于体素划分和成像值重建的基础上，本方法采用体素划分后，对路径长度和散射角度等参数进行大体素化合并，并结合角度盖帽法，对体素内成像值进行重建，以提高缪子成像的清晰度。本发明还提供了采用上述成像方法的成像系统。

为了实现上述目的，本发明提供了如下技术方案：

缪子散射的成像方法，包括以下步骤：

(a)将缪子经过的成像区域的三维空间划分为体素；

(b)在每次散射事件中，

(1)获取缪子经过成像区域最邻近点位置，

(2)根据缪子经过成像区域的入射径迹和出射径迹，计算散射角度初值；将散射角度初值与散射角度阈值进行比对：

当散射角度初值大于散射角度阈值时，散射角度值统计为散射角度阈值的平方；

当散射角度初值不大于散射角度阈值时，散射角度值统计为散射角度初值其本身的平方；

(3)计算并统计从入射点到最邻近点所经过每个体素的每一段路径长度，以及计算并统计从最邻近点到出射点所经过的每个体素的每一段路径长度；

(c)根据所统计的路径长度和散射角度值进行大体素化合并，所述大体素化合并为对中心小体素临近网格内的体素数据进行合并；

(d)根据大体素化合并后的散射角度值和路径长度，计算得到成像密度和相应的成像值，并进行图像重建。

在本发明的一个实施例中，在获取最邻近点位置后，计算得到最邻近点所在的体素坐标。

在本发明的一个实施例中，在将成像区域的三维空间进行体素划分时，在x、y、z方向上分别进行均匀划分，从而划分为呈小长方体的体素。

在本发明的一个实施例中，所述散射角度阈值的典型值为收集并计算得到的散射角度初值统计标准差。

在本发明的一个实施例中，所述计算散射角度初值并与散射角度阈值进行比对包括步骤如下：

(b1)在每次散射事件中，探测缪子经过体素时入射径迹对应的入射直线，以及出射径迹对应的出射直线；

(b2)根据入射直线和出射直线，计算每一个散射角度初值；

(b3)统计散射角度初值，并计算散射角度初值集合的统计标准差，作为散射角度阈值；

(b4)遍历所有散射角度初值，并与散射角度阈值进行比对；

(b5)以此散射角度阈值的典型值成像后，将散射角度阈值以30％增大或缩小，再次成像，如此重复数次，以选出最后成像结果。

在本发明的一个实施例中，在获取和统计路径长度、散射事件和散射角度值时：

分配路径数组L，用于统计并存放缪子经过成像区域每个体素时，从入射点到最邻近点所经过每个体素的每一段路径长度，以及从最邻近点到出射点所经过的每个体素的每一段路径长度；

分配事件数组N

分配角度数组Θ，用于统计并存放散射角度值。

在本发明的一个实施例中，在大体素化合并过程中，路径数组L、事件数组N

其中，路径数组L的第一次大体素化合并为：

对于整个N

并用L

在对路径数组L的第一次大体素化合并后，再进行第二次大体素化合并：

对每个维度进行分离，在x,y,z方向依次进行大体素合并操作：

并用

在本发明的一个实施例中，所述成像密度为所述角度数组Θ除以路径数组L，数组的除法为逐元素相除。

在本发明的一个实施例中，根据事件数组N

本发明还提供了一种缪子散射的成像系统，包括；

第一组探测器，其位于成像物体的一侧，用于测量得到缪子入射成像区域的入射径迹数据；

第二组探测器，其用于成像物体的另一侧，用于测量得到缪子从成像区域出射的出射径迹数据；

存储器，用于存储入射径迹数据、出射径迹数据、入射路径数组L、事件数组N

处理器，其用于接收入射径迹信息和出射轨迹信息，并计算和统计缪子经过体素的路径长度、散射事件和散射角度信息，并根据上述信息获取成像值，进行图像重建。

基于上述的技术方案，本发明取得的技术效果为：

本发明提供的缪子散射的成像方法，先按规则对体素进行细分，计算并统计路径长度、散射事件和散射角度值等数据，之后再将中心小体素临近网格内的数据进行大体素化合并，作为中心小体素的成像值，以增加数据量，降低数据的涨落，提高成像质量。

本发明的成像方法在进行体素内成像值重建时，在大体素化合并的基础上，结合了散射角度盖帽的方法，以去除大于散射角度阈值的散射角度初值；即在大于散射角度阈值的散射事件中，其对于成像值的贡献，最多只计算相当于散射角度为某个盖帽值(散射角度阈值)对成像值的贡献，这样减少大角度事件对于成像值稳定性的破坏作用。

附图说明

图1为本发明的缪子散射的成像方法所使用的探测系统示意图。

图2为本发明的缪子散射的成像方法的流程图。

图3为本发明的缪子经过成像区域的示意图。

图4为本发明的缪子散射的成像系统的结构框图。

图5为采用本发明的缪子散射成像方法与传统的PoCA的方法对U型钨块的成像对比图。

具体实施方式

为了便于理解本发明，下面将结合附图和具体的实施例对本发明进行更全面的描述。附图中给出了本发明的较佳实施方式。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施方式。相反地，提供这些实施方式的目的是使对本发明的公开内容理解的更加透彻全面。

需要说明的是，当元件被称为“固定于”另一个元件，它可以直接在另一个元件上或者也可以存在居中的元件。当一个元件被认为是“连接”另一个元件，它可以是直接连接到另一个元件或者可能同时存在居中元件。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施方式的目的，不是旨在于限制本发明。

实施例1

图1为本实施例的缪子散射的成像方法所使用的探测系统示意图，图2为本实施例的缪子散射的成像方法的整体流程图，结合参考图1和图2，本实施例提供了一种缪子散射的成像方法，包括以下步骤：

(a)将缪子经过的成像区域的三维空间划分为体素；

(b)在每次散射事件中，

(1)获取缪子经过成像区域最邻近点位置，

(2)根据缪子经过成像区域的入射径迹和出射径迹，计算散射角度初值；将散射角度初值与散射角度阈值进行比对：

当散射角度初值大于散射角度阈值时，散射角度值统计为散射角度阈值的平方；

当散射角度初值不大于散射角度阈值时，散射角度值统计为散射角度初值其本身的平方；

(3)计算并统计从入射点到最邻近点所经过每个体素的每一段路径长度，以及计算并统计从最邻近点到出射点所经过的每个体素的每一段路径长度；

(c)根据所统计的路径长度和散射角度值进行大体素化合并，所述大体素化合并为对中心小体素临近网格内的体素数据进行合并；

(d)根据大体素化合并后的散射角度值和路径长度，计算得到成像密度和相应的成像值，并进行图像重建。具体来说，在成像物体的两侧分别设置了一组探测器，每组探测器具有两个探测器单元，如图1所示，位于成像物体一侧的两个探测器单元上的两个点p

(1)在成像过程中，需要获取和统计路径长度、散射事件和散射角度值等参数。因此，本实施例分配了路径数组L、事件数组N

其中，路径数组L和角度数组Θ为浮点型数组，事件数组N

如图2的整体流程图所示，先将成像物体区域的三维空间细致划分为多个体素；具体方法是，在x方向的成像范围(x

在y方向和z方向上按x方向的均分方式进行同样操作，从而将成像区域划分成的每一个小长方体，这些小长方体便称为体素。在一些实施例中，小长方体的各边可相等。

成像区域的三维空间内，可划分为N

需要说明的是，成像物体的成像区域范围要比所有体素的集合要略大一些，即成像区域范围除包含被测(成像)物体外，还需要留出一定的边缘。

对于成像区域三维空间内的N

本实施例直接用三维数组为例，路径数组L的坐标从0开始。在分布完数组之后，将路径数组L所有元素初始化为零；事件数组N

(2)接下来，按以下方法获取和计算得到缪子经过成像区域每个体素时的最邻近点(Point of closest Approach,PoCA)：

缪子的入射径迹为p

其中，

再求出p

与q

取p

由PoCA的位置坐标(x,y,z)计算其所属的体素位置坐标(i,j,k)，

以上为公式组(1)。

体素位置坐标(i,j,k)存储于路径数组L中。

(3)在获取缪子经过成像区域每个体素时的最邻近点(PoCA)位置坐标，以及对应的每个体素位置坐标后，需要计算并统计从缪子入射点到最邻近点所经过每个体素的每一段路径长度，以及计算并统计从最邻近点(PoCA)到缪子出射点所经过的每个体素的每一段路径长度。在本实施例中，将经过成像区域的缪子路径假定为经过PoCA的折线。

具体来说，如图3为的缪子经过成像区域的示意图所示：

(3.1)获取p

计算p

具体计算方法为：

先采用参数方程标记直线p

当t＝0时，对应起点为p

当t＝1时，对应终点p

求出

对{t

其中，点

经过以上处理，计算并统计得到从缪子入射点到最邻近点所经过每个体素的每一段路径长度。

(3.2)获取q

采用与(3.1)“计算p

(4)请结合上述(2)关于“获取和计算得到缪子经过成像区域每个体素时的最邻近点(Point of closest Approach,PoCA)”的描述：

缪子的入射径迹为p

其中，

计算入射直线p

将计算得到的散射角度初值θ统计，并计算散射角度初值集合的统计标准差，作为散射角度阈值c(也可称为角度盖帽值)的典型值。

在实际处理操作过程中，c的重要参考值是散射角度初值的典型值，比如可统计标准差，如本实施例所处理那样。由于散射角度阈值c过大，会导致对涨落的限制效果不好，而散射角度阈值c过小，又会导致成像值之间的区分度不好。因此，在一些实施例中，对于散射角度阈值c的选取可由散射角度初值的统计标准差开始，在其周围如四分之一至四倍范围内多测量几个值，以获得最好的成像效果。

在角度数组Θ[i,j,k]中，遍历所有散射角度初值θ，并与散射角度阈值c进行比对。

当散射角度初值θ大于散射角度阈值c时，最终统计入角度数组Θ的散射角度值α，为散射角度阈值c的平方，即将散射角度值α＝c

当散射角度初值θ不大于散射角度阈值c时，最终统计入角度数组Θ的散射角度值α，为散射角度初值θ的平方，即将散射角度值α＝θ

以上为本实施例采用的角度盖帽法，可限制θ过大对于角度数组Θ[i,j,k]的破坏性作用。

需要说明的是，在上述(2)关于“获取和计算得到缪子经过成像区域每个体素时的最邻近点(Point of closest Approach,PoCA)”过程中，由PoCA的位置坐标(x,y,z)计算其所属的体素位置坐标(i,j,k)后，对记录事件次数的事件数组N

(5)还需要说明的是，对于所有事件以及数组累加的处理，可采取并行的方式。在处理过程中，本实施例记录每个体素中对于每个事件信息的路径数组L、事件数组N

待处理完所有事件后，将所有子任务的L、N

(6)在处理完所有每次事件后，接下来对得到的各个数组进行大体素化处理以得到最后的图像。

(6.1)对路径数组L进行第一次第大体素化合并。在整个N

并用L

(6.2)本实施例在对路径数组L的第一次大体素化合并后，再进行第二次大体素化合并(优化)：在路径数组L的第一次大体素化过程中，一共有(N

在第二次大体素合并中，将此每个维度进行分离，在x,y,z方向(顺序不必固定)依次进行大体素合并操作：

并用

(6.3)在路径数组第二次大体素化合并过程中，求得

因此可以将这一步骤的元素访问次数由N

将三个步骤的计算量相加，合并后的总计算量为：

6N

在N

在上述进行每个维度的一维合并时，如对x维度进行合并操作，则y和z维度各元素相互不依赖，应当并行计算，以提高速度。上述的大体素合并过程，可适合GPU并行计算，当运用于实时成像时，将起到优秀的提速效果。

(6.4)需要说明的是，对事件数组N

(6.5)根据大体素化合并后的路径数组和角度数组，计算得到成像密度和相应的成像值，并进行图像重建。

具体来说，成像密度λ＝Θ/L，角度数组与路径数组之间为逐元素相除。

(7)结合事件数组N

对上述成像值去掉的依据是缪子散射成像装置的构造导致在成像区域边缘探测到的事件数会较成像中心区域更低，且操作人员会把待成像物体放置在中心区域。

在边缘如果出现少数大角度散射事件(散射角度初值大于散射角度阈值)，会造成成像值可能非常大，大于中心区域成像值的最大值，从而使中心区域的成像的对比度降低，影响成像质量。

因此，根据事件总数、成像区域大小、大体素尺寸等因素的不同，最小的事件次数从数次至数十次以上，合适的值要满足两个标准：(一)此值足够大，以使作图后成像区域边缘无物质区域没有明显的高亮度点成像值即可；(二)此值不可过大，使得中心成像区域遭遇到了去除。满足此条件即可，不必过于追求精确。在一些实施例中，对于事件次数为数千次的三维成像，可以从5开始看成像是否合适而调大此值至满足标准(一)。

实施例2

图4为本实施例的缪子散射的成像系统的结构框图，如图4所示，本实施例提供了一种缪子散射的成像系统，包括；第一组探测器，其位于成像物体的一侧，用于测量得到缪子入射成像区域的入射径迹数据；第二组探测器，其用于成像物体的另一侧，用于测量得到缪子从成像区域出射的出射径迹数据；存储器，用于存储入射径迹数据、出射径迹数据、入射路径数组L、事件数组N

实施例3

采用实施例1的缪子散射的成像方法，以及实施例2的缪子散射的成像系统，在N

再按照实施例1提供的方法进行计算，合并后在z(或y、x)方向所有体素信息都相同，实际为二维截面图。

此外，也可以先将各数组的z(或y、x)方向先作求和合并，生成二维数组，再对二维数组进行大体素合并操作，以进一步节省计算资源。

图5为采用实施例1的缪子散射成像方法与传统的PoCA的方法对U型钨块的成像对比图，如图5所示，左图为采用传统PoCA方法对U型钨块的成像图5A，右侧为采用实施例1的大体素结合角度盖帽法对U型钨块的成像图5B。可以看出，传统PoCA方法的不同成像值随空间的涨落很大，有个别超高亮度点影响整体成像效果，在有物体区域也存在低亮度点，在无物体区域也存在高亮度点，最终导致呈现物体轮廓不清晰，图像质量不佳。

而本实施例所采用的缪子散射成像方法，可有效提高成像的清晰度，并降低成像值在每个像素内的涨落和空间不同位置的涨落，同时降低对于原始入射轨迹与出射轨迹的准确度要求，图像质量明显优于传统PoCA方法的成像质量。

以上内容仅仅为本发明的结构所作的举例和说明，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些显而易见的替换形式均属于本发明的保护范围。