凌同华，张 胜，易志强，李品钰

(长沙理工大学 土木与建筑学院，长沙 410004)

声发射技术的研究对象已由金属为主扩展到岩石、混凝土、复合材料等固体材料。作为一种无损检测手段，声发射技术已广泛地应用于地下采矿、石油化工、材料试验、民用工程、航天和航空等众多领域。

对声发射信号的分析处理是声发射技术能够得到更好地发展和推广的重要前提，岩石声发射信号具有瞬态性和多样性的特点，属于典型的非平稳信号。用傅里叶频谱、小波(包)分析处理声发射信号已得到广泛应用。刘希灵［1］编制快速傅里叶变换频谱分析程序对冲击荷载作用下岩石声发射信号进行频谱分析，得出冲击荷载作用下岩石声发射信号采集的时间短，信号的能量、信号强度、绝对能量在某一个波击数处有很大的值以及低频较显著、很集中，且低频和其他谐波成分的振幅都很大。王余刚［2］提取了混凝土材料不同破坏阶段的全波形声发射信号并分析了其频谱特性。康玉梅［3］基于小波变换，研究了声发射信号时频能量特征提取的原理和方法，利用时频能量分析技术来实现声发射信号的时延估计。赵奎［4］采用小波包频带分解方法，对岩石声发射Kaiser点信号的能量分布特征进行了研究，分析了砂岩声发射信号不同频带的能量分布规律。凌同华［5］利用小波包分析技术对冲击荷载作用下岩石声发射信号的能量分布特征进行了研究，重点讨论了冲击荷载作用下不同岩石对声发射信号频带能量分布的影响。聂鹏［6］基于声发射信号特点和小波包分解理论对不平稳信号特征提取的优势，提出了一种利用声发射信号的能量变化来监测刀具磨损状态的方法。赵元喜［7］根据滚动轴承声发射信号在各频段的能量分布与轴承的故障类型相关性，利用谐波小波包将不同故障滚动轴承的声发射信号分解到不同频段，进而将各频段的能量组成特征向量输入BP神经网络，通过神经网络判别滚动轴承的故障类型。

HHT(Hilbert-Huang Transform)法是近年提出的适合处理稳态信号和非平稳信号的有效方法，它由EMD方法和Hilbert变换两部分组成，其核心是EMD方法。相对小波(包)等分析方法而言，EMD(Empirical Mode Decomposition)方法是一种更具适应性的时频局部化分析方法，它没有固定的基函数，是自适应的［8－10］。本文针对冲击荷载作用下岩石声发射信号的特点，用EMD方法进行分析，探讨了冲击荷载作用下岩石声发射信号的能量分布特征。

1 EMD方法的原理

EMD方法依据信号本身的局部特征信息进行自适应地分解，得到一系列具有不同特征时间尺度的IMF(Intrinsic Mode Function)分量。IMF分量必须满足下面两个条件:① 在整个时间序列内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个;②在任意时刻，由局部极大值点、局部极小值点分别形成上、下包络线的平均值为零。

采用EMD方法通过以下步骤对岩石声发射信号x(t)进行分解［11－12］:

(1)确定信号所有的局部极值点，然后用三次样条曲线将所有的局部极大值点连接起来形成上包络线;

(2)用三次样条曲线将所有的局部极小值点连接起来形成下包络线，上、下包络线应该包络所有的数据点;

(3)上、下包络线的平均值记为m1(t)，求出:

在理想状况下，如果h1(t)满足IMF分量的两个基本条件，那么h1(t)就是x(t)的第1个IMF分量;

(4)如果h1(t)不满足IMF分量的两个基本条件，把h1(t)作为原始数据，重复步骤(1)～步骤(3)，得到上、下包络线的平均值m11(t)，再判断h11(t)=h1(t)－m11(t)是否满足IMF分量的两个基本条件，若不满足，则重复循环k次，得到 h1k(t)=h1(k－1)(t)－m1k(t)，使得h1k(t)满足IMF分量的两个基本条件。记c1(t)=h1k(t)，则c1(t)为信号x(t)的第1个IMF分量;

(5)从原始信号x(t)中分解出第1个IMF分量c1(t)后，从x(t)中减去c1(t)，得到剩余值序列r1(t):

将r1(t)作为原始数据重复步骤(1)～步骤(4)，得到x(t)的第2个IMF分量c2(t)，以此类推可得第3、第4个直至第n个IMF分量，即为c3(t)，c4(t)，…，cn(t)，最后剩下原始信号的余项rn(t)。

这样，原始信号x(t)可分解为若干IMF分量和一个余项rn(t)的和，即:

步骤(5)中的停止条件被称为分解过程的停止准则，它可以是如下两种条件之一:① 当最后一个IMF分量cn(t)或剩余分量rn(t)比预期值小时便停止;②当剩余分量rn(t)变成单调函数，从中不能再分解出IMF分量为止。

2 岩石声发射信号能量的表征

岩石声发射信号经EMD分解后的各IMF分量分别代表了一组特征尺度下的平稳信号，各平稳信号能量的变化就表征了岩石声发射信号的特征情况。

基于EMD方法的岩石声发射信号能量特征提取步骤如下［13－15］:

(1)对岩石声发射信号进行EMD分解，求得各IMF分量ci(t);

(2)求各IMF分量的总能量Ei:

(3)以能量为元素构造一个特征向量T:

当能量较大时，Ei(i=1，2，…，n)通常是一个较大的数值，在数据分析上会带来一些不方便的地方。由此，可以对特征向量T进行改进，即对向量进行归一化处理。令:

向量T'即为归一化后的向量。

3 冲击荷载作用下岩石声发射信号的采集

本次试验采用PCI－2声发射仪，其中有中心响应频率为250(500)kHz的PICO型谐振式窄频带传感器、2/4/6型前置放大器及AE–win声发射采集软件。试验选用长径比为1:2的三种不同的岩石试件(花岗岩、石灰岩、矽卡岩)，在杆径为50 mm的SHPB(Split Hopkinson Pressure Bar)试验系统上分别对它们进行冲击荷载作用下声发射试验。表1为试件的物理力学实验指标，其中花岗岩有轻微破碎，节理不发育，抗压强度高，表面硬度大，化学稳定性好，耐久性强;石灰岩有轻微破碎，节理裂隙不发育，具有良好的加工性、磨光性和很好的胶结性能;矽卡岩呈深灰色，节理不发育，块状构造，相对密度较大，且三种岩石试件完整性较好。每个试件在加工过程中进行标记，如“y1”对应花岗岩试件，“y2”对应石灰岩试件，“y3”对应矽卡岩试件。

表1 试件的物理力学实验指标Tab.1 Specimen of physical and mechanical test indexes

将声发射探头(传感器)的接触面涂抹少量的耦合剂(如黄油等)，其目的是获得高质量的信噪比及减少背景干扰信号，然后用胶带绑在试件表面居中位置，最后将岩石试件夹在SHPB试验装置入射杆和透射杆之间。试验中花岗岩、矽卡岩的采样频率为4×107Hz，而石灰岩采样频率为5×106Hz，为便于分析，被分析信号的采样频率必须一致，对石灰岩信号进行重新采样(重新采样并不改变信号的本质特征)，重新采样的频率为4×107Hz。试验时采集了多组信号，现取每次试验声发射起始点的信号进行分析。图1为冲击荷载作用下岩石声发射信号的幅值时程曲线。

图1 冲击荷载作用下岩石声发射信号的幅值时程曲线Fig.1 The amplitude versus time curves of rock acoustic emission signals under impact loading

4 声发射信号的EMD分解及能量分布特征分析

4.1 预处理

实测冲击荷载作用下岩石声发射信号除具有反映有关岩石本身的特征信息外，还包含大量的背景噪声。因此，在提取特征向量前有必要对信号进行预处理，以突出特征信息，提高信号的信噪比。

4.2 EMD分解层

EMD方法能对非平稳、非线性信号进行平稳化、线性化处理，并在分解的过程中保留原始信号的固有特性。因此，在Matlab语言平台上采用EMD方法对采集的声发射信号进行分解得到n个IMF分量c1(t)，c2(t)，…，cn(t)，每一个IMF分量都包含了不同的特征尺度信息，这样通过EMD分解，原始信号的特征完全可以由这n个IMF分量c1(t)，c2(t)，…，cn(t)来表征。考虑到原始信号的特征信号主要集中在高频带，因此可以只选取前8个IMF分量作为进一步的研究对象。岩石声发射信号y1，y2，y3的EMD分解结果见图2～图4。

图2 岩石声发射信号y1的EMD分解结果Fig.2 The EMD results of rock acoustic emission signal y1

图3 岩石声发射信号y2的EMD分解结果Fig.3 The EMD results of rock acoustic emission signal y2

图4 岩石声发射信号y3的EMD分解结果Fig.4 The EMD results of rock acoustic emission signal y3

从图中可以看出，第1个IMF分量是原始信号中分解出的时间尺度最短、频率最高的分量，代表信号中的高频成分，且其振幅大，表明其所占能量大，然后依次分解出余下各IMF分量。随着分解的进行，所得各IMF分量时间尺度越来越长、频率越来越低，直到分解出时间尺度最长、频率最低的最后一个IMF分量。

EMD方法把冲击荷载作用下岩石声发射信号分解成8个IMF分量，不同的IMF分量包含不同的特征时间尺度，可以使信号的特征信息在不同的分辨率下显示出来，表明EMD分解中分辨率是自适应的，相对小波的多分辨率分析，EMD方法在分解信号时更加简单。EMD分解中没有固定的基函数，每个IMF分量的提取都是由信号本性所决定的，分解出的各IMF分量与原始信号相似。在数字信号处理中常用相关系数来反映两个信号的相似程度［16］，各IMF分量与原始信号相关系数计算结果见表2。

表2 IMF分量与原始信号的相关系数Tab.2 Correlation coefficients between the IMF components and the original signals

4.3 岩石声发射信号的能量分布特征

冲击荷载作用下岩石声发射信号经EMD方法分解后得到8个IMF分量，按式(4)～式(7)分别提取每个IMF分量的特征向量T'。计算结果见表3、图5。

表3 岩石声发射信号的能量分布Tab.3 Energy distribution of rock acoustic emission signals

(1)从表3和图5可以看出，3种岩石声发射信号经EMD分解出的前4个IMF分量的能量占该信号总能量的比例分别为 93.96%，85.43%，85.50%，表明冲击荷载作用下岩石声发射信号的能量在频域上分布比较广泛，但绝大部分能量集中在前4个IMF分量内，且分布不均匀。

(2)为了显示岩石声发射信号各IMF分量相关系数与能量分布的关系，以各IMF分量为横坐标，将表2和表3中的数值绘制成曲线，如图6所示。从图6中可以看出各IMF分量相关系数与能量分布存在很大相关性，两者的相关系数分别为0.94，0.84，0.96，表明各 IMF分量相关系数与能量分布的相关程度高，反映了EMD方法是依据信号本身的固有特性进行自适应地分解。

图5 岩石声发射信号的能量分布Fig.5 The energy distribution of rock acoustic emission signals

图6 各IMF分量的相关系数与能量分布Fig.6 Correlation coefficients and energy distribution of IMF components

图7 岩石声发射信号的功率谱Fig.7 The power spectrum of rock acoustic emission signals

图8 岩石声发射信号y1的IMF分量(c1～c3)的功率谱Fig.8 The power spectrum of IMF components(c1～ c3)of rock acoustic emission signal y1

(3)冲击荷载作用下岩石声发射信号的频谱丰富，且绝大部分分布在500 kHz以下。从图7中可知3种岩石声发射信号在0～500 kHz的能量占该信号总能量的比例分别为 99.56%，97.42%，85.32%，并可分成多个子震频带，如图8中分量c1，c2及c3各自的主频率所示，这一点与图7中声发射信号y1的频谱基本一致。与此同时，图7中花岗岩声发射信号的优势频率较石灰岩、矽卡岩更集中，这一特征同参考文献［5］是一致的，但EMD方法没有固定的基函数，具有自适应性，使信号分析更加灵活多变;而小波包分析中不同的小波基具有不同的特性，对信号的分析效果也不相同，对同一个信号采用不同的小波基得到的结果基本上也没有可比性。

(4)从表1、图6和图7还可以看出，随着岩石单轴压缩强度增加，冲击荷载作用下岩石声发射信号分解出的各IMF分量相关系数与能量分布的相关程度越高;随着岩石的密度、纵波波速、弹性模量的降低，冲击荷载作用下岩石声发射信号的优势频率越来越集中，且其优势频率有往低频发展的趋势。

5 结论

(1)EMD方法把冲击荷载作用下岩石声发射信号分解为有限个不同特征时间尺度的IMF分量，其信号能量主要分布在前4个IMF分量内，且分布不均匀;

(2)EMD方法分解出的IMF分量大都有清晰的物理意义，且各IMF分量的频谱与原始信号的频谱基本一致，表明EMD方法分解出的各IMF分量集中了原始信号中最显著的特征信息;

(3)随着岩石的密度、纵波波速、弹性模量的降低，冲击荷载作用下岩石声发射信号的优势频率越来越集中，且其优势频率有往低频发展的趋势。相比小波法分析［5］，EMD分析不需预设基函数，具有自适应和高效的特点。

［1］刘希灵.基于激光三维探测的空区稳定性分析及安全预警的研究［D］.长沙:中南大学，2008.

［2］王余刚，骆 英，柳祖亭.全波形声发射技术用于混凝土材料损伤监测研究［J］.岩石力学与工程学报，2005，24(5):803－807.

［3］康玉梅，朱万成，白 泉，等.基于小波变换时频能量分析技术的岩石声发射信号时延估计［J］.岩石力学与工程学报，2010，29(5):1010 －1016.

［4］赵 奎，王更峰，王晓军，等.岩石声发射Kaiser点信号频带能量分布和分形特征研究［J］.岩土力学，2008，29(11):3082－3088.

［5］凌同华，廖艳程，张 胜.冲击荷载下岩石声发射信号能量特征的小波包分析［J］.振动与冲击，2010，29(10):127－130.

［6］聂 鹏，王东磊，徐 涛，等.频带能量特征法在声发射刀具磨损监测系统中的应用［J］.工具技术，2009，43(2):24－27.

［7］赵元喜，胥永刚，高立新，等.基于谐波小波包和BP神经网络的滚动轴承声发射故障模式识别技术［J］.振动与冲击，2010，29(10):162 －165.

［8］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proc.R.Soc.Lond.A，1998，454:903 －995.

［9］Huang N E，Shen Z，Long S R.A new view of nonlinear water waves:the hilbert spectrum［J］.Annual Review of Fluid Mechanics，1999，31:417 －457.

［10］ Ma S，Zhang R.Empirical mode decomposition of the 1994 northridge earthquake and its interpretation for seismic source mechanism［A］.The 10th International Conference on Soil Dynamics and Earthquake Engineering.USA，Philadelphia:［s.n.］，2001:7 －10.

［11］于德介，程军圣，杨 宇.机械故障诊断的Hilbert－Huang变换方法［M］.北京:科学出版社，2006.

［12］李夕兵，凌同华，张义平.爆破震动信号分析理论与技术［M］.北京:科学出版社，2009.

［13］杨 宇，于德介，程军圣.基于EMD与神经网络的滚动轴承故障诊断方法［J］.振动与冲击，2005，24(1):85－88.

［14］胡昌华，李国华，周 涛.基于MATLAB7.X的系统分析与设计－小波分析(第三版)［M］.西安:西安电子科技大学出版社，2008.

［15］凌同华，李夕兵.单段爆破振动信号频带能量分布特征的小波包分析［J］.振动与冲击，2007，26(5):41－43.

［16］楼顺天，刘小东，李博涵.基于MATLAB7.X的系统分析与设计－信号处理(第二版)［M］.西安:西安电子科技大学出版社，2005.