吴梦羽 1 毛淑娟 2 宁杰远 1,† 唐启家 3 蒋一然 1
1.北京大学地球与空间科学学院, 北京 100871; 2. Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge 02139; 3.中国地质大学(武汉)地球物理与空间信息学院, 武汉 430074
摘要 在模板识别和地震相对定位方法的基础上, 提出一种探测断层结构的地震学方法。借用为了发现更多相似地震而在模板识别方法基础上提出的弱模板识别方法, 探测与模板地震有空间位置差异的地震, 然后根据波形互相关得到的到时差确定新识别出的地震相对于模板地震的空间分布。由于该方法对相似度要求较低, 且因计算量小而可以扫描更长时间的地震记录, 所以可以探测到更多能够比较精准地刻画断层细结构的微地震。利用渤海湾地区40天的地震波数据, 初步证明了这一方法的可行性。
关键词 模板识别方法(MFT); 断层结构; 地震定位; 渤海湾; 微地震
对断层精细结构的探测既是地震物理预测的基础 [1] , 也是统计预测的基础 [2 - 3] 。构造地质学家在断层的浅部构造研究方面做了卓有成效的工作, 对地震预测工作有很大帮助 [4 - 5] 。但是, 这些研究限于地表, 难以对地震高发的 10 km 左右深度断层的精细结构进行确定。地震定位方法在很大程度上解决了这一问题, 能够给出断层的深部结构, 有助于加深对断裂带动力学的认识 [1,6 - 7] 。为了刻画断层的精细结构, 需要识别出更多(绝对或相对)定位准确的地震事件。
近年来, 模板识别技术成功地用于识别因信噪比低而容易被人工识别遗漏的事件 [8 - 15] , 可以有效地检测到比普通地震目录中记录的事件低一个震级的新事件 [9 - 10,16] 。虽然传统的模板识别技术在研究地震活动性方面发挥了重要作用 [8,16] , 但由于很少关注目标地震与已知模板地震空间位置的不同, 所以尚未用于探测断层的精细结构。在模板识别技术的基础上, Zhang等 [17] 利用格点搜索法, 寻求在各台站的观测波形与模板地震波形互相关函数叠加出最大值的目标地震空间位置, 可以多识别9%的地震。
本文借鉴Shelly等 [8] 为了识别更多的低频地震(LFE)而提出的弱识别方法(weak detection), 放大单个台站单个分量上的互相关函数峰值在叠加中的影响范围, 从一个时间点放大到一个时间段, 从而加强目标地震与模板地震相似度不高时的识别效果, 并利用目标地震与模板地震震相的到时差, 对目标地震进行定位, 实现基于模板识别技术有效地探测断层精细结构的目标。我们将这一方法应用到渤海湾地区的微地震检测和定位中, 初步验证了方法的可行性。
模板识别技术(matched filter technique, MFT)广泛应用于研究触发地震 [12] 、非火山震颤 [8] 、地震监测 [13] 等地球物理学相关问题。最初, 模板识别技术是为寻找相似地震而提出的 [8,16,18] 。由于发震位置接近, 震源机制解接近, 能量接近, 因此同一观测台站记录的地震波形应该有一定程度的相似性 [19] , 所以该方法是将已知地震(模板地震, template)的波形与同一地震台站的连续波形记录进行互相关计算, 并对不同台站的互相关结果进行匹配叠加, 找出与已知波形相似的地震事件(图1)。
传统模板识别技术假设目标地震与模板地震在位置上一致, 因而采用模板的到时信息作为对齐的基准。对于很多发生在模板地震附近的地震, 尽管在处理单个台站资料时能得到较高的互相关值, 但在叠加过程中由于没有对齐而不能将目标地震有效地识别出来(图2(a))。
弱模板识别(Weak MFT, WMFT)方法极大地降低了对目标地震与模板地震之间波形相似程度的要求 [8] 。该方法将单一台站互相关极大值附近一段时间内的互相关值也设为极大值(扩展)。将极大值扩展后, 即使目标地震与模版地震的震相时差不一致, 不同台站之间的互相关极大值仍能通过叠加得到增强。这样, 原先在叠加中不能得到加强的信号就易于加强而被识别出来, 从而找到与模板不在同一位置, 但震相到时相似的目标地震(图2(b))。弱模板识别方法的计算方案包含模板识别的互相关计算, 但匹配叠加时利用的是修改后的互相关结果。同时, 存储的是匹配叠加超过阈值的结果以及与之相对应的短时间窗内未经修改的互相关结果(用于相对定位)。这样, 既能多识别地震, 又可以节约计算资源。
为了测试弱模板识别方法的地震检测能力和我们据此编制的计算程序的正确性, 首先重复Tang 等 [20] 在台湾地区的工作。将文献[20]中15071631. P10 作为模板, 对有观测资料的 11 个台站 2010 年3月7日的连续数据进行扫描。分别采用模板识别方法和弱模板识别方法进行匹配叠加。利用弱模板识别方法, 除找到所有模板识别方法的地震外, 新发现 4 个可能的地震事件(表1)。新发现的地震事件的观测波形及与模板地震波形的对比见图3, 可以看出, 弱模板识别技术可以有效地提升模板的识别能力。有时, 弱模板识别技术识别出的地震与模板地震的相似度不高, 这在寻找重复地震时是要避免的, 但在探测断层的精细结构时却求之不得。不仅如此,弱模板识别技术的地震识别能力也明显强于Zhang等 [17] 的方法。他们的方法能多识别9%的地震, 而弱模板识别方法能多识别55%的地震。Zhang等 [17] 的方法虽然允许目标地震与模板地震的空间位置有所偏离, 但由于计算量的限制, 格点搜索的网格不能太小, 导致大量事件被遗漏。
(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
新识别出的地震只有少数台站有观测记录, 震相不清晰, 但根据已经得到的波形互相关结果, 可以得到其震相与模板地震震相的精确到时差。由于同一模板地震能识别出很多地震, 这些地震和模板地震经常发生于同一断层, 因此这些地震相对于模板地震的空间分布可以用来刻画断层的几何形态。
表1 模板识别与弱模板识别方法的地震识别能力对比(以台湾地区的一个模板为例)
Table 1 Comparison of effectivity between MFT and WMFT(from study case of Taiwan)
说明: 加粗数字为弱模板方法多识别出的 4 个新地震事件。
本文据此提出设想: 利用互相关结果, 对新识别的地震相对于模板地震进行定位, 然后根据同一模板地震识别出的地震相对于模板地震的空间构形来研究断层的精细结构。
对于浅源地震, 在距离较近的台站, P 波与 S波的走时与台站至地震的距离基本上呈线性关系。因此有
灰色曲线为4个新地震事件到时前后的波形记录, 黑色曲线为模板地震记录
图3 台湾地区新发现地震的波形
Fig. 3 Waveforms of newly detected events in Taiwan area
(1) 其中, i 代表不同的震相种类; 
和 
为地震的类 S波和类 P 波的走时; 
为地震与台站的距离; 
和 
为修正参数, 只取决于震相类别; 
为常数, 也只取决于震相类别。
对于模板地震E和新找到的地震事件E′, 容易推导出
(2) 其中, 
, 
和 
为新地震事件 E′ 对应的参数。式(2)进一步变换为
。(3)
当台站距离模板地震较近时, 我们选取的主要是上行S波和P波, 此时 
; 当台站距离较远时, 和 值很小。在这两种条件下均有
。 (4) 由于 
及 
可以由互相关精确计算, 根据式(4), 我们就可以对新事件的位置进行初步定位。实际使用中, 以最速降线法求上述方程的最小二乘解, 即寻找新事件位置 
, 满足在所有台站上理论预测与对应的实际观测量之间的残差平方和最小:
。 (5)
我们对此定位方法进行理论测试。按照上述台湾地区11个台站设定台站位置, 模板位置为所有台站的中心点, 并以此为坐标原点, 深度设为10km。对3个新识别的地震位置, 我们采用TauP程序 [21] 计算所有需要的走时信息, 然后进行重新定位。重定位结果基本上与实际位置相符, 结果见表 2。
为了检验本文提出方法的可行性, 我们在渤海湾地区应用该技术进行地震的识别和相对定位。由于这一地区有大面积海域及近海沉积平原, 台站分布非常稀疏, 地震记录的信噪比很低, 所以需要利用长时间的观测数据才能得到相对准确的构造学结论。限于现有的计算设备, 本文仅使用短时间的观测数据来检验方法的可行性。
选用国家地震台网在渤海湾附近的固定台站2011年1月1日至 5月31日的连续记录进行测试, 数据采样率为50Hz。寻找新地震的目标时间段为2011年日本 3.11地震前后各 20 天, 其余数据用于模板事件波形的截取(图 4)。
采用文献[22‒23]中的方式, 对数据进行去平均值、滤波等预处理。由于我们的互相关运算只在相同台站的不同时间段之间进行, 各个台站仪器响应的不同不会影响结果, 因此未对数据进行去仪器响应处理, 以避免反卷积带来的不确定性。
在处理渤海湾地区的资料时, 我们借鉴 Wang 等 [12] 的经验对地震波能量进行分析, 以便确定滤波参数。随机选取部分模板地震进行时频分析(time-frequency analysis), 发现高频部分能量仍然很高, 10Hz 以上部分的能量不可忽略; 噪音部分在15~20Hz 之间也存在能量集中的现象(图5)。为了尽可能地保留模板地震的主要能量, 同时压制噪音能量, 选取15Hz 作为滤波的上限频率。在挑选到时的过程中, 我们发现部分台站存在明显的频率接近 1 Hz的背景噪音, 为了避免这些噪音对计算结果的干扰, 将滤波频率下限设为 2 Hz。
模板地震来自国家地震台网中心提供的中国地震台网目录。我们以 TauP 程序计算的理论到时为基准, 在邻近时间段内手动选取到时信息。没有直接采用理论到时的原因主要有两个: 1)理论速度模型与实际结构有差别, 理论到时可能不准确; 2)需要了解数据质量, 并剔除大量由于地震震级过小、距离过远而没有记录到有效数据的台站。
在各个台站的垂直分量记录中选取“P波”到时, 在两个水平分量记录中选取“S 波”到时 [24] , 并以这些到时所在分量上附近的波形记录作为模板波形。需要注意的是, 这里选取的“S 波”与“P 波”不一定是直达波, 也有可能是 Pn 或 sS 一类的转换波。虽然因新地震也会存在这些波形, 识别地震时并不关心具体震相, 但在地震相对定位时要进行区分。
在选择模板波形长度方面, 由于用途不同, 不同的研究之间存在较大的差别。Yao等 [25] 采用到时前2 s至后10 s的窗长, Slinkard等 [26] 采用Lg波到时前5 s至后15 s的数据作为模板波形。本文主要关注体波震相, 所以与 Tang 等 [23] 相同, 截取地震波到时前后各2 s, 共4 s时间窗内的波形记录。
图4 渤海地区台站以及所使用的模板地震分布
Fig. 4 Distribution of stations and templates used in Bohai Bay area
表2 相对定位测试结果
Table 2 Test of relative location
上图为高通滤波(1 Hz)后的地震归一化速度波形, 下图为时频分析得到的能量分布结果
图5 模板地震波形与时频分析结果
Fig. 5 Waveform of template event and the result of time-frequency analyze
尽管 Junek 等 [27] 认为 MFT 对新地震的探测能力可以达到模板地震事件附近数十公里的级别, 但为了以减小噪音叠加形成高于阈值假象的可能性,我们参照文献[8], 仍将极大值延展的窗长度选取为0.4s, 相当于只探测 3~5km 范围内的新地震事件。
在选取阈值方面, 不同研究者根据不同的数据、地震类型以及研究目的, 给出各种各样的判据 [16,28] 。我们按照 Tang 等 [23] 的研究, 初步选取叠加序列的6 倍标准差作为阈值。
为使互相关函数满足高斯分布的统计特性, 保障所选阈值的有效性, 我们仅对单道互相关系数较大、信噪比较高的部分进行延展。
为了确定互相关系数与信噪比的关系, 我们选取一段模板地震记录叠加上滤波后的随机噪音, 再计算其与原先波形的互相关系数。在噪音振幅与地震信号相当, 地震在刚刚可以被识别出的临界情况下, 互相关系数为 0.4~0.5。因此, 我们只对互相关系数大于0.45的部分进行延展, 确保被延展信号的信噪比。
我们采用 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个事件, 明显多于传统方法。由于二者阈值选取相同, 因此可以看出改进后的识别方法能更加有效地识别出潜在的地震。
表3 弱模板识别与传统模板识别的效果对比
Table 3 Comparison between effect of MWFT and MFT
为了提高用本文方法找到的新事件的可信度, 尽可能去除噪音部分对模板识别的干扰, 我们在 6倍标准差的基础上设定进一步的阈值标准: 1)阈值标准设定为 0.27 的 CC 值与 6 倍标准差中的较大者, 对大部分事件近似于至少 10 倍标准差; 2)对每个识别出的事件, 舍弃部分 CC 值较低的台站以降低其负面效果, 但仍要求具有 3 个以上的可用台站。按照上述阈值标准, 共得到 478 个可能的潜在事件(图 7(c))。
灰色曲线为新地震到时前后的波形记录, 黑色曲线为模板地震记录
图6 模板地震与新地震波形对比
Fig. 6 Waveforms of the template event and the newly detected event
去除定位可信度较低的事件后, 得到所有互相关系数大于 0.27 的事件的定位结果(图 8(a))。作为对比, 我们同时采用 TomoDD 方法进行定位(图 8 (b))。可以看出, 在地震分布较分散的情况下, 与TomoDD 相比, 本文的相对定位方法能够对更多的地震进行定位, 适用性更强。本文着重于介绍方法, 要真正研究该地区断层的精细结构, 后续工作中需要更长时间的观测资料。但是, 从图 8(a)已经可以看出张渤带地震条带细结构与整个地震条带走向不一致的特征以及郯庐断裂带有与其走向斜交的断层细结构, 这些断层细结构特征与地表断层的分布特征是一致的。
(a)所有识别出的地震的互相关系数分布; (b)去除台站少于3个的事件; (c)去除部分互相关系数较低的台站, 并提高阈值到0.27
图7 新识别地震互相关系数
Fig. 7 CC value distribution of new detected events
(a) 本文不依赖于速度结构的相对定位方法的定位结果; (b) TomoDD方法的定位结果。圆点为模板地震, 五角星为识别出的新地震
图8 新地震的空间分布
Fig. 8 Distribution of new detected events
本文利用弱模板识别方法, 对与模板地震有空间位置差异的未知地震进行探测。探测到这些与模板地震相似度较低的地震事件后, 进一步利用由波形互相关得出的到时差, 可以确定新识别出的地震事件相对于模板地震的空间分布。这种方法能够探测到相似度不高的地震, 且由于计算量相对较小, 能够扫描更长时间的地震记录, 因此可以探测到更多能比较精准地刻画断层细结构的地震, 从而更好地把握断层附近应力场的演化趋势, 提升统计地震学的预测能力。
在渤海湾地区的检验结果初步证明了本文方法的可行性, 可以检测到更多的地震事件, 并得到合理的定位结果。即将开始的《中国地震科学台阵观测——华北地区东部》野外观测项目将在该地区布设大量的观测台站。未来, 利用该计划提供的高质量数据, 有望使用本文方法研究渤海湾地区断层的精细结构。
本文方法也可用于研究其他地区断层的精细结构。利用本文方法, 有望揭示全新的全球断层精细结构图像, 为研究全球地球动力学演化及地震预测提供亟需的基础资料。
致谢 中国地震局地球物理研究所国家数字测震台网数据备份中心提供地震波形数据, 在此表示感谢。
参考文献
[1]Cai C, Yu CQ, Tao K, et al. Spatial distribution and focal mechanism solutions of the Wenchuan earth-quake series: results and implications. Earthquake Science, 2011, 24(1): 115‒126
[2]Jiang M M, Zhou S Y, Chen Y S, et al. A new mul-tidimensional stress release statistical model based on coseismic stress transfer, Geophys J Int, 2011, 187: 1479‒1494
[3]贾科, 周仕勇. 西南地区1900—1970年历史地震震源参数推断及结果的不确定性分析. 地球物理学报, 2012, 55(9): 2948‒2962
[4]Xu C, Xu X, Shyu J B H. Database and spatial distribution of landslides triggered by the Lushan, China Mw 6.6 earthquake of 20 April 2013. Geomor-phology, 2015, 248: 77‒92
[5]Korup O. Geomorphic implications of fault zone weakening: slope instability along the Alpine Fault, South Westland to Fiordland. New Zealand Journal of Geology and Geophysics, 2004, 47(2): 257‒267
[6]Lou X, Cai C, Yu C, et al. Intermediate-depth earth-quakes beneath the Pamir-Hindu Kush region: evi-dence for collision between two opposite subduction zones. Earthquake Science, 2009, 22(6): 659‒665
[7]Zhou S Y, Xu Z H, Han J, et al. Analysis on the master event method and precise location of 1997 Jiashi strong earthquake swarm in western China. Acta Seismologica Sinica, 1999, 12(3): 285‒291
[8]Shelly D R, Beroza G C, Ide S. Non-volcanic tremor and low-frequency earthquake swarms. Nature, 2007, 446: 305‒307
[9]Schaff D P. Broad-scale applicability of correlation detectors to China seismicity. Geophysical Research Letters, 2009, 36(11): 192‒200
[10]Schaff D P, Waldhauser F. One magnitude unit re-duction in detection threshold by cross correlation applied to Parkfield (California) and China seismi-city. Bulletin of the Seismological Society of America, 2010, 100(6): 3224‒3238
[11]Gonzalez-Huizar H, Velasco A A, Peng Z, et al. Re-mote triggered seismicity caused by the 2011, M9.0 Tohoku-Oki, Japan earthquake. Geophysical Research Letters, 2012, 39(10): 10302
[12]Wang W, Meng X, Peng Z, et al. Increasing back-ground seismicity and dynamic triggering behaviors with nearby mining activities around Fangshan Pluton in Beijing, China. Journal of Geophysical Research: Solid Earth, 2015, 120(8): 5624‒5638
[13]Zhang M, Wen L. An effective method for small event detection: match and locate (M&L). Geophy-sical Journal International, 2015, 200(3): 1523‒1537
[14]Junek W N, Kværna T, Pirli M, et al. Inferring aftershock sequence properties and tectonic structure using empirical signal detectors. Pure and Applied Geophysics, 2015, 172(2): 359‒373
[15]Slinkard M, Heck S, Schaff D, et al. Detection of the Wenchuan aftershock sequence using waveform corre-lation with a composite regional network. Bulletin of the Seismological Society of America, 2016, 106(4): 1371–1379
[16]Gibbons S J, Ringdal F. The detection of low mag-nitude seismic events using array-based waveform correlation. Geophysical Journal International, 2006, 165(1): 149‒166
[17]Zhang M, Wen L. An effective method for small event detection: match and locate (M&L). Geophy-sical Journal International, 2015, 200(3): 1523‒1537
[18]Zhang J, Richards P G, Schaff D P. Wide-scale de-tection of earthquake waveform doublets and further evidence for inner core super-rotation. Geophysical Journal International, 2008, 174(3): 993–1006
[19]Helmberger D V. Theory and application of synthetic seismograms. Earthquakes: Observation, Theory and Interpretation, 1983, 37: 173‒222
[20]Tang C C, Peng Z, Lin C H, et al. Statistical proper-ties of low-frequency earthquakes triggered by large earthquakes in southern Taiwan. Earth and Planetary Science Letters, 2013, 373: 1‒7
[21]Crotwell H P, Owens T J, Ritsema J. The TauP Toolkit: flexible seismic travel-time and ray-path utilities. Seismological Research Letters, 1999, 70(2): 154‒160
[22]Peng Z, Zhao P. Migration of early aftershocks fol-lowing the 2004 Parkfield earthquake. Nature Geo-science, 2009, 2(12): 877‒881
[23]Tang C C, Lin C H, Peng Z. Spatial-temporal evolu-tion of early aftershocks following the 2010 ML 6.4 Jiashian earthquake in southern Taiwan. Geophysical Journal International, 2014, 199(3): 1772‒1783
[24]Aso N, Ohta K, Ide S. Volcanic-like low-frequency earthquakes beneath Osaka Bay in the absence of a volcano. Geophysical Research Letters, 2011, 38(8): L08303
[25]Yao D, Peng Z, Meng X. Remotely triggered earth-quakes in South-Central Tibet following the 2004 Mw 9.1 Sumatra and 2005 Mw 8.6 Nias earthquakes. Geophysical Journal International, 2015, 201(2): 543‒551
[26]Slinkard M, Schaff D, Mikhailova N, et al. Multista-tion validation of waveform correlation techniques as applied to broad regional monitoring. Bulletin of the Seismological Society of America, 2014, 104(6): 2768‒ 2781
[27]Junek W N, Kværna T, Pirli M, et al. Inferring aftershock sequence properties and tectonic structure using empirical signal detectors. Pure and Applied Geophysics, 2015, 172(2): 359‒373
[28]Schaff D P, Richards P G. Improvements in magni-tude precision, using the statistics of relative amp-litudes measured by cross correlation. Geophysical Journal International, 2014, 197(1): 335‒350
Fault Structure Detection by Matched Filter Technique
WU Mengyu 1 , MAO Shujuan 2 , NING Jieyuan 1,† , TANG Qijia 3 , JIANG Yiran 1
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
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
doi: 10.13209/j.0479-8023.2017.099
中图分类号 P315
收稿日期: 2017-05-04;
修回日期: 2017-05-24;
网络出版日期: 2017-05-26
†通信作者 , E-mail: njy@pku.edu.cn
† Corresponding author , E-mail: njy@pku.edu.cn
内蒙古自治区 2016 年度科技重大专项“重点地区地震预测预警技术研究开发与推广示范”和中国地震科学台阵探测项目(DQJB16A0307)资助