基于流域邊界分析的山脊線提取方法研究

黃萌萌，李 燦，郭鳳嬌

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

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

1 研究原理與方法

1.1 流域分割

1.1.1 流域及其組成

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

圖1 流域的組成示意圖

1.1.2 流域邊界搜索

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

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

1.2 山脊線的提取原理

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

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

1.2.2 山脊線精細(xì)化處理

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

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

表1 圖層L的主要字段

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

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

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

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

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

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

2.2 實(shí)驗(yàn)結(jié)果對(duì)比分析

2.2.1 河網(wǎng)的提取

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

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

2.2.2 流域邊界的提取

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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