★ FEATURED ARTICLE 最近在整理多通道信号分析流程时又把“多元经验模式分解”这条技术路线完整地走了一遍。说实话这个方向一直被不少做数据处理的朋友忽视觉得用普通傅里叶变换或者一路EMD分解就够了。但当你手里的数据变成了多传感器同步采集、多通道互相耦合的时候单通道方法很快就会撞到天花板。我这次从模拟数据生成、算法实现、参数调到可视化全流程走下来最大的体会是多元经验模式分解不是“把单通道EMD复制几份”它是在高维空间里重新定义了“经验模式”这件事。这篇文章适合三类人看一是做振动、声学、生理或状态监测数据分析的工程师二是刚开始接触经验模式分解想搞明白MEMD和EMD到底差在哪的学生三是已经在用EMD/EEMD但被多通道模态不对齐问题折磨过的人。我会把原理、代码、踩坑和可视化方案都摊开讲尽量让所有内容可以直接复现。1. 先说痛点单通道EMD在多源数据面前为什么不够用1.1 通道独立分解带来的模态错位经验模式分解EMD的核心思想很简单把一个信号自适应地分解成若干个本征模态函数IMF从高频到低频逐层剥离。它不需要预设基函数完全由数据本身驱动。这个特性让它在分析非线性、非平稳信号时比傅里叶变换灵活得多。但问题在于经典EMD是单通道算法一次只能处理一个一维时间序列。多通道数据到来后很多人第一反应是“那我对每个通道分别做EMD不就行了”。我最初也是这么干的直到在跨通道分析时发现了一个让人头大的问题通道A的第2个IMF是50Hz附近的分量通道B的第2个IMF却是15Hz附近的分量两个IMF的阶数对应上了但物理含义完全不同根本没法直接做后续的联合分析。这种现象本质上是EMD的“数据驱动”特性造成的。EMD在分解时只看当前通道自身的极值分布和包络结构不同通道的信号形态稍有差异分解出的IMF带宽和中心频率就会有偏移。结果就是你得到了一堆IMF但它们的“模式层”并不对齐跨通道比较、聚类、相关性分析全都没法做。1.2 共享模式被拆散跨通道融合信息丢失更隐蔽的问题是如果多个通道本来共享同一个物理源的模式单通道EMD很容易把这个共享模式在不同通道中拆进不同的IMF甚至拆散成几段。我举个例子某设备振动监测里三个测点都能感受到由同一转轴引起的振动分量这个分量的幅值、相位在不同测点不同但频率是一致的。对每个通道单独做EMD时这个共享分量在通道1可能完整落在第3个IMF里在通道2却因为局部干扰被硬生生拆成了两个IMF在通道3则可能因为掩膜了趋势项直接混进了残差。这样一来你在做多通道联合分析时根本没有办法把“属于同一个源的模式”抽出来。想用多通道信息做故障定位或者源分离第一步就卡住了。1.3 独立分解的复现性与一致性难题还有一个经常被忽视的问题对同一组数据如果你把通道顺序换一下或者只截取其中一部分数据重新分解各通道IMF的编号和频率分布可能都会变。单通道EMD本身就对端点效应和极值分布敏感独立分解等于把这种不稳定性放大到了整个通道集合。我在一个多通道实验数据上验证过同样的通道1数据单独做EMD和放在一组三通道数据里独立做EMD前3个IMF的中心频率差异能达到15%左右。如果你需要用分解结果做阈值判断或者搭建自动化流程这种不稳定性会让你不得不反复调参难以交付可靠的处理链路。这些问题综合下来指向一个核心需求我们需要一种能把整个多通道数据当成一个整体来分解的方法让分解出的每个IMF天然就带有“多通道”的维度所有通道在同一模式下对齐。这也是多元经验模式分解存在的根本理由。2. 多元经验模式分解的运作机制从方向向量到多元包络2.1 把多通道数据看成一个高维空间里的曲线要理解MEMD先做一个视角转换。单通道信号是一维时间序列它可以被看成一条在时间轴上起伏的曲线。如果是三通道同步采样那么每个时刻的三个值可以组成一个三维空间里的点整段数据就是一条在三维空间里随时间延伸的曲线。推广到C个通道就是一条在C维空间里的参数曲线。EMD里的“包络”概念在一维情况下是通过连接局部极大值和极小值得到的。到了C维空间我们需要构造的是一个“多元包络”它不再是上下两条线而是在C维空间里围绕数据曲线的一个包络面。这个包络面怎么构造MEMD给出的方案是先从多个方向观察这条曲线再把观察结果综合起来。2.2 方向向量采样、投影与多元包络构建MEMD的完整循环可以分为这样几个步骤。第一步在C维单位球面上均匀采样一组方向向量。这些方向向量可以理解为不同角度的“观察方向”数量一般取16到64个不等具体和通道数有关。第二步将当前残余信号沿着每个方向向量做投影得到一条一维投影序列。这条投影序列相当于“从某个方向看过去数据曲线的侧面轮廓”。第三步对每条投影序列找局部极值点。这些极值点代表的是数据曲线在该方向投影上的“凸起”和“凹陷”位置。第四步把不同方向投影得到的极值点全部映射回原始C维空间在时间轴上做多维插值。常用的是径向基函数插值也有用薄板样条的做法。插值完成后会得到一组平滑的包络曲面对每个通道而言就有了一条包络线。第五步计算所有通道包络的均值用原始信号减去这个包络均值完成一次筛分迭代。这个减除过程和一维EMD类似是逐次逼近本征模态函数的过程。这里最关键的区别在于包络均值不是某个通道单独算出来的而是所有方向投影、所有通道的极值信息共同作用的结果。也就是说每个通道的IMF迭代过程都“知道”其他通道发生了什么。这正是MEMD能实现多通道模式对齐的根本原因。2.3 停止准则与残差的递归分解筛分迭代不能无限进行下去。如果一直减下去最终每个IMF都会变成一个纯调频信号失去物理意义。MEMD里常用的停止准则有两类一类是包络均值的相对标准差小于某个阈值另一类是限制筛分迭代的最大次数。我在工程实践中倾向于同时用两个条件先到先停这样既能保证IMF的稳定性又不会过筛。第一轮迭代收敛后得到的是第一个多元IMF它包含所有通道中最高频的公共模式成分。然后把这个IMF从原始数据中减去对残余信号重新执行上述过程得到第二个IMF。如此递归下去直到残余信号不再含有可分解的模式为止。最终你会得到一组多通道IMF它们的每个IMF都包含所有通道在同一时刻的数据切片。2.4 一个生活化类比多人观察同一物体如果觉得上面的过程不太好想象可以用一个类比。假设你面前有一个形状不规则的物体你想画出它的轮廓。一个人从正面看会漏掉侧面细节一个人从侧面看也会漏掉背面信息。但如果你找一群人站在不同方向同时观察把每个人的轮廓信息汇总起来再综合描出物体的完整轮廓这个过程就非常接近MEMD的多元包络思想。不同之处在于MEMD是动态的物体形状随时间变化观察者们每时每刻都在记录轮廓变化最后生成的不是一张静态轮廓图而是一组随时间演化的模式分量。正因为有了多方向观察任何单一方向被噪声干扰导致漏检的极值都可能被其他方向的投影信息补回来这也让MEMD在抗噪和模式分离上比单通道EMD更稳。3. 动手前的数据准备模拟多通道信号与工程采集注意事项3.1 用已知成分的模拟数据验证算法开始写分解代码之前我强烈建议先构造一组“答案已知”的模拟数据。否则你根本不知道算法输出的IMF是对是错。我这次用的是三通道数据采样率1000Hz时长4秒。三个通道共用一个10Hz的正弦分量用来模拟“共享源”。通道1额外加了一个60Hz的高频正弦通道2加了一个35Hz的方波通道3加了一个缓慢漂移的趋势项最后每个通道都叠加白噪声。这样设计的目的很明确我希望分解出的某个IMF能把10Hz公共分量完整提出来另外几个IMF分别对应各通道的独立特征。实际生成数据的代码大致是这样import numpy as np fs 1000 t np.linspace(0, 4, fs * 4, endpointFalse) common 0.8 * np.sin(2 * np.pi * 10 * t) comp_ch1 0.5 * np.sin(2 * np.pi * 60 * t) comp_ch2 0.4 * np.sign(np.sin(2 * np.pi * 35 * t)) comp_ch3 np.linspace(0.6, -0.4, len(t)) rng np.random.default_rng(42) noise lambda: 0.15 * rng.standard_normal(len(t)) x1 common comp_ch1 noise() x2 common comp_ch2 noise() x3 common comp_ch3 noise() X np.stack([x1, x2, x3], axis0)这样得到的X形状是3行、4000列行是通道列是时间点。后面所有MEMD流程都基于这种[通道数, 采样点数]的矩阵布局。模拟数据的好处是分解完你能直接对照IMF和真值知道哪些分量被正确分离、哪些地方出现了模式混叠这是真实采集数据给不了的验证条件。3.2 真实采集数据的三条红线模拟数据做好之后再说说真实采集数据的要求。第一所有通道必须是同一时钟采样长度必须一致。如果各通道分别采集再拼接时间基准对不齐MEMD分解出的公共模式毫无意义。处理前先检查时间戳必要时做重采样对齐。第二采集前做好抗混叠滤波。采样率至少要高于你关心的最高频率成分的2.56倍这是工程监测里比较稳妥的经验值不是理论下限。否则高频成分混叠进低频带会让MEMD把虚假模式也当成真实IMF分解出来而且很难在事后补救。第三数据进入分解前先做中心化。MEMD本身不要求零均值但如果通道之间存在明显直流偏置差异会导致方向向量投影结果失真。我习惯先把每个通道减去自身均值再做分解。如果趋势项太强我会先用高通滤波或者多项式去除大趋势只把剩余部分交给MEMD处理。4. Python落地从零编写的MEMD核心流程与实测4.1 现有库的现实情况别指望拿来即用Python生态里单通道EMD的实现不少但完整支持多通道MEMD的开源库其实很稀缺。很多库提供的是EMD、EEMD甚至CEEMDAN但你要找“输入一个多通道矩阵输出一组对齐IMF”的现成接口大概率要失望。我的做法是自己实现核心流程。一方面是为了不受库的限制另一方面是只有把筛分循环看懂你才能真正理解方向向量数量、停止阈值这些参数在调什么。如果只是调库出了问题根本无从下手。4.2 核心筛分循环方向向量、投影极值、包络均值下面是我写的简化版MEMD核心循环。严格来说完整的MEMD要使用多维插值来构造高维包络代码会非常长。我这里用“逐个方向投影做包络再取平均”的近似方式兼顾可读性和实用性工程上已经能获得不错的效果。from scipy.interpolate import CubicSpline def sample_direction_vectors(n_channels, n_dirs, seed0): rng np.random.default_rng(seed) v rng.normal(size(n_dirs, n_channels)) return v / np.linalg.norm(v, axis1, keepdimsTrue) def find_local_extrema(s): d np.sign(np.diff(s)) max_idx np.where((d[:-1] 0) (d[1:] 0))[0] 1 min_idx np.where((d[:-1] 0) (d[1:] 0))[0] 1 return max_idx, min_idx def envelope_mean_memd(H, vectors): n_ch, n_len H.shape env_sum np.zeros((n_ch, n_len)) for v in vectors: proj H.T v t np.arange(n_len) for sign in (1, -1): idx, _ find_local_extrema(sign * proj) if len(idx) 2: continue cs CubicSpline(t[idx], (sign * proj)[idx]) envelope sign * cs(t) env_sum envelope return env_sum / (2 * len(vectors)) def memd_sift(X, n_dirs16, max_sift200, sd_threshold0.01): n_ch, n_len X.shape vectors sample_direction_vectors(n_ch, n_dirs) residue X.copy() imfs [] while True: h residue.copy() for _ in range(max_sift): m envelope_mean_memd(h, vectors) h h - m sd np.sum((m ** 2)) / np.sum((h ** 2)) if sd sd_threshold: break imfs.append(h) residue residue - h n_extrema 0 for v in vectors[:4]: proj (residue.T v) n_extrema len(find_local_extrema(proj)[0]) if n_extrema 2: break imfs np.stack(imfs, axis0) # shape [n_imfs, n_ch, n_len] return imfs这段代码的核心在envelope_mean_memd函数里对每个方向向量先把信号投影成一维序列找极值再做三次样条插值得到包络所有方向的包络汇总后取平均作为多元包络均值的近似。sd的计算是包络均值的能量占当前信号能量的比例小于阈值就认为筛分收敛。用前面生成的模拟数据跑一遍我在实测里得到了比较理想的结果10Hz公共分量基本集中到了某个中频IMF上60Hz正弦和35Hz方波分别被拆进了不同的高频IMF趋势项则残留在最后的残余分量里。虽然35Hz方波因为自身谐波丰富被拆成了相邻两个IMF但整体模式分离效果远超三个通道分别做EMD的结果。4.3 为什么先投影、再极值、再包络很多人刚接触时有一个疑问为什么不能直接在原始多通道数据的每个通道上找极值、做包络那样不是更直接吗原因在于如果你只在单个通道上找极值和做包络你的包络信息和其它通道没有任何联系又回到了单通道EMD的老路。MEMD的精髓是“同时观察所有通道在某个方向上的整体变化”。沿着方向向量投影后每个时刻的投影值实际是所有通道的加权组合极值点反映的是这个组合形态的转折。不同方向的投影极值综合起来才能刻画出高维数据曲线在不同方向上的起伏最后构造出的包络均值才带有全局信息。这也解释了为什么方向向量的数量和质量很重要。方向太少等于只有几个人从少数角度观察物体轮廓自然不完整方向太多计算量剧增而且相邻方向提供的信息高度冗余边际收益递减。5. 分解结果的四层可视化时域波形、筛分过程、时频谱与相关性5.1 多元IMF波形排列图先看大体分解效果拿到一组IMF后我第一步永远是画“通道×IMF”的排列图。横轴是时间纵轴把原始信号和各个IMF分开展示。三通道数据会得到3个原始信号、若干组3通道IMF还有最后的残差。画这种图不需要什么特殊技巧用matplotlib的subplots就能搞定。重点看两点一是每个IMF的频率是否由高到低排列二是同阶IMF下三个通道的波形是否在时间上对齐。如果两个相邻IMF在某个时间段都出现了明显的相同振荡说明存在模式混叠如果某个通道的IMF突然出现大面积端点大幅摆动说明端点效应已经影响到了这个分量。从这张图里你能直观判断出分解是否值得继续往下做而不是直接跳到定量分析。5.2 筛分过程可视化把迭代收敛过程摆出来这是很多人会跳过的一步但我认为非常值得做。在筛分循环里每一次迭代都会计算出一个包络均值你可以把这些包络均值按迭代次数堆叠成一张图看看它是如何逐次逼近零的。实际操作时我在memd_sift函数里加了一个history参数每轮迭代把包络均值存下来。然后在分解完成后把某一通道的原始信号、第1次迭代的包络均值、第5次、第20次迭代的包络均值画在一起。观察结果是前几次迭代包络均值幅度很大主要是在剥离最高频波动后面迭代的包络均值幅度迅速衰减说明高频成分已经被充分剥离再继续迭代容易把真实模式的能量也减掉。这张图能帮你判断停止准则是否合理。如果包络均值衰减得很慢说明方向向量数量可能不够如果衰减得太快说明筛分阈值太松IMF还混着不少噪声成分。5.3 希尔伯特时频谱与边际谱分解只是第一步可视化才是让结果“说话”的关键环节。对每组IMF做希尔伯特变换可以得到瞬时幅度和瞬时频率然后把它们放在时间–频率平面上就得到希尔伯特时频谱。代码思路如下from scipy.signal import hilbert def hilbert_spectrum(imf_group, fs): analytic hilbert(imf_group, axis-1) inst_amp np.abs(analytic) phase np.unwrap(np.angle(analytic)) inst_freq np.diff(phase, axis-1) * fs / (2 * np.pi) return inst_amp[:, :, :-1], inst_freq画图时用pcolormesh横轴时间纵轴频率颜色代表瞬时幅度。需要注意的是瞬时频率对噪声非常敏感我建议先对IMF做轻微的平滑滤波或者只显示幅度高于某个阈值的点否则图上会被噪点铺满看不出结构。边际谱则是对瞬时幅度在时间上求和得到的包络谱它能反映每个频率成分在整个时间区间上的能量贡献。和传统傅里叶幅值谱相比希尔伯特边际谱更擅长处理频率随时间变化的非平稳信号。我是把两种谱叠在一起对比看的如果两类谱在主要频峰上一致说明分解结果是可信的如果出现明显差异多半是某个IMF出现了过分解或欠分解。5.4 通道相关性热力图检验模式是否对齐最后一个可视化是相关性热力图。对每一阶IMF计算三个通道之间的相关系数得到一个相关系数矩阵。这个矩阵能直接回答前面提到的“模式对齐”问题。具体方法对第k阶IMF取它在三通道上的数据切片两两计算Pearson相关系数然后画成heatmap。我实测模拟数据时的结果很有意思公共10Hz分量所在的IMF三通道相关系数都在0.98以上而包含各通道独立分量的IMF相关系数明显偏低。这个对比说明MEMD确实把共享模式集中到了同一个IMF里跨通道相关性被显著强化了。实际工程中这套逻辑可以进一步发展成“模式指纹”每个IMF的三通道相关系数矩阵、中心频率、能量占比组合起来就能作为判断设备状态、识别故障类型的特征向量。6. 实战复盘多通道振动数据中冲击特征的提取6.1 场景与数据设计下面用一个虚构但贴近实际的状态监测案例来走一遍完整流程。场景是某旋转设备实验台三个加速度传感器分别布置在三个测点采样率25600Hz采集时长0.5秒。设备转频约25Hz故障特征表现为每转一次的周期性冲击冲击激发设备结构产生约1230Hz的高频衰减振荡。三通道都能感知到冲击但冲击的幅值和相位略有差异且叠加了较强的背景噪声和转频谐波。这种数据的特点是单通道里冲击特征混在高频谐波和噪声里直接看原始波形很难辨认用固定带通滤波器又需要知道共振频带而共振频带往往随工况漂移。MEMD的价值在这里体现得很明显它不需要预设频带完全根据数据自适应地把冲击模式和其他模式拆开。6.2 分解流程与可视化按前面代码跑分解三通道数据会得到6到8个IMF。在我的实测中高频冲击的衰减振荡集中在前1到2个IMF里转频和谐波落在第3到第5个IMF最低频的残余分量代表的是整体趋势和偏置。关键操作是对包含冲击的IMF做希尔伯特时频谱。在时频谱上你能看到一系列等时间间隔出现的竖直能量条纹间隔刚好对应转频周期25Hz。这些竖条纹就是冲击在时频图上的足迹比在原始时域波形里看要清晰得多。接着对该IMF的包络信号做一次幅值谱就能得到调制包络谱。冲击的重复频率正好是转频在包络谱里会看到25Hz及其整数倍处的明显谱峰。如果再对包络谱做包络分析还可以提取出冲击的阻尼衰减特征用于进一步估计故障严重程度。6.3 结果解读的注意事项我必须强调一点MEMD在这里是特征提取工具不是判决工具。它帮助你干净地分离出冲击相关模式但具体这个冲击代表什么故障还需要结合设备结构、转速、历史数据综合判断。不同实验台、不同传感器布局分解出的IMF阶数和内容都会有差异不要照搬别人的结论。还有一个容易被忽略的细节如果你的冲击很微弱在MEMD分解后可能被拆到多个相邻IMF里而不是集中在一个IMF。这时不要急着把单个IMF丢掉可以尝试把相邻的两个高频IMF叠加起来再做包络分析有时能把微弱冲击重新“拼”出来。这算是实操中一个比较实用的小技巧。7. 参数调优与踩坑记录方向向量数量、端点效应与模态混叠7.1 方向向量数量经验范围与选择逻辑方向向量数量是MEMD最重要的参数也是最容易让人困惑的参数。数量太少包络均值不够平滑分解出的IMF会抖动数量太多计算量成倍增加但包络信息逐渐冗余效果提升有限。我的经验公式是对于C通道数据方向向量数量在3C到5C之间比较合适。三通道数据用10到20个方向已经足够我通常直接取16。如果分解结果在重复运行时不稳定优先增加方向向量数量而不是去改筛分停止阈值。通道数达到八个以上时我会把方向数增加到32或48同时做好接受较长运行时间的心理准备。7.2 停止准则过筛与欠筛的平衡筛分停止准则直接影响IMF的纯度。阈值设得太严筛分迭代次数太多会把信号的真实能量一点点“刮”走阈值设得太松IMF里会残留大量相邻频带成分导致模态混叠。我常用的组合是包络均值相对标准差阈值取0.01最大筛分迭代次数200次。如果发现某个IMF明显呈现出“噪声包络”的样子即傅里叶谱非常平坦说明它被过筛了应降低最大迭代次数。如果相邻IMF的中心频率离得很近疑似欠筛则把阈值改到0.001并提高最大迭代次数到400。7.3 端点效应短序列的噩梦与缓解办法EMD族算法都有端点效应问题MEMD也逃不掉。数据两端缺乏极值包络会在端点处发散导致IMF两端出现大幅摆动的假振荡。数据越长端点效应影响比例越小数据只有几百个采样点时整个两端都可能被污染。我的处理办法是镜像延拓分解前把每段数据两端各扩展一段镜像序列分解完成后裁掉。具体思路是取数据末尾一小段翻转后接在尾部开头同理。这样做能在不引入过多虚假频率的前提下缓解端点发散。如果数据实在太短另一个务实做法是放弃两端的分析结果只取中间稳定区段。7.4 模态混叠当相邻IMF开始“串门”模态混叠指的是本应属于一个IMF的多种频率成分被拆散到不同IMF里或者不同时间段的相同频率成分在两个IMF之间交替出现。这种现象在MEMD中依然存在只是比单通道EMD轻一些。我的实测经验是噪声幅值越大模态混叠越明显。缓解办法有两类。一类是在分解前对数据进行预处理把明显的高频噪声带滤掉再分解。另一类是引入噪声辅助的思想向多通道数据中加入低幅值的白噪声重复多次分解再把对应阶次的IMF做平均。这个过程在文献中称为噪声辅助多元经验模式分解思路和单通道EEMD一脉相承。需要注意的是注入噪声的幅值不要超过原始信号标准差的0.2倍否则会引入虚假模式。7.5 算力不够怎么办MEMD的计算量远大于单通道EMD。通道数多、方向向量多、数据长、筛分迭代上限高这几个因素叠加起来跑一组数据等半天的情况完全可能出现。我在处理高通道数数据时有两个习惯一是对数据做PCA降维后再分解保留主要能量对应的前几个主成分能大幅降低MEMD的通道维度二是把分解任务拆到多进程里跑各个通道组互相独立并行化收益很明显。写这段的时候我特意回忆了一下自己在这条路上踩过的坑。最初用现成库时我根本不理解输入数据要不要归一化结果方向向量的采样严重受高幅值通道支配分解结果乱成一团。后来认真看了几次筛分循环的输出才明白包络均值是跨通道加权的不把通道幅值拉到一个量级高幅值通道会“带偏”所有方向的投影极值。现在我的流程里永远包含一步标准化到零均值、单位标准差分解完再乘回原始幅值。这个细节不算高深但非常影响最终结果的质量。分解结束后我还会做一次“账目核对”把全部IMF和残差逐点相加重建原始信号看重建误差的数量级。如果重建误差在1e-10以下说明分解过程没有丢信息业务逻辑上可以放心使用。这个方法简单有效建议你每次跑完MEMD都顺手做一次比我见过的任何“质量指标”都更让人安心。
阅读完成 · 觉得有帮助？
