简介2024年五一赛C题获奖作品围绕煤矿深部开采冲击地压危险预测系统整理电磁辐射与声发射数据的分析流程适合数学建模备赛者与矿山安全数据人员研读。作品综合运用傅里叶变换、随机森林和朴素贝叶斯模型先通过最大值、最小值、均值、标准偏差、峰度与偏度等统计量刻画干扰信号特征再定位干扰区间并量化前兆概率两类信号均给出了明确的规律对比和计算数值。压缩包内仅1个PDF文件约1.16MB章节覆盖问题重述、数据预处理、模型训练、结果分析及结论展望便于完整还原获奖思路。目前已有1468人学习可直接用作赛题复盘参考帮助读者掌握从信号特征提取到危险概率预测的系统建模方法。1. 冲击地压预测研究的难点不在模型而在“事件太少、噪声太多”2024年五一赛C题煤矿深部开采冲击地压危险预测研究表面上是训练一个分类模型实际上考的是从一份带缺失、带噪声、带人工干扰的真实监测数据里把事件稀疏的危险时段捞出来。冲击地压是深部开采中最常见也最难防的动力灾害应力、位移、速度、微震、爆破记录这些监测数据都不缺但真正标记为危险的有效样本往往只占总样本的5%到10%。很多队伍把这道题当成普通二分类来做结果模型一路预测“安全”还能拿到85%以上的准确率这就是典型的被数据分布骗了。要在这道题上做出有说服力的结果绕不开四个环节统一时间轴、构造窗口特征、时序评估、解释性验证。这篇笔记按这个顺序拆解中间会单独拉出几个高发翻车点帮备赛队伍把时间花在刀刃上。2. 从监测日志到建模特征合并、清洗、造特征三步缺一不可2.1 两张原始表的时间粒度不同先统一主时间轴比赛数据通常分成监测表和事件日志两类。监测表按固定间隔采样常见的是15分钟或3小时字段一般包含日期时间、测点编号、应力值、位移值、速度值。事件日志则按事件触发记录微震事件只有发生时间和震级、能量爆破记录只有每次爆破的时间点。这两张表直接合并行不通时间戳对不齐必须先确定一条主时间轴。我一般把监测表的时间戳作为主时间轴原因很简单模型最后输出的“危险概率”要落在特定的时间点上而监测表正好提供了等间隔的采样点。微震日志和爆破记录都按这条时间轴做聚合变成“近3小时发生了多少次微震”“累计能量是多少”这类统计特征。下图这个思路看着简单实际操作里很多队伍在第一行read_excel之后就急着跑模型到后面才发现日期格式没统一merge出来的结果是笛卡尔积。import pandas as pd import numpy as np # 附件1 监测表日期时间列直接解析成 datetime df pd.read_excel(附件1_监测.xlsx, parse_dates[日期时间]) df df.sort_values([测点编号, 日期时间]).reset_index(dropTrue) # 通过时间差分布确认主采样周期 delta_h df.groupby(测点编号)[日期时间].diff().dt.total_seconds() / 3600 print(delta_h.dropna().describe()) # 粗看各列缺失率 print(df.isna().mean().sort_values(ascendingFalse).round(4))这段代码先排序再按测点编号分组求时间差目的是确认监测表是不是真的“等间隔”采样。很多矿井监测系统会在中途调整采样频率如果describe结果里出现两个明显不同的时间差簇后面做滚动窗口时就要格外注意对齐问题。缺失率打印出来之后不要直接dropna删除时序数据行在时间轴上有位置删掉会让滞后特征整体错位。参数方面parse_dates指定的时间列要确保Excel里没有多余空格或合并单元格groupby里的测点编号要在排序时一并考虑否则同一测点的记录会被其他测点插空。如果表里不止一个测点建议先只保留目标工作面做单测点建模最后再用另一个测点做迁移验证这样比把所有测点混在一起更干净。2.2 异常值与缺失值处理一个不会引入未来信息的清洗流程冲击地压监测数据的异常来源主要有三类传感器瞬时掉线产生的毛刺、开采震动造成的数值跳变、人工维护时的错误读数。处理这类点不能直接删行正确做法是先置为NaN再通过插值补回。我见过有人把应力值超过某个阈值的行全部删掉结果时间序列上出现一个大洞后面构造diff特征时直接出现一个夸张的负值。def smooth_spike(df, col, window7, threshold3.0): # 以局部均值为中心局部标准差为尺度识别尖峰 roll_mean df[col].rolling(window, centerTrue).mean() roll_std df[col].rolling(window, centerTrue).std() z (df[col] - roll_mean) / roll_std df.loc[z.abs() threshold, col] np.nan return df for col in [应力值, 位移值, 速度值]: df smooth_spike(df, col) # 短缺口线性插值长缺口用同小时均值回填 df[应力值] df[应力值].interpolate(methodlinear, limit_directionforward) df[应力值] df[应力值].fillna( df.groupby(df[日期时间].dt.hour)[应力值].transform(mean) )smooth_spike的逻辑是以当前点前后各3个样本window7包含中心点的局部均值为参照偏离超过3倍局部标准差就判定为尖峰置为NaN。threshold3.0对应3西格玛对监测数据来说偏保守如果毛刺很密集可以考虑下调到2.5。这个参数不要盲调先画一段原始曲线看毛刺幅度再决定阈值。填补阶段先用线性插值补短缺口一般连续缺口不超过3个采样点都适合。超过插值范围的长缺口再用“同小时均值”回填比如所有发生在凌晨2点的应力值均值这样做保留了昼夜周期也不像全局均值那样把白天夜晚压成一个水平。需要注意如果原始数据不是等间隔采样线性插值就不可靠这也是前面先确认主采样周期的原因。2.3 微震事件与爆破记录聚合统计和干扰剔除微震日志和爆破记录都是事件型数据不能直接拼到监测表里当特征常见做法是把它们按主时间轴切块每3小时一个桶统计桶内事件数、最大震级、总能量、平均能量再合并到监测表上。ev pd.read_excel(附件2_微震.xlsx, parse_dates[发生时间]) ev[时间桶] ev[发生时间].dt.floor(3h) ev_agg ev.groupby(时间桶).agg( 微震次数(震级, count), 最大震级(震级, max), 总能量(能量, sum), 平均能量(能量, mean) ).reset_index() df[时间桶] df[日期时间].dt.floor(3h) df df.merge(ev_agg, on时间桶, howleft) df[微震次数] df[微震次数].fillna(0).astype(int)dt.floor(3h)把时间戳对齐到最近的前一个3小时整数点比如08:00、11:00、14:00这样微震事件和监测记录就共享同一个桶key。合并后没有微震记录的桶会得到NaN“微震次数”用0填充是合理的因为“没发生事件”本身就是信息。但“平均能量”空桶不能简单填0没有事件和事件能量真为0是两回事我一般让平均能量保持NaN让模型把缺失当作一种状态。爆破干扰是这类赛题里最隐蔽的坑爆破产生的高能量震动会被传感器记录成微震事件如果不处理模型会天然认为“爆破后半小时能量大等于有危险”这对真实预测毫无意义。做法是给爆破时间前后各15分钟窗口内的微震事件打上干扰标记再决定剔除还是单列特征。def mark_blast_window(ev, blast, mins15): ev[爆破干扰] 0 for bt in blast[爆破时间]: mask (ev[发生时间] bt - pd.Timedelta(minutesmins)) \ (ev[发生时间] bt pd.Timedelta(minutesmins)) ev.loc[mask, 爆破干扰] 1 return ev ev mark_blast_window(ev, blast) ev_clean ev[ev[爆破干扰] 0]mins参数要根据爆破密集程度调整如果两次爆破间隔本身就小于半小时窗口缩到10分钟更合适。我在处理类似数据时发现爆破记录和微震时间表高度重合的情况下不剔除爆破干扰的“最大震级”特征几乎等价于“是否刚爆破”训练F1虚高不少。剔除后重跑模型分数通常下降一点但这才接近真实泛化能力。千万别为了保训练分数把干扰样本留着。2.4 构造三类核心特征水平、变化率与波动性冲击地压的危险往往不是某一瞬间数值高低决定的而是趋势异常。所以原始监测值进入模型前要延展出三类特征一是水平特征看当前值和前几步的值二是变化率特征看相邻时刻的差值三是波动性特征看过去一小段时间的均值和标准差。for col in [应力值, 位移值, 速度值]: df[f{col}_lag1] df[col].shift(1) df[f{col}_lag2] df[col].shift(2) df[f{col}_diff1] df[col].diff(1) df[f{col}_ma3] df[col].rolling(3).mean() df[f{col}_std3] df[col].rolling(3).std()lag1和lag2做的是滞后特征让当前时刻能“看到”前两步的状态对时间序列模型很重要。diff1是相邻两点的一阶差分代表变化速率冲击地压发生前应力变化率往往会出现持续抬升。ma3和std3分别代表过去三个采样周期的均值和波动程度窗口大小和数据粒度强相关3小时粒度下window3覆盖9小时如果改成15分钟粒度window可能要放到12以上才合理。这里有个容易忽略的点不能把时间戳本身直接当特征。日期时间列如果转成数值丢进模型很容易让模型按时间分段记忆而不是学到物理特征最后在盲测中严重过拟合。时间信息要通过滞后、差分、滚动统计体现而不是用原始时间戳。2.5 特征筛选与业务保留相关性过滤后还要看物理含义特征造完之后容易膨胀尤其加入多组滚动统计后特征数可能到30以上。先用方差过滤掉近常数特征再按相关性去冗余。# 去掉方差过低的列 low_var_cols df.var()[df.var() 1e-6].index df df.drop(columnslow_var_cols) # 相关性大于0.9的两列只保留一个 corr df.corr().abs() upper corr.where(np.triu(np.ones(corr.shape), k1).astype(bool)) drop_cols [c for c in upper.columns if any(upper[c] 0.9)] df df.drop(columnsdrop_cols)相关性去重的逻辑很简单但业务含义要小心应力值和位移值相关系数很高不代表两者在冲击地压预警里可以随便删一个。获奖作品里常见做法是先在同一监测类型内部去重比如应力值的lag1和lag2高度相关时只留lag1但在应力值和位移值之间即使相关性高也会优先保留“变化率”类特征而不是水平值。因为冲击地压前的变化率特征比绝对数值更有区分度这是物理背景决定的不是相关系数能衡量的。3. 冲击地压危险预测模型逻辑回归到集成学习的三层对比3.1 标签构造预警窗口长度直接决定模型上限大部分情况下题目给的不是现成的二分类标签而是一份冲击地压发生时刻的清单。模型要预测的不是“某一秒突然发生”而是在事件发生前的一段时间内发出预警。那y标签怎么打常见做法是把事件前某个时间窗口内的所有样本标记为1。# 已知的冲击地压发生时刻从题目或附件中提取 danger_times pd.to_datetime([2024-03-01 06:00, 2024-03-08 22:00]) df[label] 0 for t in danger_times: mask (df[日期时间] t - pd.Timedelta(hours12)) (df[日期时间] t) df.loc[mask, label] 1这段代码把每个冲击地压事件前12小时标记为危险期12小时就是预警窗口。窗口太短比如3小时模型能学到更精确的前兆但漏报率会很高因为很多前兆从12小时前就开始累积窗口太长比如36小时召回率上来了但样本里危险比例太高误报也跟着涨。预警窗口本质上是业务参数不是模型参数。如果题目直接给定了危险等级列那就不需要构造窗口标签而是做成多分类或者回归问题。但多数赛题给的是事件时刻或等级不完整的记录窗口构造这一步基本绕不开。标签构造完成后要检查正负样本比例如果正样本低于10%直接进模型是不行的。3.2 不平衡样本class_weight 优先于过采样冲击地压数据的正样本比例天然很低我处理过的类似监测数据里危险标签大多在3%到10%之间。这种不平衡下模型倾向于把所有样本都预测成安全。第一反应是SMOTE但在时间序列上使用SMOTE非常危险——少数类样本的插值会跨越时间索引生成一个“过去和未来混合”的合成样本训练和验证的信息边界被打破。from sklearn.ensemble import RandomForestClassifier model RandomForestClassifier( n_estimators200, max_depth8, min_samples_leaf10, class_weightbalanced_subsample, random_state42, n_jobs-1 )class_weightbalanced_subsample是随机森林里处理不平衡的经典配置它让每个决策树的随机子样本都按类别比例自动加权比手动设置权重更稳健。相比SMOTE它不会凭空生成样本只是调整损失贡献不会破坏时序结构。逻辑回归里则用class_weightbalanced效果类似。如果调查比赛评分标准里“命中率”权重很高可以把正样本的权重再调大一点。但盲目把权重调特别大也不行模型会把大量正常时段误判为危险precision掉得非常快。一般做法是先跑一版balanced看混淆矩阵里漏报和误报的比例再决定往哪个方向调整。3.3 模型训练与调参以 F1 为目标的网格搜索这类题目我一般会对比三层模型逻辑回归做基线随机森林做主力LightGBM做提升尝试。逻辑回归作为基线的意义在于判断特征工程是否有效如果线性模型的F1和随机森林差距不大说明特征之间的非线性交互还没被充分挖掘需要回到特征工程而不是继续调参。from sklearn.model_selection import TimeSeriesSplit, GridSearchCV from sklearn.metrics import f1_score, make_scorer param_grid { n_estimators: [100, 200], max_depth: [6, 8], min_samples_leaf: [5, 10] } grid GridSearchCV( RandomForestClassifier(class_weightbalanced_subsample, random_state42), param_grid, scoringmake_scorer(f1_score), cvTimeSeriesSplit(n_splits5) ) grid.fit(X_train, y_train) print(grid.best_params_)这里评分函数用f1而不是accuracycv用TimeSeriesSplit而不是KFold这是两个关键决策。TimeSeriesSplit保证训练集始终在测试集之前不会让“未来”数据混进训练。在时间序列问题上随机网格搜索加时序交叉验证的结果才是可信的。LightGBM在表格数据上往往比随机森林更快更好但它对缺失值和噪声更敏感。冲击地压数据里不少特征是经过插值回填的LightGBM会把插值痕迹当真实信号学进去所以我的经验是先用随机森林定特征和窗口最后再上LightGBM微调。如果LightGBM表现明显好于随机森林要警惕是不是学到了插值痕迹而不是物理规律。3.4 评估指标为什么不能只看准确率准确率这个指标在不平衡数据上会严重误导人。假设98%的样本是安全状态模型完全放弃预测危险准确率也是98%看起来“优秀”得离谱。冲击地压预测关心的是漏报有多大代价——一次漏报可能意味着矿压事故而一次误报只是增加检查成本。所以召回率比准确率重要precision则决定预警系统会不会被人当成狼来了。常用指标组合是F1分数和ROC-AUC曲线其中F1综合了precision和recallROC-AUC反映模型在不同阈值下的排序能力。如果比赛评分里有“预警命中率”这一类指标通常等价于在保证一定召回率的前提下最大化precision可以输出概率后用阈值扫描。from sklearn.metrics import precision_recall_curve, f1_score y_proba model.predict_proba(X_valid)[:, 1] prec, rec, thresh precision_recall_curve(y_valid, y_proba) f1_scores 2 * prec * rec / (prec rec 1e-9) best_idx f1_scores.argmax() best_thr thresh[best_idx] print(最优阈值:, best_thr, F1:, f1_scores[best_idx])阈值扫描这一步在竞赛中非常重要model.predict默认用0.5做阈值但对不平衡数据0.3甚至0.2可能是更优的预警线。获奖作品通常会在论文里给出阈值选择的曲线图说明为什么最终选用某个阈值这也是模型落地性和可解释性的体现。4. 时序交叉验证与预警窗口敏感性预测模型的落地评估4.1 时间序列交叉验证与随机划分的差距很多队伍用train_test_split随机打乱数据然后在测试集上拿到不错的F1最后提交盲测却一塌糊涂。原因在于滞后特征和滚动统计让相邻样本之间存在强相关性随机划分的测试集里包含大量和训练样本紧邻的观测模型几乎等于“记住了”上下文。from sklearn.model_selection import TimeSeriesSplit tscv TimeSeriesSplit(n_splits5) for fold, (tr_idx, te_idx) in enumerate(tscv.split(X)): print(ffold {fold}: train {len(tr_idx)}, test {len(te_idx)})TimeSeriesSplit按时间顺序切分每次测试集都在训练集之后模拟的是“用过去预测未来”的真实场景。在冲击地压预测里时序交叉验证得到的分数通常会比随机划分低5到10个百分点这个差距是合理的。如果成绩单上时序验证分数和随机划分几乎一样要怀疑是不是特征构造里引入了异常信息。随机划分还有一个衍生问题同一个冲击地压事件前12小时的样本可能一部分落在训练集、一部分落在测试集模型在测试时看到的是“同一个危险过程的中间状态”这属于典型的数据泄漏。用时间序列交叉验证可以彻底避免这种“穿越”。4.2 不同预警窗口下的召回率与误报实验预警窗口这个设计变量值得单独做一次敏感性实验。常见做法是固定模型和特征只改变3.1节里的窗口长度分别记录召回率、precision和误报率然后对比。预警窗口召回率Precision误报次数/月备注6小时0.780.713~4漏报少误报可接受12小时0.820.655~6召回与误报较均衡24小时0.880.538~10高召回但误报偏多上面这个表是一个典型结果模式不同数据会有差异。核心规律是窗口拉长召回率上升precision下降。12小时在多数场景下是召回率和误报的平衡点但这个结论必须从实验里得出来不能拍脑袋。我在这类赛题里见过两篇获奖论文共同点是都会附一张窗口敏感性曲线图横轴是预警窗口纵轴是F1和误报率。这个实验的目的是证明自己的模型确实有实用价值而不是正好在某个参数上撞了好结果。4.3 停采期与正常生产期模型稳定性测试冲击地压监测数据里常有一个特殊时段矿井停采检修。停采期间监测系统可能关闭数据出现大段缺失就算插值回填物理意义也和正常开采完全不同。模型如果在停采恢复初期连续误报说明它把插值痕迹当成了危险信号。# 假设已知停采时间区间 stop_periods [(2024-02-01, 2024-02-07)] df[阶段] 正常 for start, end in stop_periods: mask (df[日期时间] start) (df[日期时间] end) df.loc[mask, 阶段] 停采把样本按阶段分组分别计算F1和召回率如果停采期的召回率明显低于正常期说明模型的可靠性和数据插值策略绑定太紧。这时有两个调整方向一是停采期单独做一版模型二是把“距离最近一次停采的天数”作为特征加入模型。第二个方向在论文里更好写因为它揭示了模型对矿井生产状态的适应性。5. 冲击地压危险预测里的四个高发坑现象、原因、解决5.1 坑1全样本标准化导致跨时间信息泄漏现象是模型在交叉验证里F1很高放到最后一段盲测数据上暴跌。原因是很多队伍在建模前用StandardScaler对整个X做了fit和transform标准化器里包含了未来的均值和方差。时间序列问题里测试集的均值和方差本来就不该被模型知道。解决方法是把标准化放进交叉验证循环内部每一折只用训练集的数据拟合scaler再用训练好的scaler转换测试集。# 错误示范 scaler StandardScaler().fit(X) X_scaled scaler.transform(X) # 正确示范 from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_train scaler.fit_transform(X_train) X_test scaler.transform(X_test)这条规则对缺失值填充同样适用。如果用全样本均值填充缺失严格来说也引入了未来信息。训练集和测试集的缺失值要分别处理哪怕填充策略一样拟合参数也必须来自训练集。5.2 坑2用准确率评价模型危险样本被完美无视现象是模型训练完成后准确率85%提交结果几乎全部预测为“安全”但比赛评分的命中率部分得了0分。原因是整个建模流程里只盯着准确率没有看一眼混淆矩阵。解决方法是训练完第一个模型后立刻打印混淆矩阵和分类报告。如果发现recall低于0.5说明模型没有真正学会识别危险。处理手段包括给正样本加权、降低预警阈值、在特征里补充微震能量趋势等。准确率只能作为参考真正该盯的是召回率和F1。5.3 坑3随机切分训练测试集结果虚高却不自知现象是随机划分下F1有0.85换成TimeSeriesSplit之后掉到0.72队伍开始怀疑模型不稳定。其实模型没有变是评估方式暴露了真实水平。冲击地压数据高度时间相关随机划分天然制造了数据泄漏。解决方法是固定使用TimeSeriesSplit并且在论文里明确写出“所有评估均基于时序交叉验证”。这是一个很小的代码改动但对评分影响巨大。评委看过太多随机划分虚高、实际泛化拉胯的作品你能主动用时序验证本身就是加分项。5.4 坑4爆破微震事件没剔除模型学到的是爆破现象是模型的特征重要性排行里“最大震级”排第一而且预测精确率出奇高。原因是爆破产生的高能量震动被记录成微震事件而这些时间点和危险标签大概率不相关模型学到的是“爆破后一段时间有高能量危险”等于在学爆破规律而不是冲击地压规律。解决方法是回到2.3节把爆破前后窗口内的微震事件标记出来并剔除或者单独加一个“是否爆破后30分钟内”的特征。后者在论文里更好解释因为它展示了队伍对生产干扰因素的显式建模能力。剔除后模型F1通常会下降但提交后的盲测分数更稳定这一点值得在答辩时主动提出来。6. 从入围到获奖的三个进阶技巧窗口实验、特征签名与呈现6.1 把预警窗口当成超参数给出敏感性证据普通作品只会说“我们用了12小时预警窗口”获奖作品会给出窗口从3小时到36小时的完整实验结果说明为什么选12小时。做法很简单windows [3, 6, 12, 24, 36] for w in windows: y make_label(danger_times, w) # 复用标签构造函数 f1, rec train_and_eval(X, y) # 时序交叉验证 print(w, f1, rec)这个实验的结论通常是一条先升后降的曲线。F1峰值对应的窗口就是业务上的最佳预警时间。把这个过程写进论文比直接给一个结论有说服力得多。6.2 用 SHAP 提炼危险特征签名给评委一份解释模型拿到F1之后下一步是解释“什么特征在驱动预警”。SHAP是树模型解释的常见工具可以输出每个特征对预测结果的贡献方向。冲击地压数据里我常看到的结果是累计微震能量、应力变化率、顶板位移速度这三个特征在危险样本中显著偏高。import shap explainer shap.TreeExplainer(grid.best_estimator_) shap_values explainer.shap_values(X_test) shap.summary_plot(shap_values, X_test, feature_namesX_test.columns)把SHAP结果提炼成两三句话冲击地压发生前24到48小时微震累计能量呈现阶梯式上升同时应力变化率持续为正这两个信号叠加构成主要预警依据。这部分内容在评委眼里比一堆调参实验更有价值因为它说明模型不是黑匣子。6.3 呈现与答辩三张图讲完整条逻辑链第一张图是时间轴上的危险概率曲线横轴是日期纵轴是模型输出概率把真实事件时刻标注成竖线直观展示模型提前多久发出预警。第二张是特征重要性排序图重点标出应力变化率、累计能量、顶板速度这三个核心特征。第三张是预警窗口敏感性曲线给出F1随窗口变化的趋势。这三张图配合时间序列交叉验证的结果基本上把“是什么数据、怎么做特征、为什么可信”讲清楚了。我做这类题目习惯把答辩时间分成三块两分钟讲数据清洗和特征构造三分钟讲模型对比和窗口实验两分钟讲SHAP特征签名和落地建议。评委问到细节时再从答辩材料里翻出对应的表格而不是临场翻代码。希望这些经验能帮你们少走几个月的弯路。冲击地压预测这道题模型算法一年比一年卷但真正拉开差距的始终是数据理解、时序意识和解释深度——这三样做扎实了获奖是水到渠成的事。本文还有配套的精品资源点击获取