线性回归预测PM2.5：从特征工程到模型部署的完整指南

简介面向机器学习入门者与数据分析学习者基于线性回归的PM2.5预测系统完整源码包围绕空气污染物浓度预测任务提供从数据读取、预处理到模型训练与结果导出的全流程Python实现。压缩包共19个文件以12个csv数据文件为主涵盖训练集、测试集及预测输出另含5个png结果可视化图、1个npy模型参数文件及1个主程序py脚本整体大小2.14MB结构清晰便于直接运行与二次开发。资源中代码对繁体中文big5编码的csv进行了正确处理并演示了矩阵运算、最小二乘法参数求解、特征拼接等关键步骤适合课程设计、毕业设计或竞赛练手。目前已有593人浏览学习可作为理解线性回归实际应用的直观参考帮助读者快速复现实验并对比预测效果。1. 用线性回归预测 PM2.5先分清能预测和预测得准PM2.5 浓度预测在环境监测和健康预警场景里一直有真实需求而它恰恰是线性回归算法教学案例里最贴地的一个。PM2.5 与气象变量温度、气压、湿度、风速之间存在可解释的线性相关趋势用线性回归可以在几分钟内得到一条可用的预测基线。这不是玄学而是能把特征选择、数据切分、残差诊断、正则化这些环节全部落在真实数据上验证。这里要展开的是一套 Python 源码级方案以公开的北京 PM2.5 数据集为输入用 scikit-learn 的 LinearRegression 建模把评估、调优、模型持久化串成可离线运行的预测系统。它适合正在做课程设计、刚入坑机器学习回归任务或者想把预测系统从概念落成可运行代码的工程师。2. 数据先行PM2.5 预测的数据集选择与预处理2.1 公开数据集怎么找UCI 北京 PM2.5 数据集做 PM2.5 预测第一步不是写模型而是拿到一份字段完整的观测数据。常见做法是使用 UCI Machine Learning Repository 上的 Beijing PM2.5 Data Set它记录了 2010 年 1 月 1 日到 2014 年 12 月 31 日北京某监测站点逐小时的气象与污染物浓度。数据总量约 4.3 万条覆盖四年完整季节周期对线性回归建模来说这个数据量既不会让训练发散也不会小到过拟合不可控。字段含义类型在模型里的角色year / month / day / hour时间戳四要素int构造时间特征的基础pm2.5PM2.5 浓度μg/m³float目标变量 yDEWP露点温度℃float数值特征TEMP温度℃float数值特征PRES气压hPafloat数值特征cbwd组合风向str类别特征需要编码Iws累计风速m/sfloat数值特征Is / Ir积雪 / 降雨时长小时int数值特征稀疏这个数据集的价值在于它同时包含连续特征、类别特征和时间序列特性。其中 cbwd 是四类风向NE、NW、SE 和代表静风的 cvcv 类观测往往对应高浓度时段在建模时要保留这个原始分类不要合并。需要提醒的是UCI 版本早期记录里有 pm2.5 的缺失值这是真实环境监测数据的特点直接喂给 sklearn 会报错必须在预处理里消化掉。2.2 缺失值与时间戳处理别让 NaN 吃掉模型的 R²拿到 CSV 之后我一般先做三件事把年月日时拼成时间索引、统计缺失分布、查看风向类别的样本量。代码如下import pandas as pd import numpy as np df pd.read_csv(PRSA_data_2010.1.1-2014.12.31.csv) # 拼接标准时间索引便于排序和构造滞后特征 df[datetime] pd.to_datetime(df[[year, month, day, hour]]) df df.sort_values(datetime).reset_index(dropTrue) print(df.isnull().sum()) # 查看各列缺失数量 print(df[cbwd].value_counts()) # 查看风向类别分布pd.to_datetime 接收的数据帧必须包含 year、month、day、hour 四列输出一个可直接排序的 datetime 列sort_values 按时间重排让后续 shift 操作严格沿时间轴进行。isnull().sum() 能快速暴露缺失集中区如果 pm2.5 缺失行出现在数据开头直接用 dropna 删掉即可如果缺失出现在中间时段用前向填充 ffill 更符合物理逻辑——PM2.5 浓度属于缓变过程前一小时的浓度比全天均值更接近当前真实值。提示不要用 SimpleImputer 做全局均值填充处理时间序列。均值会把浓度曲线人为拉平制造出训练集里不存在的时间跳变最终削弱的恰恰是线性回归对污染过程的还原能力。对 cbwd 这类类别特征最简单的方案是 pandas 的 get_dummies 做独热编码。这里要注意 drop_first 参数四类风向独热后会生成 4 列但类别间存在共线性我习惯设置 drop_firstTrue只保留 3 个哑变量列避免特征矩阵里出现完全相关的列导致正规方程不可逆。2.3 特征工程滞后特征与周期变量的组织方式线性回归对特征的组织方式非常敏感。PM2.5 预测里最常用的先验是过去的浓度影响现在因此要提前把浓度序列的滞后项做进特征矩阵# 按小时粒度构造多个滞后项捕捉短时惯性与日周期 for lag in [1, 3, 6, 12, 24]: df[fpm2.5_lag{lag}] df[pm2.5].shift(lag) df df.dropna().reset_index(dropTrue)shift(lag) 把 pm2.5 整列下移 lag 行使第 t 行的滞后特征恰好等于 t-lag 时刻的真实浓度。选择 1、3、6、12、24 这几个步长对应的分别是小时级惯性、半日趋势和 24 小时日循环滞后 1 小时的系数通常会显著大于其他滞后项这说明 PM2.5 浓度具有很强的短时记忆。dropna 必须放在最后否则前 24 行因 shift 产生的 NaN 会进入训练集连带污染损失函数的计算。时间字段 hour 不能直接作为数值特征喂给模型。第 23 时与第 0 时数值上相差 23物理上只相隔 1 小时线性回归会把这种虚假的数值距离当成真实规律。改进做法是周期编码# 周期编码把小时映射到单位圆上的两点保留周期性 df[hour_sin] np.sin(2 * np.pi * df[hour] / 24) df[hour_cos] np.cos(2 * np.pi * df[hour] / 24) df[weekend] (df[datetime].dt.weekday 5).astype(int)sin/cos 编码把 hour 转换成单位圆上的两个坐标模型可以学到午夜到清晨浓度爬升、午后随扩散下降这类以 24 小时为周期的模式。hour_sin 和 hour_cos 需要同时保留只用其中一个会把早晚高峰的信息折叠到一起。weekend 是简单的星期特征因为周末的交通排放与工作日差异明显这个特征在包含移动源污染的模型里通常有正向贡献。到这里特征矩阵已经包含数值特征、周期特征、哑变量和滞后特征可以作为第 3 章建模的输入。3. 用 sklearn 跑通线性回归最小代码与必调参数3.1 训练集 / 测试集划分时序数据不能随机打乱PM2.5 数据是典型的时间序列。如果调用 train_test_split 默认的随机切分训练集和测试集会混入同一个时间段样本模型相当于提前看到了未来测试集上的 R² 虚高但没有部署参考价值。正确做法是按时间顺序切分feature_cols [DEWP, TEMP, PRES, Iws, Is, Ir, pm2.5_lag1, pm2.5_lag3, pm2.5_lag6, pm2.5_lag12, pm2.5_lag24, hour_sin, hour_cos, weekend] X df[feature_cols].values y df[pm2.5].values cut int(len(X) * 0.8) # 前 80% 时间段训练后 20% 时间段验证 X_train, X_test X[:cut], X[cut:] y_train, y_test y[:cut], y[cut:]这里刻意不用 sklearn 的 train_test_split避免有人顺手把 shuffleTrue 也带进来。按 cut 切分后训练集全部位于测试集之前这样测试集评估的是模型对未知未来的预测能力和真实部署场景保持一致。切分比例 8:2 在这个数据量下比较合理训练集约 3.4 万条足够稳定估计 14 个特征的回归系数。3.2 LinearRegression 的核心参数与拟合逻辑scikit-learn 的 LinearRegression 是机器学习回归任务里接口最精简的回归器之一拟合与评估的最小代码可以压缩到一个脚本里from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score # 默认 fit_interceptTrue模型带偏置项 model LinearRegression(fit_interceptTrue, positiveFalse) model.fit(X_train, y_train) y_pred model.predict(X_test) print(R2:, r2_score(y_test, y_pred)) print(RMSE:, np.sqrt(mean_squared_error(y_test, y_pred))) print(MAE:, mean_absolute_error(y_test, y_pred))fit_intercept 控制是否拟合偏置项 b在 PM2.5 场景里必须保留 True因为浓度存在一个非零的背景水平不拟合截距会把这部分压力强行压到某个特征系数上扭曲系数含义。positiveFalse 表示允许系数符号为负这个默认值不需要改——风速对浓度的作用在静稳天气和强风天气下表现不同强制非负反而违背真实关系。训练过程实际求解的是线性回归的正规方程闭式解数据量约 3.4 万行、特征维度 14闭式解的矩阵求逆耗时在毫秒级完全不需要像神经网络那样迭代调学习率。这也是线性回归作为 PM2.5 基线模型的一大优势可以在几秒内完成全流程验证把精力留给特征工程和残差诊断。3.3 评估指标与特征标准化经验三分钟看评估结果。R²、RMSE、MAE 三个指标最好一起看它们的含义和行为差异在 PM2.5 场景里非常明显。指标计算公式对错误的敏感度PM2.5 场景解读R²1 - SS_res/SS_tot反映整体解释力0.85 说明特征对浓度变化的解释力较好0.6 说明存在明显的模式没学到RMSEsqrt(mean((y-y_pred)²))对离群点敏感平方放大污染爆表时段预测偏差会被放大适合做安全预警MAEmean(y-y_pred)特征标准化这个变量在 LinearRegression 场景下常被误解。如果模型是裸的LinearRegression标准化完全不影响预测结果因为回归系数会自动缩放以抵消特征量纲差异预测值不变。但一旦引入 Ridge 或 Lasso正则化项对系数整体大小施加惩罚量纲大的特征如 PRES 气压值上千会被无差别压低这时必须先标准化再训练才能保证每个特征受到的惩罚强度一致。from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline from sklearn.linear_model import Ridge # 标准化 岭回归组合避免量纲影响正则化 model_ridge make_pipeline( StandardScaler(), Ridge(alpha1.0) ) model_ridge.fit(X_train, y_train)make_pipeline 的作用是把标准化和回归器串成一个整体对象fit 时先做标准化再做岭回归predict 时会自动重复同样的变换避免在测试集上单独调用 scaler.transform 时忘记使用训练集得到的均值和方差。alpha1.0 是最常见的起调值之后可以按 0.1、1、10 的对数刻度搜索观察验证集 R² 的峰值落在哪个区间。4. 线性回归模型诊断与特征优化让 R² 从 0.6 走到 0.94.1 残差分析线性回归失效的信号基线模型跑完第一个动作永远是看残差。残差是 y_test 减去 y_pred如果残差对预测值呈现出漏斗形或弯曲的分布说明模型存在系统性的低估或高估公式层面没问题也掩盖不了。import matplotlib.pyplot as plt residual y_test - y_pred # 残差图横轴预测值纵轴残差观察是否随机散落 plt.scatter(y_pred, residual, s4, alpha0.4) plt.axhline(y0, colorred, linestyle--) plt.xlabel(Predicted PM2.5 (μg/m³)) plt.ylabel(Residual (μg/m³)) plt.savefig(residual.png, dpi120)PM2.5 数据残差图几乎必然会展示出两种特征一是预测值超过 300 时残差发散真实高浓度时段被明显低估二是低浓度段残差贴近零轴说明模型对清洁天气拟合得不错。第一个现象的物理背景是重污染通常由静稳天气加高湿度共同触发线性回归只能学到这些变量的平均效应无法单独刻画极端组合条件。这类问题靠换模型不如靠补特征和调滞后窗口。4.2 用系数可解释性做特征筛选线性回归的可解释性在预测类任务里是稀缺品。训练完成后把系数输出成表格按绝对值排序能直接看出模型在依赖哪几个特征coef_df pd.DataFrame({ feature: feature_cols, coef: model.coef_ }).sort_values(coef, keyabs, ascendingFalse) print(coef_df.head(10))实际训练结果里pm2.5_lag1 的系数通常远高于其他特征这符合 PM2.5 的短时记忆效应风速 Iws 的系数显著为负也符合风大浓度低的物理直觉。如果某个特征系数绝对值非常小且训练集和验证集上的表现没有明显差异就可以把它从 feature_cols 中移除。特征数量从 14 降到 10 后模型复杂度下降残差分布往往更稳定。4.3 多项式扩展线性回归的另一种升级路径当四个气象变量与浓度之间不只是简单线性关系时可以引入交互项和平方项让线性回归获得非线性表达能力。这一步常见做法是直接上 PolynomialFeatures。from sklearn.preprocessing import PolynomialFeatures # degree2 生成一次项、平方项与两两交互项 poly PolynomialFeatures(degree2, include_biasFalse, interaction_onlyFalse) X_poly poly.fit_transform(X_train) X_test_poly poly.transform(X_test) model_poly LinearRegression() model_poly.fit(X_poly, y_train) y_pred_poly model_poly.predict(X_test_poly)degree2 时特征维度会从 14 膨胀到约 100 个原因是公式 n_features * (n_features 1) / 2 涵盖了所有平方项和交叉项。维度大幅升高后LinearRegression 的方差会上升容易在测试集上出现 R² 不升反降。这时的正确组合是第 3 章提到的 RidgeL2 正则化能约束交互项系数的范数让模型保留多项式表达能力的同时不剧烈过拟合。interaction_onlyTrue 则只生成交互项不给平方项适合已有明确先验温度与湿度的乘积影响浓度的场景。一个必须警惕的坑PolynomialFeatures 必须先在训练集上 fit再对测试集单独 transform不能对全量数据做扩展后再切分。因为拟合过程会记录每列的均值、方差和特征名称测试集信息一旦混入 fit 阶段评估结果就会虚高。这一点和 StandardScaler 的使用规则完全一致。4.4 用 TimeSeriesSplit 验证跨时段稳定性单次 8:2 切分只能说明模型在某一特定未来时段的表现。要验证模型在跨季节、跨年份是否稳定应该用时间序列交叉验证。sklearn 的 TimeSeriesSplit 专门解决这个问题from sklearn.model_selection import TimeSeriesSplit, cross_val_score # 5 折滚动验证每次用前 k-1 折训练第 k 折验证 tscv TimeSeriesSplit(n_splits5) scores cross_val_score(model, X, y, cvtscv, scoringr2) print(scores, scores.mean())TimeSeriesSplit 的折切分永远保持训练集在验证集之前n_splits5 会把全量数据切成 5 个连续时段依次做 5 次训练与验证。输出结果如果显示出前几折 R² 高、后几折 R² 低说明模型对近期数据的拟合能力在衰退原因是污染物组成和气象条件在年份间有缓慢漂移。反过来的趋势说明数据尾部存在极端污染事件模型预测能力被某个时段拉低。5. 交付离线预测系统模型持久化与定时预测脚本5.1 joblib 保存与加载模型的正确姿势模型调好之后训练脚本和预测脚本要分离。sklearn 官方推荐用 joblib 持久化模型对象它比 pickle 对大对象序列化更高效import joblib # 保存训练好的模型到本地文件 joblib.dump(model, pm25_lr_model.joblib) # 加载并直接用于预测 loaded_model joblib.load(pm25_lr_model.joblib)joblib.dump 的第一个参数是模型对象第二个参数是保存路径。加载后的模型不需要重新 fitpredict 方法与原先完全一致。需要提醒的是保存时顺手把 feature_cols 也序列化成一个 JSON 文件或者 joblib 文件加载预测时通过字段对比校验输入顺序能避免因后续新增特征而导致的预测时列错位问题。5.2 一个可直接排期的预测脚本骨架预测系统的最终形态通常是一个每小时被 cron 调起的脚本。它读取最近一个小时的观测值构造特征输出下一小时浓度预测并把结果追加写入日志import pandas as pd import numpy as np import joblib # 加载模型全局只加载一次 model joblib.load(pm25_lr_model.joblib) def build_features(obs): # 输入为一条最新观测记录输出与训练特征顺序一致 return [obs[DEWP], obs[TEMP], obs[PRES], obs[Iws], obs[Is], obs[Ir], obs[pm2.5_lag1], obs[pm2.5_lag3], obs[pm2.5_lag6], obs[pm2.5_lag12], obs[pm2.5_lag24], np.sin(2 * np.pi * obs[hour] / 24), np.cos(2 * np.pi * obs[hour] / 24), int(obs[weekday] 5)] latest pd.read_csv(latest_weather.csv).iloc[-1] features build_features(latest) pred model.predict([features])[0] # 追加写入预测结果供后续自检 with open(predictions.log, a) as f: f.write(f{latest[datetime]},{pred:.1f}\n)build_features 函数里的字段顺序必须与训练时的 feature_cols 完全一致这是预测脚本最容易出错的地方。latest_weather.csv 由上游的数据采集服务写入脚本只读最后一行不关心采集细节这样的设计把数据流向拆成采集 → 预测 → 落盘三段任何一段都可以独立替换。追加写日志的方式简单但可靠之后可以用同一组历史数据画时间序列曲线。注意cron 每小时调用脚本时需要保证滞后特征使用的历史浓度是已观测到的真实值而不是上一轮的预测值。使用预测值做滞后输入误差会在小时尺度上发生累积也就是常说的反馈误差漂移。5.3 预测结果自检用最后 7 天滑动 MAE 判断模型漂移系统部署后我保留一个固定位置的最近 7 天窗口每天把模型输出的预测值与真实监测值做一次对比算滑动平均绝对误差超过阈值就告警。原因是环境监测传感器会缓慢漂移气象输入分布也会随季节改变模型第一天表现好不代表第二个月仍然稳定。df_eval pd.read_csv(predictions.log, names[datetime, pred]) df_real df_recent[[datetime, actual]] merged pd.merge(df_eval, df_real, ondatetime, howinner) merged[err] (merged[pred] - merged[actual]).abs() # 输出最近 7 天168 小时的平均绝对误差 print(merged.tail(7 * 24)[err].mean())当结果高于 35 μg/m³先把最近几天大风、静稳等天气事件标出来核对如果误差集中出现在固定时段优先怀疑数据采集端的时间戳对齐问题而不是模型本身。用这样一个轻量自检脚本离线预测系统就具备了基本的可观测性后续再往里面加新的污染源特征也能用同样的滑动窗口做 A/B 对比评估新特征的真实增益。本文还有配套的精品资源点击获取