張藝衡，李曉峰，2，劉小玲，楊智宇

（1.四川大學(xué) 計算機(jī)學(xué)院，四川 成都610064；2.四川大學(xué) 視覺合成圖形圖像技術(shù)國防重點學(xué)科實驗室，四川 成都610064）

0 引 言

對于可交互水波等流體現(xiàn)象的模擬，目前比較典型的方法是數(shù)值求解納維斯托克方程。該方法雖然能夠模擬出具有高真實感的水波，但其數(shù)值求解過程需要大量的計算資源，不適用于大規(guī)模場景中的可交互水波實時模擬。另外，基于波動方程數(shù)值求解的水波模擬也是目前應(yīng)用比較廣泛的算法之一，與基于納維斯托克方程方法類似的是，其數(shù)據(jù)求解的過程依然會產(chǎn)生龐大的計算開銷，算法的繪制效率偏低。這2類方程相較而言，波動方程更加簡單，對于流體的表面模擬也能達(dá)到比較真實的效果。因此，本文將波動方程作為算法的理論基礎(chǔ)，并對該方程進(jìn)行近似建模，構(gòu)建水面高度場，并結(jié)合粒子系統(tǒng)進(jìn)行波浪的實時更新，以實現(xiàn)具有高效率和高真實感的可交互水波實時模擬。

1 相關(guān)工作

在水體模擬的早期研究中，主要是針對水面形態(tài)及運(yùn)動狀態(tài)直接進(jìn)行顯示建模，這類方法一般采用三角函數(shù)、噪音函數(shù)和傅里葉合成公式來進(jìn)行水表面的模擬建模。Schneider等人利用Perlin噪音函數(shù)實現(xiàn)了基于GPU的實時水波模擬。Tessendorf提出采用傅里葉合成法進(jìn)行海浪運(yùn)動的模擬，并且該方法成為目前模擬環(huán)境波的主流算法之一。但是，該算法不能模擬水面浮游固體與水面之間的交互運(yùn)動。為此，Tessendorf提出基于該算法的改進(jìn)算法，實現(xiàn)了浮游固體對水面的單向交互。水面高度場構(gòu)建法也是水體模擬的主流算法之一，這類算法主要以淺水方程作為模擬的理論基礎(chǔ)，出現(xiàn)了許多研究成果。T.R.Hagen等人采用半拉格朗日求解法求解淺水方程構(gòu)造水面高度圖，但該方法在效果逼真度方面的損耗較大。為此，B.M.Kim等人對該算法進(jìn)行了改進(jìn)，取得了一定程度的算法逼真度的提高。Daniel等人提出了一種基于淺水方程的自適應(yīng)網(wǎng)格模擬算法，實現(xiàn)了大規(guī)模水面的模擬［1］。Hyokwang等人采用二維SPH實現(xiàn)了水面高度圖的創(chuàng)建［2］。Cem Yukel提出WP（wave particles）概念，將水波離散成滿足波動方程的波浪粒子，以構(gòu)建水面高度場［3］。Hilko等人對WP算法進(jìn)行了擴(kuò)展，實現(xiàn)了水流效果的模擬［4］。另外，歐拉網(wǎng)格法在水波模擬方面也出現(xiàn)了許多研究成果。Enright等人提出了基于歐拉網(wǎng)格的PLSM算法，Schroeder介紹了對該算法的一些改進(jìn)，在一定程度上降低了計算復(fù)雜度［5］。Andrew等人采用超網(wǎng)格單元 （tall grid cells）思想，將N－S方程與高度場相結(jié)合，減少了單元網(wǎng)格的數(shù)量，實現(xiàn)了固體與水之間的交互［6］。四面體網(wǎng)格在歐拉法流體模擬中也得到了應(yīng)用［7－9］。除此之外，拉格朗日粒子法在水模擬方面也得到了非常廣泛的應(yīng)用，如經(jīng)典SPH流體模擬法。Becker等人基于SPH算法，并對其進(jìn)行了改進(jìn)，實現(xiàn)了體積近似恒定的流體模擬［10］。Solenthaler等人則提出了預(yù)期修正SPH算法 （PCISPH），保持穩(wěn)定時間步長和迭代求解，使粒子系統(tǒng)趨向于不可壓縮狀態(tài)［11］。Clavet等人實現(xiàn)了基于拉格朗日粒子的粘性流體的實時交互模擬［12］。拉格朗日粒子法易于表達(dá)和實現(xiàn)，相對于歐拉網(wǎng)格法，粒子運(yùn)動具有自由和隨機(jī)性，不需要對空間進(jìn)行有限網(wǎng)格的定義。但拉格朗日法隨著粒子數(shù)量的增加，其計算量往往會隨之增大。

2 算法思想及設(shè)計

2.1 高度場離散化

對于水面波浪的模擬，本文采用表面高度場的方法進(jìn)行。對于一個完整的水面高度場，通常包含了數(shù)量不定的多個波浪高度信息，對整體高度場進(jìn)行直接構(gòu)建的方法可行性不高。因此，本文采用高度場離散化的思想來構(gòu)建水面高度場。即將整體高度場以波浪為單位進(jìn)行劃分，更進(jìn)一步將單個波浪劃分成多個波段，形成多個小的局部高度場，從而將問題簡化為對局部高度場的獨立構(gòu)建。如圖1所示，對于左圖中的完整波浪，本文將其離散為右圖中方框標(biāo)記的局部高度場，并對局部高度場進(jìn)行一定規(guī)則的排列控制，以形成左圖所示的完整平滑波浪高度場。

圖1 高度場離散化

2.2 局部高度場建模

本文采用離散化的構(gòu)建思想，將完整水波離散為多個離散小型水波波段，并結(jié)合粒子系統(tǒng)思想，每個波段對應(yīng)一個粒子。通過更新粒子屬性的方式，達(dá)到模擬跟蹤水波運(yùn)動形態(tài)的目的。在現(xiàn)實生活中，水波由橫波和縱波2個部分組成，分別對應(yīng)水波的垂直偏移和水平偏移。與此類似，在本文中，每個局部波段對應(yīng)一個局部高度場，由垂直偏移和水平偏移兩部分組成，分別對應(yīng)一個局部水平偏移函數(shù)和一個局部垂直偏移函數(shù)。

假設(shè)局部垂直偏移函數(shù)為D（X，t），則該函數(shù)表示在t時刻，水面上X位置處的垂直偏移，即相對于水面水平位置的高度值。假設(shè)水面高度場為函數(shù) H（x，t），則 H（x，t）的計算公式如下所示

式中：h0——水面的水平高度，φz（X，t）——t時刻，p（x，y）位置處由水波引起的垂直偏移，Di——第i個粒子表示的垂直偏移。對于D（X，t），本文提出一種徑向定義方式，即將每個局部高度場構(gòu)造為底面為圓形的弧面，每個弧面以波浪圓心為中心，隨著波浪的傳播進(jìn)行徑向運(yùn)動，以確保各個弧面之間的相對獨立，以充分利用粒子系統(tǒng)進(jìn)行整體波形及運(yùn)動狀態(tài)的模擬。具體構(gòu)造方式如下所示

式中：ai——局部高度場對應(yīng)波段的振幅，ri——弧面底部圓面的半徑，｜x－xi（t）｜——圓面上某點與圓心之間的距離。通過這種方式構(gòu)造的局部高度場，以一定的間隔進(jìn)行弧形排列，即能構(gòu)造出一個完整逼真的波形高度場。

局部垂直偏移函數(shù)反應(yīng)了波形在垂直方向上的形態(tài)，局部水平偏移函數(shù)則模擬了波浪運(yùn)動過程中在水平方向上的變化，如下所示

2.3 波粒子系統(tǒng)

本文算法的關(guān)鍵思想在于離散水面高度場，通過控制各個離散高度場的高度變化及運(yùn)動方向控制，以達(dá)到模擬真實水波的形態(tài)變化及運(yùn)動規(guī)律。對此，本文引入粒子系統(tǒng)，進(jìn)行水波運(yùn)動變化的實時跟蹤。即每個局部高度場對應(yīng)一個粒子，粒子的屬性存儲了各個高度場的信息，包括產(chǎn)生時間，產(chǎn)生位置，擴(kuò)散角，傳播方向以及振幅。隨著水波的傳播擴(kuò)散，通過更新波粒子的屬性來模擬水波的運(yùn)動規(guī)律及形態(tài)變化。但是，隨著水波的擴(kuò)散，波粒子進(jìn)行徑向運(yùn)動，粒子之間的距離也隨之增大，當(dāng)距離達(dá)到一定閾值時，水波波面會出現(xiàn)相對分離的現(xiàn)象，即出現(xiàn)失真的情況，如圖1右圖所示。為了避免這樣的失真現(xiàn)象，隨著水波的擴(kuò)散，本文通過分裂波粒子的方式，在粒子之間添加新粒子從而控制相鄰粒子之間的間隔距離，如圖2所示。為了提高算法效率，本文提出一種基于母粒子的更新方法。

圖2 波粒子分裂

2.4 基于母粒子的粒子更新

在本文算法中，一個完整的水波由多個波粒子所代表的圓底弧面，以一定的均勻間隔排列形成。即同一個水波上的波粒子，具有相等的振幅及運(yùn)動速度。因此，同一水波上的波粒子也具有同時分裂和同時消亡的特征。據(jù)此，本文采用二維的方式組織粒子，即每個水波對應(yīng)一個母粒子，在進(jìn)行粒子系統(tǒng)更新時，只需要檢測母粒子，當(dāng)母粒子達(dá)到分裂或消亡的條件時，則分裂或消亡母粒子下的所有子粒子。除此之外，粒子振幅的更新也只需要針對母粒子進(jìn)行實時更新計算，所有子粒子與母粒子具有相等的振幅，從而避免了遍歷更新所有的波粒子，極大地降低了計算量，實現(xiàn)高效率的粒子系統(tǒng)更新。至于粒子的實時坐標(biāo)位置，本文通過對粒子屬性進(jìn)行巧妙的設(shè)計，通過設(shè)置粒子的產(chǎn)生時間，運(yùn)動速率及擴(kuò)散角，使其可充分利用可編程著色語言GLSL，在GPU上進(jìn)行高效快速的并行計算。避免了在CPU上對每個粒子進(jìn)行坐標(biāo)位置的實時計算，利用GPU強(qiáng)大的并行計算能力達(dá)到高效的計算效率。

2.5 水的光學(xué)特性模擬

為了實現(xiàn)具有高真實感的水面效果，本文基于水面的基本光照模擬，引入環(huán)境貼圖技術(shù)。基于基本光照物理模型，水面光照模型如圖3所示。水面基本的光照效果主要是折射和反射，且其計算方法基于真實水面的折射和反射物理定律。

在現(xiàn)實中，水面對周圍環(huán)境具有非常強(qiáng)的反射性，對水面反射性的模擬對水面的真實感具有非常重要的作用。在虛擬仿真領(lǐng)域，對水面的反射效果模擬最為精確且經(jīng)典的技術(shù)是光線跟蹤技術(shù)，但是該技術(shù)的計算開銷非常大，不適用于實時繪制系統(tǒng)中。為此，本文采用立方體環(huán)境貼圖技術(shù)進(jìn)行模擬。該技術(shù)不僅可以充分利用GLSL的可編程屬性高效地模擬出高質(zhì)量的光照效果和環(huán)境貼圖，并且基于該技術(shù)可以比較方便地生成遠(yuǎn)反射效果。

3 算法實現(xiàn)

第2小節(jié)介紹了本文的算法思想，本節(jié)將介紹波粒子系統(tǒng)與高度場結(jié)合模擬水波的算法實現(xiàn)過程。本文的波粒子系統(tǒng)主要負(fù)責(zé)根據(jù)當(dāng)前水面上的交互信息產(chǎn)生新粒子，并實時更新系統(tǒng)中波粒子的位置，處理粒子傳播過程中的分裂和反彈。反彈主要針對波浪遇到水體邊界時的處理。對于粒子系統(tǒng)的迭代更新，傳統(tǒng)的基于CPU的更新方式實現(xiàn)簡單，但效率很低。為了充分利用GPU強(qiáng)大的并行計算能力，本文對粒子屬性進(jìn)行了巧妙的設(shè)計，使得粒子的更新處理可以全部在GPU上進(jìn)行，獲得了非常高效的計算效率。

3.1 母波粒子的產(chǎn)生

對于本文基于波粒子系統(tǒng)的水波，水波從產(chǎn)生到擴(kuò)散，最后到消失的過程實際上對應(yīng)了波粒子的生成，分裂和消失的過程。若水面上某處發(fā)生了交互，將對應(yīng)一個波粒子的生成，其屬性包括：生成時間、生成位置、擴(kuò)散角、傳播方向和振幅。通過這幾個屬性的設(shè)計，使得粒子的更新可以完全通過GPU來進(jìn)行，具體方式將在3.2小節(jié)進(jìn)行介紹。

3.2 波粒子更新

波粒子的運(yùn)動基于二維空間且粒子之間相對獨立，不存在交互。除此之外，沒有其它外力作用在該粒子系統(tǒng)上，因此波粒子具有相對恒定的運(yùn)動速率，并且除了與水體邊界發(fā)生的碰撞反彈之外，粒子的運(yùn)動方向保持不變。基于本文對粒子屬性的設(shè)計，不需要每幀對波粒子的位置進(jìn)行更新存儲，而是通過GLSL實時計算獲得，如下所示

式中：v——波粒子的運(yùn)動速率， ——波粒子的傳播方向。由于粒子數(shù)量隨著波浪的傳播呈指數(shù)級增長，甚至導(dǎo)致其更新存儲所需的CPU緩存大于系統(tǒng)的CPU緩存，以致算法異常退出。通過利用GPU實時并行計算的方式，避免了這部分內(nèi)存的讀寫開銷，而消耗小部分式 （4）所需的乘法和加法開銷，確保了算法的高效性和穩(wěn)定性。

在現(xiàn)實中，水波在傳播的過程中，由于運(yùn)動擴(kuò)散和阻尼導(dǎo)致其振幅的不斷降低直至降低到0。為此，本文采用式（5）來進(jìn)行模擬

式中：δ——阻尼系數(shù)，t——當(dāng)前時刻，t0——波浪粒子的產(chǎn)生時刻。

3.3 波粒子分裂

隨著水波的傳播，同一水波上相鄰波粒子之間的間隔增大，為了確保波形的平滑連續(xù)，本文對相鄰波粒子之間間距大于波粒子半徑1／2的粒子進(jìn)行分裂處理，分裂出2個新粒子，分布在其左右兩邊。而粒子一旦產(chǎn)生，其分裂時間可以通過公式進(jìn)行計算，如下所示

因此，水波傳播過程中，不需要計算粒子之間的間隔，避免了粒子之間的關(guān)聯(lián)，從而將復(fù)雜問題轉(zhuǎn)化成了簡單的時間檢測。除此之外，本文基于粒子的二維粒子系統(tǒng)，極大地減少了檢測次數(shù)，達(dá)到了高效的檢測效率。

在波粒子分裂為3個波粒子的過程中，新粒子振幅為分裂粒子振幅的1／3，同時，擴(kuò)散角縮小為原擴(kuò)散角的1／3，如圖4所示。即保證了波浪的波陣面的平滑連續(xù)，也模擬了現(xiàn)實中水波振幅隨著傳播而降低的現(xiàn)象。

3.4 交互水面繪制

交互水面的繪制則結(jié)合波粒子系統(tǒng)對水面的形態(tài)采用GLSL完成最終的繪制。本文設(shè)計了4個GPU自定義繪制管線來完成，即4對GLSL。第一個繪制管線根據(jù)CPU傳入的粒子信息創(chuàng)建粒子信息圖，如圖5所示，該圖存儲粒子的振幅及傳播方向信息，并傳入第2個繪制管線中；在第二個繪制管線中，結(jié)合第1個管線產(chǎn)生的粒子信息圖，采用式 （2）對信息圖中的振幅在X軸方向上進(jìn)行過濾；在第3個繪制管線中，則對信息圖在Y軸方向上進(jìn)行過濾；第4個繪制管線則結(jié)合兩個方向上過濾出來的高度圖，結(jié)合式 （3）進(jìn)行圓形平滑過濾，使得高度場平滑連續(xù)，滿足波形特征。

圖5 粒子信息

4 實驗結(jié)果及分析

本文采用基于母粒子的波浪粒子系統(tǒng)，并且利用GPU進(jìn)行波浪粒子屬性的更新，取得了非常高效的繪制效率，如表1和圖6所示。表1是一個幀率表，其中，行表頭表示水波數(shù)量，列表頭表示水波同時分裂的次數(shù)，為了便于統(tǒng)計和控制，測試水波都是同時產(chǎn)生的，因此分裂也是同時進(jìn)行。本文實驗采用的屏幕渲染尺寸為1440×900，水面網(wǎng)格尺寸為1024×1024。

表1 不同波浪數(shù)量和分裂次數(shù)下的運(yùn)行幀率／（ftp／s）

效果圖如圖7和圖8所示。從圖中可以看出，本文算法實現(xiàn)的水面具有較高的真實感，且水波波形也具有非常高的逼真度。

5 結(jié)束語

本文設(shè)計了一個基于GPU的實時可交互模擬算法，對水面高度場進(jìn)行離散化處理，基于現(xiàn)實中水波的組成，設(shè)計了局部垂直和局部水平偏移函數(shù)，進(jìn)行水波波形的模擬。同時，通過離散化處理后，結(jié)合粒子系統(tǒng)，實現(xiàn)了水波的動態(tài)實時更新，而基于母粒子的波粒子系統(tǒng)的提出和應(yīng)用以及波粒子屬性的合理設(shè)計，極大地降低了遍歷開銷，并將粒子的更新全部放到GPU上進(jìn)行，實現(xiàn)了高效地更新計算。為了提高水面的真實感，本文結(jié)合水的光學(xué)特性，基于水面光照模型，對水面進(jìn)行了光照計算，并采用立方體貼圖技術(shù)模擬水面對周圍環(huán)境的反射效果。實驗結(jié)果表明，本文算法實現(xiàn)了高效穩(wěn)定的可交互水波實時模擬。

圖7 鼠標(biāo)滑動觸發(fā)的水波效果

圖8 小船在水面運(yùn)動導(dǎo)致的水波

［1］Kallin D.Real－time large scale fluids for games ［C］／／SIGRAD，2008：31－35.

［2］Lee H，Han S.Solving the shallow water equations using 2D SPH particles for interactive applications ［J］.The Visual Computer，2010，26 （6－8）：865－872.

［3］Yuksel C，House D H，Keyser J.Wave particles ［J］.ACM Transactions on Graphics（TOG），2007，26 （3）：991－998.

［4］Cords H.Moving with the flow：Wave particles in flowing liquids［C］／／International Conference in Central Europe on Computer Raphics，Visualization and Computer Vision，2008：145－152.

［5］Craig Schroeder.Semi－implicit surface tension formulation with a Lagrangian surface mesh on an Eulerian simulation grid ［J］.Journal of Computational Physics，2012，231 （4）：2092－2095.

［6］Andrew Selle，Ronald Fedkiw，ByungMoon Kim，et al.An unconditionally stable MacCormack method ［J］.Journal of Scientific Computing，2008，2008 （35）：350－371.

［7］Zhang Yizhong，Wang Huamin，Wang Shuai，et al.A deformable surface model for real－time water drop animation ［J］.IEEE Transactions on Visualization and Computer Graphics，2012，18 （8）：1281－1289.

［8］Chentanez N，F(xiàn)eldman B E，Labelle F，et al.Liquid simulation on lattice－based tetrahedral meshes ［C］／／Proceedings of the ACM SIGGRAPH／Eurographics Symposium on Computer Animation.Eurographics Association，2007：219－228.

［9］Nipun Kwatra，Chris Wojtan，Mark Carlson，et al.Fluid simulation with articulated bodies ［J］.IEEE Transactions on Visualization and Computer Graphics，2010，16 （1）：70－80.

［10］Becker M，Teschner M.Weakly compressible SPH for free surface flows ［C］／／Proceedings of the ACM SIGGRAPH／Eurographics Symposium on Computer Animation.Eurographics Association，2007：209－217.

［11］Solenthaler B，Pajarola R.Predictive－corrective incompressible SPH ［C］.ACM Transactions on Graphics.ACM，2009，28 （3）：401－406.

［12］NVIDIA.Physx page ［EB／OL］.［2010－12－15］.http：／／www.nvidia.com／object／physx＿new.html.