基于GIS的隧道洞口边坡的稳定性分析

魏新江,张 阳,陈浙江,韩同春

(1. 浙江大学 建筑工程学院，浙江 杭州 310058；2. 浙江大学 城市学院土木工程学系，浙江 杭州 310015；3. 金华市公路管理局，浙江 金华 321013； 4.金华市交通投资集团有限公司， 浙江 金华 321015)

1 概述

在我国高等级公路的建设中，由于隧道可以克服高程障碍，缩短路线长度，减小坡度和曲率，从而提高技术标准，隧道往往成为选线最佳方案，尤其是在地形复杂的山区公路中，从而出现了大量隧道。大量的既有隧道，随着时间的增加会产生很多问题。其中隧道洞口段围岩埋深浅、稳定性差、地质构造条件复杂，且容易受到偏压地形影响，导致洞口段易受到滑坡、剥落、崩塌等边坡变形破坏问题，对公路的安全通行产生影响，对隧道洞口段灾变研究显得十分重要[1-3]。2016年4月25日，S210浦江段杭口岭隧道右洞顶部山体发现一系列张拉裂缝，裂缝最大宽度达2 m，可见深度约6 m，延伸长度达到了30 m，给隧道右洞造成了安全隐患[4]；2012年9月18日大广高速江西龙南境内一个在建隧道洞口发生塌方事故，致使16名正在施工的工人被困在隧道内；2016年7月2日晚11时50分左右，重庆奉溪高速杨家湾隧道奉节端洞口发生滑坡，左右线洞口被滑落的泥石封住，造成交通中断。这些事故阻塞交通，往往造成巨大的经济损失，甚至毁坏车辆，威胁人民的生命安全。因此，有针对性地加强隧道洞口部位的灾变机理分析，制定有效的加固方案，具有重要的意义。

地理信息系统(简称GIS)，是在计算机硬、软件系统支持下，对整个或部分地球表层空间中的有关地理分布数据进行采集、储存、管理、运算、分析、显示和描述的技术系统。针对GIS发展至今所拥有的各项强大功能，许多学者开始在岩土各个领域中利用GIS技术，特别是辅助地质体的建模、工程勘察数据的管理、地质结构的再现和边坡稳定性评价等方向，取得了不少成果[5-7]。由于在空间分析方面的明显优势，GIS在复杂边坡的稳定性分析方面应用更是越来越广泛。早在20世纪90年代，美国、意大利和澳大利亚等国家就已经开始在边坡的地质灾害监测中使用GIS技术[8-10]。相对来说我国的起步较晚，但国内的学者凭借近10 a快速发展的拟真技术和图形算法，也取得了许多先进的研究成果。如谢谟文[11]等基于4个采用了GIS栅格数据的边坡稳定三维极限平衡分析模型，开发用于边坡稳定性评价的GIS 扩展模块。王纯祥[12]等基于GIS和数值模型相结合的方法，采用两步骤预测和评价日本某地区的滑坡和泥石流灾害。赵春宏[13]等利用GIS建立的三维极限平衡分析模型，不仅算出边坡的三维安全系数，还得到了滑动力的空间分布规律。

这些在复杂边坡中应用GIS技术的研究成果，主要是借助GIS空间数据分析能力来完成边坡的稳定性评价，但仍存在不足，主要是模型理论的普适性差，同样的建模方法不一定适用于其他的边坡，特别是对地形起伏大、土层构成复杂的边坡，本文针对隧道洞口边坡，借助于GIS进行三维边坡建模，进行稳定性分析，以提高隧道洞口边坡数值模拟的易用性和建模效率。

2 分析方案

2.1 地理信息系统

GIS发展到现在，出现了大量的专业性软件，功能强大且价格昂贵。与价格昂贵的商用软件相比，开源GIS软件除界面不友好外，功能和模块一致，同时其开源性，便于用户进行二次开发。本文采用开源GIS中的GRASS，即Geographic Resources Analysis Support System，属于开源地理空间基金会的一个正式项目，生态发展良好，更新频繁，用户可获取源代码，根据具体需求增删功能，已在世界各地的学术和商业领域得到广泛应用[14]。现阶段地表及地层的监测数据主要以离散点的形式表示，需要人工推断或插值拟合才能得到场地内其他位置的信息。GRASSS内置了多个理论完善、易于拓展的插值模块，建模过程中可印证不同插值方法的效果，选取最符合实际的插值理论。对于边坡而言，降雨及地下水的情况对稳定性也有着较大的影响，但在水文的三维模型方面，现阶段岩土范畴内的数值模拟软件尚处于发展阶段，若引入完善的GIS水文分析模块，能更好地模拟实际情况。因此，选择GRASS代替数值模拟软件完成三维边坡建模的工作，可降低建模难度，提高最终模型的精度。

2.2 稳定性分析框架

在整个技术方案中，GIS是实现边坡建模和模型管理的主要工具，具有强大的地理数据处理能力，并带有一定的可视化功能，能对已建模型进行剖切、划分、叠加等操作。数值模拟软件则是专注于计算，对已建模型进行深入的分析，完成如边坡稳定性分析、有限元渗流分析、有限元应力分析、地震动力响应的仿真模拟。当二者耦合，既可充分发挥GIS的空间数据分析能力，又可利用岩土领域内数值模拟软件在边坡计算方面的经验，提高边坡工程中数值模拟方法的易用性。但在耦合过程中，存在有数据格式不一致，图形界面操作过多，数据传递耗时长，以及原始数据、模型配置高度离散化等关键难题，需要通过文件结构解析和程序语言等其他手段解决。

本文主要利用Python语言解决GIS与数值模拟软件之间的数据传递和格式转换等关键问题。选择Python作为接口工具主要有两大优势：一是GRASS的核心和图形界面系统是采用Python语言编写，可在shell窗口中调用GRASS的核件[14]。一旦设置完关键参数，用户可在后台完成GRASS的调用，删改和微调模型时可减少耗时长的图形界面操作。二是Python具有大量拓展功能，拥有科学计算、数据库链接、图形绘制等组件，以及与某些格式对应的转换模块，能够以一门语言来完成大部分程序耦合工作。但作为一门动态解释型的语言，Python的编译过程长，多线程计算优化差，进行三维建模计算时效率不满足实际应用要求。因此，Python可承担特定剖面的地质数据计算，主要的建模流程需由GRASS完成。在实际应用中，Python调用的模块主要是将点集和面数据转换成标准dxf文件的dxfwrite模块，和用于预览图像和绘制图表的matplotlib模块。Python还能调用SQL Server数据库封存每个模型及剖面的数据，调用视觉化工具函式库VTK渲染模型的三维图像。

数值模拟软件，如Abaqus、Plaxis、SolidWorks、Geostudio等，能读取inp、dxf、xml、json等以文本格式存储的模型文件。所以在耦合GIS与数值模拟软件的过程中，Python可利用正则表达式和数据库等工具，解析相应文件格式，实现inp、dxf、xml、json等模型文件的自动读取和写入，减少耦合中图形界面的繁琐操作，节省大量数据模拟软件在正常使用中的建模耗时。至此，如图1的方案总体结构所示，Python通过组件的拓展，以粘合剂的方式实现GIS与数值模拟软件无缝耦合的技术方案，大大增加了岩土实际工程中数值模拟方法的易用性和建模效率。

图1 方案总体结构图Figure 1 Overall framework of the scheme

3 实现步骤

3.1 数据导入

实现三维数据建模的第一步就是要导入场地的数据，包括地表的高程数据以及地层的分布数据。目前岩土工程勘察中地形图测绘多采用GPS与全站仪相结合的方法，这些数据通常包含在工程地质平面图和包含土层信息的钻孔图中呈现，格式多为dxf和dwg。在边坡建模的过程中，各项原始的勘探数据必须经过处理才能批量处理。针对dxf格式的图纸，可利用方案自主编写的程序，直接读取高程控制点的高程值以及等高线的多段线数据。针对dwg格式的图纸，可利用AutoCAD的数据提取功能，批量输出控制点的x、y坐标及高程值。

3.2 模型建立

根据地表和地层的数据，建立相应的三维模型。在实际工程中的边坡模型，一般场地范围较小，而高程变化范围大，唯一能揭露土层信息的钻孔数量通常是有限的。在这种信息量不足，难以组成一定密度的三角网的情况下，栅格结构模型可通过定义场地内所有位置的地表及土层分界线高程，降低复杂边坡的模拟难度，还可与遥感数据联动，本文采用栅格数据结构建立边坡三维模型。

在GIS中，栅格数据结构的数字高程模型是以像元为基本单位，索引参照矩阵，采用行列号进行定位，以西北角作为起点，从西往东、从北向南依次增大。导入并栅格化高程控制点后，用插值或模拟等方法，给场地内空值像元赋值。空间模型的插值方法有很多，常用的有最邻近点插值法、距离反比加权法、最小曲率插值法、Kriging法和径向基函数插值法等[16]。本文采用径向基函数插值法。

除了地表模型、每个土层的分界面和地下水模型外，一个完整的边坡模型还应包含各模型之间完整的拓扑关系和对应位置，见图2。导入原始数据并完成每个部分的建模后，GRASS可按照土层年代、现场具体情况等因素将各个模型进行组合，输出最终的三维边坡模型。

3.3 模型输出及稳定性分析

模型建立后，就可以进行边坡稳定性分析。为此需要在模型中选择计算剖面，将此剖面的数据几何模型数据、材料参数、拓扑关系等信息导入数值模拟软件，完成边坡的仿真分析和模拟计算。

图2 三维边坡模型示意图Figure 2 Schematic of three-dimensional slope model

4 工程实例

浙江金华杭口岭隧道位于S210省道浦江段杭坪至浦江城区方向K52+550 m处，是连接浦江县城与西北乡镇以及桐庐、建德等县市的交通要道，日交通量达到2万辆。在公路安全巡查中发现洞口隧道南向左洞顶部山体发现了系列张拉裂缝，危及道路行车安全。险情发生后，勘察部门对洞口附近进行了补充勘探。勘察人员利用GPS和全站仪，测量了860个控制点高程，并结合初步测量成果及裂缝位置，在场地内打了6个钻孔，以便探明土层构成。勘探提供的地质资料有地质平面图、地质断面图及6张地质钻孔图，希望用这些数据建立边坡模型，为紧接下来的边坡处理方案提供支撑。此边坡位于隧道入口和道路附近，人工干预多，地形变化大，土层分布不均匀，难以找到具有代表性的断面，与二维模型相比，建立三维的边坡模型更能反映实际情况。为此进行了场地内的三维建模，并利用三维模型实现场地稳定性分析工作，为加固方案的设计的提供支撑。

4.1 边坡表面数据输入建模

4.1.1数据提取

将场地内的高程点数据导入到GIS系统中，是建模首先要完成的工作。在抢险阶段，提供有地质平面图，图中包含有场地内控制点的平面坐标及高程。为此我们采用AutoCAD的数据提取功能，将每个控制点的坐标及高程值等相应信息以统一格式提取，过程如下：在打开工程地质平面图后，点击插入，提取数据，出现数据提取对话框，点击下一步后，弹出“将数据提取另存为”对话框。给文件取名后，按“保存”，后在“数据提取-选择对象”对话框中取消勾选你不需要的对象，接着单击下一步，在“数据提取-选择特性”中勾选你需要的特性，需要注意的是，对于GIS建模仅需要X和Y坐标及高程值，其他数据不需要选择，如图3所示。点击下一步后，最后生成的数据以EXCEL表格的形式保存。

将提取出来的数据通过GRASS数据导入命令导入GRASS，生成土层表面的栅格图。栅格图是将地理空间分割成有规则的网格，在每个网格上给出相应的属性来表示地理实体，在项目中每个格栅点有一个空间坐标。在生成过程中，可以指定格栅的精度，格栅精度与采集到的数据量，以及问题的类型有关，论文中指定格栅精度为0.5 m×0.5 m。

4.1.2数据插值

数据导入生成的地表格栅图，由于数据量的限制，一般情况下格栅图形内会有大量空值，也就是有一部分网格没有被赋值。 这时就需要进行数据插值，插值的过程实际上就是要从已有的数据中找到一个函数关系式，这个关系式可以很好地逼近已有的空间数据，然后再根据这个函数关系式推求出区域范围内其他未知点的值。插值方法的选取直接影响到边坡模型的精度和模拟结果的可信性，其中散点数据拟合的众多插值理论中，空间模型运用最广泛的是以Gauss分布函数为基函数的Kriging法，及以样条函数为基函数的径向基函数插值法[18]。本文使用RST(Regularized Spline with Tension)径向基函数法对散乱的高程控制点进行插值，以建立完整的三维模型。RST径向基函数插值法采用完全规则样条函数(CRS)作为基底函数，通过限制基函数的光滑半范数取值，可得到二元问题中有明确形式的插值公式：

图3 数据提取Figure 3 Data extraction

(1)

式中：常数a1为选定的趋势函数T(x，y)；N为采样点的数量；ρj=(φrj/2)2，其中rj为带插值点到已知点距离，φ为张力参数；E1为指数积分函数；CE为欧拉常数。

当基函数已知，就可解出式(1)中待求解系数a1、{λj}，求解方程组为：

z(xi，yi)=

(i=1，2，…，N)

(3)

(4)

(5)

式中:ω0/ωj为每个样本点带有权重的光滑因子，用于消除对控制点所在格栅插值时可能产生的误差。

根据上述原理和GRASS内自带的径向基函数插值模块，笔者编译了RST径向基函数法的GRASS建模用插件，只需选定合适的张力参数，即可迅速完成场地的建模，建成的地表模型如图4所示。

图4 地表模型图Figure 4 Surface model of slope

4.1.3模型验证

为了验证建成的地表模型与实际的地形间的相近性，将地质平面图中的等高线[图5(a)]和用GRASS求解已建模型的等高线图[图5(b)]进行对比，可以看到，已建边坡三维模型的地形整体走势与实际地貌一致，山谷线和山脊线的分布也相同。

(a) 实际等高线

(b) 模型等高线

分析图5可以看出来，当小范围内的坡度变化较小时，模拟值与监测值十分接近，但当坡度发生骤变时，误差会相应增大，绝对误差大于1 m的4个样本点均是位于坡度骤变的地方。以上验证和分析结果说明，GIS在进行模型的插值拟合时，为保证表面整体平滑度，会牺牲部分坡度变化剧烈地区的模拟精度。但如图5的中部地区，样本点附近的控制点数量较多时，即使坡度变化大，该点的误差也在平均范围内。因此，为减少此类拟合误差，勘察阶段应在高程骤变区域的边缘，如断崖的顶部与底部、公路沿线两旁等地方适当增加高程控制点，提高边坡模型精度。

4.2 土层建模

场地内共有6个钻孔，一共揭露了5个土层。从地表向下，土层的名称及参数见表1。

表1 土层材料参数Table 1 Physico-mechanical parameters of stratum土层名称密度ρ/(g·cm-3)粘聚力c/MPa内摩擦角φ/( °)含碎石粉质黏土1.8528.117.5强风化泥岩1.9515.628.0含碎石块石粉质黏土1.8936.720.3全风化晶屑凝灰岩2.303.636.0中风化晶屑凝灰岩2.504.239.0

按顺序对不同土层分别建立相应的模型，首先是表层土含碎石粉质黏土层。由于场地内仅有6个钻孔，导致用于土层插值运算的点很少，土层建模质量差。为此假定原有的860个地表的控制点为虚拟钻孔，先用原有钻孔的数据进行插值计算，获得每个虚拟钻孔的土层信息。利用经加密的钻孔数据，可求得每个空值像元的具体值。建立起来的含碎石黏土层的维模型如图6所示。

图6 含碎石粉质黏土层模型Figure 6 Model of silty clay layer with gravel

应用相同的方法建立其他土层的模型，将各个模型合并到一起，形成边坡模型，如图7所示。可以看到强风化泥岩和含碎石粉质黏土在场地内是局部分布的，其中强风化泥岩可以在模型北侧看到。

图7 边坡模型Figure 7 Three-dimensional slope model

4.3 边坡稳定性分析

进行边坡稳定性分析，需要对相应的剖面进行稳定性计算。在整个场地范围内进行稳定性分析，工作量巨大，也没有必要。对于边坡稳定分析来说，场地范围内在土层条件变化不大的情况下，坡角越大，发生边坡失稳的可能性越高，这里借助GRASS GIS的坡角计算功能对场地进行识别，划出坡角较大的区域，也就是可能失稳的范围进行计算。

在GRASS GIS中采用拟合曲面法来拟合局部的地形表面，采用3×3的网格来计算中心网格点Z0的坡度，如图8所示。

图8 坡度计算图Figure 8 Chart for slope calculation

中心网格点的坡度按下式计算：

(6)

式中:φx、φy分别是x和y方向的坡度，计算式如下:

(7)

根据上述计算过程，在GRASS系统中计算得到的坡角分布如图9所示。

图9 坡角轮廓线图Figure 9 Contour map of slope angle

从图9可以看出，场地内的最大坡角可达70°左右，主要分布在洞口外侧道路东侧，为此取隧道洞口道东侧边坡进行稳定性分析，共取5个断面，位置如图10所示。

图10 剖面位置图Figure 10 Scheme of profile locations

为分析这些剖面位置的边坡稳定性，根据起止坐标，在GRASS GIS中获得相应剖面的地表高程、土层分界面高程及土层信息等数据，其中剖面3在图形界面下所建计算模型，见图11。

图11 计算剖面Figure 11 Calculation section

为完成剖面的稳定性计算，需要将剖面数据导入到稳定性分析软件Geostudio中，手工操作需要进行繁琐的工作。为增加计算的效率，论文将剖面数据用编好的接口转换为xml文件，再用Geostudio直接打开用坐标剖切边坡，到完成完整的计算模型，全过程所耗时间极短，在Geostudio中只需要指定滑移面进入范围和退出范围外，其他步骤皆可在后台完成。到这一步，本文的技术方案快速实现了此边坡模型精度高、结果可靠的稳定性分析，剖面的计算结果见表2。

表2 剖面安全系数汇总Table 2 Summary of section safety factor剖面号安全系数Bishop法Janbu法Morgenstern-Price法11.2841.2131.26221.0321.0211.02530.9350.9130.92540.9880.9770.99551.1251.1121.119

采用了不同计算方法计算剖面的稳定性安全系数，从计算结果来看，Bishop法计算得到的安全系数最大，Morgenstern-Price法结果次之，Janbu法结果最小。同一剖面中，3种方法的临界滑移面分布相同或相近，但相较于Bishop法和Janbu法，Morgenstern-Price法的计算结果能满足水平力、竖向力和力矩的平衡，安全系数计算结果更可靠。因此在这部分山体的整体稳定性评价中，可以将Morgenstern-Price法作为主要的稳定性分析方法，余下方法作为参考。除了能对特定位置进行分析，本实例所建的边坡模型还可深度挖掘已有数据，如求解雨水流动路径、计算准确的土方量等功能，可以进一步开发。

5 结论

a.本文以隧道洞口地段边坡为研究背景，基于GIS和Geostudio等软件的耦合，建立了一套全过程通用性高、适用性强的三维边坡建模和模拟计算的技术方案，为具体工程提供了实体模型和理论计算方面的支持。通过杭口岭隧道工程边坡实例的证明，此方案在运行过程中建模部分理论完整，满足复杂边坡建模的各项要求，在边坡分析部分计算效率高，分析方法灵活性强，还利用自行编写的程序，真正实现从控制点数据处理，到边坡三维建模，到稳定性分析的无缝连接。此外，全过程的快速完成，也说明此方案拥有为出现险情的边坡提供模拟计算支持的能力。

b.利用Python语言编写的一系列插件和程序，能实现地理信息系统GRASS GIS和稳定性分析软件Geostudio间的数据传递，实现系统间的高效连接。

c.利用GRASS GIS的边坡坡角分析功能，自动寻找危险区域，为有针对性稳定性分析提供指导。

d.论文通过三维边坡的建模计算，结合实例，模型的分析建议在前期勘察阶段，应在地形骤变区域的边缘，如断崖的顶部与底部、公路沿线等位置适当增加控制点，以提高后期模型的精度。

e.本文在三维边坡模型的数值计算方面开展了一些探索性工作，计算模型仅是在只在2.5维阶段进行稳定性分析，并没有建立可用的真三维极限平衡模型,需要在以后的研究中继续完善。