我先把话放在前面如果你是被标题里的“非下采样小波包”六个字劝退的人这篇文章就是写给你看的。搞设备故障诊断这几年我试过时域指标、频域包络、EMD分解、传统小波包最后真正在工程现场扛住考验的反而是听起来最“绕口”的非下采样小波包分析NSWPT。尤其是在MATLAB r2021b环境下用自带的最大重叠离散小波包变换modwpt函数几行代码就能把轴承故障特征从强噪声里捞出来效果比传统小波包稳得多。这个项目的核心任务其实就一句话用MATLAB r2021b实现非下采样小波包分解对轴承振动信号做时频分析提取故障特征频率从而判断轴承是外圈坏了、内圈坏了还是滚动体磨损。适合正在做故障诊断课题的研究生、设备状态监测的工程师以及所有被“轴承故障特征提取”折磨过的人。下面我把整个原理、选型、代码、踩坑过程全部拆开讲。1. 项目概述与技术背景1.1 为什么要盯住轴承故障诊断不放轴承是旋转机械里最容易坏的部件没有之一。电机、泵、齿轮箱、风机只要带转子的设备轴承一旦出问题轻则振动超标停机重则整条产线瘫痪。我见过太多现场案例操作工听到异响了还在硬扛结果轴承保持架碎裂转子直接扫膛维修成本翻了几倍。所以轴承故障诊断不是“锦上添花”的学术研究而是实打实的设备维护刚需。轴承故障诊断的核心挑战在于故障特征信号非常微弱而且被强烈的背景噪声淹没。正常运转时设备本身的振动、齿轮啮合冲击、流体扰动都会叠加在传感器信号上。轴承早期故障比如剥落坑直径零点几毫米引起的冲击能量极其有限在时域波形上几乎看不出来。这时就需要用信号处理方法把隐藏在宽带信号里的周期性冲击成分剥离出来。传统做法是包络分析Envelope Analysis先对信号做带通滤波再取包络谱找特征频率。这个方法在固定转速、单一故障的实验室环境下很好用但到了现场转速波动、多振源耦合、传递路径复杂包络谱经常被噪声带偏。这也是我后来转向非下采样小波包分析的根本原因——我们需要一种更精细的“频率放大镜”能把信号按频带无损地拆开再逐个分析每个频带的能量和冲击特征。1.2 为什么选择非下采样小波包而不是传统小波包很多人在故障诊断里用的是小波变换或小波包变换。传统小波包WPT确实比普通小波分析更进一步它能把高频段也细分不像离散小波变换只对低频逐层分解高频细节丢失得很厉害。但传统小波包有一个致命缺陷每层分解后都做了下采样抽取导致分解结果不具有平移不变性。意思是信号稍微平移几个采样点小波包系数的分布就会发生明显变化。对轴承故障信号来说平移不变性太重要了。因为故障冲击是瞬态的它的相位位置直接决定了包络特征。传统小波包在下采样时可能把冲击样本恰好抽掉或者在重构时出现频率混叠导致特征频率被扭曲。我早期用wpdec做轴承外圈故障分析同一个数据集截取不同起始点算出来的特征频带能量能差出30%以上这就是下采样的锅。非下采样小波包NSWPT的思路完全不一样分解过程中不对系数做下采样而是对滤波器进行上采样膨胀。每一层分解后子带信号长度和原信号完全一致总数据量逐层翻倍冗余分解换来的是严格的平移不变性和无混叠分解。MATLAB r2021b中的modwpt函数就是这种思想的工程实现。实测下来同样一组轴承数据非下采样小波包提取的特征频带稳定性比传统小波包高很多故障识别的可重复性完全不在一个量级。2. 核心原理深度拆解2.1 从小波包到非下采样小波包的演进逻辑要理解非下采样小波包得先搞清楚小波包到底在做什么。小波变换的本质是用一组带通滤波器把信号分解为不同频带的成分。一维离散小波变换DWT每一层只对低频近似分量继续分解高频细节分量不动导致高频分辨率差。小波包变换WPT改变了这个策略它同时分解低频和高频分支每一层都产生两个子带所以最终能在整个频域上形成均匀的二叉树结构。比如采样率是12kHz对信号做4层小波包分解会得到2的4次方即16个子带每个子带的理论带宽是12000/2^4 750Hz。这个频带划分非常精细轴承故障引起的共振频带通常耦合在几千赫兹的固有频率上可以被精准锁定位在某几个子带里。传统小波包的实现方式是Mallat算法每经一层分解就做一次二抽取每个子带采样率减半。这样做的好处是计算效率高、内存占用小但代价就是平移敏感和混叠。非下采样小波包把抽取环节去掉了取而代之的是对分解滤波器做隔点插零即膨胀/上采样。第j层分解时滤波器膨胀系数为2^(j-1)所以信号长度始终不变子带数量逐层翻倍。这个设计理念上的差异带来了三个直接影响第一平移不变性。信号在时域上无论怎么平移分解系数的能量分布保持一致。对故障诊断来说这意味着特征提取不依赖于截取窗口的起点结果更加稳定可靠。第二频率局部性更好。没有下采样就没有混叠效应各子带之间频率独立清晰不会再出现传统小波包中相邻频带互相污染的问题。第三计算代价上升。每一层的数据量都是上一层的两倍因为各分支都保持原长度4层分解的系数总量是原信号的16倍。好在MATLAB r2021b的modwpt基于快速滤波算法实现纯CPU处理1分钟长度的振动数据只需几秒工程上完全可接受。2.2 非下采样小波包的数学描述与算法实现细节从数学上看非下采样小波包分解可以表示为一个迭代滤波过程。设原始信号为 (x[n])分解层数为 (J)第 (j) 层有 (2^j) 个子带。定义两组基本滤波器尺度低通滤波器 (h[n]) 和小波高通滤波器 (g[n])它们来自所选小波基比如Daubechies系列的db4。在第 (j) 层滤波器被膨胀为 (h_j[n]) 和 (g_j[n])其中 (h_j[n]) 在相邻非零抽头之间插入 (2^{j-1}) 个零。第 (j1) 层的第 (2k) 个子带由第 (j) 层的第 (k) 个子带经低通滤波得到第 (2k1) 个子带则由高通滤波得到。由于没有下采样每个子带输出长度都是 (N)原始信号长度。用MATLAB的modwpt函数时内部还会对滤波器做归一化处理确保分解后各子带能量之和等于原信号能量Parseval定理成立。这一点对能量特征提取非常关键因为我们需要拿各子带能量占比作为故障判别指标能量守恒是前提。MATLAB r2021b里modwpt的完整调用格式是[WPT, wfreq] modwpt(x, lev, wname)其中x是输入信号lev是分解层数wname是小波基名称如db4、sym5、fk18返回的WPT是一个2^lev × N的矩阵wfreq是对应每个子带的中心频率单位是归一化频率需要乘以采样率换算成Hz。这里有一个容易被忽略的细节modwpt的输出子带排序方式和传统小波包wpdec不同。传统小波包按频率从低到高排列节点索引而modwpt的排列顺序遵循小波包树的自然顺序实际频率并不是严格单调排列。用wfreq输出可以直接拿到每个子带的中心频率避免排序混淆。2.3 为什么选择MATLAB r2021b这套组合选定非下采样小波包之后很容易陷入一个纠结用MATLAB哪个版本合适用Python自己写还是用MATLAB我个人的建议很明确——直接用MATLAB r2021b没必要在早期阶段自己造轮子。r2021b在信号处理领域的实用性很突出。它内置了Wavelet Toolbox 5.4版本modwpt函数经过多轮迭代算法稳定性和边界处理都做得很成熟。而且这个版本对内存管理做了优化处理百万数据点级别的信号时基本不卡顿。对比老版本r2021b的modwpt支持更多小波基选择尤其是Fejer-Korovkin族滤波器fk18、fk22它们在频带分离度上优于经典的db系列对轴承故障特征提取效果更好。再有一点r2021b的App Designer和脚本交互体验比之前的版本更好调试代码时能看到实时的变量变化和信号图形。实际做故障诊断往往需要反复试验不同的分解层数和小波基交互流畅度直接决定效率。我试过在Python里用PyWavelets实现非下采样小波包代码量更长边界处理还要自己操心并不比MATLAB省事。所以这篇文章的定位就很清楚了基于MATLAB r2021b工程版用现成工具箱完成任务把精力集中在信号分析和故障判断上而不是抠底层算法实现。3. 诊断流程与方案设计3.1 从原始振动信号到故障特征的整体路径轴承故障诊断不是“一个函数搞定一切”的魔法而是一套完整的数据处理链。我在实际项目中总结的流程是信号采集 → 预处理 → 非下采样小波包分解 → 特征提取 → 特征筛选 → 故障识别 → 结果验证。每一步都有必须注意的坑少一步都可能前功尽弃。首先是信号采集。加速度传感器贴在最靠近轴承座的位置采样率的选择很讲究。轴承故障特征频率比如外圈BPFO、内圈BPFI、滚动体BSF一般在几十赫兹到几百赫兹之间但它们激发的共振频带往往高达几千赫兹。为了捕捉到调制在共振频带上的冲击成分采样率至少要高于共振频率的2倍。我通常设采样率为12kHz或20kHz这样既能覆盖主要共振频带又不会产生海量数据。预处理环节里有一个高频操作必须做去除直流分量和高频噪声。去直流直接用x x - mean(x)即可高频噪声可以用带通滤波但要注意滤波器会导致信号相位失真。如果你打算做包络谱分析零相位滤波器filtfilt要比普通filter函数更合适。我在预处理时还会检查信号里有没有冲击性干扰比如对敲产生的毛刺这些非轴承故障引起的瞬态冲击会把后续的时频分析彻底带歪。3.2 非下采样小波包分解的层数与小波基选择策略分解层数和小波基选择是整个流程里最依赖经验的两个参数。层数太少频率分辨率不够故障特征频率和相邻频带混在一起分不出来层数太多计算量大增而且子带变窄后每个子带内的能量减少反而降低信噪比。经验准则是让子带宽度和故障特征频率的量级匹配。以采样率12kHz为例做4层分解得到16个子带每个子带宽750Hz适合捕捉轴承故障激发的宽频共振包络做5层分解得到32个子带每个子带宽375Hz如果共振频带更集中这个分辨率更有优势。我通常先用4层跑一遍看能量分布如果相邻几个子带能量都很高说明频带没拆干净再尝试提升到5层。小波基的选择上db4和sym5是通用首选fk18Fejer-Korovkin长度18是进阶选项。Daubechies小波系的特点是紧支撑性好db4的滤波器长度适中既不会因为太短导致频带泄漏也不会因为太长增加计算量适合轴承数据的通用分析。fk18的频带分离能力更强对共振频带边界更清晰但过度依赖精确的频带定位如果信号里的共振耦合比较弱效果反而不好。我可以分享一个筛选小技巧用同一组故障数据分别用不同小波基分解比较故障特征频带子带的“包络谱峰值因子”峰值频率处的幅值除以包络谱均方根值数值越大说明特征越突出。选峰值因子最高的小波基作为该场景的分析基准。3.3 特征提取与故障识别策略分解完成后每个子带都得到一组长度等于原始信号的系数序列。接下来要回答的问题是这些系数怎么变成“故障判断”的依据最直接的特征是子带能量。正常轴承运转时振动能量主要分布在低频段轴承出现故障后冲击激发了轴承系统的固有频率能量会聚集到某个中高频子带。所以各子带能量占比的变化本身就是一种有效的故障指示。再进一步子带包络谱的特征频率幅值也可以作为特征——在外圈故障特征频率处出现明显谱峰就能确诊外圈故障。我常用的特征向量包括三部分各子带归一化能量、各子带包络谱特征频率幅值、子带信息熵。子带信息熵衡量的是子带系数的稀疏性——冲击性越强熵值越小纯噪声信号熵值接近最大值。这三个维度拼出的特征向量可以很好地区分外圈故障、内圈故障和滚动体故障。故障识别这一步可以用简单的阈值判断也可以用机器学习分类器SVM、随机森林等。实验室研究中常用分类器来评估特征有效性但在工程现场我更推荐先用频带能量和包络谱特征做人工判读积累一段时间的数据后再训练分类器实现自动诊断。原因很现实现场工况复杂没有足够的标注数据之前分类器容易过拟合不如专家规则可靠。4. MATLAB r2021b代码实现全流程4.1 环境准备与数据组织在开始写代码之前先确认你的MATLAB r2021b安装了Wavelet Toolbox和Signal Processing Toolbox。检查办法很简单% 检查工具箱版本 ver(wavelet) ver(signal)如果没有安装在命令行执行matlab.addons.toolbox.installToolbox或者到附加功能里补装。这个步骤虽然基础但我在换电脑、装新版环境时经常忘记导致后续一堆函数报Undefined function错误。数据准备方面如果你手头有实验台采集的振动数据比如凯斯西储大学轴承数据集直接加载即可。如果是现场采集的原始数据需要做一点格式适配我习惯把振动信号存成CSV或MAT文件字段包含信号数组和采样率。% 数据加载示例 load(bearing_data.mat); % 包含 x: 振动信号序列, fs: 采样率(Hz) x x(:); % 转为列向量 fs fs; % 例如 12000如果暂时没有真实轴承故障数据可以用仿真信号调试流程。用一个简化的轴承外圈故障仿真信号做代码测试fs 12000; % 采样率 T 1; % 时长1秒 t (0:fs*T-1)/fs; fr 30; % 转频 30Hz BPFO 108.7; % 外圈故障特征频率示例值 n_impact floor(T * BPFO); % 生成脉冲序列间隔为 1/BPFO impulse zeros(size(t)); for idx 1:n_impact pos round(idx / BPFO * fs) 1; if pos length(t) impulse(pos) 1; end end % 模拟系统共振 wn_res 2500; % 共振频率 2500Hz xi 0.05; % 阻尼比 [den, num] ord2(wn_res, xi); res lsim(num*wn_res^2, den, impulse, t); % 加入噪声 x_sim res 0.3 * randn(size(t)); x x_sim;这段仿真代码的思路是在轴承故障特征频率处构造冲击序列经过一个二阶共振系统模拟轴承固有频率的响应最后叠加噪声。如果你有真实数据把这部分的代码替换掉即可。4.2 核心实现非下采样小波包分解与特征提取预处理和分解部分核心代码非常紧凑%% 预处理 x x - mean(x); % 去直流 x x / max(abs(x)); % 幅值归一化可选 %% 非下采样小波包分解 level 4; wname fk18; % 首选fk18对比时可改用db4 [WPT, wfreq] modwpt(x, level, wname); nBands size(WPT, 1); % 子带数 2^level % 查看子带中心频率换算为Hz freq_axis wfreq * fs; disp(各子带中心频率(Hz):); disp(freq_axis.);分解结果WPT是一个16行的矩阵每一行对应一个子带的时域波形。到这里我们已经把原始信号按频带精细地拆开了。接下来是特征提取模块我在这里同时计算子带能量和包络谱特征频率幅值%% 特征提取 % 子带能量及归一化比例 bandEnergy sum(WPT.^2, 2); energyRatio bandEnergy / sum(bandEnergy); % 对各子带分别做包络谱分析找特征频率处的幅值 BPFO 108.7; % 根据轴承参数计算 target_freq BPFO; % 若分析内圈或滚动体替换此值 env_amp_band zeros(nBands, 1); for k 1:nBands % 对子带信号取包络Hilbert变换求解析信号幅值 env abs(hilbert(WPT(k, :))); % 去掉包络的直流分量再做FFT env env - mean(env); N length(env); f_env (0:N-1) * fs / N; env_fft abs(fft(env)); % 在特征频率附近找最大峰值 band_idx find(f_env target_freq - 3 f_env target_freq 3); if ~isempty(band_idx) env_amp_band(k) max(env_fft(band_idx)); end end % 子带信息熵 bandEntropy zeros(nBands, 1); for k 1:nBands coef WPT(k, :); p coef.^2 / sum(coef.^2); bandEntropy(k) -sum(p .* log2(p eps)); end % 组合特征向量 featureVector [energyRatio; env_amp_band; bandEntropy];这段代码里包络谱分析用了Hilbert变换来提取信号的包络然后在包络谱的目标特征频率附近找一个窄带范围内的最大峰值。选择窄带搜索而不是取单一频率点的幅值是为了应对实际信号中转速微波动导致的频率偏移。4.3 故障状态识别与可视化得到特征向量之后故障识别的逻辑相对直接。以轴承外圈故障为例典型特征是与BPFO相关的子带通常是共振频带对应的子带索引能量占比突出且该子带包络谱在BPFO处有显著谱峰。内圈故障则往往伴随明显的边带成分转频调制BPFI。一个简单的规则判断如下% 定位能量占比最高的三个子带 [~, topIdx] sort(energyRatio, descend); top3 topIdx(1:3); % 判断是否在目标特征频率附近出现明显峰值 % 计算前三能量子带的包络谱峰值因子峰值/均方根 peakFactor zeros(3, 1); for i 1:3 k top3(i); env abs(hilbert(WPT(k, :))) - mean(abs(hilbert(WPT(k, :)))); env_fft abs(fft(env)); peakFactor(i) max(env_fft) / (mean(env_fft(2:end)) eps); end % 如果峰值因子大于设定阈值判定存在故障 if max(peakFactor) 5 disp(检测到周期性冲击疑似轴承故障); else disp(未检测到明显故障特征); end实际工程中只用这个简单规则还不够我会把特征向量保存下来配合历史数据做趋势分析。比如每天提取一次特征向量观察能量比值随时间的漂移一旦某个子带能量占比连续几天升高就可以提前预警这才是状态监测的正确姿势。可视化也是项目交付的必要环节。我习惯把原始波形、子带能量柱状图、关键子带包络谱三张图放在一张figure里展示figure(Position, [100 100 1200 800]) subplot(3,1,1); plot(t, x); xlabel(时间(s)); ylabel(幅值); title(原始振动信号); grid on; subplot(3,1,2); bar(1:nBands, energyRatio); xlabel(子带索引); ylabel(归一化能量); title(非下采样小波包子带能量分布); grid on; subplot(3,1,3); % 选择能量最高的子带绘制包络谱 kBest top3(1); env abs(hilbert(WPT(kBest, :))) - mean(abs(hilbert(WPT(kBest, :)))); N length(env); f_env (0:N-1) * fs / N; env_fft abs(fft(env)); plot(f_env(1:N/2), env_fft(1:N/2)); xlabel(频率(Hz)); ylabel(幅值); title(sprintf(子带%d包络谱, kBest)); xlim([0 1000]); grid on;这张图基本可以拿来做诊断报告了。子带能量分布图一眼就能看出能量集中到哪个频带包络谱上目标特征频率处的谱峰则直接指向故障类型。5. 常见问题与排查技巧5.1 问题速查表实际调试过程中我踩过的坑和帮别人解决的报错整理成一个速查表按症状分类症状可能原因解决办法modwpt报错Invalid wavelet小波基名称拼写错误或工具箱未安装检查waveletfamilies命令确认可用小波基WPT矩阵维度与预期不符输入信号不是列向量统一使用x x(:)转为列向量分解速度极慢信号过长且层数过多可先降采样到2kHz再分析或降低分解层数能量分布无明显规律信号包含强非平稳干扰先做带通滤波预处理或增加分解层数包络谱无特征峰值Hilbert变换受噪声影响对子带先做带通滤波再做包络分析多次运行结果不一致数据截取起点不同统一从同一采样点开始或加窗函数fk18分解结果和db4差异大频带分隔特性不同用已知故障数据验证选择最优小波基5.2 场景化避坑经验汇总第一个大坑是分解层数贪多。新手容易觉得层数越多分辨率越高一上来就设6层、7层。但层数高意味着子带数指数增长6层已有64个子带每个子带的实际物理带宽只有几十赫兹轴承故障特征频率附近的能量被分割到多个相邻子带里每个子带的能量比例都不突出反而“稀释”了故障特征。我的经验是先用4层跑通流程确认该方案能捕捉到目标频带再根据实际频带宽度做微调。第二个坑是忽视仿真信号和真实信号的差异。仿真信号是理想化的共振频率、特征频率、信噪比都是设定好的算法流程很容易跑出漂亮结果。但真实轴承信号里滚动体滑移会导致特征频率有微小波动设备启停阶段转速变化会引入频率漂移这些都会让仿真中设定的参数失效。所以验证流程时最好留一部分真实数据做盲测而不是只看仿真结果。第三个坑是边界效应。任何基于滤波的信号分解方法都会在信号起始端产生边界失真。非下采样小波包分解也不例外起始端一段数据的系数可能与中部有很大偏差。如果你截取的信号段短边界区域占比大特征提取结果就会失真。解决办法是截取数据时预留一段“安全边距”分析时丢弃边界两侧各几百个采样点只保留中间稳定区域。第四个坑是把能量特征当唯一判据。子带能量分布对故障确实敏感但它也容易受到载荷变化、转速波动、传感器安装位置的影响。同样的外圈故障轻载和重载工况下的能量分布可能差别很大。所以我从第一版方案开始就坚持“能量分布包络谱特征频率子带熵”三特征融合不把鸡蛋放在一个篮子里。6. 总结与个人经验项目做到最后回头看整个技术路线我最大的体会是方法本身不难难的是对一个方法吃透边界条件。非下采样小波包分析在轴承故障诊断中的核心价值是把信号按照严格的频带邻域拆开既保留了时间分辨率又避免了传统小波包的混叠问题。但它的上限也由这些边界条件决定——分解层数要匹配信号长度和频率范围小波基要匹配信号形态特征提取要匹配故障性质这些都需要在具体数据上反复试错没有一劳永逸的万能参数。再说一个容易被忽视的小技巧做非下采样小波包分析前永远先看一眼原始信号的频谱。如果频谱上根本看不到明显的共振峰或者信号被随机噪声完全覆盖那再多的分解层数也提炼不出故障特征。先用功率谱密度估计pwelch函数大致判断信号是否有可分析的价值能省下后面大量的试错时间。这个项目用到的modwpt方法还有一个后续扩展方向把子带能量分布作为时变特征对不同时间段滑动提取特征向量形成二维的“频带-时间-能量”热力图可以观察故障特征随时间的演化趋势。如果再配合分帧处理甚至可以做变转速工况下的阶次跟踪分析。对于想深入做下去的朋友这个方向值得花时间试试。最后送大家一句话项目做到最后拼的不是算法多新颖而是每个参数你都踩过坑、知道它为什么这么设、换一个场景要怎么调。这些经验书上看不来代码跑出来。
基于MATLAB r2021b非下采样小波包的轴承故障诊断实战
我先把话放在前面如果你是被标题里的“非下采样小波包”六个字劝退的人这篇文章就是写给你看的。搞设备故障诊断这几年我试过时域指标、频域包络、EMD分解、传统小波包最后真正在工程现场扛住考验的反而是听起来最“绕口”的非下采样小波包分析NSWPT。尤其是在MATLAB r2021b环境下用自带的最大重叠离散小波包变换modwpt函数几行代码就能把轴承故障特征从强噪声里捞出来效果比传统小波包稳得多。这个项目的核心任务其实就一句话用MATLAB r2021b实现非下采样小波包分解对轴承振动信号做时频分析提取故障特征频率从而判断轴承是外圈坏了、内圈坏了还是滚动体磨损。适合正在做故障诊断课题的研究生、设备状态监测的工程师以及所有被“轴承故障特征提取”折磨过的人。下面我把整个原理、选型、代码、踩坑过程全部拆开讲。1. 项目概述与技术背景1.1 为什么要盯住轴承故障诊断不放轴承是旋转机械里最容易坏的部件没有之一。电机、泵、齿轮箱、风机只要带转子的设备轴承一旦出问题轻则振动超标停机重则整条产线瘫痪。我见过太多现场案例操作工听到异响了还在硬扛结果轴承保持架碎裂转子直接扫膛维修成本翻了几倍。所以轴承故障诊断不是“锦上添花”的学术研究而是实打实的设备维护刚需。轴承故障诊断的核心挑战在于故障特征信号非常微弱而且被强烈的背景噪声淹没。正常运转时设备本身的振动、齿轮啮合冲击、流体扰动都会叠加在传感器信号上。轴承早期故障比如剥落坑直径零点几毫米引起的冲击能量极其有限在时域波形上几乎看不出来。这时就需要用信号处理方法把隐藏在宽带信号里的周期性冲击成分剥离出来。传统做法是包络分析Envelope Analysis先对信号做带通滤波再取包络谱找特征频率。这个方法在固定转速、单一故障的实验室环境下很好用但到了现场转速波动、多振源耦合、传递路径复杂包络谱经常被噪声带偏。这也是我后来转向非下采样小波包分析的根本原因——我们需要一种更精细的“频率放大镜”能把信号按频带无损地拆开再逐个分析每个频带的能量和冲击特征。1.2 为什么选择非下采样小波包而不是传统小波包很多人在故障诊断里用的是小波变换或小波包变换。传统小波包WPT确实比普通小波分析更进一步它能把高频段也细分不像离散小波变换只对低频逐层分解高频细节丢失得很厉害。但传统小波包有一个致命缺陷每层分解后都做了下采样抽取导致分解结果不具有平移不变性。意思是信号稍微平移几个采样点小波包系数的分布就会发生明显变化。对轴承故障信号来说平移不变性太重要了。因为故障冲击是瞬态的它的相位位置直接决定了包络特征。传统小波包在下采样时可能把冲击样本恰好抽掉或者在重构时出现频率混叠导致特征频率被扭曲。我早期用wpdec做轴承外圈故障分析同一个数据集截取不同起始点算出来的特征频带能量能差出30%以上这就是下采样的锅。非下采样小波包NSWPT的思路完全不一样分解过程中不对系数做下采样而是对滤波器进行上采样膨胀。每一层分解后子带信号长度和原信号完全一致总数据量逐层翻倍冗余分解换来的是严格的平移不变性和无混叠分解。MATLAB r2021b中的modwpt函数就是这种思想的工程实现。实测下来同样一组轴承数据非下采样小波包提取的特征频带稳定性比传统小波包高很多故障识别的可重复性完全不在一个量级。2. 核心原理深度拆解2.1 从小波包到非下采样小波包的演进逻辑要理解非下采样小波包得先搞清楚小波包到底在做什么。小波变换的本质是用一组带通滤波器把信号分解为不同频带的成分。一维离散小波变换DWT每一层只对低频近似分量继续分解高频细节分量不动导致高频分辨率差。小波包变换WPT改变了这个策略它同时分解低频和高频分支每一层都产生两个子带所以最终能在整个频域上形成均匀的二叉树结构。比如采样率是12kHz对信号做4层小波包分解会得到2的4次方即16个子带每个子带的理论带宽是12000/2^4 750Hz。这个频带划分非常精细轴承故障引起的共振频带通常耦合在几千赫兹的固有频率上可以被精准锁定位在某几个子带里。传统小波包的实现方式是Mallat算法每经一层分解就做一次二抽取每个子带采样率减半。这样做的好处是计算效率高、内存占用小但代价就是平移敏感和混叠。非下采样小波包把抽取环节去掉了取而代之的是对分解滤波器做隔点插零即膨胀/上采样。第j层分解时滤波器膨胀系数为2^(j-1)所以信号长度始终不变子带数量逐层翻倍。这个设计理念上的差异带来了三个直接影响第一平移不变性。信号在时域上无论怎么平移分解系数的能量分布保持一致。对故障诊断来说这意味着特征提取不依赖于截取窗口的起点结果更加稳定可靠。第二频率局部性更好。没有下采样就没有混叠效应各子带之间频率独立清晰不会再出现传统小波包中相邻频带互相污染的问题。第三计算代价上升。每一层的数据量都是上一层的两倍因为各分支都保持原长度4层分解的系数总量是原信号的16倍。好在MATLAB r2021b的modwpt基于快速滤波算法实现纯CPU处理1分钟长度的振动数据只需几秒工程上完全可接受。2.2 非下采样小波包的数学描述与算法实现细节从数学上看非下采样小波包分解可以表示为一个迭代滤波过程。设原始信号为 (x[n])分解层数为 (J)第 (j) 层有 (2^j) 个子带。定义两组基本滤波器尺度低通滤波器 (h[n]) 和小波高通滤波器 (g[n])它们来自所选小波基比如Daubechies系列的db4。在第 (j) 层滤波器被膨胀为 (h_j[n]) 和 (g_j[n])其中 (h_j[n]) 在相邻非零抽头之间插入 (2^{j-1}) 个零。第 (j1) 层的第 (2k) 个子带由第 (j) 层的第 (k) 个子带经低通滤波得到第 (2k1) 个子带则由高通滤波得到。由于没有下采样每个子带输出长度都是 (N)原始信号长度。用MATLAB的modwpt函数时内部还会对滤波器做归一化处理确保分解后各子带能量之和等于原信号能量Parseval定理成立。这一点对能量特征提取非常关键因为我们需要拿各子带能量占比作为故障判别指标能量守恒是前提。MATLAB r2021b里modwpt的完整调用格式是[WPT, wfreq] modwpt(x, lev, wname)其中x是输入信号lev是分解层数wname是小波基名称如db4、sym5、fk18返回的WPT是一个2^lev × N的矩阵wfreq是对应每个子带的中心频率单位是归一化频率需要乘以采样率换算成Hz。这里有一个容易被忽略的细节modwpt的输出子带排序方式和传统小波包wpdec不同。传统小波包按频率从低到高排列节点索引而modwpt的排列顺序遵循小波包树的自然顺序实际频率并不是严格单调排列。用wfreq输出可以直接拿到每个子带的中心频率避免排序混淆。2.3 为什么选择MATLAB r2021b这套组合选定非下采样小波包之后很容易陷入一个纠结用MATLAB哪个版本合适用Python自己写还是用MATLAB我个人的建议很明确——直接用MATLAB r2021b没必要在早期阶段自己造轮子。r2021b在信号处理领域的实用性很突出。它内置了Wavelet Toolbox 5.4版本modwpt函数经过多轮迭代算法稳定性和边界处理都做得很成熟。而且这个版本对内存管理做了优化处理百万数据点级别的信号时基本不卡顿。对比老版本r2021b的modwpt支持更多小波基选择尤其是Fejer-Korovkin族滤波器fk18、fk22它们在频带分离度上优于经典的db系列对轴承故障特征提取效果更好。再有一点r2021b的App Designer和脚本交互体验比之前的版本更好调试代码时能看到实时的变量变化和信号图形。实际做故障诊断往往需要反复试验不同的分解层数和小波基交互流畅度直接决定效率。我试过在Python里用PyWavelets实现非下采样小波包代码量更长边界处理还要自己操心并不比MATLAB省事。所以这篇文章的定位就很清楚了基于MATLAB r2021b工程版用现成工具箱完成任务把精力集中在信号分析和故障判断上而不是抠底层算法实现。3. 诊断流程与方案设计3.1 从原始振动信号到故障特征的整体路径轴承故障诊断不是“一个函数搞定一切”的魔法而是一套完整的数据处理链。我在实际项目中总结的流程是信号采集 → 预处理 → 非下采样小波包分解 → 特征提取 → 特征筛选 → 故障识别 → 结果验证。每一步都有必须注意的坑少一步都可能前功尽弃。首先是信号采集。加速度传感器贴在最靠近轴承座的位置采样率的选择很讲究。轴承故障特征频率比如外圈BPFO、内圈BPFI、滚动体BSF一般在几十赫兹到几百赫兹之间但它们激发的共振频带往往高达几千赫兹。为了捕捉到调制在共振频带上的冲击成分采样率至少要高于共振频率的2倍。我通常设采样率为12kHz或20kHz这样既能覆盖主要共振频带又不会产生海量数据。预处理环节里有一个高频操作必须做去除直流分量和高频噪声。去直流直接用x x - mean(x)即可高频噪声可以用带通滤波但要注意滤波器会导致信号相位失真。如果你打算做包络谱分析零相位滤波器filtfilt要比普通filter函数更合适。我在预处理时还会检查信号里有没有冲击性干扰比如对敲产生的毛刺这些非轴承故障引起的瞬态冲击会把后续的时频分析彻底带歪。3.2 非下采样小波包分解的层数与小波基选择策略分解层数和小波基选择是整个流程里最依赖经验的两个参数。层数太少频率分辨率不够故障特征频率和相邻频带混在一起分不出来层数太多计算量大增而且子带变窄后每个子带内的能量减少反而降低信噪比。经验准则是让子带宽度和故障特征频率的量级匹配。以采样率12kHz为例做4层分解得到16个子带每个子带宽750Hz适合捕捉轴承故障激发的宽频共振包络做5层分解得到32个子带每个子带宽375Hz如果共振频带更集中这个分辨率更有优势。我通常先用4层跑一遍看能量分布如果相邻几个子带能量都很高说明频带没拆干净再尝试提升到5层。小波基的选择上db4和sym5是通用首选fk18Fejer-Korovkin长度18是进阶选项。Daubechies小波系的特点是紧支撑性好db4的滤波器长度适中既不会因为太短导致频带泄漏也不会因为太长增加计算量适合轴承数据的通用分析。fk18的频带分离能力更强对共振频带边界更清晰但过度依赖精确的频带定位如果信号里的共振耦合比较弱效果反而不好。我可以分享一个筛选小技巧用同一组故障数据分别用不同小波基分解比较故障特征频带子带的“包络谱峰值因子”峰值频率处的幅值除以包络谱均方根值数值越大说明特征越突出。选峰值因子最高的小波基作为该场景的分析基准。3.3 特征提取与故障识别策略分解完成后每个子带都得到一组长度等于原始信号的系数序列。接下来要回答的问题是这些系数怎么变成“故障判断”的依据最直接的特征是子带能量。正常轴承运转时振动能量主要分布在低频段轴承出现故障后冲击激发了轴承系统的固有频率能量会聚集到某个中高频子带。所以各子带能量占比的变化本身就是一种有效的故障指示。再进一步子带包络谱的特征频率幅值也可以作为特征——在外圈故障特征频率处出现明显谱峰就能确诊外圈故障。我常用的特征向量包括三部分各子带归一化能量、各子带包络谱特征频率幅值、子带信息熵。子带信息熵衡量的是子带系数的稀疏性——冲击性越强熵值越小纯噪声信号熵值接近最大值。这三个维度拼出的特征向量可以很好地区分外圈故障、内圈故障和滚动体故障。故障识别这一步可以用简单的阈值判断也可以用机器学习分类器SVM、随机森林等。实验室研究中常用分类器来评估特征有效性但在工程现场我更推荐先用频带能量和包络谱特征做人工判读积累一段时间的数据后再训练分类器实现自动诊断。原因很现实现场工况复杂没有足够的标注数据之前分类器容易过拟合不如专家规则可靠。4. MATLAB r2021b代码实现全流程4.1 环境准备与数据组织在开始写代码之前先确认你的MATLAB r2021b安装了Wavelet Toolbox和Signal Processing Toolbox。检查办法很简单% 检查工具箱版本 ver(wavelet) ver(signal)如果没有安装在命令行执行matlab.addons.toolbox.installToolbox或者到附加功能里补装。这个步骤虽然基础但我在换电脑、装新版环境时经常忘记导致后续一堆函数报Undefined function错误。数据准备方面如果你手头有实验台采集的振动数据比如凯斯西储大学轴承数据集直接加载即可。如果是现场采集的原始数据需要做一点格式适配我习惯把振动信号存成CSV或MAT文件字段包含信号数组和采样率。% 数据加载示例 load(bearing_data.mat); % 包含 x: 振动信号序列, fs: 采样率(Hz) x x(:); % 转为列向量 fs fs; % 例如 12000如果暂时没有真实轴承故障数据可以用仿真信号调试流程。用一个简化的轴承外圈故障仿真信号做代码测试fs 12000; % 采样率 T 1; % 时长1秒 t (0:fs*T-1)/fs; fr 30; % 转频 30Hz BPFO 108.7; % 外圈故障特征频率示例值 n_impact floor(T * BPFO); % 生成脉冲序列间隔为 1/BPFO impulse zeros(size(t)); for idx 1:n_impact pos round(idx / BPFO * fs) 1; if pos length(t) impulse(pos) 1; end end % 模拟系统共振 wn_res 2500; % 共振频率 2500Hz xi 0.05; % 阻尼比 [den, num] ord2(wn_res, xi); res lsim(num*wn_res^2, den, impulse, t); % 加入噪声 x_sim res 0.3 * randn(size(t)); x x_sim;这段仿真代码的思路是在轴承故障特征频率处构造冲击序列经过一个二阶共振系统模拟轴承固有频率的响应最后叠加噪声。如果你有真实数据把这部分的代码替换掉即可。4.2 核心实现非下采样小波包分解与特征提取预处理和分解部分核心代码非常紧凑%% 预处理 x x - mean(x); % 去直流 x x / max(abs(x)); % 幅值归一化可选 %% 非下采样小波包分解 level 4; wname fk18; % 首选fk18对比时可改用db4 [WPT, wfreq] modwpt(x, level, wname); nBands size(WPT, 1); % 子带数 2^level % 查看子带中心频率换算为Hz freq_axis wfreq * fs; disp(各子带中心频率(Hz):); disp(freq_axis.);分解结果WPT是一个16行的矩阵每一行对应一个子带的时域波形。到这里我们已经把原始信号按频带精细地拆开了。接下来是特征提取模块我在这里同时计算子带能量和包络谱特征频率幅值%% 特征提取 % 子带能量及归一化比例 bandEnergy sum(WPT.^2, 2); energyRatio bandEnergy / sum(bandEnergy); % 对各子带分别做包络谱分析找特征频率处的幅值 BPFO 108.7; % 根据轴承参数计算 target_freq BPFO; % 若分析内圈或滚动体替换此值 env_amp_band zeros(nBands, 1); for k 1:nBands % 对子带信号取包络Hilbert变换求解析信号幅值 env abs(hilbert(WPT(k, :))); % 去掉包络的直流分量再做FFT env env - mean(env); N length(env); f_env (0:N-1) * fs / N; env_fft abs(fft(env)); % 在特征频率附近找最大峰值 band_idx find(f_env target_freq - 3 f_env target_freq 3); if ~isempty(band_idx) env_amp_band(k) max(env_fft(band_idx)); end end % 子带信息熵 bandEntropy zeros(nBands, 1); for k 1:nBands coef WPT(k, :); p coef.^2 / sum(coef.^2); bandEntropy(k) -sum(p .* log2(p eps)); end % 组合特征向量 featureVector [energyRatio; env_amp_band; bandEntropy];这段代码里包络谱分析用了Hilbert变换来提取信号的包络然后在包络谱的目标特征频率附近找一个窄带范围内的最大峰值。选择窄带搜索而不是取单一频率点的幅值是为了应对实际信号中转速微波动导致的频率偏移。4.3 故障状态识别与可视化得到特征向量之后故障识别的逻辑相对直接。以轴承外圈故障为例典型特征是与BPFO相关的子带通常是共振频带对应的子带索引能量占比突出且该子带包络谱在BPFO处有显著谱峰。内圈故障则往往伴随明显的边带成分转频调制BPFI。一个简单的规则判断如下% 定位能量占比最高的三个子带 [~, topIdx] sort(energyRatio, descend); top3 topIdx(1:3); % 判断是否在目标特征频率附近出现明显峰值 % 计算前三能量子带的包络谱峰值因子峰值/均方根 peakFactor zeros(3, 1); for i 1:3 k top3(i); env abs(hilbert(WPT(k, :))) - mean(abs(hilbert(WPT(k, :)))); env_fft abs(fft(env)); peakFactor(i) max(env_fft) / (mean(env_fft(2:end)) eps); end % 如果峰值因子大于设定阈值判定存在故障 if max(peakFactor) 5 disp(检测到周期性冲击疑似轴承故障); else disp(未检测到明显故障特征); end实际工程中只用这个简单规则还不够我会把特征向量保存下来配合历史数据做趋势分析。比如每天提取一次特征向量观察能量比值随时间的漂移一旦某个子带能量占比连续几天升高就可以提前预警这才是状态监测的正确姿势。可视化也是项目交付的必要环节。我习惯把原始波形、子带能量柱状图、关键子带包络谱三张图放在一张figure里展示figure(Position, [100 100 1200 800]) subplot(3,1,1); plot(t, x); xlabel(时间(s)); ylabel(幅值); title(原始振动信号); grid on; subplot(3,1,2); bar(1:nBands, energyRatio); xlabel(子带索引); ylabel(归一化能量); title(非下采样小波包子带能量分布); grid on; subplot(3,1,3); % 选择能量最高的子带绘制包络谱 kBest top3(1); env abs(hilbert(WPT(kBest, :))) - mean(abs(hilbert(WPT(kBest, :)))); N length(env); f_env (0:N-1) * fs / N; env_fft abs(fft(env)); plot(f_env(1:N/2), env_fft(1:N/2)); xlabel(频率(Hz)); ylabel(幅值); title(sprintf(子带%d包络谱, kBest)); xlim([0 1000]); grid on;这张图基本可以拿来做诊断报告了。子带能量分布图一眼就能看出能量集中到哪个频带包络谱上目标特征频率处的谱峰则直接指向故障类型。5. 常见问题与排查技巧5.1 问题速查表实际调试过程中我踩过的坑和帮别人解决的报错整理成一个速查表按症状分类症状可能原因解决办法modwpt报错Invalid wavelet小波基名称拼写错误或工具箱未安装检查waveletfamilies命令确认可用小波基WPT矩阵维度与预期不符输入信号不是列向量统一使用x x(:)转为列向量分解速度极慢信号过长且层数过多可先降采样到2kHz再分析或降低分解层数能量分布无明显规律信号包含强非平稳干扰先做带通滤波预处理或增加分解层数包络谱无特征峰值Hilbert变换受噪声影响对子带先做带通滤波再做包络分析多次运行结果不一致数据截取起点不同统一从同一采样点开始或加窗函数fk18分解结果和db4差异大频带分隔特性不同用已知故障数据验证选择最优小波基5.2 场景化避坑经验汇总第一个大坑是分解层数贪多。新手容易觉得层数越多分辨率越高一上来就设6层、7层。但层数高意味着子带数指数增长6层已有64个子带每个子带的实际物理带宽只有几十赫兹轴承故障特征频率附近的能量被分割到多个相邻子带里每个子带的能量比例都不突出反而“稀释”了故障特征。我的经验是先用4层跑通流程确认该方案能捕捉到目标频带再根据实际频带宽度做微调。第二个坑是忽视仿真信号和真实信号的差异。仿真信号是理想化的共振频率、特征频率、信噪比都是设定好的算法流程很容易跑出漂亮结果。但真实轴承信号里滚动体滑移会导致特征频率有微小波动设备启停阶段转速变化会引入频率漂移这些都会让仿真中设定的参数失效。所以验证流程时最好留一部分真实数据做盲测而不是只看仿真结果。第三个坑是边界效应。任何基于滤波的信号分解方法都会在信号起始端产生边界失真。非下采样小波包分解也不例外起始端一段数据的系数可能与中部有很大偏差。如果你截取的信号段短边界区域占比大特征提取结果就会失真。解决办法是截取数据时预留一段“安全边距”分析时丢弃边界两侧各几百个采样点只保留中间稳定区域。第四个坑是把能量特征当唯一判据。子带能量分布对故障确实敏感但它也容易受到载荷变化、转速波动、传感器安装位置的影响。同样的外圈故障轻载和重载工况下的能量分布可能差别很大。所以我从第一版方案开始就坚持“能量分布包络谱特征频率子带熵”三特征融合不把鸡蛋放在一个篮子里。6. 总结与个人经验项目做到最后回头看整个技术路线我最大的体会是方法本身不难难的是对一个方法吃透边界条件。非下采样小波包分析在轴承故障诊断中的核心价值是把信号按照严格的频带邻域拆开既保留了时间分辨率又避免了传统小波包的混叠问题。但它的上限也由这些边界条件决定——分解层数要匹配信号长度和频率范围小波基要匹配信号形态特征提取要匹配故障性质这些都需要在具体数据上反复试错没有一劳永逸的万能参数。再说一个容易被忽视的小技巧做非下采样小波包分析前永远先看一眼原始信号的频谱。如果频谱上根本看不到明显的共振峰或者信号被随机噪声完全覆盖那再多的分解层数也提炼不出故障特征。先用功率谱密度估计pwelch函数大致判断信号是否有可分析的价值能省下后面大量的试错时间。这个项目用到的modwpt方法还有一个后续扩展方向把子带能量分布作为时变特征对不同时间段滑动提取特征向量形成二维的“频带-时间-能量”热力图可以观察故障特征随时间的演化趋势。如果再配合分帧处理甚至可以做变转速工况下的阶次跟踪分析。对于想深入做下去的朋友这个方向值得花时间试试。最后送大家一句话项目做到最后拼的不是算法多新颖而是每个参数你都踩过坑、知道它为什么这么设、换一个场景要怎么调。这些经验书上看不来代码跑出来。
