1. 项目背景与核心挑战去年参加天府杯数学建模竞赛的经历现在回想起来依然觉得收获颇丰。我们团队当时选的是A题关于仪器故障智能诊断。这个题目乍一看感觉像是传统工业领域的问题但组委会给的数据集和问题描述一下子就把我们拉到了数据科学和智能算法的前沿。题目要求我们基于给定的传感器时序数据构建一个能够自动、精准识别多种潜在故障模式的智能诊断系统。这不仅仅是套用一个现成的分类模型那么简单它涉及到信号处理、特征工程、模型选择、结果可解释性等一系列环环相扣的挑战。当时我们面临的核心痛点非常明确第一数据是典型的多维时间序列包含了振动、温度、压力等多种传感器在不同时间点的读数噪声大、维度高直接喂给模型效果肯定不好。第二故障模式并非独立发生早期故障信号极其微弱容易被噪声淹没如何从海量数据中提取出对故障敏感的特征是诊断准确率提升的关键。第三赛题不仅要求诊断出故障类型还希望我们对故障的严重程度或发展阶段做出评估这要求模型具备一定的回归或排序能力。最终我们团队通过一套融合了信号处理、传统机器学习与深度学习的混合策略成功解决了这些问题并拿到了一等奖。这篇文章我就把我们的解题思路、技术细节以及用Python实现的核心代码毫无保留地分享出来。无论你是正在备战数学建模竞赛还是对工业预测性维护、时序数据分析感兴趣相信都能从中获得直接的启发和可复用的代码。2. 解题总览从问题定义到技术路线图面对“仪器故障智能诊断”这样一个开放性问题第一步也是最关键的一步就是明确我们要解决的具体是什么问题。组委会提供的数据通常是一个包含多个csv文件的数据包每个文件可能对应一台设备、一段时间内的运行数据或者不同故障模式下的样本。列通常包括时间戳、若干传感器通道如acc_x,acc_y,temp,pressure等以及一个label或fault_type列在训练集中。我们的目标可以拆解为三个层次故障检测判断设备在某个时间窗口内是否发生了故障二分类正常 vs 异常。故障识别如果发生故障具体是哪种类型多分类如轴承内圈故障、外圈故障、齿轮磨损等。故障程度评估量化故障的严重性回归或有序分类如轻微、中等、严重。技术路线上我们没有押宝单一模型而是设计了一个分阶段的流水线Pipeline这样既能保证基础模型的稳健性又能利用深度模型挖掘深层特征。整体流程如下原始时序数据 - 数据预处理与清洗 - 时域/频域/时频域特征提取 - 特征选择 - (路径A)传统机器学习模型 - (路径B)深度学习模型 - 模型融合与决策 - 结果输出与可视化这个双路径设计是我们的核心策略。路径A传统机器学习依赖精心设计的特征工程模型如XGBoost、LightGBM解释性强训练快。路径B深度学习如1D-CNN、LSTM能自动学习特征对原始数据中的复杂模式捕捉能力更强。两者优势互补通过加权投票或堆叠Stacking方式融合最终诊断的鲁棒性和准确性显著提升。3. 数据预处理为模型提供“干净”的燃料原始工业传感器数据几乎不可能是完美无缺的。直接建模等于让模型在噪音中学习事倍功半。我们的预处理步骤主要解决以下四个问题3.1 缺失值与异常值处理传感器可能短暂失灵产生NaN或明显超出物理量程的异常值。import pandas as pd import numpy as np def handle_missing_and_outliers(df, sensor_columns): 处理缺失值和基于标准差/分位数的异常值。 df: 包含传感器数据的DataFrame sensor_columns: 传感器列名的列表 df_filled df.copy() # 1. 缺失值处理对于时间序列用前后时刻的均值填充更合理 for col in sensor_columns: df_filled[col] df_filled[col].interpolate(methodlinear) # 线性插值 # 如果开头或结尾还有NaN用最近的有效值填充 df_filled[col] df_filled[col].fillna(methodbfill).fillna(methodffill) # 2. 异常值处理使用基于IQR四分位距的方法 for col in sensor_columns: Q1 df_filled[col].quantile(0.25) Q3 df_filled[col].quantile(0.75) IQR Q3 - Q1 lower_bound Q1 - 1.5 * IQR upper_bound Q3 1.5 * IQR # 将超出边界的值替换为边界值或视为缺失值再填充 df_filled[col] np.where((df_filled[col] lower_bound) | (df_filled[col] upper_bound), np.nan, df_filled[col]) # 再次填充因异常值替换产生的NaN df_filled[col] df_filled[col].interpolate(methodlinear).fillna(methodbfill).fillna(methodffill) return df_filled注意对于高频振动信号简单的插值可能会引入虚假频率成分。在要求极高的场景下需要考虑更专业的信号处理方法如基于模型预测的插值。但在数学建模的有限时间内IQR插值是稳健且高效的选择。3.2 数据标准化与平滑不同传感器量纲和量级差异巨大例如加速度单位是g温度是摄氏度。必须进行标准化防止量级大的特征主导模型。我们通常使用StandardScaler减去均值除以标准差因为它能保留数据的分布形状对后续的PCA等线性变换友好。同时为了抑制高频噪声可以对信号进行滑动平均滤波。from sklearn.preprocessing import StandardScaler def normalize_and_smooth(df, sensor_columns, window_size5): 标准化并应用简单的移动平均平滑。 window_size: 滑动窗口大小需为奇数。 df_processed df.copy() scaler StandardScaler() # 先平滑再标准化顺序有时有影响可根据实验调整 for col in sensor_columns: # 滑动平均平滑 df_processed[col] df_processed[col].rolling(windowwindow_size, centerTrue, min_periods1).mean() # 标准化 df_processed[sensor_columns] scaler.fit_transform(df_processed[sensor_columns]) # 保存scaler用于后续的测试数据转换 return df_processed, scaler3.3 样本切片与标签对齐原始数据是长序列但我们需要将其切割成固定长度的时间窗口作为模型的一个个“样本”。这里的关键是标签对齐一个时间窗口对应一个故障标签。通常我们假设在一个短时间窗口内故障类型是稳定的。采用滑动窗口方法进行切片并可以设置重叠overlap以增加样本量。def create_samples(data_sequence, labels, window_size, step_size): 将长序列切割成重叠的时间窗口样本。 data_sequence: 形状为 (n_timesteps, n_features) 的传感器数据数组 labels: 形状为 (n_timesteps,) 的标签数组每个时间点一个标签 window_size: 窗口长度时间步数 step_size: 滑动步长 X, y [], [] n_samples len(data_sequence) for start in range(0, n_samples - window_size 1, step_size): end start window_size X.append(data_sequence[start:end]) # 取窗口内最主要的标签作为该样本的标签对于分类问题 window_labels labels[start:end] from scipy.stats import mode label, _ mode(window_labels, keepdimsFalse) y.append(label) return np.array(X), np.array(y)实操心得window_size和step_size是超参数。window_size要足够长以包含故障特征周期可通过分析故障频率初步估算但太长会混入不同状态的信息。step_size小于window_size会产生重叠样本能有效增加数据量防止切割时恰好切掉关键特征但也会引入样本相关性。我们通常设置重叠率为50%。4. 特征工程从原始信号中“榨取”信息这是传统机器学习路径路径A的成败关键。好的特征应该对故障敏感同时对工况变化如转速、负载相对鲁棒。我们从三个域进行特征提取4.1 时域特征直接从时间序列的幅值统计信息中提取计算简单物理意义明确。import numpy as np from scipy import stats def extract_time_domain_features(signal): 提取单个传感器通道在一个时间窗口内的时域特征。 features {} features[mean] np.mean(signal) features[std] np.std(signal) features[rms] np.sqrt(np.mean(signal**2)) # 均方根值反映能量 features[peak] np.max(np.abs(signal)) # 峰值 features[skewness] stats.skew(signal) # 偏度衡量分布不对称性 features[kurtosis] stats.kurtosis(signal) # 峰度衡量分布尖锐程度 features[crest_factor] features[peak] / features[rms] if features[rms] ! 0 else 0 # 峰值因子 features[clearance_factor] features[peak] / (np.mean(np.sqrt(np.abs(signal)))**2) if np.mean(np.sqrt(np.abs(signal))) ! 0 else 0 # 裕度因子 # 还可以增加波形因子、脉冲因子等 return features4.2 频域特征故障常常在振动信号的频谱中表现出特定的频率成分如轴承的故障特征频率。通过快速傅里叶变换FFT将信号转换到频域。from scipy.fft import fft, fftfreq def extract_freq_domain_features(signal, sampling_rate): 提取频域特征。 signal: 时间窗口信号 sampling_rate: 采样频率 (Hz) n len(signal) yf fft(signal) # 取单边频谱 yf_abs 2.0/n * np.abs(yf[:n//2]) xf fftfreq(n, 1/sampling_rate)[:n//2] features {} features[dominant_freq] xf[np.argmax(yf_abs)] # 主频 features[dominant_amp] np.max(yf_abs) # 主频幅值 # 计算频谱重心、均方频率、频率方差等 features[spectral_centroid] np.sum(xf * yf_abs) / np.sum(yf_abs) if np.sum(yf_abs) ! 0 else 0 features[spectral_rms] np.sqrt(np.sum((xf**2) * yf_abs) / np.sum(yf_abs)) if np.sum(yf_abs) ! 0 else 0 # 可以计算特定频带如故障特征频率附近的能量占比 return features4.3 时频域特征对于非平稳信号即统计特性随时间变化的信号单纯的频域分析会丢失时间信息。短时傅里叶变换STFT或小波变换能提供联合时频信息。我们常用小波包变换WPT因为它能对高频部分进行更精细的分解适合提取故障引起的瞬态冲击特征。import pywt # 需要安装PyWavelets def extract_wavelet_features(signal, waveletdb4, level3): 进行小波包分解并计算各节点子频带的能量作为特征。 wp pywt.WaveletPacket(datasignal, waveletwavelet, modesymmetric, maxlevellevel) # 获取第level层所有节点的名称如 aaa, aad, ada, ... nodes [node.path for node in wp.get_level(level, natural)] energy_features [] for node_name in nodes: node_coeffs wp[node_name].data node_energy np.sum(node_coeffs**2) energy_features.append(node_energy) # 通常将能量归一化构成能量分布向量 total_energy np.sum(energy_features) energy_features_norm [e/total_energy for e in energy_features] if total_energy ! 0 else energy_features return energy_features_norm将所有传感器通道、所有域的特征拼接起来会得到一个高维特征向量。接下来必须进行特征选择去除冗余和无关特征防止“维数灾难”。我们使用了基于树模型如XGBoost的特征重要性排序结合递归特征消除RFE来选择Top-N个最重要的特征。5. 模型构建双路径融合策略5.1 路径A基于特征工程的机器学习模型我们选择了LightGBM作为主力模型。它训练速度快对类别不平衡有一定处理能力并且能输出特征重要性与我们的特征工程流程完美契合。import lightgbm as lgb from sklearn.model_selection import train_test_split, StratifiedKFold from sklearn.metrics import accuracy_score, classification_report, confusion_matrix def train_lightgbm(X_features, y, paramsNone): X_features: 特征工程后得到的特征矩阵 (n_samples, n_features) y: 标签 if params is None: params { objective: multiclass, # 多分类 num_class: len(np.unique(y)), metric: multi_logloss, boosting_type: gbdt, num_leaves: 31, learning_rate: 0.05, feature_fraction: 0.9, bagging_fraction: 0.8, bagging_freq: 5, verbose: -1, seed: 42 } # 划分训练集和验证集 X_train, X_val, y_train, y_val train_test_split(X_features, y, test_size0.2, stratifyy, random_state42) # 创建Dataset train_data lgb.Dataset(X_train, labely_train) val_data lgb.Dataset(X_val, labely_val, referencetrain_data) # 训练使用早停法防止过拟合 model lgb.train(params, train_data, valid_sets[val_data], num_boost_round1000, callbacks[lgb.early_stopping(stopping_rounds50), lgb.log_evaluation(period100)]) # 验证集评估 y_pred model.predict(X_val, num_iterationmodel.best_iteration) y_pred_class np.argmax(y_pred, axis1) print(fValidation Accuracy: {accuracy_score(y_val, y_pred_class):.4f}) print(classification_report(y_val, y_pred_class)) # 可视化特征重要性 lgb.plot_importance(model, max_num_features20, figsize(10,6)) return model5.2 路径B基于原始信号的深度学习模型我们设计了一个结合1D-CNN和LSTM的混合网络。CNN擅长提取局部空间特征如振动信号中的冲击波形LSTM擅长捕捉时间依赖关系。模型直接输入标准化后的原始时序窗口数据(window_size, n_sensors)。import tensorflow as tf from tensorflow.keras import layers, models, callbacks def build_hybrid_cnn_lstm(input_shape, num_classes): 构建1D-CNN LSTM混合模型。 input_shape: (window_size, n_sensors) model models.Sequential([ # 第一部分1D-CNN 提取局部特征 layers.Input(shapeinput_shape), layers.Conv1D(filters64, kernel_size3, activationrelu, paddingsame), layers.BatchNormalization(), layers.MaxPooling1D(pool_size2), layers.Conv1D(filters128, kernel_size3, activationrelu, paddingsame), layers.BatchNormalization(), layers.MaxPooling1D(pool_size2), layers.Dropout(0.3), # 第二部分LSTM 捕捉时序依赖 # 将CNN输出的序列输入到LSTM。return_sequencesTrue表示输出每个时间步的状态。 layers.LSTM(units64, return_sequencesTrue), layers.Dropout(0.3), layers.LSTM(units32), layers.Dropout(0.3), # 第三部分全连接层分类 layers.Dense(units64, activationrelu), layers.Dense(unitsnum_classes, activationsoftmax) ]) model.compile(optimizertf.keras.optimizers.Adam(learning_rate0.001), losssparse_categorical_crossentropy, metrics[accuracy]) model.summary() return model # 训练深度学习模型 def train_deep_model(model, X_train_seq, y_train, X_val_seq, y_val, epochs50): X_train_seq: 原始序列样本形状 (n_samples, window_size, n_sensors) early_stopping callbacks.EarlyStopping(monitorval_loss, patience10, restore_best_weightsTrue) reduce_lr callbacks.ReduceLROnPlateau(monitorval_loss, factor0.5, patience5, min_lr1e-6) history model.fit(X_train_seq, y_train, validation_data(X_val_seq, y_val), epochsepochs, batch_size32, callbacks[early_stopping, reduce_lr], verbose1) return model, history踩坑实录直接训练这个混合网络很容易过拟合尤其是在数据量有限的情况下。我们采用了强力的正则化策略除了网络结构中的Dropout和BatchNorm还在数据上做了随机缩放、添加高斯噪声、时间轴轻微扭曲等数据增强显著提升了模型的泛化能力。另外LSTM层对输入数据的标准化非常敏感务必确保输入数据已标准化。5.3 模型融合112的策略我们采用了加权投票法进行融合。两个模型在验证集上的准确率作为其权重的基础。def weighted_ensemble_predict(model_lgb, model_dl, X_feat, X_seq, weightsNone): 加权投票融合。 model_lgb: LightGBM模型输入特征工程后的数据X_feat model_dl: 深度学习模型输入原始序列数据X_seq weights: 两个模型的权重列表如 [0.4, 0.6]。默认为None则根据验证集准确率自动计算。 proba_lgb model_lgb.predict(X_feat, num_iterationmodel_lgb.best_iteration) # 已经是概率形式 proba_dl model_dl.predict(X_seq) if weights is None: # 这里假设我们已经有了两个模型在某个验证集上的准确率 acc_lgb, acc_dl # 例如acc_lgb 0.92, acc_dl 0.94 acc_lgb, acc_dl 0.92, 0.94 total_acc acc_lgb acc_dl weights [acc_lgb/total_acc, acc_dl/total_acc] # 加权平均概率 weighted_proba weights[0] * proba_lgb weights[1] * proba_dl final_pred np.argmax(weighted_proba, axis1) return final_pred, weighted_proba融合后我们在测试集上的准确率比单一的最佳模型通常是深度学习模型提升了约1-2个百分点更重要的是对于某些单一模型容易混淆的故障类别融合模型的判断更加稳定。6. 故障严重程度评估与结果可视化对于故障程度评估我们将其建模为一个**有序分类Ordinal Regression**问题而不是简单的多分类或回归。因为“轻微”、“中等”、“严重”之间存在明确的顺序关系。我们使用了“序数逻辑回归”的思想将其转化为多个二分类问题例如模型1区分“无/轻微” vs “中等/严重”模型2区分“无/轻微/中等” vs “严重”或者直接使用支持有序分类的损失函数如CORAL损失函数在神经网络中的实现。结果可视化对于诊断系统的可解释性至关重要。我们主要做了以下几类图混淆矩阵热力图清晰展示模型在各类别上的混淆情况。特征重要性条形图从LightGBM模型获取告诉我们哪些传感器、哪些特征对诊断贡献最大这对于后续的传感器优化布置有指导意义。t-SNE/PCA降维图将高维特征或深度学习模型最后一层隐藏层的输出降到2维或3维进行可视化观察不同故障类别的样本在特征空间是否能够被良好区分。关键传感器信号对比图将正常状态和不同故障状态下的关键传感器如振动最大的那个原始信号或频谱图画在一起直观展示故障特征。import matplotlib.pyplot as plt import seaborn as sns from sklearn.manifold import TSNE def visualize_tsne(features, labels, titlet-SNE Visualization of Features): 使用t-SNE对高维特征进行降维可视化。 tsne TSNE(n_components2, random_state42, perplexity30) features_2d tsne.fit_transform(features) plt.figure(figsize(10,8)) scatter plt.scatter(features_2d[:,0], features_2d[:,1], clabels, cmaptab20, alpha0.7, s10) plt.colorbar(scatter) plt.title(title) plt.xlabel(t-SNE Component 1) plt.ylabel(t-SNE Component 2) plt.tight_layout() plt.show()7. 参赛总结与可复现性建议回顾整个项目拿到一等奖的关键在于系统性的问题拆解和务实的技术选型。我们没有追求最花哨的模型而是确保数据预处理、特征工程、基础模型训练每个环节都扎实可靠最后用融合策略提升天花板。有几个特别重要的点想分享关于数据数学建模竞赛给的数据往往“不完美”可能存在标签噪声、传感器漂移等问题。我们花了近三分之一的时间在数据探索和清洗上这是后续所有工作的基石。可视化每一类故障的典型信号波形和频谱能建立直观认识甚至能发现数据本身可能存在的问题。关于特征时域、频域、时频域特征各有千秋。对于周期性明显的故障如轴承频域特征非常有效对于瞬态冲击故障如齿轮断齿小波包能量特征可能更好。不要盲目堆砌特征一定要结合特征重要性分析进行筛选。关于模型LightGBM这类树模型对特征工程的质量要求高但训练快、调参相对简单、解释性强非常适合作为基线模型和提供特征重要性。深度学习模型潜力大但依赖大量数据和高超的调参技巧防止过拟合。双路径并行的策略让我们在有限时间内既能有一个稳健的保底方案又能冲击更高的性能。关于代码在竞赛中代码的可复现性和模块化至关重要。我们将整个流程封装成多个函数和类数据加载、预处理、特征提取、模型训练、评估可视化使得调整参数、更换模型、交叉验证变得非常方便。最终提交的论文中清晰的流程图和核心代码片段也是加分项。如果你想在自己的项目或未来的竞赛中复现这套方法我的建议是从理解数据开始画出数据分布听听“数据的声音”。先搭建一个简单的基线系统比如只用时域特征LightGBM确保整个Pipeline能跑通。迭代优化在此基础上逐步加入频域特征、尝试深度学习模型、调整融合策略。每次只改变一个变量评估其效果。重视验证策略使用分层K折交叉验证来更稳健地评估模型性能避免因为数据划分的偶然性导致过拟合。这个项目让我深刻体会到解决一个复杂的工程问题往往不是靠一个“银弹”算法而是靠对问题的深刻理解、扎实的基础工作以及将多种工具巧妙组合的系统性思维。希望这份详细的总结和代码能为你打开一扇门助你在智能诊断或相关的数据科学道路上走得更远。