基于三维重建技术的坡面细沟侵蚀演变过程研究

吴淑芳，刘勃洋，雷 琪，孙立全，郭慧莉，冯 喆，钱 阔

（西北农林科技大学旱区农业水土工程教育部重点实验室，杨凌 712100）

0 引 言

黄土高原丘陵沟壑区微地貌形态以沟缘线为界，划分为由梁峁顶面和斜坡构成的坡面及由切沟、冲沟、干沟等构成的沟道[1]。沟缘线以上的坡面在强降雨下常常伴随着溅蚀、片蚀、细沟侵蚀、浅沟侵蚀、切沟侵蚀等并存的水力侵蚀过程，其过程错综复杂。这类坡面既是流域侵蚀产沙的主要源地，又是控制水土流失、恢复与重建生态环境的基本治理单元[2]，对于其坡面侵蚀现象与侵蚀规律的探究，历年来都是黄土高原环境整治中理论性与实践性均很强的重要科学任务。但由于丘陵区坡面的复杂性，研究手段的薄弱，大多研究主要集中于坡面系统产沙量的定位观测[3]、某一侵蚀方式（如细沟、浅沟）下的侵蚀发育过程[4-5]、坡沟径流泥沙输移[6]、坡沟系统侵蚀产沙模型研究[7-8]。尤其作为严重沟头溯源侵蚀和水流汇集发源地的“梁峁坡面”，对其侵蚀演变的时空过程描述还未曾有过详细报道。实质上，这一特殊断面—梁峁坡面产流产沙对沟缘线以下的坡面侵蚀和沟道侵蚀有着重大影响，陈浩等[9]结合观测资料和模拟试验证明上方来水来沙显著强化整个坡沟系统产流产沙。

近年来，随着科学技术的发展，一些较为先进、可靠的侵蚀监测方法逐步出现，如 REE示踪法[10]、KINEROS2模型法[11]、机载和地面激光雷达、GPS测量和三维激光扫描技术等，使得细沟侵蚀过程研究从传统测量迈入精准测量阶段[12-13]。张鹏等[14]利用高精度GPS、三维激光扫描仪和测针板3种测量方法动态监测坡面细沟侵蚀形态演变过程，以此分析3种方法动态监测沟蚀演变过程的可行性；霍云云等[15]利用三维激光扫描仪和基于元胞自动机CA-Rill模型对黄土坡面细沟侵蚀动态发育过程进行了图形化描述和计算机可视化模拟；和继军等[16]进行室内人工模拟降雨试验，其借助三维激光扫描仪对坡面微地形进行扫描，对塿土和黄绵土的坡面细沟发育过程和水沙关系进行了研究总结。三维激光扫描技术精度较高，但原始点云数据量巨大，设备仪器价格昂贵，获得的数据缺乏表面颜色分布。随着计算机视觉和计算机图形学的发展，基于照片的三维重建技术逐渐受到大众的关注[17]，它基于数码相机及无人机等设备拍摄的照片完成三维模型和场景的重建，成本低廉，不受地点、空间限制，具有一定的优越性[18-19]，可结合环境条件灵活运用，对土壤侵蚀过程研究具有广泛的应用前景。

本研究通过选取黄土高原丘陵沟壑区纸坊沟流域典型梁峁坡面作为研究对象，构建实体模型，进行人工模拟降雨试验，利用基于三维重建技术的PhotoScan软件研究坡面侵蚀演变过程及形态特征，同时验证三维重建技术在土壤侵蚀研究领域的适用性。

1 材料与方法

1.1 试验材料及装置

本试验所用土壤取自陕北坡耕地耕作层（0～20 cm），土壤类型为黄绵土，将供试土样风干，用激光粒度分析仪测其机械组成，黏粒质量分数为102.1 g/kg，粉砂粒质量分数为 545.4 g/kg，砂粒质量分数为352.5 g/kg，按照国际制土壤质地划分标准判定土壤为粉砂质土壤。表1为供试土壤机械组成。

表1 供试土壤机械组成Table 1 Particle distribution of experimental soil

降雨模拟试验在西北农林科技大学中国旱区农业节水研究院进行，试验采用侧喷式人工降雨设备，共3组，喷头组合形式为对喷，试验前进行降雨率定，保证率定雨强偏差在设计雨强5%范围以内，降雨均匀度在80%以上。收集研究区域降雨资料及侵蚀产沙资料，根据研究区多年自然降雨数据，取暴雨事件中常见的最大 30 min降雨强度作为设计雨强[20]，约为90 mm/h。试验土槽为变坡土槽，规格：长12 m、宽3 m、高5.3 m。土槽下端设集流装置，用于收集径流、采集泥沙样品。

1.2 试验布设

由于黄土高原地区细沟侵蚀大多集中于 10°～30°的裸露坡耕地上[21]，同时参考黄土丘陵沟壑区坡面坡度分级数据，将梁峁坡面从坡顶到坡底分成 5°、10°、15°、20°、25°、35°。为保证降雨试验过程中土壤水分充分下渗，在试验土槽底部装填15 cm厚的石子和砂砾，在其上方装填陕北坡耕地黄绵土。装填时为保持土的自然结构状态不过筛不研磨，称取一定质量土样分层装填，每层装土10 cm，边填边夯实，使容重控制在1.25 g/cm3左右，以确保降雨试验过程下垫面的均匀性。将构建好的实体模型久置室外以更接近于野外实地条件。实体模型见图1所示。

图1 梁峁坡面模型Fig.1 Liang-mao slope

本研究共进行 6次人工模拟降雨试验，每次试验前均对雨强进行率定，以确保降雨均匀度和降雨强度达到试验要求。为保证试验下垫面地表状况均一，试验前一天用30 mm/h降雨强度进行前期降雨，直至坡面出现产流为止停止降雨。次降雨历时依据细沟发育过程中的演变特征确定，即当坡面出现明显微地形变化时停止降雨。

1.3 研究方法

控制点布设影响重建坡面三维模型的精度，照片的采集关系到模型能否建成。控制点用于优化相机位置和方位数据，得到更好的重建效果。采集照片前，在坡面两侧边界处沿坡长方向每隔1 m布设标靶作为控制点，共布设24个控制点，并用全站仪精确测量每个控制点的三维坐标。利用数码相机进行照片采集，拍摄时尽量与坡面保持垂直，镜头与坡面距离约1.5 m左右，保持一个较小的俯仰角，减少建模误差。设置相机获取图像格式为TIFF无压缩照片，照片分辨率为6 000×4 000，采集照片所用镜头为定焦50 mm f/1.8式标准镜头，定焦镜头因为只需对一个焦段的成像进行纠正与优化，所以往往很少会出现畸变现象。

该研究使用Agisoft公司开发的PhotoScan软件实现坡面侵蚀过程的三维重建。该软件中集成运动恢复结构原理、多视图立体视觉原理以及生成不规则三角网DEM（digital elevation model）计算模型等，可以实现对数码照片建立稀疏点云模型、密集点云模型以及生成DEM等过程。由于建模是依靠跟踪照片中拍摄主体上的定位点来确定最终三维模型的点集，因此要求拍摄相邻两张照片的角度相差不能过大，照片重叠率在50%～60%以上，以防跟踪点丢失。照片拍摄位置的选取及轨迹、角度的选定根据试验环境和条件合理规划，如图2所示。

图2 拍摄位置及角度Fig.2 Photographing position and angle

将采集的照片先通过运动恢复结构原理生成稀疏点云，在这个阶段，PhotoScan会搜索照片上的公共点并将它们匹配，同时会找到每张照片的相机位置，并优化相机校准参数，结果形成一个稀疏点云模型和一组摄像机位置。再通过多视觉立体视图原理由稀疏点云模型基于估计的相机位置和图片本身构建密集点云模型。数据处理过程中将不参与计算以及存在模糊处无效点云数据进行剔除，以减少试验误差，建模效果如图 3所示。试验土槽面积大约36m2，生成稀疏点云的个数与密集点云的个数依据对坡面所拍摄照片的数目以及坡面的复杂程度有所差异，据统计，6次试验稀疏点云密度在 25 000～45 000 points/m2范围内，密集点云密度在25 000 000～35 000 000 points/m2范围内。PhotoScan生成密集点云模型后，可直接在软件内部计算获得 DEM 数据及正射影像，该试验中所获取DEM数据分辨率约为0.2 mm/pix，但为便于数据的处理，在 ArcGIS中将数据插值为0.01 m×0.01 m级分辨率数据进行分析。

图3 基于PhotoScan的三维重建Fig.3 3D reconstruction based on PhotoScan

坡面开始产流后每隔2 min接取径流泥沙样，之后将样品静置 24 h以上，倒掉上层清液，将剩余泥沙样在105 ℃烘箱内烘干并称其质量，以上所获取数据用于计算降雨各阶段径流量和含沙量；利用PhotoScan软件对每次试验坡面进行建模并计算坡面 DEM（分辨率为 0.01 m×0.01 m），将试验前后2次数字高程模型DEM相减，获取坡面高程变化量，将高程变化量导出到excel表格，继而分别与对应单元格面积相乘求出各单元格侵蚀体积，将所有单元格侵蚀体积求和即得总侵蚀体积，并结合土壤容重获取该时间段的计算侵蚀量。

1.4 侵蚀形态参数

梁峁坡面细沟侵蚀演变过程数字化、图形化及侵蚀形态参数的获取（图4、图5）。通过所获取数码照片在PhotoScan中三维建模获得坡面三维模型，提取梁峁坡面侵蚀横截面形态，以沟深≥1cm侵蚀深度作为确定侵蚀沟的标准，在三维模型以及实体模型中分别量取坡面细沟最大沟深（Hmax）、最大沟长（Lmax）以及平均沟宽（D）。通过此软件还可获取坡面 DEM，以此作为数据源，在ArcGIS中提取细沟沟网用于计算侵蚀形态参数。沟网的提取过程是基于水系特征提取的基本原理，ArcGIS中提取水系相关命令被集成在Hydrology分析模块中，数据处理时可直接调用相应函数进行计算，实现细沟沟网的提取与分析。

在此基础上计算细沟平面密度，计算方法如下：细沟平面密度（μ）：坡面所有细沟平面面积之和与坡面总面积的比值。

式中A0为坡面总面积，m2；Ai（i为自然数）为坡面每条细沟的平面面积，m2。

图4 细沟沟网提取Fig.4 Rill net extraction

图5 细沟平面提取Fig.5 Rill plane extraction

2 结果与讨论

2.1 坡面细沟侵蚀演变过程

图 6为利用 PhotoScan软件获得的 6场试验坡面DEM，其展现了连续6场间歇性降雨细沟形态的发育过程，其中不同高程用不同的颜色加以区分。随着降雨历时的推移，细沟侵蚀过程主要经历片蚀、细沟雏形、细沟沟网及小切沟4个阶段。

图6 不同侵蚀试验阶段坡面DEMFig.6 DEM of different erosion stages

降雨初期，雨滴的击打以及超渗产流所产生的薄层水流、股流相互作用，把表土中的细小颗粒带走，并在坡面留下鳞片状凹地，坡面主要发生溅蚀、片蚀，由图6a可看到在坡面3～9 m之间形成一系列串珠状跌坑，跌坑宽度多为5～9 cm，平均约为7 cm，深度多为1～4 cm，平均约为2 cm；分析图6b可知，第2次降雨后坡面片蚀的数量和面积不断增加。随着径流的汇集，径流侵蚀能力逐渐增强，在某一处径流侵蚀能力大于土壤抗蚀能力时，促使该处跌水演变为下切沟头，并有向上溯源侵蚀的趋势，这标志着细沟侵蚀的开始，梁峁坡面内主要侵蚀方式由面状侵蚀转变为细沟侵蚀，由面蚀所产生的微小跌坑在径流作用下长、宽、深均不断增大，最大分别达到266、7.6、13.8 cm；由图6c、6d可知，初始形成的细沟在降雨与径流的作用下沟尾下切及沟头溯源侵蚀，使细沟沿坡向延伸并向两侧拓宽，降雨历时 58 min（降雨量约90 mm左右），在原有细沟基础上坡面新增诸多小细沟，坡面上的断续细沟出现分叉、合并及互相连通现象，细沟网雏形基本显现，细沟向长发展较向宽发展迅速，测得最大沟长增加至376 cm，最深处达16.2 cm；图6e、6f为细沟沟网的进一步发育，最大沟深达36 cm，小切沟形成并停止降雨试验。从图6f中来看，梁峁坡面自坡底到坡顶上共发育出 4条明显的主细沟，均伴随着分叉、合并等现象。对于细沟侵蚀过程前人已经做出过大量研究，无论是何种微地貌，大都经历溅蚀、片蚀、细沟及细沟网阶段[22-24]，但是前人所做的研究仅仅是定性描述，本研究通过三维建模，图像化、数字化侵蚀过程，更能直观的展现土壤侵蚀过程，而且进一步对建模数据进行挖掘，提取侵蚀参数、计算侵蚀量等。

前期降雨试验后，在坡面顶部鲜有细沟出现，细沟主要集中于梁峁坡面中下部，这主要是由于坡面顶部坡度及汇水面积均较小，水流剪切力不足以剥离表层土壤；坡面中下段，径流量和流速均较顶部显著增加，使得产生细沟机率增大；同时，在径流的流程上如果不出现小跌水，细沟也很难形成，He等[25]研究了裸坡的细沟侵蚀发育过程，也认为坡面首先在溅蚀的作用下出现下跌的小水坑，且水坑主要沿坡面中部及下部平行分布，逐渐各个跌坑顺着坡面相互贯通，为细沟出现创造条件。

2.2 侵蚀形态参数的变化特征

引入细沟侵蚀形态参数细沟沟深、沟长、沟宽、细沟密度等，量化坡面细沟网发育特征，反映细沟侵蚀沿垂向以及水平方向上的发育强度。对整个坡面细沟侵蚀过程中不同时期的侵蚀形态参数进行统计，结果如表2。在前10 min的降雨，坡面主要发生溅蚀、片蚀，并未出现细沟，跌坎儿深约为 0.04 m，该阶段由以溅蚀为主迅速过渡到溅蚀、片蚀共存；10 min后即第2次降雨试验细沟雏形形成，表现为径流引起的部分片蚀区域贯通，最大沟深增加至0.05m，细沟沟长、沟宽、沟密度分别为1.51 m、0.058 6 m、1.100 m/m2；降雨历时20～40 min过程中，坡面出现多条断续细沟，最大沟深和最长沟长分别增加至0.138 3 m和2.66 m；第4场降雨试验后，降雨量达80 mm以上，细沟沟网阶段基本形成，沟网密度增大至1.999 m/m2，最大沟深为0.162 m；在之后的降雨试验中，最长沟长、最大沟深及沟网密度仍有不同程度增加，坡面愈加破碎，小切沟形成，最大沟长及最大沟深较细沟形成时增大 3倍以上。细沟宽度、沟深均随降雨历时的增加而增加，且细沟深度的增加速率大于细沟宽度的增加速率，细沟的宽深比均随降雨历时的增加而减小。覃超[26]做了基于立体摄影技术的黄土坡面细沟侵蚀发育过程量化研究，也得出相似结论，其将细沟分成主次级，得出细沟长度随降雨的变化同时受宽度和深度制约，主细沟和次级细沟除长度外深度和宽度变化相同。Zhang等[27]研究细沟侵蚀过程也表明细沟密度随降雨历时逐渐增加，其还研究了不同雨强作用下细沟密度的变化情况，表明在90 mm/h雨强作用下细沟密度增加幅度最为剧烈。

表2 细沟侵蚀形态参数Table 2 Rill erosion form parameter

2.3 基于三维重建技术-PhotoScan软件精度评价

利用PhotoScan软件间接计算侵蚀体积，并结合土壤容重获取该时间段的计算侵蚀量。表3展示了6次降雨试验期间坡面次降雨侵蚀量，将计算侵蚀量与实测侵蚀量进行对比，用以评价PhotoScan三维重建技术获取侵蚀量的精度。

表3 模拟值与实测值相对误差表Table 3 Relative error table between simulated and experimental values

将每场试验实测总侵蚀量与通过 PhotoScan建模间接计算的总侵蚀量对比，实测值与计算值之间的相对偏差见表3。试验值和模拟值的相对偏差在5.20%～20.82%范围内。第 1场降雨试验，坡面主要发生溅蚀与面蚀，地表松散、细小颗粒较多，随着径流发生，泥沙样中含沙量较高。然而此阶段基于三维重建技术-PhotoScan软件计算的侵蚀量是基于地形变化而得，坡面形态变化微小，故试验偏差较大导致误差较大（20.82%）。随着侵蚀过程的加剧，各次试验总侵蚀量迅速增加，侵蚀量计算值与实测值绝对偏差逐渐缩小，故其他 5次试验的相对偏差减小，均小于等于 10%。与肖海等[28]使用三维激光扫描仪对坡面侵蚀研究所得结果相比，PhotoScan软件三维重建模型精度总体要略低于三维激光扫描技术，李俊利等[17]利用照片重建技术生成坡面侵蚀沟三维模型得出通过重建模型间接计算所得侵蚀量与三维激光扫描仪相差8%左右，但其研究对象对特定的侵蚀沟，小区面积较小，因而外部环境因素所带来的影响也相对较小，基于照片重建和三维激光扫描方法计算的侵蚀量相对误差更为接近。

综合上述分析，实测值和计算值的数据比较接近，其偏差在合理范围内，PhotoScan软件对所获试验照片进行三维重建具有一定的可信度。同三维激光扫描技术相比，其也可展现一定的优越性。三维激光扫描技术能够测量精度为毫米级的物体的三维信息，尤其对小区尺度细沟侵蚀的观测更加试用。张姣等[29]应用三维激光扫描技术对试验土槽中的细沟侵蚀量进行了测量，在 15 m2试验小区上获取了 130 000个数据点，点云密度达8 667 points/m2，用时0.7 h。基于数码照片重建技术较三维激光扫描技术能够更加真实地反映实物的形状、明暗以及纹理等细节。李俊利等[17]应用照片重建技术对小区尺度上沟蚀演化过程三维建模并进行分析，其在VisualSFM 软件中完成特征点的检测与匹配以及稀疏点云模型，又使用CMVS 及PMVS2 工具进行稠密重建，最终在 Meshlab 中实现可视化，获取点云数据密度在4 800～6 200 points/m2范围内。本研究中，稀疏点云、密集点云以及DEM数据生成均在Agisoft PhotoScan中进行，只需简单设置相应的参数，无需掌握过多的计算机视觉相关理论就可以完成点云密度更加密集的三维模型。除小区尺度外，对通过无人机对某一区域进行照片采集，三维重建技术还可以完成对某一区域的数字测量，在土壤侵蚀监测和水土保持决策中应用前景也更加广泛。影响三维重建技术精度因素有多种，最根本的在于原始数据获取，获取优质量的数码照片可提高建模精度，对于野外试验，拍照时易受到阳光直射、阴影等影响，阴天试验更为适宜；在试验过程中，收集到的照片重叠率、分辨率越高，所得到的模型效果也就越好，但是相应的处理时间会更长，对计算机性能的要求也会提高，所以不建议使用分辨率过高的照片，照片的重叠率控制在50%～60%[30]即可。

3 结 论

本文以黄土高原丘陵沟壑区梁峁坡面为研究对象，利用人工模拟降雨试验，结合PhotoScan软件（三维重建技术）和ArcGIS软件，对试验坡面进行三维建模并提取相应的侵蚀形态参数，研究了坡面细沟侵蚀演变过程，同时将试验结果与同等条件下 ArcGIS软件模拟结果进行了对比，验证了PhotoScan软件的精度。主要得到以下结论：

1）利用PhotoScan软件对坡面侵蚀发育过程进行三维重建。其细沟形态发育过程表现为：溅蚀、面蚀、细沟雏形、细沟网形成及小切沟出现 4个阶段，每一阶段都有较为明显的特征。

2）分析降雨试验过程中各侵蚀形态参数的变化特征可知：细沟平面密度、最大沟深、最长沟长、平均沟宽分别由 1.1 m/m2、0.04 m、1.51 m、0.058 6 m 增长至2.991 m/m2、0.368 4 m、11.3 m、0.121 4 m，坡面细沟侵蚀不断发展并不断加强，定量地展现了坡面细沟侵蚀的演变过程。

3）利用人工模拟降雨试验所得实测总侵蚀量数据对PhotoScan软件三维建模间接计算的总侵蚀量进行验证，试验值和模拟值的相对误差在5.20%～20.82%内，除第1次试验相对误差较大（20.82%）外，其他 5次试验的相对误差均小于等于10%。利用PhotoScan软件对所获试验照片进行三维重建具有一定的可信度，可以作为土壤侵蚀研究的一种手段。