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

(武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院， 武漢 430079)

?

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

(武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院， 武漢 430079)

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

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

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

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

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

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

1.1 提取主要山脊線

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

圖1 山頂點(diǎn)和鞍部點(diǎn)Fig.1 Hill and saddle points

圖2 山頂點(diǎn)的坐標(biāo)系Fig.2 Coordinate system of hill point

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

圖3 連接山脊點(diǎn)的坐標(biāo)系Fig.3 Coordinate system of connection point

圖4 連接點(diǎn)的八鄰域Fig.4 Eight neighborhoods

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

表1 山脊線連接點(diǎn)的優(yōu)先級匯總表

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

1.2 次要山脊線的提取

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

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

2 實驗與分析

2.1 實驗過程與數(shù)據(jù)

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

① DemMatrixPre=PretreatDemdata(DemMatrix);

② HillSaddleArray=LocateHillSaddle(DemMatrixPre);

③ MajorRidgeArray=ConnectMajorridge(HillSaddleArray);

④ MiniorRidgeArray=ConnectMiniorridge();

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

2.2 對比與分析

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

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

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

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

3 結(jié)束語

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

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

[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的巖溶地區(qū)地形地貌特征提取與分析[J]. 華中師范大學(xué)學(xué)報(自然科學(xué)版). 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] 郭明武， 吳 凡. 基于規(guī)則格網(wǎng)DEM自動提取地性線的一種簡便方法[J]. 測繪信息與工程, 2005(2): 33-34．

[6] 黃培之. 提取山脊線和山谷線的一種新方法[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版). 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] 朱 慶， 趙 杰， 鐘 正， 等. 基于規(guī)則格網(wǎng)DEM的地形特征提取算法[J]. 測繪學(xué)報, 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.

國家自然科學(xué)基金(41371428；41201473);國家基礎(chǔ)科學(xué)人才培養(yǎng)基金(J1103409).

1000-1190(2015)06-0979-05

P283.8; TU198+.5

A

*通訊聯(lián)系人. E-mail: 0020377@whu.edu.cn.