劉鵬程， 黃 欣*， 郭明武

(1.華中師范大學地理過程分析與模擬湖北省重點實驗室， 武漢 430079；2.華中師范大學城市與環(huán)境科學學院， 武漢 430079； 3.武漢市測繪研究院， 武漢 430022)

“粒子系統(tǒng)”是在1983年由William T Reeves初次提出的[1]，用以描述自然界中一些運動不存在規(guī)律性的事物，目前已成功模擬了爆炸、焰火等動態(tài)運動效果[2].粒子系統(tǒng)的基本原理是將模擬的事物定義為由無數(shù)個隨機產生且在其有限的生命周期內做不規(guī)則運動的粒子構成[3].每個粒子需經歷“產生”“運動”以及“消失”三個階段，粒子運動到下一位置生命周期減1，每個時刻均有新粒子的產生和舊粒子的消失.通過在某一時刻 “運動”粒子的速度及方向等信息可以展示事物的整體特征和動態(tài)變化.粒子系統(tǒng)較為充分地呈現(xiàn)事物運動的隨機性與動態(tài)性，逼真地模擬事物的變化[2].

在模擬不規(guī)則運動物體的研究領域，粒子系統(tǒng)是學術界公認的最成功的圖像模擬生成算法[4]，具備重大的研究意義與廣闊的發(fā)展前景.國內外的專家學者使用粒子系統(tǒng)的方法成功模擬出了很多的自然事物，并不斷地對模擬方法進行改進，以追求更高的真實和實時性.

1985年，Reeves和Blau[5]提出“volume filling”單元用以模擬隨著時間的變化，形狀發(fā)生改變但整體基本不變的事物，如隨風飄動的花草樹葉.1990年，Sims[6]首次實現(xiàn)了使用計算機成功繪制粒子并呈現(xiàn)出粒子運動的動畫效果，提出了粒子系統(tǒng)的并行繪制算法.該算法實現(xiàn)了電影“Particle Dream”中火焰、爆炸、瀑布等影視效果.Richard和David[7]使用粒子系統(tǒng)對彈性物體表面進行建模.此外，Unbescheiden和Trembilski[8]從云的物理原理出發(fā)，將粒子系統(tǒng)與紋理映射技術相結合，建立了云的模型.萬華根等[9]基于粒子系統(tǒng)和流體動力學模擬出了噴泉中水的實時流動效果.謝劍斌等[10]通過粒子系統(tǒng)模擬了雨點和雪花降落過程.Kolb等[11]采用 GPU 進行粒子系統(tǒng)的模擬，使得大規(guī)模進行粒子的模擬變成了可能.付帥等[12]利用粒子系統(tǒng)對全球表層流場進行仿真表達.Reynolds[13]使用粒子系統(tǒng)模擬了動物群體的運動.此外，Goss[14]通過粒子系統(tǒng)實時地模擬了船舶行駛時的船跡.

數(shù)字高程模型(digital elevation model ，DEM)是格網數(shù)據(jù)點的二維平面坐標(X，Y)和高程值(Z)的數(shù)據(jù)集[15]，用來表達地表的起伏和特征[16].目前，DEM的可視化表達方法主要在一維剖面、二維平面以及三維立體等空間維度展開，使用如等高線法[17]、陰影法[18]、明暗等高線法[19]、地形剖面線法[20]等靜態(tài)表達方法表達地形地貌.此類方法一般使用一幅圖形或圖像來對地形起伏狀況進行表達，其信息的獲取與規(guī)律的挖掘主要取決于用圖者分析的能力.

本文嘗試利用粒子系統(tǒng)模擬水體在地表上的運動來形象描述地表的高低起伏狀況，其中粒子的流向和流速可以用以表示地表單元的坡向和坡度.這為地形的表達提供了一種新的動態(tài)可視化方法，能夠很好地將局部和全局的地形地貌特征展現(xiàn)在觀察者面前，并實現(xiàn)交互式、動態(tài)性的地形瀏覽.

流模式是指利用粒子系統(tǒng)的方法構建水流模型，用以模擬水流在地表上的流動進而實現(xiàn)對地形起伏情況的動態(tài)可視化表達.其中，粒子系統(tǒng)將水流視為是無數(shù)隨機分布的粒子集合，每個粒子都有其固有的屬性，如初始位置、流動的方向和速度、顏色、生命周期等.每個粒子都能在水流場中根據(jù)其初始狀態(tài)不斷地自主進行不規(guī)則運動，模擬出接近真實流動形態(tài)的動態(tài)水流，實現(xiàn)對DEM的動態(tài)可視化表達.基于粒子系統(tǒng)的水流動態(tài)模擬主要流程包括以下幾步.

1) 水流場數(shù)據(jù)轉化為二維數(shù)組數(shù)據(jù).由于粒子的生成位置是隨機且大量的，所以為了便于快速地計算水流場中任何一個粒子在當前格網的水流方向以及速度值，需要將水流場數(shù)據(jù)轉化為二維數(shù)組數(shù)據(jù).二維數(shù)組中的每一項均存儲當前格網的橫向水流速度(u)和縱向水流速度(v).

2) 生成若干個粒子，并計算粒子的運動位置.二維數(shù)組數(shù)據(jù)相當于基礎數(shù)據(jù)，數(shù)據(jù)量有限，但是在大多數(shù)情況下，粒子的數(shù)量要遠大于該數(shù)據(jù)量.所以想要更為詳細地展現(xiàn)研究區(qū)范圍的水流情況，即地勢的高低起伏，就需要在網格上創(chuàng)建較多的粒子來進行動態(tài)可視化顯示.本文主要采用兩種方法來確定粒子的位置，一種是利用Math.random()函數(shù)在水流場范圍內生成大量粒子，使粒子能夠在整個研究區(qū)域內隨機分布.另一種方法是通過查看粒子的流向走勢可以獲取地形的特征，如山脊線、山谷線等，分析總結粒子在各類地形中的流動規(guī)律對地形特征進行提取，進而實現(xiàn)在重要特征地形上粒子的自動放置.這樣既能反映研究區(qū)的總體地形特征，又能逼近局部關鍵地形.

3) 計算粒子流動位置并進行繪制.多次調用animate方法，用于計算每個粒子的當前位置、水流的速度及方向、下一次運動要抵達的位置.每次使用animate方法，粒子的生命周期參數(shù)就會減少1，直到其值變化為0.當粒子生命周期為0或者粒子移動到水流速度為0處，會使得粒子消失進而導致水流場中的水流粒子越來越少.為解決此類問題，需對粒子生命周期進行重新賦值，每個時刻均有粒子的產生和消失，確保水流場內一直有水流在流動.

在Canvas上進行當前運動粒子所構成水流的繪制，實現(xiàn)水流的動態(tài)可視化.選擇合適的地理底圖，在相同的坐標系統(tǒng)下實現(xiàn)粒子與底圖的位置匹配.每次進行地圖的縮放操作時，計算出粒子的數(shù)量及運動狀態(tài)，根據(jù)當前地圖所在范圍進行水流的重新繪制.同時，還能實現(xiàn)在不同大小的窗口中自適應顯示粒子的流動.

4) 將重要位置標注在底圖上.將基于DEM數(shù)據(jù)提取得到研究區(qū)主要河流、湖泊、匯流口等特征位置的矢量數(shù)據(jù)轉化為GeoJSON格式，加載到地圖底圖上進行表示，可以增強地形起伏的可視化效果，提高美觀性.

1水流場數(shù)據(jù)的生成

粒子的流動是非周期性的，水流的強弱和方向可以使用流速和流向來表示，故可將水流數(shù)據(jù)定義為向量場數(shù)據(jù)，在向量場中每個點對應的數(shù)據(jù)均為一個向量，某點的向量方向代表水流在該位置的流向，向量大小代表流速.在高低起伏的地面單元上粒子傾向于運動，其運動的驅動力源自地表起伏所形成的各個坡度，且驅動力的方向與坡向一致.將驅動力分解為u和v兩個正交方向，構建二維水流場數(shù)據(jù).其具體步驟如下.

1.1構建水流模型

水流場數(shù)據(jù)的生成需要考慮到地面坡度對水流大小及方向的影響.使用D8算法[21](deterministic eight-node)來構建水流模型.D8算法假設某中心柵格單元在3×3窗口中有8種水流的可能方向，即其所在位置的八鄰域方向.使用公式(1)分別計算中心柵格單元與八鄰域的格網單元的高程差與距離的比值，即坡度.最大坡度值所對應的方向為中心柵格單元的水流方向.

(1)

在計算中心網格單元與其周圍八個相鄰柵格單元的高程差時，存在兩種情況：若高程差均小于或等于0，則該中心單元無水流流出，粒子速度為0；否則，中心單元的水流流向坡度最大值所對應的網格單元，粒子的初速度值為坡度的倍數(shù).

圖1展現(xiàn)了對DEM數(shù)據(jù)進行逐網格遍歷，計算粒子流動方向的過程：在3×3窗口中計算得到中心網格單元的最大坡度值，并確定水流流動方向，為了方便后序粒子分量的內插求解，將柵格數(shù)據(jù)轉化為具有方向屬性的點數(shù)據(jù).

圖1 水流方向的計算過程Fig.1 The process of calculating flow directions

1.2計算水流速度

各個網格單元的水流速度與最大坡度成正比，設置高差縮放因子對單元的水流速度進行調整從而優(yōu)化粒子流的效果.如圖2所示，以中心柵格點為原點構建uov直角坐標系，坐標系的象限與水流可能流動方向的八鄰域相對應.假定正東(E)方向為u軸的正方向，正北(N)方向為v軸的正方向.使用公式(2)、(3)將中心單元的水流向量速度分解為u、v分量，分別記錄水流粒子的速度與方向.

圖2 uov坐標系Fig.2 uov coordinate system

u=k×slope×sinθ，

(2)

v=k×slope×cosθ，

(3)

式中，k為縮放因子，用以調整研究區(qū)高程差；slope為各中心網格單元在八鄰域方向坡度的最大值.角度θ為水流方向與v軸的夾角，θ=n×45°(n=0，1，…，7)，0、1、2、3、4、5、6、7分別對應中心網格單元在N、NE、E、SE、S、SW、W及NW的鄰域方向.

1.3生成水流場數(shù)據(jù)

將基于DEM數(shù)據(jù)進行逐網格計算得到的水流粒子速度u、v分量，與研究區(qū)地理坐標范圍、網格的行列數(shù)和空間分辨率等信息相結合，生成JSON格式的水流場數(shù)據(jù).圖3展現(xiàn)了JSON文件的基本格式，主要包括兩個數(shù)據(jù)塊，分別用于存儲u、v分量的相關數(shù)據(jù)，每個數(shù)據(jù)塊均包含{“header”：{ }，“data”：[ ]}的結構.

“header”結構存儲一些關于水流場的屬性參數(shù)信息，一般JSON數(shù)據(jù)中兩個數(shù)據(jù)塊的“header”內容基本一致.“data”結構存儲水流速度的數(shù)據(jù)，每個值對應研究范圍內的一個網格點.第一個“data”存儲的是各網格點的橫向水流速度，即u分量，第二個“data”存儲的是各網格點的縱向水流速度，即v分量.兩個“data”中的值對應進行向量運算，即可計算得到該點的水流向量值.

2地形起伏動態(tài)模擬

2.1雙線性插值法計算水流分量

在粒子系統(tǒng)中，水流是由無數(shù)個隨機分布的粒子所構成的集合，每個粒子都能在水流場中根據(jù)其初始狀態(tài)不斷地自主進行運動.粒子系統(tǒng)能夠模擬出接近真實流動形態(tài)的動態(tài)水流，其驅動力是地面的勢能場，人們可以通過粒子流的運動來感知地形的起伏，進而實現(xiàn)對DEM的動態(tài)可視化表達.

由于水流場數(shù)據(jù)文件記錄的是中心網格點的水流速度，而隨機生成的粒子，其空間位置不可能完全位于中心網格點上.故粒子所在位置的速度需要依據(jù)粒子相鄰網格點的速度值進行插值運算得到.本文使用雙線性插值算法(Bilinear Interpolation)，來獲取該粒子u、v分量的近似值，并計算出該水流粒子下一次運動所抵達的位置.

圖3 JSON數(shù)據(jù)格式Fig.3 JSON data format

以插值計算u分量為例，圖4中點p1=(x1y1)，p2=(x2y1)，p3=(x2y2)，p4=(x1y2)為待插值粒子p=(xy)相鄰的四個網格點，五個點的平面坐標及四個相鄰網格點對應的u分量值均已知.使用粒子p所在位置的周圍4個相鄰點對應的分量值進行加權插值估計p點對應的u分量值.

圖4 雙線性插值示意圖Fig.4 The schematic map of bilinear interpolation

首先，使用公式(4)及(5)在x方向上進行線性插值，計算Q1、Q2點所對應的u分量值：

(4)

(5)

其中，Q1和Q2的平面坐標分別為(xy1)及(xy2).

再基于計算得到的Q1、Q2點所對應的u分量值，使用公式(6)在y方向上進行線性插值，即可得到粒子p所在位置的u分量值.

(6)

為了簡化計算，每次計算粒子對應速度分量時，將p1、p2、p3、p4坐標設置為(0，0)、(1，0)、(0，1)及(1，1)，并相應調整p點的坐標.使用上述相同的方法，對p點對應的v分量進行計算，即可獲得粒子所在位置對應的u、v分量值.

根據(jù)粒子當前所在位置和計算得到的u、v分量，即可計算得到粒子下一次運動到的位置，依次計算下去.每計算1次，粒子的生命周期參數(shù)就會減少1，直到其值變化為0.當粒子生命周期為0或者粒子移動到水流速度為0處，粒子將消失.

2.2粒子位置布設及屬性說明

粒子初始位置的布設使用兩種方法，其一使用隨機函數(shù)在數(shù)據(jù)范圍內隨機生成，粒子的數(shù)量根據(jù)水流場數(shù)據(jù)的分布規(guī)律而變化，也可以通過調節(jié)參數(shù)來進行控制.其二通過查看粒子的流向走勢可以獲取地形的特征，實現(xiàn)在一些重要地形特征點或者特征線(如山脊線和山谷線等)粒子的自動放置，專門針對這些特征點進行水流的動態(tài)模擬，充分展現(xiàn)局部地貌骨架線.

隨機生成的粒子屬性需要進行初始化，粒子屬性包括粒子的初始位置、流動的方向和速度、顏色、樣式及生命周期等，表1中列出了一些粒子屬性參數(shù)和相關描述.設置不同的粒子屬性參數(shù)值，會產生不同的動力學效應，水流粒子的運動也會呈現(xiàn)出不同的動態(tài)效果.在實際的應用中，可以通過調整粒子流動速度的比例因子即公式(4)、(5)中的k值，對區(qū)域整體的粒子流動速度進行夸張，使粒子流動的整體效果能夠更好地反映地表起伏和特征.

表1 水流粒子的屬性參數(shù)說明Tab.1 The parameter description of water flow particle attributes

2.3粒子流的多尺度表達

通過設置水流場的網格大小、粒子的密度以及在地形特征線上基于不同尺度放置粒子，可以實現(xiàn)不同空間分辨率地形粒子流的動態(tài)可視化表達.當對地圖進行縮放操作時，平臺將根據(jù)當前屏幕窗口范圍進行計算，重新繪制粒子并自動切換與當前比例相對應的底圖.

為了獲得更好的地形起伏效果，有時會在底圖上覆蓋不同空間分辨率的粒子流數(shù)據(jù)，并且分別設置不同的粒子屬性(大小、顏色等)，控制顯示級別，在不同的屏幕視窗下展現(xiàn)不同級別的地形特征，進而達到地形粒子流的多尺度動態(tài)表達.

原則上應使非懸掛質量盡量小，這樣才能將非懸掛質量對振動特性(車身固有頻率)的影響降至最小。同時非懸掛質量的減小也降低了慣量所產生的沖擊負荷，大大改善懸架的響應特性，明顯提高乘坐舒適性。另外懸掛質量和彈簧的特征參數(shù)(彈簧剛度)，確定了車身固有頻率，懸掛質量較大或彈簧較軟，車身固有頻率就低一些，而彈簧的行程(振幅)就大。懸掛質量較小或彈簧較硬，車身固有頻率就高一些，而彈簧的行程(振幅)就小。試驗結果表明如果車身固有頻率超過1.5Hz時，會使乘坐舒適性惡化，超過5Hz，車輛行駛會強烈振動。

3實驗研究

實驗中所使用的數(shù)據(jù)是從地理空間數(shù)據(jù)云網站(http：//www.gscloud.cn/)下載的ASTER GDEMV2，武漢市共覆蓋五景.該DEM數(shù)據(jù)空間坐標系為GCS＿WGS＿1984，空間分辨率為30 m.

使用leaflet-velocity.js插件搭建地形動態(tài)可視化平臺，首先以0.005°、0.005°的分辨率對原始DEM數(shù)據(jù)進行重采樣，進行拼接、裁剪等預處理操作，生成水流場數(shù)據(jù)，載入平臺中實現(xiàn)水從高到低的動態(tài)流動效果.地理底圖在天地圖的地形渲染圖以及ArcGIS庫中的

World＿Street＿Map地圖中進行切換.將基于研究區(qū)提取得到的邊界線、河流、湖泊、匯流口等特征位置數(shù)據(jù)加載在底圖之上.設置合適的粒子屬性參數(shù)，水流效果如圖5所示.最終的動態(tài)可視化效果能夠呈現(xiàn)出研究區(qū)的地形特征及地勢高低.

圖5 粒子系統(tǒng)模擬水流界面Fig.5 The website interface of water flow simulated by particle system

在地形渲染地圖下，粒子具有向河流等低洼地區(qū)流動的整體趨勢，水流的流動方向及聚集情況能夠較好的體現(xiàn)出河段的位置及范圍.如圖6、圖7展現(xiàn)了武漢市黃陂區(qū)的灄水河、新洲區(qū)的舉水河周邊粒子流動的整體情況.此外，從地圖上還可以發(fā)現(xiàn)到研究區(qū)內的大量湖泊(如東湖、湯遜湖、魯湖等)的粒子流動特征：湖泊周邊的水流往湖區(qū)流動，湖內幾乎無水流.結合以上特征和在底圖上標注的重要特征位置，可以直觀地觀察得到研究區(qū)內的河流與湖泊所處的位置及周邊地形情況.

圖6 水流模擬下的灄水河Fig.6 Sheshui River based on water flow simulation

圖7 水流模擬下的舉水河Fig.7 Jushui River based on water flow simulation

基于原始DEM數(shù)據(jù)，設置8 000為匯水閾值，共提取武漢市河流647條，將其加載到底圖上，結合上述粒子流動的整體趨勢，通過視覺比較來評價該方法的合理性.結果表明，338條主要河流與周邊粒子運動效果一致，準確率達到52.24%.

由于使用匯水閾值方法提取的河流存在寬度一致的誤差且涵蓋了部分湖泊等問題，因此出現(xiàn)水流粒子在距河流一段距離內才產生聚集的情況.結合粒子表達河流及湖泊的特點，設置一定的河流誤差閾值(50 m)進行多次篩選，有609條河流均能夠被表達出來，準確率能夠達到94.13%，具體比較結果如表2所示.

表2 河流多次篩選結果Tab.2 Results of multiple River screening

一部分長度較短的河流不能被很好地展現(xiàn)出來，主要是由于短河段周邊的低洼點造成的.此外，將河流數(shù)據(jù)按長度進行升序排序，前75%均能在水流模擬系統(tǒng)中被表達出來.大量的粒子在匯流口處匯集，這也驗證了本文方法的準確性.

基于原始DEM數(shù)據(jù)提取山脊線、山谷線，加載在底圖之上.綜合地圖底圖，采用人工識別的方法，發(fā)現(xiàn)使用粒子系統(tǒng)支持下地形起伏動態(tài)表達方法能夠展現(xiàn)85%的山脊線及山谷線，較好地呈現(xiàn)出武漢市的地形起伏狀況.

武漢市黃陂區(qū)院基寺水庫附近有兩列東南走向的群山，如圖8黃線所示，水流從分水線向四周流動，最終匯聚到山谷線.研究區(qū)內的木蘭山、磨山、珞珈山等山區(qū)均呈現(xiàn)出這樣的特點.

圖8 山體水流動態(tài)模擬Fig.8 The dynamic simulation of flow on mountain

將原始DEM數(shù)據(jù)按照0.005°、0.005°、0.01°、0.01°及0.025°、0.025°的分辨率進行重采樣，生成三個不同尺度的水流場數(shù)據(jù).將粒子的顏色分別設置為紅色、綠色和藍色，其余屬性均設置為相同且適宜.在同一區(qū)域，單個水流向量場數(shù)據(jù)最終呈現(xiàn)的水流效果如圖9所示.

圖9 水流模擬示意圖Fig.9 The diagram of water flow simulation

在同一區(qū)域，進行不同空間分辨率粒子流的疊加，比較表明，較高的空間分辨率DEM數(shù)據(jù)生成的水流圖與實際地形更加接近和貼合，水流也更加接近真實.在相同尺度下， DEM數(shù)據(jù)的空間分辨率越高，即重采樣間隔越小，水流粒子在其生命周期內所流動的距離越短，粒子越稀疏，粒子向低地匯聚的次數(shù)越多，流動旋渦越小.而空間分辨率越低的DEM數(shù)據(jù)，水流粒子流動的距離越長，流動旋渦覆蓋面積越大，只能對DEM起伏呈現(xiàn)出整體特征，而忽略了地形的局部特征.

此外，本文選取武漢市北部的一個矩形范圍區(qū)域的DEM數(shù)據(jù)，生成水流場數(shù)據(jù)并實現(xiàn)水粒子的動態(tài)流動.其水流流動的效果仍然呈現(xiàn)出從高地向湖泊、河流、山谷等低洼地匯聚的規(guī)律.

圖10 矩形范圍水流模擬示意圖Fig.10 The diagram of water flow simulation in rectangular range

同時，根據(jù)該矩形區(qū)域DEM數(shù)據(jù)生成的水流場數(shù)據(jù)，嘗試調整表1中列舉的水流粒子各個屬性參數(shù).實驗發(fā)現(xiàn)，只有選擇合適的粒子屬性參數(shù)，如velocityScale、lineWidth、particleMultiplier、frameRate、colorScale等，才能更好地表現(xiàn)研究區(qū)地形的真實起伏，更具美感.通過反復調整屬性參數(shù)進行實驗，發(fā)現(xiàn)將屬性參數(shù)velocityScale設置為0.001、線寬設置為1、particleMultiplier設置為1/300、frameRate設置為15時，可以得到更好的地形流模式表達.圖11展示的是僅更改某些粒子屬性參數(shù)后的水流模擬效果圖.

圖11 不同粒子屬性參數(shù)的水流模擬圖Fig.11 The water flow simulation diagrams of different particle attribute parameters

4結語

本文使用粒子系統(tǒng)的水流模擬的方法，間接地表示研究區(qū)的高程信息.由于水流粒子具有隨機性和動態(tài)性，粒子系統(tǒng)平臺能夠呈現(xiàn)出對真實地形的動態(tài)模擬效果.該DEM表達方法更加美觀且符合人類視覺感受，且用戶的交互性強，用圖者可進行交互式查看和比較某一范圍內的地物相對高程，從局部和整體的角度把握區(qū)域的地形起伏狀況.并可以任意切換和疊加基于不同的空間分辨率數(shù)據(jù)模擬得到的粒子流數(shù)據(jù)圖層，從水流的流動規(guī)律中可以挖掘更深層次的信息.

與DEM靜態(tài)可視化方法相比，本文的方法更具有動態(tài)性，可以與其他方法相結合使用.此外，水流動態(tài)效果與粒子參數(shù)及數(shù)據(jù)空間分辨率有一定的關系，需要進行相關的實驗比較，才能得出適合研究區(qū)的數(shù)據(jù)空間分辨率及參數(shù)信息，更好地呈現(xiàn)出研究區(qū)的地表起伏狀況.

然而本文的DEM動態(tài)可視化方法僅停留在二維層面，因此對水流模型的精度和地形的三維動態(tài)可視化等方面還需要進一步展開研究，為DEM 可視化表達提供更好的技術支撐.