微地震监测技术作为获得地层裂缝分布的一种较为有效方法, 已经在油气田开发与煤矿开采等领域得到了广泛应用^[1]。由于微地震波能量微弱, 检波器采集到的微地震信号存在大量随机噪声, 因而后续的初至拾取、震源定位及裂缝解释工作会受到很大影响^[2-3], 因此, 需要对微地震信号进行滤波、去噪等处理, 将微地震有效信号从噪声中分离出来。

传统的傅里叶变换是信号去噪的有效方法, 但对于具有非平稳信号特点的微地震信号去噪效果不佳^[4-5]。连续小波变换在低频部分具有较高的频率分辨率, 在高频部分具有较高的时间分辨率, 能够在时频域很好地区分非平稳信号的突变部分, 但去噪结果受到Heisenberg不确定性和小波基函数选取的影响^[6-9]。经验模态分解(EMD)是一种信号自适应分解方法, 已经在生物医学的信号处理、机械设备的检测与故障诊断、地震勘探的信号处理^[10-14]等领域得到了广泛的应用, 但其稳定性会受到模态混叠现象的影响^[15-16]。

针对EMD中出现的模态混叠现象, WU等^[15]在信号中加入功率谱密度均匀分布的白噪声, 使其在不同尺度上具有连续性的特征, 给出了总体平均经验模态分解方法; YEH等^[16]通过加入正、负成对的辅助噪声, 提出了互补集合经验模态分解方法。上述方法能在一定程度上抑制模态混叠现象, 但计算的复杂程度会显著上升。

本文根据微地震信号具有随机性、非平稳性与时频耦合的特点, 兼顾EMD出现的模态混叠问题, 提出了基于EMD互信息熵与同步压缩变换(SST)的微地震信号去噪方法。模型测试和实际资料处理验证了方法的有效性与实用性。

1 方法原理

在实际微地震信号去噪中, 有效信号往往位于低频部分, 所以通常把高频成分的固有模态函数(IMF)分量作为噪声直接去除, 重构原信号剩余的IMF分量即可实现压制噪声的目的^{[11, 17-18]}。但是, 简单地舍弃高频成分会丢失部分有效信息, 造成去噪效果不佳, 因此需要发展一种有效剔除混叠噪声并保留有效信号的方法。

1.1 IMF互信息熵成分辨识

EMD算法将时间序列信号自适应地分解成一组含有不同尺度的IMF分量, 能够较好地处理非平稳信号, 它不依赖于基函数的选取。EMD得到的IMF必须满足以下两个条件^[11-16]:①IMF中极值点的个数与过零点的个数相等或最多相差一个; ②在任意数值点处, 由局部极大值与极小值确定的包络线均值为0。

经过EMD后, 微地震信号s(t)可以表示为：

$ s\left( t \right) = \sum\limits_{i = 1}^n {{S_{{\rm{IM}}{{\rm{F}}_i}}}\left( t \right) + {r_n}\left( t \right)} $

(1)

式中:S_IMFi(t)表示频率从高到低排列的IMF分量; r_n(t)表示残差分量, 代表信号的平均趋势。将公式(1)改写为：

$ s\left( t \right) = H\left( t \right) + L\left( t \right) $

(2)

式中: $H\left( t \right) = \sum\limits_{i = 1}^{k - 1} {{S_{{\rm{IM}}{{\rm{F}}_i}}}} \left( t \right)$ 表示高频成分的IMF分量组合; $L\left( t \right) = \sum\limits_{i = k}^n {{S_{{\rm{IM}}{{\rm{F}}_i}}}} \left( t \right) + {r_n}\left( t \right)$ 表示低频成分的IMF分量与残差分量的组合。

互信息熵用于衡量两个随机变量之间的统计依存性^[18-21], 其表达式为：

$ I\left( {X,Y} \right) = \sum\limits_{i = 1}^r {\sum\limits_{j = 1}^s {p\left( {{x_i},{y_i}} \right){\rm{lb}}\frac{{p\left( {{x_i},{y_i}} \right)}}{{p\left( {{x_i}} \right)p\left( {{y_i}} \right)}}} } $

(3)

式中:p(x_i, y_j)为联合概率分布; p(x_i)和p(y_j)为边缘概率分布; X和Y表示不同的IMF分量; r和s分别是X和Y的符号数量。

假设高频成分为噪声干扰, 低频成分为有效信号, 利用IMF分量之间的互信息熵关系可以辨识出高频与低频成分的分界。由EMD得到的IMF分量特征可知：高频噪声对每一个IMF分量的依存性逐渐降低, 而低频有效信号对每一个IMF分量的依存性逐渐增强。因此, 可以假设高频成分与低频成分是部分相互统计独立的。由互信息熵的特点可知：两个相互独立的随机变量之间的互信息熵应等于0^{[18, 22]}。所以, 在计算每个相邻IMF分量之间的互信息熵时会出现局部极小值, 只需搜索第一个局部极小值即得到高频成分与低频成分的分界。于是可得以下搜索目标函数：

$ k = {\rm{first}}\left\{ {\mathop {\min }\limits_{1 \le i \le n - 1} \left[ {I\left( {{S_{{\rm{IM}}{{\rm{F}}_i}}},{S_{{\rm{IM}}{{\rm{F}}_{i + 1}}}}} \right)} \right]} \right\} $

(4)

式中:k为高频成分与低频成分分界处的IMF分量序号。

1.2 基于同步压缩变换的混叠噪声压制改进算法

SST算法是在小波变换的基础上发展起来的一种时频域重排算法, 但与原始重排算法不同的是, 其在提高时频分辨率的同时, 也能够支持信号的重构^[23-27]。由于SST能够获得较高频率精度的时频谱, 使得各频率曲线之间不存在交叉项, 因此可以较好地解决频率的混叠问题。利用SST实现高频成分混叠噪声压制的步骤如下：

1) 对高频成分H(t)进行连续小波变换, 得到小波变换系数W(ω_l, b), 并计算其瞬时频率ω(a, b)。

2) 利用公式(5)计算自适应阈值γ, 进行小波阈值去噪 $\left( {\hat a = \Delta t/{2^{\hat n}}} \right)$ , 尺度因子 ${\hat n}$ 取32或者64最佳^[23-24]。

$ \gamma = \frac{{\sqrt {2\lg n} \cdot {\rm{median}}\left( {\left| {W\left( {\hat a,b} \right)} \right|} \right)}}{{0.6745}} $

(5)

式中:n表示信号采样点数。

3) 利用计算得到的瞬时频率建立(b, a)→(b, ω(a, b))之间的映射关系, 对时间-尺度平面的能量进行重新分配, 将其转化为时间-频率平面, 使得中心频率ω_l相邻区间[ω_l－Δω/2, ω_l+Δω/2]的值被压缩到ω_l上^[28], 获得SST系数T(ω_l, b), 其离散形式为

$ {T_{\rm{s}}}\left( {{\omega _l},b} \right) = \frac{1}{{\Delta \omega }}\sum\limits_{{a_k}:\left| {\omega \left( {a,b} \right) - {w_l}} \right| \le \Delta \omega /2} {{W_{\rm{s}}}\left( {{a_k},b} \right){a^{ - 3/2}}\Delta {a_k}} $

(6)

4) 利用SST获得的高分辨率时频谱确定有效信号的时间与频率范围, 再对其进行SST重构, 获得信号x(t), 实现对高频成分混叠噪声的压制。

$ x\left( t \right) = 2{\mathop{\rm Re}\nolimits} \left\{ {\sum\limits_l {\left[ {{T_{\rm{s}}}\left( {{\omega _l},b} \right)\Delta \omega } \right]} } \right\}/C $

(7)

式中: $C = \int_0^\infty {\frac{{{\mathit{\Psi }^*}\left( \xi \right)}}{\xi }} {\rm{d}}\xi ,{\mathit{\Psi }^*}$ 表示母小波函数的共轭傅里叶变换。

将高频成分混叠噪声压制后的信号x(t)与剩余的低频成分有效信号进行重构, 即可实现对微地震信号噪声的去除。

2 测试与应用 2.1 模型测试

为了验证本文去噪方法的可行性, 采用如图 1所示的井中微地震观测系统, 利用有限差分正演方法模拟了一个井中16级检波器接收的微地震记录(图 2), 数据采样间隔为0.5ms, 采样点数为1000。将图 2模拟的微地震记录加入高斯白噪声, 加噪后数据的信噪比为0.3, 如图 3所示, 从中抽取一道数据进行去噪处理。图 4为抽取的单道记录加噪前、后波形及其频谱, 对单道加噪信号进行自适应EMD, 得到8个IMF分量以及对应的频谱, 如图 5所示。可以看出每一个IMF分量都有不同的振幅, 而且分解的频率按照从高到低的顺序排列。

利用互信息熵公式(3)计算相邻IMF分量之间的互信息量, 结果如表 1所示。按照本文1.1节所给的高频部分与低频部分分界点判定方法, 可知IMF34为出现的第一个极小值, 所以IMF4即为加噪信号的高频部分与低频部分分界点, 图 5中IMF4的波形以及对应的频谱也验证了判定结果的正确性。图 6为将IMF4之后的分量叠加获得的低频部分信号和SST的结果, 从中可以看出低频部分的噪声含量较低, 可以将低频部分当作有效信号进行处理。

利用1.2节给出的SST混叠噪声压制方法, 对IMF4之前的分量叠加组成的高频部分进行同步压缩变换, 去除混叠噪声, 如图 7所示。图 7a为IMF4之前的分量叠加组成的高频部分信号, 从中可以看出噪声基本上分布在高频部分。图 7b和图 7c分别为高频部分短时傅里叶变换(SIFT)和SST时频谱, 可以看出, SIFT时频谱模糊, 不能很好地分辨出有效信号, 而SST时频谱具有较高的频率精度, 各时频曲线间不存在交叉项, 提高了有效信号脊线的分辨率, 也证明了SST能够较好地改善频率混叠现象。图 7d为SST自适应阈值去噪后的时频谱, 可以看出高频部分的噪声得到了很好的压制。图 7e为对高频部分有效信号进行提取后得到的有效信号, 将其与剩余的低频部分进行重构, 即可得到降噪后的信号。

图 8对比了单道加噪信号直接舍弃高频分量的去噪结果和本文方法的去噪结果, 图 9对比了图 3所示的微地震加噪记录直接舍弃高频分量的去噪结果和本文方法的去噪结果。可以看出, 直接舍弃高频分量去噪的方法虽然能够压制噪声, 但同时使有效信号频率出现部分缺失, 并且还混叠了部分噪声, 造成了有效信号的波形畸变与衰减; 而本文方法不仅很好地压制了噪声, 而且使有效信号的波形得到了较好保持。定义去噪后信号与原信号的差值为含有的噪声, 信号与噪声能量的比值为信噪比R_SN^{[27, 29]}：

$ {R_{{\rm{SN}}}} = \frac{{\left\| M \right\|_2^2}}{{\left\| {S - M} \right\|_2^2}} $

(8)

式中:M表示原信号; S为去噪后信号。由公式(8)计算直接舍弃高频分量方法去噪后的信噪比为2.43, 本文方法去噪后的信噪比为6.67, 定量说明了本文方法去噪的有效性。

不断地改变含噪信号的信噪比, 观察不同噪声水平下两种方法的去噪效果, 并计算去噪后的信噪比, 结果如表 2所示。可以看出, 在不同噪声水平下, 两种方法都能提高含噪信号的信噪比, 但是本文方法具有更好的去噪效果。

2.2 应用实例

本文截取某油田压裂的一段井中微地震有效事件数据进行了去噪效果分析, 如图 10所示。原始数据由16级三分量检波器接收, 采样间隔为0.25ms, 采样点数为1000。图 10a中1~16, 17~32, 33~48道分别为数据的z, x, y分量, 由于数据的信噪比较低, 噪声几乎淹没了有用的微地震信号, 无法清晰地识别出P波、S波。图 10b为对原始数据进行EMD后直接去除高频成分的结果, 可以看出, 虽然高频噪声被滤除, 但是微地震有效信号同相轴仍然不够清晰, 有效信号波形发生了畸变。图 10c为利用SST方法^[27]提取的结果, 可以看出, SST方法虽然较好地提取了有效信号, 但是混叠了部分噪声, 说明噪声压制还不够彻底。图 10d为利用本文方法去噪后的结果, 原始数据的噪声得到了很好压制, 有效信号的波形也得到了很好保持, 为后续的初至拾取、震源定位及裂缝解释工作奠定了基础。

3 结束语

本文针对微地震信号具有随机性、非平稳性与时频耦合的特点, 以及EMD存在的模态混叠问题, 利用同步压缩变换能提高信号时频分辨率, 同时支持信号重构的优点, 提出了基于EMD互信息熵与同步压缩变换的微地震信号去噪方法。不同噪声强度的理论模型测试和实际资料分析表明, 该方法能够很好地去除信号中的混叠噪声, 较好地将有效信号从噪声中提取出来, 提高资料信噪比, 为后续的初至拾取、震源定位及裂缝解释等工作奠定了基础。