2019, Vol. 40
震群发生期间，许多地震波形相互重叠，难以区分，较小地震震相难以识别。若能自动识别地震事件，对于完善地震目录、深入研究地震序列的发展过程、全面分析地震活动性等将具有重要意义。目前，在地震自动识别方法中，基于互相关技术的模板识别匹配滤波方法（Matched Filter Technique，MFT）可在信噪比较低情况下提取微弱信号，得到广泛应用（Gibbons et al，2006；Shally et al，2007；Peng et al，2009），如：Gibbons等（2006）利用模板匹配滤波技术，在单通道记录上分辨出Kara Sea地震主震后4小时发生的1个微地震；Peng等（2009）将该方法运用到2004年帕克菲尔德6.0级地震，检测到地震发生后3天内大量余震事件，数量为地震目录的11倍，清晰刻画出余震的迁移特性；谭培毅等（2014）检测出河北蔚县震群、涿鹿震群地震目录中被遗漏的小微震，利用完善的地震目录重新探讨了震群发震机理；李璐（2016）完善了2015年尼泊尔地震目录，全面分析地震活动性；王宁等（2017）成功检测出唐山震群部分遗漏地震。
河北省测震台网现接收河北省属53个数字遥测地震台站及周边地区共166个地震台站提供的实时连续波形数据。随着地震台站密度的增加，地震监测能力有了较大提高，为更好地研究震群活动奠定了坚实基础。2017年以来，河北束鹿断陷盆地及附近地区地震活动有所增强，相继发生1月3日隆尧ML 3.8地震和9月4日临城ML 4.4地震。本文选取2017年9月4日—6日临城ML 4.4震群中ML≥1.0地震作为模板，采用模板匹配滤波方法，通过波形互相关检测遗漏的地震事件，并根据互相关结果，标定遗漏地震事件的P波、S波到时，估计震中位置和震级，结合地震精定位结果和震源机制解，初步分析临城地震序列可能的发震构造。
1 地质背景2017年9月4日起，河北省临城县发生系列地震活动（图 1）。此次临城ML 4.4震群发生于华北平原中部束鹿断陷盆地附近。束鹿断陷盆地为典型的半地堑，新河断裂为其主边界断裂，主体凹陷内部发育NE和NW至近EW向2组次级断裂，其中NE向次级断裂为数条斜列分布的正断裂，大多位于凹陷西部斜坡地带，NW至近EW向次级断裂亦为正断裂（徐杰等，1988）。断陷盆地附近多有中强震发生，1966年3月22日河北宁晋7.2级地震即发生于盆地内部。
|
图 1 观测报告给出的临城地震序列震中分布 Fig.1 The epicenter distribution of the Lincheng earthquake swarm from the observation report |
据《河北省测震台网地震观测报告》，2017年9月4日—6日河北省临城县共发生地震370次，其中ML 0.9以下地震339次，ML1.0—1.9地震24次，ML 2.0—2.9地震6次，ML 3.0以上地震1次，最大地震ML 4.4（图 1）。排除重叠地震，选取临城震群中震级较大且震相清晰的ML 1.0以上地震事件作为模板地震，共挑选出21条符合条件的模板地震，结果见表 1。
|
表 1 选取的模板地震 Table 1 Selected template events |
按照台站距离震中较近、包围性较好、记录的地震波形信噪比较高、波形清晰的特点，选取临城震群不同方位的LIC台、ZAH台、YUS台、XIY台、HNS台、XTT台等6个台站，采用模板匹配滤波技术（MFT），对连续波形进行互相关扫描，识别与模板事件波形具有一定相似度的地震事件。台站分布见图 2。
|
图 2 观测台站分布 Fig.2 Distribution of the stations |
使用选取的21条模板地震，逐一扫描9月4日—6日地震记录连续波形，计算模板事件波形与对应分量连续波形的互相关系数（Cross-Correlation coefficient），计算公式如下
| $ {\rm{CC}} = \frac{{\sum {_{{t_0}}^{{t_1}}\left\{ {\left[ {X\left(t \right) - \overline X } \right] \times \left[ {Y\left(t \right) - \overline Y } \right]} \right\}} }}{{\sqrt {\sum {_{{t_0}}^{{t_1}}{{\left[ {X\left(t \right) - \overline X } \right]}^2} \times \sum {_{{t_0}}^{{t_1}}{{\left[ {Y\left(t \right) - \overline Y } \right]}^2}} } } }} $ | (1) |
式中，t0、t1分别为用于计算互相关系数窗口的初始与结束时刻，X(t)、Y(t)分别为模板事件波形和连续波形的时间序列，X、Y分别为模板波形与连续波形的时间序列平均值。为提高计算速度，连续波形和波形模板分别经过重采样，采样间隔由0.01 s变为0.05 s，扫描窗长为模板波形长度，计算时每次在连续波形上移动1采样点（0.05 s），最终得到该分量上的一道互相关波形。
将三分向CC取平均，当平均互相关波形上某一时刻的互相关系数高于该阈值时，认为是一个新检测事件。阈值选取需要权衡，一般选择绝对中位差（Median Absolute Deviation，MAD）的9倍（Peng et al，2009）或12倍（Yao et al，2015）作为阈值。但是，若地震台站分布稀疏、波形信噪比较低，可根据需要适当提高阈值（如15倍MAD），以减少模板事件与噪声信号间的错误检测，由于临城地区台站密集度较高，文中阈值选择为9倍中位差（Peng et al，2009）。MAD表达式为
| $ {\rm{MAD}} = {\rm{Median}}\left({\left| {{X_i}} \right| - \overline X } \right) $ | (2) |
式中，Xi为第i个互相关系数序列，X为其平均值。给出模板地震Eq904072636连续波形的互相关扫描结果（LIC、YUS、ZAH台站叠加），结果见图 3，其中：扫描时间段为2017年9月4日04:00—08:00，阈值为0.261，互相关序列图中最大值为1.0（模板地震）。图中实心圆为目录地震，空心圆为地震事件检测结果，共检测到疑似遗漏地震4个。
|
图 3 基于匹配滤波技术的波形互相关扫描结果 Fig.3 Cross-correlation scanning result based on matched filter |
得到疑似遗漏地震事件后，使用原始连续波形和模板波形，同样进行互相关扫描，以确认遗漏地震，并搜索遗漏事件的P波和S波到时。在距震中最近的3个地震台站至少检测出一个P波或S波震相，则认为此次事件为遗漏地震（谭毅培，2014）。以模板地震Eq0904030811临城台（LIC）记录波形为例，扫描得到遗漏地震发震时刻为2017年9月4日06:02:54.3，利用垂向波形检测P波，水平向波形检测S波，在连续波形垂直向截取P波到时前2 s至到时后2 s数据，水平向分别截取S波到时前2 s至到时后2 s数据，扫描得到互相关序列中互相关系数最大值所在位置对应时刻即为遗漏事件的P波、S波到时，其中P波互相关系数最大值为0.64，S波互相关系数最大值为0.87，结果见图 4。
|
图 4 遗漏地震的P、S波到时 Fig.4 P- and S-wave arrival time of overlapped events |
若遗漏地震未被任何台站检测到P、S震相，则认为是误检测地震。利用匹配滤波技术波形互相关方法，共检测得到21个疑似遗漏地震，经验证，其中2个事件为某些干扰引起的波形异常，则此次临城震群共有19个遗漏地震。采用遗漏事件水平向波形S波段最大振幅与模板地震的振幅比估计遗漏事件震级。每个遗漏事件可更换几个模板地震计算震级，同一个模板地震也可计算若干遗漏事件震级，互相关系数越大，震级计算结果越可靠。通过计算，得到9月4日—6日遗漏地震震级范围为ML 0.0—1.2，仅9月4日06:02:54.5地震震级为ML 1.2，其余地震震级范围为ML 0.0—0.8，结果见表 2。
|
|
表 2 检测到的遗漏地震发震时刻和震级 Table 2 The origin time and magnitude of the detected events |
采用不同方法，对2017年9月4日—6日临城震群震中进行精定位，其中：①对于临城震群地震目录事件，采用HypoDD精定位方法进行重定位，精定位使用一维P波速度结构参数见表 3，其中vP/vS波速比设为1.73（蔡玲玲等，2017）；②对于有4个以上台站检测出震相的遗漏地震，采用双差定位方法给定震中位置；③其他遗漏地震事件震中，震中置于互相关系数最大的模板地震震中位置（Peng，2009）。
|
|
表 3 地震精定位所使用速度模型 Table 3 Velocity structured in earthquake relocation |
目录地震与遗漏地震事件的精定位结果见图 5(a)，可见精定位后震中呈近东西状条带分布，见图中AA′剖面，遗漏地震的加入强化了震中位置的条带化分布特征，显示临城震群发震构造走向可能近EW向。采用双差定位方法，计算震群地震震源深度，沿AA′剖面深度分布见图 5(b)，可见震源深度成丛集分布，优势分布在5.5—7 km深度范围内，表明震群主要发生在地壳中上部。由图 5可知，震群地震发生在NE向元氏断裂（F1）和近EW向隆尧断裂（F4）交汇部位附近，距隆尧断裂较远，初步推测发震构造可能是其他近EW向活动断层。
|
图 5 震群精定位分布 (a)目录地震和遗漏地震精定位结果；(b)震群沿AA′剖面深度分布 Fig.5 Relocation distribution of the earthquake swarm |
利用P波初动法，求取本次震群中ML 4.4地震（最大地震）震源机制解，其中：节面Ⅰ走向103°，倾角83°，滑动角-22°；节面Ⅱ走向196°，倾角68°，滑动角-173°；P轴方位角58°，仰角20°；T轴方位角151°，仰角10°；错动方式为走滑型，推测发震构造近EW向，由此判断节面Ⅱ为真实发震破裂面。
5 结论与讨论综上所述，可以得到以下结论：①基于模板识别匹配滤波方法，检测得到2017年9月4日—6日临城震群19条遗漏地震；②通过波形互相关震相检测技术，检测遗漏地震事件直达P波、S波到时，得到震中结果；③采用遗漏事件水平向波形S波段最大振幅与模板地震的振幅比，估计遗漏事件震级，发现9月4日06:02:54.5地震事件ML为1.2，其余地震震级范围为ML 0.0—0.8，地震目录得到进一步完善；④据临城震群精定位结果，结合最大地震震源机制解，认为此次震群主要发震构造可能为近EW向活动断层。
利用模板识别匹配滤波方法检测遗漏地震，提高了地震目录完整性，但由于震群中震级较小的地震波形信噪比较低，难以检测到遗漏的全部地震。另外，以完善后的地震目录为依据，结合地震精定位结果与震源机制、地质构造特征，可为发震构造可行性判断提供参考。
蔡玲玲, 李冬圣, 孟立朋, 等. 2017. 京西北地区地震重定位分析[J]. 地震地磁观测与研究, 38(3): 47-52. DOI:10.3969/j.issn.1003-3246.2017.03.009 |
李璐.台阵处理和模板匹配滤波技术在微弱地震信号检测中的应用研究[D].北京: 中国地震局地球物理研究所, 2016.
|
谭毅培, 曹井泉, 卞真付, 等. 2014. 2013年8月河北蔚县小震群遗漏地震检测与发震构造分析[J]. 地震学报, 36(6): 1022-1031. |
谭毅培, 曹井泉, 刘文兵, 等. 2014. 2013年3月涿鹿微震群遗漏地震事件检测和发震构造分析[J]. 地球物理学报, 57(6): 1847-1856. DOI:10.6038/cjg20140616 |
王宁, 杨春利, 李雪英, 等. 2017. 遗漏震相检测技术在2016年9月唐山震群中的应用[J]. 华北地震科学, 35(1): 38-43. DOI:10.3969/j.issn.1003-1375.2017.01.006 |
徐杰, 方仲景. 1988. 1966年邢台7.2级地震的构造背景和发震构造[J]. 地震地质, 10(4): 51-59. |
Gibbons S J, Ringdal F. 2006. The detection of low magnitude seismic events using array-based waveform correlation[J]. Geophys J Int, 165(1): 149-166. DOI:10.1111/j.1365-246X.2006.02865.x |
Peng Z G, Zhao P. 2009. Migration of early aftershocks following the 2004 Parkfield earthquake[J]. Nature Geoscience, 2(12): 877-881. DOI:10.1038/ngeo697 |
Shally D R, Beroza G C, Ide S. 2007. Non-volcanic tremor and low-frequency earthquake swarms[J]. Nature, 446(7133): 305-307. DOI:10.1038/nature05666 |
Yao D D, Peng Z G, Meng X F. 2015. Remotely triggered earthquakes in South-Central Tibet following the 2004 MW 9.1 Sumatra and 2005 MW 8.6 Nias earthquakes[J]. Geophys J Int, 201(2): 543-551. DOI:10.1093/gji/ggv037 |