萬(wàn) 燕, 朱甜甜, 姚 礪, 馮子涵

(東華大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，上海 201620)

隨著數(shù)字化影業(yè)的蓬勃發(fā)展，3D游戲、虛擬試衣和動(dòng)畫(huà)電影迅速發(fā)展，成為了熱門(mén)的新興產(chǎn)業(yè)。動(dòng)畫(huà)中虛擬角色建模是一個(gè)重要的環(huán)節(jié)，其中，服飾層的建模工作量占80%以上，要使虛擬角色具有令人滿(mǎn)意的真實(shí)感，布料動(dòng)畫(huà)是必不可少的。因此，有關(guān)布料模擬一直是研究熱點(diǎn)。

布料動(dòng)畫(huà)[1]主要是在計(jì)算機(jī)中建立布料的虛擬模型，產(chǎn)生逼真的形變，實(shí)現(xiàn)高度真實(shí)感和實(shí)時(shí)的動(dòng)態(tài)效果。早期的布料建模主要采用幾何方法，但該方法通常只產(chǎn)生布料的靜態(tài)變形狀態(tài)，很難實(shí)現(xiàn)真實(shí)連貫的復(fù)雜動(dòng)畫(huà)。隨著計(jì)算機(jī)圖形學(xué)技術(shù)的發(fā)展，采用物理方法[2]來(lái)模擬布料受到越來(lái)越多專(zhuān)家和學(xué)者的青睞。物理方法將布料視為連續(xù)的彈性介質(zhì)(如有限元模型[3])，或者由相互作用的非連續(xù)粒子組成(如質(zhì)點(diǎn)-彈簧模型[4])。盡管不同的物理建模方法對(duì)于布料的表達(dá)方式不盡相同，但是大部分是根據(jù)牛頓運(yùn)動(dòng)定律計(jì)算布料的基本形態(tài)?；拘螒B(tài)在數(shù)學(xué)上的表現(xiàn)形式為動(dòng)力學(xué)方程，針對(duì)不同的動(dòng)力學(xué)方程進(jìn)行求解，需要對(duì)模型進(jìn)行改進(jìn)。

其中，質(zhì)點(diǎn)-彈簧模型為表現(xiàn)柔性材料動(dòng)態(tài)特性的建模提供了一種簡(jiǎn)單而實(shí)用的方法，這些柔性材料包括布料[5]、頭發(fā)[6]和可變形固體[7]。然而，與其他用于彈性物體的建模方法一樣，獲得真實(shí)的柔性材料行為通常需要本構(gòu)參數(shù)，這些參數(shù)與剛性系統(tǒng)的數(shù)值計(jì)算方法有關(guān)。數(shù)值計(jì)算方法主要有顯式時(shí)間積分和隱式時(shí)間積分。顯式時(shí)間積分方法具有快速性，但是當(dāng)應(yīng)用于剛性系統(tǒng)時(shí)，其會(huì)出現(xiàn)不穩(wěn)定的問(wèn)題。Baraff等[8]使用隱式時(shí)間積分方法進(jìn)行數(shù)值計(jì)算，該方法能夠保證穩(wěn)定性，但是需要求解大的非線(xiàn)性方程組。由于大的非線(xiàn)性方程組的計(jì)算量大，這限制了隱式時(shí)間積分方法在實(shí)時(shí)應(yīng)用中的使用。另外，就模擬精度而言，隱式積分引入過(guò)多的近似，由此造成的過(guò)度數(shù)值阻尼問(wèn)題會(huì)嚴(yán)重影響模擬質(zhì)量。

牛頓法以其出色的收斂特性而聞名。當(dāng)?shù)銐蚪咏顑?yōu)解時(shí)，牛頓方法表現(xiàn)出二次收斂性，這優(yōu)于塊坐標(biāo)下降法的線(xiàn)性收斂性。但是，牛頓法的每次迭代都需要重新求解系統(tǒng)海森矩陣的逆矩陣，這增加了系統(tǒng)的計(jì)算成本。因此，在系統(tǒng)有限的計(jì)算成本下，根據(jù)獲得穩(wěn)定、可視化的布料效果實(shí)際需求，只使用牛頓法一次迭代下的解，這是未完全收斂下的數(shù)值解，距離實(shí)際的精確解還存在較大誤差。

Müller等[9]提出的位置動(dòng)力學(xué)(position based dynamics， PBD)法從約束投影的角度來(lái)模擬各種物體的材料屬性。由于其健壯性、簡(jiǎn)單性和快速性，PBD在游戲和動(dòng)畫(huà)行業(yè)中得到廣泛應(yīng)用。PBD與傳統(tǒng)的彈性模型和本文的彈簧投影不同，其使用了與標(biāo)準(zhǔn)胡克模型不兼容的參數(shù)，因此，模擬材料的剛度會(huì)受到迭代次數(shù)的影響。

針對(duì)以上問(wèn)題，本文提出了一種應(yīng)用于質(zhì)點(diǎn)-彈簧系統(tǒng)的隱式快速解算器，該解算器重新優(yōu)化了Martin等[10]提出的隱式歐拉時(shí)間積分能量最小值表達(dá)式。相對(duì)于使用傳統(tǒng)的牛頓法，本文通過(guò)引入輔助變量(彈簧方向變量)，重新推導(dǎo)隱式歐拉時(shí)間積分能量最小值表達(dá)式，將原有的非線(xiàn)性方程組轉(zhuǎn)化為若干個(gè)小的獨(dú)立的非線(xiàn)性方程，之后使用塊坐標(biāo)下降法求解最佳彈簧方向，即系統(tǒng)能量最小化時(shí)彈簧的方向向量，進(jìn)而根據(jù)最佳彈簧方向向量求得整個(gè)系統(tǒng)的解。由于分解后的系統(tǒng)矩陣是相互獨(dú)立的，因此，可以使用預(yù)先計(jì)算的稀疏Cholesky因子分解矩陣，從而加快解算速度。

1 質(zhì)點(diǎn)-彈簧模型系統(tǒng)

構(gòu)建能夠表達(dá)布料變形特性的力學(xué)模型是服飾動(dòng)畫(huà)首先要考慮的基礎(chǔ)問(wèn)題。質(zhì)點(diǎn)-彈簧模型是一種形式簡(jiǎn)單、適合表現(xiàn)柔性材料動(dòng)態(tài)的模型，目前在布料動(dòng)畫(huà)中得到廣泛應(yīng)用。

首先，假定布料是3D中m個(gè)離散質(zhì)點(diǎn)組成的系統(tǒng)，其還包括一個(gè)離散的時(shí)間集合，形如t1，t2, …，tn，時(shí)間步長(zhǎng)為h(通常設(shè)h=1/30 s)，假定Ρn∈3m作為在時(shí)間tn下質(zhì)點(diǎn)系統(tǒng)的位置，同時(shí)系統(tǒng)遵守牛頓動(dòng)量守恒定律，作用于質(zhì)點(diǎn)上的力由非線(xiàn)性函數(shù)表示為f:3m→3m，f(Ρn)是作用在質(zhì)點(diǎn)上的力向量。質(zhì)點(diǎn)的位置只與這個(gè)力向量有關(guān)。假定力是守恒力，而守恒力在物理上是指該力對(duì)物體所作的功，僅與物體的位置變化量有關(guān)，而與路徑無(wú)關(guān)。常見(jiàn)的守恒力有彈簧力和重力。這里該力表示為f=-E，其中E：3m→為勢(shì)能函數(shù)(非線(xiàn)性且不收斂)。質(zhì)點(diǎn)-彈簧模型的最終目標(biāo)是計(jì)算出彈簧系統(tǒng)質(zhì)點(diǎn)的各個(gè)狀態(tài)Ρ1，Ρ2，…，Ρn。

給定質(zhì)量矩陣M∈3m×3m, 根據(jù)隱式歐拉積分：

Ρn+1=Ρn+hvn+1

(1)

vn+1=vn+hM-1f(Ρn+1)

(2)

式中：vn和vn+1分別為tn和tn+1時(shí)刻的質(zhì)點(diǎn)速度，f(Ρn+1)表示作用在質(zhì)點(diǎn)Ρn+1上的力向量。根據(jù)式(1)和(2)，質(zhì)點(diǎn)速度的表達(dá)式為

hvn=Pn-Pn-1

(3)

hvn+1=Pn+1-Pn

(4)

聯(lián)合式(3)和(4)以消除式(2)中的速度因子，推導(dǎo)出式(5)。

Ρn+1-2Ρn+Ρn-1=h2M-1f(Ρn+1)

(5)

式(5)是牛頓第二定律的離散表達(dá)版本(通常其表示為F=ma)，如果Ρn和Ρn-1的狀態(tài)已知，那么，由式(5)可以求出Ρn+1的狀態(tài)。

關(guān)于非線(xiàn)性系統(tǒng)式(5)的經(jīng)典解法是Baraff等[8]提出的在已知狀態(tài)下力的系統(tǒng)線(xiàn)性化，如式(6)所示。

f(Pn+1)≈f(Pn)+(f|Pn)(Pn+1-Pn)

(6)

M(x-y)=h2f(x)

(7)

由式(7)推導(dǎo)出新的變形隱式積分函數(shù)，如式(8)所示。

(8)

2 彈簧投影過(guò)程

在質(zhì)點(diǎn)-彈簧模型中，傳統(tǒng)胡克定律下的彈性函數(shù)是非線(xiàn)性非凸的，求解整個(gè)大的系統(tǒng)耗時(shí)長(zhǎng)、成本高。但是，將大的系統(tǒng)分解為許多小的可以解決的凸問(wèn)題，不僅易于求解，而且易于并行化，可加快求解速度，這也是本文基于彈簧投影算法的主要思想。因此，這里重新定義能夠使用塊坐標(biāo)下降法的勢(shì)能函數(shù)E的方程，將整個(gè)仿真問(wèn)題轉(zhuǎn)換為求解隱式歐拉能量方程最小值問(wèn)題。其中，引入輔助變量將整個(gè)大的非凸非線(xiàn)性系統(tǒng)分解為若干個(gè)小而獨(dú)立的凸問(wèn)題，接著用塊坐標(biāo)下降法求解新的能量方程。

在質(zhì)點(diǎn)-彈簧系統(tǒng)中，勢(shì)能函數(shù)的主要組成部分是慣性勢(shì)能和彈性勢(shì)能。慣性勢(shì)能是線(xiàn)性方程，比較容易求解，因此本文主要是對(duì)彈性勢(shì)能方程進(jìn)行改變。

根據(jù)胡克定律，彈簧彈性勢(shì)能定義為

(9)

式中：p1,p2為彈簧兩個(gè)端點(diǎn)，p1，p2∈3；r為彈簧原始長(zhǎng)度，r≥0；k為彈簧的硬度參數(shù)，k≥0。

即使是如此簡(jiǎn)單的彈性系統(tǒng)，卻是一個(gè)非凸函數(shù)，非線(xiàn)性和非凸函數(shù)特性使得基于牛頓法的解算代價(jià)比較高。非線(xiàn)性使得系統(tǒng)每次都要重新計(jì)算海森矩陣，而非凸使得海森矩陣是非正定的，這將導(dǎo)致計(jì)算無(wú)法進(jìn)行。從運(yùn)行效率的需求上來(lái)看，理想的系統(tǒng)是一個(gè)凸函數(shù)且有恒定矩陣，或者可以分解成許多小而獨(dú)立的非凸問(wèn)題，能夠并行運(yùn)行。為此引入帶約束的輔助變量d表示彈簧原長(zhǎng)。重新定義勢(shì)能函數(shù)的關(guān)鍵在于，式(9)表示的彈性勢(shì)能函數(shù)是帶約束的能量最優(yōu)化值問(wèn)題。

對(duì)于每個(gè)彈簧的端點(diǎn)p1，p2∈3并且r≥0：

(10)

式(10)的證明如下所述。

證明定義p12∶=p1-p2，考慮到‖d‖=r，重寫(xiě)式(10)的左邊變量為

(‖p12‖-r)2=(‖p12‖-‖d‖)2

根據(jù)逆三角不等式，得出：

(‖p12‖-‖d‖)2≤(‖p12‖-d)2

至此證明了式(10)。

(11)

將所有的彈簧勢(shì)能加在一起，獲得整個(gè)質(zhì)點(diǎn)彈簧系統(tǒng)的彈性勢(shì)能，如式(12)所示。

(12)

經(jīng)過(guò)整理得到式(13)。

(13)

將式(13)轉(zhuǎn)換成線(xiàn)性矩陣之后，得到式(14)。

(14)

式中：s為彈簧總數(shù)；i1,i2∈{1, 2, 3, …,m}為第i個(gè)彈簧端點(diǎn)的編號(hào)；向量X=(p1, …,pm)；矩陣L∈3m×3m，J∈3m×3s的定義如式(15)所示。

(15)

式中：Ai∈m為彈簧i的關(guān)聯(lián)向量，即該彈簧質(zhì)點(diǎn)與其他質(zhì)點(diǎn)是否組成彈簧，如果組成彈簧，則值為1，否則為-1，例如，Ai,i1=1,Ai,i2=-1；Si∈s為彈簧i的方向，即彈簧i在原長(zhǎng)狀態(tài)下的旋轉(zhuǎn)量；I3∈3×3為單位矩陣；?為克羅內(nèi)克積。事實(shí)上，矩陣L是質(zhì)點(diǎn)-彈簧系統(tǒng)圖的剛度權(quán)重圖。

接下來(lái)，加入外力(重力、交互力、碰撞響應(yīng)力)fext∈3m，那么系統(tǒng)的勢(shì)能函數(shù)為

(16)

式中：U={(d1,…,d2)∈2s:||di||=ri}為彈簧原長(zhǎng)狀態(tài)的方向向量，將式(16)代入式(8)，則得式(17)。

(17)

式(17)中的最小值就是隱式歐拉時(shí)間步長(zhǎng)的精確解。

3 阻尼力和碰撞響應(yīng)

碰撞檢測(cè)和響應(yīng)處理是布料模擬的關(guān)鍵和難點(diǎn)問(wèn)題，其解決方法的優(yōu)劣直接影響布料仿真的實(shí)時(shí)性和精確性。

碰撞力也是外力的一部分，可以將其直接加入fext項(xiàng)，進(jìn)而在全局步驟中參與計(jì)算。

在加入阻尼力和碰撞力之后，隱式歐拉方程為

(18)

本文使用塊坐標(biāo)下降法[12]求解式(18)的最小值優(yōu)化問(wèn)題。

本文解算器算法如下：

(1) 初始化X∶=y

(2) 循環(huán)：

1) 局部步驟：給定X，找到對(duì)于所有彈簧最優(yōu)的方向得到d。

2) 全局步驟：固定d，求解X：(M+h2L)X=My+h2Jd+h2fext。

(3) 直到收斂或者超出時(shí)間限制。

從對(duì)X的初始猜測(cè)(這里使用y)開(kāi)始，首先調(diào)整X并計(jì)算d(局部步長(zhǎng))的最佳方向值，其次，調(diào)整d并計(jì)算X的最佳值(全局步長(zhǎng))，重復(fù)此過(guò)程，直到達(dá)到最大迭代次數(shù)為止。與以前方法不同，本文在求解d的局部步驟中不需要奇異值分解，而僅需要向量歸一化。在求解X的全局步驟(固定d)中，需要解決凸二次最小化問(wèn)題。實(shí)際上，因?yàn)長(zhǎng)是對(duì)稱(chēng)的且為正半定的，所以系統(tǒng)矩陣M+h2L是對(duì)稱(chēng)正定。而且隨著迭代開(kāi)始，粒子質(zhì)量、彈簧剛度和連通性保持不變，系統(tǒng)矩陣將保持不變。因此，本文預(yù)先計(jì)算了系統(tǒng)稀疏矩陣的Cholesky因式分子，這將使得線(xiàn)性系統(tǒng)求解速度非?？?。

加入輔助彈簧后彈簧狀態(tài)在整個(gè)系統(tǒng)中經(jīng)歷3個(gè)狀態(tài)之間的轉(zhuǎn)換如圖1所示。圖1(b)局部狀態(tài)對(duì)應(yīng)解算器的局部步驟，可抽象為彈簧處于原長(zhǎng)狀態(tài)，稱(chēng)為最佳彈簧方向，此時(shí)系統(tǒng)彈性勢(shì)能最小，系統(tǒng)處于穩(wěn)定狀態(tài)；圖1(c)全局狀態(tài)對(duì)應(yīng)解算器的全局步驟，根據(jù)局部狀態(tài)下的彈簧方向向量，可以求出質(zhì)點(diǎn)的位置向量，也就是系統(tǒng)的全局解。

4 試驗(yàn)與結(jié)果分析

在布料實(shí)時(shí)模擬中，通常使用一個(gè)常量時(shí)間步長(zhǎng)h確定目標(biāo)幀速率。Baraff等[8]使用了半隱式積分法，該方法在h=1/120、1/60 s時(shí)比較正常，而h=1/30 s時(shí)系統(tǒng)矩陣的解將可能不收斂。為了避免這個(gè)問(wèn)題，Baraff推薦使用自適應(yīng)時(shí)間步長(zhǎng)，但這將導(dǎo)致系統(tǒng)運(yùn)行時(shí)間大大增加，系統(tǒng)代價(jià)過(guò)大。

一般而言，在實(shí)時(shí)系統(tǒng)中，h=1/30 s時(shí)仿真效果比較能讓人接受，因此本文采用時(shí)間步長(zhǎng)h=1/30 s。

本文方法與牛頓法針對(duì)一塊2 500個(gè)質(zhì)點(diǎn)的布料進(jìn)行模擬數(shù)值解的比較，結(jié)果如圖3所示，其中，相對(duì)誤差定義如式(19)。

(19)

式中：x0為初始狀態(tài)值；xi為當(dāng)前狀態(tài)值；x*為精確解的值。

由圖2(a)可知，本文方法整體呈現(xiàn)線(xiàn)性收斂，牛頓法呈現(xiàn)二次收斂[13]。但是，在有限的計(jì)算成本下，牛頓法只迭代1次，在這種情況下，本文方法只迭代一次的相對(duì)誤差小于牛頓法(圖2(a)中的A點(diǎn))。也就是說(shuō)，在迭代的初始階段，本文方法的表現(xiàn)優(yōu)于牛頓法。由圖2(b)可知，本文方法在迭代100次以?xún)?nèi)數(shù)值解的相對(duì)誤差低于牛頓法的相對(duì)誤差(如圖2(b)中的C點(diǎn))。表1列出了本文算法在A點(diǎn)、B點(diǎn)、C點(diǎn)的時(shí)間效率和相對(duì)誤差，其分別代表布料動(dòng)畫(huà)在迭代1、10、100次時(shí)的運(yùn)算結(jié)果。牛頓算法的1次迭代相對(duì)誤差和迭代時(shí)間分別為0.250和181 ms。由表1可知，本文算法的相對(duì)誤差更小，并且運(yùn)行時(shí)間更短。

表1 本文算法的時(shí)間效率及相對(duì)誤差

雖然本文方法的整體收斂速度無(wú)法趕上牛頓法的整體二次收斂，但是，本文方法在在第一階段(阻尼階段)優(yōu)于牛頓法。當(dāng)?shù)膞i值接近于x*時(shí)，將本方法的結(jié)果作為牛頓法的一個(gè)起點(diǎn)，可以加快計(jì)算速度。

當(dāng)時(shí)間步長(zhǎng)=1/30 s時(shí)，使用本文方法和牛頓法模擬同一布料的動(dòng)畫(huà)效果和幀速率表現(xiàn)，如圖3所示。由圖3可知：本文方法的一次迭代已經(jīng)產(chǎn)生穩(wěn)定合理的模擬(見(jiàn)圖3(a))，但是褶皺細(xì)節(jié)相對(duì)不明顯；10次迭代在速度與質(zhì)量上得到比較好的平衡，如圖3(b)所示。

在幀速率表現(xiàn)方面，由圖3可知，本文方法迭代1次和10次的模擬幀速率約為30幀/s，滿(mǎn)足實(shí)時(shí)動(dòng)畫(huà)要求，而牛頓法的1次迭代的幀速率只有7.6幀/s，無(wú)法應(yīng)用在實(shí)時(shí)應(yīng)用中，只能離線(xiàn)渲染。一般而言，幀速率越大，動(dòng)畫(huà)畫(huà)面看起來(lái)越流暢。由上述分析可知，本文方法在長(zhǎng)時(shí)間步的布料模擬系統(tǒng)中表現(xiàn)良好。

當(dāng)模擬的布料有2 500個(gè)質(zhì)點(diǎn)，共有4 802個(gè)三角形，時(shí)間步長(zhǎng)=1/30 s，在不同迭代次數(shù)下，PBD法的布料模擬效果如圖4所示。由圖4可以觀(guān)察到，在相同的時(shí)間步長(zhǎng)下，基于PBD的布料模擬迭代次數(shù)為1時(shí)，布料出現(xiàn)過(guò)拉伸現(xiàn)象，隨著迭代次數(shù)成倍增加，布料變得更加“僵硬”。這是由于PBD在增加迭代次數(shù)時(shí)會(huì)增加系統(tǒng)的有效剛度，導(dǎo)致剛度參數(shù)與標(biāo)準(zhǔn)胡克模型不兼容。但是本文方法不存在這個(gè)問(wèn)題，即隨著迭代次數(shù)成倍增加并沒(méi)有增加系統(tǒng)的有效剛度，實(shí)際表現(xiàn)為布料模擬效果不會(huì)過(guò)度受迭代次數(shù)的影響。

在碰撞處理方面，本文設(shè)計(jì)了一個(gè)窗簾模型布料與小球交互的試驗(yàn)，以測(cè)試本文的方法在接觸和碰撞情況下的行為，實(shí)際結(jié)果如圖5所示。

由圖5可以看到布料自碰撞和小球交互碰撞的效果，布料的垂墜感和褶皺基本符合現(xiàn)實(shí)中布料的物理特性。由此證實(shí)本文方法的有效性和實(shí)用性。

5 結(jié) 語(yǔ)

本文基于質(zhì)點(diǎn)-彈簧投影模型的布料仿真算法提出了一種改進(jìn)的模型求解非線(xiàn)性數(shù)值方法，用于隱式歐拉時(shí)間步長(zhǎng)的質(zhì)點(diǎn)-彈簧系統(tǒng)。該算法在布料模擬的動(dòng)力學(xué)方程中引入輔助彈簧方向變量，將原有的非線(xiàn)性方程組轉(zhuǎn)化為若干個(gè)小的非凸方程，然后，利用塊坐標(biāo)下降法來(lái)求解每一個(gè)非凸方程，從而實(shí)現(xiàn)快速的迭代。雖然本文方法的整體收斂性不如牛頓法，但是，在牛頓法只迭代一次的情況下，本文方法的數(shù)值解相對(duì)于精確值的誤差比牛頓法要小，相對(duì)誤差曲線(xiàn)斜率更大，即在布料模擬初始階段，本文方法的求解速度更快，并且數(shù)值解更接近精確值。本文的方法非常適用于對(duì)實(shí)時(shí)性要求高的系統(tǒng)，也可以將本文的方法應(yīng)用于系統(tǒng)的初始仿真階段，后續(xù)可以繼續(xù)采用牛頓法用于仿真，有利于提高系統(tǒng)的整體仿真速度。同時(shí)，該方法也克服了PBD方法的布料模擬剛度受約束投影迭代次數(shù)影響的缺點(diǎn)，具備良好的實(shí)時(shí)性和穩(wěn)定性。