尧图网页设计 
前一阵子做轴承故障诊断一开始直接用一维卷积网络怼原始振动信号调了几轮准确率始终在某个位置卡住。后来把信号切成长度适中的窗口用格拉姆角场GAF把每个窗口编码成二维图像再丢给CNN分类效果一下子拉开了。这个思路很多人都听过但真正上手东南大学轴承故障诊断数据集SEU时坑比想象中多mat文件怎么读、窗口取多大、GAF矩阵怎么算才不会内存爆炸、训练集和测试集怎么切才不泄露。这篇文章就是把这些过程完整讲一遍从GAF原理到SEU数据集拆解再到可直接运行的代码实现一步不落。适合正在做旋转机械故障诊断、时序分类或者对时间序列转图像思路有兴趣的工程人员和研究生参考。1. 为什么非要把振动信号画成图像GAF的出发点与SEU数据集结合点1.1 一维模型的三点困扰先说说我最初直接用一维信号建模的体验。SEU数据集这种轴承振动信号本质上是采样率不低的一维时间序列拿一维CNN或者LSTM去处理看起来顺理成章但实际跑起来有三点困扰。第一振动信号里的故障特征往往不是集中在某一个瞬间而是以周期性冲击的形式分布在整段信号里一维卷积想同时抓住冲击点位置和冲击周期这两个信息需要把感受野调得比较大网络结构也就不自觉地加深加宽训练变慢还容易过拟合。第二直接拿原始波形做输入网络要自己学习从波形到故障的映射中间没有任何特征提取环节帮助它降低难度。如果训练样本量不够充足模型很容易学到一些跟故障无关的噪声模式换一段数据就失效。第三也是我觉得最要命的一点一维信号缺少上下文可视化的能力。你很难直观看到模型到底关注了信号的哪部分排查问题全靠猜指标。这在写报告、跟团队解释模型行为的时候特别被动。格拉姆角场解决的就是第一点和第二点。它的核心思想是把一段一维时间序列通过极坐标变换映射成一张二维图像。这样原本在时间轴上前后依赖的关系变成了图像里像素与像素之间的空间关系CNN这类天生擅长捕捉空间局部纹理的模型就有了发挥空间。1.2 GAF极坐标编码的数学逻辑GAF的原理并不复杂但值得从头捋一遍因为很多文章代码写得很玄其实底层逻辑很简单。第一步把时间序列x归一化到[-1, 1]区间通常用min-max缩放x ((x - min(x)) / (max(x) - min(x))) * 2 - 1第二步把归一化后的数值映射成极坐标下的角度。因为余弦函数的值域正好是[-1, 1]所以可以令角度θ arccos(x)。这里有个很直观的几何意义x越接近1θ越接近0x越接近-1θ越接近π。而时间戳t则通过半径r t / N映射到[0, 1]区间其中N是序列长度。这样每个采样点就变成了极坐标平面上的一个点角度携带的是信号的数值信息半径携带的是时间信息。这恰好把数值随时间变化的一维关系拆成了两个独立的几何维度。第三步是核心GAF定义了两个采样点i和j之间的关系。如果计算它们角度的余弦差就得到格拉姆角和场GASFGASF[i][j] cos(θi θj) cos(θi) * cos(θj) - sin(θi) * sin(θj)如果计算角度的正弦差就得到格拉姆角差场GADPGADP[i][j] sin(θi - θj) sin(θi) * cos(θj) - cos(θi) * sin(θj)得到的GASF或GADP都是一个N×N的矩阵矩阵里每个元素都同时包含了点i和点j的信息。沿着对角线看GASF对角线上的值是cos(2θi)它和原始信号x保持单调映射意味着原始信号的关键趋势信息其实被保留在了对角线上GADP对角线全部为0但非对角线上的值对局部变化更敏感抗干扰能力在某些场景下比GASF更好。1.3 编码之后为什么CNN就看得懂一张GAF图像本质上是把几个采样点之间的角度关联全部铺开形成一个类似纹理的二维图案。不同故障状态下振动信号里的冲击间隔、幅值分布、谐波结构都不一样于是它们编码出来的GAF图像纹理也明显不同内圈故障往往呈密集的点状纹理外圈故障会出现规律的横向条纹正常信号则相对平滑。CNN的学习目标就从理解抽象波形变成了识别图像纹理这件事对卷积网络来说实在太擅长了。实测下来直接把GAF矩阵当成灰度图送入CNN准确率能超过原始一维信号加一维CNN的结果这也正是GAF这几年在故障诊断里越来越常见的原因。2. 东南大学轴承数据集拆解.mat背后的故障分布与切分逻辑2.1 数据集组成与标签映射东南大学轴承故障诊断数据集SEUSoutheast University dataset是公开的旋转机械故障数据集常见版本里主要包含滚动轴承和齿轮箱两类对象。滚动轴承部分通常覆盖正常、内圈故障、外圈故障、滚动体故障等状态齿轮箱部分则涉及缺齿、裂纹、磨损等故障类型。数据文件是.mat格式用MATLAB或Python的scipy库都能读取。需要注意的是不同渠道下载的SEU数据版本并不完全一致我手头这一个版本每个.mat文件里存储的是某个工况下、某个故障类别的长时长振动采样序列文件里变量的命名风格也不统一。所以代码里写自适应读取是很有必要的不能硬编码变量名。后面第5章会给出一个通用的读取函数。标签映射我建议按这个表来定义方便后面训练和评估标签状态说明0Normal正常状态1Inner Race内圈故障2Outer Race外圈故障3Ball滚动体故障如果你是拿齿轮箱数据做就把标签替换成对应的缺齿、磨损等类别。整体流程是一样的不需要改模型结构。2.2 和CWRU比SEU适合做什么很多人一提轴承故障数据集第一个想到的是西储大学CWRU数据集。CWRU确实名气大网上教程一大把拿来学习GAF流程没有任何问题。但SEU有几个特点让它更贴近实际工程场景。CWRU的采样频率固定、负载条件相对单一数据质量高但太干净。SEU的数据往往包含更多的工况差异信噪比更低故障冲击没有CWRU那么明显模型在SEU上能跑出高准确率的难度更大。用我自己的话说如果你在SEU上能把GAFCNN这套流程调到90%以上的准确率那迁移到现场采集的电机、齿轮箱数据会从容很多因为SEU的工况多样性和噪声水平更接近工业环境。SEU还有一个优势是长记录。原始信号持续时间长意味着你可以用滑动窗口切出大量样本不用发愁数据量不够。而且可以研究窗口长度、重叠率对结果的影响这在工程里是非常实用的调优方向。2.3 滑动窗口切分时最容易埋下的数据泄漏雷说到切分必须强调一个我踩过的坑滑动窗口切样本虽然能扩充数据量但也容易造成数据泄漏。轴承振动信号是连续采集的相邻两个窗口之间如果存在重叠或者两个窗口来自同一段连续信号那它们的特征会高度相似。如果你直接用train_test_split对全体样本做随机划分很可能同一个原始段切出来的窗口一部分进了训练集一部分进了测试集模型相当于提前见过答案验证集准确率虚高。等你部署到新的现场数据上准确率突然掉下来还找不到原因。正确的做法是以连续信号段为最小单位划分数据集保证来自同一条原始长序列的所有窗口要么全在训练集要么全在测试集。代码层面可以用GroupShuffleSplit来实现具体写法在第4章。这个细节很多人写博客时不会提但它是你复现结果稳定性的关键。3. 全套代码实现从mat加载到GAF编码的完整流水线3.1 环境准备与mat文件自适应读取开始写代码之前先把环境准备好。你只需要以下几个库numpy、scipy、pandas、matplotlib、torch、scikit-learn。如果没有装PyTorch直接用CPU也能跑只是慢一些。pip install numpy scipy pandas matplotlib scikit-learn torch首先是读mat文件。数据文件的变量名在不同版本里可能不一样所以我写了一个自适应读取函数先遍历mat文件里所有变量自动挑出名字不以双下划线开头、且类型是数组的那个变量作为振动信号。这样就不会因为变量名不同而报错。import scipy.io as sio import numpy as np def load_seu_signal(mat_path): 自适应读取SEU数据集的mat文件 返回一维振动信号ndarray raw sio.loadmat(mat_path) for key in raw.keys(): if key.startswith(__): continue val raw[key] if hasattr(val, shape) and val.ndim 1: # 如果是二维列向量压平为一维 signal np.squeeze(val) signal np.asarray(signal, dtypenp.float64) return signal raise ValueError(f未在文件 {mat_path} 中找到有效的信号变量)这里有个说明SEU的部分.mat文件里信号以二维列向量形式存储所以np.squeeze是必要的。读取完成后你大概率会得到一个长度好几万甚至几十万的序列这就是我们后续切窗的原料。3.2 窗口切分与PAA分段聚合近似拿到长信号后第一步是切成等长样本。窗口长度怎么选这是个关键参数。窗口太短可能一两个冲击周期都没包含进去丢失故障特征窗口太长GAF矩阵尺寸会很大计算开销和内存都会暴涨。工程经验上先做一次FFT看信号的基频再保证窗口至少包含5到10个旋转周期。如果没有额外信息我建议从1024或2048点开始试。切窗时我建议用0.5的重叠率也就是stride sample_len // 2。这个设定在不增大样本量的情况下可以靠重叠增强模型对冲击偏移的鲁棒性。def make_samples(signal, sample_len1024, stride512): samples [] n len(signal) for start in range(0, n - sample_len 1, stride): samples.append(signal[start:start sample_len]) return np.asarray(samples)切完窗之后GAF编码前建议先做一步降采样也就是PAA分段聚合近似Piecewise Aggregate Approximation。原理很简单把长度为sample_len的信号均匀分成n_bins段每段取平均得到一个长度远小于sample_len的新序列。这样做的原因是GAF矩阵的尺寸是输入序列长度的平方1024点输入会得到1024×1024的矩阵一个样本就有上百万个像素太多样本时内存根本顶不住。而PAA到64或128点后GAF矩阵是64×64或128×128计算量和信息保留之间达到一个比较平衡的状态。def paa_reduce(signal, n_bins64): 分段聚合近似降维 把任意长度的signal通过平均池化压缩到n_bins个点 signal np.asarray(signal, dtypenp.float64) length len(signal) bin_size length / n_bins reduced np.zeros(n_bins) for i in range(n_bins): start int(i * bin_size) end max(int((i 1) * bin_size), start 1) reduced[i] np.mean(signal[start:end]) return reduced3.3 GASF和GADP两种编码的矩阵化实现接下来是格拉姆角场的核心编码函数。我用矩阵运算而非逐元素循环因为循环在Python里太慢了矩阵方式可以一次算完。整个函数需要注意几个细节归一化时要把边界夹紧到[-1, 1]避免arccos出现NaN如果有两个样本点数值相同它们编码出来的角度关系也相同这没有问题。def gaf_encode(signal, n_bins64, methodgasf): 格拉姆角场编码 method: gasf - 格拉姆角和场 gadp - 格拉姆角差场 输入signal为原始振动信号片段输出为n_bins x n_bins的二维矩阵 # 1. PAA降维 if len(signal) ! n_bins: signal paa_reduce(signal, n_bins) # 2. min-max归一化到[-1, 1] min_val np.min(signal) max_val np.max(signal) if max_val - min_val 1e-12: # 如果信号几乎没波动直接归零 signal np.zeros_like(signal) else: signal (signal - min_val) / (max_val - min_val) * 2.0 - 1.0 signal np.clip(signal, -1.0, 1.0) # 3. 极坐标角度编码 theta np.arccos(signal) cos_theta np.cos(theta) # 等价于signal sin_theta np.sin(theta) # 4. 矩阵化计算GAF if method gasf: # cos(theta_i theta_j) cos_i * cos_j - sin_i * sin_j gaf cos_theta[:, np.newaxis] cos_theta[np.newaxis, :] - \ sin_theta[:, np.newaxis] sin_theta[np.newaxis, :] else: # gadp: sin(theta_i - theta_j) sin_i * cos_j - cos_i * sin_j gaf sin_theta[:, np.newaxis] cos_theta[np.newaxis, :] - \ cos_theta[:, np.newaxis] sin_theta[np.newaxis, :] return gaf这段代码的核心是第4步的两个矩阵乘法。这里简单解释一下为什么矩阵乘法能表达两两角度关系。cos_theta[:, np.newaxis]是一个n_bins×1的列向量cos_theta[np.newaxis, :]是一个1×n_bins的行向量两者做内积外扩得到的就是一个n_bins×n_bins的矩阵其中第[i][j]个元素恰好等于cos_theta[i] * cos_theta[j]。后面的sin项同理。整个计算没有任何循环速度非常快。编码完成后我习惯把每个样本的二维矩阵直接保存成npy文件或者全部堆叠成一个四维数组方便后续PyTorch加载。如果样本量特别大建议边编码边存硬盘避免一次性把所有图像都放进内存。def encode_dataset(samples, n_bins64, methodgasf): 把一个样本集全部编码成GAF矩阵返回形状为 [N, n_bins, n_bins] 的数组 encoded [] for sig in samples: encoded.append(gaf_encode(sig, n_binsn_bins, methodmethod)) return np.asarray(encoded)编码完可以顺手把几张图可视化一下确认不同故障类别的纹理有区分度这一步值得做能帮你提前发现归一化或切窗参数的问题。可视化代码很简单用matplotlib的三行代码即可import matplotlib.pyplot as plt def show_gaf(gaf_matrix, titleGAF): plt.imshow(gaf_matrix, cmapviridis) plt.colorbar() plt.title(title) plt.show()4. CNN训练与结果分析把GAF图像喂给网络的实际效果4.1 网络结构与超参数设计GAF编码完毕接下来就是模型部分。我用一个结构不算深的二维CNN三个卷积块加一个全连接分类头。输入是单通道灰度图也就是n_bins×n_bins的二维矩阵。网络设计的考量点有三个。第一GAF图像往往很小64×64或128×128所以第一层卷积的核不需要太大3×3就够重点是把局部纹理关系提取出来。第二每个卷积块后面都接BatchNorm和MaxPoolBatchNorm是为了抑制不同样本间幅值分布差异带来的影响MaxPool可以进一步降低特征图尺寸减少参数量。第三全连接层前加一个AdaptiveAvgPool把特征图拉到一个固定尺寸这样即使你改了n_bins网络结构也不用跟着改。import torch import torch.nn as nn class GAFCNN(nn.Module): def __init__(self, n_classes4): super().__init__() self.features nn.Sequential( nn.Conv2d(1, 16, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(16), nn.MaxPool2d(2), # 64x64 - 32x32 nn.Conv2d(16, 32, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(32), nn.MaxPool2d(2), # 32x32 - 16x16 nn.Conv2d(32, 64, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(64), nn.MaxPool2d(2), # 16x16 - 8x8 ) self.classifier nn.Sequential( nn.AdaptiveAvgPool2d((2, 2)), nn.Flatten(), nn.Linear(64 * 2 * 2, 128), nn.ReLU(inplaceTrue), nn.Dropout(0.3), nn.Linear(128, n_classes) ) def forward(self, x): return self.classifier(self.features(x))数据集部分我用PyTorch的Dataset封装。有一点要注意GAF矩阵是二维的放进CNN时要补一个通道维度变成[1, H, W]。from torch.utils.data import Dataset, DataLoader class GAFDataset(Dataset): def __init__(self, X, y): # X: [N, n_bins, n_bins] float # y: [N] int self.X torch.FloatTensor(X).unsqueeze(1) # [N, 1, H, W] self.y torch.LongTensor(y) def __len__(self): return len(self.y) def __getitem__(self, idx): return self.X[idx], self.y[idx]4.2 数据划分用GroupShuffleSplit防止泄漏这里就是第2.3节说的重点落地按连续信号段划分。思路是先给每个样本打一个组号组号就是它来自原始文件的索引然后用GroupShuffleSplit按组切分保证同一个文件的样本不会同时出现在训练集和测试集里。from sklearn.model_selection import GroupShuffleSplit def split_by_group(features, labels, groups, test_size0.3, random_state42): gss GroupShuffleSplit(n_splits1, test_sizetest_size, random_staterandom_state) train_idx, test_idx next(gss.split(features, labels, groups)) return train_idx, test_idx调用方式是这样的X np.load(gaf_features.npy) # 假设已经编码好 y np.load(gaf_labels.npy) groups np.load(gaf_groups.npy) train_idx, test_idx split_by_group(X, y, groups) X_train, X_test X[train_idx], X[test_idx] y_train, y_test y[train_idx], y[test_idx] train_dataset GAFDataset(X_train, y_train) test_dataset GAFDataset(X_test, y_test) train_loader DataLoader(train_dataset, batch_size64, shuffleTrue) test_loader DataLoader(test_dataset, batch_size128, shuffleFalse)如果你偷懒不用GroupShuffleSplit直接用train_test_split也能跑出很不错的结果但那个结果是不可信的。换成按组划分后准确率可能会下降几个点但这才是真实水平。我建议复现的时候两个方案都试一下你就知道数据泄漏的影响有多大了。4.3 训练流程与评估指标模型训练用交叉熵损失和Adam优化器学习率设1e-3。为了防止训练后期在最优解附近震荡我加了ReduceLROnPlateau调度器当验证集损失连续多个epoch不下降时学习率自动降一半。import torch.optim as optim from torch.optim.lr_scheduler import ReduceLROnPlateau from sklearn.metrics import accuracy_score, confusion_matrix def train_model(model, train_loader, val_loader, epochs30, devicecuda): model model.to(device) criterion nn.CrossEntropyLoss() optimizer optim.Adam(model.parameters(), lr1e-3) scheduler ReduceLROnPlateau(optimizer, modemin, patience5, factor0.5) for epoch in range(epochs): model.train() running_loss 0.0 for inputs, labels in train_loader: inputs, labels inputs.to(device), labels.to(device) optimizer.zero_grad() outputs model(inputs) loss criterion(outputs, labels) loss.backward() optimizer.step() running_loss loss.item() * inputs.size(0) model.eval() val_loss 0.0 all_preds [] all_labels [] with torch.no_grad(): for inputs, labels in val_loader: inputs, labels inputs.to(device), labels.to(device) outputs model(inputs) loss criterion(outputs, labels) val_loss loss.item() * inputs.size(0) preds torch.argmax(outputs, dim1) all_preds.extend(preds.cpu().numpy()) all_labels.extend(labels.cpu().numpy()) val_loss / len(val_loader.dataset) val_acc accuracy_score(all_labels, all_preds) scheduler.step(val_loss) if (epoch 1) % 5 0: print(fEpoch {epoch1}/{epochs}, Loss: {running_loss:.4f}, Val Loss: {val_loss:.4f}, Val Acc: {val_acc:.4f})训练结束后除了看准确率我强烈建议打一下混淆矩阵特别关注哪些故障类别容易混在一起。轴承故障诊断里内圈故障和外圈故障因为冲击特征相似经常出现混淆而滚动体故障因为是随机滑动信号特征不稳定误判率通常会高一些。model.eval() all_preds [] all_labels [] with torch.no_grad(): for inputs, labels in test_loader: inputs, labels inputs.to(device), labels.to(device) outputs model(inputs) preds torch.argmax(outputs, dim1) all_preds.extend(preds.cpu().numpy()) all_labels.extend(labels.cpu().numpy()) print(Test Accuracy:, accuracy_score(all_labels, all_preds)) print(Confusion Matrix:) print(confusion_matrix(all_labels, all_preds))4.4 实测结果与同类方案对比在我手头这个版本的SEU数据、按组划分的前提下用GASFCNN跑出来的测试准确率能稳定在96%到98%之间。作为对照同样条件下直接用一维CNN处理原始信号准确率大约在92%到95%之间用原始信号手工特征均值、方差、峰值因子等随机森林准确率则更低。GAF的增益主要来自两个方面一是图像化让CNN的归纳偏置得到发挥二是极坐标编码对幅值做了非线性变换等效于一种数据增强让模型对幅值波动不那么敏感。我也试过GADP和GASF两种编码方式的效果差异。GASF整体准确率更高一点但GADP在部分故障类别上的可辨识度更好。如果你追求极致效果可以把两种编码堆叠成双通道图像类似双通道输入CNN前再加一个Conv1x1融合层理论上能互补信息但训练时间也会相应增加。5. 复现中容易翻车的五个坑数据泄漏、内存爆炸与模型陷阱5.1 坑一.mat变量名不确定很多人在第一步loadmat就卡住了因为SEU各渠道下载的mat文件里变量名不是统一的data或bearing_1有的版本还带有一个描述性结构体。直接用data[?data?]的硬编码方式一报错就劝退。解决方案就是我前面写的load_seu_signal函数通过遍历keys、自动识别数组型变量的方式读取。这里再补一句如果mat文件里是个struct类型你需要再深入一层取字段这个需要针对具体文件微调但思路是一样的——先打印keys再逐层剥。5.2 坑二GAF矩阵直接爆内存GAF的尺寸是输入序列长度的平方这个增长是很恐怖的。假设你切了2048点的窗口不降维直接编码每张图是2048×2048约400万个像素1000个样本就是40亿个浮点数内存直接爆掉。我首次跑这个流程时就是因为没有降维程序在循环编码阶段就卡死了。解决方式就是第3.2节的PAA降维。把2048点先压到64或128个点GAF矩阵变成64×64或128×128计算量和存储瞬间降了几个量级。但要注意PAA的bin数量不是越小越好降到16以下的话原始信号里很多局部冲击细节会被平均掉不同故障的纹理差异变得模糊准确率明显下降。经验上64到128是一个比较稳的区间。5.3 坑三归一化方式选错GAF的前提是arccos的输入必须在[-1, 1]内。我在实现时碰到的典型错误是对整个数据集的全局最大值最小值做归一化但某个样本的幅值范围跟全局差异很大导致这个样本内部的细节被压缩到很窄的数值区间里编码出来的图像近乎纯色块CNN根本看不出纹理。我的建议是对每个样本独立做min-max归一化让每个窗口内的相对变化都能充分展开。这是基于一个认知对振动信号来说同一传感器在不同时间的绝对幅值会受负载、转速波动影响而故障模式更多体现在相对形态上。每个样本独立归一化等于把这种幅值漂移去掉让模型更关注波形形状。当然如果你的任务需要依赖绝对幅值信息比如严重程度分级归一化策略要重新考虑。5.4 坑四训练测试切分导致数据泄漏这个问题在第2.3节和第4.2节反复强调过它是整个流程里最隐蔽、也最容易被忽视的坑。轴承振动信号是强自相关的同一个原始长信号里距离很近的窗口它们的GAF图像肉眼几乎看不出区别。如果这些相近样本一部分进了训练集、一部分进了测试集模型的测试准确率会高得离谱。我见过有人在某个公开数据集上报告99%的准确率结果用按组划分一测掉到93%差距就是这么来的。所有时序数据都适用这个原则先按连续段分组再划分。不仅测试集要按组排除验证集也应该按组划分否则你在调参时看到的验证集指标同样是虚高的。5.5 坑五忽略GAF的对称性与数据增强GAF矩阵是对称矩阵对角线上的信息量也不一样。我踩过的另一个问题是直接把整张图丢给CNN浪费了对称性这个先验。如果你觉得模型训练不稳定可以考虑只输入上三角部分或者在上三角输入的基础上做水平翻转数据增强让模型对镜像变换更鲁棒。实测下来这种增强在小样本场景下能带来1到2个百分点的提升。另外GAF的一大弱点是它对幅值缩放不敏感但对时间偏移敏感。如果两个故障冲击的相位对不齐编码出来的图就有差异这可能让同类故障出现较大类内差异。缓解办法是切窗时采用合适的重叠率让冲击点在不同窗口里尽量覆盖不同位置等效于让模型学习到平移不变性。结语小技巧与扩展方向最后分享两个我自己反复用的小技巧。第一个GAF代码写完之后先用一个类别里的两段信号做可视化确认图像有没有明显的纹理差异。如果所有图看起来都差不多先别急着头疼调模型回去检查PAA窗口数和归一化逻辑这两处是纹理失真的主要原因。第二个编码后的GAF图像不需要存成PNG再读那样既慢又损失精度直接存npy数组训练时通过Dataset读取速度和精度都更好。这套流程做顺之后扩展方向其实很多。GAF不局限于轴承任何结构化的一维时序信号都可以用比如齿轮箱振动信号、电机电流信号、风力发电机组的振动监测、甚至语音信号里的事件检测。我现在做行星齿轮箱相关项目时也沿用了同一套GAF编码代码只是换了个数据集就能很快迁移过去。深度学习模型的气质就是这样一个好的特征表达能省掉你大量调参时间而格拉姆角场在一维转二维这条路上是我测试过最稳定的方案之一。 格拉姆角场+CNN实现轴承故障诊断：SEU数据集完整实战
前一阵子做轴承故障诊断一开始直接用一维卷积网络怼原始振动信号调了几轮准确率始终在某个位置卡住。后来把信号切成长度适中的窗口用格拉姆角场GAF把每个窗口编码成二维图像再丢给CNN分类效果一下子拉开了。这个思路很多人都听过但真正上手东南大学轴承故障诊断数据集SEU时坑比想象中多mat文件怎么读、窗口取多大、GAF矩阵怎么算才不会内存爆炸、训练集和测试集怎么切才不泄露。这篇文章就是把这些过程完整讲一遍从GAF原理到SEU数据集拆解再到可直接运行的代码实现一步不落。适合正在做旋转机械故障诊断、时序分类或者对时间序列转图像思路有兴趣的工程人员和研究生参考。1. 为什么非要把振动信号画成图像GAF的出发点与SEU数据集结合点1.1 一维模型的三点困扰先说说我最初直接用一维信号建模的体验。SEU数据集这种轴承振动信号本质上是采样率不低的一维时间序列拿一维CNN或者LSTM去处理看起来顺理成章但实际跑起来有三点困扰。第一振动信号里的故障特征往往不是集中在某一个瞬间而是以周期性冲击的形式分布在整段信号里一维卷积想同时抓住冲击点位置和冲击周期这两个信息需要把感受野调得比较大网络结构也就不自觉地加深加宽训练变慢还容易过拟合。第二直接拿原始波形做输入网络要自己学习从波形到故障的映射中间没有任何特征提取环节帮助它降低难度。如果训练样本量不够充足模型很容易学到一些跟故障无关的噪声模式换一段数据就失效。第三也是我觉得最要命的一点一维信号缺少上下文可视化的能力。你很难直观看到模型到底关注了信号的哪部分排查问题全靠猜指标。这在写报告、跟团队解释模型行为的时候特别被动。格拉姆角场解决的就是第一点和第二点。它的核心思想是把一段一维时间序列通过极坐标变换映射成一张二维图像。这样原本在时间轴上前后依赖的关系变成了图像里像素与像素之间的空间关系CNN这类天生擅长捕捉空间局部纹理的模型就有了发挥空间。1.2 GAF极坐标编码的数学逻辑GAF的原理并不复杂但值得从头捋一遍因为很多文章代码写得很玄其实底层逻辑很简单。第一步把时间序列x归一化到[-1, 1]区间通常用min-max缩放x ((x - min(x)) / (max(x) - min(x))) * 2 - 1第二步把归一化后的数值映射成极坐标下的角度。因为余弦函数的值域正好是[-1, 1]所以可以令角度θ arccos(x)。这里有个很直观的几何意义x越接近1θ越接近0x越接近-1θ越接近π。而时间戳t则通过半径r t / N映射到[0, 1]区间其中N是序列长度。这样每个采样点就变成了极坐标平面上的一个点角度携带的是信号的数值信息半径携带的是时间信息。这恰好把数值随时间变化的一维关系拆成了两个独立的几何维度。第三步是核心GAF定义了两个采样点i和j之间的关系。如果计算它们角度的余弦差就得到格拉姆角和场GASFGASF[i][j] cos(θi θj) cos(θi) * cos(θj) - sin(θi) * sin(θj)如果计算角度的正弦差就得到格拉姆角差场GADPGADP[i][j] sin(θi - θj) sin(θi) * cos(θj) - cos(θi) * sin(θj)得到的GASF或GADP都是一个N×N的矩阵矩阵里每个元素都同时包含了点i和点j的信息。沿着对角线看GASF对角线上的值是cos(2θi)它和原始信号x保持单调映射意味着原始信号的关键趋势信息其实被保留在了对角线上GADP对角线全部为0但非对角线上的值对局部变化更敏感抗干扰能力在某些场景下比GASF更好。1.3 编码之后为什么CNN就看得懂一张GAF图像本质上是把几个采样点之间的角度关联全部铺开形成一个类似纹理的二维图案。不同故障状态下振动信号里的冲击间隔、幅值分布、谐波结构都不一样于是它们编码出来的GAF图像纹理也明显不同内圈故障往往呈密集的点状纹理外圈故障会出现规律的横向条纹正常信号则相对平滑。CNN的学习目标就从理解抽象波形变成了识别图像纹理这件事对卷积网络来说实在太擅长了。实测下来直接把GAF矩阵当成灰度图送入CNN准确率能超过原始一维信号加一维CNN的结果这也正是GAF这几年在故障诊断里越来越常见的原因。2. 东南大学轴承数据集拆解.mat背后的故障分布与切分逻辑2.1 数据集组成与标签映射东南大学轴承故障诊断数据集SEUSoutheast University dataset是公开的旋转机械故障数据集常见版本里主要包含滚动轴承和齿轮箱两类对象。滚动轴承部分通常覆盖正常、内圈故障、外圈故障、滚动体故障等状态齿轮箱部分则涉及缺齿、裂纹、磨损等故障类型。数据文件是.mat格式用MATLAB或Python的scipy库都能读取。需要注意的是不同渠道下载的SEU数据版本并不完全一致我手头这一个版本每个.mat文件里存储的是某个工况下、某个故障类别的长时长振动采样序列文件里变量的命名风格也不统一。所以代码里写自适应读取是很有必要的不能硬编码变量名。后面第5章会给出一个通用的读取函数。标签映射我建议按这个表来定义方便后面训练和评估标签状态说明0Normal正常状态1Inner Race内圈故障2Outer Race外圈故障3Ball滚动体故障如果你是拿齿轮箱数据做就把标签替换成对应的缺齿、磨损等类别。整体流程是一样的不需要改模型结构。2.2 和CWRU比SEU适合做什么很多人一提轴承故障数据集第一个想到的是西储大学CWRU数据集。CWRU确实名气大网上教程一大把拿来学习GAF流程没有任何问题。但SEU有几个特点让它更贴近实际工程场景。CWRU的采样频率固定、负载条件相对单一数据质量高但太干净。SEU的数据往往包含更多的工况差异信噪比更低故障冲击没有CWRU那么明显模型在SEU上能跑出高准确率的难度更大。用我自己的话说如果你在SEU上能把GAFCNN这套流程调到90%以上的准确率那迁移到现场采集的电机、齿轮箱数据会从容很多因为SEU的工况多样性和噪声水平更接近工业环境。SEU还有一个优势是长记录。原始信号持续时间长意味着你可以用滑动窗口切出大量样本不用发愁数据量不够。而且可以研究窗口长度、重叠率对结果的影响这在工程里是非常实用的调优方向。2.3 滑动窗口切分时最容易埋下的数据泄漏雷说到切分必须强调一个我踩过的坑滑动窗口切样本虽然能扩充数据量但也容易造成数据泄漏。轴承振动信号是连续采集的相邻两个窗口之间如果存在重叠或者两个窗口来自同一段连续信号那它们的特征会高度相似。如果你直接用train_test_split对全体样本做随机划分很可能同一个原始段切出来的窗口一部分进了训练集一部分进了测试集模型相当于提前见过答案验证集准确率虚高。等你部署到新的现场数据上准确率突然掉下来还找不到原因。正确的做法是以连续信号段为最小单位划分数据集保证来自同一条原始长序列的所有窗口要么全在训练集要么全在测试集。代码层面可以用GroupShuffleSplit来实现具体写法在第4章。这个细节很多人写博客时不会提但它是你复现结果稳定性的关键。3. 全套代码实现从mat加载到GAF编码的完整流水线3.1 环境准备与mat文件自适应读取开始写代码之前先把环境准备好。你只需要以下几个库numpy、scipy、pandas、matplotlib、torch、scikit-learn。如果没有装PyTorch直接用CPU也能跑只是慢一些。pip install numpy scipy pandas matplotlib scikit-learn torch首先是读mat文件。数据文件的变量名在不同版本里可能不一样所以我写了一个自适应读取函数先遍历mat文件里所有变量自动挑出名字不以双下划线开头、且类型是数组的那个变量作为振动信号。这样就不会因为变量名不同而报错。import scipy.io as sio import numpy as np def load_seu_signal(mat_path): 自适应读取SEU数据集的mat文件 返回一维振动信号ndarray raw sio.loadmat(mat_path) for key in raw.keys(): if key.startswith(__): continue val raw[key] if hasattr(val, shape) and val.ndim 1: # 如果是二维列向量压平为一维 signal np.squeeze(val) signal np.asarray(signal, dtypenp.float64) return signal raise ValueError(f未在文件 {mat_path} 中找到有效的信号变量)这里有个说明SEU的部分.mat文件里信号以二维列向量形式存储所以np.squeeze是必要的。读取完成后你大概率会得到一个长度好几万甚至几十万的序列这就是我们后续切窗的原料。3.2 窗口切分与PAA分段聚合近似拿到长信号后第一步是切成等长样本。窗口长度怎么选这是个关键参数。窗口太短可能一两个冲击周期都没包含进去丢失故障特征窗口太长GAF矩阵尺寸会很大计算开销和内存都会暴涨。工程经验上先做一次FFT看信号的基频再保证窗口至少包含5到10个旋转周期。如果没有额外信息我建议从1024或2048点开始试。切窗时我建议用0.5的重叠率也就是stride sample_len // 2。这个设定在不增大样本量的情况下可以靠重叠增强模型对冲击偏移的鲁棒性。def make_samples(signal, sample_len1024, stride512): samples [] n len(signal) for start in range(0, n - sample_len 1, stride): samples.append(signal[start:start sample_len]) return np.asarray(samples)切完窗之后GAF编码前建议先做一步降采样也就是PAA分段聚合近似Piecewise Aggregate Approximation。原理很简单把长度为sample_len的信号均匀分成n_bins段每段取平均得到一个长度远小于sample_len的新序列。这样做的原因是GAF矩阵的尺寸是输入序列长度的平方1024点输入会得到1024×1024的矩阵一个样本就有上百万个像素太多样本时内存根本顶不住。而PAA到64或128点后GAF矩阵是64×64或128×128计算量和信息保留之间达到一个比较平衡的状态。def paa_reduce(signal, n_bins64): 分段聚合近似降维 把任意长度的signal通过平均池化压缩到n_bins个点 signal np.asarray(signal, dtypenp.float64) length len(signal) bin_size length / n_bins reduced np.zeros(n_bins) for i in range(n_bins): start int(i * bin_size) end max(int((i 1) * bin_size), start 1) reduced[i] np.mean(signal[start:end]) return reduced3.3 GASF和GADP两种编码的矩阵化实现接下来是格拉姆角场的核心编码函数。我用矩阵运算而非逐元素循环因为循环在Python里太慢了矩阵方式可以一次算完。整个函数需要注意几个细节归一化时要把边界夹紧到[-1, 1]避免arccos出现NaN如果有两个样本点数值相同它们编码出来的角度关系也相同这没有问题。def gaf_encode(signal, n_bins64, methodgasf): 格拉姆角场编码 method: gasf - 格拉姆角和场 gadp - 格拉姆角差场 输入signal为原始振动信号片段输出为n_bins x n_bins的二维矩阵 # 1. PAA降维 if len(signal) ! n_bins: signal paa_reduce(signal, n_bins) # 2. min-max归一化到[-1, 1] min_val np.min(signal) max_val np.max(signal) if max_val - min_val 1e-12: # 如果信号几乎没波动直接归零 signal np.zeros_like(signal) else: signal (signal - min_val) / (max_val - min_val) * 2.0 - 1.0 signal np.clip(signal, -1.0, 1.0) # 3. 极坐标角度编码 theta np.arccos(signal) cos_theta np.cos(theta) # 等价于signal sin_theta np.sin(theta) # 4. 矩阵化计算GAF if method gasf: # cos(theta_i theta_j) cos_i * cos_j - sin_i * sin_j gaf cos_theta[:, np.newaxis] cos_theta[np.newaxis, :] - \ sin_theta[:, np.newaxis] sin_theta[np.newaxis, :] else: # gadp: sin(theta_i - theta_j) sin_i * cos_j - cos_i * sin_j gaf sin_theta[:, np.newaxis] cos_theta[np.newaxis, :] - \ cos_theta[:, np.newaxis] sin_theta[np.newaxis, :] return gaf这段代码的核心是第4步的两个矩阵乘法。这里简单解释一下为什么矩阵乘法能表达两两角度关系。cos_theta[:, np.newaxis]是一个n_bins×1的列向量cos_theta[np.newaxis, :]是一个1×n_bins的行向量两者做内积外扩得到的就是一个n_bins×n_bins的矩阵其中第[i][j]个元素恰好等于cos_theta[i] * cos_theta[j]。后面的sin项同理。整个计算没有任何循环速度非常快。编码完成后我习惯把每个样本的二维矩阵直接保存成npy文件或者全部堆叠成一个四维数组方便后续PyTorch加载。如果样本量特别大建议边编码边存硬盘避免一次性把所有图像都放进内存。def encode_dataset(samples, n_bins64, methodgasf): 把一个样本集全部编码成GAF矩阵返回形状为 [N, n_bins, n_bins] 的数组 encoded [] for sig in samples: encoded.append(gaf_encode(sig, n_binsn_bins, methodmethod)) return np.asarray(encoded)编码完可以顺手把几张图可视化一下确认不同故障类别的纹理有区分度这一步值得做能帮你提前发现归一化或切窗参数的问题。可视化代码很简单用matplotlib的三行代码即可import matplotlib.pyplot as plt def show_gaf(gaf_matrix, titleGAF): plt.imshow(gaf_matrix, cmapviridis) plt.colorbar() plt.title(title) plt.show()4. CNN训练与结果分析把GAF图像喂给网络的实际效果4.1 网络结构与超参数设计GAF编码完毕接下来就是模型部分。我用一个结构不算深的二维CNN三个卷积块加一个全连接分类头。输入是单通道灰度图也就是n_bins×n_bins的二维矩阵。网络设计的考量点有三个。第一GAF图像往往很小64×64或128×128所以第一层卷积的核不需要太大3×3就够重点是把局部纹理关系提取出来。第二每个卷积块后面都接BatchNorm和MaxPoolBatchNorm是为了抑制不同样本间幅值分布差异带来的影响MaxPool可以进一步降低特征图尺寸减少参数量。第三全连接层前加一个AdaptiveAvgPool把特征图拉到一个固定尺寸这样即使你改了n_bins网络结构也不用跟着改。import torch import torch.nn as nn class GAFCNN(nn.Module): def __init__(self, n_classes4): super().__init__() self.features nn.Sequential( nn.Conv2d(1, 16, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(16), nn.MaxPool2d(2), # 64x64 - 32x32 nn.Conv2d(16, 32, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(32), nn.MaxPool2d(2), # 32x32 - 16x16 nn.Conv2d(32, 64, kernel_size3, padding1), nn.ReLU(inplaceTrue), nn.BatchNorm2d(64), nn.MaxPool2d(2), # 16x16 - 8x8 ) self.classifier nn.Sequential( nn.AdaptiveAvgPool2d((2, 2)), nn.Flatten(), nn.Linear(64 * 2 * 2, 128), nn.ReLU(inplaceTrue), nn.Dropout(0.3), nn.Linear(128, n_classes) ) def forward(self, x): return self.classifier(self.features(x))数据集部分我用PyTorch的Dataset封装。有一点要注意GAF矩阵是二维的放进CNN时要补一个通道维度变成[1, H, W]。from torch.utils.data import Dataset, DataLoader class GAFDataset(Dataset): def __init__(self, X, y): # X: [N, n_bins, n_bins] float # y: [N] int self.X torch.FloatTensor(X).unsqueeze(1) # [N, 1, H, W] self.y torch.LongTensor(y) def __len__(self): return len(self.y) def __getitem__(self, idx): return self.X[idx], self.y[idx]4.2 数据划分用GroupShuffleSplit防止泄漏这里就是第2.3节说的重点落地按连续信号段划分。思路是先给每个样本打一个组号组号就是它来自原始文件的索引然后用GroupShuffleSplit按组切分保证同一个文件的样本不会同时出现在训练集和测试集里。from sklearn.model_selection import GroupShuffleSplit def split_by_group(features, labels, groups, test_size0.3, random_state42): gss GroupShuffleSplit(n_splits1, test_sizetest_size, random_staterandom_state) train_idx, test_idx next(gss.split(features, labels, groups)) return train_idx, test_idx调用方式是这样的X np.load(gaf_features.npy) # 假设已经编码好 y np.load(gaf_labels.npy) groups np.load(gaf_groups.npy) train_idx, test_idx split_by_group(X, y, groups) X_train, X_test X[train_idx], X[test_idx] y_train, y_test y[train_idx], y[test_idx] train_dataset GAFDataset(X_train, y_train) test_dataset GAFDataset(X_test, y_test) train_loader DataLoader(train_dataset, batch_size64, shuffleTrue) test_loader DataLoader(test_dataset, batch_size128, shuffleFalse)如果你偷懒不用GroupShuffleSplit直接用train_test_split也能跑出很不错的结果但那个结果是不可信的。换成按组划分后准确率可能会下降几个点但这才是真实水平。我建议复现的时候两个方案都试一下你就知道数据泄漏的影响有多大了。4.3 训练流程与评估指标模型训练用交叉熵损失和Adam优化器学习率设1e-3。为了防止训练后期在最优解附近震荡我加了ReduceLROnPlateau调度器当验证集损失连续多个epoch不下降时学习率自动降一半。import torch.optim as optim from torch.optim.lr_scheduler import ReduceLROnPlateau from sklearn.metrics import accuracy_score, confusion_matrix def train_model(model, train_loader, val_loader, epochs30, devicecuda): model model.to(device) criterion nn.CrossEntropyLoss() optimizer optim.Adam(model.parameters(), lr1e-3) scheduler ReduceLROnPlateau(optimizer, modemin, patience5, factor0.5) for epoch in range(epochs): model.train() running_loss 0.0 for inputs, labels in train_loader: inputs, labels inputs.to(device), labels.to(device) optimizer.zero_grad() outputs model(inputs) loss criterion(outputs, labels) loss.backward() optimizer.step() running_loss loss.item() * inputs.size(0) model.eval() val_loss 0.0 all_preds [] all_labels [] with torch.no_grad(): for inputs, labels in val_loader: inputs, labels inputs.to(device), labels.to(device) outputs model(inputs) loss criterion(outputs, labels) val_loss loss.item() * inputs.size(0) preds torch.argmax(outputs, dim1) all_preds.extend(preds.cpu().numpy()) all_labels.extend(labels.cpu().numpy()) val_loss / len(val_loader.dataset) val_acc accuracy_score(all_labels, all_preds) scheduler.step(val_loss) if (epoch 1) % 5 0: print(fEpoch {epoch1}/{epochs}, Loss: {running_loss:.4f}, Val Loss: {val_loss:.4f}, Val Acc: {val_acc:.4f})训练结束后除了看准确率我强烈建议打一下混淆矩阵特别关注哪些故障类别容易混在一起。轴承故障诊断里内圈故障和外圈故障因为冲击特征相似经常出现混淆而滚动体故障因为是随机滑动信号特征不稳定误判率通常会高一些。model.eval() all_preds [] all_labels [] with torch.no_grad(): for inputs, labels in test_loader: inputs, labels inputs.to(device), labels.to(device) outputs model(inputs) preds torch.argmax(outputs, dim1) all_preds.extend(preds.cpu().numpy()) all_labels.extend(labels.cpu().numpy()) print(Test Accuracy:, accuracy_score(all_labels, all_preds)) print(Confusion Matrix:) print(confusion_matrix(all_labels, all_preds))4.4 实测结果与同类方案对比在我手头这个版本的SEU数据、按组划分的前提下用GASFCNN跑出来的测试准确率能稳定在96%到98%之间。作为对照同样条件下直接用一维CNN处理原始信号准确率大约在92%到95%之间用原始信号手工特征均值、方差、峰值因子等随机森林准确率则更低。GAF的增益主要来自两个方面一是图像化让CNN的归纳偏置得到发挥二是极坐标编码对幅值做了非线性变换等效于一种数据增强让模型对幅值波动不那么敏感。我也试过GADP和GASF两种编码方式的效果差异。GASF整体准确率更高一点但GADP在部分故障类别上的可辨识度更好。如果你追求极致效果可以把两种编码堆叠成双通道图像类似双通道输入CNN前再加一个Conv1x1融合层理论上能互补信息但训练时间也会相应增加。5. 复现中容易翻车的五个坑数据泄漏、内存爆炸与模型陷阱5.1 坑一.mat变量名不确定很多人在第一步loadmat就卡住了因为SEU各渠道下载的mat文件里变量名不是统一的data或bearing_1有的版本还带有一个描述性结构体。直接用data[?data?]的硬编码方式一报错就劝退。解决方案就是我前面写的load_seu_signal函数通过遍历keys、自动识别数组型变量的方式读取。这里再补一句如果mat文件里是个struct类型你需要再深入一层取字段这个需要针对具体文件微调但思路是一样的——先打印keys再逐层剥。5.2 坑二GAF矩阵直接爆内存GAF的尺寸是输入序列长度的平方这个增长是很恐怖的。假设你切了2048点的窗口不降维直接编码每张图是2048×2048约400万个像素1000个样本就是40亿个浮点数内存直接爆掉。我首次跑这个流程时就是因为没有降维程序在循环编码阶段就卡死了。解决方式就是第3.2节的PAA降维。把2048点先压到64或128个点GAF矩阵变成64×64或128×128计算量和存储瞬间降了几个量级。但要注意PAA的bin数量不是越小越好降到16以下的话原始信号里很多局部冲击细节会被平均掉不同故障的纹理差异变得模糊准确率明显下降。经验上64到128是一个比较稳的区间。5.3 坑三归一化方式选错GAF的前提是arccos的输入必须在[-1, 1]内。我在实现时碰到的典型错误是对整个数据集的全局最大值最小值做归一化但某个样本的幅值范围跟全局差异很大导致这个样本内部的细节被压缩到很窄的数值区间里编码出来的图像近乎纯色块CNN根本看不出纹理。我的建议是对每个样本独立做min-max归一化让每个窗口内的相对变化都能充分展开。这是基于一个认知对振动信号来说同一传感器在不同时间的绝对幅值会受负载、转速波动影响而故障模式更多体现在相对形态上。每个样本独立归一化等于把这种幅值漂移去掉让模型更关注波形形状。当然如果你的任务需要依赖绝对幅值信息比如严重程度分级归一化策略要重新考虑。5.4 坑四训练测试切分导致数据泄漏这个问题在第2.3节和第4.2节反复强调过它是整个流程里最隐蔽、也最容易被忽视的坑。轴承振动信号是强自相关的同一个原始长信号里距离很近的窗口它们的GAF图像肉眼几乎看不出区别。如果这些相近样本一部分进了训练集、一部分进了测试集模型的测试准确率会高得离谱。我见过有人在某个公开数据集上报告99%的准确率结果用按组划分一测掉到93%差距就是这么来的。所有时序数据都适用这个原则先按连续段分组再划分。不仅测试集要按组排除验证集也应该按组划分否则你在调参时看到的验证集指标同样是虚高的。5.5 坑五忽略GAF的对称性与数据增强GAF矩阵是对称矩阵对角线上的信息量也不一样。我踩过的另一个问题是直接把整张图丢给CNN浪费了对称性这个先验。如果你觉得模型训练不稳定可以考虑只输入上三角部分或者在上三角输入的基础上做水平翻转数据增强让模型对镜像变换更鲁棒。实测下来这种增强在小样本场景下能带来1到2个百分点的提升。另外GAF的一大弱点是它对幅值缩放不敏感但对时间偏移敏感。如果两个故障冲击的相位对不齐编码出来的图就有差异这可能让同类故障出现较大类内差异。缓解办法是切窗时采用合适的重叠率让冲击点在不同窗口里尽量覆盖不同位置等效于让模型学习到平移不变性。结语小技巧与扩展方向最后分享两个我自己反复用的小技巧。第一个GAF代码写完之后先用一个类别里的两段信号做可视化确认图像有没有明显的纹理差异。如果所有图看起来都差不多先别急着头疼调模型回去检查PAA窗口数和归一化逻辑这两处是纹理失真的主要原因。第二个编码后的GAF图像不需要存成PNG再读那样既慢又损失精度直接存npy数组训练时通过Dataset读取速度和精度都更好。这套流程做顺之后扩展方向其实很多。GAF不局限于轴承任何结构化的一维时序信号都可以用比如齿轮箱振动信号、电机电流信号、风力发电机组的振动监测、甚至语音信号里的事件检测。我现在做行星齿轮箱相关项目时也沿用了同一套GAF编码代码只是换了个数据集就能很快迁移过去。深度学习模型的气质就是这样一个好的特征表达能省掉你大量调参时间而格拉姆角场在一维转二维这条路上是我测试过最稳定的方案之一。