在我国广西、贵州、云南等西南省份以及广东、江西等省份的部分地区存在较为广泛的灰岩地层。灰岩是一类可溶解的碳酸盐类岩石，在流水的作用下易产生溶解、分崩和破碎并逐渐形成溶洞(穴)、岩溶管道和破碎带等溶蚀现象。水流是灰岩产生溶蚀作用的关键因素，与季节性、周期性的地表径流相比，丰富发达且持久的天然地下水动力场是灰岩溶蚀作用的主要贡献者^[1]。因此，灰岩地层的地质条件与风化程度在时间和空间上存在显著的变异性，地表的地质条件往往与内部地层之间存在较大的差异，且这种差异性随着时间的推移将变得更加复杂多变。理论研究和工程实践表明^[2]，正是由于上述显著的时空变异性导致灰岩边坡与花岗岩、大理岩、白云岩等难溶解的岩体边坡相比，在施工扰动、降雨冲刷等因素的影响下更易失稳下滑，危及施工安全，阻碍施工进度。对此，在施工阶段对其稳定性进行预测评估成为灰岩边坡工程设计与施工中的重要内容。

目前，理论估算法、现场监测法、土工模型试验以及数值模拟分析法是边坡工程稳定性评估预测的常规方法。在工程实践中，通常会采用多种方法进行预测以实现相互验证及提高预测准确性的目的，如现场监测法常与数值模拟分析法相互配合应用。数值模拟分析法涵盖多种二级方法，如：有限元法、有限差分法、离散元法、边界元法和混合元法，其基本思路均为需将岩土体进行网格单元化处理并通过力学参数赋值的方式进行抗力设置^[3]。因此，在进行数值模拟分析时，需准确掌握岩土体的物理力学参数，这是保证模拟分析结果准确性的关键步骤。

为了保证设计、施工的时效性和力学参数获取的便捷性，通常会依据相关规范或者根据工程类比对岩土体的力学参数进行赋值^[4]。不可否认的是，这种做法虽然简单、快速，但有一个难以避免的缺陷，即对技术人员本身的个体主观性和经验性具有极大的依赖性，即便面对相同的岩土体，不同的技术人员也极有可能会选择完全不同的力学参数，导致分析结果离散化甚至完全失效，这显然难以乃至无法为施工过程提供有效的参考信息。此外，天然状态下的真实岩土体，其富水性、破碎性在三维空间里都是不均匀的，地表特征与深部特征之间存在差异, 且规模越大差异越明显，这在灰岩溶蚀地区尤为典型。在这种条件下，简单笼统地根据地表特征难以准确地对岩土体进行级别归类和参数赋值，也很难保证计算结果的有效性及为工程实践提供指导信息。

为了应对这一不足，一种以岩土体实际位移、应力为原始依据并通过一定的数学模型对力学参数进行反算的方法被提出，也就是反分析法。岩土体的实际位移、应力和数学模型是反分析法的两个关键客观要素，在分析过程中可较好地避免个体主观性和经验性的影响，近年来得到了广泛应用，如：凌同华等^[5]以里岩垄隧道YK11+150断面为实际依托，采用监控量测获取了围岩的位移取值, 利用一种改进的粒子群算法和BP神经网络建立了待反演参数与实测位移值之间的非线性关系，并对围岩的力学参数进行了反演，正交试验的结果表明其相对误差较小；徐中华等^[6]依托上海长峰大酒店基坑工程，依据基坑的实测变形，综合采用UCODE反分析软件和ABAQUS有限元分析在二维平面内对基坑土层的竖向抗力单参数m值进行了反演，并以上海地区20个基坑工程为实例进行验证，获得了上海地区典型土层的m值及其取值范围；Caudal等^[7]以加拿大温石棉矿东墙的尾矿边坡为实际依托，综合采用遥感与数值模拟技术进行了滑坡预测，并采用改进权重系数与传递函数的方法对集群应力反分析模型实施了优化，同时对尾矿大型岩块的深部应力和三维位移分别实施了反求与模拟预测；张继勋等^[8]依托某水电站引水隧洞工程，基于遗传算法基本模型编写了改进程序, 对ABAQUS商业软件进行了二次开发并将其应用于地下洞室的围岩力学参数反分析，发现围岩的反演参数大于设计值，利用反演参数进行计算时位移值小于初始参数的计算值；Mahmoudi等^[9]为了提高卡拉季(位于伊朗德黑兰以西36 km处)地铁隧道围岩参数的估算精度及为隧道维护提供指导，利用FLAC^3D软件中的备选单变量反分析程序，基于隧道前期监测获得的位移值对围岩的弹性模量进行了反分析，实践表明基于该参数可将初期喷射混凝土的厚度减少43%；为了根据有限的现场测量数据对初始应力进行反分析，Pu等^[10]利用多输出决策树回归器(DTR)建立了初始地应力场与其影响因子之间的关系模型，使用系数优化纠正子模型策略建立并计算获得了400个DTR全尺寸有限元模型训练样本，结果表明实测值与反分析值的相关系数r达到了0.92，证明了采用DTR及神经网络存储全局初始地应力场的初步可行性；纳启财等^[11]以姜路岭隧道为依托，采用室内直剪压缩试验获取了围岩的天然密度和含水率，并借助数值反分析的手段获得了围岩的弹性模量及其增量取值。毫无疑问，上述研究对于丰富、完善反分析方法与技术具有重要的理论参考价值和实践意义。

然而，诸如上述反分析研究，多数主要集中于算法模型，如函数改进、迭代简化、算法编制等模型优化，忽视了岩土体本身物理性质的时空变异性，且多以单参数或两参数反求为主，忽视了多参数之间的相互影响。这显然不符合天然岩土体的真实条件，也就限制了研究成果的进一步推广与应用。

为了对上述不足进行完善改进，基于理论分析，以广西宜州-河池高速公路No.2标段某灰岩边坡为依托，应用探地雷达实测和BMP90位移反分析模型对灰岩边坡岩体的力学参数进行反求，并据此采用MIDAS-GTS有限元数值模拟软件对边坡的稳定性进行分析，以期为溶蚀区灰岩边坡的力学参数求取与数值模拟计算提供参考。

宜州-河池高速公路是广西壮族自治区高速公路网络规划“四纵六横三支线”的重要组成部分，是连接桂中-桂西北的重要通道。其设计时速为120 km/h，荷载等级为公路-I级，处于黔中高原向广西盆地过渡地带。岩溶特征显著，以丘陵、峰丛洼地为主，植被发育。谷地狭长且边坡较陡，最大标高和切割深度分别达到530 m和200 m。沿线基岩出露，地表多见岩溶管道和落水洞，亦可见石芽和出露钟乳石，局部被地表残坡积黏土覆盖。第四系岩性为残坡积层(Q₄^el+dl)粉质黏土、黏土和碎石；基岩为石炭系上统马平组白云质灰岩、石炭系中统可溶碳酸盐岩、石炭系下统马平组强风化白云质灰岩和泥盆系上统白云质灰岩。

所在地四季分明，自然降水雨量充沛，雨季集中在5—8月份。根据宜州市水文站提供的水文统计资料显示，年平均降水量为1 470 mm。天然水动力场发达，沿线河流属珠江流域西江水系，含龙江、大环江、小环江；地下水类型主要为第四系松散层孔隙水潜水、基岩裂隙水，对混凝土无腐蚀性。

依托边坡位于宜州市德胜镇境内，紧邻测垌隧道，桩号里程范围为K49+105~136，沿路线走向的长度约31 m，属于路堑岩质边坡。边坡整体稍陡峭，坡外角约118°~130°，坡顶至坡底的高程为11.2 m，坡面风化程度不均，坡顶和左侧1/3坡面以强-全风化为主，表层多为全风化土块夹杂灰岩碎石，完整性、稳定性差；坡面中部-右侧整体以中风化为主，但岩体不完整，有明显的节理、裂隙，呈块-块碎状结构；可见该边坡的整体完整性一般且稳定性较差。

2020年5月12—20日当地连续降雨，地表湿润。5月22日上午，边坡中部与右侧部分表层坡体产生滑塌，导致应急停车道旁的边沟被掩埋及部分隔离护栏变形损坏。为了防止落石危险，应急部门于当天下午对边坡实施了简单挂网处理。图 1所示为滑坡现场，可见落石已侵入应急车道，所幸未对过往车辆造成不利影响。因此，尽快对整个坡体的稳定性进行预测评估是接下来需要完成的重要工作，以便为后续处治提供相关参考。

图 1 滑坡现场 Fig. 1 Landslide site

图选项

通过坡体表层地质素描和现场踏勘可知，该边坡的浅层完整性和风化程度显著不均。对于整个边坡浅层坡体而言，全风化、强风化和中风化岩体的比例大致相当，无论将其视为何种级别岩体均不合理。此外，该边坡属于溶蚀区灰岩边坡，当地岩溶作用特征显著，天然的地下水动力场彼此辐射相连，坡体内部易发育溶洞(穴)、破碎带等不良地质体，仅依靠表层地质调查难以发现此类不良地质体, 从而影响边坡的整体性和稳定性预测结果。

为此，摒弃以往对整体坡体进行一次性定级归类的做法，首先进行坡面地质素描并采用探地雷达对边坡内部的富水性、破碎性进行探测，然后再依据《工程岩体分级标准》 (GB/T 50218—2014)的建议对边坡进行分块与单元级别划分。

一般而言，当岩体的级别得到确定时，可以便捷地根据规范进行力学参数的赋值。但这忽视了岩体物理特征的时间-空间变异性，没有考虑力学参数随时间-空间的动态变化性，更多的在于依据工程类比和地区经验实施静态赋值，这显然不符合工程实际。

根据岩体本构关系可知，块体的形变与位移是弹性抗力与效应二者之间的平衡结果，而岩体的抗力又与其自身的抗剪强度指标、泊松比等参数具有直接关系^[12]。此外，由于接触与挤压，块体之间存在位移和应力的传递，从而产生位移和应力的叠加效应。工程实践表明，这种位移和应力是可以通过在岩块上安设具有一定锚固长度的可产生协调变形的柔性构件进行测量估算的，一旦获得了各个级别岩块的位移或者应力取值，就可以通过岩体本构关系对岩块的力学参数进行反求，最终得到不同岩块的力学参数。

为此，基于前述探地雷达探测及岩块单元划分，引入位移反分析技术对完整边坡不同单元岩块的力学参数进行反演。将得到的力学参数应用于有限元数值模拟计算，从而对边坡的稳定性进行预测。最后，采用现场监控量测对模拟分析结果进行对比验证。技术路线如图 2所示。

图 2 技术路线 Fig. 2 Technical route

图选项

探地雷达是近年来在土建领域被广泛采用的一种时效性强、设备简易、操作便捷的电磁勘探技术，其有效探测深度可达20~35 m。探地雷达系统由3部分组成，即：控制电脑、连接电缆和发射天线。探测时，首先在控制主机中设定好发射参数，然后发射天线将根据该参数向目标地层发射具有初始频率、相位、能量和振幅的入射电磁波。

与地震波等机械波不同，电磁波本质上属于交替变换的电场-磁场所形成的统一场在空间中的传播^[13]。天然状态下的真实岩土体是电介质(既有导电性又有绝缘性)，土体内部不同组成物之间的接触面属于电性质突变带。当入射电磁波到达该突变带时将产生反射和折射，且这种反射和折射在土层内部将不断循环直至能量衰减至无法继续传播。理论与试验研究表明^[14]，土层的物理特征会影响其内部组成物的电性质，据此可通过分析反射波的波形特征对内部土层的地质条件实施预测。

在有效探测深度范围内，探地雷达不仅对地层中的大断层、大空洞等大目标体具有较高的分辨率，还对节理、裂隙、小型溶洞(穴)、软弱夹层等具有良好的识别度，对自由水也具有灵敏度高、可判性好的优势，是短距离无损探测的代表^[15]。收-发分离式探地雷达的探测原理如图 3所示。

图 3 探地雷达探测原理示意图 Fig. 3 Schematic diagram of detection principle of GPR

图选项

在前期地质勘察资料的基础上，采用青岛电波所生产的LTD-2500型探地雷达实施现场探测，相关探测参数见表 1。需要说明的是，探地雷达的探测深度随天线中心频率(75 MHz~2.5 GHz)的提高而降低，常规地质预报选用的中心频率通常为100 MHz，此时其探测深度可达20~35 m且同时能满足分辨率的要求，并能完整地对边坡内部的地质情况进行探测，故天线频率选用100 MHz。

由于边坡陡峭且高度较高，实测时无法用人力完成探测，故需要事先在铲车的铲斗上焊接好钢架操作平台，然后选配2名探测人员进入平台并手持雷达天线且要求紧贴坡面。铲车司机需严格控制沿道路纵向的移动速度从而使天线在坡面上缓慢移动以形成水平向测线、测点以保证探测质量。由于不易在竖向控制铲斗的下放速度且容易导致天线脱离坡面，因而仅沿坡面水平向设置测线和测点。

测线的起点、终点和测点个数视现场操作难易而定(如坡面是否凹凸不平)，间距为1~2 m。为了完整地对边坡进行探测，除坡面外还需在边坡周围3~5 m范围内布置测点以扩大探测范围。

探测完毕后，采用自带的IDSP6.0探地雷达分析软件对探测文件进行背景去除、反褶积、信号增益、0点校正等处理。剔除无效文件后，获得线测图像12副(如图 4(a)所示)和点测图像37副(如图 4(b)所示)。综合分析前期地质勘察资料和探测图像的波形特征，对边坡内部的地质条件进行预测，再结合浅层坡体的地质条件综合对边坡进行单元划分。

图 4 探测图像 Fig. 4 Detection image

图选项

统计表明，该边坡共划分为15个单元块，其中包括6个V级单元岩块，位于坡顶和坡体左侧；7个IV级单元岩块，位于边坡中部和右侧坡体上部；2个III级单元岩块，位于右侧坡体的中部和下部。至此，边坡单元分块完毕，后续还需在不同的单元岩块上安设位移计以继续完成单元岩块的力学参数反分析。此外，通过图 4可知，回波反射强烈，线测图像波形杂乱，同相轴断裂不连续，且从图像中段开始存在“雪花状”特征，表明内部岩体多数较破碎，完整性差且潮湿富水^[16]。

天然真实岩体由微小单元岩块和结构面共同构成，单个岩块自身与结构面共同形成相互约束与协调变形，最终在宏观上表现出位移与应力特征，且可以通过安设一定的测试元件进行测量。表面应变计、表面压力盒等面接触式元件虽然也可测量位移和应力，但通常更适合测量接触面之间的接触应力与相互位移，不适用于对岩体内部全尺度范围内的位移进行测量。因此，为了获取不同单元块的全尺度位移值和即时力学参数，首先需通过钻孔在不同单元块上安设可与岩块产生协调变形的位移计(型号：WYJ-100，量测范围0~100 mm)，然后基于位移本构关系和位移反分析模型反求不同单元岩块的力学参数。

事实上，采用人工神经网络开展岩土体位移反分析的研究与实践相继出现且具有一定的可行性，但对有效数据样本的容量、样本训练时间具有较高的要求，时效性与便捷性较差，更适用于非紧急情况^[17]。对此，采用一种对样本数据容量依赖性更低、效率更高的“试凑式”位移反分析模型——BMP90位移反分析模型对不同岩块的力学参数进行反求。

BMP90位移反分析模型是李世辉等人根据边界元法基本理论研制开发的一款适用于各类型工程岩土体的位移反分析算法，早期多应用于隧道工程岩体，后来逐渐在各类型岩土体中得到了应用和完善。其基本原理为：通过不断地调整参数取值使实测值与拟定值之间的差值满足某容许范围，此时拟定值视为可用，如式(1)所示^[18]：

(1)

式中，R为实测值与拟定值之间的差值的容许范围；K₁，K₂分别为拟定值和实测值。

因此，基于上述不同单元岩块的位移实测数据，通过该反分析程序来获取单元岩块的弹性模量E、泊松比μ等力学参数，并将其应用于后续的有限元数值模拟计算中。

根据前述方案，完成了位移计的安设并开始实施位移测量，由于当边坡内部1 m范围内发育有溶洞(穴)时无法安设位移计，因此避开内部溶洞(穴)位置后共安设13根位移计。位移计的安设示意如图 5所示，钻孔直径为76 mm；反分析结果见表 2。

图 5 位移计安设示意图 Fig. 5 Schematic diagram of displacement gauge installation

图选项

由此可见，不同级别单元岩块的力学参数变化规律基本上与其理论本构关系一致，且即便是相同级别的岩块，其力学参数亦不相同甚至差别较大，这较为真实地反映了不同岩块本身所具有的即时力学特征，与根据规范对整个边坡进行统一性赋值相比具有显著的优势。

为了对边坡的下滑稳定性进行预测，采用MIDAS-GTS有限元数值模拟软件建立边坡的网格模型并对位移进行计算，同时采用现场位移监测对计算结果进行验证对比。

以边坡的起止里程桩号为参照建立网格模型，横向宽度为31 m，竖向高度为11.2 m，坡顶往外侧的范围取5 m。根据边界约束条件设置准则，在模型左右两侧设置X方向的约束，在模型底部设置Y方向的约束，默认边坡表面为自由面且考虑基岩裂隙水渗流。与三边形网格相比，四边形网格单元具有更好的力传递效率，因而采用四边形单元生成网格^[19]。岩体内部力和位移的传递本质上依赖于块体接触与相互运动，因而考虑到边坡在三维空间里的地质条件差异性，模型中的网格单元不宜采用等距离设置，而应按照1∶1.25的比例对V，IV级单元岩块的网格进行加密处理，从而提高此类岩体内部节点对位移和力的传递效率以及提高运算精度。

为了节省运算存储空间与加快运算效率，每一次迭代都设置应力释放系数^[20]。假定边坡岩体为弹塑性体，遵守M-C应力准则。为了充分实现岩体内部应力的重新分布与平衡，视应力分配达到平衡为基本原则由软件自行运算，因而不设置固定的迭代次数^[21]。网格模型见图 6，图 6(a)为添加基岩裂隙水渗流压力后显示的网格模型，为了直观显示已隐藏了网格线, 图 6(b)为未隐藏网格线的模型图。

图 6 网格模型 Fig. 6 Grid model

图选项

计算时，各单元岩块分别采用两种力学参数，一种是由前述反分析所得参数，一种是依据《工程岩体分级标准》(GB/T 50218—2014)的建议进行赋值。重度γ取22 kN/m³。具体而言，首先利用BLOCK-SET命令，将每一个单元块设置成隔离体，由此保证独立的隔离体之间具有互不相同的力学参数；然后再对每一个单元块的力学参数进行依次逐个输入；当全部彼此相互独立的单元块的参数都输入完毕时，整体边坡不同级别单元便获得了各自独立的力学参数，最终避免规范建议赋值所存在的盲目性和笼统性。

网格模型的位移云图显示的是整个边坡的位移状态，无法显示不同岩块的位移变化趋势，对于分析单元岩块的位移特征而言没有意义。对此，采用GTS内置的节点提取技术对不同单元岩块的位移值进行提取。为了加快提取速度，每个单元岩块提取一个节点，取不同力学参数时的提取结果见表 3。

由此可见，基于位移反分析所得力学参数的计算值与实测值更为接近，这说明与规范建议赋值相比，通过雷达实测、单元分块及反分析获取的力学参数更加可靠。此外，表中多数数值模拟计算的位移值小于实测值，这是由于网格边界条件限定了岩体的变形范围和方向，而真实岩体是三向变形的。

本研究在理论分析的基础上，以广西壮族自治区典型溶蚀区灰岩边坡为实际依托，采用探地雷达实测和位移反分析获取边坡不同单元岩块的力学参数，并利用有限元数值模拟计算和现场量测对边坡的位移进行了预测对比，结论如下：

(1) 灰岩边坡的溶蚀作用强烈，求取力学参数时需考虑其地质条件的典型时空变异性，宜先对边坡进行分块处理再实施单元块定级归类。

(2) 基于探地雷达实测的单元划分可较好地避免定级时的主观性和经验性影响，在此基础上采用BMP90位移反分析模型可反演得到更加真实的力学参数。

(3) 获取有效的力学参数是保证数值模拟计算结果可靠性的重要前提，今后可尝试综合利用探地雷达和其他人工智能反分析模型对岩体的力学参数进行求取，这是一种新思路、新方法，应继续加以重视。