基于机器学习和短周期密集台阵资料研究北流地震余震特征

文玺翔, 沈旭章*, 周启明

1 中山大学地球科学与工程学院, 广东省地球动力作用与地质灾害重点实验室, 广州 510275 2 南方海洋科学与工程广东省实验室(珠海), 广东珠海 519082

0 引言

2019年10月12日22时55分,广西壮族自治区北流市与广东省化州市交界发生MS5.2地震,震中位于22.16°N、110.52°E,震源深度10 km,震中烈度VI度. 此次地震是广西地区自2016年苍梧地震以来又一次5级以上的地震,致使震源区5个乡镇受灾,经济损失超过500万元.

震中地区自新生代以来处在喜马拉雅构造域与滨太平洋构造域的复合部位(如图1所示),新构造运动主要包括间歇性整体抬升为主的运动以及沿NE向断裂的断块差异运动(李冰溯等, 2019).李细光等(2007)和张培震等(2013)结合目前的震源机制解、震区发震构造以及GPS速度场等研究成果,得到该区域构造应力场的主压应力方向为NW-SE向.研究区域主要分布北东向蕉林断裂、北西向石窝断裂及近南北向的新丰断裂. 根据野外地质资料(潘建雄和黄日恒,1995),新丰断裂为逆断性质,倾向南东; 石窝断裂为左旋正断层,倾向南西,处于北西向巴马—博白断裂带的分支断裂.自2016年苍梧地震后,该地区进入另一个地震活跃期(任镇寰和罗振暖,1998).深入了解该区域发震构造对于华南地区地震活动、抗震设防以及地震预测预报工作都有重要意义.

图1 地震及台站分布图(a) 台站分布及研究区主要断层. 主图中三角形表示短周期密集台阵位置,圆圈表示余震分布,颜色深浅表示地震发生时间; 图中显示了主震和前震的震源机制解 (数据来自于中国地震局地震预测研究所和广东省地震局),并且展示了该区域三条主要断裂: 石窝断裂,新丰断裂,蕉林断裂; (b) 震源区域放大图,其中圆圈的大小对应于地震震级(震级大小分布于-0.4～2.2之间),颜色深浅代表地震发生的时间.Fig.1 Map of earthquakes and stations(a) Map of stations and the major faults in the studying area; Triangles represent the locations of the high dense short-period seismic array, circles represent the distribution of aftershocks, and different colors represent the time of the earthquake; The focal mechanism solutions are also presented (These data came from the Institute of Earthquake Forecasting (IEF) of the China Earthquake Administration (CEA) and the Guangdong Earthquake Agency); The figure shows three major faults in this region, the NW Shiwo fault, the SE Xinfeng fault, and the NE Jiaolin fault; (b) Zoomed window for the region enclosed by the box in (a), where the size of the circle indicates the magnitude of the earthquake (the magnitudes of these events ranged from -0.4 to 2.2) and the different colors represent the time of the earthquake.

目前不同研究机构结果均显示(广西壮族自治区地震局,苏珊; 中国地震局地震预测研究所; 中国地震局地球物理研究所; 中国地震台网中心,赵博等; 防灾科技学院,Seismology小组): 此次北流地震震源机制解为走滑性质,节面I走向为NWW,倾角约为70°,节面II走向NNE,倾角近乎垂直,主压应力P轴方向为NW向,同该地区的构造应力场方向一致(许忠淮等,1989; 张培震等,2013).尽管根据目前已有的震源机制解结果与野外地质勘探调查、地震烈度分布及余震空间展布可以基本确定NWW向节面与真实的发震断层石窝断裂相对应(李冰溯等,2019; 王小娜等,2020; 阎春恒等,2019).但是受观测资料所限,目前对北流地震的认识还面临以下问题: (1) 震源机制解得出该节面运动性质为右旋走滑,与野外地质调查得到的石窝断裂为左旋走滑的运动性质不同; (2) 根据定位和震源机制解结果,5.2级主震前2 s发生了4.2级前震,震中位于22.16N°,110.53E°,震源深度9.5 km,约在主震北西1 km处.这两次地震的震源机制解虽均为走滑型,节面走向却不同,因此可能是由于不同断裂作用而导致的结果(王小娜等,2020; 周斌,2019),但是缺乏足够观测资料证明; (3) 此外,后续得到的部分余震震源机制解结果与主震前震都不同(郭培兰等,2019).

大量余震事件空间分布和震源机制解是深入认识和了解发震构造的最重要信息.近年来,机器学习方法在识别微震事件中取得了较大突破,此类方法已经被应用于多个区域并取得了较有意义的丰硕成果.如Ross等(2018)使用卷积神经网络基于南加利福尼亚区域的1820000条地震波形所训练的模型,通过交叉验证,在P波到时拾取上与人工拾取平均误差在0.023 s,并且在测定初动极性方面,与人工测定精度相比有95%的准确率; 赵明等(2019)利用U型卷积神经网络在处理首都圈地震台网以及四川地震台网的数据时,Pg、Sg拾取到时上表现出了相较于传统震相自动识别方法更优的均方根误差,准确率和查全率; Wang等(2019)提出基于深度卷积神经网络(DCNNs)VGG-16结构的PickNet算法,在处理日本地区300个地震事件并与日本气象厅(JMA)提供的震相目录进行对比后,表现出了速度快、精度高以及P、S震相拾取数量更多的特点; Zhou等(2019)提出的基于卷积神经网络的DetNet与循环神经网络的PpkNet算法,在处理汶川余震的数据上表现出了地震信号检测的完备性以及震相拾取的稳定性与精确性,且在S波到时拾取上较好的表现;Mousavi等(2020)提出的EQTransfomer方法,在处理日本2016年鸟取地震的波形数据后,得到了相比人工识别震级分布更广,小震数量更多的地震目录,且P波和S波到时拾取精度相对于人工拾取误差在0.05 s左右.Zhou等(2021b)、Jiang等(2021)、廖诗荣等(2021)和苏金波等(2021)将机器学习方法应用于2021年云南漾濞6.4地震研究,获得了高精度、高完备性的余震序列目录.Wong等(2021)使用机器学习方法和波形互相关方法对四川盆地地震进行研究,构建了一个高精度地震目录, 进一步完善了该目录的震级的完整性.Zhou等(2021a)将基于深度学习的PhaseNet方法应用于四川威远页岩气区域临时台网记录的连续地震数据中,获得了相比于中国地震台网约60倍的地震事件,证明了机器学习震相拾取方法在密集台网中的适用性.

本研究将基于2019年10月12日广西北流5.2级地震震源区震后约30天的短周期密集台阵观测资料,使用EQTransformer模型(Mousavi et al.,2020),对余震进行识别及拾取,进而使用Hypoinverse和Hypodd方法对结果进行精定位,挑取信噪比较高的地震事件进行震源机制解反演.根据余震空间分布及震源机制解特征,对该区域中强地震发震构造进行探讨.

图2 地震目录中4次地震事件波形发震时刻标注在每幅图的上方,台站号位于每幅图左侧,按震中距排列.Fig.2 Waveforms of 4 earthquake events in regional earthquake catalogThe origin time is marked at the top of each subgraph, and the station numbers are located at the left, arranged according to epicentral distance.

图3 959台站24小时连续波形检测结果黄色星表示识别出来的地震事件,红色星表示多个台站共同识别到的地震事件.Fig.3 24 hours continuous records and detection results of earthquakesYellow stars indicate detected events, and red stars indicate events recorded by multi-stations.

图4 不同台站同一地震事件的波形(a)、(b) 为台站959(SP424)和台站983(SP327)记录到的同一个地震事件; (c)、(d) 为台站1025(SP209)和台站1011(SP309)记录到的同一个地震事件.Fig.4 Waveforms recorded by different stations of the same earthquake(a) and (b) indicate the same seismic events recorded by stations 959 (SP424) and 983 (SP327); (c) and (d) indicate the same seismic event recorded by stations 1025 (SP209) and 1011 (SP309).

图5 单台站余震识别结果三角形表示台站位置,不同颜色表示识别出的余震事件数目.(a) 每个台站单独识别出的微震数目; (b) 精定位之后每个台站的识别结果; (c) 精定位与地震台网共有的11个余震事件的台站识别结果.Fig.5 The detection results of single stationTriangles indicate station locations, and different colors indicate the number of aftershocks. (a) the number of events detected individually from each station; (b) the number of events detected by each station after relocation; (c) single station about the 11 events in the regional catalog during the same period.

1 数据

北流地震之后一周,中山大学地球科学与工程学院在震源区布设了由120个EPS-M6Q短周期地震仪组成的密集地震台阵,台间距3～5 km,自2019年10月18日起进行了为期约一个月的观测,在观测期间记录到了大量震级较小的余震.在密集台阵观测同期,区域地震目录共记录到了余震13次,这些地震事件在本研究布设的密集地震台阵上都有质量较好的地震事件波形记录.如图2为已有地震目录中4次震级较小的地震事件波形记录,在波形记录中可以看到清晰的P波和、S波.由于本研究中台阵观测密度远远高于区域地震台网,因此也记录到了非常多的区域地震台网没有记录到的余震事件.使用机器学习方法,可以识别出更多震级较小的余震.

2 机器学习方法及微震识别

本文使用最新深度学习神经网络算法(Mousavi et al.,2020)对北流的短周期地震数据进行处理.该方法包含由多任务结构组成的一个主编码器和三个独立的解码器,解码器和编码器主要由一维卷积层(1D convolutions)、双向和单向长短期记忆门(long-short-term memories,LSTM)、NIN网络层(Network-in-Network)、残差连接层(residual connections)、前馈层(feed-forward layers)、转换层(transformer)和自注意力层(self-attentive)组成.编码器在时域中分解地震信号,并通过在它们的时间依赖性上提取高维特征.然后,解码器利用这些信息将高级特征映射到与地震信号相关的三个概率序列: 即每个时间点是否存在地震信号、P波和S波.

相比其他方法,该方法主要优点为: (1) 效率更高且降低了人为主观性的影响; (2) 所受信噪比与检测区间的影响更小; (3) 对于内存和计算的需求更低; (4) 具有更好的泛化性.

图3展示了利用连续记录波形确定余震事件的过程.首先对输入数据进行基于三分量数据标准差的振幅归一化,并将地震事件、P波和S波的识别阈值分别设置为0.4、0.5和0.5,这样做是在保证每个台站识别及拾取结果相对可靠的情况下,尽可能增加关联地震事件的数目.本研究所使用的震相关联算法基于每个事件(包括其P波和S波到时)的识别时间,将不同台站之间识别误差不超过15 s的事件归为同一事件(Mousavi et al.,2020).在SP424台站(仪器编号959)连续24小时的观测记录中,使用机器学习方法识别出了8个疑似余震事件.进而根据其他台站的记录结果对这些疑似事件进行了核实和确定,如SP327台站(仪器编号983)同样在20时34分识别出一例事件,当某一地震事件由至少三个(如3个P或者两个P加1个S,等等)台站记录到时才会将其关联到初始地震目录中.图4为短周期密集台阵不同台站确定的余震事件波形.

使用该方法,本研究最终识别出了地震数据中记录到的大量余震事件波形,并对结果进行关联,得到有441个地震事件的地震目录(如图5所示),大约是同时间段的广东和广西联合测定得到地震事件数目的34倍.

3 地震精定位及震源机制解反演

机器学习方法识别确定出的地震事件形成了关联后的震相文件.为了对余震事件空间分布进行深入分析,我们使用Hypoinverse方法(Klein,2002)对地震事件先进行初定位.初定位使用的地壳模型如图6所示.

图6 初定位使用的地壳速度模型橙色线表示S波波速,蓝色线表示P波波速.Fig.6 Crustal velocity model used for initial locatingThe orange line shows the S-wave velocity, and the blue line shows the P-wave velocity.

为消除所选速度模型对非震源区附近的影响,确认关联后识别地震的可靠性,本研究进一步使用双差定位法(Waldhauser,2001)对余震序列进行重定位.其原理是将同一台站记录到的两个相邻地震(事件对)的观测走时差与理论走时差的残差(即“双差”) 表示为一个方程,利用最小二乘法求解所有台站和事件对组成的方程组,来获得地震事件的相对位置.该方法近年来在国内外地震序列精定位研究方面得到了较好的应用,如房立华等(2018),Yang等(2022).经过上述的精定位过程,最终得到了一个包含299个事件的地震目录.

本研究中使用台站比较密集,对很多余震事件观测的方位覆盖较好,为研究震源机制解提供了可靠保证.因此我们使用gCAP(general Cut And Paste)全矩张量反演方法求解了余震震源机制解及震源矩心深度.由于反演过程对波形记录质量和方位覆盖要求较高,为得到较为可靠稳定的震源机制解,本研究选取发生在主震附近且信噪比较高的65个余震波形进行震源机制解反演.该方法通过对P波震相及其后续震相(Pnl)和S波(或面波)赋予不同的权重,依据理论与实际波形之间的拟合误差函数,使用网格搜索法来反演地震矩张量的六个独立分量,即走向、倾角、滑动角、标量地震矩M0, 以及ISO(各向同性)参数ζ和CLVD(补偿线性偶极子)参数(Zhu and Ben-Zion,2013).

求解处理过程主要包括以下步骤: (1) 在SAC文件中写入头段信息,截取重定位后的地震事件,选用起震时刻前20 s至后100 s的数据; (2) 对截取波形得到的地震数据进行去均值、去尖以及去除仪器响应等操作,接着对数据进行重采样,并进行分量旋转; (3) 采用频率-波数法(FK)来计算(Zhu and Rivera,2002) 拟合所需要的格林函数模型,使用全球地壳模型CRUST1.0,网格中心位于110.5°E,21.5°N.反演的震源深度范围2～16 km,每0.5 km计算一个模型; (4) 将非双力偶分量ISO和CLVD约束为0,拟合求解地震的最佳双力偶节面解.反演中Pnl和S波截取波形窗长分别为25 s和40 s,二者的滤波范围为0.05～0.1 Hz,走向,倾角和滑动角的搜索步长为5°,每隔0.5 km反演一次震源机制获取最佳矩心深度.

4 结果与讨论

4.1 余震精定位结果

初定位的结果(如图7所示)显示余震事件主要发生在研究区域的三个位置,即主震及前震附近(测线AA′-BB′)、蕉林断裂附近(测线CC′-DD′),以及石窝断裂南端(测线EE′-FF′).主震附近的余震平面有近NS和EW向两个延展方向,深度多在3～11 km之间,震源位置随深度有向东移动的趋势,且大部分余震发生在较前震和主震更浅的位置,震中位置也更接近于前震; 蕉林断裂附近与石窝断裂南侧的余震,震源深度集中在5 km左右,虽然这些余震在已有的固定台站地震目录中并没有记录,且余震数目较少,对于指示孕震断层空间展布的作用较弱,但是却可以证明受北流地震的影响,该地区在至少60 km的范围内有多条次级断层处于活动状态.

重定位的结果如图8所示,重定位后的地震目录共包括299个地震事件,而同期的广东地震台网的余震目录中只有13个地震事件.主震附近余震源深度主要在4～11 km之间,与已有地震目录中的事件深度相近.但与之不同的是,该区域还有许多发生在浅部5 km以内的地震,且余震中分布仅仅表现出沿EW向的条带特征.余震震中位置多集中在前震附近,位于主震西北1～3 km范围内; 其它区域的余震,震源深度集中在3～7 km左右,水平和垂向上没有明显的延展趋势.大部分余震发生在布设台站后的15天以内,但是震中位置与震源深度同时间没有明显相关性.

由主震附近精定位结果的三维空间分布(如图9所示)可以看到微震序列随震源深度有沿NEE向移动的趋势,在8～10 km的范围内,微震分布集中,有指示断层面参数的作用(万永革等,2008),该断层走向约为NWW-SEE,倾向约在SSW-SW之间,倾角约为70°,与目前多个机构得到的主震震源机制解节面产状吻合(广西壮族自治区地震局,苏珊;中国地震局地震预测研究所;中国地震局地球物理研究所; 中国地震台网中心,赵博等;防灾科技学院,Seismology小组),因此确定该断层对应主震的孕震断层,即石窝断裂.此外,还可以由微震分布推测该地区存在一条分布于地下4～12 km的断层,走向约在E-NEE之间,倾角近90°,其产状与前震震源机制解节面相吻合(王小娜等,2020),这条断裂与石窝断裂相切,且二者走向较为接近,因此很有可能是该断层先发生错动,而后触发石窝断裂活动.

图7 初定位结果(a) 余震序列分布图,蓝色圆点表示识别的余震事件,红色圆点表示区域地震目录中的事件; (b)—(g) 测线AA′-FF′ 的余震序列纵剖面,绿色五角星表示主震和前震的位置.Fig.7 Initial locations of the aftershocks(a) Distribution of the aftershocks. The blue circle represents the picked events, the red circle represents the events in the catalog; (b)—(g) The longitudinal profile of aftershocks, and the green star represents the position of the mainshock and foreshock.

图8 重定位结果(a) 余震序列分布图,颜色表示识别的余震事件的时间,蓝色圆点表示区域地震目录中的事件; (b)—(g) 测线AA′-FF′的余震序列纵剖面,绿色五角星表示主震和前震的位置.Fig.8 The results of relocation(a) Distribution of the aftershocks, color shows the elapsed time of the recognized aftershocks relative to the main shock, the blue circle shows the events in the regional earthquake catalog; (b)—(g) The longitudinal profile of aftershocks, and the green star shows the position of the main earthquake and foreshock.

图9 主震附近重定位结果(a) 余震序列分布三维图, (b) 余震分布三视图; 不同颜色的球表示不同的发震时刻,紫色平面表示石窝断裂的空间展布,青色平面表示推测断裂的空间展布.Fig.9 3D view of relocation results near main shock(a) 3D map of aftershock distribution; (b) depth distribution of aftershocks in latitude and longitude profiles; Different colors of balls show different seismogenic moments. The purple plane shows the spatial distribution of the Shiwo fault, while the cyan plane shows the spatial distribution of the presumed fault.

图10 震源机制解结果(a) 震源机制解结果,不同颜色的震源球表示不同的震级范围; (b)—(g) 测线AA′-FF′的余震震源机制解纵剖面.Fig.10 The results of focal mechanism solution(a) The result of focal mechanism solution. Different colors of the focal spheres show different magnitude; (b)—(g) The longitudinal profile of aftershock focal mechanism solution from survey line AA′ to FF′.

图11 主震附近震源机制解结果(a) 震源机制解结果,不同颜色的震源球表示不同的震级范围; (b)—(c) 测线AA′-BB′的余震震源机制解纵剖面.Fig.11 The results of focal mechanism solution near the mainshock(a) The result of focal mechanism solution. Different colors of the focal spheres show different magnitude; (b)—(c) The longitudinal profile of aftershock focal mechanism solution from survey line AA′ to BB′.

4.2 震源机制解反演结果

挑选的65个余震事件分布及其震源机制解如图10所示,主震附近的余震(具体震源机制解结果如图11所示) 矩震级基本在2.0以内,多为走滑型,震源机制解的两个节面走向分别在NNW-NWW之间以及NNE-NEE之间,倾角接近直立,P轴方向基本与该地区主压应力方向保持一致.震级较小的地震震源深度在5～8 km,震级较大的地震震源深度集中在9～11 km.该部分余震震源机制解与前震的震源机制解类似,节面Ⅰ走向为NNW-SSE,节面Ⅱ走向为NEE-SWW,震源深度在9～10 km,位于主震NW约1～3 km,进一步证明了上述推测断层的存在.

蕉林断裂附近的余震,震源机制解多表现为走滑型以及逆冲走滑型,两组节面倾角在40°～60°,节面走向多为NW-NWW,震级多在MW1.3～1.9之间; 石窝断裂南侧地震,震源机制解多为过渡型和走滑型,节面Ⅰ倾角在70°～90°,节面Ⅱ倾角在40°～60°,两组节面走向近SN、EW向,震级较小,都在1.3以下.

由于前震和主震的孕震过程并不相同(王小娜等,2020; 周斌,2019),暗示主震的产生可能是由于上述另一条断层的错动而触发,且在主震后,两条断层同时处于活动状态,随之又产生了不同震源机制解的余震.此外,蕉林断裂北东端以及石窝断裂北端的孕震断裂同样受到此次北流地震的影响而处于活动状态,产生了诸多震级较小的地震.由于该部分地震数量较少,分布较为随机,难以确定孕震断层的具体位置及产状.

5 结论

本文基于广西北流5.2级地震震后布设在震源区120个短周期地震仪组成的密集台阵资料,使用机器学习方法识别了441个余震事件,并通过Hypoinverse和Hypodd方法进行精定位.从精定位的299个地震事件中选择65个波形清晰、有多台站记录的事件使用gCAP方法进行震源机制解反演.根据主震、前震以及余震精定位及震源机制解结果,结合震区地质构造情况和烈度情况,对广西北流5.2级地震进行发震构造分析,得到如下主要结论:

(1)震源区位于石窝断裂与新丰断裂的交汇处,可能存在相互切割的断裂.主震附近的余震主要发生在其北西约1～3 km的位置,震源深度在4～11 km以内,震源机制解多表现为走滑型,节面走向与前震更接近.结合精定位与震源机制解结果,确定主震区域主要存在两条断裂,即走向NWW-SEE向的石窝断裂,倾角约70°,以及走向NEE-SWW,倾角近90°的一条断裂,该断裂活动性更强,石窝断裂可能是伴随该断裂的影响而发生错动.

(2)蕉林断裂附近同样识别出了较多余震,这些余震震源深度在3～8 km以内,震级相对于主震区域更小,震源机制解多表现为逆冲型与走滑型.表明蕉林断裂在受北流地震的影响后处于活动状态.此外石窝断裂以南也识别出了部分余震,震级约在1.0左右,多为走滑型与过渡型,该区域并没有出露于地表的断裂,但是处在石窝断裂与蕉林断裂的交汇处,同样可能存在次级隐伏断裂.

致谢本文在处理波形数据的过程中使用到了SAC以及Obspy(Beyreuther et al.,2010) 软件,余震识别使用了S.Mostafa Mousavi等提供的模型和EQTransformer软件,地震定位使用了Fred W.Klein提供的HYPOINVERSE和Waldhauser提供的HypoDD程序,作震源机制解使用了美国圣路易斯大学朱露培教授分享的gCAP程序,本文图件绘制主要用到了 GMT(Wessel et al.,2019),在此一并表示感谢.