1.2 国内外的相关研究及进展

斜坡或边坡作为一种人类不可回避的地学环境与工程形式，总与人类的工程活动伴生，为此人类为了安全起见始终关注着边坡的稳定性。边坡在受到人类的工程活动及外界环境影响时，坡体有可能发生破坏，给人类的带来灾害。100多年来，人们对边坡变形过程、失稳形式、失稳机制、稳定评价及滑坡预测预报等进行了广泛而深入的研究，借助数学、力学及计算科学的理论和方法，试图对边坡的稳定、演化及滑坡的预测预报进行研究，并应用到工程实践中去。经过国内外无数工程地质和岩石力学、岩土工程工作者的努力，已形成了针对边坡工程的理论体系和工作方法，为人类工程建设活动奠定了一定的理论及实践基础。然而，由于岩土体赋存于一定的地质环境中，在其形成、变形和破坏的过程中，经历了一定的演化历史和各种地质作用，从而形成了自身的几何及物理力学特性，具有一定的物质组成，一定的结构和构造，并继续在特定的地质环境中遭受质的变形、破坏［4］，从而使岩土体具有极高的复杂性，加上人类认知能力有限，因而此领域的研究难度较大，还有大量的难题未解决，下面分静力问题和动力问题两个方面对边坡稳定性方面的研究进展做一综述。

1.2.1 边坡静力稳定问题的研究及进展

在边坡的静力稳定问题研究中，岩土体的物理力学参数的准确性会对边坡稳定的评价结果产生重大的影响，因此本节首先对边坡物理力学参数选取的研究进行总结，然后讨论边坡静力稳定问题的研究及进展。

1.2.1.1 边坡静力学参数选取的研究

人类所建设的工程无不希望建筑物稳定可靠且工程造价低，要做到这一点其关键性的一步是建筑物设计力学参数选取合理，参数值的高低会影响建筑物的建造尺寸，当然也就涉及工程投资。如果所选参数可靠度较低，往往会导致建筑物结构破坏，甚至有灾难性事故发生。根据世界各国的调查和统计，45%的重力坝事故是由于地质问题没有查清，岩体力学参数选取不合理造成的［5］，如法国马尔帕赛拱坝高66m，1959年溃坝，造成近500人死亡，究其原因是地质条件没有查清，岩体力学参数选取过高所致。由此可见，岩体力学参数选取是边坡稳定分析的重要内容，只有准确的力学参数才能得到符合工程实际的计算结果，岩体力学参数在工程建设中所起的作用极为重要，会直接影响到工程的造价和安全。

当前国内外岩体力学参数选取研究的总趋势是从有经验、半经验、精度较低的数值计算方法向考虑多种因素影响、计算过程复杂、精度较高代表性较强的数值计算分析法发展，尤其是计算机的使用，使这一领域的研究加快。

岩体力学参数选取常用的方法有点群中心法、优定斜率法、最小二乘法、随机-模糊法等。点群中心法由于人为因素影响过多，目前已不常采用，国内对于岩体力学参数的研究主要是从岩体力学参数本身所包含的随机性和模糊性出发，应用随机理论和模糊数学的方法，对试验所得的数据进行分析以获得更为逼近岩体力学实际参数的“真值”。李华晔等人建立了最小二乘选取c、φ值的算法并分析了300多组岩体抗剪强度试验成果［6］，成都勘测设计院改进了日本人的优定斜率法，并用于二滩、溪洛渡等工程岩体抗剪参数选取值［7］，熊文林、李华晔、黄志全等人建立了优选c、φ值的随机-模糊分析方法［8-10］，李华晔、黄志全等人把随机-模糊分析用于小浪底、溪洛渡、宝泉抽水蓄能电站等工程岩体c、φ值计算^［10^，11^］，胡小荣、谭文辉、张征等人用泛克拉格法、中心点离散法、局部平均离散法和分形理论对岩体单轴抗压、抗拉、变形模量等参数进行了分析和计算［12-15］，谭文辉将修正的地质强度指标（GSI）和非线性广义Hoek-Brown破坏准则相结合，在岩石三轴试验和现场地质调查的基础上，对岩体的宏观力学参数进行了研究［16］。

国外对于岩体力学参数选取的研究与国内大体相似，但他们在研究时间上更早些，如Brown E.T.等人对现场应力渗透系数和干密度的研究［17］，Derkiareghian和Zhu W.Q.等人用中心点离散法对单轴抗压、单轴抗拉所做的随机有限元分析等［5］。

1.2.1.2 边坡静力稳定分析方法研究

边坡稳定分析是进行边坡研究的核心内容，人们根据边坡的几何形态、破坏方式等发展了各种适用于不同滑坡类型的评价方法。

1.工程地质分析法

工程地质分析法的核心内容是对边坡进行工程地质定性评价，最早应用工程地质类比法，随着工作的深入，人们逐渐进入理性分析阶段，工程地质分析法的理论基础是地质成因演化论及20世纪60年代末、70年代初被明确提出的岩体结构的概念［18］，在这之后，大量基于该理论的研究不断进行，20世纪80年代是岩体结构理论发展的高潮期，相关专著不断出版^［17^，19-23^］。基于岩体结构理论的工程地质分析方法在边坡稳定评价中占有重要的地位，特别是对于地质条件复杂的岩质高边坡，工程地质分析法更有其独特的价值［2］。

岩体结构力学［23］的发展促进了工程地质分析法的发展，多数边坡的破坏是由岩体结构面控制的，在岩体结构理论的指导下，人们开始研究边坡的破坏模式，为定量研究边坡稳定性奠定基础。1980年李铁汉等人以滑动面的形态、数目、组合特征以及边坡岩体破坏的力学性质，将边坡变形破坏划分为5类，每类中又分为若干亚类［24］。孙玉科等人（1983，1988）总结了我国岩质边坡变形破坏的主要地质模式，提出了边坡变形的常见5大模式，即金川模式（反倾边坡）、盐池河模式（水平层状上硬下软）、葛洲坝模式（水平薄层状软硬相间）、白灰厂模式（水平厚层状软硬相间）以及塘岩光模式（顺倾薄层状结构）^［20^，24^］，这些基于岩体结构理论的成果的提出，为准确确定边坡工程地质模式，从而进行准确定量计算作出了巨大的贡献。

2.刚体极限平衡分析法

刚体极限平衡法是边坡稳定分析的重要方法之一，也是目前工程设计部门在进行简单的边坡稳定分析中最常用的方法，该方法是以莫尔-库仑抗剪强度理论为基础，建立滑坡体力或力矩平衡方程，通过一定的假定条件，减少未知量的个数，从而将边坡稳定的超静定问题转化为静定问题，然后求解方程组，得到边坡的安全系数［25］。这种方法没有像许多数值分析法一样，引入弹性力学或结构力学求解超静定问题的解法；所以理论简单，操作简洁，可用于求解介质材料均一，边界条件简单的边坡稳定问题。

条分法是刚体极限平衡法的重要方法之一，1916年由瑞典人彼德森提出，后来经过Fellenious、Taylor等人不断改进，现已成为简单边坡稳定分析的经典方法，他们假定边坡稳定问题是一个平面应力问题，滑裂面是个圆柱面，计算中不考虑条块之间的作用力，边坡稳定的安全系数是用滑裂面上全部抗滑力矩与滑动力矩之比来定义［26］。20世纪40年代以后，随着岩土力学的发展，很多学者致力于对条分法的改进，他们的努力大致有两个方面：一方面是着重探索最危险滑弧位置的规律，制作数表、曲线，以减少计算工作量；另一方面是对基本假定做修改和补充，提出新的计算方法，使计算结果更加符合工程实际。1955年Bishop对传统的瑞典圆弧法进行了重要的改进，提出了关于安全系数定义的改变，对条分法的发展起了重要的作用［27］。在这之后的几十年里，世界各国的学者纷纷致力于通过力的平衡条件来确定安全系数，如Janbu（1954）假定土条间力为水平［28］，美国陆军工程师团假定条间力的倾角等于平均坝坡的坡角［29］，这些方法实际上都是通过假定条件减少未知量的个数，在满足合理性要求的前提下求解安全系数，所以严格地讲，他们都是一些简化的方法。

1965年，随着计算机的出现和普及，早期的手工计算已经被计算机所代替，人们在计算机的帮助下，开始努力寻找新的、更为严格、精确和实用的求导边坡安全系数的方法，Morgenstern和Price提出了适用于任意形状滑裂面的严格分析方法，Spencer（1967）提出了土条侧向力相互平行的假定，建立了同时满足力和力矩平衡方程的分析方法。

以上的计算方法从严格的角度上来说，对于均质土坡分析较为理想，而对于一些介质和边界条件复杂的岩质边坡是不适用的，Sarma（1979）提出了一种重要的刚体极限平衡方法，解决了很多长期困扰科学家的问题，他认为边坡岩体只有沿着一个理想的平面或圆弧滑动，才能成为一个理想的刚体运动，否则岩体必先破坏成为多块可相互滑动的块体之后才可能发生滑动。由于这种方法对于岩质边坡分析能力的增强，后来被广泛地应用于边坡的分析中，并被称为Sarma法^［30^，31^］。

1983年，陈祖煜对Morgenstern-Price法作了重要改进［32］，完整地导出了静力平衡微分方程的解，提出求解安全系数的解析方法，从根本上解决了数值分析的收敛问题，进一步推动了刚体极限平衡法的发展。

目前，刚体极限平衡方法已经从二维发展到了三维，关于边坡稳定分析的三维极限平衡方法，已有很多文献介绍，Duncan曾经总结了24篇文献资料，列举了这些方法的特点和局限性［33］。

3.数值计算方法

由于刚体极限平衡法自身固有的局限性，在对边坡进行稳定分析时，需对滑坡边界条件大大地进行简化，计算中所选用的各种计算参数往往是确定的且呈线性变化，对于由复杂介质和边界组成的滑坡体，如果进行这样的简化处理，将不能客观地反映工程实际的真实性，且使计算结果有很大的误差。实际上，不仅滑坡的各种计算参数是不确定且随机的，而且边坡系统本身就是一个不平衡、不稳定、充满复杂性的动态系统，其与外界环境不断地有物质、能量和信息交换。鉴于此，科学家们为了寻求能够反应工程实际情况的方法，进行了各种尝试，在数值计算方法上，随着计算机的普及和发展，出现了一批以弹性力学、结构力学为基础的数值计算方法：FDM（有限差分法）、FEM（有限单元法）、DEM（离散单元法）、DDA（不连续变形分析）、FLAC（快速拉格朗日插值）、NNM（流形元方法）等方法都是数值计算方法飞速发展的产物。

在处理复杂边界和介质的计算上，数值分析方法由于其强大的灵活性和通用性，在工程技术的各个领域发挥了极大的作用，现已成为广大工程技术人员的得力计算工具。

4.边坡稳定分析的系统方法

边坡稳定分析的系统方法是近年来兴起的一类评价方法，其主要理论依据是边坡变形破坏行为所表现出的非线性，这些非线性主要表现在：

（1）岩石材料的无序分布，地应力的随时空分布。

（2）边坡系统是高度复杂、规模宏大的系统，影响边坡稳定的因素很多，而通常条件下，难以确定边坡的地质条件和环境信息。岩体的变形、损伤、破坏及演化过程包含了相互耦合的多种非线性过程，传统的确定性分析理论和力学方法无法描述如此复杂的力学行为。

（3）岩体变形经过塑性、断裂和破坏过程中，会在系统中出现分叉、突变等非线性复杂力学行为。

（4）岩石的变形、损伤和破坏过程是一个动态的非线性不可逆过程。

复杂非线性系统科学使我们对事物的研究实现了由静态到动态、由局部到整体、由线性到非线性、由简单到复杂的认知走向［34］。岩体工程问题的极端复杂性是其他建筑工程问题所不能比拟的，因而在经典的理论计算方面遇到了严重困难，而20世纪70年代出现的复杂非线性系统科学为解决这个问题开辟新的研究领域，谢和平等人（1996）认为：“岩土工程失稳的研究要取得突破性发展，迫切需要引进非线性科学研究的原理和方法［35］”。郑颖等人（1996）认为，把岩体的破坏与远离平衡条件下的非线性动力系统理论联系起来，可能成为21世纪岩石或岩体理论的突破口［36］。基于上述思想，为了更准确地描述边坡非线性复杂动力学过程，科学家们将混沌理论中的耗散结构、协同分岔、灰色系统、突变理论、分形理论以及神经网络系统等非线性科学引入边坡稳定分析中，秦四清等人结合大量的位移观测资料，利用滑坡系统动力行为中的多组状态变量，采取非线性动力学方法反演边坡系统的非线性动力学方程，对滑坡系统的动态行为进行了有益的尝试［37］；周萃英等人从耗散结构理论和协同论的观点出发，研究边坡系统破坏的自组织行为取得一定的成果［38］；李造鼎等人以能量突变为材料的破坏特征提出了灰色能量突变模型预测边坡稳定分析的方法［39］。

非线性理论近年来取得了突飞猛进的发展，目前已形成一整套边坡稳定分析的非线性解决方案，但这些研究目前仍然处于理论研究阶段，具体的工程适用效果仍有待检验。

5.可靠性理论与边坡可靠性稳定分析的研究进展

可靠度分析方法是一门正在迅速发展的新学科，被誉为近代工程技术的重要发展，无论从国民经济建设，还是社会效益方面，可靠度研究都起到了非常重大的作用［40］。自20世纪70年代初应用于边坡工程领域以来，由于国内外岩土力学界的重视，可靠度分析方法已经得到长足发展，并成为边坡工程研究的重要发展方向［41］。

由于岩土介质的特殊性，在工程设计、施工、适用过程中具有种种影响工程安全、适用、耐久的不确定性，这些不确定性大致有以下几个方面。

（1）岩土体参数的离散性、随机性。

由于自然界中的岩土体是非均匀、非连续、各向异性的地质体，岩土体力学参数受矿物成分、结构构造、水等因素的控制，因此在实际中岩土体力学参数难以通过测量获得准确值。所以整个工程系统就具有了随机性、不确定性，它是在一定范围内变化的随机变量。因此，在工程设计中采用确定性方法进行建筑物安全性能评估是不准确的。研究这类问题的方法有概率论、数理统计和随机过程。

（2）安全系数的模糊性。

在传统的工程设计方法中，安全系数是最为常用的评价建筑物或者工程结构安全的标准，它常被理解为材料的破坏强度对于使用中的最大应力之比，并简单地认为一个较高的安全系数就能保证结构的安全［40］。而在岩土工程中，把安全系数定为破坏应力与最大应力之比，只适用于保持弹性状态直到破坏的时为止的岩土体，因此通常所说的安全系数不能给出结构可靠度的直接信息，常规的安全系数的模糊性可能引导工程技术人员得到一个虚假的安全感，从而影响人们对工程结构安全的认识，导致工程事故的发生［42］。

将概率论和数理统计方法应用于结构可靠度分析的最早尝试可以追溯到20世纪20年代［43］，尽管早期的研究工作富有创造性，但囿于当时的科技发展水平和显示需求，基于可靠度的结构分析方法并未引起社会的足够重视。第二次世界大战期间及随后的岁月里，有关机电设备、船舶、压力容器、飞行装置和海上石油钻井平台等复杂结构，在设计使用寿命期限内，在规定的荷载条件与环境条件下不能预期正常工作的事例不断增多和日趋严重，说明以安全系数设计法为代表的传统设计方法对环境条件和结构特性的确定性假设是不适当的，必须以概率论的观点出发，对有关的设计参量进行统计分析，研究他们的分布规律和相关特性，从而制订出一整套新的、符合实际情况的结构设计规范[40]。

1947年苏联的尔然尼钦提出了一次二阶矩的基本概念，标志着结构安全度的研究进入了应用概率论和数理统计学的阶段［40］。

美国的Freudenthal在20世纪40年代开创了美国安全度的研究工作，它的研究工作与同年代苏联的尔然尼钦的工作，有某些相似之处，他在1951年提出，破坏概率的选择，应使结构建造费与期望的破坏损失费的总和最小。A.H-S.Ang发展了Freudenthal的工作，提出了广义可靠性概率法［40］。Cornell于20世纪60年代后期提出了与尔然尼钦相同的一次二阶矩法及安全指标β的概念，对安全度的实际应用作出了贡献［44］。此后，结构可靠度的研究取得了长足的发展，加拿大的Hasofer和Lind对Cornell的算法进行了修正，提出了H-L算法，随后Rackwitz和Fiessler将H-L算法的适用条件推广为任意随机变量构成的安全裕量方程，提出了R-F算法^［43^，44^］。由于R-F算法良好的普适性，目前已被国际结构安全性联合会（JCSS）所采纳，并正式命名为JC算法，此后关于结构可靠度的研究都是围绕着对JC算法的改进展开的，正是由于Hasofer、Lind、Rackwitz和Fiessler等人的卓越贡献，20世纪80年代前后，与有限元法、计算机应用技术和随机网络分析理论的迅猛发展潮流相呼应，以一次二阶矩方法为基础的结构可靠性分析理论和应用技术开始了由元件级水平向系统级水平的实质性过渡［43］。从而出现了一种称为分项安全系数设计法的细节可靠性设计方法，我国制订的有关港口工程、水利水电工程和铁路工程的结构可靠性设计规范^［45^，46^］都是以分项安全系数设计法为基础的。

可靠度理论应用于边坡工程是从20世纪70年代后期开始的^［41^，42］，许多学者对边坡可靠度的研究作出了卓有成效的贡献。Kraft、McMahon、Matsuo、Kuroda以及Harr曾对土质边坡作过典型研究论述了概率方法在边坡中的应用；NgMyen和Chowdhury报告了矿山排土场的可靠性分析；Cambou等人根据线性一次逼近理论，采用随机有限元方法进行可靠性分析^［25^，41^，47-51^］。在我国，边坡可靠性研究工作开展较晚，时至1983年，攀钢石灰石矿边坡可靠性分析与经济分析研究课题才作为第一个这个领域的研究成果通过冶金部部级鉴定［41］。在这之后，国内工程地质与岩土工程专家对可靠度理论在边坡工程中的应用给予了足够的重视，科研和学术成果纷纷发表，1990年，日本学者松尾稔关于岩土工程可靠性设计的理论和实践的专著《地基工程学》被翻译出版，该书首次系统阐述了土坡可靠性设计的理论和方法［52］，在这之后，高谦、王思敬（1991）利用可靠度理论对红水河龙滩水电站船闸高边坡的可靠性进行了研究，将可靠度理论系统地应用到岩质高边坡的优化设计和决策上［53］。1993年，我国学者祝玉学的《边坡可靠性分析》一书的出版，标志着我国对于边坡可靠性理论方面的研究进入了一个崭新的阶段［41］。此后，国内外众多学者均开展了这方面的研究，相继取得了一定的进展。

在可靠度分析中，有许多种计算方法，由于出发点或解决问题的角度不同，使各种方法存在一定的差异。可以预计，可靠性分析方法的推广和应用必将进一步推动我国边坡工程技术的进步［54］。

目前，边坡可靠性分析所常用的方法有：

（1）蒙特卡洛模拟法。

蒙特卡洛（Monte Carlo）法又称随机模拟法或统计试验法，是一种依据统计抽样理论，利用电子计算机研究随机变量的数值计算方法。在目前可靠度计算中，是相对精确的方法。它源于第二次世界大战期间，Von Neuman和Ulam在对裂变物质的中子随机扩散进行直接模拟中，以世界闻名的赌城蒙特卡洛作为该项研究的秘密称呼而得名。

由于影响边坡稳定性的不确定因素很多，而且由于边坡工程浩大，不适宜也不可能在实际作业过程中，以高昂代价对边坡形成、破坏过程进行直接的实验研究。面对这样理论上和实验上的双重困难，采用蒙特卡洛模拟方法，通过构造模型和进行模拟试验，把复杂的空间条件和时间过程简化成集总参数，进而求得安全系数或安全储备的概率分布及分布参数的特征值均值、标准差、破坏概率等，以此作为边坡工程问题的近似解，是极为适宜的［54］。其基本思想是，若已知状态变量的概率分布，根据边坡的极限状态条件g（X₁，X₂，…，X_n）=0，利用蒙特卡洛方法产生符合状态变量概率分布的一组随机数x₁，x₂，…，x_n，并代入状态函数g（X₁，X₂，…，X_n）=0计算得到状态函数的一个随机数，如此用同样的方法产生N个状态函数的随机数。如果在N个状态函数的随机数中有M个小于或等于1（以安全系数表示边坡状态）或小于或等于零（以安全储备表示边坡状态），当N足够大时，根据大数定律，此时的频率已近似于概率，因而可得边坡的破坏概率：

也可由已得的N个g（x）来求均值u_g和标准差σ_g，从而得到可靠指标β（β=u_g/σ_g）。

最近5年国内学术界对蒙特卡洛法应用于边坡工程的研究主要集中在改进其收敛速度上，杨林德（1999）提出通过重构响应面来提高蒙特卡洛模拟法的计算效率［55］；蔡为新（2000）在直接蒙特卡洛模拟方法中采用了Neumann展式，并采用谱分析方法，在模拟过程中使用FFT技术来减少计算时间［56］；张社荣（1999）利用Bayes方法对断层抗剪强度参数进行分析，较好地解决了小样本试验数据问题，从而减少蒙特卡洛模拟次数，提高了计算效率［54］；闫强刚、冯紫良（1999）根据NTM的均匀分布理论，将NT-net引入岩土工程可靠度分析的蒙特卡洛方法中，加速了蒙特卡洛方法的收敛速度［57］。

（2）一次二阶矩法。

蒙特卡洛模拟法需要预知各基本状态变量的分布形式和参数。当基本变量X_i（i=1，2，…，M）的概率密度未知，或者在概率密度函数复杂不易求其分布参数的积分时，可利用Taylor级数展开后忽略两次以上的项，只考虑它们的一阶原点矩（即均值）和二阶中心矩（即方差）这两个特征参数，近似地计算状态函数的均值和方差，求得可靠指标和破坏概率，故也称作一次二阶矩法^［40^，41^，44^，47^，48^，52^］。

对于任何一个给定的边坡，其状态函数为n个基本变量的函数，即Z=g（X）=g（X₁，X₂，…，X_n），而每种状态就是一组随机变量的实现，即Z_i=g（X₁，X₂，…，X_n）。换句话说，变量Z_i是n维基本变量空间Ω的一个点。因此，极限状态面可定义在n维基本变量空间Ω，极限状态面即用极限状态方程：Z=g（X）=g（X₁，X₂，…，X_n）=0刻画的曲面，这个面上所有点的基本状态变量的组合结果Z=g（X）=0，它把基本变量空间划分为破坏域Ω_f和安全域Ω_s。破坏域包含导致破坏［即g（X）≤0］的一组基本变量的实现，安全域包含导致安全［即g（X）＞0］的一组基本变量的实现。

这种方法的基本思路是把极限状态面变换到标准化正态分布空间U，如图1-1所示，假设极限状态面近似为一个超平面，该超平面为极限状态面距坐标原点最近的点的切平面，而从原点到该切点的距离则等于可靠指标β。这样，问题得以简化，就仅仅需要验算极限状态面上一个点。

对于非线性极限状态方程，如果对状态函数在平均值处用Taylor级数展开即用平均值赋值，称为中心点法。中心点法会带来较大的误差，因为由平均值所确定的坐标点位于安全区，而不是在极限状态曲面上。为此，而将Taylor展开式的点选在极限状态曲面上，并且只在破坏概率最大的这个点上，这个赋值点称为验算点，经过这样改进后的一次二阶矩法称为验算点法。

在边坡稳定性的可靠性分析方法中，最常用的方法是一次二阶矩法（First-Order Reliability Method，FORM）中的验算点法，经过20年的发展，它已成功应用于我国的多个工程的优化设计中，我国的尖山铁矿边坡设计［41］、二滩水电站龙山滑坡分析［58］、三峡库区滑坡［59］、清江水布垭坝址大岩淌滑坡［60］等大型工程项目均采用一次二阶矩法对边坡的稳定性进行了评价。

（3）统计矩法。

统计矩法是20世纪80年代初开始引入边坡可靠性分析的一种方法，它的基本数学工具是Rosenblueth E.于1975年提出的统计矩点估计方法，故又称为Rosenblueth法。这是一种近似方法，当单个状态变量的概率分布为未知时，只要利用它们的均值和方差，就可以求得状态函数（安全系数或安全储备）的一阶矩（均值）、二阶矩（方差）以及三、四阶矩，从而求得边坡的可靠指标，且在状态函数值的假定概率分布下求得破坏概率P_f，因而便于应用［61-63］。我国龙滩水电站进水口边坡倾倒稳定分析采用的就是统计矩法［64］。

（4）随机有限元法。

随机有限元法是可靠指标法与数值法相耦合的方法，包括一次线性逼近法和迭代验算法。自20世纪60年代以来，随着有限单元法在边坡工程分析、评价和设计中的深入发展，对岩体力学模型、影响应力及变形因素的定量化提出了更高的要求，期望为计算分析提供更合理的模型和更精确的数据。然而，由于岩体组成、结构、强度、荷载条件以及力学状态的不确定性和变化性，岩体的力学参数并非常数，而是一组随机变量，这样，用确定性分析方法不可能真实、确定地描述边坡的本征属性。因为再精确的计算数据，也只是相对的经验性确认，其分析结果自然会带有一定程度的局限性。因此，近几年来，国际上开始研究一个新的领域，即探索随机模型与数值解法的耦合，以寻求数值方法的概率解答，逐步发展成随机有限元法，并正在走向工程应用^［40^，41^，44^，52^］。

边坡随机有限元分析的实质，是在常规有限元方法基础上，改输入常数为随机变量{X}，如弹性模量（E）、泊松系数（μ）、容重（γ）等。因为单元应力{σ}和结点位移{δ}均为随机变量{X}的函数，只要随机变量存在均值和方差，则可以求得{σ}和{δ}的统计上的均值和方差。岩体强度分量黏结力（c）和摩擦角（φ）亦为随机变量。依照Mohr-Coulomb破坏准则，根据剪应力与抗剪强度的关系，可以建立起抗剪稳定性的极限状态方程Z=g（X）=g（X₁，X₂，…，X_n）=0。在由n个随机变量组成的n维状态空间中，Z=0构成一个极限状态曲面，而该曲面至变量空间坐标原点的最短距离即为单元的局部可靠指标。在标准正态空间，由式P_f=1-φ（β）可求得单元的局部破坏概率。显然，在有确定的潜在破坏面情况下，根据破坏面所穿过的所有单元的力学参数和几何参数，可以求得边坡总体的可靠指标和破坏概率。

随机有限元方法集有限单元法和可靠性分析的优点于一身，适应性强，考虑问题全面，精度较高，可解决非线性问题，随着计算机技术和有限元算法技术的提高，其工作量庞大和计算时间长的缺点基本已被克服，所以虽然发展历史不长，但随机有限元分析已经成为可靠度分析的一个重要方法之一^［63^，65^］。

1.2.2 边坡动力稳定问题的研究及进展

由于地震力等动力荷载的特殊性，地震荷载作用下边坡的稳定性分析方法与边坡的静力分析方法有很大的区别，其主要研究内容为：

（1）研究地震力的作用方式，即研究地震力如何作用于边坡，并讨论地震力的计算。

（2）研究地震荷载作用下边坡的稳定性，并判断边坡失稳的可能性。

（3）研究地震作用下边坡发生失稳的位置，讨论边坡破坏面的形状和位置。

（4）如果边坡在地震作用下失稳，还需研究边坡发生失稳的结果，即研究地震作用下边坡的永久变形问题。

（5）研究边坡在地震作用过程中的行为响应，从而研究边坡的动力失稳机制。

在以上几个地震作用边坡的稳定分析研究的主要方面，研究边坡在地震作用过程中的行为响应，进而进行地震力的计算，并判断边坡的破坏面的形状和位置是进行其他研究的前提，而边坡的稳定状态判定和永久变形（位移）的计算则是研究的重点内容。所以地震作用下边坡稳定的研究的根本是对边坡在地震力作用下破坏过程的研究。

地震力对边坡的作用表现为：地震力引起的惯性力和因循环退化引起的剪应力降低。因此，地震作用下边坡失稳可分为：惯性失稳（Inertial instability）和衰减失稳（Weakening instability）［66］。根据边坡的失稳原因，目前对于地震边坡稳定性的分析方法主要基于极限平衡理论和应力—变形分析［67］。

对于边坡惯性失稳，常用的方法有从边坡静力稳定分析方法发展起来拟静力法（Pseudo-Static Analysis），该方法的优势是原理和计算简单，可以很方便得到边坡的安全系数，但计算精度较差；基于滑块理论的Newmark滑块分析法（Newmark Sliding Block Analysis），该方法现已成为求解边坡在动力作用下产生的永久位移的经典方法；在Newmark方法的基础上，Makdisi和Seed对进行改进和发展提出了Makdisi-Seed法；另一类求解边坡惯性失稳问题的方法是数值计算方法，最常见的是动力有限单元法，即通过数值模拟的方法计算边坡的应力及位移，但这一类方法不能对边坡的稳定程度进行有效评价。

对于衰减失稳的边坡常采用流动破坏分析法（Flow Failure Analysis）和变形破坏分析法（Deformation Failure Analysis）。由于边坡的衰减失稳不是本书的研究重点，本书仅作简要介绍。

1.2.2.1 拟静力法（Pseudo-Static Analysis）

拟静力法是人类最早使用的对地震作用下边坡稳定性研究的方法，从19世纪20年代开始，就已用于结构地震稳定分析［66］。其实质是将地震产生的惯性力看作水平方向和垂直方向不变的加速度与坡体自重的乘积，施加于潜在的不稳定滑动体的重心，并将加速度的方向指向使边坡失稳的方向，然后根据极限平衡理论，把所有作用在滑动体的力和力矩进行分解，建立滑动体的静力平衡和力矩平衡方程，进而求解滑动体的安全系数。简单说，就是用振动加速度值所确定相当于静荷载的地震力即拟静力将作用短暂瞬间的荷载同长期作用的荷载等效起来［68］。安全系数的大小与边坡的材料性质（c，φ）、破坏面的形状和位置以及地震惯性力（k_h，k_v）的大小直接相关。在拟静力分析中，边坡材料特性值可通过现场或室内试验测定：而破坏面的形状和位置常根据边坡的地质条件通过经验或工程类比法来确定，且通常简化为直线形、圆形或阶梯形等；对地震力大小的计算，即水平和垂直拟静力系数k_h和k_v的研究较多。

太沙基（Terzaghi，1950）首次提出：对一般性地震k_h=0.1；对破坏性地震k_h=0.2；对灾难性地震k_h=0.5［69］。Seed（1979）对10个地震区国家的14座大坝列出了拟静力设计准则：若安全系数为1.0～1.5，则拟静力系数应为0.1～0.12［70］，他同时指出由柔性土（不会产生高孔隙水压力，循环加载强度降低15%以上的土）建成的大坝在小于0.75g最大加速度作用下的变形小于k_h=0.1（M=6.5）到k_h=0.15（M=8.25）的作用下的变形，即允许13%～20%的峰值加速度作为拟静力加速度［71］。Marcuson（1981）考虑大坝对地震动的特性影响，建议大坝的适宜地震加速度系数应取最大加速度的1/2～1/3［72］。Hynes-Griffin和Franklin（1984）采用Newmark滑块分析计算得出安全系数大于1的土坝采用k_h=0.5a_max/g不会产生危险的变形［69］。

目前没有直接而快速的方法来选择设计拟地震加速度系数，其取值均建立在不稳定体所受的实际加速度基础上（考虑地震动特性影响），取考虑某个系数的实际加速度峰值［66］。

由于拟静力法将地震等动力作用等效为静力作用，可以采用常规的静力稳定分析方法进行分析，极大简化了对滑动体在地震力作用下的破坏过程，计算和理论均较为简单，所以一经提出就被广大的工程设计部门采用，并编写入工程设计规范^［45^，46^，73^］，并在地震边坡分析中得到广泛应用。各国科学家在拟静法原理的基础上，努力寻找更为合理的评价地震作用边坡稳定性的方法，如文献［74］用数值计算方法计算潜在滑动面上的正应力分布，并用此正应力通过显式求解满足所有的静力平衡方程的最小安全系数，制定了用于评价地震作用下简单边坡的设计表，该表格在非地震条件下与文献［75］的表格相同；文献［76］对具有软弱层的垃圾场采用波传播理论和拟静力法进行了分析，提出了垃圾场地地震稳定性评价过程；文献［77］在正常固结土边坡地震分析中采用拟静力法，应用参数分析确定了不同抗剪强度内摩擦角的安全系数并考虑了土体抗剪强度降低的影响；文献［78］将拟静力法用于沿节理面滑动的岩体地震稳定性分析中，进行了地震和稳定分析和永久位移计算；文献［79］将拟静力法应用于加固边坡的地震稳定分析中，利用极限平衡理论考虑不同破坏模式，提出了加固力大小计算公式及与地震力相关的屈服强度表达式。

虽然拟静力由于简单易用，得到世界上广大工程技术人员的认可，但正是由于拟静力法的简单，用振动加速度值所确定相当于静荷载的地震力即拟静力是把作用短暂瞬间的荷载同长期作用的荷载等效起来，在计算的过程中，为了将边坡这样一个超静定系统转化为静定系统，拟静力法采用了较多的假设条件和简化条件，从而夸大了动态力的作用；另外，由于拟静力法求解的方程是边坡的力矩或力的平衡方程，本质上与边坡稳定分析的刚体极限平衡没有区别，故在边坡的刚体极限平衡分析方法存在的问题在地震边坡的拟静力分析普遍存在，如刚体假定、滑动面假定，条间法向和切向作用力大小和方向的假定等，同时由于拟静力法分析的是地震力的作用，故除与刚体平衡分析法相同的假定外，还有一些动力分析的假设条件，如加速度假定（认为边坡的地基加速度和边坡的加速度一致）、力不变假定（认为施加在边坡滑动体中心处的地震惯性力是一个不变的量）以及失稳方式假定（认为安全系数小于1时，边坡发生失稳），这些假定往往夸大了动态力的作用，导致计算结果的失真；此外地震时振动荷载可引起山体的重复加载和松弛变形，也不是拟静力法所能完全模拟的［68］。

拟静力法存在的上述缺陷，已有多位学者予以证明：Wu X.Y.（1991）对黏性土坡采用拟静力法和参数分析，发现该法的不足：在某些条件下，对给定圆心的圆弧滑动面没有最小安全系数，最小安全系数随圆弧滑动面半径增加而单调降低；在某些条件下，滑动体的最小安全系数不能确定，因此得不到临界圆弧滑动面获得到的安全系数过低［80］。H.I.Ling（1997）研究得出水平和垂直加速度很大时垂直加速度对稳定性和位移有重大影响［81］。另外需要注意的是边坡发生破坏的形式不外乎两种：一种是变形，另一种是失稳；而拟静力法考虑的只有失稳一种，通常情况下，由于岩土体的弹塑性性质，地震荷载产生的弹性变形在地震结束后即恢复，但塑性变形在地震后不能恢复，会造成边坡的永久变形或残余变形，拟静力法不能考虑边坡的永久变形对稳定性的影响，是其最大的不足之处，Seed（1973）曾指出，即使边坡的安全系数暂时小于1，也不一定会导致边坡的失稳，而只会产生一定的永久变形，边坡的动力响应行为由边坡的变形控制［82］。

拟静力法另一重大不足是不能考虑地震所引起的砂土液化对边坡稳定的影响。地震液化是进行地震动力分析不可缺少的重要内容之一，地震液化所引起的地面沉陷是地震引起的土体的竖向残余变形，它是衡量工程结构是否丧失功能的主要指标之一，而且是提出抗震对策措施的重要依据。地震引起的永久性地基沉降会造成建筑物的严重破坏，1957年墨西哥城地震和1986年Miyagioki地震就是典型的例子。震陷大部分是由液化所引起的，因此对液化震陷预估具有很重要的工程意义。

地震荷载作用下所产生的沉降大体可分为两类：①动力荷载作用下，土体的不排水残余应变引起的瞬间变形；②地震期间，土体中产生的超孔隙水压力的消散所引起的再固结变形［83］。

1.2.2.2 Newmark滑块位移法（Newmark Sliding Block Analysis）

Newmark滑块位移法是以Newmark（1965）提出的屈服加速度的概念为基础的［84］，Newmark指出：边坡的稳定性取决于边坡的变形，变形由应力时程变化决定，瞬时惯性力如果足够大，导致边坡的安全系数暂时小于1，将会导致边坡产生永久位移，但这种变形并不是随时间持续增加的，当地震加速度降低或反向施加于边坡时，边坡的变形就会停止，所以边坡的永久位移实际上是边坡在地震过程中的累积位移，如果在地震过程中，岩土体的强度没有显著降低，边坡将不会发生进一步的严重变形［84］。

Newmark的这一理论澄清了传统拟静力法认为决定边坡稳定性的是边坡的最小安全系数，边坡是否失稳决定于边坡内部最大应力的错误看法，而是将地震力作为短暂的变向荷载，考虑了地震引起的边坡变形，并提出了屈服加速度和累积位移的概念，推动了地震荷载作用下边坡稳定性分析方法的发展。

Newmark将边坡的潜在滑动体与位于斜坡上的物体类比，认为当安全系数小于1时，潜在滑动体不再保持平衡，破坏体在不平衡力的作用下加速下滑，并利用此类比发展与任何地震动相关的边坡永久位移的预测方法。该模型尽管很近似，但比传统拟静力分析提供更多信息，比有限单元法实用性更强。

Newmark标准滑块分析基于如下假定：潜在滑动体是完全刚性的；滑块与滑面间是完全塑性应力—应变关系；破坏面是平面或圆弧面（图1-2）。其基本方法是将超过可能滑动体屈服加速度的加速度响应作两次积分，即可估算边坡的有限滑动位移（图1-3）。

位于坝坡上的滑块在未滑动之前，安全系数F_s&gt；1。在地震荷载下，当其处于临界滑动状态时，即图1-2中滑块的下滑力等于摩擦力时，F_s=1。此时的地震加速度成为屈服加速度a_y。目前一般沿用静力法或拟静力法加以计算，即将地震惯性力作为地震荷载，因此屈服加速度也可用屈服地震系数或称加速度系数k_y来表示。滑块开始滑动时，F_s&lt；1，此时滑块产生速度和位移。如果地震有多个脉冲，第一个脉冲的加速度超过k_y后，滑动开始并产生速度和位移。当这个脉冲的加速度减小至小于k_y时，滑动速度减小，直至停止滑动。下一个地震脉冲又重复这一过程，如图1-3所示。向下滑动的位移被认为是不可恢复的，所以位移是逐次累加的，此即所谓永久变形。但是当地震脉冲产生的加速度是反向的，并且超过反向的屈服加速度时，所产生的残余可以是反向的，即负的。正负位移可以抵消。一次地震的脉冲次数是有限的，所产生的残余位移可根据地震加速度而变化。在确定了屈服加速度后，可通过积分求得各次脉冲的滑动残余位移及其累积量，直至地震结束。因此，这样产生的滑动位移的总量是有限的，不可能像静力失稳那样出现很大的位移或造成严重破坏。

滑动位移的计算可归纳如下几个步骤。

1.确定屈服加速度a_y

屈服加速度a_y的含义是使坝身沿着某一可能滑动面的滑动安全系数等于1时的加速度，它与重力加速度的比值成为屈服地震系数k_y。屈服地震系数k_y值与坝身几何尺寸、土料不排水强度、可能滑动体的位置等因素有关。为了计算k_y，需先确定土的屈服强度。当超过这一屈服强度时，即产生永久塑性变形，其大小与应力循环次数及频率等因素有关。

从图1-4（a）可以看出，当应力（包括静应力和周期应力）达到不排水强度80%时，100次循环作用几乎不产生永久变形。当应力达到95%不排水强度时［图1-4（b）］，10次应力循环即可产生很大的永久变形。所以，对于填筑材料，在缺乏试验资料的情况下，屈服强度可取静不排水强度的80%，或查相关的试验资料，如图1-5所示。根据确定的屈服强度，即可用一般的圆弧滑动法算出滑动面上的安全系数恰好等于1的屈服加速度a_y或屈服地震系数k_y。

2.确定平均地震系数k_av（t）

平均地震系数k_av（t）定义为作用于滑动面上剪切力的水平分量与可能滑动面上重力W之比，即

平均地震加速度a_av（t）为平均地震系数k_av（t）与重力加速度g的乘积，即a_av（t）=k_av（t）g。

F（t）可用剪切楔法近似确定，也可用有限元法确定。用有限元法确定各单元的应力时程曲线后，可进一步计算滑动面上各单元的正应力和剪应力，如图1-2所示。任意时刻在滑动面上的水平力为

式中：n为滑动面上的单元数。

平均地震系数k_av（t）取决于可能滑动体在坝体中的部位。如果以可能滑动体的滑出点与坝顶的距离y表示它在坝体中的位置，则平均地震系数应写成k_av（y，t）。以k_max，y表示k_av（y，t）的最大值，成为滑动深度为y的最大平均系数，以k_max，0表示y=0（即坝顶的地震系数）的最大值。则y/H与k_max，y/k_max，0的关系如图1-6所示。坝顶地震系数k_max，0可用有限元法通过求取坝顶加速度而获得，进而又图1-6上的曲线推求各种深度y的k_max，y。

3.确定有限滑动位移

Newmark认为，当平均加速度在可能滑动体中产生的惯性力的方向与滑动面上静剪切力的水平投影方向相同时，如果k_av（t）≤k_y，滑动体不产生滑动；如果k_av（t）≤k_y，滑动体产生滑动，滑动方向与静剪切力方向相同。滑动加速度的水平分量a_x（t）等于平均地震系数k_av（t）与屈服地震系数k_y之差乘以重力加速度g，即

在整个振动过程中，只有正半周或负半周（即惯性力方向与滑动面上静剪切力水平投影方向的相同半周）才能使滑动体滑动。每次滑动的水平位移按下式计算，即

在整个地震期间，总的水平位移δ应是每次滑动水平位移δ_i之和，即

式中：n为整个地震期间的滑动次数。

最大的总水平位移δ_max（y）通过试算确定。假定不同的y即可求得相应的δ（y）。

自Newmark提出这一方法的基本思想后，大量研究者就这一方法在地震波形对计算结果的影响方面、算法扩充方面及实际应用方面展开了研究，下面总结一下国内外学者对这一方法的研究进展。

由于Newmark滑块分析采用的地震波形为简单波形，研究者发现滑块在矩形波、正弦波和三角波作用下的永久位移与周期平方成比例，在很多学者研究了边坡位移与a_y/a_max之间的关系后，得到了很多近似表达式。采用滑块分析方法分析真实地震动作用下边坡的永久位移，发现a_y/a_max大于0.5时永久位移与a_y/a_max间曲线形状与正弦波和三角波相同［85-87］。Makdisi和Seed（1978）利用Chopra（1966）的方法得到了平均加速度，并利用Newmark分析法计算了土坝和坝堤在地震作用下的永久位移，在通过有限元分析结果的简化假定和此类结构的剪切楔分析，提出了一种简化预测永久变形的方法，通过几座实际和虚拟大坝在不同震级的实际地震动和人工地震动的作用分析，绘制了不同震级下的永久位移［标准化位移u/（a_maxT₀）］与a_y/a_max的关系图［86］。Ambraseys和Menu（1988）发现当a_y/a_max值较低时，其曲线形状是否考虑受非边坡位移的影响［87］。

由于用峰值加速度来描述地震动有一定限制，一些学者还采用其他地震动参数来进行边坡位移预测，Crespellani T（1998）引入地震破坏趋势因子（PD）来考虑控制边坡稳定的主要影响因素，通过分析310条实际水平地震波，提出了水平震动的Newmark刚性滑块与PD间关系的经验公式［88］。

一些研究者对Newmark分析法进行了扩充，Giovanni Biondi（2000）在饱和无黏性土土坡中采用Newmark分析法并考虑了地震引起的孔隙水压力和剪切强度降低对位移的影响，应用此模型利用谐波和实际地震波进行了参数分析，得出各参数的影响。位移反应显示饱和无黏性土坡地震稳定分析若忽略孔隙水压力会低估其地震反应［89］。

Newmark分析法在很多实际边坡中得到了应用。Taylor Larry R（1995）对有黏性土界面的加州某固体垃圾场的覆盖层进行地震稳定分析，根据界面剪切实验得到的屈服加速度，采用Newmark分析法对土与界面上屈服位移进行了计算，对比不同覆盖物的计算屈服加速度和变形分析结果，获得了不同覆盖物的预期结果并进行了比较［75］。Kanthaasamy等人将离心模型试验、完全动力有限元分析程序、传统的拟静力分析和Newmark分析法有效组合，提出了安全经济的设计方法。

自Newmark提出滑块位移的基本算法后，一些研究者就这一方法的具体内容提出了各种算法，Makdisi和Seed的简化方法是其中最有代表性的一个。Makdisi和Seed的简化方法适用于填筑密实的，由压实的黏性土构筑的土石坝。其屈服加速度可按常规圆弧滑动条分法通过拟静力分析求的。由于该方法是建立在许多实例计算结果之上，其优点是简单方便，但由于所用的地震记录、坝体材料的动力性质指标以及坝高范围、坝坡比等资料均十分有限，故而有一定局限性。

国内关于Newmark分析法的研究也很多，其中具有代表性的有：王思敬（1977）基于Newmark分析法提出的边坡块体滑动动力学方法［90］；王思敬、张菊明（1982）通过试验，提出运动起始摩擦力和运动摩擦力的概念，并在振动台上测得花岗岩光滑节理面的动摩擦系数和运动速度的关系［91］，在此基础上王思敬、薛守义（1989）、张菊明（1994）又分别推导了楔形体情形下和层状山体情形下的三维动力反应方程式^［92^，93^］。

近年来还开展了有关Newmark分析法可靠性的试验研究。Wartman Joseph（1998）为评价Newmark分析方法的应用范围和准确性，在振动台上进行了边坡物理模型的强震实验，发现边坡中出现了较大的平动/转动和较小的破坏运动，应用Newmark分析方法进行能够反分析得出：①测量位移介于采用峰值和参与剪切强度计算结果中间；②真正的变形时间超过堤坝地震反应结束时的时间；③地表运动时剪切面位移为滑动体内部应变之和［94］。

1.2.2.3 数值分析法（Numerical Analysis Method）

近年来，随着计算机的逐步发展，以往不可能靠手工完成的工作现在已经能够依靠计算机的强大计算能力解决，在边坡稳定分析的领域也是如此，不论是静力问题还是动力问题，目前都已有较为成熟的计算软件，而且计算方法各异，如有限单元法（FEM）、离散单元法（DEM）、差分方法（FLAC）、边界元（BEM）方法等。

将有限单元法应用于地震作用下土体的反应分析始于20世纪60年代^［95^，96^］，有限单元法又称为整体变形分析，它是将地震前后坝体及其地基假定为连续体，按照连续介质的有关理论来进行考虑，对地震力作用下边坡有限元网格中各单元永久应变进行积分得到边坡的永久位移。

用有限元法计算体积变形所产生的永久位移时，常用的方法有初应力法［83］、应变趋势法、刚度折减法（修正模量法）、等效节点法和考虑土体非线性非弹性应力—应变特性的非线性方法［66］。

初应力法就是把每一时段各单元产生的孔隙水压力增量Δu_g看作初应力，然后把它化作等效的节点荷载，再引入到静力平衡方程的荷载项内，解出节点位移，即为永久位移。Seed（1973）通过对1971年San Fernando地震期间大坝失稳的调查，从线性和等效线性有限单元分析的结果发展了一种用于估计地震引发的边坡变形方法［82］；Lee（1974）和Serff（1976）提出了一种新的用于计算边坡永久位移的方法，此方法在计算应变趋势时对土的刚度进行了折减，可用于线性和非线性模型估计水平方向和垂直方向的位移，但仍是一种近似的方法［97］。

土体动力分析有限单元法的总体思路和静力情况基本一样，不过由于荷载和时间有关，相应的位移、应变和应力都是时间的函数，因此在建立单元体的力学特性时，除静力作用外还需考虑动荷载及惯性力的阻尼作用，在引入这些量的影响后，就可类似静力有限单元分析过程建立单元体和连续体的动力方程，然后采用适当的计算方法求解［2］。

由于有限单元法不但可以应用总应力法，而且还可以以有效应力法为基础，可以考虑复杂地形、土的非线性、非均质性、弹塑性及土中孔隙水等复杂条件对地震期间边坡稳定的影响，能够深入分析土的自振特性及土体各部分的动力反应，因此有限单元法已成为边坡动力响应分析的重要方法之一。

应该注意到的是，虽然有限单元法是进行动力响应分析的一个很好的方法，但由于有限单元法自身基于连续介质，对于均质土体来说其分析结果可靠，但对于材料介质连续性很差的岩质边坡来说，有限单元的分析就不那么准确了，而且有限单元法不能分析对边坡稳定有重大影响的节理面、裂隙面上的动力响应情况也是其重大的缺陷。针对这种情况，科学家们发展了新的数值计算，以解决该问题，目前在岩土工程领域较为流行的美国Itasca公司的UDEC和Flac软件是其中的代表之作。

Flac软件采用的是有限差分的拉格朗日算法，能解决岩土工程中常见的非线性大变形问题，该方法遵循连续介质的假定，利用差分格式，按时步积分求解，随着构型的变化不断更新坐标，允许介质产生较大的变形。

UDEC软件采用的是与有限单元法和有限差分有很大区别的离散单元法（Discrete Element Method，DEM），是Cundall P.A.（1971）提出的^［98^，99^］，于20世纪80年代引入我国（王泳嘉等，1986）［100］。目前这一方法已在数值模拟与工程应用方面取得了长足的发展。二维与三维可变形离散单元法亦已问世，并在岩石工程与岩石力学中得到日益广泛的应用。值得注意的是，由石根华与Goodman提出的DDA方法（Discontinuous Deformation Analysis）也属离散元模型的一类，它在求解块体变形与用应变模态组合来代替Cundall的有限差分格式。与此同时Williams与Mustoe则用正交模态组合来近似岩块的变性特征，仅需用几个低阶模态就能描绘岩块的复杂变性。为了模拟颗粒体的力学行为，Ghaboussi与Barbosa提出了一种特殊的刚体离散元。在离散元与其他数值计算方法的耦合技术方面，Lorig等人提出了一种离散元与边界元的耦合模型模拟裂隙岩体与完整岩体的组合系统。

离散单元法的一个突出功能是它在反映岩块之间的接触面的位移、分离与倾翻等大位移的同时，又能计算岩块内部的变形与应力分布。因此，任何一种岩体材料的本构模式都可引入到模型中，例如弹性、弹塑性、黏弹性或断裂等均可考虑。该方法的另一个优点是它利用显式时间差分解法（动态松弛法）求解动力平衡方程，这一方法用于求解非线性大位移与动力稳定问题具有先天的本能，用于求解岩质边坡对地震荷载作用的响应是目前最为适合的方法。

在求解动力问题方面，Lemos发展了离散元动力分析程序，可用于计算岩石工程的地震响应问题以及混凝土大坝—地基动力相互作用问题，后者还引入了考虑辐射阻尼的投射边界；Dowding等人利用离散元与有限元的耦合系统分析了地下洞室与围岩介质的动力相互作用，其中洞室附近的节理岩体用刚体离散元模拟，远场完整岩石则用有限元离散，计算结果还与模型试验做了比较。此外，Iwashita与Hakuno研究了岩石陡坡的动力稳定过程，Meguro与Hakuno利用离散元模型研究混凝土结构的地震破坏过程，也取得了一定的进展。

在我国，离散元方法的研究与应用开始于20世纪80年代中期，王泳嘉［98］首次在我国应用这一方法于节理岩体的数值分析中，研究了放矿的数值模拟与自然崩落机制；并应用离散单元法与实验研究结合指出边坡失稳过程中应力核存在的现象［101］。魏群研究了椭圆颗粒的离散单元法［102］，并进行了模型试验验证。随后，王光纶、张楚汉提出了利用离散元分析岩质边坡的失稳判别标准和安全系数定义［103］。王光纶、张楚汉、彭冈用刚块试验验证了刚体离散元的滑移原理［104］。鲁军、张楚汉则完成三维离散元的自动剖分系统并初步应用于拱坝坝肩稳定分析［105］。李世海等人利用离散单元法对三峡永久船闸高边坡进行了数值模拟［106］，赵坚等人用UDEC模拟爆炸波在节理岩体中的传播［107］，谭文辉等人用离散单元法对静力边坡演化过程中非线性现象进行了研究［108］。目前利用这一数值方法研究岩石工程与岩石力学的论著和文献日渐增多，文献［109-112］从各个方面论述了离散单元法的应用，在此不再赘述。

此外其他的数值分析方法如非连续变形方法（DDA）、流形元法、边界元法以及它们之间的耦合关系，限于篇幅不再介绍。

1.2.2.4 概率分析方法（Probability Analysis Method）

在地震边坡稳定分析中存在很多不确定因素，如输入地震动的随机性，边坡材料特性的随机性。和边坡的可靠性分析一样，这些随机性和模糊性对计算结果有很大影响，在分析时有必要考虑这些随机性，因此出现了地震边坡概率分析方法。Halatchev Rossen A.（1992）提出了一种用于堤坝和边坡地震稳定性分析的概率方法，该方法建立在Sarma解的基础上，考虑了地震力的水平和垂直分量，即地震力具有任意倾角；土体剪切强度参数假定为正态分布，采用Monte-Carlo模拟；破坏概率由地震系数确定；将地震系数视为一个随机量，并将其应用于露天矿边坡分析中［113］。AI-Homoud A.S.（2000）利用安全系数和边坡的临界位移提出了地震力作用下土体边坡和堤坝的概率三维稳定性分析模型，考虑了如下的不确定性因素：实验室和现场所测的剪切强度参数存在大差异的不确定性；发生地震和地震所引起的加速度的随机性。同行还考虑了空间上变化以及土参数间的关系。提出了5种基于极限值下的非超越概率的地震位移概率模型，并提出了基于安全系数的三维动力边坡稳定分析的概率模型，编写了PTDDSSA计算程序，将这些模型应用于著名滑坡Selest Landslide中，发现震源距离和震级对地震引起的位移、破坏概率、二维和三维安全系数影响很大［114］。

1.2.2.5 衰减失稳流动破坏分析（Flow Failure Analysis）

流动破坏发生于剪切强度低于维持边坡平衡的静力剪切强度时，发生迅速而无任何警告且产生巨大变形。该分析方法分两步，首先确定是否发生流动破坏，其次是确定流动破坏的范围。

潜在流动滑块不稳定性常采用传统的边坡稳定性分析方法来评价，采用地震结束时的强度，若稳定分析中得出流动破坏可能发生，则流动破坏变形分析可确定破坏影响区范围，忽略流动滑体滑动前的小变形，流动滑体的变形可通过极限平衡、流体力学和应力—变形分析来粗略估计。Lucia和Iverson在这方面做了大量的工作［67］，另有有限元计算程序TARA-3FL能够在单元初始液化时将边坡单元的强度降至残余强度［115］。

1.2.2.6 衰减失稳变形破坏分析（Deformation Failure Analysis）

衰减失稳变形破坏发生于地震引起的剪切应力暂时超过土体剪切强度时，类似于惯性破坏，变形破坏产生一系列永久变形“脉冲”，在地震结束时即停止，变形破坏比流动破坏产生的变形要小，但也可能产生大的破坏，侧向扩展式变形破坏最普遍的形式。近几年来，很多研究者发展了多种方法评价变形破坏产生的永久变形，但因变形破坏机理复杂，预测变形主要凭经验。Hamda［67］，Byrne［116］，Bartlett［67］在这方面做了许多工作。