地球物理学报  2021, Vol. 64 Issue (2): 546-554   PDF    
极坐标系有限差分中起伏地表边界条件处理
徐剑侠1, 张伟2, 陈晓非2     
1. 中国科学技术大学地球和空间科学学院, 合肥 230026;
2. 南方科技大学地球与空间科学系, 广东深圳 518055
摘要:有限差分算法是地震学中重要的算法,在直角坐标系下同位网格有限差分中使用牵引力镜像方法,可以高效准确地处理起伏地表边界条件.当研究区域-全球尺度问题时需要考虑地球曲率影响,此时选择极坐标系更加直观方便,但已有方法无法在极坐标系下准确计算起伏地表影响.本文在极坐标系有限差分中引入贴体网格和牵引力镜像方法处理起伏地表边界条件,并在多个算例中验证算法正确性和适用范围,证明牵引力镜像法在极坐标系有限差分中有效.
关键词: 有限差分      极坐标系      起伏地表      牵引力镜像法     
Finite difference in the polar coordinate system with irregular topography
XU JianXia1, ZHANG Wei2, CHEN XiaoFei2     
1. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. Department of earth and space science, Southern University of Science and Technology, Shenzhen Guangdong 518055, China
Abstract: The finite difference method (FDM) is an important numerical tool in seismological research. The application of the traction free surface boundary condition is a key issue in the FDM. For the Cartesian coordinate system FDM, when irregular topography exists, the traction image method can satisfy the traction free condition with high accuracy and high efficiency. For seismic wave modeling on regional or global scales, the effect of Earth's curvature must be considered. In this case, the polar coordinator system FDM is more convenience, but the existing method cannot calculate the effect of topography accurately. In this work we introduce a boundary confirm grid and a traction image method into the polar coordinate system FDM to solve this problem. Several numerical tests confirm the validity of our work.
Keywords: Finite difference method    Polar coordinate system    Irregular topography    Traction image method    
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)

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

(2)

其中ρ=ρ(r, θ)是介质密度,是质点运动速度,是二阶应力张量,是二阶应变张量,是四阶弹性常数张量,此处给出的是在各向同性介质中的表达式,λμ是拉梅系数(尹祥础,2011).当把(1)式代入(2)式之后,经过化简可以得到规则地表下均匀介质中的弹性动力学方程组.

(3)

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

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

(4)

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

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

(5)

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

(6)

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

(7)

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

自由地表处边界条件为牵引力T等于0,即T=n·σ=0,其中n为自由地表处垂直方向单位法向量.牵引力镜像方法假定牵引力关于自由地表反对称,从而保证自由地表处牵引力为0.对于规则网格,,此时牵引力镜像方法退化为应力镜像法.对于非均匀网格,自由地表处垂直方向单位矢量可以写成:

(8)

将(8)式代入牵引力公式后,自由地表处牵引力为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,根据地表起伏进行一定的拉伸或压缩.

图 25中,每幅图的最上方为地表附近网格划分示意图,*代表震源,并用S标识震中位置,A、B则为两台站,方便在震相图中对应识别,空心正三角形为密集布置台站位置,实心倒三角为稀疏布置台站,图中网格示意图为计算网格抽稀后结果,以清晰显示.为方便直观对比波形,我们将极坐标系下波形结果,全部转换到相应直角坐标系下显示.每幅图的中间为密集台站模拟结果对比,以Sta Vx和Sta Vz分别标识水平方向和垂直方向震动,可以用于检验本方法波形同相轴的一致性.每幅图的最下方为波形显示比例放大3倍后的稀疏台站的波形对比,以Zoom Vx和Zoom Vz分别标识水平方向和垂直方向震动,可以精细检验本文方法的振幅和相位.在波形对比图中,红色波形为采用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)采用的球心处理方法,可以为本文方法完善提供重要参考.

致谢  感谢南方科技大学张振国老师对本文完成提供的参考意见和帮助.
References
Alterman Z, Karal F C Jr. 1968. Propagation of elastic waves in layered media by finite difference methods. Bulletin of the Seismological Society of America, 58(1): 367-398.
Appel D, Petersson N A. 2009. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in Computational Physics, 5(1): 84-107.
Bayliss A, Jordan K E, LeMesurier B J, et al. 1986. A fourth-order accurate finite-difference scheme for the computation of elastic waves. Bulletin of the Seismological Society of America, 76(4): 1115-1132.
Cerjan C, Kosloff D, Kosloff R, et al. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics, 50(4): 705-708. DOI:10.1190/1.1441945
Fornberg B. 1988. The pseudospectral method:accurate representation of interfaces in elastic wave calculations. Geophysics, 53(5): 625-637. DOI:10.1190/1.1442497
Hestholm S, Ruud B. 1994. 2D finite-difference elastic wave modelling including surface topography. Geophysical Prospecting, 42(5): 371-390. DOI:10.1111/j.1365-2478.1994.tb00216.x
Hestholm S, Ruud B. 1998. 3-D finite-difference elastic wave modeling including surface topography. Geophysics, 63(2): 613-622. DOI:10.1190/1.1444360
Hixon R. 1998. Evaluation of a high-accuracy MacCormack-Type scheme using benchmark problems. Journal of Computational Acoustics, 6(3): 291-305. DOI:10.1142/S0218396X9800020X
Hixon R, Turkel E. 2000. Compact implicit MacCormack-type schemes with high accuracy. Journal of Computational Physics, 158(1): 51-70. DOI:10.1006/jcph.1999.6406
Hu F Q, Hussaini M Y, Manthey J L. 1996. Low-dissipation and low-dispersion runge-kutta schemes for computational acoustics. Journal of Computational Physics, 124(1): 177-191. DOI:10.1006/jcph.1996.0052
Igel H, Weber M. 1995. SH-wave propagation in the whole mantle using high-order finite differences. Geophysical Research Letters, 22(6): 731-734. DOI:10.1029/95GL00312
Käser M, Igel H, Sambridge M, et al. 2001. A comparative study of explicit differential operators on arbitrary grids. Journal of Computational Acoustics, 9(3): 1111-1125. DOI:10.1142/S0218396X01000838
Lan H Q, Zhang Z J. 2011. Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface. Bulletin of the Seismological Society of America, 101(3): 1354-1370. DOI:10.1785/0120100194
Lan H Q, Liu J, Bai Z M. 2011. Wave-field simulation in VTI media with irregular free surface. Chinese Journal of Geophysics (in Chinese), 54(8): 2072-2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
Lan H Q, Zhang Z J. 2012. Seismic wavefield modeling in media with fluid-filled fractures and surface topography. Applied Geophysics, 9(3): 301-312. DOI:10.1007/s11770-012-0341-5
Li D Z, Helmberger D, Clayton R W, et al. 2014. Global synthetic seismograms using a 2-D finite-difference method. Geophysical Journal International, 197(2): 1166-1183. DOI:10.1093/gji/ggu050
Liu X B, Chen J Y, Zhao Z C, et al. 2018. Simulating seismic wave propagation in viscoelastic media with an irregular free surface. Pure and Applied Geophysics, 175(10): 3419-3439. DOI:10.1007/s00024-018-1879-9
Lombard B, Piraux J, Gélis C, et al. 2008. Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves. Geophysical Journal International, 172(1): 252-261. DOI:10.1111/j.1365-246X.2007.03620.x
Mulder W A. 2017. A simple finite-difference scheme for handling topography with the second-order wave equation. Geophysics, 82(3): T111-T120. DOI:10.1190/geo2016-0212.1
Mulder W A, Huiskes M J. 2017. A simple finite-difference scheme for handling topography with the first-order wave equation. Geophysical Journal International, 210(1): 482-499. DOI:10.1093/gji/ggx178
Nilsson S, Petersson N A, Sj green B, et al. 2007. Stable difference approximations for the elastic wave equation in second order formulation. SIAM Journal on Numerical Analysis, 45(5): 1902-1936. DOI:10.1137/060663520
Rao Y, Wang Y H. 2013. Seismic waveform simulation with pseudo-orthogonal grids for irregular topographic models. Geophysical Journal International, 194(3): 1778-1788. DOI:10.1093/gji/ggt190
Rao Y, Wang Y H. 2018. Seismic waveform simulation for models with fluctuating interfaces. Scientific Reports, 8(1): 3098. DOI:10.1038/s41598-018-20992-z
Robertsson J O A. 1996. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics, 61(6): 1921-1934. DOI:10.1190/1.1444107
Ruud B, Hestholm S. 2001. 2D surface topography boundary conditions in seismic wave modelling. Geophysical Prospecting, 49(4): 445-460. DOI:10.1046/j.1365-2478.2001.00268.x
Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics, 69(2): 583-591. DOI:10.1190/1.1707078
Solano C A P, Donno D, Chauris H. 2016. Finite-difference strategy for elastic wave modelling on curved staggered grids. Computational Geosciences, 20(1): 245-264. DOI:10.1007/s10596-016-9561-8
Tessmer E, Kosloff D, Behle A. 1992. Elastic wave propagation simulation in the presence of surface topography. Geophysical Journal International, 108(2): 621-632. DOI:10.1111/j.1365-246X.1992.tb04641.x
Tornberg A K, Engquist B. 2004. Numerical approximations of singular source terms in differential equations. Journal of Computational Physics, 200(2): 462-488. DOI:10.1016/j.jcp.2004.04.011
Toyokuni G, Takenaka H. 2006. FDM computation of seismic wavefield for an axisymmetric earth with a moment tensor point source. Earth, Planets and Space, 58(8): e29-e32. DOI:10.1186/BF03352593
Virieux J. 1986. P-SV wave propagation in heterogeneous media; velocity-stress finite-difference method. Geophysics, 51(4): 889-901. DOI:10.1190/1.1442147
Wang X M, Zhang H L. 2004. Modeling of elastic wave propagation on a curved free surface using an improved finite-difference algorithm. Science in China Series G:Physics, Mechanics and Astronomy, 47(5): 633-648. DOI:10.1360/03yw0253
Wang Y B, Takenaka H, Furumura T. 2001. Modelling seismic wave propagation in a two dimensional cylindrical whole earth model using the pseudospectral method. Geophysical Journal International, 145(3): 689-708. DOI:10.1046/j.1365-246x.2001.01413.x
Wang Y B, Takenaka H, Jiang X H, et al. 2013. Modelling two-dimensional global seismic wave propagation in a laterally heterogeneous whole-moon model. Geophysical Journal International, 192(3): 1271-1287. DOI:10.1093/gji/ggs094
Yin X C. 2011. Solid Mechanics (in Chinese). 2nd ed. Beijing: Seismological Press.
Zhang W, Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. Geophysical Journal International, 167(1): 337-353. DOI:10.1111/j.1365-246X.2006.03113.x
Zhang W, Shen Y, Zhao L. 2012a. Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method. Geophysical Journal International, 188(3): 1359-1381. DOI:10.1111/j.1365-246X.2011.05331.x
Zhang W, Zhang Z G, Chen X F. 2012b. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids. Geophysical Journal International, 190(1): 358-378. DOI:10.1111/j.1365-246X.2012.05472.x
兰海强, 刘佳, 白志明. 2011. VTI介质起伏地表地震波场模拟. 地球物理学报, 54(8): 2072-2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
尹祥础. 2011. 固体力学. 2版. 北京: 地震出版社.