基于二维超声图像的肝脏运动跟踪方法研究

庞 杰，周著黄，吴水才

（北京工业大学环境与生命学部生物医学工程系，北京 100124）

0 引言

目前，原发性肝肿瘤是我国第4 位常见恶性肿瘤及第2 位肿瘤致死病因，严重威胁我国人民的生命和健康[1]。肝切除、肝移植和影像引导下的热消融是目前肝肿瘤最有效的治疗方法[2]。由于肝肿瘤早期发现困难，确诊时已是晚期，因此常常错过最佳治疗机会，再加上没有足够的肝源等原因，只有小部分患者可以进行肝移植和肝切除治疗。热消融术因其侵入性小、普适性好、经济有效等原因，在肝肿瘤治疗领域逐步推广。

微波/射频热消融治疗是肝肿瘤的主要治疗方法之一。虽然已有不少研究人员开发了引导手术、立体定向导航、机器人辅助穿刺[3-5]等热消融辅助技术和方法，但其中大部分的模型都是通过术前CT 图像建立，并且假设肝脏处于静止状态。虽然在术中患者处于麻醉状态，但人体呼吸运动引起的肿瘤位置改变给热消融治疗带来了很大困难[6-9]。为了解决呼吸运动对肿瘤定位造成的影响，早期研究人员提出了外部标记法[10-11]来应对。该方法通过在患者身上放置传感器，并对患者进行CT 扫描，记录体外传感器的位置和对应的组织位置，通过建模从体外标志点的位置求取组织位置。但该方法有一些缺点，模型需要在术前建立，而且手术过程中该模型可能会发生偏移。而人体内部器官运动复杂，外部传感器不能完全对应，模型产生的轻微误差可能导致手术失败[12]。

许多医学成像技术被应用于临床以辅助医生治疗，包括平面X 射线、MRI 和超声。虽然X 射线可以实时使用，但由于有电离辐射，不能长时间跟踪胸腹部器官的运动。MRI 虽然可以提供高分辨力和很好的组织对比度，但其设备占地面积较大，移动不便且很难实时产生图像，不适合跟踪胸腹部器官的运动。随着计算机技术的发展，超声图像采集质量不断提高，且超声所具有的实时采集图像和无电离辐射的优势，使其成为非侵入式治疗目标跟踪的一个很好的选择。Hallack 等[13]使用logDemons 非线性配准，并通过感兴趣点周围血管大小自适应选取感兴趣区域（region of interest，ROI）对感兴趣点进行跟踪，其平均跟踪误差为1.21 mm。缺点是感兴趣点在相邻两帧中位移大小和ROI 边长大小近似时容易失败，且其算法运算速度较慢，平均帧率为12 帧/s。Bharadwaj等[14]对孪生神经网络进行适当改进，在超声图像中跟踪目标，其平均跟踪误差为1.60 mm，在感兴趣点接近超声图像边界、被遮挡以及对比度不高的情况下容易跟踪失败、产生较大误差。该研究通过添加模板匹配模块和线性卡尔曼滤波模块进行调整，使得误差更低。以上方法表明，需要更高的跟踪精度以及可以实时运行的跟踪算法来对肝脏运动进行跟踪，从而减少消融手术过程中呼吸运动带来的误差影响。

为此，本文提出一种中值流[15]、模板匹配、特征点匹配融合的跟踪方法，该方法通过超声图像跟踪肝脏运动并进行补偿，减少呼吸运动带来的影响。跟踪结果表明该方法较为准确，且运行速度快，能够满足实时要求，可以辅助手术导航系统校正呼吸运动带来的影响，提高手术治疗效果。

1 方法

1.1 实验数据与平台

本研究在Visual Studio 2017 上进行，采用CLUST（Challenge on Liver Ultrasound Tracking）2015 2D[16-19]数据集进行测试和验证。该数据集目前包含24 个训练集和39 个测试集。其中有4 种类型超声图像序列，不同类型的超声图像对比度、图像质量、感兴趣点位置分布不同，每个超声图像类型至少包含4位受试者，每位受试者有1～5 个感兴趣点，感兴趣点通常为血管中心。将本方法的跟踪结果提交至CLUST 2015 官方，根据官方返回跟踪误差大小判断本方法的跟踪精度，并进行误差分析与可行性分析。

1.2 算法流程

在CLUST 2015 竞赛中，跟踪不同类型的超声任务中，参数是固定的或者自动调节的，本文设定所有参数都是固定的。由于CLUST 2015 所给出的感兴趣点的坐标为亚像素精度，所以首先选用中值流算法，因为其可以跟踪亚像素精度的点。中值流作为短时跟踪算法，在长时间跟踪过程中会发生漂移并且积累误差，所以需要进行校正。在跟踪过程中，如果中值流算法成功跟踪到目标，那就在目标附近进行模板匹配，通过模板匹配来校正中值流算法跟踪过程的漂移；如果中值流算法未能跟踪到目标，则需使用ORB[20]（oriented fast and rotated brief）特征点进行匹配，判断ROI 的位移，重新选定ROI 继续进行跟踪，如果连续跟踪失败，则说明中值流算法不适用该情况下的跟踪，需使用KCF[21]（kernelized correlation filters）算法代替中值流算法进行跟踪。跟踪算法流程图如图1 所示。

图1 跟踪算法流程图

1.2.1 中值流算法

中值流算法是基于传统光流法的改进[22-23]。光流法是利用图像序列中像素的变化寻找前一帧图像和当前帧图像间的对应关系，进而得到两帧图像间物体运动状态的一种方法。光流法具有2 个必要的前提假设：

（1）同一个物体在图像中对应的像素亮度不变。由于光流法是根据像素亮度寻找两帧图像中目标的运动关系，因此如果像素亮度发生改变，那么将无法在两帧图像中实现同一个像素或物体的匹配。

（2）要求两帧图像必须具有较小的运动。光流法只在原图像点附近搜索对应的像素点，因此两帧图像像素位置不能有较大的移动。

假设在时间t，点（x，y）处的像素灰度值为I（x，y，t），在时间t+Δt，位置（x，y）移动到位置（x+Δx，y+Δy），此时的像素灰度值为I（x+Δx，y+Δy，t+Δt）。根据灰度一致假设，即像素的灰度在瞬时运动中应该保持一致，因此有

将公式（1）右侧项泰勒一阶展开，则有：

忽略高阶无穷小量，两边同时除以dt，化简为

稀疏光流法（Lucas-Kanade，L-K）是假设领域内其他像素的位移矢量和中心像素相同，然后通过加权最小二乘法对领域内的所有像素点进行计算。而中值流算法是基于稀疏光流法来进行改进的，稀疏光流法只计算一幅图像中ROI 部分特征点的光流强度。

中值流算法是通过度量前向、后向之间的误差来完成自我跟踪精度的提升，即在两帧相邻图像it、it+1中，先在图像it中选取Harris 角点Hi，通过稀疏光流算法计算出在图像it+1对应的角点集Hi+1，再对Hi+1反向使用稀疏光流算法计算出对应的角点集Hi'，如果角点集Hi和Hi'中对应角点位移差距过大则滤去，如图2 所示。左图中选取2 个Harris 角点（点1、点2）进行前向-反向光流计算。点1 在左右2 张图中均可成功跟踪；点2 在右图中被遮挡导致跟踪错误，反向跟踪并没有返回到点2，而是反向跟踪到了另一个点3，点3 与点2 间距离过大，说明跟踪过程发生了漂移。因此点2 在中值流算法跟踪过程中被滤去，从而解决遮挡、漂移的问题。

图2 前向-反向追踪图

1.2.2 模板匹配

模板匹配是一种最原始、最基本的模式识别方法，主要研究某一特定对象物的图案位于图像的什么地方，进而识别对象物，这就是一个匹配问题，是图像处理中最基本、最常用的匹配方法。模板就是一幅已知的小图像，而模板匹配就是在一幅大图像中搜寻目标，已知该图中有要找的目标，且该目标与模板有相同的尺寸、方向和图像元素，通过一定的算法可以在图中找到目标，确定其坐标位置。

模板匹配的度量方法是归一化相关系数匹配算法。其计算公式为

式中，w、h 分别代表模板的长和宽；I（x，y）代表输入大图像中（x，y）处的像素值大小；T（x'，y'）代表模板小图像（x'，y'）处的像素值；R（x，y）代表输入的大图像中（x，y）处的相似性度量结果。R（x，y）的值为[-1，1]，-1 代表负相关，0 代表不相关，1 代表正相关。对于本实验而言，计算结果越接近1 代表与原始模板越匹配，效果越好。

对于在超声边界附近的感兴趣点，判断模板匹配的搜索范围是否全部在超声图像内部，如果不在则不能使用模板匹配校正，因为要防止超声图像外部的信息对跟踪过程造成干扰。除此之外，在跟踪过程中，模板不能更新，应一直使用初始模板，防止跟踪过程中发生漂移。

1.2.3 ORB 特征点的检测和匹配

ORB 特征点是因已有特征点检测算法[如SIFT（scale invariant feature transform）[24]]，不能满足需求而开发出来的。SIFT[24]特征点在1999 年首次提出，并在2004 年得到进一步完善。SIFT 特征点在光照、噪声、缩放和旋转等干扰下仍具有良好的稳定性，在对象识别、图像拼接等领域应用出色。但是由于其运算速度过慢，无法应用到实时系统中，最主要的是它受到专利保护，学者不能随意使用。后续学者对SIFT算法进行了改进和优化，出现了如SURF[25]（speeded up robust features）、ORB 等更优的算法，其中ORB 特征点具有和SIFT 相似的匹配性能，受图像噪声影响较小，并且能够实时应用。ORB 特征点由FAST[26]角点和BRIEF[27]描述子组成，首先通过FAST 角点确定图像中与周围像素存在明显差异的像素点作为角点，之后计算每个角点的BRIEF 描述子，从而确定唯一的ORB 特征点。FAST 角点的核心思想是如果在某个灰度值较小的区域存在一个灰度值明显变大的像素点，那么该像素点在这个区域中具有明显的特征，可以作为特征点。FAST 角点计算过程如下：

第一步：选择某个像素点作为中心点p，其像素值为Ip。

第二步：设置判定FAST 角点的像素阈值，例如Tp=0.2×Ip。

第三步：比较中心点p 的像素值与半径为3 的圆周上所有像素的像素值，如果存在N 个像素的像素值大于Ip+Tp或者小于Ip-Tp，那么中心点p 为FAST角点。

第四步：遍历所有像素点，重复上述步骤，计算图中的FAST 角点。

FAST 角点的选择如图3 所示。

图3 FAST 角点（绿色）与周围像素的灰度值比对

BRIEF 描述子用于描述特征点周围像素灰度值的变化趋势，如果2 个图中具有相同的描述子，那么认为2 个特征点是同一个特征点。ORB 特征点匹配是通过计算2 个描述子之间的汉明距离来判断。图4 为CLUST 2015 数据集中相邻两帧ORB 特征点的匹配结果。

图4 ORB 特征点匹配示意图

1.2.4 KCF 算法

KCF 跟踪算法是一种核相关滤波算法，其核心思想就是扩充负样本数量以增强跟踪器的性能，而扩充负样本的方法就是采用循环矩阵的方法进行构造，只有基样本为正样本，其余均为负样本。

该算法采用岭回归的方法训练跟踪器。训练的目标是找到一个函数f（z）=wTz 最小化样本xi和回归目标yi的平方误差，其中λ 为正则化参数，w 的计算公式如下：

式中，X 为循环矩阵；I 为单位矩阵；y 为回归矩阵。

循环矩阵X 如公式（6）所示：

然后对循环矩阵X 进行离散傅里叶变换对角化，可变为

1.2.5 算法融合

对于感兴趣点的跟踪是基于其周围区域的跟踪，在初始帧以感兴趣点的中点选取一个40×40 像素大小的正方形，将该区域作为参考模板和ROI，对该区域使用中值流算法跟踪。如果跟踪成功，在后续图像所跟踪到的ROI 周围±2 像素进行模板匹配，如果两者跟踪到的区域中点的欧氏距离大于3 像素，则以中值流算法作为结果，否则以模板匹配作为结果，且以模板匹配结果重新初始化中值流算法ROI；如果中值流算法跟踪失败，则用ORB 特征点匹配的方法校正，如果ORB 特征点匹配计算出对应特征点的位移大于7 则去除，通过求取ORB 特征点平均位移计算出ROI 的位移大小，重新对ROI 区域初始化。在使用ORB 特征点校正中值流算法跟踪失败时，如果ROI 周围没有超过10 个ORB 特征点，那么只能使用模板匹配来跟踪。如果连续超过30 帧图像都不能使用中值流算法成功跟踪以及ORB 特征点匹配矫正失败，这种情况代表该任务无法使用中值流算法来完成，那么则使用KCF 跟踪。同样地，使用模板匹配进行长时目标跟踪误差校正。由于呼吸运动具有周期性，可以通过历史运动轨迹校正跟踪错误。记录每连续200 帧图像的跟踪结果，对于后续所跟踪到的感兴趣点判断其是否在运动范围±3 像素，如果没有超出范围，则认为正确，否则在运动范围内进行模板匹配，重新初始化ROI。

2 实验结果

对测试集中某一感兴趣点进行测试，分别使用没有校正的中值流算法和校正后的算法进行比对，其结果如图5 所示。从图5 可以看出，校正后的误差明显减小，证明校正使算法跟踪精确度明显改善。

图5 算法校正前后对比

通过CLUST 2015 竞赛官方的评判，本方法的跟踪误差为（1.29±1.43）mm（95%置信区间：3.59 mm），详见表1。

表1 CLUST 2015 各系列数据跟踪误差统计 单位：mm

由于所有任务的ROI 大小都是40×40 像素，所以跟踪速率相对稳定，平均帧率大于22 帧/s，基本满足实时要求。在某些情况下，帧率可能会突然降低，可能是由于感兴趣点突然出现大位移或者需要额外计算，例如重置位置、对ORB 特征点进行匹配等。

在所有跟踪任务中存在部分任务跟踪误差较大的情况，具体见表2。

表2 跟踪任务误差较大的情况统计 单位：mm

造成这些任务误差较大的原因可能有很多，具体如下：MED-07-3_1 误差大的可能原因是该点代表的血管过大，大于中值流算法40×40 像素的ROI，导致中值流算法不能很好地跟踪该区域，以及某些情况下会出现大位移也同样导致中值流算法产生误差；MED-11_2 误差大的原因是感兴趣点十分靠近超声图像边界，导致在跟踪过程中中值流算法和模板匹配应用效果较差；CIL-03_2 可能是由于部分遮挡导致目标被遮挡时跟踪误差较大；ETH-13-1_1可能是由于图像对比度较低以及感兴趣点所代表血管过小导致的；MED-06-1_4 也可能同样是因为感兴趣点所代表的血管过小导致；MED-07-4_1 与MED-07-3_1 导致误差过大的原因也应该是相同的，都是感兴趣点所代表血管过大导致的；MED-11_1 和MED-12_2 可能因为感兴趣点周围强度受到干扰，导致跟踪效果较差。根据所有跟踪结果的平均误差可知，图像对比度越高，处于超声图像中间、无大位移和遮挡的情况下跟踪精度越高。

3 讨论

使用手术导航系统进行肝肿瘤射频消融时，呼吸运动会引起肝脏的运动和偏移，进而导致误差。为了解决这一难题，本文方法主要基于中值流算法、模板匹配、ORB 特征点匹配和KCF 组合对二维超声图像中的感兴趣点进行跟踪。CLUST 2015 数据集由24个训练数据集和39 个测试数据集组成，分别具有53 个和85 个带标注的目标。本方法平均跟踪误差为（1.29±1.43）mm（95%置信区间：3.59 mm），具有一定的鲁棒性，快速准确，且不需要对每位患者进行训练。由于算法的运行时间主要由ROI 大小确定，且ROI 大小为固定值（40×40 像素），所以本方法的运行时间很稳定，平均帧率大于22 帧/s。Bharadwaj 等[14]使用孪生神经网络和模板匹配在超声图像中跟踪目标，添加模板匹配模块和线性卡尔曼滤波模块校正跟踪过程中的漂移，其平均跟踪误差为1.60 mm，缺点是在感兴趣点接近超声图像边界、遮挡以及对比度不高的情况下容易跟踪失败、产生较大误差。Hallack 等[13]使用logDemons 非线性配准以及使用SIFT 特征点匹配校正跟踪过程中的漂移和误差，其平均跟踪误差为1.21 mm，缺点是SIFT 特征点计算过程很慢，更适合离线操作，算法运算速度较慢，平均帧率为12 帧/s，且感兴趣点在相邻两帧中位移大小与其所选中ROI 大小近似时容易失败。本文使用中值流算法以及模板匹配和ORB 特征点匹配进行超声肝脏的跟踪，跟踪过程中的漂移和误差可以较好地跟踪肝脏运动，且ORB 特征点匹配具有很好的实时性，所以时间复杂度低、计算速度快，可以满足实时需求。

对于后续的研究，一是研究一种可以自动选择ROI 大小的跟踪算法，使得ROI 内所携带的血管信息和血管周围的信息比例保持在一个适当的范围，从而解决血管过大或过小导致跟踪误差变大的问题；二是针对遮挡、感兴趣点在超声边界以及对比度较差的情况，使用supports[28]方法来解决，该方法可以通过关联点计算目标点的位置；三是对于时间延迟问题[29-31]的补偿，由于数据传输、计算以及消融针进入人体内均需要一定的时间，如果不考虑这段进针时间或不进行时间补偿，消融针到达肿瘤上的目标点时会导致消融针进针位置错误。因此，需要对该时间延迟进行补偿，可以使用支持向量回归、卡尔曼滤波等方法来解决。

4 结论

本文设计了一种基于中值流的方法在超声图像中跟踪肝脏运动，包括跟踪、漂移校正、跟踪失败的ROI 重置。通过提交本文所提出的跟踪方法的运行结果至CLUST 2015，验证了本方法的跟踪准确性。本方法跟踪精度较高、运行速度较快，可以满足热消融手术中对肿瘤跟踪的临床需求。