周梟楠，周貴云

（1.電子科技大學(xué)，四川 成都611731）

基于DEM數(shù)據(jù)的SPH水體流動模擬

周梟楠1，周貴云1

（1.電子科技大學(xué)，四川 成都611731）

對基于DEM數(shù)據(jù)的SPH水體流動模擬技術(shù)進行研究。對整個模擬流程進行重點介紹，并根據(jù)DEM數(shù)據(jù)的特點，對SPH的邊界條件處理方法進行修改。最后借助OpenGL的顯示能力，在2種不同類型的地形上對水體流動進行實時計算和模擬。實驗表明，對邊界條件處理方法的修改是合理的，基于DEM數(shù)據(jù)的SPH水體流動模擬結(jié)果令人滿意。

光滑粒子流體動力學(xué)；數(shù)字高程圖像；OpenGL；邊界條件

1 基本概念介紹

DEM是一種用X、Y、Z坐標描述地表形態(tài)的數(shù)據(jù)模型。DEM具有數(shù)據(jù)精度高，易于獲取和便于分析等一系列優(yōu)點，被廣泛應(yīng)用在測繪、地質(zhì)、地理、軍事等領(lǐng)域[1-4]。利用DEM數(shù)據(jù)不但可以讓地形場景可視化，還可以配合對應(yīng)區(qū)域的遙感影像進行紋理貼圖，從而使得地形場景更為逼真。

SPH（smoothed particle hydrodynamic）即光滑粒子流體動力學(xué)，是一種用來模擬水體流動的拉格朗日無網(wǎng)格粒子法，最初由Gingold和Monaghan[5]、Lucy[6]于1977年分別提出。SPH算法的原理是將連續(xù)系統(tǒng)離散化為一堆相互作用的單獨粒子，同時給這些粒子賦予諸如質(zhì)量、速度、粘度等各種物理特性，而這樣做就使得這些粒子具有真實粒子的材料特性。這樣一來，研究者只需通過對這一堆相互作用的粒子集合的運動狀況、相關(guān)特性進行觀察和研究，便可以得到整個系統(tǒng)的運動狀態(tài)和相關(guān)特性。

2 模擬流程

本文是將DEM數(shù)據(jù)與SPH算法相結(jié)合，來實現(xiàn)基于DEM數(shù)據(jù)的SPH水體流動模擬。圖1為整個模擬過程流程圖。

2.1 數(shù)據(jù)預(yù)處理

從“國際科學(xué)數(shù)據(jù)服務(wù)平臺”下載選定區(qū)域的DEM數(shù)據(jù)和與對應(yīng)區(qū)域的LandSat7第3、4、5波段數(shù)據(jù)。由于DEM數(shù)據(jù)所表示的區(qū)域與LandSat7數(shù)據(jù)所表示的區(qū)域并不一定完全重合，故需要對數(shù)據(jù)進行裁剪和拼接，保證2種數(shù)據(jù)所表示的區(qū)域重合。數(shù)據(jù)處理完成后，將處理后的DEM數(shù)據(jù)以ASCII的格式進行保存，而將融合的LandSat7數(shù)據(jù)以BMP格式進行保存。

圖1 水體流動模擬流程圖

2.2 數(shù)據(jù)導(dǎo)入

這一步需要將保存后的DEM數(shù)據(jù)和LandSat7數(shù)據(jù)導(dǎo)入系統(tǒng)內(nèi)存中，供后面進行相關(guān)計算和結(jié)果輸出顯示。本文以O(shè)penGL為技術(shù)基礎(chǔ)來構(gòu)建整個模擬系統(tǒng)，它本身只提供了點、線、多邊形等簡單圖像的繪制，并未直接提供對ASCII格式文件的讀取和顯示功能。因此，需要依靠人工編寫程序?qū)ζ溥M行讀取和繪制?？紤]到DEM數(shù)據(jù)規(guī)則網(wǎng)格的特點，同時方便后面進行紋理貼圖和與SPH法相結(jié)合進行水體流動模擬，故決定將DEM數(shù)據(jù)以圖2的方式進行保存。

圖2 DEM數(shù)據(jù)保存示意圖

考慮到DEM數(shù)據(jù)的規(guī)則性，本文以一個小網(wǎng)格為單位，從左上角開始，從左到右，從上到下，依次讀取并保存點坐標。當讀取到DEM數(shù)據(jù)中任意一個小網(wǎng)格，在存儲的時候，將其按圖2中的樣子，分割成2個小三角形，并以逆時針的順序存儲這2個三角形中的頂點坐標，同時也將該三角形所對應(yīng)的編號進行存儲。而同時，也將該點所對應(yīng)的紋理坐標進行存儲。頂點的紋理坐標計算公式為（i/rows，j/cols），其中i，j分別為該頂點對應(yīng)的行列號，而rows和cols為整個紋理貼圖的行列總數(shù)。得到任意頂點的紋理坐標后，因為在第一步已經(jīng)將DEM數(shù)據(jù)和紋理數(shù)據(jù)進行了匹配，而通過這個紋理坐標就可以在紋理數(shù)據(jù)中找到對應(yīng)點的紋理，從而對DEM地形進行正確的紋理貼圖。

2.3 水體流動模擬計算

結(jié)合DEM數(shù)據(jù)自身的特點，本文使用一種簡潔的方法來處理SPH算法中的邊界問題。

2.3.1 SPH公式

前文提到，SPH算法是將連續(xù)系統(tǒng)離散為一堆相互作用的粒子集合，通過研究每一個粒子的運動情況，最終來模擬整個系統(tǒng)的運動情況。因此，如何得到每一個粒子的運動狀態(tài)就顯得尤為重要。本文采用了文獻[7]中提到的SPH計算方法，對于任意粒子，其加速度計算公式如下：

式中，g?為重力加速度，在通常情況下為常數(shù)；分別為水流粒子受到的壓力和速度；h為粒子的光滑核半徑；ρi為粒子的密度；m為粒子的質(zhì)量；粒子的密度、質(zhì)量和重力加速度均為常數(shù)。根據(jù)系統(tǒng)的連續(xù)性，在初始條件已知的情況下，根據(jù)上式，可以得到粒子在任意時刻的加速度當?shù)玫郊铀俣群?，根?jù)蛙跳積分法，便可以得到粒子所在的位置。當?shù)玫秸麄€系統(tǒng)任意時刻所有粒子的運動狀態(tài)，也就得到了整個系統(tǒng)在該時刻的運動狀態(tài)。

2.3.2 邊界條件處理

為了讓SPH算法能夠正確地模擬水體在真實地形上的流動，且不會產(chǎn)生粒子物理性穿越邊界或是與真實運動情況不符合的現(xiàn)象，根據(jù)DEM數(shù)據(jù)的特點，本文采用了邊界碰撞法來處理邊界問題，如圖3。

判斷某一點與邊界的關(guān)系時，考慮到DEM數(shù)據(jù)的特點，即在知道整個地形網(wǎng)格的起始坐標和矩形網(wǎng)格寬度的情況下，可以很容易得到該點落在哪一個矩形網(wǎng)格里面。如圖3中左圖所示。將該網(wǎng)格與該點同時向z軸方向進行投影，只需要判斷投影后的點落在哪一個投影三角形內(nèi)，即可得到該點與哪一個三角網(wǎng)格最接近。判斷該粒子到這個三角網(wǎng)格的垂直距離 ，如圖3右圖所示。當該點到三角網(wǎng)格的垂直距離小于設(shè)定值時，即認為該點過于接近邊界，需要讓其進行反彈。此時，根據(jù)預(yù)先設(shè)定好的參數(shù)來修正該粒子的加速度，從而保證其無法穿越邊界。2.4 結(jié)果輸出顯示

圖 3 邊界處理示意圖

使用上面方法得到系統(tǒng)在任意時刻任意點的運動狀態(tài)后，只需要進行簡單的循環(huán)顯示，便可以將該系統(tǒng)在任意時刻的運動狀態(tài)展現(xiàn)出來。本文采用OpenGL進行三維顯示。整個過程均采用實時計算和實時演示。

3 水流推進演示

為了觀察SPH方法在DEM數(shù)據(jù)地形上的流動效果，同時驗證前文提出的邊界處理方法是否正確，本文采用2種不同的地形：人工生成地形和DEM數(shù)據(jù)生成地形，在這2種地形上分別模擬水體流動。整個模擬過程借助OpenGL進行實時計算和實時演示。

本文實驗的硬件平臺配置如下：中央處理器為 I3-2120；內(nèi)存大小為4 G；顯卡為AMD 7400。軟件平臺為如下：WIN7操作系統(tǒng)；VS2010 編譯工具；OpenGL版本為4．1。

實驗1：采用人工生成的地形。首先將用于模擬水體流動的SPH粒子群放在一個長方形的水槽左邊，在水槽的中心設(shè)置一個矩形障礙物。當粒子群開始流動后，粒子群沿著水槽向右側(cè)開始流動。從圖4中左圖可以看出，當粒子遇到障礙物后，并未出現(xiàn)粒子穿越邊界的情況。同時粒子會正確地繞開障礙物并向前推進，與實際情況相符合。測試中粒子總數(shù)為4 114個，刷新率為45 fps左右。可見在小范圍場景下能得到令人滿意的顯示效果。

圖4 實驗效果圖

實驗2：采用DEM數(shù)據(jù)生成的地形。首先從網(wǎng)上下載DEM數(shù)據(jù)，并將其裁剪成256×256的網(wǎng)格,然后下載對應(yīng)區(qū)域LandSat7的3、4、5波段數(shù)據(jù)，然后融合、裁剪并導(dǎo)出。

搭建好地形后，便可以將用于模擬水體的粒子群布置到地形上進行水體流動模擬。最終結(jié)果如圖5右。從圖上來看，在SPH算法的驅(qū)動下代表水體的深藍色粒子群沿著山谷向前推進，模擬效果正確。測試中粒子的總數(shù)為10 648個，刷新率在20 fps左右。雖然顯示效果不是特別流暢，但依然在可以接受的范圍之內(nèi)。

[1] 梁欣,安如,張發(fā)明,等． 庫岸滑坡地質(zhì)災(zāi)害三維演變動態(tài)顯示方法[J]．地理空間信息, 2013,11(2):10-14

[2] 張姝慧．求解淺水波方程的光滑粒子流體動力學(xué)[D]．合肥:安徽大學(xué),2007

[3] 林昊． 基于光滑粒子流體動力學(xué)的OpenGL可視化[D]． 合肥:安徽大學(xué)，2011

[4] 趙慶祝,張清,寧川．基于OpenGL的DEM地形可視化和虛擬漫游系統(tǒng)[J]．計算機系統(tǒng)應(yīng)用,2006(5):66-69

[5] Gingold R A, Monaghan J J． Smoothed Particle Hydrodynamics: Theory and Application to Non-spherical Stars[J]． Monthly Notices of the Royal Astronomical Society, 1977(181): 375-389

[6] Lucy L B． Numerical Approach to Testing the Fission Hypothesis[J]．Astronomical Journal, 1977(82): 1 031-1 027

[7] SPH算法簡介（二）． 粒子受力分析[EB/OL] ．Http://www．starming．com/index．php?action=plugin&v=wave&tpl=union&a c=viewgrouppost&gid=34694&tid=12522，2011-03-31

P208

B

1672-4623（2015）02-0086-02

10.3969/j.issn.1672-4623.2015.02.032

周梟楠，碩士，主要研究方向為地理信息系統(tǒng)理論及應(yīng)用。

2014-02-12。

項目來源：國家自然科學(xué)基金資助項目（41371399)。