基于无人机的滑坡地形快速重建与稳定性分析

张 欢，巨能攀，陆 渊，万 勋，蹇志权

（1.成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059；2.四川公路桥梁建设集团有限公司勘察设计分公司，四川 成都 610041）

斜坡的调查与稳定性评价一直是地质灾害研究领域的重点，传统工程地质分析方法更多地进行定性分析，所获得的结果与研究人员的业务水平和经验密切相关[1]；而物理试验方法耗时长且难以达到理想效果。采用轻型无人机航拍数据精细建模，再通过数值模拟进行定量分析，最后与工程地质分析方法相互验证，是一条高效、实用、可靠的途径。

地质灾害调查是地质灾害应急处置的首要和基础环节[2]。为了促进滑坡调查评价从经验走向科学，从感性判断走向理性量化，进而有效地为滑坡应急处置提供科学依据，必须快速查明其地质条件、潜在破坏面的形态位置及变形机制等特征[3]。传统地质调查手段需要技术人员深入灾害隐患区，耗费大量的人工和时间成本。卫星遥感对地观测精度、时效和频度上存在不足，同时噪点也造成很大的信息损失和干扰[4]。

近年来，无人机航摄系统以其机动灵活、响应迅速、成本低廉、适应性强、适用于高危地区探测等优点，被越来越广泛地应用于地质灾害调查之中[5]。无人机航空摄影测量技术可以通过携带CMOS（互补金属氧化物半导体）或CCD（电荷耦合器件）感光传感器的数码相机获取地质体不同方向的多幅二维数字影像，利用实际空间坐标系和数字影像平面坐标之间进行透视变换，通过计算机立体视觉技术，匹配影像中的大量同名点，以此进行空中三角测量并获得数码相机的空间参数，然后通过相关算法生成三维点云数据，最终生成三维网格模型[6−9]。传统的建模方法数据采集效率低、生产成本高、建模耗时长、建模精度不理想，无法满足地质灾害应急调查需求[10]。李永树等[11]的研究表明：与卫星遥感数据相比，无人机航摄影像形成的三维可视化模型能提供更详细、更丰富的几何和语义信息。王国洲等[12]通过无人机航摄数据的校正与拼接，生成DEM（数字高程模型）、DOM（正射影像图）等可视化三维数据，为地质灾害的调查评估及快速处理提供数据支撑。赵星涛等[13]基于ArcGIS软件将矿区下沉等值线与正射影像叠加，对矿山灾害影响区内的建（构）筑物和主要工程设施进行沉降分析。目前无人机数据采集、实景建模等技术日趋成熟，但是缺少利用无人机系统快速调查评价滑坡的方法体系。

本文将无人机技术应用于册亨平庆组滑坡调查过程，使用航摄数据建立三维有限差分模型进行数值分析。通过传统工程地质调查、无人机航片解译和数值模拟三者结合，快速、准确地调查评价滑坡的发育状况，系统地总结出一套以轻型无人机建模进行定量评价为核心的滑坡应急调查、评价方法。

1 无人机遥感系统

无人机系统搭载的高分辨率数码相机，不仅能够垂直拍摄获取正射影像，还能获取超低空多角度、高精度影像，满足三维建模的需求。无人机遥感系统是一个具有高智能化、作业能力强和低成本的航摄平台。系统由空中、地面和数据处理三部分组成（图1）。地面部分包括地面监视系统、地面控制系统以及数据传输显示系统。空中部分的主要功能是航线控制、姿态控制、数据采集和数据传输。飞行控制系统是空中部分的核心，该系统通过主控单元将GPS 模块、舵机、IMU（惯性测量单元）和遥控接收器等设备接入，从而实现无人机自主飞行功能。数据处理部分包括影像数据预览系统、影像数据处理系统，作用是对影像数据进行加工，以提取有效信息。

图1 无人机遥感系统Fig.1 Unmanned aerial vehicle remote sensing

2 无人机航摄建模流程

2.1 无人机航摄数据获取

影像数据的获取主要包括资料收集、定点勘查、参数设置和无人机航拍等流程。先规划出测区范围，然后对特定飞行作业区域进行参数设置和航线规划。航线规划是在航飞前按照应用要求、飞行作业区特点、飞行器和任务设备性能参数等因素，规划出飞行区域的航线与拍摄点。为了保证模型的细节以及数据的精确性，规划航线时需注意尽量提高相片重叠度[14−15]，调查区内各处至少需要有三幅不同位置的影像的重叠部分。为了确保现场航摄数据获取充分，应对采集数据进行低精度快速检测[16]，Pix4Dmapper（航测数据处理系统）提供了相应的功能。

2.2 数据预处理

影像数据处理包括影像的预处理、影像匹配、自动空中三角测量等技术。由于透镜固有的特性（凹透镜发散、凸透镜汇聚）无人机搭载的相机拍摄的影像存在一定程度的变形和失真，需要根据相机的畸变参数对影像进行畸变校正[17]。此外，部分影像由于天气差异和曝光程度造成了光线、色泽上不一致，为了提高影像匹配准确度并使拼接后的DOM 影像色彩自然均匀，还需对影像进行羽化和重曝光处理[18−19]。

2.3 三维地质模型建立

初步处理后的航摄照片通过Pix4Dmapper 生成加密点云文件和正射影像。然后通过Global mapper（地图绘制工具）滤波处理点云文件，将地面点数据与非地面点数据区分开并使用Polyworks（3D 测量软件平台）进行抽稀。最后将处理后的点云文件以及地质信息通过Rhinoceros（犀牛软件）结合建立出三维地质模型。

将无人机获取的含经纬度及高程信息的高精度正射影像导入Pix4Dmapper 进行处理。处理过程包括：（1）建立英文路径的工作目录。（2）新建英文名称和路径的项目，选择航拍项目。（3）选择图像设置图片属性，照片分区分组处理以降低单次运算量和对电脑配置的要求，区域处理完成后进行整合[20]。点云坐标的解算需要相机的内部参数以及相机校准参数，软件自动识别坐标体系和相机的型号参数，以及照片内的GPS 高程和经纬度数据。若采用该软件校正畸变，需要手动输入相机畸变参数。（4）本地一键化处理，先通过Pix4D 的AAT（自动空三测量）和BBA（平差数据）来计算原始图像的真实位置及参数，然后在空中密云中增加3D 密云点。输出的所有点都包含地理坐标信息，最终可以得到高分辨率的DSM（数字表面模型）和DOM 图像（.tif）[21−22]和点云文件（.las）。

由于点云中包含大量非地面点，接下来需要将点云数据通过Global mapper 滤波处理将地面点数据与非地面点数据区分开，得到反映真实坡体的空间模型即DEM[23]。处理过程如下：（1）加载点云（.las），将图形裁剪为矩形区域。（2）选择“自动分类地面点”筛选地面点。软件采用形态学滤波方法，需要根据现场地形条件、建筑和高植被特征设置最大高差和预期地形坡度等几何参数。反复修改参数，直到滤出需要的地面点，每次处理必须勾选“启动时重置现有的地面点为未分类”。（3）选择“筛选雷达数据”勾选“ground”分离出地面点。（4）使用3D 视图和RGB 附色视图检视处理成果。（5）选择“导出矢量/雷达格式”输出为txt文件。

将筛选后的数据导入Polyworks 并按照“工具”→“数据对象”→“采样”→“应用”对数据进行不断抽稀，降低内存的同时保留高精度的地形数据，抽稀程度根据计算机的运算能力和需求的精度确定，抽稀后输出为点云数据。

最后通过Rhinoceros 软件逆向建模功能建立坡体三维计算模型。流程如下：（1）导入点云数据（.txt），使用嵌面工具构建三维坡面。（2）执行“_MeshPatch”建立三角形网格。（3）执行“_Drape”建立模型表面。（4）曲面挤出成体，采用“布朗分割”裁减掉多余部分，分割后使用“平面洞加盖”封闭模型底面。（5）根据前期勘探资料或现场实测剖面获取的地层信息构建地质结构。（6）使用“grid”插件输出为FLAC3D三维计算模型。

2.4 成果解译

地质灾害在遥感影像上呈现的形态、色调、阴影、纹理等均与周围背景存在一定的区别[24]。因此，采集的遥感影像通过Smart3D 建立实景模型并解译，可以对目标区域内的地质灾害点和地质灾害隐患点进行系统全面的调查，查明其空间形态及细部变化特征。针对解译结果进行野外复核，验证解译结果的准确性，完善各灾害点的属性信息，建立地质灾害信息数据库，方便地质灾害的动态监测和实时更新，为灾害治理提供依据。

3 应用研究实例

本文通过对多个灾害点研究的总结，以及对册亨平庆组滑坡调查的思考，梳理出基于无人机航摄的滑坡应急处置流程（图2）。

图2 基于无人机航摄的滑坡应急处置流程图Fig.2 Flow chart of landslide emergency disposal based on UAV aerial photography

流程的核心是无人机航摄获取的高精度数据，然后通过Smart3D Capture 构建三维模型解译出坡体的裂缝发育和整体变形情况。同时建立三维有限差分模型进行数值模拟，得到的位移应力云图可以有效地指导边坡支护结构的布设形式、布设位置和支护参数设置，模拟出的变形破坏过程可以增强对边坡破坏机理的认识。由于滑坡的危险性，传统手段一般难以获得灾害体全面的信息，可作为分析变形破坏机理和验证数值模拟成果的必要辅助手段。将传统工程地质定性分析、三维模型解译、数值计算定量评价三种手段相互结合、相互验证能快速提高调查人员对灾害体的认识，为后期灾害治理提供科学有力的依据。

3.1 册亨平庆组滑坡发育特征

研究区位于云贵高原东部脊状斜坡南侧，向广西丘陵倾斜的斜坡地带，地貌单元上属黔桂边境的中低山区，地形上主要为山峦斜坡、河谷、冲沟。斜坡高程900～1 300 m，坡度20°～30°，两侧被山脊所夹持，总体地势北高南低，为典型的堆积体斜坡。平面图如图3所示，工程地质剖面图如图4所示。主要地层岩性特征如下：

图3 剖面位置及研究区域航迹图Fig.3 Track and profile location map of study area

图4 册亨平庆组滑坡工程地质剖面图Fig.4 Engineering geological section of the Pingqing landslide in Ceheng

（1）三叠系中统许满组（T2xm）：岩性主要为灰黑色薄层至中厚层黏土质粉砂岩，夹粉砂质黏土岩，岩层产状335°∠40°。发育两组节理面：J1节理面产状为75°∠59°；J2节理面产状为95°∠84°。

3.2 工程地质条件分析结果

滑坡体东西两侧分别为山脊和冲沟，常年受地表水冲刷使该坡体地形孤立。同时坡体前缘临空，为坡体的滑动以及裂缝的发育提供了较好的地形条件。综合现场情况判断该滑坡破坏模式为降雨条件下由坡体自重应力作用引发的蠕滑-拉裂式滑坡，第四系崩坡积层沿基覆界面滑动。由理论分析得出册亨平庆组滑坡的变形破坏过程分为三个阶段：

（1）后缘拉裂缝形成阶段：在降雨条件下，雨水从覆盖层下渗，在基覆界面受阻，雨水长期作用增大坡体下滑力，当下滑力大于抗滑力时，斜坡蠕滑变形形成拉裂缝。裂缝的形成又为雨水下渗提供了有利条件，通过这种相互强化效应加剧斜坡变形。

（2）前缘鼓胀阶段：坡体蠕滑引起斜坡堆积体应力重新分布，斜坡重心向前缘移动，斜坡前缘位置形成鼓胀区域，区域内产生大量的鼓胀裂隙。

（3）整体破坏阶段：持续降雨时，降雨持续沿后方裂缝下渗至基覆界面，滑面性质劣化，堆积体容重增加，潜在滑面进一步贯通，当堆积体抗滑阻力不足以抗衡下滑推力后，坡体沿前缘剪出口剪出。

3.3 无人机航摄建模及成果解译

（1）数据获取

本次航摄采用大疆MACIVE 2 悬翼无人机。飞行区域面积约0.34 km2，航摄时长约205 min。所用无人机和相机的主要参数见表1。

表1 无人机和相机主要参数Table 1 Main parameters of UAV and camera

规划航线时航向重叠度设置为80%，旁向重叠度设置为70%，研究区域航迹图见图3，整个飞行任务进行了6 个架次，在区域内获得了122 幅正射影像以及454 幅倾斜影像，其影像地面分辨率为0.05～0.10 m（高处低处分辨率不同）。

（2）数据处理

影像数据处理时，采用Pix4Dmapper 软件对正射影像数据进行处理，处理时长约6 h，获取的密点云数据生成测区的正射影像（图5）。采用三维建模软件Smart3D capture（实景建模大师）对影像数据进行三维实景模型建模，建模耗时约73 h，建立的三维模型如图6所示。

图5 数字正射影像Fig.5 Digital orthophoto map

图6 三维模型图（基于Smart3D）Fig.6 3D model map（based on Smart3D）

（3）航摄成果解译

结合裂缝展布规律及地形地貌，判断此滑坡后缘

以现有裂缝为界，前缘以鼓胀区域为界，右侧以地形陡缓相交为界，左侧以斜坡现有冲沟为界，平面形态呈“舌形”，主滑方向190°，纵向长约480 m，横向宽度平均约165 m，滑坡面积7.92×104m2，滑体平均厚8 m，体积约63.4×104m3。灾害体整体分布如图7所示。

图7 灾害体解译结果Fig.7 Results of disaster interpretation

根据影像解译出滑坡后缘产生多条裂缝，其中5 条裂缝规模较大，裂缝宽度5～50 cm 不等，如图7（a）所示。其中裂缝1 长约150 m，延伸方向近79°；裂缝2 长约28 m，延伸方向近91°；裂缝3 长约30 m，延伸方向近136°；裂缝4 长约64 m，延伸方向近46°；裂缝5 长约80 m，延伸方向近42°，局部变形较大的地方出现下挫现象，在主裂缝带四周发育较多的细小次裂缝。滑坡体中部右侧居民区出现1 处小规模垮塌现象，垮塌高度近10 m，垮塌方向近165°，如图7（b）所示。滑坡前部农田发育较多裂缝与垮塌，裂缝宽10～60 cm，局部出现下挫，下挫高度10 cm，裂缝可见深度30～40 cm，裂缝延伸方向40°，如图7（d）所示。滑坡前缘位置出现隆起区域，可见2 个局部垮塌，垮塌方向分别为161°和153°，区域内可见较多近剪出方向的鼓胀裂缝，如图7（e）所示。

根据解译结果结合工程地质分析判断：滑坡目前已进入前缘鼓胀阶段，如果降雨持续劣化滑面及滑体性质，滑坡将进入整体破坏阶段。

（4）无人机航摄数据精细建模

无人机航摄数据处理生成的DEM图像如图8所示。

图8 数字高程模型Fig.8 Digital elevation model

最后通过Rhinocerosceros 软件建立三维地质模型。模型X轴方向长574 m，Y轴方向长760 m，Z轴高差500 余m，建模耗时约1.6 h。模型力学边界采用X、Y方向双向约束，底面（Z方向）垂直方向约束，地表为自由面。计算采用莫尔-库伦屈服条件的弹塑性模型，将计算模型划分为253 495 个单元和46 356 个节点。截取滑坡主剖面进行位移、应力应变分析，模型如图9所示。

图9 FLAC3D 计算模型及剖面选取Fig.9 FLAC3D calculation model and section selection

结合现场情况及钻孔揭露情况，对场地地层进行简化，主要采用黏土岩、粉质黏土两种岩性，力学参数取值参考相邻工程，参数如表2所示。整个计算过程耗时约6.5 h。

表2 物理力学参数Table 2 Values of physical and mechanical parameters

3.4 三维数值计算结果分析

（1）天然工况计算结果及分析

应力场计算结果显示，滑坡前缘最大主应力值约71.4 kPa。滑坡后壁与滑坡体交界处基覆界面附近的局部区域，出现了明显的应力分异。图10 是滑坡在1 000时步时的位移变形云图。与主应力分布差异相关的，斜坡后缘出现了最大主应变，增量集中，沿基覆界线向下出现了量值递减。上覆堆积体可能沿该“潜在变形带”自坡顶位置出现蠕滑。坡体中上部最先产生位移增量，坡顶位置位移达3.67 cm 左右，该位移量值、变形位置与前期无人机航拍监测到的后缘拉张裂缝基本吻合。随着时步推进，位移量值逐渐增大，后缘位移达11.19 cm 左右，最后收敛最大位移达到57.62 cm左右。同时坡体变形范围向前扩展，因此在坡体中部位置的公路路基也产生沿坡体向下的位移，当达到一定量值时导致路面产生横向裂缝，这与勘察到的路面沿路线方向的裂缝一致。

图10 天然工况下位移分布特征Fig.10 Displacement distribution characteristics under the natural conditions

（2）暴雨工况计算结果及分析

暴雨工况下，最大主应力总体上除了量值有所增大外，在分布上未出现明显的调整。最小主应力，在上覆堆积体浅表部附近出现了拉应力条带，量值较小。这是由于表层土体受人类活动扰动，土体孔隙率高、力学性能差，易出现位移变形。如图11所示：滑坡后缘首先出现位移，量值为4.23 cm 左右，随着时步推移，变形逐渐增大且具有向下扩散的趋势，变形量最大值达到84.47 cm。

图11 暴雨工况下位移分布特征Fig.11 Characteristics of displacement distribution under the rainstorm condition

计算收敛后，斜坡拉应力主要集中在堆积体表部（图12）。剪应力具有由后缘逐渐向坡脚贯通的趋势，在坡脚位置主要受剪应力作用。目前通过无人机多期影像比较，发现前缘部分区域有鼓胀隆起现象，并未产生剪切破坏。据此分析目前坡体处于前缘鼓胀阶段，与前文航摄成果解译结果一致。随着时间推移，剪切位移累积，剪应变量逐渐增大，在暴雨作用下，该斜坡将由“时效变形阶段”进入“累进性破坏阶段”，最终整体滑动。

图12 暴雨工况下塑性区分布特征Fig.12 Distribution characteristics of plastic zone under rainstorm condition

4 结论

（1）无人机航摄系统具备机动灵活、响应迅速、成本低廉、适应性强、适用于高危地区探测等优点，使研究人员能在安全区内快速获得地质灾害发生地区高精度的影像数据。通过对航片和Smart3D 建立的三维地质模型快速解译可以准确、全面地分析滑坡的整体发展情况。

（2）将无人机获取的数据以及各控制点的空间坐标信息导入Pix4Dmapper、Global mapper、Polyworks、Rhinoceros 等软件进行一系列处理，建立三维计算模型的方法比传统建模方式具有简单实用、快速的特点。该模型进行数值分析的结果与现场调查结果和解译结果基本吻合，能满足地质灾害应急处置的要求。

（3）本文总结出以轻型无人机航摄为核心的滑坡应急调查、定量评价流程，将传统工程地质分析、遥感解译、数值模拟结果相互验证进行综合评价，提高了应急调查的精细程度，弥补了以往过多依赖技术人员业务水平的局限性，为后续应急处置提供了有力的依据及数据支撑。

（4）册亨平庆组滑坡为降雨条件下坡体自重应力作用引发的蠕滑-拉裂式滑坡。调查时滑坡已启动，地面裂缝有加速变形趋势，坡体稳定性较差，在降雨条件下具有发生大规模滑动的可能性。