1. 项目概述从数学建模竞赛到工业实践去年带队参加天府杯数学建模竞赛A题“仪器故障智能诊断技术”让我们团队拿了个一等奖。赛题本身很有意思它把一个典型的工业预测性维护问题包装成了一个数据驱动的建模任务。简单来说就是给你一堆仪器运行时的传感器时序数据比如振动、温度、压力其中有些数据对应正常状态有些则对应已知的故障类型要求你构建一个模型不仅能识别出仪器是否“病了”还得精准诊断出它得的是哪种“病”。这其实就是工业界常说的“故障诊断”或“健康管理”PHM的核心。赛后很多同学和业内的朋友都来问具体的思路和代码索性就把我们当时的解决方案、踩过的坑以及用Python实现的核心代码整理成这篇干货。无论你是想学习如何将数学建模竞赛的题目落地为代码还是对故障诊断这个方向感兴趣希望了解从数据到模型的全流程这篇文章都能给你一个清晰的参考。我们会避开复杂的数学公式推导重点讲清楚为什么这么设计、代码怎么写、以及实际调参时的心得。2. 解题核心思路与方案选型面对“仪器故障智能诊断”这类问题解题路径通常很清晰数据预处理 - 特征工程 - 模型构建 - 评估优化。但难点在于每个环节都有大量细节和选择不同的选择直接决定了最终成绩的上限。2.1 问题本质与建模路径拆解题目提供的通常是时间序列数据例如来自多个传感器的振动信号。故障诊断本质上是一个模式识别或分类问题。我们需要从看似杂乱无章的波形数据中提取出能够表征不同健康状态正常、故障A、故障B...的“指纹”特征然后用一个分类器去学习这些特征与状态标签之间的映射关系。我们的建模路径分为四步信号预处理原始传感器数据往往包含噪声、趋势项等干扰。第一步是“清洗”数据为后续分析打下基础。我们采用了去趋势、滤波如低通滤波去除高频噪声、标准化等方法。特征提取这是最关键的一步。直接从原始波形数据分类效果通常很差。我们需要计算一系列时域、频域、时频域特征。例如时域的特征有均值、方差、峭度、峰值因子等频域特征可以通过快速傅里叶变换FFT得到如重心频率、均方频率等对于非平稳信号小波变换能提供更丰富的时频联合信息。特征选择/降维提取的特征维度可能很高且存在冗余。直接喂给模型会导致计算量大且可能过拟合。我们使用了递归特征消除RFE结合随机森林进行特征重要性排序并采用主成分分析PCA进行降维在保留大部分信息的同时减少特征数量。分类模型构建与集成单一模型有时会陷入瓶颈。我们采用了模型集成的思路以梯度提升树如XGBoost/LightGBM作为主力模型因其对表格型特征处理能力强、精度高。同时用支持向量机SVM和多层感知机MLP构建了一个简单的 stacking 集成模型进一步提升泛化能力。注意竞赛中特征工程和模型集成往往是拉开差距的关键。一个精心设计的特征其提升效果可能远优于换一个更复杂的模型。2.2 为什么选择“特征工程集成学习”路线在有限的时间和算力下竞赛通常只有3-4天这条路线是最务实、最高效的。可解释性与物理意义时域、频域特征有明确的物理含义如峭度对冲击信号敏感这比端到端的深度学习黑箱模型更容易让评委理解和接受也便于在论文中阐述。对数据量要求相对宽松与需要海量数据的深度学习相比基于特征的传统机器学习方法在中等规模的数据集上就能取得很好效果非常适合竞赛场景。集成学习提升稳健性单一模型可能在某些类别的故障上表现不佳。通过集成多个差异化的基学习器如树模型、SVM、神经网络可以平滑误差提高最终诊断的稳定性和准确率这在面对未知的测试集时尤为重要。快速迭代验证特征工程和传统机器学习模型的训练/预测速度很快允许我们在短时间内尝试多种特征组合和模型参数快速验证想法。3. 核心模块Python代码实现与详解接下来我将分模块展示核心代码并附上详细的注释和操作意图说明。我们假设数据已经以CSV格式加载其中包含多列传感器数据sensor_1,sensor_2, ... 和一列标签label0表示正常1,2,...表示不同故障。3.1 数据预处理与特征提取模块这个模块负责将原始信号转化为特征矩阵。import numpy as np import pandas as pd from scipy import signal, stats from scipy.fft import fft from sklearn.preprocessing import StandardScaler import pywt # 需要安装PyWavelets class FeatureExtractor: 信号特征提取器 提取时域、频域及时频域特征 def __init__(self, sampling_freq25600): self.sampling_freq sampling_freq # 采样频率需根据实际数据设置 def extract_time_domain_features(self, x): 提取时域统计特征 features {} x np.array(x).flatten() features[mean] np.mean(x) features[std] np.std(x) features[peak] np.max(np.abs(x)) features[rms] np.sqrt(np.mean(x**2)) # 均方根值表征能量 features[skewness] stats.skew(x) # 偏度衡量分布不对称性 features[kurtosis] stats.kurtosis(x) # 峭度对冲击信号敏感是故障诊断黄金特征 features[crest_factor] features[peak] / features[rms] if features[rms]!0 else 0 # 峰值因子 features[shape_factor] features[rms] / np.mean(np.abs(x)) if np.mean(np.abs(x))!0 else 0 # 波形因子 # 更多特征如脉冲因子、裕度因子等可根据需要添加 return features def extract_freq_domain_features(self, x): 提取频域特征 features {} n len(x) # 计算FFT yf fft(x) # 取单边频谱 yf_abs 2.0/n * np.abs(yf[:n//2]) freqs np.fft.fftfreq(n, 1/self.sampling_freq)[:n//2] # 计算频谱重心、均方频率等 features[freq_center] np.sum(freqs * yf_abs) / np.sum(yf_abs) if np.sum(yf_abs)!0 else 0 features[freq_rms] np.sqrt(np.sum((freqs**2) * yf_abs) / np.sum(yf_abs)) if np.sum(yf_abs)!0 else 0 # 寻找主频幅值最大的频率成分 if len(yf_abs) 0: features[main_freq] freqs[np.argmax(yf_abs)] else: features[main_freq] 0 return features def extract_wavelet_features(self, x, waveletdb4, level3): 提取小波包分解能量特征一种时频域特征 features {} # 进行小波包分解 wp pywt.WaveletPacket(datax, waveletwavelet, modesymmetric, maxlevellevel) # 获取最后一层所有节点的能量 nodes [node.path for node in wp.get_level(level, natural)] energy_list [] for node_path in nodes: node wp[node_path] coeff node.data energy np.sum(coeff**2) energy_list.append(energy) total_energy np.sum(energy_list) # 将能量占比作为特征 for i, energy in enumerate(energy_list): features[fwavelet_energy_ratio_{i}] energy / total_energy if total_energy ! 0 else 0 return features def extract_all_features(self, signal_data): 提取所有特征 all_features {} # 时域 all_features.update(self.extract_time_domain_features(signal_data)) # 频域 all_features.update(self.extract_freq_domain_features(signal_data)) # 时频域小波 all_features.update(self.extract_wavelet_features(signal_data)) return all_features # 使用示例 def process_raw_data(df, sensor_columns, window_length1024, overlap0.5): 将原始长序列数据通过滑动窗口分割并提取每个窗口的特征 df: 原始DataFrame sensor_columns: 传感器列名列表 window_length: 窗口长度 overlap: 重叠率0.5表示50%重叠 feature_extractor FeatureExtractor(sampling_freq25600) # 根据实际数据设置 all_samples_features [] all_samples_labels [] step int(window_length * (1 - overlap)) for idx, row in df.iterrows(): label row[label] for col in sensor_columns: signal_series row[col] # 假设每个传感器的数据是一个长序列这里简化为处理一个数组 # 实际中可能需要先确保数据是等长的序列 if isinstance(signal_series, (list, np.ndarray)): signal_len len(signal_series) for start in range(0, signal_len - window_length 1, step): end start window_length window_data signal_series[start:end] # 对每个窗口提取特征 features feature_extractor.extract_all_features(window_data) features[sensor] col # 可加入传感器标识作为特征或后续处理依据 all_samples_features.append(features) all_samples_labels.append(label) else: # 如果数据已经是标量特征则直接使用需调整流程 pass # 转换为DataFrame features_df pd.DataFrame(all_samples_features) labels_series pd.Series(all_samples_labels) return features_df, labels_series操作意图与关键点解析滑动窗口仪器数据是连续的故障可能在任何时刻发生。滑动窗口将长序列切分成多个短样本既能增加训练数据量也能让模型学习局部模式。50%的重叠是常用设置能避免丢失跨窗口的故障信息。峭度Kurtosis这是故障诊断尤其是轴承、齿轮故障诊断中的明星特征。正常信号通常近似高斯分布峭度接近3。当出现局部损伤如点蚀、裂纹产生周期性冲击时信号分布会变得“尖峰厚尾”峭度值显著增大。小波能量特征傅里叶变换只告诉我们信号包含哪些频率但不知道这些频率何时出现。小波变换能同时提供时间和频率信息。我们通过计算小波包分解后各频带能量的占比来刻画信号能量在不同时频区域的分布变化这对非平稳故障信号非常有效。特征标准化在将所有特征合并成DataFrame后必须进行标准化如Z-score标准化消除不同特征量纲和数量级的影响。这一步通常在特征提取后模型训练前进行。from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from sklearn.feature_selection import RFECV from sklearn.ensemble import RandomForestClassifier def feature_engineering_pipeline(features_df, labels_series, n_components0.95): 特征后处理流水线标准化、特征选择、降维 # 1. 处理缺失值如果有 features_df.fillna(features_df.mean(), inplaceTrue) # 2. 标准化 scaler StandardScaler() features_scaled scaler.fit_transform(features_df) # 3. 递归特征消除 (RFE) 进行特征选择 # 使用一个计算效率较高的基模型进行初步筛选 estimator RandomForestClassifier(n_estimators50, random_state42, n_jobs-1) selector RFECV(estimator, step10, cv5, scoringaccuracy, n_jobs-1, min_features_to_select20) selector.fit(features_scaled, labels_series) features_selected features_scaled[:, selector.support_] print(f原始特征数: {features_scaled.shape[1]}, 筛选后特征数: {features_selected.shape[1]}) # 4. 主成分分析 (PCA) 降维 pca PCA(n_componentsn_components, random_state42) # 保留95%的方差 features_pca pca.fit_transform(features_selected) print(fPCA降维后特征数: {features_pca.shape[1]}, 累计方差贡献率: {sum(pca.explained_variance_ratio_):.4f}) return features_pca, labels_series.values, scaler, selector, pca3.2 模型构建、训练与集成模块特征准备好后就是构建模型。我们采用XGBoost作为主力并结合SVM和MLP进行Stacking集成。import xgboost as xgb from sklearn.svm import SVC from sklearn.neural_network import MLPClassifier from sklearn.ensemble import StackingClassifier from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold from sklearn.metrics import classification_report, confusion_matrix, accuracy_score import warnings warnings.filterwarnings(ignore) def build_and_train_models(X, y, test_size0.2, random_state42): 构建并训练单一模型及集成模型 # 划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split( X, y, test_sizetest_size, random_staterandom_state, stratifyy ) # 1. 定义基学习器 # XGBoost 参数是调优后的结果后面会讲调参过程 xgb_clf xgb.XGBClassifier( n_estimators200, max_depth6, learning_rate0.05, subsample0.8, colsample_bytree0.8, objectivemulti:softprob, random_staterandom_state, n_jobs-1, use_label_encoderFalse, eval_metricmlogloss ) # SVM 对特征尺度敏感我们的数据已经标准化了 svm_clf SVC(kernelrbf, C10, gammascale, probabilityTrue, random_staterandom_state) # MLP 神经网络 mlp_clf MLPClassifier( hidden_layer_sizes(128, 64), activationrelu, solveradam, alpha0.0001, batch_size256, learning_rateadaptive, max_iter500, random_staterandom_state, early_stoppingTrue ) # 2. 训练并评估单一模型 print( 单一模型性能 ) models {XGBoost: xgb_clf, SVM: svm_clf, MLP: mlp_clf} for name, model in models.items(): model.fit(X_train, y_train) y_pred model.predict(X_test) acc accuracy_score(y_test, y_pred) print(f{name} 测试集准确率: {acc:.4f}) # 可以详细打印分类报告 # print(classification_report(y_test, y_pred)) # 3. 构建Stacking集成模型 # 第一层基学习器 estimators [ (xgb, xgb_clf), (svm, svm_clf), (mlp, mlp_clf) ] # 第二层元学习器使用简单的逻辑回归 stacking_clf StackingClassifier( estimatorsestimators, final_estimatorLogisticRegression(max_iter1000, random_staterandom_state, n_jobs-1), cv5, # 使用5折交叉验证产生元特征 n_jobs-1 ) print(\n 训练Stacking集成模型 ) stacking_clf.fit(X_train, y_train) y_pred_stack stacking_clf.predict(X_test) acc_stack accuracy_score(y_test, y_pred_stack) print(fStacking集成模型 测试集准确率: {acc_stack:.4f}) print(\n 集成模型分类报告 ) print(classification_report(y_test, y_pred_stack)) # 返回训练好的模型和测试集用于后续分析 return { models: models, stacking_model: stacking_clf, X_test: X_test, y_test: y_test, X_train: X_train, y_train: y_train }实操心得Stacking的CV策略StackingClassifier中的cv参数至关重要。它指定了如何生成元特征即第一层模型的预测概率。这里设置cv5意味着使用5折交叉验证对于训练集的每一折用其他四折训练基模型并在本折上预测这样得到的预测值作为元特征可以防止信息泄露比直接用全量数据训练基模型再预测更稳健。元学习器的选择第二层模型元学习器不宜太复杂。因为第一层模型已经提供了很强的特征即各类别的预测概率一个简单的线性模型如逻辑回归往往就能很好地完成融合任务且不容易过拟合。类别不平衡处理如果故障样本数量差异很大需要在模型层面进行处理。对于XGBoost可以设置scale_pos_weight参数更通用的做法是在训练时使用class_weightbalanced如果模型支持或使用过采样/欠采样技术如SMOTE。我们在比赛中发现对少数故障类别进行适度的过采样能有效提升其召回率。4. 模型调优与超参数搜索实战模型框架搭好了但默认参数往往不是最优的。调参是提升模型性能的最后一道也是效果最显著的工序之一。我们以XGBoost为例展示如何使用网格搜索GridSearchCV和贝叶斯优化BayesianOptimization进行调参。4.1 网格搜索调参示例from sklearn.model_selection import GridSearchCV def tune_xgb_with_gridsearch(X_train, y_train): 使用网格搜索优化XGBoost参数 # 定义基础模型 xgb_model xgb.XGBClassifier(objectivemulti:softprob, random_state42, use_label_encoderFalse, n_jobs-1) # 定义参数网格 param_grid { n_estimators: [100, 200, 300], max_depth: [3, 6, 9], learning_rate: [0.01, 0.05, 0.1], subsample: [0.7, 0.8, 0.9], colsample_bytree: [0.7, 0.8, 0.9], gamma: [0, 0.1, 0.2] # 控制节点分裂所需的最小损失下降 } # 创建GridSearchCV对象 # 使用分层K折交叉验证确保每折类别分布一致 kfold StratifiedKFold(n_splits5, shuffleTrue, random_state42) grid_search GridSearchCV( estimatorxgb_model, param_gridparam_grid, scoringaccuracy, cvkfold, verbose1, n_jobs-1 ) # 执行搜索 print(开始网格搜索...) grid_search.fit(X_train, y_train) # 输出最佳参数和分数 print(f最佳参数: {grid_search.best_params_}) print(f最佳交叉验证分数: {grid_search.best_score_:.4f}) return grid_search.best_estimator_4.2 贝叶斯优化调参更高效网格搜索在参数组合多时计算量爆炸。贝叶斯优化通过建立目标函数模型性能与超参数之间的概率模型智能地选择下一个待评估的参数点能用更少的迭代找到更优解。from bayes_opt import BayesianOptimization from sklearn.model_selection import cross_val_score def xgb_cv(n_estimators, max_depth, lr, subsample, colsample_bytree, gamma, X, y): 为贝叶斯优化定义的交叉验证目标函数 # 将传入的浮点参数转换为合适的类型 params { n_estimators: int(n_estimators), max_depth: int(max_depth), learning_rate: max(lr, 0.01), # 确保学习率为正 subsample: max(min(subsample, 0.99), 0.1), # 限制在合理范围 colsample_bytree: max(min(colsample_bytree, 0.99), 0.1), gamma: max(gamma, 0), objective: multi:softprob, random_state: 42, use_label_encoder: False, eval_metric: mlogloss, n_jobs: -1 } model xgb.XGBClassifier(**params) # 使用3折交叉验证评估 cv_scores cross_val_score(model, X, y, cvStratifiedKFold(3, shuffleTrue, random_state42), scoringaccuracy, n_jobs-1) return cv_scores.mean() # 返回平均准确率 def tune_xgb_with_bayesian(X_train, y_train, init_points5, n_iter20): 使用贝叶斯优化优化XGBoost参数 # 定义参数边界 pbounds { n_estimators: (50, 400), max_depth: (3, 10), lr: (0.005, 0.2), # 学习率 subsample: (0.5, 1.0), colsample_bytree: (0.5, 1.0), gamma: (0, 0.5) } # 创建优化器 optimizer BayesianOptimization( flambda n_estimators, max_depth, lr, subsample, colsample_bytree, gamma: xgb_cv(n_estimators, max_depth, lr, subsample, colsample_bytree, gamma, X_train, y_train), pboundspbounds, random_state42, verbose2 ) # 执行优化 print(开始贝叶斯优化...) optimizer.maximize(init_pointsinit_points, n_itern_iter) # 获取最佳参数 best_params optimizer.max[params] # 对需要整数的参数进行处理 best_params[n_estimators] int(best_params[n_estimators]) best_params[max_depth] int(best_params[max_depth]) best_params[learning_rate] best_params.pop(lr) # 键名转换 print(f贝叶斯优化最佳参数: {best_params}) print(f最佳目标值 (CV Accuracy): {optimizer.max[target]:.4f}) # 用最佳参数训练最终模型 final_model xgb.XGBClassifier(**best_params, objectivemulti:softprob, random_state42, use_label_encoderFalse, n_jobs-1) final_model.fit(X_train, y_train) return final_model调参经验与避坑指南先粗后精不要一开始就在很细的网格上搜索。先用大范围、少步长确定大致最优区间如learning_rate: [0.01, 0.1, 0.3]再在最优值附近进行精细搜索。理解参数含义learning_rate(eta)学习率控制每棵树的权重。越小越稳健但需要更多的树n_estimators来达到相同效果。通常设为0.01-0.2。max_depth树的最大深度控制模型复杂度。越深越容易过拟合通常3-10之间。subsample和colsample_bytree行采样和列采样比例类似于随机森林用于增强模型多样性和防止过拟合。gamma节点分裂所需的最小损失下降。值越大算法越保守。使用交叉验证绝对不要用测试集来调参必须使用训练集内部的交叉验证分数来评估参数好坏否则会严重高估模型性能数据泄露。贝叶斯优化 vs 网格搜索对于超过3个参数的情况强烈推荐贝叶斯优化。它通常能用网格搜索1/10甚至更少的模型训练次数找到同等或更优的参数组合极大节省时间。早停法Early Stopping在训练XGBoost时可以设置early_stopping_rounds当验证集性能在若干轮内不再提升时停止训练防止过拟合并节省时间。这在n_estimators设得较大时特别有用。5. 故障诊断结果可视化与模型解释模型训练好后不能只看一个准确率数字就完事。我们需要深入分析模型在哪里做对了在哪里做错了以及它做出判断的依据是什么。这对于竞赛论文的撰写和实际工业应用都至关重要。5.1 混淆矩阵与分类报告可视化import matplotlib.pyplot as plt import seaborn as sns from sklearn.metrics import confusion_matrix, classification_report def evaluate_and_visualize(model, X_test, y_test, class_namesNone): 评估模型并可视化混淆矩阵和特征重要性 y_pred model.predict(X_test) y_pred_proba model.predict_proba(X_test) if hasattr(model, predict_proba) else None # 1. 混淆矩阵 cm confusion_matrix(y_test, y_pred) plt.figure(figsize(10, 8)) sns.heatmap(cm, annotTrue, fmtd, cmapBlues, xticklabelsclass_names if class_names else np.unique(y_test), yticklabelsclass_names if class_names else np.unique(y_test)) plt.title(Confusion Matrix) plt.ylabel(True Label) plt.xlabel(Predicted Label) plt.tight_layout() plt.show() # 2. 打印详细分类报告 print(Detailed Classification Report:) print(classification_report(y_test, y_pred, target_namesclass_names)) # 3. 特征重要性 (如果是树模型) if hasattr(model, feature_importances_): importances model.feature_importances_ indices np.argsort(importances)[::-1] plt.figure(figsize(12, 6)) plt.title(Feature Importances) plt.bar(range(20), importances[indices[:20]], aligncenter) # 显示前20个重要特征 plt.xticks(range(20), indices[:20], rotation90) plt.xlim([-1, 20]) plt.tight_layout() plt.show() # 可以关联回原始特征名分析哪些特征最有用 # top_feature_names [feature_names[i] for i in indices[:20]] # 4. 对于集成模型可以查看各基学习器的贡献以Stacking为例 if hasattr(model, named_estimators_): print(\n Stacking基学习器在测试集上的独立表现 ) for name, estimator in model.named_estimators_.items(): if hasattr(estimator, predict): acc accuracy_score(y_test, estimator.predict(X_test)) print(f{name}: {acc:.4f})5.2 使用SHAP进行模型解释高级对于重要的项目或论文我们还需要知道模型为什么做出某个预测。SHAP是一种强大的模型解释工具可以量化每个特征对于单个预测结果的贡献。import shap # 注意SHAP计算可能较慢建议对测试集的一个子集进行解释 def explain_model_with_shap(model, X_sample, feature_namesNone): 使用SHAP解释模型预测 X_sample: 用于解释的样本数据通常是测试集的一部分 # 创建SHAP解释器 if isinstance(model, xgb.XGBClassifier): explainer shap.TreeExplainer(model) elif isinstance(model, (SVC, MLPClassifier, LogisticRegression)): # 对于线性/核模型使用KernelExplainer或LinearExplainer # 注意KernelExplainer非常慢 explainer shap.KernelExplainer(model.predict_proba, X_sample[:100]) # 使用少量背景数据 else: print(fModel type {type(model)} not fully supported for quick SHAP explanation.) return # 计算SHAP值 shap_values explainer.shap_values(X_sample) # 可视化摘要图 shap.summary_plot(shap_values, X_sample, feature_namesfeature_names, plot_typebar) shap.summary_plot(shap_values, X_sample, feature_namesfeature_names) # 可视化单个预测的解释 # 例如解释测试集第一个样本 sample_idx 0 shap.force_plot(explainer.expected_value, shap_values[sample_idx], X_sample[sample_idx], feature_namesfeature_names)结果分析要点混淆矩阵能清晰看出模型容易混淆哪些故障类型。例如如果故障A和故障B经常被互相误判可能说明我们提取的特征对这两种故障的区分度不够需要设计更具判别性的特征。分类报告关注精确率Precision、召回率Recall和F1-score。在故障诊断中召回率往往更重要——我们宁愿误报将正常判为故障也不愿漏报将故障判为正常。对于少数类故障要特别关注其召回率。特征重要性如果发现时域峭度、某频带小波能量等特征排名靠前这与我们的领域知识冲击信号导致峭度升高是吻合的可以作为论文中的有力论据。SHAP值它能告诉我们对于某一个具体的故障样本是哪些特征值例如极高的峭度、某个频率成分的能量突增导致了模型将其判定为当前故障类型。这种局部可解释性在向领域专家解释模型决策时极具说服力。6. 项目部署与持续改进思路竞赛结束拿到一等奖不是终点。如何将这套方案转化为一个可用的故障诊断系统这里分享一些简单的工程化思路。6.1 构建简易实时诊断流水线我们可以将训练好的模型、特征提取器、标准化器、PCA转换器等打包成一个Pipeline对象并保存为文件便于部署。import joblib from sklearn.pipeline import Pipeline def build_and_save_pipeline(best_model, scaler, selector, pca, feature_extractor_params): 构建并保存完整的诊断流水线 # 注意由于特征提取器是我们自定义的类需要将其关键参数保存并在部署时重新实例化 # 这里我们保存特征提取器的配置参数 pipeline_components { best_model: best_model, scaler: scaler, feature_selector: selector, pca: pca, feature_extractor_params: feature_extractor_params # 保存采样频率等参数 } # 保存到文件 joblib.dump(pipeline_components, fault_diagnosis_pipeline.pkl) print(Pipeline saved to fault_diagnosis_pipeline.pkl) def load_and_predict(pipeline_path, new_raw_signal, sensor_namevibration): 加载流水线并对新数据进行预测 new_raw_signal: 一维数组新的传感器读数 # 加载 components joblib.load(pipeline_path) model components[best_model] scaler components[scaler] selector components[feature_selector] pca components[pca] fe_params components[feature_extractor_params] # 1. 特征提取 (模拟线上过程) fe FeatureExtractor(**fe_params) # 假设新信号也需要滑动窗口处理这里简化为对整个信号提取特征或取最后一个窗口 # 实际部署中需要与训练时一致的窗口化和特征提取逻辑 features_dict fe.extract_all_features(new_raw_signal[-1024:]) # 取最后1024个点 features_vec np.array(list(features_dict.values())).reshape(1, -1) # 2. 特征预处理 (与训练时一致) features_scaled scaler.transform(features_vec) features_selected features_scaled[:, selector.support_] features_final pca.transform(features_selected) # 3. 预测 prediction model.predict(features_final) prediction_proba model.predict_proba(features_final) return prediction[0], prediction_proba[0] # 模拟使用 # best_model stacking_clf # 假设这是我们训练好的最佳模型 # build_and_save_pipeline(best_model, scaler, selector, pca, {sampling_freq: 25600}) # ... 在另一个地方 ... # pred_label, pred_proba load_and_predict(fault_diagnosis_pipeline.pkl, new_signal_data)6.2 持续改进与进阶方向深度学习探索对于数据量足够大的场景可以尝试端到端的深度学习模型如1D卷积神经网络1D-CNN或长短时记忆网络LSTM。CNN能自动学习信号中的局部特征LSTM擅长捕捉时序依赖关系。可以将原始信号直接输入网络省去复杂的手工特征工程。我们的经验是在数据量大的情况下深度学习模型上限可能更高但可解释性较差。迁移学习如果目标设备的故障数据很少可以考虑使用其他类似设备或公开数据集如著名的CWRU轴承数据集上预训练的模型进行微调Fine-tuning。这在小样本学习场景下非常有效。在线学习与模型更新工业现场的数据是持续产生的。可以设计一个在线学习系统当有新的、经过专家确认的故障数据到来时能够以较小的代价更新模型使其适应设备可能的缓慢退化或工况变化。不确定性量化模型给出“故障A”的预测时我们还应知道它有多“确信”。除了预测概率可以使用蒙特卡洛Dropout对于深度学习模型或集成模型预测的方差来量化预测的不确定性。对于低置信度的预测可以触发人工复核提高系统可靠性。结合机理模型纯粹的数据驱动模型有时会做出违反物理规律的预测。将领域知识如设备的物理方程、故障演化机理融入模型形成物理信息神经网络PINN或采用知识图谱约束模型的输出是当前的研究前沿能显著提升模型的泛化能力和可信度。在整个项目实现过程中最大的体会是故障诊断是一个交叉领域既需要扎实的机器学习功底也需要对诊断对象仪器、设备的物理背景有基本了解。最好的特征往往来自于对信号和故障机理的深刻洞察而不是盲目地堆砌算法。我们的代码提供了一个完整、可复现的框架但当你面对具体问题时一定要花时间去分析数据、理解业务这样才能设计出最有针对性的解决方案而不是仅仅当一个“调参侠”。 基于特征工程与集成学习的工业故障诊断：从数学建模到Python实战
1. 项目概述从数学建模竞赛到工业实践去年带队参加天府杯数学建模竞赛A题“仪器故障智能诊断技术”让我们团队拿了个一等奖。赛题本身很有意思它把一个典型的工业预测性维护问题包装成了一个数据驱动的建模任务。简单来说就是给你一堆仪器运行时的传感器时序数据比如振动、温度、压力其中有些数据对应正常状态有些则对应已知的故障类型要求你构建一个模型不仅能识别出仪器是否“病了”还得精准诊断出它得的是哪种“病”。这其实就是工业界常说的“故障诊断”或“健康管理”PHM的核心。赛后很多同学和业内的朋友都来问具体的思路和代码索性就把我们当时的解决方案、踩过的坑以及用Python实现的核心代码整理成这篇干货。无论你是想学习如何将数学建模竞赛的题目落地为代码还是对故障诊断这个方向感兴趣希望了解从数据到模型的全流程这篇文章都能给你一个清晰的参考。我们会避开复杂的数学公式推导重点讲清楚为什么这么设计、代码怎么写、以及实际调参时的心得。2. 解题核心思路与方案选型面对“仪器故障智能诊断”这类问题解题路径通常很清晰数据预处理 - 特征工程 - 模型构建 - 评估优化。但难点在于每个环节都有大量细节和选择不同的选择直接决定了最终成绩的上限。2.1 问题本质与建模路径拆解题目提供的通常是时间序列数据例如来自多个传感器的振动信号。故障诊断本质上是一个模式识别或分类问题。我们需要从看似杂乱无章的波形数据中提取出能够表征不同健康状态正常、故障A、故障B...的“指纹”特征然后用一个分类器去学习这些特征与状态标签之间的映射关系。我们的建模路径分为四步信号预处理原始传感器数据往往包含噪声、趋势项等干扰。第一步是“清洗”数据为后续分析打下基础。我们采用了去趋势、滤波如低通滤波去除高频噪声、标准化等方法。特征提取这是最关键的一步。直接从原始波形数据分类效果通常很差。我们需要计算一系列时域、频域、时频域特征。例如时域的特征有均值、方差、峭度、峰值因子等频域特征可以通过快速傅里叶变换FFT得到如重心频率、均方频率等对于非平稳信号小波变换能提供更丰富的时频联合信息。特征选择/降维提取的特征维度可能很高且存在冗余。直接喂给模型会导致计算量大且可能过拟合。我们使用了递归特征消除RFE结合随机森林进行特征重要性排序并采用主成分分析PCA进行降维在保留大部分信息的同时减少特征数量。分类模型构建与集成单一模型有时会陷入瓶颈。我们采用了模型集成的思路以梯度提升树如XGBoost/LightGBM作为主力模型因其对表格型特征处理能力强、精度高。同时用支持向量机SVM和多层感知机MLP构建了一个简单的 stacking 集成模型进一步提升泛化能力。注意竞赛中特征工程和模型集成往往是拉开差距的关键。一个精心设计的特征其提升效果可能远优于换一个更复杂的模型。2.2 为什么选择“特征工程集成学习”路线在有限的时间和算力下竞赛通常只有3-4天这条路线是最务实、最高效的。可解释性与物理意义时域、频域特征有明确的物理含义如峭度对冲击信号敏感这比端到端的深度学习黑箱模型更容易让评委理解和接受也便于在论文中阐述。对数据量要求相对宽松与需要海量数据的深度学习相比基于特征的传统机器学习方法在中等规模的数据集上就能取得很好效果非常适合竞赛场景。集成学习提升稳健性单一模型可能在某些类别的故障上表现不佳。通过集成多个差异化的基学习器如树模型、SVM、神经网络可以平滑误差提高最终诊断的稳定性和准确率这在面对未知的测试集时尤为重要。快速迭代验证特征工程和传统机器学习模型的训练/预测速度很快允许我们在短时间内尝试多种特征组合和模型参数快速验证想法。3. 核心模块Python代码实现与详解接下来我将分模块展示核心代码并附上详细的注释和操作意图说明。我们假设数据已经以CSV格式加载其中包含多列传感器数据sensor_1,sensor_2, ... 和一列标签label0表示正常1,2,...表示不同故障。3.1 数据预处理与特征提取模块这个模块负责将原始信号转化为特征矩阵。import numpy as np import pandas as pd from scipy import signal, stats from scipy.fft import fft from sklearn.preprocessing import StandardScaler import pywt # 需要安装PyWavelets class FeatureExtractor: 信号特征提取器 提取时域、频域及时频域特征 def __init__(self, sampling_freq25600): self.sampling_freq sampling_freq # 采样频率需根据实际数据设置 def extract_time_domain_features(self, x): 提取时域统计特征 features {} x np.array(x).flatten() features[mean] np.mean(x) features[std] np.std(x) features[peak] np.max(np.abs(x)) features[rms] np.sqrt(np.mean(x**2)) # 均方根值表征能量 features[skewness] stats.skew(x) # 偏度衡量分布不对称性 features[kurtosis] stats.kurtosis(x) # 峭度对冲击信号敏感是故障诊断黄金特征 features[crest_factor] features[peak] / features[rms] if features[rms]!0 else 0 # 峰值因子 features[shape_factor] features[rms] / np.mean(np.abs(x)) if np.mean(np.abs(x))!0 else 0 # 波形因子 # 更多特征如脉冲因子、裕度因子等可根据需要添加 return features def extract_freq_domain_features(self, x): 提取频域特征 features {} n len(x) # 计算FFT yf fft(x) # 取单边频谱 yf_abs 2.0/n * np.abs(yf[:n//2]) freqs np.fft.fftfreq(n, 1/self.sampling_freq)[:n//2] # 计算频谱重心、均方频率等 features[freq_center] np.sum(freqs * yf_abs) / np.sum(yf_abs) if np.sum(yf_abs)!0 else 0 features[freq_rms] np.sqrt(np.sum((freqs**2) * yf_abs) / np.sum(yf_abs)) if np.sum(yf_abs)!0 else 0 # 寻找主频幅值最大的频率成分 if len(yf_abs) 0: features[main_freq] freqs[np.argmax(yf_abs)] else: features[main_freq] 0 return features def extract_wavelet_features(self, x, waveletdb4, level3): 提取小波包分解能量特征一种时频域特征 features {} # 进行小波包分解 wp pywt.WaveletPacket(datax, waveletwavelet, modesymmetric, maxlevellevel) # 获取最后一层所有节点的能量 nodes [node.path for node in wp.get_level(level, natural)] energy_list [] for node_path in nodes: node wp[node_path] coeff node.data energy np.sum(coeff**2) energy_list.append(energy) total_energy np.sum(energy_list) # 将能量占比作为特征 for i, energy in enumerate(energy_list): features[fwavelet_energy_ratio_{i}] energy / total_energy if total_energy ! 0 else 0 return features def extract_all_features(self, signal_data): 提取所有特征 all_features {} # 时域 all_features.update(self.extract_time_domain_features(signal_data)) # 频域 all_features.update(self.extract_freq_domain_features(signal_data)) # 时频域小波 all_features.update(self.extract_wavelet_features(signal_data)) return all_features # 使用示例 def process_raw_data(df, sensor_columns, window_length1024, overlap0.5): 将原始长序列数据通过滑动窗口分割并提取每个窗口的特征 df: 原始DataFrame sensor_columns: 传感器列名列表 window_length: 窗口长度 overlap: 重叠率0.5表示50%重叠 feature_extractor FeatureExtractor(sampling_freq25600) # 根据实际数据设置 all_samples_features [] all_samples_labels [] step int(window_length * (1 - overlap)) for idx, row in df.iterrows(): label row[label] for col in sensor_columns: signal_series row[col] # 假设每个传感器的数据是一个长序列这里简化为处理一个数组 # 实际中可能需要先确保数据是等长的序列 if isinstance(signal_series, (list, np.ndarray)): signal_len len(signal_series) for start in range(0, signal_len - window_length 1, step): end start window_length window_data signal_series[start:end] # 对每个窗口提取特征 features feature_extractor.extract_all_features(window_data) features[sensor] col # 可加入传感器标识作为特征或后续处理依据 all_samples_features.append(features) all_samples_labels.append(label) else: # 如果数据已经是标量特征则直接使用需调整流程 pass # 转换为DataFrame features_df pd.DataFrame(all_samples_features) labels_series pd.Series(all_samples_labels) return features_df, labels_series操作意图与关键点解析滑动窗口仪器数据是连续的故障可能在任何时刻发生。滑动窗口将长序列切分成多个短样本既能增加训练数据量也能让模型学习局部模式。50%的重叠是常用设置能避免丢失跨窗口的故障信息。峭度Kurtosis这是故障诊断尤其是轴承、齿轮故障诊断中的明星特征。正常信号通常近似高斯分布峭度接近3。当出现局部损伤如点蚀、裂纹产生周期性冲击时信号分布会变得“尖峰厚尾”峭度值显著增大。小波能量特征傅里叶变换只告诉我们信号包含哪些频率但不知道这些频率何时出现。小波变换能同时提供时间和频率信息。我们通过计算小波包分解后各频带能量的占比来刻画信号能量在不同时频区域的分布变化这对非平稳故障信号非常有效。特征标准化在将所有特征合并成DataFrame后必须进行标准化如Z-score标准化消除不同特征量纲和数量级的影响。这一步通常在特征提取后模型训练前进行。from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from sklearn.feature_selection import RFECV from sklearn.ensemble import RandomForestClassifier def feature_engineering_pipeline(features_df, labels_series, n_components0.95): 特征后处理流水线标准化、特征选择、降维 # 1. 处理缺失值如果有 features_df.fillna(features_df.mean(), inplaceTrue) # 2. 标准化 scaler StandardScaler() features_scaled scaler.fit_transform(features_df) # 3. 递归特征消除 (RFE) 进行特征选择 # 使用一个计算效率较高的基模型进行初步筛选 estimator RandomForestClassifier(n_estimators50, random_state42, n_jobs-1) selector RFECV(estimator, step10, cv5, scoringaccuracy, n_jobs-1, min_features_to_select20) selector.fit(features_scaled, labels_series) features_selected features_scaled[:, selector.support_] print(f原始特征数: {features_scaled.shape[1]}, 筛选后特征数: {features_selected.shape[1]}) # 4. 主成分分析 (PCA) 降维 pca PCA(n_componentsn_components, random_state42) # 保留95%的方差 features_pca pca.fit_transform(features_selected) print(fPCA降维后特征数: {features_pca.shape[1]}, 累计方差贡献率: {sum(pca.explained_variance_ratio_):.4f}) return features_pca, labels_series.values, scaler, selector, pca3.2 模型构建、训练与集成模块特征准备好后就是构建模型。我们采用XGBoost作为主力并结合SVM和MLP进行Stacking集成。import xgboost as xgb from sklearn.svm import SVC from sklearn.neural_network import MLPClassifier from sklearn.ensemble import StackingClassifier from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold from sklearn.metrics import classification_report, confusion_matrix, accuracy_score import warnings warnings.filterwarnings(ignore) def build_and_train_models(X, y, test_size0.2, random_state42): 构建并训练单一模型及集成模型 # 划分训练集和测试集 X_train, X_test, y_train, y_test train_test_split( X, y, test_sizetest_size, random_staterandom_state, stratifyy ) # 1. 定义基学习器 # XGBoost 参数是调优后的结果后面会讲调参过程 xgb_clf xgb.XGBClassifier( n_estimators200, max_depth6, learning_rate0.05, subsample0.8, colsample_bytree0.8, objectivemulti:softprob, random_staterandom_state, n_jobs-1, use_label_encoderFalse, eval_metricmlogloss ) # SVM 对特征尺度敏感我们的数据已经标准化了 svm_clf SVC(kernelrbf, C10, gammascale, probabilityTrue, random_staterandom_state) # MLP 神经网络 mlp_clf MLPClassifier( hidden_layer_sizes(128, 64), activationrelu, solveradam, alpha0.0001, batch_size256, learning_rateadaptive, max_iter500, random_staterandom_state, early_stoppingTrue ) # 2. 训练并评估单一模型 print( 单一模型性能 ) models {XGBoost: xgb_clf, SVM: svm_clf, MLP: mlp_clf} for name, model in models.items(): model.fit(X_train, y_train) y_pred model.predict(X_test) acc accuracy_score(y_test, y_pred) print(f{name} 测试集准确率: {acc:.4f}) # 可以详细打印分类报告 # print(classification_report(y_test, y_pred)) # 3. 构建Stacking集成模型 # 第一层基学习器 estimators [ (xgb, xgb_clf), (svm, svm_clf), (mlp, mlp_clf) ] # 第二层元学习器使用简单的逻辑回归 stacking_clf StackingClassifier( estimatorsestimators, final_estimatorLogisticRegression(max_iter1000, random_staterandom_state, n_jobs-1), cv5, # 使用5折交叉验证产生元特征 n_jobs-1 ) print(\n 训练Stacking集成模型 ) stacking_clf.fit(X_train, y_train) y_pred_stack stacking_clf.predict(X_test) acc_stack accuracy_score(y_test, y_pred_stack) print(fStacking集成模型 测试集准确率: {acc_stack:.4f}) print(\n 集成模型分类报告 ) print(classification_report(y_test, y_pred_stack)) # 返回训练好的模型和测试集用于后续分析 return { models: models, stacking_model: stacking_clf, X_test: X_test, y_test: y_test, X_train: X_train, y_train: y_train }实操心得Stacking的CV策略StackingClassifier中的cv参数至关重要。它指定了如何生成元特征即第一层模型的预测概率。这里设置cv5意味着使用5折交叉验证对于训练集的每一折用其他四折训练基模型并在本折上预测这样得到的预测值作为元特征可以防止信息泄露比直接用全量数据训练基模型再预测更稳健。元学习器的选择第二层模型元学习器不宜太复杂。因为第一层模型已经提供了很强的特征即各类别的预测概率一个简单的线性模型如逻辑回归往往就能很好地完成融合任务且不容易过拟合。类别不平衡处理如果故障样本数量差异很大需要在模型层面进行处理。对于XGBoost可以设置scale_pos_weight参数更通用的做法是在训练时使用class_weightbalanced如果模型支持或使用过采样/欠采样技术如SMOTE。我们在比赛中发现对少数故障类别进行适度的过采样能有效提升其召回率。4. 模型调优与超参数搜索实战模型框架搭好了但默认参数往往不是最优的。调参是提升模型性能的最后一道也是效果最显著的工序之一。我们以XGBoost为例展示如何使用网格搜索GridSearchCV和贝叶斯优化BayesianOptimization进行调参。4.1 网格搜索调参示例from sklearn.model_selection import GridSearchCV def tune_xgb_with_gridsearch(X_train, y_train): 使用网格搜索优化XGBoost参数 # 定义基础模型 xgb_model xgb.XGBClassifier(objectivemulti:softprob, random_state42, use_label_encoderFalse, n_jobs-1) # 定义参数网格 param_grid { n_estimators: [100, 200, 300], max_depth: [3, 6, 9], learning_rate: [0.01, 0.05, 0.1], subsample: [0.7, 0.8, 0.9], colsample_bytree: [0.7, 0.8, 0.9], gamma: [0, 0.1, 0.2] # 控制节点分裂所需的最小损失下降 } # 创建GridSearchCV对象 # 使用分层K折交叉验证确保每折类别分布一致 kfold StratifiedKFold(n_splits5, shuffleTrue, random_state42) grid_search GridSearchCV( estimatorxgb_model, param_gridparam_grid, scoringaccuracy, cvkfold, verbose1, n_jobs-1 ) # 执行搜索 print(开始网格搜索...) grid_search.fit(X_train, y_train) # 输出最佳参数和分数 print(f最佳参数: {grid_search.best_params_}) print(f最佳交叉验证分数: {grid_search.best_score_:.4f}) return grid_search.best_estimator_4.2 贝叶斯优化调参更高效网格搜索在参数组合多时计算量爆炸。贝叶斯优化通过建立目标函数模型性能与超参数之间的概率模型智能地选择下一个待评估的参数点能用更少的迭代找到更优解。from bayes_opt import BayesianOptimization from sklearn.model_selection import cross_val_score def xgb_cv(n_estimators, max_depth, lr, subsample, colsample_bytree, gamma, X, y): 为贝叶斯优化定义的交叉验证目标函数 # 将传入的浮点参数转换为合适的类型 params { n_estimators: int(n_estimators), max_depth: int(max_depth), learning_rate: max(lr, 0.01), # 确保学习率为正 subsample: max(min(subsample, 0.99), 0.1), # 限制在合理范围 colsample_bytree: max(min(colsample_bytree, 0.99), 0.1), gamma: max(gamma, 0), objective: multi:softprob, random_state: 42, use_label_encoder: False, eval_metric: mlogloss, n_jobs: -1 } model xgb.XGBClassifier(**params) # 使用3折交叉验证评估 cv_scores cross_val_score(model, X, y, cvStratifiedKFold(3, shuffleTrue, random_state42), scoringaccuracy, n_jobs-1) return cv_scores.mean() # 返回平均准确率 def tune_xgb_with_bayesian(X_train, y_train, init_points5, n_iter20): 使用贝叶斯优化优化XGBoost参数 # 定义参数边界 pbounds { n_estimators: (50, 400), max_depth: (3, 10), lr: (0.005, 0.2), # 学习率 subsample: (0.5, 1.0), colsample_bytree: (0.5, 1.0), gamma: (0, 0.5) } # 创建优化器 optimizer BayesianOptimization( flambda n_estimators, max_depth, lr, subsample, colsample_bytree, gamma: xgb_cv(n_estimators, max_depth, lr, subsample, colsample_bytree, gamma, X_train, y_train), pboundspbounds, random_state42, verbose2 ) # 执行优化 print(开始贝叶斯优化...) optimizer.maximize(init_pointsinit_points, n_itern_iter) # 获取最佳参数 best_params optimizer.max[params] # 对需要整数的参数进行处理 best_params[n_estimators] int(best_params[n_estimators]) best_params[max_depth] int(best_params[max_depth]) best_params[learning_rate] best_params.pop(lr) # 键名转换 print(f贝叶斯优化最佳参数: {best_params}) print(f最佳目标值 (CV Accuracy): {optimizer.max[target]:.4f}) # 用最佳参数训练最终模型 final_model xgb.XGBClassifier(**best_params, objectivemulti:softprob, random_state42, use_label_encoderFalse, n_jobs-1) final_model.fit(X_train, y_train) return final_model调参经验与避坑指南先粗后精不要一开始就在很细的网格上搜索。先用大范围、少步长确定大致最优区间如learning_rate: [0.01, 0.1, 0.3]再在最优值附近进行精细搜索。理解参数含义learning_rate(eta)学习率控制每棵树的权重。越小越稳健但需要更多的树n_estimators来达到相同效果。通常设为0.01-0.2。max_depth树的最大深度控制模型复杂度。越深越容易过拟合通常3-10之间。subsample和colsample_bytree行采样和列采样比例类似于随机森林用于增强模型多样性和防止过拟合。gamma节点分裂所需的最小损失下降。值越大算法越保守。使用交叉验证绝对不要用测试集来调参必须使用训练集内部的交叉验证分数来评估参数好坏否则会严重高估模型性能数据泄露。贝叶斯优化 vs 网格搜索对于超过3个参数的情况强烈推荐贝叶斯优化。它通常能用网格搜索1/10甚至更少的模型训练次数找到同等或更优的参数组合极大节省时间。早停法Early Stopping在训练XGBoost时可以设置early_stopping_rounds当验证集性能在若干轮内不再提升时停止训练防止过拟合并节省时间。这在n_estimators设得较大时特别有用。5. 故障诊断结果可视化与模型解释模型训练好后不能只看一个准确率数字就完事。我们需要深入分析模型在哪里做对了在哪里做错了以及它做出判断的依据是什么。这对于竞赛论文的撰写和实际工业应用都至关重要。5.1 混淆矩阵与分类报告可视化import matplotlib.pyplot as plt import seaborn as sns from sklearn.metrics import confusion_matrix, classification_report def evaluate_and_visualize(model, X_test, y_test, class_namesNone): 评估模型并可视化混淆矩阵和特征重要性 y_pred model.predict(X_test) y_pred_proba model.predict_proba(X_test) if hasattr(model, predict_proba) else None # 1. 混淆矩阵 cm confusion_matrix(y_test, y_pred) plt.figure(figsize(10, 8)) sns.heatmap(cm, annotTrue, fmtd, cmapBlues, xticklabelsclass_names if class_names else np.unique(y_test), yticklabelsclass_names if class_names else np.unique(y_test)) plt.title(Confusion Matrix) plt.ylabel(True Label) plt.xlabel(Predicted Label) plt.tight_layout() plt.show() # 2. 打印详细分类报告 print(Detailed Classification Report:) print(classification_report(y_test, y_pred, target_namesclass_names)) # 3. 特征重要性 (如果是树模型) if hasattr(model, feature_importances_): importances model.feature_importances_ indices np.argsort(importances)[::-1] plt.figure(figsize(12, 6)) plt.title(Feature Importances) plt.bar(range(20), importances[indices[:20]], aligncenter) # 显示前20个重要特征 plt.xticks(range(20), indices[:20], rotation90) plt.xlim([-1, 20]) plt.tight_layout() plt.show() # 可以关联回原始特征名分析哪些特征最有用 # top_feature_names [feature_names[i] for i in indices[:20]] # 4. 对于集成模型可以查看各基学习器的贡献以Stacking为例 if hasattr(model, named_estimators_): print(\n Stacking基学习器在测试集上的独立表现 ) for name, estimator in model.named_estimators_.items(): if hasattr(estimator, predict): acc accuracy_score(y_test, estimator.predict(X_test)) print(f{name}: {acc:.4f})5.2 使用SHAP进行模型解释高级对于重要的项目或论文我们还需要知道模型为什么做出某个预测。SHAP是一种强大的模型解释工具可以量化每个特征对于单个预测结果的贡献。import shap # 注意SHAP计算可能较慢建议对测试集的一个子集进行解释 def explain_model_with_shap(model, X_sample, feature_namesNone): 使用SHAP解释模型预测 X_sample: 用于解释的样本数据通常是测试集的一部分 # 创建SHAP解释器 if isinstance(model, xgb.XGBClassifier): explainer shap.TreeExplainer(model) elif isinstance(model, (SVC, MLPClassifier, LogisticRegression)): # 对于线性/核模型使用KernelExplainer或LinearExplainer # 注意KernelExplainer非常慢 explainer shap.KernelExplainer(model.predict_proba, X_sample[:100]) # 使用少量背景数据 else: print(fModel type {type(model)} not fully supported for quick SHAP explanation.) return # 计算SHAP值 shap_values explainer.shap_values(X_sample) # 可视化摘要图 shap.summary_plot(shap_values, X_sample, feature_namesfeature_names, plot_typebar) shap.summary_plot(shap_values, X_sample, feature_namesfeature_names) # 可视化单个预测的解释 # 例如解释测试集第一个样本 sample_idx 0 shap.force_plot(explainer.expected_value, shap_values[sample_idx], X_sample[sample_idx], feature_namesfeature_names)结果分析要点混淆矩阵能清晰看出模型容易混淆哪些故障类型。例如如果故障A和故障B经常被互相误判可能说明我们提取的特征对这两种故障的区分度不够需要设计更具判别性的特征。分类报告关注精确率Precision、召回率Recall和F1-score。在故障诊断中召回率往往更重要——我们宁愿误报将正常判为故障也不愿漏报将故障判为正常。对于少数类故障要特别关注其召回率。特征重要性如果发现时域峭度、某频带小波能量等特征排名靠前这与我们的领域知识冲击信号导致峭度升高是吻合的可以作为论文中的有力论据。SHAP值它能告诉我们对于某一个具体的故障样本是哪些特征值例如极高的峭度、某个频率成分的能量突增导致了模型将其判定为当前故障类型。这种局部可解释性在向领域专家解释模型决策时极具说服力。6. 项目部署与持续改进思路竞赛结束拿到一等奖不是终点。如何将这套方案转化为一个可用的故障诊断系统这里分享一些简单的工程化思路。6.1 构建简易实时诊断流水线我们可以将训练好的模型、特征提取器、标准化器、PCA转换器等打包成一个Pipeline对象并保存为文件便于部署。import joblib from sklearn.pipeline import Pipeline def build_and_save_pipeline(best_model, scaler, selector, pca, feature_extractor_params): 构建并保存完整的诊断流水线 # 注意由于特征提取器是我们自定义的类需要将其关键参数保存并在部署时重新实例化 # 这里我们保存特征提取器的配置参数 pipeline_components { best_model: best_model, scaler: scaler, feature_selector: selector, pca: pca, feature_extractor_params: feature_extractor_params # 保存采样频率等参数 } # 保存到文件 joblib.dump(pipeline_components, fault_diagnosis_pipeline.pkl) print(Pipeline saved to fault_diagnosis_pipeline.pkl) def load_and_predict(pipeline_path, new_raw_signal, sensor_namevibration): 加载流水线并对新数据进行预测 new_raw_signal: 一维数组新的传感器读数 # 加载 components joblib.load(pipeline_path) model components[best_model] scaler components[scaler] selector components[feature_selector] pca components[pca] fe_params components[feature_extractor_params] # 1. 特征提取 (模拟线上过程) fe FeatureExtractor(**fe_params) # 假设新信号也需要滑动窗口处理这里简化为对整个信号提取特征或取最后一个窗口 # 实际部署中需要与训练时一致的窗口化和特征提取逻辑 features_dict fe.extract_all_features(new_raw_signal[-1024:]) # 取最后1024个点 features_vec np.array(list(features_dict.values())).reshape(1, -1) # 2. 特征预处理 (与训练时一致) features_scaled scaler.transform(features_vec) features_selected features_scaled[:, selector.support_] features_final pca.transform(features_selected) # 3. 预测 prediction model.predict(features_final) prediction_proba model.predict_proba(features_final) return prediction[0], prediction_proba[0] # 模拟使用 # best_model stacking_clf # 假设这是我们训练好的最佳模型 # build_and_save_pipeline(best_model, scaler, selector, pca, {sampling_freq: 25600}) # ... 在另一个地方 ... # pred_label, pred_proba load_and_predict(fault_diagnosis_pipeline.pkl, new_signal_data)6.2 持续改进与进阶方向深度学习探索对于数据量足够大的场景可以尝试端到端的深度学习模型如1D卷积神经网络1D-CNN或长短时记忆网络LSTM。CNN能自动学习信号中的局部特征LSTM擅长捕捉时序依赖关系。可以将原始信号直接输入网络省去复杂的手工特征工程。我们的经验是在数据量大的情况下深度学习模型上限可能更高但可解释性较差。迁移学习如果目标设备的故障数据很少可以考虑使用其他类似设备或公开数据集如著名的CWRU轴承数据集上预训练的模型进行微调Fine-tuning。这在小样本学习场景下非常有效。在线学习与模型更新工业现场的数据是持续产生的。可以设计一个在线学习系统当有新的、经过专家确认的故障数据到来时能够以较小的代价更新模型使其适应设备可能的缓慢退化或工况变化。不确定性量化模型给出“故障A”的预测时我们还应知道它有多“确信”。除了预测概率可以使用蒙特卡洛Dropout对于深度学习模型或集成模型预测的方差来量化预测的不确定性。对于低置信度的预测可以触发人工复核提高系统可靠性。结合机理模型纯粹的数据驱动模型有时会做出违反物理规律的预测。将领域知识如设备的物理方程、故障演化机理融入模型形成物理信息神经网络PINN或采用知识图谱约束模型的输出是当前的研究前沿能显著提升模型的泛化能力和可信度。在整个项目实现过程中最大的体会是故障诊断是一个交叉领域既需要扎实的机器学习功底也需要对诊断对象仪器、设备的物理背景有基本了解。最好的特征往往来自于对信号和故障机理的深刻洞察而不是盲目地堆砌算法。我们的代码提供了一个完整、可复现的框架但当你面对具体问题时一定要花时间去分析数据、理解业务这样才能设计出最有针对性的解决方案而不是仅仅当一个“调参侠”。