0.   引言

全球地震活动强度自21世纪以来明显提高，而中国又处于环太平洋地震带与欧亚地震带的交汇部位，中国地震局官方数据显示，自2000年以来中国发生7级及以上地震11次，地震灾害造成的损失触目惊心。中国西北地区黄土分布广泛且与高烈度地震区高度重合。天然黄土的震陷性及动力易损性使得黄土地区地震灾害十分严重，因此黄土的动力特性及动力响应破坏问题引起了人们高度重视^[1]。

黄土地区沟壑纵横，强震引起地表裂缝遍布和重大滑坡灾害，地震边坡稳定问题日趋突出，基于平衡理论的拟静力法、Newmark滑块分析法及有限元法等边坡稳定分析方法得以发展^[2]。继海原地震之后，众多学者对地震黄土滑坡形成机制进行了研究，按机理将黄土滑坡分为液化型、剪切型和震陷型^[1]。

黄土具有高孔隙比的架空结构，震陷是地震作用下黄土联接和排列架空结构崩塌引起土体急剧变形和强度大幅损失的力学现象，黄土的震陷量评价是判定工程是否产生灾害的重要依据^[3]。目前，对于黄土震陷的定量预测方法主要分为3类，一类是基于室内试验得到残余应变经验式来预测场地震陷变形^[4]，另一类是通过振动模型试验评价场地震陷特性。曹继东通过一系列振动台试验提出填海土的震陷经验公式，并对场地震陷性进行了评价^[5]。第3类是基于有限元计算方法对场地进行动力模拟。陈正汉等从动三轴试验、振动台试验以及有限元计算3个方面综合评价了厦门填海土的地基震陷情况^[6]。

地震是一种不常见且无法重现的自然灾害，尤其针对震陷的研究很难找到实例^[7]。动力离心模型试验能够更好的还原岩土工程的动力响应问题，是一种有效的研究手段。甘肃省80%地区被不同厚度黄土所覆盖，绝大部分为丘陵地形，土体含水量低，自然坡度平缓，易发生震陷型边坡失稳^[8]。本文以黄土丘陵区震陷型黄土滑坡为背景，研究了边坡的动力响应，以边坡震陷预测为目的，应用动单剪试验条件下的黄土震陷系数经验公式，进行有限差分非线性动力分析，提出一种新型估算黄土边坡震陷量的方法，并将计算结果与动力离心模型试验结果进行了对比，验证了该计算方法的合理性和可行性。

1.   边坡动力分析模型及参数确定

1.1   试验方案

试验采用中国水利水电科学研究院岩土所LXJ-4-450 g·t型土工离心机及配套的R500B型振动台^[9]。模型箱为单向层状剪切模型箱，该箱体在振动方向可自由滑动，用于模拟动力作用下土体的剪切变形，在非振动方向上则为固定边界。有研究表明，层状剪切模型箱可更好的模拟天然场地的自由振动^[10]。为了研究黄土边坡动力响应的一般规律，考虑兰州当地震陷型黄土层厚度及模型箱的限制，将几何相似比尺确定为1/20，并设计出原状黄土边坡概化模型（图 1）。试验前将原状黄土切削成设计尺寸，用直径20 mm、长800 mm的特制微型长柄洛阳铲钻孔，并在预定位置埋设加速度传感器（量程50g，线性度≤5%），模型箱上部安装钢支架及激光位移传感器（量程35 mm，精度0.01 mm）^[11]。

图  1  动力离心试验模型设计及传感器布置

Figure  1.  Design and arrangement of sensors for dynamic centrifugal model tests

下载: 全尺寸图片 幻灯片

试验用土取自兰州震陷型黄土场地的天然原状黄土，其物性指标见表 1。输入地震波为单向水平地震波，采用紫坪铺站台实测的汶川地震加速度时程曲线（图 2）。试验时首先将离心机加速运转达到20g并保持稳定，待试样固结变形稳定后将地震波加速度幅值调整为0.3g和0.5g并进行动力输入。

表  1  兰州黄土基本物理性质指标

Table  1.  Basic physical properties of Lanzhou loess

天然密度ρ/(g·cm-3)	含水率w/%	干密度ρd/(g·cm-3)	液限wL/%	塑限wP/%	塑性指数IP/%
1.38	8.2	1.27	28.2	17.2	11.0

下载: 导出CSV 

| 显示表格

图  2  汶川地震紫坪铺站台实测地震波（PGA≈0.2g）

Figure  2.  Measured seismic waves at Zipingpu platform of Wenchuan Earthquake (PGA≈0.2g)

下载: 全尺寸图片 幻灯片

1.2   数值模型及边界条件

为便于动力数值模拟与离心模型试验间的对比验证，按照相似比将边坡模型放大为原型场地尺寸建模（20倍），采用六面体网格和楔形体网格进行剖分，如图 3所示。

图  3  黄土边坡有限差分数值模型

Figure  3.  Finite difference numerical model for loess slope

下载: 全尺寸图片 幻灯片

FLAC^3D采用完全非线性的分析方法，基于显式差分，使用周围真实密度得出网格节点的集中质量来求解全部运动方程，动力计算可采用任意本构模型，其参数与静本构模型参数一致^[12]。本文静、动力计算皆采用Mohr-Coulomb模型，体积模量和剪切模量可根据弹性模量和泊松比计算得到，其余参数由常规土工试验得到，表 2汇总了模型计算参数。静力计算采用位移边界条件；动力计算时模型底部为黏滞边界，侧面为自由场边界；顶部为自由面（图 3）。黏滞边界在模拟模型底部静止的同时可以吸收入射波，自由场边界则提供了与无限场地相同的效果以保证动荷载在向上传播过程中在边界上不会扭曲。计算阻尼采用局部阻尼，按照经验，岩土体的临界阻尼比在2%~5%，本文黄土临界阻尼比取5%，计算得到局部阻尼系数为0.15。计算时记录到模型底部中心监测点的加速度响应与输入地震波一致，证明动力输入正确。

表  2  黄土边坡数值计算参数

Table  2.  Numerical parameters of loess slope

体积模量/MPa	剪切模量/MPa	泊松比	密度ρ(g·cm-3)	黏聚力c/kPa	内摩擦角φ/(°)
27.5	12.7	0.3	1.38	30	26

下载: 导出CSV 

| 显示表格

2.   地震作用下黄土边坡动力响应规律

2.1   加速度响应规律

为描述黄土边坡的加速度响应规律，本节以监测点响应加速度峰值与输入地震波的加速度峰值之比作为加速度放大系数进行分析。

加速度时程曲线得到不同加速度幅值的地震波作用下边坡坡面的加速度放大系数如图 4所示。动力数值模拟和离心模型试验结果均表明，沿坡面自下而上各个监测点的加速度放大系数呈递增的趋势，表明沿竖向的放大效应。坡脚（A15）与加速度输入位置接近，放大系数最小，在1.1左右，坡肩（A7）的放大系数最大。试验测得在加速度幅值为0.3g的地震波作用下，坡肩最大加速度放大系数为3.38，0.5g作用下的放大系数反而减小，为3.10。可见，边坡对较小地震的动力放大效应更加显著。

图  4  黄土边坡沿坡面加速度响应规律

Figure  4.  Time-history curves of horizontal displacement of slope surface at monitoring points

下载: 全尺寸图片 幻灯片

图 5为沿坡肩铅垂线上不同高度的加速度响应，加速度放大系数随高度的增加逐渐增大，坡体内监测点放大系数小于相同高度的坡面监测点放大系数，临空面的放大效应显著。数值模拟结果与试验结果所得加速度放大系数曲线较为接近，相互验证了黄土边坡动力响应的合理性。在坡中以上区域，模型试验所得加速度放大系数略大于数值计算，这是由于动力离心模型试验所用原状黄土具有天然结构性，当土体含水率较低时，边坡整体刚度较大，能更好的传递地震波，动力放大效应更强。坡肩放大系数可达到3以上。

图  5  黄土边坡沿坡肩铅垂线加速度响应规律

Figure  5.  Acceleration response laws of loess slope along vertical line of slope shoulder

下载: 全尺寸图片 幻灯片

2.2   位移分布规律

数值计算记录了幅值为0.3g的地震荷载下坡面不同高度监测点的水平位移时程曲线（图 6）。地震作用下边坡的水平位移发展可分为4个阶段：①弹性振动阶段：在0~3 s内，地震加速度幅值较小，水平位移在0附近波动；②缓慢增长阶段：在3~14 s内，加速度幅值随时间增大，塑性位移缓慢累积；③加速增长阶段：在第14秒时地震波加速度达到峰值，塑性位移加速发展；④稳定阶段：28 s后大幅值地震波已结束，位移增长减缓并趋于稳定。该图反映出在地震过程中土体位移呈动态波动变化且位移发展总体趋势逐渐累积增长。

图  6  坡面监测点水平位移时程曲线

Figure  6.  Time-history curves of horizontal displacement of slope surface at monitoring points

下载: 全尺寸图片 幻灯片

根据不同高度位移曲线可知，A15处于坡脚位置，水平方向上波动位移始终处于弹性阶段，无塑性累积位移；A6为坡中下侧监测点，处于滑移体下侧，永久水平位移最大；A11为坡中上侧监测点，处于滑移体上侧，水平位移略小于A6处；A7为坡肩位置，其变形主要为竖向压密变形，水平位移较小。由振动结束后边坡的水平位移等值线分布（图 7）可知边坡的水平位移主要集中在坡体下部，与静力状态下边坡的破坏模式类似。

图  7  水平位移等值线分布

Figure  7.  Contours of horizontal displacement

下载: 全尺寸图片 幻灯片

3.   黄土边坡震陷变形分析

3.1   黄土震陷系数经验公式与震陷量估算方法

常规的数值分析方法难以考虑黄土的结构性，很难模拟黄土的震陷变形，本文基于有限差分非线性动力反应分析，引入动单剪试验条件下的黄土震陷系数经验公式，对黄土边坡的震陷变形进行计算。

西安理工大学自主研发的新型立方体铰接机构动单剪仪，具有加载条件更接近实际地震作用，消除接触板与试样间的摩擦效应，剪应力分布更均匀等优点^[13]。吴飞洁用该仪器对不同含水率原状黄土进行动力试验，提出了将压缩应力${\sigma _{\text{v}}}$、含水率w、等效剪应变幅值${\gamma _{\text{d}}}$以及等效振次N作为影响的黄土震陷系数经验公式^[14]。本文引用吴飞洁的震陷公式，根据试验结果，震陷系数${\delta _{\text{d}}}$与剪应变幅值呈二次多项式函数关系，与等效振次N呈对数函数关系，与含水率呈线性关系，考虑多种影响因素，将震陷系数公式用式（1）表示，其中α和β均为土性参数。通过分析参数与压缩应力的关系，参数α在不同压缩应力下变化不大，可取定值-0.012，参数β与$({\sigma _{\text{v}}}/{p_{\text{a}}})$呈线性关系，其中${p_{\text{a}}}$为标准大气压。

将式（2）代入式（1）：

土单元的压缩应力为土体中的土压力，即${\sigma _{\text{v}}} =$ $\rho gh$，代入式（3）并分离变量h即可得到震陷系数与深度h的关系如式（4）所示，其中A和B整理成式（5）和式（6）。

根据动单剪试验条件下的黄土震陷系数经验公式提出了黄土场地震陷量估算方法，具体步骤如下：

（1）将待计算的天然黄土场地用等间距的水平线和铅垂线划分为若干土单元。

（2）利用FLAC3D计算黄土场地的动力响应，并记录各土单元的动剪应变时程曲线，根据等效荷载法得到相应的等效剪应变${\gamma _{\text{d}}}$（取最大剪应变的0.65倍）和等效振次（汶川地震为8.0级地震，等效振次取30次）。

（3）将各土单元参数代入式（4）计算其震陷系数，乘以该土单元高度得到该土单元的震陷量。

（4）将具体土单元以及该土单元以下处于同一条铅垂线上土单元的震陷量进行求和，得到该土单元的累积震陷位移。

（5）根据线性插值法绘出地层的震陷变形等值线，并进行场地震陷性评价。

3.2   黄土边坡震陷量计算分析

黄土边坡的震陷变形可简化为平面应变问题，本节对输入地震加速度峰值为0.3g作用下黄土边坡的震陷量进行计算。

本文黄土边坡为高度10 m的单一均质土体，为方便计算，首先用间距1 m的水平线和铅垂线将其划分为边长1 m正方形土单元，按照上述计算方法得到每个土单元的累积震陷量，进而得到整个黄土边坡的震陷变形。以监测点A7所在土单元为例，根据数值计算得到的剪应变时程曲线（图 8）可知最大剪应变${\gamma _{{\text{max}}}}$为6.43%，取等效剪应变${\gamma _{\text{d}}}$为4.18%，等效振次N取30次，深度h为1 m，含水率w为8%，将其代入式（4）可得该层的震陷系数为4.25%，土单元高度为1 m，则该土单元震陷量为4.25 cm。同理可得同一铅垂线上土单元的震陷系数及震陷累积变形如图 9所示。按照同样方法计算所有土单元的累积震陷变形，并用线性插值法得到黄土边坡的震陷变形等值线图（图 10）。

图  8  监测点动剪应变时程曲线（A7）

Figure  8.  Time-history curve of dynamic shear strain of monitoring points

下载: 全尺寸图片 幻灯片

图  9  不同高度震陷系数及震陷累积位移

Figure  9.  Seismic subsidence coefficients and cumulative displacements at different elevations

下载: 全尺寸图片 幻灯片

图  10  黄土边坡震陷变形等值线图

Figure  10.  Contours of seismic subsidence deformation of loess slope

下载: 全尺寸图片 幻灯片

根据计算结果，图 9（a）表明土体震陷系数随高程增高呈增大速率先快后慢的对数性增大，含水率与等效振次认为常数，震陷系数仅与土压力和动剪应变有关。位于较深土层的土单元压缩应力较大，动剪应变较小，深层土的震陷系数主要受压缩应力的影响；对于较浅土层，土体动剪应变响应较大，土的震陷系数较大，但土单元压缩应力较小，震陷系数增长速率减小，说明浅层土震陷系数受动剪应变影响较大。综合上述分析，动剪应变决定土的定震陷系数的大小，土压力决定土的震陷系数的变化速率。图 9（b）为震陷量沿高度的累积过程，震陷量是土层厚度乘以震陷系数，即震陷性黄土层厚度是震陷量的关键因素。

根据黄土边坡震陷变形等值线图（图 10）可知，随着边坡高度增加，累积震陷变形逐渐增大，直至坡顶达到峰值。同一水平位置，越靠近临空面的土体压缩应力越小，震陷量也越小，靠近坡顶等值线间距越小，表明坡体上部的震陷压缩量大于下部。

图 11为动力离心模型试验监测到的竖向变形发展过程。地震作用下黄土边坡滑移变形较小时，可近似认为竖向监测变形都是由震陷变形导致。在地震历时12 s前，地震波加速度幅值较小，边坡处于弹性振动阶段；12 s到14 s之间，竖向位移开始发展；自14 s开始，加速度幅值达到峰值后，竖向位移加速增长；23 s后竖向位移增长随之减缓，并逐渐趋于稳定，但此时的地震加速度幅值相对较大，说明震陷变形在达到一定程度后并不再随时间持续发展。坡顶最终沉降为18 mm，坡面上部为7 mm，坡面下部和坡底的竖向位移近乎为0，按照相似比还原后坡顶的最大沉降为0.36 m，坡面上部为0.14 m。而通过计算得到的坡顶最大震陷变形为0.32 m，坡面上部为0.14 m，坡脚和坡底仅有0.01 m的沉降，可见通过震陷系数经验公式所得的边坡震陷变形与试验结果基本一致，验证了该计算方法的合理性和可行性。

图  11  地震作用下竖向位移时程曲线（0.3g）

Figure  11.  Time-history curves of vertical displacement under earthquakes

下载: 全尺寸图片 幻灯片

《震害调查及地震损失评定工作指南》中规定：震陷量小于2 cm时为基本完整；2~4 cm时为轻微破坏；4~8 cm时为中等破坏；8~40 cm时为严重破坏；大于40 cm为失稳破坏^[1]。根据评定标准，该边坡在加速度幅值为0.3g的水平地震波作用下坡顶形成严重震陷破坏，沿坡面向下逐渐过渡为中等破坏和轻微破坏，坡脚和坡底则基本完整。

3.3   震陷型黄土边坡破坏形式

地震作用下黄土边坡的总变形是水平滑移变形与竖向震陷变形双向耦合的结果，在架空孔隙结构、含水率、土体应力等综合因素的影响下，水平方向上土体位移向临空面一侧发展，坡中下侧鼓起明显，呈现出滑动的趋势；竖直方向上原状孔隙结构破坏后土颗粒填充大孔隙，导致土体产生竖向的永久塌陷，土层越厚，下陷位移越大。

离心机试验表明（图 12）地震荷载作用下黄土边坡的塑性变形主要为震陷压密，坡体内部形成潜在滑移面，边坡整体水平滑动位移较小，坡面和坡顶张拉裂隙发育显著，根据动力离心模型试验及数值计算结果绘制了黄土边坡震陷失稳的变形破坏特征（图 13）。

图  12  离心机试验黄土边坡破坏特征

Figure  12.  Failure characteristics of loess slope in centrifugal model tests

下载: 全尺寸图片 幻灯片

图  13  震陷型黄土边坡破坏特征

Figure  13.  Failure characteristics of loess slope of seismic subsidence

下载: 全尺寸图片 幻灯片

将震陷前后边坡进行对比，天然黄土垂直节理发育，坡顶震陷显著；同时，坡面向临空面一侧鼓出，导致坡顶张拉裂隙发育较深，震陷变形不均导致大量水平错位裂缝的发展，部分土体脱落并在坡脚堆积，坡体内部形成潜在的滑动面，但未发生滑坡。

在黄土覆盖层较厚的丘陵地区，自然斜坡角度平缓，含水率较低，易发生震陷型边坡失稳。张振中等对1995年永登地震引起的边坡震陷破坏进行了调查，该地震仅5.8级，但由于结构性黄土的动力易损性，山坡裂隙遍布，坡顶土体酥裂，震陷沉降不均导致坡面形成错位阶梯，土体仍附在坡面尚未下滑^[15]，而震后边坡可能由于降雨入渗、再次地震等原因诱发滑坡。

4.   结论

黄土震陷是地震动引起的附加沉陷。本文针对黄土边坡进行了非线性动力计算，通过动单剪试验的黄土震陷系数经验公式，提出了黄土震陷量的估算方法，并根据动力离心模型试验和动力数值分析讨论了震陷型黄土边坡的变形破坏形式，主要得到以下4点结论。

（1）黄土边坡对地震动具有放大效应。加速度放大系数沿高度逐渐增大，且同一水平高速坡体内部的放大系数小于坡面，表现出临空面效应，加速度幅值为0.5g作用下的放大系数小于0.3g。

（2）地震作用下黄土边坡的水平变形与静力条件下的边坡类似。含水率较低时天然黄土的结构性较强，整体刚度较大，坡面向临空面一侧隆起，坡内形成潜在滑动面，滑坡趋势明显但未发生大距离滑动。

（3）采用动单剪试验下的黄土震陷系数经验公式对边坡震陷量进行预测，其结果与动力离心模型试验十分接近，验证了该计算方法的合理性。结果表明，黄土边坡的震陷变形与震陷型黄土层厚度关系较大，震陷变形等值线近似呈水平线，靠近临空面时震陷量减小。随着高程的增大，土层震陷系数逐渐增大。

（4）结合动力离心模型试验和数值计算结果讨论了震陷型黄土边坡的破坏形式。边坡水平向临空面一侧变形的同时，震陷引起的竖向沉降变形显著，同时坡顶张拉裂隙和坡面错位裂隙大量发育，震陷沉降不均导致坡面形成错位阶梯。