2017, Vol. 60
2. 滨海土木工程结构与安全教育部重点实验室, 天津 300072;
3. 中国人民武装警察部队后勤学院建筑工程系, 天津 300304
2. Key Laboratory of Coastal Structures in Civil Engineering and Safety of Ministry of Education, Tianjin 300072, China;
3. Department of Civil Engineering, Logistics College of Chinese People's Armed Police Force, Tianjin 300304, China
多次的震害调查和地震观测均表明，凸起地形对地震波的传播有着显著的影响，凸起表面的地震动会被显著放大.如1985年Central Chile地震导致了大量的结构破坏，其中山脊顶部建筑破坏尤为严重，有的几乎不可修复 (Çelebi, 1987)；1994年美国北岭地震中Tarzana山顶记录到的地震加速度高达1.78 g，其中部分原因来自于山体地震效应 (Spudich et al., 1996).鉴于许多居民建筑及重要特殊结构 (通讯及输电塔等) 建造于山体等凸起地形之上，研究凸起地形对地震动的放大作用在理论和工程上都有着重要的意义.
为探索凸起地形对地震动的影响机理，诸多学者针对该问题进行了理论研究.解析方法方面：Yuan和Men (1992)首次采用波函数展开法给出了均匀半空间中半圆凸起对SH波的散射解析解；此后，袁晓铭和廖振鹏 (1996)又进一步提出辅助函数方法求解了浅圆弧凸起对SH波的散射问题；Lee等 (2006)采用Cosine Half-Range展开方法给出了半圆凸起对SH波的散射解析解；Tsaur (2011)采用Mathieu函数展开法给出了矮半椭圆凸起对SH波散射解析解；Liang和Fu (2011)采用正弦和余弦级数展开方法研究了半椭圆凸起对SH波散射问题，文中首次给出了扁长凸起高频散射结果；梁建文等 (2006)采用波函数展开法，并借助辅助函数思想，给出了平面SV波入射下半圆凸起散射解析解.解析解的波函数需严格满足波动方程和边界条件，因此研究限于均匀半空间中简单形状凸起对弹性波的散射.数值方法方面：Boore (1972)和Simth (1975)分别采用有限差分法和有限元法求解了三角凸起对SH波的散射；Bouchon (1973)采用Aki-Larner方法研究了余弦型凸起对SH波的散射问题；Sills (1978)、Bard (1982)和Sánchez-Sesma等 (1982)则分别采用积分方程法、离散波数法和边界元方法研究了高斯型凸起对SH波的散射；Fu (2005)也分别采用几种近似理论研究了高斯型凸起对SH波的散射；梁建文和巴振宁 (2008)采用间接边界元方法研究了层状半空间中凸起地形对SH波的散射；Pedersen等 (1994)、Takenaka等 (1996)和Zhang等 (1998)分别采用间接边界元方法、离散波数法和边界元-有限元耦合方法研究了二维凸起地形对弹性波的三维散射；Sánchez-Sesma (1983)采用边界方法研究了三维凸起对垂直入射P波的散射；Bouchon等 (1996)以及Sohrabi-Bidar和Kamalian (2013)分别采用离散波数和边界元方法研究了三维高斯型山体对垂直入射P和SV波的散射.近年来，高性能计算机的出现，使得数值方法得到了极大发展.
值得指出的是，以上文献仅以单一凸起地形为模型进行了研究，然而凸起 (山体) 多连绵分布，研究多个凸起对地震动的影响更为符合实际情况.目前，还鲜有针对多个凸起的研究，据作者所知，刘晶波 (1996)采用有限元法研究了三个顺次排列的山梁对P、SV和Rayleigh波的二维散射，研究表明多个山体的存在使得地面运动进一步放大，持时延长；Liu等 (2010)研究了三角凸起与半圆凸起相连地形对SH波的散射问题；韩峰等 (2013)利用复变函数和多极坐标方法，研究了两个等腰三角凸起与半圆凹陷相连地形对SH波的散射问题；杜永军等 (2009)采用“分区和契合”技术，研究了两个等腰三角凸起对SH波的放大效应；Yang等 (2014)采用波函数展开和复变函数法研究了两个不等边三角形凸起与半圆凹陷相连地形对SH波的散射问题.针对任意多个 (无限多个) 凸起情况的精确求解在数学处理上则仍有一定的困难，若直接针对所研究区域，截取一定范围内的凸起进行求解，而放松其他凸起的边界条件，会因截断边界引入误差，同时计算量和存储量也会非常大，有时甚至难以求解 (如较低的频率、较小的凸起间距和较低的介质阻尼).为此，本文在梁建文和巴振宁 (2008)研究的基础上，提出一种新的以层状半空间中周期分布斜线荷载动力格林函数为基本解的间接边界元方法，并具体求解了周期分布凸起地形对平面SH波的散射问题.研究的意义在于提出一种用于求解周期局部地形对平面SH波散射的间接边界元方法，同时阐明任意多个 (周期) 凸起间的地震相互作用.文中对提出的方法进行了介绍，并以均匀半空间和基岩上单一土层中周期分布凸起为例，进行了数值计算分析，讨论了凸起间距、凸起截面形状、凸起介质材料、土层、波入射角度和入射频率对动力响应的影响，得到了一些结论.
2 模型与计算方法模型如图 1所示，周期分布任意截面形状凸起地形 (标号为凸起“-∞”、…、“-2”、“-1”、“0”、“1”、“2”、…、“∞”) 位于层状半空间中，各凸起自由表面分别为S′-∞、…、S′-2、S′-1、S′0、S′1、S′2、…、S′∞，各凸起与半空间交界面分别为S-∞、…、S-2、S-1、S0、S1、S2、…、S∞，凸起沿x轴等间距周期分布，各凸起截面形状保持不变，凸起间距离为L.层状半空间由N层水平土层和其下的基岩半空间组成，各水平土层的质量密度为ρjL，复剪切波速为cs, jL*=cs, jL[1+2i sgn (ω)ζjL]，阻尼比为ζjL(j=1~N)，基岩半空间的质量密度为ρR，复剪切波速为csR*=csR[1+2i sgn (ω)ζR]，阻尼比为ζR，凸起介质的质量密度为ρH，复剪切波速为csH*=csH[1+2i sgn (ω)ζH]，阻尼比为ζH.入射波为平面SH波，入射圆频率为ω，入射方向与x轴的夹角为ψ.
|
图 1 层状半空间中周期分布凸起地形平面SH波入射模型 Fig. 1 Model of scattering of plane SH-waves by periodic hills in a layered half-space |
图 1所示的层状半空间中周期分布凸起对SH波的散射问题，属无限边界散射问题.若直接针对所研究的区域，截取一定范围内的凸起地形进行单元的离散和求解，会因放松其他凸起地形的边界条件而引入误差，同时计算量和存储量也较大，有时甚至无法完成求解.但值得注意的是，图 1所示周期分布凸起地形在平面SH波入射下，各凸起地形周围波场幅值完全相同，仅因所在x轴的位置不同而在时域内存在一个时间差 (在频域内存在一个相位差)，因此只要其中任意一个凸起地形的边界条件满足，其余凸起地形的边界条件也自动满足.利用这一特征，我们可仅取一个凸起进行边界单元离散和求解，即可求得问题的解.
具体求解时，选取标号为“0”的凸起进行求解.然后针对求解区域，如图 2所示，采用Yuan和Men (1992)提出的“分区契合”方法将求解区域分解为凸起闭合域ΩH和凸起外部开口域ΩL.假定ΩH内仅存在散射波场，而ΩL内存在自由波场和散射波场.ΩL内的散射波场是所有不规则边界引起散射波场之和，如上所述，考虑到各不规则边界产生的散射波场完全相同，在频域内仅因位置的不同而存在一个相位差，因而本文提出了层状半空间中周期分布斜线荷载动力格林函数，进而通过在凸起地形与半空间交界面S0上的各离散斜线单元上施加虚拟周期分布荷载产生的动力响应来模拟.ΩL的内自由波场则可通过“直接刚度法”求得.ΩH内的散射波场完全独立，不受其他凸起闭合域的影响，因而本文采用Wolf (1985)提出的斜线荷载动力格林函数 (非周期)，通过在边界S0和S′0上的各离散斜线单元上施加虚拟均布荷载产生动力响应来模拟.
|
图 2 求解模型分区示意图 Fig. 2 Decompose the model into a closed hill region and an opened layered half-space region |
本文采用Wolf (1985)给出的“直接刚度法”求解自由波场.所采用的土层和半空间的平面外精确动力刚度矩阵分别如下式所示：
|
(1) |
|
(2) |
其中，
集整各土层刚度矩阵
|
(3) |
其中，{Q1, Q2, …QN}T为作用在各土层交界面上的外荷载幅值向量，{v1, v2, …vN}T为各土层交界面上的位移幅值向量.
对于由基岩面入射的平面SH波，若将控制点选在基岩露头，则外荷载幅值向量{Q1, Q2, …QN}T中，其他元素均为零，仅有
|
(4) |
其中，v0为基岩露头位移幅值，可取为v0=1.0.将式 (4) 代入式 (3) 可求得各土层交界面处位移幅值，进而可求得各土层内上行和下行波的幅值系数，最后可求得各土层内任意位置的位移和应力幅值，也即自由波场.以上计算均在频域中进行.关于自由场更为详细的求解内容可参考Wolf (1985).
2.2 格林函数及散射波场构造本文将散射波场分解为凸起闭合域ΩH和凸起外部开口域ΩL内的散射波场.对于ΩH内的散射波场，考虑到其独立性 (各凸起间无相互影响)，采用在凸起边界S0和S′0各离散单元上施加均布荷载产生的动力响应来模拟 (梁建文和巴振宁, 2008；Wolf, 1985).对于ΩL内的散射波场，虽然各凸起地形间相互影响，但考虑到层状半空间中周期分布凸起地形在平面SH波入射下，每个凸起周围的散射波场完全相同，仅在时域内存在一个时间差 (频域内存在一个相位差).依据此规律，本文提出了层状半空间中周期分布斜线荷载动力格林影响函数，采用在凸起与层状半空间交界面S0各离散单元上施加周期分布荷载产生的动力响应来模拟.
如图 3所示，层状半空间中周期分布斜线荷载可表示为
|
(5) |
其中，qe是等效均布荷载，qn(x, z, t)为标号“n”的斜线荷载 (出平面，方向沿y轴方向)，δ为狄拉克函数，θ为斜线与水平方向的夹角，L表示各斜线荷载间的距离，z表示距荷载上边界的距离.特别的对于标号为“n”和“0”的斜线荷载可表示为
|
(6) |
|
(7) |
|
图 3 层状半空间中周期分布斜线荷载动力函数示意图 Fig. 3 Periodically distributed loads acting on inclined lines over part of a layer |
由图 1以及式 (6) 和 (7) 可得
|
(8) |
其中，nL/ca可认为是第“n”个斜线荷载与第“0”个斜线荷载产生动力响应传播的时间差.由图 1，ca=csR/cosψ为SH波沿x轴传播的视速度.将式 (8) 代入式 (5) 得
|
(9) |
为了求解方便，引入傅里叶变换及其逆变换
|
(10) |
|
(11) |
采用傅里叶变换方法推导周期分布斜线荷载的动力格林函数，对式 (9) 进行关于t和x的傅里叶变换得
|
(12) |
其中，q0(k, z, ω)为第“0”个斜线荷载在频率和波数域内的幅值，k′=ω/ca.设f0(k, z, ω) 为第“0”个斜线荷载在频率和波数域中产生的动力响应，则由式 (12) 可得第“n”个斜线荷载在频率和波数域中产生的动力响应为
|
(13) |
由式 (13) 可得所有斜线荷载 (周期分布斜线荷载) 产生的动力响应为
|
(14) |
以上计算都是在频率和波数域内进行的，采用傅里叶逆变换可求得频率和空间域内的动力响应：
|
(15) |
由式 (15) 可知，只要求得f0(k, z, ω)，便可通过式 (15) 求得周期分布斜线荷载动力格林函数.
本文采用Wolf (1985)给出的方法求解f0(k, z, ω).其主要思路为：首先假定作用荷载q (k, z, ω)的土层固定，求得层内反应和固定断面反力；然后将反力以相反方向作用到总体系上，求得反力产生动力响应；最后叠加层内反应和固定端面解，求得总的响应.以下简要介绍f0(k, z, ω) 的求解过程.
以v (k, z, ω)eiωt-ikx表示的动力平衡方程为
|
(16) |
其中，μ*=μ(1+2iζ) 是复剪切模量，ζ是介质阻尼比，ρ是质量密度.方程 (16) 的位移特解 (以上标“p”表示) 为
|
(17) |
其中，qy0为均布荷载密度，
|
(18) |
为使作用荷载土层上下端面固定，特解还必须加上与负的vj1p和vj2p相应的齐次解 (以上标“h”表示).固定土层内位移齐解为
|
(19) |
外荷载 (齐解固定端面反力) 可采用层刚度矩阵SSH, jL以及vj1h和vj2h求得
|
(20) |
叠加特解固定端面反力和齐解固定端面反力并取负号，可得作用到整个层状半空间上的外荷载幅值
|
(21) |
将式 (21) 代入式 (3)，采用直接刚度法可求得固定端面反力产生的动力响应.然后叠加固定端面反力动力响应和固定荷载作用层内动力响应，可求得波数域内总响应f0(k, z, ω).其中，层内动力响应包括式 (17) 给出的特解和式 (19) 给出的齐解.将f0(k, z, ω) 代入 (14) 便可求得频域内周期分布斜线荷载动力格林函数.值得指出的是，式 (14) 仅考虑第“0”项求得的结果，即为频域内层状半空间中斜线荷载 (单一斜线荷载) 动力格林函数.
求得周期分布斜线荷载动力格林函数以及斜线荷载动力格林函数后，则凸起闭合区域ΩH内的位移和牵引力以及凸起外部开口域ΩL中的位移和牵引力可分别表示为
|
(22) |
|
(23) |
|
(24) |
|
(25) |
其中，[g1vH(x, z, ω)]和[g1tH(x, z, ω)]分别为域ΩH内位移和沿y方向牵引力斜线荷载动力格林影响函数.[g2vL(x, z, ω)]和[g2tL(x, z, ω)]分别为域ΩL内位移和沿y方向牵引力周期分布斜线荷载动力格林影响函数.{p1}T为求解ΩH内散射波场而施加边界S0和S′0上的均布斜线荷载向量，{p2}T为求解ΩL内散射波场而施加在交界面S0上的周期分布斜线荷载向量.上标“H”则表示与ΩH相应的位移和牵引力，上标“L”表示与ΩL相应的位移和牵引力，下标“1”代表斜线荷载动力格林函数，下标“2”代表周期分布斜线荷载动力格林函数.
2.3 边界条件及问题求解以标号为“0”的凸起地形为研究对象，进行边界单元的离散和求解.标号为“0”的凸起地形自由表面S′0上零牵引力条件，交界面S0上的位移和牵引力连续条件可表示为
|
(26) |
|
(27) |
|
(28) |
其中，vfL(s0) 和tfL(s0) 为交界面S0上的自由场位移和牵引力，可由式 (3) 求得.[W (s)]为权函数，由梁建文和巴振宁 (2008)，可取为单位矩阵，使零牵引力条件在每个边界单元上独立满足，将式 (22)、(23)、(24) 和 (25) 代入式 (26)、(27) 和 (28) 得
|
(29) |
|
(30) |
|
(31) |
其中，
|
(32) |
|
(33a) |
|
(33b) |
|
(34) |
|
(35) |
由式 (29)、(30) 和 (31) 可求得荷载向量{p1}和{p2}，结合式 (23) 和式 (25)，可得凸起内外的位移幅值为
|
(36) |
|
(37) |
本文给出了一种用于求解层状半空间中周期凸起地形对SH波散射问题的新方法，方法仅需针对一个凸起地形进行边界单元的离散和求解，便可求得问题的解.该方法的精度取决于两个方面：凸起地形边界单元的离散和式 (15) 给出周期分布斜线荷载动力格林函数的精度.其中，式 (15) 给出格林函数的精度又取决于波数域到空间域积分的精度和求和公式的精度.
对于凸起地形边界单元的离散，经验算表明，每波长采用12个单元离散，即可达到满意的离散精度，本文在后面实际计算中采用每波长16个斜线单元离散来保证边界离散的精度.对于式 (15) 中波数域到空间域的积分，本文采用分段高斯积分来完成，要保证其精度必须采用足够大的波数积分区间和足够小的波数积分间隔.最大积分波数与入射频率和数值nL有关，针对不同的频率和不同nL，本文采用的最大积分波数为kmax=300-500，最小积分的间隔为Δk=0.00005-0.0001.
对于式 (15) 中的求和公式，要保证其精度必须采用足够大的项数Nmax.经研究发现，Nmax的取值与入射频率η、凸起间距离L和介质的阻尼比ζ有关，较小的频率η、较小L和较小的ζ，需要较大的Nmax.以均匀半空间中均质半圆凸起地形为例，表 1和表 2分别给出了间距L/a=4.0和L/a=8.0两种情况，SH波入射角度为ψ=45°、入射频率分别为η=ωa/πcs=0.5和2.0时 (a为凸起半径，cs为材料剪切波速)，“0”号凸起地形附近位置点的位移幅值随Nmax的变化，不考虑阻尼 (实际取ζ=0.001).从表 1和表 2可以看出，无论是较低频率 (η=0.5) 还是较高频率 (η=2.0)，无论是较小的凸起间距 (L/a=4.0) 还是较大的凸起间距 (L/a=8.0)，位移幅值随着求和项数的增大均有着良好的收敛性.考虑到后面频率的计算范围为η=0.0~6.0，针对L/a=4.0和8.0两种情况，后面计算中取Nmax=20~240来满足求和公式的精度.
|
表 1 间距L/a=4.0周期凸起地形各位置点的位移幅值随项数Nmax收敛情况 Table 1 Convergence of displacement amplitudes versus Nmax at some positions for L/a=4.0 |
|
|
表 2 间距L/a=8.0周期凸起地形各位置点位移幅值随项数Nmax收敛情况 Table 2 Convergence of displacement amplitudes versus Nmax at some positions for L/a=8.0 |
Liang和Fu (2011)采用正余弦级数展开方法求解了半椭圆凸起对SH波的散射问题，文中首次给出了扁长凸起结果.本文通过与Liang和Fu (2011)解析结果的比较来验证方法的正确性.计算中，扁长椭圆凸起高度与半宽之比为b/a=10.0，无量纲频率η=ωa/πcs=1.0，其中a为半圆凸起半径，cs为剪切波速，SH波入射角度分别取ψ=0°、30°、60°和90°(其中，ψ=0°表示水平入射).且经验算，本文方法在入射频率为η=1.0时，将凸起间距离取为L/a=160.0，可退化到单一凸起结果.从图 4可以看出，本文的结果 (右侧) 与Liang和Fu (2011)给出结果 (左侧) 非常吻合，从而证明了本文方法的正确性.
|
图 4 本文地表位移幅值与梁建文和付佳 (2011)扁长凸起 (b/a=10.0) 结果的比较 (a) 梁建文和付佳的结果 (2011); (b) 本文方法的结果. Fig. 4 Comparison of surface computed results with those of Liang and Fu (2011) for the prolate hill (b/a=10.0) (a) Results of Liang and Fu (2011); (b) Results of the present method. |
为研究凸起间距对凸起附近地表位移幅值的影响，以均匀半空间中均质周期凸起地形 (凸起材料与半空间材料相同) 为例，图 5给出了凸起间距分别为L/a=4.0、8.0和∞时 (其中L/a=∞代表单一凸起地形)，“0”号凸起附近位移的幅值.定义无量纲频率为η=ωa/πcs，其中a为凸起半径，cs为剪切波速，不考虑阻尼影响 (实际取阻尼比ζ=0.001，以消除波数积分中的奇异性问题).图 5计算中分别取无量纲频率η=0.5、1.0和2.0，SH波入射角度分别取ψ=5°、45°和90°(ψ=90°表示垂直入射).无量纲位移幅值为|v/v0|，其中v0可认为是输入SH波的幅值.
|
图 5 “0”号凸起附近位移幅值 (均匀半空间，不同凸起间距) Fig. 5 Surface displacement amplitudes around 0th hill (different distance between hills in a uniform half-space) |
从图 5中可以看出，周期凸起地形 (L/a=4.0和8.0) 与单一凸起地形 (L/a=∞) 对SH波的散射有着显著的差别.单一凸起情况：SH波掠入射 (ψ=5°) 时，由于凸起地形波入射一侧角点的较强散射作用，凸起左侧位移幅值振荡明显，且大于凸起右侧位移幅值；位移幅值的最大值多出现在波垂直入射 (ψ=90°) 情况；整体上，随着入射频率的逐渐增大，位移幅值逐渐增大，且在高频 (η=2.0) 时，凸起表面的位移幅值明显大于凸起两侧水平地表的位移幅值.周期凸起情况：SH波掠入射 (ψ=5°) 时，由于前方凸起地形对波的耗散效应，使得位移幅值明显小于单一凸起地形情况，同时凸起地形对SH波的“放大效应”也明显减弱，凸起左右两侧的位移幅值差异减小，这种现象在高频率时尤为明显；SH波斜入射 (ψ=45°) 和垂直入射 (ψ=90°) 时，由于凸起地形间的动力相互作用，则使得位移幅值大于单一凸起地形情况；位移幅值的最大值仍出现在波垂直入射情况，但整体上凸起表面的位移幅值与其两侧的位移幅值变的较为接近；随着入射频率的逐渐增大，位移幅值反而逐渐减小.
比较凸起间距离L/a=4.0和L/a=8.0情况发现，周期凸起地形间距离对凸起附近位移幅值亦有着显著的影响.SH波掠入射 (θ=5°) 时，随着凸起间距离的减小，凸起附近位移幅值进一步减小，凸起左右两侧的位移幅值变得更为接近.SH波斜入射 (θ=45°) 和垂直入射 (θ=90°) 时，随着凸起间距离的减小，凸起间动力相互作用进一步增强，凸起两侧位移幅值逐渐增大.入射频率对位移幅值也有着重要影响，整体上，频率为η=0.5和1.0时，间距为L/a=8.0与L/a=4.0对应的位移幅值差别较大，而η=2.0时两者对应的位移幅值差别较小.另外，对于间距为L/a=8.0的周期凸起地形，对所有入射角度，x/a=-4.0和x/a=4.0处位移幅值相等，而对于间距L/a=4.0的周期凸起地形，位移幅值在区间 (-4.0≤x/a≤0.0) 和区间 (0.0≤x/a≤4.0) 两段上重复出现，这些现象也进一步证明了方法的正确性.
5.2 凸起地形材料对位移幅值的影响为研究凸起地形的材料对地表位移幅值的影响，以均匀半空间中异质 (凸起材料与半空间材料不同) 周期凸起地形为例，图 6给出了凸起与均匀半空间剪切波速比值分别为csH/csR=0.5、1.0和2.0(其中csH/csR=1.0为均质凸起) 时，“0”号凸起附近的位移幅值.凸起与均匀半空间密度比ρH/ρR=1.0，不考虑材料阻尼 (实际取阻尼比ζ=0.001).定义无量纲频率η=ωa/πcsR.图 6计算中分别取无量纲频率η=0.5、1.0和2.0，SH波入射角度分别取ψ=5°、45°和90°(ψ=90°表示垂直入射)，凸起间距L/a=8.0.无量纲位移幅值为|v/v0|，其中v0仍可认为是输入SH波的幅值.
|
图 6 “0”号凸起附近位移幅值 (均匀半空间，不同硬度材料凸起) Fig. 6 Surface displacement amplitudes around 0th hill (different material of the hill in a uniform half-space) |
从图中可以看出，凸起材料对凸起附近的地表位移幅值有着显著的影响.整体上随着凸起介质逐渐变硬 (剪切波速比越大)，凸起附近地表位移幅值逐渐减小，尤其是凸起表面的位移幅值，且这种现象在SH波以较高频率 (η=2.0) 入射时尤为明显.同样，随着凸起介质逐渐变硬，各凸起间的动力相互作用也逐渐减弱，表现为凸起地表两侧的地表位移幅值逐渐减小，但这种情况在SH波以η=0.5和ψ=90°入射时，却正好相反.另外，从图中还可以看出，当凸起地形材料强度较低 (csH/csR=0.5) 时，凸起地表面位移幅值的最大值，较多情况出现在了凸起两侧位置 (山腰)，而不是出现在凸起地形的顶部.
5.3 凸起截面形状对位移幅值的影响为研究凸起截面形状对周期凸起附近地表位移幅值的影响，以均匀半空间中均质周期半圆、等腰梯形和等腰三角形凸起为例，图 7给出了三种截面情况，“0”号凸起附近的位移幅值.图 7计算中分别取无量纲频率η=0.5、1.0和2.0，SH波入射角度分别取ψ=5°、45°和90°(ψ=90°表示垂直入射)，凸起间距L/a=8.0.半圆、等腰三角和等腰梯形凸起具有相同高度以及相同的底面宽度，且等腰梯形的上底宽度为下底宽度的0.5倍.图 7中无量纲频率和无量纲位移幅值的定义方式均同图 5.
|
图 7 “0”号凸起附近位移幅值 (均匀半空间，不同凸起截面形状) Fig. 7 Surface displacement amplitudes around 0th hill (different cross section of hill in a uniform half-space) |
从图中可以看出，凸起的截面形状对周期凸起附近的位移幅值亦有一定程度的影响，且影响又依赖于波的入射频率和入射角度.整体上，三角凸起表面的位移幅值最大，半圆凸起表面位移次之，梯形凸起表面位移幅值最小，表现为顶点越尖锐，顶部位移幅值越大，且这种情况在波垂直入射 (ψ=90°) 时尤为明显.另外，频率较低 (η=0.5) 时，三种截面形状凸起对应的位移幅值差异很小，随着频率的增大，三种截面情况对应的位移幅值差异逐渐增大，说明凸起的截面形状对高频入射波更为敏感.
5.4 层状场地中周期凸起位移幅值为研究土层对周期凸起附近地表位移幅值的影响，以基岩上单一土层中均质周期凸起地形 (凸起地形材料与土层材料相同) 为例，图 8给出了“0”号凸起附近位移幅值.基岩介质由剪切波速csR、质量密度ρR和阻尼比ζR确定，土层介质由剪切波速csL、质量密度ρL和阻尼比ζL确定.图 8中所采用计算参数如下：基岩与土层剪切波速比为csR/csL=2.0，基岩与土层密度比为ρR/ρL=1.0，基岩阻尼比为ζR=0.02，土层阻尼比为ζL=0.05，土层厚度分别为H/a=0.5、1.0和∞(其中H/a=∞代表均匀半空间情况)，凸起间距L/a=8.0，无量纲频率分别为η=0.5、1.0和2.0，SH波入射角度分别为ψ=5°、45°和90°(ψ=90°表示垂直入射).无量纲位移幅值为|v/v0|，其中v0为输入SH波的幅值.
|
图 8 “0”号凸起附近位移幅值 (基岩上单一土层，不同土层厚度) Fig. 8 Surface displacement amplitudes around 0th hill (different thickness of layer in a single layered half-space) |
比较图 8均匀半空间 (H/a=∞) 与图 5均匀半空间结果，发现介质的材料阻尼对周期凸起地形附近的动力响应有着显著的影响.有阻尼半空间中 (图 8) 各凸起地形间的动力相互作用明显小于无阻尼半空间中 (图 5) 凸起地形间的动力相互作用，且这种现象在高频 (τ=2.0) 时尤为明显.另外，介质的阻尼也使得凸起附近位移幅值显著降低.
比较层状半空间 (H/a=0.5和1.0) 与均匀半空间 (H/a=∞) 结果发现，层状半空间与均匀半空间结果存在显著差异.层状半空间情况，SH波垂直入射 (ψ=90°) 和斜入射 (ψ=45°) 时凸起附近位移幅值较大，掠入射 (ψ=5°) 时位移幅值较小，然而均匀半空间中位移幅值则在掠入射时较大，垂直入射时较小.这是因为对应给定入射频率和入射角度，层状半空间中凸起附近动力响应由层状场地的自身动力特性 (层状场地自由波场) 和凸起截面形状 (凸起形成的散射波场) 共同决定，而均匀半空间情况则仅取决于凸起截面形状 (凸起形成的散射波场).
土层厚度的改变直接导致了自由波场的改变，进而也导致了散射波场的改变，使得不同厚度基岩上单一土层中凸起地形附近位移幅值显著不同.SH波垂直入射 (ψ=90°) 且入射频率为η=1.0时，H/a=0.5对应的位移幅值显著大于H/a=1.0和H/a=∞情况对应的位移幅值，这是因为η=1.0接近于土层厚度H/a=0.5场地的第一共振频率.同样，由于η=0.5接近于土层厚度H/a=1.0场地的第一共振频率，SH波垂直入射 (ψ=90°) 时，H/a=0.5情况对应的地表位移幅值显著大于H/a=1.0和H/a=∞对应位移幅值.
5.5 周期凸起时域结果
为研究周期凸起地形在平面SH波入射下的时域响应，仍以均匀半空间均质半圆凸起地形 (凸起材料与半空间材料相同) 为例，进行了Ricker时程输入下凸起附近表面位移幅值的求解.输入Ricker时程为
首先，比较图中周期凸起 (L/a=4.0和8.0) 和单一凸起 (L/a=∞) 的结果，发现两者间的时域响应有着明显的差异.周期凸起情况较单一凸起情况，凸起附近波的传播更为复杂且波的持续时间更长.以SH波垂直入射 (ψ=90°) 为例进行说明.单一凸起情况，波的传播相对简单，“0”号凸起右侧区间 (1.0≤x/a≤4.0) 的观测点，直达波D1到达后，SR01、SR02和SR03等散射波相继到达，“0”号凸起左侧区间 (-4.0≤x/a≤-1.0) 的观测点亦有着相同的规律，散射波分别为SL01、SL02和SL03.图 9中，字母D1表示直达波，字母“SL0”和“SR0”分别表示从“0”号凸起左右角点发出的散射波.“0”号凸起表面区间 (-1.0≤x/a≤1.0) 的观测点，散射波在凸起两角点间往复传播.
|
图 9 “0”号凸起附近周期凸起时域结果 (L/a=4.0、8.0和∞) Fig. 9 Synthetic seismograms around the 0th hill (L/a=4.0, 8.0 and ∞) |
周期凸起地形情况，波的传播变得复杂，相邻凸起角点发出的散射波也被观察到.以L/a=8.0且SH波垂直入射 (ψ=90°) 为例进行说明.位于“0”号凸起右侧区间 (1.0≤x/a≤4.0) 的观测点，直达波D1到达后，SR01、SL11、SR02、SL12等散射波相继达到，位于“0”号凸起左侧区间 (-4.0≤x/a≤-1.0) 的观测点亦有着相同的规律，散射波分别为SL01、SR-11、SL02和SR-12.字母“SL1”和“SR-1”分别表示从“1”号凸起左侧角点和“-1”号凸起右侧角点发出的散射波.位于“0”号凸起表面区间 (1.0≤x/a≤1.0) 的观测点，往复于“0”号凸起左右角点间的散射波得到加强，这是由于“1”号凸起左侧角点和“-1”号凸起右侧角点传来散射波的参与造成的.
相对于L/a=8.0情况，随着凸起间距离的减小，L/a=4.0时波动变得更为复杂且持续时间进一步延长.SH波斜入射 (ψ=45°) 和垂直入射 (ψ=90°) 时，L/a=4.0情况较L/a=8.0情况，位移幅值更大且振荡更为剧烈，表现为凸起间更强的动力相互作用.SH波掠入射 (ψ=5°) 时，L/a=4.0情况较L/a=8.0情况，位移幅值明显要小，表现为“0”号凸起前方凸起对波更强的耗散作用.以上分析也说明了图 9中的时域结果与图 5中的频率结果规律是一致的.
在周期凸起L/a=8.0情况下，不考虑时间延迟时，在x/a=-4.0的波形和在x/a=4.0的波形完全一致.在周期凸起L/a=4.0情况下，不考虑时间延迟的情况下，区间-2.0≤x/a≤2.0内的波形与区间2.0≤x/a≤6.0内的波形完全相同.注意到x/a=-4.0和x/a=4.0两点间距离恰为L/a=8.0，而区间-2.0≤x/a≤2.0与区间2.0≤x/a≤6.0内对应点间的距离恰为L/a=4.0.正是依据以上特征，提出了本文的求解方法.同时，由式 (8) 和图 9可以看出，随着入射角度的逐渐增大，地表视速度逐渐增大，观测点间的时间差逐渐减小，当波垂直入射 (ψ=90°) 时，时间差为零.
6 结论本文在层状半空间平面外精确动力刚度矩阵和斜线荷载动力格林函数的基础上，给出了层状半空间中周期分布斜线荷载动力格林函数，并进一步建立了求解层状半空中周期分布地形平面外散射问题的间接边界元方法.验证了方法的正确性并以均匀半空间和基岩上单一土层中周期分布凸起地形为例进行了数值计算分析，得到了以下主要结论.
(1) 本文给出了一种用于求解层状半空间中周期分布地形平面外散射问题的新方法，方法仅需针对一个凸起地形进行边界单元离散和求解，显著降低了求解自由度和计算时间，避免了无限边界截断引入的误差，同时该方法具有与解析法一样物理概念清晰的优点，且经验证该方法具有较高的精度.
(2) 周期凸起地形与单一凸起地形对平面SH的散射有着显著的差别.周期凸起地形在SH波掠入射时，由于前方凸起对波的耗散作用，使得周期凸起地形位移幅值显著小于单一凸起情况，SH波斜入射和垂直入射时，由于凸起地形间的动力相互作用，则使得周期凸起地形位移幅值大于单一凸起情况；凸起地形间距离对凸起附近位移幅值亦有着显著的影响.随着凸起间距离的减小，凸起间的动力相互作用增强，同时SH波掠入射时，前方凸起对波的耗散作用也越强.
| Bard P Y. 1982. Diffracted waves and displacement field over two-dimensional elevated topographies. Geophys. J. Int., 71(3): 731-760. DOI:10.1111/j.1365-246X.1982.tb02795.x | |
| Boore D M. 1972. A note on the effect of simple topography on seismic SH waves. Bull. Seismol. Soc. Am., 62(1): 275-284. | |
| Bouchon M. 1973. Effect of topography on surface motion. Bull. Seismol. Soc. Am., 63(2): 615-632. | |
| Bouchon M, Schultz C A, Toksöz M N. 1996. Effect of three-dimensional topography on seismic motion. J. Geophys. Res., 101(B3): 5835-5846. DOI:10.1029/95JB02629 | |
| Çelebi M. 1987. Topographical and geological amplifications determined from strong-motion and aftershock records of the 3 March 1985 Chile earthquake. Bull. Seismol. Soc. Am., 77(4): 1147-1167. | |
| Du Y J, Liu D K, Guo F, et al. 2009. Surface displacement function of two isosceles-triangular hills during incidence of SH-wave. Journal of Earthquake Engineering and Engineering Vibration (in Chinese), 29(3): 1-8. | |
| Fu L Y. 2005. Rough surface scattering:comparison of various approximation theories for 2D SH waves. Bull. Seismol. Soc. Am., 95(2): 646-663. DOI:10.1785/0120040081 | |
| Han F, Wang G Z, Chen H. 2013. Research on scattering of SH-waves on multiple hills and canyons. Applied Mathematics and Mechanics (in Chinese), 34(4): 355-363. | |
| Lee V W, Luo H, Liang J W. 2006. Antiplane (SH) waves diffraction by a semicircular cylindrical hill revisited:An improved analytic wave series solution. J. Eng. Mech., 132(10): 1106-1114. DOI:10.1061/(ASCE)0733-9399(2006)132:10(1106) | |
| Liang J W, Zhang Y S, Lee V W. 2006. Surface motion of a semi-elliptical hill for incident plane SV waves:analytical solution. Acta Seismologica Sinica (in Chinese), 28(3): 238-249. | |
| Liang J W, Ba Z N. 2008. Surface motion of a hill in layered half-space subjected to incident plane SH waves. Journal of Earthquake Engineering and Engineering Vibration (in Chinese), 28(1): 1-10. | |
| Liang J W, Fu J. 2011. Surface motion of a semi-elliptical hill for incident plane SH waves. Earthquake Science, 24(5): 447-462. DOI:10.1007/s11589-011-0807-1 | |
| Liu G, Chen H T, Liu D K, et al. 2010. Surface motion of a half-space with triangular and semicircular hills under incident SH waves. Bull. Seismol. Soc. Am., 100(3): 1306-1319. DOI:10.1785/0120090273 | |
| Liu J B. 1996. The effect of local irregular topography on seismic ground motion. Acta Seismologica Sinica (in Chinese), 18(2): 239-245. | |
| Pedersen H A, Sánchez-Sesma F J, Campillo M. 1994. Three-dimensional scattering by two-dimensional topographies. Bull. Seismol. Soc. Am., 84(4): 1169-1183. | |
| Sánchez-Sesma F J. 1983. Diffraction of elastic waves by three-dimensional surface irregularities. Bull.Seismol.Soc.Am., 73(6A): 1621-1636. | |
| Sánchez-Sesma F J, Herrrra I, Avilés J. 1982. A boundary method for elastic wave diffraction:Application to scattering of SH waves by surface irregularities. Bull. Seismol. Soc. Am., 72(2): 473-490. | |
| Sills L B. 1978. Scattering of horizontally-polarized shear waves by surface irregularities. Geophys. J. Int., 54(2): 319-348. DOI:10.1111/j.1365-246X.1978.tb04263.x | |
| Smith W D. 1975. The application of finite element analysis to body wave propagation problems. Geophys. J. Int., 42(2): 747-768. | |
| Sohrabi-Bidar A, Kamalian M. 2013. Effects of three-dimensionality on seismic response of Gaussian-shaped hills for simple incident pulses. Soil Dyn. Earthq. Eng., 52: 1-12. DOI:10.1016/j.soildyn.2013.04.009 | |
| Spudich P, Hellweg M, Lee W H K. 1996. Directional topographic site response at Tarzana observed in aftershocks of the 1994 Northridge, California, Earthquake:Implications for mainshock motions. Bull. Seismol. Soc. Am., 86(1B): S193-S208. | |
| Takenaka H, Kennett B L N, Fujiwara H. 1996. Effect of 2-D topography on the 3-D seismic wavefield using a 2. 5-D discrete wavenumber-boundary integral equation method. Geophys. J. Int., 124(3): 741-755. | |
| Tsaur D H. 2011. Scattering and focusing of SH waves by a lower semielliptic convex topography. Bull. Seismol. Soc. Am., 101(5): 2212-2219. DOI:10.1785/0120100324 | |
| Wolf J P. Dynamic Soil-Structure Interaction.Englewood Cliffs: Prentice-Hall, 1985. | |
| Yang Z L, Xu H N, Hei B P, et al. 2014. Antiplane response of two scalene triangular hills and a semi-cylindrical canyon by incident SH-waves. Earthquake Engineering and Engineering Vibration, 13(4): 569-581. DOI:10.1007/s11803-014-0264-7 | |
| Yuan X M, Men F L. 1992. Scattering of plane SH waves by a semi cylindrical hill. Earthq. Eng. Struct. Dyn., 21(12): 1091-1098. DOI:10.1002/(ISSN)1096-9845 | |
| Yuan X M, Liao Z P. 1996. Scattering of plane SH waves by a cylindrical hill of circular-arc cross-section. Earthquake Engineering and Engineering Vibration (in Chinese), 16(2): 1-13. | |
| Zhang B, Papageorgiou A S, Tassoulas J L. 1998. A hybrid numerical technique, combining the finite-element and boundary-element methods, for modeling the 3D response of 2D scatterers. Bull. Seismol. Soc. Am., 88(4): 1036-1350. | |
| 杜永军, 刘殿魁, 郭凤, 等. 2009. 双等腰三角形凸起地形在SH波入射时的地表位移函数. 地震工程与工程振动, 29(3): 1–8. | |
| 韩峰, 王光政, 陈翰. 2013. SH波对多个凸起与凹陷相连地形的散射问题研究. 应用数学和力学, 34(4): 355–363. | |
| 梁建文, 张彦帅, LeeV W. 2006. 平面SV波入射下半圆凸起地形地表运动解析解. 地震学报, 28(3): 238–249. | |
| 梁建文, 巴振宁. 2008. 弹性层状半空间中凸起地形对入射平面SH波的放大作用. 地震工程与工程振动, 28(1): 1–10. | |
| 刘晶波. 1996. 局部不规则地形对地震地面运动的影响. 地震学报, 18(2): 239–245. | |
| 袁晓铭, 廖振鹏. 1996. 任意圆弧形凸起地形对平面SH波的散射. 地震工程与工程振动, 16(2): 1–13. | |