基于二维超声-三维CT 配准的肾脏穿刺定位方法可行性研究

郑小威 雷 隆 赵保亮 李世博 武雨琦 王祥卫 何玉成 胡 颖

1(中国科学院深圳先进技术研究院 深圳 518055)

2(中国科学院大学 北京 100039)

3(深圳大学总医院 深圳 518055)

4(香港中文大学 香港 999077)

1 引 言

近年来，医学影像引导的经皮肾穿刺手术因其创伤性小、并发症少、保留肾功能等优点，已经成为了经皮肾镜取石术[1]中建立手术通道的重要手段。清晰、实时的医学成像技术可以提高术中穿刺精度，减少重复穿刺次数[2]。目前用于穿刺手术导航的医学影像技术主要有电子计算机断层扫描(CT)成像以及超声成像。由于呼吸、组织变形以及穿刺手术过程中的不确定因素，基于术前高分辨率 CT 影像的穿刺规划难以在术中准确执行[3]。而超声图像由于实时性强，检查成本低等优点成为了术中常用的影像引导手段，但较低的空间分辨率以及狭小的视野限制了其在靶点定位中的作用[4]。

为了突破传统医学影像技术的局限，研究人员提出了术中超声-术前 CT 图像配准融合的方法，结合 CT 图像高空间分辨率以及超声图像实时性强的优势，可以清晰、实时地定位靶点位置。Brooks 等[5]、Yan 等[6]和 Xu 等[7]的研究表明，多切片到体数据的多模态配准方法在计算量小的同时还能够获得接近体数据到体数据的配准融合效果。针对具体的手术场景，研究人员们提出了多种二维超声-三维 CT 的配准方法。

一般的二维超声-三维 CT 配准框架首先建立二维超声图像与 CT 体数据中切片的匹配关系，然后进行二维超声到 CT 切片的配准，最终确定二维超声与 CT 体数据的空间位置关系。在匹配关系的确定上，Cifor 等[8]通过术前手动设置的方式确定了体数据中与二维超声对应的切片，但基于这种方法，超声探头只可做绕图像平面法线方向的旋转扫描，而且只能用于患者屏息的情况。Ramalhinho 等[9]同样在术前通过手动设置的方式初始化匹配关系，并在术中基于光学定位传感器跟踪超声探头运动，进而通过物理空间和图像空间之间的坐标转换关系对匹配关系进行实时更新。但该方法没有考虑到患者呼吸运动的影响，且要求定位传感器与超声探头之间不能有遮挡。Bhattacharjee 等[10]基于血管的边缘方向特征计算二维超声与 CT 切片的对应关系，尽管该方法考虑到了患者的呼吸运动，但是只适用于血管特征丰富的器官部位。

二维超声-CT 切片配准的挑战在于，相比于 CT 图像稠密的图像信息，超声图像中所包含的信息量是稀疏的，这会导致基于灰度的配准方法难以适用。针对这种稀疏图像到稠密图像的配准问题，基于图像特征的配准是一个鲁棒性的方法。Haque 等[11]提出了一种基于血管结构特征和梯度特征的从粗到精的配准方法，当血管特征足够清晰时，配准结果令人满意。Bhattacharjee 等[10]提出了一种基于器官边缘特征以及其内部血管特征的图像配准方法。Pohlman 等[12]在肿瘤消融治疗手术的应用场景下，提出了利用消融探针的几何特征进行配准的方法。然而上述基于迭代的配准方法在实时性方面仍难以满足临床需求。目前，基于深度学习的图像配准方法逐步成为了配准领域的研究热点，并且在实时性方面也展现出了较大优势[13]。Balakrishnan 等[14]提出了一种名为 VoxelMorph 的基于无监督学习的图像配准框架。Hu 等[15]提出了一种用于多模态图像配准的弱监督卷积神经网络框架。然而，现有的深度学习配准框架对于稀疏图像特征的学习仍较为困难，在超声-CT 配准精度上仍需进一步提高。

综上所述，针对自由呼吸下经皮肾穿刺靶点的导航定位问题，本文提出了一种基于术中二维超声-术前三维 CT 实时配准的穿刺靶点定位方法。特别地，提出了一种新的基于肾脏轮廓特征的二维超声-三维 CT 由粗到精的快速配准融合方法。另外，通过超声探头标定以及空间坐标转换，将融合图像注册到手术空间，实现靶点的实时定位。

2 方法与材料

2.1 方法概述

本文提出了一种基于术中二维超声-术前三维 CT 实时配准的穿刺靶点定位方法，整体框架如图 1 所示。术前，对超声探头进行标定，采集患者腹部三维 CT 体数据并从中提取肾脏轮廓特征；术中，实时采集超声图像并从中提取肾脏轮廓特征，通过二维超声图像与 CT 体数据中切片的快速匹配方法，结合二维图像配准方法，确定超声图像与 CT 体数据的空间位置关系，进而实现二维超声-三维 CT 的融合；最后，通过空间坐标转换，将融合图像注册到手术空间，实现靶点的实时定位。

2.2 肾脏轮廓特征提取

肾脏的像素灰度以及轮廓边缘梯度是提取肾脏轮廓特征的关键信息。手术前采用基于肾脏像素灰度信息的方法提取 CT 切片中的肾脏轮廓特征。首先，基于 CT 图像窗口技术增强 CT 切片肾脏部分的对比度；其次，提取 CT 切片中的感兴趣区域(Region of Interest，ROI)，避免无关区域对目标的干扰；然后，采用高斯滤波去除图像噪声；最后，对 CT 切片进行二值化处理，得到 CT 肾脏轮廓特征掩码。

图 1 基于二维超声-三维 CT 配准的肾脏穿刺定位方法整体框架Fig. 1 Overall framework of renal puncture positioning method based on 2D US-3D CT registration

由于噪点、阴影和信号衰减产生的伪影等因素，超声图像分割是一项较为困难的任务。此外，采集过程中的超声图像成像质量对图像分割也有着较大的影响。因此本文用基于肾脏轮廓的边缘梯度信息方法提取术中二维超声图像中的肾脏轮廓特征。首先，采用中值滤波去除图像噪声并提取感兴趣区域；其次，采用 Canny 算子检测轮廓边缘特征；然后，基于拉普拉斯梯度算子提取边缘特征；最后，采用形态学操作连接断开的部分，对超声图像进行二值化处理，得到超声肾脏轮廓特征掩码。

2.3 二维超声-三维 CT 配准框架

健康状态下人体肾脏器官的形态稳定，不会跟随呼吸而发生明显的膨胀、收缩等形变[16-18]，因此本文只考虑肾脏在自由呼吸运动中的刚体运动。本文提出的二维超声-三维 CT 配准框架如下：(1)进行二维超声图像与 CT 体数据中切片的快速粗匹配，筛选出若干候选 CT 切片；(2)进行二维图像配准，得到二维超声与候选 CT 切片之间的刚性变换矩阵，并以此将候选 CT 切片变换到超声图像坐标系下；(3)通过计算变换后的候选 CT 切片与二维超声之间的图像相似性，确定二维超声与 CT 体数据中切片的精确匹配关系，同时确定二维超声与 CT 体数据之间的空间位置关系。

2.3.1 二维超声-CT 切片粗匹配

一般来说，基于图像特征相似性的切片匹配方法，都必须先通过配准计算对齐两幅图像之间的线性或非线性偏差，才能有效地计算两幅图像之间的特征相似性。本文提出了一种新的基于肾脏轮廓面积特征的切片快速粗匹配方法，利用轮廓面积特征具备旋转、平移不变性的特点，通过计算二维超声-CT 切片之间轮廓面积特征的相似性，初步确定二维超声与 CT 体数据切片之间的匹配关系。

基于术前三维 CT 体数据，可以得到肾脏 CT 切片及其对应轮廓面积的先验知识，然后通过多项式拟合的方法得到关于 CT 切片序号与肾脏轮廓面积的拟合曲线，具体如图 2 所示。术中从实时采集的超声图像中提取出肾脏轮廓特征并计算其面积，即可根据术前得到的拟合曲线粗略地筛选出与二维超声图像匹配的若干候选 CT 切片。

图 2 肾脏 CT 轮廓面积特征拟合曲线Fig. 2 Fitting curve of renal CT characteristic area

2.3.2 二维超声-三维 CT 精配准

为了获得术中二维超声图像与 CT 切片之间的刚性变换关系，需要进行 2D 图像刚性配准。由于肾脏器官内部血管特征信息稀疏，基于肾脏轮廓特征计算二维超声到 CT 切片的刚性变换矩阵是一个鲁棒性的方法。给定超声图像 IUS和若干候选 CT 切片 ICT，二者之间的刚性变换 通过对公式(1)进行迭代优化获得，基于梯度下降的优化器根据目标函数 Loss 的梯度不断更新配准参数，直至目标函数收敛至最小值或达到最大迭代次数[19]。

经配准计算得到二维超声图像与候选 CT 切片之间的刚性变换矩阵后，将候选 CT 切片变换到超声图像坐标系中，采用骰子系数(Dice Coefficient)测量肾脏轮廓特征的相似性，如公式(2)所示：

候选 CT 切片中 Dice 值最高的即为与二维超声图像对应的最佳 CT 切片，由此确定了二维超声图像与 CT 体数据中切片的精确匹配关系，同时确定了超声图像与 CT 体数据之间的空间位置关系。

2.4 超声探头标定及穿刺定位方法

为了确定超声图像坐标系(后续简称为超声坐标系)与超声探头坐标系(后续简称为探头坐标系)之间的空间转换关系，首先对探头进行标定。本文提出了一种基于改进的空间单 N 线模型(图3(a))的超声探头标定方法，标定模型由直径 0.5 mm 的尼龙线以及高精度 3D 打印的外壳组成。N 线由两根平行线以及斜线组成，其中斜线位于平行线上方，仅需单个 N 线即可求取二维超声图像中 x、y 轴方向上的像素间距。

像素间距 sx、sy的计算方法如图 3(b)所示。已知标定模型上各点之间的距离，模型坐标系以 A 点为原点、AB 为 x 轴、AC 为 y 轴、AI 为 z 轴，则 EF 在模型坐标系 z 轴方向上的投影与 HD 相等；使用机械臂把持探头确保探头平面 EFG 平行于 AB 与 CD，则模型坐标系中 EG 的长度与 AB、CD 相等。在 N 线横截面的超声图像中，可得 EG 之间的像素行数 u 以及 EF 之间的像素列数 v，如图 3(d)所示。从而计算出超声图像 x、y 方向上的像素间距 sx、sy，如公式(5)所示：

图 3 空间单 N 线模型Fig. 3 Spatial single N-wire phantom

3 结 果

3.1 超声探头标定精度评估

本文采用光学定位跟踪系统(Polaris Vicra，Northern Digital Inc，Canada)作为定位传感器，包括双目摄像头以及若干个无源光学靶点。其中的光学定位探针经过针尖标定(标定误差为 0.026 mm)后可直接读取针尖在传感器坐标系中的坐标。用于标定的超声探头为超声诊断系统(Affiniti 30，Philips Inc，Holland)的凸阵探头(C6-2)。

超声探头标定数据采集流程如下：(1)超声探头由六轴机械臂(UR5)通过探头夹具把持，机械臂末端-超声探头-光学靶点之间均采用刚性连接的方式固定以保证三者之间的相对位姿固定不变，同时确保在扫描过程中，探头的 x、y、z 轴分别与模型的 y、z、x 轴平行；(2)探头沿模型 y 轴方向扫描盛满水的 N 线模型，采集 N 线的横断面超声图像，同时记录每张图像对应的探头靶点位姿，具体如图 4 所示。

每次标定实验采集 35 组数据，其中 20 组数据用于求解标定矩阵，15 组数据用于评估标定精度。评估方法为将特征点 F 在图像坐标系中的坐标值转换到模型坐标系中，计算其与模型坐标系中点 F 的真实坐标值之间的距离 D(公式(12))，同时计算距离 D 的均方根误差 RMSE，如公式(13)所示：

其中，n 为特征点数量。

本文总共进行了 5 次标定实验，计算所得的标定精度如表 1，平均标定误差为 0.998 mm。

表 1 超声探头标定点均方根误差Table 1 Root-mean-square error of ultrasonic probe calibration

图 4 标定数据采集场景Fig. 4 Calibration data collection scenario

3.2 二维超声-三维 CT 配准

3.2.1 数据集

本文采用人体腹部模型(057A，CIRS，USA)进行实验验证(图 5)。体模的三维 CT 体数据图像尺寸为 512×512×123，切片层厚为 1.25 mm，像素间距为 0.676 mm/0.676 mm。二维超声图像尺寸为 1 024×768，像素间距为 0.265 mm/0.265 mm。

图 5 人体腹部模型Fig. 5 Human abdomen model

为了评估二维超声到三维 CT 体数据的切片匹配精度，需要获取二维超声与 CT 切片的真实匹配关系。在超声图像采集实验中，首先，六轴机械臂把持经过标定的超声探头，沿探头 x 轴方向采集二维肾脏超声图像，同时记录固定于探头的光学靶点位姿。其次，使用经过标定的光学定位探针读取若干个黏贴在体模底部的直径为 1 mm 的小钢珠坐标。然后，通过空间坐标转换，计算传感器坐标系下超声图像与体模底部特征点之间在探头扫描方向上的距离，结合三维 CT 体数据的切片层厚，得到与超声图像真实对应的 CT 切片序号。最终获取了 5 组真实对应的体模肾脏二维超声图像-CT 切片图像对。另外，为了评估二维超声与 CT 切片的配准精度，实验过程请临床医生在三维 CT 体数据与二维超声图像上分别标注了 5 组对应的解剖标记点。

3.2.2 配准精度评估

表 2 二维超声-三维 CT 配准误差Table 2 2D US-3D CT registration erro

其中，n 为标记点数量。

二维超声-三维 CT 配准实验在切片匹配以及二维超声-CT 切片的配准误差如表 2 所示。所提出的方法正确地计算出了二维超声-三维 CT 的切片匹配关系，二维超声-CT 切片平均标记点配准误差为 0.709 mm，平均配准运算时间为 1.15 s。二维超声-CT 切片的融合图像以及其与 CT 体数据之间的空间位置关系分别如图 6 所示。

图 6 二维超声-三维 CT 融合图像Fig. 6 Fusion image of 2D US-3D CT

3.3 穿刺定位方法精度

为了评估本文所提定位方法的定位误差，将体模底部小钢珠作为特征点，获取特征点分别位于传感器坐标系、CT 图像坐标系中的坐标 PSensor、PCT。将特征点在 CT 图像坐标系中的坐标转换到传感器坐标系中，计算其与传感器坐标系中对应点的坐标值之间的距离，作为定位误差 PE，计算公式如(15)所示：

其中，n 为特征点个数；

最终选取了 5 个特征点，其定位误差如表 3 所示，平均定位误差为 2.265 mm。

表 3 定位误差Table 3 Positioning error

4 讨论与分析

针对自由呼吸下经皮肾穿刺靶点的导航定位问题，本文提出了一种基于术中二维超声-术前三维 CT 实时配准的穿刺靶点定位方法。其中，本文提出的二维超声-三维 CT 由粗到精的快速配准融合方法，不需要术前手动设置，即可在术中切片粗匹配阶段，基于肾脏轮廓特征高效率、大幅度地缩小匹配搜索范围。在精配准阶段，基于肾脏轮廓特征进行二维图像配准，得到二维超声与候选 CT 切片之间的刚性变换矩阵。同时，以此将候选 CT 切片变换到超声图像坐标系中，通过计算变换后的候选 CT 切片与二维超声之间的图像相似性，确定二维超声图像与 CT 体数据中切片的精确匹配关系，确保了切片匹配的容错率。此外，由于自由呼吸下肾脏轮廓面积特征以及肾脏轮廓特征基本稳定不变，因此整个配准过程将不受呼吸运动影响，理论上本文方法适用于自由呼吸运动下的配准。

在超声探头标定实验中，本文提出了一种改进的空间单 N 线标定模型。相比于多 N 线模型，所提出模型结构简单，并充分利用超声图像中的像素间距信息，减少对标定模型精度的依赖。在标定数据采集实验中，使用机械臂替代人手把持超声探头，减少了实验系统误差以及人为因素对于标定精度的影响，最终的标定精度为 0.998 mm。Lin 等[21]提出了一种基于 N 形楔子模型的标定方法获得了 1.730 mm 的标定精度；Carbajal 等[22]基于 N 线标定模型获得了 1.400 mm 的标定精度。与上述方法相比，本文的标定方法取得了更为准确的标定结果。在基于体模验证的配准方法中，Gomes-Fonseca 等[23]提出了一种基于器官表面的 CT-超声配准方法，使用猪肾模型进行实验验证，标记点配准误差为 1.260 mm。Zhang 等[24]提出了一种用于腹腔镜肾部分切除的视频增强现实导航的无标记自动形变注册框架，使用人体肾脏模型，目标配准误差为 1.280 mm。与上述方法相比，本文提出的二维超声-三维 CT 配准方法，标记点配准误差为 0.709 mm，取得了较好的结果。在最终的定位精度评估实验中，体模肾脏中的病变靶点直径约为 7.500 mm，本文的平均定位误差为 2.265 mm，满足体模肾脏穿刺的定位需求。

本文工作尚存在以下不足：(1)基于特征的配准方法，其局限性在于配准精度受特征分割精度影响较大；(2)当前所提出的方法只适用于超声图像所在平面与 CT 切片所在平面相互平行的情况；(3)所提出的方法仅在体模上取得了较好的实验结果。未来将进一步研究任意平面下二维超声-三维 CT 的配准问题，并在临床数据上进行实验验证。

5 结 论

针对自由呼吸下经皮肾穿刺靶点的导航定位问题，本文提出了一种基于术中二维超声-术前三维 CT 实时配准的穿刺靶点定位方法，将术中二维超声与术前三维 CT 配准后的融合图像注册到手术空间，提供导航信息。经实验评估，探头标定精度为 0.998 mm，配准误差为 0.709 mm，平均配准运算时间为 1.15 s，平均定位误差为 2.265 mm。本文方法在定位精度以及运算时间方面均取得了较好的结果，基本满足经皮肾穿刺手术对于导航定位的需求，这将有效提高穿刺手术的成功率。