4D-CT在肺癌放射治疗中发挥着重要的作用，它提供了一个全面的高精度放射治疗呼吸运动表征，有助于跟踪肿瘤运动，实施精确放射治疗，减少对正常组织的损伤^[1]。肺4D-CT数据的获取，通常是根据床位和肺体积，将多个自由呼吸3D-CT数据段排序而得。然而，由于CT固有的高剂量照射^[2-3]，故往往只能降低沿纵向（Z轴方向）的采样来减少4D-CT扫描时间以求降低辐射量，从而导致肺4D-CT图像层间分辨率远低于层内分辨率，造成数据显著的各向异性。因此，在对肺各相位3D数据进行冠矢状面观察时，为了获得正确比例的图像，需要根据3D数据的层间分辨率和层内分辨率的比例，沿Z轴方向进行插值放大。常用的插值方法为最近邻或双线性插值法，这些方法都会导致图像的模糊，尤其当层间分辨率与层内分辨率比例差别比较大时。因此，研究如何提高肺4D-CT图像分辨率有着实际意义。

图像超分辨率重建技术是一种提高图像分辨率的有效方法。其原理是利用多帧关于同一场景的有相互位移的低分辨率降质图像来重建高分辨率图像^[4]。本文的主要目的是提高肺4D-CT数据多平面显示图像的质量。肺4D-CT数据提供随呼吸运动变化的肺部低分辨率图像序列，故多平面间不同相位对应图像可以认为是一系列低分辨率图像“帧”，这些图像间信息是互补的，根据这个特征我们可采用超分辨率技术重建清晰的高分辨率图像。超分辨率重建方法已有多年的研究。其可大体划分为两类：频域方法和空域方法^[5]。频域方法是基于傅里叶变换的方法，将图像变换到频域内以消除低分辨率图像的频谱混叠。其算法简单，运行速度快，但是难以处理图像噪声问题，而且只能处理全局整体运动的降质模型^[6]。因此，该方法有很大的局限性，并不常用。目前的研究主要集中在空域方法方面。常用空域超分辨率重建方法有凸集投影法（POCS）、迭代反投影法（IBP）、最大后验法（MAP）等。他们各有优缺点，其中迭代反投影法因为理论简单、收敛快而得到了广泛应用。在文献[7]中，我们提出基于运动估计的迭代反投影的超分辨率重建方法，实现了肺4D-CT数据冠（矢）状面图像分辨率增强。在此方法中，图像间运动场的估计是影响重建精度和速度的主要因素。由于我们采用的是全搜索运动估计法^[8]估计不同帧图像之间的运动场，因此该方法的缺点是速度慢，且只能是固定搜索步长，这样既影响重建速度，也影响重建精度。

本文中我们结合三步搜索法^[9]和光流场方法^[10-12]，实现了快速的图像间亚像素运动估计，而后，采用迭代反投影（IBP）算法^[11]，重建出高分辨率的肺4D-CT图像。实验结果表明，本文基于快速亚像素运动估计的超分辨率技术的重建结果显著优于传统的双线性插值法；并且，本文算法重建结果较文献[7]提出的基于完全搜索块匹配运动估计的超分辨率重建结果好，计算时间也大幅度缩减。

1 方法 1.1 不同相位间肺4D-CT图像运动估计

本文将多相位低分辨率图像作为输入“帧”，首先估计各帧图像间的运动场。具体包含两个步骤：（1）采用三步搜索法，搜索图像块之间的整数像素位移；（2）采用光流法搜索亚像素位移。整个过程叙述如下。

1.1.1 三步搜索法（three step search，TSS）

TSS采用的是一种由粗到细的搜索模式，从搜索窗中心点开始，按一定步长取周围8个点构成每次搜索的点群，然后按照匹配准则进行匹配计算，找到误差最小的匹配块中心点。匹配准则本文仍选择简便有效的绝对误差和（SAD）来求解两图像间的初始运动矢量，SAD定义如下：

$ S\left({{\nu _x}, {\nu _y}} \right)=\sum\nolimits_{i=0}^{N - 1} {\sum\nolimits_{j=0}^{N - 1} {\left| {{I_t}\left({p+i, q+j} \right)- {I_{t - {\Delta _t}}}\left({p+{\nu _x}+i, q+{\nu _y}+j} \right)} \right|} } $

(1)

式中，图像块的大小为N×N，左上角的坐标为(p，q)；运动矢量为(ν_x，ν_y)；I_t(i，j)和I_{t -△t}(i，j)分别为当前帧和参考帧在像素(i，j)处的值。在搜索区域内使ν_x，ν_y)值最小的运动矢量即为当前块的最优运动矢量。

TSS的步骤为：第一步，首先确定一个中心点，设定最大搜索长度，以最大搜索长度的1/2作为步长，将中心点及周围距离相同步长的8个检测点根据匹配准则，找到最小块误差点，如果最小块误差点位于原中心点，则算法结束，否则，进行下一步；第二步，步长减半，在上一步确定的最小块误差点及其周围相同步长的8个检测点中找到最小误差匹配块的中心点；第三部，重复第一部和第二步，直到步长达到搜素精度要求，即得到最佳匹配点。其原理示意图如图 1所示^[9]。

1.1.2 光流法

采用三步搜索法进行初始运动估计之后，所得运动矢量精度为整数，其为图像运动矢量的粗略估计，因此本文还将采用光流法对初始运动矢量进行精度优化，精确到亚像素。

假设运动前后两帧连续图像分别为f(x，y)和g(x，y)，根据经典光流简化模型^[12]，可得式（2）：

$ g\left({x, y} \right)=f\left({x+\Delta {x_s}, y+\Delta {y_s}} \right)\approx f\left({x, y} \right)+\Delta {x_s}\frac{\partial }{{\partial x}}f\left({x, y} \right)+\Delta {y_s}\frac{\partial }{{\partial y}}f\left({x, y} \right) $

(2)

该式近似为一阶泰勒级数展开式。

对（2）式进行最小化求解可得最优位移矢量：

$ \mathop {\min \;\;imize\Phi \left({\Delta {x_s}, \Delta {y_s}} \right)}\limits_{\Delta {x_s}, \Delta {y_s}} $

(3)

其中：

$ \Phi \left({\Delta {x_s}, \Delta {y_s}} \right)={\sum\limits_{x, y} {\left({g\left({x, y} \right)- f\left({x, y} \right)- \Delta {x_s}\frac{\partial }{{\partial x}}f\left({x, y} \right)- \Delta {y_s}\frac{\partial }{{\partial y}}f\left({x, y} \right)} \right)} ^2} $

(4)

式（4）可看作线性最小二乘求解问题，可将目标函数的偏导数置0得到最优值△x _s，△y_s。因此，可有如下方程：

$ \frac{{\partial \Phi }}{{\partial \Delta {x_s}}}=0和\frac{{\partial \Phi }}{{\partial \Delta {y_s}}}=0 $

(5)

求解方程（5），可得最优解△x _s，△y_s。

设三步搜索法得到的初始运动矢量为△x_I，△yI, 则最终得到的总运动矢量△x，△y为：

$ \Delta x=\Delta {x_I}+\Delta {x_s}, \Delta y=\Delta {y_I}+\Delta {y_s} $

(6)

1.2 迭代反投影（IBP）超分辨率重建

IBP算法是由Irani和Peleg ^[13]率先提出，其基本思想是：首先，假设一个初始高分辨率图像H ⁽⁰⁾，然后，根据退化模型^[14]模拟成像过程得到一个低分辨图像的集合{L_k⁽ⁿ⁾}，k=1，...，K, K表示原始序列低分辨率图像的数量，n为迭代次数。如果H ⁽⁰⁾是正确的高分辨率图像，则得到的这一集合就应该和已知的低分辨率图像序列集合{ L_k}相等。如果H⁽⁰⁾不是理想的高分辨率图像，则{L_k}与{L_k⁽ⁿ⁾}将产生差值{L_k-L_k⁽ⁿ⁾}。将此差值返回给H ⁽⁰⁾并更新之，得到一个相对较接近的正确结果H ⁽¹⁾。反复迭代该过程直至收敛，则重建结束，误差计算为：

$ {e^{\left(n \right)}}=\sqrt {\frac{1}{K}\sum\nolimits_{k=1}^K {\left\| {{L_k} - L_k^{\left(n \right)}} \right\|_2^2} } $

(7)

在第n次迭代过程中，L_k⁽ⁿ⁾可表示为：

$ L_k^{\left(n \right)}=\left({{T_k}\left({{H^{\left(n \right)}}} \right)h} \right)\downarrow s $

(8)

式中，T_k表示从H到L_k的二维几何变换，且可逆；↓s是下采样算子，h为高斯模糊算子。

重建的高分辨率图像的更新过程表示为：

$ {H^{\left({n+1} \right)}}={H^{\left(n \right)}}+\frac{1}{K}\sum\nolimits_{k=1}^K {T_{^k}^{ - 1}\left({\left({\left({{L_k} - L_k^{\left(n \right)}} \right)\uparrow s} \right)p} \right)} $

(9)

式中，↑ s表示上采样算子；p表示背投影算子，可设定为h的逆。

本文中选定一幅需要重建的低分辨率肺4D-CT图像，插值重建为初始高分辨率图像H ⁽⁰⁾，其他相位图像作为已知低分辨率图像序列，根据获得的运动场，采用IBP方法重建出高分辨率4D-CT图像。

1.3 超分辨率重建结果量化评价

本实验采用一个公共可用的数据集^[15]，该数据集由10组肺4D-CT数据组成，每组数据包含10个相位。各组图像数据层内分辨率为0.97~1.16 mm，层间分辨率均为2.5 mm。实验针对10组数据，在不同相位中选取冠矢状面低分辨率图像，然后分别采用双线性插值，文献[5]算法和本文算法进行图像重建。

本文采用图像平均梯度^[16]这个量化指标来评价重建结果。平均梯度能够反映出图像细微反差的程度，平均梯度越大，表明影像越清晰、反差越好、细节保持的越好。其定义为：

$ \nabla \bar f=\frac{1}{{\left({M - 1} \right)\left({N - 1} \right)}}{\rm{ \times }}\sum\nolimits_{i=1}^{M - 1} {\sum\nolimits_{j=1}^{N - 1} {\sqrt {\frac{{\nabla _i^2f\left({i, j} \right)+\nabla _i^2f\left({i, j} \right)}}{2}} } } $

(10)

式中，f(i，j)，▽_i f(i，j)和▽_j f(i，j)分别为像素点灰度以及其在行、列方向上的梯度；M和N分别为图像的行数和列数。

2 结果 2.1 超分辨率重建结果视觉评价

图 2显示了1例典型的冠（矢）状面图像重建结果。从左至右，分别为双线性插值，基于全搜索运动估计超分辨率重建（文献[5]）和本文基于快速亚像素运动估计超分辨率的重建结果。第1、3行分别为冠状面和矢状面图像，为了更清楚地进行效果对比，我们分别将红色方框内的区域进行了放大，置于原图下方显示。从实验结果可看出，本文算法获得的超分辨率重建结果比双线性插值法有很大改进，图像的分辨率有较明显的提高，从局部放大图像来看，肺实质中的血管和周边组织的亮度和清晰度更有明显增强。同时，本文结果相比于基于全搜索的超分辨率重建结果，清晰度亦有所改善。这是因为全搜索只能设定固定步长，为了节约时间，步长不能设置太小，故影响了它的运动估计精度。另外，在运算时间方面，基于全搜索的超分辨重建时间大约40 min，而本文方法只需约2 min。时间缩短了近20倍，重建效率明显改进。

2.2 超分辨率重建结果量化评价

平均梯度既反映了图像中的微小细节反差与纹理变化特征，也反映了图像的清晰度。本文根据图像平均梯度量化指标，针对10组数据重建的图像进行了量化评价。根据式（10）对10组图像的双线性插值，基于全搜索运动估计重建（文献[5]）和本文方法的重建结果分别进行平均梯度计算，结果如表 1所示。可见，本文算法较双线性插值方法的平均梯度值显著增大，且与基于全搜索重建结果相比，亦有一定的提高，图像更清晰，细节保持更好。

2.3 参数讨论

我们选取大小为16*16的图像块进行匹配运算，其主要原因是：在图像块搜索过程中，常用图像块大小有8*8像素、16*16像素和32*32像素等。若选取大的图像块，将会明显增加搜索计算。如：32*32的块与16*16的块相比，其块相似度匹配的运算量将增加3倍。而如果图像块过小，由于肺纹理比较稀疏，有时会导致块包含的结构特征较少，影响匹配搜索精度。因此，本文设定图像块大小为16*16像素，同时考虑肺的运动幅度不大（像素位移一般为10个像素以内），故我们将搜索范围设定为16*16像素，也就是在目标像素点周围16*16像素的范围内搜索。

在公式（4）计算中，数值解法将影响结果。我们采用实验图像，比较了5点差分和中心差分方法这两种较典型方法对运动矢量估计的影响，并用均方根误差（RMSE）评价两种方法估计的运动场差异，其结果为RMSE=0.1504。这意味着使用两种差分方法，得到的运动场差异较小。这是因为，5点差分法是结合前后5个像素点的灰度进行差分计算，其结果有平滑作用。但因为我们在运动估计前，首先对图像进行了平滑预处理，因此，采用5点差分与采用中心差分方法，其结果差异不大是合理的。本文选择了相对简单的中心差分法求解。

IBP方法作为一种迭代处理方法，Irani等^[13]证明了它的收敛性。他们指出：算法收敛与初始高分辨率图像估计无关，且能较快速收敛（通常在5次迭代以内）。同时，背投影算子p是影响数值稳定性和收敛速度的一个主要因素。p取h的逆能获得稳定的解，但收敛较慢，如果p选择脉冲函数，则可快速收敛，但数值稳定性可能会有所降低。文中，我们设定p为h的逆，采用数据集1中的图像进行了测试，其迭代误差曲线如图 3所示。该曲线呈单调下降趋势，验证了算法的收敛性。

3 结论

本文提出了一种基于快速亚像素运动估计的超分辨率重建算法。与传统的双线性插值方法相比，本文算法可使图像清晰度显著增强。同时，本文算法要优于基于全搜索运动估计的重建算法。由于运动矢量估计速度和精度的优化，本文算法不仅使最后重建图像更清晰，而且计算时间较全搜索算法缩减了近20倍，是很有意义的改进。从视觉效果和定量评估两方面结果综合说明，本文算法能较快速重建出结构更清晰的肺部图像。