简介凯斯西储大学CWRU轴承故障数据集的时频分析文档系统演示短时傅里叶变换STFT与连续小波变换CWT两种主流时频分析工具在故障诊断中的实际应用适合机械健康监测、故障预测方向的研究生和工程技术人员作为入门参考。资源包为单个docx文件大小约1.05MB内容涵盖数据集实验台构成、驱动端/风扇端/基座三个测点的振动信号特点并通过正常信号与0.021英寸内圈、滚珠、外圈三类故障信号的对比试验逐步展示STFT在16、32、64尺度下的时频图变化以及morl、cmor1-1、cmor1.5-2、cgau8四种小波函数在128尺度下的辨识效果最终给出参数选择结论与可复现的Python代码。文档已有1582人学习浏览对希望掌握时频分析参数调优思路、快速复现经典数据集实验的初学者尤为实用。1. 凯斯西储大学轴承故障数据集为什么信号时频分析绕不开它把凯斯西储大学轴承故障数据集的信号时频分析做扎实绝大多数做设备诊断、故障识别和可靠性运维的人都要过这一关。这套数据集来自电机驱动端的振动加速度信号故障类型包含内圈、外圈、滚动体损伤采样率常用 12 kHz 和 48 kHz给了不同负载和转速工况足够用来验证一套诊断算法能不能从振动里把故障特征频率挖出来。很多人一开始只看时域幅值或直接做 FFT结果发现转速一变、负载一变谱线就乱了这时候才意识到必须上时频分析。本文按我平时调这套数据的顺序把方法选型、代码路径、关键参数和踩坑点讲清楚新手照着走能跑通熟手也能在参数边界上做一些检查。2. 信号时频分析的方法选型短时傅里叶、小波与包络谱各自解决什么2.1 先想清楚振动信号的“非平稳”在故障里意味着什么轴承故障信号不是简单的稳态正弦叠加。滚动体经过缺陷点时会产生周期性的冲击这个冲击会激起轴承座、壳体甚至传感器的结构共振共振频率通常在大几千赫兹到上万赫兹而冲击重复频率却只有几十赫兹到一百多赫兹。更麻烦的是电机转速不是恒定的。凯斯西储大学数据虽然看起来是固定转速工况但 0 到 3 HP 负载下电机转速分别约 1797、1772、1750、1730 rpm转频本身就偏移实际运行中还有轻微波动导致故障频率的边带被“抹宽”。这种情况用纯频域做整段信号的 FFT能看出大致峰值但很难同时表达“哪段时间故障冲击强、哪段时间弱”。时频分析的意义就在这里它把一维时间信号展开到二维平面上横轴是时间纵轴是频率颜色代表该时刻该频率的能量强弱。我自己的经验是拿到一段凯斯西储大学数据后不要把“做时频分析”当成一个必须执行的仪式而是先问自己三个问题要观察的是瞬态冲击还是稳定周期振动转速是否可能有轻微波动最终交付物是时频图、特征频率数值还是给分类模型用的特征矩阵这三个问题直接决定你该选短时傅里叶变换、小波变换还是包络谱。顺序搞反了后面调参基本靠玄学。2.2 三种常用时频方法的适用边界和典型参数短时傅里叶变换是最容易上手的方法它把信号切成很多小段对每一段做 FFT再把结果按时间堆叠成谱图。它的代价是时间分辨率和频率分辨率互相制约窗越长频率越准但时间颗粒越粗。小波变换在这一点上更灵活它在高频段用短窗、低频段用长窗适合捕捉像滚动体冲击这类瞬态特征。包络谱并不是严格意义上的时频分布但它和时频分析配在一起非常管用对带通滤波后的信号做 Hilbert 变换取包络再对包络做频谱能把淹没在共振带里的故障特征频率翻出来。方法输出最适合的场景需要盯紧的参数短时傅里叶变换时间-频率-能量谱图转速缓变、长期趋势观察窗长、重叠率、窗函数连续小波变换变分辨率时频图冲击瞬态、突变特征提取小波基、尺度范围、采样周期包络谱一维频谱故障特征频率定位与量化带通范围、滤波器阶数、谱估计点数选型上我先看最终要落地的场景。如果只是给一张“故障存在吗”的图谱短时傅里叶变换最短平快如果要做分类模型特征一般用小波包或者时频图降维如果要给维修人员一个“哪个元件坏了”的判断包络谱往往比彩色图谱更直白。多数故障诊断论文里喜欢把所有方法画一遍但实际项目里没必要。我一般先用短时傅里叶把整个工况扫一眼确定共振频带再用包络谱锁定特征频率最后用小波确认冲击位置。2.3 凯斯西储大学数据里用 FFT 看故障特征频率前先要算对理论值直接画频谱之前先算理论故障特征频率。凯斯西储大学使用的 6205 深沟球轴承参数在公开文件里可以看到但在代码里我不会把系数敲死而是用包含节圆直径、滚动体直径、滚动体数量和接触角的公式计算方便换数据集时直接复用。import numpy as np from scipy.io import loadmat # 读取一个典型文件12kHz 驱动端 0.007 英寸外圈故障负载 0 mat loadmat(12k_Drive_End_B007_0.mat) print([k for k in mat.keys() if DE_time in k]) # 取驱动端振动信号转成 float 一维数组 x mat[X097_DE_time].flatten() fs 12000 # 取同文件里的转速变量避免用名义转频 rpm_key [k for k in mat.keys() if RPM in k] rpm mat[rpm_key[0]][0, 0] if rpm_key else 1797.0 fr rpm / 60.0 # 6205 轴承几何参数单位只要和 Pd/Bd 匹配即可 n_balls 9 Bd 0.3126 # 滚动体直径 Pd 1.537 # 节圆直径 alpha 0.0 # 接触角取 0 度近似值 ftf fr / 2.0 * (1.0 - (Bd / Pd) * np.cos(alpha)) bpfo n_balls * ftf bpfi n_balls * (fr - ftf) bsf fr * (Pd / (2.0 * Bd)) * (1.0 - ((Bd / Pd) * np.cos(alpha)) ** 2) print(f转频 fr {fr:.2f} Hz) print(f外圈 BPFO {bpfo:.2f} Hz) print(f内圈 BPFI {bpfi:.2f} Hz) print(f滚动体 BSF {bsf:.2f} Hz)这段代码必须先跑通因为后面所有时频图上的竖线、横线都来自这几个数值。注意转速变量不一定存在如果文件里没有 RPM就用对应负载表格里的名义值例如负载 0 时用 1797 rpm。单位方面Bd 和 Pd 只要比例一致用英寸、毫米都行因为公式里只出现比值。算出理论值之后再画频谱就不容易把旁边某个不知道来源的谱峰当成故障频率。3. 用 Python 跑通一次完整的信号时频分析从 .mat 文件到彩色谱图3.1 读取凯斯西储大学 .mat 文件的最小代码凯斯西储大学数据集的常用文件是 MATLAB .mat 格式。优先用 scipy 读取因为它不依赖 MATLAB 环境参数也简单。读取后先打印变量名不要直接硬编码某个信号名因为同一种故障在不同文件里的变量名可能不同。from scipy.io import loadmat import numpy as np mat loadmat(12k_Drive_End_IR007_0.mat) keys [k for k in mat.keys() if not k.startswith(__)] print(keys) # 找到驱动端时域信号变量 de_keys [k for k in keys if k.endswith(DE_time)] print(de_keys) x mat[de_keys[0]].flatten() fs 12000 # 如果信号长度明显不对先检查是不是 48kHz 采样 print(len(x), len(x) / fs)这里的DE_time表示驱动端加速度时间序列FE_time是风扇端BA_time是基座或其他测点。文件内部的变量名通常带数字例如X097_DE_time、X105_DE_time数字对应不同故障尺寸和工况。我在实际中经常遇到有人直接拿mat[X097_DE_time]但文件里其实没有这个变量然后报 KeyError。所以打印keys这一步不要省。len(x) / fs能算出信号时长正常文件一般是几秒钟到几十秒如果结果零点几秒说明可能选错了采样率。3.2 短时傅里叶变换的窗口、重叠与谱图输出短时傅里叶变换用 scipy 的stft就能实现。这里最烦的是窗口参数它同时影响图上能分辨的频率间隔和每个时刻的最小持续时间。下面的代码用 1024 点窗、75% 重叠率在 12 kHz 采样率下频率分辨率约为 11.7 Hz足够区分转频及其倍频也不会把冲击抹得太平。from scipy.signal import stft import matplotlib.pyplot as plt # 截取 1 秒避免一次性画太长导致细节丢失 x_seg x[:fs] nperseg 1024 noverlap int(nperseg * 0.75) f, t, Zxx stft( x_seg, fsfs, npersegnperseg, noverlapnoverlap, windowhann, boundaryNone, ) plt.figure(figsize(10, 5)) plt.pcolormesh( t, f, 20.0 * np.log10(np.abs(Zxx) 1e-12), shadinggouraud, cmapturbo, ) plt.ylim(0, 3000) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.colorbar(labeldB) plt.title(CWRU Bearing Vibration STFT) plt.tight_layout() plt.savefig(cwru_stft.png, dpi150)stft返回三个值频率轴f、时间轴t、复数矩阵Zxx。Zxx的形状是频率点数时间点数颜色值我习惯用 dB 而不是线性幅值否则故障冲击的大幅值会把背景全压成蓝色特征根本看不见。ylim(0, 3000)是我在 12 kHz 数据里看轴承特征频率时的常用范围因为 BPFO、BPFI 一般都在 200 Hz 以内但冲击激起的共振边带可能到一两千赫兹如果看 48 kHz 数据共振频带会更高这时需要把纵轴放宽到 5000 Hz 或以上。3.3 小波时频图与边际谱的对照连续小波变换不需要固定窗长它对瞬态冲击的定位更好。PyWavelets 的cwt用起来简单但必须传sampling_period1/fs否则返回的频率轴不是真实赫兹而是伪频率很多人就在这里翻车。import pywt # 取前 0.3 秒约 3600 点突出冲击过程 x_short x[: int(0.3 * fs)] scales np.arange(1, 128) coefs, freqs pywt.cwt( x_short, scales, morl, sampling_period1.0 / fs, ) plt.figure(figsize(10, 5)) plt.imshow( np.abs(coefs), extent[0, x_short.size / fs, freqs[-1], freqs[0]], aspectauto, cmapturbo, ) plt.ylim(0, 3000) plt.xlabel(Time (s)) plt.ylabel(Frequency (Hz)) plt.colorbar(label|CWT coeff|) plt.savefig(cwru_cwt.png, dpi150)pywt.cwt的coefs是二维数组第 0 维对应尺度第 1 维对应时间点。映射到真实频率时extent里的纵轴顺序必须从大频率到小频率因为 imshow 默认把数组第 0 行画在最上面而freqs通常随尺度递增而递减所以我写了[freqs[-1], freqs[0]]。这里的morl是 Morlet 小波适合提取振荡型冲击如果信号里的冲击更窄可以换cmor或gaus但要重新看频率轴。scales选 1 到 128 是经验起点太小会漏掉低频调制太大会把谱图拉得很长且计算量明显上升。4. 参数怎么设才不翻车采样率、窗长、频率范围才是关键4.1 采样率决定你能看到多高频率的故障特征凯斯西储大学数据有 12 kHz 和 48 kHz 两档很多人把这当成“数据越多越好”用 48 kHz 做整个分析结果计算慢、图还糊。其实采样率决定了奈奎斯特频率也就是你能无混叠看到的最高频率。12 kHz 数据能覆盖到 6 kHz足够看 6205 轴承的常见共振和特征频率48 kHz 能看到更高频的局部共振但未必每段信号都值得用那么宽的带宽。分析时我一般先明确如果是确认 BPFO、BPFI 这类几百赫兹的重复频率12 kHz 就够如果要找共振频带或做深度学习模型再考虑 48 kHz。遇到 48 kHz 数据想降采样时不要直接signal[::4]那样会引入混叠必须先用低通滤波器把 6 kHz 以上能量滤掉。from scipy.signal import butter, sosfiltfilt # 从 48kHz 降到 12kHz先防混叠 fs_high 48000 fs_target 12000 btype butter(4, fs_target / (fs_high / 2.0), btypelow, outputsos) x_lp sosfiltfilt(btype, x_high) x_12k x_lp[::4]这段代码里的butter的截止频率写成fs_target / (fs_high / 2.0)也就是 12000 除以 24000 等于 0.5再乘以奈奎斯特频率后得到截止点设在 6000 Hz。注意顺序必须先滤波再抽取否则 6000 Hz 以上的能量会折叠到低频区域把故障特征频率淹没。实际项目中我曾用 48 kHz 数据不加滤波直接降采样结果谱图上 100 Hz 附近出现一堆来历不明的谱峰后来才发现是驱动端一个 5 kHz 共振被混叠回来了。4.2 窗口长度与频率分辨率的取舍短时傅里叶变换的窗口长度直接决定谱图质量。很多教程只说“窗口越长频率分辨率越高”但没告诉你频率分辨率和时间分辨率是跷跷板。以 12 kHz 采样率为例nperseg 对应的频率分辨率是fs / nperseg时间窗长度是nperseg / fs。nperseg频率分辨率单个窗口时长适合场景256约 46.9 Hz约 21 ms冲击定位清晰但特征频率容易糊1024约 11.7 Hz约 85 ms平衡好我常用的起点4096约 2.9 Hz约 341 ms频率精细但瞬态被抹平所以当你发现时频图上故障特征频率附近一片模糊时不要急着换小波先看自己的窗口是不是选得太长反过来如果想确认一次冲击发生在哪个精确时刻窗口要比两段冲击间隔更短。重叠率只影响谱图平滑度和计算量不改变频率分辨率。重叠 50% 是底线75% 带来的平滑效果明显计算开销增加但通常可接受。高频段和低频段需要的信息不同我会先用 1024 点窗跑全图再用 256 点窗放大冲击段两张图配合看。4.3 转频、故障特征频率与边带的对应关系轴承故障信号本质上是调制信号故障冲击以 BPFO 或 BPFI 的周期重复它激起的共振振动被这个周期调制所以在共振频带附近会看到围绕中心频率的一簇边带边带间隔等于故障特征频率。只对整个 0 到 5000 Hz 范围画 STFT很多时候看起来像一片红色条带根本看不出边带。正确做法是先确定共振频带再做带通滤波最后取包络谱。from scipy.signal import butter, sosfiltfilt, hilbert from scipy.signal import periodogram lo, hi 2000.0, 4500.0 sos butter( 4, [lo / (fs / 2.0), hi / (fs / 2.0)], btypeband, outputsos, ) y_band sosfiltfilt(sos, x) env np.abs(hilbert(y_band)) f_env, psd_env periodogram(env, fs, nfft8192) plt.figure(figsize(10, 4)) plt.semilogy(f_env[f_env 500], psd_env[f_env 500]) plt.axvline(bpfo, colorr, linestyle--, labelfBPFO{bpfo:.1f}Hz) plt.xlabel(Frequency (Hz)) plt.ylabel(Envelope PSD) plt.legend() plt.savefig(cwru_envelope.png, dpi150)butter里的截止频率必须除以奈奎斯特频率所以是lo / (fs/2)而不是直接写 2000。带宽范围怎么定我一般回到 STFT 图上看哪个频带的能量随故障冲击明显跳起而不是拍脑袋写死。hilbert的作用是构造解析信号取模得到包络periodogram的nfft决定包络谱的频率分辨率8192 点足够在 500 Hz 内分辨出 0.5 Hz 量级的间隔。得到包络谱后再和理论 BPFO 对比峰值对上了这个故障基本就能锁定。注意带通范围不是标准答案。同一只轴承在不同负载、不同测点下共振频率会移动固定 2 kHz 到 4 kHz 只是起点。每次分析前用 STFT 看一遍共振带再定滤波上下限。5. 避坑凯斯西储大学数据用错的四个高频问题5.1 现象直接 FFT 找不到外圈故障特征频率很多人拿一段外圈故障信号直接做快速傅里叶变换希望在频谱图上看到 BPFO 对应的尖峰结果只看到一堆转频谐波BPFO 附近很小。原因是原始振动信号里转频和它的高次谐波能量远大于故障冲击调制出的边带FFT 是整段时间的平均会把故障冲击的统计特征抹掉。解决方法是不要看原始 FFT改用带通滤波加包络谱。先按第 4 节的流程取共振频带再做 Hilbert 解调故障特征频率在包络谱里通常非常突出。如果包络谱也找不到再检查转速不要用 60 Hz 或 50 Hz 工频去套而是用数据里的 RPM 变量或从振动信号里测量转频。5.2 现象loadmat 报错或者读出来的变量名不对遇到ValueError: Unknown mat file type时大概率是文件不是 MATLAB v5 格式而是 v7.3 格式。v7.3 本质是 HDF5scipy.io.loadmat读不了。解决办法是改读文件而不是硬扛。凯斯西储大学官方旧文件一般不是 v7.3但网上有人重新整理过的副本可能已经是 HDF5。用 h5py 读取时变量名会变成X097_DE_time加上外层数组需要做一层[:]索引。import h5py with h5py.File(cwru_v73.mat, r) as f: de f[X097_DE_time][:].flatten() fs 12000这里f[X097_DE_time][:]读出来的仍然可能是形状(1, N)或(N, 1)所以.flatten()必须保留。另一个隐藏坑是 MATLAB 里字符串变量在 h5py 中会变成引用不要用list(f.keys())直接推断字段名先把所有 key 打印出来再逐个核对。读出来之后不要立刻相信变量名上的数字一定要确认信号长度和采样率是否匹配。5.3 现象时频图里冲击看不见只有平行边带如果 STFT 图上看到一组等间距的平行亮线却找不到局部的冲击亮点一般是因为窗太长时间分辨率不够冲击被平均掉了。平行亮线通常是转频及其谐波不是故障冲击。解决方法是缩小nperseg比如从 1024 改到 256同时把绘图色标改成 dB。如果还是看不到再看是不是信号里工频分量太强把小分量压没了。可以用滤波先去掉转频谐波再画时频图。具体做法是把信号做低通、剔除 300 Hz 以下的能量或者用scipy.signal.detrend去掉趋势项都不行就做一次预白化用 AR 模型残差替代原始信号。诊断冲击类故障最怕的是看图只看颜色深浅不关心时间轴上有没有重复的冲击节拍。5.4 现象小波变换画出来的横纵轴不是期望的“时间-频率”pywt.cwt返回的第二个值是伪频率需要通过sampling_period转换。漏掉这个参数时纵轴显示的是尺度或伪频率无法和理论 BPFO 对应。还有一个常见错误是extent顺序写反导致频率轴高低颠倒。解决方法是先把freqs打印出来确认最大频率和最小频率落在哪个范围然后按我第 3 节代码里的写法用freqs[-1]作为图顶freqs[0]作为图底。如果发现频率轴范围不合理一般需要调整scales的上下限。小波基选择也影响结果morl对振荡信号友好gaus更偏瞬态。不要一张图走天下先看故障信号的冲击形态再换基函数。6. 进阶用故障特征频率做能量带量化再把时频图落到文档里时频分析不能只停留在“画一张好看图”。我一般还会把时频矩阵压缩成一个可比较的故障指数方便放进数据库或技术报告。做法是取 STFT 结果里故障特征频率附近一个频带的能量沿时间方向求中位数再和正常工况的同一个带对比。这样既保留了时频分析的优势又给出一个可量化的阈值。def band_energy(f, Zxx, f_target, bw30.0): lo f_target - bw / 2.0 hi f_target bw / 2.0 idx np.where((f lo) (f hi))[0] return np.mean(np.abs(Zxx[idx, :]) ** 2, axis0) energy_fault band_energy(f, Zxx_fault, bpfo, bw30.0) energy_normal band_energy(f, Zxx_normal, bpfo, bw30.0) ratio np.median(energy_fault) / np.median(energy_normal) print(fBPFO band energy ratio: {ratio:.2f} dB)这段代码里Zxx形状是(频率点, 时间点)np.abs(Zxx[idx, :]) ** 2得到每个时间窗内该频带的瞬时能量mean(axis0)把所有频率点压缩成一个时间序列最后用中位数而不是均值因为冲击信号里偶尔有离群大值均值容易被拉偏。bw30.0表示以故障特征频率为中心上下各 15 Hz这个带宽要大于频率分辨率但又不能大到把相邻倍频包进来。如果做不同工况对比必须保证每段信号的采样率、窗长和带通滤波设置完全一致否则能量比值没有可比性。验证时我再拿帕德劳恩轴承故障数据集交叉跑一遍同一套代码。帕德劳恩轴承故障数据集的试验台结构、轴承几何参数和故障来源都和凯斯西储大学不一样它不是 CWRU 的简单复刻把同一套时频分析脚本换到帕德劳恩数据上如果能量带比值依然能区分正常和故障说明这套方法不是只对某个数据集有效。要注意的是帕德劳恩数据的轴承几何参数并不一定和 6205 相同理论特征频率不能用 CWRU 的那组数字硬套必须先查清滚动体数、节圆直径和滚动体直径再代进公式。最后把分析整理成 docx 文档时我会在每张时频图下面写明采样率、窗长、重叠率、带通滤波范围和基准段文件名没有这些参数任何一张图都只是没法复现的黑匣子。用自己的代码把 CWRU 这套流程跑顺之后再去碰其他轴承数据集会轻松不少希望帮到你。本文还有配套的精品资源点击获取