虚拟手术中带摩擦的碰撞力计算及其物理仿真研究

欧惠棠 李晋芳 莫建清

(广东工业大学机电工程学院 广州 510006)

引言

虚拟医疗是应用虚拟现实(virutalreatility,VR)技术,构建虚拟的人体模型器官以及手术，提高虚拟实境的真实感。借助于虚拟外设，可以使人们更逼真地学习医疗知识以及治病救人的虚拟现实应用，在未来有着非常乐观的前景。其中，手术训练与虚拟现实相结合的虚拟手术是现代手术训练的新趋势。虚拟手术仿真利用医学生物信息和计算机仿真技术以及其他仿真模拟器械，为医疗工作者提供一个虚拟且逼真的手术场景，使其借助这些虚拟场景进行手术计划训练和研究，积累经验，掌握技巧，以便在实际的手术过程中完成复杂的手术操作，为病患减少手术后的痛苦和手术风险[1]。通过这一手段，无疑可使训练更真实、准确、可靠。

国内外对于虚拟手术仿真中碰撞力的研究，着重于变形体形变技术的实现上，如Montgomery等提出的现在非常流行的弹簧-质点网络(spring-mass networks)[2]，Felippa等提出的基于共旋转模型的有限元方法(FEM)[3]，吴涓等的一种基于沿径向方向分割为呈同心圆分布的弹簧-质点模型及实时力学响应算法[4]，陈卫东等提出的虚拟体弹簧的变形算法[5]等。这些研究仅仅着重于变形体内部力的计算上，并没有系统地给出碰撞检测、碰撞力的计算、时间积分等整个仿真流程。国内的研究仅仅给出当模型发生形变时内力与粒子群位移之间的关系，却没有给出与手术器械模型发生碰撞时如何计算出碰撞力并使变形体发生变形。国外Hasegawa等提出惩罚方法来求解碰撞接触力[6]，但会导致一系列不稳定性问题。Anitescu等提出k边形金字塔算法[7]，但时间复杂度较大，导致运算时间过长。

本研究提出虚拟手术训练中带摩擦的碰撞力计算及其仿真方法，阐述粒子系统的建立、碰撞检测得到碰撞力方向、隐式欧拉积分、碰撞约束方程(线性互补问题)、摩擦约束方程的构建、求解方程得出碰撞力大小这一系列的完整流程。该方法适用于刚体和变形体(弹簧质点网络和有限元)物理模型仿真，不仅实现了碰撞功能，而且还添加了摩擦特性，能够在一定接触点数目的场景下满足交互的真实性、稳定性和实时性的要求。

在计算机图形学中，物理仿真通过粒子系统[8]实现。粒子系统可以建立刚体和变形体等物理对象，并能实现其物理行为。本文所述的碰撞研究对刚体和变形体均适用。对于碰撞物理仿真的流程，已有相关的研究。下面通过约束动力学构建碰撞与摩擦相关约束方程并求解，然后通过实验验证是否计算正确和性能是否符合要求。

1 材料和方法

1.1 材料

采用自行设计的虚拟手术仿真系统，加入带摩擦的碰撞力计算算法，进行蓝钳模型和半月板模型的碰撞测试。实验环境为Win10 64 位系统，处理器为Intel(R) Core(TM) i5-3470CPU @ 3.20 GHz，显卡为NVIDIA Quadro 600。先分别通过蓝钳与半月板模型的单接触面和双接触面碰撞进行10次实验，验证所计算出的碰撞力方向和响应时间是否达到要求，再验证100次测试实验中，单接触面接触碰撞从1个接触点到45个接触点同时碰撞的平均碰撞时间。

1.2 方法

计算模型的碰撞接触力的基本仿真流程如图1所示。

图1 碰撞仿真流程Fig.1 Schematic diagram of collision simulation

为了计算接触力，必须先通过碰撞检测，检测出潜在的碰撞对象，当模型发生碰撞后，精确计算出碰撞点的位置和接触力的方向(大小未知)，然后构建相应的约束方程，再利用相应的数值方法进行求解。求解后，通过图形界面验证计算出的反馈力的方向和大小是否符合要求，通过20次的碰撞测试，并作出碰撞力的曲线，验证接触力的响应时间和稳定性是否符合要求。通过细分半月板模型的网格并通过100次的实验，验证单接触面接触从1个接触点到50个接触点同时接触的平均碰撞时间。

1.2.1碰撞检测

如图2所示，通过碰撞检测，可以得出两个即将碰撞的模型的潜在碰撞点P和Q。可以任意选取法向量，本研究选取模型D1上P点对应的朝内的法向量，记为n。

图2 两个物体之间的碰撞Fig.2 The collision between two objects

定义两个模型在P点上的空隙为

(1)

fn21(P)+fn12(Q)=0

(2)

由Signorini法则[9]，对于碰撞行为，δn(P)与fn21(P)非零并且具有互补关系，记为

0≤δn(P)⊥fn21(P)≥0

(3)

互补关系是指δn(P)与fn21(P)中必有一个为零[10]。

1.2.2插值

在计算机图形系统中，模型表面为三角面片构成的网格。模型之间的碰撞化归为两个三角形之间的标准碰撞，如图3所示。碰撞检测出来的碰撞点可能位于三角形的内部或边线上，但由于在粒子系统中，粒子通常是三角面片的顶点，因此必须把碰撞点的碰撞力等效到三角形的顶点上。

图3 两个三角形碰撞的全部5种情况Fig.3 All five cases of the collision between two triangles

根据解析几何相关知识，三角形所在平面的任意点都能表示为顶点的加权平均值，其中的权称为重心坐标。在图4中，对于点P，有

P=αA1+βB1+γC1

(4)

式中，α、β、γ为重心坐标分量。

式(4)中两边取微位移，有

u(P)=αu(A1)+βu(B1)+γu(C1)

(5)

图4 碰撞点在三角形内部Fig.4 The contact point located in the interior of atriangle

由虚功原理，有

uT(P)fP=uT(A1)fA1+uT(B1)fB1+uT(C1)fC1

(6)

从而

uT(A1)(fA1-αfP)+uT(B1)(fB1-βfP)+

uT(C1)(fC1-γfP)=0

(7)

得出

fA1=αfPfB1=βfPfC1=γfP

(8)

1.2.3构建线性互补问题

在图2中，设D1、D2中所有粒子的位置分别为q1、q2，三角形碰撞面片的顶点与q1、q2有一定的映射关系，称为接触蒙皮(contact skinning)，不详细讨论。接触蒙皮能使得两个模型之间三角面片顶点的约束转化为模型之间粒子群的约束。因此，δn(P)≥0可表示为

Ψ(q1,q2)≥0

(9)

称为单面约束，Ψ仅与位置有关，称为完整约束[11]。

粒子系统的物理仿真实质上是求解牛顿第二定律常微分方程(ODE)有

(10)

式中，fext为外部力，fint为内部力，fcon为约束力。

令fext+fint=f(q,v)，fcon=HTλ，其中q为粒子群的位置，v为粒子群的速度，H为约束力的方向，λ为约束力的大小。

使用隐式欧拉方法来求解常微分方程，有

vt+h-vt=hat+h

(11)

qt+h-qt=hvt+h

(12)

(13)

式中，t为上一个时间步的时刻，h为隐式欧拉方法时间步的步长。

(M-hB-h2K)dv=hf(xt,vt)+h2Kvt+HTλh

(14)

式中，dv=vt+h-vt。

令A=M-hB-h2K,b=hf(qt,vt)+h2Kvt，新的λ=λh。则有Adv=b+HTλ，从而对两个粒子群q1、q2之间的约束，有

(15)

若取λ=0，求得的解为没有约束的运动状态，称为自由运动解，为dvfree1和dvfree2。

事实上，dv=dvfree+dvcor，因此，

(16)

最后得到

(17)

对于线性互补问题，其求解过程可由图5直观表示。

1.2.4建立摩擦约束

图5 线性互补问题求解。(a)未发生碰撞；(b)发生碰撞Fig.5 Linear complementar yproblem solving.(a)Collision does not occur； (b)Collision occurs

对于摩擦约束，仅知道摩擦力位于所在碰撞点的切平面上，而在切平面的哪个方向未知，需要求解。令摩擦力为ft，滑动摩擦位移为δt，是一个二维向量。

碰撞约束方程与摩擦约束方程合并，得到总约束方程为

(18)

其中，接触约束方程为线性互补问题，有

0≤fn⊥δn≥0

(19)

对于摩擦约束方程，要满足干摩擦法则(Coulomb’s law)，有

对于摩擦约束，其求解可由图6直观表示。

图6 摩擦约束求解。(a)静摩擦；(b)动摩擦Fig.6 Friction constrain tsolving.(a) static friction； (b) dynamical friction

1.2.5建立总约束方程

对于接触面约束，总的约束方程为式(18)。对于复杂的手术交互，除了单接触面碰撞之外，还会涉及双接触面碰撞，如软组织夹取交互。若双接触面碰撞仍按照单接触面碰撞的情况处理，则一个接触面碰撞响应会影响另一个接触面的碰撞响应，导致两个接触面不能稳定地产生碰撞压力。灵活地改变约束方程，能使约束方程仍然适用于双接触面的情形。

如果碰撞检测出软组织的同时与两个接触面发生碰撞，则约束方程发生改变。如图7所示，图中三角面片代表接触区域。(a)中软组织模型与夹子模型的一个面接触，此时采用单接触面约束求解方法。(b)中软组织模型与夹子模型的两个面同时接触，则原来的单接触面状态转变为双接触面状态，使用双接触面约束求解方法。

图7 接触状态的转变。(a)单接触面接触；(b)双接触面接触Fig.7 Change of contact status.(a)Single contact surface collision； (b)Double contact surface collision

软组织模型碰撞点所在的三角面片顶点与夹子碰撞面建立新的线性互补，有

0≤δn⊥(ξfn+∑ξadjfnadj)≥0

(22)

建立基于顶点领域的摩擦方程，有

δt=0⟹‖ft‖<μ‖ξfn+∑ξadjfnadj‖

(23)

(24)

式中，fnadj为顶点的邻接顶点的法向约束力，ξ为权重。

当fn≤0，ξ=0；当fn>0，ξ=fn/N(N为顶点及其邻接顶点的数目)。

由于新的约束方程允许fn≤0，并且把三角面片顶点及其邻接顶点法向力的综合作用考虑进来，可以有效提高交互的稳定性，并且实质上更符合于现实的接触情况。

从而双接触面接触总的约束方程为

(25)

式中，m为接触的三角面片顶点的数目。

1.2.6求解总约束方程

对于式(25)，使用类高斯-赛德尔算法(Gauss-Seidel-like algorithm)求解。在粒子仿真计算力学领域[13]中，该算法能保证较好的收敛性。先对接触逐个求解，然后不停地迭代，直到数值收敛。对于m个瞬时摩擦接触，先考虑其中一个摩擦接触α，对其进行求解，有

(26)

求出fα，其他接触也以此进行计算，然后不断迭代求出新的fα，直到数值收敛。

(27)

本研究仅给出双接触面接触约束求解算法，单接触面接触求解算法与双接触面接触类似，并且要简单，由于篇幅所限，这是不再赘述。由于双接触面接触要求查询三角面片顶点的邻接顶点，故模型的网格要求为流形网格，并采用半边数据结构，可以使查询操作的时间复杂度为O(1)。算法如下：

Input:(δfree)3 m×1,[W]3 m×3m

Output:(f)3 m×1

设置阈值ε1、ε2、ε3

k=0

fori=1…mdo

(fi(0))(3×1)=0

Λmin/max=eig([Wii]tt) 每个([Wii]tt)均有两个特征值

end

repeat

k=k+1

foreachi=1…mdo

(δ)(3×1)=(δfree i)(3×1)

end

foreachj=1…i-1do

(δ)(3×1)+=[Wij]3×3(fj(k))3×1

end

foreachj=i…mdo

(δ)(3×1)+=[Wij](3×3)(fj(k-1))(3×1)

end

if(fi)nadj>ε1then

(fi(k))t=(fi(k-1))t-(δ)t/Λi

else

(fi(k))=0

end

foreachi=1…mdo

if(fi)n>ε3thenξ=1

elseξ=0

end

(fi)nadj=ξ(fi)n

((fik)nadj)l×1=adjacent Vertices Normal Forces(i)

foreachs=1…ldo

if(fis)nadj>ε3thenξs=1

elseξs=0

end

(fi)nadj+=ξs(fis)nadj

end

(fi)nadj/=l+1

if(fi)nadj≤ε3

(fi(k))t=0

continue

end

if‖(fi(k))t‖>μ(fi)nadjthen

(fi(k))t×=μ(fi)nadj/‖(fi(k))t‖

end

end

1.2.7实验方法

在碰撞试验中，类高斯-赛德尔算法收敛阈值ε1、ε2、ε3取0.001，摩擦系数μ取0.5。半月板模型设定杨氏模量为59 MPa，泊松比为0.45[14]。实验方案为先分别进行10次单接触面和双接触面接触碰撞实验。单接触面接触时，蓝钳模型简单地碰撞半月板模型；双接触面接触时，蓝钳模型夹取碰撞模型。碰撞时，分析碰撞力的大小和方向、响应时间、稳定时候的起伏是否符合要求。然后，再细分半月板模型的网格，使得碰撞时的接触点增加，进行100次实验，计算1个接触点到45个接触点同时接触的平均运算时间，其运算时间的记录方法如图8所示。

图8 碰撞时间测试流程Fig.8 Diagram of the collision time test

1.2.8统计学分析

单接触面接触中，以膝关节模型的正方向(髌骨方向)为基准，分别在竖直和水平方向从-90°、45°、0°、45°、90°各进行1次碰撞实验，计算平均的响应时间和稳定时的起伏。在双接触面接触中，从水平方向-90°、45°、0°、45°、90°进行10次夹取试验，计算平均的响应时间和稳定时的起伏。在多接触点平均碰撞力运算时间的测试中，为了能使足够的接触点接触，接触的方向任意选取，然后计算平均计算时间。

2 结果

图9为笔者构建的膝关节镜手术场景，在场景中，蓝钳模型与半月板模型发生碰撞。碰撞接触力通过上述方法得到计算并绘制，由图可知接触力的方向符合实际情况。

图9 场景中碰撞力的绘制Fig.9 Contact force rendering in the scene

图10为10次单接触面碰撞实验中其中一次的碰撞力的曲线。由图可知，模型在1和2 s的时候发生了碰撞，x、y、z方向力的大小变化趋势无显著差异，开始碰撞时，接触力急剧上升，随着半月板模型的变形行为趋于稳定，碰撞力也趋于稳定。在2 s的时候，蓝钳模型在z轴负方向前进0.2个单位长度，稳定后碰撞力增加。在10次的实验中，模型碰撞的平均响应时间为0.02 s，接触力稳定后的平均起伏范围为±0.02个单位。

图10 碰撞力的分析Fig.10 Analysis of the contact force

图11 通过摩擦力夹取变形体Fig.11 Grasping a deformable object by friction

图11为半月板模型夹取的效果。在不同的时间段，其状态从左到右、从上到下分别为：准备夹取阶段、夹取阶段、向右拖动阶段、向下拖动阶段、释放阶段。

图12为10次夹取实验中的一次蓝钳下颚碰撞力曲线。在6.12 s的时候，蓝钳夹取半月板，接触力在y方向发生突变，随后趋于稳定；在8.78 s的时候，蓝钳拉动半月板运动，摩擦使接触力在z轴方向上发生突变；在10.78 s的时候，蓝钳释放半月板，接触力回复为零。在10次的实验中，模型碰撞的平均响应时间为0.02 s，接触力稳定后的平均起伏范围为±0.02个单位。

图13为100次单接触面碰撞下1～45个接触点同时接触时的平均运算时间散点。运算时间的测试结果显示，1个接触的运算时间约为0.9 ms，11个接触点时平均每个接触点的运算时间为1.9 ms。随着接触点数目的增长，运算时间也急剧上升。到44个接触点时，平均每个接触的运算时间为8 ms，实时性有所降低。

图13 不同接触点数目的运算时间Fig.13 Computation time of various numbers of contact points

3 讨论

单接触面接触的实验结果表明，单接触面接触过程的运算速度快，计算出来的接触力稳定，能很好地达到虚拟手术中的实时性和真实性要求。

双接触面接触的实验结果表明，在双接触面的情况下能计算出接触力和摩擦力，稳定性能得到保证，得到良好的摩擦效果，仍能满足虚拟手术中的实时性和真实性要求。

在多接触点同时接触的情形下，碰撞力的计算时间与接触点的数量有直接的关系，会随着接触点的增加而增长。在多模型的场景中，碰撞的发生会较多，从而碰撞点增多，碰撞的时间也越长。但是，若模型不发生碰撞，那么模型数量影响的仅仅是下一个时间步粒子群位置的计算时间，接触计算不会发生。

在本次多接触点接触实验中，采用类高斯-赛德尔算法求解线性互补问题。类高斯-赛德尔算法的时间复杂度为O(m2)，m为接触点的数目，显然求解的时间随着接触点的增加而增加，与本次实验的结果基本吻合。这表明，随着接触点的增加，仍然会使总的求解时间上升，从而影响实时性。国外有一种求解带摩擦接触的k边形金字塔算法[7]，其时间复杂度为O(k×m2)，m为接触点的数目，k为近似摩擦圆锥底面的多边形边数。显然，类高斯-赛德尔算法更具优势。

本研究表明，手术交互的方式复杂而多样，即便是摩擦接触交互，对不同的接触方式也不能使用单一不变的方法。本研究使用了基于约束的方法来处理摩擦接触交互问题：对于单接触面接触，使用单面约束和摩擦约束，简单而有效；对于双接触面接触，则原本的约束方程不再适用，改为基于三角面片顶点及其邻接顶点的单面约束和摩擦约束。双接触面接触的情形比单接触面的情况要复杂一些，但能提高交互的稳定性，也没有跳出约束方程这一方法的框架。由此表明，使用基于约束方程的方法，能够应对复杂多变的交互，即使交互的形式不同，仍能通过适当地改变约束方程来实现。

4 结论

本研究提出了虚拟手术中带摩擦的接触力的计算方法。首先，建立计算机图形学物理仿真中最为常用的粒子系统，通过碰撞检测出潜在的碰撞对象，精确计算碰撞接触点和接触力的方向，总体的常微分方程通过隐式欧拉方法进行时间积分；其次，建立碰撞接触和摩擦接触的约束方程，并通过类高斯-赛德尔方法求解；再次，分别进行10次单接触面和双接触面接触碰撞实验，碰撞时分析碰撞力的大小和方向、响应时间、稳定时候的起伏是否符合要求；最后，再细分半月板模型的网格，使得碰撞时的接触点增加，进行100次实验，计算1～45个接触点同时接触的平均运算时间。

实验结果表明，通过这一系列方法，在建立的变形体夹取场景中，能绘制出符合实际情况的力的大小和方向。该方法在单接触面和双接触面接触的情况下，满足虚拟手术的真实性和实时性要求。在多接触点接触的情况下，随着接触点数目的增长，运算时间也急剧上升，然而在11个接触点的情况下，仍能满足实时性要求。求解约束方程所用的高斯-赛德尔方法具有更低的时间复杂度O(m2)，而且有很好的灵活性，通过稍加修改即可用于双接触面接触的情形，并且能很好地保持其性能。不过，在接触点的数目上升到一定程度后，还是会对仿真的实时性产生影响。

未来将会与自行设计的力反馈设备进行虚耦合，建立一整套力反馈交互系统。