基于改进EMD的地下金属矿山围岩声发射信号去噪声处理

罗周全¹，左红艳^{1, 2}，吴超¹，王益伟¹

(1. 中南大学 资源与安全工程学院，湖南 长沙，410083；

2. 湖南涉外经济学院 商学院，湖南 长沙，410205)

摘要：为提高地下金属矿山围岩声发射信号检测精度，提出一种可用于地下金属矿山围岩声发射信号去噪声处理的改进EMD方法。该方法首先将地下金属矿山围岩声发射信号进行经验模态分解，然后，通过确定高通滤波器、低通滤波器和带通滤波器所组成的滤波器组的滤波等级而选择相应滤波器进行滤波，最后，将剩余IMF分量进行重构得到去噪声处理后的地下金属矿山围岩声发射信号。研究结果表明：该改进EMD方法的去噪声误差较小，重构后的地下金属矿山围岩声发射能率信号能较好地反映地下金属矿山围岩声发射能率信号的真实趋势。

关键词：地下金属矿山；声发射信号；经验模态分解方法；本征模函数；去噪声处理

中图分类号：TF811           文献标志码：A         文章编号：1672-7207(2013)11-4694-08

Denoising processing on acoustic emission signals from underground metal mine rock by using improved EMD method

LUO Zhouquan¹, ZUO Hongyan^{1, 2}, WU Chao¹, WANG Yiwei¹

(1. School of Resource and Safety Engineering, Central South University, Changsha 410083, China;

2. Commercial College, Hunan International Economics University, Changsha 410205, China)

Abstract: In order to improve measurement accuracy of acoustic emission signals from underground metal mine rock, an improved EMD method was put forward, which can be used to remove the noises from acoustic emission signals from underground metal mine rock. Noise removing process on the acoustic emission signals from underground metal mine rock by use of the improved EMD method could be expressed as follows: Firstly, acoustic emission signals from underground metal mine rock were decomposed into a series of IMFs by using empirical mode decomposition. Secondly, appropriate filter for filtering IMFs was selected when filter grade of filter group including high-pass filter, low-pass filter and band-pass filter was defined. Finally, remainder IMFs were used to reconstruct acoustic emission signals which noises from underground metal mine rock were removed. The results show that noise removing error based on the improved EMD method is small, and reconstructed acoustic emission energy rate signals can present true trend of acoustic emission energy rate signals from underground metal mine rock.

Key words: underground metal mine; acoustic emission signals; empirical mode decomposition; intrinsic mode function; denoising processing

矿山开采对地下围岩造成采动损害，改变了地下岩体的平衡状态^[1-2]，采动影响导致地下工程发生冲击地压是常见的灾害事故。由于冲击地压的演化和形成是一个非线性过程，且非线性是导致冲击地压的根本原因^[3-6]，因此，不能用传统的确定论观点来研究矿山冒顶这一非线性动力学过程^[7-10]。采用非线性动力学的理论和方法来研究冲击地压灾害的孕育、演化和发生机制，并建立相应的预测方法，是从本质上认识和解决冲击地压灾害的正确途径^[1]。在地下金属矿山围岩声发射信号实际测试过程中，采集到的围岩声发射信号将不可避免地受到白噪声、局部高频化噪声等干扰，导致获取的围岩声发射信号出现不同程度的失真^[11]。傅立叶滤波算法^[12]可滤除高斯白噪声，但同时也会滤除围岩声发射信号中十分重要的高频信息。样条拟合方法^[13]能较好地抑制高斯白噪声，但同样会将无关信息引入地下金属矿山围岩声发射信号时间序列等信号。小波变换极大模算法的主要缺点是准确度不高^[14-15]，而经验模态分解(empirical mode decomposition, EMD) 方法是由Huang等^[16-18]发展的一种新的数据分析方法。该方法能有效地去除高斯白噪声对采集信号的干扰，而且不删除采集信号中的有用信息，不引入无关信息，消噪效果好，但经验模态分解方法中需要依靠主观经验选择合适的连续本征模函数(intrinsic mode function，IMF)进行组合重构信号，其滤波信噪比必然会受到主观经验的制约。为此，本文作者提出一种地下金属矿山围岩声发射信号EMD新型滤波器，可通过确定高通滤波器、低通滤波器和带通滤波器所组成的滤波器组的滤波等级k而选择相应滤波器进行滤波重构信号的改进EMD方法，为地下金属矿山围岩声发射信号提取提供技术保障。

1  基于改进EMD的地下金属矿山围岩声发射信号去噪声方法

1.1  地下金属矿山围岩声发射信号的经验模态分解条件

地下金属矿山围岩声发射信号IMF分量要满足以下条件：(1) 整个数据序列的极大、极小值数目与过零点数目相等或最多相差1个；(2) 地下金属矿山围岩声发射信号数据序列的任意一点由极大值所确定的包络与由极小值所确定的包络均值始终为0。

这2个条件使分解得到的所有地下金属矿山围岩声发射信号IMF分量是窄带信号，而且这种分解应满足以下假设：(1) 地下金属矿山围岩声发射信号至少有1个极大值和1个极小值；(2) 特征时间尺度由极值间的时间间隔确定；(3) 若数据中没有极值点，而只有拐点，则可通过一阶或多阶微分得到极值点，再通过分解、积分的方法获得IMF分量。

1.2  地下金属矿山围岩声发射信号EMD端点问题处理

基于三次样条插值的特点，当原始的地下金属矿山围岩声发射信号的两端点不是极值点时，若能够根据端点以内极值点序列的规律得到地下金属矿山围岩声发射信号时间序列在端点处的近似值，则可以防止对极值点进行样条插值得到的包络线出现极大的摆动。取出地下金属矿山围岩声发射信号极值点序列最左端1/3的极值点，根据地下金属矿山围岩声发射信号数据的间距均值和左端点的幅值，确定左端点需增加的极值点位置和幅值。同理，可确定右端点需增加的极值点位置和幅值，当所构成的新的极大值和极小值数据集的最大间距小于原始的地下金属矿山围岩声发射信号时，以近似的左端点处增加的极值点为起始点，向左进行数据对称延拓；以近似的右端点处增加的极值点为起始点，向右进行数据对称延拓。将得到的新的数据作为一个整体进行经验模态分解。需要指出的是所取的结果只是中间部分。虽然该方法只是求出左、右端的近似值，对近似极值点进行对称延拓，但延拓的目的不是为了给出准确的原序列以外的数据，而是提供一种条件，使得包络完全由端点以内的数据确定。

1.3  地下金属矿山围岩声发射信号EMD过程

对于包含N个地下金属矿山围岩声发射信号而言，其EMD过程具体步骤如下。

(1) 找出地下金属矿山围岩声发射信号所有的局部极值点，用三次样条线将所有的局部极大值点和局部极小值点连接起来形成上、下包络线，上、下包络线应包含所有的数据点。2条包络线的平均值记为m₁(t)，数据X(t)与m₁(t)之差为h₁(t)：

h₁(t)=X(t)-m₁(t)               (1)

式中：t=1, 2, …, N。若h₁(t)满足IMF的定义，则h₁(t)是1个IMF，h₁(t)就是X(t)的第1个分量。

(2) 若h₁(t)不满足IMF的定义，则将h₁(t)作为原始数据，重复以上步骤，得到：

h₁₁(t) = h₁(t)-m₁₁(t)            (2)

式中：m₁₁(t)为h₁(t)的上、下包络线的平均值。判断h₁₁(t)是否满足IMF的定义，若不满足，则重新循环k次，得到h_1k(t) = h_1(k-1)(t)-m_1k(t)使m_1k(t)满足IMF的定义，记c₁(t)= h_1k(t)。

(3) 将c₁从数据X(t)中分离出来得到：

r₁(t) =X(t)-c₁(t)            (3)

再把r₁(t)作为新的原始数据，重复以上步骤，就得到第2个满足IMF的分量c₂(t)。重复循环，得到：

r_j-1(t)-c_j(t)=r_j(t)；j=2, 3, 4, …, n      (4)

为判断所处理地下金属矿山围岩声发射信号不再含有IMF分量，一般采取地下金属矿山围岩声发射信号IMF分量结束循环条件。地下金属矿山围岩声发射信号IMF分量满足条件(2)过于苛刻，会删除掉具有物理意义的幅度波动。因此，为保证每一个地下金属矿山围岩声发射信号IMF分量具有幅度和频率调制的物理意义，将条件(2)转化为比较容易实现的数量标准。该标准由下式给出：

        (5)

式中：m_1k (t)为IMF分量提取模块中本次循环过程中求得的平均包络；m_1(k-1)(t)为上次循环过程中求得的平均包络；0，…，T为平均包络线所包含的所有时刻。理想的应该在0.2~0.3之间。满足上述2个条件的IMF分量，既是进一步进行Hilbert变换的基础，又保证了每个分量蕴含必要的物理意义。

直到r_n(t)成为一个单调函数不能从中提取满足地下金属矿山围岩声发射信号IMF的分量时，循环结束。

这样，原始数据可以由地下金属矿山围岩声发射信号IMF分量和最后残量之和表示为：

                 (6)

所以，可以把地下金属矿山围岩声发射信号X(t)分解为n个地下金属矿山围岩声发射信号IMF和1个地下金属矿山围岩声发射信号残量r_n(t)之和。其中：分量c₁(t)，c₂(t)，…，c_n(t)分别包含了信号从高到低不同频率段的成分，而且不是等带宽的，EMD方法是一个自适应的信号分解方法，这是它优于小波包分析法之处。

(4) 估计IMF系数的小波阈值方法。对于噪声和声发射信号在IMF成分混叠的情况，直接利用尺度滤波方法不能除去噪声。这时，可采用小波变换中的软硬阈值折中法对每一个IMF成分作阈值处理。则阈值处理后的IMF系数可变表示为：

硬阈值：

         (7)

式中：为阈值，；为噪声的方差；N为噪声和声发射信号在IMF成分混叠信号的长度。

软阈值：

    (8)

1.4  地下金属矿山围岩声发射信号EMD新型滤波器设计

从地下金属矿山围岩声发射信号IMF筛选的过程看，EMD分解过程总是不断提取地下金属矿山围岩声发射信号中的高频部分，即c_j(t)总是包含比c_j-1(t)更低的频率部分。可见：通过选择不同连续的IMF进行线性组合可以得到不同频段的信号，地下金属矿山围岩声发射信号EMD滤波方法可以等效为确定高通滤波器、低通滤波器和带通滤波器所组成的滤波器组的滤波等级k，并选择相应滤波器进行滤波的方法。因此，地下金属矿山围岩声发射信号高通滤波器表示为

              (9)

其中：k₁为高通滤波器参数，k₁=n/(Q+1)；Q为地下金属矿山围岩声发射信号品质因子，可根据地下金属矿山围岩声发射信号特征在1~3范围内取值。

地下金属矿山围岩声发射信号低通滤波器表示为

         (10)

式中：k₂为低通滤波器参数，k₂=n-n/(Q+1)。地下金属矿山围岩声发射信号带通滤波器表示为

             (11)

设地下金属矿山围岩声发射信号由反映围岩动态特征的期望信号s(t)和白噪声e(t)组成，则地下金属矿山围岩声发射信号X(t)可表示为

X(t)= s(t)+e(t)             (12)

地下金属矿山围岩声发射信号X(t)的变化趋势主要取决于期望信号s(t)的变化，而噪声e(t)表现为在期望信号s(t)上的波动，s(t)的变化相对e(t)是一低频信号。在EMD分解过程中，期望信号s(t)的能量多数集中在k级低频部分的IMF中，而噪声e(t)能量则多数集中在k-1级表示相对高频部分的IMF中，所以，地下金属矿山围岩声发射信号重构信号X_k(t)的信噪比f[X_k(t)]应高于地下金属矿山围岩声发射信号重构信号的信噪比f[X_k-1(t)]，即

f [X_k(t)]＞f [X_k-1(t)]            (13)

为此，采用地下金属矿山围岩声发射信号IMF分量的相邻2个重构信号的均方误差σ(k)来检测地下金属矿山围岩声发射信号重构信号信噪比的提高程度，即

   (14)

当地下金属矿山围岩声发射信号相邻2个重构信号的均方误差σ(k)达到最大值时，地下金属矿山围岩声发射信号重构信号X_k(t)具有最大的信噪比，因此，重构信号X_k(t)的等级k可用下式确定：

       (15)

式(14)和(15)表明地下金属矿山围岩声发射信号期望信号s(t)的能量主要集中在一定时间尺度内波动幅度最大的IMF中，因此，式(13)~(15)提供了确定地下金属矿山围岩声发射信号重构信号X_k(t)的IMF等级k的方法。η通过确定k即可按照式(9)~(11)重构地下金属矿山围岩声发射信号的输出信号，从而实现对地下金属矿山围岩声发射信号的去噪声处理。

1.5  基于改进EMD的去噪声方法仿真实例

为验证基于改进EMD的地下金属矿山围岩声发射信号去噪声方法的有效性，采用Matlab2008生成如图1所示的近似正弦函数动态信号y，同时对图1所示的近似正弦函数动态信号添加随机信号y而生成如图2所示的含噪声近似正弦函数动态信号。

图1  近似正弦函数y动态信号

Fig.1  Dynamic signals of approximate sine function

图2  含噪声近似正弦函数y动态信号

Fig.2  Dynamic signals with noises of approximate sine function

分别采用EMD方法、小波包分解方法和本文提出的改进EMD方法对图2所示的含噪声近似正弦函数动态信号进行去噪声处理后得到的动态信号与图1所示信号的相对误差如图3所示。

由图3可见：小波包分解重构信号的最大相对误差为7.88%，最小相对误差为1.90%，平均相对误差为4.83%；EMD重构信号最大相对误差为6.90%，最小相对误差为3.10%，平均相对误差为4.45%。而本文提出的改进EMD重构信号最大相对误差为4.11%，最小相对误差为1.34%，平均相对误差为2.95%。可见：本文提出的改进EMD方法的去噪声误差明显小于小波包分解重构方法和EMD方法的去噪声误差。

图3  3种去噪处理方法的相对误差η比较

Fig.3  Comparison to relative errors based on three denoising processing methods

2  地下金属矿山围岩声发射信号去噪声处理实例

2.1  地下金属矿山围岩声发射能率E信号经验模态分解和Hilbert谱

以某南方地下铅锌矿围岩声发射信号为例对基于改进EMD的地下金属矿山围岩声发射信号去噪声方法的应用效果进行实例验证。在某南方地下铅锌矿采场测取500个含噪声的地下金属矿山围岩声发射能率E信号，如图4所示。

图5所示为地下金属矿山围岩声发射能率E信号通过经验模态分解得到如的IMF1~IMF8分量。从图5可见：IMF1~IMF8依次表示本征模函数由低阶向高阶发展变化，即依次表示高频段成分占主要成分的地下金属矿山围岩声发射能率E信号向低频段成分占主要部分的地下金属矿山围岩声发射能率E信号自适应分解变换，逐次把地下金属矿山围岩声发射能率E信号的频率分量从高到低分解出来。由于经验模态分解方法本身的特点，地下金属矿山围岩声发射能率E信号的经验模态分解基底函数是自适应的，因此，不会像传统的信号分析方法那样受到先验基底函数的影响，所得到的地下金属矿山围岩声发射能率E信号的IMF1~IMF8分量是信号直接和真实的反映，其IMF1~IMF8分量的Hilbert谱如图6所示。

图4  含噪声的地下金属矿山围岩声发射信号

Fig.4  Acoustic emission signals with noises from underground metal mine rock

2.2  地下金属矿山围岩声发射能率E信号重构信号的确定

由图6所示的地下金属矿山围岩声发射能率E信号的IMF1~IMF 8分量的Hilbert谱可知：地下金属矿山围岩声发射能率E信号的能量主要集中在IMF7和IMF8上，其频率集中在0.05~4.50 Hz，故对地下金属矿山围岩声发射能率E信号的IMF分量进行重构时应该舍去IMF1，IMF2，IMF3 ，IMF4，IMF5与IMF6，只对IMF7和IMF8进行重构。

图5  地下金属矿山围岩声发射信号IMF分量

Fig.5  IMF components of acoustic emission signals from underground metal mine rock

图6  地下金属矿山围岩声发射能率E信号的Hilbert谱

Fig.6  Hilbert spectrum of acoustic emission energy rate signals from underground metal mine rock

此外，采用式(14)和式(15)可得地下金属矿山围岩声发射信号IMF分量的相邻2个重构信号信噪比提高等级k的变化如图7所示。由图7可知：地下金属矿山围岩声发射能率E信号的相邻2个重构信号的均方误差σ(k)达到最大值时所对应的期望信号能量主要集中在第7级IMF及其后级IMF中。地下金属矿山围岩声发射信号重构信号X_k(t)的IMF等级k的取值为7。故对地下金属矿山围岩声发射能率E信号的IMF分量进行重构时也应该舍去IMF1，IMF2，IMF3，IMF4，IMF5与IMF6，而只对IMF7和IMF8进行重构。可见：采用Hilbert谱对地下金属矿山围岩声发射能率E信号的IMF分量进行重构和采用地下金属矿山围岩声发射信号重构信号X_k(t)的IMF等级k的方法进行重构是一致的，但更为快捷和有效。

按照式(10)重构的地下金属矿山围岩声发射能率E信号如图8所示。从图8可以看出：经改进EMD方法去噪处理后，重构后的地下金属矿山围岩声发射能率E信号波形变得光滑，也很好地保留了尖峰和突变部分，为较明显的阻尼振荡过程，能较好地反映地下金属矿山围岩声发射能率E信号的真实趋势。

图7  相邻2个重构信号信噪比的变化

Fig.7  Change of signal to noise ratio for adjacent two reconstructed signals

图8  地下金属矿山围岩声发射能率E信号IMF分量重构

Fig.8  IMF components reconstruction of acoustic emission energy rate signals from underground metal mine rock

3  结论

(1) 针对经验模态分解重构信号的滤波信噪比高低受主观经验制约问题。提出一种地下金属矿山围岩声发射信号EMD新型滤波器，可通过确定高通滤波器、低通滤波器和带通滤波器所组成的滤波器组的滤波等级k而选择相应滤波器进行滤波重构信号的改进EMD方法。改进EMD方法的去噪声误差较小。

(2) 利用改进EMD方法对地下金属矿山围岩声发射信号进行了去噪声处理，重构后的地下金属矿山围岩声发射能率E信号波形变得光滑，也很好地保留了尖峰和突变部分，为较为明显的阻尼振荡过程，能较好地反映地下金属矿山围岩声发射能率E信号的真实趋势。

参考文献：

[1] 左红艳. 地下金属矿山开采安全机理辨析及灾害智能预测研究[J]. 长沙: 中南大学资源与安全工程学院, 2012: 95-100.

ZUO Hongyan. Analysis on mining safety mechanism of underground metal mine and its disaster intelligent prediction[J]. Changsha: School of Resource and Safety Engineering. Central South University, 2012: 95-100.

[2] Palei S K, Das S K. Sensitivity analysis of support safety factor for predicting the effects of contributing parameters on roof falls in underground coal mines[J]. International Journal of Coal Geology, 2008, 75(4): 241-247.

[3] Palei S K, Das S K.Logistic regression model for prediction of roof fall risks in bord and pillar workings in coal mines: An approach[J]. Safety Science, 2009, 47(1): 88-96.

[4] LIU Hongtao, MA Nianjie, GONG Peinuan, et al. The technology research about hidden danger identification of tunnel roof fall[J]. Procedia Engineering, 2011, 26: 1220-1224.

[5] Phillipson S E. Texture, mineralogy, and rock strength in horizontal stress-related coal mine roof falls[J]. International Journal of Coal Geology, 2008, 75(3): 175-184.

[6] Maiti J, Khanzode V V. Development of a relative risk model for roof and side fall fatal accidents in underground coal mines in India[J]. Safety Science, 2009, 47(8): 1068-1076.

[7] Bobick T G, Jr McKenzie E A, Kau T Y. Evaluation of guardrail systems for preventing falls through roof and floor holes[J]. Journal of Safety Research, 2010, 41(3): 203-211.

[8] Bertoncini C A, Hinders Mark K. Fuzzy classification of roof fall predictors in microseismic monitoring[J]. Measurement, 2010, 43(10): 1690-1701.

[9] ZHENG Minggui, CAI Sijing. An intelligent system for calculating the scale of rational, enlarged production of an underground non-ferrous metal mine[J]. Journal of China University of Mining and Technology, 2008, 18(2): 214-219.

[10] Asad M W, Dimitrakopoulos R. Optimal production scale of open pit mining operations with uncertain metal supply and long-term stockpiles[J]. Resources Policy, 2012, 37(1): 81-89.

[11] 左红艳, 罗周全, 王益伟, 等. 基于模糊自适应变权重算法的采场冒顶函数链神经网络预报[J]. 中国有色金属学报, 2011, 21(4): 894-900.

ZUO Hongyan, LUO Zhouquan, WANG Yiwei, et al. Functional link neural network to predict roof caving based on method of fuzzy adaptive variable weight[J]. The Chinese Journal of Nonferrous Metals, 2011, 21(4): 894-900.

[12] Alsdorf D. Noise reduction in seismic data using Fourier correlation coefficient filtering[J]. Geophysics, 1997, 62(5): 1617-1627.

[13] 吴先良, 焦丹, 王良之, 等. 从目标瞬态响应中提取极点的样条拟合及有理逼近方法[J]. 中国科学技术大学学报, 1996, 26(4): 528-533.

WU Xianliang, JIAO Dan, WANG Liangzhi, et al. A method for singularity extraction from targets transient response with spline functions and rational approximation[J]. Journal of China University of Science and Technology, 1996, 26(4): 528-533.

[14] 焦丹, 徐善驾, 吴先良, 等. 采用频域紧支集正交小波基消除瞬态散射回波中的高斯白噪声干扰[J]. 电子学报, 1999, 27(6): 120-122.

JIAO Dan, XU Shanjia, WU Xiangliang, et al. Recovery of signal from transient scattered response contaminated by gaussian white noise based on orthogonal bases of compactly supported wavelets in frequency domain[J]. Acta Electronica Sinica, 1999, 27(6): 120-122.

[15] Jr O M, Domingues M O, Costa A M D, et al. Wavelet analysis applied to magnetograms: Singularity detections related to geomagnetic storms[J]. Journal of Atmospheric and Solar-Terrestrial Physics, 2005, 67(17/18): 1827-1836.

[16] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proc R soc Lond A, 1998, 454: 903-995.

[17] 鄂加强, 王春华, 龚金科, 等. 铜火法冶炼热动力学系统实测数据EMD处理[J]. 中国有色金属学报, 2008, 18(5): 946-951.

E Jiaqiang, WANG Chunhua, GONG Jinke, et al. Process on measurement data from copper pyrometallurgical heat dynamical system by using of EMD method[J]. The Chinese Journal of Nonferrous Metals, 2008, 18(5): 946-951.

[18] Echeverria J C. Application of empirical mode decomposition to heart rate variability analysis[J]. Medical & Biological Engineering & Computing, 2001(39): 471-479.

(编辑  陈灿华)

收稿日期：2013-02-27；修回日期：2013-04-13

基金项目：国家自然科学基金资助项目(51274250)；国家“十二五”科技支撑计划项目(2012BAK09B02-05)

通信作者：左红艳(1975-)，女，湖南湘潭人，博士，从事矿产安全研究；电话：15084931748；E-mail: zuohongyan18@126.com