做滚动轴承故障诊断的同行应该都有过这种体验：故障特征明明存在，包络谱里却糊成一团，尤其是转速波动和背景噪声大的时候。这两年，稀疏表示在机械故障诊断里被频繁提起，常见的做法是拿小波字典或傅里叶字典做稀疏分解，但实际用下来，固定字典对轴承冲击信号的适应性并不理想。我去年在一个设备健康管理项目里换成了学习字典的思路，配合加权稀疏表示做信号重构，对早期故障的识别效果比普通稀疏重构稳定不少。这篇文章就把这套方案完整拆开来说：为什么要用学习字典而不是固定字典，加权稀疏表示的“权重”到底加在哪里、怎么加，以及可以直接跑的Python代码和调参时容易踩的坑。核心流程是字典学习、迭代重加权L1求解、包络谱分析三段，你把它拼到自己的数据上就能用。

1. 为什么要用稀疏表示做轴承故障诊断

1.1 轴承故障信号到底长什么样

滚动轴承发生局部损伤时（比如外圈剥落、内圈点蚀），滚子滚过缺陷位置会产生一个周期性的冲击。这个冲击本身不是一段简单的正弦波，它的频率成分很宽，会激励起轴承座和传感器安装位置的高频固有共振。所以我们实测到的故障信号，宏观上可以看成“周期性冲击 + 高频共振调制 + 背景噪声”，而且早期故障的冲击幅值非常小，信噪比经常是负的。

经典的包络谱分析（希尔伯特解调 + FFT）能解决一部分问题：它先带通滤波把高频共振分量保留下来，再通过包络解调把低频的冲击重复频率找出来。但它的前提是“带通滤波频带选择要合理”。频带选宽了，噪声进来太多；选窄了，冲击能量又被削掉。另一个麻烦是，冲击周期在转速波动下不是严格不变的，包络谱的谱线会被拉宽，特征频率处的幅值下降。这时候，稀疏表示的思路就体现出优势：它不再假设信号是“单频带”的，而是假设信号可以被少量“原子”线性组合出来，通过系数稀疏性把冲击成分从噪声里挑出来。

打个比方，包络谱像一把固定的梳子，梳齿间距是固定的，遇到稍微跑偏的冲击就梳不顺；稀疏表示则是先把头发里的“毛刺”都识别出来，再用少数几根主线把主要形态重构出来，对信号本身的畸变有更强的容忍度。

1.2 固定字典的局限与学习字典的直觉

早期做稀疏分解，大家习惯用傅里叶字典、小波字典、经验小波字典。这些字典有明确的数学定义，用起来也很方便。但它们有一个共同问题：原子形态是人为预设的。傅里叶原子是正弦波，适合稳态信号；小波原子虽然能刻画瞬态，但小波基函数、分解层数、小波参数都需要按经验去试。实际轴承冲击的形态受传递路径、转速波动、载荷变化影响，同一台设备在不同工况下冲击衰减波形都不一样，固定字典很难全部匹配。

学习字典则完全不同。它不预设原子形状，而是直接从振动数据里“学”出一组能最佳稀疏表示这批数据的原子。K-SVD是最经典的学习算法，它的思想是：先用OMP做稀疏编码，再用SVD逐个更新字典原子，让字典越来越贴近数据本身的结构。后面在线字典学习（比如sklearn里的DictionaryLearning）又做了优化，用随机块坐标下降，一次处理一批样本，速度和内存占用都更适合长信号。

我用下来最直观的感受是：学习字典训练完成后，很多原子确实长得很像“被共振调制的衰减震荡”，这就是机器自己总结出的冲击响应形态。相比小波字典，它在重构信号时不需要手工调“用哪个小波基”，模型的自适应性好很多。

1.3 为什么还要加一个“权重”

字典学习把字典准备好了，接下来就是稀疏编码。最正统的目标是求最小L0范数，即让非零系数尽量少，但L0问题是个NP难问题，实际都用L1范数做凸松弛。L1正则的毛病在于，它对所有系数一视同仁：无论你是来自真实冲击的大系数，还是来自噪声的小系数，惩罚力度都一样。结果就是稀疏重构出来的冲击幅值容易被压小，而噪声里的一些小峰又没有被剔除干净。

加权稀疏表示的动机就是打破这种“一视同仁”。我们给每个系数一个权重，大系数对应的权重小，允许它保留较大幅值；小系数对应的权重大，进一步把它压向零。用数学语言说，目标函数从 min ||x||₁ 变成 min ||Wx||₁，其中W是对角权重矩阵。权重的选择有很多种方案，工程上最常用的是迭代重加权L1：初始权重全为1，解完一次L1后，用当前系数幅值的倒数来更新权重，然后继续迭代。通常3到5轮之后，系数分布就稳定了，冲击对应的大系数被“保下来”，噪声对应的杂散小系数被“按下去”。

一个更生活的类比：普通的L1像是一把均匀口径的筛子，所有的沙砾都按同一个标准过滤；加权L1则是你观察完第一轮筛出的结果后，主动调小筛孔去过滤那些细沙，同时放宽对大颗粒的阻碍，几轮之后留在筛面上的就是真正有价值的“冲击块”。

2. 学习字典与加权稀疏表示的数学原理

2.1 稀疏表示的建模过程和字典学习目标

先设立统一的记号。假设我们从振动信号里截取一段长度为m的帧，用向量y∈R^m表示。我们希望找到一个过完备字典D∈R^(m×n)（n > m，即原子数量大于信号维度），使得y≈Dx，并且系数向量x尽量稀疏。加上噪声后，基本模型是：

y = Dx + e，e为噪声项。

稀疏编码的优化目标可以写成：

min ||x||_0 subject to ||y - Dx||_2^2 ≤ ε

由于L0范数不可处理，工程上替换为：

min 0.5 * ||y - Dx||_2^2 + λ * ||x||_1

这里λ是正则化参数，控制拟合误差和稀疏性之间的平衡。λ越大，系数越稀疏，但重构误差也越大；λ太小，稀疏性失效，退化成最小二乘。

而字典学习就是同时求字典D和稀疏系数X，使得整个训练样本集Y能被稀疏表示：

min 0.5 * ||Y - DX||_F^2 + λ * Σ||x_i||_1，附带对字典原子的归一化约束。

K-SVD的思路分两步交替迭代：第一步固定字典，用OMP或基追踪求每个样本的稀疏系数；第二步固定系数矩阵，逐个原子做SVD更新——每次只更新一个原子和它对应的系数行，这样可以显著降低计算复杂度。在线字典学习则是每次抽一个小批量样本，用一个可学习的步长去更新字典，更适合信号很长、样本量很大的场景。

2.2 加权L1的数学形式和变量替换技巧

加上权重矩阵W之后，稀疏编码目标变成：

min 0.5 * ||y - Dx||_2^2 + λ * ||Wx||_1

其中W是对角矩阵，第i个对角元素w_i表示对系数x_i的惩罚权重。如何高效求解这个问题？我常用一个简单的变量替换技巧：设z_i = w_i * x_i（等价于z = Wx），则Dx可以写成：

Dx = Σ d_i * x_i = Σ (d_i / w_i) * z_i

也就是说，如果把字典D的每一列除以对应的权重w_i，得到D_w，那么原问题变成标准Lasso形式：

min 0.5 * ||y - D_w * z||_2^2 + λ * ||z||_1

求解出z之后，反变换x_i = z_i / w_i就得到原始系数。sklearn里的Lasso、坐标下降法、近端梯度法都可以直接用来求解这个标准Lasso问题。这个技巧很实用，不用额外去实现加权L1专用求解器。

迭代重加权L1的完整流程如下：

1. 初始化w_i = 1（所有i）
2. 按上面的变量替换求解标准Lasso，得到x
3. 更新权重：w_i = 1 / (|x_i| + ε)，ε是一个很小的正数，防止除零
4. 重复步骤2-3，直到权重变化较小或达到预设迭代次数

我实践下来，迭代次数设3-5次就够。多次迭代的目的不是追求数学上的收敛到全局最优，而是让权重能根据第一次稀疏重构的结果做一次“筛选”——把真正可能属于冲击的大系数暴露出来，把疑似噪声的小系数压下去。迭代次数超过6次之后收益很小，反而可能因为权重过于极端，把某些真实冲击系数也压没了。

2.3 权重设计与实际效果的关系

权重设计的核心矛盾在于“容忍大系数”和“抑制小系数”之间如何取舍。如果w更新得太激进，比如ε取得太小，权重动态范围会非常大，可能导致某些原子虽然匹配冲击、但第一次迭代时系数没有足够大，后续权重被压得极大，结果该原子被彻底冻结，重构信号丢失部分冲击细节。我的经验是把ε取在1e-6到1e-4之间，并且对权重做裁剪：

w = np.clip(1.0 / (np.abs(x) + 1e-6), 1e-3, 1e4)

上限1e4防止权重爆炸，下限1e-3防止所有权重都变成零。实际调试时，如果重构信号出现“脉冲丢失”，优先怀疑ε太小或λ太大；如果重构信号留下大量毛刺，则说明λ太小，稀疏性不够。

3. 可直接运行的代码与诊断流程

3.1 环境准备与仿真信号生成

我用Python实现整套流程，依赖numpy、scipy、scikit-learn和matplotlib。安装很简单：

pip install numpy scipy scikit-learn matplotlib

先构造一个外圈故障的仿真信号，方便复现和调试。采样率取12000 Hz，转频29.5 Hz，外圈故障特征频率BPFO约107.3 Hz，冲击激起的共振频率1600 Hz，衰减系数0.02。这个仿真信号的核心是用周期冲击激励一个二阶衰减振荡，再叠加高斯白噪声。

import numpy as np from scipy import signal from scipy.fft import fft, fftfreq from scipy.signal import hilbert from sklearn.decomposition import DictionaryLearning from sklearn.linear_model import Lasso import matplotlib.pyplot as plt fs = 12000 # 采样率 fr = 29.5 # 转频 fn = 1600.0 # 共振频率 BPFO = 107.3 # 外圈故障特征频率 zeta = 0.02 # 衰减阻尼比 T = 4.0 # 信号时长 N = int(fs * T) t = np.arange(N) / fs np.random.seed(42) # 冲击时刻，带微小随机抖动模拟转速波动 impulse_times = np.arange(0, T, 1.0 / BPFO) impulse_times = impulse_times + np.random.normal(0, 0.2 / BPFO, size=impulse_times.shape) x = np.zeros(N) A = 1.0 for ti in impulse_times: if ti >= T: continue idx = np.where((t >= ti) & (t < ti + 0.03))[0] rel = t[idx] - ti x[idx] += A * np.exp(-zeta * 2 * np.pi * fn * rel) * np.sin(2 * np.pi * fn * rel) # 加高斯白噪声 noise_std = 0.15 x_noisy = x + noise_std * np.random.randn(N)

这里冲击时刻加了随机抖动，是为了模拟真实工况下的转速波动和滑动。故障特征频率在包络谱上不会是一条理想的谱线，而是有一定宽度的峰，更接近真实情况。

3.2 滑动分帧与字典训练

训练学习字典之前，需要把原始振动信号切成一帧一帧的样本。窗口长度选256，覆盖两个左右的冲击周期（冲击周期约112个采样点），这样每帧既能包含完整的冲击波形，又不会把信号切得太碎。重叠率取50%，也就是步长128，保证相邻帧的连续性。

def sliding_windows(sig, win=256, step=128): n = len(sig) frames = [] starts = [] for s in range(0, n - win + 1, step): frames.append(sig[s:s + win]) starts.append(s) return np.array(frames), np.array(starts) Y, idxs = sliding_windows(x_noisy, 256, 128) # 随机抽取一部分帧训练字典 np.random.seed(0) train_idx = np.random.choice(len(Y), min(2000, len(Y)), replace=False) Y_train = Y[train_idx] dl = DictionaryLearning( n_components=128, alpha=0.5, max_iter=50, transform_algorithm='lasso_lars', fit_algorithm='lars', random_state=0 ) dl.fit(Y_train) D = dl.components_.T # 转置成 m x n_atoms，每一列是一个原子 D = D / np.linalg.norm(D, axis=0, keepdims=True) # 原子归一化

DictionaryLearning里原子数量和alpha这两个参数最直观。原子数量太少，字典表达能力不够；原子数量太多，训练时间长且容易过拟合噪声。先按128起步，后续再看重构效果调整。

3.3 加权稀疏编码与信号重构

字典训练好之后，对每一帧信号做加权稀疏编码。我写了一个可复用的函数，内部用迭代重加权L1加变量替换求解：

def weighted_l1_solve(y, D, lam=0.02, n_iter=5, eps=1e-6): n_atoms = D.shape[1] w = np.ones(n_atoms) x = np.zeros(n_atoms) for _ in range(n_iter): Dw = D / w[None, :] # 每一列除以权重 lasso = Lasso(alpha=lam, fit_intercept=False, max_iter=2000) lasso.fit(Dw, y) z = lasso.coef_ x = z / w w = np.clip(1.0 / (np.abs(x) + eps), 1e-3, 1e4) return x

这里的核心逻辑在前面已经解释过：设z = Wx，则字典列需要除以权重。第一次迭代w全为1，等价于普通L1；从第二次开始，权重会依据前一次系数调整。调用这个函数对所有帧重构：

recon = np.zeros(N) weight_sum = np.zeros(N) win = 256 for s in idxs: seg = x_noisy[s:s + win] x_sparse = weighted_l1_solve(seg, D, lam=0.02) seg_rec = D @ x_sparse recon[s:s + win] += seg_rec weight_sum[s:s + win] += 1 recon /= weight_sum # 重叠区取平均

实际信号有4秒，滑窗总数约370帧，每帧内部要跑5次Lasso，所以整段会跑几秒到十几秒，属于正常的计算量。工程上如果不做在线处理，这么几秒完全可以接受；如果信号长达几十秒，建议先做降采样或只对疑似故障段做重构。

3.4 希尔伯特包络谱分析与故障判定

重构完成后，用希尔伯特变换求包络，再对包络做FFT，得到包络谱。特征频率处出现明显峰值，是诊断外圈故障的关键依据。

analytic_sig = hilbert(recon) env = np.abs(analytic_sig) Env_spectrum = np.abs(fft(env))[:N // 2] freqs = fftfreq(N, 1 / fs)[:N // 2] # 限制只看到500Hz以内，对特征频率的倍频看得更清楚 mask = freqs < 500 plt.figure(figsize=(10, 4)) plt.plot(freqs[mask], Env_spectrum[mask]) plt.axvline(BPFO, color='r', linestyle='--', label='BPFO') plt.axvline(2 * BPFO, color='g', linestyle='--', label='2xBPFO') plt.xlabel('Frequency / Hz') plt.ylabel('Amplitude') plt.legend() plt.show()

诊断判据通常有三条：一是BPFO处有主导峰；二是BPFO的2倍频、3倍频谐波清晰可见；三是边频带结构相对干净。如果是内圈故障，BPFI处会有峰值，并且由于内圈随轴旋转，故障冲击的幅值受到转频调制，所以BPFI附近还会有fr间隔的边频带。滚动体故障的BSF则更隐蔽，因为冲击位置在滚动体和内外圈之间变动，能量更分散，包络谱的峰往往比较弱，这时候稀疏重构的优势会更明显——它能从噪声里把微弱的周期冲击“捞”出来。

为了对比普通L1和加权L1，可以把上面代码中n_iter设成1，也就是只跑一次标准Lasso，再算包络谱。我实测的效果是：标准L1在BPFO处的峰幅值明显偏低，谐波几乎看不出第3倍频；加权L1在BPFO处的幅值更高，2倍频和3倍频都能清晰分辨，噪声底也低了一截。特别是信噪比降到5dB以下时，标准L1的重构结果还会出现冲击幅值被严重压扁的问题，故障特征频率容易淹没在噪声底中，加权L1则好很多。

4. 参数调优、常见问题与工程经验

4.1 关键参数速查表

这套方法涉及参数比较多，我整理成了一份速查表，按经验给出发起值和调整方向。

一个比较隐蔽的坑是：lambda和字典训练alpha是有相关性的。字典训练时用了L1稀疏约束，得到的原子“整体形态”偏稀疏；测试阶段再用大的lambda，会过度压缩系数，导致重构信号能量偏低。我习惯先固定训练alpha=0.5，然后测试阶段在0.01到0.03之间微调lambda，这样参数空间比较小。

4.2 特征频率算不准怎么办

最让我头疼的不是稀疏重构本身，而是轴承参数不完全准确的时候，特征频率理论值算出来和实际峰值对不上。比如轴承型号标注的接触角、滚动体直径是名义值，实际安装后会有偏差；转速也不是恒定值，直接拿名义转速去算BPFO，偏差可能超过2%。这时候如果只盯着理论值找峰，很可能误判。

我的排查流程是：先按名义参数算出理论特征频率，然后在理论值附近±3%范围内找局部峰值，同时看它的2倍频、3倍频是否也存在。如果只有基频没有谐波，大概率是某个随机冲击的干扰峰；如果基频、2倍频、3倍频都有，并且间隔一致，基本可以认定为同类故障特征。对于内圈故障，还要观察特征频率周围有没有转频间隔的边频带。边频带出现得越密集，说明调制越严重。为了快速检验，可以在代码里写一个简单的谐波匹配函数：对目标频率的1~4倍频窗口内取局部最大幅值，累加作为“故障证据分”，分数超过阈值再报警，这样可以减少人眼判图的负担。

4.3 常见问题与调试记录

我在实际项目中遇到过几类典型问题，这里直接按“现象—原因—办法”列出来。

第一类：包络谱里BPFO处没有峰，但原始时域波形里能隐约看到周期冲击。原因是字典学到的原子趋于“小波状”，而测试帧的冲击相位和原子形状没有对齐，导致稀疏系数分散。解决办法是把窗口长度调小到覆盖1~2个冲击周期，让原子更像单个冲击，而不是一小段多个冲击；同时增加重叠率到75%，让冲击在不同帧里都有相对完整的出现机会。

第二类：重构信号出现很多伪冲击，波形看起来毛刺很多。这往往是lambda太小，稀疏性不够，噪声也被当作有效成分重构出来。先把lambda调到0.05量级，如果毛刺消失但冲击幅值也变小，再回调一点。调节过程中可以参考重构信号与原始信号的相关系数：相关系数太高说明可能连噪声一起重构了，太低说明丢失了有效成分。

第三类：不同段信号跑出来结果波动很大。原因往往是字典训练样本随机抽取导致的原子不固定。解决办法是固定随机种子（代码里的random_state和np.random.seed），并且训练样本量尽量超过2000帧。样本太少时，字典会偏向某一小段信号的形态，泛化能力差。

第四类：加权迭代后期把某些真实冲击系数压没了。这是权重更新过于激进的典型表现。我会检查更新后的权重分布，如果很大一部分权重接近上限1e4，说明epsilon太小。把epsilon从1e-6改到1e-4，权重分布会温和很多，重构出的冲击也更完整。

4.4 在线部署时的效率优化与个人经验

如果要把这套流程放到在线系统里，我建议分两步走。第一步离线训练字典，输入不同类型的历史故障数据，学出一个通用字典；第二步在线只做稀疏编码和包络谱分析，不重复训练字典。在线编码一个256点的帧，跑5次Lasso大约需要几十毫秒，基本可以满足每秒几十个帧的处理速度。如果实时性要求更高，可以用近端梯度法或者ISTA替代sklearn的Lasso，速度会再快一些，代码也完全可控。

我个人在实际操作中还有一个体会：稀疏加权重构和传统包络谱不是替代关系，而是互补关系。传统包络谱速度快、物理含义明确，适合作为在线监控的主通道；稀疏加权重构计算量大，但对早期微弱故障更敏感，适合在报警后作为“确认手段”或者离线故障复判工具。先快后准，这套组合在项目里比单独使用任何一种方法都稳。做完一次完整的轴承故障诊断流程后，你再看那些“包络谱一塌糊涂”的数据，就不会急着去调滤波器参数了——把信号交给学会了冲击形态的字典，它往往能给出意想不到的清晰结果。