说实话做了几年振动信号分析我在轴承故障诊断这件事上有过一个比较强烈的感受常规小波分析不是不能用而是总觉得差那么一口气。明明频谱图上转频、倍频都在可故障特征频率就是不够干净提取得模棱两可。直到我把思路从“拿现成小波去套信号”转成“给信号专门设计一种脉冲小波原子再用稀疏表示去匹配”整个流程才豁然开朗。这篇文章就以此为切入点完整记录我在MATLAB里如何基于模型的脉冲小波model-based impulsive wavelet构造稀疏字典、编写OMP分解程序并且把它应用到轴承故障诊断中的过程。内容包含原理拆解、代码实现、仿真验证、实测数据推断以及我在这个过程中踩过的几个值得注意的坑。适合正在做故障诊断研究的学生也适合刚接触稀疏表示、想在一个具体场景里落地这项技术的工程师参考。1. 轴承故障脉冲为什么选通用小波总差一口气信号形态与基函数错位1.1 轴承故障信号的真实长相不是正弦波而是一连串“敲击声”滚动轴承出现局部损伤外圈点蚀、内圈剥落、滚动体缺陷时滚动体每经过损伤位置就会产生一次冲击。这种冲击的持续时间极短通常只有几毫秒甚至更短但它会立刻激起轴承座和传感器安装系统的固有振动。所以传感器拾取到的信号数学上看是这样一个东西在某个时刻 t τ 处突然出现一个幅值跳变随后信号以系统固有频率 f_n 做振荡同时振荡幅度按指数规律衰减大约经过几个周期就被阻尼吃掉。然后隔一段时间下一次冲击又来了。写成表达式就是典型的“阻尼正弦振荡”形态其有效部分只在冲击发生后的很短窗口内存在其余时间几乎为零。这个物理特征决定了它与连续平稳的正弦波、甚至与常见的连续小波基函数有本质区别。1.2 通用小波库的错位在哪里尺度匹配不等于形态匹配很多人一开始会选择Daubechies小波、Symlets小波或者Morlet小波去做连续小波变换直觉是“小波擅长捕捉突变信号”。这话没错小波变换确实有良好的时频局部化能力。但问题出在基函数形态上Daubechies小波虽然有紧支撑但它没有内在的“阻尼振荡”结构用它去分解一段脉冲信号需要非常多层的细节系数才能凑出那一下冲击Morlet小波是复指数乘高斯窗它的形态是“中间大、两边小”的连续振荡和故障脉冲的“一下子敲下去然后快速收敛”仍然不一样。结果就是你用小波阈值去噪阈值压低了脉冲细节也跟着被抹平阈值放宽了噪声又被保留。始终处于一个两难状态。本质上不是算法出了问题而是基函数和信号形态没有对齐。1.3 “基于模型的脉冲小波”想解决的三个核心矛盾既然通用基函数不匹配自然的思路是为这一类故障脉冲“量身定做”一种基函数也就是基于模型的脉冲小波。这个思路可以解决的矛盾有三个形态矛盾基函数本身就是一个“带阻尼的短促振荡”波形与真实故障脉冲高度相似。稀疏表示追求用尽量少的基函数逼近信号当基函数与信号足够“像”时几个原子就能精准刻画。参数矛盾真实脉冲不是固定不变的。转速变了、故障位置变了、系统固有频率变了脉冲的衰减快慢和振荡频率都会变。参数化原子通过调节时移、阻尼系数、固有频率三个参数让字典能够覆盖多种脉冲形态。可解释性矛盾小波分解系数是一大片分布哪个系数对应哪次冲击并不直观而稀疏表示的结果是一组明确选中的原子每个原子自带时移位置和参数直接告诉你“哪个时刻出现了什么形态的冲击”。也就是说这个方案本质上不是为了做漂亮的时频图而是为了把“一长段含噪信号”压缩成“少数几个带物理意义的脉冲原子参数”这对后续的故障频率识别非常有利。2. 参数化脉冲小波原子与字典设计从一行公式到千万级原子库2.1 原子表达式的推导逻辑三个参数控制一个波形在文献中脉冲小波原子有多种表达形式。不同文献的系数略有差别但核心结构一致都可以统一写成[ \psi_{\tau,\beta,f_n}(t) \begin{cases} A \cdot e^{-\beta (t-\tau)} \cdot \sin\big(2\pi f_n (t-\tau)\big), t \ge \tau \ 0, t \tau \end{cases} ]其中时移 τ决定冲击发生的时刻。对应滚动体经过损伤位置的时间点是后续计算故障特征频率的关键参数。阻尼系数 β决定振荡衰减的快慢。它和系统的阻尼比、接触刚度有关。β越大脉冲衰减越快波形越“尖”β越小波形拖得越长越接近连续振荡。固有频率 f_n决定振荡频率。实际中通常对应轴承座或传感器安装结构的谐振频率多在数百Hz到数千Hz范围。参数 A 是幅值。构造字典时每个原子会做L2归一化所以 A 本身不影响原子选择只影响后续系数重构。这看起来只是一个简单的衰减正弦但正是这个“简单”保证了它与故障物理过程的对齐。它不是随便选一个波形去拟合而是模型认为结构对冲击的响应就是这样衰减振荡的。2.2 离散化策略时移、频率、阻尼三个网格怎么设置把连续参数变成离散字典需要给三个维度划分网格。这里直接给出一套可用的工程化配置时移网格按采样间隔滑动。采样率 fs 12000 Hz信号时长 1 秒则时移 τ n/fs其中 n 1, 2, ..., 12000。时移网格的密度就等于采样点数。频率网格根据信号频谱和系统谐振情况圈定范围。假设固有频率在 1000~4000 Hz按步长 100 Hz 划档频率维度就有 31 个候选值。阻尼网格我没有把 β 划得过密。一般取 5~10 档即可例如 β 300, 600, 1000, 1500, 2000, 3000。原因后面会讲阻尼对波形的影响是渐变的网格太密只会让相邻原子高度相似对稀疏分解的稳定性没有帮助。如果直接做全排列字典原子总数 12000时移× 31频率× 6阻尼≈ 223万个原子。每个原子长度 12000 点以double存储的话光存字典就超过 200 GB。这是很多人初次尝试时马上碰到的硬墙。2.3 字典规模与内存问题不能硬存要按“子字典”结构组织我采用的方案是不直接存储完整大矩阵而是按频率和阻尼分组每组只存一个基础波形再用循环相关或分块卷积的方式去快速计算正交匹配追踪所需要的内积。具体来说设某一组频率 f_n、阻尼 β 对应的基础波形为 g(t)其定义域在 [0, L/fs] 上L 是波形有效长度。所有不同时移的原子其实就是这个基础波形的循环移位版本。字典与残差的内积 (\langle r, \psi_{\tau,\beta,f_n}\rangle) 本质上就是残差 r 与基础波形 g 的互相关函数。MATLAB里可以用xcorr或conv来做。计算时内存占用只跟一个基础向量和一个残差向量有关223万个原子变成一个十几万点的相关结果内存压力完全可控。这个方法的核心思想足够实用不是把字典建出来而是把字典“算出来”。下表是我在实验中使用的字典参数配置参考参数维度网格范围档数说明时移0 ~ T12000等于采样点数密度无需更高频率1000~4000 Hz16步长 200 Hz先粗搜再细搜阻尼300~30006对数均匀分布更合理信号长度0.5~1 秒-每次分解的输入信号段注意频率网格我只取了 16 档。这是因为大网格用于粗定位确定优势频率和阻尼后可以在该点附近做一次小范围的细化。两步搜索法能有效避免“字典太大算不动”和“字典太稀疏错过匹配”的两难。2.4 原子归一化一个容易被忽略但极其关键的步骤字典中的每个原子在做匹配之前必须进行L2归一化也就是把每个基本波形的能量归一。这一点的重要性在于如果不归一化幅值更大的原子在计算内积时天然占优OMP算法会优先选择它们而不是选择与残差波形最相似的那个原子。在实际实现时归一化的位置要小心。因为基础波形在循环移位过程中能量不会改变所以可以对每个基础波形先完成L2归一化再把移位后的原子用于相关计算。这样最后稀疏系数的幅值就代表了该原子在重构信号中的能量占比而不是原始幅值。3. OMP稀疏分解的MATLAB落地手写求解器与计算量控制3.1 为什么选择OMP而不是直接套凸优化稀疏表示单靠字典不够关键还得求解。常见路径有两条一条是凸优化类把稀疏问题写成L1范数最小化用迭代软阈值或内点法去解另一条是贪婪类以匹配追踪为代表每步挑最相关的原子然后逐步逼近信号。在MATLAB里做轴承故障诊断我推荐先用OMPOrthogonal Matching Pursuit正交匹配追踪。理由很简单凸优化方法需要额外工具箱且参数设置正则系数λ、迭代容差对结果影响很大OMP的每一步都有明确的几何含义迭代次数可以直接映射为“用多少个脉冲原子来解释这段信号”物理可解释性更强。3.2 OMP的核心逻辑与最小可运行代码OMP的算法本质是残差初始化为原始信号在字典中找与当前残差内积绝对值最大的那个原子把选中的原子加入支撑集合用最小二乘把信号投影到支撑集合张成的子空间上更新残差重复直到达到预设稀疏度或残差能量低于阈值。这里用pinv做最小二乘投影示例代码function [supp, coef] omp_solve(y, D, K) % y : Nx1 信号 % D : NxM 字典若内存不足请改用相关函数计算 % K : 最大原子数 r y; supp zeros(K,1); coef zeros(K,1); for ii 1:K corr D * r; [~, idx] max(abs(corr)); supp(ii) idx; A D(:, supp(1:ii)); coef(1:ii) pinv(A) * y; r y - A * coef(1:ii); if norm(r) / norm(y) 1e-6 supp supp(1:ii); coef coef(1:ii); break; end end end把上面的D换成“按需生成相关函数”的版本就是我在实验中真正使用的求解器框架。对于小规模的验证性实验这个代码已经够用。3.3 大规模字典下的OMP加速技巧一个直接的改进是把相关运算分解到子字典层。比如把字典按频率-阻尼分组每组计算一次相关向量 (\rho_{\beta,f_n} \text{xcorr}(r, g_{\beta,f_n}))然后取峰值。这样内外循环的维度从“千万级原子”降为“几十组参数×每次一组相关”速度能提高一到两个数量级。第二个改进是增量最小二乘。每轮迭代时支撑集合只是增加一个原子不需要重新对全部支撑原子求一遍pinv。用QR分解或Cholesky分解对 A^T A 进行更新复杂度从 (O(K^3)) 降到 (O(K^2))。当K达到几十甚至上百时这个差异非常明显。再一个比较实用的习惯是每次不要把整段几分钟的信号一次性丢给OMP。先按转轴周期切段每段做一次稀疏分解。这样既能减少单次计算规模又方便后续统计脉冲间隔。3.4 迭代次数K怎么定有时候残差下降曲线自己会说话K值太小脉冲细节不够K太大算法开始把噪声也当作有效原子。工程上我常用两个办法来确定合适的稀疏度肘部法则画一条“迭代次数-残差能量比”的曲线观察残差下降从陡峭变平缓的位置那个拐点就是有效稀疏度。自适应终止条件如果新加入原子后残差能量下降比例小于某个值例如1%就停止迭代。% 自适应终止示意 threshold 0.01; while iter K_max % 计算最优候选原子 % 获得加入后的残差能量 if (res_old - res_new) / res_old threshold break; end end一般来说一段1秒、采样率12000 Hz的轴承信号有效稀疏度在20~60之间就能得到很好的重构效果。超过100通常意味着信号段中包含了太多非脉冲成分或噪声干扰。4. 仿真验证与实测诊断全流程从“能看到稀疏系数”到“能报故障频率”4.1 构造已知答案的故障仿真信号先把整个链条放在一个可控环境里验证。我构造一个外圈故障仿真信号采样率 (f_s12000) Hz轴转频 (f_r29.5) Hz外圈故障特征频率约为 (f_{BPFO}107.8) Hz固有振荡频率取 2400 Hz阻尼系数 β1200。每次冲击产生一个衰减振荡原子再叠加高斯噪声信噪比约 6 dB。fs 12000; t (0:fs*1-1)/fs; fr 29.5; f_bpfo 107.8; fn 2400; beta 1200; y zeros(size(t)); for k 0:floor(f_bpfo * 1) tau k / f_bpfo randn*0.0002; % 加入微小随机滑移 idx find(t tau, 1); if isempty(idx), continue; end tt t(max(idx,1):end) - tau; n min(length(tt), round(fs * 0.02)); seg exp(-beta * tt(1:n)) .* sin(2*pi*fn*tt(1:n)); y(idx:idxn-1) y(idx:idxn-1) seg; end y y 0.5 * randn(size(t));这段代码不复杂但它忠实地反映了故障信号的物理结构确定周期、衰减振荡、随机滑移和噪声。有了这个已知答案后面每一步都能验证。4.2 OMP分解后如何看结果不要只盯着稀疏系数幅值对仿真信号做OMP分解选出的原子按参数排序。最容易犯的错误是直接看稀疏系数幅值列表并宣称“找到故障了”。实际上更有信息量的两个东西是被选中原子的时移序列。把这些τ从小到大排列计算相邻间隔。理想情况下间隔应该集中在一个故障特征周期附近约9.27ms。重构脉冲分量。用选中的原子和系数重构出纯脉冲信号再做Hilbert包络谱。在频谱上故障特征频率及其倍频会形成明显峰线。我两个方法都做了交叉验证后才下结论。4.3 从时移序列提取故障频率的两条路线路线一脉冲间隔统计法。对时移序列做一阶差分去掉极端异常值绘制直方图。峰的横坐标对应故障冲击周期。路线二重构信号包络谱法。用稀疏系数重构出脉冲成分做Hilbert变换取包络再FFT。峰值对应的频率即故障特征频率。下面这个表格列出了两种方法的适用场景方法优势不足适用场景时移间隔直方图直接利用稀疏表示的物理语义直观对滑移和漏检敏感信号信噪比高、脉冲间隔稳定重构包络谱抗噪性更强频域呈现稳定需要对字典顺序正确重构信噪比低、需要自动判读实测中我更喜欢两者配合先用间隔法快速确认周期再用包络谱去验证。4.4 把流程迁移到实测数据分段处理与转速波动应对经历过仿真验证后我拿某试验平台的数据跑了一遍。流程整理如下把原始振动信号按1秒一段做切分每段分别做OMP稀疏分解记录选中原子的时移和参数把多段结果的时移间隔做聚类用直方图寻找主导间隔对主导间隔换算成频率与轴承参数计算的理论故障特征频率对比。实测中有一个非常值得注意的现象由于转速存在轻微波动理论上的107.8 Hz在实际数据中并不会严格锁定在某一个值。它会在约103~112 Hz之间游动。如果字典或间隔统计只盯着单一频率做匹配很容易漏诊。更好的做法是在一个合理带宽比如理论特征频率的±6%~±8%内搜索峰值这比我一开始固定频率的做法稳健得多。5. 这四个坑我替你们先踩了字典归一化、迭代停止、共振重叠与算力妥协5.1 不归一化导致的选择偏差我最早实现OMP时为了让代码跑得快省略了原子归一化结果选出来的原子全都集中在幅值最大的那段信号附近而不是真正能解释故障脉冲的形态。加上L2归一化后相同参数下选原子顺序立刻合理。归一化注意守恒关系原子归一化后OMP给出的系数已是相对值重构信号时要把系数乘回原子的原始范数。我在项目里把原子原始范数单独存了一列向量防止重构幅度失真。5.2 迭代停止条件不能只看“次数”固定用K50跑所有信号很省事但不同信噪比下最优K差别很大。噪声强时K过大稀疏系数开始逐个解释噪声毛刺时移间隔变得杂乱噪声弱时K偏小少数弱冲击可能漏选间隔统计出现空洞。推荐做法是同时观察两条线残差能量曲线和时移间隔稳定性曲线。每次加入原子后看新原子加入是否让时移间隔直方图峰更集中如果更散说明已经进入噪声拟合区可以停了。这个方法有点“玄学”但工程中确实好用。5.3 阻尼小与固有频率接近时的多原子冗余问题当阻尼系数β很小振荡拖得很长或者两个不同频率的原子都跟残差有较高相关时OMP会倾向于用两三个相邻原子“合力”表示同一个物理脉冲。这会导致时移间隔出现虚假的短间隔比如明明故障周期是9.27ms结果相邻原子时移间隔只有1.2ms干扰了直方图判读。解决办法是在后处理中做原子合并把时移差小于某个阈值例如 0.5/ f_n的原子合并成一个簇簇的等效时移取能量加权平均。这个操作对提取故障周期非常重要。5.4 算力不够时的妥协方案分段、池化、两步细搜在普通笔记本上处理几十秒的实测信号直接用完整字典做OMP不可行。我的工程降级方案是第一步用原始信号段和粗网格字典频率步长200Hz、阻尼4档做OMP找出主导频率和阻尼范围第二步把字典细化到刚才锁定的小范围用细化后的子字典重新OMP第三步信号切段后并行跑多个分段。这样实际运行时间从不可忍变成分钟级结果稳定性还更好。如果是纯教学验证也可以进一步降低采样率或信号长度保证一个上午就能跑通全流程。至于最终的收敛判断我认为应该结合信号本身的物理意义当你看到OMP选出的时移序列其间隔分布有一个明显的主峰而这个主峰对应的频率与轴承特征频率误差小于3%时就可以很有把握报告故障了。不要把稀疏表示当作黑盒它的每一条输出都应该是可以回到物理模型去解释的。我在实际测试中发现把粗搜-细搜两步法加上原子合并处理后即使是信噪比低到4 dB左右的信号也能稳定识别出外圈故障特征频率。这个结果比我原先用连续小波变换包络谱的做法要干净得多。这就是“让基函数适应物理模型而不是让物理模型迁就现成基函数”的回报。