CN104346820A - 一种x光双能ct重建方法 - Google Patents

Abstract

本发明涉及一种基于自先验信息的X光双能CT重建方法，能够利用数据内在的信息提供先验模型，从而获得高质量的重建图像。本发明的X光双能CT重建方法具有如下步骤：（a）进行能谱标定和双能查找表建立；（b）利用双能CT成像系统的探测器采集得到双能CT成像系统的高能数据 p _H以及低能数据 p _L；（c）根据所得到的所述高能数据 p _H以及所述低能数据 p _L得到比例图像 r ₁以及 r ₂的投影图像 R ₁以及 R ₂；（d）利用第一局部平滑限制条件重建比例图像 r ₂，由此，得到电子密度图；（e）利用第二局部平滑限制条件重建比例图像 r ₁，由此，得到等效原子序数图像。在本发明中，有效地利用数据固有信息，能够在保持分辨率的情况下有效地抑制双能重建图像的噪声。

Description

一种X光双能CT重建方法

技术领域

本发明涉及断层影像（CT）重建的方法，特别涉及一种基于自先验信息的X光双能CT重建方法。

背景技术

断层影像(CT)的图像对比度与扫描所用的X光源能谱分布有很大关系。传统CT使用具有能谱分布的一个光源进行成像，有时会出现信息模糊致使两种不同材料在CT图像上完全相同，而双能CT使用两个不同分布的能谱对物体进行成像，能够消除单能谱情况下的信息模糊。双能X光断层成像技术（CT）利用物质在不同能量下的衰减差异能够获得关于物体的多重物理特性参数的分布信息，例如，电子密度分布、等效原子序数分布、多个能量下的单能衰减图像。这样，双能X光CT能够被用来进行传统CT的射线硬化校正、临床高对比度能谱图像获取、工业和安检的特殊和危险物品检测等。与传统的X光CT成像技术相比，双能CT在其成像功能上的突破对医学诊断技术、无损检测、以及安全检查等应用都具有重大意义，因此，在近年来得到越来越广泛的关注。此外，双能X光CT的重建方法是当前的一个研究热点。

就目前来说，双能CT重建主要有如下三种途径。（1）后处理方法：从低能数据和高能数据分别重建衰减系数分布图像，然后，对这两个衰减系数分布图像进行综合计算，由此，能够获得单能图像或者与能量无关的物理量（例如，原子序数、电子密度）的分布图像；（2）前处理方法：从低能数据和高能数据来解析能量相关和能量无关信号（也就是所谓的双能分解），解析后的信号属于投影域，然后，利用传统的CT重建方法对解析后的信号进行重建；（3）综合迭代法；对低能数据和高能数据直接使用迭代方法进行重建。目前，被广泛采用的是前处理方法，这是因为前处理方法比后处理方法更加准确，能够更好地去除X光宽谱的影响，另一方面，前处理方法相对于综合迭代法计算量较少。

目前，关于双能分解，有基材料分解和双效应分解这两种分解方式（例如，参照非专利文献1）。但是，在双效应分解方法中，通常等效原子序数重建图像的信噪比差，相对于此，电子密度图像具有比较好的信噪比（例如，参照非专利文献2）。此外，在双能CT重建中，能够利用蒙特卡罗（MonteCarlo）方法或者实验方法估计双能CT系统的能谱数据，此外，也可以建立查找表（例如，参照非专利文献2、3）。

现有技术文献

非专利文献1：Y. Xing， L. Zhang， X. Duan， J. Cheng， Z. Chen， "A reconstruction method for dual high-energy CT with Mev X-rays，" IEEE Trans Nucl. Sci. vol. 58， no. 2， pp537-546， 2011.；

非专利文献2：张国伟，双能X射线成像算法及应用研究[D]，北京：清华大学工程物理系，2008；

非专利文献3：G. Zhang， Z. Chen， L. Zhang， and J. Cheng， Exact Reconstruction for Dual Energy Computed Tomography Using an H-L Curve Method， 2006 IEEE Nuclear Science Symposium Conference Record， pp.  M14-462， 2006。

此外，在双能CT重建中，两个主要物理特性参数是等效原子序数和电子密度，由于在双能分解的过程中带有很强的不平衡性，因此，导致双能CT重建图像的噪声被放大，特别是等效原子序数分布的噪声被放大。

发明内容

本发明是为了解决上述课题而提出的，其目的在于提供一种双能CT重建方法，能够利用数据内在的信息提供先验模型，从而获得高质量的重建图像。

本发明提供一种X光双能CT重建方法，其特征在于，具有如下步骤：

（a）利用双能CT成像系统的探测器采集得到双能CT成像系统的高能数据p_H以及低能数据p_L；

（b）根据所得到的所述高能数据p_H以及所述低能数据p_L得到比例图像r₁以及r₂的投影图像R₁以及R₂；

（c）利用第一局部平滑限制条件重建比例图像r₂，并且，得到分解系数a₁；

（d）利用第二局部平滑限制条件重建比例图像r₁，并且，得到分解系数a₂。

此外，在本发明的X光双能CT重建方法中，其特征在于，

比例图像r₁以及r₂以如下的式（4）定义，

           （4），

投影图像R₁以及R₂以如下的式（6）定义，

                   （6），

其中，，，H是投影矩阵，a₁和a₂为分解系数，ε为一个小值常系数向量，ω₁和ω₂为能够任意选择的向量，

在所述步骤（c）中，以作为第一局部平滑限制条件的下式（9）重建r₂，

 s.t.                  (9)，

在所述步骤（d）中，以作为第二局部平滑限制条件的下式（8）重建r₁，

 s.t.                     (8)。

此外，在本发明的X光双能CT重建方法中，其特征在于，

在所述步骤（c）中，根据所述高能数据p_H重建高能情况下的有效线衰减系数μ_H，并且，选择ω₂ =μ_H，以作为第一局部平滑限制条件的式（9）重建r₂。

此外，在本发明的X光双能CT重建方法中，其特征在于，

在上述步骤（d）中，取ω₁ = a₂并且以作为第二局部平滑限制条件的式（8）重建r₁。

此外，在本发明的X光双能CT重建方法中，其特征在于，

还具有如下步骤：利用双效应分解，根据a₂＝r₂×diag（ω₂+ε）以及ρ_e＝2a₂得到电子密度图ρ_e。

此外，在本发明的X光双能CT重建方法中，其特征在于，

还具有如下步骤：利用双效应分解，根据来得到等效原子序数图Z^eff。

此外，在本发明的X光双能CT重建方法中，其特征在于，

还具有如下步骤：利用双效应分解，根据a₂＝r₂×diag（ω₂+ε）以及ρ_e＝2a₂得到电子密度图ρ_e。

此外，在本发明的X光双能CT重建方法中，其特征在于，

还具有如下步骤：利用双效应分解，根据来得到等效原子序数图Z^eff。

此外，在本发明的X光双能CT重建方法中，其特征在于，

在所述步骤（d）中，根据所述高能数据p_H重建高能情况下的有效线衰减系数μ_H，取ω₁ =μ_H并且以作为第二局部平滑限制条件的下式（8）重建r₁，

 s.t.                    (8)，

然后，利用a₁＝r₁×diag（ω₁+ε）以及来得到等效原子序数图Z^eff。

此外，在本发明的X光双能CT重建方法中，其特征在于，

在所述步骤（c）和步骤（d）中，利用ART+TV法对r₁和r₂进行重建。

此外，在本发明的X光双能CT重建方法中，其特征在于，

在所述步骤（c）和步骤（d）中，利用split Bregman方法对r₁和r₂进行重建。

根据本申请发明，与现有技术相比，具有如下效果：（1）通过有效地利用数据固有信息（例如，选择ω₁和ω₂为μ_H或μ_L，或者选择ω₁为a₂），从而能够在保持分辨率的情况下有效地抑制双能重建图像的噪声；（2）使用先验模型的方式建立重建方式，从而能够方便地进行算法设计；（3）该方法不局限于扫描方式，也适用于扇束、锥束，圆轨道、螺旋轨道等不同扫描方式，能够通过该先验方法增加迭代重建的鲁棒性；（4）与本领域内已有的双效应分解方法相比，直接重建系数的比值能得到更稳定的结果。

附图说明

图1是本发明的以双效应分解为例的基于自先验信息的双能CT重建方法的流程图。

图2是本发明的以双效应分解为例的基于自先验信息的双能CT重建方法重建得到的图像，（a）为电子密度重建图，（b）为等效原子序数重建图。

具体实施方式

关于自先验，这是发明人提出的一种说法，因为在重建如下所述的比例图像r₁和r₂的过程中所使用的先验模型能够从数据本身得到，例如，之前按照常规CT方法重建单一高能数据所得到的线衰减系数结果，或者重建a₁（分解系数）时所使用的ω₁ = a₂（a₂也是分解系数），即，利用a₂的结构信息作为先验信息。此外，在本申请发明中，处理数据本身获得单一能量的衰减图像重建或者a₂的重建，将它们放在先验模型里，所以，称之为利用数据内在的信息，由此，能够得到获得高质量的重建图像。

以下，参照附图对本发明的实施方式进行说明。

首先，进行能谱标定和双能查找表建立。然后，利用双能CT成像系统的探测器采集得到双能CT成像系统的高能数据和低能数据，此处，设双能CT的高能数据和低能数据分别为p_H和p_L，其中，P_H以及P_L以如下的式（1）和式（2）表示，

                 （1），

                 （2），

其中，上述式（1）和式（2）中的、分别为归一化的高能和低能的能谱分布，这些能谱分布能够利用多种方式来产生，包括例如两层三明治式探测器的伪双能、快速切换的Ｘ光机产生的两个能量谱或者使用两个Ｘ光机所得到的能谱分布等。此外，μ(E)是物体的线衰减系数，H是投影矩阵。按照传统的CT重建方法，能够得到高能以及低能下的有效线衰减系数分布μ_H和μ_L的估计和，并且，和能够作为先验来使用。

目前，在本领域中，双能分解包括基材料分解和双效应分解这两种分解方式，并且，这两种分解都能够表示为如下的式（3），

                            （3），

其中，上述式（3）中的和在基材料分解和双效应分解中表现为不同的预先确定的函数形式，此外，a₁和a₂为分解系数。

此外，在本发明的方法中引入两个新的向量，即，比例图像r₁和r₂，r₁和r₂以如下的式（4）表示，

                              （4），

在上述式（4）中，diag()表示对角线矩阵，并且，其对角线上的元素为括号内的向量的值，ε为一个小值常系数向量，以避免分母出现零值，ω₁和ω₂为能够任意选择的两个向量。此外，定义，，综合上述的式（1）～（4），能够得到如下的式（5），

           （5）。

此外，令

                     （6），其中，R₁和R₂分别是比例图像r₁和r₂的投影图像，和分别为对应的投影算子。

接着，利用式（6）将式（5）简化，得到如下的式（7），

                （7）。

此处，能够将上述式（7）称为除式先验定义下的数据。此外，对于每一对采集得到的高能数据和低能数据，式（7）形成一对非线性二元方程组，通过解此方程组，由此，能够根据(p_H, p_L)得到(R₁, R₂)，或者，也能够使用类似于非专利文献2中所述的方法，使用建立查找表根据已知的数据对(p_H, p_L)得到(R₁, R₂)。剩下的问题就是由(R₁, R₂)重建(r₁, r₂)。根据上述的式（6）可知，(r₁, r₂)的重建能够使用任何的传统CT重建方法完成。但是，在利用传统的双能CT重建方法重建(r₁, r₂)的情况下，由于在双能分解的过程中带有很强的不平衡性，因此，导致双能CT重建图像的噪声被放大，特别是等效原子序数分布的噪声被放大。

相对于此，本发明方法的特征在于通过限定r₁和r₂的局部平滑（piece-wise smooth），即其数学表示可稀疏化，使a₁与ω₁具有相似性、a₂和ω₂具有相似性，能够在重建中利用这个性质改善重建图像的质量。a₁与ω₁都是图像，它们具有相似性指的是图像的结构具有相似性，例如，a₁平滑的地方ω₁也平滑，a₁有边缘的地方ω₁可能也有边缘。此外，由于ω₁和ω₂与μ_H或μ_L之间的结构相似性，所以，能够选择ω₁和ω₂为μ_H或μ_L（其中，μ_H或μ_L能够分别从p_H和p_L按照传统单能CT重建得到）或者其他的与ω₁和ω₂具有相似性的先验图像（例如，在本发明中，在双效应分解方法中能够使用a₂）。此外，在取ω₁和ω₂为均匀常数向量的情况下，此方法退化为常规的双能重建方法。此外，在本发明中，为了改善重建图像的质量而使用以下的条件即式（8）（第二局部平滑限制条件）以及式（9）（第一局部平滑限制条件）的方式来实现对r₁和r₂的局部平滑限制，

 s.t.                     (8)，

 s.t.                     (9)，

在上述的式（8）以及式（9）中，表示r的梯度的p阶范数，此处，通过令其最小化，从而实现r的局部平滑。

此外，在上述的式（8）和式（9）中，除了下标以外完全一致，所以，r₁和r₂的重建能够相互独立地用同样的方法实现。但是，另一方面，从双能CT来说，为了优化重建图像的质量，能够根据实际情况从双能CT的特点出发分别优化选择ω₁和ω₂。例如，对双效应分解方法（参照非专利文献2）来说，通常等效原子序数重建图像的信噪比比较差，而电子密度图像具有较好的信噪比。此外，基于双效应分解方法，等效原子序数Z^eff和电子密度ρ_e的计算公式为如下的式（10），

             （10），

其中λ是双效应分解中光电效应随能量的变化关系的参数，能够选择ω₁=a₂。这样，得到如下的式（11），

            （11）。

这样，能够更好地控制等效原子序数重建图像的噪声水平，这是因为能够通过限制r₁的局部平滑而使等效原子序数的值相对稳定。

此处，在本发明中，具体地以基于双效应分解的双能CT重建为例，使用本发明的方法给出如下的具体实施例。

图1是以双效应分解为例的基于自先验信息的双能CT重建方法的流程图，图2是以双效应分解为例的基于自先验信息的双能CT重建方法重建得到的图像，（a）为电子密度重建图，（b）为等效原子序数重建图。

在本实施例中，如图1所示，首先，进行能谱标定和双能查找表建立。关于建立查找表，使用确定的已知的材料做成不同厚度，从而(R₁, R₂)已知，把它们放在双能CT系统采集数据，获得(p_H, p_L)，把这些数据做成表格，一般取大于两种材料，取几十个不同厚度。更细的描述可以参考本领域内已发表的文献（例如非专利文献2）。

然后，利用双能CT成像系统的探测器采集得到双能CT成像系统的高能数据p_H以及低能数据p_L。

接着，通过查找表或者根据上述的二元方程组（7）得到比例图像r₁和r₂的投影图像R₁和R₂的值。

接着，使用传统的单能CT重建方法由p _H重建μ_H，并且，选择ω₂=μ_H，即，此处将μ_H的结构信息作为先验信息。

接着，根据上述的算式来得到，然后，根据上述的式（9）（第一局部平滑限制条件）即 s.t. ，重建r₂。

此外，在本领域内有多种方法可供选择以实现上述优化问题，一般通过迭代的方式完成。在本申请中以ART+TV方法为例进行说明。

1) 初始化r₂ = 1。

n = 1， ;

2) for :

ART迭代: for :

;  其中，i为射线索引号，[H] _i 为系统矩阵的第i行向量，N_iter为总的迭代次数。

3) 对向量的每个元素进行正则化约束：

，其中，N_ART为ART的跌代次数。

4) 对向量进行全变分最小化迭代:

， 

全变分最速下降法：for ， ，:

， ;

;

5) ; 返回到2）开始下一次迭代，其中，N_TV为TV最小化的跌代次数。

此外，如上所述，在本申请发明中，除了上述的ART+TV方法以外，也可以利用例如split Bregman方法等其他方法进行r₂的重建。

接着，根据上述的式（4）以及式（10）得到分解系数a₂以及电子密度ρ_e，即，计算，由此，能够获得电子密度图，图2（a）是如上得到的电子密度重建图，如图2（a）所示那样，电子密度图的边缘清晰。

接着，取ω₁＝a₂，从而根据式得到，即，此处将a₂的结构信息作为先验信息。

接着，根据上述的式（8）（第二局部平滑限制条件）即 s.t. ，重建r₁。并且，重建r₁的具体实现与上述的重建r₂的步骤相同。

接着，根据式（11）得到等效原子序数，由此，能够获得等效原子序数图像。图2（b）是如上得到的等效原子序数重建图，如图2（b）所示那样，等效原子序数图像在同一物质的局部区域的噪声均匀性得到很大提高，无异常点出现，能够很好地重建细条状物质。

此外，关于r₁的重建，也可以与上述不同，取ω₁＝μ_H（即，将μ_H的结构信息作为先验信息）并且以作为第二局部平滑限制条件的式（8）重建r₁，然后，利用a₁＝r₁×diag（ω₁+ε）以及来得到等效原子序数图Z^eff。这样，也与上述同样地，能够得到个高质量的等效原子序数图像。

如上述那样，以基于双效应分解的双能CT重建为例进行了说明，但是并不限于此，本发明的方法也能够应用于基材料分解。

如上所述，在本发明中，通过有效地利用数据固有信息，从而能够在保持分辨率的情况下有效地抑制双能重建图像的噪声，并且，使用了先验模型的方式建立重建方式，从而能够方便地进行算法设计，此外，本发明方法不局限于扫描方式，也适用于扇束、锥束，圆轨道、螺旋轨道等不同扫描方式，能够通过该先验方法增加迭代重建的鲁棒性，此外，本发明的X光双能CT重建方法与本领域内已有的双效应分解方法相比，直接重建系数的比值能得到更稳定的结果。

Claims (11)

1. 一种X光双能CT重建方法，其特征在于，具有如下步骤：

（a）利用双能CT成像系统的探测器采集得到双能CT成像系统的高能数据p_H以及低能数据p_L；

（b）根据所得到的所述高能数据p_H以及所述低能数据p_L得到比例图像r₁以及r₂的投影图像R₁以及R₂；

（c）利用第一局部平滑限制条件重建比例图像r₂，并且，得到分解系数a₁；

（d）利用第二局部平滑限制条件重建比例图像r₁，并且，得到分解系数a₂。

2. 如权利要求1所述的X光双能CT重建方法，其特征在于，

比例图像r₁以及r₂以如下的式（4）定义，

           （4），

投影图像R₁以及R₂以如下的式（6）定义，

                   （6），

其中，，，H是投影矩阵，a₁和a₂为分解系数，ε为一个小值常系数向量，ω₁和ω₂为能够任意选择的向量，

在所述步骤（c）中，以作为第一局部平滑限制条件的下式（9）重建r₂，

 s.t.                  (9)，

在所述步骤（d）中，以作为第二局部平滑限制条件的下式（8）重建r₁，

 s.t.                     (8)。

3. 如权利要求2所述的X光双能CT重建方法，其特征在于，

在所述步骤（c）中，根据所述高能数据p_H重建高能情况下的有效线衰减系数μ_H，并且，选择ω₂ ＝ μ_H，以作为第一局部平滑限制条件的式（9）重建r₂。

4. 如权利要求2或3所述的X光双能CT重建方法，其特征在于，

在上述步骤（d）中，取ω₁ ＝ a₂并且以作为第二局部平滑限制条件的式（8）重建r₁。

5. 如权利要求2或3所述的X光双能CT重建方法，其特征在于，

还具有如下步骤：利用双效应分解，根据a₂＝r₂×diag（ω₂+ε）以及ρ_e＝2a₂得到电子密度图ρ_e。

6. 如权利要求2或3所述的X光双能CT重建方法，其特征在于，

还具有如下步骤：利用双效应分解，根据来得到等效原子序数图Z^eff。

7. 如权利要求4所述的X光双能CT重建方法，其特征在于，

还具有如下步骤：利用双效应分解，根据a₂＝r₂×diag（ω₂+ε）以及ρ_e＝2a₂得到电子密度图ρ_e。

8. 如权利要求4所述的X光双能CT重建方法，其特征在于，

还具有如下步骤：利用双效应分解，根据来得到等效原子序数图Z^eff。

9. 如权利要求2所述的X光双能CT重建方法，其特征在于，

在所述步骤（d）中，根据所述高能数据p_H重建高能情况下的有效线衰减系数μ_H，取ω₁ ＝ μ_H并且以作为第二局部平滑限制条件的下式（8）重建r₁，

 s.t.                     (8)，

然后，利用a₁＝r₁×diag（ω₁+ε）以及来得到等效原子序数图Z^eff。

10. 如权利要求2所述的X光双能CT重建方法，其特征在于，

在所述步骤（c）和步骤（d）中，利用ART+TV法对r₁和r₂进行重建。

11. 如权利要求2所述的X光双能CT重建方法，其特征在于，

在所述步骤（c）和步骤（d）中，利用split Bregman方法对r₁和r₂进行重建。