一种基于深度神经网络的CT重建方法与流程

本发明属于计算成像技术领域，尤其涉及一种基于深度神经网络的ct重建方法。

背景技术：

计算机断层扫描(x-raycomputedtomography，简称ct)技术已经被广泛地用于医疗诊断、安全检测以及工业检测等领域。传统的ct需要围绕目标物体或样本，从所有探测角度进行完整的照射和测量，并通过数值算法，如滤波反投影算法(filteredback-projectionalgorithm，简称fbp)对目标物体或样本进行重建，这将不可避免地导致较高的辐射暴露、测量时间和检测成本。在医疗诊断或生物医学成像等领域，大量的辐射将会严重影响病人的健康或对生物样本组织产生不利影响，例如细胞和神经的损伤等。为了在兼顾图像重构质量的条件下，尽量降低辐射暴露，可以采用空间结构照明对物体进行照射，即在目标物体与x射线光源之间添加一个透射式的编码模板，用于阻挡部分x射线，形成结构照明。同时，为了降低辐射剂量、缩短检测时间，可以从少数几个角度对目标物体进行照射，获得非完备测量，而非进行密集视角下的完备测量。但是，如果大幅降低x射线的辐射剂量，往往会导致ct图像重建质量的下降。

除此之外，现有的ct重建方法运算量大，且耗时较长。另外，现有的x射线编码层析成像方法为了提高图像重构质量，需要对编码模板图案进行优化。但是现有的x射线编码模板优化方法的计算复杂度较高，且无法实现编码模板与物体重建之间的协同优化与设计，因此影响了ct图像重建的精度和速度。

综上所述，如何在降低x射线辐射剂量的同时，提高编码模板优化和ct图像重建过程的效率，同时进一步提高ct图像重建精度，是本领域一个亟待解决的关键问题。

技术实现要素：

为解决上述问题，本发明提供一种一种基于深度神经网络的ct重建方法，能够快速估计出目标物体或样本的影像，极大地提高运算效率。

一种基于深度神经网络的ct重建方法，包括以下步骤：

s1：基于ct成像系统分别获取b个目标物体的真实衰减系数向量fⁱ，且i＝1,2,…,b；同时，分别将各目标物体作为当前目标物体，执行衰减系数向量估计模型获取操作，得到各目标物体的衰减系数向量估计值其中，所述衰减系数向量估计模型获取操作具体为：

s11：构建ct成像系统的成像模型的稀疏表示如下：

其中，yp表示ct成像系统中第p个x射线光源经各自对应的编码模板调制后，再对当前目标物体进行扫描而在探测器上得到的成像模型，cp为第p个x射线光源对应的编码模板的模板矩阵对角化后的孔径矩阵，hp为第p个x射线光源对应的感知矩阵，ψ为稀疏基，ap＝cphpψ为对角化后的第p个x射线光源对应的稀疏感知矩阵，为当前目标物体的理论最优稀疏表示系数，p＝1,2,…,p，p为x射线光源的数量；

s12：获取当前目标物体的理论最优稀疏表示系数的估计值的迭代公式如下：

其中，为第n+1次迭代得到的估计值，为第n次迭代得到的估计值，σβ为非线性操作函数，w和z均为辅助变量，且w＝i-a^ta，z＝a^t，i为单位矩阵，t表示转置，

s13：将辅助变量w和z分别等效为第一卷积核φ与第二卷积核γ，将每一次迭代等效为一个卷积层，构建估计值基于深度神经网络的迭代公式如下：

其中，θj+1为第j+1层卷积层输出的稀疏表示系数矩阵，θj为第j层卷积层输出的稀疏表示系数矩阵，σβ为非线性操作函数，φj为第j层卷积层对应的第一卷积核，γj为第j层卷积层对应的第二卷积核，j＝1,2,...,k，k为深度神经网络的卷积层的设定数量；y为各x射线光源对应的成像模型y的矩阵形式；

s14：构建用于获取当前目标物体的衰减系数向量估计值的解码器，具体为：

其中，θk为θk转换而成的向量形式，且θk为第k层卷积层输出的稀疏表示系数矩阵；

s2：基于各目标物体的真实衰减系数向量fⁱ与衰减系数向量估计值构建深度神经网络的目标函数e：

其中，ε为设定的收敛因子，为第i个目标物体对应的深度神经网络第k层卷积层输出的稀疏表示系数向量；

s3：根据目标函数e优化变量αpl、φj、γj以及ψ，具体为：分别将变量αpl、φj以及γj以及ψ作为变量ω代入下式迭代maxloopp次，将第maxloopp次得到的变量αpl、φj以及γj以及ψ作为最优解：

其中，为目标函数e对变量ω的梯度，ωm为第m次迭代时得到的变量ω的值，ωm+1为第m+1次迭代时得到的变量ω的值，m＝1,2,…,maxloopp，step为设定步长；αpl为与孔径矩阵cp相关的辅助参数，且两者满足如下关系：

其中，cpl为第p个x射线光源所对应的孔径矩阵cp对角线上的第l个元素；

s4：根据步骤s3优化得到的αpl得到ct成像系统中的各x射线光源对应的编码模板的最优模板矩阵；同时，将优化得到的φj以及γj代入步骤s13中的迭代公式，得到最优深度神经网络；将优化得到的ψ代入步骤s14中的估计值的计算公式，得到最优解码器；

s5：采用物理方法获取各最优模板矩阵对应的最优编码模板，并将各最优编码模板作为ct成像系统中的各x射线光源对应的编码模板；

s6：采用装有最优编码模板的ct成像系统对目标物体进行层析摄影后，从探测器上得到成像量测值y*，然后将成像量测值y*输入最优深度神经网络中，得到目标物体对应的稀疏表示系数θ*，再将稀疏表示系数θ*输入最优解码器中，得到目标物体对应的衰减系数向量f*，并将衰减系数向量f*作为目标物体最终的ct重构图像。

进一步地，当前目标物体的理论最优稀疏表示系数的估计值的获取方法为：

根据ct成像系统的成像模型构建当前目标物体的稀疏表示系数优化函数：

其中，θ为当前目标物体的稀疏表示系数，为二范数的平方，||·||0为零范数，s表示稀疏表示系数中非零元素数量的最大值；

采用迭代硬阈值算法求解所述稀疏表示系数优化函数，得到当前目标物体的理论最优稀疏表示系数的估计值

有益效果：

本发明提供一种基于深度神经网络的ct重建方法，首先将传统的梯度迭代算法进行展开和截断，仅将前k次迭代的过程等效为深度神经网络的各层卷积层，而不是基于经验模型或迁移方法构建神经网络结构，因此本发明依靠物理模型就可以给出结构优良的网络结构，同时可以对网络参数进行初始化；

其次，本发明将x射线编码模板，稀疏基，以及网络中的其它参数放在同一个框架，即同一个迭代公式下进行协同优化，属于端到端的非监督式训练方法，无需对训练样本进行标记，因此可以极大地降低神经网络的训练成本和工作量，缩短训练时间，同时也可以避免标记错误所导致的训练误差或训练失败，由此提高了优化自由度，更有利于获得网络参数优化问题的全局最优解，能够有效提高了ct系统的成像精度；

最后，采用本发明优化得到的x射线编码孔径可以极大的减少暴露给病人的辐射剂量，使得ct系统更加安全可靠；同时，本发明通过训练过程对编码模板图案进行优化，可以在降低x射线辐射剂量的同时，提高ct图像的重构质量。

附图说明

图1为压缩x射线层析摄影合成系统示意图；

图2为本发明涉及的iht算法的迭代过程示意图；

图3为本发明中mcnn网络结构示意图和解码器示意图；

图4为用于训练mcnn网络的训练数据以及测试数据；

图5为采用mcnn优化后的编码模板对测试数据中的物体进行x射线层析摄影，再利用训练好的mcnn网络重构图像的结果；

图6为采用mcnn对x射线编码模板优化的结果及其与未优化的随机编码模板图案对比；

图7为采用传统方法优化后的编码模板对测试数据中的物体进行x射线层析摄影，再利用梯度投影稀疏重构(gradientprojectionforsparsereconstruction，简称gpsr)算法获得的物体重建图像；

图8为采用传统的利用梯度迭代方法获得的编码孔径优化结果及其与未优化的随机编码模板图案对比。

具体实施方式

为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。

本发明的原理：现有的基于梯度的ct重建方法需要通过大量的循环迭代过程，使得重建的ct图像结果尽可能接近于物体真实值。然而，基于梯度的ct重建方法运算量大，且耗时较长。本发明根据ct系统的物理成像模型，以iht算法为例，将传统的求解目标物体稀疏系数的迭代过程进行展开和截断，进而构建模型驱动的卷积神经网络(model-drivenconvolutionneuralnetwork，简称mcnn)网络，每一步迭代过程看作mcnn网络的一层结构。mcnn网络分为两个模块，分别记为“mcnn模块1”和“mcnn模块2”。mcnn模块1将编码模板视为mcnn网络中的二值神经网络层，它的输入为目标物体或样本，输出为x射线编码ct成像系统的压缩投影测量数据，即mcnn模块1模拟了ct系统的成像过程。mcnn模块2的输入为x射线编码ct成像系统的测量数据，输出为ct重建的物体在某个稀疏基下进行展开后的稀疏系数。综合两个模块后，mcnn整体网络的输入为目标物体或样本，输出为ct重建的物体在某个稀疏基下进行展开后的稀疏系数。另外，通过构建一种端到端的非监督式网络训练方法，避免了对训练样本的标记任务，可以极大地简化mcnn的训练过程，缩短训练时间，同时也可以避免标记错误所导致的训练误差或训练失败。该方法在mcnn之后添加解码器，解码器为一个全连接层，用于模拟物体图像的稀疏基。解码器的输入为mcnn网络的输出数据，即目标物体或样本的稀疏系数，解码器的输出为ct系统的重建图像。利用反向传播算法对首尾相连的mcnn网络和解码器进行优化，从而构建了一种可以协同优化x射线编码模板图案、物体图像稀疏基以及网络中跟其它各项参数的框架。很大程度上提高了信息的利用率，使得优化结果取得全局最优解，进而在降低辐射剂量的同时提高物体图像的重构质量和重建速度。

一种基于深度神经网络的ct重建方法，包括以下步骤：

s1：基于ct成像系统分别获取b个目标物体的真实衰减系数向量fⁱ，且i＝1,2,…,b；同时，分别将各目标物体作为当前目标物体，执行衰减系数向量估计模型获取操作，得到各目标物体的衰减系数向量估计值其中，所述衰减系数向量估计模型获取操作具体为：

s11：构建ct成像系统的成像模型的稀疏表示如下：

其中，yp表示ct成像系统中第p个x射线光源经各自对应的编码模板调制后，再对当前目标物体进行扫描而在探测器上得到的成像模型，cp为第p个x射线光源对应的编码模板的模板矩阵对角化后的孔径矩阵，hp为第p个x射线光源对应的感知矩阵，ψ为稀疏基，ap＝cphpψ为对角化后的第p个x射线光源对应的稀疏感知矩阵，为当前目标物体的理论最优稀疏表示系数，p＝1,2,…,p，p为x射线光源的数量。

需要说明的是，ct成像系统的成像模型的稀疏表示的推导过程如下：

s11a：如图1所示，本发明的ct成像系统为：将目标物体或样本栅格化，若物体为三维的，则栅格化为nx×ny×nz的图形，其中nx，ny和nz分别为物体沿着x，y和z轴的维度，若物体为二维的，则栅格化为nx×ny的图形，其中nx，ny分别为物体沿着x，y轴的维度；令表示为经过逐点扫描后得到的目标物体或样本的衰减系数向量。然后，将探测器上的测量值栅格化，设表示为探测器上测量值经过逐行扫描后得到的向量。

s11b：成像模型的离散形式可以简述为：

yp＝hpf

其中，hp为系统矩阵；

s11c：在目标物体与x射线光源间添加一个部分透光的编码模板，其中编码孔径的像素与探测器的像素一一对应，这样成像模型可以简述为：

yp＝cphpf

其中孔径矩阵cp的对角线元素为编码孔径的像素值，则对角化的孔径矩阵cp选出了hp中对应于透光像素的行。根据压缩感知理论，信号可以在某一组基(称为稀疏基)上进行稀疏表示，则上式可以等价写为：

上述公式即为mcnn模块1的表达式。mcnn模块1的输入为目标物体或样本，输出为x射线编码ct成像系统的测量数据，即mcnn模块1模拟了ct系统的成像过程。

s12：获取当前目标物体的理论最优稀疏表示系数的估计值的迭代公式如下：

其中，为第n+1次迭代得到的估计值，为第n次迭代得到的估计值，σβ为非线性操作函数，w和z均为辅助变量，且w＝i-a^ta，z＝a^t，i为单位矩阵，t表示转置；

需要说明的是，当前目标物体的理论最优稀疏表示系数的估计值的获取方法为：

根据ct成像系统的成像模型构建当前目标物体的稀疏表示系数优化函数：

其中，θ为当前目标物体的稀疏表示系数，为二范数的平方，||·||0为零范数，s表示稀疏表示系数中非零元素数量的最大值；

进一步地，估计值的迭代公式的推导过程为：

利用迭代算法可以求解稀疏表示系数优化函数中的优化问题，如图2所示，本发明以迭代硬阈值(iterativehardthresholding，简称iht)算法为例，说明mcnn网络的推导过程。需要注意的是，使用其它的迭代方法也可以用于求解稀疏表示系数优化函数中的问题，并由此推导出mcnn的网络结构。不同算法所推导出来的具体网络结构可能存在差异。令θ⁽ⁿ⁾表示稀疏系数在第n次迭代中的估计值，则第(n+1)次迭代中，变量θ由iht算法迭代更新为：

θ⁽ⁿ⁺¹⁾＝ωs[θ⁽ⁿ⁾-a^t(aθ⁽ⁿ⁾-y)]

其中ωs[·]为非线性的硬阈值函数。令w＝i-a^ta，z＝a^t，则经过化简后上式可表述为：

θ⁽ⁿ⁺¹⁾＝σβ[wθ⁽ⁿ⁾+zy]

s13：将辅助变量w和z分别等效为第一卷积核φ与第二卷积核γ，将每一次迭代等效为一个卷积层，构建估计值基于深度神经网络的迭代公式如下：

其中，θj+1为第j+1层卷积层输出的稀疏表示系数矩阵，θj为第j层卷积层输出的稀疏表示系数矩阵，φj为第j层卷积层对应的第一卷积核，γj为第j层卷积层对应的第二卷积核，j＝1,2,...,k，k为深度神经网络的卷积层的设定数量；y为各x射线光源对应的成像模型y的矩阵形式；

需要说明的是，本发明对步骤s12中的迭代过程进行展开，并截取前k步迭代，每一步迭代过程可以视为mcnn网络中的一层结构，mcnn网络中的上一层输出作为下一层的输入。非线性函数作为mcnn网络的激励层，线性函数可以转化为卷积操作，等效为mcnn网络中的卷积层。

mcnn模块2的输入为x射线编码ct成像系统的测量数据，输出为ct重建的物体在某个稀疏基下进行展开后的稀疏系数。通过上述转换，mcnn模块2网络每层的迭代可以转化为估计值基于深度神经网络的迭代公式如上。

s14：构建用于获取当前目标物体的衰减系数向量估计值的解码器，具体为：

其中，θk为θk转换而成的向量形式，且θk为第k层卷积层输出的稀疏表示系数矩阵；

也就是说，如图3所示，将mcnn网络的最后一层输出θk转化为其向量表示形式θk，在其后方接入解码器，利用mcnn网络提取出的目标物体或样本的稀疏系数，通过解码运算，快速估计出待测物体的ct重构图像。解码器的输出为

其中稀疏基ψ被当作是解码器中的全连接层，通过mcnn网络和解码器相连接可以协同优化编码模板、物体图像的稀疏基，以及网络中的其它各项参数。

此外，需要说明的是，基于ct成像系统分别获取b个目标物体的真实衰减系数向量fⁱ可以是采用传统方法模拟仿真得到，也可以是根据传统方法采用未经过优化的ct成像系统对b个目标物体实际成像得到。

s2：基于各目标物体的真实衰减系数向量fⁱ与衰减系数向量估计值构建深度神经网络的目标函数e：

其中，fⁱ为第i个目标物体的真实衰减系数向量，为第i个目标物体的衰减系数向量估计值，ε为设定的收敛因子，为第i个目标物体的稀疏表示系数向量；

需要说明的是，ε是为了使网络加速训练，通过在线寻找等方式得到的参数，使得网络训练过程中的初始误差最小，有利于训练过程的快速收敛。

s3：根据目标函数e优化变量αpl、φj、γj以及ψ，具体为：分别将变量αpl、φj以及γj以及ψ作为变量ω代入下式迭代maxloopp次，将第maxloopp次得到的变量αpl、φj以及γj以及ψ作为最优解：

其中，为目标函数e对变量ω的梯度，ωm为第m次迭代时得到的变量ω的值，ωm+1为第m+1次迭代时得到的变量ω的值，m＝1,2,…,maxloopp，step为设定步长；αpl为与孔径矩阵cp相关的辅助参数，且两者满足如下关系：

其中，cpl为第p个x射线光源所对应的孔径矩阵cp对角线上的第l个元素。

需要说明的是，由于x射线编码模板的像素与探测器像素一一对应，因此其可以被当作一个点到点的二值神经网络层，它的网络权重被限制为0或1，0表示对应的编码模板的像素阻挡x射线，1表示对应像素透过x射线。为了在后续的网络训练过程中对编码模板进行优化，需要将二值的编码模板像素连续化，例如采用下列函数：

也就是说，本发明采用αpl表示孔径矩阵对角线元素连续化后的辅助参数。在网络训练过程中，利用基于梯度的算法优化参数αpl，然后让其通过硬阈值函数操作，就可以根据上式计算出二值的编码模板图案。

由此可见，本发明采用一种端到端的非监督训练方法，来联合优化物体图像稀疏基、x射线编码模板图案，以及mcnn网络中的其它各项参数，为了加速训练过程，本发明采用批量训练方式，每次网络的前向和反向过程都有b个数据参与训练。然后，基于每一批训练样本，采用反向传播算法，对mcnn中的各项参数进行maxloopp轮迭代更新。

s4：根据步骤s3优化得到的αpl得到ct成像系统中的各x射线光源对应的编码模板的最优模板矩阵；同时，将优化得到的φj以及γj代入步骤s13中的迭代公式，得到最优深度神经网络；将优化得到的ψ代入步骤s14中的估计值的计算公式，得到最优解码器；

s5：采用物理方法获取各最优模板矩阵对应的最优编码模板，并将各最优编码模板作为ct成像系统中的各x射线光源对应的编码模板；

s6：采用装有最优编码模板的ct成像系统对目标物体进行层析摄影后，从探测器上得到成像量测值y*，然后将成像量测值y*输入最优深度神经网络中，得到目标物体对应的稀疏表示系数θ*，再将稀疏表示系数θ*输入最优解码器中，得到目标物体对应的衰减系数向量f*，并将衰减系数向量f*作为目标物体最终的ct重构图像。

也就是说，本发明网络训练结束后，将优化后得到的二值编码模板图案提取出来，用物理方法制造优化所得的编码模板的硬件设备，并安装于ct系统中作为真实实验系统中的编码模板。将mcnn模块1与mcnn模块2进行分离。当需要计算某个目标物体的重构图像时，将目标物体放在ct系统的物体空间。通过x射线层析摄影后，获得测量数据。将得到的测量数据输入到训练好的mcnn模块2中，可得到重构物体图像的稀疏系数。再利用优化后的稀疏基进行解码运算，即可快速重构出目标物体的图像。

下面结合附图对本发明的ct重建效果进行说明。

图4(a)为用于训练mcnn的训练数据集，(b)为测试数据集。图5(a)-(d)、(e)-(h)、(i)-(l)、(m)-(p)分别为采用mcnn优化后的编码模板对测试数据中的四个物体进行x射线层析摄影，再利用训练好的mcnn网络重构物体图像的结果，利用峰值信噪比(peaksignal-to-noiseratio，简称psnr)来评估图像质量，其平均psnr值分别为54.77，44.28，54.86，51.71。图6(d)为初始的编码模板图案，(e)为采用mcnn优化后的编码模板图案，(f)为编码模板的变化图案。图7(a)-(d)、(e)-(h)、(i)-(l)、(m)-(p)分别为采用传统方法优化后的编码模板对测试数据中的四个物体进行x射线层析摄影，再利用gpsr算法重构得到的物体图像，其平均psnr值分别为36.41，37.10，36.41，36.75。图8(d)为初始的编码模板图案，(e)为采用传统的基于梯度迭代方法优化后的编码模板图案，(f)为编码模板的变化图案。对比图5和图7可知，采用本发明中的mcnn方法能够快速计算出目标物体或样本的ct重建结果，并有效提高了物体图像的重建质量。

综上所述，本发明的目的是提供一种基于深度神经网络的快速ct重建方法。该方法使用mcnn网络对ct图像进行重建。首先将传统的梯度迭代算法进行展开和截断，进而构建mcnn网络结构。mcnn网络分为两个模块，分别记为“mcnn模块1”和“mcnn模块2”。mcnn模块1的输入为目标物体或样本，输出为x射线编码ct成像系统的测量数据，即mcnn模块1模拟了ct系统的成像过程。mcnn模块2的输入为x射线编码ct成像系统的测量数据，输出为ct重建的物体在某个稀疏基下进行展开后的稀疏系数。综合两个模块后，mcnn整体网络的输入为目标物体或样本，输出为ct重建的物体在某个稀疏基下进行展开后的稀疏系数。与传统的基于经验或迁移学习构建神经网络的方法不同，mcnn方法通过物理模型推导出网络结构，可以给出网络结构的系统性设计方法，获得结构优良的神经网络，并可以实现网络参数的系统性初始化。

其次，本发明将x射线编码模板视为mcnn网络中的一层结构，与网络的其它参数进行协同优化。由于编码模板的像素与探测器的像素一一对应，因此编码模板可被当作一个点到点的二值神经网络层，它的网络权重被限制为0或1，0表示对应的编码孔径像素阻挡x射线，而1表示对应像素透过x射线。

另外，本发明提供了一种端到端的非监督式网络训练方法，避免了对大量训练样本进行标记的任务。本发明的训练方法将mcnn网络看作对目标物体或样本的编码过程，mcnn整体网络的输入为目标物体或样本图像，输出为目标物体或样本在某组稀疏基上展开的稀疏系数。在mcnn之后添加解码器，解码器为一个全连接层，用于模拟物体图像的稀疏基。解码器的输入为mcnn网络的输出数据，即目标物体或样本的稀疏系数，解码器的输出为物体或样本的重建图像。在端到端的无监督训练过程中，将mcnn网络与解码器首尾相连，采用反向传播算法对编码模板图案、mcnn网络参数和解码器参数进行联合优化，使得mcnn网络的输入数据与解码器的输出数据之间的误差最小化。通过这种设计，使得x射线编码模板，物体图像稀疏基以及网络中的其它各项参数在同一个框架下进行协同优化，很大程度上提高了信息的利用率，使得优化结果取得全局最优解，进而提高物体图像的重构质量。

当网络训练结束后，可用物理方法制造优化所得的编码模板的硬件设备，并安装于ct系统中。对于任意的待测物体，通过带有编码模板的ct系统，可以获得压缩投影测量数据。将mcnn模块1与mcnn模块2进行分离，将投影测量数据输入mcnn模块2，可得到重构物体图像的稀疏系数，再利用优化后的稀疏基进行解码运算，即可快速重构出待测物体的图像。相比现有的迭代式ct重构算法，本发明中的方法能够明显降低计算复杂度，缩短x射线编码模板优化和ct重建图像所需的计算时间。另外，本发明提出的端到端非监督式训练方法，无需对训练样本进行标记，因此可以极大地简化mcnn的训练过程，缩短训练时间，也可以避免标记错误所导致的训练误差或训练失败。同时，该训练方法还可以同时优化网络参数和稀疏基，提高优化自由度，有利于进一步提高物体图像的重构精度。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。