一种基于X射线源阵列成像的投影图像恢复方法与流程

本发明涉及一种X射线成像方法，具体涉及一种基于X射线源阵列成像的投影图像恢复方法。

背景技术：

随着X射线成像在医学、安检、无损检测、工业探伤等领域的应用越来越广泛，人们对于降低X射线成像的剂量、提高成像分辨率和成像速度的要求也愈加迫切。传统的X射线成像受制于单点源X射线锥角的限制，需要射线源与被探测物体有足够远的距离才能覆盖待检测区域，加之对机械扫描速度的追求，使得X射线成像系统复杂、体积较大、成本高昂，应用范围受到了局限。其中，成像系统体积较大的问题还导致在对部件进行检测或对病人进行临床检查时，难以避免地使得成像目标区域之外的部分受到X射线辐射，可能对待测部件或人体产生不必要的损伤。

当前，X射线源正开始朝着冷阴极和平板化的方向发展，其可以实现一个器件内单独寻址大量微X射线源阵列，同时，由于X射线源能够与被测目标紧密耦合，避免对相关区域的辐射，有望依此降低X射线剂量、减少辐射泄露，并缩小成像系统体积。自2001年日本名古屋工业大学最先报导了以碳纳米管(CNT)作为电子源的X射线管以来，微焦X射线源及其在动态成像和分布式X光源CT系统的研究成为热点，其中美国北卡罗来纳大学研制的分布式冷阴极X射线管乳腺CT原理型样机已经进入临床验证阶段。而到了2015年，中山大学已报导了较大面积的氧化锌纳米冷阴极平板X射线源，并实现了小于25微米的静态成像。等等一系列研究成果都为利用X射线源阵列探索新的成像模式与图像恢复方法奠定了基础。

技术实现要素：

为了解决现有技术中的问题，本发明提出一种基于X射线源阵列成像的投影图像恢复方法，能够利用X射线源阵列微焦点和多射线源的特点，实现快速高分辨率X射线成像，同时大大缩小X射线成像系统体积，其中，适应待检测物体形状大小来组装成像模块使得X射线成像能满足更多的应用场景的要求，辐射泄露降低，而利用锥束射线空间重排的图像恢复方法，结合分时和同时扫描能够实现快速高分辨率X射线成像，具有很好的实用前景。

为了实现以上目的，本发明所采用的技术方案为：包括以下步骤：

1)使用X射线源单元组成线阵或面阵的X射线成像扫描结构，依据被检测对象的形状规格来组合射线源单元与平板探测器模块；

2)使用X射线源阵列上的多个射线源单元分时或同时发出锥束射线穿过被检测物体，平板探测器模块接收被检测物的锥束投影图像，得到投影数据；

3)根据投影数据建立被检测物的分时或同时投影物理模型；

4)利用锥束射线空间插值重排法在分时或同时投影物理模型下，将锥束投影图像恢复为平行X光束投影图像，完成投影图像恢复。

所述步骤2)中得到投影数据的具体过程为：首先X射线源阵列上的Q个射线源单元，依照顺序先后发出锥束X射线穿过被检测物体，其中第q个射线源在探测器单元p上得到光子强度I′_pq，然后利用对应路径空扫数据I_pq，依据Beer定理，分别得到Q个源中射线穿过物体路径上衰减系数的线积分投影，其中第q个射线源在探测器单元p上得到投影值为y_pq。

所述步骤2)中得到投影数据的具体过程为：首先将X射线源阵列上的Q个射线源单元中发出锥束在空间上无重合部分的分为同一组，得到K个射线源组，依照组号先后发出射线穿过被检测物体，其中第k个射线源组在探测器单元p上得到光子强度I′_kq，然后利用对应路径空扫数据I_pq，依据Beer定理，得到K个射线源组中射线穿过物体路径上衰减系数的线积分投影，其中第q个射线源在探测器单元p上得到投影值为y_kq。

所述步骤2)中得到投影数据的过程为：X射线源阵列上的Q个射线源单元，同时发出锥束穿过被检测物体，探测器单元p上得到光子强度y_p，其为各射线源发出射线穿过不同衰减路径之后到达p处光子强度之和。

所述步骤3)中分时投影物理模型建立过程：射线源阵列上的Q个射线源单元依次发出射线，结合空扫描数据得到Q个线积分投影图，或者将射线源阵列上的Q个射线源单元分为锥束覆盖空间不重合的K个射线源组，每组依次发出射线，结合对应空扫描数据，根据Beer定理得到K个线积分投影图；然后根据下列公式(1)建立分时投影物理模型：

其中，q表示射线源的位置，f∈R^N表示待成像物体的线性衰减系数，a_pqn表示由射线源q和探测器p形成的X线光束与体素f_n之间交集的体积，n_pq表示线积分投影图上对应点处的噪声分量。

所述步骤3)中同时投影物理模型建立过程：射线源阵列上的Q个射线源单元同时发出射线，平板探测器接收穿过被检测物体的X射线强度，然后根据下列公式(2)建立同时投影物理模型：

其中，q表示射线源的位置，f∈R^N表示待成像物体的线性衰减系数，a_pqn表示由射线源q和探测器p形成的X线光束与体素f_n之间交集的体积，I_pq表示射线源q对应于探测器p方向的入射光子强度，s_p与n_p分别表示位置为p的探测器接收到的散射分量和噪声分量。

所述步骤4)中空间插值重排法的具体过程如下：

(a)定义从点源q射线到探测器单元p的射线途径上，成像物体f的线性衰减系数线积分为：

即为函数f关于路径pq的Radon变换；

(b)定义虚拟射线源阵列，其所在平面与实际射线源阵列平面重合，平行于探测器平面，相距为d，虚拟平板射线源上存在与探测器单元数目P相同的射线源单元，各自发出垂直于对应探测器的强度为I′的笔形射线束，穿过成像物体f被衰减后只到达对应探测器单元，即虚拟射线源p′发出的射线穿过物体，在对应探测器单元p处接收到的虚拟点源信号为：

其中，b_p′pn为笔形束射线p’p与体素f_n相交的体积，同时定义点源p′到探测器p形成的笔形束射线穿过成像物体，在穿过路径上对该虚拟射线进行衰减的线性衰减系数总和为：

由于位置p′和p是一一对应的，式(5)可表示为：

(c)定义从g′_p计算g_pq的线性变换矩阵w_pq，为虚拟投影g′_i向路径pq投影g_pq的线性加权系数，于是有：

(d)对于分时投影物理模型，有：

y_pq＝w_pqig′_i+n_pq (8)

建立优化目标函数，有：

求解目标函数，得到分时投影模型下基于空间插值重排法的图像恢复结果；

对于同时投影物理模型，采用矩阵D∈R^P×(P×Q)表示I_pq产生的效果，有：

建立优化目标函数，有：

式中，λ表示权重因子，表示散射信号的梯度，β表示正则化因子，R表示正则项，Y为投影信号，S为散射信号；

求解目标函数，得到同时投影模型下基于空间插值重排法的图像恢复结果。

所述步骤4)中已知射线源单元q，探测器单元p，取p、q线的中点M，沿M作阵列射线源平面和探测器平面的垂线，对应射线路径p′p为虚拟射线源到对应探测器单元的笔形束，由公式(3)得从点源q到探测器单元p的射线途径上成像物体f的线性衰减系数线积分g_pq，考虑pq间的欧氏距离为‖p-q‖，探测器平面与射线源阵列平面间距d，则相应虚拟射线路径对应的衰减系数线积分为：

g_pq表示为：

所述步骤d)中采用加速迭代法求解目标函数：令Ψ(x)表示目标函数，x为求解目标，t表示加速因子，z表示加速项，初始化重建目标及加速项x⁽⁰⁾＝z⁽⁰⁾，t₀＝1；

根据上述公式，并设定收敛条件，迭代求解目标函数得到图像恢复结果。

所述步骤d)中采用基于泰勒展开法求解目标函数，采用下列公式：

式中g′₀为求解目标g′某一次迭代的更新值，w_pq为对于射线路径pq，从理想投影恢复实际投影值的变换矩阵，对于同时投影的情况，依照公式(17)的近似方法，选择迭代更新机制，使g′₀不断接近求解目标，得到图像恢复结果。

与现有技术相比，本发明使用X射线源单元组成线阵或面阵的X射线成像扫描结构，依据被检测对象的形状规格来组合射线源单元与平板探测器模块，使用X射线源阵列上的多个射线源单元分时或同时发出锥束射线扫描待检测物体，每个射线源发出的射线仅覆盖待检测物体的一部分，所有射线源发出的射线在空间上完整覆盖整个待检测区域，探测器接收透过被检测物体的X射线，得到分散或混叠的投影图像，本发明能够依据待检测物体大小调整成像系统的体积，利用阵列射线源的结构显著减小射线源与待检测物体的距离，适应待检测物体的不同形状规格，降低系统的设计复杂度及成本，丰富系统的应用场景，避免X射线对不相关区域的辐射。使用锥束射线空间重排的方法在分时或同时投影模型下，将锥束投影图像恢复为理想的透过物体到达每一个探测器单元的平行X光束投影，利用锥束射线空间重排的图像恢复方法，结合分时和同时的扫描方法能够实现快速高分辨率X射线成像，能够使得X射线成像适应更多的应用场景，具有很好的实用前景。

进一步，采用虚拟射线源投影到实际阵列源射线投影的简化线性变换方法可以在误差允许的范围内，降低空间插值重排的运算量，提高计算速度，实现快速成像。

进一步，采用加速迭代法，设定收敛条件，迭代求解目标函数，能够较快的收敛得到投影恢复图像。

进一步，采用基于泰勒展开法对同时投影模型下的图像恢复目标函数，从理想投影恢复实际投影值的变换矩阵，对于同时投影的情况，选择适当的迭代更新机制，使g′₀不断接近求解目标，可以在误差允许的范围内，简化目标函数的求导计算，减少计算量，较快的收敛得到投影恢复图像。

附图说明

图1是本发明提出的使用射线源阵列的X射线成像模式的实施例的立体结构示意图，图中，1、射线源阵列组合模块；2、射线源单元；3、被检测物放置区间；4、被检测物体；5、射线阻挡模块；6、平板探测器组合模块；

图2是本发明提出的通用性的射线源阵列与平板探测器组合模块立体结构示意图，图中，7、8、9铰接结构；

图3是本发明提出的使用射线源阵列锥束射线空间重排的方法将投影图像恢复为理想平行光束投影的效果示意图，图中，10、虚拟射线源阵列；11、虚拟射线源单元；

图4是本发明提出的使用射线源阵列上的多个射线源或射线源组分时投影流程图；

图5是本发明提出的使用射线源阵列上多个射线源同时投影流程图；

图6是本发明提出的对射线源阵列锥束投影进行空间重排的方法将投影图像恢复为理想平行光束投影的算法流程图；

图7是本发明提出的一种虚拟射线源投影到实际阵列源射线投影的简化线性变换方法，其为本发明提出的基于空间插值重排的图像恢复方法的一种特例；

图8是本发明提出的一种加速迭代求解目标函数的图像恢复方法，其为在本发明提出的目标函数下求解方法的一种特例；

图9a～9d是本发明的一组范例所使用的模体及其理想投影，其中各标记如下：9a～9c分别为所用模体的中心横断面、冠断面和矢断面图；9d为垂直模体矢断面方向理想投影图；

图10是本发明在图9a～9d示模体下，使用分时投影方法获得投影数据，建立分时投影模型图像恢复目标函数，使用简化的空间插值重排方法，利用图8示加速迭代求解方法获得的一组图像恢复结果；

图11是本发明在图9a～9d示模体下，使用同时投影方法获得投影数据，建立同时投影模型图像恢复目标函数，使用简化的空间插值重排方法，并对图像恢复目标函数使用基于泰勒展开的近似方法，利用图8示加速迭代求解方法获得的一组图像恢复结果。

具体实施方式

下面结合具体的实施例和说明书附图对本发明作进一步的解释说明。

本发明使用X射线源阵列的X射线成像扫描结构，依据被检测对象的形状规格来组合射线源单元与平板探测器模块，使用X射线源阵列上的多个射线源分别或同时发出锥束射线得到投影图像的扫描方案，使用在分时或同时投影模型下平板射线源锥束射线空间重排的方法将投影图像恢复为理想平行光束投影。以上成像方法、适形结构设计方法及图像恢复方法，能够缩小X射线成像系统体积，适应更多的应用场景，降低辐射泄露，实现快速高分辨率X射线成像。

参见图1，使用图1所示使用X射线源阵列的X射线成像模式时，根据待检测物体大小规格选择图2所示基本单元模块进行连接，使得各射线源单元发出锥束所覆盖的空间可以容纳待检测物，由电脑主机控制调控电路，应用图4或图5所示分时和同时投影方案对物体进行扫描，平板探测器接收X射线穿过物体后的投影数据，完成投影后输入电脑主机计算单元，利用图6所示锥束投影恢复算法将投影数据恢复为理想平行光束投影，经过图像后处理，显示输出。其中，对求解目标迭代在迭代中重排为锥束射线投影时可使用图7所示一种虚拟射线源投影到实际阵列源射线投影的简化线性变换方法，对求解目标迭代更新时可使用如图8所示一种加速迭代求解目标函数的图像恢复方法。

下面结合图1与图4详述分时扫描方案获得投影的过程：

(a)X射线源阵列上的Q个射线源单元，在不放置被探测物体时(空扫描)，依照顺序先后发出锥束，其中第q个射线源在探测器单元p上得到光子强度I_pq，在放置被检测物体时，依照顺序先后发出锥束，其中第q个射线源在探测器单元p上得到光子强度I′_pq，依据Beer定理，可得到Q个源中射线穿过物体路径上衰减系数的线积分投影，其中第q个射线源在探测器单元p上得到投影值为y_pq；

(b)将X射线源阵列上的Q个射线源单元中发出锥束在空间上无重合部分的分为同一组，得到K个射线源组，在不放置被探测物体时(空扫描)，依照组号先后发出射线，其中第k个射线源组在探测器单元p上得到光子强度I_pk，在放置被检测物体时，依照组号先后发出射线，其中第k个射线源组在探测器单元p上得到光子强度I′_kq，依据Beer定理，得到K个射线源组中射线穿过物体路径上衰减系数的线积分投影，其中第q个射线源在探测器单元p上得到投影值为y_kq。

上述即使用分时扫描方案获得投影的过程，获得的投影图像由一系列来自不同射线源或不同射线源组的投影图像，每个射线源发出的锥束射线仅穿过被检测物体的一部分，所有射线源发出的锥束在空间上完全覆盖完整的被检测部分，显然与传统的需要被检测物与射线源保持足够远距离才能成像的模式不同，应用图1示结构可以明显的缩小成像系统体积，从而降低辐射泄露，扩大X成像的应用场景。同时X射线源阵列分组扫描的方式使得投影图像的获取更加快速，减小了被检测物体可能的移动为成像结果带来的影响。在医疗诊断或小动物X射线成像的用途中，分时方法的投影频率可根据如心跳呼吸等生物信号来控制，计算机得到触发信号后，判断探测器应该曝光，若满足条件，则获取一次投影，这种控制方式能进一步消除生理运动所造成的运动伪影。

下面结合图1与图5详述同时扫描方案获得投影的过程：X射线源阵列上的Q个射线源单元，在不放置被探测物体时(空扫描)，依照顺序先后发出锥束，其中第q个射线源在探测器单元p上得到光子强度I_pq，在放置被检测物体时，同时发出锥束，探测器单元p上得到光子强度y_p，其为各射线源发出射线穿过不同衰减路径之后到达p处光子强度之和。在上述扫描模式下，得到的y_p与目前单X线点源成像系统不同，除去散射和噪声等成分，y_p是一个由多点源产生的多条射线积分的总和，是一个模糊混叠的X线成像，而单点源系统的成像只与一条射线积分相关，是一个有相对清晰结构的图像。使用同时扫描方案可以只一次扫描快速的得到被检测物的X射线投影图像，极大的提升了成像的速度，避免了成像过程中被检测物体移动对成像质量的影响，但是由于图像的混叠模糊，需要进一步对投影图像进行图像恢复，在下文中将详细介绍。

下面结合图3、图4与图6详述将分时扫描方案获得的锥束投影恢复为理想平行光束投影的过程：

(c)X射线源阵列上的Q个射线源单元依次发出射线，或将射线源阵列上的Q个射线源单元分为锥束覆盖空间不重合的K个射线源组，结合对应空扫描数据得到被检测对象的线积分投影图，以上两种扫描模式都称为分时扫描，均假设为式(1)分时投影模型，

其中，q表示射线源的位置，f∈R^N表示待成像物体的线性衰减系数，a_pqn表示由射线源q和探测器p形成的X线光束与体素f_n之间交集的体积，n_pq表示线积分投影图上对应点处的噪声分量；同时，我们定义从点源q发出射线到探测器单元p的射线路径上，成像物体f对射线衰减的线性衰减系数总和为：

即函数f关于路径pq的Radon变换。

(d)由图3所示，虚拟X射线源阵列所在平面与实际射线源阵列平面重合，平行于相距为d的探测器平面，虚拟射线源阵列上存在与探测器单元数目P相同的射线源单元，各自发出垂直于对应探测器的强度为I′的笔形射线束，穿过成像物体f被衰减后只到达对应探测器单元。虚拟射线源p′发出的射线穿过物体，在对应探测器单元p处接收到的虚拟点源信号可表示为：

其中b_p′pn为笔形束射线p’p与体素f_n相交的体积。同时定义点源p′到探测器p形成的笔形束射线穿过成像物体，在穿过路径上对该虚拟射线进行衰减的线性衰减系数总和为：

由于位置p′和p是一一对应的，式(4)可表示为：

(c)定义从g′计算g_pq的线性变换矩阵w_pqi，为虚拟投影g′_i向路径pq投影g_pq的线性加权系数，于是有：

联合式(1)有：

建立优化目标函数，有：

(d)对目标函数使用如下方法求解：

初始g′₀为全零向量，k＝0，

设置适当的松弛因子ω与TV参数，迭代求解即可获得分时投影模型下的投影恢复图像。此外，使用图8所示的一种加速迭代求解目标函数的图像恢复方法也可以用于求解目标函数，其中加速收敛技巧的使用可以更快的实现投影图像的恢复。

下面结合图3、图5与图6详述将同时扫描方案获得的锥束投影恢复为理想平行光束投影的过程：

(e)X射线源阵列上的Q个射线源单元同时发出射线，探测器接收透过被检测物体的X射线信号，得到一幅光子数投影图，可建立式(11)同时投影物理模型：

其中q表示射线源的位置，f∈R^N表示待成像物体的线性衰减系数，a_pqn表示由射线源q和探测器p形成的X线光束与体素f_n之间交集的体积，I_pq表示射线源q对应于探测器p方向的入射光子强度，s_p与n_p分别表示位置为p的探测器接收到的散射分量和噪声分量；同时，与式(2)相同，定义g_pq为从点源q发出射线到探测器单元p的射线路径上，成像物体f对射线衰减的线性衰减系数总和。如图4所示，定义虚拟平板射线源和实际平板射线源的位置关系，同式(5)定义g′为笔形束射线穿过成像物体，在穿过路径上对该虚拟射线进行衰减的线性衰减系数总和。因而，结合式(11)，有：

建立优化目标函数，有：

式中，λ表示权重因子，表示散射信号的梯度，β表示正则化因子，R表示正则项，Y为投影信号，S为散射信号；

(f)不考虑散射信号，目标函数可改写为：

对在g′₀处进行一阶泰勒展开得：

用式(15)替换目标函数中对应部分得修正的目标函数，利用牛顿迭代算法，对其分别求一阶导和二阶导，同时加入TV正则，迭代求解目标函数，得到同时投影模型下的投影恢复图像。此外，使用图8所示的一种加速迭代求解目标函数的图像恢复方法也可以用于求解目标函数，其中加速收敛技巧的使用可以更快的实现投影图像恢复。

下面结合图3、图7详述一种虚拟射线源投影到实际阵列源射线投影的简化线性变换方法：

已知射线源单元q，探测器单元p，取p、q线的中点M，沿M作阵列射线源平面和探测器平面的垂线，对应射线路径p′p为虚拟射线源到对应探测器单元的笔形束。由公式(3)可得从点源q到探测器单元p的射线途径上成像物体f的线性衰减系数线积分g_pq，考虑pq间的欧氏距离为‖p-q‖，探测器平面与射线源阵列平面间距d，则相应虚拟射线路径对应的衰减系数线积分为：

g_pq可表示为：

使用虚拟射线源投影到实际阵列源射线投影的简化线性变换方法可以在误差允许的范围内，降低空间插值重排的运算量，提高计算速度，实现快速成像。

下面结合图8详述本发明提出的一种加速迭代求解目标函数的图像恢复方法：

令Ψ(x)表示目标函数，x为求解目标，t表示加速因子，z表示加速项，初始化重建目标及加速项x⁽⁰⁾＝z⁽⁰⁾，t₀＝1；

依照上述方法设定收敛条件，迭代求解目标函数，可以较快的收敛得到投影恢复图像。

使用以上所述的两种投影模型以及与对应的图像恢复方法，作为范例，选取一组特定参数进行试验:平板射线源上共Q＝9×9个射线源(相距d＝5.5mm)，使用射线源分时和同时扫描两种模式，每个锥束半锥角θ＝4.57°，射线源阵列与平板射线源间距100mm，中间放置的成像模体，其像素大小为50×256×256，每个像素0.25mm，探测器分辨率为P＝640×640，每个探测器单元长宽均为0.1mm，比较本发明提出的投影图像恢复方法与理想的平行束投影图像的相对误差，以此来评估算法的有效性。实验用模体及理想投影如图9a～9d所示。

定义重建图像的相对误差为：

其中，m,n为图像长宽，I为真实目标图像，I′为重建目标图像。

实验结果图像由图10、图11所示，数据分析如下表所示：

使用本发明提出的成像模式在分时投影与同时投影下基本可以达到令人满意的投影图像恢复结果，其中分时模式具有成像效果好的特点，而同时扫描模式具有成像速度快的特点。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。