1 引言（Introduction）

世界卫生组织（WHO）发布的2020年全球癌症统计报告^[1]显示，在全部癌症种类中肺癌的全球死亡人数位列第一，对人类健康构成了严重威胁。在我国，每年肺癌发病约78.1万例，死亡约62.6万例，发病率、死亡率均位列全球第一^[2-3]。肺穿刺术最常见的并发症是气胸、出血，其中气胸发生率2.4%~60.0%（平均20%），出血发生率5.0% ~16.9%^[4]。若穿刺针与胸膜切面不垂直、穿刺时多次经胸膜穿刺、穿刺路径跨肺间裂或肺大疱、穿刺距离过长，容易导致肺和胸膜损伤，增加气胸发生率；若穿刺时损伤血管、穿刺次数过多，容易导致血管破裂引发出血^[5-7]。CT影像引导下的经皮肺穿刺活检作为胸部疾病诊断和鉴别的常规方法，在实际临床应用中仍然存在多方面的问题^[4]。由于无法实时采集肺穿刺活检针与病灶的位置关系，患者常常会经历反复穿刺和CT扫描确认才能成功获取活体组织，这增大了术后并发症的发生概率以及患者受到的放射性辐射剂量^[8-9]。肺穿刺手术的操作很大程度上依赖于经验^[10-11]，受直接影响的包括：手术操作时间、对患者身体的损害程度以及手术成功率。尤其是对于经验不足的医生来说，很可能会因为路径选择不当或未精准到达靶点而导致对患者的反复穿刺，加大患者术后并发症的风险，甚至可能因为没有取到适当的肿块组织导致取样失败。此外，多次的CT扫描确认也使患者遭受反复的电离辐射。

合理的肺穿刺手术路径规划能提高医生手术效率，减小患者受辐射的剂量，降低手术并发症风险，使肺穿刺手术更快速、更精确、更安全。目前穿刺手术路径规划方案主要是针对于腹部^[12-17]。Seitel^[12]利用器官硬性约束、3个软性约束对入刺路径进行筛选，并引入帕累托最优方法计算最优路径，分析了10例失败（伴有并发症）的临床案例；Schuman^[13-15]团队补充了部分约束条件，并提出了以肿瘤为原点向患者皮肤做圆柱形投影，利用投影过程中的水平、垂直偏转角进行皮肤点离散采样，需预留2 mm以上的边界误差；刘少丽等^[16-17]采用球形坐标，以靶心为中心分析可行穿刺区域，但该算法偏向理论，没有经过临床数据验证。

由于肺穿刺手术有其不同于其他器官穿刺手术的特有临床准则，有必要根据这些准则进行肺穿刺手术路径规划的研究。基于RRT（快速扩展随机树）的3维闭环针转向算法^[18]、卡尔曼滤波算法^[19]、结合了贪婪启发策略和可达引导策略的改进RRT算法^[20]、粒子群算法^[21]等一系列肺穿刺手术路径规划算法都是基于柔性针的，对医生操作的要求高且受实际力学影响难以预测路径，不常用于微创介入治疗中。硬性针穿刺具有相对确定性强、控制简单等优点，在临床中使用更为广泛。东北大学团队针对硬性针肺穿刺进行了一系列研究，牛珂楠^[22]对骨骼、肺、气管、肺部血管进行分割并重建，帮助医生观察病情并设计手术路径；陈宇^[23]提出了建立风险模型的思路，将血管与气管作为不可穿刺的重要器官；庄锦峰^[24]提出3个硬性条件与5个软性条件，细化了穿刺模型；鲍楠^[25]基于6组数据提出一种半自动路径规划方案，其中叶间裂、肿瘤分割任务需要医生手动进行，并提供多条筛选后的路径由医生分析选择。

本文利用图像分割、3维重建技术、手术导航规划技术突破现有肺穿刺手术在2维层面上展开路径设计的局限，设计了一套肺穿刺手术路径规划系统。收集来自国家三级甲等综合医院四川省自贡市第一人民医院的25套带有肺部肿块的胸部CT平扫数据，共计4871张术前胸部CT影像数据，分割出目标肿块以及患者的骨骼、心脏、血管、气管、叶间裂等重要解剖结构，并利用VTK（Visualization Toolkit）可视化工具库的面绘制算法将分割处理后的胸部CT图像展示在图像空间中。将肺穿刺手术中路径规划的风险因素抽象化为3个影响手术成败的约束性条件和6个衡量手术规划路径优劣的目标性条件，根据约束性条件筛除不符合临床要求的穿刺点，将剩余穿刺点聚类为“簇”，再利用目标性条件对“簇”的中心点进行帕累托优化和风险加权，得出最优穿刺路径。

本文工作的创新点如下：相较于其他针穿刺路径规划^{[17, 24-33]}，增大了样本数据量，采用25套带有肺部肿块的胸部CT平扫数据，避免数据量小导致的过拟合现象；使用传统学习与深度学习算法，自动对器官进行分割，减小规划前期工作量；提出自适应AGNES算法，将邻近的剩余可穿刺皮肤点凝聚成“簇”，进一步缩小可穿刺路径数；结合帕累托前沿与风险权重解决多目标最优问题，减少人为干扰。此外，为进一步提高计算效率，改进皮肤入刺角物理模型，将皮肤平面拟合的范围邻域改为$ k $近邻点，再对皮肤邻近点进行超定方程求解，减小平面拟合过程的计算量；改进路径到重要器官距离的物理模型，由膨胀算法改为根据路径方向向量与路径投影计算路径到重要器官距离，降低不同体素间距参数的CT图像距离计算的误差。

2 肺穿刺路径规划系统流程（System process of the lung biopsy path planning）

肺部穿刺手术通过胸部CT扫描检测肿块的大小和具体位置，使用针刺穿透法取肺部的局部组织进行活检，以确认体内的肿瘤是否为恶性^[34-35]。术前手术规划直接影响手术成功率，现阶段的肺部穿刺路径都是医生通过观察患者2维断层图像，根据经验手动进行规划，更注重路径长度的相关评价，很难同时考虑到2维图像上非直观的临床需求，如路径延长线与重要器官的距离、胸膜穿透角度等。

术前规划阶段有2个关键任务：明确胸腔解剖结构和规划穿刺路径。肺穿刺手术路径需要避开肺外关键结构（如肋骨、脊椎、心脏、主动脉等）和肺内组织（如肺内血管、气管、叶间裂等）^[36]，精准定位肿瘤位置，选择最短穿刺路径。因此，需要根据CT影像，对胸部皮肤区域、胸部骨骼、心脏、主动脉、肺部实质及叶间裂、肺部血管、气管、肺部肿瘤进行精准的分割，结合VTK面绘制算法立体显示技术实现胸部解剖结构分割结果的3维可视化，并根据手术约束、目标条件实现肺穿刺手术的最优路径规划，为医生提供重复性和一致性的可视化辅助，指导临床医师操作。肺穿刺路径规划系统流程如图 1所示。

3 胸部图像多器官分割流程（Multi-organ segmentation process of chest images）

在肺穿刺过程中，穿刺针对人体的伤害是不可避免的。如果规划路径穿过了重要器官或组织，可能会对人体造成极大的伤害，导致气胸、出血和咯血等并发症^[37]。为了降低手术并发症风险，需要尽量避免经过重要器官和骨骼。正确定位胸部重要器官组织是肺穿刺路径规划的基础，直接影响肺穿刺效果。需要避开的重要器官包括：肋骨、肩胛骨等骨组织，心脏及血管、气管、叶间裂。此外，作为肺穿刺路径的起点和终点，皮肤区域和肺部肿块区域也需要单独进行分割。综合考虑重要组织器官结构关系，为减小后续图像处理的计算量，将各器官分割流程设计为如图 2所示的过程，各器官分割结果如图 3所示。

胸部区域分割。通过k均值聚类^[38]将肺部周边骨骼、脂肪区域聚集为一类，肺部区域与周围背景聚集为另一类。利用形态学操作去除图像中的微小粒状噪声以及背景中的衣物，提取最大连通域，得到肺部周边骨骼、脂肪区域，对获得的最大连通域进行孔洞填充，去除中间肺部区域，实现胸腔区域分割，检测二值图像的边缘轮廓^[39]作为所需皮肤区域。

骨骼/心脏分割。根据CT图像的特点，骨骼的HU值（Hounsfield unit）在300以上，但心脏和主动脉的HU值与一些骨骼的HU值比较接近，且与周围组织HU值的差别明显，基于胸腔区域分割结果采用固定阈值^[40]将骨骼、心脏及主动脉分割出来。

气管分割^[41]。指定一个气管种子像素作为生长的起点，使用区域生长法将种子点周围邻域的像素点和种子点进行对比，将具有相似性质的点合并起来继续向外生长，分割气管组织。

肺叶分割。将胸部区域分割结果作为图像与处理结果输入深度学习U-Net网络，利用U-Net（LTRCLobes_R231）模型^[42]，融合包括肿瘤区域的左右肺分割结果与不包括肿瘤区域的肺叶分割结果，得到包含肿瘤区域的肺叶分割结果。

肺内血管分割^[43]。错综复杂的血管网络一般在图像中会呈现出线性的树状结构，黑塞矩阵的特征值能够很好地描述这种管状信息，使用基于海森矩阵的3维管状识别滤波器对管状结构进行边缘检测增强滤波器处理，再和肺叶分割结果进行“与”运算，将增强后的血管归一化并进行阈值分割，提取出肺血管组织。

肿块分割。将肺叶分割结果作为图像预处理结果输入深度学习3D U-Net网络^[44]，用3维卷积操作替换U-Net网络中的2维卷积，对肺部肿块进行分割。

4 肺穿刺术最优路径规划算法（Optimal path planning algorithm for lung biopsy）

经皮穿刺活检在肺部肿块的诊断、分期和治疗中起着十分重要的作用，选择手术路径也应该遵循一定的临床准则。通常情况下，需要先对病人进行胸部CT扫描，以确定肺部肿块的大小、位置、与周围重要解剖结构的关系，制定肺穿刺方案及路径^[45]。

图 4为肺穿刺手术环境的平面示意图。穿刺路径在避开重要脏器和肋骨、肩胛骨等骨性结构的前提下，还应避开肺大泡、大血管、气管和叶间裂，使病灶与胸膜穿刺点的距离尽可能短，同时尽可能减小通过正常肺组织的长度^[4]。此外，在实际手术操作过程中，最便捷、安全的是垂直穿刺。因此，在肺穿刺过程中需要重点考虑3个因素：穿刺长度、重要器官避障、入刺角度。

基于肺穿刺术需要遵循的临床准则，建立肺穿刺术风险模型，最优路径规划流程如图 5所示。穿刺针可到达深度、重要器官避障、入刺角度3个约束性条件是手术必须遵循的临床准则，胸壁厚度、肺内长度、路径与重要器官距离、路径延长线与重要器官距离、皮肤入刺角、胸膜入刺角6个目标性条件用于衡量手术路径优劣。通过3个约束性条件，可以将不符合临床应用的路径删除；但是剩余路径的手术风险仍然有高低之分，再在剩余的路径中利用6个目标性条件对其风险进行定量化分析，筛选出肺穿刺的最优路径。

4.1 约束性条件 4.1.1 穿刺深度

肺穿刺活检手术中穿刺深度受穿刺针长度限制，在进行路径规划时需要考虑穿刺点到肿瘤区域的长度应小于穿刺针最大可达深度，如图 6所示，阴影部分为不可穿刺区域。在实际临床中使用的肺穿刺活检针为ACHIEVE^®同轴活检针，产品型号为CA1815，由带有针座的穿刺针与同轴导引针共同构成。穿刺针的最大总长度为15 cm，通过按钮控制弹簧相继释放针芯和切割套管抽取活体组织，但考虑临床操作难度，实际中穿刺路径最长为14 cm，因此筛除穿刺深度大于14 cm的路径。

设皮肤入刺点$ P(x_{P}, y_{P}, z_{P}) $与肿瘤中心$ O(x_{0}, y_{0}, $ $ z_{0}) $之间的距离$ d $为穿刺深度，其公式为

$ \begin{align} d = \sqrt{ (x_{0} - x_{P})^{2} + (y_{0} - y_{P})^{2} + (z_{0} - z_{P})^{2}} \end{align} $

(1)

当$ P $点和$ O $点之间的距离$ d $小于或等于穿刺针可达深度14 cm时，满足穿刺深度临床需求，即穿刺针从该点入刺可达到肿瘤区域，否则不满足。

4.1.2 重要器官避障

穿刺最优路径考虑的因素较多，其中最关键的是减小相邻组织的损伤。如图 7所示，阴影部分为不可穿刺区域，若穿刺路径经过肺叶间裂、气管，则极易引发气胸；若穿刺路径损伤血管，则会导致血管破裂，严重的可能会引起出血、咯血等并发症；骨骼组织由于具有不可穿透性，一定程度上阻碍了穿刺操作，也需要进行避障。

假设皮肤点的合集为$ S $，皮肤入刺点$ P(x_{P}, y_{P}, $ $ z_{P}) $与肿瘤中心$ O(x_{0}, y_{0}, z_{0}) $的连线即穿刺路径为$ PO $，其方向向量$ \mathit{\boldsymbol{J}}= (m, n, p) $：

$ \begin{align} m & = \frac{x_{0} - x_{P}}{\sqrt{ (x_{0} - x_{P})^{2} + (y_{0} - y_{P})^{2} + (z_{0} - z_{P})^{2}}} \\ n & = \frac{y_{0} - y_{P}}{\sqrt{ (x_{0} - x_{P})^{2} + (y_{0} - y_{P})^{2} + (z_{0} - z_{P})^{2}}}\\ p & = \frac{z_{0} - z_{P}}{\sqrt{ (x_{0} - x_{P})^{2} + (y_{0} - y_{P})^{2} + (z_{0} - z_{P})^{2}}} \end{align} $

(2)

从$ P $点沿方向向量$ \mathit{\boldsymbol{J}} $搜索下一体素是否为重要器官，如果是重要器官，则返回碰撞，不满足临床需求；若直至搜寻到肿瘤中心$ O $，仍未发生碰撞，则证明该路径不会经过重要器官，满足重要器官避障临床需求。

4.1.3 皮肤入刺角

在进行路径规划时，穿刺手术的实际操作难度十分重要，穿刺针与皮肤的夹角是影响穿刺可实施性的一个难点，其物理模型示意图如图 8所示。因此，在选择穿刺路径时，需要计算穿刺路径与皮肤的夹角。当穿刺针垂直于入刺点处皮肤平面时，为最佳穿刺角90$ ^{\circ} $，理论上最便于操作。当小于20$ ^{\circ} $时，针容易打滑，导致穿刺不准确，所以在规划路径时应确保穿刺角度大于20$ ^{\circ} $。

首先，根据CT图像的体素间距参数约为0.75 mm$ \times $0.75 mm$ \times $1.50 mm，获取每个皮肤体素的8$ \times $8$ \times $4邻域，即空间体积为6.00 mm$ \times $6.00 mm $ \times $6.00 mm邻域皮肤点，再从邻域皮肤点中找到距离入刺点$ P $最近的$ k $个数据点作为该点的$ k $近邻点，直接求解超定方程对入刺点$ P $及其$ k $近邻进行平面拟合。假设拟合后的平面方程为$ z = ax + by + c $，由离散点拟合平面，实际上就是求解超定方程：

$ \begin{align} \begin{bmatrix} x_{1} & y_{1} & 1 \\ x_{2} & y_{2} & 1 \\ x_{3} & y_{3} & 1 \\ \vdots & \vdots & \vdots \\ x_{k+1} & y_{k+1} & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} z_{1} \\ z_{2} \\ z_{3} \\ \vdots \\ z_{k+1} \end{bmatrix} \end{align} $

(3)

上述方程也可以表示为$ \mathit{\boldsymbol{AX}} = \mathit{\boldsymbol{\gamma}} $。由于$ \mathit{\boldsymbol{A}} $是一个$ (k+1) \times 3 $的矩阵，首先在等号两边分别乘以$ \mathit{\boldsymbol{A}} $的转置矩阵$ \mathit{\boldsymbol{A}}^{\rm T} $，使系数矩阵变为$ 3 \times 3 $的方阵，然后乘以系数矩阵的逆矩阵解得：

$ \begin{align} {\mathit{\boldsymbol{X}}} = (\mathit{\boldsymbol{A}}\mathit{\boldsymbol{A}}^{\rm T}) ^{- 1}\mathit{\boldsymbol{A}}^{\rm T}\mathit{\boldsymbol{\gamma}} \end{align} $

(4)

该平面法向量为$ \mathit{\boldsymbol{I}}(a, b, - 1) $，路径$ PO $的方向向量为$ \mathit{\boldsymbol{J}}(m, n, p) $，可以求得经过皮肤入刺点$ P $的穿刺路径与该点处皮肤平面法向量的夹角$ \alpha $为

$ \begin{align} \alpha = \arccos \frac{ma + nb + pc}{\sqrt{m^{2} + n^{2} + p^{2}} \times \sqrt{a^{2} + b^{2} + (- 1)^{2}}} \end{align} $

(5)

当夹角$ \alpha > 90^{\circ} $时，穿刺路径$ PO $与皮肤之间的夹角$ \beta = \alpha - 90^{\circ} $；当夹角$ \alpha < 90^{\circ} $时，穿刺路径$ PO $与皮肤之间的夹角$ \beta = 90^{\circ} - \alpha $。

若穿刺路径$ PO $与皮肤间的夹角$ \beta $不小于20$ ^{\circ} $，则满足皮肤入刺角临床需求；否则，不满足。

4.2 目标性条件 4.2.1 胸壁厚度风险子目标函数

在肺穿刺手术中，路径所经过的胸壁厚度（CWT）越厚，对患者的损伤越大^[46-47]。因此在设计手术路径时，需要尽量减小穿刺针在穿刺过程中经过的胸壁厚度。设路径$ PO _{i} $对应的皮肤入刺点为$ P(x_{Pi}, y_{Pi}, z_{Pi}) $，则皮肤－胸膜长度为

$ \begin{align} d^{\rm CWT}_{i} = \sqrt{(x_{0} - x_{Pi})^{2} + (y_{0} - y_{Pi})^{2} + (z_{0} - z_{Pi})^{2}} \end{align} $

(6)

为了更为准确地对指标进行评判，需要对该指标进行归一化，归一化的公式如下：

$ \begin{align} L^{\rm{CWT}}(i) = \frac{d^{\rm CWT}_{i} - d^{\rm CWT}_{\min}}{d^{\rm CWT}_{\max} - d^{\rm CWT}_{\min}} \end{align} $

(7)

其中$ L^{\rm{CWT}}(i) $为穿刺路径所经过胸壁厚度的归一化结果，其值越低，经过的胸壁厚度越薄，对应的手术路径风险越低。$ d^{\rm SPL}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径所经过胸壁厚度的最大值，$ d^{\rm SPL}_{\min} $是对应的最小值。$ d^{\rm SPL}_{i} $是当前穿刺路径所经过的胸膜厚度。

4.2.2 肺内长度风险子目标函数

肺穿刺手术中需要尽量减小穿刺针所经过的肺内长度$ d^{\rm PLL} $（puncture length in lung），以减少手术并发症^[48-49]。设路径$ PO _{i} $与肺实质的交点为$ Q(x_{Qi}, y_{Qi}, z_{Qi}) $，则穿刺路径的肺内长度为

$ \begin{align} d^{\rm PLL}_{i} = \sqrt{(x_{0} - x_{Qi})^{2} + (y_{0} - y_{Qi})^{2} + (z_{0} - z_{Qi})^{2}} \end{align} $

(8)

其归一化公式如下：

$ \begin{align} L^{\rm PLL}(i) = \dfrac{d^{\rm PLL}_{i} - d^{\rm PLL}_{\min}}{d^{\rm PLL}_{\max} - d^{\rm PLL}_{\min}} \end{align} $

(9)

其中$ L^{\rm PLL}(i) $为肺内距离的归一化结果，该子目标函数值越低，肺内距离越短，对应的手术路径风险越低。$ d^{\rm PLL}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径的肺内最大长度，$ d^{\rm PLL}_{\min} $是对应的最小长度。$ d^{\rm PLL}_{i} $是当前路径的肺内长度。

4.2.3 重要器官距离风险子目标函数

在肺穿刺手术中，由于呼吸会带动器官移动，在设计穿刺路径时需要尽量远离重要器官，降低穿刺到重要器官的可能性^{[11, 50-51]}。路径到重要器官的距离$ D^{\rm DVO} $（distance to vital organ）的物理模型如图 9所示，计算路径到重要器官的距离可简化为计算重要器官上的点$ M(x_{M}, y_{M}, z_{M} $) 到路径所在直线$ PO $的距离。

首先，排除掉与路径夹角大于90$ ^{\circ} $的重要器官点，即重要器官点$ M $在直线$ OP $上的投影应在线段$ OP $上。根据向量的平行四边形法则，可以将问题转化为求向量$ \mathit{\boldsymbol{e}} $的模。$ \overrightarrow{OM} $为向量$ \mathit{\boldsymbol{a}} $，$ \overrightarrow{OP} $为向量$ \mathit{\boldsymbol{b}} $，向量$ \mathit{\boldsymbol{c}} $为向量$ \mathit{\boldsymbol{a}} $在向量$ \mathit{\boldsymbol{b}} $上的投影，向量$ \mathit{\boldsymbol{e}} $为$ M $点到$ OP $的垂线，根据向量的内积公式$ \mathit{\boldsymbol{a}} \cdot \mathit{\boldsymbol{b}} = | \mathit{\boldsymbol{a}} || \mathit{\boldsymbol{b}} |\cos\theta $可求得向量$ \mathit{\boldsymbol{c}} $的模：

$ \begin{align} | \mathit{\boldsymbol{c}} | = | \mathit{\boldsymbol{a}} |\cos\theta = \frac{| \mathit{\boldsymbol{a}} || \mathit{\boldsymbol{b}} |\cos\theta}{| \mathit{\boldsymbol{b}} |} = \frac{\mathit{\boldsymbol{a}} \cdot \mathit{\boldsymbol{b}}}{| \mathit{\boldsymbol{b}} |} \end{align} $

(10)

因为$ \mathit{\boldsymbol{c}} $与$ \mathit{\boldsymbol{b}} $的方向相同，所以$ \mathit{\boldsymbol{c}} $的单位向量与$ \mathit{\boldsymbol{b}} $的单位向量相同，向量$ \mathit{\boldsymbol{b}} $的单位向量为$ \mathit{\boldsymbol{J}}(m, n, p) $，整理可计算得到向量$ \mathit{\boldsymbol{c}} $：

$ \begin{align} \mathit{\boldsymbol{c}} = \frac{\mathit{\boldsymbol{a}} \cdot \mathit{\boldsymbol{b}}}{| \mathit{\boldsymbol{b}} |} \mathit{\boldsymbol{J}} \end{align} $

(11)

向量$ \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{a}} - \mathit{\boldsymbol{c}} $，再取$ \mathit{\boldsymbol{e}} $的模即可得到点到直线的距离。获取其最小值即为路径到重要器官的距离$ {D^{\rm DVO}} $，对其进行归一化得：

$ \begin{align} D^{\rm DVO}(i) = 1 - \frac{D^{\rm DVO}_{i} - D^{\rm DVO}_{\min}}{D^{\rm DVO}_{\max} - D^{\rm DVO}_{\min}} \end{align} $

(12)

其中$ D^{\rm DVO}(i) $为路径与重要器官距离的归一化结果，该子目标函数值越小，路径距离重要器官越远，对应手术路径风险越低。$ D^{\rm DVO}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径中路径到重要器官的最大距离，$ D^{\rm DVO}_{\min} $是对应的最小距离。$ D^{\rm DVO}_{i} $是当前路径到重要器官的距离。

4.2.4 路径延长线与重要器官距离风险子目标函数

路径规划时，需要在路径的延长线方向上预留一定距离（一般为2 cm），减小手术风险^{[11, 51]}。从$ O $点沿其方向向量$ {{\mathit{\boldsymbol{J}}}}(m, n, p) $搜索下一体素是否为重要器官，如果是重要器官，则返回发生碰撞时$ O $点与该体素点之间的距离作为路径延长线与重要器官的距离$ D^{\rm EDVO} $（extended distance to vital organ）；若直至搜寻到图像边缘，仍未发生碰撞，则证明该路径延长线不会经过重要器官，返回$ O $点与该图像边缘点的距离。对其进行归一化得：

$ \begin{align} D^{\rm EDVO}(i) = 1 - \frac{ D^{\rm EDVO}_{i} - D^{\rm EDVO}_{\min}}{ D^{\rm EDVO}_{\max} - D^{\rm EDVO}_{\min}} \end{align} $

(13)

其中$ D^{\rm EDVO}(i) $为路径延长线与重要器官之间距离的归一化结果，该子目标函数值越小，路径延长线到重要器官的距离越远，对应手术路径风险越低。$ D^{\rm EDVO}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径中路径延长线到重要器官的最大距离，$ D^{\rm EDVO}_{\min} $是对应的最小距离。$ D^{\rm EDVO}_{i} $是当前路径的延长线到重要器官的距离。

4.2.5 皮肤入刺角风险子目标函数

肺穿刺手术中最理想的是垂直于皮肤入刺，此时穿刺针最容易掌控，同时也更方便取材^[48]。设路径$ PO _{i} $对应的皮肤入刺角为$ \theta^{\rm PAS} $（puncture angle of skin），其归一化结果如下：

$ \begin{align} \theta^{\rm PAS}(i) = 1 - \frac{\theta^{\rm PAS}_{i} - \theta^{\rm PAS}_{\min}}{\theta^{\rm PAS}_{\max} - \theta^{\rm PAS}_{\min}} \end{align} $

(14)

其中$ \theta^{\rm PAS}(i) $为皮肤入刺角的归一化结果，该子目标函数值越小，皮肤入刺角越接近90$ ^{\circ} $，对应手术路径风险越低。$ \theta^{\rm PAS}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径的最大皮肤入刺角，$ \theta^{\rm PAS}_{\min} $是对应的最小皮肤入刺角。$ \theta^{\rm PAS}_{i} $是当前路径的皮肤入刺角。

4.2.6 胸膜入刺角风险子目标函数

穿刺针与胸膜切面不垂直时，容易引发气胸，因此需要尽量使穿刺路径与胸膜切面垂直^{[11, 51]}。设路径$ PO _{i} $对应的胸膜入刺角为$ \theta^{\rm PAP} $（puncture angle of pleura），其归一化结果如下：

$ \begin{align} {\theta^{\rm PAP}}(i) = 1 - \frac{\theta^{\rm PAP}_{i} - \theta^{\rm PAP}_{\min}}{\theta^{\rm PAP}_{\max} - \theta^{\rm PAP}_{\min}} \end{align} $

(15)

其中$ {\theta^{\rm PAP}}(i) $为胸膜入刺角的归一化结果，该子目标函数值越低，胸膜入刺角越接近90$ ^{\circ} $，对应手术路径风险越低。$ \theta^{\rm PAP}_{\max} $是所有经过约束性条件筛选的剩余穿刺路径的最大胸膜入刺角，$ \theta^{\rm PAP}_{\min} $是对应的最小胸膜入刺角。$ \theta^{\rm PAP}_{i} $是当前路径的胸膜入刺角。

4.3 最优路径筛选

经过3种约束性条件筛选后，剩余路径数量看似比较多，但实际上路径对应的入刺点主要集中在几个区域中。对本文所采用数据集而言，凝聚层次聚类算法与分散性聚类算法（例如k均值法）、密度聚类算法（例如DBSCAN算法，density-based spatial clustering of applications with noise）能够得到高质量的聚类结果，但$ m $个数据需要存储$ m^{2} $邻近度，空间复杂度为$ O(m^{2}) $，对计算量和存储空间需求大。为进一步缩小其范围，提出一种自适应凝聚层次聚类算法，其伪代码如算法1所示。

算法1：自适应凝聚层次聚类算法

输入: 皮肤入刺点样本集$\{\mathit{\boldsymbol{x}}_{1} , \mathit{\boldsymbol{x}}_{2} , \cdots, \mathit{\boldsymbol{x}}_{m}\}$、聚类簇距离度量函数${d}$ 输出: 簇中心点$\{\mathit{\boldsymbol{y}}_{1} , \mathit{\boldsymbol{y}}_{2} , \cdots, \mathit{\boldsymbol{y}}_{k}\}$ 1: for $j = 1, 2, \cdots, m$ do 2:     计算目标Calinski-Harabasz值$s(j)$ 3:     取最大的$s(j)$对应的$j$值为最佳聚类数$k$ 4: end for 5: for $ j = 1, 2, \cdots, m $ do 6:     $C_{j} = \{ \mathit{\boldsymbol{x}}_{j}\}$ 7: end for 8: for $i = 1, 2, \cdots, m $ do 9:     for $j = 1, 2, \cdots, m $ do 10:         $ {M} ( i, j ) = {d}( C_{i} , C_{j})$ 11:         $ {M} (j, i)= {M}( i, j )$ 12:     end for 13: end for 14: 设置当前聚类簇个数：$q = m$ 15: while $q > k$ do 16:     找出距离最近的两个聚类簇$C_{i^{*}}$和$C_{j^{*}}$ 17:     合并$C_{i^{*}}$和$C_{j^{*}}$：$C_{i^{*}}=C_{i^{*}}\cup C_{j^{*}}$ 18:     for $j =j^{*}+ 1, j^{*}+ 2, \cdots, q$ do 19:         将聚类簇$C_{j}$重新编号为$C_{j - 1}$ 20:     end for 21:         删除距离矩阵$M$的第$j^{*}$行和第$j^{*}$列 22:     for $j = 1, 2, \cdots, q - 1$ do 23:         $ {M}(i^{*}, j)= {d}(C_{i^{*}} , C_{j})$ 24:         $ {M}(j, i^{*})= {M}(i^{*}, j)$ 25:     end for 26:     $q = q - 1$ 27: end while 28: for $j = 1, 2, \cdots, k$ do 29:     计算簇$C_{j}$内各点与其他点的距离总和$D$ 30:     取最小距离$D_{\min}$对应的点为簇中心点${y}_{j}$ 31: end for

首先计算剩余路径点的Calinski-Harabasz分数，其计算公式如下：

$ \begin{align} s(k) = \frac{{\rm tr\; } ({{\mathit{\boldsymbol{B}}}}_{k})}{{\rm tr\; } ({{\mathit{\boldsymbol{W}}}}_{k})}\frac{m - k}{k - 1} \end{align} $

(16)

其中$ m $为样本总量，$ k $为分类的数量（簇数），$ {{\mathit{\boldsymbol{B}}}}_{k} $是簇间协方差矩阵，$ {{\mathit{\boldsymbol{W}}}}_{k} $是簇内协方差矩阵，$ {\rm tr} $是矩阵的迹（矩阵主对角线上各个元素的总和）。也就是说，类别内部数据的协方差越小越好，类别之间的协方差越大越好；Calinski-Harabasz分数越高，聚类的结果就越好。获取最优聚类结果，将属于一个区域的剩余点归为一“簇”，区域中每个点$ ( x, y, z ) $与其他$ n - 1 $个点$ (x_{i}, y_{i}, z_{i}) $距离的总和可表示为

$ \begin{align} D = \sum\limits_{i = 1}^{n - 1}\sqrt{(x - x_{i})^{2} + (y - y_{i})^{2} + (z - z_{i})^{2}} \end{align} $

(17)

其中距离最小$ D_{\min} $所对应的点作为“簇”中剩余点的最中心穿刺点，代表该“簇”。

一般情况下，与临床需求相关的6个子目标风险函数无法同时达到最优，任一个目标函数值的提高都可能以牺牲其他目标函数值为代价。为了尽量让这些目标在区域内达到最佳状态，需要对“簇”聚类结果利用帕累托前沿进行多目标优化求解。

2维帕累托最优点示意图如图 10所示，在变量空间中，找不到其他解能够优于前沿点（任意子目标函数评分都优于前沿点），前沿点的合集为该2维变量空间中的帕累托前沿，即帕累托最优解集。6个风险子目标两两获取帕累托前沿，计算每个点在帕累托前沿中出现的次数，出现次数最多的穿刺点为衡量所有子目标后的最优穿刺路径。对于帕累托前沿次数相同的穿刺点，按其重要程度赋以适当的权系数，将其乘积之和作为新的目标函数，求得总风险系数$ R(i) $：

$ \begin{align} R(i) = \sum\limits_{n = 1}^{6} \omega_{n}R_{n}(i) \end{align} $

(18)

其中$ R_{n}(i) $为路径的子目标函数值，$ \omega_{n} $为其对应的权重。

5 实验结果与分析（Experimental results and analysis）

本节采用医生盲评验证肺穿刺最优路径规划算法。选取国家三级甲等综合医院四川省自贡市第一人民医院放射科进行诊治的25例肺肿块患者的胸部CT平扫数据为研究对象，图像层厚1.5 mm，每层图像的矩阵尺度为512$ \times $512。所有患者中男11例，女14例，病灶直径12~36 mm，平均21.64 $ \pm $ 5.28 mm；年龄42~88岁，平均63.04 $ \pm $ 8.29岁。其中病灶位于左肺7例，右肺18例。纳入标准：(1) 临床资料健全者；(2) 均经病理学检查及影像学检查确定存在肺部占位性病变；(3) 经纤维支气管镜、脱落细胞学、细菌学及肺部影像学检查均未能明确病变的性质，需给予进一步肺穿刺检查治疗；(4) 血常规、凝血各指标均处于正常范围内。排除标准：(1) 合并严重心、肝、肾器质性疾病；(2) 合并严重免疫功能缺陷者；(3) 严重咳血、肺气肿及肺动脉高压者。

临床中医生对绝大部分病例的病理诊断是正确无误的，可认为是疾病判断的“金标准”^[52-54]。手术路径规划的算法评价方案大致可分为3个大类：手术实际相关参数评价^{[24, 55-59]}、仿真实验评价^[18-21]、医生注释结果对比评价^{[12-15, 25]}。其中硬性针穿刺算法^[12-25]主要由医生通过评分、排序对路径优劣进行评定。本文采用差异实验评估方法对路径规划结果进行评估，由四川省自贡市第一人民医院放射科专业医师进行评价，将盲评获取的路径排序结果作为“金标准”。此时，医生对这些路径的来源并不知情，可以避免主观因素带来的影响。对路径的评估包括2个部分：评估所有规划路径是否符合手术临床需求，对所有路径优先级进行排序。

肺穿刺路径规划系统使用的台式机内存为32.0 GB，操作系统为Windows 10，处理器为Intel Xeon CPU E5-2623 v3，主频为3.00 GHz，在PyCharm下进行开发。

5.1 肺穿刺路径规划结果与分析

利用本文提出的肺穿刺路径规划方法对25套胸部CT平扫图像进行最优路径求解，以最优路径结果在临床“金标准”中的名次作为评估指标，结果如表 1所示。

每组数据扫描获得的层数、患者的身形不同，皮肤点的数量、规划所得路径的总数也有所不同。由医生盲评结果可看出，大部分规划路径符合手术要求。大部分不符合要求的路径是因为穿刺路径未避开腋窝、前胸组织中的重要血管、神经、淋巴结群，这些血管、神经组织在胸部CT平扫图像中很难进行观察，后续可采用增强CT图像进行改进；有一条路径经过肺大泡组织，容易引发气胸，也不符合路径规划临床需求。计算得到的最优路径在“金标准”中排名均为前3，证实了本文提出的最优路径筛选方法的合理性。此外，医生盲评指出，该算法在复杂环境下的路径规划表现更好，这主要是由于复杂环境下主要考虑路径对骨骼的规避，系统可以在3维图像中解决2维图像难以合理设计的进针路线规避骨骼的问题。选择穿刺路径时需考虑穿刺针入刺角度以避免楔形进针、令路径经过的肺内长度尽量地短、使路径尽量平行于血管。以上这些问题都增加了路径规划的难度，本文将这些问题作为未来研究的内容。

以患者1数据为例，利用本文算法计算得到的最优路径解集的3维面绘制显示结果如图 11所示（俯卧位）。蓝色路径为规划出的所有路径结果，绿色路径为规划得到的最佳路径，鼠标选中后路径变为紫色，并在界面右下角弹出路径编号以及对应的6个目标性条件值。

该患者皮肤点数为256 520，即未经筛选路径总数为256 520，依次经过穿刺深度、重要器官避障、皮肤入刺角3种约束性条件筛选后剩余的穿刺点数为68 721、1 033、690，再经过自适应“簇”的聚类后，剩余的穿刺点数进一步缩小为8，其中1条路径（图中最右侧路径）未避开分布有重要血管、神经、淋巴结群的腋窝区域（胸部平扫CT中难以观察），不符合临床需求。

最优路径6个目标性条件风险子目标函数情况如图 12所示。每个柱状图由下至上风险递增，横线表示最优路径对应风险子目标函数值。

本文算法规划出的最优路径在规划出的所有路径中排序第2，“金标准”最优路径为图 11中紫色路径，两者对应6个目标性条件值，如表 2所示。

当胸膜厚度$ d^{\rm CWT} $与肺内长度$ d^{\rm PLL} $越小，路径与重要器官距离$ D^{\rm DVO} $、路径延长线与重要器官距离$ D^{\rm EDVO} $、皮肤入刺角$ \theta^{\rm PAS} $、胸膜入刺角$ \theta^{\rm PAP} $越大时，对应路径风险越小，路径越优。算法规划路径在肺内长度$ d^{\rm PLL} $、路径延长线与重要器官距离$ D^{\rm EDVO} $、皮肤入刺角$ \theta^{\rm PAS} $、胸膜入刺角$ \theta^{\rm PAP} $四个目标性条件方面优于“金标准”最优路径，但其入刺角同时包含纵向、横向2个方向，对医生操作要求更高，优先级低于横向角度入刺的“金标准”最优路径。

5.2 各约束性条件筛选强度分析

根据手术必须遵循的临床需求来筛选所有路径，筛选主要依据包括：穿刺深度不能大于穿刺针长度；不能穿刺到重要器官导致严重出血、气胸等并发症；穿刺的入刺角不能太小，避免增大手术实际操作难度。每个约束条件对路径筛选的强度不同，对25组患者数据分别计算每种约束性条件筛除的路径占路径总数的百分比，结果见图 13。

穿刺深度、重要器官避障、穿刺角度的平均筛选强度依次为71.05%、97.61%、20.21%，顺序执行这3个约束性条件时的平均筛选总强度为98.66%。其中穿刺深度、穿刺角度受患者体型、图像采集层数的影响，筛选强度的浮动较大；重要器官避障是最严格的约束性条件，筛选强度最高，这主要是由于骨骼、叶间裂限制了可行的穿刺区域，能够筛除大部分的皮肤入刺点。

以患者1的数据为例，图 14展示了各约束性条件对穿刺点的筛选效果。符合3种约束条件的手术路径各不相同，经过3个约束条件的筛选后，入刺点基本集中在一个大区域中，这主要是为减小气胸发生概率，穿刺时需要尽量避开肺叶间裂入刺。此外，穿刺路径需要避开骨骼组织（尤其是肋骨组织），入刺针可进入的区域变小，入刺点分布呈条状。

5.3 各约束性条件执行顺序分析

由于每个约束条件的定义不同，其计算方法与时间也不尽相同，3个约束条件之间的执行顺序也会影响整个约束筛选的计算效率。单独执行每种约束条件并记录其执行时间，根据其各自的执行时间确定3个约束条件的最佳执行顺序，并以此顺序执行各约束得到最终的整体执行时间。计算25组患者数据，各约束性条件的执行时间如图 15所示。

由于每组数据扫描获得的层数、患者的身形不同，皮肤点的数量与路径的总数也不同，因此进行约束性条件计算的时间相应地也会不同。穿刺深度、重要器官避障、入刺角度的计算时间最长，这是由于计算入刺角度时需要先计算入刺点的$ k $近邻后再进行平面拟合，根据拟合平面法向量$ {{\mathit{\boldsymbol{I}}}}(a, b, - 1) $与穿刺路径方向向量$ {{\mathit{\boldsymbol{J}}}}(m, n, p) $的夹角进行判定，计算量大，耗时最长。因此设计3种约束的执行顺序为穿刺深度、重要器官避障、入刺角度。单独执行这3种约束时穿刺角度的平均筛选时间依次为1.92 s、71.81 s、325.26 s，顺序执行3种约束时的平均筛选时间为38.93 s，远小于3个约束各自分别执行的时间总和。

6 结论（Conclusion）

精确的术前路径规划是手术成功的关键因素，可以提高诊断的敏感性和特异性，缩短手术时长，降低病人的不适感，减轻对病人的伤害^[27-28]。计算机辅助手术路径规划能够预测手术风险、提供术前重要参数、减少人为因素干扰、帮助外科医生更好地适应复杂的手术环境，具有较高的可靠性和准确性^[60-62]，逐渐被医学界认可并应用于腹腔^[9-27]、骨科^[63-65]、神经外科^[28-30]等领域。

本文对肺穿刺手术术前路径规划算法进行了研究，初步建立了肺穿刺手术路径规划系统，能够根据CT影像数据对目标肿块及其周围的骨骼、血管、气管等重要解剖结构进行图像分割和3维重建，显示在图像空间上，制定穿刺路径最优准则，从穿刺的进针位置、深度、角度对手术风险量化，完成合理的肺穿刺手术路径规划。实验结果表明，对25套胸部CT平扫图像利用肺穿刺规划方法进行最优路径求解，计算得到的最优路径在“金标准”中排名均为前3，满足肺穿刺手术的安全性要求。规划所得最优路径可展示为3维绘制模型，以直观、生动的方式展示在医生面前，作为术前计划应用于肺穿刺活检手术中。同时，定量化的手术路径风险指标可为医生提供参考，智能化地辅助和引导医生进行肺穿刺手术。

后续工作中，将致力于验证规划方法在实际应用中存在的不足、进一步完善肺穿刺手术路径规划风险模型；探索呼吸运动导致的器官边界误差，减小呼吸运动对手术导航精度的影响；研究肺穿刺手术导航实时定位方法，并将最优路径规划所得结果（包含入刺点、入刺角度等）通过激光标线定位等方式展示出来，实现肺穿刺手术路径规划的临床应用。此外，为降低入刺点、穿刺角度等人为误差导致手术风险的概率，将考虑进一步研究肺穿刺手术导航机器人，实现全自动手术操作。