微小地震虽然破坏性不大，却是研究大地震和整个地震带应力场的基础^[1]，其特点是频率高、信噪比低、数量多，若能实现微地震事件的自动识别和初至震相的自动拾取，将对海量微地震数据的自动处理有重要意义^[2]。Stevenson^[3]首次提出应用长短时窗平均能量比值法(STA/LTA)识取地震初至震相，其后发展出各种基于该方法的时间域、频率域、偏振分析、高阶统计量、模式识别等综合分析震相拾取方法^[4]；Gelchinsky等^[5]提出用相关法自动检测走时曲线初至到时；Schaff等^[6]提出基于波形互相关的重复地震检测方法，并在此基础上发展了一种模板匹配技术，用于检测波形相似的地震^[7]。由于模板匹配技术能够有效抑制背景噪声、较大地震面波和震荡等造成的低频信号干扰^[8]，在地震自动识别领域被广泛应用^[9-11]。与STA/LTA方法相比，波形互相关法在波形相似度较高的情况下能够检测出更多的地震^[10]，但当信噪比不高、地震记录信号与噪声频带重叠时，其检测效果并不理想。

小波相似性分析是研究时间序列多尺度相似性的一种新的时频分析方法。该方法基于多尺度分析框架并考虑局部相位，可较好地再现2个序列时间尺度上的相似性，同时也保持了局部的细节特征，进而揭示2个序列间多尺度的相互作用关系^[12]。Johnson等^[13]首次应用小波相似性算法分析地球物理钻孔数据，并对不同情况下数据的相似性进行讨论。本文提出一种基于小波互相关技术的微地震波形自动检测方法，利用复Morlet小波在特定尺度上对地震波形进行变换，分析其相位角相关性，设置合理的相似度阈值，自动检测微地震。

1 基本原理

由于Morlet小波采用时域窗面积最小的高斯窗函数，时域、频域局部化性能都很好，对称性也较好，故本文选用该小波作为母小波进行地震信号检测。地震信号f(x)的连续小波变换CWT^[14]为：

$ \mathrm{CWT}(u, s)=\int_{-\infty}^{\infty} f(t) \frac{1}{|s|^{0.5}} \psi^{*}\left(\frac{t-u}{s}\right) \mathrm{d} t $

(1)

式中，s为小波分析尺度，u为位移，ψ为小波母函数，*表示共轭复数。用复Morlet小波作为小波母函数，定义如下^[15]：

$ \psi(x)=\frac{1}{\pi f_{b}} \mathrm{e}^{2 \pi \mathrm{i} f_{c} x} \mathrm{e}^{-x^{2} / f_{b}} $

(2)

式中，f_b为复Morlet小波的带宽，f_c为复Morlet小波的中心频率。

对于地震数据X(t)，由于同一地震在同一台站内应具有相同的频率特征，故可以利用连续小波在某些特定尺度上的相位角相似性来检测获得一段时间内的所有地震事件。

定义2个信号的连续交叉小波变换^[16]为：

$ \mathrm{CWT}_{1, 2}=\mathrm{CWT}_{1} \times \mathrm{CWT}_{2}^{*} $

(3)

式中，CWT_{1, 2}的幅值A为|CWT_{1, 2}|，相位角θ为$\arctan \frac{I\left(\mathrm{CW} \mathrm{T}_{1, 2}\right)}{R\left(\mathrm{CWT}_{1, 2}\right)}$，－π≤θ≤π。用θ的余弦值表示2组地震数据的相似度，即

$ S=\cos ^{n} \theta $

(4)

式中，n为大于0的奇数。当S=1时表示数据正相关，S＝0时表示数据不相关，S=－1时表示数据负相关。

当信号的背景噪声很大时，相似性的变化对噪声较为敏感，所以考虑信号的幅值信息，2组数据的相似度又可表示为：

$ S=A \cos ^{n} \theta $

(5)

2 模拟信号的小波相关性分析

假设有2个时间序列x₁、x₂分别为：

$ x_{1}=A \cos \omega_{1} t+A R \times \operatorname{rand} n(1, n) $

$ x_{2}=B \cos \omega_{2} t+B R \times \operatorname{rand} n(1, n) $

对x₁和x₂进行小波相关性分析，波形及分析结果见图 1。图 1(a)为x₁数据，是一组10 Hz正弦信号掺杂50 Hz的正弦信号；图 1(b)为x₂数据，是一组50 Hz的正弦信号；图 1(c)、1(d)分别为x₁、x₂数据的复Morlet连续小波变换结果；图 1(e)为2组数据进行Morlet小波相位角相关性分析的结果。从图中可以看出，在2组信号频率相同的地方表现出明显的相关性。

由于小波分析的多尺度特性，小波相关分析具备明显的尺度特征，即可以在不同尺度下进行相关性比较。这给相关性分析带来了很好的频率指向性，从而能够进行更精细的数据相关性分析，因此可以使用统计学方法进行2组波形相似度的定量计算。

假定x₁为基准波形，试图从序列x₂中寻找与x₁最相似的数据段时，可以通过计算x₁在不同尺度下的小波系数的能量权重，将该权重引入相似度中。取数据序列中所有数据的平均值，可以标定2组数据序列的定量相似性，表示为：

$ S_{\mathrm{avg}}=\frac{\sum\limits_{n=1}^{N} K(n) \times S(n)}{(N \times M)} $

(6)

式中，K(n)为x₁在各尺度下的小波系数的能量权重，S(n)为x₁与x₂进行相关性分析的结果，N为x₁与x₂进行相关性分析的每个尺度下的系数数量，M为x₁与x₂进行相关性分析的尺度。

在x₂逐点取得与x₁同长度数据，使用式(1)~(6)计算信号x₂中任意位置信号片段与信号x₁的相关性，得到相关性系数集合，在集合中寻找满足离群条件的数据，则可以找到信号x₂中与模板信号x₁相关的信号片段。这个过程可以用绝对离差中位数(median absolute deviation，MAD)计算。

对所有数据进行MAD计算，并约定阈值Thr，则可以基本确定序列x₂中与x₁较为相似的数据段位置，即

$ \mathrm{MAD}=\operatorname{median}\left(\left|X_{i}-\operatorname{median}(X)\right|\right) $

(7)

$ \mathrm{Thr}>\kappa \mathrm{MAD} $

(8)

继续上述实验，在序列x₂中寻找与x₁相似性的定量分析，结果见图 2。图 2(a)、2(b)分别为用式(4)、(5)对x₁和x₂进行小波相位角相关性定量分析的结果。从图 2可看出，该方法能直观地确定相关数据的具体位置，从而自动检测强地震后的余震或震群。

3 实际地震的小波相关性分析

2018-05-01 11:14新疆昌吉发生M_L4.8地震(43.61°N，86.62°E)，其后发生一系列余震。主震与余震通常有共同的发震构造和互相关联的震源机制，波形也具有很好的相似性^[17]。为验证本文方法的可靠性，选取离震中最近的STZ、LHG和HTB台记录的2018-05-01 16:00~16:30三分量波形为待检测数据，以信噪比高、记录清晰的地震波形为模板地震。首先用2~8 Hz四阶巴特沃斯(Butterworth)带通滤波器对地震波形进行去均值及归一化处理，然后对模板地震三分量数据进行复Morlet连续小波计算，计算其主要能量集中的尺度范围，并对这个尺度范围内的小波数据进行相位角计算，得出相位角序列。根据此结果截取待检测数据中相似的地震数据，对每一次截取的数据进行复Morlet小波变换，在能量集中的尺度范围内用式(4)或式(5)、式(6)和式(7)进行小波相位角相关性计算。计算MAD时采用三分量数据小波相位角相关系数的平均值，阈值选取根据实际地震不同而有所不同，本文选取40倍MAD作为判别地震的阈值Thr。地震检测结果见表 1，30 min内共检测出11个地震，与中国地震台网中心地震编目结果基本一致，无漏检、错检。比较P波到时，STZ台相差约1 s，LHG台相差约为0.05 s，HTB台相差约为0.26 s。图 3为2018-05-01 16:20:50~16:26:10的地震检测情况，图 3(a)、3(c)及3(e)为在连续波形上检测到的地震P波初至的位置，图 3(b)、3(d)及3(f)为检测波形示意图。STZ台和LHG台记录的波形信噪比高，用式(4)进行计算；HTB台记录波形背景噪声大，用式(5)进行计算，所有结果均能清晰地检测到地震。由此可见，将该方法用于地震序列中微地震检测具有一定的可行性。

4 结语

本文利用复Morlet小波相位角互相关变换分析数据的相关性，通过统计学方法对数据相似度进行定量计算，获得原始数据中与模板数据较为相似的数据段初始位置，进而将该方法用于地震序列微地震检测。对2018-05-01新疆昌吉M_L4.8地震后30 min的波形进行检测测试，共发现11个微地震事件，并给出其检测台站的P波到时结果。与中国地震台网中心地震编目结果进行对比发现，无漏检、错检，P波到时结果相差基本在1 s之内，认为将该方法用于地震序列微地震检测具有一定的可行性。

本文所用方法的一个显著特点是使用三分量地震波形进行小波变换，利用小波相位角互相关检测地震序列，同时标定检测出的地震序列P波到时。与前人用互相关检测地震的方法^{[8-9, 11]}相比，该方法应用小波变换将时间序列分解到时间频率内，从而得出时间系列显著的相位变化特征。由于小波变换在时频两域都有表征局部特性的能力，对地震噪声有很强的压制作用，即使地震记录信号与噪声频带重叠，也能通过变换尺度突显检测信号的瞬态或奇异点。

本文所用方法选取的模板地震为信噪比大于4、震级M_L≥1.5、至少有3个以上的台站清晰记录，若震级太小，则能记录到该地震的台站不多，不具有通用性，如第1个模板地震。作为方法检验仅选择了30 min的数据，用到3个模板地震，对于地震序列较多、待检测波形时间更长的情况，采用该原则选择的模板地震会有很多，若进行实际操作，工作量将非常庞大。在以后的工作中，笔者会选择几个小时的数据进行试算，若有N个地震都被检测到，则确定该地震为模板地震，N的确定要保证不漏检、错检和多检地震。不可否认，在检测过程中，若存在与任何模板地震波形均不互相关的地震，该地震将难以被检测到。

目前，该方法仅用于检测相似微地震的P波初至，对于形成完整的地震序列目录还存在一定的局限性，今后会尝试用该方法检测微地震的S波初至。

致谢: 感谢新疆维吾尔自治区地震局阿克苏中心地震台龚固斌对相关公式给予的指导。