摘要
针对传统规则格网数字高程模型地形特征线提取方法中存在阈值难以定量调控、连接方式无法自适应调整以及地形特征线类型不完整的问题，将滚动球变换模型应用于不规则三角网数字水深模型（triangulated irregular network digital depth model，TIN-DDM）地形特征线的自动提取，在构建临界滚动球半径关联的地形特征点定量识别判定准则基础上，引入地形形态边界点概念，采用逆向工程的建模思路，建立了以剖分单元为基础的地形特征线自动提取模型，结合地形类型判定准则的多尺度表达特性及顾及水深数值的地形特征优化模型，提出了可多尺度表达且类型完整的地形特征线自动提取方法。试验结果表明：相比于经典地表流水模拟方法，该方法可实现完整、连续且细分的TIN-DDM地形特征线自动提取及多尺度表达，且提取的地形特征线具有更高的地形重构精度。
关键词
Abstract
In view of the problems in the traditional regular grid digital elevation model terrain feature line extraction method, such as the threshold is difficult to be adjusted quantitatively, the connection mode cannot be adjusted adaptively, and the type of terrain feature line is incomplete, an automatic extraction method of TIN-DDM (triangulated irregular network digital depth model) terrain feature line based on rolling ball transformation model was proposed. Based on the judgment criterion of a quantitative identification for terrain feature points associated with the critical rolling sphere radius, the concept of terrain shape boundary points was introduced, and the automatic extraction model of terrain feature lines based on subdivision unit was established according to the modeling idea of reverse engineering. Combined with the multi-scale expression characteristics of the terrain type judgment criteria and the terrain feature optimization model considering water depth values, an automatic extraction method of terrain feature lines that can be expressed at multiple scales and complete types was proposed. The experimental results show that compared with the classical surface water simulation method, this method can automatically extract complete, continuous, subdivided and multi-scale TIN-DDM terrain feature lines, and the generated terrain feature lines have higher terrain reconstruction accuracy.
地形特征线是数字水深模型（digital depth model，DDM）的结构化框架，是建立局部地形特征分析与整体地形形态控制的关联纽带。一方面，地形特征线可实现诸如山脊、山谷、鞍部、平台等地形特征的识别与判定，可为无定位条件下的潜艇、水下潜航器等武器平台提供地形匹配导航信息支持；另一方面，地形特征线可实现骨架地形与细部地形的有效区分，已成为DDM数据压缩及地形重构的关键条件约束。此外，由于山脊线与山谷线等不同类型地形特征线对水流作用的显著差异性，其在水系提取及流域分割等水文工程领域同样具有广泛应用[1-3]。因此，面向DDM的地形特征线精确提取方法研究具有重要的理论和现实意义。
现阶段针对地形特征线的提取研究方法有很多[4-17]，主要基于依高程数据构建的数字高程模型（digital elevation model，DEM）。根据数据组织方式不同，DEM通常分为规则格网数字高程模型（regular grid digital elevation model，GRID-DEM）和不规则三角网数字高程模型（triangulated irregular network digital depth model，TIN-DEM）。其中由于GRID-DEM具有数据结构简单、易于计算机存储与管理等优点，针对GRID-DEM地形特征线提取的研究相对深入，主要有图像处理[4]、曲面几何分析[5-8]、地表流水模拟[10-14]等方法。图像处理方法通过对固定窗口内局部最高点和最低点的位置标记实现了GRID-DEM地形特征点的快速判定，但该方法无法提供地形特征点间有效连接方式；曲面几何分析方法通过局部地形的几何形态特征分析选取地形特征点，但仍存在地形特征点连接的无序特性，且由于其严重依赖于选取方向，故地形特征线提取的完整性难以保证；地表流水模拟方法利用物理水流模型中地形特征分析由局部至整体扩展的特性，通过固定窗口格网点流水方向计算和整体格网汇水量的累积计算，实现了依汇水量累积阈值的地形特征点自动判定及依流水方向的地形特征点自适应连接，被广泛应用于提取DEM分水线与合水线的工程实践。
然而，需要指出的是：一方面，地表流水模拟方法受汇水量累积阈值的影响较大，地形特征线（山脊线和山谷线）提取结果的唯一性难以保证；另一方面，对照地形特征线定义[2]，地表流水模拟方法所提取的地形特征线缺乏山脊线与山谷线之间过渡线——坡折线的提取，难以构建完整的DEM地形结构化框架，导致地表流水模拟方法在DEM数据压缩、地形重建、地貌识别等重点关注地形保真性的应用领域相对受限。此外，由于TIN-DEM数据结构的复杂性及拓扑关系维护的低效性，故专门针对TIN-DEM地形特征线提取的研究较少。但是相比于经过内插处理的GRID-DEM，TIN-DEM直接以原始采样点作为模型数据，其内地形特征点的准确性相对较高且地形特征线与原始地形的贴合程度更加紧密，更适合作为地形特征线提取的基础。
为克服GRID-DEM地形特征线提取方法中阈值无法自动调控及地形特征线连接难以自适应调整的应用局限、解决不规则三角网数字水深模型（triangulated irregular network digital depth model，TIN-DDM）完整地形特征线提取的难点问题，本文将滚动球变换模型应用于TIN-DDM地形特征线的提取过程，分析了采样点地形曲率与临界滚动球半径识别地形特征点的相关性，建立了依临界滚动球半径的地形特征点判定准则；通过引入地形形态边界点，对非地形特征点邻接聚类构造面向地形特征线提取的剖分单元，设计了基于自然邻域的采样点地形特征属性与所属剖分单元的判定准则，建立了临界滚动球半径定量控制的地形特征线自动提取模型；通过滚动球变换的多尺度地形形态判定原理，建立了联合地形特征属性与地形形态尺度变换的地形特征点优化提取模型；通过临界滚动球半径与水深数值识别特性的差异性分析，建立了临界滚动球半径与深度信息联合提取完整且细分的山脊线、山谷线及坡折线模型。通过基于三维地形的目视解译、与经典的地表流水模拟方法对比、效率分析，验证了该方法的可行性与有效性。
1 基于滚动球变换的地形特征划分
1.1 滚动球变换概念及多尺度地形形态判定原理
滚动球变换是滚动圆变换在三维空间上的应用扩展，其物理含义可理解为半径为r的光滑球体沿地形表面滚动，在其球心轨迹处形成新增痕迹面的过程[18-20]。滚动球变换有正向和负向之分，分别表示光滑球体沿地形表面上方和下方的滚动过程。其物理表述如图1所示，黑色实线表示地形表面的纵向剖面，虚线表示滚动球球心轨迹处新增痕迹面的纵向剖面。通过滚动球痕迹面与原始地形表面的特征比对分析，可明确正向滚动变换保凸减凹及负向滚动变换保凹减凸的几何特性[19]。
图1滚动球变换的物理表述
Fig.1Physical representation of rolling ball transformation
为解决连续尺度变换条件下TIN-DDM地形形态的划分问题，前人将滚动球变换的几何特性运用至TIN-DDM地形形态判定的过程中，通过对滚动球接触点状态变化与滚动球半径间的数值关系分析，提出了以临界滚动球半径为分割边界的TIN-DDM采样点Pi地形类型判定准则[21]。即：
(1)
式中：R表示空间尺度阈值（即预先设置的滚动球半径）；r′表示正向临界滚动球半径；r″表示负向临界滚动球半径；Q（Pi）表示采样点Pi的地形类型属性，其中，1表示正向地形点，-1表示负向地形点，0表示平坦地形点，2表示正向嵌套内负向地形点，-2表示负向嵌套内正向地形点。然后通过对TIN-DDM中采样点Pi地形类型属性的邻接聚类分析，实现了任意尺度条件下TIN-DDM地形形态的连续尺度划分。
1.2 依临界滚动球半径的地形特征点判定准则
相比于通过正负向临界滚动球半径数值比较，定性描述局部地形表面起伏状况的地形形态，地形特征则侧重于地形形态范围中任意采样点特征数值（如地形曲率等）的定量比较。为探明采样点临界滚动球半径与地形特征数值之间的关联机理，将以地形曲率为参照对象，分析两类参数数值变化规律的相关性。如图2所示为TIN-DDM中任意负向地形剖面，其中：A~G为负向地形采样点；沿着采样点A至D（G至D）的方向，各采样点的地形曲率逐渐增大，在D（负向地形特征点）处达到极大值；负向地形剖面上各采样点以过D的垂线对称；各采样点的临界滚动球剖面以红色虚线表示，且球心位于采样点D的垂线上。
图2依临界滚动球半径的地形特征识别特性示意图
Fig.2Schematic diagram of terrain feature recognition characteristics according to critical rolling ball radius
由图2可知，沿着采样点A至D（G至D）方向，相比于数值逐渐增大的地形曲率，各采样点的临界滚动球半径数值逐渐减小，且地形曲率的最大采样点D为A至D（G至D）方向上临界滚动球半径的最小点。由此，临界滚动球半径与地形曲率的数值变化趋势相反，但在突变地形（地形特征点）处的数值变化规律一致。相比于地形曲率，TIN-DDM中任意采样点Pi被定量识别为地形特征点的判定准则为：
(2)
式中，表示采样点Pi的临界滚动球半径，Ωpou（Pi）表示某一剖面上采样点Pi的邻域范围，Tj表示Ωpou（Pi）范围内采样点，表示采样点Tj的临界滚动球半径，min表示取极小值。式（2）的含义可理解为：若在采样点Pi的某一剖面邻域范围Ωpou（Pi）内，采样点Pi的临界滚动球半径在某一剖面邻域范围Ωpou（Pi）内取极小值，则采样点Pi被判定为地形特征点。
2 临界滚动球半径定量控制的地形特征线自动提取模型
滚动球变换的地形特征点判定准则减少了地形特征点提取中的阈值设定环节，进而摒除算法中的人工干预、增强算法自动化水平。然而，该准则本质上是传统断面极值法在滚动球变换中的特殊应用，受制于断面极值法的构建原理，难以克服地形特征点判定对方向的依赖性问题；更为关键的是，由于分析对象为单一地形形态，相比于地表流水模拟方法中依流水方向的地形特征点自适应连接方式，其在地形特征点的判定过程中无法建立地形特征点之间的有效关联。为此，本节采用逆向工程的建模思路，在引入地形形态边界点概念的基础上[18]，通过地形形态边界点与非地形特征点的邻接聚类，构造基于地形形态边界点剖分单元。由于地形形态边界点剖分单元的邻接共边特性，地形特征线与地形形态边界点剖分单元边界等价。
2.1 依邻接聚类分析的地形特征线逆向工程构造原理
引入地形形态边界点的概念，结合滚动球变换的多尺度地形形态判定原理，相邻地形特征点之间的地形形态及与之对应的采样点临界滚动球半径数值之间呈相关性变化规律，以负向地形为例，负向地形特征点与负向地形形态边界点之间的采样点随着地形形态由负向至正向变化，其正向临界滚动球半径逐渐增大，如图3所示。
图3滚动球变换的地形弯曲识别特性示意图
Fig.3Schematic diagram of terrain bending recognition characteristics of rolling ball transformation
图3中，负向地形形态边界点C与负向地形特征点B（D）之间采样点的正向临界滚动球半径逐渐增大，故在B、C、D所构成的正向地形中，采样点C为正向临界滚动球半径极大值点，并以采样点C为边界实现了相邻负向地形（A、B、C所构负向地形与B、C、D所构正向地形）的有效分割。基于上述分析，将TIN-DDM中任意采样点Pi判定为地形特征点的等价条件为该点实现了相邻地形形态的有效分割，且采样点Pi的邻域范围Ω（Pi）边界上必然存在两个或以上的采样点，满足其在各自邻域范围内的临界滚动球半径为最大的条件。即：
(3)
式中，、分别表示采样点Pi邻域范围边界ΩB（Pi）上的任意两个采样点Pk、Pl的临界滚动球半径，max表示取极大值，ΩB表示任意采样点的邻域范围边界。式（3）含义为：在ΩB（Pk）、ΩB（Pl）邻域范围内满足临界滚动球半径取极大值的采样点Pk、Pl为采样点Pi所在地形的边界点，且ΩB（Pk）、ΩB（Pl）相交于采样点Pi，则采样点Pi实现了地形形态边界点Pk、Pl所在地形形态的有效分割。在式（3）的基础上，本文可选用逆向工程建模思路，即通过地形特征点判断与地形形态边界点搜索的问题转化，采取从地形形态边界点至地形特征点的逆向搜寻，可在有效分割地形形态的同时实现采样点邻域范围的唯一确定，避免通过剖面判定地形特征点所存在的方向依赖性问题，地形形态边界点剖分单元如图4所示。
图4中，地形形态边界点及其周边的非地形特征点共同构成地形形态边界点剖分单元。类比于图3中的地形剖面，图4中采样点A、C为相邻地形形态边界点，B为地形特征点，地形形态边界点A的剖分单元Ω（A）与C的剖分单元Ω（C）相邻，依据地形形态的连续渐变特性，Ω（A）与Ω（C）的公共边即为地形特征线，且地形特征点B位于A、C连线的剖面上。基于上述分析，通过逆向工程的建模思路，将滚动球变换的地形特征线自动提取转换至地形形态边界点的判定及其剖分单元的构造，以实现地形特征点有效判定与连接。
图4地形形态边界点剖分单元示意图
Fig.4Schematic diagram of boundary point sub division unit of terrain form
2.2 面向地形特征线网络构造的采样点特征属性判定及所属剖分单元划分
由地形形态边界点剖分单元与地形特征线的相关性分析可知：在地形形态边界点至其剖分单元边界线（即地形特征线）的任意方向上，临界滚动球半径均呈现逐渐减小的数值特性。因此，一方面，对于由n个采样点Pi构成的集合Λ={P1，P2，···，Pi，···Pn|i∈（1，n）}，按照临界滚动球半径数值由大到小的顺序对各采样点进行顺序整理，形成新的采样点集ΛNEW={P′1，P′2，···，P′j，···，P′n|j∈（1，n）}。另一方面，对任意采样点被判定为地形形态边界点的邻域范围可简化为该采样点所在的自然邻域。故，对于任意采样点P′j，若其临界滚动球半径在其自然邻域Ω0（P′j）内为最大，则将P′j判定为地形形态边界点，即：
(4)
在判定地形形态边界点的基础上，依据临界滚动球半径逐渐减小的采样点顺序，对与地形形态边界点邻接的采样点进行聚类构造剖分单元，直至当前地形形态边界点剖分单元与其他地形形态边界点剖分单元相邻接。记地形形态边界点集为Ψ={P″1，P″2，···，P″k，···，P″m|k∈（1，m）}（m为点集ΛNEW中地形形态边界点的数量），Ψ中任意地形形态边界点P″k所构剖分单元为Ω（P″k）（初始条件下Ω（P″k）中仅包含地形形态边界点P″k）。为实现每个剖分单元在空间邻域上的逐步扩充，依据点集ΛNEW中的采样点序，判断采样点P′j与剖分单元内采样点的邻接状态，若采样点P′j仅与剖分单元Ω（P″k）内的采样点邻接（两者均处于同一Delaunay三角形边上），则采样点P′j为剖分单元Ω（P″k）的非地形特征点并更新剖分单元Ω（P″k），从而可将Ω（P″k）内非地形特征点与地形形态边界点P″k进行邻接聚类。随着剖分单元范围的逐渐扩大，任意采样点必然存在与多个剖分单元内其他采样点邻接的情形。若采样点P′j和相异地形形态边界点P″k、P″l所构剖分单元Ω（P″k）与Ω（P″l）内的采样点同时邻接，则采样点P′j为剖分单元Ω（P″k）与Ω（P″l）公共边界线上的点，此时采样点P′j在空间分布上既属于剖分单元Ω（P″k）也属于剖分单元Ω（P″l）。
基于上述分析，地形形态边界点剖分单元的邻接状态分析及地形特征线提取的关键在于对任意采样点所属剖分单元的判断。记任意采样点P′j所属剖分单元的索引号构成集合Index（P′j）（初始为空），依据点集ΛNEW中的采样点序，判定采样点P′j所属剖分单元Ω（P″k）并更新索引号集合Index（P′j）。若P′j自然邻域Ω0（P′j）内所有采样点Ty同属于剖分单元Ω（P″k），即各采样点Ty索引号集合Index（Ty）的交集为{k}，则更新索引号集合Index（P′j）为{k}，并判定P′j为非地形特征点。若采样点P′j自然邻域Ω0（P′j）内各采样点Ty不属于同一剖分单元Ω（P″k），即各采样点Ty索引号集合Index（Ty）的交集为空，则采样点P′j为地形特征点，其索引号集合Index（P′j）为Ω0（P′j）内各采样点Ty索引号集合Index（Ty）的并集。则非地形特征点与地形特征点的判定条件如下：
(5)
结合滚动球变换的地形特征线构造原理，TIN-DDM中的地形特征线可由点集Ψ中任意地形形态边界点P″k所构剖分单元的边界ΩB（P″k）间接获得。即：
(6)
式中，ΓBall表示地形特征线。由于滚动球变换有正负向之分，依正向临界滚动球半径提取的地形特征线为负向地形特征线（记为ΓNeq-Ball）；反之，依负向临界滚动球半径提取的地形特征线为正向地形特征线（记为ΓPos-Ball）。由此，地形特征线提取的具体流程可由式（4）~（6）构成，此流程主要涉及临界滚动球半径的大小比较与剖分单元索引号集合Index（Ty）的交集、并集运算，完全可实现地形特征线的自动化提取。
3 TIN-DDM地形特征线多尺度表达及补充
3.1 基于多尺度地形形态判定原理的地形特征线优化提取模型
受地形形态边界点剖分单元邻接特性的影响，基于滚动球变换提取的地形特征线在几何形态上表现为网络状结构，这与地表流水模拟方法（传统地形特征线提取方式）差异较大。上述现象产生的原因在于：基于临界滚动球半径数值变化趋势分析的地形特征线提取方式未充分考虑采样点特征属性的强度数值分析与评价，所提取的地形特征点无法满足不同应用领域中对于地形特征线的多尺度表达需求。
由于地形表达的多尺度特性，小尺度下的局部地形形态与大尺度下的整体地形形态呈现明显的层次嵌套关系，因此局部地形与整体地形形态的判定及对应地形特征点的识别结论关联性不强。TIN-DDM的地形特征线分布如图5所示，依据地形形态的连续尺度划分[21]，在设置空间尺度阈值为R的条件下，对原始TIN-DDM进行正向、负向地形类型划分。为便于讨论，图5中显示的地形特征点及其连线特指负向地形。
图5TIN-DDM地形特征线分布
Fig.5Distribution of TIN-DDM terrain feature lines
图5中，剖分单元Ω（A）、Ω（B）及Ω（C）实现了基于原始TIN-DDM地形形态的有效分割，采样点D、E由于在其局部剖面地形形态范围内对应的临界滚动球半径数值为极小，表现为处于各剖分单元边界上的负向地形特征点。但在给定空间尺度阈值R的条件下，由于局部地形的邻接及与整体地形的嵌套关系，对剖分单元Ω（A）、Ω（B）及Ω（C）进行了重新聚类与分割，导致原始TIN-DDM中局部地形范围内的负向地形特征点D、E在不同尺度条件下的地形范围中分属于不同性质的地形（D、E分别属于正向和负向地形）。依据地形特征多尺度表达特性分析，当前小尺度下的地形特征点特征属性强度及所在地形性质与较大尺度变换下的无关[21]。因此，相比于采样点E，采样点D的负向地形特征属性随着空间尺度的变换及地形嵌套关系的化简发生了改变，其负向地形特征属性相对较弱。
考虑到嵌套地形的本质为正向、负向及平坦地形的邻接聚类，且其性质的判定与其内各类地形特征属性的强弱及空间尺度相关。因此，TIN-DDM中的负向嵌套地形点的负向地形特征属性较强，其在一定空间尺度下等价于负向地形内采样点；反之，亦然。基于上述分析，在空间尺度阈值为R条件下地形特征线为：
(7)
式中：Φpositive（Φnegative）表示在空间尺度阈值R下由正向（负向）地形点与正向（负向）地形嵌套点邻接聚类构建的正向（负向）地形范围；ΓPos-Ball-R表示位于正向地形范围Φpositive内的正向地形特征线ΓPos-Ball；ΓNeg-Ball-R表示位于负向地形范围Φnegative内的负向地形特征线ΓNeg-Ball。
3.2 考虑采样点水深数值的地形特征线补充与细分模型
滚动球变换的地形特征识别关键在于基于地形形态变化的临界滚动球半径数值分析，这与地表流水模拟方法（传统的基于地形起伏变化的水深数值分析方式）存在一定差异，尤其在水文分析等应用领域，缺乏纵向深度数值分析的地形特征识别难以保证流域分析的准确性。以下以正向地形为例，分析两种分析方式在识别正向地形特征点方面的差异性，三种典型正向地形形态走势剖面如图6所示。
利用滚动球变换的地形特征识别特性可实现图6（a）中山脊点和图6（b）、（c）中坡折点的定量识别。而图6（b）中的山脊点由于其地形曲率（临界滚动球半径）为其局部区域范围内的最小值（最大值），该识别特性难以实现此类地形特征点的有效识别。参照断面极值法提取地形特征的论述[5]，基于水深数值分析的地形特征识别特性可实现图6（b）中特殊山脊点的有效判定，但存在无法提取图6（b）、（c）中坡折点的应用局限。综上所述，临界滚动球半径数值分析方式难以有效识别高地上的山脊点，且未对已识别的山脊点、坡折点进行细致分类；而水深数值分析方式尽管可实现山脊点的有效判定，但无法进行坡折点的识别。因此，若在临界滚动球半径数值分析方式基础上结合水深数值分析方式，则可有效识别和细分所有的正向地形特征点。
图6三种典型正向地形形态走势剖面图
Fig.6Trend profile of three typical positive terrain forms
由于依深度信息的地形特征点判定条件与依临界滚动球半径的地形特征点判定条件（式（2））类似，即：对于TIN-DDM中任意采样点Pi，采样点Pi的水深在其剖面邻域范围Ω（Pi）内取极小值（极大值），则Pi被判定为山谷线（山脊线）上的地形特征点。可通过借鉴滚动球变换的地形特征线提取思路，并将临界滚动球半径替换为水深数值，实现TIN-DDM中的山脊线（记为ΓRidge）、山谷线（记为ΓValley）的完整提取和顾及空间尺度阈值R的山脊线（记为ΓRidge-R）、山谷线（记为ΓValley-R）的提取。由此，可获得TIN-DDM中完整细分且顾及空间尺度阈值R的地形特征线（记为ΓFeature-R），即：
(8)
式中：ΓPos-Slope-R（ΓNeg-Slope-R）表示在空间尺度阈值R条件下的正向（负向）地形坡折线。式（8）表示：在空间尺度阈值R的条件下，地形特征线ΓFeature-R由山脊线ΓRidge-R、山谷线ΓValley-R、正向地形坡折线ΓPos-Slope-R和负向地形坡折线ΓNeg-Slope-R构成。其中：正向地形坡折线ΓPos-Slope-R由依负向临界滚动球半径提取的正向地形特征线ΓPos-Ball-R中除山脊线ΓRidge-R外的地形特征线构成；负向地形坡折线ΓNeg-Slope-R由依正向临界滚动球半径提取的负向地形特征线ΓNeg-Ball-R中除山谷线ΓValley-R外的地形特征线构成。
4 试验结果与分析
为验证本文所提方法的可行性和有效性，本节采用MATLAB实现本文方法。试验所采用的数据为深度范围在10~90 m的多波束水深数据，共包含12 774个水深点，水深点平均间隔为18 m。试验环境为Inter（R）core（TM）i5处理器，主频为2.5 GHz，内存为8 GB。
4.1 目视解译评价
为验证本文方法在地形特征线提取方面可实现特征属性强度的有效划分，且能表征地形形态变化与地形起伏变化，采用目视解译的方式对试验结果进行评价。利用本文方法从试验数据中提取完整的正向地形特征线网络，以200 m、500 m、800 m和1 100 m作为空间尺度阈值R在该网络下提取正向地形特征线，结果如图7所示，其中：红色线条表示山脊线；黄色线条表示正向地形坡折线；A区域为海底沟壑区域；B区域为平缓山顶区域。图8、9分别为A、B区域在不同R下的正向地形特征线提取图像。
图7不同空间尺度阈值下的正向地形特征线提取图像
Fig.7Positive terrain feature line selection images under different spatial scale thresholds
图8不同空间尺度阈值下A区域的正向地形特征线提取图像
Fig.8Positive terrain feature line selection image of region A under different spatial scale thresholds
试验结果表明：①在图7~9中，随着空间尺度阈值的增大，本文方法提取的正向地形特征线数量逐渐增多，地形特征线从复杂地形向平坦地形不断扩张。即随着空间尺度阈值的增大，区分地形特征线特征属性强度的标准逐渐下降。因此，基于多尺度地形形态判定原理的地形特征线提取可实现地形特征属性强度的有效划分。②本文方法在沟壑分布较密的区域A中，提取的红色地形特征线和黄色地形特征线分明，并分别与沟壑中的山脊线、明显的坡折处相重合；在山脊不明显的区域B中，提取的红色地形特征线和黄色地形特征线区分明显，并分别与平缓山顶处的山脊线以及平缓山顶的坡折线相重合。因此，基于临界滚动球半径和水深数值所提取的山脊线、正向地形坡折线真实、可靠且细分，能明显表征正向地形内地形起伏变化与地形形态变化。
图9不同空间尺度阈值下B区域的正向地形特征线提取图像
Fig.9Positive terrain feature line selection image of region B under different spatial scale thresholds
4.2 对比分析试验
4.2.1 地形特征线提取
为验证本文所提方法（以下简称方法Ⅰ）在水文分析领域方面的优势，在利用ArcGIS软件实现地表流水模拟方法[12]（以下简称方法Ⅱ）的基础上，通过对方法Ⅰ和方法Ⅱ所提取山脊线的差异性分析，实现两类方法在地形特征线提取完整性和连续性方面的量化评估。为最大程度保证方法Ⅱ提取地形特征线的完整和连续性，参照前人关于流水累积阈值的设置方法[11]，试验中将流水累积阈值设置为90。同时，考虑到方法Ⅱ的构建原理与空间尺度无关，故方法Ⅱ提取的地形特征线无多尺度变换的概念，可认为是规则格网数字水深模型（regular grid digital depth model，GRID-DDM）的初始地形特征线。为保证试验的有效性，方法Ⅰ同样不设定空间尺度阈值，以保证所提取TIN-DDM地形特征线网络的完整性。图10为两种方法提取的山脊线与5m等深线及各自DDM渲染图叠加效果图，其中，黄色虚线框区域表示海底沟壑构成条带状地形上的山脊区域，绿色虚线框表示海底沟壑构成的面状区域。图11为图10中绿色虚线框区域内两种方法提取结果的放大图，其中的黄色虚线框区域为两种方法提取山脊线的部分差异位置。
试验结果表明：相比于图10（a），图10（b）中黄色虚线框区域内的山脊线断裂状态明显；在图11中黄色虚线框标识的复杂区域（依据图例与等深线可判断该区域地形为处于两侧明显沟壑山脊间的马鞍状地形），相比于图11（a），图11（b）中出现黄色虚线框内马鞍状地形处山脊线无法提取的情况。造成上述现象的原因在于：由于方法Ⅱ中汇水累积量与高程密切相关，高程值较低的山脊线上的汇水量较大，故易被误判为非山脊点，因此山脊线断裂状态明显；而由于方法Ⅰ采用与阈值无关的特征数值趋势变化分析方式识别山脊点，通过构建地形形态边界点剖分单元，提取其连续边界线作为山脊线，故提取的山脊线完整且连续。
图10方法Ⅰ与方法Ⅱ的提取结果
Fig.10Extraction results of method Ⅰ and method Ⅱ
图11在海底沟壑构成的复杂区域内方法Ⅰ与方法Ⅱ的提取差异对比图
Fig.11Comparison of extraction differences between method Ⅰ and method Ⅱ in the complex area composed of seabed gullies
4.2.2 地形重构
为验证方法Ⅰ提取的地形特征线在地形重构方面的优势，分别利用无空间尺度阈值约束的方法Ⅰ与方法Ⅱ，提取DDM中的地形特征线，并进行地形重构的精度分析。其中基于方法Ⅰ提取的正向地形特征点数量为6 193个（包含3 421个山脊点和2 772个正向坡折点）；负向地形特征点数量为2 921个（包含1 501个山谷点和1 420个负向坡折点）。图12所示为地形重构前后的等深线（等深距为5 m）的叠合图像。由于两种方法提取地形特征线所基于的模型不同，故两种方法重构地形的参照DDM分别为原始TIN-DDM以及经格网插值后的GRID-DDM。
试验结果表明：相比于方法Ⅱ，在地形复杂区域A、B内，基于方法Ⅰ的重构等深线与原始等深线贴合程度更高，即经方法Ⅰ提取的地形特征线在地形重构的准确性方面优势明显。更进一步地，以复杂区域B为例，分析对比两种方法提取的地形特征线形态差异。如图13~15所示，分别为区域B内基于两种方法进行地形重构前后的等深线叠合图、地形特征线平面图和方法Ⅰ提取地形特征线的三维图。图13~15中的黑色虚线椭圆区域为基于两种方法进行地形重构前后的等深线差异较大区域。
图12重构前后模型的等深线叠合
Fig.12Contour superposition map of the model before and after reconstruction
试验结果表明：相比于方法Ⅱ，方法Ⅰ在黑色虚线椭圆区域内提取的地形特征线相对完整，体现在方法Ⅰ中黄色线条（坡折线）对于红色线条（山脊线）的形态补充。因此，方法Ⅰ提取的地形特征线在地形重构的准确性方面优势在于其对坡折线的有效识别。为定量评价方法Ⅰ在地形重构方面的准确性优势，对基于两种方法进行地形重构后的DDM与原始DDM进行相关差异性数值的统计分析。具体包括：重构DDM与原始DDM在垂直方向的最小差值（Hmin）、最大差值（Hmax）、重构DDM及原始DDM 的整体体积差与投影面积比值（即平均深度差（Haverage））、重构DDM高于原始DDM的体积差（V1）、重构DDM低于原始DDM的体积差（V2），统计结果如表1所示。
图13不同方法重构前后等深线叠合图
图13Contour superposition map before and after reconstruction by different methods
图14不同方法的地形特征线平面图
Fig.14Terrain feature line plan of different methods
图15方法Ⅰ地形特征线三维图
Fig.15Method Ⅰ 3D map of terrain feature line
表1方法Ⅰ与方法Ⅱ进行地形重构后的DDM与原始DDM差异性数值统计
Tab.1 Statistics of difference between DDM and original DDM after terrain reconstruction by method Ⅰ and method Ⅱ
试验结果表明：相比于方法Ⅱ，方法Ⅰ除在Hmin与方法Ⅱ保持相同之外，其他各项差异性数值均表现较小，定量验证了方法Ⅰ在地形重构方面的准确性优势。尤其是Hmax与Haverage的相对数值比较更是明显反映出，相比于局部地形重构精度，坡折线对于整体地形重构精度的影响更大。此外，作为代表负向地形特征线形态差异的统计量V1，相对于正向地形特征线形态差异的统计量V2的数值差异，不仅直接反映出方法Ⅰ中重构地形与原始地形的贴合程度优势，也间接验证了方法Ⅰ中提取的正负向地形特征点的数值差异。
4.2.3 效率分析
在DDM数据压缩、地形重建、地貌识别、水文分析等相关应用领域中，地形特征线的提取及尺度变换效率十分重要。对于采样点数量为n的点集ΛNEW，计算所有采样点的临界滚动球半径时间复杂度为O（n2）[21-22]。在假设P′j自然邻域Ω0（P′j）内采样点Ty数量为a的前提下，本文方法中地形特征线网络提取的关键在于P′j所属剖分单元的判定及基于剖分单元边界线的地形特征线提取。对于任意P′j，其所属剖分单元的判定需遍历Ω0（P′j）内所有采样点Ty，并进行采样点Ty剖分单元索引号集合Index（Ty）的交集和并集运算，共需进行2a次运算；而基于剖分单元边界线的地形特征线提取同样需要遍历Ω0（P′j）内所有采样点Ty，即进行a次运算。在此基础上，对于在给定空间尺度阈值R下与P′j地形特征属性相一致的地形范围Φ特征线的提取，需遍历Ω0（P′j）内所有采样点Ty，共需进行a次运算。结合坡折线提取的a次运算，并考虑到依临界滚动球半径和水深数值提取正负向地形特征线的4次点集ΛNEW遍历，故本文方法的时间复杂度为O（20a×n）≈O（n）。
由此结合数据预处理过程，本文方法提取地形特征线的总体时间复杂度为O（n2+n）≈O（n2），且方法Ⅱ的时间复杂度同为O（n2）[14]。需要指出的是：尽管两种方法的时间复杂度相同，但在地形特征线的多次提取或多尺度表达应用中，方法Ⅰ的时间复杂度可快速降至O（n），试验结果如表2所示，其中，为获取TIN-DDM中采样点的临界滚动球半径数值，需要对TIN-DDM数据进行预先处理，预先处理共耗时112.34 s。
表2试验耗时统计
Tab.2 Statistics of test time
试验结果表明：①相比于方法Ⅱ，尽管方法Ⅰ需进行前期数据的预处理，且耗时占比较大，但结合预先处理时间后整体效率相对较高。前期预处理中计算所得数据链与临界滚动球半径可用于地形形态快速划分[21]与缓冲面快速构建[22]等领域。因此，在多种相关领域使用中，方法Ⅰ前期预处理所占比重将明显下降。②进行一次数据预处理后，方法I在不同空间尺度阈值R下地形特征线提取效率均衡，且相比于算法Ⅱ的效率优势明显。
5 结论
针对传统地形特征线提取方法存在的阈值调控及地形特征线连接问题，结合不同应用领域的不同尺度的应用需求，本文提出了基于滚动球变换的TIN-DDM地形特征线自动提取方法。
试验分析表明：相比于地表流水模拟方法，本文方法可快速实现山脊线、山谷线、坡折线的多尺度提取及有效分类，并且所提取的地形特征线在水文分析与地形重构方面具有明显优势。需要指出的是，本文在获取临界滚动球半径的前期数据预处理过程中，未对采样点建立空间格网索引和并行运算，致使算法前期数据预处理耗时较大。在多尺度提取中选用滚动球半径为参照尺度，未与传统尺度建立对应关系。由此，下一步将在TIN-DDM数据分块处理、算法并行运算、尺度对应等方面加以研究。
