首先从原始振动信号中加载数据并进行预处理，包括去直流、标准化和信号分段。然后，采用微分几何方法从每个信号段中提取几何特征，包括曲率、挠率、测地距离、Frenet标架、黎曼曲率张量和流形拓扑特征等，这些特征描述了信号在相空间中的局部和全局几何结构。
同时通过相空间重构将一维时间序列转换为高维相空间中的轨迹矩阵，并计算其协方差矩阵作为对称正定矩阵，这些矩阵构成黎曼流形上的点。接着，使用黎曼图嵌入技术将这些高维SPD矩阵降维，并通过对数映射将它们投影到黎曼均值的切空间中，得到黎曼流形特征向量。
随后将微分几何特征与黎曼特征按维度对齐并融合，形成综合的特征表示。最后，使用随机森林分类器对融合特征进行训练和测试，实现对正常状态、滚珠故障、外圈故障和内圈故障四类轴承状态的精确分类，并通过多种可视化方法展示特征分布、分类结果和特征重要性。
整个算法从振动信号数据开始，首先进行数据预处理，包括去除直流分量、标准化和将长信号分割为固定长度的样本段。然后分两路并行提取特征：第一路采用微分几何方法，计算每个信号段的曲率特征来表征局部弯曲程度，通过相空间重构分析信号在高维空间中的几何结构，提取流形特征描述信号的拓扑性质，计算测地距离反映信号在嵌入空间中的内在距离，构建Frenet标架分析信号的运动几何特性，并从协方差矩阵中提取黎曼曲率张量特征；第二路采用黎曼流形方法，通过相空间重构构建轨迹矩阵并计算协方差矩阵作为对称正定矩阵，这些矩阵形成黎曼流形，接着使用图嵌入技术将这些高维流形点降维，再通过黎曼均值和对数映射将流形点投影到切空间中得到向量表示。然后将两路特征分别进行标准化和降维处理，按维度对齐后进行特征融合。使用融合后的特征训练随机森林分类器，学习不同故障类型的判别模式。在测试阶段，对新样本提取同样的特征并进行分类预测。最后通过混淆矩阵、特征分布可视化、预测置信度分析和特征重要性评估等方式全面展示算法性能，并将训练好的模型保存以备后续使用，整个流程实现了从原始振动信号到故障类别判别的端到端智能诊断。
import numpy as np import pandas as pd from scipy.linalg import logm, expm, eigvals, eig, sqrtm from scipy.spatial.distance import pdist, squareform from scipy.signal import savgol_filter 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 from scipy.spatial import ConvexHull # 设置样式 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 DifferentialGeometryFeatureExtractor: """微分几何特征提取器""" 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): """计算曲率""" # 使用Savitzky-Golay滤波器计算导数 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 # 处理NaN和无穷大值 curvature = np.nan_to_num(curvature, nan=0.0, posinf=0.0, neginf=0.0) return curvature except Exception as e: # 如果滤波失败，使用简单差分 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/1986781299621524163
工学博士，担任《Mechanical System and Signal Processing》审稿专家，担任
《中国电机工程学报》优秀审稿专家，《控制与决策》，《系统工程与电子技术》，《电力系统保护与控制》，《宇航学报》等EI期刊审稿专家。
擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。
