探测器故障下医用CT图像重建的容错算法开发

董 建，毕丹阳，杨耿煌，秦转萍

（1.天津职业技术师范大学天津市信息传感与智能控制重点实验室，天津 300222；2.唐山工业职业技术学院自动化工程系，唐山 063299）

计算机断层成像（computed tomography，CT）技术具有适用面广、分辨率高、对比度好、价格相对低廉等优势，已成为使用最广泛的无创临床成像技术之一[1-2]。自1979年我国引入第一台医用CT 以来，CT 设备已经在我国各级医院工作了40年，早期投入使用的设备正相继进入设备老化和故障期[3-4]。X 线探测器是CT设备的核心部件，也是故障高发部件。探测器故障会导致CT 图像出现环状伪影，严重影响正常组织结构和病变区域的刻画[5-6]。而现阶段我国CT 研发及设备维护技术团队已经无力应对各医疗机构集中出现的设备故障问题，研发具有容错机制的CT 图像重建算法已经迫在眉睫。据本团队了解，还未发现国内做相关研究的报告，本文提出的探测器故障下的CT 图像重建的容错算法开发尚属解决此类问题的首例研究。

本文提出一种在未知探测器受损数量和具体位置信息的故障条件下，依然能够精准重建的CT 图像重建算法。传统的迭代式重建算法采用正向投影与实测投影数据的最小二乘差，即L2 范数差作为最小化评价函数[7-8]。众所周知，L2 范数对实测投影数据的异常值尤其敏感，这决定了此传统方案无法避免CT 图像中环状伪影的出现。而一种能有效预测异常值位置并将其精准滤除的方案，即采用L1 范数为最小化评价函数。基于这一想法，开发了一种类似于代数式重建法（algebraic reconstruction technique，ART）的行处理型（row-action-type）迭代式重建算法[9-10]，该算法由凸优化领域的近端分裂理论（proximal splitting）推导构建而成[11-12]。此外，通过向L1 范数追加全变分（total variation，TV）惩罚项[13-16]，新算法在去除环状伪影方面更具鲁棒性。在经典Shepp-Logan 模型重建实验以及实际临床CT 图像的重建实验中，新算法较传统L2 范数重建算法表现出显著优越性。

1 算法介绍

1.1 问题定义

CT 设备获取人体断面图像的扫描过程如图1所示。

图1 CT 设备获取人体断面图像的扫描过程

重建一张人体断层图像，CT 设备需要对该断面进行1 000～1 500 次X 射线扫描。在X 线的某一通路上，令X 线源发出的X 线的强度为I0，探测器接收到的被人体部分吸收后的X 线的强度为I1。在不考虑噪声影响的前提下，I1与I0存在如下等式关系

式中：b 为投影数据；γ 为探测器对X 线的灵敏度。探测器正常工作条件下γ = 1，若探测器出现故障，则γ→0。此时投影数据b 为

即探测器故障情况下，测得的投影数据值远小于探测器正常工作条件下的测得值。重建图像中环状伪影的产生机理阐述如下：在图像重建的反投影这一关键步骤中，受损害的投影数据使得该通路上的各像素由于得不到有效恢复而取值大幅减小。此外，根据重建规则，反投影运算会在各角度方向依次进行，进而使得无法恢复正常取值的像素在某一路径进行重合。180°的反投影运算使得该路径逐渐演变为半圆形结构，而且重合次数足够多的像素取值偏大，重合次数少的像素取值偏小。该半圆形结构即对应重建图像中由相邻的亮区和暗区所组成的环状伪影。完整、1 个和2 个探测器故障的投影数据与相应滤波反投影算法（filtered back-projection，FBP）的重建图像如图2 所示。

图2 完整、1 个和2 个探测器故障的投影数据与相应FBP 算法的重建图像

探测器无故障和随机性出现1 个和2 个探测器故障时，由一般商业CT 搭载的FBP 算法重建了经典Shepp-Logan 模型。由图2 可知，完好投影数据可重建出理想图像，探测器故障可引起投影数据的某一列数据异常，进而导致重建图像中的环状伪影。伪影出现个数与出现故障探测器个数一致。

CT 图像重建即利用重建算法将投影数据转化为CT 图像的过程。迭代式重建算法可将问题定义为求解多元线性方程组 Ax = b 的问题，此处 x =（x1，x2，…，xJ）T为J 维未知量即待重建的CT 图像，b=（b1，b2，…，bI）T为I 维已知量即实测投影数据，A={aij}为I×J 维系统矩阵[17-18]。本文研究I≥J 的超定方程组求解问题。迭代式重建领域的标准评价函数为L2 范数差

本研究将传统L2 范数置换为L1 范数，如式（4），并通过将f（x）最小化而推导出行处理型迭代式重建算法。

2 种范数差展开式表示为

1.2 近端分裂理论

在实施算法推导前，有必要先对其基础内容的近端分裂理论及其对算法推导的重大贡献做简要介绍。首先考虑以下凸优化问题

假定式中f（x）为可能不可微分的下半连续凸函数。近端算子法（proximal operator）为求解此类最优化问题提供了良好解决方案。对应于f（x）的近端算子表达式为

式中：α 为步长参数。

近端算子具备2 种优秀属性。

（1）它可针对任何一个下半连续凸函数进行定义，即使该函数不可微（如L1 范数）；

（2）在步长参数α ＞0 的前提下，满足x=proxαf（g）的定点x 与使凸函数f（x）取最小值时的变量x 相一致。

这2 个性质决定了可以应用如下迭代式求得凸函数f（x）取得最小值时的变量x

式中：k 为迭代次数。

该迭代算法叫做近端最小化算法，它为不可微凸函数的最优化问题提供了强有力的解决方案。

凸函数f（x）可拆分成若干分量函数fi（x），即其中fi（x）（i=1，2，…，I）均为可能不可微但满足下半连续的凸函数。此时，假定凸函数f（x）的近端算子求解难度巨大，而其各个分量函数fi（x）的近端算子能够由低计算成本求得。近端分裂理论能够为该条件下的凸优化问题提供解决方案，它的处理方式为在求解由x(k)到x(k+1)的过程中，进一步将问题分解为按照i=1，2，…，I 的顺序进行各分量函数所对应的近端算子的计算，即

当 i=I 时，x(k+1)=x(k，I+1)。

式中：k 为主循环变量；i 为子循环变量。

步长参数α(k)值会随k 值改变。为了使所求解序列x(k)收敛于函数f（x）取最小值时的变量x，α(k)需满足如下条件[12]

近端分裂理论为本文L1 范数重建算法的构建提供了理论依据。

1.3 L1范数重建算法描述

式中：ai为系统矩阵A 所对应的第i 行分量。进而求解各分量函数所对应的近端算子，即求解以下最优化问题

该近端算子的计算可由拉格朗日乘子技术简单求得。首先，引入替代变量则以上问题转变为约束极小化问题

对应于式（13）的拉格朗日函数可定义为

式中：λ 为拉格朗日乘子，也称作对偶变量。

对应于式（14）的对偶问题可定义为

通过对该对偶问题的求解，可得到未知解x 的更新公式为

同时可得到对偶变量的表达式及取值范围

可将λ 分情况表示为

将式（18）代入式（16），即可得到 x(k,i)的更新解x(k,i+1)=proxα(k)fi（x(k,i)）。由于步长参数α(k)需满足随循环变量k 而递减的条件，故将其设定为

式中：α0＞ 0，є ＞ 0 为提前给定的正数常量。

至此，L1 范数重建算法已经构建完成。

1.4 L1范数重建算法

步骤0设定步长参量α0和є 的初始值，为变量x 取初值 x(0，1)，以 k=0，1，2，…的顺序执行步骤 1 到步骤3；

步骤1计算步长参数

步骤2以i=1，2，…，I 的顺序执行以下操作：

①计算λ

②更新x

步骤3x(k+1,1)=x(k,I+1)

为了得到更好的图像重建效果，进一步向L1 范数的评价函数中追加了全变分惩罚项‖x‖TV，定义如下

式中：xm，n为在一幅二维图像中坐标为（m，n）的像素的像素值。全变分模型是一个依靠梯度下降流对图像进行平滑的各向异性的模型，能对图像起到良好的去噪和平滑作用[13-16]。因此，L1-TV 重建的最小化评价函数为

式中：β 为可调节全变分作用强弱的超参数。

2 实验结果及分析

2.1 少数探测器故障下的图像重建

本文对经典Shepp-Logan 数字模型进行重建实验。该数字模型所含像素数为256×256，像素值的分布范围为[0，4.5]。图像重建领域的常用做法是获取与重建图像的像素数等量的投影数据量，因此本文以[0°，180°]和 180°/256 的角度间隔进行投影数据采集，投影数据的像素数同样为256×256。

假定出故障的探测器个数为1 个和2 个，并且假定出现故障的探测器的位置具有随机性。分别用L2范数重建算法、L1 范数重建算法和L1-TV 重建算法对含有异常值的投影数据进行重建实验。经实验验证，对于本实验的Shepp-Logan 数据模型，在各参数取以下数值时可得到最佳重建效果。在各项实验中α0=0.03，є=50，L1-TV 实验中 β =260。各项实验的迭代次数k 设定为100。在不同探测器故障个数条件下各算法的实验结果如图3所示。其中，图3（a）、（c）、（e）为 1 个探测器故障时的重建结果，图3（b）、（d）、（f）为2 个探测器故障时的重建结果。由结果图像可知，L2范数重建的图像仍含有严重的环状伪影，并且在半圆环伪影的起点和终点处二次产生了直线型放射状伪影。L2 范数的重建结果与FBP 算法的重建结果的图像效果相似。L1 范数的重建效果较L2 范数的重建效果有了显著提高，环状伪影得以有效去除。然而，伪影处的黑色瘢痕仍有残留，并且图像全体可见微弱的直线型伪影。L1-TV 算法得到了最好的重建效果，图像中的环状伪影得以完全去除，并且各组织结构的灰度得到了有效的平滑。

2.2 多数探测器故障下的图像重建

本文进一步设计了条件严苛的验证实验，即假定探测器的故障个数为9 个。各项实验中的参数取值不变。9 个探测器故障条件下各算法的无噪声影响和有噪声影响的重建结果如图4所示。图4（a）、（c）、（e）为在无噪声影响下用各种范数进行图像重建的实验结果。由结果图像可知，严重的环状伪影分布于整幅L2范数的重建图像，并且在环状伪影的起点和终点处存在显著的直线型伪影。图像的某些细节信息被伪影覆盖，并且高灰度值的环状伪影降低了图像原有信息的对比度。在L1 范数的重建图像中，高灰度值的环状伪影被有效去除，并且在环状伪影的起点和终点处的直线型伪影也被完全去除。图像的细节信息得到了较好的保存。然而，在L1 的结果图像中，仍有暗黑色环状伪影瘢痕残留。L1-TV 算法得到了最佳实验效果，伪影被完全去除，图像细节保存完好，并且各组织结构均匀一致。鉴于CT 图像采集过程中会受到电子噪声的影响，为了使该实验更接近实际运行机制，在投影数据中加入了均值为0、方差值为6.0 的高斯噪声和发射光子数为5.0×105条件下的泊松噪声。各项实验中其余参数取值不变，为控制统计噪声将循环次数k设定为 200。图4（b）、（d）、（f）为在 2 种噪声影响下用各种范数进行图像重建的实验结果。在L2 范数的重建图像中充满了亮暗并存的环状伪影以及统计噪声，环状伪影与物体边缘交界处可见部分组织结构变形。L1 范数重建充分去除了环状伪影的明亮区域和部分统计噪声。L2 范数重建中的部分组织结构变形得以恢复。L1-TV 算法有效去除了环状伪影和统计噪声，然而由于为达到去除统计噪声目的而进行的过度平滑使得图像微小细节结构有所消失。

图3 不同数量探测器故障条件下各算法的实验结果

图4 9 个探测器故障条件下各算法的无噪声影响和有噪声影响的重建结果

2.3 临床CT图像的重建

为进一步验证算法有效性，选取了临床中常用的上腹部CT 图像进行了重建实验。该图像的像素数为512×512，像素间距为0.625 mm，因此图像的采集范围为32 cm×32 cm。各像素的取值为对应区域的X 线衰减系数值，像素值分布范围为[0，0.387]。采集的投影数据像素数为512×512。各实验参数设置如下：α0=0.002，є=50，L1-TV 实验中 β =8 600（β 值越大，图像的平滑效果越小）。含3 个探测器故障的上腹部CT 断层图像的重建结果如图5 所示，实际图像的重建效果与上述数字模型重建效果保持一致。L2 范数重建结果中有明显的环状伪影出现，而L1 范数重建有效去除了环状伪影的明亮区域，L1-TV 重建得到了良好质量的重建图像。

图5 含3 个探测器故障的上腹部CT 断层图像的重建结果

实验结果表明，在单个及多数探测器故障条件下L2 范数重建不能消除图像中的环状伪影，而L1 范数重建能有效去除高灰度值的环状伪影，L1-TV 重建算法能够得到良好的图像重建效果。在受高斯噪声和泊松噪声影响的条件下，L1 算法依然展现出对环状伪影去除的有效性。深究3 种算法的不同发现，L1 算法较L2 算法具有显著优越性的原因在于图像更新过程具备异常值滤除程序。式（16）中，λ 的取值被限定在的取值被限定在一定合理区间，该过程即决定了若投影数据bi存在异常值则会被滤除的运行机制。而L1-TV 算法较L1 算法具有优越性，原因在于全变分惩罚项能够对图像起到有效的平滑和去噪作用。

3 结 语

本文研究了在医用CT 设备的核心部件X 线探测器部分损坏的条件下，进行具有容错机制的图像重建算法的开发问题。以经典Shepp-Logan 数字模型作为重建对象，设计了单个及多个X 线探测器出现故障情况下的重建实验。此外，还讨论了在电子噪声影响下的重建情况以及实际的上腹部CT 图像的重建情况。分别对比了传统L2 范数重建算法和本文提出的L1范数重建和L1-TV 重建算法。实验结果表明，L2 范数重建不具备去除探测器故障引起的环状伪影的功能，其重建效果接近于当前医用商业CT 所搭载的FBP 算法的重建效果。而L1 范数重建算法较传统L2 范数重建在去除探测器故障引起的环状伪影方面具有显著优越性，能有效去除CT 图像中光亮的环状伪影，然而该算法在去除暗黑色环状伪影方面还存在不足。增强算法L1-TV 算法有效弥补了L1 算法的不足，具备完全去除环状伪影和有效保护图像细节和对比度的能力。本文提出的探测器故障下医用CT 图像重建的容错算法具有较强的容错能力，能够实现多数探测器同时出现故障的情况下的精准图像重建。因此，该技术为我国各临床机构中已经进入老化和故障期的CT 器械能够继续服役提供了有效的技术支撑。