CN110415307A - 一种基于张量补全的多能ct成像方法、装置及其存储设备 - Google Patents

Abstract

本发明公开了一种基于张量补全的多能CT成像方法、装置及其存储设备，该方法及系统，首先需使用FDK算法，对得到的所述每段窄束能谱的投影值，分别进行处理后，得到每个能量段的重建图像；然后，对所得的每个能量段的重建图像建模为3阶张量，其中，建立张量核范数和全变分正则化的最小化模型，提高各能量段重建图像的精度；最后，将建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。实施本发明的有益效果是，通过基于GATE的多能CT模拟系统采集数据，通过将CT问题固有的多维性质与张量结合起来，更加精确的重建了CT扫描的图像。

Description

一种基于张量补全的多能CT成像方法、装置及其存储设备

技术领域

本发明涉及断层影像(CT)成像系统以及成像方法，更具体地说，涉及一种基于张量补全的多能CT成像方法及系统。

背景技术

CT成像基本原理是用X线束对人体检查部位一定厚度的层面进行扫描，由探测器接收透过该层面的X线，转变为可见光后，由光电转换器转变为电信号，再经模拟/数字转换器(analog/digital converter)转为数字信号，输入计算机处理。

为了进一步的深入研究，人们运用虚拟仿真技术，借助模拟软件，如：Geant4、SimSet等实现CT成像系统的仿真模拟。当前基于光子计数探测器的CT成像(PhotonCounting Computed Tomography，PCCT)仿真技术主要包括解析仿真算法和蒙特卡洛仿真方法。上述技术的研究均基于GATE(Geant4Application for Tomographic Emission)仿真软件，此仿真软件具有准确的物理模型，友好的用户界面，灵活的使用方法等特点。本发明根据光子计数探测器的能谱滤波分离成像原理，搭建了基于GATE的多能CT虚拟平台。

对于多能谱CT问题，最近，在基于低秩、稀疏分解或者低秩假设的张量或者矩阵的恢复已经做了大量的工作。利用CT问题固有的张量性质，可以使用更广泛的工具来分析这些结构。如今，张量分解工具如CP分解、Tucker分解已经用于化学分析、神经科学、计算机视觉和数据挖掘。高阶奇异值分解(High Order Singular Value Decomposition，HOSVD)分解也可以看作多维奇异值分解，已经应用于图像处理应用，如面部识别。尽管这些工具处理多维数据很有效，但是要找到这些解，需要解决一个复杂的非凸优化问题，该问题同时具有较差的收敛性。此外，对于CP和Tucker方法，需要已知成分个数，因此考虑使用一种替代方法进行进一步的分析。

发明内容

本发明要解决的技术问题在于，针对现有技术的缺陷，提供一种基于张量补全的多能CT成像方法及系统。

本发明解决其技术问题所采用的技术方案是：构造一种基于张量补全的多能CT成像方法，具有如下步骤：

S1、构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

S2、确定一定电压下的连续X射线能谱，并对所述连续X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到步骤S1构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

S3、使用FDK算法，对步骤S2得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

S4、将步骤S3所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

S5、在能量维度下，将步骤S4建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。

进一步的，步骤S1中，所述光源类型为锥束，其中，锥束角度为6.8°，使用gamma粒子，每次投影粒子数为2×10⁶；探测器为50*50mm²的面探测器，由200×200个0.25×0.25×1mm³的硅探测元组成。

进一步的，步骤S4中，基于张量核范数和全变分正则化建立的最小化模型为：

s.t.Z＝bcirc(χ)

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||_*表示求解“·”的矩阵核范数。

进一步的，步骤S5中，在使用加权融合算法进行优化加权之前，首先需提取各能量段重建图像的噪声和对比噪声比信息，所述噪声水平σ和对比噪声比CNR的数学表达式分别为：

其中，j＝1,2,…,M,M为像素个数；和为第j个能量段重建图像的像素值及其均值，x_c和x_b分别为第一、第二衰减系数材料对应的像素值。

进一步的，在提取出各能量段重建图像的噪声和对比噪声比后，利用对比噪声方差比来表示每个能量段的加权权重因子w_n：

其中，k为能量段的总数，表示第n个能量段的噪声方差；C_n＝|x_c,n-x_b,n|为第n个能量段重建图像中，第一、第二衰减系数材料的像素差异值。

进一步的，在步骤S5中，根据所得的每个能量段的加权权重因子w_n，利用加权融合算法对各能量段的重建图像进行优化加权，得到最终的成像图像x：

其中，k为能量段总数，ω_k为第k个能量段的重建图像的权重，x^k为第k个能量段重建的图像，x为各能量段加权融合后形成的最终成像图像。

进一步的，针对步骤S3得到的各个能量段的重建图像，使用加权融合算法对所得的每个重建图像进行优化加权，获得对应的第一全能谱图；将所述第一全能谱图与步骤S5中得到的最终成像图像进行对比，进一步判断当前张量补全的算法的计算精度。

一种存储设备，所述存储设备存储指令及数据用于实现上述任意一项提出的一种多能CT成像方法。

一种基于张量补全的多能CT成像装置，具有以下模块：

多能CT模拟系统构建模块，用于构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

X射线能谱抽样模块，用于确定一定电压下的连续X射线能谱，并对所述连续的X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到多能CT模拟系统构建模块构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

重建图像构建模块，用于使用FDK算法，对X射线能谱抽样模块得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

张量建模模块，用于将重建图像构建模块所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

图像成像模块，用于在能量维度下，将张量建模模块建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。

进一步的，张量建模模块中，基于张量核范数和全变分正则化建立的最小化模型为：

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||_*表示求解“·”的矩阵核范数。

实施本发明的一种基于张量补全的多能CT成像方法及系统，通过基于GATE的多能CT模拟系统采集数据，通过将CT问题固有的多维性质与张量结合起来，更加精确的重建了CT扫描的图像。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明公开的多能CT成像方法流程图；

图2是本发明公开的多能CT成像系统结构图；

图3为SpectrumGUI生成的0-125keV能谱；

图4为投影过程仿真模型；

图5为模体横截面模型；

图6为模体参数设置；

图7为康普顿散射过程；

图8为3阶张量展开示意图；

图9为张量纤维示意图；

图10为张量切片示意图；

图11为FDK重建图像直接加权融合形成的最终图像；

图12为FDK重建图像建模为3阶张量后，经过张量补全运算后加权融合形成的最终图像。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

请参考图1，其为本发明公开的多能CT成像方法流程图，一种基于张量补全的多能CT成像方法，具有如下步骤：

S1、构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

本实施例中，通过搭建投影过程仿真模型，模拟得到每段窄束能谱的投影值，请参考图4，其为投影过程仿真模型，整个模型可以分为三个部分：

1、光源模块即x射线源部分，在此部分设置粒子类型为gamma、每次投影发射粒子数量2×10⁶、x射线能谱以及光源类型为锥形束，锥束角度为6.8°；

2、模体，参考图5所示的模体横截面示意图设置模体，整个模体包含5个圆柱体，材料分别为：水、铝、玻璃、尖晶石以及PVC，具体参数请参考图6；

3、探测器，所述探测器为50×50mm²的阵面探测器，由200×200个0.25×0.25×1mm³的硅探测元组成；

此外，当x射线穿透物体时，除了吸收和反应作用，部分光子还会发生散射过程，如图7所示；其中，工业CT能量范围内散射的主要成分是康普顿散射，因此，为了更加真实地模拟CT成像系统，在模拟过程中加入了康普顿散射过程。

本实施例中，采用狄拉克电子理论，得到康普顿散射微分截面如下：

式中：r₀为经典电子半径，α为以静止电子能量为单位的入射光子能量，θ为散射角，dσ/dΩ为单位角的散射微分截面。

由于在GATE的模拟系统中，Physical模块主要是用来设置粒子运输过程中与物质发生的相互作用和一些物理过程。这其中包括有：核素的衰减、电磁场作用、康普顿效应、光电效应、光的反射和折射以及粒子的追踪传输等物理过程。同时，GATE的模拟系统中对于电磁场作用过程中提供了两个模拟数据包，分别为：标准能量电磁过程(StandardEnergyElectromagnetic Process，SPE)和低能量电磁过程(Low EnergyElectromagnetic Process，LPE)，这两个数据包是由大量实验数据所组成。因此，可以根据模拟所得的数据与实际采集的数据进行很好的比较。其中，SEP所处理的粒子是能量在10keV以上所发生的光电效应和康普顿效应；而LEP则用于模拟粒子在10keV以下的电磁过程。

S2、确定125keV电压下的连续X射线能谱，并对所述连续X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到步骤S1构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

S3、取各能量段的投影数据，使用FDK算法，对步骤S2得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

FDK算法实际上是二维扇形滤波反投影算法的三维扩展，它包括投影数据的预加权、一维滤波和反投影几个步骤，重建的具体过程为：

1、对投影数据加权：

物体上待重建点位f(x,y,z)的射线源在探测器平面坐标为(Y,Z)，p_β(Y,Z)表示物体点的投影数据，R为射线源S绕中心轴Z的旋转半径。修正的权系数的几何意义是任一射线与中心射线之间夹角的余弦。

2、对二维投影数据，进行滤波，其滤波函数为h(Y)；

3、沿射线方向进行反投影重建，获得密度图像f(x,y,z)：

式中，U为加权因子，β为射线源与x轴正方向的夹角(即旋转角度)，dβ为在中心平面上的扇束的投影角增量。

S4、将步骤S3所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

由于在利使用量补全的算法前，首先需要要定义张量(用χ表示)，具体的操作步骤为：

1、从张量中抽取向量，具体为在矩阵中固定其中任意一个维度，可以得到对应的行或者列。类似于矩阵操作，固定其他维度，只保留一个维度变化(如图9所示)。

2、在张量中保留两个维度变化，将其他维度固定之后，可以得到一个矩阵，所得的这个矩阵即为张量的切片(如图10所示)；其中，用X^(k)表示χ的第k个前面切片。

3、将一个k阶张量的元素重新排列(χ_(l)，即对张量mode-l的纤维重新排列)，得到一个矩阵的过程(如图10所示)，张量元素展开为矩阵元素(i_l，j)时，元素对应位置为：

用x表示张量χ的向量化，x_l表示张量χ_(l)的向量化，则X和X_(l)的关系可以表示为：

其中，是转置矩阵，由于x和x₁是相等的，则P₁为单位矩阵。

4、计算fold、unfold以及bcirc，其中unfold为张量前面切片的重排列为：

fold是unfold的逆操作为：fold(unfold(χ))＝χ；X^(k)为张量χ的第k个前片。

5、使用X^(k)组成块循环矩阵有如下定义：

由于bcirc(·)是一个线性操作，可以写作一系列转置矩阵，其中，利用x_c表示bcirc(χ)的向量形式，则对应的有：

其中x为张量χ的向量化，P_c为转置矩阵，P_c，i为P_c的子矩阵，P_c，i根据bcirc(χ)的第一列分块矩阵重排列x。

6、根据块循环矩阵的块对角性质，对任意块循环矩阵有如下性质：

其中，和为单位矩阵，是归一化的离散傅里叶变换矩阵,表示张量χ的mode-3纤维经过快速傅里叶变换后的张量的第n个前片。

基于上述张量操作的定义，由于矩阵核范数的单一不变性：

基于张量核范数和全变分正则化建立其最小化模型为：

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||_*表示求解“·”的矩阵核范数。

针对上述最小化模型，本实施例中使用的是二维离散TV：

其中： 表示张量χ的第k个前片的第i行第j列元素。

本实施例中，定义最小化模型的增广拉格朗日形式为：

其中,Y是对偶变量，η＞0作为惩罚项，＜.＞为内积，K₁，K₂∈R^M×N，定义如下：

利用ADMM算法求解该模型得到：

Yⁿ⁺¹：＝Yⁿ+η(χⁿ⁺¹-Zⁿ⁺¹).；

其中，使用循环块矩阵转置操作，可以将其写作：

利用前面讲到的转置操作，当前χ的更新可以写作：

对于Z的更新，进一步可以写作：

按照上述操作，对建模成的张量使用本发明提出的最小化模型优化后，在能量维度优化加权融合形成的重建图像如图12所示。

S5、在能量维度下，将步骤S4建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像；其中：

在使用加权融合算法进行优化加权之前，首先需提取各能量段重建图像的噪声和对比噪声比信息，所述噪声水平σ和对比噪声比CNR的数学表达式分别为：

其中，j＝1,2,…,M,M为像素个数；和为第j个能量段重建图像的像素值及其均值，x_c和x_b分别为第一、第二衰减系数材料的像素值；

多能谱重建图像的对比噪声比可由各个能量段的图像对比噪声比表示，即：

式中，为第k能量段的噪声方差，C_k为第k能量段图像中两种种材料间平均衰减系数的绝对差值，本发明中使用的是旋转铝和水。

对上述式子中第n(n＝1，2，3，…，K)个能量段图像权重w_n求其偏导数，可得到对比噪声比最大时的归一化权重w_n为：

其中，k为能量段的总数，表示第n个能量段的噪声方差；C_n＝|x_c,n-x_b,n|为第n个能量段重建图像中，第一、第二衰减系数材料的像素差异值；

根据所得的每个能量段的加权权重因子w_n，利用加权融合算法对各能量段的重建图像进行优化加权，得到最终的成像图像x：

其中，k为能量段总数，ω_k为第k个能量段的重建图像的权重，x^k为第k个能量段重建的图像，x为各能量段加权融合后形成的最终成像图像。

将各能量段FDK算法重建图像，按照上述步骤优化加权融合，形成的图像效果如图11所示。

请参考图2，其为本发明公开的多能CT成像系统结构图，该系统具体包括以下模块：

多能CT模拟系统构建模块L1用于构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

X射线能谱抽样模块L2用于确定一定电压下的连续X射线能谱，并对所述连续的X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到多能CT模拟系统构建模块构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

重建图像构建模块L3用于使用FDK算法，对X射线能谱抽样模块得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

张量建模模块L4用于将重建图像构建模块所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

图像成像模块L5用于在能量维度下，将张量建模模块建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。

本实施例中，所述多能CT成像装置的张量建模模块中，基于张量核范数和全变分正则化建立的最小化模型为：

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||_*表示求解“·”的矩阵核范数。

作为一个优选的实施例，本发明公开了一种存储设备，所述存储设备用于存储指令及数据用于实现上述任一一种多能CT成像方法实现最终的图像成像。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (10)

1.一种基于张量补全的多能CT成像方法，其特征在于，具有如下步骤：

S1、构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

S2、确定一定电压下的连续X射线能谱，并对所述连续X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到步骤S1构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

S3、使用FDK算法，对步骤S2得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

S4、将步骤S3所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

S5、在能量维度下，将步骤S4建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。

2.根据权利要求1所述的多能CT成像方法，其特征在于，步骤S1中，所述光源类型为锥束，其中，锥束角度为6.8°，使用gamma粒子，每次投影粒子数为2×106^；探测器为50*50mm²的面探测器，由200×200个0.25×0.25×1mm³的硅探测元组成。

3.根据权利要求1所述的多能CT成像方法，其特征在于，步骤S4中，基于张量核范数和全变分正则化建立的最小化模型为：

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||*表示求解“·”的矩阵核范数。

4.根据权利要求1所述的多能CT成像方法，其特征在于，步骤S5中，在使用加权融合算法进行优化加权之前，首先需提取各能量段重建图像的噪声和对比噪声比信息，所述噪声水平σ和对比噪声比CNR的数学表达式分别为：

其中，j＝1,2,…,M,M为像素个数；和为第j个能量段重建图像的像素值及其均值，x_c和x_b分别为第一、第二衰减系数材料对应的像素值。

5.根据权利要求4所述的多能CT成像方法，其特征在于，在提取出各能量段重建图像的噪声和对比噪声比后，利用对比噪声方差比来表示每个能量段的加权权重因子w_n：

其中，k为能量段的总数，表示第n个能量段的噪声方差；C_n＝|x_c,n-x_b,n|为第n个能量段重建图像中，第一、第二衰减系数材料的像素差异值。

6.根据权利要求5所述的多能CT成像方法，其特征在于，在步骤S5中，根据所得的每个能量段的加权权重因子w_n，利用加权融合算法对各能量段的重建图像进行优化加权，得到最终的成像图像x：

其中，k为能量段总数，ω_k为第k个能量段的重建图像的权重，x^k为第k个能量段重建的图像，x为各能量段加权融合后形成的最终成像图像。

7.根据权利要求1所述的多能CT成像方法，其特征在于，针对步骤S3得到的各个能量段的重建图像，使用加权融合算法对所得的每个重建图像进行优化加权，获得对应的第一全能谱图；将所述第一全能谱图与步骤S5中得到的最终成像图像进行对比，进一步判断当前张量补全的算法的计算精度。

8.一种存储设备，其特征包括：所述存储设备存储指令及数据用于实现权利要求1～7所述的一种多能CT成像方法。

9.一种基于张量补全的多能CT成像装置，其特征在于，具有以下模块：

多能CT模拟系统构建模块，用于构建基于GATE的多能CT模拟系统；所述多能CT模拟系统包括光源模块、模体和探测器；所述光源模块用于发射X射线到模体，其中，发射到模体上的X射线将进一步穿过模体，通过探测器对其进行接收并探测当前X射线的强度，通过所述X射线的强度生成对应的投影值；

X射线能谱抽样模块，用于确定一定电压下的连续X射线能谱，并对所述连续的X射线谱抽样得到若干段相邻的窄束能谱，其中，将每段窄束能谱输入到多能CT模拟系统构建模块构建的多能CT模拟系统中，通过所述探测器模拟得到每段窄束能谱对应的投影值；

重建图像构建模块，用于使用FDK算法，对X射线能谱抽样模块得到的所述每段窄束能谱的投影值，分别进行预加权、一维滤波和反投影处理后，得到每个能量段的重建图像；

张量建模模块，用于将重建图像构建模块所得的每个能量段的重建图像建模为3阶张量，建立张量核范数和全变分正则化的最小化模型，其中，所述最小化模型用于修复在投影过程中由于存在外界因素的干扰而引起失真的数据，提高各能量段重建图像的精度；

图像成像模块，用于在能量维度下，将张量建模模块建模所得张量中的每一个切片按加权融合算法进行优化加权，得到最终的成像图像。

10.根据权利要求9所述的多能CT成像装置，其特征在于，张量建模模块中，基于张量核范数和全变分正则化建立的最小化模型为：

其中，χ为各能量段重建图像构建的张量，bcirc(χ)为张量χ的块循环矩阵，TV(χ)为对张量χ进行全变分正则化；||·||_*表示求解“·”的矩阵核范数。