★ FEATURED ARTICLE 简介本资源为2015年第十二届五一数学建模联赛B题优秀论文完整文档面向参加数学建模竞赛的高校学生及指导教师尤其适合研究空气污染扩散建模与评价体系的读者参考。压缩包内仅含1个doc文件约1.5MB完整收录了该获奖论文的承诺书、编号专用页、摘要、问题重述、问题分析及模型求解全过程。论文围绕京津冀地区空气污染问题依次构建了基于层次分析法的空气质量评价模型、结合因子分析与动态加权的污染源识别方法、修正高斯烟羽扩散模型以及灰色预测模型并针对单污染源与多污染源场景给出浓度分布与空气质量等级求解结果最后提出可行性治理建议。目前已有329人学习下载读者可从中获取完整的赛题解题思路、建模框架、公式推导与论文写作范式适合作为数模竞赛备赛与空气污染建模研究的参考范例。1. 空气污染问题研究一份数模赛题为什么值得用代码重做一遍2015年第十二届五一数模联赛B题给了一份空气质量监测数据要求参赛队在三天内完成污染特征分析、影响因素建模和治理建议。十年过去这类题目的价值反而更清晰了——它几乎是数据建模全流程的最小闭环数据清洗、特征工程、回归预测、结果解释一个都不少。很多做数据分析的同行第一次真正理解“缺失值处理会直接改变结论”就是在这种赛题里翻的车。如果你手头有类似的空气质量监测数据或者想找一个完整练手项目把 pandas、sklearn、matplotlib 串起来这份题目的结构比大多数教程都扎实。下面我按自己重做这类题目的实际路径把每一步拆开讲清楚。2. 拿到数据先别建模空气污染数据的清洗与特征构造2.1 先搞清楚数据长什么样再动手这类赛题的数据通常以 Excel 或 CSV 形式给出字段一般包括监测点编号、日期时间、PM2.5、PM10、SO2、NO2、CO、O3 等污染物浓度以及温度、湿度、风速、风向等气象要素。不同监测点的记录频率可能不一致有的是逐小时有的是逐日。拿到数据后的第一件事不是画图而是用几行代码把数据的基本面貌摸清楚。import pandas as pd import numpy as np # 读取数据注意编码问题中文列名常见gbk或utf-8 df pd.read_csv(air_quality.csv, encodingutf-8) # 基本信息行数、列数、每列数据类型 print(df.shape) print(df.dtypes) # 缺失情况一览 missing df.isnull().sum() missing_pct (missing / len(df) * 100).round(2) print(pd.DataFrame({缺失数: missing, 缺失占比%: missing_pct})) # 时间列解析这是后续所有时序分析的基础 df[datetime] pd.to_datetime(df[datetime], errorscoerce) df df.sort_values(datetime).reset_index(dropTrue) # 数值列描述统计快速看有没有离谱的极值 print(df.describe().T[[mean, std, min, max]])这段代码做了四件事确认数据规模、检查缺失分布、解析时间列、看数值范围。逻辑很直白但每一步都有讲究。errorscoerce是为了把无法解析的时间字符串变成 NaT 而不是直接报错方便后续统一处理。describe()里如果某个污染物的 max 值比 mean 大几十倍大概率存在异常值或者单位不统一的问题。参数方面编码格式需要根据实际文件调整。如果utf-8报错就换gbk这是中文数据集的常见坑。时间列的列名也可能是“时间”“监测时间”“date”等需要按实际字段名替换。2.2 缺失值不是填个数就完事空气污染数据的缺失有很强的规律性夜间某些监测点可能不记录设备故障会导致连续多小时缺失极端天气下传感器可能直接掉线。不同缺失模式对应不同的处理策略。# 先看缺失是不是连续的 df[pm25_missing] df[PM2.5].isnull().astype(int) # 按小时统计缺失率看是否有时间聚集性 df[hour] df[datetime].dt.hour hourly_missing df.groupby(hour)[pm25_missing].mean() print(hourly_missing) # 短缺口连续3小时用线性插值 # 长缺口连续3小时标记后不插值避免引入虚假数据 df[PM2.5_filled] df[PM2.5].interpolate(methodlinear, limit3) # 对长缺口保留NaN后续建模时用样本筛选排除 long_gap_mask df[PM2.5].isnull() df[PM2.5_filled].isnull() print(f长缺口样本数: {long_gap_mask.sum()}, 占比: {long_gap_mask.mean()*100:.1f}%)这里的核心判断是插值只在缺口很短的时候才可靠。连续缺失超过 3 小时线性插值等于在编数据。我一般会把长缺口样本单独标记建模时要么排除要么用模型预测填充并在结果中注明。limit3这个参数控制最大插值长度可以根据数据记录频率调整——如果是逐小时数据3 小时以内的缺口插值还算合理如果是逐日数据limit 应该设为 1。2.3 特征构造决定模型上限原始字段直接扔进模型也能跑但效果通常差一截。空气污染有明显的时间周期性和气象依赖性把这些先验知识变成特征比换模型管用得多。# 时间特征 df[month] df[datetime].dt.month df[dayofweek] df[datetime].dt.dayofweek df[is_weekend] (df[dayofweek] 5).astype(int) df[hour_sin] np.sin(2 * np.pi * df[hour] / 24) df[hour_cos] np.cos(2 * np.pi * df[hour] / 24) # 滞后特征前一小时、前24小时的污染物浓度 for lag in [1, 3, 24]: df[fPM2.5_lag{lag}] df[PM2.5_filled].shift(lag) # 滚动统计过去6小时均值 df[PM2.5_roll6_mean] df[PM2.5_filled].rolling(window6, min_periods3).mean() df[PM2.5_roll6_std] df[PM2.5_filled].rolling(window6, min_periods3).std() # 气象交互特征风速低湿度高是污染累积的典型条件 df[wind_humidity] df[wind_speed] / (df[humidity] 1)时间特征用 sin/cos 编码是为了保留周期性——直接给 0-23 的整数模型会认为 23 点和 0 点差很远实际上它们相邻。滞后特征和滚动统计是时序预测的核心shift(1)表示用前一小时的值预测当前值这在做预测任务时必须注意不能引入未来信息。min_periods3保证滚动窗口内至少有三个有效值才计算避免前几行全是 NaN。3. 建模路线怎么选从线性回归到树模型的取舍3.1 先跑一个线性基线别一上来就上深度学习很多同行拿到数据直接上 LSTM 或者 Transformer结果发现效果还不如线性回归。原因很简单数据量不够、特征工程没做到位、过拟合严重。我的习惯是先跑一个带正则的线性模型作为基线确认特征方向对不对再决定要不要上复杂模型。from sklearn.linear_model import Ridge from sklearn.model_selection import TimeSeriesSplit from sklearn.preprocessing import StandardScaler from sklearn.metrics import mean_absolute_error, r2_score from sklearn.pipeline import Pipeline # 构造特征矩阵和目标变量 feature_cols [hour_sin, hour_cos, month, is_weekend, PM2.5_lag1, PM2.5_lag3, PM2.5_lag24, PM2.5_roll6_mean, PM2.5_roll6_std, wind_humidity, temperature, wind_speed, humidity] # 排除长缺口样本 model_df df.dropna(subsetfeature_cols [PM2.5_filled]).copy() X model_df[feature_cols].values y model_df[PM2.5_filled].values # 时序交叉验证不能用随机KFold tscv TimeSeriesSplit(n_splits5) ridge_pipe Pipeline([ (scaler, StandardScaler()), (ridge, Ridge(alpha1.0)) ]) scores [] for train_idx, test_idx in tscv.split(X): X_train, X_test X[train_idx], X[test_idx] y_train, y_test y[train_idx], y[test_idx] ridge_pipe.fit(X_train, y_train) pred ridge_pipe.predict(X_test) scores.append({ MAE: mean_absolute_error(y_test, pred), R2: r2_score(y_test, pred) }) import pandas as pd print(pd.DataFrame(scores).mean())这里最关键的选择是TimeSeriesSplit而不是KFold。时序数据用随机切分会导致未来信息泄露——训练集里混入了测试集之后的数据评估结果虚高。TimeSeriesSplit保证每次训练集都在测试集之前模拟真实预测场景。alpha1.0是 Ridge 的正则强度值越大对系数压缩越狠如果特征之间共线性严重可以适当调大。3.2 树模型什么时候比线性模型强线性模型假设特征和目标之间是线性关系但空气污染和气象因素之间明显不是——风速对污染物浓度的影响在低风速区间很敏感高风速区间就趋于平缓。这种非线性关系用树模型能自动捕捉。import lightgbm as lgb from sklearn.model_selection import TimeSeriesSplit # 用同样的特征矩阵 lgb_params { objective: regression, metric: mae, learning_rate: 0.05, num_leaves: 31, max_depth: 6, min_child_samples: 20, subsample: 0.8, colsample_bytree: 0.8, reg_alpha: 0.1, reg_lambda: 0.1, verbose: -1 } tscv TimeSeriesSplit(n_splits5) lgb_scores [] for train_idx, test_idx in tscv.split(X): X_train, X_test X[train_idx], X[test_idx] y_train, y_test y[train_idx], y[test_idx] dtrain lgb.Dataset(X_train, labely_train) dval lgb.Dataset(X_test, labely_test, referencedtrain) model lgb.train( lgb_params, dtrain, num_boost_round500, valid_sets[dval], callbacks[lgb.early_stopping(30), lgb.log_evaluation(0)] ) pred model.predict(X_test, num_iterationmodel.best_iteration) lgb_scores.append({ MAE: mean_absolute_error(y_test, pred), R2: r2_score(y_test, pred), best_iter: model.best_iteration }) print(pd.DataFrame(lgb_scores).mean())参数里几个关键项learning_rate0.05配合early_stopping(30)是防止过拟合的标准组合学习率低就需要更多轮次但泛化更好。num_leaves31控制树的复杂度值越大模型越容易记住训练数据。min_child_samples20保证每个叶子节点至少有 20 个样本对小数据集尤其重要。subsample和colsample_bytree都是行采样和列采样增加模型多样性。如果 LightGBM 的 MAE 比 Ridge 低 10% 以上说明数据里确实有非线性关系值得挖掘。如果差距很小直接用 Ridge 更省事解释性也更好。3.3 特征重要性怎么看才不误导树模型可以输出特征重要性但默认的 split 次数统计会偏向高基数特征。更可靠的方式是看 gain 或者用 permutation importance。# 用最后一折的模型看特征重要性 importance pd.DataFrame({ feature: feature_cols, gain: model.feature_importance(importance_typegain), split: model.feature_importance(importance_typesplit) }).sort_values(gain, ascendingFalse) print(importance.head(10)) # permutation importance 更稳健 from sklearn.inspection import permutation_importance # 需要把lgb模型包装成sklearn接口或者手动实现 # 这里用最后一折的测试集做permutation result permutation_importance( model, X_test, y_test, n_repeats10, random_state42, scoringneg_mean_absolute_error ) perm_imp pd.DataFrame({ feature: feature_cols, importance_mean: result.importances_mean, importance_std: result.importances_std }).sort_values(importance_mean, ascendingFalse) print(perm_imp.head(10))gain 重要性统计的是特征在所有分裂中带来的损失下降总和比 split 次数更能反映实际贡献。permutation importance 则是打乱某个特征的值看模型效果下降多少更接近“这个特征对预测有多重要”的直觉。两种方法结合看如果某个特征在两种方法里都排前列那它的重要性就比较可信。4. 避坑与排查空气污染建模里最容易翻车的五个地方4.1 用随机切分做时序验证指标虚高现象交叉验证 R² 达到 0.95上线后实际预测误差翻倍。原因train_test_split或KFold默认随机打乱训练集里混入了测试集时间点之后的数据。空气污染有强自相关性今天的数据和明天的数据高度相关模型相当于“偷看”了答案。解决所有时序建模必须用TimeSeriesSplit或者手动按时间切分。如果数据有多个监测点还要保证同一时间点的数据不会同时出现在训练集和测试集。4.2 缺失值插值引入虚假规律现象插值后模型在缺失段附近的预测异常准确但实际部署时这些时段误差很大。原因线性插值在连续缺失段会生成一条直线模型学到了这条直线的模式但真实数据并不是直线变化。解决限制插值长度超过阈值的缺口保留 NaN 并在建模时排除。如果必须填充用模型预测填充比线性插值更合理但要在结果中注明填充比例。4.3 特征泄露用了未来才知道的信息现象模型在训练集和测试集上表现都很好但实际预测时完全不可用。原因构造特征时不小心用了未来数据。比如计算“当日均值”时用了全天数据来预测当天上午的浓度或者滚动窗口没有加shift。解决构造任何统计特征时问自己一句“这个值在预测时刻真的能拿到吗”。滚动统计要shift(1)后再 rolling确保只用历史信息。4.4 气象特征单位不统一现象风速特征的系数方向反了风速越大污染越重。原因不同监测点的风速单位可能不一致有的用 m/s有的用 km/h甚至有的用风力等级。单位不统一导致模型学到的关系混乱。解决建模前统一所有物理量的单位风速统一转 m/s温度统一摄氏度浓度统一 μg/m³。转换公式要写在数据预处理脚本里不要手动改。4.5 过拟合赛题评分标准现象在本地验证集上调到最优提交后排名靠后。原因赛题评分通常用独立测试集和本地验证集分布不同。过度调参会让模型在本地验证集上过拟合。解决保留一个完全独立的验证集只在最后评估时用一次。调参用交叉验证但不要反复用同一个验证集调参。模型复杂度控制在合理范围树模型深度不要超过 8 层。5. 从赛题到落地把分析结果变成可复用的监测工具5.1 用滚动预测模拟真实部署赛题只要求给出分析结论但实际工作中需要的是持续预测能力。滚动预测是检验模型实用性的最好方式每次只用历史数据预测下一时刻然后把真实值加入历史继续预测。def rolling_forecast(model, df, feature_cols, target_col, start_idx, steps24): 滚动预测每次预测一个时间点用真实值更新历史 start_idx: 从第几个样本开始预测 steps: 预测多少步 predictions [] actuals [] for i in range(start_idx, min(start_idx steps, len(df))): # 用当前样本之前的全部数据构造特征 history df.iloc[:i].copy() # 重新计算滞后和滚动特征实际部署时这些是实时更新的 for lag in [1, 3, 24]: history[fPM2.5_lag{lag}] history[PM2.5_filled].shift(lag) history[PM2.5_roll6_mean] history[PM2.5_filled].rolling(6, min_periods3).mean() history[PM2.5_roll6_std] history[PM2.5_filled].rolling(6, min_periods3).std() # 取最后一行作为预测输入 last_row history[feature_cols].iloc[-1:].values if np.isnan(last_row).any(): continue pred model.predict(last_row)[0] predictions.append(pred) actuals.append(df[target_col].iloc[i]) mae mean_absolute_error(actuals, predictions) return predictions, actuals, mae # 用最后24小时做滚动预测 preds, acts, mae rolling_forecast(model, model_df, feature_cols, PM2.5_filled, start_idxlen(model_df)-24, steps24) print(f滚动预测MAE: {mae:.2f})这段代码的关键在于每次预测都重新计算特征模拟真实场景中数据逐步到达的过程。如果滚动预测的 MAE 比交叉验证的 MAE 高很多说明模型对历史数据的依赖过强实际部署时效果会打折扣。5.2 阈值预警比精确预测更实用实际监测场景中决策者关心的不是“明天 PM2.5 是多少”而是“明天会不会超标”。把回归问题转成分类问题用预测值是否超过阈值来评估往往更有业务价值。from sklearn.metrics import classification_report, confusion_matrix # 以75μg/m³为阈值常见日均二级标准 threshold 75 # 用滚动预测结果做分类评估 pred_labels (np.array(preds) threshold).astype(int) actual_labels (np.array(acts) threshold).astype(int) print(confusion_matrix(actual_labels, pred_labels)) print(classification_report(actual_labels, pred_labels, target_names[未超标, 超标]))分类评估的好处是直接对应业务动作预测超标就触发预警预测未超标就正常监测。阈值可以根据实际标准调整不同污染物用不同阈值。如果漏报实际超标但预测未超标代价高可以适当降低分类阈值牺牲一点精确率换召回率。5.3 模型更新频率怎么定空气污染的影响因素会随季节变化冬季取暖排放增加夏季光化学反应增强。一个在冬季数据上训练的模型到夏季可能完全失效。我的习惯是用最近 3 个月的数据训练模型每两周重新训练一次。如果数据量足够可以按季节分别建模。验证方式是看滚动预测的误差是否随时间增大如果连续一周 MAE 上升超过 20%就该重新训练了。# 简单的模型衰减监控 def monitor_model_decay(model, df, feature_cols, target_col, window7): 监控最近window天的预测误差判断是否需要重训 recent_mae [] for i in range(len(df) - window, len(df)): history df.iloc[:i].copy() for lag in [1, 3, 24]: history[fPM2.5_lag{lag}] history[PM2.5_filled].shift(lag) history[PM2.5_roll6_mean] history[PM2.5_filled].rolling(6, min_periods3).mean() last_row history[feature_cols].iloc[-1:].values if np.isnan(last_row).any(): continue pred model.predict(last_row)[0] recent_mae.append(abs(pred - df[target_col].iloc[i])) avg_mae np.mean(recent_mae) print(f最近{window}天平均MAE: {avg_mae:.2f}) return avg_mae # 如果avg_mae比训练时高50%以上触发重训这套监控逻辑不复杂但能避免模型“悄悄失效”。很多团队模型上线后就不管了等到业务方反馈预测不准才去排查中间可能已经积累了几周的误差。5.4 一个容易被忽略的技巧分位数回归看极端值普通回归预测的是条件均值但空气污染治理最关心的是极端高值。用分位数回归可以给出不同分位数的预测区间比如 90 分位数预测值代表“有 90% 把握不超过这个浓度”。import lightgbm as lgb # 训练一个90分位数的模型 lgb_quantile_params { objective: quantile, alpha: 0.9, # 90分位数 metric: quantile, learning_rate: 0.05, num_leaves: 31, verbose: -1 } dtrain_q lgb.Dataset(X_train, labely_train) model_q lgb.train(lgb_quantile_params, dtrain_q, num_boost_round300) # 预测90分位数 pred_q90 model_q.predict(X_test) # 对比均值预测 pred_mean model.predict(X_test) # 看极端高值样本的覆盖情况 high_pollution_mask y_test np.percentile(y_test, 90) print(f极端高值样本中90分位预测覆盖比例: f{(pred_q90[high_pollution_mask] y_test[high_pollution_mask]).mean():.2f})分位数回归的alpha参数控制分位数0.9 就是 90 分位。这个模型预测值通常比均值模型高因为要覆盖更多极端情况。实际使用时可以同时部署均值模型和分位数模型均值用于日常播报分位数用于预警决策。这套流程从数据清洗到滚动预测再到模型监控基本覆盖了空气污染数据分析的完整链路。我自己的习惯是每接一个新数据集先把清洗和特征构造的代码写成函数后面换数据只需要改列名映射。模型部分保持简单优先用 Ridge 和 LightGBM 两个基线效果不够再考虑更复杂的方案。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
