运用地理信息系统新技术进行滑坡稳定性三维评价和滑动过程模拟研究.pdf

1 运用地理信息系统新技术进行滑坡稳定性三维评价运用地理信息系统新技术进行滑坡稳定性三维评价 ①① 和滑动过程模拟研究和滑动过程模拟研究 Mowen Xie 1 Tetsturo Esaki 1 Guoyun Zhou1 Yasuhiro Mitani 1 著 张晓娟 2 译 罗靖筠2 校 朱汝烈2 复校 （1 Environmental System Insitute,Kyushu University,Hakozaki 6-10-1,Higashi Ku,Fukuoka,Japan; 2中国地质调查局水文地质工程地质技术方法研究所，河北保定，071051） [摘 要摘 要] 本文在传统的边坡稳定性三维分析模型的基础上,提出了一个全新的基于 GIS 的边 坡稳定性三维栅格分析模型。 在这个模型中， 假定初始滑动面就是椭球底面, 采用 蒙特卡洛（Monte-Carlo）随机模拟方法，在求取最小安全系数法的同时,确定出最 危险滑动面。运用 GIS 栅格模型和 GIS 数据模拟滑坡滑动过程时，滑坡体将沿主 滑方向滑动，直到其安全系数上升到 1 为止。所有的计算均可通过一个称为三维 边坡地理信息系统（3DSLOPGIS）的计算程序来完成，该程序主要利用 GIS 的空 间数据处理分析功能。 [关键词关键词] 确定性模型 地理信息系统GIS 蒙特卡洛 （Monte-Carlo）模拟 滑动模拟 三 维边坡稳定性。 1 引言 1 引言 滑坡不稳定性和风险评价不但已成为地学家和工程专家们感兴趣的主要课题， 同时也成 了世界各地政府部门和管理者关注的焦点。据统计世界上每年约有 600 人葬身于滑坡灾害 中。在许多发展中国家，自然灾害所带来的经济损失，占总国民生产总值的 12。 近年来， 由于地理信息系统具有强大的空间数据处理功能， 被广泛运用于自然灾害评价 领域。GIS 是由硬件和软件组成的系统，它可以实现数据采集、输入、操作、转换、可视化、 组合、质疑、分析、建模和输出等过程。GIS 对空间数据具有强大的分析和处理功能。同时， 基于 GIS 的地质技术分析模型，可以简便而有效的分析滑坡稳定性。目前已经被广泛地用 于土木工程和地质工程中， 进行边坡稳定性的分析。 我们通常认为一个传统的模型无论是对 均质滑坡还是非均质滑动都是适用的。 稳定性指数是被广泛应用的、 基于岩土工程模型和物 理力学参数的安全系数。安全系数的计算需要几何数据、剪切强度数据及孔隙水压力数据， ① 译自 Environment Grology,200343503512 2 正确的结果取决于可靠的数据和恰当的模型。尽管输入的数据会较大程度地影响安全系数， 但一个可靠的确定性模型对于取得可靠结果则更为重要。确定性计算可在 GIS 系统内执行， 也可利用其它程序完成。若使用其它程序计算，则 GIS 只作为一个空间数据库用来存储、 显示、更新输入数据。此方法主要优点是利用外部模型计算可以节约时间；而其缺陷是对从 外部模型获得的数据进行转化时较为复杂。 因为每一个程序都有其自己的数据格式和数据结 构，数据转换成为一个主要的问题。有些程序的输入模块只允许人工输入数据。只有当这些 程序所默认的数据格式都是 ASCII 码时，数据转换才可直接进行。运用外部模型的另一个 缺点是计算结果通常不是按 GIS 的空间分布模式来表达，而是以点或线的形式表述的。因 此，改变这种计算结果的形式也是个主要的问题。 用来计算安全系数稳定性模型的边坡是二维或三维的。 因为一个地区包括很多边坡， 而 且必须分别对每个边坡做分析，所以利用这些模型计算安全系数的空间分布非常花费时间。 要克服数据转换的困难，可以利用 GIS 内部确定性计算模型来实现。然而这一方法也有缺 点，那就是由于应用复杂算法，迭代过程及在常规二维 GIS 中的三维体积等复杂局限性， 使得只有简单的模型能较容易实现。当前，只有基于 GIS 的无限边坡模型能分别计算出每 个像元的安全系数。 研究表明， 只有当越来越多的成熟的三维模型和 GIS 系统得到使用后， 才能彻底解决这类问题。 从近来对 GIS 用于边坡稳定性分析的调查中发现，大部分研究者潜心于运用统计学方 法来确定边坡破坏与影响因素之间的关系。尽管 GIS 能对区域数据进行了准备和处理，但 是只有极少量的研究者运用了 GIS 的集成功能和边坡稳定性的确定性模型。 即使在很短的距离范围内，边坡破坏在空间上都有其不同的几何结构。因而，运用三 维模型分析边坡稳定性是合理的。 从二十世纪七十年代中期以来， 三维稳定性模型的发展和 运用日益受到关注。在地质力学的著作中提到了几个三维分析方法。 上面提到的大部分方法都用到了柱状图法。这些方法将柱体之间的作用力，或者说作 为三维安全系数计算的假定前提，都忽略不计。因为所有与斜坡相关的 GIS 数据都可转成 栅格数据，所以这些基于三维模型的柱体，就可能用来借助于使用 GIS 栅格数据，进行三 维稳定性的计算。然而，长期以来大家习惯采用人尽皆知的“一维模型”“无限斜坡”模 型， 来描述滑动面与地面平行的长期天然边坡的潜在危险性。 这样的模型仅仅可以用于浅层 斜坡失稳分析和一些存在深层滑坡的区域性研究。 由于算法复杂、步骤重复和三维数据在二维 GIS 中难于表达，早期的文献中并没有提 及三维确定模型的应用。为了克服 GIS 数据的外部转换和 GIS 内部算法复杂等困难，此次 3 研究中，在 GIS 软件组件（a GIS component）中使用了 Visual Basic 程序。三维因子的计算 和滑动过程的模拟由计算机内的三维边坡地理信息系统（3-DSLOPGIS）的计算程序完成。 在这个系统中，GIS 组件（ESRI 公司生产的 MapObjects2.1）可以完成所需的 GIS 功能，就 像普通的 GIS 软件一样，它可以有效的管理和分析所有与滑动相关的数据。所有用来计算 三维斜坡安全系数的数据都采用 GIS 的数据格式（例如矢量和栅格数据层） ，因此，没必要 在 GIS 数据格式和其它程序的数据格式之间进行数据转换；同时，复杂算法和三维问题的 交互程序也可以理想的实现。 在此次研究中， 将基于 GIS 栅格数据和基于柱状图的三维边坡稳定性分析模型相结合 （Hovland，1977） ，演绎了一个新的基于 GIS 栅格的三维确定性分析模型。 运用蒙特卡洛随机模拟方法求最小安全系数值，从而确定临界滑动条件。假定基本滑 动面是一椭球体的较低部分， 临界滑动则受不同地层受力情况和不连续界面状况的影响而变 化。客观事物的这种变化引出最小三维安全系数。 如果滑坡的三维安全系数小于 1，滑坡就有滑动的危险，那么评估滑坡灾害的规模和 影响范围是非常重要的。因此，在此研究中，采用基于 GIS 三维栅格数据模型和 GIS 栅格 数据来模拟滑坡滑动过程的目的，就是评估滑坡危险性和预测其影响范围。 2 基于基于 GIS 的三维模型的三维模型 利用 GIS 的空间分析功能，所有与三维安全系数计算有关的输入数据（如高程、倾向、 斜坡、地下水、地层、滑动面和力学参数等）都有其对应的栅格元，而所有与斜坡相关的数 据都是栅格化的。 当这些数据输入到确定的边坡稳定性模型中时， 就可计算出一个安全系数 值。在 Hovland 模型的基础上，下面详细介绍基于 GIS 的三维模型，在这个模型中，考虑 了孔隙地下水压力，所有输入数据都能简单地转换成栅格数据。 图 1 是具有潜在滑动面的滑体的三维几何示意图。滑坡的稳定性与地质岩层、地貌、 地质力学参数和水动力条件有关。 4 图 2 所示是土壤（或岩石）小柱状研究体物质的离散性。所有与滑坡相关的数据都可 用如图 2 所示的柱状三维可视图来表示。假定每一个柱体单元的垂面均为无摩擦面（柱 体单元的垂面不受其它边界影响， 或其影响可忽略不计） ， 三维安全系数可用公式 （1） 表示 F D−3 ∑∑ ∑∑ JI JI W cA θ φθ sin tancos （1） 式中F D−3 为三维斜坡安全系数，W 为一个柱体的重量，A 为滑动面面积，c 为内聚 力，φ 为内摩擦角，θ为滑动面的角度，而 J、I 为在斜坡破坏范围栅格内的行列数和柱体 数。如果没有 GIS，则基于柱体模型的三维安全系数的计算将是冗长且耗时的工作，数据的 更新和增加也极其不便。然而，在 GIS 中，通过运用 GIS 空间数据处理与分析功能，整个 研究区的边坡稳定性相关数据，可用如图 3 所示的矢量图层来描述；而对于每一层，则可通 过 GIS 空间数据处理与分析功能得到栅格数据，其像元大小可根据精度需要而定。 5 现在，将斜坡破坏划分为基于栅格数据的柱体。参考图 2，诸如地表、地层、地下水、 裂缝和滑动面之类的空间数据均可从栅格数据层中得到。因为与斜坡相关的数据量非常大， 所以不能高效的管理所有的栅格数据集。因此，在三维边坡地理信息系统中，有一个专门储 存这些栅格数据的点数据库，其中，有一个属性表用来链接所有与滑动相关的数据。每个栅 格柱状图的中心点设置点类型，其它区域则设置与滑坡相关的一些数据（例如地面高程、 地层和裂缝的高程、地下水、滑动面的深度等等） 。表 1 所示即是属性表的一个实例 表 1 点数据库的实例描述 编 号 类 型 地面高度 （ m） 坡 向 （ ○ 坡 度 （ ○ ） 地层 1 高程 （ m） 地层 2 高程 （ m） 破坏滑动面深度 （ m） 地 下 水 头 （ m） 23号 点 栅 格 中 心 184 . 90 249 .77 6 . 04 地 形 参 数 170 . 00 164 . 22 139 .76 168. 80 地 质 和 水 动 力 资 料 另一方面，为了控制滑坡边界和有效管理空间数据并进行分析，滑坡的边界线被定义 为多边形类型文件。 基于这种点数据库，公式 1 可以改成基于 GIS 的方程。这里所有的阻力和滑力都是沿 着滑动方向的，而不必如 Hovland 的模型所用的 Y 轴方向。在本研究中，假定斜坡区域的 主要倾斜方向为可能滑动方向。根据图 4，滑动表面面积可由公式（2）得到。 Aabsinα 2 从图 4 推导出如下公式 cos/.tan Aspgdgcθ 3 yzxz θθαsinsincos 4 6 接着，x 和 y 轴的倾角推导如下 sintantan.costantanAspAsp xzyz θθθθ 5 记 xz cellsizeaθcos/和 yz cellsizebθcos/，则一个栅格柱状图的滑动表面面积为         − yzxz yzxz cellsizeA θθ θθ coscos sinsin1 22 2 6 滑坡范围主滑动方向的倾角计算公式如下 costantanAvrAspAsp Avr −θθ 7 至此，三维边坡水平滑动方向安全系数可以用下面的公式计算 [] ∑∑ ∑∑ − −− − JI AvrAvrjiji JI Avrjijiji D zZ uzZcA F θθγ φθγ cossin costancos 3 8 这里，对于每个栅格，Zji，zji分别为地表高程和滑动面高程，uji为在滑动面上的孔隙 水压力，而 γ为单位重量。 为了检验基于栅格的 GIS 三维稳定分析模型，我们运用这个模型做了一个实例计算。 实例问题为一个均质的粘土滑坡，具有球形滑动面，其它各种参数如图 5 所示。在图 5 中， c 为内聚力，φ为摩擦角，R 为瞬时摩擦力，γ为土的单位重量。运用封闭式（closed-） 算法得出三维安全系数为⒈402。运用 CLARA 模型算得安全系数为⒈422。同样的问题运 用三维边坡模型算得三维安全系数范围为⒈386 到⒈472，它取决于用于被分离的边坡柱体 的数量。 7 运用基于 GIS 栅格的三维稳定分析模型（图 5） ，并将格网尺寸定为 0.5m 时，算得三 维安全系数为 1.386；而当格网尺寸为 0.6m 时，算得安全系数为 1.388。很明显，与封闭式 算法相比，基于栅格模型的 GIS 可有效的用于三维边坡稳定性评估。 3 确定临界滑动表面和蒙特卡洛模拟确定临界滑动表面和蒙特卡洛模拟 滑动面只能通过岩土工程调查来确定，由于地质调查的费用比较昂贵，因此滑动面通 常是很难确定的。因此，边坡稳定性评价对临界滑动面的确定是非常重要的。 为了判定三维临界滑动情况， 利用蒙特卡洛随机模拟方法来计算三维安全系数最小值。 假定最初的滑动面是一个椭球体的较低部分， 边坡表面则根据不同地层受力情况和不连续界 面条件而改变。最终得到危险滑动的表面，同时可得到相关三维安全系数的最小值。 4 椭圆坐标转换椭圆坐标转换 假定最初的滑动面是一椭球体的较低部分， 椭球体的倾斜方向设置为与研究区主要的倾 斜方向一致；将椭圆的倾角基本上设定得近乎于研究区起伏变化的倾角。其主要倾向为α， 主倾角为β，它们是由边坡破坏区域主要栅格像元的值确定的。假定倾向和倾角属正常分 布，则将主要的倾向α，倾角代入分布模型中。 ββααcos.sin.cos.sin 2211 cscs 9 运用公式（10）和（11）完成坐标转换。图 6 显示了坐标转换过程。 8           − − −          − − −           − − −                     0 0 0 2 2 21 1 21 21 1 21 0 0 0 321 .0.. zz yy xx c s ss c cs sc s cc z y x z y x BBB z y x 10                               −−           0 0 0 2 21 21 1 1 2 21 21 0z y x z y x c ss sc c s s cs cc z y x 11 式中，x,y,z 为全球大地坐标， ,,zyx为当地坐标， 000 ,,zyx为椭球体中心点坐标。 5 Z 值的确定和滑动面的倾斜度值的确定和滑动面的倾斜度 滑动面上“B”点的 Z 值是通过中直线 AB 和椭圆的由公式（12）的结果确定的。 （见 图 7） 9 1 cossin 2 2 2 2 2 2 0 00 − − c z b y a x yy Slope zz Slope xx （12） 对于每个栅格像元，滑动面的倾向和倾角可通过下面的公式计算得出，像元（j，i）的 斜度可以通过图 8 中点 14 的 Z 值来确定。点 14 的值由公式（13） ， （14） ， （15）算出， 滑动面的倾向和倾角由公式（16）算出。 [] [] [] []4/1.1. 1. 1. 4/. 11. 11.. 4/. 11. 11.. 4/ 1.1. 1. 1. 4 3 2 1 −− −−− −−− i jZijZijZi jZZ ijZijZi jZi jZZ ijZijZi jZi jZZ i jZijZijZi jZZ （13） 4/ 3 4/ 3 4/ 3 4/ 3 43214 43213 43212 43211 ZZZZZ ZZZZZ ZZZZZ ZZZZZ − − − − （14） 10 13 12 ZZb ZZa z z − − （15） z z zz b a d ba − 0 22 tan tan β θ （16） 这里，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） , , , maxmin maxmin maxmin ccc bbb aaa ∈ ∈ ∈ （17） 假定 a，b，c 都均匀分布，则蒙特卡洛模拟的随机变量由公式（18）和（19）来算出。 在[0，1]范围内平均分布的随机变量可通过全等乘积方法得出 myr mModayy ii ii / 1 − （18） 式中，γi为在[0，1]范围内平均分布的随机变量。在[a，b]范围内平均分布的随机变量 可由公式（19）计算得出。 aabrx ii − （19） 式中，χi为在[a，b]范围内平均分布的随机变量。 椭球体的倾角设定为平均分布的一个随机变量。平均分布范围为主倾角及其在一个确 定的波动范围之内变化的变量。 7 计算三维安全系数最小值的过程计算三维安全系数最小值的过程 整个研究区（或边坡破坏范围）可以被均分为若干小矩形栅网，如同基于栅格的 GIS 11 一样。 关于此基于栅格的三维边坡稳定性分析的数值计算， 前面提到的所有的计算过程都可 以通过 Visual Basic （利用 GIS 组件）来完成。这个软件叫三维边坡地理信息系统，是运用 Visual Basic 6.0 和 ESRI 公司生产的 MapObjects 2.1 开发的。MapObjects 作为 GIS 的一个组件，用来对 GIS 数 据进行组织和空间分析。计算三维安全系数的过程如图 9 所示。 在这个过程中，数据模块的功能用来获得所有与边坡相关的地质、地貌、水动力学数据 和地质力学参数； 随机变量参数模块用来随机选择蒙特卡洛模拟的实验滑动面； 三维边坡稳 定性模块可用于计算三维安全系数；而危险滑动面及其安全系数可以通过一些实验计算得 出。在图 9 中可以看到关于 GIS 空间分析功能的所有模块可以通过 GIS 组件来实现。因为 一个 GIS 组件是在三维边坡地理信息系统系统中完成的，所以可以有效的计算三维安全系 数；同时利用与边坡相关的 GIS 数据，所有的相关数据和结果可以在三维边坡地理信息系 统系统中实现可视化。 实例剖面如图 10 所示。 在这个实例中考虑的因素有四种地层、 地下水和破坏面 （断层） ； 物理和力学参数如表 2 所示。 12 表 2 研究实例的物理和地质力学参数 地层 C（kN/m2） 0φ / 3 mkNγ 1 2 3 4 破坏面（断层） 9.8 58.8 49.8 64 9.8 30 25 30 35 20 22 24 26 27 -- 为确定临界滑动面，对蒙特卡洛随机计算次数进行了实验，总共计算次数达到了 1,000 次。每次实验计算的三维安全系数最小值的结果如图 11 所示。图中明显的显示了在实验中 计算 300 次后，可得到的安全系数最小值。这 300 次实验的结果见图 12，这些计算结果差 别不太大，其最小值为 1.34，最大值是 1.68。这个临界滑动的研究程序是建立在最小安全系 数的计算基础之上的。 而最小安全系数的计算取决于随机选择的参量。 有关这一临界滑动实 例的三维可视图见图 13。通过三维模型与二维模型结果的比较，用 Janbu 法确定临界滑动 面时，使用的是图 10 所示的二维模型和表 2 所列的参数，通过这种二维模型计算出的安全 系数为 1.18，这要比用三维模型计算出结果的极小值（1.346）略小一点。 13 8 滑坡滑动过程模拟滑坡滑动过程模拟 基于 GIS 栅格三维边坡稳定性分析模型和 GIS 栅格数据，对滑坡滑动过程进行了模拟， 直到三维安全系数大于 1 为止。滑动方向按滑动面的主滑方向确定。图 14 中展示了由滑面 确定的八个滑动方向。例如，若滑面方向的倾角在 22.5゜～67.5゜之间，则滑坡将要滑动的 方向恰在该图的右上方（即“5”方向） 。 滑坡滑动过程的模拟流程见图 15。首先，要计算滑坡初始状态时的三维安全系数，以 确定其滑动的可能性。若其安全系数小于 1，则接着进行下一步滑动过程模拟。先沿着由滑 面主倾向确定的滑动方向移动滑坡多边形；接着，在新的滑坡多边形范围内，分步（每一步 14 等于一个栅格大小）计算每一个栅格的 DEM 和滑动的变化，并再次计算下一步滑动的新滑 动方向。并在新的 DEM 数据和滑动多边形范围的基础上，计算出新的三维安全系数。如果 三维安全系数仍然小于 1，则进行以下的新滑动步骤模拟。 在这种滑动模拟模型中，假定滑动面内摩擦角不改变，但除了在初始三维边坡安全系 数的计算过程之外，假定滑动面没有内聚力（即内聚力为零） 。 仍然用同样的实例（如图 5 所示） ，用不同的两种动力学参数进行滑坡滑动过程模拟 情况 1 302 /23.11./4mKNmKNcγφ 情况 2 302 /23.5 .10./6mKNmKNcγφ 第一种情况下， 初始边坡安全系数为 0.82， 在进行７步滑动之后， 滑坡体开始趋于平稳， 其安全系数是 1.04。部分滑动步骤剖面及三维视图变化如图 16 所示。在此图中，DEM 的改 变及滑坡体移动过程一目了然。运用三维边坡地理信息系统，也可将可视滑动过程表现为 GIS 地图和剖面图的形式。滑坡体沿水平方向的最终滑动距离为 3.0m。 15 16 第二种情况下，滑坡体将一直向下滑动到平坦地区，水平方向滑动距离为 14m。滑坡体 最后停止滑动位置的三维展视图如图 17 所示。 ９ 讨论和结论９ 讨论和结论 在三维边坡稳定性柱状分析模型的基础上，开发了一个全新的基于三维确定性模型的 GIS 栅格，并且通过一个问题实例证实了其正确性。在三维边坡稳定性分析模型中，假定其 初始滑面为一椭圆面； 其三维临界滑面， 是利用蒙特卡洛随机模拟求取最小三维安全系数而 确定的。基于 GIS 的栅格三维模型，滑坡滑动过程模拟用于判断滑坡灾害和预测滑动距离。 已开发了作为计算程序软件的三维边坡地理信息系统， 它足以完成一切有关三维边坡问题的 计算，其中的 GIS 组件用于实现 GIS 的空间分析功能和有效数据的管理。因其具有空间分 析、数据管理和与边坡相关的综合数据的 GIS 可视化等优点，所以三维边坡稳定性问题已 经比较易于研究。自打全新的基于 GIS 栅格三维边坡稳定性分析模型问世，就为惯于使用 传统数学方法研究边坡稳定性的工作者拓展了一个新的研究领域和数据库方法。