基于梯田信息的地形湿度指数表达研究

张 鹏, 赵牡丹, 郗家琪, 吴宇鑫

(西北大学 城市与环境学院, 陕西 西安 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.