☰
经验模态分解(EMD)原理与Python实践:非平稳信号的自适应时频分析
2026/10/13 7:56:53 网站建设 项目流程

前阵子做一个振动信号分析任务,数据里既有恒定转速下的周期性冲击,又有一小段载荷突变带来的非平稳过渡过程。傅里叶谱能告诉我冲击频率在哪里,但我最想知道的那个“从低频过渡到高频的具体时间点”却完全看不到。把信号切片做短时傅里叶变换,勉强能看到趋势,可频率分辨率又掉得厉害。最后解决问题的正是经验模态分解(EMD):它把原始信号按自身时间尺度自适应地拆成若干本征模态函数(IMF),再结合希尔伯特变换做瞬时频率、边际谱和相关分析,整个流程才算真正有了一个可复现的完整示例。

这篇内容不是教科书式复述,而是从实际使用者的视角,把EMD的来龙去脉、筛分原理、完整Python示例、以及我在实践中反复踩过的坑一次性讲清楚。适合正在做信号处理、故障诊断、气象数据分析或金融时间序列建模的工程师和研究人员,尤其适合那些发现常规频域分析不够用、想引入自适应分解方法的读者。

1. EMD出场之前:傅里叶和小波的“别扭之处”

1.1 平稳性假设:最常被忽略的那个前提

很多人做频谱分析时第一反应是直接上傅里叶变换,却忽略了一个关键前提:傅里叶变换默认信号是平稳的,至少是分段平稳的。平稳意味着统计特性不随时间变化,频率成分始终一致。但实际采集到的信号往往不是这样:机械设备启停机过程中的振动信号、风电塔筒在变风速下的响应、脑电信号中的瞬态波形,都是典型的非平稳、非线性信号。

对于非平稳信号,傅里叶变换能给出一个“平均意义上的频谱”,告诉你整个时间段内哪些频率占主导,却无法回答“这个频率成分是什么时候出现的”。比如一个信号前半段是5Hz正弦波,后半段变成20Hz正弦波,傅里叶谱上会同时出现5Hz和20Hz的谱峰,但你完全看不出两者的时间顺序。短时傅里叶变换通过加窗勉强解决了这个问题,却又引入了新的矛盾:窗越窄时间分辨率越高,频率分辨率就越差;窗越宽频率分辨率好,但时间信息又变得模糊。这个不确定性原理是窗函数方法绕不过去的坎。

1.2 小波基函数:选择困难症的根源

小波变换比短时傅里叶进步了一大截,它通过平移和伸缩基函数来同时获得时间与频率信息。但小波分析有一个让人头疼的地方:必须提前选定小波基函数和分解层数。选db4还是sym8?分解到第几层合适?这些选择会直接影响分析结果,而且不同的基函数对同一段信号的分解结果差异可能很大。实际做项目时,换一种小波基就等于换了一种分析视角,结论的稳定性很难保证。

这还不是最麻烦的。小波变换本质上是将信号与固定的基函数做内积,一旦基函数内在的振荡形态和你的信号特征不匹配,分解结果就会失真。换句话说,小波依然是“用别人的尺子量你的信号”,尺子本身的刻度是固定的。对信号先验信息不足的场景,这种依赖人为设定的方法经常让人无从下手。

1.3 EMD究竟带来了什么变化

经验模态分解(Empirical Mode Decomposition, EMD)的出现,把“选基函数”这个环节整个拿掉了。它不需要预先指定任何基函数,而是从信号自身出发,基于局部极值的时间尺度来提取本征模态函数(Intrinsic Mode Function, IMF)。每个IMF都是一个单一振荡模式的近似表示,可以窄带地反映信号在不同时间尺度上的内在特征。

这种自适应性带来的直接好处是:同一套算法、同一组参数,可以不加修改地处理振动信号、脑电信号、气象观测序列甚至金融价格序列。你不用先把信号分类再决定用哪把“尺子”,EMD把分解的过程交给了数据和算法本身。当然,这不是说EMD没有任何参数需要设置,但相比小波基函数的选择,EMD的参数更多体现在停止准则和端点处理策略上,对结果的影响方式也完全不同。这一点后面在实操章节里我会展开讲。

2. 拆解筛分过程:从包络线到本征模态函数

2.1 什么是合格的IMF:两个硬性条件

要理解EMD,必须先理解IMF的定义。IMF需要同时满足两个条件:

  • 在整个数据序列中,极值点(局部极大值和局部极小值)的数量与过零点的数量相等,或至多相差1;
  • 在任意时间点,由局部极大值定义的上包络线和由局部极小值定义的下包络线的均值必须为0。

这两个条件的目的,是让IMF可以被看作一个近似窄带的振荡信号,从而可以直接对其做希尔伯特变换来获得有物理意义的瞬时频率。需要注意的是,第二个条件中的“均值为0”在实际计算中只能近似满足,通常用上下包络均值足够接近0来判断,而不是绝对等于0。

打个比方:把复杂信号想象成一杯由果汁、果肉和冰粒组成的混合饮品,IMF就相当于从中分层提取出的不同成分,每一层都有相对单纯的形态。上下包络均值接近0,本质上是在要求每一个IMF都是关于时间轴对称的振荡,这样后续求瞬时频率时才不会出现负值或剧烈跳变。

2.2 筛分是怎么一步步进行的

EMD的核心算法是“筛分”(sifting),大致流程如下:

  1. 找出原始信号x(t)的全部局部极大值点和局部极小值点;
  2. 用三次样条曲线分别拟合极大值点和极小值点,形成上包络线u(t)和下包络线l(t);
  3. 计算两条包络线的均值m(t) = (u(t)+l(t))/2;
  4. 从原始信号中减去均值包络,得到候选分量h(t) = x(t) - m(t);
  5. 检查h(t)是否满足IMF的两个条件;如果不满足,就把h(t)当作新的原始信号,重复步骤1到4;
  6. 当h(t)满足IMF条件后,把它记作第一个IMF,记为imf1;
  7. 用原始信号减去imf1,得到残差r(t),对r(t)继续重复上述筛分过程,依次提取imf2、imf3……直到残差变成单调函数或小于预设阈值为止。

第4步是整个方法最直观的地方:每一次减去的“均值包络”近似代表信号中的慢变成分,剩下的h(t)则保留了一个相对快速的振荡模式。反复迭代的目的,是为了让h(t)逐渐逼近一个对称的、单一振荡特征的分量。

2.3 停止准则:筛到什么时候算完

筛分次数不是越多越好。如果一直“筛”下去,分量会逐渐趋于纯正弦波,原始信号的幅度和形态信息会被过度抹平,甚至产生没有物理意义的分量。因此必须设置停止准则。最经典的做法是计算连续两次筛分结果的标准差:

S D = sum(|h_{k-1}(t) - h_k(t)|^2) / sum(h_{k-1}(t)^2)

当SD落在0.2到0.3之间时中止当前IMF的提取。这个阈值既保证了包络均值的局部对称性足够好,又避免筛分过度。此外还有基于筛分次数上限的停机策略,比如限制单个IMF的最大迭代次数为100次,防止程序死循环。

在实际使用Python中的EMD库时,这些停止准则通常已经封装成了默认参数,但理解背后的逻辑仍然重要。因为当你遇到分解结果异常“破碎”或“过度平滑”时,首先应该调整的就是筛分停止阈值,而不是去怀疑库本身。

3. 可直接复跑的Python完整示例

3.1 安装依赖与合成测试信号

这里我用Python生态里比较常用的PyEMD库来演示,安装命令很简单:

pip install EMD

这个库实际导入时的模块名是PyEMD,所以代码里写的是from PyEMD import EMD。除了EMD库本身,还需要numpy、scipy和matplotlib,这几个用pip一并装上就行。

先构造一个非平稳多分量测试信号:一个5Hz基频、一个25Hz谐波分量,再叠加一段只在后半段出现的瞬态冲击。这个信号本身并不复杂,但已经超出了简单傅里叶分析能完美处理的范畴。

import numpy as np import matplotlib.pyplot as plt fs = 1000 t = np.linspace(0, 2, 2 * fs) signal = ( np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 25 * t) + 0.2 * np.sin(2 * np.pi * 55 * t) ) # 后半段加入瞬态冲击 signal[int(fs * 1.0):] += 0.6 * np.sin(2 * np.pi * 80 * t[int(fs * 1.0):])

采样率设为1000Hz,数据长度2秒。为什么用这个采样率?因为信号最高的60Hz分量加上瞬态成分后,200Hz以内的频率成分都值得关注,1000Hz采样率足以保证奈奎斯特频率远高于分析需求,不会有混叠风险。

3.2 EMD分解主流程与画图

把信号喂给EMD之后,会得到一个二维数组,每一行就是一个IMF分量。

from PyEMD import EMD emd = EMD() imfs = emd(signal) fig, axes = plt.subplots(imfs.shape[0] + 1, 1, figsize=(10, 8), sharex=True) axes[0].plot(t, signal, 'k') axes[0].set_ylabel('Original') for i, imf in enumerate(imfs): axes[i + 1].plot(t, imf, 'k') axes[i + 1].set_ylabel(f'IMF{i + 1}') plt.xlabel('Time [s]') plt.tight_layout() plt.show()

分解完成后,你会发现前几个IMF分别对应信号的快速振荡成分,越靠后的IMF振荡越慢,最后一个分量往往是趋势项或残差。这里有一个经验性的观察方法:IMF1通常对应最高频的细节分量,最后几个IMF对应低频趋势,这个规律在大多数实测信号中都成立,但不要机械地认为IMF1就一定是最重要的分量,它的物理意义需要结合信号特征去判断。

3.3 结果解读:每一个IMF对应什么

判断IMF对应什么物理成分,最直接的方法是看它的频谱或直接观察时域波形。上面这个测试信号中,正常应该能分解出80Hz瞬态冲击对应的IMF、55Hz和25Hz分量对应的IMF、5Hz基频对应的IMF,以及一个接近0的趋势项。

注意,如果是理想合成信号,EMD可以分得很干净;但真实信号总伴随噪声和幅值波动,某些频率接近的分量可能会被分到同一个IMF里。此时不要急着下结论,而应结合希尔伯特谱或滤波手段进一步确认。你还可以计算每个IMF与原始信号的相关系数,如果某个IMF和原始信号的相关系数极低且波形看起来杂乱无章,那它大概率是分解产生的伪分量。

3.4 继续到希尔伯特谱和边际谱

EMD的典型伙伴是希尔伯特变换。对每个IMF做希尔伯特变换,可以得到解析信号,进而提取瞬时幅度和瞬时频率:

from scipy.signal import hilbert analytic = hilbert(imfs[0]) inst_amp = np.abs(analytic) inst_phase = np.unwrap(np.angle(analytic)) inst_freq = np.diff(inst_phase) * fs / (2 * np.pi)

对这组瞬时频率序列做统计,可以得到瞬时频率随时间变化的曲线。把所有IMF的瞬时频率和瞬时幅度画在同一张图上,就构成了希尔伯特谱,它比普通时频图更适合表达非平稳信号的时间-频率-能量分布。边际谱则是把希尔伯特谱沿时间轴做积分,得到类似“频谱”的表示,但它保留了信号局部变化信息,比傅里叶频谱更锐利。

4. 相关分析的价值:瞬时频率、边际谱和特征提取

4.1 瞬时频率的计算逻辑

很多人第一次接触“瞬时频率”这个概念会有点别扭:傅里叶分析里频率是一个全局量,它和正弦波的周期一一对应;瞬时频率却允许频率在每一个时间点都被定义,它来自解析信号相位的导数。

对一个实信号x(t),希尔伯特变换构造解析信号z(t) = x(t) + j·H[x(t)],其中H表示希尔伯特变换。解析信号的相位θ(t)对时间求导,再除以2π,就得到瞬时频率。IMF之所以适合做这个计算,是因为它的上下包络均值对称、振荡模式接近单一正弦,相位曲线单调增长,求出的瞬时频率才有明确的物理意义。对非IMF信号直接求瞬时频率,会出现大幅跳变甚至负频率,这就是为什么Hilbert-Huang变换要求先做EMD再算瞬时频率,顺序不能颠倒。

4.2 边际谱和傅里叶谱到底哪里不一样

傅里叶谱描述的是信号中“存在哪些频率成分”,边际谱描述的则是“这些频率成分在时频平面上累积了多少能量”。两者在平稳纯正弦信号上的结果很接近,但在非平稳信号上差异明显:边际谱的频率分辨率不受窗函数限制,不会出现傅里叶谱中常见的能量泄漏拖尾。

实际操作中,我经常同时画出傅里叶谱和边际谱来对照。当发现某个谱峰在傅里叶谱中很模糊、但在边际谱中非常尖锐时,往往就是非平稳特征在起作用。比如一个频率随时间周期性变化的信号,傅里叶谱会把它摊平成一个宽峰,而边际谱能把它集中到实际的频率范围附近,这对诊断变速工况下的机械故障非常有帮助。

4.3 从IMF到故障特征的提取路径

相关分析这个词在信号处理里通常有两层含义:一是统计意义上的相关性计算,比如两个IMF之间的相关系数、IMF与原始信号的互相关函数;二是希尔伯特谱分析所揭示的“频率、时间和能量之间的关联关系”。实际做故障诊断时,我推荐的路径是:

  1. 对原始信号做EMD分解,得到若干IMF;
  2. 分别对关键IMF做包络谱或希尔伯特谱分析,提取对应的特征频率;
  3. 计算各IMF的瞬时能量和样本熵,作为刻画信号复杂度变化的特征;
  4. 用相关系数或互信息筛选出和原始信号高度相关的IMF,丢弃伪分量。

例如滚动轴承出现外圈故障时,IMF中会出现以转频为间隔的边频带,希尔伯特谱上可以看到周期性的能量脉冲。如果只看原始信号的傅里叶谱,这些边频可能被淹没在噪声中;但经过EMD分解后,故障成分集中到了少数几个IMF中,再做包络谱分析,特征频率就非常清晰了。

5. 实战中会被反复坑的三件事

5.1 端点效应:首尾污染为何避不开,又如何缓解

EMD筛分过程依赖三次样条包络,而样条插值在数据序列两端之外没有约束点,上包络线和下包络线在端点附近会剧烈摆动,直接影响第一个IMF和最后一个分量的形态。这种端点效应最典型的表现是:分解结果的头部和尾部出现明显的大幅度扭曲,甚至冒出原本不存在的振荡模式。

缓解方法主要有三种:信号延拓、端点镜像、特征波延拓。最常用的是镜像延拓,把端点附近的波作为镜子,在两端各延拓出一定长度的数据,再做样条插值,分解完成后再把延拓部分裁掉。有些EMD库已经内置了端点处理策略,但默认效果不一定适合你的数据。我通常会在预处理阶段自己用最大值、最小值端点做一次最小振铃延拓,保证两端各保留5到10个极值点。

即便做了延拓,端点效应也不可能完全消除。因此,分析结果时最好避开首尾各10%的时间区段,尤其当你看的是瞬时频率曲线时,一定要忽略前几个点和最后几个点的异常跳变,不要让端点的毛刺干扰你对整体趋势的判断。

5.2 模态混叠:接近频率分量被搅成一团

模态混叠是EMD最经典也很难绕开的缺陷。当原始信号中包含两个频率非常接近的分量,或者高频分量与低频分量幅度差距悬殊时,EMD可能会把本应分开的成分揉进同一个IMF,导致分解结果语义不清。

这个问题的典型触发场景是原始信号中叠加了间断的强干扰,比如振动信号中混入了一段突发的电脉冲。脉冲信号的能量会在很宽的频带内扩散,EMD的筛分过程会把间歇性高频成分和连续高频成分纠缠在一起。解决思路通常是加噪声辅助:EEMD(集合经验模态分解)通过多次添加不同实现的白噪声再做平均,利用噪声的均匀频谱特性把不同尺度“拉开”,最后平均消除噪声影响。CEEMDAN则更进一步,引入了自适应噪声并在每个阶段添加定制的噪声分量。

EEMD和CEEMDAN的实现同样可以在PyEMD库中找到。使用时要注意噪声幅值和集合次数的配合:如果噪声幅值太小,起不到抗混叠作用;太大又会让分解结果失真。一般经验是噪声标准差设为原始信号标准差的0.1到0.3倍,集合次数取100到300次,视数据长度和复杂度而定。

5.3 参数与预处理:采样率、去均值、噪声辅助的取舍

数据处理中常被忽略的一步是去均值。如果信号带有直流分量,EMD会很容易把直流趋势分解为一个单独的IMF,这本身没什么问题,但它会占用一个筛分通道,还可能让后续IMF的幅度估计出现偏差。我习惯先把信号减去均值,再取方差归一化,分解完成后再按实际物理单位还原。

采样率对EMD的影响也需要留意。过低的采样率会让极值点稀疏,包络拟合失真;过高的采样率则会显著增加计算量。EMD的时间复杂度并不低,它是逐级迭代的过程,数据量翻倍带来的耗时增长接近几倍而不是线性关系。对长序列数据,建议先做降采样或分段处理,而不是直接把几十万点的数据一次丢给EMD。

筛分参数上,默认设置通常会表现得不错,但当你发现高频IMF出现很多毛刺、或者低频IMF出现明显分段畸变时,要优先检查两点:一是筛分停止阈值是否过严,导致过度筛选;二是端点处理策略是否与信号形态匹配。在我看来,EMD的参数调试更像是一门经验手艺,针对同一类信号积累几轮调参经验之后,效果会稳定得多。

6. 把EMD接入更大系统:特征工程与工程化约束

6.1 IMF相关性筛选:剔除伪装成分量的噪声

EMD分解得到的IMF并不是每一个都有用。由于端点效应、过筛分和噪声干扰,分解结果中经常出现一些能量虽然不小、但和原始信号物理过程毫无关系的伪分量。用一个简单有效的过滤器:对每个IMF,计算它与原始信号的皮尔逊相关系数,以及它与相邻IMF的相关系数。如果某个IMF与原始信号相关系数低于0.1,且它与其他IMF之间没有明显的频谱重叠,就可以果断丢弃。

还有一种更实用的做法是计算每个IMF的互相关函数峰值和时延。在分析两个通道信号时,可以分别对两个通道做EMD分解,然后按频率范围配对IMF,观察互相关函数的峰值大小和时延方向,判断两个通道之间的耦合关系。这种“先分解、再相关”的思路,在处理机械振动多测点数据和脑电通道间联动分析时特别有效。

6.2 把IMF特征喂给机器学习模型

EMD和机器学习的结合是近年来很火的方向。通常做法是:先对每条样本做EMD,然后分别计算前几个IMF的时域统计特征、频域峰值频率、能量占比、样本熵或排列熵,组成一个完整的特征向量,再送入分类器。

这里有一个值得注意的小技巧:不要直接截取固定数目的IMF,然后对齐特征,因为不同样本分解出的IMF个数可能不一样。我在实际项目中的做法是,统一只取前K个IMF(通常取3到5个)做特征提取,原因是高频IMF携带的个性化特征最明显,而低频趋势项在不同样本之间差异很大,反而会干扰分类模型。如果K个IMF不够,可以用IMF累加能量占比超过90%来确定有效分量数,再把特征向量统一到固定维度。

6.3 数据长度、计算耗时与边缘端部署的现实问题

EMD在离线分析场景很好用,但如果你想把EMD算法搬上实时监测系统,要提前做好心理准备。一次完整的EMD分解耗时随数据长度和筛分次数显著增加,对一个10000点的信号可能耗时几十毫秒到几百毫秒,这还要看硬件性能。实时性要求高的场景,通常只能做分段滑窗分解,而且窗口长度不能太长,可以选择每次只分解最近2048点,并允许一定的端点效应。

还有一个工程化选择:EMD算法本身没有闭式解,每一步都依赖包络拟合和均值迭代,很难直接塞进单片机等资源受限设备。在嵌入式侧,更常见的是先用带通滤波器或小波包做粗筛,再对筛选出的关键子带做EMD分析。这样既利用了EMD的自适应分解能力,又把计算开销控制在了可接受范围内。

我个人在实际使用中的习惯是:面对测量信号,先用傅里叶谱做快速预览,再用EMD处理那些“明显非平稳”的部分。EMD不是万能的,它解决的是“自适应分解”和“时频聚焦”的问题,而不是所有信号分析问题。当你掌握了对IMF做物理解释的能力之后,会发现这套方法在陌生信号上带来的洞察,很多时候是其他频域方法很难替代的。最后建议你把自己常用的信号类型积累成标准测试集,每次调整参数后都跑一遍回归,这样心里才有底。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询