基于DEM的三维微条柱法计算滑移型塌岸安全系数的研究

王小东，马潇菲，袁广祥

(华北水利水电大学， 河南 郑州 450011)

目前边坡稳定性分析，一般简化为二维平面应变问题并采用二维极限平衡法来处理，二维条分法是对边坡实际情况简化后采用的方法，实际上，任何边坡失稳现象均是在现实的三维空间中发生的，因此进行三维计算与分析是研究的方向，即将边坡置于三维空间之中，在传统的二维条分法的基础上进行拓展，二维条分法中的土条也相应地拓展为三维条柱法中的条柱，通过分析每个条柱上的受力状况，根据极限平衡理论，最终计算得到三维安全系数。三维极限平衡法发展至今，主要的理论模型有三维简化毕肖普法、三维简化简布法、三维斯宾塞法以及三维Morgenstern-Price法等[1-6]。三维极限平衡法计算结果与二维极限平衡法的计算相比，大多偏于安全[7]，也正因此，三维极限平衡法需要进一步研究、发展和完善。二维极限平衡分析法的计算机实现程序已相对成熟，Simon等[8]基于excel的VBA开发模式，融合了三种极限平衡计算公式，实现安全系数的计算；一些商业化的软件如Geoslope等广泛应用于边坡稳定性计算，并被广大的研究者和工程师所接受。在三维极限平衡分析法的计算机程序实现中，Reid等[9]在研究构造火山侧翼边坡三维稳定性分析中，采用三维毕肖普法与DEM(Digital Elevation Model，简称DEM)相结合，进行边坡稳定性评价；姜清辉等[10]通过多层DEM构造三维地质模型，基于Visual C++可视化开发平台，结合OpenGL，开发了三维边坡稳定性极限平衡分析软件slope3D。简化毕肖普法是毕肖普提出的边坡稳定性计算方法，在工程中得到广泛的应用，其计算方法简单且计算结果有较高的精度。Hungr于1987年在简化毕肖普法的基础上，将该方法在三维上进行了直接的拓展，即三维毕肖普法。其假设条件与二维状态时完全一样[1，11]。此后，许多研究者也做过类似的研究[12-13]，就三维毕肖普法而言，其基本假设和理论基础均是一致的。

本文在Reid等[9]对构造火山侧翼进行三维稳定性研究所采用的理论和方法的基础上，结合高分辨率DEM，探讨三维毕肖普法的实现，其优势在于结合了地理信息系统(Geographical Information System，简称GIS)软件平台，能够快速实现三维模型的构建，使得三维模型构建与三维安全系数的计算更加的自动化和一体化。

1 三维简化毕肖普法的基本原理

本研究基于高分辨DEM数据，选择在二维条分法的应用领域里被工程所接受的毕肖普法进行扩展，即三维毕肖普法，进行三维安全系数的计算。图1展示在三维状态下条柱的全部受力状况[14]。

基于图1中条柱的受力状况，结合Reid等对构造火山侧翼进行三维稳定性研究所采用的理论和方法，即满足每个条柱沿竖直方向即Z轴方向的力的平衡，并根据绕X轴的整体力矩平衡求解安全系数，计算三维安全系数的公式为：

图1 条柱受力状况Fig.1 Forces acting on individual columns

(1)

其中，R为底滑面抗滑力的力臂；W为条柱重力；Ac为滑动面与条柱相交部分的面积；c为粘聚力；φ为内摩擦角；k为水平地震力加速度系数，其作用点为条柱中点，力臂为e，无地震力作用时，k=0；γu为孔隙压力比；θ为失稳条柱底部的真倾角；α为底滑面沿滑动方向的视倾角。

其中，mα的计算中包括安全系数F，因此，三维毕肖普法的计算也同样需要迭代运算才能完成，mα采用公式(2)表示：

ma=cosq+tanjsina/F

(2)

重力W的计算与二维不同，即

W=Vγt

(3)

其中，体积V通过计算微条柱的体积得到，通过计算平截头六面体的体积近似得到，即

V=1/6Δx(S0+4S1+S2)

(4)

该方法由于未考虑x轴和y轴的整体力的平衡，所以根据整体力矩平衡条件求解安全系数时与坐标轴的位置有关，所以该方法比较适合于滑裂面为弧形面的情况。

2 基于DEM的三维毕肖普法实现的方法和步骤

本文提出了基于GIS的三维极限平衡法的实现思路和方法，即参照Reid等的研究[9]，基于所获取的2.5 m×2.5 m DEM，采用GIS组件开发模式，完成三维安全系数的计算。

应用三维微条柱法实现三维稳定性分析的具体原理及步骤如下：

(1)确定滑移型塌岸边界。通过野外调查，结合遥感解译，确定涉水滑坡堆积体范围，以矢量多边形的形式在二维空间中表达，如图2(a)所示。

(2)确定失稳滑动的方向。滑移型塌岸失稳滑动方向是进行塌岸体三维条柱划分的基础。在塌岸体详细调查的基础上，对塌岸可能的失稳方向作出判断。

(3)三维微条柱划分。以DEM的空间分辨率为参考，对所要进行稳定性计算的三维岸坡进行微条柱划分，即沿滑动方向和与滑动方向垂直的方向生成一系列垂直相交的直线，覆盖整个塌岸体在二维平面的投影范围，并与塌岸体边界线进行相交运算，最终完成对整个塌岸体的条柱划分，如图2(b)所示。

图2 条柱划分过程Fig.2 The division process of individual columns

(4)滑动面的选择与确定。首先确定滑移形塌岸体岸坡性质，即判断塌岸岸段土体为粘性土还是无粘性土，这样便可确定塌岸体可能的失稳滑动模式，为选择试算滑动面提供依据。

(5)确定滑动角。滑动面确定之后，滑动角也基本能够确定，对于滑动面为平面的情况，如图3(a)所示的二维剖面，每一个条柱的滑动角即为滑动面的倾角；对于滑动面为球形或椭球形时，如图3(b)所示的二维剖面，每一个条柱的滑动角是不同的，即沿着过条柱底部中心点的切平面与水平面的夹角；对于组合滑动面的情况，如图3(c)所示的二维剖面，即不同的滑动面区域采用相应的滑动角；当然，除了图3给出的这三种情况外，滑动面还有很多种情况，但均可以通过平面与平面、平面与球面的组合方式得到。

图3 滑坡面形状Fig.3 The shape of Landslide surface

(6)土层划分与参数载入。由于考虑到水库蓄水对滑移型岸坡的影响，因此，即使针对均质岸坡而言，也被蓄水位线分割为水上和水下两个部分，因此要在竖直方向上对土体进行分层，针对不同土层采用不同的物理力学参数，按土层顶部和底部高程确定不同土层采用不同的参数，通过键盘输入或文件导入的方式载入程序，参与运算。

通过(1)—(6)步可完成对塌岸体三维地质模型的构建，下面则针对极限平衡法进行计算和分析。

(7)依据三维简化毕肖普法计算安全系数。为了计算安全系数，首先按照以上所述(1) —(6)步构建三维塌岸体地质模型。

(8)结果输出与表现。采用三维毕肖普法计算安全系数以及对应的最危险滑动面，并在三维空间上进行展示。

3 算法实现与实例分析

本研究以Visual Studio为可视化开发平台， 嵌入ArcGIS Engine(简称AE)地理信息系统二次开发组件，以高分辨率DEM数据为基础，基于三维微条柱法实现的基本方法和步骤，以三维简化毕肖普法为基础，完成了塌岸体三维地质模型的构建和三维安全系数的计算。以金沙江溪洛渡库区老木沟滑移型塌岸体为例，进行岸坡稳定性评价，过程如下：

(1)基于二维毕肖普法搜索得到临界滑动面

沿滑动方向绘制主剖面线，参照瑞典条分法实现过程中构建试算滑动面的方法，构建试算圆弧滑动面，采用毕肖普法获得最小安全系数对应的最危险滑动面的位置，并记下圆心坐标和半径，如图4所示。这一过程和瑞典条分法的实现过程一致。

图4 二维毕肖普法获得临界滑动面Fig.4 Critical slip surface obtained by 2D Bishop method

(2)基于AE组件二次开发构建球形滑动面

得到二维状态下的临界滑动面所在的圆弧的圆心和半径之后，就可通过该圆心和半径构造三维球形滑动面，其过程在AE二次开发中实现时，最终的目标是获得沿与主剖面线垂直方向按一定距离间隔平移后得到的一组剖面线在球形滑动面上投影点的三维坐标，如图5所示。

图5显示了在AE地图容器中构造球面上点的过程，图5(a)中的间隔距离即为过该剖面线所在球体截面圆与过主剖面线截面圆之间的垂直距离，将其投影到过主剖面线所在的截面圆上，如图5(a)所示，可得到一系列同心圆，这些圆的半径可由球半径和与过主剖面所在截面圆的距离通过勾股定理计算得到，如图5(b)所示。经过这一过程之后，便构造出一组同心圆弧，将其离散成点之后，计算圆弧滑动面上点的高程，就得到了该点的三维坐标，最终得到一组球面上的三维坐标点，插值即可生成三维球形滑动面。

图5 三维球形滑动面的构造过程示意图Fig.5 Construction process of 3D spherical sliding surface

图6以老木沟堆积体为例，展示了通过AE二次开发构建三维滑动面的过程，如图所示，第一步：由二维毕肖普法计算得到主剖面上二维临界圆弧滑动面；第二步：由二维滑动面构造三维球形滑动面的一组三维点坐标，三维球形滑动面在平面上的投影为一椭圆；第三步，由三维点坐标通过插值生成三维球形滑面DEM。

图6 构造三维滑动面并生成滑动面DEM表面Fig.6 Constructing 3D sliding surface and generating DEM surface of sliding surface

(3)基于三维毕肖普法计算得到三维安全系数

三维滑动面构建完成之后即可按照三维毕肖普法的计算公式，对岸坡所在区域的多边形进行微条柱划分，如图7所示，将塌岸多边形沿滑动方向进行正方形网格剖分，分别投影在原始坡面DEM和三维滑动面DEM之上，以便于进行每个条柱上的计算，根据公式(1)进行三维安全系数的计算，与二维毕肖普法一样，三维毕肖普法的实现同样是一个迭代的过程，即采用公式(1)和公式(2)完成这一过程。

图7 网格与坡面和滑动面的关系Fig.7 The relationship between mesh and original slope DEM and spherical sliding surface DEM

对于式(1)—式(4)中所包含θ，α和V等参数，在AE组件开发中可采用下列方法确定：

真倾角θ的确定：真倾角θ与真倾向有关，是指条柱底部与三维球形滑动面相交的面的中心所在的切平面与水平面的夹角，从几何意义上讲，是真倾向与其水平面上投影线的夹角，真倾向可通过对三维滑动面DEM进行坡向分析得到，如图8中坡向图所示。真倾角则可通过对三维滑动面DEM进行坡度分析得到，如图8中坡度图所示，在二次开发进行三维毕肖普法计算时，基于条柱底面中心点位置坐标和坡度图获得该条柱底面的真倾角θ。

图8 三维球形滑动面坡向图和坡度图Fig.8 Aspect map and slope map of 3D spherical sliding surface

视倾角α的确定：视倾角α与真倾向和视倾向有关，对某个条柱底面而言，存在多个视倾角，但对于三维毕肖普法计算安全系数而言，视倾角α是指沿岸坡主滑方向的视倾角，即将主滑方向作为视倾向，通过视倾向与真倾向之间的夹角ω，如图9所示，可通过公式tanθ=tanα·cosω，计算得到视倾角α。

图9 真倾角与视倾角之间的关系Fig.9 The relationship between true dip angle and apparent dip angle

最终，以金沙江溪洛渡库区老木沟滑移型塌岸体为例进行三维安全系数的计算，其在天然状态下的抗剪强度指标为：c=25 kPa，φ=40°，γ=20 KN/m3，以25 m间隔构建条柱，通过三维毕肖普法计算得到的三维安全系数为1.48，而采用二维毕肖普法计算得到安全系数为1.04。由此可见，三维毕肖普的计算结果偏于安全。引起这种差异的原因较多，其中一个原因是由滑动面的形状造成的，球形滑动面在现实中很少见到，图7也明确显示出球形滑动面在平面的投影范围没有完全覆盖实际调查中所确定的塌岸体的范围，因此，在理论上是一种可行的方法，在实际中是一种近似的方法。

4 结论

本文通过AE二次开发的方式，将三维毕肖普法引入GIS软件平台，利用GIS平台的数据处理与空间分析功能，耦合极限平衡理论计算边坡安全系数的方法，实现了三维圆弧滑动面的构造和三维安全系数的计算，主要结论如下：

1)以高分辨率DEM为结合点，耦合地理信息系统分析方法与三维极限平衡法，进行安全系数的计算。高分辨率DEM提供了详细的地表信息，十分适合边坡的三维模型构建，本文将三维毕肖普法引入GIS软件平台，实现了三维安全系数的计算，为单体稳定性评价与区域稳定性评价的结合提供了新的思路和方法。

2)提高了边坡三维模型构建与三维安全系数的计算的自动化、一体化程度。AE二次开发的方式提高了三维地质模型构建的效率，即通过图切剖面线方式构建球形滑动面，提供了边坡三维模型构建的自动化程度，并将GIS中坡度、坡向等空间分析方法与三维极限平衡法相结合，实现了三维安全系数计算的一体化。

3)三维滑动面的构造需要进一步的优化。应用本论文方法，以金沙江溪洛渡库区老木沟滑移型塌岸体为例，进行三维稳定系数的计算，与二维计算结果相比，结果偏于安全。主要的原因在于将滑动面视为球形滑动面，是一种理想的状态，实际的滑动面大多是非圆弧的，因此需要在圆弧滑动面的基础上进行优化，可在程序中加入控制性结构面或控制性节点，使其与实际情况更加接近。