GPU 并行加速下的弱模板匹配及其在川滇地区的应用

2021-07-14蒋一然宁杰远李春来

华北地震科学 2021年2期

蒋一然，宁杰远,2*，李春来

（1.北京大学地球与空间科学学院，北京 100871；2.河北红山地球物理国家野外科学观测研究站，河北邢台 054000；3.中国地震局地球物理研究所，北京 100081）

0 引言

模板匹配方法是检测微小地震、慢滑移等活动的重要方法[1-4]。对于位置接近、震源机制解相似的2 个地震，其到同一台站的路径相近，记录也相似。模板匹配方法就是使用其中一个地震的波形作为模板，与连续记录互相关，当互相关值较高时则认为可能发生了相似地震。单台的结果容易受随机噪声的干扰产生误识别，因而将多个台站的互相关结果叠加起来以减少误识别的比例。同时，由于噪声的影响，能量较小的震相与模板波形的互相关值虽然会略高于平均水平，但并不显著。借助多台叠加的方法，随机噪声引起的互相关值波动被平均，来自微小地震的高互相关值因为时间上的一致性被加强，使得小能量的微小地震信号变得显著。这也使得模板匹配方法可以检测到单台上能量与噪音水平相当，甚至略低于噪音水平的微小地震。

与模板位置不同的微小地震的识别具有重要的意义，不仅能够识别到更大范围内的地震，还能用于断层微结构研究。但是，当微小地震与模板位置有一定差距时，即使波形相似度较高，但由于波场到达台站的走时与模板地震不同，不同台站的互相关值峰值时间不一致，也不能叠加起来，因此无法被检测到。一种改进的思路[5]是在将模板地震的位置在原始位置附近进行格点式的挪动，根据与原始位置的位置差和参考的速度结构对时间移动的大小做微小调整，从而提高模板检测微震的空间范围。在这种方法中，格点移动的范围、移动的步长都将显著地影响计算量，阻碍高精度定位；参考模型选取得合适与否，也将对检测能力造成较大影响。另一种思路则是扩宽单台上高互相关值峰的宽度，使得与模板走时差在一定范围内的互相关峰值也能够叠加起来[1,6]。这种扩宽单台互相关峰值后叠加的方法被称为弱模板匹配方法（weak matched filter technique，WMFT）。弱模板匹配方法在不显著增加计算量的情况下提高了单个模板探测微小地震的半径，因而是一种较好平衡计算代价和检测能力的微小地震识别方法。

对于模板匹配而言，模板多可以提高检测微小地震的能力。对于缺乏人工目录和到时标准的流动台网，模板相对缺乏。而在有大量地震目录和波形积累的固定台网中，人工标注的部分也往往是信噪比较高的部分，对应的地震震级较高。一般来讲，震级差距较大的地震，震源时间函数差距也较大，台站的波形记录相似度也较低。震级较高的地震所能检测到的微小地震的震级也会偏高。近年来，大量以深度学习方法为基础的地震自动检测和震相拾取方法被提出[7-11]，大大减少了获取地震目录和震相走时的成本，并提高了检测结果的一致性。借助深度学习的震相自动拾取算法，可以获得震级下限更小的地震目录作为模板，从而探测震级更小的微小地震。

由于每个模板都需要和连续记录进行互相关，即使是已经大大减少运算代价的弱模板匹配，仍然需要进行大量的运算。互相关运算对于每段波形的扫描结果是独立的，可以很好地并行化。使用torch[12]软件包实现了弱模板匹配中互相关、归一化、拓展互相关峰值3 个主要计算部分的GPU 并行，极大地提高了运算速度，使得在相同算力的情况下，可以使用更多的模板，从而检测出更多的微小地震。

本文将并行优化后的弱模板匹配方法应用到川滇地区，以基于深度学习的APP++震相拾取方法[11]拾取到的地震[13]作为模板，扫描了2019 年6 月17 日四川长宁地震前后各1 个月共60 天的连续波形，检测到81 704 个微小地震，并筛选出7 618 个可靠的微小地震。在由双差地震层析成像[14-15]确定的模板地震位置的基础上，利用主事件法对模板探测到的微小地震位置进行定位，并分析了地震序列的时空分布特征。

1 方法和数据

1.1 GPU 并行下的弱模板匹配

包括弱模板匹配在内的模板匹配方法面临的主要困难之一是庞大的计算量，弱模板方法因计算量更大，问题尤为突出。提高模板匹配的计算效率，就可以在相同的时间内，使用更多的模板扫描更多的连续记录，从而检测更多的微小地震，从而得到有意义的研究结果。为了尽可能地提高模板匹配方法的运行效率，本文具体分析了弱模板匹配方法的计算过程和主要的运算代价部分，借助深度学习软件包torch 中的相关函数，实现了弱模板计算的GPU 并行化。

弱模板匹配方法的第一步是将模板地震的每一个震相波形与相应台站的连续数据互相关，得到互相关值的时间序列。具体地讲，对于一个模板，其第i个波形记录wi(t)与相应台站的联系记录ri(t)做互相关得到的互相关时间序列corri(t)如下：

式中：t表示时间；dt是时间采样间隔；Ti为用于互相关的模板波形的起始时刻；N为用到与互相关的采样点的个数；To为模板地震的发震时刻。这里使用了(Ti-To)对不同走时的震相的互相关结果进行了对齐，分母用于将互相关结果归一化到0～1 之间。

深度学习算法中，一般不关心卷积核中参数的时序，因而一般用互相关操作来进行等价卷积，以简化程序的实现。具体到torch 软件包中，用一维的卷积操作conv1d（实际为互相关）来实现在GPU中的并行计算。值得注意的是，这里的卷积操作并没有对每一段进行归一化。因此，作为归一化因子的分母需要单独计算。分母的计算由两部分构成，一部分是模板波形的平方和，另外一部分是每时刻用于互相关的连续波形的平方和。第一部分对于每个时刻是不变的，计算量较小；第二部分是随时间变化的，计算量与归一化因子的分子相当。在此将第二部分的计算转化为互相关计算：

式中：1(t)表示一个值恒为1 的时间序列，分母第二部分的计算相当于ri的平方与这个序列做互相关，也可以使用conv1d 函数予以实现。在北京大学实验室的工作站上，对使用CPU（AMD 2990wx）和GPU（1 080 Ti）计算互相关值的平均用时进行了统计（图1），结果显示使用GPU 的平均耗时约为CPU的1/40。需要特别指出的是，如果使用频率域互相关的方法，运算速度会更快。但是，由于连续记录时间长度长、波形能量变化大，采用频率域计算的方法，会面临大数吃小数的数值问题，对于互相关值的计算精度影响很大。在测试过程中，conv1d 函数在CPU 上进行时域互相关时仅调用了2～3 个核心，并行度不高。如果提高CPU 上的并行化程度，能够在一定程度上提高计算效率，但GPU 仍然有较大的效率优势。

图1 CPU-GPU 互相关效率测试

在得到归一化后的互相关值后，弱模板匹配方法要求需要对每一道互相关值的峰值进行拓宽，得到新的互相关值序列。具体地讲，根据直接互相关序列的均值m和标准差s，设定单台的互相关值阈值threshold为：

式中：mul为峰值需要大于平均值的最低标准差倍数。当互相关值大于阈值时，视为需要拓宽的峰值。这里倍数mul的选取，需要保证能够正确拓宽微小地震的高互相关值，又不过多地拓宽随机噪声引起的互相关高值。经过实验，一般取mul为4。

从实际操作上讲，这个拓宽的过程和深度学习中的最大值池化过程maxPooling 类似。最大值池化是选取时间点附近一定时间窗内的最大值代替原有值，即：

式中：twin表示选取的用于池化的时间窗大小，等于需要拓宽峰的时间长短。时窗长度变长可以增加模板的检测范围，但同时也会增加误识别比例，其选取需要平衡这二者。带阈值的峰的拓展中与最大值池化略有不同的地方在于，当最大值小于阈值时，就不使用最大值对原有的值进行替代。因此，最终得到的拓宽后的互相关值序列corriw(t)为：

这样，就可以借助torch 软件包中的maxPooling函数实现并行的峰的拓宽。结合上面的处理方法，就可以将弱模板主要计算部分全部通过GPU 并行加以实现。在实际扫描过程中，一般先把一天的连续记录加载到内存中，然后与多个模板进行互相关。考虑到数据在内存和显存之间交换会耗费相当的时间，预先将波形数据加载到显存中。一般来讲，因显存的大小有限，无法同时将一天的连续数据和全部的模板波形加载到显存中。由于每个台站的连续记录会和很多模板在这个台站上的波形进行互相关，而每一次互相关中连续数据的数据量大于模板波形的数据量，所以可以优先将连续记录加载到显存中。

之后，将不同模板波形的互相关结果进行平均，得到总的互相关叠加结果corrall(t)：

式中：I表示所有使用到的模板波形数量。根据corrall(t)的均值M和标准差S设定检测微小地震的阈值Threshold为：

式中：Mul是需要给定的微小地震互相关值高于平均值的标准差倍数，一般取Mul等于6，得到初步的微小地震目录后再进行分析和质量控制。

1.2 滤波频段的选取

在进行互相关之前需要对数据进行滤波，对微小地震波形主要能量所处频段外的其他信号进行压制。在模板匹配中，滤波频段的选择还应该考虑使得相似地震的互相关值更容易具有较高的互相关值。有两个完全相同的波场，都记为wi(t)，两者之间做互相关。在（1）式的定义下，互相关值在T0处取得最大值1。但是由于离散采样，实际记录中两个波场的采样点会有一定差异，设其为dt′∈，dt为数据的采样间隔。此时实际只能获得T0+dt′处的互相关值：

将wi(Ti+n×dt)和wi(Ti+n×dt+dt′)作傅立叶展开：

在dt′的偏移下，同一频率的信号会有一个相差，越高频的部分，相差越大，越容易造成互相关值的降低。通过压制高频部分的能量，可以在一定程度上减少采样点差异对互相关值的影响。

本文中使用数值试验的方式来测试在采用不同的滤波频段对互相关值的影响。先生成1 000 Hz采样的随机信号当作连续波场的近似，并将波场降采样到50 Hz，生成两段采样点相差0.01 s 的波形，并且进行0.5～20 Hz、0.5～12 Hz、0.5～8 Hz 三个不同频段的带通滤波。将在同一频段滤波但有采样点偏差的记录做互相关（图2）。结果可以看出，当波形包括较高频段的能量时，即使是相同连续波场下采样的波形互相关值也显著低于1。随着带通滤波最大频率的下降，互相关值逐渐上升。当取0.5～8 Hz 时，互相关值趋于稳定。因此我们选取0.5～8 Hz作为弱模板匹配的互相关频率。

图2 采样点有偏差时相同波场的互相关结果在不同带通滤波下的变化测试

1.3 微小地震的相对定位

在检测到微小地震后，在每个台站的直接互相关结果对应时间点附近选取互相关值较高的点作为微小地震在该台站上的震相到时，并以叠加的互相关值峰值位置的时刻作为对发震时刻的初步估计。之后采用主事件法，计算小地震相对于模板地震的位置。用x表示模板地震的位置，x′表示微小地震的位置，ti和τio分别表示模板地震和微小地震的第i个波形的走时。这里 τio是由微小地震初步估计的发震时刻计算的。当对微小地震的发震时刻调整 Δτ时，获得新发震时刻下的走时 τi：

考虑到2 个地震事件的位置十分接近，于是根据模板地震走时及其对震源位置的一阶导数来估计微地震的走时作为理论到时 τit：

式中：I表示使用到的震相波形的总数；wi是根据每个震相波形的互相关值设定的权重。在此设定为当互相关值大于0.4 时，权重wi为互相关值的平方，否则设定为千分之一。由于模板地震一般具有较高的能量，走时和绝对位置更容易确定。而微小地震能量较低，绝对位置确定较困难。这里依赖于模板地震，计算微小地震相对于模板地震的相对位置，则更为容易。

1.4 模板地震与连续数据

蒋一然等[11]基于U-net++[18]网络改进了APP 方法，提出了基于台阵策略APP++震相自动拾取算法，并使用该方法扫描了川滇地区的固定台网（116 个）和小江断裂带附近的高密度流动台阵（51 个）2014—2019 年间6 年的连续波形数据（台站分布如图3），共检测到73 291 个地震，537 554 个P 波记录，471 459 个S 波记录。在此基础上，蒋一然等[13]挑选了其中的21 160 个地震，利用双差地震层析成像方法，获得了该地区浅部的三维速度结构和17 240 个地震的精确空间位置。本文使用这17 240个地震作为研究模板地震，并依据其精定位结果对微小地震进行相对定位。需要特别指出的是，蒋一然等[13]在进行双差地震层析成像的时候，使用了由互相关提取的走时差作为约束，而本文中对微小地震定位也是使用了互相关提取的走时差，两者具有一致性，因而不同模板扫描得到的微小地震的位置间也是可比较的。

图3 APP++方法同时间段内识别到的地震

下面选取与模板地震同一台阵的2019 年6 月17 日长宁地震前后一个月共60 天的连续波形，使用GPU 并行下的弱模板匹配技术扫描其中的微小地震。

2 弱模板匹配结果

2.1 初步结果和质量控制

在50 Hz 采样率下，设定带通滤波频带为0.5～8 Hz，震相波形截取为震相到时前3 s 和后4 s；模板地震至少含有6 个震相波形，如果超过40 个震相波形，则取相应台站震中距最小的40 个震相波形；对于单台互相关结果，当某些时刻的互相关值大于整体平均值4 倍标准差时，将其后面0.45 的互相关值全部置为该互相关值；取叠加结果中高于平均值6 倍标准差的峰值为检测到的微小地震。一共检测到81 704 个地震，图4 展示的是其中的2 个例子。

图4 两个弱模板匹配图

进一步统计了相同时间段内，在该区域由弱模板匹配和APP++方法拾取到的地震震级频次分布（图5a）。可以看到，APP++的结果在完备震级之上基本满足G-R关系[19]，而弱模板匹配的结果则在震级小于2 时，数量明显增多。经分析，这是由于设计的阈值过低，导致一部分噪声被误识别为地震，从而使得地震目录的G-R关系被破坏。通过统计不同Mul取值下，弱模板匹配检测到地震数目的变化（图5b），发现在8～12 之间，地震数目的变化和Mul取值变化间满足很好的指数关系。而当Mul小于8 时，误识别比例增加，地震数目变化不再符合之前的指数关系。故将Mul设定为8，重新统计了地震震级频次分布（图5a），发现提高阈值后的震级分布基本满足G-R关系。和APP++方法比较，质量控制后的弱模板匹配方法的地震检测能力比APP++方法提高约0.5 个震级，数目约为APP++方法的3 倍（APP++拾取到2 779 个地震分布如图3，弱模板拾取到7 618 个地震分布如图6）。

图5 地震数目统计图

图6 筛选后弱模板匹配检测到的地震分布图

2.2 地震序列

图6 中，地震的震中位置与当地的断层位置有很好的一致性。图中剖面水平线汇聚的位置为长宁地震主震及其余震发生的区域。选择红色矩形框内的地震，统计其发震时刻和震级分布（图7）。图7a～7b 中分别展示了APP++方法拾取到的地震序列和弱模板匹配方法拾取到的地震序列。APP++地震序列中，可以看到有4 个左右时间上的丛聚序列（近垂直的细条）；而在弱模板匹配目录中，借助更多小震发现较多时间上的丛聚序列。另外，由于弱模板匹配识别到的地震事件较多，对于地震活动性的分辨力也会更高。由弱模板匹配方法的目录中可以发现，在长宁地震发生前有一个较为明显的平静期，这可能与地震临震状态的动力学过程有关。借助弱模板匹配方法识别到的地震目录可以更好地研究地震活动性及其相应的构造过程。

图7 2019 年6 月17 日长宁地震前后地震震级随时间的变化

2.3 地震的空间位置分布

以图6 中的4 条线段为剖面，展示地震在垂直剖面上的位置（图8）。A—A′剖面120～200 km 之间和C—C′剖面170～250 km 之间，地震丛集大致为南北两个中心，主要集中在20 km 以内，更深部的地震没有体现出较为明显的取向；B—B′的120～200 km和D—D′剖面200～300 km，两个丛集的垂直投影相互重叠不易分辨，深部的地震体现出一定的线性特征，可能与该位置的断层结构有关；D—D′剖面150 km处出现的地震丛集体出现明显的取向，反映地震在相应断层发生的位置分布。获得的微小地震空间位置分布可以进一步用于对断层几何形态、地震活动性及地震活动性变化的研究。