2020年MW7.7加勒比海地震：反投影成像确定的一次超剪切事件

曹博男, 盖增喜

北京大学地球与空间科学学院地球物理系， 北京 100871

0 引言

2020年1月28日19∶ 10∶ 24 (UTC)，加勒比海发生一次MW7.7走滑型强震，震中位于古巴南部、牙买加西北部海域.这次地震是近一个世纪以来该地区发生的最大地震.尽管地震规模很大，但它发生在远离人口稠密地区的近海区域，并未造成太大的破坏.根据美国地质调查局(United States Geological Survey, USGS)给出的结果，震中位置为19.419° N，78.756° W，震源深度为 14.9 km，位于该区域莫霍面(最深为大约13 km)下方的上地幔顶部(参考Crust1.0模型，Laske et al., 2013).此次地震发生于北美板块与加勒比板块的交界处(如图1所示)，为左旋走滑破裂.在该地区，地质构造背景呈现多样性和复杂性.加勒比板块沿着北美板块的南缘以大约19 mm·a-1的相对速度向东移动(Molnar and Sykes, 1969; Leroy et al., 1996).主震和余震位于Cayman海沟附近的Oriente断层带上，该断层带分布于大洋中脊东侧，为转换断层，断层南侧的大洋海底较深，是典型的拉张盆地构造.根据前人对该区域发生的一系列历史地震的研究，由加勒比板块和北美板块之间的相对运动所积累的几乎所有的应力，都是在相对较少但较强烈的地震中沿着板块边缘的地震带释放出来，而在这些地震中发生于海洋的地震绝大多数为走滑型地震(Perrot et al., 1997).海洋走滑断层往往是世界上最直、结构最简单的断层，其构造有利于超剪切破裂的发生，且超剪切破裂通常发生于走滑断层的直断层段(Wang et al., 2012).利用反投影(BP)方法，可以对本次走滑型地震的破裂过程进行成像研究，确定地震的运动学参数，为研究地震破裂机制提供观测约束，有助于进一步研究北加勒比板块边界区域的构造和运动学机制.

图1 震源区域地质构造概况图中红色震源机制球显示主震震源的位置、震级和断层滑动信息，黑色震源机制球表示该区域1990年以来M6.0以上的历史地震，灰色圆形表示主震发生后三天内M4.0以上的余震，黑色实线表示断层，白色箭头指示板块相对运动方向.Fig.1 Tectonic overview of focal areaThe red beach ball shows the location, magnitude and focal mechanism of the main shock. The black beach balls show historical earthquakes with a magnitude of M6.0 or higher in the region since 1990. The gray circles denote the aftershocks of M4.0 and above within three days after the main shock. The black solid lines represent the faults, and white arrows indicate the relative plate motion.

地震发生后，快速准确地确定地震震源破裂过程，获取地震能量释放的时空分布情况，对于减轻地震灾害、快速救援响应具有重要意义.区别于震源破裂研究领域中另一种应用广泛的有限断层反演方法(Kikuchi and Kanamori, 1991)，反投影方法具有结果相对稳定，不依赖断层几何和破裂速度等先验信息，无需计算格林函数等优点，可以较快地获得断层破裂的一系列运动学参数(Ishii et al., 2005；Zhang and Ge, 2010；李琦等，2019).2005年，Ishii等首次将反投影方法应用于地震破裂研究，成功得到了2004年苏门答腊安达曼大地震的破裂过程(Ishii et al., 2005).此后，反投影分析就成为了一种常用的地震破裂成像方法，广泛地应用于许多大地震的破裂过程成像，如2004年苏门答腊安达曼大地震(Krüger and Ohrnberger, 2005; Ishii et al., 2007)、2008年中国汶川MW7.9地震(Xu et al., 2009; Zhang and Ge, 2010)、2010年智利地震(Koper et al., 2012; Palo et al., 2014)、2011年日本东北MW9.0大地震(Koper et al., 2011; Meng et al., 2011)、2015年尼泊尔大地震(刘志鹏和盖增喜, 2015; Fan and Shearer, 2015)、2016年MW7.8新西兰地震(刘志鹏等，2018)等.

地震发生之后美国地质调查局(USGS)给出了利用全球地震台网宽频带地震数据得到的有限断层反演结果(https:∥earthquake.usgs.gov/earthquakes/eventpage/us60007idc/finite-fault)，结果显示断层面的滑移分布长约220 km，与地震后沿断层的余震分布相吻合，最大滑移为～24 m，分布距震中约80 km.美国地震学联合研究会(Incorporated Research Institutions for Seismology, IRIS)给出了分别利用北美大陆、欧洲和澳大利亚台阵地震数据得到的反投影破裂结果(http:∥ds.iris.edu/spud/backprojection).由于北美台阵的震中距过小，得到的结果容易受到其他震相的影响从而降低分辨率，而欧洲台阵和澳大利亚台阵的台站分布相对不均匀，且数据质量较差，成像结果略粗糙，通过这些结果无法获得相对准确的地震破裂的运动学参数.因此，本文使用震中距和方位角分布更加有利的美国阿拉斯加台阵记录到的宽频带远场P波垂直分量数据，基于一维层状地球结构模型——AK135模型计算射线走时，利用反投影方法分析本次地震的破裂模式、破裂方向、持续时间、延伸长度、破裂速率和破裂能量释放等运动学参数.

1 研究方法

传统波形聚束反投影方法的核心思想是将格林函数简化为一个时间域内的只与从震中到台站的走时相关、与振幅无关的时移函数(Ishii et al., 2005; Yao et al., 2012).将一次地震所释放的能量离散为许多子事件，通过震源区网格化将格点视为可能的备选子事件，将台站记录的地震图根据从格点到台站的走时进行相移，采用一定的时间窗和滑动步长对叠加后的波形进行截取，将截取波形反向传播回震源区域，并根据聚束能量分布进行成像，以获得地震动态破裂的时空特征.

在实践中，利用短时窗多道互相关(VanDecar and Crosson, 1990)将P波初至对齐，以去除台阵中各台站下方的横向不均匀性的影响，并且将最初破裂点设在震源位置.基于对齐后的波形和主事件重定位方法的概念，对于给定的标号为n的备选子事件源(震源位置标号为o)，其叠加震源波形可以表示为

(1)

(2)

挑选出某一时窗中所有格点上聚束能量的极大值所在的位置即为该时窗内的震源位置.滑动时窗确定的一系列震源即破裂子事件的能量辐射源，同样指示破裂前沿，从而可以获得地震的时空破裂过程.

2 数据处理与参数设置

2.1 台站筛选与数据获取

选取美国阿拉斯加台网的远场P波垂直分量地震数据对2020年1月28日MW7.7级加勒比海地震的破裂过程进行反投影成像分析，数据下载自IRIS数据中心(http:∥ds.iris.edu/wilber3/find_stations，最后访问于2020年1月30日).如图2所示，台阵震中距范围为55～75°，方位角范围为330～340°，断层走向与台站方位角之差大约为95～105°，该方位角范围的台站记录的P波震相清晰且振幅较强，从而能够更好地通过反投影得到破裂过程(Zhang and Ge, 2010).在该震中距范围，平均PP-P大约130 s，PcP-P大约30 s.但是PcP震相频率较低、幅度较小，对高频反投影分析影响不大(杜海林等，2009).

图2 震中与所选台站分布黑色三角形表示台站，黑色五角星表示震中.Fig.2 Distribution of epicenter and selected stationsThe black triangles denote the stations and the black star represents the epicenter.

2.2 台阵响应函数

台阵响应函数(ARF)通常考虑在小孔径或中孔径台站阵列上平面波传播的情况，可以将台阵的几何形状对分辨率的影响量化和可视化(Rost and Thomas, 2002, 2009)，直观地给出由于台阵分布的有限性和离散性，所引起的能量在慢度域中的扩散或泄漏(Xu et al., 2009).本文采用修改和简化后的公式计算台阵响应函数：

(3)

式中Δt(θ,φ)j=t(θ,φ)j-t(hypocenter)j,θ和φ分别为经度和纬度，f为频率，t(θ,φ)j为从震源区某点(θ,φ)到第j个台站的走时，而t(hypocenter)j为从震源到第j个台站的走时.台阵响应函数为正实数，取值范围为0～1，在震源处有最大值，理想情况下，它是一个以震源为中心的脉冲函数.利用式(3)分别计算了图2所示台阵的中心频率为0.5 Hz和1 Hz的台阵响应函数，为了在图上表示得更清晰直观，将得到的能量换算为以分贝表示，即

arf(θ,φ,f)=10log10ARF(θ,φ,f)(dB).

(4)

根据式(4)，可以换算50%的最大能量对应于约-3.01 dB，图3中等值线间距约为2 dB.结果显示(图3)，最大能量的确位于震中位置，但在图示区域范围内，还存在震中以外的个别局部极大值，称其为旁瓣效应.同时，观察到，中心频率越高，能量的分布越集中，考察大于-3 dB(约50%最大能量)的等值线范围，在1 Hz(图3a)和0.5 Hz(图3b)台阵响应函数图中分别对应直径约为50 km和约100 km的区域，这意味着台阵对高频信号具有更高的空间分辨率.

图3 (a)1 Hz， (b)0.5 Hz台阵响应函数白色五角星表示震中，色标表示以分贝(dB)为单位的能量强度.Fig.3 (a) 1 Hz, (b) 0.5 Hz array response functionsThe white star denotes the epicenter, the color scale represents the energy intensity in dB.

2.3 初至P波互相关对齐

为了消除地球介质的横向非均匀性对成像位置的影响，需要对地震数据的初至P波通过波形互相关进行对齐(Ishii et al., 2005; Yao et al., 2012; Zhang and Ge, 2010).采用两步互相关方法以减少数据中噪声的影响(Zhang and Ge, 2010)，首先对原始数据进行0.05～1 Hz的带通滤波，提取波形的低频成分，利用波形多道互相关将所有台站数据与选取的参考台站波形的P波初至对齐，对齐时窗为7 s.在低频成分大致对齐的基础上，再对波形进行0.5～2 Hz的带通滤波，对波形的高频成分进行更精确的对齐，基于对齐的效果进行数据筛选，保留互相关系数大于0.6的波形数据.

经过筛选后，保留了160条原始波形，再对其分别进行0.2～1 Hz和0.5～2 Hz的带通滤波，滤波后的波形如图4所示.可以看出，两个频段的波形，在P波初至时刻，均存在清晰明显几乎完全对齐的同相轴.

图4 对齐筛选后(a)0.5～2 Hz波形数据， (b)0.2～1 Hz波形数据Fig.4 (a) 0.5～2 Hz filtered, (b) 0.2～1 Hz filtered waveform data after aligning and selecting

2.4 震源区网格划分与参数设置

根据USGS给出的震中位置和震源深度，将震源区域(18.419—20.419° N，77.756—82.756° W)划分为0.1°×0.1°的网格，共21×51个网格点作为可能的破裂点.由于相对走时差对深度的变化不敏感(Xu et al., 2009)，整个破裂过程在位于震源深度14.9 km的水平面上进行成像.为了得到此次地震的时空破裂过程，分别利用0.5～2 Hz的高频数据(图4a)和0.2～1 Hz的低频数据(图4b)进行反投影成像.反投影成像研究一般采用该频段地震波(Ishii et al., 2005; Zhang and Ge, 2010; Meng et al., 2011;)，是因为波形频率太高，会导致信号相关性较差，而波形频率太低，子事件定位的分辨率较低.采用一个长10 s的滑动时窗，滑动步长为1 s，持续滑动120 s，对数据进行分段截取和叠加成像.考虑到数据量充足且台站分布较均匀，此次研究使用线性叠加的聚束成像方法，从而避免对原始波形的改动.

3 成像结果与讨论

3.1 结果筛选、破裂时间

利用上述参数设置，分别对两个不同频段的波形数据进行反投影分析，得到了一系列能量辐射源的时空分布.假设断层面不会重复破裂，破裂传播是不可逆的，即破裂不会在某时刻向震中方向反向传播.根据这一原则，筛去逆向传播距离较大且分布散乱的能量辐射源.低频段成像结果在大约76 s之后，出现了较大尺度的逆向传播，能量辐射源位置分散且跳动剧烈，可能受到低频PcP震相的干扰，故只截取保留大约前76 s的低频段成像结果；对于高频段成像结果，在筛去个别逆向传播的能量辐射源后，我们确定地震破裂的持续时间大约为98 s.

3.2 破裂方向、破裂模式

进一步得到了波形聚束叠加不同时窗内的能量快照(图5).高频能量(图5a)空间分辨率和定位精度相对更高，第40 s时窗中聚束能量存在两个极值点，根据前面给出的筛选原则确定第二大极值点为此时窗的破裂点.可以看出，破裂从震中开始向偏西方向单向传播，能量辐射源的位置具有较好的连续性，个别能量辐射源位置快速变化，破裂速度存在波动.低频能量(图5b)空间分辨率和定位精度相对较低，同样也可以看到与高频段类似的西向单侧破裂过程.

图5 (a) 0.5～2 Hz高频， (b) 0.2～1 Hz低频能量快照从发震时刻起始，每隔10 s绘制一次叠加能量快照.白色五角星表示震中，白色十字表示该时窗内破裂能量峰值的位置，即能量辐射源的位置.Fig.5 Snapshots of the (a) 0.5～2 Hz high frequency, (b) 0.2～1 Hz low frequency energyPlotted every 10 s from the beginning of the earthquake. White crosses indicate the positions of the local maximum, which are also regarded as the locations of the energy radiators.

综合两个频段所有的能量辐射源，得到破裂过程反投影成像结果(图6).高频段成像结果(图6a)进一步直观地显示，本次地震的破裂从震中处开始，向偏西方向单向近线性传播，与转换断层的走向基本一致.破裂延伸至大洋中脊附近，推测断层很可能完全破裂并受到大洋中脊的阻挡而停止.低频段成像结果(图6b)可以得到同样的破裂方向、破裂模式，并观察到相同区域内能量辐射源沿传播方向跳动的现象.

图6 (a)高频和(b)低频能量辐射源时空分布黄色五角星表示震源位置；黑色实线表示断层；红色圆形代表能量辐射源，其圆心标示能量辐射源即破裂点的位置，圆的大小与归一化能量成正比，颜色深浅参照色标表示破裂的相对发生时刻.Fig.6 Rupture process of (a) high frequency, (b) low frequencyThe yellow star and black solid lines represent the epicenter and the faults, respectively. The circles show the locations of the energy radiator while sizes of the circles are proportional to the normalized energy and color indicates the elapsed time relative to the origin time according to the color bar.

3.3 破裂距离、破裂速度

分析不同频段能量辐射源距震中的距离与时间的关系，进而可以获得破裂速度.对于高频能量辐射源(图7a)，此次地震的破裂规模约为250 km，平均破裂速率约为2.55 km·s-1.考察瞬时破裂速率，发现破裂过程主要包含了六个速度特征显著的阶段，包括三个低速阶段：Ⅰ.0～40 s，Ⅲ.47～54 s，Ⅴ.62～91 s和三个高速阶段：Ⅱ.40～47 s，Ⅳ.54～62 s，Ⅵ.91～98 s.其中第Ⅰ阶段破裂速率约为1.5 km·s-1，第Ⅲ，Ⅴ阶段破裂传播几乎停滞，第Ⅳ阶段破裂速率约为6 km·s-1，第Ⅵ阶段破裂速率约为3 km·s-1.对于低频能量辐射源(图7b)，前76 s的破裂规模约为200 km，平均破裂速率约为2.6 km·s-1，与高频段基本相同.破裂过程主要包括三个速度特征显著的阶段，分别为两个低速阶段：Ⅰ.0～60 s，Ⅲ.68～76 s和一个高速阶段：Ⅱ.60～68 s.其中第Ⅰ，Ⅲ阶段破裂速率约为1.5 km·s-1，第Ⅱ阶段破裂速率约为6 km·s-1.对比图7a和图7b，发现不同频段的不同破裂阶段在时间和空间上基本对应，速度变化也十分相似.值得注意的是，高频第Ⅳ阶段和低频第Ⅱ阶段的破裂速度均达到约6 km·s-1，超过了该区域上地幔顶部的剪切波速度(约4.5 km·s-1，参考Crust1.0模型，Laske et al., 2013)，为超剪切破裂，这可以解释图5，图6中不同频段的反投影结果在破裂中期均出现的跳跃现象.王墩等的研究结果表明浅源走滑型大地震普遍具有接近或超过区域剪切波速度的平均破裂速度，且在开始时表现出相对缓慢的破裂速度，随后加速，并有相当一部分阶段具有超剪切破裂速度(Wang et al., 2016).我们在此次地震破裂速率的分析结果(图7)，尤其是高频段结果(图7a)中，观察到破裂从缓慢到加速直至发生超剪切的过程，该现象支持王墩等的论断.但本次地震破裂发生于莫霍面以下，该深度剪切波速度显著高于普通地壳的剪切波速度，这可能导致本次地震的平均破裂速率没有达到区域剪切波速度，而只在个别阶段发生了超剪切破裂.

图7 (a)高频段结果及(b)低频段结果破裂速率分析横轴表示各能量辐射源的破裂时刻，纵轴表示各能量辐射源距震中的距离，不同斜率的射线分别代表不同速率.Fig.7 Rupture speed analysis of (a) high frequency results and (b) low frequency resultsThe horizontal and vertical axis respectively represents the rupture moment and the distance to the epicenter of each energy radiator. The rays with different slopes represent different speed.

3.4 破裂能量分析

分析破裂过程中辐射能量随时间的变化，可以发现本次地震的高频能量辐射(图8a)的峰值(归一化能量大于0.8)主要存在于三个区间，分别为：ⅰ .5～15 s，ⅱ.21～30 s和ⅲ.66～79 s，低频能量辐射(图8b)的峰值主要存在于区间ⅰ .14～23 s.对比图8a和图8b，发现高频第ⅰ ，ⅱ区间与低频第ⅰ区间的能量辐射峰值基本对应，这是因为高频段分辨率更高，对能量释放的刻画更精细.对比图7和图8，不难发现：低速阶段与高能区间相吻合，说明该阶段内破裂传播缓慢，破裂能量得以在局部区域集中释放；高速阶段也恰好对应于能量较弱的区间，在该阶段内由于发生了高速且相对匀速的超剪切破裂，导致其能量强度较弱.图8中破裂能量的谷值所在区间，对应于图6，图7中破裂在时空上出现的间断部分，由于我们根据筛选原则删去了这一部分的能量辐射源，导致对此阶段的破裂无法进行有效观测.

图8 (a)高频及(b)低频破裂能量变化横轴表示各能量辐射源的破裂时刻，纵轴表示各能量辐射源的归一化能量，经平滑得到的曲线刻画了的破裂能量的释放过程.Fig.8 Rupture energy change of (a) high frequency and (b) low frequency resultsThe horizontal and vertical axis respectively represents the rupture moment and the normalized energy of each energy radiator.

3.5 结果检验与破裂速度校正

为了确认反投影结果在时间和空间的准确性，分别选取了高频结果中三个能量峰值和超剪切破裂起始、终止时刻对应的能量辐射源，根据反投影得到的破裂位置和破裂时刻，计算该破裂点在波形上对应的理论时间.我们认为能量辐射源对应的波形时间应该与波形振幅较大的部分吻合.如图9a所示，三个能量峰值对应的能量辐射源的波形时间恰好位于波峰处，吻合程度良好，说明通过反投影成像得到的能量辐射源的破裂时刻和破裂位置是准确可信的.由于反投影方法采用一个滑动窗确定破裂时间，超剪切破裂发生前后对应的两个能量辐射源的波形时间与波峰或波谷没有完全对应，因此在该波形时间为中心时刻的时窗(图9a，9b，9c，窗长为10 s)内，重新叠加波形(图9d，9e)并将叠加波形能量极值对应的时刻作为校正后的能量辐射源波形时间，将二者的差值作为时间校正量，修正之前得到的破裂时刻.经过校正后，高频结果得到的超剪切破裂阶段为53～62 s，对应的破裂速度仍然约为6 km·s-1，表明此时段具有较为准确和稳定的破裂速度.

图9 (a)按震中距排列的用于反投影成像的波形数据；对齐后的由超剪切破裂起始(b)、终止(c)时刻所在时窗截取的波形；叠加后的由超剪切破裂起始(d)、终止(e)时刻所在时窗截取的波形(a)中黑色点标示P波初至时刻，蓝色、绿色和青色点分别标示三个高频能量峰值对应的辐射源的理论波形时间，红色、粉色点和波段分别标示超剪切破裂起始、终止时刻对应的理论波形时间和以该时刻为时窗中心所截取的波形.(d，e)中红色和青色点分别标示叠加波形的最大值和最小值.Fig.9 (a) Waveform data arranged by epicentral distance for backprojection imaging; the waveform intercepted by the window at the beginning (b) and ending (c) of the supershear rupture after alignment; the stacked waveform intercepted by the window at the time of starting (d) and ending (e) of the supershear rupture(a) The black dots denote the initial P wave; blue, green, and cyan dots, respectively, mark three peaks corresponds to the theory waveform time of high frequency energy radiator; red, pink dots, and the wave fragments respectively label the theory waveform time of super shear rupture initiation and termination and the waveform intercepted by the window which center is the time. (d)(e) The red and cyan dots denote the maximum and minimum value of the stacked waveform respectively.

4 结论

本文基于AK135地球速度模型，使用美国阿拉斯加台网的远场P波垂直分量宽频带地震数据，对2020 年1月28日MW7.7加勒比海地震的破裂过程进行反投影成像分析.结果显示，该地震沿西向单侧近线性破裂，破裂规模约为250 km，破裂时长约为98 s，平均破裂速率约为2.55 km·s-1.瞬时破裂速率存在明显波动，在大约53～62 s达到约6 km·s-1，符合超剪切破裂的特征.地震在5～15 s，21～30 s和66～79 s存在三个高频能量释放主要峰值，相应时段破裂速率较低，破裂传播甚至近乎停滞，使破裂能量得以在局部区域集中释放.低频段反投影结果在大约76 s后可能受到其他震相干扰，通过保留的部分结果得出的破裂方向、破裂模式、能量峰值和速度区间与高频段结果基本一致.本文的反投影结果位于沿板块边缘的Oriente断层带上并沿走向分布，与主震后余震的空间分布情况和USGS给出的有限断层反演结果相吻合.

Oriente断层的走滑速率约19 mm·a-1(Molnar and Sykes, 1969; Leroy et al., 1996)，且该断层近一个世纪以来未发生过7级以上强震(Perrot et al., 1997)，长期积累的形变蕴含了强大的能量，为形成超剪切破裂提供了动力支持.此次地震的高频能量辐射源延伸至断层末端，受大洋中脊阻挡而终止，转换断层从震中处开始几乎完全破裂，使加勒比板块和北美板块之间由于相对运动所长期积累的应力得以沿着该海洋走滑断层的直断层快速释放.发震断层特殊的空间几何特征和长期积累的应变能，是能够形成此次超剪切破裂的基础.

致谢本文采用的地震信息来自USGS官方网站，地震波数据来自IRIS数据中心，部分图件利用Generic Mapping Tools(GMT)软件绘制．衷心感谢审稿专家及编辑在审稿过程中提出的建设性意见和修改建议.