CWRU轴承故障诊断：STFT与CWT时频分析实战指南

简介这份资源围绕凯斯西储大学CWRU轴承故障数据集展开信号时频分析面向从事机械设备故障诊断、健康监测与振动信号处理的研究生、工程师及科研人员帮助读者理解如何从实测振动数据中提取故障特征。压缩包内仅含1个docx文档约1.05MB以图文与代码结合的方式组织内容。文档先介绍CWRU数据集的实验台构成包括1.5KW电机、扭矩传感器与功率测试计并区分驱动端、风扇端和基座三处加速度计信号所反映的物理特征随后以正常信号与0.021英寸内圈、滚珠、外圈故障信号为对象演示短时傅里叶变换在重叠比例0.5、尺度16/32/64下的时频图差异说明尺度32的选取依据再对比morl、cmor1-1、cmor1.5-2、cgau8等小波函数最终确定cmor1.5-2并比较尺度32至256的连续小波变换效果。目前已有1582人学习适合作为时频分析入门与故障分类特征工程的实操参考。1. 凯斯西储大学轴承故障数据集为什么时频分析是绕不开的第一道坎如果你手头有一份凯斯西储大学CWRU轴承故障数据打算做故障诊断大概率会先画个时域波形看看。然后你会发现正常轴承和故障轴承的时域波形肉眼看上去差别并不大——尤其是内圈故障和滚动体故障波形毛刺混在一起光看幅值根本分不清谁是谁。这不是你数据读错了而是时域信号本身把频率信息藏起来了。凯斯西储大学轴承故障数据集之所以成为旋转机械故障诊断领域被反复使用的基准核心原因是它的实验条件标注清晰故障类型内圈、外圈、滚动体、故障直径0.007、0.014、0.021 英寸、负载03 马力都有明确记录。但标注清晰不等于信号好分。轴承故障的冲击成分是瞬态的、非平稳的傅里叶变换把整个时段的频率成分平均掉之后冲击发生的时刻信息就丢了。时频分析要解决的正是这个问题把一维时间信号映射到时间-频率二维平面让冲击在哪个时刻、哪个频带出现变得可见。这篇内容面向的是已经拿到 CWRU 数据、准备做特征提取或故障分类的工程师和研究生。我会从数据加载开始把短时傅里叶变换STFT和连续小波变换CWT两条最常用的时频分析路径走通给出可复现的代码、参数设置逻辑以及我在实际处理 CWRU 数据时踩过的坑。不涉及深度学习模型搭建只聚焦时频分析这一步怎么做扎实。2. 先把 CWRU 数据读对文件结构、采样率与通道选择2.1 CWRU 数据的目录结构和命名规则CWRU 数据集在官网以 .mat 文件分发命名规则通常是「故障类型 故障直径 负载」的组合。比如IR007_0.mat表示内圈故障、直径 0.007 英寸、负载 0 马力B014_1.mat表示滚动体故障、直径 0.014 英寸、负载 1 马力正常基线数据一般命名为Normal_0.mat到Normal_3.mat。驱动端加速度计数据存在DE_time变量里风扇端存在FE_time基座存在BA_time。做时频分析我一般只用驱动端DE_time。原因很直接驱动端离轴承座最近故障冲击的信噪比最高风扇端信号经过传递路径衰减后冲击成分容易被淹没。如果你要做多通道融合那是另一个话题单通道时频分析从 DE 通道起步最稳。import scipy.io as sio import numpy as np # 读取 CWRU 驱动端数据 def load_cwru_de(filepath): 读取 CWRU .mat 文件中的驱动端加速度信号 filepath: .mat 文件路径 返回: 一维 numpy 数组 (DE_time) mat_data sio.loadmat(filepath) # 变量名通常是 DE_time但不同版本可能带后缀 de_key [k for k in mat_data.keys() if DE_time in k] if not de_key: raise KeyError(f未找到 DE_time 变量现有键: {list(mat_data.keys())}) signal mat_data[de_key[0]].flatten() return signal signal load_cwru_de(IR007_0.mat) print(f信号长度: {len(signal)}, 采样点数: {signal.shape})这段代码的关键在变量名匹配。CWRU 官网下载的 .mat 文件里变量名有时是X105_DE_time这种带编号的形式直接写mat_data[DE_time]会报 KeyError。用列表推导式模糊匹配DE_time是更稳妥的做法。flatten()把列向量压成一维后续做时频变换时不用再关心原始维度。2.2 采样率确认与重采样决策CWRU 驱动端数据的采样率是 12 kHz12,000 Hz部分文件是 48 kHz。这个信息不在 .mat 文件里需要根据你下载的数据版本确认。采样率搞错的后果很严重STFT 的频率轴会整体偏移CWT 的尺度-频率对应关系也会错最后得到的时频图看着有模有样但频率刻度全是错的。我一般会在代码里显式定义采样率而不是从数据里猜FS 12000 # CWRU 驱动端常用采样率48kHz 数据需改为 48000 # 如果要做降采样先做抗混叠滤波 from scipy.signal import decimate def downsample_signal(signal, fs_original, fs_target): 降采样带抗混叠滤波 decimate 的 q 参数必须是整数 q fs_original // fs_target if q 2: return signal, fs_original signal_ds decimate(signal, q, ftypeiir, zero_phaseTrue) return signal_ds, fs_original // q signal_ds, fs_new downsample_signal(signal, FS, 3000) print(f降采样后采样率: {fs_new}, 点数: {len(signal_ds)})降采样不是必须的但如果你后续要做 CWT 且尺度范围设得比较大原始 12 kHz 数据计算量会很大。降到 3 kHz 通常足够覆盖轴承故障的特征频率——CWRU 轴承故障特征频率一般在几百赫兹到 2 kHz 之间3 kHz 采样率满足奈奎斯特条件。decimate函数自带 IIR 抗混叠滤波比直接切片signal[::q]安全得多后者会把高频噪声混叠到低频时频图上会出现虚假的低频能量带。注意降采样前一定确认目标采样率大于你关心的最高故障特征频率的两倍。CWRU 内圈故障特征频率在负载 0 时约 160 Hz 左右但冲击的宽带成分可以延伸到 5 kHz 以上做包络分析时降采样要谨慎。3. 短时傅里叶变换窗长、重叠与频率分辨率的三角权衡3.1 STFT 在 CWRU 信号上的参数选择逻辑STFT 的核心思想是用一个滑动窗把长信号切成短段对每段做 FFT再把结果按时间排列。窗长决定了时间分辨率和频率分辨率的平衡窗越长频率分辨率越高但时间定位越模糊窗越短时间定位越准但频率分辨率下降。这不是玄学是海森堡不确定性原理在信号处理里的直接体现。对 CWRU 轴承故障信号我一般从窗长 256 点起步12 kHz 采样率下约 21 ms重叠 50%75%。故障冲击的持续时间通常在几毫秒量级256 点窗能覆盖一个完整的冲击衰减过程同时频率分辨率约 47 Hz12000/256足以区分故障特征频率及其谐波。from scipy.signal import stft import matplotlib.pyplot as plt def compute_stft(signal, fs, nperseg256, noverlapNone): 计算 STFT nperseg: 窗长采样点数 noverlap: 重叠点数默认 75% 重叠 if noverlap is None: noverlap int(nperseg * 0.75) f, t, Zxx stft(signal, fsfs, npersegnperseg, noverlapnoverlap, windowhann, detrendFalse, return_onesidedTrue) return f, t, Zxx f, t, Zxx compute_stft(signal, FS, nperseg256) magnitude np.abs(Zxx) plt.figure(figsize(10, 4)) plt.pcolormesh(t, f, 20*np.log10(magnitude 1e-10), shadinggouraud) plt.ylabel(Frequency (Hz)) plt.xlabel(Time (s)) plt.title(STFT Magnitude (dB)) plt.colorbar(labeldB) plt.ylim(0, 3000) plt.tight_layout() plt.show()windowhann是默认选择汉宁窗的旁瓣衰减比矩形窗好能减少频谱泄漏。detrendFalse表示不做去趋势因为 CWRU 信号已经做过预处理再去趋势可能把故障冲击的低频成分也去掉。return_onesidedTrue对实数信号只返回正频率半边省一半计算量。3.2 从 STFT 图里读出故障特征频率轴和实际特征频率的对应STFT 图出来之后关键是从图上找到故障特征频率的踪迹。CWRU 内圈故障特征频率BPFI的计算公式是BPFI (n/2) × fr × (1 (d/D) × cosθ)其中 n 是滚动体个数fr 是轴转频d 是滚动体直径D 是节径θ 是接触角。CWRU 实验台的具体参数n9d0.3126 英寸D1.537 英寸θ0。负载 0 时轴转频约 29.95 Hz算出来 BPFI 约 162 Hz。在 STFT 图上你应该能在 162 Hz 及其倍频324 Hz、486 Hz…附近看到能量集中。如果看不到先检查采样率是否设对再检查窗长是否太长导致频率分辨率不够——窗长 256 点时频率分辨率 47 Hz162 Hz 和 324 Hz 能分开但如果窗长只有 64 点分辨率降到 187 Hz这两个峰就糊在一起了。# 提取 STFT 在特定频率切片上的时间演化 def extract_freq_slice(f, t, Zxx, target_freq, bandwidth20): 提取目标频率附近的时间-幅值曲线 target_freq: 目标频率 (Hz) bandwidth: 频率带宽 (Hz) freq_mask (f target_freq - bandwidth) (f target_freq bandwidth) slice_mag np.mean(np.abs(Zxx[freq_mask, :]), axis0) return t, slice_mag t_slice, mag_slice extract_freq_slice(f, t, Zxx, target_freq162, bandwidth20) # mag_slice 的周期性峰值对应冲击重复频率这段代码把 STFT 矩阵在频率轴上做切片得到目标频率附近能量随时间的变化。如果轴承有内圈故障mag_slice会呈现周期性的峰值峰峰间隔对应冲击重复周期。这个周期信息在时域波形上很难直接量出来但在时频平面上变得直观。3.3 STFT 的局限窗长固定带来的分辨率困境STFT 最大的问题是一次变换只能用一种窗长。轴承故障信号里冲击成分是高频瞬态需要短窗来定位而故障特征频率的谐波结构是低频持续成分需要长窗来分辨。固定窗长意味着你必须在两者之间做取舍。我的经验是如果目标是检测冲击是否存在窗长取 128256 点如果目标是精确测量特征频率窗长取 5121024 点。不要试图用一个窗长解决所有问题那只会两头不讨好。这也是为什么后面要引入 CWT——小波变换的多分辨率特性天然适合处理这种「高频短时、低频长时」的信号。4. 连续小波变换尺度参数怎么选、复小波为什么更适合轴承信号4.1 CWT 的尺度-频率映射关系与参数初始化CWT 用一组尺度可变的母小波去卷信号尺度小对应高频、时间分辨率高尺度大对应低频、频率分辨率高。这个特性正好匹配轴承故障信号的结构冲击瞬间用小尺度捕捉谐波结构用大尺度分辨。在 Python 里做 CWT常用PyWavelets库。母小波选复 Morlet 小波cmor因为轴承故障信号是调制信号复小波的实部和虚部能同时给出幅值和相位信息后续做包络分析更方便。import pywt def compute_cwt(signal, fs, scalesNone, waveletcmor1.5-1.0): 计算连续小波变换 scales: 尺度数组None 时自动生成 wavelet: 母小波名称cmorB-C 中 B 是带宽参数C 是中心频率 if scales is None: # 尺度范围对应频率范围f (center_freq * fs) / scale # cmor1.5-1.0 的中心频率约 1.0 Hz归一化 scales np.arange(1, 128) coefficients, frequencies pywt.cwt(signal, scales, wavelet, sampling_period1.0/fs) return coefficients, frequencies coeffs, freqs compute_cwt(signal, FS, scalesnp.arange(1, 256)) print(fCWT 系数矩阵形状: {coeffs.shape}) print(f频率范围: {freqs[0]:.1f} Hz ~ {freqs[-1]:.1f} Hz)cmor1.5-1.0里的1.5是带宽参数1.0是中心频率。带宽参数越大小波在频率域的带宽越宽时间分辨率越低。对轴承冲击信号我一般用cmor1.5-1.0或cmor1.0-1.0前者频率聚集性稍好后者时间定位更锐利。scales从 1 到 256 覆盖的频率范围取决于采样率12 kHz 采样率下scale1 对应约 12 kHzscale256 对应约 47 Hz。如果你只关心 03 kHz把 scales 上限设到 256 就够了再大只会增加计算量。4.2 用 CWT 系数矩阵做故障可视化的完整流程CWT 系数矩阵是复数矩阵取模之后得到幅值矩阵画出来就是时频图。和 STFT 图相比CWT 图在低频段频率分辨率更高在高频段时间分辨率更高冲击的起始和结束时刻更清晰。import matplotlib.pyplot as plt def plot_cwt_scalogram(coeffs, freqs, fs, titleCWT Scalogram): 绘制 CWT 尺度图 coeffs: CWT 系数矩阵 (len(scales), len(signal)) freqs: 对应的频率数组 magnitude np.abs(coeffs) # 转成 dB 尺度动态范围设 60 dB mag_db 20 * np.log10(magnitude 1e-10) vmax np.max(mag_db) vmin vmax - 60 plt.figure(figsize(10, 5)) plt.pcolormesh(np.arange(len(coeffs[0]))/fs, freqs, mag_db, shadinggouraud, vminvmin, vmaxvmax, cmapjet) plt.ylabel(Frequency (Hz)) plt.xlabel(Time (s)) plt.yscale(log) plt.title(title) plt.colorbar(labelMagnitude (dB)) plt.ylim(10, fs/2) plt.tight_layout() plt.show() plot_cwt_scalogram(coeffs, freqs, FS, titleCWT Scalogram - IR007)vmin和vmax的设置很关键。如果不设matplotlib 会自动按数据最小最大值拉伸噪声会被放大成和冲击一样的亮度图上一片花。设 60 dB 动态范围是经验值冲击能量比背景噪声高 4060 dB这个范围能把冲击突出来同时压住噪声。plt.yscale(log)让低频段展开因为轴承故障特征频率集中在低频线性频率轴会把低频区域压得很扁。4.3 CWT 和 STFT 在 CWRU 数据上的效果对比同一段内圈故障信号STFT 和 CWT 给出的信息侧重不同。STFT 图在 162 Hz 附近能看到周期性的能量团但每个能量团的时间宽度被窗长模糊了CWT 图在同样频率位置能看到更锐利的冲击时刻而且冲击的衰减过程从高频到低频的能量下降更清楚。这不是说 CWT 一定比 STFT 好。STFT 的计算量小参数少适合做实时监测的初步筛查CWT 计算量大参数多母小波类型、带宽、中心频率、尺度范围但多分辨率特性让它在分析非平稳信号时信息量更大。我的做法是先用 STFT 快速扫一遍确认故障特征频率的大致位置再用 CWT 在目标频段做精细分析。提示CWT 的尺度数组不要用默认值。pywt.cwt如果不传 scales会用一个默认的尺度范围那个范围不一定覆盖你关心的频率。显式传入np.arange(1, 256)或根据目标频率反算尺度结果才可控。5. 避坑与排查CWRU 时频分析中五个让我返工的问题5.1 现象时频图上出现一条横贯整个时间轴的水平亮线原因这通常是工频干扰或传感器固有共振。CWRU 实验台在 60 Hz美国电网频率附近有干扰如果时频图在 60 Hz 或 120 Hz 有一条不随时间变化的亮线基本可以确定是电源干扰。另一个可能是加速度计的谐振频率CWRU 使用的加速度计谐振频率在 2030 kHz但如果你降采样前没做抗混叠滤波谐振能量会混叠到低频。解决先检查信号频谱确认干扰频率。如果是 60 Hz 工频用陷波滤波器滤掉如果是传感器谐振降采样前加抗混叠滤波。不要试图在时频图上「忽略」这条线它会掩盖故障特征频率处的微弱能量。5.2 现象STFT 图上的频率刻度和理论计算的 BPFI 对不上原因采样率设错了。CWRU 数据有 12 kHz 和 48 kHz 两个版本如果你下载的是 48 kHz 数据但代码里写FS12000所有频率都会偏小 4 倍。另一个可能是 .mat 文件里的信号经过了降采样但你没注意到。解决用信号本身的特性反推采样率。CWRU 数据在 0 Hz 附近有明显的轴转频成分负载 0 时约 29.95 Hz负载 1 时约 29.5 Hz。在频谱上找到这个峰用理论值反推采样率。如果对不上检查数据来源。5.3 现象CWT 系数矩阵全是 NaN 或 Inf原因pywt.cwt对信号中的 NaN 或 Inf 零容忍。CWRU 原始数据一般不会有 NaN但如果你做了滤波或降采样边界处理不当可能引入 NaN。另一个可能是信号幅值过大导致浮点溢出但这种情况少见。解决在 CWT 之前加一行signal np.nan_to_num(signal, nan0.0, posinf0.0, neginf0.0)。同时检查滤波器的边界处理scipy.signal.filtfilt比lfilter更安全因为它做了前后向滤波相位不失真。5.4 现象不同负载下的时频图差异巨大无法用同一套参数分析原因负载变化会改变轴转频进而改变故障特征频率。负载 0 时轴转频约 29.95 Hz负载 3 时约 28.5 HzBPFI 从 162 Hz 变到 154 Hz。如果你用固定频率切片去提取特征负载 3 的数据会偏。解决要么按负载分别设置目标频率要么做阶次分析order analysis把频率轴转成轴转频的倍数。阶次分析需要转速信号CWRU 数据里没有直接的转速通道但可以从轴转频成分估算。我一般对每个负载单独算 BPFI不跨负载用同一套参数。5.5 现象时频图看着很漂亮但分类器准确率上不去原因时频图是给人看的分类器需要的是数值特征。直接把时频图当图像喂给 CNN 是一种做法但如果你用的是传统机器学习需要从时频矩阵里提取统计量比如目标频带的能量均值、方差、峰值因子、小波系数的熵。光有图没有特征分类器学不到东西。解决从 STFT 或 CWT 矩阵里提取频带能量比、时频熵、小波包能量等特征。我一般会在 BPFI、BPFO、BSF 及其谐波位置各取一个频带算每个频带的能量占比再加上时频图的整体熵组成特征向量。这套特征在 CWRU 上的分类准确率通常能到 95% 以上比直接用时域统计量高 10 个百分点左右。6. 从时频图到特征向量一个可复用的 CWRU 特征提取函数时频分析本身不是终点终点是拿到能喂给分类器的特征。我把自己反复用的一套特征提取流程整理成函数输入是原始信号和采样率输出是一个固定长度的特征向量。这套流程在 CWRU 的 10 类故障正常 3 种故障类型 × 3 种直径上做过验证配合简单的 SVM 就能得到不错的分类结果。from scipy.stats import kurtosis, skew from scipy.signal import stft import numpy as np def extract_tf_features(signal, fs, nperseg256): 从时频分析结果中提取特征向量 返回: 特征字典 features {} # 1. STFT 频带能量比 f, t, Zxx stft(signal, fsfs, npersegnperseg, noverlapint(nperseg*0.75)) mag np.abs(Zxx) total_energy np.sum(mag**2) 1e-10 # 定义四个频带0-500, 500-1500, 1500-3000, 3000-6000 Hz bands [(0, 500), (500, 1500), (1500, 3000), (3000, 6000)] for i, (lo, hi) in enumerate(bands): mask (f lo) (f hi) band_energy np.sum(mag[mask, :]**2) features[fband_energy_ratio_{i}] band_energy / total_energy # 2. 时频熵 p mag**2 / total_energy features[tf_entropy] -np.sum(p * np.log2(p 1e-10)) # 3. 目标频率切片的统计量 # 以 162 Hz 为例内圈故障特征频率 target_freq 162 freq_mask (f target_freq - 20) (f target_freq 20) if np.any(freq_mask): slice_mag np.mean(mag[freq_mask, :], axis0) features[bpfi_slice_kurtosis] kurtosis(slice_mag) features[bpfi_slice_skew] skew(slice_mag) features[bpfi_slice_peak] np.max(slice_mag) / (np.mean(slice_mag) 1e-10) # 4. 时域统计量作为补充 features[time_kurtosis] kurtosis(signal) features[time_rms] np.sqrt(np.mean(signal**2)) features[time_peak] np.max(np.abs(signal)) return features # 使用示例 feat extract_tf_features(signal, FS) for k, v in feat.items(): print(f{k}: {v:.4f})这段代码的逻辑分四层。第一层是频带能量比把 STFT 幅值平方后按频带求和再除以总能量得到四个比值反映能量在不同频段的分布。轴承故障的冲击能量集中在高频段正常信号能量集中在低频段这个比值有区分度。第二层是时频熵衡量时频平面上能量分布的均匀程度故障冲击会让能量集中熵值下降。第三层是目标频率切片的统计量峰度反映冲击的尖锐程度偏度反映不对称性峰值因子反映冲击的相对强度。第四层是时域统计量作为兜底特征。参数方面nperseg256是起点如果你的信号采样率是 48 kHz窗长要相应加大到 1024 点保持时间窗宽度在 20 ms 左右。频带划分不是固定的根据你的故障特征频率调整如果 BPFI 在 162 Hz第一个频带 0-500 Hz 能覆盖它如果 BPFO 在 100 Hz 左右也在第一个频带里。目标频率切片那里target_freq要按实际故障类型改内圈用 BPFI外圈用 BPFO滚动体用 BSF。这套特征向量长度在 1215 之间取决于频带数量和是否包含目标频率切片。我一般会把所有样本的特征向量拼成矩阵做标准化后用 SVM 或随机森林分类。在 CWRU 的 0 负载数据上10 类故障的分类准确率能到 97% 左右跨负载用 0 负载训练1/2/3 负载测试会降到 85%90%这是 CWRU 数据集的已知难点不是特征提取的问题。最后说一个我自己的习惯每次做完时频分析我会把原始信号、STFT 图、CWT 图、特征向量存成一个字典用joblib.dump保存。这样后面换分类器或调参数时不用重新算时频变换省时间。时频分析的计算量不小尤其是 CWT12 kHz 采样率下 10 秒信号做 256 个尺度的 CWT 要好几秒缓存下来能省不少事。希望帮到你。本文还有配套的精品资源点击获取