简介基于变分模态分解VMD、排列熵与极限学习机ELM的滚动轴承故障诊断Python实现面向机械设备健康管理研究人员、工业运维人员及相关专业学生着力解决振动信号特征提取与故障模式识别的实际需求可用于教学实验和工程验证。资源包共407个文件压缩后约4.8MB其中404个txt为不同工况下的轴承振动信号数据2个Python脚本分别覆盖VMD分解与排列熵计算、ELM模型训练与分类另有1个csv文件存放特征或标签整体目录清晰便于定位和复用。目前已有633人学习下载。代码结构简单、可直接运行完整演示了从原始信号到VMD模态分解、排列熵特征提取、ELM分类诊断的端到端流程并给出了数据准备、算法调用与结果输出的具体实现读者可在此基础上更换数据集、调整参数将方法迁移到其他旋转机械的故障诊断场景中。1. 从振动信号到故障标签VMD、排列熵和ELM这条链路为什么值得跑通拿到一段轴承振动加速度信号直接做FFT往往只能看到一片宽带噪声故障特征频率被转频和随机冲击淹没。更反直觉的一点是真正有助于分类的信息不一定在频谱峰值里而在信号复杂度的细微变化中。这个项目把VMD变分模态分解、排列熵和ELM极限学习机串成一条完整的Python故障诊断流水线先对振动信号做自适应分解把非平稳信号摊成多个窄带模态再用排列熵把每个模态压缩成一个复杂度标量最后用ELM做多分类。整个流程对CWRU这类公开轴承数据集非常友好从读入txt到输出准确率只需要几十秒。它适合设备健康管理方向的算法工程师做基线方案也适合研究生拿来做故障诊断课程的完整课设下面按实际运行顺序把每个环节拆开讲。2. VMD模态分解带宽约束优化与调参落地的第一步2.1 为什么批处理场景下VMD优于EMD经验模态分解EMD的递归筛选方式存在两个老问题一是筛分过程中极值点包络的插值误差会累积导致模态混叠同一个物理频率成分被拆到多个IMF里二是端点处的样条拟合容易发散产生虚假的振荡成分。VMD换了个思路不再筛而是把一个约束变分问题的求解作为分解手段假设原始信号可以表示为K个带限本征模态函数之和每个模态围绕各自的中心频率分布优化目标是最小化所有模态的估计带宽之和同时要求模态总和等于原始信号。求解过程中引入二次惩罚因子alpha和拉格朗日乘子用交替方向乘子法迭代更新模态和中心频率直到收敛判据满足。这个框架带来的直接好处是VMD对采样率变化和噪声更稳健分解结果对模态数K有明确依赖且迭代收敛比EMD的筛分循环更容易控制。它在故障诊断中的典型用法是把高频冲击成分和低频调制成分分开让后续特征提取只看故障相关的模态。代价是K需要预先给定alpha也影响模态的带宽。这两个参数在2.3节里单独展开。在项目对应的Python实现中核心分解功能由vmdpy库提供特征提取脚本只负责把原始振动数据送入VMD再对返回的模态矩阵做排列熵计算。2.2 vmd-pailieshang.py的信号加载与分解调用项目里的vmd-pailieshang.py承担特征提取入口的角色。它的输入是原始振动信号的txt文件输出是排列熵特征表pailieshao.csv。原始数据文件如inner0.txt、inner60.txt每行一个采样点的加速度幅值属于典型的时序信号存储格式。数据量较大时通常截取一段分析窗再分解避免整段信号直接参与VMD导致计算时间不可控。import numpy as np from vmdpy import VMD # 读取CWRU内圈故障振动信号只取前2048个采样点 signal np.loadtxt(inner60.txt, dtypenp.float64)[:2048] # VMD关键参数 alpha 2000 # 惩罚因子控制模态带宽 tau 0 # 噪声容忍度0表示严格保真 K 4 # 模态分解个数 DC 0 # 0表示不把第一个模态当作直流分量 init 1 # 中心频率初始化方式1为均匀初始化 tol 1e-7 # 迭代收敛阈值 u, u_hat, omega VMD(signal, alpha, tau, K, DC, init, tol) print(模态矩阵形状:, u.shape) print(各模态最终中心频率:, omega[-1, :])这段代码执行后u的形状是(K, N)每一行是一个带限模态分量omega[-1, :]保存的是最后一次迭代收敛时各模态的中心频率。中心频率的用途不只是观察它还能帮你判断分解是否合理比如两个模态的中心频率挤在一起说明K设大了。alpha取2000是工程中较稳妥的起点它让低频调制成分和高频冲击成分在能量上保持相对均衡tau取0意味着不使用噪声容忍项适合信噪比尚可的实验室公开数据集。DC置0是因为轴承振动信号不存在需要单独分离的直流漂移。tol控制迭代停止精度1e-7已经足够。2.3 模态数K、惩罚因子alpha和数据窗长的联合调整K和alpha不是独立起作用的。alpha固定为2000时K从2增加到8会出现三个阶段K过小时不同频率成分挤压在同一个模态里排列熵偏低且彼此接近K适当时各模态中心频率分散故障冲击集中在某一两个模态上K过大时会产生虚假模态把正常噪声也拆成独立的分量。可以打印最终中心频率来辅助判断例如K4时如果得到一组类似125、580、2100、3800 Hz的中心频率分布且相邻间隔没有重叠一般可以认为分解合理。reconstruct np.sum(u, axis0) rmse np.sqrt(np.mean((signal - reconstruct) ** 2)) / np.std(signal) print(fK{K}, alpha{alpha}, 归一化重组误差{rmse:.6f})重组误差是另一个直接指标理论上分解加和等于原信号误差应趋近于0。工程上我一般允许误差在1e-3量级如果误差偏大优先检查alpha是否设置过小。下面是本地跑的一组参考结果不同数据源会有波动但趋势可复用。K值中心频率分布Hz归一化重组误差结论倾向2320, 24502.1e-4模态数不足低频段信息被合并4130, 620, 1900, 36004.8e-4分布合理推荐695, 420, 1200, 2400, 4100, 68003.6e-3高频出现虚假模态数据窗长方面2048点是常用选择。窗太长会摊平故障冲击的局部性窗太短则VMD对频率分辨不足模态中心频率不稳定。做批量样本时把同一段信号按50%重叠切窗能让样本数量成倍增加且不引入额外采集成本。3. 排列熵特征提取从模态波形到特征向量的压缩3.1 序模式统计如何量化信号复杂度排列熵是一种不依赖信号幅值绝对大小的复杂度指标。它把时间序列转换成一组符号序列对每个位置取连续m个点按数值大小排列得到一个顺序索引比如三个点的排序结果是(1,2,0)表示第二个点值最小、第三个点其次、第一个点最大。这些排列模式出现的概率分布就是信号局部序结构的统计描述均匀随机信号中所有排列模式趋近等概率熵值高规则周期信号中只有少数几个排列模式反复出现熵值低。轴承故障对振动信号的影响会改变这种序结构早期故障产生周期性冲击使信号在局部时段内更有规律排列熵下降严重故障引入大量非线性调制排列熵回升。因此排列熵并不直接对应某种故障频率而是作为区分健康状态的特征量进入分类器。计算时两个参数很关键嵌入维度m和时延t。m过小只有m!种模式区分度不够m过大会出现很多未充分采样的模式统计噪声上升。振动信号通常取m4到5t取1或2。3.2 排列熵函数实现要点import numpy as np from itertools import permutations def permutation_entropy(signal, order3, delay1, normalizeTrue): 计算一维信号的排列熵 order: 嵌入维度对应排列长度 delay: 采样间隔默认1表示逐点比较 normalize: 是否除以 log(order!) 归一化 n len(signal) patterns list(permutations(range(order))) pattern_to_idx {p: i for i, p in enumerate(patterns)} cnt np.zeros(len(patterns)) for i in range(n - (order - 1) * delay): window signal[i:i order * delay:delay] order_idx tuple(np.argsort(window)) cnt[pattern_to_idx[order_idx]] 1 prob cnt / cnt.sum() prob prob[prob 0] pe -np.sum(prob * np.log(prob)) if normalize: pe / np.log(len(patterns)) return pe这个实现里有三个容易踩坑的细节。第一np.argsort默认返回索引而索引元组正是排列模式的编码方式不需要额外排序比较第二窗口点数为order * delay但排序只针对order个点它比直接取连续order个点能捕捉更大时间跨度的结构第三概率向量里可能包含零计算Shannon熵之前必须剔除否则log(0)会产生警告并把结果变成nan。归一化除以log(order!)把输出压到[0,1]区间不同维度下的特征值才有可比性。这段代码在2048点数据上运行耗时极短适合放进批量特征提取循环。3.3 特征矩阵组装与csv落地特征工程不能只算一个模态要把VMD分解出的每个模态都计算排列熵形成一个与模态数等长的特征向量。对第i个样本做K个模态分解得到的特征行是[PE_1, PE_2, ......, PE_K]。在vmd-pailieshang.py中把内圈故障不同损伤程度的inner0.txt、inner3.txt、inner25.txt等文件依次读取每段信号经过分解、熵计算后追加到特征矩阵最后在矩阵末尾拼上一个整数标签列写入pailieshao.csv。文件来源PE_IMF1PE_IMF2PE_IMF3PE_IMF4标签inner0.txt0.8120.9030.4710.5580inner25.txt0.7840.8870.4360.6021inner60.txt0.7590.9120.4020.5882实际组装时要注意样本与标签的对应关系。一种常见做法是把同一个txt文件先做分段每2048点一段得到N个样本后再统一打标签而不是对整个大文件只取一个样本。这样能显著增加训练数据量ELM才有足够的样本学习类间差异。pailieshao.csv的格式本质上就是特征列最后一列标签的二维表格用np.savetxt或pandas.to_csv都能落盘后续ELM.py直接读取该csv即可。4. ELM分类模型隐藏层随机映射与快速最小二乘训练4.1 单隐藏层前馈网络的随机映射原理极限学习机与BP网络的关键区别在于隐藏层参数不需要反向传播。输入权重W和偏置b在训练前随机生成一旦固定就不再更新隐藏层输出矩阵H只由输入X和固定的(W,b)决定输出层权重beta通过解一个线性最小二乘问题获得。这样的设计让训练过程退化为矩阵运算速度比迭代式梯度下降快一个数量级以上同时避免了BP容易陷入局部极小点的问题。在轴承故障诊断这种典型的小规模分类任务中样本量通常在几百到几千之间特征维度等于VMD模态数ELM的结构优势非常明显。它不需要漫长的调参训练隐藏层节点数在几十到一两百之间就能达到足够好的分类效果。ELM也有代价隐藏层随机映射意味着同样的数据、同样的节点数两次独立运行的结果可能不同因此实际工程中需要通过固定随机种子或适当增大节点数来稳定输出。4.2 ELM.py核心代码与C正则项import numpy as np class ELM: 极限学习机分类器输出层用带L2正则的最小二乘求解 def __init__(self, n_hidden80, C0.01, random_state42): self.n_hidden n_hidden # 隐藏层神经元个数 self.C C # L2正则化系数 self.rng np.random.RandomState(random_state) self.W None # 输入层到隐藏层权重 self.b None # 隐藏层偏置 self.beta None # 隐藏层到输出层权重 def fit(self, X, y): n_samples X.shape[0] # 随机初始化输入权重和偏置 self.W self.rng.uniform(-1, 1, (X.shape[1], self.n_hidden)) self.b self.rng.uniform(-1, 1, (1, self.n_hidden)) # 隐藏层输出矩阵激活函数用双曲正切 H np.tanh(np.dot(X, self.W) self.b) # 加入正则项的最小二乘(H^T H I/C)^{-1} H^T y beta np.linalg.solve(H.T H np.eye(self.n_hidden) / self.C, H.T y) self.beta beta return self def predict(self, X): H np.tanh(np.dot(X, self.W) self.b) return np.argmax(H self.beta, axis1)代码逻辑分三步。第一步用均匀分布在[-1,1]内随机生成输入权重和偏置所有随机性都受random_state控制重复运行时结果可复现。第二步计算隐藏层输出矩阵H激活函数选择双曲正切tanh是因为它对输入的正负号敏感比sigmoid更容易表达振动特征中的正负冲击差异。第三步用np.linalg.solve直接解带正则项的线性方程组。加入I/C的意义在于当隐藏层节点数大于样本数时H^T H可能奇异正则项让它保持可逆同时抑制过大的beta值。C越小正则越强分类边界越平滑。4.3 样本划分、混淆矩阵与评估指标在训练前需要对pailieshao.csv做两件事把特征列和标签列拆开并把原始整数标签转成one-hot矩阵。ELM的输出层节点数等于故障类别数预测时取输出向量中最大值的下标作为类别。以下是一个完整训练与验证的流程片段。import pandas as pd import numpy as np from sklearn.model_selection import train_test_split from sklearn.preprocessing import OneHotEncoder # ELM.py 主流程 data pd.read_csv(pailieshao.csv, headerNone) X data.iloc[:, :-1].values # 前K列为排列熵特征 y data.iloc[:, -1].values.reshape(-1, 1) # 最后一列为标签 enc OneHotEncoder(sparse_outputFalse) y_onehot enc.fit_transform(y) X_train, X_test, y_train, y_test train_test_split( X, y_onehot, test_size0.3, random_state42, stratifyy ) model ELM(n_hidden100, C0.01, random_state42) model.fit(X_train, y_train) y_pred model.predict(X_test) y_test_label np.argmax(y_test, axis1)stratifyy很重要它能保证划分后训练集和测试集中各类别比例与原始数据一致避免某一类故障样本全部集中到训练集、测试集缺失这一类。ELM的隐藏层节点数不是越大越好节点数过少拟合不足节点数过多时在几十个样本的小数据集上会产生轻微过拟合。经验上K4维特征、几百个样本时从50开始逐步增大到200观察测试集准确率的变化选择准确率不再提升的节点数作为最终配置。评估时除了总体准确率建议看一眼混淆矩阵。真实\预测正常内圈轻故障内圈重故障正常3200内圈轻故障1292内圈重故障0133如果错分集中在相邻损伤等级之间说明特征在当前K和alpha配置下区分度不够优先回到VMD调参而不是继续改ELM。如果错分均匀分布则考虑增加样本数或降低C值增强正则。5. 用inner系列文件做验证时容易被忽略的细节最后一章的实战重点放在可复现性上。ELM的随机初始化带来一个隐患即使数据完全不变两次运行也可能得到不同准确率。项目里inner0.txt到inner62.txt多个文件提供的是不同工况下的振动信号验证时需要注意三个问题。第一样本打乱顺序。实验脚本如果直接按文件顺序组装特征矩阵前80%都是同一类故障测试集又集中于另一类故障准确率会虚高。把样本ID或特征行先整体随机打乱再做划分能让评估结果更有意义。可以参考上面代码中train_test_split的用法加上random_state保证每次运行都保持同样的划分结果。第二固定随机种子。在创建ELM实例后如果要做对比实验建议固定random_state否则无法判断准确率变化来自参数调整还是随机波动。更进一步可以用多次运行取平均的做法把同一训练过程跑10次每次用不同的random_state记录准确率均值与标准差标准差低于2%才认为模型稳定。accs [] for seed in range(10): model ELM(n_hidden80, C0.01, random_stateseed) model.fit(X_train, y_train) accs.append((model.predict(X_test) y_test_label).mean()) print(f准确率均值 {np.mean(accs):.4f}标准差 {np.std(accs):.4f})第三样本量不足时优先调整分段策略。如果每个inner文件只取一段2048点信号得到的样本数少于50个ELM很难区分四个故障等级之间的细微差异。把每个文件滑窗切段步长取512点一个文件能生成十几个样本整体训练集就达到百级别。切段后排列熵特征的均值与方差可能变化此时需要检查pailieshao.csv中同一标签下特征值是否相对集中如果某一类特征方差过大说明窗长或K不合适。最后提醒一点故障诊断的终点不是训练集准确率而是留出集上对未知工况的泛化表现验证脚本应保留独立的测试集不做任何手工筛选。本文还有配套的精品资源点击获取