0 引言

地形起伏对于强地面运动和区域震相的解释有重要影响.在过去的几十年中，地形散射效应的研究一直是地震学研究的热点问题之一(例如Boore, 1972, 1973; Bouchon，1973; Geli et al., 1988; Sánchez-Sesma and Campillo, 1993; Bouchon et al., 1996).在声学领域中，起伏表面会引起散射衰减以及振幅和相位扰动，其散射强度取决于表面粗糙度与入射波长的相对关系(Hu et al., 2009).自从Rayleigh的经典著作发表以来，粗糙表面散射级数逼近方法有了很多的进展(例如，Gilbert and Knopoff, 1960；Kennett，1972；Hudson et al., 1973；Snieder，1986).这些散射近似涉及的级数收敛问题取决于散射算子逼近的函数类型和展开项数(Yu and Fu, 2012).Fu(2005)对各种粗糙表面散射逼近方法及其收敛特性进行了综述.

边界积分方程方法由于可以灵活处理粗糙表面而被广泛用作粗糙表面散射的精确模拟工具(例如，Sánchez-Sesma et al., 1982；Dravinski，1983；Kawase，1988；Sánchez-Sesma and Campillo, 1991；Fu et al., 2002；Fu，2003).Cao等研究了入射SH波对各种地形的散射(2004)，评估了Aki-Larner方法(Aki and Larner, 1970；Bard，1982)，Bouchon-Campillo方法(Bouchon，1985；Campillo and Bouchon, 1985; Fu and Bouchon, 2004)和全球广义反射/透射矩阵方法(Chen, 1990, 1995, 1996)的准确性、计算效率和稳定性.这些全波形数值方法由于需要对每个频率进行矩阵求逆，在高频情况下，需要的计算量巨大.一种简单的方法是直接利用边界积分方程离散系数矩阵来分析崎岖地表地震散射的复杂性(Chen et al., 2017; Liu et al., 2017)，但是这些方法不能计算波形数据.Fu和Wu(2001)及Ge等(2005)通过基于边界积分的波场连接技术，将一个大的计算区域分解成多个小区域，从而提高了全波形数值模拟的计算效率.

Zhou和Chen(2006)通过对Bouchon-Campillo方法(BC)进行分析，根据格林函数的特性提出了一种新的局域离散波数方法，该方法只需直接求解起伏地表部分的位移，极大提高BC方法的计算效率.由于BC方法在水平方向均匀采样，对于高倾角地形，其计算结果精度不足(Zhou and Chen, 2006).边界元方法沿起伏边界离散，可以很好处理高倾角问题.因此，本文将采用类似Zhou和Chen(2006)的方式进行处理，提出了一种有效模拟起伏地形散射的局域边界元法.本方法对起伏地形散射的边界积分进行局域化处理，即将起伏地形分解成水平部分和起伏部分，水平部分的位移可以表示成起伏部分位移解的线性组合，因此只需求解起伏部分位移的线性方程.我们首先对二维SH波的地表散射问题的全波形边界元方法进行简单回顾，然后根据格林函数的特点，提出了在不损失任何精度的情况下，计算地表散射的局域边界元方法，并通过与解析解，以及全波边界元方法的比较，验证了方法的正确性.

1 方法介绍 1.1 边界元方法简介

二维稳态SH波入射到起伏表面时，其位移u(x)在满足自由边界条件t=0时，其满足的边界积分方程为

(1)

其中系数C(x)在地表以下为1，在地表Γ上等于1/2，而表示边界Γ上的点ξ上外法线方向的方向导数.G(x, ξ)是均匀介质中位移的格林函数，有

(2)

其中, 为波数, H₀⁽¹⁾(kr)零阶第一类汉克尔函数，r=|x－ξ|.

在边界元方法中，我们首先将自由边界Γ离散成NE个单元，第i个单元用表示，每个单元的两端为节点，共NP个节点(NP=NE+1).经过边界离散化后，方程(1)变为

(3)

线性单元边界元中假设每个单元上的位移为线性分布，位于节点e₁和e₂之间的单元Γ_e上的位移可以利用定义域为η∈[－1, 1]上的形函数ϕ₁(η)=以及节点e₁和e₂上的位移u(e₁)和u(e₂)表示，

(4)

则有

(5)

边界积分方程(1)变为

(6)

.令

(7)

方程(7)可以写成矩阵形式

(8)

其中C为对角矩阵，H元素由(7)式给出, u为边界Γ上SH波的位移.传统的边界元方法，直接求解方程(8)来得到地表的位移，方法涉及到矩阵求逆，当计算模型增大时，计算时间指数增加.下面我们将根据格林函数的特点，提出一种计算地表散射问题的局域边界元方法.

1.2 局域边界元方法

如图 1所示，我们将边界部分分为水平边界和起伏边界两部分，位于水平边界上的位移定义为u_F，位于起伏边界上的位移定义为u_C，(8)式则可以写为

(9)

对于SH波，水平边界上，由于垂直于r, 所以有

(10)

所以可知在水平界面上，H_FF=0，同时，水平边界上的位移u_F满足

(11)

整理可得

(12)

即水平部分的位移场可以由其入射波场和起伏表面的散射波场完全确定，因此我们只需求解起伏表面上的波场，就可以求得整个表面区域的波场.

对于起伏部分的节点，其位移场满足

(13)

将(12)式代入(13)式，可得

(14)

整理可得

(15)

再将(15)式中的结果代入(12)式中，就可以得到整个地表的地震波场.

2 方法验证

首先，我们通过计算不同频率平面SH波垂直入射到半圆柱形峡谷时的地表位移来检验新方法的有效性.Trifunac(1973)给出了这一经典模型的解析解，从而为各种数值方法提供了参考.将本研究中该模型在垂直入射SH波作用下的计算结果与解析解(Trifunac, 1973)比较，如图 2所示.图中显示了无量纲频率η=0.5、1.5、2.5时本文的计算结果与解析解的比较，dx/a=0.05，L/a＝20，η=2a /λ，λ是入射平面SH波的波长.细实线为解析解的计算结果，粗虚线为本文的计算结果，比较结果表明，本方法所得结果与解析解吻合得很好，从而验证了方法的正确性.

3 初步应用和效率分析 3.1 理论模型

作为方法的应用测试，我们设计了一个楔形的地形，如图 3a所示，这一地形的右侧垂直断崖，对于均匀离散方法(如BC方法)存在巨大的挑战，我们将模型沿地表起伏地形每20 m一个单元，共419个单元，震源时间函数为中心频率f=10 Hz的Ricker子波，我们分别用传统边界元方法和局域边界元方法进行了计算，分别选取位于不同位置的8个台站的数据进行对比.这8个台站中，其中1号和176号位于水平的界面处，26、126、151三个台站位于拐角处，51，76，101三个台站位于倾斜的斜坡上.对比结果展示在图 3b中，可以看到由于在局域边界元方法没有做任何假设，因此其计算结果与边界元方法完全相同.

图 4展示了整个表面上所有台站的位移，从图中可以看到，对与垂直入射的平面波，在几个转折点都产生了散射波场.B点产生的散射沿边界传播到两侧；由于C、D两点距离较近，在两点之间产生了多次往复的散射，因此在山顶和地面之间可以产生震荡.这一结果在频率-空间域的波场图中也可以清楚地看到，如图 5所示，在B-C、C-D之间，随着频率的增加，都可以看到明显的地震波的往复震荡. 图 6分别给出了f=2.5, 5.0, 7.5, 10.0, 12.5 Hz时地表振幅随观测点位置的变化，图中实线表示地表位移的实部，虚线表示地表位移的虚部，从波形图中可以更加清楚地看到转折点之间的震荡变化.

本例计算时所采用的模型，水平部分的单元和起伏部分的单元数量相同，均为320个节点.传统的边界元方法每个频率中需要求解n×n个线性方程组，若采用高斯消去法求解矩阵，需要次乘法运算和次加法运算，而采用局域边界元方法，只需进行一次大小为的矩阵乘法，一次矩阵加法和一次的高斯消去法，其中矩阵乘法的计算复杂度约为O(n³/8), 矩阵加法的计算复杂度为O(n²/4)，高斯消元法为O(n³/12)，因此总计算量约为O(5n³/24)，约为传统边界元方法的5/16，显著提高了计算效率.

3.2 实际应用

2019年12月26日18时36分，湖北省发生了M_L4.9地震，震中位置为北纬30.87°，东经113.4°，震源深度10 km.为了考察不规则起伏地形地表起伏对地震波散射的影响，我们利用USGS提供的全球高分辨率高程数据GMTED2010 (https://www.usgs.gov/centers/eros/science/usgs-eros-archive-digital-elevation-global-multi-resolution-terrain-elevation?qt-science_center_objects=0#qt-science_center_objects)沿着北纬30.83°在东经90°—120°之间提取出一条地形数据，数据分辨率为30 s(约等于0.8 km).地形数据如图 7a所中蓝色线所示，为了利用本文给出的局域边界元方法，我们对模型的两端进行了简化，变成了水平界面(如图 7a中红线所示)，分别用传统的边界元方法和局域边界元法对两个模型计算SH波的位移.图 7b展示了两种方法得到的18个台站上得到的波形，其中台站位置如图 7a中的三角形所示，分别位于模型设置的地表上.震源位置取在x=1400 km的位置，震源深度10 km，震源时间函数为中心频率1 Hz的Ricker子波.为了展示散射波形，我们对波形进行了放大处理，并按照直达波的峰值进行了对齐.可以清楚地看到，由于地表不规则地形的存在，在直达波之后存在长时间的散射尾波.两组结果对比表明，在模型的中部两种方法得到的结果几乎完全一致，说明在两端模型的改变对散射波场影响很小.即使是模型的最右侧的三个台站，台站高程存在一定差别，计算得到的散射波形也非常接近，说明在地形相对平缓的区域，我们可以用水平地表代替起伏地表.而在最左侧的台站，在直达波到达后0~50 s之间，两种模型计算得到的地震波形非常接近，而50 s之后原始模型计算结果和简化模型中的波形存在较大差异.这是因为50 s之前的尾波是震源到台站路径上的前向散射产生的; 而50 s之后，原始模型中尾波是由左侧复杂的起伏地形产生的后向散射波，简化模型由于忽略了左端的起伏地形，后向散射波较弱.

4 结论与讨论

本文提出了一种用于研究不规则地形地震波散射问题的有效方法——局域边界元法.在新方法中，我们根据边界元方法中水平地表单元间格林函数为零的特点，将自由边界分成水平部分和起伏部分，解析地推导出水平部分位移与起伏部分位移之间的线性关系，只需直接求解不规则部分中的位移极大地减少了未知数，提高了边界元方法的计算效率.

通过与解析方法对半圆柱形盆地模型SH波场散射结果的比较，验证了新方法的正确性.局域边界元方法对楔形模型地震波场的模拟结果与边界元法完全相同.在水平地形与起伏地形单元数大致相等的情况下，计算效率可以提高近70%.很多情况下，由于需要考虑介质边缘的边界效应，我们需要扩大模型，当模型增大时，新方法的计算量增加很少.因此，局域边界元方法有望作为模拟不规则地形地震波散射问题的有效工具.

值得指出的是，由于我们的方法利用了水平界面上矢径方向和法线方向相互正交而形成的单元间格林函数为零这一特性对边界元方法进行了简化，提高了计算效率.因此在忽略地球曲率的情况下，用来计算地表散射问题.当问题的尺度大到必须考虑地球的曲率时，Sun和Kennet(2016)将传统的层析成像方法应用到了球坐标系下，避免了展平近似带来的误差.但是在球坐标系下，即使时平坦的球面，格林函数为零的特性也不能满足，因此目前还不能将此方法推广到球坐标系下.将边界元方法应用到球坐标系下，也是未来可行的改进方向.

致谢  在本文的修改过程中，审稿专家对结果论证和应用实例方面都提出了许多建设性的意见，使作者受益良多，在此对审稿专家的辛勤工作表示衷心的感谢.