CWRU轴承故障时频分析：STFT窗口与CWT小波参数选择实战

简介凯斯西储大学CWRU轴承故障数据集是故障诊断领域公认的经典基准数据这份文档以正常、内圈、滚珠、外圈四类振动信号为对象完整演示短时傅里叶变换STFT与连续小波变换CWT的时频分析全过程。文档先介绍1.5KW电机实验台构成以及驱动端、风扇端、基座三类加速度计数据各自的物理含义再逐步讲解信号读取、滑窗取数和参数设置。STFT部分比较重叠比例0.5下16、32、64三种尺度的时频图说明大尺度频率分辨率高、小尺度时间分辨率高的权衡CWT部分对比morl、cmor1-1、cmor1.5-2、cgau8四种小波在尺度128下的效果并进一步考察32、64、128、256尺度对低频特征的响应最终给出选择尺度32和cmor1.5-2的明确依据。文末附可运行的Python代码可直接复现分析并迁移到其他故障数据资源共1个docx文档压缩包约1.05MB图文代码对应、结构清晰。目前已有1582人学习下载适合机械故障诊断、设备健康监测方向的研究生和工程师快速建立时频分析思路。1. 凯斯西储大学轴承故障数据集的时频分析从原始振动信号到能用的时频图做轴承故障诊断的人第一份拿来练手的数据集十有八九是凯斯西储大学CWRU的轴承数据。这套数据驱动端、风扇端、基座三路振动信号齐整故障类型覆盖内圈、滚珠、外圈还带不同故障直径和工作负载算是故障诊断领域的“免检产品”。但数据拿回来后怎么把振动信号变成能说明问题的时频图很多人卡在短时傅里叶变换窗口长度选多少、连续小波变换小波函数用哪个这类细节上。下面用 0.021 英寸内圈、滚珠、外圈故障数据和正常信号做一组完整对比把 STFT 和 CWT 的参数选择逻辑讲透代码拿过去就能跑。适合正在入门轴承故障诊断、或者想快速复现一篇时频分析论文结果的从业者。2. CWRU 数据集选型三个测点与故障样本怎么影响你的对比实验2.1 实验台构成与三个加速度计位置CWRU 实验台的核心部件是一个 1.5kW2 马力电机、一个扭矩传感器/编码器、一个功率测试计外加电子控制器。轴承安装在电机驱动端和风扇端通过在电机壳体的驱动端、风扇端和基座三个位置放置加速度计能同时拿到三路振动信号。这三路信号不是简单的“同一信号不同音量”它们反映的物理过程差别很大。驱动端加速度计测的是电机驱动端的振动信号主要受电机转子旋转和传动系统激励影响适合看轴承故障、齿轮啮合故障这类与旋转轴直接相关的缺陷。风扇端数据受风扇叶片旋转和风扇系统激励主导能反映风扇叶片失衡、风扇端轴承故障。基座数据测的是整个电机系统的振动传递对电机整体不平衡、底座松动这类“全局问题”更敏感。测点加速度计位置主导激励源典型适用故障驱动端 DE电机壳体驱动端转子旋转、传动系统轴承内圈/外圈/滚珠故障、齿轮啮合故障风扇端 FE电机壳体风扇端风扇叶片旋转风扇叶片失衡、风扇端轴承故障基座 BA电机底座整机结构振动电机不平衡、底座松动做时频分析对比实验时我一般优先用驱动端数据原因是驱动端轴承故障的振动能量传递路径最短信噪比最高特征最容易在时频图上显现出来。风扇端数据受风扇背景噪声干扰比较明显基座数据则会把电机自身的电磁振动也叠进来对初学者判断故障特征不够友好。2.2 MAT 文件的读取与信号提取CWRU 数据集发布的是 MAT 格式文件每个文件对应一种工况下的振动记录。读取时要特别注意loadmat读出来的不是数组而是一个 Python 字典真正的振动数据藏在以DE_time、FE_time、BA_time结尾的键里。import numpy as np from scipy.io import loadmat # 读取 MAT 文件data 是字典格式可用 type(data) 确认 data1 loadmat(0_0.mat) # 正常信号 data2 loadmat(21_1.mat) # 0.021 英寸 内圈故障 data3 loadmat(21_2.mat) # 0.021 英寸 滚珠故障 data4 loadmat(21_3.mat) # 0.021 英寸 外圈故障 # DE - drive end accelerometer data 驱动端加速度数据 data_list1 data1[X097_DE_time].reshape(-1) data_list2 data2[X209_DE_time].reshape(-1) data_list3 data3[X222_DE_time].reshape(-1) data_list4 data4[X234_DE_time].reshape(-1) # 划窗取值大多数论文实验窗口大小取 1024 data_list1 data_list1[0:1024] data_list2 data_list2[0:1024] data_list3 data_list3[0:1024] data_list4 data_list4[0:1024]代码逻辑不复杂但有两个关键点一是键名中的X097、X209这类编号对应具体的故障位置和故障尺寸做实验前必须核对清楚否则会把内圈数据当成滚珠数据去分析二是reshape(-1)把二维的 MATLAB 列向量压成一维数组后续stft和pywt.cwt才能直接处理。我一般会顺手打印一下data2.keys()确认字段名因为不同版本的数据集文件命名规则不完全一致直接抄网上的键名偶尔会踩空。2.3 划窗取 1024 点还是全段处理原代码里划窗取了前 1024 个点这是故障诊断论文里最常见的处理方式。原因有两个一是 CWRU 数据在 12kHz 采样率下1024 点对应约 85ms 的振动时长足够覆盖轴承故障特征频率的好几个周期二是后续做训练集时通常会把长信号切成大量短样本1024 点是一个兼容 STFT 窗口设置和 CNN 输入尺寸的折中选择。如果做纯时频分析可视化取 2048 点或 4096 点会更舒服频率分辨率更高时频图上的故障频带更清晰。但如果是为了后续做分类模型1024 点已经够用还能让数据量翻倍。这个取舍取决于你最终要什么看趋势用长窗做样本增强用短窗。3. 短时傅里叶变换窗口长度才是 STFT 的“命门”3.1 STFT 的时间-频率分辨率矛盾短时傅里叶变换的思路很直接信号整体是非平稳的那就用固定时间窗去截取让窗内信号近似平稳再做傅里叶变换。窗口在时间轴上滑动就得到了时间-频率二维分布。核心公式里那个时间窗g(t-α)的宽度直接决定分析效果。这里有个绕不开的矛盾窗口越长频率分辨率越高但时间分辨率越差瞬态冲击在时频图上会被抹平窗口越短时间定位越准但频率分辨率下降相邻故障特征频率可能糊成一团。这就是为什么 STFT 的参数调试看起来像“玄学”——没有绝对最优只有针对具体故障类型的相对最优。3.2 窗口长度 16、32、64、128 的对比实验用 0.021 英寸内圈故障数据做测试重叠比例固定为 0.5窗口长度分别取 16、32、64、128得到四张时频图对比。现象很直观窗口长度频率分辨率时间分辨率内圈故障特征表现16差好频带很宽故障特征淹没在噪声中32一般较好故障频带轮廓开始清晰瞬态冲击可见64较好一般故障频带集中但时间轴上冲击被拉宽128好差频率定位准但时变特征丢失严重实际看时频图窗口长度 16 时整张图几乎看不出内圈故障的周期性冲击窗口长度 64 和 128 时故障特征频率虽然集中但冲击事件在时间轴上变得模糊无法区分是连续磨损还是间歇性冲击。最终选窗口 32是频率分辨率和时间分辨率的折中方案——内圈故障的特征频率在 12kHz 采样率下约为 162Hz 左右窗口 32 足以把这个频率和附近的边带分开同时又能保留每次滚珠通过缺陷时产生的冲击细节。3.3 STFT 核心代码与参数说明from scipy.signal import stft import matplotlib.pyplot as plt import numpy as np # 第一组参数窗口长度 32重叠比例 0.5 window_size 32 overlap 0.5 overlap_samples int(window_size * overlap) # 重叠样本数 16 frequencies1, times1, magnitude1 stft(data_list2, npersegwindow_size, noverlapoverlap_samples) # 第二组参数窗口长度 64 window_size 64 overlap_samples int(window_size * overlap) frequencies2, times2, magnitude2 stft(data_list2, npersegwindow_size, noverlapoverlap_samples) # 第三组参数窗口长度 128 window_size 128 overlap_samples int(window_size * overlap) frequencies3, times3, magnitude3 stft(data_list2, npersegwindow_size, noverlapoverlap_samples) # 第四组参数窗口长度 256 window_size 256 overlap_samples int(window_size * overlap) frequencies4, times4, magnitude4 stft(data_list2, npersegwindow_size, noverlapoverlap_samples) plt.figure(figsize(20, 10), dpi100) plt.subplot(2, 2, 1) plt.pcolormesh(times1, frequencies1, np.abs(magnitude1), shadinggouraud) plt.title(窗口 32-内圈) plt.subplot(2, 2, 2) plt.pcolormesh(times2, frequencies2, np.abs(magnitude2), shadinggouraud) plt.title(窗口 64-内圈) plt.subplot(2, 2, 3) plt.pcolormesh(times3, frequencies3, np.abs(magnitude3), shadinggouraud) plt.title(窗口 128-内圈) plt.subplot(2, 2, 4) plt.pcolormesh(times4, frequencies4, np.abs(magnitude4), shadinggouraud) plt.title(窗口 256-内圈) plt.show()nperseg是 STFT 里的关键参数它决定每个时间窗包含多少个采样点对应上文说的窗口长度。noverlap是相邻窗口的重叠样本数设置为窗口长度的一半也就是重叠比例 0.5能让时频图在时间轴上更平滑避免窗口边界处出现明显接缝。shadinggouraud是pcolormesh的平滑着色模式让时频图的色块过渡更自然方便肉眼观察能量分布。需要说明的是原代码里plt.title写的是“尺度 16”而实际nperseg是 32这个命名错位不影响计算结果但对比实验记录时最好统一口径否则后面写论文容易把自己绕进去。4. 连续小波变换为什么最终选了 cmor1.5-2 而不是 morl4.1 复小波比实小波更适合振动信号连续小波变换和 STFT 最大的区别是STFT 用固定宽度的窗CWT 用可伸缩的小波基函数高频段自动用窄窗、低频段自动用宽窗这就是小波变换的“双尺度”特性。但小波函数的选择直接决定分析结果没有通用最优解。Morlet 小波是复值小波在频率域和时间域都有较好的局部化性质适合处理振动信号这类非平稳信号。cmorComplex Morlet是它的变种通过在复指数上叠加高斯包络实现cgauComplex Gaussian是复数高斯小波更适合近似高斯形状的信号。实小波会把相位信息丢在时频图外复小波能同时保留幅值和相位而振动故障诊断更关心幅值能量分布所以复小波是更稳的选择。4.2 四种小波在尺度 128 下的对比实验尺度长度固定为 128用内圈故障数据分别测试cgau8、morl、cmor1-1、cmor1.5-2四种小波。从时频图看morl的时频图能量分布比较散故障特征频带的边界模糊cmor1-1的带宽参数是 1频带偏窄但能量集中度一般cgau8和cmor1.5-2的故障辨识度明显更高——内圈故障的特征频率及其倍频在时频图上形成清晰的能量脊。这里有个小波参数的概念要先理清cmor1.5-2中1.5 是带宽参数2 是中心频率参数。带宽参数越大小波在频率域的支撑范围越宽频率分辨率越低中心频率参数越大小波振荡越快。对 CWRU 内圈故障这种特征频率明确的信号cmor1.5-2的带宽和中心频率组合出来的时频图故障频带与噪声底之间的对比度最好最终实验选定它做进一步分析。4.3 CWT 尺度序列构造与核心代码import pywt import numpy as np import matplotlib.pyplot as plt # 采样频率CWRU 驱动端数据通常为 12kHz这里做演示用 1024 需按实际数据修改 sampling_rate 1024 sampling_period 1.0 / sampling_rate totalscal 128 # 小波对比组 wavename1 cgau8 fc1 pywt.central_frequency(wavename1) cparam1 2 * fc1 * totalscal scales1 cparam1 / np.arange(totalscal, 0, -1) wavename2 morl fc2 pywt.central_frequency(wavename2) cparam2 2 * fc2 * totalscal scales2 cparam2 / np.arange(totalscal, 0, -1) wavename3 cmor1-1 fc3 pywt.central_frequency(wavename3) cparam3 2 * fc3 * totalscal scales3 cparam3 / np.arange(totalscal, 0, -1) wavename4 cmor1.5-2 fc4 pywt.central_frequency(wavename4) cparam4 2 * fc4 * totalscal scales4 cparam4 / np.arange(totalscal, 0, -1) # 连续小波变换返回系数矩阵和对应的频率序列 coefficients1, frequencies1 pywt.cwt(data_list2, scales1, wavename1, sampling_period) coefficients2, frequencies2 pywt.cwt(data_list2, scales2, wavename2, sampling_period) coefficients3, frequencies3 pywt.cwt(data_list2, scales3, wavename3, sampling_period) coefficients4, frequencies4 pywt.cwt(data_list2, scales4, wavename4, sampling_period) # 取系数矩阵绝对值作为时频幅值 amp1 abs(coefficients1) amp2 abs(coefficients2) amp3 abs(coefficients3) amp4 abs(coefficients4) # 生成时间轴 t np.linspace(0, 1.0 / sampling_rate, sampling_rate, endpointFalse) plt.figure(figsize(20, 10)) plt.subplot(2, 2, 1) plt.contourf(t, frequencies1, amp1, cmapjet) plt.title(内圈-cgau8) plt.subplot(2, 2, 2) plt.contourf(t, frequencies2, amp2, cmapjet) plt.title(内圈-morl) plt.subplot(2, 2, 3) plt.contourf(t, frequencies3, amp3, cmapjet) plt.title(内圈-cmor1-1) plt.subplot(2, 2, 4) plt.contourf(t, frequencies4, amp4, cmapjet) plt.title(内圈-cmor1.5-2) plt.show()尺度序列构造是这段代码的核心。pywt.central_frequency(wavename)返回小波函数的中心频率cparam 2 * fc * totalscal是尺度常数项。np.arange(totalscal, 0, -1)生成从 128 递减到 1 的整数序列cparam除以这个序列后scales数组从小到大排列——前面是小尺度对应高频后面是大尺度对应低频。这样做的好处是让pywt.cwt返回的frequencies从高频到低频规则排列画出来的时频图频率轴是顺的。sampling_period必须显式传参它是采样周期采样率的倒数直接决定frequencies序列的物理单位。如果这里漏定义代码会在运行时直接报NameError这是新手最容易翻车的地方。5. CWRU 数据集实战避坑五个让人翻车的细节5.1 loadmat 读出来的字典结构不等于数组现象loadmat(21_1.mat)后直接print(data)发现输出一个巨大的 dict找不到振动数据在哪甚至data[X209_DE_time]报 KeyError。原因MAT 文件被读成字典结构键名是 MATLAB 变量名而且不同版本数据集的键名可能不同。部分文件还带__header__、__version__、__globals__这类元数据键新手容易被带偏。解决先执行data.keys()把所有键列出来确认振动字段名后再取数。字段名通常以DE_time、FE_time、BA_time结尾取出来后再reshape(-1)转成一维数组。我自己的习惯是写一小段脚本打印所有键名做一次“数据体检”再进入分析流程。5.2 采样频率拿错导致时间轴和频率轴全错现象时频图画出来后频率轴范围和理论值差了一个数量级时间轴长度也不对。原因原代码里sampling_rate 1024是演示值但 CWRU 驱动端数据的实际采样率通常是 12kHz部分工况是 48kHz。频率轴跟着采样率走采样率错整个时频图的物理意义就废了。解决分析前先确认所读文件对应的采样率。CWRU 官方文档里驱动端 12kHz 采样是主力配置风扇端和基座端数据也是 12kHz。如果是 48kHz 版本的数据文件sampling_rate改成 48000t np.linspace(0, 1.0/sampling_rate, sampling_rate, endpointFalse)里的长度也跟着改不能照抄 1024 的配置。5.3 sampling_period 未定义导致 CWT 直接报错现象运行 CWT 代码时抛出NameError: name sampling_period is not defined或者frequencies返回全空。原因pywt.cwt需要接收采样周期参数sampling_period但很多网上流传的代码片段只写pywt.cwt(data, scales, wavelet)漏传这个参数或者定义了sampling_rate却忘了定义sampling_period。解决在 CWT 代码块前补一行sampling_period 1.0 / sampling_rate并把完整参数传给pywt.cwt。这是个一秒钟能修好的问题但第一次跑的时候确实会卡很久因为报错信息不会直接告诉你是哪个参数缺失。5.4 0.021 英寸滚珠故障的区分度最低现象四种小波函数对比中滚珠故障的时频图总是比内圈和外圈模糊能量分布散特征频率不明显。原因滚珠故障的振动激励路径最长缺陷经过承载区时冲击能量被滚动体分散加上保持架的随机滑动故障特征频率的周期性比内圈、外圈弱很多。这不是代码问题是物理本质决定的。解决做滚珠故障分析时不要只对比内圈数据的参数必须单独用滚珠数据做一组尺度对比实验如果多个小波函数下滚珠故障仍然模糊优先检查尺度长度是否足够试着把totalscal从 128 提到 256让低频段有更细的划分。如果手里有帕德劳恩轴承故障数据集也可以拿来做交叉验证那个数据集的滚珠故障样本更丰富对比结果更有说服力。5.5 文件命名和故障位置对应关系容易记错现象分析做完才发现代码里标注的“内圈”数据实际上是“滚珠”数据整套对比实验白跑。原因CWRU 文件名编号规则不直观21_1、21_2、21_3并不是简单的故障类型顺序MAT 文件内部键名里的X209、X222、X234也各有对应关系。看网上代码直接复制容易把映射关系张冠李戴。解决开始实验前先做一次“文件清单核对”把文件名、内部键名、故障类型、故障直径做成一张对照表确认无误后再执行分析。下面是这次实验的对照关系文件名内部字段名故障直径故障类型0_0.matX097_DE_time正常无故障21_1.matX209_DE_time0.021 英寸内圈21_2.matX222_DE_time0.021 英寸滚珠21_3.matX234_DE_time0.021 英寸外圈6. 把时频图变成可用的故障特征批量提取与保存技巧6.1 从时频矩阵到能量特征时频图看着直观但喂给分类模型前必须转成数值特征。一个实用做法是沿时间轴对时频幅值矩阵做均值池化得到“平均频谱”再取特征频段的能量占比作为判别指标。# 基于 cmor1.5-2 小波系数提取能量特征 coeff, freqs pywt.cwt(data_list2, scales4, cmor1.5-2, sampling_period) amp abs(coeff) # 沿时间轴做均值池化得到每个频率点的平均幅值 mean_amp np.mean(amp, axis1) # 计算频段能量取频率在 100-300Hz 范围内的能量占比 band_mask (freqs 100) (freqs 300) band_energy np.sum(mean_amp[band_mask] ** 2) total_energy np.sum(mean_amp ** 2) energy_ratio band_energy / total_energy print(f特征频段能量占比: {energy_ratio:.4f})这里的逻辑是故障特征频率及其倍频集中在特定频段能量占比越高故障特征越明显。mean_amp把二维时频矩阵压缩成一维频谱band_mask圈定特征频段energy_ratio就是后续分类的输入特征。这个特征比直接用原始振动信号做统计量稳定得多因为它天然去掉了时间轴的相位影响。6.2 批量保存时频图沉淀自己的可视化模板做完一轮实验后最该做的是把代码固化成模板。下面这段代码可以把四种故障的 CWT 时频图批量保存到本地下次换数据集只需要改文件路径和字段名。import os import pywt import numpy as np from scipy.io import loadmat import matplotlib.pyplot as plt output_dir cwt_results os.makedirs(output_dir, exist_okTrue) cases [ {file: 0_0.mat, key: X097_DE_time, label: normal}, {file: 21_1.mat, key: X209_DE_time, label: inner}, {file: 21_2.mat, key: X222_DE_time, label: ball}, {file: 21_3.mat, key: X234_DE_time, label: outer}, ] for case in cases: data loadmat(case[file])[case[key]].reshape(-1)[:1024] coeff, freqs pywt.cwt(data, scales4, cmor1.5-2, sampling_period) plt.figure(figsize(12, 4)) plt.contourf(t, freqs, abs(coeff), cmapjet) plt.title(case[label]) plt.colorbar() plt.savefig(os.path.join(output_dir, f{case[label]}_cwt.png), dpi150, bbox_inchestight) plt.close()bbox_inchestight会裁掉多余白边dpi150保证图片清晰度够写报告用。plt.close()是必须的否则循环里会累积出大量内存中的 figure 对象跑几十个文件就可能卡死。这套流程跑通之后我基本告别了每次实验都从零调参的状态。从那以后做任何数据集我都会强制先跑一遍“数据体检 文件核对 参数验证”三步再进入正式的时频分析对比。CWRU 数据的坑不算多但每个坑都踩过一次之后你会在其他数据集的时频分析里少走很多弯路。希望帮到你。本文还有配套的精品资源点击获取