基于GIS新技术的滑坡稳定性三维评价及滑动过程模拟
谢1江崎哲郎1郭宇N周1三谷康弘1。
2翻译罗2校对朱如烈2校对
(1九州大学环境系统研究所，日本福冈东Ku箱崎6-10-1；2中国地质调查局水文地质工程地质技术方法研究所，河北保定，071051)
在传统边坡稳定性三维分析模型的基础上，提出了一种基于GIS的边坡稳定性三维网格分析模型。该模型假设初始滑动面为椭球底，在计算最小安全系数法的同时，采用蒙特卡罗随机模拟法确定最危险滑动面。当用GIS网格模型和GIS数据模拟滑坡滑动过程时，滑坡会沿着主滑动方向滑动，直到其安全系数上升到1。所有的计算都可以通过一个名为三维边坡地理信息系统(3DSLOPGIS)的计算程序来完成，该程序主要利用了GIS的空间数据处理和分析功能。
确定性模型；地理信息系统；蒙特卡罗模拟；滑动模拟；三维边坡稳定性
1简介
滑坡失稳和风险评估不仅成为地质学家和工程专家感兴趣的主要课题，也成为世界各国政府部门和管理者关注的焦点。据统计，全球每年约有600人死于山体滑坡。在许多发展中国家，自然灾害造成的经济损失占国民生产总值的1% ~ 2%。
近年来，地理信息系统(GIS)以其强大的空间数据处理功能在自然灾害评估领域得到了广泛应用。GIS是由硬件和软件组成的系统，可以实现数据的采集、输入、操作、转换、可视化、组合、查询、分析、建模和输出。GIS具有强大的空间数据分析和处理功能。同时，基于GIS的地质技术分析模型可以简单有效地分析滑坡稳定性。目前，边坡稳定性分析已广泛应用于土木工程和地质工程中。
我们通常认为传统模型既适用于均质滑坡，也适用于非均质滑坡。稳定性指标是一种广泛使用的基于岩土工程模型和物理力学参数的安全系数。安全系数的计算需要几何数据、抗剪强度数据和孔隙水压力数据，正确的结果依赖于可靠的数据和合适的模型。虽然输入数据会极大地影响安全系数，但可靠的确定性模型对于获得可靠的结果更为重要。确定性计算可以在GIS系统中进行，也可以通过其他程序进行。如果使用其他程序进行计算，GIS只是作为一个空间数据库来存储、显示和更新输入数据。这种方法的主要优点是利用外部模型计算可以节省时间；但是，它的缺点是对从外部模型获得的数据进行转换很复杂。因为每个程序都有自己的数据格式和数据结构，所以数据转换成为一个主要问题。有些程序的输入模块只允许手工输入数据。只有当这些程序的默认数据格式为ASCII码时，才能直接进行数据转换。使用外部模型的另一个缺点是计算结果通常以点或线的形式表示，而不是GIS的空间分布方式。因此，改变这个计算结果的表达方式也是一个主要问题。
用于计算安全系数稳定性模型的边坡是二维或三维的。由于一个地区包括许多边坡，每个边坡必须单独分析，使用这些模型计算安全系数的空间分布需要花费大量时间。要克服数据转换的困难，可以利用GIS内部确定性计算模型来实现。但这种方法也有缺点，即由于常规二维GIS中复杂算法、迭代过程和三维体的复杂限制，只能容易地实现简单模型。目前只有基于GIS的无限边坡模型可以单独计算每个像元的安全系数。研究表明，只有使用越来越成熟的三维模型和GIS系统，才能彻底解决此类问题。
从近年来GIS在边坡稳定性分析中的应用研究中发现，大多数研究者致力于用统计方法来确定边坡失稳与影响因素之间的关系。虽然GIS可以准备和处理区域数据，但只有少数研究人员使用GIS的集成功能和边坡稳定性的确定性模型。
即使在很短的距离内，边坡破坏也有其不同的空间几何结构。因此，采用三维模型分析边坡稳定性是合理的。自20世纪70年代中期以来，三维稳定模型的发展和应用越来越受到重视。地质力学著作中提到了几种三维分析方法。
上面提到的大多数方法都使用直方图方法。这些方法忽略了圆柱体之间的力，或者作为三维安全系数计算的假设。因为所有与斜坡相关的GIS数据都可以转换为栅格数据，所以这些基于三维模型的柱可以用于通过使用GIS栅格数据来计算三维稳定性。然而，长期以来，人们习惯用著名的“一维模型”——“无限边坡”模型来描述滑动面平行于地面的长期自然边坡的潜在危险。该模型只能用于浅层边坡失稳分析和一些深层滑坡的区域研究。
由于2D地理信息系统中算法的复杂性、步骤的重复性和表达三维数据的困难性，三维确定性模型的应用在早期文献中并未提及。为了克服GIS数据外部转换的困难和GIS内部算法的复杂性，本研究在GIS软件组件中使用了Visual Basic程序。三维因子的计算和滑动过程的模拟由计算机中的三维边坡地理信息系统(3-DSOPGIS)的计算程序完成。在该系统中，GIS组件(ESRI生产的MapObjects2.1)可以完成所需的GIS功能，就像普通的GIS软件一样，可以有效地管理和分析所有与滑动相关的数据。用于计算三维边坡安全系数的数据全部采用GIS数据格式(如矢量和栅格数据层)，因此不需要在GIS数据格式和其他程序的数据格式之间进行数据转换；同时，复杂算法和三维问题的交互程序也可以理想地实现。
本研究将GIS网格数据与基于直方图的三维边坡稳定性分析模型(Hovland，1977)相结合，推导出一种新的基于GIS网格的三维确定性分析模型。
采用蒙特卡罗随机模拟法寻找最小安全系数，从而确定临界滑动条件。假设基本滑动面为椭球体的下部，临界滑动受不同地层和不连续界面条件的影响。这种客观事物的变化导致了最小三维安全系数。
如果滑坡的三维安全系数小于1，则滑坡处于危险状态，因此评价滑坡灾害的规模和影响范围至关重要。因此，在本研究中，利用基于GIS的三维网格数据模型和GIS网格数据模拟滑坡滑动过程的目的是评价滑坡风险和预测其影响范围。
2基于GIS的三维模型
利用GIS的空间分析功能，输入所有与计算三维安全系数有关的数据(如高程、倾角、坡度、地下水、地层、滑动面和力学参数等。)都有其对应的网格元素，所有与坡度相关的数据都被栅格化。当这些数据输入到确定的边坡稳定性模型中时，可以计算安全系数值。在Hovland模型的基础上，下面详细介绍基于GIS的三维模型。在该模型中，考虑到孔隙地下水压力，所有输入数据可以简单地转换成栅格数据。
图1是具有潜在滑动面的滑动体的三维几何示意图。滑坡的稳定性与地质地层、地貌、地质力学参数和水动力条件有关。
图1边坡崩塌的三维景观
图2显示了用于研究该物质的土壤(或岩石)小柱的离散性。所有与滑坡有关的数据都可以用一个柱状的三维视图来表示，如图2所示。假设每个圆柱单元的垂直面为无摩擦面(圆柱单元的垂直面不受其他边界的影响，或者其影响可以忽略)，三维安全系数可以用公式(1)表示:
地质灾害调查与监测技术方法论文集
其中:F3-D为三维边坡的安全系数，W为一根柱子的重量，A为滑动面的面积，C为粘聚力，φ为内摩擦角，θ为滑动面的角度，J和I为边坡破坏范围网格内的行数和列数。如果没有GIS，基于圆柱模型的三维安全系数计算将是一项冗长费时的工作，更新和增加数据也极其不便。而在GIS中，利用GIS的空间数据处理和分析功能，整个研究区的边坡稳定性相关数据可以用矢量图层来描述，如图3所示。对于每个图层，可以通过GIS空间数据处理和分析功能获取栅格数据，并根据精度要求确定像素大小。
图2滑动面和三维棚直方图
现在，基于栅格数据将边坡破坏划分为圆柱体。参考图2，诸如地表、地层、地下水、裂缝和滑动面的空间数据都可以从栅格数据层获得。因为与坡度相关的数据量非常大，所以不可能高效地管理所有栅格数据集。因此，在三维斜坡地理信息系统中，有一个点数据库专门用来存储这些栅格数据，其中有一个属性表来链接所有与滑动相关的数据。每个网格直方图的中心点设置为点类型，其他区域设置为与滑坡相关的一些数据(如地面高程、地层和裂隙的高程、地下水、滑动面深度等。).表1是属性表的一个示例。
图3用于边坡稳定性分析的GIS图层
表1点数据库示例描述
另一方面，为了控制滑坡边界，有效地管理和分析空间数据，将滑坡边界定义为多边形文件。
基于这个点数据库，可以将公式1变成基于GIS的方程。这里所有的阻力和滑动力都是沿着滑动方向的，不一定像Hovland模型中使用的那样沿着Y轴。在本研究中，假设斜坡区域的主要倾斜方向是可能的滑动方向。根据图4，滑动表面积可以通过公式(2)获得。
地质灾害调查与监测技术方法论文集
下面的公式来自图4:
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
然后，X轴和Y轴的倾斜角推导如下:
地质灾害调查与监测技术方法论文集
写出α = cellsize/cos θ xz，b=cellsize/cosθyz，则一个网格直方图的滑动表面积为:
地质灾害调查与监测技术方法论文集
滑坡范围内主滑动方向倾角的计算公式如下:
地质灾害调查与监测技术方法论文集
此时，三维边坡在水平滑动方向的安全系数可通过以下公式计算:
地质灾害调查与监测技术方法论文集
图4三维安全系数推导公式的网格直方图。
这里，对于每个网格，ZJI和ZJI是表面高程和滑动面高程，uji是滑动面上的孔隙水压力，γ′是单位重量。
为了检验基于网格的GIS三维稳定性分析模型，我们利用该模型做了实例计算。示例问题是具有球形滑动面的均质粘土滑坡，其他参数如图5所示。图5中，c为粘聚力，φ为摩擦角，r为瞬时摩擦力，γ为土的容重。采用闭合形式算法，三维安全系数为1.402。克拉拉模型计算的安全系数为1.422。同样的问题，利用三维边坡模型计算的三维安全系数在1.386 ~ 1.472之间，取决于分离边坡所用的柱数。
图5示例问题验证
采用基于GIS网格的三维稳定性分析模型(图5)，设置网格尺寸为0.5m，计算三维安全系数为1.386；当网格尺寸为0.6m时，计算的安全系数为1.388。显然，与封闭算法相比，基于格网模型的GIS可以有效地用于三维边坡稳定性评价。
3临界滑动面的确定和蒙特卡罗模拟
滑动面只能通过岩土工程勘察来确定，而由于地质勘察费用较高，通常很难确定滑动面。因此，边坡稳定性评价对确定临界滑动面非常重要。
为了确定三维临界滑动情况，采用蒙特卡罗随机模拟法计算三维安全系数的最小值。假设初始滑动面是一个椭球体的下部，那么边坡面就根据不同地层的应力和不连续的界面条件而变化。最后得到危险滑动面，并得到相关三维安全系数的最小值。
4椭圆坐标变换
假设初始滑动面为椭球体的下部，椭球体的倾斜方向设置为与研究区的主倾斜方向一致；椭圆的倾斜度基本上被设置为接近研究区域中波动的倾斜度。主倾角为α，主倾角为β，由边坡破坏区域的主网格像素值决定。假设倾角和倾角属于正态分布，将主倾角α和倾角β代入分布模型:
地质灾害调查与监测技术方法论文集
利用公式(10)和(11)完成坐标转换。图6示出了坐标转换过程。
图6坐标转换过程
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
其中:x，y，z为全局大地坐标，局部坐标，x0，y0，z0为椭球中心点坐标。
滑动面5 Z值和倾角的确定
滑动面上“b”点的z值由公式(12)根据直线AB和椭圆的计算结果确定(见图7)。
地质灾害调查与监测技术方法论文集
对于每一个网格像素，滑动面的倾斜度和倾角可由下式计算，像素(j，I)的倾斜度可由图8中点1 ~ 4的z值确定。点1 ~ 4的值由公式(13)(14)(15)计算，滑动面的倾斜度和倾角由公式(16)计算。
图7确定了滑动表面上的z值。
图8滑动倾角的计算
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
这里，Z(j，I)是像素(j，I)的Z值，θ是倾斜角，β0是相对于X轴的倾斜度。在GIS中，倾角是与Y轴的夹角。因此，最高点为点3时，倾向为90°-β0°；最高点为点4时，倾向为90+β0；最高点为点2时，趋势为270-β0；最高点1时，趋势为270+β 0。
6随机模拟
为了确定临界滑动面，通常采用蒙特卡罗模拟法来选择三维边坡稳定性分析的变量。这些变量是椭球的中心点、几何参数和倾角。作为研究区域的中心，首先需要确定椭球的中心点，然后在一定范围内随机选取。
椭球的几何参数A、B、C由用户在一定范围内随机设定，确定的范围如公式(17)所示:
地质灾害调查与监测技术方法论文集
假设A、B、C为均匀分布，蒙特卡罗模拟的随机变量由公式(18)和(19)计算。
在[0，1]范围内均匀分布的随机变量可以通过同余乘积法得到:
地质灾害调查与监测技术方法论文集
地质灾害调查与监测技术方法论文集
其中:ri是平均分布在[0，1]范围内的随机变量。均匀分布在[a，b]范围内的随机变量可以用公式(19)计算。
地质灾害调查与监测技术方法论文集
其中:xi是均匀分布在[a，b]范围内的随机变量。
椭球倾角设为平均分布的随机变量。平均分布范围是主倾角及其在一定波动范围内变化的变量。
7三维安全系数最小值的计算过程
整个研究区域(或边坡破坏范围)可以划分为若干个小的矩形网格，就像基于网格的GIS一样。关于基于网格的三维边坡稳定性分析的数值计算，所有的计算过程都可以通过前面提到的Visual Basic(使用GIS组件)来完成。该软件名为三维边坡地理信息系统，是用Visual Basic 6.0和ESRI公司的MapObjects 2.1开发的。作为GIS的一个组成部分，MapObjects用于组织和分析GIS数据。计算三维安全系数的过程如图9所示。
图9三维安全系数最小值的计算过程
在这个过程中，数据模块的作用是获取所有与边坡有关的地质、地貌、水动力数据和地质力学参数；随机可变参数模块用于随机选择蒙特卡罗模拟的实验滑动面；三维边坡稳定性模块可用于计算三维安全系数；危险滑动面及其安全系数可以通过一些试验来计算。从图9可以看出，所有关于GIS空间分析功能的模块都可以通过GIS组件来实现。由于在三维边坡地理信息系统中完成了GIS组件，可以有效地计算三维安全系数；同时，利用与边坡相关的GIS数据，所有相关的数据和结果都可以在三维边坡地理信息系统中可视化。
示例剖面如图10所示。本例考虑的因素有:4地层、地下水和滑动面；其物理机械参数见表2。
表2研究实例的物理和地质力学参数
图10断层面研究实例
图11计算次数和最小三维安全系数实验
为了确定临界滑动面，进行了蒙特卡罗随机计算次数试验，总* * *计算次数达到1000次。每次实验计算的最小三维安全系数结果如图11所示。图中清楚地显示了经过300次实验计算后得到的最小安全系数。这300次实验的结果如图12所示，这些计算结果相差不大，最小值为1.34，最大值为1.68。临界滑动的研究方案是基于最小安全系数的计算。最小安全系数的计算结果取决于参数的随机选择。参见图13，获取该临界滑动示例的三维视图。通过比较三维模型和二维模型的结果，发现图10所示的二维模型和表2中所列的参数是用Janbu法确定临界滑动面的。这个二维模型计算的安全系数为1.18，略小于三维模型计算的最小值(1.346)。
图12三维安全系数分布曲线
滑坡滑动过程模拟
基于GIS栅格三维边坡稳定性分析模型和GIS栅格数据，模拟滑坡滑动过程，直至三维安全系数大于1。滑动方向根据滑动面的主滑动方向确定。图14显示了由滑动面确定的八个滑动方向。比如滑动面方向的倾角在22.5° ~ 67.5°之间，滑坡将要滑动的方向正好在图的右上方(即“5”方向)。
图13临界滑动面三维展开图
图14滑动面主要滑动趋势
图15滑坡滑动过程模拟框图。
滑坡滑动过程模拟流程见图15。首先要计算滑坡在初始状态下的三维安全系数，以确定其滑动的可能性。如果安全系数小于1，则进行下一次滑动过程模拟。首先，沿滑动面主倾角确定的滑动方向移动滑坡多边形；然后在新的滑坡多边形范围内，逐步计算每个网格的DEM和滑动变化(每一步等于一个网格的大小)，再次计算下一个滑坡的新的滑动方向。基于新的DEM数据和滑动多边形范围，计算新的三维安全系数。如果三维安全系数仍然小于1，则执行以下新的滑动步长模拟。
在这种滑动模拟模型中，假设滑动面的内摩擦角不变，但假设除了在初始三维边坡安全系数的计算过程中，滑动面没有内聚力(即内聚力为零)。
还是用同一个例子(如图5所示)，用两个不同的动力学参数模拟滑坡滑动过程:
案例1:c=4kN/m2，φ = 110，y = 23kn/m3。
情况二:c = 6kn/m2，φ = 10.5，γ = 23kn/m3。
第一种情况，初始边坡安全系数为0.82，经过7步滑动，滑坡体趋于稳定，其安全系数为1.04。图16显示了一些滑动步骤的截面和三维视图变化。在该图中，DEM的变化和滑坡的移动过程一目了然。利用三维边坡地理信息系统，还可以将可视化的滑动过程表示为GIS地图和剖面图。滑坡沿水平方向的最终滑动距离为3.0m
图16不同滑动阶段地表和剖面的三维视图
第二种情况，滑坡会一路滑下到平地，水平滑动距离为14m。滑坡最后停止滑动位置的三维视图如图17所示。
图17滑坡最后停止位置
9讨论和结论
在边坡稳定性三维矿柱分析模型的基础上，提出了一种新的基于GIS网格的三维确定性模型，并通过实例验证了其正确性。在三维边坡稳定性分析模型中，假设其初始滑动面为椭球体；采用蒙特卡罗随机模拟法确定其三维临界滑动面，得到最小三维安全系数。基于GIS的网格三维模型，模拟滑坡滑动过程，用于判断滑坡灾害和预测滑动距离。开发了三维边坡地理信息系统作为计算程序软件，可以完成三维边坡问题的所有计算，并利用GIS组件实现GIS的空间分析功能和有效数据的管理。三维边坡稳定性问题因其具有空间分析、数据管理和综合边坡数据GIS可视化的优势，研究起来相对容易。一种全新的基于GIS网格的三维边坡稳定性分析模型问世以来，为习惯于用传统数学方法研究边坡稳定性的工作者拓展了一个新的研究领域和数据库方法。