0 引言

从微地震数据中提取微地震事件，可实时监测裂缝的空间位置与尺度，为后续压裂参数的调整与控制提供数据支撑，微地震监测对油气开发具有重要的指导意义^[1]。

当前微地震监测面临的主要问题是资料信噪比较低^[2-4]，甚至有效信号被噪声完全掩盖。这对以初至拾取为基础的震源定位、压裂预测、机制分析等后续工作造成了严重影响。因此，针对低信噪比或超低信噪比的微地震资料，选择合理有效的初至拾取方法是微地震资料处理的关键。

通常，初至拾取以信号与噪声在振幅、频率、偏振等方面的差异为基础。常采用长短时窗平均能量比(STA/LTA)法^[5-7]、基于自回归模型的Aksike信息准则(AIC)法^[8]、相关法^[9-12]、偏振分析法^[13]、二分法^[14]、分维相关法、时频分析法、神经网络法^[15]、独立分量法^[16]、基于高阶统计量的PAI-S/K法^[17]以及基于小波分解的高阶统计初至拾取方法^[18]。Coppens^[19]对不同时窗内不同炮检距道集的能量进行统计相关分析，提出一种自动定位初至点的方法，有效提高了拾取精度；邹红星^[20]提出一种基于自回归模型的判读法，该方法具有一定自适应能力，可对不同噪声影响程度的接受信息进行初至拾取。贾瑞生等^[21]、蔡剑华等^[22]基于经验模态分解(EMD)和本征模态函数(IMF)重构对低信噪比信号去噪，然后应用正交性指标与能量保存度指标对传统EMD中的端点效应进行压制，从而实现震相初至的拾取，其应用效果优于常规基于AIC准则的拾取方法。谭玉阳等^[23]提出了基于信号信噪比构造不同检测函数的初至拾取算法(SLPEA)，通过综合考虑有效信号与环境噪声间的多种特征差异，提高拾取过程的抗噪能力，但时窗的选择依据并不明确，需依赖经验或反复调试。

在前人的基础上，本文基于三分量间有效成分的特征相似性，采用CEEMDAN算法对各分量信号进行分解，然后用各阶IMF构建三分量IMF矩阵，再进行主成分分析(PCA)，提取主成分特征，达到对低信噪比微地震资料有效波初至拾取的目的。通过对不同信噪比合成信号与实测信号的处理，证实该方法具有很强的抗噪能力，适用于有效信号被噪声掩盖的低信噪比微地震资料。

1 方法原理 1.1 CEEMDAN算法原理

EMD可将信号x(t)分解为一系列IMF，但信号需要满足两个条件：①极值点数与过零点数相等或相差不超过一个；②上、下包络线的均值在各处均为0。

对于实际非平稳时变信号，常规EMD会受到模态混叠的影响。通过引入正态分布的白噪声，形成噪声辅助信号处理方法(NADA)，即为常用的EEMD^[24]。而本文采用基于改进EEMD方法——CEEMDAN算法，通过在原始信号上加入有限方差约束的多组独立同分布自适应白噪声，得到含噪信号的集合。该改进方法可进一步降低迭代次数，压缩频率混叠区域，提高收敛性能，对于非平稳信号不同频率成分具有更高的分辨能力。

对信号x(n)的CEEMDAN的实现步骤如下。

(1) 生成含噪信号集

$ {x^i}\left( n \right) = x\left( n \right) + {w^i}\left( n \right) $

(1)

式中wⁱ(n)(i=1，2，…，I)为满足高斯分布的噪声，I为集合样本数；

(2) 对xⁱ(n)进行EMD，得到各样本的一阶IMF F₁ⁱ，再计算其均值作为x(n)的一阶IMF

$ {{\tilde F}_1}\left( n \right) = \frac{1}{I}\sum\limits_{i = 1}^I {F_1^i} $

(2)

(3) 计算一阶残差量

$ {r_1}\left( n \right) = x\left( n \right) - {{\tilde F}_1}\left( n \right) $

(3)

(4) 计算二阶IMF

$ {{\tilde F}_2}\left( n \right) = \frac{1}{I}\sum\limits_{i = 1}^I {{E_1}\left\{ {{r_1}\left( n \right) + {\varepsilon _1}{E_1}\left[ {{w^i}\left( n \right)} \right]} \right\}} $

(4)

式中：E_j(·)表示信号的j阶IMF；ε_j为控制白噪能量的参数，Wu等^[24]认为对于主频较高的信号，取值应较小，反之较大，本文采用较小常量进行试算。

(5) 对于k(k=2，3，…，K，K为设置的最高IMF阶次)阶分量，计算k阶残差

$ {r_k}\left( n \right) = {r_{k - 1}}\left( n \right) - {{\tilde F}_k}\left( n \right) $

(5)

(6) 计算k+1阶分量

$ {{\tilde F}_{k + 1}}\left( n \right) = \frac{1}{I}\sum\limits_{i = 1}^I {{E_1}\left\{ {{r_k}\left( n \right) + {\varepsilon _k}{E_k}\left[ {{w^i}\left( n \right)} \right]} \right\}} $

(6)

(7) 重复步骤(5)、步骤(6)直到残差不可再分解或达到最高IMF阶次。

最终残差满足

$ R\left( n \right) = x\left( n \right) - \sum\limits_{k = 1}^K {{{\tilde F}_k}\left( n \right)} $

(7)

信号可表示为

$ x\left( n \right) = \sum\limits_{k = 1}^K {{{\tilde F}_k}\left( n \right)} + R\left( n \right) $

(8)

与传统EEMD方法相比，CEEMDAN算法最终得到的各阶IMF分量对原始信号频率成分的表征不同。如图 1所示，对3~6阶IMF分量进行频谱对比，可见传统EEMD各分量间能量差异较大，频率混置区较大。而CEEMDAN方法得到的各IMF分量间能量均衡，频率混叠区较窄，对于非平稳信号，不同频率成分具有更高的分辨能力。通过对心电图信号(图 1a)的测试考察其计算性能，该信号具有非平衡性、非线性与随机性，具有丰富的频率成分，适合用于时频分析方法测试。该测试信号的各阶分量的迭代次数箱形图如图 2所示，可见CEEMDAN各阶分量的迭代次数均显著降低。对于主要承载有效信息的3阶至6阶IMF分量，CEEMDAN方法的计算量约为EEMD方法的1/3。

图 3为实测微地震信号及其频谱，信号采样间隔为0.25ms。图 4为实际微地震数据传统EMD方法分解的6~9阶IMF及其频谱。各阶IMF在靠近端点处产生较大幅度的振荡，并且在6阶和7阶IMF对应的频谱上均显示了两个明显峰值(图 4a、图 4b中的蓝色区域)，明显包含两组频率成分，表明存在明显的模态混叠现象^[25-26]。

图 5为实际微地震数据CEEMDAN方法分解的6~9阶IMF分量及其频谱，端点附近的异常振荡得到有效压制，且各阶IMF对应的频谱均无明显的多峰现象，模态混叠现象得到有效消除。这是因为CEEMDAN方法利用高斯白噪声频率均匀分布的统计特征，向原始信号中加入不同的白噪声，使信号在不同尺度上具有连续性，并通过在分解时的每一阶段添加特定白噪声，并计算一个唯一残差以得到每个IMF分量，使原始信号的模态频谱分离效果更好，端点效果的影响显著降低。

1.2 主成分分析与特征提取

通常采用三分量检波器采集微地震数据，有效信号在各分量上能量不同，但初至点一致、有效成分波形高度相似，故可通过对三分量信号先进行PCA，再进行特征提取。微地震信号中常出现的噪声为近似满足独立高斯分布的自然噪声与低频震荡。通过CEEMDAN分解，不同阶次的IMF承载原始信号中的不同信息。由于原信号中噪声能量强于有效成分，故在各阶IMF中，仍存在噪声能量，可采用PCA对三分量IMF矩阵提取主成分。三分量数据的k阶IMF矩阵可表示为

$ \mathit{\boldsymbol{M}}\left( k \right) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{D}}_x^k}\\ {\mathit{\boldsymbol{D}}_y^k}\\ {\mathit{\boldsymbol{D}}_z^k} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {D_x^k\left( 1 \right)}&{D_x^k\left( 2 \right)}& \cdots &{D_x^k\left( N \right)}\\ {D_y^k\left( 1 \right)}&{D_y^k\left( 2 \right)}& \cdots &{D_y^k\left( N \right)}\\ {D_z^k\left( 1 \right)}&{D_z^k\left( 2 \right)}& \cdots &{D_z^k\left( N \right)} \end{array}} \right] $

(9)

式中：D_x^k、D_y^k、D_z^k分别为x、y、z分量的k阶IMF；N为各阶IMF样点数。计算构建的三分量IMF矩阵M(k)的特征向量W，可测得矩阵在各特征向量上的能量分配

$ {E^k} = \mathit{\boldsymbol{M}}\left( k \right) \cdot \mathit{\boldsymbol{W}} $

(10)

其中第一主成分对应的特征向量，应该满足

$ {\mathit{\boldsymbol{W}}_1} = \arg {\max _{\left\| \mathit{\boldsymbol{W}} \right\| = 1}}\left\{ {\sum\limits_k {\left[ {\mathit{\boldsymbol{M}}\left( k \right) \cdot \mathit{\boldsymbol{W}}} \right]} } \right\} $

(11)

写成矩阵形式

$ {\mathit{\boldsymbol{W}}_1} = \arg {\max _{\left\| \mathit{\boldsymbol{W}} \right\| = 1}}\left[ {{{\left( {{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{W}}} \right)}^{ - 1}}{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\left( k \right)\mathit{\boldsymbol{M}}\left( k \right)\mathit{\boldsymbol{W}}} \right] $

(12)

W₁使第一主成分具有最高能量，第一主成分能量可表示为

$ E_1^k = \mathit{\boldsymbol{M}}\left( k \right) \cdot {\mathit{\boldsymbol{W}}_1} $

(13)

将求取的主成分从原矩阵M(k)中剔除，得到剩余矩阵$ \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over M} }}$

$ \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over M} }} = \mathit{\boldsymbol{M}}\left( k \right) - \mathit{\boldsymbol{M}}\left( k \right){\mathit{\boldsymbol{W}}_1}\mathit{\boldsymbol{W}}_1^{\rm{T}} $

(14)

令M(k)=$ \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over M} }}$，$ \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}$_i=W₁(i=1，2，…，S)，重复式(11)~式(14)，逐次从IMF分量矩阵M(k)中减去当前最大主成分分量，直到求出第S个主成分，并依次可求取各阶主成分能量

$ E_S^k = \mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over M} }} \cdot {{\mathit{\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over W} }}}_i} $

(15)

通过对相邻三道三分量微地震记录的各阶IMF进行PCA，统计各主成分所占的能量百分比(图 6)，可见4、5、6阶IMF的第一主成分具有明显的能量优势。其原因在于有效成分能量主要落于这三阶IMF的第一主成分中，而其他各阶IMF中，前三项主成分能量差异不大，主要由高能量的噪声构成。通过将高能量噪声等特征较弱的成分分配到各阶IMF中，而特征明显的有效成分能量则集中分配于前几项主成分中，剔除后几项主成分，可对有效成分进行提取。

1.3 IMF主成分加权重构

对三分量数据进行CEEMD后，对各阶IMF进行PCA，选取目标主成分，最后加权重构。通过CEEMD分解得到的各阶IMF，反映出信号的不同成分。以图 7含噪信号为例，可以看出，1、2阶IMF主要反映信号的大体趋势，7~9阶IMF则主要承载高频信息。在低信噪比情况下，高频信息主要为噪声，那么有效信号的时域特征则主要由中间的3~6阶IMF反映。因此，重构时，可将1、2阶和7~9阶IMF的权重减少，将3~6阶IMF的权重加大，以便更加突出有效信号的特征。

对于重构时各级主成分权重系数的选取，本文的思路如下。

(1) 归一化处理：将各阶IMF主成分单独重构回时域波形，对时域波形进行振幅归一化处理；

(2) 计算各阶时域波形与原始信号所提取的子波间的相关系数；

(3) 以相关系数绝对值作为各阶IMF主成分的重构权值。

1.4 处理流程

针对低信噪比条件下的微地震信号，首先对原始信号进行CEEMDAN；然后基于同一检波器中三个分量间的有效成分相关性，在IMF域内从高能量噪声中提取出有效成分特征，并剔除其余成分；最后将剩余的IMF主成分加权重构，得到体现有效信号特征的时域波形。具体处理流程如下：

(1) 对三分量微地震资料进行CEEMDAN；

(2) 对三分量中各阶IMF分量进行PCA，以75%为门槛值，剔除次要成分；

(3) 按不同权值进行IMF主成分重构，得到主成分在时域内波形。

2 模型测试与实际应用

采用两组模型对本文方法进行可行性测试与适应性测试。模型Ⅰ采用合成信号与白噪声，测试该方法对噪声掩盖下的具有内在固有特征成分检测能力。模型Ⅱ采用实际有效信号，叠加白噪声，通过调整噪声能量，测试该方法适用的信噪比下限。通过对2000组信号测试，对比基于本文方法与常规基于AIC信息准则的长短视时窗(STA/LTA-AIC)方法的拾取误差。最后对具有不同品质的低信噪比实测资料进行了处理。

2.1 模型Ⅰ

设计模型Ⅰ合成信号为

$ x\left( t \right) = \lambda y\left( t \right) \otimes \left[ {\exp \left( { - 4.5t} \right)\sin \left( {100{\rm{ \mathsf{ π} }}t} \right)} \right] + w\left( t \right) $

(16)

式中：y(t)为标准Ricker子波；“⊗”为褶积运算符；w(t)为白噪声；λ为振幅因子，控制合成信号的信噪比。初至时间设置为500ms，采样间隔为0.25ms，用式(16)生成能量随时间呈指数衰减的正弦信号与标准Ricker子波的褶积作为有效信号(图 8中红线)，调整λ得到一系列不同信噪比(SNR)的测试信号(图 8中黑线)。

应用本文方法对测试信号进行处理，结果如图 9所示。图 9a中由于信噪比过低，在500~650ms可以看到持续的低能量伪响应，随着SNR的提高，该伪响应的相对能量降低。经过本文方法处理后，初至事件的响应清晰度与分辨率都有显著提升，应用常规方法在该结果上均能得到高精度的初至拾取结果。需要指明，有效事件的检测响应，需要一致地出现在三分量信号中的至少两个分量上，否则不应被视为微地震事件的有效响应。如图 9e中所示，在约210ms处出现单一响应，而在其他分量上则并没有同时出现同步响应，因此将其作为伪响应剔除。该模型测试结果表明，本文方法最低可适用于SNR=-20dB的信号。当SNR低于-20dB，有效信号因能量过低而无法识别。

2.2 模型Ⅱ

江汉油田水力压裂实验中采集到高信噪比微地震信号(图 10a)，采样频率为4000Hz，有效频带为50~220Hz。利用常规方法(如STA/LTA-AIC等)拾取的P波初至为287.50ms，S波初至为527.25ms。加入不同能量的白噪声，获得模型Ⅱ不同SNR的合成数据(图 10b~图 10f)，当SNR降至-6.9787dB时(图 10c)，STA/LTA方法无法有效拾取初至。不同SNR的测试信号频谱如图 11所示，随着噪声能量的加强，有效频段的能量优势逐渐消失。本文方法对不同SNR的合成数据(图 10b~图 10f)的处理结果如图 10g~图 10k所示，可以看出，即使SNR降到-19.0199dB，处理后仍可对微地震事件进行有效检测。

通过对2000组不同信噪比信号的测试，本文方法与常规STA/LTA-AIC法的拾取精度误差统计结果如图 12所示。可见本文方法的拾取精度在SNR为-2~-20dB时，远优于常规STA/LTA-AIC方法。

2.3 实际资料应用

实测井下压裂微地震数据(图 13a)采样频率为4000Hz，部分记录数据品质较差，常规方法拾取效果不理想。图 13a中蓝色显示的记录上可肉眼观测到疑似P波与S波初至波峰，约在1080ms与1130ms处(激发时刻不为0)，在红色显示的纪录上仅在1230ms附近可见单一疑似初至波峰，而在黑色显示的纪录上，无任何突变波形体现。本文方法处理结果如图 13b所示，蓝、红色记录均可见明确的检测响应。而对于黑色记录，本文方法没有得到有效的处理结果，原因在于有效成分能量过低，处理后仍被噪声残留成分所淹没。在蓝、红色记录的检测结果基础上，采用常规STA/LTA-AIC方法拾取初至，以图 13b中不同底色(蓝—黄、黄—灰)的分界线标识P波、S波初始时刻。

3 结论

本文针对噪声能量较大、有效信号被淹没的初至拾取难题，提出在CEEMDAN的基础上，利用微地震记录中三分量信号有效成分具有高度相似性，对三分量信号的各阶IMF进行主成分分析，保留有效成分的特征，剔除干扰成分，有效实现了低信噪比资料中的初至特征检测。实际资料的处理结果表明该方法对于在有效信号被噪声掩盖的情况下仍能有效、精确地拾取初至。