地球物理学报  2015, Vol. 58 Issue (5): 1675-1691   PDF    
上月壳中的散射引起月震尾波的数值模拟研究
姜祥华1, 王彦宾1, 古村孝志2    
1. 北京大学地球与空间科学学院地球物理系, 北京 100871;
2. 日本东京大学地震研究所, 东京 113-0032
摘要:在阿波罗月震记录中普遍存在着强烈持久的尾波信号,这样的波形特征无法用均匀分层月球模型解释.一个普遍被接受的解释是月震尾波由月球浅层结构对月震波的散射引起.我们采用基于交错网格的伪谱和有限差分混合方法模拟研究非均匀上月壳对月震波的散射效应,在此基础上解释月震尾波的形成机制,并估计出上月壳速度扰动的强度.我们发现,在均匀分层模型基础上,进一步考虑上月壳中的非均匀结构对月震波的散射效应,能有效地解释月震信号中强烈持久的尾波.我们认为月震尾波可能是由上月壳中的低波速、低衰减和散射这三个因素的共同作用所引起.采用不同的扰动标准差模拟上月壳的非均匀性,并比较模拟波形与真实月震图的相似程度,我们发现上月壳中速度扰动的标准差应该在3%到5%之间,很可能接近于3%.
关键词散射     月震尾波     数值模拟     计算地震学     地震波传播    
Numerical simulation of lunar seismic coda caused by scattering in upper Moon crust
JIANG Xiang-Hua1, WANG Yan-Bin1, Takashi Furumura2    
1. Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871, China;
2. Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan
Abstract: The Apollo lunar seismic data generally exhibits a long and strong coda, which masks seismic phases other than initial P and S waves and causes difficulties in analyses of secondary waves. Recent numerical simulation of 2D global seismic wave propagation in a whole-Moon model showed that lunar seismic coda couldn't be produced by a layered Moon model. A widely accepted explanation is that lunar seismic coda is resulted from scattering in a near-surface layer, which could be several to more than twenty kilometers thick. This study presents numerical simulations of scattering effects caused by small-scale heterogeneity in upper Moon crust in order to explain the coda.
Wave equations in velocity-stress form are solved in a 2D cylindrical coordinate system. The spatial derivatives are calculated by a hybrid pseudospectral and finite difference method defined on staggered grid. The background layered model is based on recent published seismic velocity and attenuation models of the Moon. In order to simulate the heterogeneity of upper Moon crust, von-Karman auto correlation function is used to generate velocity fluctuations, superimposed on the background velocity.
We firstly check the effectiveness of scattering in upper Moon crust to reproduce lunar seismic coda. Wavefields are calculated for an upper Moon crust model without scattering and one with scattering. Calculations are performed for a 71 km shallow moonquake and an 867 km deep moonquake. Comparison of the wavefield snapshots generated by the two models shows that wave energy can last a very long time in the upper Moon crust due to scattering in this layer. Comparison of waveforms shows that scattering in upper Moon crust can effectively produce long-duration coda. Comparison with Apollo seismograms shows that the duration and strength of coda caused by scattering in upper Moon crust can be well consistent with observations. The mechanism which produces lunar seismic coda is tentatively discussed. We propose that scattering, low attenuation and low velocity in the upper Moon crust may all contribute to the lunar seismic coda.
Effects of fluctuant strength on synthetic coda are discussed. Different standard deviations of velocity fluctuations are set for upper Moon crust: 1%, 3%, 5%, 7% and 9%. For both shallow and deep moonquakes, wavefields are simulated for each model and the characteristic decay times of synthetic coda are calculated at four epicentral distances: 22°, 45°, 67.5° and 90°. The characteristic decay times firstly increase with standard deviation of velocity fluctuations and then decrease. The characteristic decay times reach the largest value where standard deviation equals to 3%. Finally, similarities between synthetics and Apollo lunar seismograms are measured via correlation coefficient in order to estimate the most realistic strength of velocity fluctuations. We propose that the standard deviation of velocity fluctuations should be between 3% and 5%, and it's likely close to 3%.
Key words: Scattering     Lunar seismic coda     Numerical simulation     Computational seismology     Wave propagation    
1 引言

月球的内部结构和物质组成对于研究月球和地-月系统的起源与演化历史具有十分重要的科学意义,同时对探月飞行器的轨道控制和探测器的着陆有影响(欧阳自远,2005占伟和李斐,2007).地震学方法是研究地球内部结构最有效方法之一,它对月球内部结构的研究也有着非常重要的作用(姜明明和艾印双,2010).

最早的月震学研究始于20世纪的美国阿波罗登月计划(1969—1977).该计划中所布设的四台月震仪共记录了超过12500个月震事件(Nakamura,2005).这些月震资料为到目前为止的月震学研究提供了数据上的支持.和典型的地震观测波形相比,月震波形表现出非常大的差异.月震信号中普遍存在着很强的尾波特征,并且尾波的持续时间很长,振幅衰减很缓慢(Dainty et al., 1974).强烈的尾波信号掩盖了各种后续震相,使得识别P波和S波以外的次达震相非常困难.

在地震学研究中,通常认为尾波是由介质小尺度(相对波长而言)的非均匀结构对波的散射所引起(Sato and Fehler, 1998).对于月震波强烈而持久的尾波特征,早期研究者指出可能是由月球介质尤其是浅层介质对月震波的强烈散射所引起(Latham et al., 1970Nakamura,1977).研究月球浅层介质对月震波的散射效应非常重要.它不但有助于弄清月震尾波信号的成因,也能提供关于月球浅层速度结构更精细的约束(在更小尺度上的速度变化),还有助于理解月壳的形成及演化历史.

最近,Blanchette-Guertin等(2012a)基于已有的月震资料对月震波信号的尾波衰减特征进行了系统研究,讨论了尾波衰减时间随事件类型、频率以及震中距的变化.他们指出散射是在整个月球的近表浅层中都存在的现象,因为所有类型和深度的事件在不同震中距的记录都展现出强烈的尾波特征.Latham等(1972)指出散射层的可能厚度为几公里至20 km.Dainty等(1974)将扩散理论应用于月震波的散射模拟,通过拟合登月舱冲击信号记录的能 量包络函数,推断出散射层的厚度范围为14~25 km.

前人的研究为我们探寻月震信号尾波的成因提供了方向,但也留下了不少有待探明的问题.比如,散射层中速度结构非均匀的剧烈程度为多大?以及尾波衰减如何随速度非均匀的剧烈程度变化?这些问题正是本文的研究动机之所在.

作为对观测资料的补充,数值模拟技术在研究地震波传播的物理过程以及全球地震学中起着非常重要的作用(Tromp,2007).对月震波而言,由于强烈尾波掩盖了各种后续震相,因此通过数值方法模拟月震波的传播过程,对于理解和解释月震波形特征就显得更为重要.纵观国内外,目前已有少量关于月震波的数值模拟研究.Lawrence和Johnson(2010)将地震声子法(Shearer and Earle, 2004)应 用于全月球模型,模拟了表面粗风化层(megaregolith)及月壳非均匀性对月震波散射的效应,发现三维非均匀性对月震波形的散射特征有重要贡献.Blanchette-Guertin等(2012b2013)改进了地震声子法,在全月球模型中既考虑近月表的强烈散射层,也包含全月球的背景散射,通过全波形模拟以及和观测资料的对比,讨论了散射层厚度、不均匀尺度以及固有衰减对月震波形的影响.Wang等(2013)采用伪谱和有限差分混合方法模拟了全月球模型中月震波的传播过程,从理论上预测了各种震相以及它们的到时,并指出近月表低速层引起的多重反射对月震波波形有很大影响.现有的关于月震波散射的模拟全都采用声子法,该方法是一种基于射线理论的间接方法,优点是计算效率高,但不如直接求解波动方程的模拟方法对波场信息的包含完整和准确.而Wang等的月震波模拟,虽然基于波动方程的直接求解,但未能将散射层的影响包含到模拟中,所以不能有效解释尾波特征.因此,还需要一个将散射效应直接包含到波动方程中进行求解的模拟工作.本研究是对这一空白的很好填补.另外,Lawrence,Johnson,Blanchette-Guertin等采用地震声子法的月震波散射模拟,仍没有给出关于散射层中速度结构非均匀的剧烈程度的估计.

本文采用之前Wang等(2013)所用的伪谱和有限差分混合方法模拟月球浅层非均匀结构对月震波的散射,以解释月震记录的尾波特征.结合前人对散射层厚度的估计,在模拟中选取上月壳(15 km)作为散射层,即只考虑上月壳速度结构的非均匀性对月震波的影响.在本方法能模拟的频率下,将检验上月壳散射对月震尾波形成的有效性.还将探讨尾波衰减随上月壳速度扰动幅度的变化.最后,试图通过模拟结果与实际记录的比较,给出对上月壳速度扰动标准差的估计. 2 算法与模型 2.1 基于交错网格的伪谱和有限差分混合方法

考虑P波和SV波在柱坐标系下的传播,并设所有变量不随z变化,可导出二维柱坐标系下的P-SV波动方程.方程选取速度应力形式:

其中,vr和vθ分别是径向和横向的质点运动速度,ρ是介质密度,fr和fθ表示径向和横向上的体力,σrrσσθθ表示应力分量,λ和μ是拉梅常数.

为了求解上述方程,采用如图 1所示的交错网格在r和θ方向对方程组进行离散化.应力和速度的分量位于节点上不同位置,相互间错开半个节点间距.在θ方向,采用相同的节点间距Δθ.在r方向上,采用相同的节点间距Δr.

图 1 二维柱坐标系下交错网格的节点分布(Wang et al., 2013)Fig. 1 Distribution of grid points for the staggered grid scheme in 2D cylindrical coordinates(Wang et al., 2013)

采用基于交错网格的伪谱和有限差分混合方法来求解方程组(1)和(2).在横向,对θ的偏导数采用傅里叶伪谱法计算(Witte and Richards, 1990)

径向上,对r的偏导数采用交错网格的四阶有限差分格式(Levander, 1988)

对时间的偏导数采用与Furumura等(1998)的文章中一样的有限差分格式.

求解区域为通过月球大圆选取的一个深1000 km 在水平方向延伸120°的扇形剖面,如图 2所示.这样的求解区域,有四个边界条件需要处理.在月球的表 面施加自由边界条件,即令表面处的牵引力为0.在 扇形区域的左右两边以及底边,采用Cerjan等(1985)给出的由20个节点组成的吸收边界条件,以减小人为带来的反射.方程组(1)中的体力项对应二维柱坐标系下的线源,通过下面的公式由震源矩张量分量给出(Wang et al., 2001)

其中r0和θ0为震源中心坐标,Mrr,M和Mθθ为震源矩张量分量.式中的δ函数采用Herrmann伪δ函数近似(Herrmann,1979),这是一个保持单位面积不变的钟形函数.介质的滞弹性衰减通过在每一时间步的计算中对每个网格节点乘以由Graves(1996)所给出的衰减系数实现

其中,A(r,θ)表示滞弹性衰减系数,f0为参考频率,Δt为时间步长,Q(r,θ)为品质因数.
图 2 模拟中的区域模型Fig. 2 Regional model in our simulation
2.2 月球分层结构模型

Weber等(2011)最近利用台阵技术重新分析了阿波罗月震波形资料,首次发现了月核产生的反射波和转换波,并据此提出了月核和核幔过渡带的速度结构模型.我们采用他们给出的分层速度和密度结构模型作为模拟的背景模型(不包含表面沉积 层),如图 3a所示(为了模拟散射,还需要在此模型基础上叠加随机速度扰动,参见2.3节).Garcia等(2011)对阿波罗月震波形资料进行了更细致的再分析处理,分辨出了从月核反射的SH波信号,通过联合反演月震与月球重力测量数据,给出了一个非常初步的参考月球模型(VPREMOON).参考他们的 模型,我们得到了模拟中所用的反映滞弹性衰减的Q值模型,如图 3b.在我们的分层模型中,上月壳厚度为15 km,下月壳 厚度为25 km,从40 km到1000 km深为月幔.上月 壳的P波速度为3.2 km·s-1,S波速 度为1.8 km·s-1;下月壳的P波速度为5.5 km·s-1,S波速度为3.2 km·s-1;月幔中P波速度的变化范围在7.6~8.5 km·s-1之间,S波速度的变化范围在4.4~4.5 km·s-1之间.由于月球的结构、物质组成、热的历史和演化历史与地球不同,所以月球中的Q值和地球中的Q值有很大差异(Khan et al., 2006).例如,本模型中月壳中的Q值高至几千,而地壳中的Q值大小一般为几百,月幔中的Q值也整体比地 幔中的Q值(Dziewonski and Anderson, 1981)要高.

图 3(a)波速和密度的分层模型,(b)分层Q值模型Fig. 3(a)Velocity and density structure in a layered model,(b)Q values in a layered model
2.3 随机速度扰动模型与自相关函数

为了模拟上月壳中的非均匀结构对月震波的散射,还需要在2.2节中给出的分层模型基础上叠加上一个随机速度扰动.月壳由于自形成之初持续遭受各种天体撞击,造成整个表面浅层分布着许多破碎体(McGetchin et al., 1973Thompson et al., 2009).这样的破碎体分布会在结构上引起近似于随机的速度扰动,正是这样的随机速度扰动造成了月震波的散射.在地球介质对地震波散射的研究中,人们普遍采用统计学的方法去描述介质的随机速度扰 动,并获得了巨大成功(例如Kennett and Furumura, 2013). 对于速度随位置变化的非均匀介质,可以将速度拆分成平均背景值项v0和扰动项δv之和:

这里 y 是位置矢量,ξ(y)表示速度的相对扰动大小.用〈…〉表示整体平均值,有

对ξ(y)做自相关,得自相关函数(ACF):

这里,x表示空间延迟矢量.当介质各向同性时,自相关函数仅是r=x的函数.自相关函数是对介质的非均匀空间尺度和扰动幅度大小的统计学测量.相对扰动的幅度通过相对扰动的方差来衡量:

ε2表示扰动的方差,ε表示扰动的标准差.扰动的空间变化特征由关联长度a决定.自相关函数的傅里叶变换叫做功率谱密度函数(PSDF):

k是波数矢量(Sato and Fehler, 1998).

常用的自相关函数有三种,分别是高斯自相关函数(Gaussian),指数自相关函数(Exponential),以及冯 · 卡尔曼自相关函数(von Karman).表 1给出了二维情形的三种自相关函数以及它们相应的功率谱密度函数.

表 1 二维情形下的三种自相关函数以及功率谱密度函数Table 1 Three types of correlation functions and power spectral density functions in a two dimensional case

通过自相关函数生成随机速度扰动的步骤如下:(1)产生一个白噪声,进行傅里叶变换得到波数域的白噪声谱;(2)用自相关函数的功率谱密度函数的平方根在波数域对上述白噪声谱滤波;(3)进行傅里叶反变换得随机速度扰动(Frankel and Clayton, 1986).在我们的模拟中,采用冯 · 卡尔曼自相关函数来产生上月壳的随机速度扰动.图 4是用冯 · 卡尔曼自相关函数产生的随机速度扰动的一个示例.

图 4 由冯·卡尔曼自相关函数生成的随机速度扰动Fig. 4 R and om velocity fluctuations produced by von Karman autocorrelation function
3 上月壳散射产生月震尾波的有效性

本节将检验上月壳的散射是否足以产生像真实月震记录那样强烈又持久的尾波.我们将分别模拟均匀分层模型中的波场和上月壳含有散射时的波场,并对比两种情形的波场快照以及波形记录,以弄清上月壳散射对月震波传播的影响.最后,将对比模拟波形和真实记录的波形,以探明上月壳散射是否能在记录中造成足够强的尾波.

月球内部发生的月震按震源深度的分布不同分为浅震和深震.浅源月震的深度通常在50~220 km以内,震级在4.8到5.5级之间,最大地震矩为 3×1014Nm.深源月震具有丛集性,深度通常在700~1200 km之间,震级据估计最大不超过3级,最大地震矩为5×1013Nm(Goins et al., 1981).浅源月震很可能具有构造起源,与地球中的板内地震比较类似(Nakamura et al., 1979).深源月震据推测可能是由地球施加给月球的潮汐引力所引起,震源机制尚不能完全确定.不过,Nakamura(1978)Koyama和Nakamura(1980)对于深震A1集群记录的P波和S波振幅以及极性的分析,表明深震的机制可能比较接近于水平断层.在我们的模拟中,浅震的深度为71 km,深震的深度为867 km,分别对应后面用于比较的两个实际月震记录的震源深度(参见表 3).源都放置在扇形上20°的位置,如图 2中所示.对于震源机制,采用二维的双力偶线源.浅震取Mrr=-1.0,Mθθ=1.0,M=Mθr=0,相应于45°倾角的走滑断层;深震取Mrr=0,Mθθ=0,M=Mθr=1.0,相应于水平滑动断层.

结合计算机性能和计算时间的考虑,对于图 2所示的扇形区域,将θ方向离散为2048个节点,r 方向离散为1000个节点.月球的半径为1737.1 km,这样横向的节点间距将从自由表面的1.78 km变到底部边界的0.76 km.径向节点间距保持1 km不变.对于伪谱法,这样的横向网格划分能计算的最小波长为横向最大节点间距的2.5倍(Wang et al., 2013).λmin=4.45 km是模拟中能计算的最小波长,结合模型的最小波速,我们能计算的最大频率为fmax=Vminmin=0.4Hz.对于震源时间函数,采用 宽度为2.5 s的Herrmann伪δ函数.为确保计算稳定,时间步长 Δt应满足稳定性条件Δt≤α×Δmin/Vmax(Furumura et al., 1998).采用的时间步长为Δt=0.015s.为了有效地观察尾波特征,需要确保模拟月震图的持续时间足够长.我们模拟了长达3000 s的月震波传播过程,共200000个时间步.

为了考察上月壳散射对月震波尾波形成的影响,对上月壳分别设置了一个均匀模型和一个扰动模型(表 2),用于对月震波场进行对比模拟研究.在模型ModelA0中,整个上月壳为均匀的一层,月震波不会发生散射.在模型ModelA1中,上月壳包含了关联长度为2.0 km,扰动标准差为3%的随机速度扰动(P波和S波的速度扰动一样).这样的关联长度和扰动标准差设置,也是在地球岩石圈中模拟地震波的散射常见的设置(Furumura and Kennett, 20052008).对于关联长度a的选取我们经过了如下考量:首先,目前尚没人指出对于上月壳的不均匀应该采用多大的关联长度;其次,由于长期遭受天体撞击,上月壳的破碎可能比较严重,非均匀的尺度应该不会超过地球岩石圈的非均匀尺度,不宜选择比2.0 km更大的关联长度;再次,考虑到模拟中的节点间距,选择太小的关联长度意义不大;最后,根据Aki和Richards(1980)的讨论,当关联长度a和波长λ相当时,散射效应最强,考虑到模拟中可计算的最小波长为4.45 km,取a=2.0 km也是一个较为合适的选择.同样,目前也没有关于上月壳中速度扰动幅度的估计.所以对于扰动标准差ε,我们选取了3%这一在地震波散射模拟中较为常见的值,用于模拟上月壳对月震波的散射.至于真实情况下,上月壳中速度扰动的标准差应该为多大,将在第4节给出具体的估计.

表 2 两个上月壳模型:ModelA0为均匀模型,ModelA1为扰动模型Table 2 Two models for upper Moon crust: ModelA0 is a uniform model,ModelA1 is a fluctuant model

3.1 上月壳均匀模型和扰动模型的波场快照对比

图 5展示了位于71 km深的浅源月震所激发的月震波在210 s,600 s和1200 s的波场快照.P波和S波的贡献分别由波场的散度和旋度计算而得(Furumura et al., 1998).

图 5 深度为71 km的月震在210 s,600 s,1200 s的波场快照
左边一列对应上月壳均匀模型(ModelA0),右边一列对应上月壳扰动模型(ModelA1).P波和S波分别由红色和绿色表示.
Fig. 5 Wavefield snapshots at 210 s,600 s,1200s for a 71 km deep moonquake
Left panels are for ModelA0 with a uniform upper Moon crust,right panels are for ModelA1 with velocity fluctuations in upper Moon crust. P and S waves are shown in red and green colors.

在210 s,对于均匀上月壳模型ModelA0,图 5a中P波和S波的波前具有十分清晰的轮廓,低速上月壳多次反射所引起的混响波也具有规则的轮廓.而对于上月壳含随机速度扰动的模型ModelA1,图 5b中波前轮廓的边界有些模糊,特别是从上月壳表面反射回来的S波(第四道绿色波前),轮廓边界十分模糊.在这之后的混响S波则完全看不到,取而代之的是上月壳速度扰动所造成的散射S波能量,在图 5b中显示为大量不规则分布的绿色散点.而对于P波的混响波,仍能看到一道道的轮廓,只是比图 5a中的轮廓模糊,同时我们还注意到P波后面的散射波能量(红色散点)并不像S波的散射波能量那样强烈.这说明上月壳对P波的散射要比对S波的散射弱,也就是说,月震波的尾波主要由S波的散射引起.另外,由于月震波在上月壳发生了散射,图 5b中月壳里的波场能量要比图 5a中强,特别是震源上方那部分区域,差别十分明显.在600 s,图 5c中扇形左半部分月壳已看不见什么能量,而图 5d中整个上月壳里仍然存有明显的波场能量,这归因于上月壳随机速度扰动对月震波造成的散射效应.在1200 s,图 5e中几乎已看不到能量,而图 5f中整个上月壳仍然广泛分布着可见的能量.以上关于浅源月震的波场快照表明,上月壳速度结构的随机扰动会对月震波造成持续强烈的散射,所引起的散射波能量将 在上月壳之中长时间存在,一直持续到1200 s以后.

图 6是867 km深的深源月震所激发的月震波在390 s,600 s和1200 s的波场快照.在390 s,由于上月壳非均匀结构对月震波的散射,图 6b中的波前轮廓不如图 6a中清晰,并且在P波和S波后面跟随着许多红色散点和绿色散点,它们分别是散射P波和散射S波.与图 6a相比,第三道S波(从月球表面反射回来的下行S波)后面的混响S波被散射S波取代,P波的部分混响波也被散射P波取代.与浅源月震的波场(图 5a,5b)类似,由于上月壳中散射的影响,图 6b中月壳里的波场能量要明显强于图 6a.在600 s,图 6c中月壳里的波场能量已经很微弱,而图 6d中整个月壳里仍存在很强的波场能量.在1200 s,图 6e中已完全看不到波场能量,而图 6f中仍然能够看到波场能量存在于上月壳中.与浅源月震十分类似,对于深源月震模拟出的波场快照也表明,上月壳里的散射会使月震波持续强烈地存在于上月壳之中,一直到1200 s以后.

图 6 深度为867 km的月震在390 s,600 s,1200 s的波场快照 左边一列对应上月壳均匀模型(ModelA0),右边一列对应上月壳扰动模型(ModelA1).P波和S波分别由红色和绿色表示. Fig. 6 Wavefield snapshots at 390 s,600 s,1200s for an 867 km deep moonquake Left panels are for ModelA0 with a uniform upper Moon crust,right panels are for ModelA1 with velocity fluctuations in upper Moon crust. P and S waves are shown in red and green colors.

对比模型ModelA0和模型ModelA1,浅震和深震的模拟结果都表明,当上月壳速度结构存在随机扰动变化时,月震波会在里面发生持续强烈的散 射,使得上月壳中长时间充斥着散射波能量,这在月震记录上将表现为强烈持久的尾波(参见图 7图 8). 3.2上月壳均匀模型和扰动模型的模拟月震图对比

图 7图 8分别显示了浅源月震和深源月震的合成位移月震波形以及P波和S波的到时.我们模拟的区域是三维月球模型的一个扇形截面,对应的模拟波形包含水平和垂直两个分量,其中水平分量与实际记录的径向分量(R分量)相对应,垂直分量与实际记录的垂直分量(Z分量)相对应.在模拟中,波的振幅近似地遵循1/ √ r的几何扩散衰减,r是接收点到震源的距离.为了更清晰地研究诸如散射等非几何扩散因素对波形的影响,我们给每条模拟记录乘上因子√ r以简单地补偿几何扩散.模拟中的震源是一个二维线源,在辐射效果上不同于三维点源,为此需要进行二维线源到三维点源的辐射校正.参考Helmberger和Vidale(1988)提出的针对二维合成地震图的校正公式,以及Wang等(20012013)和Li等(2014)所采用的公式形式,我们采用下面的校正公式:

图 7 深度为71 km的月震在几个震中距的模拟位移月震波形
上面一行图对应模型ModelA0,下面一行图对应模型ModelA1.左边一列图为R分量,右边一列图为Z分量.
Fig. 7 Synthetic displacement seismograms at several epicentral distances for a 71 km deep moonquake
Upper panels are for ModelA0,lower panels are for ModelA1. Left panels show radial components,right panels show vertical components.

图 7是位于月表下71 km的浅源月震在不同模型下的模拟位移月震波形.每张图上分别记录了震中距为11.25°,22.5°,33.75°,45.0°,56.25°,67.5°,78.75°,90.0°处的模拟月震波形.上下两行图显示了上月壳均匀模型(ModelA0)和上月壳为散射层模型(ModelA1)在波形记录上的巨大差别.对于上月壳均匀模型,不论是R分量还是Z分量的记 录,S震相后面信号的振幅衰减都比较迅速,在1300 s 之后,各个震中距所记录到的信号振幅都近乎为零.而当上月壳为散射层时,不论R分量和Z分量的记录,S震相后面都出现了强烈又持久的尾波信号.尾波在强弱上与S震相处于同一水平,尾波的振幅衰减十分缓慢,尾波信号可一直持续到3000 s.图 8是位于月表下867 km的深源月震在模型ModleA0和模型ModelA1下的模拟位移月震波形.可以看到上月壳散射对月震波形的影响也十分明显.对于模型ModleA0,R和Z分量的记录在约900 s之后就近乎为零,而对于模型ModelA1,尾波信号一直持续到3000 s,并且十分强烈.相比于浅源月震而言,对于深源月震,散射所引起的尾波效应甚至更强.在图 8中我们看到,散射所引起的尾波信号的振幅有时比直达S波的振幅还大,比如在45.0°,56.25°,67.5°,78.75°这几个震中距所记录的月震波形.这说明上月壳非均匀速度结构所引起的散射对月震波形记录的影响是可以十分巨大的.

图 8 深度为867 km的月震在几个震中距的模拟位移月震波形
上面一行图对应模型ModelA0,下面一行图对应模型ModelA1.左边一列图为R分量,右边一列图为Z分量.
Fig. 8 Synthetic displacement seismograms at several epicentral distances for an 867 km deep moonquake
Upper panels are for ModelA0,lower panels are for ModelA1. Left panels show radial components,right panels show vertical components.

以上对浅源和深源月震的模拟波形的分析,一致表明上月壳散射可在月震信号记录中产生长时间持续的强烈尾波.这一点与图 5图 6的波场快照所展示的结果也是一致的.

3.3 上月壳均匀模型和扰动模型的模拟月震记录与真实月震记录对比

本节将模拟月震图与阿波罗月震图进行比较.从Wang等(2013)所采用的月震数据中挑选出浅源事件SH1和深源事件A1在台站S12的记录作为用于比较的观测数据,分别命名为SH1S12和A1S12.图 9展示了事件SH1,事件A1,以及台站S12在月球上的位置.表 3列出了所选记录的震源深度、震中距、方位角三个参数.每组记录都包含了X,Y,Z三个长周期分量(LPX,LPY和LPZ).

表 3 用于比较的阿波罗月震数据的一些参数Table 3 Parameters of Apollo seismic data used in the comparison

图 9 浅源月震(SH1),深源月震(A1)和记录台站(S12)在月球上的分布Fig. 9 Distribution of shallow moonquake(SH1),deep moonquake(A1) and station(S12)on the Moon

在进行比较之前,我们遵循Bulow等(2005)Garcia等(2011)Wang等(2013)所采用的步骤对原始的观测数据进行处理.首先,对数据进行从0.01 Hz到3.31 Hz(奈奎斯特频率)的带通滤波.然 后,采用稳健滑动中位数去尖峰算法(robust running median despiking algorithm)去掉每条记录中的过大尖峰(可能由噪声和热异常引起).又由于X和Y两 个方向的传感器可能具有不同的增益,所以接着采用Garcia等(2011)所估计的X/Y增益比对Y分量记录进行校正.最后,将X分量和经过校正的Y分量旋转为径向(R)分量和切向(T)分量以便和模拟月震图比较.考虑到能够计算的最大频率为0.4 Hz,在进行比较之前,对于观测数据和模拟月震记录,统一选取从0.05 Hz到0.4 Hz的频率范围进行带通滤波.模拟月震记录在滤波之前都已采用公式(13)进行过震源辐射的校正.

图 10展示的是两个模型的模拟月震图与阿波罗月震图的比较,所提取的模拟月震图与观测数据具有相同的震中距.对比三道波形,我们看到,简单的均匀分层模型(ModelA0)不足以解释阿波罗月震记录中的尾波特征,而把非均匀结构对月震波的散射效应包含进上月壳中后(ModelA1),模拟月震图中出现了与阿波罗月震记录十分相似的尾波.不论是浅源月震还是深源月震,模拟月震图尾波和阿波罗月震图尾波的强弱和持续时间都处在同一水平,尾波振幅的衰减趋势也比较一致.这些都表明上月壳非均匀结构对月震波的散射是解释月震信号尾波形成的有效机制.

图 10 浅源和深源月震的阿波罗月震图与模拟月震图比较Fig. 10 Comparison between Apollo seismograms and synthetics for both shallow and deep moonquakes

本节探讨了上月壳非均匀结构的散射对于在月震波形中引起尾波的有效性.通过模拟我们发现,在均匀分层模型基础上,进一步考虑上月壳非均匀结构在月震波传播过程中所引起的散射效应,能有效地解释月震信号的尾波特征.下面,让我们试着探讨一下尾波形成的物理机制.首先,上月壳波速很低(模拟中上月壳的S波速度为1.8 km·s-1,下月壳的S波速度为3.2 km·s-1),波阻抗低,月震波在这一层中的混响(Wang et al., 2013)有效地保存了波的能量.其次,上月壳的滞弹性衰减非常低或者说Q值非常大(模拟中上月壳的Q值为6750),使得月震波的能量在其中的衰减缓慢,保证了月震波在上月壳中长时间持续的可能性.最后,也最为重要的是,在前两个原因的基础上,散射层对月震波造成的多次散射使月震波不断散焦和聚焦同时,延长了散射月震波在上月壳内的传播时间,使得在直达波之后不断地有散射波到达接收点,从而在月震信号记录中引起长时间持续的强烈尾波信号.综合以上分析,上月壳中的低波速,高Q值和散射这三个因素的共同作用造就了月震记录中强烈持久的尾波信号. 4 速度扰动幅度对尾波的影响

采用自相关函数描述上月壳速度结构的非均匀性,不同的速度扰动标准差(ε)会对应不一样的散射效果,相应的月震尾波信号也会呈现出不同的波形特征.本节讨论速度随机扰动的标准差对月震尾波信号的影响.为此设置了五个不同的非均匀上月壳模型,分别对应于五个不同的速度扰动标准差ε=1%,3%,5%,7%,9%,如表 4所示.其中,各个模型的关联长度a都设置为2.0 km,如在第3节中所述这是一个相当合理的值.采用与第3节中一样的震源设置、节点间距和时间步长分别对上述五个模型进行模拟.对由这五个模型计算出的模拟波形,采用适当的数学方法进行分析,可获知扰动标准差ε对尾波形状的影响.4.1节研究尾波衰减随扰动标准差ε的变化,4.2节将查看不同扰动标准差下的模拟月震图与真实月震记录的相似度,挑选出相似度最高的模型,并进一步给出关于上月壳中真实的速度扰动幅度的估计.

表 4 不同标准差的上月壳扰动模型Table 4 Fluctuant models of upper Moon crust with different st and ard deviations

4.1 尾波的特征衰减时间随扰动标准差ε的变化

尾波的衰减可简单地用下面的解析式描述(Shearer,2007):

其中,A(t)表示t时刻的振幅,A0表示初始振幅,τ叫做特征衰减时间,是尾波振幅衰减到原来的1/e所需的时间.τ值的大小反映了尾波衰减的快慢,τ越大,尾波衰减越慢,τ越小,尾波衰减越快.为了获得尾波的特征衰减时间τ,需要用(14)式所示的指数衰减函数去拟合地震图中的尾波部分.

采用如下步骤计算模拟月震图的特征衰减时间:(1)采用希尔伯特变换求取模拟月震图的包络函数;(2)采用宽为300 s的滑动时间窗平滑包络函数;(3)采用(14)式所示的指数衰减函数对平滑后的包络函数进行最小二乘拟合,计算出特征衰减时间τ(Blanchette-Guertin et al., 2012a).拟合的起点通过手动的方式选取,考虑到指数衰减并不总是从最大振幅处开始.拟合的终点选取采用如下的标准:先求出尾波在2000~2750 s这段时间内的最小振幅Amin,再求出整段待拟合尾波的最大振幅Amax,然后选取Amax/3和2Amin中的最小者作为终点标准,尾波振幅首次下降到这一最小值的地方即为拟合的终点(可确保振幅衰减至1/e以下,使拟合时间窗的长度大于τ).上述关于拟合起点和终点的选取基于一定的经验,其中终点的选取过程可通过程序自动 实现.图 11摘录了求取特征衰减时间τ的大体步骤.

图 11 求取特征衰减时间的步骤展示Fig. 11 Example of steps for calculating the characteristic decay time

图 12显示了22.5°,45.0°,67.5°,90.0°这几个震中距处的特征衰减时间τ随扰动标准差ε的变化.这四幅图中的τ值大致呈现出如下变化趋势:当ε小于3%时,τε的增大而增大;在ε大于3%时,τε的增大而减小.上述变化趋势对于浅源月震较为明显,对于深源月震也大体符合.ε值的大小反映了随机速度扰动的剧烈程度.当ε较小时,速度扰动较小,散射相对也较弱,特别地,当ε→0时,介质变为均匀介质,散射将不复存在.当ε较大时,速度扰动大,散射也将变得强烈.τε先增大后减小表明,要形成长时间持续的尾波需要适中的散射强度,散射太弱或太强都会导致尾波快速衰减.若散射过弱(ε<3%),则不足以在记录中产生足够强的尾波,尾波振幅很难保持在一定大小之上;若散射过强,波在非均匀结构中的聚散焦过于强烈,也会缩短尾波的持续时间.模拟结果表明,当上月壳中速度扰动的标准差在3%左右时,最利于在月震记录中引起长时间持续的尾波,这一规律在浅源月震中尤其明显.

图 12 几个震中距处的特征衰减时间随扰动标准差的变化
上面一行对应71 km深的月震,下面一行对应867 km深的月震.左右两列分别对应于R和Z两个分量.
Fig. 12 Dependence of the characteristic time on st and ard deviation of velocity fluctuations at several epicentral distances
Upper panels are for a 71 km deep moonquake,lower panels are for an 867 km deep moonquake. Left and right panels are for radial and vertical components,respectively.
4.2 确定上月壳中速度扰动的标准差

本节将通过比较不同扰动标准差下的模拟月震图与阿波罗月震图的相似度,确定出上月壳中最有可能的速度扰动大小.

图 13展示了五个模型的模拟月震图与阿波罗月震图的比较,所选的阿波罗月震记录以及数据的处理过程与第3节中相同.上下两行图分别对应浅源和深源月震,左右两列图分别是R分量和Z分量.在每张图中,最上面一道波形为真实月震数据,下面五道波形为模拟月震图,从上到下分别对应模 型ModelB1,ModelB2,ModelB3,ModelB4,ModelB5,相应的扰动标准差分别为1%,3%,5%,7%,9%.为了便于比较波形,每个图中的模拟月震图都已按振幅的标准差进行了归一化.对于模型ModelB1,上 月壳散射所引起的尾波相对较弱,虽然在原有P波和S波的后面带上了尾波信号,但尾波信号的振幅 明显比S波的振幅小.到模型ModelB2时,尾波信号在整体上明显增强,振幅达到与S波振幅同一水平,尾波信号紧跟在S震相后面十分平缓地衰减.从ModelB2到ModelB5,四个图中尾波的衰减都依次变快,其中ModelB5的尾波衰减尤为迅速.这五个模型的模拟月震图在波形上呈现出的变化趋势,与4.1节中基于特征衰减时间的分析所得的结果完全一致.从波形上看,模型ModelB2产生的模拟月震图与真实月震图最为接近,表明上月壳中速度扰动的标准差可能比较接近于3%.

图 13 阿波罗月震图与不同速度扰动标准差模型的模拟月震图比较Fig. 13 Comparison between Apollo seismograms and synthetics for models with different st and ard deviations of velocity fluctuations

图 14给出了图 13中所示月震图平滑后的包络线,平滑所采用的时间窗宽度为100 s.包络线上每个点所处的高度反映了该时刻尾波振幅的大小.包络线的比较也显示出与波形相一致的结果.在四张图中,模型ModelB1的月震图包络线的振幅在约1500 s之后就变得很小,反映到月震图上就是尾波信号在1500 s之后很弱.从ModelB2到ModelB5,每张图中尾波一段的包络线下降速度都依次变快,进一步表明了尾波衰减从ModelB2到ModelB5依次变快这一趋势.五个模型中,只有ModelB2的月震图包络线在尾波一段衰减最为平缓,并且振幅能长时间保持在不太低的水平,这两点 与真实月震图的包络线最为接近.包络线的比较也表明上月壳中速度扰动的标准差可能比较接近于3%.

下面,采用互相关系数来量化模拟月震图与真实月震图的相似程度,以便精确地挑出最佳模型.用x(t)y(t)分别表示两个月震信号,它们的互相关系数corr定义为:

该式的离散形式为:

xi和yi分别表示对应的两个离散信号序列.

我们先统一采用0.2 s的时间间隔(15000个样本点)对真实月震图和模拟月震图的包络线重采样,然后采用公式(16)计算互相关系数.对每个模型的每条模拟月震图,都可以算出一个互相关系数.图 15显示了浅源月震和深源月震的两个分量的互相关系数R随扰动标准差ε的变化情况.图中的四条线均显示,ε=3%(ModelB2)时的模拟月震图包络线与真实月震图包络线之间的互相关系数最大,都在约93%以上,ε=5%(ModelB3)时的互相关系数次之,其余三个扰动标准差下的互相关系数都比这两者的小.如果连续的看,互相关系数最大的地方应该在3%到5%之间并且更靠近3%的地方.这表明上月壳中速度扰动的标准差应该在3%到5%之间,而且很可能就在3%附近.

图 14 图 13中所有月震图的包络线 每条包络线都已采用100 s的时间窗进行了平滑. Fig. 14 Envelope functions of seismograms shown in Fig. 13 Each envelope was smoothed using a 100 s running window.

图 15 互相关系数随扰动标准差的变化
红色数据点对应71 km浅源月震,绿色数据点对应867 km 深源月震.圆圈代表R分量的互相关系数,正方形代表Z 分量的互相关系数.
Fig. 15 The change of correlation coefficient with st and ard deviation of velocity fluctuations
Red markers are for shallow moonquake(71 km),green markers are for deep moonquake(867 km). Circles represent radial components,squares represent vertical components.

5 结论与讨论

本文采用伪谱和有限差分混合方法成功模拟了上月壳中的散射在月震记录中引起尾波的现象.与前人在月震波散射模拟中所采用的基于射线理论的声子法相比,我们将上月壳中的速度扰动直接包含到波动方程当中,计算出的波形能包含更加完整和准确的信息.诸如多次散射波、界面反射波和转换波、衍射波、波的聚散焦等各种机制引发的效应都能被有效地包含到模拟波形当中.因此,我们的模拟具有很高的可靠性.

我们发现上月壳中速度扰动所引发的散射能有效地在月震记录中引起长时间持续的尾波.不论是浅源月震还是深源月震,在均匀分层模型的基础上,进一步考虑上月壳中的小尺度非均匀结构对月震波的散射,都能有效地解释在月震记录中广泛存在的尾波.对于尾波形成的物理机制,结合上月壳中的波速以及Q值进行分析,我们认为月震记录中的尾波可能是由上月壳中的低波速、高Q值和散射这三个因素的共同作用所引起.接着,我们研究了上月壳中的速度扰动幅度对尾波衰减快慢的影响.我们发现,随着扰动标准差的逐渐增大,尾波的衰减呈现出先变慢后变快的变化规律.从采用1%,3%,5%,7%,9%这几个标准差模拟所得到的结果看,当上月壳中速度扰动的标准差为3%时,月震记录中的尾波衰减最慢,对应了更长的尾波持续时间.最后,为了给出对上月壳中真实的速度扰动标准差的估计,从波形、包络线和互相关系数这三个角度比较了不同扰动标准差下的模拟月震图与真实月震图的相似程度.从比较的结果看,上月壳中速度扰动的标准差应该在3%到5%之间,而且在3%附近的可能性比较大,因为3%的扰动标准差所模拟出的月震图与真实月震图最接近.

本文通过模拟上月壳中的速度扰动所造成的月震波散射,解释了月震记录中的尾波的成因,并给出了关于上月壳速度扰动的标准差的合理估计.月震波的散射,不但受速度扰动的幅度影响,也受速度非均匀的尺度影响(在自相关函数中由关联长度a控制).在本文的研究中,我们结合网格节点间距、上月壳的构造情况以及散射理论选择出了一个较为合适的非均匀尺度用于模拟当中.这样的选择具有一定的合理性,不过,非均匀尺度对于月震波散射的影响仍然是一个在未来的模拟工作中值得研究的问题.

我们的模拟能计算的最高频率为0.4 Hz,真实的月震长周期记录的频率可达2 Hz.未来采用并行模拟将更多的高频成分包含到计算结果中,可在更大的频率范围内进一步研究月震波的散射在月震尾波的形成中所起的作用.模拟中,求解区域为通过大圆从三维月球模型中截取的一个扇形平面.与三维月球模型相比,二维扇形模型不能包含来自平面之外的结构所引起的散射和聚散焦等效应.将来在计算机硬件条件成熟的情况下,采用三维月球模型来模拟月震波的散射可能会得到更好的结果.

致 谢 感谢北京大学地球与空间科学学院地球物理系王桥同学为本研究的计算所提供的帮助.
参考文献
[1] Aki K, Richards P G. 1980. Quantitative Seismology Volume Ⅱ. San Francisco: W.H. Freeman.
[2] Blanchette-Guertin J F, Johnson C L, Lawrence J F. 2012a. Investigation of scattering in lunar seismic coda. J. Geophys. Res.: Planets (1991—2012), 117(E6), doi: 10.1029/2011JE004042.
[3] Blanchette-Guertin J F, Johnson C L, Lawrence J F. 2012b. Modeling seismic waveforms in a highly scattering Moon.// 43rd Lunar and Planetary Science Conference, Abstract 1473, Lunar and Planetary Institute, Houston.
[4] Blanchette-Guertin J F, Johnson C L, Lawrence J F. 2013. Effect of variable scatterer length-scales and frequency dependent attenuation on the decay of lunar seismic coda.// 44th Lunar and Planetary Science Conference, Abstract 1234. Houston: Lunar and Planetary Institute.
[5] Bulow R C, Johnson D L, Shearer P M. 2005. New events discovered in the Apollo lunar seismic data. J. Geophys. Res., 110, E10003, doi: 10.1029/2005JE002414.
[6] 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.
[7] Dainty A M, Toksöz M N, Anderson K R, et al. 1974. Seismic scattering and shallow structure of the moon in oceanus procellarum. The Moon, 91(1-2): 11-29.
[8] Dziewonski A M, Anderson D L. 1981. Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4): 297-356.
[9] Frankel A, Clayton R W. 1986. Finite difference simulations of seismic scattering: Implications for the propagation of short-period seismic waves in the crust and models of crustal heterogeneity. J. Geophys. Res.: Solid Earth (1978—2012), 91(B6): 6465-6489.
[10] Furumura T, Kennett B L N, Furumura M. 1998. Seismic wavefield calculation for laterally heterogeneous whole earth models using the pseudospectral method. Geophys. J. Int., 135(3): 845-860.
[11] Furumura T, Kennett B L N. 2005. Subduction zone guided waves and the heterogeneity structure of the subducted plate: Intensity anomalies in northern Japan. J. Geophys. Res., 110, B10302, doi: 10.1029/2004JB003486.
[12] Furumura T, Kennett B L N. 2008. A scattering waveguide in the heterogeneous subducting plate.// Sato H, Fehler M eds. Scattering of Short-Period Seismic Waves in Earth Heterogeneity, Advances in Geophysics. Vol.50, Chap. 7. Elsevier, 195-217.
[13] Garcia R F, Gagnepain-Beyneix J, Chevrot S, et al. 2011. Very preliminary reference Moon model. Physics of the Earth and Planetary Interiors, 188(1-2): 96-113, doi: 10.1016/j.pepi.2011.06.015.
[14] Goins N R, Dainty A M, Toksöz M N. 1981. Seismic energy release of the Moon. J. Geophys. Res., 86(B1): 378-388.
[15] Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bulletin of the Seismological Society of America, 86(4): 1091-1106.
[16] Helmberger D V, Vidale J E. 1988. Modeling strong motions produced by earthquakes with two-dimensional numerical codes. Bulletin of the Seismological Society of America, 78(1): 109-121.
[17] Herrmann R B. 1979. SH-wave generation by dislocation sources—a numerical study. Bulletin of the Seismological Society of America, 69(1): 1-15.
[18] Jiang M M, Ai Y S. 2010. Moonquakes and lunar interior. Geochimica (in Chinese), 39(1): 15-24.
[19] Kennett B L N, Furumura T. 2013. High-frequency Po/So guided waves in the oceanic lithosphere: I—long-distance propagation. Geophys. J. Int., 195(3): 1862-1877.
[20] Khan A, Maclennan J, Taylor S R, et al. 2006. Are the Earth and the Moon compositionally alike? Inferences on lunar composition and implications for lunar origin and evolution from geophysical modeling. J. Geophys. Res., 111, E05005, doi: 10.1029/2005JE002608.
[21] Koyama J, Nakamura Y. 1980. Focal mechanism of deep moonquakes.// Bedini S A ed. Lunar and Planetary Science Conference Proceedings, Vol. 11. 1855-1865.
[22] Latham G, Ewing M, Dorman J, et al. 1972. Moonquakes and lunar tectonism. The Moon, 4(3-4): 373-382.
[23] Latham G V, Ewing M, Press F, et al. 1970. Passive seismic experiment. Science, 167(3918): 455-457.
[24] Lawrence J F, Johnson C L. 2010. Synthetic seismograms with high-frequency scattering for the Moon.// 41st Lunar and Planetary Science Conference, Abstract 2701. Houston: Lunar and Planetary Institute.
[25] Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425-1436.
[26] Li D Z, Helmberger D, Clayton R W, et al. 2014. Global synthetic seismograms using a 2-D finite-difference method. .Geophys. J. Int, 197(2): 1166-1183.
[27] McGetchin T R, Settle M, Head J W. 1973. Radial thickness variation in impact crater ejecta: Implications for lunar basin deposits. Earth Planet. Sci. Lett., 20(2): 226-236.
[28] Nakamura Y. 1977. Seismic energy transmission in an intensively scattering environment. J. Geophys. Res., 43(1-2): 389-399.
[29] Nakamura Y. 1978. A1 moonquakes-source distribution and mechanism.// Lunar and Planetary Science Conference Proceedings, Vol.9.3589-3607.
[30] Nakamura Y. 2005. Farside deep moonquakes and the deep interior of the Moon. J. Geophys. Res., 110(E1): E0100, doi: 10.1029/204JE002332.
[31] Nakamura Y, Latham G V, Dorman H J, et al. 1979. Shallow moonquakes-depth, distribution and implications as to the present state of the lunar interior.// Lunar and Planetary Science Conference Proceedings, Vol. 10. 2299-2309.
[32] Ouyang Z Y. 2005. Introduction to Lunar Science (in Chinese). Beijing: China Astronautic Publishing House.
[33] Sato H, Fehler M C. 1998. Seismic Wave Propagation and Scattering in the Heterogeneous Earth. New York: Springer-Verlag.
[34] Shearer P M, Earle P S. 2004. The global short-period wavefield modelled with a Monte Carlo seismic phonon method. Geophys. J. Int., 158(3): 1103-1117.
[35] Shearer P M. 2007. Deep Earth structure—Seismic scattering in the deep Earth.// Seismology and the Structure of the Earth, Treatise Geophys., Vol. 1. Amsterdam: Elsevier, 695-729.
[36] Thompson T W, Campbell B A, Ghent R R, et al. 2009. Rugged crater ejecta as a guide to megaregolith thickness in the southern nearside of the Moon. Geology, 37(7): 655-658, doi: 10.1130/G25565A.1.
[37] Tromp J. 2007. Forward modeling and synthetic seismograms: 3D
[38] numerical methods.// Romanowicz B, Dziewonski A eds.Treatise on Geophysics. Amsterdam: Elsevier, 191-217.
[39] Wang Y B, Takenaka H, Furumura T. 2001. Modelling seismic wave propagation in a two dimensional cylindrical whole earth model using the pseudospectral method. Geophys. J. Int., 145(3): 689-708.
[40] 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. Geophys. J. Int., 192(3): 1271-1287, doi:10.1093/gji/ggs094.
[41] Weber R C, Lin P Y, Garnero E J, et al. 2011. Seismic detection of the lunar core. Science, 331(6015): 309-312.
[42] Witte D C, Richards P G. 1990. The pseudospectral method for simulating wave propagation. Computational Acoustics, 3: 1-18.
[43] Zhan W, Li F. 2007. The inner structure of the moon. Progress in Geophysics (in Chinese), 22(3): 737-742.
[44] 姜明明,艾印双. 2010. 月震与月球内部构造. 地球化学,39(1):15-24.
[45] 欧阳自远. 2005. 月球科学概论. 北京:中国宇航出版社.
[46] 占伟,李斐. 2007. 月球内部构造研究综述. 地球物理学进展,22(3):737-742.