西南艰险山区具备地形地貌复杂多变、地震频发及气象条件恶劣等特点。错综复杂的地形地质条件、巨大的风化差异和频发的地震动导致高位危岩落石、崩塌和滑坡灾害频发，严重威胁铁路建设和行车安全^[1-3]。顺层岩质边坡是西南山区分布最为广泛的边坡类型之一，在风化作用下，岩质边坡内部存在大量的非连续结构面，包括层面、节理及软弱夹层等，导致边坡的动力响应与破坏模式也更为复杂^[4-5]。

针对顺层岩质边坡的动力特性与破坏模式的研究方法通常主要包括现场调查、理论分析、数值仿真计算及模型试验等^[6-9]。在研究含大量非连续结构面岩体在地震作用下的响应时，通常采用振动台模型试验或数值仿真的方法。振动台试验方面，董金玉^[10]、Yang^[11]等基于振动台试验研究了顺层、反倾层状岩质边坡的动力响应和破坏模式，得出不同倾向岩质边坡加速度放大系数的变化规律。许彬等^[12]基于振动台试验研究了水库诱发地震时消落带岩体劣化岩质边坡的动力响应与破坏模式，并建立了微震和强震作用下的损伤模型。吴多华等^[13]研究了地震波频率、幅值和持时对边坡动力响应的影响，得出当输入波频率小于自振频率时，边坡的水平向PGA放大系数随输入波频率的增大而增大；当输入波频率大于自振频率时，放大效应减弱。数值仿真方面，刘树林等^[14]在振动台试验的基础上，结合UDEC离散元方法，研究了不同倾角的顺层岩质边坡在微振作用下的动力响应。Feng等^[15]通过DDA非连续介质方法，研究了顺层、反倾层状岩质边坡在横波和纵波作用下的动力响应规律。胡训健等^[16]基于PFC2D颗粒流软件，研究了水平厚层状岩质边坡内部岩桥长度及节理间距对其动力响应与破坏过程的影响。李鹏等^[17]通过UDEC分析了软弱夹层倾角、埋深、弹性模量、厚度等影响因素对岩质边坡动力响应的影响。

已有研究成果大多是针对某种特定岩质边坡的动力特性，或是对节理的物理力学参数展开敏感性研究，对于岩体破碎程度对顺层岩质动力响应的影响研究较少。西南山区存在大量碎裂状的边坡，局部风化程度大，岩体破碎严重，结构面多，而内部风化程度低，岩体破碎程度低，结构面少。为揭示岩体破碎程度不同的顺层岩质边坡的动力响应和破坏模式，本文通过层理和垂直节理切割边坡块体的大小表征岩体的破碎程度，并通过UDEC离散元法，对2维不同节理参数的岩质边坡的地震响应规律进行研究。研究结果可为碎裂状顺层岩质边坡在地震动力作用下的动力响应特性及加固防护提供参考。

1.   工程背景

1.1   工程地质条件

研究工点为西南高陡山区某边坡，两侧山体陡峭，地形险峻，为高山峡谷地貌，地面高程为2 021～3 800 m。实地调查显示，该边坡为顺层岩质边坡，坡度为50°～60°，岩性以片麻岩为主，岩体多为弱风化，局部强风化（强风化区域较小规模的结构面较多，岩体更为破碎），地质剖面如图1所示。岩体中裂隙发育，当遇到突发地震时，可能产生较大范围的破坏。

图  1  边坡地质剖面

Fig.  1  Geological profile of the slope

下载: 全尺寸图片

1.2   结构面分布

在隧道进口端布设一平硐，设计长度为70 m，据现场调查及平硐内测试和统计，平硐内共有结构面147条。节理面玫瑰花图（图2）显示，平洞揭露的结构面以陡倾角为主，主要优势节理组倾向为W10°～30°N或E10°～30°S，其次为倾向E0°～20°N的节理。与边坡表面揭示的结构面优势组对比可见，表面结构面组与内部结构面组发育规律差别较大。

图  2  节理面的玫瑰花图

Fig.  2  Rose chart of jointed surface

下载: 全尺寸图片

1.3   赤平极射投影分析

为进一步研究隧道进口的结构块体特征，应用结构面优势分组的3组平均产状，结合隧道进口岸坡的临空面产状N22°E/37°SE和N22°E/70°SE，进行赤平极射投影分析，如图3所示。由图3可知：3组结构面产状中，J1走向与岸坡面走向近正交，只能起到切割岩体的作用；J3为缓倾结构面，其走向与岸坡面走向也近于正交并偏向于倾向坡内；唯一与岸坡面走向近于一致的J2由于其属于中陡倾角，介于自然坡面倾角间，且切割其他两组节理而形成卸荷裂隙。故根据结构面优势分组平均产状与岸坡临空面的赤平投影分析，局部结构块体可能不稳定，易形成危岩落石。在地震等不利荷载作用下，更易发生崩塌落石灾害。

图  3  赤平投影图

Fig.  3  Equatorial projection

下载: 全尺寸图片

2.   计算模型及验证

2.1   离散元模型的建立

据勘察资料可知，该边坡强风化岩体被结构面切割得较为破碎，层面、节理面对边坡整体稳定性的影响较大。为模拟岩体破碎程度对顺层岩质边坡动力响应的影响，对模型进行一定程度的简化，利用层理和节理切割边坡块体的大小表征岩质边坡的破碎程度，块体越小，表示边坡破碎程度越大^[18]。

综合考虑计算速度和计算结果精度，以及已有算例经验^[18-19]，模型范围确定如下：边坡高度为160 m，顶部宽度为52 m，底部宽度为250 m，单台阶高50 m，坡度设计为60°，层理倾角为60°。为监测模型在坡面及坡体内部的动力响应规律，在坡表布置8个监测点（M1～M8），在M6位置处沿铅锤方向向下布置6个监测点，记录各点处的加速度响应及坡表处的位移时程。离散元计算模型和监测点布置如图4所示。图4中：A表示弱风化岩体区域层理间距，B表示弱风化岩体区域节理间距，C表示强风化岩体区域层理间距，D表示强风化岩体区域节理间距。

图  4  模型及监测点示意图

Fig.  4  Schematic diagram of the model and monitoring points

下载: 全尺寸图片

2.2   本构模型及计算参数设置

在静力分析时，将岩体视为可连续变形体，初始平衡阶段采用Elastic本构模型。在动力分析时，岩体采用Mohr–Coulomb本构模型，分析岩体的塑性变形，其中，弹性模量和泊松比应为岩体的动弹性模量与动泊松比。模型的层理、垂直节理采用接触–库伦滑移本构模型（AERA）。根据工程地质勘察报告，研究区域岩石力学参数见表1。为使得计算结果可以更为精确地反应模型的动力响应，根据已有计算经验，结构面参数根据文献[19]选取，具体见表2。

2.3   边界条件及阻尼设置

在静力计算时，坡顶和坡面为自由边界，模型的左右两侧和底面约束法向位移。在地震动力计算时，边界条件的设置对计算结果有很大影响^[20-21]。为防止输入地震波的反射和失真，将模型的左右两侧变为自由场边界，高度与模型两侧垂直高度相同；在模型底部设置黏滞边界。

图5为模型在动力分析时的边界设置，计算过程中黏滞边界参数与片麻岩参数相同。

图  5  边坡动力模型边界

Fig.  5  Boundary of dynamic slope model

下载: 全尺寸图片

由于模型底部为黏滞边界，故地震波无法以加速度时程方式直接施加，需要将加速度时程积分成速度时程，再将速度时程利用式（1）～（2）转化为应力时程：

式中： ${\sigma }_{{\rm{n}}}$ 为施加在边界上的法向应力荷载，MPa； ${\sigma }_{{\rm{s}}}$ 为施加在边界上的切向应力荷载，MPa； ${v}_{{\rm{n}}}$ 、 ${v}_{{\rm{s}}}$ 分别为输入的法向振动速度与切向振动速度，m/s； $ \rho $ 为岩体密度； ${C}_{{\rm{p}}}$ 、 ${C}_{{\rm{s}}}$ 分别为地震纵波波速和横波波速。

由于局部阻尼计算速度快、效率高，在研究中选用局部阻尼进行动力计算。参考相关文献资料，经过试算分析，阻尼系数取0.15^[19]。

2.4   输入地震波

为对比不同地震波对顺层岩质边坡动力响应的影响，选取的输入地震波类型为EL–Centro地震波和Kobe地震波（简称EL波和Kobe波），对地震波进行归一化处理，两个地震波持续时间均为20 s，地震波的加速度时程曲线及傅里叶谱如图6所示。

图  6  输入波形的加速度时程曲线及傅里叶谱（0.1g）

Fig.  6  Acceleration time history curve and Fourier spectrum of input waveform (0.1g)

下载: 全尺寸图片

查阅《中国地震动参数区划图》，工况类别为Ⅱ类，基本地震动峰值加速度为0.2g，基本烈度地震下设计地面地震最大加速度为0.6g。故本次数值计算地震波地面峰值加速度（PGA）选为0.2g、0.4g和0.6g。

2.5   工况设计

为得出不同节理参数对试验结果的影响，通常情况下需使用控制变量法进行大量计算，受离散元UDEC计算量与计算时间的限制，需要对计算工况进行优化设计。正交试验法是根据正交性从全面性试验中选取部分具备代表性的点进行试验，其结果分析较为方便，更易评估各影响因子间的相互关系。

结合数值计算模型，针对弱风化岩体区域的层理间距A、垂直节理间距B，以及强风化岩体区域层理间距C、垂直节理间距D这4个因素，选取2种变化水平进行研究，根据已有边坡计算实例和工程经验，设计8组试验，不同工况的节理参数取值见表3。

2.6   模型验证

为验证本文数值方法的合理性，对该数值方法进行验证。肖克强等^[22]基于离散元UDEC分析了地震荷载作用下坡高、坡脚等影响因素对顺层岩质边坡位移的影响规律。因此，首先建立与文献[22]中相同的模型，岩层倾角β分别取20°、40°及45°，并将计算结果和文献[22]进行对比，结果如图7所示。由图7可知，当加载频率为2 Hz、振幅为1.25 m/s²、持时为10 s的正弦波时，模型的位移响应规律和文献[22]保持一致，并且位移幅值基本相同，这说明本文所采用的数值模拟计算方法比较合理，得出的结果具有合理性。

图  7  本文位移计算结果同文献[22]对比

Fig.  7  Comparison of displacement calculation results with Reference[22]

下载: 全尺寸图片

3.   模型动力响应分析

3.1   坡体位移

坡体模型位移分析时，实测曲线包含震荡和基线两部分，震荡部分表征坡面位移在地震动作用下的实际响应，基线部分表征坡面位移的变化趋势^[23]。选取地震结束后坡面最终位移（基线部分的最终位移）进行分析，并结合位移云图研究其失稳演化模式，测点位置由相对高程h·H^-1表示，其中，h表示测点到底部的垂直距离，H表示斜坡高度。以工况1为例，图8为该工况下模型的坡面位移。

图  8  模型坡面位移（工况1）

Fig.  8  Model slope displacement (condition 1 )

下载: 全尺寸图片

由图8可知：输入不同地震荷载的坡面位移规律大致相同，从坡底至坡顶呈现先增大再减少的趋势，位移响应最大处位于测点M4处，即强风化岩体的下部；输入地震波的幅值越高，坡面的最终位移越大。同一震级下，加载Kobe波时边坡坡面的位移响应大于加载EL波时的位移响应。

3.2   坡体加速度响应

定义PGA（peak ground acceleration）放大系数为数值模型测点加速度峰值与输入加速度峰值的比值^[24]。本文主要从坡表（M1～M8）及坡体内部（M6～M14）分析水平向PGA放大系数的变化规律。以工况1为例，模型各监测点的PGA放大系数如图9所示。

图  9  PGA放大系数变化规律（工况1）

Fig.  9  Variation law of PGA amplification coefficient (condition 1 )

下载: 全尺寸图片

由图9（a）可以看出，PGA放大系数沿坡表呈现先增大后减小的规律，表现出“高程放大效应”，这与杨长卫^[25]和付晓^[26]等所得结论基本相同。值得注意的是：在下方岩体破碎程度低的区域，随着高程的增加，PGA放大系数有所增加，但增加幅度较小；上方岩体破碎程度高的区域，PGA放大系数随着高程的增加，呈现先增大后减小的趋势，且增加幅度较大，PGA响应的最大处位于强风化区域坡面的中上部。这种现象可能是因为与强风化岩体区域相比，下方岩体节理密度低，完整性相对较好，在进行地震波加载时，下方岩体区域整体的PGA放大系数相对较小；而强风化岩体区域节理密度高，完整性较差，导致其PGA放大系数变大。同时可以发现，随着输入地震波幅值的增加，模型的PGA放大系数逐渐减小。这可能是因为随着输入地震波幅值的增加，边坡的整体刚度降低，塑性位移增加，模型对于能量的耗散作用逐渐增加。对比同一震级下EL波与Kobe波的模拟结果可知：加载EL波时，最大加速度响应发生在0.86倍高程处；而Kobe波的最大加速度响应发生在0.74倍高程处。且相同震级下，Kobe波的PGA最大值大于EL波的最大值。这是因为加速度响应与输入地震波的频谱特性有关，也与边坡的形状、自振频率及振动模态有关^[27]。由图9（b）可以看出，坡体内部PGA放大系数随高程增加而增加，符合“高程放大效应”。值得注意的是：强风化区域岩体的PGA响应强于弱风化区域；坡体内部PGA放大系数随输入地震波幅值的增加有所降低，且输入Kobe波的PGA响应比输入EL波更为剧烈。

3.3   FFT频谱

为分析边坡模型在不同地震波下的频率响应差异规律，对强风化岩体区域（M4、M6、M8）和弱风化岩体区域（M1、M2、M3）的加速度时程曲线做快速傅里叶变换，不同区域的FFT谱如图10和11所示，其中，f为各监测点傅里叶谱主频。由图10和11可知，不同类型和震级的地震波对边坡模型的傅里叶谱的分布规律有所影响。具体分析：在弱风化岩体区域，随着地震波幅值增加，输入不同类型地震波时，模型的主频f₁、f₂变化不大，f₃变化较为显著，表明弱风化区域M1和M2位置处的整体性较好，而M3处高程更大，在地震动作用下更易发生损伤，从而影响地震波的传递；与EL波相比，输入Kobe波时模型响应的FFT幅值更大。与弱风化岩体区域相比，强风化岩体区域FFT幅值整体偏小，这是因为随着地震波逐渐向上传递，至强风化岩体区域时，由于大量结构面的存在，导致地震波能量耗散增加，进而导致FFT幅值降低。对比M3和M4位置处的FFT幅值谱可知，地震波在经过强弱风化交界位置的结构面时会消耗部分能量，导致其FFT幅值降低，且该局部效应的影响大于高程效应。

图  10  弱风化岩体区域坡面监测点FFT谱

Fig.  10  FFT spectrum of slope monitoring points in weakly weathered rock mass region

下载: 全尺寸图片

图  11  强风化岩体区域坡面监测点FFT谱

Fig.  11  FFT spectrum of slope monitoring points in strongly weathered rock area

下载: 全尺寸图片

坡面主频所对应的FFT幅值随地震波输入震级（PAG）和高度的变化如图12所示。由图12可知：在加载Kobe波时，边坡模型的FFT幅值响应比加载EL波时更大；不同的地震波加载时，FFT的幅值均出现在测点M3处；FFT幅值总体呈现随着高程的增加先增大再减小，最后再增大的规律。

图  12  坡面FFT幅值随输入震级和高程变化

Fig.  12  Change of FFT amplitude of slope with input PGA and elevation

下载: 全尺寸图片

3.4   破坏模式

基于离散元UDEC的动力仿真分析过程中，可以清晰观察到模型的裂缝发展，进而全面掌握模型的失稳演化过程^[18]。为揭示边坡模型在地震动力作用下的失稳破坏模式，通过提取模型破坏发展的不同阶段，揭示其在地震作用下的失稳演化过程。由于不同工况条件下模型的破坏模式大致相同，选取具备代表性的工况1所对应的模型破坏现象进行具体分析。根据模型破坏模式特征的变化，将边坡破坏过程主要分为4个阶段，每个阶段均采取位移云图与位移矢量图相结合的方式进行分析，为清晰展现模型裂隙的发展，对部分坡体进行局部放大。加载0.6gEL波时，工况1模型的破坏过程如图13所示。

图  13  边坡失稳演化过程

Fig.  13  Instability evolution process of the slope

下载: 全尺寸图片

由图13（a）可知：在地震波作用下，坡表的强风化岩体区域产生初始震裂，且部分破碎块体向临空面抛出，形成落石；同时，强风化岩体区域下部节理部分张开，该区域块体向临空面产生滑移趋势。由图13（b）可知：随着地震波的持续输入，部分落石沿坡表向下加速滚动，强风化岩体区域节理裂缝迅速发展，同时下部弱风化岩体区域节理裂缝初步孕育。由图13（c）可知：强风化岩体区域的中上部节理裂缝逐渐变大，形成贯通结构面，且滑动面后缘部位产生裂缝，同时下部弱风化岩体节理裂缝进一步扩展。模型的最终破坏如图13（d）所示，强风化岩体区域中上部产生多组贯通结构面，并且整体产生明显的下滑趋势，落石在坡脚处形成堆积体。

综上所述，碎裂状顺层岩质边坡的破坏模式为裂隙拓展–强风化岩体滑移型破坏。可认为在地震作用下模型的破坏主要是由于节理裂缝的扩展和贯通导致上方强风化岩体区域向临空面产生整体性滑移失稳。在地震过程中，模型下方的弱风化岩体区域整体性保持较好，只是局部节理裂缝发展，未发现整体的大规模滑塌、崩塌破坏。因此，在治理破碎程度较大的边坡时，应重点关注岩体破碎程度较大的区域。

4.   岩质边坡稳定性影响因素敏感性分析

4.1   极差分析法

为揭示不同因素的敏感性，针对不同工况的计算结果采用极差分析法进行分析。首先，求出各个因素在不同水平上的计算结果的平均值；然后，对其求极差。当某列因素能够体现该列因素的水平变动时，试验结果产生相应变动，极差越大则该因素对试验结果影响越大，也即：极差越大，敏感度越高；反之，敏感度越低。极差计算公式为：

式中， $ {R}_{i} $ 为第i列因素的极差， ${\overline{K}_{ki}}$ 为第i列因素水平k的计算结果的平均值。

4.2   试验指标的选取

在评估边坡稳定性时，多数是基于其安全系数的变化。但由于地震作用下边坡的安全系数在不断地动态折减，离散元UDEC中虽然能通过命令流实现安全系数的求解，但是耗时较长。由于模型中每次输入的地震波加载时间相同，因而，与计算动态安全系数相比，通过采用边坡加载完成后的最终位移和PGA放大系数作为指标评估边坡的动力稳定性更为简便、有效。

4.3   敏感性结果分析

4.3.1   位移响应

选取输入0.4gEL波时坡表M4处的最终位移结果进行分析。不同工况下模型M4处的最终位移结果见表4，各因素极差分析结果见表5。

为了能够更加直观地表达A、B、C、D这4个参数的敏感性，利用柱状图表示A、B、C、D各因素对边坡最终位移的影响程度，结果如图14所示。由图14可以看出，极差大小为 ${R}_{A} > {R}_{B} > {R}_{D} > {R}_{C}$ 。该结果表明，对此边坡位移影响最大的因素为弱风化岩体区域的层理距离，其次依次为弱风化岩体区域的次级节理间距、强风化岩体区域垂直节理间距、强风化岩体区域层理间距。

图  14  影响因素对位移敏感性排序

Fig.  14  Sensitivity ranking of influencing factors to displacement

下载: 全尺寸图片

4.3.2   加速度响应

根据加速度响应分析结果，选取坡表加速度响应最大的M7处的PGA放大系数分析。不同工况下模型M7处的PGA放大系数见表6，极差分析结果如图15所示。

图  15  影响因素对PGA放大系数敏感性排序

Fig.  15  Sensitivity ranking of influencing factors to PGA amplification factor

下载: 全尺寸图片

由图15可知，各个因素对PGA放大系数的敏感性排序为 $ {R}_{D} > {R}_{A} > {R}_{C} > {R}_{B} $ 。该结果表明，对PGA放大系数影响最大的因素为强风化岩体区域的次级节理间距，其次依次为弱风化岩体区域的层理距离、强风化岩体区域层面间距、弱风化岩体区域的层理距离。

5.   结　论

基于离散元UDEC，研究了碎裂状顺层岩质边坡在地震作用下的动力响应与破坏模式；揭示了不同地震波加载下，边坡模型的PGA放大系数、位移及FFT频谱的响应规律与边坡破坏过程；基于正交试验法，研究了不同节理参数的敏感性。主要得到以下结论：

1）不同地震荷载作用下，坡面位移规律大致相同，从坡底至坡顶呈先增加再减小的规律，地震波幅值越高，坡面的最终位移越大，坡面位移响应最大处位于M4处。在震级相同的情况下，与EL波相比，Kobe波的坡面位移响应更剧烈。

2）边坡模型的PGA放大系数明显呈现“高程放大效应”。地震波的类型对PGA有一定的影响，输入EL波时，PGA响应的最大位置在M7处；输入Kobe波时，PGA响应的最大位置在M6处。并且，地震波幅值相同时，Kobe波的PGA响应比EL波的更剧烈。

3）地震波类型和震级影响FFT频谱分布规律，总体上FFT幅值随着高程的增加呈先增大再减小，再增加的规律；强风化岩体区域的FFT幅值整体偏小；地震波在通过强弱风化交界位置的结构面时，会导致局部FFT幅值降低。

4）与EL波相比，加载Kobe波时，边坡模型的FFT幅值更大。不同地震波作用下FFT幅值的最大位置均位于模型的M3处。

5）碎裂状边坡的破坏过程主要为：坡表震裂→形成落石→强风化岩体下部节理张开→弱风化岩体节理孕育→裂缝发展→强风化岩体区域形成贯通结构面→强风化岩体区域沿结构面下滑失稳破坏→坡脚形成堆积体。

6）采用正交试验法和极差分析法所得的节理参数对位移和加速度响应的敏感性不同，针对位移的敏感性顺序为A>D>B>C，针对加速度响应的敏感性顺序为D>A>C>B。

由此可见，与传统的岩质边坡相比，碎裂状岩质边坡在加速度、位移及FFT频谱等方面的响应均有所区别。本文仅讨论了节理参数对碎裂状岩质边坡的敏感性，对于其他因素的敏感性应在今后继续研究。