王 丽,孙 宇,王瑞强,齐晓志
(1.河北工程大学机械与装备工程学院,河北 邯郸 056038;2.中国科学院深圳先进技术研究院,广东 深圳 518055;3.哈尔滨工业大学(深圳)机电工程与自动化学院,广东 深圳 518055)
腰椎管狭窄症(lumbarspinalstenosis,LSS)是骨科常见疾病,治疗LSS有切除减压与有限减压两类术式,切除减压会破坏腰椎结构影响脊柱的稳定性引起术后并发症,椎板有限减压能保留较多脊柱骨及韧带结构减少术后并发症。
目前机器人辅助外科手术广泛应用临床中,如全椎板切除、人工椎间盘置换术等[1-2],首个用于骨科手术的机器人是Curexo Technology公司的Robodoc系统,该系统根据术前扫描的CT影像定制磨削方案,并基于末端力传感器实现安全控制[3]。目前有代表性的是以色列Mazor公司推出的SpineAssist系统[4],该系统由6自由度的小型并联操作器、T型支撑架以及导航系统构成,并行操纵器在半主动模式下操作进行定位术中的手术工具。中科院深圳先进技术研究院开发的脊柱手术机器人系统由图像导航系统、控制系统和5自由度机械臂组成[5-8],并选用猪脊椎骨为样本进行磨削实验磨削误差为0.15mm,满足脊柱手术临床需求[5]。
针对椎板减压术设计的手术机器人系统,第一节介绍机器人系统结构;第二节对机器人正逆运动学分析并求解机器人的工作空间;第三节介绍基于医学影像构建轨迹规划的实际临床环境,采用Isomap与遗传算法实现机器人磨削轨迹规划;并采用A*算法实现椎板换位过程中的轨迹规划,避免机器人在三维空间运动时与非目标区域发生碰撞。
手术机器人系统由手术机器人控制系统和图像导航系统组成,图像引导手术导航是进行精准手术的关键技术,其关键在于基于3D/2D配准技术及PSU注册块实现术前CT图像与术中XRay影像、术中X-Ray影像与导航系统间的配准,这将保证医生在影像中的规划数据应用于机器人的轨迹控制,不仅可以提高手术操作的安全性,且降低医生受辐射的程度。机器人磨削环境,如图1所示。术前将机器人固定在可调浮架上,浮架前端连接有骨钳夹,通过骨钳夹与棘突的固定来保证浮架跟随脊柱在呼吸过程中浮动,降低呼吸运动导致的机器人磨削误差。
图1 椎板减压机器人工作环境Fig.1 Working Environment of Lamina Decompression Surgical Robot
五自由度的椎板减压手术机器人包括机械臂和进给装置两部分,如图2所示。由水平面内的旋转关节(1st)和两个垂直平面内的旋转关节(2nd,3rd)构成位置调节部分,姿态调节部分(4th、5th)主要是为了调整骨钻在两个方向上的角度,使骨钻能够对空间复杂区域磨削,所设计椎板减压手术机器人不仅结构简单,且适用于大多数脊柱外科手术。
图2 机器人结构示意图Fig.2 The CAD Model of the Universal Robotic Holder
手术机器人由五自由度机械臂和进给装置组成,本节基于D-H参数法建立机器人运动学分析模型,如图3所示。机器人各关节变量的D-H参数表,如表1所示。
图3 手术机器人运动学分析模型Fig.3 The Kinematics Model of the Robot
表1 机器人D-H参数列表Tab.1 The D-H Parameters of the Robot
机器人逆运动学是已知末端连杆的位姿求关节变θ1~θ5,将机器人的运动方程写为:
同理可得:
椎板减压手术机器人工作空间需满足手术区域位姿需求,如图4(a)所示。红色方块代表手术区域,文章根据机器人运动学获得机器人工作空间,如图4所示。各关节角运动范围及杆长参数不小于上文所提要求,如图显示工作空间约为半径为400mm的半圆球完全覆盖手术区域满足要求。
图4 手术区域和机器人工作空间Fig.4 Surgical Area and Robot Workspace
椎板减压术中对机器人任务规划的决策主体为医生,文章所提规划技术主要是给医生提供合适的交互方式选取三维影像椎板上下表面区域为机器人磨削提供数据支持。通过对医学影像处理提取部分椎体用于三维重建,如图5所示。实操中医生可用鼠标调整重建图像中长方体包围盒的尺寸与位姿使其合理地包围手术区域,并使包围盒的上下表面与椎板上下表面相对应,如图5(a)所示。在椎板表面搜索时对包围盒的上下表面划分矩形网格,网格交点设定成照射方向垂直于包围盒表面的光源,计算光线与椎板表面的交点得到两组离散点集来描述手术区域,如图5(b)所示。并建立虚拟约束曲面对机器人运动空间约束。
图5 图像处理Fig.5 Image Processing
采用水平集方法分割脊柱CT图像,医生根据需要在影像中选择目标点约束全局搜索范围,再采用水平集方法搜索椎体边界。因椎体轮廓在相邻切片间具有相似性,可将当前影像边界作为相邻影像的初始边界继续搜索以减少迭代次数,直至获得用于三维重建的椎体数据。
根据文献[9]提出的两相位水平集法,构建水平集函数的泛函φ避免因医学图像的不均匀引起分割算法搜索重叠区域降低目标区域识别度。
式中:εΩ—定义在图像域Ω上的能量函数数据项,M1与M2—Ω被分割出的两个不相交域的隶属函数,定义M1(φ)=H(φ)及M2(φ)=1-H(φ),其中H为Heaviside函数,ci为局部K均值聚类函数,其核函数为截断高斯函数。
式(5)为文献[9]中变分水平集能量函数Γ,式右边第二项用于计算零水平集轮廓长度作为曲线演化的平滑约束项;第三项用于控制水平集梯度稳定在1附近作为平滑性保持的惩罚项,其中p为潜在函数,定义为p(s)=(1/2)(s-1)2,水平集分割图,如图6所示。将分割出的空间数据进行中值滤波处理,并结合移动立方体法进行三维重建,同时算出重建后长方体包围盒的八个顶点坐标向量P0=(p1,p2,…,p8)。医生通过鼠标交互操作使包围盒产生平移、旋转与尺寸变换,分别为Um、Rm及Sm,顶点坐标向量可以表达为Pm=Rm(P0+Sm)+Um。
图6 水平集分割图Fig.6 Level Set Partitioning Graph
由医生通过调整上节生成的包围盒形态与位姿参数主导磨削区域规划,仅将待磨削区域包含进包围盒内,在保证有足够数据描述椎板形态前提下对包围盒表面数据降采样,以每个采样点作为光源进行光线投射获得用于描述椎板上下表面形态的点云集合Qu与Qd,使用配准算法实现规划数据向机器人坐标系映射。
式中:Nu—面向椎板上表面Mu⊃Qu的包围盒表面Bu上离散矩形网格的最大交点个数;Su—Bu在影像空间中的长宽比;eu—垂直于Bu且指向包围盒的内部的单位法向量;Nw、Nl—Bu离散后宽度方向的网格个数与长度方向的网格个数。
在三维重建时为减小计算量仅对能反映外部表面特征的数据重建,故实际体数据模型边缘厚度仅为1个像素单位。文章采用BSP-Tree的方式加速对Qu的搜索防止搜索步长过大或过小影响实际应用。
在对机器人磨削路径规划时,要保证机器人运动平稳性和轨迹曲线连续性,由于椎板磨削区域较小,要求机器人有较高的操作精度和人机协同能力,由于机器人磨削位置与深度因人而异,需要结合术中影像信息自主生成合理轨迹,为避免机器人末端位移的不确定性,其轨迹应严格限制在笛卡尔空间中。
在图像规划中获得的采样点分布于一个潜在的二维流形中,为降低轨迹计算的复杂度,采用Isomap法对环绕空间内的坐标数据SM∈R3降维,以获取能描述其内在非线性结构特征的内蕴坐标数据SN∈R2。通过降维获得椎板表面采样数据在二维空间的分布情况,同时尽可能保留成对距离信息。文章采用遗传算法求解磨削轨迹,使其路径能经过且仅经过一次所有的采样点且满足路径最短。若机器人实际运行的轨迹中出现图7(a)中较大折角的情况,会导致换向加速度过大造成系统不稳定,采用三次样条曲线对轨迹进行插值,以保证轨迹光滑性。将降维数据点一对一映射至原空间后,经过空间中三次样条插值的轨迹能较好地反映椎板表面形态,路径简单、完整且换向过程中不会产生冲击,如图7(b)所示。
图7 轨迹曲线Fig.7 Trajectory Locus
在完成磨削轨迹规划的同时采用A*算法实现椎板换位过程中的轨迹规划,避免机器人在三维空间中运动时与非目标区域发生碰撞。A*算法在人工智能中是一种典型的启发式搜索算法具有较高的目的性和效率,三维空间中A*算法的估计函数,如式(8)所示。其中,f(i)是目标点到节点i的估计值;g(i)是初始节点n到节点的实际代价;h(i)是节点n到目标节点最佳路径的估计代价。
利用Matlab对CT扫描图像进行处理,利用A*算法对椎板换位的路径规划进行验证,算法验证结果,一个区域为提取的路径通道,另一个线条为规划的路径,A、B分别为路径起始点和终止点,如图8所示。
图8 算法实施结果Fig.8 Algorithm Implementation Results
针对椎板减压术研制了一款手术机器人,该机器人由五自由度机械臂和进给机构组成。对机器人运动学分析并验证机器人工作空间满足手术要求;其次以磨削轨迹规划为目的对医学影像进行分割与重建,并采用Isomap与遗传算法对机器人的磨削轨迹进行规划;同时采用A*算法实现椎板换位过程中的轨迹规划,避免机器人在三维空间中运动时与非目标区域发生碰撞。