汶川大地震地震波传播的谱元法数值模拟研究

根据张勇等[5]提供的主震矩张量解, 结合本研究 中所采用的坐标系 , 换算得到复合源数值模拟中的 三个主要子事件的地震矩张量解 ( 单位为 dyne cm), 分别为
⎡ M rr ⎢ M1 = ⎢ M θ r ⎢ Mφr ⎣ M rθ M θθ M φθ M rφ ⎤ ⎥ M θφ ⎥ M φφ ⎥ ⎦
[17,18]
1
谱元法及数值模型
谱元法的基本思想是在传统有限元方法的基础
上, 每一个单元内使用谱方法或者伪谱法 , 选取正交 多项式为基的插值基函数 , 在各个单元上通过采用 高阶插值基函数 , 提高插值精度, 并进而提高数值解 的收敛速度 . 利用谱元法求解偏微分方程的步骤类 似于有限元方法 , 在网格划分后 , 进行单元处理时 , 谱元法是将解近似的表示成正交多项式的有限级数 展开式 , 单元基函数一般是用截断的正交多项式展 开来确定 [27,28]. 谱元法实质上是高阶有限元方法 , 一 般选择阶数为 4~10 的插值多项式表示为基函数[26,29]. 本文中我们利用的谱元法结合拉格朗日插值基函数 及 GLL 积分, 保证插值点与积分点在同一位置, 从 而使得生成的单元质量矩阵为对角分布的特性 , 同 时采用了显式求解方法 , 大大简化了算法实现过程 , 且易于实现并行计算 , 使得数值模拟方法变得快捷 , 计 算 时 间 缩 短 , 计 算 效 率 进 一 步 提 高 . Komatitsch 等 [24,26,29] 给出了波动方程理论及利用谱元法数值模 拟地震波传播的详细理论基础. 我们的数值模型采用 AK135 模型[30,31], 将整个 球体划分为 144 层, 考虑了地球介质的衰减特性、地 球椭率效应及地表地形 . 其中对于介质的衰减特性 处理是基于标准线性体 [24,26], 地球椭率根据 Clairaut 方程而得出[21,32], 地形数据为 ETOPO5[33]. 数值模型 的网格单元是 27 结点六面体单元, 将地心作为一个 单独的六面体单元 , 沿此六面体的六个面分别向地 表延伸建立整个模型, 如图 1 所示. 从地心到地表逐 层划分网格, 并且在 ICB, 670 km 及 Moho 面三处分 别进行了网格加密,在 ICB 边界处网格首次进行加密, 加密后的网格数为加密前的一倍, 670 km 处网格数再 次加倍 , Moho 面处网格数第三次加倍 , 以满足接近 地表区域需要更高分辨率的要求. 图 1 为数值模型. 图 1(a)显示了地心、外地核及 地幔, 整个球体由六个大小相等的块体组成; 图 1(b) 显示了整个数值模型形态以及地球地形和震中位置. 本次研究中, 使用的程序是 Jeroen Tromp 等人提
关于汶川地震破裂过程的研究结果(图 2(a))来更真实 的再现汶川地震震源的时空特性.
N ijk = N (ε i ,η j , ξ k ) = li (ε )l j (η )lk (ξ ) .
震源项为一双力偶力源, 其一般可以表示为[32]:
f = − M ⋅∇δ ( X −Βιβλιοθήκη BaiduX s ) S (t ),
摘要
汶川大地震震中位于震害活跃的龙门山断裂带上, 其强度超过 1976 年唐山大地震.
关键词
汶川大地震 谱元法 数值模拟 全球地震波 并行计算
地震激发的地震波可以用来研究震源破裂过程及地球内部的层次结构. 从弹性波理论基础出 发, 采用 AK135 理论地球模型, 考虑地表地形、 地球介质衰减及地球椭率等特性, 利用谱元法 结合高性能并行计算, 分别对人工点源和复合源所激发的汶川大地震地震波的全球传播过程 进行数值模拟重现. 数值模拟结果显示了地震波在地球表面的传播形态. 通过比较发现, 复合 源数值模拟结果较点源数值模拟结果更能体现汶川地震震源破裂过程的时空特性. 另外, 将 两种震源的数值模拟结果分别与实际观测台站记录波形资料进行拟合对比, 进一步定性地认 识到多个点源组成的复合源的数值模拟结果较单个点源的数值模拟结果与实际观测资料的拟 合程度更好, 证明了汶川大地震破裂过程为多阶段破裂组成的一个复合破裂过程.
.
该方法计算精度高 , 但不能直接处理起伏自由界面 , 容易受到模型结构限制 . 有限元法适宜于模拟任意 地质构造, 可以任意三角形逼近地层界面 , 保证复杂 地层形态模拟的逼真性 [19], 但是低阶有限元法对高 频和短波长信号的模拟效果不好 , 需要提高网格分 辨率来弥补, 计算量大
[20]
. 采用高阶有限元法也会因
395
严珍珍等: 汶川大地震地震波传播的谱元法数值模拟研究
⎡ 1.523 × 1028 ⎢ = ⎢ −0.065 × 1028 ⎢ 28 ⎢ ⎣ 0.327 × 10
⎡ M rr ⎢ M 2 = ⎢ Mθ r ⎢ Mφr ⎣ M rθ M θθ M φθ
② ①
②
① 中国科学院地质与地球物理研究所 , 北京 100029; ② 中国科学院研究生院计算地球动力学实验室 , 北京 100049 * 联系人 , E-mail: hzhang@gucas.ac.cn 收稿日期 : 2008-11-30; 接受日期 : 2009-03-16 国家重点基础研究发展计划 (编号 : 2004CB418406)、国家自然科学基金项目 (批准号 : 40774049, 40474038)、中国科学院网络中心 (编号 : INF105-SCE-02-12)资助
[2~4]
茂汶北迭溪 , 烈度达到Ⅹ度. 中国地震研究专家及地 质灾害研究专家就汶川大地震进行了深刻地探讨和 研究 , 并详细解释了这次地震的破坏性强于唐山地 震的主要原因. 张勇等[5]利用远场台站垂向波形资料 反演得到主震及部分余震的震源机制解 , 解释了地 震破裂时间过程 ; 王卫民等[6]结合了汶川地震区域的 地质构造背景, 基于有限断层震源模型, 利用实际地 震波形记录及 GPS 等观测资料反演得到了地震震源 破裂过程; 张培震等[7]对地震发震断裂的滑动速率、 复发周期及其构造成因进行了研究 , 认为此次大地 震是低滑动速率、 长复发周期和高破坏强度的巨大地 震, 是一种新的地震类型 . 国内外很多专家也对汶川 地震震源机制等特性进行了理论研究 , 并对 GPS 观 测数据进行了分析研究 . 如用层析成像方法反演龙 门山断裂带地下结构变化 , 地震引起的应力变化以
(1)
(2)
其中, M 为地震矩张量, S(t)为震源时间函数, δ (X−Xs) 为狄垃克函数(Dirac Delta)分布, Xs 为震中位置. 谱元法计算中一般是输入 (2)式中的点源地震矩 张量参数, 但考虑到汶川地震的破裂长度 , 我们尝试 分别以复合源和点源作为震源项进行数值模拟. 我们的复合源模型参数依据张勇等人提出的震 源破裂过程模型 , 认为汶川大地震由三个相继发生 的子事件组成. 三个子事件的发震时间、破裂时间函 数和大小均根据张勇等 [5] 的模型参数 ( 图 2(a) 所示 ) 给出 . 复合源需要的地震矩张量解及震源时间函数
[21,24~26]
, 使得数值模拟全球地震波传播变得可
能 . 本文利用谱元法数值模拟方法, 以弹性波理论为
394
中国科学 D 辑: 地球科学
2009 年 第 39 卷 第 4 期
图1
数值模型
▲为震中位置
供的 SPECFEM3D_GLOBE 谱元法程序. 六面体网格 单元数约为 1205600 个, 网格节点数约为 81517845 个, 时间步长为 0.26 s, 地球表面节点间平均距离约 为 0.141°, 需要计算机进程约为 150 个, 模拟地震波 传播时间约为 200 min, 计算时间约为 30 h, 数值模 拟精度能够到达最短周期约为 27 s. 虽然六面体单元在处理复杂几何形状上显得不 如四面体好 , 在谱元法中仍需要采用四边形单元或 是六面体单元, 这样在单元处理时, 三维基函数可以 表示为一维的张量积[34]. 如基于权重为 1 的正交多项 式的三维拉格朗日插值形函数可以进一步表示为
2008 年 5 月 12 日 14 时 28 分, 四川省汶川县 (31.49°N, 104.11°E)发生 Ms8.0 级大地震 . 这是中国 1949 年以来破坏性最强、波及范围最大的一次地震, 也是 1976 年唐山大地震以来, 国内最为严重的一次 地震灾害, 地震的强度远远超过了唐山大地震, 造成 了巨大的人员和经济财产损失 . 汶川特大地震激发 的地震波传播特性为科学研究本次地震提供了丰富 的资料. 中国及世界各国都对这次大地震进行了迅速报 道和研究[1]. 震后研究结果显示汶川大地震发生在龙 门山断裂带上 , 造成了四川盆地西北缘沉降以及断 层以西部分山区的抬升 . 该区域地质构造复杂
( 图 2(b) 所示 ) 参考张勇等研究结果 . 其中 , t1=0.0 s, t2=20.8 s, t3=49.6 s, hd1=5.2 s, hd2=7.8 s, hd3=5.2 s, h1=1.9, h2=6.4, h3=2.4. 我们希望能够利用张勇等[5]
图 2 地震震源时间函数 [5] (a)及本文数值模拟采用的复合 源时间函数 (b)
Runge 现象产生虚假波[21]等偏差. 谱元法可以看作是 广义有限元方法的一种 , 是解偏微分方程的一种有 效的数值方法 , 它结合了有限元法处理边界结构的 灵活性和谱方法的高精度快速收敛等优点成为地震 波数值模拟的重要工具
[22]
. Maday 等
[23]
将拉格朗日插
值函数应用于谱元法中, 并与 Gauss-Lobatto- Legendre(GLL)积分相结合, 从而保证对角质量矩阵的产生, 大大简化了算法实现和存储量 . 在达到相同精度前 提下, 谱元法能够通过采用较稀疏单元划分, 很好的 解决了有限元方法因单元内高阶插值带来的 Runge 现象 , 并且易与大规模高性能并行计算数值模拟相 结合
中国科学 D 辑:地球科学 www.scichina.com
2009 年 第 39 卷 第 4 期: 393 ~ 402
《中国科学》杂志社
SCIENCE IN CHINA PRESS
earth.scichina.com
汶川大地震地震波传播的谱元法数值模拟研究
严珍珍
①②
, 张怀 *, 杨长春 , 石耀霖
基础 , 基于 AK135 理论地球模型 , 对汶川大地震激 发地震波传播进行数值模拟研究 , 并将数值模拟结 果与实际观测资料进行了拟合对比分析.
.
随着地震波理论在天然地震和勘探地震中的应 用 , 数值模拟地震波传播已成为研究复杂地下结构 的有效方法 , 可以用来探究不同介质构造条件下地 震波场传播特性 . 数值模拟地震波场可以利用数值 的方法 , 重现地震波产生和传播过程中的关键性信 息 . 地震波属于机械波方程, 其数值模拟方法一直在 地震波研究中占有重要地位, 主要包括有限差分法、 有限元法、伪谱法及谱元法等. 有限差分方法是最早 被用于地震波传播数值模拟研究中的 [10,11],能够成功 的应用到非对称结构模型中 , 有效的实现了对横向 不均匀性介质地球模型的地震波传播数值模拟 [12]. 该方法算法简单 , 计算速度快 , 占用内存小 , 但只适 应于相对简单的地质模型 [13]. 伪谱方法又称配置点 法 , 基 于 快 速 Fourier 变 换 、 Chebyshev 变 换 及 Legendre 变换, 被广泛的用于数值模拟区域地震波传 播问题 [14~16] 和全球结构的地震波传播问题 [15], 通过 求解二维柱坐标弹性波动方程结合 Fourier 伪谱法实 现了基于整个地球模型的地震波传播数值模拟 , 也 能 够有效 的考 虑了地 球介 质的横 向不 均匀性
393
.
中国地震信息网资料显示, 有地震记载以来, 震中附 近 200 km 范围内发生过 8 次 7 级以上地震, 最大的 是 1933 年四川茂汶北迭溪 7.5 级地震, 极震区为四川
严珍珍等: 汶川大地震地震波传播的谱元法数值模拟研究
及对今后四川盆地地震灾害进行了初步的危险性评 估
[8,9]