用MC算法優(yōu)化基于SPH的海浪粒子模型

劉文龍，張 靜

（1.中國地震局第一監(jiān)測中心，天津 300180；2.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266510）

用MC算法優(yōu)化基于SPH的海浪粒子模型

劉文龍1，張 靜2

（1.中國地震局第一監(jiān)測中心，天津 300180；2.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266510）

結(jié)合風(fēng)浪場復(fù)雜的動力學(xué)、時空特性以及對目前海浪仿真方法的研究總結(jié)，針對海浪仿真模擬中存在的真實感不足，難以模擬海浪破碎等問題，根據(jù)流體的物理特性，采用光滑粒子流體動力學(xué)（SPH）的方法實現(xiàn)了海浪的動態(tài)演變仿真。然后針對基于SPH方法的海浪粒子模型存在的離散性和真實感較差的問題，采用一種基于移動立方體法（MC）的海浪運動的自由表面抽取算法，完成了海浪場的表面建模，優(yōu)化了基于粒子系統(tǒng)的海浪模擬效果。實驗表明，用MC算法優(yōu)化的方法既滿足了海浪動態(tài)演變的仿真要求，又使水體模型更加連貫，更接近真實水體效果，優(yōu)化效果較為明顯，是綜合利用體繪制與面繪制法實現(xiàn)三維可視化模擬的一次成功的探索。

海浪模擬；光滑粒子流體動力學(xué)；移動立方體法

海浪仿真是近些年來計算機圖形學(xué)領(lǐng)域研究的熱點之一，它在虛擬現(xiàn)實應(yīng)用、計算機游戲以及電影制作中起著重要的作用[1]。另外，海浪的仿真模擬是構(gòu)建海洋環(huán)境信息可視化平臺的一個重要組成部分，平臺的建立可顯著提高我國海洋信息的可視化水平。

海浪具有復(fù)雜的動力學(xué)、時空特性，對其三維模擬有著特殊的要求，這使得海浪建模及繪制成為海浪仿真模擬中的難點。海浪建模的本質(zhì)是海浪的真實感繪制，它是海洋環(huán)境仿真的基礎(chǔ)，在很大程度上決定了仿真效果的好壞[2]。針對海浪復(fù)雜的特性，本文在基于粒子的拉格朗日法（SPH方法）建立的海浪粒子模型基礎(chǔ)上采用移動立方體（MC算法）算法對海浪粒子表面建模，在海浪動態(tài)演變的連貫性方面做出改進，形成更為逼真的海浪動態(tài)演變效果。

1 基于SPH的海浪建模

SPH方法是一種由一組粒子代替流體獲得流體動力學(xué)公式的近似數(shù)學(xué)解決方案的方法[3]，該方法的基本思想是將連續(xù)的流體(或固體)用相互作用的質(zhì)點組來描述[4]，各個物質(zhì)點上承載各種物理量，包括質(zhì)量、速度等，通過求解質(zhì)點組的動力學(xué)方程和跟蹤每個質(zhì)點的運動軌道，求得整個系統(tǒng)的力學(xué)行為。

SPH方法是一種純Lagrange方法，能避免Euler描述中歐拉網(wǎng)格與材料的界面問題，特別適合于求解高速碰撞等動態(tài)大變形問題[5]。因此，SPH是一種適合于三維流體的模擬方法。采用SPH方法進行海浪建模時必須首先將其粒子化。這些海浪粒子在仿真初始的時候給定，所有的粒子都使用統(tǒng)一的光滑核半徑，而且都具有質(zhì)量、密度、速度、加速度、位置等屬性，這些屬性在模擬仿真過程中會不斷地發(fā)生變化。

本文在.NET環(huán)境下利用OpenGL設(shè)計并實現(xiàn)一個粒子系統(tǒng)。SPH海浪運動模擬主要包括以下3個方面分析建立海浪運動模型。

（1）初始化粒子的屬性，生成基本的海浪粒子類；

（2）分析海浪粒子受力情況，計算粒子速度、加速度等屬性的變化，模擬海浪的運動狀態(tài)；

（3）計算粒子位置，判斷是否會與周圍環(huán)境障礙物發(fā)生碰撞；若發(fā)生了，重新計算粒子加速度、速度等變量，以確定粒子碰撞后的位置。圖1為渲染效果圖。

圖1 SPH渲染效果

2 基于MC的粒子群表面建模

由于SPH方法的局限性，從基于SPH方法的海浪粒子系統(tǒng)的模擬中可以看出海浪模型只是一堆離散的水粒子，流動水體的連貫性和真實感較差。針對此問題，提出了一種基于移動立方體法（MC算法）的海浪運動的自由表面抽取算法，對海浪場的表面建模。

MC算法是應(yīng)用最廣泛的等值面方法,其研究主要集中在算法的改進和算法在不同領(lǐng)域的應(yīng)用上。MC算法在三維流體可視化中的應(yīng)用己成為計算機圖形學(xué)中一個重要研究熱點[6]。MC算法的不斷完善和發(fā)展,使得提供直觀、逼真而且能夠包含原始信息中隱含的豐富的三維信息成為可能。在混合流和多相流的模擬中,不同流體間的界面繪制極其重要,利用Ma方法來進行繪制將成為一個廣泛的應(yīng)用方法。將計算流體力學(xué)的數(shù)值方法MC方法結(jié)合起來完成對流體三維可視化的模擬對分析結(jié)果的準(zhǔn)確性和正確性有深遠意義,并促進了應(yīng)用領(lǐng)域的發(fā)展[7]。

密度是海浪粒子的一個重要屬性，通過比較密度值的大小可確定海浪與大氣環(huán)境之間的交界面，即海浪的自由運動表面，這在海浪仿真模擬中起著重要的作用。本文采用MC算法從海浪密度場中提取海浪表面，主要包括下面幾個方面的內(nèi)容。

2.1 基于立體網(wǎng)格的密度場數(shù)據(jù)采樣

根據(jù)粒子的空間分布建立三維空間網(wǎng)格，將立體網(wǎng)格的節(jié)點作為海浪場的采樣點，生成規(guī)則分布的海浪密度場數(shù)據(jù)。網(wǎng)格節(jié)點的密度值可以通過SPH密度插值公式來計算。SPH方法核心是核函數(shù)，它表示在一定的光滑長度范圍內(nèi)其它臨近粒子質(zhì)點對研究粒子影響程度的權(quán)函數(shù)。假設(shè)流體中某點r，在光滑核半徑h范圍內(nèi)受數(shù)個粒子影響，其位置分別是r0,r1,r2…,rj，該點位置處某項屬性A的累加值就可以用式（1）來表示。

式中：Aj是要累加的某種屬性，mj和ρj是周圍粒子的質(zhì)量和密度；rj是該粒子的位置，h是光滑核半徑；函數(shù)W為光滑核函數(shù)。

根據(jù)式（1），用密度ρ代替式中的A，可以得到密度的計算公式：

網(wǎng)格節(jié)點的密度值計算使用的光滑核函數(shù)為Poly6函數(shù)，所有網(wǎng)格的質(zhì)量相同且都為m，在三維情況下，ri處的密度計算公式可表示為：

2.2 基于MC算法的密度等值面提取

MC算法的基本思想是將整個三維數(shù)據(jù)所在空間按照一定的規(guī)則劃分成一個個體元，立方體元8個頂點都擁有自己的函數(shù)值，根據(jù)給定的閾值對體元12條邊進行插值，構(gòu)造體元內(nèi)部的三角面片，連接所有體元的三角面片完成等值面的提取。

在離散的三維空間數(shù)據(jù)場中，通過每次讀取連續(xù)兩張切片數(shù)據(jù)形成一個層，兩張切片上下相對應(yīng)的8個點構(gòu)成一個立方體體元，其中體元各頂點具有自己的屬性值。這樣整個三維空間可以視為由多個擁有頂點值的規(guī)則體元構(gòu)成，如圖2所示為三維規(guī)則數(shù)據(jù)場與體元。在MC算法中，需要給定一個等值面的閾值，將閾值依次與體元的8個頂點作比較，根據(jù)比較結(jié)果將各頂點標(biāo)記為內(nèi)部和外部兩種狀態(tài)，可以理解為二進制值0和1。如果一條邊的兩個頂點分別屬于這兩種情況，則說明這條邊與等值面有交點。

圖2 三維空間數(shù)據(jù)場及體元

由于每個體元的8個頂點分別可能具有兩種狀態(tài)，因此每個體元進行比較后有28=256種情況。經(jīng)過分析，這256種情況存在對稱性。如果將一個體元頂點的兩種狀態(tài)互換后，等值面的連接是不變的（如圖3（a）），根據(jù)體元的這種互補對稱性，256種不同情況的種類將減少到32種。另外，體元8個頂點存在旋轉(zhuǎn)對稱性，即旋轉(zhuǎn)立方體體元后，許多情況是同構(gòu)的（如圖3（b）），這樣將不同情況進一步組合，可以減少到16種情況，這16種情況包括所有頂點均在等值面內(nèi)與均在等值面外兩種，這兩種情況下體元12條邊與等值面均無交點，拓?fù)浣Y(jié)構(gòu)一致，可以合并為一種情況。

圖3 互補對稱性與旋轉(zhuǎn)對稱性

經(jīng)過合并后15種情況如圖4所示，把這15種基本的情況經(jīng)旋轉(zhuǎn)，倒置后可以得到完整的256種情況。

圖4 體元插值的15種基本情況

基于上面的分析，在實現(xiàn)方法上，可以由立方體體元的8個頂點的內(nèi)部和外部兩種情況得到一個0～255之間的索引值，如圖5所示，MC方法用一個字節(jié)的空間構(gòu)造了一個體元的狀態(tài)表，表中每一位代表相應(yīng)頂點的狀態(tài)（0或者1）。根據(jù)這個索引表和體元的互補對稱性與正反對稱性，可以知道體元屬于圖4中的哪一種情況，以及等值面與體元哪條邊相交。

圖5 體元頂點狀態(tài)表

式中：V代表所求交點的位置坐標(biāo)；V0，V1代表體元邊界的兩個頂點坐標(biāo)；C代表給定的閾值，Value0與Value1代表體元邊界兩個頂點的函數(shù)值。

采用MC算法從海浪密度場中抽取海浪表面，首先須確定海浪表面的密度值，以此密度值作為MC算法中插值用的閾值。然后根據(jù)MC算法原理，將每個體元8個頂點的密度值與水面密度值進行比較，確定出與閾值有交點的體元邊界；通過線性插值，計算閾值與體元邊界交點的坐標(biāo)值，并利用中心差分法求體元頂點處的法向量，然后同樣利用線性插值方法求出三角面片頂點處的法向量；最后根據(jù)三角面片的頂點坐標(biāo)及法向量繪制等值面，完成海浪表面的提取。

首先，根據(jù)體元與等值面相交的256種情況構(gòu)建狀態(tài)表。此表包括索引、指向256種相交邊情況的指針以及三角剖分模式。其中，索引指的是標(biāo)記8個頂點狀態(tài)的有序二進制編碼。

然后，確定海浪粒子群的空間范圍，構(gòu)造立方體體元。通過遍歷所有海浪粒子獲取粒子三維坐標(biāo)位置的最大最小值，然后設(shè)定立方體體元的邊長，劃分體元并為各體元8個頂點賦坐標(biāo)值。其中，體

在確定等值面與體元的相交狀態(tài)后，需要計算構(gòu)造的三角面片頂點的位置。當(dāng)三維離散數(shù)據(jù)場的密度較高，即體元很小時，可以假設(shè)函數(shù)值沿著體元的邊界呈線性變化。所以，等值面與體元邊的交點可以通過線性插值來求得，如（4）式所示。元的存儲結(jié)構(gòu)為：

struct gridcell

{

//存儲體元八個頂點坐標(biāo)的Vector3DF類型數(shù)組

Vector3DF v[8];

//存儲體元八個頂點法向量值的Vector3DF類型數(shù)組

Vector3DF n[8];

//存儲體元八個頂點密度值的double類型數(shù)組

double val[8];

};

第三步，在劃分體元后，計算體元各頂點的函數(shù)值，即密度值。密度值的計算可通過式（3）密度計算公式得到，具體的計算方法同海浪粒子的密度計算。

第四步，比較體元每個頂點的密度值與給定的等值面值，若頂點密度值大于給定閾值，賦值為0，否則賦值為1，這樣會構(gòu)造一個表示頂點與閾值之間關(guān)系的二進制值。比如：10111111表示頂點V6的值大于給定的等值面值。

第五步，查找構(gòu)造的狀態(tài)表，得到一個代表與體元邊界有交點的二進制值。例如第四步中根據(jù)頂點的二進制值查到狀態(tài)表得到的代表相交邊狀態(tài)的二進制值為：010001100000，它表示與邊5、6、10有交點。見立方體體元示意圖6。

圖6 立方體體元示意圖

第六步，通過線性插值方法，計算體元邊界與海浪密度等值面的交點坐標(biāo)；利用中心差分方法，計算體元各頂點處的法向，然后通過線性插值方法，求得所構(gòu)造的三角面片各頂點處的法向量。

最后，利用OpenGL中的glBegin(GL_TRIANGLES)三角面片繪制方法，根據(jù)插值得到的各三角面片的頂點坐標(biāo)值以及法向量繪制等值面圖像。

3 模擬效果對比

采用MC算法優(yōu)化的基于SPH方法的海浪流動模型實現(xiàn)方法綜合了SPH方法適合與求解高速碰撞等動態(tài)大變形問題的優(yōu)點，彌補了該方法粒子性強、連貫性差的缺點，相比單一的SPH方法構(gòu)建的海浪流動模擬效果連貫性更好，更符合流動的海浪的狀態(tài)，如圖7所示，SPH法較好地模擬了海浪碰撞效果，但是渲染效果粒子感較強，沒有形成連貫的水體模型；采用MC法優(yōu)化的方法既保留了海浪碰撞效果，又使水體模型更加連貫，更接近真實水體效果，優(yōu)化效果較為明顯。

圖7 SPH與MC優(yōu)化后的渲染效果對比

4 總結(jié)

結(jié)合海浪基本特征，利用物理模型與粒子系統(tǒng)相結(jié)合的方法—SPH方法完成了對海浪流動的模擬，然后用MC算法對海浪粒子模型進行了表面建模，使海浪模型更具有連貫性，效果更為逼真。將流體力學(xué)數(shù)值計算法與MC法結(jié)合起來完成了流動海浪的三維可視化模擬，是綜合利用體繪制與面繪制法實現(xiàn)三維可視化模擬的一次成功的探索。

存在問題：MC算法改進的基于SPH的海浪模擬并沒有脫離粒子系統(tǒng)，因此在場景較大，粒子數(shù)目增多后，模擬的實時性會較差。這可能需要考慮到GPU加速等技術(shù)。所以該模擬系統(tǒng)還有很多方面需要進一步研究與改進。

參考文獻：

[1]張靜.基于粒子系統(tǒng)的海浪動態(tài)演變模擬仿真研究[D].青島:山東科技大學(xué),2012.

[2]丁紹潔.虛擬海洋環(huán)境生成及場景特效研究[D].哈爾濱:哈爾濱工程大學(xué),2008.

[3]高峰.基于SPH的化學(xué)溶液傾倒過程仿真[D].長春:東北師范大學(xué),2010.

[4]鐘子春.實時流體交互性模擬算法的研究與實現(xiàn)[D].成都：電子科技大學(xué),2009

[5]陸慧蓮.基于LS-DYNA和HyperMesh的某型飛機垂尾前緣鳥撞分析[J].航空工程進展，2013,4(4):498-502.

[6]葉再春.MC算法研究及在三維流體可視化模擬中的應(yīng)用[D].蘇州:蘇州大學(xué),2009.

[7]矯春海.基于改進MC算法的醫(yī)學(xué)圖像三維重建研究[J].微型機與應(yīng)用,2011,30(3):39-45.

Optimization of the Wave Particle Mode Based on SPH by the MC Algorithm

LIU Wen-long1,ZHANG Jing2
1.First Crust Deformation Monitoring and Application Center,CEA,Tianjin 300180,China; 2.College of Geosciences and Technology,Shandong University of Science and Technology,Qingdao 266510,Shandong Province, China

Combined with the complex dynamics,spatial and temporal characteristics of the wave field and research of wave simulation methods,and in order to solve the problem of the lack of realistic rendering in wave simulation,and difficulty to simulate wave breaking and other issues,this paper realizes the simulation of dynamic evolution of the ocean waves by adopting the SPH method based on the physical characteristics of the fluid.In view of the discreteness and the poor realistic rendering in wave particle model based on the smoothed particle hydrodynamics method,this paper uses a free surface extraction algorithm based on the marching cubes algorithm,and completes surface modeling of the wave field,while optimizing the ocean wave effect based on particle system.The results can meet the requirements of simulating dynamic evolution of ocean waves with optimized MC algorithm,and the water model is more fluent and closer to the effects of real water bodies.The results proposed in this paper is a successful exploration of the comprehensive utilization of volume rendering and surface rendering method to realize 3D visual simulation.

wave simulation;SPH;marching cubes method

P731；TP391.9

A

1003-2029（2017）02-0041-05

10.3969/j.issn.1003-2029.2017.02.007

2016-05-19

劉文龍(1987-)男，碩士，主要研究方向為地理信息系統(tǒng)。E-mail：liuwenlon1987@126.com