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

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

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

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

1 研究區(qū)概況與數(shù)據(jù)源

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

2 研究方法

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

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

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

圖1坡式梯田斷面示意圖

由圖1可以得知：

d=H/Ltanα

(1)

b=cotβ(H-Ttanα)

(2)

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

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

3 結(jié)果與分析

3.1 坡度因子分析

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

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

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

3.2 單位匯水面積因子分析

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

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

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

3.3 地形濕度指數(shù)分析

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

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

表1 地形濕度指數(shù)統(tǒng)計分布特征值

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

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

4 討論與結(jié)論

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

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

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

不同DEM數(shù)據(jù)源對于地形濕度指數(shù)等地形特征要素的準確表達所產(chǎn)生的干擾是一個值得關注的點，特別是在局部區(qū)域高精度的地形分析中，對類似梯田這種突變地形的特征信息進行表達時選擇DEM數(shù)據(jù)應該慎重。同時，在后來的梯田地形表達研究中兼顧地形特征與水文等屬性特征的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] 鄧慧平,李秀彬.地形指數(shù)的物理意義分析[J].地理科學進展,2002, 21(2):103-110.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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