ACTA Scientiarum Naturalium Universitatis Pekinensis
全火星模型中震波传播特征的数值模拟研究
邓迪 肖万博 王彦宾†
北京大学地球与空间科学学院地球物理学系, 北京100871; † 通信作者, E-mail: ybwang@pku.edu.cn
摘要 采用基于交错网格的傅里叶伪谱与有限差分混合方法, 求解弹性波动方程, 根据地球化学分析得到的两个火星理论结构模型, 模拟二维全火星模型中P-SV波和SH波的传播过程。根据理论地震图和波场快照,讨论全火星模型中震波的传播过程以及各种震相的产生和演变, 分析模型内部火星壳厚度以及火星核幔边界深度对震波传播的影响。结果表明, 在低速火星壳内部多重反射波及转换波的相干叠加会形成很强的波列,其特征受火星壳厚度的影响较大, 在切向分量上可以更清晰地观测到核幔边界的反射震相。关键词 全火星模型; 震波传播; 数值模拟; 伪谱与有限差分混合方法
第 56 卷 第 4 期2020 年 7 月
器于 年 月 日登陆火星上的乌托邦平原, 着陆器上的地震仪连续工作19个月, 但除在第80个火星日记录到一个可能的火星震事件外, 没有其他发现[3]。该事件发生时, 没有记录到风数据, 因此不能排除风噪声的干扰。没有记录到其他事件的原因可能是地震仪在远震体波的频带宽度中灵敏性不足[4],也可能是由于仪器的位置导致其对风噪声的高灵敏度[5‒6]。陨石撞击是引起行星内部震波传播的另一个潜在来源。在阿波罗计划1969 — 1977年对月球的8年观测中, 4个台站观测到13000多个月震事件, 其中包含约1800多个陨石撞击事件[7]。在月球上, 陨石直接落在月球表面并产生月震信号, 由于缺乏大气层, 因此仅通过表面位移就可以检测到所有的撞击。但是, 火星的大气层能非常有效地防止陨石撞击, 到达火星表面的10 kg以下陨石的数量仅为入射到大气层陨石数量的10%, 因此在火星上检测到陨石撞击的数量很少[3]。对火星核结构的研究目前主要通过重力观测进行。Yoder等[8]通过火星全球探勘者号(Mars Global Surveyor)的无线电跟踪分析, 测得火星太阳引力下的潮汐变形勒夫数k2。观察得到的k2 值为 0.153± 0.017, 排除了火星核为固体铁芯的可能性, 至少外核是液体。之后, 研究人员对2001年火星奥德赛号(Mars Odyssey)、2003年机遇号火星漫游车(Mars Exploration Rover Opportunity)和 2005年火星勘测轨道飞行器计划(Mars Reconnaissance Orbiter)发回的跟踪数据进行一系列的分析, 确定k2值为 0.169± 0.006[9‒10]。
Sohl等[11]提出火星内部结构的两个模型, 他们在流体静力平衡和封闭传热的条件下进行求解, 得到质量、静水压力、重力、温度和热流密度的径向分布, 还获得火星模型的分层密度和火星震波速结构。两个模型分别是满足地球物理约束(即极惯性矩达到最大可能值)的模型A以及与来自火星的SNC陨石的地球化学约束一致的模型B。
2018 年 11月 26日, 洞察号火星地球物理探测器降落在埃律西昂平原, 并在火星表面部署地震实验仪器(Seismic Experiment for Internal Structure, SEIS), 用于研究火星内部结构。地震仪由长周期三轴宽频带仪和三轴短周期仪组成, 覆盖 0.01~50 Hz的频率范围。与海盗号地震仪相比, SEIS检测火星震的分辨率有所提升。另一项与海盗号不同的重
630大改进是, 地震仪通过机械臂直接部署在火星表面,并通过高效的隔热和隔风保护, 使其免受温度和风变化的影响。基于现有的火星知识, 预计 SEIS 每年可检测到几十次震中距40°以内、矩震级大于3级的火星震[12]。2019 年4月 23 日, 法国国家空间研究中心(Centre National d’etudes Spatiales, CNES)宣布洞察号首次在火星上检测到发生于2019 年4月6日的火星震信号。目前, 科学家仍在研究该火星震的数据, 以期确认其震源信息[13]。在洞察号发射之前, 针对可能观测到的火星震数据, 人们利用地震学方法开展了相关研究, 如火星震的活动性、火星震定位、震波传播、火星浅层结构和壳幔结构等, 为洞察号的设计与发射以及火星震观测数据的分析奠定了基础[14‒20]。地震波数值模拟研究在地震学中起到非常重要的作用, 有助于理解复杂地震波的传播过程和解释
[21]地震观测数据。Furumura 等 将傅里叶伪谱法与有限差分法相结合, 对 1999年台湾集集地震进行三
[22]维强震地面运动的数值模拟。马德堂等 运用傅里叶伪谱法与有限元法的混合方法, 对弹性波动方
[23]程进行求解。魏星等 运用傅里叶伪谱法与四阶精度有限差分方法的混合方法, 基于傅里叶伪谱法精度高、内存消耗少的前提, 发挥四阶精度有限差分方法对人工边界和自由界面的处理优势, 对二维
[24]模型进行弹性波场的数值模拟计算。秦艳芳等将魏星等的方法推广到三维非均匀介质地震波传播模拟计算中, 在确保精度和效率的情况下, 成功地实现对三维沉积盆地模型的并行模拟计算。Wang
[25]等 和 Jiang 等[26]应用交错网格傅里叶伪谱法与有限差分法的混合方法求解弹性动力学方程, 分别模拟计算全月球模型中P-SV波和SH波的传播。本文对二维全火星模型中火星震波传播进行数值模拟, 针对已有的理论火星结构模型, 模拟火星震波传播的理论地震图和波场快照, 讨论火星壳厚度和核幔边界深度对火星震波传播的影响, 可为解释未来可能获取的火星震数据以及火星结构和火星震发生机制的研究提供方法参考。
1基于交错网格的傅里叶伪谱与有限差分混合方法和全火星模型
对于各向同性的弹性介质, 在柱坐标系(r, θ, z)下, 假设z方向上所有变量保持不变, 在(r, θ)平面内, 以速度和应力形式表示的P-SV波和SH波的二
邓迪等 全火星模型中震波传播特征的数值模拟研究
在横向上, 通过傅里叶伪谱法计算空间变量对的偏导数[27]: i l k d 1 k j i( l k)e N 1 2 d 2 2 l 0 lj F ( l k)ei2 , N其中, j=1, 2, …, N–1, f (j∆θ)表示在 θ方向上离散的空间变量, F(l∆k)表示其相应的波数域的傅里叶变换, N表示空间离散的网格节点数, ∆θ和∆k分别表示在θ方向上和波数域的离散间隔。在半径方向上, 空间变量 f(m∆r)对 r的偏导数由交错网格的四阶精度有限差分格式计算[28]: 3 3 f m r f m r 1 2 2 fm r 24 r 1 1 f m r f m r 9 2 2 8 r 。(4)在模型区域内共需要处理两个边界条件, 通过满足表面处牵引力为零的自由边界条件, 将自由表面引入模型中。对于模型的中心区域, 采用 Cerjan [29]等 提出的带有20个网格点的“缓冲区”的吸收边界条件。式(1)中的震源体力采用Wang 等[30‒31]给出的二维柱坐标的线源形式: 1 fr ( r , , t ) Mrr ( t ) ( r r 0) ( 0) r r M ( t )1 ( r r) r 0 rr0 d dr f ( ), 0 (5a)
第 56 卷 第 4 期2020 年 7 月 ( r r ) ( ) 0 0 r 0) ( ), 0 (5b)
( r r 0) ( 0)
r 0) ( ),
0其中, (r0, Mrr, Mθθ, Mrθ, Mzr 和Mzθ分别为震源矩张量的分量; δ(r)和 δ(θ)表示震源体力的空间分布方式, 采用 Herrmann[32]提出的伪 δ函数进行近似计算。
本文选取以往研究中广泛采用的Sohl 等[11]提出的火星内部结构参考模型A和模型B, 两个模型的火星壳厚度分别为110和 250 km, 火星幔上层与火星幔下层由橄榄石‒尖晶石过渡区分隔开, 火星幔下层又分为β尖晶石层和γ尖晶石层两层。因此,火星内部分为5 层, 由浅至深分别是火星壳、火星幔上层(α橄榄石层)、β尖晶石层、γ尖晶石层和火星核, 各层的P波速度、S波速度和密度分布如图2
[33]所示。本文的Q值模型采用 Toyokuni 等 给出的值作为参考, 即液体核心内的Q值为 1000000, 其他区域的 Q值均为 600。(5c)由于火星震的震源机制没有任何观测资料的约束, 本文在模拟过程中将典型的剪切地震震源作为参考。本研究的震源机制为二维双力偶线源, 取Mrr= –1.0, Mθθ=1.0, Mrθ=mθr=0, 相当于一个45°倾角的倾滑断层。模型的θ方向离散为2048个网格点,范围为 0~2π; r方向离散为512个网格点, 范围为0~3390 km。因此, r方向的网格间距为6.62 km, θ 方向的网格间距最小为0.02 km, 最大为10.40 km, 位
[34]于火星模型表面。根据Knapmeyer 等 构建的包含 13681个火星震的目录, 火星震事件的震源深度最大可达100 km, 大多数事件的震源深度小于60 km, 因此本文模拟的震源深度为60 km。震源时间函数的宽度为10 s, 计算中使用的时间间隔∆t由稳定性条件确定, 即受模型中最小网格间距∆min与最大波速Vmax之比的约束。时间步长∆t满足
min t
Vmax其中, α为常系数0.26。中心附近的∆min为 0.02 km, Vmax 为 10.5 km/s, 因此得出的∆t 为 0.0004 s。该值在实际计算中过小, 因此需要设法加大∆t。本文参考 Wang等[30]提出的波数域滤波算法, 解决邻近球心 r=0处横向网格间距过小的问题。在计算横向导数时应用一个平滑因子, 即应用低通滤波器滤除高频波数成分, 由此得到的∆t<0.049, 比滤波前的时间步长提高约100倍。本研究取∆t值为0.04 s, 计算的总步数为75000, 能得到3000 s的理论地震图。
在研究切向分量波传播时, 由于该分量上只有S波, 且 S波在火星液核内不传播, 因此对模拟参数做了调整, 将火星核内的网格点剔除。θ方向离散为2048个网格点, 范围为0~2π; r方向离散为550个网格点, 范围为 1190~3390 km, 因此r方向的网格间距为5.0 km, 同时r方向的网格能完整覆盖SH波的传播路径。
2.1
,
模型A的理论地震图与波场快照
(6)
图3展示模型A在火星表面的合成理论地震图以及主要震相的理论射线到达时间, 可以看到, 由于火星核的存在, P波传播到核幔边界时速度会明显减小, 并改变路径方向, 因此存在着直达P波的影区, 影区范围的大小与火星核半径有关。3个分量上均能观测到互相对应的直达S 波, 但只有径向和垂向分量能观测到直达P 波, 同时能观测到P波
变得更强。420 s 时, 来自核幔边界的直达SV波的反射波和转换P波向上传播, 由于P波影区的存在,核幔边界附近可以看到弯曲的绕射P波。720 s 时,可以看到整个波场存在各种“Y”字形P波和SV波的多重反射波, P波已经传播至火星核的另一面, 产生转换的SV波。
图 5展示SH波在全火星模型A中传播的位移波场快照, 能够很清晰地看到直达SH波、多重反射 SS 和 SSS震相以及传播到核幔边界反射的SCS震相。180 s 时, 直达SH波和经过表面反射的SH波向下传播, 波场非常清晰。420 s时, 直达SH波到达核幔边界, 并产生反射 SH 波向上传播, SCS和sscs震相非常清晰。630 s时, 由于影区的存在, 核幔边界附近产生明显的绕射波。750 s时, 反射SH波已传回地表, 能看到火星壳内传播的勒夫面波和低速火星壳内由于多重反射产生的很长的波列。
2.2模型B的理论地震图与波形快照
图6展示模型B在火星表面的合成理论地震图以及各个主要震相的理论射线到达时间。与模型A相比, 模型B的火星核半径更大, 因此直达P波的影区范围也更大。切向上, 由于火星核半径更大,直达SH波的传播距离更短, SCS震相的到时也比模型A更早。同时, 由于模型B的火星壳更厚, 低速区内SH波的传播距离更长, 因此 SCS 与 sscs震相的到时相差更大。二次反射的SCS2震相和三次反
633
射的 SCS3震相等也是如此。由于切向上只有SH波, 径向和垂向上P波和SV波会发生震相转换, 因此在切向上更容易估算波的传播距离和传播时间。
634
Fig. 8
[4] [5]