基于无人机和Rockfall Analyst的崩塌落石特征分析与运动学模拟——以察雅县崩塌落石为例

刘福臻，李旭德，王军朝，刘建康，陈 龙，张佳佳

(1.西南石油大学 土木工程与测绘学院，四川 成都 610500; 2.中国地质调查局地质灾害防治技术中心，四川 成都 611734)

崩塌落石作为常见的地质灾害，是指斜坡上的部分危岩体受到力的作用后发生滚落的现象，是边坡勘察防护中的重点[1]。引发崩塌落石的条件简单，容易崩落，截至2018年昌都地区就共有2 192处地质灾害及隐患点，崩塌共计647处，占总数的29.52%，其中小型崩塌落石共计403处[2]。目前针对崩塌落石灾害，在形成机制、运动学模拟、易发性等方面的研究成果较多，调查手段也越来越丰富。

董秀军等[3]在2006年就提出将三维激光扫描技术应用于高陡边坡地质调查测量；王凤艳等[4]也早在2008年将近景摄影测量技术用在岩体结构面调查。但是二者对调查环境的要求较高，复杂环境下工作难以开展。随着小型无人机技术的快速发展，凭借其高精度、灵活机动、实用性强，能克服危险环境影响，可提供多角度、高分辨率影像等独特优点，近几年已逐渐应用于各种灾害地质调查的研究[5-13]。

影响滚石运动轨迹和规律的因素多且情况复杂[14]。在国内已有众多学者基于数值模拟软件研究影响滚石运动的参数和运动规律，常用的是二维模拟软件Rockfall和三维模拟软件Rockyfor3D。Rockfall是一个基于统计的软件程序，计算时将块石的崩落过程当作一个质点，在人为规定好的剖面线上让其坠落、碰撞、滚落，应用极为广泛[15-18]。但其崩塌路径是斜坡上的特定剖面，具有一定的人为控制性。Rockyfor3D作为三维模拟软件，在森林地区和铁路沿线的落石灾害评估上都有应用[19-22]。该软件在模拟计算时会考虑单个块石的形状和尺寸大小，通常块石形状精确数据的获取具有一定的难度[23]。

Rockfall Analyst是基于地理信息系统，将分布式建模和物理落石过程相结合的三维模拟软件。模拟时采用集中质量法，使用矢量和栅格数据，能在ArcGIS平台有效呈现落石轨迹、飞行弹跳高度以及速度等的分布，在山区落石灾害预测以及风险评价方面都有应用[24-25]。相比Rockyfor3D和Rockfall，该软件不需要统计落石大小、形状的参数，也无需人为规定落石运动的剖面方向，整体考虑的岩体参数较少，又减少了人为干扰因素，因此更适用于本文研究区环境下崩塌落石的运动学模拟。

察雅县为G349国道起始路段，沿线多为高山峡谷地段，每逢雨季沿线崩塌落石灾害频发，而且随机性大、物源位置分布高、环境复杂，传统的人工手持罗盘等工具的接触式调查无法完成，由于公路上车流量较多，采用近景摄影测量和三维激光扫描技术方法不仅极为不便，而且也无法全面调查评估灾害特征，不能满足调查需求。复杂的地形环境为灾害调查评估和防治带来了较大的困难，严重威胁着过往车辆安全，调查分析研究区崩塌落石灾害特征及威胁范围，可为该处落石灾害防治以及沿线灾害评估提供技术借鉴。因此，本文借助现场调查和无人机三维倾斜摄影测量技术，分析了研究区崩塌落石物源区岩体特征、堆积区特征以及岩体稳定性。获取了研究区高精度DEM，结合基于ArcGIS的三维模拟软件Rockfall Analyst，对物源区块石进行了运动学模拟。

1 研究区概况

昌都市位于青藏高原东部，地处横断山脉以东高山深谷区，海拔3 000～6 378 m，平均海拔3 500 m，地势北高南低。在大地构造上， 区域上处于巨型青藏滇缅构造， 浅层地壳结构特征为板片拼合结构[2]。区域基本地震加速度为0.2g，地震烈度为VII度，50年内超越该烈度值的可能性为10%[26]。察雅县位于昌都市南东侧，县城属于宽谷向峡谷过度的地貌，主河麦曲从县城穿过。县城南侧斜坡基岩出露少，整体为滑坡群，县域内地层以侏罗系汪布组和三叠系夺盖拉组为主(图1)。

图1 崩塌点周边地质图Fig.1 Geological map of the surrounding of the rockfall

本文研究区位于县城北西侧的一处岩质斜坡，坡顶高程3 187 m，坡脚处高程3 094 m。该斜坡基岩以紫红色粉砂质泥岩为主，夹少量浅灰绿色粉砂岩，差异性风化严重。崩塌落石点所在坡体岩体破碎，结构面发育，曾发生过多次小型崩塌落石，坡脚形成堆积区，目前仍存在多处危险源，物源分布位置高，地势陡峭，最低一处物源与坡脚公路的高差为15 m，传统的人工调查方式、近景摄影测量和三维激光扫描技术方法均无法满足调查需求。高位物源块石一旦诱发崩落，对坡脚G349国道、乡道和果园产生着极大的威胁，如图2所示。

图2 崩塌落石位置概况Fig.2 Location of collapse and rockfall

2 研究方法

2.1 无人机三维倾斜摄影测量

相比传统的摄影测量，无人机倾斜摄影可以从多个不同角度拍摄采集数据，作业机动灵活，受地形环境影响小，能够完成传统的人工手段无法进行的高陡边坡危岩体调查工作，构建的三维模型可以真实反映地物特征，实现非接触式测量。无人机倾斜摄影测量的灾害调查流程主要由调查区资料准备、野外拍摄和内业影像数据处理和分析结果现场验证四部分组成。流程图如图3所示。调查区资料包括遥感影像、地形资料、历史灾害调查资料等，野外拍摄包括布像控点和飞行拍摄，内业影像处理主要为影像三维重建。

图3 调查流程Fig.3 Investigation process

2.2 Rockfall Analyst运动学模拟

灾害模拟分析和评估灾害危险性对于灾害防治具有实在意义，常用的分析评估方法为物理力学模型分析法[27]。力学数值模拟模型在牛顿运动定律和碰撞理论的基础上，通过大量试验分析总结出影响落石的特征参数，利用运动学和动力学完成滚石运动轨迹、速度等的模拟[28,29]。本文采用基于ArcGIS的三维数值模拟软件Rockfall Analyst来模拟研究区崩塌落石的动力学。该模型采用了集中质量法，落石动力学不受块石的质量和形状的影响，分析时只考虑震源位置和初始方向的内在随机性。岩石和边坡之间的相互作用由两个恢复系数(RN、RT)和一个摩擦系数特定描述，在模拟过程中，这3个系数是不变的[30]。

该软件可以模拟落石飞行、碰撞反弹和滚动这3个主要的运动过程。落石运动的路径、高度和速度通过抛物线方程来模拟。位移和速度如式(1)(2)所示：

(1)

(2)

其中，g为重力加速度(取9.8m/s2)，X0，Y0，Z0是滚石在三维空间中的初始位置，Vx0，Vy0，Vz0是滚石在x、y和z方向的初始速度。

碰撞和反弹阶段首先是要确定岩石飞行结束时的撞击点，即飞行滚石和栅格表面的交点，然后由质点碰撞反弹模型(图4)表示弹跳速度矢量。弹跳速度矢量根据恢复系数(RN、RT)计算。RN、RT可以在野外现场用高速摄影机搭配解析软件测得[31]，数学表达式如式(3)，弹跳速度矢量由式(4)定义：

(3)

(4)

式中，VDip，VTrend，VN为完整的弹跳速度矢量，VDip为倾角方向上的速度分量，VTrend为法向方向上的速度分量，RN、RT为恢复系数，范围在[0，1]。

滚石沿坡面运动过程中，根据牛顿运动定律对加速和速度变换等模拟计算，其受力分析(图5)和计算公式如下。

m石a=Gsinθ-F阻.

(5)

V末=V初-at.

(6)

3 研究结果

3.1 崩塌落石特征

(1)物源区岩体特征

根据建模生成的报告可知模型平面精度3 cm，高程精度6 cm，满足灾害调查精度要求。从三维模型观察到研究区共有3处物源区，从低到高依次划分为A区、B区和C区，如图6(a)所示。根据模型统计出A区物源方量为710 m3，B区为1 530 m3，C区为450 m3。3处岩体整体受3组主要结构面J1、J2、J3控制，被切割成块状(图6(b))。其中A区块度最大，出露岩面也最新，B区次之，C区最小。物源B区岩体表面X形节理极其发育(图6(c))，可见3处裂缝，由三维模型统计了长度，其中受J3节理控制的两处裂缝分别长1.5 m和0.75 m，J2控制的一条裂缝长2.5 m(图6(d))。

图6 岩体结构面Fig.6 Structural plane of rock mass

除了3组主控结构面，岩体还发育部分次生结构面。由于结构面强度低于岩块强度，不同结构面之间的组合形式产生了多种破坏方式，而且结构面容易积水，高寒地带冰雪冻融交替作用使得结构面间隙扩张，导致进一步加重了崩塌落石的发生。

软、硬岩风化速度不同而形成空腔导致岩体失稳是一种常见的岩体破坏特征。物源区岩体岩性整体为粉砂质泥岩夹粉砂岩，呈“软-硬-软”夹层状。由于粉砂岩的抗风化能力比粉砂质泥岩要强，在受到暴雨冲刷和高原冰雪冻融侵蚀破坏后，粉砂质泥岩风化剥蚀速度加快，从而形成空腔(图7)，粉砂岩失去支撑力，同时又被结构面切割成不稳定块体，在自身重力和外力作用到一定程度时，危岩块体失稳崩落。

(2)岩体结构面特征

岩体结构面是岩体的主要特征，对岩体稳定性起控制作用。岩体结构面在三维点云数据中一般都表现为规则的平面，可以通过提取平面上不在同一条直线上的构成范围较大，且内角均匀的任意n个特征点的三维坐标(n≥3)来计算，倾向和倾角误差可以控制在5°以内[32-33]。

为了验证由三维模型获取的岩体产状的精度，从坡脚选取5处平整的面(图8)并提取每个面的4组三维坐标，基于Matlab语言编写计算程序算出每个结构面的3组产状，求出平均值，然后用手工罗盘实测每个面的产状，将2种方式算得的产状值进行对比，如表1所示。结果发现倾向最大误差11°，最小误差7°，平均误差8.6°，倾角误差最大为6°，最小0°，平均误差2.2°，均符合地质灾害调查要求。

表1 岩体产状验证Table 1 Verification of rock mass occurrence

验证了由三维模型获取的岩体产状精度符合调查要求后，通过提取J1、J2、J3每组结构面的四组三维坐标，同样基于Matlab语言编写的计算程序算出每个结构面的3组产状，求出平均值，并统计结构面间距，如表2所示。

表2 物源区岩体产状Table 2 Occurrence of rock mass in provenance area

(3)落石堆积特征

调查发现堆积区块石都是零散崩落堆积在坡脚。统计了粒径大于1 m的块石共有15处，最大粒径2.3 m。堆积最远的一处块石粒径2.2 m，距离物源A区斜距51 m，距G349国道3.5 m，坡脚处未做任何防治措施，如图9所示。从三维模型统计的结构面间距来看，3组结构面切割形成的块石中，粒径大于1 m的均在物源A区。3处物源区中A区有明显崩落的痕迹，参考上述图6(b)，并且现场可以观察到A区出露的岩面较新，无明显风化现象，再结合块石堆积的位置，综合判断坡脚堆积的大块石均来自于物源A区。确定现有的堆积大块石来自于A区，若B区和C区的块石失稳崩落，能威胁多大范围，需要进行数值模拟分析判断。

图9 崩塌堆积区Fig.9 Partial accumulation area of rockfalls

3.2 岩体稳定性分析

将J1、J2、J3的产状和斜坡坡向、坡度(J4)绘制到极射赤平投影图上，如图10所示。其中J1和J2组合交线的交点位于J4的外侧，因此交线倾角小于坡角，且交线倾向与坡向相同，所以切割体不稳定；J1、J3的组合交线倾向和J2、J3组合交线倾向方向与坡向相反，切割体相对稳定。影响岩体稳定性的除了主控结构面之外，部分次生结构面也可导致切割体失稳。基于HOEK楔形体模型计算出危岩稳定系数，在降雨充水状况下为0.99，干燥状况下为2.2，由此说明在干燥状况下稳定，在降雨充水状况下不稳定，稳定系数参考依据如表3所示。危岩块石失稳后的运动路径复杂且难以确定，需要用数值模拟软件进行分析预测。

图10 赤平投影图Fig.10 Stereographic projection

表3 稳定系数Table 3 Stability coefficient

3.3 落石数值模拟

(1)参数获取

利用Rockfall Analyst对落石模拟分析时，DEM精度、恢复系数(RN、RT)和内摩擦角(Deg)的取值以及落石源点的设置对模拟结果的准确性影响较大。高精度DEM由无人机成果获取，落石源点根据高精度三维模型和正射影像结合现场调查确定。确定研究区恢复系数和内摩擦角时根据历史崩塌落石的堆积位置，采用参数反演的方法确定：先根据参数取值范围确定大概数值，然后进行试算模拟，最后将模拟结果与现场调查的落石堆积位置作对比，最终确定该区域斜坡下垫面的恢复系数和内摩擦角[19]。

通过调整参数，模拟A区落石生成的运动轨迹如图11所示，模拟的运动轨迹最远到达坡脚，没有危及到G349国道。对比发现模拟的轨迹终点与现场调查的块石堆积体分布情况基本吻合(图12)，由此确定生成该轨迹的下垫面参数作为其它物源区块石的模拟参数，如表4所示。在对A区落石模拟时，除了对已崩落的块石进行反演模拟，还模拟了大量潜在崩落的块石，生成的轨迹没有对G349国道造成破坏，因此不再对A区潜在崩落块石的速度和弹跳高度进行模拟分析。

表4 下垫面参数Table 4 Parameters of underlying surface

(2)运动特征分析

根据模拟得出的参数对物源B区和C区块石进行模拟，最终得到落石的运动轨迹如图13(a)、(b)和图14所示，弹跳高度的栅格图如图13(c)、(d)所示。

青岛工学院工程管理类本科的教学模式同样存在重理论轻实践、授课方式单一、课程内容狭窄等问题，本研究主要从以下三个方面探讨工程管理专业教学模式创新[5]：

图13 B、C区落石和弹跳高度轨迹Fig.13 Trajectory of rockfall and bounce height in areas B and C

从图13(a)、(b)可以看出，两处物源区落石都危及到了G349国道和乡道。其中B区西侧落石最远堆积到了国道外侧，东侧落石大部分堆积到乡道，没有危及到果园；C区西侧落石基本都堆积到斜坡中部，对国道没有大的威胁，东侧部分落石不仅危及省道和果园，最远堆积到国道外侧公共卫生场所附近。弹跳高度栅格图如图13(c)、(d)所示，B区落石运动过程中最大弹跳高度为7 m，C区落石最大弹跳高度为6 m。结合落石的运动三维轨迹图(图14)可以看出，B区堆积到国道外侧的落石都是从坡脚处以弧线运动轨迹飞跃过去，C区部分落石直接从斜坡中部弹起飞跃乡道进入果园。结合现场调查，乡道一侧因修路切坡形成一高4 m的陡坎，因此模拟的弹跳高度符合实际。

图14 落石三维轨迹Fig.14 Three dimensional trajectory of rockfall

两处物源区落石的运动速度栅格图如图15、16所示。B区落石速度最大为9 m/s，C区落石速度最大为15 m/s。结合野外调查将B区块石的平均质量设为800 kg，C区块石的平均质量设为200 kg，算得B区落石的最大冲击动能为32.4 KJ，C区落石的最大冲击动能为22.5KJ，如表5所示。

表5 模拟值Table 5 Analog value

3.4 防治建议

结合现场调查和运动学模拟来看，B、C两区落石对公路和果园产生的威胁大，应及时做防治措施。在靠近G349国道这一侧，大部分落石到坡脚处速度为0，如果修筑挡墙，受到落石的冲击力最小，但是无法避免部分飞跃的落石，仍然会危及国道；如果在落石飞跃之前做防治，虽受到落石冲击力最大，但防护效果最好。在乡道这一侧，如果在坡脚处修挡墙能防止B区的落石危及乡道，但是C区的落石将直接飞入果园，因此建议修筑被动网才能做到最有效的防治。

4 结论与讨论

根据现场调查结合无人机倾斜摄影测量技术，实现了对昌都市察雅县一处崩塌落石的精细化调查，并利用Rockfall Analyst软件完成了危岩块石威胁范围的数值模拟。主要结论如下。

(1)分析了物源区结构面的基本信息，得到三组结构面产状：J1：217°∠71°，J2：82°∠38°，J3：315°∠63°，并统计了结构面间距。物源区岩体变形破坏特征主要表现在不同结构面组合形式上的破坏和地层岩性上的差异性分化形成的破坏。

(2)统计了堆积区粒径大于1 m的块石，为数值模拟过程中的参数反演确定了块石来源。通过参数反演确定了研究区下垫面的参数，反演结果和现场进行调查分析的结果基本吻合，保证了后期模拟结果的准确性。

(3)从三维模拟结果来看，物源B区和C区的落石均能危及到G349国道、乡道和果园。B区落石距地面的最大弹跳高度为7 m，C区为6 m；B区落石的最大动能为32.4 KJ，C区为22.5 KJ。通过模拟确定了滚石能够威胁的范围、弹跳高度以及运动速度，有效评估了灾害影响范围。将运动轨迹以三维的形式展示，不仅消除了二维数值模拟软件人为控制因素的影响，而且能直观反映落石运动过程中弹跳和飞行的空间变化，有利于防护位置的选择。通过计算落石威胁到公路时的最大冲击动能，在防护强度的选择上提供了数据依据。

致谢：感谢中国地质科学院探矿工艺研究所提供的实习机会，以及各位老师、师兄给与的指导与帮助！感谢兰恒星老师团队无偿提供的Rockfall Analyst软件，感谢仉义星师兄在软件使用过程中的无私指导！