爆破振动信号分析中模态混叠和虚假分量消除的改进方法

2019-09-17徐文龙冯丹丹

振动与冲击 2019年17期

李清，徐文龙，张迪，李娜，冯丹丹

(1.中国矿业大学(北京) 力学与建筑工程学院，北京 100083; 2.广州大学土木工程学院，广州 510006)

爆破振动信号是爆破设计参数的外在表现。通过分析信号的峰值振速、主频、能量在不同时间和不同频带下的分布特征能有效的控制爆破振动强度，减小爆破危害。以傅氏变换为基础的信号变换方法主要针对平稳信号分析，显然对于爆破振动信号分析存在局限性。Huang等[1]提出了HHT变换，并广泛应用于爆破振动等非平稳信号的分析。HHT变换分成两步，首先对信号进行EMD分解，得到一组频率由高到低的模态函数IMF，然后对IMF进行Hilbert变换，得到信号的时频图谱、能量谱等。但是EMD中存在一些需要改进的缺陷，如模态混叠和虚假分量。即任一IMF里含有明显不相同的时间尺度或同一范围的时间尺度出现在不同的IMF中。

针对模态混叠问题，Yang等[2]分析得到，在间断信号处包络线产生的过冲是产生模态混叠的重要原因。胡爱军等[3]采用非线性形态滤波器，通过布尔运算和相应的加减计算来分离信号内的脉冲分量，保留原始信号的局部特性，减少信号中的脉冲噪声来抑制模态混叠现象。高云超等[4]采用对混合分量一阶差分的方法对分量进行再分解，对分解得到的新的待验证IMF分量进行累计求和，消除模态混叠现象。周颖涛等[5]通过改进EMD筛分停止阀值，增加单次筛分，使包络平均值小于原始幅值抑制信号分解的的模态混叠现象。郑近德等[6]检测各IMF分量排列熵，剔除前几阶异常分量，将剩余分量再次进行EMD分解抑制模态混叠现象。Ryan等[7]提出了对原始信号添加合适的掩膜信号来抑制模态混叠。赵玲等[8]、凌立鹏等[9]分别对Ryan的方法进行了改进，认为由于EMD能量泄露的影响,仅通过第一个IMF分量的包络幅值和瞬时频率来定义s(t)的幅值和频率对模态混叠的抑制效果较差，改进的方法中将表征能量泄露分量的IMF的幅值与频率和IMF1相加获得新的s(t)。当信号内噪声成分较多时，掩膜信号的确定十分困难，不便于操作。秦品乐等[10]对信号进行先期小波阀值消噪，消除信号内的突变成分，降低包络线的波动。

关于虚假分量的判断，黄迪山[11]根据EMD分解的完备性，利用能量守恒准则，认为当EMD分解产生虚假分量时，各分量的均方值和要大于原始信号的均方值。将第一阶本征模态分量与虚假分量相加，消除虚假分量的同时更新第一阶模态分量，得到最终的各阶分量。胡爱军[12]基于互信息理论，将各IMF分量与原始信号进行互信息计算判断虚假分量。李纪永等[13]通过改进的SVD对原始信号降噪，对降噪信号进行EMD分解，利用频谱比值法判断虚假分量。此外，判断虚假分量方法还有K-L散度法[14]，相关系数法[15-16]等。但对于虚假分量如何处理，目前关于这方面的研究并不多见。

本文以CEEMD为基础，分解含有间断信号和脉冲信号的仿真信号，分析并讨论CEEMD中存在的模态混叠和虚假分量问题，提出了一种改进的解相关CEEMD方法，有效地消除了仿真信号和爆破振动信号的模态混叠和虚假分量现象。

1 信号分解的CEEMD方法

1.1 CEEMD基本原理

CEEMD是针对EMD和EEMD (Ensemble Empirical Model Decomposition)在抑制信号分解中模态混叠问题的不足提出的改进方法。由于EMD仅通过信号局部的极值来不断的构造上下包络线，因此即使非常小的信号波动(如幅值很小的间歇成分)，也会造成EMD分解产生新的极值点，甚至新的IMF。这样产生的IMF即不满足唯一性。采用EMD分解无法解决模式混合问题。Wu等[17]对不同长度的白噪声做了大量的EMD分解，并提出了噪声辅助EMD分解的方式，即EEMD。EEMD可以在一定程度上抑制模态混叠的现象。EEMD是解决模式混合的一大创新，然而越来越多的研究指出EEMD存在着不足，比如残留噪声的叠加会导致严重的重构误差。为此，Yeh等[18]提出了噪声互补的改进算法CEEMD，主要步骤为：

对x(t)添加白噪声时，是成对的添加组，即

(1)

(2)

最终得到的每个模态分量为

(3)

IMFj是对x(t)做CEEMD分解所得到的固态模式函数，它一般应该有M=log2N-1个分量，即c1(t)、c2(t)、…、cM(t)，N是x(t)离散化以后的长度。

1.2 CEEMD的模态混叠和虚假分量现象

虽然CEEMD是一种较好的改进算法，可以很好的抑制信号分解过程中存在的模态混叠现象。但根据以往的文献和作者构造的仿真信号分解结果表明，当信号内噪声种类较多或信号突变点较多时，CEEMD的分解效果并不理想，同时在分解的过程中可能产生多个虚假分量。构造仿真信号D(t)说明CEEMD的分解效果。

D(t)=cos(2π·70t)+0.3cos(2π·30t)+

v1(t)+u(t)

(4)

式中：v1(t)是200 Hz高频间歇振荡信号，振荡区间为t=0.1～0.15 s、0.3～0.35 s，其余时间范围幅值为0，两个余弦信号频率分别为70 Hz和30 Hz。u(t)为一脉冲信号，表达式如式(5)

(5)

CEEMD分解中IMF=ceemd (Y,Nstd,NE,TNM)，其中Nstd是噪声标准偏差，取0.18，NE是集总次数，研究确定最佳集总平均次数为300次。TNM是分解获得的分量个数，取12。得到的IMF分量见图1。共分解出12个分量，仅分析前10个主分量，最后一个是余项rt(n)。

为抑制端点效应，CEEMD分解中采用了张迪[19]提出的适用于爆破振动信号的端点处理方法。从图1看到，各IMF分量在端点均未出现大的扰动。CEEMD将信号按频率从高到低依次分解出来，但分解过程中出现了模态混叠和过分解现象。IMF3和IMF4均为间歇振荡信号；IMF5为70 Hz的余弦分量，分解效果较好，信号内未出现模态混叠现象；IMF6为30 Hz的余弦分量，其中出现了模态混叠现象，对IMF6求频谱发现其混入了脉冲信号成分，且能量较小。由于分解过程中出现能量泄露，IMF5和IMF6的幅值比原始信号要小。IMF7和IMF8为脉冲分量，IMF9是否为脉冲分量需要做进一步的判断。

图1 D(t)经CEEMD分解得到的前10个IMF分量

对于易造成模态混叠的间歇信号和脉冲信号，容易出现两个甚至多个虚假分量，需要对各IMF分量进行虚假分量的判断。采用求取分量与原始分量的互相关系数和频谱判断图1中的虚假分量。IMF3、IMF4与原始分量的互相关系数分别为0.62和0.753，IMF7、IMF8和IMF9与原始分量的互相关系数分别为0.881、0.546和0.158。求取5个IMF分量的频谱，IMF3和IMF4的主频率均为200 Hz，并混有少量的70 Hz的成分，且IMF3的频谱幅值较小；IMF7和IMF8的主频率均在10～13 Hz，IMF8的频谱幅值较低，而IMF9的主频率为4 Hz。判断可知，IMF4和IMF7和原始单分量最为接近，为主要IMF分量，IMF3和IMF8与原始分量存在一定程度相关性，但能量相对较小，判断这两个分量为“伪”虚假分量。IMF9从互相关系数和主频率上和原始分量都不同，为真实虚假分量。

1.3 CEEMD中虚假分量的处理

目前没有虚假分量的统一处理方法，文中对分解得到的间歇信号和脉冲信号的虚假分量分别进行合成与直接剔除处理，并与原始单分量对比分析，合成结果及其频谱如图2和图3所示。对比结果如表1。

Fig.2 The original intermittency signal, mixture intermittency signal, intermittency signal by decomposed and

Fig.3 The original pulse signal, mixture pulse signal, pulse signal by decomposed and their frequency spectrum

从图2和图3和表1看出，间歇振荡信号的频率较高，最先被分离，因此合成后的IMF分量的幅值和原始单分量相差不大，能量损失较小。而由于端点效应和间歇信号内的奇异点的干扰，低频IMF的分解误差越来越大导致IMF7和IMF8的合成分量的能量损失较多，最大幅值约为原始分量的2/3。将“伪”虚假分量与主分量合并更有利于消除模态混叠，而且减轻了能量损失。

基于以上分析，由于爆破振动信号本身也是一种脉冲信号，而且含有多种高频噪声，CEEMD的分解效果不会特别理想。

2 改进的解相关CEEMD方法

上述分析可知，当信号内噪声种类和信号突变点较多时，仅用CEEMD分解还是会产生模态混叠和过分解现象，本质原因是分解得到的IMF不满足EMD的原始定义，即各个IMF没有完全正交，存在信息耦合问题。从以下两方面解决IMF正交性问题。

(1) 减少原始信号的奇异点、噪声的干扰，即对信号进行去噪处理。

(2) 将分解得到的IMF分量进行后处理，使得到的最终IMF分量满足正交性定义。

肖瑛等[20]提出了对EMD分解后进行解相关操作抑制各分量间的模态混叠问题。研究发现，对频率接近的平稳信号，采用解相关的EMD算法可以较好抑制模态混叠现象。但在分解含脉冲噪声和间歇信号噪声时会出现严重的模态混叠问题。对解相关的EMD方法进行改进，使其能更有效的分离出爆破振动信号中的间歇振荡信号和脉冲信号成分。

运用解相关方法处理CEEMD分解得到的各IMF分量的算法如下。

信号经EMD(CEEMD)分解后，得到n个模态分量dn(t)(n=1,2,…)。任意两相邻分量的相关系数rc为(以d1(t)和d2(t)为例)

(6)

式中:上标T代表信号的转置。rc·d2(t)为d1(t)和d2(t)中的共同部分，将这部分从d1(t)中剔除，即得到与d2(t)中不相关的部分dv1(t)

dv1(t)=d1(t)-rcd2(t)

(7)

由HHT分解特性可知，IMF分量的频率和幅值都随着阶数的增大逐渐减小。因此通过式(6)得到的rc不保证恒小于1，必须对两个相邻分量做归一化处理，即每个分量内所有数据均除以数列中的最大值，得到归一化分量。归一化后的相关系数rn与rc的关系为

(8)

n1和n2分别为d1(t)和d2(t)的归一化系数，归一化后得到的|rn|≤1，rn值越小，说明两个分量独立性越强；反之则表明分量间的共同部分越多。在归一化条件下，式(7)改写为

dv1(t)=d1(t)-rnd2(t)

(9)

若不对信号进行端点处理，由于端点处的飞翼现象，信号最大值可能出现在端点附近，此时归一化处理得到的归一化信号与实际信号的误差相差较大，因此必须对信号进行端点处理后再进行归一化计算。采用张迪提出的适合于爆破振动信号的端点处理方法。

基于上述算法，提出了改进的解相关CEEMD方法，步骤如下：

(1) 对原始信号端点处理，并进行CEEMD分解。

(2) 利用互相关系数和频谱判断虚假分量，将虚假分量和主分量相加，得到新的合成分量。

(3) 将噪声分量和含突变成分的前几阶分量去除，通过互相关系数法判断第一个主分量，从这一个分量开始依次设为d1(t)、d2(t)、d3(t)…。

(4) 求取d1(t)和d2(t)的最大幅值并做归一化处理，保证相关系数|rc|≤1。根据式(6)对d1(t)和d2(t)做解相关，求取两个分量的相关系数，并根据式(9)得到d1(t)与d2(t)的不相关分量dv1(t)。

(5) 得到的解相关后的信号余量dr(t)为

dr(t)=d1(t)+d2(t)-dv(t)

(10)

对dr(t)再次进行CEEMD分解，根据步骤(2)和(3)得到新的模态分量d11(t)、d22(t)、d33(t)…。将dv1(t)与d22(t)继续做解相关处理，得到不相关分量dv2(t)。

(6) 重复(2)～(5)，直到混叠部分消除为止。其余的IMF分量参照前两个分量的计算方法逐一分离。

由于CEEMD存在能量泄露，得到的IMF分量小于原始分量，因此分解次数并不是越多越好。针对不同信号类型和模态混叠程度，选择合理的计算次数。由于CEEMD的分解效果比EMD有着较大的提升，不需要对每个IMF分量都进行解相关计算。对得到的各个IMF分量求取频谱，仅对包含多个频带成分的IMF分量进行解相关操作，有利于对工程信号处理的实时性。

3 仿真分析

采用改进的解相关CEEMD方法对D(t)求解，观察其对模态混叠和虚假分量问题的处理效果，得到的4个IMF分量，采用傅里叶变换求取频谱见图4。

(a) IMF分量

(b) 频谱

间歇振荡信号的振荡部分是由高频余弦分量构成，本质上也是平稳信号，因此使用解相关计算后的分解效果较好，从频谱图中看到，IMF1和IMF2内已经没有模态混叠现象，由于使用的解相关次数较少，两个信号的能量泄露也较小。对于IMF3和IMF4，脉冲信号属于非平稳信号，分解较困难。两个分量经多次解相关和CEEMD计算后，能量泄露较IMF1和IMF2多，而且IMF3和IMF4内还存在少量模态混叠现象，但模态混叠部分能量较小，在频谱图中几乎无法观察到。对比图4和图1，经解相关计算后得到的30 Hz余弦分量和脉冲分量与仅通过CEEMD分解得到的两个分量相比，模态混叠现象得到极大改善，能量泄露也较小，因此信号经CEEMD分解后再进行解相关操作，分解的效果最好。

为了更直观的评价改进的解相关计算对消除模态混叠和虚假分量问题的效果，对D(t)分别进行EMD、CEEMD和解相关CEEMD方法的分解，根据EMD、EEMD和CEEMD定义，得到的各IMF分量不可避免的存在一定程度的能量泄露，计算各IMF分量与原始单分量的能量差值也可以判断抑制效果。对分解结果和原始分量计算相关系数ρi和能量均方根差值γm。

原始信号X(t)的均方根值SO

(11)

式中，N为X(t)的采样点数。若X(t)经CEEMD分解得到k个IMF，趋势项rn(t)为第k+1个IMF分量，分解后各分量的均方根值SR为

(12)

ci(t)为分解得到的各阶IMF。定义γm为能量均方根差值，即分解前后的均方根差值比值

(13)

γm值越小，说明得到的IMF分量越接近原始信号，分解效果越好；反之效果越差。计算结果见表2，其中D-CEEMD指端点处理、解相关计算和CEEMD联合处理方法。

表2 IMF分量和原始分量的相关系数和能量均方根差值

Tab.2 The correlation coefficient and the difference value of energy root mean square between IMFs components and the original components

评价方法模态分量EMDCEEMDD-CEEMDIMF10.0890.760.929相关IMF20.9240.9990.998系数IMF30.6170.8890.985IMF40.7380.8730.978IMF10.860.0940.033能量IMF20.030.1260.116差值法IMF30.720.0980.122IMF40.160.5030.103

EMD分解中，4个信号主要分布在后3个IMF分量中，计算的ρi和γm没有规律性。从CEEMD和改进的解相关CEEMD方法得到的ρi和γm看到，改进方法得到的4个分量的精度明显高于只用CEEMD分解的计算精度，相关系数均在0.92以上，且利用能量均方根差值得到的结果均比CEEMD的小。从能量差值上也可以看到，两种方法得到的IMF的γm值逐渐增大，这是由于分解次数的增加，能量泄露也逐渐增大导致。

4 爆破振动信号分析

图5是深圳某地铁隧道爆破振动信号，采用普通毫秒延期雷管起爆。图中看出，爆破前期，由于围岩夹制作用和较差自由面条件，质点振动速度往往较大，前期段位微差间隔较小，相邻段位之间波形会发生混叠，随后期段位间隔增大，段与段之间基本可看成独立爆破。

使用CEEMD方法和改进的解相关CEEMD方法处理上述爆破振动信号，CEEMD分解得到12个IMF分量，列出前10个分量，最后一个为趋势项。对CEEMD方法和改进的CEEMD方法处理得到的波形使用傅里叶变换求频谱图，结果如图6和图7所示。

图5 典型的爆破振动信号

如图6，前两个IMF分量频率最高、但所占能量最小，为现场环境中高频噪声，IMF3～IMF7分量频率逐渐降低，但相比IMF1，幅值很高，包含信号大部分能量，为原始信号的主频段。频谱图中看出IMF3～IMF7之间模态混叠现象最严重，求取各分量与原始信号相关系数，IMF3～IMF10分别为0.76、0.79、0.39、0.24、0.15、0. 01、0.001 5、0.003 6。对比相关系数和频谱图可见IMF3～IMF7为主要分量，IMF8为IMF7的“伪”虚假分量，IMF9为真实虚假分量。使用改进的解相关CEEMD方法，将虚假分量与主分量相加并对各分量进行处理后，剔除噪声分量和趋势项，共得到5个IMF分量。如图7，可以清楚看出，信号内按频率从高到低主要包含100 Hz、77 Hz、62 Hz、43 Hz、25 Hz、11 Hz和6 Hz共7种频率成分。除了100 Hz和62 Hz的成分未分量出来外，其它5个IMF中的主频率能量均较高，非频率成分得到较好抑制。改进方法得到的各IMF振速峰值分别为0.11，0.05，0.021，0.013，0.007。主频基本分布在100 Hz以下，较之CEEMD得到的分量峰值(0.1，0.04，0.018，0.008 7，0.004 7)均有一定程度提高，对比频谱图中主频对应的峰值亦是如此，可见改进的CEEMD方法在处理模态混叠的同时可以减轻一部分能量损失。