Numerical Modeling of Global Seismic Wave Propagation in the Whole Mars Models

DENG Di, XIAO Wanbo, WANG Yanbin†

Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871; † Corresponding author, E-mail: ybwang@pku.edu.cn

北京大学学报(自然科学版) 1976 9 3

维波动方程为 vr 1  1   ( r )  r rr t r r r  r v 1  1   r 2   t r( ) r r2  r   vz  1    zr z    zr t r r  对于各向同性的弹性介质, 其本构关系为  v   v , rr r   (  2 )  v  r t  r r  v    vr  2  v       , r  t r  r      vr , v  r    v   t  rr     v    , zr z  t r    vz z   t r 其中, vr, vθ 和 vz分别为垂向、径向和切向的速度分量, fr, fθ 和 fz 分别为垂向、径向和切向的体力分量, σrr, σrθ, σθθ, σzr 和 σzθ为各应力张量的分量, λ 和 μ为拉梅常数。本文中网格划分采用如图1所示的交错网格,应力和速度分量在不同位置离散化, 彼此之间有半个网格间距。半径方向采用等间距网格离散化, 对于每个固定的半径, 横向网格点的数量是相同的,横向网格的间距随着深度的增加而逐渐减小。采用基于交错网格的有限差分与傅里叶伪谱的混合方法, 对波动方程组进行求解。  , f  , f z, f r ,

北京大学学报(自然科学版) f

( r ,  , t ) M ( t )1 r r r 1  ( r rr0 1   Mzr ( t )  r r 1 Mz ( t )  ( r rr0 θ0)为震源中心位置;

M ( t )

Fig. 3

的多重反射波PP和 PPP等。径向和垂向分量上均能看到传播经过火星核的PKP震相, 只能看到很弱的入射到核幔边界发生反射的震相, 如 PCP, PCS和SCP等。切向分量上波形单一, 只能观测到S 波,多重反射波SS 和 SSS最为清晰明显, 同时能观测到 SCS, sscs, SCS2 和 SCS3等核幔边界反射震相。由于模型A中火星壳厚度较小, 因此SCS和 sscs 震相的到时相差不大。在直达S波和SS波之间, 可以看到很强的低速火星壳内部产生的多重反射震相。同时可以看到, 3个分量的理论地震图上面波均比较清晰明显, 径向和垂向分量上的振幅较大, 切向分量上的勒夫面波没有另外两个分量上的瑞利面波强。随着面波的传播, 波列越来越长, 频散现象越来越明显。

图 4展示震源深度为60 km的火星震激发的P波和SV波在全火星模型A中传播的波场快照, 可以看到, 120 s时直达P波和SV波在火星幔上层向下传播, 具有清晰的波前, 同时能看到来自地表的直达P波和SV波的反射波和转换波, pp和 sp震相的波前依次在P波之后向下传播。在火星壳中, 速度比较慢的面波开始产生。在低速的火星壳中, 由于P波和SV波的多重反射和转换作用, 形成直达波之后的较长波列。360 s 时, 直达SV波到达核幔边界, 可以看到比较清晰的P和pp震相的核幔边界反射波, 面波和直达波之后多重反射和转换波的波列

第 56 卷 第 4 期2020 年 7 月

邓迪等 全火星模型中震波传播特征的数值模拟研究

图 8展示SH波在全火星模型B中传播的位移波场快照。270 s 时, 直达SH波到达火星幔上层与下层的分界面。由于模型B的火星壳更厚, 因此向下传播先经壳幔边界反射, 再经火星表面二次反射

第 56 卷 第 4 期

SH 2020 年 7 月ss

照, 可知模型A具有持续时间更长、频散更强的面波, 这与 Toyokuni等[33]的结论相符, 即模型差异严重地影响面波的波列, 较薄的火星壳模型中面波的波列更长。

本文数值模拟的震源时间函数的主周期为10 s,在以后的研究中, 可以考虑不同的主周期, 进一步研究火星壳对不同成分面波传播的影响。模拟过程中采用的震源为二维线源, 与实际的三维点源相比,由于存在波形和几何扩散的差异, 二者的模拟结果有差异, 可以利用二维线源到三维点源结果的转换进行近似的校正[30]。另外, 后续研究中可以尝试对火星震波在三维空间中的传播进行数值模拟, 并且讨论横向非均匀火星壳和地形变化对火星震波传播的影响。

邓迪等 全火星模型中震波传播特征的数值模拟研究

波动方程. 地球科学与环境学报, 2004, 26(1): 61‒64魏星, 王彦宾, 陈晓非. 模拟地震波场的伪谱和高阶有限差分混合方法. 地震学报, 2010, 32(4): 392‒400秦艳芳, 王彦宾. 地震波传播的三维伪谱和高阶有限差分混合方法并行模拟. 地震学报, 2012, 34(2): 147‒156 Wang Y, Takenaka H, Jiang X, et al. Modelling twodimensional global seismic wave propagation in a laterally heterogeneous whole-moon model. Geophysical Journal International, 2013, 192(3): 1271‒1287 Jiang X, Wang Y, Qin Y, et al. Global Sh-wave propagation in a 2D whole Moon model using the parallel hybrid PSM/FDM method. Earthquake Science, 2015, 28(3): 163‒174 Witte D C, Richards P G. The pseudospectral method for simulating wave propagation. Computational Acoustics, 1990, 3: 1‒18 Levander A R. Fourth-order finite-difference P-SV seismograms. Geophysics, 1988, 53(11): 1425‒1436 Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics, 1985, 50(4): 705‒708 Wang Y, Takenaka H, Furumura T. Modelling seismic wave propagation in a two-dimensional cylindrical whole-earth model using the pseudospectral method. Geophysical Journal International, 2001, 145(3): 689‒708 Wang Y, Takenaka H. Sh-wavefield simulation for a laterally heterogeneous whole-earth model using the pseudospectral method. Science China Earth Sciences, 2011, 54(12): 1940‒1947 Herrmann R B. Sh-wave generation by dislocation source — a numerical study. Bulletin of the Seismological Society of America, 1979, 69: 1‒15 Toyokuni G, Ishihara Y, Takenaka H. Preliminary modeling of global seismic wave propagation in the whole Mars // 42nd Lunar and Planetary Science Conference. Woodlands, 2011: Abstract no. 1631 Knapmeyer M, Oberst J, Hauber E, et al. Working models for spatial distribution and level of Mars’ seismicity. Journal of Geophysical Research, 2006, 111: 11006

1. 中国电子科学研究院, 北京 100041; 2. 中国地震局地球物理研究所, 北京 100081; 3. 中国地质科学院地质研究所,北京 100037; 4. 北京大学地球与空间科学学院, 北京 100871; † 通信作者, E-mail: wangwt@cea-igp.ac.cn

摘要 利用地震观测台阵中相邻台站间背景噪声互相关函数随时间的偏移的特性, 考虑参考互相关函数、噪声源变化、反演方法和数据误差的影响, 给出基于背景噪声互相关函数的密集台阵时钟偏差分析及其误差评价方法。对盐源短周期地震观测密集台阵的单台时钟偏差分析结果表明, 基于背景噪声互相关函数的时钟偏差分析方法可以给出连续的单台时钟偏差, 并能够筛选出有明显时钟偏差的台站及偏差出现的时段。该台阵的209个台站中, 有17个台站在观测期内出现大于1 s的时钟偏差, 可能与仪器数据采集装置的硬件或软件故障有关。关键词 地震观测密集台阵; 背景噪声互相关函数; 时钟偏差分析; 噪声源分布

田原等 基于背景噪声互相关函数的地震观测密集台阵时钟偏差分析方法

传感器内部时钟进行授时。在仪器运转正常且GPS信号连接稳定时, 系统的时间精度可以达到50 ns,能够保证数据授时的精度。授时系统在两种情况下可能产生较大的偏差: 1) GPS信号连接失败, 导致数据采集系统只能使用传感器内部时钟, 此时记录的波形可能产生连续的时间漂移, 直至GPS信号重新连接成功[1‒2]; 2) 由于采集或授时系统的硬件或软件出现故障, 数据采集系统的时间可能会出现较

[3‒4]大的变化 。无论哪种情况产生的时间偏差, 数据使用者都无法获得准确的时间校正信息。

在无法获知准确时间信息的情况下, 可以通过记录的波形数据对台站时钟进行简单的分析。利用波形的时钟偏差分析方法主要有两类。1) 通过地震震相在各台站的到时, 估计台站的时钟偏差。例如, Anchieta等[5]利用海底地震仪与临近海岛上的陆地宽频带地震仪的P波到时差, 估计海底地震仪的时钟偏差。此类方法依赖于地震震相的清晰度和出现频率, 因此仅适用于地震频发区域, 对于地震稀少的地区, 不能保证时钟的不间断校正。2) 利用台站间的相干噪声信号, 检测时钟偏差。此类方法可用于检测小孔径台阵、海底台阵及区域台网中的台站时钟偏差。例如, Stehly等[3]利用地球微震频带噪声信号(5~10 s和 10~20 s), 对固定台站的时间偏差进行连续检测; Sens-schönfelder[1]利用两个小孔径火山监测台阵的每日互相关函数相对于全时段叠加互相关函数的偏移量, 给出每个检波器的单台时间日变化, 并指出该台阵的时间漂移主要由GPS授时失败期间的内部时钟漂移产生; Xia等[4]利用几内亚湾的26 s固定噪声源信号, 分析全球范围的台网时钟偏差; 王俊等[6]对江苏省的固定台网进行连续的时钟分析, 并进行时钟偏差的误差分析。

随着台阵技术的进一步发展, 人们开始尝试布设密集的短时观测台阵(台间距数十米至数公里,观测周期十几天至数月), 利用噪声中的高频成分,提取高频(5 s以下)面波乃至体波数据, 以期研究局部小尺度浅层地壳乃至近地表的精细结构[7‒10]。相对于中长周期台阵, 短周期密集台阵具有台间距小、总台站数大和总观测时间短等特点, 数据对时间精度要求较高, 因此需要一种能够同时检测观测时段内所有台站时间精度的快速检测方法。

本文采用类似Sens-schönfelder[1]提出的基于台阵间噪声互相关函数的单台时间偏差分析方法, 针对短周期密集观测台阵的特点, 给出适用于它的台

阵时间偏差分析方法, 并对盐源短周期地震观测密集台阵(简称盐源台阵)进行时间偏差分析, 探讨参考互相关函数计算误差、矩阵优化方案、数据读取误差及噪声源方向改变等因素对分析结果产生的影响, 最终给出盐源台阵所有台站的时钟偏差以及偏差的误差估计。

1 方法原理1.1 基本原理

在理想状态下, 噪声源均匀地分布, 两个台站之间的背景噪声互相关函数代表两台站间的经验格林函数, 与两台站间的路径结构有关, 不随时间变化[11‒12]。以下因素可能会导致互相关函数的变化: 1) 台站间路径结构发生变化; 2) 某个台站的时钟出现偏差, 或数据采集器出现相位偏差; 3) 噪声源分布发生变化。可将某一时段内台站间互相关函数相对于参考互相关函数的时间偏差 δt 用如下式[3]表达:

δt = φ(t)+ D(t)+ ε(t), (1)其中, φ(t)代表路径结构变化导致的时间偏差, D(t)代表由台站时钟偏差或相位误差导致的时间偏差, ε(t)表示由噪声源分布变化产生的时间偏差。短周期密集台阵的观测时间通常小于3个月, 可以认为观测期内台站间的路径结构几乎不随时间变化, 因此相对于δt, φ(t)可以忽略不计。此时, 台阵中任意两台站i与j间的时钟(或相位)偏差可以表示为

Dij(t)=∆i(t)–∆j(t)= δtij – εij(t), (2)其中, ∆i(t)表示第 i个台站的单台时钟偏差, εij(t)表征由于噪声源分布的变化在两台站间产生的时间偏移。若噪声源在空间分布均匀, 或空间分布虽不均匀但随时间变化不大时, 相对于δt, ε(t)的影响可忽略, 此时δt几乎全部由两个台站各自的单台时间偏差贡献。若噪声源的优势方向在一段时间内发生明显的改变, 则ε(t)的影响不可忽略, 需要想办法予以修正。使用较长观测时间(1个月以上)的互相关函数叠加, 可以改善噪声源的分布情况, 进而减少ε(t)的影响, 这种做法多应用于对长时间观测的固定台站进行的时钟偏差分析中[2‒3,6], 但会在一定程度上降低计算结果的时间分辨率。本文利用密集观测台阵的特性, 拟合噪声源变化在不同方向的台站间产生的时间变化, 并加以去除, 尽可能地减小噪声源分布变化的影响, 同时保证计算结果的时间分

639

第 56 卷 第 4 期2020 年 7 月

辨率。

对于整个台阵中的共M个台站, 通过台阵中所有K个台站对互相关函数的时间偏移量(式(2)), 构建求解单台时钟偏差的线性方程组:

引入绝对时间偏移量后,

1.5噪声源分布改变的影响修正

ij

利用背景噪声互相关函数进行时钟偏差分析的一个必要条件是背景噪声源均匀分布, 或在观测期内的分布情况没有明显变化。若噪声源分布不均匀且分布方式发生改变, 噪声互相关函数会在时间轴上产生小幅的移动。对于较低频带内距离较近的两个台站, 噪声源分布变化的影响类似传播方向改变的平面波在台站处产生的到时变化, 此时台站对互相关函数相对其参考互相关函数的时间变化量与噪声源的分布以及台站对连线的方位角有关。对密集台阵而言, 台阵内存在大量方位角及台间距相似的台站对, 这些相似的台站对对噪声源改变应该有相似的响应, 可以利用这一特性来估算噪声源改变在台站间互相关函数上产生的时间响应εij(t)。在盐源台阵每个方位角相似的台站对组 I 中, 我们将每天组内台站对相对时间偏移量 的平均值 作为该

ij组内所有台站对的噪声源时间响应 。此时, 式(5)可改写为

1.6

ij i

ij j

ij

式(2)改写为ref ij

ij ij

ef

I

单台时钟偏差计算及误差分析

(6)进行参考互相关函数修正和噪声源响应修正后, 得到各个台站对的时钟偏差Dij ( t), 继而可以利用式(3)的线性方程组来计算台阵每一天的单台时

n钟偏差 。由于式(3)中的系数矩阵永远为不满秩

i矩阵, 因此在求解线性反演问题前, 要先对矩阵进行优化。本文选择的优化方案为: 选出一个台站作为参考台站, 设定其单台时钟偏差为零, 此时线性方程组变为典型的超定方程组, 可以利用最小二乘法[14]求出稳定解。选取的参考台站k应满足以下条件: 1) 台站在观测期内数据没有间断; 2) 与之相关的所有台站对的参考互相关函数没有明显偏离“正确时间”, 即与参考台站相关的参考互相关函数的

ref极值偏移量dtkj 均小于整个台阵dtref标准差的一半; 3) 台站的时钟变化最小, 即与参考台站相关的所有台站对时钟偏差Dij(t)的标准差最小。单台时钟偏差分析结果的误差由两部分组成: 1) 线性方程组求解过程中引入的误差; 2) 数据误差造成的结果误差, 即台站对时钟偏差Dij(t)的误差对单台时钟偏差的影响。在线性方程组的求解过程作为本文应用示例的短周期密集台阵数据来自由中国地质科学院和中国地震局地球物理勘探中心布设于中国西南部盐源盆地及其周边山区的短周期密集台阵, 台站总数为209 个, 以约2 km的台间距密集地布设于约 40 km×60 km的菱形区域内(图1)。台阵使用的仪器全部为EPS-2型短周期三分量地震仪, 频带宽度为 0.2~150 Hz, 采样频率为 200 Hz。台阵观测时间为 2017 年 6 月 4日至 2017 年 7月 10 日, 单台最长观测时间为36 天, 最短18天。

[15]本文采用 Bensen 等 的背景噪声处理方法,对背景噪声记录数据进行预处理。首先, 将台站原始记录的采样频率降至50 Hz, 并对每天记录完整度大于80%的数据进行去均值和去趋势处理, 再进行去除仪器响应的操作。然后, 对每个长度为1天

641

北京大学学报(自然科学版) 20

图 3模型 A中火星表面的位移理论地震图Synthetic displacement seismograms at Mars’ surface for Model A