歐惠棠 李晉芳 莫建清

(廣東工業(yè)大學機電工程學院 廣州 510006)

引言

虛擬醫(yī)療是應用虛擬現(xiàn)實(virutalreatility,VR)技術,構建虛擬的人體模型器官以及手術，提高虛擬實境的真實感。借助于虛擬外設，可以使人們更逼真地學習醫(yī)療知識以及治病救人的虛擬現(xiàn)實應用，在未來有著非常樂觀的前景。其中，手術訓練與虛擬現(xiàn)實相結合的虛擬手術是現(xiàn)代手術訓練的新趨勢。虛擬手術仿真利用醫(yī)學生物信息和計算機仿真技術以及其他仿真模擬器械，為醫(yī)療工作者提供一個虛擬且逼真的手術場景，使其借助這些虛擬場景進行手術計劃訓練和研究，積累經(jīng)驗，掌握技巧，以便在實際的手術過程中完成復雜的手術操作，為病患減少手術后的痛苦和手術風險[1]。通過這一手段，無疑可使訓練更真實、準確、可靠。

國內(nèi)外對于虛擬手術仿真中碰撞力的研究，著重于變形體形變技術的實現(xiàn)上，如Montgomery等提出的現(xiàn)在非常流行的彈簧-質(zhì)點網(wǎng)絡(spring-mass networks)[2]，F(xiàn)elippa等提出的基于共旋轉(zhuǎn)模型的有限元方法(FEM)[3]，吳涓等的一種基于沿徑向方向分割為呈同心圓分布的彈簧-質(zhì)點模型及實時力學響應算法[4]，陳衛(wèi)東等提出的虛擬體彈簧的變形算法[5]等。這些研究僅僅著重于變形體內(nèi)部力的計算上，并沒有系統(tǒng)地給出碰撞檢測、碰撞力的計算、時間積分等整個仿真流程。國內(nèi)的研究僅僅給出當模型發(fā)生形變時內(nèi)力與粒子群位移之間的關系，卻沒有給出與手術器械模型發(fā)生碰撞時如何計算出碰撞力并使變形體發(fā)生變形。國外Hasegawa等提出懲罰方法來求解碰撞接觸力[6]，但會導致一系列不穩(wěn)定性問題。Anitescu等提出k邊形金字塔算法[7]，但時間復雜度較大，導致運算時間過長。

本研究提出虛擬手術訓練中帶摩擦的碰撞力計算及其仿真方法，闡述粒子系統(tǒng)的建立、碰撞檢測得到碰撞力方向、隱式歐拉積分、碰撞約束方程(線性互補問題)、摩擦約束方程的構建、求解方程得出碰撞力大小這一系列的完整流程。該方法適用于剛體和變形體(彈簧質(zhì)點網(wǎng)絡和有限元)物理模型仿真，不僅實現(xiàn)了碰撞功能，而且還添加了摩擦特性，能夠在一定接觸點數(shù)目的場景下滿足交互的真實性、穩(wěn)定性和實時性的要求。

在計算機圖形學中，物理仿真通過粒子系統(tǒng)[8]實現(xiàn)。粒子系統(tǒng)可以建立剛體和變形體等物理對象，并能實現(xiàn)其物理行為。本文所述的碰撞研究對剛體和變形體均適用。對于碰撞物理仿真的流程，已有相關的研究。下面通過約束動力學構建碰撞與摩擦相關約束方程并求解，然后通過實驗驗證是否計算正確和性能是否符合要求。

1 材料和方法

1.1 材料

采用自行設計的虛擬手術仿真系統(tǒng)，加入帶摩擦的碰撞力計算算法，進行藍鉗模型和半月板模型的碰撞測試。實驗環(huán)境為Win10 64 位系統(tǒng)，處理器為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

為了計算接觸力，必須先通過碰撞檢測，檢測出潛在的碰撞對象，當模型發(fā)生碰撞后，精確計算出碰撞點的位置和接觸力的方向(大小未知)，然后構建相應的約束方程，再利用相應的數(shù)值方法進行求解。求解后，通過圖形界面驗證計算出的反饋力的方向和大小是否符合要求，通過20次的碰撞測試，并作出碰撞力的曲線，驗證接觸力的響應時間和穩(wěn)定性是否符合要求。通過細分半月板模型的網(wǎng)格并通過100次的實驗，驗證單接觸面接觸從1個接觸點到50個接觸點同時接觸的平均碰撞時間。

1.2.1碰撞檢測

如圖2所示，通過碰撞檢測，可以得出兩個即將碰撞的模型的潛在碰撞點P和Q?？梢匀我膺x取法向量，本研究選取模型D1上P點對應的朝內(nèi)的法向量，記為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插值

在計算機圖形系統(tǒng)中，模型表面為三角面片構成的網(wǎng)格。模型之間的碰撞化歸為兩個三角形之間的標準碰撞，如圖3所示。碰撞檢測出來的碰撞點可能位于三角形的內(nèi)部或邊線上，但由于在粒子系統(tǒng)中，粒子通常是三角面片的頂點，因此必須把碰撞點的碰撞力等效到三角形的頂點上。

圖3 兩個三角形碰撞的全部5種情況Fig.3 All five cases of the collision between two triangles

根據(jù)解析幾何相關知識，三角形所在平面的任意點都能表示為頂點的加權平均值，其中的權稱為重心坐標。在圖4中，對于點P，有

P=αA1+βB1+γC1

(4)

式中，α、β、γ為重心坐標分量。

式(4)中兩邊取微位移，有

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

(5)

圖4 碰撞點在三角形內(nèi)部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)，不詳細討論。接觸蒙皮能使得兩個模型之間三角面片頂點的約束轉(zhuǎn)化為模型之間粒子群的約束。因此，δn(P)≥0可表示為

Ψ(q1,q2)≥0

(9)

稱為單面約束，Ψ僅與位置有關，稱為完整約束[11]。

粒子系統(tǒng)的物理仿真實質(zhì)上是求解牛頓第二定律常微分方程(ODE)有

(10)

式中，fext為外部力，fint為內(nèi)部力，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，求得的解為沒有約束的運動狀態(tài)，稱為自由運動解，為dvfree1和dvfree2。

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

(16)

最后得到

(17)

對于線性互補問題，其求解過程可由圖5直觀表示。

1.2.4建立摩擦約束

圖5 線性互補問題求解。(a)未發(fā)生碰撞；(b)發(fā)生碰撞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)。對于復雜的手術交互，除了單接觸面碰撞之外，還會涉及雙接觸面碰撞，如軟組織夾取交互。若雙接觸面碰撞仍按照單接觸面碰撞的情況處理，則一個接觸面碰撞響應會影響另一個接觸面的碰撞響應，導致兩個接觸面不能穩(wěn)定地產(chǎn)生碰撞壓力。靈活地改變約束方程，能使約束方程仍然適用于雙接觸面的情形。

如果碰撞檢測出軟組織的同時與兩個接觸面發(fā)生碰撞，則約束方程發(fā)生改變。如圖7所示，圖中三角面片代表接觸區(qū)域。(a)中軟組織模型與夾子模型的一個面接觸，此時采用單接觸面約束求解方法。(b)中軟組織模型與夾子模型的兩個面同時接觸，則原來的單接觸面狀態(tài)轉(zhuǎn)變?yōu)殡p接觸面狀態(tài)，使用雙接觸面約束求解方法。

圖7 接觸狀態(tài)的轉(zhuǎn)變。(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為頂點及其鄰接頂點的數(shù)目)。

由于新的約束方程允許fn≤0，并且把三角面片頂點及其鄰接頂點法向力的綜合作用考慮進來，可以有效提高交互的穩(wěn)定性，并且實質(zhì)上更符合于現(xiàn)實的接觸情況。

從而雙接觸面接觸總的約束方程為

(25)

式中，m為接觸的三角面片頂點的數(shù)目。

1.2.6求解總約束方程

對于式(25)，使用類高斯-賽德爾算法(Gauss-Seidel-like algorithm)求解。在粒子仿真計算力學領域[13]中，該算法能保證較好的收斂性。先對接觸逐個求解，然后不停地迭代，直到數(shù)值收斂。對于m個瞬時摩擦接觸，先考慮其中一個摩擦接觸α，對其進行求解，有

(26)

求出fα，其他接觸也以此進行計算，然后不斷迭代求出新的fα，直到數(shù)值收斂。

(27)

本研究僅給出雙接觸面接觸約束求解算法，單接觸面接觸求解算法與雙接觸面接觸類似，并且要簡單，由于篇幅所限，這是不再贅述。由于雙接觸面接觸要求查詢?nèi)敲嫫旤c的鄰接頂點，故模型的網(wǎng)格要求為流形網(wǎng)格，并采用半邊數(shù)據(jù)結構，可以使查詢操作的時間復雜度為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，摩擦系數(shù)μ取0.5。半月板模型設定楊氏模量為59 MPa，泊松比為0.45[14]。實驗方案為先分別進行10次單接觸面和雙接觸面接觸碰撞實驗。單接觸面接觸時，藍鉗模型簡單地碰撞半月板模型；雙接觸面接觸時，藍鉗模型夾取碰撞模型。碰撞時，分析碰撞力的大小和方向、響應時間、穩(wěn)定時候的起伏是否符合要求。然后，再細分半月板模型的網(wǎng)格，使得碰撞時的接觸點增加，進行100次實驗，計算1個接觸點到45個接觸點同時接觸的平均運算時間，其運算時間的記錄方法如圖8所示。

圖8 碰撞時間測試流程Fig.8 Diagram of the collision time test

1.2.8統(tǒng)計學分析

單接觸面接觸中，以膝關節(jié)模型的正方向(髕骨方向)為基準，分別在豎直和水平方向從-90°、45°、0°、45°、90°各進行1次碰撞實驗，計算平均的響應時間和穩(wěn)定時的起伏。在雙接觸面接觸中，從水平方向-90°、45°、0°、45°、90°進行10次夾取試驗，計算平均的響應時間和穩(wěn)定時的起伏。在多接觸點平均碰撞力運算時間的測試中，為了能使足夠的接觸點接觸，接觸的方向任意選取，然后計算平均計算時間。

2 結果

圖9為筆者構建的膝關節(jié)鏡手術場景，在場景中，藍鉗模型與半月板模型發(fā)生碰撞。碰撞接觸力通過上述方法得到計算并繪制，由圖可知接觸力的方向符合實際情況。

圖9 場景中碰撞力的繪制Fig.9 Contact force rendering in the scene

圖10為10次單接觸面碰撞實驗中其中一次的碰撞力的曲線。由圖可知，模型在1和2 s的時候發(fā)生了碰撞，x、y、z方向力的大小變化趨勢無顯著差異，開始碰撞時，接觸力急劇上升，隨著半月板模型的變形行為趨于穩(wěn)定，碰撞力也趨于穩(wěn)定。在2 s的時候，藍鉗模型在z軸負方向前進0.2個單位長度，穩(wěn)定后碰撞力增加。在10次的實驗中，模型碰撞的平均響應時間為0.02 s，接觸力穩(wěn)定后的平均起伏范圍為±0.02個單位。

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

圖11 通過摩擦力夾取變形體Fig.11 Grasping a deformable object by friction

圖11為半月板模型夾取的效果。在不同的時間段，其狀態(tài)從左到右、從上到下分別為：準備夾取階段、夾取階段、向右拖動階段、向下拖動階段、釋放階段。

圖12為10次夾取實驗中的一次藍鉗下顎碰撞力曲線。在6.12 s的時候，藍鉗夾取半月板，接觸力在y方向發(fā)生突變，隨后趨于穩(wěn)定；在8.78 s的時候，藍鉗拉動半月板運動，摩擦使接觸力在z軸方向上發(fā)生突變；在10.78 s的時候，藍鉗釋放半月板，接觸力回復為零。在10次的實驗中，模型碰撞的平均響應時間為0.02 s，接觸力穩(wěn)定后的平均起伏范圍為±0.02個單位。

圖13為100次單接觸面碰撞下1～45個接觸點同時接觸時的平均運算時間散點。運算時間的測試結果顯示，1個接觸的運算時間約為0.9 ms，11個接觸點時平均每個接觸點的運算時間為1.9 ms。隨著接觸點數(shù)目的增長，運算時間也急劇上升。到44個接觸點時，平均每個接觸的運算時間為8 ms，實時性有所降低。

圖13 不同接觸點數(shù)目的運算時間Fig.13 Computation time of various numbers of contact points

3 討論

單接觸面接觸的實驗結果表明，單接觸面接觸過程的運算速度快，計算出來的接觸力穩(wěn)定，能很好地達到虛擬手術中的實時性和真實性要求。

雙接觸面接觸的實驗結果表明，在雙接觸面的情況下能計算出接觸力和摩擦力，穩(wěn)定性能得到保證，得到良好的摩擦效果，仍能滿足虛擬手術中的實時性和真實性要求。

在多接觸點同時接觸的情形下，碰撞力的計算時間與接觸點的數(shù)量有直接的關系，會隨著接觸點的增加而增長。在多模型的場景中，碰撞的發(fā)生會較多，從而碰撞點增多，碰撞的時間也越長。但是，若模型不發(fā)生碰撞，那么模型數(shù)量影響的僅僅是下一個時間步粒子群位置的計算時間，接觸計算不會發(fā)生。

在本次多接觸點接觸實驗中，采用類高斯-賽德爾算法求解線性互補問題。類高斯-賽德爾算法的時間復雜度為O(m2)，m為接觸點的數(shù)目，顯然求解的時間隨著接觸點的增加而增加，與本次實驗的結果基本吻合。這表明，隨著接觸點的增加，仍然會使總的求解時間上升，從而影響實時性。國外有一種求解帶摩擦接觸的k邊形金字塔算法[7]，其時間復雜度為O(k×m2)，m為接觸點的數(shù)目，k為近似摩擦圓錐底面的多邊形邊數(shù)。顯然，類高斯-賽德爾算法更具優(yōu)勢。

本研究表明，手術交互的方式復雜而多樣，即便是摩擦接觸交互，對不同的接觸方式也不能使用單一不變的方法。本研究使用了基于約束的方法來處理摩擦接觸交互問題：對于單接觸面接觸，使用單面約束和摩擦約束，簡單而有效；對于雙接觸面接觸，則原本的約束方程不再適用，改為基于三角面片頂點及其鄰接頂點的單面約束和摩擦約束。雙接觸面接觸的情形比單接觸面的情況要復雜一些，但能提高交互的穩(wěn)定性，也沒有跳出約束方程這一方法的框架。由此表明，使用基于約束方程的方法，能夠應對復雜多變的交互，即使交互的形式不同，仍能通過適當?shù)馗淖兗s束方程來實現(xiàn)。

4 結論

本研究提出了虛擬手術中帶摩擦的接觸力的計算方法。首先，建立計算機圖形學物理仿真中最為常用的粒子系統(tǒng)，通過碰撞檢測出潛在的碰撞對象，精確計算碰撞接觸點和接觸力的方向，總體的常微分方程通過隱式歐拉方法進行時間積分；其次，建立碰撞接觸和摩擦接觸的約束方程，并通過類高斯-賽德爾方法求解；再次，分別進行10次單接觸面和雙接觸面接觸碰撞實驗，碰撞時分析碰撞力的大小和方向、響應時間、穩(wěn)定時候的起伏是否符合要求；最后，再細分半月板模型的網(wǎng)格，使得碰撞時的接觸點增加，進行100次實驗，計算1～45個接觸點同時接觸的平均運算時間。

實驗結果表明，通過這一系列方法，在建立的變形體夾取場景中，能繪制出符合實際情況的力的大小和方向。該方法在單接觸面和雙接觸面接觸的情況下，滿足虛擬手術的真實性和實時性要求。在多接觸點接觸的情況下，隨著接觸點數(shù)目的增長，運算時間也急劇上升，然而在11個接觸點的情況下，仍能滿足實時性要求。求解約束方程所用的高斯-賽德爾方法具有更低的時間復雜度O(m2)，而且有很好的靈活性，通過稍加修改即可用于雙接觸面接觸的情形，并且能很好地保持其性能。不過，在接觸點的數(shù)目上升到一定程度后，還是會對仿真的實時性產(chǎn)生影響。

未來將會與自行設計的力反饋設備進行虛耦合，建立一整套力反饋交互系統(tǒng)。