引言

滑坡的稳定性分析与评价是工程地质中的一项经典工作，传统的分析方法一直以二维极限平衡法为主。但目前随着计算机软、硬件技术的飞速发展，基于真实“三维滑裂面”的坡体稳定性分析已成为共识。在已知或设定三维滑裂面的情况下，如何能够同时兼顾严密性和实用性，高效地确定坡体的稳定性是一个新的问题。国外以Hungr为代表的一些学者 (Hungr，1987；1994；Hungr等，1989；) 曾将经典的二维Janbu法和Bishop法扩展至三维，在满足三维极限平衡条件下，给出了一套较具代表性的实用化方法。在国内，许多学者 (冯树仁等，1999；陈祖煜等，2001；李同录等，2003；张均锋等，2005；陈胜宏等，2005；张常亮等，2010) 亦从不同的方面就三维滑裂面稳定性系数的计算方法进行了实用化的改进或扩展。例如，冯树仁等 (1999)通过忽略条块间垂向剪力，给出了一种类似于三维Janbu法的计算方法，可利用垂直方向和滑动方向力的平衡条件求得稳定性系数；陈祖煜等 (2001)将经典的二维Morgenstern-Spencer法扩展至三维；李同录等 (2003)提出了三维简化的Sarma法等等。

本文介绍作者所开发的一套实用型黄土坡体“三维最危险滑裂面”搜索和稳定性评价软件系统的核心思想，即采用Monte Carlo随机搜索法与遗传算法相结合的优化方法，高效生成一系列接近 (或包含)“最危险滑裂面”的三维滑裂面，并以Hungr法所确定的稳定性系数最小为筛选依据，搜索确定任意形态黄土坡体的“三维最危险滑裂面”，并基于已知的“三维最危险滑裂面”，进一步考虑各种可能的参数变化，如地震烈度、土参数分层差异、坡体雨雪载荷、坡角开挖等因素，进行更加严密的稳定性分析和评价。

1.   黄土滑坡滑裂面的几何特征

黄土滑坡的滑裂面在几何特征及发生机制上与岩质滑坡存在着本质的差异。岩质滑坡的滑裂面通常发生在先存的软弱结构面，如层理、节理、裂隙、断层面、不整合接触面等，因此，可以说岩质坡体的滑裂面是岩体本身固有的，是否出现滑动现象仅取决于滑裂面上的滑动条件。黄土滑坡情况却有所不同，尽管其滑裂面的产生也受一些先存的竖向结构性节理的影响，但在总体上，其实际滑裂面不是坡体内部本身先存的，而是在许多因素的综合作用下，坡体具备了滑动条件临时产生的。就物理学的角度而言，黄土边坡的失稳，其本质是坡体的一部分沿某一滑裂面所受“下滑力”与“抗滑力”失衡，导致该部分的宏观位移。在一定条件下，最有可能使坡体失稳的所谓“最危险滑裂面”，应是这样一个几何面：沿其发生的每个微量滑移所耗总能量与整个坡体的相应的势能降之比最小。而在实际工程中，“最危险滑裂面”可以理解为“总阻滑力”与“总促滑力”之比 (即稳定性系数F_s) 最小的滑裂面。因此，在理论上，我们可以以稳定性系数最小为标准，搜索出一定条件下的“最危险滑裂面”。

关于黄土坡体最危险滑裂面的几何形态，传统上往往将其过度地简化为二维上的圆弧面 (或椭圆弧面) 以及三维上的马蹄形面。实际上，无论是理论推测还是野外的实测勘察均表明，黄土滑坡滑裂面的大小和形状受坡体形态、黄土介质的均匀性状况、黄土竖向结构性节理的发育状况、黄土坡体基岩基底状况等等一系列因素的影响和控制，往往并非简单的马蹄形，更为普遍的情形是整体大致呈具有陡立后缘和左右臂、滑坡床具有一定宽度的马蹄形、局部复杂而非规则的光滑曲面。因此，我们采用计算机模拟生成一系列可能黄土滑裂面时，无需施加规则或简单曲面的假设约束。

2.   黄土边坡“三维最危险滑裂面”的搜索确定

对于任意黄土坡体，除确定其坡体表面形态外，还需探明其坡体基岩基底的具体状况。在此基础上，我们可首先在实际滑裂面可能出现的大致范围内，按照Monte Carlo法随机地模拟生成一系列“可能的”三维滑裂面，这些滑裂面虽然形态、位置和尺寸随机各异，但每一条在几何形态上应基本合乎上述黄土滑坡滑裂面的主要特征，即整体大致呈马蹄形、局部非规则的光滑曲面。

对于每一个三维滑裂面，我们采用国际上使用较广的Hungr法计算其对应的稳定性系数。理论上，只要采用Monte Carlo随机搜索法生成的“三维滑裂面”的数量足够大，则通过稳定性系数“小则优”的准则，总可筛选出“最危险滑裂面”(图 1)。但是，由于采用这种单纯的随机搜索，每一条滑裂面彼此独立生成，“最危险滑裂面”是随机“碰”到的，故效率非常低。为此，我们在Monte Carlo随机搜索的基础上引入了遗传算法的基本思想，按以下思路实现高效搜索：

图 1  Monte Carlo随机搜索编程思路

Figure 1.  Flowchart for Monte Carlo random search programming method

下载: 全尺寸图片 幻灯片

(1) 采用三维滑裂面Monte Carlo随机模拟系统，在具体黄土坡体的合理范围内随机搜索出N个相对最优的“候选滑裂面”，作为遗传算法的第一代滑裂面“种群”。

(2) 将每个三维滑裂面用结构统一的“三维网格点阵”(三维矩阵) 进行归一化的数学表述。这里所谓“归一化”是指同一坡体不同大小和形态的三维滑裂面均具有完全相同的平面二维网格节点 (X_i，Y_i)，彼此差异只体现为第三维数值 (Z_i) 各不相同，由此实现了用统一的数据结构表述不同滑裂面的空间形态 (图 2)，其中，Z为斜坡的高度，Y为长度，X为宽度。

图 2  黄土坡体任意三维滑裂面的“归一化”数学表述示意

Figure 2.  A normalized mathematical representation of any 3D slip surface on a loess slope

下载: 全尺寸图片 幻灯片

(3) 将三维网格点阵看作“滑裂面染色体”。由于各染色体结构一致，故可借鉴生物遗传进化的过程，让滑裂面染色体随机变异和相互杂交，不断产生新的进化改良型滑裂面。具体而言，随机变异即相当于让单个滑裂面矩阵的第三维数值Z_i做局部随机变化，进行自我突变改良。相互杂交即让不同滑裂面某些对应部位进行彼此互换或加权平均等，进行取长补短的改良。

(4) 计算每一滑裂面的稳定系数，以稳定系数“小则优”为原则，优胜劣汰，进化出N个新一代种群。

(5) 不断重复上述遗传过程，“滑裂面种群”将最终进化 (收敛) 为一组稳定系数最小的“最危险滑裂面”。

(6) 为了防止“滑裂面种群”收敛于“局部最优”，将采用独立搜索、多次比较结果以及加入扰动因子迫使滑裂面种群跳出局部极值的措施 (图 3)。

图 3  Monte Carlo随机搜索法与遗传算法相结合的搜索编程思路

Figure 3.  Flowchart of programming method for the combination of Monte Carlo random search and genetic algorithm

下载: 全尺寸图片 幻灯片

3.   “三维滑裂面”稳定系数计算方法

本系统在滑裂面搜索过程中采用目前国际上广泛使用的Hungr法计算三维滑裂面稳定性系数。该方法继承了Bishop简化法 (Bishop，1955) 的迭代算法，其特点在于考虑了“土柱”上的正压力和剪切力，详细计算原理见参考文献 (Hungr，1987；1994；1997；Hungr等，1989)，此处不再赘述。下面简单介绍其关键性的计算原理。

假设三维滑裂面以上的滑动体被规则地切割为一系列垂直土柱，则其受力状况如图 4所示。

图 4  单个土柱的受力情况示意图

Figure 4.  Schematic diagram of the stress situation of a single soil column

下载: 全尺寸图片 幻灯片

设W为土柱的重量，a为地震的水平加速度，X_y为垂直剪应力，E_y为正应力，整个稳定性系数F的计算过程如下：

(1) 假设λ值。

(2) 对每一行土柱，用公式 (1) 计算a_c。a_c对于土柱的任意给定行是个常量。

(3) 用公式 (2) 计算正应力E_y：

(4) 用公式 (3) 计算剪应力X_y：

垂直剪切力X_y和正压力E_y之间的关系由土柱间的力函数给出。此函数由Morgenstern Price给定：

λ是常量。土柱间的力函数对于Spencer法 (Spencer，1967) 是个常数或者对于Morgenstern-Price法 (Morgenstern等，1965) 是一个从滑坡上缘头到坡底距离的半正弦函数。

(5) 按照Bishop的简化法迭代计算，将剪切力X_y加在土柱的重力之上。

(6) 将整个滑动体看作是一个整体，为了满足水平力的平衡，不断地迭代变化λ值，直到满足公式 (4)。从而可确定出稳定系数F。

这里$\sum {{F_n}} $是所有土柱水平力F_n的总和。$\sum {{a_c}W} $是所有土柱的a_c系数和自重W乘积的总和。

4.   实例应用

本文选取甘肃省兰州市皋兰县城东南约10km处一个黄土坡体作为实例，该坡体地处黄土梁峁与基岩交界处，属自然斜坡 (图 5)。坡体中轴线坡面长约70m，宽约40m，平均坡度约25°，最陡处约50°。坡体上覆黄土层为第四系 (Q₃) 马兰黄土，下伏花岗岩风化壳顶板，坡底黄土层厚度≥12m，坡顶黄土层厚度≥30m，为典型黄土边坡 (图 6)，H、L分别为斜坡的坡体高度和水平长度。基于坡体样品的有关黄土力学参数见表 1。

图 5  黄土高边坡实例现场

Figure 5.  The field site of a high loess slope used for the experimental study

下载: 全尺寸图片 幻灯片

图 6  黄土高边坡实测剖面

Figure 6.  Profile of the high loess slope

下载: 全尺寸图片 幻灯片

表 1  皋兰县黄土高边坡土力学参数及稳定性预测结果

Table 1.  The soil mechanical parameters of the high loess slope in Gaolan county and the predicted result for its stability

取样环境	土样号	容重ρ/g·cm-3	含水量ω/%	振次N	粘聚力C/kPa	内摩角Φ/°	Fs(不同地震系数)
							0.1	0.15	0.2
兰州皋兰破坏应变3%	L-1	1.45	7.9	10	36.6	22.4	1.66	1.49
				20	35.2	21.8			1.3
	L-2	1.44	10.63	10	32.6	21.6	1.52	1.36
				20	30.0	21.4			1.16
	L-3	1.62	17.86	10	22.4	19.2	1.05	0.94
				20	21.6	18.7			0.82
	L-4	1.72	23.40	10	21.0	17.0	0.92	0.82
				20	19.0	16.5			0.69

下载: 导出CSV 

| 显示表格

图 7是采用黄土坡体“三维最危险滑裂面”搜索系统所确定的该坡体的三维潜在最危险滑裂面。表 1列出了该滑裂面所对应的不同土力学参数和地震加速度参数下的稳定性系数。显然，通过尝试性地变换各参数的量值，可以从数值模拟的角度判定每种参数变化对坡体稳定系数的影响程度和敏感性，从而为进行坡体稳定性分析和坡体安全性防护措施的制定提供一定的参考依据。

图 7  “三维最危险滑裂面”搜索系统所确定的马蹄型三维潜在最危险滑裂面

Figure 7.  The most possible 3D "horseshoe-like" slip surface predicted by our integrated search system

下载: 全尺寸图片 幻灯片

5.   结论

对于任意给定的黄土坡体，当我们采用结构统一的“三维网格点阵”(三维矩阵) 对其各种可能的三维滑裂面进行归一化的数学表述后，即可方便地将Monte Carlo随机搜索法与遗传算法相结合，实现“三维最危险滑裂面”的高效搜索。其中，对于每个模拟生成的三维滑裂面的稳定性系数计算，可采用兼具严密性和实用性特点的Hungr法。基于搜索确定的三维最危险滑裂面，能够更好地分析和评估在各种可能的地震烈度、土参数变化条件及各种防护条件下黄土坡体的稳定性状况，对我国广大西北地区黄土滑坡灾害的预测和防治具有较好的参考价值。