基于无人机倾斜摄影的白格堰塞区三维重建

2020-10-18包腾飞郑东健胡雨菡谈家诚

水利水电科技进展 2020年5期

朱征，包腾飞,3，郑东健，龚健，胡雨菡，谈家诚

(1.河海大学水文水资源与水利工程科学国家重点实验室,江苏南京 210098;2.河海大学水利水电学院,江苏南京 210098; 3.三峡大学水利与环境学院，湖北宜昌 443002)

我国西南地区滑坡、泥石流、崩塌等地质灾害多发，运动土石冲入河道就会截断河流，形成堰塞湖。仅2008—2017年，有记录的滑坡堰塞坝超过100座[1]。该地区历史上曾发生过多起大型堰塞堵江事件，如1786年大渡河[2]、1933年叠溪[3]、1967年唐古栋[4]、2000年易贡[5-6]、2008年唐家山[7]、2010年舟曲[8]、2014年红石岩[9]。堰塞湖形成迅速，牵涉范围广，具有很高的危险性。2018年10和11月发生的白格堰塞事件，位于四川省白玉县则巴村和西藏自治区江达县白格村交界处(31.08°N, 98.71°E)。白格堰塞湖是同一地点连续两次大规模高位滑坡阻断金沙江干流形成的堰塞湖[10]，造成上游多个乡被淹，下游多座桥梁被冲毁，四川、西藏及云南多地紧急转移安置数万人，引起了社会的广泛关注。

堰塞体往往结构松散、稳定性差，极易受外界环境影响发生形态变化，存在随时溃决的风险。堰塞体一旦失稳，其蓄积的大量水体将冲击沿线数百公里，造成的破坏难以估量。作为了解灾情和科学决策的重要依据，快速、准确地获取现场信息必不可少。然而，受灾地区一般伴有通信与电力中断、交通瘫痪、基础设施遭受重大破坏等紧急情况，堰塞灾害的事发地也存在交通不便、环境原始、人员物资输送困难等问题。这些因素制约了常规测绘手段的开展。

无人机倾斜摄影作为一种新兴的测绘技术，在灾害救援[11]、建筑[12]、城市规划[13]、测绘[14]、考古[15]、水利规划[16]等领域已经有许多应用。Kirnbauer等[17]利用倾斜摄影量化研究了高山雪线的变化。Doneus[18]使用航空影像进行倾斜摄影成图，从而还原了Teurnia古城遗址。Mergili等[19]以Piz Cengalo-Bondo滑坡事件的倾斜摄影重建成果为基础，对该灾害事件进行了后分析。Tonkin等[20]使用地面控制点提升了无人机倾斜摄影重建精度，并在一处野外边坡开展了验证试验。Harwin等[21]利用无人机采集海边悬崖的图像，并以重建点云实现了亚分米级侵蚀监测。Eltner等[22]在野外搭设多目相机系统，监测了雨水冲刷引起的土壤流失。在滑坡灾害研究方面，帅向华等[23]将倾斜摄影技术用于在云南鲁甸地震引发的滑坡测绘；许建华等[24]利用无人机倾斜摄影对地震裂缝、滑坡等几何信息进行量测，从而为烈度评估提供参考；许强等[25-26]将倾斜摄影用于白格堰塞堵江事件的灾害调查；王晓刚等[27]利用旋翼无人机倾斜摄影技术在茂县高位滑坡灾害进行了应急测绘。陈诚等[28]采用无人机进行了表面流场测量。

倾斜摄影在众多领域的运用中表现出该技术具有测绘效率高、人力资源占用较少、成果呈现方式多元化等优点，具备用于堰塞现场的非接触式应急测绘的潜力。

本文介绍了倾斜摄影的建模原理，阐述了无人机倾斜摄影的步骤要点，并在白格堰塞现场进行了无人机飞行作业。重建了白格堰塞现场的三维数字模型，并生成了数字高程模型(DEM)、数字地表模型(DSM)、数字正射图像(DOM)。无人机倾斜摄影三维重建技术在白格堰塞现场测绘方面表现出了良好的适用性。

1 倾斜摄影的建模原理

倾斜摄影是运用无人机云台摄像机，以一定倾斜角度获取地面物体的图像。通过设定拍摄间隔，保证图像间存在重叠。利用图像间的重叠部分进行空中三角测量计算，还原相机姿态并获得稀疏点云。再以稀疏点云为基础进行稠密重建得到稠密点云，从而实现倾斜摄影三维重建。

1.1 空中三角测量

空中三角测量(aerial triangulation)包含4个部分：图像特征检测、图像特征点匹配、点云初值计算、误差优化。

图像特征检测指的是在图像中找出具有一定代表性的像素点，并且用描述算子对该点进行描述。在众多算法中，SIFT(scale-invariant feature transform)算法[29]得到广泛应用。该方法对旋转、尺度缩放、亮度变化保持不变性，同时对视角变化、仿射变换、噪声也保持一定程度的稳定性。

特征点匹配将图像间的具有一定相似度的特征点进行配对。如图1所示，倾斜摄影图像间存在重叠，重叠区域内存在大量的相似特征点。相似度计算主要有两类方法，一类是按候选点描述算子间的空间距离按阈值过滤，另一类则是采用聚类算法寻找同类点。由于匹配时不可避免地存在一些错误匹配，这需要应用双重检验机制来提高准确性，如随机采样一致性(RNSAC)[30]算法。

图1 基于重叠图像的特征检测原理

点云初值计算是利用特征点间的配对关系进行的特征点三维坐标计算的过程。根据经典的小孔成像模型，像素点、实物点、焦点之间满足式(1)的关系。在此成像模型下，实物点坐标和像素点坐标之间呈线性关系。当以多视角拍摄物体时，相机焦点、实物点、像素点间存在一定关系，即外极几何约束。如图2所示，两张图像拍摄于不同位置，对应的焦点位置为F1、F2，其中，像素点A1、A2指向同一实物点A，那么在极线平面F1AF2上存在关系式(2)。

图2 基于两视图的实物点坐标计算原理

式中：K为内部参数矩阵；f为相机焦距；u、v为相机主点位置；R、t分别为相机成像坐标系与世界坐标间的旋转矩阵和平移矩阵；X、Y、Z为实物点坐标；x、y为像素点坐标。

(2)

式中：A1、A2分别为像素点A1、A2的坐标(x1、y1)和(x2、y2)；K1、K2分别为两相机的内部参数矩阵;E为本质矩阵。

式(2)中仅本质矩阵E中存有未知量，应用8组以上的特征点对即可求解E全部未知量。再通过SVD分解得到相机间相对姿态，以此代入式(1)可计算出每个特征点的三维坐标。匹配点组数越多，该过程的准确度越高。

实际成像过程与理论存在差异，如镜头畸变、投影模型非线性、特征点的检测与匹配错误。这些由各处积累的误差会降低实物点三维坐标计算的精准度。目前，得到广泛运用的思想是最小化整体重投影误差平方(式(3))。重投影误差平方和最小化的实现，普遍采用的是光束平差法(bundle adjustment)。这是一种启发式的阻尼高斯牛顿法，在几何视觉中广泛使用，具有较高的求解速度和稳健性。

(3)

式中：P为实物点向像素点的投影矩阵；Xi为实物点坐标；xi为像素点坐标。

1.2 稠密重建

空中三角测量通过对图像集特征点的三维投影计算，重建了特征点对应的实物点云。然而，特征点数量相对全部像素点而言是稀疏的，因此空中三角测量得到的点云也是天然稀疏的，包含信息丰富度较低。

多视影像密集匹配(multiple view stereo)以稀疏点云为基础，以一定方法寻找具有一致性的像素点加密点云，从而增强点云的可利用度。点云扩散法是较为常用的稠密重建方法，主要思想是以初始点云为基础按照邻近关系逐步填充未重建像素。在此过程中会产生部分飞离点，以深度值与邻近区域一致程度可以剔除这些飞离点。

经过多年的发展，倾斜摄影及三维重建技术的理论水平和应用能力随着研究的深入不断提升。摄像器材、数字成像技术、三维重建算法和计算机性能的发展使得三维重建更加容易实现，其成图精度更高、速度更快、成图规模更大。

2 倾斜摄影三维重建的基本步骤

2.1 原始数据采集

原始数据指倾斜摄影作业采集到的图像数据、POS(position and orientation system)数据(位置、姿态)和相机参数。图像由机载相机拍摄，POS数据由机载传感器记录，相机参数(径向畸变、成像中心点等)由厂商标定。

图像数据使用无人机机载相机获取，采集的图像应当清晰、重叠度高、视角丰富。倾斜摄影作业前应拟定好飞行线路、摄像重叠度、拍摄数量等参数，选择光照充足、少云正光的时段进行飞行作业，避免云朵移动造成被摄物体阴阳变化及斜光照射产生的大背光区。

2.2 数据筛选

采集到图像后，有必要对图像进行筛查。一方面，检查采集成果是否达到计划的数量、范围，重点关注地形起伏、细节复杂、光影交错的区域，若不满足要求，应重飞补拍；另一方面，筛查图像质量，防止低质量图像增大重建误差，应舍弃失焦、杂物遮挡的图像，有杂点、移动物的图像可作蒙盖处理。

2.3 重建计算

图像集经过预处理后，首先进行空中三角测量计算。空中三角测量通过图像间的像素特征匹配生成大量连接点，利用连接点的匹配信息计算获得连接点的三维位置(稀疏点云)。POS数据和图像畸变的误差会随迭代得到修正。然后，采用多视影像密集匹配技术，利用像素网格匹配对稀疏点云进行加密计算，生成稠密点云。由于稠密点云数据量极大，通常采用切块降低计算机的压力。切块也有利于计算机群的平行处理。随后，稠密点云经过优化处理降低数据冗余，通过三角网格化最终获得不规则三角网模型。利用不规则三角网模型与像素间的匹配关系，将原始图像纹理贴附到不规则三角网模型上，可得到带纹理的三维真实感模型。

2.4 成果生成

倾斜摄影三维重建的基本成果是三维网格模型(3D mesh)，三维网格模型由不规则三角网构成。对三维网格模型进行光滑处理可生成不同精度的数字高程模型(DEM)；三维网格模型经材质附着可得到数字地表模型(DSM)；数字地表模型经正射投影可生成数字正射图像(DOM)。

这些成果可以被运用到广泛的领域，数字高程模型可用于地理信息系统(GIS)、正射投影影像可用于测绘成图(Mapping)、数字地表模型可用于建筑信息模型(BIM)。

3 研究区域和数据获取

3.1 白格堰塞湖形成过程

金沙江是长江上的重要支流，位于青藏高原东缘。金沙江落差大、流速快，两岸边坡陡峻，周围地形沟壑交错，山顶与谷底的落差可达1 km。频繁的地质活动和强烈的风化作用使沿线的滑坡、崩塌事件非常多，由此引发的堰塞灾害是沿线居民生产生活安全的重大威胁。图3为卫星记录的白格堰塞的初期征兆、堰塞规模、溃后状态等光学影像信息。

图3 白格堰塞事件全阶段卫星光学图像

2018年10月11日，第1次滑坡发生，滑坡体积2 400万m2、堆积高度85 m；1天后，堰塞湖蓄水达2.9亿m3，堰塞体自溃留下左岸残积体[31-32]；11月3日，第1次滑坡区域原位发生第2次滑坡，下滑土石350万m2，冲压第1次滑坡残积体形成135 m高的堰塞体，数天内水位快速上涨，抢险人员采取人工开挖导流槽的方式主动降低坝体高度；11月12日堰塞湖蓄水达5.24亿m3，水位达到人工导流槽高程2 956.40 m并开始泄流，3天后金沙江水位基本稳定，险情解除。这起连续滑坡堰塞事件造成上游江达县波罗乡、白玉县金沙乡等先后被淹；导致下游318国道金沙江大桥损毁严重，云南迪庆藏族自治州、丽江市部分道路中断、学校被淹、农田被冲毁。

白格堰塞事件发生突然、滑坡体量大、周边环境原始，造成的经济损失巨大、威胁人口众多，具备大型堰塞灾害的典型特征，对研究无人机倾斜摄影技术在堰塞现场的适用性有很高的参考价值。

3.2 数据获取

采用无人机对堰塞体展开了倾斜摄影作业，使用的无人机型号是DJI Phantom 4 RTK，该设备质量1.3 kg、单次续航能力约30 min，综合GPS、北斗和伽利略全球卫星导航系统(global navigation satellite system，GNSS)，机载相机型号为DJI FC6310R(1 英寸 CMOS；5 472像素×3 648像素分辨率；3轴云台稳定系统)。

如图4所示，无人机飞行两个架次分别重点采集右岸泄流槽和左岸残积体数据。为保证三维数据生产成果质量,航线采用大重叠率:航向重叠率为80%,旁向重叠率为70%。这两架次飞行共获取1 497幅JPEG格式图像(29.9 G像素)及每张图像的位置信息、姿态信息。如图5所示(图中箭头指向为水流方向)，图像涵盖了堰塞区的重点部位，如泄流槽右岸边坡、河谷左岸残积体、河谷右岸滑坡面等。

图4 白格堰塞体无人机倾斜摄影作业

图5 白格堰塞区现场无人机航拍影像

4 基于图像的三维重建处理

三维重建需要图像数据、POS数据(位置、姿态)、相机参数。完成现场采集后，导出无人机内置GNSS中记录的位置数据、云台稳定器中记录的姿态数据以及相机拍摄的图像数据。

无人机GNSS记录的位置数据存在一定误差，这部分误差可在空中三角测量计算过程中由光束平差法迭代优化，调整后的位置数据精度得到提高。姿态信息(俯仰角、横滚角、偏航角)由云台稳定器记录，云台稳定器的硬件性能保证了角度抖动量低于0.02°，具备高精度。相机参数在无人机出厂时经过校正，其径向畸变、横向畸变、焦距、成像中心点等参数较准确。

图像经过初步处理后，进行空中三角测量计算。以GNSS定位最大误差为调整限值对位置数据采取平差调整策略，姿态数据采用云台稳定器记录数据原值，相机参数采用了出厂校正值。如图6所示，空中三角测量计算总计生成442 947个连接点(平均每张图像包含1 223个连接点)，覆盖面积为417万m2，像素点的地面分辨率平均值为5.79 cm，重建均方根误差(RMS)为0.62个像素点。

图6 白格堰塞区倾斜摄影三维重建

完成空中三角测量计算后，进行三维场景重建。白格堰塞区的场景范围较大，以堰塞体为建模中心(31.083 8°N，98.715°E)设定建模区域(1 796 m×2 322 m×573 m)。三维场景重建对计算机内存的需求较高(310 GB)，因此采用自适应切块方法。自适应切块是以稀疏点云密度分度预估稠密点云的密度分布，将高密区域分成小瓦片，从而保证单个瓦片的内存需求在硬件容量内。三维重建成果如图7所示。其中数字高程模型描述的是堰塞区地面高程信息，可用于堰塞体覆盖面积、体积、堆积高度、坡度等特征的计算，作为现场作业规划的直接依据；数字地表模型是在数字高程模型的基础上，进一步涵盖了除地面以外的其他地表信息的高程，丰富的信息用于判别堰塞体的结构分布、材料类型、颗粒级配等，可为堰塞体稳定评估提供参考。

图7 白格堰塞区倾斜摄影建模三维成果

本节三维场景重建处理所采用的工具为三维图像建模软件ContextCapture和PhotoScan。

5 成果应用

5.1 白格堰塞区建模成图

无人机倾斜摄影获取现场图像数据、POS数据，经过空中三角测量计算、三维场景重建，获得堰塞区三维网格模型。以此为基础可生成白格堰塞区的数字高程模型(DEM)、数字地表模型(DSM)、数字正射图像(DOM)等形式的成果如图8所示。数字正射图像是利用数字高程模型对无人机图像进行投影差改正、镶嵌、剪裁生成的数字正射影像数据集，具有精度高、信息丰富、直观真实等特点，可作为减灾规划决策的参考依据(图8(a))。

图8 白格堰塞区倾斜摄影建模平面成果

5.2 堰塞体形态测算

结合倾斜摄影重建成果，在数字三维模型中进行坐标读取、距离量测、坡比计算等工作，获得了堰塞体形态、结构等方面的信息，为救灾减灾提供服务。

人工泄流槽初始尺寸为：坎底高程2 952.52 m、长220 m、顶宽42 m、底宽3 m、最大深度15 m[32]。经测量，冲刷扩大后泄流槽左岸边坡角度为35°～38°、右岸边坡角度为32°～36°，与初始边坡角度35°基本一致；槽底高程约为2 865.50 m，较初始降低85 m；槽底宽度约为80 m，约为初始底宽的26倍；槽底长度约720 m，约为初始底长的3.5倍；槽顶宽度约为300 m，约为初始顶宽的7倍。

堆积体最高点位于左岸边坡(高程2 974.93 m)，距泄流槽左坡顶距离约160 m。以最大纵剖面为基准，距上游坡脚水平距离约230 m(23°～25°)，距下游坡脚水平距离约500 m(11°～14°)。