基于不同流向算法的地形濕度指數(shù)模型差異性分析

唐瀑漾 羅明良,2 白雷超 劉晉紅

（1.西華師范大學地理科學學院,四川 南充 637009；2.西華師范大學嘉陵江流域研究所,四川 南充 637009）

土壤流失極難恢復，土壤侵蝕是制約地區(qū)可持續(xù)發(fā)展最嚴重的生態(tài)環(huán)境問題之一[1]。溝道侵蝕的特征、發(fā)育過程及形成機理研究一直是土壤侵蝕與溝蝕地貌學領域的熱點[2]。目前基于數(shù)字高程模型（DEM）的應用越來越豐富，針對溝道識別算法模型的探究也越來越深入，相關研究主要集中于復合地形指數(shù)（CTI）、坡面面積（SA）、地形濕度指數(shù)（TWI）和坡面面積冪（SAP）4種地形指數(shù)模型的優(yōu)化和運用[3]。有研究表明，TWI在實際運用中受到流向算法和數(shù)字高程模型（DEM）分辨率的影響較大[4-5]。然而，已有文獻對于TWI模型的運用，少有同時考慮到其受不同流向算法和分辨率影響。因此，本研究以元謀干熱河谷為研究區(qū)，基于高精度DEM為實驗樣本，運用4種流向算法（Rho8、D8、Dinf和MFD-md）定量化研究不同算法之間的差異性，進一步揭示不同流向算法的適用性。

1 研究區(qū)概況

研究區(qū)域位于云南省楚雄彝族自治州北部的元謀縣黃瓜園鎮(zhèn)小雷宰村禿禿脊梁子，經(jīng)緯度范圍在101°35'E～102°069'E，25°23'N～26°06'N，海拔980～1 400 m。該區(qū)域內(nèi)土壤類型以燥紅土和變性土為主，并且土壤與巖石膠結(jié)性差，易受侵蝕[6-7]。實驗樣區(qū)流域面積1 875.39 m2，在該區(qū)域內(nèi)部發(fā)育著由溝蝕作用形成的大量侵蝕溝，本文選取其中一條溝道（見圖1），該溝道水沙輸送情況明顯，坡面存在大量細溝。

圖1 實驗樣區(qū)圖

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)采集及處理

數(shù)據(jù)主要來源于對元謀干熱河谷的典型溝道區(qū)域的野外踏勘，基于D200飛馬無人機進行激光點云數(shù)據(jù)的實測收集 。無人機的模塊選擇為D-LiDAR300模塊（測量精度為5 cm）、D-OP320傾斜模塊，同時用實時動態(tài)差分儀（real-time kinematic,RTK）進行校正。通過對0.05 m高分辨率的DEM數(shù)據(jù)處理，用CUBIC三次卷積插值重采樣方法得到0.05 m、0.1 m、0.2 m分辨率的三組DEM數(shù)據(jù)。利用ArcGIS和SAGA軟件對數(shù)據(jù)進行運算，最終得到不同流向算法在不同分辨率下的TWI示意圖。

2.2 地形濕度指數(shù)模型

地形濕度指數(shù)（topographic wetness index,TWI）是基于DEM對徑流路徑長度、產(chǎn)流面積等的定量描述，能夠用來方便地模擬理想狀態(tài)下流域內(nèi)每點土壤水分的干濕狀況的模型[8]。其應用最為廣泛的計算公式如下：

2.3 流向算法

流向算法可分為單流向算法（single flow direction algorithm,SFD）和多流向算法（multiple flow direction algorithm,MFD）[9]，其中單流向算法認為每一格網(wǎng)單元本身產(chǎn)生的流量及其上游流量都會流向其下游唯一的相鄰單元；多流向算法認為水流分布具有分散性質(zhì)，每一格網(wǎng)單元的流量并不是全部流入某一下游格網(wǎng)，而是流入比其低的所有或部分下游格網(wǎng)單元[10]。本研究中選擇了單流向算法中常見的Rho8、D8，以及多流向算法中的Dinf、MFD-md，如表1所示。

表1 本文選用的流向算法及其原理

3 結(jié)果與分析

3.1 不同分辨率DEM下的TWI

本研究針對不同分辨率的數(shù)據(jù)進行處理，利用4種流向算法取得禿禿脊侵蝕溝樣區(qū)的TWI的分布圖，如圖2所示。

圖2 禿禿脊侵蝕溝樣區(qū)TWI細節(jié)及整體分布圖

從圖2可以看出，由單流向算法到多流向算法，TWI分布圖展示趨勢由線到面，由細小溝變?yōu)榧殰蠀R水趨勢面。在相同分辨率下，單流向算法所得TWI分布圖中水流路徑細密，其中D8算法出現(xiàn)較多平行水流路徑，Rho8算法下水流路徑散亂無章；多流向算法所得TWI分布圖呈平滑的趨勢面，可以展示整體匯水趨勢情況，且簡化了部分局部信息，表達溝道時有一定主次。

分辨率差異使得溝道的表達效果不同。利用高精度的DEM計算得到的TWI分布圖可以表現(xiàn)出地表細微變化，并且TWI細節(jié)圖可以大致展示有植株區(qū)域（見圖3）。對于完整的較長的百米侵蝕溝而言，分辨率下降之后TWI分布圖能更明顯地表現(xiàn)出溝的走勢（如圖2中0.2 m DEM整體圖）。

圖3 TWI對植株的展示（以0.05 m DEM下Rho8算法為例）

3.2 TWI結(jié)果的分布特征

從表2可以看出，同一流向算法下的TWI平均值隨分辨率的減小而增大。柵格變大的同時，算法預測分配到下游柵格的水量變大，形成徑流的可能性也變大。相同分辨率下不同流向算法的極差接近，但分辨率下降極差明顯下降。相同算法下，分辨率越高，標準差越大，這說明在高分辨率DEM下的TWI值分布更分散，分辨率會顯著影響TWI值的離散程度。

表2 禿禿脊侵蝕溝樣區(qū)TWI分布特征

單流向算法中D8算法的標準差小于Rho8算法，多流向算法中Dinf算法的標準差小于MFD-md算法，這是因為Rho8和MFD-md算法隨機性更強，造成更明顯的離散性。4種算法均為過度峰度、正偏態(tài)，峰度、偏度均隨著分辨率的降低而增大。總體來說，不同算法下TWI值的空間統(tǒng)計特征差異小，不同分辨率下TWI值的空間統(tǒng)計特征差異更顯著。

3.3 TWI的變化趨勢

選擇Akima樣條曲線分析侵蝕溝的TWI頻率和累積頻率，得到了圖4。不同分辨率下的頻率變化趨勢是大致相同的。累積頻率在TWI值為1左右的時候出現(xiàn)趨于平緩的轉(zhuǎn)折，頻率在TWI值為0的時候出現(xiàn)下降轉(zhuǎn)折，說明大于0值的柵格隨著TWI值的增大，柵格數(shù)量減少，這符合匯流能力強的區(qū)域占侵蝕溝區(qū)域小部分的規(guī)律。通過分析TWI值在不同地形單元上的差異，按其分布規(guī)律可以將地形分成2個部分：非溝道區(qū)域和溝道區(qū)域。通過目視解譯圖2可以看出TWI的高值主要集中在溝道區(qū)域，而TWI的低值都出現(xiàn)在非溝道區(qū)域。結(jié)合圖3對植被的表達和圖4的TWI頻率和累積頻率可以考慮在TWI值大于0的情況下，該區(qū)域開始有匯流能力。這為利用TWI識別溝道的閾值界定提供了參考。

圖4 不同分辨率下4種徑流算法的沖溝TWI的分布頻率圖

4 結(jié)論

分辨率變化引起的差異在于：分辨率下降，柵格變大，算法預測出分配到下游柵格的水量變大，形成徑流的可能性變大。這說明，不同大小等級的溝道在進行TWI計算的時候需要考慮其尺度效應。

算法變化引起的差異在于：（1）單流向算法中D8算法會得到較多的平行水流路徑；Rho8算法所得水流路徑散亂。（2）多流向算法都簡化了部分局部信息，相較于單流向算法有一定的主次層序，能更好展示出研究區(qū)整體匯水趨勢情況。這說明，識別坡面細溝時選單流向算法進行表達，識別切溝、沖溝等溝道走向時選多流向算法進行表達更加合理。

本研究僅選擇了元謀的一條侵蝕溝樣區(qū)，實驗范圍有限，需要進一步擴大實驗樣區(qū)以找到普適性規(guī)律。雖然TWI模型可以利用地形對土壤濕度空間分布的控制來提供關于侵蝕溝潛在位置的信息，從而能夠識別溝道，但實際侵蝕溝的存在可能受到其他因素（土壤類型、土地利用和管理以及實際降雨特征等）的影響。例如，有植被的區(qū)域可以更好地維持土壤濕度，增加入滲，減少徑流和土壤侵蝕潛力，所以有植被的區(qū)域?qū)δＰ彤a(chǎn)生影響。如果能將這些因素結(jié)合到TWI模型進行計算，該模型評估潛在溝道的準確性將會得到更大的提升，從而為相關保護規(guī)劃和管理工作提供幫助。