滚动轴承出问题的时候整条产线都不会等你。故障从“轻微”到“报废”可能就隔着一个检修周期而这窗口期里振动信号里的故障特征往往弱到被背景噪声完全淹没。我做过一段时间的设备故障诊断最头疼的就是理论上都知道内圈外圈的故障频率是算得出来的但拿实际振动数据一跑特征峰就是出不来甚至被强干扰吃掉。这个项目“基于学习字典与加权稀疏表示的滚动轴承稀疏表示故障诊断”要解决的就是这个问题。核心思路一句话不靠固定的小波基或者频谱分析硬套而让字典从数据本身“学”出故障冲击的形态再用加权策略把噪声碎片压下去最后从包络谱里把故障特征峰顶出来。它不是一个纯理论推导演示而是带完整可运行代码的工程级方案采样率、字典参数、特征频率计算、阈值判断都给你铺好了适合做设备状态监测的工程师、机械故障诊断方向的研究生以及所有对稀疏表示好奇但不知道怎么落地的信号处理学习者。1. 方案的底层逻辑为什么是“学习字典”而不是小波包1.1 稀疏表示先想清楚“稀疏”到底在求什么稀疏表示的本质很简单一段信号里能起决定作用的分量其实很少剩下的大部分是冗余、噪声、或者与目标无关的干扰。就像一锅汤真正决定味道的只有几片肉其他都是清水。数学上我们找一个字典矩阵D它由一堆原子列向量构成然后把信号x拆成x ≈ Dα其中系数向量α大部分是零或者接近零的“稀疏”结构。对这个项目来说滚动轴承故障时产生的周期性冲击在时域上是很有特征的瞬态响应——衰减振荡的一串小尖峰。这种信号的特点是“少数几个冲击事件占主导”与稀疏表示的前提天然吻合。理论上只要字典D里存在与冲击形态相似的原子噪声再大稀疏求解也倾向于用几个大系数去逼近冲击而不是用一堆小系数去拟合噪声。但这里有个关键点如果我们不好好设计字典稀疏表示的效果会大打折扣甚至直接失效。这也是为什么众多论文里对“固定基字典”的质疑这么深。1.2 为什么非要用“学习字典”而不是固定的傅里叶/小波字典常用的固定字典有很多傅里叶基、离散余弦变换DCT基、小波基、曲波基等等。它们的好处是数学性质明确算起来快坏处是“一刀切”。轴承故障冲击的形态会随着故障尺寸、转速高低、负载大小、传递路径远近而变。同一台设备、同一个轴承外圈早期点蚀和中后期剥落的冲击波形差异很大换一台设备固有频率不同冲击的振荡频率和衰减快慢又是另一个样。拿一套固定写死的小波字典去套所有工况大概率出现“匹配不精、系数不稀疏”的窘境。学习字典的思路是数据驱动的把运行状态的实测信号切成一段段样本用K-SVD这类算法去学出一组“针对这台设备”的原子。K-SVD这个名字念出来有点吓人但它的思路其实通俗K代表字典里有K个原子SVD则是原子更新时用的矩阵分解工具。整个流程就是交替做两件事固定字典求稀疏系数再固定系数逐列更新字典原子循环往复直到字典能更好地表达训练信号。实践下来学习字典最大的好处是学到的原子里往往真的会出现“长得像故障冲击”的原子一段快速衰减的类指数振荡波形紧接着一小段静默或者是周期性间隔的低频调制形态。这种原子是我们在预定义字典里很难手动构造出来的而它一旦出现后续的稀疏编码就会非常有指向性故障特征的“可分离度”直接上一个台阶。1.3 “加权”到底加在哪里又是怎么起作用的很多人看到“加权稀疏表示”会困惑觉得稀疏表示本身已经带了约束为什么还要加权。我的理解是普通稀疏表示是“一视同仁地追求系数个数少”但实际信号里不同原子对故障诊断的价值是不一样的。有些系数是代表故障冲击的我们想保住它有些系数是噪声偶然碰到的“伪匹配”我们想压低它。这就需要根据系数本身的强度或者原子与待检信号的关联度给不同的分量赋予不同权重。项目里我用的是迭代重加权思路先用一次普通的OMP稀疏编码得到初始系数然后根据系数绝对值构造权重——系数大的地方说明是真信号权重给得小一点约束放松些让真冲击能保留系数小的地方多半是噪声碎片权重给大一点进一步把它往零压。如此迭代两三轮相当于给每个原子“按贡献分蛋糕”而不是稀里糊涂地平均分配稀疏预算。这个操作在强噪声场景的效果相当明显不加权时包络谱里总有若干假峰晃动加权之后谱线干净很多。2. 整体流程设计与关键参数2.1 从振动信号到诊断结论的完整链路整个项目跑通的路子是一条非常典型的“降噪-重构-解调-判决”链路。我先把它完整列出来后面每步再拆开讲读取振动信号按采样率和工况记录去除趋势项必要时做带通滤波。对信号做滑窗分帧每帧长度取512或1024点相邻帧重叠50%~75%。对每帧做去均值与方差归一化然后送入字典学习阶段。用K-SVD或MOD从训练帧中学习字典原子数通常取128~256。对目标信号分段执行加权稀疏编码得到稀疏系数。用字典原子与稀疏系数重构“去噪冲击信号”。对重构信号做Hilbert变换求包络再做FFT得到包络谱。在包络谱里寻找轴承特征频率及其谐波处的谱峰与理论值对比后给出诊断结论。这套流程里最容易犯的错是把前面的分帧和归一化当“杂活”随手做。实际上字典学习很敏感如果信号里带着直流偏置字典学出来的第一个原子可能不是冲击而是“一条直线”如果各帧振幅差异太大字典会把训练集中振幅最大的几帧反复拟合其他帧的信息被无视。所以预处理我建议一律先中心化再除以帧内均方根让字典专注学“形态”而不是“幅值”。2.2 核心参数定多少心里要有谱我给出的是经过大量实验验证的经验区间同时也是自查清单参数推荐范围说明窗口长度 N512~1024 点必须至少覆盖1~3个故障冲击间隔太短学不到完整冲击太长冲击被淹没字典原子数 K128~256满足K远大于稀疏度即可K过大会让训练变慢且出现大量冗余原子稀疏度 L4~8太小漏掉冲击太大会连噪声一起编进重构信号字典学习迭代次数10~20太少字典没学成型太多对有限样本容易过拟合分帧重叠率50%~75%提高数据使用效率减少重构时帧边界断层加权迭代轮数2~3更多轮收益很有限基本收敛为什么窗口长度建议至少1~3个冲击周期因为稀疏表示虽然不强制信号严格周期但我们后续要靠包络谱判断故障频率如果窗内连两个完整冲击都装不下重构出来的冲击序列在时域上就没有周期性包络谱里自然看不见清晰的故障频率峰。稀疏度这块我得啰嗦一句它不是越大越好。稀疏度太大相当于允许字典用很多原子去逐点逼近信号结果噪声的“细节”也被当成有效成分重构出来包络谱虽然还是有故障峰但旁瓣和杂峰会明显增多容易误判。我在实测里常见的情况是同样一段信号L5时谱峰干干净净L15时理论频率处确实有峰但旁边跟着三四个与转频相关的边带反而干扰判断。2.3 先算好理论特征频率再拿结果比对整个诊断流程离不开“靶子”轴承故障特征频率。以滚动轴承最常见的深沟球轴承为例几个典型特征频率的计算关系简化形式如下外圈故障频率BPFO ≈ n × fr × (1 - d/D × cosα) / 2内圈故障频率BPFI ≈ n × fr × (1 d/D × cosα) / 2滚动体故障频率BSF ≈ fr × D/d × (1 - (d/D × cosα)^2)其中n是滚动体数目fr是转频d是滚动体直径D是节径α是接触角一般工程估算时α可以按20度左右带入。举个例子某型轴承滚动体数12节径65mm滚动体直径15mm轴转速1800rpm转频30Hz。算下来外圈BPFO约为0.5 × 12 × 30 × (1 - 15/65) ≈ 138.5Hz内圈BPFI约为0.5 × 12 × 30 × (1 15/65) ≈ 221.5Hz。有了这个数值在包络谱里直接看这两个频率位置有没有峰峰高相对本底噪声是否有3到5倍以上的优势是否有2阶甚至3阶谐波就可以给诊断结论了。特别提醒理论频率是按当前实际转速算的不是按电机铭牌上的额定转速。变频驱动的设备转速可能随时在飘拿额定值硬套结果必然对不上谱峰。所以现场用的时候务必先通过转频峰或者键相设备获取真实转速再代入公式。3. 代码实现与踩坑实录3.1 工程框架一个类把诊断流程串起来这个项目我把整个诊断流程封装成了一个类核心方法包括分帧、字典学习、稀疏编码、加权重构、包络谱分析和诊断判决。类的好处是便于复用你拿到新设备数据改一下初始化参数和特征频率就能直接跑。下面贴出核心骨架代码用Python写依赖numpy、scipy和scikit-learn其中字典学习我直接用了scikit-learn的DictionaryLearning接口省得自己手搓K-SVD的细节。import numpy as np from scipy.signal import hilbert from sklearn.decomposition import DictionaryLearning def frame_signal(signal, win_len, overlap): step int(win_len * (1 - overlap)) n_frames (len(signal) - win_len) // step 1 frames np.stack([signal[i * step : i * step win_len] for i in range(n_frames)]) # 去均值 RMS归一化消除直流和幅值差异 frames frames - frames.mean(axis1, keepdimsTrue) frames frames / (np.linalg.norm(frames, axis1, keepdimsTrue) 1e-12) return frames def omp(block, D, max_atoms, residual_threshold1e-5): 正交匹配追踪逐次挑选与残差最相关的原子 residual block.copy() n_atoms D.shape[1] coef np.zeros(n_atoms) support [] for _ in range(max_atoms): proj np.abs(D.T residual) idx np.argmax(proj) if proj[idx] residual_threshold: break support.append(idx) Ds D[:, support] coef_s np.linalg.lstsq(Ds, block, rcondNone)[0] residual block - Ds coef_s coef[support] coef_s return coef def weighted_ista(y, D, W, lamb0.001, max_iter80): 迭代软阈值算法求解加权L1稀疏表示 L np.linalg.norm(D, 2) ** 2 DtD D.T D Dty D.T y alpha np.zeros(D.shape[1]) for _ in range(max_iter): grad DtD alpha - Dty alpha alpha - grad / L alpha np.sign(alpha) * np.maximum( np.abs(alpha) - lamb * W / L, 0) return alpha class SparseBearingDiagnoser: def __init__(self, n_atoms128, sparse_level5, win_len512, overlap0.5, fs12000): self.n_atoms n_atoms self.sparse_level sparse_level self.win_len win_len self.overlap overlap self.fs fs self.D None def train_dictionary(self, signal, max_iter15): frames frame_signal(signal, self.win_len, self.overlap) dl DictionaryLearning( n_componentsself.n_atoms, max_itermax_iter, transform_algorithmomp, transform_n_nonzero_coefsself.sparse_level, positive_codeFalse, random_state42 ) self.D dl.fit(frames).components_.T return self.D def reconstruct(self, signal): frames frame_signal(signal, self.win_len, self.overlap) rec_frames np.zeros_like(frames) for i, frame in enumerate(frames): coef omp(frame, self.D, self.sparse_level) # 用系数绝对值做权重强项保留、弱项收缩迭代两次 for _ in range(2): W 1.0 / (np.abs(coef) 1e-6) coef weighted_ista(frame, self.D, W) rec_frames[i] self.D coef # 重叠相加还原为完整信号 step int(self.win_len * (1 - self.overlap)) rec np.zeros(len(signal)) count np.zeros(len(signal)) for i, frame in enumerate(rec_frames): s i * step rec[s : s self.win_len] frame count[s : s self.win_len] 1 return rec / (count 1e-12) def envelope_spectrum(self, signal): analytic hilbert(signal) env np.abs(analytic) spec np.abs(np.fft.rfft(env)) ** 2 freq np.fft.rfftfreq(len(signal), 1.0 / self.fs) return freq, spec这段代码的注释栏已经把关键逻辑标出来了。有一点要特别说明weighted_ista里的lamb参数并没有做复杂的自适应实测下来0.001到0.01这个量级对归一化后的信号基本都能用。如果你换的数据幅值差异很大记得先对信号做全局归一化否则lamb的标定会失效。3.2 用仿真信号先验证流程是否走得通跑新方法之前我强烈建议先用仿真信号打通全流程别一上来就丢真实数据。为什么因为真实现场信号复杂影响因素多算法如果效果不好你分不清是“算法本身有问题”还是“数据预处理没做干净”。仿真信号的好处是我们提前知道故障频率在哪里答案在手再来看算法的输出对不对得上这就把变量控制住了。我用这样的方式构造了一段“外圈故障”仿真信号采样率12000Hz转频30HzBPFO约138.5Hz每个冲击是衰减振荡波形冲击重复周期约7.2ms再叠加高斯白噪声信噪比压到5dB左右。这样做的效果是原始波形基本看不清冲击包络谱直接做也是一团噪声但经过字典学习和加权稀疏重构之后包络谱里在138~139Hz处会顶出明显的峰并且能看到二次谐波。验证逻辑很简单重构后包络谱峰值频率与理论BPFO的偏差小于1Hz且谱峰超过底部噪声3倍说明流程自洽。这个步骤也帮我筛掉了不少参数问题第一次设置窗口长度256点时冲击间隔约86点一个窗口只能装下三个冲击字典学出来的原子几乎全是噪声形态把窗口加到512点后原子形态才像样。后面我调试任何新数据都先跑一遍仿真确认参数没问题再换数据。3.3 实操中最容易翻车的5个地方这几点都是我实际摸爬滚打总结出来的文档里一般不会写。第一分帧不做去均值和归一化字典学出来的原子是直流和幅值“老大”的形状而不是冲击形状。标准做法是帧内减均值、除以帧范数或RMS。第二原子数拍脑袋设512训练时间翻了几倍效果却和128基本没区别。学习字典的原子数应略大于稀疏度即可128到256点已经覆盖绝大多数轴承诊断场景再大纯属浪费算力。第三OMP的停止条件只设稀疏度不设残差阈值。低信噪比时允许的原子数还没用完残差已经全是噪声了算法照样继续匹配结果那几颗“额外原子”完全是在拟合噪声。建议同时设残差阈值残差降到一定水平就直接停。第四Hilbert包络谱没做平方结果峰是出来了但旁边多了明显的调制边频带很容易和滚动体故障的特征混淆。我建议取平方包络谱也就是先求包络再平方再FFT这样特征频率处的能量更集中。第五理论特征频率拿计算器按完就算忘了和实际转速核对。现场设备转速和额定转速差个5%非常常见138Hz的理论峰实际飘到145Hz你死盯着138Hz自然看不到东西。诊断前务必通过转频峰确认实际fr。4. 常见问题排查从“算法能跑”到“诊断可信”4.1 包络谱里就是看不到特征峰先查这四步这是被问得最多的问题代码跑通了也没有报错但是重构后包络谱里找不到理论上的BPFO或BPFI峰。我的排查顺序是这样固定的先查理论频率本身有没有算错。接触角、滚动体直径、节径、滚动体数目这些几何参数少一个都对不上。有的工程师直接拿轴承型号手册里的参数代入却忽略了手册给出的是节径还是外径这一步最容易差出百分之十几。再查采样率和窗长。采样率太低导致冲击振荡的高频成分已经被混叠字典学出来的原子形态是畸形的。看窗长是不是覆盖不了两个冲击周期窗口太短会让字典把“冲击一段噪声”学成一个原子重构出来的信号周期性很差。然后查预处理。原始信号里的转频成分和齿轮啮合频率如果能量太强会把故障冲击“压”下去字典学出来的是轴频弦波形态。这种情况我建议先做谱峭度或者带通滤波粗略锁定故障共振频带后再进字典。最后才怀疑算法参数稀疏度太小、迭代次数不足、字典初始化不合适。4.2 重构后的信号像低通滤波冲击被磨平了这个问题出过好几回现象是重构信号波形很光滑衰减振荡的尖锐冲击变成了圆润的小鼓包。原因有两个一是字典里缺少“高频瞬态型”的原子K-SVD学到的基本是低频慢变形态二是稀疏度开太小一个冲击就需要两三个原子配合才能拟合出快速衰减形态给少了凑不出来。解决办法第一换更长的窗口或更多训练样本让字典见过足够多的高频冲击形态第二把稀疏度上调一档从4调到6或8第三检查OMP的残差阈值设置残差阈值太紧匹配被迫提前终止高频细节就丢了。4.3 变工况、大数据量这套方案还顶得住吗很多设备的转速不是恒定的特别是风机和泵类负载转速会缓慢漂移。这时候固定时间域的稀疏表示效果会打折因为等时间间隔采样下冲击之间的间隔是变化的字典学的“周期冲击原子”就不够准了。工业上的通行做法是先做角域重采样用键相或编码器信号把等时间间隔的振动数据重采样到等角度间隔把变转速问题变成准恒定转速问题再做稀疏表示。如果设备没有键相也可以用转速计信号离线估算瞬时转速再做重采样。这一步做完前面整套流程不用大改。大数据量方面振动数据动辄几百万点逐帧OMP加ISTA确实不算快。实测下来128原子的字典、512点窗口、75%重叠处理一分钟的12k采样数据大约需要几十秒。在线监测场景如果嫌慢可以把重叠率降到50%或者先用谱峭度把关注频带缩小只对带通滤波后的信号做稀疏计算量能降一个量级。4.4 问题速查表现象可能原因解决路径包络谱无理论特征峰理论频率算错/转速不准核实轴承几何参数实时测量转速重构信号“过平滑”字典缺少瞬态原子/稀疏度过低增加训练样本上调L至6~8谱峰有大量边带平方包络没做/包络中有调制改用平方包络谱抑制低速调制边带字典原子全是噪声形态分帧没归一化/窗长过短加去均值与RMS归一化加窗长不同工况结果不稳转速波动/负载变化做角域重采样或分工况单独建字典计算太慢字典过大/重叠率过高减小K降重叠率先带通再稀疏5. 从实验到现场这套方案还能怎么演进5.1 与包络谱、谱峭度的组合打法单靠稀疏表示包打天下是不现实的工程上我更推荐组合拳先用谱峭度或者快速峭度图在原始信号里找到包含故障冲击的共振频带做带通滤波把明显的无关干扰先剔掉再用学习字典对滤波后的信号做稀疏重构进一步把残留噪声压掉最后做平方包络谱在窄带里找特征频率。每一步都在做“减法”最后留给判决的就非常干净了。我在现场项目里实测过单纯谱峭度包络谱在信噪比低时容易出假峰但加上稀疏重构后误报明显减少。这套组合的另一个好处是大幅降低字典学习的压力。带通滤波之后的信号形态相对简单字典可以用更少的原子学得更准速度也跟着上来了。5.2 在线诊断与预警阈值的设定离线分析跑通之后很多人的下一步就是做成在线监测。我的建议是不要拿单帧结果直接报警容易误触发。更稳的做法是滑窗连续诊断把包络谱里特征频率处的幅值和时间做趋势记录只有连续多帧比如5帧以上都超过本底噪声3倍才触发预警。阈值不是拍脑袋定的最好用正常运行数据先跑一遍统计出特征频率处的本底波动范围再乘以安全系数当作报警门槛。5.3 我个人的使用习惯与换参建议做了这么多轮实验我自己的默认参数基本固定为原子数K128稀疏度L5窗口长度512重叠率75%字典迭代15次加权轮数2遍。这套设置在不同采样率和不同转速下都表现出很好的稳定性。新手一开始不需要到处调参先把这套默认值跑通再根据自己设备的信号特点微调。换到新数据源时我的建议永远是先仿真后实测先带通后稀疏先单段后连续。严格按照这个次序来能少走一半弯路。这套方案的潜力还不止于轴承齿轮箱的局部断齿、往复机的阀体漏气只要故障形态是周期性瞬态冲击本质上都可以套这套“学习字典加权稀疏包络谱”的框架只是特征频率的计算公式要跟着传动链重新推一遍。诊断方法的价值不在于代码本身跑得多快而在于它对每一种真实故障形态的“理解”有多深而学习字典恰恰给了我们这种理解能力。