地球物理学报  2019, Vol. 62 Issue (11): 4339-4352   PDF    
基于辛-谱元-FK混合方法的保结构远震波场模拟
李冰非1,2,3, 董兴朋4, 李小凡1,2, 司洁戈5     
1. 中国科学院地质与地球物理研究所, 北京 100029;
2. 中国科学院大学, 北京 100049;
3. 中国地震灾害防御中心, 北京 100029;
4. 清华大学数学科学系, 北京 100084;
5. 中国电子科技集团公司第三研究所, 北京 100015
摘要:远震全波形层析成像能获得研究区域下方岩石圈乃至地幔过渡带高分辨率速度结构,是研究地球深部构造与动力学过程的有效工具.该类方法需以高精度及长时程远震波场正演模拟为基础,这为设计高精度长时程稳定的正演算法带来了挑战.在此背景之下,本文提出了一种适用于远震波场模拟的保结构算法.该方法采用谱元法(SEM)对研究区域进行空间离散,在不考虑耗散项情况下,将空间离散后的常微分方程变换为哈密顿系统形式,采用保辛分部龙格-库塔方法数值求解.在三级保辛分部龙格-库塔算法基础上添加额外空间离散项,得到修正辛算法.本文将该时间-空间全离散形式称为修正辛-谱元法(SSEM),并将SSEM算法与频率波数域(FK)方法结合,发展了可模拟高频远震波场在局域模型内传播的SSEM-FK混合方法.该方法结合了FK方法模拟层状介质中平面波传播的高效性和SSEM计算复杂介质中弹性波传播的精确性.数值实验表明,SSEM-FK能够准确模拟高频远震波场在研究区域内的传播,结合该方法在计算效率上的优势,可为高效、高精度的远震全波形层析成像打下基础.
关键词: 辛算法      弹性波方程      远震波场模拟      谱元法     
Structure-preserving modeling of teleseismic wavefield using symplectic SEM-FK hybrid method
LI BingFei1,2,3, DONG XingPeng4, LI XiaoFan1,2, SI JieGe5     
1. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. China Earthquake Disaster Prevention Center, Beijing 100029, China;
4. Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China;
5. The 3rdResearch Institute of CETC, Beijing 100015, China
Abstract: Teleseismic full-waveform tomography can obtain high resolution velocity structure of upper mantle, and even the structure of mantle transition zone beneath the study area. It has been a vital tool in probing the Earth's internal structure and dynamic process. This method is based on effective forward simulation of seismic wavefield, which requires forward-modeling methods with high accuracy and long-term stability. Under this background, we present a structure-preserving scheme based on symplectic geometric theory, which is specially tailored for teleseismic wave modeling. The spectral element method (SEM) is used to discretize the study area. Without the consideration of the dissipation term, the discretized ordinary differential equation is transformed into a Hamiltonian system, which can be solved numerically by the symplectic-preserving partitioned Runge-Kutta (PRK) method. After adding an extra spatial discretization term into the PRK scheme, we developed a new symplectic scheme with fourth-order temporal accuracy. We call the spatial-temporal discretized scheme symplectic SEM (SSEM). We combine SSEM with FK method, which calculates plane wave propagation in one-dimensional layered media. After that, we construct a hybrid method (SSEM-FK) to efficiently calculate teleseismic elastic wave propagation in regional heterogeneous model. Numerical experiments show that SSEM-FK can accurately simulate the propagation of high-frequency teleseismic wavefield in the study area. The advantage of this method in computational efficiency forms the foundation for high-efficient and high-accurate teleseismic full-waveform tomography.
Keywords: Symplectic scheme    Elastic wave equation    Teleseismic wavefield modeling    Spectral element method    
0 引言

地震层析成像是研究地球内部结构、揭示地球深部动力学过程的重要手段(Rawlinson et al., 2010; Liu and Gu, 2012).由于大部分地震发生在地壳,区域地震层析成像较难获得岩石圈深部构造信息(Tian et al., 2009).经过地幔向上传播被研究区域台站记录得到的远震波形可有效用于研究地壳和上地幔构造(Rondenay, 2009; Liu et al., 2017a, 2018), 如:接收函数可利用远震地震波在速度不连续介质中发生的转换波研究地下的速度间断面(Langston, 1977; Kind et al., 2012),但接收函数难以反演岩石圈横向非均匀性(Langston, 1977).采用广义Radon变换方法,Bostock等(2001)使用地震台站接收的远震尾波数据获得了研究区域下方的弹性参数扰动.尽管如此,他们的方法仍基于射线理论,难以对地下复杂介质高分辨成像(Tape et al., 2009).地壳-上地幔的高分辨成像迫切需要发展基于伴随原理的远震波形成像方法(Tong et al., 2014a).

能有效模拟高频远震波形的数值模拟方法是高分辨率远震波形成像的关键.对于远震波场在区域模型中传播的数值模拟,需要模拟的高频地震波形(0.2~2 Hz)时长往往超过200 s (Liu et al., 2017a, 2017b),因此有必要构造高效且具有长时间数值稳定能力的地震波模拟算法(Nissen-Meyer et al., 2008; Chen, 2009a; Li et al., 2012; Cai et al., 2016, 2017a, 2017b; Zhang et al., 2009).地震波场数值模拟算法涉及到时间离散和空间离散.对于地震波运动方程时间离散,多采用二阶差分格式(Virieux, 1984, 1986; Bohlen and Wittkamp, 2016; Wang et al., 2015).二阶差分方法虽然简单容易实现、计算速度快,但数值频散较为严重(Liu et al., 2014; Tan and Huang, 2014).除二阶差分格式之外,另一种常用的时间离散格式为显式Newmark算法(Komatitsch and Vilotte, 1998; Liu et al., 2017c),该格式计算效率高,且在空间离散格式精度足够高的情形下,具有保能量的优点(Liu et al., 2017c),但数值实验表明该方法在用于地震波长时间模拟时积累误差及数值频散仍较大(Liu et al, 2014a, b).

在无耗散的假设下,地震波方程可表示为Hamilton系统(Qin and Zhang, 1990).在数值求解Hamilton系统时,辛算法可保持Hamilton系统辛结构守恒(Feng and Qin, 2010).对于辛算法在波动方程的求解方面,Qin和Zhang(1990)首先做了尝试,其在空间上采用有限差分离散,利用显式辛算法对一维声波方程时间离散.随后,孙耿(1997)证明了在空间离散算子对称的条件下,离散方程的辛守恒性将独立于空间离散.在隐式辛算法设计方面,罗明秋等(2001)构造了基于Padé逼近的二阶隐式辛格式,并用于声波方程的数值求解.隐式辛格式需要求解大规模线性系统,计算量较大,计算效率较低.针对计算效率低的问题,罗明秋等(2003)提出了基于查表法的改进算法.但相比于显式算法,改进算法计算效率仍然较低(Qin and Zhang, 1990; Li et al., 2011).显式算法的高效性吸引了Chen(2007)的注意,Chen(2009a)在系统对比了三种高阶时间格式之后,得出Lax-Wendroff格式时间误差累积较大的根本原因在于Lax-Wendroff格式仅为二阶精度的辛算法,高阶辛算法误差压制能力明显高于低阶辛算法.最近,Liu等(2017c)研究表明中心差分方法与一阶辛算法等价,其研究也说明了中心差分难以运用于高精度、长时间地震波场模拟.关于辛算法与SEM的结合得益于Nissen-Meyer等(2008)的研究,其研究表明四阶辛算法较适用于模拟大尺度模型中的地震波传播.

在常规辛算法中,高阶格式(如四阶)虽然其精度较高,但计算量较大(Chen, 2009a),低阶格式(如二阶)计算量虽较小,但精度较差(刘少林等, 2013b, 2015).相对而言,三阶格式计算量较为适中,且精度也较高(Li et al., 2011).为了进一步提高辛格式的计算效率,Liu等(2015)首次提出适用于地震波场模拟的修正辛算法.其在空间采用有限差分(FDM)离散的条件下,讨论了修正辛算法应用于弹性波模拟的计算效率,研究表明修正辛算法在数值频散压制和稳定性方面优于常规辛算法.前人发展的修正辛算法精度仅为三阶,且选择的空间离散算子FDM难以适用于复杂地质模型.鉴于此,本文构造了四阶修正辛算法,在提高精度的同时亦可降低计算量;本文采用谱元法对地震波方程空间离散,得到适用于复杂模型中地震波模拟的辛谱元法(SSEM).随后,本文采用FK方法构造远震入射边界条件(Liu et al., 2017a),最后得到适用于保结构远震波场模拟的SSEM-FK方法.本文利用SSEM-FK模拟远震波场在区域模型中的传播,并模拟了华北岩石圈模型中远震波场传播特征.

1 数值求解弹性波方程的SSEM

区域尺度以及全球尺度地震层析成像往往基于射线理论(Zhao et al., 2012).为了提高成像分辨率,研究者逐渐转向采用伴随层析成像方法(Tape et al., 2009; Tong, et al., 2014b).伴随层析成像方法需要迭代求解波动方程,因此需要巨大的计算资源.对于高分辨地震成像问题而言,发展适用于地震波模拟的高效算法至关重要.

在不考虑外力的情形下,各向同性介质中弹性波方程写为

(1)

其中u为位移,c为四阶弹性张量,f为外力.对于大尺度地震波传播与层析成像问题, 较为流行空间离散方式是SEM(Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 2002; De Basabe and Sen Mrinal, 2007; 汪文帅等, 2012).为了适用于复杂模型中地震波场模拟,本文亦选择SEM对弹性波运动方程空间离散.在不考虑外力的情形下,离散空间域后的微分方程组可写为(Komatitsch and Vilotte, 1998)

(2)

其中M为全局质量矩阵,K为全局刚度矩阵,U表示全局位移向量,引入变量(2)式等价为

(3)

其中L =- M-1 K为空间导数算子,V为速度向量.(3)式可表示为Hamilton形式

(4)

其中哈密顿函数H

(5)

(5) 式表明(4)式为可分的哈密顿系统.应用辛PRK方法(Liu et al., 2017b)可得如下方程:

(6)

其中Δt为时间步长,s为辛格式级数,UnVn的上标表示第n时间步,UiVi为中间变量,aijbiAijBi为辛系数.若(6)式为辛格式,则数值离散以后方程辛性守恒(Feng and Qin, 2010),即

(7)

其中∧表示外积(Feng and Qin, 2010).根据(6)和(7)式,辛系数应满足如下等式:

(8)

根据(6)和(8)式,显式三级(s=3) PRK格式可以写为如下简洁形式(Liu et al., 2017b):

(9)

可以运用根数理论(Feng and Qin, 2010)或者泰勒展开得到(9)式的阶条件方程组.Liu等(2014a)证明了显式二阶辛PRK格式只能达到二阶精度.为了提升二阶辛PRK格式的精度,常用的两种方法有频率误差最小化方法(刘少林等, 2015)和截断误差最小化方法(McLachlan and Atela, 1992; Liu et al., 2014a).然而二阶精度难以满足高精度地震波场模拟的要求(Liu et al., 2017b).为了获得高阶时间精度,本文引入空间导数项用于修正原格式,并获得更高阶精度,由(9)式可得

(10)

在(10)式中,我们加入了Δt3 LV1项修正原格式的精度.为了使格式(10)具有三阶时间精度,本文通过Taylor展开得到系数B1, b1, B2, b2, B3B4的限定方程(具体见附录A).三阶辛系数满足的阶条件为

(11)

上式中可用B3表示其他系数,即

(12)

为获得对称系数,我们设b1=b2,于是得到以下系数:

(13)

可以证明(13)式中的系数可以保证(10)式满足辛格式(见附录B).事实上在对称系数(13)式的条件下(10)式具有四阶精度.值得注意的是,(13)式中的系数全为正实数,从理论上可以证明收敛性及稳定性.传统辛PRK和RKN方法当精度为三阶或以上时总是出现负的辛系数(McLachlan and Atela, 1992).负系数在时间积分上是不稳定的(Chin and Anisimov, 2006; Liu et al., 2017b),而修正辛算法中辛系数均为正数.本文中(10)式适用于弹性波动方程,对于带耗散的地震波运动波动方程,可以算子分裂为一个守恒系统和一个耗散系统(Cai et al., 2017b).守恒系统可用传统辛格式或本文修正辛算法求解,耗散系统可以解析求解, 最终带有耗散项的系统可以应用Strang方法将两部分解结合得到带衰减的地震波运动方程的数值解(Cai et al., 2017b).

2 算法理论分析

精度和稳定性是评价数值方法的重要指标.本文应用Fourier分析方法(Yang et al., 2006)和矩阵本征值理论(Liu et al., 2015)分析SSEM格式的稳定性和频散.为了对比的需要,本文引入常用四种方法对比.将本文SSEM记为M1,二阶辛PRK-SEM记为M2(刘少林等, 2013a, b),三级辛PRK- SEM记为M3(Li et al., 2011),四阶Lax-Wendroff Correction (Deablain, 1986)记为M4,二阶中心差分方法记为M5(Di Bartolo et al., 2012).事实上M2等价于无衰减的显式Newmark方法(Liu et al., 2017c),故M1不再与经典Newmark方法(Hashamder et al., 2011)对比.

2.1 稳定性分析

n时刻第j个网格结点上的谐波解为

(14)

其中Δt为时间步长,ωnum为频率,MN为任意常向量,下标j表示空间离散节点.将(14)式代入(10)式,第(n+1)时刻数值解可以写为.转换矩阵A见附录B.由(14)式,得到如下特征问题:

(15)

稳定性条件通过求解不等式ρ(A) < 1来获得,其中ρ(A)表示A的谱半径, ρ(A)=max |λ|, λA的特征值.由于矩阵A较为复杂,稳定性条件的解析表达式难以解析获得.在这里,我们将矩阵A转换为A的相似矩阵D (Liu et al., 2015),即

(16)

其中相似矩阵D可表示为

(17)

其中I为单位矩阵,Σ为空间离散矩阵.显然,ρ(D)=ρ(A).若用γ表示矩阵D的特征值,λ表示Σ的特征值,则γλ满足如下方程:

(18)

ρ(D) < 1可等价为

(19)

(19) 式可等价简化为

(20)

从(19)式中的多项式可知时间离散和空间离散耦合,(20)式中时间格式确定后,λj仅与空间离散有关,(20)式中的数值12又称为时间放大因子(temporal amplification factor, TAF) (Liu et al., 2015).表 1列出了五种时间格式的TAF.从表 1可以看出,M1和M4方法TAF值最大,M2和M5方法TAF值最小.

表 1 五种时间格式的TAF值 Table 1 TAF values of five temporal schemes

根据(20)式,λj与空间离散有关.无限空间下求矩阵L的特征值较难,空间采用规则的正方形单元可将SEM空间离散的无限特征值问题化简为有限维问题(De Basabe and Sen, 2007; Liu et al., 2017a, 2017b).根据Liu等(2015)的研究,(10)式的稳定性条件可写为

(21)

其中VP表示P波速度,r=VP/VS为纵横波速度比,h0=h/N(N阶插值的平均空间间隔),ab为待定常数.表 2给出N=1时ab的值.表 3给出N≥2时b的值(此时a为零).

表 2N=1时M1—M3的稳定性参数 Table 2 Stability parameters for M1—M3 with N=1
表 3 空间为不同插值阶数时M1—M3的稳定性参数 Table 3 Stability parameters for M1—M3 with N≥2

由于a只与空间离散相关,当SEM空间插值阶数为1时,表 2a的值都为1.方法的稳定性强弱取决于b的值.从表 2表 3可知M1的b值最大,故该方法的稳定性优于其他两种方法的.由表 3可知,随着阶数N的增加,所有方法的稳定性都降低,但M1的稳定性依然最强.

2.2 频散分析

数值频散是影响数值计算精度的重要因素,本小节中讨论M1,M2和M3三种方法的数值频散.谱元法在空间上往往采用高阶插值多项式以提高计算精度(Cohen, 2002),De Basabe和Sen Mrinal(2007)证明了高阶谱元法的弱数值频散特性.为了分析仅由时间离散所引入的数值频散,本文空间采用高阶插值(6阶)来压制空间离散而引入的数值频散.在相同频率下,P波的波长要比S波的大,因此在相同空间网格的情形下,一个波长内P波的采样点数要大于S波的,同时值得注意的是,在相同采样的条件下S波的频散往往大于P波的(De Basabe and Sen Mrinal, 2007; 刘少林等, 2014).因此在频散分析中,S波的频散能更好地评价数值方法的整体数值频散.

(15) 式中exp(iωnumΔt)=γ,其中γA的特征值.可推导M1的数值频散关系为

(22)

其中VSnumVS分别为S波的数值波速和物理波速,α=VPΔt/h为库朗数;采样率G定义为G=h0/λ(λ表示波长),γχ的关系式为χt2γ/h02.图 1为四个不同库朗数条件下泊松介质中S波数值频散曲线.当G>1/(2N)时,特征值变化极不规则,因此本文只给出G < 1/(2N)时的数值频散曲线(Liu et al., 2017b).

图 1显示0≤G≤1/N时的频散曲线(N=6).由图 1可知,M2的数值频散最为严重,相对速度误差远大于M1和M3.为更好地对比M1和M3的数值频散, 在图 1(a—d)中分别将误差放大5000, 5000, 5000和4000倍.从图 1可看出尽管M1有更高阶精度(时间精度为4阶),但其数值频散要稍大于M3的.这些对比结果表明:数值频散大小由增长矩阵A的迹决定,而数值精度取决于增长矩阵A的非对角元素,频散和精度之间并没有必然的联系.当M1和M3数值频散放大4000或5000倍时才显示微弱区别,即说明M1和M3几乎具有相同的数值频散特性.

图 1 M1—M3及参考格式的S波数值频散曲线 谱元法采用6阶空间插值.纵横波速度比为图(a—d)分别对应库朗数为0.2, 0.4, 0.5和0.6.图(a—c)中将除M2外所有方法误差放大5000倍.图(d)中将除M2外所有方法误差放大4000倍. Fig. 1 Relative velocity errors of S-wave numerical dispersions SEM with sixth-order interpolation is used in the space. The P- to S-wave velocity ratio is Plots (a—d) correspond to the four Courant numbers 0.2, 0.4, 0.5 and 0.6, respectively. For visualization, the results of all the methods in (a—c), except for that of M2, are multiplied by 5000; the results of all the methods in plot (d), except for that of M2, are multiplied by 4000.

值得注意是,图 1中由于参考算法不含时间域频散, 其相对速度误差应该低于三种辛格式, 但在图 1a中采样率大于0.06时M3数值频散最小,这种现象的原因是空间和时间离散耦合,在某些情况下时空数值频散可以相互抵消.由于采样率最大仅为1/12,所以库朗数为0.6时M2和M3依然稳定.

3 FK-SEM混合方法的数值传递与散射波吸收

自从Aki等的开创性工作以来(Aki and Lee, 1976),地震层析成像已成为研究地球内部结构的主要方法.根据震源是否在研究区域内,地震层析成像法可分为区域内层析成像和远震层析成像,对于地震活动较少的区域可利用远震数据进行成像(Monteiller et al., 2013, 2015).当远震波入射到局部研究区域时, 远震波可以近似为平面波(Tong et al., 2014a).本节介绍SEM方法和FK方法在边界结合处的数值传递原理.如图 2所示,远震激发的波场以平面波形式入射到研究区域下方,入射角为θ,灰色阴影区T是我们要研究的目标区域,其嵌入在水平背景模型之中.区域T内是非均匀介质,可能包含起伏的不连续面(如Moho面)或高/低速的速度异常体等.远震平面波在一维水平层状区域中的传播可以采用FK方法计算(Haskell, 1953),然而当远震平面波传播到区域T时,由于T的非均匀性,FK方法不再有效(Tong et al., 2014a),因此需要用数值方法如SEM来模拟远震波场在区域T中的传播.在两种算法的交界处,需要进行数值传递,该过程在数据交换域中进行.数值交换域(图 3)具有如下的两个作用:

图 2 远震波场以平面波形式入射到研究区域下方,入射角为θ.FK方法得到水平层状介质中的半解析解,研究区域T内结构复杂(如包含Moho起伏面或高/低速异常体等) Fig. 2 The teleseismic wavefield is incident to the lower part of the study area, and the incident angle is θ. The FK method gives the semi-analytical solution in the horizontal layered medium. The complex structures, including Moho undulation or high/low velocity zone, are in the study domain T
图 3 目标研究区域和PML区域之间的数据交换区 Fig. 3 Data-exchange area between study domain and PML domain

(1) 在每个时刻,用FK方法计算得到的远震波场为SEM方法计算的区域T提供边界条件;

(2) 将区域T内由于非均匀性产生的散射波向外传播,并用吸收边界条件吸收.

我们记目标区域边界最外层为Bc, B,PML与目标区域相邻的最外层为Bp, B,Bc, B处由SEM方法计算得到的总波场为UT,FK方法得到的层状介质中的半解析解场为UFK,则散射波场为:Us= UTUFK,该散射波场需要在Bp, B处应用PML条件吸收,记散射波场Us经过PML吸收后的剩余波场为Us,则在Bc, B处的总波场为

4 数值实验

本节设计三个数值实验用于验证本文构造的SSEM-FK在远震波场模拟时的有效性.

4.1 一维地壳-地幔模型

首先设计一维地壳-地幔模型验证方法的精度.模型大小为100 km×60 km,具体参数设置见表 4.P波入射角设为θ=15°.由于震源和介质速度模型都会影响地震记录, 因此采用准确的震源子波对远震波场模拟极其关键.可以对实际地震数据应用反褶积技术来获取震源时间函数(Rondenay, 2009; Kind et al., 2012).对于合成理论地震图,本文采用截止频率为2 Hz的高斯子波.将模型剖分为960个边长为2500 m的正方形单元,空间采用六阶插值.在初始t=0时刻P波还未到达计算区域(它的波前与莫霍面交于x=-350 km).在t=7.6 s时刻, 平面波的波前到达计算区域的左下角.图 4为三种辛格式与SEM结合模拟一维层状模型中远震传播16 s和20 s时速度水平分量的波场快照.入射P波首先在莫霍面分为上行和下行P波和S波.随后,上行P波在地球表面发生反射和转换.这些反射波或转换波可能经莫霍面反射后再次传播至地表.位于地表的接收点会记录到直达P波以及地表和莫霍面界面产生的透射波、转换波和反射波等多次波.

表 4 一维地壳-地幔模型参数 Table 4 Material parameters for the 1-D crust-mantle model
图 4 三种方法模拟一维模型中远震波场t=16 s和20 s时速度水平分量的波场快照 Fig. 4 Snapshots of velocity wavefields computed by the three symplectic schemes in 1-D crust-mantle model. The number in each plot indicates the time instants

为了验证方法的准确性,设置接收器位于R(60 km, 0 km), 图 5为采用FK和M1-M3模拟远震波传播80 s得到的合成地震记录.从图 5可以看出M1—M3都可以较正确模拟远震波传播.图中显示的远震震相可以借助射线理论来标注(Tong et al., 2014a).图 5(a, b)振幅最强震相是直达Pp震相,之后第二个到达的是Ps震相,接下来到达的震相包括PpPmpPmp震相等.为了检验这些方法的长时间计算能力,图 5c5d分别给出了图 5a5b中的合成波形与解析波形之间的数值残差.可以看出,M2合成的波形误差最大,而M1与M3误差较小,且两者接近.

图 5 采用FK和M1—M3模拟远震波传播80 s得到的合成地震记录对比 (a)接收器水平速度分量;(b)接收器垂直速度分量;(c)和(d)分别为图(a)和(b)矩形框中的波形放大图.(a)和(b)中远震震相a:Pp; b: Ps; c: PpPmP; d: PpPms, PsPmP, PpSmp和e: PsSmp, PsPms, PpSms.(c)和(d)中数字位置显示M2误差远大于M1和M3. Fig. 5 Synthetic seismograms computed by FK and M1—M3 (a) Horizontal-velocity components at the receiver R; (b) Vertical-velocity components at the receiver R; (c) The error of the waveforms generated by M1—M3 in the rectangular box in (a); (d) The error of the waveforms generated by M1—M3 in the rectangular box in (b). The teleseismic phases in (a) and (b) are labeled with letters. a: Pp; b: Ps; c: PpPmP; d: PpPms, PsPmP, PpSmp; and e: PsSmp, PsPms, PpSms.The numbers in (c) and (d) indicate that the error of M2 are much larger than those of M1 and M3.
4.2 含起伏界面的地壳-地幔模型

真实地球模型中Moho不连续面并不水平,而是存在起伏,为了测试M1方法对于起伏Moho面的有效性,我们在4.1节模型基础上,在Moho面上加入两个弧形起伏(如图 6所示),使得模型更接近实际地下起伏结构.起伏界面高差约2 km.震源时间函数与上一算例相同,P波入射角仍设为θ=15°.

图 6 含起伏莫霍面的地壳地幔模型 Fig. 6 Crust-mantle model with undulating Moho:configuration and parameters

选取三个时间步长: 0.01 s, 0.013 s, 0.015 s用于检验上述三种辛方法对于含起伏Moho面模型的有效性.M1,M2和M3能取到的最大时间步长分别为0.015 s、0.01 s和0.013 s.模拟结果如图 7所示.在t=15 s时刻, 波前刚到达莫霍面左侧的起伏处,此时的波前面相比于水平莫霍面波前面(图 4)的轻微起伏说明莫霍面左侧的起伏构造加快了地震波传播.在t=16.5 s时刻, 不仅可以看到波前小的起伏, 还可以看到明显的莫霍面反射波.在t=18 s时刻, 波前到达右侧起伏的莫霍面.在t=19.5 s, 莫霍面的反射波传到模型底部, 可以看到散射波几乎被完全吸收了.值得注意的是人工边界未产生明显虚假反射, 说明本文远震模拟混合方法的吸收边界条件能很好地吸收模型内部产生的散射波.

图 7 三种方法模拟二维复杂模型中地震波速度波场水平分量快照 Fig. 7 Snapshots of the horizontal velocity wavefield computed by the three symplectic schemes in the 2-D complex model

随着时间步长的增大, M2和M3面临失稳.其中当时间步长取0.013 s时M2失稳, 当时间步长取0.015 s时M3失稳.当M2和M3失稳时,M1仍然稳定,说明M1的稳定性优于M2和M3.数值模拟时选择强稳定性的算法意味着可以选择较大的时间步长来提高计算效率.M1的强稳定性对于远震波数值模拟具有明显优势.

4.3 中国华北地区岩石圈模型

华北克拉通(NCC)是在早中生代被破坏的古代克拉通,是最明显和最典型的破坏的克拉通.东部华北克拉通跨华北造山带和西部华北克拉通是该地区的三大构造单元.华北克拉通自形成以来至中生代早期一直保持相对稳定.然而,华北克拉通发生了大规模的构造变形和岩浆活动,特别是在华北克拉通的东部,原克拉通的结构和性质受到严重破坏.岩石圈减薄是华北克拉通破坏的重要标志(Gao et al., 2009).根据密集地震台阵数据获得的最新地震成像结果(Chen, 2009b; Zheng et al., 2006; Zhao et al., 2012), 不同区域华北克拉通的地壳和岩石圈厚度存在明显的差异.克拉通东部广泛分布着较薄的地壳(< 35 km)和岩石圈(60~100 km),而位于克拉通西部的鄂尔多斯盆地岩石圈厚度约为200 km.这一巨大差异表明,华北克拉通的破坏主要发生在太行山以东.

我们基于最新的地震学结果建立了一个深度超过华北克拉通岩石圈-软流圈边界(LAB)的正演模型.根据Chen(2009b)我们建立了一个分层模型,由地壳、上地幔、起伏的莫霍面和LAB组成(图 8),介质参数参考iasp91全球模型.模型跨越华北克拉通的东部和西部连接区域,且LAB有一个快速变化区,最大的深度差异达30km.这一变化符合华北克拉通独特的地质特征:东部被破坏,西部保持稳定.设P波入射角为θ=15°.为了与实际地震数据频带范围一致,我们将截断频率设置为0.5 Hz, 以5 km的间隔共放置41个接收器.

图 8 华北克拉通分层模型 Fig. 8 Layered model for the Central North China Craton

图 9给出了M1得到的速度垂直分量地震记录.一般来说,Pp震相最先到达并具有最强的能量;随后是Ps震相,其能量弱于Pp震相;最后是来自LAB的转换波,将其记为Ps1, 能量最弱.模型右侧处Ps1震相较早到达,对应该处下方较薄的岩石圈.Ps震相也表现出相同的特征,但不如Ps1震相明显.进一步地,我们可以利用正演模拟得到的理论波形计算该区域的理论接收函数,并与实际接收函数进行对比,这可以用来验证由实际数据得到的该区域下方不连续面位置是否可靠.

图 9 M1得到的速度垂直分量地震记录 Fig. 9 Velocity seismograms of vertical component generated by the M1 scheme
5 讨论

本文基于辛几何理论提出了一种适合远震波场模拟的SSEM-FK混合方法.在常规辛算法中加入作用于速度场空间离散项,得到了一种时间四阶精度的修正辛格式.相较于传统时间离散算法, 本文方法稳定性更强, 更能有效地压制数值频散.理论分析和数值实验表明修正辛-谱元法在计算精度、效率以及稳定性等方面表现优异,适合高精度和长时程的远震波场模拟.虽然本文给出了二维远震弹性波场模拟, 但这一格式也可直接推广至三维情形.

发展适用于远震波场模拟的高效、高精度数值算法是远震伴随层析成像的基础.本文格式采用FK方法快速构造入射边界条件,利用强稳定的SSEM模拟区域模型中远震波场传播,可在保持高精度的前提下,通过采用大时间步长来节省计算时间.本文提出的SSEM-FK混合算法为高效、高精度的远震全波形层析成像提供了一种有效的正演模拟工具.

附录A 阶条件

根据泰勒展开得到

(A1)

根据本文格式(10)式,用第n步的UV表示的第n+1步的UV如下:

(A2)

将上式系数与泰勒展开对比,得到如下关系式:

(A3)

其中等式B1+B2+B3=1未列出因为其可以通过(A3)的其他等式推出.

附录B 辛条件

为了使格式(10)是辛的,格式(10)须满足辛条件(7)式.我们将(13)式阶条件的解代入格式(10),得到两个相邻时间步的UV的关系式为

(B1)

其中矩阵.矩阵A的元素为

(B2)

其中I表示单位矩阵.根据Hairer(2005), (7)式等价为

(B3)

其中J

(B4)

根据(B3)和(B4), 辛条件为

(B5)

aij代入(B5)等式成立,表明满足阶条件时辛条件自动成立.

References
Aki K, Lee W H K. 1976. Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes:1. A homogeneous initial model. Journal of Geophysical Research, 81(23): 4381-4399.
Bohlen T, Wittkamp F. 2016. Three-dimensional viscoelastic time-domain finite-difference seismic modelling using the staggered Adams-Bashforth time integrator. Geophysical Journal International, 204(3): 1781-1788. DOI:10.1093/gji/ggv546
Bostock M G, Rondenay S, Shragge J. 2001. Multiparameter two-dimensional inversion of scattered teleseismic body waves 1. Theory for oblique incidence. Journal of Geophysical Research:Solid Earth, 106(12): 30771-30782.
Cai W J, Zhang H, Wang Y S. 2016. Novel Symplectic Discrete Singular Convolution Method for Hamiltonian PDEs. Communications in Computational Physics, 19(5): 1375-1396. DOI:10.4208/cicp.scpde14.32s
Cai W J, Sun Y J, Wang Y S, et al. 2017a. Local discontinuous Galerkin methods based on the multisymplectic formulation for two kinds of Hamiltonian PDEs. International Journal of Computer Mathematics, 95(1): 114-143. DOI:10.1080/00207160.2017.1335866
Cai W J, Zhang H, Wang Y S. 2017b. Modelling damped acoustic waves by a dissipation-preserving conformal symplectic method. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Science, 473(2199): 20160798. DOI:10.1098/rspa.2016.0798
Chen J B. 2007. High-order time discretizations in seismic modeling. Geophysics, 72(5): SM115-SM122. DOI:10.1190/1.2750424
Chen J B. 2009a. Lax-Wendroff and Nyström methods for seismic modeling. Geophysics-Prospecting, 57(6): 931-941. DOI:10.1111/j.1365-2478.2009.00802.x
Chen L. 2009b. Lithospheric structure variations between the eastern and central North China Craton from S- and P-receiver function migration. Physics of the Earth and Planetary Interior, 173(3-4): 216-227. DOI:10.1016/j.pepi.2008.11.011
Chin S A, Anisimov P. 2006. Gradient symplectic algorithms for solving the radial Schrodinger equation. The Journal of Chemical Physics, 124(5): 054106. DOI:10.1063/1.2150831
Cohen G. 2002. Higher-Order Numerical Methods for Transient Wave Equation. Berlin, Heidelberg: Springer-Verlag.
De Basabe J D, Sen Mrinal K. 2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equation. Geophysics, 72(6): T81-T95. DOI:10.1190/1.2785046
Deblain M A. 1986. The application of high-order differencing to the shear wave equation. Geophysics, 51(1): 54-66. DOI:10.1190/1.1442040
Di Bartolo L, Dors C, Mansur W J. 2012. A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation. Geophysics, 77(5): T187-T199. DOI:10.1190/geo2011-0345.1
Feng K, Qin M Z. 2010. Symplectic Geometric Algorithms for Hamiltonian Systems. Berlin Heidelberg: Springer-Verlag.
Gao S, Zhang J F, Xu W L, et al. 2009. Delamination and destruction of the North China Craton. Chinese Science Bulletin, 54(14): 3367-3378.
Hashamder H, Ibrahim Z, Jameel M. 2011. Finite element analysis of nonlinear structures with Newmark method. International Journal of Physical Sciences, 6(6): 1395-1403.
Haskell N A. 1953. The dispersion of surface waves on multilayered media. Bulletin of the Seismological Society of America, 43(1): 17-34.
Hairer E, Lubich C, Wanner G. 2005. Geometric Numerical Integration I, second ed., Springer, Berlin, Heidelberg. New York, 2005. https://link.springer.com/book/10.1007%2F3-540-30666-8
Kind R, Yuan X H, Kumar P. 2012. Seismic receiver functions and the lithosphere-asthenosphere boundary. Tectonophysics, 536-537: 25-43. DOI:10.1016/j.tecto.2012.03.005
Komatitsch D, Vilotte J. 1998. The spectral element method:An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America, 88(2): 368-392.
Komatitsch D, Tromp J. 2002. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation. Geophysical Journal International, 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x
Langston C A. 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. Bulletin of the Seismological Society of America, 67(3): 713-724.
Li X F, Li Y Q, Zhang M G, et al. 2011. Scalar seismic-wave equation modeling by a multisymplectic discrete singular convolution differentiator method. Bulletin of the Seismological Society of America, 101(4): 1710-1718. DOI:10.1785/0120100266
Li X F, Wang W S, Lu M W, et al. 2012. Structure-preserving modelling of elastic waves:a symplectic discrete singular convolution differentiator method. Geophysical Journal International, 188(3): 1382-1392. DOI:10.1111/j.1365-246X.2011.05344.x
Liu H F, Dai N X, Niu F L, et al. 2014. An explicit time evolution method for acoustic wave propagation. Geophysics, 79(3): T117-T124. DOI:10.1190/geo2013-0073.1
Liu Q, Gu Y J. 2012. Seismic imaging:From classical to adjoint tomography. Tectonophysics, 566-567: 31-66. DOI:10.1016/j.tecto.2012.07.006
Liu S L, Li X F, Wang W S, et a. 2013. Optimal symplectic scheme and generalized discrete convolutional differentiator for seismic wave modeling. Chinese Journal of Geophysics (in Chinese), 56(7): 2452-2462. DOI:10.6038/cjg20130731
Liu S L, Li X F, Liu Y S, et al. 2013. Symplectic RKN schemes for seismic scalar wave simulations. Chinese Journal of Geophysics (in Chinese), 56(12): 4197-4205. DOI:10.6038/cjg20131222
Liu S L, Li X F, Wang W S, et al. 2014a. A new kind of optimal second-order symplectic scheme for seismic wave simulations. Science China Earth Sciences, 57(4): 751-758. DOI:10.1007/s11430-013-4805-0
Liu S L, Li X F, Wang W S, et al. 2014b. A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling. Journal of Geophysics and Engineering, 11(5): 055009. DOI:10.1088/1742-2132/11/5/055009
Liu S L, Li X F, Wang W S, et al. 2015. A modified symplectic scheme for seismic wave modeling. Journal of Applied Geophysics, 116: 110-120. DOI:10.1016/j.jappgeo.2015.03.007
Liu S L, Li X F, Wang W S, et al. 2015. A symplectic RKN scheme for solving elastic wave equations. Chinese Journal of Geophysics (in Chinese), 58(4): 1355-1366. DOI:10.6038/cjg20150422
Liu S L, Yang D H, Dong X P, et al. 2017a. Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling. Solid Earth, 8(5): 969-986. DOI:10.5194/se-8-969-2017
Liu S L, Yang D H, Ma J. 2017b. A modified symplectic PRK scheme for seismic wave modeling. Computers & Geosciences, 99: 28-36.
Liu S L, Yang D H, Lang C, et al. 2017c. Modified symplectic schemes with nearly-analytic discrete operators for acoustic wave simulations. Computer Physics Communications, 213: 52-63. DOI:10.1016/j.cpc.2016.12.002
Liu S L, Suardi I, Yang D H, et al. 2018. Teleseismic traveltime tomography of northern Sumatra. Geophysical Research Letters, 45(24): 13231-13239. DOI:10.1029/2018GL078610
Luo M Q, Liu H, Li Y M. 2001. Seismic wave modeling with implicit symplectic method based on spectral factorization on helix. Chinese Journal of Geophysics (in Chinese), 44(3): 379-388.
Luo M Q, Zhu G T, Liu H, et al. 2003. A hybrid method of matrix inversion suited for 3D implicit prestack depth migration. Chinese Journal of Geophysics (in Chinese), 46(5): 684-689.
McLachlan R, Atela P. 1992. The accuracy of symplectic integrators. Nonlinearity, 5(2): 541-562.
Monteiller V, Chevrot S, Komatitsch D, et al. 2013. A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model. Geophysical Journal International, 192(1): 230-247. DOI:10.1093/gji/ggs006
Monteiller V, Chevrot S, Komatitsch D, et al. 2015. Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method. Geophysical Journal International, 202(2): 811-827. DOI:10.1093/gji/ggv189
Nissen-Meyer T, Fournier A, Dahlen F A. 2008. A 2-D spectral-element method for computing spherical-earth seismograms-II. Waves in solid-fluid media. Geophysical Journal International, 174(3): 873-888.
Qin M Z, Zhang M Q. 1990. Multi-stage symplectic schemes of two kinds of Hamiltonian systems for wave equations. Computers & Mathematics with Applications, 19(10): 51-62.
Rawlinson N, Pozgay S, Fishwick S. 2010. Seismic tomography:a window into deep Earth. Physics of the Earth and Planetary Interiors, 178(3-4): 101-135. DOI:10.1016/j.pepi.2009.10.002
Rondenay S. 2009. Upper mantle imaging with array recordings of converted and scattered teleseismic waves. Surveys in Geophysics, 30(4-5): 377-405. DOI:10.1007/s10712-009-9071-5
Sun G. 1997. A class of explicitly symplectic schemes for wave equations. Mathematica Numerica Sinica, (1): 1-10.
Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation. Geophysical Journal International, 197(2): 1250-1267. DOI:10.1093/gji/ggu077
Tape C, Liu Q Y, Maggi A, et al. 2009. Adjoint tomography of the southern California crust. Science, 325(5943): 988-992. DOI:10.1126/science.1175298
Tian Y, Zhao D P, Sun R M, et al. 2009. Seismic imaging of the crust and upper mantle beneath the North China Craton. Physics of the Earth and Planetary Interiors, 172(3-4): 169-182. DOI:10.1016/j.pepi.2008.09.002
Tong P, Chen C W, Komatitsch D, et al. 2014a. High-resolution seismic array imaging based on an SEM-FK hybrid method. Geophysical Journal International, 197(1): 369-395. DOI:10.1093/gji/ggt508
Tong P, Zhao D, Yang D, et al. 2014b. Wave-equation-based travel-time seismic tomography-Part 1:Method. Solid Earth, 5(2): 1151-1168. DOI:10.5194/se-5-1151-2014
Virieux J. 1984. SH-wave propagation in heterogeneous media; velocity-stress finite-difference method. Geophysics, 49(11): 1933-1942. DOI:10.1190/1.1441605
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 M X, Xu S. 2015. Finite-difference time dispersion transforms for wave propagation. Geophysics, 80(6): WD19-WD25. DOI:10.1190/geo2015-0059.1
Wang W S, Li X F, Lu M W, et al. 2012. Structure-preserving modeling for seismic wavefields based upon a multisymplectic spectral element method. Chinese Journal of Geophysics (in Chinese), 55(10): 3427-3439. DOI:10.6038/j.issn.0001-5733.2012.10.026
Yang D H, Peng J M, Lu M, et al. 2006. Optimal nearly analytic discrete approximation to the scalar wave equation. Bulletin of the Seismological Society of America, 96(3): 1114-1130. DOI:10.1785/0120050080
Zhang H, Zhou Y Z, Wu Z L, et al. 2009. Finite element analysis of seismic wave propagation characteristics in Fuzhou basin. Chinese Journal of Geophysics, 52(3): 604-614. DOI:10.1002/cjg2.1382
Zhao L, Allen R M, Zheng T Y, et al. 2012. High-resolution body wave tomography models of the upper mantle beneath eastern China and the adjacent areas. Geochemistry, Geophysics, Geosystems, 13(6): Q06007. DOI:10.1029/2012GC004119
Zheng T Y, Chen L, Zhao L, et al. 2006. Crust-mantle structure difference across the gravity gradient zone in North China Craton:Seismic image of the thinned continental crust. Physics of the Earth and Planetary Interiors, 159(1-2): 43-58. DOI:10.1016/j.pepi.2006.05.004
刘少林, 李小凡, 汪文帅, 等. 2013a. 最优化辛格式广义离散奇异核褶积微分算子地震波场模拟. 地球物理学报, 56(7): 2452-2462. DOI:10.6038/cjg20130731
刘少林, 李小凡, 刘有山, 等. 2013b. 求解声波方程的辛RKN格式. 地球物理学报, 56(12): 4197-4205. DOI:10.6038/cjg20131222
刘少林, 李小凡, 汪文帅, 等. 2015. 求解弹性波方程的辛RKN格式. 地球物理学报, 58(4): 1355-1366. DOI:10.6038/cjg20150422
罗明秋, 刘洪, 李幼铭. 2001. 基于螺旋线上谱因式分解的地震波场隐式辛算法. 地球物理学报, 44(3): 379-388. DOI:10.3321/j.issn:0001-5733.2001.03.010
罗明秋, 朱国同, 刘洪, 等. 2003. 适用于三维隐式叠前深度偏移中矩阵求逆的混合算法. 地球物理学报, 46(5): 684-689. DOI:10.3321/j.issn:0001-5733.2003.05.016
孙耿. 1997. 波动方程的一类显式辛格式. 计算数学, (1): 1-10.
汪文帅, 李小凡, 鲁明文, 等. 2012. 基于多辛结构谱元法的保结构地震波场模拟. 地球物理学报, 55(10): 3427-3439. DOI:10.6038/j.issn.0001-5733.2012.10.026