基于无人机影像的危岩体识别及公路地震风险研究

李培锋,王 晖,吴雨辰,李斯涛,李 春

(1.云南玉临高速公路建设有限责任公司,云南 昆明 677000;2.云南省交通运输厅工程质量监督局,云南 昆明 650206;3.上海交通大学船舶海洋与建筑工程学院,上海 200240)

0 引言

崩塌落石是指具有失稳征兆的危岩体在重力或地震、降雨等外力作用下突然脱离母体向下崩落的现象[1],往往会对坡底的公路、建筑等造成冲击和破坏[2]。近年来,随着我国西南山区等地(如云南省)大力发展公路等基础设施建设,危岩崩塌的影响更加凸显,不仅造成生命财产损失,还会影响工程建设,威胁行车安全。尤其是地震引发的崩塌落石,具有很强的突发性和随机性,更增大了防治的困难。

为有效预防和减轻地震引发的崩塌落石灾害,对危岩体进行调查并分析其运动特征十分必要。但传统的人工勘察方式往往耗费大量人力物力,且难以在复杂地质环境下开展,因此无人机摄影测量技术近几年来广泛应用于地质灾害调查。徐画等[3]通过无人机点云数据实现了危岩体结构面产状的数字化。贾曙光等[4]采用轻小型的单镜头无人机进行倾斜摄影,最终拟合得到岩体出露结构面参数。黄海宁等[5]用无人机进行高陡边坡危岩体调查,得到结构面产状和落石方量信息。目前多数无人机影像解译都是采用肉眼识别的方法统计危岩体参数,当危岩体数量较多时效率不高,本文基于无人机影像,使用图像识别算法,能实现危岩体的自动识别和参数提取,极大提高识别效率。

目前对落石运动特性的分析可分为理论计算[6]、现场实验和数值模拟[7],其中数值模拟方法由于其简便性和灵活性,成为最常用的分析手段。数值模拟方法包括二维数值模拟软件Rockfall[7-10]、三维落石模拟软件Rockfall Analyst[11-12]、非连续变形分析方法DDA[2,13]和离散元方法[14]等。其中Rockfall常用于计算落石运动距离、弹跳高度及动能变化,但由于其在模拟过程中不能施加自重外的荷载,因此无法考虑地震作用[13]。本文结合地震能经验公式,将地震荷载作为初始速度赋予落石,使Rockfall可以用来计算地震作用下的落石运动距离。

本文针对危岩带分布典型的八代村山体,结合无人机影像识别技术和二维数值模拟软件,通过图像识别对危岩体进行定位并快速获取尺寸参数;基于获得的参数,使用Rocfall软件计算地震作用下危岩体的运动距离,从而确定落石影响范围;最后评估坡脚公路受崩塌落石的影响风险,并给出防护建议。

1 研究区概况

研究区域为位于云南省程海断裂带沿线的八代村,如图1所示。程海断裂带北起挖家坪,向南经金官、永胜、程海、宾川,止于弥渡盆地西北,呈倒S状,总体走向近NS,倾向W,全长约200 km,主干断裂分为东、西两支。断裂带沿线地层岩性受影响明显,岩体破碎,裂隙发育,泥石流频发。程海断裂带东缘自全新世以来发生过4次MS≥7.5的古地震。近年来,断裂带周边发生过2001年永胜6.0级、2007年宁洱6.4级、2009年姚安6.0级、2014年鲁甸6.5级、2019年永胜4.9级、2021年漾濞6.4级地震等多次地震。目前程海断裂带还具有较强的活动性,并且在未来一定时间内还可能发生6级左右的中强震[15]。

图1 研究区域Fig.1 Research area

八代村为冲洪积扇地形,海拔高处为物源区,低处为洪积区,地表有明显的冲刷沟,是每年雨季强降水导致的泥石流造成的(图2)。在八代村冲洪积扇地形表面,为防止泥石流对山脚村庄造成破坏,在泥石流路径上修建有多道拦石坝。同时在作为物源区的山体表面,散落着大量泥石流遗留下来的落石(图3)。地表落石形成危岩体,在扰动的情况下虽然处于稳定状态,但在地震动作用下极易发生滚动,导致崩塌滚石灾害,对坡脚的公路及村庄造成威胁。

图2 八代村地形Fig.2 Geography of Badai Village

图3 危岩体示意图Fig.3 Schematic diagram of dangerous rock mass

云南省大永高速公路k97+020(八代村段)通过巨型冲积扇,位于物源区坡脚,如果坡体表面的危岩体发生崩塌,极有可能对高速公路产生影响。尤其是在地震发生后,落石对公路的破坏会严重阻碍人员疏散及应急救援工作的开展,因此需要在考虑地震作用的情况下,对物源区表面危岩体可能对公路造成的风险进行分析。

2 基于无人机影像的危岩体识别

2.1 数据采集

危岩体的影像数据使用瑞士生产的eBee无人机进行采集。利用 Pix4Dmapper 软件对所拍摄的图像进行全自动处理,生成高质量、空间参照式的二维和三维图像,以及数字表面模型(Digital Surface Model,DSM)。生成的三维点云可达到分米级建模的要求,拼接后得到的全域影像如图4所示。

图4 八代村无人机影像及研究区域Fig.4 UAV image of Badai Village and the research area

2.2 图像处理

(1) 研究区域选取

危岩体物源区集中于八代村西北山坡,因此在无人机影像中选取左上的坡体部分作为研究区域,(图4)。

(2) 颜色空间转换

无人机航拍所获得的照片是基于RGB颜色空间的。同RGB颜色空间的设备依赖性相比,Lab是一种设备无关的颜色空间,也更加容易调整。在Lab颜色空间中,L代表亮度,a和b是两个颜色通道,分别代表从绿色到红色与从蓝色到黄色的分量[16]。在本研究中,统计Lab值与RGB值的分布能够发现Lab空间呈现的各元素特征更为明显,因此首先将无人机影像从RGB空间转换到Lab空间,以便于后续图像识别处理。

(3) 滤波去噪

在进一步处理前需进行滤波去噪,去除干扰信号并提升图像质量[17]。均值滤波是图像滤波最常用的手段,其优点在于计算速度较快,并且能够去除尖锐的噪声,使得图像中的噪声像素的灰度值更加平滑。本研究采用均值滤波的方法,保留图像的完整性。

(4) 裸土与植被元素去除

在无人机影像中,将像元元素划分为危岩体、植被、裸土三类,目标危岩体仅占整幅图像中极小一部分,而绝大多数都是裸土与植被等无关元素。因此本文提取目标元素时采用基于阈值分割的无关元素去除法。

常用的阈值分割方法包括p-分位数法、大津法(OTSU法)、直方图双峰法、手动法和迭代法等[18]。在本研究区域中,元素种类较少,手动法可根据实际图像Lab值频数分布人为对阈值进行反复调整,以获取最佳分割效果,且操作简单、易于实现,故选择手动法进行阈值分割。图5为无人机影像的Lab值频数分布。

图5 无人机影像Lab值频数分布Fig.5 Frequency distribution of Lab value in the UAV image

结合手动查询对比,各元素的a值差异较为明显[图5(b)],裸土元素a值多在11以上,植被元素a值多在0以下,而目标危岩体a值基本分布在0～11区间,因此初步筛选出a值位于0～11区间的像元。余下部分主要包含植被与危岩体,其中植被L值大多高于30,且b值大多高于20,因此去除上述部分,即可大致得到只有危岩体的图像,如图6(b)所示。但在该步骤中,由于危岩体的阴影部分与植被Lab值相近,存在去噪过度的状况,从而导致了识别误差。

图6 图像处理Fig.6 Image processing

(5) 二值化与区域去除

将图像二值化处理,以便于明确危岩体位置,如图6(c)所示。二值图像中存在一定的大小区域干扰,通过人工判别,过小的区域往往是与危岩颜色相近的土地上稀疏的杂草,过大的区域基本上是与危岩颜色相近的拦石坝。去除其中面积过大或过小的区域,图像分辨率为40 px,小区域面积取10 px,大区域面积取40 px,所得部分即为最终目标,如图6(d)所示。

2.3 识别结果

将图6(d)与初始图像重叠,得到能够直观呈现危岩体的识别结果(图7)。如图7所示,共识别出危岩体32处,而通过人工判断能够确定剩余未识别的危岩体共10处,说明所提方法的识别精度达到76.2%,具有良好的识别效果。

图7 识别结果Fig.7 Recognition result

对于未识别情况,主要是由于所取阈值未能完全区分危岩本体与裸土、植被部分所致。在阈值选取中,裸土与植被的a、b值不可避免地和危岩部分存在交集,因此只能尽量择优取值。

由于本研究采用的是无人机正射摄影,无法获取危岩体高度。在对每处危岩区域的识别中,长边长度识别是相对准确的,而短边长度的识别则受到危岩体阴影的影响,故以长边长度作为唯一可靠的变量。以长边作为直径,将危岩形体拟合为球体,岩体密度取现场调查得到的2 600 kg/m3,结合图像识别得到的直径,计算其面积和体积,用于后续的数值模拟计算。统计危岩体参数如表1所列。

表1 图像识别得到的参数Table 1 Parameters obtained by image recognition

3 基于Rocfall的地震作用下危岩体运动距离分析

Rocfall是一个结合统计原理,对危岩体沿边坡剖面崩塌坠落过程进行模拟的二维数值计算软件[11]，利用其能够分析得到崩塌落石的速度变化、运动距离、弹跳高度等参数。目前许多研究使用Rocfall分析危岩体崩塌破坏的运动特性[10,19],取得了良好的模拟结果,能够为危岩体崩塌的工程治理提供合理参考。

3.1 Rocfall计算原理

Rocfall软件中,崩塌落石沿给定的坡面运动,运动过程遵循牛顿运动定律,并做出如下假定[13]:(1)边坡坡面由若干折线段组成,忽略坡面微小突起或凹陷对落石运动的影响;(2)落石与坡面为刚体碰撞,在滚落过程中不考虑二者的变形;(3)不考虑空气阻力对落石的影响。

软件通过定义摩擦系数控制落石运动速度的变化,通过引入法向恢复系数(Rn)和切向恢复系数(Rt)衡量落石和坡面碰撞前后的能量损失[13],恢复系数定义为碰撞前后速度的比值。计算时摩擦系数和恢复系数的具体取值是根据坡面物质组成和粗糙程度确定的经验值。除摩擦系数和恢复系数外,软件的计算还需设定危岩体的形状和质量、初始状态、坡面形态等参数,并结合统计原理考虑这些参数的不确定性[20]。

由于Rocfall软件无法在模拟过程中对落石施加除重力作用之外的荷载[13],因此通常无法用来计算地震荷载作用下危岩体的崩塌过程。针对这一点,本研究结合前人提出的经验公式,计算得到地震作用传递给危岩体的总能量,并将这一能量转化为初速度施加到危岩体块石上,从而实现Rocfall对地震作用下崩塌落石运动的模拟。

3.2 地震荷载的施加

在程海断裂带附近,永胜县历史地震尤其频繁,20世纪以来发生过2001年MS6.0和2019年MS4.9等中强震。永胜县所在的程海—宾川断裂带自晚更新世以来一直具有较强的活动性,历史上发生过多次强震,如1515年断裂带北段活动引发的MS7.5永胜红石崖地震,是滇西北地区有记录以来震级最大的地震。此外,1515—2001年间程海—宾川断裂带共发生了9次MS5.0以上的强震[15],小震则更为频发,因此可以推断该地区未来仍有极高的地震风险。选取20世纪以来永胜县内震级最强的MS6.0永胜地震(震源深度10 km)作为本次地震荷载,考虑最极端的情况,假设震中在研究区域1 km之内,即所有危岩体与震源的空间距离都可近似看作10 km。

Gutenberg于1995年根据历史地震数据提出了计算地震波总能量的经验公式[21]:

lgE0=1.5M+11.8

(1)

式中:E0为地震波的总能量;M为地震震级。

假设地震的主要能量是由垂直传播的SH波传递,并且波的能量按球形辐射,则可以用下式表示地表底层单位面积获得的地震能[22]:

(2)

式中:EIP为地表底层获得的地震能量;A为获得能量的面积;R为地表底层距震源的空间距离。

为了计算的便捷,式(2)忽略了地震波从震源传到坡底过程中的衰减[21],因此其计算出的坡底获得的地震能量大于实际值。将该能量取值代入后续计算,所得的危岩体运动距离同样将大于实际距离,所以可以认为本研究的分析结果是偏于安全的。最后假定所有传递到坡体的能量全部耗散,不存在从坡体顶层反射回底层的波,则可以用下式计算真正从地表底层传递到坡体顶层的能量[22]:

(3)

式中:EEQ为坡体顶层获得的经过耗散的地震能;α为阻抗比,是坡体顶层和底层土体密度和波速的乘积之比,通常按经验值取0.3。

由此,就可以在已知地震震级和震源到危岩体空间距离的情况下,根据各危岩体的面积计算其获得的地震能。假定地震对危岩体的能量全部转化为初始动能,由动能公式即可换算得到危岩体初始的水平速度，这样就可以通过在Rocfall中赋予初始状态参数,模拟地震作用下危岩体的运动路径。

3.3 基于Rocfall软件的二维数值模拟

使用Rocfall软件进行数值模拟,首先需定义危岩体崩塌所在的边坡剖面。由于危岩体数目较多,且分布在物源区的不同位置,因此选取物源区内的4处主要沟谷作为典型剖面(图8)。将识别得到的32处危岩体分别划分至邻近剖面,根据危岩体的高程在剖面对应高度上设置运动起点,然后将图像识别得到的危岩体质量以及地震传递能量计算得到的危岩体初始水平速度作为初始运动参数，赋予对应的危岩体。

图8 剖面选取Fig.8 The selection of profiles

在设定坡面的摩擦系数和恢复系数时,据现场勘察结果确定坡面主要为少量植被覆盖的土壤,因此参考前人文献，设置边坡的切向恢复系数、法向恢复系数、摩擦角分别为0.32、0.83和30°[19]。

考虑参数的不确定性,利用Rocfall软件对每个危岩体进行50次运动轨迹模拟[10]，选取50次模拟结果中危岩体运动距离的众数作为最终的运动距离。图9为剖面1中高程1 539.09 m的危岩体利用Rocfall分析得到的运动轨迹。以图9为例,可以看出危岩体由于地震作用获得较大的水平初速度,首先做抛物线运动,与坡面发生碰撞后危岩体的水平速度迅速减小,然后与坡面发生多次碰撞弹跳,最后整体沿坡面滚动直至停止。统计32处危岩体的运动距离,并绘于图10。结果显示,危岩体运动距离在8.3～75.2 m间,大部分位于40 m附近。

图9 运动过程模拟Fig.9 Movement simulation

图10 运动距离统计折线图Fig.10 Line chart of movement distance

为了分析坡度和岩性因素对危岩体运动距离的影响,使用控制变量法,将同一危岩体放置在不同坡度(10°、15°、20°和25°)和岩性的坡面上,使用Rocfall对危岩体运动进行模拟。3种不同岩性类型的参数取值列于表2,模拟结果列于表3。从表3中可以发现,随着坡度从10°增加至25°,危岩体运动距离单调增加,特别是坡度超过20°时,运动距离的增加非常显著。当坡度相同时,从覆盖植被的地表到清洁基岩,危岩体运动距离呈明显的增加趋势。这是因为随着坡面粗糙程度的减小,危岩体滚落过程中损耗的能量更少,恢复系数越来越大,最后计算得到的运动距离自然更大。总结上述规律,可以认为坡度和岩性都对运动距离有较大影响,尤其是当坡度超过20°以及岩性非常光滑时,运动距离可以达到相当大的值，不过在实际工程中,这种理想状况通常不会出现。

表2 边坡基本参数设置Table 2 Basic parameters of the slope

表3 不同参数下危岩体运动距离计算结果Table 3 Calculation results of movement distance of dangerous rock mass under different parameters

如图11 所示,在地理信息平台ArcGIS中,以危岩体所在位置为中心,以运动距离为半径,绘制圆形缓冲区,即可作为危岩体在地震作用下可能的影响范围。可以看出,代表影响范围的红色区域和大永高速公路没有重叠,但是在最南端存在较为接近的情况。测量得到危岩体在地震作用下的影响范围和公路的最近距离为57 m,存在一定的安全风险。由于危岩体在坡脚基本上处于滚动状态,不会产生飞跃的落石,因此公路与山脚的相邻路段应该修筑挡墙作为防护手段。

图11 危岩体影响范围Fig.11 Influence range of the dangerous rock mass

4 结论

针对具有典型地貌特征和大量危岩体分布的八代村,结合现场调查和无人机影像识别技术,实现了对危岩体的智能化识别;在假设物源区附近有强震发生的情况下,使用Rockfall二维数值模拟软件计算危岩体的运动距离,分析危岩体对坡脚高速公路的影响,所得结论如下:

(1) 使用图像识别方法,基于Lab色彩空间对无人机影像进行分析,成功剔除了草地和土壤后共识别出32处危岩体分布,与肉眼识别和现场调查进行对比,精度达到76.2%。通过自动化的图像识别方法,实现了危岩体位置及尺寸参数的快速获取。

(2) 使用地震能量计算的经验公式,通过赋予危岩体水平初速度的方式施加地震荷载,结合图像识别得到的危岩体尺寸，在Rocfall软件中计算得到地震作用下危岩体的运动距离,结果显示在MS6.5地震作用下危岩体的运动距离为8.3～75.2 m。

(3) 结合危岩体的位置和运动距离,确定危岩体的运动范围,发现其影响区域和坡脚的大永高速公路没有重叠,但与公路直线距离最短仅为57 m,因此判断出公路的危险段,并给出了相应的防护建议。