变速工况轴承故障诊断：倒谱预白化与平方包络谱组合方案

先聊一个很多做设备状态监测的朋友都绕不过去的场景：现场振动数据一测一大把，可一遇到变速工况，传统频谱分析就抓瞎——转频在飘，故障特征频率也跟着飘，固定频点上一会儿有峰一会儿没峰，整个谱图看着就像一团乱麻。我早年做风电齿轮箱和轧机轴承诊断时，被这种"时有时无"的特征搞得头大过不少次。

后来项目里逐步用上一套组合方案，效果非常稳定：带通滤波 + 倒谱预白化 + 平方包络谱，英文缩写就是CPW（Cepstral Pre-whitening）加SES（Squared 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 整体流程设计

完整处理流程可以分成六步：

1. 读取加速度振动信号，设定采样率 fs；
2. 对信号做快速峭度图分析或根据先验知识确定带通滤波频带 [f_lo, f_hi]；
3. 设计带通滤波器，对原始信号滤波，得到带通信号 x_bp；
4. 对 x_bp 做倒谱预白化，得到预白化信号 x_w；
5. 对 x_w 做Hilbert变换求包络，取平方，做FFT得到平方包络谱；
6. 在平方包络谱中定位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_samples+2: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 参数怎么定（实测经验）

参数选择是这套方法落地最容易出问题的地方。我根据自己的实际测试，总结一张参数参考表：

倒谱截止倒频率的物理含义要强调一下：它决定了你"抹掉"多慢的调制成分。设得太小，很多有用的冲击调制信息也被抹掉了；设得太大，齿轮啮合产生的边带成分清理不干净。实际调试时，我会把cut_queffrency从0.002s到0.05s扫一遍，看哪种参数下平方包络谱中特征频率峰值最突出，然后定下来。

另外，对于变速工况，转折点是"是否需要阶次跟踪"。这套方案本身不依赖角域重采样，因此即使转速有波动也可以直接用。但如果转速波动特别剧烈（比如快速升降速），建议先用转速计信号做阶次跟踪，把时域信号映射到角域后再做CPW+SES。这样的话，特征频率变成"阶次"，就不受转速影响，效果更稳。不过这是另一个话题，这里点到为止。

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的工况，直接用CPW+SES后特征频率峰几乎不可见。我的解决办法有两个，都很务实：

一是分段处理。把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_samples+1 : 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('分段平方包络谱随时间变化');

这个方法本质上是"短时分段+CPW+SES"，相当于在时频平面上追踪特征频率轨迹。我没有用角域重采样那么重的操作，但实测对转速范围不大（比如1:3以内）的场景足够用。

4.4 一个容易被忽视的采样率与滤波器问题

还有一个项目里踩到的坑带通滤波器在Matlab里设计时用的采样率如果和实际数据采样率不一致，会导致滤波频带偏移，整个分析结果全错。有个朋友拿我的代码去跑，带通设了3k~8k，结果他的数据采样率是25.6kHz，我用代码里写的51.2kHz采样率直接滤波，频带实际偏移到了6k~16k，包络谱里特征频率完全找不到。

所以每次跑数据前，第一件事就是把采样率fs从原始数据里读出来，不要写死在代码里。这个看似基础的问题在实际协作中出现的频率远比你想象的高。

4.5 常见问题速查表

5. 变速工况下的扩展用法：从单点诊断到趋势监测

这套方法除了做单次数据诊断，我在实际项目里还发现了一个特别好用的扩展——配合转速波动范围内的连续监测，做故障趋势预判。

现场的滚动轴承故障是一个渐进过程：早期只是微小的剥落，冲击能量很弱，常规指标（RMS、峭度）几乎不敏感；但CPW+SES的平方包络谱中，特征频率峰值已经可以稳定出现。我可以在现场采集一段数据，跑完处理流程后，把特征频率处的幅值作为"健康指数"，按天记录，形成趋势曲线。这个曲线对轴承早期退化非常敏感，往往比振动总值提前几周甚至几个月发出预警。

具体实施上有个小技巧：如果现场有转速计，可以把转速波动区间等分成几个小范围（比如600~800rpm、800~1000rpm），每个范围单独建立基线趋势。因为不同转速下轴承冲击能量本身就有差异，不做转速归一化就把趋势混在一起，容易误报。

如果现场没有转速计，也没关系，CPW+SES本身不依赖转速，直接用固定频段的SES峰值做趋势，虽然不如分段细致，但也能看个大概。对于大多数现场来说，"大概准确"的早期预警比"精确但不实用"的分析方案有价值得多。

6. 一些个人体会和下一步想试的方向

这套"带通滤波+倒谱预白化+平方包络谱"的方案，从原理验证到现场应用，我前后折腾了小半年。最初看文献接触CPW时，总觉得它是个"黑魔法"——名字唬人，公式抽象。真正把倒谱域高通滤波那段代码写出来、跑通数据后，才理解它的本质就是一个自适应的频谱整形工具。

我个人在实操中的体会是，这套方法最值得称道的地方不是它有多"先进"，而是它足够"皮实"：不需要转速信号，不需要复杂的参数寻优，几个参数按照物理含义调整就能出稳定结果。把它的定位想清楚，它就是包络分析的老路加了一个"预处理模块"，但正是这个模块，解决了变速工况下确定性成分掩盖轴承冲击的核心痛点。

踩过几次坑之后总结一句话：带通滤波负责聚焦，倒谱预白化负责清理，平方包络谱负责暴露。每一环都有自己明确的职责，调试时哪一环效果不对就单独检查哪一环，思路非常清晰。

后面我还打算把这个流程整合进一个自动化分析脚本，批量处理历史数据，用平方包络谱特征峰值自动生成设备健康报告，不再需要人工逐条看图。如果你也在做类似方向，欢迎交流数据和处理经验。