李晉芳，韋光揚，何漢武，蔡嘉鴻，陳基榮

（1. 廣東工業(yè)大學 機電工程學院，廣東 廣州 510006；2. 廣東工貿(mào)職業(yè)技術(shù)學院 機電工程學院，廣東 廣州 510510）

虛擬現(xiàn)實技術(shù)(Virtual Reality，VR)的應(yīng)用集中于虛擬手術(shù)、遠程醫(yī)療系統(tǒng)、醫(yī)學教育等方面。其具有許多優(yōu)點，可以提高訓練效率，有效降低手術(shù)訓練成本，減少一些當代教學實踐帶來的負面影響[1-2]。

在虛擬手術(shù)的眾多研究中，軟組織形變建模是其核心。常見的軟組織有皮膚、肌肉、神經(jīng)、血管、黏膜等。雖然它們在結(jié)構(gòu)上差異明顯，但是力學特點相似，都涉及復(fù)雜生物力學[3]，在力的作用下會發(fā)生形變，表現(xiàn)出明顯的非線性。

牙齦是指緊貼于牙頸周圍及鄰近的牙槽骨上淡紅色的結(jié)構(gòu)，由復(fù)層扁平上皮及固有層組成，是一種典型的軟組織。牙齦軟組織血管豐富，呈淡紅色，堅韌而有彈性，直接與骨膜緊密相連，是口腔黏膜的一部分，一般根據(jù)牙齦的厚度可分為薄型和厚型兩種。虛擬口腔手術(shù)中牙齦軟組織的建模仿真研究主要抓住牙齦的這幾個特點制定相應(yīng)的技術(shù)方案。

有限元法和質(zhì)點彈簧模型是當前軟組織建模的兩種重要的建模方法[4-5]。有限元法在擁有較少的材料參數(shù)時也可以描述某一個物理系統(tǒng)，然而運算量較大，實現(xiàn)起來較困難，耗時較長[6-7]，尤其是對復(fù)雜的場景，所耗費的時間和計算機內(nèi)存是相當大的，因此它不適用于虛擬切割這類實時性要求高的仿真。

相比有限元法，質(zhì)點彈簧模型由于無須連續(xù)參數(shù)化，既可以用于靜態(tài)分析，也可以用于動態(tài)分析，建模簡單、計算速度較快、去除和增加點的操作較容易實現(xiàn)[8]。另外，如果整個模型頂點都被計算會造成計算量巨大和不符合實際，離碰撞點遠的軟組織區(qū)域變化不大，尤其是牙齦軟組織。

因此本文提出一種基于質(zhì)點彈簧模型的適用于牙齦軟組織的算法，算法只對碰撞點附近的網(wǎng)格頂點進行更新計算，能夠在保證軟組織形變符合要求的前提下減少計算量，提升求解精度和系統(tǒng)穩(wěn)定性。此外，本文通過仿真模擬驗證所提出算法的處理效果。

1 牙齦軟組織的物理建模

1.1 質(zhì)點彈簧模型

質(zhì)點彈簧模型是一種常用的軟組織物理建模方法。主要思想是把仿真對象用質(zhì)點離散化，質(zhì)點之間用符合線性彈性模型的彈簧連接而成，質(zhì)點除受彈簧的彈力作用外，彈簧產(chǎn)生的阻尼力與速度成比例關(guān)系[9]。

質(zhì)點彈簧模型如圖1所示，網(wǎng)格的邊界和節(jié)點由質(zhì)點和彈簧組成，一個彈簧連接兩個質(zhì)點，而一個質(zhì)點可以連接多個彈簧，這些彈簧是線性的，遵循胡克定律[10-11]。質(zhì)點在受到外力時，產(chǎn)生的外力作用傳遞到相鄰質(zhì)點，相鄰質(zhì)點接著傳下去帶動其相鄰質(zhì)點運動，變形由質(zhì)點運動產(chǎn)生。軟組織形變仿真需要利用數(shù)值解算的方法求解形變的動力學模型[12]。

圖 1 質(zhì)點彈簧模型Fig.1 Mass-spring model

質(zhì)點彈簧模型運算效率高且易于編程實現(xiàn)，能很好地滿足虛擬口腔手術(shù)的實時性要求[13]，展示牙齦軟組織的粘彈性、塑性和非線性等特性的存在。

1.2 一般的數(shù)值求解方法及對比

在基質(zhì)點彈簧模型的虛擬口腔手術(shù)系統(tǒng)中，形變仿真的時間步長直接決定了虛擬手術(shù)的實時性[14]，選取合適的時間步長是保證牙齦軟組織形變精度和實時性的必要前提。求解質(zhì)點彈簧模型常用的數(shù)值求解方法主要有顯式歐拉算法、隱式歐拉算法、Runge-Kutta算法等。

為驗證各算法的精度分析，本文通過MATLAB以步長h=0.05，分別使用顯式歐拉算法、隱式歐拉算法、Runge-Kutta算法求解式(1)的常微分方程，給定自變量x，求解值y，同時計算理論值，綜合繪制出曲線圖如圖2所示。

圖 2 3種算法精度比較圖Fig.2 Comparisons of three algorithms in precision

從圖2中可以看出，Runge-Kutta算法求解結(jié)果幾乎與理論曲線重合，它的精度最高，而顯式歐拉算法精度相比Runge-Kutta算法和隱式歐拉算法最低，隱式歐拉算法居中。

根據(jù)以上的分析，可知Runge-Kutta算法、隱式歐拉算法都具有較好的穩(wěn)定性，精度較顯式歐拉算法高。但是它們共同的缺點是計算量較大，不適合實時性要求高的虛擬手術(shù)系統(tǒng)，所以大部分的虛擬手術(shù)都采用顯式歐拉算法來求解質(zhì)點彈簧模型，從而進行軟組織形變的模擬。顯式歐拉算法一般只是一階收斂，精度不高，誤差較大，為實現(xiàn)逼真的模擬效果，時間步長需要設(shè)置得很小，這會導致整個形變過程延長，并存在進退性沖擊波[15-16]。

為了提高顯式歐拉算法的計算精度，文獻[16]提出了一種改進的歐拉算法來求解質(zhì)點彈簧系統(tǒng)的速度與位移矢量，以提高求解精度及降低求解步長的影響。運用顯式歐拉算法求解速度v，并直接使用v的結(jié)果，用隱式的方法求解位移x。表達式如式(2)所示。其中，m為質(zhì)點質(zhì)量，變量上標t表示本時刻該變量的值，t+h表示下一時刻該變量的值，f 為系統(tǒng)所受的力。

化簡式(2)得位移表達式，如式(3)所示：

式(3)表明速度會影響系統(tǒng)所受的力f，從而影響位移x的求解精度。

這種算法相較于顯式歐拉算法提高了計算精度和穩(wěn)定性，但由于速度v采用顯式歐拉算法求解，所以求解速度的精度還不夠高，且受步長h的影響還比較大。

2 中點平均算法

2.1 算法基礎(chǔ)

中點法常用于求解常微分方程，即應(yīng)用于求解質(zhì)點彈簧模型。采用本時刻 ft與下一時刻 ft+h之間的中間時刻t +h/2 處的 ft+h/2來求解速度vt+h。它所需的時間步長小于Runge-Kutta算法，精度比顯式歐拉算法要高，但低于Runge-Kutta算法，具有較快的計算速度和較好的實時性。

平均法采用 vt和vt+h的平均值來求解 xt+h，表達式如式(4)所示。它的計算速度相比顯式歐拉算法慢，但比Runge-Kutta算法要快，具有較好的計算精度和穩(wěn)定性。

將式(4)中的vt+h代入位移x的表達式中，化簡位移表達式，如式(5)所示：

基于這兩種算法，同時考慮到形變的實時性和精度，本文提出一種改進的算法——中點平均算法，即算法融合中點法和平均法的特點，采用中點法對速度v進行迭代求解，位移x則采用平均法進行迭代求解。雖然減緩了計算速度，但能夠綜合兩者的優(yōu)勢，減小誤差及減小受步長的影響，提升求解的精度。算法流程如圖3所示。

圖 3 中點平均算法流程圖Fig.3 Midpoint-averaging algorithm step flow chart

2.2 算法后處理

在算法基礎(chǔ)上，對提出的中點平均算法進行后處理，并提出一個比重參數(shù)α。大致思想是：在采用中點平均算法求解質(zhì)點彈簧模型前，對于求解速度v，以步長h=0.000 1的Runge-Kutta算法的求解結(jié)果作為基準，不斷對式(6)所示的速度表達式進行求解并與基準比較，直到找到一個合適的α，最后用此中點平均算法求解虛擬口腔手術(shù)系統(tǒng)中的質(zhì)點彈簧模型，方法流程如圖4所示。

式(6)中α分別取值1/2，3/8，5/8，來驗證α值的大小對求解精度有著不同程度的影響。

圖 4 提升求解精度的方法流程圖Fig.4 Method sketch for improving solution precision

以步長h=0.000 1的Runge-Kutta算法求解的速度為理論曲線。在算法的對比中，步長大小的選取會影響比較結(jié)果的直觀程度。因此這里選取步長較大的2個值，h為0.1和0.4來求解質(zhì)點運動狀態(tài)。分別繪制出質(zhì)點速度v的理論曲線和α=1/2，α=3/8，α=5/8時的求解速度v曲線，如圖5的(a)、(b)所示。若某α下求解的速度v曲線與理論曲線更接近，則證明該α值調(diào)整的中點平均算法具有更高的求解精度。

圖 5 不同α對應(yīng)的速度曲線Fig.5 Speed curves corresponding to different α

從圖5結(jié)果中可看出，步長h=0.1，α=1/2時取得最好的求解精度；步長h=0.4，α=3/8時取得最好的求解精度。

因此，在原有算法基礎(chǔ)上，在確定的步長下，可以通過求解一個適合的α來提升求解精度，并且不增加計算量。

在Unity3D中建立基于質(zhì)點彈簧模型的虛擬牙齦軟組織模型，并應(yīng)用后處理的中點平均算法進行求解，圖6為牙齦軟組織切割前后的形變效果圖。

圖 6 牙齦軟組織切割效果Fig.6 Cutting effect of gingival tissue

3 算法與歐拉算法的比較

中點平均算法與改進歐拉算法可通過式（5）和式（3）進行比較。改進歐拉算法采用t時刻的力ft計算下一時刻的位移 xt+h，而中點平均算法采用t+h/2 時刻的力 ft+h/2來 求解下一時刻的位移 xt+h，可知中點平均算法具有更高的計算精度，并且中點平均算法沒有用到顯式歐拉算法，使它對求解步長的影響不大。所以對于計算精度與受步長的影響，中點平均算法都具有更好的表現(xiàn)。

在MATLAB中，以步長h為0.000 1的Runge-Kutta算法求解的位移和作用力曲線作為理論值曲線，分別采用顯式歐拉算法、改進歐拉算法及中點平均算法對質(zhì)點進行狀態(tài)求解，作出求解結(jié)果曲線。從而驗證中點平均算法的求解精度優(yōu)于顯式歐拉算法和改進歐拉算法，對比如圖7所示。

圖 7 不同數(shù)值求解方法求解質(zhì)點運動狀態(tài)對比圖Fig.7 Comparisons of particle motion state by different numerical solutions

從圖7可知，中點平均算法具有較好的求解精度，幾乎與理論線重合。改進歐拉算法相較顯式歐拉算法有更高的求解精度，但沒有中點平均算法高。

這里進一步探討步長的選取分別對顯式歐拉算法、改進歐拉算法以及中點平均算法的影響。在MATLAB中分別用這3種算法對不同步長h的質(zhì)點運動狀態(tài)進行求解。

在圖8中繪制出顯式歐拉算法、改進歐拉算法及中點平均算法求解的質(zhì)點運動狀態(tài)曲線圖。

圖 8 3種算法求解的質(zhì)點運動狀態(tài)Fig.8 Motion state of point solved by three algorithms

對比圖8的(a)、(b)、(c)可以發(fā)現(xiàn)，顯式歐拉算法受步長h的影響最大。改進歐拉算法受步長影響相對于顯式歐拉算法更小，但由于其求解速度的時候采用的是顯式歐拉算法，所以它受步長h影響依然較大。而中點平均算法采用中點法求解速度，然后再用平均法求解位移，這樣就一定程度上減小了受步長的影響。從圖8(c)可以看出，當h取0.02或h取0.1時，所求的結(jié)果幾乎與理論結(jié)果重合；當步長h取0.4這樣較大的值時，中點平均算法相比顯式歐拉算法和改進歐拉算法依舊有更高的求解精度，這充分說明了中點平均算法能很好地滿足虛擬口腔手術(shù)的實時性和真實性要求。

這里對虛擬口腔手術(shù)器械碰撞牙齦軟組織進行模擬，在Unity3D中分別應(yīng)用中點平均算法、顯式歐拉算法和改進歐拉算法建立的質(zhì)點彈簧模型模擬牙齦軟組織的形變，圖9為牙齦軟組織碰撞前后的形變效果圖。

圖 9 牙齦軟組織碰撞效果Fig.9 Collision effect of gingival tissue

從圖9可以看到，在相同的環(huán)境和作用力下，3種算法處理的軟組織形變都滿足虛擬口腔手術(shù)的真實性要求。中點平均算法處理的牙齦形變效果相較另外兩種算法要更平滑，這也說明了中點平均算法相比具有更高的求解精度。

4 結(jié)論

本文提出了一種中點平均算法并對其進行后處理，使其對基于質(zhì)點彈簧模型的牙齦軟組織的求解精度在不增加計算量的前提下進一步提升。隨后通過在MATLAB中分別采用顯式歐拉算法、改進歐拉算法和中點平均算法這3種算法求解質(zhì)點彈簧模型中質(zhì)點的運動狀態(tài)，驗證了中點平均算法在三者中具有更高的求解精度，以及受時間步長的影響更小，滿足虛擬口腔手術(shù)的實時性要求。在Unity3D中分別應(yīng)用這3種算法對牙齦軟組織模擬仿真，中點平均算法的形變效果也更為理想。盡管本文提出并進行處理的中點平均算法已展示其優(yōu)越性，但算法尚需繼續(xù)優(yōu)化，使之得到更廣泛的應(yīng)用。