利用模板匹配技术识别大岗山水库库区爆破事件

杜瑶 阮祥 戴仕贵 邵玉平 胥珍 余洋洋 罗勇

（中国成都 610041 四川省地震局）

0 引言

地震事件的识别与拾取在地震数据处理过程中非常重要，根据识别手段的不同，可划分为人工经验分析和计算机自动识别。随着地震台网的日益密集及专用地震台网的不断布设，对水库库区、页岩气开采区、石油钻探区等工程场地的地震监测逐步开展起来，监测数据越来越多，人工拾取地震事件的效率和精度已远不能满足地震数据处理的需求。同时，随着微小震级地震事件研究的不断深入，对如何高效率、高精度、低错误率地识别微小地震事件提出了更高的要求。在专业地震台网监测的水库库区及附近，地震类型较复杂，构造地震及人工地震活动等均可能会在监测区内出现。由于人工爆破波形与天然地震波形具有极高的相似性（尹欣欣等，2020），因此，如何有效地将不同类型的地震加以区分，从波形上分辨天然地震、人工爆破，提高地震定位、编目的质量，对于蓄水前后水库库区及附近地震活动研究，以及分析预报工作中正确判定人工爆破等异常都具有较重要的实际意义。

国内外针对地震类别的识别已开展了一些相关研究。张春贺等（2006）在小爆破识别与研究的基础上，通过周期—频度谱、波形特征等对爆破事件的自动识别进行了研究，自动识别概率达到了77%以上；崔鑫等（2016）通过对天然地震和人工爆破的频谱分析认为，天然地震和人工爆破的记录具有明显不同的时频特征；尹欣欣等（2020）利用波形互相关的模板匹配方法对爆破及天然地震事件进行震相识别，对于人工爆破，准确率为80%，对于天然地震，准确率为100%。在日常编目工作中发现，人工识别时震级较小的近震及爆破极易混淆，因此，建立一种能够及时准确识别天然地震和人工爆破的方法，把人工爆破与震级较小的地方震区别开来，以得到真实的区域地震活动变化很有必要（郑秀芬等，2006；包淑娴，2011）。

大岗山水电站是在四川境内大渡河干流近期开发的大型水电工程之一。电站正常蓄水位 1 130 m，最大坝高210 m，总库容7.42 亿m3，电站装机容量2 600 MW。电站于2005 年开工建设，2015 年建成。大岗山水库于2015 年5 月28 日导流底孔下闸蓄水，这标志着水库正式开始蓄水发电，该水库属典型的高山峡谷型高坝大库容水库。

大岗山水库库区位于青藏高原东南缘，地属川西高原区。库区所在区域地处以鲜水河—安宁河—小江断裂为界的川滇断块与大凉山断块结合部位SN 向构造带北段，该处为SN 向与NW、NE 向等多组构造的交汇复合部位，区域构造上位于由金坪断裂和磨西断裂所围限的黄草山断块的西侧边缘，库区内的主要构造形迹以近SN、NNW、NW 向的褶皱和断裂为主。主要断裂有大渡河断裂、磨西断裂、大发断裂、安宁河断裂等及规模较小的断层，与库水直接接触的岩性主要为花岗岩，占库区的大部分库段，在水库田湾支流有泥盆系、二叠系的一套碎屑岩出露，在库中段有三叠系碎屑岩出露。从工程地震条件来看，库区属于断裂块状岩体亚型，可能诱震强度为中等或强烈地震（《中国水力发电工程》编审委员会，2000）。区域内有较强构造地震背景，其中，1975 年1 月15 日大坝西边的康定6.2 级地震震中距坝址约31 km，是历史上距库区最近且影响最大的地震。与低烈度地区水库诱发地震不同，高烈度地区水库在蓄水前就存在明显地震活动（程万正，2013），因此，在水库蓄水后如何判断是否库水加卸载诱发的地震活动则需要开展更细致的工作。2014 年1 月、2015 年4 月、2015 年10 月等大岗山水库蓄水前后库区及附近均发生过大量微小地震，并且地震类型较复杂，有天然地震、工业爆破及蓄水后诱发地震等。

1 方法及数据

2012 年12 月大岗山水库地震监测台网软硬件设备安装完成，2013 年2 月完成调试运行，同年3 月开始考核运行，同时开始监测及观测数据的处理。库区共布设了8 个高增益短周期速度记录地震仪，数据采集器为EDAS-24GN3 型，地震计为RSFS-1A 型，仪器观测频带为2 s—40 Hz，布局孔径NS20 km×EW40 km，台站沿库区均匀展布，监测范围包括可能诱发地震的重点监视区段，监测震级下限为ML0.5，震中定位误差≤1 km（阮祥等，2017）。根据大岗山水库地震监测台网分布，并结合其外围地区的区域地质构造环境和历史地震情况，将大岗山水库库区地震活动监测研究区域定为29.2°—29.9°N、101.9°—102.5°E。

采用Zhang 等（2015）对微震同时进行检测及定位的方法，即匹配定位方法（Matching and location，简称M&L 方法），并在GeoTaos 软件平台下进行模板扫描与定位。该方法需要精确的模板地震事件，同时更注重微弱小震事件的拾取与定位。M&L 方法计算流程（图2）及原理如下 （Zhang et al，2015）。

图2 M&L 方法模板扫描流程Fig.2 Flow diagram of M&L method

图1 研究区域位置及地震台站分布阴影部分为研究区Fig.1 Map showing study area and stations

根据震群序列重新定位后将震级较大的地震作为模板事件，以模板事件位置为中心在经度、纬度、深度等3 个方向上进行网格化，计算每个可能的微震位置（格点）与模板事件位置之间的参考震相（S 波）在同一地震台站上的走时差，然后按此走时差对所有台站分量记录的模板事件参考震相（S 震相）与相应台站分量的连续数据进行滑动互相关，并计算互相关波形叠加后的平均相关系数CC和信噪比SNR。当平均相关系数和信噪比大于阈值时认为检测到了1 个微震，并把其位置确定在最大相关系数的格点位置（曾宪伟等，2017）。其中，由模板事件与检测出的微震事件之间的位置差异所导致的走时差计算公式如下（Wen，2006；Wen et al，2010；Zhang，2013）

其中，dDk为模板事件与可能的微震事件之间相对差异所导致的台站k震中距的差异；dh为2个事件间深度的相对变化；分别为地震相位P 的走时分别对模板震中距D（水平慢度）和模板深度h（垂直慢度）的导数，利用一维速度模型，分别计算每个台站及其相关地震相位的。走时差Δt(k,P)的计算对参考模型依赖较小。

研究中获取了部分爆破作业的准确时间数据，并根据这些资料在地震台网目录获取事件，从连续波形中挑选了信噪比高、震相清晰的5 次爆破事件（图3、表1）。利用这5 次已知的爆破事件作为模板事件，进行去均值、去线性趋势等预处理以及1—10 Hz的带通滤波，连续波形重采样率至50 Hz，因爆破事件样本及记录到爆破的台站数量较少，故直接使用观测报告中定位的震相文件，没有再进行重新定位。由于震中距越大的台站记录波形特征越不明显，近台震相波形特征明显，因此模板选取近台记录的爆破事件（图4）。

图3 研究区域爆破事件分布Fig.3 Distribution of epicenters of blasting events in the study area

图4 爆破模板事件波形（a）HCP 台2014-03-10 T 18:44:56.1；（b）CAK 台2014-03-12 T 18:18:36.4；（c）BLG 台2014-03-24 T 18:01:43.5；（d）CAK 台2014-04-09 T 18:10:39.5；（e）CAK 台2014-04-28 T 18:07:50.6Fig.4 Template waveforms of blasting events

表1 爆破事件统计Table 1 Statistics of blasting events

选取大岗山水库蓄水前后2014 年3 月1 日至2015 年12 月31 日库区地震活动监测研究区域（29.2°—29.9°N，101.9°—102.5°E）8 个台站（XIX、AJW、ZMC、NTW、HCP、BLG、CAK、XMC）三分量地震仪记录的连续波形，利用M&L 方法检测定位。

由于近震爆破与天然地震在波形上有差异，前者P 波能量强于S 波，更容易从噪声中检测出来，因此选取P 波及后5 s 的时间作为匹配相位和匹配窗口，水平、垂直搜索距离分别为2 km，水平搜索网格数为20，垂直搜索网格数为10，模板相位信噪比阈值为10，综合考虑匹配的爆破数量及质量，确定3 个及以上台站同时记录到时检测相关系数阈值为0.3，只有2 个台站记录到时为0.8，仅1 个台站记录到时为0.95，所采用的四川地区一维平均速度模型见表2。图5 分别选取了与模板事件2、4、5 同一台站（CAK 台）记录的事件波形匹配出的事件，2 个不同波形之间的相关系数分别为0.628 5、0.773 4、0.409 7。从图5 还可见，扫描出的事件波形与模板事件均有震中距为数十千米的小吨位爆破的特征：①纵波在垂直分向的初动向上；②垂直分向上呈振幅较大的脉冲型；③有较明显的反射和折射波等。同一台站记录的波形相似度极高，可确定为同类型事件。

图5 M&L 方法检测出的事件Fig.5 Comparison of the template seismograms (blakc traces) with signals detected in the continuous waveforms (red traces) by using M&L Method

表2 四川地区一维平均速度模型Table 2 One-dimensional average velocity model for Sichuan area

2 结果及分析

利用M&L 方法，剔除因断记等因素造成的相关系数＞1 的误匹配结果，再经过人工复合，最终获得23 次疑爆记录，其中，5 次疑爆事件与已知爆破事件波形位置吻合且相关系数为1.0，可以确定为爆破。图6 为爆破模板事件及扫描出的爆破（疑爆）事件分布。在图6 中划分了A、B、C、D、E 等5 个事件集中区域，其中，A 区域有疑爆事件12 次（含爆破1 次）；B 区域有疑爆事件2 次（含爆破1 次）；C 区域有爆破事件1 次；D 区域有疑爆事件4 次（含爆破1 次）；E 区域有疑爆事件4 次（含爆破1 次）。匹配结果为已知爆破模板事件个数的4.6 倍，有效地提高了爆破事件检出率，同时将匹配识别出的爆破事件形成目录与原分析编目报告中的事件进行比对检验，区别并剔除原报告中人工误判的爆破事件，提高报告准确性。将匹配出来的疑爆事件及距离最近的台站（CAK）记录到的NS 向波形分别在Google Earth 上投影（图7），由于疑爆事件发生时间均为2014—2015 年，而Google Earth 的影像则是2020 年8 月，所以存在部分爆破环境与当时实际环境不一致的可能，并且可能受速度模型的影响，事件定位也会有一点偏离，但偏离量很小，基本在0.5 km 之内。由图7 可见，扫描出的爆破（疑爆）事件附近存在采石场或蓄水前修建大坝的场地，这与爆破实际情况较符合，通过模板搜索出的爆破（疑爆）事件详见表3。

图6 爆破模板及扫描出的爆破（疑爆）事件分布Fig.6 Distribution of templates of blasting events and detected blasting or suspicious blasting events

图7 匹配出来的疑爆事件及CAK 台记录到的NS 向波形在Google Earth 上的投影Fig.7 Projection of detected blasting or suspicious blasting events and seismograms recorded by CAK station in NS direction at satellite map

表3 通过模板扫描出的疑爆及爆破事件Table 3 Catalogue of suspicious blasting and blasting events detected by template

对表3 中的发震时刻在2014 年3 月1 日至2015 年12 月31 日人工编目观测报告中进行搜索，找到并剔除了观测报告中人工误分的爆破及疑爆事件8 次（图8），表明可以利用M&L 模板匹配方法来有效区别不同类型的地震事件。

图8 人工编目观测报告中被误判为爆破事件的地震震相报告Fig.8 Seismic phase report of misjudged blasting event in the original observation reports

3 结论

选取大岗山水库库区及附近波形记录较好的5 次爆破作为模板事件，采用M&L 方法对2014 年3 月至2015 年12 月连续波形进行检测与定位，并进行人工复核，最终获得23 次爆破（或疑爆）事件，这些事件均呈现出震中距为数十千米范围内小吨位爆破的特征：①纵波在垂直分向的初动向上；②垂直分向上呈振幅较大的脉冲型；③有较明显的反射和折射波等。扫描出的爆破（或疑爆）事件与模板事件有较高的相关性，相关系数达0.773 4。将部分扫描出的事件投影至Google Earth，从影像图上可见附近存在采石场或蓄水前修建大坝的场地，这与爆破实际情况较符合，重定位偏离基本在0.5 km 之内。将识别出来的爆破（或疑爆）目录与人工编目观测报告进行对比，剔除了8 次人工误判的爆破事件，有效地提高了观测报告的精度。大岗山水库库区及附近在水库蓄水前存在大量的人工活动，而真实的地震活动背景才对库区蓄水前后地震活动的研究具有参考意义。微小地震与爆破在人工分析处理中容易混淆，而利用人工地震活动作为模板事件在连续波形中扫描定位的方法，可有效地剔除人工地震活动事件。建议建立非天然地震数据库，利用匹配定位等方法区分人工活动地震事件，在地震监测工作中复检工程地震台网监测报告，落实因人工活动形成的震群异常。此外，因爆破等人工活动地震事件震源位置精确，可将其用于反演速度结构等研究工作中。

日本产业技术综合研究所雷兴林研究员提供Geotaos 软件，并给予了指导与帮助，在此表示感谢。