X射线CT为医学诊断和治疗提供了人体主要的解剖和病理信息^[1]。近年来，越来越多的临床医生通过CT图像对患者的病症进行更准确的诊断，甚至借助CT图像的引导对病灶进行穿刺活检，对人体管道、体腔或器官组织内病理性积液、血肿、或脓肿进行穿刺引流达到治疗的效果。虽然CT技术的发展提高了人类的健康水平，但是高剂量和累积剂量的X射线照射会对人体造成不可逆转的放射性损害，甚至会诱发癌症和遗传性疾病^[2-5]，这蕴藏着低剂量CT(LDCT)技术研究和发展的新契机^[6-8]。目前，在临床上降低CT诊断辐射剂量最行之有效的方法是通过调节mA(毫安)或mAs(毫安秒)的设置来降低X射线管电流，但是X射线管电流的降低势必导致重建后的LDCT图像被斑点噪声和条纹伪影所覆盖，严重退化的LDCT图像已不能满足临床诊断需求，甚至会导致漏诊、误诊，这与医学CT技术发展的初衷相背而驰。

为使获取的LDCT图像满足临床诊断需求，有效提高LDCT图像的质量成为当前的研究热点^[9-13]。学者们探索了许多策略^[14-16]来提高LDCT图像的质量，目前直接对重建的LDCT图像实现消除噪声和抑制条形伪影的后处理方法是不依赖投影数据的，易于集成到大多数CT系统中，为LDCT成像技术的临床应用发展注入新的活力。鉴于难以对LDCT图像中分布不均匀的噪声和条形伪影进行统计建模，目前大多数后处理技术往往依赖于更多的启发式(heuristic)方法。近20年来，变分法在图像去噪领域独树一帜，具有良好的保边去噪性能^[17-19]。基于此，研究者们将变分法应用到医学成像领域并取得新成效。受文献[17]启发，本文提出了一种基于可变阶变分模型的后处理技术，用于改善LDCT图像的质量，该模型根据图像的局部特征自适应地选择正则化算子的阶数，在抑制斑点噪声和条形伪影的同时可以保留具有临床诊断意义的结构特征。

1 可变阶变分模型

基于变分法的图像降噪是求解关于图像的能量泛函最优解的问题，而图像降噪是一个不适定的反问题，所以常常运用正则化理论将不适定问题转变为适定性问题来求最优解。

首先提出如下凸结合变分正则化模型：

式中：Ω为图像的定义域；u为降噪后的低剂量图像；u₀为带噪声和伪条纹的低剂量图像；▽为梯度算子；λ为正则化参数；θ∈[0, 1]为变分阶数控制参数；等式右边第1项为保真项，第2项为对图像具有平滑作用的正则项。此外，本文模型中Hessian矩阵(即)记为：▽²u=。显然，合理地选择θ，本文模型有以下几种情况：

1) 当θ=1时，本文模型可以改写为

此时模型退化为正则化算子为一阶变分的总变分(TV)模型，该模型具有良好的保边性能，但会产生阶梯效应。

2) 当θ=0时，本文模型可以改写为

此时模型退化为正则化算子为二阶变分的有界Hessian(BH)模型，该模型具有良好的平滑性能，但是保边能力欠佳。

3) 当θ∈(0, 1)时，本文模型类似于TVBH模型，兼顾了一阶变分和二阶变分，是对TV模型和BH模型的融合。

综上所述，参数θ的选择决定着本文模型的滤波形式与滤波性能。由于参数常量θ通常通过大量实验来试凑或凭经验取得最佳值，是对全局内容的粗略评估，忽略了图像的局部特征。

接下来考量图像的局部特征, 用边缘扩散函数θ(z)来代替常量θ改善本文模型的自适应性，提出了可变阶变分模型，优化的模型形式为

式中：k为边缘阈值参数。

为了更好地保持医学图像中所包含的边缘细节信息，关于边缘扩散函数θ(z)的特征检测因子z的构造，融合了梯度和局部熵2种特征检测算子。

1) 图像的梯度

梯度表征了图像灰度值的变化大小和方向，因此常用梯度模值来区分图像的边缘区域和非边缘区域。边缘区域梯度较大，平坦区域梯度较小，但是有些细节信息的梯度与平坦区域相差不大，噪声点处的梯度甚至比边缘的梯度还要大，这样梯度边缘检测算子可能会对图像中细节丰富的弱边缘区域和强噪声点造成误判，导致处理后的图像细节信息损失或者去噪不彻底。

2) 图像的局部熵

局部熵表征了图像中局部区域的像素灰度值变化的剧烈程度，从而可以反映图像中所包含信息的丰富程度。一副大小为M×N的灰度图像f的熵值定义为

式中:f(i, j)为位于图像(i, j)点处像素的灰度值；p_{i, j}表示(i, j)点处像素的灰度值在大小为M×N的局部邻域中的分布概率；H为图像的局部熵。

根据局部熵可以有效地判定图像的局部特性，在灰度分布复杂的边缘细节区域，局部熵值较大；在灰度分布均匀的平坦区域，局部熵值较小。此外，局部熵具有较强的抗噪能力，独立噪声点对其影响很小。因此，局部熵可被广泛应用于图像处理中。

根据梯度和局部熵的特性，新的特征检测因子z设计为，u_σ为对图像进行参数为σ(通常取为0.5)的高斯平滑。从式(5)可知，在给定阈值参数k的边缘扩散函数的引导下，本文模型可自适应地选择变分阶数。在强边缘处，θ(z)≈1，此时正则化算子为一阶变分，能够在滤除噪声时较好地保护边缘纹理等细节信息；在平坦区域处，θ(z)≈0，此时正则化算子为二阶变分，能够很好地平滑噪声；在弱边缘处，0＜θ(z)＜1，此时正则化算子为混合阶变分，实现较好的保边去噪效果。

2 可变阶变分模型的数值解法 2.1 离散化

为了提高本文提出的可变阶变分模型的运算速率，类似文献[17]采用基于快速傅里叶变换(FFT)的分裂Bregman算法求解。

首先对连续的一阶、二阶和四阶微分算子及散度算子进行离散化。此外，为了进一步提高分裂Bregman算法计算的速度，采用周期性边界条件使FFT适用于分裂Bregman算法。设Ω为尺寸大小为M×N的二维灰度图像区域，图像列和行方向的坐标分别用x和y表示。在像素点(i, j)处沿着坐标x和y方向的一阶向前差分记为

一阶向后差分记为

离散的二阶差分记为

离散的四阶差分形式如式(18)~式(20)所示，其中关于边界的有限差分的离散形式可参考式(12)~式(15)。

离散化的梯度、Hessian矩阵为

最后对于任意p=(p₁ p₂)的一阶散度div，任意的二阶散度div²的离散形式分别记为

2.2 本文模型的分裂Bregman算法

为了求得本文提出的可变阶变分模型的能量泛函的最优值，首先引入辅助变量w、v，Bregman迭代参数b、d和惩罚参数θ₁、θ₂。将本文提出模型的能量泛函转换为以下形式：

式中:w=(w₁ w₂)为一个与函数u的梯度相关的二维函数且为一个与函数u的Hessian矩阵相关的四维向量函数且；b=(b₁ b₂)代表约束的Bregman迭代参数；代表约束的Bregman迭代参数。

接下来固定变量(w, v; b, d)，得到关于u的欧拉方程为

根据式(21)~式(24)可得式(26)的离散形式为

式中:G_{i, j}为式(26)等号右侧部分，离散形式为

在周期边界的条件下，对等式(27)两端分别进行离散傅里叶变换，变换形式为

式中: 为傅里叶变换符号。

式(29)关于离散频率r和s的等价形式为

通过对式(29)进行傅里叶逆变换，求得u的解析解：

式中: 为傅里叶逆变换符号；为复数的实部。

此外，关于变量w的更新，通过固定(u, w; b, d)，获得与w相关的欧拉方程：

式(33)关于w的解析解可以用如下二维广义软阈值公式求得：

采用同样的方法更新变量v，其解析解为

最后更新Bregman迭代参数b、d。

综上所述, 本文提出的可变阶变分模型的基于FFT的分裂Bregman算法的求解流程可以归结如下：

步骤1     初始化：u₀(u₀设定为待改善的LDCT图像)。

步骤2    初始化：(w, v; b, d)=0。

步骤3    设定参数λ、θ₁及θ₂。

步骤4    由d)更新变量u、w及v。

步骤5    更新Bregman迭代参数。

步骤6    更新Bregman迭代参数。

步骤7    若满足迭代终止条件： 10^－6时(uⁿ和u^n－1分别为图像第n次和第n－1次迭代滤波后的结果)，则停止迭代。反之，则转至步骤4，直到满足终止条件。

3 实验结果与分析

为了验证本文提出的可变阶变分模型对医用LDCT图像的去噪性能及在临床应用上的价值，选择实际CT扫描的胸腔体模、临床数据作为实验对象，其尺寸均为512 mm×512 mm。图 1(a)为胸腔体模HDCT图像，其通过一台多排探测器西门子16 CT扫描仪在管电流为240 mAs、管电压为120 kVp(千伏峰值)的高剂量扫描协议下实际扫描一个人体胸腔的仿真体模，获得经反投影重建(Filtered Back Projection，FBP)算法重建后的图像，再利用大尺度邻域非局部伪影抑制(AS-LNLM)算法对重建后的图像质量进行进一步的改善而获得^[20]。图 1(b)为胸腔体模LDCT图像，其是通过将管电流降为30 mAs，在此低剂量的条件下扫描胸腔仿真体模的同一部位，经FBP重建后获得。图 2(a)为胸腔体模HDCT图像的局部放大部分, 图 2(b)为胸腔体模LDCT图像的局部放大部分。图 3(a)为临床胸腔LDCT数据，该数据是在管电流为60 mAs、管电压为120 kVp的低剂量扫描协议下得到的FBP重建后的DICOM格式的低剂量图像。图 4(a)和图 5(a)分别为临床腹腔和盆腔LDCT数据，这2个数据是将辐射剂量降为标准剂量的25%的低剂量扫描协议下得到的FBP重建后的DICOM格式的低剂量图像。

实验中，分别采用TVBH模型^[17]、EWSO模型^[18]、NLM模型^[21]和本文模型对胸腔体模LDCT图像和临床LDCT图像进行处理。本文模型的边缘扩散函数θ(z)中，k是一个重要参数，其决定着要保持边缘的对比度并影响扩散过程。更具体地说，如果k值过大，图像的边缘被过度平滑；如果k值过小，图像中的一些噪声和伪影难以滤除。经大量实验验证，当k=10时，本文模型降噪效果最佳。对于胸腔体模，本文模型的参数设定为：λ=30，θ₁=θ₂=5。对于临床胸腔数据，本文模型的参数设定为:λ=80，θ₁=θ₂=10。对于临床腹腔和盆腔数据，本文模型的参数设定为：λ=50，θ₁=θ₂=10。为了客观说明本文模型的性能，其他对比模型的参数均按取得最优降噪效果设定。图 1~图 5展示了各模型的降噪结果。

观察图 1可知，本文模型对胸腔体模LDCT图像降噪效果在视觉上整体优于其他3种模型。为了更清晰地展示各模型的降噪效果，对图 1进行了局部放大(见图 2)。图 2(c)为TVBH模型的降噪效果，该模型有效去除了噪声和条形伪影，但是过度平滑LDCT图像的一些细节信息会造成临床诊断的误诊。图 2(d)为EWSO模型的处理结果，该模型虽然可以抑制大部分噪声和条形伪影，但是把LDCT图像的组织结构附近的伪影误判为结构信息保留了下来，不能满足临床诊断需求。图 2(e)所示NLM模型的结果较好地保留了边缘信息，但是丢失了一些细小结构，同时也引入了一些新的细小伪影，亦不能满足临床诊断需求。图 2(f)为本文模型的去噪结果，相对于其他3种对比模型，本文模型在抑制LDCT图像的条形伪影和边缘保持方面具有明显优势，其处理结果与HDCT图像最为相似，具有更加理想的视觉效果，为临床诊断带来便利。临床医学图像含有大量的弱边缘信息，相对仿真体模而言，临床医学LDCT图像的伪影去除难度更大。为了进一步验证本文模型在临床上的有效性，图 3~图 5给出了各模型对临床数据LDCT图像的处理结果。根据对比结果可知，TVBH模型很好地滤除了条形伪影，但是严重破坏了LDCT图像中组织结构；EWSO模型和NLM模型的处理结果均含有部分伪影，视觉效果较差；本文模型可以有效滤除LDCT图像中的条形伪影，而且较好地保留了LDCT图像中的边缘组织结构信息，比其他3种模型具有更佳的视觉效果。

此外，为了验证上述实验的有效性和合理性，针对图 1(a)所示的实际胸腔体模采用峰值信噪比(Peak Signal to Noise Ratio, PSNR)、平均结构相似度(Mean Structural Similarity Index Measure, MSSIM)和降噪时间作为衡量各模型去噪性能的客观性能评价指标，如表 1所示。显然，本文模型相比于其他3种模型取得了较高的PSNR值和MSSIM值，这与上述实验的视觉效果对比图是吻合的，说明本文模型处理结果的降噪效果最优，且与对应的HDCT图像的相似度最高。另外，本文模型的降噪时间略高于TVBH模型和EWSO模型，明显低于NLM模型，总的来说耗时较短，便于应用于临床中。针对图 3~图 5中的(a)图所示的临床数据缺乏高质量的参考图像，采用信噪比^[22](Signal to Noise Ratio, SNR)对局部感兴趣区域(ROI)进行质量评价，如表 2所示。可以看出，本文模型对临床LDCT图像处理结果的定量评价均优于其他模型，进一步说明了本文模型处理临床数据的有效性。综上所述，本文模型是一个综合性能较优的去噪模型。

4 结论

本文针对LDCT图像存在严重的斑点噪声和条纹伪影问题，提出了一种可变阶变分模型降噪模型。该模型兼顾了一阶变分模型和二阶变分模型的优点，根据图像的梯度模值和局部熵融合的特征检测因子将图像划分为平坦区域和细节区域，在不同的区域自适应地选择变分阶数，从而实现不同的扩散模式。同时采用基于FFT的分裂Bregman算法求解所提出的变分模型，提高了运算速度。仿真实验表明，本文所提出的可变阶变分模型不但有效地去除了LDCT图像中存在的噪声和条纹伪影，而且很好地保留了图像的边缘细节，明显改善了LDCT图像的质量，获得了与HDCT图相近的视觉效果，适用于临床诊断，具有实际应用价值。