格拉姆角场与轴承故障诊断：从时序信号到图像识别的数据预处理实战

1. 项目概述：从代码到数据，理解故障诊断的基石

拿到一份名为“格拉姆角场东南大学轴承故障诊断代码解读——数据集解读”的代码，很多朋友可能会直接一头扎进模型构建和训练的部分，急切地想看到诊断准确率。但根据我多年的工业数据分析经验，这恰恰是新手最容易踩坑的地方。一个故障诊断项目的成败，在模型跑起来之前，就已经被数据决定了七八成。这份代码的标题将“数据集解读”放在后半部分，而我认为，它应该是我们打开任何类似项目时，第一个、也是最需要花时间吃透的环节。

格拉姆角场（Gramian Angular Field, GAF）是一种将一维时间序列转换为二维图像矩阵的编码方法，它通过保留时间序列的绝对时序关系和数值信息，为后续使用成熟的图像分类模型（如CNN）处理振动信号铺平了道路。而“东南大学轴承数据集”则是在机械故障诊断领域一个非常经典且公开的基准数据集。这个项目本质上，就是利用GAF技术，将东南大学轴承的振动时序信号“翻译”成图像，然后用计算机视觉的方法来识别轴承的健康状态和故障类型。

所以，在动手修改任何一行模型代码之前，我们必须彻底搞清楚：我们喂给模型的是什么“粮食”？这些“粮食”是怎么从原始的振动信号加工而来的？数据里有没有“杂质”或“偏见”？只有把数据集这第一道关把好了，后续的模型调优、结果分析才有意义。否则，很可能出现模型在训练集上表现完美，一到实际场景就“翻车”的情况。接下来，我就带大家深入这个项目的“后厨”，看看数据是如何被准备和处理的。

2. 核心数据集：东南大学轴承数据深度解析

在解读任何相关代码前，我们必须先独立于代码，理解数据本身的来源、结构和物理意义。这是避免被代码实现带偏、形成自己判断力的关键。

2.1 数据来源与采集背景

东南大学（SEU）的轴承数据集是在实验室环境下，通过转子实验台采集的。实验台通常包含电机、转轴、支撑轴承、加载装置等部分，通过在健康轴承和预设故障的轴承上安装振动加速度传感器来收集数据。故障类型通常包括内圈故障、外圈故障、滚动体故障，并且每种故障会有不同尺寸（如0.007英寸，0.014英寸，0.021英寸）的模拟损伤。

这个数据集之所以经典，是因为它：

1. 工况相对可控：负载、转速通常是固定的，这减少了变量，便于初学者聚焦于故障特征本身。
2. 故障模式典型：涵盖了旋转机械中最常见的几种轴承故障类型。
3. 数据格式规整：通常以.mat（MATLAB数据文件）或文本文件形式提供，采样频率、数据长度标注清晰。

注意：实验室数据与现场数据存在“鸿沟”。实验室数据信噪比高，故障特征明显；而现场数据受背景噪声、工况波动、多源耦合振动影响巨大。因此，在实验室数据集上表现优异的模型，直接部署到工厂可能需要大幅调整。

2.2 数据结构与文件组织

通常，下载到的SEU数据集文件夹结构如下：

SEU_Bearing_Dataset/ ├── 正常/ │ ├── normal_1.mat │ ├── normal_2.mat │ └── ... ├── 内圈故障/ │ ├── IR007_1.mat │ ├── IR007_2.mat │ ├── IR014_1.mat │ └── ... ├── 外圈故障/ │ └── ... └── 滚动体故障/ └── ...

每个.mat文件里，通常存储着一个或多个通道的振动加速度时序数据。关键参数需要从数据说明或代码中提取：

• 采样频率（Fs）：如12kHz， 24kHz。这决定了信号能捕获的最高频率（奈奎斯特频率为Fs/2）。
• 数据长度：每个文件可能包含几十万甚至上百万个数据点，代表一段连续采样的振动信号。
• 转速与负载：这些信息对于理解故障特征频率至关重要，但数据集有时未必直接附带，需要从实验描述中查找。

2.3 故障的物理特征与在信号中的体现

轴承的故障会在振动信号中产生周期性冲击，其频率由故障类型和几何参数决定，称为故障特征频率。

• 内圈故障频率（BPFI）：与转频相关，频率较高，且由于载荷方向变化，振幅会有调制现象。
• 外圈故障频率（BPFO）：频率相对固定，振幅稳定。
• 滚动体故障频率（BSF）：频率通常低于内外圈故障。
• 保持架故障频率（FTF）：频率最低。

在时域波形上，健康信号相对平稳，而故障信号会出现明显的、周期性的冲击脉冲。在频域（通过傅里叶变换），我们可以在相应的故障特征频率及其倍频处看到突出的谱线。

理解这些物理背景至关重要，因为GAF将一维信号转为图像后，这些时域和频域的特征会以某种纹理、形状或亮度的模式体现在图像中。我们的目标，就是让CNN模型学会识别这些与特定故障对应的图像模式。

3. 格拉姆角场原理与代码实现拆解

理解了“原材料”（原始振动数据）后，我们来看“烹饪方法”——格拉姆角场。代码中实现GAF的部分是核心，我们需要明白每一步的数学意义和工程考量。

3.1 GAF转换的核心步骤

GAF主要分为两种：格拉姆角和场（GASF）和格拉姆角差场（GADF）。项目代码中通常使用其中一种或两者结合。其转换流程可分解为以下四步，我结合代码中可能出现的函数进行解释：

第一步：数据归一化将原始振动信号X = [x1, x2, ..., xn] 缩放到区间[-1, 1]或[0, 1]。

# 常见代码片段示例 from sklearn.preprocessing import MinMaxScaler scaler = MinMaxScaler(feature_range=(-1, 1)) X_normalized = scaler.fit_transform(X.reshape(-1, 1)).flatten()

为什么必须归一化？因为GAF基于角度计算，而归一化到[-1,1]区间后，数据点可以映射到单位圆上的余弦值。如果数据量纲不统一（例如不同通道、不同实验的数据），绝对值大的信号会主导角度计算，导致信息失真。

第二步：将数值转换为角度通过反余弦函数，将归一化后的值映射为角度（弧度制）。

import numpy as np phi = np.arccos(X_normalized) # 此时 phi 在 [0, pi] 区间内

这一步是GAF的精华。每个数据点不再是一个孤立的振幅值，而是单位圆上的一个点，由其与横轴的夹角φ来表征。时序信息被巧妙地编码进了这个角度序列中。

第三步：计算格拉姆矩阵（核心）这是生成二维图像的关键。格拉姆矩阵的元素由每两个点之间的三角和差关系构成。

• 格拉姆角和场（GASF）：计算角度之和的余弦。它更侧重于捕捉信号之间的“和”关系，反映的是整体相关性。

  # 伪代码逻辑 GASF = np.cos(phi_i + phi_j) # 其中 i, j 遍历所有数据点 # 实际代码利用向量化操作，避免低效循环

• 格拉姆角差场（GADF）：计算角度之差的正弦。它更侧重于捕捉信号之间的“差”或相对变化关系。

  GADF = np.sin(phi_i - phi_j)

生成的GASF或GADF矩阵是一个n x n的对称矩阵（对于GASF）或反对称矩阵（对于GADF），其中n是输入时序片段的长度。这个矩阵就是我们要的“图像”。

第四步：图像化与裁剪生成的n x n矩阵可能很大（例如，1000x1000）。直接作为CNN输入可能计算量过大。因此，代码中通常会有以下操作：

1. 降采样：在计算GAF前，先对长时序信号进行切片或降采样，使n控制在一个合理大小（如224，适配ImageNet预训练模型）。
2. 图像缩放：生成GAF矩阵后，使用cv2.resize或PIL.Image.resize将其缩放到统一尺寸（如224x224）。
3. 伪彩色映射：GAF矩阵是单通道的（每个像素一个值）。为了适配通常输入为3通道的CNN（如ResNet），需要将其转换为“伪彩色”图像。常见方法是使用cv2.applyColorMap（如COLORMAP_JET）或简单地将同一矩阵复制到三个通道。

3.2 代码中的关键参数与选择

在解读data_preprocessing.py或generate_gaf_images.py这类文件时，要重点关注以下参数：

• 切片长度（segment_length）：如1024个点。这决定了每张“图像”代表多长时间的振动信号。太短可能包含不了一个完整的故障冲击周期，太长则图像分辨率过高且可能混合多种状态。经验上，这个长度应能覆盖至少2-3个故障特征周期。例如，转速为1800 rpm（30 Hz），BPFI约为5.4倍转频即162 Hz，周期约6.2毫秒。若采样频率为12 kHz，则一个周期约74个点。因此，切片长度取512或1024点是合理的。
• 重叠率（overlap_rate）：如0.5（50%）。滑动窗口切分信号时，相邻切片之间的重叠比例。提高重叠率可以增加生成图像的数量，缓解数据量不足的问题，但也会引入更强的样本相关性，可能影响模型泛化能力评估的准确性。
• GAF类型选择：只用GASF，还是GADF，或者将两者合并为双通道图像？不同的故障特征在不同场中的表现可能不同。我个人的经验是，对于轴承的周期性冲击故障，GASF往往能更好地保留冲击的时序相关性，效果更稳定。可以尝试融合，但会增加模型输入通道和计算量。
• 缩放方法与插值算法：将GAF矩阵缩放到目标尺寸（如224x224）时，cv2.INTER_LINEAR（双线性插值）是常用选择。应避免使用INTER_NEAREST（最近邻插值），因为它可能在图像中引入块状伪影，破坏连续的特征模式。

4. 数据预处理流程全链路实操

现在，我们把数据集和GAF原理串联起来，看一个完整的、可复现的数据预处理流水线应该如何构建。这是项目代码的核心骨架。

4.1 步骤一：原始数据加载与探查

在写任何处理代码之前，先用Jupyter Notebook或脚本进行数据探查。

import scipy.io as sio import numpy as np import matplotlib.pyplot as plt # 1. 加载一个.mat文件示例 data_dict = sio.loadmat(‘path/to/IR007_1.mat‘) # 打印所有键，查看数据结构 print(data_dict.keys()) # 通常振动数据在 ‘data‘, ‘vibration‘, ‘X‘ 等键下 vibration_signal = data_dict[‘X‘].flatten() # 假设键名为‘X‘，并转换为一维数组 # 2. 绘制时域波形 plt.figure(figsize=(12, 4)) plt.plot(vibration_signal[:5000]) # 只看前5000个点 plt.title(‘Raw Vibration Signal (Time Domain)‘) plt.xlabel(‘Sample Points‘) plt.ylabel(‘Amplitude‘) plt.grid(True) plt.show() # 3. 计算并绘制频谱（快速傅里叶变换） from scipy.fft import fft, fftfreq Fs = 12000 # 假设采样频率为12kHz N = len(vibration_signal) yf = fft(vibration_signal) xf = fftfreq(N, 1/Fs)[:N//2] # 取正频率部分 plt.figure(figsize=(12,4)) plt.plot(xf, 2.0/N * np.abs(yf[0:N//2])) plt.title(‘Frequency Spectrum‘) plt.xlabel(‘Frequency (Hz)‘) plt.ylabel(‘Magnitude‘) plt.grid(True) plt.xlim([0, Fs/2]) # 显示到奈奎斯特频率 plt.show()

这个探查步骤能帮你确认信号质量，观察是否有明显的故障冲击，并验证采样频率。

4.2 步骤二：数据切片与标签生成

这是为后续GAF转换准备输入片段和对应标签。

def segment_signal(signal, label, segment_length, overlap_ratio): """ 将一维信号切分成固定长度的片段，并分配标签。 参数: signal: 一维振动信号数组。 label: 该信号对应的整数型标签（如0:正常，1:内圈故障...）。 segment_length: 每个片段的长度。 overlap_ratio: 重叠率，0-1之间。 返回: segments: 片段列表，形状为 (num_segments, segment_length)。 labels: 标签列表，形状为 (num_segments,)。 """ segments = [] labels = [] step = int(segment_length * (1 - overlap_ratio)) if step == 0: step = 1 # 避免死循环 num_segments = (len(signal) - segment_length) // step + 1 for i in range(num_segments): start = i * step end = start + segment_length segment = signal[start:end] # 可选：这里可以添加片段能量检查，过滤掉能量过低的无效片段 segments.append(segment) labels.append(label) return np.array(segments), np.array(labels) # 遍历所有数据文件夹，收集所有片段和标签 all_segments = [] all_labels = [] class_folders = {‘normal‘: 0, ‘IR‘: 1, ‘OR‘: 2, ‘Ball‘: 3} # 示例映射 for class_name, label_id in class_folders.items(): folder_path = os.path.join(‘dataset‘, class_name) for file_name in os.listdir(folder_path): if file_name.endswith(‘.mat‘): file_path = os.path.join(folder_path, file_name) signal = load_signal_from_mat(file_path) # 自定义加载函数 segments, labels = segment_signal(signal, label_id, segment_length=1024, overlap_ratio=0.5) all_segments.extend(segments) all_labels.extend(labels) # 转换为NumPy数组 all_segments = np.array(all_segments) all_labels = np.array(all_labels) print(f“总片段数： {all_segments.shape[0]}， 片段长度： {all_segments.shape[1]}， 类别数： {len(np.unique(all_labels))}“)

4.3 步骤三：GAF图像批量生成

将上一步得到的所有信号片段批量转换为GAF图像。

from pyts.image import GramianAngularField # 可以使用pyts库，也可以自己实现 import cv2 # 方法1：使用pyts库（推荐，稳定且高效） def generate_gaf_images_pyts(segments, image_size=224, method=‘summation‘): """ 使用pyts库生成GAF图像。 method: ‘summation‘ (GASF) 或 ‘difference‘ (GADF) """ gasf = GramianAngularField(image_size=image_size, method=method) # pyts要求输入形状为 (n_samples, n_timestamps) images_gasf = gasf.fit_transform(segments) # 输出形状 (n_samples, image_size, image_size) # 将值域从[-1,1]或[0,1]映射到[0, 255]的uint8，并应用伪彩色 images_uint8 = ((images_gasf + 1) * 127.5).astype(np.uint8) # 假设值域为[-1,1] colored_images = [] for img in images_uint8: colored = cv2.applyColorMap(img, cv2.COLORMAP_JET) colored_images.append(colored) return np.array(colored_images) # 形状 (n_samples, image_size, image_size, 3) # 方法2：手动实现（更灵活，便于理解原理） def gramian_angular_field(series, method=‘summation‘): """手动计算单一样本的GAF矩阵""" # 归一化 min_val, max_val = series.min(), series.max() scaled_series = (2 * (series - min_val) / (max_val - min_val)) - 1 # 归一化到[-1,1] scaled_series = np.clip(scaled_series, -1, 1) # 防止反余弦计算溢出 # 计算角度 phi = np.arccos(scaled_series) # 计算格拉姆矩阵 if method == ‘summation‘: # GASF = cos(φ_i + φ_j) cos_sum = np.cos(np.add.outer(phi, phi)) return cos_sum elif method == ‘difference‘: # GADF = sin(φ_i - φ_j) sin_diff = np.sin(np.subtract.outer(phi, -phi)) # 注意符号处理 return sin_diff # 批量处理 image_size = 224 gaf_images = [] for segment in all_segments[:100]: # 示例：先处理100个 gaf_matrix = gramian_angular_field(segment) # 缩放 resized_matrix = cv2.resize(gaf_matrix, (image_size, image_size), interpolation=cv2.INTER_LINEAR) # 伪彩色和归一化到[0,255] normalized = ((resized_matrix + 1) * 127.5).astype(np.uint8) colored = cv2.applyColorMap(normalized, cv2.COLORMAP_VIRIDIS) gaf_images.append(colored) gaf_images = np.array(gaf_images)

4.4 步骤四：数据集划分与保存

将生成的图像数据集划分为训练集、验证集和测试集，并保存为文件（如TFRecord或直接NumPy数组+标签），方便模型加载。

from sklearn.model_selection import train_test_split import pickle # 划分数据集（先划分索引，避免数据混乱） indices = np.arange(len(gaf_images)) X_train_idx, X_temp_idx, y_train_idx, y_temp_idx = train_test_split( indices, all_labels[:len(gaf_images)], test_size=0.3, random_state=42, stratify=all_labels[:len(gaf_images)] ) X_val_idx, X_test_idx, y_val_idx, y_test_idx = train_test_split( X_temp_idx, y_temp_idx, test_size=0.5, random_state=42, stratify=y_temp_idx ) # 根据索引获取数据 X_train, y_train = gaf_images[X_train_idx], all_labels[X_train_idx] X_val, y_val = gaf_images[X_val_idx], all_labels[X_val_idx] X_test, y_test = gaf_images[X_test_idx], all_labels[X_test_idx] print(f“训练集： {X_train.shape}， 验证集： {X_val.shape}， 测试集： {X_test.shape}“) # 保存数据集 save_dict = { ‘X_train‘: X_train, ‘y_train‘: y_train, ‘X_val‘: X_val, ‘y_val‘: y_val, ‘X_test‘: X_test, ‘y_test‘: y_test, ‘label_names‘: {0: ‘正常‘, 1: ‘内圈故障‘, 2: ‘外圈故障‘, 3: ‘滚动体故障‘} } with open(‘seu_bearing_gaf_dataset.pkl‘, ‘wb‘) as f: pickle.dump(save_dict, f) print(“数据集已保存为 ‘seu_bearing_gaf_dataset.pkl‘“)

5. 关键注意事项与避坑指南

在实际操作中，有几个细节如果不注意，很容易导致模型效果不佳或结论错误。

5.1 数据泄露问题

这是时序数据划分中最常见的坑。绝对不能用随机打乱后再划分的方法！

• 错误做法：将all_segments和all_labels用sklearn.model_selection.train_test_split直接随机划分。因为重叠切片的存在，同一个原始样本的不同切片可能被分到了训练集和测试集，导致模型通过“记忆邻居”就能做出正确预测，严重高估泛化能力。
• 正确做法：按“样本源文件”划分。即，将所有.mat文件列表打乱，然后按比例分配给训练集、验证集和测试集。确保来自同一个原始数据文件的所有切片，都只出现在同一个集合中。这样才能模拟现实场景：模型用一批机器历史数据训练，去预测另一批机器未来的状态。

5.2 类别不平衡处理

轴承故障数据中，正常状态的数据往往远多于各种故障状态的数据（因为故障是偶发事件）。直接训练会导致模型偏向于预测“正常”类。

• 对策：在数据加载或训练时进行处理。
  1. 过采样：对少数类样本进行复制或使用SMOTE（需谨慎，对于图像数据，简单的复制可能导致过拟合）。
  2. 欠采样：随机丢弃一部分多数类样本，可能损失有用信息。
  3. 类别权重：在损失函数中为少数类赋予更高的权重。这是最常用且有效的方法。在PyTorch或TensorFlow中，可以方便地设置class_weight参数。

  # 以sklearn为例计算类别权重 from sklearn.utils.class_weight import compute_class_weight class_weights = compute_class_weight(‘balanced‘, classes=np.unique(y_train), y=y_train) # 在训练时，将class_weights传递给损失函数

5.3 GAF图像的可视化与检查

生成GAF图像后，一定要抽样可视化，检查转换是否合理。

import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 4, figsize=(16, 8)) for i in range(4): # 展示4个类别，每个类别2个样本 # 选取某个类别的样本索引 class_idx = np.where(y_train == i)[0] sample_idx = class_idx[0] axes[0, i].imshow(X_train[sample_idx]) axes[0, i].set_title(f‘{label_names[i]} - Sample 1‘) axes[0, i].axis(‘off‘) sample_idx = class_idx[1] axes[1, i].imshow(X_train[sample_idx]) axes[1, i].set_title(f‘{label_names[i]} - Sample 2‘) axes[1, i].axis(‘off‘) plt.tight_layout() plt.show()

你需要观察：不同类别的图像在纹理、颜色分布上是否有肉眼可辨的差异？同一类别的不同样本是否具有相似的模式？如果看起来都是杂乱无章的噪声，可能需要检查GAF转换的参数（如归一化是否正确）或考虑原始信号是否信噪比过低。

5.4 计算资源与效率优化

GAF计算是计算密集型操作，尤其是当数据量大、切片长度长时。

• 向量化操作：务必使用NumPy的广播和向量化函数（如np.add.outer,np.subtract.outer）来替代Python循环，速度可提升成百上千倍。
• 分批处理与存储：不要试图一次性将所有数据转为GAF图像并加载进内存。应该设计一个生成器（Generator）或使用tf.data.Dataset/torch.utils.data.DataLoader的预处理管道，在需要时动态生成批次数据，或者预处理后存储为.tfrecord或.h5格式，便于流式读取。
• 并行处理：利用multiprocessing库或joblib对多个数据文件或片段进行并行GAF转换。

6. 进阶思考：从SEU数据集到工业实际应用

基于SEU数据集和GAF的方法在实验室环境下可以取得很高的准确率（>99%），但这离真正的工业应用还有距离。在解读代码、复现结果之后，我们更应该思考以下几个问题，这也是项目价值能否延伸的关键：

6.1 模型学到了什么？——可解释性分析

当你的CNN模型在测试集上达到高精度后，不要满足于此。我们需要知道模型是根据图像的哪些部分做出判断的，这有助于验证模型是否真的学到了有物理意义的故障特征，而不是一些无关的伪影。

• 梯度加权类激活映射（Grad-CAM）：这是最常用的可视化方法。它可以生成一个热力图，叠加在原始GAF图像上，显示哪些区域对模型决策的贡献最大。
• 分析：观察正常样本和故障样本的热力图差异。对于故障样本，高亮区域是否集中在图像的某些特定对角线或区块？这些区域可能对应着原始振动信号中周期性冲击发生的时间点。如果热力图总是集中在图像边缘或无关区域，那模型的决策依据可能是可疑的。

6.2 如何处理变工况与噪声？——数据增强与域适应

SEU数据集的工况是固定的。但现实中，设备的转速、负载会变化，环境噪声也更大。

• 数据增强：在生成GAF图像后或之前，可以引入增强策略。
  • 时域增强：对原始信号添加随机缩放、轻微抖动、添加高斯白噪声。
  • 图像域增强：对GAF图像使用标准的图像增强，如随机旋转（小角度）、水平/垂直翻转（需谨慎，GAF图像具有对称性，翻转可能改变物理意义）、亮度对比度微调。
• 域适应（Domain Adaptation）：如果目标工况与SEU数据集差异巨大，可以考虑使用迁移学习或域自适应算法（如DANN），利用SEU的带标签数据和目标域的无标签数据，让模型学会工况不变的特征。

6.3 除了GAF，还有哪些编码方式？

GAF不是唯一将时序信号转为图像的方法。了解其他方法有助于你在不同场景下做出选择。

• 马尔可夫变迁场（Markov Transition Field, MTF）：将时序数据离散化后，计算状态间转移概率的矩阵。对刻画信号的状态变迁规律有效。
• 递归图（Recurrence Plots, RP）：展示信号在相空间中哪些时刻回到相似状态。对非线性、非平稳信号分析有优势。
• 连续小波变换时频谱图（CWT Scalogram）：这是更经典的方法，能同时提供时间和频率信息。对于轴承故障，冲击成分会在特定频带产生“条纹”，非常直观。

如何选择？一个实用的建议是：从时频谱图（CWT）开始。因为它最符合振动分析工程师的看图习惯，特征物理意义明确。如果追求极致的端到端自动化，再尝试GAF、MTF等编码方法，并与CWT进行效果对比。

6.4 部署考量：轻量化与实时性

实验室代码往往不考虑推理速度。但工业部署要求模型轻量化、推理快。

• 模型压缩：考虑使用MobileNet、EfficientNet等轻量级CNN backbone，或者对训练好的模型进行剪枝、量化。
• 推理流水线优化：将GAF转换和模型推理集成到C++或嵌入式环境中，使用ONNX Runtime、TensorRT等推理引擎进行加速。
• 边缘部署：对于实时监测，可以考虑在设备边缘（如工控机、带算力的网关）直接运行轻量级模型，只将报警结果或特征上传到云端，而非原始振动数据。

解读“格拉姆角场东南大学轴承故障诊断代码”的最终目的，绝不是为了在SEU数据集上刷出一个漂亮的数字。而是通过这个标准的“实验室样板”，掌握“数据理解 -> 特征工程（GAF）-> 模型构建 -> 评估分析”的完整方法论链条。当你吃透了数据集的每一个字节，理解了GAF矩阵中每一个像素的来源与含义，你才真正拥有了将这套方法迁移到其他设备、其他故障诊断场景的能力。数据是地基，特征工程是骨架，模型只是血肉。把地基打牢，骨架搭正，项目才能立得住，走得远。