上海及邻区微震检测与定位

张雅楠， 李红谊*， 张盛中， 李炎臻， 黄雅芬， 钟卫星

1 中国地质大学(北京)地球物理与信息技术学院， 北京 100083 2 上海佘山地球物理国家野外科学观测研究站， 上海 201602 3 中国地质大学(北京)信息网络中心， 北京 100083 4 上海市地震局佘山地震监测中心站， 上海 201602

0 引言

地震目录作为地震危险性分析和地震活动性研究的关键基础资料，其完整性以及地震事件定位的精准性会对后续如地震活动性分析等研究的可靠性产生直接影响(朱艾斓等，2005；王同利等，2019；Ross et al.，2019).由于微小地震能量弱，信噪比低，因此日常的地震目录中，对微小地震的记录往往会出现比较严重的遗漏.对微震活动的检测，不仅可以用来间接监测大的自然灾害，如火山喷发、诱发地震等(Kato et al.，2012；Kiser and Ishii，2013；高保彬等，2014；余国锋等，2018)；同时对微小地震时空分布特征的研究，能够为我们深入认知大地震的发震机理和断层带结构等(Mori and Hartzell，1990；Wu et al.，2017)提供宝贵的信息.

微震检测过程中，由于中强地震后的尾波干扰、台站的分布不合理、以及背景噪声较强等原因，微弱的地震信号很难被准确的识别，导致微震事件的遗漏.本文拟利用GPU-M&L技术(Liu et al.，2020a)对上海及附近地区进行微震检测与识别，以期提升地震目录的完备性.现有的地震事件识别方法中，较常用到的有两种方法(Zhang and Wen，2015；Grigoli et al.，2018)：一种是基于波形的地震事件检测方法(Kao and Shan，2004；Shelly et al.，2006；Drew et al.，2013)，另一种方法是基于震相拾取的地震事件检测方法(Allen，1982；Withers et al.，1998).其中，机器学习算法被广泛应用于震相拾取工作(Ross et al.，2018；于子叶等，2018；Wang et al.，2019；Zhu and Beroza，2019；Liu et al.，2020b)，并表现出高效率、高精度的优势，在地震目录构建中展现出了很好的应用前景(Huang et al.，2020；Wang et al.，2020；Tan et al.，2021)；基于波形互相关的匹配滤波技术近十年被广泛应用于微小地震的检测，该方法源于滑动互相关检测技术，是一种在数据信噪比较低的情况下仍能提取到微弱信号的有效方法.

地震记录到的信号存在各种噪声的干扰，当地震信号能量较弱，或者与噪声信号量级相当时，我们通常很难从噪声信号里识别地震信号.因此，对地震资料去噪是数据处理工作中的一项重要的内容.本文采用深度神经网络构建的DeepDenoiser方法(Zhu et al.，2018)来进行去噪.该方法在有效地剔除环境噪声对记录信号干扰的同时，可有效地将地震信号与噪声分离，并且在低信噪比的条件下对不同形式的信号有较好的去噪效果.

地震在时与空维度上的分布可充分发掘出活动断层深部构造及运动特征，同时地震精定位已被广泛应用在研究活动构造和隐伏断层过程中.地震震相拾取的精准程度对后续事件定位处理有着极大地影响，在地震数据处理过程中扮演着重要的角色，其中包含了信号的检测与到时估计和震相识别三步骤.本文采用深度学习算法PhaseNet(Zhu and Beroza，2019)进行P、S震相自动识别和到时拾取，之后基于机器学习获取的震相拾取结果，利用双差定位法(Waldhauser and Ellsworth，2000)对检测到的事件进行精定位.该方法消除了震源至地震台站共同的传播路径效应，有效减少了由于对地壳结构了解不精细而引起的定位误差.

本文的研究区域位于30.3°N—32°N，119.4°E—122.35°E，在海域一侧，台站分布较少(黄耘等，2008).同时，由于该地区地处中国人口众多的长三角地区，背景环境噪声较大，许多微弱的地震信号埋没在背景噪声中，不易被距离较远的台站发现，若依照传统地震事件检测方法，诸多微小地震大概率将发生遗漏.而上海地区位于扬子地块与华夏地块之间的钦杭结合带北东向延伸处，从全球构造的角度，该地区位于太平洋板块与欧亚板块相互作用的影响范围，是受到板块构造运动及板内构造变形影响的地区之一(Waldhauser，2009；于海英等，2021).研究区内的昆山—嘉定断裂沿江苏唯亭、昆山南、蓬朗及上海外岗一线呈北东东—东西向延伸，江苏境内唯亭—蓬朗段长约40 km，为安角断凹的北界(顾澎涛和王尧舜，1993；张浩等，2021).安角断凹为一长期继承性发育的箕状断陷盆地，由于昆山—嘉定断裂的持续活动，致使盆地沉积中心向北迁移.历史上断凹内地震频繁，1731年昆山以南的5级地震可能与该断裂有关(王卫平，1997).

本文收集了中国数字化地震台网所提供的上海地区及邻区13个台站记录的2011—2020年共10年的连续地震资料.首先基于已有的地震目录利用GPU-M&L方法进行地震检测与识别，然后通过DeepDenoiser方法对检测到的地震事件进行去噪处理，进一步确认检测到的地震事件，并利用PhaseNet进行震相拾取，最后使用双差定位构建高分辨地震目录，为后续研究发震断层构造等提供了数据支撑.

1 数据与方法

1.1 台站和数据

本研究资料来源为上海佘山地球物理国家野外科学观测研究站提供的中国数字化地震台网连续数据和地震目录，资料范围为30.3°N—32°N，119.4°E—122.35°E，涉及上海及其邻近区域等共13个台站，2011—2020年之间共10年的连续记录地震资料.图1给出了中国地震台网中心目录(简称台网目录)给出的研究区域内146个地震事件以及地震台站与断裂带的分布.

图1 研究区内的台站与断层分布图蓝色三角形代表台站，黑色虚线代表断层；黄色方框代表安角断凹；黄色字体JS、ZJ、SH分别代表该区域位于江苏省，浙江省，上海市；F1，南通—上海断裂；F2，昆山—嘉定断裂；F3，千灯—黄渡断裂；F4，枫泾—川沙断裂带；F5，湖州—嘉善—平湖断裂；F6，天马山断裂.Fig.1 Map of stations and faults in the study area The blue triangles represent stations, the black dashed lines represent fault; the yellow square denotes Anjiao sunken; the yellow fonts JS, ZJ and SH represent the area located in Jiangsu Province, Zhejiang Province, and Shanghai City, respectively; F1, Nantong-Shanghai fault; F2, Kunshan-Jiading fault; F3, Qiandeng-Huangdu fault; F4, Fengjing-Chuansha fault zone; F5, Huzhou-Jiashan-Pinghu fault; F6, Tianma mountain fault.

1.2 基于图像加速器改进的微震检测与定位

传统模板匹配技术(Gibbons and Ringdal，2006；Shelly et al.，2006；Peng and Zhao，2009;李璐等，2017)通过对模板地震波形与连续波形进行互相关叠加，来实现对微小地震的检测.但该方法的实施仅针对紧邻参考模板的地震且无法对地震事件进行定位.Zhang和Wen(2015)通过在相关波形叠加前考虑微震与模板之间的位置差，开发了一种对微震事件同时进行检测和定位的方法，该方法可以检测到更小震级的地震事件以及那些距离模板事件较远的微震.但由于该方法需要计算模板与微震每个可能位置之间的走时差，以实现对模板周围的三维空间进行搜索，因此运行起来非常耗时.为了提高匹配定位的计算效率和探测能力，Liu等(2020a)发展了GPU-M&L方法，该方法针对不同的台站分量记录，综合考虑几何扩散和地震波形的信噪比来赋予不同的权重系数，在利用GPU加速的同时，提高地震事件的探测能力.具体处理步骤如下：

(1)连续波形预处理：首先对连续波形进行降采样，采样率从0.01降到0.05，在保证扫描效果良好的情况下，减小运行内存且大幅提升运行速度；然后对连续波形做2～8 Hz带通滤波；设置每天的连续波形数据的零时刻作为参考，扫描长度取决于扫描时连续波形记录中的最短时间长度.

(2)模板事件挑选原则：1)模板事件的三分量能够被所有的台站记录到；2)地震事件有清晰的P波、S波震相，信号无缺失，且信噪比(SNR)大于3.0.依据以上挑选原则，我们从台网目录中选取了136个地震事件作为扫描的模板事件，利用PhaseNet标记模板事件波形到时，模板扫描的时间窗口为模板波形S波到时的前1 s至后5 s，扫描后计算平均滑动互相关系数.

(3)地震事件的检测：通过参考已有的研究以及对本文中数据的分析后，我们设定当平均互相关值超过9倍的绝对偏差中值(Median Absolute Deviation， MAD)或者探测事件的平均滑动互相关系数(CC)大于预设的经验探测阈值(经验阈值一般为CC=0.19)时，则将该事件认定为一个潜在地震事件.

我们将检测得到的地震事件目录称为检测目录，图2展示了一个检测到的地震事件的实例，模板事件为2012年10月31日的ML0.8地震事件，检测到的事件为2012年12月19日ML0.64地震事件.

图2 利用GPU-M&L检测到的一个地震事件样例(a) 黑色实线：检测标准为MAD的9倍，红点：检测到的事件； (b) 对应(a)所示窗的CC值分布直方图；(c) 模板波形(红色)与检测到的地震事件波形(灰色)比较.Fig.2 An example of seismic event detected by GPU-M&L method(a) The black solid line indicates that the detection standard is 9 times MAD, and the red dot represents the detected event; (b) Histogram of CC value distribution corresponding to the window shown in (a); (c) Comparison of the template waveform (red) and the detected seismic event waveform (gray).

1.3 深度去噪

以往的研究通常为利用模板进行微震检测后，通过波形包络线等方法对检测到的地震事件进行简单的辨识(Liu et al., 2019).由于检测到的微震事件波形的低信噪比，通常很大一部分检测到的地震事件震相信息很难提取，从而影响微震事件的精定位.为了进一步对GPU-M&L技术检测到的地震事件进行确认，本文采用DeepDenoiser来进一步辨认探测到的地震事件.将带有噪声的地震事件数据输入到已经训练好的神经网络模型中，输出两个单独的模块，一个用于记录地震信号时频特征，另一个用于记录地震信号波形特征.我们将GPU-M&L方法检测到的事件数据作为输入，传输到神经网络模型中，该模型把输入数据变换到相应的时间域和频率域.频率域模块给出了信号与噪声的时频分布，即不同时间和频率的能量强度，通过时频分析，我们可以清晰的将信号与噪声分离出来；时间域模块则分别展示了分离后的信号与噪声波形.

1.4 震相拾取与双差定位

本文采用PhaseNet来对地震事件的P波和S波到时进行拾取，通过将原始波形输入到卷积神经网络中，逐层自动地从三分量地震记录中识别P波与S波震相.HypoDD基于地震对的震相相对到时差对地震事件进行初始相对定位校正，从而获得精度较高的地震事件相对位置.因此，利用HypoDD相对定位过程中，精准的震相识取这个步骤最为关键(赵明等，2021).本研究中，我们将利用机器学习拾取到的震相作为双差定位的输入.在定位过程中，考虑P波震相走时精度高于S波震相走时，因此我们赋予P波震相1.0的权重，S波0.5的权重.此外，我们设定事件对与台站之间的最大距离为300 km；事件可以组成的最大事件对为10个；事件对之间的间距小于40 km；地震对的最大距离为25 km，一个地震事件最多与25个最近的地震事件组成地震对.

2 结果分析

2.1 检测地震事件

完整的地震目录有助于对地震活动性的准确分析.本文利用从台网目录中挑选的136个高信噪比的地震事件作为模板(图3)，采用GPU-M&L技术对研究区10年连续地震资料进行扫描，共检测到824个地震事件，约为台网目录提供地震数目的5.5倍.通过检测事件与地震目录数量对比(图9b)，可以看出大多数遗漏的事件震级偏小，以ML2.0以下的微震为主.

图3 利用GPU-M&L检测模板事件分布图红色实心圆表示模板事件，实心圆的大小代表事件的震级大小.Fig.3 Distribution of template events used by GPU-M&L The red dots represent template events, and the size of the solid circles represent the magnitude.

2.2 检测事件去噪处理

我们利用DeepDenoiser对扫描的地震事件进行去噪处理以进一步确认扫描到的地震事件.图4展示了我们检测到的一个ML0.5的地震事件的时频分布特征图，以及分离后时间域的地震信号与噪声信号.图4a2展示了利用DeepDenoiser去噪后的地震信号频谱特征，图4b2给出了去噪后相应于图4a2的地震信号的时间序列.可以看到，经过DeepDenoiser去噪后可以提取到清晰的地震波形信息，不仅显著提高了信噪比，同时还为我们后续的震相拾取工作提供了良好的数据基础.通过DeepDenoiser去噪，我们最终从824个检测到的地震事件中确认了333个高信噪比的地震事件.

图4 利用DeepDenoiser分离ML0.5的事件波形噪声与信号(a) 时频域； (b) 时间域. (a1, b1) 带有噪声的信号； (a2, b2) 剔除噪声后的地震信号； (a3, b3) 剔除的噪声信号.Fig.4 Separation of noise and earthquake signal from a ML0.5 event by DeepDenoiser(a) Time-frequency domain; (b) Time domain. (a1, b1) Raw data; (a2, b2) Earthquake signal extracted from raw data; (a3, b3) Noise extracted from raw data.

图5给出了台网目录记录的2011年2月25日ML0.3地震经过2～8 Hz带通滤波后的地震事件波形和经过DeepDenoiser去噪后提取的地震波形记录.从图中可以看到，经过DeepDenoiser去噪，地震信号的信噪比大幅度提高了，同时也为我们后续的震相拾取提供了高质量的波形数据.

图5 利用DeepDenoiser对一个地震事件Z分量去噪前后的波形对比(a) 去噪前的地震波形； (b) 去噪后的地震波形.图中T1, T2分别代表根据理论模型计算得到的S波和P波到时.Fig.5 Comparison of vertical component waveform of an earthquake before and after DeepDenoiser The seismic event waveform (a) before and (b) after DeepDenoiser. T1 and T2 represent the theoretical S wave and P wave arrival time, respectively.

2.3 P波、S波震相拾取

利用基于机器学习的震相拾取技术，我们对333个检测到的地震事件和136个模板地震事件进行了震相拾取，最后分别得到了13个台站记录到的2540和2122个高质量的P波和S波到时数据(如图6所示).

图6 定位后P波，S波走时分布图红色实点与黑色实点分别代表P波与S波.Fig.6 Travel time distribution of P wave and S wave after relocationRed dots and black dots represent P wave and S wave, respectively.

2.4 双差定位后的地震分布特征

高质量的震相到时数据为高精度的地震定位提供了重要的保障.本文参考温燕林等(2021)的上海地区速度结构模型进行精定位，速度模型结构如表1所示.图7给出了检测到的2012年12月19日一个ML0.73地震事件在双差定位前后的震中位置和震相到时数据对比.从图中可以清晰的看到，通过双差精定位后，图7c给出了更合理的按震中距排列的P波与S波震相到时.总体来说，上海地区地震活动性弱，小震分布较稀疏，邻近的地震事件波形有着较高的相似性，通过采用HypoDD我们对479个事件进行精定位，其中可重定位的地震数为412个，重定位后RMS的平均值为0.157 s.图8给出了精定位后地震事件的震中分布图，平均定位误差在E-W，N-S，U-D方向分别为0.035，0.044，0.044 km.

表1 上海地区速度结构模型Table 1 Velocity structure model in Shanghai area

图7 地震事件定位前后震相到时对比示例图(a) 红色实心圆代表地震目录事件，白色五角星代表定位前的事件位置，紫色五角星代表定位后的事件位置； (b) 定位前按震中距排列的震相到时; (c) 定位后按震中距排列的震相到时，蓝色短线与红色短线分别代表P波与S波到时.Fig.7 An event example of phase arrivals comparisions before and after relocation(a) The red dots represent earthquakes given by the routine catalog, and the white and purple stars represent events before and after relocation, respectively. The seismic event before relocation (b) and after relocation (c) waveform are sorted by epicenter distances; Blue short line and red short line mark P and S wave arrival time respectively.

图8 重定位后的地震事件震中分布Fig.8 The distribution of epicenters of seismic events after relocation

精定位结果显示，上海地区及邻区的地震事件分布相对较为分散，以小震活动为主，其中：(1)部分地震沿枫泾—川沙断裂分布，该断裂形成较早，燕山期断裂重新活动；(2)部分地震事件集中分布于安角断凹.安角断凹以千灯—黄渡断裂和昆山—嘉定断裂为界，受周边断裂构造控制呈北东东向展布，为中新生代隐伏断陷盆地.历史上断凹内地震事件相对较频繁.我们检测到的地震事件在该隐伏断陷盆地也较为集中，主要分布于近东西走向的昆山—嘉定断裂与千灯—黄渡断裂之间，推测可能在这两条断裂之间存在一条北北东走向的隐伏断层；(3)部分小震事件沿北西向南通—上海断裂呈线性分布；(4)通过对太湖地区的微震事件发生时间进行分析，推测可能是由于2018年夏天建设太湖隧道时的爆破或引起的微震.

2.5 检测事件的震级-频度分布

Mc(Magnitude of completeness)为地震目录的最小完备震级，是地震台网监测能力的重要参数之一，通常将其定义为所检测能力覆盖范围内的所有地震事件中检测到的最小震级.

图9给出了本研究台网目录、检测事件与去噪后得到的地震事件震级与频度关系并分别给出了其不同震级频段对应的事件数量对比.从图中可以清晰的看到，经过模板检测及去噪后，最小完备震级由台网目录的Mc1.0降为最终检测目录的Mc0.8，有效提升了地震目录在小震级范围的完备性.

图9 原始地震目录、检测及去噪后地震目录的震级-频率关系图与震级-数量分布图(a) 蓝色空心圆、橙色空心圆和灰色空心圆分别代表台网目录、去噪后目录与检测目录；实心三角形代表相应目录的最小完备震级；(b) 蓝色柱体、橙色柱体和灰色柱体分别代表台网目录、去噪后目录与检测目录事件位于不同震级区间的数量分布.Fig.9 Magnitude-frequency and magnitude-number diagrams of local catalog, catalog after detection and catalog after denosing(a) The blue circle, orange circle and grey circle represent the China Earthquake Network Center (CENC) catalogue, denoised catalogue and detection catalogue, respectively; the black triangle represents the minimum complete magnitude of the corresponding catalogue； (b) The blue bar, orange bar and grey bar mean the number of events in different magnitude within CENC catalogue, denoised catalogue and detection catalogue, respectively.

3 讨论与结论

GPU-M&L方法通过计算模板波形和连续地震信号之间波形的相关性，将其互相关进行叠加，并根据信噪比给不同的台站分量赋予不同的权重系数，显著提高微震事件探测能力；同时采取GPU内核进行加速，提高了检测的速度(Liu et al.，2020a).图10为一个模板事件与检测到的多个地震事件之间的互相关系数随震级与距离变化图，由图中可以看到，对于与模板事件间距离越小且震级越大的检测事件，对应的互相关系数越大，即事件波形的相似度越高.但由于该方法对微震事件的检测基于阈值统计的局限性，尽管可以检测到更多的微震事件，但同时有可能产生如误检等病态检测.

图10 模板事件与检测事件的互相关系数、震级与距离关系图Fig.10 The relationship of cross-correlation coefficients, magnitude and distance between the template event and detected events

DeepDenoiser去噪与PhaseNet震相拾取均采取深度学习对地震信号进行模拟与学习，经过对美国北加州地区海量地震数据的训练，Deepdonoiser所提炼的模型同时从时间和相位信息中学习，能有效提取出地震信号；Phasenet采用迁移训练的模型对到时的拾取也更为精准.这两类模型所采用的深度卷积神经网络利用Relu非线性激活函数，有助于提高网络的泛化能力.PhaseNet已经被广泛应用于意大利及中国江苏、云南等地区地震序列的快速检测与精定位的研究，并取得良好的效果(Huang et al.，2020；Cianetti et al.，2021；廖诗荣等，2021)，证实此类模型有着较好的泛化性.

我们首先采用GPU-M&L技术进行微震事件的检测，然后利用DeepDenoiser方法对检测到的事件进行进一步识别，将识别到的微震事件通过震相拾取算法PhaseNet标记震相到时，最后采用HypoDD算法对微震事件进行进一步精定位.我们将这套流程首次应用于上海地区的微震检测，取得了显著效果，也验证该套方法的技术可行性.尽管上海地区地震活动性较弱，近十年记录到的地震事件有限且分布稀疏，我们仍然检测到了比原来地震目录要多出2倍多的地震事件.

本文主要关注微震的检测和定位，对震级的估算并不是关注的重点，因此我们简单比对了去噪前后估算的震级，去噪后的地震事件震级较去噪前整体偏小.采用去噪前的信号计算震级可能会由于噪声将地震信号淹没，导致计算的过程中实际振幅偏大，使得震级计算结果相应偏大.

综上所述，本文利用GPU-M&L技术，从中国地震台网中心提供的146个地震事件目录中挑选了136个地震事件作为模板事件，对上海及邻区13个台站连续10年的地震记录进行扫描，共识别出824个地震事件.然后我们通过深度去噪方法将信号与噪声分离，对去噪后的记录的频率和振幅特性进行分析，从而进一步对识别出的地震事件进行确认.检测和去噪后的地震目录的完备震级由台网目录的Mc1.0降为Mc0.8，基于此，我们利用PhaseNet对去噪后的333个地震事件和136个模板事件进行了震相的拾取.最后基于机器学习拾取的震相数据，我们利用双差定位法对479个地震事件进行了精定位.精定位后的结果显示：上海地区整体地震活动性较弱，地震事件空间分布比较零星，以小震为主，重定位后地震事件部分集中于安角断凹，部分沿北东向枫泾—川沙隐伏断裂分布，部分沿北西向的南通—上海断裂呈线性分布.本研究构建的从GPU-M&L检测与定位、DeepDenoiser深度去噪、PhaseNet震相拾取到HypoDD双差定位技术路线也为其他地区开展更高效、高质的微震检测与定位提供了一套可参考的流程.我们的微震检测结果也为研究上海及邻区地震风险性评估、地震灾害和地震活动性等提供了重要的数据支撑.

致谢感谢上海市地震局为本文提供数据支持.感谢香港中文大学刘敏博士开源的GPU-M&L代码(GitHub-MinLiu19/GPU-MatchLocate1.0)，加州理工学院朱尉强博士的DeepDenoiser及PhaseNet代码(GitHub-wayneweiqiang/DeepDenoiser，GitHub-wayneweiqiang/PhaseNet).