基于多点源的W-phase反演及其在2004年苏门答腊MW9.1大地震中的应用

彭宇洋, 王 墩,2

(1. 中国地质大学(武汉)地球科学学院, 湖北 武汉 430074;2. 中国地质大学湖北巴东地质灾害国家野外科学观测研究站, 湖北 武汉 430074)

0 引言

一般而言,在大震发生后,因数据及相关参数选择的不确定性,地震震源参数快速测定过程中往往力求简单、准确测定地震基本参数,来服务地震减灾及快速响应。因此,在震后地震参数反演中往往将地震震源当作一个点源。但是在这种情况下,通常只能得到基本的地震破裂参数(如震级、震中/质心位置)。如果地震是“多段式”破裂,那么不同破裂阶段的地震信息将在单个点源设定的反演结果中不能区分。这虽然对震级大小测定没有什么太大影响,但可能使得有些信息被掩盖。如2008年汶川大地震,该地震面波震级8.0,破裂尺度超过300 km,在地震发生初期,地震应急救援力量主要聚焦于震中位置的汶川附近,忽视了其他地区。震后地震调查发现,同在一个断层面的北川破坏同样严重。究其原因主要是在地震应急中往往把地震设定为点源,默认震中破坏最大。但实际上地震是沿线或者沿面破裂的,最开始破裂的地方不一定是地震灾害最严重的地方。

现在发展起来的地震波形拟合方法可以获得精细的断层面位错分布,但其往往需要更多的人工设定参数。目前较为成熟的获得断层面滑移分布的方法为有限断层反演方法,但其准确结果往往在地震震后数天至一个月产生,远不能满足大震应急响应的需要。那么如何在时效性和结果精细准确之间实现一个“两全其美”的方法就显得尤为重要。

1 震源破裂反演方法简介

实时地震研究中使用的方法有传统地震学、测地学和最新提出的重力波[1-2]方法。这里主要介绍与本研究相关的地震学方法和地震波观测数据。按使用地震观测台站远近不同,实时地震研究可大致分为基于近场地震的观测数据和全球地震波观测数据两个主要研究范畴。基于近场地震观测数据的实时地震学研究主要包括利用P波和S波的地震预警技术及系统架构等[3-9]或利用近场或区域地震波数据反演矩张量[10-11]或搜索与之最匹配的理论波形从而求得矩张量[12]。基于全球地震波观测数据的研究主要包括W-phase反演、反投影技术、有限断层反演等[13-23]。因数据易获取性等原因,这里我们主要聚焦基于全球地震波观测数据的实时地震学方法。

基于全球地震波观测数据的实时地震学研究方法主要有基于W-phase的矩张量反演、反投影成像、有限断层反演技术和全波形拟合反演技术。这些方法所需数据类型、时效性、优缺点列于表1。因全波形拟合反演费时较多(如GCMT等),常不用于实时地震学研究,故这里略去。

表1 基于全球地震波数据的实时地震学常用研究方法简介

W-phase因其到时快(介于P、S震相之间)、不易限幅、结果相对稳定、反演速度快等特点,在大震震后应急响应及次生灾害(特别是海啸风险评估)中得到了广泛应用[15,24]。该技术在反演中利用滑动时窗和格点搜索来获得矩张量质心的最佳位置和时刻,所以能较为准确地获得基于点源模型的矩张量解,为海啸风险模拟计算提供科学依据[25]。然而,这种基于简单点源模型的假定,在提高结果稳定性的同时,也失去了对复杂地震事件震源破裂时空过程的解析度,使反演结果产生偏差[26]。

反投影方法自2005年Ishii 等[13]和Kruger等[17]在《Nature》杂志上同期发表利用该方法成功成像2004年苏门答腊岛MW9.1大地震震源破裂后,在地震学界掀起了应用热潮,获得了对震源破裂参数及其他特征的进一步认识[27-31]。该方法无需设定断层模型,只用通过简单的波形偏移叠加,就可获得稳定可靠的震源空间迁移分布图。但是该方法不使用格林函数和反演,得到的结果只能反映破裂空间迁移信息,其物理意义尚未明晰。

有限断层反演技术是三种方法中最成熟的一种。其在20世纪70—80年代就开始被用来反演断层面位错及子事件[16,32-33],并逐渐作为地震震源研究的一种标准手段得到广泛应用[18-22,34-35]。通过有限断层反演可以获得矩震级、断层面位错分布、破裂速度等重要震源信息,为海啸预警、震后应急响应和地震物理等研究提供了基础数据。在震后快速(或自动)有限断层反演中,由于借助尺度律和初始震级来确定模型参数,常常存在模型设置不合理导致结果偏差较大或者稳定性较差。因此,本工作在现有的W-phase反演方法基础上[15,22,25],按照有限断层原理发展了多点源W-phase方法。通过余震空间展布获得粗略的断层空间特征,在震源区域设定多点源,用W-phase反演方法反演各点源的震源机制及地震矩张量。并将其应用到2004年苏门答腊岛MW9.1大地震中。其结果与该处俯冲带海沟滑移方向以及历史地震震源机制特征都有较好的一致性。

2 基于多点源的W-phase反演方法

在介绍本文所使用的W-phase多点源方法之前,首先介绍一下传统的多点源方法。Kikuchi等[16,36-37]开发了一种迭代算法来反演复杂破裂过程的多点源解,其算法的基本操作总结如下:

(1) 依据先验信息进行准备操作。对已知的地震断层面布设网格,并计算每个网格到台站的格林函数,并将格林函数与假定的震源时间函数进行卷积。

(2) 通过已知的格林函数与台站波形,对每个网格点进行反演求解,得到误差最小的点源结果。

(3) 从观测值中减去步骤(2)点源结果与格林函数的合成值,然后将剩余波形作为观测数据输入到下一次迭代。重复这个过程,得到N个点源,直到达到可接受的整体波形拟合或没有进一步合理的点源解。

如果将传统方法运用到W-phase反演中,相当于对式(1)的反复迭代运算,每一次得到一组最佳地震矩张量后都从台站中将这个地震矩张量正演形成的波形剔除掉,剩余波形再进行一次反演。这样做有许多优点,比如能够按照能量释放从大到小依次提取相关震源。但是也有不足之处:先提取出来的震源可能为了尽可能地减小误差而包含了后续震源的地震矩张量、越到后续提取出的震源结果可能越来越不可靠、整体拟合的误差与点源解的个数难以平衡等等。

(1)

(2)

与传统方法不同的是,本文中离散的震源点是一次性选好,而不是在网格中搜索,反演结果是一步到位,不是迭代得到。这样做的好处是不仅仅节省了运算时间,也避免了传统迭代法中后续提取出来的震源结果可靠性不佳的影响。

虽然这样做有许多优点,但其缺点是需要一次性配置好所有的震源点。这里我们首先依据余震空间展布信息获得的粗略断层空间特征来等纬度布设多点源。这样就可以通过W-phase反演获得设定多点源的地震矩张量,从而得到其震源机制及相对释放能量大小,为定量判定破裂空间展布特征及破裂尺度提供科学数据。

3 2004年MW9.1苏门答腊大地震结果及分析

3.1 地震基本信息

2004年12月26日苏门答腊MW9.1地震是一次毁灭性的地震,引发了巨大的海啸,其破裂时间持续了8～10 min,破裂长度延伸了1 300～1 500 km,矩震级 (MW) 为9.1。这使它成为21世纪目前最大的地震,仅比1960年的智利和1964年的阿拉斯加地震小,它是现代宽带地震仪和全球定位系统 (GPS) 设备记录的第一次大型逆冲断层破裂。在其主震破裂后的三个小时内,发生了45次大型余震,这些余震分布在整个大型逆冲断层面和西安达曼断层沿线。其中7次余震在2 h内从主震处由南向北连续单发[38]。海底的垂向运动引发了一场巨大的海啸,淹没了苏门答腊岛北部和孟加拉湾周围的沿海地区,造成了20多万人死亡。

该地震发生在大洋岩石圈逆冲到大陆边缘和岛弧下方的位置,大型逆冲断层型地震基本都发生在这里。根据古登堡-里氏震级的频率关系,此处小地震比大地震多得多,很少有超过MW8以上的地震。然而在极少数情况下,部分地震事件接近或超过MW9,这样的大型逆冲断层地震涉及到断层宽度(150～300 km)、长度(500～1 500 km)很大的运动,释放了俯冲带中相对板块运动的大部分应变能,经常产生巨大的海啸,这大大放大了它们的破坏潜力。在2004年之前,人们认为大型逆冲区地震最可能发生在相对年轻的、以高收敛速度俯冲的洋壳区域[39];1960年智利MW9.5地震和1964年阿拉斯加MW9.2地震都发生了这种情况。然而,苏门答腊地震迫使人们重新评估这些理论[40-42],因为它发生在一个具有较老洋壳和相对收敛缓慢的地区,在其破裂带北部体现的尤为明显。

本次地震发生的具体位置在苏门答腊断层的北部边缘,印度洋板块正在俯冲到苏门答腊大陆架下面[43]。震中附近的相对板块运动约为50 mm/a,但不垂直于海沟,沿苏门答腊逆冲断层的印度—澳大利亚板块边界变形带在10°S～5°N之间广泛发育,造成此处板块构造环境格外复杂[44-45]。

3.2 地震震源机制解反演结果

基于前文中所提到的多点源方法,此次地震多点源W-phase反演依据余震空间展布特征(图1)布设多点源,使用了震中距5°～85°范围内全球地震台网(GSN)记录的地震波形(图2)来反演本次地震的震源机制解,为了测试和验证方法的可靠性,本研究分别设定1、2、3、4、5、6个点源。基于6个点源模型的波形拟合结果见图3。

红色五角星:2004年苏门答腊MW9.1大地震震中;黑色震源机制:震后3个月5级以上地震GCMT结果;红色小圆圈:震后3个月4级以上余震分布位置,共计2 456个(数据源自USGS)。图1 余震分布图Fig.1 Distribution of aftershokes

图2 本次地震所选取的台站分布Fig.2 Distribution of stations selected for the study of Sumatra earthquake

多点源的选取按照余震空间展布得到的粗略断层空间特征进行等纬度选择:在断层点为6个时按照2°间隔进行布设;在断层点为5个时按照2.5°间隔进行布设;在断层点为4个时按照2.8°间隔进行布设;在断层点为3个或2个时按照4°间隔进行布设。其时窗的设置见图4,从上至下依次为1点源到6点源的时窗设置。纵坐标为地震矩释放率(即对地震矩Mo求取关于时间的导数),横坐标为时间。

反演结果如图5所示。当设定为单点源时[图5(a)],震源机制为低倾角断层,与哥伦比亚大学GCMT结果吻合性较好。 当为2个点源时[图5(b)],震源机制断层节面走向在震中附近点源时为北西方向,远离震中点源时为北北西,与该处中强余震震源机制走向一致。 3～5点源结果[图5(c)～(e)]展现的特征比2个点源时更加细致,总体没有太大变化,但波形拟合误差更小。6个点源时[图5(f)],其震源机制断层节面走向的偏转与俯冲带海沟滑移方向表现出高度一致性,也与周缘的中强震震源机制一致,证实了该方法的可靠性及准确性。

在传统W-phase反演中,由于使用的是单一点源,反演结果只能得到一个震源机制解,它获得的是发震断层面整体特征。地震震级越大,破裂尺度往往更长,破裂时间更久,此时基于单点源模型反演并不合适,较难获取复杂发震断层面空间展布形态的具体特征;在空间上其能量释放特征也不清晰,震级测定相对多点源往往偏小(表2)。

以左上角台站为例:IC:台网名,QIZ:台站名,LHZ:宽频带垂直向,φ:方位角,Δ:震中距,α:台站南北分量与正北方向夹角。黑色波形是实际观测波形,红色波形为理论预测波形,两个红点之间是W-phase震相时窗。通过计算实际观测波形与理论预测波形的拟合误差均方根值,当设定为单点源时值为0.417 94 mm,当6点源时值为0.330 16 mm图3 台站分布和基于6个点源模型的波形拟合结果Fig.3 Stations distribution and waveform fitting results based on six point source models

图4 多点源W-phase时窗设置图Fig.4 Time windows diagram of multiple point W-phase source

从基于6个点源模型的结果分析表明,本次苏门答腊地震的断层展布从破裂开始(最南边震源点)到破裂结束(最北边震源点)的断层面都是低倾角断层,但是断层节面倾向由东北逐渐向东东南方向转变,这与海沟处断层的滑移方向相吻合(图5)。本次地震最终得到的震级为9.0,能量分布在第2个点达到最大,后面3个点能量释放相对较少,能量主要集中在9°N以南(表2)。这与现有研究的反演结果相吻合(图6、图7)。

红色线条:断层以及俯冲带运动方向;红色震源机制解:多点源W-phase反演结果;浅蓝色的震源机制解:GCMT此次地震结果;浅绿色的震源机制解:所有该处大于6级以上地震并与本次地震成因相同的震源机制解图5 多点源W-phase反演结果图Fig.5 Results of multiple point W-Phase source inversion

4 讨论

对于W-phase多点源方法,虽然其结果展现出与其他方法结果的高度一致性,但是由于其较为依赖人工设定的点源位置与持续时间,尚不能做到完全自动、实时地运行,这是它将来有待进一步发展的地方。如果我们根据震后余震序列展示的主震破裂断层空间展布来设置多点源,那么其时效性将会从传统W-phase方法的震后30 min得到结果延后到震后几个小时甚至几天才能得到结果。

表2 多点源W-phase相关参数的设置及其结果的具体展示Table 2 Setting of relevant parameters of multiple point W-Phase source and specific display of their results

黑色五角星为震源位置;灰色圆圈是GCMT点源结果,颜色变化表示了滑动量的大小[46-51]。结果源自瑞士地震服务 (SED) (http://www.seismo.ethz.ch/),由Martin Mai汇编得到图6 有限断层反演结果图Fig.6 Finite-fault inversion results

但是如果基于快速反投影方法所得到的最大破裂点位置来进行点源位置设置,那么它将和传统的W-phase方法一样迅速,因为快速反投影方法的结果大概震后13 min加上震源破裂持续时间内产出[53],而传统W-phase结果大概在震后30 min产出。这意味着在等待W-phase波列传到所需最远台站之前,我们就能获得该地震的多点源破裂特征,而且本程序在高性能计算服务器环境下运行,运算时间基本在数十秒内完成。多点源W-phase程序将完全可以实现实时自动反演,从而获得大震破裂空间能量释放特征。

W-phase多点源反演方法表明:此方法可以快速测定破裂过程复杂的大震,获得其变化的震源机制,这是其他基于点源模型的传统方法所不能达到的。特别是对于俯冲带地震,快速了解地震震源机制特别重要。因为走滑型地震一般不产生大海啸,而倾滑型地震产生强烈海啸风险,如果某俯冲带附近先发生了一个走滑地震,然后随之触发了海沟垂向分量的运动,那么本算法将能快速判定逆冲事件,为海啸风险评估及海啸预警提供可靠科学数据。

除此以外,多点源W-phase反演的结果也可以更好协助快速有限断层反演工作。在有限断层反演中需要准确的震源机制解(一般由W-phase反演提供)来布设断层面,现在利用多点源W-phase反演能够给多断层有限断层反演提供更加准确的参数,实现更精细的断层面位错反演。

对于本算法,除了与快速反投影方法相结合以外,也可以与空间格点搜索方法相结合:在获得快速反投影得到的破裂点后,在每个破裂点四周布设加密网格点,并对每个格点组合都运算一遍其不拟合度,找出目标范围内的最佳格点组合。这样将会对快速反投影得到的破裂点更加包容,提升本算法的稳健性。

5 总结

基于传统W-phase方法,本研究开发了多点源W-phase反演方法。在通过结合余震分布结果来等纬度设定多点源后能够获得更加精细的断层空间展布及能量释放特征,为深入认识地震物理过程提供了新手段。

将多点源W-phase反演方法应用到2004年苏门答腊MW9.1大地震后,其结果较好地反映出了该地震1 200 km长破裂的分段震源机制变化及能量释放特征。从基于6个点源模型的结果进行分析表明,本次苏门答腊地震从破裂开始到破裂结束,断层面都是低倾角断层,断层面倾向由东北逐渐向东东南方向转变,能量释放在6°N达到最大,主要集中在9°N以南,这与已知的有限断层方法、此处历史地震的结果以及海沟处断层的滑移方向都有高度的一致性。

除此以外,本方法还可以基于快速反投影方法所得到的最大破裂点位置来实时配置程序,这将使本文结果时效性与传统W-phase方法时效性一致,在震后约30 min即可产出可靠的断层空间展布及能量释放特征,使多点源W-phase反演方法有应用于海啸预警及地震应急的潜力。