黃萌萌，李 燦，郭鳳嬌

（1.重慶市勘察規(guī)劃設計有限公司，重慶 401121；2.重慶欣榮土地房屋勘測技術研究所有限責任公司，重慶 401121；3.重慶市巴南區(qū)規(guī)劃和自然資源局，重慶 401000）

地形中沿山脊走向布設的路線即為山脊線。山脊線與山谷線作為地形骨架線，共同決定著地形的高低起伏。在地形表達中，山脊線和山谷線對其他地理要素起著控制性的作用。因此，山脊線常被作為專題進行研究，山脊線的提取和分級方法則是地形綜合、制圖表達等領域的重要研究內(nèi)容。目前，根據(jù)基礎數(shù)據(jù)，山脊線的提取方法可分為基于等高線[1]的方法和基于DEM分析的方法?；贒EM分析的方法主要包括地表徑流模型（反地形）[2-3]和基于流域分析[4]的提取方法，前者已有較成熟的算法，并集成于GIS軟件中，可快速自動化提取山脊線，但存在無法精確提取平坦地區(qū)山脊線，且出現(xiàn)平行效應的問題；后者是由郭萬欽[4]等提出的一種利用流域邊界和坡向差提取山脊線的方法，能較精準地提取山脊線，有效解決平坦地區(qū)的平行效應，但該方法需要設定匯流累積量、坡向差、坡向距離等6個閾值參數(shù)，人工干預性較大、自動化程度較低。針對上述問題，本文提出了一種基于流域邊界與河網(wǎng)空間關系分析的山脊線提取方法，可有效避免平坦地區(qū)的平行效應問題，解決了閾值參數(shù)過多的問題，提高了山脊線提取的自動化效率。

1 研究原理與方法

1.1 流域分割

1.1.1 流域及其組成

流域是指流經(jīng)其中的水流或其他物質(zhì)從不同的方向流過公共的出水口，這些水流或其他物質(zhì)流經(jīng)的由分水線所包圍而形成的閉合集水區(qū)域[5-6]。從尺度上看，不同的水系等級都有自己的流域，流域本身又可按照水系級別細分為多個子流域。相鄰流域之間公用的邊界為分水嶺，也是集水區(qū)域的邊界線。水流或其他物質(zhì)流經(jīng)出水口的路徑形成河網(wǎng)，從水流方向定向的特點來看，河網(wǎng)的形態(tài)為樹結構，樹的根節(jié)點是出水口，即為該流域內(nèi)的最低點[7-9]。流域的組成如圖1所示。

圖1 流域的組成示意圖

1.1.2 流域邊界搜索

流域分割是以出水點（該集水區(qū)域內(nèi)的最低點）為基礎，結合水流方向矩陣，逆流而上搜索該出水點上游所有流經(jīng)該點的柵格單元，直到該集水區(qū)域內(nèi)所有柵格單元都確定了位置，即搜索到指定出水口點的流域邊界[10]。為了確保山脊線網(wǎng)絡與河網(wǎng)的匹配性，可采用河流鏈接數(shù)據(jù)（Stream Link）作為出水口數(shù)據(jù)，其中隱含著每條河網(wǎng)弧段的聯(lián)結信息，包括河段的起點和終點，河段的終點即為該匯水區(qū)域出水口所在的位置[11-12]。同時，河網(wǎng)提取采用的閾值也決定了該河段所處匯水區(qū)域的詳略程度，即提取河網(wǎng)的閾值越大，子流域面積也越大。

如圖2所示，通過D8算法計算無洼地DEM每個柵格單元的水流方向，并計算流入每個柵格單元的匯流總量，獲得匯流累積量矩陣；圖2d為將柵格單元(4,4)指定為出水口點，逆向搜索所有流經(jīng)該點的柵格單元（圖中的藍色單元），所有藍色單元形成了出水口(4,4)的流域邊界。

1.2 山脊線的提取原理

1.2.1 流域邊界到山脊線的轉(zhuǎn)換

流域邊界除了分水嶺界線外，還有一部分是位于負地形的邊界線。因此，可利用流域邊界線與山脊線之間的包含關系，剔除位于平坦河谷地帶等負地形的流域邊界，保留的則是需要提取的山脊線網(wǎng)絡。從流域邊界與河網(wǎng)的空間位置關系來看，位于負地形的流域邊界與河網(wǎng)相交，可通過以下思路粗略提取山脊線：首先將柵格的流域邊界轉(zhuǎn)換為矢量數(shù)據(jù)，并由面圖層轉(zhuǎn)換為線圖層；然后合并所有的線要素，并在節(jié)點（線與線的交點）處打斷；最后通過河網(wǎng)空間位置查詢，選取位于負地形的流域邊界線要素，刪除該部分要素，形成粗略的山脊線網(wǎng)絡。

1.2.2 山脊線精細化處理

1）閉合環(huán)檢查。受DEM分辨率、流域提取算法或其他不明原因的影響，刪除負地形的流域邊界后，山脊線網(wǎng)絡內(nèi)部尚存在閉合環(huán)。為提高山脊線的使用效率，確保山脊線網(wǎng)絡的拓撲正確性，應對網(wǎng)絡內(nèi)部的閉合環(huán)進行破環(huán)處理。對刪除位于負地形的邊界后剩下的流域邊界（簡稱圖層L）進行要素轉(zhuǎn)面操作，可將所有的閉合圖形轉(zhuǎn)為面要素，未形成閉合圖形的線要素則不能轉(zhuǎn)面成功，形成面圖層P。面圖層P中的每個面要素則是待處理的環(huán)，對每個面要素的RID進行編號，RID是每個環(huán)的唯一標識碼。

2）閉合環(huán)分析。由于面圖層P與圖層L中的部分要素存在共線的關系，將面圖層P的唯一編號RID通過空間鏈接賦值到圖層L，空間位置關系為：SHARE_A_LINE_SEGMENT_WITH。若鏈接要素中具有與目標要素共線的要素，則匹配這些要素，賦值到圖層L的RID字段，并將空間鏈接成功要素的BZ字段賦值為環(huán)要素，其余未鏈接成功的賦值為非環(huán)要素，RID為0。圖層L的主要字段如表1所示。

表1 圖層L的主要字段

根據(jù)環(huán)要素與非環(huán)要素的空間位置關系，可將環(huán)分為兩種類型：①多要素環(huán)，由多條要素組成，每個節(jié)點由環(huán)與非環(huán)要素相連接，如圖3a所示，多邊形（環(huán)）ABCDE由邊AB、BC、CD、DE、EA組成，每個 節(jié)點均是環(huán)與非環(huán)要素的交點；②單要素環(huán)，由一條邊形成，單環(huán)與非環(huán)要素只有一個交點，如圖3c所示，多邊形與非環(huán)要素相交于V0點，V0既是邊界的起點，也是終點。

圖3 多要素環(huán)和單要素環(huán)破環(huán)前后示意圖

3）破環(huán)方法。通過對圖層L的RID字段進行頻數(shù)分析，可識別位于同一環(huán)的各線要素，RID相同的要素為同一環(huán)的邊界，RID的頻數(shù)可區(qū)分環(huán)的類型，即

為保持破環(huán)后，山脊線網(wǎng)絡整體的連通性、完整性和拓撲正確性，可分別對單要素環(huán)和多要素環(huán)進行以下處理：①針對多要素環(huán)由多條線要素構成的幾何特征，可對位于邊界上的各線要素進行對比分析，選取長度最長的要素（DE），并刪除；②針對單要素環(huán)由一條線要素組成的幾何特征，線要素的起點和終點重合，可從起點開始，計算每個節(jié)點至起點的歐氏距離，標記距離最長的節(jié)點P，刪除該節(jié)點的下一個節(jié)點至終點的所有節(jié)點，點P為終點。

2 實驗結果與分析

2.1 實驗數(shù)據(jù)

實驗區(qū)位于108.051°～109.493 E、35.635°～ 37.206 N之間，海拔高度為780～1 800 m。實驗數(shù)據(jù)來源于ASTER GDEM 30 m分辨率的DEM數(shù)據(jù)。本文提出的數(shù)據(jù)處理方法基于ArcGIS DeskTop10.2提供的工具箱以及ArcGIS Engine10.2，采用C#語言在Microsoft Visual Studio 2012開發(fā)平臺下實現(xiàn)。

2.2 實驗結果對比分析

2.2.1 河網(wǎng)的提取

首先計算實驗區(qū)的水流方向，利用ArcGIS Toolbox里的水文分析工具，雙擊Flow Direction工具，選擇實驗區(qū)域的DEM數(shù)據(jù)，生成水流方向圖層FlowDir；然后采用Sink工具對FlowDir進行洼地檢測，若存在洼地，則利用Fill工具對原始DEM進行填洼處理，若無洼地，則采用Flow Accumulation計算匯流累積量數(shù)據(jù)，生成FlowAcc數(shù)據(jù)。本文通過不斷試驗和利用現(xiàn)有地形圖等其他數(shù)據(jù)輔助檢驗的方法來確定能滿足研究需要且符合研究區(qū)域地形地貌條件的合適的閾值。為滿足不同尺度的河網(wǎng)和山脊線數(shù)據(jù)要求，選取了3組閾值（閾值T分別設置為500、200、100）進行實驗，結果如圖4所示。通過Raster Calculator工具，將匯流累計量大于等于閾值的屬性值設定為1，這些柵格就是河網(wǎng)的潛在位置；將小于閾值的屬性值設定為No Data，形成StreamNet柵格數(shù)據(jù)。采用Stream to Feature完成對河網(wǎng)的柵矢轉(zhuǎn)換，形成Stream矢量圖層。

圖4 不同閾值的河網(wǎng)提取結果

2.2.2 流域邊界的提取

首先雙擊Hydrology工具集中的Stream Link工具，在Input Stream Raster文本框中選擇StreamNet，在Input Flow Direction Raster文本框中選擇fdirfill，輸出數(shù)據(jù)為Stream Link；然后雙擊Watershed工具，打開集水區(qū)域（貢獻區(qū)域）計算的對話框，分別在水流方向數(shù)據(jù)和出水口數(shù)據(jù)的輸入文本框中選擇fdirfill和Stream Link數(shù)據(jù)，輸出數(shù)據(jù)為Watershed，并利用Raster to Polygon工具將Watershed轉(zhuǎn)換為矢量圖層。由圖5可知，流域的邊界能很好地概括山脊線的走向，并與相同閾值的河網(wǎng)數(shù)據(jù)形成良好的銜接關系。

圖5 基于不同閾值河網(wǎng)提取的流域邊界

2.2.3 流域邊界與河網(wǎng)空間分析

首先選擇閾值相匹配的河網(wǎng)S t r e a m圖層和Watershed圖層；然后通過Feature to Polyline將Watershed轉(zhuǎn)換為線圖層；最后將目標圖層設置為Watershed，源要素設置為Stream，空間位置關系選擇與源要素相交，通過空間位置查詢選取Watershed中位于負地形的流域邊界，并將其刪除，剩下的流域邊界則是粗略的山脊線Ridge_R。受DEM分辨率、流域邊界提取算法等因素的影響，提取的山脊線內(nèi)部存在部分小環(huán)，如圖6所示。

圖6 山脊線內(nèi)部存在的閉合環(huán)

為了獲取科學合理、實用性強的山脊線，需對粗略提取的山脊線網(wǎng)絡進行精細化處理。首先采用Feature to Polygon工具將Ridge_R圖層轉(zhuǎn)為面圖層，確定Ridge_R內(nèi)部的閉合環(huán)，并對所有環(huán)要素的RID字段編碼，可直接賦值為要素的FID+1；然后利用空間鏈接將RID鏈接回圖層Ridge_R，通過頻數(shù)分析確定閉合環(huán)的類型，并利用前文敘述的環(huán)處理思路，通過代碼對山脊線內(nèi)部環(huán)進行破環(huán)處理。由圖7可知，與傳統(tǒng)的基于反地形D8算法的提取方法相比，本文方法的實驗結果能有效地與地形相吻合，且在山頂平坦位置或洼地位置避免了因填洼出現(xiàn)的平行效應；同時，在參數(shù)設置和自動化提取方面，較參考文獻[4]也有一定的簡化。

圖7 D8算法與本文方法提取山脊線的對比圖

本文利用基于流域邊界分析的山脊線提取方法，采用3組閾值分析了提取的山脊線與相同密度的山谷線的銜接情況，如圖8所示，可以看出，山脊線與相同密度的山谷線的銜接較好，能構建出不同層次的地形網(wǎng)絡骨架，為有效劃分地形層次結構提供了實踐基礎。

圖8 本文算法提取的山脊線與相同閾值的山谷線銜接情況

3 結 語

針對山脊線提取存在的平行效應、自動化程度較低的問題，本文提出了基于流域邊界分析的山脊線提取方法。該方法能夠得到連續(xù)的山脊線，且與DEM可見的山脊線走向保持一致，避免了脊線段斷續(xù)和平坦地段出現(xiàn)平行效應等問題；同閾值山脊線與河網(wǎng)能相互有效銜接，形成可概括地形的骨架線網(wǎng)絡；不同匯流閾值的河網(wǎng)可獲得不同面積大小的流域邊界，從而形成詳略程度不一的山脊線網(wǎng)絡，為地形綜合研究提供了新的思路。目前還沒有科學有效的山脊線提取評價機制，山脊線的提取還應配合行之有效的分級方法，才能更好地為地形綜合和水文分析提供理論基礎。今后筆者將針對山脊線的分級做進一步研究。