2021地理设计组二等奖：DEM的构建及地形特征线的提取研究

作品简介

摘要

格网DEM是数字地面模型的重要数据表现形式，基于格网DEM提取的山脊线、山谷线是地形骨架结构的表达与数字地形分析的重要方面。本文通过比较距离倒数乘方法、自然邻域法、克里金法、样条函数法、趋势面分析法5种空间内插方法在本溪关门山地区的适用性，选出最优法用于格网DEM的构建。基于格网DEM，对比运用“水文分析法”和“曲率判别法”两种方法提取山脊、山谷线的结果，并研究两种提取算法的融合效果。以GeoScene Pro建模和脚本的形式分别实现了DEM构建、等高线曲率计算，山脊线、山谷线提取；制作了本溪关门山地区数字高程模型图、等值线专题图、山脊线山谷线专题图等。提出了坡度、平面曲率与等高线曲率的转换方法；对水文分析法和曲率判别法做出了改进，并提出了基于水文分析法和曲率判别法的融合算法，进一步提高了山脊线山谷线提取的准确率。融合算法兼具二者之长，提取结果较为全面，整体性较好。

0 引言

数字高程模型（Digital Elevation Model，DEM）的数据结构简单、容易管理、便于存贮，许多国家的DEM数据都是以规则格网的数据矩阵形式提供，因此，基于规则格网DEM提取等高线、山脊线、山谷线一直是研究人员的研究方向[1]。山脊线、山谷线是地形骨架结构的重要组成部分[2]，常应用于地貌表达以及三维地形分析。

建立格网DEM的空间内插方法可以分为分块内插法、整体内插法、逐点内插法3类。其建立一般包括数据取样、数据内插和数据精度分析等步骤[3]。整体内插法的原理简单，适用于总体预测，但不适用于复杂地形，对于复杂地形，常使用分块内插法[4,5]。逐点内插法是应用最为广泛的方法，包括距离倒数乘方法、克里金法。当样本点分布均匀且数量较多时均取得较好的结果。当实际地形的坡度较大时，采用距离倒数乘方法较为适宜[6]。

提取山脊、山谷线常用的算法有3类，一是几何分析法，它利用地表几何形态的实际规律，通过一定的算法（如曲率判别法、断面极值法和骨架化法）提取地形特征点，其中断面极值法与曲率判别法均是利用极值提取地形特征点。断面极值法一般采用两个方向的断面，并不能提取出所有地形特征点[7-9]。骨架化法假设山脊、山谷线两侧地形对称，提取中心轴线，与实际不符，因此，很大程度的近似性，所得结果差强人意[10,11]。二是图像处理法，它将DEM数据作为灰度图像处理，通过数学形态学的“序贯减薄”运算得到骨架线[12]，该方法可利用移动窗口法提取特征点[13]，但是将特征点连成线非常困难。三是水文分析法，它利用地表水流的客观规律对山脊、山谷线进行提取，常用D8流向法计算流向[14]，在带有汇聚型河流的多山地形能生成较好结果。但在平坦地区不易确定流向，不适用于含泛滥平原和湿地的多变地形[15]。

本文对格网DEM构建常用的五种算法（距离倒数乘方法、自然邻域法、趋势面分析法、克里金法、样条函数法）进行了实现，并分析它们在本溪关门山地区的适用性与准确性，选出最适合的插值方法“距离倒数乘方法”和“自然邻域法”用于格网DEM的构建。基于格网DEM，利用“水文分析法”和“曲率判别法”两种方法提取山脊线、山谷线，将提取结果进行对比分析，并研究两种山脊、山谷线提取算法的融合，实现了融合算法提取山脊、山谷线。

1 格网DEM构建及山脊、山谷线提取的算法原理与技术路线

1.1 技术路线

本文的主要任务是针对本溪关门山地区，采用不同的空间插值方法、山脊、山谷线提取方法，通过GeoScene建模和二次开发的形式分别完成格网DEM的构建、等高线的绘制、等高线曲率的计算、山脊线和山谷线的自动化提取。同时分析不同插值算法的适宜性，以及不同山脊、山谷线提取算法的优缺点。技术路线图（图1）如下：

第一步：格网DEM的构建，利用五种空间插值算法进行栅格插值，并通过交叉验证选出最合适的方法用于本溪关门山地区格网DEM的构建。

第二步：等高线的自动绘制，基于格网DEM数据，生成等高线。

第三步，等高线曲率的计算，基于等高线数据，计算等高线曲率。

第四步：曲率判别法的实现，基于等高线曲率，完成山谷线及山脊线的提取。

第五步：水文分析法的实现，基于水文分析法，完成山谷线及山脊线的提取。

第六步：水文分析法与曲率判别法提取结果的对比分析以及两种算法的融合。

图1 总体技术路线图

2 空间插值算法实现与分析

本文实现了五种空间插值算法构建格网DEM，包括距离倒数乘方法、克里金法、样条函数法、趋势面分析、自然邻域法。在此基础上，本文对插值结果进行对比分析。

在进行插值之前，首先将样本点分为两部分，其一为训练点，其二为测试点，样本点共有37660个，训练点为1000个，测试点为36660个，其中训练点用于检测插值质量，测试点用于格网DEM的构建。

在插值成功之后，将5种插值结果进行3次对比，即将样本点随机分配三次，进行独立试验，互不干扰。每次对比的测试点与训练点的分配比例相同，但分布各不相同，均是随机的。

根据测试点的5种插值算法插值结果，将训练点所在位置的值分别提取至训练点属性表，共分别提取5次。之后便可计算每种插值算法在每个训练点处的偏差。基于此属性表，进行统计分析，分别计算均值、最大值、最小值、标准差、方差。统计表（表1）如下：

表1 三次插值结果对比

图2 距离倒数乘方法（a）、克里金法（b）、样条函数法（c）、趋势面法（d）和自然邻域法（e）的插值结果

综合表1以及图2可以看出，本实验区（本溪关门山地区）应用的5种插值算法中，距离倒数乘方法、自然邻域法的均方根误差为4.5330和3.3393，均值为-0.1622和-0.0738，插值结果最好，最能反映真实地形。克里金法和样条函数法稍差，均方根误差为5.6032和5.7489，均值分别为-0.2934、-0.2381，插值结果不如距离倒数乘方法和样条函数法，但差距不大，也能反映真实地形，最后是趋势面分析法，该方法插值结果极差，标准差123.7280，偏差较大，只能反映整体趋势，不能反映真实地形，不适用于本试验区。

3 山脊、山谷线提取模型的建立

3.1 曲率计算模型

等高线曲率与GeoScene中的曲率有根本上的不同，它是指地表面上通过某一点的等高面与地表的交线在该点的曲率，反映的是地形表面在水平方向的凹凸性[19]。

图3 等高线曲率(C1)、平面曲率(C2) 和剖面曲率(C3)

如图3所示，地形表面S的曲面方程为

，点P为曲面S上任意一点，平面α为过点P的水平面，α与S的交线为过点P的等高线；为等高线在点P处的切向量，为P点的坡向向量。平面β、γ 分别过 、垂直于水平面XOY并与曲面S相交于曲线。曲面S：的各阶偏导数用以下符号表示：

…（3）

在GeoScene 中提供了平面曲率的计算工具，该曲率是垂直于坡度方向的曲率。从图3中的数学关系来看，平面曲率就是曲线

在P点处的二阶导数，它的计算公式为： 

…（4）

在GeoScene 中，用最大平均值法计算坡度值。GeoScene中的坡度 可利用式（3）中的变量进行表达，其计算公式如下：

…（5）

地形曲面S：

的任意一条等高线方程可表示为（c 为表示任意高程值的常数），由曲率计算公式可得等高线曲率计算公式为：

…（6）

通过综合公式（4）、（5）、（6），可以发现平面曲率、等值线曲率、坡度三者之间具有一定的关系，通过平面曲率与坡度便可计算等值线曲率栅格。

本文推导出的公式如下：

…（7）

曲率判别法的实现思路是：首先获取等值线曲率栅格，其次基于曲率栅格获取地形的特征点，通过一定的分类阈值提取山脊点与山谷点。模型图如附录A所示，基于该模型，山谷线提取结果如图6（a1）所示，山脊线的提取结果如图6（a2）所示。

3.2 水文分析模型

水文分析法的基本原理是：对于山脊线而言，它也是分水线，分水线所在的栅格区域在一定程度上不存在水流，即流量累积栅格的对应值为零值。通过地表径流模拟计算可以得到流量累积栅格，并提取零累积值区域，即分水线所在区域，也就得到了山脊线[20]；其次，若以某一相对高程面为对称基准面，求原始地形的对称地形，则山谷线与山脊线也位置互换。因此，可以通过本文建立的反地形提取模型，求出对称地形（反地形），那么，正地形山谷线就可以在反地形中利用提取山脊线的方法进行提取[21]。

1）已填洼DEM

水文分析法中已填洼DEM是指不存在小洼地，大多数小洼地是DEM生成过程中带来的数据错误所致，如果未进行填充，则生成的水系网络可能会呈现不连续性，常用的除去洼地的方法是把洼地像元值加高至周围最低像元值[22]。填洼操作不仅能去除小洼地，同时也可以去除错误的峰，峰是一种伪像元，其高程高于所预期的高程，是反地形中的小洼地。

2）反地形DEM

传统的水文分析法提取山谷线过程中，反地形的计算是基于原始DEM数据的，本文利用已填洼DEM数据提取反地形，从而消除反地形中的峰和小洼地，进而消除山谷线提取结果的局部错误值。反地形的计算公式为：

…（8）

式中Abs是绝对值运算，Dem代表已填洼DEM数据，

分别代表原始DEM的最大高程值和最小高程值。

模型中山脊、山谷线提取阈值的确定和反地形的计算极为重要。对汇流累积量零值数据进行重分类时，属性值越接近1越可能是山脊线或山谷线的位置，分界阈值的确定极为重要。本文通过综合根据原始DEM建立的山体阴影（地貌晕渲图），提取的等值线，3D显示的卫星影像三方面的因素，辅助判读，经过反复验证与试验，得到山谷线提取的最佳阈值为0.5102，山脊线提取的阈值为0.5614。

模型图如附录B所示，山谷线提取结果如图6（b1）所示，山脊线的提取结果如图6（b2）所示。

3.3 融合算法模型

融合算法实现思路是对水文分析法和曲率判别法的一种集成，利用两种算法进行加权评价，分值越高，代表两种算法均认可提取结果的准确，分值越低则与之相反。具体的实现过程是。

首先利用本文创建的等值线曲率计算模型或脚本，获取等值线曲率；同时，截取水文分析模型中的上半部分，获取蓄积栅格数据(流量)和正反地形。

然后基于等值线曲率，利用条件函数与阈值获取初始山脊线曲率、初始山谷线曲率；基于水文分析获取初始山脊线流量，山谷线流量

再者，基于上一步的初始山脊线、山谷线，使用分位数法将以上四个初始数据各自分为八个等级，等级越高，代表正确的可能性越大。

最后，将拉伸过的初始山脊线曲率与初始山脊线流量数据进行相加，将结果拉伸到0-1之间，越接近于1，结果是山脊线的可能性越大，山谷线与之相同。

模型图如附录C所示，山谷线提取结果如图6（c1）所示，山脊线的提取结果如图6（c2）所示。

4 山脊线山谷线提取结果的对比分析

4.1 水文分析法与曲率判别法的对比分析

本文采用了两种山脊线山谷线的提取算法，一是水文分析法，二是曲率判别法，两种算法各有自己的优缺点，本文分别对两种算法进行了建模实现与开发实现，并对两种算法的提取山脊线山谷线的结果分别进行对比分析。

首先建立了对比分析模型，将水文分析法与曲率判别法提取山脊线山谷线的结果进行分析，然后在本模型中将该对比分析功能进行了实现，对比分析模型建立的原理主要是地图代数运算与栅格局域运算，得到的分析结果栅格具有四个值，即为四个类型，其中

“0”值代表两种算法均为提取的区域，即并不是山脊线/山谷线区域；

“1”值代表曲率分析法提取的单独区域，此区域水文分析法未提取出山脊线/山谷线；

“2”值代表水文分析法提取的单独区域，此区域曲率分析法未提取出山脊线/山谷线；

图4 两种算法提取山脊线的结果对比

图5 两种算法提取山谷线的结果对比

利用水文分析法与曲率判别法分别提取山脊线的结果（栅格）对比如图4所示，提取山脊线的像元数如表5所示。利用水文分析法与曲率判别法分别提取山谷线的结果（栅格）对比如图5所示，提取山谷线的像元数如表6所示。

表5 两种算法提取山脊线的像元数的对比

值

说明

像元数

占比

0

其它区域

75729

89.43%

1

单独曲率判别法提取

2956

3.49%

2

单独水文分析法提取

2214

2.61%

3

公共区域

3777

4.46%

表6 两种算法提取山谷线的像元数的对比

值

说明

像元数

占比

0

其它区域

75769

89.48%

1

单独曲率判别法提取

2116

2.50%

2

单独水文分析法提取

3015

3.56%

3

公共区域

3775

4.46%

1）两种算法提取山脊线结果的对比分析

综合图4和表5的对比结果，可以看出两种算法提取山脊线均取得了很好的效果，在整体上都是正确的。相比而言，曲率判别法提取的范围较为全面，但连续性较差，水文分析法提取山脊线的连续性较好，但局部区域未能提取。

山脊线的提取结果受限于多方面的因素，首先是DEM的质量，其次是填洼、流向、流量、坡度、曲率等中间数据生成算法的影响，以流向计算为例，在GeoScene中采用的是D8流向法，即八方向法，该方法并非最优算法，但却是相对而言效率高精度高的算法。最后是人为干涉因素，包括各类阈值的设定，需要人为辅助判断。

2）两种算法提取山谷线结果的对比分析

综合图5和表6的对比结果，可以看出水文分析法与曲率判别法在山谷线的提取方面各有优缺点，首先，二者提取的结果都是正确的，但都有不足之处。水文分析法提取山谷线的连续性效果较好，在主山谷线处提取结果突出，正确率较高，且少有断裂现象发生。但是在分支处的提取结果较差，部分局部山谷线未得到提取。曲率判别法提取山谷线的整体性较好，局部的山谷线也得到了提取，但连续性较差，有断裂情况发生，局部存在错误点，需要进行众数滤波算法处理错误值。

山谷线的提取结果同样受限于多方面的因素，与山脊线的影响因素相同，首先是DEM的质量，其次是流向、流量提取算法，最后是各类阈值的设定。

4.2 融合算法与水文分析法、曲率判别法的对比分析

通过表5和表6的中的数据来看，水文分析法和曲率判别法单独提取区域占总栅格区域的2.50%-3.56%，将该部分数据通过图4和图5进行分析对比，可以看出：水文分析法提取的单独部分主要位于主体特征线，曲率判别法单独提取的部分大部分是水文分析法未提取的局部山脊线，因此均是正确的。总的来说，二者各有长短，优势互补。水文分析法的主体提取结果较好，连续性强，断裂线少，但局部提取结果差；而曲率判别法局部提取结果较好，但连续性较差；因此将二者算法进行融合，得到的结果相对而言是更为全面而准确的。

本文研究了两种算法的融合，利用水文分析提取的山脊、山谷线结果辅助曲率判别法的阈值设定，并将两种算法的提取结果进行叠加，利用众数滤波剔除局部错误值。两种方法融合后的提取结果如图6（c1）和图6（c2）所示，与图4和图5进行对比，可以明显看出，融合算法的提取结果较水文分析法和曲率判别法有所提高，无论是整体还是局部的山脊线与山谷线都得到了提取。三种结果的对比图如下：

图6 曲率判别法提取的山谷线（a1）、山脊线（a2）、水文分析法提取的山谷线（b1）、山脊线（b2）与融合算法提取的山谷线（c1）、山脊线（c2）

5 结论与展望

5.1 结论

本文基于GeoScene建模和脚本，实现了DEM构建、等高线绘制、等高线曲率计算、山脊线、山谷线的提取等功能。主要研究结果： 

1）完成了格网DEM的构建，包括5种空间插值算法：距离倒数乘方法、薄板样条函数法，趋势面分析法，自然邻域法、克里金法。五种算法中，自然邻域法和距离倒数乘方法对本溪关门山地区的适宜性最佳。

2）实现了山脊、山谷线的自动提取，完成了水文分析法与曲率判别法两种算法对山脊线与山谷线的自动提取。提取结果显示水文分析法连续性较好，但存在局部特征线未提取，曲率判别法整体性较好，提取结果较为全面，但存在特征线断裂情况。

3）研究了水文分析法与曲率判别法的融合，在两种算法的基础上实现了融合算法对山脊线与山谷线的自动提取。并对水文分析法、曲率判别法以及二者融合算法提取结果进行了对比分析，建立了对比分析模型，将分析结果量化显示。研究发现融合算法兼具二者之长，提取的山脊线与山谷线不仅连续性好，而且更加完整、全面。

4）提出了基于坡度与平面曲率计算栅格等高线曲率的方法，并实现了栅格形式等高线曲率的计算，同时也实现了矢量形式等高线曲率的计算。通过对比发现，栅格形式的等高线曲率计算速度较快，且精度与准确性更高。

5）制作了本溪关门山地区数字高程模型图、等高线专题图、山脊线专题图、山谷线专题图。

5.2 展望

本文对栅格空间插值方法、等高线的曲率计算方法、山脊线与山谷线的提取方法做了初步探索，存在以下问题有待进一步研究：在山脊线山谷线的提取方面，仍需验证更多的提取算法，并进行对比分析；水文分析法与等高线曲率判别法的融合策略有待进一步的加深。在今后的学习和工作中，需要不断探索，吸取更多经验，综合多领域的技术，继续加深研究，提高算法的效率与准确性，进一步优化“格网DEM的构建及山脊、山谷线提取模型”。

参考文献

[1] 刘淑琼,邹时林.基于格网DEM的地形特征线提取方法比较[J].测绘与空间地理信息,2015,38(2):85-86.

LIU Shuqiong, ZOU Shilin.Comparison of Terrain Feature Line Extraction Methods Based on Grid DEM[J].Surveying and Spatial Geography Information,2015,38(2):85-86.

[2] 茆德柱.TIN模型的构建方法研究:[D].南京:河海大学,2007.

Zhai Dezhu. Research on the construction method of TIN model: [D]. Nanjing: Hohai University, 2007.

[3] 黄杏元,马劲松.地理信息系统概论[M]. 北京: 高等教育出版社, 2001.

Huang Xingyuan, Ma Jinsong. Introduction to Geographic Information System [M]. Beijing: Higher Education Press, 2001.

[4] 周兴华,姚艺强,赵吉先.DEM内插方法与精度评定[J].测绘科学,2005,30(5):86-88.

ZHOU Xinghua, YAO Yiqiang, ZHAO Jixian. DEM interpolation method and accuracy evaluation[J]. Science of Surveying and Mapping, 2005, 30(5): 86-88.

[5] 冯桂.数字高程模型内插中等高线保形方法的研究:[D].武汉:武汉大学,2001.

Feng Gui. Research on the interpolation method of medium-high line in digital elevation model: [D]. Wuhan: Wuhan University, 2001.

[6] 焦晓双.机载激光雷达点云滤波算法与DEM内插方法研究:[D].太原:太原理工大学,2018.

Jiao Xiaoshuang. Research on airborne lidar filtering algorithm and DEM interpolation method for airborne laser radar: [D]. Taiyuan: Taiyuan University of Technology, 2018.

[7] Thibault D,Gold C M.Terrain Reconstruction from Contours by Skeleton Construction[J].GeoInformatica,2000,4(4):349-373.

[8] Moore I D,O"Loughlin E M,Burch G J.A contour‐based topographic model for hydrological and ecological applications[J].Earth Surface Processes and Landforms,1988,13(4):305-320.

[9] 黄培之,刘泽慧.地形断面高程极值法的理论研究[J].测绘通报,2005,(04):11-13.

HUANG Peizhi, LIU Zehui. Theoretical Study on the Elevation Extreme Method of Topographic Section[J]. Bulletin of Surveying and Mapping, 2005, (04): 11-13.

[10] 黄培之.提取山脊线和山谷线的一种新方法[J].武汉大学学报(信息科学版),2001,(03):247-252.

Huang Peizhi. A New Method for Extracting Ridge Lines and Valley Lines[J]. Journal of Wuhan University (Information Science Edition), 2001, (03): 247-252.

[11] 崔先国,汝绪伟,曲均浩.一种山脊线和山谷线提取的新方法[A].见：中国地理信息系统协会. 中国地理信息系统协会第八届年会论文集[C],中国地理信息系统协会:中国地理信息系统协会,2004.4.

Cui Xianguo, Qu Xuwei, Qu Junhao. A new method for extracting ridgeline and valley lines [A]. See: China Geographic Information System Association. Proceedings of the 8th Annual Meeting of China Geographic Information System Association [C ], China Geographic Information System Association: China Geographic Information System Association, 2004.4.

[12] 胡启明.地形特征线提取方法的研究进展[J].测绘与空间地理信息,2017,40(4):47-50+54.

Hu Qiming.Research progress on terrain feature line extraction method[J].Surveying and Spatial Geography Information,2017,40(4):47-50+54.

[13] Peucker T K,Douglas D H.Detection of surface-specific points by local parallel processing of discrete terrain elevation data[J].Computer Graphics and Image Processing,1975,4(4):375-387.

[14] O'callaghan J F.The extraction of drainage networks from digital elevation data[J].COMPUTER VISION, GRAPHICS, AND IMAGE PROCESSING,1984,28(3):323-344.

[15] Liang C,Mackay D S.A general model of watershed extraction and representation using globally optimal flow paths and up-slope contributing areas[J].International Journal of Geographical Information Science,2000,14(4):337-358.

[16] 贾蓉,刘春.基于机载LiDAR点云数据的复杂城市区域数字地面模型提取[J].遥感信息,2009,(05):3-7.

Jia Rong, Liu Chun. Extraction of Digital Urban Ground Model in Complex Urban Area Based on Airborne LiDAR Point Cloud Data[J]. Remote Sensing Information, 2009, (05): 3-7.

[17] Philip G M,Watson D F.A precise method for determining contoured surfaces[J].The APPEA Journal,1982,22(1):205-212.

[18] 陈娟,罗明良,张斌,(等).反距离加权插值参数变化对元谋冲沟DEM构建的影响[J].中国水土保持科学,2015,13(1):29-34.

CHEN Juan,LUO Mingliang,ZHANG Bin,(et al).Effects of inverse distance weighted interpolation parameters on DEM construction of Yuanmou gully[J].Chinese Journal of Soil and Water Conservation,2015,13(1):29-34.

[19] 张亮.矢量等高线曲率的计算方法研究:[D].陕西杨凌:中国科学院研究生院(教育部水土保持与生态环境研究中心),2014.

Zhang Liang. Research on the calculation method of vector contour curvature: [D]. Shaanxi Yangling: Graduate School of Chinese Academy of Sciences (Center for Soil and Water Conservation and Ecological Environment, Ministry of Education), 2014.

[20] 吴成永.艾比湖水底地形遥感研究:[D].乌鲁木齐:新疆师范大学,2010.

Wu Chengyong. Remote sensing research on the underwater topography of Aibi Lake: [D]. Urumqi: Xinjiang Normal University, 2010.

[21] 栗丹.基于DEM自动提取山脊线、山谷线方法研究[J].环球人文地理,2015,(2):77-77.

Li Dan. Research on automatic extraction of ridgeline and valley line based on DEM[J].Global Human Geography, 2015, (2): 77-77.

[22] Tarboton D G,Bras R L,Rodriguez-Iturbe I.On the extraction of channel networks from digital elevation data[J].Hydrological processes,1991,5(1):81-100.

[23] 牟乃夏,王海银,李丹,(等).GeoScene Engine 地理信息系统开发教程[M]. 北京: 测绘出版社, 2015.

Mou naixia, Wang Haiyin, Li Dan, (etc.). GeoScene Engine Geographic Information System Development Course [M]. Beijing: Surveying and Mapping Press, 2015.

[24] Chang K-T.地理信息系统导论[M]. 陈健飞等,译. 北京: 科学出版社, 2010.

Chang K-T. Introduction to Geographic Information System [M]. Chen Jianfei et al., Beijing. Science Press, 2010.