SfM摄影测量法对黄土高原典型切沟的测量精度评价

陈薪伊,史扬子,杨 扬*,黄婷婷,刘瑛娜,张光辉

(1 北京师范大学 地理科学学部 地表过程与资源生态国家重点实验室,北京 100875;2 北京师范大学 地理科学学部 地理学院,北京 100875;3 水利部 水土保持监测中心,北京 100053)

切沟是指普通农具无法横跨耕作的侵蚀沟,通常切破犁底层,深入母质,宽度和深度均超过50 cm[1-2]。切沟侵蚀是指切沟发生发展过程中导致的侵蚀。作为水力侵蚀最严重的阶段,切沟侵蚀是全球河流泥沙的主要来源之一,其产沙量可占流域产沙总量的10%～94%[3]。在我国黄土高原地区,切沟侵蚀产沙量可占小流域总产沙量的50%以上,部分地区甚至可以达到80%以上[4-5]。选择准确、高效的测量方法对切沟及切沟侵蚀进行监测,对于揭示切沟侵蚀时空分布规律、阐明切沟侵蚀机理、构建和完善切沟侵蚀模型具有重要意义。

近年来,以激光扫描法(laser scanning,LS)、实时动态差分(real-time kinematic,RTK)GPS和摄影测量法为典型代表的现代土壤侵蚀测量技术开始逐渐应用于切沟及切沟侵蚀监测,均取得了良好效果[6-9]。LS利用激光扫描设备发射激光,根据其返回时间推算所测地物的三维坐标,进而生成DEM。该方法测量精度极高[10],结果常被用作验证其他方法的标准数据[11];然而,LS所用设备昂贵,对操作人员的专业要求较高,目前应用十分有限[12]。RTK GPS是随GPS定位精度不断突破发展而来的一项全新测量技术,兼具定位精度高(厘米级)、测量耗时短等优势[13-14];但是,该方法的实际应用受到卫星信号限制,成本较高且要求操作人员到达实地测量点,从而在一定程度上限制了其广泛应用。摄影测量法可准确测量地表形态变化,早在1980年代末便被用于切沟侵蚀研究[12]。传统的立体摄影测量法需专业人员设计专门拍摄路线并利用笨重、昂贵的量测相机对切沟进行连续拍照,总体应用有限。随着计算器视觉技术的快速发展,尤其是尺度不变特征转换等图像特征检测和匹配技术的出现,基于运动恢复结构(structure from motion,SfM)的摄影测量方法应运而生[11]。该方法利用一部普通相机从多角度对切沟进行连续拍照,然后借助若干已知地面控制点坐标便可重建所测地物的三维模型。尽管在切沟监测中的应用较晚[8],但该方法因具有便携、快速、成本低和操作简单等优点,发展势头迅猛[12,15-16]。

目前,SfM摄影测量法在国内切沟监测中的应用还较少,沟壑纵横、切沟侵蚀严重的黄土高原也不例外。本研究在黄土高原六道沟小流域选取规模不等、形态各异的7条典型切沟,以LS和RTK GPS的测量结果为基准,评估SfM摄影测量法对切沟和切沟侵蚀的测量精度,以期为SfM摄影测量法在黄土高原的实践提供参考。

1 材料与方法

1.1 研究区概况

研究区为陕西省神木市的六道沟小流域(110°21′～110°23′ E,38°56′～38°51′ N),总面积6.9 km2,主沟道长4.2 km。该区地处毛乌素沙地边缘,属水蚀、风蚀交错带(图1a);气候为中温带半干旱气候,多年平均降水量为437.4 mm,主要集中在夏季;植被为沙生植物与黄土区植物交错出现;区内地面组成物质主要为第四纪黄土沉积物,土壤类型包括沙黄土、黄绵土、风沙土及坝地淤土,土壤结构松散,易遭受侵蚀[17-18]。

图1 六道沟小流域(a)及所选切沟(b)的地理位置Fig.1 Geographical locations of Liudaogou catchment (a) and gullies selected (b) for investigation

2018年,在六道沟小流域开展实地调查,选取不同规模和形态的代表性切沟共7条(图1b),于植被尚未完全萌发的初春(4月底、5月初)及植被枯萎凋落的秋末(10月)进行测量。在春季,利用SfM摄影测量法对1～6号切沟(记为1S、2S、3S、4S、5S、6S)进行测量。由于三维激光扫描及RTK GPS设备均为外借,时间所限,仅测量了其中部分切沟;前者对1～4号及6号共5条切沟进行了测量、后者只对规模最大及中等的3号和5号切沟进行了测量。在秋季,同时利用SfM摄影测量法和RTK GPS对5号和7号切沟(记为5F、7F)进行测量。在7条切沟中,仅5号切沟在春、秋季均利用RTK GPS进行了测量,通过对比可计算该切沟在当年雨季的侵蚀体积。将该计算结果作为参考基准,评估SfM摄影测量法对切沟侵蚀的测量精度。

1.2 野外切沟测量

1.2.1 SfM摄影测量法

进行SfM摄影测量之前,在切沟沟缘、沟底的局部地势平坦且无植被遮挡处布设十字形地面控制点(ground control point,GCP),各GCP均匀分布,在沟缘转折处适当加密。利用Leica D510激光测距仪测量相邻GCP之间的空间直线距离及两点连线与水平面之间的夹角,在此基础上计算两点之间的水平投影距离和相对高程,再推算各GCP的空间直角相对坐标[19]。

各切沟利用Canon EOS 70D相机进行拍摄,相机焦距设置为18 mm。为确保后期三维重建的顺利进行,相邻照片的重叠率至少为60%[20]。拍摄照片数量因切沟规模而异,介于1 320～11 600张(表1)。将各切沟照片导入软件Agisoft Metashape 1.8.2中,将清晰度级别不足的照片删除后进行三维重建,依次生成稀疏点云、稠密点云和0.01 m×0.01 m分辨率的DEM[21-22]。

表1 各切沟拍摄照片数量Tab.1 Number of images taken for each gully

1.2.2 LS

利用激光扫描法进行野外实地测量时,先在所测区域内选择宽阔、无地物与植被遮挡的位置布设棱镜靶标,利用Leica Nova MS60三维激光扫描全站仪进行地理坐标校正,至少重复3次。校正结束后,利用该全站仪扫描切沟不同部位,并在地形剧烈变化处适当增加扫描次数。将扫描得到的点云数据导入ArcGIS 10.2,基于不规则三角网模型(triangulated irregular network,TIN)构建DEM。为与SfM摄影测量法结果进行对比,利用重采样生成0.01 m×0.01 m分辨率的DEM。

1.2.3 RTK GPS

利用中海达iRTK2经典版RTK GPS以大致50 cm的间距测量各切沟的沟缘、沟坡及沟底,并在地形剧烈变化处适当增加测量点。与LS类似,将点云数据导入ArcGIS 10.2中构建DEM,并重采样至0.01 m × 0.01 m分辨率。

1.3 切沟相对DEM计算

3种测量方法采用的坐标系不同,得到的切沟DEM无法直接进行比较。RTK GPS采用大地坐标系;LS和SfM摄影测量法的点云数据虽然均采用独立空间直角坐标系,但二者的坐标原点并不完全一致。因此,先使用ArcGIS 10.2中的Georeferencing工具,以LS或RTK GPS得到的数据为基准,对SfM摄影测量法所得DEM进行配准。初春的植被虽较稀疏,也会对切沟沟缘的识别造成一定的影响。在春季,使用3种方法对3号切沟进行了测量,计算其相对DEM之前,以RTK GPS数据为基准对LS和SfM摄影测量法所测数据进行配准;接下来,结合目视解译和局部坡度变化识别不同测量方法得到的各切沟沟缘,利用TIN生成各切沟沟缘面DEM并重采样至0.01 m × 0.01 m分辨率,将各切沟DEM减去相应的沟缘面DEM得到相对DEM后再行比较。在相对DEM中,沟底的高程为负值,且绝对值大;沟缘的高程绝对值小,在理想状况下应接近0。然而,沟缘经矢量勾绘后转化为栅格并重采样得到,在与切沟DEM做相减计算时往往会出现对应栅格高程不等于0的情况。

1.4 切沟形态参数提取与对比

切沟各形态参数通过沟缘面DEM和切沟相对DEM在ArcGIS 10.2中加以提取。基于沟缘面DEM可提取切沟长度、周长、面积和周长-面积比[23],基于切沟相对DEM可得到平均表面宽度、平均深度和体积。利用相对误差(E)评估SfM摄影测量法对各切沟形态参数的测量精度[9]:

(1)

式中:mi为SfM摄影测量法所测第i条切沟的形态参数;Mi为LS或RTK GPS所测第i条切沟的形态参数。

2 结果与分析

2.1 切沟测量精度

2.1.1 SfM与LS所测切沟DEM对比

图2对比了春季SfM与LS所测各切沟的相对DEM。1号(1S)、2号(2S)、3号(3S)、4号(4S)、6号(6S)切沟的高程平均差值分别为-0.077、0.120、-0.162、0.436和0.373 m。

图2 SfM(a～e)与LS(f～j)所测1S、2S、3S、4S和6S相对DEMFig.2 Relative DEM of 1S, 2S, 3S, 4S and 6S measured by SfM (a～e) and LS (f～j)注:网络版为彩图。

然而,从高程差值的累积频率分布(图3)与空间分布(图4)来看,2种方法对1S和3S的测量差异总体较小,相对高程差值主要集中在-0.5～0.5 m,对应栅格占比均在85%左右。SfM较LS总体高估了4S和6S的高程,各栅格高程差值多为正值,主要集中在0～1.5 m范围内,对应栅格占比分别为70.2%和67.8%,介于-0.5～0.5 m的栅格占比分别为41.3%和43.2%。至于2S,尽管2种方法所测平均高程差值较小,但空间变异较大,相对高程差值的分布较分散(图3)。高程差值主要集中在-2～2 m,对应栅格占比为70.0%,介于-0.5～0.5 m的栅格仅占22.0%。

图3 SfM与LS所测不同切沟相对高程差值累积频率分布Fig.3 Cumulative frequency distribution of relative elevation differences between SfM and LS for different gullies

图4 SfM与LS所测1S(a)、2S(b)、3S(c)、4S(d)和6S(e)相对高程差值空间分布Fig.4 Spatial distribution of relative elevation differences of 1S (a), 2S (b),3S (c), 4S (d) and 6S (e) between SfM and LS注:网络版为彩图。

2.1.2 SfM与RTK GPS所测切沟DEM对比

利用SfM和RTK GPS测量春季3号(3S)和5号(5S)切沟,以及秋季5号(5F)和7号(7F)切沟的相对DEM(图5)发现,2种方法所测各切沟相对高程的平均差值不大,依次为-0.033、-0.096、-0.128和-0.075 m。从累积频率分布(图6)与空间分布(图7)来看,SfM与RTK GPS所测5S、5F和7F的高程差值均主要集中在-0.5～0.5 m,对应栅格占比分别为95.6%、96.3%和94.7%,超过±1.0 m的栅格仅占0.1%、1.2%和0.3%。相比之下,SfM和RTK GPS对3S的测量差异略大,高程差值介于-0.5～0.5 m的栅格较少,对应占比为85.3%,甚至还有3.8%的栅格高程差值超过±1 m。

图5 SfM(a～d)与RTK GPS(e～h)所测3S、5S、5F、7F相对DEMFig.5 Relative DEM of 3S, 5S, 5F and 7F measured by SfM (a～d) and RTK GPS (e～h)注:网络版为彩图。

图6 SfM与RTK GPS所测不同切沟相对高程差值累积频率分布Fig.6 Cumulative frequency distribution of relative elevation differences between SfM and RTK GPS for different gullies

图7 SfM与RTK GPS所测3S(a)、5S(b)、5F(c)、7F(d)相对高程差值空间分布Fig.7 Spatial distribution of relative elevation differences of 3S (a),5S (b), 5F (c) and 7F (d) between SfM and RTK GPS注:网络版为彩图。

2.1.3 切沟形态参数对比

SfM、LS和RTK GPS所测各切沟形态参数如表2所示。配对t检验结果显示,无论与LS还是与RTK GPS相比,SfM所测各切沟形态参数均无显著差异(P>0.05)。

表2 SfM、LS和RTK GPS所测各切沟形态参数Tab.2 Morphological parameters of each gully measured by SfM, LS and RTK GPS

相比LS,SfM所测各切沟长度和周长精度最高,相对误差均不超过±6%,平均相对误差仅分别为2.7%和1.9%(表3)。SfM对切沟平均表面宽度的测量精度也较高,除6S的相对误差高达11.8%外,各切沟的相对误差均不超过±7%。相比之下,SfM对切沟深度的测量误差较大。例如,SfM较LS大面积高估了4S和6S的高程(图4d、4e),其对平均深度的低估百分比分别高达20.2%和32.6%。对于其余3条切沟,SfM对平均深度的测量误差介于-4%～13%。尽管SfM所测各切沟面积与LS结果相近,相对误差均不超过±7%;由于对深度的测量误差,除1S和2S外,SfM对3S、4S和6S的体积测量误差均超过±10%,分别为12.5%、-20.7%和—37.4%。两种方法测得的各切沟周长-面积比也十分接近(表2),但由于LS所测周长-面积比本身较小,部分切沟的相对误差超过±10%(表3)。

表3 SfM所测不同切沟形态参数的相对误差Tab.3 Relative errors of different gully morphological parameters measured by SfM 单位:%

对同时利用RTK GPS和SfM 测量的3S、5S、5F和7F各形态参数进行统计分析发现,2种方法所测结果总体相差不大(表2)。与SfM和LS的对比类似,SfM对各切沟长度、平均表面宽度、周长、面积和周长-面积比的测量误差较小,均不超过±8%,平均相对误差依次为2.9%、1.3%、-1.9%、1.4%和-3.2%;对深度和体积的测量误差较大,但也不超过±15%,平均相对误差分别为7.9%和9.0%(表3)。

2.1.4 切沟测量精度影响因素

与LS和RTK GPS相比,SfM对切沟的测量误差主要与植被覆盖、切沟沟缘曲折程度、光照条件以及沟坡形态有关。尽管开展野外调查的初春和秋末植被较为稀疏,但其也会对LS和SfM的测量结果造成干扰,影响切沟沟缘的识别和DEM重建,使2种方法所测DEM无法完全重合、相对高程出现差异(图4)。相对而言,RTK GPS直接对地表进行测量,基本不受植被和枯落物干扰[12]。当SfM与其进行对比时,往往会高估地表高程。例如,3S沟头附近的沟底发育了1条小沟,但由于沟底植被干扰,SfM难以对其进行完整识别和测量,导致SfM在此所测高程偏大(图8)。

图8 SfM与RTK GPS所测3S相对高程差值空间分布(a)及沟底小沟实地拍摄照片(b)Fig.8 Spatial distribution of relative elevation differences of 3S between SfMand RTK GPS (a) and a photo of a small channel at the bottom of 3S (b)注:网络版为彩图。

总体而言,切沟的周长-面积比越大,沟缘越曲折,SfM的测量误差相应越大。例如,在SfM和LS同时开展监测的5条切沟中,6S的沟缘最为曲折,LS所测周长-面积比高达1.19 m-1,SfM所测面积和体积误差均最大(表2)。此外,6S因拍摄过多照片造成的信息冗余(表1),可能也是造成SfM测量误差的原因之一[24]。而在与RTK GPS的对比中,7F的周长-面积比最大(表2),SfM所测切沟长度、平均表面宽度和面积误差均最大。

SfM对切沟的测量精度还受到光照条件的显著影响。例如,SfM和LS所测4S的高程差值呈西低东高的态势(图4d),这主要是因为SfM测量当天下午光线充足,斜射的阳光只可到达东岸,无法进入较深的沟底。东岸的细节充足曝光,相机与所摄地物的距离较实际偏小,进行三维重建的三角测量计算时得到的高程偏大(沟深被低估),与LS结果相减为正值。相比之下,西岸光线较暗,物体远近对比较不明显,相机与地物的距离较实际偏大,SfM后期计算得到的高程偏小(沟深被高估),与LS相减为负值。可见,SfM的应用受限于天气情况,最好在阴天或多云天气进行,确保所测地物各部位的光线基本一致[19,25]。

沟坡形态主要影响SfM与RTK GPS对切沟的测量差异。对于沟坡上某些坡度较大的部位,RTK GPS无法直接进行测量,只能通过沟缘和沟底的测量结果插值获得。如果沟坡并非直形坡,而呈凸形、凹形等不规则形状,那么SfM与RTK GPS所测高程就会出现较大差异。SfM与RTK GPS共同测量的各切沟沟坡均出现了局部坡度超过60°的陡峭凸出部位(图9),SfM所测相对高程明显较RTK GPS大(图7)。可见,相比RTK GPS,SfM对地形陡峭、无法到达地物存在明显的测量优势[12]。然而,SfM对地物的三维重建高度依赖GCP的坐标信息,对于黄土高原某些沟深坡陡的切沟(如2S)来说,GCP只能布设在沟岸和沟底,可能会在一定程度上影响SfM对切沟的测量精度(图3、图4b)[19]。

2.2 切沟侵蚀测量精度

2.2.1 SfM与RTK GPS所测切沟侵蚀对比

将SfM和RTK GPS所测5F的DEM减去5S的DEM,得到5号切沟在当年雨季的土壤侵蚀空间分布(图10)。2种方法都发现该切沟的沟头处发生侵蚀,沟尾处出现堆积,但SfM测得的切沟沟底高程变化明显较RTK GPS复杂。此外,SfM所测切沟沟缘尤其是东侧和北侧沟缘的高程变化较RTK GPS所测结果大,且主要表现为高程下降。从切沟形态参数(表4)来看,相比RTK GPS,SfM所测切沟平均表面宽度和体积变化较大,但周长变化明显偏小。对于切沟长度、平均深度和面积变化,2种方法所测结果相差不大。

表4 SfM和RTK GPS所测5号切沟形态参数变化Tab.4 Changes of morphological parameters of gully No.5 measured by SfM and RTK GPS

图10 SfM(a)与RTK GPS(b)所测5号切沟侵蚀空间分布Fig.10 Spatial distribution of soil erosion measured by SfM (a) and RTK GPS(b) for gully No. 5注:网络版为彩图。

2.2.2 切沟侵蚀测量精度影响因素

与切沟监测类似,SfM对切沟侵蚀的测量精度主要与植被及其枯落物覆盖有关。春季萌发的植被会对沟缘线造成遮挡,当使用SfM测量时,解译得到的沟缘高程偏大且易向切沟内侧偏移,尤其是植被较为丰富的东侧和北侧沟缘;到了秋季,植被枯萎、叶片凋落,切沟沟缘出露。因此,相比基本不受植被干扰的RTK GPS,SfM所测切沟沟缘高程和平均表面宽度变化均较大(表4),因2种方法所测切沟面积变化相差不大,因而SfM所测切沟体积变化相应较大(表4)。然而,雨季过后,SfM所测切沟周长仅增加0.21 m,远小于RTK GPS所测结果。当使用RTK GPS进行测量时,5号切沟在雨季前后的周长-面积比基本相当,可见侵蚀并未显著改变沟缘的曲折程度。2种方法测得的切沟周长差异可能主要还是与沟缘植被有关。雨季前,由于沟缘植被干扰,SfM所测切沟沟缘较曲折,周长-面积比为0.66 m-1,周长与RTK GPS测量结果相差不大(表2)。雨季过后,尽管侵蚀使切沟有所扩张,但植被枯萎凋落后沟缘曲折程度降低,SfM所测周长-面积比降至0.62 m-1(表2),该方法测量得到的周长增加有限(表4)。因设备所限,本研究仅对1条切沟的雨季侵蚀进行了系统测量。未来有必要增加切沟测量条数,尤其是雨季期间土壤侵蚀明显的切沟,系统开展SfM对切沟侵蚀的测量精度评估。

2.2.3 切沟侵蚀影响因素

从SfM和RTK GPS所测5号切沟的相对DEM(图10)和各形态参数变化(表4)可以发现,该切沟在雨季前后总体变化不大。一方面,该切沟接近小流域分水岭,汇水面积小,侵蚀动力不大;另一方面,该切沟沟缘处的植被生长较好,沟底的枯落物也较丰富,在一定程度上可保护切沟,削弱降雨和径流侵蚀。

3 结论

本文以黄土高原六道沟小流域内7条规模与形态各异的典型切沟为研究对象,通过与LS和RTK GPS对比,评估了SfM摄影测量法对切沟和切沟侵蚀的测量精度。研究结果表明,SfM摄影测量法可作为一种简单、便捷的高精度方法应用于黄土高原的切沟和切沟侵蚀测量。与LS和RTK GPS相比,SfM摄影测量法所测各切沟相对DEM总体差异不大,平均高程差值介于-0.162～0.436 m;各切沟长度、表面宽度、周长和面积的相对误差基本不超过±7%;深度和体积的相对误差略大,介于-37.4%～18.0%。SfM对切沟的测量精度主要与植被覆盖、切沟沟缘曲折程度、光照条件以及沟坡形态有关。因设备用时所限,仅评估了其中1条切沟的雨季侵蚀,发现RTK GPS和SfM摄影测量法所测土壤侵蚀空间分布与切沟形态参数变化均十分接近。

致谢:诚挚感谢西华师范大学张斌教授、硕士研究生郑吉林和吴曦在野外调查工作中的协助。