基于GIS的边坡三维稳定性计算

杨 伟，丁伯阳，潘晓东，朱益军，陶海冰

（1.浙江工业大学 岩土工程研究所，浙江 杭州 310032；2.浙江省交通规划设计研究院，浙江 杭州 310014；3.杭州师范大学 理学院，浙江 杭州 310036）

基于GIS的边坡三维稳定性计算

杨 伟1，丁伯阳1，潘晓东1，朱益军2，陶海冰3

（1.浙江工业大学 岩土工程研究所，浙江 杭州 310032；2.浙江省交通规划设计研究院，浙江 杭州 310014；3.杭州师范大学 理学院，浙江 杭州 310036）

根据高速公路高边坡的钻孔资料，利用Arc Map建立进行三维边坡稳定性分析的TIN数据层和栅格数据层，将栅格数据模型的各层对应栅格块视为一个滑坡单元体，通过读取相应栅格值进行稳定性计算.针对顺层滑坡的特点，采用ArcGIS集成的VBA开发环境，提出了一种三维折线法分析边坡的稳定性.该方法通过对栅格数据的重采样，可以计算任意滑动方向三维边坡的安全系数，通过选择最小安全系数来确定边坡的主滑动方向.结合某高速公路大型岩质滑坡的工程实例，通过比较三维折线法计算结果与某典型剖面的二维不平衡推力法计算结果，讨论该方法对工程滑坡治理的指导意义.

三维边坡；边坡稳定；GIS；滑动方向；折线法

GIS（地理信息系统）软件在土木工程中的应用越来越广泛，在公路边坡这种涉及大量空间信息存储、处理和分析的工程实践中，更加体现出其优越性.目前，滑坡稳定性分析常用二维的极限平衡法，对某一可能滑体，采用计算几个不同剖面的稳定系数来综合评价.虽然出现了很多三维极限平衡法，多以柱体单元为计算单元，由于计算方法繁琐，在二维GIS中较难实现，通常只借助GIS作为空间信息存储和处理的工具，通过数据的转换实现在其他软件中的三维稳定性计算.对于边坡三维问题，谢漠文结合GIS栅格数据和4个边坡稳定极限平衡模型开发了一个用于边坡三维稳定分析的GIS扩展模块，提供了一种便利的数据处理方法和分析方法［1］.

利用ArcGIS的VBA开发环境，采用工程中广泛应用的折线法，进行滑坡三维稳定性计算，可以方便地确定滑坡体的三维稳定系数，并对滑坡体的假定的可能滑动方向进行试算，找出最有可能的滑动方向.

1 滑坡体模型的建立

滑坡体的DEM模型和各土层的TIN数据、栅格数据，可以根据钻孔资料，通过合适的插值方法得到.各土层的上下关系通过栅格法进行空间拓扑关系的处理［2］.ArcGIS提供的空间插值方法有反距离权重插值、样条法插值和克里格插值.无论选择哪种插值方法，样本点越多，样本点分布越广，插值结果越接近实际值［3］.而在实际工程中，钻孔点的分布可能达不到计算精度的要求，这时可以引入虚拟钻孔来提高精度.

顺层滑坡的滑动面，通常可以通过钻孔资料得到破碎滑动带来确定，或分别对各地质层层底作为潜在滑面，分别做稳定计算而得到.图1为边坡DEM模型，图2为以破碎带为滑动面得到的滑坡体.

图1 边坡DEM模型Fig.1 A slope digital elevation model

2 滑坡体三维稳定计算

图2 滑坡体Fig.2 A landslide mass

对滑坡体进行柱状单元划分，由于滑坡体基于栅格数据，事实上可以看成是以栅格大小为面积，以地面到滑面距离为高的柱状体.与滑坡体稳定性分析相关的数据，如：各土层底面高度，水位线高度，滑动面倾角，滑动面倾向等，存储在相应的栅格数据中［4］.在ArcGIS的VBA开发环境中，可以在滑坡区域内提取相应的栅格数据进行计算.

假定滑体的主滑动方向为X轴，Y轴垂直滑动方向.每一个栅格单元都与X轴垂直.栅格柱体间的条间力可分解为水平向的剪力Q、沿X轴方向的推力Ei和Y轴方向的挤压力Eh，对于每一个滑体的受力情况如图3所示.

图3 栅格柱状体受力图Fig.3 Force on a cell column

采用折线法计算，忽略滑块水平方向的剪力的作用，假定Q＝0，滑块只传递主滑动方向的推力，不考虑挤压力的影响，则基于栅格数据的滑体向主滑动方向的推力［5－6］为

其中：ψi－1＝cos（θxzi－1－θxzi）－sin（θxzi－1－θxzi）tgφi／Fs3D 为传递系数；θA为方位角 ，θxz为滑面与X轴的夹角；Wi为各滑块重力；Ai为各滑面面积.可以计算出X轴滑动方向的整体的稳定系数为

其中：Ri＝（W cosθi·sinφi＋c Aicosφi）sinθA；Ti＝W sinθi·sinθA.当所有滑体的φi值相同，则：Ri＝（W cosθi·tgφi＋c Ai）sinθA；Ti＝W sinθi·sinθA.式中各项的值可从栅格数据中取值得到.

3 计算机程序实现

ArcGIS下的VBA开发环境中，提供了众多栅格数据的接口和方法.对于基于同一坐标的土层栅格数据集合，可以直接读取栅格值进入计算程序，不必进行数据转换就在Arc Map中实现了边坡三维稳定性计算.与边坡计算有关的各土层的栅格数据可以从基于钻孔数据而来的TIN数据得到.这样得到的栅格数据层在边界处会有栅格不重合的现象，即在同一坐标点上，有可能只有地面高程数据，而没有滑动面高程数据.在计算程序中，边界上无法同时取到地面数据和滑动面数据的点不读入数据.为了计算简便，对土层参数取加权平均值，作为整个栅格柱体的参数，将栅格柱体作为单一土层考虑.

求解安全系数时，按以下步骤：1）首先假定安全系数K＝1，读取第一行第一列的栅格数据，计算得到R1，T1；2）依次读取下一列的栅格数据，计算对应的传递系数ψi，得到的R′i＝Ri－1ψi＋Ri，T′i＝Ti－1ψi＋Ti继续往下传递；3）逐行、逐列计算，按式（2）得到新的安全系数K′；4）当｜K′－K｜＞0.01，则K＝K′，重新进行计算，直到K值满足要求.计算框架图如图4所示.

图4 稳定系数计算框图Fig.4 Calculating diagram of stability cofficent of landslide

4 滑动方向的确定

确定滑坡滑动方向的传统方法是根据地表和滑动面形态，通过经验方法估计确定，缺乏定量基础，基于栅格数据的边坡稳定性计算公式中引入了倾角、倾向等参数，则为滑坡体滑动方向的确定提供了条件.选取不同的滑动方向计算出对应的边坡的安全系数，则安全系数最小的某个滑动方向，为最可能的滑动方向［7］.利用ArcGIS可以较方便的确定主滑动方向.在ArcGIS下，栅格数据可以通过Rotate工具重采样得到新的栅格数据，提供的重采样方法有三种：最邻近采样（NEAREST），双线性（BILINEAR）和三次卷积采样（CUBIC）.最邻近采样精度较低，对旋转后得到的栅格数据进行坡度分析和坡向分析误差较大.为提高内插精度，选用三次卷积采样，取待计算点周围相邻的16个点，可先在某一方向上内插，如先在X方向上，每四个值依次内插四次，再根据四次的计算结果在Y方上内插，最终得到内插结果.对旋转后的滑动面数据进行坡度分析和坡向分析，得到旋转后的坡度角和方位角的栅格数据层.此时得到的栅格单元与原坐标系X轴垂直，可直接套用上式计算X轴滑动方向的安全系数，不必进行坐标变换.利用上面提到的折线法计算程序，假定某一方向为主滑动方向，只需通过栅格数据的旋转，即可以计算出不同主滑动方向的安全系数.

5 工程实例

杭金衢K103高边坡是一个大型的岩质滑坡，该滑体体积约为160万m3，属于工程施工诱发的、潜伏型路堑滑坡.根据地质分析和定量评价，该滑坡治理前长期处于缓慢蠕滑状态，经历公路开挖后4年多缓慢蠕动变形，2005年春季，滑坡体后缘及两侧均发现多条裂缝，坡脚挡墙出现拉裂、挡墙沉降缝错位、墙面与坡面局部鼓起等现象，相应路段高速公路路面也出现了不同程度的裂缝开裂与路面鼓包现象.滑坡方向大致为正东向，与公路延伸方向近于正交.滑坡体的厚度一般在15～40 m之间，局部大于40 m.滑坡岩体破碎，风化强烈，呈碎块石夹泥状，节理面风化后大多呈泥状，潜在滑面主要沿破碎带发展，滑面分布次生夹泥，强度低.根据现场地质勘察结论，该滑坡滑面前缓后陡，目前还处于蠕滑变形阶段，按照规范要求，对其取稳定系数K＝0.98，利用Bishop法对该滑面进行二维的反演分析，得到滑面参数的加权平均值为：c＝20.5，φ＝14.8.

利用高程数据生成地面TIN数据，再将其转换为栅格数据.利用钻孔中的软弱夹层位置数据与滑坡裂隙编辑数据生成滑面TIN，也将其转换为栅格数据.对滑动面进行坡度分析和坡向分析，得到滑动面倾角栅格图和滑动面方位角栅格图，见图5.利用VBA开发的折线法计算程序计算安全系数，当c＝20.5，φ＝14.8时，正东方向的安全系数为：K＝1.152.对倾向进行统计分析，多分布在60°～120°范围内，在此范围内假定不同主滑动方向所得到安全系数如表1所示.

图5 坡度角和方位角Fig.5 Slope angle and aspect angle

表1 安全系数表Table 1 Slope Stability coefficient on different direction

由此推断，最危险滑动方向在倾向为110°，相对应的安全系数为1.096.如果在计算程序中考虑地下水的作用，安全系数可能会进一步降低.该结果与实际边坡的情况较吻合，有一定的参考价值.

用不平衡推力法对滑坡进行稳定性计算，该法在水利部门、铁路部门广泛应用.选取滑坡上某一典型剖面，用二维不平衡推力法计算安全系数，通过调整安全系数的大小，最终使滑坡出口处的剩余下滑推力为0，此时得到的安全系数为滑坡的稳定系数.同时计算地下水位变化时，滑坡的稳定系数变化情况［8］，见表2.

表2 地下水位变化时滑坡安全系数表Table 2 Slope stability coefficient on different groundwater level

可见，在完全疏干的情况下，该典型剖面的安全系数取值为1.26，与计算程序得到的滑动方向正东10°范围内的结果比较接近，但并非最小值.而最危险滑动方向的确定，对缓慢发展的工程滑坡治理具有一定的指导意义.

6 结 论

基于GIS栅格层数据模型能将各层对应栅格块视为一个滑坡单元体，极大的方便了边坡三维稳定性的计算.在ArcGIS环境下，结合VBA二次开发技术，能直接读取ArcGIS中栅格数据的值进行稳定性计算，避免了复杂的数据转换工作.基于栅格数据的三维折线法边坡稳定性计算公式，引入了倾角、倾向等参数，对定量确定滑动方向提供了条件，对不同滑动方向计算所得到的安全系数最小的方向为最可能的滑动方向.不同滑动方向的安全系数相差最大的可到30%，可见滑动方向的选取具有一定的工程意义.边坡三维稳定性计算模型对不同土层的参数取用的是加权平均值，相当于作为单一土层考虑，且未考虑地下水作用，有待改进.确定滑动方向时需要对较大范围进行试算，计算量较大.可以进一步探究边坡的倾向与安全系数之间的解析式，来确定可能的滑动方向，这是下一步要继续深入的工作.

［1］谢漠文，江崎哲郎，周国云.基于边坡单元的三维滑坡灾害评价的 GIS方法［J］.岩石力学与工程学报，2003，22（6）：969-976.

［2］陶海冰，蒋蓓.基于数据的栅格法模拟三维土层研究［J］.杭州师范大学学报，2007，7（6）：469-473.

［3］吴秀芹.ArcGIS 9地理信息系统应用与实践［M］.北京：清华大学出版社，2007.

［4］谢漠文，蔡美峰，江崎哲郎.基于GIS边坡稳定三维极限平衡方法的开发及应用［J］.岩土力学，2006，27（1）：117-122.

［5］陈斗勇.折线形滑面边坡安全系数的一种新的求算法［J］.水文地质工程地质，1991，18（6）：46-47.

［6］殷宗泽.土工原理［M］.北京：中国水利水电出版社，2007.

［7］李先华，管群，杨武年，等.数字地形图上滑坡动态方向稳定性系数的计算与滑坡滑动方向的确定［J］.岩石力学与工程学报，1999，18（3）：308-311.

［8］方鹏飞，朱益军，朱向荣.杭金衢高速公路K103滑坡治理［J］.工程勘察，2009（1）：31-35.

Three-dimensional slope stability analysis based on GIS

YANG Wei1，DING Bo-yang1，PAN Xiao-dong1，ZHU Yi-jun2，TAO Hai-bin3
（1.Geotechnical Engineering Institute，Zhejiang University of Technology，Hangzhou 310032，China；2.Institute of Communications，Planning，Design ＆Research，Hangzhou 310014，China；3.College of Science，Hangzhou Nomal University，Hangzhou 310036，China）

According to the borehole date of Highway，the TIN data layer and grid data layer for three-dimensional slope stability analysis were carried out by use of Arc Map.The value of each gird cell of different layers are picked up as slide block，it makes convenience in the slope stability analysis.Considering the property of bedding landslide，the three-dimensional slope stability analysis program based on VBA of ArcGIS is put forward in this paper.The program can calculate the Stability factor of any directions by resampling of raster date，thus the real sliding direction can be determined by the minimumValue of the Stability factor.Taking one of the large rocky landslide as an example，the results of the three-dimensional with the two-dimensional method of a typical profile of the slope are compared.

three-dimensional slope；slope stability analysis；GIS；sliding direction；line method

TU457

A

1006-4303（2012）01-0092-04

2010-10-12

浙江省科技厅基金资助项目（2009C33047）

杨 伟（1982—），男，浙江金华人，硕士研究生，主要从事岩土工程数字化研究，E-mail：airball150336＠126.com.

（

刘 岩）