引言

大量震害调查、强震观测及理论研究表明（Borcherdt，1970；Boore，1972；Ashford等，1997），地形条件对地震动的影响较大，不容忽视。通过整理震害出现的地形资料，依据地貌几何形态将地形分为河谷、山地、阶地、斜坡、陡坎，山脊、山梁，台塬、盆地。大量学者研究了不同地形对地震动特性的影响，如梁建文等（2003）给出了凹陷地形SV波入射下的解析解；巴振宁等（2018）采用间接边界元法研究了SV波入射下沉积谷地对地震动的放大效应，并分析了有无覆盖层及多个谷地对结果的影响；刘中宪等（2019）基于快递多极子间接边界元法，针对城市区域沉积盆地进行了三维地震动响应模拟，结果表明入射频率、盆地形状和盆地内部位置点对沉积盆地地表位移放大作用均有不同程度的影响；周国良等（2012）和孙纬宇等（2019）研究了河谷地形对地震动分布特性的影响，结果表明河谷地形对地震动分布具有重要影响。沿河谷斜坡顶到谷底的各观测点地震动峰值加速度存在显著差异，地震动位移存在差动现象。袁晓铭等（1996）给出了任意圆弧形凸起地形对平面SH波散射的解析解；蒋涵等（2015）采用三维谱元法模拟地震波传播，研究了芦山地区地震动频谱特征；李英民等（2019）和郝明辉等（2019）分别采用动态子结构法、有限元有限差分法研究了凸起地形对地震动特性的影响，并分析了坡角、坡高、入射角度等因素对计算结果的影响。近年来，大量城市建筑建于山坡上，坡地地形对地震动的放大效应不容忽视（殷跃平，2008）。《建筑抗震设计规范》（GB 50011—2010）（中华人民共和国住房和城乡建设部等，2010）针对坡地地形给出了设计地震动影响系数的调整；李英民等（2010）采用脉冲输入研究了岩质坡地地形对反应谱特性的影响；李家祥等（2021）研究了SH波入射下斜坡地形对地震动特性的影响，定量分析了SH波垂直入射下二维斜坡地形地震动响应与斜坡角度、土层厚度及介质阻抗比的关系。但上述研究均采用脉冲输入，而实际地震动是非常复杂的，且受地震波频散和干涉等因素影响，实际地震动脉冲之间存在复杂的相位关系。因此，基于脉冲输入的地形对地震动峰值的放大倍数无法直接用于评估地形对实际地震动峰值的放大效应（郝明辉等，2014）。此外，除峰值、频谱特性外，阿里亚斯烈度也是表征地震滑坡稳定性的重要地震动特性参数（刘甲美等，2015）。因此，本文针对坡地地形，基于透射人工边界的有限元有限差分方法，采用人工合成地震动，研究SV波垂直入射下，坡地地形对地震动峰值、反应谱、阿里亚斯烈度等的放大效应，并定量分析坡高、坡角对计算结果的影响，从而为坡地地形上建筑物选址及地震动参数的确定提供理论参考。

1.   模型及分析方法

局部坡地地形模型如图1所示。坡地地表包括顶部平台、斜坡面和坡底段。利用人工边界从无限大场地中切取包含局部坡起在内的有限元计算区域，脚点与左、右人工边界的垂直距离取为5H/tanα（H为坡高，α为坡角），与底人工边界的垂直距离为5H，可满足计算精度要求，减小人工边界的影响。采用四边形有限元网格离散，网格尺寸为1 m×1 m，以获取对高频成分可靠的数值模拟结果。介质假定为均匀、各向同性、黏弹性，纵波波速表示为c_p，横波波速表示为c_s，阻尼比表示为ζ。本文主要考虑岩石介质，其剪切波速取为1400 m/s，压缩波速取为2 000 m/s，人工边界区介质波速取为1200 m/s，仅考虑平面SV波垂直入射情况。

图 1  局部坡地模型

Figure 1.  A model with local terrain

下载: 全尺寸图片 幻灯片

为消除输入地震动随机性对分析结果的影响，本文采用多组输入地震动计算得到的反应谱谱比平均值，研究坡地地形对地震动反应谱的放大效应。针对第$ i $个输入地震动样本$ {a_i}（t） $，利用程序TPG2D求解图1所示凸起的地震反应，得到不同输出点地震反应加速度时间过程$ a_i^{\text{r}}（t,\vec x） $，其绝对加速度反应谱为$ S_{{\text{a}},i}^{\text{r}}（T,\vec x）$，$ \vec x $为输出点空间坐标，$ \vec x = {[x,y]^{\text{T}}} $。

程序TPG2D假定$ {a_i}（t） $为图1所示模型左下角点处入射波加速度，以此为基础，可确定自由场（即不存在凸起的弹性半空间）地震反应，其加速度反应记为$ a_i^{\text{f}}（t,\vec x） $，反应谱记为$ S_{{\text{a}},i}^{\text{f}}（T,\vec x） $。相应地，定义式（1）所示反应谱谱比用于描述局部凸起对地震动的放大效应：

已有研究表明（郝明辉等，2014），利用Hanning窗对3个输入地震动样本计算所得平均反应谱谱比曲线进行平滑，所得结果与大样本数量所得平均值相近，可有效消除输入地震动随机性的影响，同时降低计算工作量：

式中，$ {\text{Hann}}[f\left( x \right)] $表示利用Hanning窗对函数$ f\left( x \right) $进行平滑处理。下文研究的反应谱谱比是由式（2）定义的反应谱谱比。后续计算结果均采用3个输入地震动样本计算所得结果的平均值表征。

基于郝明辉等（2014）的研究，采用人工合成的地震动时程作为模型输入，研究凸起地形对地震动反应谱的放大效应，其中1个输入地震动样本加速度、速度和位移时程曲线如图2所示。

图 2  输入地震动时程曲线

Figure 2.  Acceleration, velocity and displacement curves of artificial ground motion samples

下载: 全尺寸图片 幻灯片

2.   计算方法

2.1   集中质量显式有限元内点计算方法

集中质量显式有限元的实质是从当前时刻节点运动方程推求下一时刻节点运动，无须进行刚度矩阵、质量矩阵、阻尼矩阵集成，其右端项的形成仅需在单元一级水平上根据每个单元对有效荷载向量的贡献累加，整个计算基本在单元一级水平上进行，仅需很小的高速存贮区，计算效率较高。当一系列单元刚度矩阵、质量矩阵、阻尼矩阵相同时，无须重复计算，效率更高。

对于模型内部节点（图1中2～7号点，假定其总体编号为1），由动力平衡条件可建立运动方程为：

式中，${N_e}$为包含节点1的所有单元个数，${N_n}$为第$e$个单元节点数，$i$为节点1在单元$e$中的局部编号，$ \hat{n} $为单元$e$中第$n$个节点总体编号，$M_{in}^{\left( e \right)}$和$K_{in}^{\left( e \right)}$分别为单元$e$质量矩阵和刚度矩阵，$\alpha $和$\beta $为Rayleigh阻尼系数，$\left\{ {{{\ddot u}_j}\left( t \right)} \right\}$、$\left\{ {{{\dot u}_j}\left( t \right)} \right\}$和$\left\{ {{u_j}\left( t \right)} \right\}$分别为加速度、速度和位移，$\{ {R_1}\left( t \right)\} $为节点等效集中荷载，${M_1}$为节点1等效集中质量，且有：

利用李小军等（1992）提出的显式差分格式求解方程（1），可得到节点1运动数值解为：

式中，$\left\{ {{\dot{{u}}}_{{i}}^{{p}}} \right\}$、$\left\{ {{{u}}_{{i}}^{{p}}} \right\}$和$\left\{ {{{R}}_{{i}}^{{p}}} \right\}$分别为$ {{p}}\Delta {\text{t}} $时刻节点$i$速度、位移和等效集中动力荷载，$ \Delta {\text{t}} $为离散时间间隔。

利用式（5）可显式求解内部节点动力响应。

2.2   透射人工边界

对于人工边界节点（图1节点1、7，假定其总体编号为0），基于透射人工边界原理，可得动力响应计算公式为：

式中，$ {{N}} $为透射阶数（本文取为2阶），$\left\{ {{{u}}_0^{{{p}} + 1}} \right\}$为人工边界节点在$({{p}} + 1)\Delta {{t}}$时刻的位移，$\left\{ {{{u}}_{{j}}^{{{p}} + 1 - {{j}}}} \right\}$为计算点${{x}} = - {{j}}{{{c}}_{\text{a}}}\Delta {{t}}$在$({{p}} + 1 - {{j}})\Delta {{t}}$时刻的位移，${{{c}}_{\text{a}}}$为人工波速值，${{C}}_{{j}}^{{N}}$为二项式系数：

2.3   程序验证

采用上述基于透射人工边界的显式有限元计算方法，开发了并行化Fortran数值计算程序TPG2D，并移植到北京工业云计算平台上运行。为验证该程序计算结果的正确性，计算了弹性半空间在SV波和P波垂直入射下的动力响应，并与理论解进行对比。

从二维半无限空间中截取$ 6\;{\text{m}} \times 50\;{\text{m}} $的有限范围，单元网格尺寸为$ 1\;{\text{m}} \times 1\;{\text{m}} $，材料弹性模量$ E = $$ 2.4 \times {10^7}\;{\text{Pa}} $，泊松比$ \mu = 0.2 $，质量密度$ \rho = 1\;000\;{\text{kg}}/{{\text{m}}^3} $。在底部垂直向上入射水平方向单位脉冲剪切位移波和竖直方向单位脉冲压缩位移波：

式中，$ f = 4\;{\text{Hz}} $。

图3所示为程序TPG2D计算得到的有限元模型底部、中部、顶部水平向和竖直向位移时程曲线与理论解的比较。入射波由底部向上传播，在自由地表发生反射，自由地表处位移幅值为输入波幅值的2倍。由图3可知，程序TPG2D计算结果与理论解吻合较好，证明了本文采用计算方法的正确性和有效性。

图 3  数值解与理论解比较

Figure 3.  Comparison of viscous-elastic boundary calculated results with the theoretical solution

下载: 全尺寸图片 幻灯片

3.   坡地地形地表不同位置放大系数

以坡高50 m、坡角45°的坡地地形为例，分析地形效应对地震动反应谱特性的影响。坡底段、斜坡段、平台段不同空间点反应谱谱比曲线如图4～6所示。由图4～6可知，在远离斜坡的坡底段反应谱谱比基本在1附近波动，越靠近斜坡放大效应越明显（节点2）；对于斜坡段，各点反应谱谱比最大值沿坡高逐渐增大，而坡脚点对地震动反应谱呈缩小效应（蓝色实线），斜坡顶点在大部分频段对地震动起放大作用，反应谱谱比最大值达1.24（红色实线）；对于斜坡平台段，平台节点在大部分频段反应谱谱比均>1，对地震动均起放大作用。

图 4  坡底段地表点反应谱谱比曲线

Figure 4.  Spectral ratio chart of surface point on bottom

下载: 全尺寸图片 幻灯片

图 5  斜坡段地表点反应谱谱比曲线

Figure 5.  Spectral ratio chart of surface point on the slope

下载: 全尺寸图片 幻灯片

图 6  平台段地表点反应谱谱比曲线

Figure 6.  Spectral ratio chart of the surface point on the platform

下载: 全尺寸图片 幻灯片

4.   坡高对计算结果的影响

坡角α取为45°，分析坡高变化对斜坡顶点反应谱谱比的影响。入射地震波采用垂直入射的SV波。参与计算的凸起地形坡高H分别取20.0、30.0、40.0、50.0、60.0、70.0、80.0、90.0、100.0、110.0、130.0、150.0、170.0、190.0、210.0 m。斜坡顶点反应谱谱比随坡高变化曲线如图7所示，由图7可知，不同坡高反应谱谱比曲线形状基本相同。反应谱谱比峰值点对应的特征周期随着坡高的增加呈增长趋势，峰值点对应的反应谱谱比随着坡高的增加呈增大趋势。在本文研究的坡高范围内，周期>1.0 s的反应谱谱比随着坡高的增加呈增大趋势，最大值达1.3，且随着周期的增大，反应谱谱比逐渐趋于1.0，说明这一坡高范围对较长周期地震动成分的影响逐渐减弱，在反应谱上表现为斜坡顶点地震动反应谱与自由地表地震动反应谱相差较小。周期为0.1～0.3 s时，随着坡高的增加，反应谱谱比呈减小趋势。在其他周期范围内，反应谱谱比随坡高的变化不再呈单调变化趋势，变化规律较复杂。

图 7  斜坡顶点反应谱谱比随坡高变化曲线

Figure 7.  Spectral ratio curve of the vertex with different heights

下载: 全尺寸图片 幻灯片

图8所示为坡脚点反应谱谱比随坡高变化的曲线，由图8可知，不同坡高反应谱谱比曲线在大部分周期点对地震动的放大作用均不明显，基本呈减弱趋势，反应谱谱比最小值为0.5。在本文研究的坡高范围内，所有周期范围内反应谱谱比随坡高的变化不再呈单调变化趋势，变化规律较复杂。

图 8  斜坡脚点反应谱谱比随坡高变化曲线

Figure 8.  Spectral ratio curve of the Slope toe with different heights

下载: 全尺寸图片 幻灯片

图9所示为斜坡顶点不同反应谱控制周期点处反应谱谱比随坡高的变化曲线，由图9可知，长周期范围（蓝线）反应谱谱比随坡高变化较明显，且反应谱谱比随坡高的增加呈增大趋势，最大值达1.25左右；当坡高>160 m后，反应谱谱比随坡高的增加呈减小趋势；短周期范围（黑线）随坡高变化反应谱谱比基本呈减弱趋势，当坡高>60 m后，峰值加速度基本随着坡高的增加而减弱，且当坡高>130 m后峰值加速度放大系数基本保持不变；峰值速度随着坡高的增加呈增大趋势，PGV放大倍数最大值达1.2左右，PGD放大倍数基本在1左右。

图 9  不同周期点处斜坡顶点反应谱谱比随坡高变化曲线

Figure 9.  Variation curve of slope peak spectral ratio with height at different periodic points

下载: 全尺寸图片 幻灯片

图10所示为斜坡顶点阿里亚斯烈度比值随坡高变化曲线，由图10可知，对于宽度为50 m、坡角为45°的计算模型，当坡高<30 m时，斜坡顶点阿里亚斯烈度比值随着坡高的增加呈增大趋势；当坡高>30 m时，斜坡顶点阿里亚斯烈度比值随着坡高的增加呈减小趋势。

图 10  斜坡顶点阿里亚斯烈度比值随坡高变化曲线

Figure 10.  Curve of arias intensity ratio with height at the top of slope

下载: 全尺寸图片 幻灯片

图11所示为节点1～7（坡高50 m，坡角45°）阿里亚斯烈度比值，由图11可知，地形对平台中点和平台顶点阿里亚斯烈度的放大作用较大，最大值达1.2；地形对脚点、斜坡中点阿里亚斯烈度起减小作用；地形对自由地表点阿里亚斯烈度的放大系数基本为1，这与地表点反应谱谱比放大规律基本一致。

图 11  地表点阿里亚斯烈度比值

Figure 11.  Arias intensity ratio map of surface points

下载: 全尺寸图片 幻灯片

5.   坡角对地震动特性的影响

将图1所示凸起地形模型宽度B取为50 m，坡高H取为50 m，研究坡角变化对地震动反应谱特性及阿里亚斯烈度的影响。入射地震波为垂直入射的SV波，参与计算的坡角分别取20°、25°、30°、35°、40°、45°、50°、55°、60°。

不同坡角对应的斜坡顶点反应谱谱比曲线如图12所示，由图12可知，相比凸起台地地形，不同坡角反应谱谱比曲线变化规律不同（郝明辉等，2019），反应谱谱比最大值为1.28左右，说明坡角变化对反应谱某频段上的最大反应的影响较大。在坡高不变的情况下，随着坡角的增大，高频段（0.01～0.03 s）及低频段（0.3～6.0 s）放大倍数增大，中间频段（0.03～0.3 s）在不同周期点处放大倍数与坡角之间的变化规律较复杂。

图 12  不同坡角对应的斜坡顶点的反应谱谱比曲线

Figure 12.  Spectral ratio curve of the vertex with different slopes

下载: 全尺寸图片 幻灯片

图13所示为不同反应谱控制周期点处反应谱谱比随坡角变化曲线，由图13可知，长周期范围（蓝线）反应谱谱比随着坡角的增加呈增大趋势，反应谱谱比最大值为1.1左右；短周期范围（黑线）反应谱谱比随着坡角的变化较复杂，在坡角较大的情况下，反应谱谱比随着坡角的增加呈增大趋势；PGA受坡角的影响较大，反应谱谱比最大值达1.3左右；PGV、PGD随坡角变化不明显，PGV放大倍数基本为1.1左右，PGD放大倍数基本为1左右。

图 13  不同周期点处反应谱谱比随坡角变化曲线

Figure 13.  Variation curve of slope peak spectral ratio with slope at different periodic points

下载: 全尺寸图片 幻灯片

图14所示为斜坡顶点阿里亚斯烈度随坡角变化曲线，由图14可知，对于宽度50 m、坡角45°计算模型，斜坡顶点阿里亚斯烈度比值在坡角>30°后基本随着坡角的增加呈增大趋势，最大值达1.38。

图 14  斜坡顶点阿里亚斯烈度比值随坡角变化曲线

Figure 14.  Curve of arias intensity ratio with slope at the top of slope

下载: 全尺寸图片 幻灯片

6.   结论

本文基于有限元有限差分法，结合人工透射边界理论，研究了局部坡地地形对地震动特性的影响，分析了坡高、坡角对地形放大效应的影响，得出以下结论：

（1）斜坡顶点地震动受地形效应的影响较大，在零阻尼条件下，反应谱谱比曲线呈双峰特点，脚点在大部分周期点对地震动呈缩小效应。

（2）在坡角一定的情况下，斜坡顶点反应谱谱比峰值点对应的特征周期随着坡高的增加呈增长趋势，对应的谱比随着坡高的增加呈增大趋势。斜坡顶点阿里亚斯烈度比值在斜坡坡高大于一定值后随着坡高的增加呈减小趋势，与地震动高频成分受地形影响的规律基本一致。

（3）当坡高不变时，斜坡顶点反应谱谱比放大倍数与坡角之间的变化规律较复杂。此外，斜坡顶点阿里亚斯烈度比值在坡角大于一定值后基本随着坡角的增加呈增大趋势。