一種在DEM平地中的河網(wǎng)水流方向估算方法

董良，張宏鳴，王巖澤，詹淇茹，樊曉

(西北農(nóng)林科技大學(xué) 信息工程學(xué)院，陜西 楊凌 712100)

0 引言

在區(qū)域尺度研究中，基于數(shù)字高程模型(digital elevation model，DEM)獲取河網(wǎng)、匯水面積、坡度等相關(guān)地形特征時，都需要提取水流方向[1]。通常，原始DEM經(jīng)填洼處理[2]，加之地形本身的特征，平地面積會占有很大的比例。由于平地高程一致，快速、有效獲取其流向一直是比較棘手的問題[3]。在國內(nèi)外關(guān)于流向的計算方法研究中，O’Callaghan等[4]提出的坡面流模擬方法D8(deterministic eight-neighbors)算法應(yīng)用最廣，即在3像素×3像素的DEM 柵格上，計算中心柵格與各相鄰柵格間的距離權(quán)重高程差，差值最大的柵格為中心柵格的流出柵格。而在平地流向提取時，無法明確水流方向。

關(guān)于平地流向的計算，一些研究者先后提出了各自的方法。Jenson等[5]提出的方法(J&D)在后來許多研究中得以應(yīng)用，也被集成到ArcGIS等軟件中。該方法根據(jù)平地出口點(diǎn)向平地內(nèi)部迭代確定平地柵格流向，每次迭代需遍歷所有DEM柵格，處理效率較低，且結(jié)果易產(chǎn)生“平行流”問題[6]。Martz等[7]采用高程增量疊加的方式修改平地DEM，確定研究區(qū)域內(nèi)的水流方向。之后，Garbrecht等[8]又改進(jìn)了高程增量的計算方法，G&M方法有效避免了平行流問題。

Barnes等[9]基于G&M方法將平地的增高處理分為3個部分，先計算平地與出口的距離，再計算平地與邊界的距離，最后疊加這2種距離，得到平地高程增量。在此基礎(chǔ)上，又有研究者做了進(jìn)一步推進(jìn)，如使用鏈碼矩陣(CCM)[10]和距離變換(DT)[11]方法來確定平地區(qū)域的流向。Liu等[12]采用3種不同的高程增量修正平地，第一級和第二級增量分別用來確定平地出口點(diǎn)和平地中心的流向，第三級增量通過重新修正區(qū)域內(nèi)高程值來確定流向。盡管減少了DEM增量，但反復(fù)修正DEM值，效率還需進(jìn)一步提升。此外，修正DEM的方法也會破壞地形形態(tài)，影響地形特征的提取質(zhì)量[13]。

國內(nèi)很多研究者還提出了其他方法來確定平地流向。Kong等[14]認(rèn)為在輸入的等高線資料中，根據(jù)地形圖資料增加能夠反映地表起伏的高程信息，避免平地和洼地出現(xiàn)，但這會使地形地貌發(fā)生形變，影響地形特征的提取質(zhì)量。Zhang等[15]采用徑向基函數(shù)插值法確定每個柵格的高程增量，利用數(shù)字信道網(wǎng)絡(luò)結(jié)合周圍地形的高程來確定平地流向。Pan等[16]通過線性插值的方法來修正DEM值，從而用D8算法計算平地流向。趙美玲[17]提出在平坦區(qū)域結(jié)合實際水系，修正主干河道高程，使整個區(qū)域水流有序向流域出口匯流。Yu等[18]改進(jìn)了Wang等[19]方法中由于填洼重新生成平地的問題，并且通過修改平地DEM值計算其流向。此外，Wu等[20]在綜合了FD8算法[21]和Rho8算法[22]的基礎(chǔ)上，從多流向和地形隨機(jī)性因素2個方面對流向算法進(jìn)行了改進(jìn)。Zhang等[23]提出了結(jié)合非平地區(qū)域流累積值的MW方法來確定流動方向。夏譽(yù)玲等[24]提出的混合流向算法將數(shù)字地形分類后，對緩坡和鞍部區(qū)域采用多流向算法確定徑流方向，并進(jìn)行水量分配，能大幅減少非自然的平行徑流。

上述方法大多參考平地周圍高程值，獲得的流向幾乎一致，易產(chǎn)生平行的流線，無法提取清晰而獨(dú)立的河網(wǎng)。而對平地采用逐級抬高增加高程來構(gòu)造微坡面的方法，使得平地流向處理效率不高，并且DEM的微小變形可能在平地周圍出現(xiàn)洼地，對結(jié)果產(chǎn)生較大影響[13]。

本文在前期研究的基礎(chǔ)上[23]，提出以平地柵格的初始匯流累積量為參考，采用最短路徑方法和就近原則方法來計算平地的水流方向。本研究直接對平地進(jìn)行流向判斷，不需對DEM進(jìn)行修改，通過與J&D方法、Barnes方法對比，用所提出方法提取流向結(jié)果合理并效率較高，具有較好的運(yùn)行效果和可行性，是對數(shù)字地形分析方法的有益嘗試。

1 平地處理的方法介紹

1.1 平地處理的基本思想

在數(shù)字流域特征提取中，形成數(shù)字河網(wǎng)的關(guān)鍵環(huán)節(jié)之一是獲取每個DEM柵格的單位匯水面積[25]。單位匯水面積是指其他柵格直接或間接流向一個柵格的匯流累積量(即累計流量)。通過計算非平地區(qū)域的匯水面積，平地最外圈的柵格會有初始累計流量，而平地內(nèi)部柵格的初始累計流量則為0，這些平地區(qū)域便是河網(wǎng)提取中存在的斷點(diǎn)。

一般在流域中，上游柵格的累計流量按照流向趨勢逐漸向下游累積，逐漸增大，最終會形成一個數(shù)字河網(wǎng)。因此，在平地上游存在一個在平地區(qū)域中初始累計流量最大的平地柵格，將其視為平地水流的一個入口點(diǎn)；在平地下游有一個高程值大于周圍非平地柵格高程的平地柵格，即有流向的平地柵格，將其視為平地的一個出口點(diǎn)。

為提取平地區(qū)域內(nèi)的流向，可將其入口點(diǎn)和出口點(diǎn)之間的最短路徑視為主流，從而確定主流路徑上各柵格的流向。同樣，將平地上游中初始累計流量次大的平地柵格作為支流入口點(diǎn)，計算該入口點(diǎn)到達(dá)平地中已分配流向柵格之間的最短路徑，確定該路徑上平地柵格的流向，該方法稱為最短路徑方法。需要注意的是，在某一平地區(qū)域中支流的個數(shù)不宜太多，否則會出現(xiàn)區(qū)域內(nèi)流向散亂的情況，無法獲得清晰的河網(wǎng)。而對于平地中未分配流向的柵格，或初始累計流量小于規(guī)定閾值的小流量柵格，使用就近原則方法，將其流向收斂于附近已有流向的平地柵格。

本文將前期計算DEM數(shù)據(jù)的平地柵格的初始匯流累積量和確定的平地出口點(diǎn)，作為使用以上2種方法的基本前提條件。由出口點(diǎn)找到待處理的平地，由初始累計流量的大小來確定使用最短路徑方法或者就近原則方法。

1.2 平地處理的過程

平地處理總體流程如圖1所示。圖2表示某平地區(qū)域及周邊非平地柵格的高程值。在平地處理前，需要使用D8算法確定非平地的流向，并找到平地的出口點(diǎn)，即平地下游邊緣區(qū)域周邊柵格高程值，低于該平地柵格高程值的平地點(diǎn)(圖3)，再由非平地流向計算其累計流量，便得到平地邊緣柵格的初始累計流量。規(guī)定區(qū)分流量大小的閾值為8倍柵格長度，將平地中初始累計流量較大的柵格，依次作為平地區(qū)域的主流和不超過3條支流的入口點(diǎn)(圖4)。對于平地中剩余的小流量柵格，采用就近原則方法，賦予其附近流域的流向。

圖4 平地初始累計流量

圖3 D8算法確定平地出口

圖2 平地區(qū)域及周邊柵格高程值

圖1 平地處理總體流程

將這塊平地中的所有柵格均處理完，再以相同的方法依次處理所有的平地。

1.3 大流量柵格最短路徑算法

平地區(qū)域主流和支流的最短路徑是一個簡化的迷宮求解[26]的過程。從平地入口點(diǎn)開始按照寬度優(yōu)先搜索的方法走到出口點(diǎn)，然后回溯確定路徑，并分配給路徑上各平地柵格的流向。所求最短路徑首先滿足包含最少步數(shù)的條件。

假設(shè)某2點(diǎn)間的路徑由m條長度為1的直邊和n條長度為1的斜邊組成，該路徑長度為L1；每減少k步斜邊則增加2k步直邊，得到路徑L2；則有表達(dá)式如式(1)、式(2)所示。

(1)

(2)

圖7 支流路徑確定過程

圖6 主流路徑確定過程

圖5 平地主流支流的路徑

2)按照寬度優(yōu)先搜索方法依次遍歷每個節(jié)點(diǎn)，并計入輔助隊列中。

3)按照步驟1)方法確定當(dāng)前節(jié)點(diǎn)的子節(jié)點(diǎn)，以步驟2)方法不斷遍歷至平地的出口點(diǎn)(即有流向的平地柵格)后停止搜索。接下來從出口點(diǎn)回溯，尋找當(dāng)前節(jié)點(diǎn)的父節(jié)點(diǎn)，并賦予其水流方向。

4)找到平地的主流之后，使用上述最短路徑算法，從初始累計流量次大的平地柵格搜索到距離最近的主流路徑上或出口點(diǎn)，確定支流及路徑上的柵格流向。

1.4 小流量柵格就近原則算法

通過最短路徑算法確定平地區(qū)域中主流和支流路徑上部分柵格的流向之后，剩余未分配流向的都是初始累計流量較小的柵格。這些小流量的柵格在匯水過程中，不能形成獨(dú)立的水流。為避免區(qū)域內(nèi)出現(xiàn)散流、亂流現(xiàn)象，提出就近原則算法，將剩余未分配流向的平地柵格，收斂于附近的主流或者支流上。就近原則算法的具體步驟描述如下。

1)設(shè)置一個數(shù)組依次存放平地中已有流向的柵格。

2)遍歷該數(shù)組中有流向平地柵格的8鄰域，如果鄰域中有尚未分配流向的柵格，則令其流向8鄰域中間的平地柵格，并將此平地柵格加入該數(shù)組中。

3)按照步驟2)方法依次遍歷數(shù)組中有流向的柵格，來確定其他未分配流向柵格的流向。當(dāng)數(shù)組中有流向平地柵格數(shù)等于該平地總柵格個數(shù)時，處理完畢。圖8(a)、圖8(b)、圖8(c)為就近原則方法的示意過程。

圖8 無流向柵格就近原則

2 平地流向處理的實驗與結(jié)果分析

2.1 實驗數(shù)據(jù)和環(huán)境

使用4個地區(qū)不同規(guī)模的數(shù)據(jù)進(jìn)行測試和分析。第1個是陜西省綏德縣韭園溝流域，海拔高度851～1 103 m，平均高度980 m，占地總面積100 km2。第2個是陜西省安塞縣縣南溝流域，流域面積44 km2，海拔1 100～1 400 m，平均高度1 200 m。這2個區(qū)域地貌溝壑縱橫，是典型的黃土丘陵溝壑區(qū)，侵蝕較嚴(yán)重。第3個是USGS官網(wǎng)上提供的SRTM1數(shù)據(jù)，澳大利亞中部(23.00°S，131.00°E～20.00°S，134.00°E)的9塊SRTM數(shù)據(jù)，該區(qū)域中平地面積所占比重很大。第4個是貝加爾湖以東(47.00°N，112.00°E～56.00°N，121.00°E)的100塊SRTM數(shù)據(jù)。

本算法在Visual Studio 2017集成開發(fā)環(huán)境下采用C#語言實現(xiàn)。運(yùn)行于Windows10(64位)操作系統(tǒng)，CPU為Interl(R)Core(TM)i7，內(nèi)存8 GB。為了和具有代表性的J&D方法、Barnes方法的計算結(jié)果進(jìn)行比較，要確保程序運(yùn)行環(huán)境一致。Barnes算法從其提供的下載站點(diǎn)獲取，同本算法一樣，每組DEM數(shù)據(jù)測試5次，取得平均運(yùn)行時間，輸入相同的DEM數(shù)據(jù)，輸出計算后的柵格流向存儲在一個文本文件中。將結(jié)果都導(dǎo)入ArcGIS 10.2軟件中，用水文分析工具處理得到可視化的圖層。J&D方法已集成在ArcGIS軟件的水文分析步驟中，直接使用工具完成。

2.2 平地流向處理結(jié)果

實驗中使用的10 m分辨率的韭園溝流域數(shù)據(jù)，局部溝谷存在連續(xù)的平坦區(qū)域，地形特征最具代表性。在平地處理前，直接對該地區(qū)進(jìn)行流向和匯水的計算，在設(shè)置河網(wǎng)提取的閾值為500后，得到有斷點(diǎn)的河網(wǎng)(圖9)。經(jīng)本算法處理平地后，確定平地區(qū)域的流向，進(jìn)而計算匯水，可提取到完整的河網(wǎng)(圖10)。

圖9 韭園溝流域局部有斷點(diǎn)的河網(wǎng)

圖10 韭園溝流域局部完整的河網(wǎng)

2.3 同現(xiàn)有方法的對比和分析

本文和J&D方法、Barnes方法在處理結(jié)果的合理性和效率方面做出了比較和分析。為驗證本文方法的普適性，從澳大利亞中部地區(qū)一塊SRTM數(shù)據(jù)中分割出500像素×500像素的樣區(qū)進(jìn)行詳細(xì)分析對比。該樣區(qū)地形起伏較小，樣區(qū)內(nèi)存在較多平坦區(qū)域和緩坡(圖11)。根據(jù)本文方法，對該樣區(qū)進(jìn)行流向計算，然后設(shè)置閾值提取河網(wǎng)(圖12)。在相同河網(wǎng)閾值下，對利用J&D方法得到的結(jié)果(圖13)以及Barnes方法得到的結(jié)果(圖14)進(jìn)行了分析和比較。

圖11 澳大利亞中部地區(qū)500像素×500像素的樣區(qū)

圖12 本文方法對樣區(qū)的河網(wǎng)提取結(jié)果

圖13 J&D方法對樣區(qū)的河網(wǎng)提取結(jié)果

圖14 Barnes方法對樣區(qū)的河網(wǎng)提取結(jié)果

J&D方法主要基于坡面流模擬D8算法，根據(jù)平地出口點(diǎn)向平地內(nèi)部迭代來確定平地流向。而D8算法在實際執(zhí)行過程中會規(guī)定一個8鄰域的遍歷順序，最終確定的所有平地柵格的流向往往受遍歷順序的影響，所得平地流向方向相同，易產(chǎn)生大面積的平行偽河道。Barnes方法修正平地DEM數(shù)據(jù)，為平地設(shè)置高程增量，將其改為微坡面，然后依舊用D8算法確定流向。整個修正過程分為3步：①計算平地內(nèi)部每個柵格與出口的距離；②計算每個柵格到邊界的距離；③將2種距離進(jìn)行累加，確定最終的平地高程增量。將平坦的區(qū)域重構(gòu)為周邊略高中心低、出口邊更低的一個“微型溝谷匯水面”。這種方式在徑流匯集處會得到不錯的效果，減少了部分的平行徑流的出現(xiàn)，但在平地的中心區(qū)域仍會有小面積的平行偽河道。

本文方法的提取結(jié)果其河網(wǎng)總體走勢和上述2種方法的結(jié)果基本一致，而設(shè)置相同河網(wǎng)閾值后，得到的河網(wǎng)流線清晰，并無平行流產(chǎn)生。這是由于平地流向算法首先確定主流和支流最短路徑，保證了水流路線的整體趨勢，而對于小流量柵格的流向，采用就近原則收斂于主流支流，不會出現(xiàn)流向散亂的情況，總體分配合理。

在算法運(yùn)行效率方面，本研究設(shè)計的算法需要前期計算DEM數(shù)據(jù)的非平地區(qū)域的匯水面積，同時找到所有平地的出口點(diǎn)以及對應(yīng)的平地。在設(shè)計實現(xiàn)時，將計算DEM柵格數(shù)據(jù)的匯水面積和計算平地區(qū)域水流方向相結(jié)合起來，這里只分析平地處理部分。對于Barnes方法，相應(yīng)地，也僅分析平地處理過程的效率。對不同地區(qū)、不同規(guī)模的數(shù)據(jù)，算法具體運(yùn)行時間如表1所示。

通過使用不同特征的DEM數(shù)據(jù)進(jìn)行測試，將Barnes方法和本文方法對比。從表1中分析發(fā)現(xiàn)，對于同一地區(qū)的不同分辨率的數(shù)據(jù)，分辨率下降時，DEM數(shù)據(jù)存儲量減少，對應(yīng)的平地數(shù)量相對減少，本文方法對所有測試數(shù)據(jù)的處理時間均小于Barnes方法，效率較高。Barnes方法對平地的處理過程主要在于修正DEM數(shù)據(jù)，需要計算2次高程增量再疊加，每次計算還需要遍歷所有平地柵格；而本文方法直接進(jìn)行計算，遍歷一次平地柵格即可。在測試數(shù)據(jù)比較龐大或區(qū)域中平地所占面積比較大時，Barnes方法對內(nèi)存空間的要求更高，這也是增加計算量產(chǎn)生的結(jié)果。本文方法計算簡單，可行性較好。

表1 不同數(shù)據(jù)、不同規(guī)模下的平地處理效率對比

3 結(jié)束語

DEM數(shù)據(jù)預(yù)處理是水文應(yīng)用中的一個重要環(huán)節(jié)，只有合理確定水流方向，才能得到精度較高的河網(wǎng)水系及數(shù)字流域特征，為分布式水文模擬提供基礎(chǔ)數(shù)據(jù)。基于DEM的水系提取，在地形起伏較大的山區(qū)等非平緩地區(qū)應(yīng)用精度較高，但在平緩區(qū)域或人類活動影響較大的平原區(qū)還有待進(jìn)一步探索，這也將該領(lǐng)域的研究推向熱潮。

借鑒前人研究，在對流域分布式侵蝕學(xué)平地處理認(rèn)識的基礎(chǔ)上，進(jìn)行理論分析和研究，設(shè)計并實現(xiàn)了平地水流方向處理的算法。本文主要與J&D方法、Barnes方法進(jìn)行對比。相比而言，本算法得到河網(wǎng)的結(jié)果與以上2種方法基本一致，且?guī)缀醪淮嬖诓环蠈嶋H規(guī)律的平行流的問題。在運(yùn)行時間上，Barnes方法通過修正平地DEM數(shù)據(jù)，2次計算平地的高程增量，再進(jìn)一步將2種距離累加確定最終的平地高程增量，然后用D8算法確定修正平地后的微坡面的流向，這樣不斷重復(fù)計算和迭代，計算量很大，效率較低。而經(jīng)過本算法的計算，可以明顯看到平地中的主流，以及其他支流匯聚到主流之上，也解決了平地區(qū)域中部分水流方向散亂的問題，整體結(jié)果比較合理，無平行流產(chǎn)生。在運(yùn)行時間上，本算法無須修正平地DEM數(shù)據(jù)，通過前期計算DEM數(shù)據(jù)的非平地區(qū)域的累計流量，同時找到平地的出口點(diǎn)和部分入口點(diǎn)，直接對其進(jìn)行計算，運(yùn)行時間短，效率較高。

綜上所述，本處理方法利用DEM對平緩地區(qū)的數(shù)字流域構(gòu)建，是平地河網(wǎng)提取方法的有益嘗試，也是對數(shù)字地形分析技術(shù)的補(bǔ)充。