简介这套基于HHT-ELM的故障诊断项目实例面向具备信号处理与机器学习基础、熟悉MATLAB的科研人员与工程师聚焦旋转机械早期故障识别与智能运维场景。项目以希尔伯特–黄变换HHT为核心通过EMD自适应分解提取IMF分量结合希尔伯特变换构造瞬时频率与幅值时频特征再交由极限学习机ELM完成高效分类建模有效应对非平稳强噪声信号下的故障辨识难题。文档从信号采集、预处理、特征提取、模型训练、参数调优到GUI设计部署均有完整阐述并配有可直接运行的MATLAB代码、混淆矩阵等可视化方案及可解释性分析。资源包包含1个docx文档大小109KB内容预览显示其目录涵盖项目背景、挑战与解决方案、模型架构、代码示例等模块方便读者按图索骥。目前已有48人学习适合用于工程验证、科研对比及教学演示。1. 从“看到波形”到“认出故障”为什么 HHT 比传统频谱更先发现问题旋转机械的振动信号里轴承点蚀、齿轮裂纹这类早期故障往往只表现为短暂的冲击脉冲和局部调制常规 FFT 频谱会把这种瞬态能量平均到整个频段上等你能在频谱图上看出边频带时故障通常已经到了中后期。这也是为什么工业现场越来越倾向于用希尔伯特–黄变换HHT这类时频分析方法做前端特征提取——它不假设信号平稳先用经验模态分解EMD把信号自适应地拆成若干个本征模态函数IMF再做希尔伯特变换得到瞬时频率和瞬时幅值把“什么时候出现异常冲击”这件事显式地表达出来。但 HHT 只解决了特征表达的问题分类还要靠模型。极限学习机ELM在这里的价值是训练极快隐层权重随机生成、输出权重闭式解求解不需要反向传播迭代特别适合故障样本维度高、需要频繁更新模型的场景。这套 HHT-ELM 组合在滚动轴承、齿轮箱、电机等设备的故障诊断中很有代表性本文按“信号分解 → 时频特征构造 → ELM 分类建模 → 调参与部署”的完整链路展开所有代码基于 MATLAB 环境可直接对照项目实例复现。无论是做科研验证还是搭建实际的状态监测系统这条技术路线都值得完整走一遍。2. EMD 分解与希尔伯特分析把非平稳信号拆成能用的特征2.1 EMD 的自适应分解逻辑与 IMF 判据经验模态分解的核心假设是任何复杂信号都可以看成若干个本征模态函数与一个残差项的叠加。IMF 需要满足两个条件一是整个数据段内极值点数目与过零点数目相等或至多相差一个二是任意位置处由局部极大值定义的上包络和由局部极小值定义的下包络均值为零。分解过程是一个迭代筛选sifting算法找出原始信号x(t)的全部局部极值点用三次样条插值分别拟合上包络u(t)和下包络l(t)计算包络均值m(t) (u(t)l(t))/2提取细节分量h(t) x(t) - m(t)检查h(t)是否满足 IMF 条件不满足则将其作为新输入重复以上步骤满足则记为第一个 IMF用原信号减去该 IMF对残差继续执行上述过程直到残差为单调函数或幅值小于阈值。function [imfs, residual] emd_decompose(x, max_imf) % 简化版 EMD实际工程建议直接调用 MATLAB 自带 emd() 函数 imfs []; residual x(:); for k 1:max_imf h residual; sd 1; % 内循环筛选通过包络均值不断逼近 IMF 条件 while sd 0.1 size(unique(h), 1) 3 env_upper interp1(find(islocalmax(h)), h(islocalmax(h)), 1:length(h), spline); env_lower interp1(find(islocalmin(h)), h(islocalmin(h)), 1:length(h), spline); env_mean (env_upper env_lower) / 2; h_new h - env_mean; sd sum((h - h_new).^2) / sum(h.^2 eps); % 标准差判据 h h_new; end imfs [imfs; h]; residual residual - h; if length(unique(residual)) 3 break; end end end这段代码里最关键的是sd这个停止条件它计算筛选前后信号的归一化能量差工程上通常取 0.1 到 0.3 之间。值取太小会导致过度筛选IMF 变成纯正弦波失去物理意义取太大则筛得不干净IMF 里还混着其他尺度成分。另外unique(h)的判断是为了处理信号退化成常数的边界情况。2.2 Hilbert 变换与瞬时特征提取每一个 IMF 都是窄带信号对它做 Hilbert 变换可以得到解析信号function [inst_freq, inst_amp] hilbert_features(imf, fs) analytic hilbert(imf); % 解析信号 inst_amp abs(analytic); % 瞬时幅值包络 inst_phase unwrap(angle(analytic)); % 瞬时相位unwrap 消除跳变 inst_freq diff(inst_phase) / (2 * pi) * fs; % 瞬时频率单位 Hz inst_freq(end 1) inst_freq(end); % 对齐长度 end注意diff之后序列会缩短一个点这里用重复最后一个值的方式对齐实际应用中更推荐对瞬时相位做多项式拟合后再求导能有效抑制端点处的频率抖动。提取出的瞬时幅值能反映故障冲击的强度变化瞬时频率则能体现调制特征——比如齿轮断齿时啮合频率附近会出现以转频为间隔的边频带而这些在瞬时频率曲线上会表现为周期性的频率牵引。2.3 为什么 HHT 比小波更适合故障冲击信号小波变换需要预先选择小波基函数和分解层数选错了基特征表达就会偏差而 HHT 的基函数是从信号自身自适应生成的不依赖先验假设。对于轴承故障那种持续时间短、衰减快的冲击响应HHT 能在时间轴上准确定位冲击发生的时刻同时给出对应的瞬时频率成分。提示EMD 的端点效应是绕不开的问题。信号两端在包络拟合时缺少极值点约束容易产生发散。处理方法是做信号延拓常见做法是在两端各延长 3~5 个周期分解完成后裁剪掉延拓部分。MATLAB 自带的emd()函数内部已处理此事自研实现时务必考虑。3. 特征向量构造从 IMF 到能喂给分类器的数值表达3.1 时域、频域与时频特征的组合策略EMD 分解后通常得到 5~10 个 IMF如果把所有 IMF 的完整波形都作为特征维度太高且冗余严重。工程上常用的做法是从每个 IMF 中提取统计量再结合 Hilbert 边际谱构造紧凑的特征向量。我一般会从以下三个层面取特征时域层每个 IMF 的均方根值、峭度、峰值因子、脉冲因子。峭度对冲击型故障非常敏感正常轴承的峭度在 3 附近出现点蚀后会明显升高频域层对每个 IMF 做 FFT提取重心频率、频率方差、主频幅值占比时频层对 Hilbert 边际谱计算重心频率和能量集中度边际谱是瞬时幅值平方对时间积分的结果反映各频率成分在整个信号持续期内的累积能量分布。function feature_vec extract_features(imfs, fs) num_imf size(imfs, 1); feat []; for i 1:num_imf imf imfs(i, :); % 时域特征 rms_val sqrt(mean(imf.^2)); kurt_val kurtosis(imf); % 峭度默认已中心化 peak_val max(abs(imf)); feat [feat, rms_val, kurt_val, peak_val / (rms_val eps)]; % 频域特征 fft_amp abs(fft(imf)); f (0:length(imf)-1) * fs / length(imf); centroid sum(f .* fft_amp) / sum(fft_amp eps); feat [feat, centroid, std(fft_amp)]; end % 取前 6 个 IMF 的特征不足则补零控制维度稳定 model_dim 6 * 5; if length(feat) model_dim feat [feat, zeros(1, model_dim - length(feat))]; else feat feat(1:model_dim); end feature_vec feat; end这段代码里特征维度的上限控制很关键。不同样本 EMD 分解出的 IMF 数量不一致如果不截断或补零特征向量的长度就不统一ELM 的输入节点数没法固定。取前 6 个 IMF 是基于能量占比的考虑——一般前 3~4 个 IMF 已携带 90% 以上的信号能量后续 IMF 主要是低频残余和趋势项区分度有限。3.2 特征标准化与样本划分特征提取完成后需要标准化。ELM 的输入层权重是随机生成的如果特征量纲差异过大输出权重的求解会受大数值特征主导小数值但判别力强的特征会被稀释。function [X_train, X_test, y_train, y_test] prepare_data(all_feat, all_label, ratio) % all_feat: N x M 特征矩阵all_label: N x 1 类别标签 rng(42); % 固定随机种子保证可复现 idx randperm(size(all_feat, 1)); num_train round(length(idx) * ratio); train_idx idx(1:num_train); test_idx idx(num_train1:end); X_train_raw all_feat(train_idx, :); X_test_raw all_feat(test_idx, :); % 标准化参数只从训练集计算防止测试集信息泄露 mu mean(X_train_raw); sigma std(X_train_raw); sigma(sigma 0) 1; % 常数列不参与缩放 X_train (X_train_raw - mu) ./ sigma; X_test (X_test_raw - mu) ./ sigma; y_train all_label(train_idx); y_test all_label(test_idx); end标准化参数只用训练集计算这一点很多人会忽略。如果在标准化时用了全量数据的均值和方差测试集的信息就已经混进了训练流程评估结果会偏乐观现场部署时性能会打折扣。特征类别具体指标对故障类型的敏感性时域统计峭度、峰值因子、脉冲因子冲击类故障轴承点蚀、齿轮断齿频域统计重心频率、频率方差磨损类故障频谱重心偏移时频分析边际谱能量集中度、IMF 能量占比调制类故障齿轮啮合边频带4. ELM 分类器原理与 MATLAB 实现随机映射背后的数学逻辑4.1 ELM 的核心思想不需要迭代的训练ELM 的结构是单隐层前馈网络。输入层到隐层的权重W和偏置b随机生成后固定不动隐层输出通过激活函数计算然后求解输出层权重β。网络输出可以写成矩阵形式H · β T其中H是隐层输出矩阵维度为N x LN 个样本L 个隐层节点T是目标标签矩阵。由于H是已知的β可以通过最小二乘直接求解β H† · T这里的H†是H的 Moore-Penrose 广义逆。整个训练过程没有反向传播没有梯度下降因此训练速度比 BP 网络快一两个数量级。function model elm_train(X, y, hidden_num, act_type) % X: 训练特征矩阵 N x My: 标签 N x 1 % 标签转 one-hot 编码 classes unique(y); num_class length(classes); T zeros(length(y), num_class); for i 1:length(y) T(i, find(classes y(i))) 1; end % 随机生成输入权重和偏置 rng(7); input_w rand(size(X, 2), hidden_num) * 2 - 1; % [-1, 1] 均匀分布 bias rand(1, hidden_num) * 2 - 1; % 计算隐层输出矩阵 H X * input_w repmat(bias, size(X, 1), 1); switch act_type case sigmoid H 1 ./ (1 exp(-H)); case relu H max(0, H); case tanh H tanh(H); end % 广义逆求解输出权重加入正则化项防止过拟合 lambda 1e-3; beta (H * H lambda * eye(hidden_num)) \ H * T; model.beta beta; model.input_w input_w; model.bias bias; model.act_type act_type; model.classes classes; end求解β时用的公式是(HH λI)⁻¹HT这相当于在最小二乘目标函数里加了一个 L2 正则项λ||β||²。工程上λ取1e-3到1e-5比较稳妥。正则化在故障诊断场景里很重要——故障样本往往数量有限隐层节点数又比较多不加正则化的最小二乘解会过度拟合训练集测试集准确率明显下降。4.2 隐层节点数的选择逻辑隐层节点数L是 ELM 唯一需要认真调的超参数。理论上隐层节点越多模型的拟合能力越强但超过某个临界点后泛化性能会下降。常见做法是做一个简单的扫描实验从 10 到 500 按步长 10 递增分别训练并记录验证集准确率选择准确率曲线平台期的节点数。function best_node search_hidden_node(X_train, y_train, X_val, y_val) node_list 10:10:500; acc_list zeros(size(node_list)); for i 1:length(node_list) model elm_train(X_train, y_train, node_list(i), sigmoid); pred elm_predict(model, X_val); acc_list(i) mean(pred y_val); end [~, best_idx] max(acc_list); best_node node_list(best_idx); % 绘制曲线辅助观察 figure; plot(node_list, acc_list, -o); xlabel(Hidden Nodes); ylabel(Validation Accuracy); grid on; end这里要注意验证集和测试集必须分开。如果用测试集来选隐层节点数测试集就变成了验证集最终评估的准确率会有偏。4.3 预测与概率输出function [pred_label, prob] elm_predict(model, X) H X * model.input_w repmat(model.bias, size(X, 1), 1); switch model.act_type case sigmoid H 1 ./ (1 exp(-H)); case relu H max(0, H); case tanh H tanh(H); end output H * model.beta; % N x num_class 的连续值 % softmax 转概率 exp_out exp(output - max(output, [], 2)); prob exp_out ./ sum(exp_out, 2); [~, idx] max(prob, [], 2); pred_label model.classes(idx); end代码里exp(output - max(output))是数值稳定的 softmax 写法直接对output做 exp 可能溢出。输出的prob可以用来做置信度筛选——当最高类别的概率低于某个阈值比如 0.7时标记为“不确定”而不是强行归类这在工业现场比单纯追求准确率更实用。5. 完整故障诊断流程整合从原始信号到评估报告5.1 模拟信号构造与数据生成项目的第一步是构造模拟数据来验证整个技术路线。这里模拟四种状态正常、轴承内圈故障、外圈故障、滚动体故障。每种状态的振动信号模型不同例如内圈故障会在转频调制下产生周期性冲击function data generate_fault_signal(fault_type, fs, duration) % fs: 采样率 Hz, duration: 信号时长 s t 0:1/fs:duration-1/fs; n length(t); % 工频转频与故障特征频率 fr 25; % 转频 25 Hz对应 1500 RPM f_fault 95; % 外圈故障特征频率 BPFO 示例 rng(fault_type * 10); switch fault_type case 1 % 正常 data 0.5 * sin(2*pi*fr*t) 0.1 * randn(1, n); case 2 % 内圈故障 BPFI impact 0.8 * exp(-50 * mod(t, 1/f_fault)) .* sin(2*pi*2000*mod(t, 1/f_fault)); data 0.5*sin(2*pi*fr*t) impact .* (1 0.5*sin(2*pi*fr*t)) 0.1*randn(1, n); case 3 % 外圈故障 BPFO impact 0.9 * exp(-60 * mod(t, 1/f_fault)) .* sin(2*pi*1800*mod(t, 1/f_fault)); data 0.5*sin(2*pi*fr*t) impact 0.1*randn(1, n); case 4 % 滚动体故障 BSF impact 0.7 * exp(-40 * mod(t, 1/60)) .* sin(2*pi*1500*mod(t, 1/60)); data 0.5*sin(2*pi*fr*t) impact .* (1 0.8*sin(2*pi*12*t)) 0.1*randn(1, n); end end将每类故障生成 50 段样本每段 1 秒构造出总共 200 个样本的数据集。实际项目中数据来源是传感器采集但验证算法链路时人工构造信号能精确控制故障特征便于对比不同参数下的表现差异。5.2 特征提取批量处理与数据集构建对每个样本执行 EMD 分解、特征提取最终得到200 x 30的特征矩阵每个样本 6 个 IMF × 5 个统计量。然后用prepare_data函数按 7:3 划分训练集和测试集。5.3 多维度评估与可视化ELM 训练完成后需要从多个维度评估分类效果。准确率只能反映整体表现对故障诊断来说每类故障的召回率同样重要——把“内圈故障”误判成“正常”意味着漏报一次真实故障代价远比把“正常”误判成“故障”高。function report evaluate_model(pred_label, true_label, class_names) % 混淆矩阵 cm confusionmat(true_label, pred_label); num_class length(class_names); % 精确率、召回率、F1 precision zeros(num_class, 1); recall zeros(num_class, 1); f1 zeros(num_class, 1); for i 1:num_class tp cm(i, i); fp sum(cm(:, i)) - tp; fn sum(cm(i, :)) - tp; precision(i) tp / (tp fp eps); recall(i) tp / (tp fn eps); f1(i) 2 * precision(i) * recall(i) / (precision(i) recall(i) eps); end report table(class_names, precision, recall, f1); % 绘图 figure; confusionchart(cm, class_names); title(ELM 故障分类混淆矩阵); end在 MATLAB 中confusionchart可以直接生成带颜色映射的混淆矩阵图数值格上的颜色深浅代表样本数量级。热力图对识别“哪些类别容易混淆”非常直观比如内圈故障和外圈故障之间的特征很接近时混淆矩阵上对应位置的数值会偏高提醒我们可能需要增加额外特征来区分这两类。5.4 交叉验证与模型稳定性验证单次划分训练集和测试集的评估结果波动较大特别是样本数量少的时候。项目里我建议至少做 5 折交叉验证。基本做法是将数据分成 5 份轮流取其中 4 份训练、1 份验证最终报告 5 次结果的平均值和标准差。标准差能反映模型对数据划分的敏感性——如果标准差超过 2 个百分点说明特征或模型参数还不够稳定。评估指标公式故障诊断中的意义准确率TPTN / 总样本整体分类正确比例精确率TP / (TPFP)报警可靠性虚警率越低越好召回率TP / (TPFN)漏报率控制关键设备必须高F1 值2PR / (PR)精确率与召回率的平衡度量宏平均各类指标取算术平均小样本类别不被大类别淹没5.5 GUI 设计要点从脚本到可交付工具项目附带的 MATLAB GUI 把以上流程全部封装进界面。GUI 的框架通常是左侧控制面板放置“载入数据”“特征提取”“模型训练”“预测评估”四个按钮右侧展示区划分数据概览表格、特征分布图窗、训练状态信息区和评估结果展示区。GUI 回调函数的核心是状态机管理。例如点击“特征提取”按钮之前必须先检查数据是否已经载入点击“模型训练”之前必须确认特征矩阵已经生成。用一组标志位isDataLoaded、isFeatureExtracted来控制按钮的可用性能避免误操作导致程序崩溃。6. 排查与实战进阶把 HHT-ELM 模型用到真实设备信号上真实工业信号远比模拟信号复杂项目从实验室走向现场时通常会碰到几个具体问题这里给出对应的处理策略。第一是端点效应。EMD 分解时信号两端发散会导致 IMF 在边界处失真。常见做法是信号延拓——在两端各延长数个周期后分解再裁剪掉延拓部分。MATLAB 自带的emd()函数内部有端点处理机制但用自定义 EMD 实现时这一步不能省。另一种思路是改用集合经验模态分解EEMD通过多次叠加白噪声取平均来抑制模态混叠代价是计算量增加 10 倍以上。第二是模态混叠。当故障特征频率与信号中的其他频率成分过于接近时相邻 IMF 之间会出现频率交叉。排查方法是逐个 IMF 做频谱分析如果发现两个 IMF 的主频重叠说明分解不干净。实际项目中增加 EMD 筛选迭代次数、或者对信号先做带通滤波保留关注频段都能改善模态混叠问题。第三是降低虚警率的置信度策略。ELM 的 softmax 输出可以转化为置信度现场系统可以设置双阈值置信度高于 0.9 直接告警介于 0.7 和 0.9 之间标记为“关注”低 0.7 则归为“不确定”等待下一轮采样再判断。这样不会频繁误报也不会漏掉真正的早期故障。关于超参数调优ELM 比深度网络简单很多。隐层节点数通过网格扫描找到准确率平台区间的下限值即可激活函数在sigmoid和tanh之间二选一。relu在 ELM 中效果不稳定因为随机生成的输入权重可能导致大量隐层节点输出恒为零反而损失表达能力。故障诊断项目的价值在于把信号处理方法和分类模型组合成一条完整可用的技术链路HHT-ELM 只是其中一种被验证过有效的方式。当你把代码跑通、把混淆矩阵看明白、把 GUI 界面组装完成之后这套“分解-特征-分类”的框架完全可以迁移到电流信号分析、声发射检测等其他场景中。动手改一改特征参数换一组真实数据你会比只看代码更快理解这套组合到底在解决什么问题。本文还有配套的精品资源点击获取