空气质量预测实战:从数据预处理到SHAP模型解释
简介一份基于机器学习的空气质量预测数据挖掘实战资料适合机器学习初学者、数据挖掘课程学生以及需要快速搭建预测项目的开发者。项目以空气质量污染数据集为对象借助Python环境及Jupyter笔记本完成数据读取、缺失值清洗、特征分析、模型训练与结果评估完整覆盖了数据挖掘的常规流程。资源共包含3个文件分别是ipynb代码笔记本、html结果报告和csv原始数据文件压缩包整体约1.39MB轻量便捷目前已有128人浏览学习可作为入门参考。使用时可以对照代码理解数据预处理和建模思路也可以直接查看HTML报告中的图表与结论即使不运行程序也能掌握项目脉络。数据集与代码配套提供方便替换数据或调整参数进行迁移练习对完成类似预测任务很有帮助。1. 空气质量预测模型数据挖掘实战里最值得复现的一类项目做数据挖掘最怕的不是算法不会调而是拿到一份数据集不知道从哪下手。这份空气质量预测模型资源正好卡在这个痛点上数据集、完整分析代码、可交互的 HTML 报告三样都齐拿到就能跑通一整套“数据预处理 → 特征工程 → 建模评估 → 结果解释”的流程。对正在学机器学习但缺实战项目练手的人或者做环境数据分析想参考建模套路的人来说它是一个能直接抄作业的样本。核心是那份 updated_pollution_dataset.csv 和配套的 分析2.ipynb跑完你会看到数据挖掘在回归预测问题上到底是怎么一步步落地的。这比单看算法理论有用得多。2. 数据预处理与探索性分析先把 CSV 里的“脏东西”揪出来拿到 updated_pollution_dataset.csv 之后别急着建模第一步永远是搞清楚数据长什么样。很多人在 Kaggle 上跑过类似项目上来就 pd.read_csv() 然后直接切特征训练结果模型效果差还不知道为什么。常见原因是没检查数据类型、缺失值和离群点这三个问题会一路传导到模型输出。2.1 先看数据类型和缺失值pandas 的 dtype 和 isna 组合拳我一般拿到 CSV 先跑一遍 info() 和 isna().sum()这两行代码能快速暴露大部分问题。空气质量数据集里最常见的坑是时间列被读成字符串污染物浓度列混入 NA 之类的不规范缺失标记或者 PM2.5 等列里出现 -999 这类哨兵值。import pandas as pd import numpy as np df pd.read_csv(updated_pollution_dataset.csv, encodingutf-8) print(df.info()) print(df.isna().sum()) # 把常见的哨兵值替换为 NaN避免脏数据进入统计量 sentinel_values [-999, -99, 9999] df df.replace(sentinel_values, np.nan) # 数值列统一转 float时间列后面单独处理 for col in [PM2.5, PM10, SO2, NO2, CO, O3]: if col in df.columns: df[col] pd.to_numeric(df[col], errorscoerce)逻辑说明info() 看每列的非空计数和数据类型isna().sum() 看缺失分布replace 这一步是为了把数据源里用来占位的极端值还原成真正的缺失。errorscoerce 保证脏字符串被转成 NaN 而不是直接抛异常中断流程。参数说明哨兵值列表里的 -999、-99、9999 是环境监测数据里比较常见的占位写法实际用的时候先跑一次 df.describe() 看最小值和最大值如果出现明显不合理的数值就加进哨兵列表。errorscoerce在 pandas 里表示转换失败时填入 NaN真遇到数据里有空格或中文逗号的情况也不会让程序崩掉。2.2 时间字段解析与索引设置pd.to_datetime 与 UTC 偏移空气质量的时序属性很强时间索引没设对后面构造滞后特征、滚动窗口全都会错位。很多新手在这里只把时间列解析成 datetime 就完事但实际上还要考虑时间频率是否均匀、时区是否需要统一。# 假设时间列名是 date格式类似 2023-01-01 08:00:00 df[date] pd.to_datetime(df[date], format%Y-%m-%d %H:%M:%S) # 检查时间间隔是否均匀有重采样需求时先归一到小时粒度 df df.set_index(date).sort_index() print(df.index.is_monotonic_increasing) hourly df.resample(1H).mean() print(hourly.head())逻辑说明把时间列转成 datetime 并设为索引后时间序列的排序、重采样、滑动窗口才能按时间语义执行。is_monotonic_increasing 检查索引是否按时间递增数据乱序时后续的 shift 会错位。参数说明format 参数按实际 CSV 里时间字符串的写法定常见的是 %Y-%m-%d %H:%M:%S如果数据里只有日期没有小时就要考虑 resample 的粒度是否要改成 1D。resample(1H).mean() 把非整点数据归到整点这一步做不做取决于原始采样频率如果原始数据本身就是小时级直接跳过重采样。2.3 离群点与分布初探describe、箱线图、对数变换的判断污染物浓度数据经常是长尾分布个别极端日比如沙尘暴当天的 PM10 飙升到上千会把均值拉高、把线性模型的系数带偏。处理离群点不能无脑删除得先看业务合理性。import matplotlib.pyplot as plt desc df[[PM2.5, PM10, SO2, NO2, CO, O3]].describe() print(desc) # 用 IQR 方法找出疑似离群点但只记录不直接删 Q1 df[PM2.5].quantile(0.25) Q3 df[PM2.5].quantile(0.75) IQR Q3 - Q1 outlier_mask (df[PM2.5] Q1 - 1.5 * IQR) | (df[PM2.5] Q3 1.5 * IQR) print(PM2.5 疑似离群点数量:, outlier_mask.sum()) # 看分布形态右偏严重的话考虑对数变换 df[PM2.5_log] np.log1p(df[PM2.5])逻辑说明describe 先看均值和中位数差异均值远大于中位数说明右偏。IQR 方法只做标记不自动删数据因为空气质量数据里高浓度值可能是真实污染事件删掉反而丢失信息。对数变换压缩长尾让模型更容易拟合。参数说明1.5 倍 IQR 是标准做法但对空气质量这种重尾分布可以放宽到 3 倍 IQR。np.log1p 是 log(x1)专门处理存在 0 值的数据避免 log(0) 报错。3. 特征工程与数据集划分把时间、滞后、类别转成模型能吃的形状数据挖掘项目的成败很大程度在特征工程。空气质量预测除了原始污染物浓度时间特征和滞后特征往往能带来显著的精度提升。这份资源里的分析 2.ipynb 在特征工程上做了不少文章下面把关键步骤拆开讲。3.1 时间特征构造小时、星期、月份、节假日空气污染呈现明显的周期性早晚高峰 PM2.5 升高、周末车流量变化导致污染物浓度波动、冬季采暖季和夏季臭氧高发期各有规律。这些周期信号不喂给模型模型就只能靠原始浓度硬猜。df[hour] df.index.hour df[weekday] df.index.weekday df[month] df.index.month df[is_weekend] (df.index.weekday 5).astype(int) # 周期特征用正弦余弦编码避免 23 点和 0 点被模型当作距离很远 df[hour_sin] np.sin(2 * np.pi * df[hour] / 24) df[hour_cos] np.cos(2 * np.pi * df[hour] / 24) print(df[[hour, hour_sin, hour_cos]].head())逻辑说明小时这种循环特征直接当数值喂给树模型还行但喂给线性回归会出现“23 点 和 0 点差 22 个小时”的荒谬距离。正弦余弦编码把小时映射到单位圆上保留周期性。参数说明24 是小时周期weekday 的周期是 7month 的周期是 12做编码时把周期分母换成对应值即可。is_weekend 是二值特征有些数据集里还有节假日列有的话直接拼上。3.2 滞后特征与滑动窗口shift 和 rolling 的边界问题空气质量有很强的自相关性今天的 PM2.5 和昨天、前天高度相关。构造滞后特征就是把这个“前因后果”显式地交给模型。但这一步最容易翻车用 shift 之前必须先按时间排序否则滞后的是乱序数据。df df.sort_index() for lag in [1, 3, 6, 24]: df[fPM2.5_lag_{lag}] df[PM2.5].shift(lag) # 24 小时滑动平均代表历史污染水平 df[PM2.5_rolling_24] df[PM2.5].rolling(window24).mean() # 滑动窗口的标准差刻画波动程度 df[PM2.5_rolling_std_24] df[PM2.5].rolling(window24).std() # 删除因为 shift 和 rolling 产生的 NaN 行避免模型训练时报错 df_clean df.dropna() print(df_clean.shape)逻辑说明shift(lag) 把过去第 lag 个时刻的值拉到现在这一行rolling(window24).mean() 计算过去 24 小时的均值。滞后特征让模型能用历史值做回归滑动统计量捕捉趋势和波动。参数说明lag 的选择要看数据粒度和业务周期小时级数据常用 1、3、6、24日级数据常用 1、2、7。rolling window 设 24 对应小时级数据的完整一天。dropna 必须放在最后因为 shift 和 rolling 会引入前几个时刻的 NaN。3.3 类别特征编码与特征选择LabelEncoder、卡方检验空气质量数据里除了数值型污染物浓度通常还有风向、天气状况等类别特征。把字符串类别直接喂给 sklearn 模型会报错需要先编码。风向是环形的和小时特征同理最佳做法是把角度拆成 sin 和 cos 两个分量。from sklearn.preprocessing import LabelEncoder # 风向如果是 16 方位字符串先映射成角度再拆分量 wind_dir_map {N: 0, NNE: 22.5, NE: 45, ENE: 67.5, E: 90, ESE: 112.5, SE: 135, SSE: 157.5, S: 180, SSW: 202.5, SW: 225, WSW: 247.5, W: 270, WNW: 292.5, NW: 315, NNW: 337.5} df[wind_deg] df[wind_dir].map(wind_dir_map).astype(float) df[wind_sin] np.sin(np.deg2rad(df[wind_deg])) df[wind_cos] np.cos(np.deg2rad(df[wind_deg])) # 天气描述这种纯类别列用 LabelEncoder le LabelEncoder() df[weather_encoded] le.fit_transform(df[weather])逻辑说明风向的 0 度和 360 度是同一个方向直接编码成数字会切断周期关系。映射成角度再拆成 sin/cos 后模型能理解风向的连续性。LabelEncoder 适合树模型喂线性回归建议用 OneHotEncoder。参数说明wind_dir_map 里的角度对应 16 方位标准气象定义数据是 8 方位就只需保留 N、NE、E、SE、S、SW、W、NW 八个键。weather 列如果类别太多可以先分组再编码比如把“晴”、“多云”归为“晴好天气”把“雨”、“雪”归为“降水天气”。4. 模型构建与训练线性回归、随机森林、XGBoost 的选型与参数数据干净、特征做完了才轮到算法登场。很多人一上来就跑 XGBoost结果发现效果还不如线性回归就是因为没有先用简单模型做基线。基线模型的意义是给后续复杂模型设一个及格线连及格线都超不过的模型直接淘汰。4.1 数据划分时间顺序切分 vs 随机切分这一步决定是否泄漏时序预测任务不能直接调用 train_test_split 默认的随机切分。随机切分会把未来数据混进训练集模型在验证集上的分数虚高一旦上线做真实预测立刻打回原形。这份数据是时间序列要么按时间顺序切分要么用 TimeSeriesSplit 做交叉验证。from sklearn.model_selection import train_test_split feature_cols [c for c in df_clean.columns if c not in [PM2.5, PM10, date, weather]] X df_clean[feature_cols] y df_clean[PM2.5] # 按时间顺序切分前 80% 训练后 20% 验证 split_idx int(len(X) * 0.8) X_train, X_test X.iloc[:split_idx], X.iloc[split_idx:] y_train, y_test y.iloc[:split_idx], y.iloc[split_idx:] print(ftrain 时间范围: {X_train.index.min()} ~ {X_train.index.max()}) print(ftest 时间范围: {X_test.index.min()} ~ {X_test.index.max()})逻辑说明iloc 按行位置切分保持了时间的先后顺序。这样做是模拟真实场景——用过去预测未来。打印时间范围是确认切分点经常有人切完才发现测试集时间比训练集还早。参数说明0.8 的比例在数据量足够时合理数据量少于一万条可以调到 0.85。如果数据有明显的季节性周期切分点尽量卡在季节边界附近避免训练集只含夏天、测试集只含冬天。4.2 基线模型线性回归与小规模标准化线性回归是判断特征工程做得是否合格的标尺。如果特征和 PM2.5 之间确实存在线性关系线性回归的 R² 应该不会太差。数值特征量纲差异大时先做标准化否则梯度下降类模型收敛慢。from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LinearRegression from sklearn.pipeline import Pipeline # 只用数值特征含类别编码后的列 numeric_cols [c for c in feature_cols if c not in [weather_encoded]] pipeline_lr Pipeline([ (scaler, StandardScaler()), (lr, LinearRegression()) ]) pipeline_lr.fit(X_train[numeric_cols], y_train) train_score pipeline_lr.score(X_train[numeric_cols], y_train) test_score pipeline_lr.score(X_test[numeric_cols], y_test) print(f线性回归 train R²: {train_score:.4f}, test R²: {test_score:.4f})逻辑说明Pipeline 把标准化和线性回归打包fit 的时候先对训练集计算均值和标准差再把同样的变换应用到测试集。这样避免用测试集的统计量做标准化而造成信息泄漏。参数说明StandardScaler 默认把数据转成均值 0、方差 1对量纲差几个数量级的特征比如 CO 浓度是 0.xPM2.5 是 100尤其重要。如果换用 RandomForest 这类树模型标准化反而没必要树模型对单调变换不敏感。4.3 树模型与集成模型随机森林和 XGBoost 的关键参数当线性回归做不动的时候就该上树模型了。随机森林适合做特征重要性的初步评估XGBoost 在梯度提升框架下后发制人。这两个模型在网格调参时参数策略完全不同下面给的参数组合是我的常用起点不是最优解但胜在不容易过拟合。from sklearn.ensemble import RandomForestRegressor from xgboost import XGBRegressor rf RandomForestRegressor( n_estimators300, max_depth12, min_samples_leaf3, random_state42 ) rf.fit(X_train[feature_cols], y_train) rf_test_score rf.score(X_test[feature_cols], y_test) print(f随机森林 test R²: {rf_test_score:.4f}) xgb XGBRegressor( n_estimators300, learning_rate0.05, max_depth6, subsample0.8, colsample_bytree0.8, reg_lambda1.0, random_state42 ) xgb.fit(X_train[feature_cols], y_train) xgb_test_score xgb.score(X_test[feature_cols], y_test) print(fXGBoost test R²: {xgb_test_score:.4f})逻辑说明RandomForest 的 min_samples_leaf3 防止树在训练集上记得太死max_depth12 给足深度但留了缓冲。XGBoost 里 learning_rate 调低后必须配合更多棵树subsample 和 colsample_bytree 是防过拟合的双保险。参数说明n_estimators300 在小样本数据集上是够用的数据量超过十万条可以加到 500 以上。reg_lambda1.0 是 XGBoost 的 L2 正则默认值就是这个过拟合严重时往上调到 2~5。random_state 固定是为了让结果可复现。5. 评估与调参MAE、R² 之外交叉验证和 GridSearchCV 的常见坑模型训练完不能只看 R² 就收工。空气质量预测这种连续值回归任务需要看多维度指标才能判断模型是“整体差”还是“极端值差”。这个章节既是评估方法也是整篇最值得看的避坑内容。5.1 评估指标怎么选MAE、RMSE、R² 的组合判断R² 高不代表预测准它只代表模型解释了方差的比例。MAE 反映平均绝对误差RMSE 对大的误差更敏感。空气质量预测里如果存在偶尔飙升的重度污染日RMSE 会远大于 MAE说明模型极端值拟合差。from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score import numpy as np y_pred xgb.predict(X_test[feature_cols]) mae mean_absolute_error(y_test, y_pred) rmse np.sqrt(mean_squared_error(y_test, y_pred)) r2 r2_score(y_test, y_pred) print(fXGBoost MAE: {mae:.2f}, RMSE: {rmse:.2f}, R²: {r2:.4f}) # 看误差分布预测值明显偏低还是偏高 error y_test - y_pred print(误差均值:, error.mean(), 误差中位数:, np.median(error))逻辑说明MAE 和 RMSE 同时看如果 RMSE 明显大于 MAE说明个别样本的误差很大模型在拖后腿。误差均值偏离 0 越远模型系统性偏差越严重均值接近 0 但中位数偏负则说明多数样本被低估。参数说明误差的单位和 PM2.5 一致如果 PM2.5 均值在 80 左右MAE 在 15 以内算不错。不同城市污染水平不同这个数字要结合数据分布看别拿别的城市的绝对数值对比。5.2 交叉验证与 GridSearchCV时间序列里不能直接 KFold时序数据做交叉验证有一个经典的翻车点直接用 KFold 随机打乱数据让未来的数据进入训练集导致验证分数虚高。正确的做法是 TimeSeriesSplit它保证每一折的训练数据都在验证数据之前。from sklearn.model_selection import TimeSeriesSplit, GridSearchCV tscv TimeSeriesSplit(n_splits5) param_grid { max_depth: [6, 8, 10], learning_rate: [0.03, 0.05], reg_lambda: [1.0, 2.0] } grid_search GridSearchCV( estimatorXGBRegressor(n_estimators200, subsample0.8, colsample_bytree0.8, random_state42), param_gridparam_grid, cvtscv, scoringneg_mean_absolute_error, n_jobs-1 ) grid_search.fit(X_train[feature_cols], y_train) print(best params:, grid_search.best_params_) print(best score (负 MAE):, grid_search.best_score_)逻辑说明TimeSeriesSplit(n_splits5) 把数据切成 5 份第一折用第 1 份训练、第 2 份验证第二折用第 12 份训练、第 3 份验证以此类推。scoring 用负 MAE 是因为 GridSearchCV 里分数越大越好而 MAE 越小越好取负号就对齐了。参数说明n_splits5 对小时级数据来说偏少可以调成 8 或 10但折数越多训练越慢。param_grid 里的组合 3×2×212 组每组 5 折共 60 次训练数据量大时把 n_jobs 调成 -1 充分利用多核。5.3 空气质量预测避坑手册四个真实翻车现场这一节把我在类似项目里踩过的坑集中写出来每一条都是现象、原因、解决三段式照着排查能省半天时间。翻车现场一预测结果整体比真实值偏低一个固定数值。现象误差均值是 -10 左右MAE 不高但系统性偏差明显。 原因目标变量做了对数变换之后忘了在预测结果上做指数还原或者标准化时用了 fit_transform 的全局统计量导致训练集信息泄漏。 解决检查数据处理流程里是否有 np.log1p / np.expm1 的对应操作有变换就必须有反变换。标准化必须在 train_test_split 之后做而且只 fit 训练集。翻车现场二时间索引没有排序就做 shift。现象滞后特征全乱套模型训练分数很高但验证分数极差。 原因read_csv 之后没有 set_index 和 sort_indexshift(1) 拿到的“前一天”实际上是文件里上一行。 解决在构造任何滞后特征之前先执行 df df.sort_index() 并打印 df.index.is_monotonic_increasing 确认。翻车现场三GridSearchCV 的验证分数比实际测试分数高一大截。现象网格搜索里 best_score 是 0.92测试集 R² 只有 0.75。 原因用了默认 KFold 而不是 TimeSeriesSplit未来的数据混进了训练折。或者特征里包含目标变量的滞后值而滞后值在测试集拆分时因为索引错位泄漏到了训练集。 解决重跑交叉验证确认 cv 参数是 TimeSeriesSplit。检查特征列里是否有目标变量的 shift 值这些特征本身就是带泄漏的评估时一定要保证训练和测试在时间上完全分开。翻车现场四PM10 的预测效果比 PM2.5 差很多模型似乎“失灵”。现象同一套特征和模型PM2.5 的 R² 有 0.85PM10 只有 0.6。 原因PM10 受沙尘、扬尘等随机事件影响更大自相关性本来就弱滞后特征携带的信息量少。这不是模型坏了是任务本身的预测难度不同。 解决不要用同一套特征硬套所有污染物单独做一次特征筛选给 PM10 增加风速、风向的交互项或者把训练目标改成对数变换后的 PM10 再看指标。6. 进阶验证滚动预测、残差分析与 SHAP 解释模型预测模型在测试集上分数达标只是第一步真实部署时要面对的是连续多步预测和模型可解释性。这个章节讲三个我常用的验证技巧能帮你判断模型是真能打还是只会在测试集上“刷题”。6.1 滚动预测验证模拟真实上线逻辑单次预测测试集和“每天预测未来 24 小时”是两种任务。滚动预测的做法是用前 N 小时的数据预测下一小时然后把真实值并入历史再进行下一步预测。这个过程更贴近空气质量的真实预报场景。# 取最后 240 小时做滚动验证 history X_train[feature_cols].tail(240).copy() history_y y_train.tail(240).copy() roll_preds [] for t in range(24 * 7, len(X_test[feature_cols])): if t len(X_test[feature_cols]): x_input X_test[feature_cols].iloc[t:t1] pred xgb.predict(x_input)[0] roll_preds.append(pred) # 注意这是简化示范真实滚动验证要逐步更新滞后特征逻辑说明简化版滚动预测直接对测试集逐点预测真正的滚动验证需要在每步预测后用真实值回去更新滞后特征列代码复杂度高很多但验证结果更接近部署场景。参数说明24 * 7 表示预热一周的数据才开始预测时间跨度可以根据实际需要调整。这种验证方式得到的 MAE 会略高于单次预测属于正常现象别看到分数下降就觉得模型坏了。6.2 残差分析找模型系统性盲区残差是真实值减预测值把残差按时间画出来或者对特征维度聚合能看出模型在哪些场景下稳定失效。residual y_test - xgb.predict(X_test[feature_cols]) # 按小时聚合残差 resid_by_hour pd.DataFrame({hour: X_test.index.hour, residual: residual}) print(resid_by_hour.groupby(hour)[residual].mean()) # 残差绝对值最大的 10 个样本 idx_top_err residual.abs().nlargest(10).index print(X_test.loc[idx_top_err, [PM2.5_lag_1, hour]])逻辑说明残差按小时聚合能看出来早高峰或夜间是不是系统性地预测偏大或偏小这是特征工程改进的线索。找出最大误差样本看它们的特征值往往能发现模型对某些极端天气或突变工况没有感知。参数说明groupby 的 hour 列可以在构造特征时用 X_test.index.hour 生成。nlargest(10) 取前 10 个误差最大的索引核心目的是剖析样本不是筛异常值别顺手把她们删了。6.3 SHAP 解释模型给预测结果一个说法机器学习模型的黑匣子问题在环境领域尤其敏感业务方会问“为什么预测明天 PM2.5 是 120”。SHAP 值能给出每个特征对单次预测的贡献量是目前做模型可解释性最直接的工具。import shap explainer shap.TreeExplainer(xgb) shap_values explainer.shap_values(X_test[feature_cols].iloc[:100]) shap.summary_plot(shap_values, X_test[feature_cols].iloc[:100])逻辑说明TreeExplainer 专为树模型设计速度比通用 KernelExplainer 快几个数量级。summary_plot 展示每个特征在所有样本上的 SHAP 值分布直观看到 PM2.5_lag_1 这类滞后特征对预测的拉动方向。参数说明iloc[:100] 是限制样本数SHAP 计算量大样本多的时候先取子集看态势。如果项目报告里需要特征重要性排名也可以导出 shap_values 后用 mean(|shap|) 排序比内置的 feature_importances_ 更稳定。做数据挖掘项目我一直坚持“测试集分数只是中间产物滚动验证和残差分析才决定模型能不能用”。从那以后我每次跑完一个模型都会强制走一遍滚动预测和 SHAP 解释这两步流程哪怕结果不理想也比模型在测试集上虚高强得多。这套流程配合这份资源里的 updated_pollution_dataset.csv 和 分析2.ipynb完整跑一遍能帮你把从数据预处理到模型解释的每个环节都打通。希望帮到你。本文还有配套的精品资源点击获取