李冰非， 李小凡， 李峰， 龚飞

1 中国地震灾害防御中心， 北京 100029 2 中国科学院地质与地球物理研究所， 北京 100029 3 中国地质大学地球物理与空间信息学院， 武汉 430074

0 引言

地球自由振荡是由大地震在有限尺度的地球内部激发的三维驻波，即整体地球的固有振荡(Benioff et al., 1961).地球自由振荡的固有频率和频谱形态与地球内部结构密切相关(Lapwood and Usami, 1981)，包括地球内部横向非均匀性及径向非均匀性(Tsuboi et al., 1985; Tsuboi,1995)、地球内部各向异性(Woodhouse et al., 1986; Tromp, 1993,1995)、介质非完全弹性(Masters and Widmer, 1995)、以及介质对波的衰减(Davies and Bunge, 2001; Alterman et al., 1974; Pollitz et al., 1987; Park et al., 2005)等性质.由于应用驻波方法研究地球内部结构具有对采样点要求不高的优势，地球自由振荡数值模拟研究可以方便、经济地重建全球多尺度三维地球内部结构及获得地球内部某些物理场的三维特征参数分布.

谱元法是近年来大尺度地震波及地球自由振荡数值模拟的重要方法之一(Komatitsch et al., 2002, 2004; Komatitsch and Tromp, 2002; Tsuboi et al., 2004; Yan et al., 2014).谱元法将有限元算法和谱方法相结合，兼顾了有限元方法的灵活性和谱方法的高精度(Komatitsch and Vilotte, 1998; Liu et al., 2017a).谱元法是弹性波方程在空间近似的高阶变分方法，Komatitsch和Tromp(1999)在谱元单元上采用Gauss-Lobatto-Legendre(GLL)积分并结合高阶精度的拉格朗日插值基函数，得到具有对角化的质量矩阵的显式时间域微分方程组，且数值解更好地趋近于理论解，模拟精度和计算效率得以提高；而后又将谱元方法应用于全球地震波传播模拟研究(2002).Tromp等(2010)在多年研究的基础上，将数值方法和高性能并行计算相结合，研制谱元法程序包SPECFEM3D_GLOBE，可以用来精细剖分地球模型、模拟三维地球模型或区域模型中地震波传播.赵明等(2013)拓展了一种谱元-简正振型耦合方法适用于各向异性介质中地震波传播数值模拟.近二十年来，研究者们将超大规模高性能并行计算数值模拟方法应用于全球地震波传播研究，将计算结果与实际观测的超长波记录对比研究，获得了一系列较为令人满意的结果(刘少林等, 2021).特别应该提到的是，严珍珍等(2008a)应用谱元法结合弹性波理论及高性能并行计算对由特大地震激发的地球自由振荡进行了成功的大规模数值模拟，从而证实了这类方法的可行性，为由大地震激发地球自由振荡的数值模拟研究奠定了基础.

尽管理论上时间域数值方法均可处理诸如地球内部横向非均匀性及地震激发地球自由振荡之全过程等问题，但现有地球自由振荡时间域数值模拟技术均基于非保结构计算，其长时程计算中的高积累误差难以避免，很难精确处理地球自由振荡模拟这一长时程追踪问题.近年来已有声波及弹性波保结构模拟方法问世(罗明秋等, 2001; Li et al., 2012; 汪文帅等, 2012; 陈世仲等, 2016; Liu et al., 2017b,c; 李冰非等, 2019)；Li 等(2012)及Gao等(2016)进一步将逆时离散变换与三阶辛相结合，获得了长时程的低频散数值模拟，这类保结构方法非常适合处理有关保守动力学体系的问题.若将地球近似为弹性球体，用上述保结构方法模拟其自由振荡问题将是很好的选择；然而，地球本身是一具有衰减性质的非完全弹性球体，真实地球介质极其复杂, 除了各向异性以外,可能还有多相多孔性以及黏性(Kuznetsov, 1997; Kustowski et al., 2008a)，其自由振荡(驻波)问题涉及非保守动力学系统问题.针对这一问题，Li等(2015)发展了阻尼声波保结构辛算法相关理论，数理逻辑的合理性和方法的可行性均得到验证.Cai等(2017)基于同样的模型方程提出了另一种耗散保结构模拟方法，对于带耗散的地震波运动波动方程，可以分为一个守恒系统和一个耗散系统来求解.为更精确地模拟衰减模式下的地球自由振荡，需要发展新的、高效驻波保结构数值方法，研究适用于衰减介质的保弥散数值模拟方法.

本文将保结构计算延拓至非保守动力学体系，发展了一种更切合实际的、适用于处理地球自由振荡衰减问题的高精度、高效率模拟计算方法，以2012年印度尼西亚苏门答腊西北海域8.6级大地震观测记录为数据源，进行了衰减模式下的全球自由振荡数值实验，研究地球深部(如地核)弥散衰减对地球自由振荡的影响和各种衰减模式对地球自由振荡的影响.

1 基于辛-谱元法的地球自由振荡保弥散衰减数值模拟

1.1 地球自由振荡波动方程

地球的一般弹性运动方程为

(1)

其中ρ为密度，u为质点位移，ω为地球角转动矢量，f表示外力和重力的扰动，λ，μ为Lame常数(严珍珍等，2008b).

为便于表达，我们对地球的一般弹性运动方程进行简化，忽略大气变化因素，忽略自转、重力、地球椭率以及横向不均匀性，考虑简化的地球弹性运动方程为：

(2)

其中F表示震源.数值模拟地球自由振荡，本质上也是求满足适当初边值条件的地球振动的常系数偏微分方程组的数值解.

1.2 基于辛-谱元法的地球自由振荡保弥散衰减数值模拟

用谱元法将方程组(2)离散化后可写为：

(3)

其中M为全局质量矩阵，K为全局刚度矩阵，U表示全局位移向量.

对于本文研究的问题，弥散衰减介质中弹性波方程可写为：

(4)

其中u为位移矢量，σ为应力张量，ε为应变张量，c为四阶弹性张量，D为弥散阻尼系数，f为外力.谱元法离散空间域后的微分方程组可写为：

(5)

将式(5)应用阻尼项保结构辛格式(Li et al., 2015)进行时间离散，可得如下基于谱元法的保弥散衰减波动方程的时间格式：

(6)

其中V为速度向量. Δt为时间步长，c2，d1，d2为辛系数，U(n)和V(n)的上标表示第n时间步，U1和V1为中间变量.

本文将采用基于(6)式的辛格式-谱元法开展衰减模式下的地球自由振荡数值模拟研究，并探究地核中不同衰减模式及其对地球自由振荡的影响.

2 数值实验

本节通过数值模拟实验验证本文提出的方法的有效性.采用2012年4月11日印度尼西亚苏门答腊西北海域8.6级地震(2.327°N，93.063°E)作为震源，图1所示为本次地震事件震源位置和Global Centroid-Moment-Tensor (CMT) Project提供的8.6级地震的震源球以及全球地震台网(Global Seismographic Network, GSN)公布的152个台站所在位置(https:∥www.iris.edu/hq/programs/gsn).数值模拟使用Computational Infrastructure for Geodynamics(CIG)提供的SPECFEM3D-GLOBE谱元法程序包(经本文保弥散-辛格式升级).

图1 地震事件震源和全球152个台站位置示意图(https:∥www.iris.edu/hq/programs/gsn)Fig.1 Epicenter, 152 stations on globe and earthquake mechanism

2.1 长时程数值模拟实验

本节分别采用保弥散-辛格式和SPECFEM3D-GLOBE谱元法程序包中所采用的Newmark算法模拟同一模型，考虑地形、椭率、重力及自转等因素，模拟地震波传播600 min.采用S362ANI地球模型(Kustowski et al., 2008a,b)，并将数值模拟结果进行同参数处理.图2所示为NNA台站Z分量地震记录波形对比，从图中波形放大图可以看出，模拟时长约388 min时Newmark算法的结果开始显现出高频数值频散(波形曲线出现高频毛刺)，本文算法的波形记录则比Newmark算法在长时间计算后更加稳定，显示出在模拟长时程地震波传播时的优越性.采用同样的时间步长(dt=0.1425 s)，在同样的计算节点条件下，本文保弥散-辛格式模拟耗时68 h 23 m 24 s，Newmark算法模拟耗时71 h 05 m 10 s，本文保弥散-辛格式模拟长时程地震波传播时更加稳定，同时不增加计算量甚至计算效率略有提高.本实验仅简要展示低阶辛格式处理大尺度地震波传播问题的有效性及抑制数值频散的能力，同类辛格式的数值频散前人已有详尽的分析(Li et al., 2015).

2.2 内外核弥散阻尼相等时地球模型自由振荡数值实验

本小节设计两组实验，其整体均采用S362ANI地球模型，并考虑地形、椭率、重力及自转等因素.第一组实验采用本文保弥散衰减方法模拟外核和内核弥散衰减模型(设内核和外核弥散阻尼D=-0.1 s-1),和采用SPECFEM3D-GLOBE谱元法程序包自带的Newmark算法模拟无弥散衰减的模型，模拟2012年印度尼西亚北苏门答腊8.6级地震发生后地震波传播600 min (dt=0.1425 s，模拟252700步).将数值模拟结果与实际地震观测数据进行对比分析，本文采用的实际地震观测数据来源为IRIS(Incorporated Research Institutions for Seismology).根据地球自由振荡的甚低频特性，我们将两种数值模拟算法的结果和实际地震观测结果(IRIS下载)采用0.0039 Hz的低通滤波，选择Chirp-z变换算法求取频谱进行对比，并引入地球自由振荡理论值(简振模)作为参考.由于外核中不存在环型振荡，本文仅研究地核弥散阻尼对球型振荡的影响.且随着简正模球谐阶l的增大，自由振荡所发生的深度越来越浅(Jeans, 1923)，故针对地核的研究我们选择频率较低的低阶简正振型来分析.

图3所示为本文数值模拟算法的结果与无衰减(弥散衰减或黏滞衰减)的数值模拟结果和实际地震观测结果频谱图，对比所用的简正模采用的是由一维地球模型得到的自由振荡简正模(图中虚线)，可作为近似参考或约束.数值模拟结果频谱图显示，相较于无衰减模型，本文实验在地核加入了弥散阻尼后，低阶的简正振型0S2、0S3、…、0S9和3S1振型的振幅衰减明显，模拟结果更加接近实测数据.而0S10之后的振型未见振幅的衰减.这一结果与方明(1991)进行的对球型振荡的深度与球谐阶之间的关系的数值实算结果一致.下面，我们在模拟研究地核中不同的弥散衰减对地球自由振荡的影响时，选择0S2、0S3、…、0S9和3S1振型的频谱来进行对比分析.

图3 引入地核弥散阻尼模型和无衰减模型以及实测数据(IRIS)频谱图对比Fig.3 Spectrum comparison of models with and without the core dispersive attenuation and IRIS data

第二组实验将内核和外核弥散阻尼均设为D=-0.25 s-1，其他参数设置与上一组相同.图4所示为不同弥散阻尼模型数值模拟实验、无衰减(弥散衰减或黏滞衰减)数值模拟结果和实际地震观测结果频谱图对比，可看出地核弥散阻尼加大后，不影响地球自由振荡的频率，即本征频率对应性良好，只影响振幅衰减.由于地球衰减和仪器响应等因素，实际地震台阵记录到的长时间传播的地震波信号能量会大大衰减，超长周期信号较难识别，故实测数据的甚低频段能量非常低，可以识别出的自由振荡振型仅为0S5、0S6、0S7、0S8、0S9振型，在这些振型处本文模拟的地核自由振荡振幅更加接近实测数据.由于计算所用服务器和计算稳定性的约束，数值模拟无法达到台网地震仪的采样频率.实测数据来源为IRIS，地震仪器的采样率为50 sps (sample per second)，采样间隔为0.05 s.本文方法采用的时间步长为0.1425 s.因此数值模拟结果的有效信息远少于实测数据的有效信息，信息量约相当于实测数据的三分之一.但即使是在信息量相对较少的情况下，数值模拟的结果依然是可靠的，说明了本文算法具有高效求解全球尺度地震波长时程传播问题及地球自由振荡问题的优势.

图4 地核弥散阻尼D=-0.25 s-1模型(蓝)、D=-0.1 s-1模型(黑)和无衰减模型以及实测数据频谱图对比Fig.4 Spectrum comparison of models with core dispersive attenuation D=-0.25 s-1 (blue)，D=-0.1 s-1 (black) and model without the core dispersive attenuation and IRIS data

在此基础上我们试验了仅增大内核弥散阻尼，实验结果表明，当内核弥散阻尼系数增高至大于外核弥散阻尼系数时，部分振型频率和振幅偏离实际观测数据.

2.3 仅增大外核弥散阻尼时自由振荡数值实验

本小节进行仅增大外核弥散阻尼时自由振荡的数值实验.给出两组模型实验结果，分别为将外核弥散阻尼设为DOC=-0.9 s-1、内核弥散阻尼设为DIC=-0.3 s-1，以及将外核弥散阻尼设为DOC=-1.1 s-1、内核弥散阻尼设为DIC=-0.1 s-1.由图5可见，外核弥散阻尼较大时模型模拟结果最接近实际地震观测数据， 0S5、3S1、0S6、0S7、0S8、0S9振型的模拟结果与观测数据在本征频率对应性、幅值关系等方面拟合较好，说明外核弥散衰减对地球自由振荡的影响大于内核弥散衰减.图6为地核敏感振型3S1、3S2频谱分析结果，从图中可以看出，外核弥散阻尼DOC=-1.1 s-1时3S1、3S2振型频谱幅值衰减，更加接近实测数据.表1和表2分别给出外核弥散阻尼DOC=-0.9 s-1模型和DOC=-1.1 s-1模型的数值模拟结果和实测数据的频谱分析结果，两个数值模型识别出的球型振荡0S3到0S10振型本征频率值相同，与实测频率的误差如表1所示；表2可以看出外核弥散阻尼DOC=-1.1 s-1模型比DOC=-0.9 s-1模型幅值均减小，更加接近实测数据幅值.在频谱的甚低频段，频率的可识别度对地球自由振荡的模拟研究至关重要.因此，本文采用的地核弥散阻尼等效模型可作为一种模拟甚低频地球自由振荡的有效模型.

图5 不同地核弥散阻尼DIC=-0.3 s-1& DOC=-0.9 s-1模型(蓝)、DIC=-0.1 s-1& DOC=-1.1 s-1模型(红)和无衰减模型以及实测数据频谱图对比Fig.5 Spectrum comparison of models with core dispersive attenuation DIC=-0.3 s-1& DOC=-0.9 s-1 (blue)，DIC=-0.1 s-1& DOC=-1.1 s-1 (red) and model without the core dispersive attenuation and IRIS data

图6 不同地核弥散阻尼DIC=-0.3 s-1& DOC=-0.9 s-1模型(红)、DIC=-0.1 s-1& DOC=-1.1 s-1模型(蓝)模拟结果地核敏感振型(3S1、3S2)频谱图对比Fig.6 Core sensitive normal modes (3S1，3S2) spectrum comparison of models with core dispersive attenuation DIC=-0.3 s-1& DOC=-0.9 s-1 (red) and DIC=-0.1 s-1& DOC=-1.1 s-1 (blue)

表1 数值模拟结果与实测数据频谱分析——频率Table 1 Spectrum analysis of simulation and observation data—frequency

表2 数值模拟结果与实测数据频谱分析——幅值Table 2 Spectrum analysis of simulation and observation data—amplitude

继续增大外核弥散阻尼设为DOC=-1.5 s-1，内核弥散阻尼设为DIC=-0.1 s-1, 模拟结果如图7所示，部分振型的频率难以分辨，说明此模型外核弥散阻尼过大，与实际地球内部结构不相符.

图7 不同地核弥散阻尼DIC=-0.1 s-1& DOC=-1.5 s-1模型(蓝)、DIC=-0.3 s-1& DOC=-0.9 s-1模型(红)和无衰减模型以及实测数据频谱图对比Fig.7 Spectrum comparison of models with core dispersive attenuation DIC=-0.1 s-1& DOC=-1.5 s-1 (blue)，DIC=-0.3 s-1& DOC=-0.9 s-1 (red) and model without the core dispersive attenuation and IRIS data

综上，采用内核弥散阻尼为0.1 s-1、外核弥散阻尼为1.1 s-1的地核衰减模型时，地球自由振荡数值模拟结果最接近实际观测记录.以上数值模拟实验结果表明，对甚低频段的地球自由振荡而言，弥散衰减为地核自由振荡衰减的重要成分.对于低频段，主波长远远大于等效非均匀尺度, 引起较强的Rayleigh散射，使得等效散射衰减(弥散衰减)较黏滞衰减更为明显；对于高频段，由于主波长远远小于等效非均匀尺度，其散射主要为前向散射(沿传播方向)，其等效散射衰减则远弱于高频情况下的黏滞衰减.前期数值实验结果表明，当λ/a≫1(Rayleigh散射模态，其中λ为波长，a为等效平均自由程)或λ/a～1(Mie散射模态)时，弥散衰减为衰减之主要成分；反之黏滞衰减占主导地位(李冰非，2019).就地球自由振荡与地球内部结构而言，大多为λ/a≫1情形，因而弥散衰减为衰减的主要成分.

数值实验结果还显示对于地球自由振荡超低频段外核弥散衰减的影响大于内核弥散衰减，这种现象可能是由外核和内核非均匀尺度的差异以及外核和内核体积的差异造成.

3 总结

本文基于地震台网实际观测记录，得到了更加适用于真实地球自由振荡数值模拟的地核等效弥散衰减系数，但鉴于计算资源的限制，本文选择的模型单元剖分不够精细，未来预计进一步增加剖分精细程度；本文通过数值实验模拟给出地核等效弥散阻尼，但限于模型单元剖分精细程度不够及未考虑介质各向异性等因素，其精度评估亦有难度.未来计划对本文等效地核弥散衰减开展更深入的研究.

对于存在衰减效应的地球自由振荡问题，本文发展了具有保弥散衰减性质、适用于低频问题的数值模拟方法，将其应用于全球尺度地球自由振荡衰减效应数值模拟，发展了基于辛格式-谱元法的保弥散数值模拟方法，提供了一种更切合实际的、适用于处理地球自由振荡衰减问题的高精度、高效率模拟计算方法.数值模拟实验结果表明，在地核区域，对甚低频段的地球自由振荡而言，外核中的自由振荡衰减主要源自弥散衰减.本文的研究方法提供了一种可通过自由振荡这样的甚低频振荡研究探索地核非均匀性的有效途径.

致谢非常感谢中国科学院地质与地球物理研究所计算模拟实验室为本文研究提供高性能计算集群，非常感谢IRIS提供地震数据以及CIG提供SPECFEM3D-GLOBE程序．