1、第 32 卷第 4 期2022 年 12 月安徽地质Geology of AnhuiDec.2022Vol.32 No.4文章编号：1005-6157（2022）04-0引言数字高程模型（digital elevation model，DEM）的数据结构简单、容易管理、便于存贮，许多国家的DEM数据都是以规则格网的数据矩阵形式提供的，因此，基于规则格网DEM提取等高线、山脊线、山谷线一直是研究人员的研究方向1。山脊线、山谷线是地形骨架结构的重要组成部分2，常应用于地貌表达以及三维地形分析。建立格网DEM的空间内插方法可以分为分块内插法、整体内插法、逐点内插法3类。其建立一般包括数据取样、数据内
2、插和数据精度分析等步骤3。整体内插法的原理简单，适用于总体预测，但不适用于复杂地形，对于复杂地形，常使用分块内插法4-5。逐点内插法是应用最为广泛的方法，包括距离倒数乘方法、克里金法，在样本点分布均匀且数量较多时均可取得较好的结果。当实际地形的坡度较大时，采用距离倒数乘方法较为适宜6。提取山脊线、山谷线常用的算法有3类。一是几何分析法，它利用地表几何形态的实际规律，通过一定的算法（如曲率判别法、断面极值法和骨架化法）提取地形特征点，其中断面极值法与曲率判别法均是利用极值提取地形特征点。断面极值法一般采用两个方向的断面，并不能提取出所有地形特征点7-9。骨架化法假设山脊线、山谷线两侧地形对称，提
3、取中心轴线，与实际不符，有很大程度的近似性，因此所得结果不尽如人意10-11。二是图像处理法，它将DEM数据作为灰度图像处理，通过数学形态学的“序贯减薄”运算得到骨架线12，该方法可利用移动窗口法提取特征点，但是将特征点连成线非常困难。三是水文分析法，它利用地表水流的客观规律对山脊线、山谷线进行提取，常用D8流向法计算流向，在带有汇聚型河流的多山地形能生成较好的结果。但在平坦地区不易确定流向，不适用于含泛滥平原和湿地的多变地形。对格网DEM构建常用的5种算法，即距离倒数乘方法、自然邻域法、趋势面分析法、克里金法、样条函数法进行研究，并分析它们在本溪关门山地区的适用性与准确性，选出最适合的插值方
4、法距离倒数乘方法和自然邻域法用于格网DEM的构建。基于格网DEM，利用水文分析法和曲率判别法两种方法提取山脊线、山谷线，将提取结果进行对比分析，并研究和实现了两种山脊线、山谷线提取算法的融合。收稿日期：2022-3-19作者简介：李胜天（1982），男，江西萍乡人，高级工程师，主要从事不动产测绘和航空摄影测量、地理信息工程、数据处理等工作。E-mail：摘要：格网DEM是数字地面模型的重要数据表现形式，基于格网DEM提取的山脊线、山谷线是地形骨架结构的表达与数字地形分析的重要方面。通过比较距离倒数乘方法、自然邻域法、克里金法、薄板样条函数法、趋势面分析法5种空间内插方法在本溪关门山地区的适用性
5、，选出最优法用于格网DEM的构建。基于格网DEM，对比运用水文分析法和曲率判别法两种方法提取山脊线、山谷线的结果，并研究两种提取算法的融合效果。以GeoScene Pro建模和脚本的形式分别实现了DEM构建、等高线曲率计算、山脊线与山谷线提取；制作了本溪关门山地区数字高程模型图、等值线专题图、山脊线与山谷线专题图等。提出了坡度、平面曲率与等高线曲率的转换方法；对水文分析法和曲率判别法做出了改进，并提出了基于水文分析法和曲率判别法的融合算法，进一步提高了山脊线、山谷线提取的准确率。融合算法兼具二者之长，提取结果较为全面，整体性较好。关键词：DEM；数字高程模型；空间差值；地形特征中图分类号：P2
6、58文献标志码：ADEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEMDEM 构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形
7、特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究构建及地形特征线提取研究李胜天，徐贵兴（江西省地质局地理信息工程大队，江西南昌330001）365-6安徽地质2022年1格网DEM构建及山脊线、山谷线提取方法与技术路线本溪关门山地区，采用不同的空间插值方法提取山脊线、山谷线，通过GeoScene Pro建模和二次开发的形
8、式，分别完成格网DEM的构建、等高线的绘制、等高线曲率的计算、山脊线和山谷线的自动化提取。同时分析不同插值算法的适宜性，以及不同山脊线、山谷线提取算法的优缺点。其技术路线步骤如下：第一步：格网DEM的构建，利用5种空间插值算法进行栅格插值，并通过交叉验证选出最合适的方法用于本溪关门山地区格网DEM的构建。第二步：等高线的自动绘制，基于格网 DEM 数据，生成等高线。第三步：等高线曲率的计算，基于等高线数据，计算等高线曲率。第四步：曲率判别法的实现，基于等高线曲率，完成山谷线及山脊线的提取。第五步：水文分析法的实现，基于水文分析法，完成山谷线及山脊线的提取。第六步：水文分析法与曲率判别法提取结果
9、的对比分析以及两种算法的融合。2空间插值算法实现与分析空间插值算法构建格网DEM有5种方法，包括距离倒数乘方法、克里金法、薄板样条函数法、趋势面分析、自然邻域法。在此基础上，对插值结果进行对比分析。在进行插值之前，首先将样本点分为两部分，其一为训练点，其二为测试点，样本点共有37 660个，训练点为1 000个，测试点为36 660个，其中训练点用于检测插值质量，测试点用于格网DEM的构建。在插值成功之后，将5种插值结果进行3次对比，即将样本点随机分配3次，进行独立试验，互不干扰。每次对比的测试点与训练点的分配比例相同，但分布各不相同，均是随机的。根据测试点的5种插值算法的插值计算结果，将训练
10、点所在位置的值分别提取至训练点属性表，共分别提取5次。之后便可计算每种插值算法在每个训练点处的偏差。基于此属性进行统计分析，分别计算均值、最大值、最小值、标准差、方差。统计表见表1。从表1、图1可以看出，本实验区（本溪关门山地区）应用的 5 种插值算法中，距离倒数乘方法、自然邻域法的均方根误差为 4.533 0 和 3.339 3，均值为-0.162 2和-0.073 8，插值结果最好，最能反映真实地形。克里金法和样条函数法稍差，均方根误差为5.603 2和5.748 9，均值分别为-0.293 4、-0.238 1，插值结果不如距离倒数乘方法和样条函数法，但差距不大，也能反映真实地形。趋势面
11、分析法插值结果极差，标准差为 123.728 0，偏差较大，只能反映整体趋势，不能反映真实地形，不适用于DEM的建立。3山脊线、山谷线提取模型的建立3.1 曲率计算模型等高线曲率与GeoScene Pro中的曲率有根本上表1三次插值结果对比Table parison of the results of cubic interpolation条目均方根误差标准值均值最大值最小值距离倒数乘方法4.533 07.712 9-0.162 243.701 5-67.685 6自然邻域法3.339 35.739 8-0.073 837.830 7-108.141 8克里金法5.603 29.487 5-0
12、.293 442.132 5-69.869 9薄板样条函数法5.748 99.680 0-0.238 184.812 5-484.689 3趋势面法74.013 7123.728 0-1.442 0291.429 0-419.482 7（a）距离倒数乘方法；（b）克里金法；（c）样条函数法；（d）趋势面法；（e）自然邻域法图1构建网格的5种空间插值算法插值结果Figure 1.Interpolation results of five spatial interpolation algorithms for grid construction366第32卷第4期的不同，等高线曲率是指地表面上通
13、过某一点的等高面与地表的交线在地形表面水平方向上的凹凸性。地形表面S的曲面方程为z=f(x,y)，点P为曲面S上任意一点，平面为过点P的水平面，与S的交线为c1过点P的等高线；a为等高线在c1点P的切向量，b为点P的坡向向量。ba平面、，分别过a、b垂直于水平面XOY并与曲面S相交于曲线c2、c3。曲面S：z=f(x,y)的各阶偏导数用以下符号表示：p=fx，q=fy，r=2fx2，t=2fy2，s=2fxy（1）GeoScene Pro中提供了平面曲率的计算工具，该曲率是垂直于坡度方向的曲率。从图2中的数学关系来看，平面曲率就是曲线c2在点P的二阶导数，它的计算公式为：c_plan=-q2r
14、-2pqs+p2tp2+q2（2）在 GeoScene Pro 中，用最大平均值法计算坡度值。GeoScene Pro中的坡度可利用式（1）中的变量进行表达，其计算公式如下：slope_Degree=57.29578 arctan()p2+q2（3）地形曲面S：z=f()x,y的任意一条等高线方程可表示为f()x,y=c（c为表示任意高程值的常数），由曲率计算公式k=|y|/()1+y23可得等高线曲率计算公式为：c_cont=-q2r-2pqs+p2t()p2+q23 2（4）通过综合公式（2）（4），可以发现平面曲率、等值线曲率、坡度三者之间具有一定的关系，通过平面曲率与坡度便可计算等值线
15、曲率栅格。推导出的公式如下：c_cont=c_plantan(Slope_Degree/57.29578)（5）曲率判别法的实现思路是：首先获取等值线曲率栅格，其次基于曲率栅格获取地形的特征点，通过一定的分类阈值提取山脊点与山谷点。3.2 水文分析模型水文分析法的基本原理是：对于山脊线而言，它也是分水线，分水线所在的栅格区域在一定程度上不存在水流，即流量累积栅格的对应值为零。通过地表径流模拟计算可以得到流量累积栅格，并提取零累积值区域，即分水线所在区域，也就得到了山脊线13。若以某一相对高程面为对称基准面，求原始地形的对称地形，则山谷线与山脊线位置也互换。因此，可以通过建立的反地形提取模型，求
16、出对称地形（反地形），那么，正地形山谷线就可以在反地形中利用提取山脊线的方法进行提取14。（1）已填洼DEM。水文分析法中已填洼DEM是指不存在小洼地，大多数小洼地是DEM生成过程中带来的数据错误所致，如果未进行填充，则生成的水系网络可能会呈现不连续性。常用的除去洼地的方法是把洼地像元值加高至周围最低像元值15。填洼操作不仅能去除小洼地，同时也可以去除错误的峰。峰是一种伪像元，其高程高于所预期的高程，是反地形中的小洼地。（2）反地形 DEM。传统的水文分析法在提取山谷线过程中，反地形的计算是基于原始DEM数据的，利用已填洼DEM数据提取反地形，从而消除反地形中的峰和小洼地，进而消除山谷线提取结果的局部错误值。反地形的计算公式为：FanDem=Abs(Dem-Dmax-Dmin)（6）式中：Abs 为绝对值运算，Dem为已填洼 DEM 数据；Dmax与Dmin分别为原始DEM的最大高程值和最小高程值。模型中山脊线、山谷线提取阈值的确定和反地形的计算极为重要。对汇流累积量零值数据进行重分类时，属性值越接近1越可能是山脊线或山谷线的位置，分界阈值的确定尤为重要。根据原始DEM建立的山体阴影（