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

(1.河海大學(xué)水文水資源學(xué)院，江蘇南京210098； 2.西藏自治區(qū)水文水資源勘測(cè)局，西藏拉薩850000； 3.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇南京210098)

0 引 言

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

表2 流向算法原理

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

1 材料與方法

1.1 研究材料

為了更好地研究流向算法在不同山丘區(qū)的應(yīng)用效果，本文分別在青藏高原、東南丘陵、黃土高原選取了姜灣研究區(qū)、奪底溝研究區(qū)和杜家灣研究區(qū)(見(jiàn)圖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ù)通過(guò)航天飛機(jī)搭載雷達(dá)進(jìn)行測(cè)繪， V4.1版本采用了新型的插值方法，讓數(shù)據(jù)的準(zhǔn)確度有了較大的提升。本文使用的全部數(shù)據(jù)下載自美國(guó)地質(zhì)調(diào)查局網(wǎng)站(http://earthexplorer. usgs.gov/)。

1.2 流向算法

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

1.3 研究方法

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

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

1.3.1累計(jì)頻率分布

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

1.3.2相對(duì)差系數(shù)

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

(1)

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

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

1.3.3Google Earth實(shí)際效果對(duì)比

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

1.3.4運(yùn)行時(shí)間對(duì)比

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

2 結(jié)果與分析

2.1 運(yùn)行成果對(duì)比

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

2.2 累積頻率分布分析

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

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

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

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

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

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

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

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

2.4 Google Earth實(shí)際效果對(duì)比

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

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

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

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

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

2.5 運(yùn)行時(shí)間對(duì)比

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

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

3 結(jié) 論

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

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

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

(3)rwflood算法與其他算法接近無(wú)關(guān)，對(duì)真實(shí)河道還原效果最差。

表6 流向算法運(yùn)行時(shí)間 ms

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