2012年以来新疆于田强震前震序列的时空演化模式

闫坤， 王伟君， 冯磊， 彭菲， 寇华东， 袁爱璟， 付广裕

1 中国地震局地球物理研究所， 北京 100081 2 中国地震局地震预测研究所， 北京 100036 3 新疆维吾尔自治区地震局， 乌鲁木齐 830011 4 中国地质大学(北京)， 北京 100083

0 引言

前震现象在全球板块边界、板内断裂等活动构造带中广泛存在(Jones, 1984; Dodge et al., 1995; Marzocchi and Zhuang, 2011; Guilhem and Nadeau, 2012; Kato et al.,2012; Ruiz et al., 2014; Schurr et al., 2014; Huang et al., 2020).据统计，在1966—1996年间中国大陆发生的强震中，约5%具有前震(Chen et al., 1999).强震前的地震活动性图像，如前震(foreshock)，震前平静(preseismic quiescence)、前兆性震群(precursory swarms)和地震空区(gap/doughnut patterns)，可以为断层带的应力状态或地震预测提供一些重要的线索(Kanamori, 1981).

研究前震时-空演化规律有助于提高我们对地震发生机理的认识(Gomberg, 2018).目前普遍认为地震成核可能存在断层预滑(Dodge et al., 1996；Bouchon et al., 2011；Kato et al., 2012, 2016)、级联触发(Helmstetter and Sornette, 2003；Marzocchi and Zhuang, 2011；Ellsworth and Bulut, 2018；Yoon et al., 2019)或两种机制综合驱动(McLaskey, 2019; Yao et al., 2020)三种模式.除了岩石实验(Latour et al., 2013; McLaskey, 2019)、数值模拟(Dieterich, 1992; Noda et al., 2013; Kazemian et al., 2015)、统计学分析(Kamogawa et al., 2019)以及地表连续GPS观测(Ruiz et al., 2017; Radiguet et al., 2016; Socquet et al., 2017)等手段，通过追踪前震的时-空演化能够获得临震前岩石破裂和应力状态等重要信息，是探究地震成核过程的关键.

研究前震活动规律需要完整的地震目录和较高的定位精度(Sánchez-Reyes et al., 2021; Mignan, 2014)，稀疏的地震监测台网成为前震研究的主要障碍.依赖人工识别震相的检测方法在低信噪比的情况下容易遗漏地震，而基于波形相关的模板匹配法(Shelly et al., 2007; Peng and Zhao, 2009)能够自动化对地震进行检测，降低人为因素影响，显著提高目录的完整性，从而更好地展现前震活动的迁移规律(如Bouchon et al.,2011; Kato et al., 2012; Walter et al., 2015; Huang et al., 2020)，获得更准确的地震活动性参数(Huang et al., 2020)

自1997年以来，巴颜喀拉块体周缘强震频发，陆续发生了1997年玛尼M7.5、2001年昆仑山口西M8.1、2008年于田M7.3、2008年汶川M8.0、2010年玉树M7.1、2013年庐山M7.0、2017年九寨沟M6.6以及2021年玛多M7.4等地震(图1a)，为强震过程研究提供了丰富的震例.特别是巴颜喀拉块体西边界的于田地区相继出现4次6级以上的地震(其中2008年和2014年为7.3级地震)，反映出巴颜喀拉块体西部正在进行强烈的应力加载和调整.

图1 研究区概况

近些年来于田几次强震活动总体上分布在阿尔金断裂西南尾端弧形过渡带中.全长约2000 km的阿尔金断裂是调节印度板块向北挤压应变(Li et al., 2016)的一条大型左旋走滑断裂，其东部滑动速率约20 mm·a-1，西端滑动速率约30 mm·a-1(Tapponnier et al., 2001).于田西部区域在拉张变形作用下NE向走滑断裂逐渐演化为帚状排布的多条NEE左旋走滑羽列状分支断裂，包括北部的康西瓦断裂，南部的阿什库勒断裂和硝尔库勒断裂等，并由SN方向上的正断层连接，形成了硝尔库勒、南硝尔库勒和东琼木孜塔格等活跃的拉分盆地(徐锡伟等，2011；宋春燕等，2015；李海兵等, 2015; 武艳强等，2020; 袁兆德等，2021).在这些拉分盆地及其周缘容易闭锁并累积应变，从而孕生破坏性地震(宋春燕等, 2015).

2008年于田M7.3地震开启了阿尔金断裂西端一系列强震的序幕，这些地震都有比较明显的前震活动，如2014年于田M7.3地震前有一次M5.4前震，2020年于田M6.4地震前有M4.7前震，可以为地震成核过程研究提供难得的资料.但是当地地震监测台网稀疏，初始目录中的前震活动较不完备，难以刻画出其时-空迁移特征以及与主震的关联.因此，本研究将通过波形模板匹配技术和单台定位方法，完善该地区的地震目录和定位，再以此为基础，深入分析2012年以来3次强震序列的时-空迁移规律，探讨该区域的强震发生模式.

1 数据和方法

1.1 数据

本文以2012年以来于田三次强震的震源区为研究区(35.5°N—37°N，82.1°E—83.4°E)；该区域地震台站稀疏，最近的于田台从2008年下半年才开始正式运行，是2014年于田M7.3主震150 km范围内唯一的固定台站(图1).该台站为三分量宽频地震仪，采样率为100 Hz，于2012年8月27日对仪器的南北方位进行了重新标定，标定前有-3.1°的误差，校正后误差小于0.1°.由于台网稀疏，很多小震(特别是前震)无法用常规方法进行检测或定位，因此考虑只利用于田台进行事件波形模板匹配检测和单台定位；重新检测的时间为2011年6月1日—2020年12月31日.

1.2 微震波形模板匹配检测

数据处理和检测流程参照Wang等(2015).从2012至2020年2月中国地震台网的697个地震目录中，选择P波信噪比≥4的事件制作波形模板.对模板和连续波形进行1～8 Hz带通滤波后，用模板P波和S波到时前1 s至后4 s的波形分别与一天长度的连续波形垂直分量和水平分量计算相关系数，并逆推至发震起始点进行三分量相关系数的叠加；然后以0.01 s为步长，滑动计算得到一天长度的互相关曲线，从中计算绝对偏差中位数(MAD).

单台检测受随机噪声干扰影响较大，为了降低误检率，需要提高检测阈值标准；但过高的阈值，会降低地震的检出率.经过测试，我们采用17倍MAD为阈值.由于震中距较大，地震尾波较长，为了避免重复检测，我们设置紧临事件的时间间隔要大于等于6 s.新检测事件的位置和模板相同，其震级通过新事件与模板之间的最大振幅对数比确定：

(1)

其中M、M0分别表示待测事件与模板事件的震级；A、A0分别为待测事件与模板事件的最大震振幅.

全部新检测的事件波形将被截取出来进一步做目视检查，并将其中误判为地震事件的瞬态噪声或仪器标定脉冲信号剔除，保证结果的可靠性.

完成检测后，对地震目录事件波形进行互相关分析，评估其相似性并进行分类.流程为首先截取地震事件P波后20 s长度的垂直向波形，进行1～8 Hz带通滤波后再计算互相关矩阵.根据最小距离层次聚类算法(single linkage, Thompson et al., 2017)生成平均相关系数的分层聚类树，最后用0.9和0.95阈值对事件分组.此外，人工读取P波S波到时，约束发震位置.其中波形相似度大于0.95，S-P波相对到时差小于0.1 s的地震对认为是相似地震.

1.3b值计算

b值大小与应力状态关系密切，通常应力的集中会导致b值下降(Spada et al., 2013; Scholz, 2015)，因此可用于评估孕震断层的应力状态.地震目录完备性对b值计算有一定的影响(Zhou et al., 2018)，相对完备的新检测地震目录能够更可靠地获得b值动态变化.我们从2011年6月1日开始，以250个事件作为计算窗口，10个事件为滑动步长，计算b值曲线(Gulia and Wiemer, 2019).为计算每个窗口的b值，先用最大曲率法计算出完备震级，统计大于完备震级的事件，如满足50个以上事件，则用最大似然法(Aki, 1965)计算：

(2)

其中Mc为最小完备震级，Ma为大于完备震级地震的平均震级，由于震级不是连续的，其最小宽度Δ为0.1，上式对震级的非连续性进行了校正.b值标准偏差为(Wiemer, 2001)

(3)

其中M为样本地震震级，Mmean表示平均震级，var表示方差.

1.4 事件单台定位

由于研究区台站分布稀疏而且台站距离震源超过100 km，信噪比低，可拾取的震相数据有限，很多地震无法用常规定位方法确定地震位置.这种情况下，可以利用单台三分量波形，通过计算地震波入射的背方位角和S-P到时差，单台确定地震位置.由于离震中最近的台站接收到的地震波形质量相对较高，单台定位可以获得更多的地震定位；尽管难以约束深度信息，但足于辨别出地震在平面空间的相对位置，有利于识别地震的空间迁移特征.

于田台离震源区最近，有大量事件的波形能够可靠地获得S-P到时差和地震背方位角.事件到台站的距离L可以用下式估算：

L=(tS-tP)·VP(VP/VS-1)，

(4)

其中tP、tS分别为P波和S波走时；VP、VS分别为P波和S波速度.

震源区附近的人工地震测深结果(李秋生等，2001)将地壳速度结构分为6层，每层的顶界面深度分别是0、12、20、30、40、59 km，对应的层速度分别为5.4、6.0、6.3、6.6、6.9、8.1 km·s-1.根据附近台站的接收函数h-k扫描结果，这个区域S/P波速比为1.73(刘文学等，2011).因为于田地震主要分布在15 km以上(房立华等, 2015; 罗钧等, 2021)，地表至地下20 km的平均VP为5.48 km·s-1，因此L=7.507×(tS-tP).深度对震中距有一定的影响，但震中距越大，其影响越小.在70 km震中距，10 km深度差异导致的震中距误差约为0.8 km，而该地区双差定位方法的水平误差约0.8 km(房立华等, 2015)，贝叶斯绝对定位误差约1.2 km(张广伟等, 2014).

背方位角主要通过入射P波进行估算.估算方法有多种，本文主要使用协方差矩阵偏振分析方法(Vidale, 1986).将三分量波形旋转后，其径向u(t)、切向v(t)和垂向w(t)组成的协方差矩阵C(t)，代表了地震波在三维空间上的能量分布：

(5)

其特征方程为

(6)

λi为特征方程的解.调整特征向量(xi,yi,zi)在三维特征空间中的方向α(0°～180°)，当特征向量实部的长度X(下式)达到最大时，该方向的极性也最大.

(7)

x0、y0表示特征值λ0最大时特征向量(x0,y0,z0)的水平元素，cisα代表cosα+isinα，Re为向量实部.地震事件的背方位角Φ如下：

(8)

协方差矩阵偏振分析方法自由参数只有计算偏振的滑动时间窗长，在本研究中设为4 s.为保证结果的稳定性，我们首先对事件三分量事件波形做0.01～8 Hz滤波，然后自P波初至时刻起以0.01 s步长滑动进行极性分析，取前10个结果的平均值作为该事件的背方位角.波形信噪比对计算结果的稳定性会有一定的影响(图2)，我们取P波前1～2 s和P波后0.5～2.5 s时间窗计算信噪比，重点分析信噪比大于5的事件.图2也对比给出了三分量波形短时互相关方法(Roberts et al., 1989)的计算结果，可见偏振分析方法得到的背方位角结果一致性更高.由三分量波形计算背方位角存在180°的不确定性，但考虑到于田2012年以来三次强震序列基本上发生在于田台的东侧(图1b)，大于180°的背方位角应该投影回东侧.

图2 不同方法计算的背方位角结果

2 结果

2.1 微震检测结果

对0级以上检测事件进行复检并剔除误检测后，最终获得了2011年6月—2020年12月研究区内21411个新事件，是台网目录的约15倍，它们的M-t图如图4a示.采用最大曲率法得到的最小完备震级Mc为1.3(图3b)，相比于台网目录的2.8(图3a)，地震目录完备性显著提升.

检测结果表明台网目录有较多的遗漏，包含一定的3～4级地震.图3c是部分遗漏事件在于田台的波形，它们都有明显的事件特征，P和S波震相明显、信噪比高.重新检测后前震的活动图像得到了较好的补齐，如2012年M6.3地震前7个月，新发现了23个地震事件；又如台网目录中2014年M5.4震后至M7.3主震之间台网目录只记录到7次事件，而新目录中出现多达103个地震；此外，我们发现2020年M4.7前震与M6.4主震之间并非平静，还存在7次事件，最高震级为2.4.总的来说，新的目录可以更好地还原前震过程的细节(图4a).

图3 台网目录(a)和检测目录(b)的震级频度分布, 红色实线表示完备震级拟合曲线; (c)台网目录中遗漏的部分地震波形

研究区的背景b值大约为0.7(图4b)，b值曲线显示在2012、2014和2020年三次强震事件前均出现b值的显著降低(低于0.5，降幅大于35%)，临震前b值又逐渐回升，震后恢复至正常水平(图4b).由公式(2)可知，目录的最小完备震级对b值计算有一定的影响，台网目录虽然能够大致得到b值的变化趋势(如图4b)，但结果的形态和检测目录有显著差别；特别是主震前的b值变化，如2014和2020年主震前b值的异常幅度和持续时间，检测目录结果似乎更为可信；因此提高目录完备性对前震特征研究是有必要的.

图4 (a) 检测目录震级时间分布与累积频度, 绿色阴影表示b值出现明显变化的时段; (b) b值随时间变化曲线，红实线和蓝虚线分别表示检测目录和台网目录结果.由于计算窗口地震数量不满足计算要求，部分时间没有b值

震前b值起伏变化在很多强震中都有类似的现象，如2009年L′Aquila地震、2011年Tohoku地震、2016年Pawnee地震以及2019年Ridgecrest地震(Sugan et al., 2014; Tamaribuchi et al., 2018; Walter et al., 2017; Gulia et al., 2020).震前断层相对闭锁，应力集中并趋同，此时小震占比减小，b值降低；震后断层破裂、应力迅速释放，从而b值上升(如Tamaribuchi et al., 2018).因此，地震活动和b值都可作为断层应力状态的指示器(Walter et al., 2017).

2.2 波形相似性

对大于0级的21411个检测事件进行了互相关计算，得到了事件对的互相关值.三次主震前后一段时间的地震序列互相关矩阵如图5f所示，前震序列表现出高相似性，余震序列的相似性较差.2012年的前震序列与2014年和2020年的前震序列之间有较高的互相关系数(>0.6)，主要是因为这些强震最初的活跃区域都在于田台东南角阿尔金断裂北侧，有相似的震源参数.

2014年于田M7.3地震前震序列中互相关系数大于0.95的事件丛集有5组(图5a—5d).2020年M6.4前震序列中互相关系数大于0.95的有3组，分别出现在1月28日，3月1日和5月2日(图5e).这些丛集中地震的平均S-P相对到时差约0.01 s，应该属于相似或重复地震.2014年M5.4前震和M7.3主震之间的相似地震的复现间隔在1～6 h，活动持续到主震前5小时，之后进入震前的平静期.这些相似地震在时间上为主震前活跃，并没有跨越主震；丛集之间平均S-P相对到时差约0.1 s，空间上是相对集中的.2020年M6.4震前出现的相似地震均为时间间隔只有几个小时的双震型(图5e)，但三组之间波形差异较大，可能处于不同的构造环境.

图5 (a)—(e)为于田台记录的互相关系数大于0.95的丛集波形, 蓝色和黑色波形分别为模板和新检测事件; (f) 三次强震序列波形的互相关系数矩阵

3 讨论

于田地区台网稀疏，台网目录的地震定位不确定性较大.前人对2014年M7.3地震与2020年M6.4地震序列进行了相对定位研究(如房立华等, 2015; 罗钧等, 2021)，但是这些结果主要是针对余震，前震研究较少，而且不同研究得到的地震定位差异较大(如张广伟等, 2014; 唐明帅等, 2016).本文通过模板检测和单台定位，能够获得更为完整、而且空间位置有所约束的前震活动图像.

3.1 2012年M 6.3强震序列的时空迁移

2012年M6.3主震的震源机制解为正断型(图1)，单台定位结果表明其发生在南硝尔库勒盆地东缘的隐伏断裂上，与前人定位结果相比(刘建明等, 2016)，两个结果相对于于田台的背方位角相似，但单台定位结果到于田台的震中距更大(图6d).主震前的地震分布在主震NNE方向硝尔库勒盆地北侧边界康西瓦断裂上，距离于田台站50～65 km(图6b,d).由两组地震构成，第一组发生在震前150～200天，最大震级2.7，第二组发生于震前50天，最大震级2.3，两组地震间隔约100天.第二组地震活动结束后经历了20天的地震平静后发生主震(图6b).整体上地震活动的间隔期缩短，同时地震位置从硝尔库勒盆地内部转移至盆地东南缘构造，有显著的空间转移.但发生主震的隐伏断裂上无明显前震活动.

图6 2012年M 6.3地震序列距离于田台的震中距

余震序列的背方位角集中在140～150°(图6d)，到于田台的震中距与主震相比减小约15 km(图6a,c)，说明余震向北部的盆地核部迁移，但余震和主震前的早期地震活动在空间上有较大的差异(图6c)，它们不是在一条断裂带上，这也说明2012年主震的破裂可能没有交汇到阿尔金和阿什库断裂.

3.2 2014年M 7.3强震序列的时空迁移

2014年M7.3地震的前震序列有显著的阶段性迁移现象.早期地震活动散布在阿什库勒断裂以南(图7d)，距离于田台站的距离分别为75～85 km和45～60 km的两个区域(图7b)，从主震前210天持续到主震前15天.经过15天的平静期后地震活动转移并集中到南硝尔库勒断裂北段，主震前31 h发生M5.4前震(图7d)，后续地震沿南硝尔库勒断裂分布，并在主震发生前5 h再次进入地震平静期(图7a)，随后M7.3主震在阿什库勒断裂上成核(图7d).震后地质调查研究结果表明南硝尔库勒和阿什库勒断裂上均有地表破裂(如袁兆德等, 2021; 李海兵等, 2015)，我们的结果说明南硝尔库勒断裂北段在震前已经沿断裂自西向东开始活动(图7d)，主震发生后地震活动带继续向东扩张.

图7 2014年M 7.3地震序列距离于田台的震中距

在前震序列中，我们共找到5组地震相似性高于0.95的地震，按照各组第一个地震的发震时间排序(如图7d)，最早的相似地震开始于M5.4震后14 min，持续到主震前5 h(图7a)，各组内最大事件震级分别为2.3、2.7、1.5、2.2和1.9.以上5组相似地震位于南硝尔库勒断裂北段(组3位于阿什库勒-硝尔库勒断裂)，主震东南缘(图7d)，无明显的迁移.

3.3 2020年M 6.3强震序列的时空迁移

2020年M6.3地震是一次NNE走向正断型地震，罗钧等(2021)结果显示主震位于琼木孜塔格峰附近(图8d).我们的定位结果与之不同，本文通过背方位角和震中距计算，认为主震应该位于南硝尔库勒盆地东缘断裂的南端(图8d).首先M4.7前震的背方位角约179°(图8e)，说明起始破裂区在盆地南缘，另外余震的能量也主要来自于田台170°方向(图8d)，因此2020年于田地震序列发生在南侧的可能性更大.南硝尔库勒盆地东西两侧相对的正断构造在盆地南部交汇，也更易造成应变累积进而促进地震的发生.

图8 2020年M 6.4地震序距离于田台的震中距

2020年于田M6.4主震前约170天的地震活动分布在距离于田台50～60 km的硝尔库勒盆地内(图8c,d)，背方位角小于150°(图8d).这些事件主要位于2012—2014年余震区内，总体上有向西南迁移的趋势.震前36min，地震活动转移至南硝尔库勒盆地，在其西缘产生M4.7前震，后续地震集中在这一区域直至主震在盆地东缘成核(图8).

3.4 于田三次强震的前震模式

于田三次强震的前震过程具有一些相似的模式特征.表现为在主震前半年至一年，b值曲线开始逐渐下降，然后在临震前抬升(图4b).这一期间的地震活动分散在主震区外围构造上(主要为硝尔库勒盆地周缘断裂)，震级小、频度低.临震前，地震活动强度增加，并集中到主震发震断裂的相邻构造上(图7d,8d)，最终引发跨断层的级联破裂(如Sieh et al., 1993; Xu et al., 2006; 袁兆德等, 2021).以上现象表明，于田几次强震的产生是区域内多条断裂组成的构造系统上应力的传递、转移以及形变局部化的结果(如Kato and Ben-Zion, 2021).

在2014年于田地震序列中M5.4前震的后续事件未出现明显的空间扩散，它们集中在主震南侧的硝尔库勒断裂北段，强度逐渐减弱，在主震前5 h几乎停止(图7a).这一阶段的前震活动特征反映出孕震断层上应力的调整和累积.单台互相关得到的5组高相似性地震丛集，有较短的复现周期，导致这一现象的原因可能是流体活动导致的小尺度无震瞬变(如Huang et al., 2020).前震在空间上的聚集以及波形高一致性现象在海城地震、孟连地震中均有发现(Chen et al., 1999)，探究其背后的物理机制有助于推动地震预测理论的发展.

于田震源区普遍的前震现象可能与该地区的构造特征有关.巴颜喀拉块体西边界的应力场呈东西向拉张，可能导致北东向断裂上摩擦系数的降低，此外西昆仑至东昆仑一带的地热流体活动也能够促进前震的形成.

4 结论

本文通过对2012以来三个于田强震震源区进行单台的波形模板匹配扫描和地震事件的单台定位，检测出大量遗漏地震事件(约为台网地震目录数量的15倍，完备震级从M2.8下降为M1.3)，同时有效地约束了这些事件的相对位置，完善了它们的活动图像.这些工作不仅有利于我们对前震时空演化规律的研究，也提高了我们对发震构造的认识；如我们认为2020年M6.4地震应该发生在南硝尔库勒盆地南缘而不是之前普遍定位的琼木孜塔格峰附近.

于田主震前的地震活动，总体上是比较活跃的，它们的b值随时间的变化基本表现为长时间(几个月至1年)的持续下降，降幅可达35%，在临震前有明显的回升，较好地反映震前应力的调整和微破裂的发生，有一定的前兆预报意义.它们的时空活动图像存在外围向主震位置逐渐迁移的模式；这种现象在2014年M7.3地震前最为明显，前震的波形相似性高，刻画出前震在主震发生区域外围的移动、停滞和复发，反映了地震成核过程，有助于未来我们对成核机制的深入研究.

青藏高原作为印度—欧亚板块碰撞的产物，构造复杂，有大量强震活动，是研究地震机理的天然试验场.然而，这些地区的地震台站分布非常稀疏，严重地阻碍着这些宝贵震例的地震过程深入研究.本文在于田地区的初步研究，说明利用新技术可以较大程度地改善稀疏台站条件下地震活动图像的重建，这对于青藏高原或其他观测条件欠佳的地区具有一定的参考意义.

致谢感谢新疆地震局陈亮和深圳防灾减灾技术研究院王宝柱对本研究的帮助；感谢两位审稿专家提出的宝贵建议.文中波形模板检测工作基于中国地震局地震预测研究所集群系统.GMT6(Wessel et al., 2019)被用于文中图件绘制；GSIMO(Thompson et al., 2016)和ZMAP程序(Wiemer, 2001)分别被用于波形聚类分析和b值计算.