应用自适应带宽信号的BS-EMD混叠消除

蒋永华， 焦卫东， 李荣强， 唐 超， 郑佳佳， 蔡建程

(浙江师范大学 精密机械研究所，浙江 金华 321004)

在工程实际中，旋转机械的故障振动信号大部分是非平稳、非线性信号。对于非平稳信号，采用传统的Fourier变换无法满足对此类信号的分析需求[1]。EMD[2-3]是美籍华人Huang提出的一种非平稳信号自适应分解方法[4]，其本质是对信号进行平稳化处理，把复杂的信号分解成有限个本征模态分量(Intrinsic Mode Function, IMF)[5]。它是一种无需任何先验知识的自适应时频分析方法，能够很好地体现非平稳信号的局部特性，已应用于信号检测、桥梁故障监测、医学信号监测等众多工程领域[6]。但是，EMD方法仍然存在许多需要改进的地方，如端点效应、模态混叠等问题。模态混叠现象最早由Huang通过对含有间断信号的EMD分解时发现，它严重影响了EMD的应用效果，阻碍了EMD在工程中的推广应用。

针对模态混叠问题，Huang提出通过预先设定尺度(极值点间距离)上限实现间断判别来抑制模态混叠的思想，但是并没给出详细的处理方法，而且需要预先了解信号的基本特点，难以实现；马文朋等[7]针对EMD在工程应用中存在的端点效应和模式混叠问题，提出了一种改进的EMD方法;Rilling等[8]提出局部EMD方法，在一定程度上改善了模态混叠，但在窗函数的宽度选择上会带入主观因素；郑源等[9]提出了一种改进的EMD-ICA方法，以IMF的能量波动系数为指标设定阈值，判断是否发生模态混叠，排除虚假IMF；禹丹江等[10]通过设置间断频率来避免模态混叠，使每一个IMF表示结构的某一阶固有模态，但是间断频率作为筛分过程的一个附件条件会对分解结果造成未知的偏差；Ryan等[11]提出了掩膜信号法，该方法对于小频率比混合信号引起的模态混叠有较好效果，但在频谱图上仍能看到混叠分量存在；对此，赵玲等提出一种改进的掩膜信号法，以能量为基础对掩膜信号的选择进行了改进，但是对于不同的混合信号，掩膜信号的幅值和频率确定具有一定难度；Wu等[12]提出了总体平均经验模态分解(EEMD)方法，对间断信号引起的模态混叠具有很好的抑制效果，但是对于小频率比信号产生的模态混叠效果不理想。另外，EEMD需要进行几十至上百次的EMD运算，效率很低。同时，为了使EEMD实现理想的分解结果，预先需要知道信号的信噪比；肖瑛等[13]提出解相关EMD方法，通过在EMD过程中嵌入解相关操作来解决模态混叠问题；陈建国等[14]采用独立分量分析(ICA)来抑制模态混叠；汤宝平等[15]对信号进行形态滤波后再进行峭度最大化的ICA来实现混叠成分的分离，但是ICA的幅度不确定性会限制对信号能量信息的分析判断；郑近德等[16]针对EMD存在的包络误差，模态混叠，端点效应等缺陷，提出了一种改进的HHT方法(IHHT)；汤宝平等采用形态奇异值分解对信号进行预处理来消除模态混叠和边界效应，提高EMD分解质量；曹莹等提出了基于形态滤波预处理与端点延拓相结合的方法来抑制模态混叠；胡爱军等提出了一种高频谐波加入法(HFHA-EMD)，通过加入高频谐波后再进行EMD分解来消除模态混叠。高频谐波频率根据采样频率选择在分析频率的上限附近，幅值取接近原始信号的幅值，据此来确定高频谐波的幅值和频率具有一定难度。本文利用EMD的二进带通滤波特性和总是先分离高频分量的特点，基于高频谐波加入法原理并对高频谐波信号的构造进行改进，根据被分析信号的BS-EMD得到的第一个IMF分量确定带宽限制频率和带宽限制幅值，从而构造出自适应带宽限制信号。通过在原始信号中添加自适应带宽限制信号来改变EMD带通滤波器的中心频率，再进行BS-EMD分解消除模态混叠。

1 EMD模态混叠现象及产生机理

EMD的分解过程依赖信号本身包含的变化信息，是一种自适应的信号分解方法。但是，在某些情况下，EMD得到的IMF会产生模态混叠现象。模态混叠具体有以下两种表现形式：①单个IMF中包含不同尺度或频率的多个信号；②同一尺度或频率的信号被分解到多个不同的IMF中。表现为相邻两个IMF波形混叠，相互影响，难以辨别，致使后续的时频分布混淆，也导致混叠的IMF缺乏足够的物理意义。

产生模态混叠的原因主要包括以下几个方面：①混合信号的组合分量的频率过于接近(不满足条件f1/f2<0.5)；②信号中含有小幅值高频间断信号；③信号中存在脉冲干扰和噪声。

对于混合信号的组合分量的频率过于接近引起的模态混叠现象。究其原因，是因为信号所包含的两个分量频率太过接近时，EMD很容易将这个信号解释成由一定调制程度的IMF组成，这就偏离了原信号的本质，因此，EMD不能将它们完全分解到两个IMF中，造成模态混叠。

对于信号中含有异常成分(小幅值高频间断信号)引起的模态混叠现象。究其原因，是因为异常事件(如小幅值高频间断信号)的存在造成了局部极值点的异常分布，为保证包络线的光滑性，三次样条插值产生的包络不可避免会产生失真，比如出现过冲和欠冲现象，从而使得分解结果中包含了信号和异常事件，造成模态混叠。可以说，在间断信号处包络线的过冲和欠冲现象是产生模态混叠的重要原因之一。

对于信号中存在脉冲干扰和噪声引起的模态混叠现象，作者已做了相关研究，限于篇幅，不再赘述，详见汤宝平等的研究结论。

2 模态混叠消除方法

2.1 BS-EMD

曲线拟合是EMD中的一个关键问题，直接影响到EMD结果。三次样条插值算法容易产生过冲和欠冲现象。而在间断信号处包络线的过冲和欠冲现象则是产生模态混叠的重要原因之一。B样条曲线不仅保持了三次样条曲线的优点，而且还具有良好的局部性和连续性。此外，B样条插值直接拟合均值曲线，无需计算上、下包络线，因此减少了计算量，加快了拟合均值曲线的速度，提高了EMD分解的时效性。因此，在EMD的曲线拟合方法上，B样条插值算法是一种更为理想的选择。故本文研究基于三次B样条曲线插值算法的EMD，简称BS-EMD。BS-EMD分解过程如下：

(1) 求得信号x(t)的所有极值点(不区分极大值点和极小值点)；

(4) 计算x(t)与m(t)的差， 即x(t)-m(t)=h1(t)；

(5) 判断是否为第一个IMF分量c1(t)， 否则，重复以上步骤，直至得到第一个IMF分量；

(6)x(t)减去c1(t)作为新的原信号重复以上步骤n次， 可依次得到x(t)的第n个IMF分量和第n阶残余分量rn(t)；

(7) 当rn(t)成为一个单调函数不能再从中提取IMF分量时，循环结束，整个分解过程也终止。

用仿真算例说明BS-EMD优于三次样条EMD。采用两个频率接近的余弦信号叠加生成的信号，如下

x1(t)=cos(2π10t)+cos(2π14t)

(1)

信号由幅值为1、频率分别为f1=10 Hz和f2=14 Hz的两个余弦信号组成。采样频率为500 Hz，采样点数为500点。

显然，信号不满足条件f1/f2<0.5。 其时域波形如图1所示。 采用基于三次样条的EMD对x1(t)进行分解， 结果如图2所示， 其Hilbert谱如图3所示。

图2 仿真信号x1(t)三次样条EMD的分解结果Fig.2 Decomposition result of x1(t) by cubic spline based

图3 仿真信号x1(t)三次样条EMD的Hilbert谱Fig.3 Hilbert spectrum of x1(t) by cubic spline based EMD

从图2可知，分解出了4个IMF分量和一个余项。其中，IMF1和IMF2分别为14 Hz和10 Hz的两个余弦信号。但是，EMD没有将两个余弦分量完全分解到IMF1和IMF2中，10 Hz的IMF2分量的一部分信息被分解到了14 Hz的IMF1分量中。因此，EMD没有将它们完全分解到两个IMF中，分解结果无论是幅值还是频率都有误差，产生了模态混叠。从图3中可看出，两个信号分量的频率波形很大，且相互混叠，很难准确地分辨出两个余弦信号的频率，即模态混叠致使时频分布混淆。此外谱线的两端出现了明显的端点效应。

究其原因，是因为信号所包含的两个分量频率太过接近时，EMD很容易将其解释成由一定调制程度的IMF组成，这就偏离了原信号的本质，因此，EMD不能将它们完全分解到两个IMF中，造成模态混叠。

对仿真信号x1(t)采用BS-EMD方法进行分解，结果如图4所示，其Hilbert谱如图5所示。

图4 仿真信号x1(t)的BS-EMD分解结果Fig.4 Decomposition result of x1(t) by BS-EMD

图5 仿真信号x1(t) BS-EMD的Hilbert谱Fig.5 Hilbert spectrum of x1(t) by BS-EMD

对比图2与图4可知，尽管BS-EMD也存在着模态混叠现象，但仍可看出BS-EMD的结果要优于三次样条EMD。三次样条EMD的IMF1的幅值达到了2.5且波形相较于余弦信号而言更接近于调制信号，同时，IMF2的幅值却不足0.5且波形两端存在着明显的发散现象，此外，虚假分量IMF3的幅值也较大。相比较而言，BS-EMD的IMF1和IMF2更接近于真实的余弦信号，且IMF2端点处无明显发散现象。这是因为B样条插值算法是直接拟合均值曲线，避免了三次样条插值算法的过冲与欠冲现象，减小了EMD的分解误差，提高了EMD分解精度。

进一步地，对比图3和图5可知。虽然BS-EMD的两个信号的频率也存在着一定的波动，但是相比较三次样条EMD而言波动要小的多，尤其是第二个频率分量。从图5中已能够准确地区分出两个余弦信号的频率，即10 Hz和14 Hz。此外，谱线两端的端点效应也没三次样条那么明显。因此，基于BS-EMD的频率分辨能力也要优于三次样条EMD。虽然BS-EMD的Hilbert谱可以区分出两个余弦分量的频率，但是谱线并不是两条直线，而是具有一定的波动性，这是由于模态混叠造成的。

促销活动，主要是为了更多地促进本公司产品的销售，而在一定时间内对产品的营销。对于南通鹏越纺织有限公司而言，人员推销、直复营销等方式，是其使用较为频繁的促销方式。一方面，通过促销人员对潜在客户和直接客户的沟通交流，去推卖产品；另一方面，南通鹏越纺织有限公司的直复营销种类较多，其中电话营销是其运用最多的营销方式。通过电话的对接，与客户进行多次的明确的沟通交流，按照客户要求确定数量和产品，生产包装，并启动物流送达客户手中。

2.2 自适应带宽限制信号法

EMD分解是从高频到低频逐次分解的过程，含高频成分的IMF总是最先被分解出来。因此，本文利用EMD分解的二进带通滤波特性和总是先分离高频分量的特点，借鉴胡爱军等的高频谐波加入法(HFHA-EMD)原理，采用在原始信号中添加带宽限制信号的方法，来改变EMD带通滤波器中心频率。由于带宽限制信号频率较高，因此，EMD分解总是最先把这个带宽限制信号的所有信息分解出来，包含在第一阶IMF分量中，由于带宽限制信号已知，因此，将第一个IMF减去带宽限制信号便可得到原始信号真正的第一个IMF分量。若原始信号中含有异常成分，则第一阶IMF即为带宽限制信号和异常成分，而原始信号中的第一个真实分量将基本不会被分解到这阶IMF中，因此，可直接将这阶IMF去掉。

因此，如何构造一个合理的带宽限制信号是关键。胡爱军等建议构造高频谐波时，其频率根据采样频率选择在分析频率的上限附近，其幅值取接近原始信号的幅值，据此来确定高频谐波的幅值和频率具有一定难度。此外，这种方法构造的高频谐波信号没有充分利用被分析信号自身的特点，因此也不具有自适应性。本文对其进行改进，根据被分析信号自身的特点进行自适应构造带宽限制信号。具体构造方法如下：

(1) 对信号x(t)进行BS-EMD分解， 得到第一阶IMF分量c1(t)和一个余项r1(t)， 即：x(t)=c1(t)+r1(t)；

(2) 对c1(t)做Hilbert变换，根据式(2)计算平均频率favr：

(2)

式中：a1(t)和f1(t)分别为c1(t)瞬时幅值和瞬时频率。

(3) 根据式(3)得到带宽限制信号频率fbr：

fbr=αfavr

(3)

式中：α为带宽限制频率系数。

(4) 根据式(4)计算c1(t)的平均幅值aavr：

(4)

abr=βaavr

(5)

式中：β为带宽限制幅值系数。

(6) 根据fbr和abr构造带宽限制信号：

xbr(t)=abrcos(2πfbrt)

(6)

由于带宽限制信号频率fbr和带宽限制信号幅值abr都是根据被分析信号的第一个IMF分量c1(t)计算得到， 而EMD本身是一种自适应的信号分解方法，因此， 由fbr和abr构造的带宽限制信号也具有自适应性，称之为自适应带宽限制信号。

带宽限制频率系数α的取值，对于频率相近信号之间相互作用引起的模态混叠问题， 信号c1(t)的平均频率favr介于几个真实分量频率之间， 因此系数α可适当取大一些，使得最终构造的带宽限制信号频率高于信号自身最高频率分量的频率，实验分析表明α=1.4为宜(此值附近适当调整取值，确定最佳取值，以下类同)；对于异常成分(小幅值高频间断信号)引起的模态混叠问题，由于异常成分是高频信号(通常情况下都是比信号大好多倍)，使得信号c1(t)的平均频率favr通常比较高，因此，α可适当取小一些，取α=0.5为宜。带宽限制幅值系数β的取值与α类似，但是，实验表明EMD方法对于带宽限制信号的幅值取值的敏感性远不及频率的取值，因此，取β=1.6为宜。

3 仿真分析与比较

(1) 仿真分析一：对含频率相近分量的信号分解

考察3.1章节中的仿真信号x1(t)。首先根据对分析信号x1(t)进行BS-EMD分解，再构造自适应带宽限制信号。由式(2)和式(4)可得favr=13.34 Hz、aavr=1.22。 取α=1.4、β=1.6则fbr=19、abr=2。因此，自适应带宽限制信号xbr(t)=2cos(2π19t)， 如图6所示。采用本文所提的基于自适应带宽限制信号的BS-EMD方法方法对x1(t)进行分解， 其BS-EMD分解结果和Hilbert时频谱分别如图7和8所示。

图6 自适应带宽限制信号的时域波形Fig.6 The waveform of adaptive bandwidth constrained

从图7和图8可看出，两个分量已经被完全分解到IMF1和IMF2中，频率分别为14 Hz和10 Hz，幅值范围均为[-1,1]，可知分解结果准确地反应了原始信号的特征。对比图5和图6可以说明，基于自适应带宽限制信号的BS-EMD方法可以有效抑制混合信号的组合分量的频率过于接近引起的模态混叠现象。

图7 基于自适应带宽限制信号的BS-EMD分解结果Fig.7 Decomposition result of x2(t) by the proposed method

图8 基于自适应带宽限制信号的BS-EMD的Hilbert谱Fig.8 Hilbert spectrum of x2(t) by the proposed method

(2) 仿真分析二：对含异常成分(小幅值高频间断信号)的信号分解

采用参考文献[4]中的示例信号，由下式给出

x2(t)=sin(2π8t)+i(t)

(7)

式中：i(t)为幅值0.1， 频率160 Hz的三段间断正弦信号(分别为[287,337]、[537,587]、[787,837])。

采样频率为2 000 Hz，采样点数为1 200点，其时域波形如图9所示。参考胡爱军等给出采用三次样条EMD的x2(t)分解结果，限于篇幅，本文不再重复给出。但在后续章节的其他内容中会给出相应的处理结果。

图9 仿真信号x2(t)的时域波形Fig.9 The waveform of simulation

由式(2)和式(4)可得favr=268.33 Hz、aavr=0.58。 取α=0.5、β=1.6则fbr=134、abr=0.9。 构造自适应带宽限制信号xbr(t)=0.9cos(2π134t)， 如图10所示。对x2(t)进行分解，BS-EMD结果如图11所示。

图10 自适应带宽限制信号的时域波形Fig.10 The waveform of adaptive bandwidth constrained signal

图11 基于自适应带宽限制信号的BS-EMD分解结果Fig.11 Decomposition result of x2(t) by the proposed method

从图11可知，间断信号和带宽限制信号都已被分解到IMF1分量中，且间断信号完全被带宽限制信号所“淹没”，而8 Hz的低频正弦信号则基本被完全分解到IMF2分量中。将IMF1分量减去图11的带宽限制信号，则可以得到原始信号中的间断正弦信号，结果如图12所示。因此，此仿真分析说明基于自适应带宽限制信号的BS-EMD方法可以有效抑制异常成分(小幅值高频间断信号)引起的模态混叠。

图12 从IMF1中提取出的间断正弦信号Fig.12 Intermittent sinusoidal signal obtained from IMF1

对比胡爱军等的仿真计算(图4)可知，两种方法在抑制异常成分(小幅值高频间断信号)引起的模态混叠方面都取得了较为理想的效果。为了不失一般性，胡爱军等采用的仿真信号的三个间断信号偏离了正弦信号的正峰位置。而为了便于与EMD和BS-EMD比较，也限于篇幅，本文没有采用与胡爱军等相同的信号进行比较，应该说本文所用信号的情况更加复杂恶劣，也更易引起模态混叠现象。另外，本文在构造自适应带宽限制信号上更加明确可行，还充分利用了被分析信号自身的特点，因此也更具自适应性。

胡爱军等将HFHA-EMD与EEMD进行了比较，证明了其比EEMD更有优势，因此，本文不再与EEMD做类似比较。此外，胡爱军等没有考虑到混合信号的组合分量的频率过于接近引起的模态混叠现象，而本文则对此展开了较为翔实的论证，因此，本文所做工作与胡爱军等的研究工作互为补充、互为验证，也为EMD模态混叠问题做了一个较为全面的分析和讨论。

4 实际信号分析与比较

现以一实际信号来验证本文方法的有效性。图13为使用加速度传感器采集的某旋转机械转子测量面缺陷信号，采样频率为2 000 Hz，工频为46.01 Hz。

图13 旋转机械转子测量面缺陷信号Fig.13 Signal of the Rotor with measurement surface defects

图14为该信号的三次样条EMD分解结果，从图中可以看出，信号时域波形含有脉冲、噪声等异常成分， IMF3与IMF4之间存在着明显的调制特征，表明出现了模态混叠现象。图15为IMF3与IMF4两个分量的功率谱，图16(a)中，IMF3分量中不仅出现了二倍频，还出现了功率较大的三倍频和四倍频；图15(b)中，IMF4分量除了工频成分及边频外，还出现了二倍频，模态混叠现象严重。EMD分解结果不能体现设备的真实状态。

图14 转子信号的三次样条EMD的分解结果Fig.14 Decomposition result of the rotor signal based on EMD

利用本文方法进行分析，由式(2)和式(4)得该转子测量面缺陷信号favr=488.73 Hz、aavr=0.68。 取α=0.6(0.5附近适当调整取值，确定最佳取值为0.6)、β=1.6(与α取值同理)， 则fbr=293、abr=1.1。 构造自适应带宽限制信号xbr(t)=1.1cos(2π293t)， 分解结果如图16所示。可以看出，各IMF分量没有明显的调制现象，表明模态混叠现象得到了有效抑制。相较于三次样条EMD得到的9个IMF分量，本文方法仅有7个分量，说明本文方法能更好的消除虚假分量。

图15 转子信号三次样条EMD分解结果的功率谱Fig.15 Power spectrum of EMD of the rotor signal

图16 转子信号基于自适应带宽限制信号的BS-EMD结果Fig.16 Decomposition result of the rotor signal x2(t)by the proposed method

取α=0.5、β=1.6(仅改变α取值，减小α)，则fbr=244、abr=1.1。 构造自适应带宽限制信号xbr(t)=1.1cos(2π244t)， 分解结果如图18所示。

相较于图16结果可知，分解质量没有图16高，IMF分量存在一定的调制现象，模态混叠现象虽得到了一定抑制但效果不够理想。此外，分解得到8个分量，虽比三次样条EMD得到的9个分量要好，但是比图16的7个分量要差。同样，对IMF3和IMF4两个分量做功率谱分析，结果如图19所示。图19(a)中，IMF3分量除了含有二倍频成分，还出现了工频与三倍频分量；图19(b)中，IMF4分量出现了一定的边频。总体而言，比图15所示的三次样条EMD结果好，但是不如图17所示的最佳参数BS-EMD结果好。

图17 转子信号基于本文方法EMD分解结果的功率谱Fig.17 Power spectrum of decomposition result of the rotorsignal based on the proposed method

图18 转子信号基于自适应带宽限制信号的BS-EMD结果Fig.18 Decomposition result of the rotor signal by the proposedmethod with different adaptive bandwidth constrained signal

图19 转子信号基于本文方法EMD分解结果的功率谱Fig.19 Power spectrum of decomposition result of the rotorsignal based on the proposed method

取α=0.7、β=1.6(仅改变α取值，增大α)，则fbr=342、abr=1.1。 构造自适应带宽限制信号xbr(t)=1.1cos(2π342t)， 分解结果如图20所示。IMF3和IMF4两个分量的功率谱如图21所示。对比分析可知，取α=0.7与α=0.5时情况类似，优于三次样条EMD方法效果，但是不如本文取值结果。

β的取值与α类似，限于篇幅，不再赘述。

图20 转子信号基于自适应带宽限制信号的BS-EMD结果Fig.20 Decomposition result of the rotor signal by the proposedmethod with different adaptive bandwidth constrained signal

图21 转子信号基于本文方法EMD分解结果的功率谱Fig.21 Power spectrum of decomposition result of therotor signal based on the proposed method

5 结 论

(1) 针对EMD的模态混叠现象，本文对模态混叠的具体表现形式、引起模态混叠的主要原因进行了归纳，对模态混叠的产生机理进行了分析。

(2) 借鉴高频谐波加入法原理，提出了基于自适应带宽限制信号的BS-EMD方法。充分利用被分析信号自身的特点，构造出自适应带宽限制信号。仿真分析和实际含复杂异常事件的转子故障信号分析验证了所提方法的有效性和可行性。

(3) 与EMD、BS-EMD的对比分析验证了本文方法的优越性。与HFHA-EMD的对比仿真分析表明，两种方法都可以有效消除模态混叠现象，但是本文方法在构造自适应带宽限制信号上更加明确可行，也更具自适应性。