
简介2024年五一数学建模竞赛C题获奖作品聚焦煤矿深部开采冲击地压危险预测问题完整呈现基于电磁辐射与声发射数据的数学建模全过程。论文综合运用傅里叶变换进行信号频域分析构建随机森林模型识别干扰信号所在时间区间并借助朴素贝叶斯模型计算前兆特征出现的概率。内容包含最大值、最小值、均值、标准偏差、峰度、偏度等统计特征的细致对比以及从问题重述、模型假设、数据预处理到结果分析的全流程逻辑。整份资料为PDF格式论文共1个文件压缩包约1.16MB已有1468人学习下载。适合参加数学建模竞赛的学生、从事煤矿动力灾害预测的工程技术人员阅读也可为机器学习在工业数据中的应用提供参考从中可获得完整建模思路、特征工程细节和论文写作规范。1. 冲击地压危险预测在五一赛C题里到底考什么时窗、漏报与现场可用性深部开采的冲击地压很多搞算法的人第一次听会以为是个地质灾害预测题一上来就堆模型。实际上把这题拆开看它考的是两件事一件是怎么把带时间戳的监测数据正确切成“用过去一段预测未来一段”的样本另一件是当危险样本很少的时候怎么让模型记住危险而不是牺牲准确率。前一件错了后面所有模型都白搭后一件没做好作品看起来分数很高现场却根本不敢用。这道赛题给的原始数据通常来自井下监测系统包含应力、微震能量、声发射计数、钻屑量等字段并且是按时间排列的。大多数第一次参赛的队伍会用处理静态表格的思路去做每个时刻的特征预测每个时刻的危险等级随机划分交叉验证模型很快就收敛到很高的分数。但这种方案在现场等于没有提前量。下面按“拆任务、造特征、选模型、避坑、验证”的顺序把一套能落地的建模方案完整拆开讲适合正在备赛的队也适合做矿压大数据预警的工程师参考。2. 先把赛题拆成能建模的样本目标时窗、滑窗切分与评估指标2.1 用过去一段预测未来一段才是“预警”赛题里的“预测”不是识别当前状态。对冲击地压来说报警必须发生在事故前。机械地用当前时刻的特征预测当前时刻风险会把问题退化成相关性问题不能给现场留出处置时间。我一般拿到数据后会先画两条时间窗特征窗口是t之前的一段时间目标窗口是t之后的一段时间。如果一条记录代表一个班次比较自然的方式是特征窗口取过去12个班次目标窗口取未来4个班次中间特意留出间隔。这样构造出来的是有监督学习样本过去的数据是X未来的风险状态是y。import pandas as pd import numpy as np def make_samples(df, feature_steps12, gap_steps1, target_steps4, step1, risk_level2): 把原始监测时间序列切成监督学习样本。 df: 按时间排序的监测数据必须包含 risk 列和所需传感器字段 feature_steps: 特征窗口覆盖多少个时间点 gap_steps : 特征窗口结束到目标窗口开始之间的间隔 target_steps: 目标窗口覆盖多少个时间点 step : 每隔多少个时间点取一个样本控制样本重叠 risk_level : 当目标窗口内最大风险等级大于等于该值时视为正样本 rows [] for i in range(0, len(df) - feature_steps - gap_steps - target_steps, step): feature_part df.iloc[i:i feature_steps] future_part df.iloc[i feature_steps gap_steps: i feature_steps gap_steps target_steps] label int((future_part[risk] risk_level).any()) rows.append({ feature: feature_part[[stress, ae_energy, microseismic_count]].values, label: label, time: df.iloc[i feature_steps gap_steps][time] }) return rows这段代码里最容易被忽略的是gap_steps。很多队伍切完窗口就直接训练训练集和验证集之间实际上只差一个步长导致交叉验证分数虚高加入间隔后特征和目标之间有一段“不知道”的时间这更接近真实预警场景。step参数是控制重叠的如果只想快速验证就设1若要评估真实泛化性能我会调到目标窗口长度以上让相邻样本尽量不共享太多观测值。特征部分我用.values直接存成数组等全部样本构造完之后再展开成二维表格。注意不要把未来窗口的计算结果合入特征部分比如把未来风险的均值放进去当特征这种错误在代码里非常隐蔽后面专门讲检查方法。2.2 评估指标先算不偏向多数类的F1准确率只当参考冲击地压危险数据有个显著特点无危险和弱危险的记录占多数真正的中等以上危险记录往往只占百分之几到十几。这时候只看准确率会骗人全部预测成无危险准确率可能到95%以上但一次漏报在矿上就是事故。如果官方没有额外说明评分规则我一般以少数类的F1为主召回率作为次要指标。from sklearn.metrics import accuracy_score, f1_score, recall_score, confusion_matrix def evaluate_model(y_true, y_prob, threshold0.5): y_pred (y_prob threshold).astype(int) cm confusion_matrix(y_true, y_pred) metrics { acc: accuracy_score(y_true, y_pred), f1: f1_score(y_true, y_pred, zero_division0), recall: recall_score(y_true, y_pred, pos_label1, zero_division0) } return cm, metrics打印出的混淆矩阵比一个总分更重要。我习惯先看右上角的漏报数如果漏报很多说明模型在有用信号上还没有学会。准确率只在报告里作为整体参考不作为调参依据。对于预警类任务漏一次和误报十次的代价完全不一样所以召回率要优先保住。2.3 先跑逻辑回归基线再用时间序列交叉验证基线模型不是用来拿奖的是用来锚定后续所有改进的值。如果先进模型连逻辑回归都打不过说明特征构造有问题而不是模型不够深。这里必须用TimeSeriesSplit而不是普通的KFold因为监测数据一旦打乱顺序训练集和验证集就会出现“知识穿越”。from sklearn.model_selection import TimeSeriesSplit from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LogisticRegression from sklearn.metrics import f1_score def run_baseline(X, y, n_splits5): tscv TimeSeriesSplit(n_splitsn_splits) fold_f1 [] for train_idx, test_idx in tscv.split(X): X_tr, X_te X.iloc[train_idx], X.iloc[test_idx] y_tr, y_te y.iloc[train_idx], y.iloc[test_idx] model make_pipeline( StandardScaler(), LogisticRegression(max_iter1000, class_weightbalanced) ) model.fit(X_tr, y_tr) fold_f1.append(f1_score(y_te, model.predict(X_te), zero_division0)) return np.mean(fold_f1), np.std(fold_f1)class_weightbalanced会按照类别比例自动放大少数类样本的惩罚权重这是最简单的一种代价敏感学习。逻辑回归在这里的作用是看线性可分能力如果F1接近于零说明原始特征在时间维度上的表达能力不足下一步就应该构造趋势和变化率而不是急着调树模型。如果在多工作面数据上做验证还要按工作面ID分组切分避免同一个工作面的相邻时段同时落在训练集和验证集里。到这一步为止基线评估基准已经确定后续所有特征和模型改动都可以跟这个数去对比。3. 从原始监测数据到特征矩阵清洗、时序特征与信息泄漏检查3.1 清洗不是全填缺失值而是区分短缺和长缺井下传感器数据拿到手第一眼往往是缺得七零八落。常见现象是某个传感器在某一时刻跳成0或者连续几个小时都没上报。这时候如果用均值填充全表等于把异常当成了正常如果直接删行时间序列又会断裂。我一般先分两类处理。短时长缺失用线性插值因为它前后值仍然有趋势关系长时长缺失直接在该字段置为NaN让后续树模型自己处理缺失值。物理量的边界准则也要跟上应力、能量这些量不可能为负出现负值按跳数处理瞬时值如果超过历史正常范围的数倍我会直接截断掉不参与特征计算。def clean_sensor(series, max_gap5, max_ratio10.0): s series.copy().astype(float) # 连续缺失超过 max_gap 的段保留为 NaN s s.interpolate(methodlinear, limitmax_gap) # 负值按异常跳数处理 s[s 0] np.nan # 超过滚动中位数若干倍的尖峰值也按缺失处理 rolling_med s.rolling(24, min_periods6).median() spike_mask (s / rolling_med) max_ratio s[spike_mask] np.nan s s.interpolate(methodlinear, limitmax_gap) return s线性插值在前后值都正常时是安全的但如果一段数据整体漂移插值结果可能把趋势做平所以max_gap不要设太大。max_ratio判断用的是滚动中位数而不是均值因为中位数对突变不敏感单个尖峰不会把基线整体抬高。清洗之后我会顺手做一个统计表把每个字段的缺失比例、跳数数量、截断数量记下来。这不是给代码跑分用的是为了在报告里回答评委的提问为什么这个字段在某些时间点没有参与建模。3.2 特征构造统计量、趋势量、事件量三组分开建监测字段不是越多越好而是每组特征都要说明它对应哪种物理意义。我习惯把特征分成三组再拼接。第一组是窗口统计量包括均值、标准差、最大值、最小值。均值反映应力水平标准差反映波动性这些都是矿压工程师一眼能看懂的。第二组是趋势量包括窗口内的线性斜率、末时刻相对窗口均值的偏离比例表达的是“这段时间在变好还是变坏”。第三组是事件量包括微震事件次数、累计能量、声发射计数直接和冲击前兆的频次变化挂钩。def build_features(df, target_cols, window12): df df.copy() for col in target_cols: roll df[col].rolling(window, min_periodsmax(3, window // 3)) df[f{col}_mean] roll.mean() df[f{col}_std] roll.std() df[f{col}_max] roll.max() df[f{col}_trend] df[col].diff(window) / (roll.mean() 1e-6) # 事件量微震频次、能量累计 if micro_count in df: df[micro_count_win] df[micro_count].rolling(window, min_periods1).sum() if micro_energy in df: df[micro_energy_cum] df[micro_energy].rolling(window, min_periods1).sum() return df注意min_periods不能设得太小否则窗口前几个点的标准差没有意义但也不能设满否则序列前一段全变成NaN树模型虽然能处理逻辑回归却会把这些行全部放弃。trend里除以均值是为了抵消不同工作面背景值的差别这组相对量在时间上更可比。对冲击地压问题我还会额外算一个类似地震学里b值的指标也就是微震事件能量分布的拟合斜率。它表达的是“小事件和大事件的比值”如果小事件减少、大事件增加往往是应力集中的信号。这个特征不具有强线性关系但放进树模型里经常能带来稳定提升。矿压领域常用的“能量释放率”也可以算成窗口累计能量除以时间长度单位是焦耳每小时表达的是冲击前能量积累的快慢。3.3 信息泄漏自检往特征里塞一个“未来标签”跑一遍特征工程做多了最怕的不是模型不行而是不知道什么时候把未来信息放进了特征。最常见的放法是计算滚动均值时窗口边界没有对齐把目标窗口的起始时刻也包了进去又或者把整份数据标准化时用了全量的均值和方差导致验证集的分布被训练集污染。我有个土办法专门用来检查滑窗有没有漏未来先在数据里加一个fake_future特征数值直接等于下一时段的标签然后照常建模。如果模型在这个假特征上投入了很高的重要性并且测试分数明显上升基本可以判断整个特征管道里存在跨时间的信息穿透。def leak_check(X, y, fake_feature_indexNone): 检查一个位置已知的未来特征能否被模型利用。 如果模型对假特征给出高重要性特征管道很可能存在泄漏。 from sklearn.ensemble import RandomForestClassifier X2 X.copy() if fake_feature_index is None: fake_feature_index X2.shape[1] X2[fake_future] 0 X2.iloc[:, fake_feature_index] y.shift(-1).fillna(0) model RandomForestClassifier(n_estimators100, random_state42) model.fit(X2, y) importance model.feature_importances_ fake_importance importance[fake_feature_index] ratio fake_importance / (importance[:-1].mean() 1e-6) return ratio这个自检要分开做两次第一次把假特征加进整条特征管道第二次把假特征放在原始数据列里从零开始构造。如果前一次重要性低、后一次重要性高问题就出在特征构造函数如果两次都高说明假特征混进了时间对齐关系要立刻检查切窗代码。另外还有一种不算泄漏但会高估的写法在切分前对全部特征做标准化用全样本的均值和方差把特征压到同一量纲。验证集信息在建模阶段就已经被看到了虽然不至于出现夸张的高分但会给阈值选择带来偏差。正确做法是先切分再只对训练部分fit再transform验证部分。信息泄漏是这类赛题里最隐蔽的扣分项切窗方式严谨了后续训练集、验证集和测试集的划分才会真正有意义。4. 不平衡数据下的模型选型与调参顺序从加权树模型到概率集成4.1 三条不平衡处理路线重采样、代价敏感、阈值移动冲击地压危险预测的正样本占比往往很低模型如果一味追求最小误差会把所有样本都判成无危险。处理不平衡有三条常见路线它们解决的问题不一样。重采样是通过增加正样本数量或减少负样本数量来平衡比例常见的是SMOTE。代价敏感是在损失函数里给少数类更高的权重比如class_weight和scale_pos_weight。阈值移动是训练时不刻意平衡但在预测阶段把判定为正类的阈值从0.5下调或上调直接控制漏报和误报的取舍。我的习惯是先做代价敏感因为改动最小、不容易产生样本分布扭曲。如果正样本实在太少再用SMOTE并且只在训练折内执行。这里最容易被忽略的是很多人直接对全体数据做了SMOTE才切分结果验证集里出现了和训练集重复的合成样本分数直接被高估。from sklearn.metrics import precision_recall_curve def pick_threshold(y_true, y_prob, min_recall0.8): prec, rec, ths precision_recall_curve(y_true, y_prob) candidates [(p, r, t) for p, r, t in zip(prec, rec, ths) if r min_recall] best max(candidates, keylambda x: x[0]) return best[2], best[0], best[1]pick_threshold的作用是把阈值选择变成一个显式环节而不是默认0.5。对冲击地压预警漏报的代价比误报高我会把min_recall设成0.8或更高再在满足召回率的前提下选精度最高的阈值。这样做的价值在于将来要给现场解释模型表现时可以直接说“我保证危险样本能报出80%在这个前提下误报率是多少”。4.2 梯度提升树的参数先定树深再定采样比最后定权重表格类数据上梯度提升树通常比神经网络稳尤其是样本量几千到几万这个量级。以LightGBM为例我一般不会一上来就套默认参数而是按下面的顺序来调。import lightgbm as lgb model lgb.LGBMClassifier( objectivebinary, learning_rate0.05, n_estimators600, num_leaves15, max_depth4, min_child_samples30, subsample0.8, colsample_bytree0.8, scale_pos_weight(y_train 0).sum() / max((y_train 1).sum(), 1), random_state42 )num_leaves和max_depth决定树模型的复杂度在这类噪声多、样本少的数据上树太深会记住局部波动。min_child_samples设到30以上是防止某些叶节点只含少数几条记录。scale_pos_weight直接按负样本与正样本的比值填不需要太细先让它把少数类召回率提上来。迭代次数交给早停来控制而不是拍一个600。把训练集拆出一个时间切片做验证只要连续50轮验证分数不再上升就停止。这样即使后面改了特征也不需要重新手工调n_estimators。调参顺序也有讲究。先固定learning_rate0.05去调num_leaves和max_depth这两项对过拟合影响最大然后调min_child_samples和采样比最后才动scale_pos_weight。一次只动一个参数记录下每个组合在时间序列交叉验证上的F1避免多个参数同时调整时分不清收益来自哪里。4.3 集成用概率加权而不是硬投票权重在验证折里搜多个模型做集成是赛题里常见加分项但硬投票对概率阈值敏感不如直接在预测概率上做加权平均。权重可以用验证集网格搜索来定不要用手拍。做法是先把两三个模型的预测概率保存下来然后搜索权重组合。def search_ensemble_weights(prob_dict, y_val): names list(prob_dict.keys()) best None best_score 0.0 for a in np.linspace(0, 1, 11): for b in np.linspace(0, 1 - a, 11): c 1 - a - b prob a * prob_dict[names[0]] b * prob_dict[names[1]] c * prob_dict[names[2]] score f1_score(y_val, (prob 0.5).astype(int), zero_division0) if score best_score: best_score score best (a, b, c) return best, best_score网格搜索直接看F1可以用来找最优权重。它有三个注意点。第一权重必须在验证集上搜索不能在训练集上搜索否则会过拟合训练集。第二参与集成的模型预测概率要先经过统一的概率校准否则一个模型的概率整体偏高会把权重带偏。第三搜索得到的权重不要保留太多小数尽量取0.1或0.2的粒度避免部署时出现没有意义的精度。这种概率加权的另一个好处是还能继续接阈值移动因为输出的还是概率。现场如果想降低误报只需要把最后的判定阈值调高不需要重新训练模型。集成带来的提升通常在5个百分点以内如果某个单模型F1只有0.3另一个有0.6集成后最多也就是0.6附近救不了模型本身的短板。因此做集成前先把单模型调到一个合理水平。5. 冲击地压预测的五个常见坑与排查清单先过时间关再过平衡关不管特征工程做得多细模型选得多花最后被问倒的往往是一些看起来很小的地方。下面这几个坑是我在类似监测数据项目里反复见过的每一条都会直接影响最终结果的可靠性。5.1 随机切分交叉验证导致分数虚高现象同一个模型用随机KFold跑出来的F1是0.9换成时间序列切分之后掉到0.5。原因监测数据是时间序列相邻记录之间高度相似。随机切分会把同一时间段的样本同时分进训练集和验证集模型实际上见过验证集附近的数据。解决统一用TimeSeriesSplit或按时间粒度分块验证。分块时还要保证同一工作面的数据不会跨块否则仍然存在邻近泄漏。如果赛题数据本身包含多个工作面请把工作面ID也作为分层的依据别让模型通过位置信息偷看验证结果。5.2 滑窗重叠过大样本独立性被破坏现象在时间序列切分下F1依然远高于预期但把预测结果画回时间轴后发现模型只是把上一个预测结果复制了一遍。原因滑窗时step1样本间只差一条记录重叠率很高。解决把step调到目标窗口长度以上或者训练集和验证集之间多留一段缓冲区域。前面make_samples里的gap_steps就是为这个准备的验证时要检查训练折和验证折之间的时间间隔是否足够。5.3 SMOTE在全体数据上重采样验证集被污染现象SMOTE之后F1提升明显但用新一批数据测试时性能掉得厉害。原因先对全体数据做了过采样再切分训练集和验证集验证集里含有由训练集生成的合成样本的复制或近似样本评估结果不再公正。解决只在训练折内做SMOTE最好的办法是用imblearn.pipeline把重采样和训练模型包在一起保证每次交叉验证都只对训练部分过采样。5.4 只盯准确率把危险样本全部漏掉现象准确率97%混淆矩阵显示正样本召回率只有0.1。原因负样本占比过高模型把全部样本判为负类准确率依然很高。解决调参过程用F1、召回率、精确率做主要指标准确率放到报告最后展示。对冲击地压这种安全类问题漏报的后果大于误报要优先保证正样本召回率。5.5 把预测概率0.5直接当报警阈值形成大量误报现象模型训练时F1不错部署后每天报警几十次现场很快失去信任。原因训练目标只关心分类正确率没有把误报成本纳入目标函数0.5这个阈值本身就是人为约定。解决用precision_recall_curve在验证集上选择报警阈值。面向现场使用时把阈值设到让精确率不低于某个可接受水平同时把对应的漏报率写清楚。这个阈值必须在报告中写明评委或矿方最关心的是“漏报率多少、误报率多少”。这五关过完模型就具备了一个能摆上台面的基础。剩下的工作不在于继续堆模型而是验证模型换了时间段之后还稳不稳。6. 让模型多扛一层验证滚动回测、概率校准和物理反查6.1 滚动回测模拟边采边测的真实推进过程静态交叉验证再严谨也只是一次性切分。现场更关心的是如果从某一天开始每周用最新数据重训一次模型后面两个月预警准不准。滚动回测正是解决这个问题。def rolling_backtest(X, y, model, train_size400, test_size50): start 0 scores [] while start train_size test_size len(X): idx_tr np.arange(start, start train_size) idx_te np.arange(start train_size, start train_size test_size) model.fit(X.iloc[idx_tr], y.iloc[idx_tr]) prob model.predict_proba(X.iloc[idx_te])[:, 1] scores.append(f1_score(y.iloc[idx_te], prob 0.5, zero_division0)) start test_size return scores滚动回测的特点是每个测试块都严格在训练块之后和真实部署一致。它能看出的问题主要是特征是否跨时间退化。如果第一段时间评分高、后面骤降很可能是因为井下生产条件变化模型需要按周期重训。6.2 概率校准输出分档而不是单一预测值如果评估方要求输出危险等级而模型给你的是概率直接设0.5阈值转换会丢失很多信息。我习惯在验证集上查看概率分布再按分位点做区间而不是按0.25、0.5这种固定比例切。一个具体做法是把验证集里预测概率按从小到大排序取0.2和0.6两个分位点切出三个区间对应弱、中、强风险。这样分档的好处是每个区间的样本占比相对合理。但要注意分位点会因为数据分布漂移而失效上线之后还需要定期重新校准避免长期使用一个旧分位点。6.3 物理反查把预测结果拉回时间轴看能否对上冲击事件最后一步也是最容易被人忽视的一步把模型输出的概率序列和实际冲击地压事件画在同一个时间轴上。如果事故前没有看到概率抬升证明模型的预警价值不足如果事故前概率明显走高再把误报的时段列出来去查这几个时段的应力状态是否有特殊性。例如有段时间应力并不高但模型给出高概率这时候多半是特征里混入了其他工作面的信号或者传感器存在长期偏差。这个反查不需要复杂的代码一张带事件标记的时间序列图就够。真正麻烦的是把“预测结果能否解释”作为模型选型的一个硬性条件。我在几个类似项目里吃过亏模型调得很高但没有任何人能解释它为什么发出高报警最后现场根本不敢用。所以后来我落地的习惯是先画时间序列确认预测趋势再做概率分档最后才写报告。获奖作品和普通作品的分水岭也在这个地方复杂公式谁都会套但能做到时间切分干净、漏报误报有过交代、预测结果和物理背景对得上才会让评委相信这套方法不是只在训练集上好看。希望这个流程帮到你。本文还有配套的精品资源点击获取