由于微地震事件信号通常能量较弱, 加之传播过程中能量被地层吸收衰减, 造成微地震事件淹没在噪声中, 对后续的微地震事件精确定位产生了不利影响, 因此提高微地震资料的信噪比是微地震数据处理和解释的首要任务^[1-2]。随着定位精度和震源机制全矩张量反演需求的增加, 对去噪技术的要求也逐步提高^[3]。传统压制随机噪声的方法可分为空间域及变换域两类, 前者主要包括均值滤波、中值滤波及各向异性扩散滤波等, 后者主要包括傅里叶变换域滤波方法及基于小波变换、曲波变换等的阈值去噪方法等^[4-5]。地面微地震资料具有强噪声、弱有效信号的特点, 常规的去噪方法往往难以获得较好的去噪效果^[6]。针对经验模态分解(EMD)方法的不稳定问题和模态混叠现象^[7-9], 总体经验模态分解(EEMD)利用高斯白噪声频谱均匀分布的统计特性, 在原始信号中加入不同的白噪声, 使得信号在不同尺度上具有连续性, 但是该方法计算效率不高^[10-11]; 完备总体经验模态分解方法(CEEMD)通过加入正、负成对的辅助噪声, 能够有效消除重构信号中的残余辅助噪声, 而且计算效率较EEMD有所提高^[12], 但是信号重构精度会有所降低; 而TORRES等^[13]提出的改进的完备总体经验模态分解方法(ICEEMD)能够精确重构原始信号, 有效减少虚假模态和模态中的噪声, 同时也能降低计算成本。

本文研究了基于ICEEMD的有效信号提取方法, 将非平稳信号分解为一系列相对平稳的固有模态函数, 定义了一种自适应间隔阈值, 用来去除固有模态中的噪声成分, 通过重构数据得到去噪后的结果, 同时利用分解后的信号与Hilbert变换进行初至检测。最后用数值模拟数据及实际微地震资料验证了方法的有效性。

ICEEMD方法将待处理微地震信号看作初始数据, 在分解的每一阶段添加一个特定白噪声, 计算一个固有模态函数(IMF)分量, 然后得到一个残差, 进而得到每个IMF分量, ICEEMD能自适应地将一个复杂信号分解为一系列IMF分量, 且该分量满足从高频到低频的顺序分布。ICEEMD克服了EMD原有的模态混叠等缺点, 保持了EMD的完备性, 能够进行原始信号的精确重构, 且有更好的收敛性^[14-15]。

ICEEMD具体步骤如下。

1) 通过EMD实现加入不同噪声后信号的分解计算, 获得第一个IMF分量, $I_{\mathrm{MF}_{1}}=\frac{1}{m} \sum\limits_{i=1}^{m} E_{1}\left(x^{i}\right)$。

2) 计算一级残差r₁=x－I_MF1。

3) 然后通过EMD实现r₁+ε₁E₁(ωⁱ), 得到第二个IMF分量$I_{\mathrm{MF}_{2}}=\frac{1}{m} \sum\limits_{i=1}^{m} E_{1}\left[r_{1}+\varepsilon_{1} E_{1}\left(\omega^{i}\right)\right]$, 获得二级残差(r₂), r₂= $r_{1}-I_{\mathrm{MF}_{2}}=r_{1}-\frac{1}{m} \sum\limits_{i=1}^{m} E_{1}$[r₁+ε₁E₁(ωⁱ)]。

4) 对于k=3, 4, …, K, 计算第k个IMF分量$I_{\mathrm{MF}_{k}}=\frac{1}{m} \sum\limits_{i=1}^{m} E_{1}\left[r_{k-1}+\varepsilon_{k-1} E_{k-1}\left(\omega^{i}\right)\right]$。

5) 计算k级残差r_k=r_k－1－I_MFk。

6) 重复步骤4)和步骤5)直至残差不能被分解。

上述各式中, x表示输入的原始信号; r_k表示k级残差; 算子E_j(·)表示对给定信号通过EMD求得第j个模态分量; wⁱ为单位方差均值为0的高斯白噪声, i=1, 2, …, m, m表示噪声组总数; xⁱ=x+wⁱ为加入不同噪声后的信号; ε_k表示噪声标准差, 在每一个分解阶段表示一个常数, 允许在每个阶段选择信噪比。

图 1分别给出了原始信号时域波形(全文需要分解的原始信号不做特殊说明均为各图中的c0)及经过EMD、EEMD、ICEEMD分解后的IMF分量。通过对比可以看出, EMD在分解信号时出现模态混叠, 不能将各个模态分量有效分解出来, 而EEMD与ICEEMD分解效果相近, 均可以将不同模态分量有效地分解出来。

图 1(Fig. 1)

图 1 原始信号时域波形及分解后的IMF分量 a EMD分解; b EEMD分解; c ICEEMD分解

对主要包含高斯噪声的IMF分量可以通过高阶统计量(HOS)和Kurtosis准则进行检测和去除, 留下混合噪声和信号的分量。包含N个数据的第i个IMF的Kurtosis峰度(K_i)可通过下式计算^[16]:

式中:σ_IMFi和μ_IMFi分别是I_MFi的标准差和平均值, i=1, 2, …, K。区分高斯分布和非高斯分布的HOS准则为:

式中:α是置信水平, RAVIER等^[17]给出了优化值α=90%。K_i若满足(2)式即为高斯分布, 需要在信号中滤除。该过程类似自动带通滤波, 去除地震信号频带相对较低和较高的频率成分, 通过此步骤可以去除高能量相干噪声。

阈值去噪方法的关键是阈值和阈值函数的选取。传统的阈值函数主要有硬阈值函数和软阈值函数, 若选取的阈值过大, 虽然去噪彻底, 但同时造成了有效信息的丢失; 若选取的阈值过小, 则导致去噪不彻底^[18-19]。本文采用了自适应间隔阈值:

式中:h(z_i)是输入信号相邻零点之间的采样间隔; h(r_i)是这个间隔中的局部极值; i=1, 2, …, I; 0 < γ < λ, 0≤α≤1, γ是截断值, 截断值以下数值设为0;α决定了阈值函数的形状。该阈值函数类似于硬阈值, 但门槛值是平稳过渡的, 因此它可以视为硬阈值和软阈值的组合。如果该高斯噪声的方差为σ², 则对于J个采样点的优化阈值为^[20]:

图 2a为含噪声合成信号时域波形及分解后的IMF分量; 图 2b为去噪后信号时域波形及分解后的IMF分量; 图 3a为去噪前的时频分布; 图 3b为去噪后的时频分布。合成记录中含有大量噪声, 信噪比很低, 对比去噪前后波形可看出, 噪声被基本去除。从分解后的模态分量可以看出, 噪声分量得到去除; 从时频分析结果也可以看出合成记录中的噪声得到有效去除, 信噪比大大提高。

图 2(Fig. 2)

图 2 含噪声合成信号时域波形及ICEEMD分解后的IMF分量(a)及去噪后信号时域波形及ICEEMD分解后的IMF分量(b)

图 3(Fig. 3)

图 3 含噪声合成信号去噪效果分析 a去噪前时频分布; b去噪后时频分布

图 4a为含噪声微地震信号时域波形及分解后的IMF分量; 图 4b为去噪后信号时域波形及分解后的IMF分量; 图 5a为去噪前的时频分布; 图 5b为去噪后的时频分布。所选微地震信号信噪比较低, 无法有效识别初至。对比去噪前后波形可见, 噪声被基本去除, 在分解后的各个分量中, 噪声分别得到有效去除; 从时频分析结果也可以看出去噪效果良好, 噪声基本得到去除, 信噪比得到极大提高。

图 4(Fig. 4)

图 4 含噪声微地震信号时域波形及分解后的IMF分量(a)及去噪后时域波形及分解后的IMF分量(b)

图 5(Fig. 5)

图 5 含噪声微地震信号去噪效果分析 a去噪前时频分布; b去噪后时频分布

微地震记录的特点是频率高、信噪比低, 因此微地震事件的自动识别和初至到时拾取对实现海量微地震数据的自动处理具有重要意义。微地震事件的自动识别和初至到时拾取主要根据地震记录中地震事件到达前后质点振动性质的差别构建特征函数来实现^[22-23]。刘玉海等^[24]利用相邻道有效信号相关性好、与随机噪声相关性差的特性来识别有效信号, 压制随机噪声。谭玉阳等^[25]针对低信噪比事件容易被遗漏的情况, 提出利用多道相似系数来检测微地震事件, 在实际资料的应用中具有较高的准确率。本文结合ICEEMD与Hilbert变换进行信号初至检测, 利用ICEEMD将信号分解为一系列IMF分量, 然后通过Hilbert变换将信号变换到频率域, 噪声与信号可以更好地分开, 进而可以更加有效地进行微地震有效信号检测。

对分解出的分量, 利用Hilbert变换获得每个点的振幅(H_i):

式中:n=1, 2, …, N, N表示分解出的IMF分量个数; A_MP表示取振幅; i=1, 2, …, I, I表示信号大小。

然后计算持续能量比(E_R1):

式中:L是τ测试点附近的能量求和窗口的长度。定义一个关于高能量地震信号的参数(E_R2):

在获得了E_R2(τ)后, 设定一个阈值, 通过给定的阈值找到E_R2的局部极大值, 大于该阈值则说明存在微地震初至信息。阈值被设置为最大峰值的分数, 定义了检测灵敏度。

图 6为选取的某地面微地震资料和利用传统的滤波方法和本文方法分别处理后的结果以及所对应的时频分布。其中图 6a是实际地震数据, 其信噪比比较低, 无法确定准确到达时间; 图 6b和图 6c分别显示了利用带通滤波与本文方法去噪后的结果。带通滤波可以去除某一频带外的噪声, 与有效信号频率接近的噪声无法完全去除; 而本文方法去除了大部分噪声, 提高了信号的信噪比, 使所选事件的信噪比更高。图 6d到图 6f给出了利用短时傅里叶变换将图 6a到图 6c对应信号变换到时频域的结果。对比时频分析结果可见, 带通滤波仅去除了特定频率范围之外的噪声, 而本文方法去除了大部分噪声, 信噪比得到了提高。图 7为去噪后的微地震数据以及对应的特征函数, 由图 7可以看出, 去噪后的有效信号(图 7a)与其特征函数(图 7b)的峰值呈一一对应关系, 可以有效地检测微地震信号初至。

图 6(Fig. 6)

图 6 对实际微地震资料利用带通滤波与本文方法去噪效果对比 a实际微地震数据; b带通滤波去噪结果; c本文方法去噪结果; d, e, f为对应a, b, c的时频分布

图 7(Fig. 7)

图 7 去噪后微地震数据(a)及对应的特征函数(b)

本文研究了基于ICEEMD的微地震噪声压制与初至检测方法, 理论模型和实际数据的应用取得了很好的效果, 得出的主要结论有:

1) ICEEMD方法可以将复杂信号有效地分解为一系列固有模态函数, 克服了EMD方法存在的模态混叠问题。阈值去噪中阈值的大小会影响去噪效果与有效信号的保护, 针对传统阈值大小不好确定的问题定义了自适应间隔阈值, 可以看做是硬阈值与软阈值的组合, 与ICEEMD结合的去噪方法可以有效地压制噪声, 数值模拟资料与实际微地震资料的应用结果表明该方法在微地震信号噪声压制方面效果明显, 提高了数据信噪比, 有较好的应用价值;

2) ICEEMD与Hilbert变换相结合的信号检测方法可以显著提高微地震事件的探测精度, 利于后续震源定位及反演。