用Python和Matlab搞定东南大学齿轮箱数据集：从数据读取到故障分类实战
东南大学齿轮箱数据集实战：Python与Matlab双视角解析
引言：机械故障诊断的数据科学之旅
机械故障诊断领域正经历着从传统经验判断向数据驱动决策的转变。东南大学齿轮箱数据集作为国内知名的开放研究资源，包含了轴承和齿轮两大类共10种故障状态，为算法验证提供了宝贵素材。但面对8通道振动信号、多工况CSV文件时，许多研究者常陷入"数据在手，无从下手"的困境。
本文将采用工程化思维，以Python生态为主力工具链，Matlab为辅助参考，构建从原始振动信号到分类模型的完整分析流程。不同于简单的代码展示，我们会重点探讨：
- 多源异构振动信号的标准化处理方法
- 时频域特征工程的构建逻辑
- 分类模型在机械信号中的适配技巧
- 两种编程环境的协同工作模式
无论您是刚开始接触机械信号分析的在校学生，还是需要快速验证算法的工程人员，这套方法论都能帮助您建立系统化的分析框架。
1. 数据环境搭建与预处理
1.1 跨平台数据加载方案
东南大学数据集采用CSV格式存储，但包含特殊配置信息头。我们推荐以下跨语言处理方案：
Python方案（Pandas优先）
import pandas as pd def load_gear_data(filepath): # 跳过前两行配置信息 raw_data = pd.read_csv(filepath, skiprows=2, header=None) # 8通道信号标准化命名 columns = ['motor_vib', 'planet_x', 'planet_y', 'planet_z', 'motor_torque', 'reducer_x', 'reducer_y', 'reducer_z'] return raw_data.set_axis(columns, axis=1)Matlab对照方案
function data = loadGearData(filename) opts = detectImportOptions(filename); opts.DataLines = [3 Inf]; data = readtable(filename, opts); data.Properties.VariableNames = {'motor_vib', 'planet_x',... 'planet_y', 'planet_z', 'motor_torque',... 'reducer_x', 'reducer_y', 'reducer_z'}; end注意：原始数据采样频率为5120Hz，建议在加载阶段就记录该元数据，后续分析会频繁用到这个关键参数。
1.2 多文件批量处理框架
面对数十个CSV文件，需要建立系统化的批处理框架：
from pathlib import Path def build_dataset(root_dir): dataset = {} for csv_file in Path(root_dir).glob('*.csv'): # 从文件名解析故障类型和工况 parts = csv_file.stem.split('_') fault_type = parts[0] condition = parts[1] if len(parts)>1 else 'default' # 统一数据存储结构 if fault_type not in dataset: dataset[fault_type] = {} dataset[fault_type][condition] = load_gear_data(csv_file) return dataset关键处理技巧：
- 使用文件命名规范自动提取元信息
- 构建嵌套字典保存原始信号
- 保留完整的信号通道关系
2. 振动信号可视化与探索
2.1 多通道信号对比分析
典型的多通道信号可视化方案：
import matplotlib.pyplot as plt def plot_multi_channel(data, start=0, duration=0.1): samples = int(duration * 5120) # 计算采样点数 subset = data.iloc[start:start+samples] fig, axes = plt.subplots(4, 2, figsize=(12, 8)) for idx, col in enumerate(data.columns): ax = axes[idx//2, idx%2] ax.plot(subset[col]) ax.set_title(f'Channel {col}') plt.tight_layout() return fig信号观察要点：
- 电机振动（motor_vib）通常包含丰富的高频成分
- 行星齿轮XYZ三向振动呈现明显相关性
- 扭矩信号往往具有不同的量纲需要单独处理
2.2 时频域联合分析
结合快速傅里叶变换（FFT）和短时傅里叶变换（STFT）的综合观察：
from scipy import signal def time_freq_analysis(channel_data, fs=5120): # 时域特征 time_features = { 'peak': channel_data.max(), 'rms': np.sqrt(np.mean(channel_data**2)), 'kurtosis': channel_data.kurtosis() } # 频域特征 freqs, psd = signal.welch(channel_data, fs) dominant_freq = freqs[np.argmax(psd)] return {'time': time_features, 'main_freq': dominant_freq}提示：齿轮故障常在啮合频率及其谐波处表现出异常，轴承故障则多在外圈/内圈特征频率附近出现峰值。
3. 特征工程构建策略
3.1 时域特征组合
推荐提取的时域特征集：
| 特征类别 | 具体指标 | 物理意义 |
|---|---|---|
| 幅值特征 | 峰值、峰峰值、RMS | 振动能量水平 |
| 波形指标 | 偏度、峭度、波形因子 | 信号分布特性 |
| 脉冲指标 | 脉冲因子、裕度因子 | 冲击成分强度 |
| 相关性特征 | 通道间互相关系数 | 系统耦合程度 |
Python实现示例：
def extract_time_features(data_frame): features = {} for col in data_frame.columns: x = data_frame[col].values features.update({ f'{col}_peak': np.max(x), f'{col}_kurtosis': stats.kurtosis(x), f'{col}_crest': np.max(x)/np.sqrt(np.mean(x**2)) }) return features3.2 频域特征提取
关键频域特征包括：
频谱重心：反映能量分布位置
def spectral_centroid(psd, freqs): return np.sum(freqs*psd) / np.sum(psd)频带能量比：划分5个典型频段计算能量占比
% Matlab实现 function ratio = band_energy_ratio(psd, f_edges) total = sum(psd); ratio = []; for i = 1:length(f_edges)-1 mask = (f >= f_edges(i)) & (f < f_edges(i+1)); ratio(end+1) = sum(psd(mask)) / total; end end谐波成分检测：针对齿轮特征频率的谐波分析
4. 故障分类模型构建
4.1 特征选择与数据集构建
典型特征矩阵构建流程：
from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler # 构建特征矩阵 def build_feature_matrix(dataset): X, y = [], [] for fault_type, conditions in dataset.items(): for condition, data in conditions.items(): # 合并时域频域特征 features = extract_time_features(data) features.update(extract_freq_features(data)) X.append(list(features.values())) y.append(fault_type) return np.array(X), np.array(y) # 标准化管道 preprocessor = make_pipeline( StandardScaler(), PCA(n_components=0.95) )4.2 模型训练与评估
推荐模型对比表：
| 模型类型 | 优点 | 注意事项 |
|---|---|---|
| SVM | 小样本效果好，泛化能力强 | 需要仔细调参 |
| 随机森林 | 自动特征选择，抗噪声 | 可能过拟合 |
| 1D CNN | 自动特征学习，端到端 | 需要大量数据 |
| XGBoost | 处理不平衡数据优秀 | 超参数较多 |
Python实现示例：
from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score # 随机森林实现 model = RandomForestClassifier( n_estimators=200, max_depth=10, class_weight='balanced' ) cv_scores = cross_val_score( model, X_preprocessed, y, cv=5, scoring='accuracy' ) print(f"Mean CV Accuracy: {cv_scores.mean():.2%}")4.3 模型解释与故障分析
SHAP值分析示例：
import shap # 训练最终模型 model.fit(X_train, y_train) # 解释模型预测 explainer = shap.TreeExplainer(model) shap_values = explainer.shap_values(X_test) # 可视化特定样本分析 shap.force_plot( explainer.expected_value[0], shap_values[0][0,:], X_test[0,:], feature_names=feature_names )关键分析维度：
- 各通道振动信号的贡献度排序
- 时域与频域特征的相对重要性
- 不同故障类型的判别特征差异
5. 工程实践中的进阶技巧
5.1 数据增强策略
机械信号特有的数据增强方法：
时域窗口滑动：通过重叠采样增加样本量
def sliding_window(data, window_size, step): n_samples = len(data) return [data[i:i+window_size] for i in range(0, n_samples-window_size, step)]噪声注入：添加符合实际工况的高斯噪声
def add_industrial_noise(signal, snr_db=20): rms_signal = np.sqrt(np.mean(signal**2)) noise = np.random.normal(0, rms_signal/(10**(snr_db/20)), len(signal)) return signal + noise
5.2 跨语言协作模式
Python与Matlab混合编程方案：
数据交换格式：
- 使用HDF5(.h5)实现大型数组交换
- 通过MAT文件接口
scipy.io.savemat/loadmat
实时通信：
% Matlab作为计算引擎 eng = py.matlab.engine.start_matlab() result = eng.sqrt(4.0)性能分工：
- Matlab处理信号处理专用算法
- Python负责机器学习流程
5.3 实际部署考量
边缘设备部署优化策略：
特征压缩：使用PCA保留95%方差
from sklearn.decomposition import PCA pca = PCA(n_components=0.95) X_compressed = pca.fit_transform(X)模型轻量化：
- 决策树剪枝
- 神经网络量化
实时性保障：
- 固定长度滑动窗口
- 增量特征计算
常见问题解决方案
Q1 数据量不足怎么办？
- 采用迁移学习，复用ImageNet预训练模型
- 使用生成对抗网络(GAN)进行数据增强
Q2 不同工况数据分布不一致？
- 添加工况标签作为额外特征
- 采用领域自适应(Domain Adaptation)技术
Q3 如何选择采样长度？
- 齿轮分析至少包含5个完整啮合周期
- 轴承故障建议0.5-1秒时长
Q4 特征重要性如何评估？
- 排列重要性(Permutation Importance)
- 基于SHAP值的全局解释
延伸方向与资源推荐
进阶学习路径：
- 时频分析：小波变换、Hilbert-Huang变换
- 深度学习：1D CNN、LSTM、Transformer
- 异常检测：One-Class SVM、Autoencoder
参考工具库：
- Python: PyBearings, SciPy Signal
- Matlab: Predictive Maintenance Toolbox
- 开源项目: SKF @ptitude Analyst
公开数据集扩展：
- NASA轴承数据集
- 凯斯西储大学数据集
- Paderborn大学轴承数据集