基于等高線的地形特征提取算法研究

任振娜

（武警指揮學院，天津 300250）

[關鍵字] 等高線數據；Delaunay三角網；地形特征；提取

1 引言

地形特征是指對于描述地形形態(tài)有著特別意義的地形表面上的點、線、面，它構成了地形變化起伏的骨架。其中，地形特征點包括：山峰點、谷底點、鞍部點等，地形特征線包括：山脊線、山谷線等。地形特征是地形地表在空間分布、延展的具有骨架化控制作用的重要結構化信息，在地圖制圖、地貌形態(tài)識別、水文分析等領域具有重要的支持作用。由于歷史的原因，現有的地形數據更多的是等高線，DEM建立的數據來源也主要是等高線，這樣如何基于等高線建立DEM模型，如何基于建立的模型提取地形結構特征便成為一個感興趣的問題。

總體思路：首先利用等高線數據構建不規(guī)則三角網DEM模型，然后通過對三角網進行分析，提取地形結構特征[1]。而本文是在基于等高線一次性構建Delaunay三角網模型的基礎上，進一步分析了模型的特點，研究了山脊線和山谷線的提取算法。

2 基于等高線一次性構建Delaunay三角網模型算法描述

2.1 算法中用到的概念

定義1：約束邊，本文中稱由等高線上的相鄰兩點確定的線段為約束邊。

定義2：約束邊的起始點和終止點，設約束邊是由等高線上相鄰兩個點M，N逆時針相連而成，則稱M點是約束邊的起始點，N點是約束邊的終止點。

定義3：基邊，以某一邊為基準邊，在點集中選取第三點與其構成一個三角形，該邊就稱為基邊。

定義4：左右三角形，以約束邊的起始點→終止點的方向為基準方向，約束邊與位于基準方向右側的點組成的三角形稱為右三角形，該點稱為該條約束邊的右側第三點，約束邊與位于基準方向左側的點組成的三角形稱為左三角形，該點稱為該條約束邊的左側第三點。

定義5：平三角形，是指三角形的三個頂點都在同一條等高線上。

2.2 基于等高線構建Delaunay三角網算法[2]

算法的具體過程描述如下：

步驟1：將每一段等高線上的點數據以其X，Y坐標，進行逆時針排列，相鄰兩點構成一條約束邊。記錄點數據到點數據鏈表,記錄約束邊數據到邊鏈表中，并將每一條等高線數據作為一個子元素記錄到等高線數據鏈表。

步驟2：根據點數據鏈表創(chuàng)建點數據格網索引。

步驟3：置標志位于邊鏈表的表尾。

步驟4：在等高線鏈表中，依次將每一個子元素中的每一條約束邊作為基邊，應用夾角最大準則生成左右三角形。更新點數據信息，記錄新生成的邊到邊鏈表，記錄新生成的三角形到三角形鏈表，更新基邊信息。

步驟5：從標志位的下一條邊開始(即第一條非約束線段)，依次取出邊作為基邊，向左應用夾角最大準則生成新的三角網，更新點數據信息，記錄新生成的邊到邊鏈表，記錄三角形信息，更新基邊信息。

步驟6：從標志位的下一條邊開始，依次對同時存在左右三角形的非約束邊，判斷該邊與其左右側第三點組成的邊是否相交，如果相交，則應用簡化的LOP優(yōu)化算法對三角網進行優(yōu)化處理。

3 地形特征提取算法

3.1 特征點提取算法

山脊線和山谷線上的點在等高線上的特征表現為等高線局部曲率最大點，也就是等高線彎曲變化的特征點。曲線特征點的提取算法有多種。其中Split法是一種較常用的曲線特征點提取算法。該法從原理上講屬于整體算法，它所提供的曲線特征點能夠保證曲線變形在規(guī)定的限差之內。

Split法的基本思想是[3,4]：先用曲線的最左邊和最右邊的兩個點作為起始點(對于閉合曲線)將閉合曲線分為兩部分，對于非閉合曲線選擇其兩個端點作為起始點。起始點確定后，順序計算曲線上位于兩個起始點之間的每一個點距兩個起始點連線的垂距，并找出最大垂距點。若該點處等高線張角小于給定的閾值，則該點為特征點。它將原曲線分為兩部分，對每一部分確定新的起始點，即用該點分別與原兩個起始點構成兩對新的起始點，用相同的方法對這兩段曲線找出各自的特征點。Split法的原理可用圖1表示。

圖1 Split法的原理

圖2 識別山脊點山谷點原理圖

3.2 山脊點和山谷點的識別算法

應用Split法，即可找出特征點。找出特征點后，還要對其進行分析判斷，找出哪些特征點是山脊點，哪些特征點是山谷點。

最直觀的判斷山脊點山谷點的方法如圖2所示，即計算特征點C處等高線張角范圍內某點D的高程，并與C點高程比較，如果D點的高程大于C點高程，則C點為山脊點，反之，C點為山谷點。此種方法看似簡單，但計算D點高程算法復雜，計算量大，影響判斷效率。本文建立的Delaunay三角網模型是以等高線線段為約束邊構建的，因此，模型中每一個平三角形三個點均落于一條等高線上；每一個非平三角形三個點分別位于兩條等高線上。因此，在識別山脊點和山谷點時，主要步驟如下：

步驟1：從等高線鏈表中，依次選出第n條線段ln和第n+1條線段ln+1[5,6]，以及兩條線段的三個端點dn、dn+1、dn+2。

步驟2：計算以ln為起始方向線，順時針到ln+1的角度值。

步驟3：當角度值>195°（195°是閾值）時：根據ln找到其右三角形的第三個點，當第三個點的高程值小于dn+1的高程值，則點dn+1是山谷點，更新點信息；當第三個點的高程值大于dn+1的高程值，則點dn+1是山脊點，更新點信息；當第三個點的高程值等于dn+1的高程值，則轉入步驟5。

步驟4：當角度值<165°（165°是閾值）時：根據ln找到其右三角形的第三個點，當第三個點的高程值小于dn+1的高程值，則點dn+1是山脊點，更新點信息；當第三個點的高程值大于dn+1的高程值，則點dn+1是山谷點，更新點信息；當第三個點的高程值等于dn+1的高程值，則轉入步驟6。

步驟5：根據ln找到其左三角形的第三個點，如果第三個點的高程值大于dn+1的高程值，則點dn+1是山谷點，更新點信息；如果第三個點的高程值小于dn+1的高程值，則點dn+1是山脊點，更新點信息；如果第三個點的高程值等于dn+1的高程值，則轉入步驟7。

步驟6：根據ln找到其左三角形的第三個點，如果第三個點的高程值大于dn+1的高程值，則點dn+1是山脊點，更新點信息；如果第三個點的高程值小于dn+1的高程值，則點dn+1是山谷點，更新點信息；如果第三個點的高程值等于dn+1的高程值，則轉入步驟8。

步驟7：從ln開始，圍繞dn+1依次找其右三角形，直至右三角形的某一點高程值不等于dn+1的高程值為止。當該點高程值小于dn+1的高程值，則點dn+1是山谷點，更新點信息；當該點高程值大于dn+1的高程值，則點dn+1是山脊點，更新點信息。

步驟8：從ln開始，圍繞dn+1依次找其右三角形，直至右三角形的某一點高程值不等于dn+1的高程值為止。當該點高程值小于dn+1的高程值，則點dn+1是山脊點，更新點信息；當該點高程值大于dn+1的高程值，則點dn+1是山谷點，更新點信息。

3.3 山脊線和山谷線的生成

3.3.1 連接山谷線

連接山谷線是按照由高向低逐條等高線來搜索的，其搜索過程如下：

步驟1：先從高程最高的等高線起，找出沒有向下連接的山谷點dn，以及與其相連接的兩個點dn-1、dn+1和兩條線段ln、ln+1。

步驟2：以圍繞dn，搜索邊ls從ln開始，向右搜索三角形的第三點d3，當d3不是山谷點時，ls=d3dn，繼續(xù)向右搜索三角形的第三點，直至d3是山谷點或者ls=ln+1時為止。

步驟3：當存在山谷點d3時,將dn與d3相連接。

步驟4：重復1、2、3步，直至所有的山谷點都向下搜索不到新的山谷點為止。

3.3.2 連接山脊線

連接山脊線則是由低向高逐條等高線來搜索的，其搜索過程如下：

步驟1：先從高程最低的等高線起，找出沒有向上連接的山脊點dn，以及與其相連接的兩個點dn-1、dn+1和兩條線段ln、ln+1。

步驟2：以圍繞dn，搜索邊ls從ln開始，向左搜索三角形的第三點d3，當d3不是山脊點時，ls=d3dn，繼續(xù)向左搜索三角形的第三點，直至d3是山脊點或者ls=ln+1時為止。

步驟3：當存在山脊點d3時,將dn與d3相連接。

步驟4：重復1、2、3步，直至所有的山脊點都向上搜索不到新的山脊點為止。

3.3.3 說明

連接山谷線時，在某點向下找不到可以連接的山谷點的原因可能是此山谷線到達了封閉的盆地底部或谷底線流入湖泊或海洋，另一種可能是等高線在這一部分分布太稀，數字化點太少而不能足夠詳細地描述地貌形態(tài)；連接山脊線時，在某點向上找不到可以連接的山脊點的原因可能是此山脊線到達了山頂；在某點搜索到的山谷線(山脊線)的下一點是別的山谷線(山脊線)上的點的原因是兩條山谷線(山脊線)匯合在此點。

4 實驗實例

為了檢驗本文算法的正確性和有效性，筆者在VisualC++編程環(huán)境中對上述算法進行了實踐。本實驗采用的數據是某一地區(qū)的147條數字化地形等高線。圖3為原始等高線數據，圖4為基于等高線數據構建的局部Delaunay三角網，圖5為提取的山谷線，圖6為提取的山脊線。從圖4-6所示的實驗結果可以看出，用本文所提出的算法提取出的山脊線和山谷線是較為準確的，與實際地形也是相符合的。

圖3 等高線數據

圖4 Delaunay三角網模型

圖5 提取的山谷線

圖6 提取的山脊線