刘鹏程， 黄 欣*， 郭明武

(1.华中师范大学地理过程分析与模拟湖北省重点实验室， 武汉 430079；2.华中师范大学城市与环境科学学院， 武汉 430079； 3.武汉市测绘研究院， 武汉 430022)

“粒子系统”是在1983年由William T Reeves初次提出的[1]，用以描述自然界中一些运动不存在规律性的事物，目前已成功模拟了爆炸、焰火等动态运动效果[2].粒子系统的基本原理是将模拟的事物定义为由无数个随机产生且在其有限的生命周期内做不规则运动的粒子构成[3].每个粒子需经历“产生”“运动”以及“消失”三个阶段，粒子运动到下一位置生命周期减1，每个时刻均有新粒子的产生和旧粒子的消失.通过在某一时刻 “运动”粒子的速度及方向等信息可以展示事物的整体特征和动态变化.粒子系统较为充分地呈现事物运动的随机性与动态性，逼真地模拟事物的变化[2].

在模拟不规则运动物体的研究领域，粒子系统是学术界公认的最成功的图像模拟生成算法[4]，具备重大的研究意义与广阔的发展前景.国内外的专家学者使用粒子系统的方法成功模拟出了很多的自然事物，并不断地对模拟方法进行改进，以追求更高的真实和实时性.

1985年，Reeves和Blau[5]提出“volume filling”单元用以模拟随着时间的变化，形状发生改变但整体基本不变的事物，如随风飘动的花草树叶.1990年，Sims[6]首次实现了使用计算机成功绘制粒子并呈现出粒子运动的动画效果，提出了粒子系统的并行绘制算法.该算法实现了电影“Particle Dream”中火焰、爆炸、瀑布等影视效果.Richard和David[7]使用粒子系统对弹性物体表面进行建模.此外，Unbescheiden和Trembilski[8]从云的物理原理出发，将粒子系统与纹理映射技术相结合，建立了云的模型.万华根等[9]基于粒子系统和流体动力学模拟出了喷泉中水的实时流动效果.谢剑斌等[10]通过粒子系统模拟了雨点和雪花降落过程.Kolb等[11]采用 GPU 进行粒子系统的模拟，使得大规模进行粒子的模拟变成了可能.付帅等[12]利用粒子系统对全球表层流场进行仿真表达.Reynolds[13]使用粒子系统模拟了动物群体的运动.此外，Goss[14]通过粒子系统实时地模拟了船舶行驶时的船迹.

数字高程模型(digital elevation model ，DEM)是格网数据点的二维平面坐标(X，Y)和高程值(Z)的数据集[15]，用来表达地表的起伏和特征[16].目前，DEM的可视化表达方法主要在一维剖面、二维平面以及三维立体等空间维度展开，使用如等高线法[17]、阴影法[18]、明暗等高线法[19]、地形剖面线法[20]等静态表达方法表达地形地貌.此类方法一般使用一幅图形或图像来对地形起伏状况进行表达，其信息的获取与规律的挖掘主要取决于用图者分析的能力.

本文尝试利用粒子系统模拟水体在地表上的运动来形象描述地表的高低起伏状况，其中粒子的流向和流速可以用以表示地表单元的坡向和坡度.这为地形的表达提供了一种新的动态可视化方法，能够很好地将局部和全局的地形地貌特征展现在观察者面前，并实现交互式、动态性的地形浏览.

流模式是指利用粒子系统的方法构建水流模型，用以模拟水流在地表上的流动进而实现对地形起伏情况的动态可视化表达.其中，粒子系统将水流视为是无数随机分布的粒子集合，每个粒子都有其固有的属性，如初始位置、流动的方向和速度、颜色、生命周期等.每个粒子都能在水流场中根据其初始状态不断地自主进行不规则运动，模拟出接近真实流动形态的动态水流，实现对DEM的动态可视化表达.基于粒子系统的水流动态模拟主要流程包括以下几步.

1) 水流场数据转化为二维数组数据.由于粒子的生成位置是随机且大量的，所以为了便于快速地计算水流场中任何一个粒子在当前格网的水流方向以及速度值，需要将水流场数据转化为二维数组数据.二维数组中的每一项均存储当前格网的横向水流速度(u)和纵向水流速度(v).

2) 生成若干个粒子，并计算粒子的运动位置.二维数组数据相当于基础数据，数据量有限，但是在大多数情况下，粒子的数量要远大于该数据量.所以想要更为详细地展现研究区范围的水流情况，即地势的高低起伏，就需要在网格上创建较多的粒子来进行动态可视化显示.本文主要采用两种方法来确定粒子的位置，一种是利用Math.random()函数在水流场范围内生成大量粒子，使粒子能够在整个研究区域内随机分布.另一种方法是通过查看粒子的流向走势可以获取地形的特征，如山脊线、山谷线等，分析总结粒子在各类地形中的流动规律对地形特征进行提取，进而实现在重要特征地形上粒子的自动放置.这样既能反映研究区的总体地形特征，又能逼近局部关键地形.

3) 计算粒子流动位置并进行绘制.多次调用animate方法，用于计算每个粒子的当前位置、水流的速度及方向、下一次运动要抵达的位置.每次使用animate方法，粒子的生命周期参数就会减少1，直到其值变化为0.当粒子生命周期为0或者粒子移动到水流速度为0处，会使得粒子消失进而导致水流场中的水流粒子越来越少.为解决此类问题，需对粒子生命周期进行重新赋值，每个时刻均有粒子的产生和消失，确保水流场内一直有水流在流动.

在Canvas上进行当前运动粒子所构成水流的绘制，实现水流的动态可视化.选择合适的地理底图，在相同的坐标系统下实现粒子与底图的位置匹配.每次进行地图的缩放操作时，计算出粒子的数量及运动状态，根据当前地图所在范围进行水流的重新绘制.同时，还能实现在不同大小的窗口中自适应显示粒子的流动.

4) 将重要位置标注在底图上.将基于DEM数据提取得到研究区主要河流、湖泊、汇流口等特征位置的矢量数据转化为GeoJSON格式，加载到地图底图上进行表示，可以增强地形起伏的可视化效果，提高美观性.

1水流场数据的生成

粒子的流动是非周期性的，水流的强弱和方向可以使用流速和流向来表示，故可将水流数据定义为向量场数据，在向量场中每个点对应的数据均为一个向量，某点的向量方向代表水流在该位置的流向，向量大小代表流速.在高低起伏的地面单元上粒子倾向于运动，其运动的驱动力源自地表起伏所形成的各个坡度，且驱动力的方向与坡向一致.将驱动力分解为u和v两个正交方向，构建二维水流场数据.其具体步骤如下.

1.1构建水流模型

水流场数据的生成需要考虑到地面坡度对水流大小及方向的影响.使用D8算法[21](deterministic eight-node)来构建水流模型.D8算法假设某中心栅格单元在3×3窗口中有8种水流的可能方向，即其所在位置的八邻域方向.使用公式(1)分别计算中心栅格单元与八邻域的格网单元的高程差与距离的比值，即坡度.最大坡度值所对应的方向为中心栅格单元的水流方向.

(1)

在计算中心网格单元与其周围八个相邻栅格单元的高程差时，存在两种情况：若高程差均小于或等于0，则该中心单元无水流流出，粒子速度为0；否则，中心单元的水流流向坡度最大值所对应的网格单元，粒子的初速度值为坡度的倍数.

图1展现了对DEM数据进行逐网格遍历，计算粒子流动方向的过程：在3×3窗口中计算得到中心网格单元的最大坡度值，并确定水流流动方向，为了方便后序粒子分量的内插求解，将栅格数据转化为具有方向属性的点数据.

图1 水流方向的计算过程Fig.1 The process of calculating flow directions

1.2计算水流速度

各个网格单元的水流速度与最大坡度成正比，设置高差缩放因子对单元的水流速度进行调整从而优化粒子流的效果.如图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为缩放因子，用以调整研究区高程差；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生成水流场数据

将基于DEM数据进行逐网格计算得到的水流粒子速度u、v分量，与研究区地理坐标范围、网格的行列数和空间分辨率等信息相结合，生成JSON格式的水流场数据.图3展现了JSON文件的基本格式，主要包括两个数据块，分别用于存储u、v分量的相关数据，每个数据块均包含{“header”：{ }，“data”：[ ]}的结构.

“header”结构存储一些关于水流场的属性参数信息，一般JSON数据中两个数据块的“header”内容基本一致.“data”结构存储水流速度的数据，每个值对应研究范围内的一个网格点.第一个“data”存储的是各网格点的横向水流速度，即u分量，第二个“data”存储的是各网格点的纵向水流速度，即v分量.两个“data”中的值对应进行向量运算，即可计算得到该点的水流向量值.

2地形起伏动态模拟

2.1双线性插值法计算水流分量

在粒子系统中，水流是由无数个随机分布的粒子所构成的集合，每个粒子都能在水流场中根据其初始状态不断地自主进行运动.粒子系统能够模拟出接近真实流动形态的动态水流，其驱动力是地面的势能场，人们可以通过粒子流的运动来感知地形的起伏，进而实现对DEM的动态可视化表达.

由于水流场数据文件记录的是中心网格点的水流速度，而随机生成的粒子，其空间位置不可能完全位于中心网格点上.故粒子所在位置的速度需要依据粒子相邻网格点的速度值进行插值运算得到.本文使用双线性插值算法(Bilinear Interpolation)，来获取该粒子u、v分量的近似值，并计算出该水流粒子下一次运动所抵达的位置.

图3 JSON数据格式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分量值.

根据粒子当前所在位置和计算得到的u、v分量，即可计算得到粒子下一次运动到的位置，依次计算下去.每计算1次，粒子的生命周期参数就会减少1，直到其值变化为0.当粒子生命周期为0或者粒子移动到水流速度为0处，粒子将消失.

2.2粒子位置布设及属性说明

粒子初始位置的布设使用两种方法，其一使用随机函数在数据范围内随机生成，粒子的数量根据水流场数据的分布规律而变化，也可以通过调节参数来进行控制.其二通过查看粒子的流向走势可以获取地形的特征，实现在一些重要地形特征点或者特征线(如山脊线和山谷线等)粒子的自动放置，专门针对这些特征点进行水流的动态模拟，充分展现局部地貌骨架线.

随机生成的粒子属性需要进行初始化，粒子属性包括粒子的初始位置、流动的方向和速度、颜色、样式及生命周期等，表1中列出了一些粒子属性参数和相关描述.设置不同的粒子属性参数值，会产生不同的动力学效应，水流粒子的运动也会呈现出不同的动态效果.在实际的应用中，可以通过调整粒子流动速度的比例因子即公式(4)、(5)中的k值，对区域整体的粒子流动速度进行夸张，使粒子流动的整体效果能够更好地反映地表起伏和特征.

表1 水流粒子的属性参数说明Tab.1 The parameter description of water flow particle attributes

2.3粒子流的多尺度表达

通过设置水流场的网格大小、粒子的密度以及在地形特征线上基于不同尺度放置粒子，可以实现不同空间分辨率地形粒子流的动态可视化表达.当对地图进行缩放操作时，平台将根据当前屏幕窗口范围进行计算，重新绘制粒子并自动切换与当前比例相对应的底图.

为了获得更好的地形起伏效果，有时会在底图上覆盖不同空间分辨率的粒子流数据，并且分别设置不同的粒子属性(大小、颜色等)，控制显示级别，在不同的屏幕视窗下展现不同级别的地形特征，进而达到地形粒子流的多尺度动态表达.

原则上应使非悬挂质量尽量小，这样才能将非悬挂质量对振动特性(车身固有频率)的影响降至最小。同时非悬挂质量的减小也降低了惯量所产生的冲击负荷，大大改善悬架的响应特性，明显提高乘坐舒适性。另外悬挂质量和弹簧的特征参数(弹簧刚度)，确定了车身固有频率，悬挂质量较大或弹簧较软，车身固有频率就低一些，而弹簧的行程(振幅)就大。悬挂质量较小或弹簧较硬，车身固有频率就高一些，而弹簧的行程(振幅)就小。试验结果表明如果车身固有频率超过1.5Hz时，会使乘坐舒适性恶化，超过5Hz，车辆行驶会强烈振动。

3实验研究

实验中所使用的数据是从地理空间数据云网站(http：//www.gscloud.cn/)下载的ASTER GDEMV2，武汉市共覆盖五景.该DEM数据空间坐标系为GCS＿WGS＿1984，空间分辨率为30 m.

使用leaflet-velocity.js插件搭建地形动态可视化平台，首先以0.005°、0.005°的分辨率对原始DEM数据进行重采样，进行拼接、裁剪等预处理操作，生成水流场数据，载入平台中实现水从高到低的动态流动效果.地理底图在天地图的地形渲染图以及ArcGIS库中的

World＿Street＿Map地图中进行切换.将基于研究区提取得到的边界线、河流、湖泊、汇流口等特征位置数据加载在底图之上.设置合适的粒子属性参数，水流效果如图5所示.最终的动态可视化效果能够呈现出研究区的地形特征及地势高低.

图5 粒子系统模拟水流界面Fig.5 The website interface of water flow simulated by particle system

在地形渲染地图下，粒子具有向河流等低洼地区流动的整体趋势，水流的流动方向及聚集情况能够较好的体现出河段的位置及范围.如图6、图7展现了武汉市黄陂区的滠水河、新洲区的举水河周边粒子流动的整体情况.此外，从地图上还可以发现到研究区内的大量湖泊(如东湖、汤逊湖、鲁湖等)的粒子流动特征：湖泊周边的水流往湖区流动，湖内几乎无水流.结合以上特征和在底图上标注的重要特征位置，可以直观地观察得到研究区内的河流与湖泊所处的位置及周边地形情况.

图6 水流模拟下的滠水河Fig.6 Sheshui River based on water flow simulation

图7 水流模拟下的举水河Fig.7 Jushui River based on water flow simulation

基于原始DEM数据，设置8 000为汇水阈值，共提取武汉市河流647条，将其加载到底图上，结合上述粒子流动的整体趋势，通过视觉比较来评价该方法的合理性.结果表明，338条主要河流与周边粒子运动效果一致，准确率达到52.24%.

由于使用汇水阈值方法提取的河流存在宽度一致的误差且涵盖了部分湖泊等问题，因此出现水流粒子在距河流一段距离内才产生聚集的情况.结合粒子表达河流及湖泊的特点，设置一定的河流误差阈值(50 m)进行多次筛选，有609条河流均能够被表达出来，准确率能够达到94.13%，具体比较结果如表2所示.

表2 河流多次筛选结果Tab.2 Results of multiple River screening

一部分长度较短的河流不能被很好地展现出来，主要是由于短河段周边的低洼点造成的.此外，将河流数据按长度进行升序排序，前75%均能在水流模拟系统中被表达出来.大量的粒子在汇流口处汇集，这也验证了本文方法的准确性.

基于原始DEM数据提取山脊线、山谷线，加载在底图之上.综合地图底图，采用人工识别的方法，发现使用粒子系统支持下地形起伏动态表达方法能够展现85%的山脊线及山谷线，较好地呈现出武汉市的地形起伏状况.

武汉市黄陂区院基寺水库附近有两列东南走向的群山，如图8黄线所示，水流从分水线向四周流动，最终汇聚到山谷线.研究区内的木兰山、磨山、珞珈山等山区均呈现出这样的特点.

图8 山体水流动态模拟Fig.8 The dynamic simulation of flow on mountain

将原始DEM数据按照0.005°、0.005°、0.01°、0.01°及0.025°、0.025°的分辨率进行重采样，生成三个不同尺度的水流场数据.将粒子的颜色分别设置为红色、绿色和蓝色，其余属性均设置为相同且适宜.在同一区域，单个水流向量场数据最终呈现的水流效果如图9所示.

图9 水流模拟示意图Fig.9 The diagram of water flow simulation

在同一区域，进行不同空间分辨率粒子流的叠加，比较表明，较高的空间分辨率DEM数据生成的水流图与实际地形更加接近和贴合，水流也更加接近真实.在相同尺度下， DEM数据的空间分辨率越高，即重采样间隔越小，水流粒子在其生命周期内所流动的距离越短，粒子越稀疏，粒子向低地汇聚的次数越多，流动旋涡越小.而空间分辨率越低的DEM数据，水流粒子流动的距离越长，流动旋涡覆盖面积越大，只能对DEM起伏呈现出整体特征，而忽略了地形的局部特征.

此外，本文选取武汉市北部的一个矩形范围区域的DEM数据，生成水流场数据并实现水粒子的动态流动.其水流流动的效果仍然呈现出从高地向湖泊、河流、山谷等低洼地汇聚的规律.

图10 矩形范围水流模拟示意图Fig.10 The diagram of water flow simulation in rectangular range

同时，根据该矩形区域DEM数据生成的水流场数据，尝试调整表1中列举的水流粒子各个属性参数.实验发现，只有选择合适的粒子属性参数，如velocityScale、lineWidth、particleMultiplier、frameRate、colorScale等，才能更好地表现研究区地形的真实起伏，更具美感.通过反复调整属性参数进行实验，发现将属性参数velocityScale设置为0.001、线宽设置为1、particleMultiplier设置为1/300、frameRate设置为15时，可以得到更好的地形流模式表达.图11展示的是仅更改某些粒子属性参数后的水流模拟效果图.

图11 不同粒子属性参数的水流模拟图Fig.11 The water flow simulation diagrams of different particle attribute parameters

4结语

本文使用粒子系统的水流模拟的方法，间接地表示研究区的高程信息.由于水流粒子具有随机性和动态性，粒子系统平台能够呈现出对真实地形的动态模拟效果.该DEM表达方法更加美观且符合人类视觉感受，且用户的交互性强，用图者可进行交互式查看和比较某一范围内的地物相对高程，从局部和整体的角度把握区域的地形起伏状况.并可以任意切换和叠加基于不同的空间分辨率数据模拟得到的粒子流数据图层，从水流的流动规律中可以挖掘更深层次的信息.

与DEM静态可视化方法相比，本文的方法更具有动态性，可以与其他方法相结合使用.此外，水流动态效果与粒子参数及数据空间分辨率有一定的关系，需要进行相关的实验比较，才能得出适合研究区的数据空间分辨率及参数信息，更好地呈现出研究区的地表起伏状况.

然而本文的DEM动态可视化方法仅停留在二维层面，因此对水流模型的精度和地形的三维动态可视化等方面还需要进一步展开研究，为DEM 可视化表达提供更好的技术支撑.