基于几何流形学习的轴承故障诊断Python实现

1. 项目概述

在工业设备维护领域，机械故障诊断一直是个极具挑战性的课题。传统方法往往依赖于简单的时频域分析，难以捕捉复杂工况下设备状态的细微变化。今天要分享的这个项目，我们采用了一种融合高维几何流形学习和最优传输理论的创新方法，通过Python实现了对轴承故障的精准诊断。

这套系统的核心思路是将振动信号映射到高维几何空间，从多个角度提取表征故障特征的几何信息。与常规方法相比，这种几何视角能够更全面地捕捉信号的拓扑结构和动态演化特性。在实际测试中，我们的模型达到了惊人的100%分类准确率，而且特征重要性分析显示，前13维特征就贡献了90%的鉴别力。

2. 核心原理与技术路线

2.1 整体架构设计

系统采用模块化设计，主要包含以下几个关键环节：

1. 数据预处理：对原始振动信号进行去直流、标准化和分段处理
2. 几何特征提取：采用五种不同的几何分析方法并行提取特征
3. 黎曼流形处理：通过相空间重构获取信号的流形结构特征
4. 特征融合与降维：将多源特征标准化后融合并进行降维处理
5. 模型训练与评估：使用随机森林分类器进行训练和性能评估

这种架构的优势在于：

• 多角度特征提取确保不遗漏任何有价值的诊断信息
• 几何方法对噪声和干扰具有更好的鲁棒性
• 模块化设计便于针对不同设备类型进行调整

2.2 关键技术原理

2.2.1 谱几何分析方法

基于图拉普拉斯算子，我们将信号视为图结构数据，通过谱分析提取全局连接特性。具体步骤包括：

1. 构建信号的邻接矩阵，反映各采样点间的相似性
2. 计算归一化拉普拉斯矩阵：L = D^(-1/2)(D-W)D^(-1/2)
3. 求解特征值和特征向量，得到谱特征值、谱熵等指标

这种方法特别适合捕捉信号中的周期性故障特征，比如轴承的周期性冲击。

2.2.2 曲率流方法

曲率流模拟了热扩散过程，可以追踪信号局部几何特性的动态演化：

1. 计算信号的一阶和二阶导数
2. 使用公式 κ = |f''|/(1+f'²)^(3/2) 计算曲率
3. 通过曲率随时间的变化分析信号的动态特性

曲率流对信号的局部突变非常敏感，适合检测早期微弱故障。

2.2.3 李群方法

将信号片段视为李群元素（特别是旋转矩阵），通过李代数分析其群结构特征：

1. 将信号分段映射到SO(n)群
2. 通过指数映射和对数映射在群和代数空间转换
3. 计算交换子等群运算特征

这种方法保留了信号的相位信息，对某些类型的故障特别有效。

3. 详细实现步骤

3.1 数据准备与预处理

python复制class DataLoader:
    def __init__(self, data_dir='data', sampling_freq=12000):
        self.data_dir = data_dir
        self.sampling_freq = sampling_freq
        self.label_map = {
            '98raw.txt': 0,    # 正常状态
            '106raw.txt': 1,   # 滚珠故障
            '131raw.txt': 2,   # 外圈故障
            '119raw.txt': 3    # 内圈故障
        }
    
    def load_data(self, segment_length=2048, overlap=0.5):
        X, y = [], []
        for filename in os.listdir(self.data_dir):
            if filename.endswith('.txt'):
                filepath = os.path.join(self.data_dir, filename)
                label = self.label_map[filename]
                
                data = pd.read_csv(filepath, header=None, names=['vibration'])
                signal = data['vibration'].values
                
                # 预处理
                signal = signal - np.mean(signal)  # 去直流
                signal = (signal - np.mean(signal)) / np.std(signal)  # 标准化
                
                # 分段
                step = int(segment_length * (1 - overlap))
                for i in range(0, len(signal) - segment_length + 1, step):
                    segment = signal[i:i + segment_length]
                    X.append(segment)
                    y.append(label)
        
        return np.array(X), np.array(y)

关键参数选择考量：

• segment_length=2048：足够捕获故障特征又不至于过长
• overlap=0.5：确保不会遗漏任何重要事件
• 标准化处理：消除不同传感器和工况下的量纲差异

3.2 几何特征提取实现

python复制class AdvancedGeometryFeatureExtractor:
    def calculate_curvature(self, signal):
        """计算曲率特征"""
        try:
            dy = savgol_filter(signal, self.window_size, self.polyorder, deriv=1)
            d2y = savgol_filter(signal, self.window_size, self.polyorder, deriv=2)
            curvature = np.abs(d2y) / (1 + dy**2)**1.5
            return np.nan_to_num(curvature)
        except:
            dy = np.gradient(signal)
            d2y = np.gradient(dy)
            curvature = np.abs(d2y) / (1 + dy**2)**1.5
            return np.nan_to_num(curvature)
    
    def spectral_analysis(self, signal):
        """谱几何分析"""
        # 构建相似性矩阵
        dist_matrix = squareform(pdist(signal.reshape(-1,1)))
        W = np.exp(-dist_matrix**2 / (2 * np.median(dist_matrix)**2))
        
        # 计算拉普拉斯矩阵
        D = np.diag(np.sum(W, axis=1))
        L = np.eye(len(signal)) - np.dot(np.dot(np.linalg.inv(D), W), np.linalg.inv(D))
        
        # 特征分解
        eigvals = np.linalg.eigvalsh(L)
        eigvals = np.sort(eigvals)[1:]  # 去掉0特征值
        
        # 提取特征
        features = {
            'spectral_entropy': -np.sum(eigvals * np.log(eigvals + 1e-10)),
            'spectral_energy': np.sum(eigvals**2),
            'spectral_slope': (eigvals[-1] - eigvals[0]) / len(eigvals)
        }
        return features

实际应用中发现几个关键点：

1. 曲率计算时，savgol_filter比直接梯度更稳定
2. 谱分析中，高斯核的带宽选择对结果影响很大
3. 特征值计算需要处理数值稳定性问题

3.3 黎曼流形处理

python复制def riemannian_features(signal, embed_dim=5):
    """黎曼流形特征提取"""
    # 相空间重构
    tau = 10  # 延迟时间
    embedded = np.array([signal[i:i+embed_dim*tau:tau] 
                        for i in range(len(signal)-embed_dim*tau)])
    
    # 计算协方差矩阵
    cov_matrices = np.array([np.cov(sig.T) for sig in embedded])
    
    # 黎曼均值计算
    mean_matrix = np.mean(cov_matrices, axis=0)
    
    # 切空间投影
    tangent_vectors = []
    for C in cov_matrices:
        log_map = logm(np.dot(np.linalg.inv(sqrtm(mean_matrix)), 
                             np.dot(C, np.linalg.inv(sqrtm(mean_matrix)))))
        tangent_vectors.append(log_map[np.triu_indices(embed_dim)])
    
    return np.array(tangent_vectors).flatten()

注意事项：

1. 延迟时间τ的选择很关键，通常用自相关函数确定
2. 嵌入维度要足够大以捕获动态特性，但也不能太大
3. 矩阵对数运算需要处理奇异值问题

4. 特征融合与模型训练

4.1 特征融合策略

python复制def feature_fusion(geo_features, riemann_features):
    """特征融合"""
    # 标准化
    geo_scaled = StandardScaler().fit_transform(geo_features)
    riemann_scaled = StandardScaler().fit_transform(riemann_features)
    
    # PCA降维
    pca_geo = PCA(n_components=0.95)
    geo_reduced = pca_geo.fit_transform(geo_scaled)
    
    pca_riemann = PCA(n_components=0.95)
    riemann_reduced = pca_riemann.fit_transform(riemann_scaled)
    
    # 特征拼接
    fused_features = np.hstack([geo_reduced, riemann_reduced])
    return fused_features

融合过程中的经验：

1. 先分别标准化再融合，避免量纲差异
2. PCA保留95%方差是个不错的起点
3. 特征拼接前要确保维度匹配

4.2 随机森林模型训练

python复制def train_model(X, y):
    """模型训练"""
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.2, stratify=y)
    
    # 参数网格
    param_grid = {
        'n_estimators': [100, 200],
        'max_depth': [10, 20, None],
        'min_samples_split': [2, 5],
        'min_samples_leaf': [1, 2]
    }
    
    # 网格搜索
    grid = GridSearchCV(RandomForestClassifier(),
                       param_grid,
                       cv=5,
                       scoring='accuracy',
                       n_jobs=-1)
    grid.fit(X_train, y_train)
    
    # 评估
    best_model = grid.best_estimator_
    y_pred = best_model.predict(X_test)
    print(classification_report(y_test, y_pred))
    
    return best_model

调参心得：

1. 随机森林对max_depth比较敏感
2. 增加n_estimators总能提升性能，但要考虑计算成本
3. 类别不平衡时需要调整class_weight

5. 结果分析与可视化

5.1 性能评估指标

我们在测试集上获得了以下结果：

这种完美分类结果在实际应用中很少见，可能的原因包括：

1. 数据集本身区分度很好
2. 特征提取非常有效
3. 故障模式比较典型

5.2 特征重要性分析

前13个最重要特征的贡献度：

可视化代码示例：

python复制def plot_feature_importance(model, feature_names):
    """绘制特征重要性"""
    importances = model.feature_importances_
    indices = np.argsort(importances)[::-1]
    
    plt.figure(figsize=(10,6))
    plt.title("Feature Importance")
    plt.bar(range(len(indices)), importances[indices], align='center')
    plt.xticks(range(len(indices)), [feature_names[i] for i in indices], rotation=90)
    plt.tight_layout()
    plt.show()

5.3 故障诊断可视化

我们提供了多种可视化方式帮助理解诊断结果：

1. 原始信号波形与曲率变化对比
2. 特征空间的t-SNE降维投影
3. 混淆矩阵与分类置信度
4. 特征重要性排序

这些可视化不仅验证了模型效果，也为故障机理分析提供了直观依据。

6. 工程实践建议

在实际部署这套系统时，有几个关键点需要注意：

1. 数据采集质量至关重要，采样频率要足够高
2. 不同设备可能需要调整特征提取参数
3. 在线应用时要考虑计算效率
4. 定期用新数据重新训练模型

对于想要复现这个项目的工程师，我建议：

• 先从公开数据集（如CWRU轴承数据）开始
• 逐步实现各个特征提取模块
• 重点理解每种几何特征的实际物理意义
• 可视化中间结果帮助调试

这个项目的完整代码已经整理成模块化结构，便于集成到现有监测系统中。对于特定应用场景，可能需要调整以下参数：

• 信号分段长度
• 各种几何特征的权重
• 分类器的超参数

在实际应用中，这套方法不仅适用于轴承故障诊断，经过适当调整后也可用于齿轮箱、电机等其他旋转机械的故障诊断。关键在于理解设备的故障机理，选择合适的几何表征方法。