SfM摄影测量法测定切沟的精度评价

2020-04-27史扬子黄婷婷罗建勇刘宝元

水土保持研究 2020年1期

史扬子, 黄婷婷, 罗建勇, 杨扬, 刘宝元

(1.北京师范大学地理科学学部地表过程与资源生态国家重点实验室, 北京100875;2.中国科学院水利部水土保持研究所黄土高原土壤侵蚀与旱地农业国家重点实验室, 陕西杨凌 712100)

切沟是指普通农具无法横过耕作的侵蚀沟，宽度和深度一般均超过50 cm[1-2]。切沟侵蚀是在切沟发生发展过程中造成的侵蚀，主要包括沟头溯源侵蚀、沟底下切侵蚀和沟岸侧向侵蚀3种方式[2-3]。切沟侵蚀是全球河流泥沙的重要来源，其产沙量可占流域产沙总量的10%～94%[4]。为了深入研究切沟侵蚀规律，指导流域及区域水土流失综合治理，亟需建立一套简便、快速且行之有效的切沟监测方法。

切沟测量方法可大致分为传统方法和现代方法两类。传统方法以断面测量法为代表，即沿切沟纵剖面将其划分成若干小段，利用卷尺分别测量各段的横断面面积及纵向长度，进而获得切沟总体积。断面测量法原理简单、成本低廉、操作便捷，是切沟监测的常用方法[5-7]，但在切沟分段与测量时受卷尺精度与人为因素的影响；且只能反映切沟的整体变化状况，无法获知切沟的空间发展方向[5,8]。现代方法以三维激光扫描法、GPS(Global Positioning System)测量法和摄影测量法为典型代表。三维激光扫描法是指通过激光脉冲采集目标物体的点三维信息，经后期处理转换为三维模型的方法，因其具有高精度和高效率的优势，被广泛认为是目前最精确的切沟测量方法[9-11]。但该方法对设备要求高，资金投入大，且测量人员需具有一定的专业测量知识[9]，往往难以大面积推广。GPS测量法近年来发展迅速，其中，实时动态差分(Real-time Kinematic，RTK)GPS作为GPS测量技术发展的一项全新突破，兼具定位精度高(cm级)与作业速度快的优势，被成功应用于切沟监测[12-14]。并且，鉴于其在精度方面固有的优越性，RTK GPS的测量结果也可作为基准数据校验其他监测方法[14]。然而，该方法的应用受卫星状况限制以及障碍物的影响，且具有成本高、对测量人员的仪器操作水平要求较高等不可避免的缺点[13,15]。

摄影测量法可高效快速地获得地表形态变化，于20世纪60年代开始应用于土壤侵蚀监测[16]。彼时的摄影测量方法具有较大的局限性，如需要成本高、笨重的量测相机拍摄照片，需要昂贵的软件重建三维点云模型，后期处理数据所需时间长，受人为影响较多等[10]。近年来，随着计算器视觉技术的高速发展，基于运动恢复结构(Structure from Motion，SfM)的摄影测量方法应运而生。该方法从运动摄像机拍摄的多幅二维图像中估计摄像机运动并重建所摄地物的三维场景结构[17-19]，然后借助地面控制点(Ground Control Point，GCP)的坐标信息获取被测地物的确切形态[20]。因具有快速、简单、成本低等优势，SfM摄影测量法自20世纪80年代发展以来便广泛应用于包括测绘、城市建设、考古研究与文物保护在内的各学科领域，但其在地学的应用起步较晚[19]。近10 a来才开始逐渐应用于地貌与侵蚀测量[9,21]。

目前，SfM摄影测量法在国内切沟监测中应用较少，且多集中于较小尺度的室内模拟侵蚀沟测量[22]，缺乏野外切沟测量实践。鉴于三维激光扫描仪高昂的成本，本研究选取河北省滦平县两间房村的一条典型切沟，以RTK GPS的测量结果为基准，分析断面测量法与SfM摄影测量法对该切沟的测量精度，评价SfM摄影测量法在切沟测量中的适用性。并在此基础上，设计不同的GCP布设方案，探讨其对切沟测量精度的影响，为SfM摄影测量法的实施提供参考依据。

1 研究方法

1.1 测量切沟概况

目标切沟沟长约15 m，宽2～4 m，深0.5～1.5 m，横断面多呈V形，位于河北省东北部承德市滦平县两间房村(图1)，属中温带向暖温带过渡的半干旱半湿润大陆性季风气候，四季分明，冬长夏短。多年平均降雨量545 mm，多年平均气温7.7 ℃。2018年4月，对该切沟开展测量。同时，在沟缘与沟底采集土壤表层样品，测得其平均容重为1.42 g/cm3，有机碳平均含量为1.04%，砂粒平均含量70.54%，粉粒23.39%，黏粒6.07%，质地属砂壤土[23]。

1.2 RTK GPS测量法

采用水平测量精度为2.5 mm±2 ppm，垂直测量精度为5 mm±2 ppm的中海达iRTK2经典版RTK GPS，以大致40 cm的间距测量目标切沟沟缘、沟坡与沟底各点的地理坐标。同时，在地形变化剧烈处，适当增加测量点。目标切沟共测量点数248个，将其对应坐标导入ArcGIS 10.2中重建DEM，并提取与计算切沟长度、平均宽度、平均深度、周长、面积和体积等参数。

图1 目标切沟实地照片

1.3 断面测量法

根据切沟横断面形状(近似三角形或梯形)与面积将其划分为7段，使用卷尺测量每个横断面的各边边长，用以计算各横断面的面积；同时利用卷尺测量切沟各段的纵向长度，计算各段体积，加和得到切沟总体积[5]：

(1)

式中：Si为第i个横截面的面积(m2)；Hi为第i段切沟的纵向长度(m)；n为切沟分段总数，此处取值为7；V为切沟总体积(m3)。

1.4 SfM摄影测量法

1.4.1 切沟照片获取与处理切沟照片利用Canon EOS 70 D相机获取。拍摄过程中，相机焦距设置为18 mm，并保证相邻两张照片的重叠率为60%以上[24]。共拍摄切沟照片678张，导入至SfM专业软件Agisoft Photoscan中。为防止如模糊等拍摄质量差的照片参与后期三维重建而影响数据精度，通过该软件评估照片中最聚焦部分的清晰度级别，将该级别低于0.5的照片删除[25]。最终，共有635张有效照片参与目标切沟的三维重建。

1.4.2 SfM数据处理将有效照片导入至Agisoft Photoscan进行目标切沟的三维重建时，先通过尺度不变特征转换算法(scale-invariant feature transform，SIFT)[17]完成特征点的提取与描述，再使用随机采样一致性算法(random sample and consensus，RANSAC)[26]过滤最近邻匹配(nearest neighbor，NN)[17]产生的错误匹配点，接着采用光束法平差(bundle adjustment，BA)[27]进行非线性优化，从而建立目标切沟的稀疏点云。在此基础上，利用面片的多视图立体视觉算法(patch-based multi-view stereo，PMVS)[27]重建稠密点云，进而生成目标切沟的DEM。最后将DEM导入ArcGIS 10.2中，提取与计算切沟相关参数。

1.4.3 GCP布设与坐标计算在切沟沟缘与沟底依据以下3条原则布设十字形GCP：(1) GCP平均间隔大致为3 m，均匀分布；(2) 沟缘转折处优先或适当加密；(3) 选取局部地势平坦、易固定、没有植物遮挡的位置。目标切沟共布设GCP 18个，其中沟缘14个、沟底4个。

利用中海达iRTK2经典版RTK GPS测量各GCP的大地坐标。为了减轻SfM摄影测量法的成本与野外负担，本研究同时使用小巧轻便的Leica D510激光测距仪(Laser Range Finder，LRF)测量GCP两两之间的空间直线距离及两点所在直线与水平面间的夹角，并由此计算GCP之间的高程差及其空间直线距离在水平面的投影距离，进而推算各GCP的空间直角相对坐标。方法如下：将各GCP依次连线构成三角网，取水平面为水平方向，自定义坐标原点(0，0)与x，y轴方向；取垂直于水平面方向为z轴方向，高程增加方向为正，并自定义z轴0点，构建空间直角坐标系。在水平面上，根据三角网各边的投影距离计算各GCP的平面坐标(x，y)；在z轴方向，根据GCP之间的高程差计算各GCP的坐标z值。

为了校验LRF对GCP的测量精度，将RTK GPS所测各GCP的大地坐标转换为空间直角坐标后求取GCP两两之间的空间直线距离与高程差，与LRF的结果进行对比。

1.4.4 GCP布设方案比较 GCP的布设、测量与坐标计算复杂繁冗。因此，在原有18个GCP的基础上，对沟缘和沟底的GCP数量进行不同程度的抽稀。具体原则为沟缘GCP不抽稀、均匀抽取1/2，均匀抽取1/3；沟底不抽稀、均匀抽取1/2，全部缺失，共9种组合(图2)。将不同布设方案的控制点信息输入SfM专业软件Agisoft Photoscan 中，重建目标切沟的三维形态，比较各方案的重建效果。

图2 不同布设方案的GCP空间分布

2 结果与分析

2.1 LRF的GCP测量精度

图3对比了LRF与RTK GPS所测GCP空间直线距离与高程差。结果表明，无论是空间直线距离还是高程差，LRF的测量结果与RTK GPS十分接近，对应散点均靠近1∶1线分布。此外，两种方法得到的空间直线距离相关系数为0.993，高程差相关系数为0.999，均在0.01的置信水平显著。可见，LRF可有效测量GCP的空间分布与相对位置，为获取GCP坐标与重建切沟三维形态提供了另外一种可能方法。相比RTK GPS，LRF便宜、轻便，测量过程高效易行，可极大地减少SfM方法测量切沟时的设备费用，减轻野外作业负担。

2.2 SfM摄影测量法的切沟测量精度

2.2.1 RTK GPS、断面测量法、SfM摄影测量法所测切沟参数比较表1为RTK GPS、断面测量法与SfM摄影测量法测量得到的目标切沟主要形态参数。其中，RTK GPS通过其所测各点的坐标信息重建DEM(图4A)后提取切沟长度、平均宽度、平均深度、体积、面积与周长；断面测量法主要测量切沟各段的横断面面积与其代表长度，用以估算切沟长度、平均宽度、平均深度与体积，并未测量切沟面积与周长；SfM摄影测量法分别利用RTK GPS所测的大地坐标(SfMRTK GPS)与LRF所测的独立空间直角坐标(SfMLRF)重建DEM(图4B，C)，提取得到对应的两组结果。

除平均宽度外，断面测量法相比RTK GPS测得的各参数误差均较大。长度、体积分别高估了15.75%与37.28%，平均深度低估了17.12%。因为平均深度的影响，宽深比较RTK GPS高估了29.91%。这主要是因为断面测量法法概化了切沟形态，忽视了沟内细节，且切沟测量横断面的位置选取与卷尺读数受人为因素影响较大。相比于尹佳宜等[5]在东北黑土区切沟体积测量中得到的～10%的误差，本研究的断面测量法误差较大，这可能与选取切沟的规模有关。前者测量了3条切沟，体积均超过500 m3，在所用卷尺精度一致的情况下，总体积越大，误差百分比相应越小。

相比之下，SfM方法对切沟的测量精度高得多。由于SfMRTK GPS与SfMLRF均基于重建切沟DEM的坡度转折提取其沟缘线，所得沟缘线轮廓一致(图4B，C)，对应的长度、平均宽度、周长与面积也相等。其中，长度相比RTK GPS仅低估了0.56%，平均宽度与面积分别高估了8.26%与1.72%。就体积而言，SfMRTK GPS得到的结果为26.85 m3，相对RTK GPS的误差为2.40%；SfMLRF的结果为26.28 m3，误差仅为0.23%。无论基于何种GCP坐标，SfM摄影测量法均表现出较高的切沟体积测量精度，极大地优于断面测量法(表1)。然而，SfM摄影测量法提取的切沟周长却远大于RTK GPS的结果，前者比后者长8.10 m，增幅达23.63%。这是由于RTK GPS得到的是密度较低的点数据，仅可体现沟缘的大致走向(图4A)；而SfM法基于重建的三维点云得到DEM，由此提取的坡度信息可使沟缘走向变得清晰，细节得以充分体现，经矢量化的沟缘线更加细致，周长更长(图4B，C)。

图3 LRF与RTK GPS所测GCP两两之间空间直线距离与高程差对比

表1 RTK GPS、断面测量法与SfM摄影测量法测定的

图4 不同方法所测切沟DEM

2.2.2 SfM摄影测量法与RTK GPS所测切沟DEM对比虽然SfMRTK GPS与RTK GPS均基于大地坐标开展DEM的重建与计算，但SfMLRF采用独立空间直角坐标系，无法与RTK GPS所测切沟DEM进行直接比较。因此，将所测切沟DEM与对应沟缘面DEM作相减运算，获得切沟相对DEM后进行SfM方法与RTK GPS的比较。

图5为SfMRTK GPS，SfMLRF所得切沟相对DEM与RTK GPS所测结果的差值。显然，无论基于何种坐标系，SfM方法得到的切沟相对DEM比之RTK GPS结果总体呈西北高、东南低的态势。这可能是因为试验当天的切沟照片拍摄于上午，切沟西北岸光线充足，深处细节曝光充足，进行三维重建的三角测量计算时，相机与所摄地物的距离较实际距离偏小，高程偏大，与GPS结果相减呈正值。而切沟东南岸光线较暗，物体远近对比较不明显，相机与地物的距离较实际距离偏大，后期计算的高程偏低，与GPS结果相减呈负值。可见，摄影测量法的使用一定程度上受限于天气条件，需确保所测地物所有部分的光线条件基本一致[28]。

图5 (A)SfMRTK GPS、(B)SfMLRF与RTK GPS所测切沟相对DEM差值

此外，基于两种方法得到的目标切沟相对DEM，对各栅格高程作配对t检验，结果显示，SfMRTK GPS，SfMLRF均与RTK GPS存在显著差异(p<0.01)，但总体差异很小，前者与RTK GPS的平均高程差值仅为-0.02 m，后者与RTK GPS的平均高程差值略高，也仅为-0.04 m。进一步对比DEM各栅格的高程差值也发现，SfMRTK GPS与RTK GPS所测切沟相对DEM高程差值均不超过1.00 m，且主要集中于-0.2～0.2 m范围内，对应栅格数量占比73.72%；SfMLRF与RTK GPS的DEM相对差值不超过0.9 m，同样集中于-0.2～0.2 m范围内，对应占比74.58%(图6)。这一结果说明SfM摄影测量法与RTK GPS所测切沟DEM总体差别不大，是一种高精度的切沟测量方法。同时，LRF也可有效替代RTK GPS应用于GCP的测量及后期切沟三维形态的重建。

2.3 不同GCP布设方案对SfM摄影测量法测量精度的影响

图7为不同沟缘与沟底GCP抽稀方案所得切沟DEM与基于所有18个GCPs(未抽稀)所得DEM的差值。由于LRF所测GCP两两之间的空间距离与高程差与RTK GPS结果无显著差异，且由此获得的切沟DEM与RTK GPS结果差别也不大，此处各方案所取GCP坐标为LRF测得的独立空间直角坐标。总体而言，不同程度的沟缘与沟底GCP抽稀对生成的切沟DEM影响不大。尽管GCP沟底不抽稀、沟缘均匀抽取1/2大范围地低估了切沟高程(图7D)，GCP沟底全部缺失、沟缘均匀抽取1/3大面积高估了切沟高程(图7I)，但其低估与高估程度均不超过0.15 m。

图6 SfMRTK GPS与SfMLRF相对RTK GPS所测切沟相对DEM差值累积频率百分比

基于各抽稀方案与未抽稀方案得到的切沟DEM，对各栅格高程作配对t检验，结果显示，不同方式的GCP抽稀均显著改变了各栅格高程(p<0.01)，但总体幅度不大。各抽稀方案与未抽稀情况的DEM栅格平均高程差值介于-0.03～0.01 m范围内。进一步统计DEM各栅格的高程差发现，各抽稀方案所得切沟DEM相比未抽稀结果的高程差值主要集中于-0.05～0.05 m范围内，对应栅格数量占比均超过70%(图8)。此外，各GCP抽稀方案得到的目标切沟体积介于26.43～27.10 m3，相比未抽稀方案的误差介于0.56%～3.13%(表2)。这一结果说明SfM摄影测量法基于不同GCP布设方案得到的切沟DEM总体误差较小，均不超过5%。因此，即使GCP数量低至4个，也可满足大多数情况下的切沟体积测量精度。

图7 SfMLRF基于不同GCP布设方案所得切沟DEM差值

图8 SfMLRF基于不同GCP布设方案所得切沟DEM差值的累积频率百分比

一般来说，GCP数量越多，所摄地物的三维形态越精确[20,29]。并且，GCP是否均匀分布、是否在地势剧烈处进行加密，也会影响SfM方法的测量精度[30]。但本研究中切沟DEM或体积精度并未随GCP数量，特别是沟底GCP数量明显变化。这可能与目标切沟沟底较浅、规模较小有关。一方面，目标切沟沟宽约3 m，长约15 m，少量的GCP便可满足精度需求；另一方面，沟缘最高处与沟底最低处高程差不足7 m(图4)，在沟缘拍摄照片时也可将沟底拍摄在内，沟底的尺度信息可在后期的三维重建过程中通过提取照片中的其他地物特征信息获得，沟底有无GCP，GCP数量多少对切沟DEM重建影响较小。但当研究对象为如黄土高原地区广泛分布的体积、高差较大的切沟时，还需进一步探究GCP布设方法对切沟测量的影响。