Fault Structure Detection by Matched Filter Technique

WU Mengyu1, MAO Shujuan2, NING Jieyuan1,†, TANG Qijia3, JIANG Yiran1

1. School of Earth and Space Sciences, Peking University, Beijing 100871; 2. Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge 02139; 3. Institute of Geophysics and Geomaties, China University of Geoscience, Wuhan 430074; † Corresponding author, E-mail: njy@pku.edu.cn

Abstract Based on matched filter technique (MFT) and earthquake relative location technique, a new method is proposed to detect fault structure. The weak matched filter technique (WMFT) is modified from MFT to find more similar earthquakes, which is now used to detected events different from known templates in location. After new events are detected, their locations relative to templates can be given by using time differences from cross correlation. Due to low requirement of the WMFT in wave similarity and its low computational cost, more events can be detected and might help depict the structure of faults clearly. The feasibility of this method is primarily proved by the study of Bohai Bay area using records of 40 days. Key words matched filter technique (MFT); fault structure; earthquake location; Bohai Bay; micro seismicity

对断层精细结构的探测既是地震物理预测的基础[1],也是统计预测的基础[23]。构造地质学家在断层的浅部构造研究方面做了卓有成效的工作,对地[45]震预测工作有很大帮助。但是,这些研究限于地表,难以对地震高发的10 km左右深度断层的精细结构进行确定。地震定位方法在很大程度上解决了这一问题,能够给出断层的深部结构,有助于加深对断裂带动力学的认识[1,67]。为了刻画断层的精

细结构, 需要识别出更多(绝对或相对)定位准确的地震事件。近年来, 模板识别技术成功地用于识别因信噪比低而容易被人工识别遗漏的事件[815], 可以有效地检测到比普通地震目录中记录的事件低一个震级的新事件[910,16]。虽然传统的模板识别技术在研究地震活动性方面发挥了重要作用[8,16], 但由于很少关注目标地震与已知模板地震空间位置的不同, 所

以尚未用于探测断层的精细结构。在模板识别技术的基础上, Zhang等[17]利用格点搜索法, 寻求在各台站的观测波形与模板地震波形互相关函数叠加出最大值的目标地震空间位置, 可以多识别9%的地震。

[8]本文借鉴Shelly等 为了识别更多的低频地震(Lfe)而提出的弱识别方法(weak detection), 放大单个台站单个分量上的互相关函数峰值在叠加中的影响范围, 从一个时间点放大到一个时间段, 从而加强目标地震与模板地震相似度不高时的识别效果,并利用目标地震与模板地震震相的到时差, 对目标地震进行定位, 实现基于模板识别技术有效地探测断层精细结构的目标。我们将这一方法应用到渤海湾地区的微地震检测和定位中, 初步验证了方法的可行性。

1 方法介绍1.1 模板识别技术简介

模板识别技术(matched filter technique, MFT)广泛应用于研究触发地震[12]、非火山震颤[8]、地震监测[13]等地球物理学相关问题。最初, 模板识别技术是为寻找相似地震而提出的[8,16,18]。由于发震位置接近, 震源机制解接近, 能量接近, 因此同一观测台站记录的地震波形应该有一定程度的相似性[19],所以该方法是将已知地震(模板地震, template)的波形与同一地震台站的连续波形记录进行互相关计算, 并对不同台站的互相关结果进行匹配叠加, 找 出与已知波形相似的地震事件(图 1)。

1.2 弱模板识别技术

传统模板识别技术假设目标地震与模板地震在位置上一致, 因而采用模板的到时信息作为对齐的基准。对于很多发生在模板地震附近的地震, 尽管在处理单个台站资料时能得到较高的互相关值, 但在叠加过程中由于没有对齐而不能将目标地震有效地识别出来(图 2(a))。

弱模板识别(Weak MFT, WMFT)方法极大地降低了对目标地震与模板地震之间波形相似程度的要求[8]。该方法将单一台站互相关极大值附近一段时间内的互相关值也设为极大值(扩展)。将极大值扩展后, 即使目标地震与模版地震的震相时差不一致,不同台站之间的互相关极大值仍能通过叠加得到增强。这样, 原先在叠加中不能得到加强的信号就易于加强而被识别出来, 从而找到与模板不在同一位置, 但震相到时相似的目标地震(图2(b))。弱模板识别方法的计算方案包含模板识别的互相关计算,但匹配叠加时利用的是修改后的互相关结果。同时, 存储的是匹配叠加超过阈值的结果以及与之相对应的短时间窗内未经修改的互相关结果(用于相对定位)。这样, 既能多识别地震, 又可以节约计算资源。

1.3 两种模板识别技术地震识别能力的对比

为了测试弱模板识别方法的地震检测能力和我

们据此编制的计算程序的正确性, 首先重复Tang [20]等在台湾地区的工作。将文献[20]中15071631. P10 作为模板, 对有观测资料的 11 个台站2010年3月7日的连续数据进行扫描。分别采用模板识别方法和弱模板识别方法进行匹配叠加。利用弱模板识别方法,除找到所有模板识别方法的地震外,新发现4个可能的地震事件(表1)。新发现的地震事件的观测波形及与模板地震波形的对比见图3,可以看出,弱模板识别技术可以有效地提升模板的识别能力。有时,弱模板识别技术识别出的地震与模板地震的相似度不高,这在寻找重复地震时是要避免的,但在探测断层的精细结构时却求之不得。不仅如此,弱模板识别技术的地震识别能力也明显强于Zhang [17]等的方法。他们的方法能多识别9%的地震, 而[17]弱模板识别方法能多识别55%的地震。Zhang等的方法虽然允许目标地震与模板地震的空间位置有所偏离,但由于计算量的限制,格点搜索的网格不能太小,导致大量事件被遗漏。

1.4 对新识别地震的相对定位

新识别出的地震只有少数台站有观测记录, 震相不清晰, 但根据已经得到的波形互相关结果, 可以得到其震相与模板地震震相的精确到时差。由于同一模板地震能识别出很多地震, 这些地震和模板地震经常发生于同一断层, 因此这些地震相对于模板地震的空间分布可以用来刻画断层的几何形态。 本文据此提出设想: 利用互相关结果, 对新识别的地震相对于模板地震进行定位, 然后根据同一模板地震识别出的地震相对于模板地震的空间构形来研究断层的精细结构。

对于浅源地震, 在距离较近的台站, P 波与 S波的走时与台站至地震的距离基本上呈线性关系。因此有

当台站距离模板地震较近时, 我们选取的主要是上行 S波和 P 波, 此时∆ Ti  ∆Li  0; 当台站距离较远时, ∆Ti 和 ∆Li值很小。在这两种条件下均有 由于T S  T S及T P  T P 可以由互相关精确计算,根据式(4), 我们就可以对新事件的位置进行初步定位。实际使用中, 以最速降线法求上述方程的最小二乘解, 即寻找新事件位置( x  , y  , z) , 满足在所有台站上理论预测与对应的实际观测量之间的残差平方和最小:

我们对此定位方法进行理论测试。按照上述台湾地区11个台站设定台站位置, 模板位置为所有台站的中心点, 并以此为坐标原点, 深度设为10 km。对 3个新识别的地震位置, 我们采用Taup程序[21]计算所有需要的走时信息, 然后进行重新定位。重定位结果基本上与实际位置相符, 结果见表2。

2 利用渤海湾地区的观测结果检验新方法

为了检验本文提出方法的可行性, 我们在渤海湾地区应用该技术进行地震的识别和相对定位。由于这一地区有大面积海域及近海沉积平原, 台站分布非常稀疏, 地震记录的信噪比很低, 所以需要利用长时间的观测数据才能得到相对准确的构造学结论。限于现有的计算设备, 本文仅使用短时间的观测数据来检验方法的可行性。

2.1 数据处理

选用国家地震台网在渤海湾附近的固定台站2011年1月1日至5月31日的连续记录进行测试,数据采样率为50 Hz。寻找新地震的目标时间段为2011年日本3.11地震前后各 20 天, 其余数据用于采用文献[22‒23]中的方式,模板事件波形的截取(图4)。

对数据进行去平均值、滤波等预处理。由于我们的互相关运算只在相同台站的不同时间段之间进行,各个台站仪器响应的不同不会影响结果,因此未对数据进行去仪器响应处理,以避免反卷积带来的不确定性。

在处理渤海湾地区的资料时, 我们借鉴Wang等[12]的经验对地震波能量进行分析,以便确定滤波参数。随机选取部分模板地震进行时频分析(timefrequency analysis), 发现高频部分能量仍然很高, 10 Hz 以上部分的能量不可忽略; 噪音部分在15~20 Hz 之间也存在能量集中的现象(图5)。为了尽可能地保留模板地震的主要能量,同时压制噪音能量,选取15 Hz 作为滤波的上限频率。在挑选到时的过 程中, 我们发现部分台站存在明显的频率接近1 Hz的背景噪音, 为了避免这些噪音对计算结果的干扰,将滤波频率下限设为 2 Hz。

2.2 模板获取

模板地震来自国家地震台网中心提供的中国地震台网目录。我们以 Taup 程序计算的理论到时为基准, 在邻近时间段内手动选取到时信息。没有直接采用理论到时的原因主要有两个: 1) 理论速度模型与实际结构有差别, 理论到时可能不准确; 2) 需要了解数据质量, 并剔除大量由于地震震级过小、距离过远而没有记录到有效数据的台站。

在各个台站的垂直分量记录中选取“P 波”到时,在两个水平分量记录中选取“S 波”到时[24], 并以这些到时所在分量上附近的波形记录作为模板波形。需要注意的是, 这里选取的“S波”与“P 波”不一定是直达波, 也有可能是 Pn 或 ss 一类的转换波。虽然因新地震也会存在这些波形, 识别地震时并不关心具体震相, 但在地震相对定位时要进行区分。

在选择模板波形长度方面, 由于用途不同, 不同的研究之间存在较大的差别。Yao 等[25]采用到时

前 2s 至后 10 s 的窗长, Slinkard 等[26]采用 Lg波到时前5 s 至后 15 s的数据作为模板波形。本文主要[23]关注体波震相,所以与 Tang 等 相同,截取地震波到时前后各2 s, 共 4s时间窗内的波形记录。

2.3 匹配叠加参数的选取

尽管 Junek 等[27]认为 MFT 对新地震的探测能力可以达到模板地震事件附近数十公里的级别, 但为了以减小噪音叠加形成高于阈值假象的可能性,我们参照文献[8], 仍将极大值延展的窗长度选取为0.4 s, 相当于只探测3~5 km 范围内的新地震事件。在选取阈值方面, 不同研究者根据不同的数据、地震类型以及研究目的, 给出各种各样的判据[16,28]。我们按照 Tang 等[23]的研究, 初步选取叠加序列的6倍标准差作为阈值。

为使互相关函数满足高斯分布的统计特性, 保障所选阈值的有效性, 我们仅对单道互相关系数较大、信噪比较高的部分进行延展。

为了确定互相关系数与信噪比的关系, 我们选取一段模板地震记录叠加上滤波后的随机噪音, 再计算其与原先波形的互相关系数。在噪音振幅与地震信号相当, 地震在刚刚可以被识别出的临界情况 下, 互相关系数为 0.4~0.5。因此, 我们只对互相关系数大于 0.45 的部分进行延展, 确保被延展信号的信噪比。

2.4 弱模板识别方法的效果

我们采用 2011 年 1 月 1 日至 5 月 31 日共 277个模板地震, 分别用 MFT 和 WMFT 方法, 对 2011年 3 月 11 日日本大地震前后各 20 天内可能存在的地震事件进行检测, 结果见表 3。可以看出, 弱模板识别方法的识别能力大大优于传统模板识别方法。通过直接的波形对比(图 6), 可以证实本文提出的方法是有效的。将弱模板识别方法识别出的新事件按互相关值从高到低排列, 并对比我们使用的阈值(图 7(a))。部分新事件仅对应 1~2 个台站, 数量过少, 无法完成定位等工作, 因此将这一部分舍弃(图 7(b))。从整体上看, 阈值的分布较为稳定, 表明渤海地区噪音与地震波互相关的6倍标准差基本上处在CC 值0.15~0.25 之间。去除少于 3个台站参与识别的结果后, 改进的模板识别方法共得到 2969 个事件, 明显多于传统方法。由于二者阈值选取相同, 因此可以看出改进后的识别方法能更加有效地识别出潜在

的地震。为了提高用本文方法找到的新事件的可信度,尽可能去除噪音部分对模板识别的干扰, 我们在6倍标准差的基础上设定进一步的阈值标准: 1) 阈值标准设定为 0.27 的 CC 值与 6 倍标准差中的较大者, 对大部分事件近似于至少 10 倍标准差; 2) 对每个识别出的事件, 舍弃部分 CC 值较低的台站以降低其负面效果, 但仍要求具有 3 个以上的可用台站。按照上述阈值标准, 共得到 478 个可能的潜在事件(图7(c))。

去除定位可信度较低的事件后, 得到所有互相关系数大于 0.27 的事件的定位结果(图 8(a))。作为对比, 我们同时采用 Tomodd 方法进行定位(图 8 (b))。可以看出, 在地震分布较分散的情况下, 与

(a) 上图为单一台站的单道互相关结果, CC (cross correlation coefficient)表示归一化的互相关系数; 下图的黑色曲线为已知地震模板的归一化速度波形, 灰色曲线为目标的归一化速度波形, 模板窗口的中心对应上图中最大互相关系数出现的位置, 横轴的具体时间由模板地震与目标地震的记录时刻差决定。(b) 灰色曲线为 3 个台站(CHN1, SGS, WTP)东西方向的互相关函数,黑色曲线为按模板到时对齐后 3个台站东西分量的互相关函数以及对齐的互相关函数叠加后结果 图 1采用模版识别技术探测地震的实例Fig. 1 An example of earthquake detection using MFT

(a) 同图 1, 按照模板到时不能有效地对齐不同台站的互相关信号; (b) 按照弱模板识别方式, 将(a)中 3个分量的记录进行扩展后再叠加图 2 传统模板识别与弱模板识别的比较Fig. 2 Comparison between traditional MFT and the WMFT

灰色曲线为 4个新地震事件到时前后的波形记录, 黑色曲线为模板地震记录图 3台湾地区新发现地震的波形Fig. 3 Waveforms of newly detected events in Taiwan area

图 4渤海地区台站以及所使用的模板地震分布Fig. 4 Distribution of stations and templates used in Bohai Bay area

上图为高通滤波(1 Hz)后的地震归一化速度波形, 下图为时频分析得到的能量分布结果图 5模板地震波形与时频分析结果Fig. 5 Waveform of template event and the result of time-frequency analyze

(a) 所有识别出的地震的互相关系数分布; (b) 去除台站少于 3 个的事件; (c) 去除部分互相关系数较低的台站, 并提高阈值到 0.27图 7新识别地震互相关系数Fig. 7 CC value distribution of new detected events

灰色曲线为新地震到时前后的波形记录, 黑色曲线为模板地震记录图 6模板地震与新地震波形对比Fig. 6 Waveforms of the template event and the newly detected event