基于麦夸特算法的径向基无网格软组织数据的力反馈模型建模方法与流程

本发明涉及虚拟手术软组织力反馈模型建模，用于实现在虚拟手术中生物软组织力反馈模型与操作者的精确性实时性交互。
背景技术：
：一个理想的生物软组织模型能够准确的将软组织位移情况及受力大小即时反馈给操作者，这就要求生物软组织模型具有高精确性及实时性。而在虚拟手术中，高效的软组织力反馈物理模型建立对虚拟手术系统的成功起着至关重要的作用。目前，基于物理驱动的软组织模型日渐成熟，大致可分为两类：一种是基于网络结构，主要有质点-弹簧模型和有限元模型。质点-弹簧模型(MS)的优点是原理简单、计算量比较小和实时性好，缺点是形变具有局部性质，并且对模型中弹簧刚度系数设置合理的材料参数是没有理论依据，很容易造成软组织在形变仿真的时候，出现失真问题。有限元模型(FEM)的优点是计算精度高，但是由于求解过程中需要组装刚度矩阵，计算量比较大，如果进行大形变或者切割的时候，网格模型的拓扑结构发生复杂的变化，使得形变实时性差；另一种是无网格结构，主要方法有粒子系统、有限球方法和无网格方法。HieberSE等基于线弹性力学模型，利用光滑粒子流(SPH)无网格法，不需要网格拓扑结构，实时性高，可以满足软组织大范围的形变，但存在稳定性差以及不具备粘弹性等问题，该方法常用于快速粒子流或气体等情况的仿真。比如，流血、爆炸和烟雾等。SuzukiS的填充球模型，是一种基于几何驱动的无网格法，不包含软组织的生物特性，不能满足虚拟手术仿真的高精确度要求。无网格法是机械工程领域中的数值方法，近年来发展迅速，出现了如无网格伽辽金法(EFG)、径向基函数法(RBF)和多尺度重构核粒子方法(MRKP)等多种无网格法。它采用一组离散的场节点来构造问题域和边界，通过插值来生成每个场节点的形函数，在不需要任何点与点之间拓扑关系的情况下近似逼近问题域内未知的场函数，最后求解离散系统方程，得到每个场节点的位移、应力和应变。无网格法常用于求解悬臂梁等机械固体受力后的形变，在软组织形变方面并没有形成一套成熟的理论系统。由于无网格法本身存在的特性使得其可以应用与软组织形变仿真中。其优点：第一，自适应性强，采用点云结构模型，不需要点与点之间复杂的拓扑结构，可适用于大形变或者切割。第二，建立模型只需要节点的信息。第三，计算结果光滑连续，适合力反馈设备。但是由于无网格存在计算效率低下的问题，所以在实时性方面并不是很理想。综上所述可知，基于网络结构的质点-弹簧模型(MS)结构简单实时性好，但其在形变上具有局部性质，精确性不足。有限元计算精度高但计算量大形变实时性差。而无网格结构具有非常高的精确性及自适应性满足虚拟手术中反馈精确要求，但是无网格法依然面临计算效率低下，在实时性方面有所缺陷。技术实现要素：本发明的目的是提出一种基于麦夸特算法的径向基无网格软组织数据的力反馈模型建模方法。针对现有生物软组织物理模型缺陷，考虑生物软组织各向异性特性，由径向基无网格法一次性获得以软组织受力点为中心，面上有限个等间距点受力位移数据，使用麦夸特法拟合软组织节点力-位移模型，再以点至线，以线带面，建立生物软组织力反馈模型。在保留径向基无网格法精确性的同时，又大大简化计算量，保证了交互的高实时性，实现虚拟手术中生物软组织力反馈模型与操作者的精确性实时性交互。本发明是通过以下技术方案实现的。生物软组织第hαβ节点有如下系统：以受力点O为坐标原点，其中x、y、z分别为hαβ节点横轴方向位移量、纵轴方向位移量及垂直水平面方向位移量，F为受力点O反馈力大小，为随机误差项。利用径向基无网格法获取以受力点O为中心，间距为R的H个矩形点阵节点的软组织输出力F(k)和输入位移x(k),y(k),z(k)数据。(在现实手术中生物组织受力形变的组织直径范围在1～4cm左右)。α×β＝H考虑生物软组织力与形变的特点可将系统写成：其中x、y、z分别为hαβ节点横轴方向位移量、纵轴方向位移量及垂直水平面方向位移量，p为系统待估参数。开始取α＝1，β＝1。步骤(1)由径向基无网格法获得各节点软组织力与形变输入输出数据，采样当前第hαβ节点输入输出数据。获取L组软组织数据。步骤(2)为避免过度拟合，利用交叉验证法先确定hαβ节点系统自变量阶数na、nb、nc，而确定了自变量阶数的同时也确定了参数p的个数m。步骤(3)设置待估参数p初值设定阻尼因子d＝10ud(0)(d(0)>0)。步骤(4)设定迭代次数u＝0，进行第一次迭代，利用L组软组织数据、参数初值p(0)及阻尼因子初值d(0)通过麦夸特算法参数估计解得p值，将求得参数p代入系统原函数式，计算残差平方和J(0)。步骤(5)将步骤(4)计算获得参数p赋予p(0)，即p(0)＝p，令u＝1即d＝10d(0)进行第二次迭代，再执行步骤(4)解得新的参数p并代入系统原函数式，计算新的残差平方和J(1)。比较J(1)与J(0)，增加迭代次数，直至J(1)<J(0)为止。步骤(6)继续增加阻尼因子d的数量级即u＝u+1，反复执行步骤(5)，直至新计算获得参数p与上一次计算参数即p(0)之差绝对值|Φ|≤δ允许最小误差为止。此时hαβ节点系统参数估计完成。步骤(7)令α＝α+1回到步骤(1)重新执行，即对节点矩阵纵向建模，当时，步骤(7)结束。步骤(8)令α＝1，β＝β+1回到步骤(1)重新执行，即对节点矩阵横向建模，当时，步骤(8)结束。至此建立了H个矩阵节点输入输出系统，即建立以受力点O为中心面积为HR2生物软组织力反馈矩阵模型。本发明的优点：以麦夸特法快速拟合径向基无网格法获得的软组织输入输出数据，建立生物软组织力反馈模型，该方法结构简洁、计算方便，在保留径向基无网格法精确性的同时，又大大简化计算量，保证了交互的高实时性，实现虚拟手术系统软组织力反馈模型高效建模。附图说明图1为发明流程框图。图2为力反馈示意图。具体实施方式本发明将通过以下实施例作进一步说明。生物软组织第hαβ节点有如下系统：以受力点O为坐标原点，其中x、y、z分别为hαβ节点横轴方向位移量、纵轴方向位移量及垂直水平面方向位移量，F为受力点O反馈力大小，为随机误差项。利用径向基无网格法获取以受力点O为中心，间距为R的H个矩形点阵节点的软组织输出力F(k)和输入位移x(k),y(k),z(k)数据。(在现实手术中生物组织受力形变的组织直径范围在1～4cm左右)。α×β＝H考虑生物软组织力与形变的特点可将系统写成：令原系统函数式可以写成其中x、y、z分别为hαβ节点横轴方向位移量、纵轴方向位移量及垂直水平面方向位移量，p为系统待估参数。通过径向基无网格算法提供第hαβ节点的M组软组织力与形变数据：(xi1,xi2,...,xina;yi1,yi2,...,yinb;zi1,zi2,...,zinc;Fi),i=1,2,...,M]]>开始取α＝1，β＝1。步骤1采样当前由径向基无网格法获得的生物软组织第hαβ节点输入输出数据。获取L组软组织数据：(xi1,xi2,...,xina;yi1,yi2,...,yinb;zi1,zi2,...,zinc;Fi),i=1,2,...,L]]>步骤2为避免过度拟合，先确定hαβ节点系统自变量阶数na、nb，而确定了自变量阶数的同时也确定了参数p的个数。确定自变量阶数具体过程描述如下：(a)利用交叉验证法，有L个样本，对于na阶水平方向位移自变量x，设纵轴位移、垂直水平位移自变量y，z为常量C、D，对x做L次拟合，将第i(i＝1,2,…L)次观测值剔除，计算残差平方和,寻找到使其最小的阶数A即其中b-i,a为删掉第i个观测之后的估计值。(b)设定横向位移x、垂直水平位移z为常量C、D，同理可求得nb。(b)设定水平位移自变量x、y为常量C、D，同理可求得nc。(d)初步描述系统为：F=p1+p2x1+p3x2+...+pna+1xna+pna+2y1+pna+3y2+...+pna+nb+1ynb+pna+nb+nc+1znc---(1)]]>步骤3设置待估参数p初值设定阻尼因子d＝10ud(0)(d(0)>0)；设定迭代次数u＝0即d＝d(0)，对于阻尼因子d＝0即高斯-牛顿法，在高斯-牛顿法中参数初值选择的合适与否,将决定迭代计算的工作量,甚至迭代成功与否。即使有较好的选择初值的方法,也难免出现发散的情况，而麦夸尔特法由于阻尼因子的引入，在一般的初值选择条件下都可收敛于所求结果。步骤4设定迭代次数u＝0，进行第一次迭代，利用L组软组织数据、参数初值p(0)及阻尼因子初值d(0)通过麦夸特算法参数估计解得p值，将求得参数p代入系统原函数式，计算残差平方和J(0)。麦夸特求解待估参数具体过程描述如下：(a)将f(xi,yi,zi,p)在p(0)处按泰勒级数展开并略去二次及二次以上的项得：f(xi,yi,zi,p)≈f(xi,yi,zi,p)+∂f(xi,yi,zi,p)∂p1|p=p(0)(p1-p1(0))+∂f(xi,yi,zi,p)∂p2|p=p(0)(p2-p2(0))+...∂f(xi,yi,zi,p)∂pm|p=p(0)(pm-pm(0))---(2)]]>(b)式(2)中除待估参数(p1,p2,…pm)外，都是已知。对此使用最小二乘残差平方原理：J=ΣiL{Fi-[f(xi,yi,zi,p(0))+Σj=1m∂f(xi,yi,zi,p(0))∂pj|p=p(0)(pj-pj(0))]}2+dΣjm(pj-pj(0))2---(3)]]>其中d≥0为阻尼因子，当d＝0时就是牛顿-高斯法，他对迭代初值的选择要比麦夸特法更加严格。(c)通过将残差平方J分别对p1,p2,…pm一阶偏导设为0，以达到将其最小化的目的：∂J∂pk=2Σi=1L[Fi-f(xi,yi,zi,p(0))+Σj=1m∂f(xi,yi,zi,p(0))∂pj|p=p(0)(pj-pj(0))]·∂f(xi,yi,zi,p(0))∂pk|p=p(0)+2d(pk-pk(0))=0,k=(1,2,...,m)---(4)]]>(d)令ajk=∂f(xi,yi,zi,p)∂pk|p=p(0)·Σi=1L∂f(xi,yi,zi,p)∂pj|p=p(0)ajy=∂f(xi,yi,zi,p)∂pj|p=p(0)·Σi=1L(f(xi,yi,zi,p)-Fi)j=1,2,...,m;k=1,2,...,m---(5)]]>可将式(4)写成：(a11+d)(p1-p1(0))+a12(p2-p2(0))+...+a1m(pm-pm(0))=α1ya21(p1-p1(0))+(a22+d)(p2-p2(0))+...+a2m(pm-pm(0))=α2y...am1(p1-p1(0))+am2(p2-p2(0))+...+(amm+d)(pm-pm(0))=amy---(6)]]>(f)解得P=p1p2···pm=p1(0)p2(0)···pm(0)+a11+da12...a1ma21a22+d...a2m···am1am2...amm+d---(7)]]>步骤5将步骤4计算获得参数p赋予p(0)，即p(0)＝p，令u＝1即d＝10d(0)进行第二次迭代，再执行步骤4解得新的参数p并代入系统原函数式，计算新的残差平方和J(1)。比较J(1)与J(0)，增加迭代次数，直至J(1)<J(0)为止。具体过程描述如下：比较：若J(1)<J(0)，步骤5结束；若J(1)≥J(0)，取u＝2即d＝100d(0)，再执行步骤4重解p并计算新残差平方J(1)，若新计算的J(1)<J(0)步骤5结束，反之则再取u＝3即d＝1000d(0)再执行步骤4…，如此反复迭代直至J(1)<J(0)为止。在式(6)中，aij为定值，所以当d越大时Φ＝|p-p(0)|越小，残差平方和J收敛，而d相对初值d(0)增加的越大即u越大迭代的次数越多。选择的界限是看残差平方和是否下降。步骤6继续增加阻尼因子d的数量级即u＝u+1，反复执行步骤5，直至新计算获得参数p与上一次计算参数即p(0)之差绝对值|Φ|≤δ允许最小误差为止，使得估计参数波动范围忽略不计。此时hαβ节点系统参数估计完成。步骤7令α＝α+1回到步骤1重新执行，即对节点矩阵纵向建模，当时，步骤7结束。步骤8令α＝1，β＝β+1回到步骤1重新执行，即对节点矩阵横向建模，当时，步骤8结束。至此建立了H个矩阵节点输入输出系统，即建立以受力点O为中心面积为HR2生物软组织力反馈矩阵模型。当前第1页1 2 3