基于动态优化的X射线CT低剂量重建

王 康，赵 琦，李 铭

(1.南京理工大学 电子工程与光电技术学院，江苏 南京 210094；2.中国科学院 苏州生物医学工程技术研究所，江苏 苏州 215163)

1 引 言

计算机断层成像(Computed Tomography,CT)是通过无损方式获取人体内部解剖信息的一种医学影像技术。CT扫描的普及，给临床诊断带来了极大方便，但也引发了受检者的辐射损伤。由于X射线辐射剂量的累加效应，接受过量的X光照射可能会显著增加病人诱发癌症和白血病等疾病的概率。因此，国际辐射防护委员会建议，在满足临床诊断需求的情况下，CT辐射剂量越低越好[1-2]。

低剂量CT成像技术是减少辐射损伤的有效方法。临床应用中的CT低剂量扫描成像可大致分为两类：降低管电流或管电压和降低采样数[3-4]。前者多适用于投影数据采集完备的螺旋CT设备，后者则更适用于投影数据欠采样的CBCT(Cone-beam CT,CBCT)设备。目前，国内主流的CT设备仍采用滤波反投影法(Filtered Back-projection,FBP)进行重建，采用CT低剂量扫描协议将导致重建图像中伪影和噪声的增加以及对比度的下降，直接影响医生对病变组织的诊断分析。因此，如何保证低剂量扫描下的图像质量是CT成像领域的热点问题之一，而研究迭代重建算法[5-14]则是抑制噪声和伪影，提高CT重建图像质量的关键。Idris等人[5]首先采用复合泊松分布模型对投影数据的噪声特性进行建模，并提出基于这一模型的统计迭代重建算法(Statistical Iterative Reconstruction，SIR)。Lasio等[6]分析了基于泊松分布模型的统计迭代重建算法在实际CT系统中的性能，并验证了泊松分布模型对于实际CT设备采集投影数据能取得较好的重建效果。在此基础上，为了进一步提升重建效果，需改进目标函数中的先验约束项。Zhang等[7]和Wu等[8]将L1范数字典学习(DL)正则化项引入统计迭代重建模型(SIR)，减轻了L2最小化带来的过平滑效果并保留更多图像细节。Gou等[9]测试了一种基于Lp字典学习(DL)稀疏约束的新颖重建算法，结论是较低的p值会呈现出更好的重建性能。基于字典学习的CT低剂量重建对于图像细节有所改善，但是其重建时间较长，且由于低剂量投影数据会导致迭代前期存在伪影，这些伪影会被字典所表示从而影响最终的重建结果。Xu等[10]则结合全变分最小化(Total Variation，TV)，提出了SIR-TV算法。该算法能较好地抑制噪声和伪影，但对于稀疏性表征差的重建目标，重建图像中易出现块状伪影。此后，Sidky等[11]又提出分数阶全变分算法(TpV)。该算法通过调整p值可有效改进待重建目标的稀疏性表征，但是TpV算法易导致重建图像中出现散点噪声。Lei等[12]对TV结构进行改进，通过引入自适应的加权权重，提出了自适应变权全变分算法(Adaptively Reweighted TV,ARWTV)。ARWTV算法对于噪声和伪影有很好的抑制效果，但重建结果中弱对比度组织边缘模糊。近年来，随着深度学习的发展，许多学者将基于深度学习的方法应用于医学领域[15-23]。与传统方法相比，深度学习的图像质量有所提高，但这需要大量的训练数据，且重建质量取决于数据集的质量。目前大多数深度学习中使用的数据是通过噪声插入来模拟的，这可能无法反映低剂量CT扫描的实际噪声分布。相反，传统的基于压缩感知(Compressed Sensing ,CS)的重建方法不需要大量的训练数据，并且可以很容易地集成到传统的迭代重建框架中。

针对CT低剂量扫描下出现的噪声和伪影问题，文中提出一种基于分数阶全变分动态优化(TpV dynamic optimization,TpV-DO)的CT低剂量成像算法。首先通过一族双曲正切函数集构造TpV的动态复合函数模型，以增强待重建目标的稀疏性表征；其次，通过引入动态优化设计以实现对噪声和伪影的抑制以及对弱对比度组织边缘的保护；最后，结合统计迭代重建框架，使用交替优化算法实现对目标函数的优化求解。采用数值模型和动物扫描获取的投影数据对新算法进行了实验，并与FBP算法、TpV算法、ARWTV算法进行了比较。

2 算法及理论分析

2.1 统计迭代成像模型

基于泊松分布模型来模拟投影数据采集过程，将统计迭代模型与正则项结合的目标函数可以表示为[10]：

(1)

2.2 基于TpV的动态优化模型描述

在2D图像空间中，基于TpV模型的稀疏性变换可以用双索引的方式表示如下[11]：

(2)

其中：

(3)

其中：j=(m-1)×W+n，m=1,…,H，n=1,…,W，H和W分别代表重建图像x的宽度和高度；p为弹性的范数模型控制参数，通过p值调整可以实现范数模型在L0范数与L2范数之间选择。

压缩感知理论指出[24]：重建信号的稀疏性表示越好，越有利于精确地恢复原始信号。为了应用TpV模型增强重建图像的稀疏表征，并建立动态优化模型。现考虑如下一族双曲正切函数集，其具体表示形式如下：

(4)

其对应的范数求和形式如下：

(5)

其中：σ为动态调整参数，用于控制图像表征和相邻迭代过程的近似性。

如图1所示，当σ参数取值较大时，对重建图像的稀疏性表征较差，重建结果容易导致过平滑，图像纹理细节模糊；当σ参数取值较小时，对重建图像的稀疏性表征较好，重建结果容易陷入局部极值，图像噪声水平显著。基于上述双曲正切函数集构造TpV动态复合函数模型，其具体形式如下：

图1 双曲正切函数曲线动态演化图Fig.1 Dynamic evolution of hyperbolic tangent function curve

(6)

其中：p和σ均为动态优化参数，对式(6)分析，当p趋近于2，σ取较大值时，R(x)接近于L2范数，对于伪影和噪声抑制效果较好；当σ趋近于0，p取较小值时，R(x)接近于理想的L0范数，对图像纹理细节的保护效果较好。因此通过对p和σ的合理选择从而在保持图像边缘和细节的同时能有效克服TV正则项存在的阶梯效应问题。

为了使上述的正则化估计同时表现出较好的纹理保护和噪声抑制性能，本文提出动态最优化过程，即初始迭代优化时，采用较大的p和σ值以更好地抑制伪影和噪声，接下来的优化过程中依次降低p和σ值以更好地保存图像纹理信息。

实验中，为了保证算法初始时具有较好的伪影和噪声抑制效果，初始设定：σ0=0.8，p0=1.2；为了保证相邻两次迭代的近似性，相邻迭代的递减因子Δσ和Δp分别设定为0.016和0.01；同时为了避免算法收敛于局部极值，实验中设定：σmin≥0.1，pmin≥0.8。

2.2 交替最小化DO-TpV重建算法

如式(1)为典型的统计重建目标函数，该目标函数由对数似然数据保真项和正则化项两部分构成，本文采用交替最小化的方法对其进行优化求解，具体算法流程如下：

步骤一、最小化对数似然数据保真项：采用典型的可分离抛物面算法[10]进行优化，其迭代更新公式为：

(7)

式中：t表示迭代次数，x0为初始迭代图像。

步骤二、式(6)的动态最优化求解：采用梯度下降法对其进行优化。其中式(6)的梯度公式可以表示为：

(8)

其中：TpV模型的梯度公式可表示为：

(9)

式中：ζ为扰动约束因子，通常取ζ=1×10-8。

基于上述的梯度公式，采用梯度下降法进行迭代更新，其计算公式为：

(10)

上述算法具体实施中，首先进行参数x、σ、p的初始化，接下来重复执行步骤一与步骤二直至满足预先设定的迭代终止条件。下文实验中将迭代终止条件设定为t不大于预先设定的迭代总次数Tmax或相邻两次迭代误差小于设定阈值err，如果t>Tmax，则迭代终止；否则，t=t+1 ，更新σ和p值，并启动下一轮迭代。本文提出的基于分数阶全变分动态优化算法(DO-TpV)的CT低剂量成像重建流程如下所示：

输入:A,l^,β,σ,p,err,初始化x0,σ0,p0,Δσ,Δp,t=1;

输出:重建图像x.While(t‰Tmaxor‖xt-xt-1‖/‖xt‖≥err){步骤一:优化SIR过程xtj=xt-1j-∑Ndi=1ai,jyi[Axt-1]i-^li ∑Ndi=1ai,jyi∑Jk=1ai,k ,j=1,...,J;步骤二:分数阶全变分动态模型优化(TpV-DO)d=‖xt-xt-1‖2;dx=∂R(x)∂xj;xtj=xtj-β·d·dx;步骤三:参数更新t=t+1;ifσ≥0.1,p≥0.8σ=σ-Δσ,p=p-Δp.}

3 实验与结果

为了验证本文算法在噪声抑制和软组织细节保护方面的有效性，分别采用数值仿真实验和动物扫描实验对本文提出的CT低剂量成像算法(DO-TpV)进行实验验证，并对其与解析重建算法(FBP)、TpV算法、字典学习(DL)、ARWTV算法进行了比较。其中，算法的优化参数设定：本文所提算法采用2.2节所述方式；TpV算法正则化参数β=0.2，动态优化参数p=0.8～1.2；ARWTV算法正则化参数β=0.2，变全次数设定为5；字典学习(DL)算法正则化参数β=0.2，图像块尺寸设定为7×7，过完备字典的稀疏度设定为4。

3.1 数值模拟实验

本部分实验采用数值模型生成圆轨迹扇形束扫描投影数据，投影参数设定如下：X射线源到旋转中心的距离为54.1 cm，探测器到旋转中心的距离为40.8 cm，CT图像的像素数为512×512，尺寸大小为20 cm×20 cm，探测器数为642，相邻探测器间距为0.672 mm。扫描范围为0°～360°，步长为1°和2°，以生成180°和360°两种采样角度数下的投影数据。实验中采用Patrick等[10]提出的Poisson分布模拟光子探测的随机过程，其中：实验总光子数设定为8.0×105。通过对探测器记录信号进行对数化处理后，获得含噪声的投影数据。分别使用FBP算法、TpV算法、字典学习(DL)算法、ARWTV算法和DO-TpV算法进行CT图像重建。图2和图3分别给出了数值模型在180°和360°两种采样数下的实验结果，窗口显示范围[0.185,0.205]。

图3 数值模型实验结果(360采样数)Fig.3 Numerical model results (360 views)

如图2和图3所示，在低剂量扫描条件下，FBP算法重建图像中出现严重的伪影和噪声，且软组织细节完全被噪声和伪影遮盖，总体图像质量极差。对比解析算法(FBP)结果，4种迭代算法重建图像的质量有了明显提升。对比4种迭代算法结果，ARWTV算法取得了很好的伪影和噪声抑制效果，但重建图像软组织细节模糊严重，如图2(e)红色箭头所示的放大区域。相比ARWTV算法，字典学习算法(DL)和TpV算法在软组织细节保护方面有了显著的提升，但字典学习(DL)算法的部分软组织产生变形模糊，如图2(d)红色箭头所示，且字典学习(DL)从整体看尤其在180°较低采样数下会存在条纹状伪影，而对于TpV算法部分软组织细节结构不清晰，如图2(c)红色箭头所示的放大区域。由图2和图3可以看出，DO-TpV算法在软组织细节保护方面表现最好，且取得较好的伪影和噪声抑制效果，总体图像质量的提升显著。

为了评价不同算法的重建结果对高对比度边缘细节的保护能力，本文提取图2红色实线所示位置，并采用剖面密度曲线显示TpV算法、字典学习(DL)算法、ARWTV算法和DO-TpV算法的重建结果与金标准图像的接近程度。图4实验结果表明，与其他算法的重建结果相比，DO-TpV算法的重建结果更好地拟合了金标准图像的密度曲线。

图4 重建结果的剖面灰度值比较图(图2)Fig.4 Line profiles of reconstructed results of Fig.2

本文采用归一化平均绝对偏差(Normalized mean absolute deviation,NMAD)和信噪比(Signal noise ratio,SNR)2个标准定量评价重建结果，NMAD和SNR的定义如下：

(11)

(12)

由表 2 可以看出，DO-TpV算法的NMAD小于FBP算法、TpV算法、字典学习(DL)算法、ARWTV算法的NMAD，表明DO-TpV算法重建的图像与理想图像更接近。在180°视角和360°视角的数值模型实验中，DO-TpV算法重建图像的信噪比比FBP算法、TpV算法、字典学习(DL)算法、ARWTV算法要高，表明DO-TpV算法能更好地保护图像纹理信息和抑制图像噪声，重建图像的质量更高。

表2 数值模型定量分析Tab.2 Quantitative evaluation of the numerical results

3.2 动物数据实验

为了进一步验证DO-TpV算法的性能，采用真实CT系统扫描小鼠来生成投影数据。系统扫描参数设定如下：X光管电压50 kVp，X光管电流1 mA，曝光时间0.466 9 s，X射线源到扫描中心距离为22.188 cm，X射线源到探测器中心的距离为65.85 cm，平板探测器阵列数1 536×880，探测单元尺寸0.15 mm×0.15 mm，扫描间隔为2°，共采集180个等间距投影图像。分别使用FBP算法、TpV算法、字典学习(DL)算法、ARWTV算法和DO-TpV算法对中心片层投影数据进行重建实验。重建像素大小为512×512，重建像素单元尺寸为0.1 mm×0.1 mm。图5给出小鼠扫描数据重建结果，窗口显示范围为[-0.09，1.19]。

图5 小鼠扫描数据重建结果Fig.5 Reconstructed results of scanned mouse datasets

从图5可以看出，低剂量扫描协议下，解析算法(FBP)重建图像中出现严重的伪影和噪声，且软组织边缘细节极其模糊，如图5(b)放大区域所示。对比解析算法重建结果，4种迭代算法重建结果在抑制伪影和噪声，保护软组织结构方面取得了明显改善。对比4种迭代算法结果，ARWTV算法在抑制伪影和噪声方面效果显著，但重建图像部分细节模糊不清晰，如图5(c)放大区域红色箭头所示位置。相比ARWTV算法，TpV算法能够重建出图像的边缘细节，但重建结果中存在一些散点噪声，如图5(c)放大区域所示。而对于字典学习(DL)方法，伪影和噪声得到了一定程度的抑制，但结果仍然存在由于过完备字典而导致的不规则伪影和噪声，如图5(d)放大区域所示。对比ARWTV算法和TpV算法的重建图像，DO-TpV算法重建的图像细节更为清晰，且图像的噪声水平也显著降低，总体图像质量最佳。

4 结 论

本文提出一种基于动态优化的CT低剂量成像算法。该算法通过双曲正切函数集构造分数阶变分的动态复合函数模型，使基于分数阶变分的重建逐步接近理想的L0范数重建。通过动态最优化过程逐次增强当前重建图像的稀疏表征，改善了DO-TpV算法在图像纹理细节保护和噪声抑制方面的性能，提升了重建图像的质量。数值模型实验结果表明，在180个采样角度下，文中算法重建图像的信噪比分别比FBP算法、TpV算法、字典学习(DL)方法、ARWTV算法高出29.51，8.03，9.15，6.81 dB。动物数据实验结果表明，DO-TpV算法重建结果有效抑制了噪声和伪影，清晰重建出小鼠软组织细节，极大地提高了低剂量采集数据重建图像的质量。