进入某一 层组织后，操作者会感受到阻力的增减变化，这对操作者来说是一个提示，有经验的操作者 可据此判断针尖所处位置。本系统可仿真穿刺过程中阻力的变化，从而为学员提供逼真的 训练效果。本系统可训练硬脊膜外腔减压和推注麻醉药液进行硬脊膜外麻醉，抽取脑脊液 等手术项目。
[0041] -个成功的穿刺包括以下阶段：
[0042] ①针尚未接触皮肤并逐渐刺破皮肤，此时针尖受到来自皮肤的弹性力 stiffness；
[0043] ②针尖刺破皮肤之后，针尖透过皮肤，来自皮肤的弹性力消失；针尖继续前进而且 未接触黄韧带之前，这个阶段针体逐渐穿过脂肪和肌肉组织，此过程中由于脂肪和肌肉之 间没有质地坚韧的筋膜，针尖持续刺透软组织，针尖一直受到切割力cutting的作用；
[0044] ③针尖到达黄韧带并逐渐突破刺穿黄韧带，此过程中针尖受到来自黄韧带的弹性 力stiffness;
[0045] ④针尖刺穿黄韧带以后并继续前进，针尖受到切割力cutting;
[0046] ⑤针尖到达硬脑膜并逐渐突破刺穿硬脑膜，此过程中针尖受到来自硬脑膜的弹性 力stiffness;
[0047] ⑥针尖刺穿硬脑膜以后到达脑脊液，此时针尖没有阻力。
[0048] 本发明是一种支持力觉交互的腰椎穿刺虚拟仿真与训练系统，包括：计算机主机 和力觉交互设备，在计算机主机上连接有显示器，所述主机的软件系统至少包括：具有皮 肤、软组织和骨骼的人体躯干部分的几何模型，经过模型约减后的躯干有限元模型，力觉交 互模块、躯干有限元模型变形计算模块、虚拟穿刺教学与训练模块，所述软件系统通过USB 接口与力觉交互设备连接，用于将力觉设备检测到的使用者做出的动作传输至计算机系 统，同时将系统计算得到的穿刺针作用力传输至力觉交互设备，最终由使用者感知。系统总 体架构如图2所示。
[0049] 一、包含脊柱骨骼和软组织的人体躯干的几何模型的构建，主要包括以下步骤：
[0050] 1)获取穿刺体位状态下人体腰部和部分胸部的CT图像。实际穿刺操作时，病人侧 卧位，背部与床板垂直，头向前胸屈曲，两手抱膝使其紧贴腹部，躯干呈弓形。其目的是使得 脊椎尽量后突以增宽脊椎间隙，便于进针。保持这样的体位，以2毫米为间隔，对图3中下 腰部的区域进行CT扫描。
[0051] 2)图像分割：采用主动轮廓线跟踪和马尔可夫随机场的图像分割算法对骨骼和 软组织的图像数据进行分割，对CT断层扫描数据采用人工识别的方法进行精细、准确的图 像分割。
[0052] 3)三维建模与矢量化：首先，采用计算机三维重建算法对分割数据进行三维重 建，建立面片模型，并进行平滑和简化处理。然后，对面片模型进行重采样，获取点云数据， 由点云数据构建自由曲面表示的三维数字模型，完成三维模型的矢量化。
[0053] 4)三维腰部躯干模型：系统所需几何模型包括：L1-L5腰椎（椎体、椎弓板、棘 突）、骨盆、骶骨、骶裂孔、骶角、棘上韧带、黄韧带、硬脊膜与珠网膜，以及由上述组织形成的 珠网膜下腔、硬膜外腔、骶后上棘、胸椎棘突、腰椎棘突。此外，也构建了皮肤和肌肉的几何 模型。如图3所示，其中背部区域需要CT扫描的部分。
[0054] 二、躯干有限元模型的建立和有限元模型的约减
[0055]1、下腰部躯干有限元模型的构建
[0056] 穿刺仿真涉及软组织的变形仿真。本系统采用有限元模型仿真软组织的变形。由 于系统使用力觉交互，而力觉交互通常要求变形仿真计算具有1000Hz以上的计算刷新率， 从而保证力觉交互的稳定性。为此，本发明采用了有限元分析、模态分析和模型约减技术来 达到这一要求。
[0057] 采用六面体单元，将躯干几何模型在三维空间中均匀划分，得到分配有不同物理 属性的有限元模型。单元的杨氏模量如下：皮质骨12000Mpa，松质骨lOOMpa，椎间盘纤维环 500MPa，椎间盘髓核4Mpa，韧带lOMpa，肌肉8Mpa，皮肤为6Mpa。
[0058] 2、有限元模型约减
[0059] 有限元模型的形变通常采用位移向量uGR3n(n为模型网格顶点的数目）表示，通 常n非常大，求解位移u的过程十分耗时。为了能实时对躯干有限元模型进行仿真求解，我 们用qGR1X3n)来代替位移向量u，q表示约减子空间的位移向量，r是约减空间的维 度，远远小于3n，通常为20以内。原始空间中的位移向量u和约减空间的位移向量q之间 的关系由方程（1)给出：
[0060] u = Uq (1)
[0061] 其中UGR3#是约减空间位移基矩阵（r彡3n)，在仿真过程中不发生变化。对有 限元模型进行约减，是实时仿真之前通过预计算获得U，从而大幅度减小有限元运动方程的 求解规模。
[0062] 获得U的方法如下：
[0063] ①采用商业有限元软件，对躯干的有限元模型进行模态分析，模态分析的输出是 cpGR30 ;
[0064] ②只保留前r阶（i= 1，2，…，r)的模态，用以构建线性模态的基矩阵 U= …屮r] ￡R3n>r，r取值 12〇
[0065] 基矩阵U就是模型约减基矩阵。
[0066] 三、躯干有限元模型变形计算模块
[0067]1、约减后的有限元模型的运动方程
[0068] 采用质量加权正交化方法使得U质量正交，则UtMU=I。可得到约减运动方程 (2)：
[0070] 方程⑵中Lt (q)表示约减内力，为约减空间中的外力，可以推导出约减刚 度矩阵表示如下：
[0072] 其中1是约减空间的刚度矩阵。方程⑵描述了可变形模型在约减空间中的运动 q=q(t)GRr0
[0073] 躯干有限元模型在原始空间中的运动也就可以通过u(t) =Uq(t)来获得。
[0074] 2、积分方法
[0075]采用隐式NewMark积分方法在约减空间中对约减动力学方程进行仿真。方程（4)、 (5)、（6)表示了隐式NewMark积分方法，其中A和y为参数，取值为A= 1/2,y= 1/4， 这样的取值可以获得较好的仿真稳定性。
[0079]方程（4)、（5)、（6)中q为约减空间中的位移向量，t为当前时刻，St为积分步长， M为质量矩阵，D为阻尼矩阵，R为变形造成的内力矩阵，f?为外力矩阵。矩阵的维数都为r。
[0080] 3、重建形变向量u
[0081] 方程（4)、（5)、（6)积分可以得到约减空间中的形变q，然而，必须根据q重建原始 空间的形变向量uG1^才可以更新有限元模型，实现躯干有限元模型的变形。
[0082] 重建形变向量u只需根据式（1)进行矩阵乘法即可。然而，由于本系统有限元规 模较大，实际的U的规模也很大，因此采用基于GPU的矩阵乘法来进行加速。
[0083] 基于CUDA实现了基于GPU的加速。采用了CUBLAS库方法来实现u=Uq，对大规 模矩阵与矢量的乘法运算问题进行加速。具体过程如下：
[0084] ①在显存中为矩阵分配存储空间；
[0085] ②将矩阵U、q从内存拷贝到显存中；
[0086] ③调用CUBLAS函数实现矩阵乘法运算u=Uq，得到u;
[0087] ④将结果u从显存拷贝回内存中。
[0088] 4、有限元模型的局部更新和采样
[0089] 躯干有限元模型的网格和顶点数量很大，没有必要更新和重建所有的网格和顶点 的形变。主要对两类情况进行局部更新和采样：
[0090] 为满足图形绘制的要求，只需要更新模型外表面的顶点的变形和位置。找出模型 外表面对应的有限元网格顶点，数量为ns，从约减空间基矩阵U中提取与外表面顶点对应 的行，并丢弃其余的行，得到外表面顶点的子空间基矩阵Usurf_eR3ns'外表面顶点的位移 为UsurfaceUsurfaced?其中（UsurfaceER)。
[0091] 对可能发生碰撞、干涉的局部区域（比如穿刺针附件的网格和顶点）进行更新，为 碰撞检测模块提供更新的几何信息。首先，找出有限元模型中局部兴趣网格顶点，数量为 ni，从子空间基矩阵U中提取与外表面顶点对应的行，并丢弃其余的行，得到局部兴趣顶点 的子空间基矩阵Uint_steRn1'然后，可以求出外表面顶点的位移为Uintarast=Uint_tq，其 中(uinteresteR3m) 〇
[0092] 四、力觉交互模块
[0093] 此模块主要包含三个子模块：碰撞检测子模块、碰撞响应子模块，力觉渲染子模 块。
[0094] 碰撞检测子模块的功能是，根据虚拟穿刺针和躯干各部分的相对位置进行碰撞检 测。
[0095] 碰撞响应子模块：根据碰撞检测的结果，计算出虚拟穿刺针受到的阻力。虚拟穿刺 针受到的阻力与虚拟穿刺针施加给躯体软组织上的力互为作用力与反作用力，把虚拟穿刺 针施加给躯体软组织上的力作为输入，传送给躯干有限元模型变形计算模块，从而计算出 躯干有限元模型变形后的形状和位移。
[0096] 力觉渲染子模块的功能是：接受碰撞检测子模块与碰撞响应子模块发送来的力觉 信息并将其发送至力觉交互设备，从而让操作者感受到虚拟穿刺产生的力反馈；同时，读取 力觉交互设备的操作手柄的位姿，并根据此位姿信息更新虚拟环境中的虚拟穿刺针的位姿 状态。
[0097] 1、碰撞检测子模块