WO2015010632A1

WO2015010632A1 - 一种x光双能ct重建方法

Info

Publication number: WO2015010632A1
Application number: PCT/CN2014/082900
Authority: WO
Inventors: 刑宇翔; 沈乐; 张丽; 陈志强; 李亮
Original assignee: 清华大学; 同方威视技术股份有限公司
Priority date: 2013-07-26
Filing date: 2014-07-24
Publication date: 2015-01-29
Also published as: HK1204124A1; CN104346820B; CN104346820A; EP2843623A2; US20150030225A1; US9513235B2; EP2843623B1; PL2843623T3; EP2843623A3

Abstract

一种基于自先验信息的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) 进行能谱标定和双能查找表建立；

(b)利用双能 CT成像系统的探测器采集得到双能 CT成像系统的高能数据 p_H以及低能数据

(c) 根据所得到的所述高能数据 p_H以及所述低能数据寻到比例图像以及 r₂的投影图像以及 R_{2 ;}

(d)利用第一局部平滑限制条件重建比例图像 r₂，由此，得到电子密度图； (e) 利用第二局部平滑限制条件重建比例图像 1^，由此，得到等效原子序数图像。

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

比例图像！^以及 r₂以如下的式（4) 定义，

' 1

ι ≡ diag

)

(4)，

、

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

R_1≡H;_ri

(6)，

R₂≡ H;r₂ 其中， H;

+ s)， H;≡Hdiag(co₂ + ε)， H是投影矩阵， a Pa₂为分解系数， ε为一个小值常系数向量， co^Pco₂为能够任意选择的向量，在所述步骤（d) 中，根据所述高能数据 p_H重建高能情况下的有效线衰减系数 μ_Η，选择 ω₂=μ_Η，并且，以作为第一局部平滑限制条件的下式（9) 重建 r₂， mm Vi I s.t. R₂≡H;r₂ (9)，然后，利用 a₂ = r₂xdiag (ωι+ε) 以及 p_e=2a₂得到电子密度图 p_e，在上述步骤（e)中，取 co^a并且以作为第二局部平滑限制条件的下式（8) 重建 r ,

(8)，

然后，利用 Z^rff «_ri^,A«4来得到等效原子序数图 Z^rff 。

此外，在本公开的 X光双能 CT重建方法中，其特征在于，在上述步骤（e)中，取 _ωι = μ_Η并且以作为第二局部平滑限制条件的下式（8)

mm II S.t. R,≡ Η'^ (8)， i、 λ - 然后，利用 a^r diag ( ωι+ε) 以及 Z^eff = diag , A « 4来得到等效原子序数图 Z^rff 。

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

在上述步骤（d) 和步骤（e) 中，利用 ART+TV法对！^和！^进行重建。此外，在本公开的 X光双能 CT重建方法中，其特征在于，

在上述步骤（d) 和步骤（e) 中，利用 split Bregman方法对！^和 r₂进行重建。

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

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

图 2是本公开的以双效应分解为例的基于自先验信息的双能 CT重建方法重建得到的图像，（a) 为电子密度重建图，（b) 为等效原子序数重建图。具体实施方式

关于自先验，这是发明人提出的一种说法，因为在重建如下所述的比例图像！^和！^的过程中所使用的先验模型能够从数据本身得到，例如，之前按照常规

CT方法重建单一高能数据所得到的线衰减系数结果，或者重建₃₁ (分解系数）时所使用的 co a ( a₂也是分解系数）， BP , 利用 a₂的结构信息作为先验信息。此外，在本申请公开中，处理数据本身获得单一能量的衰减图像重建或者 a₂的重建，将它们放在先验模型里，所以，称之为利用数据内在的信息，由此，能够得到获得高质量的重建图像。

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

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

- In εχρ(-Ημ( ) Κ = p_H (1)，

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

^，并且， ^和^能够作为先验来使用。目前，在本领域中，双能分解包括基材料分解和双效应分解这两种分解方式，并且，这两种分解都能够表示为如下的式（3)， μ(Ε) = _{Άι ί}(Ε) + ^ ₂(Ε) (3)，其中，上述式（3) 中的和 (E)在基材料分解和双效应分解中表现为不同的预先确定的函数形式，此外， a^Pa₂为分解系数。此外，在本公开的方法中引入两个新的向量，即，比例图像！^和！^，！^和 r₂以如下的式（4) 表示，

' 1 )

ι ≡ diag ― a_:

+ - εε )

(4)， 1

r₂≡ diag

ω, + ε

在上述式（4) 中， diag()表示对角线矩阵，并且，其对角线上的元素为括号内的向量的值， ε为一个小值常系数向量，以避免分母出现零值， co^Pco₂为能够任意选择的两个向量。此外，定义11;_≡11(113_§(( ₁ + ， H'₂≡Hdiag(_W2 + _£) , 综合上述的式 (1) 〜（4)，能够得到如下的式（5)，

-In fw_H (E) exp (― Η;ι·_ι(¾(£·) - H'₂r₂^₂ (E))dE = p_H

^J. (5) 。

- In j w_L (E) exp (-Η _ΓΙ^ (E) - H'₂r₂^₂ (E))dE = p_L 此外，令

^{Rl ≡} (6), R₂≡ H;r₂ 其中， 1^和分别是比例图像！^和！^的投影图像，和 H'₂分别为对应的投影算子。

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

(7)。

— lnjv_L( )_eXp(— R^^)— R ₂( ) ffi = p_L

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

相对于此，本公开方法的特征在于通过限定 r^Pr₂的局部平滑（piece-wise smooth), 即其数学表示可稀疏化，使 ^与^具有相似性、 a Pco具有相似性，能够在重建中利用这个性质改善重建图像的质量。 ^与^都是图像，它们具有相似性指的是图像的结构具有相似性，例如，平滑的地方 cq ^p , ^有边缘的地方 0¾可^也有边缘。此外，由于 to Pto₂与 μ_Η或之间的结构相似性，所以，能够选择₍0₁和₍0₂为 ₁₁或 ₁^ (其中， ₁₁或 ₁ ^能够分别从

p^安照传统单能 CT 重建得到）或者其他的与 ^和^具有相似性的先验图像（例如，在本公开中，在双效应分解方法中能够使用 a₂ )。此外，在取 co^P co₂为均匀常数向量的情况下，此方法退化为常规的双能重建方法。此外，在本公开中，为了改善重建图像的质量而使用以下的条件即式（8) (第二局部平滑限制条件）以及式（9) (第一局部平滑限制条件）的方式来实现对！^和！"的局部平滑限制， mm II L s.t. R_1≡H _ri (8), min|| S.t. R₂≡ H'₂r₂ (9)，在上述的式（8) 以及式（9) 中， ||V 表示 r的梯度的 p阶范数，此处，通过令其最小化，从而实现 r的局部平滑。

此外，在上述的式（8) 和式（9) 中，除了下标以外完全一致，所以， r P

1"₂的重建能够相互独立地用同样的方法实现。但是，另一方面，从双能 CT来说，为了优化重建图像的质量，能够根据实际情况从双能 CT的特点出发分别优化选择 co^P co₂。例如，对双效应分解方法（参照非专利文献 2) 来说，通常等效原子序数重建图像的信噪比比较差，而电子密度图像具有较好的信噪比。此外，基于双效应分解 ^eff和电子密度 p_e的计算公式为如下的式（ 10)，

( 10)，

p = 2a 其中 λ是双效应分解中光电效应随能量的变化关系的参数，能够选择 co^ a^ 这样，得到如下的式（11 )，

这样，能够更好地控制等效原子序数重建图像的噪声水平，这是因为能够通过限制！^的局部平滑而使等效原子序数的值相对稳定。此处，在本公开中，具体地以基于双效应分解的双能 CT重建为例，使用本公开的方法给出如下的具体实施例。图 1 是以双效应分解为例的基于自先验信息的双能 CT 重建方法的流程图，图 2是以双效应分解为例的基于自先验信息的双能 CT重建方法重建得到的图像，（a) 为电子密度重建图，（b) 为等效原子序数重建图。

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

然后，利用双能 CT成像系统的探测器采集得到双能 CT成像系统的高能数据p_H以及低能数据接着，通过查找表或者根据上述的二元方程组（7) 得到比例图像！^和 r₂ 的投影图像 R P R₂的值。接着，使用传统的单能 CT重建方法由重建 μ_Η，并且，选择 ω₂ = μ_Η，即，此处将 μ_Η的结构信息作为先验信息。接着，根据上述的算式 H'_2≡Hd_iag((0_{2 + S})来得到 Η'₂，然后，根据上述的式 ( 9) (第一局部平滑限制条件）即 s.t. R₂ = H'₂r₂，重建 r₂。

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

1) 初始化 r₂ = l。

2) for η=2, . .., N_lter：

ART迭代： for m=2, ■■；：其中，！'为射线索引号， [_H]i为系

3) 对向量 '^Nart的每个元素进行正则化约束： r-^N- = max(r-^N-,0)，其中， NAR_T为 ART的跌代次数。

4) 对向量 r₂"'^N»_T进行全变分最小化迭代： d" = |i ¹ _ i ^Nart ， r₂ ^{n l} = r₂ ^{n Nart} 全变分最速下降法： for fe=2, ...,N_TV， α=0.2， ε= 10"⁸:

5) r₂ ' =r₂ ^TV _; 返回到 2)开始下一次迭代，其中， Ν_ΤΝ为 TV最小化的跌代次数。

此外，如上所述，在本申请公开中，除了上述的 ART+TV方法以外，也可以利用例如 split Bregman方法等其他方法进行 r₂的重建。接着，根据上述的式（4) 以及式（10) 得到分解系数以及电子密度 p_e，即，计算^^^ ^^+^), _Pe=2a₂ , 由此，能够获得电子密度图，图 2 (a)是如上得到的电子密度重建图，如图 2 (a) 所示那样，电子密度图的边缘清晰。

接着，取 _ωι=ίΐ₂，从而根据式！！^！^ ^+ 得到!^， BP,此处将的结构信息作为先验信息。

接着，根据上述的式（8) (第二局部平滑限制条件）即 ηπηΙνηΙΙ s.t. R^H^ , 重建 η。并且，重建的具体实现与上述的重建 r₂的步骤相同。接着，根据式（11)得到等效原子序数 Ζ »!·^¹, 由此，能够获得等效原子序数图像。图 2 (b)是如上得到的等效原子序数重建图，如图 2 (b)所示那样，等效原子序数图像在同一物质的局部区域的噪声均匀性得到很大提高，无异常点出现，能够很好地重建细条状物质。

此外，关于！^的重建，也可以与上述不同，取 _ωι = μ_Η (BP, 将 ₁₁的结构信息作为先验信息）并且以作为第二局部平滑限制条件的式（8) 重建，然后，利 i、 λ - 用 a₁ = r₁xdiag ( ωι+ε) 以及 Z^eff = diag , A « 4来得到等效原子序数图

Z^rff 。这样，也与上述同样地，能够得到个高质量的等效原子序数图像。

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

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

Claims

权利要求书

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

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

(b) 根据所得到的所述高能数据 p_H以及所述低能数据寻到比例图像以及 r₂的投影图像以及 R_{2 ;}

(c)利用第一局部平滑限制条件重建比例图像 r₂，并且，得到分解系数 a_1;

(d)利用第二局部平滑限制条件重建比例图像 1^，并且，得到分解系数。

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

比例图像！^以及 r₂以如下的式（4) 定义，

' 1 )

ι ≡ diag ―

+ - εε I

(4)，

1

r₂≡ diag

ω + ε 投影图像1^以及 R₂以如下的式（6) 定义，

(6)，其中， H;

+ s)， H;≡ Hdiag(w₂ + ε)， Η是投影矩阵， ₁和 ₂为分解系数， ε为一个小值常系数向量， co^Pco₂为能够任意选择的向量，

在所述步骤（c) 中，以作为第一局部平滑限制条件的下式（9) 重建 r₂， mm V] I s.t. R₂≡H;r₂ (9)，

在所述步骤（d) 中，以作为第二局部平滑限制条件的下式（8) 重建， minllVrJI s.t. R,≡ H _: (8)

η 。

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

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

4. 如权利要求 2或 3所述的 X光双能 CT重建方法，其中，在上述步骤（d) 中，取 co^ a并且以作为第二局部平滑限制条件的式（8)

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

还具有如下步骤：利用双效应分解，根据 a₂ = r₂xdiag (ω₂+ε) 以及 p_e=2a₂ 得到电子密度图 p_e。

6. 如权利要求 2或 3所述的 X光双能 CT重建方法，其中，还具有如下步骤：利用双效应分解，根据 Z^rff « _ri^, A « 4来得到等效原子序数图 Z^rff 。

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

8. 如权利要求 4所述的 X光双能 CT重建方法，其中，还具有如下步骤：利用双效应分解，根据 Ζ » !·^¹, A « 4来得到等效原子序数图 Z^rff 。

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

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

mm L s.t. R_1≡H _ri (8),

1

然后，利用 a₁ = r₁xdiag (ωι+ε) 以及 Z^eff = diag , A « 4来得到等效原子序数图 Z^rff 。

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

在所述步骤（c) 和步骤（d) 中，利用 ART+TV法对！^和！^进行重建。

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

在所述步骤（c) 和步骤（d) 中，利用 spht Bregman方法对！^和！^进行重建