先聊一个很多做设备状态监测的朋友都绕不过去的场景现场振动数据一测一大把可一遇到变速工况传统频谱分析就抓瞎——转频在飘故障特征频率也跟着飘固定频点上一会儿有峰一会儿没峰整个谱图看着就像一团乱麻。我早年做风电齿轮箱和轧机轴承诊断时被这种时有时无的特征搞得头大过不少次。后来项目里逐步用上一套组合方案效果非常稳定带通滤波 倒谱预白化 平方包络谱英文缩写就是CPWCepstral Pre-whitening加SESSquared Envelope Spectrum。这套东西不挑转速不需要你额外装编码器也不用先做阶次跟踪直接对原始振动信号动手就能把淹没在齿轮啮合、结构共振、随机噪声里的轴承微弱冲击给捞出来。这篇文章就把我实际跑通的完整思路、原理、Matlab代码和踩过的坑一次讲清楚适合正在做轴承故障诊断研究、设备状态监测系统开发或者看论文看到CPW但不知道怎么落地的朋友。1. 项目背景变速工况下轴承故障检测为什么难1.1 轴承故障检测的基本逻辑滚动轴承故障诊断的底层逻辑很简单轴承内外圈、滚动体出现局部损伤剥落、点蚀、裂纹后滚动体经过损伤点时会产生周期性冲击。这种冲击会激起轴承座和传感器的结构共振形成一个高频衰减振荡。诊断这件事本质上就是想办法把这个微弱的周期冲击从振动信号里识别出来并算出它对应的特征频率。轴承故障特征频率有四个固定公式这里先摆出来后面代码要用外圈故障频率BPFO( BPFO \frac{n}{2} f_r \left(1 - \frac{d}{D} \cos \alpha\right) )内圈故障频率BPFI( BPFI \frac{n}{2} f_r \left(1 \frac{d}{D} \cos \alpha\right) )滚动体故障频率BSF( BSF \frac{D}{2d} f_r \left(1 - \left(\frac{d}{D} \cos \alpha\right)^2\right) )保持架故障频率FTF( FTF \frac{1}{2} f_r \left(1 - \frac{d}{D} \cos \alpha\right) )其中 ( n ) 是滚动体个数( d ) 是滚动体直径( D ) 是节圆直径( \alpha ) 是接触角( f_r ) 是转频。在恒定转速下这些频率是固定的直接用FFT幅值谱找峰值就行这是教科书里最基础的套路。但在现场转速总是波动的设备启停机、负载变化、调速运行转频一会儿30Hz一会儿50Hz故障特征频率成比例地跟着变。固定频率的谱峰被抹平在很宽的频带里甚至完全湮没在噪声中传统的包络谱分析效果很不理想。1.2 传统方法在变速工况下的失效原因变速工况下传统方法失效的原因可以拆成三层来看第一层转速变化导致特征频率非平稳。故障特征频率是转频的固定比例倍数转频一变特征频率就跟着变。做长时FFT时能量被分散到多个频点上峰值幅度被稀释信噪比下降。第二层振动信号里干扰成分复杂。齿轮啮合频率及其谐波、边带、结构共振、随机冲击混在一起。尤其齿轮箱里齿轮啮合能量通常比轴承冲击能量大一个数量级以上。变速时齿轮啮合频率也在变和轴承微弱信号搅在一起很难靠简单的频带选择分离。第三层经典的包络谱方法依赖共振频带选择。传统做法是先用带通滤波器把某个共振频带选出来再做包络谱。可是现场没有先验知识不知道轴承冲击激起了哪个共振峰。选宽了噪声多选窄了漏掉冲击成分。在变速工况下共振频带虽然相对固定但特征频率是漂移的如果滤波器中心频率和带宽没选好包络谱照样一塌糊涂。1.3 倒谱预白化在变速工况下的价值定位面对这些难题学术圈和工业界给出过很多解法阶次跟踪、时频分析、同步平均、盲解卷积MED、MCKD、谱峭度等。各有各的适用范围。**倒谱预白化CPW**属于近年来逐渐被重视的方法它在变速工况下有一个非常关键的优势——信号处理过程完全不依赖转速信息。CPW的思路是把振动信号在频域的对数幅值谱转换到倒谱域然后把倒谱中代表慢变结构分量比如齿轮啮合、共振调制、传输路径效应的低倒频率成分置零或削弱再还原到时域。这一顿操作之后信号里的确定性周期成分齿轮成分、轴频谐波被大幅压制残留的主要是非确定性成分和轴承冲击成分相当于对信号做了白化让后续的包络谱更容易凸显轴承故障特征。这种方法论上的含义是无论转速怎么变只要轴承有冲击CPW都能把冲击从确定性干扰中剥离出来。它不像阶次跟踪那样需要转速计也不像盲解卷积那样需要迭代优化和参数调优计算效率高、鲁棒性好特别适合工业现场数据。2. 核心原理拆解三条链路串起来看这套方案本质上是三个处理阶段的串联先用带通滤波把最能代表轴承冲击的共振频带抓出来再用倒谱预白化把残余的确定性干扰按在倒谱域里干掉最后用平方包络谱把冲击的周期性特征搬到频谱上找特征频率。下面逐个拆开讲原理和设计理由。2.1 倒谱预白化到底做了什么要理解CPW先得理解倒谱。倒谱的概念最早是Bogert等人1963年提出的用于回波检测。它的核心操作是对信号的功率谱取对数再对这个对数谱做一次逆傅里叶变换。得到的新域叫倒谱域横轴的单位叫倒频率quefrency量纲是时间。为什么要取对数再做变换因为对数可以把频谱中相乘的关系变成相加的关系。振动信号里的很多干扰比如传输路径效应、结构共振调制在频域表现为频谱包络的慢变起伏而轴承冲击则表现为频谱中均匀分布的精细谐波结构。这两者在原频谱中是相乘的但在对数谱中就变成了叠加。再做一次傅里叶变换它们在倒谱域中就分布在不同的倒频率区间低倒频率靠近0代表频谱的慢变包络对应结构共振、传输路径等确定性调制高倒频率代表频谱的快速波动对应周期性冲击的精细谐波结构。所以CPW的做法非常直观把倒谱中低倒频率的成分清零高通滤波然后把剩余的高倒频率成分还原成新的频谱幅值再结合原始相位做逆傅里叶变换得到预白化后的时域信号。这个过程相当于把频谱的形状抹平了把信号变成近似白噪声激励下的响应而白化后剩下的就是非确定性的、类冲击的成分。我个人的理解类比是倒谱域里做高通滤波就像照片处理里减去低频背景光只保留高频纹理细节。轴承冲击就是纹理细节齿轮啮合和共振就是背景光。2.2 带通滤波的作用与频带选择在CPW之前加一个带通滤波工程上的意义非常大。第一减少无关频率成分提高白化的针对性。原始振动信号频带往往很宽可能从几赫兹到几十千赫兹。如果直接对全频段做CPW那些远离轴承共振频带的噪声会在白化过程中被放大干扰后续分析。先带通滤波把注意力集中到轴承冲击占主导的频带上CPW才能更有效地凸出冲击。第二降低信号的非平稳性干扰。低频段的轴频、齿轮啮合频率及其谐波分布非常密集而且随转速变化剧烈。带通滤波可以把这些低频大能量成分挡在门外让后续处理聚焦于高频冲击成分。关键问题是频带怎么选。工程上有几种思路如果目标是轴承局部损伤用谱峭度Spectral Kurtosis或快速峭度图Fast Kurtogram自动寻找使峭度最大的频带中心频率和带宽。峭度大的频带意味着冲击性最强这是最常用的方法。如果项目里已知结构共振频带比如做过锤击试验或历史数据的谱分析直接采用先验频带。如果轴转速稳定或变化不大可以用边带能量最大化的方式选频带但在大范围变速时不太靠谱。我实测下来的经验是对于滚动轴承共振频带通常集中在2kHz到20kHz之间带宽一般选原始采样率能覆盖的1/4到1/8左右。不过最稳的还是先用Fast Kurtogram扫一遍把Kurtogram峰值对应的频带拿来做带通滤波。在Matlab里可以直接用kurtogram函数第三方工具箱比如Antoni的Fast Kurtogram或者自己写短时傅里叶变换算峭度。2.3 平方包络谱为什么能看见故障特征预白化之后信号已经干净了很多但仍然是时域信号冲击的周期性还没有直接呈现。**平方包络谱SES**就是把这层周期性挖出来的关键一步。包络谱的基本原理是通过Hilbert变换构造解析信号取模得到包络信号。包络信号的物理意义是原始信号的瞬时幅值变化轴承故障冲击引起的幅值调制正好体现在包络中。平方向操作的原因有两个一个是数学上的——经典窄带包络谱也就是对包络直接做FFT得到的幅值谱在某些情况下会出现谱线幅度和特征频率幅值不成线性关系的问题。而平方包络谱在数学推导中和谱相关分析有密切联系理论上对非平稳冲击信号的检测性能更好也更容易在平方包络谱中看到故障特征频率的二次谐波、三次谐波等。另一个是工程上的——对包络做平方运算等价于加强了冲击成分的幅值调制深度让周期冲击的能量更集中到特征频率及其谐波上信噪比更高。我做仿真和实测数据对比时同样的数据平方包络谱里特征频率峰值的突出程度明显优于直接包络谱。所以这三个环节的组合逻辑是带通滤波负责划重点倒谱预白化负责打掩护的抹掉平方包络谱负责亮出真身。三条链路环环相扣少了哪一环效果都会打折扣。3. Matlab实现流程与关键代码理论说清楚之后上干货。这一节给出我实际在Matlab里跑通的完整流程包括参数如何选取、每个环节为什么要这样写、哪些地方容易踩坑。3.1 整体流程设计完整处理流程可以分成六步读取加速度振动信号设定采样率 fs对信号做快速峭度图分析或根据先验知识确定带通滤波频带 [f_lo, f_hi]设计带通滤波器对原始信号滤波得到带通信号 x_bp对 x_bp 做倒谱预白化得到预白化信号 x_w对 x_w 做Hilbert变换求包络取平方做FFT得到平方包络谱在平方包络谱中定位BPFO、BPFI、BSF、FTF及其谐波判断故障类型。流程图用文字描述就是这样原始信号 → [带通滤波] → [倒谱预白化] → [平方包络谱] → 特征频率识别。字不多条理很清楚。3.2 关键函数与Matlab代码实现下面按环节给出核心代码。由于实际项目的数据差异很大这里用示意代码配合注释关键是体现实现思路。第一步带通滤波% 参数定义 fs 25600; % 采样率单位Hz现场采集卡一般2.56k~100k N length(x); % 信号长度 % 带通频带选择——这里以快速峭度图提取结果为例 % 假设kurtogram计算得到最佳频带为 [fc - Bw/2, fc Bw/2] fc 6000; % 中心频率单位Hz Bw 4000; % 带宽单位Hz f_lo fc - Bw/2; f_hi fc Bw/2; % 设计带通滤波器 bpFilt designfilt(bandpassiir, ... FilterOrder, 8, ... HalfPowerFrequency1, f_lo, ... HalfPowerFrequency2, f_hi, ... SampleRate, fs, ... DesignMethod, butter); % 零相位滤波避免相位畸变影响包络分析 x_bp filtfilt(bpFilt, x);这里我特意用filtfilt替代filter因为零相位滤波不会让包络产生时延畸变。对于后续的包络分析相位一致性很重要这一点后面在避坑部分会再强调。第二步倒谱预白化function x_w cepstral_prewhitening(x, cut_queffrency) % CEPSTRAL_PREWHITENING 倒谱预白化 % 输入 % x - 输入时域信号建议先做过带通滤波 % cut_queffrency - 倒频率截止值单位s小于该值的成分被置零 % 输出 % x_w - 预白化后的时域信号 N length(x); % 加汉宁窗减少频谱泄漏可选但推荐 win hanning(N, periodic); xw x(:) .* win; % 计算FFT X fft(xw); X_amp abs(X); X_phase angle(X); % 避免幅值为0导致对数无穷大加一个极小值 eps_val 1e-12; log_amp log(X_amp eps_val); % 对对数幅值做IFFT得到倒谱 cep real(ifft(log_amp)); % 倒谱域高通将小于截止倒频率的成分置零 % 注意倒谱的横轴是时间0对应直流N-1对应(N-1)/fs % 截止倒频率对应样本数需要换算 cut_samples round(cut_queffrency * fs); cep_cut cep; cep_cut(1:cut_samples) 0; cep_cut(end-cut_samples2:end) 0; % 对称部分也置零 % 由处理后的倒谱重建对数幅值谱 log_amp_white real(fft(cep_cut)); % 重建白化后的幅值谱归一化处理 amp_white exp(log_amp_white); % 白化过程将原频谱幅值除以重建的慢变幅值保留相位 % 这样等同于把频谱“抹平” X_white (X_amp ./ (amp_white eps_val)) .* exp(1j * X_phase); % 逆变换回时域取实部 x_w real(ifft(X_white)); % 由于之前加了窗这里做简单幅度修正窗函数能量补偿 win_energy sum(win.^2) / N; x_w x_w / win_energy; % 去除直流分量 x_w x_w - mean(x_w); end这里需要解释几个关键细节。第一倒谱域的高通操作截止值cut_queffrency的单位是秒需要乘以采样率换算成样本数。对于旋转机械通常设置截止倒频率为 1/(2*最大特征频率) 或根据齿轮啮合周期确定。我平时的起步值是0.01秒对应100Hz的周期然后看结果调整。第二因为倒谱是对称的置零时要同时处理正半轴和负半轴对应位置否则重建出来的频谱会不对称。第三重建时用原幅值除以估计的慢变幅值就实现了白化的效果——原来幅值高的地方被压低原来幅值低的地方被抬升频谱趋于平坦。第三步平方包络谱% 对预白化信号求包络 analytic_signal hilbert(x_w); envelope abs(analytic_signal); % 平方包络 envelope_sq envelope.^2; % 去除直流后计算FFT N length(envelope_sq); win_e hanning(N, periodic); env_win (envelope_sq - mean(envelope_sq)) .* win_e; SES abs(fft(env_win)) * 2 / N; % 单边谱 % 频率轴 f_axis (0:N/2-1) * fs / N; % 绘制平方包络谱 figure; plot(f_axis, SES(1:N/2)); xlabel(频率 (Hz)); ylabel(幅值); title(平方包络谱 SES); xlim([0, min(fs/2, 1000)]); % 关注低频段高频段一般是噪声这里的几个细节要提醒包络平方后直流分量很大作图前先去直流不然低频段会被巨大的直流旁瓣淹没加窗是为了减少谱泄漏但窗口会引入旁瓣对幅度精度要求高的话也可以不加窗看你的场景。3.3 参数怎么定实测经验参数选择是这套方法落地最容易出问题的地方。我根据自己的实际测试总结一张参数参考表参数建议范围/取值说明带通滤波中心频率 fc2kHz ~ 15kHz具体用快速峭度图或共振分析确定不要盲目定带通滤波带宽 Bw1kHz ~ 8kHz太窄容易丢失冲击成分太宽噪声多带通滤波器阶数6~10阶 Butterworth阶数太高相位畸变严重太低过渡带太宽倒谱截止倒频率 cut_queffrency0.005s ~ 0.02s需大于最大分析频率周期的2倍我一般先取0.01s再调信号时长至少包含50~100次故障冲击例如BPFO80Hz采集时长至少2~3s才能有足够分辨率采样率 fs至少是最高关心频率的8~10倍处理20kHz共振频带建议不低于50kHz采样倒谱截止倒频率的物理含义要强调一下它决定了你抹掉多慢的调制成分。设得太小很多有用的冲击调制信息也被抹掉了设得太大齿轮啮合产生的边带成分清理不干净。实际调试时我会把cut_queffrency从0.002s到0.05s扫一遍看哪种参数下平方包络谱中特征频率峰值最突出然后定下来。另外对于变速工况转折点是是否需要阶次跟踪。这套方案本身不依赖角域重采样因此即使转速有波动也可以直接用。但如果转速波动特别剧烈比如快速升降速建议先用转速计信号做阶次跟踪把时域信号映射到角域后再做CPWSES。这样的话特征频率变成阶次就不受转速影响效果更稳。不过这是另一个话题这里点到为止。4. 实操过程中的典型问题与排查实录这套方案我在多个数据集上试过包括公开数据如CWRU轴承数据、实验室实测数据和现场工业数据。效果整体很好但过程中也踩过不少坑列几个典型问题供参考。4.1 预白化把冲击信号也削掉了第一次跑通CPW时我遇到最大的问题是预白化之后的信号变得很平连轴承冲击都看不见了平方包络谱里也找不到特征频率。排查后发现问题出在倒谱截止值设得太小。我一开始设cut_queffrency 0.001也就是把周期小于1ms的调制成分删掉。但轴承故障冲击引起的周期调制通常在几个毫秒到几十毫秒量级这个设置把冲击对应的倒谱成分一起抹掉了。解决方法是把截止值调大到0.008~0.02s并对比不同截止值下的包络谱峰值。我还发现可以用一个简单规则先算出你关心的最小故障特征频率比如最低转速下的BPFO然后让截止倒频率大于 1/(2*f_min)这样冲击成分一定能保留。4.2 带通频带选错导致包络谱一塌糊涂有几次测试CPW和SES参数都没问题但包络谱里完全找不到特征频率。回来后重新检查发现问题在带通滤波频带选择上。直接用Kurtogram自动搜索时算法选中的最优频带可能是齿轮啮合谐波集中频带而不是轴承共振频带。这种频带里虽然峭度高但冲击性来自齿轮振动而非轴承故障。后来我总结了两条经验第一Kurtogram的结果要人工复核。把选择的频带信号画出来看时域波形确认它有典型的轴承冲击衰减波形幅值突然增大然后快速衰减而不是连续的齿轮调制波形。第二可以结合包络谱的物理意义来验证选频带后先对原始信号做包络谱看看有没有可疑的峰值结构尤其看BPFO附近的边带或谐波如果完全没有任何疑似轴承特征的结构说明频带很可能选错了。4.3 变速工况下包络谱特征频率漂移怎么办虽然CPW不依赖转速但如果你直接对一段转速变化明显的长信号做SES特征频率会因为转速漂移被抹宽成一个弧形的频率带峰值幅度下降看起来还是不明显。我实测的数据里有一段从600rpm加速到1800rpm的工况直接用CPWSES后特征频率峰几乎不可见。我的解决办法有两个都很务实一是分段处理。把30秒信号切成1~2秒的小段每段内部转速变化比较小SES谱线比较集中然后按时间顺序把每段的SES谱图拼成谱图瀑布图可以直观看到特征线的变化轨迹。二是结合短时包络谱/时频图。在matlab里可以用spectrogram对预白化信号的包络做短时傅里叶变换横轴时间、纵轴频率、颜色表示幅值特征频率随转速变化的轨迹一目了然。这个做出来发给现场工程师比任何数值指标都直观。segment_dur 1.0; % 每段时长1秒 seg_samples round(segment_dur * fs); num_seg floor(N / seg_samples); SES_stack zeros(num_seg, seg_samples/2); for i 1:num_seg seg x_w((i-1)*seg_samples1 : i*seg_samples); % 对该段做平方包络谱 % ... SES_stack(i, :) SES_seg; end % 绘制瀑布图 figure; imagesc((0:num_seg-1)*segment_dur, f_axis(1:seg_samples/2), log(SES_stack)); xlabel(时间 (s)); ylabel(频率 (Hz)); title(分段平方包络谱随时间变化);这个方法本质上是短时分段CPWSES相当于在时频平面上追踪特征频率轨迹。我没有用角域重采样那么重的操作但实测对转速范围不大比如1:3以内的场景足够用。4.4 一个容易被忽视的采样率与滤波器问题还有一个项目里踩到的坑带通滤波器在Matlab里设计时用的采样率如果和实际数据采样率不一致会导致滤波频带偏移整个分析结果全错。有个朋友拿我的代码去跑带通设了3k~8k结果他的数据采样率是25.6kHz我用代码里写的51.2kHz采样率直接滤波频带实际偏移到了6k~16k包络谱里特征频率完全找不到。所以每次跑数据前第一件事就是把采样率fs从原始数据里读出来不要写死在代码里。这个看似基础的问题在实际协作中出现的频率远比你想象的高。4.5 常见问题速查表现象可能原因解决方案预白化后信号幅值极小像被削平倒谱截止值过小冲击被误删增大 cut_queffrency调到0.01s以上包络谱中找不到任何轴承特征频率带通频带选择错误复核快速峭度图结果手动确认冲击波形特征频率峰值过度展宽转速变化明显分段做短时平方包络谱或在角域重采样谱图中出现在转频附近的密集边带齿轮调制未除干净减小 cut_queffrency把慢变调制成分清理更彻底低频段出现巨大直流分量包络平方后未去直流SES计算前减去均值两组数据结果差异大且不稳定信号长度不足频率分辨率低延长数据采集时间确保最少覆盖几十次冲击周期5. 变速工况下的扩展用法从单点诊断到趋势监测这套方法除了做单次数据诊断我在实际项目里还发现了一个特别好用的扩展——配合转速波动范围内的连续监测做故障趋势预判。现场的滚动轴承故障是一个渐进过程早期只是微小的剥落冲击能量很弱常规指标RMS、峭度几乎不敏感但CPWSES的平方包络谱中特征频率峰值已经可以稳定出现。我可以在现场采集一段数据跑完处理流程后把特征频率处的幅值作为健康指数按天记录形成趋势曲线。这个曲线对轴承早期退化非常敏感往往比振动总值提前几周甚至几个月发出预警。具体实施上有个小技巧如果现场有转速计可以把转速波动区间等分成几个小范围比如600~800rpm、800~1000rpm每个范围单独建立基线趋势。因为不同转速下轴承冲击能量本身就有差异不做转速归一化就把趋势混在一起容易误报。如果现场没有转速计也没关系CPWSES本身不依赖转速直接用固定频段的SES峰值做趋势虽然不如分段细致但也能看个大概。对于大多数现场来说大概准确的早期预警比精确但不实用的分析方案有价值得多。6. 一些个人体会和下一步想试的方向这套带通滤波倒谱预白化平方包络谱的方案从原理验证到现场应用我前后折腾了小半年。最初看文献接触CPW时总觉得它是个黑魔法——名字唬人公式抽象。真正把倒谱域高通滤波那段代码写出来、跑通数据后才理解它的本质就是一个自适应的频谱整形工具。我个人在实操中的体会是这套方法最值得称道的地方不是它有多先进而是它足够皮实不需要转速信号不需要复杂的参数寻优几个参数按照物理含义调整就能出稳定结果。把它的定位想清楚它就是包络分析的老路加了一个预处理模块但正是这个模块解决了变速工况下确定性成分掩盖轴承冲击的核心痛点。踩过几次坑之后总结一句话带通滤波负责聚焦倒谱预白化负责清理平方包络谱负责暴露。每一环都有自己明确的职责调试时哪一环效果不对就单独检查哪一环思路非常清晰。后面我还打算把这个流程整合进一个自动化分析脚本批量处理历史数据用平方包络谱特征峰值自动生成设备健康报告不再需要人工逐条看图。如果你也在做类似方向欢迎交流数据和处理经验。 
长河君 · 2026/9/11 14:38:03 · 阅读 1.2k 变速工况轴承故障诊断：倒谱预白化与平方包络谱组合方案
先聊一个很多做设备状态监测的朋友都绕不过去的场景现场振动数据一测一大把可一遇到变速工况传统频谱分析就抓瞎——转频在飘故障特征频率也跟着飘固定频点上一会儿有峰一会儿没峰整个谱图看着就像一团乱麻。我早年做风电齿轮箱和轧机轴承诊断时被这种时有时无的特征搞得头大过不少次。后来项目里逐步用上一套组合方案效果非常稳定带通滤波 倒谱预白化 平方包络谱英文缩写就是CPWCepstral Pre-whitening加SESSquared Envelope Spectrum。这套东西不挑转速不需要你额外装编码器也不用先做阶次跟踪直接对原始振动信号动手就能把淹没在齿轮啮合、结构共振、随机噪声里的轴承微弱冲击给捞出来。这篇文章就把我实际跑通的完整思路、原理、Matlab代码和踩过的坑一次讲清楚适合正在做轴承故障诊断研究、设备状态监测系统开发或者看论文看到CPW但不知道怎么落地的朋友。1. 项目背景变速工况下轴承故障检测为什么难1.1 轴承故障检测的基本逻辑滚动轴承故障诊断的底层逻辑很简单轴承内外圈、滚动体出现局部损伤剥落、点蚀、裂纹后滚动体经过损伤点时会产生周期性冲击。这种冲击会激起轴承座和传感器的结构共振形成一个高频衰减振荡。诊断这件事本质上就是想办法把这个微弱的周期冲击从振动信号里识别出来并算出它对应的特征频率。轴承故障特征频率有四个固定公式这里先摆出来后面代码要用外圈故障频率BPFO( BPFO \frac{n}{2} f_r \left(1 - \frac{d}{D} \cos \alpha\right) )内圈故障频率BPFI( BPFI \frac{n}{2} f_r \left(1 \frac{d}{D} \cos \alpha\right) )滚动体故障频率BSF( BSF \frac{D}{2d} f_r \left(1 - \left(\frac{d}{D} \cos \alpha\right)^2\right) )保持架故障频率FTF( FTF \frac{1}{2} f_r \left(1 - \frac{d}{D} \cos \alpha\right) )其中 ( n ) 是滚动体个数( d ) 是滚动体直径( D ) 是节圆直径( \alpha ) 是接触角( f_r ) 是转频。在恒定转速下这些频率是固定的直接用FFT幅值谱找峰值就行这是教科书里最基础的套路。但在现场转速总是波动的设备启停机、负载变化、调速运行转频一会儿30Hz一会儿50Hz故障特征频率成比例地跟着变。固定频率的谱峰被抹平在很宽的频带里甚至完全湮没在噪声中传统的包络谱分析效果很不理想。1.2 传统方法在变速工况下的失效原因变速工况下传统方法失效的原因可以拆成三层来看第一层转速变化导致特征频率非平稳。故障特征频率是转频的固定比例倍数转频一变特征频率就跟着变。做长时FFT时能量被分散到多个频点上峰值幅度被稀释信噪比下降。第二层振动信号里干扰成分复杂。齿轮啮合频率及其谐波、边带、结构共振、随机冲击混在一起。尤其齿轮箱里齿轮啮合能量通常比轴承冲击能量大一个数量级以上。变速时齿轮啮合频率也在变和轴承微弱信号搅在一起很难靠简单的频带选择分离。第三层经典的包络谱方法依赖共振频带选择。传统做法是先用带通滤波器把某个共振频带选出来再做包络谱。可是现场没有先验知识不知道轴承冲击激起了哪个共振峰。选宽了噪声多选窄了漏掉冲击成分。在变速工况下共振频带虽然相对固定但特征频率是漂移的如果滤波器中心频率和带宽没选好包络谱照样一塌糊涂。1.3 倒谱预白化在变速工况下的价值定位面对这些难题学术圈和工业界给出过很多解法阶次跟踪、时频分析、同步平均、盲解卷积MED、MCKD、谱峭度等。各有各的适用范围。**倒谱预白化CPW**属于近年来逐渐被重视的方法它在变速工况下有一个非常关键的优势——信号处理过程完全不依赖转速信息。CPW的思路是把振动信号在频域的对数幅值谱转换到倒谱域然后把倒谱中代表慢变结构分量比如齿轮啮合、共振调制、传输路径效应的低倒频率成分置零或削弱再还原到时域。这一顿操作之后信号里的确定性周期成分齿轮成分、轴频谐波被大幅压制残留的主要是非确定性成分和轴承冲击成分相当于对信号做了白化让后续的包络谱更容易凸显轴承故障特征。这种方法论上的含义是无论转速怎么变只要轴承有冲击CPW都能把冲击从确定性干扰中剥离出来。它不像阶次跟踪那样需要转速计也不像盲解卷积那样需要迭代优化和参数调优计算效率高、鲁棒性好特别适合工业现场数据。2. 核心原理拆解三条链路串起来看这套方案本质上是三个处理阶段的串联先用带通滤波把最能代表轴承冲击的共振频带抓出来再用倒谱预白化把残余的确定性干扰按在倒谱域里干掉最后用平方包络谱把冲击的周期性特征搬到频谱上找特征频率。下面逐个拆开讲原理和设计理由。2.1 倒谱预白化到底做了什么要理解CPW先得理解倒谱。倒谱的概念最早是Bogert等人1963年提出的用于回波检测。它的核心操作是对信号的功率谱取对数再对这个对数谱做一次逆傅里叶变换。得到的新域叫倒谱域横轴的单位叫倒频率quefrency量纲是时间。为什么要取对数再做变换因为对数可以把频谱中相乘的关系变成相加的关系。振动信号里的很多干扰比如传输路径效应、结构共振调制在频域表现为频谱包络的慢变起伏而轴承冲击则表现为频谱中均匀分布的精细谐波结构。这两者在原频谱中是相乘的但在对数谱中就变成了叠加。再做一次傅里叶变换它们在倒谱域中就分布在不同的倒频率区间低倒频率靠近0代表频谱的慢变包络对应结构共振、传输路径等确定性调制高倒频率代表频谱的快速波动对应周期性冲击的精细谐波结构。所以CPW的做法非常直观把倒谱中低倒频率的成分清零高通滤波然后把剩余的高倒频率成分还原成新的频谱幅值再结合原始相位做逆傅里叶变换得到预白化后的时域信号。这个过程相当于把频谱的形状抹平了把信号变成近似白噪声激励下的响应而白化后剩下的就是非确定性的、类冲击的成分。我个人的理解类比是倒谱域里做高通滤波就像照片处理里减去低频背景光只保留高频纹理细节。轴承冲击就是纹理细节齿轮啮合和共振就是背景光。2.2 带通滤波的作用与频带选择在CPW之前加一个带通滤波工程上的意义非常大。第一减少无关频率成分提高白化的针对性。原始振动信号频带往往很宽可能从几赫兹到几十千赫兹。如果直接对全频段做CPW那些远离轴承共振频带的噪声会在白化过程中被放大干扰后续分析。先带通滤波把注意力集中到轴承冲击占主导的频带上CPW才能更有效地凸出冲击。第二降低信号的非平稳性干扰。低频段的轴频、齿轮啮合频率及其谐波分布非常密集而且随转速变化剧烈。带通滤波可以把这些低频大能量成分挡在门外让后续处理聚焦于高频冲击成分。关键问题是频带怎么选。工程上有几种思路如果目标是轴承局部损伤用谱峭度Spectral Kurtosis或快速峭度图Fast Kurtogram自动寻找使峭度最大的频带中心频率和带宽。峭度大的频带意味着冲击性最强这是最常用的方法。如果项目里已知结构共振频带比如做过锤击试验或历史数据的谱分析直接采用先验频带。如果轴转速稳定或变化不大可以用边带能量最大化的方式选频带但在大范围变速时不太靠谱。我实测下来的经验是对于滚动轴承共振频带通常集中在2kHz到20kHz之间带宽一般选原始采样率能覆盖的1/4到1/8左右。不过最稳的还是先用Fast Kurtogram扫一遍把Kurtogram峰值对应的频带拿来做带通滤波。在Matlab里可以直接用kurtogram函数第三方工具箱比如Antoni的Fast Kurtogram或者自己写短时傅里叶变换算峭度。2.3 平方包络谱为什么能看见故障特征预白化之后信号已经干净了很多但仍然是时域信号冲击的周期性还没有直接呈现。**平方包络谱SES**就是把这层周期性挖出来的关键一步。包络谱的基本原理是通过Hilbert变换构造解析信号取模得到包络信号。包络信号的物理意义是原始信号的瞬时幅值变化轴承故障冲击引起的幅值调制正好体现在包络中。平方向操作的原因有两个一个是数学上的——经典窄带包络谱也就是对包络直接做FFT得到的幅值谱在某些情况下会出现谱线幅度和特征频率幅值不成线性关系的问题。而平方包络谱在数学推导中和谱相关分析有密切联系理论上对非平稳冲击信号的检测性能更好也更容易在平方包络谱中看到故障特征频率的二次谐波、三次谐波等。另一个是工程上的——对包络做平方运算等价于加强了冲击成分的幅值调制深度让周期冲击的能量更集中到特征频率及其谐波上信噪比更高。我做仿真和实测数据对比时同样的数据平方包络谱里特征频率峰值的突出程度明显优于直接包络谱。所以这三个环节的组合逻辑是带通滤波负责划重点倒谱预白化负责打掩护的抹掉平方包络谱负责亮出真身。三条链路环环相扣少了哪一环效果都会打折扣。3. Matlab实现流程与关键代码理论说清楚之后上干货。这一节给出我实际在Matlab里跑通的完整流程包括参数如何选取、每个环节为什么要这样写、哪些地方容易踩坑。3.1 整体流程设计完整处理流程可以分成六步读取加速度振动信号设定采样率 fs对信号做快速峭度图分析或根据先验知识确定带通滤波频带 [f_lo, f_hi]设计带通滤波器对原始信号滤波得到带通信号 x_bp对 x_bp 做倒谱预白化得到预白化信号 x_w对 x_w 做Hilbert变换求包络取平方做FFT得到平方包络谱在平方包络谱中定位BPFO、BPFI、BSF、FTF及其谐波判断故障类型。流程图用文字描述就是这样原始信号 → [带通滤波] → [倒谱预白化] → [平方包络谱] → 特征频率识别。字不多条理很清楚。3.2 关键函数与Matlab代码实现下面按环节给出核心代码。由于实际项目的数据差异很大这里用示意代码配合注释关键是体现实现思路。第一步带通滤波% 参数定义 fs 25600; % 采样率单位Hz现场采集卡一般2.56k~100k N length(x); % 信号长度 % 带通频带选择——这里以快速峭度图提取结果为例 % 假设kurtogram计算得到最佳频带为 [fc - Bw/2, fc Bw/2] fc 6000; % 中心频率单位Hz Bw 4000; % 带宽单位Hz f_lo fc - Bw/2; f_hi fc Bw/2; % 设计带通滤波器 bpFilt designfilt(bandpassiir, ... FilterOrder, 8, ... HalfPowerFrequency1, f_lo, ... HalfPowerFrequency2, f_hi, ... SampleRate, fs, ... DesignMethod, butter); % 零相位滤波避免相位畸变影响包络分析 x_bp filtfilt(bpFilt, x);这里我特意用filtfilt替代filter因为零相位滤波不会让包络产生时延畸变。对于后续的包络分析相位一致性很重要这一点后面在避坑部分会再强调。第二步倒谱预白化function x_w cepstral_prewhitening(x, cut_queffrency) % CEPSTRAL_PREWHITENING 倒谱预白化 % 输入 % x - 输入时域信号建议先做过带通滤波 % cut_queffrency - 倒频率截止值单位s小于该值的成分被置零 % 输出 % x_w - 预白化后的时域信号 N length(x); % 加汉宁窗减少频谱泄漏可选但推荐 win hanning(N, periodic); xw x(:) .* win; % 计算FFT X fft(xw); X_amp abs(X); X_phase angle(X); % 避免幅值为0导致对数无穷大加一个极小值 eps_val 1e-12; log_amp log(X_amp eps_val); % 对对数幅值做IFFT得到倒谱 cep real(ifft(log_amp)); % 倒谱域高通将小于截止倒频率的成分置零 % 注意倒谱的横轴是时间0对应直流N-1对应(N-1)/fs % 截止倒频率对应样本数需要换算 cut_samples round(cut_queffrency * fs); cep_cut cep; cep_cut(1:cut_samples) 0; cep_cut(end-cut_samples2:end) 0; % 对称部分也置零 % 由处理后的倒谱重建对数幅值谱 log_amp_white real(fft(cep_cut)); % 重建白化后的幅值谱归一化处理 amp_white exp(log_amp_white); % 白化过程将原频谱幅值除以重建的慢变幅值保留相位 % 这样等同于把频谱“抹平” X_white (X_amp ./ (amp_white eps_val)) .* exp(1j * X_phase); % 逆变换回时域取实部 x_w real(ifft(X_white)); % 由于之前加了窗这里做简单幅度修正窗函数能量补偿 win_energy sum(win.^2) / N; x_w x_w / win_energy; % 去除直流分量 x_w x_w - mean(x_w); end这里需要解释几个关键细节。第一倒谱域的高通操作截止值cut_queffrency的单位是秒需要乘以采样率换算成样本数。对于旋转机械通常设置截止倒频率为 1/(2*最大特征频率) 或根据齿轮啮合周期确定。我平时的起步值是0.01秒对应100Hz的周期然后看结果调整。第二因为倒谱是对称的置零时要同时处理正半轴和负半轴对应位置否则重建出来的频谱会不对称。第三重建时用原幅值除以估计的慢变幅值就实现了白化的效果——原来幅值高的地方被压低原来幅值低的地方被抬升频谱趋于平坦。第三步平方包络谱% 对预白化信号求包络 analytic_signal hilbert(x_w); envelope abs(analytic_signal); % 平方包络 envelope_sq envelope.^2; % 去除直流后计算FFT N length(envelope_sq); win_e hanning(N, periodic); env_win (envelope_sq - mean(envelope_sq)) .* win_e; SES abs(fft(env_win)) * 2 / N; % 单边谱 % 频率轴 f_axis (0:N/2-1) * fs / N; % 绘制平方包络谱 figure; plot(f_axis, SES(1:N/2)); xlabel(频率 (Hz)); ylabel(幅值); title(平方包络谱 SES); xlim([0, min(fs/2, 1000)]); % 关注低频段高频段一般是噪声这里的几个细节要提醒包络平方后直流分量很大作图前先去直流不然低频段会被巨大的直流旁瓣淹没加窗是为了减少谱泄漏但窗口会引入旁瓣对幅度精度要求高的话也可以不加窗看你的场景。3.3 参数怎么定实测经验参数选择是这套方法落地最容易出问题的地方。我根据自己的实际测试总结一张参数参考表参数建议范围/取值说明带通滤波中心频率 fc2kHz ~ 15kHz具体用快速峭度图或共振分析确定不要盲目定带通滤波带宽 Bw1kHz ~ 8kHz太窄容易丢失冲击成分太宽噪声多带通滤波器阶数6~10阶 Butterworth阶数太高相位畸变严重太低过渡带太宽倒谱截止倒频率 cut_queffrency0.005s ~ 0.02s需大于最大分析频率周期的2倍我一般先取0.01s再调信号时长至少包含50~100次故障冲击例如BPFO80Hz采集时长至少2~3s才能有足够分辨率采样率 fs至少是最高关心频率的8~10倍处理20kHz共振频带建议不低于50kHz采样倒谱截止倒频率的物理含义要强调一下它决定了你抹掉多慢的调制成分。设得太小很多有用的冲击调制信息也被抹掉了设得太大齿轮啮合产生的边带成分清理不干净。实际调试时我会把cut_queffrency从0.002s到0.05s扫一遍看哪种参数下平方包络谱中特征频率峰值最突出然后定下来。另外对于变速工况转折点是是否需要阶次跟踪。这套方案本身不依赖角域重采样因此即使转速有波动也可以直接用。但如果转速波动特别剧烈比如快速升降速建议先用转速计信号做阶次跟踪把时域信号映射到角域后再做CPWSES。这样的话特征频率变成阶次就不受转速影响效果更稳。不过这是另一个话题这里点到为止。4. 实操过程中的典型问题与排查实录这套方案我在多个数据集上试过包括公开数据如CWRU轴承数据、实验室实测数据和现场工业数据。效果整体很好但过程中也踩过不少坑列几个典型问题供参考。4.1 预白化把冲击信号也削掉了第一次跑通CPW时我遇到最大的问题是预白化之后的信号变得很平连轴承冲击都看不见了平方包络谱里也找不到特征频率。排查后发现问题出在倒谱截止值设得太小。我一开始设cut_queffrency 0.001也就是把周期小于1ms的调制成分删掉。但轴承故障冲击引起的周期调制通常在几个毫秒到几十毫秒量级这个设置把冲击对应的倒谱成分一起抹掉了。解决方法是把截止值调大到0.008~0.02s并对比不同截止值下的包络谱峰值。我还发现可以用一个简单规则先算出你关心的最小故障特征频率比如最低转速下的BPFO然后让截止倒频率大于 1/(2*f_min)这样冲击成分一定能保留。4.2 带通频带选错导致包络谱一塌糊涂有几次测试CPW和SES参数都没问题但包络谱里完全找不到特征频率。回来后重新检查发现问题在带通滤波频带选择上。直接用Kurtogram自动搜索时算法选中的最优频带可能是齿轮啮合谐波集中频带而不是轴承共振频带。这种频带里虽然峭度高但冲击性来自齿轮振动而非轴承故障。后来我总结了两条经验第一Kurtogram的结果要人工复核。把选择的频带信号画出来看时域波形确认它有典型的轴承冲击衰减波形幅值突然增大然后快速衰减而不是连续的齿轮调制波形。第二可以结合包络谱的物理意义来验证选频带后先对原始信号做包络谱看看有没有可疑的峰值结构尤其看BPFO附近的边带或谐波如果完全没有任何疑似轴承特征的结构说明频带很可能选错了。4.3 变速工况下包络谱特征频率漂移怎么办虽然CPW不依赖转速但如果你直接对一段转速变化明显的长信号做SES特征频率会因为转速漂移被抹宽成一个弧形的频率带峰值幅度下降看起来还是不明显。我实测的数据里有一段从600rpm加速到1800rpm的工况直接用CPWSES后特征频率峰几乎不可见。我的解决办法有两个都很务实一是分段处理。把30秒信号切成1~2秒的小段每段内部转速变化比较小SES谱线比较集中然后按时间顺序把每段的SES谱图拼成谱图瀑布图可以直观看到特征线的变化轨迹。二是结合短时包络谱/时频图。在matlab里可以用spectrogram对预白化信号的包络做短时傅里叶变换横轴时间、纵轴频率、颜色表示幅值特征频率随转速变化的轨迹一目了然。这个做出来发给现场工程师比任何数值指标都直观。segment_dur 1.0; % 每段时长1秒 seg_samples round(segment_dur * fs); num_seg floor(N / seg_samples); SES_stack zeros(num_seg, seg_samples/2); for i 1:num_seg seg x_w((i-1)*seg_samples1 : i*seg_samples); % 对该段做平方包络谱 % ... SES_stack(i, :) SES_seg; end % 绘制瀑布图 figure; imagesc((0:num_seg-1)*segment_dur, f_axis(1:seg_samples/2), log(SES_stack)); xlabel(时间 (s)); ylabel(频率 (Hz)); title(分段平方包络谱随时间变化);这个方法本质上是短时分段CPWSES相当于在时频平面上追踪特征频率轨迹。我没有用角域重采样那么重的操作但实测对转速范围不大比如1:3以内的场景足够用。4.4 一个容易被忽视的采样率与滤波器问题还有一个项目里踩到的坑带通滤波器在Matlab里设计时用的采样率如果和实际数据采样率不一致会导致滤波频带偏移整个分析结果全错。有个朋友拿我的代码去跑带通设了3k~8k结果他的数据采样率是25.6kHz我用代码里写的51.2kHz采样率直接滤波频带实际偏移到了6k~16k包络谱里特征频率完全找不到。所以每次跑数据前第一件事就是把采样率fs从原始数据里读出来不要写死在代码里。这个看似基础的问题在实际协作中出现的频率远比你想象的高。4.5 常见问题速查表现象可能原因解决方案预白化后信号幅值极小像被削平倒谱截止值过小冲击被误删增大 cut_queffrency调到0.01s以上包络谱中找不到任何轴承特征频率带通频带选择错误复核快速峭度图结果手动确认冲击波形特征频率峰值过度展宽转速变化明显分段做短时平方包络谱或在角域重采样谱图中出现在转频附近的密集边带齿轮调制未除干净减小 cut_queffrency把慢变调制成分清理更彻底低频段出现巨大直流分量包络平方后未去直流SES计算前减去均值两组数据结果差异大且不稳定信号长度不足频率分辨率低延长数据采集时间确保最少覆盖几十次冲击周期5. 变速工况下的扩展用法从单点诊断到趋势监测这套方法除了做单次数据诊断我在实际项目里还发现了一个特别好用的扩展——配合转速波动范围内的连续监测做故障趋势预判。现场的滚动轴承故障是一个渐进过程早期只是微小的剥落冲击能量很弱常规指标RMS、峭度几乎不敏感但CPWSES的平方包络谱中特征频率峰值已经可以稳定出现。我可以在现场采集一段数据跑完处理流程后把特征频率处的幅值作为健康指数按天记录形成趋势曲线。这个曲线对轴承早期退化非常敏感往往比振动总值提前几周甚至几个月发出预警。具体实施上有个小技巧如果现场有转速计可以把转速波动区间等分成几个小范围比如600~800rpm、800~1000rpm每个范围单独建立基线趋势。因为不同转速下轴承冲击能量本身就有差异不做转速归一化就把趋势混在一起容易误报。如果现场没有转速计也没关系CPWSES本身不依赖转速直接用固定频段的SES峰值做趋势虽然不如分段细致但也能看个大概。对于大多数现场来说大概准确的早期预警比精确但不实用的分析方案有价值得多。6. 一些个人体会和下一步想试的方向这套带通滤波倒谱预白化平方包络谱的方案从原理验证到现场应用我前后折腾了小半年。最初看文献接触CPW时总觉得它是个黑魔法——名字唬人公式抽象。真正把倒谱域高通滤波那段代码写出来、跑通数据后才理解它的本质就是一个自适应的频谱整形工具。我个人在实操中的体会是这套方法最值得称道的地方不是它有多先进而是它足够皮实不需要转速信号不需要复杂的参数寻优几个参数按照物理含义调整就能出稳定结果。把它的定位想清楚它就是包络分析的老路加了一个预处理模块但正是这个模块解决了变速工况下确定性成分掩盖轴承冲击的核心痛点。踩过几次坑之后总结一句话带通滤波负责聚焦倒谱预白化负责清理平方包络谱负责暴露。每一环都有自己明确的职责调试时哪一环效果不对就单独检查哪一环思路非常清晰。后面我还打算把这个流程整合进一个自动化分析脚本批量处理历史数据用平方包络谱特征峰值自动生成设备健康报告不再需要人工逐条看图。如果你也在做类似方向欢迎交流数据和处理经验。