面向东川复杂山地泥石流沟谷三维地形建模及特征分析的无人机遥感探测应用研究

2021-07-02李绕波

中国地质灾害与防治学报 2021年3期

毕瑞，甘淑,2，李绕波，胡琳

（1.昆明理工大学国土资源工程学院，云南昆明 650093；2.云南省高校高原山地空间信息测绘技术应用工程研究中心，云南昆明 650093）

0 引言

我国云南地区海拔高、地形地势复杂、自然环境多变，采用GNSS[1]、光学遥感卫星[2]、雷达遥感[3]、地面三维激光扫描等传统监测方法能实现对灾害区的监测，但是存在以下主要问题[4−6]：①卫星过境周期较长，无法实现短时间周期监测；②受自然环境因素干扰较大；③影像数据分辨率较低，对小型、单一灾害区域无法进行精确灾害监测。因此，对于灾害监测来说，高效、精确地获取灾害区的数据，才能准确地对灾害特征进行定性、定量的分析。

目前，无人机遥感技术应用于地质灾害监测成为了国内外学者的热点研究方向。无人机遥感技术具有低空作业、成本低、机动性强、多维获取高精度影像数据等特点[7−9]。众多学者利用无人机遥感技术实现了高效、精确获取灾害区不同灾害类型的影像数据。结合地面控制点，影像数据进行绝对定向后，获取高精度的数字地表模型（DSM）、数字正射影像（DOM）、数字高程模型（DEM）和高密集匹配点云数据[10−12]。通过对无人机影像数据成果的分析，获取灾害地区的特征变化信息，实现准确地地质灾害监测。黄海峰等[13]针对2014年9月初三峡库区的一处滑坡，生成地面分辨率为4.25 cm的DOM 和三维模型，分析得到滑坡整体下移已处于欠稳定状态。王庆国等[14]面向云南澜沧江某水电站滑坡监测，提出了将地面激光点云数据和航空影像相结合的滑坡监测方法，通过点云数据配准，实现地面激光点云数据和航空影像数据的融合统一，基于地表模型的整体变形分析和基于剖面的采样分析，对滑坡区域的变形规律进行研究。SOTIRIS 等[15]利用无人机获取2015年希腊莱夫卡达奥基诺斯滑坡影像，可以更早地对滑坡变化进行评估。DARREN 等[16]使用无人机收集了一个时间序列的高分辨率图像，建立滑坡地表数字模型（DSM），在时间序列上测量了滑坡特定区域的体积变化，用于滑坡动力学评估。总的来说，无人机遥感技术用于地质灾害的监测日益广泛，众多研究表明无人机遥感技术能够快速、高效、精确地实现对灾害区域的评估。但是，针对泥石流沟谷地区的无人机遥感技术应用研究较少，针对数据采集、数据处理和成果分析等应用方面的研究还有待加强。

本文以云南省东川区小江流域中游大白泥河沟谷为研究对象，获取该区域一段泥石流沟谷的无人机遥感数据。首先，对影像数据进行相对定向，获取具有相对位置的点云数据；其次，与同期地面三维激光扫描数据结合，选取明显的点作为控制点，实现无人机影像数据的绝对定向；然后，结合Smart3D 影像数据处理，获取绝对定向后的三维模型、DOM、DSM 和密集匹配点云，通过对密集匹配点云处理，得到DEM 数据；最后，利用得到的数据成果，结合ArcGIS 和Cloud Compare 中的相关分析模块，实现对该段泥石流沟谷地区的地形特征分析。针对提出的结合地面三维激光扫描数据实现影像绝对定向的技术和相关地形特征分析方法进行探索。

1 研究方法

1.1 研究区概括

东川区小江为云南省东北部金沙江右岸支流，起源于寻甸回族彝族自治县，全长141.9 km，流域面积约为344.3 km2，自然落差约2 860 m，多年平均流量约为51 km3/s，泥石流沟分布密度为37.8 条/1 000 km2[17]。流域内遍布泥石流冲沟，植被覆盖度低，泥石流灾害频发，是典型的暴雨泥石流区。

大白泥河位于小江流域的中下游地区，流域面积18.1 km2，主沟长11.8 km，沟床纵比降10.66%[18]。流域内地势险峻，沟谷两侧分布陡峭高山，沟谷地区滑坡体、崩塌体分布较多。雨季水量丰富，是一条具有代表性的暴雨性泥石流沟谷。本文研究区选取大白泥河沟谷地区的一段区域，具体研究区位置见图1。

图1 研究区位置Fig.1 Location of study area

1.2 技术路线

针对选取的泥石流沟谷地区实验区，选择泥石流频发的夏季，实现对泥石流沟谷地区的监测应用分析，研究中采用的主要技术路线见图2。

图2 技术路线Fig.2 Technical route

首先，针对该地形、地势起伏较大的复杂地区，通过实地踏勘，利用大疆DJIGSPro 航线规划软件进行航线规划，获得清晰、完整的泥石流沟谷影像数据；其次，结合无人机自带的POS 数据，利用Smart3D 软件建立具有POS 坐标系的无人机影像点云数据，结合同一时期地面三维激光扫描仪获取的具有准确坐标、完整的点云数据成果，利用Cloud Compare 软件，对两点云数据进行粗配准，通过目视判读，选取三维激光点云中的具有明显特征的点位作为无人机影像的地面控制点，实现对无人机影像数据的绝对定向；之后，利用Smart3D对绝对定向后的无人机数据处理，生成具有颜色信息的三维模型以及相关数据成果；然后，采用PhotoScan 中不规则三角网渐进加密点云滤波，去除无人机点云中植被、建筑物等非地面点，获得泥石流沟谷的地面点，在ArcGIS 中，采用反距离权重(IDW)插值方法，生成0.5m 分辨率的DEM 数据；最后，结合三维模型、DOM和DEM 对泥石流沟谷地区进行可视化表达，实现对研究区定性、定量的地形特征分析。

2 面向研究区的无人机倾斜摄影影像数据采集与空间定向处理

2.1 无人机数据采集

采用大疆Phantom 4 Pro 四旋固定翼无人机进行数据采集，无人机的基本参数见表1。考虑到大白泥河沟谷地区地势复杂，陡峭高山分布于中心沟谷两侧，中心沟谷最低点海拔约1 380 m，两侧高山最高点海拔约2 265m，高差约885m，为了确保飞行安全，航线规划中避免过度靠近山体，特将无人机飞行高度设置为300m，地面分辨率为7.0 cm/px。按照《CH/Z 3005 低空数字航空摄影测量规范》要求，航向和旁向重叠度均设置为70%，保证影像数据重叠度的完整。

表1 无人机基本参数Table 1 Basic parameters of UAV

此外，基于大疆平台下的DJIGSPro 软件实现对研究区的航线规划，采用单镜头下设置不同角度而采集地面数据的倾斜摄影模式，即规划飞行3 架次，且对各架次分别获取正视、左视和右视影像数据，共计获得影像274 张。具体飞行参数和航线布设见表2、图3。

表2 飞行参数Tab le 2 Flight param eters

图3 无人机航线布设Fig.3 UAV routes are laid out

2.2 基于地面激光三维扫描的无人机影像定向处理

无人机获取的影像数据中包含了每一张照片的经度、纬度、海拔和RGB 颜色信息，但无人机搭载的是GPS/GLONASS 双模式定位系统，获取POS 数据在平面方向上能满足单点测量的误差精度要求，但在高程方向上误差较大，使得整体点位精度较低。单纯利用POS 数据，能实现对无人机影像进行相对定向，但无法进行绝对定向。

泥石流沟谷地区具有突发性，对于范围较大的区域无法实现全面的地面控制点布设，实现影像的绝对定向。因此，提出结合已有的同期地面三维激光扫描仪获取的泥石流沟谷地区点云数据处理成果实现影像的绝对定向。为了后续得到精度较高的无人机影像数据成果，将POS 数据下无人机点云与高精度三维激光扫描仪获取的点云进行粗配准、利用明显的地物特征点云坐标值实现无人机影像的绝对定位提高影像数据成果精度。

2.2.1 无人机点云数据粗配准定向处理

地面采用MAPTEK I-Site8200 地面三维激光扫描仪采集点云数据。利用仪器内置的GPS 和罗盘采集测站点的坐标和方位值，实现各测站点的点云坐标统一（WGS84，UTM 48N）。本文中使用的是经过内业处理后得到的地面点云数据成果。

利用Smart3D 专业无人机影像数据处理软件，对无人机影像进行空中三角测量，由于未采集地面控制点，基于自身坐标系（WGS84）以及采集的POS 数据，实现对影像的相对定向。采用运动恢复结构算法（Structure from Motion,SfM），通过相机的移动来确定目标的空间和几何关系[19]结合多视角立体算法（MultiView Stereo,MVS）,得到相对定向后的无人机高密集匹配点云数据。

无人机获取的点云数据具有位置信息和RGB 颜色信息，地面三维激光点云数据有位置信息，没有RGB 颜色信息，单独利用地面三维激光点云数据选取明显地物特征点存在一定的困难。并且，未配准的点云数据，存在水平和垂直方向上的偏差，垂直方向偏差大于水平方向。因此，将二者点云数据加载至Cloud Compare 软件中，通过软件中粗配准工具，以地面三维激光采集的点云数据为参考点云，将无人机点云数据与地面三维激光采集的点云数据进行配准。选取3 对以上明显的地物点进行粗配准，利用透视变换向量，计算一定的旋转和平移变换矩阵，将两组点云数据统一到同一坐标系下。粗配准精度并不能满足后续成果的精度，但可以减少后续在选取地面三维激光采集的点云提取坐标值的工作量，利于快速选取坐标点。两个点云数据粗配准处理前后结果对比见图4。

图4 点云数据粗配准Fig.4 Coarse registration of point cloud data

2.2.2 无人机影像精准匹配定向处理

通过粗配准后，从地面三维激光采集的点云上选取分布于沟谷地区的6 个明显地物特征点作为无人机影像绝对定向的地面控制点，控制点以“XK+点号”进行命名（图5）。在Smart3D 进行控制点刺点，实现影像的绝对定向，获取与地面三维激光采集的点云有相同比例尺且坐标系（WGS84 坐标系）一致的高精度数据成果。

图5 控制点布设Fig.5 Control points arrangement

3 面向泥石流沟谷的地形建模及其特征分析

3.1 研究对象三维建模及DEM 处理

选取地面控制点对无人机影像进行绝对定向后，基于无人机影像，通过SfM 算法，获得精度较高的密集匹配点云，基于点云数据，得到地表三维模型、DOM、DSM和无人机点云数据。相关成果见图6。

图6 无人机影像数据成果Fig.6 UAV images data results

基于无人机影像，生成的DSM 数据，其中不仅包含了地面点，而且也包含了地物、植被、建筑等信息，利用DSM 进行地形分析会存在较大的误差，因此，需要对无人机点云数据进行去噪处理，去除非地面点数据，通过插值得到只含有地面点的DEM 数据。

结合研究区的实地情况，可以看到，对于泥石流沟谷地区，少量的建筑物和分布于沟谷和两侧山坡的植被是最主要的非地面点数据。由于该地区地形复杂，山体坡度较大，因此，尝试利用PhotoScan 软件中的不规则三角网渐进加密的方法去除非地面点数据。通过设置像元尺寸大小，对整个点云区域进行划分，搜索每个划分区域中的高程最小值，视为初始地面点，再设置最大角度值、最大距离值，通过对整个点云区域迭代计算，小于设定阈值则作为地面点。去噪后得到的点云数据见图7。

图7 无人机点云去噪处理Fig.7 UAV point cloud denoising

得到地面点云数据后，利用ArcGIS 10.6 中的“LAS 数据集转栅格”工具，采用反距离权重插值方法，构建0.5m 分辨率的泥石流沟谷地区DEM 数据。研究区最小高程值为1 352.92m、最大高程值为1 655.54m，相对高差302.62m。由于采用地面三维激光点云数据作为无人机影像的地面控制点，选取一定数量的点位，得到高程值，对比有无控制点生成DEM 的精度，无人机DEM 成果及精度对比见图8。

图8 无人机DEM 及精度对比Fig.8 UAV DEM and accuracy compariso

3.2 泥石流沟谷灾害地形特征分析

泥石流的发生通常受众多自然因素和人为因素的影响，本文针对地形地貌要素进行分析。首先，通过DOM 影像直观的对泥石流沟谷地区中的灾体类型进行识别，获取灾体边界线；其次，利用DEM 数据进行坡度、坡向分析、剖面线分析，对泥石流沟谷地区的地形进行分析。

3.2.1 泥石流沟灾害体识别

结合DOM 影像，目视解译得到，研究区主要分布白色状沟谷，沟谷两侧分布滑坡体和小型泥石流冲沟，结合三维模型充分进行识别（图9）。统计分析得到，3 个滑坡体的边界长为251.12，334.07，421.17m；整个泥石流沟谷边界长为3 351.71 m；小型冲沟边界长为1 070.50m。

图9 灾害体目视解译Fig.9 Disaster bodies visual interpretation

3.2.2 泥石流沟坡度、坡向分析

通过DEM 数据分析得到，研究区相对高差302.62m，高差较大，对泥石流的发生、松散的土石体的搬运提供了位移能量。坡度也是泥石流的发生重要的因素之一，地形陡峻，且多以碎石、土石体分布，受暴雨冲刷，极易发生滑坡和崩塌地质灾害，为泥石流的发生提供了充足的固体物质。通过对研究区进行坡度分析，有利于分析不同坡度上固体物质的分布情况。对DEM 进行坡度分析、重分类后，坡度分析结果见图10（a）。将坡度（α）结果重分为7 类，统计得到：0°<α≤5°总面积为97 161.00m2，占研究区总面积9.72%；5°<α≤8°总面积为138 110.75m2，占研究区总面积13.81%；8°<α≤12°总面积为133 275.50m2，占研究区总面积13.33%；12°<α≤25°总面积为225 721.75m2，占研究区总面积22.57%；25°<α≤35°总面积为132 558.00 m2，占研究区总面积13.26%；35°<α≤55°总面积为223 236.00m2，占研究区总面积22.33%；α>55°总面积为49 830.75m2，占研究区总面积4.98%。

坡向是决定地表接收阳光和分配太阳辐射能量的重要地形因子，太阳辐射能量的分布影响泥石流沟地区植被的覆盖程度，从而引起泥石流灾害的发生。通过坡向分析，将坡向分为10 类，结果见图10（b）。分析得到，该泥石流沟谷地区大部分坡向为东、东北、东南和西北方向，结合DOM 分析，大部分滑坡体位于东北方向，是常发生泥石流的地区。

图10 坡度、坡向分析图Fig.10 Slope slope direction analysis chart

3.2.3 泥石流沟剖面线分析

结合剖面线分析泥石流沟谷横向和纵向的高程变化和地形起伏变化特征。获取泥石流沟谷横向和纵向的剖面线，横向剖面线为泥石流主沟西南方向至东北方向A（a→b），纵向剖面线为泥石流沟谷中部区域B（c→d），剖面线见图11。分析得到，泥石流沟谷主沟横向高程值逐渐降低，流通区高程起伏变化较大，外部堆积区高程值变化较小。泥石流沟谷纵向剖面线地形起伏较大，剖面线两侧高程变化较大，结合三维模型分析，两侧均为山体分布，随着距离的增加，沟谷地区的高程逐渐增加，350～450m 距离内，出现两个峰值，主要为泥石流固体物质流动冲击形成的特征地形。