简介这是一篇数学建模省级优秀奖论文《城市空气质量评估及预测》面向数学建模竞赛参赛者、环境数据分析学习者。资源以层次分析法完成10个城市空气污染严重程度排名并利用指数平滑法预测成都11月空气质量同时分析首要污染物及影响空气污染的主要因素可作为建模方法落地、论文结构安排与竞赛获奖规范的完整参考。包体为1个doc文档约846KB正文含摘要、问题重述、基本假设、模型建立与求解、结果分析与建议等完整章节且给出Matlab求特征值、一致性检验及Excel统计污染天数的具体思路便于读者对照复现。目前已有885人学习下载适合正在备赛数学建模、需要同类城市环境问题案例或想学习层次分析法与指数平滑法应用的读者。1. 城市空气质量评估与预测为什么数据清洗比模型选型更决定奖项归属城市空气质量评估与预测这道数模题数据是公开的、模型是现成的为什么还有一大堆队伍连省级奖项都摸不到我做过参赛队员也帮老师看过卷最深的体会是这题比的不是谁的模型更花哨而是谁先把数据收拾明白。拿到“城市空气质量评估及预测”这个题第一步不是去翻论文选模型而是先搞清楚监测站点的数据长什么样哪个字段是浓度、哪个字段是时间、缺了多少值、两套数据的时间对不对得上。数据管线搭完评估和预测基本就是套公式的事。这篇笔记按这个顺序把从数据预处理、评估建模、预测建模到论文呈现的完整链路过一遍顺带把那些让人翻车的细节都讲透。2. 数据清洗才是主力战场站点数据、气象数据怎么对齐才不会跑偏一个典型数模题从发布到交卷只有三天左右大多数队伍第一天全耗在数据上了。城市空气质量评估与预测这个题数据通常来自两类渠道环保部门公布的污染物浓度监测数据以及气象站发布的逐日或逐时气象数据。前者至少包含 PM2.5、PM10、SO₂、NO₂、CO、O₃ 六项常规污染物后者包含温度、湿度、气压、风速、风向。两份数据的发布时间粒度不一致是后面所有预测误差的第一个来源。2.1 数据源和字段表先搞清楚手里有什么拿到数据后我做的第一件事不是跑代码而是建一张字段核对表。监测数据常见字段是站点编号、时间、六项污染物浓度和 AQI气象数据常见字段是时间、温度、相对湿度、气压、风向、风速。很多队伍上来就用 Excel 直接打开看一眼就开始建模结果做预测时才发现气象时间比监测时间晚了 8 个小时整个模型的预测曲线整体偏移评分直接崩掉。数据类别常见字段时间粒度典型来源污染物监测station, time, pm25, pm10, so2, no2, co, o3, aqi逐小时或逐日环境监测站点日报、开放数据平台气象观测time, temp, rh, pressure, wind_dir, wind_speed逐小时或逐 3 小时国家气象信息中心、气象数据共享平台辅助数据city, lon, lat, population静态统计年鉴、城市公开数据我一般会把所有文件读进来后先打印dtypes和head()核对时间字段是不是解析成了 datetime。这一步看起来基础但能省掉后面大量排查时间。文件名里的 .doc 如果是别人的成稿里面通常还带了站点经纬度、人口密度这些静态数据这些都可以作为评估模型的补充特征但前提是先把数据读完、字段对齐。2.2 缺失值处理别信均值填补先看缺失模式缺失值处理是整道题最容易写得像“玄学”的地方。常见做法是df.fillna(df.mean())一把梭但这么做有两个问题一是当缺失集中在某几天时均值填补会把那几天的波动拉平预测模型学不到真实特征二是如果某些站点长期缺测均值填补会把站点差异抹掉评估排名就失真了。我一般会先按站点分组计算缺失比例缺失少的用线性插值缺失多的站点直接剔除或标注为不可用。import pandas as pd import numpy as np df pd.read_csv(air_quality.csv, parse_dates[time]) df df.sort_values([station, time]).reset_index(dropTrue) # 按站点检查各指标的缺失比例 miss_rate df.groupby(station)[pm25].apply( lambda s: s.isna().sum() / len(s) ) print(miss_rate[miss_rate 0.1])这段代码先按站点和时间排序再计算每个站点 PM2.5 的缺失比例。isna().sum() / len(s)算的是缺失占比阈值我一般取 0.1。缺失超过 10% 的站点后续做时序预测时风险很大宁可去掉也不硬补。# 连续缺失不多时用线性插值limit 控制最大连续插值数 df[pm25] df[pm25].interpolate(methodlinear, limit12) df[pm10] df[pm10].interpolate(methodlinear, limit12) # 对极端异常值做窗口过滤超过前后 5 天中位数 5 倍的点视为异常 window_med df.groupby(station)[pm25].transform( lambda s: s.rolling(5, centerTrue).median() ) df.loc[df[pm25] window_med * 5, pm25] np.nan df[pm25] df[pm25].interpolate(methodlinear, limit12)interpolate(methodlinear)是线性插值适合短时间连续缺失limit12表示最多连续补 12 个点超过就不再补避免长段缺失被虚假填充。异常值过滤这步容易被忽略监测仪器偶发故障会产生一个“浓度 2000”之类的离谱点如果不处理后续模型计算会直接带偏。用滚动中位数做阈值比用均值更稳因为中位数不受极端值本身影响。2.3 时间对齐UTC 和北京时间的 8 小时坑这个坑我踩过不止一次。气象数据很多用 UTC 时间记录监测数据基本是北京时间两个时间差 8 小时直接合并会让当天的气象因子对到第二天的污染物浓度上。做预测时这个错位会藏在模型系数里看不出来但一到解释结论就被评委问倒。# 气象数据统一到北京时间 met[time_bj] ( pd.to_datetime(met[time], utcTrue) .dt.tz_convert(Asia/Shanghai) .dt.tz_localize(None) ) # 监测数据本身是北京时间直接解析 aqi[time_bj] pd.to_datetime(aqi[time]) # 以站点时间为主键做左连接保留监测数据的全部时间点 merged pd.merge( aqi, met, on[station, time_bj], howleft ) print(merged.isna().sum())这里用utcTrue先把字符串解析成带时区的时间戳再用tz_convert转到北京时间最后用tz_localize(None)把时区信息去掉变成普通的 pandas 时间对象。这样merge才按主键对齐。合并后打印缺失值能发现大多数缺失集中在第一天和最后一天这是因为两份数据的起止日期不完全一致属于正常现象不需要补。时间对齐做完后面所有的评估和预测才有基础。这个环节没有可跳跃的捷径多花半小时做对齐能省掉后面两小时的排查。3. 空气质量评估从单因子指数到主成分综合得分这一步做扎实排名就立住了评估部分的目标是回答“哪个城市空气质量更好”或者“某城市空气质量在时间上怎么变化”。很多队伍直接把官方 AQI 拿过来一排名就交差但评阅老师想看到的是你理解 AQI 的计算逻辑并且能提出自己的综合评估方法。3.1 官方 AQI 计算逻辑先走一遍中国空气质量指数 AQI 依据 HJ 633-2012 计算核心是把六项污染物浓度各自折算成 0 到 500 的单项评价指数再取最大值。单项指数 IAQI 的计算公式是分段线性插值IAQI (IAQI_hi - IAQI_lo) / (C_hi - C_lo) * (C - C_lo) IAQI_lo其中 C 是实测浓度C_lo 和 C_hi 是该项污染物在标准表中对应的浓度区间边界IAQI_lo 和 IAQI_hi 是对应的指数区间边界。六项污染物中CO 的单位是 mg/m³其他五项是 μg/m³单位混用是这里最常见的翻车点。污染物浓度限值24小时均值对应 IAQI100 边界PM2.575 μg/m³75PM10150 μg/m³150SO₂150 μg/m³150NO₂80 μg/m³80CO4 mg/m³4O₃160 μg/m³8小时滑动160我建议即使题目没要求算 AQI也先把这套官方计算跑通拿它作为后续自建模型的基准线。AOI 的最大值取法有一个明确缺陷它只看污染最严重的单项其他污染项的信息全被丢掉了。这就是为什么数学建模题里还要做综合评估。3.2 主成分分析作为综合评价工具的适用边界主成分分析 PCA 的作用是把六项高度相关的污染物指标压缩成两到三个互不相关的综合指标再按方差贡献率加权得到总得分。选它的理由是六项污染物浓度彼此强相关比如 PM2.5 和 PM10 的相关系数经常超过 0.8直接用原始指标做加权平均会产生信息冗余而 PCA 可以自动去相关权重由数据本身决定不需要主观设定。但 PCA 有适用前提变量之间要存在一定的线性相关性。如果六项指标相关性很弱PCA 压缩后第一主成分方差贡献率可能不到 40%那就不适合用这时可以换熵权法。熵权法通过计算各指标的信息熵来确定权重指标差异越大权重越高也是完全客观的适合在论文里和 PCA 做对比。3.3 sklearn 实现 PCA标准化、降维、计算综合得分的完整流程from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA features [pm25, pm10, so2, no2, co, o3] X merged[features].dropna().values # 标准化让六项指标处在同一量纲 scaler StandardScaler() X_std scaler.fit_transform(X) # 先跑全量 PCA 看方差贡献率 pca_full PCA(n_components6, random_state42) pca_full.fit(X_std) print(pca_full.explained_variance_ratio_)先输出六个主成分的方差贡献率观察前几个主成分累积是否达到 80% 以上。如果前两个主成分就累积超过 80%说明城市间的空气质量差异主要由两类因素主导这时再降维才有意义。# 正式降维到 2 个主成分计算综合得分 pca PCA(n_components2, random_state42) scores pca.fit_transform(X_std) # 按方差贡献率加权得到综合得分 w pca.explained_variance_ratio_ composite scores[:, 0] * w[0] scores[:, 1] * w[1] rank np.argsort(composite)[::-1] # 得分从高到低排序random_state42是固定随机种子保证每次运行结果一致这在高分论文里是必要条件。explained_variance_ratio_返回每个主成分的解释方差占比作为权重时表示该主成分对原始信息的保留程度。最后算出的 composite 值越低表示污染越轻、空气质量越好排序时注意方向别把排名写反。做完 PCA 之后必须看载荷矩阵pca.components_它会告诉你每个主成分和哪些原始指标强相关。一个容易解释的模式是第一主成分在 PM2.5、PM10、CO 上载荷高代表燃煤和机动车排放的一次污染第二主成分在 O₃ 和 NO₂ 上载荷高代表光化学二次污染。这样写进论文评阅老师会觉得你们不是拿 sklearn 跑一遍交差而是真正理解了指标含义。3.4 熵权法做对照让评估结论更有说服力PCA 的权重是“方差最大”导向熵权法的权重是“信息量”导向。做法是先把各指标做归一化然后计算每个指标的熵值 E权重 W (1 - E) / Σ(1 - E)最后用原始浓度乘权重求和得到综合得分。它的优势是完全不依赖线性相关性假设缺点是没有考虑指标之间的重叠信息。论文里用两种方法各算一遍城市排名然后用 Spearman 秩相关系数检验排名一致性这是省级优秀奖论文里的常见加分操作。4. 预测模型ARIMA 做底线、LSTM 做上限、GM(1,1) 救小样本预测部分的常见问题是“只做一种模型不做对比”。评阅老师看过的卷子中至少一半队伍只跑了一个 ARIMA 就写完了。获奖论文的普遍做法是用两种以上模型做预测并给出误差对比让评委看到你们了解每种模型的适用边界。4.1 模型选型预判表先想清楚再做模型适用场景最短样本量输出主要坑ARIMA线性、平稳或差分后平稳的序列至少 30 个时间点点预测区间定阶靠猜过差分LSTM非线性、长周期依赖至少几百个样本点预测数据泄漏、训练不稳定GM(1,1)样本极少、指数趋势8-20 个点点预测只适合短期、光滑序列我做预测的顺序永远是先跑 ARIMA 作为 baseline再跑 LSTM 看能不能更好最后用小样本场景补一个灰色预测。这样模型的三个层次都覆盖了论文里“模型对比”这一节就不会空。4.2 ARIMA从差分到定阶的完整代码ARIMA 的数学逻辑是对非平稳序列做差分去掉趋势得到一个平稳序列AR 项捕捉序列对自身历史值的依赖MA 项捕捉历史白噪声的影响。参数 (p, d, q) 里的 d 是差分阶数p 是自回归阶数q 是移动平均阶数。from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.arima.model import ARIMA series merged.set_index(time_bj)[pm25].dropna() series series.resample(D).mean() # 统一为日均值 # 单位根检验p 值大于 0.05 说明序列非平稳 adf_result adfuller(series.dropna()) print(ADF p-value:, adf_result[1]) # 差分一次后再次检验 diff1 series.diff().dropna() adf_result_diff adfuller(diff1) print(Diff1 ADF p-value:, adf_result_diff[1])如果原始序列的 ADF 检验 p 值大于 0.05说明序列存在单位根、非平稳需要差分。d 的取值从 0 开始逐个试每次差分后重新做检验直到 p 值小于 0.05 就停止。不要看到 p 值不小就一直差分差分次数过多会把序列里的有效信息差没。# 通过 ACF 和 PACF 图观察截尾/拖尾特征辅助定阶 from statsmodels.graphics.tsaplots import plot_acf, plot_pacf plot_acf(diff1, lags30) plot_pacf(diff1, lags30)ACF 和 PACF 的判断方法比较依赖经验PACF 在某个阶数后突然截尾对应 AR 的 pACF 截尾对应 MA 的 q。如果看不出明显截尾直接用 AIC 准则在 (0,1,0) 到 (3,1,3) 之间网格搜索选 AIC 最小的组合这是最省事也最不容易被质疑的做法。import itertools best_aic, best_order float(inf), None for p in range(0, 4): for q in range(0, 4): try: model ARIMA(series, order(p, 1, q)) fit model.fit() if fit.aic best_aic: best_aic, best_order fit.aic, (p, 1, q) except Exception: continue print(best order:, best_order, AIC:, best_aic) final_model ARIMA(series, orderbest_order).fit() forecast final_model.forecast(steps7) print(forecast)AIC是赤池信息准则值越小表示模型在拟合优度和复杂度之间平衡得越好。网格搜索时 p 和 q 一般不超过 3因为城市空气质量数据里超过三阶的自回归关系已经非常罕见。forecast(steps7)表示预测未来 7 天的 PM2.5 日均值。ARIMA 的输出是数值列表论文里要配上置信区间final_model.get_forecast(steps7).conf_int()可以拿到上下界。4.3 LSTM建窗和归一化是决定成败的两件事LSTM 做时间序列预测的直觉逻辑是用过去 N 天的数据作为窗口预测第 N1 天的值。它比 ARIMA 强的地方在于能捕捉非线性关系但代价是需要大量数据而且对数据预处理方式极其敏感。import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense from sklearn.preprocessing import MinMaxScaler data series.dropna().values.reshape(-1, 1) scaler MinMaxScaler(feature_range(0, 1)) data_scaled scaler.fit_transform(data) # 构建滑动窗口用前 7 天预测第 8 天 def make_window(data, n7): X, y [], [] for i in range(len(data) - n): X.append(data[i:in, 0]) y.append(data[in, 0]) return np.array(X), np.array(y) X, y make_window(data_scaled, n7) X X.reshape((X.shape[0], X.shape[1], 1)) # 前 80% 训练后 20% 测试测试集不参与训练 split int(len(X) * 0.8) X_train, X_test X[:split], X[split:] y_train, y_test y[:split], y[split:] model Sequential([ LSTM(32, activationrelu, input_shape(7, 1)), Dense(1) ]) model.compile(optimizeradam, lossmse) history model.fit(X_train, y_train, epochs50, batch_size16, verbose0) # 反归一化回真实浓度单位 y_pred model.predict(X_test) y_pred scaler.inverse_transform(y_pred) y_test_orig scaler.inverse_transform(y_test.reshape(-1, 1))这段代码有几个关键参数n7表示窗口长度也就是用过去一周预测下一天。窗口不是越大越好过大的窗口不仅增加计算量还会把久远的历史噪声带进来。LSTM(32)是隐藏单元数城市空气质量这种单变量序列 16 到 64 都够用。epochs50是训练轮数建议保存训练过程的 loss 曲线看到 loss 不再明显下降就停止。batch_size16是每批样本数样本总量小的时候 batch 设小一点有利于收敛。这里最隐蔽的坑在make_window函数里如果构建 X 时不小心把 y 也包含进了窗口测试集精度会高得不正常肉眼看着完美但一到真实预测就没法用这就是典型的数据泄漏。检查方法很简单打印 X[i] 的最后一个值和 y[i] 的值正常情况下 X[i] 最后一位是第 i6 天的数据y[i] 是第 i7 天的数据两者不是同一个时间点。4.4 GM(1,1)小样本场景下的灰色预测有些城市站点只有十几天的有效数据ARIMA 定阶都不够LSTM 更是没法训练这时灰色预测 GM(1,1) 是唯一还能给出预测结果的方案。它不对原始序列建模而是对原始序列的一次累加生成序列建模。核心公式是x_next (x0 - b/a) * (1 - exp(-a)) * exp(-a * k)其中 a 是发展系数b 是灰作用量通过最小二乘法从累加序列中估计出来。GM(1,1) 只适合短期预测一般用来预测 3 到 5 个点超过这个范围预测值会迅速退化。这个模型的价值不在精度而在于证明你“把所有能试的方法都试过了”在评阅逻辑里属于完整性加分项。4.5 预测误差评估只看 RMSE 不够方向准确率被多数人忽略from sklearn.metrics import mean_absolute_error, mean_squared_error rmse np.sqrt(mean_squared_error(y_test_orig, y_pred)) mae mean_absolute_error(y_test_orig, y_pred) # 方向准确率真实趋势和预测趋势同号的占比 dir_real np.sign(np.diff(y_test_orig[:, 0])) dir_pred np.sign(np.diff(y_pred[:, 0])) dir_acc np.mean(dir_real dir_pred) print(RMSE:, rmse, MAE:, mae, Direction Acc:, dir_acc)RMSE 对误差大的点敏感MAE 对异常点更稳健两者一起看能避免单一指标被个别极端值主导。方向准确率衡量的是“明天涨还是跌”判断得对不对这个指标在很多预测类赛题里比 RMSE 更贴近实际关心的问题。建议论文里把三种模型的这三个指标做成一张对比表这比贴一堆曲线图更能让评委快速抓住结论。5. 数模评审最看重的细节5 条把省奖拉回参赛奖的常见翻车现场三天比赛时间赶出来的论文翻车点高度集中在数据逻辑和模型细节上而不是模型的数学推导本身。评委的耐心有限下面这几个问题只要出现在论文里奖项基本就和你们告别了。5.1 AQI 计算把 24 小时均值浓度当成日均值来折算现象算出来的 IAQI 比官方公布值低很多或高很多站点排名和公开数据对不上。 原因六项污染物中有几项在标准里用的是 24 小时平均值比如 PM2.5 日均值限值是 75 μg/m³但如果手里是逐小时数据直接拿某个时刻的小时浓度去套日均值限值折算指数会忽高忽低。还有一种是把 O₃ 的 8 小时滑动平均浓度直接当成小时浓度去算。 解决先明确手头数据的统计口径逐小时数据先按天聚合只有日均序列才能套用日均限值表。写 AQI 计算函数时把单位也一并处理CO 是 mg/m³其他是 μg/m³代码里统一注意单位。5.2 主成分分析没有标准化第一主成分被 PM10 绑架现象PCA 跑完后第一主成分的方差贡献率极高但载荷几乎全部集中在 PM10 上其他五个指标的载荷接近零。 原因PM2.5、PM10 的数值范围是几十到几百而 CO 的数值范围是 0.5 到 4数值范围大的变量在协方差矩阵里天然占据主导地位。没有做标准化就相当于给 PM10 加了超高权重PCA 的意义就丢了。 解决在PCA.fit之前必须用StandardScaler对数据做标准化。判断标准检查pca.components_矩阵看每个主成分的载荷是否覆盖了多个原始指标如果某个指标独占一个主成分说明预处理有问题的可能性大。5.3 ARIMA 反复差分直到信号被差没了现象ARIMA 模型 AIC 越来越低但预测结果几乎是一条水平直线或者预测值剧烈震荡。 原因对已经平稳的序列继续差分会把序列中的真实波动逐步抹掉模型最后学到的只剩随机噪声。常见做法是看到 ADF 检验 p 值大于 0.05 就一直差分差分到第三阶还在继续试序列早就退化成了纯噪声。 解决d 的取值从 0 开始每加一阶都重新做 ADF 检验p 值小于 0.05 就停。差分后把序列画出来看目测是否平稳两个条件同时满足才进入下一步。d 最大一般取到 2超过 2 就要怀疑数据本身有问题。5.4 LSTM 里悄悄混入了未来数据现象测试集 R² 高过 0.95预测曲线基本贴合真实值看似完美但换成新数据就完全失效。 原因构建窗口时把目标时间点的数据也放进了输入窗口。比如用前 7 天预测第 8 天但窗口取的是第 1 到第 8 天第 8 天的真实值已经作为输入喂给了模型。另一种常见做法是用当天的气象数据预测当天的 PM2.5在建模阶段能成立但真实场景里当天气象数据是未知的这就属于顺向数据泄漏。 解决构建窗口时严格区分特征和目标make_window函数里把X切到in之前y从in开始取值。测试时拿到真实值再反归一化做误差评估但 LSTM 模型本身只用历史数据预测未来。5.5 时间戳合并错了预测曲线整体滞后 8 小时或 1 天现象模型效果显示趋势对得上但每天的预测峰谷都比真实值晚 8 小时RMSE 数值不大图一拉出来就穿帮。 原因气象数据用 UTC 记录监测数据是北京时间或者两份数据虽然都是北京时间但一个记录的是“日期开头零点”另一个记录的是“日期结尾 24 点”差了一整天。 解决在 merge 之前先打印两份数据的时间范围逐条核对某一天的具体时间值。统一用pd.to_datetime(..., utcTrue).dt.tz_convert(Asia/Shanghai)转时区合并后做一次交叉验证画出某站点 PM2.5 和温度两条曲线如果污染物浓度峰值比温度峰值明显滞后半天到一天这是符合大气化学常识的如果完全同步说明时间对齐大概率已经错了。6. 用灵敏度分析把论文厚度撑起来三个必做的验证实验到了这个阶段模型已经跑通、论文初稿也成形了但还差一道工序灵敏度分析。评阅老师看到“模型对不同参数设置是否稳健”这一节会明显更认可工作量的完整性。这里我一般做三个实验全部不依赖额外数据源只花一个晚上就能完成。第一个实验是参数扰动分析。对 PCA 里取主成分个数分别设为 2、3、4重算城市排名计算各组排名之间的 Spearman 秩相关系数。如果相关系数都在 0.9 以上说明评估结果对参数不敏感可以写“综合评估排名受主成分个数影响极小”。同理对 ARIMA 的 p、q 在最优值附近 ±1 扰动看 7 天预测的 RMSE 变化幅度变化不超过 15% 就认为模型稳定。第二个实验是训练/测试切分比例验证。将 LSTM 的训练集比例从 70% 改成 75%、80%、85% 各跑一次记录 RMSE。如果性能波动很大说明样本量不足以支撑 LSTM这时论文里要明确写“LSTM 在样本量较小时表现不稳定实际场景建议优先选用 ARIMA”这个诚实的结论比硬撑 LSTM 更得评委好感。第三个实验是真预测的可行性检验。把最后一段时间完整藏起来比如最后 10 天不参与训练用前 60 天训练、预测后 10 天再把预测结果和真实值对齐画图。这模拟的是真实应用场景历史数据只到昨天今天及未来的值就是靠模型推出来的。这个图表放论文里标题写“基于历史窗口的向前滚动验证”视觉说服力远强于普通的训练集测试集对比图。数据可视化的规范也要一并检查坐标轴必须带物理单位和量级μg/m³ 或 mg/m³、折线图必须标明用哪个模型、预测区间用浅色填充带表示、对比表格里 RMSE 和 MAE 保留两位小数而不是一长串。评委翻论文的时候先看图后看公式图脏乱差公式再漂亮也没人愿意细看。我自己做这个题最大的教训是别把时间花在调一个更复杂的模型上把 PCA 载荷解释清楚、把时间对齐的口径写明白、把灵敏度分析补全这三件事对奖项的贡献比任何模型技巧都大。数据管线的每一步都留下了脚印评阅老师能顺着你的脚印走完整条逻辑链他就愿意相信你的结果是可靠的。希望帮到你。本文还有配套的精品资源点击获取
