CN109685865B - 适合直线扫描轨迹的锥束断层重建方法 - Google Patents

Abstract

本发明提供一种适合直线扫描轨迹的锥束断层重建方法，其特点是X射线源和检测器间无相对运动，与待检测物体间存在均匀的直线运动，系统运动过程中，等时间间隔获得沿直线扫描轨迹的投影数据；对待检测物体进行切片式网格划分，计算切片式网格的全变差系数矩阵和投影权重贡献系数矩阵；对切片网格全变差的梯度进行开方运算，以开方后的全变差为目标函数，X射线投影方程为约束，将断层重建问题转化为凸二次规划；进一步将凸二次规划求解问题转化为增广拉格朗日函数的优化求解问题，利用神经网络进行迭代求解。该发明的方法具有扫描机械结构简单、可获得较清晰物体的三维影像，一定程度上解决物体互相遮挡的问题等优点。

Description

适合直线扫描轨迹的锥束断层重建方法

技术领域

本发明属于无损检测领域，具体来讲公开了一种适合直线扫描轨迹的锥束CT重建方法，适用于X射源和检测器沿直线运动的计算机断层成像系统。

背景技术

计算机断层成像(Computed Tomography，简称CT)是通过无损方式获取物体内部结构信息的一种技术手段，已被广泛地应用到医疗诊断、工业无损探伤等领域。近年来，国内外时有暴力恐怖袭击事件的发生，对诸如机场、地铁站、火车站、海关等的安检系统提出越来越高的要求。目前阶段，基于安检系统主要使用双能X射线安检机和双视角X射线安检系统。然而，无论是那种方式，计算机系统所显示是具有一定立体感的透视图像，不是真正意义上的断层成像，不能完全有效地分辨出行李中重叠遮挡的问题。

断层图像首选的射源扫描轨迹是圆轨迹。在医学领域内，通过复合检查床的直线运动，获得了满足锥束精确重建条件的螺旋扫描轨迹，从而获得人体内部比较精确的重建结果。直线扫描轨迹不满足断层图像精确重建的扫描轨迹条件，理论上采用适合完备轨迹的锥束重建算法，无法获得精确的重建结果。安检系统中可以采用医学领域内螺旋扫描的形式进行危险物品的检测，从而可以获得待检查物体内部的精确影像。然而，这样的成本较高，检测的速度也会受到影响，造成人流堵塞。从这两个角度而言，直线轨迹运动是最适合安检系统。本发明利用全变差作为正则项能够在欠采样的条件下很好地重建图像，提出了将重建图像的全变差构造矩阵表达形式，将X射线衰减规律用一组线性方程组来表示，从而将断层重建算法转换为带有等式约束的凸二次规划优化求解问题，该算法适合X射源直线扫描轨迹，并且在投影角度方向上允许稀疏采样。

发明内容

本发明的目的是针对现有技术的不足而提供一种适合直线扫描轨迹的锥束断层重建方法，其特点是X射线源和检测器间无相对运动，与待检测物体间存在均匀的直线运动，系统运动过程中，等时间间隔获得沿直线扫描轨迹的投影数据；对待检测物体进行切片式网格划分，计算切片式网格的全变差系数矩阵和投影权重贡献系数矩阵；对切片网格全变差的梯度进行开方运算，以开方后的全变差为目标函数，X射线投影方程为约束，将断层重建问题转化为凸二次规划；进一步将凸二次规划求解问题转化为增广拉格朗日函数的优化求解问题，利用神经网络进行迭代求解。该发明的方法具有扫描机械结构简单、可获得较清晰物体的三维影像，一定程度上解决物体互相遮挡的问题等优点。

本发明的目的由以下技术措施实现：一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；重建方法的步骤包括：根据检测系统结构参数以及运动参数，对待检测物体进行切片式网格划分，计算切片式网格的全变差系数矩阵和投影权重贡献系数矩阵，运算结果存储到存储器中；待检测物体以均匀速度与检测系统形成直线扫描轨迹，在等间隔时间内，检测系统采集射线源发射的X射线，并转换为投影数据；复合全变差系数阵、投影权重系数阵和投影数据，建立凸二次规划目标函数及约束条件，利用增广拉格朗日法将优化问题转化为无约束的优化问题；利用神经网络算法对无约束优化问题进行求解，当网路达到平衡状态后，切片式网格的断层图像，通过反插值获得被检测物体的断层图像，实现直线扫描轨迹的锥束重建。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的检测系统是指由X射线射源和检测器组成的检测系统，射源和检测器的中心在一条水平直线上。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的直线扫描轨迹是指X射线射源和检测器在工作过程中无相对运动，两者与被测物体存在相对直线运动，射源在等间隔时间内发射X射线。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的全变差系数矩阵是指对被测物体依检测器间距进行切片网格划分，每个网格的像素值代表被检测物体的衰减系数，所有网格点与其相邻的两个网格点的衰减系数差的平方和进行累加，二次项的系数所构成的稀疏矩阵正定矩阵。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的投影权重贡献系数阵是指建立从射源到检测器每个阵元的X射线投影方程，投影方程的系数是切片网格对投影射线的权重贡献值，投影方程的系数构成了投影权重贡献系数矩阵。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的投影数据是指X射线入射强度与接受强度的比值的对数值。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的凸二次规划目标函数是全变差系数矩阵与各重建网格的衰减系数构成满足二次规划的正定条件的矩阵表达形式。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的约束条件指投影权重贡献系数矩阵与各重建网格的衰减系数构成线性方程组。

根据本发明所述一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的梯度神经网络算法是指无约束的优化问题的优化变量进行梯度表达，利用神经网路算法进行迭代求解。

与现有技术相比，本发明具有如下优点：本发明适合直线扫描轨迹的锥束重建成像方法，该技术方案采用的检测系统与被检测物体间只有直线运动，扫描结构简单，适合长物体断层重建，特别适合X射线安检系统；重建物体沿射源与检测器水平阵元方向所在的平面进行划分，加快了重建速度；修改图像全变差的形式，将锥束断层重建问题转化为无约束的优化问题，利用梯度神经网络进行优化问题求解，优化问题的原形满足正定的凸二次规划，理论上可以获得全局最优解，实现了直线扫描轨迹的锥束重建，一定程度上解决了透视成像中的物体重叠问题。

附图说明

图1直线扫描轨迹的锥束断层重建方法的流程图；

图2直线扫描重建系统结构和运动示意图；

图3被检物体的网格划分示意图。

其中，(a)被检物体的等间距网格划分示意图；(b)依检测器间距的非等距切片划分示意图；(c)非等距的切片网格与被检物体的等间距网格关系示意图；

图4切片网格的全变差系数分析过程图；

其中，(a)切片网格元素的全变差计算关系图；(b)网格元素与相邻元素的二次项系数分布图；

图5确定切片网格对投影数据的权重贡献值示意图；

其中，(a)均匀介质中X射线衰减规律；(b)非均匀介质中X射线衰减规律；(c)投影射线的权重贡献值示意图；(d)t₁时刻投影射线l的权重贡献值示意图；(e)t₂时刻投影射线l的权重贡献值示意图；

图6神经网络求解优化求解的架构流程图。

具体实施方式

以下通过实施例对本发明进行具体的描述，有必要在此指出的是本实施例只用于对本发明进行进一步说明，不能理解为对发明保护范围的限制，该领域的技术熟练人员可以根据上述本发明的内容作出的一些非本质的改进和调整。

下面对照附图详细描述本发明的实施例。

图1示出了本发明直线扫描轨迹的锥束断层重建方法的流程图。其扫描和重建的具体流程如下：

根据检测系统结构参数以及运动参数，对待检测物体进行切片式网格划分，计算切片式网格的全变差系数矩阵和投影权重贡献系数矩阵，运算结果存储到存储器中；

待检测物体以均匀速度与检测系统直线扫描轨迹，在等间隔时间内，直线扫描系统采集射线源发射的X射线，并转换为投影数据；

复合全变差系数阵、投影权重系数阵和投影数据，建立凸二次规划目标函数及约束条件，利用增广拉格朗日法将优化问题转化为无约束的优化问题；

利用神经网络算法对无约束优化问题进行求解，当网路达到平衡状态后，切片式网格的断层图像，通过反插值获得被检测物体的断层图像，实现直线扫描轨迹的锥束重建。

图2示出了直线扫描重建系统结构和运动示意图。X射线射源和检测器组成的检测系统，检测器是二维平板检测器，射源和检测器的中心在一条水平直线上。X射线射源和检测器在工作过程中无相对运动，两者与被测物体存在相对直线运动，如图中箭头所示。射源在等间隔时间内发射X射线，检测器接受X射线信号，该信号经一定运算可获得直线扫描轨迹的投影数据。通过检测系统和被测物体的相对运动，可以获得不同角度的投影数据，用于直线扫描轨迹的锥束重建。

图3示出了被检物体的网格划分示意图。在计算机系统中，通常用二维数组或矩阵来描述一幅图像。

图3(a)示出了被检物体的等间距网格划分示意图。通常情况下，物体是沿三维方向进行等间距的网格划分。因为要计算网格的全变差，获取全变差的系数矩阵。矩阵的维数是原网格数量的平方数量。如果直接用于三维等间距网格，数据量过大，不利于问题求解。

图3(b)示出了依检测器间距的非等距切片划分示意图。射源与位于同一水平线的阵元的上下边界在三维空间中形成了很薄的立体空间。被检测物体被该立体空间截切成非常薄的切片。该切片倾斜于射源与检测器水平中线所在的平面，沿竖直方向切片的间距不等。非等距的切片在随被测物体与检测系统相对直线运动中始终位于同一平面。

图3(c)示出了非等距的切片网格与被检物体的等间距网格关系示意图。对倾斜非等距的切片进行网格划分。保证网格在水平上的投影与被测物体的等间距水平网格，这样的划分保证在基于全变差的重建过程仅发生在切片网格上。重建结束后可通过反插值操作获得被测物体的断层重建图像。

图4示出了切片网格的全变差系数分析过程图。在一幅图像中，全变差可离散为以下的近似公式：

其中f表示图像的像素值，该公式表达了图像的梯度方向的变化。

图4(a)示出了切片网格元素的全变差计算关系图。全变差公式表明全变差是网格元素与相邻元素作差的平方和的平方根之和。重建的目标是切片的衰减系数，因此用衰减系数代替图像的像素值。同时为了得到方便进行求解，去除根号操作后，相应的表达式为：

为后续描述方便，上式称为切片网格的衰减系数平方全变差。

图4(b)示出了网格元素与相邻元素的二次项系数分布图。对衰减系数平方全变差进行全展开合并同类项，该式是仅含有二次项的多项式。而二次项是具有局部属性的，表达了每个切片网格元素最多与其相邻的四个元素存在关联，在切片的四个顶点和边界点因相邻元素个数不同而发生退化。不失一般性，以切片网格非边界点为例，说明衰减系数平方全变差中二次项系数分布的含义。中心点处数字4，表明含有

的系数为4；左侧处数字-2，表明含有μ_i-1,jμ_i,j的系数为-2；右侧处数字-2，表明含有μ_i,jμ_i+1,j的系数为-2；上侧处数字-2，表明含有μ_i,j-1μ_i,j的系数为-2；上侧处数字-2，表明含有μ_i,jμ_i,j+1的系数为-2。

得到二次项的系数分布后，很容易用矩阵形式表达衰减系数平方全变差，其公式为：

式中U表示切片网格的衰减系数形成的列向量；

是二次项系数矩阵，因衰减系数平方全变差是平方和，显然无论U取何值，/>

所以/>

一定是正定矩阵，并且该矩阵为稀疏矩阵，仅在五条对角线存在非零值。除主对角线上的竖直会在2,3,4有规律变化外，其余四条次对角线上的值均为-1。若U的维数为n，四条次对角线的索引值分别为-n-1，-1,1，n+1。负数表达次对角线在主对角线的下方，正数表达次对角线在主对角线的上方。显然，Q是对称正定矩阵。

为方便优化求解，令

并且Q也是对称正定矩阵，则衰减系数平方全变差公式变为：/>

图5确定切片网格对投影数据的权重贡献值示意图。X射线衰减规律服从Beer定理。

图5(a)示出了均匀介质中X射线衰减规律。在衰减系数为μ的均匀介质中，入射强度为I₀的射线穿过介质Δx距离后，强度变为I，其公式为：

图5(b)示出了非均匀介质中可对介质内部进行适当等间距网格划分。各段网格内的衰减系数分别为μ_i,i＝1,2,3,……n，网格的厚度记为Δx_i，则根据Beer定理，X射线强度的衰减关系可表示为：

上式即表示X射线在非均匀介质中沿直线路径传播时的线性衰减规律。在CT系统中，投影数据是X射线的入射强度I₀与检测器测量的强度I之间的比值经过对数运算后得到的。

图5(c)示出了投影射线的权重贡献值示意图。当物体经过切片网格划分后，投影数据表示为衰减系数和距离的有限积分和形式。每个网格中衰减系数被认为是常数，距离Δx_i被看成X射线在切片网格中沿直线路径传播时对投影射线的贡献值，为后续表达方便记为ω_i,j，投影数据记为P，有限积分和也称之为投影方程。

图5(d)和(e)示出了t₁和t₂两个时刻同一投影射线上的权重贡献值示意图。断层成像利用X射线对物体从不同的角度进行照射，将检测器探测的信号转化为物体的一系列投影数据，经过计算机重建后得到物体所得到衰减系数μ(x)分布。这就需要在不同方向上采集了若干条投影数据。通过不同时间的同一条投影射线穿过物体的网格不同，获得不同的投影方程，其表达如下：

其中脚标k上式可进一步改写成矩阵形式：

WU＝P

上式称之为投影方程组，至此获得了投影权重贡献系数矩阵W。

图6示出了神经网络求解优化求解的架构流程图。二次规划是非线性优化中的一种特殊情形，它的目标函数是二次实函数，约束函数都是线性函数。从上面的推导过程中获得了衰减系数平方全变差的矩阵表达以及描述X射线衰减规律的投影方程组。将这两项写成二次规划的形式可以获得：

其中，U是含有n个元素的列向量，Q∈R^nxn是对称正定矩阵，W∈R^mxn的矩阵，P是含有n个元素的列向量，WU＝P是满足行满秩的线性方程组。上式满足凸二次规划的求解问题，目标函数为衰减系数平方全变差，等式约束为X射线衰减规律的投影方程组，通过求解方法，可获得最优解U^*。

求解有约束的凸规划时，关键是优化问题转化为无约束问题，通过构造增广拉格朗日函数，求其最小值，其表达式为：

式中，λ^T＝(λ₁,λ₂,...,λ_m)是含有m个元素的行向量，被称为拉格朗日乘子。ψ是较大的常数；为加快收敛速度，减轻解的振荡，引入附加惩罚项

根据优化理论，获得最优解的必要条件为:

上式中，

表示增广拉格朗日函数对U方向导数；/>

增广拉格朗日函数对λ方向导数，两个表达式具体为：

为使增广拉格朗日函数取得局部极小值，必须向函数上当前点对应梯度的反方向的规定步长进行迭代搜索，即函数沿着导数的方向解的收敛速度最快，U和λ的递推公式为：

λ(k+1)＝λ(k)-η[WU(k)-P]

上式中，0＜ζ＜1和0＜η＜1表示搜索的步长。一旦有变量的迭代关系，很容易利用神经网络进行求解。神经网络的架构流程图如图6所示，其中输入量是U和λ。当λ的改变量趋于0时，网络达到平衡，此时的输入量即为增广拉格朗日的最优解，也即获得切片网格的衰减系数分布。

获得所有切片网格的衰减系数分布后，被测物体的等距三维网格每个点向检测器投影，以确定切片网格的位置。因为划分切片网格时，保证水平方向网格划分与等距三维网格划分一致，确定切片网格的位置后，只要进行竖直方向的发插值操作就可以获得被测物体的重建图像。

通过上述步骤即可所获得直线扫描轨迹上被测物体的重建图像。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (4)

1.一种适合直线扫描轨迹的锥束断层重建方法，其特征在于；重建方法的步骤包括：

步骤1：根据检测系统结构参数以及运动参数，对待检测物体进行切片式网格划分，计算切片式网格的全变差系数矩阵和投影权重贡献系数矩阵，运算结果存储到存储器中；

所述的全变差系数矩阵是指对被测物体依检测器间距进行切片网格划分，每个网格的像素值代表被检测物体的衰减系数，所有网格点与其相邻的两个网格点的衰减系数差的平方和进行累加，二次项的系数所构成的稀疏矩阵正定矩阵；

所述的投影权重贡献系数阵是指建立从射源到检测器每个阵元的X射线投影方程，投影方程的系数是切片网格对投影射线的权重贡献值，投影方程的系数构成了投影权重贡献系数矩阵；

步骤2：待检测物体以均匀速度与检测系统直线扫描轨迹，在等间隔时间内，直线扫描系统采集射线源发射的X射线，并转换为投影数据；

所述的投影数据是指X射线入射强度与接受强度的比值的对数值；

步骤3：复合全变差系数阵、投影权重系数阵和投影数据，建立凸二次规划目标函数及约束条件，利用增广拉格朗日法将优化问题转化为无约束的优化问题；

所述的凸二次规划目标函数是全变差系数矩阵与各重建网格的衰减系数构成满足二次规划的正定条件的矩阵表达形式；

所述的约束条件指投影权重贡献系数矩阵与各重建网格的衰减系数构成线性方程组；

步骤4：利用神经网络算法对无约束优化问题进行求解，当网路达到平衡状态后，切片式网格的断层图像，通过反插值获得被检测物体的断层图像，实现直线扫描轨迹的锥束重建；

网络达到平衡，此时的输入量即为增广拉格朗日的最优解，也即获得切片网格的衰减系数分布；

获得所有切片网格的衰减系数分布后，被测物体的等距三维网格每个点向检测器投影，以确定切片网格的位置；因为划分切片网格时，保证水平方向网格划分与等距三维网格划分一致，确定切片网格的位置后，只要进行竖直方向的发插值操作就可以获得被测物体的重建图像。

2.根据权利要求1所述的适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的检测系统是指由X射线射源和检测器组成的检测系统，射源和检测器的中心在一条水平直线上。

3.根据权利要求1所述的适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的直线扫描轨迹是指X射线射源和检测器在工作过程中无相对运动，两者与被测物体存在相对直线运动，射源在等间隔时间内发射X射线。

4.根据权利要求1所述的适合直线扫描轨迹的锥束断层重建方法，其特征在于；所述的神经网络算法是指无约束的优化问题的优化变量进行梯度表达，利用神经网路算法进行迭代求解。

Images