DDEEEMMM,构建及地形特征线提取研究
李胜天，徐贵兴
（江西省地质局地理信息工程大队，江西南昌 330001）
数字高程模型（digital elevation model，DEM）的数据结构简单、容易管理、便于存贮，许多国家的DEM 数据都是以规则格网的数据矩阵形式提供的，因此，基于规则格网DEM 提取等高线、山脊线、山谷线一直是研究人员的研究方向[1]。山脊线、山谷线是地形骨架结构的重要组成部分[2]，常应用于地貌表达以及三维地形分析。建立格网DEM 的空间内插方法可以分为分块内插法、整体内插法、逐点内插法3 类。其建立一般包括数据取样、数据内插和数据精度分析等步骤[3]。整体内插法的原理简单，适用于总体预测，但不适用于复杂地形，对于复杂地形，常使用分块内插法[4-5]。逐点内插法是应用最为广泛的方法，包括距离倒数乘方法、克里金法，在样本点分布均匀且数量较多时均可取得较好的结果。当实际地形的坡度较大时，采用距离倒数乘方法较为适宜[6]。提取山脊线、山谷线常用的算法有3 类。一是几何分析法，它利用地表几何形态的实际规律，通过一定的算法（如曲率判别法、断面极值法和骨架化法）提取地形特征点，其中断面极值法与曲率判别法均是利用极值提取地形特征点。断面极值法一般采用两个方向的断面，并不能提取出所有地形特征点[7-9]。骨架化法假设山脊线、山谷线两侧地形对称，提取中心轴线，与实际不符，有很大程度的近似性，因此所得结果不尽如人意[10-11]。二是图像处理法，它将DEM 数据作为灰度图像处理，通过数学形态学的“序贯减薄”运算得到骨架线[12]，该方法可利用移动窗口法提取特征点，但是将特征点连成线非常困难。三是水文分析法，它利用地表水流的客观规律对山脊线、山谷线进行提取，常用D8 流向法计算流向，在带有汇聚型河流的多山地形能生成较好的结果。但在平坦地区不易确定流向，不适用于含泛滥平原和湿地的多变地形。对格网DEM 构建常用的5 种算法，即距离倒数乘方法、自然邻域法、趋势面分析法、克里金法、样条函数法进行研究，并分析它们在本溪关门山地区的适用性与准确性，选出最适合的插值方法——距离倒数乘方法和自然邻域法用于格网DEM 的构建。基于格网DEM，利用水文分析法和曲率判别法两种方法提取山脊线、山谷线，将提取结果进行对比分析，并研究和实现了两种山脊线、山谷线提取算法的融合。
本溪关门山地区，采用不同的空间插值方法提取山脊线、山谷线，通过GeoScene Pro 建模和二次开发的形式，分别完成格网DEM 的构建、等高线的绘制、等高线曲率的计算、山脊线和山谷线的自动化提取。同时分析不同插值算法的适宜性，以及不同山脊线、山谷线提取算法的优缺点。其技术路线步骤如下：
第一步：格网DEM 的构建，利用5 种空间插值算法进行栅格插值，并通过交叉验证选出最合适的方法用于本溪关门山地区格网DEM的构建。
第二步：等高线的自动绘制，基于格网DEM 数据，生成等高线。
第三步：等高线曲率的计算，基于等高线数据，计算等高线曲率。
第四步：曲率判别法的实现，基于等高线曲率，完成山谷线及山脊线的提取。
第五步：水文分析法的实现，基于水文分析法，完成山谷线及山脊线的提取。
第六步：水文分析法与曲率判别法提取结果的对比分析以及两种算法的融合。
空间插值算法构建格网DEM有5种方法，包括距离倒数乘方法、克里金法、薄板样条函数法、趋势面分析、自然邻域法。在此基础上，对插值结果进行对比分析。在进行插值之前，首先将样本点分为两部分，其一为训练点，其二为测试点，样本点共有37 660 个，训练点为1 000 个，测试点为36 660 个，其中训练点用于检测插值质量，测试点用于格网DEM 的构建。在插值成功之后，将5种插值结果进行3次对比，即将样本点随机分配3次，进行独立试验，互不干扰。每次对比的测试点与训练点的分配比例相同，但分布各不相同，均是随机的。根据测试点的5 种插值算法的插值计算结果，将训练点所在位置的值分别提取至训练点属性表，共分别提取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，插值结果不如距离倒数乘方法和样条函数法，但差距不大，也能反映真实地形。趋势面分析法插值结果极差，标准差为123.728 0，偏差较大，只能反映整体趋势，不能反映真实地形，不适用于DEM的建立。
图1 构建网格的5种空间插值算法插值结果Figure 1. Interpolation results of five spatial interpolation algorithms for grid construction
表1 三次插值结果对比Table 1. comparison of the results of cubic interpolation
3.1 曲率计算模型
等高线曲率与GeoScene Pro 中的曲率有根本上的不同，等高线曲率是指地表面上通过某一点的等高面与地表的交线在地形表面水平方向上的凹凸性。
地形表面S的曲面方程为z=f(x,y)，点P为曲面S上任意一点，平面α为过点P的水平面，α与S的交线为c1过点P的等高线；
a为等高线在c1点P的切向量，b为点P的坡向向量。b⊥a平面β、γ，分别过a、b垂直于水平面XOY并与曲面S相交于曲线c2、c3。曲面S：z=f(x,y)的各阶偏导数用以下符号表示：
GeoScene Pro 中提供了平面曲率的计算工具，该曲率是垂直于坡度方向的曲率。从图2中的数学关系来看，平面曲率就是曲线c2在点P的二阶导数，它的计算公式为：
在GeoScene Pro 中，用最大平均值法计算坡度值。GeoScene Pro中的坡度可利用式（1）中的变量进行表达，其计算公式如下：
地形曲面S：z=f(x,y)的任意一条等高线方程可表示为f(x,y)=c（c为表示任意高程值的常数），由曲率计算公式可得等高线曲率计算公式为：
通过综合公式（2）~（4），可以发现平面曲率、等值线曲率、坡度三者之间具有一定的关系，通过平面曲率与坡度便可计算等值线曲率栅格。
推导出的公式如下：
曲率判别法的实现思路是：首先获取等值线曲率栅格，其次基于曲率栅格获取地形的特征点，通过一定的分类阈值提取山脊点与山谷点。
3.2 水文分析模型
水文分析法的基本原理是：对于山脊线而言，它也是分水线，分水线所在的栅格区域在一定程度上不存在水流，即流量累积栅格的对应值为零。通过地表径流模拟计算可以得到流量累积栅格，并提取零累积值区域，即分水线所在区域，也就得到了山脊线[13]。若以某一相对高程面为对称基准面，求原始地形的对称地形，则山谷线与山脊线位置也互换。因此，可以通过建立的反地形提取模型，求出对称地形（反地形），那么，正地形山谷线就可以在反地形中利用提取山脊线的方法进行提取[14]。
（1）已填洼DEM。水文分析法中已填洼DEM 是指不存在小洼地，大多数小洼地是DEM 生成过程中带来的数据错误所致，如果未进行填充，则生成的水系网络可能会呈现不连续性。常用的除去洼地的方法是把洼地像元值加高至周围最低像元值[15]。填洼操作不仅能去除小洼地，同时也可以去除错误的峰。峰是一种伪像元，其高程高于所预期的高程，是反地形中的小洼地。
（2）反地形DEM。传统的水文分析法在提取山谷线过程中，反地形的计算是基于原始DEM 数据的，利用已填洼DEM 数据提取反地形，从而消除反地形中的峰和小洼地，进而消除山谷线提取结果的局部错误值。反地形的计算公式为：
式中：Abs为绝对值运算，Dem为已填洼DEM 数据；
Dmax与Dmin分别为原始DEM 的最大高程值和最小高程值。模型中山脊线、山谷线提取阈值的确定和反地形的计算极为重要。对汇流累积量零值数据进行重分类时，属性值越接近1 越可能是山脊线或山谷线的位置，分界阈值的确定尤为重要。根据原始DEM 建立的山体阴影（地貌晕渲图），提取等值线，通过综合3D 显示的卫星影像三方面的因素，辅助判读，经过反复验证与试验，得到山谷线提取的最佳阈值为0.510 2，山脊线提取的阈值为0.561 4。
3.3 融合算法模型
融合算法实现思路是对水文分析法和曲率判别法的一种集成，利用两种算法进行加权评价，分值越高，代表两种算法均认可提取结果的准确，分值越低，则与之相反。具体的实现过程是：首先利用创建的等值线曲率计算模型或脚本，获取等值线曲率；
同时，截取水文分析模型中的上半部分，获取蓄积栅格数据（流量）和正反地形；
然后基于等值线曲率，利用条件函数与阈值获取初始山脊线曲率、初始山谷线曲率。基于水文分析获取初始山脊线流量、山谷线流量；
再者，基于上一步的初始山脊线、山谷线，使用分位数法将以上初始数据各自分为8个等级，等级越高，代表正确的可能性越大。最后，将拉伸过的初始山脊线曲率与初始山脊线流量数据进行相加，将结果拉伸到0～1，越接近于1，结果是山脊线的可能性越大。山谷线与之相同。
4.1 水文分析法与曲率判别法的对比分析
采用两种山脊线、山谷线的提取算法，一是水文分析法，二是曲率判别法。运用两种算法进行了建模实现与开发实现，并对两种算法的提取山脊线、山谷线的结果分别进行对比分析，两种算法各有优缺点。首先建立了对比分析模型，将水文分析法与曲率判别法提取山脊线、山谷线的结果进行分析，然后在本模型中将该对比分析功能进行了实现。对比分析模型建立的原理主要是地图代数运算与栅格局域运算，得到的分析结果栅格具有4 个值，即为4 个类型，其中“0”值代表两种算法均为提取的区域，而并不是山脊线、山谷线区域；
“1”值代表曲率分析法提取的单独区域，此区域水文分析法未提取出山脊线、山谷线；
“2”值代表水文分析法提取的单独区域，此区域曲率分析法未提取出山脊线、山谷线。
利用水文分析法与曲率判别法分别提取山脊线的结果（栅格）对比如图2 所示，提取山脊线的像元数如表2所示。利用水文分析法与曲率判别法分别提取山谷线的结果（栅格）对比如图3 所示，提取山谷线的像元数如表3所示。
表2 两种算法提取山脊线的像元数对比Table 2. Comparison of pixel elements of ridge lines extracted by two algorithms
表3 两种算法提取山谷线的像元数对比Table 3. Comparison of pixel elements of valley lines extracted by two algorithms
图2 两种算法提取山脊线的结果对比Figure 2. Comparison of ridge lines extracted by two algorithms
图3 两种算法提取山谷线的结果对比Figure 3. Comparison of valley lines extracted by two algorithms
4.1.1 两种算法提取山脊线结果的对比分析
综合图2 和表2 的对比结果，可以看出两种算法提取山脊线均取得了很好的效果，在整体上都是正确的。相比而言，曲率判别法提取的范围较为全面，但连续性较差，水文分析法提取山脊线的连续性较好，但局部区域未能提取。山脊线的提取结果受限于多方面的因素，首先是DEM 的质量，其次是填洼、流向、流量、坡度、曲率等中间数据生成算法的影响。以流向计算为例，在GeoScene Pro中采用的是D8流向法，即八方向法，该方法并非最优算法，但相对而言却是效率高、精度高的算法。最后是人为干涉因素，包括各类阈值的设定，需要人为辅助判断。
4.1.2 两种算法提取山谷线结果的对比分析
综合图3 和表3 的对比结果，可以看出水文分析法与曲率判别法在山谷线的提取方面各有优缺点，首先，二者提取的结果都是正确的，但都有不足之处。水文分析法提取山谷线的连续性效果较好，在主山谷线处提取结果突出，正确率较高，且少有断裂现象发生，但是在分支处的提取结果较差，部分局部山谷线未得到提取。曲率判别法提取山谷线的整体性较好，局部的山谷线也得到了提取，但连续性较差，有断裂情况发生，局部存在错误点，需要进行众数滤波算法处理错误值。山谷线的提取结果同样受限于多方面的因素，与山脊线的影响因素相同，首先是DEM的质量，其次是流向、流量提取算法，最后是各类阈值的设定。
4.2 融合算法与水文分析法、曲率判别法的对比分析
通过表2 和表3 中的数据来看，水文分析法和曲率判别法单独提取区域占总栅格区域的2.50%～3.56%，将该部分数据通过图3 和图4 进行分析对比，可以看出：水文分析法提取的单独部分主要位于主体特征线，曲率判别法单独提取的部分大部分是水文分析法未提取的局部山脊线，因此均是正确的。总的来说，二者各有长短，优势互补。水文分析法的主体提取结果较好，连续性强，断裂线少，但局部提取结果差；
而曲率判别法局部提取结果较好，但连续性较差。因此将二者算法进行融合，得到的结果相对而言是更为全面而准确的。研究了两种算法的融合，利用水文分析提取的山脊线、山谷线结果辅助曲率判别法的阈值设定，并将两种算法的提取结果进行叠加，利用众数滤波剔除局部错误值。由两种方法融合后提取结果可以明显看出，融合算法的提取结果较水文分析法和曲率判别法有所提高，无论是整体还是局部的山脊线与山谷线都得到了提取。3 种结果的对比如图4所示。
图4 曲率判别法、水文分析法、融合算法分别提取的山谷线、山脊线对比图Figure 4. Comparison of valley lines and ridgelines extracted by curvature discrimination method, hydrologic analysis method and fusion algorithm
5.1 结论
基于GeoScene Pro 建模和脚本，实现了DEM 构建、等高线绘制、等高线曲率计算、山脊线与山谷线的提取等功能。主要研究结果：
（1）完成了格网DEM 的构建，包括5 种空间插值算法：距离倒数乘方法、薄板样条函数法、趋势面分析法、自然邻域法、克里金法。5 种算法中，自然邻域法和距离倒数乘方法对本溪关门山地区的适宜性最佳。
（2）实现了山脊线、山谷线的自动提取，完成了水文分析法与曲率判别法两种算法对山脊线与山谷线的自动提取。提取结果显示水文分析法连续性较好，但存在局部特征线未提取，曲率判别法整体性较好，提取结果较为全面，但存在特征线断裂情况。
（3）研究了水文分析法与曲率判别法的融合，在两种算法的基础上实现了融合算法对山脊线与山谷线的自动提取。并对水文分析法、曲率判别法以及二者融合算法提取结果进行了对比分析，建立了对比分析模型，将分析结果量化显示。研究发现融合算法兼具二者之长，提取的山脊线与山谷线不仅连续性好，而且更加完整、全面。
（4）提出了基于坡度与平面曲率计算栅格等高线曲率的方法，并实现了栅格形式等高线曲率的计算，同时也实现了矢量形式等高线曲率的计算。通过对比发现，栅格形式的等高线曲率计算速度较快，且精度与准确性更高。
5.2 展望
对栅格空间插值方法、等高线的曲率计算方法、山脊线与山谷线的提取方法做了初步探索，存在以下问题有待进一步研究：在山脊线、山谷线的提取方面，仍需验证更多的提取算法，并进行对比分析；
水文分析法与等高线曲率判别法的融合策略有待进一步的加深。在今后的学习和工作中，需要不断探索，吸取更多经验，综合多领域的技术，继续加深研究，提高算法的效率与准确性，进一步优化“格网DEM 的构建及山脊线、山谷线提取模型”。