李晉芳，韋光揚(yáng)，何漢武，蔡嘉鴻，陳基榮

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

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

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

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

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

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

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

1 牙齦軟組織的物理建模

1.1 質(zhì)點(diǎn)彈簧模型

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

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

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

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

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

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

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

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

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

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

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

化簡(jiǎn)式(2)得位移表達(dá)式，如式(3)所示：

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

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

2 中點(diǎn)平均算法

2.1 算法基礎(chǔ)

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

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

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

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

圖 3 中點(diǎn)平均算法流程圖Fig.3 Midpoint-averaging algorithm step flow chart

2.2 算法后處理

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

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

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

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

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

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

因此，在原有算法基礎(chǔ)上，在確定的步長(zhǎng)下，可以通過(guò)求解一個(gè)適合的α來(lái)提升求解精度，并且不增加計(jì)算量。

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

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

3 算法與歐拉算法的比較

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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