首页 / 专利库 / 视听技术与设备 / 均值漂移过程 / 基于带电粒子测量的统计层析重建

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

阅读:31发布:2020-12-18

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

如下定义了其中请求保护专有财产或权利的发明实施例。具有这样描 述的、请求保护的发明为:
1.一种从带电粒子层析数据中探测物体体积的方法,所述带电粒子层析 数据从所述物体体积中获得,所述方法包含:
(a)获得与穿过物体体积的带电粒子的散射度和估计动量相对应的预 定带电粒子层析数据;
(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布, 所述概率分布基于统计多重散射模型;
(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实 质最大似然性估计;以及
(d)基于所述实质最大似然性估计输出重建的物体体积散射密度。
2.如权利要求1所述的方法,包含:
基于所述重建的物体体积散射密度,对所述物体体积中的目标物体的(1) 存在和(2)类型中的至少一个做出决定。
3.如权利要求1所述的方法,其中,提供用于期望最大化(ML/EM) 算法的带电粒子散射的概率分布包含:
(g)基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布;
(h)基于所述2D概率分布获得所述带电粒子的3D概率分布;
(i)获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散 射的概率分布;以及
(j)基于多重散射的定义和所述带电粒子的散射和动量测量,获得多重 散射的概率分布。
4.如权利要求3所述的方法,其中,基于均质物体的预先定义的散射密 度获得带电粒子的2D概率分布包含:
(k)确定材料的散射密度作为穿过所述材料的单位深度的所述带电粒子 的期望均方散射;
(l)基于高斯模型来近似带电粒子散射分布;以及
(m)基于相关联2D
高斯分布来近似带电粒子射线散射和位移。
5.如权利要求3所述的方法,其中,基于所述2D概率分布获得所述带 电粒子的3D概率分布包含:
添加坐标并定义三维路径长度;
计算3D位移;以及
定义3D协方差矩阵。
6.如权利要求5所述的方法,其中,获得穿过经由基函数特性化的非均 质物体体积的多个带电粒子的散射的概率分布包含:
建立代表离散散射密度估计的基函数的3D网格;
将每个单独的μ介子的散射/位移协方差矩阵确定为射线路径和散射密 度估计的函数;以及
将多个带电粒子的概率分布确定为单独的带电粒子概率的乘积。
7.如权利要求1所述的方法,其中,使用所述期望最大化(ML/EM) 算法确定所述物体体积散射密度的实质最大似然性估计包含:
采集对每一带电粒子的散射和动量的测量;
估计每一带电粒子与所述统计散射模型的每一个基函数的相互作用的几 何结构;
对每一带电粒子基函数对,确定权重矩阵:Wij;
利用猜测初始化每一基函数中的散射密度;以及
迭代地求解物体体积容量的近似最大似然性解;
其中,在预定数目的迭代处或当所述解变得小于预定容限值时,停止所 述迭代过程。
8.如权利要求1所述的方法,其中,使用所述期望最大化(ML/EM) 算法确定所述物体体积散射密度的实质最大似然性估计包含:
采集对每一带电粒子i=1到M的散射和动量的测量:(Δθx,Δθy,Δx,Δy,pr2)i;
估计每一μ介子与每一体素j=1到N的相互作用的几何结构:(L,T)ij;
对每一带电粒子体素对,将所述权重矩阵Wij计算为
W ij L ij L ij 2 / 2 + L ij T ij L ij 2 / 2 + L ij T ij L ij 3 / 3 + L ij 2 T ij + L ij T ij 2 ,
在每一体素中初始化散射密度的猜测λj,old;以及
对全部体素使用停止判据过程来设置λj,old=λj,new。
9.如权利要求1所述的方法,其中,所述期望最大化(ML/EM)算法 包括均值更新规则或中值更新规则。
10.如权利要求1所述的方法,其中,所述带电粒子层析数据包含宇宙 射线μ介子层析数据。
11.一种从带电粒子层析数据中探测物体体积的、计算机实施的方法, 所述带电粒子层析数据从所述物体体积中获得,所述方法包含:
(a)获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的带 电粒子层析数据;
(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布, 所述概率分布基于统计多重散射模型;
(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实 质最大似然性估计;以及
(d)输出重建的物体体积散射密度。
12.如权利要求11所述的方法,进一步包含基于所述重建的物体散射体 积密度来做出决定。
13.如权利要求11所述的方法,其中,提供用于期望最大化(ML/EM) 算法的带电粒子散射的概率分布包含:
(g)基于均质物体的预先定义的散射密度获得带电粒子的2D概率分布;
(h)基于所述2D概率分布获得所述带电粒子的3D概率分布;
(i)获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的散 射的概率分布;以及
(j)基于多重散射的定义和所述带电粒子的散射和动量测量,获得多重 散射的概率分布。
14.如权利要求13所述的方法,其中,所述带电粒子层析数据包含宇宙 射线带电粒子层析数据。
15.如权利要求14所述的方法,其中,基于均质物体的预先定义的散射 密度获得带电粒子的2D概率分布包含:
确定材料的散射密度作为穿过该材料的单位深度的具有标称动量 p0=3GeV的宇宙射线带电粒子的期望均方散射;
使用高斯模型来近似宇宙射线带电粒子散射分布;以及
通过相关联2D高斯分布来近似宇宙射线散射和位移;
其中,基于所述2D概率分布获得所述带电粒子的3D概率分布包含:
添加坐标并定义三维路径长度;
计算3D位移;以及
定义3D协方差矩阵,并且
其中,获得穿过经由基函数特性化的非均质物体体积的多个带电粒子的 散射的概率分布包含:
建立代表离散散射密度估计的基函数的3D网格;
将每个单独的宇宙射线带电粒子的散射/位移协方差矩阵确定为射 线路径和散射密度估计的函数;以及
将多个宇宙射线带电粒子的概率分布确定为单独的带电粒子的概率 的乘积。
16.如权利要求14所述的方法,其中,使用所述期望最大化(ML/EM) 算法确定物体体积散射密度的实质最大似然性估计包含:
采集每一个宇宙射线带电粒子的散射和动量的测量;
估计每一带电粒子与所述统计多重散射模型的每一个基函数的相互作用 的几何结构;
对每一带电粒子基函数对,确定权重矩阵:Wij;
利用猜测初始化每一基函数中的散射密度;以及
迭代地求解物体体积容量的近似最大似然性解;
其中,在预定数目的迭代处或当所述解变得小于预定容限值时,停止所 述迭代过程。
17.如权利要求10所述的方法,其中,使用所述期望最大化(ML/EM) 算法确定物体体积散射密度的实质最大似然性估计包含:
采集每一带电粒子i=1到M的散射和动量的测量:(Δθx,Δθy,Δx,Δy,pr2)i;
估计每一带电粒子与每一体素j=1到N的相互作用的几何结构:(L,T)ij;
对每一带电粒子体素对,将所述权重矩阵Wij计算为
W ij L ij L ij 2 / 2 + L ij T ij L ij 2 / 2 + L ij T ij L ij 3 / 3 + L ij 2 T ij + L ij T ij 2 ,
在每一体素中初始化散射密度的猜测λj,old;以及
对全部体素使用停止判据过程来设置λj,old=λj,new。
18.如权利要求17所述的方法,其中,所述停止判据过程包含: 对每一带电粒子,使用 W ij L ij L ij 2 / 2 + L ij T ij L ij 2 / 2 + L ij T ij L ij 3 / 3 + L ij 2 T ij + L ij T ij 2 , 来计算,并取逆 矩阵,
对每一带电粒子体素对,使用
S ij ( n ) = p r , i - 2 Tr ( A ij - 1 Σ H ij - A ij - 1 Σ D i H ij T Σ D i - 1 Σ D i H ij ) + p r , i - 2 D i T Σ D i - 1 W ij Σ D i - 1 D i ( p r , i 2 λ j ( n ) ) 2
= 2 λ j ( n ) + ( D i T Σ D i - 1 W ij Σ D i - 1 D i - Tr ( Σ D i - 1 W ij ) ) p r , i 2 ( λ j ( n ) ) 2 ,
来计算条件期望项:Sij,
其中,在最后步骤中Tr(AB)=Tr(BA);
并且,在更新规则中使用
S ij ( n ) = S ij , x ( n ) + S ij , y ( n ) 2 ,
来合并x和y坐标散射数据两者。
19.如权利要求18所述的方法,其中,所述带电粒子包含μ介子。
20.如权利要求18所述的方法,其中,所述ML/EM算法包括定义为 λ j ( n + 1 ) = 1 2 M j Σ i : L ij 0 S ij ( n ) . 的均值更新规则,或定义为 λ j ( n + 1 ) = 1 2 median i : L ij 0 ( S ij ( n ) ) , 的中值 更新规则。
21.一种计算机程序产品,包含:存储指令的计算机可用的数据载体, 当被计算机执行时,所述指令使所述计算机执行用于带电粒子层析数据的物 体体积密度剖面的统计重建的方法,所述方法包含:
(a)获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的预 定带电粒子层析数据;
(b)提供用于期望最大化(ML/EM)算法的带电粒子散射的概率分布, 所述概率分布基于统计多重散射模型;
(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实 质最大似然性估计;以及
(d)输出重建的物体体积散射密度。
22.一种用于经由穿过物体体积的带电粒子探测物体体积的探测系统, 包含:
第一组位置敏感探测器,其位于物体体积的第一侧,以测量向所述物体 体积入射的带电粒子的位置和角度;
第二组位置敏感探测器,其位于与所述第一侧相对的、所述物体体积的 第二侧,以测量离开所述物体体积出射的带电粒子的位置和角度;和
信号处理单元,其接收来自所述第一组位置敏感探测器的测量信号和来 自所述第二组位置敏感探测器的测量信号的数据,其中,所述信号处理单元 处理所接收的数据以便产生所述物体体积内体积散射密度分布的统计重建。
23.如权利要求22所述的系统,其中,所述信号处理单元被配置来:
(a)获得与穿过物体体积的带电粒子的散射角度和估计动量相对应的带 电粒子层析数据;
(b)基于统计多重散射模型提供带电粒子散射密度的概率分布;
(c)使用所述期望最大化(ML/EM)算法确定物体体积散射密度的实 质最大似然性估计;以及
(d)基于所述实质最大似然性估计输出重建的物体体积散射密度。
24.如权利要求22所述的系统,其中:
所述第一组和第二组粒子探测器中的每一个包括漂移管,该漂移管被安 排以允许在第一方向中的至少三个带电粒子位置测量和在不同于第一方向的 第二方向中的至少三个带电粒子位置测量。
25.如权利要求22所述的系统,其中,所述带电粒子为入射到所述物体 体积的自然宇宙射线μ介子,并且所述信号处理单元被配置来基于所述物体 体积内体积散射密度分布的统计重建来指示目标物体是否存在于所述物体体 积中。

说明书全文

临时申请的交叉引用

本PCT申请要求(1)申请号为60/855,064、标题为“SYSTEMS, METHODS AND APPARATUS FOR PARTICLE DETECTION AND ANALYSIS 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 ray muon)的散射的宇宙射线层析。来自外层空间的稳定粒子大部分为质子,这 些粒子持续不断地轰击地球。所述粒子与在高层大气中的原子相互作用,以 产生包括许多短寿命的π介子(pion)的粒子的簇射(shower),而所述π介 子衰变并产生较长寿命的μ介子。μ介子主要通过库仑与物质相互作用, 而不具有核相互作用,并与电子相比更不容易辐射。它们只通过电磁相互作 用缓慢地丢失能量。因此,许多μ介子作为高穿透性的带电辐射到达地球表 面。在海平面的μ介子通量为大约每分钟每平方厘米一个μ介子。
当μ介子运动穿过材料的时候,亚原子粒子的电荷的库仑散射干扰它的 轨迹(trajectory)。虽然总偏离量取决于几种材料性质,但是主导参数为原子 核的原子序数Z和材料密度
需要提供一种改进的方法和系统,用于从穿过体积(volume)的μ介子 或其他带电粒子重建感兴趣的体积。

发明内容

以下发明概述被提供来帮助对涉及对探测诸如像μ介子那样的带电粒子 的粒子的技术、装置和系统的一些技术特征,以及对带电粒子层析数据的物 体体积散射密度剖面(object volume scattering density profile)的统计重建的 理解,而不是为了全面的描述。对本发明各方面全面了解可以通过将整个说 明书、权利要求附图摘要作为一个整体来获得。
本发明的上述方面以及其他目标和优势现在能够像这里所描述的那样达 成。
根据一个方面,描述一种探测系统用于经由穿过物体体积的带电粒子来 探测物体体积。所述系统包括:第一组位置敏感探测器,其位于物体体积的 第一侧,来测量向物体体积入射的带电粒子的位置和度;第二组位置敏感 探测器,其位于与第一侧相对的物体体积的第二侧,来测量离开物体体积出 射的带电粒子的位置和角度;和信号处理单元,其接收来自第一组位置敏感 探测器的测量信号和来自第二组位置敏感探测器的测量信号的数据。信号处 理单元处理所接收的数据,以产生物体体积内体积散射密度分布的统计重建。
信号处理单元能够被配置来:(a)获得相应于穿过物体体积的带电粒子 的散射角度和估计动量(estimated momenta)的带电粒子层析数据;(b)基 于统计多重散射模型(multiple scattering model),来提供带电粒子散射密度 的概率分布;(c)使用期望最大化(expectation maximization)(ML/EM)算 法,来确定物体体积散射密度的实质最大似然性估计(substantially maximum likelihood 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网格(3D grid);将每个单独的μ介子的散射/位移的协方差矩阵确定为射线路径和散射 密度估计的函数;以及将多个带电粒子的概率分布确定为单独的带电粒子的 概率的乘积。
使用期望最大化(ML/EM)算法来确定物体体积散射密度的实质最大似 然性估计能够包含:采集每一个带电粒子的散射和动量的测量;估计每一带 电粒子与统计散射模型的每一个基函数的相互作用的几何结构(geometry of interaction);对每一带电粒子基函数对,确定权重矩阵:Wij;以猜测(guess) 来初始化每一基函数中的散射密度;以及迭代地求解(solving)物体体积容 量(contents)的近似最大似然性解(solution);其中在迭代的预定数目处、 或当所述解变得小于预定的容限值时,停止迭代过程。
使用期望最大化(ML/EM)算法来确定物体体积的散射密度的实质最大 似然性估计能够包含:采集每一带电粒子i=1到M的散射和动量的测量 (Δθx,Δθy,Δx,Δy,pr2)i;估计每一μ介子与每一体素(voxel)j=1到N的相互作 用的几何结构:(L,T)ij;对每一带电粒子体素对,将权重矩阵Wij计算为 W ij L ij L ij 2 / 2 + L ij T ij L ij 2 / 2 + L ij T ij L ij 3 / 3 + L ij 2 T ij + L ij T ij 2 , 在每一体素中初始化散射密度的猜测λ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-premise computer)中实现信号处理单元的层析处理部分。可替换地,可以在远程计 算机中实现信号处理单元的层析处理部分,而该远程计算机连接到诸如专用 网络或公用网络(诸如互联网)的计算机网络上。
在图1的说明性实施例中,带电粒子是宇宙射线μ介子或其他宇宙射线 带电粒子,而位置敏感探测器10为充满用于感测(sensing)带电粒子的工作 气体的漂移单元。例如,可以用漂移管来实现漂移单元,所述漂移管具有沿 着每一漂移管的纵向延伸的中心阳极线。然而,能够使用不同于漂移单元的 位置敏感传感器来探测不同于μ介子的带电粒子。此外,能够由不同于宇宙 射线的源生成带电粒子。例如,μ介子能够从加速器中作为低强度射束(low intensity 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 Δθ ( Δθ ) 1 2 π σ θ exp ( - Δθ 2 2 σ θ 2 ) , 公式(1)
虽然真实的分布与高斯型相比具有更重或更大的尾部。可以按照材料的性质 来近似地表达分布的宽度。如在S.Eidelman等人的、Phys.Lett,vol.B592,p.1, 2004的“Review of particle physics”中所评论的那样,许多研究者已给出作 为各种材料性质的函数的用于散射的经验表达式,通过参考将其公开内容并 入这里。一个特别简单的形式为
σ θ 15 MeV βcp H L rad . 公式(2)
这里,p为以MeV/c计的粒子动量,H为材料深度,并且Lrad为材料的辐射长 度,βc为速度(c为光速),并且使用i、β=1的近似。辐射长度随原子序数 和材料密度的增加而减小。我们建立标称(nominal)μ介子动量,p0,并将 具有辐射长度Lrad的材料的散射密度定义为:
λ ( L rad ) ( 15 p 0 ) 2 1 L rad . 公式(3)
因而,材料的散射密度λ代表穿过那个材料的单位深度的、具有标称的 动量的μ介子的均方散射角度。例如,对于某些材料,以平方毫弧度每厘米 计的散射密度的值为:对铝大约为3,对大约为14,而对大约为78。因 而,穿过具有散射密度λ和深度H的材料的、具有动量p的μ介子的散射的 方差为
σ θ 2 = λH ( p 0 p ) 2 . 公式(4)

p r 2 = ( p 0 / p ) 2 , 公式(5)
得出
σ θ 2 = λ Hp r 2 . 公式(6)
将位移Δx与散射角度Δθ相关联(correlated)。一并考虑,散射角度和位 移提供了指出在较大体积中的局部散射贡献者(contributors)的位置的信息, 如图1中路径的“弯折”(kinks)所指出的。可以将散射角度和位移的分布特 性化为具有零均值的联合高斯型(jointly Gaussian),并且
σ Δx = H 3 σ Δθ , 公式(7)
ρ ΔθΔx = 3 2 . 公式(8)
我们可以将协方差矩阵表达为
Σ σ Δθ 2 σ ΔθΔx σ ΔθΔx σ Δx 2 = λ H H 2 / 2 H 2 / 2 H 3 / 3 p r 2 . 公式(9)

A H H 2 / 2 H 2 / 2 H 3 / 3 , 公式(10)
得出
Σ = λA p r 2 . 公式(11)
注意到前面提到的,能够如在图15B的流程图中所概述的那样来描述根 据一个实施例来获得单个μ介子散射的、2D的散射的概率分布(处理步骤 110)。如处理步骤150所示,按照公式(3),将材料的散射密度定义为p0=3 GeV/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是有用的。穿过所述层到所投影的(非散射的)点(xp,yp)的μ介 子路径的直线延伸(即,3-D路径长度)为
L = H 1 + tan 2 θ x , 0 + tan 2 θ y , 0 HL xy . 公式(12)
将出射的μ介子的x位置和角度定义为(x1,θx,1),然后令
Δθx=θx,1-θx,0.          公式(13)
虽然测量的x位移可以计算为xm=x1-xp,但是我们必须将这个测量旋转到正 交于射线路径的平面,并为3-D路径长度进行调整。将位移定义为
Δ x = ( x 1 - x p ) cos ( θ x , 0 ) L xy cos ( Δ θ x + θ x , 0 ) cos ( Δ θ x ) , 公式(14)
其中,中间的两项考虑3-D路径长度,并且最后一项将测量投影到适当 的朝向。
最后,将协方差权重重新定义为
A L L 2 / 2 L 2 / 2 L 3 / 3 . 公式(15)
然后,需要对散射和位移以相似的方式继续进行,并且公式(11)定义 对这两者以及坐标散射的协方差矩阵。在两个正交的、水平坐标中独立进行 散射测量。为了简化记号,我们发展(develop)对仅仅一个坐标的分析。稍 后将讨论组合来自两个坐标的信息。我们必须注意,该模型对“小的”散射 角度和位移是有效的(valid)。模型推导中忽略的二阶项对于大角度散射可能 变得显著。
注意到前面提到的,根据一个实施例,如在图15C的流程中所概述的那 样,获得延伸到3D的统计模型。首先,添加y坐标,并且定义三维路径长 度(处理步骤160、公式(12))。接下来,在处理步骤161中,按照公式(13-14) 来计算3D位移。最后,按照公式(15)表达3D协方差矩阵(处理步骤162)。
对材料的非均质体积,为了重建的目的,按照利用系数{v1,...,vj,vN}的3-D 基函数{φ1,...,φj,φN}的线性组合来代表密度剖面,即,
λ ( x , y , z ) = Σ j v j φ j ( x , y , z ) . 公式(16)
虽然对于基函数存在许多选择,但是此处我们的注意力在于矩形的3D体素。 λj被用来表示第j个基函数的系数,即,第j个体素中的散射密度。考虑图4, 示出三个层(或体素),以射线穿过所述栈(stack),传达(delivering)观测 的信息Δθ和Δx。第j个体素中“隐藏的”散射和位移分别表示为Δθj和Δxj。 再次,在图中散射的量被夸大。我们可以通过以下表达式将所观测的与隐藏 的数据联系起来
Δθ=Δθ1+Δθ2+Δθ3,              公式(17)
Δx=Δx1+L2tan(Δθ1)+Δx2+
L3tan(Δθ1+Δθ2)+Δx3
≈Δx1+Δx2+Δx3+T1Δθ1+T2Δθ2.      公式(18)
此处,我们在第二个公式中依赖小角度散射的假设,并将Tj定义为从第j个体 素的出口点到从重建体积的出口点的3-D射线路径长度。更一般地讲,对于 穿过一组体素的射线,
Δθ = Σ j Δθ j , 公式(19)
Δx = Σ j ( Δx j + T j Δ θ j ) . 公式(20)
最后,我们可以表达对第i个射线的、合计的散射/位移的协方差,通过首先 注意到,对于第j个体素
Σ ij = λ j A ij p r , i 2 , 公式(21)
其中
A ij L ij L ij 2 / 2 L ij 2 / 2 L ij 3 / 3 公式(22)
并且,Lij为第i个射线穿过第j个体素的路径长度,对未被射线“撞击”的体 素定义为零。组合公式(19)至(22),我们可以写出
Σ i = p r , i 2 Σ j N λ j W ij . 公式(23)
此处,N为体素总数,并且我们基于对要素(elements)的简单但冗长的计算, 将权重矩阵定义为
W ij L ij L ij 2 / 2 + L ij T ij L ij 2 / 2 + L ij T ij L ij 3 / 3 + L ij 2 T ij + L ij T ij 2 , 公式(24)
对未知μ介子路径进行一些假设,以便估计穿过体素的射线路径长度。 参考图5,所述近似以计算进入的和出射的轨道的最接近点(point of closest approach,PoCA)(xca,yca)开始。然后,将PoCA的入口(entry)到出口点被 连接,以便估计体素路径长度。
最后,定义数据矢量
D i Δ θ i Δx i 公式(25)
并且令D表示来自M个μ介子的全部测量。我们将给定密度剖面λ的数据似 然性写为
P ( D | λ ) = Π i M P ( D i | λ ) 公式(26)
其中利用因子
P ( D i | λ ) = 1 2 π | Σ i | 1 / 2 exp ( - 1 2 D i T Σ i - 1 D i ) . 公式(27)
注意到前面提到的,如图15D中所概述的那样,根据一个实施例,能够 获得多个μ介子穿过非均质材料的散射和位移的概率分布。首先,建立体素 的3D网格(或者其他基函数)(处理步骤170)。然后,在处理步骤171中, 按照公式(23,24)计算每一μ介子的散射/位移的协方差矩阵。最后,如处 理步骤172所示,按照公式(25-27),给定射线路径和体素散射密度,计算 全部μ介子的总体概率分布。
已描述多重散射的统计模型,现在将参考对实验效果的模型的延伸(处 理步骤113)。真实的μ介子探测器呈现有限的位置分辨率。由从拟合到多个 位置测量的轨道推导的角度和位置来特征化进入的和出射的μ介子的轨道。 因而,测量误差传播到构成μ介子层析的数据集(dataset)的散射角度和位 移测量。由RMS误差ep特性化给定探测器的精度。对于探测器的特定布置, 可以基于误差如何传播来定义误差矩阵
E e Δθ 2 e ΔθΔx e ΔθΔx e Δx 2 公式(28)。
在迭代的重建方法中,这样的误差是相对容易处理的。在我们的情况下,我 们可以通过补充公式(23)的协方差矩阵来考虑探测器误差
Σ i = E + p r , i 2 Σ j N ( λ j W ij ) . 公式(29)
以这种方式,减少噪声,否则所述噪声将由于探测器误差而出现在重建 中。对探测器误差的更加精确的模型应该考虑动量依赖性,因为跟踪误差 (tracking error)的一个源是探测器自身中的散射,并且所述散射随着粒子动 量增加而减小。如果单独的(individual)μ介子动量的估计是可用的,那 么能够对每一射线估计误差矩阵如从公式(2)所明显看出的那样, 多重库仑散射的宽度取决于粒子动量。通过在公式(5)中引入因子pr2来考虑 不同的μ介子动量。实践中,μ介子动量不是精确地已知的,但是可以从已 知散射体中的散射测量(如宇宙射线μ介子的已知光谱)中估计单独μ介子 的动量的估计。此处,假设我们有对每一μ介子的很好的估计
可以使用期望最大化算法来确定物体体积密度的最大似然性估计(方法 100的处理步骤103)。EM算法依赖按照“完整”数据(即,观测数据加上 隐藏数据)来表达“不完整”数据的似然性。在我们的应用中,观测数据 D={Di:1≤i≤M}是测量的散射。隐藏数据H={Hij: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 DLR = E H | D , λ ( n ) [ log ( P ( D , H | λ ) ) ] . 公式(30)
对于H给定和参数矢量λ(n)的条件分布(conditional distribution),给定 参数矢量λ,这个函数是观测的和未观测的数据两者的对数似然性(log likelihood)的期望值。由下列的两个步骤组成算法的每一迭代。
E步骤:估计或特性化P(H|D,λ(n)),隐藏数据的条件概率。
M步骤:最大化辅助函数Q,该辅助函数Q是关于E步骤中特性化的分布的 期望值。
在我们的情况下,由于隐藏数据唯一地确定所观测的数据,其中通过使用更 简单的辅助函数
Q ( λ ; λ ( n ) ) = E H | D , λ ( n ) [ log ( P ( H | λ ) ) ] 公式(31)
我们获得可以通过使用QDLR来获得的相同的序列估计从参数估计λ(n) 中,由如下公式,所述算法的迭代产生新的估计λ(n+1)
λ(n+1)=arg maxλ(Q(λ;λ(n)).        公式(32)
我们首先注意到单个μ介子穿过单个体素的散射的概率分布简单地从单 层的统计模型得出。
P ( H ij | λ ) = 1 2 π | Σ ij | 1 / 2 exp ( - 1 2 H ij T Σ ij - 1 H ij ) , 公式(33)
其中, Σ ij = λ j A ij p r , i 2 , 定义在公式(21)中。由于在每一体素中的散射的无 条件分布独立于在其他体素中的散射,因此合计的隐藏数据集的概率为每一 要素的概率的乘积。从而,可以将对数似然性写为
log ( P ( H | λ ) ) = Σ j N Σ i : L ij 0 ( - log λ j - H ij T A ij - 1 H ij 2 λ j p r , i 2 ) + C , 公式(34)
其中C代表不包含λ的项。考虑条件期望,我们将Q函数写为
Q ( λ ; λ ( n ) ) = C + Σ j N Q j ( λ j ; λ j ( n ) ) 公式(35)
其中利用求和项(summands)
Q j ( λ j ; λ j ( n ) ) = - M j log λ j - 1 2 λ j Σ i : L ij 0 S ij ( n ) . 公式(36)
这里,Mj是其Lij≠0的射线的数目(即,撞击第j个体素的射线的数目),并 且Sij(n)定义为
S ij ( n ) E H | D , λ ( n ) [ p r , i - 2 H ij T A ij - 1 H ij ] . 公式(37)
将公式(36)相对于λj的导数(derivative)设置为零,我们找到以下最 大化辅助函数(M步骤)的迭代公式
λ j ( n + 1 ) = 2 2 M j Σ i : L ij 0 S ij ( n ) . 公式(38)
Sij的二次型形式保证λ(n+1)的正性(positivity)。其保留以计算条件期望Sij。令 X表示随机变量Hij|Di。二次型XTA-1X的期望值为
E [ X T A - 1 X ] = Tr ( A - 1 Σ X ) + μ X T A - 1 μ X , 公式(39)
其中μX和∑X分别为X的均值和协方差。由于Di线性地依赖于Hij,因此它们 为联合高斯型。给定Di的Hij的条件分布也是高斯型,来自多变量(multivariate) 分布理论的结果。使用所述理论以及Hij与Di每个都具有零均值的事实,我们 找到
μ X = Σ D i H ij T Σ D i - 1 D i , 公式(40)
Σ X = Σ H ij - Σ D i H ij T Σ D i - 1 Σ D i H ij . 公式(41)
这里,由公式(29)给定观测数据的协方差并且可以经由公式(21)表 达隐藏数据要素的协方差我们能够执行简单的(虽然为冗长的)矩阵计 算以示出(而非明确地写出观测的隐藏数据的协方差)
Σ D i H ij A ij - 1 Σ D i H ij T = W ij ( p r , i 2 λ j ) 2 . 公式(42)
将来自公式(39)至(42)的结果替代到公式(37)中,我们找到
S ij ( n ) = p r , i - 2 Tr ( A ij - 1 Σ H ij - A ij - 1 Σ D i H ij T Σ D i - 1 Σ D i H ij ) + p r , i - 2 D i T Σ D i - 1 W ij Σ D i - 1 D i ( p r , i 2 λ j ( n ) ) 2
= 2 λ j ( n ) + ( D i T Σ D i - 1 W ij Σ D i - 1 D i - Tr ( Σ D i - 1 W ij ) ) p r , i 2 ( λ j ( n ) ) 2 , 公式(43)
其中,我们已在上一个步骤中使用Tr(AB)=Tr(BA)。
最后,为了并入x和y坐标散射数据两者,我们简单地在更新的公式(38) 中使用平均
S ij ( n ) = S ij , x ( n ) + S ij , y ( n ) 2 , 公式(44)
注意到的是公式(38)代表正在撞击体素的射线的均值。此后,将把该 公式的使用称作为均值方法。下面将示出,这个更新公式的可替换形式在减 少由于外露层(outlier)μ介子数据造成的噪声中是有用的。由所改变的更新 公式来定义算法的中值方法:
λ j ( n + 1 ) = 1 2 media n i : L ij 0 ( S ij ( n ) ) , 公式(45)
注意到前面提到的,在图15E中概述了使用期望最大化算法的、物体体 积的估计密度剖面的最大化似然性的过程(方法100的处理步骤103)。如处 理步骤180所示,对每一μ介子i=1到M采集散射和动量的测量: (Δθx,Δθy,Δx,Δy,pr2)i。估计每一μ介子与每一体素j=1至N的相互作用的几何 结构:(L,T)ij(处理步骤181)。对于每一μ介子体素对,使用公式(24)来 计算权重矩阵:Wij(处理步骤182)。以猜测λj,old来初始化每一体素中的散射 密度(处理步骤183)。如下面的那样,在处理步骤184中指示停止判据。对 于每一μ介子,使用公式(29)并且取逆矩阵(inverse taken)来计算对 每一μ介子体素对,使用公式(43、44)计算条件期望项:Sij,使用公式(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.8mrad2/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%)。给定在模拟和反演模型(inversion model)之间的匹配,这个结果使反演算法和实施有效。
重建表现为与物体场景相同。相应于三个物体中的每一个的、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),或为了限制本发明的范 围。按照上面的教导,许多修改和改变都是可能的,而不背离以下权利要求 书的范围。可以预期的是本发明的使用会涉及具有不同的特性的组成部分。
高效检索全球专利

专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台,是提供专利分析,专利查询,专利检索等数据服务功能的知识产权数据服务商。

我们的产品包含105个国家的1.26亿组数据,免费查、免费专利分析。

申请试用

分析报告

专利汇分析报告产品可以对行业情报数据进行梳理分析,涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度,系统通过AI智能系统对图表进行解读,只需1分钟,一键生成行业专利分析报告。

申请试用

QQ群二维码
意见反馈