天府杯A题仪器故障诊断全流程baseline:从SG去噪、时域特征到K-means与MLP
简介这份资源是2022年第二届天府杯全国大学生数学建模竞赛A题的完整参赛论文作者为参赛队#A485成员面向备战数学建模竞赛的高校学生及对仪器故障智能诊断感兴趣的机器学习初学者。论文围绕基于图像识别与多模型建模的故障检测技术展开系统覆盖信号去噪、时域特征提取、无监督与有监督分类建模以及迁移学习图像分类等核心环节可作为赛题复盘与建模思路参考。资源包内含1个PDF文件大小约1.29MB完整呈现从问题重述、模型假设到各问求解的论文结构。目前已有274人学习下载。读者可从中获取滑动平均、Savitzky-Golay滤波与小波变换的去噪对比方案10类时域特征的提取与相关性分析K-means、DBSCAN、Spectral聚类及SVM、Xgboost、MLP等八种有监督算法的实验设计以及AlexNet与ResNet-18微调策略在故障图像分类中的应用对掌握多模型建模与迁移学习流程具有实际参考价值。1. 天府杯A题仪器故障诊断一份能跑通全流程的建模资源去年带学弟复盘天府杯赛题时他盯着附件里那堆一维信号数据发懵——A类按行存、B类按列存连读进来都费劲更别说后面去噪、提特征、跑分类了。这份2022年第二届天府杯A题作品参赛队号#A485恰好把这条链路走通了从txt矩阵化、Savitzky-Golay去噪、10类时域特征提取到K-means无监督聚类、8种有监督模型对比最后还把一维信号reshape成64×64图像喂给AlexNet和ResNet-18做迁移学习。适合正在做信号类故障诊断、想找一套完整baseline的建模选手也适合想看清特征工程到底比原始数据强多少的从业者。它不只是一篇论文附件里还有表1-1到表4的支撑材料和实验代码能直接对着复现。2. 数据预处理与去噪从txt到矩阵再到SG滤波2.1 为什么第一步必须做数据矩阵化拿到附件先别急着上模型。A类数据按行存储、B类按列存储这个差异如果不处理后面所有实验都是在错误的数据结构上跑。论文里的做法是把两类数据统一转成100×4096的矩阵——100个样本每个样本4096个采样点。这个维度信息很关键4096点意味着一维信号长度固定后面做reshape成64×64图像时刚好能整除说明出题方在设计数据时就埋了图像化的伏笔。常见做法是用numpy的loadtxt或pandas读入后判断shape如果列数远大于行数就转置。我一般会先打印前5行和前5列确认存储方向别信文件名里的暗示。import numpy as np def load_signal_data(filepath, expected_samples100, expected_length4096): 读取故障信号txt自动判断行列存储方向并转为 (100, 4096) 矩阵 raw np.loadtxt(filepath) # 如果行数远小于列数说明是按列存储需要转置 if raw.shape[0] raw.shape[1]: raw raw.T # 校验维度是否符合预期 assert raw.shape (expected_samples, expected_length), \ f维度异常: {raw.shape}, 期望 ({expected_samples}, {expected_length}) return raw # A类按行存储B类按列存储统一处理 data_A load_signal_data(attachment_A.txt) data_B load_signal_data(attachment_B.txt) print(fA类: {data_A.shape}, B类: {data_B.shape})逻辑说明loadtxt默认按空白字符分隔读入返回二维数组。判断shape[0] shape[1]是因为按列存储时列数4096会远大于行数100。assert这行别省——我见过有人转置反了还跑了半天模型结果准确率一直在50%晃悠排查了两小时才发现数据是躺着的。参数说明expected_samples和expected_length根据附件实际维度调整论文里是100和4096。如果你的数据是其他维度改这两个参数即可但要注意后续reshape图像时4096必须能开方。2.2 五种去噪方法怎么选SG滤波胜出的数据依据去噪这步最容易翻车的地方是凭感觉选方法。论文对比了滑动平均、SG滤波、小波固定阈值、小波软阈值、小波硬阈值五种策略用了6个指标来评判MSE、SSE、RMSE是题目给定的SNR、平滑度R、互相关系数ρ是自行补充的。这个补充思路值得学——题目只给3个指标时自己加指标要覆盖不同维度SNR看信号保留程度R看平滑程度ρ看波形相关性。看A类第一个样本的数据SG滤波的MSE是3.62e-7SNR达到54.94ρ为1.0000小波固定阈值MSE是1.05e-6SNR是50.30。B类样本差距更明显SG滤波SNR高达81.73而软阈值只有25.05。SG滤波在两类数据上都拿到了最高SNR和最低MSE这不是偶然——SG滤波基于局部多项式最小二乘拟合对震荡型信号既能去噪又能保留峰值特征而小波阈值法在阈值选取不当时容易把故障冲击成分也滤掉。from scipy.signal import savgol_filter from scipy.stats import pearsonr def evaluate_denoise(original, denoised): 计算6项去噪评价指标 n len(original) mse np.mean((original - denoised) ** 2) sse np.sum((original - denoised) ** 2) rmse np.sqrt(mse) # 信噪比 snr 10 * np.log10(np.sum(original**2) / np.sum((original - denoised)**2)) # 平滑度去噪后差分的方差根 / 原始差分的方差根 r np.sqrt(np.sum(np.diff(denoised)**2) / np.sum(np.diff(original)**2)) # 互相关系数 rho, _ pearsonr(original, denoised) return {MSE: mse, SSE: sse, RMSE: rmse, SNR: snr, R: r, rho: rho} # SG滤波窗口长度取11多项式阶数取3 denoised savgol_filter(data_A[0], window_length11, polyorder3) metrics evaluate_denoise(data_A[0], denoised) for k, v in metrics.items(): print(f{k}: {v:.6f})逻辑说明savgol_filter的window_length必须是奇数polyorder必须小于window_length。论文没写具体参数但按4096点信号和故障冲击特征窗口11、阶数3是常见起点。pearsonr返回相关系数和p值这里只取相关系数。参数说明window_length越大去噪越强但可能过度平滑建议从7开始试到21polyorder一般取2到4阶数越高越能保留信号细节但去噪能力下降。如果SNR上不去先调窗口长度别急着换方法。提示去噪后一定要把数据存下来后面特征提取和无监督实验都要用。别每次重新算SG滤波在4096点上跑100个样本也要几秒。3. 时域特征提取10类判据的物理含义与计算3.1 为什么选这10类特征而不是别的论文选了均值、标准差、方根幅值、均方根、偏度、峭度、峰值因子、裕度因子、波形因子、脉冲指数。这10个不是随便凑的——它们覆盖了信号的一阶统计量均值、二阶统计量标准差、均方根、三阶统计量偏度反映分布不对称性、四阶统计量峭度反映冲击成分以及无量纲指标峰值因子、裕度因子等。对故障诊断来说峭度和脉冲指数对冲击性故障特别敏感而均方根反映能量水平。从论文给出的特征值表能看出区分度A类数据峭度0.084B类数据峭度9.132差了100多倍脉冲指数A类是3.757B类是11.426。这两个特征单独拿出来做阈值分类都能有不错的效果。但注意特征不是越多越好——论文在问题五里专门讨论了特征数量和质量的权衡后面会细说。from scipy.stats import skew, kurtosis def extract_time_features(signal): 提取10类时域特征判据 n len(signal) abs_signal np.abs(signal) mean_val np.mean(signal) std_val np.std(signal) # 方根幅值 sqrt_amp (np.mean(np.sqrt(abs_signal))) ** 2 # 均方根 rms np.sqrt(np.mean(signal ** 2)) # 偏度 skewness skew(signal) # 峭度 kurt kurtosis(signal) # 峰值因子 peak np.max(abs_signal) crest_factor peak / rms # 裕度因子 margin_factor peak / sqrt_amp # 波形因子 shape_factor rms / np.mean(abs_signal) # 脉冲指数 impulse_factor peak / np.mean(abs_signal) return { 均值: mean_val, 标准差: std_val, 方根幅值: sqrt_amp, 均方根: rms, 偏度: skewness, 峭度: kurt, 峰值因子: crest_factor, 裕度因子: margin_factor, 波形因子: shape_factor, 脉冲指数: impulse_factor } # 对去噪后的全部样本提取特征 features_A [extract_time_features(s) for s in denoised_A] features_B [extract_time_features(s) for s in denoised_B]逻辑说明skew和kurtosis来自scipy.stats注意kurtosis默认返回超额峭度正态分布为0论文里的峭度值0.084和9.132应该也是超额峭度。方根幅值的计算是先开方再平方别写成先平方再开方。参数说明这些特征都是无量纲或固定量纲不需要额外归一化就能跨样本比较。但如果要喂给SVM或KNN建议还是做标准化因为均值和均方根的量纲差异可能影响距离计算。3.2 特征分布分析与相关性热力图提取完特征别直接扔进模型。论文先对A类和B类数据的各特征分布做了分析又画了特征间相关性热力图。这一步的价值在于如果两个特征相关系数超过0.95说明信息冗余可以考虑删一个。比如均方根和方根幅值在平稳信号里相关性很高但在冲击信号里可能分开。我一般会先算特征矩阵的相关系数矩阵把相关系数大于0.9的特征对列出来人工判断是否保留。论文里没写具体删了哪些但从最终10个特征都保留来看它们之间的相关性应该都在可接受范围内。import pandas as pd import seaborn as sns # 构建特征DataFrame df_A pd.DataFrame(features_A) df_B pd.DataFrame(features_B) # 计算相关系数矩阵 corr_A df_A.corr() # 找出高相关特征对 high_corr [(i, j, corr_A.loc[i, j]) for i in corr_A.columns for j in corr_A.columns if i j and abs(corr_A.loc[i, j]) 0.9] print(A类高相关特征对:, high_corr)逻辑说明corr()默认算Pearson相关系数。i j避免重复对和自相关。如果高相关特征对很多说明特征集需要精简。参数说明阈值0.9是经验值也可以设0.85更严格。但别设太低否则会把互补特征也删掉。4. 无监督与有监督模型从K-means到MLP的选型逻辑4.1 无监督聚类K-means为什么能赢DBSCAN和Spectral问题三要求无监督二分类准确率均值90%以上、标准差10以内。论文试了K-means、DBSCAN、Spectral三种算法还用t-SNE做了可视化。结论是K-means配合提取的特征数据效果最好。这个结论有数据支撑DBSCAN对参数eps和min_samples极其敏感故障信号在高维特征空间里密度不均匀很难找到一组参数同时适配A类和B类Spectral聚类依赖相似度矩阵的构建在样本量只有100时图结构的稳定性不够。K-means虽然假设簇是凸形的但经过特征提取后A类和B类在特征空间里已经分得很开了凸形假设基本满足。论文还做了一个关键实验对比原始数据输入和特征数据输入对聚类效果的影响。t-SNE可视化显示用特征数据时两类分得很开用原始4096维数据时混在一起。这直接证明了特征提取的必要性。from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score, adjusted_rand_score from sklearn.preprocessing import StandardScaler # 合并A、B特征构建无标签数据集 X np.vstack([df_A.values, df_B.values]) # 标准化 X_scaled StandardScaler().fit_transform(X) # K-means二分类 kmeans KMeans(n_clusters2, n_init10, random_state42) labels kmeans.fit_predict(X_scaled) # 评价指标 sil silhouette_score(X_scaled, labels) # 伪标签映射与A类第一个样本最近的簇中心为A类 center_dists np.linalg.norm(kmeans.cluster_centers_ - X_scaled[0], axis1) a_cluster np.argmin(center_dists) pred (labels a_cluster).astype(int) true np.array([0]*len(df_A) [1]*len(df_B)) ari adjusted_rand_score(true, pred) print(f轮廓系数: {sil:.4f}, 调整兰德系数: {ari:.4f})逻辑说明n_init10表示跑10次不同初始化取最优避免局部最优。伪标签映射是无监督分类的关键一步——聚类输出的簇编号是随机的需要根据已知样本A类第一个来对齐。参数说明n_clusters2对应二分类任务。StandardScaler对K-means很重要因为K-means基于欧氏距离量纲差异会主导距离计算。4.2 有监督分类8种模型对比与MLP的胜出问题四要求准确率均值95%以上、标准差5以内。论文试了SVM、Xgboost、RF、MLP、DT、GaussianNB、KNN、Adaboost八种模型还对比了20%、50%、80%三种训练样本比例。看80%训练样本的结果所有8个模型的准确率都是1.0000F1、Precision、AUC全是1.0。这说明在特征提取到位的前提下这个二分类任务对模型选择不敏感。但论文最终选MLP我猜是因为MLP在20%训练样本时就已经达到1.0准确率而其他模型在20%时还有波动——MLP的样本效率最高。更值得看的是实验设置4用原始4096维数据直接输入80%训练样本。结果SVM准确率只有0.60KNN只有0.35GaussianNB和Adaboost倒是1.0但耗时明显增加。这个对比太有说服力了——特征提取把4096维降到10维不仅准确率上去了计算量还降了两个数量级。from sklearn.neural_network import MLPClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import accuracy_score, f1_score, roc_auc_score # 构建有标签数据集 X np.vstack([df_A.values, df_B.values]) y np.array([0]*len(df_A) [1]*len(df_B)) X_train, X_test, y_train, y_test train_test_split( X, y, train_size0.8, random_state42, stratifyy ) # MLP建模 mlp MLPClassifier(hidden_layer_sizes(64, 32), max_iter500, random_state42) mlp.fit(X_train, y_train) y_pred mlp.predict(X_test) y_prob mlp.predict_proba(X_test)[:, 1] print(f准确率: {accuracy_score(y_test, y_pred):.4f}) print(fF1: {f1_score(y_test, y_pred):.4f}) print(fAUC: {roc_auc_score(y_test, y_prob):.4f})逻辑说明stratifyy保证训练集和测试集的类别比例一致样本量小时这步很重要。hidden_layer_sizes(64, 32)是两层隐藏层节点数按输入特征10维来配别搞太大否则过拟合。参数说明max_iter500是迭代上限如果收敛警告出现说明需要增加。random_state固定后结果可复现论文里跑了100次测试实际使用时可以换不同随机种子看稳定性。注意论文里MLP在20%训练样本时准确率1.0这在小样本下要警惕过拟合。建议用交叉验证而不是单次划分来评估尤其是样本量只有100时。5. 避坑与排查复现这套流程时最容易翻车的5个点5.1 现象去噪后SNR很高但分类准确率反而下降原因过度去噪把故障冲击成分也滤掉了。SG滤波窗口设太大比如31以上时信号峰值会被平滑掉而峰值恰恰是区分A、B类的关键特征。解决去噪后重新算一遍峭度和脉冲指数如果这两个特征在去噪后区分度下降说明去噪过头了。把窗口长度调小到7或9再试。5.2 现象K-means聚类结果每次跑都不一样原因K-means对初始簇中心敏感n_init默认值在不同sklearn版本里不一样老版本是10新版本可能改成auto。解决显式设置n_init10或更高并固定random_state。如果还是不稳定说明特征空间里两类边界模糊需要回到特征提取那步检查。5.3 现象MLP训练时出现ConvergenceWarning原因max_iter不够或学习率不合适。论文里MLP在80%样本时准确率1.0但没写迭代次数实际跑的时候可能500次不够。解决先设max_iter1000如果还警告就调learning_rate_init从0.001降到0.0001。别直接上Adam小样本下SGD加动量往往更稳。5.4 现象一维信号reshape成64×64图像后可视化是一片黑原因原始数据量纲是[-1,1]直接乘以255后大部分像素值接近0或255对比度极低。解决论文里的做法是先平移缩放到[0,1]再乘255。但更稳的做法是做归一化(x - x.min()) / (x.max() - x.min()) * 255这样每张图都能充分利用灰度范围。5.5 现象迁移学习时AlexNet和ResNet-18的准确率远低于MLP原因预训练模型是在ImageNet上训的输入是3通道224×224而你的数据是单通道64×64。直接喂进去要么报错要么特征提取层学不到东西。解决把单通道复制成3通道用transforms.Resize(224)上采样。微调时只训最后的全连接层前面卷积层冻结。论文里用了SGD和二分类交叉熵batch_size16这些参数在小样本下是合理的。6. 进阶技巧一维信号转图像与迁移学习的参数调优把4096点一维信号reshape成64×64图像这个操作乍看只是换个数据形状实际上改变了模型看数据的方式。一维卷积关注的是局部时序模式二维卷积关注的是空间纹理。论文的可视化结果显示A类和B类在图像上差异明显——A类纹理更均匀B类有明显的亮斑聚集。这意味着图像分类模型能捕捉到一维模型忽略的空间结构信息。但迁移学习不是把预训练模型拿来就能用。我踩过的坑是直接微调全部层结果过拟合到训练集测试集准确率比MLP还低。后来改成只训最后三层前面冻结准确率才上来。论文里提到的预训练加微调策略具体操作时要注意学习率要设小因为预训练权重已经很好大学习率会破坏它们。import torch import torch.nn as nn from torchvision import models, transforms # 数据预处理单通道转3通道上采样到224 transform transforms.Compose([ transforms.Resize(224), transforms.Grayscale(num_output_channels3), transforms.ToTensor(), transforms.Normalize(mean[0.485, 0.456, 0.406], std[0.229, 0.224, 0.225]) ]) # 加载预训练ResNet-18 model models.resnet18(pretrainedTrue) # 冻结前面所有层 for param in model.parameters(): param.requires_grad False # 替换最后的全连接层为二分类 model.fc nn.Linear(model.fc.in_features, 2) # 只训练fc层 optimizer torch.optim.SGD(model.fc.parameters(), lr0.001, momentum0.9) criterion nn.CrossEntropyLoss()逻辑说明requires_grad False冻结卷积层只更新全连接层。model.fc是ResNet最后的分类头替换成2输出对应二分类。Normalize用的ImageNet均值和标准差因为预训练权重就是在这个分布上学的。参数说明lr0.001比从头训练的学习率小一个数量级因为微调不需要大更新。momentum0.9是SGD的标配。如果显存不够把batch_size从16降到8但别降到4以下否则BatchNorm统计量不准。论文在模型改进部分提到可以用Transformer做Patch序列建模这个思路在2022年还算前沿现在已经是标配了。如果要把这套方案用到实际故障诊断场景我的习惯是先跑通SG滤波10特征MLP这条baseline拿到一个可解释的准确率再尝试图像化迁移学习看能不能提升。如果baseline已经95%以上图像化带来的提升可能不值得增加的复杂度。从那以后我每次做信号类项目都强制先跑一遍特征提取传统模型的baseline再决定要不要上深度学习。希望帮到你。本文还有配套的精品资源点击获取