极坐标系有限差分中起伏地表边界条件处理

徐剑侠， 张伟， 陈晓非*

1 中国科学技术大学地球和空间科学学院， 合肥 230026 2 南方科技大学地球与空间科学系， 广东深圳 518055

0 引言

有限差分算法是地震学中广泛应用的重要算法(Alterman and Karal, 1968; Bayliss et al., 1986; Saenger and Bohlen, 2004)，可以简洁高效地模拟二维/三维非均匀模型中地震波波场，尤其是在Virieux(1986)提出了基于应力-速度的一阶差分格式后得到了广泛的应用.

地形起伏会明显影响地震波场的传播.为了在有限差分中精确处理起伏地表，需要解决网格划分和自由地表边界条件引入两个问题.对于网格划分，现有工作已经研究比较完善，主要有规则网格，非结构网格和垂向变换/贴体网格等方法.规则网格容易生成和使用，最简单的方法是在加密网格后采用台阶状网格近似描绘起伏地表(Robertsson, 1996; Wang and Zhang, 2004)；非结构网格适用性好，较易生成网格，但是对精度有一定的影响(Käser et al., 2001)；垂向变换/贴体网格思路相近，依照实际地表起伏划分曲线网格，使用比较方便，适用于同位网格和交错网格(Hestholm and Ruud, 1994, 1998; Zhang and Chen, 2006)，但是对网格变形范围具有一定的限制(Tessmer et al., 1992; Hestholm and Ruud, 1994).自由地表边界条件近年来具有较大发展.对于速度分量法向导数的需求，多通过牵引力为零边界条件，转化为速度分量水平导数计算.对于应力分量法向导数计算，主要实现方法大体有：单边差分方法实施方便，通过在地表格点上采用单边差分计算，可以满足全局二阶空间差分精度，并被成功拓展到复杂介质(Nilsson et al., 2007; Appelö and Petersson, 2009; Lan and Zhang, 2011,2012; 兰海强等, 2011; Liu et al., 2018)；应力镜像法是地表边界条件在水平规则地表下的简化表述，因此直接用于起伏地表存在一定近似(例如Robertsson, 1996)，部分工作也考虑了起伏地表法向对牵引力的影响，和牵引力镜像法物理思路一致(Solano et al., 2016)；牵引力镜像法考虑了起伏地表法向量变化，通过牵引力关于地表反对称表述自由地表边界条件，对于自由地表边界条件表述更加准确，在起伏地表边界满足四阶空间差分精度(Zhang and Chen, 2006; Zhang et al., 2012b)；Taylor多项式方法适用于规则网格，通过引入Taylor多项式拟合外推虚拟点内波场信息，同时考虑了自由地表法向影响(Lombard et al., 2008)；在地球物理勘探领域不关注地表反射，也有工作将吸收层作用于起伏地表(Rao and Wang, 2013, 2018)；还有工作通过样条或者插值多项式外推物理量，进而处理自由地表边界条件(Mulder, 2017; Mulder and Huiskes, 2017)等.上述方法已经可以在笛卡尔坐标系框架下较好地处理起伏地表边界条件，其中基于贴体网格牵引力镜像法(Zhang and Chen, 2006)在贴体网格中将牵引力关于地表反对称，可以简洁清晰地表征自由地表物理边界条件，同时具有较高空间差分精度.

当研究区域-全球尺度的地震波传播问题时，必须要考虑地球曲率的影响.虽然可以通过网格变形(Appelö and Petersson, 2009)或者展平变换(Li et al., 2014)等方法处理，但是选用极坐标系(二维)或球坐标系(三维)处理起来更为简洁方便(Igel and Weber, 1995; Toyokuni and Takenaka, 2006; Zhang et al., 2012a; Wang et al., 2013).

目前在极/球坐标系中，有限差分法不能较好处理起伏地表边界条件的，所以我们希望可以将牵引力镜像法(Zhang and Chen, 2006)推广到极/球坐标系有限差分中.作为第一步工作，本文集中于将牵引力镜像法推广到二维极坐标系有限差分中处理起伏地表边界条件.

1 数学原理

1.1 规则地表情况下极坐标系方程

(1)

运动学方程和广义胡克定律可以写为

σ=C:ε,

(2)

(3)

1.2 起伏地表情况下极坐标系方程

当存在起伏地形情况下，离散网格和实际地形不一致将引起的虚假波形，因此采用贴体网格贴合实际地形的有限差分(Fornberg, 1988; Zhang and Chen, 2006)具有更高精度.参照前人工作(Zhang et al., 2012a)，以ζ为横向坐标η为径向坐标，引入贴合地形的贴体网格ζ=ζ(r,θ),η=η(r,θ)，将物理空间中不规则的贴体网格映射到计算空间中的规则网格，参见图 1.应用链式求导法则，(1)式中的哈密顿算子和极坐标基矢量关系变为(4)式形式：

(4)

通过(4)式中正交基矢量的微分关系，可以得到坐标基矢量在贴体网格坐标系中的空间导数，参见 (5)式：

图1 物理网格-计算网格示意图通过ζ=ζ(r,θ),η=η(r,θ)将不规则的物理网格(r,θ)映射为规则的计算网格(η,ζ).Fig.1 Physical-computational spaceThe operators ζ=ζ(r,θ),η=η(r,θ) mapping from irregular physical grid (r,θ) to regular computational grid (η,ζ).

(5)

通过r和θ的正交关系，可以得到贴体网格的空间空间变形系数矩阵：

(6)

将(4)、(5)、(6)式代入(2)式，经过化简可以得到贴体网格下弹性动力学方程组：

(7)

1.3 极坐标系下贴体网格中牵引力镜像法

(8)

将(8)式代入牵引力公式后，自由地表处牵引力为0的条件可以写为

Tη=nη·σ=0,

(9)

在对(9)式进行时间域求导后，将(7)式中应力变量时间导数表达式代入，整理后可以得到速度分量在η方向导数使用速度分量、速度分量在ξ方向导数表达的具体形式，参见(10)式.由于体贴网格中所有变量在ξ方向的导数全部可解，所以自由地表处速度分量在η方向导数可以由(10)式计算.

(10)

由于自由地表边界条件此时仅有牵引力关于地表反对称的假设，所以需要把(7)式中所有应力分量η方向的导数通过牵引力及其镜像表达.将(9)式中牵引力表达式对η求导，通过使用(5)式和(6)式化简整理可得(11)式.至此，(7)式中自由地表处所有应力分量η方向导数，可以用组合的形式，由应力分量、牵引力分量及牵引力分量η方向导数表示.

(11)

至此，在自由地表处(7)式中所涉及的自由地表边界条件施加完成.

1.4 其他数值计算处理

仿照Zhang和Chen(2006)，本研究所用空间差分模板采用DRP/opt MacCormack格式(Hixon, 1998)，时间步积分采用4-6阶Runge-Kutta积分(Hu et al.,1996).

对于自由地表下临近格点，如果差分模板跨过自由地表，此时速度/应力分量η方向导数需要额外处理.其中应力分量η方向导数可以继续沿用(11)式，但是速度分量η方向导数无法采用(10)式求解，此处沿用Zhang和Chen(2006)处理方法，采用紧致差分格式(Hixon and Turkel, 2000)，在保证4阶精度情况下计算相关格点速度分量η方向导数.

本算法中震源引入方法与Zhang和Chen(2006)一致，空间上的分布为3个格点半宽度的光滑高斯形，二维分布采用两个一维函数乘积构成以提高精度(Tornberg and Engquist, 2004).

本算法在非自由地表的吸收边界处，采用e指数衰减层处理(Cerjan et al., 1985).

2 算例对比

我们需要计算一系列算例和参考结果对比，验证本算法正确性.考虑到Zhang和Chen(2006)可以在直角坐标系下很好地处理地表起伏，所以我们将其方法结果作为参照解.通过选用相同介质模型和贴体网格，分别使用直角坐标系和极坐标系进行计算，在归一化后对二者进行对比检验.

对于DRP/opt MacCormack格式，传播20到30个波长的距离上稳定性条件为每波长内网格点6～8个(Hixon, 1998).本文选取横波内每波长网格点25，CFL约0.3，以提高计算精度.

我们选择模型为均匀介质，纵波波速为3.0 km·s-1，横波波速为1.2 km·s-1，密度为1800 kg·m-3；震源为爆炸源，震源深度为0.5 km，震源时间函数为Ricker子波，中心频率为0.96 Hz，时间延迟为1.625 s；计算时间步长为0.005 s；空间网格采用贴体网格划分，在横向/θ方向网格间距不大于0.05 km.对地表起伏光滑变化的情况，垂向/R方向采用垂向变换网格(Ruud and Hestholm, 2001)，网格初始间距为0.05 km，根据地表起伏进行一定的拉伸或压缩.

在图2—5中，每幅图的最上方为地表附近网格划分示意图，*代表震源，并用S标识震中位置，A、B则为两台站，方便在震相图中对应识别，空心正三角形为密集布置台站位置，实心倒三角为稀疏布置台站，图中网格示意图为计算网格抽稀后结果，以清晰显示.为方便直观对比波形，我们将极坐标系下波形结果，全部转换到相应直角坐标系下显示.每幅图的中间为密集台站模拟结果对比，以StaVx和StaVz分别标识水平方向和垂直方向震动，可以用于检验本方法波形同相轴的一致性.每幅图的最下方为波形显示比例放大3倍后的稀疏台站的波形对比，以ZoomVx和ZoomVz分别标识水平方向和垂直方向震动，可以精细检验本文方法的振幅和相位.在波形对比图中，红色波形为采用Zhang和Chen(2006)方法在直角坐标系下计算的结果，黑色波形为采用本文结果在极坐标系下的结果.在波形图中，不同台站波形，按照台站沿自由表面到震中距离排列；所有极坐标系下计算，均设定震源正上方处地表半径为100 km，计算网格为扇形，不包括极坐标系原点.

图2 极坐标下规则地表波形对比Fig.2 Waveform comparison in polar coordinate system with regular topography

图3 极坐标系下起伏地表情况波形对比Fig.3 Waveform comparison in polar coordinate system with irregular topography

图4 水平地表情况波形对比Fig.4 Waveform comparison under horizontal surface

图5 水平地表叠加光滑起伏地表情况波形对比Fig.5 Waveform comparison under horizontal surface with irregular topography

首先我们对比规则地表和仅存在局部起伏地表两种情况.图2为本文算法规则地表下退化结果的对比；图3为在规则地表上叠加光滑地形起伏时结果对比，所加山谷深度和山峰高度均为2 km.可以看到在这两种情况下本文算法和参照结果吻合，因此本方法适用于局部光滑地形起伏情况下自由地表边界条件处理.

考虑到本文方法在极坐标系下使用贴体网格，同样可以处理直角坐标系下水平地表情况，此时的水平地表在极坐标系中可以认为是全局光滑变形的起伏地表，对比结果参见图 4.图 5对比了在极坐标系中全局光滑变形基础上，叠加小尺度变形的波形对比结果.本文方法在这些对比中和参照结果吻合，因此本方法适用于全局光滑地形起伏情况下自由地表边界条件处理，同时也适用于全局光滑地形起伏和局部地形起伏耦合情况.

上述各算例中，本文方法计算波形和参考算法结果高度一致，说明本文方法准确、有效，可以精确处理极坐标系有限差分方法中起伏地表处自由边界条件.

3 总结与展望

本文基于极坐标系贴体网格有限差分，引入贴体网格和牵引力镜像方法处理起伏地表边界条件，并和直角坐标系下已有方法对比，算法得到了较全面检验，证明本文方法正确、有效.

本文针对二维极坐标系，利用链式求导法则展开哈密顿算子到贴体网格坐标系中，并利用矢量正交性质化简方程.这些方法全部可以方便地拓展到三维球坐标系，可以预期拓展后的方法对更加精确模拟实际地震波场在地球中传播具有重要意义.

本文研究关注点为起伏地表边界条件，模拟区域均为扇形区域，由于极坐标系固有的周期性，在施加循环边界条件后，本工作模拟区域可以方便地从扇形区域扩展到完整圆环区域；现阶段本文方法未包含极坐标系原点，因此无法模拟穿过原点及附近区域的相关震相，例如天然地震中的PKIKP震相，也需要进一步完善，Wang等(2001)采用的球心处理方法，可以为本文方法完善提供重要参考.

致谢感谢南方科技大学张振国老师对本文完成提供的参考意见和帮助.