线性回归实现PM2.5预测:从数据清洗到模型部署
简介:这份压缩包是基于线性回归的PM2.5预测系统Python源码,面向机器学习初学者、数据竞赛参与者及课程设计开发者,提供一套完整可运行的空气质量预测项目方案。包内共19个文件,包含12个csv数据文件(训练集、测试集及中间处理结果)、1个py主程序、1个npy模型权重文件以及5张过程截图,压缩包整体仅2.14MB,便于快速下载与本地复现。代码涵盖数据读取、预处理、矩阵运算、线性回归模型训练与预测等关键环节,对理解台湾地区空气质量监测数据的清洗方式与特征拼接有一定帮助。目前已有593人学习下载,适合希望以真实数据集练习线性回归、参考模型保存与结果输出流程的读者。资源附带的图片可辅助查看数据分布和预测效果,整体结构简洁,能节省从零搭建环境的时间。
1. 线性回归做 PM2.5 预测:为什么这个老算法依然值得跑通一遍
拿到一份名为“基于线性回归的PM2.5预测系统python源码.rar”的压缩包,你大概率带着两个问题:这套系统到底能不能跑出能看的预测结果,以及线性回归这种最基础的机器学习算法,凭什么被单独打包成一个项目。先说结论:PM2.5 预测的本质不是“预测明天是 50 还是 80”,而是从历史观测数据里找出哪些气象条件和污染物浓度对未来的 PM2.5 有稳定影响。线性回归在这个任务上不是一个炫技的选择,而是一个性价比极高的起点——数据量不到几千条也能训练,特征列基本不用做复杂变换,训练时间以秒计,而且每个特征的系数可以直接解释成“风速每提高 1m/s,PM2.5 大约下降多少微克”。这正好是新手从传统统计分析跨进机器学习的第一步,也是很多量化环境数据分析项目的底层基线。
这个方案适合三类人:第一类是刚接触 Python 机器学习、想用真实数据跑通全流程的初学者;第二类是环境科学背景、需要快速验证气象与污染物相关性的研究人员;第三类是想给现有监测系统加一个“简单可用”预测模块的开发者。你不需要 GPU,不需要深度学习的黑匣子,一个带 pandas 和 scikit-learn 的 Python 环境就够了。下面我从数据准备开始,把这套系统的完整落地路径、参数选择和踩坑记录都拆给你看。
2. 数据准备与特征构造:决定线性回归上限的是表,不是算法
2.1 拿到历史监测数据后,第一步不是训练而是清洗
任何 PM2.5 预测项目的第一关都是数据。常见做法是使用城市空气质量监测站的历史数据,字段通常包括日期时间、PM2.5 浓度、PM10、SO2、NO2、CO、O3,以及温度、气压、风速、风向、降雨量等气象变量。很多公开数据集是按小时采样的,一天 24 条,一个月大约 720 条,一年的数据足够支撑一个像样的线性回归模型。
清洗阶段最常见的需求是处理缺失值。监测设备断电、校准、通讯故障都会产生空行。线性回归对缺失值很敏感,因为训练时每一行特征必须完整。我的习惯是先看缺失比例再决定策略:如果单列缺失低于 5%,用中位数填充即可;如果超过 20%,那这一列基本是废的,直接删掉。注意不要用均值填充,因为 PM2.5 和气象数据往往有偏态分布,均值会被极端值拉高,中位数更稳健。另外,时间字段要醒目地解析成 datetime 类型,否则后面做滞后特征时会反复踩坑。
import pandas as pd
df = pd.read_csv("pm25_data.csv", parse_dates=["datetime"])
print("原始数据形状:", df.shape)
print("缺失值统计:\n", df.isnull().sum())
# 删除缺失比例过高的列
df = df.dropna(thresh=len(df) * 0.8, axis=1)
# 数值列用中位数填充
num_cols = df.select_dtypes(include=["float64", "int64"]).columns
df[num_cols] = df[num_cols].fillna(df[num_cols].median())
# 检查时间是否连续,是否存在重复时间戳
print("重复时间戳数量:", df["datetime"].duplicated().sum())
这段代码做了三件事:先看数据形状和缺失分布,再按缺失阈值砍列,最后用中位数填充保证每一行都能进入训练。参数上 thresh=len(df) * 0.8 意味着某列必须有至少 80% 的非空值才保留,对于小时级空气质量数据,这个阈值通常比 0.8 更严一点更稳妥,我一般建议提到 0.85。重复时间戳必须处理,否则滞后特征会串行错乱。
2.2 线性回归不是只能学“当前值”,滞后特征才是预测的关键
如果你直接拿当前时刻的风速、温度去预测同一时刻的 PM2.5,那做的是“拟合”不是“预测”。预测系统的核心在于用已知信息推断未来。因此要构造滞后特征(lag features):比如用 t-1 小时、t-2 小时、t-24 小时的 PM2.5 浓度和气象数据,来预测 t+1 小时的 PM2.5。为什么滞后特征有效?因为污染物浓度有很强的自相关性——今天晚上的 PM2.5 数值和昨晚同时段相关性很高,而且气象变化对污染物扩散的影响不是瞬时的,往往滞后几小时才体现出来。
def make_lag_features(df, target_col="PM2.5", lags=[1, 2, 3, 24]):
df = df.sort_values("datetime").reset_index(drop=True)
for lag in lags:
df[f"{target_col}_lag_{lag}"] = df[target_col].shift(lag)
df[f"wind_speed_lag_{lag}"] = df["wind_speed"].shift(lag)
df[f"temperature_lag_{lag}"] = df["temperature"].shift(lag)
df = df.dropna().reset_index(drop=True)
return df
df_feat = make_lag_features(df)
print("构造特征后的数据形状:", df_feat.shape)
print(df_feat.head(10)[["datetime", "PM2.5", "PM2.5_lag_1", "PM2.5_lag_24"]].to_string())
shift 是 pandas 里做滞后特征的利器。 shift(1) 表示取上一时刻的值, shift(24) 在小时数据里表示昨天同一时刻,天然捕获了日周期。注意 dropna 的位置:生成滞后特征后,最前面的 24 行会因为缺少历史值变成 NaN,直接丢掉。还有一个细节:必须按时间排序后再 shift,否则滞后关系全乱了。如果你的数据是小时级但有间断,建议先重采样到小时网格再做 shift,保证时间索引对齐。
除了滞后特征,时间编码也是线性回归需要的。把“小时”“星期几”“月份”直接当整数塞给线性回归不太合适,因为线性回归认为 23 和 0 之差是 23,但 23 点与 0 点在周期性上是相邻的。常见做法是把小时拆成两个连续特征: sin(2*pi*hour/24) 和 cos(2*pi*hour/24) 。同理,月份和星期也可以用这种方式。这样模型能学到“夜间 vs 白天”“工作日 vs 周末”的平滑差异,而不会强行把 23 点和 0 点切到两端。
import numpy as np
hour = df_feat["datetime"].dt.hour.values
month = df_feat["datetime"].dt.month.values
weekday = df_feat["datetime"].dt.weekday.values
df_feat["hour_sin"] = np.sin(2 * np.pi * hour / 24)
df_feat["hour_cos"] = np.cos(2 * np.pi * hour / 24)
df_feat["month_sin"] = np.sin(2 * np.pi * month / 12)
df_feat["month_cos"] = np.cos(2 * np.pi * month / 12)
这是把周期性时间变量转换为线性回归可解释形式的常用手段,和热词里反复出现的“线性回归模型”“python 入门教程”几乎绑定出现。很多初学者在这里会漏掉 np.pi 的浮点精度问题,导致计算被截断成整数,但 pandas 的时间字段默认是 int64,而 hour / 24 在 Python 3 里是浮点除法,所以一般不出错。如果你的 Python 环境是 2.7,那就得手动转 float。现在的机器学习项目基本都跑在 3.8 以上,这个问题遇到的概率已经很小了。
3. 用 scikit-learn 训练多元线性回归:从特征矩阵到模型落盘
3.1 划分训练集和测试集必须按时间切,不能随机打乱
线性回归训练本身不复杂,但数据划分的方式直接决定模型评估是否可信。很多教程习惯用 train_test_split 的默认随机切分,这在 PM2.5 预测场景是禁忌。因为时间序列数据有着强烈的时间自相关性,随机打乱会把未来的数据混进训练集,模型相当于“偷看答案”,在测试集上表现虚高,部署到真实环境立刻失效。正确做法是按时间顺序切分:比如用前 80% 的时间段做训练,后 20% 做验证。更进一步,应该留出最后一周作为最终测试集,因为短期预测的验证方式模拟了“用昨天预测今天”的真实过程。
from sklearn.model_selection import train_test_split
feature_cols = [c for c in df_feat.columns if c not in ["datetime", "PM2.5"]]
X = df_feat[feature_cols].values
y = df_feat["PM2.5"].values
cutoff = int(len(X) * 0.8)
X_train, X_test = X[:cutoff], X[cutoff:]
y_train, y_test = y[:cutoff], y[cutoff:]
print(f"训练集大小: {X_train.shape}, 测试集大小: {X_test.shape}")
cutoff 是时间序列切分的关键参数。如果你有 8760 条小时数据(一年), int(8760 * 0.8) = 7008 ,训练集是前 7008 条,测试集是最后 1752 条,对应完整的一个冬季或春季,能真实反映模型在陌生时段的预测能力。注意不要把 train_test_split 默认的 shuffle=True 用在这里,否则时间顺序就毁了。如果坚持要用 sklearn 的函数,必须显式传入 shuffle=False 。
3.2 标准化的必要性和 model 训练
线性回归对特征尺度比较敏感。如果风速范围是 [0, 20],而 PM2.5 的滞后值范围是 [0, 500],那么大数值的特征会在损失函数中占据更高权重,导致梯度更新偏向该特征。用 StandardScaler 做标准化后,所有特征均值归 0,方差归 1,模型才能公平学习。注意一个常见错误:标准化器必须只用训练集拟合,然后用同一个 scaler 转换测试集,不能对全部数据一起 fit。否则测试集信息泄露到训练过程,评估还是不准。
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
model = LinearRegression()
model.fit(X_train_scaled, y_train)
print("模型截距:", model.intercept_)
print("特征系数:", dict(zip(feature_cols, model.coef_)))
fit_transform 用于训练集,transform 用于测试集,这是 sklearn 的标准契约。 model.coef_ 保存每个特征的权重,权重绝对值越大表示该特征对预测结果影响越强。比如 PM2.5_lag_1 的系数通常是所有特征里最大的,因为 t-1 时刻的浓度对 t+1 时刻有着极强的正相关。查看系数是理解线性回归可解释性的最佳途径,这也是线性回归区别于随机森林或 XGBoost 的地方——每个系数的业务含义都可以直接读出来。在执行这段代码之前,建议先打印 X_train_scaled.mean() 和 X_train_scaled.std() ,确认标准化后的均值接近 0、标准差接近 1,这是排查特征处理是否正确的快速方法。
3.3 训练完成后的保存与预测调用
训练好的模型和标准化器都需要保存,否则每次预测都要重新训练一遍。常见做法是用 joblib.dump 把 model 和 scaler 分别落盘成 .pkl 文件。预测阶段只需要加载这两个文件,把新实时数据按相同顺序构造特征并标准化,然后调用 model.predict 。注意特征顺序必须与训练时完全一致,所以最好把 feature_cols 也保存下来,加载后检查列名对齐。
import joblib
joblib.dump(model, "pm25_model.pkl")
joblib.dump(scaler, "scaler.pkl")
joblib.dump(feature_cols, "feature_cols.pkl")
# 预测新数据时
loaded_model = joblib.load("pm25_model.pkl")
loaded_scaler = joblib.load("scaler.pkl")
loaded_features = joblib.load("feature_cols.pkl")
# 假设 new_df 是经过同样滞后处理的实时数据
new_X = new_df[loaded_features].values
new_X_scaled = loaded_scaler.transform(new_X)
pred = loaded_model.predict(new_X_scaled)
print(f"预测 PM2.5 浓度: {pred[0]:.2f} ug/m3")
用 joblib 保存 sklearn 模型是标准做法,比 pickle 更高效,尤其是模型包含大数组时。这里要特别注意:scaler 也要保存并一起加载,很多项目只保存模型不保存 scaler,导致在线预测时数据没被标准化直接塞进模型,输出数值完全不可信。另一个坑是特征顺序:如果在线特征拼接顺序与训练时不同,预测结果会错乱,所以保存 feature_cols 并在加载后和实际数据的列顺序做一次断言校验。
4. 模型评估与调参:R2、MAE 与特征筛选的真实含义
4.1 回归任务要看的四个指标,以及它们的业务边界
线性回归的评估指标不能只看决定系数 R2。在 PM2.5 预测场景中,R2 超过 0.9 并不罕见,尤其是滞后特征把 t-1 时刻的浓度直接带入模型后,R2 会显得很高。但 R2 高不代表预测实用,因为你需要知道误差到底有多大。常用的互补指标是 MAE(平均绝对误差,单位微克/立方米)和 RMSE(均方根误差)。MAE 更像“普通情况下的预期偏差”,RMSE 会放大异常点的影响。如果 RMSE 明显大于 MAE,说明模型在大浓度峰值时段表现较差——这正是重污染过程,而重污染恰恰是预测系统最应该抓准的场景。
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
y_pred = model.predict(X_test_scaled)
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(f"MAE: {mae:.2f} ug/m3")
print(f"RMSE: {rmse:.2f} ug/m3")
print(f"R2: {r2:.3f}")
写评估代码很简单,但读指标需要业务常识。如果 MAE 在 15-25 微克/立方米左右,对小时级 PM2.5 预测来说是一个可以接受的区间。空气质量指数(AQI)的分段大致是 0-35 为优,35-75 为良。如果你的 MAE 已经超过 35,那预测结果连“良/优”这种粗粒度区分都做不准,这时不应继续调线性回归参数,而是要回到特征工程层面——增加滞后阶数、引入风向类别编码或加入降雨量。这就是我要强调的:R2 只是拟合优度,MAE 和 RMSE 才是业务决策依据。
4.2 回归系数矩阵与特征相关性排查
做完一轮训练后,把特征列名和系数打印成表格,可以直观发现哪些特征是无效的,哪些特征是共线性的。共线性是线性回归的一个典型问题:当两个特征高度相关时,模型会把权重在这两个特征之间随意分配,导致单个系数的解释性崩塌。PM2.5 数据里典型的共线性来源是 PM10 和 PM2.5 本身强相关,以及 PM2.5_lag_1 与 PM2.5_lag_2 之间接近 0.95 的相关系数。
corr_matrix = df_feat[feature_cols].corr()
high_corr_pairs = []
for i in range(len(corr_matrix.columns)):
for j in range(i):
if abs(corr_matrix.iloc[i, j]) > 0.85:
col_i = corr_matrix.columns[i]
col_j = corr_matrix.columns[j]
high_corr_pairs.append((col_i, col_j, round(corr_matrix.iloc[i, j], 3)))
print("相关系数大于0.85的特征对:")
for pair in high_corr_pairs:
print(pair)
这是一个快速排查共线性的方法。阈值 0.85 是经验值,超过这个值通常需要考虑删掉其中一个或做 PCA。删除策略优先保留业务含义更直接的特征。比如保留 PM2.5_lag_1 而删掉 PM2.5_lag_2,因为滞后 1 小时的业务意义更清晰;PM10 和 PM2.5 同时出现时,如果目标是预测 PM2.5,删除 PM10 反而让系数更稳定。如果删完特征后测试集 MAE 没明显变差,那就大胆删。线性回归的优势就在这里:它给你一个统计学工具,用简单计算换取模型的可解释性和稳定性。
4.3 用岭回归替代普通最小二乘:一个参数解决过拟合
线性回归的普通最小二乘解在某些情况下会震荡——当训练集存在近似共线性时,系数会变大甚至变号,导致模型在测试集上不稳定。解决办法是加 L2 正则化,也就是岭回归(Ridge Regression)。它的损失函数在均方误差的基础上加了一项惩罚项,惩罚力度由 alpha 控制。alpha 越大,系数越被压缩向 0,模型越平滑。在 PM2.5 预测中,岭回归通常比普通线性回归更能抵抗滞后特征之间的共线性。
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
alphas = [0.1, 1.0, 10.0, 100.0]
param_grid = {"alpha": alphas}
ridge = Ridge()
grid = GridSearchCV(ridge, param_grid, cv=5, scoring="neg_mean_absolute_error", n_jobs=-1)
grid.fit(X_train_scaled, y_train)
best_alpha = grid.best_params_["alpha"]
print("最优 alpha:", best_alpha)
best_ridge = grid.best_estimator_
y_pred_ridge = best_ridge.predict(X_test_scaled)
print(f"岭回归 MAE: {mean_absolute_error(y_test, y_pred_ridge):.2f} ug/m3")
GridSearchCV 是 sklearn 里做超参数搜索的标准工具。这里用 5 折交叉验证搜索最合适的 alpha,评估指标是负 MAE——sklearn 里所有“越大越好”的评分指标,如果是误差类指标就要取负值。使用交叉验证而非测试集来选参,能避免把测试集信息泄漏进超参数选择过程。注意,这里交叉验证的折切分也在时间序列上存在泄漏问题,但因为在训练集内部,风险可控,实际项目更严格会使用 TimeSeriesSplit 。如果你想在时间和效果之间取平衡,先跑默认 KFold 看趋势,再用 TimeSeriesSplit 做最终确认也说得过去。
5. 避坑指南:PM2.5 线性回归项目最常见的五个翻车现场
5.1 现象:预测结果全部趋近于一个常数,且这个常数约等于 y 的均值
原因:滞后特征没有生效,或者特征列中净剩的只有时间编码,而历史 PM2.5 值没有进模型。线性回归学不到复杂非线性关系,如果输入特征和目标之间没有强关联,它只能输出均值作为最“安全”的预测。特征矩阵如果全部是常数或接近常数,模型学出来 intercept 等于目标均值,coef 全为 0。这通常发生在 shift 之后忘记 dropna,导致所有滞后特征全为 NaN,pandas 会自动把这些行丢掉,剩下可用的样本量极小。
解决:检查 df_feat.shape 与原始 df 的差异,确认滞后特征确实生成了非空值。打印 df_feat[["PM2.5_lag_1"]] 的前几行,非空值数量应接近总行数减去最大滞后阶数。如果特征量很少,建议无条件把 PM2.5_lag_1 和 PM2.5_lag_2 作为必选特征加入。
5.2 现象:训练集 R2 是 0.98,测试集 R2 掉到 0.3
原因:随机切分了时间序列数据,或者滞后特征的阶数取得太大(如 lag 48、72),模型大量记忆了训练集中短期波动的细节,而测试集是完全不同的时间段。这是时间序列预测里最常见的“轨道平移”翻车。
解决:检查是否用了默认 train_test_split 。如果是,改成按 cutoff 顺序切分,训练集只保留起点到 cutoff 之间的样本。如果已经是顺序切分但测试集 R2 仍然极低,检查训练集和测试集的时间跨度——比如训练集是春夏,测试集是冬季,PM2.5 浓度在冬季会因为供暖而系统性抬高,线性回归很难外推到没见过的量级。这时应把整个数据集按“包含完整四季”的方式划分,而不是简单按 80/20 切。
5.3 现象:模型系数中有极大正值和极大负值,符号与业务常识相反
原因:共线性。典型如 PM10 和 PM2.5 同时作为特征,模型难以区分二者贡献,系数在两者之间分摊并出现符号翻转。业务上 PM10 升高必然伴随 PM2.5 升高,但系数可能显示 PM10 升、PM2.5 降,这完全失真。
解决:计算相关系数矩阵,把绝对值高于 0.85 的特征对找出来。保留业务意义更直接的 PM2.5 滞后值,删除 PM10。如果删掉后 RMSE 反而小幅上升,这是正常现象——模型牺牲了一点点精度换取系数可解释性,对这个项目来说是值得的。
5.4 现象:测试集 MAE 尚可,但重污染时段(PM2.5 > 200)的预测严重偏低
原因:线性回归的损失函数是均方误差,模型会优先减少常见浓度区间的误差,对于少数极端高值,即使预测偏差较大,在总平方误差中的占比也不足以让模型调整到准确捕捉峰值。这在重污染预测中是天然弱势。
解决:一方面可以引入特征来提升峰值识别,比如增加污染过程持续时间特征(连续几小时 PM2.5 超过 150 的累计值),另一方面可以换用分位数回归或梯度提升树,但那就偏离了线性回归的主题。线性回归能做得更实际的是:把预测结果按阈值分层,输出时对高值区间加一个经验修正系数。这不是严谨做法,但工程上可以先顶一阵,后续再升级模型。
5.5 现象:在线预测时加载模型返回的结果与离线测试完全对不上
原因:在线预测时没有使用同一个 scaler,或者特征列顺序不一致,或者新数据的时间字段未对齐导致滞后特征计算错位。最常见的是只加载了模型文件,而忘了加载 scaler。
解决:落盘时统一保存三件套——model、scaler、feature_cols,加载后先断言 new_df.columns.tolist() == loaded_features ,再做 transform。把在线预测逻辑写成函数,每次调用前自检特征顺序,这样即使后面有人改了训练脚本也不会悄悄破坏线上逻辑。
6. 进阶技巧:把线性回归预测结果用残差分析再校准一轮
简单线性回归落地后,多数人会直接上线,但我想多说一个见效快的技巧:残差分析。所谓残差,就是真实值减预测值的差。把残差当成新的分析对象,按时间、按 PM2.5 浓度区间、按风速区间去统计残差分布。比如,你可能会发现残差在风速低于 2m/s 时普遍为正(即预测偏低、真实高),在风速高于 5m/s 时普遍为负。这说明模型没有充分捕捉“静稳天气下的污染物累积效应”。针对这个现象,可以构造一个交互特征:风速与 PM2.5_lag_1 的乘积,让模型学习“高污染且低风速时预测应该更高”。
df_feat["wind_pm25_interaction"] = df_feat["wind_speed"] * df_feat["PM2.5_lag_1"]
# 重新训练
feature_cols.append("wind_pm25_interaction")
X = df_feat[feature_cols].values
X_train, X_test = X[:cutoff], X[cutoff:]
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)
model.fit(X_train_scaled, y_train)
y_pred_new = model.predict(X_test_scaled)
print(f"加入交互特征后 MAE: {mean_absolute_error(y_test, y_pred_new):.2f} ug/m3")
交互特征就是两个原始特征相乘,它给线性回归引入了“协同效应”的建模能力。如果你发现加入后 MAE 下降了 3 微克以上,就值得保留;如果只降了 0.5,则可能是噪声掺入,建议删掉。另外一个更实用的落地技巧是输出概率区间:线性回归给的是点预测,但你可以用历史残差的标准差构造一个近似置信区间,预测值加减 1.65 倍标准差就是 90% 的波动范围。这对环境决策者是有用的——他们不只是想知道“明天 PM2.5 是多少”,更想知道“最坏可能到多少”。这两步做下来,你的线性回归预测系统已经不是玩具,而是一个具备在线预测、误差解释、置信区间输出能力的完整模块,作为毕业设计、量化环境分析项目、或企业环境监测的可解释基线,都足以交差。我自己在这个方向上迭代过三轮,最深的感觉是:线性回归不是被替代的老古董,它是让你把数据和业务之间关系彻底想清楚的那面镜子。希望帮到你。
更多推荐
所有评论(0)