ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

城市空气质量评估与预测:数学建模中的数据分析与机器学习方法

城市空气质量评估与预测:数学建模中的数据分析与机器学习方法 简介城市空气质量评估及预测是一份获省级优秀奖的数学建模竞赛论文面向数学建模参赛者、环境数据分析及城市规划相关人群重点解决多城市空气质量对比排序和未来状况预测问题。资源共1个文件为完整doc文档大小846KB包含摘要、问题提出、基本假设、模型建立与求解全过程其中层次分析法部分详细展示了判断矩阵构造、Matlab求最大特征值与一致性检验、组合权重计算指数平滑法部分则给出了成都11月空气质量预测模型及平滑系数选择思路。已有885人学习下载。读者通过该文档既能学习完整建模流程也能掌握用Excel统计污染天数、用Matlab求解权重等实操技能还可借鉴其围绕首要污染物探究成因并提出环保建议的展开方式对备战数学建模竞赛或开展环境数据研究均有直接参考价值。1. 数学建模里的空气质量评估与预测先搞清楚这道题在考什么城市空气质量评估及预测是数学建模竞赛里出现频率极高的一类题目。表面看是环境问题实际考的是三件事能不能把监测数据清洗成可用样本能不能用合理的指标体系给城市空气质量排序或定级以及能不能把时间序列预测做到误差可接受。当年我拿省级优秀奖的这篇论文核心思路其实很朴素污染浓度是连续监测数据空气质量等级是由六项污染物浓度换算出来的离散标签两者之间的映射关系适合用统计模型描述而浓度随时间的变化则适合用时间序列和机器学习模型捕捉。适合谁一句话说清手里有一份城市逐日空气质量数据、想参加数模竞赛或者做环境数据分析的人这篇文章的评估思路和预测流程可以直接复现。2. 数据预处理与指标构建评估和预测之前先把数据底子打牢很多人拿到题目直接套模型结果评估排序和预测全都不理想。问题往往不在模型本身而是数据底子没打牢。城市空气质量数据通常来自国控监测站点包含 SO₂、NO₂、CO、O₃、PM10、PM2.5 六项污染物的小时浓度或日均浓度外加气象数据如温度、湿度、风速、风向、气压。省级优秀奖和国赛省一之间差的往往不是模型复杂度而是对 AQI 计算口径和数据质量的把控。2.1 AQI 计算口径六项污染物与 IAQI 分段线性插值评估城市空气质量第一步不是把所有污染物浓度加起来算平均而是先算单项空气质量指数 IAQI。国家标准把每项污染物的浓度区间和对应的 IAQI 分段绑在一起区间内按线性插值计算。不同污染物的分段阈值不同很多人栽在 O₃ 的统计口径上O₃ 用的是日最大 8 小时滑动平均浓度不是全天平均也不是小时最大值。# 计算单项IAQI以PM2.5为例 def calc_iaqi_pm25(c_pm25): # (浓度上限, 浓度下限, IAQI上限, IAQI下限) breakpoints [ (35, 0, 50, 0), (75, 35, 100, 50), (115, 75, 150, 100), (150, 115, 200, 150), (250, 150, 300, 200), (350, 250, 400, 300), (500, 350, 500, 400), ] for bp_high, bp_low, iaqi_high, iaqi_low in breakpoints: if c_pm25 bp_high: return (iaqi_high - iaqi_low) / (bp_high - bp_low) * (c_pm25 - bp_low) iaqi_low return 500 # 示例PM2.5浓度为82微克/立方米落在75~115区间 print(calc_iaqi_pm25(82))这段代码的逻辑是逐段判断浓度落在哪个区间然后按比例插值。注意 PM2.5 的浓度超过 500 时直接取 500这是国标明确规定的封顶值。实际项目中我会把六项污染物都写成类似函数然后取最大值作为当日 AQI。这里有个细节AQI 不是六项 IAQI 的平均值而是最大值意味着当天空气质量等级由最差的那项污染物决定。做评估时把 AQI 分级映射成 1 到 6 的等级标签常见的分段是优、良、轻度污染、中度污染、重度污染、严重污染。2.2 数据清洗的隐藏工作量缺失值、异常值和日期对齐监测数据里最常见的三类问题站点停机导致的小时数据缺失、仪器校准产生的负值或极端离群值、以及不同数据源日期格式不一致。省级优秀奖的评审不会看你的数据清洗代码但数据清洗直接影响模型分数。缺失值处理我一般按缺失比例分层处理单日缺失不超过 20% 的用前后两天均值插补连续多日缺失的站点直接剔除该时段不硬插负浓度值直接置为空再插补。import pandas as pd import numpy as np # 假设df是逐日数据列包含[date, PM2.5, PM10, O3, NO2, SO2, CO] df[date] pd.to_datetime(df[date]) df df.sort_values(date).reset_index(dropTrue) # 负值置空 for col in [PM2.5, PM10, O3, NO2, SO2, CO]: df.loc[df[col] 0, col] np.nan # 缺失率高于30%的站点列直接丢弃其余用前后均值插补 missing_ratio df.isnull().mean() valid_cols missing_ratio[missing_ratio 0.3].index.tolist() df df[valid_cols].interpolate(methodlinear, limit_directionboth)逻辑说明先把负值当作缺失处理避免把仪器故障数据当成真实浓度然后按列计算缺失率超过 30% 的列说明该站点或污染物数据质量太差插补反而会引入系统性偏差最后用线性插补填充小比例缺失。参数说明limit_directionboth 表示序列开头和结尾的缺失值也会被填充如果不加这个参数开头缺失会保留 NaN后续建模直接报错或丢弃整行。2.3 特征工程的边界气象要素与滞后特征怎么选评估用哪些指标、预测用哪些特征两者要分开设计。评估层面除了六项污染物浓度还会加入气象因素的负面影响比如静稳天气指数预测层面特征要包含历史浓度、气象预报值和日历特征。一个常见做法是构造滞后特征用前一天的 PM2.5、前两天的 PM2.5、当天的风速和湿度去预测当天 PM2.5。滞后特征的窗口不是越大越好对逐日数据3 到 7 天的滞后窗口足够超过 14 天基本是噪音。# 构造滞后特征和日历特征 for lag in [1, 2, 3, 7]: df[fPM2.5_lag{lag}] df[PM2.5].shift(lag) df[fPM10_lag{lag}] df[PM10].shift(lag) df[month] df[date].dt.month df[weekday] df[date].dt.weekday df[is_weekend] (df[weekday] 5).astype(int) # 删除前7行滞后特征为NaN的行 df df.dropna().reset_index(dropTrue)特征不是越多越好特别是加入风速和湿度后模型会出现一些奇怪的非线性效应。风速低且湿度高时 PM2.5 容易累积这两者的交互项有时比单独的风速和湿度更有效。参数说明shift(1) 是前一天shift(3) 是前三天weekday 从 0 到 65 和 6 对应周末。如果数据集跨年还要考虑春节这种对污染排放有显著影响的时间节点简单做法是加一个 is_spring_festival 布尔特征虽然粗糙但在预测精度上通常有可感知的提升。3. 空气质量评估模型从熵权法到综合指数评出可解释的排序评估城市空气质量最忌讳的就是拿一个黑匣子模型算出个总分然后直接排序。评委希望看到的是可解释的、有物理意义的评估过程。常见做法是先用熵权法或 AHP 确定六项污染物和气象指标的权重再用加权综合指数计算每个城市的得分最后按得分划分等级。省级优秀奖的论文里这块一般会配一张权重表、一张综合得分排序表。3.1 选型判断用 AHP 还是熵权法给指标定权重AHP 是主观赋权法依赖专家打分适合指标少、物理意义清晰的场景熵权法是客观赋权法根据数据本身的离散程度定权重指标差异越大权重越高。空气质量评估我更倾向熵权法原因很简单污染物浓度数据本身就包含污染源结构信息数据波动大的污染物往往对应污染事件频发权重应该更高AHP 的打分矩阵很容易被评委质疑主观性。熵权法的计算逻辑分三步归一化、算信息熵、算权重。归一化要区分正向指标和负向指标PM2.5 这类浓度越高越差的指标是负向的归一化时要用最大值减当前值再除以极差。def entropy_weight(df_norm): # df_norm: 已经归一化到[0,1]的指标矩阵行是样本列是指标 n, m df_norm.shape # 计算每个样本在指标j下的比重 p df_norm / df_norm.sum(axis0) # 计算信息熵 k 1 / np.log(n) e -k * (p * np.log(p 1e-12)).sum(axis0) # 计算差异系数和权重 d 1 - e w d / d.sum() return w逻辑说明p 是每个样本在某个指标上的占比信息熵越小说明该指标的数据越分散区分度越高权重应该越大。参数说明np.log 里加了 1e-12 是为了防止 p 为 0 时取对数报错。实际项目中熵权法算出来的权重可能和直觉不符比如某城市 CO 浓度全年都很低且平稳熵权法会给它很小的权重这其实是合理的因为它对空气质量排序几乎没有区分度。3.2 主成分降维与综合得分保留多少成分合适有些题目会额外要求做污染源结构分析这时可以把六项污染物浓度做主成分分析。主成分分析的直接产出是各主成分的载荷矩阵和方差解释率能看出哪几种污染物经常同步超标。比如第一主成分通常在 PM2.5、PM10、NO₂ 上有高载荷对应燃煤和机动车源第二主成分可能只在 O₃ 上有高载荷对应光化学污染。from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler # X为六项污染物浓度 scaler StandardScaler() X_scaled scaler.fit_transform(X) pca PCA(n_components0.85) # 保留方差解释率85% X_pca pca.fit_transform(X_scaled) print(pca.explained_variance_ratio_) print(pca.components_)参数说明n_components0.85 表示自动选择使累计方差解释率达到 85% 的主成分个数这个阈值是我个人经验对六项污染物通常能取到 2 到 3 个主成分太低会丢失信息太高则失去了降维的意义。注意要在标准化之后做 PCA否则浓度量级大的 PM10 会主导主成分方向。输出 components_ 矩阵时留意载荷值载荷的正负和大小是解释污染源结构的关键。3.3 综合指数与等级划分算出分数后怎么落到结论综合指数法把熵权法权重和标准化后的指标相乘再求和得到一个 0 到 1 之间的得分。分数本身没有绝对意义关键是排序和分级。这里有一个很多队伍会犯的错直接用线性加权算综合指数但 AQI 本身是非线性映射六项污染物浓度到 AQI 是分段线性直接用原始浓度线性加权会低估或高估某些城市。# 综合指数得分score_i sum(w_j * x_norm_ij) score (X_norm * w).sum(axis1) # 按分位数划分等级注意这里用的是样本分位数而非国标 q25, q50, q75 np.percentile(score, [25, 50, 75]) grade np.where(score q75, 较差, np.where(score q50, 中等, np.where(score q25, 较好, 好)))按样本分位数分级的好处是不同城市之间有可比性坏处是如果某年整体污染很重即使得分不高也可能被评为较差。这就是为什么最终定级最好结合 AQI 国标等级而不是纯看分位数。我的做法是两种都算分位数等级用来横向对比城市国标 AQI 等级用来描述城市自身的达标天数两者在论文里对应两张不同的表。4. 预测模型从 ARIMA 到随机森林逐日浓度预测的完整流程评估是打分预测是看未来。逐日空气质量预测的难点不是模型选得有多高级而是数据的时间依赖性和特征的有效性。常见路径是先跑一个时间序列基线模型 ARIMA再上机器学习模型如随机森林或 XGBoost最后做残差修正。如果数据量够大且想追求精度可以再加入 LSTM但对逐日数据、一两年的训练样本树模型往往比深度模型更稳。4.1 ARIMA 基线定阶、训练与预测窗口ARIMA 的作用是给后续机器学习模型提供一个对比基线。如果机器学习的精度还打不过简单的 ARIMA说明特征工程存在问题。ARIMA 定阶用 ACF 和 PACF 图结合 AIC 准则逐日 PM2.5 序列通常表现出一阶自相关和 7 日周期的季节性所以先尝试 ARIMA(1,1,1) 并对比加入季节项的结果。from statsmodels.tsa.arima.model import ARIMA # y为PM2.5日均浓度序列 model ARIMA(y, order(1, 1, 1)) model_fit model.fit() # 预测未来7天 forecast model_fit.forecast(steps7) print(model_fit.summary())逻辑说明order 里的三个参数分别对应自回归阶数 p、差分阶数 d、移动平均阶数 q。对非平稳的浓度序列先做一阶差分让它平稳所以 d 取 1。预测未来 7 天时ARIMA 只能基于历史值外推它看不到气象预报数据所以精度上限不高但作为基线足够了。参数说明如果 AIC 显示更低阶的 ARIMA(0,1,1) 或 ARIMA(2,1,2) 更好就换不要迷信默认参数。4.2 随机森林与 XGBoost输入输出设计与参数范围树模型的优势是能同时吃进历史浓度、气象因素和日历特征而且不需要对特征做标准化这对混合类型特征很友好。用随机森林做逐日预测输入特征用前面构造的滞后浓度、风速、湿度、温度、月份、是否周末输出是当天的 PM2.5 浓度。from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split from sklearn.metrics import mean_absolute_error, r2_score features [PM2.5_lag1, PM2.5_lag2, PM2.5_lag3, PM10_lag1, wind_speed, humidity, temperature, month, is_weekend] X df[features] y df[PM2.5] X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, shuffleFalse) rf RandomForestRegressor( n_estimators300, max_depth10, min_samples_leaf3, random_state42 ) rf.fit(X_train, y_train) y_pred rf.predict(X_test) print(MAE:, mean_absolute_error(y_test, y_pred)) print(R2:, r2_score(y_test, y_pred))逻辑说明shuffleFalse 非常关键时间序列数据不能随机打乱再切分否则训练集里包含测试集前后的信息测试集指标会虚高这就是典型的数据泄漏。参数说明n_estimators 取 300 以上可以让误差收敛max_depth10 防止树过深导致过拟合min_samples_leaf3 限制叶子节点最少样本数。这三个参数是我在实际项目里最常用的默认组合如果特征数很多可以调大 max_depth如果数据噪声很大可以调大 min_samples_leaf。4.3 滑窗验证与多步预测验证集不能偷看未来用 train_test_split 切分时间序列仍然不够严谨因为单次切分的测试集只有一段验证结果受这段时间的天气事件影响很大。我一般会做滑窗验证从训练集末尾取多个连续时间段作为验证集模型在每个验证集上重新训练并评估最后取平均误差。这里以滚动预测未来 7 天为例逐日预测比一次性预测 7 天更稳。def rolling_evaluate(model, X, y, train_days300, val_days7): preds, trues [], [] for start in range(train_days, len(X) - val_days, val_days): X_train X.iloc[start - train_days:start] y_train y.iloc[start - train_days:start] X_val X.iloc[start:start val_days] y_val y.iloc[start:start val_days] model.fit(X_train, y_train) preds.extend(model.predict(X_val)) trues.extend(y_val) return mean_absolute_error(trues, preds)逻辑说明start 不断向后滑动每次用前 300 天训练、后 7 天验证这样能模拟真实使用场景中模型持续更新的过程。参数说明train_days 和 val_days 按数据总长度调整一年数据可以设 train_days250、val_days7。如果数据量太少滑窗会切出很多短训练集这时可以重叠窗口但每段训练集至少要包含两个完整季节周期。5. 避坑城市空气质量评估与预测中最常见的五个翻车点做这类题模型算法其实都是现成的真正拉开差距的是谁能避开那些隐蔽的坑。以下五条都是实战中反复遇到的每条我都按现象、原因、解决的顺序写清楚其中前两条几乎每个队伍都会踩。5.1 现象一预测的 PM2.5 浓度出现负值现象模型预测结果里有负值比如 -12 μg/m³这在物理上不可能直接拉低均方根误差指标评委一眼就能看出问题。原因线性回归或某些集成模型没有输出范围约束当输入特征落在训练数据范围之外时预测值可能越过零线。解决对模型输出做后处理裁剪小于 0 的按 0 处理或者改用对数值做训练目标即对浓度取对数再预测最后取指数还原。后者效果更好因为浓度数据本身呈右偏分布对数变换能让模型更关注低浓度区间的精度。5.2 现象二训练集 R² 很高测试集一塌糊涂现象训练集 R² 达到 0.95测试集 R² 只有 0.4误差曲线在验证段突然失真。原因两种可能一是树模型过深导致过拟合训练集噪声二是时间序列切分时用了 shuffleTrue造成数据泄漏。解决先检查切分代码确认 shuffleFalse再把 max_depth 从默认的 None 降到 10 左右加上 min_samples_leaf 约束。如果过拟合仍然严重就用特征重要性排序删掉滞后 7 天以上的特征这些特征在训练集里可能碰巧与目标相关但本质是噪声。5.3 现象三沙尘暴和节假日让误差暴增现象模型整体误差尚可但三四月的某一天预测误差超过 200 μg/m³检查发现当天有沙尘暴春节期间误差也明显偏大。原因模型训练数据里沙尘暴样本极少模型根本没有见过 PM10 浓度飙到数千的场景春节期间排放结构变化大工厂停工、机动车减少源结构完全不同于工作日。解决在特征里加 is_dust_event 和 is_holiday 两个标志变量有历史标记的话直接加没有标记的话用 PM10/PM2.5 比值大于某个阈值识别沙尘事件。同时可以考虑对极端事件单独建模或直接剔除预测常规天气时精度更稳定。5.4 现象四PM10 和 PM2.5 高相关导致多重共线性现象随机森林的特征重要性显示 PM10_lag1 和 PM2.5_lag1 几乎平分重要性两个特征高度相关权重被分散模型解释困难。原因PM10 和 PM2.5 同源排放浓度曲线走势高度一致相关系数通常超过 0.8放在一起会产生共线性问题。对线性模型来说这会导致系数不稳定对树模型来说是特征重要性被稀释。解决评估模型中把 PM10 和 PM2.5 分成两类或者用 PM10-PM2.5 的差值作为粗颗粒物指标预测模型中只保留其中一个作为滞后特征优先保留 PM2.5因为它是预测目标本身的历史值。5.5 现象五AQI 与污染物浓度混用导致评估结论矛盾现象论文里评估用 AQI 等级预测用污染物浓度两个模块的结论对不上比如某城市 AQI 排名靠前但 PM2.5 浓度排名靠后。原因AQI 由六项污染物中的最差值决定一个城市可能臭氧很高但 PM2.5 很低评估排名反映的是短板效应预测排名反映的是单项浓度水平。解决在论文里明确区分两个对象评估结论基于综合指数相对排名预测结论基于 PM2.5 或 AQI 绝对数值在框架图里画出两条线。这样不仅逻辑自洽而且评委能看到你在用不同视角回答不同问题。6. 用模型融合收尾Stacking 与残差修正拿稳精度分预测模块做到随机森林已经是很多队伍的终点但要拿到省级优秀奖级别的精度还得加一层模型融合。模型融合在竞赛里是常规动作在校赛和省赛里反而是很多人不敢碰的加分项。本节给出一个可复现的融合方案以及让图表看起来更专业的两个验收指标。6.1 Stacking 框架让随机森林和 XGBoost 互相纠错用随机森林和 XGBoost 作为基模型以它们的预测结果作为新特征训练一个线性回归作为元模型。基模型之间相关性越低融合效果越好两个树模型相关性偏高但 XGBoost 对噪声的鲁棒性略好随机森林对特征交互的捕捉略好互为补充。from xgboost import XGBRegressor from sklearn.linear_model import LinearRegression from sklearn.model_selection import KFold # 基模型定义 base_models [ RandomForestRegressor(n_estimators300, max_depth10, min_samples_leaf3, random_state42), XGBRegressor(n_estimators300, max_depth4, learning_rate0.05, reg_lambda1.0) ] # 5折交叉生成元特征防止基模型在训练集上过拟合导致元模型失效 kf KFold(n_splits5, shuffleFalse) meta_features np.zeros((len(X_train), len(base_models))) for i, model in enumerate(base_models): for train_idx, val_idx in kf.split(X_train): fold_model model.__class__(**model.get_params()) fold_model.fit(X_train.iloc[train_idx], y_train.iloc[train_idx]) meta_features[val_idx, i] fold_model.predict(X_train.iloc[val_idx]) # 元模型 meta_model LinearRegression() meta_model.fit(meta_features, y_train) # 测试集预测 test_meta np.column_stack([m.predict(X_test) for m in base_models]) final_pred meta_model.predict(test_meta)逻辑说明基模型在 KFold 的每一折上重新训练并预测验证折这样元模型的输入特征不会包含基模型“见过”的样本避免元模型学到基模型的记忆性错误。参数说明shuffleFalse 保持时间顺序n_splits5 是常用默认数据多可以到 10数据少就 3。XGBoost 的 reg_lambda1.0 是 L2 正则防止它在小样本上过拟合。6.2 残差修正用误差序列再学一轮融合模型在污染事件日的残差通常有正相关性即前一天低估了后一天往往也低估。利用这一点把融合模型的预测残差作为目标用简单的线性回归或随机森林学习残差与气象特征的关系最后把残差预测值加回融合预测。这一步在空气质量这类强自相关的时间序列里效果明显。# 第一轮预测 pred_train meta_model.predict(meta_features) residual y_train - pred_train # 用同一批特征预测残差 resid_model RandomForestRegressor(n_estimators200, max_depth5, min_samples_leaf5) resid_model.fit(X_train, residual) # 最终预测 融合预测 残差预测 resid_pred resid_model.predict(X_test) final_adjusted final_pred resid_pred逻辑说明残差模型学习的是系统误差比如湿度高时模型系统性低估 PM2.5 累积速度这种误差模式在气象特征上是有规律的。参数说明max_depth 控制在 5 以下因为残差里的可学习信号比原始浓度弱得多树太深会把残差里的噪声也学进去反而放大误差。6.3 验证与呈现误差分解图和特征重要性的正确用法模型做完之后论文里最容易加分也最容易被忽视的是误差的可视化呈现。画时间序列预测对比图时不要只画预测线和真实线要画出误差柱状图并且用颜色标注误差超过 50 μg/m³ 的日期旁边注明对应天气事件。特征重要性图不要用默认的散点要按数值排序画水平条形图并在每个特征名后标注物理含义比如 PM2.5_lag1 标注为“前一天 PM2.5 浓度残留”。验收时我自己会按三个指标卡标准R² 达到 0.75 以上MAE 控制在 20 μg/m³ 以内7 天预测中第 1 天误差不超过第 7 天误差的一半。如果第 1 天误差就很大说明滞后特征构造有问题不是模型问题。这是我从那次省级优秀奖之后一直沿用的习惯先看分时段误差走势图再决定改特征还是改模型而不是一上来就换算法。希望帮到你。本文还有配套的精品资源点击获取
返回列表