线性回归预测PM2.5：从特征工程到模型评估的完整实践指南

简介：面向机器学习初学者的PM2.5预测大作业项目，基于合肥地区历史空气质量月均值数据，使用线性回归模型完成建模与预测，涵盖矩阵运算及梯度下降公式的具体实现。资源共21个文件，压缩包约2.58MB，其中12个CSV为训练与测试数据集，3个Python脚本承载数据预处理、模型训练和预测评估流程，另有npy格式模型参数、README说明文档及结果图像。项目以完整代码结合真实数据的形式，展示了从数据读取、特征构造到梯度下降求解、结果可视化的全流程，适合课程设计、期末作业或入门练习参考。目前已有280人学习浏览，文件结构清晰，便于按需取用各模块。

1. 线性回归预测 PM2.5：为什么课程大作业都爱选它

交机器学期末大作业的人常陷入“模型越高级越好”的误区，可在答辩现场被追问参数时，黑匣子最容易露馅。基于线性回归的 PM2.5 预测是个反直觉的好选择：模型简单到能手推公式，PM2.5 数据里又塞满了缺失值、异方差和时间自相关，正够把数据预处理、特征工程、模型评估这些核心知识点完整过一遍。这篇文章按我实际处理这类大作业的顺序展开：先理特征，再建模型，然后看评估，最后把我踩过的坑原样摆出来，帮你少绕几圈。适合正在写课程大作业、又不想只靠跑通别人源码交差的人。

2. 数据准备从哪下手：把 PM2.5 的“记忆”变成特征

2.1 字段审查：先确认手里是真实浓度还是空气质量指数

数据是起点，但我见过不少同学在起点就翻车——把 PM2.5 的 AQI 分指数当成浓度值直接喂给模型。常见的大作业数据源有两类：一类是 UCI 的 Air Quality 数据集，字段包括 CO、NO2、O3 等气体浓度以及温度湿度；另一类是城市环保站点公布的小时级历史数据，字段里除了污染物浓度还有气压、风速、风向、降雨。拿到表先别急着建模，用下面这段代码看清楚结构：

import pandas as pd df = pd.read_csv("air_quality_hourly.csv", parse_dates=["date"]) df = df.sort_values("date").set_index("date") print(df.head()) print(df.describe()) print(df.isna().sum())

head() 看的是字段名和取值量级，describe() 看每个列的均值、极值和分位数，isna().sum() 统计缺失规模。以小时级数据为例，一列正常的 PM2.5 浓度，均值应该在 30～80 µg/m³，最大值往往能到 300 甚至更高；如果你看到均值只有 20 上下、最大值不到 100，就要怀疑它是不是已经被换算成了 AQI 分指数，或者是某次数据清洗时做过归一化。这种错误藏得最深，因为它不报错，只表现为模型怎么调都学不好。

2.2 缺失值处理：别让 dropna 把夜间样本全扔了

空气质量站点数据有个特点：夜间经常出现连续几小时的校准停测，冬季设备结露也会造成长时间缺失。很多现成源码里直接写一行 df.dropna()，我现在一看到这种写法就头大。你可以自己做个统计：如果原始数据有三万行，dropna 之后只剩下两万行，丢失的几乎都是夜间和降雨时段，模型学到的就是“白天+晴天偏好”。正确的做法是按缺失成因分开处理。

# 污染物传感器短时间缺失，用时间插值 df["pm2_5"] = df["pm2_5"].interpolate(method="time") # 气象字段零星缺失，用前后向填充 df["temp"] = df["temp"].ffill().bfill() df["humidity"] = df["humidity"].ffill().bfill() df["pressure"] = df["pressure"].ffill().bfill() df["wind_speed"] = df["wind_speed"].ffill().bfill() print(df.isna().sum())

interpolate(method="time") 会用前后时间点按时间间隔比例插值，适合连续缺失不超过几个小时的缺口；ffill().bfill() 适合把某一小时偶发的空值用最近邻时刻补上。之后如果目标列还有残留缺失，再去按行剔除。这样处理完再统计 isna()，应该只剩插值补不了的零头。如果你拿到的数据缺了整整一天，建议直接把那一整段标为无效剔除，因为任何插值在长时间缺口上都是猜。

提示：pm2_5 是代码里的目标列名，阅读时可以直接把它当作 PM2.5 浓度的缩写。后面构造的 pm2_5_lag24、pm2_5_roll7 同理。

2.3 时间特征和滞后特征：线性回归也能“记住”昨天

线性回归没有循环结构，但它可以靠特征构造出“记忆”。PM2.5 浓度有很强的自相关：今天的浓度和 24 小时前同一小时的浓度高度相关；冬季采暖季的浓度又和过去一周的累积趋势有关。日期时间列本身不能直接进模型，所以我们要把它拆成小时、月份、星期几，再构造滞后特征。

df["hour"] = df.index.hour df["month"] = df.index.month df["dayofweek"] = df.index.dayofweek # 24小时前同刻浓度，体现日周期自相关 df["pm2_5_lag24"] = df["pm2_5"].shift(24) # 过去一周（7*24小时）平均浓度，体现累积污染 df["pm2_5_roll7"] = df["pm2_5"].rolling(7 * 24).mean() df = df.dropna() print(df.shape)

这里 shift(24) 的前提是数据为小时粒度，每条记录代表一小时；如果你的数据是 3 小时间隔，24 小时前就要写成 shift(8)。滚动均值 rolling(7 * 24) 的窗口是 168 个小时，也就是完整一周。窗口越长，历史信息越平滑，但越滞后；窗口过短又挡不住日周期波动，一周是我在污染数据上的默认值。构造完滞后特征后，最初 7 天因为缺少滚动窗口会变成空值，这里 dropna() 扔掉的是边界样本，而不是中间数据，可以放心。

2.4 切分数据：时间序列的测试集不能随机抽样

这一点怎么做都不嫌强调。很多网上流传的源码会把数据直接丢给 train_test_split，这在时间序列任务里是严重的概念错误：随机抽样会把未来某一周的数据混进训练集，模型等于“偷看了答案”。正确的做法是按时间顺序切分。

train = df.loc[df.index < "2020-01-01"] test = df.loc[df.index >= "2020-01-01"] print(train.shape, test.shape)

切分点按你自己的数据时间范围选择，原则是训练集在前、测试集在后，而且最好让测试集覆盖完整的一周以上，包含工作日和周末。污染数据有明显的星期效应：周一到周五通勤排放高，周末低，测试集如果只有三天，很可能恰好全是晴天或全是周末，评估出来的 R2 会虚高或虚低。做完这一步，数据侧的工作才算收尾。

3. 建模与调参：一条 Pipeline 走完线性回归到正则化

3.1 先跑通最小可用模型：别让“高级模型”变成黑匣子

拿到现成大作业源码时，最常见的翻车方式是把项目里一堆模型逐个跑一遍，哪个分数高交哪个，结果答辩时说不清任何参数。我的习惯是反着来：先跑最朴素的多元线性回归，把它当作基线，再往上做正则化和特征增强。下面构造 Pipeline，把特征标准化和线性回归打包在一起：

from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LinearRegression features = ["temp", "pressure", "humidity", "wind_speed", "co", "hour", "month", "dayofweek", "pm2_5_lag24", "pm2_5_roll7"] X_train = train[features] y_train = train["pm2_5"] pipe_lr = Pipeline([ ("scale", StandardScaler()), ("lr", LinearRegression()) ]) pipe_lr.fit(X_train, y_train) print("train R2:", pipe_lr.score(X_train, y_train)) print("test R2:", pipe_lr.score(test[features], test["pm2_5"]))

Pipeline 的意义不只是看着整洁。它保证 fit 时在训练集上计算标准化参数，predict 时用同一套参数变换测试集，避免你在手动处理时不小心用测试集的均值、方差去缩放训练集，那属于标签泄漏的一种隐蔽形态。LinearRegression 在 sklearn 里用的是最小二乘闭式解，不需要梯度下降，但如果数据量到了百万级，Pipeline 以后想换 SGDRegressor 也是一行替换的事。第一次跑出来测试 R2 在 0.3～0.6 之间都正常，别急着嫌低，先继续往下加正则化。

跑通以后做一个快速自检：如果 train R2 接近 1 而 test R2 为负，优先怀疑特征泄漏或目标列引用错误；如果两者都低到 0.2 以下，回头看特征列表里有没有列名拼错，以及 2.2 节缺失值补完以后样本量还够不够。这一步能把后面调参的无效劳动省掉一半。

3.2 岭回归和拉索回归：什么时候用、alpha 怎么选

特征一到十个以上，共线性问题就来了。pm2_5_lag24 和 pm2_5_roll7 本身就是同一列数据的函数，相关性极高，普通最小二乘的解会变得不稳定：系数绝对值巨大，正负号飘忽。岭回归（Ridge）通过 L2 惩罚限制系数平方和，拉索回归（Lasso）通过 L1 惩罚把不相关特征系数压到零。如果你做的是课程大作业，我的建议是两个都跑一遍，GridSearchCV 选超参数，顺便在报告里写对比。

from sklearn.model_selection import GridSearchCV from sklearn.linear_model import Ridge, Lasso param_grid = {"ridge__alpha": [0.01, 0.1, 1.0, 10.0, 100.0]} pipe_ridge = Pipeline([ ("scale", StandardScaler()), ("ridge", Ridge()) ]) grid_ridge = GridSearchCV(pipe_ridge, param_grid, cv=5, scoring="r2") grid_ridge.fit(X_train, y_train) print("best alpha:", grid_ridge.best_params_["ridge__alpha"]) print("best CV R2:", grid_ridge.best_score_) param_grid_lasso = {"lasso__alpha": [0.01, 0.1, 1.0, 10.0, 100.0]} pipe_lasso = Pipeline([ ("scale", StandardScaler()), ("lasso", Lasso()) ]) grid_lasso = GridSearchCV(pipe_lasso, param_grid_lasso, cv=5, scoring="r2") grid_lasso.fit(X_train, y_train) print("best alpha:", grid_lasso.best_params_["lasso__alpha"]) print("best CV R2:", grid_lasso.best_score_)

alpha 是惩罚强度：alpha 越小，模型越接近普通最小二乘，容易保留噪声；alpha 越大，系数被压得越狠，模型整体偏向均值，R2 下降。我这里给的网格是 0.01 到 100，覆盖了从接近最小二乘到强正则化的范围。需要说明的是，默认的 cv=5 会把数据随机打乱再分五份，这对时间序列并不是理想做法，所以上面只用于选超参数，最终评估仍然用按时间切出的 test 集。交叉验证的 R2 只能说明相对好坏，不代表最终测试表现。

3.3 系数解读：答辩时最值钱的一段输出

线性回归相比树模型最大的优势，是它把贡献摊在系数上。模型训练完，把系数按特征打印出来，你能直接看到 PM2.5 的预测主要被谁推高、被谁压低：

best_ridge = grid_ridge.best_estimator_ coefs = best_ridge.named_steps["ridge"].coef_ for name, coef in zip(features, coefs): print(f"{name:>16}: {coef:+.4f}")

注意这里的特征都经过了标准化，系数反映的是“该特征每变化一个标准差时，PM2.5 预测值变化多少个 µg/m³”，所以系数之间可以直接比较大小。常见的结果是：pm2_5_lag24 的系数通常在 0.5～0.8 之间，是最大的贡献者；温度系数为负，说明冷天 PM2.5 更容易堆积；风速系数为负，表示大气扩散条件好时浓度下降。如果你打印出来发现滞后特征系数是负的，或者 month 系数大到离谱，那多半是第 2 章的缺失值处理或者切分方式出了岔子，回头检查比继续调参更有效。

4. 模型评估与诊断：R2、RMSE、残差图怎么看才算数

4.1 评估指标：R2 之外，RMSE 才是能向老师解释的白话

单看 R2 很容易被误导。R2 是相对指标，测试集样本分布如果比较集中，R2 会虚高；如果出现了几次极端污染事件，R2 又会剧烈下跌。PM2.5 浓度是有单位的，它的预测误差应该用同单位的 RMSE（均方根误差）和 MAE（平均绝对误差）来说话。RMSE 对大误差施加平方惩罚，更适合空气质量这种“高估低估都影响决策”的场景。一行代码同时算三个：

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score y_true = test["pm2_5"] y_pred = best_ridge.predict(test[features]) rmse = mean_squared_error(y_true, y_pred, squared=False) mae = mean_absolute_error(y_true, y_pred) r2 = r2_score(y_true, y_pred) print(f"RMSE: {rmse:.2f} µg/m³") print(f"MAE: {mae:.2f} µg/m³") print(f"R2: {r2:.3f}")

补充一个版本注意点：sklearn 1.4 以上把 mean_squared_error 的 squared=False 标记为废弃，控制台会刷警告，建议改用from sklearn.metrics import root_mean_squared_error，一行拿到 RMSE。MAE 更稳健但不会惩罚极端误差，RMSE 对极端污染事件更敏感，两个指标一对比就能判断模型到底是被普通天数拖累，还是被少数重污染天拉垮。

一次重污染天的极端误差对 RMSE 的影响非常大：预测 300 而实际 150，单个样本就贡献了 150 的误差，平方后相当于几十个正常误差叠加。所以 RMSE 高通常意味着模型在极端情况上反应滞后，这比 R2 更能指出改进方向。

4.2 残差图：线性模型在一张图里现原形

残差是真实值减去预测值。把残差和预测值画成散点图，是最快、最省事、不需要任何统计知识的模型体检方式，它比任何玄学调参都直观。这种做法本质上就是一种“机器学习检测”——检测模型有没有系统性偏差：

import matplotlib.pyplot as plt residual = y_true.values - y_pred plt.figure(figsize=(8, 5)) plt.scatter(y_pred, residual, s=10, alpha=0.5) plt.axhline(0, color="red", linestyle="--") plt.xlabel("predicted PM2.5 (µg/m³)") plt.ylabel("residual (µg/m³)") plt.show()

正常的线性拟合残差应该围绕 0 均匀分布，没有明显形状。PM2.5 数据上典型的残差图会长这样：预测值小时残差集中在 ±20 以内，预测值越大，残差往上下两个方向散得越开，形成喇叭形。这说明数据存在异方差性——高浓度时段模型的不确定性天然更大。这不是 bug，是 PM2.5 本身受实时排放、局部燃烧源影响的噪声数据特征。写报告时你甚至可以专门用一张图证明“线性回归的局限”，这比强行解释高 R2 更有说服力。

4.3 如果测试 R2 只有 0.4，是模型失败了吗

不是。线性回归预测 PM2.5，在只给气象和前值特征的情况下，R2 落在 0.3～0.6 是常态。污染物浓度不是纯线性过程：光化学反应、二次气溶胶生成、突发沙尘事件都是强非线性因素。所以线性模型的“低分”反而是合理结果。它真正的问题不在于 R2 低，而在于残差图如果出现明显的弯曲弧线，说明“线性关系”这个前提本身不对，这时候才考虑对目标做对数变换、加入非线性特征，或者直接换树模型。

4.4 评估结论怎么写：三句话撑起课程报告

课程报告不需要堆数字，但需要把数字读出来。我一般按三句话组织评估结论：第一句给水平，例如“测试集 RMSE 为 38 µg/m³，R2 为 0.42，模型能解释四成左右的浓度变化”；第二句指出失真样本，例如“残差图显示高浓度段预测偏差明显扩大，模型对重污染日反应不足”；第三句给下一步，例如“加入 48 小时滞后特征和风速突变标志后，预计能改善拐点滞后”。这样写，答辩老师听到的是一个完整的建模闭环，而不是一个孤立的分数。

5. 避坑与排查：PM2.5 预测中最容易翻车的五个现场

这一章把我在大作业里几乎必然会踩的坑按“现象-原因-解决”整理出来。每一条都不是网上搜不到的常识，是实际跑数据时才见得着的鬼。

5.1 测试 R2 虚高到 0.8 以上：先查切分方式

现象：模型训练完，测试集 R2 高达 0.8 以上，可画出预测曲线后发现“拟合”得不像话，连拐点都提前对齐了。原因：使用了 train_test_split 随机切分，未来时间段样本混进了训练集，模型间接通过同一天的样本复制了答案。这在大作业源码里非常常见，因为随机切分是多数人的默认动作。解决：改成按时间顺序切分，先时间后标签；交叉验证改用 TimeSeriesSplit。检查方法：打印 test 里第一条样本的时间，如果它不晚于训练集最后一条样本的时间，就是典型泄漏。

5.2 补完缺失后样本量从 3 万变成 1 万：别无脑 dropna

现象：构造滞后和滚动特征后，dropna 把数据行删掉近半，模型 R2 下降，系数符号也变怪。原因：shift(24) 和 rolling(168) 本身会产生大量空值，早期样本全被清空；如果原始缺失又没处理，dropna 会连锁带走大片行。解决：先按 2.2 节补齐原始缺失，再构造特征，最后只 drop 由特征窗口引入的缺失；如果样本仍然紧张，可以把滚动窗口从 7 天缩短为 3 天，减少边界损失。这个坑的特征是“代码能跑、数据行数神秘变少”，所以每次 dropna 前先记一下 df.shape。

5.3 alpha 怎么调都无效：特征缩放的账没算

现象：对 Ridge 做 GridSearchCV，alpha 从 0.01 试到 100，交叉验证分数纹丝不动，或者全取边界值。原因：特征没做标准化，co、wind_speed、month 数量级完全不同，L2 惩罚实际只作用在大尺度特征上，名义上是均匀惩罚，实际变成了选择性的随机惩罚。解决：把 StandardScaler 放进 Pipeline，确保 alpha 作用在无量纲特征上；打印标准化后的特征标准差，确认都接近 1 再调参。如果你手动建模而不是用 Pipeline，记得用同一个 scaler 变换测试集，不能重新 fit。

5.4 预测结果全是同一个数值：目标列被覆盖了

现象：模型能正常训练和预测，但 y_pred 全是一个常量，比如 55.7。原因：最常见的是在构造滞后特征时误把目标列覆盖了，典型代码是df["pm2_5"] = df["pm2_5"].shift(24)，把目标列整体往后移了一格，于是所有真实目标都指向下一行，模型只能学到一个均值。解决：打印 y_train.head(10) 和 y_train.mean()，对比原始数据列；用 df.columns 确认没有列名覆盖；再检查拟合后 intercept 是否近似等于 y_train.mean()。如果一个模型输出恒为均值，说明它完全没有学到任何特征信息，数据管道一定在某个环节换了列。

5.5 GridSearchCV 的风度分比时间切分还高：KFold 不适合时序

现象：交叉验证时 R2 高达 0.7，换到按时间切的 test 集却只有 0.4，两个分数打架。原因：GridSearchCV 默认用 KFold 随机打乱数据，前面说过，这等于把未来样本泄漏进训练过程；交叉验证分数虚高，真正的时间外评估才暴露水平。解决：交叉验证用 TimeSeriesSplit，或者在选完超参数后，完全不看交叉验证分数，只看时间切分测试集的结果。给你个血泪经验：任何时间序列模型报告“五折交叉验证 R2 0.8”，如果没说清楚切分方式，都值得先怀疑一下。

排查顺序建议：发现分数不合理时，先查切分方式，再查目标列有没有被覆盖，然后查特征缩放，最后才怀疑模型本身。按这个顺序走，绝大多数“怎么调都不对”的问题在十分钟内能定位。

6. 再往前走一步：用滞后特征、交叉验证和曲线可视化把报告撑住

6.1 把滞后窗口拉开，模型会更稳

在 2.3 节的基础上，再补两天前和一周前的特征：

df["pm2_5_lag48"] = df["pm2_5"].shift(48) df["pm2_5_lag168"] = df["pm2_5"].shift(168) df = df.dropna()

特征变多后，记得把 features 列表更新，再走一遍第 3 章的 Pipeline。线性回归吃不了几百个特征，但这类手选滞后特征完全在它能力范围内。加入一周前的特征后，模型在冬季采暖期的预测通常会更稳。如果发现 lag168 的系数接近 0，那是 Lasso 在帮你做特征选择，可以直接删掉，不是 bug。

6.2 交叉验证换成 TimeSeriesSplit，评估结果更有说服力

课程报告里写“我用五折交叉验证”很常见，但时间序列任务要接一句“用的是时间序列切分”。TimeSeriesSplit 是 sklearn 内置的时序交叉验证，每折都保持时间顺序：

from sklearn.model_selection import TimeSeriesSplit import numpy as np tscv = TimeSeriesSplit(n_splits=5) cv_scores = [] for train_idx, val_idx in tscv.split(X_train): X_cv_train, X_cv_val = X_train.iloc[train_idx], X_train.iloc[val_idx] y_cv_train = y_train.iloc[train_idx] y_cv_val = y_train.iloc[val_idx] pipe_lr.fit(X_cv_train, y_cv_train) cv_scores.append(pipe_lr.score(X_cv_val, y_cv_val)) print("CV R2 mean: %.3f ± std: %.3f" % (np.mean(cv_scores), np.std(cv_scores)))

每折训练集都在验证集之前，分数波动能反映模型在不同季节的稳定性，比单次时间切分更可信。

6.3 预测曲线叠上去，看模型在哪里反应慢

数字指标再全，不如一张曲线直观：

plt.figure(figsize=(10, 4)) plt.plot(test.index[:200], y_true.values[:200], label="true") plt.plot(test.index[:200], y_pred[:200], label="predicted") plt.legend() plt.xticks(rotation=45) plt.ylabel("PM2.5 (µg/m³)") plt.show()

截取测试集前 200 个小时，你会看到模型在普通日子贴合得不错，但在浓度骤升或骤降的拐点会滞后半拍到一天。这是线性回归的固有惯性，也说明下一步改进要么加入气象骤变特征，要么换树模型。这些观察写进报告，比直接贴一个 R2 更有说服力。

我最初做这个题时，第一次跑出 0.85 的 R2 高兴了半天，后来改用时间顺序切分，分数掉到 0.4，才知道之前是窃取未来的假象。从那以后，我每次拿到时间序列数据，第一件事就是看训练集最后一行和测试集第一行之间有没有时间空隙。这篇笔记里写的流程，从字段审查到残差图、时序交叉验证，都是我延续至今的习惯，希望帮到你。

本文还有配套的精品资源，点击获取