基于梯田信息的地形濕度指數表達研究

張 鵬, 趙牡丹, 郗家琪, 吳宇鑫

(西北大學 城市與環境學院, 陜西 西安 710100)

地形濕度指數(topographic wetness index, TWI)是1979年由Beven和Kirkby提出的用于反映土壤飽和缺水量空間分布的參數[1]，是一種基于數字高程模型(DEM)的對流路徑長度、產流面積以及土壤徑流產生能力的重要量化指標[2-3]，廣泛應用于土壤、水文、地貌等領域[4-6]。它通常表達為單位匯水面積因子SCA(specific catchment area)與坡度因子β的正切值比值的復合函數，即ln(SCA/tanβ)[7]。在黃土高原地區，水土流失和干旱缺水問題一直是制約其生態社會發展的重要因素之一，關于該地區土壤含水量的研究歷來是一個比較受關注的話題[8-10]。梯田是黃土高原地區耕地最為主要的組成部分之一，也是最為重要的水土保持措施之一，梯田在改善土壤侵蝕的同時，也造成了地表自然形態的顯著變化，這對于梯田范圍內的區域地形特征、水文特征等產生了顯著的影響[11-14]。因此，分析梯田微地形對小區域范圍內的地形濕度指數表達及其空間分布狀態的影響呈現一種怎樣的形態是一個值得探究的內容。另一方面來講，基于梯田的地形濕度指數表達研究對于黃土高原地區水土保持、徑流路徑分析等研究有著一定的指導作用。

1 研究區概況與數據源

選取的研究區域位于延安市南部3 km處的燕溝流域，該流域是黃土高原中部丘陵溝壑區第Ⅱ副區下的一個子流域，屬半濕潤半干旱氣候過渡帶。地處北緯36°28′—36°32′，東經109°20′—109°35′，隸屬延河支流，流域總面積約為47 km2，流向為東南—西北走向。燕溝年均氣溫為9.8 ℃，年均降雨量為558.4 mm，6—9月為集中多雨季節，占全年70%以上，降雨年際變化較大。燕溝流域中的梯田面積約為總流域面積的8%左右，土壤以黃綿土為主，流域植被覆蓋率較低，土壤侵蝕量大，屬于強度水土流失區域。本研究以燕溝流域內某處梯田為試驗樣區，其面積約為0.8 km2。選用的基礎數據包括：燕溝流域0.5 m Wordview-3遙感影像、5 m分辨率DEM數據以及利用瑞士徠卡HDS 8800三維激光掃描儀掃描研究區獲取的激光點云數據，其原始點云分辨率為12 mm/100 m。

2 研究方法

以梯田區域為研究對象，基于真實田坎法構建研究區1 m分辨率的梯田DEM；同時以激光點云數據為基礎，對其點云去燥、平滑濾波以及重采樣后通過插值生成高精度1 m DEM來高保真反映真實梯田地形。從坡度、單位匯水面積兩方面對5 m DEM數據、基于真實田坎法構建出的梯田DEM、點云1 m DEM數據進行對比分析，探索梯田信息的缺失與否對地形濕度指數表達的定量影響，同時以點云數據插值生成的1 m DEM為參考，探究基于真實田坎法構建出的梯田DEM與真實梯田地表的差異及其差別在地形濕度指數表達上的體現。

對于梯田DEM的構建，當前被廣泛應用的方法大致可歸為3類：基于快速構建法[15]、基于真實田坎法以及基于野外實測數據的梯田DEM構建方法[16]。快速構建法簡單便捷，但構建出的梯田地形與真實地表存在較大偏差；野外實測數據構建出的DEM精度很高，但對數據要求較高；而基于真實田坎的構建方法在實現對梯田地表真實模擬的同時，不存在太多數據要求。因此，本文選擇基于真實田坎法構建梯田DEM來實現對梯田信息的表達，由于樣區實際地形為坡式梯田，因此構建梯田DEM以坡式梯田為前提，坡式梯田斷面示意圖(數學模型)如圖1所示。

注：α表示梯田田面傾角；β表示田坎坡度；L表示梯田田面水平投影寬度；H表示上下兩田面之間的高程差；b,d分別表示田坎偏移線相對田坎在水平方向與垂直方向偏移的距離和高差。

圖1坡式梯田斷面示意圖

由圖1可以得知：

d=H/Ltanα

(1)

b=cotβ(H-Ttanα)

(2)

基于真實田坎法構建梯田DEM過程主要有4步：①通過高分辨率遙感影像繪制每塊梯田田面所對應的田坎線(臺沿線)；②在DEM基礎上通過掩膜的方式獲取每塊田面的高程平均值，并將其看作為對應田坎線的高程值；③根據試驗樣區地形情況確定梯田基本參數(α和β)，利用田坎線與田坎偏移線的數學關系對田坎線水平偏移距離b得到田坎偏移線，并高程做差d得到田坎偏移線高程值；最后，利用多組田坎線與田坎偏移線通過構TIN的方式實現梯田DEM的構建。經過對真實梯田樣區進行細致考量以及借鑒前人研究成果后，選取田面坡度構建參數分別為：田面坡度α=3°，田坎坡度β=70°。圖2分別為研究區的5 m DEM、基于真實田坎法構建出的1 m梯田DEM以及基于點云數據插值構建的高精度1 m DEM。為分析方便，文中將5 m DEM稱為原始DEM，將基于真實田坎法構建出的梯田DEM稱為T-DEM，將點云數據生成的1 m DEM稱為H-DEM。

圖2 原始DEM(a)，T-DEM(b)，H-DEM(c)

3 結果與分析

3.1 坡度因子分析

在地形表面分析中，坡度(slope)是最能直接體現地形起伏和高程變化劇烈程度的地形因子之一，在土壤侵蝕、地形水流模擬分析等方面，坡度因子同樣是影響土壤抗侵蝕能力以及水流路徑的關鍵因素，在地表地形濕度指數研究中，坡度更是土壤出水能力的表征。為分析梯田對地形坡度的影響，對3種不同DEM進行坡度分析，選擇三階差分法來計算地表坡度因子[17]。

結果表明，DEM上梯田信息的表達與否對于地表坡度的影響非常明顯，原始DEM坡度只具有宏觀的坡度分布特征，全區域坡度分布范圍在0°～47°之內，大部分區域坡度波動在15°～35°之間，有極少部分區域坡度在8°以下。基于點云數據構建出的H-DEM則真實地還原了梯田地形的坡度信息，所有區域的坡度分布在0°～71°之內，田面的坡度大都分布在0°～15°范圍值之間，坡度在60°～71°之間的區域基本都是梯田田坎的位置，連續兩田面過渡處的高程落差是形成高坡度值的根源。基于真實田坎法構建的T-DEM坡度范圍分布在0°～82°之間，坡度在45°～82°內的區域基本是田坎。從整體看，T-DEM在坡度空間分布上同樣能夠清晰地反映出梯田地形下的坡度特征，只是與H-DEM相比：雖然DEM田面和田坎區分很清晰，但過渡區域稍顯生硬，紋理特征表達不如后者地形表達真實自然；構建出的DEM梯田田面過度平滑，失去了自然地形本身的非絕對平整性，而H-DEM所表達出的坡度更顯隨機自然；T-DEM在坡度表達上有兩極化的現象，將田坎位置的坡度增大化、田面坡度減緩化。坡度頻率分布(圖3)中可以更明了地體現出3種DEM提取出坡度的差別：原始DEM在坡度表達上填低削高，坡度范圍向坡度中值區5°～35°范圍匯聚；H-DEM坡度高中低值分布相對均勻，只是在0°～3°分段處頻率分布最高，這是由梯田的獨特地形特征所決定的；T-DEM同樣能夠反映出梯田特征對坡度表達的影響，但是在0°～3°范圍相較其他兩者頻率分布相當高，同時在>45°范圍也同樣略有偏高，而在中值區域頻率較低，地形表達存在兩極化現象，這是由構建梯田DEM過程中設立統一的田面坡度以及田坎坡度造成的。

圖3 不同DEM坡度分級頻率分布

3.2 單位匯水面積因子分析

單位匯水面積(specific catchment area，SCA)是指單位等高線長度上的上游匯水面積或者單位等高線上的徑流面積，描述了地表土壤的匯水能力，是各種地貌結構和水文模型如地形濕度指數、水流強度指數等的重要參數，廣泛應用于地貌結構研究、土壤水分空間分析、流域網絡分析等研究中[18]。對于單位匯水面積而言，流向(flow dirrection)是決定其分布及量化值的主要因子，因此流向算法的選擇對于單位匯水面積的計算至關重要。梯田地形表面平滑渾圓，田面坡度較為平緩，水流特征以漫散徑流為主，適宜用多流向算法來模擬其流向特征[19]，本研究選取的多流向算法為D-Infinity算法[20]，利用David Tarboton團隊合作開發的TauDEM Tools工具集可實現基于D-Infinity多流向算法對DEM單位匯水面積因子的提取。

從圖4來看，基于D-Infinity算法提取出的不同DEM單位匯水面積值差異不是特別明顯，3種數據源提取的單位匯水面積值大都集中于1～100，其中原始DEM計算出的SCA值有97.1%位于0～136，T-DEM和H-DEM分別為96.3%和95.1%，頻率分布非常接近，但在0～5的低值范圍內三者差距稍大。同時原始DEM，T-DEM和H-DEM三者單位匯水面積的平均值分別為32.05，33.81，35.31，原始DEM平均值稍低是由于單位匯水面積的尺度效應引起的。由于單位匯水面積的量算值與水流累積量矩陣密切相關，而水流累積是從上游逐步累積至下游的一個遞進式的過程，梯田區域地形較為平整且屬于上游處，因而從中提取的單位匯水面積大多都集中在0～100，在田面上由水流沖刷出來的細小溝壑的單位匯水面積值稍大一些。因此，不同DEM對于單位匯水面積因子的提取影響效應并不是很明顯。

圖4 不同DEM單位匯水面積頻率分布

3.3 地形濕度指數分析

地形濕度指數是地形坡度與單位匯水面積因子值的復合函數，利用公式TWI=ln(SCA/tanβ)即可實現對地形濕度指數的表達。需要注意的是不論從數學角度還是實際地形角度出發，tanβ均不能為0值，因此在計算中對于DEM坡度值為0的柵格需要給一個坡度增量α，本研究將其取值為0.000 000 000 1°，同樣地單位匯水面積取值最小為1個柵格單元長度(根據單位匯水面積定義確定)，在ArcGIS軟件的柵格計算器工具中實現代碼如下：

TWI=ln((Con("flowacc.tif"==0,m,"flowacc.tif"))/Tan(Con("slope-DEM"<=0.000 000 000 1,0.000 000 000 1,"slope-DEM")*3.141 592 6/180))。

其中TWI為地形濕度指數，flowacc.tif為單位匯水面積，slope-DEM為坡度，m為柵格單元長度。計算出不同DEM表達的地形濕度指數之后，對其結果進行統計如表1所示。結果表明3種不同DEM對于地形濕度指數的提取結果差異較大，T-DEM的地形濕度指數平均值明顯大于其他兩者，而且原始DEM與T-DEM不同程度地拉伸了TWI的波動范圍。

表1 地形濕度指數統計分布特征值

圖5為地形濕度指數值的頻率分布。從圖5上可以看出：3種DEM提取出的TWI都存在數值跳躍的現象，5 m DEM TWI值的跳躍更加劇烈一些，1 m DEM則最為連續。同時對比不同DEM的TWI提取結果可以發現，5 m分辨率DEM對梯田信息的模擬失真及對坡度信息的表達偏差使得梯田區域TWI值表達不連續且存在斷續的“高值面塊區域”，而且理論上在田坎坡度陡轉處TWI值應較小、梯田田面應是TWI的高值處，這些在5 m分辨率DEM上均表現不出來，因而其對土壤水分的模擬存在不準確性。1 m分辨率DEM能對梯田區域實現高保真模擬，田面處TWI值能清晰自然地反映出流水特征信息，可以很好地映證出梯田的保水保肥作用。基于真實田坎構建的梯田DEM表達出的TWI在梯田田面處均為高值區段、田坎處為低值區段，較為準確地反映出了田面田坎TWI值的分布狀況，但因構建時對田面信息的理想化，使得田面上TWI值普遍偏高，同時流水特征信息表達存在一定偏差。

圖5 不同DEM提取TWI頻率分布

4 討論與結論

本研究對顧及梯田的地形濕度指數表達進行了探討，以地形濕度指數的主要影響因子—坡度和單位匯水面積因子對不同數據源地形濕度指數的表達差異進行了深入分析比較。綜合對比3種數據源對地形濕度指數表達的差異性發現：

(1) 5 m DEM因其自身數據對梯田信息表達的缺失，不論是在坡度還是在TWI表達上都存在一定的信息偏差和缺失，基本無法區分出梯田地形中的田面和田坎特征信息，雖然在TWI表達的統計分布特征值描述上與1 m DEM差異不太明顯，但是在TWI空間分布以及細節表達上差異較大。

(2) 基于真實田坎法構建出的DEM能夠較為精確地表達出梯田的位置、形態以及細節特征，能很好地區分出田面田坎的TWI差異性。但與高精度1 m DEM相比，構建數學方法的局限性使得田面上各點的高程以及坡度變化均帶有較明顯的機械性，構建結果偏理想化，并伴有兩極化的趨勢，丟失了真實田面內部的地形變化信息，導致田面的流水特征信息表達存在偏差，因而計算出的地形濕度指數并不能很真實地反映出梯田的土壤水分分布特征。

不同DEM數據源對于地形濕度指數等地形特征要素的準確表達所產生的干擾是一個值得關注的點，特別是在局部區域高精度的地形分析中，對類似梯田這種突變地形的特征信息進行表達時選擇DEM數據應該慎重。同時，在后來的梯田地形表達研究中兼顧地形特征與水文等屬性特征的DEM表達方法的建立是一個需待探究的問題。

[1] Beven K J, Kirkby M J. A physically based, variable contributing area model of basin hydr-ology/Un modèleà base physique de zone d’appel variable de l'hydrologie du bassin versant[J]. Hydrological Sciences Bulletin, 1979,24(1):43-69.

[2] 鄧慧平,李秀彬.地形指數的物理意義分析[J].地理科學進展,2002, 21(2):103-110.

[3] 王洪明,楊勤科,姚志宏.小流域尺度土壤水分與地形濕度指數的相關性分析[J].水土保持通報,2009, 29(4):110-113.

[4] 凌峰,杜耘,肖飛,等.分布式TOPMODEL模型在清江流域降雨徑流模擬中的應用[J].長江流域資源與環境,2010,19(1):48-53.

[5] 楊琳,朱阿興,秦承志,等.基于典型點的目的性采樣設計方法及其在土壤制圖中的應用[J].地理科學進展,2010,29(3):279-286.

[6] 張彩霞,楊勤科,李銳.基于DEM的地形濕度指數及其應用研究進展[J].地理科學進展,2005,24(6):116-123.

[7] 張鍍光，王克林，陳洪松，等.基于DEM的地形指數提取方法及應用[J].長江流域資源與環境，2005,14(6):715-719.

[8] 王洪明.基于DEM和實測數據的小流域土壤水分模擬[D].西安:西北大學,2009.

[9] 王信增.延河流域土壤水分狀態及尺度效應[D].陜西 楊凌:西北農林科技大學,2012.

[10] 姚志宏, 楊勤科,王春梅,等.基于GIS的黃土丘陵區小流域土壤水分模擬[J].草地學報,2011,19(3):525-530.

[11] 馮瑤.基于梯田DEM的水流路徑模擬[D].西安:西北大學,2015.

[12] 劉芬.黃土高原梯田DEM地形特征研究[D].西安:西北大學,2015.

[13] 王翊人,趙牡丹,馮園,等.梯田對土壤侵蝕地形因子擾動特征研究[J].山東農業大學學報:自然科學版,2017,48(1):46-51.

[14] 趙衛東.顧及梯田地形的數字高程模型研究[D].南京:南京師范大學,2011.

[15] 祝士杰,湯國安,張維,等.梯田DEM快速構建方法研究[J].測繪通報,2011 (4):68-70.

[16] 李慧.梯田DEM構建方法研究[D].西安:西北大學,2014.

[17] 陳楠,王欽敏,湯國安,等.6種坡度提取算法的應用范圍分析:以在黃土丘陵溝壑區的研究為例[J].測繪信息與工程,2006,31(4):20-22.

[18] 湯國安,李發源,楊昕,等.黃土高原數字地形分析探索與實踐[M].北京:科學出版社,2015.

[19] 龔秒.基于DEM的地形濕度指數不確定性研究[D].南京:南京師范大學,2015.

[20] 鄔倫,汪大明,張毅.基于DEM的水流方向算法研究[J].中國圖象圖形學報,2006,11(7):998-1003.