CN101606082B - 基于带电粒子测量的统计层析重建 - Google Patents

Abstract

用于带电粒子探测的系统、装置、计算机程序产品和方法包括带电粒子层析数据的物体体积散射密度剖面的统计重建，以使用统计多重散射模型确定带电粒子散射的概率分布，并使用期望最大化(ML/EM)算法确定物体体积散射密度的实质最大似然性估计，以重建物体体积散射密度。能够从所重建的体积散射密度剖面中识别占据感兴趣的体积的物体的存在和/或类型。带电粒子层析数据可以是来自扫描包裹、集装箱、车辆或货物的μ介子跟踪器的宇宙射线μ介子层析数据。能够使用可在计算机上执行的计算机程序实现这样的方法。

Description

基于带电粒子测量的统计层析重建

临时申请的交叉引用

本PCT申请要求(1)申请号为60/855,064、标题为“SYSTEMS，METHODS AND APPARATUS FOR PARTICLE DETECTION ANDANALYSIS AND FIELD DEPLOYMENT OF THE SAME”、并于2006年10月27日提交的美国临时专利申请和(2)申请号为11/771,169、标题为“RADIATION PORTAL MONITOR SYSTEM AND METHOD”、并于2007年6月29日提交的美国临时专利申请的优先权。

通过引用将以上两件申请的公开内容并入这里。

关于联邦权利的声明

本发明是在由美国能源部授予的、合同号为DE-AC52-06NA25396的政府支持之下进行的。在本发明中政府具有一定的权利。

技术领域

实施例涉及粒子探测、分析和控制领域，并且更加具体地而非排他地说，涉及方法和系统，所述方法和系统用于分析来自具有多个位置敏感探测器的带电粒子探测系统的数据，以及用于重建穿过带电粒子探测系统的、诸如宇宙射线μ介子的带电粒子的轨迹。

背景技术

带电粒子层析(charged particle tomography)基于带电粒子的散射(scattering)。带电粒子层析的一种形式是依赖宇宙射线μ介子(cosmic raymuon)的散射的宇宙射线层析。来自外层空间的稳定粒子大部分为质子，这些粒子持续不断地轰击地球。所述粒子与在高层大气中的原子相互作用，以产生包括许多短寿命的π介子(pion)的粒子的簇射(shower)，而所述π介子衰变并产生较长寿命的μ介子。μ介子主要通过库仑力与物质相互作用，而不具有核相互作用，并与电子相比更不容易辐射。它们只通过电磁相互作用缓慢地丢失能量。因此，许多μ介子作为高穿透性的带电辐射到达地球表面。在海平面的μ介子通量为大约每分钟每平方厘米一个μ介子。

当μ介子运动穿过材料的时候，亚原子粒子的电荷的库仑散射干扰它的轨迹(trajectory)。虽然总偏离量取决于几种材料性质，但是主导参数为原子核的原子序数Z和材料密度。

需要提供一种改进的方法和系统，用于从穿过体积(volume)的μ介子或其他带电粒子重建感兴趣的体积。

发明内容

以下发明概述被提供来帮助对涉及对探测诸如像μ介子那样的带电粒子的粒子的技术、装置和系统的一些技术特征，以及对带电粒子层析数据的物体体积散射密度剖面(object volume scattering density profile)的统计重建的理解，而不是为了全面的描述。对本发明各方面全面了解可以通过将整个说明书、权利要求、附图和摘要作为一个整体来获得。

本发明的上述方面以及其他目标和优势现在能够像这里所描述的那样达成。

根据一个方面，描述一种探测系统用于经由穿过物体体积的带电粒子来探测物体体积。所述系统包括：第一组位置敏感探测器，其位于物体体积的第一侧，来测量向物体体积入射的带电粒子的位置和角度；第二组位置敏感探测器，其位于与第一侧相对的物体体积的第二侧，来测量离开物体体积出射的带电粒子的位置和角度；和信号处理单元，其接收来自第一组位置敏感探测器的测量信号和来自第二组位置敏感探测器的测量信号的数据。信号处理单元处理所接收的数据，以产生物体体积内体积散射密度分布的统计重建。

信号处理单元能够被配置来：(a)获得相应于穿过物体体积的带电粒子的散射角度和估计动量(estimated momenta)的带电粒子层析数据；(b)基于统计多重散射模型(multiple scattering model)，来提供带电粒子散射密度的概率分布；(c)使用期望最大化(expectation maximization)(ML/EM)算法，来确定物体体积散射密度的实质最大似然性估计(substantially maximumlikelihood estimate)；以及(d)基于实质最大似然性估计，来输出重建的物体体积散射密度。

根据另一个方面，一种用于在从物体体积中获得的带电粒子层析数据中探测物体体积的方法包含：(a)获得相应于穿过物体体积的带电粒子的散射角度和估计动量的预定的带电粒子层析数据；(b)基于统计多重散射模型，提供带电粒子散射的概率分布；(c)使用期望最大化(ML/EM)算法，确定物体体积散射密度的实质最大似然性估计；(d)输出所重建的物体体积散射密度；以及，如果必要的话，(e)基于所重建的物体体积散射密度来做出决定。

本方法允许用户从所重建的体积散射密度剖面中识别占据(occupying)感兴趣的体积的物体的存在和/或类型。各种应用包括用于各种国土安全检查应用的宇宙射线μ介子层析，在所述应用中，车辆或货物能够通过μ介子跟踪器来扫描。

带电粒子层析数据能够包含从由宇宙射线或一些其他源生成的、诸如μ介子的带电粒子所采集的层析数据。

基于所重建的物体体积散射密度来做出决定(making a decision)能够包含：基于所重建的物体体积散射密度，来对物体体积中的目标物体的(1)存在和(2)类型中的至少一个做出决定。

提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布能够包含：(g)基于均质(homogenous)物体的预先定义的散射密度来获得带电粒子的2D概率分布；(h)基于2D概率分布来获得带电粒子的3D概率分布；(i)获得穿过经由基函数特性化的(characterized via basis functions)非均质(non-homogenous)物体体积的多个带电粒子的散射的概率分布；和(j)基于多重散射的定义和带电粒子的散射和动量测量，来确定所述多重散射(multiple scattering)的概率分布。

基于均质物体的预先定义的散射密度来获得带电粒子的2D概率分布能够包含：(k)确定材料的散射密度作为穿过材料的单位深度(unit depth)的带电粒子的期望均方散射；(l)基于高斯模型来近似(approximating)带电粒子的散射分布；以及(m)基于相关联(correlated)2D高斯分布来近似带电粒子射线的散射和位移。

基于2D概率分布来获得带电粒子的3D概率分布能够包含：添加坐标并定义三维路径长度；计算3D位移；以及定义3D协方差矩阵。

获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布能够包含：建立代表离散散射密度估计的基函数的3D网格(3Dgrid)；将每个单独的μ介子的散射/位移的协方差矩阵确定为射线路径和散射密度估计的函数；以及将多个带电粒子的概率分布确定为单独的带电粒子的概率的乘积。

使用期望最大化(ML/EM)算法来确定物体体积散射密度的实质最大似然性估计能够包含：采集每一个带电粒子的散射和动量的测量；估计每一带电粒子与统计散射模型的每一个基函数的相互作用的几何结构(geometry ofinteraction)；对每一带电粒子基函数对，确定权重矩阵：W_ij；以猜测(guess)来初始化每一基函数中的散射密度；以及迭代地求解(solving)物体体积容量(contents)的近似最大似然性解(solution)；其中在迭代的预定数目处、或当所述解变得小于预定的容限值时，停止迭代过程。

使用期望最大化(ML/EM)算法来确定物体体积的散射密度的实质最大似然性估计能够包含：采集每一带电粒子i＝1到M的散射和动量的测量(Δθ_x，Δθ_y，Δx，Δ_y，p_r ²)_i；估计每一μ介子与每一体素(voxel)j＝1到N的相互作用的几何结构：(L，T)_ij；对每一带电粒子体素对，将权重矩阵Wij计算为 $W_{\mathrm{ij}}\≡\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{T}_{\mathrm{ij}}^2\end{array}\right],$ 在每一体素中初始化散射密度的猜测λ_j，old；以及使用停止判据过程(a stopping criteria process)来对全部体素设置λ_j，old＝λ_j，new。

期望最大化(ML/EM)算法能够包括均值更新规则或中值更新规则。

根据再一个方面，用于在从物体体积获得的带电粒子层析数据中探测物体体积的计算机实施的方法包含：(a)获得相应于穿过物体体积的带电粒子的散射角度和估计动量的带电粒子层析数据；(b)提供用于期望最大化(ML/EM)算法的带电粒子散射密度的概率分布，而所述概率分布基于统计多重散射模型；(c)使用期望最大化(ML/EM)算法来确定物体体积散射密度的实质最大似然性估计；以及(d)输出所重建的物体体积散射密度。能够基于所重建的物体散射体积密度来做出决定。

根据再一个方面，计算机程序产品包含计算机可用的数据载体存储指令，而当由计算机执行时，所述指令导致计算机执行带电粒子层析数据的物体体积密度剖面的统计重建的方法，所述方法包含：(a)获得相应于穿过物体体积的带电粒子的散射角度和估计动量的预定的带电粒子层析数据；(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布，而所述概率分布基于统计多重散射模型；(c)使用期望最大化(ML/EM)算法来确定物体体积密度的实质最大似然性估计；以及(d)输出所重建的物体体积散射密度。

附图说明

在附图中，贯穿分离的视图的相同参考标号指称相同或功能相似的元件，并且附图并入说明书中、形成说明书的一部分，所述附图进一步说明了本发明，并与本发明的详细描述一起，用来解释本发明的原理。

图1说明了μ介子层析概念的一个例子；

图2说明了被用于确定库仑散射的散射和位移的二维投影；

图3说明了被用于调整3D散射模型的散射和位移的二维投影的参数；

图4说明了穿过多层材料的散射；

图5说明了为了在图3中所示的投影的路径长度计算而使用最接近点(point of closest approach)；

图6说明了用于μ介子层析的体积散射密度剖面的统计重建的自动化系统；

图7说明了模拟的物体的透视图；

图8说明了模拟的物体的俯视图；

图9说明了具有假设两条线在PoCA点处连接所估计的路径长度的高斯散射模拟的重建；

图10说明了使用有非高斯的尾部(tails)的散射的、模拟数据的重建；

图11说明了经由中值方法使用有非高斯的尾部的散射的、模拟数据的重建；

图12说明了在模拟的客运厢车(passenger van)中的主要物体；

图13说明了经由均值方法的客运厢车的一分钟的模拟μ介子曝光的重建；

图14说明了经由中值方法的客运厢车的场景的重建；

图15说明了根据一个实施例的流程图，该流程图一般地概述了用于带电粒子层析的体积散射密度剖面的统计重建的方法；

图15A说明了根据一个实施例的流程图，该流程图概述了使用多重统计散射模拟来估计穿过物体体积的μ介子的散射的概率分布的过程的例子；

图15B说明了根据一个实施例的流程图，该流程图概述了基于物体的预先定义的散射密度来估计单个μ介子在2D中散射的期望概率分布的过程的例子；

图15C说明了根据一个实施例的流程图，该流程图概述了将统计模型扩展到3D的过程的例子；

图15D说明了根据一个实施例的流程图，该流程图概述了确定穿过非均质材料的多个μ介子的散射和位移的概率分布的过程的例子；以及

图15E说明了根据一个实施例的流程图，该流程图概述了使用期望最大化算法来最大化物体体积的估计密度剖面的似然性的过程的例子。

具体实施方式

在这些非限制的例子中所讨论的特定值和配置可以变化，并且仅仅是为了说明本发明的至少一个实施例而提出，而不是为了限制本发明的范围。

在这个申请中所描述的技术特征能够用来构建各种粒子探测系统。例如，用于探测作为带电粒子的μ介子的粒子探测系统可以包括：物体存放区(object holding area)，用于放置将被检查的物体；第一组位置敏感μ介子探测器，位于物体存放区的第一侧，来测量向物体存放区入射的μ介子的位置和角度；第二组位置敏感μ介子探测器，位于与第一侧相对的、物体存放区的第二侧，来测量离开物体存放区出射的μ介子的位置和角度；以及信号处理单元，其可以包括，例如，微处理器，来接收来自第一组位置敏感μ介子探测器的进入μ介子的测量信号和来自第二组位置敏感μ介子探测器的出射的μ介子的测量信号的数据。作为一个例子，第一组和第二组粒子探测器中的每一个能够被实现为包括漂移管(drift tubes)，所述漂移管被安排以允许在第一方向中的至少三个带电粒子位置测量和在不同于第一方向的第二方向中的至少三个带电粒子位置测量。信号处理单元被配置以基于所测量的μ介子进入和出射的位置和角度，来分析由μ介子在物体存放区之内的材料中的散射所导致的μ介子的散射行为，以获得在物体存放区之内的散射中心(scattering centers)的层析剖面或空间分布。所获得的散射中心的层析剖面或空间分布能够被用来揭示在物体存放区中的一个或更多物体的存在或不存在，所述物体诸如具有高原子序数的材料，包括核材料或核装置。每一位置敏感μ介子探测器可以以各种配置来实现，包括漂移单元(drift cells)，如充满能够被μ介子电离(ionize)的气体的漂移管。这样的系统可以用来利用自然的宇宙射线μ介子作为用于探测在物体存放区中的一个或更多物体的μ介子的源。

在一个实施方式中，根据说明的实施例的用于带电粒子层析的体积散射密度剖面的统计重建的方法和系统提供了一种方案，其中基于运动穿过物体的宇宙射线带电粒子的散射来重建物体的图像或模型。

与组成更普通的物体的材料(诸如水、塑料、铝和钢)相比，特殊核材料(SNM)和可作为良好的伽马射线屏蔽的材料(诸如铅和钨)会更加强烈地影响轨迹。对于宇宙射线带电粒子，特别是宇宙射线μ介子，每一μ介子都携带有关于其所穿透的物体的信息，并且通过测量多个μ介子的散射，就能够探查这些物体的性质。具体来讲，人们能够探测在更典型的低Z(low-Z)和中等Z(medium-Z)物质中的高Z(high-Z)物体。

为了解释根据说明的实施例的用于带电粒子层析的体积密度剖面的统计重建的各种技术特征，将首先参考μ介子层析概念，在图1中说明了μ介子层析概念的例子。

位置敏感探测器组10被配置在将被成像的(imaged)物体体积11之上和之下，以提供全部进入的和出射的带电粒子轨道(track)12的位置和角度(由带箭头的实线示出)。被安排在将被成像的体积之上的两组或更多组位置敏感探测器10提供进入的带电粒子轨道的位置和角度。这些探测器在两个正交的或非正交的坐标中测量带电粒子的位置。另一组位置敏感探测器10记录出射的带电粒子的位置和角度。侧面探测器(未示出)可以用来探测更加水平朝向(orientate)的带电粒子的轨道。从一致的(coincident)进入的和出射的测量中计算每一带电粒子轨道的散射角度。从在探测器自身中、或在被放置在两组位置敏感探测器平面之间的已知性质的散射体(scatterers)的层中发生的轻微的散射中，估计带电粒子动量。

位置敏感带电粒子探测器的一个例子为被充满工作气体(operating gas)的漂移管。漂移管可以是圆柱形管，并被充满诸如氩-异丁烷(Argon-Isobutane)的探测器气体，来探测诸如μ介子的宇宙射线的带电粒子。将大约+2-3千伏的正的HV(High Voltage，高电压)施加到沿着圆柱形管的长度延伸的中心阳极线且以管的外表面接地，以使得存在高压静电场。当带电粒子与气体原子相互作用时，沿着穿过管的弦长(chord)的直线从所述原子中释放出许多电子。静电场使得电子“串”(string)向带正电的阳极线漂移，阳极线由数据采集电子仪器的TDCS(time-to-digital converter，时间-数字转换器)电子地读出。每组探测器可以是多个漂移管，所述多个漂移管被安排以允许在第一方向中的至少三个带电粒子位置测量，和在不同于第一方向、并可以与第一方向正交的第二方向中的至少三个带电粒子位置测量。

在图1内的系统中提供了信号处理单元，例如，计算机，来接收由在物体体积之上的探测器测量的进入的μ介子、以及由在物体体积之下的探测器测量的出射的μ介子的信号的数据。所述信号处理单元被配置来基于测量的μ介子的进入的和出射的位置和角度来分析由体积中的散射所导致的μ介子的散射行为，以获得在体积之内的散射中心的层析剖面或空间分布。所获得的在体积之内的散射中心的层析剖面或空间分布能够揭示体积中物体的存在或不存在。在某些实施方式中，可以在体积的侧面上实现附加的漂移管探测器，以形成盒或四面结构(four sided structure)，而包裹、车辆或货物集装箱能够进入其中，以便用系统进行扫描。因而，能够使用宇宙射线μ介子的多重散射来在正常货物的背景中选择性地探测高z材料(z-material)。有利地，本技术为无源的，不在背景上释放任何辐射剂量，并且对高z致密物体为选择性的。可以在位于与探测器10相同的位置的本地计算机(on-premisecomputer)中实现信号处理单元的层析处理部分。可替换地，可以在远程计算机中实现信号处理单元的层析处理部分，而该远程计算机连接到诸如专用网络或公用网络(诸如互联网)的计算机网络上。

在图1的说明性实施例中，带电粒子是宇宙射线μ介子或其他宇宙射线带电粒子，而位置敏感探测器10为充满用于感测(sensing)带电粒子的工作气体的漂移单元。例如，可以用漂移管来实现漂移单元，所述漂移管具有沿着每一漂移管的纵向延伸的中心阳极线。然而，能够使用不同于漂移单元的位置敏感传感器来探测不同于μ介子的带电粒子。此外，能够由不同于宇宙射线的源生成带电粒子。例如，μ介子能够从加速器中作为低强度射束(lowintensity beam)生成。

穿透致密物体(黑色轨道)的μ介子比穿过空气(灰色轨道)的μ介子更加强烈(stronger)地散射。从多个轨道测量，能够重建全部材料的物体几何结构和电子密度。穿过体积的μ介子以取决于μ介子所穿过的材料的方式散射。

由图1中系统的处理单元检查的体积(例如，包裹、集装箱或车辆)中对宇宙射线μ介子的测量的处理可以包括：重建μ介子穿过体积的轨迹，基于来自在体积的每一侧上的探测器的信号测量进入的μ介子的动量，并且确定体积的散射密度的空间分布。这些和其他处理结果可以用来构建层析剖面和测量体积的各种性质，如探测目标物体。

例如，穿过具有一组漂移管的探测器10的带电粒子的轨迹的重建能够包括：(a)获得代表由带电粒子撞击的漂移单元的标识符和相应的撞击时间的撞击信号；(b)将被识别为与穿过所述探测器的特定带电粒子的轨道相关联的时间吻合的(in-time)多个漂移单元撞击编组(grouping)；(c)初始地估计所述特定带电粒子撞击漂移单元的时刻的时间零点值(time zero value)；(d)基于时间零点值的估计、漂移时间转换数据和撞击的时间来确定漂移半径；(e)将线性的轨道拟合(fitting)到相应于特定时间零点值的漂移半径；以及(f)搜索和选择与为特定带电粒子执行的最好轨道拟合相关联的时间零点值，并计算时间零点和轨道参数中的误差。这种基于时间零点拟合的轨道重建提供了穿过带电粒子探测器的带电粒子的重建的线性轨迹，而不必使用快速探测器(如有闪烁器桨(scintillator paddles)的光电倍增管)或一些其他的快速探测器，所述其他的快速探测器将μ介子对本装置的穿过探测到最近的几个纳秒，以便提供时间零点。

又例如，由于基于来自图1中的探测器10的信号测量进入的或出射的μ介子的动量的处理能够包括：(a)配置多个位置敏感探测器来散射穿过那里的带电粒子；(b)测量位置敏感探测器中带电粒子的散射，其中测量散射包含获得散射的带电粒子的至少三个位置的测量；(c)从位置的测量来确定带电粒子的至少一个轨迹；以及(d)从所述至少一个轨迹确定带电粒子的至少一个动量测量。本技术可以用来基于带电粒子的轨迹确定带电粒子的动量，而带电粒子的轨迹由位置敏感探测器自身中带电粒子的散射确定，而没有使用探测器中附加的金属板。

下面提供对带电粒子层析数据的物体体积散射密度剖面的统计重建的示例系统和方法的细节。

在图6内的框图中说明了根据一个实施例的、带电粒子层析的体积散射密度剖面的统计重建的自动化系统的例子。自动化系统50具有被适配和安排以接收带电粒子层析数据54的控制器51。例如，带电粒子层析数据可以是μ介子层析数据，所述μ介子层析数据使用图1的带电粒子探测器1，或可替换地，使用任何具有被配置以使能对穿过体积的带电粒子的跟踪的位置敏感探测器的其他带电粒子探测器，而从μ介子的测量中确定。结果，μ介子或其他带电粒子层析数据可以被用来提取或确定穿过物体体积的μ介子或其他带电粒子的散射角度和估计动量。

自动化系统50包括被存储在控制器上的统计重建器模块56。重建器模块56负责统计重建μ介子或其他带电粒子的层析的体积散射密度剖面。模块56可以实现为软件模块或硬件模块。

在图6的自动化系统50的说明性实施例中，使用一个或更多诸如计算机(PC)的基于可操作地链接的计算机处理器单元(CPU)的系统、或诸如基于数字信号处理器的系统的其他基于微处理器的系统来形成控制器51。控制器可以是单一的标准计算机，但是为了达到实时结果，控制器典型地包括在数目上足够提供达到实时结果所需处理能力的、并行处理计算机的机群(farm)(未示出)。例如，控制器能够包括，比方说，20个CPU。μ介子探测器的扫描体积越大，并且所期望的分辨率越好，则需要越大的处理计算机机群(processing computer farm)。

操作系统在控制器51上运行，并且所述操作系统可以为商业上可获得的或开放源码的操作系统，包括但不限于来自Apple、Windows、Linux或Unix或其他将来可能开发的操作系统。由于操作系统的指令以及应用或程序被存储在诸如硬盘驱动器的存储设备中。而且，在自动化系统50中，跟踪重建器模块56是存储指令的计算机可用的数据载体形式的软件，当由控制器执行时，所述指令使得控制器执行根据说明的实施例的、用于带电粒子层析的体积散射密度剖面的统计重建的方法。模块能够如图6中所示安装在控制器本地，或经由耦接到控制器的网络从遥远位置运行。本领域技术人员应当理解有多种实现这种模块的方式。

如果需要，自动化系统50还包括显示器58，其可操作地耦接到控制器51，由于向用户显示由系统重建的物体密度剖面的图像或数据。如果需要，用户接口(未示出)可以可操作地连接到处理系统，以允许操作人员来操作处理系统。

本领域技术人员将理解，图6的说明仅仅表示自动化系统50的实施例的一个例子，并且实施例不限于此。例如，可以将重建器模块功能中的一些或全部实现为诸如模拟或数字电路的硬件，而不使用微处理器。

参考图15，其说明了根据一个实施例的、一般地概述用于带电粒子层析的体积密度剖面的统计重建的方法的流程图。如处理步骤101所示，方法100通过获得相应于穿过物体体积的带电粒子轨道位置、散射角度和估计动量的预定的带电粒子层析数据来初始化。例如，能够从图1的探测器获得预定的带电粒子层析数据。之后，如处理步骤102所示，基于多重统计散射模型，提供由散射密度的空间分布(将在以下定义)所代表的、多个带电粒子穿过物体体积的散射的概率分布。然后，如处理步骤103所示，使用期望最大化算法来确定物体体积散射密度剖面的最大似然性估计。然后，将所重建的体积散射密度剖面输出以做出决定(处理步骤104)。如处理步骤105所示，做出决定的过程为可选的，并可以是用于识别占据体积的物体的存在和/或类型的过程。做出决定的过程能够涉及代表物体体积的重建的密度剖面的图像的人工解读和/或通过附加的算法的自动化决定。

实施例的方法和自动化系统允许基于由许多带电粒子提供的数据执行感兴趣的体积的离散层析重建。迭代的期望最大化(EM)算法的实例被用于寻找物体的密度剖面的最大似然性估计。实施例的方法和系统允许用户从重建的体积密度剖面中识别占据感兴趣的体积的物体的存在和/或类型。各种应用包括用于各种国土安全检查应用的宇宙射线μ介子层析，其中能够由μ介子跟踪器(tracker)扫描车辆或货物。使用说明性实施例的方法和自动化系统，作为结果的μ介子层析数据可用于来重建和显示车辆或货物的密度剖面，以允许识别任何构成威胁的物体。

在将最大似然性用于医学图像重建，特别地，对于PET和SPECT重建的同时，几个重要的不同妨碍为那些应用开发的标准方法的使用。首先，测量信号-散射角度-是随机的，并且具有等于零的均值(mean)和由穿透的材料的性质定义的标准偏差(deviation)。其次，宇宙射线μ介子并非来自规定好的离散的方向或角度，而是具有在天顶角度周围并几乎延伸到地平线角度的宽广、连续的角度分布。最后，μ介子的轨迹不是直的；轨迹的弯曲使得我们能够找到强烈散射物体的粗略位置。EM算法是灵活和计算有效率的，并且能够说明它对复杂的几何结构的应用。

现在将根据一个实施例来描述处理步骤102至104，其中数据是从测量的穿过体积的μ介子的、图1的探测器获得的宇宙射线μ介子层析数据。

在图15A的流程图中概述了根据一个实施例的、使用多重统计散射模型来提供μ介子穿过物体体积的散射的估计概率分布的过程(处理步骤102)。如处理步骤110至113所示，所述过程具有四个主要的组成部分。首先，估计基于均质物体的预先定义的散射密度的、单个μ介子的2D概率分布(处理步骤110)。然后，将2D分布扩展到3D(处理步骤111)。接下来，在处理步骤112中，使用体素基函数来表达非均质物体体积，并且表达对给定体素化(voxelized)散射密度的、多个μ介子的散射的概率分布。最后，将概率分布表达式扩展到有限精度的μ介子散射和动量测量(处理步骤113)。

使用多重散射统计模型来实现处理步骤110至113，而将首先参考单层均质材料中的散射、然后参考非均质材料中的散射来描述所述模型。

如在图2中所说明的那样，穿过材料的宇宙射线μ介子经历多重库仑散射(multiple Coulomb scattering)，图2说明被用于描述多重库仑散射的散射和位移的二维投影。在这个以及其他附图中，为了说明的目的，极大地夸大散射的幅度。与入射的μ介子的朝向和位置对比，可以通过散射角度和位移来特性化出射的μ介子轨道。典型的散射角度为几十毫弧度(milliradian)(1毫弧度≈0.06度)，并且比几度更大的散射角度为非常罕见的。可以将散射角度的中央98％的分布近似为零均值高斯型(zero-mean Gaussian)。

$f_{\mathrm{\Δ\θ}}\left(\mathrm{\Δ\θ}\right)\≅\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{\θ}}}\exp (-\frac{{\mathrm{\Δ\θ}}^2}{{2\mathrm{\σ}}_{\mathrm{\θ}}^2}),$ 公式(1)

虽然真实的分布与高斯型相比具有更重或更大的尾部。可以按照材料的性质来近似地表达分布的宽度。如在S.Eidelman等人的、Phys.Lett，vol.B592，p.1，2004的“Review of particle physics”中所评论的那样，许多研究者已给出作为各种材料性质的函数的用于散射的经验表达式，通过参考将其公开内容并入这里。一个特别简单的形式为

${\mathrm{\σ}}_{\mathrm{\θ}}\≅\frac{15\mathrm{MeV}}{\mathrm{\βcp}}\sqrt{\frac{H}{L_{\mathrm{rad}}}}.$ 公式(2)

这里，p为以MeV/c计的粒子动量，H为材料深度，并且L_rad为材料的辐射长度，βc为速度(c为光速)，并且使用i、β＝1的近似。辐射长度随原子序数和材料密度的增加而减小。我们建立标称(nominal)μ介子动量，p₀，并将具有辐射长度L_rad的材料的散射密度定义为：

$\mathrm{\λ}\left(L_{\mathrm{rad}}\right)\≡{\left(\frac{15}{p_0}\right)}^2\frac{1}{L_{\mathrm{rad}}}.$ 公式(3)

因而，材料的散射密度λ代表穿过那个材料的单位深度的、具有标称的动量的μ介子的均方散射角度。例如，对于某些材料，以平方毫弧度每厘米计的散射密度的值为：对铝大约为3，对铁大约为14，而对铀大约为78。因而，穿过具有散射密度λ和深度H的材料的、具有动量p的μ介子的散射的方差为

${\mathrm{\σ}}_{\mathrm{\θ}}^2=\mathrm{\λH}{\left(\frac{p_0}{p}\right)}^2.$ 公式(4)

令

$p_r^2={(p_0/p)}^2,$ 公式(5)

得出

${\mathrm{\σ}}_{\mathrm{\θ}}^2=\mathrm{\λ}{\mathrm{Hp}}_r^2.$ 公式(6)

将位移Δx与散射角度Δθ相关联(correlated)。一并考虑，散射角度和位移提供了指出在较大体积中的局部散射贡献者(contributors)的位置的信息，如图1中路径的“弯折”(kinks)所指出的。可以将散射角度和位移的分布特性化为具有零均值的联合高斯型(jointly Gaussian)，并且

${\mathrm{\σ}}_{\mathrm{\Δx}}=\frac{H}{\sqrt{3}}{\mathrm{\σ}}_{\mathrm{\Δ\θ}},$ 公式(7)

${\mathrm{\ρ}}_{\mathrm{\Δ\θ\Δx}}=\frac{\sqrt{3}}{2}.$ 公式(8)

我们可以将协方差矩阵表达为

$\mathrm{\Σ}\≡\left[\begin{array}{cc}{\mathrm{\σ}}_{\mathrm{\Δ\θ}}^2& {\mathrm{\σ}}_{\mathrm{\Δ\θ\Δx}}\\ {\mathrm{\σ}}_{\mathrm{\Δ\θ\Δx}}& {\mathrm{\σ}}_{\mathrm{\Δx}}^2\end{array}\right]=\mathrm{\λ}\left[\begin{array}{cc}H& H^2/2\\ H^2/2& H^3/3\end{array}\right]p_r^2.$ 公式(9)

令

$A\≡\left[\begin{array}{cc}H& H^2/2\\ H^2/2& H^3/3\end{array}\right],$ 公式(10)

得出

$\mathrm{\Σ}=\mathrm{\λA}{p}_r^2.$ 公式(11)

注意到前面提到的，能够如在图15B的流程图中所概述的那样来描述根据一个实施例来获得单个μ介子散射的、2D的散射的概率分布(处理步骤110)。如处理步骤150所示，按照公式(3)，将材料的散射密度定义为p₀＝3GeV/c的μ介子穿过所述材料的单位深度的期望均方散射。然后，如处理步骤151、公式(1，5，6)中所示，对RMS散射进行高斯近似。最后，如处理步骤152所示，经由公式(10，11)所总结的，由零均值相关联的2D高斯分布来近似射线的散射和位移分布。

在三维中，考虑正交于x的y坐标来特性化散射，并且参考散射角度Δθ_x和Δθ_y、以及位移Δx和Δy。到x和y平面的偏离量为独立的和同样分布的(见Eidelman等)。上面的发展(development)基于朝向正交于入射的μ介子的方向的坐标系统。在3-D模型中，我们必须考虑3-D路径长度，并将位移测量投影到正交于入射的μ介子路径的平面。在说明用于调整3-D散射模型的参数的图3中，说明了以与垂直方向成θ_x，0的投影角度入射的μ介子。

为了帮助理解这个3-D散射，想象从页面朝外的正交y坐标中相关联的投影角度θ_y，0是有用的。穿过所述层到所投影的(非散射的)点(x_p，y_p)的μ介子路径的直线延伸(即，3-D路径长度)为

$L=H\sqrt{1+{\tan }^2{\mathrm{\θ}}_{x,0}+{\tan }^2{\mathrm{\θ}}_{y,0}}\≡{\mathrm{HL}}_{\mathrm{xy}}.$ 公式(12)

将出射的μ介子的x位置和角度定义为(x₁，θ_x，1)，然后令

Δθ_x＝θ_x，1-θ_x，0.          公式(13)

虽然测量的x位移可以计算为x_m＝x₁-x_p，但是我们必须将这个测量旋转到正交于射线路径的平面，并为3-D路径长度进行调整。将位移定义为

${\mathrm{\Δ}}_x=(x_1-x_p)\cos \left({\mathrm{\θ}}_{x,0}\right)L_{\mathrm{xy}}\frac{\cos (\mathrm{\Δ}{\mathrm{\θ}}_x+{\mathrm{\θ}}_{x,0})}{\cos \left(\mathrm{\Δ}{\mathrm{\θ}}_x\right)},$ 公式(14)

其中，中间的两项考虑3-D路径长度，并且最后一项将测量投影到适当的朝向。

最后，将协方差权重重新定义为

$A\≡\left[\begin{array}{cc}L& L^2/2\\ L^2/2& L^3/3\end{array}\right].$ 公式(15)

然后，需要对散射和位移以相似的方式继续进行，并且公式(11)定义对这两者以及坐标散射的协方差矩阵。在两个正交的、水平坐标中独立进行散射测量。为了简化记号，我们发展(develop)对仅仅一个坐标的分析。稍后将讨论组合来自两个坐标的信息。我们必须注意，该模型对“小的”散射角度和位移是有效的(valid)。模型推导中忽略的二阶项对于大角度散射可能变得显著。

注意到前面提到的，根据一个实施例，如在图15C的流程中所概述的那样，获得延伸到3D的统计模型。首先，添加y坐标，并且定义三维路径长度(处理步骤160、公式(12))。接下来，在处理步骤161中，按照公式(13-14)来计算3D位移。最后，按照公式(15)表达3D协方差矩阵(处理步骤162)。

对材料的非均质体积，为了重建的目的，按照利用系数{v₁，...，v_j，v_N}的3-D基函数{φ₁，...，φ_j，φ_N}的线性组合来代表密度剖面，即，

$\mathrm{\λ}(x,y,z)=\underset{j}{\mathrm{\Σ}}{v}_j{\mathrm{\φ}}_j(x,y,z).$ 公式(16)

虽然对于基函数存在许多选择，但是此处我们的注意力在于矩形的3D体素。λ_j被用来表示第j个基函数的系数，即，第j个体素中的散射密度。考虑图4，示出三个层(或体素)，以射线穿过所述栈(stack)，传达(delivering)观测的信息Δθ和Δx。第j个体素中“隐藏的”散射和位移分别表示为Δθ_j和Δx_j。再次，在图中散射的量被夸大。我们可以通过以下表达式将所观测的与隐藏的数据联系起来

Δθ＝Δθ₁+Δθ₂+Δθ₃，              公式(17)

Δx＝Δx₁+L₂tan(Δθ₁)+Δx₂+

L₃tan(Δθ₁+Δθ₂)+Δx₃

≈Δx₁+Δx₂+Δx₃+T₁Δθ₁+T₂Δθ₂.      公式(18)

此处，我们在第二个公式中依赖小角度散射的假设，并将T_j定义为从第j个体素的出口点到从重建体积的出口点的3-D射线路径长度。更一般地讲，对于穿过一组体素的射线，

$\mathrm{\Δ\θ}=\underset{j\∈\ℵ}{\mathrm{\Σ}}{\mathrm{\Δ\θ}}_j,$ 公式(19)

$\mathrm{\Δx}=\underset{j\∈\ℵ}{\mathrm{\Σ}}({\mathrm{\Δx}}_j+T_j\mathrm{\Δ}{\mathrm{\θ}}_j).$ 公式(20)

最后，我们可以表达对第i个射线的、合计的散射/位移的协方差，通过首先注意到，对于第j个体素

${\mathrm{\Σ}}_{\mathrm{ij}}={\mathrm{\λ}}_j{A}_{\mathrm{ij}}{p}_{r,i}^2,$ 公式(21)

其中

$A_{\mathrm{ij}}\≡\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2\\ L_{\mathrm{ij}}^2/2& L_{\mathrm{ij}}^3/3\end{array}\right]$ 公式(22)

并且，L_ij为第i个射线穿过第j个体素的路径长度，对未被射线“撞击”的体素定义为零。组合公式(19)至(22)，我们可以写出

${\mathrm{\Σ}}_i=p_{r,i}^2\underset{j\≤N}{\mathrm{\Σ}}{\mathrm{\λ}}_j{W}_{\mathrm{ij}}.$ 公式(23)

此处，N为体素总数，并且我们基于对要素(elements)的简单但冗长的计算，将权重矩阵定义为

$W_{\mathrm{ij}}\≡\left[\begin{array}{cc}{L}_{\mathrm{ij}}& L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}\\ L_{\mathrm{ij}}^2/2+L_{\mathrm{ij}}{T}_{\mathrm{ij}}& L_{\mathrm{ij}}^3/3+L_{\mathrm{ij}}^2{T}_{\mathrm{ij}}+L_{\mathrm{ij}}{T}_{\mathrm{ij}}^2\end{array}\right],$ 公式(24)

对未知μ介子路径进行一些假设，以便估计穿过体素的射线路径长度。参考图5，所述近似以计算进入的和出射的轨道的最接近点(point of closestapproach，PoCA)(x_ca，y_ca)开始。然后，将PoCA的入口(entry)到出口点被连接，以便估计体素路径长度。

最后，定义数据矢量

$D_i\≡\left[\begin{array}{c}\mathrm{\Δ}{\mathrm{\θ}}_i\\ {\mathrm{\Δx}}_i\end{array}\right]$ 公式(25)

并且令D表示来自M个μ介子的全部测量。我们将给定密度剖面λ的数据似然性写为

$P\left(D\right|\mathrm{\λ})=\underset{i\≤M}{\mathrm{\Π}}P\left(D_i\right|\mathrm{\λ})$ 公式(26)

其中利用因子

$P\left(D_i\right|\mathrm{\λ})=\frac{1}{2\mathrm{\π}{\left|{\mathrm{\Σ}}_i\right|}^{1/2}}\exp (-\frac{1}{2}{D}_i^T{\mathrm{\Σ}}_i^{-1}{D}_i).$ 公式(27)

注意到前面提到的，如图15D中所概述的那样，根据一个实施例，能够获得多个μ介子穿过非均质材料的散射和位移的概率分布。首先，建立体素的3D网格(或者其他基函数)(处理步骤170)。然后，在处理步骤171中，按照公式(23，24)计算每一μ介子的散射/位移的协方差矩阵。最后，如处理步骤172所示，按照公式(25-27)，给定射线路径和体素散射密度，计算全部μ介子的总体概率分布。

已描述多重散射的统计模型，现在将参考对实验效果的模型的延伸(处理步骤113)。真实的μ介子探测器呈现有限的位置分辨率。由从拟合到多个位置测量的轨道推导的角度和位置来特征化进入的和出射的μ介子的轨道。因而，测量误差传播到构成μ介子层析的数据集(dataset)的散射角度和位移测量。由RMS误差e_p特性化给定探测器的精度。对于探测器的特定布置，可以基于误差如何传播来定义误差矩阵

$E\≡\left[\begin{array}{cc}{e}_{\mathrm{\Δ\θ}}^2& e_{\mathrm{\Δ\θ\Δx}}\\ e_{\mathrm{\Δ\θ\Δx}}& e_{\mathrm{\Δx}}^2\end{array}\right]$ 公式(28)。

在迭代的重建方法中，这样的误差是相对容易处理的。在我们的情况下，我们可以通过补充公式(23)的协方差矩阵来考虑探测器误差

${\mathrm{\Σ}}_i=E+p_{r,i}^2\underset{j\≤N}{\mathrm{\Σ}}\left({\mathrm{\λ}}_j{W}_{\mathrm{ij}}\right).$ 公式(29)

以这种方式，减少噪声，否则所述噪声将由于探测器误差而出现在重建中。对探测器误差的更加精确的模型应该考虑动量依赖性，因为跟踪误差(tracking error)的一个源是探测器自身中的散射，并且所述散射随着粒子动量增加而减小。如果单独的(individual)μ介子动量的估计是可用的，那么能够对每一射线估计误差矩阵如从公式(2)所明显看出的那样，多重库仑散射的宽度取决于粒子动量。通过在公式(5)中引入因子p_r ²来考虑不同的μ介子动量。实践中，μ介子动量不是精确地已知的，但是可以从已知散射体中的散射测量(如宇宙射线μ介子的已知光谱)中估计单独μ介子的动量的估计。此处，假设我们有对每一μ介子的很好的估计

可以使用期望最大化算法来确定物体体积密度的最大似然性估计(方法100的处理步骤103)。EM算法依赖按照“完整”数据(即，观测数据加上隐藏数据)来表达“不完整”数据的似然性。在我们的应用中，观测数据D＝{D_i：1≤i≤M}是测量的散射。隐藏数据H＝{H_ij：1≤i≤M&1≤j≤N}为穿过第j个体素的第i个μ介子的散射角度和位移。在A.Dempster，N.Laird和D.Rubin的“Maximum likelihood from incomplete data via the EM algorithm”(J.Roy.Statist.Soc.B，vol.39，pp.1-78，1977)中(通过引用将其公开内容并入这里)，按照下列的辅助函数来描述算法：

$Q_{\mathrm{DLR}}=E_{H|D,{\mathrm{\λ}}^{\left(n\right)}}[log\left(P(D,H|\mathrm{\λ})\right)].$ 公式(30)

对于H给定和参数矢量λ⁽ⁿ⁾的条件分布(conditional distribution)，给定参数矢量λ，这个函数是观测的和未观测的数据两者的对数似然性(loglikelihood)的期望值。由下列的两个步骤组成算法的每一迭代。

E步骤：估计或特性化P(H|D，λ⁽ⁿ⁾)，隐藏数据的条件概率。

M步骤：最大化辅助函数Q，该辅助函数Q是关于E步骤中特性化的分布的期望值。

在我们的情况下，由于隐藏数据唯一地确定所观测的数据，其中通过使用更简单的辅助函数

$Q(\mathrm{\λ};{\mathrm{\λ}}^{\left(n\right)})=E_{H|D,{\mathrm{\λ}}^{\left(n\right)}}\left[\log \left(P\left(H\right|\mathrm{\λ})\right)\right]$ 公式(31)

我们获得可以通过使用Q_DLR来获得的相同的序列估计从参数估计λ⁽ⁿ⁾中，由如下公式，所述算法的迭代产生新的估计λ⁽ⁿ⁺¹⁾

λ⁽ⁿ⁺¹⁾＝arg max_λ(Q(λ；λ⁽ⁿ⁾).        公式(32)

我们首先注意到单个μ介子穿过单个体素的散射的概率分布简单地从单层的统计模型得出。

$P\left(H_{\mathrm{ij}}\right|\mathrm{\λ})=\frac{1}{2\mathrm{\π}{\left|{\mathrm{\Σ}}_{\mathrm{ij}}\right|}^{1/2}}\exp (-\frac{1}{2}{H}_{\mathrm{ij}}^T{\mathrm{\Σ}}_{\mathrm{ij}}^{-1}{H}_{\mathrm{ij}}),$ 公式(33)

其中， ${\mathrm{\Σ}}_{\mathrm{ij}}={\mathrm{\λ}}_j{A}_{\mathrm{ij}}{p}_{r,i}^2,$ 定义在公式(21)中。由于在每一体素中的散射的无条件分布独立于在其他体素中的散射，因此合计的隐藏数据集的概率为每一要素的概率的乘积。从而，可以将对数似然性写为

$\log \left(P\left(H\right|\mathrm{\λ})\right)=\underset{j\≤N}{\mathrm{\Σ}}\underset{i:L_{\mathrm{ij}}\≠0}{\mathrm{\Σ}}(-\log {\mathrm{\λ}}_j-\frac{H_{\mathrm{ij}}^T{A_{\mathrm{ij}}}^{-1}{H}_{\mathrm{ij}}}{{2\mathrm{\λ}}_j{p}_{r,i}^2})+C,$ 公式(34)

其中C代表不包含λ的项。考虑条件期望，我们将Q函数写为

$Q(\mathrm{\λ};{\mathrm{\λ}}^{\left(n\right)})=C+\underset{j\≤N}{\mathrm{\Σ}}{Q}_j({\mathrm{\λ}}_j;{\mathrm{\λ}}_j^{\left(n\right)})$ 公式(35)

其中利用求和项(summands)

$Q_j({\mathrm{\λ}}_j;{\mathrm{\λ}}_j^{\left(n\right)})=-M_j\log {\mathrm{\λ}}_j-\frac{1}{2{\mathrm{\λ}}_j}\underset{i:L_{\mathrm{ij}}\≠0}{\mathrm{\Σ}}{S}_{\mathrm{ij}}^{\left(n\right)}.$ 公式(36)

这里，M_j是其L_ij≠0的射线的数目(即，撞击第j个体素的射线的数目)，并且S_ij ⁽ⁿ⁾定义为

$S_{\mathrm{ij}}^{\left(n\right)}\≡E_{H|D,{\mathrm{\λ}}^{\left(n\right)}}\left[p_{r,i}^{-2}{H}_{\mathrm{ij}}^T{A}_{\mathrm{ij}}^{-1}{H}_{\mathrm{ij}}\right].$ 公式(37)

将公式(36)相对于λ_j的导数(derivative)设置为零，我们找到以下最大化辅助函数(M步骤)的迭代公式

${\mathrm{\λ}}_j^{(n+1)}=\frac{2}{{2M}_j}\underset{i:L_{\mathrm{ij}}\≠0}{\mathrm{\Σ}}{S}_{\mathrm{ij}}^{\left(n\right)}.$ 公式(38)

S_ij的二次型形式保证λ⁽ⁿ⁺¹⁾的正性(positivity)。其保留以计算条件期望S_ij。令X表示随机变量H_ij|D_i。二次型X^TA^-1X的期望值为

$E\left[X^T{A}^{-1}X\right]=\mathrm{Tr}\left(A^{-1}{\mathrm{\Σ}}_X\right)+{\mathrm{\μ}}_X^T{A}^{-1}{\mathrm{\μ}}_X,$ 公式(39)

其中μ_X和∑_X分别为X的均值和协方差。由于D_i线性地依赖于H_ij，因此它们为联合高斯型。给定D_i的H_ij的条件分布也是高斯型，来自多变量(multivariate)分布理论的结果。使用所述理论以及H_ij与D_i每个都具有零均值的事实，我们找到

${\mathrm{\μ}}_X={\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\Σ}}_{D_i}^{-1}{D}_i,$ 公式(40)

${\mathrm{\Σ}}_X={\mathrm{\Σ}}_{H_{\mathrm{ij}}}-{\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\Σ}}_{D_i}^{-1}{\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}.}$ 公式(41)

这里，由公式(29)给定观测数据的协方差并且可以经由公式(21)表达隐藏数据要素的协方差我们能够执行简单的(虽然为冗长的)矩阵计算以示出(而非明确地写出观测的隐藏数据的协方差)

${\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}}{A}_{\mathrm{ij}}^{-1}{\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}}^T=W_{\mathrm{ij}}{\left(p_{r,i}^2{\mathrm{\λ}}_j\right)}^2.$ 公式(42)

将来自公式(39)至(42)的结果替代到公式(37)中，我们找到

$S_{\mathrm{ij}}^{\left(n\right)}=p_{r,i}^{-2}\mathrm{Tr}(A_{\mathrm{ij}}^{-1}{\mathrm{\Σ}}_{H_{\mathrm{ij}}}-A_{\mathrm{ij}}^{-1}{\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}}^T{\mathrm{\Σ}}_{D_i}^{-1}{\mathrm{\Σ}}_{D_i{H}_{\mathrm{ij}}})+p_{r,i}^{-2}{D}_i^T{\mathrm{\Σ}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\Σ}}_{D_i}^{-1}{D}_i{\left(p_{r,i}^2{\mathrm{\λ}}_j^{\left(n\right)}\right)}^2$

$=2{\mathrm{\λ}}_j^{\left(n\right)}+(D_i^T{\mathrm{\Σ}}_{D_i}^{-1}{W}_{\mathrm{ij}}{\mathrm{\Σ}}_{D_i}^{-1}{D}_i-\mathrm{Tr}\left({\mathrm{\Σ}}_{D_i}^{-1}{W}_{\mathrm{ij}}\right))p_{r,i}^2{\left({\mathrm{\λ}}_j^{\left(n\right)}\right)}^2,$ 公式(43)

其中，我们已在上一个步骤中使用Tr(AB)＝Tr(BA)。

最后，为了并入x和y坐标散射数据两者，我们简单地在更新的公式(38)中使用平均

$S_{\mathrm{ij}}^{\left(n\right)}=\frac{S_{\mathrm{ij},x}^{\left(n\right)}+S_{\mathrm{ij},y}^{\left(n\right)}}{2},$ 公式(44)

注意到的是公式(38)代表正在撞击体素的射线的均值。此后，将把该公式的使用称作为均值方法。下面将示出，这个更新公式的可替换形式在减少由于外露层(outlier)μ介子数据造成的噪声中是有用的。由所改变的更新公式来定义算法的中值方法：

${\mathrm{\λ}}_j^{(n+1)}=\frac{1}{2}\mathrm{media}{n}_{i:L_{\mathrm{ij}}\≠0}\left(S_{\mathrm{ij}}^{\left(n\right)}\right),$ 公式(45)

注意到前面提到的，在图15E中概述了使用期望最大化算法的、物体体积的估计密度剖面的最大化似然性的过程(方法100的处理步骤103)。如处理步骤180所示，对每一μ介子i＝1到M采集散射和动量的测量：(Δθ_x，Δθ_y，Δx，Δ_y，p_r ²)_i。估计每一μ介子与每一体素j＝1至N的相互作用的几何结构：(L，T)_ij(处理步骤181)。对于每一μ介子体素对，使用公式(24)来计算权重矩阵：W_ij(处理步骤182)。以猜测λ_j，old来初始化每一体素中的散射密度(处理步骤183)。如下面的那样，在处理步骤184中指示停止判据。对于每一μ介子，使用公式(29)并且取逆矩阵(inverse taken)来计算对每一μ介子体素对，使用公式(43、44)计算条件期望项：S_ij，使用公式(38)或公式(45)计算λ_j，old，并且对全部体素设定λ_j，old＝λ_j，new。

现在将参考数值例子，以便进一步说明方法100。模拟与图1中所示出的相似的设置。作为最初的有效性测试，使用被设计来严密地匹配多重统计散射模型的简单模拟。单个探测器平面(而不是如图中所示出的3个)有2×2平方米的大小，并且在顶部和底部的探测器之间的垂直距离为1.1米。这些探测器理想地记录μ介子的位置和角度。使用简化的μ介子光谱，其中μ介子的动量均匀分布于500至10000MeV/c。粒子在上部探测器平面、以均匀跨越(spanning)垂直的投影角度进入体积。

根据处理步骤110至113模拟μ介子的多重散射和位移。如图7和8中所显现的那样，将物体放置在体积的中心1.1×1.1×1.1立方米部分中。模拟三个10×10×10立方厘米的立方体的材料钨(W)、铁(Fe)和铝(Al)，分别具有散射密度71.5、14.2和2.8mrad²/cm。模拟假设相应于大约10分钟的曝光的400000个μ介子在上部探测器栈入射。所述μ介子中的大约160000个未碰到(missed)下部探测器平面，留下240000用于重建。假设理想地知道每一μ介子的动量，5×5×5立方厘米的体素大小被用于重建，并实现上面所描述的均值方法。模拟以被充满空气的体积开始，并将算法运行100次迭代(足够块特征(block features)的收敛)。结果出现在图9中。相应于三个物体中的每一个的、8个体素的重建值的平均，对于(W，Fe，Al)块分别为(74.0，14.7，2.7)。组成每一立方体的8个体素的分数散布(fractional spread)(rms/mean)为(12.6％，13.2％，12.1％)。给定在模拟和反演模型(inversionmodel)之间的匹配，这个结果使反演算法和实施有效。

重建表现为与物体场景相同。相应于三个物体中的每一个的、8个体素的重建值的平均，对于(W，Fe，Al)块分别为(74.0，14.7，2.7)。组成每一立方体的8个体素的分数散布(rms/mean)为(12.6％，13.2％，12.1％)。给定在模拟和反演模型之间的匹配，这个结果使反演算法和实施有效。

其次，使用GEANT4 Monte Carlo程序包来重新模拟相同的场景。能够在J.Allison的“Geant4 developments and applications”(IEEE Trans.Nucl.Sci.，vol.53，no.1，pp.270-278，Feb.2006)的出版物中找到GEANT4的细节，通过引用将其公开内容并入这里。GEANT4实现了对多重散射的更加完整的、精确的和有效的模型。这个模型包括散射分布的中心高斯部分的宽度的更加精准的计算、重尾(heavy tails)的实现以及μ介子在穿过材料时的能量损失的模拟。还使用μ介子事件发生器，其复制宇宙射线μ介子的海平面角度和动量分布。假设模拟中探测器为理想的，理想地知道每一μ介子的动量，不包括宇宙射线电子或轨道二级粒子。结果出现在图10中。相应于(W，Fe，Al)块的体素值的平均分别为(674.4，63.4，5.4)。

体素值太高，并且中间以及低Z区域中的几个的错误分类(missclassification)是显然的。通过将全部体素值大约除以4以产生中间Z体素的正确的平均体素值归一化(normalizing)重建，并不产生对高和低Z体素的正确的值，也不消除全部错误分类。这个效应的原因是一小部分μ介子以高斯模型未能很好描述的方式散射。要求散射的投影角度分布的中心98％很好地近似为高斯型。全部μ介子的大约2％散射到相对于此处所描述的统计模型很大的角度，即，远大于高斯分布出现的。由于散射角度的平方确定拟合，因此效应会是显著的。落在所述尾部的μ介子散射产生太大的散射密度估计。

而且，诸如在图1的仪器内的μ介子的衰变的其他过程或显著的探测器误差会被错误地记录为非常大的角度散射事件(虽然在我们的模拟中这些源并没有出现)。这能够在体积中的任何地方发生，并且倾向于产生具有不合理地大的散射密度的单一体素。这样的事件应该消除，因为它们给出SNM的假阳性(false-positive)指示。

为了使EM算法容忍非高斯型的数据，可以由公式(45)代替均值更新规则公式(38)，即，使用中值方法。

在图11中示出了使用中值方法的结果。对(W，Fe，Al)区域的体素平均分别为(79.2，14.2，2.1)，而分数散布分别为(21.5％，26.3％，23.2％)。很清楚，使用中值更新规则改进了反演算法的鲁棒性。

在图12中示出了重建密度剖面的更加实际的例子，其说明了客运厢车的详细的GEANT4模拟。去掉厢车车体的主要组成部分的说明出现在图12中。在说明的中心的红色的块代表钨的10×10×10立方厘米的实心的块，作为高Z威胁物体的代表物(proxy)。在这种情况下，我们使用位于场景的四个长边的、模拟的探测器平面，以便利用到更多水平朝向的μ介子。根据均值方法和中值方法，模拟一分钟的宇宙射线μ介子曝光，并且执行使用5×5×5立方厘米大小体素的数据的重建。图13和14分别示出使用均值EM方法和中值EM方法进行的重建的可视化(visualization)。非高斯数据的效果在这个场景的均值方法重建中是明显的，表现为散射在图像上的更暗的点。在中值方法的重建中，这些假象全部消失，并且厢车的更加致密的组成部分(引擎、电池、传动系统)示出为(低Z)或(中Z)，而威胁物体突出地更暗(高Z)。中值方法的使用产生了倾向于对抗由非高斯的散射分布和其他异常事件引起的假阳性具有鲁棒性的结果。

这里阐述的实施例和例子用来最好地解释本发明及其实际应用，并且由此使得本领域技术人员能够制作和利用本发明。然而，本领域技术人员应当认识到前述描述和例子仅仅是为了说明和示例的目的。

本发明的其他改变和修改(如重建过程的调整)对本领域技术人员将是清楚的，并且涵盖这样的改变和修改是所附权利要求书的意图。

所阐述的描述并不打算为穷举的(exhaustive)，或为了限制本发明的范围。按照上面的教导，许多修改和改变都是可能的，而不背离以下权利要求书的范围。可以预期的是本发明的使用会涉及具有不同的特性的组成部分。

Claims (17)

1.一种从带电粒子层析数据中探测物体体积的方法，所述带电粒子层析数据从所述物体体积中获得，所述方法包含： 

(a)获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的预定带电粒子层析数据； 

(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布，所述概率分布基于统计多重散射模型； 

(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实质最大似然性估计；以及 

(d)基于所述实质最大似然性估计输出重建的物体体积散射密度； 

其中，提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布包含： 

(g)基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布； 

(h)基于所述2D概率分布获得所述带电粒子的3D概率分布； 

(i)获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布；以及 

(j)基于多重散射的定义和所述带电粒子的散射和动量测量，获得多重散射的概率分布。 

2.如权利要求1所述的方法，包含： 

基于所述重建的物体体积散射密度，对所述物体体积中的目标物体的(1)存在和(2)类型中的至少一个做出决定。 

3.如权利要求1所述的方法，其中，基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布包含： 

(k)确定材料的散射密度作为穿过所述材料的单位深度的所述带电粒子的期望均方散射； 

(1)基于高斯模型来近似带电粒子散射分布；以及 

(m)基于相关联2D 

高斯分布来近似带电粒子射线散射和位移。 

4.如权利要求1所述的方法，其中，基于所述2D概率分布获得所述带电粒子的3D概率分布包含： 

添加坐标并定义三维路径长度； 

计算3D位移；以及 

定义3D协方差矩阵。 

5.如权利要求4所述的方法，其中，获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布包含： 

建立代表离散散射密度估计的基函数的3D网格； 

将每个单独的μ介子的散射/位移协方差矩阵确定为射线路径和散射密度估计的函数；以及 

将多个带电粒子的概率分布确定为单独的带电粒子概率的乘积。 

6.如权利要求1所述的方法，其中，使用所述期望最大化(ML/EM)算法确定所述物体体积散射密度的实质最大似然性估计包含： 

采集对每一带电粒子的散射和动量的测量； 

估计每一带电粒子与所述统计多重散射模型的每一个基函数的相互作用的几何结构； 

对每一带电粒子基函数对，确定权重矩阵：W_ij； 

利用猜测初始化每一基函数中的散射密度；以及 

迭代地求解物体体积容量的近似最大似然性解； 

其中，在预定数目的迭代处或当所述解变得小于预定容限值时，停止所述迭代过程。 

7.如权利要求1所述的方法，其中，所述期望最大化(ML/EM)算法包括均值更新规则或中值更新规则。 

8.如权利要求1所述的方法，其中，所述带电粒子层析数据包含宇宙射线μ介子层析数据。 

9.一种从带电粒子层析数据中探测物体体积的、计算机实施的方法，所述带电粒子层析数据从所述物体体积中获得，所述方法包含： 

(a)获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的带电粒子层析数据； 

(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布，所述概率分布基于统计多重散射模型； 

(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实质最大似然性估计；以及 

(d)输出重建的物体体积散射密度； 

其中，提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布包含： 

(g)基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布； 

(h)基于所述2D概率分布获得所述带电粒子的3D概率分布； 

(i)获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布；以及 

(j)基于多重散射的定义和所述带电粒子的散射和动量测量，获得多重散射的概率分布。 

10.如权利要求9所述的方法，进一步包含基于所述重建的物体体积散射密度来做出决定。 

11.如权利要求9所述的方法，其中，所述带电粒子层析数据包含宇宙射线带电粒子层析数据。 

12.如权利要求11所述的方法，其中，基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布包含： 

确定材料的散射密度作为穿过该材料的单位深度的具有标称动量p₀＝3GeV的宇宙射线带电粒子的期望均方散射； 

使用高斯模型来近似宇宙射线带电粒子散射分布；以及 

通过相关联2D高斯分布来近似宇宙射线散射和位移； 

其中，基于所述2D概率分布获得所述带电粒子的3D概率分布包含： 

添加坐标并定义三维路径长度； 

计算3D位移；以及 

定义3D协方差矩阵，并且 

其中，获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布包含： 

建立代表离散散射密度估计的基函数的3D网格； 

将每个单独的宇宙射线带电粒子的散射/位移协方差矩阵确定为射线路径和散射密度估计的函数；以及 

将多个宇宙射线带电粒子的概率分布确定为单独的带电粒子的概率 的乘积。 

13.如权利要求11所述的方法，其中，使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实质最大似然性估计包含： 

采集每一个宇宙射线带电粒子的散射和动量的测量； 

估计每一带电粒子与所述统计多重散射模型的每一个基函数的相互作用的几何结构； 

对每一带电粒子基函数对，确定权重矩阵：W_ij； 

利用猜测初始化每一基函数中的散射密度；以及 

迭代地求解物体体积容量的近似最大似然性解； 

其中，在预定数目的迭代处或当所述解变得小于预定容限值时，停止所述迭代过程。 

14.如权利要求9所述的方法，其中，所述带电粒子包含μ介子。 

15.一种从带电粒子层析数据中探测物体体积的系统，所述带电粒子层析数据从所述物体体积中获得，所述系统包含： 

(a)用于获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的预定带电粒子层析数据的装置； 

(b)用于提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布的装置，所述概率分布基于统计多重散射模型； 

(c)用于使用期望最大化(ML/EM)算法确定物体体积散射密度的实质最大似然性估计的装置；以及 

(d)用于基于所述实质最大似然性估计输出重建的物体体积散射密度的装置， 

其中，用于提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布的装置包含： 

(g)用于基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布的装置； 

(h)用于基于所述2D概率分布获得所述带电粒子的3D概率分布的装置； 

(i)用于获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散射的概率分布的装置；以及 

(j)用于基于多重散射的定义和所述带电粒子的散射和动量测量，获得 多重散射的概率分布的装置。 

16.如权利要求15所述的系统，包含： 

用于基于所述重建的物体体积散射密度，对所述物体体积中的目标物体的(1)存在和(2)类型中的至少一个做出决定的装置。 

17.如权利要求15所述的系统，其中，所述带电粒子为入射到所述物体体积的自然宇宙射线μ介子。