ACTA Scientiarum Naturalium Universitatis Pekinensis
Utilizing Back-projection Method Based on 3-D Global Tomography Model to Investigate Mw 7.8 New Zealand South Island Earthquake
LIU Zhipeng, SONG Chao, GE Zengxi†
School of Earth and Space Sciences, Peking University, Beijing 100871; † Corresponding author, E-mail: zge@pku.edu.cn
Abstract Based on a 3-D global velocity structure model, the authors used teleseismic P wave data from Asian and South American array to image the rupture process of 2016 Mw 7.8 New Zealand earthquake via backprojection analysis. The results show that the rupture is a unilateral one with northeast direction, extending to the ocean. The rupture speed is about 1.65 km/s. There are two phases dominated by high frequency power radiation, occurring during 2040 s and 6080 s, respectively. The second phase is the major one, whose distribution of low frequency power radiation is consistent with the centroid location of the event. The high frequency back projection result of the South American data is better correlated with the peak ground acceleration result. According to the comparison and analysis of the Asian and South American results, it could be inferred that in order to obtain more detailed rupture information of high frequency, the data of array deployed in the region with lower 3-D heterogeneity should be adopted in back projection analysis to enhance the coherency of waveforms. Key words 2016 New Zealand earthquake; rupture process; 3-D earth model; back projection method
2016 年 11 月 13 日 11:02 (UTC), 新西兰南岛东北端发生一起 Mw 7.8 级大地震, 造成数十人受伤,两人死亡, 并摧毁当地十余座建筑。地震还引发多起山体滑坡等次生灾害, 造成很大的经济损失。根
据美国地震学联合研究会(Incorporated Research Institutions for Seismology, IRIS)给出的结果, 震中位置为 42.7245°S, 173.0647°E, 震源深度为 22 km。该地震发生的位置处于Puysegu俯冲带与Hikurangi
俯冲带之间的 Alpine 断层系[1]。这两个俯冲带的俯冲方向完全相反, Alpine 断层系中的断层以转换断层为主, 断层性质的多样性导致该地区地质背景结构极为复杂。因此, 有必要对该地震的破裂过程进行深入研究。
在地震学研究中, 理解大地震的破裂过程可以帮助我们研究破裂机制, 并为震后救援提供快速响应。2005 年, 反投影方法被首次使用, 从此广泛地应用于研究大地震的动态破裂过程[27]。之前的反投影方法研究中, 通常采用一维层状地球结构模型来计算射线走时, 例如利用 PREM[8]、IASP91[9]和AK135[10]来计算 P 波的走时。一维模型的反投影方法仅在下列两种情况下适用: 1) 破裂的区域较小; 2) 射线穿过部分的速度结构与一维模型的结构极其相似。因此, 当破裂范围大至成百上千公里或射线路径穿过部分的速度结构与一维模型相差较大时, 这样的假设就会失效。
近年来, 得益于宽频带数字地震台站的布置和反演技术的进步, 地震学家得到越来越精确的全球三维速度结构模型。与传统的一维模型相比, 这些模型可以更精确地计算走时。因此我们或许可以将三维模型结合到反投影分析中来研究大地震, 从而获得比使用一维模型分辨率更高的结果。本文采用一个全球三维速度结构模型—— LLNL-G3DV3 模型[11], 作为参考模型来进行三维射线追踪, 得到走时, 从而实现反投影分析。本文分别使用来自南美洲和亚洲的台站数据, 对 2016 年新西兰地震破裂过程进行反投影成像和分析。
1 基于三维地球模型的反投影方法
与地震勘探中经常使用的逆时偏移方法类似,传统波形聚束反投影方法的核心思想是将格林函数简化为一个时间域内的只与从震中到台站的走时相关、与振幅无关的时移函数[2,12]。为了得到地震源区能量释放的时空过程, 将地震台站记录到的波形反向传播到源区, 并根据聚束能量(normalized power)进行成像。具体过程为,首先将震源区域网格化为格点, 作为备选的子事件源。然后, 为了找到某一时刻子源所在的位置, 将所有台站记录到的波形在时间轴上根据走时进行移动, 序号为 n 的备选子事件源上叠加的波形可以表示成1 sn ( T ) vi (T ), (1) M M i 1 n ,i 其中, M 是台站的总数, vi是归一化后台站 i记录的原始波形, T 是破裂时间, τn,i代表从格点 n 到台站i的绝对走时。实际应用中,经常使用 P波初至开始很短一段时间内的波形来进行多道互相关[13],从而对齐同一台阵的数据。这样做类似于主事件定位法的思想,可以减少三维不均匀性在震中位置的影响,将初始破裂点成像在震源位置,并确定其余子源相对于震源点的位置。最后,对齐后的数据可以写成ui ( T ) vi ( T ), (2) o ,i τo,i 是从震中到台站 i 的绝对走时, 实际应用中可以cal写成基于一个参考速度结构模型的理论走时 和o ,i走时误差项δτo,i,其中误差项代表地球实际结构与参考模型的差异对走时的影响。因此,式(2)可以分解为
cal u ( T ) v ( T )。(3) i i o ,i o ,i当进行实际反投影成像时,误差项 δτo,i可以通过互相关估计。值得注意的是,直接通过互相关得到的时移ξτo,i只是一个相对值。所以有 Co , (4) o ,i o , i这里Co是一个常数, 用来代表 δτo,i 与τo,i的差异。虽然无法得到Co的具体值,但是不会影响反投影的结果。因此,对齐后的波形可以写成cal (5) u ( T ) v ( T )。i i o ,i o ,i基于对齐的波形和主事件定位法的概念,将式(1)变形为
1 M ( )] 。 (6) sn ( T ) ui [T M n ,i o ,i i 1与式(2)同理, 式(6)可以分解为1 s ( T ) u [ T ( cal ) ( cal )] , (7) M n i 1 i n ,i n ,i o ,i o ,i M δτn,i 是从格点 n 到台站 i 的理论走时与绝对走时的差值。如果根据式(7)中各项的物理意义进行重新分组, 则有
1 M sn ( T ) ui (T ), (8) cal n , i o ,i n , i o ,i M i 1其中, cal cal cal 是一阶项(理论走时差别n , i o ,i n ,i o ,i项); 是二阶项(误差差别项)。n , i o , i n , i o ,i虽然式(8)在形式上比式(1)复杂, 但是更容易实现。
在传统反投影方法中, 一阶项 是基于一cal n , i o ,i维各向同性层状地球模型计算的, 二阶项
n , i o ,i经常被忽略, 因此, 震中所在的格点到各个台站的走时误差被用于所有格点。但是, 这样的假设只对震中附近的格点有效。随着可能的子源位置向远离震中的方向移动, 假设就会逐渐失效。尤其是使用不那么精确的模型时, 的误差将会累积至cal n , i o ,i , 对反投影结果会有负面影响。很多借助
n , i o ,i余震信息进行改进的反投影方法被引入地震学研究中[1416]。与传统的反投影方法相比, 在这些研究中地震的破裂过程得到校正, 并展现了更高的空间分辨率。这些依赖于余震的方法将对余震数据进行互相关得到的时移等相关信息用在余震震中所在的格点上, 然后通过不同的插值方法得到所有离散化的源区格点二阶项的值。但是, 这种校正方法有两个局限性: 1) 通过互相关得到的时移并不是绝对值,而是一个相对值, 所以根据式(4), 中会有
n , i o ,i一个位置的常数, 这个常数对每个格点是不同的; 2) 由于余震的定位精度可能较低, 因此会导致插值的结果有偏差, 从而引起二次误差, 影响反投影的结果。以上两点导致这些依赖余震的反投影校正方法在实际应用中并不稳定。使用更多的余震信息来校正反投影结果的初衷是, 可以借助余震信息得到一个更精确的走时表,用于反投影的实现。如果使用余震信息, 这个走时表可以将每一个子源到台站的射线路径上三维不均匀性考虑进去, 从而得到更精确的走时。实际上,三维地球速度结构模型就是一个基于大量的全球地震数据和观测的综合的结果。这让我们想到, 或许可以采用三维地球模型直接计算走时。当然, 这种方法等效于使用前震信息, 而不是某个地震的余震信息。这种方法的优点是, 如果三维模型是更可信的, 那么一阶项 将会更精确, 所以二阶项cal n , i o ,i 会被最小化, 并且不会引入由于使用余震
n , i o ,i信息而导致的误差。除此之外, 由于无需等到余震发生, 因此可以更快地得到结果。本研究使用的三维全球 P 波速度结构模型——
[11] LLNL-G3DV3 模型, 是由 Simmons 等 反演得到的。这个模型基于一个非球面的地球结构, 并以层状镶嵌的格点为框架, 包含地壳和上地幔的多层起
伏界面。为了得到这个模型, 共使用全球 7000 个台站来自 1 万个地震事件的总计约 280 万条的直达P 波记录。所有的地震事件都通过名为 Bayesloc 的多事件定位算法[1718]重新精确定位, 并作为得到模型的输入参数, 这显然提高了地震定位的精度。由于使用了大量的质量较好的数据, 因此我们使用这个模型来进行射线追踪, 从而计算三维模型下 P 波的走时, 进行反投影分析。我们利用一个专门支持提取该模型信息的电脑程序 LLNL-EARTH3D①来提高后续计算过程的效率。
2 数据处理流程
为了使用三维模型反投影方法对2016年 7.8 级新西兰地震的同震破裂过程进行高频成像, 我们从IRIS 的数据中心(http://ds.iris.edu/wilber3/find_stations/ 5197722)获取来自南美洲和亚洲的两个台阵数据,分别进行计算。这些台阵位于震中位置的不同方位(图 1)。首先, 所有的地震记录都经过 0.5~2 Hz 的带通滤波, 并且进行归一化处理。然后, 所有的波形都通过多道互相关, 并且根据理论到时进行对齐(图2)。我们采用的震中位置为 42.73°S, 173.07°E,深度为 22 km。震源区(42.23°—44.73°s, 172.57°— 175.07°E)被分隔为 0.1°×0.1°的格点, 这表示源区被分隔为25×25个格点。由于相对走时差对深度的变化不敏感[19], 整个破裂过程是在震源深度为 22 km的水平面上进行成像的。为了得到此次地震的时空破裂过程, 我们采用一个长为 10s 的滑动窗, 滑动步长为 1 s, 对数据进行分段化处理。考虑到数据量充足且台站分布较好, 我们使用线性叠加的聚束成像方法, 从而避免对原始波形的改动。
3 结果
基于三维模型, 分别得到使用南美洲数据和亚洲数据发震后 90s 的反投影破裂过程(图 3)。选择截取至90 s 的原因是, 90s 之后的聚焦结果较差, 无法观测到稳定聚焦的破裂点。两组数据的结果有共同的特征, 同时也有一些差异。首先, 两组结果都证明, 本次地震的破裂从震中处开始, 向东北方向传播, 这与该地区 Alpine 断层系中断层的走向基本上一致。此外, 在高频能量辐射的过程中, 两组数
据的结果均表明出破裂过程中存在两个能量释放的峰值, 如图 4所示。第一个峰值发生在初始破裂后20~40 s, 第二个峰值发生在初始破裂后60~80 s。从两组结果的差异来看, 南美洲数据的结果破裂尺度较大, 约为140 km, 并且一直传播到海岸线; 亚洲数据的结果破裂尺度较小, 约为100 km。分别将两组结果的高频辐射源位置投影在通过最小二乘法拟合的各自的破裂方向上, 利用得到的投影距离来研究破裂速度(图5)。南美洲数据得到的破裂速度较大, 为 1.65 km/s; 亚洲数据得到的破裂速度较小,为 1.28 km/s。
4 讨论
区别于前人的方法, 本文使用三维地球模型,因此有必要对该三维模型对于本次事件的分辨率进 行检测。检测的方法为, 选取本次地震的一个余震(http://ds.iris.edu/spud/eventplot/13298018), 将其波形视为主震事件后的连续波形, 对其相对于主震震中的矩心位置分别使用一维模型和三维模型进行反投影(图 6)。需要注意的是, 由于余震的绝对位置定位存在误差, 因此目录位置仅作为参考, 而不作为判断反投影位置是否准确的依据。这里使用的判断依据为: 使用不同台阵数据得到的反投影余震位置基于哪个模型更为一致。可以看到, 使用一维模型时, 不同数据成像的位置相差较大。当使用三维模型后, 南美洲和亚洲数据所确定的矩心位置明显靠拢, 展现出更高的分辨率。因此证明, 使用三维模型可以提升该区域的反投影成像精度。但是,由于该余震距震中位置较近, 而远离震中的区域(如图 3中破裂尾端区域)没有适用震级的余震发生,
因此, 我们无法对远离震中区域成像的分辨率进行测试。为了探究两组台阵数据得到的结果在破裂的后半段有差异的原因, 分别使用低频数据(0.05~1 Hz)对 50s 之后的破裂(即第二个能量释放峰值的时间段), 使用 40 s 的时间窗进行聚束能量分析(图 7)。 可以发现, 两组数据低频能量的峰值基本上与该地震的矩心位置一致(图 8)。由于矩心位置也是使用低频数据进行矩张量反演得到的, 且低频部分60s前的能量极低[20], 因此低频反投影方法与矩张量反演方法得到的矩心位置一致, 说明这两组台阵中台站位置的设定均可以对低频数据进行相对准确的成像。一般来讲, 由于数据相关程度较高, 低频数据的能量较为集中, 受三维非均匀性影响较小, 因此比高频数据更容易成像, 但分辨率较低。两组台阵
均能对低频能量位置进行成像, 说明台阵中台站的位置分布情况较好。但是, 对于破裂末期, 南美洲的高频数据反投影仍然可以成像, 而亚洲数据的高频结果则停滞在矩心位置附近。这可能是由于亚洲台站大多位于海岛, 台站下方三维不均匀性较强,虽通过三维模型有一定程度的修正, 但其影响仍然无法消除, 所以对破裂末期的高频部分聚焦效果较差, 其结果中破裂前沿停滞在极大能量释放区域。相反, 南美洲台站位于相对稳定的南美洲大陆, 三维非均匀性相对较弱, 因此通过三维模型修正后,高频部分更容易相干聚焦, 对于破裂末期的成像分辨率也更高。此外, 图 3 中, 南美洲数据的结果与通过地表观测得到的可能破裂的断层分布也更为一致[21]。因此, 由南美洲数据得到的结果能更好地展示本次地震的破裂细节特征。
为了更好地理解高频源分布的意义, 将反投影得到的位置与地面峰值加速度(peak ground acceleration, PGA)进行对比(图 9)。可以发现, 南美洲数据的结果与 PGA 的分布更为一致, 破裂均传播到东北端海中的断层。由于 PGA 是与加速度相关的物理量, 因此能更好地展示高频部分的信息, 所以从理论上讲, PGA与高频反投影结果应该更为耦合。这表明, 我们可以通过 PGA 的结果, 对反投影分析的结果进行一些有效性分析。同时, PGA 与高频反投影结果互相补充, 可以对地震高频特征分析提供帮助。
前面的分析说明, 与亚洲数据相比, 南美洲数据给出的结果能更完整地展示本次地震的高频破裂特征。因此, 我们认为南美洲数据给出的 1.65 km/s的破裂速度更为准确。亚洲数据给出的破裂速度较 小, 是因为末期的破裂停滞在距离震中约100 km处, 破裂前沿未向远处移动。1.65 km/s 的破裂速度
[20]与 Hollinsworth 等 基于澳洲数据给出的 1.5~2.0 km/s 的速度值区间较为一致, 略为大于Zhang 等[22]使用多台阵联合反投影给出的1.4~1.6 km/s 的破裂速度。
与传统一维模型反投影方法相比, 本文使用的三维模型反投影方法的最大改进是对走时的估计更为准确, 并依此得到空间分辨率更高的地震破裂过程反投影结果。为了更直观地展示使用三维模型后对走时估计准确性的提升, 我们在地震图上分别标出使用三维模型和一维模型得到的理论走时(图10)。可以看到, 对于南美洲数据, 一维模型估计的走时明显地系统性偏小, 三维模型的估计更为准确;对于亚洲数据, 一维模型计算的走时比三维模型稍大, 但由于亚洲数据信噪比较低, 很难直观地分辨P 波初至, 因此无法准确地判断哪个模型估计得更准确。值得注意的是, 根据式(4), 由走时预测的系统性偏差所得到的常数并不会影响反投影的结果,因此我们分别计算基于三维模型和一维模型的理论走时与实际走时的差值, 并进行去平均处理以消除系统性偏差, 得到理论走时与实际走时的相对误差(图 11)。可以观察到, 使用三维模型的结果明显更为集中, 走时相对误差值接近 0, 而一维模型的结果涨落较大, 特别是对于亚洲数据, 走时计算结果更为离散。
5 结论
本文分别使用南美洲和亚洲的数据, 用三维地球模型对 2016 年 Mw 7.8 级新西兰地震进行反投
影分析。结果显示, 本次事件为东北方向的单边破裂, 破裂长度约为140 km, 延伸至海中, 破裂速度约为 1.65 km/s。该地震的高频能量释放有两个阶段, 分别为 20~40 s 和 60~80 s, 其中, 第二个阶段为能量释放的主要阶段。该阶段的低频能量聚束分布与该地震的矩心位置较为一致。本文讨论了亚洲数据成像结果不如南美洲数据理想的原因, 认为使 用三维模型可以在一定程度上减弱地球介质横向不均匀性对反投影结果的影响, 但是受限于三维模型的分辨率, 该影响无法完全修正, 特别是在台站分布区域的结构差异较大、区域性非均匀性较强的情况下, 修正幅度有限。这对选择用于反投影分析的数据提出了要求, 即尽可能选取位于稳定的横向不均匀性较弱区域的台阵。此外, 本文对比了使用不
同模型时理论走时与实际走时的差别, 发现使用三维模型可以明显地提升走时估计的精确度。
对余震位置的反投影成像测试说明, 与一维模型相比, 三维模型的使用可以有效地提高空间分辨率。同时, 与使用余震信息对不均匀性进行修正的方法相比, 三维模型反投影无需等待足够数量的用于校正的余震的发生(数天甚至数周), 并且可以避免因余震定位不准而引入的二次误差, 因此可在震后第一时间给出更准确的反投影结果, 达到快速响应的目的, 为震后救灾提供帮助。随着观测数据的不断累积, 三维模型不断的迭代升级, 分辨率也会获得提升。因此, 三维地球模型可以在今后的应用中为反投影分析提供更多帮助。
参考文献
[1] Shi X H, Wang Y, Liu J, et al. How complex is the 309‒311 2016 Mw 7.8 Kaikoura earthquake, South Island, New Zealand?. Science Bulletin, 2017, 62(5): [2] Ishii M, Shearer P M, Houston H, et al. Extent, duration and speed of the 2004 Sumatra-andaman 933‒936 earthquake imaged by the Hi-net array. Nature, 2005, 435: [3] Meng L, Inbal A, Ampuero J P. A window into the Tohoku‐oki complexity of the dynamic rupture of the 2011 Mw 9 earthquake. Geophysical Research Letters, 2011, 38(7): L00G07 [4] Zhang H, Ge Z. Tracking the rupture of the 2008
Wenchuan earthquake by using the relative backprojection 2551‒2560 method. Bulletin of the Seismological Society of America, 2010, 100(5B): [5] Krüger F, Ohrnberger M. Tracking the rupture of the 937‒939 Mw = 9.3 Sumatra earthquake over 1,150 km at teleseismic distance. Nature, 2005, 435: [6] Walker K T, Shearer P M. Illuminating the near‐sonic rupture velocities of the intracontinental Kokoxili Mw 7.8 and Denali fault Mw 7.9 strike‐slip earthquakes with global P wave back projection ima-ging. Journal 1205‒1222 of Geophysical Research Solid Earth, 2009, 114(2): [7] Satriano C, Kiraly E, Bernard P, et al. The 2012 Mw 8.6 Sumatra earthquake: Evidence of westward sequential seismic ruptures associated to the reactivation of a N-S ocean fabric. Geophysical Research Letters, 2012, 39(15): L15302 [8] Dziewonski A M, Anderson D L. preliminary reference 297‒356 earth model. Physics of the Earth & Planetary Interiors, 1981, 25(4): [9] Kennett B L N. Traveltimes for global earthquake 429‒465 location and phase identific. Geophysical Journal International, 1991, 105(2): [10] Kennett B L N, Engdahl E R, Buland R. Constraints 108‒ on seismic velocities in the Earth from traveltimes. Geophysical Journal International, 1995, 122(1): 124 [11] Simmons N A, Myers S C, Johannesson G, et al. LLNL-G3DV3: global P wave tomography model for improved regional and teleseismic travel time prediction. 189‒200 Journal of Geophysical Research Atmospheres, 2012, 117(B10): [12] Yao H, Shearer P M, Gerstoft P. Subevent location and rupture imaging using iterative backprojection for 1152‒1168 the 2011 Tohoku Mw 9.0 earthquake. Geophysical Journal International, 2012, 190(2): [13] Vandecar J, Crosson R. Determination of teleseismic relative phase arrival times using multi-channel crosscorrelation 150‒169 and least squares. Bulletin of the Seismological Society of America, 1990, 80(1): [14] Meng L, Zhang A, Yagi Y. Improving back projection
physics‐based (2): imaging earthquake. 628‒636 approach: with Geophysical a novel a case Research study of Letters, the aftershock 2015 2016, Gorkha calibration 43 [15] Ishii wave Andaman M, imaging Shearer and 28 of P M, the March Houston 26 2005 December H, Sumatra et al. 2004 Teleseismic earthquake Sumatra‐ P ruptures using the Hi‐net array. Journal of Geophysical 429‒430 Research Atmospheres, 2007, 112(B11): [16] Palo M, Tilmann F, Krüger F, et al. High-frequency seismic radiation from Maule earthquake (Mw 8.8, 2010 February 27) inferred from high-resolution back projection 2014, analysis. 199(2): 1058‒1077 Geophysical Journal International, [17] Myers S C, Johannesson G, Hanley W. Incorporation of probabilistic seismic phase labels into a Bayesian multiple-event seismic locator. Geophysical Journal International, 2009, 177(1): 193–204 [18] Myers S C, Johannesson G, Hanley W. A Bayesian hierarchical method for multiple‐event seismic location. Geophysical Journal International, 2007, 171: [19] Xu 1049‒1063 Y, Koper K D, Sufri O, et al. Rupture imaging of the Mw 7.9 12 May 2008 Wenchuan earthquake from back projection of teleseismic P waves. Geochemistry, Geophysics, Geosystems, 2009, 10(4): Q04006 [20] Hollinsworth J, Ye L, Avouac J P. Dynamically triggered slip on a splay fault in the Mw 7.8, 2016 Kaikoura (New Zealand) earthquake. Geophysical Research Letters, 2017, 44(8): 3517–3525 [21] Bradley B A, Razafindrakoto H N, Polak V. Groundmotion observations from the 14 November 2016 Mw 7.8 Kaikoura, New Zealand, earthquake and insights from broadband simulations. 740‒756 Seismological Research Letters, 2017, 88(3): [22] Zhang H, Koper K D, Pankow K, et al. Imaging the 2016 Mw 7.8 Kaikoura, New Zealand earthquake with teleseismic P waves: a cascading rupture across multiple faults. Geophysical Research Letters, 2017, 44: doi: 10.1002/2017GL073461