高速列车轴承故障诊断：多尺度加权融合模型与迁移学习可解释性分析（附代码）

时速350公里的复兴号高铁，轮对轴承以每秒近30转的速度高速旋转，承受着数吨的载荷。轴承上任何微小的损伤，若没能被及时捕捉，都可能像滚雪球一样演变成严重事故。这不仅是机械问题，更是一个典型的数据科学挑战——我们拥有海量的实验室台架数据（源域），却极度缺乏真实运行中的故障数据（目标域）。如何将实验室里的“知识”有效迁移到实际列车上？这正是本文要讲述的故事。

本文改编自实际完成的咨询项目，核心思路是：通过多尺度加权GOPAE（全局序模式注意力熵）与迁移学习技术，解决高速列车轴承故障诊断中“数据丰富但场景不足”的痛点。我们构建了一套从数据筛选、特征提取、源域建模到跨域迁移诊断的完整解决方案，并特别关注模型的可解释性，让每一个诊断结果都能追溯到物理机理。

引言

想象这样一个场景：一列高速列车正载着上千名旅客穿越华北平原，轴承的微小损伤如果不能被及时察觉，后果不堪设想。实验室里我们能模拟各种故障，但实际运营中故障样本极其稀缺——这正是工业智能诊断的核心难题。

1 背景与问题

高速列车已成为中国客运的主流方式，轴承作为走行部的关键旋转部件，其安全状态直接关系行车安全。实际运营中，轴承故障数据极其稀缺，而实验室台架数据丰富但工况不同。如何利用台架数据（源域）辅助诊断实际列车轴承（目标域），是工业智能诊断的核心难题。

图：迁移学习分析图

2 数据筛选与特征提取

2.1 总体思路

整个建模过程遵循“数据分析与特征提取 → 源域模型构建 → 跨域迁移诊断 → 模型可解释性分析”的逻辑链条，思路流程如下图：

图：问题总体思路流程图

2.2 数据筛选策略

基于故障机理，从161个源域文件中筛选出12个代表性文件，涵盖正常(N)、外圈故障(OR)、内圈故障(IR)、滚动体故障(B)四种状态，优先选择驱动端(DE)数据，并重采样至32kHz统一频率。

下图展示了筛选后的数据分布：

图：源域数据筛选与重采样统计分析

图：源域数据筛选与预处理

2.3 特征提取方法

我们构建了“时域-频域-动态复杂性”三维特征体系：

• 时域特征：均值、标准差、RMS、峰值、峭度系数、峰度、偏度、形状因子等，共9维。
• 频域特征：谱心、谱方差、低频比(<2kHz)、中频比(2-8kHz)、高频比(>8kHz)，经转频归一化，共5维。
• 动态特征：多尺度加权GOPAE（全局序模式注意力熵），捕捉故障冲击的周期性规律，共5维（尺度1~5），后加权融合为1维。

GOPAE的核心思想是将信号转换为序模式，计算每种模式出现间隔的熵，熵值低表示信号规律性强（如故障冲击），熵值高表示随机性强（如正常状态）。

图：GOPAE模型分析图

图：多尺度加权GOPAE模型分析图

图：GOPAE核心函数实现图

图：多尺度权重计算与特征融合图

图：GOPAE多尺度权重分布图

图：频谱分析图

2.4 特征分析与选择

通过F检验评估各特征的判别能力，筛选出F值最高的6维特征：加权GOPAE、峭度系数、峰值因子、高频能量比、RMS、谱心。这些特征在四类状态间表现出显著差异。

图：源域各故障类型代表性信号波形时域特征分析

图：源域各故障类型功率谱密度分析

图：各故障类型多尺度GOPAE特征对比

图：各故障类型序模式分布特征

图：关键时域特征分布箱线图

图：频域特征分布特性分析

图：多尺度GOPAE特征雷达图

图：特征相关性矩阵分析图

图：特征分布箱线图

3 源域诊断模型构建

3.1 数据预处理

图：数据预处理流程图

3.2 模型设计

采用集成学习策略，融合三种基础模型：

• SVM (支持向量机)：使用RBF核，处理非线性特征，捕捉故障细微差异。算法流程如图：

图：SVM算法流程图

• 随机森林 (Random Forest)：通过多棵决策树投票，降低过拟合，擅长处理统计特征。算法流程如图：

图：随机森林算法流程图

• 梯度提升 (Gradient Boosting)：迭代优化，聚焦难分样本，提升对相似故障的区分能力。

通过软投票机制融合三模型概率，最终输出类别。

图：集成策略流程图

图：特征协同流程图

3.3 训练与评估

将筛选后的6维特征输入集成模型，采用5折分层交叉验证。模型在测试集上准确率达0.98，宏平均F1分数0.97，表明模型具有良好的泛化能力。

混淆矩阵显示，各故障类别均被准确识别，无明显误判：

图：混淆矩阵图

图：多分类分析图

特征重要性分析表明，加权GOPAE和高频比贡献最大，与故障机理高度吻合。

4 跨域迁移诊断

4.1 目标域数据适配

目标域为16个无标签列车轴承振动文件（编号A-P），采样频率32kHz，转速约600rpm。将其下采样至12kHz（与源域一致），采用相同窗口（4096点，重叠率0.5）提取特征，并复用源域的特征选择与标准化参数。域偏移度量MMD（最大均值差异）从0.28降至0.06，表明特征分布已很好对齐。

4.2 迁移诊断策略

采用“特征迁移+模型迁移”混合策略：

• 复用源域集成模型的结构与权重。
• 对目标域每个文件的多个窗口样本进行预测，通过多数投票确定文件最终类别。
• 计算平均预测概率作为置信度，分为四级（高≥0.9，中0.75~0.9，低0.6~0.75，不可靠<0.6）。

图：迁移诊断模型流程图

4.3 诊断结果

16个文件的分类结果：外圈故障(OR) 4个，内圈故障(IR) 5个，滚动体故障(B) 4个，正常(N) 3个。平均置信度0.79，标准差0.11，无极端低置信结果，满足工业应用要求。

5 模型可解释性分析

5.1 事前可解释性：特征与机理绑定

每个特征都有明确的物理意义：

• RMS：振动能量，故障时升高。
• 峭度：冲击程度，故障时显著增大。
• 高频比：故障激发高频成分，外圈故障最明显。
• GOPAE：序模式间隔规律性，故障时熵值降低。

这些特征与轴承故障机理（BPFO、BPFI等特征频率）直接对应，使模型决策有据可依。

5.2 过程可解释性：域分布可视化

通过箱线图对比源域与目标域的特征分布，观察迁移前后分布重叠度的变化，验证特征适配效果。

5.3 事后可解释性：逐例溯源

对每个目标域文件，输出诊断报告，包含：

• 投票结果与置信度。
• 关键特征值（如高频比=0.45，GOPAE=0.5）。
• LIME局部解释：各特征对决策的贡献度。
• 频谱图标注故障特征频率，与机理对照。

排列重要性分析显示，高频比和GOPAE是最重要的两个特征，打乱它们会导致准确率下降超过20%。

6 总结与推广

本文提出了一套完整的高速列车轴承故障诊断方案，创新点在于：

1. 多尺度加权GOPAE：将全局序模式注意力熵与多尺度分析结合，有效量化故障的动态复杂性，并通过权重融合提升跨工况鲁棒性。
2. 集成迁移学习：融合SVM、随机森林、梯度提升的软投票模型，结合特征迁移策略，解决了目标域无标签样本的诊断难题。
3. 全链路可解释性：从事前、过程、事后三个维度，将诊断结果与物理机理深度关联，增强了模型的信任度。

该方案已成功应用于某轨道交通企业的轴承监测系统，后续可推广至风电、航空发动机等关键旋转机械的故障诊断中。

附录：核心代码（节选）

以下为数据预处理与特征提取的关键代码（已调整变量名以适应本文逻辑）。完整代码及数据可在项目仓库中找到。

# -*- coding: utf-8 -*-
import numpy as np
import os
from scipy.io import loadmat
from scipy.signal import welch
from scipy.stats import kurtosis, skew
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier, VotingClassifier
from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold
from sklearn.metrics import classification_report, accuracy_score, confusion_matrix
from sklearn.preprocessing import RobustScaler
from sklearn.feature_selection import SelectKBest, f_classif
import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
from collections import defaultdict, Counter
from itertools import permutations
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
# ========== GOPAE相关函数 ==========
def compute_ordinal_patterns(series, embed_dim=3, time_delay=1):
    """将时间序列转换为序模式列表和位置索引"""
    patterns = []
    indices = []
    for i in range(len(series) - (embed_dim - 1) * time_delay):
        window = series[i:i + time_delay * embed_dim:time_delay]
        ranks = tuple(np.argsort(window))
        patterns.append(ranks)
        indices.append(i)
    return patterns, indices
def get_pattern_positions(patterns):
    """记录每种序模式出现的位置"""
    pos_dict = defaultdict(list)
    for idx, pat in enumerate(patterns):
        pos_dict[pat].append(idx)
    return pos_dict
def calc_attention_entropy(positions):
    """计算给定位置序列的注意力熵"""
    if len(positions) < 2:
        return 0.0
    intervals = [positions[i+1] - positions[i] for i in range(len(positions)-1)]
    freq = Counter(intervals)
    probs = np.array([count / len(intervals) for count in freq.values()])
    return -np.sum(probs * np.log2(probs))
def gopae_single_scale(ts, embed_dim=3, time_delay=1):
    """计算单个尺度的GOPAE值"""
    patterns, _ = compute_ordinal_patterns(ts, embed_dim, time_delay)
    pos_dict = get_pattern_positions(patterns)
    total = len(patterns)
    all_patterns = list(permutations(range(embed_dim)))
    gopae_val = 0.0
    for pat in all_patterns:
        pos = pos_dict.get(pat, [])
        ae = calc_attention_entropy(pos)
        prob = len(pos) / total if total > 0 else 0
        gopae_val += prob * ae
    return gopae_val
def multiscale_gopae(signal, scales=[1,2,3,4,5], embed_dim=3, time_delay=1, per_scale=True):
    """多尺度GOPAE特征提取"""
    if len(signal) < embed_dim * time_delay:
        return np.zeros(len(scales)) if per_scale else 0.0
    # 多尺度分解与计算的详细实现（完整代码见项目仓库）
    # ......
    return np.array(gopae_vals) if per_scale else np.mean(gopae_vals)
# ========== 时域特征提取 ==========
def extract_time_features(signal):
    """提取9维时域统计特征"""
    mean_val = np.mean(signal)
    std_val = np.std(signal)
    var_val = np.var(signal)
    rms_val = np.sqrt(np.mean(signal**2))
    peak_val = np.max(np.abs(signal))
    crest_factor = peak_val / rms_val if rms_val != 0 else 0
    kurt_val = kurtosis(signal)
    skew_val = skew(signal)
    shape_factor = rms_val / np.mean(np.abs(signal)) if np.mean(np.abs(signal)) != 0 else 0
    return [mean_val, std_val, var_val, rms_val, peak_val, crest_factor, kurt_val, skew_val, shape_factor]
# ========== 频域特征提取 ==========
def extract_freq_features(signal, fs, rpm):
    """提取5维频域特征，并进行转频归一化"""
    freqs, psd = welch(signal, fs, nperseg=min(1024, len(signal)))
    # 谱心、谱方差、频带能量比等计算（完整代码见项目仓库）
    # ......
    return freq_features
# ========== 主处理流程 ==========
def process_source_file(file_path, window_size=4096, overlap=0.5, scales=[1,2,3,4,5]):
    """处理单个源域文件，生成样本特征矩阵和标签"""
    data = loadmat(file_path)
    signal = data['DE_time'].flatten()
    label = extract_label(file_path)
    samples = []
    step = int(window_size * (1 - overlap))
    for start in range(0, len(signal) - window_size, step):
        segment = signal[start:start+window_size]
        time_feat = extract_time_features(segment)
        freq_feat = extract_freq_features(segment, fs=12000, rpm=1797)
        gopae_feat = multiscale_gopae(segment, scales=scales, per_scale=True)
        combined = np.concatenate([time_feat, freq_feat, gopae_feat])
        samples.append(combined)
    return np.array(samples), label
# 实际使用时，需遍历所有文件并整合特征矩阵
# ......