一种针对虚拟微创血管介入手术的导丝建模方法
【专利摘要】本发明公开了一种针对虚拟微创血管手术的导丝建模方法，该方法包括：步骤1：对导丝进行离散化表示，得到导丝模型；步骤2：对所述导丝模型的参数进行初始化；步骤3：基于初始化的导丝模型参数，获得所述导丝模型的能量值；步骤4：基于所述步骤3得到的能量获得相应的力和力矩；步骤5：利用拉格朗日乘子式法实现所述导丝模型的不可拉伸约束，并计算出由该不可拉伸约束产生的约束力；步骤6：对所述步骤4和步骤5求得的力和力矩进行求和，并结合所述步骤4和步骤5求得的力和力矩，对所述导丝模型的参数进行更新；步骤7：循环调用步骤3～6来进行导丝的动态仿真。实验证明，本发明导丝模型能够逼真地、实时地模拟真实导丝的物理形变。
【专利说明】
一种针对虚拟微创血管介入手术的导丝建模方法
技术领域
[0001] 本发明涉及虚拟微创血管介入手术仿真领域，具体涉及一种针对虚拟微创血管介 入手术的导丝建模方法。
【背景技术】
[0002] 心血管疾病是全球致死率最高的疾病，微创血管介入手术是治疗该疾病的主要方 法，具有创伤小、恢复快等优势。但是，微创血管介入手术流程复杂，需要医生具有丰富的手 术操作经验。目前，传统手术培训主要在手术室内基于活体动物、尸体进行手术培训，培训 者很难随时随地进行手术训练。虚拟微创血管介入手术培训系统为培训者提供一个沉浸感 很强的虚拟培训场景，培训者可随时随地与虚拟环境交互，不断提高手术技能，提高医生培 训的效率，缩短培训周期。
[0003] 导丝建模作为虚拟血管介入手术培训系统的核心，一直是虚拟血管介入手术的研 究重点，当前存在很多建模方法，但都存在实时性和逼真度无法同时满足的缺点。基于此， 本发明提出了一种针对虚拟微创血管介入手术的导丝建模方法，该方法具有很好的实时性 和逼真度。
【发明内容】
[0004] 本发明实施的主要目的在于提供一个针对虚拟微创血管介入手术的导丝建模方 法，该方法基于弹性棒理论和拉格朗日乘子式法对导丝建模，能够同时满足实时性和逼真 度的需求。
[0005] 为实现上述目的，本发明提供了以下技术方案：
[0006] -种针对虚拟微创血管介入手术的导丝建模方法，所述方法包括以下步骤：
[0007] 步骤1:对导丝进行离散化表示，得到导丝模型；
[0008] 步骤2:对所述导丝模型的参数进行初始化；
[0009]步骤3:基于初始化的导丝模型参数，获得所述导丝模型的能量值；
[0010]步骤4:基于所述步骤3得到的能量获得相应的力和力矩；
[0011] 步骤5:利用拉格朗日乘子式法实现所述导丝模型的不可拉伸约束，并计算出由该 不可拉伸约束产生的约束力；
[0012] 步骤6:对所述步骤4和步骤5求得的力和力矩进行求和，并结合所述步骤4和步骤5 求得的力和力矩，对所述导丝模型的参数进行更新；
[0013] 步骤7:循环调用步骤3~6来进行导丝的动态仿真。
[0014] 可选地，所述导丝模型通过N个空间控制节点和N-1个四元数来表示，其中，N个空 间控制节点由对导丝离散化得到，N-1个四元数表示N个空间控制节点形成的N-1个中心线 段的右手正交坐标系。
[0015]可选地，所述空间控制节点的数量大于等于3,相邻空间控制节点之间的距离相同 或不同。
[0016] 可选地，所述步骤2中，需要进行初始化的导丝模型参数包括导丝的空间控制节点 坐标、四元数坐标、杨氏模量、剪切模量、导丝质量、导丝半径、导丝长度、导丝固有弯曲扭转 属性、重力大小和能量公式常量参数中的一种或多种。
[0017] 可选地，导丝模型头部的杨氏模量小于导丝模型导丝体部分的杨氏模量。
[0018] 可选地，所述步骤3中，所述能量值包括弯曲势能￡#、内摩擦耗能、粘滞力耗能 %和向量值约束补偿能Ce。
[0019] 可选地，所述步骤4中，首先将每个能量分别分解为最小单元，然后对每个最小单 元分别计算力和力矩，然后对每个能量对应的最小单元力和力矩分别进行求和得到每个能 量对应的力和力矩，最后将所有能量的力和力矩分别进行求和得到所有能量的合力以及力 矩之和。
[0020] 可选地，弯曲势能的最小单元为两个相邻四元数向量；内摩擦耗能的最小单元为 两个相邻空间控制节点；粘滞力耗能的小单元为两个相邻四元数向量；向量值约束补偿能 的最小单元为两个空间控制节点和空间控制节点间的一个四元数向量。
[0021] 可选地，所述步骤5进一步包括以下步骤：
[0022] 步骤51，通过下式得到拉格朗日乘子式系数λ:
[0024] 其中，W是空间控制节点的质量矩阵的逆，J是不可拉伸约束向量的雅克比矩阵，/ 表示雅克比矩阵的时间导数，％为空间控制节点的向量集合，衫ρ表示％的时间导数，f代表 能量产生的力与导丝所受外力之和，kdPkd分别代表弹簧和阻尼系数，C为不可拉伸约束的 向量形式，.亡表示约束的时间导数；
[0025] 步骤52,通过下式得到由不可拉伸约束产生的约束力Fif:
[0027]可选地，所述步骤6中，基于以下公式对导丝模型的参数进行更新：
[0033]其中，h代表时间步长，代表t时刻空间控制节点i的速度，nu表示空间控制节点 的质量，产表示t时刻导丝的受力，0表示t时刻的空间控制节点的位置，<4代表四元数j在 时亥叶时的角速度，I代表惯性张量，%#表示七+1!时刻的非单位四元数，//f表示t时刻的四 元数矩阵，0表示代表t时刻的单位四元数，II.qUI表示的绝对值，代表t时刻的 四元数j的临时力矩，其可以通过下式获得：
[0035] 其中，?表示力矩的中间变量，hJ"表示四元数矩阵的转置，T j表示求得的力矩之 和。
[0036] 与现有技术相比，本发明具有以下优点：
[0037] 1、本发明基于弹性棒理论对导丝进行建模，能够有效地模拟导丝的真实物理形 变，从而保证了手术器械模型的逼真度。
[0038] 2、本发明基于拉格朗日乘子式法保证了导丝模型不可拉伸的物理特性，进一步提 高了导丝模型的稳定性和逼真度。
[0039] 3、本发明通过导丝能量直接获得相应力和力矩，用于后面的导丝状态迭代，提高 了模型的计算速度，保证了模型的实时性。
【附图说明】
[0040] 图1为根据本发明一实施例的针对微创血管介入手术的导丝建模方法的流程图。
[0041] 图2为导丝的离散化表示。
[0042] 图3(a)_(b)为导丝模型在推拉操作下的误差变化情况。
[0043] 图4(a)_(d)为导丝模型在动脉弓的仿真效果。
[0044] 图5为导丝模型在冠脉分支中的仿真效果。
【具体实施方式】
[0045] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照 附图，对本发明进一步详细说明。
[0046] 图1为根据本发明一实施例的针对微创血管介入手术的导丝建模方法的流程图， 如图1所示，所述方法包括以下步骤：
[0047]步骤1:为了能够让计算机对导丝进行仿真，需要首先对导丝进行离散化表示，得 到导丝模型；
[0048]如图2所示，导丝可被离散成一系列的空间控制节点，如图2中所示的节点Gi，所述 空间控制节点主要存储该节点的X，y，z坐标，用来表示该空间控制节点在三维世界中的具 体位置。
[0049]在所述导丝模型中，将导丝离散成N个空间控制节点，这N个空间控制节点构成了 N-1个中心线段，如图2中的线段GwGi，为了实现导丝模型的弯曲和扭转操作，对每一个中 心线段附一个右手正交坐标系，如图2中的（Cl( j)，C2( j)，C3( j))，其中，Cl( j)、C2( j)、C3( j) 分别表示中心线段j的x，y，z轴。基于该坐标系，可以实现导丝的弯曲和扭转效果。在所述导 丝模型中，空间控制节点和右手正交坐标系通过向量值约芽
实现关联，其 中，G'表示导丝的切向量，所述向量值约束用来保证右手正交坐标系的C3轴与导丝的切向 量平行。由于真实导丝具有不可拉伸的物理特性，所以本发明中实现的导丝模型也具有不 可拉伸属性，于是，就可以假设在模型的迭代过程中，导丝长度不可拉伸，这样就可以将导 丝的切向量表示成
_这里，G' i为第i个切向量的离散化表示，li表示中心线段i 的长度。在本发明一实施例中，对于所述导丝模型，为了提高该导丝模型的计算速度，通过 四元数Q(j) = (qi (j)，q2 (j)，q3 (j)，q4(j))来表示右手正交坐标系，所述四元数与（ci (j)，C2 (j)，c3( j))之间的转换关系如下所示：
[0053]其中，91〇)^2〇)，阳〇)^4(」）分别代表第」个四元数的四个分量数值。于是，整 个导丝模型可以通过N个空间控制节点和N-1个四元数来表示。
[0054]步骤2:对所述导丝模型的参数进行初始化，以使所述导丝模型能够更加逼真地模 拟真实导丝；
[0055] 在得到导丝模型之后，需要对所述导丝模型的参数进行初始化。在该导丝模型中， 需要初始化的参数包括导丝的空间控制节点坐标、四元数坐标、杨氏模量、剪切模量、导丝 质量、导丝半径、导丝长度、导丝固有弯曲扭转属性、重力大小和能量公式中的常量参数。其 中，导丝的空间控制节点和四元数信息用来初始化导丝模型在空间中的位置和初始状态； 杨氏模量用来表示导丝用来抵制弯曲的能力大小，杨氏模量越大，导丝越能够提供较大的 恢复力。考虑到真实导丝分为导丝头部和导丝体，两个部分具有不同的杨氏模量，因此在对 导丝模型杨氏模量进行配置的过程中，导丝模型头部的杨氏模量应该小于导丝模型导丝体 部分的杨氏模量，这样导丝能够更容易的弯曲，也能够更容易模拟真实导丝的形变;剪切模 量用来表示导丝抵制扭转的能力大小，剪切模量越大，导丝在发生扭转时，提供的扭转恢复 力越大；导丝固有弯曲扭转属性用来控制导丝在稳定状况下，导丝所表现出来的弯曲扭转 效果，可以通过该参数改变导丝的外在弯曲形状。
[0056] 所述导丝模型中，所述空间控制节点的数量必须大于等于3,另外，相邻空间控制 节点之间的距离可以相同也可以不同。
[0057]步骤3:基于初始化的导丝模型参数，获得所述导丝模型的能量值，该步骤中，需要 求取的能量值主要包括弯曲势能￡#、内摩擦耗能β|、粘滞力耗能和向量值约束补偿能 Ce；
[0058]其中，弯曲势能用来表示导丝模型的弯曲扭转程度，可以通过下式获得：
[0060] 在这个公式中，L代表导丝的长度，Sk代表导丝的弯曲扭转应变，81和8 2代表导丝在 切面方向的弯曲应变，S3代表导丝在轴向的扭转应变大小，玲代表导丝的固有弯曲和扭转 属性，分别对应导丝在切面的两个弯曲固有应变和轴向的扭转固有应变，Kkk表示刚度张量 的主对角线元素。在本实施例的导丝模型中，假设导丝模型具有统一的横截面，于是导丝模 型的刚度张量就可以忽略非对角线元素，其刚度张量K可以通过其对角线元素 Kn、K22、K33来
表示，
，其中Ε为步骤1中提到的杨氏模量，G为剪切 模量，r为导丝半径。
[0061] 内摩擦耗能用来维持导丝模型的稳定性，其主要用来模拟空间控制节点间的摩 擦力，可以通过下式获得：
[0063] 在这里，γν代表内摩擦系数，Vrel代表空间控制节点相对平移速度。
[0064] 粘滞力耗能￡^.同样是用来维持导丝模型的稳定性，其主要用来模拟四元数间的粘 滞力，其可以通过下式获得：
[0066]其中，代表粘滞力系数，ω〇代表角速度，ω%表示相对角速度。
[0067]向量值约束补偿能用来维持向量值约束，在本实施例导丝模型中，基于补偿法 来实现该约束，可以通过下式获得：
[0069]其中，κ代表补偿系数。
[0070]步骤4:基于所述步骤3得到的能量获得相应的力和力矩，具体地，先将得到的每个 能量进行分解，对分解的部分能量分别进行力和力矩的求解，然后再将每个能量对应的力 和力矩进行求和统计；
[0071]在该实施例导丝模型中，力和力矩可以通过下式获得：
[0073]其中，η代表导丝空间控制节点和四元数向量的坐标集合，办表示η的时间导数， rin代表由这些能量产生的力和力矩。考虑到微分运算为线性运算，所以在对该公式进行求 解的过程中，可以分别对每个能量进行求解，然后再将其结果进行相加。在这里，r ιη包含力 和力矩Tin。
[0074] 具体地，在通过上式进行力和力矩求解的过程中，对上述四个能量分别计算相应 的力和力矩，计算过程中，将其它的能量赋值为零，并且在计算每个能量相应的力和力矩的 过程中，将每个能量拆分成导丝模型的最小单元，对每个最小单元计算出来相应的力和力 矩之后，再将其分别加到能量相应的全局力和力矩中，最后将所有能量的力和力矩分别进 行求和得到最终能量合力以及合力矩。其中，所述导丝模型的最小单元针对不同的能量有 不同的最小单元。其中，对于弯曲势能，其最小单元为两个相邻四元数向量;对于内摩擦耗 能，其最小单元为两个相邻空间控制节点;对于粘滞力耗能，其最小单元为两个相邻四元数 向量;对于向量值约束补偿能，其最小单元为两个空间控制节点和空间控制节点间的一个 四元数向量。
[0075] 步骤5:利用拉格朗日乘子式法实现所述导丝模型的不可拉伸约束，并计算出由该 不可拉伸约束产生的约束力；
[0076]真实的导丝具有不可拉伸性，为了实现该物理特性，本发明基于拉格朗日乘子式 法来进行实现，导丝的不可拉伸性是指导丝的长度不会随着操作而改变，转化到导丝模型 中就是指相邻的空间控制节点之间的距离保持不变。
[0077]在本实施例导丝模型中，将不可拉伸约束表示成向量形式C，在该导丝模型中，由 不可拉伸约束产生的内力Fit可以表示成JTA，其中，λ是乘子式系数，J是不可拉伸约束向 量的雅克比矩阵，可以通过|，％为空间控制节点的向量集合。在这里，λ可以通过 下式获得：
[0079] 其中，W是空间控制节点的质量矩阵的逆，j表示雅克比矩阵的时间导数，表示％ 的时间导数，f代表1^5与导丝所受外力之和，k4Pk d分别代表弹簧和阻尼系数，亡表示约束 的时间导数。通过求解该方程可以获得乘子式系数λ，之后可以通过/λ获得由不可拉伸约 束产生的约束力Fg e。
[0080] 步骤6:对所述步骤4和步骤5求得的力和力矩进行求和，并结合所述步骤4和步骤5 求得的力和力矩，通过欧拉显式方法对所述导丝模型的参数进行更新；
[0081] 通过步骤4和步骤5,可以获得所述导丝模型迭代所需要的全局合力F和合力矩T， 其中合力矩T即为步骤4中求得的Τιη，合力F包括f和Ft，其中所述导丝模型的参数更新主 要包括导丝的空间控制节点和四元数的信息更新。基于以下公式，可以对导丝模型的参数 进行更新：
[0087]其中，h代表时间步长，代表t时刻空间控制节点i的速度，nu表示空间控制节点 的质量，产表示t时刻导丝的受力，即t时刻求得的合力F，Gf:表示t时刻的空间控制节点的位 置，代表四元数j在时刻t时的角速度，I代表惯性张量，泛表示t+h时刻的非单位四元 数，$表示t时刻的四元数矩阵，表示代表t时刻的单位四元数，|| ||表示的绝 对值，代表t时刻的四元数j的临时力矩，其可以通过以下式子获得：
[0089] 其中，A表示力矩的中间变量，i/f表示四元数矩阵的转置，L表示求得的力矩和。
[0090] 步骤7:循环调用步骤3~6来进行导丝的动态仿真。
[0091] 如图3(a)_(b)所示，为本发明导丝模型在不同时间步长情况下，进行推拉操作后 最大误差的变化情况，图3(a)为本发明导丝模型在时间步长为lms的情况下，进行推拉操作 后最大误差的变化情况，图3(b)为本发明导丝模型在时间步长为0.1ms的情况下，进行推拉 操作后最大误差的变化情况。其中，每次迭代过程中，推拉操作的距离为1_，误差为2mm左 右，并且该误差能够在0.5秒以内减少到0.01mm以下。由图3(a)-(b)可知，本发明的导丝模 型具有良好的稳定性。
[0092] 图4(a)_(d)和图5展示了本发明导丝模型在动脉弓和冠脉分支中的形变效果，其 中，图4(a)显示导丝即将进入动脉弓，形变较小，为67mm，图4(b)显示导丝开始进入动脉弓， 形变变大，为174mm，图4(c)显示导丝进入动脉弓，开始产生较大形变，为277mm，图4(d)为导 丝完全进入动脉弓，产生180度大形变，为352mm。由图4(a)-(d)和图5可知，本发明导丝模型 具有很好的稳定性以及逼真的形变效果。
[0093] 以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详 细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡 在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保 护范围之内。
【主权项】
1. 一种针对虚拟微创血管手术的导丝建模方法，其特征在于，所述方法包括W下步骤： 步骤1:对导丝进行离散化表示，得到导丝模型； 步骤2:对所述导丝模型的参数进行初始化； 步骤3:基于初始化的导丝模型参数，获得所述导丝模型的能量值； 步骤4:基于所述步骤3得到的能量获得相应的力和力矩； 步骤5:利用拉格朗日乘子式法实现所述导丝模型的不可拉伸约束，并计算出由该不可 拉伸约束产生的约束力； 步骤6:对所述步骤4和步骤5求得的力和力矩进行求和，并结合所述步骤4和步骤5求得 的力和力矩，对所述导丝模型的参数进行更新； 步骤7:循环调用步骤3~6来进行导丝的动态仿真。2. 根据权利要求1所述的方法，其特征在于，所述导丝模型通过N个空间控制节点和N-1 个四元数来表示，其中，N个空间控制节点由对导丝离散化得到，N-1个四元数表示N个空间 控制节点形成的N-1个中屯、线段的右手正交坐标系。3. 根据权利要求2所述的方法，其特征在于，所述空间控制节点的数量大于等于3,相邻 空间控制节点之间的距离相同或不同。4. 根据权利要求1所述的方法，其特征在于，所述步骤2中，需要进行初始化的导丝模型 参数包括导丝的空间控制节点坐标、四元数坐标、杨氏模量、剪切模量、导丝质量、导丝半 径、导丝长度、导丝固有弯曲扭转属性、重力大小和能量公式常量参数中的一种或多种。5. 根据权利要求4所述的方法，其特征在于，导丝模型头部的杨氏模量小于导丝模型导 丝体部分的杨氏模量。6. 根据权利要求1所述的方法，其特征在于，所述步骤3中，所述能量值包括弯曲势能 反韓、内摩擦耗能％、粘滞力耗能C訂日向量值约束补偿能Ce。7. 根据权利要求1所述的方法，其特征在于，所述步骤4中，首先将每个能量分别分解为 最小单元，然后对每个最小单元分别计算力和力矩，然后对每个能量对应的最小单元力和 力矩分别进行求和得到每个能量对应的力和力矩，最后将所有能量的力和力矩分别进行求 和得到所有能量的合力W及力矩之和。8. 根据权利要求7所述的方法，其特征在于，弯曲势能的最小单元为两个相邻四元数向 量；内摩擦耗能的最小单元为两个相邻空间控制节点；粘滞力耗能的小单元为两个相邻四 元数向量;向量值约束补偿能的最小单元为两个空间控制节点和空间控制节点间的一个四 元数向量。9. 根据权利要求1所述的方法，其特征在于，所述步骤5进一步包括W下步骤： 步骤51，通过下式得到拉格朗日乘子式系数λ:其中，W是空间控制节点的质量矩阵的逆，J是不可拉伸约束向量的雅克比矩阵，/表示 雅克比矩阵的时间导数，％为空间控制节点的向量集合，嗦Ρ表示％的时间导数，f代表能量 产生的力与导丝所受外力之和，ks和kd分别代表弹黃和阻尼系数，C为不可拉伸约束的向量 形式，€表示约束的时间导数； 步骤52,通过下式得到由不可拉伸约束产生的约束力Flije;10.根据权利要求1所述的方法，其特征在于，所述步骤6中，基于W下公式对导丝模型 的参数进行更新：其中，h代表时间步长，代表t时刻空间控制节点i的速度，nil表示空间控制节点的质 量，F嗦示t时刻导丝的受力，Gf表示t时刻的空间控制节点的位置，代表四元数j在时刻 t时的角速度，I代表惯性张量，Qi + h表示t+h时刻的非单位四元数，妇I表示t时刻的四元数 矩阵，0表示代表t时刻的单位四元数，I暢+i表示巧+'，.的绝对值，巧代表t时刻的四元 数j的临时力矩，其可W通过下式获得：其中，表示力矩的中间变量，表示四元数矩阵的转置，Tj表示求得的力矩之和。
【文档编号】G06F19/00GK106096265SQ201610404708
【公开日】2016年11月9日
【申请日】2016年6月8日
【发明人】谢晓亮, 侯增广, 高占杰, 边桂彬, 奉振球, 郝剑龙, 王莉, 周小虎, 程笑冉
【申请人】中国科学院自动化研究所