栗 业

（河北省地质矿产勘查开发局第二地质大队，河北 唐山 063000）

1 研究背景概述

1.1 研究背景

数字高程模型（Digital Elevation Model， DEM）是依据高程数据对地形曲面的数字化模拟，是地表形态数字化的量化模型。DEM是地理空间地形分析的核心数据，是新一代的地形图，同时也是GIS数据库中重要的组成部分。DEM与传统地形图相比有着巨大优势。首先，它可以在计算机中直接应用、计算，为各种辅助设计系统利用；其次，DEM是以数字形式存储，便于数据的录入、修改、更新、复制和管理，同时，可以方便地转换成其他形式的材料文件和产品；最后，DEM中含有传统地形图中所无法容纳和表达的垂直分布地物信息，即高程信息，所以DEM数据可以满足众多行业的发展需求。

1.2 研究意义

描述地貌形态的骨架线、控制线称为地性线，包括山脊线、山谷线、鞍部、凹地等线性物体。山谷线和山脊线构成了地势起伏变化的特征线，对地形地貌研究有重要意义。DEM数字高程模型中表达地形地貌的特征线，即地性线，其中包含了大量的地性线信息。基于DEM提取地性线的研究在区域水文分析、土木工程建设、地质勘探选定资源靶区、遥感影像自动配准等诸多领域具有重要的实践应用价值和科学意义。

2 国内外研究现状

多年来，国内外专家学者在数字高程模型中提取隐含着的山脊线和山谷线信息，及其如何将其快速自动提取，并应用于水文分析、地形分析、3D视域等多个研究领域方面，取得了一定的研究成果。

1984年，Yeoli提出了一个两步算法。首先，使用断面识别方法确定特征明显的点，存储于XMIN、YMIN、HMIN这3个矩阵中；然后从上述3个矩阵中去除已有山谷线上的点，用剩余的点中最大的点开始一条新山谷线的搜索。1991年，John Fairfield等在地表流水模型的基础上，推出了随机水流模型。该模型基本解决了DEM网格问题，但数据量成几十倍的增加，计算耗时多，且只适用于区域较小且边界较规则的范围；1995年，Amnon Meisels等提出了骨架法。该算法是提取等值线上曲率值较大的点，然后把它们连接起来作为特征线。这种算法所需的参数少，缺点是不能在较大凹陷地区使用。

近年来，国内学者也在此领域取得了重要突破。张维军应用道格拉斯-帕克提取山脊和山谷的特征点，并结合提取区域地形大致划分和汇水线的几何分析方法， 用概略的分水线和汇水线对特征点进行识别、归类、顺序提取各条山脊线和山谷线。姚慧敏等阐述了利用数学形态学的流域分割算法提取地性线的原理， 采用矢量栅格混合算法提取地性线的具体步骤， 最后实验结果表明， 其提出的算法是可行的， 提取的地性线有较好的连续性。黄培之研究了特定实验区域的山脊线和山谷线的几何特性或物理特性的单一方面设计的提取山脊线和山谷线的算法，随后又提出了一种基于地形表面流水分析与等高线几何分析结合的提取地性线的方法。该方法把等高线几何分析的方法与地形表面流水模拟分析的方法有机地结合起来， 能够克服各自所具有的弊端。

3 DEM与地性线的数学表达

3.1 DEM的数学表达及数据结构

DEM数字高程模型是地表形态的数字化表达。从狭义角度说DEM是区域地表正高的数字化表达模型；广义上来讲，它是时空中的地理对象表面正高高度的数字化表达模型。

数学意义上的DEM是表示区域A上的三维向量有限序列，函数表达为：

式中，（Xi，Yi）是平面坐标，Zi是（Xi，Yi）对应的高程。

目前，DEM主要有4种数据结构：①不规则离散点数据结构。不规则离散点是最简单的组织形式，它由不规则的离散点表示地表形态。这些点是通过测量获得的独立原始特征点数据。②规则格网（Grid）数据结构。在垂直、水平两个方向上以相等间隔逐行逐列记录每个网格的高程值称为规则格网DEM。为方便行列号和平面坐标的转换，规则网格中还标识了西南角坐标值、网格间距等。为压缩存储空间，常将规则网格每个单元值减去区域内最小高程值，再乘以一个常数是规则网格值全变为整数。规则格网DEM通常包含记录西南角起点坐标、坐标类型、行列数的头文件和记录行列高程矩阵的数据体。规则格网的编码结构有行程编码、块状编码、四叉树编码。③不规则三角网（TIN）数据结构。不规则三角网就是一系列互不重叠互不相交的三角形连接在一起表示地表形态。不规则三角网一般采用链表存储，由结点和三角形顶点两组记录组成。TIN是典型的矢量拓扑结构。TIN克服了规则格网的缺点，数据量精简，具有可变的分辨率，能表达复杂的地形。但该数据结构、算法都较复杂，管理困难，计算耗时长，适用于高精度小范围的地形建模。④等高线数据结构。等高线是最常用的地形表达形式，在地形图上高程相等的点所连成的闭合曲线即为等高线。它直观地反映了高程、坡度、坡向、山脉走向等基本地貌。等高线以坐标点对序列形式存储，在利用等高线重建地表模型时容易丢失细节信息，而且单条等高线无法反映地貌。

3.2 地性线的数学表达

3.2.1 山脊线

在实际地形地貌中条带状隆起的顶部线即山脊线。曲面上二阶方向导数取得极小值，且一阶导数为0的点是山脊点。山脊线则是相邻山脊点的连线。

假设Ta和Tb正负型反斜坡，Lab=Ta∩Tb是两斜坡的交线。在交线上存在一点P，其高程为Hp，则对任意iLab，存在CMi=Bi∩Ta∩Tb，使 得Hp-Hi≥0，得 出T=Ta∪Tb为山脊，Lab为山脊线或称分水线，Bi为过点i的平面，CMi为山脊线等高线。如果地形曲面T的表达式为Z = A00+A10X+ A01Y+A11XY+A20X2+A02Y2+……（A20<0；A02<0且A11≠0） 时，判断点P是否为山脊上一点，需满足任意下面条件：坡度α=0；最大曲率Cmax>0；且最小曲率Cmin=0或坡度α>0；且断面曲率Cs>0。

3.2.2 山谷线

山脊线和山谷线是两个完全对应相反的概念。曲面上二阶方向导数取得极大值，且一阶导数为0的点是山谷点。山谷线则是相邻山脊点的连线。山谷线收集两侧的水，属于汇水区域，易形成河流。山谷线是纯粹的局部地貌概念，因为它在端点处的一小段范围内可能不会形成水系。由于它的局部性质，它比河流更适合，更便于计算机提取。所有我们近似地把它当成河流的一部分。

假设Ta和Tb正负型反斜坡，Lab=Ta∩Tb是两斜坡的交线。在交线上存在一点P，其高程为Hp，则对任意iLab，存在CMi=Bi∩Ta∩Tb，使 得Hp-Hi≤0，得 出T=TaTb为 谷底，Lab为山谷线或称汇水线，Bi为过点i的平面，CMi为山谷线等高线。如果地形曲面T的表达式为Z=A00+A10X+ A01Y+A11XY+A20X2+A02Y2+……（A20>0； A02>0且A11≠0） 时，判断点P是否为山谷上一点，需满足任意下面条件：坡度α=0；最大曲率Cmax=0；且最小曲率Cmin<0或坡度α>0；且断面曲率Cs<0。

4 基于DEM提取地性线的数学算法

4.1 基于图像处理的算法

DEM实质是一种栅格数据形式，网格大小即为像素大小，高程值为灰度值。DEM数据可以按数字图像处理的方法来设计算子。在DEM中高程值越大，灰度越大，山脊线上的点比其他地方的点都要亮。按照各种滤波算子就能进行边缘提取，将地形特征线提取出来。Chakrevavanieh提出用高斯算子从DEM上自动提取山脊线和山谷线上的点，并将其用于规则格网的数据压缩。唐心红等在提取地性线前利用小波分析进行多尺度表示，在多尺度DEM中提取。余生晨利用傅里叶变换在频域中计算山脊线和山谷线，但该方法计算量大。基于图像处理的技术还有以下两个方面的难点：排除DEM中的噪声影响；所提取特征点连接成线的算法问题。

4.2 基于地表几何形态的算法

断面极值法是地表几何形态分析的基础算法。该算法的基本思路是：在DEM垂直和水平横截面上的极值，极大值就是分水点，极小值是汇水点；然后根据相应的准则划分极值点归属形成地性线。地表几何形态分析算法由于只判断两个断面，提取过程中会出现遗漏现象。而且阈值过大会丢失特征点，形成的地性线也就断裂较多，长度较短，阈值过小，特征点会出现过度提取现象，地性线冗余。

Fowler改进该算法，提出用4栅格扫面标记山脊点山谷点，这4栅格单元分别为（i， j）、（i， j+1）、（i+1， j）、（i+1， j+1），相比于8邻域算法效率高而且不会遗漏山脊点和山谷点，缺陷是如果DEM采样点密集会提取多余的特征点。Jenson认为任意一栅格与相邻8单元格形成4对位置，至少一对高程都比该点高程值大，则该点可能是山谷点，如果至少有一对高程点比该点小，该点就可能为山脊点。这种改进算法易于编程实现，但产生地性线的断裂较多。针对特征点丢失现象，黄培之提出再增加一组对角线上断面的方法来克服这一缺陷，研究实验证明该法有效，但同时增加了特征点噪声。

4.3 基于地形表面流水物理模型分析算法

根据水往低处流，最终顺谷而流，汇集成水系网的自然规律，Chakrevavanieh率先提出模拟地表流水物理模型。D8法是该模型中的最常用算法：假设只有一个单一的网格的水流进八个相邻的网格单元。水流方向以最陡坡度法来确定，即在3×3的DEM栅格上，计算中心栅格与各相邻栅格间的距离权落差，取距离权落差最大的栅格为中心栅格的流出栅格。以数值表示每个单元的流向，变化范围是1~255。然后汇流栅格图由流向栅格图生成。汇流栅格上每个单元的值代表上游汇流区汇入该单元的流入路径数较大者。从流域汇流栅格图中可以轻易地提取流域各种参数特征值。在实际应用时需要解决以下问题：①水道起始点位置的确定；②伪负地形的识别及填充；③凹陷与平坦处水流方向的设定。而且基于地表水流模拟方法存在两个缺陷：①洼地的填充和平坦地区水流方向的设定十分麻烦；②该方法所计算的汇水量与高程有关，计算的结果为高程值大的区域汇水量小，高程值小的区域汇水量大，由于水流遵从高处连续性的向下流动的自然特征，这就使得在一个河道水系或山谷线上的相对位置较高的某一点因其汇水量较小而被遗漏，而再非山谷处处于低位置的某点因其汇水量较大而被误判为山谷线上的点。在后来的研究中，许多学者提出利用平滑DEM，垫高填平凹陷区，标识洼地集水区的方法解决上述问题，并取得了良好效果。

4.4 基于地形表面几何形态分析和地表流水物理模型分析的算法

由于地表几何形态分析和地表流水分析模型的算法有一些不足之处，所以专家们提出了将两种模式结合提取地形线。这种算法的基本思路是首先按照地表流水模型在稀疏的DEM中概略提取地性线；然后对周围地域进行几何地理形态分析，在精确定位地性线。该方法的关键是求出概略地性线与DEM网格的交点，在该点周围的小范围内进行几何分析，找出正交方向地形断面上的高程极值点，正是这一点为精确点。这种综合算法在特征点提取阈值的选取、DEM网格大小的确定方面还有待进一步研究。

4.5 平面曲率与坡形组合法

求取DEM平面曲率和地表的正负地形，正地形上曲率的大值是山脊，负地形上曲率的大值是山谷。这种方法中平面曲率的大小可以用来调教提取宽度，方法简便效果好，但提取的图形是栅格形式，在转为矢量线性方面还需进一步研究改进。

5 结语

综上所述，基于DEM提取地性线是项复杂的工作，目前提取算法相对成熟，但不同算法各有利弊，针对实际工作的不同地域地形和收集到的DEM数据质量，应选择最佳的算法进行提取。