梯田DEM构建方法精度比较研究

张 鹏,赵牡丹,王翊人

(西北大学,城市与环境学院,西安 710100)

梯田是黄土高原地区最为重要的水土保持措施之一，不仅能够极大程度地减缓土壤侵蚀过程，而且具有很好的保水保肥作用，可提高坡耕地的土壤农业生长潜力[1-3]。梯田在黄土高原地区分布广泛，包括水平梯田和坡式梯田大约400万hm2[4]，其修筑过程一定程度上改变了自然地表微形态，使得坡面坡度变缓、坡长被截断、水流路径发生了变化，形成了突变地形。目前而言，DEM一般是以现有地形图通过内插方法生成的，而中低分辨率DEM大多局限于反映光滑连续的自然地表，对自然陡坎、梯田等突变地形存在描述失真，因而促成了顾及梯田地形DEM构建方法的产生。杨蕾等[5]从遥感影像上对梯田田坎进行目视解译，通过加权取平均值法对田坎的高程进行赋值，实现了水平梯田的三维表达。周春寅等[6]利用高程增量法、偏移线法构建TIN，实现了坡式梯田的三维可视化表达。赵卫东等[7]在1∶10 000格网DEM的基础上提取了梯田的特征点、特征线以及周围边界，基于此构建了充分考虑梯田信息的Grid-TIN混合格网DEM，可以实现对梯田形态的准确描述；祝士杰，汤国安等[8]提出了梯田DEM快速构建方法，利用梯田约束特征线构建出了能够表达梯田地形的DEM；罗仪宁[9]基于野外数字测图和地形图数据获取了突变地形特征并利用ANUDEM软件将其表达到DEM中生成Hc-DEM，以实现突变地形的表达；为了能够使DEM更真实地表达梯田地形，刘芬[10]提出了基于真实田坎信息的梯田DEM构建方法，程德强等[11]在其基础上又提出了基于梯田田面的梯田DEM构建方法，实现了对梯田信息的强化表达。这些均为基于DEM进行梯田信息的表达和分析研究奠定了基础。

虽然目前在梯田DEM的构建上已经有了一定的研究成果，已能初步实现梯田DEM的构建及对梯田地形信息的强化，但不同方法构建结果的精度差异分析并没有受到太多重视，不同构建方法在实际应用中的选择上缺乏一定依据性，而这些直接影响着相关梯田DEM研究的结果精度。因此，本研究提出对不同梯田DEM构建方法进行精度评价，分析探讨不同梯田DEM构建方法的优势与不足，为顾及梯田等突变地形的DEM表达研究提供参考，同时可促进对梯田DEM构建方法的完善。

1 数据源与研究方法

1.1 研究区与数据基础

选取延安市南部燕沟流域内的某一处梯田为研究样区。燕沟流域是黄土高原中部丘陵沟壑区第Ⅱ副区下的一个子流域，属于半湿润半干旱气候过渡带，地处北纬36°28′—36°32′、东经109°20′—109°35′。燕沟平均海拔为1 198.7 m，年均降雨量为558.4 mm，流域内沟壑纵横、坡陡沟深，土壤侵蚀量大，属于强度水土流失区域。选用的基础数据包括：墨卡托投影坐标系下的研究区0.5 m分辨率Wordview-3遥感影像、5 m分辨率DEM数据以及基于激光点云数据构建的1 m分辨率DEM数据。其中激光点云数据是利用瑞士徕卡HDS8800三维激光扫描仪扫描梯田样区获取的，点云原始采样间隔为12 mm，对其进行点云去燥、平滑滤波及重采样后构TIN再生成格网DEM。5 m分辨率DEM数据是利用遥感摄影测量方式制作生成的。

1.2 研究方法

本研究以燕沟流域梯田样区为研究对象，在5 m分辨率DEM、高分遥感影像的数据支撑下，基于真实田坎信息法[10]和快速构建法[8]对样区构建顾及梯田地形的1 m分辨率DEM。以激光点云数据构建的1 m分辨率DEM作为参考，从高程、坡度、坡长等地形特征信息分析、测量误差及地形描述误差[12-14]指标评价三方面入手，对不同方法构建的梯田DEM进行地形特征分析及精度评价，探究不同方法各自的优势及不足之处。

研究选择当前最被广泛应用的方法之二：快速构建法、基于真实田坎法来构建梯田DEM以实现DEM表面对梯田信息的表达。由于样区实际地形为坡式梯田，因此本研究构建的梯田DEM以坡式梯田为前提，其断面示意图(数学模型)如图1所示。

图1 坡式梯田断面示意图

图1中α表示梯田田面倾角；β表示田坎坡度；L表示梯田田面水平投影宽度；H表示上下两田面之间的高程差；h表示田坎偏移线与田坎的高差(田坎偏移线原始高程)；b和d分别表示田坎偏移线相对田坎在水平方向与垂直方向偏移的距离和高差。由图1可以得知：

d=H-Ltanα

(1)

b=cotβ(H-Ltanα)

(2)

(3)

快速构建法根据梯田特征线提取方法的不同，分为参数构建法和简易构建法。参数构建法具体步骤有：取梯田范围内最低的田坎高程为E，以E为初始高程，田坎高度H为等高距，采用DEM自动提取等高线方法提取田坎线；再以(E+h)(公式2—3)为田坎偏移线的初始高程，H为等高距，利用提取等高线法提取台沿线偏移线；将田坎偏移线高程降低修正值Δh，其中Δh=h-Ltanα。α的大小决定梯田类型，α>0为坡式梯田，α=0为水平梯田，α<0为反坡梯田。以两组约束线高程为内插属性构TIN，转化为格网DEM。简易构建法适用于一些无法获取全部构建参数值的情况中，通过科学估计田坎高度H，再通过提取等高线法获得梯田特征线，实现快速构建梯田DEM。

刘芬[10]利用高分影像和5 m分辨率DEM获取梯田特征线的位置与高程信息，并基于此构建能反映真实梯田信息的1 m分辨率DEM。基于真实田坎法构建梯田DEM过程主要有四步：通过高分辨率遥感影像绘制每块梯田田面所对应的田坎线(台沿线)；在DEM基础上通过掩膜的方式获取每块田面的高程平均值，并将其看作为对应田坎线的高程值；根据试验样区地形情况确定梯田基本参数(α和β)，利用田坎线与田坎偏移线的数学关系对田坎线水平偏移距离b得到田坎偏移线，并高程做差d得到田坎偏移线高程值；最后，利用多组田坎线与田坎偏移线通过构TIN的方式实现梯田DEM的构建。

经过对研究区5 m分辨率DEM原始坡度进行分析、对真实梯田样区细致考量以及借鉴前人研究成果[15-16]后，选取的梯田构建参数为：田面坡度α=5°(梯田区域平均坡度)、田坎坡度β=70°。基于各参数构建样区1 m分辨率的梯田DEM，利用ArcGIS，LS_Tool[17]软件对不同方法构建出的梯田DEM进行坡度(三阶差分法)、坡长(D8算法)提取，结果如图2所示。

图2 基于不同方法构建的DEM及其地形因子提取

2 结果与分析

2.1 地形特征分析

地形属性指标是表达与研究地表形态最为有效的参数指标，也是DEM精度高低最直观的表征。由图2可知，激光点云数据构建的DEM梯田特征明显，梯田田面和田坎信息表达清晰，其他两种方法构建出的梯田DEM同样具备完整的梯田信息。因快速构建法和基于真实田坎法均是通过TIN-Grid方式构建梯田DEM的，所以表达出的地形因子不如点云数据构建出的1 m分辨率DEM自然圆滑，存在一些三角面状区域，在坡长信息的表达中较为明显。对图2中地形属性指标的特征值进行统计，结果如表1所示(表中将基于点云数据构建出的梯田DEM、基于真实田坎法构建出的梯田DEM以及快速构建法构建出的梯田DEM分别用H-DEM，T-DEM，K-DEM表示，下文亦同)。结合表1来讲，快速构建法表达出的梯田DEM对高程信息和坡度信息摆动较大，而基于真实田坎法的构建结果与1 m DEM相差较小，更接近真实值；从坡长信息特征值的对比中可以明显看出，除最大值外快速构建法表达出的坡长特征更接近1 m DEM。

表1 不同DEM地形因子特征分布值统计

图3为基于不同梯田DEM对其高程、坡度和坡长特征分布值进行统计分析的结果。高程剖面图和频率曲线图显示出T-DEM高程信息的分布与H-DEM紧密相随，其高程值摆幅也较小，各高程值段分布均匀，而K-DEM的高程信息摆动相对较大，对梯田高程值整体有抬升作用，同时对田坎和田面位置信息的表达存在微小偏差，这是由于快速构建法的梯田特征线是沿地形等高线提取的且其田面宽度和田坎高的取值均是统一值，而真实的梯田是综合考虑诸多因素建造而成的，各个田面的田面宽度和田坎高度也不尽相同，由此造成了其梯田特征线、面与真实地表的差异；在坡度信息的频率曲线图中，K-DEM和T-DEM的频率曲线摆动较大，没有H-DEM平滑自然，尤其是快速构建法。在坡度为1°～7°处存在一个较宽的波峰，其频率值远高于其他值处的频率值，说明梯田大部分区域(梯田田面)的坡度值处于这一范围内，同时3种DEM在1°～7°处的波峰并没有完全重合，与H-DEM相比，其他两者波峰均向右发生了偏移，在坡度中值区域，K-DEM和T-DEM的频率值均小于H-DEM，在两端极值区域，两者频率值又大于H-DEM，说明两种构建方法对坡度值的表达均存在极端化趋势。总的来看，T-DEM对坡度信息的表达更接近H-DEM；从坡长(自然对数)频率曲线图可以看出梯田地形主要以短坡为主，超过一半区域坡长值小于4.5 m，且短坡坡长值都比较集中。K-DEM在坡长信息的表达中整体与H-DEM更为接近，但是在短坡长值处的频率(坡长为0.5 m)存在过于偏大的现象。T-DEM提取出的坡长(自然对数)曲线除在短坡区域更接近H-DEM外，长坡值(坡长大于10 m)频率要高于其他两者，这也是其坡长平均值和标准差值大于其他两者的原因，对应地坡长累计频率曲线图中T-DEM坡长值到36 m时累计频率曲线才达到稳定，远大于其他两种DEM对应的坡长稳定值。

2.2 测量误差分析

DEM数据在不考虑粗差和系统误差时，其数据误差基本由随机误差构成，通常以有限的采样点为样本，利用中误差(RMSE)来衡量DEM数据精度。为避免中误差作为唯一指标的绝对性，研究除中误差外，另选标准差(SD)和平均绝对误差(MAE)作为DEM数据误差的评定指标，其计算公式如表2所示(设H-DEM的特征值为Z，另外两种DEM数据的特征值为z，则误差ε=Z-z，n为样本点个数)。在梯田样区内选取200个采样点，分别提取不同DEM及其地形特征因子在每个采样点处的特征值，并以点云高精度数据构建的DEM为参考，统计不同方法下3种地形特征因子的误差分布情况，见表3。

与H-DEM相比，两种构建方法对3种地形特征因子均存在一定的表达误差。对于高程和坡度因子，基于真实田坎法的3种误差指标均小于快速构建法，说明基于真实田坎法表达出的高程和坡度信息的稳定性、整体误差大小及误差极值情况均高于快速构建法，因此在这两个因子的表达上基于真实田坎法要优于快速构建法。而从坡长因子的误差描述中可看出快速构建法对坡长因子的表达更准确一些，与上文地形特征分析中的结果正好吻合。

2.3 地形描述误差分析

即使DEM高程点上的误差均为零时，有限的栅格采样点也只能实现对真实地表的近似模拟，汤国安将这种情况下构建出的模拟地形与真实地形之间的差异称为地形描述误差[14]，并基于此提出了DEM地形描述精度数学模型来量化地形描述误差。通过栅格窗口分析法即可求取地形描述误差(Et)，当分析窗口为3×3时，DEM中某一栅格单元(第i行，第j列)的Et可表示为：

Et(i,j)=H(i,j)-(H(i-1,j-1)+H(i-1,j+(1)+H(i+1,j-1)+H(i+1,j+(1) )/4

(1)

图3 不同地形特征因子分布值统计

名称中误差(RMSE)标准差(SD)平均误差(ME)公式RMSE=∑ni=1εi2nSD=∑ni=1εi-ME()2nME=∑ni=1εin

表3 不同构建方法误差分布统计

由地形描述误差数学模型可看出，当DEM分辨率为d时，提取出的地形描述误差分辨率为2 d。在ArcGIS软件中利用Python编码可实现3种梯田DEM的Et值提取，表4为不同梯田DEM地形描述误差的特征分布值统计结果。

表4 不同梯田DEM的Et特征分布值统计

从表4可以看出：H-DEM的收敛性最高，数据波动最小，也就意味着最接近地表真实值。K-DEM的Et值的平均值和标准差值最大，分别为0.009 1，0.387，说明其Et值变化起伏较大；T-DEM的Et值波动范围最广，其最大值和最小值分别为5.08和-4.16，但是平均值和标准差的数值表明其相对更稳定，与H-DEM更为接近。总之，从地形描述误差表达出的信息来看，基于真实田坎法能够在整体上更好地描述实际地表，而快速构建法在最大、最小误差上与H-DEM相近，但在整体的地形描述精度上要稍弱于基于真实田坎法。

3 结 论

不同方法构建的梯田DEM均能表达出梯田地形特征，梯田田面和田坎特征清晰，但是不同构建方法各有优势与不足：(1) 从测量误差和地形描述误差的误差指标统计和对比分析结果来看，相比快速构建法，基于真实田坎法的构建结果在各种地形因子的表达上整体误差更小，更接近梯田真实地表形态。(2) 从地形特征分析结果可知，快速构建法构建出的梯田DEM在梯田田面和田坎位置的描述上与真实梯田地表存在偏差，这是其构建方法本身决定的。基于真实田坎法则整体精度较高，对各种地形因子的描述更准确。(3) 快速构建法虽构建精度较低，但具有快速简单、省时省力的优势，在大面积的梯田DEM构建中具有优势；基于真实田坎法构建出的梯田DEM能够精确地描述原始地表形态，适用于精度较高的梯田地形研究中，但在遥感影像中梯田特征线的自动提取上存在改善的空间。同时梯田特征线、面要素(田面数、田坎线位置)和梯田参数(田面倾角、田坎坡度)是其构建精度高低的决定性因素，须谨慎确定；三维激光扫描数据对田面田坎信息的描述准确度很高，虽然目前仅限于小范围区域且其数据成本较高，但不可否认的是，三维激光数据未来必然会成为最为常用的地形数据之一。