王佳希, 邓家勇, 张 岩†, 宋晓鹏, 马 瑞
(1.北京林业大学水土保持学院,100083,北京;2.山西吉县森林生态系统国家野外科学观测研究站,100083,北京)
地形地貌是影响土壤侵蚀的重要因子。土壤侵蚀研究中,数字高程模型(digital elevation model, DEM)是重要的地形地貌基础数据,可以用于提取坡度、坡向及水系特征;计算侵蚀沟形态指标和土壤侵蚀模型中的地形因子;分析土壤侵蚀过程中的地貌变化[1-4]。DEM的制作早期只能通过实地勘测或者使用已有地形图[5]。随着测绘技术的不断发展[6-7],基于GPS[8]、卫星和航空遥感[9]、摄影测量[10]、三维激光扫描[11]等都可制作DEM。三维激光扫描所产生的数据是具有空间坐标及高程信息的点,称为点云[11]。虽然三维激光扫描具有一定的穿透作用,但由于植被的反射作用等,点云中仍会存在植被等非地面点。所以点云仍需经过滤波、插值[12-14]而最终生成高精度的DEM[13,15-19]。点云滤波方法有数学形态法滤波[20]、坡度滤波[21]、基于回波次数[13]、基于回波强度[22]等。常用的插值方法有反距离权重法(inverse distance weighted,IDW)、克里金法(Kriging)、自然邻域法(natural neighbor,NN)[14]和不规则三角网(triangulated irregular network,TIN)[8]等。
点云滤波在平坦、植被较单一地区相对简单,在崎岖、植被低矮复杂地区较为困难[22]。黄土高原丘陵沟壑区土壤侵蚀严重,地形破碎,植被多为次生演替的灌木和草本,属于点云滤波困难的地区。潘少奇等[23]将点云用TIN法插值生成DEM,与1∶1万的地形图比较,精度有所提升,但没有植被滤波处理。马鼎等[22]利用点云的激光回波强度衰减补偿模型进行了植被滤波,得到的DEM效果较好,但需要点云记录有回波强度信息。
如上,植被对于DEM的影响研究多集中于植被滤波、插值方法的开发[24-25],而关于植被季相变化对DEM的影响研究较少。植被非生长季时,枯萎凋落,激光脚点落在地面的可能性增大,DEM精度可能会提高;但是植被季相变化具体有何影响,缺乏数据支撑。笔者基于黄土高原丘陵沟壑区植物生长季 (2018年9月)、非生长季(2018年12月)2次实地测绘的点云,分析生长季与非生长季的点云以及DEM差异,讨论植被季相变化是否对其具有影响,探究提高DEM精度的手段。
陕西省绥德县桥沟小流域(E 110°17′24″~110°17′43″,N 37°29′33″~37°30′13″)位于无定河下游(图1),属于黄土高原丘陵沟壑区第一副区。气候类型为温带半干旱大陆性季风气候。土壤类型均为黄绵土亚类。植被类型以草地、灌木林地为主,有少量人工林地。植物群落多为次生演替,代表性群系有白羊草(Bothriochloaischaemum)群系、白草(Pennisetumflaccidum)-铁杆蒿(Artemisiasacrorum)群系、茭蒿(Artemisiagiraldii)群系、茵陈蒿(Artemisiacapillaris)群系和兴安胡枝子(Lespedezadavurica)群系。人工林造林主要有刺槐(Robiniapsdudoacacia)、油松(Pinustabulaeformis)、侧柏(Platycladusorientalis)、柠条(Caraganakorshinskii)等[26]。在桥沟小流域中选取A、B和C样区。样区A(1 445 m2)切沟沟头边缘平滑,沟壁较陡,沟道比降在1.08~1.25之间;以禾本科,菊科等草本为主,沟底有柠条,冠层高度约20~60 cm;植被覆盖度较低。样区B(3 177 m2)沟壁较陡,沟头边缘锋利,沟道比降在0.75~1.11之间,切沟较大;除垂直沟壁裸露外,在沟底、沟间有零星刺槐分布,其高度约2.5~3.5 m;灌草覆盖度较高,灌木冠层高度约30~50 cm。样区C(1 386 m2)切沟密集,沟道比降较大(1.2~2.0),但切沟宽度较小;以低矮灌木和草本为主,灌木主要分布在沟底,冠层高度20~30 cm;植被覆盖度中等。
A、B和C为三维激光测量样区编号。底图来自Google 影像。Sampling sites are named A,B and C for 3D laser measurement. Background images are from Google Earth images. 图1 研究区位置图Fig.1 Location of the study area
已有研究表明黄土高原植被覆盖在7—9月达到峰值,12月达到最低值[27]。本研究分别于植物生长季 (2018年9月)、非生长季(2018年12月),使用拓普康IS—IMAGING STATION型三维激光扫描仪对桥沟内切沟特征差异较大的3个样区(A、B和C)内11条切沟进行了扫描。为了保证测量精度,无障碍物遮挡,仪器架设点选择在样区对面山坡视野最佳且大致水平的位置,并用手持GPS确定架设点地理坐标及高程。2次测量使用了相同的架设点,保证2次测量之间相对误差较小(图1)。点云平面间隔为0.15 m,密度约30~50 pts/m2。所得点云不包括回波次数、强度等。另外,在2次外业时,使用激光测距仪和卷尺抽样测量部分优势种群的株高,用来辅助判断植被季相变化。
2.2.1 点云处理流程 第1步,通过人工判读除去包括飞鸟、电力线路、仪器偶然误差等造成的离群值[13]。第2步,基于Terrascan,利用坡度滤波法[21]将点云分类为地面点和非地面点(大部分为植被点),随后去除分类中的非地面点,仅保留地面点。第3步,将滤波前后点云用TIN法[8]、克里金法[28]分别插值构建DEM(栅格分辨率0.15 m)。按照上述1~3步分别处理9月、12月数据。最终分别得到2季相的4幅DEM:TIN法、滤波-TIN法、克里金法、滤波-克里法。
2.2.2 点云和DEM精度评价方法 对点云高程进行数理统计,计算点云平均高程、标准差、最大值和最小值等,对比2次点云整体差异。
采用中误差和判定系数表征DEM和点云的符合程度。中误差,即均方根误差(root mean square error,RMSE)。它是最常用的体现DEM偏移情况的指标,计算公式为:
(1)
式中:σ为中误差,m;R为点云高程,m;Z为插值高程,m;n为测得的总点数;m为第m个点。中误差越小,DEM越符合点云。
判定系数(R2)用以表征插值高程与实测高程总体吻合程度。计算公式为:
(2)
图2 A、B和C样区点云高程频数分布图Fig.2 Frequency distribution of point cloud elevation in the site A,B and C
本研究共选取7条切沟纵剖面,29条横剖面来提取沟长、沟宽和沟深,用来代表不同方法生成的DEM地貌形态差异。另外,通过三维可视化、坡度分析和切沟纵剖面图对比不同季相,不同滤波插值方法对生成DEM的影响。三维可视化可以直观判断DEM质量,地形地貌还原程度,植被覆盖特征等。坡度图可以直观地显示出DEM坡度变化,植被在其中多表现为“鼓包”,另外从坡度图中也很容易识别出沟沿线等地貌特征。本研究区地面粗糙度不大,而三维可视化图、坡度图中粗糙度较大区域,也一定程度上反映了植被表面的影响。
3.1.1 点云高程统计特征 点云高程统计分析表明,样区A、B和C生长季平均高程比非生长季分别高0.22、0.47和0.07 m。实地抽样方面,共测量优势种灌木、多年生草本3种,9和12月分别测184和189株。其平均株高分别为0.73和0.68 m。测量的9株乔木高度无明显变化。非生长季点云标准差普遍比生长季略低,且中位数、众数也有不同程度减小,说明非生长季高程分布更趋于集中,且整体偏低。从图2可看出3个样区类似的偏态分布规律,且非生长季比生长季的高值部分略少,其他部分明显增多。样区B点云频率分布特征差异最大。植被覆盖度方面,样区A约0.60,样区B约0.85,样区C约0.75,说明点云高程变化和植被覆盖度有一定关系,植被覆盖度越高,点云高程季相变化越大。
3.1.2 植被季相变化与滤波效果 3个样区12月所测点云高程平均标准差比9月有所降低,说明非生长季的点云高程离散趋势减弱。12月所测数据滤波后保留的点云更多(表1)。从点云剖面图可看出,9月所测点云中高大植被在12月几乎消失,这说明经过坡度滤波后,高大乔灌被去除同时地面点尽可能地被保留(图3)。样区A、B和C点云平均高程分别降低了0.25、0.56和0.44 m(9月)。
表1 4种滤波插值方法所得DEM精度统计
深色点云为滤波前,浅色点云为滤波。 The dark point cloud is the one before filtering, and the light point cloud is the one after filtering图3 滤波前后切沟横剖面点云分布图Fig.3 Distribution of point cloud of gully profiles before and after filtering
首先,9月所测点云经过滤波之后,点云数目大幅度下降。12月所测点云经过滤波之后,点云数目也有所下降,但幅度不及9月。尤其样区A,可能是由于其以草本为主,且植被盖度较低。克里金与TIN法所得点云高程均值几乎相同,但滤波之后高程有所降低:9月所测点云大概0.4 m,12月所测点云在0.2 m左右。TIN法与克里金法交叉验证的R2都在0.9以上,且克里金和TIN法中误差在0.30~0.50 m之间,滤波后插值中误差降低至0.20~0.10 m(表1)。TIN法中误差普遍比克里金法要小,说明TIN法所得高程更接近于实测点云,地形地貌刻画精度较高。
由图4可见, b比a插值结果更平滑,粗糙度降低,但是c,d粗糙度差异不大。4种方法中沟沿线都十分清晰,a保留的地貌和植被细节更多。a、b中的高大植被在c、d中已被去除。a3中2次测得剖面相比,中上部基本吻合,下部12月比9月明显更低,并且12月下部与中上部更为平滑,更符合实地特征。这是由于沟底的乔灌在12月落叶枯萎造成的。b3比a3相比,其不能去除9月植被且平滑了沟头坡度。经过滤波之后,c3,d3分别比a3、b3更为平滑,并且2次测得的结果基本一致,但是d3沟头部分仍有明显差异。
a为TIN法,b为克里金法,c为滤波-TIN法,d为滤波-克里金法;1为三维可视化叠加山体阴影效果图,2为坡度图,3为对应三维可视化叠加山体阴影效果图中的白色剖线。a refers to TIN, b refers to Kriging, c refers to filtering-TIN, and d refers to filtering-Kriging. 1 refers to 3D model with hill shade effect, 2 refers to slope map, and 3 refers to profiles of white lines in corresponding 3D model.图4 4种方法生成的DEM比较(以样区B为例)Fig.4 Comparison of DEMs by four processes (Taking site B as an sample)
在2次三维激光扫描之间,研究区内没有发生强烈的侵蚀性降雨,不存在剧烈侵蚀的情况。DEM中切沟形态参数变化主要是由植被变化引起的。切沟形态参数统计表明(表2),克里金法所得的12月平均沟长、沟宽和沟深都更大,相差0.05~0.35 m。而TIN法所得12月平均沟长和沟宽比9月分别减小0.72和0.26 m,但沟深更大。滤波-TIN法所得12月平均沟长、沟宽和沟深比9月分别增大0.03、0.25和0.10 m,说明具有去除植被的效果。由于经过滤波后,12月保留更多地貌特征点,即非生长季比生长季更接近于真实地貌。滤波-克里金法所得12月份的沟长和沟深分别比9月份小0.36和0.34 m,说明克里金法的平滑作用使得插值结果稍偏离实际地貌形态。
表2 4种处理方法提取的切沟形态参数比较
点云平均高程和标准差等都表明12月点云高程总体上低于9月,高程分布更为集中。点云高程高值部分在9月数量比12月多,而中低值9月比12月数量要少(图2)。这是由于非生长季植被枯萎凋零,覆盖度降低,植被平均高度降低,使激光脚点落在地面的概率增大。由于三维激光扫描仪无法测量到植株的最高点,所以点云中的植被点总小于最大植株高度[29],即不能精确反映植被高度,加之抽样株高不能完全反映植被总体情况,所以9月与12月点云平均高程差(0.10 m),略高于抽样平均株高差(0.05 m)。总体上看,植被变化与点云变化规律相符合,表明植被变化对于点云具有一定影响。
对滤波前后的点云和DEM进行了对比分析,结果表明坡度滤波法对于去除植被具有一定效果。在植被滤波完全去除植被的理想状态下,不同季相所得滤波结果应完全一致。但本文9月与12月植被滤波后的点云并不一致,样区A 2次点云平均高程相差在0.11 m左右,样区B和C相差了0.3 m左右。所以即使滤波之后,植被季相变化对点云高程的影响仍然存在,进而影响DEM精度。
前人比较克里金法与最邻近法、反距离权重法等将点云插值生成DEM的精度,结果表明克里金法效果较好[28]。本研究将克里金法和TIN法进行比较,发现TIN法比克里金法R2更高,σ更小,且σ均符合国家测绘局1∶500 DEM一级精度要求(坡度25°以上,栅格分辨率为0.5 m的,一级精度σ<0.70,坡度6°~25°,一级精度σ<0.50)。TIN法生成的DEM与点云拟合效果更好,精度更高,对滤波前后的地形变化更为敏感。而克里金法有一定的平滑效果,一方面有使平滑高大植被的效果,使其高程一定程度地降低(9月,σ=0.46),但另一方面,也会出现“牛眼”现象[28],造成与实际地貌的偏离。另外,未经滤波处理直接生成DEM时,克里金法比TIN法得到的平均沟长沟宽更大,但平均沟深更浅,而经过滤波后,克里金法比TIN法生成的沟长沟宽和沟深更小。可能是由于克里金法插值时,沟沿线附近、沟底的部分点云插值权重过低,导致沟沿线附近、沟底略缓,出现DEM平滑现象。
总之,TIN法生成的DEM沟长沟宽等参数更接近于真实的切沟地貌形态。经过滤波之后比未滤波生成DEM的切沟形态参数有所增大,说明具有一定地植被去除效果。
1)植被季相变化对于三维激光扫描具有一定影响,生长季测量的点云高程均值高于非生长季0.25 m,从而进一步影响到DEM的精度和地形地貌刻画。
2)基于滤波-TIN法生成的DEM平均沟长、沟宽和沟深,2次相比,分别增大0.03、0.25和0.10 m,结合2次实地抽样的植物株高,说明坡度滤波能一定程度上去除植被。但滤波后9月份(生长季)点云高程均值仍然比12月份(非生长季)高出0.15 m,两者并不完全一致,表明植被季相变化对DEM所造成的影响无法完全去除。
3)与克里金法相比,TIN法生成的DEM精度较高。前者存在平滑现象,植被生长季时能粗略地去除高大植被,而在非生长季时平滑沟道边缘。TIN法比克里金法的插值结果更接近实际地貌。
黄土高原丘陵沟壑区植被季相变化对DEM精度具有影响,且仅通过植被滤波算法很难完全去除。选择植被非生长季时进行测绘,并使用适合的滤波插值方法,才能有效地提高DEM的精度。