基于无人机遥感的滑坡灾害风险分析*

2023-05-25王佳颖李玉容杨国柱

工业建筑 2023年2期

杜伟王佳颖李玉容杨国柱

(国网电力空间技术有限公司, 北京 102209)

近年来我国极端天气事件时有发生,受其影响我国多地的滑坡、崩塌和泥石流等地质灾害事件频发,由此造成人员伤亡和基础设施瘫痪等重大的损失。因此,必须加强对地质灾害防治的研究,有效保障人民生命财产和国家基础设施的安全。为了切实提高国家地质灾害防治的能力,我国先后颁布了《地质灾害防治条例》《国务院关于加强地质灾害防治工作的决定》及《全国地质灾害防治“十三五”规划》等法规,计划建立健全地质灾害调查评价、监测预警、综合治理、应急防治四个体系。风险分析包括评估地质灾害发生的可能性以及灾害发生的严重程度两个部分,是地质灾害防治体系中的重要内容。[1]

滑坡是一种典型的地质灾害,具有显著的不确定性和后果损害性,是典型的风险事件。对滑坡进行风险分析和风险管理的研究是近年来地质灾害领域的研究热点和难点。[2]滑坡灾害的风险分析可分为区域和单体尺度。区域性滑坡风险分析是在大范围内对边坡进行危险性的时空预测以及对承灾体进行易损性的调查评价,最终绘制危险性区划图,为地区性的风险管理提供依据;单体滑坡风险分析则侧重于对单一边坡破坏风险的量化分析,其研究结果有助于分析滑坡的破坏机理,是在区域性滑坡风险分析基础上的细化。[3]

关于单体滑坡风险的计算,文献[4-5]提出了不同的计算方法,其中多数方法须要求解边坡发生破坏的概率及灾害导致的损失。边坡的破坏概率与岩土体参数的不确定性密切相关,可由可靠度理论求得[1]。除了岩土体参数,影响单体滑坡风险分析结果的重要因素是边坡的空间几何形态、尺寸以及失效后果的估算。通过传统的测绘方法很难迅速准确地获得复杂边坡表面的空间特征,也很难估算潜在的滑坡体积。

基于无人机遥感(UAVRS)获取边坡的空间形态特征是近年来兴起的一种新方法。UAVRS系统是指以无人驾驶飞行器为载体,通过挂载相机或激光雷达等遥感传感器,在无线通信设备或机载计算机程控系统操控下获取空间遥感信息的一套系统。UAVRS的优点包括设备成本和使用成本相对较低,使用场景灵活,能够快速响应和非接触式等。因此,UAVRS系统已经广泛应用于边坡的稳定性分析、地质勘查和变形监测等方面。文献[6]介绍了基于无人机倾斜摄影的公路边坡三维重建和灾害识别方法;文献[7-8]介绍了使用搭载单镜头的小型无人机获取高陡边坡的岩体结构面参数及对露天矿山边坡建立精细化FLAC3D数值模型;文献[9]介绍了通过固定翼无人机监测露天采坑地表变形的动态变化,为滑坡预警预报和应急救援提供数据支撑。这些研究表明无人机遥感是获取复杂坡体表面空间特征的一种新型方法。但是,如何基于无人机遥感进行滑坡灾害的定量风险分析仍鲜有报道。

因此,依托无人机遥感平台,通过激光雷达采集边坡体的形状及位置信息,利用ArcMap和Rhino软件建立边坡的几何模型并进行网格划分;将使用有限差分软件FLAC3D计算边坡的安全系数,使用一阶可靠度方法计算边坡的破坏概率;通过对FLAC3D进行二次开发,求解三维复杂边坡的失效后果。最后进行滑坡灾害的定量风险分析。

1 工程概况

皖南山区受地形、气候、地质构造和人类工程活动等因素的影响,滑坡灾害多发,人民的生命财产安全和公路、铁路、输电线路等基础设施的平稳运行都受到严重威胁。因此,选择皖南山区池州某县道旁塔基边坡作为研究对象,用以研究无人机遥感在滑坡灾害风险分析中的应用。该边坡由碎石土层及下覆的基岩两部分组成,其土性参数见表1。

表1 边坡岩土体参数Table 1 Soil parameters of slope rock

采用的无人机设备是大疆经纬M600(图1),其参数见表2所示。无人机飞行时,搭载数字绿土LIAIR 220N无人机激光雷达扫描系统,用于采集边坡体的激光点云数据。与通过无人机搭载倾斜摄影相机获取边坡体的光学影像数据相比,对激光点云数据处理后可以滤去坡体表面植被的影响,能更真实地反映复杂坡体表面的空间几何特征。

图1 大疆经纬M600无人机Fig.1 An unmanned aerial vehicle of DJI Matrice 600

表2 无人机参数Table 2 Technical parameters of UAVs

2 边坡精细化几何模型的建立

边坡的几何模型包括地表模型及地层模型,结合二者并进行网格划分可以建立接近实际情况的坡体精细化数值模拟网格模型。

2.1 边坡地表精细化几何模型的建立

许多专家学者已经在基于UAVRS采集地表激光点云数据并处理得到地表数字高程模型(DEM)方面做了很多有益的工作[6,8]。基于DEM文件建立边坡表面(简称“坡表”)精细化三维几何模型,须要用到ArcMap及Rhino软件。其中,ArcMap软件主要用于对*.tif格式的DEM文件进行数据的光滑性处理,以便消除滤波及影像处理算法误差引起的DEM畸变,得到相对平滑的坡表几何形状,并导出带有高程信息的点云数据;Rhino软件主要用于建立边坡的实体模型并划分网格。

2.1.1DEM数据的光滑性处理。

在ArcMap软件中通过“添加数据”功能导入*.tif格式的DEM文件;使用“重采样工具”功能以得到抽稀后的DEM栅格数据。图2为根据无人机采集的激光点云数据得到的边坡DEM信息。以根据DEM生成的等高线为例,其重采样前、后的等高线分别见图3。可见,抽稀后得到的等高线更明晰平滑。

在上述步骤完成之后,须使用“栅格转点”功能将处理好的DEM栅格数据转换成带有高程的点云数据,并保存为*.dwg格式的文件。值得注意的是,此时生成的点云数据并不含有*.dwg文件可以识别的高程信息,还须要在其“属性表”里面自行添加高程字段信息。最后,通过“导出至CAD”功能将文件保存为带有高程信息的*.dwg格式的栅格点云数据。

图2 边坡的数字高程信息Fig.2 Digital elevation information of slopes

a—重采样前的等高线; b—重采样后的等高线。图3 重采样前后的DEM等高线Fig.3 DEM contours before and after resampling

2.1.2边坡实体模型的建立及网格划分

在Rhino软件中导入ArcMap软件生成的*.dwg格式的栅格点云数据文件,使用Resurf插件的“从云点创建单张曲面”功能由点云建立边坡地表的“曲面”模型。使用“挤出曲面”功能由边坡的“曲面”模型建立“实体”模型,并通过“修剪”功能去除模型底部多余的部分,最后使用“将平面洞加盖”功能即可生成如图4a所示的封闭的边坡三维实体模型。

a—坡表几何模型; b—坡体几何模型。图4 坡体的精细化几何模型Fig.4 Refined geometric models of slopes

2.2 地层几何模型的建立

当边坡由成层土组成时,须要在模型中建立土层分界面的位置,即建立地层的几何模型。地层几何模型与坡表几何模型的建立方法类似,即:在ArcMap软件中建立地层界面栅格,再将其导入Rhino软件中并建立地层界面的“曲面”模型。在此基础上,还须通过Griddle插件将地表模型和地层模型组合起来形成一个整体,方法是:1)通过“NonManifoldMerge”命令将边坡实体模型与地层曲面模型相组合。2)使用“testExtendSlits”命令使组合的模型符合Rhino软件的几何规则。3)通过“ExtractSrf”命令去除组合模型的多余部分。

加入地层后的坡体几何模型见图4b。

2.3 网格划分

基于无人机遥感获取的激光点云数据建立坡体精细化几何模型后,为了对边坡进行数值模拟分析,还须要对该几何模型进行网格划分。采用FLAC3D有限差分软件进行边坡的稳定性分析,因此,须要生成FLAC3D所需的网格文件*.f3grid。该项工作可在Rhino软件中实现,步骤[10]如下:

1)使用“mesh”命令建立实体模型的表面网格模型。

2)通过“Gsuf”命令将步骤1)表面网格模型重新划分成符合FLAC3D建模要求的表面网格模型。

3)使用“Gvol”命令将步骤2)表面网格模型转换成符合FLAC3D网格规则的体网格模型,并将其导出为符合FLAC3D格式的网格文件*.f3grid。如图5所示,模型的高约147 m,宽度约为191 m,长度约为354 m,基底厚约为41 m。该网格模型包含83 483个单元和67 285个结点。

图5 Rhino中含地层信息的边坡网格模型Fig.5 A grid model of slopes with stratum information in Rhino

3 滑坡灾害风险的量化分析

风险(R)是指某种特定的危险事件发生的可能性(破坏概率Pf)与其产生的后果(C)的组合,其基本计算式为:

R=PfC

(1)

可见,在单体滑坡的定量风险分析中,两个重要步骤是计算边坡的破坏概率及失效后果,其共同求解基础是边坡稳定性的三维数值模拟。

3.1 边坡的破坏概率

在岩土工程界通常使用安全系数Fs来定量评价边坡的稳定性,并认为Fs小于1时边坡失稳破坏。一般采用可靠度分析方法来求解边坡的破坏概率,并取式(2)所示的极限状态函数:

Z=g(X)=Fs(X)-1

(2)

其中X=[X1,X2, …,Xn]

式中:X为随机变量;Fs为边坡的安全系数,它是变量X的函数。

边坡破坏概率的计算式如下:

Pf=1-Φ(β)

(3)

式中:β为边坡的可靠指标;Φ为标准正态分布的累积函数。

研究使用常规的一阶可靠度方法(FORM)求解可靠指标β,其迭代计算式见式(4):

(4)

式中:μX、σX分别为变量X的均值及均方差。

由式(2)～(4)可见,求解边坡破坏概率的关键环节是求解边坡的安全系数,它可以由FLAC3D自带的强度折减法求得。对于图6所示的边坡模型,在FLAC3D中对模型X方向边界施加沿X方向的固定约束,对模型Y方向边界施加沿Y方向的固定约束,固定模型底部边界的位移。在此基础上,即可采用solve fos命令求解边坡的安全系数。

图6 FLAC3D三维精细化边坡模型Fig.6 A 3D refined slope model constructed by FLAC3D

3.2 边坡的失效后果

边坡失效后果的确定是个非常复杂的问题。为简单起见,参考文献[5]的方法,用滑坡体的体积V来代表边坡灾害导致的损失C,即令C=V。

文献[11]指出:滑坡在宏观上的空间特征就是滑体与滑床相脱离,故可以将二者的分离界面定义为边坡的滑动面。边坡产生破坏时滑坡体和滑床均存在变形,但是在滑裂破坏处的变形相对较大。在基于强度折减法的边坡稳定性数值模拟中,虽然可以由此绘制最大剪切应变增量等值线云来直观显示边坡滑动面的位置,但难以量化表示滑坡体的体积。而从位移的角度来描述这种变形分布时,则可以认为位移等值线最为密集区域即为剪切带,并通过临界位移等值面来表征三维滑坡面。以图7所示的边坡位移等值线剖面图为例,滑动面附近的位移等值线最密集,其位移值较小;滑体部分靠近临空面,其位移则相对较大。因此,通过定义临界位移值来确定滑动面位置的同时,可以依据临界位移值区分边坡的滑体和稳定体,即:认为边坡体中位移大于该临界位移值的单元属于边坡的滑动破坏部分,所有该类单元体积的总和即为边坡的滑坡体的体积。

基于这一思路,在FLAC3D中进行了二次开发,具体步骤如下:

图7 边坡剖面的位移等值线 mFig.7 Contours of displacement for a vertical cross section of slopes

1)载入网格模型后,输入“list zone gp”,得到所有网格单元编号(id)和该单元上的结点编号,将网格信息保存为*.txt文件以便后续调用。

2)在FLAC3D中读取步骤1)得到的网格信息,并将其存入数组中。

3)遍历步骤2)数组中每个单元,并对每个单元的所有结点位移进行判断。该过程涉及的主要命令有gp_xdisp、gp_ydisp、gp_zdisp和find_gp。若单元上的每个结点位移均大于临界位移值,则将该单元的体积保存至一个数组,该过程涉及的主要命令是z_volume。

4)对存储体积的数组求和,即为滑体体积。

该方法的一个重要步骤是确定边坡的临界位移值。参考文献[11]的建议,根据特定剖面边界的剪应变增量极大值点确定滑体剪出口位置,据此选择边坡的临界位移值为0.2 m。以边坡的定值法分析为例,求得的对应的滑体体积为9.683×105m3,约占边坡总体积(7.465×106m3)的13.0%。

基于FLAC3D进行二次开发求解滑坡体积方法的优点是能方便地求解滑动面以上、坡表以下所有单元的体积。当坡表的几何形态由无人机遥感精确获取后,即可准确地确定任意形状滑坡体的体积,在此基础上可以方便地进行滑坡灾害的风险分析。

3.3 滑坡灾害的风险评价

对单体滑坡灾害进行概率分析时,须要明确研究对象稳定性的影响因素及该因素所服从的概率分布。图6所示的边坡由碎石土层及基岩组成,影响边坡稳定性的主要变量是上层土的强度参数值,记上层土的黏聚力及内摩擦角分别为c及φ,则基本变量X=[c,φ]。设基本变量服从正态分布,变异系数为0.1～0.4。在此条件下,边坡的破坏概率、失效后果及风险值与强度参数变异系数的关系分别见图8～10。图中,颜色由蓝至红表示表示数值逐渐增加。

由图8可知:随着黏聚力和内摩擦角变异系数的增大,边坡的破坏概率逐渐增加。其中,黏聚力变异系数对边坡破坏概率的影响较小,而内摩擦角变异系数对边坡破坏概率影响较大。内摩擦角变异系数的增加会导致边坡破坏概率的大幅增加;内摩擦角变异系数较大时,黏聚力变异系数的变化对边坡破坏概率的影响程度减弱。

图8 边坡破坏概率与强度参数变异系数的关系Fig.8 Relations between failure probability of slopes and variation coefficients of strength parameters

表3给出了边坡的性能等级与破坏概率的关系。[12]结合图8及表3可知:当内摩擦角变异系数为0.4时,边坡达到了危险级别,此时须对该边坡开展进一步监测。在黏聚力变异系数为0.1而内摩擦角变异系数从0.1增长到0.4的过程中,边坡的性能从平均以上级别降至危险级别,说明在工程实践中要重视对边坡内摩擦角变异系数的取值。