企业官网建设流程全解析

说实话，做了几年振动信号分析，我在轴承故障诊断这件事上有过一个比较强烈的感受：常规小波分析不是不能用，而是总觉得差那么一口气。明明频谱图上转频、倍频都在，可故障特征频率就是不够干净，提取得模棱两可。直到我把思路从“拿现成小波去套信号”转成“给信号专门设计一种脉冲小波原子，再用稀疏表示去匹配”，整个流程才豁然开朗。

这篇文章就以此为切入点，完整记录我在MATLAB里如何基于模型的脉冲小波（model-based impulsive wavelet）构造稀疏字典、编写OMP分解程序，并且把它应用到轴承故障诊断中的过程。内容包含原理拆解、代码实现、仿真验证、实测数据推断，以及我在这个过程中踩过的几个值得注意的坑。适合正在做故障诊断研究的学生，也适合刚接触稀疏表示、想在一个具体场景里落地这项技术的工程师参考。

1. 轴承故障脉冲为什么选通用小波总差一口气：信号形态与基函数错位

1.1 轴承故障信号的真实长相：不是正弦波，而是一连串“敲击声”

滚动轴承出现局部损伤（外圈点蚀、内圈剥落、滚动体缺陷）时，滚动体每经过损伤位置就会产生一次冲击。这种冲击的持续时间极短，通常只有几毫秒甚至更短，但它会立刻激起轴承座和传感器安装系统的固有振动。所以传感器拾取到的信号，数学上看是这样一个东西：

在某个时刻 t = τ 处突然出现一个幅值跳变，随后信号以系统固有频率 f_n 做振荡，同时振荡幅度按指数规律衰减：大约经过几个周期就被阻尼吃掉。然后隔一段时间，下一次冲击又来了。

写成表达式，就是典型的“阻尼正弦振荡”形态：其有效部分只在冲击发生后的很短窗口内存在，其余时间几乎为零。这个物理特征决定了它与连续平稳的正弦波、甚至与常见的连续小波基函数有本质区别。

1.2 通用小波库的错位在哪里：尺度匹配不等于形态匹配

很多人一开始会选择Daubechies小波、Symlets小波或者Morlet小波去做连续小波变换，直觉是“小波擅长捕捉突变信号”。这话没错，小波变换确实有良好的时频局部化能力。但问题出在基函数形态上：

Daubechies小波虽然有紧支撑，但它没有内在的“阻尼振荡”结构，用它去分解一段脉冲信号，需要非常多层的细节系数才能凑出那一下冲击；Morlet小波是复指数乘高斯窗，它的形态是“中间大、两边小”的连续振荡，和故障脉冲的“一下子敲下去然后快速收敛”仍然不一样。

结果就是：你用小波阈值去噪，阈值压低了，脉冲细节也跟着被抹平；阈值放宽了，噪声又被保留。始终处于一个两难状态。本质上不是算法出了问题，而是基函数和信号形态没有对齐。

1.3 “基于模型的脉冲小波”想解决的三个核心矛盾

既然通用基函数不匹配，自然的思路是为这一类故障脉冲“量身定做”一种基函数，也就是基于模型的脉冲小波。这个思路可以解决的矛盾有三个：

1. 形态矛盾：基函数本身就是一个“带阻尼的短促振荡”，波形与真实故障脉冲高度相似。稀疏表示追求用尽量少的基函数逼近信号，当基函数与信号足够“像”时，几个原子就能精准刻画。
2. 参数矛盾：真实脉冲不是固定不变的。转速变了、故障位置变了、系统固有频率变了，脉冲的衰减快慢和振荡频率都会变。参数化原子通过调节时移、阻尼系数、固有频率三个参数，让字典能够覆盖多种脉冲形态。
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万个原子变成一个十几万点的相关结果，内存压力完全可控。这个方法的核心思想足够实用：不是把字典建出来，而是把字典“算出来”。

下表是我在实验中使用的字典参数配置参考：

注意频率网格我只取了 16 档。这是因为大网格用于粗定位，确定优势频率和阻尼后，可以在该点附近做一次小范围的细化。两步搜索法能有效避免“字典太大算不动”和“字典太稀疏错过匹配”的两难。

2.4 原子归一化：一个容易被忽略但极其关键的步骤

字典中的每个原子在做匹配之前必须进行L2归一化，也就是把每个基本波形的能量归一。这一点的重要性在于：如果不归一化，幅值更大的原子在计算内积时天然占优，OMP算法会优先选择它们，而不是选择与残差波形最相似的那个原子。

在实际实现时，归一化的位置要小心。因为基础波形在循环移位过程中能量不会改变，所以可以对每个基础波形先完成L2归一化，再把移位后的原子用于相关计算。这样，最后稀疏系数的幅值就代表了该原子在重构信号中的能量占比，而不是原始幅值。

3. OMP稀疏分解的MATLAB落地：手写求解器与计算量控制

3.1 为什么选择OMP而不是直接套凸优化

稀疏表示单靠字典不够，关键还得求解。常见路径有两条：一条是凸优化类，把稀疏问题写成L1范数最小化，用迭代软阈值或内点法去解；另一条是贪婪类，以匹配追踪为代表，每步挑最相关的原子，然后逐步逼近信号。

在MATLAB里做轴承故障诊断，我推荐先用OMP（Orthogonal Matching Pursuit，正交匹配追踪）。理由很简单：凸优化方法需要额外工具箱，且参数设置（正则系数λ、迭代容差）对结果影响很大；OMP的每一步都有明确的几何含义，迭代次数可以直接映射为“用多少个脉冲原子来解释这段信号”，物理可解释性更强。

3.2 OMP的核心逻辑与最小可运行代码

OMP的算法本质是：

1. 残差初始化为原始信号；
2. 在字典中找与当前残差内积绝对值最大的那个原子；
3. 把选中的原子加入支撑集合；
4. 用最小二乘把信号投影到支撑集合张成的子空间上，更新残差；
5. 重复直到达到预设稀疏度或残差能量低于阈值。

这里用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. 肘部法则：画一条“迭代次数-残差能量比”的曲线，观察残差下降从陡峭变平缓的位置，那个拐点就是有效稀疏度。
2. 自适应终止条件：如果新加入原子后，残差能量下降比例小于某个值（例如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_s=12000) Hz，轴转频 (f_r=29.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:idx+n-1) = y(idx:idx+n-1) + seg; end y = y + 0.5 * randn(size(t));

这段代码不复杂，但它忠实地反映了故障信号的物理结构：确定周期、衰减振荡、随机滑移和噪声。有了这个已知答案，后面每一步都能验证。

4.2 OMP分解后如何看结果：不要只盯着稀疏系数幅值

对仿真信号做OMP分解，选出的原子按参数排序。最容易犯的错误是直接看稀疏系数幅值列表并宣称“找到故障了”。实际上更有信息量的两个东西是：

• 被选中原子的时移序列。把这些τ从小到大排列，计算相邻间隔。理想情况下，间隔应该集中在一个故障特征周期附近（约9.27ms）。
• 重构脉冲分量。用选中的原子和系数重构出纯脉冲信号，再做Hilbert包络谱。在频谱上，故障特征频率及其倍频会形成明显峰线。

我两个方法都做了，交叉验证后才下结论。

4.3 从时移序列提取故障频率的两条路线

路线一：脉冲间隔统计法。对时移序列做一阶差分，去掉极端异常值，绘制直方图。峰的横坐标对应故障冲击周期。

路线二：重构信号包络谱法。用稀疏系数重构出脉冲成分，做Hilbert变换取包络，再FFT。峰值对应的频率即故障特征频率。

下面这个表格列出了两种方法的适用场景：

实测中我更喜欢两者配合：先用间隔法快速确认周期，再用包络谱去验证。

4.4 把流程迁移到实测数据：分段处理与转速波动应对

经历过仿真验证后，我拿某试验平台的数据跑了一遍。流程整理如下：

1. 把原始振动信号按1秒一段做切分；
2. 每段分别做OMP稀疏分解，记录选中原子的时移和参数；
3. 把多段结果的时移间隔做聚类，用直方图寻找主导间隔；
4. 对主导间隔换算成频率，与轴承参数计算的理论故障特征频率对比。

实测中有一个非常值得注意的现象：由于转速存在轻微波动，理论上的107.8 Hz在实际数据中并不会严格锁定在某一个值。它会在约103~112 Hz之间游动。如果字典或间隔统计只盯着单一频率做匹配，很容易漏诊。更好的做法是在一个合理带宽（比如理论特征频率的±6%~±8%）内搜索峰值，这比我一开始固定频率的做法稳健得多。

5. 这四个坑我替你们先踩了：字典归一化、迭代停止、共振重叠与算力妥协

5.1 不归一化导致的选择偏差

我最早实现OMP时，为了让代码跑得快，省略了原子归一化，结果选出来的原子全都集中在幅值最大的那段信号附近，而不是真正能解释故障脉冲的形态。加上L2归一化后，相同参数下选原子顺序立刻合理。

归一化注意守恒关系：原子归一化后，OMP给出的系数已是相对值；重构信号时要把系数乘回原子的原始范数。我在项目里把原子原始范数单独存了一列向量，防止重构幅度失真。

5.2 迭代停止条件不能只看“次数”

固定用K=50跑所有信号很省事，但不同信噪比下最优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左右的信号，也能稳定识别出外圈故障特征频率。这个结果比我原先用连续小波变换+包络谱的做法要干净得多。这就是“让基函数适应物理模型，而不是让物理模型迁就现成基函数”的回报。