做旋转机械状态监测这些年,我越来越觉得,轴承故障诊断最难的不是“坏了怎么修”,而是“怎么在信号里把‘要坏’的证据找出来”。实际采集到的振动信号,绝大多数情况下就是一团乱麻:转频及其谐波能量巨大,齿轮啮合频率和边带搅在一起,环境噪声无孔不入,而早期轴承故障产生的冲击脉冲,能量微弱到经常被淹没在背景里,直接在原始频谱上看,根本找不到特征频率的影子。这就像在闹市区找一个低声说话的人,不把噪音隔开,你什么都听不见。
这个项目就是用**可调Q因子小波变换(TQWT)**来解决这个“闹市区找人”的问题,代码运行环境是MATLAB r2021b。它的核心思路不是把信号整体做一次滤波,而是把信号按振荡特性拆解成不同成分,然后把包含故障冲击的那部分挑出来做包络谱分析。整个过程跑通之后,就算信号里噪声很强、转频成分很突出,也能把轴承内圈、外圈、滚动体的故障特征频率清晰地暴露出来。这篇文章我会从原理、参数选型、代码实现到工程踩坑完整过一遍,适合正在做故障诊断课题的研究生,以及现场做设备状态监测、想把手头信号处理手段升级一下的工程师参考。
1. 项目概述与方案选型
1.1 为什么是TQWT而不是传统小波
先说清楚轴承振动信号里到底有什么。正常运行的轴承,振动以转频及其整数倍谐波为主,幅值相对稳定。一旦滚道或滚动体上出现局部缺陷,比如剥落、点蚀,每次滚过缺陷位置就会产生一个短暂的冲击脉冲。这个脉冲的特点是持续时间极短、频带很宽,并且会激励起轴承座或传感器支撑结构的固有共振,形成一段高频衰减振荡。
问题来了:转频谐波是持续的、低频的、大能量的;故障冲击是瞬态的、宽频的、小能量的;噪声则是无处不在的。传统带通滤波器需要先确定共振频带,但故障频率和共振频率往往是未知的,选错了频带等于白做。EMD这类经验分解方法看起来自适应,但模态混叠问题严重,尤其在强噪声下,第一阶IMF往往还是噪声主导。而普通小波变换的基函数固定,一旦选定小波基和分解层数,滤波特性就定死了,很难同时适配“持续谐波”和“瞬态冲击”这两种特性差别极大的成分。
TQWT不一样的地方在于,它通过调整品质因子Q来控制小波基的振荡能力。Q值低,小波更像一个尖锐的脉冲,适合匹配故障冲击;Q值高,小波更像一段持续振荡的波形,适合匹配谐波和共振衰减。这个特性让它可以直接面向“信号分离”来设计参数,而不是先猜共振频带。实际项目中我测试过,同样一段含噪内圈故障信号,用db4小波做五层分解后包络谱,特征频率处几乎看不到峰值;换成TQWT把冲击成分提取出来再做包络,特征频率和倍频清清楚楚。这个差距非常直观。
1.2 项目整体技术路线
整个诊断流程可以拆成五个环节:信号预处理、TQWT分解、形态分量分析重构、Hilbert包络谱分析、特征频率判读。预处理负责去均值、去趋势、剔除明显异常点;TQWT负责把信号按Q因子分解成多个子带;形态分量分析(MCA)负责从子带中把符合“瞬态冲击”特征的成分选出来重新组合;包络谱分析负责把重构信号的冲击周期转化为可读的频率峰值;最后一步是把峰值与理论计算的轴承故障特征频率比对,给出诊断结论。
这个路线在MATLAB r2021b里完全可行,不需要额外安装付费工具箱,基础信号处理函数加上自己写的TQWT函数就够了。r2021b的Hilbert变换、fft、filtfilt等内置函数性能足够稳定,处理8192点采样信号,TQWT分解加重构整个过程耗时在几百毫秒到一两秒之间,完全满足离线分析需求。需要注意的一点是,r2021b的wavelet toolbox内置了小波变换,但没有TQWT函数,需要自己实现或下载开源代码。文章后半部分我会给出一个完整可用的实现思路。
2. TQWT核心原理与参数详解
2.1 滤波器组结构与数学本质
TQWT由Selesnick于2011年提出,本质是一个过完备的离散小波变换,建立在双通道滤波器组迭代分解的基础上。每一层分解把信号分成低频子带和高频子带,低频子带继续送入下一层迭代,高频子带作为这一层的输出。这里的“低通”和“高通”并不是固定的,而是由两个缩放因子α和β控制,低通缩放因子α决定了低频子带保留多少频带,高通缩放因子β决定了高频子带覆盖频带的上限。
品质因子Q和冗余度r由这两个缩放因子决定:
Q = (2 - β) / β
r = β / (1 - α)
数学上看,Q因子决定了小波基在频域上的带宽与中心频率比值。Q越高,频带越窄,时域波形越接近持续振荡;Q越低,频带越宽,时域波形越尖锐、越接近冲击。冗余度r控制滤波器组在相邻层之间的频率覆盖重叠程度,r取3时是常用默认值,频率响应在频域上的覆盖比较均匀,既不会留下明显空隙,也不会过度冗余导致计算量暴增。
分解层数J决定了最高频子带对应的中心频率范围。每一层高频子带的中心频率大致按β的幂次逐层降低,如果J设置得过小,高频细节分解不充分;设置得过大,最低层的低频子带会被过度压缩,计算量也会显著增加。实际使用中J的取值要结合信号长度和采样率确定,一般让最低层低频子带的上限频率低于转频最低分析范围即可,我通常设置在15到40之间。
2.2 可调Q因子的工程含义
用一个生活化的类比来理解Q因子。想象一下敲击一个玻璃杯边缘发出的“叮”声——声音尖锐、很快消失,这就是低Q振荡,持续时间短,频带极宽;再想象持续按着琴弦发出的长音——音调稳定、余音悠长,这就是高Q振荡,频带窄、持续性强。轴承故障冲击恰好类似于敲击声,转频谐波和齿轮啮合分量则类似于琴弦长音。
TQWT的可贵之处在于,它允许你针对这两种成分分别设计不同的Q因子,然后并行分解同一段信号。实际操作中常用双Q因子策略:一个低Q通道负责提取瞬态冲击成分,比如Q等于4到6;一个高Q通道负责建模持续振荡成分,比如Q等于20到30。两个通道分别分解、分别重构,最后把高Q通道重构结果从原信号中减掉,剩下的就是包含故障冲击的低Q成分。
这个“两通道并行加相减”的思路在日常信号处理中非常实用。我从项目实践中得到的经验是,单靠一个固定Q值跑一遍TQWT,效果往往不如双通道并行;因为现场信号里干扰成分的种类太多,仅靠单一Q值的滤波器组很难把所有干扰都压制住,而双通道天然形成了一个“信号分裂器”,把谐波类成分和冲击类成分强行分开,效果立竿见影。
2.3 边界效应与子带选择
TQWT采用分块卷积和延拓处理,在信号两端仍然会引入边界效应,也就是重构信号在首尾一小段会出现幅值畸变。这个问题在轴承故障诊断中很容易被忽略,因为大家通常只盯着包络谱峰值看,但边界畸变会在特定频带引入假冲击,导致包络谱中出现虚假峰。我处理的时候,习惯把重构信号的边缘10%直接舍弃,不参与后续分析,或者用hanning窗对首尾进行渐变衰减。对一段8192点的信号来说,舍掉两端各100到200个点对频率分辨率几乎没有影响,但能显著降低误判概率。
子带选择方面,TQWT分解后每个子带对应一个频带范围。判断哪些子带属于故障冲击源,不能只看子带能量,因为噪声能量分摊到多个子带后并不起眼,但真正包含冲击的子带在时域上会呈现明显的稀疏尖峰,幅值分布具有重尾特征。我常用的筛选指标是子带信号的峰度,峰度明显高于其他子带的,优先保留;峰度接近3的子带基本可以认定是噪声或平稳谐波,直接丢弃。这个做法操作简单,效果却比纯能量筛选靠谱得多。
3. 轴承故障特征频率计算与包络谱原理
3.1 四类故障特征频率计算公式
轴承故障诊断的最终目标是在频谱图上找到与轴承几何参数和转速相关的特征频率。假设滚动体个数为n,转频为fr,滚动体直径为d,节圆直径为D,接触角为α,则四类典型故障的特征频率计算公式如下。
外圈故障特征频率(BPFO):
f_BPFO = (n × fr / 2) × (1 - d × cos α / D)
内圈故障特征频率(BPFI):
f_BPFI = (n × fr / 2) × (1 + d × cos α / D)
滚动体故障特征频率(BSF):
f_BSF = (fr × D / (2 × d)) × (1 - (d × cos α / D)^2)
保持架故障特征频率(FTF):
f_FTF = (fr / 2) × (1 - d × cos α / D)
举一个实际算例。某深沟球轴承节径D等于39.04mm,滚动体直径d等于7.94mm,滚动体个数n等于9,接触角α等于0度,轴转速为1500r/min,转频fr等于25Hz。代入公式可以得到:
f_BPFO = (9 × 25 / 2) × (1 - 7.94/39.04) = 89.61Hz
f_BPFI = (9 × 25 / 2) × (1 + 7.94/39.04) = 135.39Hz
f_BSF = 25 × 39.04 / (2 × 7.94) × (1 - (7.94/39.04)^2) = 42.38Hz
这些数值就是后续包络谱判读的“对照表”。需要注意的是,实际轴承在运行中会存在轻微打滑,测量得到的特征频率与理论值的偏差通常在1%到2%以内,超过3%的基本可以怀疑参数录入或者转速测量出了问题。
3.2 包络谱为什么比原始频谱有效
轴承故障冲击的频谱特征是宽带共振峰被故障特征频率调制,这在原始频谱上体现为共振峰两侧等间距的边带。问题在于边带间距就是特征频率,通常在几十到一百多赫兹,而边带本身幅值很小,很容易被转频谐波的泄漏边缘掩盖。直接对原始信号做FFT,低频段的谱线拥挤在一起,很难分辨出边带间隔。
Hilbert变换解调的本质是提取信号的瞬时幅值包络。对故障冲击信号来说,包络信号的频率成分就直接反映了冲击的重复频率,也就是故障特征频率。做完Hilbert包络后再做FFT,特征频率处的峰值就变得非常突出,而且谐波和边带关系一目了然。TQWT重构后的信号已经大幅抑制了噪声和转频成分,包络谱的底噪更低,特征频率峰值更尖锐,两者搭配属于典型的“1+1>2”组合。
我实际测试过多种组合方式,结论是:直接用原始特征频率带通过滤加包络谱,效果受限于带通范围选择,选偏了就漏检;TQWT重构加包络谱,对带通范围不敏感,鲁棒性明显更好。这也是为什么把这个流程作为项目核心方案。
4. MATLAB r2021b代码实现全流程
4.1 仿真信号生成与参数设置
为了验证算法流程,我习惯先用仿真信号跑通,再加真实数据。仿真信号包含三部分:转频谐波成分、周期冲击成分和高斯白噪声。转频谐波模拟轴的旋转;冲击成分模拟轴承内圈故障,冲击周期对应内圈故障特征频率,冲击本身是高频衰减振荡;白噪声模拟环境干扰。
fs = 12800; % 采样率 12.8kHz N = 8192; % 采样点数 t = (0:N-1) / fs; % 时间序列 fr = 25; % 转频 25Hz BPFI = 135.39; % 内圈故障特征频率 % 转频谐波成分 x_harmonic = 0.6 * sin(2*pi*fr*t) + 0.3 * sin(2*pi*2*fr*t) + 0.15 * sin(2*pi*3*fr*t); % 周期性冲击成分,每1/BPFI秒出一次冲击 T_imp = 1 / BPFI; imp_idx = round(1 : T_imp*fs : N); x_imp = zeros(1, N); for k = 1:length(imp_idx) idx = imp_idx(k); if idx + 100 <= N tt = (0:100) / fs; x_imp(idx:idx+100) = x_imp(idx:idx+100) + 0.5 * exp(-tt * 800) .* sin(2*pi*3000*tt); end end % 叠加噪声 rng(42); x = x_harmonic + x_imp + 0.1 * randn(1, N);这段代码有几个细节值得注意。冲击响应用指数衰减乘以3000Hz正弦波模拟,衰减系数800决定冲击持续大约5毫秒,和真实传感器共振衰减特性接近。噪声幅值0.1是经过多次测试选定的,此时冲击成分在时域上还能隐约看到,但直接做FFT已经找不到故障特征频率,符合“早期故障微弱特征”的测试目标。
4.2 TQWT分解与MCa重构实现
MATLAB r2021b没有内置TQWT函数,但实现起来并不复杂。核心是设计两个缩放因子滤波器,然后迭代分解。Selesnick官网提供了官方实现,如果不方便下载,也可以按下面的逻辑自己写一个精简版。
function [w, f] = tqwt(x, Q, r, J) beta = 2 / (Q + 1); alpha = 1 - beta / r; w = cell(1, J+1); temp = x(:)'; for j = 1:J % 设计低通和高通滤波器 [lo_d, hi_d] = tqwt_filter(alpha, beta, length(temp)); low = conv(temp, lo_d, 'same'); high = conv(temp, hi_d, 'same'); % 下采样 w{j} = high(1:2:end); temp = low(1:2:end); end w{J+1} = temp; end这里的tqwt_filter函数需要根据α和β计算频域滤波器响应,核心公式是低通响应为0到α频段保持1,α到β频段余弦渐变;高通响应为β到π频段保持1,α到β频段补充。源码实现细节较多,正式项目中直接参考Selesnick的公开代码更稳妥。
分解出J+1个子带后,形态分量分析的目标是重建包含冲击成分的信号。最简单有效的方法是对所有高频子带计算峰度,保留峰度大于3且冲击特征明显的子带,其余置零,然后逆TQWT重构。这个做法虽然不如压缩感知框架下的稀疏优化精致,但在工程上足够可靠。我最初尝试用L1范数优化求解MCa,迭代慢且参数敏感,换成峰度选带后分析速度提升了几十倍,故障识别结果几乎没有区别。
4.3 包络谱与故障判别
重构出冲击信号后,用Hilbert变换提取包络,再对包络做FFT,就得到包络谱。
env = abs(hilbert(x_reconstructed)); f = linspace(0, fs/2, floor(N/2)); spec = abs(fft(env)); spec = spec(1:floor(N/2)); spec(1:5) = 0; % 去除零频直流分量 % 查找特征频率附近的峰值 target = BPFI; search_band = 5; [~, idx_target] = min(abs(f - target)); peak_range = (idx_target-search_band):(idx_target+search_band); peak_val = max(spec(peak_range));故障判读的逻辑是:在理论特征频率加减5Hz范围内搜索峰值,如果该峰值明显高于包络谱平均底噪(比如超过3倍底噪均值),就判定该频率存在。更严格的做法是同时检查2倍频和3倍频处是否存在峰值,因为早期故障的包络谱往往在基频处并不突出,但倍频序列非常规整,多倍频参照能显著提升确认置信度。
从这个角度来看,包络谱的判读不能只看幅值最大的峰,而要看频率轴上是否存在“等间隔的峰值序列”。这就是我后文处理真实信号时最依赖的判断准则。
5. 实测案例:内圈故障信号分析全过程
5.1 实验数据背景
某次项目实测,对象是一台电机的驱动端轴承,型号6205深沟球轴承,转速1440r/min,采样率12800Hz,采样点数8192。轴承参数为节径39.04mm,滚珠直径7.94mm,滚珠数9,接触角0度。转频24Hz,计算得到内圈故障特征频率约131.6Hz,外圈故障特征频率约84.4Hz。
采集到的原始信号时域波形基本看不出明显冲击,幅值在±0.3g范围内波动,频谱在50Hz、100Hz、150Hz处有明显转频谐波峰值,在2000Hz以上有一片宽频隆起,但没有清晰边带。这种信号如果交给没有经验的人判断,很容易得出“轴承正常”的结论,但设备已经出现间歇性异响。
5.2 参数整定与分解结果
针对这段信号,我设置双通道TQWT参数:低Q通道Q=5,r=3,J=30,负责提取冲击;高Q通道Q=20,r=3,J=30,负责建模转频谐波和持续振荡。两个通道分别分解后,低Q通道重构信号和高Q通道重构信号相减,得到冲击候选成分。
观察重构信号的时域波形,可以明显看到一系列稀疏尖峰,尖峰间隔大致为7.6毫秒,对应131.6Hz的倒数。这个发现直接说明内圈故障特征频率在时域上已经暴露出来。再看各个子带的峰度分布,前10个子带的峰度值集中在4到6之间,明显高于噪声子带的峰度2到3,选带范围很明确。
这边有一个实际经验:TQWT的分解层数J不是越大越好。我曾把J从30加到50,结果低频子带被压缩得厉害,重构时引入了额外振荡,包络谱里反而多出很多虚假峰。后来我按最低层子带中心频率不低于故障特征频率2倍来反推J,效果最稳定。对这段数据,J=30已经是合理上限附近。
5.3 包络谱判读与结论
重构信号做Hilbert包络谱后,在131.7Hz位置出现明显峰值,幅值约为底噪的6倍,符合内圈故障特征频率理论值131.6Hz,偏差小于0.1%。同时,263.4Hz处出现二次谐波,395.1Hz处出现三次谐波,谐波序列清晰。更关键的是,在131.7Hz附近还能看到以转频24Hz为间隔的边带,这是内圈故障的典型特征——内圈随轴旋转,故障点相对负载区的位置周期性变化,导致冲击幅值被转频调制。边带的存在让诊断结论几乎不需要犹豫。
同步观察外圈故障特征频率84.4Hz处的包络谱幅值,只有底噪的1.2倍,不构成诊断依据。最终结论明确为驱动端轴承内圈滚道故障。后拆检发现内圈滚道有一处约1.5mm的剥落坑,与诊断结论吻合。
这个案例让我特别想强调的一点:包络谱判读时,基频峰值加倍频序列加边带特征,三者同时出现才是高置信度判据,只看到一个孤立峰,先别下结论,继续观察或者换工况复测。
6. 常见问题与排查技巧实录
6.1 典型问题速查表
以下问题是我在大量实验和现场数据中反复遇到的,整理成速查表方便大家对照。
| 症状 | 可能原因 | 解决方法 |
|---|---|---|
| 包络谱无任何特征频率峰值 | Q值过低,冲击子带被稀释 | 提高低Q通道Q值至6-8,同时检查分解层数 |
| 特征频率峰值存在但很低 | 子带选择遗漏了真实冲击频带 | 扩大峰度筛选范围,或改用L1稀疏优化重构 |
| 包络谱出现大量虚假峰 | 边界效应未处理,重构信号两端畸变 | 舍弃重构信号边缘10%再分析 |
| 特征频率与理论值偏差大于3% | 转速测量不准确或轴承存在打滑 | 用转频谐波峰值修正实际转速,重新计算特征频率 |
| 不同时间段包络谱结果差异大 | 转速波动或载荷变化导致冲击周期不稳 | 分时间段分别分析,结合工况参数综合判断 |
| TQWT计算速度过慢 | r设置过大或J设置过高 | r保持3,J根据信号长度下调至20-25 |
6.2 两个值得说的坑
第一个坑是MCa重构时过度追求“干净”导致的漏检。刚开始做这个项目时,我以为重构信号里冲击越稀疏、底噪越低越好,于是在MCA的稀疏优化里把正则化参数调得很大。结果重构出来的信号只剩两三个极大冲击,包络谱里特征频率峰值反而消失了。原因在于轴承故障冲击的幅值并不是恒定不变的,内圈故障的冲击幅值会随转频周期性调制,过强的稀疏约束把所有小幅值冲击都清零了,周期性信息随之丢失。后来我把正则化系数降低了一个量级,保留了调制周期内的多级冲击,特征频率才重新暴露出来。这里我得到的教训是:诊断目标是“周期性”,不是“稀疏性”,重构时保周期比保稀疏更重要。
第二个坑是高Q值参数对冲击信号的“平滑化”。有一段时间我把Q值从5一路调高到12,认为Q越高频率选择性越好、信噪比越高,结果冲击成分被过度平滑成一段段振荡衰减波,包络谱的峰值变得很钝。Q值本质上决定了小波基对信号形态的匹配程度,故障冲击是瞬态突变事件,Q值超过8之后,小波基更像振荡波形,反而与冲击形态失配。对不同轴承实际测试下来,Q=4到6这个区间对大多数故障冲击的匹配效果最好。参数不是越“高级”越好,匹配信号物理特性才是唯一的调参准则。
6.3 参数自动搜索思路
手动调参在处理单组数据时可行,但面对一批不同工况、不同轴承类型的数据时,每次都手动调一遍Q和J效率太低。我在项目后期做了一个简单的网格搜索:Q从3到7以1为步长,J从15到35以5为步长,对每组参数组合做TQWT分解、MCA重构、包络谱分析,以理论特征频率处的包络谱峰度作为目标函数,选最大峰度对应的参数作为最优参数。一组8192点信号大概测试40组参数组合,总耗时不到2分钟,完全可接受。
更聪明的做法是先用Kurtogram快速定位共振频带,再用共振频带的中心频率和带宽反推TQWT的Q因子初值。Kurtogram本身是快速峭度图,能找出信号中最具冲击特性的频带,虽然没有直接给出Q值,但提供了一个很好的起点。这个方法在工程中非常实用,适合作为TQWT参数自动整定的前置环节。
结尾
在设备故障诊断这条路上,TQWT不是万能的,但它确实解决了传统方法很难兼顾的一个矛盾:既要保留瞬态冲击的细节,又要尽可能压制持续谐波和噪声。从我个人的使用体验来说,TQWT更像一个“按波形形态分拣”的工具,它不要求你知道故障频率在哪,只需要你对信号中冲击成分的形态有大致判断,剩下的让参数去匹配。实际项目里,我更习惯把TQWT和包络谱当作第一道筛选手段,确认有异常后再结合温度、油液、声音等多维数据综合判断,最终才给出维修建议。如果你正在被“频谱图一团糊、特征频率找不到”的问题困扰,不妨按这个流程试一试TQWT,也许会有意料之外的收获。最后再分享一个小技巧:不管分析结果有多漂亮,诊断结论一定要留有余地,标注“建议复测确认”永远比绝对化的“确诊故障”更能经得起现场考验。