梁 喆，彭苏萍，郑 晶

(1.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室，北京 100083;2.安徽理工大学 电气与信息工程学院，安徽 淮南 232001)

微震监测技术在矿业工程、石油开采和防灾减灾监测方面已经获得了广泛的应用。微震事件发生时，采集仪器得到的是混有背景噪声的信号。带有噪声的微震信号对定位精度有很大的影响，同时也对后续提取特征参数带来误差。为了能精确定位微震事件和计算信号特征参量，必须对微震监测系统记录到的信号进行滤波、去噪。

李夕兵等［1－2］采用EMD方法对震动信号进行滤波、消噪。虽然EMD变换非常适合处理非平稳信号，但运用EMD方法对非平稳信号进行分解时，存在端点效应问题，造成整个数据序列分解后的结果严重失真，由此得到的分析结果的精度较低。

针对EMD的端点效应问题和实际应用，许多专家学者提出了不同的延拓方法，例如极值点镜像延拓、包络线延拓、局部波形延拓、波形匹配度延拓、支持向量机延拓、神经网络延拓、AR模型延拓、最大Lyapunov指数边界延拓方法，基于时间尺度的LS-SVM端点延拓以及最大相关波形延拓等［3－9］。以上延拓方法，对短序列数据都有很好的延拓精度，但是使用时需要大量的计算时间。微震信号是一种突变的非平稳信号，其自身是一个不确定的信号，没有规律可循。对微震数据处理时，数据量大，且要求实时处理，在处理、解释的过程中，检测出突变发生的位置非常重要。采用上述方法对微震信号EMD分解时的端点效应处理时导致边界误差大，造成分解后突变信号衰减或丢失。

自适应滤波具有在未知环境下良好运行并跟踪输入统计量随时间变化的能力，在没有任何关于信号和噪声的先验知识的条件下，自适应滤波器利用前一时刻已获得的滤波器参数来自动调节现时刻的滤波器参数，以适应信号和噪声未知或随机变化的统计特性。

结合 LMS自适应滤波器理论，本文对提出了LEMD(Empirical Mode Decomposition Based on LMS)端点效应处理方法。所提出的LEMD方法是利用LMS自适应滤波器根据端点以内数据所包含的信息，在数据序列两端分别延拓出一段波形，使其符合数据的自然趋势，从而产生EMD分解时边界所需要的极大值和极小值。自适应滤波器通过自动地调节现时刻滤波器参数，适应微震信号和随机噪声随时间变化的统计特性，同时根据数据段内部的信息和端点处数据的信息，实时修正延拓的结果。在整个EMD分解的过程中，实时跟踪边界信息，采用边延拓边分解的方法，不断修正EMD筛分时插值算法带来的误差，最终实现最优延拓。

将LEMD方法应用于微震信号的降噪，成功获取了微震信号定位和特征参数等信息，结果表明该方法是合理、有效的。

1 EMD分解及端点效应分析

在Hilbert-Huang中，为了计算瞬时频率，定义了本征模态函数(IMF)，它必须满足两个条件:①在整个数据段内，其极值点的个数和过零点的个数必须相等或相差一个;②在任意时刻，由局部极大值点和极小值点形成的上、下包络线平均值为零。对于不满足IMF条件的复杂信号，可以用 EMD方法对其进行分解［10－11］。

EMD分解的基本思想是:对一给定信号，先获得信号极值点，然后通过插值获得信号上下包络，计算出包络均值，信号每减去1次包络均值，便可得到1个固有模态函数。如此重复，直到将信号分解成有限个基本模式分量IMF和残余项rn(t)的组合。

通过经验模态分解，信号x(t)可表示为:

式中，imfi(t)表示具有不同频率分量的IMF分量;rn称为残余函数，代表信号的平均趋势。

对EMD分解后的不同IMF分量的取舍可以分别实现对信号的低通、带通及高通滤波的功能［12］，根据式(1)可知，忽略趋势项rn(t)，有

式中，x(t)为原始信号，X(t)为去噪后的信号。由式(2)可知，从原始信号中除去噪声信号，就完成了噪声信号的分离。

EMD是基于数据极值点进行分解的，对数据进行分析的时候，由于时间窗的选取，分析的数据往往是被截断后的数据，这就造成信号在端点处不一定都是极值点，更不可能同时为极大值点和极小值点，致使利用样条曲线插值求包络运算的时候，在端点处形成较大的摆动，出现误差，随着分解过程的进行，重复进行插值求包络线，导致误差不断累积并向内传播污染数据，从而影响分解结果，这就是 EMD方法的端点效应问题。

微震信号的主频率一般在几赫兹到一百赫兹左右，信号在传播的过程中，经过不同介质的吸收，其频率总是在低频段分布。对于处于低频段微震信号，EMD分解时其极值点之间的距离较大，即时间跨度大，端部的边缘效应更容易影响到信号内部，所以在用EMD对微震信号进行处理时必须抑制端点效应。

2 改善端点效应新方法

2.1 自适应滤波器(预测算法)基本原理［13］

在目前的研究及应用领域中，自适应滤波器可分为四类:系统识别、系统求逆、信号预测和干扰抵消。

假设输入为某一特定矢量

自适应算法即根据第k个采样周期所获得的误差e(kT)，计算下一个采样周期计算式中的各个权系数:

式中，μ为步长因子。式(5)与式(6)的算法称为LMS算法。这种算法对于每一个输入样本，只需对其进行式(7)中的两个乘法与两个加法运算，因此该算法易于用实时系统实现。

2.2 利用自适应预测算法对信号进行延拓

自适应滤波器的权系数向量迭代更新的过程，实质上是一种自学习和自训练的过程。当使用LMS自适应滤波器进行端点延拓时，如果直接进行延拓，权系数向量会影响延拓的结果，造成结果不准确。在延拓之前可以首先对滤波器进行训练，对权系数向量反复地优化，增强其在动态条件下的适应性。训练结束后将训练得到的权系数作为延拓的初始值，从端点处开始进行延拓。具体方法为:设一个数据序列内有n个数据{x(1)，x(2)，…，x(n)］，首先根据待延拓的数据点数m产生学习样本，根据需要的目标输出利用式(5)～式(7)对自适应滤波器进行训练，得到对应的权系数Wj。然后用训练好的权系数Wj作为延拓时候的初始值，将端点x(n)前的n－1个数据代入式(5)计算出端点处的第1个延拓值x^(n+1)。并以x^(n+1)为原数据序列新的端点，利用x^(n+1)为端点构成新的数据n个数据组合{x(2)，x(3)，…，x(m)，x^(n+1)］ 可以得到第2个数据序列延拓值x^(n+2)。以此类推，根据所需延拓的数据个数m可得到全部延拓序列x^(n+j)(j=1，2，…，m)。最后将新的延拓数据与原数据序列连接起来，就可以得到延拓后的全部数据序列。实际使用时数据点延拓个数可以根据需要的极值数确定，由延拓算法自适应决定。

综上分析，LEMD算法的具体步骤如下:

(1)输入数据序列;

(2)根据数据序列的长度，选取一部分数据按式(5)～式(7)对自适应滤波器训练，训练结束的自适应滤波器在数据序列的两端延拓一段波形，波形的长度自适应决定，停止条件为两边各延拓出2个极大值和2个极小值;

(3)EMD分解，得到第一层分量，保留剩余信号中间对应延拓前的部分，采用自适应滤波器再次进行延拓，停止条件和第一层相同;

(4)每层分解结束后均对剩余信号再次延拓，直到整个分解结束。

LEMD算法采用延拓的数据完成EMD分解，信号包络插值计算采用延拓后的全部数据，而在EMD分解完成后，舍弃两端的延拓数据，只保留中间对应延拓前的部分，作为最终分解结果。

为了便于分析比较，采用文献［14］中的仿真信号作为实例来分析，该仿真信号为:

原始信号和上下包络线如图1所示。为了和文献［14］中的结果对比，采样频率设为200 Hz，数据点数为N=181。图1中可见，该数据序列仅存在3个极大值点和4个极小值点，包络线不仅没有包含所有的数据，且上包络线在数据的两端出现了严重的失真。因此在采用EMD方法对信号进行分解时必须对端点进行处理。采用本文的方法，在信号的两边各延拓40个数据。首先对自适应滤波器进行训练。以右边的延拓为例:用信号的第1～160点的值按顺序分为4组，每组40个数据，将1、3组作为学习样本，2、4组作为期望输出。也就是1组作为训练输入，2组作为期望输出，通过构造的学习样本和期望输出对自适应滤波器训练，训练结束后，利用得到的权系数，将信号的第82～121点的值代入式(4)进行预测运算，得到第122点的值，再利用第83～122点的值预测得到第123点的值，以此类推预测得到第161点的值，即得到40个数据点。同样，在左边也可以得到40个数据点。延拓的波形和包络线如图2所示。

图1 原始信号和包络线Fig.1 original signal and the envelope

图2 延拓后波形和包络线Fig.2 Waveform and envelope after extension

图2 (a)是对图1所示的信号采用LEMD方法在端点处延拓后得到的信号，其中实线是对数据延拓后的信号波形，虚线为真实信号波形。虽然延拓的波形有稍许畸变，但和真实值相比，其极值的大小非常接近，误差较小。图2(b)为延拓后的包络线，可以看出本文方法有效的抑制了包络线在端点处的发散，所有信息都被包含在包络线内。

为了比较LEMD方法的特点，将文献［14］中的的几种延拓方法性能比较的结果一并列出。从表1中可见，LEMD端点延拓方法与支持向量机、BP神经网络相比，LEMD方法的学习时间仅需要0.025 8 s，误差的均方值为0.143，LEMD方法的学习时间明显比支持向量机和BP神经网络都短很多。LEMD方法延拓的信号与原始信号误差的均方值要比上述2种方法的稍大，但是本文采用的是边分解边延拓的方法，误差不会累积，因此对数据精度没有太高的影响。同时由于支持向量机不能预测太多的数据，BP网络法训练时间过大无法实时处理，而本文的LEMD方法则没有这方面的限制，在保证系统精度的条件下数据延拓的实时跟踪性强，更适合处理非平稳的突变信号。

表1 各种方法性能对比(仿真信号)Tab.1 Performance comparison of various methods(simulation signal)

采用LEMD方法对信号外延后进行分解得到的结果如图3所示，其中虚线为分解后的IMF分量，实线为真实信号中的各分量，延拓后的IFM分量和真实信号中的各分量之间的相关系数分别为0.991 1、0.995 4、0.940 3，两者非常接近。说明LEMD方法能够有效地对数据序列进行延拓。

图3 EMD分解图Fig.3 Diagrams of EMD decomposition

3 实验结果及其分析

选取某煤矿的一段现场监测微震信号数据来验证本文方法的有效性，信号采样频率为1 kHz，采用ESG的采集仪，用地震检波器拾取微震信号。现场监测是24 h监测，这里选用1个通道10 000个点的监测数据，截取中间1.5～5.596 s的一段，时域波形如图4(a)所示，图4(b)是其频谱图。

图4 原始信号和频谱Fig.4 Original signal and frequency spectrum

从图4(a)中可见，由于背景噪声的存在，从时域波形中无法获取明显地特征参数或定位信息，从图4(b)中可观测，微震信号特征频率不明显，整个频谱图中分布各种频率的信号。

对微震信号进行分析，就是从数据序列中找出或分离出能够反映定位或提取特征参数的信号。由于微震信号非平稳性强，且信号微弱，为了突出信号中的冲击成分的特征，应该先消除噪声，可以采用EMD方法对信号降噪，然后进行相应特征参数的提取。

图5 延拓效果图Fig.5 Extension effect diagram

为了对比本文方法的延拓效果，对截取的这一段数据向两端各延拓100个数据点，以右边延拓为例，在原数据从1.5～5.596 s这一段数据向右继续取100点数据，即5.596～5.696 s时间段数据，用于和延拓后的结果进行对比。效果对比如图5所示，其中虚线为原始信号，实线为基于 LEMD方法端点延拓后信号。LEMD方法在使用之前选取内部3 000点的数据对自适应滤波器进行训练，训练完成后进行数据延拓。

从图5可以看出，5.671～5.696 s时间段数据属于非平稳的突变信号，虽然延拓后的波形在幅值上和原数据有较大的误差，但延拓后的数据与原数据极值点跳变规律基本保持一致。

为验证LEMD方法抑制端点效应的优越性，以相同的分解终止条件，分别采用镜像延拓法、LEMD方法对图4的原始信号延拓并分解。图5是该信号利用镜像延拓抑制端点效应以后分解得到的各阶IMF分量;图6为采用LEMD抑制端点效应后分解得到的各阶IMF分量。

图6 采用镜像延拓后的EMD分解图Fig.6 EMD decomposition diagrams after mirror extension

图7 采用LEMD端点延拓后的分解图Fig.7 Decomposition diagrams after LEMD endpoint extension

通过图6和图7对比可以看出，利用镜像延拓的方法和LEM方法来抑制经验模态分解的端点效应后分解同样的信号，得到相同的IMF分量，都是13个IMF分量和一个残余分量。由图6可以看出，imf8～imf13分量和残余分量的右端点产生了较大的偏离，端点摆动现象很明显。而图7中经过LEMD处理端点效应以后，信号的端点效应得到了有效抑制，得到的IMF信号在端点处的偏离程度较小，尤其是信号的右端点更为明显。imf8～imf11分量、imf13分量和残余分量的右端点近似为0，分解得到的残差更为合理。综合分析比较，本文提出的LEMD算法分解精度更高，分解后的信号更加真实。

EMD分解后，对信号进行自适应低通滤波。利用EMD分解的结果重构信号，自适应的去掉信号高频成分。图8是采用镜像延拓后的EMD分解重构后信号及频谱图，图9为 LEMD方法分解重构后信号及频谱图。对比图8(b)和图9(b)可看出2种延拓方法分析结果的差异所在。尽管在图8(b)和图9(b)频谱图中均能明显地看出微震信号，但图9(b)噪声抑制的效果高于图8(b)。即基于本文方法对微震信号降噪的效果更好，对于定位信息的提取更加精确。

图8 镜像延拓后的EMD分解重构后信号及频谱Fig.8 Reconstruction signal and frequency spectrum after Mirror extension

图9 LEMD分解重构后信号及频谱Fig.9 Reconstruction signal and frequency spectrum after LEMD

4 结论

采用EMD方法对微震信号降噪，由于数据选取时的截断，造成EMD分解时引起端点效应问题，提出一种基于自适应滤波器延拓的LEMD端点延拓方法，该方法既考虑微震信号端点处的信号走向特点，又兼顾信号整体的发展趋势，仿真结果表明该方法适应性强，有效延拓距离长，延拓实时性好。将该方法应用于微震信号降噪、分析上，在抑制端点效应的同时，较好的提取出反映微震特性的信号，为微震监控与定位提供数据支持，具有较高的实用价值。

［1］李夕兵，张义平，刘志祥，等.爆破震动信号的小波分析与HHT变换［J］.爆炸与冲击，2005，25(6):528－535.LI Xi-bing，ZHANG Yi-ping，LIU Zhi-xiang，et al.Wavelet analysis and Hilbert-Huang Transform of blasting vibration signal［J］.Explosion and Shock Waves，2005，25(6):528－535.

［2］李夕兵，张义平，左宇军，等.岩石爆破振动信号的EMD滤波与消噪［J］.中南大学学报(自然科学版)，2006，37(1):150－154.LI Xi-Bing，ZHANG Yi-ping，ZUO Yu-Jun，et al.Filtering and denoising of rock blasting vibration signal with EMD［J］.Journal of Central South University (Science and Technology)，2006，37(1):150 －154.

［3］胡劲松，杨世锡.EMD方法基于AR模型预测的数据延拓与应用［J］.振动、测试与诊断，2007，27(2):116－119.HU Jing-song，YANG Shi-xi.AR model prediction-based EMD method and its application to data extension［J］.Journal of Vibration，Measurement ＆ Diagnosis，2007，27(2):116－119.

［4］邓拥军，王伟，钱成春，等.EMD方法及Hilbert变换中边界问题的处理［J］.科学通报，2001，46(3):257－263.DENG Yong-jun，WANG Wei，QIAN Cheng-chun，et al.Boundary processing technique in EMD method and Hilbert transfor m［J］.Chinese Science Bulletin，2001，46(3):257－263.

［5］张郁山，梁建文，胡聿贤.应用自回归模型处理EMD方法中的边界问题［J］.自然科学进展，2003，13(10):1054－1059.ZHANG Yu-shan，LIANG Jian-wen，HU Yu-xian.The processing of end effects in EMDmethod by autoregressive model［J］.Progresses in Nature Science，2003，13(10):1054－1059.

［6］程军圣，于德介，杨宇.基于支持矢量回归机的Hilbert-Huang变换端点效应问题的处理方法［J］.机械工程学报，2006，42(4):23 －31.CHENG Jun-sheng，YU De-jie，YANG Yu.Process method for end-effects of Hilbert-Huang transform based on support vector regression machines［J］. Chinese Journal of Mechanical Engineering，2006，42(4):23 －31.

［7］蔡艳平，李艾华，张玮，等.HHT端点效应的最大Lyapunov指数边界延拓方法［J］.仪器仪表学报，2011，32(6):1330－1336.CAI Yan-ping，LI Ai-hua，ZHANG Wei，et al.HHT end effect processing method based on maximum Lyapunov index boundary extension model［J］.Chinese Journal of Scientific Instrument，2011，32(6):1330 －1336.

［8］郭明威，倪世宏，朱家海，等.振动信号中 HHT/EMD端点延拓方法研究［J］.振动与冲击，2012，31(8):62－66.GUO Ming-wei，NI Shi-hong，ZHU Jia-hai，et al.HHT/EMD end extension method in vibration signal analysis［J］.Journal of Vibration and Shock，2012，31(8):62 －66.

［9］高强，段晨东，赵艳青，等.基于最大相关波形延拓的经验模式分解端点效应抑制方法［J］.振动与冲击，2013，32(2):62－66.GAO Qiang，DUAN Chen-dong，ZHAO Yan-qing，et al.A maximal correlation waveform extension method for end effects reduction of empirical mode decomposition［J］.Journal of Vibration and Shock，2013，32(2):62 －66.

［10］ Huang N E，Shen Z，Steven R L，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonsta-tionary time series analysis［J］.Proceedings of the.Royal Society A Mathematical，Physical and Engineering Sciences，1998，454:899 －995.

［11］ Huang N E，Wu M L C，Long SR，et al.A confidence limit for the empirical mode decomposition and Hilbert spectral analysis［J］.Proceedings of the Royal Society.Series A:Mathematical，Physical and Engineering Sciences，2003，459:2317－ 2345.

［12］谭善文，秦树人，汤宝平.Hilbert-Huang变换的滤波特性及其应用［J］.重庆大学学报，2004，27(2):9－12.TAN Shan-wen，QIN Shu-ren，TANG Bao-ping.The filtering character of Hilbert-Huang transform and its application［J］.Journal of Chongqing University，2004，27(2):9 －12.

［13］高晋占.微弱信号检测(第二版)［M］.北京:清华大学出版社.

［14］苏玉香，刘志刚，李科亮，等.一种改善EMD端点效应的新方法及其在谐波分析中的应用［J］.电工电能新技术，2008，27(2):33－37.SU Yu-xiang，LIU Zhi-gang，LI Ke-liang，et al.A new method for end effect of EMD and its application to harmonic analysis［J］.Advanced Technology of Electrical Engineering and Energy，2008，27(2):33－37.