1. 引言
我国西部因其脆弱的地质环境和强烈的地震活动往往会诱发大的崩塌和滑坡,随着经济的高速发展,越来越多的铁路、公路及国防隧道建设重心逐渐往西部高烈度山区转移,西部山区隧道工程的规模与数量日益增长,滑坡等地质灾害的影响已不可忽视[1]。边坡在西部山区十分常见,层状边坡更是分布广泛,边坡按节理倾向和坡面倾向的关系可分为顺层边坡和反倾边坡[2]。其中,倾角大于坡角的顺层边坡是最易造成严重破坏的边坡。此类边坡在自然条件下较为稳定,但在地震作用下极易失稳,形成高速远程滑坡,对生命财产安全造成极大的威胁[3]。在“5.12”汶川地震中失稳的四川省阿坝州茂县叠溪镇新磨滑坡、四川省绵阳市安县千佛岩干磨房滑坡以及四川省绵阳市北川县唐家山滑坡[4]。因此,研究高烈度地震区顺层边坡的动力响应和失稳机制具有重要的理论和现实意义。
对地震作用下顺层边坡的研究,国内外学者进行了诸多探讨。数值模拟技术能够直观的反映边坡的破坏机制和动力响应,因此成为大多数学者研究边坡的主要手段。李龙起等[5]采用三维离散元技术以汶川地震中干磨坊滑坡和水磨沟滑坡为原型开展两种软硬互层斜坡对比分析,研究不同岩层倾角陡倾顺层软硬互层斜坡在地震作用下的动力响应及失稳机理,揭示了此类结构岩体边坡的动力破坏形式。黄志全等[6]通过FLAC 3D数值模拟方法,研究分析了不同波形、加载时间的地震波作用下陡倾顺层岩质边坡的动力响应规律以及变形破坏模式。陈文倩[7]等通过数值模拟,分析了结构面参数对顺层岩质边坡动力响应影响程度。
上述研究主要集中在地震作用下顺层边坡的加速度放大效应和破坏模式方面,关于边坡的内部动力响应影响性研究并不深入,还有进一步探讨空间,并且研究主要从时域或频域方面进行分析,未能进行更综合的研究。借助希尔伯特–黄变换,处理坡体的地震信号,可获取地震信号在三维时频域的能量分布,进行更加可靠的动力响应特征分析。杨长卫等[8]利用Hilbert谱和Hilbert边际谱分析地震能量在时频域上的传播特征,对比顺层及反倾岩质边坡动力响应差异并识别出两种边坡地震损伤的主要集中位置。雷浩等[9]利用HHT方法对隧道交叉穿越边坡时隧道及边坡的频谱特性进行定量分析,获取地震作用下隧洞及边坡的敏感区域,对边坡工程具有很强的指导意义。
本文以顺层边坡为研究对象,设计完成了顺层边坡的数值模拟,研究其在地震作用下的动力响应。利用HHT方法获取边坡在加速度时频域的响应特征,分析边坡的地震动力响应特征和破坏机制。
2. 数值模拟
2.1. 模型建立
为分析顺层边坡在地震作用下的动力响应特征及失稳过程,本文同时使用有限元软件进行数值模拟,对地震作用下该边坡动力响应及失稳过程进行分析,数值模拟模型图如图1所示。边坡坡度为45˚,岩层倾角30˚,软弱面厚度为2 m。在构建数值模拟模型时,空间网格单元尺寸Δl遵循Kuhlemeyer和Lysmer提出的规则[10],以保证仿真的准确性。其公式如下:
(1)
其中Δl为网格单元尺寸,λ为输入波最高分量频率对应的波长。因此,我们将网格尺寸设置为1 m,以保证地震波在介质中的完整传播,最终划分模型共含9529个节点,9943个网格。模型坡体和软弱面都采用塑性本构模型,选取摩尔库伦屈服准则进行分析,数值模拟参数依照前人经验确定[11]。如表1所示。加速度监测点布设如图2。
Figure 1. Numerical simulation model diagram (unit: m)
图1. 数值模拟模型图(单位:m)
Figure 2. Monitoring point layout map (unit: m)
图2. 监测点布设图(单位:m)
Table 1. Numerical simulation parameters
表1. 数值模拟参数
位置 |
密度(g/cm3) |
抗拉强度(kPa) |
弹性模量(MPa) |
泊松比 |
内摩擦角(˚) |
黏聚力(kPa) |
岩体 |
2.50 |
4.00 |
80 |
0.25 |
35 |
70 |
软弱面 |
2.30 |
0.8 |
20 |
0.35 |
20 |
50 |
2.2. 加载工况
在模型建立后,对模型底部以及四周边界水平向速度进行约束。在自重应力下达到初始平衡后,设置四周为自由场边界,底部为黏滞边界,消除应力波的反射,如图3所示。
Figure 3. Numerical simulation boundary conditions (unit: m)
图3. 数值模拟边界条件(单位:m)
地震动力荷载从模型底面输入,有限元软件可通过应力边界及速度边界进行动力加载,文章采用应力边界进行动力加载,边界应力按式(1)与(2)进行转换:
σn = −(2ρCp)vn (2)
σs = −(2ρCs)vs (3)
上式σn为法向应力(MPa),σs为切向应力(MPa),ρ为介质密度(kg/m3),vn为输入的质点的铅直方向振动速度(m/s),vs为输入质点水平方向振动速度(m/s)。
在计算过程中阻尼选用局部阻尼来消除振动时产生的动能,最终设置阻尼系数为0.157。本次试验主要进行了1 g,2 g,3 g和4 g这4种不同地震峰值的地震波加载,加载地震波形选取汶川波,该地震波是实测地震波,能真实反映实际地震动特征。加载方向均为水平向,试验加载方案如表2所示,汶川波时程曲线和频谱曲线如图4所示。
Table 2. Test loading scheme
表2. 试验加载方案
编号 |
波 |
方向 |
幅值(g) |
1 |
汶川波 |
水平向 |
1 |
2 |
汶川波 |
水平向 |
2 |
3 |
汶川波 |
水平向 |
3 |
4 |
汶川波 |
水平向 |
4 |
(a) 汶川波时程曲线 (b) 汶川波频谱曲线
Figure 4. Time history curve and frequency spectrum curve of Wenchuan wave
图4. 汶川波的时程曲线和频谱曲线
3. 结果
3.1. 加速度响应规律分析
为了探究边坡产生变形失稳的主要原因,本文对加速度展开分析。因加速度可以反映地震惯性力,故边坡加速度变化及其分布规律是评价边坡地震动力响应的基本信息。在此,参照前人研究,定义各监测点动力响应的加速度峰值(AE)与原始输入波的加速度峰值(AC)的比值为加速度放大系数(ƞ),如下表示:
(4)
在此,本文主要研究地震幅值为1 g加速度响应规律,便于反应原始地震波作用下的边坡动力响应特征。数值模拟结果表明:随着高程增加,坡面的加速度放大系数呈现高程放大效应;随着水平距离增加,坡体内部的加速度放大系数随之增加,即呈现出趋表效应。如图5所示。这种现象是由于地震波在顺层边坡传播过程中在软弱面发生反射、折射和散射而产生的效应。
(a) (b)
Figure 5. Change law of acceleration amplification coefficient of slope
图5. 坡体加速度放大系数变化规律图
3.2. 变形破坏过程分析
为研究顺层边坡失稳过程,对数值模拟模型位移云图展开分析(如图6所示),发现其失稳过程如下:
当地震波振幅为1 g时,边坡在地震波作用下坡角软弱面处出现明显破坏;在地震波振幅为2 g时,随着地震波振幅的逐渐增大,坡体破坏沿着软弱面逐渐向上扩展,同时位移也进一步增大,坡体结构的稳定性进一步下降;当地震波振幅为3 g时,随着地震波振幅的逐渐增大,滑面沿着软弱面贯通,坡体发生滑动破坏;当地震波振幅为4 g时,坡体沿滑面进一步滑动,达到完全破坏。
Figure 6. Numerical simulation displacement nephogram: (a) 1 g seismic wave, (b) 2 g seismic wave, (c) 3 g seismic wave, (d) 4 g seismic wave
图6. 数值模拟位移云图:(a) 1 g地震波,(b) 2 g地震波,(c) 3 g地震波,(d) 4 g地震波
4. 讨论
由于数值模拟试验无法清晰表征地震波作用下顺层边坡内部响应,故引入希尔伯特–黄变换理论,深入研究在地震波作用下边坡内部能量变化,以此揭示边坡失稳机理。
4.1. 希尔伯特–黄变换理论
希尔伯特–黄变换(HHT)是由N. E. Huang等于1998年提出来的一种随机信号处理方法,主要包括经验模态分解(EMD)方法和Hilbert变换2个部分,经验模态分解(EMD)和Hilbert变换的结合被定义为希尔伯特–黄变换,其是一种新颖且功能强大的时频分析方法。
经验模态分解是整个希尔伯特黄变换极为关键的起始步骤。它的目的是把一个复杂的、非平稳的信号自适应地分解成多个本征模态函数(Intrinsic Mode Functions,简称IMFs),即把原始信号的多种不同频率成分拆分。首先,需要找出原始信号中的所有极大值点,然后连接极大值点形成上包络线;同理,找出所有极小值点形成下包络线,计算上下包络线的平均值,用原始信号减去这个平均值,得到中间信号,直到中间信号满足本征模态函数的条件,也就是它的极值点数目和过零点数目相等或者最多相差一个,而且在任意时刻,由局部极大值点和局部极小值点定义的上下包络线的平均值为零,最终提取出了本征模态函数。重复上述流程,持续分解,直到剩余信号变成一个单调函数或者其幅值小到可以忽略不计,完成了经验模态分解,得到多个本征模态函数。
而后,对每个本征模态函数分别进行希尔伯特变换,得到每个本征模态函数对应的瞬时频率和瞬时相位等重要时频信息。通过对各个本征模态函数的时频分析结果进行整合,最终就可以得到原始非平稳信号完整的时频分布情况,便于分析地震波信号的变换情况。
4.2. 希尔伯特谱分析
在进行Hilbert谱分析之前,首先要对地震信号的加速度时程曲线进行EMD分解,以在1 g汶川波作用下测点的加速度时程曲线为例,进行EMD处理后得到IMF分量的时程曲线,而后进行Hilbert变换。
图7为选取部分测点在幅值为1 g压缩4倍汶川波作用下的Hilbert谱结果,整体而言,六个测点在Hilbert谱时频域时间轴上高能量的分布主要集中在3 s和8 s,这是由地震波所决定的,而在Hilbert谱时频域频率轴上,六个测点在10~20 Hz处都有明显的峰值,其主要是受边坡自身特性影响。从部分测点上看,Hilbert谱在频率轴上的分布差异明显,测点A8和A9随着高程的增加,Hilbert谱能量分布上发生较大的变化,其在10~30 Hz处能量相对较大,高频部分能量开始增强。测点A9的Hilbert谱能量较大的部分频带变宽,最高达到了40 Hz,并且高频部分能量激增。
上述现象表明,坡体的高程对地震波的Hilbert谱能量有放大作用,且其对高频部分的能量放大作用效果更为明显。
Figure 7. Two-dimensional Hilbert energy distribution ((a) A0, (b) A5, (c) A6, (d) A7, (e) A8 and (f) A9) of acceleration monitoring points in time-frequency domain under 1 g seismic wave. (Color change represents the value of seismic Hilbert energy)
图7. 1 g地震波作用下加速度监测点时频域二维Hilbert能量分布((a) A0、(b) A5、(c) A6、(d) A7、(e) A8和(f) A9。色变表示地震希尔伯特能量的值)
4.3. 边坡失稳机理分析
根据数值模拟所记录边坡在地震作用下变形破坏的全过程位移云图以及塑性区图(见图8),可看出顺层边坡其宏观破坏过程:坡体软弱面是顺层边坡中的薄弱点,受加速度“高程效应”和“趋表效应”的影响,坡脚软弱面处岩体在地震波的左右摆动下最先被拉裂。随着地震幅值的增大,坡脚处拉裂缝逐步沿着软弱面向上延伸。在地震幅值达到足够大时滑面贯通,顺层边坡失稳。
其内在失稳机理为:在地震作用下顺层边坡的软弱面因抗剪强度相对较低,更容易发生剪切变形和破坏,从而成为边坡滑动的潜在滑动面。在岩土体自重以及地震力联合作用下,地震作用弱化了软弱面的摩擦力,使得原本相对稳定的岩层连接变得松动且增强了岩层所承受的剪切力,最终岩土体沿软弱面产生整体的平面滑动。顺层边坡的滑动面往往平行于软弱面走向,类似刚体沿平面进行滑动,下滑力克服了层面的抗滑力从而破坏。在失稳过程中因坡体的高程会放大地震波高频部分的Hilbert谱,故软弱面上部产生拉裂破坏,而下部因岩土体自重以及地震力联合作用使得下部剪切破坏。
Figure 8. The plastic zone of slope under the action of 4 g seismic wave
图8. 4 g地震波作用下边坡塑性区图
5. 结论
本文通过振动台试验和动力数值模拟试验,研究地震作用下陡倾顺层岩质边坡动力响应,主要得到以下结论:
1) 边坡在地震波作用下表现出了明显的“高程效应”和“趋表效应”,边坡的破坏主要产生在软弱面。随着地震波强度增大,边坡内部开始破坏,当地震波强度达到4 g后,边坡完全破坏。
2) Hilbert谱表明高程对地震波能量的传播有放大作用,尤其是高频部分。随着坡体的高程增加,Hilbert谱能量出现明显放大,且其对高频部分的能量放大作用效果更为突出。
3) 软弱面影响着坡体的稳定性是顺层边坡中的薄弱点,坡角软弱面处最先发生破坏。随着地震幅值的增大,坡体拉裂缝逐步向上延伸,滑面贯通,边坡失稳。
4) 顺层边坡的失稳机理为:软岩相比硬岩易发生变形破坏,不协调变形会在软弱面交界面处产生,促进软弱面的破坏扩展,导致顺层边坡沿着软弱面发生滑移失稳,最终产生软弱面上部拉裂破坏和下部剪切破坏。