考虑土体参数空间变异性的边坡地震滑移概率分析方法

1.本发明属于地震作用下边坡稳定性评价的技术领域，具体涉及一种考虑土体参数空间变异性的边坡地震滑移概率分析方法。


背景技术：

2.地震引发的滑坡灾害对人身安全、经济建设等可能造成严重的影响。因此，地震作用下的边坡稳定性问题受到了学术界和工程界的普遍关注。地震引起的永久位移是衡量边坡稳定性的重要参数之一，对该参数的准确预测对边坡工程抗震设计、加固等具有重要意义。另一方面，岩土体是一种经过各种地质、环境和物理化学过程作用形成的天然材料，其性质的区域性和个性较强，存在某种程度的天然变异性。现有的边坡地震滑移分析一般仅考虑地震荷载的不确定性，而对土体参数进行确定性取值，忽略了土体天然变异性可能导致的物理力学参数不确定性。近年来有研究提出逻辑树方法，将土体参数取值用离散型随机变量来描述，从而一定程度量化该不确定性。然而，该方法由于其固有的近似特性，无法表征土体参数的空间变异性并量化其对地震滑移预测的影响。因此，如何在边坡地震滑移分析中考虑土体参数的空间变异性是当前亟待解决的问题。统计学中的随机场理论为解决该问题提供了一种可行的工具。
3.考虑土体参数空间变异性往往涉及对参数的多次随机抽样(比如，蒙特卡洛模拟)以及相应的边坡滑移计算，若采用复杂的有限元计算需要耗费大量的时间和计算机资源，且不易于工程师掌握和使用。有必要发展一种能较好权衡精确程度和计算效率的地震边坡滑移估计模型。特别地，需要在计算中有效地表征土体的非线性和非均匀性特点。另外，有必要在边坡地震滑移概率分析中同时考虑地震荷载的不确定性和周期性，输出与重现期相关的滑移阈值超越概率，进而为地震作用下的边坡稳定性提供一个有效的评价指标。


技术实现要素：

4.本发明的目的在于针对现有技术的不足之处，提供一种考虑土体参数空间变异性的边坡地震滑移概率分析方法，该方法计算过程较为简便、计算效率和精确程度均较高，同时考虑了土体参数和地震动的不确定性，可以预测地震作用下边坡在某个设计基准期的滑移概率。
5.为解决上述技术问题，本发明采用如下技术方案：
6.一种考虑土体参数空间变异性的边坡地震滑移概率分析方法，包括以下步骤：
7.步骤1、确定土体参数空间分布的统计特征，构造参数随机场变量的联合概率密度函数f(x)；
8.步骤2、根据联合概率密度函数f(x)，采用蒙特卡洛模拟方法产生土体参数随机场样本x(i)，i＝1，2，...，n，i表示随机场样本编号；
9.步骤3、构建边坡稳定性分析模型，计算土体参数随机场样本x(i)对应的边坡临界加速度系数i＝1，2，...，n，并记录临界滑动体的位置；
10.步骤4、根据具体边坡抗震设计条件，挑选一组地震动加速度时间历程k＝1，2，...，m，k表示地震动编号；
11.步骤5、构建边坡临界滑动体多自由度永久位移分析模型，采用等效线性算法计算随机场样本x(i)和地震动加速度时程对应的滑动体永久位移d
(i，k)
，i＝1，2，...，n，k＝1，2，...，m；
12.步骤6、对所得的永久位移样本进行统计分析，计算在某个设计基准期td下滑动体永久位移对设计位移阈值d的超越概率pd。
13.进一步地，采用随机场理论描述土体参数的自相关结构以及采用nataf理论描述土体参数的互相关结构，根据工程试验数据或文献资料合理选取所需的参数统计特征。
14.进一步地，统计特征包括均值μ、标准差σ、水平波动范围δ
x
和竖直波动范围δy、不同参数间的皮尔逊相关系数ρ
cross，x
。
15.进一步地，土体参数包括土体重度γ、不排水抗剪强度、动力土性参数。
16.进一步地，步骤3包括以下子步骤：
17.3.1根据步骤1的随机场网格单元对边坡材料区域进行划分，并将步骤2的随机场样本x(i)赋予相应的土体参数单元；
18.3.2采用极限平衡法迭代计算水平地震力系数kh作用下边坡最小安全系数为1的临界加速度系数并记录其对应临界滑动体的位置。
19.进一步地，根据具体边坡抗震标准、所在地区和场地条件拟定地震设计反应谱，挑选一组反应谱与该设计反应谱相匹配的地震动加速度时间历程。
20.进一步地，步骤5中，每个位移数据d
(i，k)
的计算包括如下子步骤：
21.5.1把随机场样本x(i)对应的临界滑动体理想化为一个多自由度体系，构建质量矩阵m、刚度矩阵k、阻尼矩阵c；
22.5.2根据多自由度体系的运动状态，交替建立并求解在加速度时程作用下体系黏合阶段和相对滑移阶段的运动方程，计算临界滑动体各层的等效剪应变γ
eff，l
与底部的滑动位移d；
23.5.3根据各层的等效剪应变γ
eff，l
更新土体剪切模量g
l
和阻尼比λ
l
；
24.5.4重复子步骤5.1-5.3直至各层当前等效剪应变与前一次等效剪应变间的差值均小于一个预设阈值，将此时的d作为随机场样本x(i)和地震动加速度时程对应的临界滑动体永久位移d
(i，k)
。
25.进一步地，根据泊松分布假设计算地震动的发生概率，并将其与位移超越阈值d的样本比例相乘估计某个设计基准期td下滑动体永久位移超越概率pd。
26.与现有技术相比，本发明的有益效果为：
27.(1)本发明结合结构动力学原理和土动力学原理，采用等效线性多自由度体系代表边坡临界滑动体，实现了对非均匀土质边坡地震滑移距离的高效估计；此外，本发明利用泊松分布模型考虑地震动的不确定性，可预测边坡滑移距离对工程允许阈值在任意设计基准期的超越概率；至此本发明克服了传统边坡地震滑移概率分析方法对土体参数空间变异性考虑的不足，提供了一种能考虑土体参数空间变异性、计算过程较为简便、计算效率和精确度较高的边坡地震滑移概率分析方法，可应用于边坡工程抗震设计实践；
28.(2)本发明的方法基于随机场理论和nataf理论表征多土体参数的空间变异性，基于结构动力学的多自由度体系计算非均匀土体边坡的地震滑移量，并采用等效线性算法考虑土体的非线性效应；
29.(3)本发明在边坡地震滑移概率分析中同时考虑了地震荷载的不确定性和周期性，输出与重现期相关的滑移阈值超越概率，进而为地震作用下的边坡稳定性提供一个有效的评价指标，可以预测地震作用下边坡在某个设计基准期的滑移概率。
附图说明
30.图1为本发明的具体流程示意图；
31.图2为本发明中随机场变量间的相关性示意图；
32.图3为本发明中边坡地震滑移计算方法流程图；
33.图4为本发明中边坡滑动体的等效多自由度体系示意图；
34.图5为本发明实施例中边坡的剖面图；
35.图6为本发明实施例中边坡临界加速度系数样本的概率分布图；
36.图7为本发明实施例中所选地震动的反应谱；
37.图8为本发明实施例中位移超越概率计算结果图。
具体实施方式
38.下面将结合本发明实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
39.需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
40.如图1所示，本发明提供了一种考虑土体空间变异性的边坡地震滑移概率分析方法，包括以下步骤：
41.步骤1，确定土体参数空间分布的统计特征，构造参数随机场变量的联合概率密度函数f(x)；
42.本发明中考虑的土体参数包括土体重度γ、不排水抗剪强度su、动力土性参数(小应变剪切模量g
max
、剪切模量比rg＝g/g
max
、阻尼比λ)。采用随机场理论刻画土体参数的空间变异性，以一定数量的相关随机变量x代表不同位置的不确定性土体参数，根据适当的自相关函数、互相关系数来建立这些变量的联合概率密度函数。
43.具体地，在步骤1中，首先需要确定土体参数的统计特征，如均值μ(反映某种参数的基本趋势)、标准差σ(反映某种参数的离散程度)、水平波动范围δ
x
和竖直波动范围δy(反映某种参数的空间不均匀程度)、不同参数间的皮尔逊相关系数ρ
cross，x
。当工程试验数据充足时，可直接基于样本估计；反之，可根据经验和文献资料，对参数统计特征进行合理选取。
44.现以两个不确定性土体参数x1和x2为例，介绍随机场变量联合概率密度函数的构造方法。将随机场剖分为n个网格单元，参数x1和x2的随机场变量可分别表示为(x
1，1
，x
1，2
，...，x
1，n
)
t
和(x
2，1
，x
2，2
，...，x
2，n
)
t
。采用nataf理论考虑x1与x2之间的互相关性，所有随
机场变量x＝(x
1，1
，x
1，2
，...，x
1，n
，x
2，1
，x
2，2
，...，x
2，n
)
t
的联合概率密度函数f(x)可推导为：
[0045][0046]
式中，f1(x
1，j
)和f2(x
2，j
)分别为x
1，j
＝x
1，j
和x
2，j
＝x
2，j
的边缘概率密度函数(其中j＝1，2，...，n)；z＝(z
1，1
，z
1，2
，...，z
1，n
，z
2，1
，z
2，2
，...，z
2，n
)
t
为服从2n元标准正态分布的随机向量，其相关系数矩阵为r，变量间的相关性如图2所示；z＝(z
1，1
，z
1，2
，...，z
1，n
，z
2，1
，z
2，2
，...，z
2，n
)
t
为z的具体取值；φ1(z
1，j
)和φ2(z
2，j
)分别为z
1，j
＝z
1，j
和z
2，j
＝z
2，j
的边缘概率密度函数(即一元标准正态分布的概率密度函数)；φ
2n
(z；r)为2n元标准正态分布的联合概率密度函数，其表达式为：
[0047][0048]
式中，|
·
|为矩阵行列式；维度为2n
×
2n的相关系数矩阵r用分块矩阵表达为：
[0049][0050]
式中，ρ
cross，z
为正态分布空间的相关系数，需要对原始空间的土体参数互相关系数ρ
cross，x
进行转换而得，可根据相关文献(kiureghian a d，liu p l.structural reliability under incomplete probability information[j].journal of engineering mechanics，1986，112：86-114)提供的方法求解；一般情况下，ρ
cross，z
和ρ
cross，x
相差不大，可近似取ρ
cross，z
＝ρ
cross，x
；a为空间自相关系数矩阵：
[0051][0052]
式中，矩阵的每个元素可通过自相关函数计算得到，本具体实施中采用较常用的单指数型自相关函数：
[0053][0054]
式中，τ
x
(r，s)和τy(r，s)分别为第r个与第s个随机场单元中点的水平和竖直间距；至此，式(1)中所需元素可由式(2)～(5)获得，可构建随机场变量的联合概率密度函数f(x)。
[0055]
步骤2，根据步骤1构建的联合概率密度函数f(x)，采用蒙特卡洛模拟方法产生服从f(x)的土体参数随机场样本x(i)，i＝1，2，...，n，i表示随机场样本编号；
[0056]
进一步地，步骤2中，由于不确定性土体参数可能服从除正态或对数正态分布以外的其它边缘分布类型，f(x)难以用理论形式表达，且不易直接产生相应随机数。可采用nataf逆变换产生随机场样本：首先产生相关系数矩阵为r的2n元标准正态随机向量的一组样本z＝(z
1，1
，z
1，2
，...，z
1，n
，z
2，1
，z
2，2
，...，z
2，n
)
t
，可借助编程语言(如matlab中的mvnmd函
数)轻易实现；然后根据等概率变换将样本z转换为随机场样本向量x(i)＝(x
1，1
，x
1，2
，...，x
1，n
，x
2，1
，x
2，2
，...，x
2，n
)
t
，转换公式如下：
[0057][0058]
式中，φ(
·
)为标准正态变量的累计分布函数；φ-1
(
·
)为标准正态变量累计分布函数φ(
·
)的反函数；和分别为变量x
1，j
和x
2，j
累计分布函数的反函数，j＝1，2，...，n，j表示随机场单元编号。
[0059]
步骤3，构建边坡稳定性分析模型，计算土体参数随机场样本x(i)对应的边坡临界加速度系数i＝1，2，...，n，并记录临界滑动体的位置；
[0060]
具体地，步骤3还包括以下子步骤：
[0061]
3.1根据步骤1的随机场网格单元对边坡材料区域进行划分，并将步骤2的随机场样本x(i)赋予相应的土体参数单元；
[0062]
目前常用的边坡稳定性分析法有极限平衡法、有限元法和有限差分法，可通过商业软件geostudio(极限平衡法)、abaqus(有限元法)、flac(有限差分法)等实现，本发明适用任何一种边坡稳定性分析法。在具体边坡工程设计中，进行合理选取，其中以简化毕肖普法、摩根斯坦-普赖斯法等为代表的极限平衡法能实现较为精确的边坡安全系数估计，同时要求相对较少的计算成本。因此，本具体实施中采用简化毕肖普法进行安全系数的计算。由于考虑土体参数非均匀性后无法事先确定临界滑动体的位置，设置多条潜在滑动面以基本涵盖整个边坡区域；
[0063]
3.2计算水平地震力系数kh作用下的边坡最小安全系数fs，当fs＜1时减小kh值，当fs＞1时增大kh值，直至到达临界状态(即fs＝1)，将此时的地震力系数记为临界加速度系数并记录其对应临界滑动体的位置；
[0064]
循环n次执行子步骤3.1和子步骤3.2，即i＝1，2，...，n。
[0065]
步骤4，根据具体边坡抗震标准(地震动重现期ta)、所在地区和场地条件确定地震设计反应谱(目标谱)，可借助中国地震动参数区划图(gb18306-2015)实现，进而挑选一组反应谱与该目标谱相匹配的地震动加速度时间历程k＝1，2，...，m，k表示地震动编号；
[0066]
步骤5，构建边坡临界滑动体多自由度永久位移计算模型，估计随机场样本x(i)和地震动加速度时程对应的滑动体永久位移d
(i，k)
，i＝1，2，...，n，k＝1，2，...，m。采用等效线性算法考虑土体剪切模量比rg＝g/g
max
和阻尼比λ随剪应变的非线性变化；
[0067]
其中，步骤5的具体流程见图3，主要包括以下子步骤：
[0068]
5.1将随机场样本x(i)对应的临界滑动体理想化为一个如图4所示的多自由度体系，构建质量矩阵m、刚度矩阵k、阻尼矩阵c；
[0069]
将该滑动体剖分为l个水平层，各层质量在层间位置处按结点集中(取结点上下两分层质量之和的一半)。具体而言，第l个节点(位于第l-1层与l层之间)的质量可表示为：
[0070][0071]
式中，γ
l
为第l层土体的重度；g为重力加速度；h
l
为第l层的厚度。将第层余下的一半质量赋给滑动面处的节点，即结点间由无质量的剪切弹簧和粘滞阻尼器相连，其中弹簧刚度k
l
和阻尼系数g
l
分别表达为：
[0072][0073][0074]
式中，ω为多自由度体系第一阶振型的自振频率，可参考相关文献(chopra a.dynamics of structures[m].limited：pearson education，2007)求解含有体系质量矩阵m和刚度矩阵k的特征方程而得；g
l
为第l层土体的剪切模量；λ
l
为第l层土体的阻尼比。
[0075]
当第一次执行步骤5.1时，剪切模量g
l
设定为小应变剪切模量g
max
，即随机场样本x(i)中第l层位置g
max
单元的平均值，阻尼比λ
l
则在合理范围内(如0.01～0.05)假设。当在迭代过程中执行步骤5.1时，根据该层有效剪应变γ
eff，l
分别查询剪切模量比rg和阻尼比λ随剪应变变化的经验曲线估计g
l
和λ
l
，其中前者需要运算g
l
＝g
max，l
rg，后者直接获得。进一步地，可将多自由度体系的质量矩阵m、刚度矩阵k、阻尼矩阵c分别表达为：
[0076][0077][0078][0079]
5.2根据多自由度体系的运动状态，交替建立并求解加速度时程作用下体系黏合阶段和相对滑移阶段的运动方程，计算临界滑动体各层的等效剪应变γ
eff，l
与底部的滑
动位移d。
[0080]
在地震加速度的作用下，滑体并不是始终相对于底部土体而滑动，取决于该时步滑动体承受的滑动力是否超过滑动体的抗滑力。当滑动力小于抗滑力时，滑动体与底部土体处于黏合状态，不产生滑移量；反之，滑动体与底部土体处于相对滑移状态，存在滑移量的累积。
[0081]
(1)黏合阶段
[0082]
当滑动体与底部土体处于黏合状态时，滑移体系内各质点的运动方程可表达为：
[0083][0084]
式中，为维度为l(i)×
1的滑面以上节点相对加速度向量；为维度为l(i)×
1的滑面以上节点相对速度向量；u为维度为l(i)×
1的滑面以上节点相对位移向量；1为维度为l(i)×
1的单位向量。在每个时间步，求解式(13)得到节点相对加速度向量再判断底部滑动力是否超过最大静摩擦力，即：
[0085][0086]
式中，m
all
为所有节点质量之和为步骤3求得的临界加速度系数，可在此视作摩擦系数；g为重力加速度；t为矩阵转置。当式(14)不再满足时，由黏合阶段进入相对滑动阶段。
[0087]
(2)相对滑移阶段
[0088]
当滑动体与底部土体处于相对滑移状态时，滑移体系内各质点的运动方程可表达为：
[0089][0090]
式中，为此刻滑动体相对于底部土体的滑动加速度，可以通过滑动体底部平衡方程而表达为：
[0091][0092]
将式(16)代入式(15)可消掉从而得到：
[0093][0094]
因此，在相对滑移阶段，先通过式(17)求解节点相对加速度向量再代入式(16)求解底部的滑动加速度同时，对进行一次数值积分得底部的滑动速度当减小至0时，判定当前的相对滑移阶段已结束，并重新进入黏合阶段。循环上述过程直到地震加速度的最终时刻，获得底部滑动速度的时间历程对整个时间历程再进行一次数值积分并累加便可得底部滑动位移d。
[0095]
临界滑动体各层的剪应变γ
l
(t)根据求解式(13)和(17)时获得的节点相对位移u得到，可表达为：
[0096]
[0097]
式中，u
l
为第l个节点在t时刻的相对位移，即向量u的第l个元素。由此可获得所有层的剪应变时间历程，将它们最大值的0.65倍记录为当前的等效剪应变γ
eff，l
(l＝1，2，...，l(i))。
[0098]
5.3根据各层当前的等效剪应变γ
eff，l
更新土体剪切模量g
l
和阻尼比λ
l
。将γ
eff，l
代入一种同场地土体条件相匹配的剪切模量比g/g
max
和阻尼比λ随剪应变变化的经验模型估计新的剪切模量g
l
和阻尼比λ
l
。其中，g
max
为空间随机分布的土体小应变剪切模量，包含在随机场样本x(i)中。
[0099]
5.4重复子步骤5.1-5.3直至等效线性算法收敛，即各层当前等效剪应变与前一次等效剪应变间的差值均小于一个预设阈值ε：
[0100][0101]
当式(19)不满足时，记录当前有效剪应变否则，退出算法，将此时的d作为随机场样本x(i)和地震动加速度时程对应的临界滑动体永久位移d
(i，k)
。循环n
×
m次执行子步骤5.1-5.4，即i＝1，2，...，n，k＝1，2，...，m。
[0102]
步骤6，对所得的永久位移样本进行统计分析，计算在某个设计基准期td(如50年)下滑动体永久位移对设计位移阈值d的超越概率pd。
[0103]
其中，在步骤6中，假设地震发生服从泊松分布，则超越概率pd可表示为：
[0104][0105]
式中，tg为所选地震动的重现期；i[
·
]为指示函数，当d
(i，k)-d＞0时取1，当d
(i，k)-d≤0时取0。在具体边坡工程设计中，位移阈值d的选取取决于多种因素(如工程重要性程度)，越小的取值代表更为安全和安全的工程设计。因此，可考虑多个阈值重复运算式(20)绘制超越概率pd随位移阈值d变化的趋势曲线，全面反映设计边坡的抗震性能。
[0106]
为了更好地说明本发明，下面结合具体实施例作进一步说明，但不作为本发明的限定。
[0107]
一、边坡模型简介
[0108]
考虑一个两层不排水粘土边坡模型，土层分界线位于地表以下18m，边坡剖面和随机场单元划分如图5所示。边坡坡高为24m，坡角为36.9
°
。本实施例考虑4个不确定性土体参数，分别为上层的不排水抗剪强度s
u1
和小应变剪切模量g
max1
、下层的不排水强度s
u2
和小应变剪切模量g
max2
，全部服从对数正态分布，相应的统计特征见表1。两层土体参数间互不相关，其随机场分别独自模拟。将不排水抗剪强度与小应变剪切模量在正态空间的互相关系数ρ
cross，z
设为1。其它确定性参数如下：土体重度γ为19kn/m3，土体剪切模量比g/g
max
和阻尼比λ随剪应变的变化由文献(vucetic m，dobry r.effect of soil plasticity on cyclic response[j].journal of geotechnical engineering，1991，117(1)：89-107)中塑性指数等于30的经验模型表征。
[0109]
表1不确定性参数的统计特征
[0110][0111]
二、具体实施流程
[0112]
步骤1，根据表1中各不确定性参数的统计特征，分别构造各土层随机场变量的联合概率密度函数f(x
up
)和f(x
low
)；
[0113]
边坡上层有2个不确定性土体参数和486个随机场单元，对应972个随机场变量这些变量服从多元对数正态分布，但取对数后则服从多元正态分布，对应的均值向量和协方差矩阵可分别表示为：
[0114][0115][0116]
式中：
[0117][0118]aup
为上层随机场单元的自相关系数矩阵，由式(4)获得。
[0119]
因此，x
up
取值为的概率密度函数f(x
up
)可具体表示为：
[0120][0121]
边坡下层2
×
378＝756个随机场变量的联合概率密度函数的构造与上述过程类似，在此不再赘述。
[0122]
步骤2，基于两层参数随机场变量的联合概率密度函数和利用matlab软件
的随机数生成器分别产生两层土体参数随机场的样本x
up(i)
和x
low(i)
，i＝1，2，...，5000；
[0123]
步骤3，构建边坡稳定性分析模型，计算随机场样本x
up(i)
和x
low(i)
对应的边坡临界加速度系数i＝1，2，...，5000，并记录临界滑动体的位置；
[0124]
3.1将随机场样本x
up(i)
和x
low(i)
作为输入参数分别带入到边坡模型上层和下层对应的材料单元中，设置多条潜在滑动面以基本涵盖整个边坡区域，利用matlab软件进行基于简化毕肖普法的边坡稳定性安全系数计算；
[0125]
3.2对于每个潜在滑动体，采用二分法搜索使边坡安全系数为1时的临界水平地震力系数k
ch
，搜索的范围为k
ch
∈[0.01，2]。将所有潜在滑动体对应k
ch
的最小值记录为临界加速度系数降临界滑动体的位置和尺寸储存为matlab数组，供滑动位移计算使用；
[0126]
循环执行子步骤3.1和子步骤3.2，直至i＝5000。所获得的5000个kc值如图6所示。
[0127]
步骤4，根据重现期ta＝475年和ii类场地土拟定目标反应谱曲线，挑选与该反应谱匹配的30条地震动加速度时间历程k＝1，2，...，30，目标反应谱和所选地震动反应谱如图7所示。
[0128]
步骤5，构建边坡临界滑动体永久位移分析模型，计算随机场样本x(i)和地震动加速度时程下的临界滑动体永久位移d
(i，k)
，其中i＝1，2，...，5000，k＝1，2，...，30。采用等效线性算法考虑土体剪切模量和阻尼比随剪应变的非线性变化，各土层的起始阻尼比取为0.03；
[0129]
5.1把第i组随机场样本对应的临界滑动体理想化为一个多自由度体系，构建质量矩阵m、刚度矩阵k、阻尼矩阵c；
[0130]
5.2建立多自由度体系的运动方程，求解临界滑动体各层的等效剪应变与底部的滑动位移d
(0)
，l和表示该滑动体分层编号和分层数目。体系的高度取为该滑动体最低点与地面之间的距离，滑体每层的厚度设置为0.5m，与随机场单元的竖向长度相等，层数则取为该距离与厚度的商；
[0131]
5.3根据各层当前的等效剪应变更新土体剪切模量g
l
和阻尼比λ
l
。由子步骤5.1和子步骤5.2计算更新后的滑动体底部滑动位移d
(1)
和土层等效剪应变和土层等效剪应变
[0132]
5.4重复子步骤5.1-5.3直至各层与间的绝对差值均小于阈值ε＝1
×
10-4
。将此时的d
(1)
作为最终的临界滑动体位移d
(i，k)
；
[0133]
循环执行子步骤5.1-5.4，直至i＝5000，k＝30，共获得150000组位移数据。
[0134]
步骤6，对所得永久位移样本进行统计分析，计算在设计基准期td＝50年的滑动体永久位移对阈值d＝20cm的超越概率。由于满足该条件的样本有71613个，本实例边坡的位移超越概率p
20cm
可根据式(20)计算为：
[0135][0136]
为了更全面地反映设计边坡的抗震性能，考虑阈值范围d∈[10cm，100cm]并重复带入式(20)可绘制超越概率pd随位移阈值d变化的趋势曲线，如图8所示。边坡永久位移超
过某个阈值d可视为无法接受的边坡性能，所以工程师可根据某个边坡设计方案对应的位移超越概率pd判断该方案是否具备足够的安全裕度，进而决定是否采取削坡、加固等措施提高边坡抗震性能。
[0137]
以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，均应当包含在本发明的保护范围内。