简介围绕非平稳振动信号下的故障识别需求这套MATLAB项目实例面向具备一定编程基础的信号处理、工业自动化与智能运维研发人员工作1—3年的工程师可借此打通时频分析与机器学习结合的落地路径。内容从振动信号采集与预处理、STFT时频图生成与频带能量特征提取延伸到SVM多类别分类、交叉验证调参、混淆矩阵与单样本时频图联动展示并给出GUI界面的控件回调与参数调节逻辑构成可复现、可扩展的原型系统。压缩包共1个docx文件以图文文档形式承载完整程序清单、参数配置说明与代码逐段解析体积约120KB便于离线阅读与对照调试。目前已有132人学习下载。读者可直接获得可运行的工程代码框架、STFT参数对特征稳定性影响的调试思路以及为后续引入深度学习、多传感器融合预留的技术验证平台。1. 一台变转速电机上的反常现象频谱看不出毛病时频图却全是故事一台 1800 rpm 的电机驱动轴承座内圈滚道出现早期剥落。现场采一段振动信号丢进 FFT包络谱上只看到几条模糊的谱线峰值因子也只有 3.2跟健康样本的差别不到 10%。但如果把同一段信号按 256 点汉明窗、75% 重叠做一次 STFT时频图上每隔约 5.5 ms 就出现一串沿 3 kHz 共振带衰减的竖直条纹——那正是内圈故障特征频率 BPFI 在时间轴上的投影。原因不复杂。整段傅里叶变换把时间信息积分掉了非平稳信号里那些短促的冲击被平均进背景噪声而 STFT 用滑动窗把长信号切成几十毫秒的短片段逐段做变换冲击发生在哪一刻、能量集中在哪条共振带就都能在时间-频率平面上留下来。这套流程最终要落到「分类」上把时频矩阵聚合成一维特征向量交给支持向量机做多类别判别。对做旋转机械状态监测的人来说这条路的价值在于参数量少、训练成本低、结果可解释跑到边缘计算盒子上也不吃力。下面按数据生成、时频特征、模型调参、GUI 集成四段拆开讲。2. 模拟振动信号生成与预处理把故障机理写成可复现的数据源真实故障样本稀缺是这类项目的第一个坎。手头只有几段现场录的波形四类状态加起来不到 40 个样本直接训练 SVM 一定是过拟合。常见做法是先按物理机理合成一批带标签的仿真信号把整条流水线跑通、参数调稳再用实测数据替换。这一节的重点不是合成数据本身而是让合成模型和后面的特征提取对得上。2.1 四类健康状态的合成模型轴承类故障的时域波形可以用「周期冲击串 × 结构共振衰减响应」描述冲击重复频率由故障特征频率决定每次冲击激发系统共振并以指数规律衰减。转子不平衡则是转频 1X 及其谐波的能量突出。据此把四类状态分别建模状态标签时域模型要点主要特征频率健康01X 转频 微弱 2X 宽带噪声fr外圈剥落1等间隔冲击串 转频分量BPFO ≈ 4.2 fr内圈剥落2冲击串经载荷区调制BPFI ≈ 6.5 fr边带 ±fr转子不平衡31X 强烈、2X 明显、相位固定fr、2 fr采样率必须覆盖系统共振频率的 3 倍以上否则共振带能量会被折叠进低频区时频图上的竖直条纹直接消失。这里取 fs 12.8 kHz共振频率 fn 3 kHz衰减系数 ζ 0.05——ζ 太小冲击拖尾过长会糊掉相邻两次冲击的间隔。2.2 数据生成代码function [X, Y, t] genFaultDataset(nPerClass, fs, snrDb) % 生成四类健康状态的振动样本集 % 输出 X: [4*nPerClass x N] 按行存放样本; Y: 标签列向量; t: 时间轴 rng(2024); dur 0.5; % 单样本 0.5 s t (0:1/fs:dur-1/fs); N numel(t); fr 30; % 转频 30 Hz - 1800 rpm bpfo 4.2*fr; bpfi 6.5*fr; % 外圈/内圈故障特征频率 fn 3000; zeta 0.05; % 结构共振频率与衰减系数 nTot 4*nPerClass; X zeros(nTot, N); Y zeros(nTot, 1); for c 0:3 for k 1:nPerClass idx c*nPerClass k; switch c case 0 % 健康转频及微弱谐波 x 0.30*sin(2*pi*fr*t) 0.08*sin(2*pi*2*fr*t); case 1 % 外圈等间隔冲击串 x impulseTrain(t, bpfo, fn, zeta) 0.20*sin(2*pi*fr*t); case 2 % 内圈冲击串带载荷区调制 x impulseTrain(t, bpfi, fn, zeta) .* (1 0.6*cos(2*pi*fr*t)); otherwise % 不平衡1X 主导 2X 次谐波 x 1.00*sin(2*pi*fr*t) 0.35*sin(2*pi*2*fr*t pi/4); end x x(:) randn(N,1) * 10^(-snrDb/20) * rms(x); % 按目标 SNR 加高斯白噪 X(idx,:) x.; Y(idx) c; end end end function y impulseTrain(t, fImp, fn, zeta) % 周期性冲击激发的共振衰减响应 y zeros(size(t)); T 1/fImp; % 冲击间隔 for tp 0:T:t(end) tt t - tp; mask tt 0; % 只保留冲击发生之后的时段 resp zeros(size(t)); resp(mask) exp(-zeta*2*pi*fn*tt(mask)) .* sin(2*pi*fn*tt(mask)); y y resp; end end逐段说明一下case 2里乘的那个(1 0.6*cos(2*pi*fr*t))就是载荷区调制它会在时频图上给 BPFI 主频两侧拉出 ±fr 的边带这是内圈和外圈最容易混淆的地方也是后面特征设计要专门留意的点。impulseTrain里用mask掩码而不是给负时间赋inf是为了避开exp(-inf)*sin(inf)产生的 NaN——这个坑我第一次写的时候踩过训练时整批特征全是 NaN很容易误以为是归一化写错了。2.3 预处理去趋势、零相位带通与幅度归一化function Xp preprocess(X, fs) % 逐样本去趋势 - 带通 - 峰值归一化 bp designfilt(bandpassiir,FilterOrder,6, ... HalfPowerFrequency1,500,HalfPowerFrequency2,5000, ... SampleRate,fs); Xp zeros(size(X)); for i 1:size(X,1) x detrend(X(i,:)); % 去掉直流偏置和线性漂移 x filtfilt(bp, x); % 零相位滤波 Xp(i,:) x / max(abs(x)); % 峰值归一化消除通道增益差异 end end三个参数各有理由。带通下限取 500 Hz 是为了滤掉转频及其低次谐波对共振带的干扰上限 5000 Hz 留出抗混叠余量。滤波器用filtfilt而不是filter前者前后各滤波一次实现零相位冲击包络的位置不会发生群延迟偏移——对时频分析来说冲击落在哪一帧直接决定特征值延迟几毫秒就够毁掉可分性。归一化按样本峰值做而不是全局做是为了保留样本之间的能量差异同时又不受传感器增益漂移影响。3. STFT 时频特征构造窗长、重叠率与频带聚合怎么定STFT 本身的调用只有一行难的是参数怎么定、高维矩阵怎么压。这一节把这两个问题都落到具体数字上。3.1 时间分辨率与频率分辨率的定量取舍STFT 的核心矛盾是海森堡不确定性窗长越长频率分辨率越高但时间定位越模糊。工程上有个简单的估算方式频率分辨率 Δf ≈ fs / Nw × 窗因子汉明窗约 1.36时间步长 Δt (Nw - noverlap) / fs按 fs 12.8 kHz 算几组窗长 NwΔf汉明窗时间步长75% 重叠适用场景64约 272 Hz1.25 ms高转速、冲击密集重时间定位256约 68 Hz5.0 ms通用轴承诊断均衡1024约 17 Hz20 ms齿轮啮合边带、低频调制分析内圈故障的冲击间隔在 5 ms 量级用 256 点窗刚好能分辨出每一次冲击而不至于把相邻两次糊在一起。窗口重叠率取 75% 是经验值再高收益递减、计算量线性上涨。如果后续要改动窗长做对比实验记得把spectrogram的缩放模式从power换成psd否则不同窗长下的能量值不可比。3.2 频带聚合把 N×M 矩阵压成 1×K 特征时频谱直接拉平做特征会有上千维样本才几百个必然维度灾难。做法是按物理意义划出若干频带在带内做统计聚合。function fv stftFeature(x, fs, edges, Nw, ov, nfft) % 单样本 STFT 频带聚合特征 win hamming(Nw); [~, F, ~, P] spectrogram(x(:), win, ov, nfft, fs, psd); nb numel(edges) - 1; eMean zeros(1,nb); eStd zeros(1,nb); ePk zeros(1,nb); for b 1:nb idx F edges(b) F edges(b1); Pb P(idx,:); pt mean(Pb, 1); % 先沿频率轴聚合得到带内能量时间序列 eMean(b) log10(mean(pt) eps); % 频带平均能量取对数压量纲 eStd(b) log10(std(pt) eps); % 沿时间的起伏冲击周期性越强值越大 ePk(b) max(pt) / (mean(pt) eps); % 峰均比对短时冲击敏感 end fCent sum(F .* mean(P,2)) / (sum(mean(P,2)) eps); % 能量重心 kurt kurtosis(x(:)); % 时域峭度作补充 fv [eMean, eStd, ePk, fCent, kurt]; end关键点在pt mean(Pb, 1)这一步先沿频率轴把带内所有谱线平均掉得到一条随时间变化的能量曲线再做时间轴的统计。这样做的好处是噪声在频率轴上的随机起伏被平均抑制而冲击在时间轴上留下的尖峰却被保留下来——eStd和ePk正是靠这个机制对冲击敏感。频带边界按共振带划edges [500 1200 2000 2600 3200 4000 5000]其中 2600–3200 Hz 正好覆盖 3 kHz 共振带故障能量的大头都在这里面。如果换了设备先画几张典型样本的时频图看能量集中在哪里再回头调 edges这一步没法靠公式代替。3.3 用 Fisher 判别比快速筛特征特征提出来有二十多列哪几列真正有区分度用类间方差比类内方差算一下就有答案。function score fisherScore(F, y) % 逐维计算 Fisher 判别比值越大说明该维的类间可分性越强 cls unique(y); score zeros(1, size(F,2)); for j 1:size(F,2) muAll mean(F(:,j)); sb 0; sw 0; for c cls fc F(yc, j); sb sb numel(fc) * (mean(fc) - muAll)^2; % 类间散度 sw sw sum((fc - mean(fc)).^2); % 类内散度 end score(j) sb / (sw eps); end end跑完通常会看到两类结果共振带上的eMean和ePk得分最高fCent和kurt得分中等而 500–1200 Hz 那几个低频带的eStd得分普遍偏低——低频段基本是转频分量四类状态都有区分度自然差。把得分低于中位数的维度直接砍掉特征维度从 24 降到 12 左右训练速度翻倍而精度几乎不掉。4. SVM 分类器训练与调参核函数、交叉验证与类别不平衡处理特征准备完之后模型侧要处理三件事数据怎么划、核函数和超参数怎么定、类别不平衡怎么补偿。4.1 分层划分与标准化cv cvpartition(Y, HoldOut, 0.3, Stratify, true); Xtr Ftr(cv.training,:); Ytr Y(cv.training); Xte Ftr(cv.test,:); Yte Y(cv.test);Stratify, true必须开。四类样本各 100 个随机划分有可能让某一类在测试集里只剩十几个评估结果波动极大。分层划分保证每类在训练集和测试集里的比例一致。标准化交给templateSVM的Standardize参数做它会在训练时记录均值方差、预测时复用同一套参数比手工zscore再手动还原要安全。4.2 核函数选择与多类组合方式线性核适合特征已经线性可分的情况速度最快RBF 核能处理非线性边界是振动特征分类的默认选择。多类问题 MATLAB 用fitcecoc封装编码方式有三种编码含义分类器数量4 类特点onevsone两两配对6每个分类器样本均衡综合表现最稳onevsall一对一4训练快类别多时易受不平衡影响binarycomplete纠错输出码7容错性强计算量最大4 类场景下 onevsone 是性价比最高的6 个二分类器每个只用两类样本天然规避了类别不平衡对单个分类器的影响。4.3 超参数网格与自动优化惩罚系数 BoxConstraint 控制间隔与误分类的权衡KernelScale 决定 RBF 核的作用半径。手工网格搜索是这样写的boxList [0.5 2 8 32]; % 对数尺度递增 sigList [0.5 1 2 4]; bestAcc 0; bestT []; for bc boxList for ks sigList t templateSVM(KernelFunction,rbf,BoxConstraint,bc, ... KernelScale,ks,Standardize,true); m fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone); L kfoldLoss(crossval(m, KFold, 5)); % 5 折交叉验证误分类率 if 1-L bestAcc bestAcc 1-L; bestT t; end end end16 组参数各跑 5 折样本量上千的时候耗时几分钟可以接受。如果 matlab优化工具箱可用直接交给贝叶斯优化更省事opts struct(Optimizer,bayesopt,ShowPlots,false, ... MaxObjectiveEvaluations,30, ... CVPartition, cvpartition(Ytr,KFold,5)); t templateSVM(KernelFunction,rbf,Standardize,true); mdl fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone, ... OptimizeHyperparameters,{BoxConstraint,KernelScale}, ... HyperparameterOptimizationOptions, opts);注意KernelScale和Standardize一起用时优化器搜的是标准化之后的数据尺度。如果跳过Standardize手动标准化KernelScale的最优值会完全不一样两者的搜索结果不能混着用。4.4 类别不平衡与评估指标现场数据里严重故障样本往往只有几十个轻故障几百个。两个手段一是templateSVM里设Prior,uniform让各类先验概率相等二是直接给少数类加权w ones(size(Ytr)); w(Ytr 3) 2.5; % 不平衡类样本权重抬高 t templateSVM(KernelFunction,rbf,BoxConstraint,8, ... KernelScale,auto,Standardize,true); mdl fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone, Weights, w);权重取多少可以用反频率比初估w_c N_total / (K * N_c)。补偿过度会让少数类误报率上升需要在混淆矩阵上反复看。评估别只看准确率。四类样本均衡时准确率够用一旦不均衡一个把所有样本都判成多数类的模型也能拿到 70% 准确率。Yp predict(mdl, Xte); C confusionmat(Yte, Yp); acc sum(diag(C)) / sum(C(:)); % 宏平均 F1先算每类的 P/R再等权平均不被类别规模带偏 K numel(unique(Yte)); f1 zeros(K,1); for c 1:K tp C(c,c); fp sum(C(:,c)) - tp; fn sum(C(c,:)) - tp; p tp/(tpfpeps); r tp/(tpfneps); f1(c) 2*p*r/(preps); end macroF1 mean(f1);如果混淆矩阵上看到内圈和外圈互相误判特别多说明载荷区调制带来的边带特征没被抓住——回到第 3 节把 2600–3200 Hz 共振带再细分成两个子带重新提特征通常比继续调 SVM 参数有效。5. GUI 集成与 R2025b 图形规范三个容易翻车的细节把上面几个函数串起来做成界面核心是参数传递和图形调用规范。用uifigureuigridlayout搭骨架控件回调写成嵌套函数共享父函数工作区里的模型句柄比用guidata存取清爽得多。训练这种耗时操作放进parfeval后台执行用afterEach回调刷新界面否则点一次「训练模型」窗口会卡十几秒。R2025b 有三个改动会直接让代码报错% 1) colormap 需要显式指定目标 figure colormap(fig, turbo); % 旧写法 colormap(turbo) 在多窗口下目标不明确 % 2) colorbar 的属性名变了别再设 ColorbarVisible cb colorbar(ax); cb.Label.String 功率 / dB; % 用 Label 对象而不是旧的可视性开关 % 3) confusionchart 不能当普通子级颜色要设给坐标区 cm confusionchart(ax, Yte, Yp); cm.Title 分类混淆矩阵;时频图的绘制顺序也有讲究先imagesc画功率谱再set(ax,YDir,normal)让频率轴从下往上递增最后调colormap和colorbar。顺序反了会出现颜色映射与坐标区不匹配的情况图看着没问题但刻度颜色对不上。坐标轴传的是uiaxes句柄而不是gca因为 GUI 里有多个子图时gca拿到的是当前焦点轴未必是你想画的那个。最后一个实用技巧实时预测场景下别每次重新提特征。把频带索引预先算好存成结构体流式数据进来后直接查表切片persistent bandIdx if isempty(bandIdx) [~, F] spectrogram(zeros(256,1), hamming(256), 192, 512, 12800, psd); edges [500 1200 2000 2600 3200 4000 5000]; bandIdx arrayfun((b) find(Fedges(b) Fedges(b1)), ... 1:numel(edges)-1, UniformOutput, false); end这样省掉的是每个样本都对 F 做一次逻辑比较的开销在 12.8 kHz 采样、每秒出一个诊断结果的在线场景里单样本特征提取从 38 ms 降到 24 ms 左右——对需要跑几十路通道的边缘盒子来说这个差距决定了能不能用单核跑完整个产线。本文还有配套的精品资源点击获取
简介围绕非平稳振动信号下的故障识别需求这套MATLAB项目实例面向具备一定编程基础的信号处理、工业自动化与智能运维研发人员工作1—3年的工程师可借此打通时频分析与机器学习结合的落地路径。内容从振动信号采集与预处理、STFT时频图生成与频带能量特征提取延伸到SVM多类别分类、交叉验证调参、混淆矩阵与单样本时频图联动展示并给出GUI界面的控件回调与参数调节逻辑构成可复现、可扩展的原型系统。压缩包共1个docx文件以图文文档形式承载完整程序清单、参数配置说明与代码逐段解析体积约120KB便于离线阅读与对照调试。目前已有132人学习下载。读者可直接获得可运行的工程代码框架、STFT参数对特征稳定性影响的调试思路以及为后续引入深度学习、多传感器融合预留的技术验证平台。1. 一台变转速电机上的反常现象频谱看不出毛病时频图却全是故事一台 1800 rpm 的电机驱动轴承座内圈滚道出现早期剥落。现场采一段振动信号丢进 FFT包络谱上只看到几条模糊的谱线峰值因子也只有 3.2跟健康样本的差别不到 10%。但如果把同一段信号按 256 点汉明窗、75% 重叠做一次 STFT时频图上每隔约 5.5 ms 就出现一串沿 3 kHz 共振带衰减的竖直条纹——那正是内圈故障特征频率 BPFI 在时间轴上的投影。原因不复杂。整段傅里叶变换把时间信息积分掉了非平稳信号里那些短促的冲击被平均进背景噪声而 STFT 用滑动窗把长信号切成几十毫秒的短片段逐段做变换冲击发生在哪一刻、能量集中在哪条共振带就都能在时间-频率平面上留下来。这套流程最终要落到「分类」上把时频矩阵聚合成一维特征向量交给支持向量机做多类别判别。对做旋转机械状态监测的人来说这条路的价值在于参数量少、训练成本低、结果可解释跑到边缘计算盒子上也不吃力。下面按数据生成、时频特征、模型调参、GUI 集成四段拆开讲。2. 模拟振动信号生成与预处理把故障机理写成可复现的数据源真实故障样本稀缺是这类项目的第一个坎。手头只有几段现场录的波形四类状态加起来不到 40 个样本直接训练 SVM 一定是过拟合。常见做法是先按物理机理合成一批带标签的仿真信号把整条流水线跑通、参数调稳再用实测数据替换。这一节的重点不是合成数据本身而是让合成模型和后面的特征提取对得上。2.1 四类健康状态的合成模型轴承类故障的时域波形可以用「周期冲击串 × 结构共振衰减响应」描述冲击重复频率由故障特征频率决定每次冲击激发系统共振并以指数规律衰减。转子不平衡则是转频 1X 及其谐波的能量突出。据此把四类状态分别建模状态标签时域模型要点主要特征频率健康01X 转频 微弱 2X 宽带噪声fr外圈剥落1等间隔冲击串 转频分量BPFO ≈ 4.2 fr内圈剥落2冲击串经载荷区调制BPFI ≈ 6.5 fr边带 ±fr转子不平衡31X 强烈、2X 明显、相位固定fr、2 fr采样率必须覆盖系统共振频率的 3 倍以上否则共振带能量会被折叠进低频区时频图上的竖直条纹直接消失。这里取 fs 12.8 kHz共振频率 fn 3 kHz衰减系数 ζ 0.05——ζ 太小冲击拖尾过长会糊掉相邻两次冲击的间隔。2.2 数据生成代码function [X, Y, t] genFaultDataset(nPerClass, fs, snrDb) % 生成四类健康状态的振动样本集 % 输出 X: [4*nPerClass x N] 按行存放样本; Y: 标签列向量; t: 时间轴 rng(2024); dur 0.5; % 单样本 0.5 s t (0:1/fs:dur-1/fs); N numel(t); fr 30; % 转频 30 Hz - 1800 rpm bpfo 4.2*fr; bpfi 6.5*fr; % 外圈/内圈故障特征频率 fn 3000; zeta 0.05; % 结构共振频率与衰减系数 nTot 4*nPerClass; X zeros(nTot, N); Y zeros(nTot, 1); for c 0:3 for k 1:nPerClass idx c*nPerClass k; switch c case 0 % 健康转频及微弱谐波 x 0.30*sin(2*pi*fr*t) 0.08*sin(2*pi*2*fr*t); case 1 % 外圈等间隔冲击串 x impulseTrain(t, bpfo, fn, zeta) 0.20*sin(2*pi*fr*t); case 2 % 内圈冲击串带载荷区调制 x impulseTrain(t, bpfi, fn, zeta) .* (1 0.6*cos(2*pi*fr*t)); otherwise % 不平衡1X 主导 2X 次谐波 x 1.00*sin(2*pi*fr*t) 0.35*sin(2*pi*2*fr*t pi/4); end x x(:) randn(N,1) * 10^(-snrDb/20) * rms(x); % 按目标 SNR 加高斯白噪 X(idx,:) x.; Y(idx) c; end end end function y impulseTrain(t, fImp, fn, zeta) % 周期性冲击激发的共振衰减响应 y zeros(size(t)); T 1/fImp; % 冲击间隔 for tp 0:T:t(end) tt t - tp; mask tt 0; % 只保留冲击发生之后的时段 resp zeros(size(t)); resp(mask) exp(-zeta*2*pi*fn*tt(mask)) .* sin(2*pi*fn*tt(mask)); y y resp; end end逐段说明一下case 2里乘的那个(1 0.6*cos(2*pi*fr*t))就是载荷区调制它会在时频图上给 BPFI 主频两侧拉出 ±fr 的边带这是内圈和外圈最容易混淆的地方也是后面特征设计要专门留意的点。impulseTrain里用mask掩码而不是给负时间赋inf是为了避开exp(-inf)*sin(inf)产生的 NaN——这个坑我第一次写的时候踩过训练时整批特征全是 NaN很容易误以为是归一化写错了。2.3 预处理去趋势、零相位带通与幅度归一化function Xp preprocess(X, fs) % 逐样本去趋势 - 带通 - 峰值归一化 bp designfilt(bandpassiir,FilterOrder,6, ... HalfPowerFrequency1,500,HalfPowerFrequency2,5000, ... SampleRate,fs); Xp zeros(size(X)); for i 1:size(X,1) x detrend(X(i,:)); % 去掉直流偏置和线性漂移 x filtfilt(bp, x); % 零相位滤波 Xp(i,:) x / max(abs(x)); % 峰值归一化消除通道增益差异 end end三个参数各有理由。带通下限取 500 Hz 是为了滤掉转频及其低次谐波对共振带的干扰上限 5000 Hz 留出抗混叠余量。滤波器用filtfilt而不是filter前者前后各滤波一次实现零相位冲击包络的位置不会发生群延迟偏移——对时频分析来说冲击落在哪一帧直接决定特征值延迟几毫秒就够毁掉可分性。归一化按样本峰值做而不是全局做是为了保留样本之间的能量差异同时又不受传感器增益漂移影响。3. STFT 时频特征构造窗长、重叠率与频带聚合怎么定STFT 本身的调用只有一行难的是参数怎么定、高维矩阵怎么压。这一节把这两个问题都落到具体数字上。3.1 时间分辨率与频率分辨率的定量取舍STFT 的核心矛盾是海森堡不确定性窗长越长频率分辨率越高但时间定位越模糊。工程上有个简单的估算方式频率分辨率 Δf ≈ fs / Nw × 窗因子汉明窗约 1.36时间步长 Δt (Nw - noverlap) / fs按 fs 12.8 kHz 算几组窗长 NwΔf汉明窗时间步长75% 重叠适用场景64约 272 Hz1.25 ms高转速、冲击密集重时间定位256约 68 Hz5.0 ms通用轴承诊断均衡1024约 17 Hz20 ms齿轮啮合边带、低频调制分析内圈故障的冲击间隔在 5 ms 量级用 256 点窗刚好能分辨出每一次冲击而不至于把相邻两次糊在一起。窗口重叠率取 75% 是经验值再高收益递减、计算量线性上涨。如果后续要改动窗长做对比实验记得把spectrogram的缩放模式从power换成psd否则不同窗长下的能量值不可比。3.2 频带聚合把 N×M 矩阵压成 1×K 特征时频谱直接拉平做特征会有上千维样本才几百个必然维度灾难。做法是按物理意义划出若干频带在带内做统计聚合。function fv stftFeature(x, fs, edges, Nw, ov, nfft) % 单样本 STFT 频带聚合特征 win hamming(Nw); [~, F, ~, P] spectrogram(x(:), win, ov, nfft, fs, psd); nb numel(edges) - 1; eMean zeros(1,nb); eStd zeros(1,nb); ePk zeros(1,nb); for b 1:nb idx F edges(b) F edges(b1); Pb P(idx,:); pt mean(Pb, 1); % 先沿频率轴聚合得到带内能量时间序列 eMean(b) log10(mean(pt) eps); % 频带平均能量取对数压量纲 eStd(b) log10(std(pt) eps); % 沿时间的起伏冲击周期性越强值越大 ePk(b) max(pt) / (mean(pt) eps); % 峰均比对短时冲击敏感 end fCent sum(F .* mean(P,2)) / (sum(mean(P,2)) eps); % 能量重心 kurt kurtosis(x(:)); % 时域峭度作补充 fv [eMean, eStd, ePk, fCent, kurt]; end关键点在pt mean(Pb, 1)这一步先沿频率轴把带内所有谱线平均掉得到一条随时间变化的能量曲线再做时间轴的统计。这样做的好处是噪声在频率轴上的随机起伏被平均抑制而冲击在时间轴上留下的尖峰却被保留下来——eStd和ePk正是靠这个机制对冲击敏感。频带边界按共振带划edges [500 1200 2000 2600 3200 4000 5000]其中 2600–3200 Hz 正好覆盖 3 kHz 共振带故障能量的大头都在这里面。如果换了设备先画几张典型样本的时频图看能量集中在哪里再回头调 edges这一步没法靠公式代替。3.3 用 Fisher 判别比快速筛特征特征提出来有二十多列哪几列真正有区分度用类间方差比类内方差算一下就有答案。function score fisherScore(F, y) % 逐维计算 Fisher 判别比值越大说明该维的类间可分性越强 cls unique(y); score zeros(1, size(F,2)); for j 1:size(F,2) muAll mean(F(:,j)); sb 0; sw 0; for c cls fc F(yc, j); sb sb numel(fc) * (mean(fc) - muAll)^2; % 类间散度 sw sw sum((fc - mean(fc)).^2); % 类内散度 end score(j) sb / (sw eps); end end跑完通常会看到两类结果共振带上的eMean和ePk得分最高fCent和kurt得分中等而 500–1200 Hz 那几个低频带的eStd得分普遍偏低——低频段基本是转频分量四类状态都有区分度自然差。把得分低于中位数的维度直接砍掉特征维度从 24 降到 12 左右训练速度翻倍而精度几乎不掉。4. SVM 分类器训练与调参核函数、交叉验证与类别不平衡处理特征准备完之后模型侧要处理三件事数据怎么划、核函数和超参数怎么定、类别不平衡怎么补偿。4.1 分层划分与标准化cv cvpartition(Y, HoldOut, 0.3, Stratify, true); Xtr Ftr(cv.training,:); Ytr Y(cv.training); Xte Ftr(cv.test,:); Yte Y(cv.test);Stratify, true必须开。四类样本各 100 个随机划分有可能让某一类在测试集里只剩十几个评估结果波动极大。分层划分保证每类在训练集和测试集里的比例一致。标准化交给templateSVM的Standardize参数做它会在训练时记录均值方差、预测时复用同一套参数比手工zscore再手动还原要安全。4.2 核函数选择与多类组合方式线性核适合特征已经线性可分的情况速度最快RBF 核能处理非线性边界是振动特征分类的默认选择。多类问题 MATLAB 用fitcecoc封装编码方式有三种编码含义分类器数量4 类特点onevsone两两配对6每个分类器样本均衡综合表现最稳onevsall一对一4训练快类别多时易受不平衡影响binarycomplete纠错输出码7容错性强计算量最大4 类场景下 onevsone 是性价比最高的6 个二分类器每个只用两类样本天然规避了类别不平衡对单个分类器的影响。4.3 超参数网格与自动优化惩罚系数 BoxConstraint 控制间隔与误分类的权衡KernelScale 决定 RBF 核的作用半径。手工网格搜索是这样写的boxList [0.5 2 8 32]; % 对数尺度递增 sigList [0.5 1 2 4]; bestAcc 0; bestT []; for bc boxList for ks sigList t templateSVM(KernelFunction,rbf,BoxConstraint,bc, ... KernelScale,ks,Standardize,true); m fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone); L kfoldLoss(crossval(m, KFold, 5)); % 5 折交叉验证误分类率 if 1-L bestAcc bestAcc 1-L; bestT t; end end end16 组参数各跑 5 折样本量上千的时候耗时几分钟可以接受。如果 matlab优化工具箱可用直接交给贝叶斯优化更省事opts struct(Optimizer,bayesopt,ShowPlots,false, ... MaxObjectiveEvaluations,30, ... CVPartition, cvpartition(Ytr,KFold,5)); t templateSVM(KernelFunction,rbf,Standardize,true); mdl fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone, ... OptimizeHyperparameters,{BoxConstraint,KernelScale}, ... HyperparameterOptimizationOptions, opts);注意KernelScale和Standardize一起用时优化器搜的是标准化之后的数据尺度。如果跳过Standardize手动标准化KernelScale的最优值会完全不一样两者的搜索结果不能混着用。4.4 类别不平衡与评估指标现场数据里严重故障样本往往只有几十个轻故障几百个。两个手段一是templateSVM里设Prior,uniform让各类先验概率相等二是直接给少数类加权w ones(size(Ytr)); w(Ytr 3) 2.5; % 不平衡类样本权重抬高 t templateSVM(KernelFunction,rbf,BoxConstraint,8, ... KernelScale,auto,Standardize,true); mdl fitcecoc(Xtr, Ytr, Learners, t, Coding,onevsone, Weights, w);权重取多少可以用反频率比初估w_c N_total / (K * N_c)。补偿过度会让少数类误报率上升需要在混淆矩阵上反复看。评估别只看准确率。四类样本均衡时准确率够用一旦不均衡一个把所有样本都判成多数类的模型也能拿到 70% 准确率。Yp predict(mdl, Xte); C confusionmat(Yte, Yp); acc sum(diag(C)) / sum(C(:)); % 宏平均 F1先算每类的 P/R再等权平均不被类别规模带偏 K numel(unique(Yte)); f1 zeros(K,1); for c 1:K tp C(c,c); fp sum(C(:,c)) - tp; fn sum(C(c,:)) - tp; p tp/(tpfpeps); r tp/(tpfneps); f1(c) 2*p*r/(preps); end macroF1 mean(f1);如果混淆矩阵上看到内圈和外圈互相误判特别多说明载荷区调制带来的边带特征没被抓住——回到第 3 节把 2600–3200 Hz 共振带再细分成两个子带重新提特征通常比继续调 SVM 参数有效。5. GUI 集成与 R2025b 图形规范三个容易翻车的细节把上面几个函数串起来做成界面核心是参数传递和图形调用规范。用uifigureuigridlayout搭骨架控件回调写成嵌套函数共享父函数工作区里的模型句柄比用guidata存取清爽得多。训练这种耗时操作放进parfeval后台执行用afterEach回调刷新界面否则点一次「训练模型」窗口会卡十几秒。R2025b 有三个改动会直接让代码报错% 1) colormap 需要显式指定目标 figure colormap(fig, turbo); % 旧写法 colormap(turbo) 在多窗口下目标不明确 % 2) colorbar 的属性名变了别再设 ColorbarVisible cb colorbar(ax); cb.Label.String 功率 / dB; % 用 Label 对象而不是旧的可视性开关 % 3) confusionchart 不能当普通子级颜色要设给坐标区 cm confusionchart(ax, Yte, Yp); cm.Title 分类混淆矩阵;时频图的绘制顺序也有讲究先imagesc画功率谱再set(ax,YDir,normal)让频率轴从下往上递增最后调colormap和colorbar。顺序反了会出现颜色映射与坐标区不匹配的情况图看着没问题但刻度颜色对不上。坐标轴传的是uiaxes句柄而不是gca因为 GUI 里有多个子图时gca拿到的是当前焦点轴未必是你想画的那个。最后一个实用技巧实时预测场景下别每次重新提特征。把频带索引预先算好存成结构体流式数据进来后直接查表切片persistent bandIdx if isempty(bandIdx) [~, F] spectrogram(zeros(256,1), hamming(256), 192, 512, 12800, psd); edges [500 1200 2000 2600 3200 4000 5000]; bandIdx arrayfun((b) find(Fedges(b) Fedges(b1)), ... 1:numel(edges)-1, UniformOutput, false); end这样省掉的是每个样本都对 F 做一次逻辑比较的开销在 12.8 kHz 采样、每秒出一个诊断结果的在线场景里单样本特征提取从 38 ms 降到 24 ms 左右——对需要跑几十路通道的边缘盒子来说这个差距决定了能不能用单核跑完整个产线。本文还有配套的精品资源点击获取
