柵格流向算法在山丘區(qū)水系提取中的應用評價

崔寅鶴，吳鵬飛，汪銀奎，劉金濤，3

(1.河海大學水文水資源學院，江蘇南京210098； 2.西藏自治區(qū)水文水資源勘測局，西藏拉薩850000； 3.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇南京210098)

0 引 言

地表流徑描述了水、土壤、污染物等受重力影響在陸地表面的運移路線[1]，在水文建模[2]、土壤侵蝕[3]等研究領域應用廣泛。除少數(shù)基于不規(guī)則三角網(wǎng)[4]、等高線[5]等數(shù)據(jù)的算法外，目前主流的流徑提取算法普遍以柵格數(shù)字高程模型(Digital elevation model，DEM)作為數(shù)據(jù)基礎[6- 7]。使用DEM提取流徑的算法普遍包括填洼、流向計算、匯流累積3個步驟[8]。在山丘區(qū)，流徑提取的結果主要受流向算法影響。流向算法分為單流向算法和多流向算法兩類[9]。以3×3網(wǎng)格為基礎，單流向算法將本柵格產(chǎn)生的水流和上游柵格匯流至本柵格的水流作為一個整體，認為其只能流向8個與本柵格相鄰柵格中的1個；多流向算法則將水流按比例分配給多個相鄰柵格。目前，應用較多的單流向算法包括D8算法[10]、D8-LTD算法[11]、GRASS軟件rwflood模塊使用的快速算法[12]等，多流向算法包括MFD算法[13]、Dinf算法[14]、MDinf算法[15]等。其中，Dinf算法雖然僅計算得到一個流向，但在進行匯流時需將流量按比例分配給兩個相鄰單元，故此處將其定義為多流向算法。

表2 流向算法原理

目前對各種流向算法的對比分析大多聚焦于D8、MFD、Dinf等經(jīng)典算法[16- 17]，對較新的rwflood、MDinf等算法的功能評價較少。本文選取3處典型高山、丘陵地形，使用劉學軍等[17]提出的定量化對比方法分析了6種流向算法(D8、D8-LTD、rwflood、MFD、Dinf、MDinf)的差異性，并對比了這些算法對真實河道的還原效果以及計算速率，以期為相關研究的流向算法選取提供參考。

1 材料與方法

1.1 研究材料

為了更好地研究流向算法在不同山丘區(qū)的應用效果，本文分別在青藏高原、東南丘陵、黃土高原選取了姜灣研究區(qū)、奪底溝研究區(qū)和杜家灣研究區(qū)(見圖1)。其中，奪底溝研究區(qū)坡度較大，平均坡度25.58°；姜灣研究區(qū)山地坡度居中，平均坡度20.23°；杜家灣研究區(qū)為黃土丘陵，山地坡度最小，平均坡度13.23°。

圖1 研究區(qū)位置及DEM數(shù)據(jù)

本文所選用的DEM數(shù)據(jù)源為SRTM V4.1 30 m分辨率數(shù)據(jù)。SRTM數(shù)據(jù)通過航天飛機搭載雷達進行測繪， V4.1版本采用了新型的插值方法，讓數(shù)據(jù)的準確度有了較大的提升。本文使用的全部數(shù)據(jù)下載自美國地質(zhì)調(diào)查局網(wǎng)站(http://earthexplorer. usgs.gov/)。

1.2 流向算法

流向算法旨在確定水流流動、土壤和污染物輸移路徑，體現(xiàn)在DEM數(shù)據(jù)中就是要解決各柵格向下游柵格分配流量的問題。基于這個問題可以將流向計算分為單流向算法和多流向算法兩大類。

1.3 研究方法

首先對各研究區(qū)DEM數(shù)據(jù)進行填洼(rwflood算法不需要填洼)，利用填洼后數(shù)據(jù)進行流向計算和回流累積；再根據(jù)匯流累積結果設定最大流向的1/100為閾值提取水系。本文使用的填洼算法是LIU等[18]提出的算法，對平坦地區(qū)賦予1×10-5m的梯度，匯流算法為Su等提出的基于樹狀圖的新算法[19]。這兩種算法均為目前計算速度較快的算法。

本文采用了劉學軍等提出的相似度量化分析方法，將匯水面積(Total Catchment Area， TCA)選為分析對象；并且將TCA進行了對數(shù)變換以避免由TCA分布不均勻造成誤差的現(xiàn)象[17]。此外，本文還對各算法在實際應用方面的效果進行了評估。

1.3.1累計頻率分布

在分析和表示地形屬性分布特征時，頻率和頻率分布圖是一種統(tǒng)計分析方法。TCA的數(shù)量表現(xiàn)和算法的不同通常會導致使用不同算法的TCA累計頻率出現(xiàn)差異。這種差異可以直觀地反映為頻率和頻率分布曲線形狀的不同。從差異特征中便可分析得出算法的機理以及TCA的空間分布。

1.3.2相對差系數(shù)

應用不同機理的徑流算法所得到的DEM在相同單元格內(nèi)TCA值存在差異。為了便于比較這種差異，故將其以相對差系數(shù)α來表示[17]。α的定義為

(1)

表3 各研究區(qū)TCA分布

當α≤0時，則兩種算法之間沒有可比性；當0<α<1時，兩種算法存在差異且α值越小差異越大；當α=1時，兩種算法之間沒有差異。由式(1)可看出，相對差系數(shù)是一種可將不同算法進行兩兩相比的數(shù)值指標，而在計算過程中需要假定其中一種算法作為基準算法；因此，α的大小取決于基準算法的選擇。由于基準算法的不同，每一對進行比較的算法會產(chǎn)生2個α。

1.3.3Google Earth實際效果對比

使用不同算法得到的研究區(qū)河網(wǎng)圖像在保證河道連續(xù)，沒有異常匯流現(xiàn)象的同時還要考慮其結果與實際地形的吻合程度。Google Earth作為一款虛擬地球軟件，可以清晰地看到全球各地的衛(wèi)星照片。通過ArcGIS將柵格中的最大值除以100作為閾值，提取大于閾值的柵格；再將值為1的點過濾，結果即為河道；之后再將河道數(shù)據(jù)轉換為用Google Earth可以打開的格式，便可將河網(wǎng)置入衛(wèi)星照片中以比較其與實際地形的匹配程度。

1.3.4運行時間對比

不同算法的實際應用價值，其運行時間也是評價的一個重要方面。因此，本文對每種算法在填洼(r.watershed算法中不包含填洼步驟)、流向計算、匯流計算3個步驟的用時也做了統(tǒng)計和比較，以便綜合實際效果等方面對每種算法做一個應用評價。為了減小誤差，采用運行3次求取平均值的方法。

2 結果與分析

2.1 運行成果對比

D8、D8-LTD、rwflood、MFD、Dinf和MDinf算法在3個研究區(qū)的DEM上計算得到的TCA圖(見表3)，TCA值越大，灰度越暗。由對比分析可看出，同一研究區(qū)，單流向算法較多流向算法水流路徑更加細密，D8和D8-LTD兩種算法最為細致且差異較小，多流向算法的水流路徑則比較寬。之所以會有這樣的差異是因為兩類算法在向下游分配水流時的方式不同，多流向算法向下游多個甚至所有的柵格進行分配，使區(qū)域內(nèi)多數(shù)柵格獲得的流量較單流向算法多，集水面積的差別減小。對于多流向算法，由于Dinf和MDinf算法分配給下游最陡方向的流量比MFD多，故灰度圖顏色相對較淡。另外，不同算法之間提取主要水系差異性較小，但在各柵格流量數(shù)值上存在較大差異，可見在應用過程中，流向算法對于注重整體研究區(qū)特征的研究影響不大，對于注重各柵格TCA數(shù)值的研究則存在較大影響。

2.2 累積頻率分布分析

圖2為3個研究區(qū)的TCA累積頻率分布圖。在ln(A)≤10.5時，6條曲線逐漸趨于重合，當ln(A)≥10.5時，6條曲線幾乎完全重合。

圖2 3個研究區(qū)ln(A)累積頻率分布

在圖2中，單流向算法(D8、D8-LTD、rwflood)與多流向算法(MFD、Dinf、MDinf)的圖線存在比較明顯的差別，易于區(qū)分。單流向算法中匯水面積較小的柵格所占比例遠大于多流向算法。這與2.1節(jié)相同，同樣是由單流向算法和多流向算法流量分配策略的差異導致的。

MFD在相同頻率下的匯水面積最大。這與其對流量的分配方式有關。所有的下游柵格都可按比例分配到一定的上游水量；從而使得匯水面積較其他算法更大。Dinf和MDinf在圖線的走向上基本相同，數(shù)值介于單流向算法與MFD算法之間。另外，從圖2中也可看出奪底溝研究區(qū)6種算法之間的差異略大于另外兩研究區(qū)，這是由于在實際地形中，奪底溝研究區(qū)下游存在部分平坦地區(qū)，多流向算法使水流發(fā)散的特征被顯著增強。

2.3 相對差系數(shù)分析

為了分析不同地形單元[17]上6種算法所計算的TCA之間存在的差異，將地形分為坡面區(qū)域(ln(A)≤10.5)和匯流區(qū)域(ln(A)>10.5)兩部分，分別對全研究區(qū)、坡面區(qū)域、匯流區(qū)域計算相對差系數(shù)(見表4)。表3中3個數(shù)據(jù)的含義分別為全研究區(qū)計算的相對差系數(shù)，山坡區(qū)域相對差系數(shù)，匯流區(qū)域相對差系數(shù)。

從表4中可以看出，除rwflood算法外，兩種單流向算法或者兩種多流向算法進行比較時計算得到的α值較大，而單流向算法與多流向算法相互比較時得到的α值較小。由于Dinf和MDinf算法在流量分配過程中分給最陡坡度的流量極大，因此與單流向算法相似度較高。D8-LTD算法是對D8算法基于全局偏差進行校正，大多數(shù)柵格流向與D8相同，所以相似度高。另外，相互比較的兩種算法在基準算法不同的情況下得到的山坡區(qū)域α值差異不大。

Dinf算法由于實際僅計算一個流向，具備單流向算法的部分優(yōu)勢，在與單流向算法相比時α值較其他多流向算法高。rwflood算法由于流向判斷基于高程而非物理意義更合理的坡度，得到的所有相對差系數(shù)均較小，與其他算法差異大；因而，在追求精度時的流向計算中應該盡量避免使用rwflood算法。此外，MFD算法在多流向算法中α值整體較小，即與其他算法差異較大，這是由于MFD算法計算中更加強調(diào)水流彌散的影響。以上分析結論驗證了劉學軍等[17]的成果。

2.4 Google Earth實際效果對比

表5展示了將應用6種算法提取出的全研究區(qū)的河道導入進Google Earth軟件后的主要河道分布圖，以及將篩選數(shù)據(jù)的閾值適當調(diào)小后得到的局部山區(qū)的水系分布圖。單流向算法獲取的河道利用ArcGIS轉成了線，多流向算法由于可以反應河道寬度故保留了原有的柵格形式。由于不同算法提取出的河道有所重疊，故圖線不能全部顯示。

從表5中全研究區(qū)圖中可以看出，除rwflood算法存在提取的主河道過于平行、順直，與實際河道偏離較多外，其他5種算法在河谷中的差異較小，提取出的主河道基本重合且與實際河道吻合程度較高，3種多流向算法反映的河道寬度較準確。在奪底溝研究區(qū)存在平坦地區(qū)，各種算法的結果存在比較大的差異，且準確性不高，這是由于平坦地區(qū)受填洼算法選取、人類活動等影響，河道的還原難度大。對于局部山坡，單流向算法得到的水系存在較多的過直河道，rwflood算法的過直河道最多，與其他算法差異最大，甚至出現(xiàn)了多處河道翻越山脊的錯誤情形。

表4 各研究區(qū)6種算法相對差系數(shù)

表5 Google Earth全研究區(qū)及山區(qū)局部導出圖

整體而言，rwflood算法對河道的還原效果最差，多流向算法反映真實河道寬度、走向的效果優(yōu)于單流向算法，D8-LTD算法較D8算法更接近于真實河道。相同閾值下，單流向算法得到的河長大于多流向算法，MFD算法因水流發(fā)散程度最嚴重得到的河長最短。在TCA值勉強達到閾值的源頭柵格，一旦出現(xiàn)水流發(fā)散的情況，存在該柵格的所有下游柵格TCA值均低于閾值的可能，因此多流向算法提取的可視化河道在水系上游存在斷流現(xiàn)象。

2.5 運行時間對比

本文所有的程序使用Java語言編輯，由JDK 1.8.0編譯，運行電腦使用Intel i5- 6200U處理器，運行內(nèi)存4GB，每種算法的運行時間見表6。

從表6可見，D8-LTD因流向校正步驟復雜，耗時最長；另外兩種單流向算法整體運行時間較多流向算法短；rwflood算法因省去填洼步驟，總運行時間最短。運行時間的差異與算法的復雜程度相符。

3 結 論

本文選擇了3個地勢起伏較大的山丘區(qū)，使用了6種算法提取流向進行相似性、真實效果、運行時間3個方面的對比。主要結論如下：

(1)不同流向算法間存在一定的差異，但是就單流向算法和多流向算法兩個大類而言，同類算法間的相似性普遍高于不同類算法，D8和D8-LTD兩種算法的相似性極高。

(2)單流向算法存在較多的平行河道，多流向算法較單流向算法可更好地模擬出河道的彎曲和寬度，把水流分配給多個下游單元的思想更加符合水流的實際情況。在河道提取方面，驗證了D8-LTD算法得到的河道由于D8算法的事實。

(3)rwflood算法與其他算法接近無關，對真實河道還原效果最差。

表6 流向算法運行時間 ms

(4)運行時間與實際效果不可兼得。rwflood雖然河道還原效果最差，但是擁有最快的整體計算速度；而準確性占優(yōu)的D8-LTD與多流向算法，在面積增大時運行時間大幅增長。在應用過程應對運行時間與實際效果進行綜合考量。