基于无人机航测的丹霞地貌区危岩结构面识别与三维裂隙网络模型——以重庆四面山景区为例

熊开治，任志远，赵亚龙，杨忠平，张黎健

（1.重庆市地质矿产勘查开发局208水文地质工程地质队，重庆 400700；2.重庆市二零八地质环境研究院有限公司，重庆 400700；3.重庆大学土木工程学院，重庆 400045；4.重庆大学山地城镇建设与新技术教育部重点实验室，重庆 400045；5.库区环境地质灾害防治国家地方联合工程研究中心(重庆)，重庆 400045）

0 引言

丹霞地貌是我国重要的旅游景观资源，同时也是地质灾害的高发区。受地质成因、构造背景、矿物特性、气候条件等多种内、外动力作用的影响，岩坡内部节理裂隙发育，岩体呈现错综切割的“砖块状”特征，多处于欠稳定状态，地区内高位破碎危岩崩塌灾害时有发生，严重威胁游客生命财产安全[1]。已有研究表明，结构面是造成岩质边坡灾害频发的主要原因，因此认清地质背景、查明结构面空间展布特征成为打赢防灾减灾战役的首要前提[2]。然而，受地壳差异性隆升运动的影响，以缓倾厚层碎屑岩为主的丹霞红层山体中多形成垂直向延展的深大优势结构，其控制下的山体不断后退式演化并逐渐发展成孤立的高陡岩坡，沿深切沟谷分布，近距离踏勘工作难以开展；再加上陡直坡面上由差异性风化引起的“凹腔、凹槽”地质结构密集，非接触式勘察精度十分有限。

采用工程措施对地质环境进行改造能够有效的实现实现工程与环境的有机结合的[3]。出于人员安全和勘察结果精准的双重考虑，越来越多的技术手段被运用于实际工程中。比起近距离人工踏勘的高危低效，以及红外激光扫描技术的高成本、低机动性，无人机航勘手段优势突出，并随着航摄技术的不断革新，越来越受地质工作者们青睐。在现有的勘察工作中，以无人机的倾斜航测最为常见[4−5]，但此法在近直立且凹槽、凹腔等微地貌密布的高陡岩坡中明显受限：一方面，局部负角度地质单元图像难以捕捉、信息缺失严重；另一方面，受视距影响，低分辨率的航片难以识别岩坡坡面小尺度结构信息，无法还原岩坡的真实空间结构特征，基于精细化数字模型的稳定性分析工作难以开展。

为此，本文以重庆四面山景区崩塌抢险工程为依托，提出了一种基于大型陡壁式岩坡的水平贴近飞行无人机航测方法，在合理规划航线的基础上，经强大的后处理软件即可快速获得高精度表面数字模型，并通过后期解译、分类统计建立出可供岩坡稳定性分析的空间结构面三维模型。以期为此类地质结构坡体应急抢险工程的勘察工作提供参考。

1 工程背景

1.1 基本概况

2020年12月7日上午11时左右，重庆市四面山景区滴水岩顶部发生局部危岩崩塌。如图1(a)所示，崖壁危岩在长期内外营力改造作用下，稳定性持续降低，最终沿优势弱面组合脱离母岩高位崩落。受山体地质结构影响，失稳仅发生在泥质砂岩岩组中下部，上部保留岩体中节理裂隙发育、卸荷效应显著，孕灾地质条件充分，存在较大安全隐患，岩坡整体上呈受底部“凹腔”结构控制的欠稳定状态，见图1(b)。

图1 四面山崩塌概况示意图Fig.1 Schematic diagram of the collapse of Simian Mountain

为保障旅游高峰期往来游客的生命财产安全，防止崩塌凹腔所致二次灾害的发生，急需围绕危岩带的实际工程特征提出合理的应急方案，其中查明岩坡结构面的空间展布特征是工程地质勘查中尤其关键的一步，能够为短时期内剖析隐患特征、制定防护设施等提供理论参考。

1.2 工程地质条件

研究区位于景区滴水岩段公路沿线陡坡地带，具有典型的“顶平、身陡、麓缓”的丹霞地貌特征，在地质演化历史上正处于丹霞地貌幼年期，表生改造相对强烈，地貌单元长期处于“临空卸荷、应力松弛—岩体崩塌、高位势能释放—陡崖线后退并形成新的陡崖面—继续临空卸荷、应力松弛”的周期性运动状态，崩塌是岩坡系统维持阶段性应力平衡、自组织释能的重要途径。

如图2所示，区内基岩主要由主要为白垩系下窝头山组（K1w）褐红色厚层砂岩、泥质砂岩互层组成，中细粒结构，铁、钙质胶结，山顶沉积有松散堆积物薄层，持水能力有限。砂岩主要由石英、长石组成，相比之下，泥质砂岩受黏土矿物物化性质决定，其力学性能较砂岩差，整体表现出明显的差异风化特征，这也是崖表“凹腔”结构形成的最主要原因。

图2 Ⅰ−Ⅰ'方向地质剖面图Fig.2 Geological profile in Ⅰ−Ⅰ' direction

四面山位于北碚向斜南东翼，区域内无大型褶皱断层发育，新构造运动以及各种风化劣化作用形成了地层中现有构造单元，多呈平行状成组出现。与一般红层岩组相同，受间歇性隆升构造运动影响，地层中多形成正交于岩层层面的近垂直向节理，具有“间距小、延伸远”的特点，与岩层层面有机组合，对岩坡结构稳定性起主要控制作用，决定着坡体的总体演化特征。受当地湿热型气候影响，降水强度大、持续时间久，地下水补给充分，由于顶部松散层持水能力十分有限，地下水沿坡表陡壁流下形成“滴水岩”或沿节理面入渗坡内产生劣化，进一步放大了结构面对岩坡稳定性的控制作用，劣化过程主要分为水动力作用和水化学作用两类，其中后者是由于钙质胶结物或碳酸盐岩块石颗粒的溶解造成的。此外，由于景区开发，天然坡体结构在一定程度上遭受破坏，以坡脚卸荷为主的扰动工程增大了结构面的承载负担，整体稳定性降低。

综上所述，结构面是造成坡体失稳破坏的物质基础，查明结构面的空间分布方式、相互组合关系对分析岩坡变形运动特征、提出合理防治措施来说至关重要。

2 结构面空间信息获取方法

本文针对大型高陡岩坡特征，采用了基于无人机技术的快速便捷、精度高、适用性广的坡表结构面空间信息获取方法。

2.1 航摄设备及基本原理

本文采用4Pro小型无人机进行数据获取，结合Pix4dmapper图像信息处理软件辅助建模，并通过高精度三维数字模型解译出坡面结构面信息（如产状、迹长、密度等）（图3）。

主要由以下几个核心环节组成：①设定航线，拍摄获得研究对象的二维图像；②通过相机定标、特征点提取并基于Sift匹配算法和对极几何原理，得到投影矩阵；③结合二维图像空间信息与研究对象真实三维空间之间的对应关系，将二维图像中满足匹配规律的特征点映射至三维空间坐标系中，形成三维点云；④做平差、去噪处理，提高三维模型的还原度。

2.2 航线布设方案

根据坡体结构特征，同时为保证凹腔处结构单元像素点的分辨率及匹配数满足精度要求，采用全新的基于高精度建模的水平贴近飞行航线布设方案。

像片主距按照相机焦距取值9 mm[6]，1英寸CMOS传感器中像元尺寸取2.3 µm。在综合考虑工作效率和建模精度，将拍摄距离设置在10～20 m，对应图像精度可达2.5～5 mm，能够满足解译工作精度要求。

拍摄过程中，航线应与坡面走向保持平行或近平行关系，以直线型航线为优，对于某些山形空间差异显著的工程坡体，可采用分区或弧线型飞行办法。云台保持水平姿势，相机朝向正对坡面且与航向之间成正交关系。另外，对相邻航片间的图像重叠率作出新的定义，如图3所示，四面山岩坡在垂直方向上具有很强的结构分异性，精度要求较高，采用取值为85%。

图3 无人机水平拍摄基本原理图Fig.3 The basic schematic diagram of horizontal aerial survey

2.3 建立数字模型

按照前节所述方案进行飞行拍摄后，将所得航片进行合理筛选即可采用Pix4dmapper算法程序进行像素点匹配运算并将二维像素点转化为具有三维空间信息特征的点云，进一步平差、去噪，即获得可直接用于节理统计解译的精细化三维数字化点云模型，如图4所示。根据坡表水蚀痕迹可进一步确定出发生本次崩塌的危岩范围（图4虚线范围）。

图4 水平航摄三维点云模型Fig.4 Three-dimensional point cloud model of horizontal aerial survey

3 建立三维裂隙网络模型

3.1 结构面空间特征统计与分析

以裂隙带所在坡表为调查对象，参考勘查工作中常用的节理统计方法，考虑节理岩坡表面的实际坡形形态特征及结构面发育程度，最终在坡表圈定出60 m×30 m的矩形大型节理统计窗口见图1(b)，作为潜在变形区进行研究。在统计平面上，结构面的空间信息由（视向）产状、迹长、间距、延伸距离等要素组成，是结构性和随机性的统一整体，表现为确定性（如：同期节理往往具有相同的产状呈平形状分布）和非确定性（如：各节理迹长、间距等又表现出一定的随机性）的二重特征。因此，为最大限度还原坡体结构的产出特征，须在窗口节理统计变量的基础上，结合平面要素与空间三维要素（包括产状、结构面延伸范围、分布密度等）之间的对应关系，找出三维要素的空间分布规律，从而基于随机算法建立出结构面三维随机裂隙网络模型。具体可分为以下几个环节：

（1）优势分组方法

节理窗口统计法源自统计学思想，旨在通过大量的样本值实现统计变量随机分布规律的无偏估计，从而建立出具高度还原性的随机模型。假设结构面为具有空间属性的薄圆盘[7]，自然条件下，每一个结构面均具有相似的形状，差异特征主要由产状、半径大小、与统计平面交切关系的不同所引起，最后按照一定的分布密度形成裂隙网络集合。为减小统计样本之间的离散性、提高运算精度和速度，现对结构面进行优势分组，将产状相近的节理按照同期结构面近似平行的性质归为一组，组内结构面产状要素统一取平均值。如此一来，将裂隙网络简化为由一个确定性变量和三个随机性变量共同确定的随机模型。

（2）结构面优势产状统计与分析

结构面的空间位置特征通常用产状加以描述，与大地测量坐标系中的空间平面满足一一对应关系。理论上可以通过结构面上非共线三点构造出平面空间方程，从而解译出结构面产状要素：设结构面所在平面由式(2)确定；基于最小二乘法思想，方程系数通过非共线三点P1(x1,y1,z1)，P2(x2,y2,z2)，P3(x3,y3,z3)求得，如式（3）所示；最后基于大地测量坐标系与岩层产状要素之间的对应关系解算出对应倾向、倾角[8]。

基于（1）中所述优势分组原理，对图4所示高精度点云模型进行解译，发现研究区主要由4组优势产状控制下的结构面组成，节理窗口参数统计如表1所示。以岩层层面为例，与现场勘察人员测得结果（倾向300°～320°、倾角8°～9°）比较吻合，一定程度上论证了前述产状解译方法的可行性。

另外，“岩体结构控制论”认为结构面的组合方式直接决定着岩坡变形破坏特征[9−10]，则岩体破坏最有可能发生在上述4组结构面同时参与控制的岩坡结构中。图5所示为满足这一条件的可能存在的几种结构面组合形式。显然，前两种结构面组合本身具有一定的结构稳定性，岩块下落运动趋势受限，为岩坡稳定性的有利组合；相比之下，在c、d结构面组合的控制下，卸荷岩体极易沿“八字形”结构面组合发生张拉破坏，直接控制着坡体崩塌破坏的范围及凹腔的演化方式。四面山崩塌灾害就是在这种不利组合的控制下发生的（如图4所示）。这一结论也再次证实了表1所示结构面优势分组方法及产状解译结果的正确性。

图5 岩坡结构面空间组合形式Fig.5 Spatial combination form of rock slope structural plane

表1 结构面优势分组结果Table 1 Structural plane for the group advantage

需要说明的是，结合表1解译结果，岩层层面为成岩结构面，而J1是影响红层坡体渐进演化特征的主控结构面（如：陡崖临空面是追踪J1而成的），二者皆具有长大构造结构面的特征，在坡体工程尺度范围内近乎完全贯穿，圆盘假设不再适用。而其在研究范围内的空间分布规律相对稳定，因此只需对统计窗口上的平面特征进行延展推测即可获得其空间特征：岩层层面按照统计结果平均值进行取值，产状为313°∠10°、间距1 m；J1产状取平均值210°∠85°，节理间距按照“凹腔”的平均深度进行取值，为1.63 m。相比之下，J2、J3在空间上延伸范围有限，空间特征不确定性显著，必须通过（3）、（4）步骤建立出三维随机结构面系统。

（3）确定半径随机分布规律

同组节理由于与统计平面交切关系（通常用圆心至平面的距离表征）、半径长度的随机性，导致节理迹长在统计平面上也表现出一定的随机性，并服从一定的分布规律。因此，要建立随机裂隙网络模型，就必须先结合统计结果确定出迹长的随机分布规律，再根据其与空间变量之间的对应关系，导出圆心距离、半径长度的随机分布规律。节理圆盘与统计平面存在三种空间交切关系：相交、相切和相离（图6）。本文采用窗口平面上的统计结果代表整个空间上的平均分布情况。

图6 结构面-统计平面交切关系示意图Fig.6 Schematic diagram of intersecting relationship between structural plane and statistical plane

在窗口平面上，节理与平面相交，迹长t与圆心至平面距离d、节理半径r之间满足勾股定理，记迹长的一半为tc，则三者关系由式（4）确定。

可知，已知任意两个变量随机分布规律，另一变量也可求出。假设节理面与统计平面的交切位置是随机的，即d服从（0，r）上的均匀分布[11]，亦即服从（0，1）上的均匀分布。据此，将式（4）变换至式（5），可进一步得到节理半径随机分布规律：首先在（0，1）上生成的随机数S1；通过窗口统计结果得到tc的概率密度函数，从而在对应区间内产生tc的随机数S2；将二者代入式（5），最终得到r对应的概率密度函数及随机数S3，如式（6）所示。

对于J2、J3，迹长t可由节理窗口统计结果获得，解译结果如图7所示。

图7 节理迹长统计结果Fig.7 Statistical results of joint trace length

由统计结果可知，J2、J3节理迹长空间随机特征显著，近似服从对数正态分布，其概率密度函数由对数正态分布的概率密度函数衍生而来。由于拟合值与无人机解译结果吻合度较好，即认为该拟合函数能够代表J2、J3节理在空间中的分布规律。现以J2为例，基于这类分布规律，提出一种迹长随机数生成方法：首先由积分法确定节理迹长的随机数s（如式（7）所示）[12]，式中，r为（0，1）上的平均分布随机数，通常由计算机生成。

将式（7）进行分解，结果如式（8）所示

由图7拟合结果可知，随机变量x满足u0.99=4.375（ up表 示p分位数），即迹长随机数大于4.375的概率仅为0.01，为简化计算，将其视作不可能事件，则式（8）中第二积分项满足：

又由于概率密度函数峰值显著，随机数主要集中在峰值附近（峰值半径0.75 m邻域内累计密度达0.63），因此将峰值赋予s，对第二积分项作进一步简化（由于C值较小，由此造成的误差可忽略）：

将式（10）代入式（8），得到：

上式中，R 可视为（−0.025，1.036）上均匀分布的随机数。如此一来，原问题退化为对数正态分布变量随机数的确定问题。由统计学知识可知，对数正态分布变量随机数 s′由式（12）进行确定[12]：

式中，v 为标准正态分布变量随机数；r1、r2为（0，1）上均匀分布的随机数，二者相互独立。联立式（11）—（13），原问题中的迹长随机数可以表示为：

将式（14）代入式（6），可进一步生成节理半径的随机数：

式中， S1为（0，1）上的均匀分布随机数，与式（13）中的r1、r2之间相互独立。J3节理同理可求出对应迹长随机数的分布特征。

（4）确定节理分布密度

密度为表征节理空间分布特征的又一个要素。半径决定特定截面上节理的产出情况，而密度直接影响节理在空间上的密集程度及相互组合关系。对于密度的模拟，通常采用试算法：首先建立出工程尺度的完整模型模拟坡体；假设坡体中裂隙总数为N，在（3）中所述随机数S1～S3的控制下生成预测三维模型；按照与统计窗口平行的方向任取一平面作为检验窗口，若所得交切节理数 n1与 实际解译总节理数 n之间存在差值，须对N做出 调整，并重复上述步骤，直至满足误差允许值。

3.2 模型建立

为认识丹霞红层岩坡结构，准确把握危岩体控制结构面的空间展布特征，需建立符合丹霞地貌高位破碎危岩的三维随机裂隙网络模型。

结合实际工程，采用三维离散单元法（3DEC），建立出60 m×60 m×30 m的长方体模型区域[13−15]。首先根据前节所得确定性要素（产状）和非确定性要素（半径、节理圆盘圆心距统计平面的距离）构建出由J2、J3组成的随机裂隙模型；基于裂隙总数N反复试算，直至检验窗口统计结果与实际统计结果基本一致，所得结果如图8所示；最后只需将岩层层面、J1空间要素纳入建模之中，形成还原实际工程地质特征的集随机节理、长大结构面于一体的三维裂隙网络模型，结合离散元计算理论即可进行相关计算，为后续工程提供建议。

图8 随机裂隙网络模型Fig.8 Random crack network model

4 结论

以四面山滴水岩崩塌灾害为工程背景，探讨了丹霞地貌山体结构的特殊性及其在工程勘察工作中的不足，并基于此，提出了一种适用于高陡岩坡破碎危岩野外勘察的航测方法，最终结合随机数学思想建立出了三维裂隙网络模型，为认识丹霞红层岩坡结构、建立红层岩坡地质力学模型提供参考。通过以上研究，主要得到以下认识和结论：

（1）丹霞地貌岩质坡体受深切下蚀作用强烈，高陡坡面上多发育“凹腔、凹槽”等地质结构，常规的勘察方法无法同时保障勘察人员的安全性和勘察结果的精确性。

（2）基于无人机航测技术，提出了一种水平式贴近飞行航测方法，能够获取复杂岩坡结构的高精度地质信息。并结合随机数学思想，建立了还原丹霞高陡岩坡结构的三维随机裂隙网络模型。

（3）岩质坡体是随机裂隙、长大构造结构面共同控制下的复杂地质力学系统。根据各自特点分别提出基于图像信息的识别解译办法：前者基于平面统计结果向三维随机空间映射，而后者只需对统计平面特征进行延展推测即可。

（4）建立出的三维随机裂隙网络模型，可进一步结合离散元数值计算进行坡体稳定性分析，为边坡治理、灾害防治提供参考。