地震目录是区域地震活动性研究和地震危险性分析的重要基础资料.编目是否完整，震源参数是否精确，直接影响后续研究结果的可信性和科学性(Woessner and Wiemer， 2005；冯建刚等，2012).强地震发生后或震群活动时，短时间内在震中附近或局部小区域往往发生大量地震，不同地震事件波形相互交叠，震级较小的地震震相难以清晰识别，造成地震编目遗漏一定数量的地震，给震后趋势判定和发震构造分析等工作带来一定不利影响.如何在现有观测条件下检测出遗漏的地震，并测定其震中和震级，是一个值得关注且亟需解决的问题.

在地震信号自动识别技术研究中，我国各级地震台网目前主要采用基于长短时间平均(STA/LTA)(Allen，1978；Baer and Kradolfer， 1987；Earle and Shearer， 1994)并综合多种信号特征的识别方法(沈萍等，2002；刘希强等，2009).此类方法在信噪比较高时能够可靠准确地识别震相(刘希强等，2009).另一类方法利用记录波形与已知地震事件的波形互相关识别地震信号(Schaff and Richards， 2004；Schaff and Waldhauser， 2005；Gibbons and Ringdal， 2006；Yang et al., 2009)，若假设噪声为高 斯白噪声，互相关方法又称为匹配滤波技术(Matched Filter Technique) (van Trees，1968)).已有学者使用匹配滤波技术进行非火山地脉动与低频地震的关系(Shelly et al., 2007；Tang et al., 2010)、中强地震余震序列完整性和余震分布特征(Peng et al., 2007；Peng and Zhao， 2009；Schaff，2010；Meng et al., 2012)等方面的研究.本文尝试使用匹配滤波技术识别震群活动时地震目录遗漏的地震.

华北北部是我国大陆地震活动强烈的地区之一.受区域张性断陷作用影响，首都圈及邻近地区活动断裂控制的地堑和半地堑断陷盆地是破坏性地震的主要活动场所.控制盆地发育和演化的NE向至NEE向断裂以及与其相互切错的NW向断裂是主要的发震构造和孕震构造(徐锡伟等，2002).2013年3月3日至5日，在河北省涿鹿县和怀来县交界地区发生微震震群活动，震群位置如图 1a所示.震群所在延怀盆地位于北京西北约120 km，地质构造较为复杂，属于具有左旋剪切性质的张家口—渤海构造带(徐锡伟等，2002)，历史上曾发生多次6级以上地震(国家地震局震害防御司，1995)，1995年7月20日发生一次M_L4.1地震(许向彤等，1999；2001).多年来不同学者利用地质调查(冉勇康等，1992；1998；方仲景等，1993)、人工地震勘探(张先康等，1996；刘保金等，1999)、地震精定位(朱艾斓等，2005；李乐等，2007)以及中-欧合作怀来数字地震台网观测(王培德等，1997)等手段对区域内断裂构造进行了较为详细的研究.

《河北省测震台网地震观测报告》给出此次震群发生地震45次，其中M_L0.0~0.9地震32次，M_L1.0~1.7地震13次，最大震级为M_L1.7.震中位置分布见图 1b.本文选取其中震级较大的地震波形作为模板，通过互相关扫描检测震群中目录遗漏的地震，根据互相关结果标定遗漏地震事件的P波和S波到时，对震中位置和震级做出估计，补充现有的地震目录，并进一步分析此次震群活动可能的发震构造.

本文所用数据来自首都圈台网记录的三分量数字化地震波形.首先选取M_L1.0以上的地震事件，连续波形和事件波形记录经过4阶零相位Butterworth滤波器2~8 Hz滤波.根据观测报告中的到时信息截取P波到时前0.5 s到S波到时后4.0 s的波形，挑选其中三分量平均信噪比大于3的作为备选波形模板(Template).噪声能量水平由P波到时前6.0 s到2.0 s的波形计算得到.若一个地震事件的备选波形模板不少于3条，则将此事件作为模板地震(Template Event).

本研究共挑选出满足条件的8个地震事件作为模板地震，波形模板48条，见表 1所示，每条波形模 板对应一个地震事件和一个台站.由于每个台站记 录的P、S波到时差不同，每条波形模板的长度有所不同.

将每条波形模板在3月3日至5日的连续记录 波形上进行扫描，计算互相关系数(Cross-correlation，CC). 为提高计算速度，连续波形和波形模板分别经过重采样，采样间隔由0.01 s变为0.05 s.扫描窗长为波形模板长度，扫描间隔为0.05 s(一个采样点).将三分向CC取平均，通过计算序列的绝对离差中位数(Median Absolute Deviation，MAD)搜索地震. MAD表达式为MAD=median(|X_i|－X).

本文取5倍MAD作为判别地震的阈值.假设样本序列为符合高斯分布的随机过程，其标准差σ≈1.4826×MAD(Peng and Zhao， 2009)，则单个随机变量大于5倍MAD的概率为

若样本均值μ为0，则P₅≈0.0003851.我们设定至少有3条波形模板的互相关值高于阈值，后续人工鉴别过程中每条波形模板推得的发震时刻允许有±3个采样点的误差，误差范围为7个采样点.长度为3天的0.05 s采样率连续波形共有5184000 个采样点，3天内误判地震个数期望值为(P₅×7)³×5184000≈0.099，即误判地震平均每3天约为0.1个，其可能很小.因此本文涿鹿震群遗漏的地震事件在统计上是比较可信的.

图 1 涿鹿震群观测台站及震中位置示意图 a)观测台站分布图.三角形为记录到本次震群地震波信号的台站，五角星为震群震中所在位置. 其波形用于互相关扫描的台站标注了台站名(b)观测报告给出的涿鹿震群地震震中位置分布. Fig. 1 Station and epicenter distribution of Zhuolu earthquake swarm (a)Station distribution，in which triangles are stations recorded the signal in the swarm， and the star is the location of Zhuolu earthquake swarm. The station′s names are written whose waveforms are used in cross-correlation scanning.(b)Epicenter distribution obtained from observation report.

表 1(Table 1)

表 1 选取的模板地震 Table 1 Selected template events 模板地震号 发震时刻(北京时间) 震中经度(°E) 震中纬度(°N) 震源深度(km) 震级(ML) 波形模板条数 Eq030241 2013-03-03 02∶41∶34.6 115.313 40.395 8 1.2 6 Eq030251 2013-03-03 02∶51∶24.2 115.304 40.401 5 1.2 5 Eq030318 2013-03-03 03∶18∶56.5 115.311 40.397 6 1.1 6 Eq030507 2013-03-03 05∶07∶26.3 115.303 40.392 7 1.1 7 Eq030919 2013-03-03 09∶19∶24.0 115.309 40.398 6 1.5 8 Eq031227 2013-03-03 12∶27∶42.5 115.303 40.397 6 1.2 3 Eq032228 2013-03-03 22∶28∶17.5 115.311 40.400 7 1.7 9 Eq050133 2013-03-05 01∶33∶10.1 115.310 40.389 6 1.2 4

表 1 选取的模板地震 Table 1 Selected template events

图 2展示了波形模板Eq030919ZHL在3月3日1 ∶ 00到4 ∶ 00连续波形上进行互相关扫描的结果.计算MAD得到阈值取0.217，3小时内共有37个CC大于阈值的点.这些点出现的时刻为检测到的疑似地震在ZHL台的P波到时，减去模板地震 到ZHL台的P波走时，即得到疑似地震的发震时刻.

图 1

Fig. 2

图 2 互相关扫描结果示意图波形模板为Eq030919ZHL，扫描时间段为2013年3月3日1 ∶ 00至4 ∶ 00.图中虚线为互相关匹配阈值，共有37个点互相关值大于阈值.Fig. 2 An example of cross-correlation scanning result Template is Eq030919ZHL，scanning time is from 1 ∶ 00 to 4 ∶ 00 March 3rd，2013. The dashed line is cross-correlation matched threshold. There are 37 points whose cross-correlation coefficient are higher than the threshold.

由于地震信号不是完全的随机信号，互相关系数一般不是非常理想的尖脉冲，因而需要对互相关扫描得到的疑似地震进行人工鉴别.设定在同一发震时刻有3个以上不同台站对应的波形模板扫描结果存在疑似地震，则判定此发震时刻存在地震事件，本文中称为遗漏事件.

波形模板互相关扫描得到遗漏事件发震时刻称为发震时刻初始值.由发震时刻初始值和模板地震P波、S波走时推得到时称为P波、S波到时初始值.本小节使用未经重采样的波形模板和连续波形，同样使用互相关扫描搜索遗漏事件的发震时刻、P波和S波到时.

图 3为搜索过程的示意图.首先从经过滤波的连续波形上截取发震时刻初始值之前一倍于波形模板长度、之后两倍长度的三分量波形，用波形模板在其上做互相关扫描.扫描窗长为波形模板长度，扫描间隔为0.01 s.通过得到的CC峰值时刻和P波走时推断发震时刻(图 3a).

接下来分别利用垂直向和水平向波形校正P波、S波到时.鉴于手动拾取震相到时存在一定误差，为截取比较完整的P、S直达波波形，在连续波形垂直向截取P波到时初始值前1.5 s起始，长度为4.0 s波形，在波形模板垂直向截取P波到时前0.5 s起始，长度为2.0 s的波形进行互相关扫描.扫描窗长为2.0 s，扫描间隔0.01 s(图 3b)；在连续波形两个水平向分别截取S波到时初始值前2.0 s起始，长度为5.0 s波形，在波形模板水平向截取S波到时前1.0 s起始，长度为3.0 s波形进行互相关扫描.扫描窗长为3.0 s，扫描间隔0.01 s(图 3c).扫描得到的互相关序列中CC最大值位置即为遗漏事件P波、S波到时.

本文使用双差定位法(Waldhauser and Ellsworth， 2000)对目录已给出的地震事件和遗漏事件进行精定位.所使用数据包括观测报告给出的震相到时数据，互相关扫描得到的遗漏事件震相到时，以及用匹配滤波技术计算得到的P波波段互相关系数和P波相对到时差.利用遗漏事件水平向波形S波到时后4 s内最大振幅与模板地震水平向波形S波波列最大振幅之比估计遗漏事件的震级.

利用匹配滤波技术在3月3日至5日共检测到地震观测报告中遗漏的地震事件52个，其发震时刻和震级估计结果见表 2.

图 3

Fig. 3

图 3 搜索遗漏地震事件发震时刻，P波、S波到时方法示意图灰色表示连续波形，黑色表示波形模板，波形经过4阶零相位Butterworth滤波器2~8 Hz滤波；下方为互相关系数，并标注了互相关最大值CCmax.(a)三分量叠加波形；(b)垂直向P波波形；(c)水平向叠加的S波波形.Fig. 3 The process of detecting the origin time，P-wave and S-wave arrival time of overlapped events Gray curves are continuous waveforms，black curves are template waveforms，the waveforms are filtered in 2~8 Hz by 4th order zero-phase Butterworth filter. The curve below the waveforms is cross-correlation coefficient sequence， and maximum value of cross-correlation coefficient(CCmax)is indicated. The waveforms in(a)are three component added，in(b)are P-wave of vertical component， and in(c)are S-wave of two horizontal components added.

表 2(Table 2)

表 2 遗漏地震事件发震时刻与震级 Table 2 The origin time and magnitude of detected events 发震时刻(北京时间) 震级(ML) 发震时刻(北京时间) 震级(ML) 发震时刻(北京时间) 震级(ML) 2013-03-03 01∶26∶37.22 0.4 2013-03-03 07∶28∶57.81 0.4 2013-03-03 13∶48∶24.93 0.6 2013-03-03 01∶28∶58.73 0.4 2013-03-03 08∶08∶07.70 0.1 2013-03-03 22∶03∶49.23 0.3 2013-03-03 02∶46∶39.16 0.4 2013-03-03 08∶09∶36.07 0.4 2013-03-03 22∶22∶23.58 0.5 2013-03-03 02∶47∶50.84 0.3 2013-03-03 08∶17∶50.39 0.4 2013-03-03 22∶23∶54.32 0.5 2013-03-03 02∶51∶51.96 0.6 2013-03-03 08∶27∶15.45 0.1 2013-03-03 22∶40∶54.57 0.2 2013-03-03 02∶52∶15∶10 0.3 2013-03-03 08∶28∶52.01 0.9 2013-03-04 05∶40∶09.80 0.3 2013-03-03 03∶25∶33.61 0.4 2013-03-03 08∶29∶08.52 0.9 2013-03-04 08∶46∶44.81 0.5 2013-03-03 03∶45∶41.35 0.6 2013-03-03 08∶29∶13.95 0.9 2013-03-04 09∶04∶23.94 0.5 2013-03-03 03∶46∶19.38 0.3 2013-03-03 08∶29∶21.85 0.8 2013-03-04 13∶23∶24.77 0.3 2013-03-03 03∶52∶50.56 0.4 2013-03-03 08∶30∶31.90 0.4 2013-03-04 13∶35∶37.57 0.2 2013-03-03 03∶54∶57.01 0.3 2013-03-03 09∶20∶51.63 0.5 2013-03-04 13∶55∶12.58 0.5 2013-03-03 04∶15∶51.95 0.4 2013-03-03 09∶26∶02.73 0.5 2013-03-04 22∶15∶58.27 0.4 2013-03-03 04∶32∶28.87 0.3 2013-03-03 10∶02∶14.13 0.3 2013-03-05 01∶05∶40.55 0.3 2013-03-03 05∶01∶31.96 0.2 2013-03-03 10∶30∶50.82 0.5 2013-03-05 01∶26∶19.18 0.7 2013-03-03 05∶01∶58.40 0.2 2013-03-03 12∶14∶59.36 0.6 2013-03-05 01∶28∶26.01 0.3 2013-03-03 06∶09∶27.16 0.5 2013-03-03 13∶29∶48.89 0.5 2013-03-05 01∶36∶13.80 0.2 2013-03-03 06∶41∶12.50 0.4 2013-03-05 01∶36∶30.16 0.2 2013-03-05 03∶47∶29.22 0.4 2013-03-05 07∶30∶04.27 0.4

表 2 遗漏地震事件发震时刻与震级 Table 2 The origin time and magnitude of detected events

图 4展示了4条发震时刻在3月3日8 ∶ 28至8 ∶ 30之间的遗漏地震事件滤波后的波形.由于发震时刻相距很近，使得不同地震间的波形相互交叠，且XBZ、LAY、LBG和GUY四个台站的记录受到稍早前一个塌陷事件波形的影响，台网编目人员难以通过人工分析有效辨别遗漏地震事件.

图 4

Fig. 4

图 4 4个遗漏地震事件的波形示意图不同颜色代表不同事件的波形，波形经过4阶零相位Butterworth滤波器2~8 Hz滤波，振幅经过归一化.顶部的箭头位置表示发震时刻.Fig. 4 Waveforms of four detected events in four different colorsAll waveforms are filtered in 2~8 Hz by 4th order zero-phase Butterworth filter， and the amplitudes are normalized.The locations of arrow on the top show the origin times of the four events.

本文同样利用匹配滤波技术估计了地震目录中所给出地震事件的发震时刻和震级.由于网络故障导致天津台网原始连续波形受到影响，造成数据不完整，3月4日下午的3次地震未检测到，目录中给 出的其余42次地震均检测到.图 5为估计结果与地震目录中给出参数的对比统计，发震时刻差别大部分在0.3 s以内，震级差别基本在0.3级以内，显示出本文对 遗漏地震事件发震时刻和震级的估计与观测报告的结果基本一致.个别事件发震时刻差距较大，可能与 其与模板地震震中位置有差别，P波走时不一致有关.

图 5

Fig. 5

图 5 利用本文方法估计地震目录中42个地震的发震时刻和震级结果，与目录给出参数的对比统计图a)发震时刻对比，横坐标为目录中给出的发震时刻(TCata)与本文估计所得发震时刻(TMFT)之差；(b)震级对比.Fig. 5 Statistical comparison of origin time and magnitude of the events in catalog between the result detected by matched filter technical and obtained from catalog (a)Comparison of origin time，x-axis is origin time obtained from catalog(TCata)minus the one detected in this paper(TMFT);(b)Comparison of magnitude.

图 6为补充遗漏地震前后地震目录的震级-频度关系对比，震级间隔取0.1.本文检测到的遗漏地震事件最大震级为M_L0.9，因而在≥M_L1.0地震目录没有变化.结果显示，增加了遗漏地震后，M_L0.3以上震级-频度关系与未加入遗漏地震相比呈现更好的线性特征，表明M_L0.3~0.8之间地震目录的完整性有较明显的改善.

图 7给出目录记录事件和遗漏地震事件的精定位结果.平均到时残差由0.1696 s降为0.0674 s.精定位前后震中分布对比表明，精定位后的震中分布呈现更加明显的条带状分布，显示此次震群发震构造走向可能为NW—SE向.图 7b剖面AA′上的震源深度投影表现出近垂直的条带状分布，推测发震构造可能为倾角较大的断裂.遗漏事件的加入强化了震中位置分布的条带状特征，对发震构造的推测和分析工作起到了积极作用.

如图 8所示，北西走向的施庄断裂在盆地以南出露较清晰(图 8a中F₈实线部分)，地表形成陡立的断层崖，多处可见壮观的巨型断层面(图 8c).高精度浅地震反射剖面结果(刘保金等，1999)推断其 向北西方向隐伏延伸(图 8a中F₈虚线部分).涿鹿

图 6

Fig. 6

图 6 补充遗漏地震前后地震目录的震级-频度关系对比图纵坐标中N为≥M(对应横坐标值)的地震个数.十字表示加入遗漏事件，空心圆圈表示未加入遗漏事件. Fig. 6 Comparison of magnitude-frequency relationship betweenetected+catalog events and catalog events only The parameter N in Y-axis is the event number ≥M(magnitude in X-axis). Cross means detected+catalog event， and the open circle means catalog events.

震群震中位置及震中分布与施庄断裂位置较一致，且精定位结果显示发震断层为倾角较大的断裂，施庄断裂地表出露部分近直立的特征相似，因而推测此次震群的发震构造可能为北西走向的施庄断裂.

图 7

Fig. 7

图 7 涿鹿震群地震精定位结果图a)震中分布；(b)震源在AA′剖面上的投影.黑色十字为震相报告中给出的震中位置和震源深度，红色和蓝色实心圆圈分别表示目录事件和遗漏事件的精定位结果.Fig. 7 Precise relocation result of Zhuolu earthquake swarm(a)Epicenter distribution in map view;(b)Epicenter distribution in profile AA′. Black cross represents epicenter of catalog event，red and blue solid circle represents epicenter of relocated catalog events and relocated detected events respectively.

图 8

Fig. 8

图 8 区域断裂构造分布略图(a)中F1 怀涿次级盆地北缘断裂，F2 宣化盆地南缘断裂，F3 延矾次级盆地北缘断裂，F4 蔚广盆地南缘断裂，F5 孙庄子—乌龙沟断裂，F6 紫荆关断裂，F7 南口山前断裂，F8 施庄断裂，F9 黄土窑—土木断裂，F10 新保安—沙城断裂，F11 桑干河断裂.灰色小圆点为此次涿鹿震群精定位后的震中位置，YH-1为浅地震反射剖面.(b)中FJ 方家冲西八里断裂，SZ 施庄断裂，LS 狼山断裂；S表示沉积层，UC表示上地壳，MC 为中地壳.改自文献(刘保金等，1999).(c)施庄断裂地表出露断层面.Fig. 8 Regional faults construct outline(a)F1 is north margin fault of Huailai-Zhuolu subbasin，F2 is south margin fault of Xuanhua basin，F3 is north margin fault of Yanqing-Fanshan subbasin，F4 is south margin fault of Weiguang basin，F5 is Sunzhuangzi-Wulonggou fault，F6 is Zijingguan fault，F7 is Nankou piedmont fault，F8 is Shizhuang fault，F9 is Huangtuyao-Tumu fault，F10 is Xinbaoan-Shacheng fault，F11 is Sanggan River fault. Gray points are the earthquake swarm，YH-1 is the shallow seismic reflection profile.(b)FJ is Fangjiachong Xibali fault，SZ is Shizhuang fault，LS is Langshan fault; S is sediments，UC is upper crust，MC is middle crust. Modified from reference(Liu et al., 1999).(c)Fault plane outcropping of Shizhuang fault.

本文利用匹配滤波技术对2013年3月河北涿鹿微震群中目录遗漏的地震进行了检测拾取，通过检测发现了52条遗漏地震事件，并给出其震中位置和震级估计结果.通过与地震目录中给出的参数对比，认为本文方法对遗漏地震事件的发震时刻和震级的估计是基本可信的.震级-频次统计分析表明加入遗漏事件后，地震目录的完整性在M_L0.3~0.8范围内有了较明显的改善.通过对涿鹿微震群精定位结果的分析，推测此次震群的发震构造为走向NW—SE向的近直立断层，施庄断裂为发震构造的可能性较大.

本文所使用方法一个显著的特点是在利用三分量波形互相关检测遗漏地震后，进一步使用P、S波波段波形互相关标定遗漏事件的P、S到时，从而能够对遗漏事件做精定位研究.前人对中强地震余震序列遗漏地震的研究中，多采用将遗漏地震震中位置直接置于模板地震震中的方法(Peng et al., 2007；Peng and Zhao， 2009；Meng et al., 2012).此方法简单易行，但会造成加入遗漏事件后的地震目录中多个地震震源位置完全重合，且遗漏地震的加入对地震序列空间分布没有任何影响，即无法对发震构造分析等后续研究做出贡献.本文在前人研究基础上，进一步采用震相波形互相关对遗漏事件进行精定位，使得遗漏地震的拾取对发震构造的推测和分析工作起到了积极作用.

匹配滤波技术中互相关扫描计算提供了震相互相关系数和标定时差，作为震相报告中到时数据的补充带入双差定位方法计算，能够在一定程度上克服测震台网人工拾取地震震相带来的定位误差，进一步提高地震精定位的精度和可信度.

利用匹配滤波技术拾取微震震群遗漏事件，能够改善震群目录完整性，对其他研究者利用震群目录的进一步深入研究工作有促进作用.与人工识别地震信号相比，匹配滤波技术能够有效抑制背景噪声、较大地震面波和震荡等造成的低频信号干扰，从而检测出较多的遗漏事件.Schaff(2010)对比了利用波形互相关和基于STA/LTA两种识别地震的方法，认为当地震事件之间波形相似度较高时，互相关识别方法可以识别出更多的遗漏地震事件.另一方面，匹配滤波技术需要地震模板、走时资料等先验信息，且互相关扫描计算时间较长，因而至今较难应用于测震台网的实时分析工作中.随着数字化波形资料的积累和计算技术的不断发展，匹配滤波技术的应用范围将会不断扩展.

不可否认，由于震级较小的地震波形信噪比低，且存在不与任何模板地震波形高度互相关的遗漏地震，本文方法难以检测到全部目录遗漏的地震.将地震精定位结果与震源机制相结合可以更加深入全面地研究震群的发震构造，是下一步研究的方向.

感谢两位审稿专家宝贵的建设性意见，感谢中国科学院地质与地球物理研究所陈棋福研究员的指导和建议，感谢中国地震局地震预测研究所王伟君副研究员、李乐副研究员、杨峰博士和熊仁伟等有益的讨论，感谢天津市地震局地震应急信息中心为本研究工作提供计算系统的支持.本文部分图件采用GMT软件包绘图.