河南理工大学土木工程毕业论文详解

边坡工程可靠度分析方法研究摘要：此文以可靠度理论为基础，考虑土体的随机特性，提出了应用于边坡工程可靠度分析的蒙特卡罗法的基本思路，岩土工程中基本变量大多数都为互相关，通常土性参数、之间存在负相关性，因此在计算边坡稳定失效概率时要考虑其相关性。Monte-Carlo模拟法是求解失效概率和可靠度指标较为精确的方法。此文用Monte-Carlo求解相关空间中失效概率和可靠度指标的计算模拟方法，并编制了相应的程序。关键词：边坡工程；可靠度指标；相关变量；蒙特卡罗方法Abstract：Basedonthetheoryofreliabilityandtheconsiderationofstochasticpropertiesofsoil,thebasicideaabouttheMonte-Carlomethodusedtoanalyzethereliabilityofslopeengineeringispresentedinthisdissertation.Thefoundamentalvariablesingeo-engineeringarecorrelated.Itisnecessarytotaketheircorrelationintoconsiderationwhichcalculatingthereliabilityindex.Monte-Carlosimulationisanaccuratetechniqueforestimatingfailureprobabilityandreliabilityindexoftheslopestability.Inthisdissertation,Monte-Carlomethodisappliedtocalculateandfromaspaceofcorrelatedvariables.Theappliedprogramisdevelopedtocalculatethecorrespondingreliabilityindex.Keywords：slopeengineering；reliabilityindex；correlatedvariable；Monte-Carlomethod目录一、绪论 31.1边坡工程可靠度分析的必要性[1] 31.2边坡工程可靠度分析研究现状 7二、基本理论 102.1边坡工程可靠度分析的基本方法 102.2用瑞典法建立极限状态方程 132.4产生给定随机变量分布类型的随机数 162.5相关变量条件下蒙特卡罗模拟法 17三、工程实例 21四、结语与展望 23附录 26参考 30一、绪论1.1边坡工程可靠度分析的必要性[1]边坡工程可靠度分析是上世纪六十年末代发展起来的边坡工程状态的新方法。它把边坡岩体性质、荷载、地下水、破坏模式、计算模型等作为不确定性量(图1-1)，借鉴结构工程可靠度理论方法，结合边坡工程的具体情况，用可靠指标或破坏概率(或可靠概率)描述边坡工程状态的理论体系。与传统的确定性理论比较，它能更好地反映边坡工程的实际状态，正确合理地解释许多用确定性理论无法解释的工程问题，因此引起国内外有关人广泛关注。边坡工程可靠度分析之所以能够产生和发展并逐渐为广大工程技术人员所接受，是由于工程设计和分析中不可避免地要涉及到大量的不确定性因素，如岩体性质、荷载、破坏机理等。这些不确定性因素按其成因大致可以分为四大类：物理不确定性、统计不确定性、模型不确定性和人为过失造成的不确定性等。在具体的边坡工程分析中，可以概述如下：（一）岩层(土层)剖面与边界条件的不确定性在边坡分析中，首先要知道边坡中包含的各岩层及其厚度变化等情况，从理论上讲，如果有足够多的测点或钻孔信息是可以精确地描述它们的。事实上，由于岩层出露条件和钻孔工作量的局限，人们是很难精确描述它们的，必然会加入一些人为的判断因素，在有限个测点之间，进行主观推断，给出各岩层的边界。图1-2是某坝基C5构造剪切带17.0～40.0m段的厚度变化情况，可见同一剪切带的厚度变化是不规律的。岩层剖面和边界的不确定性必然导致破坏模式和计算模型等的不确定性。（二）岩土性质的固有变异性岩体和土体都是非均质的各向异性体，各点间的性质往往有较大的差异，由于试验样品不足及其它各种因素的影响，其变化程度也表现出随机性，同一组试样在相间试验条件下测定其强度，结果也表现出一定的离散性，其中摩擦角的值可能离散，粘聚力值可能离散，给定正应力下的抗剪强度值变化约为。国内某露天煤矿北帮风化岩的抗剪切试验表明，粘聚力、内摩擦系数均服从正态分布．粘聚力的变异系数(变异系数＝标推差／均值)为0.18～0.54，内摩擦系数()的变异系数为0.19～0.86。人们在进行岩土性质的变异性研究时往往不注重取样的空间位置，只要是同一岩层的试样，对其试验结果进行笼统的统计分析，得出变异性的有关指标。事实上，由于岩土体是三维空间实体，在进行取样时，无法在精确的同一取样点获得多于一块的试样，只能在别的位置再获得试样，这样统计出的岩土性质变异性结果中，更大程度上包含有空间变异性的成分。因此，近年来人们引用了随机场、地质统计学和贝叶斯马尔可夫过程理论等研究岩土性质的空间变异性，并取得了一定进展。图1-1边坡工程的不确定性图1-2c5剪切带的厚度变化（三）试验样品数量不足从理论上讲，岩层的边界和性质等并非完全是随机的，而是确实地存在于三维空间实体中的，而且至少在一定范围内这些性质和特征是相关的，可以用某些规律描述，如果进行无限多的勘察工作，情况是可以确知的。但是实践中，无论是地质工作还是试验工作都有限，而由少量的资料去推断全部情况是不可避免的，这样就不得不借助概率统计理论的工具。为方便起见，将岩层的性质和特征看作随机的，试验样品数量的影响程度也是不确定的，人们很难确知样品数量与真实试验结果之间的确切关系。（四）外加荷载大小和分布的不确定性对矿山边坡工程来说，常见的外加荷载有地应力、地下水、地震、爆破震动、运输等生产设备、降雨等，这些荷载很难用确定指标描述，它们都是随机变量，并见其分布的类型和特征参数往往又都具有不确定性。其它的边坡工程，如建筑物的基坑边坡、路基边坡及各类开挖边坡或自然边坡，同样也会遇到荷载的不确定性问题，并且荷载的不确定性往往又会导致岩土表现出的性质不确定性，如水理性质、动力性质等。（五）勘测取样方法与试验方法的误差在取样及测试过程中对试件的扰动，试验环境条件以及试验方法的不一致等，都使结果产生差异。（六）计算模型的不确定性与其他岩土工程一样，在定量分析中，边坡工程依据对岩土地质模型的深刻观察和分析来建立边坡物理力学模型。由于工程岩土的复杂性和人们分析判断能力的局限性及模拟条件和原型的差异，必然导致结果的不确定性，致使目前许多物理模拟试验结果仅能用来做机制分析或趋势分析，而难以为设计提供准确的依据。同时为了使量化计算能够实现，人们不得不对复杂的边坡破坏机理、模式及其环境条件等做出某些简化和假设，给山一些确切的函数关系，以此作为一切计算、分析和评价的基础。实际上，在简化和假设的过程中，已经引入了不确定性，因此不能说计算模型就是精确的，即使是输入相同的变量及其分布，由于采用不同的计算模型(当然是合理的)，也会得出不同的结果。出于岩土工程中的不确定性，使人们对采用传统的确定性方法研究岩土边坡问题，即用安全系数表示安全程度产生了疑问。边坡安全系数是岩土许多参数的函数，既然这些参数具有不确定性，那么把安全系数作为确定值来判断边坡工程的安全程度显然是不够合理的。实际上，在岩土工程实例中，也确实出现过安全系数大于1的边坡却发生了破坏，而安全系数小于1的边坡却依然稳定的情况。为了最大限度地消除岩土工程中的不确定性，人们试图从提高岩土测试技术和计算技术的精确度人手，这种努力可以取得一定效果。然而，局部试验的精确性、确定性并不能消除岩土性状宏观判断的随机性和模糊性，靠无限度提高单项试验的精度、规模和完善确定性计算方法是不够的，因此用较简单的测试手段，对岩土进行大数量的信息采集，应用和发展可靠度理论、方法，以提高边坡工程状态客观判断的精度，就显得必要和有意义了。在边坡及其它工程的设计和研究中，人们有意识或无意识地都考虑了岩土中不确定性因素的影响，例如在边坡临界安全系数的选取上，包含了对不确定性因素的预计。通常认为对于土坝，临界安全系数取1.25即可满足要求；对于煤渣堆或尾矿坝，或较复杂的基坑边坡，由于材料的特性可能变化无常，则临界安全系数一船取1.5对于矿山而言，非运输帮边坡的临界安全系数通常取为1.3，运输帮边坡的通常取为1.5。从理论上讲如果能够充分认识和量化影响边坡稳定的因素，对于任何边坡，其临界安全系数都应取为1。可是由于无法真正地量化所有因素，以及对输人参数和计算模型精确性的顾虑，人们不得不凭借经验和工程实例的总结、推断，将临界安全系数取为大于1的某个值，具体取值要视对工程岩土的认识程度、环境条件及工程的重要性等而定。实际上临界安全系数的取值包含了对不确定性因素的直观预计。另一方面，虽然安全系数是一个由确定性方法获得的一个确定数值，但它本身又是一个不确定性量，例如一个安全系数为1.4的边坡，并不意味着有140%的安全程度，安全系数为1.5的边坡，也不说明它较安全系数为1.4的边坡要高10%的安全程度。这就是说，安全系数大小的本身往往不能定量表示安全富裕的程度，只能说明具有高安全系数的边坡要较低安全系数的边坡更稳定一些。上述这些事实说明，对于边坡及其他岩土工程的设计和分析，单凭一个笼统的安全系数很难定量地评价工程的安全程度，也无法比较采用不同的计算方法、取不同的安全系数时，边坡工程是否具有相同的安全储备。随着现代科学技术的高速发展，边坡工程已逐渐大型化和复杂化，如人们经常会遇到高大、复杂和环境限制严格的露天矿边坡、铁路边坡、基坑边接、坝体边坡、河谷岸坡等。为了客观地考虑到岩土工程中的不确定性因素和评价边坡工程的安全程度，加之相关学科的发展和完善，推动了可靠度分析方法和理论在岩上和边坡工程中的研究和应用，目前也逐渐为广大工程技术人员及管理人员所接受。1.2边坡工程可靠度分析研究现状边坡的可靠度分析是基于对岩土性质等的不确定性认识，随着认识逐步深化，可靠度研究得以逐步发展。自上世纪六十年代末以来，国际上对岩土不确定性的理论和应用研究以及边坡工程等的可靠度研究非常活跃，研究范围也相当广泛。原则上讲，岩土工程的可靠度方法可以应用于任何岩土工程问题，但由于各种因素的影响，目前比较常用的有以下几个方面：1)岩土工程结构物的设计和优化分析，尤其是边坡、挡土墙、海洋平台和桩基的稳定分析研究较多；2)岩土结构物的性状预测，如预测地基的孔隙水压和固结特性，沉降和水平位移特性等；3)地质勘探资料的统计分析、概率地质剖面的模拟和各种参数的概率分析及其表达、岩体节理的随机模拟等；4)各种数据和资料的检验、修正、补充和更新，其中除古典的分析方法以外，贝叶斯方法在近年得到了更广泛应用。1971年在香港召开了第一届“统计学和概率论在土工结构中的应用国际会议”；1980年在澳大利亚召开了“概率方法在岩土工程中的应用学术讨论会”；1983年召开的“第六届岩土力学中的数值方法国际会议”，已将可靠度方法作为一个专门的内容讨论。1984乍在加拿大多伦多召开的第四届和1988年在瑞土治桑召介的第五届“国际滑坡会议”均由澳大利亚Wollongong大学的R．N．Chowdhury教授对过去几年来边坡工程可靠性研究的情况进行了系统的回顾与评价。其它的国际系列会议，如“岩石力学国际会议”、“土力学与基础工程国际会议”、“露天采矿岩土稳定国际会议”、“美国岩石力学年会”以及国际上有关岩土和结构的著名杂志，如“国际岩石力学与采矿科学(InternationalJournalofRockMechanicsandMiningScience)”、“岩土工程学报(JournalofGeotechnicalEngineering)”、“岩土力学(Geotechnique)”、“加拿大岩土学报(CanadianGeotechnicalJourna1)”、“结构安全(SafetyStructure)”等都经常发表一些有关岩土不确定性及边坡工程可靠性方面的文章。尤其是近几年来，有关边坡的文章中，很大一部分那是可靠度研究和应用方面的。这些文章中除了一些理论研究以外，工程实例研究也占有相当大的比重，这说明边坡工程可靠度研究已开始运用到工程实践中。近20年来，边坡工程的可靠度研究在国内发展迅速，涉及的范围也相当广泛，从山体边坡、坝体、铁路、土建到矿山边坡等许多领域。在露天矿边坡方面，1982年冶金部马鞍山矿山研究院负责完成的“攀钢灰石矿边坡可靠性分析与经济评价”是可靠性理论方法在国内矿山边坡工程实践中的首次应用，1990年12月该院负责完成的“太钢尖山铁矿露天矿边坡优化设计方法”课题将当时的边坡可靠性理论和方法进行了比较全面的研究和应用。边坡工程的可靠性研究可追溯到60年代末，在研究初期，认为岩土参数，如粘聚力，内摩擦角、容重等均是互不相关的随机变量，分布上多采用正态分布或对数正态分布，也有采用伽玛分布的。求解上多采用蒙特卡罗方法或中心点法。随着对岩土及边坡可靠度研究的深入，人们逐渐发现各岩土参数之间具有某种相关性。Nguyen(1985)指出：“许多岩上参数的值都在不同程度上，以不同的方式受到同一过程的影响，因此任一对变量之间可能是正相关或负相关的，例如在渗流问题中，孔隙度与渗透性之间是正相关的，而粘聚力与内摩擦角之间是负相关的”。“粘聚力与内摩擦角之间的相关系数通常被假定是零，然而某些文献已经指出它们是负相关的”。在可靠度分析的方法上，目前广泛采用的仍然是蒙特卡罗随机模拟方法、中心点法、验算点法和概率矩点估计法。以上概要地总结了近年来边坡工程可靠度研究的工作和成果，这些工作对于建立边坡工程可靠度分析的基本理论、推动工程实践应用起到了积极作用。但是，客观地说边坡工程可靠度分析目前仍主要处于理论研究和实用探索阶段，没有建立一套设计和分析的规范，可靠度分析往往只作为确定性方法的一种补充和参考。有关可靠度分析结果的进一步处理和分析也缺少深入的研究，还没有达到广泛实用的阶段，问题在于方法的本身有许多不成熟的地方。比如，目前对岩土工程中种种不确定件的认识和研究仍显得不够，要定量地确定它们还有许多事情要做；又比如地质勘探工作中贯穿随机的概念和岩层资料分析中的随机处理也显得十分薄弱，有的甚至基本上还未涉及，如果缺少第一手资料的随机处理，那么边坡的可靠度分析就失去了必要的基础。因此，对岩层不确定性的认识、描述和工程实践是相当必要的，而且也是一项十分艰巨的工作。其次，目前的分析方法还比较复杂，影响了它的推广应用。复杂的方法可能会有较高的精度，但经济性和可接受的程度就会下降。尽管如此，边坡工程可靠度分析的方向是正确的，它能真实地考虑到客观和人为的一些不确定性因素，并已逐步运用于工程实践，从事研究和接受这一方法思想的人员也在增加，因此具有生命力。再者，可靠度分析方法与常规的确定性方法不是互相排斥的，它们可以互为补充。事实上，常规方法长期积累的工作经验和许多分析模式、计算方法，仍是可靠度分析方法的一个基础和重要组成部分。二、基本理论2.1边坡工程可靠度分析的基本方法可靠度是个模糊的概念，为方便分析比较，提出了可靠度指标，在可靠度分析中，结构的极限状态一般由功能函数加以描述。当有个随机变量影响结构的可靠度时，功能函数为:式中为影响结构安全度的随机变量。当时，结构处于可靠状态；当时，结构处于极限状态；当时，结构处于失效状态；而可靠度指标的定义式由Cornell提出：式中：—为变量的均值；—为变量的方差。对多个相互独立正态分布的随机变量，有：可靠度计算的主要结果就是要求出这个可靠度指标。由上式可知，可靠度指标的表示很简单，但由于各随机变量分布的相关性和随机性，该指标的获得困难重重。边坡工程的可靠性是指边坡工程在规定的使用条件和设计的服务年限内，完成其设计功能的能力。露天矿边坡的可靠性可理解为，在矿山设计的服务年限内，边坡工程保证采矿、剥离、运输等作业正常进行的能力。边坡工程的可靠性通常用可靠指标或破坏概率表示。边坡工程可靠度分析的主要任务就是求解某一具体边坡的可靠指标或破坏概率，以此对边坡工程的状态进行综合评价，或根据工程的性质、所处的环境等结定边坡可接受的可靠指标或破坏概率对边坡进行优化设计。由于描述边坡状态的极限状态函数通常是输入基本随机变量的非线性函数，因此给精确求解带来一定的困难，这就产生了许多近似求解方法。目前常用的方法有蒙特卡罗(Monte-Carlo)随机模拟方法、概率矩点估计法、一次二阶矩法(F0SM，包括中心点法和验算点法)和随机有限元法等。本文将结合工程实例利用蒙特卡罗模拟法求解，后文将重点论述。其他的计算方法简单介绍一下。通过边坡的破坏模式和破坏，建立了极限状态函数以后，将边坡工程可靠度分析抽象为对极限状态函数的求解。一般的极限状态函数有：式中是个具有一定分布的基本随机变量，如粘聚力、内摩擦系数、地下水压、容重、岩层厚度。对极限状态函数求解的目的就是通过随机变量的数字特征，求出函数的数字特征。由于往往是的非线性复杂函数，所以通过的数字特征精确地求出函数的数字特征是困难的，这就产生了许多近似计算方法，并且在工程中得到了广泛应用。（1） 一次二阶矩法所谓一次二阶矩法是针对功能函数为变量的线性函数，以变量的一阶矩和二阶矩为概率特征进行可靠度计算的一种方法。对于非线性功能函数，一般在某点进行泰勒级数展开并进似地取其一次式，使功能函数线性化，然后再用一次二阶矩法计算可靠度。在可靠性设计的初期，由于各个随机变量的分布规律难以确定，而这些变量的一阶矩和二阶矩则容易得到。对于非线性的功能函数，则在均值点进行泰勒级数展开并取其一次式来计算可靠度，这种方法称为中心点法。显然，在上述计算过程中并没有考虑到变量的实际分布情况，而只考虑了他们的均值和方差，或者说，是将各个随机变量假定为正态分布或者对数正态分布变量进行计算的。理论与实践均证明，对于非线性极限状态方程，中心点法计算误差较大。后来Hasofer-Lind和Rachwitz-Fiessler等人提出改进的中心点法，被称为验算点法。这一方法是考虑了基本变量分布类型对可靠度的影响，将非线性功能函数的线性化点选为设计验算点通过当量正态化进行可靠度计算。（2）概率矩点估计法概率矩点估计法是由Emilio-Rosenblueth于1975年提出的，1981年他又对这一方法进行了完善和理论化，通常又称为Rosenblueth法。从80年代初开始，将这一方法引入到了岩土及边坡工程的可靠度分析中。它主要是根据输入随机变量的前三阶矩（均值、方差、偏态系数）来近似地描述极限状态函数的概率矩，不必预先知道输入随机变量的精确分布，应用方便。就许多岩土工程问题来说，也具有足够的精度，是一种比较实用的工程方法，近几年有许多人对其精度和实用性进行了深入讨论，并且都给予了充分肯定。事实上概率矩点估计法与蒙特卡罗法一样，他们的最大优点是不用改变边坡的极限状态函数，只需在基本输入随机变量的空间中选取一些模拟数值计算点，对极限状态函数进行简单的数值运算就可以求出破坏概率或概率矩，因此是非常方便的。概率矩点估计法的取值点大大少于蒙特卡罗方法，而且不需事先给出随机变量的分布类型，应用前景广泛。一般说来，对于个随机变量的情况，需选取计算点，当较小时（），比较经济；当较大时，选取的计算点数就会大量增加，这时可以选用蒙特卡罗方法进行模拟计算。（3）蒙特卡罗方法蒙特卡罗方法（Monte-Carlo）方法又称随机模拟方法或统计试验法，是一种依据概率统计学理论利用计算机研究随机变量的数值计算方法，它的应用非常广泛，在目前是相对精确的方法。蒙特卡罗模拟法通过随机模拟与统计试验来求解边坡工程的可靠度，多用于理论方面，以便检验其它方法的精度。该方法的关键在于：随机抽样数和随机抽样方法的确定。由概率论知，采用频率来估算概率的基本前提是随机抽样数必须足够大，否则达不到精度要求，而抽样数大了就增加了工作量。所以为克服这个矛盾就发展了很多的方法：重要抽样法、对偶抽样法、分层抽样法、条件期望值法、公共随机数法、图解渐近法等等。以上简单介绍了三种计算边坡工程可靠度或失效概率的方法。一般来说，中心点法适用于基本变量为正态分布或对数正态分布，极限状态为线性且可靠度要求不高的情况。验算点法只有在各随机变量间是独立的正态随机变量且功能函数为线性是计算结果比较精确。目前这几种方法来说，还很难说哪一种方法更精确一些，从原理上讲，蒙特卡罗方法更明确一些，而且当模拟次数增加后，其精度是可以肯定的，因此通常用这一方法来检验其它方法的精度和适用性。Hasofer-Lind和Rackwitz的可靠指标法是近20年来结构可靠性计算方面的一个主要成果，它们将可靠指标计算的解析过程转化成了几何求解过程，并且都采用了泰勒级数在设计点一次线性展开，因此它们也是近似求解方法，在工程中有较广泛的应用。它们的不足之处是必须对极限状态函数进行求导、变换等运算，当函数形式复杂时，会感到繁琐。本文选用蒙特卡罗方法求解可靠度，从蒙特卡罗方法的思路可以看出，该方法回避了可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多就可得到一个比较精确的失效概率和可靠度指标。特别在岩土体分析中，变异系数往往较大，在这种情况下，采用蒙特卡罗方法计算的可靠度更为精确，并且由于思路简单易于编制程序。事实上，选用哪一种方法进行边披的可靠度分析往往是依个人的习惯和经历．不能笼统地否认或宣扬哪一种方法，对具体问题来说，每一种方法的计算结果都可能会有一些偏差，但这种偏差往往是从事边坡工程研究人员不能克服的，而边坡破坏机理和破坏模式等方面的不完善或不合理引起的边坡可靠性差异远大于方法本身的误差，因此更应关注的是后者。2.2用瑞典法建立极限状态方程瑞典条分法是基于极限平衡原理先假定若干可能的剪切面——滑裂面，然后将滑裂面以上土体分成若干垂直土条，对作用于各土条上的力进行力与力矩的平衡分析，求出在极限平衡状态下土体稳定的安全系数，并通过一定数量的试算，找出最危险滑裂面位置及相应的(最低的)安全系数。由于土坡的稳定问题是一个高次超静定问题，必须要作出各种简化假定以减少未知量。瑞典条分法的基本假定：假定土坡稳定属于平面稳定问题；假定滑裂面为圆柱面，弧面上的滑动体为不变形的刚体；认为条块间的作用力对边坡的整体稳定性影响不大，可以忽略，或者说，假定条块两侧的作用力大小相等，方向相反且作用于同一直线上。图2-1条块受力分析如图2-1所示，把滑动土体分成若干土条，土条的两个侧面存在着条块间的作用力，对于第i条土块，受的力有重力条块侧面ac和bd作用有法向力、，切向力、。滑弧段cd的长度为，其上作用着法向力和切向力，中包括粘聚阻力和摩擦阻力，由于条块的宽度不大，和可以看成作用于弧段cd的中点。在这些力中，由于不考虑条间力作用，根据径向平衡条件，有：（2-1）设安全系数为，根据滑弧面上的极限平衡条件，有：（2-2）式中：——条块i在滑动面上的抗剪强度。整个滑动土体对圆心O取力矩平衡得：（2-3）将式（2-1）代入（2-2）再将式（2-2）代入（2-3）得如下计算：（2-4）式中：——第i条滑裂面处土体粘聚力，——第i条滑弧弧长，——第i条自重，，其中和为第i条宽度和高度，——第i条底面中点的法线与竖直线交角，——第i条内摩擦角。因此得到边坡的极限状态方程：（2-5）2.3蒙特卡罗法的基本原理[13]当已知基本变量的概率分布时，可利用适当的随机数发生器，产生符合状态变量的概率分布的一组随机数以之代入状态函数计算状态函数的一个随机数，并看他是否小于零。以同样的方式产生个状态函数的随机数据。若个状态函数的随机数中有个小于零，则当足够大时，由大数定律可知系统的失效概率为：（2-6）可靠度指标为：（2-7）2.4产生给定随机变量分布类型的随机数用Mont-Carlo法分析边坡工程可靠度问题，关键是产生已知分布变量的随机数。为了快速、高效地产生给定随机变量的随机数，通常分两步进行，即首先在区间（0，1）内产生均匀分布随机数，然后再变换成给定变量的随机数。（1） 产生（0，1）间的随机数本文用乘同余法获取随机数，算式如下：（2-8）式中，ａ，ｃ和ｍ为正整数。式（2-8）表示以ｍ为模数的同余式，为以ｍ除后得到的余数。具体计算时，可引入参数，即令　　　　　　　　　　　　　　　　（2-9）式中，为取整符号，式（2-8）成为　　　　　　　　　　　　　　　　（2-10）　将除以模数，可得标准化的随机数，即（2-11）以上是产生（0，1）间的随机数。为了判别所得随机数的优劣，一般还应对伪随机数进行统计检验，主要是检验其均匀性和独立性。相关内容可见参考文献[14]（2）产生给定随机变量分布类型的随机数下面介绍产生给定随机变量服从正态分布、对数正态分布的随机数的方法。1）正态分布设随机数和是（0，1）区间的两个均匀随机数，则可通过下列变换得到服从标准正态分布的两个随机数和（2-12）如果随机变量为一般正态分布，则其随机数和由下式计算：（2-13）这里成对产生的随机数不仅互相独立，而且服从一般正态分布。2）对数正态分布对数正态分布变量随机数产生的方法是先将均匀随机数变换为正态分布随机数，然后再转换为对数正态分布随机数。设为对数正态分布变量，其均值为，标准差为，变异系数为。因为为正态分布，所以，得其均值和标准差分别为变量的随机数产生后，则可得到对数正态分布变量的随机数为2.5相关变量条件下蒙特卡罗模拟法对于一组正态变量，设其均值为，标准差为，相关系数矩阵为（2-14）其中，为变量与的相关系数。再设与相关的标准正态随机变量为，于是有（2-15）式中，由于变量与的协方差为所以标准正态变量的协方差矩阵为（2-16）根据线性代数的理论，通过正交变换可将矩阵化为对角矩阵，其对角元素即为的特征值。令（2-17）式中，为正交矩阵，其列向量为的特征向量，向量即为一组不相关的正态变量，他的协方差矩阵为（2-18）并且有（2-19）进而得到一组不相关的标准正态变量（2-20）而（2-21）其中（2-22）以上进行了理论分析，具体计算过程总结如下：（1） 建立极限状态方程，（2） 将相关的随机变量进行独立化处理，参见上面内容，（3） 通过变换，产生，的给定具体分布类型变量的随机数，（4） 把所得变量的随机数代入已经建立的功能函数式，（5） 累计记录的次数及抽样总次数，（6） 当计算抽样次数足够时，抽样结束，失效概率为。三、工程实例现有一高7.62m、边坡角45°、边坡体容重的简单均质土体边坡，为了简化问题，忽略土的重度变异性对可靠度的影响，土体的破坏符合Mohr-Coulomb破坏准则，不考虑地下水的作用和地震作用。设只有粘聚力和内摩擦系数是随机变量，粘聚力服从正态分布，内摩擦系数也服从正态分布，粘聚力的均值，标准差，内摩擦系数的均值，标准差，随机变量和是相关的，不考虑参数的自相关性。设滑面为圆弧形，过边坡脚，滑面半径，最危险滑裂面及滑弧圆心及土坡模型见图3-1。对一般性工程，取N=5000～10000次，已完全可以满足工程评价精度要求[9]。为满足研究需要，在下面的计算中，均取N=10000次。图3-1边坡断面1） 建立边坡的极限状态函数。根据边坡的破坏模式，采用瑞典条分法，将滑体分成10个小分条，第i个小分条的重量、滑面倾角为，不计分条间的作用力。用安全余量定义边坡的极限状态函数：（3-17）式中：——第i条滑弧弧长，——第i条自重，，其中和为第i条宽度和高度，边坡的破坏概率：2）将互相关的随机变量和进行独立化处理。3）产生符合随机变量和分布的随机数。由乘同余法产生[0，1]区间内均匀分布的随机数，然后用坐标变换法产生给定均值，标准差，，标准差的正态分布随机数。4）计算极限状态函数值。将随机变量和的组随机数代入极限状态函数式（3.1.1）中，得出的个随机值。5）累计记录的次数及抽样总次数。6）当计算抽样次数足够时，抽样结束，失效概率为式中：——个的随机值中小于等于零的个数。四、结语与展望本文采用蒙特卡罗模拟法（Monte-Carlo）求解了边坡工程的破坏概率和可靠指标，考虑了粘聚力和摩擦系数的互相关性。实际上，土体材料参数不仅具有随机性，也有相关性。更加合理、确切地分析描述边坡工程质量状态，在边坡可靠度分析中应该基于岩土性质及边坡破坏的客观事实，岩土性质不是简单的纯随机变量，他们在空间分布构成了具有空间变异性的随机场，除具有经典概率统计学中随机变量的有关性质外，还具有某种空间相关性或称连续性。对其研究表述的基本方法，理论应是随机场理论，而不是简单的经典概率统计理论。边坡工程随机计算中应该考虑到岩土性质的空间变异性，并建立相应的空间变异模型，推导出边坡随机计算中所需的一些有关参数。由于成因、剥离、风化以及其它地质作用，岩土性质的空间变异总是存在的，研究某些岩体性质的空间变异必然要联系到其节理的空间随机模拟等；另外，将大部分软岩和土体视为空间随机介质，可用随机场理论进行研究。随机场理论已被进行了大量研究。近年来，随机有限元法已成为分析边坡工程可靠性的有效的方法之一。随机有限元法是随机分析理论与确定性有限元法相结合的产物。随着科学技术的发展和对自然界事物认识的深入，对工程中不确定性因素的描述经历了从随机变量、随机过程到随机场这一过程。根据对工程进行随机分析的方法与手段，随机有限元法可以分为很多种。总的来说，分为两大类型：一类是统计的方法，就是通过大量的随机抽样，对工程结构反复进行有限元计算，将得到的结果做统计分析，得到该工程结构的失效概率和可靠度指标，这种方法称为蒙特卡罗（Monte-Carlo）随机有限元法。Monte-Carlo随机有限元法需要进行大量的模拟计算，因此工作量很大，使这一方法的应用受到一定的限制。另一类是分析的方法，就是以数学、力学分析作为工具，找出结构系统（确定的或随机的）的响应与输入信号（确定的或随机的）之间的关系，并据此得到结构内力、应力或位移的统计规律，得到结构的失效概率和可靠度指标。这一类随机有限元方法比较多，按照随机分析的目的与结果不同，又可分为两种：一种是分析结构响应的统计特性及其分布规律，如摄动随机有限元法（PSFEM）、纽曼（Neumann）随机有限元法等；另一种是直接分析结构的可靠度或失效概率和可靠度指标，如验算点随机有限元法等。边坡工程可靠度分析的方向是正确的，它能真实地考虑到客观和人为的一些不确定性因素，并已逐步运用于工程实践，从事研究和接受这一方法思想的人员也在增加，因此具有生命力。可靠度分析方法的应用会越来越广泛。未来边坡稳定性研究的发展趋势体现在以下几个方面：(l)两种或多种稳定性分析方法联合并用，对边坡工程进行综合分析验证，力求得出一个更加客观、可靠、合理的评价结果。(2)由于可靠度分析方法在处理不确定性问题方面具有非常大的优势，因而可以在边坡稳定性分析上发挥重要作用。可靠度分析方法的不断完善，将是边坡稳定性分析研究的一大趋势。如失稳概率分析，模糊数学综合评价方法等可靠度分析方法研究和应用。(3)采用反分析方法来研究边坡系统的失稳。可以转化思维，从边坡系统的变形分析出发，寻求边坡的变形失稳判据并进行稳定性评价。如位移反分析法就是先以实测的位移值为依据反演求出初始的应力与参数，再反过来应用于工程实践。(4)随着数值分析方法的不断发展，不同数值方法的相互藕合成为一大发展趋势。如有限元，离散元与块体元等的相互藕合，数值解和解析解的结合，这些方法的藕合能充分发挥各自的优点，解决更复杂的边坡问题。(5)边坡失稳的发展是一个渐变的过程，渐变到一定程度就会破坏，这是一个突变过程，因而用突变理论来评价边坡工程的失稳同样具有广阔的前景。(6)三维边坡稳定性评价一直是边坡研究人员探索的问题。目前，在破坏机制评价理论等方面都没有成熟的方法。因此，三维边坡稳定性评价的研究将是边坡工程研究的重点。(7)如何建立数字边坡用于指导设计和安全评价将是未来研究内容之一。它将对准确的边坡稳定性评价提供基础，但在数据采集方面有一定难度。通过对边坡工程可靠度分析方法的阐述与未来趋势的分析，正确认识了各种边坡稳定性分析方法的优缺点及其适用范围，未来边坡稳定性的研究趋势也将向着多元化方向发展。但由于边坡工程多建立于复杂地质体中，有极其复杂多变的特性，同时又有较强的隐蔽性。因此，十分精确的对边坡工程的稳定程度作出评价是不太现实的。但随着人们对边坡工程研究的不断深人，边坡稳定性评价方法也会随之不断地完善和发展。参考资料[1]王家臣.边坡工程随机分析原理[M].北京.煤炭工业出版社,1996.[2]林忠民.工程结构可靠性设计与估计[M].北京.人民交通出版社,1990.[3]吴世伟.结构可靠度分析[M].北京.人民交通出版社,1990.[4]姜弘道,陈和群.有限单元法的程序设计[M].北京.水利电力出版社,1989.[5]洪昌华,龚晓南.相关情况下的Hasofer-Lind可靠度指标的求解[J].岩土力学,2000.[6]傅旭东,赵善锐.用蒙特卡洛方法(Mont-CarloMethod)计算岩土工程的可靠度指标[J].西安交通大学报,1996,31(2).56-60.[7]陈祖煜.土质边坡稳定分析——原理·方法·程序[M].北京.中国水利水电出版社,2003.[8]刘宁.可靠度随机有限元法及其工程应用[M].北京.中国水利水电出版社,2001.[9]祝玉学.边坡可靠性分析.冶金工业出版社，1993.[10]黄昌乾，丁恩保.边坡工程常用稳定分析方法[J].水电站设计，1999，15(l).53-58.[11]王勤成，邵敏.有限单元法基本原理和数值方法[M].北京.清华大学出版社，2002.[12]邱天.边坡稳定的不确定性研究现状与展望[J].山西建筑，2006，32(4).114-115.[13]武清玺.结构可靠性分析及随机有限元法，机械工业出版社，2005[14]徐钟济.蒙特卡罗方法.上海.上海科学技术出版社，1985附录PROGRAMrand-cov !****************产生随机数********************* PARAMETER(N=10,PI=3.1415926,N1=3) REALarray(N),ZZ(N),AZ(N1,N1),SQ(N1,N1)REAL::EO=1E6REAL::X0=1.0,X1,u1,JZ,BZC INTEGER(4)::a=2**7+1;c=1;m=2**20 PRINT*,'请输入均值:' READ*,JZ PRINT*,'请输入标准差:' READ*,BZCOPEN(1,FILE='example1.in')OPEN(2,FILE='example1.out')J=1 DOI=1,N K0=INT((a*X0+c)/m) WRITE(2,*)K0 X1=a*X0+c-m*K0 u1=X1/m X0=X1 array(J)=u1 J=J+1 ENDDO WRITE(2,1000)(array(J),J=1,N)1000FORMAT(5X,'array(N)='/,(10X,4F13.6)/)DOJ=1,N,2 IF((JZ.EQ.0.0).AND.(BZC.EQ.1.0))THEN Z0=SQRT(-2*LOG(array(J)))*cos(2*PI*array(J+1)) Z1=SQRT(-2*LOG(array(J)))*sin(2*PI*array(J+1)) ELSE Z0=SQRT(-2*LOG(array(J)))*cos(2*PI*array(J+1))*BZC+JZ Z1=SQRT(-2*LOG(array(J)))*sin(2*PI*array(J+1))*BZC+JZ ENDIF ZZ(J)=Z0 ZZ(J+1)=Z1 ENDDO WRITE(2,2000)(ZZ(J),J=1,N)2000FORMAT(5X,'ZZ(N)='/,(10X,4F13.6)/)!******************相关性变量************************** READ(1,*)((AZ(I,J),J=1,N1),I=1,N1)WRITE(2,6)((AZ(I,J),J=1,N1),I=1,N1) 6FORMAT(5X,'AZ(I,J)='/(5X,3F13.5))DO10I=1,N1 DO10J=1,N1 SQ(I,J)=0.010CONTINUEDO15I=1,N1 SQ(I,I)=1.015 CONTINUE G=0.0 DO40I=2,N1 I1=I-1 DO40J=1,I1 G=G+2.0*AZ(I,J)**240CONTINUEWRITE(2,42)G42FORMAT(10X,'G=',F13.5)S1=SQRT(G) FN=FLOAT(N1) S2=EO*S1/FN S3=S1WRITE(2,45)S245FORMAT(10X,'S2=',F13.5)L=050S3=S3/FNWRITE(2,55)S355FORMAT(10X,'S3=',F13.5)60DO130IQ=2,N1IQ1=IQ-1 DO130IP=1,IQ1 IF(ABS(AZ(IP,IQ)).LT.S3)GOTO130 L=1 V1=AZ(IP,IP) V2=AZ(IP,IQ) V3=AZ(IQ,IQ) U=0.5*(V1-V3)WRITE(2,65)V1,V2,V3,U65FORMAT(10X,'V1*V2*V3*U=',4F13.5)IF(U)70,80,9070B=-1.0GOTO10090B=1.0100G=-B*V2/SQRT(V2*V2+U*U)GOTO10580IF(V2.GT.S3)GOTO85G=1.0 GOTO10585G=-1.0105ST=G/SQRT(2.0*(1.0+SQRT(1.0-G*G)))CT=SQRT(1.0-ST*ST)WRITE(2,107)ST,CT107FORMAT(10X,'ST*CT=',2F13.5)DO110I=1,N1 G=AZ(I,IP)*CT-AZ(I,IQ)*ST AZ(I,IQ)=AZ(I,IP)*ST+AZ(I,IQ)*CT AZ(I,IP)=G G=SQ(I,IP)*CT-SQ(I,IQ)*ST SQ(I,IQ)=SQ(I,IP)*ST+SQ(I,IQ)*CT SQ(I,IP)=G110CONTINUEDO120I=1,N1 AZ(IP,I)=AZ(I,IP) AZ(IQ,I)=AZ(I,IQ)120CONTINUEG=2.0*V2*ST*CT AZ(IP,IP)=V1*CT*CT+V3*ST*ST-G AZ(IQ,IQ)=V1*ST*ST+V3*CT*CT+G AZ(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST) AZ(IQ,IP)=AZ(IP,IQ)130CONTINUEIF(L-1)150,140,150140L=0GOTO60150WRITE(2,145)((AZ(I,J),J=1,N1),I=1,N1)WRITE(2,142)((SQ(I,J),J=1,N1),I=1,N1)142FORMAT(10X,'SQ(I,J)='/(10X,3F13.5))145FORMAT(10X,'AZF(I,J)='/(10X,3F13.5))IF(S3.GT.S2)GOTO50 END