凌 晴，张 勤，张 静，瞿 伟，孔令杰，朱 丽，张金辉

1. 兰州理工大学土木工程学院，甘肃 兰州 730050; 2. 长安大学地质工程与测绘学院，陕西 西安 710054; 3. 信息工程大学地理空间信息学院，河南 郑州 450001

我国具有世界上黄土分布面积最大、黄土发育最完整且唯一正在堆积的黄土高原[1]。黄土本身具有极强水敏性和独特强度衰减性等特点，导致黄土滑坡灾害愈演愈烈，严重威胁到人们生命财产安全和生存环境[2]。为减少及有效避免黄土滑坡灾害带来的巨大危害，亟须开展滑坡有效预警研究，而滑坡稳定性评价的准确与否是滑坡预警的重要前提。受多场及多因素互馈劣化作用影响，黄土滑坡变形具有复杂性、多样性、随机性的特点[3-6]，因此准确的黄土滑坡稳定性分析与评价是滑坡灾害防灾预警研究亟待解决的重要课题之一。

现有滑坡稳定性评价方法依据评价手段的不同主要可分为大地测量领域中的位移时间序列预测法和工程地质领域中的基于可靠度理论、有限元或极限平衡理论等数值模拟法。前者以滑坡外部多源高精度形变监测信息为基础，以变形及其变化量(包括累积位移、变形速率和切线角等)[7]为预警判据来判别滑坡当前所处的稳定状态，以及未来可能的变形发展趋势，具有模型易实现、计算精度优等突出特点，已在滑坡稳定性评价中取得了广泛应用[8-18]。文献[8]首次提出滑坡蠕变三阶段理论，并将其成功应用于日本井田铁道滑坡、高汤山滑坡的预测。文献[9]运用Voight加速蠕变理论得到了滑坡特征速率曲线，并获取了滑坡预报、预警和失稳阶段的速度阈值。文献[10]基于多种机器学习耦合模型进行了滑坡位移预测研究。文献[11]基于集成模态分解-随机森林模型开展滑坡时序预测研究。文献[12]提出了一种基于神经网络的全数据驱动滑坡时序预测模型。后者在工程地质中，滑坡稳定性评价以极限平衡理论、有限元或可靠度理论为基础，以滑坡失稳概率或安全系数为失稳判据来评价滑坡稳定性[19-29]。该方法是一种没有顾及时间变化的物理力学评价模型，能够解释滑坡滑移的灾变机理，适用于滑坡静态稳定性评价。文献[19]提出了极限平衡法。文献[21]采用有限元强度折减法研究滑坡稳定性。文献[22]将滑坡视为系统，运用可靠度理论分析了边坡稳定性。文献[27]通过对复合土滑体的非稳态非饱和渗流数值模拟，探讨了复合土滑体的再激活机理。文献[28]基于有限元分析法量化了现有及潜在滑坡的危害程度。

尽管上述各类滑坡预测预警方法被不断组合与完善，使得滑坡稳定性评价结果也有显著改善，但仍存在以下不足：大地测量位移时序法实质上还是一种数学表观模型，仅能反映滑坡位移随时间的变形规律，且位移时序曲线受外部环境因素的影响常呈现多期加速阶梯振荡趋势，并没有顾及滑坡内部地质条件的影响，导致难以准确判断滑坡稳定性；工程地质数值模拟手段实质也是一种物理力学静态评价模型，仅能基于简化的边界条件及估计的物理力学参数构建模型，且并没有与滑坡外部高精度监测信息相结合，导致无法评价滑坡随时间变化的动态稳定性，进而难以实现滑坡动态预测及预警。因此，上述两类稳定性分析方法均具有一定的局限性，而将大地测量中滑坡位移演化过程与工程地质动力灾变机理进行耦合，建立一种优势互补的滑坡稳定性综合评价模型，可进一步加深对滑坡失稳动因和灾变机理的深入认知，为滑坡防灾减灾预警研究提供十分重要的理论与实际参考。

为此，本文结合大地测量监测信息、地下水、工程地质钻孔信息及灌溉资料等，提出一种融合工程地质资料和GNSS高精度监测数据约束的黄土滑坡稳定性评价方法。基于高分辨率影像、高精度DEM、地下水、地层地貌等多源异构监测数据精细建立起滑坡三维地质模型，并将滑坡外部GNSS高精度监测数据作为位移约束进一步施加到模型上。选取我国典型黄土滑坡域甘肃黑方台党川实际发生的两起黄土滑坡失稳事件进行实例分析，进一步验证了本文提出的滑坡稳定性评价方法优于现有工程地质数值模拟法。

1 模型与方法

本文首先基于滑坡GNSS高精度监测网络、高精度DEM、地下水位、高分辨率遥感影像及工程地质钻孔信息等多源异构数据，借助Arcgis-Rhinoceros-Griddle构建出复杂滑坡精细三维数值模型；然后，将滑坡GNSS高精度监测网获取的滑坡外部高精度监测信息作为模型外部边界位移约束，进一步构建起高精度监测数据约束的滑坡精细三维地质模型；最后，融合区域灌溉、孔隙水压力等工程地质资料开展滑坡稳定性综合评价研究。具体算法流程如图1所示。

图1 融合工程地质资料与GNSS高精度监测信息的甘肃黑方台党川黄土滑坡稳定性研究Fig.1 Algorithm flow of the stability evaluation of Dangchuan loess landslide in Heifangtai area of Gansu province based on integration of engineering geological data and GNSS high-precision monitoring information

1.1 滑坡数据

(1) 高精度多源监测数据。本文试验分析对象为位于甘肃省永靖县盐锅峡镇的黑方台党川黄土滑坡群[30-32]，台塬上部和下部相对高差达70 m，滑坡主体距离下方村庄等生活聚集区约400 m，正下方为水渠和农田(图2(a))。区内人类工程活动扰动密集，黄土滑坡发育，严重威胁该区人民生命财产安全。为此，党川段从2018年10月起共建立12个GNSS连续监测站(图2(b))，并在裂缝附近布设了一系列裂缝计[33-35]。此外，为监测区内地下水位，黑方台台塬共布设了大约40个钻孔，其中32个钻孔孔底位于黄土层，8个钻孔孔底位于基岩层。因此，可获得党川段滑坡群的GNSS高精度监测数据、工程地质钻孔信息。综上，本文选取党川段中部变形较大的区域，即HF05/HF06/HF07/HF09 GNSS监测点所在范围，进行附加外部高精度监测信息下融合多源异构数据的滑坡稳定性评价。

图2 研究区滑坡位置及GNSS监测点Fig.2 The location of study area and GNSS monitoring points

(2) 高分辨率DEM。为建立复杂滑坡精准三维地质模型，需要高分辨率的DEM。因此，采用无人机在研究区航行，并以WGS-84坐标系为基准得到了研究区0.5 m DEM，再结合党川滑坡HF05/HF06/HF07和HF09 GNSS监测点的位置确定三维模型具体范围，然后将黑方台0.5 m DEM在ARCGIS中进行裁剪，获取研究区0.5 m DEM(图4)。

图3 HF06/HF07/HF09/HF05 GNSS监测点E、N、U方向累积位移Fig.3 The cumulative displacements of HF06/HF07/HF09/HF05 GNSS monitoring points in the E, N and U direction

图4 黑方台台塬0.5 m DEMFig.4 0.5 m DEM for Heifangtai

1.2 基于Arcgis-Rhinoceros-Griddle复杂滑坡精准三维模型构建

滑坡稳定性确定性评价方法的关键是求解出合理的安全系数，在解决该关键问题上多采用二维稳定性评价方法。然而，当前滑坡所处地质环境较复杂时，基于简化的滑坡二维模型会导致计算结果与实际存在较大差异，进而会降低计算结果的可靠性。研究表明，基于二维简化模型获取的安全系数小于三维模型得到的结果，会低估滑坡的稳定性[26]。为此，本文针对当前研究构建模型与滑坡实际不符造成评价结果不可靠的问题，将融合多源异构数据，构建出能真实反映滑坡实际三维运动变形特性且具有现势性的滑坡精细三维地质模型。具体过程如下。

(1) 依据滑坡地理位置，运用Arcgis裁剪高精度DEM，并将其导入到Rhinoceros中得到滑坡等高线，同时生成GNSS监测点。而后运用绕轴旋转、平移、MeshPatch等命令得到如图5(a)所示的网格。

(2) 基于GNSS监测点，运用布帘曲面、曲线、从物件建立曲线、复制边框等生成垂直壁。随后在z=0处建立一条水平线，基于水平线修剪垂直壁，并将其与布帘曲面进行组合、加盖，形成一个封闭的多面体(图5(b))。

(3) 依据研究区工程地质钻孔信息、滑坡失稳具体位置设置切割线，并运用切割线将实体进行切割，得到具有地质信息的滑坡三维实体(图5(c))。基于Griddle的Gsurf、Gvol命令生成四面体网格(图5(d))，从而获取FLAC3D可识别的网格模型(.f3grid)。

图5 复杂滑坡精准三维模型构建过程Fig.5 Construction process of accurate 3D model for complex landslide

1.3 滑坡外部高精度多源监测数据约束的施加

滑坡外部GNSS监测采用的是站心坐标系，因此可直接将监测形变量作为外部约束施加到模型上，分析在外部约束下滑坡的三维变形演化过程。此外，滑坡位移一般会经历初始变形阶段、匀速变形阶段和加速变形阶段。本文在FLAC3D中应将位移分时段施加，初始位移变动较小，中期位移匀速变形，后期位移处于加速阶段，形变量随着时间步的增大呈指数型增加。根据上述实际，本文根据时间步与位移的关系进一步调整施加的速度约束值。图6为施加的变形量与时间步的关系。

图6 GNSS监测点的位移追踪结果Fig.6 Displacement tracking results of GNSS monitoring points

2 附加外部约束下融合多影响因素的滑坡稳定性评价

2.1 滑坡精细三维地质模型

基于研究区滑坡外部多源高精度监测数据及高精度DEM，可构造研究区三维地质模型(图7)。根据研究区工程地质钻孔信息、地形地貌及地层岩性并结合已有资料[36-41]，将模型共分为5层，依次为非饱和黄土层、饱和黄土层、粉质黏土层、砂砾石层和基岩，各地层参数见表1。依据区内GNSS高精度监测数据、灌溉及地下水位信息精细模型，得到具有地质信息的滑坡精准三维模型。数值模拟过程中定义各地层为连续分布的各向同性介质，采用摩尔-库伦本构模型进行计算。设定模型底部为固定边界，模型四周为垂向边界，斜坡坡面为自由边界。将计算是否收敛、剪应变增量区及塑性区是否贯通作为滑坡失稳判据。研究施加外部位移约束与未施加外部位移约束下融合工程地质资料的黄土滑坡稳定性。

图7 滑坡三维模型Fig.7 Three dimensional landslide model

表1 黑方台党川滑坡强度参数

2.2 附加外部约束下融合多影响因素的黄土滑坡稳定性评价

针对工程地质中分析滑坡稳定性时没有结合滑坡外部形变监测信息，导致模拟结果与实际不符的突出问题，本节首先将HF06/HF07 GNSS监测点位移作为先验约束施加到模型上，研究附加外部高精度监测信息下融合灌溉、孔隙水压力等工程地质资料的黄土滑坡稳定性。其次，将已失稳区域(HF06/HF07)视为空模型，并将HF09 GNSS监测点位移作为约束评价未施加约束与施加约束条件下滑坡变形。

本文运用分阶段弹塑性求解法获取滑坡的初始地应力场。依据模型分组分别设置强度参数、杨氏模型和泊松比，平衡初始应力。当位移达到10-6m左右，表明初始地应力场已经形成。令初始地下水位为黏土层上方14 m，将HF06/HF07 GNSS监测点位移作为约束条件施加到模型上，HF06/HF07 GNSS监测点于2019年10月5日失稳，与初始模型相差5个月。此外，据资料显示，研究区每年进行春灌、冬灌及中间的5次苗灌，且降雨较少但蒸发作用较强[42]，气候干燥，本文在此仅考虑灌溉对滑坡稳定性的影响。综上，依据研究区灌溉量及灌溉强度[43-44]，将模型后缘灌溉30 d，停止灌溉4个月，评价未施加和施加GNSS高精度监测数据约束的滑坡稳定性，结果如图8—图11所示。

图9为两种方案获取的斜坡塑性区发育云图。未施加外部高精度监测数据约束时，坡体在重力及灌溉水作用下位移较小，仅坡体后缘及裂缝处有较小的位移(图8)。塑性区发育较缓慢。滑坡上部受到重力作用，沿裂缝处下滑，在裂缝附近拉张作用强烈。上部土体向下挤压，导致坡体底部向下变形剪切滑床，因此剪切作用强烈。综上，坡体未出现明显塑性区贯通区，较稳定。相反，在施加外部高精度监测信息下，模型位移与实际监测位移一致，此时坡体HF06/HF07附近区域处于剪切应力状态，且剪切活性塑性区已贯通。由此可见，在施加外部高精度监测数据约束下，坡体处于不稳定状态，沿裂缝处下滑。

图8 滑坡位移云图Fig.8 Landslide displacement nephogram

图9 滑坡塑性区发育云图及剖面图Fig.9 Landslide plastic zone development nephogram and its profile

图10为两种方案下的最大剪应变增量云图及剖面图。在施加与滑坡实际失稳时位移相同的约束时，监测点HF06/HF07区域最大剪应变增量较大，呈贯通之势。除约束区域外，取坡体中央剖面，得到剖面处位移及速度云图(图11)。由图11可知，相对于未施加约束方案，在约束区域外，本文方法获取的位移值虽然较小，但速度较大，达到10-4数量级，且具有沿着裂缝向坡脚处下滑的趋势。综上所述，相对于当前工程地质的数值模拟法，本文方法获取的结果表明滑坡经过30 d灌溉后处于不稳定状态，其中HF06/HF07区域最不稳定，首先发生失稳破坏，与实际监测结果更吻合。

图10 滑坡最大剪应变增量云图及剖面图Fig.10 Landslide maximum shear strain increment and its profile

图11 HF06/07中央处位移及速度剖面图Fig.11 The profile of displacement and velocity in the central of HF06 and HF07

为进一步验证模型的可靠性，将HF06/HF07 GNSS监测点所在区域视为空模型，模拟HF09 GNSS监测点区域斜坡变形发育。依据研究区监测资料，由于HF09处发生滑坡于2021年1月27日失稳，与HF06/HF07 GNSS监测点失稳相差1年零3个月，研究区地下水位每年抬升0.3～0.4 m，因此模拟时将模型地下水位黄土层14 m抬升至14.5 m，评价未施加约束与施加约束条件下滑坡稳定性，结果如图12—图14所示。

图12 滑坡位移云图Fig.12 Landslide displacement nephogram

图13 滑坡塑性区发育云图及剖面图Fig.13 Landslide maximum shear strain increment and its profile

由图14可知，相对于现有方法，在施加与实际滑坡失稳时位移相同约束的条件下，滑坡速度云图、剪应变增量云图和塑性区发育云图均显示HF09所在区域更不稳定，且约束区域塑性区剖面云图显示约束区域活性塑性区已贯通，速度达到10-3量级，计算不收敛。上述结果均进一步证实了本文方法优于现有工程地质数值模拟法。

图14 滑坡最大剪应变增量云图及剖面图Fig.14 Landslide maximum shear strain increment and its profile

表2给出了两种方案下滑坡稳定性评价结果。可以看出，相对于常规滑坡稳定性模型，本文方法获取的结果得出HF06/HF07 GNSS监测点首先失稳，其次是HF09。综上，在施加滑坡GNSS高精度监测数据约束下，党川滑坡发育失稳次序为HF06/HF07、HF09、HF05，而在未施加外部约束条件下，失稳次序为HF06/HF07、HF05、HF09。依据实际监测资料可知，长安大学基于自主研发的北斗滑坡卫士(BeiDou landslide security,BOLS)分别于2019年10月5日成功预警了HF06/HF07 GNSS监测点的失稳和2021年7月8日HF09 GNSS监测点的滑动。上述结果证实，本文提出的融合地下水、工程地质钻孔信息及灌溉资料等工程地质资料和GNSS高精度监测数据约束的黄土滑坡稳定性方法获取的结果与滑坡的实际失稳情况相吻合，均优于已有工程地质数值模拟法。综上，本文方法将滑坡外部动态高精度监测数据与工程地质数值模拟方法耦合，明显改善了滑坡稳定性评价结果的可靠性，有效地弥补了大地测量的位移时序法无法准确判断滑坡是否失稳的不足及工程地质的数值模拟法无法评价滑坡随时间变化的动态稳定性的局限。

表2 约束与未约束滑坡稳定性模型下滑坡稳定性评价结果对比统计

2.3 附加外部约束下滑坡安全系数

为进一步表明所施加约束的有效性，运用有限元强度折减法计算融合地下水、工程地质钻孔信息及灌溉资料等工程地质资料和GNSS高精度观测数据约束的黄土滑坡安全系数云图(图15)。可以看出，本文方法获取的监测点HF06/HF07 GNSS监测点所在位置安全系数为1.075，HF09 GNSS监测点处安全系数为1.1，HF05 GNSS监测点处安全系数为1.125，稳定性HF06/HF07

图15 滑坡安全系数云图Fig.15 The nephogram for landslide factor of safety

2.4 讨 论

滑坡外部动态高精度监测信息是滑坡内部变形的重要外部表征，能够实时获取滑坡随时间变化的规律，但其无法解释滑坡发生变形的原因；而多源异构数据能够提供较准确的滑坡内部物理力学参数(如地下水、工程地质钻孔信息)以及外部形状特征(如高精度DEM)，同时依据灌溉、降雨等监测资料能够合理解释滑坡滑移灾变的机理，但其无法描述出滑坡随时间变化的动态稳定性。因此，二者融合能够实现滑坡外部形变信息与内部变形机制的有效衔接，更准确地评价滑坡变形发展过程及失稳机理。基于此，本文利用甘肃黑方台党川区域已有的工程地质资料，包括地下水、工程地质钻孔信息及灌溉资料等以及高精度多源监测数据(GNSS监测数据、裂缝计观测资料等)，通过大地测量高精度GNSS观测信息与工程地质数值耦合建模这一独特设计，建立了一种优势互补的黄土滑坡稳定性综合评价方法，基于多学科交叉知识获取了相较于已有滑坡稳定性判识更加准确的黄土滑坡稳定性评价结果，为提高黄土滑坡预警的准确性和有效防治提供了有价值的理论与实际参考。

需要指出的是，作为一种融合多源异构信息和GNSS高精度观测数据的黄土滑坡综合稳定性评价方法，构建模型所需监测数据种类较多，缺少了某一类观测数据均可能会影响到模型的稳定性评定精度。例如，如缺少高精度DEM或工程地质钻孔信息，会导致无法建立具有现势性的滑坡精细三维地质模型，进而难以准确地评价滑坡稳定性。此外，本文模型在黄土滑坡稳定性评价中的适用性较好，对于其他类型滑坡的适用性还有待进一步深入研究，这也是下一步需要重点研究的工作。

3 结 论

对于复杂黄土滑坡，现有的黄土滑坡稳定性评价方法仅从滑坡地表形变或纯数值模拟的角度分析滑坡变形发育过程及所处的稳定状态，并没有将滑坡外部形变信息与内部变形机理相融合，降低了滑坡稳定性结果的可靠性。针对此问题，本文基于甘肃黑方台黄土滑坡实例，验证了融合地下水、工程地质钻孔信息及灌溉资料等工程地质资料和GNSS高精度观测数据约束的黄土滑坡稳定性评价方法，详细分析了滑坡外部高精度监测信息约束方法对已有黄土滑坡稳定性评价结果的改善程度，得到以下结论。

(1) 构造的融合高精度DEM、地下水位资料及工程地质钻孔信息的滑坡精细三维地质模型，可为科学合理地评价滑坡稳定性提供基础支撑。

(2) 附加高精度监测信息约束的滑坡稳定性综合评价模型，相较于现有工程地质数值模拟法，具有更高的精度及可靠性。

(3) 大地测量手段和工程地质数值模拟手段的有机融合，可获取与实际滑坡发育一致的结果，即：HF06/07 GNSS监测点首先失稳，其次是HF09 GNSS监测点，最后是HF05 GNSS监测点，进而实现了滑坡外部形变信息与内部变形机制的有效衔接。

综上，本文方法获取的结果与实际监测情况符合，滑坡失稳顺序与实际滑坡发育顺序一致，稳定性评价结果更可靠。本文将融合地下水、工程地质钻孔信息及灌溉资料等工程地质资料和GNSS高精度监测信息约束的黄土滑坡稳定性评价模型成功地应用于甘肃黑方台党川黄土滑坡的建模与机理讨论中，为黄土滑坡稳定性评价提供了新的思路与视角，加深了对黄土滑坡失稳动因和灾变机理变的认识与理解。