变转速轴承故障诊断：带通滤波+倒谱预白化+平方包络谱实战

去年在某风电齿轮箱的测试台架上我遇到过一件挺头疼的事转速在600到900转之间反复横跳轴承外圈故障的特征频率在包络谱里找不着能量被摊成一片。后来把信号按角度重采样再用带通滤波加倒谱预白化最后看平方包络谱故障阶比才像钉子一样扎出来。这篇文章就围绕带通滤波、倒谱预白化、平方包络谱这套组合专门聊它们在变速条件下做轴承故障检测时为什么管用、怎么落地以及Matlab代码怎么写。内容偏实战适合搞旋转机械状态监测、风电传动链测试、机床主轴诊断的朋友尤其是那些转速不稳定、包络谱一团浆糊的工况。1. 变速条件为什么让传统包络谱失效1.1 恒定转速下的老套路固定转速下轴承故障诊断有一个几乎“无脑”的流程取一段振动信号带通滤波到轴承共振频带做Hilbert包络再FFT得到包络谱在故障特征频率处找峰。这个流程依赖一个前提——故障特征频率是固定的。比如一个外圈故障阶比为3.57的轴承转频fr10Hz时BPFO35.7Hz转速不变35.7Hz处就会有一个稳定的谱峰。你只需要把谱图拉出来对着故障频率表看一眼基本就能下结论。这种方法的有效性建立在信号平稳的假设上。包络谱本质上是时域包络信号的傅里叶变换它隐式地假设包络中的周期性成分在整个分析时间内以固定频率重复。对恒定转速的齿轮箱、电机、风机这个假设基本成立。所以过去很多老师傅都习惯用这套流程做日常巡检。1.2 转速一波动谱线就“融化”了一旦转速开始变化事情就变麻烦。假设同一只轴承转频从10Hz慢慢爬到15HzBPFO也会从35.7Hz爬到53.6Hz。如果你把这段10秒的数据直接做FFT这18Hz的跨距会把谱峰能量摊到几十条谱线上原本应该尖锐的峰变成一座小土包高度可能只有固定转速时的五分之一甚至更低。这还只是单一变化实际现场不只转速波动齿轮啮合频率、边带、轴转频谐波都在跟着转速漂移它们从“离散的谱线”变成“连续宽带干扰”跟轴承冲击在频谱上更难区分。我最早遇到这种情况时第一反应是把FFT窗口缩短让每一段近似平稳。但窗口缩短带来的问题是频率分辨率下降包络谱上一个峰可能对应5Hz的宽度很多低速轴承的故障特征频率本来就只有几十赫兹窗口短了根本分辨不出来。加窗、重叠也没什么本质改善因为问题不在谱估计而在信号本身非平稳。1.3 角域重采样把变转速问题变成定转速问题解决非平稳的标准思路是把“时间域”换成“角度域”。只要你有同步的转速信号键相脉冲、编码器或转速计就能计算出瞬时转速曲线再积分得到累计转角。把振动信号按照“每转固定采样多少个点”重新插值一遍就得到了角域信号。角域信号的特点是转轴每转一圈对应固定数量的样本因此原本随转速变化的故障特征频率在角域中变成了固定的“特征阶比”。特征阶比的物理意义是“每转发生多少次故障事件”。以外圈故障为例BPFO_order z/2 * (1 - d/D * cosα)这个式子只跟轴承几何尺寸有关跟转速无关。角域重采样之后变速工况下的轴承故障检测就退化成恒定转速下的“固定频率检测”只是频率轴的单位从Hz变成了orderscycles/revolution。这是整个方法能够成立的地基。2. 倒谱预白化的原理与在轴承诊断中的定位2.1 用生活化方式理解倒谱倒谱Cepstrum这个名字本身就是把“Spectrum”前面的字母倒过来暗示“频谱的频谱”。概念上确实有点绕但可以这样理解如果频谱里有一组等间距的梳状谱线那么这种等间距的周期结构会在倒谱中形成一个尖峰尖峰的横坐标就是谱线间距的倒数。打个比方你在一个喧闹的食堂里录一段音频频谱能看到许多规则排列的人声泛音倒谱则能告诉你这些泛音的“周期性”来自哪里。在旋转机械振动信号中最容易在倒谱里制造尖峰的就是转频谐波、齿轮啮合频率及其高阶边带。它们共同的特点是在频谱中呈现等间隔的“梳齿”结构。而轴承故障冲击不一样它激起的是结构共振频谱上是宽带大包络没有明显的等间隔梳齿。所以倒谱天然地可以把“确定性周期成分”和“冲击类非周期成分”区分开来。2.2 预白化到底“白”在哪里倒谱预白化的具体操作一句话就能说完把倒谱中除了零倒频率以外的系数全部清零然后用这个被“破坏”的倒谱重构一个幅度谱这个新幅度谱是常数——也就是“白色”谱再保留原始相位做逆变换。因为所有频点的幅度被强制相同任何频率上原有的能量优势都消失了。齿轮啮合谐波、转频谐波这些曾经在频谱中高高凸起的谱线全部被抹平。你可能会问把幅度谱全抹平了那轴承冲击也变成白色谱了还有什么信息关键在于相位谱没有动。冲击在时域中的位置、发生时刻都编码在相位里。把幅度谱白化后逆变换回时域得到的信号保留了“什么时候有冲击”这个信息却删掉了“哪些频率有稳定大能量”这个信息。对于包络谱分析来说我们最终关心的是包络的周期性也就是冲击发生的间隔而不是冲击的形状和大小。预白化恰好把影响包络周期判断的确定性谐波干扰去掉了同时把噪声也“漂白”成与信号同量级的成分再经过平方包络周期性冲击信息就更容易浮现出来。2.3 为什么顺序是“带通滤波后再预白化”标题里的顺序是带通滤波后的倒谱预白化这不是随便排的。如果对原始宽频信号直接做倒谱预白化低频转频的巨大能量和高频噪声都会被强制拉到同一量级白化后整个频带的信噪比并没有提升甚至可能把原本藏在微弱噪声里的冲击淹没掉。先做带通滤波把频带限制在轴承共振区域附近再去白化相当于在已经提高信噪比的频带内做“公平竞争”。这样冲击频带内原本相对较低的离散干扰被剔除冲击就成了主导成分。我做对比实验时发现如果不带通滤波直接白化平方包络谱的底噪会明显抬高故障阶比的峰几乎被底噪淹没而先带通到共振频带峰高与底噪的比值能提升一个数量级以上。这个顺序值得大家在实际处理时坚持。3. 从原始振动到平方包络谱的完整流程3.1 需要什么样的数据和转速信息这套方法需要两个同步数据振动加速度信号和转速信号。振动信号最好从靠近被测轴承的测点采集采样率要能覆盖轴承共振频带一般不低于20kHz。转速信号可以从键相槽、编码器或者光电转速计获取至少要能给出每转一次的脉冲信号。没有转速信号角域重采样就无从谈起。有的人问“我只有振动数据没有转速能做吗”能但要先估计瞬时转频。比较粗糙的办法是用短时傅里叶变换在低频段找转频的轨迹或者用相位解调方法从轴频谐波中提取。这类方法的问题在于转速波动大的时候估计误差也大重采样之后阶比谱会出现虚假边带。我的建议是现场如果条件允许一定要同步采集转速省得后面处理时怀疑结果。3.2 带通滤波的频带怎么选带通滤波的频带一般是轴承故障冲击激发的结构共振频带不是随便选的。最正规的方法是快速峭度图Kurtogram它会在时频图上自动找出峭度最大的频带和中心频率。Matlab有现成函数kurtogram需要信号处理或额外工具箱也可以用“滑动窗口计算包络谱人工比较不同频带效果”这种笨办法。几个实际经验带宽不要低于共振峰的半宽否则冲击的边带会被削掉带宽也不要太宽否则白化时噪声受益。我一般先看原始功率谱找高频段有没有集中的“山包”以山包为中心取2kHz到4kHz带宽。对于风电齿轮箱8k到16k是常见选择对于小型的机床主轴共振频率可能更高。总之先多试几个候选带再做决定。3.3 角域重采样的实施细节角域重采样实质上是改变信号的“钟表”坐标系从均匀时间间隔改为均匀角度间隔。实现分四步第一步从转速脉冲算出瞬时转频对转频做平滑处理第二步积分得到累计转角φ(t)第三步设定每转采样点数Ns生成等角度间隔序列φ_k第四步用插值把原始信号映射到φ_k对应的采样点。Ns的取值需要遵循采样定理角域信号的最大阶比不得高于Ns/2。如果带通滤波后的最高分析频率是16kHz而最低转速是500rpm那么最高可能出现的阶比约等于16000/(500/60)1920因此Ns至少要取3840实际取4096更稳妥。取小了会在低速段发生高频混叠包络谱里出现一堆莫名其妙的假峰。插值方法建议用三次样条线性插值会产生额外高频分量导致包络谱中出现虚假峰值。3.4 倒谱预白化与平方包络谱的组合逻辑信号流入的顺序是角域信号 → 倒谱预白化 → 平方包络谱。这里倒谱预白化的输入是已经位于角域的信号它的横轴是角度所以倒谱处理的是“角域频谱”也就是阶比谱。离散成分转频谐波、齿轮啮合阶次在阶比谱中等间隔排列倒谱预白化可以精确地抹掉它们。然后对预白化信号做Hilbert包络、平方、FFT得到平方包络谱横轴是阶比。最后对照轴承故障特征阶比表在BPFO、BPFI、BSF等位置找峰。4. Matlab代码实现与逐段解读4.1 主程序框架一句话串起整个流程我习惯把主流程写得非常短只留清晰的步骤具体运算放函数里。下面这段代码是完整流程的主干输入是原始时域信号x、时间向量t和同步转速rpm输出是平方包络谱。% 主流程: 带通滤波 - 角域重采样 - 倒谱预白化 - 平方包络谱 fs 40000; % 原始采样率, Hz t ...; % 时间向量, s x ...; % 振动加速度信号 rpm ...; % 同步转速曲线, rpm, 与t等长 % 1) 时域带通滤波 f_low 8000; f_high 16000; [b, a] butter(4, [f_low f_high]/(fs/2), bandpass); x_bf filtfilt(b, a, x); % 2) 角域重采样 samples_per_rev 4096; % 每转采样点数, 需要根据最低转速和最高分析频率计算 [x_angle, theta] resampleToAngle(x_bf, t, rpm, samples_per_rev); % 3) 倒谱预白化 xw cepstrumPrewhiten(x_angle); % 4) 平方包络谱 [ses, order] squaredEnvelopeSpectrum(xw, samples_per_rev); % 5) 画图 figure; plot(order, ses, b-, LineWidth, 1.2); xlabel(Order (cycles/revolution)); ylabel(Squared envelope spectrum); xlim([0 30]); grid on;注意带通滤波用的是butter4阶设计、filtfilt零相位滤波避免相位偏移影响后续冲击位置。如果信号很长可以分块读取但要注意重采样和倒谱都是整体操作分块会破坏连续性。4.2 倒谱预白化函数核心三行倒谱预白化的完整实现只有几行核心是“对数量幅度谱→IFFT→倒谱系数置零→FFT重建”。下面是函数输入一个向量输出白化后的向量。function xw cepstrumPrewhiten(x) X fft(x); logAmp log(abs(X) eps); cep ifft(logAmp); cep(2:end) 0; % 只保留零倒频率, 其余清零 whiteLogAmp fft(cep); % 常数对数谱 whiteAmp exp(whiteLogAmp); % 常数幅度谱 Xw whiteAmp .* exp(1i * angle(X)); xw real(ifft(Xw)); end这里有几个细节。加eps是为了防止log(0)出现NaN。cep(2:end)0意味着只留下倒谱中的直流分量也就是对数幅度谱的平均值重组成的白色幅度谱对所有频率都取同一个放大倍数。乘上原始相位再逆变换得到时域/角域的白化信号。这个函数对输入长度没有强制要求但建议长度是2的幂FFT更快。4.3 平方包络谱函数直接可用平方包络谱与普通包络谱的区别在于先对包络做平方再做FFT。平方操作相当于对包络信号进行自相关增强周期性冲击的基频和谐波会成倍突出。函数如下横轴直接输出阶比。function [ses, order] squaredEnvelopeSpectrum(x, samples_per_rev) N length(x); env abs(hilbert(x)); sq env.^2; win hann(N, periodic); S fft(sq .* win); S S(1:floor(N/2)1); ses abs(S).^2 * 4 / N; order (0:floor(N/2)) / N * samples_per_rev; endhilbert是Matlab信号处理工具箱的函数返回解析信号取模就是包络。加汉宁窗是为了减少泄漏如果你在角域信号中已经感觉到泄漏问题可以试试flattop窗。ses的归一化4/N是经验值谱峰的绝对高度没有意义关键是相对周围的底噪高度。4.4 角域重采样函数注意的坑重采样的核心是计算累计转角和插值。下面这个函数是最简实现适合转速单向递增的场合。function [x_angle, theta] resampleToAngle(x, t, rpm, samples_per_rev) fr rpm / 60; % 转每秒 phi 2 * pi * cumtrapz(t, fr); % 累计转角, rad phi phi - phi(1); total_rot phi(end) / (2 * pi); n round(total_rot * samples_per_rev); theta linspace(phi(1), phi(end), n); x_angle interp1(phi, x, theta, spline); end最大的坑是转速曲线有噪声。如果rpm抖动厉害cumtrapz之后的角度路径会有微小台阶插值会引入虚假高频成分。建议在重采样前用smooth或medfilt1对rpm做平滑。另外如果转速有瞬间反转比如往复机械phi不再是单调函数interp1会报错需要分段处理。风电齿轮箱、离心机这类单向旋转的设备一般没问题。4.5 模拟信号验证证明整套流程有效为了验证上述代码我写了一个简单的变速轴承信号生成器。它生成一段带转速波动的时域信号包含轴承外圈故障冲击、齿轮啮合阶次干扰和随机噪声。生成代码的思路是用角度域的相位来定义齿轮谐波再映射回时间域。fs 40000; t 0:1/fs:8; % 8秒数据 rpm 700 200*sin(2*pi*0.25*t); % 转速约500~900rpm fr rpm/60; phi 2*pi*cumtrapz(t, fr); % 累计转角 % 外圈故障阶比 3.57 BPFO_order 3.57; delta_phi 2*pi/BPFO_order; n_impulse floor(phi(end)/delta_phi); th interp1(phi, t, (0:n_impulse)*delta_phi, linear); % 冲击时刻 x_bearing zeros(size(t)); for k 1:length(th) idx t th(k); x_bearing(idx) x_bearing(idx) ... 0.8*exp(-600*(t(idx)-th(k))) .* sin(2*pi*12000*(t(idx)-th(k))); end x_gear 2*sin(23*phi) 0.8*sin(46*phi); % 23、46阶齿轮谐波 x_noise 0.5*randn(size(t)); x x_bearing x_gear x_noise;把这段生成代码放到主程序之前然后运行整个流程理论上会在3.57阶、7.14阶等位置看到谱峰。我自己跑下来3.57阶的峰高大约是底噪的8到10倍。如果你不做角域重采样直接用原始x做平方包络谱则几乎看不到3.57阶的峰因为23、46阶齿轮谐波和变速带来的能量扩散把低阶冲击峰盖住了。5. 实测案例变转速轴承数据的效果对比5.1 案例数据描述这里用上一节的模拟信号作为“实测”输入持续8秒转速在500到900rpm之间正弦波动。轴承外圈故障阶比3.57这意味着每转出现3.57次冲击。同时叠加了23阶和46阶的齿轮啮合干扰噪声标准差0.5。采样率40kHz每转采样点数4096带通滤波范围8k-16kHz。5.2 直接包络谱与本文方法的结果对比我把两种处理路径做了对比。路径A是传统方法直接对原始时域信号做带通滤波、Hilbert包络、FFT横轴换算成阶比假设转速已知。路径B是本文方法带通滤波-角域重采样-倒谱预白化-平方包络谱。两者都只看3.57阶附近。结果如下表所示。处理路径3.57阶谱峰高度底噪水平峰噪比近似结论A: 直接带通包络谱0.420.331.3峰被噪声覆盖难以判定B: 带通重采样预白化平方包络谱1.850.218.8特征阶比清晰可辨路径A的问题不只是峰低而是23阶齿轮干扰和转速波动带来的“宽带高底”让谱峰完全失去特异性。路径B先通过角域重采样把所有特征固定到恒定阶比再用倒谱预白化把23阶、46阶的确定成分抹平最后平方包络把3.57阶的周期性冲击放大结果自然清晰。这张表里的具体数字因参数而异但趋势是非常稳定的。5.3 为什么先带通滤波能提升稳定性在这个案例里如果我把步骤改成“先角域重采样再倒谱预白化最后带通滤波”效果会差一些。原因很简单全频段白化会把随机噪声和齿轮谐波同时拉平而8k-16k之外的噪声大量涌入3.57阶的峰会被底噪吞噬。先带通做了一次“安检”只让共振频带的信号进入后续处理这样预白化的去干扰效果才会被放大。如果你用快速峭度图选频带效果通常比固定频带更稳。6. 工程应用中的坑与经验建议6.1 倒谱预白化会放大哪些噪声预白化不是万能的。它把所有频率的幅度拉到同一水平这也就意味着原本幅度很低的随机噪声会被放大。如果噪声在整个频带上都很大或者带通滤波带宽选择过宽预白化后的信号底噪可能比原始信号还高。我个人的经验是白化前一定要有一个带通滤波步骤而且带宽宁窄勿宽。白化的价值是“去除确定性离散成分”不是“提高随机噪声环境下的通用信噪比”。此外平方包络谱的幅值会随着白化增益变化不能拿两段不同转速、不同负荷的数据直接比较幅值来判断故障恶化。需要建立基准库或者观察特征阶比幅值与整体底噪的比值趋势。6.2 带通滤波器参数别照搬很多人拿到代码后直接把8k-16k照抄过去这不对。轴承共振频带由轴承结构、传感器安装、信号传输路径共同决定。同样一只轴承装在不同设备上共振频带可能完全不同。我建议在选频带时做几个候选频带的对比看同一段数据在不同频带下哪个3.57阶峰更突出。如果没有时间做快速峭度图就根据功率谱里高频峰群手动选一个窄带试试比盲目用固定带好得多。6.3 角域重采样的转速曲线要平滑转速曲线不平滑是导致重采样后信号质量差的最常见原因。转速计脉冲间隔天然存在抖动直接转换成转速后会有高频毛刺。这些毛刺会让累计转角路径出现微小抖动插值后的角域信号会带有人为的频率调制进而在平方包络谱里产生边带。建议在重采样前对转速曲线做滑动平均或中值滤波时间常数取一个转周期左右。如果你用的是“转速计每转一个脉冲”那转速曲线的分辨率只有一转需要插值平滑到每时刻否则重采样结果会非常粗糙。6.4 诊断阈值与特征阶比的匹配在实际诊断时不要只看理论阶比处有没有峰要看三个指标特征阶比处的峰值是否超过周围底噪的3到5倍是否在2倍、3倍谐波处也有等间隔峰是否同时出现边带。轴承故障阶比会因为滚动体打滑、保持架变形等因素略有偏移允许±0.05阶的容差。另外最好同时算BPFO、BPFI、BSF、FTF四个阶比的谱峰再做组合判断。只盯一个阶比很容易被偶然干涉峰骗到。最后说句实在话这套方法解决的是变转速下离散干扰压过轴承冲击的问题但它不是银弹。轴承严重磨损、润滑不良引发的宽频噪声升高倒谱预白化是处理不了的必须结合时域峭度、温度、油液金属颗粒这些指标一起看。工具只是帮你把特征从噪声里捞出来真正下结论还是要靠对设备工况的理解。