轴承故障诊断实战：从振动信号到Python代码的完整分析流程

轴承故障诊断实战从振动信号到Python代码的完整分析流程在工业设备维护领域轴承作为旋转机械的核心部件其健康状态直接影响设备运行效率与安全性。传统的人工巡检方式已难以满足现代工业对故障预警的实时性需求而基于振动信号分析的智能诊断技术正成为工程师手中的利器。本文将带您从原始振动信号出发通过Python代码逐步拆解轴承故障诊断的全流程涵盖数据预处理、特征提取、算法实现等关键环节最终输出可落地的诊断结论。1. 数据准备与预处理轴承振动数据的质量直接影响后续分析效果。以凯斯西储大学CWRU轴承数据集为例其采样频率为10kHz包含正常、内圈故障、外圈故障、滚珠故障等多种状态。原始数据通常存在以下问题基线漂移传感器零点偏移导致信号整体偏离零轴噪声干扰环境噪声与电磁干扰叠加在有效信号上量纲差异不同测点间的幅值范围差异显著1.1 数据加载与可视化使用scipy.io加载MATLAB格式的原始数据并通过时域波形初步观察信号特征import numpy as np import matplotlib.pyplot as plt from scipy.io import loadmat # 加载数据文件 data loadmat(97.mat) signal data[X097_DE_time].reshape(-1) # 驱动端振动信号 # 绘制时域波形 plt.figure(figsize(12,4)) plt.plot(np.arange(len(signal))/10000, signal) # 时间轴转换为秒 plt.xlabel(Time (s)) plt.ylabel(Amplitude (g)) plt.title(Raw Vibration Signal) plt.grid() plt.show()1.2 关键预处理步骤零均值化处理是后续分析的必备环节消除信号直流分量对频域分析的影响def remove_mean(signal): return signal - np.mean(signal) signal_zm remove_mean(signal)滤波去噪可根据故障特征频率范围选择带通滤波器。例如滚动轴承故障特征频率通常在几百Hz到几kHz之间from scipy.signal import butter, filtfilt def bandpass_filter(signal, lowcut, highcut, fs, order4): nyq 0.5 * fs low lowcut / nyq high highcut / nyq b, a butter(order, [low, high], btypeband) return filtfilt(b, a, signal) signal_filt bandpass_filter(signal_zm, lowcut500, highcut5000, fs10000)提示预处理阶段建议保存中间结果便于回溯分析各步骤对最终诊断的影响2. 特征提取与指标计算从振动信号中提取有效的故障特征是诊断的核心环节。下表对比了常用时域特征指标的计算公式与物理意义特征名称计算公式诊断意义均方根值 (RMS)$\sqrt{\frac{1}{N}\sum_{i1}^N x_i^2}$反映信号能量水平峰值指标 (CF)$\frac{\max(x峭度 (Kurtosis)$\frac{\frac{1}{N}\sum_{i1}^N (x_i-\mu)^4}{\sigma^4}$敏感早期故障2.1 时域特征计算Python实现关键时域特征计算函数def calc_time_features(signal): features {} features[RMS] np.sqrt(np.mean(signal**2)) features[Peak] np.max(np.abs(signal)) features[Kurtosis] np.mean((signal-np.mean(signal))**4) / np.std(signal)**4 features[CrestFactor] features[Peak] / features[RMS] return features time_features calc_time_features(signal_filt)2.2 频域特征提取傅里叶变换将信号转换到频域揭示周期性故障特征from scipy.fft import fft def compute_fft(signal, fs): n len(signal) freq np.fft.fftfreq(n, d1/fs)[:n//2] fft_vals np.abs(fft(signal))[:n//2] * 2/n return freq, fft_vals freq, fft_vals compute_fft(signal_filt, fs10000)包络分析对检测轴承局部缺陷尤为有效通过希尔伯特变换提取信号包络from scipy.signal import hilbert def envelope_analysis(signal, fs): analytic_signal hilbert(signal) envelope np.abs(analytic_signal) env_freq, env_fft compute_fft(envelope, fs) return env_freq, env_fft3. 故障特征频率理论计算轴承各部件故障特征频率与几何参数相关计算公式如下内圈故障频率 (BPFI): $$f_i \frac{N_b}{2} \times f_r \times \left(1 \frac{B_d}{P_d} \cos\phi\right)$$外圈故障频率 (BPFO): $$f_o \frac{N_b}{2} \times f_r \times \left(1 - \frac{B_d}{P_d} \cos\phi\right)$$滚珠故障频率 (BSF): $$f_b \frac{P_d}{2B_d} \times f_r \times \left[1 - \left(\frac{B_d}{P_d}\cos\phi\right)^2\right]$$Python实现这些特征频率计算def bearing_frequencies(rpm, nb, bd, pd, contact_angle0): fr rpm / 60.0 # 转频(Hz) bpfi nb/2 * fr * (1 (bd/pd)*np.cos(contact_angle)) bpfo nb/2 * fr * (1 - (bd/pd)*np.cos(contact_angle)) bsf pd/(2*bd) * fr * (1 - (bd/pd * np.cos(contact_angle))**2) return {BPFI:bpfi, BPFO:bpfo, BSF:bsf, FTF:0.4*fr} # FTF为保持架频率近似值 # CWRU 6205-2RS轴承参数示例 freqs bearing_frequencies(rpm1797, nb9, bd7.94, pd39.04)4. 故障模式识别与诊断结合特征指标与频率分析结果建立诊断逻辑4.1 基于规则的门限诊断def rule_based_diagnosis(features, freqs, env_spectrum): diagnosis [] # 峭度指标判断 if features[Kurtosis] 5: diagnosis.append(存在明显冲击特征) # 包络谱峰值检测 env_peaks np.where(env_spectrum np.mean(env_spectrum)3*np.std(env_spectrum))[0] matched_freqs [] for peak in env_peaks: for name, theory_freq in freqs.items(): if abs(peak - theory_freq) 5: # 允许±5Hz误差 matched_freqs.append(name) if BPFI in matched_freqs: diagnosis.append(内圈故障特征明显) elif BPFO in matched_freqs: diagnosis.append(外圈故障特征明显) return diagnosis if diagnosis else [无明显故障特征]4.2 机器学习方法应用对于更复杂的工况可采用机器学习模型。以下示例使用随机森林分类器from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split # 特征工程组合多个时频域特征 def extract_all_features(signals, fs): features_list [] for sig in signals: # 时域特征 feat calc_time_features(sig) # 频域特征 _, fft_vals compute_fft(sig, fs) feat[FFT_peak] np.max(fft_vals) feat[FFT_mean] np.mean(fft_vals) # 包络谱特征 _, env_fft envelope_analysis(sig, fs) feat[Env_peak] np.max(env_fft) features_list.append(list(feat.values())) return np.array(features_list) # 模型训练与评估 X extract_all_features(training_signals, fs10000) y training_labels # 预准备的标签数据 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3) model RandomForestClassifier(n_estimators100) model.fit(X_train, y_train) print(f模型准确率: {model.score(X_test, y_test):.2f})5. 工程实践中的优化技巧在实际项目中以下几个技巧可显著提升诊断效果多传感器数据融合同时分析驱动端、风扇端、水平/垂直方向的振动信号def fuse_sensors(signals, methodaverage): if method average: return np.mean(signals, axis0) elif method max: return np.max(signals, axis0)变工况适配当转速或负载变化时采用阶次分析代替固定频率分析def order_analysis(signal, tachometer, fs): from scipy.interpolate import interp1d # 基于转速脉冲信号重采样 angle np.cumsum(tachometer)/fs * 360 f interp1d(angle, signal, kindlinear) angle_uniform np.linspace(0, angle[-1], len(angle)) return f(angle_uniform)趋势监测建立特征指标的历史变化曲线比单次检测更可靠def trend_monitoring(feature_history, window30): from statsmodels.tsa.seasonal import STL # 使用STL分解趋势项 stl STL(feature_history, periodwindow) res stl.fit() return res.trend在完成多个工业现场项目后发现最有效的诊断策略往往是结合时域峭度指标、包络谱分析和机器学习模型的混合方法。当设备转速波动较大时重采样技术能显著提升特征提取的稳定性。