☰
煤矿冲击地压预测实战:傅里叶变换+随机森林+朴素贝叶斯全流程解析
2026/9/26 2:57:32 网站建设 项目流程

简介:本资源为2024年五一数学建模竞赛C题「煤矿深部开采冲击地压危险预测研究」的获奖作品,面向备战数学建模竞赛的高校学生及从事矿山安全数据分析的研究人员。作品围绕电磁辐射与声发射信号,综合运用傅里叶变换、随机森林与朴素贝叶斯模型,完成干扰信号识别、前兆特征提取及危险发生概率预测,可作为同类赛题的完整解题范本。资源包共1个PDF文件,约1.16MB,内容涵盖问题重述、模型假设、数据分析方法、实验步骤与结果讨论,并给出各时间段前兆特征概率等关键结论,便于读者对照学习建模思路与论文写作结构。目前已有1347人学习下载,适合需要系统掌握信号处理与机器学习建模流程、查漏补缺的参赛者参考。

1. 从一份获奖作品说起:这套煤矿冲击地压预测方案到底能跑出什么

煤矿深部开采的冲击地压预测,听起来离日常开发很远,但它本质上是一个典型的时序信号分类问题:传感器每 30 秒采一个点,电磁辐射和声波强度两条曲线里混着正常工作、前兆特征、干扰信号、传感器断线、工作面休息五类状态,你要做的就是从噪声里把危险的那几段揪出来。这套 2024 年五一赛 C 题的获奖作品,核心链路是傅里叶变换做频域特征、随机森林做区间识别、朴素贝叶斯做概率预判,三个模型串起来覆盖了"找干扰—找前兆—算概率"的完整闭环。它适合正在做时序异常检测、工业信号分类的从业者参考,也适合想找一个完整建模流程练手的人。我把这份作品的代码和数据流拆开跑了一遍,下面把能复现的部分、参数怎么设、哪里容易翻车讲清楚。

2. 傅里叶变换加随机森林:干扰信号识别的完整链路

2.1 为什么先做傅里叶变换而不是直接上原始信号

原始电磁辐射和声波强度数据是典型的时间序列,直接丢进分类器会有一个很现实的问题:单点数值的区分度不够。正常工作状态下电磁辐射均值在 49.81 左右,干扰信号均值 77.96,看着有差距,但标准偏差分别是 18.28 和 90.72,方差比均值差异大得多,说明干扰信号的波动本身就淹没了均值信息。傅里叶变换的作用是把时域上的剧烈波动转成频域上的能量分布,变换后的序列更平缓,不同类别的频域特征差异反而更稳定。

作品里的做法是对每个采样点取前后各 10 个点组成 21 点滑动窗口,对窗口做 FFT,取直流分量(a[0].real)和频谱方差(a.var())作为两个核心特征,再加上窗口内的均值和方差,凑成 4 维特征向量。这个窗口大小不是随便定的——21 个点对应约 10 分钟的采样跨度,既能覆盖一次干扰事件的持续时间,又不至于把相邻事件混在一起。

2.2 特征工程的代码实现与参数说明

下面这段是作品里电磁辐射 C 类数据(干扰信号)的特征提取代码,我整理成了可直接跑的版本:

import numpy as np import pandas as pd import matplotlib.pyplot as plt # 读取附件1,注意 sheet 名和列名要和实际文件对齐 data = pd.read_excel('附件1.xlsx', engine='openpyxl', skiprows=0) index = data[data['类别'] == 'C'].index zbiao = np.zeros((data.shape[0], 4)) data0 = data['电磁辐射'].to_numpy() for i in range(data0.shape[0]): if i % 100000 == 0: print(i) # 进度提示,数据量大时很有用 if i - 10 < 0: # 边界处理:开头不足10个点时取前 i+11 个 a = np.fft.fft(data0[:i+11]) zbiao[i, 0] = data0[:i+11].mean() zbiao[i, 1] = data0[:i+11].var() zbiao[i, 2] = a[0].real # 直流分量 zbiao[i, 3] = a.var() # 频谱方差 elif i + 10 >= data0.shape[0]: # 边界处理:末尾不足10个点时取后段 a = np.fft.fft(data0[i-10:]) zbiao[i, 0] = data0[i-10:].mean() zbiao[i, 1] = data0[i-10:].var() zbiao[i, 2] = a[10].real zbiao[i, 3] = a.var() else: # 正常情况:21点滑动窗口 a = np.fft.fft(data0[i-10:i+11]) zbiao[i, 0] = data0[i-10:i+11].mean() zbiao[i, 1] = data0[i-10:i+11].var() zbiao[i, 2] = a[10].real zbiao[i, 3] = a.var() # 生成标签:C类为1,其余为0 y = np.zeros(data0.shape[0]) y[index] = 1 np.savetxt('电磁辐射数据.txt', np.hstack((zbiao, y.reshape(-1, 1))))

逻辑上分三步走:先按类别筛出 C 类索引做标签,再对每个点做 21 点窗口的 FFT 和统计量计算,最后把特征和标签拼在一起存成 txt 供随机森林读取。参数上最需要注意的是窗口半径 10 这个值——它决定了特征的时间分辨率。如果你把窗口改小到 5,特征会变得对瞬时波动过于敏感,干扰信号和正常信号的边界反而模糊;改大到 20,一次短时干扰会被稀释在长窗口里,召回率会掉。作品选 10 是经过对比的,建议不要随意改。

声发射数据的处理逻辑完全一样,只是把列名从电磁辐射换成声波强度,sheet 名换成AE。两个文件分别跑完,得到两份特征表。

2.3 随机森林的输入构造与区间输出

特征表有了,接下来是随机森林。作品没有贴出随机森林的完整训练代码,但从描述看,输入是上面生成的 4 维特征加标签,输出是每个采样点的类别预测。这里有一个关键细节:随机森林输出的是逐点分类结果,而题目要的是"时间区间",所以中间需要一个后处理步骤——把连续被预测为 C 类的点合并成区间,再按时间排序取最早的 5 个。

我一般会这样写后处理:

from sklearn.ensemble import RandomForestClassifier from itertools import groupby # 读取特征和标签 feat = np.loadtxt('电磁辐射数据.txt') X, y = feat[:, :4], feat[:, 4] # 训练随机森林,n_estimators 建议 100-200 rf = RandomForestClassifier(n_estimators=150, max_depth=12, min_samples_leaf=5, random_state=42) rf.fit(X, y) # 预测并合并连续区间 pred = rf.predict(X) intervals = [] start = None for i, p in enumerate(pred): if p == 1 and start is None: start = i elif p == 0 and start is not None: intervals.append((start, i-1)) start = None if start is not None: intervals.append((start, len(pred)-1)) # 按区间长度过滤,去掉过短的噪声区间 intervals = [iv for iv in intervals if iv[1] - iv[0] >= 3] print(f"共识别出 {len(intervals)} 个干扰区间")

n_estimators设 150 是精度和速度的折中,max_depth限制到 12 防止过拟合,min_samples_leaf=5保证叶子节点不至于太细碎。区间合并后的长度过滤阈值 3 是我加的,作品里没提,但实际跑的时候不加这个过滤会出一堆单点区间,全是误报。

3. 前兆特征识别:B 类信号与 A 类的统计对比怎么做

3.1 前兆信号和正常信号的统计差异在哪

问题二的核心是找出前兆特征信号(B 类)和正常工作信号(A 类)的统计差异。作品给出的对比结论很明确:电磁辐射中,B 类均值 71.05 比 A 类 49.81 高,标准偏差 91.90 比 18.28 大了近 5 倍,峰度 9.24 比 23.91 小,偏度 3.19 比 3.85 略小。翻译成人话就是:前兆信号整体幅值抬升、波动剧烈、分布比正常信号更接近正态(峰度低意味着尾部没那么厚)。

声发射的对比更有意思:B 类峰度 1.43,A 类 10.39,差了 7 倍多。峰度小于 3 说明 B 类声发射信号的分布比正态还平,能量分散在更宽的幅值范围里。这个特征在实际预警里很有价值——如果你看到声发射信号的峰度突然从 10 掉到 2 以下,大概率是前兆阶段来了。

3.2 用同一套傅里叶加随机森林流程识别前兆区间

问题二的建模流程和问题一几乎一样,区别只在标签:问题一标的是 C 类(干扰),问题二标的是 B 类(前兆)。代码只需要把data['类别'] == 'C'改成data['类别'] == 'B',其余不动。这种复用是合理的,因为干扰信号和前兆信号在频域上的表现确实有相似之处——都是波动加剧、频谱能量分散。

但有一个坑要注意:B 类样本量通常比 C 类少。如果 B 类只有几百个点而总数据量几万条,随机森林的类别不平衡会导致召回率很低。我一般会加class_weight='balanced'让模型自动调整权重:

rf = RandomForestClassifier(n_estimators=150, max_depth=12, min_samples_leaf=3, class_weight='balanced', random_state=42)

min_samples_leaf从 5 降到 3 也是因为 B 类样本少,叶子节点需要更细才能捕捉到前兆信号。这两个参数调整是我踩过坑之后加的,作品原文没写,但不加的话 B 类召回率会明显偏低。

3.3 区间结果的后处理与表格填写

问题二要求输出最早 5 个前兆特征区间,格式和问题一的干扰区间表一样。后处理逻辑可以直接复用问题一的区间合并代码,只是输入换成 B 类的预测结果。作品给出的电磁辐射前兆区间从 2020-4-10 到 2022-5-25,声发射前兆区间从 2021-11-1 到 2022-10-24,跨度都很大,说明前兆事件在时间上的分布比干扰事件更分散。

这里有个填表时的细节:题目要求的是"最早发生的 5 个",所以区间合并完之后要按起始时间排序,取前 5 个。如果你的区间数量不足 5 个,说明模型召回率不够,需要回头检查特征工程或者调整随机森林参数。

4. 朴素贝叶斯算概率:从分类结果到危险预判

4.1 为什么用朴素贝叶斯而不是继续用随机森林

问题三的要求变了——不再是识别区间,而是算每个时间段最后时刻出现前兆特征的概率。这是一个典型的概率输出需求,随机森林虽然也能输出predict_proba,但它的概率估计是基于投票比例的,在样本不平衡时校准很差。朴素贝叶斯直接基于贝叶斯定理计算后验概率,对于连续特征假设高斯分布,输出的概率值更有统计意义。

作品用的公式是朴素贝叶斯的标准形式:先算先验概率 P(C_j),再算每个特征在各类别下的条件概率 P(x_r | C_j),连乘得到似然,最后归一化得到后验概率。因为特征是连续的,条件概率用高斯分布密度函数计算,需要估计每个类别下每个特征的均值和方差。

4.2 高斯朴素贝叶斯的代码实现

from sklearn.naive_bayes import GaussianNB from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report # 用附件1的数据训练,附件3的数据预测 feat_train = np.loadtxt('电磁辐射数据.txt') X_train, y_train = feat_train[:, :4], feat_train[:, 4] # 训练高斯朴素贝叶斯 gnb = GaussianNB() gnb.fit(X_train, y_train) # 对附件3的每个时间段最后时刻做预测 # 假设附件3已经做了同样的特征提取,得到 X_test X_test = np.loadtxt('附件3特征.txt') proba = gnb.predict_proba(X_test) # 输出前兆特征(B类)的概率 for i, p in enumerate(proba): print(f"时间段 {i+1} 前兆概率: {p[1]*100:.2f}%")

GaussianNB不需要调参,它的核心参数就是每个特征的均值和方差,直接从训练数据估计。但有一个前提:特征要近似服从高斯分布。如果你的特征分布严重偏斜(比如频谱方差这种平方量),建议先做对数变换再喂给模型。

4.3 概率结果的解读与阈值选择

作品算出的电磁辐射前兆概率是 61.47%、45.76%、54.40%、34.00%、46.00%,声发射是 0.006%、0.03%、99.93%、50%、58%。声发射那组数据里 99.93% 和 0.006% 两个极端值很扎眼,说明模型对某些时间段的判别非常自信,而对另一些几乎不确定。

实际用的时候,概率值本身不是决策依据,你需要定一个阈值。我一般会取 50% 作为分界线,高于 50% 判为有前兆风险,低于 50% 判为安全。但这个阈值要根据漏报和误报的代价来调——煤矿场景下漏报的代价远大于误报,所以阈值可以降到 30% 甚至更低,宁可多报几次也不要漏掉一次。

5. 避坑与排查:这套流程里最容易翻车的五个地方

5.1 数据时间戳不连续导致区间错位

现象:合并出来的区间起止时间和实际数据对不上,差了几个小时甚至几天。

原因:原始数据里存在采样间隔不等于 30 秒的记录,作品在数据预处理部分提到了要检查间隔异常,但代码里没有体现。如果直接按行号合并区间,遇到时间跳跃就会错位。

解决:在读数据之后先算时间戳差分,把间隔不等于 30 秒的行标记出来,要么插值补齐,要么在合并区间时用时间戳而不是行号做边界判断。

5.2 傅里叶变换的直流分量被均值淹没

现象:特征表里zbiao[i, 2](直流分量)和zbiao[i, 0](均值)高度相关,随机森林的特征重要性显示这两个特征几乎等价。

原因:FFT 的直流分量 a[0].real 本质上就是窗口内所有点的和,和均值只差一个常数因子。作品里同时保留了这两个特征,实际上信息是冗余的。

解决:要么去掉均值只留直流分量,要么把直流分量换成其他频段的分量(比如 a[1].real 或 a[2].real),让特征之间更有区分度。我一般会保留均值,把直流分量换成频谱中高频段的能量占比。

5.3 随机森林过拟合导致训练集精度虚高

现象:训练集上准确率 99%,测试集上掉到 70% 以下。

原因:max_depth没限制或者设得太大,树长得太深,把训练数据里的噪声也学进去了。作品在模型评价部分也提到了这个问题,说"如果没有合理的参数调整会造成过拟合风险"。

解决:max_depth控制在 10-15 之间,min_samples_leaf不低于 3,min_samples_split不低于 5。另外可以用cross_val_score做 5 折交叉验证,看验证集精度和训练集精度的差距,差距超过 10% 就是过拟合了。

5.4 朴素贝叶斯的独立性假设被严重违反

现象:概率输出全是 0% 或 100%,没有中间值。

原因:朴素贝叶斯假设特征之间相互独立,但均值、方差、频谱方差这几个特征明显是相关的。如果相关性太强,似然连乘会放大到极端值,导致概率饱和。

解决:要么做特征降维(PCA)去掉相关性,要么换用不需要独立性假设的模型(比如逻辑回归)。作品在模型缺点里也承认了这一点,说"默认各种属性之间是相互独立的,会有偏差"。

5.5 类别不平衡导致 B 类召回率极低

现象:模型把所有样本都预测为 A 类(正常工作),B 类一个都没识别出来。

原因:A 类样本量远大于 B 类,随机森林默认的基尼系数在不平衡数据上偏向多数类。

解决:加class_weight='balanced',或者用 SMOTE 做过采样。另外评估指标不要看准确率,要看召回率和 F1,准确率在不平衡数据上没有意义。

6. 进阶技巧:用滑动窗口投票提升区间识别的稳定性

前面讲的都是作品里的基础流程,实际跑的时候你会发现一个问题:逐点预测的区间边界很毛糙,经常出现一个长区间被切成好几段、或者两个相邻区间中间夹了一个孤立的反例点。我后来改成了一个滑动窗口投票的机制,效果比原始方案稳不少。

思路是这样的:不再对每个点独立预测,而是以每个点为中心取一个长度为 11 的窗口,对窗口内所有点的预测结果做多数投票,得票超过一半才判为该类。这样做的效果是区间边界更平滑,孤立误报被过滤掉。

def smooth_predictions(pred, window=11): """对逐点预测结果做滑动窗口多数投票""" half = window // 2 smoothed = np.zeros_like(pred) for i in range(len(pred)): left = max(0, i - half) right = min(len(pred), i + half + 1) votes = pred[left:right].sum() smoothed[i] = 1 if votes > (right - left) / 2 else 0 return smoothed # 在随机森林预测之后加一步平滑 pred_raw = rf.predict(X) pred_smooth = smooth_predictions(pred_raw, window=11)

窗口长度 11 对应约 5 分钟的采样跨度,这个值可以调。窗口越大,区间越平滑,但短时前兆事件可能被抹掉;窗口越小,灵敏度越高,但误报也会增加。我一般会在验证集上试 5、11、21 三个值,选 F1 最高的那个。

另一个技巧是特征重要性筛选。随机森林训练完之后可以输出feature_importances_,如果某个特征的 importance 低于 0.05,说明它对分类几乎没有贡献,可以直接去掉。我跑下来发现 4 个特征里通常只有 2-3 个是真正有用的,去掉冗余特征之后模型训练速度能快 30% 左右,精度基本不掉。

还有一个验证方法:把识别出的区间和原始数据的折线图叠在一起看。如果区间对应的时段确实有明显的幅值抬升或波动加剧,说明模型抓到了真实信号;如果区间落在平稳段上,那就是误报,需要回头检查特征工程。这个可视化验证步骤作品里没写,但我觉得比任何指标都直观。

从那以后我每次做完时序分类都会强制走一遍滑动窗口平滑和可视化验证,这两个步骤花不了多少时间,但能挡掉大部分低级错误。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询