2019年6月16日新西兰克马德克群岛M W7.3地震震源机制快速反演与海啸波高数值模拟

史健宇，徐志国，王君成，李宏伟

（1.国家海洋环境预报中心，北京100081；2.自然资源部海啸预警中心，北京100081）

1 引言

2019年6月16日06时55分（北京时，下同），美国地质调查局（United States Geological Survey，USGS）测定新西兰克马德克群岛附近发生MW7.3地震，地震震中位于30.805°S，178.095°W，震源深度46 km，为该地区2019年发生的最大地震[1]。太平洋海啸预警中心（Pacific Tsunami Warning Centre，PTWC）在地震发生6 min后，针对震中周边地区发布了第一份海啸预警信息，预计震中附近可能发生海啸。震中附近实际水位观测数据表明，本次地震在震中引发小规模的局地海啸，震中附近的拉乌尔岛船湾潮位站观测到0.09 m海啸波高。自然资源部海啸预警中心（National Tsunami Warning Center，NTWC）也在地震发生后发布了海啸预警信息，判定此次地震可能会在震中附近引发局地海啸，但不会对我国沿岸造成影响。

本次地震震中位于克马德克群岛附近海域（见图1），处于西南太平洋地区的汤加-克马德克俯冲带。该俯冲带是环太平洋地震带的主要组成部分，从新西兰北岛沿NNE方向一直延伸到萨摩亚岛附近，是太平洋板块向西俯冲至印度-澳大利亚板块下形成的会聚型边界。该边界主要由汤加-克马德克海沟、火山岛弧以及弧后拉伸盆地组成，其中汤加-克马德克海沟深度超过10 000 m，是全球第二深的海沟[2-4]。全球定位系统（Global Positioning System，GPS）观测资料显示，西南太平洋板块沿汤加-克马德克海沟向印度-澳大利亚板块下俯冲，俯冲速率从克马德克海沟南部的60 mm/a向北逐渐增加，在汤加海沟北部最高达240 mm/a，这是全球移动最快的板块之一[5-6]。强烈的板块间相互作用，使得汤加-克马德克地区成为全球俯冲系统中地震和火山活动最活跃的地区之一，也是目前全球地球动力学研究重要的热点区域之一。USGS历史地震资料显示，从1900年起，该地区共发生MW7.0以上的地震超过50次，包括9次MW8.0以上地震，其中最大的是2018年8月19日发生在斐济附近海域的MW8.2地震。在本次地震震中附近500 km范围内，最大地震为1976年1月14日发生的MW8.0地震。历史地震活动分布表明，该地区地震深度分布范围比较广，从浅部地壳到深部700 km左右，均表现出强烈的地震活动性[7-8]。

图1 研究区域历史地震及区域构造背景图

全球历史海啸灾害数据库（National Centers for Environmental Information，NCEI）显示，汤加-克马德克俯冲带曾发生过36次海啸事件，其中最严重的是2009年9月29日发生在萨摩亚附近的MW8.1地震引发的海啸，最大海啸波幅高达22.35 m，超过100人死亡。距离此次震中500 km范围内，历史上曾发生过11次海啸事件，其中影响最大的是1917年11月16日发生的MW7.5地震引发的海啸灾害，最大海啸波幅达到2.5 m[5,9]。

为了深入理解和认识新西兰克马德克群岛地震发震机制和断层特征，进一步解释此次地震产生海啸的原因，本文首先利用W震相方法快速得到此次地震的断层面几何参数、地震矩和矩心位置等震源参数，并根据断层面反演结果和已有区域板块构造背景，识别可能的发震断层节面。然后，采用OKADA理论模型[10]，应用美国康奈尔（Cornell）大学开发的COMCOT（COrnell Multi-grid COupled Tsunami）模型进行海啸数值模拟，通过模拟结果分析海啸波能量分布特征，对于海啸可能到达的地区及波高进行预测。

2 资料与方法

2.1 震源机制反演

本研究采用全球虚拟地震台网长周期三分量地震波形数据资料（台站分布见图2），选取震中距为5°～90°的地震台站的数据自动提取W震相进行震源机制快速反演计算。W震相是指P波与S波之间斜坡状的长周期震相（100～1 000 s），最早在1992年日本记录的尼加拉瓜地震的位移波形中被辨识出，其所携带的能量可解释为P、PP、SP和S等震相的叠加[11-13]。W震相主要在上地幔传播，其群速度明显高于面波，波形振幅小，不容易受到大地震振幅裁剪的影响。因此，相较于传统的面波波形反演方法，W震相方法更适合大地震快速应急响应，为地震抗震救灾和海啸早期预警提供快速的震源参数[14]。目前，该方法已应用于美国国家地震信息中 心（National Earthquake Information Center，NEIC）、PTWC、日本气象厅（Japan Meteorological Agency，JMA）和NTWC等多家国内外相关机构[15-17]。

图2 本研究所使用台站分布（五角星表示震中位置，三角形表示台站位置）

对于W震相矩心矩张量方法反演而言，待定的震源参数由地震矩张量f=[Mrr,Mθθ,Mφφ,Mrr,Mθθ,Mφφ,Mrθ,Mrφ,Mθφ]t和 矩 心 的4个 时 空 坐 标 参 数ηc=[θc,φc,rc,τc]t（θc是矩心余纬度，φc是矩心经度，rc是矩心深度，τc是矩心时间）确定。完整的W震相矩心矩张量方法反演未知震源参数可表示为：

矩心位置和矩心矩张量为最佳点源位置以及矩张量。式中，矩心坐标ηc可以通过使W震相实际观测数据dω和W震相理论波形Sω之间二次失配函数最小的网格搜索方法估算[16]：

同时，为了提高震源机制解反演效率，需要预先构建格林函数库，我们基于简振正模叠加（Normal mode summation）方 法，采 用PREM（Preliminary Reference Earth Model）全球一维速度模型计算格林函数。

2.2 海啸数值模拟

本文利用的COMCOT海啸模式是由美国康奈尔大学开发的一个非常成熟的长波数值模式[18-20]，基于浅水波方程的有限差分模型模拟海啸产生、传播和增水全过程。该模式有以下特点：（1）模块化设计，控制简单；（2）多层网格嵌套，最多能够嵌套4层网格，每层子网格又能嵌套4块区域；（3）大小嵌套网格的比例不固定，模式能自动调整相关参数，以保持精度和效率的平衡；（4）模式的海啸源比较灵活，地震、海底滑坡、入射波或人为设定的源都可以模拟；（5）能够模拟海啸的产生、传播、漫滩和淹没等海啸全过程[21-25]。对1960年智利海啸、2003年阿尔及利亚海啸以及2004年印度洋大海啸等大量历史海啸事件模拟和复现的结果表明，COMCOT模式具有较好的准确性和适用性[18-19]。

在实际的应用过程中，海啸数值模拟结果与实际观测值之间存在着一定偏差，主要受以下因素影响：（1）地震震源参数的不确定性对海啸数值模拟结果的影响；（2）计算区域的近岸海底地形数据精度不够，不能精确反映地形特征，影响数值模拟结果的分辨率；（3）数值模拟过程没有考虑近岸地形和非线性特征对爬高过程的影响；（4）地球自转产生的科氏力、频散效应、非线性项和底摩擦项设置对模拟结果的影响[25]。

3 结果与解释

3.1 震源机制反演结果

表1中为新西兰克马德克MW7.3地震发生后，基于W震相方法自动快速反演震源机制结果，地震发生40 min内，共得到了6次反演结果。由于震中位于台站分布较为稀疏的海域地区，地震发生17 min后，由7个台站的观测波形数据计算得到初始震源机制解，其后随着更多台站接收到地震波形，参与反演的台站数也逐渐增多，反演结果趋于稳定，在地震发生38 min后得到最终结果：节面Ⅰ走向角186°/倾角37°/滑动角71°，节面Ⅱ走向角29°/倾角55°/滑动角104°，矩心位置为30.4°S，177.82°W，矩心深度50.5 km，理论地震图与观测波形见图3。这个结果与震后统计的USGS、全球矩心矩张量解（Global Centroid-Moment-Tensor，GCMT）、德国波茨坦地学研究中心（The German Research Center for Geoscience，GFZ）、法 国 全 球 地 震 台 网（GEOSCOPE）等多家国际地震机构给出的结果（见表2）相比，震源参数结果均较为接近，都揭示本次地震的震源断层以逆冲性质为主，这说明了我们反演结果的可靠性。

图3 部分台站记录到的观测波形（黑色）与W震相获得的矩张量解拟合波形（红色）对比（波形上的红点范围为W震相时间窗；五角星表示震中位置；圆形表示参与矩张量反演的台站；大圆代表对应波形的台站位置）

对于震源机制反演给出的两个共轭节面，单纯依靠波形数据，我们很难判断哪个节面为破裂主节面，即与实际震源破裂面相符。因此，需要借助震中区域的实际地质构造环境、余震分布和地震传播的方向性等多种手段来判定哪个节面为真实的破裂断层节面。

截至2019年6月24日，本次地震共发生余震超过60次，最大余震震级M=6.5。通过余震分布图我们可以看到（见图4），余震深度分布与该地区正在向下俯冲的太平洋板片深度分布基本一致，结合Hayes等[26]的Slab板片俯冲的倾角分布资料，震源附近的板片俯冲倾角约为35°，我们可以推断节面Ⅰ可能为此次地震的主节面，即走向角186°/倾角37°/滑动角71°。我们认为本次地震为发生在俯冲带的逆冲断层，可能是由正在向下俯冲的太平洋板块造成。

3.2 海啸数值模拟

在进行地震海啸的数值模拟过程中，我们假设初始海啸波由断层的突然垂向错动引起，不考虑断层破裂过程，假设海水表面的向上运动和海底位移一致，可由各向同性的弹性半空间中的断层位错公式确定。根据JMA的经验公式[27]，我们可以快速估计断层长度L=53.7 km，宽度W=26.85 km，断层的平均滑动量D=2.5 m，刚性系数μ取3×1010N/m2。根据W震相反演得到的震源机制结果和区域构造特征，我们选取走向角=186°，倾角=37°，滑动角=71°的断层面几何参数构建断层位错模型。通过COMCOT海啸模拟数值模式计算得到海啸传播最大波幅场（见图5），结果显示，最大海啸波幅为0.28 m，在震中附近，而在拉乌尔岛船湾潮位站附近的模拟结果为0.05 m，与实际观测值基本一致。本次地震发生后，在震中附近引发了局地海啸，但是海啸波幅较小，且本次震中所处的海域附近不利于人的常年居住，人为活动较少，因此，本次地震海啸并未产生破坏性影响。

图5 克马德克M W7.3地震产生海啸传播最大波幅场

4 讨论与结论

本文使用全球虚拟地震台网长周期三分量地震波形数据资料，采用基于W震相的快速震源机制反演方法，研究得到2019年6月16日新西兰克马德克群岛附近海域MW7.3地震信息。结果表明，本次地震是一次浅源逆冲型地震事件，矩心位置为30.4°S，177.82°W，矩心深度50.5 km，断层面解结果为节面Ⅰ：走向角186°/倾角37°/滑动角71°，节面Ⅱ：走向角29°/倾角55°/滑动角104°，该结果与国际上其他机构给出的研究结果基本一致。

本次地震是发生在太平洋板块和印度-澳大利亚板块之间的板间地震。由于密度大的大洋板块俯冲至密度小的大陆板块之下，上覆板块和俯冲板块之间的浅部界面由于断层闭锁-摩擦滑动，造成局部应力不断累积并最后释放而导致地震的发生。板块界面通过瞬时滑移诱发地震来释放日益积累的应力和能量。研究表明，世界上大多数震级较大的浅源地震都是以这种方式诱发的板间地震[28]。因此，我们认为本次地震发生的原因是由于密度较大的西南太平洋板块受重力牵引往密度较小的澳大利亚板块下俯冲，两个板块之间的活动因摩擦力而受阻，板片之间被闭锁，闭锁区域内应力与能量逐渐积累，当应力积累到一定程度时，板块间发生相对滑动产生地震，释放积累的能量。

基于震源机制反演得到的震源参数，采用COMCOT海啸数值模式模拟的海啸传播结果显示，本次地震事件在震中附近引发了较小的海啸，该数值模拟结果与震后实际震中附近潮位站的观测值基本一致。由于此次地震震中远离陆地和岛屿沿岸地区，并未给沿岸居民生产生活带来破坏性影响。

值得说明的是，海啸波首波到达时间和最大波高是海啸预警时需要提供的两个重要参数。海啸数值模拟是定量化海啸预警预报的有效手段之一，能够较为准确地模拟海啸传播以及近岸爬坡的全过程，为海啸预警提供参考数据以及必要支撑。在海啸数值模拟计算过程中，W震相测定的地震断层面几何参数、地震矩、矩心位置和矩心深度等震源参数，为海啸数值模拟提供了真实的海啸波生成初始条件，对提高地震海啸预警的时效性与准确性具有重要的现实意义。