2017年3月渤海地震序列微震检测与发震构造分析

2021-07-21马婷邓莉王晓山宋程谭毅培

中国地震 2021年2期

马婷邓莉王晓山宋程谭毅培

1)天津市地震局，中国地震局地震工程综合模拟与韧性抗震重点实验室，天津 300201 2)河北省地震局，石家庄 050021

0 引言

地震序列指发震时间和震中位置相对集中的地震丛集。在华北地区，有数字化地震记录的大震相对缺乏，而中小地震波形数据较为丰富，其数量优势能够在地震学研究中发挥重要作用(Brodsky，2019)。对中小地震序列发震构造的分析，成为华北地区地震活动性研究和地震危险性分析的重要基础。

渤海海域地处郯城-庐江断裂带(简称郯庐带)与张家口-渤海地震带(简称张渤带)的交汇部位，新生代主要发育张性和剪切性的断裂，断裂主要分NNE-NE向、NW向和近EW向3组(侯贵廷，2014)，地震地质构造较为复杂。渤海海域历史上强震活动频繁，史料记载1548年以来发生7级以上强震4次(国家地震局震害防御司，1995)。而近十多年华北地区有数字化地震波形记录以来，渤海海域地震活动较弱。2017年3月24日渤海海域发生ML4.4地震，为十年来本区域发生的最大地震。根据河北省测震台网地震目录，3月24—28日，渤海海域发生地震18次，除3月24日ML4.4地震外，3月27日还发生1次ML3.5地震，地震震中位置如图1 所示，ML2.0以上地震的震源参数见表1。研究此次地震序列的发震构造，有利于深入了解现今渤海海域地震活动特点，为区域地震危险性分析提供基础资料。

图 1 渤海地震序列2次ML3.0以上地震震中及本文所用台站分布红色和蓝色的圆圈分别对应ML4.4和ML3.5地震；黑色三角形表示本文使用的台站；断层数据引自邓起东等(2007)

表 1渤海地震序列中ML2.0以上地震的震源参数

地震精定位和震源机制解是分析发震构造的重要手段，能够较为精细地刻画发震构造几何特征(Long et al，2015；Yoon et al，2019)和地震破裂的错动方式(易桂喜等，2019；王小娜等，2019)。前人对渤海地区历史强震、1969年7.4级地震和区域中等地震的震源参数进行了复核，加深了对渤海海域地震活动与区域地质构造关系的认识(环文林等，1989；赵燕来等，1993；谢卓娟等，2008)。然而渤海海域连续运行的固定地震监测台站稀少，区域台网距离本次地震序列震中最近的台站震中距约为80km(图1)，对分析渤海地震活动的研究造成影响。如何利用现有的数字化地震波形，获得可靠的渤海海域地震精定位和震源机制计算结果，值得深入研究。

随着信息技术的发展，微震检测和波形互相关技术能够深入挖掘数字地震波形中携带的信息，逐渐成为地震序列分析中不可或缺的一环。常用的微震识别方法主要有3类，第一类基于地震波形与背景噪声的振幅差别，多利用长时窗与短时窗的振幅或频率等特征的变化(Allen，1978；Earle et al，1994)以及结合赤池信息准则(AIC)进行计算(刘希强等，2009)，此类方法已应用于中国地震台网日常地震自动速报业务中；第二类为模板匹配方法，其利用地震波形相似性(Peng et al，2009)或音频指纹算法(Yoon et al，2015)识别微震，近年来此方法已广泛应用于余震序列(Wu et al，2017；Li et al，2018)、前震序列(Kato et al，2012；Ruiz et al，2014；Warren-Smith et al，2018)和震群(Frank et al，2018；Xue et al，2018)的遗漏地震检测工作中；第三类基于人工智能的方法，以大量地震记录作为训练集学习地震波形的特征，进而在连续波形上搜寻这种特征从而检测微震信号。国内外学者提出了多种人工智能地震和震相检测算法(Perol et al，2018；Ross et al，2019；于子叶等，2018；赵明等，2019a、2019b；蒋一然等，2019；李健等，2020)，在地震序列分析和区域地震活动性研究中得到了应用(Ross et al，2020)。

本文以2017年3月渤海地震序列为研究对象，首先检测序列中遗漏的微震事件，其次基于地震精定位和震源机制计算结果，对序列发震构造进行分析。

1 微震检测

由于渤海地震序列地震波形相似度较高，同时序列发生在几天之内，使得波形互相关计算量不大。因而，本文采用模板匹配识别方法(Peng et al，2009)，检测渤海地震序列目录遗漏的地震事件，从而提高地震序列的完整性。

模板匹配方法以较大地震记录波形作为模板，在连续波形上进行互相关扫描，互相关达到一定阈值即为检测到地震事件。遗漏地震事件检测的计算步骤为：①选取地震目录中所有ML2.0以上地震事件及其震相报告，选择距离震中最近的3个宽频带台站CLI、BDH和LUX的三分量数字化地震波形，截取地震事件S波到时前2s至后2s的波形作为模板；②将模板波形和连续波形通过2～8Hz的带通滤波，再由每秒100个采样点降为20个采样点；③将模板波形在连续波形上进行滑动窗互相关扫描，滑动步长为1个采样点(降采样后采样点间隔为0.05s)，计算3个台站波形互相关系数的平均数；④取平均互相关系数大于9倍绝对离差中位数(MAD)，以及至少有一个台站单台互相关系数大于0.65作为阈值，挑选大于阈值的点作为检测到的地震事件，其中台网目录中未记录的地震作为遗漏地震事件；⑤利用遗漏事件水平向波形S波段最大振幅与模板地震的振幅比计算二者的震级差，从而估计遗漏事件的震级。

使用表1 中的6个ML2.0以上地震事件的波形作为模板，对3月24—28日时间段内的连续波形进行互相关扫描。将3月27日23时37分46.7秒的ML2.6地震作为模板，其在3月27—28日进行互相关扫描的结果，如图2 所示。波形互相关值为3个台的平均互相关系数，9倍MAD作为阈值(图2 中红色虚线)。图2(a)结果显示，在27日23时至28日4时检测到较多的地震事件，显示出在此时间段内多次发生地震，且其波形相似度较高。

图 2 模板匹配微震检测计算结果(a)扫描得到的3个台的平均波形互相关系数，红色虚线表示通过计算9倍MAD得到的地震检测阈值；(b)互相关系数分布统计，纵坐标与图(a)相同，统计区间步长为0.01，横坐标为每个统计区间内的互相关系数个数，取对数坐标

图 3给出了一个检测到的目录遗漏地震波形对比实例，所使用模板为3月28日3时19分55秒的ML2.3地震，截取模板S波到时前2s至后2s波形，经过2～8Hz带通滤波(图3 中红色波形)。检测到目录遗漏地震为3月27日23时37分22秒的ML1.8地震。由于其后约25s发生一次ML2.6地震，振幅远大于ML1.8地震，且震中距较远台站(如DOY、SUZ、DL2)的2个地震的记录波形有所交叠，造成人工识别ML1.8地震较为困难。

图 3 检测到的目录遗漏地震与模板地震波形对比红色波形为模板地震；黑色波形为检测到遗漏地震时间段的波形；红色波形对比出的位置为检测到遗漏地震波形的位置；由于需要清晰显示ML1.8地震波形，ML2.6地震波形有所限幅

经检测后较完整序列目录的M-t图如图 4(a)所示，共检测到目录遗漏的地震事件32个，约为台网目录给出地震数量的1.8倍，其中检测到ML2.0以上的遗漏地震2个。检测到的遗漏地震较为集中地发生在3月24日ML4.4地震后15个小时内，以及3月27日23时至28日凌晨4时的时间段内(图4(b))。

图 4 经模板匹配检测后较完整的序列目录M-t图(a)3月24—28日的M-t图，其中红色表示第一组地震序列，蓝色表示第二组地震序列，黑色表示第三组地震，实线表示台网目录中给出的地震，虚线表示模板匹配方法检测到的遗漏地震；(b)为图(a)中3月27日23时至28日4时的放大图，其中蓝色实线表示台网目录中给出的地震，蓝色虚线表示模板匹配方法检测到的遗漏地震，均属于第二组地震序列

通过模板匹配计算，我们观测到此次渤海地震序列的地震可以分为2组，一组为主余序列，一组为震群。图4(a)中红色所示的第一组地震序列为ML4.4地震的余震序列，共检测到14次地震，其地震频度和强度均随时间而衰减；蓝色所示的第二组地震序列，集中发生在3月27日23时至28日4时之间，共检测到35次地震，其中最大地震震级为ML3.5，次大地震震级为ML2.6，符合震群的定义；另外，台网目录中给出的3月27日9时7分22.5秒的ML1.6地震(图4(a)中黑色实线所示)无法被其他任何模板地震检测出，因而不属于上述2组序列，判断为一次孤立地震。此外，震群活动结束后的3月28日9时38分3.3秒仍有一个属于第一组的地震发生，因而这2组地震并非先后发生，而是震群发生在余震序列活动期间。

第一组与第二组地震每组内地震的波形互相关系数较高，而2组地震间的波形互相关系数较低，其中一组的地震无法被属于另一组的模板地震检测出。选取记录信噪比最好的BDH台垂直向波形，经过2～8Hz带通滤波后计算互相关系数，比较两组地震序列内ML1.0以上地震间波形互相关系数占比。如图5 所示，第一组序列除了模板ML4.4地震外，共有10个ML1.0以上地震，可计算得到10个互相关系数，根据互相关系数cc值分档，其中满足0.4≤cc<0.5、0.5≤cc<0.6和0.6≤cc<0.7的各1个，占比均为10%，满足0.7≤cc<0.8的5个，占比50%，满足0.8≤cc<0.9的2个，占比20%，10个互相关系数平均值为0.7033；第二组地震序列互相关系数平均值为0.8580。震群地震间的波形互相关系数平均值(0.858)高于余震序列地震间的平均系数(0.7033)，显示出震群的波形一致性相对较高。

图 5 第一组(a)和第二组(b)地震序列地震间波形互相关系数占比

2 震源机制计算与地震精定位

首先使用CAP方法(Zhu et al，1996)计算序列中2个震级最大地震的震源机制，选择渤海周边河北、山东和辽宁地震台网的宽频带地震计记录波形，体波滤波频段为0.05～0.2Hz，面波滤波频段为0.05～0.1Hz。参考基于深反射剖面对渤海湾地壳结果的研究(张成科等，2002)设定速度结构(表2)，vP/vS波速比设为1.73。

表2地震精定位和震源机制计算所用速度结构模型

ML4.4地震的震源机制计算结果如图6 所示，节面I走向51°、倾角37°、滑动角-109°，节面Ⅱ走向254°、倾角55°、滑动角-76°；最佳拟合深度约为19km，拟合矩震级为MW4.08；结果显示ML4.4地震震源机制以拉张为主。ML3.5地震的震源机制计算结果如图7 所示，节面I走向41°、倾角90°、滑动角-18°，节面Ⅱ走向131°、倾角72°、滑动角180°；最佳拟合深度约为13km，拟合矩震级为MW3.49；结果显示ML3.5地震震源机制以左旋走滑为主。震源机制计算结果给出2条正交的节面，判断哪个节面为实际发震断裂，需要结合地震序列精定位结果。

图 6 ML4.4地震震源机制计算结果(a)波形拟合；(b)震源深度拟合

图 7 ML3.5地震震源机制计算结果(a)波形拟合；(b)震源深度拟合

为了提高地震精定位结果可信度和精度，基于波形互相关对地震序列震相相对到时进行校正。本次地震序列发生在渤海中部(图1)，100km内仅有2个短周期井下台，且其记录波形质量较差。最近的宽频带基岩台(CLI台)距震中约110km，该震中距范围极易出现Pg、PmP、Pb震相交叠的现象，导致准确识别震相较为困难，因而台网给出的震相报告存在人工识别震相误差。我们采用波形互相关震相检测技术精确拾取地震间同台震相到时差，对所有目录给出的地震和检测到的遗漏地震的震相重新标定了到时。

震相校正的计算步骤为：①人工校核地震序列2组地震中最大的ML4.4和ML3.5地震的P波和S波震相到时，作为2组地震计算波形互相关的模板地震；②对于序列中的其他小地震，以微震检测计算中得到的发震时刻为基准，根据模板地震的走时计算小地震震相到时的初始值；③将模板地震和小地震记录波形经过2～8Hz带通滤波，截取模板地震各台站P波和S波到时前0.5s至后1.5s波形，以及小地震各台站P波和S波到时初始值前0.5s至后1.5s波形，做波形互相关计算；④若波形互相关函数出现明显的尖峰，则认为检测到相应震相，作为校正后的小地震震相到时，若未出现明显尖峰，则认为未检测到震相；⑤将所有基于波形互相关检测到的震相到时数据作为校正后的新震相报告。

基于校正后的新震相报告，利用双差定位法(Waldhauser et al，2000)对序列进行精定位，使用的地壳速度结构见表2。ML4.4余震序列地震震中初始位置统一设置为台网给出的ML4.4地震震中位置，初始深度设为CAP方法拟合的ML4.4地震深度19km；ML3.5震群地震震中初始位置统一设置为台网给出的ML3.5地震震中位置，初始深度设为CAP方法拟合的ML3.5地震深度13km；反演计算选用奇异值分解(SVD)方法。

最终得到的精定位结果如图 8 所示，ML4.4余震序列有5个地震可给出精定位结果(图8(a)中红色圆圈所示)，ML3.5震群有25个地震可给出精定位结果(图8(a)中蓝色圆圈所示)。其他地震因可用波形互相关检测到的震相数量较少，无法得到高精度的定位结果，故在计算中被舍去。双差定位法给出了地震序列的反演计算误差，其中ML3.5震群精定位结果震中平均误差为16.1m，最大误差为30.1m，深度平均误差为26.3m，最大误差为37.9m；ML4.4余震序列精定位结果震中平均误差为63.3m，最大误差为328.1m，深度平均误差为101.1m，最大误差为386.9m。由于反演使用奇异值分解方法，双差定位法给出的结果误差是有意义的。

图 8 渤海地震序列精定位结果(a)震中分布；(b)ML3.5震群在AA′剖面上的投影；(c)ML3.5震群在BB′剖面上的投影

结果显示，2组地震发震构造走向均为NE向，与震源机制解中节面I走向基本一致。图8(b)给出了ML3.5震群在与断层走向正交的AA′剖面上的投影，显示其发震构造为近直立的断层，与震源机制解节面I倾角90°的结果一致。ML4.4余震序列精定位的地震仅有5个，无法刻画断层几何特征，故未给出在剖面上的投影图。

3 讨论

2017年渤海地震序列发生在张渤带与郯庐带的交汇部位。张渤带是一条晚第四纪、具有相当规模和较强活动性的NW向活动构造带(徐锡伟等，2002)。郯庐带是我国东部规模最大的一条巨型走滑断裂带，该断裂走向为NNE-NE向。地震精定位结果显示渤海序列走向为NE向，其发震构造与郯庐带走向接近。郯庐带渤海湾段历史上发生过多次强震，其走向在山东段和渤海湾南段为NNE向(曹筠等，2018；顾勤平等，2020)，渤海湾北段为NE向(陈瑛等，2020)。地震序列震中所在的渤海中部正处于郯庐带走向由NNE向NE转变的过渡地带。震源机制解给出的2次最大地震NE向节面走向分别为51°和41°，其更接近于渤海湾北段的走向。

渤海海域内郯庐断裂带东西两侧的构造存在差异，断层西侧多为高速区，中强地震震源深度浅于断层东侧(汪晟等，2017)。谢卓娟等(2008)通过统计分析认为渤海海域内地震深度多在10～20km的地壳中，渤中拗陷内地震震源深度一般较深。渤海地震序列震中位置处于郯庐带西侧(图1)，CAP方法给出的2个最大地震深度为19km和13km，精定位后序列地震的深度均不超过20km，与谢卓娟等(2008)给出地震深度分布基本一致。1969年渤海海域M7.4地震的深度尚存在争议，束沛镒等(1983)通过远震波形反演得到震源深度为25km，而周蕙兰等(1989)同样使用远场波形得到初始破裂点深度为8km，不在谢卓娟等(2008)得到的中小地震震源深度分布优势区间内。随着渤海海域地震监测能力提升，渤海小震、微震震源深度的测定精度也将不断提升。

此次渤海地震序列震中位于渤中盆地内的渤中凹陷。渤中盆地为NNE向“郯庐深断裂”与NWW向“北京-塘沽-蓬莱深断裂”的交汇区域，盆地构造演化受这2条重要的区域构造带影响。渤中凹陷是盆地新生代沉降中心，充填的新生代地层达10km以上(陆克政等，1997)。渤中凹陷区断层NE、NW和EW三组方向，主要分为2类：一类以延展断层为特征的控盆构造，另一类以花状构造和产状较陡为特征的与走滑作用有关的次级断层(侯贵廷，2014)。渤西地区拗陷的主伸展断层比较低缓，倾角多为30°(童亨茂，2003)。倾角较低的伸展断裂与ML4.4地震震源机制解节面I给出的参数(倾角37°、滑动角-109°)较为吻合；倾角较陡走滑型次级断层与ML3.5地震震源机制解节面I给出的参数(倾角90°、滑动角-18°)较为吻合。

综合地震精定位和震源机制计算结果，我们推测ML4.4余震序列发震构造为渤中凹陷内NE向低倾角的伸展性正断层，ML3.5震群发震构造为NE向倾角较陡的次级走滑断层。根据陆克政等(1997)和周斌等(2000)给出的渤中盆地构造剖面图，得到渤海地震序列发震构造示意图，如图9 所示。ML4.4余震序列发生在伸展性正断层上，浅部倾角较大，深部倾角逐渐变小。ML3.5震群发生在次级走滑断层上，倾角较陡。由于中小地震破裂尺度较小，很难将其发震构造与已知的断层相互对应。虽然可以对发震构造的几何特征做出一定判断，但确定地震发生在哪条断层上仍十分困难。

图 9 渤海地震序列发震构造示意图(a)渤中盆地构造剖面图，N+Q：新第三系和第四系，Ed：东营组，Es1-2：沙一、二段，Es3：沙三段，Es4：沙四段，据陆克政等(1997)修改；(b)视角为E向，断层走向为NE向；断层1表示ML4.4余震序列发震断层，红色圆圈为精定位后的余震序列震源位置；断层2表示ML3.5震群发震断层，蓝色圆圈为精定位后的震群震源位置

各类微震检测方法均有其优势与不足，在实际应用中应根据研究对象和研究目标选择使用的方法。长短窗类方法在信噪比较高时能够可靠准确地识别震相(刘希强等，2009)，计算速度较快，但会漏掉低信噪比的信号(赵明等，2019b)。模板匹配类方法适用于检测具有相近震源位置和震源机制的地震，对于波形高度相似的地震序列较为有效(Schaff，2010)，但其检测精度取决于使用的模板，计算量随模板数量增大呈指数级增加(Skoumal et al，2014)，且无法检测到不与任何模板地震波形高度相似的地震。人工智能类算法优势在于训练好的模型具有泛化能力，可以发现训练集不包含的新特征(赵明等，2019b)，在识别准确率和效率等综合性能方面有所提升。

利用地震波形的相似性，除了能实现检测微震的功能，还可以对序列中的地震进行聚类分析(王伟涛等，2012)。本文对渤海地震序列的分析，通过模板匹配算法检测出序列可分为2组，并对2组地震分别进行精定位，发现2组地震的发震构造确实存在较大差别。若未经过波形互相关聚类，将2组地震合在一起分析，则可能无法得到可靠的发震构造分析结果。不可否认，基于人工智能的微震检测技术在检测能力和计算效率上显示出了越来越大的优势(赵明等，2019b；Perol et al，2018)，但在地震序列研究中，由于地震波性相似度较高，地震活动时间相对较短，模板匹配方法对波形相似的地震检测能力较强的优势体现较为明显，同时能避免检测出不属于地震序列的其他地震的干扰。同时，与GPU加速技术(Beaucé et al，2017)、音频指纹技术(Yoon et al，2015)等新技术相结合，在一定程度上弥补了模板匹配类方法计算效率较低的不足。通过给模板地震各台站走时加入一个微小的变化值，M&L方法(Zhang et al，2015)或弱模板识别方法(吴梦羽等，2018)可以检测到更多与模板地震“不太相似”的微震。因而，模板匹配方法仍然广泛应用于地震序列微震检测工作中(Ross et al，2019；Yao et al，2020；Sianipar，2020)。

对地震序列精定位结果误差的估计对发震构造分析可靠性有重要作用。自助法(Bootstrap)是估计误差的有力工具。对震相到时数据加入一定的随机误差(王清东等，2015)，或对震相到时、地震初始位置和速度结构同时加入误差(Hardebeck，2013；姜金钟等，2016)，通过重复定位几百次确定精定位结果的误差。Hardebeck(2013)分析了基于波形互相关校正后的相对到时误差，发现其绝大多数在0.005s内。本文采用每秒100个采样点的数据，0.005s为采样间隔的一半。因而，可以认为经过波形互相关震相检测后，每秒100个采样点数据的相对到时误差是非常微小的，对每一组地震相对定位结果的影响很小。本文2组地震的初始位置分别置于台网给出的ML4.4和ML3.5震中位置，震源深度基于CAP拟合结果。鉴于人工识别的震相到时存在一定误差，导致震中定位结果存在一定不确定性，ML4.4和ML3.5两个模板地震初始位置的不确定性会造成本文精定位结果中2组地震的相对位置存在误差。同时，速度结构的不确定性，尤其是震源区P波和S波速度的误差，对双差定位也有影响。利用Bootstrap方法定量估计模板地震初始位置和速度结果造成的相对定位结果误差，是今后研究和改进的方向。