首先从振动信号文件中加载数据，进行去除直流分量、标准化和分段预处理，形成样本集并划分训练集和测试集。
然后通过多路并行的几何特征提取方法从信号中提取丰富的几何信息：谱几何方法基于图拉普拉斯算子分析信号的全局连接性和结构特性，提取谱特征值、谱熵和谱斜率等；曲率流方法模拟曲率随时间演化的热扩散过程，捕捉动态几何变化特征；李群方法将信号片段近似为旋转矩阵元素，通过李代数和交换子计算反映群结构的几何特征；最优传输方法使用Wasserstein距离比较信号分布与均匀分布的差异，提取分布几何特征；拓扑数据分析方法计算信号的极值点分布和连通分量，近似持久同调特征。
同时，通过相空间重构将一维时间序列转换为高维相空间中的轨迹矩阵，计算协方差矩阵形成对称正定矩阵流形，利用黎曼图嵌入技术降维后投影到黎曼均值的切空间，得到黎曼流形特征向量。
将几何特征和黎曼特征标准化降维后进行特征融合，形成综合特征表示，使用随机森林分类器训练故障分类模型，最后对新样本进行预测分类，并通过多种可视化手段展示信号特性、特征分布、混淆矩阵、预测置信度和特征重要性，全面评估系统性能后保存训练好的模型。
详细算法步骤：
数据加载与预处理：从指定文件夹读取振动信号文件，根据文件名确定故障类型标签，对原始信号进行去除直流分量和标准化处理，将长信号分割成固定长度的样本段用于后续分析。
信号分段与标准化：设置适当的样本长度和重叠率，将预处理后的信号分割成多个样本，确保每个样本包含足够的信息用于特征提取，同时对样本进行标准化处理消除量纲影响。
创新性几何特征提取：采用多种先进的几何方法从信号中提取丰富的特征，包括基于图拉普拉斯算子的谱几何分析提取全局结构特征，通过曲率流模拟提取动态演化特征，利用李群理论提取旋转结构特征，运用最优传输理论提取分布差异特征，结合拓扑数据分析提取拓扑不变特征。
黎曼流形特征提取：通过相空间重构将一维时间序列转换为高维相空间中的轨迹矩阵，计算协方差矩阵作为对称正定矩阵，构建黎曼流形结构，采用图嵌入技术降维后投影到切空间得到黎曼几何特征。
特征融合与降维：将提取的多类几何特征进行标准化处理，通过主成分分析降维保留主要信息成分，与黎曼特征按维度对齐后进行拼接融合，形成综合的特征表示向量。
分类器训练与优化：使用随机森林算法对融合特征进行训练，通过网格搜索优化超参数设置，利用交叉验证评估模型性能，建立故障类型与特征之间的映射关系模型。
测试集预测与评估：对新样本提取相同的几何特征和黎曼特征，输入训练好的分类模型进行预测，计算准确率、精确率、召回率和F1分数等性能指标，评估模型的泛化能力。
可视化分析与展示：绘制原始信号波形和曲率变化图展示信号特性，通过PCA和t-SNE降维可视化特征分布，绘制混淆矩阵分析分类结果，展示预测置信度分布和特征重要性排序。
模型保存与应用：将训练好的特征提取器和分类模型保存到文件，构建完整的故障诊断系统，支持对新采集的振动信号进行在线或离线的轴承状态识别和故障诊断。
运行结果分析
从特征提取过程来看，几何特征提取模块成功从原始信号中提取了60维的几何特征，经过PCA降维后保留了13个主成分，累计解释方差达到95.2%，表明提取的特征具有很高的信息密度和区分度。
训练过程中，黎曼流形处理模块将5×5的SPD矩阵通过图嵌入降维并投影到切空间，产生了25维的黎曼特征。将13维几何特征与13维黎曼特征融合后，形成了26维的综合特征表示，这些特征充分捕捉了信号在局部和全局几何结构上的细微差异。
随机森林分类器在训练集上达到100%准确率，并在测试集上同样实现完美分类，表明模型没有出现过拟合，具有良好的泛化能力。
特征重要性分析显示，前11个特征贡献了80%的重要性，前13个特征贡献了90%的重要性，说明这些特征具有很高的鉴别力。
整个系统运行稳定，可视化分析清晰展示了信号特性、特征分布和分类结果，为轴承故障诊断提供了一种高效可靠的解决方案。
import numpy as np import pandas as pd from scipy.linalg import logm, expm, eigvals, eig, sqrtm from scipy.spatial.distance import pdist, squareform, cdist from scipy.signal import savgol_filter, welch from scipy.sparse import csgraph, diags, lil_matrix from scipy.sparse.linalg import eigs, eigsh from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import classification_report, accuracy_score, confusion_matrix, f1_score from sklearn.decomposition import PCA from sklearn.manifold import TSNE from sklearn.neighbors import NearestNeighbors import matplotlib.pyplot as plt import matplotlib import warnings warnings.filterwarnings('ignore') from tqdm import tqdm import os import seaborn as sns from scipy import integrate, interpolate, optimize from scipy.spatial import ConvexHull, Delaunay import networkx as nx from ot import emd2, sinkhorn2 # 设置样式 matplotlib.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans'] matplotlib.rcParams['axes.unicode_minus'] = False plt.style.use('seaborn-v0_8-darkgrid') # ============================== # 数据加载模块 # ============================== 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 # 内圈故障 } self.label_names = { 0: 'Normal', 1: 'Ball Fault', 2: 'Outer Race Fault', 3: 'Inner Race Fault' } self.colors = ['#2E86AB', '#A23B72', '#F18F01', '#C73E1D'] 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) X = np.array(X) y = np.array(y) print(f"Data loaded successfully") print(f"Total samples: {len(X)}") print(f"Class distribution: {np.bincount(y)}") return X, y # ============================== # 创新性几何特征提取模块 # ============================== class AdvancedGeometryFeatureExtractor: """创新性几何特征提取器""" def __init__(self, n_components=0.95, embedding_dim=3, window_size=51, polyorder=3): self.n_components = n_components self.embedding_dim = embedding_dim self.window_size = window_size self.polyorder = polyorder self.pca = None self.scaler = StandardScaler() # ==================== 基本微分几何方法 ==================== def calculate_curvature(self, signal): """计算曲率""" if len(signal) < self.window_size: window_size = min(self.window_size, len(signal) - 1) if len(signal) > 1 else 1 polyorder = min(self.polyorder, window_size - 1) else: window_size = self.window_size polyorder = self.polyorder if window_size % 2 == 0: window_size += 1 if window_size <= polyorder: window_size = polyorder + 1 if polyorder % 2 == 1 else polyorder + 2 try: dy = savgol_filter(signal, window_size, polyorder, deriv=1) d2y = savgol_filter(signal, window_size, polyorder, deriv=2) curvature = np.abs(d2y) / (1 + dy**2)**1.5 curvature = np.nan_to_num(curvature, nan=0.0, posinf=0.0, neginf=0.0) return curvature except: try: dy = np.gradient(signal) d2y = np.gradient(dy) curvature = np.abs(d2y) / (1 + dy**2)**1.5 curvature = np.nan_to_num(curvature, nan=0.0, posinf=0.0, neginf=0.0) return curvature except: return np.zeros_like(signal)参考文章：
基于高维几何流形学习和最优传输理论融合的机械故障诊断方法（Python） - 哥廷根数学学派的文章
https://zhuanlan.zhihu.com/p/1986810770059716342
工学博士，担任《Mechanical System and Signal Processing》审稿专家，担任
《中国电机工程学报》优秀审稿专家，《控制与决策》，《系统工程与电子技术》，《电力系统保护与控制》，《宇航学报》等EI期刊审稿专家。
擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
