引言

受太平洋板块、菲律宾海板块和欧亚板块的俯冲和碰撞作用，我国台湾及台湾海峡地区的地震活动较为强烈。受台湾及台湾海峡强烈的构造运动影响，福建地区构造活动较为复杂。该地区自西而东存在三条NE走向的地震构造带，分别为邵武—河源、政和—海丰和长乐—诏安地震构造带（徐锡伟等，2016；郭晓然等，2019），与之共轭的三条断裂带自南向北分别为上杭—云霄、永安—晋江和闽江断裂带（李祖宁等，2014）（图1）。福建地区地质构造较为复杂，历史上曾发生过多次M7.0以上地震，例如1600年广东南澳附近的M7.0地震、1604年12月29日福建泉州外海M7${\raise0.5ex\hbox{$\scriptstyle 1$}\kern-0.1em/\kern-0.15em\lower0.25ex\hbox{$\scriptstyle 2$}} $地震、1918年2月13日南澳附近M7.3地震（国家地震局震害防御司，1995；中国地震局震害防御司，1999）。M5.0以上破坏性地震主要分布在福建西南地区。福建地区及台湾海峡现今中小地震也较为频繁，仅2015年福建地震台网记录到的编目地震就达974条（除单台记录外），震级主要集中在M_L0.0—2.9，地震主要分布在各断裂带上。因此，研究福建地区断裂带的地震活动性，预测未来中强地震可能发生的区域，是防震减灾工作的重要内容。

图  1  福建地震台网台站、历史中强地震、2015年中小地震分布和区域断裂带分布图

F₁：滨海断裂带；F₂：长乐—诏安断裂带；F₃：政和—海丰断裂带；F₄：邵武—河源断裂带；F₅：闽江断裂带；F₆：永安—晋江断裂带；F₇：上杭—云霄断裂带，下同

Figure  1.  The distribution of seismic stations，history moderate and strong earthquakes ，mid-small earthquakes in 2015 and fault zones in Fujian area

F₁：Binhai fault zone；F₂：Changle-Zhao’an fault zone；F₃：Zhenghe-Haifeng fault zone；F₄：Shaowu-Heyuan fault zone； F₅：Minjiang fault zone；F₆：Yong’an-Jinjiang fault zone；F₇：Shanghang-Yunxiao fault zone，the same below

下载: 全尺寸图片 幻灯片

目前针对福建省和台湾海峡地质构造和地震活动性的研究十分广泛，包括：利用天然地震Pg和P震相走时联合反演地壳上地幔顶部结构（李细兵等，2019a）；震源机制解的应力分析（钟继茂，2009）；气枪激发震源的一维P波速度结构（李细兵等，2019b）和三维P波层析成像（金震等，2018）等。虽然研究者通过人工激发地震的方式获取研究区域的精细速度结构，但精确的三维速度结构模型在目前的地震定位中应用不够广泛；而人工震相识别受人为因素影响存在系统偏差，且受背景噪声干扰，从而导致地震定位精度不高、地震目录完备性欠缺，不能满足现今断裂活动性及地震活动性研究的需要。因此，高精度地震事件自动检测技术成为目前主流的研究方向。

近年来，机器学习和匹配滤波技术在获取较完整地震目录方面取得了很好的成效。随着机器学习的突破性发展，基于数据驱动的机器学习微震检测识别技术极大地扩充了地震目录数量，在国内外地震活动性、断裂精细结构及活动习性、诱发地震等研究中得到了广泛应用（Zhang et al，2019；Zhu et al，2019；Liu et al，2020）。但机器学习需要从大量训练样本中提取统计特征，除对地震样本数有一定要求外，还对计算机性能有较高要求。匹配滤波方法（matched filter technique）（Turin，1960；Gibbons，Ringdal，2006）是通过波形互相关寻找与已知模板事件相似的地震事件，可在信噪比相对较低的条件下获取更为完备的目录（Shelly et al，2007），但其检测能力仍然受限于模板事件的分布和数量，且运行效率较低。

为了规避匹配滤波方法在地震自动检测中的局限性，在计算机配置有限的条件下，我们采用基于台网约束下的长短窗特征触发式震相检测的模板事件识别方法和匹配滤波自动检测方法（phase picking，association，location，matched filter，缩写为PALM）（Zhou et al，2019）对2015年福建地震台网连续波形进行扫描和地震精定位，以期获得断裂带的地震活动性。

1.   方法和数据

1.1   方法

PALM方法是一种基于匹配滤波技术的自动目录扩展方法，利用STA/LTA （Allen，1978）拾取震相，在台网的约束下，通过震相关联算法去除单台误识别，获得完整的模板事件，再由基于匹配滤波技术的目录扩展方法，利用高性价比GPU计算地震模板事件与原始连续波形的互相关函数，提高地震识别率，尤其是在福建地区这类背景噪音相对较高的区域，可有效扩展地震目录数目，为后续地震构造背景研究提供可靠数据。PALM方法在构建2021年漾濞M_S6.4高精度地震序列时（Zhou et al，2021）与Liu等（2022）用类似方法获取的地震条目在同一数量级上，且地震分布具有较好的一致性，该方法在2019年7月4日南加州里奇克莱斯特（Ridgecrest）地震序列（Zhou et al，2022）的构建上亦取得了良好效果。PALM方法由两部分构成：PAL模板拾取和MSMS匹配滤波。PAL模板拾取部分分为三个步骤：震相拾取（phase picking）、台网约束下的震相关联（association）及地震事件的定位（location）。首先震相拾取基于STA/LTA特征的触发式震相检测和寻峰式震相拾取算法（Allen，1978，1982），利用STA/LTA长短窗能量比检测Z分量的P波峰值，再通过偏振滤波器在水平分量上将P波能量进行压制后利用STA/LTA正确识别出S波峰值，该方法可以准确地给出P波和S波震相到时。其次，震相关联算法基于台站记录在时间和空间上的关联，地震事件在各台站的记录需满足：① 在时间域，由各台站的P波和S波到时估算出的发震时刻与实际发震时刻相近；② 在空间域，台站的震中距与P波和S波到时差估算出的震中距相近。时间域的关联通过搜索发震时刻丛集实现，而空间域的关联则通过网格搜索走时残差最小的震中位置实现。最后，利用HypoInverse （Klein，2002）定位法进行地震事件定位，得到PAL模板事件目录。

MSMS匹配滤波部分是一种可在GPU上进行计算的自动匹配滤波技术方法。该方法共分为四步：① 计算事件模板与待检测连续波形的互相关函数；② 根据模板时间走时将各台站的互相关函数对齐；③ 用蒙版的方式拓展互相关函数峰值的范围，由于模板事件与被检测事件的震中位置差异、模板事件震相拾取的误差以及介质波速变化等因素，仅通过走时偏移无法完全对齐互相关峰值，造成后续叠加时相消，而通过应用峰值蒙版可以减少此类偏差，从而提高模板的检测能力；④ 将各台用蒙版拓展过的互相关函数进行线性叠加（Zhou et al，2019），多台叠加可有效压制噪声的互相关值，从而降低误识别率。基于互相关的震相拾取可以准确地确定相对走时：

式中：t_corr为互相关峰值位置确定的震相到时；t_org为被检测事件的发震时刻；t_temp为模板事件在该台站的走时。根据模板事件与被检测事件的相对走时，即可使用双差重定位算法（HypoDD）确定事件间的相对位置（Waldhauser，Ellsworth，2000）。基于PALM方法（Zhou，2022），我们可以构建从原始波形到获取精定位地震事件的整体架构。在HypoInverse绝对定位和HypoDD相对定位中，我们均使用一维福建2015模型（朱元清等，2017）作为初始速度模型。图2分别为以2015年1月27日和2015年4月19日为起始时间的不同台站垂向原始波形记录，这两个波形事件在未滤波前初始P波震相较为清晰但后续S波震相肉眼较难识别，因此均未被福建台网拾取为地震事件，但利用PALM方法可识别为模板事件，并参与后续匹配滤波算法，可获得更多地震事件。

图  2  起始时间2015-01-27 T03:26:56 （a）和2015-04-19 T04:02:47 （b）下的不同台站垂直分量原始波形记录

Figure  2.  Vertical component original waveforms in different stations started in 2015-01-27 T03:26:55 （a） and 2015-04-19 T04:02:47 （b）

下载: 全尺寸图片 幻灯片

1.2   数据

本文收集了福建台网88个台站2015年1—4月和7—12月共10个月的连续波形数据。台站分布如图1所示，中部和南部台站密度相对较大，北部台站密度较为稀疏。将收集到的连续波形数据进行预处理，并在2—40 Hz范围内进行带通滤波，得到可进行PAL计算的波形数据。由于本试验使用的GPU具有显存为8 G的M4000显卡，在MSMS自动检测过程中，需要将福建地震台网一整天所有台站的连续波形读入内存中，因此，需要将88个台站的连续波形数据降采至每秒25个数据才能读入显卡内存中，以便快速进行GPU运算。另外，为了减少运算时间和降低噪声、脉冲标定等局部影响以确保地震事件识别的准确性，设定至少有6个台站同时接收到信号才被当作一个地震事件。PALM方法可高效地从连续波形中自动检测出地震事件，配合定位程序HypoInverse和HypoDD获得高精度震源位置，从而给出构造带的地震活动特性。

2.   计算结果

在进行STA/LTA特征触发式震相检测和寻峰式震相拾取算法时，将信噪比阈值设置为12，将长短窗平均能量比阈值设置为8，S波窗长设置为20 s，以避免福建地震台网外部及台湾地区的地震被计算。由于噪声、脉冲等非天然事件的随机性，当非天然地震事件被当作模板进行后续的MSMS自动检测并与原始连续波形进行互相关时，互相关系数会远小于设置的互相关阈值0.4，不会检测出新的地震事件，即被淘汰。为了使模板事件波形质量更高，更好地服务于后续模板匹配滤波检测，我们规定保留由最少6个台站参与计算的模板事件进行后续计算，在进行PAL计算后进行目录整理，最终获得地震事件模板919个。图3所示蓝色波形为PAL方法检测到的3个事件垂直分量波形示例，即3个模板事件，图3a蓝色波形为PAL方法检测出的2015年1月27日03:26:59发生在福建南平市附近的M_L1.5地震模板，为图2a原始波形进行预处理、滤波、降采后的波形，可从该图中不同台站垂向分量上清晰地识别出P波及S波震相，该模板事件可参与后续MSMS检测；图3b蓝色波形为PAL方法检测出的2015年4月19日04:02:49发生在福建漳州的M_L1.9地震模板，为图2b原始波形进行预处理、滤波、降采后的波形，该地震可看到清晰的P波和S波震相，该模板事件可参与后续MSMS检测；图3c为2015年12月6日04:43:23台湾花莲海域M_L4.0地震，从图中能看到清晰的初始P波震相和被误识别的S波震相，由于STA/LTA在S震相平均能量比较低，故不参与后续MSMS检测。

图  3  PAL模板事件垂向分量（蓝色）和检测到的地震事件的垂向分量（红色）

（a） 模板事件为2015年1月27日M_L1.5地震，检测到的事件为2015年1月12日M_L1.4地震； （b） 模板事件为2015年4月3日M_L1.9地震，检测到的事件为2015年4月3日M_L1.4地震； （c） 模板事件为2015年12月6日误识别事件，S震相识别错误

Figure  3.  Examples for PAL template events in vertical component （blue lines） and detected events in vertical component （red lines）

（a） The template event is on 27 January 2015 M_L1.5 and detected event is on 12 January 2015 M_L1.4；（b） The template event is on 3 April 2015 M_L1.9 and the detected event is on 3 April 2015 M_L1.4；（c） Error template event on 6 December 2015 with inaccurate S phase

下载: 全尺寸图片 幻灯片

利用上述模板事件对福建台网10个月连续波形进行MSMS检测，为了降低运算量以及噪音对地震事件识别的影响，我们同样规定必须至少有6个台站同时检测到地震事件才参与互相关计算。这个规定虽然减少了脉冲和噪音对地震识别的干扰，但同时也降低了对小震级事件的识别能力，使得最小完备震级M_c升高。我们设置互相关阈值为0.4，即模板事件和连续波形进行MSMS检测时，只有互相关值超过0.4才认为检测出新的地震事件。图3a和3b中红色波形为利用模板检测到的新地震事件，分别为2015年1-1月12日03:41:39 M_L1.4地震，震中为（26.20°N，119.04°E）和2015年4月3日03:48:16 M_L1.4地震，震中为（24.27°N，118.07°E），与相应模板事件的相关系数分别为0.6和0.5，而图3c未检测到新的地震事件。模板匹配自动检测算法共获得7万1 678个事件，但由于同一个地震事件会被不同的模板检测到而被重复计数，需要对重复地震进行筛选剔除。在进行HypoDD精定位的过程中，选取一维福建2015模型（图4d）作为一维速度模型，地震与台站的最大距离设置为200 km，采用共轭梯度法进行计算，共进行3组迭代，每组迭代8次，最终定位残差为63 ms，说明本次计算选取的参数较为合理，与此同时将发震时刻在4 s以内的地震事件认为是同一地震事件，以剔除重复地震事件。HypoDD精定位完成后，地震条目缩减到2 243条，地震分布如图4所示，主要集中分布在福建省中部和南部，尤其是沿长乐—诏安断裂带和永安—晋江断裂带分布明显。

图  4  基于PALM地震目录的地震事件分布与一维福建2015速度模型

（a） 地震事件水平分布图；（b） 震源深度分布沿纬度的投影；（c） 震源深度分布沿经度的投影；（d） 一维福建2015模型

Figure  4.  Events distribution of PALM catalog and 1-D Fujian 2015 P wave velocity model

（a） Horizontal distribution of the events from PALM catalog；（b） Vertical distribution of the events along latitude；（c） Vertical distribution of the events along longitude；（d） 1-D Fujian 2015 P wave velocity model

下载: 全尺寸图片 幻灯片

3.   结果分析

3.1   PALM目录分析

由于计算过程中，S波窗长选择、互相关阈值选择及至少6个台站记录到地震事件等限制因素，PALM目录对台网外部较远的海域即台湾海峡区域的滨海断裂带和广东南澳附近海域（23°N—23.5°N，117°E—117.5°E）的地震识别较少，不能如福建编目目录一样识别到广东南澳海域和台湾地区等地的地震。在政和—海丰断裂带与永安—晋江断裂带的交界处（25.5°N，117.7°E）附近的地震事件，PALM方法也无法识别。另外上杭—云霄断裂带东南段及福建广东交界位置识别出的地震较少。PALM方法在南平市和三明市附近（26.5°N±0.5°，117.5°E±0.5°）识别到福建地震台网目录未识别到的地震235次，原始连续波形如图2a所示，模板波形如图3a所示；在长乐—诏安断裂带上靠近漳州市—厦门市—泉州市—莆田市位置识别出了更多地震事件，尤其是漳州—厦门附近断裂带上识别到421次，原始波形如图2b所示，模板波形如图3b所示。仙游序列[（25.6°N，118.7°E）附近]在2015年处于平静期，但仍有小震级地震发生，PALM目录检测到了仙游序列大部分M_L1.0以上地震。在福州盆地及周边区域[（26.0°N，119.1°E）附近]，PALM目录识别出了更多地震，可为今后研究福州盆地构造提供依据。此外，PALM目录还识别出了福建台网网缘的浙江文成12次地震。

图4除了展示PALM目录水平分布外，还展示了PALM地震目录和台网编目目录中的地震的震源深度分别在不同经度和纬度上的投影（图4b和4c）。福建台网的地震定位主要利用人机交互系统中嵌入的单纯型和hypoSAT方法，这两种方法对震源深度的限制都依赖于含有深度信息的震相，如Pn震相（首波）和PmP震相（莫霍面反射波），当台站没有接收到这类震相时对震源深度的约束较差。福建地震台网目录的震源深度主要集中在5—25 km，小于5 km的震源较少，且深度只保留整数位。主要集中在广东南澳海域的地震在沿经度和纬度投影的震源分布图中显示为有一定间隔（间隔为1 km）的短横线。PALM目录利用HypoDD方法进行地震精定位后，可获得精度更高的震源位置，尤其对震源深度有更好的约束。PALM目录震源深度主要集中在0—20 km，且大部分地震深度小于10 km，表明地震主要发生在上地壳。

在2 243条PALM地震目录中，M_L1.0以上地震为1 991条，M_L1.1地震数目最多，为346条，整体呈现高斯分布特征，主要集中在M_L1.0—1.5 （图5a）。PALM目录中的震级由各个台站S波振幅的中位数获得，由于滤波的原因，其整体震级较台网编目地震震级偏小，存在系统偏差。图5同时给出了PALM目录和福建台网编目目录不同震级的数目及累计地震数目的统计和对比。由于PALM目录的地震事件需至少6个台站记录，而福建编目目录只需要不少于4个台站记录，因此利用zmap软件（Wiemer，2001）中的最大曲率法（Wiemer，Wyss，2000）计算福建编目目录最小完备震级M_c为M_L0.8，而PALM目录最小完备震级M_c仅为M_L1.1。福建编目有更丰富的M_L3.0以上目录，累计震级-频度曲线斜率更为平缓，而PALM目录对M_L1.0—1.5目录进行了大量补充，而小于M_L1.0小震和大于M_L2.0地震数目较少，使得累计震级-频度值增长较快，其斜率较大，平均b值增大至1.5。b值平均值较大说明PALM目录不完整，因此结合福建编目目录中小于M_L1.0和大于M_L2.0的地震研究福建地区地震活动特性更可靠。

图  5  PALM地震目录（a）和福建地震台网编目目录（b）震级-频度分布及最小完备震级M_c对比

Figure  5.  Frequency-magnitude distribution of PALM catalog （a） and Fujian Earthquake Agency catalog （b） and their complete magnitudes M_c

下载: 全尺寸图片 幻灯片

3.2   福建地区断裂带地震活动性

长乐—诏安断裂带是福建地区最活跃的断裂带，呈东北走向，内部发育一系列次级断裂带，以左旋压性活动为主，走滑速率为（0.7±1.3） mm/a，福建地块呈现出顺时针运动状态（占惠等，2015）。在PALM自动检测获得的地震目录中，有755次地震发生在长乐—诏安断裂带上，占地震总数的1/3左右；而福建台网编目目录中长乐—诏安断裂带上的地震仅59次。地震分布以永安—晋江断裂带为界，南西段地震条目571条，北东段地震条目184条，呈现南西段较强、北东段相对较弱的地震活动特性（图4），与福建沿海地壳水平形变南强北弱、东强西弱（陈光等，2010）的特征相一致。为了研究震源深度分布特征，我们沿着断裂走向提取了一条震源深度分布剖面AB （图6a），并将剖面附近20 km范围内的地震进行等深度投影。从图6b可以看出，地震震源深度主要集中在0—10 km范围内，从A点到B点震源深度呈逐渐变浅的趋势，与福建沿海莫霍面深度由西南向东北方向变浅的特征（金震，2018）一致。

图  6  剖面AB和CD位置（a）及沿两剖面的震源深度分布图（b，c）

Figure  6.  Location of profile AB and profile CD （a） and events distribution with depth along the sections AB （b） and CD （c）

下载: 全尺寸图片 幻灯片

永安—晋江断裂带的走向与福建地区主应力方向一致，整体构造活动较长乐—诏安断裂带弱，但历史上依然发生过中强地震（国家地震局震害防御司，1995；中国地震局震害防御司，1999）。永安—晋江断裂带的左旋走滑速率和倾滑速率分别为（0.6±1.8） mm/a和（−0.4±0.5） mm/a （占惠等，2015），其地震动力作用与菲律宾海板块对台湾岛弧作用密切相关，动力作用波前峰大致位于泉州湾，断裂带东南段以拉张走滑运动为主（林松建，黄宗林，2013）。根据PALM目录，其整体发生地震较长乐—诏安断裂带少，主要分布在断裂带两端，地震数目为180条。图6a中剖面CD为沿着该断裂带从北西至南东提取的地震分布剖面，同样将剖面附近20 km范围内的地震等深度投影至剖面上（图6c）。从图6c中可以看出，剖面CD上的地震分布以政和—海丰断裂带为界，可分为两段：北西段地震分布更为集中，深度主要集中在10 km范围内；南东段地震分布较为稀疏，深度主要在20 km范围内，中部存在地震平静，与其动力环境具有一致性。

4.   讨论与结论

采用PALM自动检测方法对福建台网2015年10个月的连续波形数据进行自动识别与重新定位，共获得2 243条地震目录。震级主要集中在M_L1.0—1.5之间，震源深度主要集中在0—20 km范围内，且大部分地震发生在上地壳内，深度小于10 km。PALM目录对福建台网M_L1.0—1.5地震事件进行了极大的扩充，在台网的西北边缘南平、三明一带识别出台网目录丢失的235个地震事件，在福建沿海长乐—诏安断裂带多识别出696个地震事件，为地震活动性分析提供了较福建台网更为完备、定位精度更高的地震目录。但信噪比阈值和至少6个台站接收波形等参数设置仍在一定程度上限制了震级下限及网外海域地震的检测能力。

重新识别与精定位后的地震活动沿长乐—诏安断裂带与永安—晋江断裂带分布集中，且表现出明显的分段活动性：沿长乐—诏安断裂带震源深度分布由南向北逐渐变浅，与莫霍面深度的变化相一致；在永安—晋江断裂带北西段地震活动更为集中。福建地区地震活动性与水平形变南强北弱、东强西弱的特点相一致。PALM方法可以在福建地区取得较为可靠的地震目录，为后续依靠地震活动性展开的研究提供了重要依据。

审稿专家为本文提出了宝贵的修改意见和建议，福建省地震局提供了2015年台网连续波形数据，文中图件使用GMT （Wessel，Smith，1998）制作生成，作者在此一并表示感谢。