煤矿冲击地压预测：微震数据傅里叶变换与随机森林建模实战

简介这份资源是2024年五一数学建模竞赛C题「煤矿深部开采冲击地压危险预测研究」的获奖作品面向参加数学建模竞赛的学生及从事矿山安全数据分析的研究者帮助解决电磁辐射与声发射信号中干扰识别、前兆特征提取及危险概率预判等问题。压缩包内共1个PDF文件约1.16MB完整呈现赛题论文涵盖傅里叶变换、随机森林与朴素贝叶斯模型的建模思路与结果。目前已有1347人学习下载。读者可从中获取问题重述、模型假设、数据预处理与特征提取的完整流程包括干扰信号统计特征对比、最早5个干扰区间识别、前兆特征趋势分析以及各时间段前兆出现概率的计算结果适合作为竞赛复盘、论文写作与算法迁移的参考范例。1. 煤矿深部开采冲击地压危险预测从微震数据到可复现的建模链路2024年五一赛C题把「煤矿深部开采冲击地压危险预测」摆到了台面上这道题的核心不是让你背公式而是逼你把微震监测数据、声发射信号、应力时序这些真实工况下的多源数据用傅里叶变换做频域特征提取再用随机森林和朴素贝叶斯做危险等级分类。我见过太多队伍卡在第一步——数据读进来发现采样频率不统一、标签分布极度不均衡、特征量纲差异巨大然后直接上模型结果AUC看着还行但召回率惨不忍睹。这篇笔记不讲虚的按我实际做这类时序分类题的路子从数据预处理、频域特征构造、模型选型到调参排坑一步步拆开讲。适合正在做数学建模竞赛、或者手头有工业监测数据想做危险预测的从业者新手能跟着代码跑通熟手能直接拿去改参数上自己的数据集。2. 冲击地压数据到底长什么样先搞清楚你在预测什么2.1 微震监测数据的字段结构与标签定义煤矿冲击地压预测的原始数据通常来自微震监测系统每条记录至少包含时间戳、震源坐标x, y, z、能量值、震级、以及是否发生冲击地压事件的标签。五一赛C题给的数据集我印象里是分了训练集和测试集训练集有明确标签测试集需要你预测危险等级。这里第一个坑就是标签的定义方式——有的队伍直接把「是否发生冲击地压」当二分类有的按能量阈值分三级低危、中危、高危这两种做法对后续模型选择影响很大。我一般会先做一件事把标签分布画出来。如果高危样本占比不到5%那朴素贝叶斯基本不用考虑了它的先验概率会被多数类主导预测出来全是低危。随机森林虽然对不平衡数据有一定鲁棒性但也要配合类别权重或者过采样。具体操作上用pandas读进来先看value_counts()再决定要不要做SMOTE或者调整class_weight参数。import pandas as pd import numpy as np # 读取训练数据假设是CSV格式 df pd.read_csv(train_data.csv) # 查看标签分布 label_counts df[risk_level].value_counts() print(label_counts) print(f高危样本占比: {label_counts.get(high, 0) / len(df):.4f}) # 如果极度不平衡考虑加权 from sklearn.utils.class_weight import compute_class_weight classes np.unique(df[risk_level]) weights compute_class_weight(balanced, classesclasses, ydf[risk_level]) class_weight_dict dict(zip(classes, weights)) print(class_weight_dict)这段代码的逻辑是先摸清标签分布再计算类别权重。compute_class_weight的balanced模式会根据样本频率自动反比赋权后续传给随机森林的class_weight参数即可。注意如果标签是字符串要先做LabelEncoder转成数值不然sklearn会报错。2.2 采样频率不统一怎么对齐重采样与时间窗划分微震数据另一个恶心的地方是采样频率可能不一致。有的传感器是100Hz有的是200Hz直接拼在一起做傅里叶变换频率轴对不上特征全是错的。我一般会统一重采样到最低频率或者用时间窗切片每个窗口内做统计特征。窗口长度选多少冲击地压的前兆信号通常持续几秒到几十秒我一般用10秒窗、5秒滑移这样既能捕捉短时能量突变又不会因为窗口太长丢失细节。重采样用scipy.signal.resample或者pandas的resample都行但要注意如果原信号有高频成分降采样前必须加抗混叠滤波器不然频谱会混叠傅里叶变换出来的高频特征全是假的。这个坑我踩过当时模型训练集AUC 0.95测试集直接掉到0.6查了半天才发现是混叠。from scipy.signal import resample, butter, filtfilt def resample_signal(signal, original_fs, target_fs): # 先做抗混叠滤波 nyquist target_fs / 2 b, a butter(4, nyquist / (original_fs / 2), btypelow) filtered filtfilt(b, a, signal) # 再重采样 num_samples int(len(signal) * target_fs / original_fs) resampled resample(filtered, num_samples) return resampled # 假设原始信号1000Hz目标200Hz original_signal df[energy].values[:1000] resampled_signal resample_signal(original_signal, 1000, 200) print(f重采样后长度: {len(resampled_signal)})butter(4, ...)里的4是滤波器阶数阶数越高过渡带越陡但相位失真越大一般4到6够用。filtfilt是零相位滤波避免信号延迟。重采样后的长度按比例算注意取整时可能差一两个点后续做FFT时要统一长度。3. 傅里叶变换提取频域特征不只是调用fft那么简单3.1 从时域到频域哪些特征真正对冲击地压敏感傅里叶变换在这道题里的作用是把微震信号的时域波形转成频域谱然后提取主频、频谱质心、频带能量比这些特征。但不是什么频域特征都有用。冲击地压的前兆通常表现为低频能量增加、主频向低频偏移所以我会重点算三个指标主频频谱最大幅值对应的频率、低频段0-5Hz能量占总能量的比例、以及频谱质心。这三个特征物理意义明确而且对噪声没那么敏感。具体实现上用numpy.fft.fft做完变换后频率轴是np.fft.fftfreq(n, d1/fs)取前半段正频率部分。幅值谱取绝对值再归一化。注意如果信号长度不是2的幂FFT效率会低但精度没问题不用强行补零到2的幂除非你做的是实时系统对速度有要求。import numpy as np def extract_freq_features(signal, fs): n len(signal) # 去均值避免直流分量干扰 signal signal - np.mean(signal) # FFT fft_vals np.fft.fft(signal) fft_freqs np.fft.fftfreq(n, d1/fs) # 取正频率部分 positive_mask fft_freqs 0 freqs fft_freqs[positive_mask] magnitudes np.abs(fft_vals[positive_mask]) / n # 主频 dominant_freq freqs[np.argmax(magnitudes)] # 频谱质心 if np.sum(magnitudes) 0: spectral_centroid np.sum(freqs * magnitudes) / np.sum(magnitudes) else: spectral_centroid 0 # 低频能量比0-5Hz low_freq_mask freqs 5 low_energy_ratio np.sum(magnitudes[low_freq_mask]**2) / np.sum(magnitudes**2) return { dominant_freq: dominant_freq, spectral_centroid: spectral_centroid, low_energy_ratio: low_energy_ratio } # 对每个时间窗提取特征 features extract_freq_features(resampled_signal, fs200) print(features)这段代码的关键点是去均值那一步很多人忘了做结果频谱在0Hz处有个巨大峰值主频永远算出来是0。magnitudes除以n是归一化不影响主频和质心的相对值但影响能量比的绝对值。低频段选0-5Hz是根据冲击地压信号的先验知识如果你做的是其他场景这个阈值要调整。3.2 短时傅里叶变换处理非平稳信号的正确姿势微震信号是非平稳的直接做全局FFT会丢失时间信息。比如一个信号前5秒正常、后5秒能量突增全局FFT只能看到整体频谱看不出突变发生在什么时候。这时候要用短时傅里叶变换STFT加滑动窗每个窗内做FFT得到时频图。从时频图里可以提取「能量突增时刻」「频率随时间的变化率」这些特征对冲击地压预测非常有用。STFT的窗长和重叠率是两个关键参数。窗长太短频率分辨率不够窗长太长时间分辨率不够。我一般用窗长256点、重叠128点对应200Hz采样率就是1.28秒窗、0.64秒滑移。这个配置在冲击地压数据上表现比较稳。scipy.signal.stft直接返回时频矩阵取模平方得到功率谱密度。from scipy.signal import stft def extract_stft_features(signal, fs): f, t, Zxx stft(signal, fsfs, nperseg256, noverlap128) # 功率谱密度 psd np.abs(Zxx)**2 # 能量突增检测计算每个时间帧的总能量 frame_energy np.sum(psd, axis0) # 能量变化率 energy_diff np.diff(frame_energy) max_energy_jump np.max(np.abs(energy_diff)) if len(energy_diff) 0 else 0 # 频率随时间的变化主频轨迹的方差 dominant_freqs f[np.argmax(psd, axis0)] freq_variance np.var(dominant_freqs) return { max_energy_jump: max_energy_jump, freq_variance: freq_variance, mean_frame_energy: np.mean(frame_energy) } stft_features extract_stft_features(resampled_signal, fs200) print(stft_features)nperseg256是窗长noverlap128是重叠点数。np.argmax(psd, axis0)沿着频率轴找最大值索引得到每个时间帧的主频。freq_variance大说明频率不稳定可能是冲击地压前兆。max_energy_jump直接反映能量突变这个特征在随机森林里重要性通常排前三。4. 随机森林与朴素贝叶斯的选型与调参别拿一个模型硬套4.1 随机森林在冲击地压分类中的参数怎么调随机森林是我做这类工业数据分类的首选原因很简单它对特征量纲不敏感、能处理非线性关系、自带特征重要性输出、对不平衡数据有一定容忍度。但默认参数往往不是最优的几个关键参数必须调n_estimators树的数量、max_depth最大深度、min_samples_split节点分裂最小样本数、class_weight类别权重。我的调参顺序是先固定n_estimators100跑一版看基线然后调max_depth防止过拟合再调min_samples_split控制分裂粒度最后加class_weight处理不平衡。n_estimators一般200到500够用再多边际收益很低。max_depth我一般从10开始试如果训练集AUC远高于测试集就往下调。min_samples_split默认是2太容易过拟合我通常设5到20之间。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import GridSearchCV from sklearn.metrics import classification_report # 假设X_train, y_train已经准备好 rf RandomForestClassifier(random_state42, n_jobs-1) param_grid { n_estimators: [100, 200, 300], max_depth: [8, 12, 16, None], min_samples_split: [5, 10, 20], class_weight: [balanced, None] } grid_search GridSearchCV( estimatorrf, param_gridparam_grid, cv5, scoringf1_macro, n_jobs-1, verbose1 ) grid_search.fit(X_train, y_train) print(f最佳参数: {grid_search.best_params_}) print(f最佳F1: {grid_search.best_score_:.4f}) # 用最佳模型预测 best_rf grid_search.best_estimator_ y_pred best_rf.predict(X_test) print(classification_report(y_test, y_pred))scoringf1_macro是因为类别不平衡时准确率没意义F1的宏平均能反映每个类的表现。cv5是五折交叉验证如果数据量小可以改成10折。n_jobs-1用满所有CPU核心。注意class_weightbalanced和None都试一下有时候加权反而降低整体性能因为多数类被压制太狠。4.2 朴素贝叶斯什么时候能用、什么时候别碰朴素贝叶斯在这道题里不是不能用但限制很多。它的核心假设是特征条件独立而傅里叶变换提取的频域特征之间往往有相关性比如主频和频谱质心高度相关违反假设会导致概率估计偏差。另外朴素贝叶斯对连续特征通常假设高斯分布如果你的特征不服从正态分布效果会很差。我一般只在两种情况下用朴素贝叶斯一是特征维度很高但样本量小随机森林容易过拟合二是需要极快的推理速度比如实时监测系统。如果要用先做特征独立性检验把相关性高的特征剔除。另外用GaussianNB之前最好对特征做正态性检验偏度太大的做Box-Cox变换。from sklearn.naive_bayes import GaussianNB from scipy import stats # 正态性检验 for col in X_train.columns: stat, p_value stats.normaltest(X_train[col]) if p_value 0.05: print(f{col} 不服从正态分布p{p_value:.4f}) # 如果大部分特征不正态考虑先做变换 from sklearn.preprocessing import PowerTransformer pt PowerTransformer(methodbox-cox) # 注意Box-Cox要求数据为正 X_train_positive X_train - X_train.min() 1 X_train_transformed pt.fit_transform(X_train_positive) gnb GaussianNB() gnb.fit(X_train_transformed, y_train) y_pred_nb gnb.predict(pt.transform(X_test - X_train.min() 1))normaltest返回的p值小于0.05就拒绝正态假设。PowerTransformer的Box-Cox变换要求数据严格为正所以先减去最小值再加1。变换后的数据再喂给GaussianNB效果通常比原始数据好一截。但即便如此朴素贝叶斯在这道题上的F1通常比随机森林低5到10个百分点除非你的特征工程做得特别好。5. 避坑与排查那些让模型翻车的细节5.1 数据泄漏你的交叉验证可能一直在骗你现象交叉验证F1 0.92测试集F1 0.65差距巨大。原因做特征工程时用了全局数据的信息比如用整个数据集的均值和方差做标准化或者时间窗切片时训练集和测试集有重叠。解决标准化只能在训练集上fit然后transform测试集时间序列数据必须按时间顺序切分不能随机打乱。from sklearn.preprocessing import StandardScaler # 错误做法全局标准化 # scaler StandardScaler() # X_scaled scaler.fit_transform(X) # 泄漏 # 正确做法训练集fit测试集transform scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test)时间序列还要注意用TimeSeriesSplit代替KFold避免未来信息泄漏到过去。5.2 标签噪声微震数据里的误标注怎么处理现象模型在训练集上怎么调都上不去或者某些样本的预测概率一直在0.5附近。原因微震监测系统有误报有些标注为「高危」的样本实际能量并不高。解决用置信学习Confident Learning找出疑似误标注样本或者用鲁棒损失函数如Huber损失降低异常样本影响。我一般先画一个能量值vs标签的箱线图如果高危样本的能量分布和低危大量重叠说明标签有问题。5.3 特征量纲差异随机森林不在乎但朴素贝叶斯很在乎现象随机森林表现正常朴素贝叶斯F1只有0.3。原因频域特征里主频可能是几十Hz能量比是0到1之间的小数量纲差几个数量级高斯朴素贝叶斯的概率密度估计被大量纲特征主导。解决对朴素贝叶斯必须做标准化随机森林可以不做但做了也没坏处。5.4 过采样导致过拟合SMOTE不是万能药现象用了SMOTE之后训练集F1 0.98测试集F1 0.55。原因SMOTE在少数类样本之间插值生成新样本如果少数类样本本身噪声大插值出来的样本更离谱。解决改用SMOTEENN或SMOTETomek先过采样再清理边界噪声样本或者直接用class_weight代替过采样。from imblearn.combine import SMOTEENN smote_enn SMOTEENN(random_state42) X_resampled, y_resampled smote_enn.fit_resample(X_train, y_train) print(f过采样后样本数: {len(X_resampled)})5.5 傅里叶变换的频谱泄漏窗函数选错全盘皆输现象频谱上出现很多不该有的高频分量主频判断错误。原因信号截断时不是整数个周期导致频谱泄漏。解决加窗函数汉宁窗、汉明窗再做FFT。scipy.signal.get_window可以生成各种窗。from scipy.signal import get_window def fft_with_window(signal, fs, window_typehann): n len(signal) window get_window(window_type, n) windowed_signal signal * window fft_vals np.fft.fft(windowed_signal) freqs np.fft.fftfreq(n, d1/fs) return freqs[:n//2], np.abs(fft_vals[:n//2])汉宁窗适合大多数场景如果主频很接近0Hz用矩形窗反而更好因为汉宁窗会压制低频。6. 把模型推到能用的程度后处理与阈值优化模型训练完输出的是概率但实际预警需要的是明确的危险等级。默认的predict用0.5做阈值但在不平衡数据上0.5往往不是最优的。我一般会画PR曲线找F1最大的阈值点或者根据业务需求定——比如宁可误报不可漏报那就把高危类的阈值调低。from sklearn.metrics import precision_recall_curve y_proba best_rf.predict_proba(X_test)[:, 1] # 假设高危是正类 precision, recall, thresholds precision_recall_curve(y_test, y_proba) f1_scores 2 * precision * recall / (precision recall 1e-8) best_threshold thresholds[np.argmax(f1_scores)] print(f最佳阈值: {best_threshold:.4f}) # 用最佳阈值做预测 y_pred_optimized (y_proba best_threshold).astype(int)另外时间序列预测可以做滑动平均后处理把连续几个窗口都预测为高危才触发警报降低误报率。这个后处理窗口长度根据实际监测系统的响应时间定我一般用3到5个窗口。最后说一个我自己的习惯每次跑完模型把特征重要性排序打印出来如果排第一的特征是「样本ID」或者某个明显无关的列说明数据泄漏了。这个检查花不了几秒钟但能救你一命。希望帮到你。本文还有配套的精品资源点击获取