简介这份资源是2015年第十二届五一数学建模联赛B题“空气污染问题研究”的优秀论文文档面向备战数学建模竞赛的本科生及指导教师可用于赛题复盘、论文写作参考与建模思路学习。压缩包内仅1个doc文件约1.5MB完整收录承诺书、摘要、问题重述、问题分析、模型建立与求解等论文主体结构。论文围绕京津冀地区空气污染展开依次完成国标与美标空气质量指数对比、层次分析法赋权与动态加权综合评价、修正高斯烟羽扩散模型、多污染源扩散及灰色预测模型求解并给出“APEC蓝天”治理建议报告可帮助读者掌握层次分析法、变权函数、优化高斯烟羽模型等方法的实际应用。目前已有328人学习下载适合希望借鉴优秀论文框架、学习建模论文写作规范与模型推导细节的参赛者。1. 空气污染问题研究从 AQI 数据到可复现的建模流程拿到这个标题的人手里通常已经有一份城市小时级空气质量数据却在第一步就卡住先算 AQI还是先做预测。空气污染问题研究的落点一般分三层——把六项污染物折算成可比的分指数回答「今天空气算几级」用时间序列加气象因子回答「明天 PM2.5 会不会破 150」用扩散模型反推「污染从哪个方向来」。三层要的数据粒度、模型类型和评估口径完全不同混在一起做最后容易得到一个 R² 很高但没人能解释的模型。这套流程适合做数模、做环境数据看板、做城市空气质量分析的从业者也适合第一次接触时空数据的新手。下面按数据准备、时空特征、预测建模、扩散溯源的顺序推进每一步都给能直接跑的代码和能改的参数。2. 空气污染数据的清洗与 AQI 折算2.1 六项污染物的分指数与限值对照AQI 不是对浓度直接排序得到的而是先把每项污染物单独折算成 IAQI再取最大值。折算采用分段线性插值区间断点由《环境空气质量标准》二级限值确定。污染物平均时间二级限值IAQI100 对应浓度PM2.524 小时75 μg/m³75PM1024 小时150 μg/m³150SO₂24 小时150 μg/m³150NO₂24 小时80 μg/m³80CO24 小时4 mg/m³4O₃8 小时滑动160 μg/m³160IAQI 区间级别类别050一级优51100二级良101150三级轻度污染151200四级中度污染201300五级重度污染300六级严重污染单项分指数按 IAQI (IAQI_hi − IAQI_lo)/(BP_hi − BP_lo) × (C − BP_lo) IAQI_lo 计算其中 BP_hi、BP_lo 是浓度 C 落入区间的上下断点。最终 AQI 取六项 IAQI 的最大值同时把取最大值的那项记为首要污染物。这一步在竞赛题里经常被简化成只算 PM2.5一旦题目要求评价城市整体空气状况漏掉 O₃ 就会在夏季月份给出完全相反的结论。2.2 用 pandas 完成缺失值、哨兵值与异常值处理import pandas as pd import numpy as np df pd.read_csv(air_2015_hourly.csv, parse_dates[time]) df df.set_index(time).sort_index() # asfreq 补齐缺失时刻为 NaN比 resample 更诚实不会凭空聚出平均值 df df.asfreq(h) cols [pm25, pm10, so2, no2, o3, co] # co 单位 mg/m3其余 μg/m3 # 负值和明显超量程的哨兵值-999、9999先置空否则会污染整条序列 df[cols] df[cols].mask((df[cols] 0) | (df[cols] 2000)) # 连续缺口不超过 3 小时才线性插值长缺口保留 NaN df[cols] df[cols].interpolate(methodlinear, limit3, limit_areainside) q1, q3 df[pm25].quantile([0.25, 0.75]) iqr q3 - q1 df[pm25_outlier] (df[pm25] q1 - 3 * iqr) | (df[pm25] q3 3 * iqr)逻辑上先补时间轴再补数值asfreq 把稀疏的观测时刻拉成严格小时序列后面所有 shift、rolling 才有正确的时间含义。limit3 是个经验值超过三小时的缺口用插值填平会制造出虚假的平滑趋势模型会把这部分当成真实规律学走。异常判定用 3 倍四分位距而不是常见的 1.5 倍原因是污染物浓度分布强烈右偏沙尘过程和春节烟花时段的真实峰值会被 1.5 倍规则误杀这部分样本恰恰是预测模型最该学会的部分。注意AQI 发布口径使用的是 24 小时滑动平均做短期预测时如果直接拿实时浓度训练评估结果和官方播报值对不上两者必须区分开。2.3 风向风速的向量化与静风样本取舍风向是角度量直接求平均会出大错350° 和 10° 的算术平均是 180°方向完全反了。正确做法是先转成 u、v 分量在分量上做任何线性运算需要角度时再反算回去。# wd 为风的来向单位度ws 单位 m/s rad np.deg2rad(df[wd]) df[u] -df[ws] * np.sin(rad) # 气象约定的 u 分量指向东 df[v] -df[ws] * np.cos(rad) # v 分量指向北 # 静风时风向没有物理意义置零并单独打标 df[calm] df[ws] 0.5 df.loc[df[calm], [u, v]] 0.0静风样本占比在一份年度数据里常能到 10% 以上。这些时刻污染物更接近累积状态风向信息无效直接参与训练会让模型学到噪声。把 calm 当作哑变量单独进模型或者在统计分析阶段剔除效果都比强行插值风向好。3. 空气质量的时空特征分析与插值3.1 STL 分解看趋势项与日周期项from statsmodels.tsa.seasonal import STL s df[pm25].interpolate() # 分解前必须无缺失 res STL(s, period24, seasonal7, robustTrue).fit() res.trend.to_frame(trend).join(res.seasonal.to_frame(seasonal)).to_csv(stl_out.csv)period 取 24 表示按小时数据剥离日周期如果研究的是年度规律先按天聚合再取 period7 剥离周内效应更合适。seasonal 控制季节项的平滑程度必须是奇数7 对应三个点左右25 对应十几个点取值越大抽出的日周期越干净但会吃掉部分真实波动。robustTrue 用迭代重加权最小二乘拟合单个重污染日不会把整条趋势线拉偏。分解结果读法趋势项代表排放量与区域输送的长期变化量级通常在日均值的 ±20% 内季节项反映早晚高峰的双峰结构峰值一般出现在 79 点和 1922 点残差项才是真正难以解释的部分如果残差的标准差超过日均值的 40%说明还有重要变量没进模型优先检查逆温层高度和边界层风速。3.2 IDW 与普通克里金在站点插值上的取舍from pykrige.ok import OrdinaryKriging import numpy as np ok OrdinaryKriging(lon, lat, pm25_mean, variogram_modelspherical, nlags6, weightTrue) z, ss ok.execute(grid, grid_lon, grid_lat) # z 为插值场ss 为克里金方差variogram_model 有 spherical、exponential、gaussian 三种常用选择城市尺度上球形模型在变程附近衰减更自然指数模型适合长距离输送明显的区域高斯模型过于平滑容易把局地高值抹掉。nlags 是经验变异函数的分组数站点数量在 2040 之间时取 68 比较稳站点太少会导致拟合出来的块金值失真。weightTrue 让拟合按每组的样本对数量加权避免远距离分组主导结果。方法主要参数优势局限反距离加权 IDW幂次 p、搜索半径无需拟合变异函数秒级出图极值被平滑给不出误差估计普通克里金变异函数模型、块金、基台、变程附带预测方差站点少于 15 个时变异函数不稳定回归加残差插值道路密度、土地利用、地形能吸收排放因子协变量获取成本高选型经验是只需要一张趋势图就用 IDWp 取 2要给出「这个位置的浓度区间」就用克里金要做污染源排查就用回归加残差因为残差场里剩下的才是局地源信号。3.3 留一交叉验证判断插值结果能不能用errors [] for i in range(len(lon)): m np.arange(len(lon)) ! i ok OrdinaryKriging(lon[m], lat[m], pm25_mean[m], variogram_modelspherical) pred, _ ok.execute(points, np.array([lon[i]]), np.array([lat[i]])) errors.append(float(pred[0]) - float(pm25_mean[i])) rmse np.sqrt(np.mean(np.square(errors))) print(LOOCV RMSE:, round(rmse, 2), 占比:, round(rmse / pm25_mean.mean(), 3))留一法每次剔除一个站点、用其余站点预测它得到的 RMSE 是该区域插值精度的直接度量。城市尺度 PM2.5 年均值下这个值通常落在 1525 μg/m³如果占比超过 30%说明站点太稀插值网格图只能定性看分布形态不能拿来读具体数值。这份验证结果也是判断要不要引入协变量的依据只要残差插值的 RMSE 比克里金低 10% 以上就值得多花时间找土地利用数据。4. PM2.5 浓度预测建模与参数调优4.1 特征工程滞后项、滚动统计与周期编码feat df.copy() for lag in [1, 2, 3, 24, 48]: feat[fpm25_lag{lag}] feat[pm25].shift(lag) # 先 shift 再 rolling否则窗口内包含当前时刻属于标签泄漏 feat[pm25_roll24] feat[pm25].shift(1).rolling(24).mean() feat[hour_sin] np.sin(2 * np.pi * feat.index.hour / 24) feat[hour_cos] np.cos(2 * np.pi * feat.index.hour / 24) feat[week_sin] np.sin(2 * np.pi * feat.index.dayofweek / 7) feat[week_cos] np.cos(2 * np.pi * feat.index.dayofweek / 7) feat feat.dropna()滞后 1、2、3 小时捕捉短时记忆滞后 24、48 小时捕捉日周期重复性这五个特征在多数城市能把线性模型误差压掉三成。滚动均值必须写成 shift(1).rolling(24)先移到上一时刻再开窗否则当前值会出现在自己的特征里线下评估漂亮上线直接失效。小时编号用 sin/cos 编码而不是 023 的整数是因为整数编码下 23 点和 0 点的距离被算成 23而实际只差一小时。周期编码把时间映射到单位圆上距离才符合物理含义。风向同样处理u、v 分量本身就是一种连续编码方案。4.2 岭回归与随机森林的基线对比from sklearn.ensemble import RandomForestRegressor from sklearn.linear_model import Ridge from sklearn.metrics import mean_absolute_error, mean_squared_error from sklearn.preprocessing import StandardScaler split int(len(feat) * 0.8) # 时间序列必须按时间切不能 shuffle train, test feat.iloc[:split], feat.iloc[split:] drop {pm25, pm25_outlier, calm, wd} X_cols [c for c in feat.columns if c not in drop] scaler StandardScaler().fit(train[X_cols]) # 归一化统计量只能来自训练集 ridge Ridge(alpha1.0).fit(scaler.transform(train[X_cols]), train[pm25]) rf RandomForestRegressor(n_estimators300, max_depth12, min_samples_leaf3, n_jobs-1, random_state42).fit(train[X_cols], train[pm25]) for name, model, X in [(Ridge, ridge, scaler.transform(test[X_cols])), (RF, rf, test[X_cols])]: pred model.predict(X) print(name, round(mean_absolute_error(test[pm25], pred), 2), round(np.sqrt(mean_squared_error(test[pm25], pred)), 2))参数常用取值对结果的影响Ridge alpha0.110越大越抑制共线特征的权重过大导致欠拟合n_estimators200500超过 300 后误差下降不明显只增加耗时max_depth815越大越容易记住噪声峰值外推变差min_samples_leaf25提高后峰值被削平日均误差反而更小评估指标的选择上MAE 反映日常预测水平RMSE 对重污染日的偏差更敏感。两份指标差距悬殊时说明模型只在低浓度段表现好。归一化统计量必须只用训练集计算用全量数据算均值和方差再切分是最隐蔽的一类数据泄漏线下 RMSE 能凭空低 20%。4.3 LSTM 预测 PM2.5滑窗、批大小与早停设置import torch import torch.nn as nn class PMNet(nn.Module): def __init__(self, n_feat, hidden64): super().__init__() self.lstm nn.LSTM(n_feat, hidden, num_layers2, batch_firstTrue, dropout0.2) self.fc nn.Linear(hidden, 1) def forward(self, x): # x: (batch, seq_len, n_feat) out, _ self.lstm(x) return self.fc(out[:, -1, :]).squeeze(-1) model PMNet(n_featlen(X_cols)) optimizer torch.optim.Adam(model.parameters(), lr1e-3) scheduler torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience3, factor0.5) criterion nn.HuberLoss(delta1.0) # 比 MSE 更抗极端值滑窗长度取 24 或 48是日周期的整数倍能让网络看到完整的一轮变化取 72 小时在宽数据上略好但显存和训练时间成倍增长。batch_size 在 64256 之间样本不足一年时取小值。num_layers 两层足够再深在这个数据量上基本只会过拟合dropout 0.2 是常规起点。损失函数用 HuberLoss 而不是 MSE是因为重污染日的浓度能到日常值的五倍以上MSE 会给这些样本极大的梯度权重模型为了压那几个点会把整体曲线拉平。delta1.0 表示误差在一倍标准差内用平方项、之外用线性项兼顾了对峰值的适度关注。早停的 patience 设 510 个 epoch监控验证集损失而不是训练集。学习率调度器在连续 3 个 epoch 不下降时把 lr 减半通常能再挤出 3%5% 的精度。训练完成后一定要做残差自相关检验如果滞后 1 小时的残差相关系数超过 0.3说明网络只是把上一时刻的值复制了一遍并没有真的学会气象驱动关系。5. 高斯烟羽模型做污染溯源参数标定与验证技巧预测只能回答「会变成什么样」溯源要回答「从哪里来」。城市尺度上常用高斯烟羽模型做快速反演公式为 C(x, y, z) Q / (2π u σy σz) × exp(−y²/(2σy²)) × [exp(−(z−H)²/(2σz²)) exp(−(zH)²/(2σz²))]。Q 是源强单位 g/su 取烟羽高度处的有效风速H 是有效源高等于烟囱几何高度加上抬升高度σy、σz 是水平与垂直扩散参数由下风向距离和大气稳定度共同决定。稳定度典型气象条件σy (m)σz (m)A 极不稳定晴天、风速 2 m/s0.32x(10.0004x)^−1/20.24x(10.001x)^1/2C 弱不稳定晴天、风速 23 m/s0.22x(10.0004x)^−1/20.20xD 中性阴天或风速 46 m/s0.16x(10.0004x)^−1/20.14x(10.0003x)^−1/2F 稳定晴夜、风速 2 m/s0.11x(10.0004x)^−1/20.08x(10.0015x)^−1/2上表是城市下垫面的 Briggs 扩散参数形式x 为下风向距离。城市地表的粗糙度让扩散比平原更快直接套用乡村参数会系统性高估下风向浓度这是溯源误差里最大的一项来源。标定 Q 和 H 时不要试图拟合绝对浓度。拿一个监测站的逐小时序列把每小时的实测风向分成 8 个或 16 个扇形统计每个扇形的浓度均值与上风向频率得到一张污染玫瑰图浓度最高且频率最低的扇区通常就是源区方位角。这一步只用到风向概率不受源强估计误差影响鲁棒性远好于最小二乘拟合。用 scipy.optimize.minimize 联合反演 Q 和 H 时目标函数建议用预测值与观测值的对数残差平方和而不是线性残差。浓度跨越一两个数量级时线性残差会被高值点主导对数残差能把量级差异摊平。反演出的 Q 若超过区域排放清单总量的量级说明假定单源不成立需要拆成多源叠加或者改用受体模型。验证环节看预测与观测的比值分布理想情况下落在 0.52 倍区间内的样本比例应超过 60%。如果这个比值随下风向距离单调增大问题一定出在 σz 上——稳定度类别选得偏低垂直扩散被低估如果比值在近距离明显偏大、远距离转小则是有效源高 H 设得过高烟羽中心没有落到采样高度。把逐小时稳定度按风速和云量重新分级后重跑一次多数情况下这一项就能把误差压掉三分之一。本文还有配套的精品资源点击获取 从AQI数据到PM2.5预测与污染溯源：空气污染建模全流程
简介这份资源是2015年第十二届五一数学建模联赛B题“空气污染问题研究”的优秀论文文档面向备战数学建模竞赛的本科生及指导教师可用于赛题复盘、论文写作参考与建模思路学习。压缩包内仅1个doc文件约1.5MB完整收录承诺书、摘要、问题重述、问题分析、模型建立与求解等论文主体结构。论文围绕京津冀地区空气污染展开依次完成国标与美标空气质量指数对比、层次分析法赋权与动态加权综合评价、修正高斯烟羽扩散模型、多污染源扩散及灰色预测模型求解并给出“APEC蓝天”治理建议报告可帮助读者掌握层次分析法、变权函数、优化高斯烟羽模型等方法的实际应用。目前已有328人学习下载适合希望借鉴优秀论文框架、学习建模论文写作规范与模型推导细节的参赛者。1. 空气污染问题研究从 AQI 数据到可复现的建模流程拿到这个标题的人手里通常已经有一份城市小时级空气质量数据却在第一步就卡住先算 AQI还是先做预测。空气污染问题研究的落点一般分三层——把六项污染物折算成可比的分指数回答「今天空气算几级」用时间序列加气象因子回答「明天 PM2.5 会不会破 150」用扩散模型反推「污染从哪个方向来」。三层要的数据粒度、模型类型和评估口径完全不同混在一起做最后容易得到一个 R² 很高但没人能解释的模型。这套流程适合做数模、做环境数据看板、做城市空气质量分析的从业者也适合第一次接触时空数据的新手。下面按数据准备、时空特征、预测建模、扩散溯源的顺序推进每一步都给能直接跑的代码和能改的参数。2. 空气污染数据的清洗与 AQI 折算2.1 六项污染物的分指数与限值对照AQI 不是对浓度直接排序得到的而是先把每项污染物单独折算成 IAQI再取最大值。折算采用分段线性插值区间断点由《环境空气质量标准》二级限值确定。污染物平均时间二级限值IAQI100 对应浓度PM2.524 小时75 μg/m³75PM1024 小时150 μg/m³150SO₂24 小时150 μg/m³150NO₂24 小时80 μg/m³80CO24 小时4 mg/m³4O₃8 小时滑动160 μg/m³160IAQI 区间级别类别050一级优51100二级良101150三级轻度污染151200四级中度污染201300五级重度污染300六级严重污染单项分指数按 IAQI (IAQI_hi − IAQI_lo)/(BP_hi − BP_lo) × (C − BP_lo) IAQI_lo 计算其中 BP_hi、BP_lo 是浓度 C 落入区间的上下断点。最终 AQI 取六项 IAQI 的最大值同时把取最大值的那项记为首要污染物。这一步在竞赛题里经常被简化成只算 PM2.5一旦题目要求评价城市整体空气状况漏掉 O₃ 就会在夏季月份给出完全相反的结论。2.2 用 pandas 完成缺失值、哨兵值与异常值处理import pandas as pd import numpy as np df pd.read_csv(air_2015_hourly.csv, parse_dates[time]) df df.set_index(time).sort_index() # asfreq 补齐缺失时刻为 NaN比 resample 更诚实不会凭空聚出平均值 df df.asfreq(h) cols [pm25, pm10, so2, no2, o3, co] # co 单位 mg/m3其余 μg/m3 # 负值和明显超量程的哨兵值-999、9999先置空否则会污染整条序列 df[cols] df[cols].mask((df[cols] 0) | (df[cols] 2000)) # 连续缺口不超过 3 小时才线性插值长缺口保留 NaN df[cols] df[cols].interpolate(methodlinear, limit3, limit_areainside) q1, q3 df[pm25].quantile([0.25, 0.75]) iqr q3 - q1 df[pm25_outlier] (df[pm25] q1 - 3 * iqr) | (df[pm25] q3 3 * iqr)逻辑上先补时间轴再补数值asfreq 把稀疏的观测时刻拉成严格小时序列后面所有 shift、rolling 才有正确的时间含义。limit3 是个经验值超过三小时的缺口用插值填平会制造出虚假的平滑趋势模型会把这部分当成真实规律学走。异常判定用 3 倍四分位距而不是常见的 1.5 倍原因是污染物浓度分布强烈右偏沙尘过程和春节烟花时段的真实峰值会被 1.5 倍规则误杀这部分样本恰恰是预测模型最该学会的部分。注意AQI 发布口径使用的是 24 小时滑动平均做短期预测时如果直接拿实时浓度训练评估结果和官方播报值对不上两者必须区分开。2.3 风向风速的向量化与静风样本取舍风向是角度量直接求平均会出大错350° 和 10° 的算术平均是 180°方向完全反了。正确做法是先转成 u、v 分量在分量上做任何线性运算需要角度时再反算回去。# wd 为风的来向单位度ws 单位 m/s rad np.deg2rad(df[wd]) df[u] -df[ws] * np.sin(rad) # 气象约定的 u 分量指向东 df[v] -df[ws] * np.cos(rad) # v 分量指向北 # 静风时风向没有物理意义置零并单独打标 df[calm] df[ws] 0.5 df.loc[df[calm], [u, v]] 0.0静风样本占比在一份年度数据里常能到 10% 以上。这些时刻污染物更接近累积状态风向信息无效直接参与训练会让模型学到噪声。把 calm 当作哑变量单独进模型或者在统计分析阶段剔除效果都比强行插值风向好。3. 空气质量的时空特征分析与插值3.1 STL 分解看趋势项与日周期项from statsmodels.tsa.seasonal import STL s df[pm25].interpolate() # 分解前必须无缺失 res STL(s, period24, seasonal7, robustTrue).fit() res.trend.to_frame(trend).join(res.seasonal.to_frame(seasonal)).to_csv(stl_out.csv)period 取 24 表示按小时数据剥离日周期如果研究的是年度规律先按天聚合再取 period7 剥离周内效应更合适。seasonal 控制季节项的平滑程度必须是奇数7 对应三个点左右25 对应十几个点取值越大抽出的日周期越干净但会吃掉部分真实波动。robustTrue 用迭代重加权最小二乘拟合单个重污染日不会把整条趋势线拉偏。分解结果读法趋势项代表排放量与区域输送的长期变化量级通常在日均值的 ±20% 内季节项反映早晚高峰的双峰结构峰值一般出现在 79 点和 1922 点残差项才是真正难以解释的部分如果残差的标准差超过日均值的 40%说明还有重要变量没进模型优先检查逆温层高度和边界层风速。3.2 IDW 与普通克里金在站点插值上的取舍from pykrige.ok import OrdinaryKriging import numpy as np ok OrdinaryKriging(lon, lat, pm25_mean, variogram_modelspherical, nlags6, weightTrue) z, ss ok.execute(grid, grid_lon, grid_lat) # z 为插值场ss 为克里金方差variogram_model 有 spherical、exponential、gaussian 三种常用选择城市尺度上球形模型在变程附近衰减更自然指数模型适合长距离输送明显的区域高斯模型过于平滑容易把局地高值抹掉。nlags 是经验变异函数的分组数站点数量在 2040 之间时取 68 比较稳站点太少会导致拟合出来的块金值失真。weightTrue 让拟合按每组的样本对数量加权避免远距离分组主导结果。方法主要参数优势局限反距离加权 IDW幂次 p、搜索半径无需拟合变异函数秒级出图极值被平滑给不出误差估计普通克里金变异函数模型、块金、基台、变程附带预测方差站点少于 15 个时变异函数不稳定回归加残差插值道路密度、土地利用、地形能吸收排放因子协变量获取成本高选型经验是只需要一张趋势图就用 IDWp 取 2要给出「这个位置的浓度区间」就用克里金要做污染源排查就用回归加残差因为残差场里剩下的才是局地源信号。3.3 留一交叉验证判断插值结果能不能用errors [] for i in range(len(lon)): m np.arange(len(lon)) ! i ok OrdinaryKriging(lon[m], lat[m], pm25_mean[m], variogram_modelspherical) pred, _ ok.execute(points, np.array([lon[i]]), np.array([lat[i]])) errors.append(float(pred[0]) - float(pm25_mean[i])) rmse np.sqrt(np.mean(np.square(errors))) print(LOOCV RMSE:, round(rmse, 2), 占比:, round(rmse / pm25_mean.mean(), 3))留一法每次剔除一个站点、用其余站点预测它得到的 RMSE 是该区域插值精度的直接度量。城市尺度 PM2.5 年均值下这个值通常落在 1525 μg/m³如果占比超过 30%说明站点太稀插值网格图只能定性看分布形态不能拿来读具体数值。这份验证结果也是判断要不要引入协变量的依据只要残差插值的 RMSE 比克里金低 10% 以上就值得多花时间找土地利用数据。4. PM2.5 浓度预测建模与参数调优4.1 特征工程滞后项、滚动统计与周期编码feat df.copy() for lag in [1, 2, 3, 24, 48]: feat[fpm25_lag{lag}] feat[pm25].shift(lag) # 先 shift 再 rolling否则窗口内包含当前时刻属于标签泄漏 feat[pm25_roll24] feat[pm25].shift(1).rolling(24).mean() feat[hour_sin] np.sin(2 * np.pi * feat.index.hour / 24) feat[hour_cos] np.cos(2 * np.pi * feat.index.hour / 24) feat[week_sin] np.sin(2 * np.pi * feat.index.dayofweek / 7) feat[week_cos] np.cos(2 * np.pi * feat.index.dayofweek / 7) feat feat.dropna()滞后 1、2、3 小时捕捉短时记忆滞后 24、48 小时捕捉日周期重复性这五个特征在多数城市能把线性模型误差压掉三成。滚动均值必须写成 shift(1).rolling(24)先移到上一时刻再开窗否则当前值会出现在自己的特征里线下评估漂亮上线直接失效。小时编号用 sin/cos 编码而不是 023 的整数是因为整数编码下 23 点和 0 点的距离被算成 23而实际只差一小时。周期编码把时间映射到单位圆上距离才符合物理含义。风向同样处理u、v 分量本身就是一种连续编码方案。4.2 岭回归与随机森林的基线对比from sklearn.ensemble import RandomForestRegressor from sklearn.linear_model import Ridge from sklearn.metrics import mean_absolute_error, mean_squared_error from sklearn.preprocessing import StandardScaler split int(len(feat) * 0.8) # 时间序列必须按时间切不能 shuffle train, test feat.iloc[:split], feat.iloc[split:] drop {pm25, pm25_outlier, calm, wd} X_cols [c for c in feat.columns if c not in drop] scaler StandardScaler().fit(train[X_cols]) # 归一化统计量只能来自训练集 ridge Ridge(alpha1.0).fit(scaler.transform(train[X_cols]), train[pm25]) rf RandomForestRegressor(n_estimators300, max_depth12, min_samples_leaf3, n_jobs-1, random_state42).fit(train[X_cols], train[pm25]) for name, model, X in [(Ridge, ridge, scaler.transform(test[X_cols])), (RF, rf, test[X_cols])]: pred model.predict(X) print(name, round(mean_absolute_error(test[pm25], pred), 2), round(np.sqrt(mean_squared_error(test[pm25], pred)), 2))参数常用取值对结果的影响Ridge alpha0.110越大越抑制共线特征的权重过大导致欠拟合n_estimators200500超过 300 后误差下降不明显只增加耗时max_depth815越大越容易记住噪声峰值外推变差min_samples_leaf25提高后峰值被削平日均误差反而更小评估指标的选择上MAE 反映日常预测水平RMSE 对重污染日的偏差更敏感。两份指标差距悬殊时说明模型只在低浓度段表现好。归一化统计量必须只用训练集计算用全量数据算均值和方差再切分是最隐蔽的一类数据泄漏线下 RMSE 能凭空低 20%。4.3 LSTM 预测 PM2.5滑窗、批大小与早停设置import torch import torch.nn as nn class PMNet(nn.Module): def __init__(self, n_feat, hidden64): super().__init__() self.lstm nn.LSTM(n_feat, hidden, num_layers2, batch_firstTrue, dropout0.2) self.fc nn.Linear(hidden, 1) def forward(self, x): # x: (batch, seq_len, n_feat) out, _ self.lstm(x) return self.fc(out[:, -1, :]).squeeze(-1) model PMNet(n_featlen(X_cols)) optimizer torch.optim.Adam(model.parameters(), lr1e-3) scheduler torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, patience3, factor0.5) criterion nn.HuberLoss(delta1.0) # 比 MSE 更抗极端值滑窗长度取 24 或 48是日周期的整数倍能让网络看到完整的一轮变化取 72 小时在宽数据上略好但显存和训练时间成倍增长。batch_size 在 64256 之间样本不足一年时取小值。num_layers 两层足够再深在这个数据量上基本只会过拟合dropout 0.2 是常规起点。损失函数用 HuberLoss 而不是 MSE是因为重污染日的浓度能到日常值的五倍以上MSE 会给这些样本极大的梯度权重模型为了压那几个点会把整体曲线拉平。delta1.0 表示误差在一倍标准差内用平方项、之外用线性项兼顾了对峰值的适度关注。早停的 patience 设 510 个 epoch监控验证集损失而不是训练集。学习率调度器在连续 3 个 epoch 不下降时把 lr 减半通常能再挤出 3%5% 的精度。训练完成后一定要做残差自相关检验如果滞后 1 小时的残差相关系数超过 0.3说明网络只是把上一时刻的值复制了一遍并没有真的学会气象驱动关系。5. 高斯烟羽模型做污染溯源参数标定与验证技巧预测只能回答「会变成什么样」溯源要回答「从哪里来」。城市尺度上常用高斯烟羽模型做快速反演公式为 C(x, y, z) Q / (2π u σy σz) × exp(−y²/(2σy²)) × [exp(−(z−H)²/(2σz²)) exp(−(zH)²/(2σz²))]。Q 是源强单位 g/su 取烟羽高度处的有效风速H 是有效源高等于烟囱几何高度加上抬升高度σy、σz 是水平与垂直扩散参数由下风向距离和大气稳定度共同决定。稳定度典型气象条件σy (m)σz (m)A 极不稳定晴天、风速 2 m/s0.32x(10.0004x)^−1/20.24x(10.001x)^1/2C 弱不稳定晴天、风速 23 m/s0.22x(10.0004x)^−1/20.20xD 中性阴天或风速 46 m/s0.16x(10.0004x)^−1/20.14x(10.0003x)^−1/2F 稳定晴夜、风速 2 m/s0.11x(10.0004x)^−1/20.08x(10.0015x)^−1/2上表是城市下垫面的 Briggs 扩散参数形式x 为下风向距离。城市地表的粗糙度让扩散比平原更快直接套用乡村参数会系统性高估下风向浓度这是溯源误差里最大的一项来源。标定 Q 和 H 时不要试图拟合绝对浓度。拿一个监测站的逐小时序列把每小时的实测风向分成 8 个或 16 个扇形统计每个扇形的浓度均值与上风向频率得到一张污染玫瑰图浓度最高且频率最低的扇区通常就是源区方位角。这一步只用到风向概率不受源强估计误差影响鲁棒性远好于最小二乘拟合。用 scipy.optimize.minimize 联合反演 Q 和 H 时目标函数建议用预测值与观测值的对数残差平方和而不是线性残差。浓度跨越一两个数量级时线性残差会被高值点主导对数残差能把量级差异摊平。反演出的 Q 若超过区域排放清单总量的量级说明假定单源不成立需要拆成多源叠加或者改用受体模型。验证环节看预测与观测的比值分布理想情况下落在 0.52 倍区间内的样本比例应超过 60%。如果这个比值随下风向距离单调增大问题一定出在 σz 上——稳定度类别选得偏低垂直扩散被低估如果比值在近距离明显偏大、远距离转小则是有效源高 H 设得过高烟羽中心没有落到采样高度。把逐小时稳定度按风速和云量重新分级后重跑一次多数情况下这一项就能把误差压掉三分之一。本文还有配套的精品资源点击获取