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

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

0 引言

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

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

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

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

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

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

1 平地處理的方法介紹

1.1 平地處理的基本思想

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

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

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

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

1.2 平地處理的過程

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

圖4 平地初始累計流量

圖3 D8算法確定平地出口

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

圖1 平地處理總體流程

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

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

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

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

(1)

(2)

圖7 支流路徑確定過程

圖6 主流路徑確定過程

圖5 平地主流支流的路徑

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

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

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

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

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

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

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

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

圖8 無流向柵格就近原則

2 平地流向處理的實驗與結果分析

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

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

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

2.2 平地流向處理結果

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

圖9 韭園溝流域局部有斷點的河網

圖10 韭園溝流域局部完整的河網

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

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

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

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

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

圖14 Barnes方法對樣區(qū)的河網提取結果

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

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

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

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

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

3 結束語

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

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

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