用Python和Matlab处理辛辛那提IMS轴承数据：从ASCII文件到故障特征提取的完整流程

工业轴承故障诊断实战Python与Matlab双视角下的辛辛那提IMS数据处理指南在工业设备预测性维护领域轴承故障诊断一直是研究热点。辛辛那提大学IMS轴承数据集作为经典公开数据包含了轴承从健康状态到完全失效的全生命周期振动信号为算法验证提供了宝贵资源。但原始数据以非标准ASCII格式存储且文件命名存在特殊性直接使用存在诸多障碍。本文将带你用Python和Matlab双工具链系统解决从原始数据解析到时频特征提取的全流程问题。1. 数据准备与环境配置1.1 数据集获取与初步检查辛辛那提IMS数据集包含三个子集记录了三组轴承在不同工况下的振动信号。每个数据文件对应1秒的采样数据20kHz采样率共20480个点。在开始处理前建议先创建如下目录结构ims_data/ ├── raw/ # 存放原始ASCII文件 ├── processed/ # 存放处理后的结构化数据 ├── notebooks/ # Jupyter分析笔记 └── scripts/ # 处理脚本使用Python进行批量文件检查时可以运行以下代码快速统计文件数量import os def count_files(dataset_path): for root, dirs, files in os.walk(dataset_path): print(f{root}: {len(files)} files) # 示例用法 count_files(ims_data/raw/dataset1)1.2 工具库安装Python环境需要以下核心库pip install numpy pandas scipy matplotlib scikit-learnMatlab用户需要确保已安装Signal Processing ToolboxStatistics and Machine Learning Toolbox2. ASCII数据解析实战2.1 Python解析方案原始ASCII文件每行包含一个采样点的多个通道数据。使用Pandas可以高效读取import pandas as pd def read_ims_file(filepath, channels8): # 自动检测分隔符并读取 df pd.read_csv(filepath, headerNone, delim_whitespaceTrue) # 规范列名 df.columns [fchannel_{i1} for i in range(channels)] return df # 示例读取单个文件 sample_data read_ims_file(ims_data/raw/dataset1/bearing1_0001.txt)对于批量处理建议使用多进程加速from multiprocessing import Pool def process_batch(file_list): with Pool() as p: results p.map(read_ims_file, file_list) return pd.concat(results, ignore_indexTrue)2.2 Matlab解析优化针对原始博文提到的文件名冲突问题可以改进为function data readIMSFiles(folderPath) files dir(fullfile(folderPath, *.txt)); allData cell(length(files), 1); parfor i 1:length(files) filePath fullfile(files(i).folder, files(i).name); allData{i} dlmread(filePath); end data vertcat(allData{:}); end注意Matlab并行计算需要Parallel Computing Toolbox支持3. 数据质量控制与预处理3.1 异常数据检测轴承振动数据常见问题包括传感器失电导致的零值段采样异常造成的突变值通道间相位偏移Python检测示例def check_data_quality(df): # 零值检测 zero_counts (df 0).mean() # 幅值范围检测 stats df.describe().loc[[min, max, mean]] return { zero_ratio: zero_counts, statistics: stats }3.2 数据标准化方案不同通道间可能存在量纲差异推荐使用Robust Scalingfrom sklearn.preprocessing import RobustScaler scaler RobustScaler() normalized_data scaler.fit_transform(raw_data)4. 时频域特征工程4.1 时域特征提取基础时域特征包括特征名称计算公式物理意义峰值max(x均方根值sqrt(mean(x²))能量水平峭度mean(x⁴)/std(x)⁴冲击成分敏感度脉冲因子peak/RMS极端值相对水平Python实现示例from scipy.stats import kurtosis def extract_time_features(signal): features { peak: np.max(np.abs(signal)), rms: np.sqrt(np.mean(signal**2)), kurtosis: kurtosis(signal), crest_factor: np.max(np.abs(signal)) / np.sqrt(np.mean(signal**2)) } return features4.2 频域分析方法FFT分析是基础但更推荐使用包络谱分析function [envSpectrum] envelope_analysis(signal, fs) % 希尔伯特变换获取包络 analytic hilbert(signal); envelope abs(analytic); % 计算包络谱 N length(envelope); f (0:N-1)*(fs/N); envSpectrum abs(fft(envelope-mean(envelope))); end对应Python版本from scipy.signal import hilbert def envelope_spectrum(signal, fs): analytic hilbert(signal) envelope np.abs(analytic) spectrum np.abs(np.fft.fft(envelope - np.mean(envelope))) freqs np.fft.fftfreq(len(envelope), 1/fs) return freqs, spectrum5. 特征可视化与早期故障检测5.1 趋势特征可视化绘制关键特征随时间变化曲线import matplotlib.pyplot as plt def plot_trend_features(feature_df): fig, axes plt.subplots(2, 2, figsize(12, 8)) features [rms, peak, kurtosis, crest_factor] for ax, feat in zip(axes.ravel(), features): ax.plot(feature_df[feat]) ax.set_title(feat.upper()) plt.tight_layout() return fig5.2 基于统计的过程控制使用控制图检测异常def spc_control_chart(data, feature, window30): rolling_mean data[feature].rolling(window).mean() rolling_std data[feature].rolling(window).std() plt.figure(figsize(10, 4)) plt.plot(data[feature], alpha0.5) plt.plot(rolling_mean, r) plt.fill_between(data.index, rolling_mean - 3*rolling_std, rolling_mean 3*rolling_std, colorr, alpha0.1) plt.title(f{feature} SPC Chart)6. 工程实践中的经验技巧内存优化处理大规模IMS数据时建议使用Python的Dask或Matlab的memmap处理超出内存的数据按轴承通道分别存储减少单次加载数据量特征选择策略早期故障阶段重点关注高频段能量和峭度指标发展期故障跟踪1-3倍轴承特征频率的能量变化严重故障期监测整体振动水平跨平台协作建议使用HDF5格式在Python和Matlab间交换数据统一时间戳格式建议采用Unix时间戳# Python保存HDF5示例 import h5py with h5py.File(bearing_data.h5, w) as f: f.create_dataset(vibration, dataprocessed_data) f.attrs[sampling_rate] 20000% Matlab读取HDF5示例 info h5info(bearing_data.h5); data h5read(bearing_data.h5, /vibration); fs h5readatt(bearing_data.h5, /, sampling_rate);实际项目中我们发现数据集1的轴承3在文件#1450左右开始出现内圈故障特征这时峭度值会突然增大2-3个数量级而RMS值的变化相对滞后约50个文件。这种时域特征的差异组合使用能显著提高早期故障检测的准确率。