基于約束連接方向的最速上升法提取山脊線

鄧曉龍， 江文萍， 孫二洋， 張 萍， 于 池

(武漢大學 資源與環境科學學院， 武漢 430079)

?

鄧曉龍， 江文萍*， 孫二洋， 張 萍， 于 池

(武漢大學 資源與環境科學學院， 武漢 430079)

對基于格網DEM提取山脊線的方法進行了研究，針對基于格網DEM的流水物理模擬法和最速上升法提取山脊線所存在的問題，提出了基于約束連接方向最速上升法.該方法將山頂點和鞍部點作為控制點、以約束山脊線的連接方向來達到改善山脊線提取效果的目的.減少了山脊線斷裂、交叉和回轉等錯誤現象，同時提取的山脊線主次分明層次性良好，與實際地形相符合.

格網DEM; 山脊線; 約束連接方向; 流水物理模擬法; 最速上升法

地性特征線是地貌形態的骨架線，是數字形式描述和處理地貌時的控制線，它主要包括山脊線、山谷線，山脊線的自動提取在水文分析、高質量DEM生成、制圖綜合、土壤侵蝕分析、遙感數字圖像處理以及解析地貌等方面均有著十分重要的作用[1-2].由于格網DEM獲取方法的廣泛性和數據精度不斷提高的特性，基于格網DEM自動提取地性特征線的方法被普遍使用和研究[1].基于規則格網DEM數據的山脊線提取方法從算法設計原理上可以分為以下兩種：基于地形表面幾何形態分析方法[3]和基于地形表面流水物理模擬方法[4].

基于地形表面形態分析方法主要是考慮點與其領域內其他點高程之間的關系來確定該點的特征類型[3]，該方法簡單易于實現，但是容易把不是特征點的點誤判為特征點[5].許多學者對該算法進行改進，但是都不能徹底解決該問題[6-7].此外，最速上升法提取山脊線的算法復雜度低且易于實現，在單獨山頂點情況下該方法提取的山脊線與實際地形符合良好[8]，而且該方法有效地改善了應用比較廣泛的離散拉普拉斯算法(D.L.T)[9]在提取山脊線時發生的粘連的情況.但是最速上升法在多山頂點和鞍部點且交替連續出現的情況該方法考慮不足，提取山脊線在山頂點處可能發生斷裂，使得主要的山脊線發生斷裂不連續，而且該方法提取山脊線時沒有考慮與山谷線相交的可能性.

基于規則格網DEM的流水物理模擬算法自1984年提出就被廣泛的研究和應用[10-11]，該方法首先要對格網DEM進行填平洼地和確定平地區域流水方向的數據預處理.有許多改善洼地處理效果的算法，但是如何確定平地區域的流水方向一直是該方法的難題[12-13]，如果格網數據量比較大且洼地、平地比較多，則該方法提取地性特征線的算法復雜度與時間復雜度都比較高[14-15].流水物理模擬法提取山脊線有兩種辦法[10]：一是將DEM高程矩陣整體反轉，進行洼地和平地的預處理，然后提取出的特征線就是山脊線(分水線)，但是該方法同樣面臨著洼地處理與平地流向難確定的問題；二是從分水嶺中提取山脊線，但是如何連接分水嶺離散點成山脊線是主要難題，因為這些分水嶺的點只有流出的流向沒有流入的流向，所以不能像流水物理模擬法提取山谷線時按照流向連接特征線[16]. 針對最速上升法和流水物理模擬法提取山脊線存在的問題，本文對基于控制點的山脊線提取方法進行了研究.

1 基于約束連接方向的最速上升法

1.1 提取主要山脊線

1.1.1 建立坐標系 匯水流量矩陣中匯水流量值等于零的格網單元為分水嶺離散點，本模型采用了簡單且容易實現的方法從這些分水嶺離散點中提取山頂點和鞍部點[17]，依次以離散點為中心建立3×3的窗口，如果窗口中心點的高程比其周圍的八鄰域的高程都大則該中心點為山頂點，圖1中Pc2和Pc5所示的點為山頂點；如果中心點兩個相互正交的方向上，一個方向為凹陷，另一個為凸起，則該中心點為鞍部點，圖1中Pc3所代表的是鞍部點.接下來通過離散的分水嶺的點連接山頂點，連接山頂點時首先以該山頂點為原點建立二維X、Y坐標系如圖2所示，該坐標系分為第一象限、第二象限、第三象限、第四象限、X正半軸、Y正半軸、X負半軸和Y負半軸8個區域，坐標系的8個區域分別對應八個山頂點隊列.把其余山頂點根據其在該坐標系中的位置依次添加到相應區域的山頂點隊列中，并按照山頂點到原點的距離從小到大的順序排序每個隊列中的山頂點.

圖1 山頂點和鞍部點Fig.1 Hill and saddle points

圖2 山頂點的坐標系Fig.2 Coordinate system of hill point

本模型把8個區域隊列中要連接的山頂點定義為“終止山頂點”，原點處的山頂點為“起始山頂點”，連接過程中山脊線最后一個山脊點為“連接山脊點”.依次從8個區域的山頂點隊列中取出一個山頂點，先判斷取出的山頂點與“起始山頂點”是否連接，如果連接過則取出山頂點隊列中下一個進行判斷，否則創建存儲連接山脊線的點數組，并把起始山頂點加入到山脊線點數組中.取出山脊線點數組中最后一個山脊點為“連接山脊點”，如果該“連接山脊點”為“終止山頂點”則結束了一條山脊線的連接，否則以該“連接山脊點”為原點建立二維的X、Y坐標系，并劃分為16個區域如圖3所示.然后計算“終止山頂點”在連接山脊點坐標系中的位置.

圖3 連接山脊點的坐標系Fig.3 Coordinate system of connection point

圖4 連接點的八鄰域Fig.4 Eight neighborhoods

1.1.2 建立連接優先級 “連接山脊點”的八鄰域如圖4所示，“連接山脊點”八鄰域內一般都會有多個可以連接的分水嶺離散點.本文定義了這些八鄰域內的分水嶺離散點添加到山脊線的優先級為Priority，并用Px表示‘X’點的優先級.在連接山脊線時本模型實時計算山脊線最后一個連接點與“終止山頂點”的位置關系，根據 “終止山頂點”在連接山脊點坐標系中的位置來確定連接的山脊線的大致走向，本模型用山脊線最后一個連接點指向“終止山頂點”的方向向量表示山脊線的大致走向.接下來分別計算連接點指向其八領域內的分水嶺離散點的方向向量，分水嶺離散點的方向向量與山脊線大致走向的夾角和優先級成反比，夾角越小則優先級越大，夾角越大則優先級越小，如果夾角相等則優先級相等.根據優先級從大到小的順序判斷下一個可能的山脊線連接點，當兩個連接點的優先級相等則高程值大優先加入到連接山脊線的點數組中.在連接任意一條山脊線的過程中不考慮已經加入山脊線的分水嶺離散點，這樣可以有效地避免山脊線發生回轉、與自身交叉等錯誤現象.“終止山頂點”在坐標系中16個區域所對應的下一個連接點可能位置和其優先級匯總表如表1所示.

表1 山脊線連接點的優先級匯總表

1.1.3 連接鞍部點 經過上述步驟之后大部分鞍部點在山頂點連接過程中都被連接起來，如果有鞍部點沒有連接過則把該鞍部點為原點建立如圖2所示的坐標系，因為鞍部點的特點：在某一方向上在鞍部點的兩邊往往存在著兩個山頂點.本模型分別計算出鞍部點坐標系中16個區域內距離原點最近的山頂點，并選擇兩個不在同一個區域且離鞍部點最近的山頂點，以該鞍部點為起始點分別向兩個山頂點方向連接，連接方法與上述的連接山頂點方法類似，通過建立的分水嶺離散點的優先級隊列來判斷連接的下一個山脊點，如果連接過程中遇到山頂點或則與其他山脊線相交則結束連接.連接結束后一般會有兩條山脊線，這兩條山脊線有著相同的起點即沒有連接過的鞍部點，把這兩條山脊線合并成一條山脊線.

1.2 次要山脊線的提取

在連接次要的山脊線之前需要對已經連接的主要山脊線和所有的分水嶺離散點進行如下處理：連接成主要山脊線的分水嶺離散點設置為連接過標志，連接主要山脊線過程中遍歷的離散點設置為檢測過標志.在剩余的分水嶺離散點中查找沒有連接過和檢測過且高程值最小的開始連接次要山脊線，連接每一條山脊線時根據連接山脊線點數組中最后兩個連接點來約束連接方向.圖5中表示了山脊線最后兩個連接點在八鄰域內的具體情況，其中BP表示前一個連接點，EP表示最后一個連接點，深色區域表示該處的點已經被排除，連接的時候不予考慮，因為這些連接點可能導致山脊線發生回轉和交叉.在連接點的八鄰域內查看淺色區域位置上的連接點，找出其中高程值最大且要大于最后連接點的離散點，如果該點沒有連接過和沒有檢測過則直接加到連接山脊線點數組中繼續進行連接；如果該點連接過則加入到連接山脊線點數組中并且結束這條山脊線的連接開始下一條山脊線的連接；如果該點檢測過則要判斷該點八鄰域內是否存在連接過的山脊點，如果存在則把該點和連接過的山脊點依次加入到連接點數組中并結束連接，如果不存在則結束山脊線的連接.

圖5 次要山脊線的最后兩個點約束連接方向Fig.5 Automatic constraint connecting direction

2 實驗與分析

2.1 實驗過程與數據

本模型的實驗環境是Windows系統和VS2010，算法采用C++編寫，在本模型的算法中定義了若干變量，其中DemMatrix為原始DEM矩陣，DemMatrixPre是經過洼地平地處理和計算了匯水流量的DEM矩陣，HillSaddleArray是存在于DemMatrixPre中的山頂點和鞍部點，MajorRidgeArray是在山頂點和鞍部點約束下連接的主要山脊線，MiniorRidgeArray是在山脊線最后兩個點約束下連接的次要山脊線.本模型算法主要分為原始DEM數據的預處理、山頂點和鞍部點、連接主要山脊線和連接次要山脊線，偽代碼如下所示：

① DemMatrixPre=PretreatDemdata(DemMatrix);

② HillSaddleArray=LocateHillSaddle(DemMatrixPre);

③ MajorRidgeArray=ConnectMajorridge(HillSaddleArray);

④ MiniorRidgeArray=ConnectMiniorridge();

本次實驗數據是從SRTM網站上下載的90 m的格網DEM數據，提取山脊線的實驗區域為左上經緯度是(100.85°E，30.18°N)右下經緯度是(100.99°E，30.04°N)的正方形區域，該區域位于四川西南部，最大高程為4 665 m,最小高程為2 640 m，地形起伏明顯.

2.2 對比與分析

采用流水物理模擬法提取山脊線時如何連接分水嶺離散點比較困難，因為這些分水嶺離散點只有流出的流向沒有流入的流向，所以不能像連接山谷線那樣按照流向來連接.此外最速上升法算法簡單且容易實現，該算法是從沒有連接過的且高程值最小的山脊點開始連接，查看連接點的八鄰域高程值并找出高程值最大的點，如果該高程值最大的點沒有連接過且高程大于連接點，則把該點加入到連接的山脊線中，否則結束連接.但是山頂點的高程值大于其八鄰域內點的高程值，所以最速上升法連接山脊線的過程遇到山頂點必然發生斷裂，這不僅導致主要山脊線沒有連接起來，而且影響后續的山脊線之間連接處理.

為了解決分水嶺不容易連接成山脊線和最速上升法提取山脊線容易發生斷裂的問題，本文是先從分水嶺的離散點中確定山頂點和鞍部點如圖6(a)所示，其中三角形的點為山頂點圓形的點為鞍部點，在山頂點和鞍部點的連接方向的控制下連接分水嶺的離散的點成主要的山脊線如圖6(b)所示，剩余的分水嶺的離散點通過改進的最速上升法連接，最后對次要山脊線進行刪除短小山脊線和連接到主要山脊線的處理，其提取結果如圖7(a)所示，圖7(b)中的山脊線為最速上升法提取的山脊線，從圖7(b)可知最速上升法提取的山脊線比較破碎，而且提取的山脊線沒有主次之分，主要山脊線和短小的次要山脊線同樣寬度.從圖7中(a)和(b)的對比可知，本模型提取的山脊線相較最速上升法連續且主要山脊線和次要山脊線在寬度上區別明顯，同時通過暈渲圖7(a)中的等高線可知本模型提取的山脊線與山脊的走勢相符即與實際地形較好的符合.

圖6 本模型約束連接方向的控制點與提取的主要山脊線Fig.6 Control points and major lines

圖7 本模型與最速上升法提取山脊線對比圖Fig.7 Comparison with steepest ascent method

3 結束語

本文提出的約束連接方向最速上升法提取存在于流水物理模擬法確定的分水嶺中的山脊線，通過分水嶺中的山頂點、鞍部點和山脊線最后兩個點約束山脊線的大致連接方向.與最速上升法提取分水嶺中的山脊線相比，本文提取的山脊線在連續性、層次性和整體性都有所提高，同時也解決了流水物理模擬法連接分水嶺中的山脊線的難題.

值得注意的是，本文提出的約束連接方向最速上升法也存在很大的改進空間，連接主要山脊線時只通過山頂點和鞍部點約束山脊線的連接方向，控制點不應該僅限于山頂點和鞍部點還需考慮更多的特征點類型，同時在連接主要的山脊線時需要考慮更多的情況.

[1] O’CALLAGHAN J F, MARK D M. The extraction of drainage networks from digital elevation data[J]. Computer vision, graphics, and image processing.1984, 28(3): 323-344．

[2] 羅 偉， 鄢志武， 薛重生. 基于DEM的巖溶地區地形地貌特征提取與分析[J]. 華中師范大學學報(自然科學版). 2013(03): 436-439．

[3] JENSON S K. Automated derivation of hydrologic basin characteristics from digital elevation model data[C]. 1985. 301-310．

[4] JENSON S K, DOMINGUE J O. Extracting topographic structure from digital elevation data for geographic information system analysis[J]. Photogrammetric engineering and remote sensing. 1988, 54(11): 1593-1600．

[5] 郭明武， 吳 凡. 基于規則格網DEM自動提取地性線的一種簡便方法[J]. 測繪信息與工程, 2005(2): 33-34．

[6] 黃培之. 提取山脊線和山谷線的一種新方法[J]. 武漢大學學報(信息科學版). 2001(3): 247-252．

[7] YOELI P. Computer-assisted determination of the valley and ridge lines of digital terrain models[J]. International Yearbook of Cartography. 1984, 24: 197-206．

[8] KOKA S, ANADA K, NOMAKI K, et al. Ridge Detection with the Steepest Ascent Method[J]. Procedia Computer Science, 2011, 4: 216-221．

[9] KOKA S, ANADA K, NAKAYAMA Y, et al. A Comparison of Ridge Detection Methods for DEM Data[C]// IEEE, 2012. 513-517．

[10] 朱 慶， 趙 杰， 鐘 正， 等. 基于規則格網DEM的地形特征提取算法[J]. 測繪學報, 2004(1): 77-82．

[11] GARBRECHT J, MARTZ L W. The assignment of drainage direction over flat surfaces in raster digital elevation models[J]. Journal of hydrology, 1997, 193(1-4): 204-213．

[12] BAI R, LI T, HUANG Y, et al. An efficient and comprehensive method for drainage network extraction from DEM with billions of pixels using a size-balanced binary search tree[J]. Geomorphology, 2015, 238: 56-67．

[13] ZHANG H, HUANG G. Building channel networks for flat regions in digital elevation models[J]. Hydrological Processes, 2009, 23(20): 2879-2887．

[14] BARNES R, LEHMAN C, MULLA D. An efficient assignment of drainage direction over flat surfaces in raster digital elevation models[J]. Computers & Geosciences, 2014, 62: 128-135．

[15] JANA R, RESHMIDEVI T V, ARUN P S, et al. An enhanced technique in construction of the discrete drainage network from low-resolution spatial database[J]. Computers & geosciences, 2007, 33(6): 717-727．

[16] FREEMAN T G. Calculating catchment area with divergent flow based on a regular grid[J]. Computers & Geosciences, 1991, 17(3): 413-422．

[17] WOOD J. The geomorphological characterisation of digital elevation models[Z]. University of Leicester (United Kingdom), 1996.465-479.

Detecting ridge lines with the constrained connecting direction steepest ascent method

DENG Xiaolong， JIANG Wenping， SUN Eryang， ZHANG Ping， YU Chi

(School of Resource and Environmental Sciences， Wuhan University， Wuhan 430079)

In this paper， methods for extracting ridge lines from grid digital elevation models are studied. To overcome shortcomings of the overland flow simulation method and the steepest ascent method， the constrained connecting direction steepest ascent method is proposed. Hill points and saddle points are defined as control points. Connecting direction is constrained to improve extraction. The results show that the method is effective for decreasing the fragmentation， cross and rotation of extracted ridge lines. Moreover， it has better performance in the continuity and integrity of ridge lines than the two basic methods mentioned above．

Grid DEM; ridge lines; constrained connecting direction; overland flow simulation; steepest ascent method

2015-06-03.

國家自然科學基金(41371428；41201473);國家基礎科學人才培養基金(J1103409).

1000-1190(2015)06-0979-05

P283.8; TU198+.5

A

*通訊聯系人. E-mail: 0020377@whu.edu.cn.