基于流場的流動水體仿真

楊在興，芮小平，梁 明，賴 晗，陸 瑾

（1.西安科技大學 測繪科學與技術學院，陜西 西安710054；2.中國科學院研究生院 資源與環境學院，北京100049；3.中國礦業大學（北京）地球科學與測繪工程學院，北京100083；4.中國水利水電科學研究院 水環境研究所，北京100038）

0 引 言

大規模三維地形場景構建技術日趨成熟，但是大多數三維地形場景系統局限于表現靜態的地形和幾何數據，對三維地景系統中不可或缺的動態流動水體表現的較少。建設數字流域的一個重要方面就是利用數字流域水文地理信息平臺建立起水體仿真模擬系統，對全流域進行動態實時的三維仿真，為流域水污染防治、抗洪減災和大壩水庫選址等提供科學的決策支持。所以在三維地形場景上再疊加水體的仿真層，必然能使三維場景更加真實，對水文水利、水生態保護，水污染治理等領域都有非常重要的現實意義。

目前，國內對于水體仿真的研究集中在靜止水體水面的仿真，水面的效果多為波動效果，對自適應于河道邊界的自然流動水體仿真研究較少。王海玲等通過Gerstner模型生成水面模型的基本造型，利用改進的曲面熵算法對水面進行動態分區，對不同區域施加不同程度的水面擾動，通過修正紋理映射位移的方法模擬了動態水面，未對流動水面進行探討［1］。冶運濤等基于斷面控制點提取河流的邊界，應用動態紋理技術，通過實時更換系列紋理圖片，利用人眼視覺暫留效應，在虛擬仿真平臺與長河道一維水動力學模型集成的基礎上，實現了流動水體的仿真，此方法通過紋理切換實現水體的模擬，對于更大區域的水體模擬可能會面臨效率問題［2］。王道臣等在正弦波疊加的同時，利用兩個凹凸紋理實現了水面的模擬，也為對沿河道流動水面的模擬進行探討［3］。Qizhi Yu，Fabrice Neyret等研究了應用速度場驅動紋理來模擬流體，沒有對真實河流的仿真進行闡述［4］。作者在前人已有工作基礎上，重點研究并解決了自然流動水體的三維仿真問題。主要探討了應用流體力學原理的流場模型的建立，滿足泊松碟分布的塊狀紋理的渲染技術和流場驅動塊狀紋理沿河道的自適應移動等問題。在此基礎上，實現了自然河道內流動水體的仿真，并以三峽庫區為例，結合OpenSceneGraph（OSG）和GLSL著色語言建立了三峽庫區流域水體三維可視化系統。

1 流場的計算

1.1 基本原理

流場模型的建立是模擬水體的物理學基礎。流場是河道內所有點在同一時刻的速度矢量的分布狀態。如果能實時計算出河流的流場，就可以用速度場描述任意時刻水體的流動細節，在此基礎上就可以應用計算機和地理信息技術進行自然流動水體的仿真［5－6］。

本文采用的河道數據為包含河流邊界的水系網絡圖，水系網絡圖是一個有向圖，圖中包含各水系的空間分布、水系之間的連接狀況和所有水系的體積流量等信息。把這個圖映射到三維地形上，來確定河流在三維空間中的分布。

本文假設研究的流動水體為不可壓縮且流動模式為平面流動的流體，根據流體力學原理其流線的微分方程為－vdx＋udy＝0，在不可壓縮流體的平面流動中，速度場滿足不可壓縮流體的連續性方程：，必有一個函數滿足：dy＝－vdx＋udy＝0，于是函數（x，y）即為流場的流函數，u，v為流體的速度分布，在平面直角坐標系中為水流速度沿X軸和Y軸的分量［7－9］。水流速度可以用式（1）表示

由于在不可壓縮流體的平面流動中，通過兩條流線間任一曲線單位厚度的體積流量等于兩條流線的流函數之差，且與流線形狀無關。即Q＝Ψleft－Ψright。其中Q為通過兩條流線間任一曲線單位厚度的體積流量，Ψleft，Ψright分別為河道左右邊界上的流函數。

本文用到的水系網絡圖包含河流的體積流量信息。由于在同一條流線上各點的流函數為常數，在流線上dΨ＝0，即Ψ＝C（C為常數），也即每條流線對應一個常數值。如果給任意一個河流邊界流函數賦一個初值，就能夠算出其他邊界上的流函數值，然后就可以根據所有的邊界流函數對流函數進行插值，計算河道內的流函數。

1.2 對邊界流函數插值計算河道內的流函數

如圖1所示，假設要計算P點處的流函數，可以采用這樣的方法：規定一個搜索半徑s，找出經過以點P為圓心，s為半徑的圓的所有河道邊界，設為Bi，設di為點P到Bi的距離，設Ψi為河道邊界Bi上的流函數。應用線性插值法，可得

其中權重函數定義為

其中f（t）＝6t5－15t4＋10t3，f具有二階連續性，選此函數保證流函數Ψ在河道邊界處的連續性。通過式（2）以求出河道內任意點的流函數，然后就可以應用式（1）可求得同一時刻不同空間點的液體質點的速度分布，即河流流場。

圖1 流函數插值

2 流動水體的仿真

現階段對流動水體的模擬方法常采用構造條帶形三角網，動態更換紋理來實現［10－11］，這種方法隨著水體范圍的擴大，三角形數量明顯增加，程序整體效率就會受到影響，另外動態更換紋理不可避免地會使水面出現跳動。本文提出一種運用流場驅動塊狀紋理的方法來模擬水體的流動。所謂塊狀紋理，就是一個獨立的點攜帶一個紋理片所構成的紋理。應用塊狀紋理可大大降低三角形的渲染量，提高程序效率。具體方法是運用速度場驅動并約束攜帶塊狀紋理的點在河道內移動，通過河道邊界界定塊狀紋理是否活動，單個活動的紋理經過混合和渲染來重構水面，形成一種整體的自適應于河道的水體流動效果。

在可視范圍內，因為點攜帶的紋理塊有固定的大小，所以點的分布必須比較均勻。在均勻的前提下，點既不能太稀又不能太密，分布太稀，則河面會出現空洞，太密會占用太多的系統資源。本文點的分布采用了泊松碟分布法。該方法與傳統的單純隨機分布方法相比，既保持了采樣的隨機性，又使得采樣點的分布更加均勻，有效地調節了采樣間距，從而使點達到更好的分布效果。假設點間的最小距離為d，紋理塊的半徑為r（在屏幕坐標系中），為了防止流動的水面出現空洞，必須保證r大于或等于d。由于點是按照流場不斷移動的，為了保持泊松碟分布，還需要動態地插入和刪除點。由于河道形狀的不規則性，泊松碟分布采樣區域為整個屏幕可視區域，而不是可視區域內的河段，這樣采樣區域就是一個正方形，使采樣更簡單，速度更快。對于位于河道邊界之外多余的點，不參與移動、刪除和渲染，保證了系統的效率［12］。

2.1 均勻分布點的生成

前文提到應用泊松碟分布法來保證紋理均勻分布，由于傳統的泊松碟采樣算法十分復雜，很難實時地進行采樣計算，不能滿足實時生成的需要，為此，本文對泊松碟采樣方法進行了改進，提出了一種快速泊松碟采樣算法——邊界采樣算法。此算法既保持泊松分布的特性，且算法簡單，效率較高。如圖2所示，在多邊形的區域里隨機采樣生成一個泊松碟，泊松碟的半徑為r，只對可采樣的邊界進行采樣，而不是對剩下的整個區域進行采樣，圖中實線為可采樣的邊界，虛線為不可采樣的邊界，在可采樣的邊界上采樣生成泊松碟，泊松碟的邊界落到別的泊松碟中的部分，設為不可采樣的邊界，按照這樣的方法，直到完成整個區域的采樣［13］。

泊松碟不可采樣邊界的確定采用一個角度范圍，如圖3所示，假設現在要以點A的邊界為采樣邊界進行隨機采樣，點B為一個與A相鄰的點，點B相對點A的極坐標為（d，θ），則很容易得出點 A的不可采樣角度范圍為（θ－cos－1（d/2r），θ＋cos－1（d/2r））。可采樣角度范圍減去不可采樣的角度范圍就得出當前點A的可采樣邊界的角度范圍。

邊界采樣算法的具體步驟為：

（1）隨機采樣生成點集合P，無可采樣邊界點的集合為H；

（2）如果P為空，轉到（8），從P中隨機選擇一個點A，對邊界進行采樣，初始可采樣邊界角度范圍為angleRange；

（3）找出所有相鄰的點（兩點距離小于2r）構成集合Q；

（4）循環計算集合Q中的點相對點A的極坐標（d，θ）。計算不能采樣的角度范圍cutAngleRange。點A的可采樣角度范圍等于angleRange－cutAngleRange；

（5）如果angleRange等于0，把A從P中刪除，添加到H，返回到（2），否則，繼續下一步；

（6）在點A可采樣邊界上隨機采樣，生成泊松碟C，把C添加到集合P；

（7）計算點A的可采樣邊界，如果angleRange等于0，把A從P中刪除，添加到 H，回到（2），其他，回到（6）；

（8）返回集合H，結束采樣。

此算法的時間復雜度為O（N），能滿足實時采樣計算生成足夠多均勻分布點的需求。

2.2 點的驅動、刪除和插入

應用快速邊界采樣算法生成足夠點后，就可以對這些點進行操作了，通過式（1）可以計算河道內任一點的速度，設某一點的世界坐標為P，經過時間dt后的世界坐標為P＋V（P）dt，V是速度函數。通過實時計算點的速度，就能驅動點沿著河道水流方向移動。且點的移動自適應于自然的河道邊界。

點移動后，一些點可能已經移動到視景體之外，這些點已經不可見，需要把它從內存中刪除；可視范圍內的點分為兩種：活動點和不活動點。如果點移動后與其他點的距離太近，小于點之間最小距離d，為了維持泊松分布，需要將這些點設置為不活動的點，不活動的點也需要刪除，為了防止流動水體出現因突然刪除點而出現暫時的空洞，對不活動點采用漸漸淡出的渲染方式，然后將其從內存中刪除，渲染的混合因子為式（3）。

點移動或刪除后，就可能在剩余點之間產生空洞，所以需要根據剩余的點應用泊松碟分布算法不斷生成新的點，插入到已有點的中間，以填補空洞，且整體上所有點依然滿足泊松碟分布。采用2.1節的快速邊界采樣算法，應用該算法可以根據已有的點和點之間的最小距離d插入新的點，產生泊松碟分布。為了防止流動水體出現因突然插入點而產生局部的反差，插入點采用漸入的方式。混合因子為式（4）。

式中：t0——點生成時的時間，t1——點不再活動時的時間，T——漸入或淡出需要的時間。

2.3 點與紋理塊的綁定

如果給滿足泊松碟分布的每個點綁定一個紋理，紋理隨著點的進行移動、刪除和插入，點由流場驅動，紋理隨著點移動，然后對這些紋理塊進行混合，所有的紋理塊就能構成一個連續完整的流動水面。紋理塊的混合因子定義如下

其中

式（6）中x0是塊狀紋理的中心。

紋理塊在屏幕坐標系中的半徑為r，但是在世界坐標系中，塊狀紋理的大小是不一定的，與攜帶紋理的點到視點的距離有一定的比例關系。如果為每個紋理塊指定一個不同的紋理，可能會需要很多紋理，處理起來很麻煩，本文采用一塊紋理作為參考紋理，當生成新的點時，根據塊狀紋理的大小，從這一塊紋理中隨機選擇一部分作為紋理塊的紋理，構成塊狀紋理。

2.4 紋理塊的渲染

本文所采用的紋理為法線凹凸貼圖（NormalMap），如圖4所示。波動水面的渲染采用GPU可編程渲染管道，應用GLSL編寫GLSL著色器，對凹凸紋理進行渲染。渲染需要兩個著色器，頂點著色器和片元著色器。頂點著色器主要實現頂點的坐標變換，計算視線向量、光源位置、頂點反射向量、折射向量和紋理坐標；片元著色器主要實現紋理的采樣，計算光照強度，返回片元的最終顏色［14－15］。對于水面，觀察者與水面的角度越小，反射效果越明顯；角度越大，折射效果越明顯，物理學上稱為菲涅爾效應。所以在渲染水面時要根據觀察者的角度來計算反射、折射顏色以及水面顏色的混合方式。物理上菲涅爾計算比較復雜，實際應用常采用以下的計算方法

式中：F——空氣到水面的反射系數，值為0.020037，eye－Dir——視線方向，normal——光線與水面交點處的法線向量。最終得到的水面的顏色為

式中：finalColor——最終歸水的顏色；waterColor——初始水的頂點顏色；reflectionColor——水面的反射顏色；refractionColor——水面的折射顏色。

頂點著色器具體實現步驟：

（1）計算頂點經過變換后在世界坐標系中的位置；

（2）計算視線方向向量，光源位置，并進行向量歸一化；

（3）計算頂點的反射和折射向量，并進行向量的歸一化；

（4）計算紋理坐標。

片段著色器具體實現步驟：

（1）計算反射方向和折射方向的紋理采樣值；

（2）計算環境光、散射光、鏡面光；

（3）返回片元的最終顏色。

圖4 所用紋理

3 試驗結果與分析

根據本文提出的方法，作者結合OSG和GLSL著色語言構建了一個流域水體三維仿真系統。為了系統的擴展性和推廣應用前景，作者采用的試驗環境為：Windows XP SP3系統，P4CPU，1G內存，ATI顯卡，顯存256MB。用戶可以流暢的操作系統在場景內漫游，瀏覽水體的流動效果，水面細節（水流的緩急）依賴于河道形狀和方向，并具有較好的光照效果，流動水體的仿真效果如圖5所示。

圖5 流動水體可視化效果

4 結束語

水體仿真是實現數字流域的關鍵，作者根據水體的特點，提出了一套流動水體仿真的方法，它可以概括為以下幾個方面：①根據水體流動的特點，提出用流場驅動塊狀紋理運動來實現河道自適應的水體流動效果。②根據水體的邊界流函數插值得到河道內任意點的流函數。利用流體力學原理，給出了河道內任意點的速度的求解方法。③根據本文研究的需要，提出了一種快速泊松碟邊界采樣算法，此方法滿足泊松分布的優良特性，實現起來相對簡單，且具有線性時間復雜度。另外由于本文的采樣區域只是可視區域，這樣點的數量就不會很多，提高了實時采樣效率。快速泊松碟邊界采樣算法能使點攜帶的塊狀紋理分布均勻，避免多余的渲染。④采用GLSL著色語言編寫著色器，利用GPU可編程渲染管道來處理數據，對移動的塊狀紋理進行混合和渲染，使用可編程著色管線來處理數據，減少CPU的實時運算量，將CPU從復雜龐大的運算中解放出來，利用GPU進行圖形的計算，使CPU和GPU之間更協調地工作，保證渲染的質量和效率。

實踐表明，本文提出的方法對流動水體仿真具有很好的借鑒意義。由于采用流場驅動紋理來實現河流的模擬，所以借助流場可以很容易的結合具體的水體污染擴散模型對河流污染物的擴散進行模擬仿真。本文研究的河道類型為中間沒有島的河道，對于中間有島的復雜河道，還需要進一步探討。

［1］WANG Hailing，YIN Guisheng，ZHANG Jing，et al.Dynamic water surface simulation method based on improved surface entropy ［J］.Computer Engineering，2011，37（6）：24－26（in Chinese）.［王海玲，印桂生，張菁，等.基于改進曲面熵的動態水面模擬方法 ［J］.計算機工程，2011，37（6）：24－26.］

［2］YE Yuntao，ZHANG Shanghong，WANG Xingkui.Research on 3Dvisualization of flood propagationin three gorge reservoir area ［J］.Journal of System Simulation，2009，21（14）：4379－4382（in Chinese）.［冶運濤，張尚弘，王興奎.三峽庫區洪水演進三維可視化仿真研究 ［J］.系統仿真學報，2009，21（14）：4379－4382.］

［3］WANG Daochen，WAN Wanggen，TANG Jingjou，et al.Realtime rendering algorithm for water surface based on GPU ［J］.Computer Engineering，2008，34（20）：233－234（in Chinese）.［王道臣，萬旺根，唐經洲，等.基于GPU的水面實時渲染算法 ［J］.計算機工程，2008，34（20）：233－234.］

［4］INRIA.Spectrum－preserving texture advection for animated fluids：Report of INRIA ［R］.YU Qizhi，Fabrice Neyret，Eric Bruneton：INRIA，2009.

［5］WU Dijun，SUN Haiyan，HUANG Quanyi，et al.Numerical simulation of flood routing in watercourse and its visualization on GIS platform ［J］.Science of Surveying and Mapping，2009，34（6）：242－243（in Chinese）.［吳迪軍，孫海燕，黃全義，等.河道洪水演進數值模擬及其在GIS平臺上的可視化 ［J］.測繪科學，2009，34（6）：242－243.］

［6］LI Yun，FAN Ziwu，WU Shiqiang，et al.Numerical simulation and 3Dvisualization of flood propagation in large scale detention basins ［J］.Journal of Hydraulic Engineering，2005，36（10）：1158－1164（in Chinese）.［李云，范子武，吳時強，等.大型行蓄洪區洪水演進數值模擬與三維可視化技術 ［J］.水利學報，2005，36（10）：1158－1164.］

［7］ZHANG Dawei，CHENG Xiaotao，HUANG Jinchi.2Dnumerical modeling of levee－breach floods in dense urban areas ［J］.Journal of Hydraulic Engineering，2010，41（3）：272－276（in Chinese）.［張大偉，程曉陶，黃金池.建筑物密集城區潰堤水流二維數值模擬 ［J］.水利學報，2010，41（3）：272－276.］

［8］WANG Zhili，GENG Yanfen，JIN Sheng.Numerical modeling of 2Dshallow water flow with complicated geometry and topography［J］.Journal of Hydraulic Engineering，2005，36（4）：439－443（in Chinese）.［王志力，耿艷芬，金生.具有復雜計算域和地形的二維 淺 水 流 動 數 值 模 擬 ［J］. 水 利 學 報，2005，36（4）：439－443.］

［9］XIE Zuotao，ZHANG Xiaofeng，TAN Guangming，et al.Study and application of mathematical model for one－dimensional floodrouting ［J］.Engineering Journal of Wuhan University，2005，38（1）：69－72（in Chinese）.［謝作濤，張小峰，談廣鳴，等.一維洪水演進數學模型研究及應用 ［J］.武漢大學學報（工學版），2005，38（1）：69－72.］

［10］HUANG Zhijun，WANG Shuwen.Three－dimensional visualization simulation of flood routing of city based on hydrological mathematical model［J］.Science of Surveying and Mapping，2010，35（1）：88－90（in Chinese）.［黃志軍，王樹文.基于水文仿真模型與GIS的三維可視化在城市洪水演進中的模擬研究——以馬鞍山為例 ［J］.測繪科學，2010，35（1）：88－90.］

［11］LIU Quan，HU Zhigen，HUAI Wenxin，et al.Three－dimen－sional dynamic simulation of earth－rock cofferdam failure flood［J］.Journal of Hydraulic Engineering，2010，41（12）：1432－1435（in Chinese）.［劉全，胡志根，槐文信，等.土石圍堰潰堰洪水三維動態仿真模型 ［J］.水利學報，2010，41（12）：1432－1435.］

［12］DONG Bingjiang，ZHANG Xiaofeng，QIAO Wei.Real－time interactive simulation shows for dam－break flood based on OpenGL technology ［J］.Engineering Journal of Wuhan University，2009，42（2）：168－171（in Chinese）. ［董炳江，張小峰，喬偉.基于OpenGL的潰壩洪水實時交互仿真顯示［J］.武漢大學學報（工學版），2009，42（2）：168－171.］

［13］Dunbar D，Humphreys G.A spatial data structure for fast poisson－disk sample generation ［J］.ACM Transactions on Graphics 2006，25（3）：503－508.

［14］ZHOU Zhili.Theory and applications of wave［M］.Beijing：Science Press，2005（in Chinese）.［鄒志利.水波理論及其應用 ［M］.北京：科學出版社，2005.］

［15］HAN Min，CHEN Ming.Three－dimension visualization research of flood routing simulation ［J］.Computer Applications，2005，25（8）：1906－1907（in Chinese）.［韓敏，陳明.洪水演進模擬的三維可視化研究 ［J］.計算機應用，2005，25（8）：1906－1907.］