利用锥构建DEM生成算法的研究

摘　要　在综合基于离散点数模和三角网数模生成DEM两种方法 优点的基础上，提出了
一种基于离散无序点的快速生成DEM的方法。该方法有内插速度 快、精度高及占用内存
小的优点。
关键词　锥；内插；数字高程模型
分类号　P231.5
Study of Algorithm for Generating DEM Based Cone
JIANG Hongfei　ZHAN Zhenyan
(Civil Department,Changsha Railway University,154 Shao shan Road ,Changsha，
China,410075)
Abstract：DEM(digital elevation model) has a wide variety of appli cation in
GIS and CAD .It is the basic model for generating three_dimensional terrain
feature.Gene rally speaking, there are two methods for building DEM.One is
based upon dis ord ered points digital terrain model with features of fast s
peed a nd lo w precision .The other is built upon triangular digital terrain
model with the p eculiarities of slow speed and high precision.Compositing
the advantag es of the two methods,an algorithm for generating DEM with diso
rdered points is presented in the paper.When interpolating elevation, it can
create a trangle wh ich includes interpolating point and the elevation of t
he interpolating point ca n be got with the triangle. The method has the adv
antage of fast speed,high prec ision and requirement,for little memory.
Key words：cone;interpolate;digital elevation model
　　数字高程模型(digital elevation model,DEM)在地理信息系统及计算机辅助设计中
有着广泛的应用。它是生成三维立体地貌，计算地面坡度、土石方数量等数据的模型，
也是土 木工程中实现三维可视化设计的基础和前提条件。
　　文献［1］给出了一种利用Voronoi图 生成DEM的方法，该方法的核心是先利用离散
点数据生成Voronoi图，在生成DEM时找出插值 点所在的Voronoi多边形以及与该多边形
相邻的其他多边形，利用这些多边形以及所含的离 散点生成一个局部的小三角网，再利
用最小曲率原则求算插值点的高程。该方法的优点是精 度较高，其缺陷是必须先生成V
oronoi图，前处理工作量大，且保存Voronoi图需占用大量内 存。在生成DEM时还必须判
断插值点落在哪个Voronoi多边形中，这影响了该算法的计算速度 。文献［4］提出了一
种由离散无序数据生成DEM的方法，该方法采用数据分块技术、分块搜 索 技术可快速找
出与插值点相关的离散点，然后采用距离加权法进行插值，有较快的生成速 度。但这种
方法的精度不高，其原因是该方法容易引起数据的“跑偏”，即当插值点附近的 离散点
数据分布不均匀时，内插出的高程值容易失真。
　　基于此，本文提出了一种利用锥构建DEM 的算法，该方法利用格网管理离散点数据
，在内插插值点的高程时只要先确定插值点所在的 格，即可快速找出与插值点相邻的离
散点，然后利用这些离散点生成一个包含插值点的三角 形，用这个三角形来内插插值点
的高程值即可。
1　算法的基本原理
　　首先，本文所提出的内插算法是基于离散点数模的，即内插前应该先按常规方法建
立离散点数模；其次，为了提高内插精度，本算法在离散点数模中引入了地性线。
　　在内插插值点的高程时，离散点离内插点越近，对插值点高程的影响越大，因而最
理想的情 况是生成的三角形的3个顶点距内插点最近。基于这种想法，在形成内插三角
形时先选取距 插值点最近的两个点，再按某一种方法找出第三点，使这3点构成的三角
形包含插值点，然后利用该三角形作为内插三角形来内插出插值点的高程。由于参加构
建DEM的离 散点都是三维点，这里特别强调内插三角形的生成在水平面上进行，即插值
点落在某3 个离散点构成的三角形中是针对插值点和离散点在水平面上的投影而言的。

　　为了能方便阐述算法的基本原理，先给出锥的定义。
　　定义　设C是Rn中的非空集，又设x∈C，若p∈Rn,当x+p∈C时，必有x+tp∈C，t≥0
是任意实数，则把集合C称为以x为顶点 的锥。若锥C是凸集，则称为凸锥［5］。
　　显而易见，由向量V0,V1,…,Vm组成的所有非负组合的集合为：
C=xx=x0+riVi,ri≥0
它是一个以x0为顶点的凸锥。
　　特别地，对二维空间而言，以B点为顶点的
　　.两个向量BA、BC所张成的凸锥是一个以B点为顶 点，其值小于π的∠ABC。
　　如图1所示，O是插值点，A、B是距插值点O最近的两点，现在要找出一点P，△ PAB
包含O点，并使P点尽量接近O点及△PAB的形状最好。为了找出点P，可找出落在以O为顶
点，由向量OC(即-OB)及OD(即-OA)张成的凸锥R中的点， 在这些点中找出任一点P，则△
ABP必定包含点O。
图1　内插三角形的生成原理
Fig.1　The Pri nciple of Creating an Interpolating Triangle
　　证明　由图1可得：
OP=r1AO+r2BO　　　(1)
PB=-OP-BO　　　(2)
PA=PO+OA=-OP -AO　　　(3)
由式(1)、(2)、(3)得：
PO=(r1PA+r2PB)/(1+r1+r2)　　　(4)
PA=PO+OA=(1+r1)OA+ r2OB　　　(5)
AB=OB-OA　　　(6)
由式(5)及式(6)得：
AO=(AP+r2AB)/(1+r1+r2)　　　(7)
于是有下列结论：
　　1)当P点落在凸锥R中时，有r1≥0且r2≥0，故r1/(1+r1+r2)≥ 0，r2/(1+r1+r2)≥
0。由式(4)及式(7)知，点O既在以P为顶点，由向量PA及PB张成的 凸锥中，又在以A为顶
点，由向量AP及AB张成的凸锥中，故△PAB必包含 点O。
　　2)当P点不在凸锥R中时，则r1及r2至少有一个小于0，故r1/(1+r1+r2)、r2/ (1+r1
+r2)、1/(1+r1+r2)必不会同时大于0。由式(4)及式(7)知，点O不能同时既在 以 P为顶
点，由向量PA及PB张成的凸锥中，又在以A为顶点，由向量AP及AB张成的凸锥中,故此 时
△PAB必不包含点O。
　　由此有下列定理：
　　定理　对于任意给定的不在同一直线上的3点A、B、O，点O落在△PAB中的充分条件
是：P点在以O为顶点，由向量AO及BO张成的凸锥中。
　　实际上，生成三角形时，首先找出离插值点O最近的两个点A、B；然后以O为顶点，
由向量AO及BO张成一个凸锥R；其次在O点周围的点中找出点P，使P点在凸锥R中，则△
PAB即为所求的三角形。当符合条件的P点有多个时，可选取使△PAB的最小角最大的 那
个点作为内插三角形的第三个顶点，避免狭长三角形的产生。当然也可以采用别的原则
来选择内插三角形的第三个顶点。
　　显而易见，凸锥R中不一定有点存在，这样就找不到符合条件的点P；但以其他两点
为基边生 成包含插值点的三角形有可能成功，因而必须考虑这种情况。
　　如果以AB为基边时找不到P点，此时显然以O为顶点，由向量AO及BO张成的凸锥R内无
点，如图 2所示，用O点对基边所张的角来表示已搜索的范围，目前的搜索角θ=∠AOB。
由于不能形 成符合要求的三角形，而A点又是距待定点O最近的点，显然此时应在以O为
顶点、由向量BO及OA张开的凸锥中再取距 O点最近的点C，然后以AC为基边再找P点。如
果点C或点P不存 在,那么当点C不存在时，置θ=π；C点存在而P点不存在时，置θ=θ+
∠AOC。此时如果 θ< π，再在以O为顶点，由向量OB及CO张成的凸锥中找一点距O最近
的点D，以BD为基边 找点P 。如果点D或点P不存在，那么当点D不存在时，置θ=π；当
P点不存在时，置θ=θ+ ∠BOD 。此时如果θ<π，再重复上述步骤，直至找到P点或找
不到作基边的点(θ≥ π)为止。
图2　搜寻内插三角形的方法
Fig.2　The Met hod for Searching Interpolating Triangle
2　地性线的引入
　　按照前述方法，每次内插插值点的高程时都要先生成一个内插三角形，三角形是否
能够密贴地面，真实反映出地表形状将直接影响内插精度。由于在生成三角形时采用了
最大最小角法则，因而能有效地杜绝狭长三角形的产生，对一般地形而言，三角网能密
贴地面，真实反映出地形情况。但在各种断裂线及山谷线、山脊线附近，三角形将会跨
越这些地性线，此时，三角形不能密贴地面，而是悬浮在空中或穿入地下，从而使内插
精度受到影响。解决这种情况最有效的办法是在采点时引入地性线，使三角形只能在某
条地性线的一边生成，不准 其跨越任何一条地性线，这样就能保证其内插精度。
　　在引入地性线时，可用若干个有序的三维点形成的折线来表示山谷线和山脊线及其
他地性线。为了在插值时能迅速找到与插值点相关的地性线，仍然采用格网对地性线进
行管理。具体方法是记录每格中地性线的条数及坐标，为此必须将每条地性线离散到格
网中，即如果地性线不只在某一格中，求出地性线与格边界的交点并将交点及其落在格
中的端点顺序存放在该格的有关数据变量中。
　　通过上述处理后，在高程内插时，不仅可快速找出与插值点相关的离散点，还可以
找出与插值点相关的地性线。如图3所示，在内插O点的高程时，此时有两条与O点相关的
地性线。为防止生成的三角形跨越地性线，在构造内插三角形之前，按下述方法对参加
构造内插三角形的离散点进行筛选。如果要防止生成的内插三角形跨越AB，则以向量OA
、OB张成一个凸锥；如果某离散点P在该凸锥内且O、P分别位于线段AB的两侧，则该点不
参加内插三角形的构造。用每条与插值点相关的地性线对离散点进行筛选后，即能防止
生成的内插三角形跨越地性线。
图3　在离散点数模中引入地性线
Fig.3　Adding Geographical Feature Lines in
Disordered Digital Terrain Model
3　结　语
　　为了检测本文所提高程内插方法（移动三角形内插方法）的精度、速度，先利用离
散点生成离散点数模和三角网数模，然后在此基础上分别采用三角网数模高程内插方法
、离散点数模按距离加权高程内插方法及移动三角形高程内插方法对同一区域进行DEM的
生成。检测时采用的计算机主频为133MHz，硬盘6G，内存32M。参加构建数模的离散点约
为25万个，总共布置内插点5万个，结果利用三角网数模内插用了6min，离散点数模按距
离加权高程内插方法用了65s，移动三角形高程内插方法用了90s，表明移动三角形高程
内插方法具有很快的内插速度。
　　在生成DEM后，将按3种方法内插出的部分高程值（1 000个点）进行对比，发现按距
离加权 高程内插方法内插出的高程值误差在1m以上的点约占总内插点数的10%，最大的
相差5.1m；移动三角形高程内插方法的相应值为1.8%和1.72m；三角网数模内插方法的相
应值为2.3 %和1 .80m。内插值与实际地形高程值相差较大的主要原因是离散点采集密度
太稀，从而使离散点在个别地方无法真实描述出地表形状。此外，上述数据表明，移动
三角形高程内插方法有很高的内插精度，其精度甚至比三角网数模内插出的还要高。其
原因是移动三角形内插方法在生成内插三角形时，是专门针对内插点来选择三角形的顶
点的，其3个顶点是对内插点影响最大的；而三角网数模在生成内插三角形时，是从整个
三角网的角度来考虑三角网的形状的，但对于某些特 定的内插点而言，某些三角形并不
是最佳的三角形。因而在某些情况下，移动三角形内插方 法的精度反而比三角网数模的
内插精度要高。
作者简介：蒋红斐，男，32岁，讲师，现主要从事铁路线路的计算机辅助设 计研究。代
表成果：《基于复杂地形上三角网数字地面模型的建立》、《一般窗口的一种实 现方法
》等。
　　　　　E_mail:yizhang@csru.edu.cn
作者单位：蒋红斐（长沙铁道学院土木系,长沙市韶山路154号，410075）
　　　≌舱裱祝ǔど程?姥г 土木系,长沙市韶山路154号，410075）
参考文献：
[1]　Janos Kalmar.DEM_based Surface and Volume Approximation —— G eographi
cal Applications.Computer & Geosciences,1995,21(2):147～163
[2]　Tsai，Delaunay V J D.Triangulations in TIN Creation:an Overview an d a
Linear_time Algorithm.Int. J. of GIS,1993,7(6):501～524
[3]　Nagy，Wagle.Geographic Data Processing.Computing Surveys,1979，11(2):1
39～181
[4]　王永明，林行刚.一种快速DEM生成算法.计算机应用与软件，1998，15(4):28～3
3
[5]　薛嘉庆.最优化原理与方法.北京:冶金工业出版社,1983.259～266