面向L5/S1微创椎间盘射频消融术的多约束最优穿刺路径规划算法

A multi-constraint optimal puncture path planning algorithm for percutaneous interventional radiofrequency thermal fusion of the L5/S1 segments

摘要

目的 为有效减少L5/S1微创椎间盘射频消融术中因医师操作熟练度差异导致的治疗效果波动，实现对消融针穿刺路径风险的精确量化评估。 方法 基于自主研发的深度神经网络DWT-UNet，将L5/S1节段的MR图像精准分割为7个关键结构：L5椎骨，S1椎骨，髂骨，L5/S1节段的椎间盘、神经根、硬脊膜和皮肤。基于上述分割结果，3D建模穿刺路径的规划环境。针对穿刺的临床准则提出了6个硬性约束和6个软性约束，将医师经验通过层次分析算法量化为权重添加到穿刺路径风险函数中以提升对个例的适应性。通过提出的皮肤进针点采样子算法、Kambin三角投影面积子算法，结合层次分析算法，并运用光线追踪、CPU多线程处理和GPU并行计算等多种技术，计算出一条既符合临床硬性约束条件，又能使软性约束条件综合达到最优的穿刺路径。 结果 医师团队对算法规划的21例消融针穿刺路径进行主观评估，结果显示所有路径均满足临床基本要求（及格率100%），且95.24%的路径被评为优秀或良好。在与医师规划结果的对比中，算法在D_Ilium、D_S1和Depth等3个指标上有不同程度的下降（P<0.05），在D_Dura、D_L5、D_N5和A_Kambin等4个指标上有较大提升（P<0.05）。算法规划21例的平均时间为7.97±3.73 s，而医师传统规划普遍需耗时10 min以上。 结论 本研究提出的多约束最优穿刺路径规划算法为L5/S1节段的微创椎间盘射频消融术提供了一种一站式的解决方案。经过严谨的临床评价和实验验证，该算法显示出独特的优势和广阔的临床应用潜力。

Abstract

Objective To minimize variations in treatment outcomes of L5/S1 percutaneous intervertebral radiofrequency thermocoagulation (PIRFT) arising from physician proficiency and achieve precise quantitative risk assessment of the puncture paths. Methods We used a self-developed deep neural network DWT-UNet for automatic segmentation of the magnetic resonance (MR) images of the L5/S1 segments into 7 key structures: L5, S1, Ilium, Disc, N5, Dura mater, and Skin, based on which a needle insertion path planning environment was modeled. Six hard constraints and 6 soft constraints were proposed based on clinical criteria for needle insertion, and the physician's experience was quantified into weights using the analytic hierarchy process and incorporated into the risk function for needle insertion paths to enhance individual case adaptability. By leveraging the proposed skin entry point sampling sub-algorithm and Kambin's triangle projection area sub-algorithm in conjunction with the analytic hierarchy process, and employing various technologies such as ray tracing, CPU multi-threading, and GPU parallel computing, a puncture path was calculated that not only met clinical hard constraints but also optimized the overall soft constraints. Results A surgical team conducted a subjective evaluation of the 21 needle puncture paths planned by the algorithm, and all the paths met the clinical requirements, with 95.24% of them rated excellent or good. Compared with the physician's planning results, the plans generated by the algorithm showed inferior D_Ilium, D_S1, and Depth (P<0.05) but much better D_Dura, D_L5, D_N5, and A_Kambin (P<0.05). In the 21 cases, the planning time of the algorithm averaged 7.97±3.73 s, much shorter than that by the physicians (typically beyond 10 min). Conclusion The multi-constraint optimal puncture path planning algorithm offers an efficient automated solution for PIRFT of the L5/S1 segments with great potentials for clinical application.

Graphical abstract

关键词

微创椎间盘射频消融术 / L5/S1 / Kambin三角 / 路径规划

Key words

percutaneous intervertebral radiofrequency thermocoagulation / L5/S1 / Kambin's triangle / path planning

腰椎间盘突出症（LDH）^［1］是一种高发的脊柱退行性疾病，近年来其发病年龄趋于年轻化^［2］，且以L4/L5和L5/S1椎间盘突出最为常见^［3］，占所有LDH病例的95%左右^［4］。对于早期LDH患者，经皮穿刺微创手术^{［5， 6］}，特别是微创椎间盘射频消融术（PIRFT）已被证实为替代传统椎间盘切除术的有效微创疗法^［7］。然而，未经充分训练的医师在执行PIRFT时，可能面临高手术失败率、残余椎间盘突出和神经根压迫等潜在风险^［8］。特别地，L5/S1节段的穿刺手术因其椎间孔受到髂骨和横突的阻碍，难以维持与椎间隙平行的同时避免髂骨碰撞^［9］，因而显得尤为具有挑战性。此外，由于穿刺操作缺乏直观的3D视野，可能会增加误伤神经、硬脊膜与血管等重要软组织的风险^［10］。因此，针对L5/S1节段穿刺手术的高挑战性，术前精确的穿刺路径规划对于提升手术安全性与成功率具有至关重要的临床意义^［11］。

主流术前规划软件如Mimics和3D Slicer^［12］，虽在医学图像处理上表现出色，但针对PIRFT的专门术前规划模块尚属空白。这导致在模拟PIRFT穿刺时，软件功能主要限于基础3D可视化，难以满足临床PIRFT复杂精细的规划需求，且操作繁琐。

近年来，穿刺路径规划领域的算法研究取得了重要进展，刚性针和柔性针^［13］作为关键工具，分别成为研究焦点。刚性针路径规划方面，有研究将临床准则转化为硬性约束和软性约束，并将软性约束纳入目标函数中实现路径优化，其两步优化法能在一定程度上逼近全局最优解，为临床实践提供了有价值的参考^{［14， 15］}。有学者构建了半自动路径规划系统，该系统不仅丰富了约束条件，还提出了多目标优化，并首次引入了帕累托优化概念，医师无需再手动设定权重^［16］。另有学者通过投影简化3D规划为2D，并提出了一种基于规则排列和逐步调整的启发式算法在3分钟内完成规划^［17］。柔性针路径规划则涌现出基于快速扩展随机树的三维闭环针转向算法^［18］、卡尔曼滤波算法^［19］、结合了贪婪启发策略和可达引导策略的改进快速扩展随机树算法^［20］以及粒子群算法^［21］等创新方法，旨在提升规划精度。然而，柔性针在复杂的人体环境中操作繁琐，且其路径轨迹多变且难以预测^［22］，这极大地增加了精确建模与预测的难度。综合考量操作简便性、轨迹可预测性以及其在临床中的广泛应用^［11］，本研究选定刚性针作为PIRFT的穿刺工具。

在刚性针穿刺路径规划的算法演进中，一系列共通的挑战与局限性逐渐浮现，主要有4个挑战：自动精确建模的挑战，既往研究^［14-17］依赖CT图像手动分割，效率低下且易受软组织分辨率限制；L5/S1节段特异性约束条件的缺失，该区域复杂解剖结构增加了约束条件设定的难度；最优路径的决策挑战，单目标或多目标优化均存在局限性，前者可能简化权重设置^［15］，后者则增加医师决策负担^［16］；计算效率与全局最优解的平衡挑战，穷举法计算量大^［15］，基于优化的方法易陷入局部最优^{［15， 17］}。为应对上述挑战，本研究提出了一个针对L5/S1节段的PIRFT多约束最优穿刺路径规划算法，构建了一个全自动、高效且便捷的一站式PIRFT穿刺路径规划解决方案，旨在降低医师的工作负担，提高手术成功率，并减少患者并发症风险。现报道如下。

1 材料和方法

1.1 数据获取及规划环境建模

本研究数据源自中山大学附属第五医院脊柱外科，共收集21例确诊为LDH的病例MR影像数据，图像矩阵大小为320×320，层厚为0.8 mm。本研究得到了中山大学附属第五医院医学伦理委员会的批准（审批号：ChiCTR2100052883）。所有数据已经过脱敏处理，不含患者个人隐私信息。

首先将分辨率为0.75 mm×0.75 mm×0.8 mm的MR图像重采样到各向同性的0.75 mm×0.75 mm×0.75 mm分辨率。然后将图像灰度值归一化到［0，1］范围后输入进DWT-UNet网络，从而分割出L5椎骨结构（L5），S1椎骨结构（S1），髂骨结构（Ilium），L5/S1节段的椎间盘结构（Disc）、神经根结构（N5）、硬脊膜结构（Dura mater）和皮肤结构（Skin）等7个关键组织结构并进行3D重建。为了降低计算量，需要对3D模型进行适当裁剪，最后建立PIRFT规划环境（图1）。

1.2 算法整体流程

步骤1：输入L5/S1节段的MR图像；步骤2：DWT-UNet对MRI数据分割，得到L5、S1、Ilium、L5~S1节段的Disc、Dura mater、N5和Skin，并对分割后的结构进行3D重建；步骤3：在皮肤表面进行进针点采样并与医师标定的靶点形成候选穿刺路径；步骤4：选择一个进针点所形成的穿刺路径进行临床硬性约束条件检测，剔除掉不满足条件的皮肤进针点；步骤5：如果满足步骤4的所有硬性约束条件则进入临床软性约束条件计算，否则返回步骤4；步骤6：检查是否遍历完所有皮肤进针点，否则继续步骤4；步骤7：基于本研究提出的风险目标函数对满足临床硬性约束条件的皮肤进针点进行评估；步骤8：基于最优解确定皮肤最优进针点及对应穿刺路径的相关参数（图2）。

1.3 基于3D U-Net的深度神经网络分割模型DWT-UNet

相比于Transformer^［23］，基于卷积神经网络的3D U-Net因其独特的U型网络结构和跳跃连接机制，能更有效地捕获上下文信息及细粒度特征，尤为契合处理高分辨率但标签资源稀缺的3D医学图像任务^［24］。在本研究中，鉴于其高效性和实用性，将其选定为路径规划建模中MR图像分割任务的基础框架。然而，值得注意的是，3D U-Net中的下/上采样卷积操作可能会导致关键信息丢失且难以在后续步骤中恢复。为了有效应对这一挑战，本研究引入了一个3D离散小波变换（DWT）/逆变换（IDWT）模块来缓解这一问题，模型整体框架（图3）。

在3D DWT中，通过并置操作，模型的特征通道数扩展至原始输入通道数的8倍，不仅丰富了特征的层次结构，还成功捕获了多尺度的空间频率信息。同时，图像尺寸被缩减至原始输入的1/2，既减轻了计算负担，又促使模型聚焦于更高级别的特征表征。为了平衡特征维度扩张与计算效率之间的需求，额外设计了一个卷积模块，该模块将输出特征通道数从原始输入通道数的8倍缩减至2倍，确保数据形状变换与原始3D U-Net中下采样卷积层的功能性保持一致。同理，3D IDWT亦然。

本研究综合采用了交叉熵损失函数和Dice损失函数作为联合损失函数，旨在平衡像素级分类的准确性与区域一致性。该联合损失函数定义如下：

{\mathbf{L}}_{\mathbf{t}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l}}={\mathbf{L}}_{\mathbf{D}\mathbf{i}\mathbf{c}\mathbf{e}}+{\mathbf{L}}_{\mathbf{C}\mathbf{E}}（1）

{\mathit{L}}_{\mathbf{D}\mathbf{i}\mathbf{c}\mathbf{e}}=\mathrm{1}-\frac{\mathrm{2}{\sum }_{\mathit{i}=\mathrm{1}}^{\mathit{N}} {\sum }_{\mathit{j}=\mathrm{1}}^{\mathit{M}} {\mathit{y}}_{\mathit{i}\mathit{j}}{\mathit{p}}_{\mathit{i}\mathit{j}}+\mathit{ϵ}}{{\sum }_{\mathit{i}=\mathrm{1}}^{\mathit{N}} {\sum }_{\mathit{j}=\mathrm{1}}^{\mathit{M}} {\mathit{y}}_{\mathit{i}\mathit{j}}^{\mathrm{2}}+{\sum }_{\mathit{i}=\mathrm{1}}^{\mathit{N}} {\sum }_{\mathit{j}=\mathrm{1}}^{\mathit{M}} {\mathit{p}}_{\mathit{i}\mathit{j}}^{\mathrm{2}}+\mathit{ϵ}}（2）

{\mathbf{L}}_{\mathbf{C}\mathbf{E}}=-{\sum }_{\mathit{i}=\mathrm{1}}^{\mathit{N}} {\sum }_{\mathit{j}=\mathrm{1}}^{\mathit{M}} {\mathbf{y}}_{\mathbf{i}\mathbf{j}}\mathbf{l}\mathbf{o}\mathbf{g}\left({\mathbf{p}}_{\mathbf{i}\mathbf{j}}\right)（3）

其中，L_{\mathrm{D}\mathrm{i}\mathrm{c}\mathrm{e}}表示Dice损失函数，L_{\mathrm{C}\mathrm{E}}表示交叉熵损失函数，N 表示样本数量，M 表示类别数量，y_{ij}表示第i个样本是否属于第j个类别的真实标签，p_{ij}表示模型预测第 i个样本属于第j个类别的概率，ϵ是出于数值稳定性考虑加入的一个小的常数项，取1e-12。

在模型训练阶段，本研究使用显存为12 GB的NVIDIA GTX 1080Ti GPU训练，采用SGD优化器，初始学习率设置为3e-4，通过余弦退火法逐渐降低到6e-5，批处理参数设置为1，迭代次数设置为1000。采用五折交叉验证方法客观评估模型性能并选取最优参数配置。最后，为了进一步提升分割结果的准确性和临床实用性，对最优模型输出的分割结果应用空洞填充和孤立假阳性区域剔除等形态学后处理技术。

1.4 约束条件

为了确保计算机辅助手术导航的成功与安全，既往研究将穿刺消融术约束量化为硬性约束和软性约束^［14-16］。硬性约束作为手术安全性的基石，严格界定手术条件，确保不与关键结构碰撞，构建出可行解空间。而软性约束则通过精细的代价函数优化穿刺路径，为医师提供丰富的手术路径分析参数。本研究深入融合了既往研究^［14-16］在计算机辅助穿刺路径规划领域的先驱性工作，系统性地归纳并总结了针对L5/S1节段PIRFT路径规划的关键约束条件（表1A~C）。此外，本研究深入剖析了L5/S1节段PIRFT经椎间孔入路的特点，特别是Kambin三角这一关键解剖结构的独特优势，提出了特定约束条件（表1D~F）。

穿刺针进入人体后引起的组织变形建模过于复杂，故暂未将其纳入考虑范围。违反消融针穿刺硬性约束示例（图4）。

1.5 最优穿刺路径筛选

1.5.1 穿刺路径风险总目标函数

针对消融针穿刺路径规划的软性约束提出了3个风险子目标函数，并通过自开发的层次分析算法（AHP）计算器由医师团队按其重要程度赋以权重系数。穿刺路径的风险总目标函数即为其乘积之和，表示为：

{\mathrm{R}}_{\mathrm{P}}={\mathrm{W}}_{\mathrm{P}}{\sum }_{\mathrm{i}=\mathrm{1}}^{\mathrm{n}}{\mathrm{\alpha }}_{\mathrm{i}}{\mathrm{D}}_{\mathrm{i}}+{\mathrm{W}}_{\mathrm{P}}\mathrm{\gamma }{\mathrm{S}}_{\mathrm{p}}+{\mathrm{W}}_{\mathrm{p}}\mathrm{\delta }{\mathrm{L}}_{\mathrm{p}}（4）

\Sigma {\alpha }_i+r+\delta =\mathrm{1}（5）

其中，R_P为穿刺路径的总风险值，其取值范围是［0，1］，其值越靠近0，代表该路径风险越小，越靠近1，路径风险越大。W_P是一个二值数，当穿刺路径满足所有硬性约束条件时，其值为1，否则为0。D_i是第i个关键结构距离风险子目标函数，共有5个关键结构被纳入计算，{\alpha }_i为对应关键结构距离的权重系数。S_P是Kambin三角投影面积风险子目标函数，\gamma为对应的权重系数。L_P是穿刺深度风险子目标函数，\delta为对应的权重系数。

1.5.2 关键结构距离风险子目标函数

在实际穿刺手术中，规划路径需谨慎并尽可能远离关键组织结构，以应对呼吸运动、操作误差和组织间相互作用力导致的穿刺偏移，从而提升手术容错性并降低误伤风险。三角面片是3D网格模型的基本单元，穿刺路径到关键结构3D模型距离d的物理模型（图5）。

将距离计算简化为空间中一条线段P₁P₂上的点到一个三角形面片ABC上的点的最小距离，线段和平面不相交，即求：

d=min(d_{\mathrm{1}},d_{\mathrm{2}},d_{\mathrm{3}},d_{\mathrm{4}})（6）

为了将不同量纲的数据转换到同一尺度上，减少异常值对建模的影响，需要对该指标进行归一化，归一化的公式如下：

D_i\left(j\right)=\mathrm{1}-\frac{d_j^{}-d_{min}^{}}{d_{max}^{}-d_{min}^{}}（7）

其中，D_i\left(j\right)为穿刺路径到关键结构i距离的归一化结果，本研究中共考虑Ilium、L5、S1、N5、Dura mater等5个关键结构。D_i\left(j\right)函数值越小，其对应路径到关键结构i的距离越大，对应其手术路径风险越低。d_{min}^{}是经过硬性约束条件检测通过后的所有候选路径中到关键结构的最小距离，d_{max}^{}为最大距离，d_j^{}是当前穿刺路径到关键结构的距离。

1.5.3 Kambin三角投影面积风险子目标函数

腰椎间孔的Kambin三角是一个由神经根、下位椎体上终板及邻近结构构成的安全穿刺区域。该区域血管和神经根较少，允许穿刺针灵活操作，避免关键结构损伤。由于Kambin三角的3D概念在临床应用中难以精确描述，本研究将其投影至2D空间，并定义了准确的2D边界，以便进行量化分析和计算。在避障可达空间中X、Y、Z 等3个轴最小和最大位置的极限视角下Kambin三角2D投影的示意图（图6），局部放大图中红色线勾画区域为该视角下Kambin三角的2D投影。

本研究提出的计算Kambin三角投影面积子算法的步骤如下：

步骤1：创建计算Kambin三角投影面积所需的场景。场景由L5、S1、Disc、N5和Dura mater等5个关键结构的3D模型组成。场景的背景颜色置为黑色，L5颜色置为白色， S1置为红色，Disc置为蓝色，N5置为绿色，Dura mater置为红色。

步骤2：裁剪3D模型，减少计算。在以靶点为球心，半径为30的球体范围内裁剪模型。

步骤3：调整场景相机的参数。首先，将相机的投影模式置为透视投影，以模拟出真实世界中物体在人眼中的投影效果。然后，设定相机的初始位置为原点\left[\mathrm{0,0},\mathrm{0}\right]，并初始化其朝向向量V为\left[\mathrm{0,0},\mathrm{1}\right]，即默认朝向正Z轴方向。随后，为将相机对准候选穿刺路径上的靶点，需计算目标朝向向量T，该向量由皮肤穿刺点与靶点之间的连线确定。接着，通过比较初始朝向向量V与目标朝向向量T，计算必要的旋转角度，即欧拉角，以实现相机朝向的调整。

计算向量V、T之间的夹角\alpha：

\alpha =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{V\cdot T}{‖V‖\cdot ‖V‖}\right)（8）

计算以原点为固定点，V转向T的旋转轴向量\mathrm{B}：

B=V\times T（9）

根据罗德里格旋转公式计算旋转矩阵：

R=I+\mathrm{s}\mathrm{i}\mathrm{n}(a)\cdot K+(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}(a))\cdot K^{\mathrm{2}}（10）

K=\left[\begin{array}{l}\mathrm{0}-{\stackrel{̑}{B}}_z{\stackrel{̑}{B}}_y\\ {\stackrel{̑}{B}}_z\mathrm{0}-{\stackrel{̑}{B}}_x\\ -{\stackrel{̑}{B}}_z{\stackrel{̑}{B}}_x\mathrm{0}\end{array}\right]，\stackrel{̑}{B}=\frac{\mathrm{B}}{‖\mathrm{B}‖}（11）

调整相机空间位置。将靶点G设为相机的观察焦点，相机的位置与靶点的距离d设置为150，水平视野和垂直视野角度均为30°，计算相机的旋转变换矩阵Transform。

\mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}=\left[\begin{array}{cc}R& C_w+d\cdot R[:\mathrm{3,2}]\\ \mathrm{0}& \mathrm{1}\end{array}\right]（12）

C_w=R\cdot C_c+O， C_c=R^{-T}\cdot G（13）

步骤4：给L5重新着色。L5椎骨3D模型的三角面片到相机位置的距离小于靶点与相机距离的置为黑色，否则为白色。

步骤5：获得相机图像。通过渲染管线得到相机图像。

步骤6：计算相机图像中的Kambin三角投影面积。Kambin三角的投影区域（图7A）；灰度化相机图像后提取的Kambin三角区域（图7B）。

在PIRFT中，为了更准确地将Kambin三角投影面积纳入风险函数，避免异常值和不同指标绝对值差异过大的影响，需要对Kambin三角投影面积进行归一化，归一化的公式如下：

S_j=\mathrm{1}-\frac{s_j^{}-s_{min}^{}}{s_{max}^{}-s_{min}^{}}（14）

其中S_j为Kambin三角投影的归一化结果，该子目标函数值越小，Kambin三角投影面积越大，穿刺手术工作空间越大，对应的穿刺路径风险也就越低。s_{min}^{}是经过硬性约束条件检测通过后的所有候选路径中最小的Kambin三角面积，s_{max}^{}为最大面积，s_j^{}是当前穿刺路径的Kambin三角面积。

1.5.4 穿刺深度风险子目标函数

在PIRFT中，穿刺深度必须小于所选穿刺针的最大设计可达深度。皮肤上的黑色部分为可行的皮肤穿刺区域（图8）。

在实际临床中使用的射频热凝套管针产品规格为20 G，由导电针头、套管和连接管等构成。穿刺针的长度为150 mm，限定穿刺深度不得超出此值。在PIRFT中，穿刺针进入人体的部分越多，对患者的损伤越大^［25］。规划时应优先采取最小化原则，力求穿刺路径长度最短。设某一皮肤进针点P(x_{\mathrm{1}},y_{\mathrm{1}},z_{\mathrm{1}})与靶点G(x_{\mathrm{2}},y_{\mathrm{2}},z_{\mathrm{2}})形成穿刺路径PG，则穿刺深度为：

d=\sqrt[]{(x_{\mathrm{1}}-x_{\mathrm{2}}{)}^{\mathrm{2}}+(y_{\mathrm{1}}-y_{\mathrm{2}}{)}^{\mathrm{2}}+(z_{\mathrm{1}}-z_{\mathrm{2}}{)}^{\mathrm{2}}}（15）

为了更准确地将该指标纳入风险函数，避免异常值和不同指标绝对值差异过大的影响，需要对穿刺路径长度进行归一化，归一化的公式如下：

L_j=\frac{d_j^{}-d_{min}^{}}{d_{max}^{}-d_{min}^{}}（16）

其中L_j为穿刺深度的归一化结果，该子目标函数值越小，穿刺深度越小，穿刺针进入人体的部分越少，对应的穿刺路径风险也就越低。d_{min}^{}是经过硬性约束条件检测通过后的所有候选路径中最小的穿刺深度，d_{max}^{}为最大深度，d_j^{}是当前穿刺路径的深度。

1.5.5 调整风险权重的AHP计算器

AHP是一种系统化多目标决策方法，通过层次分解、判断矩阵和数学方法（如特征值分解）计算权重，实现决策过程透明化，确保结果可信稳定。首先，建立层次结构，目标层为筛选出最优穿刺路径，准则层有候选路径与L5、S1、Ilium、N5和Dura mater等5个结构的距离，候选穿刺路径深度，候选穿刺路径对应的Kambin三角投影面积。方案层为各个候选穿刺路径。然后，医师通过问卷调查决定准则层两两之间的相对重要性，以防一次性考虑多个指标之间的关系时出现疏漏。使用1~9标度法进行打分，1表示同等重要，9表示极端重要。根据两两对比的结果，构建判断矩阵（表2）。

经过计算判断矩阵的特征向量和特征值、一致性检验和综合权重计算等步骤后便可以得到准则层每个准则的权重（图9）。

1.6 皮肤进针点采样子算法

本研究运用先进的光线追踪技术，致力于精确模拟穿刺路径，即从皮肤进针点至特定靶点的直线轨迹。在穿刺路径规划的过程中，靶点的选择是基于临床医师的手术需求精确设定的。皮肤上所有可行进针点和靶点的连线构成了消融针的避障可达作业空间。本研究仅对皮肤的3D模型腹膜侧表面进行采样，并根据临床医师的丰富经验，将采样区域限定在较为精确的小范围内^{［26， 27］}。采样频率设定为5个点/mm²，从而在降低整体采样量的同时，保证了局部区域的采样密度。该算法伪代码（表3）。最后通过自靶点向候选进针点发射射线，并检测射线是否与关键解剖结构（如椎骨、神经根等）发生碰撞和是否满足其他硬性约束条件，以判断路径的可行性。

1.7 3D符号距离场的优化计算

符号距离场的概念在于描述空间中任意一点到某对象表面的最小有符号距离。有学者将距离场定义为“场内任意一点到场内任意对象上距离最近的点的距离”^［28］。在数学中可以用函数F表示符号距离场，其表达式为：

F\left(p\right)=sign\left(p\right)\mathrm{*}\underset{x\in s}{min}‖p-x‖（17）

向函数输入定义域中的任意点p都将得到函数值域所对应的表示的该点到空间中对象的有符号距离。

本研究在每条候选路径上均匀选取1000个样本点，以其到特定解剖结构的最小距离作为路径与该结构的近似距离。此方法在确保候选路径与关键结构距离准确评估的同时，避免了全局距离场计算，实现了对Ilium、L5、S1、N5和Dura mater等结构的快速有效避让。

以抽样点与关键结构的最小距离代替路径与关键结构的距离会带来误差，但是其小到足以忽略不计，证明如下：

本研究中使用的射频热凝套管针长度为150 mm，现假设一种误差最大的极端情况：穿刺深度为针的最大长度150 mm，路径上某一点P为距离某结构最近的一点但其不在采样点中，采样点中距离点P的最近点为点Q，结构上对应的点为R，问题简化为用QR代替PQ所带来的误差（图10）：

易得：

P{Q}_{max}\le \frac{\mathrm{1}}{\mathrm{2}}\mathrm{*}\frac{\mathrm{150}}{\mathrm{1000}}（18）

假设PR为1mm，可得误差：

E_{rror}=\frac{QR-PR}{PR}\le \frac{\sqrt[]{\mathrm{1}+\mathrm{0}.{\mathrm{075}}^{\mathrm{2}}}-\mathrm{1}}{\mathrm{1}}\le \mathrm{0.0028}（19）

由该证明可知，在极端情况下，所采用的抽样方法引起的误差被控制在0.28%以内。该方法通过将直线到3D模型的最小距离问题简化为点集到点集的最小距离计算，有效降低了计算复杂度。通过将原本CPU处理的复杂算法运算转化为GPU上大量简单重复的计算，将原本分钟级的距离场运算缩短至毫秒级，从而极大提升了算法在实际应用中的性能。

1.8 评价指标

1.8.1 L5/S1多目标分割效果定量评价

为评估本研究提出的DWT-UNet脊柱多目标分割方法的性能，将神经网络的分割结果与临床医师的人工标注进行对比。其中，Dice相似性系数（DSC）作为衡量图像分割质量的关键标准用于量化自动分割结果与专家标注之间的空间一致性或重叠程度。

1.8.2 路径优劣的评价指标

采用综合主观与客观的评价策略，全面评估穿刺路径规划算法的性能。客观评估指标涵盖待评路径与L5（D_L5）、S1（D_S1）、Ilium（D_Ilium）、Dura mater（D_Dura）和N5（D_N5）等5个关键结构的距离、穿刺深度（Depth）、Kambin三角面积（A_Kambin）及规划所需时间。为确保临床评价的严谨性、准确性及公正性，本研究采用了双重独立盲审的评估机制，引入2位具备高年资脊柱外科临床经验的医师进行分级评定。评定体系细化为5个明确等级：优秀、良好、中等、及格与不及格，严格遵循核心标准（表4）。同时，必须认识到专家经验的独特价值，因此在评定体系中融入了专家自由裁量权，允许医师根据丰富的临床经验对特定情况进行灵活判断。若2位初始评审医师的评定结果存在分歧，将引入第3位高资质医师参与讨论，确保最终评定结果的公正、准确与权威性。

1.9 统计学分析

采用SPSS25.0软件对实验结果进行统计学分析。计量资料符合正态分布的以均数±标准差表示，配对样本的差异性检验采用配对t检验或两相关样本秩和检验。以P<0.05为差异有统计学意义。

2 结果

2.1 DWT-UNet语义分割实验结果

在相同的实验条件下，所有结构的DWT-UNet 模型分割结果在DSC指标上均取得了最佳效果，分别为 94.91%、 93.73%、 94.88%、 91.97%、 82.74%。相较于nnFormer^［29］、DeepLabv3^［30］和U-Net++^［31］，DWT-UNet的S1结构提升了1.67%、7.72%、1.60%，Disc结构提升了0.94%、12.59%、0.73%，L5结构提升了1.86%、12.69%、1.44%，Dura mater结构提升了1.86%、12.69%、1.44%，N5结构提升了 4.30%、63.10%、3.27%（表5）。

2.2 消融实验结果

消融实验结果显示，S1、Disc、L5、Dura mater、N5等结构的DSC指标均获得了不同程度的提升，其中Disc、Dura mater和N5获得提升较大。相较于基准模型，加入DWT后，3D U-Net^［24］分别提升了0.4%、1.54%、0.80%、1.79%、3.67%，ResUNet^［32］分别提升了0.39%、1.36%、0.82%、1.53%、3.71%，V-Net^［33］分别提升了0.23%、0.47%、0.80%、1.23%、2.05%。Disc、N5和Dura mater等非骨性结构的提升较为明显，融合了DWT模块的3D U-Net网络的表现超越了ResUNet、V-Net及其相应的DWT集成版本，性能最佳（表6）。

2.3 算法规划的最优穿刺路径结果

算法对21例脊柱MR平扫图像进行了消融针穿刺路径规划（表7）。21例病例中，男6例，女15例，年龄26~64（39.81±11.50）岁。病例的初始候选进针点数34 120~45 875（40 242.38±3577.66）个；符合硬性约束的皮肤进针点数3917~22 097（12 224.62±4357.67）个。所有病例的穿刺深度在150 mm以内，平均深度（113.92±17.44）mm。在临床评价上，21个病例（100%）的算法规划路径均达到了及格标准，其中10个病例（47.62%）被评价为“优秀”， 10个病例（47.62%）被评价为“良好”，1个病例（4.76%）被评价为“中等”。

以病例1数据为例，算法规划得到的穿刺路径在路径规划建模环境中的结果显示，红色路径为医师规划的路径，绿色路径为算法规划的最优路径，绿色路径更偏向S1关节突和横突之间的通道，距离神经根和硬脊膜更远，手术工作空间更大（图11）。

2.4 算法与医师规划结果的对比结果

实验结果显示，对于D_iliac指标，算法规划的穿刺路径平均值低于医师的规划路径，降低了约20.54%［（9.44±5.10） mm vs （11.88±6.14） mm， P<0.05］。D_S1指标和Depth指标分别降低了约53.60%（P<0.001）、4.68%（P<0.001）。

对于D_Dura指标，算法的平均值稍高于医师的平均值，提升约1.22%［（17.44±2.38） mm vs （14.00±1.90） mm， P<0.05］。D_L5、D_N5、A_Kambin等3个指标分别提升了约12.78%（P<0.001）、8.30%（P<0.001）、14.91%（P<0.05，表8）。

2.5 算法在不同分割方法上的路径规划结果

算法在自动分割和手动分割数据上，对于D_Dura、D_L5、D_N5、Depth和A_Kambin等5个指标，两种分割方法的差异无统计学意义（P>0.05）。但是对于D_iliac指标，自动分割数据比手动分割数据低约20.13%［（7.54±5.54） mm vs （9.44±5.10） mm， P<0.05］。对于D_S1指标，自动分割数据比手动分割数据高约84.48%（表9）。

2.6 算法穿刺路径规划的时间结果

算法在自动分割数据上的平均规划时间相较于手动分割数据存在约0.52 s的增幅，但两者差异无统计学意义［（8.49±3.05） s vs （7.97±3.73） s，P>0.05］。算法的规划时间远短于医师普遍所需的10 min。

3 讨论

近年来，脊柱外科手术机器人在融合多项先进技术的基础上，显著提升了手术精度与效率，有效减少了人为错误，促进了术后恢复^［34］。然而，其应用的进一步发展仍受限于自动化与智能化水平的不足，导致在如路径规划等关键手术步骤中易受人为不稳定因素的干扰，为技术的持续优化提出了挑战^［35］。尽管Mimics和3D Slicer^［12］ 等主流规划软件正积极适应多样化的手术场景需求，但在针对特定手术如L5/S1的PIRFT穿刺路径的精细化规划方面，仍面临较大挑战，且用户需承担较高的学习成本以掌握其复杂操作。鉴于此，本研究深入探索并优化了L5/S1的PIRFT穿刺路径规划，同时根据临床实际需求定制开发了相应功能。

在穿刺路径规划的环境建模问题上，传统的图像分割方法（如Mimics和3D Slicer中的阈值分割、区域生长）依赖于手动设定的阈值和启发式规则，难以自适应L5/S1节段腰椎结构的解剖复杂性（如L5和S1）和多样性（如周围软组织），从而限制了分割的准确性和鲁棒性。深度学习方法，特别是卷积神经网络，通过自动化的特征学习和自适应优化^［24］，能够深入挖掘数据中的复杂模式，有效提升了分割精度，同时大幅减少了人工干预和特征工程的复杂度。随着神经网络层数的增加，卷积操作在提取深层特征的同时不可避免地导致信息损失，特别是在3D医学图像分割中。本研究在3D U-Net^［24］架构中融入了3D DWT模块，提出了DWT-UNet网络，通过多尺度保留细节信息，有效提高了对复杂解剖结构的识别能力。模型对比实验显示，DWT-UNet在各类解剖结构分割的DSC上均优于3D U-Net^［24］、U-Net++^［31］、nnFormer^［29］、DeepLabv3^［30］、V-Net^［33］和ResUNet^［32］等模型，尤其在Disc、N5和Dura mater等组织结构的分割任务中表现更为突出。消融实验验证了DWT模块在3D U-Net、V-Net和ResUNet中的提升效果，特别是在非骨性结构分割上的精度提升较为明显。综合实验结果， 3D U-Net结合DWT模块，展示了最佳的L5/S1节段分割能力，并通过该网络一站式快速精准构建符合临床需求的PIRFT路径规划环境。

在穿刺路径规划的约束条件问题上，早期研究者将临床准则转化为硬性约束和软性约束^{［14-16， 36］}。尽管上述研究与本研究存在一定的共通之处，但并未针对L5/S1节段的PIRFT进行探讨，不可生搬硬套前人的约束条件。本研究在借鉴既往研究的基础上^［14-16］，创新性地提出了6个硬性约束和6个软性约束条件，以提升路径规划的安全性。针对L5/S1节段独特的Ilium结构和Kambin三角这一关键解剖通道，本研究通过基于Kambin三角面积的软性约束进一步扩大手术安全操作空间，并开发了相应的面积计算算法。此外，结合关键结构的3D有符号距离场，有效避开骨骼密集区域和神经密集区域，从而进一步提高穿刺的容错率。鉴于穿刺方向与椎间盘之间的切角在当前分析框架下对穿刺路径的风险影响相对有限，因此在构建风险总函数时，本研究审慎地决定不将其纳入考量范围。尽管如此，在平衡考量范围的全面性与风险模型表达的简洁性方面，仍有进一步优化的空间。

在穿刺路径的最优决策问题上，或基于穷举的暴力求解全局最优解^{［14， 15］}，但计算复杂度高；或基于优化来逼近全局最优解^{［17， 37］}，但无法避免陷入局部最优解；又或将决策权交由医师，从帕累托前沿中挑选最优路径^［16］。本研究聚焦于临床需求的明确性和风险量化的精准性，通过权重系数法将多目标优化问题简化为单目标优化，进而构建了一个总风险评估函数。为适应个体的特异性，医师可以借助AHP计算器根据临床需求快速便捷的调整软性约束条件的相对重要性，实现权重系数的个例化微调。通过对比分析算法与医生规划的结果，本研究所提出的算法在规划穿刺路径时会以牺牲与Ilium和S1两个骨性结构距离的代价以换取距离N5和Dura mater等软组织距离更远，穿刺深度更小，Kambin三角投影面积更大的优化。根据医师的专业经验，这种规划倾向使得穿刺路径更靠近横突和椎弓根从而提升手术的安全性。

在变量优化和参数减少的策略上，本研究受到既往研究^［37］启发，通过合并相似路径策略优化搜索空间，改进了文献^［38］的蒙特卡洛皮肤进针点采样法，候选进针点数减半至约40 242个，实现了三角面片级的穷举求解。此外，在面对为L5、S1、Ilium、N5及Dura mater等5个关键结构构建3D距离场时近五千万次点到模型最小距离的计算挑战时，本研究利用CPU多线程优化Kambin三角面积的计算，借助GPU并行计算精准优化后的局部距离场和候选路径的风险。实验结果显示，21个病例的平均规划时间为7.97±3.73 s。

本研究算法的优势在于：首先，这是一个针对L5/S1节段PIRFT的全自动路径规划算法，实现了从病例MR图像输入到精确路径规划结果输出的自动化流程。整个流程只需医师简单标定消融靶点，缩减了忙于各种任务的外科医师在繁琐手术计划与策略制定上的时间投入^［39］，减少了人为误差从而提升手术的安全性与效率。其次，本研究提出的DWT-UNet网络，专为L5/S1节段MR图像分割任务定制优化，经过对比实验和消融实验验证，其性能在多个维度上展现出不凡实力，极大提高了路径规划效率。再次，本研究在继承前人研究的基础上^［14-16］，对临床约束进行了全面而精细的深化与拓展，创新性地提出了包含6项硬性约束与6项软性约束在内的综合体系。这些约束不仅涵盖了广泛临床病例中的关键安全因素，而且还融入了L5/S1节段独特的Ilium和Kambin三角这一关键解剖通道。最后，本研究不仅提供了基于严谨统计分析和医师团队经验得出的各风险子函数权重分配，更是从实践角度出发开发了一个AHP计算器。其摒弃了传统上仅凭经验或直觉设定权重的做法，赋予了医师根据具体临床情境灵活调整权重的能力，确保了决策制定的科学性、客观性与个性化。

本研究在PIRFT穿刺路径规划上取得一定进展，但实际应用中仍存在优化空间。深度学习算法虽展现出强大的潜力，但其全面应用尚需技术突破，特别是在处理骨性结构与非骨性结构交界处及梯度剧变区域的精准度上。尽管深度学习算法普遍擅长骨性结构分割，但是算法在不同分割方法上的路径规划结果显示D_iliac和D_S1两个指标却具有显著性差异，下一步工作将继续研究提高深度学习算法的性能。为了更好的结合脊椎外科手术机器人，将进一步地探索脊柱配准技术^［40］，融合多模态图像手术导航系统，以精准导航优化手术流程与效果，推动PIRFT穿刺路径规划向智能化、精准化迈进。

综上所述，本研究提出针对L5/S1节段的PIRFT穿刺路径规划算法具备高精度、高安全以及高度自动化的特点。该算法可以对L5/S1节段的MR图像进行精准分割与三维重建，并能够在3D可视化模型的基础上进行PIRFT的路径规划，为医师提供切实可行的穿刺方案。

参考文献

[1] 周谋望, 岳寿伟, 何成奇, 等. “腰椎间盘突出症的康复治疗” 中国专家共识[J]. 中国康复医学杂志, 2017, 32(2): 129-35.

[2] Ucar D, Duman S, Bayram Y, et al. Extruded disc herniations are experienced earlier by inactive young people in the high-tech gaming era[J]. J Med Life, 2021, 14(3): 402-7.

[3] Deyo RA, Mirza SK. CLINICAL PRACTICE. herniated lumbar intervertebral disk[J]. N Engl J Med, 2016, 374(18): 1763-72.

[4] Martin MD, Boxell CM, Malone DG. Pathophysiology of lumbar disc degeneration: a review of the literature[J]. Neurosurg Focus, 2002, 13(2): E1.

[5] Beyaz SG, Ülgen AM, Kaya B, et al. A novel combination technique: three points of epiduroscopic laser neural decompression and percutaneous laser disc decompression with the Ho: YAG laser in an MSU classification 3AB herniated disc[J]. Pain Pract, 2020, 20(5): 501-9.

[6] Wei WB, Dang SJ, Wei L, et al. Transforaminal epidural steroid injection combined with radio frequency for the treatment of lumbar disc herniation: a 2-year follow-up[J]. BMC Musculoskelet Disord, 2021, 22(1): 347.

[7] 王锁良, 吕建国, 申晓东, 等. 射频热凝靶点消融术在颈腰椎间盘突出症中的应用[J]. 西安交通大学学报(医学版), 2008, 29(2): 236-8.

[8] Wang HW, Huang B, Li CQ, et al. Learning curve for percutaneous endoscopic lumbar discectomy depending on the surgeon’s training level of minimally invasive spine surgery[J]. Clin Neurol Neurosurg, 2013, 115(10): 1987-91.

[9] Ruetten S, Komp M, Godolias G. An extreme lateral access for the surgery of lumbar disc herniations inside the spinal canal using the full-endoscopic uniportal transforaminal approach-technique and prospective results of 463 patients[J]. Spine, 2005, 30(22): 2570-8.

[10] Ju CI. Technical considerations of the transforaminal approach for lumbar disk herniation[J]. World Neurosurg, 2021, 145: 597-611.

[11] Scorza D, El Hadji S, Cortés C, et al. Surgical planning assistance in keyhole and percutaneous surgery: a systematic review[J]. Med Image Anal, 2021, 67: 101820.

[12] Kikinis R, Pieper SD, Vosburgh KG. 3D slicer: a platform for subject-specific image analysis, visualization, and clinical support[M]// Intraoperative Imaging and Image-Guided Therapy. New York: Springer, 2014: 277-89.

[13] Huang Y, Yu LF, Zhang FF. A survey on puncture models and path planning algorithms of bevel-tipped flexible needles[J]. Heliyon, 2024, 10(3): e25002.

[14] Baegert C, Villard C, Schreck P, et al. Multi-criteria trajectory planning for hepatic radiofrequency ablation[J]. Med Image Comput Comput Assist Interv, 2007, 10(Pt 2): 676-84.

[15] Baegert C, Villard C, Schreck P, et al. Trajectory optimization for the planning of percutaneous radiofrequency ablation of hepatic tumors[J]. Comput Aided Surg, 2007, 12(2): 82-90.

[16] Seitel A, Engel M, Sommer CM, et al. Computer-assisted trajectory planning for percutaneous needle insertions[J]. Med Phys, 2011, 38(6): 3246-59.

[17] Li RK, An CY, Wang SX, et al. A heuristic method for rapid and automatic radiofrequency ablation planning of liver tumors[J]. Int J Comput Assist Radiol Surg, 2023, 18(12): 2213-21.

[18] Patil S, Burgner J, Webster RJ 3rd, et al. Needle steering in 3-D via rapid replanning[J]. IEEE Trans Robot, 2014, 30(4): 853-64.

[19] Zhao X, Guo H, Ye D, et al. Comparison of Estimation and Control Methods for Flexible Needle in 2D[C]//东北大学, IEEE新加坡工业电子分会. 第28届中国控制与决策会议论文集（下）, 2016: 6.

[20] 赵燕江, 黄 磊, 杜海艳, 等. 基于改进RRT算法的套管柔性针运动规划[J]. 仪器仪表学报, 2017, 38(3): 620-8.

[21] Cai CX, Sun CS, Han Y, et al. Clinical flexible needle puncture path planning based on particle swarm optimization[J]. Comput Methods Programs Biomed, 2020, 193: 105511.

[22] Shahriari N, Georgiadis JR, Oudkerk M, et al. Hybrid control algorithm for flexible needle steering: demonstration in phantom and human cadaver[J]. PLoS One, 2018, 13(12): e0210052.

[23] Vaswani A, Shazeer N, Parmar N, et al. Attention is all you need[EB/OL]. 2017: arXiv: 1706.03762.

[24] Çiçek Ö, Abdulkadir A, Lienkamp SS, et al. 3D U-net: learning dense volumetric segmentation from sparse annotation[C]//International Conference on Medical Image Computing and Computer-Assisted Intervention. Cham: Springer, 2016: 424-32.

[25] Ma HP, Hung YF, Tsai SH, et al. Predictions of the length of lumbar puncture needles[J]. Comput Math Methods Med, 2014, 2014: 732694.

[26] Eun SS, Lee SH, Liu WC, et al. A novel preoperative trajectory evaluation method for L5-S1 transforaminal percutaneous endoscopic lumbar discectomy[J]. Spine J, 2018, 18(7): 1286-91.

[27] Tezuka F, Sakai T, Abe M, et al. Anatomical considerations of the iliac crest on percutaneous endoscopic discectomy using a transforaminal approach[J]. Spine J, 2017, 17(12): 1875-80.

[28] Jones MW, Baerentzen JA, Sramek M. 3D distance fields: a survey of techniques and applications[J]. IEEE Trans Vis Comput Graph, 2006, 12(4): 581-99.

[29] Zhou HY, Guo JS, Zhang YH, et al. nnFormer: volumetric medical image segmentation via a 3D transformer[J]. IEEE Trans Image Process, 2023, 32: 4036-45.

[30] Chen LC, Papandreou G, Schroff F, et al. Rethinking atrous convolution for semantic image segmentation[J/OL]. arXiv E Prints, 2017: arXiv:

[31] Zhou ZW, Siddiquee MMR, Tajbakhsh N, et al. UNet: redesigning skip connections to exploit multiscale features in image segmentation[J]. IEEE Trans Med Imag, 2020, 39(6): 1856-67.

[32] Diakogiannis FI, Waldner F, Caccetta P, et al. ResUNet-a: a deep learning framework for semantic segmentation of remotely sensed data[J]. ISPRS J Photogramm Remote Sens, 2020, 162: 94-114.

[33] Milletari F, Navab N, Ahmadi SA. V-net: fully convolutional neural networks for volumetric medical image segmentation[C]//2016 Fourth International Conference on 3D Vision (3DV). Stanford, CA, USA. IEEE, 2016: 565-71.

[34] Overley SC, Cho SK, Mehta AI, et al. Navigation and robotics in spinal surgery: where are we now?[J]. Neurosurgery, 2017, 80(3S): S86-99.

[35] 李玉希, 黄浚燊, 刘 婷, 等. 脊柱手术机器人在脊柱外科手术中的应用进展[J]. 机器人外科学杂志: 中英文, 2021, 2(2): 143-50.

[36] Villard C, Baegert C, Schreck P, et al. Optimal trajectories computation within regions of interest for hepatic RFA planning[J]. Med Image Comput Comput Assist Interv, 2005, 8(Pt 2): 49-56.

[37] Liang LB, Cool D, Kakani N, et al. Multiple objective planning for thermal ablation of liver tumors[J]. Int J Comput Assist Radiol Surg, 2020, 15(11): 1775-86.

[38] 刘晓熙, 姜慧研, 骆 敏. 面向肝癌消融术的多约束最优穿刺路径规划算法[J]. 东北大学学报: 自然科学版, 2023, 44(7): 922-30.

[39] Liu YJ. Potential risk of intelligent technologies in clinical orthopedics[J]. Adv Exp Med Biol, 2018, 1093: 281-8.

[40] Zhao L, Pang SM, Chen YF, et al. SpineRegNet: spine Registration Network for volumetric MR and CT image by the joint estimation of an affine-elastic deformation field[J]. Med Image Anal, 2023, 86: 102786.

基金资助

RIGHTS & PERMISSIONS

版权所有©《南方医科大学学报》编辑部2021