简介：这是一份基于线性回归的PM2.5预测系统Python源码，适合机器学习初学者、数据分析课程设计者参考。资源共19个文件，压缩包大小2.14MB，其中12个csv文件包含训练集、测试集以及经特征工程处理后的中间数据（如特征拼接表、标准答案），5个png图片直观展示数据分布或预测效果，另有model.npy保存训练得到的模型参数，主程序.py实现从数据读取、清洗、特征构造到模型训练与预测输出的完整流程。源码以繁体中文编码的空气质量观测数据为输入，逐步骤演示了big5文件读取、numpy矩阵运算、最小二乘求解及csv结果导出等关键操作，并附带预测结果样例，便于逐行对照理解线性回归在空气污染预测中的实际应用。已有594人学习下载，可据此快速复现实验，也可结合自身数据替换修改，用于课程报告或毕业设计。
1. 用线性回归做PM2.5预测：这个老模型为什么仍值得跑一遍
去年给某地级市做空气质量数据分析时，对方上来就问能不能上深度学习，我反问一句：你们现在有能用的基线吗？对方摇头。于是我先把历史监测数据喂给了线性回归，结果这个十几分钟跑完的模型，在预报次日均值上的表现直接超过了他们之前人工估算的方案。线性回归预测PM2.5的意义在于：它是所有预测系统的基线，也是检验特征工程是否有效的试金石。模型本身没有秘密，秘密全在特征和数据处理上。这套源码适合两类人——想入门机器学习、需要一个能跑通全流程的完整案例的新手，以及需要快速搭建业务基线、给领导一个可解释结论的工程师。
2. 把监测数据整理成特征表：时间戳、滞后项与缺失值处理
2.1 先看数据和字段：这套源码最常见的输入格式
PM2.5预测系统无论代码怎么写，第一步永远是数据整理。绝大多数监测平台导出的原始数据是CSV或Excel表格，每行是一个整点记录，列包含时间戳、PM2.5浓度、PM10、SO2、NO2、CO、O3以及气象要素如温度、湿度、风速、气压。这套源码里我一般会把数据文件命名为data.csv，字段名做成英文小写，避免中文列名在pandas里反复出幺蛾子。
import pandas as pd df = pd.read_csv("data.csv", parse_dates=["time"]) df = df.sort_values("time").reset_index(drop=True) print(df.head()) print(df.dtypes)逻辑说明：parse_dates把时间列解析成 datetime 类型，sort_values按时间排序。很多原始数据是乱序的，如果跳过排序，后面做滞后特征时shift会把不同日期的数据混在一起，模型立刻就废了。这一步是整个数据预处理的第一道关卡。
值得提醒的是，现场对接的数据经常混着小数点和千分位逗号，或者PM2.5浓度列里有"—"这种占位符。读取时应当加na_values=["—", "", "NULL"]把这些统一转成缺失值，而不是让pandas把它们读成字符串。源码包里如果附带了数据说明文档，第一件事永远是确认字段含义和单位，而不是急着跑模型。
2.2 构造时间与滞后特征：让线性回归「看见」昨天的风
线性回归没有记忆能力，它看到的每一行都是独立的。如果只给模型当前时刻的PM10和气象数据，模型完全不知道昨天发生了什么，预测效果会非常差。因此要做滞后特征（lag features），把过去几小时的浓度、过去24小时的平均值作为新列加入。
for lag in [1, 3, 6, 24]: df[f"pm25_lag_{lag}"] = df["pm25"].shift(lag) df["pm25_lag24_mean"] = df["pm25"].rolling(24).mean().shift(1) df["hour"] = df["time"].dt.hour df["weekday"] = df["time"].dt.weekday df = df.dropna().reset_index(drop=True)逻辑说明：shift(lag)表示取前 lag 行的数值，比如pm25_lag_1就是上一小时的PM2.5浓度。rolling(24).mean().shift(1)是过去24小时滑动平均值，注意shift(1)必不可少，否则当前时刻的浓度会被包含进去，造成数据泄漏。hour和weekday让模型能捕捉早晚高峰和周末效应。时刻、滞后浓度、滑动均值这几类特征组合起来，线性回归的效果会有明显跃升。
参数怎么调：滞后阶数的选择取决于业务场景。做小时级预测时，lag 取 1、3、6、24 基本够用；做次日日均值预测时，lag 应该取前一天相同时刻的值，并用过去7天的均值做平滑。这个参数没有固定最优解，可以先用少量组合跑一轮，看RMSE变化再决定。
2.3 缺失值与异常值处理：坏数据才决定RMSE下限
监测数据最常见的缺失发生在设备校准时段和通信中断时，表现为连续几小时为NaN。而异常值更典型：冬季某天PM2.5实测值忽然飙到800，可能是秸秆焚烧，也可能是传感器故障。线性回归对异常值极其敏感——因为它最小化的是平方误差，一个极端值就能把回归系数拉偏不少。
df["pm25"] = df["pm25"].replace(0, np.nan) df["pm25"] = df["pm25"].interpolate(method="linear", limit=6) df = df[df["pm25"] < 500]逻辑说明：第一行把0值替换成缺失，因为PM2.5浓度出现0基本是仪器归零错误；线性插值适合短时间缺失，超过6小时连续缺失建议用前后一天的同一时刻均值填充，而不是插值；df["pm25"] < 500是过滤异常大值，阈值按地区情况调整，华北地区冬季400多并不罕见，但500以上基本都是数据质量问题。处理缺失值时，我的习惯是保留一个is_missing标记列给模型，让它自己学习缺失时段与其他特征的关系，这对线性回归也是有效特征。
3. 写最小可跑的线性回归代码：sklearn 版与手写梯度下降版
3.1 sklearn 版本：十分钟跑通基线
拿到处理好的特征表之后，建模本身反而是最快的一步。这套源码里，train.py 核心代码不超过二十行。使用 scikit-learn 的LinearRegression是最靠谱的起点，它默认使用最小二乘法求解，数据量在几万行内几乎瞬时完成。
import pandas as pd from sklearn.model_selection import train_test_split from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, r2_score df = pd.read_csv("features.csv") feature_cols = [c for c in df.columns if c not in ["time", "pm25"]] X = df[feature_cols].values y = df["pm25"].values X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, shuffle=False ) model = LinearRegression() model.fit(X_train, y_train) y_pred = model.predict(X_test) rmse = mean_squared_error(y_test, y_pred, squared=False) print("RMSE:", round(rmse, 2)) print("R2:", round(r2_score(y_test, y_pred), 4))逻辑说明：这段代码的关键在shuffle=False。时间序列数据不能随机打乱后切分，否则测试集里混着训练集前后时间段的数据，模型等于开卷考试，RMSE会虚低。squared=False让mean_squared_error直接返回RMSE（均方根误差），单位与PM2.5浓度一致，比较好解释。
如果你刚装好Python环境，这一步需要先安装依赖：pip install pandas numpy scikit-learn。版本上 scikit-learn 1.0 以上即可，不需要追新。遇到ModuleNotFoundError时先确认是不是装进了别的虚拟环境，这个坑我见得太多了。
3.2 手写梯度下降版：理解学习率和收敛
sklearn 把细节都封装了，但作为工程落地，理解参数怎么调仍然绕不开底层原理。这套源码如果附带手写版本，一般会包含一个用梯度下降训练线性回归的文件，它能让读者直观看到学习率如何影响收敛。
import numpy as np class LinearRegressionGD: def __init__(self, learning_rate=0.01, epochs=1000): self.lr = learning_rate self.epochs = epochs self.weights = None self.bias = 0 def fit(self, X, y): n_samples, n_features = X.shape self.weights = np.zeros(n_features) for _ in range(self.epochs): y_pred = X @ self.weights + self.bias dw = (2 / n_samples) * X.T @ (y_pred - y) db = (2 / n_samples) * np.sum(y_pred - y) self.weights -= self.lr * dw self.bias -= self.lr * db def predict(self, X): return X @ self.weights + self.bias逻辑说明：核心是每次迭代计算损失函数对所有参数的偏导数，然后沿负梯度方向更新参数。dw是权重梯度，db是偏置梯度，learning_rate决定每一步走多远。数据量较大时，矩阵乘法X.T @ (y_pred - y)仍然是全量计算，这是批梯度下降。
参数怎么设：学习率 0.01 对标准化后的数据通常能稳定收敛；学习率 1.0 时损失函数会震荡甚至爆炸。特征没有做标准化时，不同特征的量纲差异会导致梯度方向偏斜，收敛慢。因此使用手写版本前，务必对特征做标准化，或者把学习率调到 0.0001 级别，否则loss曲线会像心电图一样乱跳。
3.3 源码工程的目录组织：train.py / predict.py / features.py
一个能交付的预测系统源码，文件组织至少要让人一眼看出「数据从哪来、特征在哪造、训练怎么跑、预测怎么调」。这套源码包里的典型结构是：
| 文件 | 职责 |
|---|---|
load_data.py | 读取原始CSV，字段清洗，返回干净的DataFrame |
features.py | 滞后特征、滑动平均、时间特征构造 |
train.py | 训练模型，输出RMSE、R²，保存模型到本地 |
predict.py | 加载模型，对最新一条数据做下一时刻预测 |
model.pkl | 训练产出的模型文件，predict 阶段反序列化使用 |
我一般会在predict.py里也放一份特征构造代码，保证预测时的输入格式与训练时完全一致。很多翻车现场就是训练时用了pm25_lag_1，预测时忘了构造这一列，导致维度对不上直接报错。如果源码包的README里没有强调这一点，你接手后要主动把它对齐。
4. 模型评估与参数调优：R²、RMSE 与三个必调参数
4.1 评价指标怎么选：RMSE 惩罚大偏差，MAE 更稳
预测系统的评价指标决定了调参方向。PM2.5预测里，最常用的是RMSE和MAE。RMSE把误差平方后求平均再开方，因此它对那些偏差特别大的预测点非常敏感。比如某天实际浓度250，模型预测150，这个100的误差在平方之后会主导整个RMSE数值。MAE则只看绝对误差的平均值，更平稳，但对大偏差不敏感。
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score rmse = mean_squared_error(y_test, y_pred, squared=False) mae = mean_absolute_error(y_test, y_pred) r2 = r2_score(y_test, y_pred) print(f"RMSE: {rmse:.2f} MAE: {mae:.2f} R2: {r2:.4f}")逻辑说明：三个指标一起打印，能看出模型误差分布的特征。如果RMSE明显大于MAE，说明存在少数预测得很离谱的样本；如果两者接近，说明误差分布比较均匀。对于PM2.5业务而言，污染峰值时段的预测偏差才是考核重点，因为重污染天的应急响应才是刚需。所以我在项目里会额外关注「实际浓度>200时模型的RMSE」，这个指标比全局RMSE更能说明问题。
评估代码里还要关注r2_score。R²接近1说明模型解释了大部分方差，但R²低不一定意味着模型没用——如果测试集正好赶上一次极端沙尘天气，R²会很难看，可模型在常规天气下其实表现不差。所以判断模型质量时，不要只看单个指标，要结合时间序列图观察预测曲线与实际曲线的跟随程度。
4.2 三个必调参数：滞后窗口、训练集比例、是否加正则
线性回归的参数不多，但这三个地方值得花时间调。首先是滞后窗口的长度。滞后项太少，模型缺少历史信息；滞后项太多，特征之间相关性变强，数值稳定性变差。我的经验是小时级预测用[1, 3, 6, 24]四档，日均值预测用[24, 48, 72]加7日滑动平均，在RMSE上会有明显差异。
第二个是训练集与测试集的比例和切分方式。时间序列模型不能用随机切分，而是按时间顺序切。常见做法是前80%训练、后20%测试，这样测试集模拟的是「用历史预测未来」的真实场景。如果你想严谨一点，可以再用最后两周数据做验证集，用于调参，避免在测试集上反复调参造成信息泄漏。
第三个是是否使用正则化。当特征数量多、共线性严重时，普通线性回归的系数方差会变大。此时把LinearRegression换成Ridge加上 L2 正则，往往能让测试集上的RMSE下降几个点。正则化系数alpha的典型值从 0.01 到 10 之间网格搜索，选验证集上RMSE最小的值，这是性价比最高的调参操作。
from sklearn.linear_model import Ridge from sklearn.model_selection import GridSearchCV param_grid = {"alpha": [0.01, 0.1, 1.0, 10.0]} ridge = Ridge() grid = GridSearchCV(ridge, param_grid, cv=5, scoring="neg_mean_squared_error") grid.fit(X_train, y_train) print("best alpha:", grid.best_params_)逻辑说明：GridSearchCV在训练集内部做交叉验证来选超参数，而不是直接用测试集评估。因为时间序列的连续性，这里cv=5用的是默认的KFold划分，严格来说会在交叉验证中造成时间错位。更严谨的做法是用TimeSeriesSplit代替默认KFold，这点后面避坑章节会细说。neg_mean_squared_error是sklearn的常规写法：sklearn的评分规则是越大越好，所以MSE要取负号。
4.3 用交叉验证验证稳定性：别让一次划分决定结果
单次划分训练集和测试集具有偶然性——测试集恰好落在重污染过程或恰好全是好天，结果天差地别。我一般会用TimeSeriesSplit做滚动验证，把数据切成多段，每段都用前面的数据训练、后面的数据测试，综合多次结果看稳定性。
from sklearn.model_selection import TimeSeriesSplit tscv = TimeSeriesSplit(n_splits=5) rmse_list = [] for train_idx, test_idx in tscv.split(X): X_tr, X_te = X[train_idx], X[test_idx] y_tr, y_te = y[train_idx], y[test_idx] m = LinearRegression() m.fit(X_tr, y_tr) pred = m.predict(X_te) rmse_list.append(mean_squared_error(y_te, pred, squared=False)) print("RMSE mean:", round(np.mean(rmse_list), 2), "std:", round(np.std(rmse_list), 2))逻辑说明：TimeSeriesSplit保证每次划分中训练集永远在测试集之前，符合时间序列预测的真实约束。输出5次RMSE的均值和标准差：均值是模型的真实水平，标准差是稳定性。如果标准差很大，说明模型在不同时间段表现波动大，需要回头检查特征是否遗漏了季节性信息。
5. PM2.5预测避坑指南：这些坑我都替你踩过
5.1 数据泄漏：把未来数据混进了特征
现象：训练集上R²高达0.98，测试集也接近0.97，模型表现好得不真实，部署到线上后预测值却一直滞后。
原因：构造滞后特征时用了当前时刻的数据算滑动平均，比如用当天的PM2.5参与计算当天的预测特征。测试阶段没问题，因为特征是从同一时刻数据里算出来的，但线上预测时你还没有当前时刻的PM2.5真实值，特征构造不完整，模型只能使用残缺输入。
解决：所有特征必须来自「预测时刻之前」的数据。滑动平均构造完要shift(1)，确保特征矩阵里第 i 行的数据只包含第 i 行以前的信息。写代码时养成一个习惯：构造完特征后，随机抽一行，手动验证一下该行特征是否都能用历史数据独立算出。
5.2 用 R² 判断模型好坏会误判
现象：模型A的R²为0.90，模型B的R²为0.82，团队选择上线A。实际运行后发现A在高浓度时段的误差巨大。
原因：R² 衡量的是模型解释的方差比例，对均值附近的拟合程度非常敏感。PM2.5数据大部分时间在100以下，模型只要把低浓度时段拟合好，R²就会很高；高浓度时段样本少，对R²贡献小，模型即使在这些时段完全失效也不影响R²。
解决：用RMSE作为调参主指标，同时单独统计重污染时段（实测值>150）的RMSE。只有这两个指标同时优化，模型才有业务价值。在汇报时也习惯把两个数字一起放出来，避免其他人被单一指标带偏。
5.3 时间序列数据不能用随机切分
现象：train_test_split默认test_size=0.2且shuffle=True，跑出来的RMSE远好于预期，模型看起来完美知道未来趋势。
原因：随机打乱后，训练集和测试集的数据来自同一时间段，测试集里包含了与训练集相邻时刻的信息。比如测试集取了下午3点的数据，训练集里可能就有下午2点的数据。时间序列存在强自相关性，模型相当于在「默写」测试集。
解决：一律使用shuffle=False，按时间先后切成前80%和后20%。更严格的方案是用TimeSeriesSplit滚动验证，并保证训练集和测试集之间留出至少一天的间隔（gap），避免滑动平均特征跨越切分边界造成泄漏。
5.4 缺失值填充错误让滞后特征「空转」
现象：数据缺失率约15%，直接用dropna()删掉所有含缺失的行，结果剩下的数据集中在秋冬季，模型在夏季表现一塌糊涂。
原因：缺失不是均匀分布的——设备检修常常集中在固定的时间段，直接删行相当于重塑了数据分布。另外，如果先做滞后特征再删行，一行缺失会导致后面多个滞后特征全变NaN，删除量比预期大得多。
解决：先填充缺失值，再做滞后特征。短时间缺失用线性插值，长时间缺失用前后同时间点的均值填充；实在不行就保留NaN，让树模型类算法处理。线性回归不能接受NaN，所以必须在特征构造前把缺失值解决干净。
5.5 极端浓度把回归系数拉偏
现象：模型各项指标看起来可以，但某次沙尘天气后重新训练，RMSE整体变大，且预测值系统性偏高。
原因：一次极端事件产生的样本点在损失函数里被平方放大，为了拟合这少数几天，模型牺牲了大量常规天气样本的精度。我在实际项目中遇到过PM2.5实测值飚到900多的记录，事后确认是传感器故障，但模型已经为此付出代价。
解决：建模前画分布图，把明显超过物理上限（比如500）的样本挑出来人工核对。确属真实污染事件时可以保留，但要用截断或加权方式降低其影响，比如把极端值截断到某个百分位数。对业务预测来说，保住常规天气的精度比死磕极端值更有价值。
6. 进阶：从单步预测到多步预测的滞后窗口策略
前面的模型预测的是下一小时的PM2.5浓度，属于单步预测，也是这套源码的核心功能。但实际业务中常常要的是未来24小时或未来3天的浓度变化曲线，这时候单纯调参解决不了，需要换预测策略。
两种常见做法：一种是递推预测，即把上一时刻的预测值当作历史值放回特征里，继续预测下一步——优点是代码简单、依赖少，缺点是误差会逐级累积，预测周期越长越不靠谱；另一种是直接多步预测，直接训练24个独立模型，分别预测未来1小时、2小时到24小时的浓度——每个模型各自构建对应的滞后特征，缺点是训练成本增大，但避免了误差累积问题，我在项目里更倾向这种方式。
def train_direct_multistep(df, horizons=[1, 6, 24]): models = {} for h in horizons: df[f"target_{h}"] = df["pm25"].shift(-h) tmp = df.dropna() X = tmp[[c for c in tmp.columns if c.startswith("feat_")]].values y = tmp[f"target_{h}"].values m = LinearRegression() m.fit(X, y) models[h] = m return models逻辑说明：对每个预测步长h，把目标列向上移动h行，即第 t 行的目标是 t+h 时刻的浓度。注意shift(-h)后尾部会多出 h 行NaN需要删掉。每个模型都用自己对应的目标列训练，预测时分别调用不同模型得到未来曲线上的离散点。
我个人的经验是，做小时级预测系统时，递推预测能用就用，省事；做决策支撑（比如是否启动重污染应急预案）时，直接多步预测的可靠性更值得信赖。这套源码拿来在实际业务中跑，最重要的事不是追求R²再高零点几，而是把特征构造、数据切分、评估口径这三件基本功做扎实。很多团队翻车，不是模型不够强，是数据管线和评估方法出了问题，再花哨的算法也救不回来。希望帮到你。
本文还有配套的精品资源，点击获取
