CN103136772B - 基于加权阿尔法散度约束的x射线低剂量ct图像重建方法 - Google Patents

Abstract

本发明公开了一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法，包括：（1）利用CT成像设备获取低剂量CT投影数据以及成像系统参数；（2）采用阿尔法散度测度作为原始含噪声的投影数据与待恢复的投影数据的距离测度，并根据获取的系统参数计算阿尔法散度测度的权重因子，构建基于加权阿尔法散度约束的投影数据恢复模型；（3）对构建的投影数据恢复模型进行目标函数求解，建立迭代算法格式；（4）对获取的低剂量CT投影数据，应用建立的迭代算法格式对投影数据恢复模型进行迭代求解；（5）对获得的恢复后投影数据进行图像重建。本发明公开的低剂量CT图像重建方法在CT图像的噪声抑制、边缘保持方面均有上佳表现。

Description

基于加权阿尔法散度约束的X射线低剂量CT图像重建方法

技术领域

本发明涉及一种医学影像的断层图像投影数据恢复方法，具体来说涉及一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法。

背景技术

虽然X射线CT已经在医学影像诊断中广泛应用，但扫描中过高的X射线剂量使用将对人体造成不可预测的伤害。因此，在保证图像质量的前提下，最大限度地降低X射线使用剂量已成为医学CT成像领域的迫切需要。

当前，降低扫描中的管电流（mA）或管电压（kVp）是实现低剂量CT成像的最便捷且最常用的方法。但是在降低的管电流（mA）或管电压（kVp）条件下采集的投影数据中含有大量的噪声，使得基于传统滤波反投影方法重建的图像会出现严重的退化现象，难以满足临床诊断需求。因此，诸多基于降低管电流或管电压的扫描协议的CT图像重建方法相继提出，以期在保证图像质量前提下大幅降低X射线辐射剂量，如基于统计模型的图像迭代重建方法，基于投影数据恢复的解析重建方法。大量研究表明，为了实现基于降低的管电流（mA）或管电压（kVp）条件下的CT图像重建，投影数据统计特性的引入是其实现的关键技术之一。由于X射线低剂量CT投影数据的噪声统计特性非常复杂，传统的基于投影数据高斯统计分布的最小二乘距离测度不能准确地描述原始的CT投影数据与恢复后的投影数据之间的距离。由于阿尔法散度测度可以准确测量两种不同分布之间的统计分布距离，因此，本专利通过引入阿尔法散度测度来刻画原始的CT投影数据与恢复后的投影数据之间的统计关系，同时引入CT探测器各探测通道上投影数据的噪声方差作为阿尔法散度测度的权重因子，提出了一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法，可以大幅提升在降低的管电流（mA）或管电压（kVp）条件下的CT重建图像质量。本专利公开方法，相比已有的基于投影数据恢复的X射线低剂量CT图像重建方法，在CT图像的噪声抑制、边缘保持方面均有上佳表现。

发明内容

本发明的目的在于提供一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法，该方法可以大幅提升在降低的管电流（mA）或管电压（kVp）条件下的CT重建图像质量。

本发明的目的可通过以下的技术措施来实现：

一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法，包括以下步骤：

（1）利用CT成像设备采用降低管电流（mA）或管电压（kVp）的扫描协议获取低剂量CT投影数据以及相应的成像系统参数；

（2）采用阿尔法散度测度作为原始含噪声的投影数据与待恢复的投影数据的距离测度，并利用步骤（1）获取的系统参数计算阿尔法散度测度的权重因子，构建基于加权阿尔法散度约束的投影数据恢复模型；

（3）对步骤（2）中构建的投影数据恢复模型进行目标函数求解，建立迭代算法格式，并设置迭代终止条件；

（4）对步骤（1）中获取的低剂量CT投影数据，应用步骤（3）中建立的迭代算法格式对步骤（2）中构建的投影数据恢复模型进行迭代求解；

（5）对步骤（4）获得的恢复后投影数据采用CT图像重建方法进行图像重建。

所述步骤（1）中获取的成像系统参数为降低管电流（mA）或管电压（kVp）的扫描协议下CT探测器各探测通道对应的投影数据噪声方差其中i表示探测器探测通道的位置，I表示所有探测器探测通道的个数。

所述步骤（2）中构建的基于加权阿尔法散度约束的投影数据恢复模型为：

$q^{*}=\arg \underset{q\≥0}{\min }{D}_{\mathrm{\α}}^w(y,q)+\mathrm{\λR}\left(q\right),$

其中 $D_{\mathrm{\α}}^w(y,q)=\frac{1}{\mathrm{\α}(1-\mathrm{\α})}{\mathrm{\Σ}}_{i=1}^I{w}_i[\mathrm{\α}{y}_i+(1-\mathrm{\α})q_i-y_i^{\mathrm{\α}}{q_i}^{1-\mathrm{\α}}]$ 表示α加权散度测度，α为实数参数，即α∈(-∞,+∞)；y={y_i,i＝1,2,…,I}表示采集到的CT投影数据；q={q_i,i＝1,2,…,I}表示待恢复的投影数据；为α散度测度的权重因子，为获取的投影数据噪声方差；R(q)为先验约束项；λ>0为正则化参数，用于刻画先验约束的强度。

所述步骤（3）中的迭代算法格式为高斯-塞德尔（Gauss-Seidel，GS）迭代形式，即： $q_i^{n+1}=\frac{y_i^{\mathrm{\α}}+\mathrm{\λ}{\mathrm{\σ}}_i^2{\left(q_i^n\right)}^{\mathrm{\α}}\underset{m\∈N_i}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{im}}{q}_m^n}{{\left(q_i^n\right)}^{\mathrm{\α}-1}\{1+\mathrm{\λ}{\mathrm{\σ}}_i^2{q}_i^n\underset{m\∈N_i}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{im}}\}},$ 其中是表示第n步迭代过程中的迭代前投影数据，是第n步迭代恢复后的投影数据。

所述步骤（4）中的迭代求解的终止条件为：相邻两次迭代恢复的投影数据之间的均方根误差（RootofMeanSquareError,RMSE）小于0.001。

所述步骤（5）中的图像重建方法可以为：滤波反投影方法或卷积反投影方法。

本发明方法相比现有方法具有以下有益效果：

1、本发明方法通过引入阿尔法散度测度来刻画原始的CT投影数据与恢复后的投影数据之间的统计关系，同时引入CT探测器各探测通道上投影数据的噪声方差作为阿尔法散度测度的权重因子，对低剂量CT投影数据进行恢复，并实现在降低的管电流（mA）或管电压（kVp）条件下的CT图像优质重建；

2、本发明方法较已有方法能够较好地保持图像分辨率和抑制图像噪声。

附图说明

图1是本发明基于加权阿尔法散度约束的X射线低剂量CT图像重建方法（WAD-QM）的流程图；

图2(a)是修定的Shepp-Logan体模图像；

图2(b)是基于Ramp滤波的滤波反投影重建图像；

图2(c)是基于Hanning窗滤波的滤波反投影重建图像；

图2(d)是基于已有的加权最小二乘惩罚约束方法（PWLS-QM）恢复的重建图像；

图2(e)是本发明公开方法的重建图像；

图3(a)是用于性能评估的数字体模图像；

图3(b)是WAD-QM和PWLS-QM重建图像的噪声-分辨率曲线。

具体实施方式

本发明公开的基于加权阿尔法散度约束的X射线低剂量CT图像重建方法的具体实施步骤如图1所示，具体如下：

1、利用CT成像设备采用降低管电流（mA）或管电压（kVp）的扫描协议获取低剂量CT投影数据以及相应的成像系统参数，射线剂量为标准剂量的1/10至1/20。上述系统参数为降低管电流（mA）或管电压（kVp）的扫描协议下CT探测器各探测通道对应的投影数据噪声方差其中i表示探测器探测通道的位置，I表示所有探测器探测通道的个数；

2、根据低剂量CT投影数据的噪声统计特性以及降低管电流（mA）或管电压（kVp）的扫描协议下CT探测器各探测通道对应的投影数据噪声方差的非均一性，引入信息论中的阿尔法散度测度作为原始含噪声的投影数据与待恢复的投影数据的距离测度，并利用步骤1获取的系统参数计算阿尔法散度测度的权重因子，构建基于加权阿尔法散度约束的投影数据恢复模型：

$q^{*}=\arg \underset{q\≥0}{\min }{D}_{\mathrm{\α}}^w(y,q)+\mathrm{\λR}\left(q\right),$

其中 $D_{\mathrm{\α}}^w(y,q)=\frac{1}{\mathrm{\α}(1-\mathrm{\α})}{\mathrm{\Σ}}_{i=1}^I{w}_i[\mathrm{\α}{y}_i+(1-\mathrm{\α})q_i-y_i^{\mathrm{\α}}{q_i}^{1-\mathrm{\α}}]$ 表示α加权散度测度，α为实数参数，即α∈(-∞,+∞)，可以取值为1.3；y={y_i,i＝1,2,…,I}表示采集到的CT投影数据；q={q_i,i＝1,2,…,I}表示待恢复的投影数据；为α散度测度的权重因子，为获取的投影数据噪声方差；R(q)为先验约束项，可采用传统的二次平板（QuatraticMembrane，QM）先验形式，即ω_m为局部小方形邻域N_i内的权重值，其中小方形邻域N_i可以取为3×3的八邻域；λ>0为正则化参数，用于刻画先验约束的强度，λ可以取值为4×10^-5。

3、对步骤2中构建的投影数据恢复模型采用高斯-塞德尔方法进行目标函数求解，建立高斯-塞德尔迭代算法格式，并设置迭代终止条件，即相邻两次迭代恢复数据之间的均方根误差（RootofMeanSquareError,RMSE）小于0.001；

4、将步骤1采集到的低剂量CT投影数据作为初始的迭代前投影数据，应用步骤3构建的迭代模型： $q_i^{n+1}=\frac{y_i^{\mathrm{\α}}+\mathrm{\λ}{\mathrm{\σ}}_i^2{\left(q_i^n\right)}^{\mathrm{\α}}\underset{m\∈N_i}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{im}}{q}_m^n}{{\left(q_i^n\right)}^{\mathrm{\α}-1}\{1+\mathrm{\λ}{\mathrm{\σ}}_i^2{q}_i^n\underset{m\∈N_i}{\mathrm{\Σ}}{\mathrm{\ω}}_{\mathrm{im}}\}},$ 进行迭代求解。其中是表示第n步迭代过程中的迭代前投影数据，是第n步迭代恢复后的投影数据。

5、判断是否满足迭代终止条件，若不满足，则将步骤4中的迭代前的投影数据更新为步骤4中迭代后的投影数据，重复步骤4-5，直至满足迭代终止条件；

6、对步骤5获得的恢复后投影数据采用滤波反投影方法或卷积反投影方法进行最终的CT图像重建。

下面对具体数据采用本发明方法重建图像来说明本发明方法的效果。

首先采用如图2(a)所示的修定的Shepp-Logan体模图像作为本发明的模拟实验对象。体模大小设为512×512，射源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器，探测器单元为1.407mm。通过转换概率矩阵K得到投影数据（sinogram），然后通过调整光子总数值模拟生成低剂量CT投影数据。对比实验中图像重建均使用传统扇形束FBP算法，其中汉宁（Hanning）窗滤波的截止频率设为奈奎斯特频率的80%。

图2(b)至图2(e)描述了本发明公开方法和其他方法的效果对比。图2(b)为低剂量数据采用Ramp滤波后的FBP重建图像，图2(c)低剂量数据采用hanning窗滤波后的FBP重建图像，图2(d)为已有的加权最小二乘惩罚约束方法（PenalizedWeightedLeast-Squares-QuatraticMembrane,PWLS-QM）的重建图像，图2(e)为本发明公开方法的重建图像。通过对比可以看出本发明公开方法保持了良好的边缘并且在高衰减区域有更佳的抑制噪声的效果。

表1不同方法的重建图像的信噪比

方法	FBP-Ramp	FBP-Hanning	PWLS-QM	WAD-QM
					信噪比（dB）	23.0939	25.1946	27.7024	30.4681

表1列出了图2(b)-(e)中所示不同方法的FBP重建图像的信噪比。对比可知，本发明公开方法在图像信噪比方面有着上佳的表现。

为了进一步比较分析本发明公开的新方法，对WAD-QM和PWLS-QM重建方法在抑制噪声及保持分辨率方面的性能进行定量分析．实验选择的性能评估数字体模数据如图3(a)所示，主要目的是研究图中实线所示的轮廓线位置（图像中垂直经过白色嵌入物中心的轮廓线）附近的噪声抑制及分辨率保持情况。

图3(b)描述了WAD-QM和PWLS-QM重建图像的噪声-分辨率曲线，可以看出本发明所提的方法相比已有的PWLS-QM方法在分辨率保持和噪声抑制方面均有一定的改善。

本发明的实施方式不限于此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (4)

1.一种基于加权阿尔法散度约束的X射线低剂量CT图像重建方法，其特征包括以下步骤：

(1)利用CT成像设备采用降低管电流或管电压的扫描协议获取低剂量CT投影数据以及相应的成像系统参数；

(2)采用阿尔法散度测度作为原始含噪声的投影数据与待恢复的投影数据的距离测度，并利用步骤(1)获取的成像系统参数计算阿尔法散度测度的权重因子，构建基于加权阿尔法散度约束的投影数据恢复模型；

(3)对步骤(2)中构建的投影数据恢复模型进行目标函数求解，建立迭代算法格式；并设置迭代终止条件；

(4)对步骤(1)中获取的低剂量CT投影数据，应用步骤(3)中建立的迭代算法格式对步骤(2)中构建的投影数据恢复模型进行迭代求解；

(5)对步骤(4)获得的恢复后投影数据采用CT图像重建方法进行图像重建；

步骤(2)中构建的基于加权阿尔法散度约束的投影数据恢复模型为：

$q^{*}=\arg \underset{q\≥0}{min}{D}_{\mathrm{\α}}^w(y,q)+\mathrm{\λ}R\left(q\right),$

其中 $D_{\mathrm{\α}}^w\left(y,q\right)=\frac{1}{\mathrm{\α}\left(1-\mathrm{\α}\right)}{\Σ}_{i=1}^I{w}_i\[{\mathrm{\αy}}_i+\left(1-\mathrm{\α}\right)q_i-y_i^{\mathrm{\α}}{q_i}^{1-\mathrm{\α}}\]$ 表示α加权散度测度，α为实数参数，即α∈(-∞,+∞)；y＝{y_i,i＝1,2,…,I}表示采集到的CT投影数据；q＝{q_i,i＝1,2,…,I}表示待恢复的投影数据；为α散度测度的权重因子，为获取的投影数据噪声方差；R(q)为先验约束项；λ＞0为正则化参数，用于刻画先验约束的强度；i表示探测器探测通道的位置，I表示所有探测器探测通道的个数。

2.根据权利要求1所述的X射线低剂量CT图像重建方法，其特征在于：所述步骤(1)中获取的成像系统参数为降低管电流(mA)或管电压(kVp)的扫描协议下CT探测器各探测通道对应的投影数据噪声方差其中i表示探测器探测通道的位置，I表示所有探测器探测通道的个数。

3.根据权利要求1所述的X射线低剂量CT图像重建方法，其特征在于：所述步骤(3)中的迭代终止条件为：相邻两次迭代恢复数据之间的均方根误差小于0.001。

4.根据权利要求1所述的X射线低剂量CT图像重建方法，其特征在于：所述步骤(5)中的重建方法为滤波反投影方法或卷积反投影方法。