基于逆向工程建模方法的危岩体稳定性数值模拟分析*

高相波 李丽慧 廖小辉 王学良 陈子干

(①中国科学院地质与地球物理研究所，中国科学院页岩气与地质工程重点实验室，北京 100029，中国) (②中国科学院地球科学研究院，北京 100029，中国) (③中国科学院大学，北京 100049，中国) (④衢州学院，衢州 324004，中国) (⑤浙江神仙居旅游集团有限公司，台州 317300，中国)

0 引 言

自然界地质体的稳定性分析是复杂的地质力学问题，地质工作者需要在野外地质调查、工程测量等工作的基础上，建立相应的地质力学模型。并采用一定的理论与方法，分析、模拟、再现地质体的灾害发生发展过程，以解决实际工程问题(Huang et al.，1999)。

岩质边坡危岩体失稳的研究一直是学者们关注的重点问题。传统的稳定性评价方法包括定性分析方法(李靖等，1998)，如自然历史分析法、工程类比法和图解法；定量分析方法(郑颖人等，2001；卢坤林等，2013；唐乐人，2014)，如极限平衡法、数值模拟方法等；以及非确定性分析方法(夏元友等，2002；唐慧群，2005；杨天鸿等，2011)，如可靠度分析法、聚类分析法、灰色系统分析法等。目前工程实践中最常用的危岩体稳定性分析方法是极限平衡法，通过计算潜在滑移破坏面上的抗滑力(矩)与滑动力(矩)之比，得到稳定性系数来判断其稳定性(陈祖煜，2003)。并在二维极限平衡法的基础上发展了三维极限平衡法(Huang et al.，2000；Chen et al.，2001；张均锋等，2005；郑宏，2007)。也有学者实现了极限平衡法的有限元计算(Griffiths et al.，1999；邵龙潭等，2001；曾亚武等，2005；于斯滢等，2013)。但是极限平衡法适用于几何形状规则和荷载作用简单的边坡稳定性评价，这与岩质边坡变形破坏受结构面控制的实际情况存在一定的差距，仅凭极限平衡法的计算结果难以揭示危岩体失稳破坏机理。无从得知岩体内部应力分布状态，因此得到的稳定性评价结果也是不准确的。

近年来，数值模拟分析作为一种重要的边坡稳定性评价方法，可以对边坡内部应力场、变形场等的分布特征进行定量化分析，在工程地质领域中得到越来越广泛的关注和应用，解决了许多理论和实际问题。目前较为流行的数值模拟方法主要包括有限单元法(郑颖人等，2002，2004)、离散单元法(Lu et al.，2014；石秋侠，2014；王颖等，2018)、有限差分法(胡卸文等，2019)、非连续变形分析法(裴觉民等，1993；Maclaughlin et al.，2001；吴建宏等，2003；邬爱清等，2008)等。由美国Itasca公司开发的三维有限差分软件FLAC3D以其在工程地质问题分析计算方面快速而稳定的求解优势，成为研究人员理想的数值模拟计算工具。然而FLAC3D软件在复杂地质模型的建立与不规则模型网格的划分等前处理模块中仍然力不从心，存在地质模型过于简化、网格划分过程经常出错以及复杂建模工作量太大、耗时过长等缺陷。进而影响了数值模拟结果的准确性和可靠性。因此有学者通过ANSYS有限元程序完成复杂地质体建模、网格划分，再利用接口程序导入FLAC3D，实现建模的直观化、快速化和自动化(廖秋林等，2005)；基于AutoCAD平台构建三维可视化模型(郑文棠等，2007)；将三维地质建模与数值模拟分析进行耦合，开发地质建模与分析系统VisualGeo，简化数值模拟的前处理工作(李明超等，2007)；基于Surfer平台提取、转换地表及岩层分界面三维地质信息，生成模型数据文件(崔芳鹏等，2008)等。

除此之外，在高陡边坡的工程地质勘察中，研究人员往往受地形因素的影响难以开展实地踏勘工作，制约了边坡岩体结构特征等信息的获取。甚至在潜在危岩体区域的地质灾害调查工作还存在人身安全的问题。因此，有研究人员利用无人机和三维激光扫描仪(董秀军等，2006；Fekete et al.，2013)等新技术代替调查人员进行地质结构的高精度测量和地质结构信息获取工作。然后利用逆向工程建模软件将获取的地形等高线和地层界面、节理裂隙面等点云数据进行CAD曲面模型重建(徐文杰等，2008；邓小龙等，2017)，可以实现三维地质模型的构建。可以看出，已有的数值建模方法基本都是在野外调查结果的基础上，通过建模软件构建地质体模型，再导入数值模拟软件中计算分析岩体稳定性，操作较为繁琐，且过程中可能丢失地质体本身的岩体结构信息。本文基于三维激光扫描技术，快速便捷地获取边坡岩体结构点云数据，在逆向工程建模软件中构建精细化数值模型，尽可能还原岩体本身地质结构特征，最后通过数值模拟软件分析岩体稳定性，并计算其失稳概率。

以浙江神仙居景区地质灾害调查为依托，选取景区标志性景点神龙瀑为研究对象，在野外地质调查的基础上，基于逆向工程建模方法利用FLAC3D有限差分数值计算软件进行危岩体稳定性分析。通过三维激光扫描技术获取神龙瀑边坡的岩体结构几何信息，利用逆向工程软件Geomagic Studio对点云数据进行处理，重构坡体的曲面模型；再利用Hypermesh进行有限单元网格划分，将曲面模型实体化和网格化；最后基于FLAC3D软件进行数值模拟分析，并根据结构面塑性破坏单元数目计算危岩体失稳概率。为当地旅游管理部门的灾害防治工作提供理论指导和依据，避免人员伤亡和财产损失。

图 1 神龙瀑边坡Fig. 1 Shenlong waterfall slope

1 研究区工程地质条件

本文研究区域神龙瀑所在景区属于中、低山区，相对高差200～600im，地形切割强烈，最大切割深度达800im，地势陡峭，局部坡度达到80°以上。常见石峰、石岭、石壁、悬崖等地貌。神龙瀑两侧为高陡边坡，中间为沟谷。东侧山体坡角50°，坡向230°，坡长20im；西侧山体边坡反倾坡角80°，坡向50°，坡长26.5im，坡宽25im。研究区主体岩性为流纹岩、流纹质晶玻屑熔结凝灰岩，局部夹沉凝灰岩。

神龙瀑两侧边坡节理裂隙十分发育，岩体较为破碎，呈碎裂状结构。受多组结构面切割影响，已形成多处危岩体，并频有块石坠落。沟谷中堆积有数百块历史崩塌产生的块石，其中数十块直径在1im以上，最大块石直径达10im。图 2给出了研究区节理走向玫瑰花图，从图中可知北东东向节理最发育，与区域内断层的走向50°基本一致，说明节理发育主要受断层构造作用控制。

图 2 研究区节理走向玫瑰花图Fig. 2 Rose diagram of joint strikes in study area

神龙瀑危岩区有地表水汇集且水量较大，两侧山体表面潮湿，局部有细小水流渗出。地下水类型为基岩裂隙水，补给来源主要为大气降水。坡体表面风化作用强烈，为中风化-强风化，局部岩块用手即可掰碎。坡面上部长有直径5～10icm的树木22棵，树木沿裂缝和节理生长，根劈作用明显。如图 3所示神龙瀑坡面植被发育情况。

图 3 神龙瀑边坡植被发育情况Fig. 3 Vegetation development on Shenlong waterfall slope

2 危岩体形成机制

神龙瀑边坡坡向55°，坡面反倾，倾角达80°，可见多处拉张裂缝，在坡面临空方向形成多处潜在危岩体。且频有掉块和落石坠于坡底。

据野外调查发现，神龙瀑边坡两侧出露宽约1im的断层破碎带，走向50°。受断层构造作用影响，控制该区岩体稳定性的主控结构面走向大致为北东东向(参见图 2研究区节理走向玫瑰花图)。经过野外结构面统计与测量，神龙瀑边坡岩体主要受3组结构面切割，包括一组层理和两组节理，岩体结构呈碎裂状。其中层理面产状225°∠9°，可见范围内迹长10im，间距0.2im；节理1产状340°∠84°，可见范围内迹长10im，间距0.5im，节理面平直，张开度0.003im，局部有泥质填充；节理2产状210°∠78°，局部较为发育，可见范围内迹长5im，间距0.3im，节理面较粗糙，张开度0.005im，局部有泥质和碎石块填充。根据岩体稳定性受岩体结构控制理论，神龙瀑边坡危岩体稳定性受上述3组结构面的控制。

在降雨、地震、风化等自然营力作用下，节理裂隙不断扩展，岩体强度逐渐弱化，极有可能诱发危岩体失稳。降雨、地震和长期风化作用是危岩体发生崩塌的主要诱因，危岩体破坏运动方式主要表现为滑移-坠落。

3 危岩体稳定性分析

有限差分数值计算软件FLAC3D(Fast Lagrangian Analysis of Continua in 3 Dimensions)即连续介质快速拉格朗日分析方法，基于显式差分法对运动方程和动力方程进行求解，且程序内置命令可以设置具有一定接触刚度的接触面以模拟断层和节理等结构面，能够快速而稳定地求解连续介质大变形、大位移问题。本文借助数值计算软件FLAC3D对神龙瀑边坡危岩体进行稳定性分析，为当地景区地质灾害防治提供参考。

3.1 逆向工程建模方法

为了得到更为精确的数值计算结果，软件本身前处理模块的建模工具显然难以胜任地质体模型的准确还原，因此采用逆向工程建模方法对神龙瀑边坡的岩体结构特征进行精细刻画。

首先利用三维激光扫描技术对神龙瀑边坡坡面进行高精度扫描。采用Leica scanstation C10三维激光扫描仪，该仪器单次测量点位精度为6imm，单次测量距离精度为4imm，标靶扫描精度为2imm，模型表面精度为2imm。将获得的点云数据在逆向建模软件Geomagic Studio中进行三维重建，经过网格清理、模型封装等处理手段，将边坡发育的植被部分点云进行降噪、清理，构造高质量NURBS曲面片模型(图 4)。可见重建模型基本保留了边坡岩体结构特征信息。

图 4 Geomagic Studio重建曲面片模型图Fig. 4 Reconstructive surface model by Geomagic Studio

然后将重建的坡面曲面片模型进行四周延伸并三维实体化，构建神龙瀑边坡三维模型(图 5)。此时封装完成的模型由于数据格式的限制尚无法直接导入FLAC3D软件中进行数值计算，且不规则坡面形态的网格划分过程也十分困难。

图 5 神龙瀑边坡三维重建模型图Fig. 5 3D reconstructive model of Shenlong waterfall slope

图 6 Hypermesh有限单元网格模型图Fig. 6 Finite element mesh model by Hypermesh

因此借助有限单元网格划分软件Hypermesh进行处理。根据现场结构面调查统计结果，将封装模型分为神龙瀑边坡、断层破碎带和节理面3类。为尽可能准确刻画模型特征，对重点考察区域的网格进行细化，以提高数值计算精度。本次建模设置岩体有限单元网格划分的尺寸为0.2im，断层和节理单元网格划分的尺寸为0.1im。该模型共包括节点数118i588个，单元数645i805个。

将Hypermesh所生成的模型节点和单元信息文件导入FLAC3D软件中，从而建立数值计算模型(图 7)。在FLAC3D中设置边界条件，模型底部约束竖向位移，四周边界约束水平向位移，模型上表面为自由边界。

表 1 模型分组及其力学参数取值Table1 Numerical model groups and its parameters

模型部位容重/kN·m-3黏聚力/MPa内摩擦角/(°)抗拉强度/MPa体积模量/GPa剪切模量/GPaknks神龙瀑边坡25370420.822.7断层破碎带251.3702.51011节理面0.12506.9e96.9e10

图 7 FLAC3D 数值计算模型图Fig. 7 Numerical model by FLAC3D

3.2 力学参数的选取

本次计算中力学参数主要依据经验值，并根据模型计算结果做一定修正。对结构面的刚度系数采用反演分析方法确定。计算模型的本构关系选择Mohr-Coulomb弹塑性模型。数值计算中使用体积模量K和剪切模量G来代替杨氏模量E和泊松比ν。模型介质分组及其力学参数取值如表 1所示。

3.3 数值模拟分析

图 8给出了神龙瀑边坡在主控结构面切割作用下的危岩体分区。图 9为对应危岩体的塑性破坏区。从塑性破坏区可以明显看出，神龙瀑边坡潜在破坏区域主要分布在上部危岩体。包括一处明显的受拉破坏危岩体和局部剪切破坏的小型掉块和落石。

图 8 神龙瀑边坡危岩体分区图Fig. 8 Unstable rock mass zones of Shenlong waterfall slope

图 9 危岩体塑性破坏区Fig. 9 Plastic damage areas of unstable rock mass

图 10 受拉破坏危岩体主应力分布剖面图Fig. 10 Profile of principal stress in tensile damaged unstable rock massa. 最小主应力；b. 最大主应力

为了进一步对边坡岩体内部应力状态进行分析，在FLAC3D中输出受拉破坏危岩体的主应力分布剖面图(图 10)。由图 10a最小主应力分布剖面图中应力分布状态可知，在自然状态下，神龙瀑边坡岩体内部的最小主应力为拉应力，最大值超过1.74iMPa。且坡体本身反倾，在节理裂隙的切割作用下形成临空危岩体，后缘裂隙贯通，危岩体出现拉应力集中区。神龙瀑坡体除上部岩体分布拉应力集中区外，下部岩体最大主应力主要为压应力(图 10b)，坡脚位置和节理裂隙接触处出现压应力集中区，最大值在25iMPa左右。

图 11 受拉破坏危岩体垂向应力分布剖面图Fig. 11 Profile of vertical stress in tensile damaged unstable rock mass

图 11给出了受拉破坏危岩体垂向应力分布剖面图。由图分析可知，坡体上部危岩体在层理和后缘裂隙切割作用下产生拉应力集中区，拉应力值超过1.6iMPa。在降雨、风化、树木根劈作用等影响下，节理裂隙进一步发育扩展，下部岩体在压应力作用下破碎失稳，进而可能导致上部危岩体发生滑移-拉裂破坏。

图 12给出了自然状态下受拉破坏危岩体位移矢量分布图，神龙瀑边坡整体位移主要分布在受拉破坏危岩体，该危岩体位移矢量向坡体临空方向发展，极有可能产生危岩体失稳，威胁沿路游客的人身安全与景区设施。

图 12 受拉破坏危岩体位移矢量分布图Fig. 12 Displacement vector in tensile damaged unstable rock mass

3.4 基于FISH语言的危岩体失稳概率计算

在野外现场调查结果的基础上，结合上文对神龙瀑边坡危岩体形成机制的工程地质分析结果，危岩体主要受3组结构面切割影响。其自身发育的节理裂隙等结构面控制了危岩体稳定性和失稳破坏方式。因此，考虑计算危岩体结构面的破坏概率作为危岩体失稳概率，评价危岩体稳定性。

在FLAC3D数值分析结果的基础上，利用FISH语言进行二次开发，利用循环语句输出神龙瀑边坡数值模型中各个危岩体接触面单元的应力状态，拉应力值超过抗拉强度的单元定义为拉张破坏单元，剪应力值超过抗剪强度的单元定义为剪切破坏单元。然后统计各个接触面的塑性破坏单元数目，进而计算出塑性破坏单元数目在接触面单元总数中的占比作为相应危岩体的破坏概率。

经计算，出现较多塑性破坏单元的部位在受拉破坏危岩体的后缘裂隙面和层理面，两处接触面破坏概率均达到20.4%。结合上文危岩体稳定性评价结果和FLAC3D的数值计算结果，可以预测神龙瀑边坡的潜在破坏区域极有可能发生在受拉破坏危岩体，应及时采取防治措施，避免景区游客人员伤亡和财产损失。

4 结 论

(1)经过现场地质调查，神龙瀑边坡岩体在节理裂隙的切割作用下呈碎裂状结构，在外力作用下极易脱离母岩，产生掉块。且由于长期风化作用，节理裂隙进一步扩展，岩体强度不断弱化，可能诱发危岩体失稳，其破坏运动方式主要表现为滑移-坠落。

(2)FLAC3D数值分析结果表明，神龙瀑坡体受节理裂隙的切割作用形成多处临空危岩体。其中受拉破坏危岩体后缘裂隙位置产生拉应力集中，其最大拉应力值(1.7iMPa)接近岩石抗拉强度(2.5iMPa)，且该危岩体存在明显的位移，极有可能发生滑移-坠落，威胁景区内游客人身安全，破坏景区公共设施。

(3)基于FLAC3D进行FISH语言二次开发，根据结构面塑性破坏状态，计算出危岩体的失稳概率。结果表明，出现较多塑性破坏单元的部位在受拉破坏危岩体的后缘裂隙面和层理面，其破坏概率均达到20%以上。结合野外地质调查结果和FLAC3D的数值模拟分析，可以预测该危岩体为潜在破坏区域，建议景区管理部门立即对神龙瀑边坡潜在危岩体进行清理或加固。