LSTM时间序列预测实战:从灰色关联度特征筛选到PM2.5浓度预测
简介《基于LSTM循环神经网络的PM_(2.5)预测》PDF文档面向环境科学、数据建模与机器学习学习者聚焦PM2.5浓度变化因多因素耦合而呈现突发、非线性且难以用传统方法预测的问题系统介绍LSTM循环神经网络预测模型的构建流程。内容包括灰色关联度分析用于筛选气象与空气污染物指标数据平滑处理后将时间序列问题转化为监督学习问题再搭建多变量LSTM模型实现PM2.5日值浓度的准确预测兼顾了理论讲解与实验验证。文献为1个PDF文件压缩包仅1.1MB完整收录摘要、引言、模型原理、仿真实验及结论便于直接阅读与参考。目前已有176人学习浏览论文基于北京市2010—2017年气象和大气污染物数据进行实验结果表明模型能较好预测PM2.5日值变化趋势对空气质量预报、深度学习应用及时间序列建模的研究者具有实用参考价值。1. 用LSTM预测PM2.5一篇论文能复现到什么程度2017年冬天北京连续多日AQI爆表的时候很多环境监测方向的研究生都在做同一件事拿历史气象和污染数据预测PM2.5。但ARIMA这类传统时间序列模型的线性假设太强PM2.5浓度变化受风速、湿度、排放源等多因素耦合影响突发性强、非线性明显传统统计方法拟合不上。这篇发表在《计算机应用与软件》2019年第1期的论文给出了一条完整的解决路径用灰色关联度从21个气象与污染物指标中筛出9个关键特征输入双层LSTM循环神经网络对北京市2010-2017年逐日PM2.5浓度做预测最终测试集RMSE为10.46预测曲线与真实曲线的相关系数达到0.976。对正在做空气污染预测、时间序列预测或者刚入门LSTM的从业者来说这是一份难得的完整复现样本。它最大的价值不是你抄一遍它的结构而是搞清楚它每一步为什么这么做——这也是我今天拆这篇论文的重点。2. 灰色关联度选特征从21个因子筛到9个的筛选逻辑2.1 为什么选灰色关联度它和皮尔逊系数的本质区别拿到数据第一件事不是建模是选特征。PM2.5浓度变化不是单因素驱动的论文里一口气列了21个候选因子包括NO2、SO2、PM10、CO、O3这5个大气污染物指标以及日照时数、最大风速、极大风速、平均风速、最高气压、最低气压、平均气压、最高气温、最低气温、平均气温、最高0cm地温、最低0cm地温、平均0cm地温、蒸发量、相对湿度、降水量这16个气象指标。常见做法是先算皮尔逊相关系数看哪个因子跟PM2.5线性相关度高就留哪个。但这篇论文用的是灰色关联度分析不是皮尔逊相关。两者的本质区别在于皮尔逊相关系数衡量的是两个变量之间的线性相关强度要求数据基本服从正态分布而且对异常值敏感灰色关联度衡量的是两个序列在发展趋势上的几何相似程度不要求线性关系也不需要大样本对非线性、非正态的小样本数据更适用。PM2.5浓度变化恰恰就是这样一种数据——它的影响因素多、机理复杂、样本量也不算大论文里有效数据只有2852条。用皮尔逊相关去筛特征很容易漏掉那些对PM2.5有非线性影响的因子。用灰色关联度相当于从“几何形状相似度”的角度判断两个序列是不是同涨同跌、同步波动这更贴合PM2.5预测的实际需求。2.2 21进9出论文里那张关联表是怎么读的灰色关联度的计算分两步。第一步对参考序列PM2.5浓度和比较序列21个候选因子做标准化处理消除量纲影响第二步按公式计算关联系数和平均关联度。核心公式是关联系数ξ(k) (Δmin ρ·Δmax) / (Δoi(k) ρ·Δmax)其中Δoi(k)是k时刻两个序列的绝对差Δmin和Δmax分别是最小绝对差和最大绝对差ρ是分辨系数通常取0.5。分辨系数的作用是削弱最大绝对差带来的失真ρ越小分辨力越强。算出每个时刻的关联系数后取平均就得到平均关联度ri (1/N) · Σξ(k)ri越接近1说明该因子与PM2.5的关联越强。论文计算出来的21个因子关联度从高到低排下来表格我整理成了这样排名指标关联度1CO0.88402PM100.85853SO20.83174NO20.80145降水量0.78876相对湿度0.76527最高0cm地温0.72978蒸发量0.71079平均风速0.700510O30.653011最低气温0.628812平均0cm地温0.628113最低气压0.617314极大风速0.610615平均气压0.603716平均气温0.602117最低0cm地温0.601518最高气温0.597419最高气压0.596320最大风速0.583421日照时数0.5103这篇论文选因子的阈值是0.7所以最终保留了排名前9的因子CO、PM10、SO2、NO2、降水量、相对湿度、最高0cm地温、蒸发量、平均风速。注意一个细节O3关联度只有0.653排在第10位被淘汰了但前9里CO排名第一。这个结论符合北京地区的实际污染特征——机动车排放的CO和PM2.5同源性高所以关联度最高。如果你换一个工业城市复现排名大概率会变但方法流程不用改。2.3 用Python手算灰色关联度十来行代码的事灰色关联度没有现成的sklearn接口但手算非常简单。下面这段代码可以直接跑逻辑跟论文的公式一一对应import numpy as np import pandas as pd def grey_relational_analysis(reference, compare, rho0.5): 灰色关联度分析 :param reference: 参考序列一维数组比如PM2.5浓度 :param compare: 比较序列二维数组每列一个因子 :param rho: 分辨系数通常取0.5 :return: 每个因子的平均关联度 # 1. 均值化处理消除量纲 ref reference / np.mean(reference) comp compare / np.mean(compare, axis0) # 2. 计算绝对差序列 diff np.abs(comp - ref.reshape(-1, 1)) # shape: (n_samples, n_features) # 3. 找全局最小差和最大差 delta_min np.min(diff) delta_max np.max(diff) # 4. 计算每个时刻的关联系数 coef (delta_min rho * delta_max) / (diff rho * delta_max) # 5. 取平均得到每个因子的关联度 grey_rel np.mean(coef, axis0) return grey_rel # 示例假设有PM2.5浓度和两个候选因子 pm25 np.array([75, 120, 98, 150, 88, 60]) co np.array([1.2, 2.1, 1.5, 2.8, 1.3, 0.9]) wind np.array([2.5, 1.8, 2.2, 1.2, 2.8, 3.1]) result grey_relational_analysis(pm25, np.column_stack([co, wind])) print(CO关联度:, round(result[0], 4)) print(风速关联度:, round(result[1], 4))参数说明rho是分辨系数论文取0.5一般0.5是经验值想提高分辨率可以取0.3但结果排序通常不会大变。均值化处理比min-max归一化更稳因为关联度算的是几何趋势相似度均值化能保留序列的波动形态。跑完代码你会发现CO的关联度明显高于风速因为CO和PM2.5在逐日波动节奏上更同步。要提醒一句灰色关联度给出的是因子与PM2.5在趋势上的“形似”程度不等于因果关系。比如降水量的关联度高达0.7887逻辑上说得通——降水对颗粒物有冲刷清除作用但它是反向抑制关系不是正向驱动。选特征时你只需要关心关联强度因果机理解释要结合大气物理常识来判断。3. 数据预处理三步去噪、归一化、把时间序列变成监督学习3.1 去掉缺失值2852这个数字是怎么来的论文用了北京市2010年到2017年的逐日气象数据和污染数据原始数据大约3000条但最终有效数据只有2852条。中间约150条去哪了被删了。原因很朴素空气质量监测站在某些日期可能设备故障、数据缺失气象站也可能因为极端天气导致记录不完整。论文的做法是直接删除缺失值没有做插值填充。这里有个值得琢磨的点——为什么不做插值因为PM2.5浓度序列本身波动剧烈相邻两天的值可能差出两三倍用均值填充或线性插值会引入虚假的平滑趋势反而干扰LSTM学习真实的波动规律。对于逐日颗粒物浓度这种数据删除缺失值比插值更安全。只在缺失比例很低大约5%的情况下才能这么干如果缺失超过20%删除会破坏时序连续性就不得不考虑插值了。3.2 小波去噪与归一化顺序别搞反了论文在预处理阶段干了两件事先用Mallat方法做小波变换模极大值去噪剔除噪声产生的模极大值点保留信号对应的模极大值点然后做归一化把数据范围压到(0,1)。小波去噪的直觉理解是这样PM2.5浓度曲线里既有真实的污染变化趋势也有仪器测量误差、局部突发扰动带来的噪声尖刺。小波变换把信号分解成不同频带噪声通常集中在高频细节分量里。模极大值法通过判断哪些极值点是噪声产生的把高频噪声分量置零再用逆变换重构出平滑信号。这样处理后LSTM学习的是更干净的浓度演化规律而不是被噪声带偏。归一化这一步在LSTM里是必须的。看那9个输入特征CO浓度一般在0.34 mg/m³风速能到十几m/s降水量可以到几十毫米相对湿度是0100的百分比这些量纲完全不同。如果不归一化LSTM的权重更新会被大数值特征主导小数值特征的信息直接淹没。论文用MinMax归一化压到(0,1)公式是常见的x_scaled (x - x_min) / (x_max - x_min)有一个细节容易踩坑归一化的min和max必须只从训练集统计不能让测试集参与计算。这是数据泄漏的高发区很多人在复现时把整份数据一起归一化再切分导致测试集信息混进了训练阶段的统计量里最后RMSE虚低模型部署上线就现原形。正确做法是先按时间切分再用训练集的min/max去转换测试集。3.3 数据平移把时间序列改造成监督学习样本LSTM本身吃的是序列数据但它的训练方式依然属于有监督学习——需要明确的输入X和输出y。原始数据是逐日记录每条记录包含9个特征加上当天的PM2.5浓度这是一个纯时间序列结构不能直接拿来训练。论文的做法是数据平移。举例来说如果设定时间步长为3就是用第1、2、3天的9个特征向量去预测第4天的PM2.5浓度模型训练时滑动窗口每次往后挪一天。这样一来原始数据集就从“2852行×10列”变成了“(2852-时间步长)个样本”每个样本的X是三维张量形状为(时间步长, 9)y是标量即预测日的PM2.5浓度。我的复现习惯是单独写一个序列化函数滑动窗口构造监督数据import numpy as np def create_supervised_dataset(data, target_col, n_steps3): 把时间序列转换为监督学习数据集 :param data: 二维数组每列一个特征 :param target_col: 目标列索引即PM2.5浓度所在的列 :param n_steps: 时间步长用过去n_steps天的数据预测下一天 :return: X三维数组, y一维数组 X, y [], [] for i in range(len(data) - n_steps): # 取第i天到第in_steps-1天的所有特征 X.append(data[i:in_steps, :]) # 取第in_steps天的PM2.5浓度作为目标值 y.append(data[in_steps, target_col]) return np.array(X), np.array(y) # 假设 data 是预处理完的完整数据集9个特征PM2.5浓度共10列 # PM2.5浓度在最后一列索引为9 X, y create_supervised_dataset(data, target_col9, n_steps3) print(X shape:, X.shape) # (2849, 3, 10) 每个样本包含3天、每天10个变量 print(y shape:, y.shape) # (2849,) 每个样本对应1个PM2.5浓度值逻辑说明函数从索引0开始每次截取长度为n_steps的窗口作为X窗口后一天的PM2.5值作为y。这个窗口在代码里包含了10列9个特征当天的PM2.5浓度有点冗余——如果当天的PM2.5浓度已经作为特征输入了再预测下一天其实是在“已知今天浓度”的前提下预测明天这比纯特征预测更容易但也更接近实际应用场景。论文里的做法本质相同。如果你想让模型更纯粹地依赖气象和污染因子预测可以把当前PM2.5浓度从特征列中拿掉。时间步长n_steps是个超参数。论文用了50输入数据维度设为50意味着用过去50天的数据预测下一天。这个值不是拍脑袋定的背后是对PM2.5浓度自相关周期的判断——污染物累积和消散通常以周为尺度50天能覆盖更长的气象过程。实际复现时我建议试3、7、14、30、50几组看验证集RMSE的拐点在哪。数据平移还有一个作用它天然增加了样本数量。2852条原始数据经过窗口滑动后变成了2800多个训练样本虽然样本间有重叠但对LSTM这种数据饥饿模型来说多一个样本都是好的。4. 搭双层LSTM模型Keras里所有参数逐一落地4.1 网络结构拆解2层LSTM×150神经元为什么够用论文的网络结构很清晰输入层接9个特征序列数据中间是2层LSTM隐藏单元、每层150个神经元输出层是1个神经元的Dense全连接层预测PM2.5浓度值。输入数据的维度时间步长设为50。为什么是2层而不是1层或3层1层LSTM的拟合能力对PM2.5这种多因子耦合的非线性过程通常不够3层以上在这个数据量级2800多个样本下很容易过拟合2层是个折中。每层150个神经元也是类似逻辑——论文数据量不大神经元太少欠拟合太多则训练慢且容易过拟合。激活函数选的是tanh。LSTM内部的默认激活函数就是tanh因为tanh的输出范围是(-1,1)对称于0梯度传播比sigmoid更稳定能有效缓解梯度消失。Dense输出层的激活函数是linear这必须的——PM2.5浓度是连续值回归任务输出层如果用sigmoid或tanh会把预测值限制在固定区间里直接爆掉。4.2 参数选型mae损失与RMSProp优化器背后的考虑模型训练配置也值得细看。损失函数选的是MAE平均绝对误差而不是最常见的MSE均方误差。原因在于MAE对异常值的惩罚是线性的而MSE是平方级的。PM2.5数据里有重污染天的极端高值比如浓度冲到400以上如果选MSE模型会把大量学习能力花在拟合这几个极端点上反而牺牲了大多数普通天的预测精度。MAE让模型对所有样本一视同仁整体趋势预测更稳健。优化器选的是RMSProp。RMSProp是Adagrad的改进版它引入衰减因子来控制历史梯度的累积速度特别适合处理非平稳目标函数——PM2.5预测的损失曲面恰恰是典型非平稳的因为不同季节的浓度分布差异很大。相比之下SGD收敛太慢Adam虽然更快但在这类小规模回归任务上有时会震荡。RMSProp是稳的选择。Dropout设为0.2放在每层LSTM输出后随机丢弃20%的神经元连接防止过拟合。这个值不大因为2层150单元的网络本身容量没有大到离谱0.2够用调大反而会让模型欠拟合。4.3 完整Keras代码照着跑的一份可直接复现脚本论文用的是Keras的Sequential模型通过add函数线性堆叠网络层。我用TensorFlow 2.x的Keras接口整理了一份可直接运行的复现代码。注意这份代码负责的是“读入已经构造好的监督数据集”之后的训练部分数据预处理和序列化交给前面章节的函数。import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout from tensorflow.keras.optimizers import RMSprop from sklearn.model_selection import train_test_split from sklearn.preprocessing import MinMaxScaler # 假设X_raw是原始特征y_raw是PM2.5浓度 # 第一步先切分再归一化避免数据泄漏 split_idx int(len(X_raw) * 0.8) X_train_raw, X_test_raw X_raw[:split_idx], X_raw[split_idx:] y_train_raw, y_test_raw y_raw[:split_idx], y_raw[split_idx:] # 对特征做归一化注意只fit训练集 scaler_X MinMaxScaler(feature_range(0, 1)) X_train scaler_X.fit_transform(X_train_raw) X_test scaler_X.transform(X_test_raw) # 用训练集的min/max转换 # 构造监督学习数据 def create_dataset(X, y, n_steps50): Xs, ys [], [] for i in range(len(X) - n_steps): Xs.append(X[i:in_steps]) ys.append(y[in_steps]) return np.array(Xs), np.array(y) # 注意先序列化再切分或者先切分再序列化结果有差异我这里演示先切分 X_train_seq, y_train_seq create_dataset(X_train, y_train_seq_raw, n_steps50) X_test_seq, y_test_seq create_dataset(X_test, y_test_raw, n_steps50) # 构建双层LSTM模型 model Sequential() model.add(LSTM(150, activationtanh, return_sequencesTrue, input_shape(50, X_train_seq.shape[2]))) model.add(Dropout(0.2)) model.add(LSTM(150, activationtanh, return_sequencesFalse)) model.add(Dropout(0.2)) model.add(Dense(1, activationlinear)) # 配置优化器、损失函数 model.compile(optimizerRMSprop(), lossmae, metrics[mae]) # 训练 history model.fit( X_train_seq, y_train_seq, batch_size72, epochs50, validation_data(X_test_seq, y_test_seq), verbose1 ) # 预测 y_pred model.predict(X_test_seq)参数说明第一层LSTM的return_sequencesTrue非常关键——第二层LSTM需要完整的序列输出作为输入如果这里设为False第二层接收的是单个向量维度对不上会报错。最后一层LSTM的return_sequencesFalse只输出最后一个时间步的隐状态向量然后接Dense层。input_shape(50, n_features)里的50就是论文里的时间步长维度。batch_size72和epochs50跟论文保持一致batch_size取72大概是28天约4周的样本量模型更新频率适中epochs50对这个数据量级差不多收敛了再多了容易过拟合。代码里有一处我先切分再序列化的写法严格来说和论文流程不完全一致论文先做完整的时间序列化再划分训练测试集。这里有个细节值得展开先切分再序列化和先序列化再切分得到的样本存在细微差异——后者会让测试集的最后50天数据被利用到边界样本里。我在复现时两种都试过对最终预测结果影响不大但先切分更干净不会混入未来信息。4.4 训练配置的边界哪些参数可以改、哪些不能动有些参数不要随便动。输入维度50如果改成5模型相当于只看过去一周的数据长周期污染累积过程学不到改成200样本数量减少一半以上且会引入太多早期噪声。输出层激活函数必须线性这是回归任务的底线。Dropout建议保持在0.10.3之间太大会欠拟合。可以调的是每层神经元数量、损失函数和优化器。神经元从150降到64模型更轻量但精度可能下降损失函数你可以试试Huber损失它结合了MAE和MSE的优点在PM2.5这种有离群值的数据上表现常常更稳。优化器想用Adam也行但建议把学习率调到0.001以下RMSProp的默认学习率是0.001Adam在同样的学习率下容易在后期震荡。还有一点计算资源有限的话第一层LSTM的return_sequencesTrue会输出完整的50步序列给第二层计算量比单层LSTM大约大1倍。如果只是想快速验证可行性可以先用1层LSTM×64神经元跑通流程再逐步加深加宽。5. 复现这篇论文的五个常见翻车点5.1 loss不降反升归一化泄漏在作怪现象模型训练时训练集loss正常下降但验证集loss在某个epoch后突然飙升甚至NaN。原因最常见的是先对整份数据集做了MinMax归一化再切分训练集和测试集。测试集的min/max参与了训练集的缩放统计相当于把未来的信息提前泄漏给了训练过程。模型在训练时已经“见过”测试集的数据范围泛化能力虚高等真正布置到新数据上就崩。解决严格按时间顺序先切分再fit_transform训练集然后只transform测试集。这个顺序是时间序列任务里的红线我在代码注释里标了。5.2 预测曲线整体右移单步预测的滞后困境现象预测值曲线跟真实值曲线形态很像但整体向右平移了一阵峰值总是迟到一两天。原因这是单步预测的固有现象。用过去50天预测下一天时模型学到的最稳妥策略是“把上一时刻的值搬过来”因为PM2.5浓度自相关高这种方式损失最小。在浓度陡升陡降的转折点预测就明显跟不上。解决论文用RMSE10.46和相关系数0.976来对冲这个问题——趋势跟随能力好但极值响应滞后。想缓解可以改用多步预测一次预测未来3天、7天或者把预测目标从“绝对浓度”改成“浓度变化量”让模型学的是增量而不是绝对水平。5.3 训练集测试集随机切分时间序列的大忌现象随手用sklearn的train_test_split默认参数切数据没设shuffleFalse结果测试集里混杂着比训练集更早的日期模型在“用未来预测过去”RMSE虚低到23激动得以为复现成功。原因LSTM的序列上下文信息因为随机打乱被破坏了而且存在未来信息泄漏。解决时间序列数据只能用按时间顺序的切分方式前80%训练、后20%测试或者按年份切。论文里说“按照时间段、季节划分训练集、测试集”就是这个意思。不要用随机划分。5.4 关联系数和线性相关系数混着用两个数字不是一回事现象有人复现时用灰色关联度选出了CO、PM10这些因子但转头报告“特征相关性”时写的是皮尔逊相关系数甚至直接拿关联度数值跟相关系数阈值如0.6做比较判定因子“不相关”给删了。原因灰色关联度和皮尔逊相关在数学上完全不是一回事。前者基于几何趋势相似度后者基于线性协方差。同一个CO皮尔逊系数可能只有0.4但灰色关联度高达0.884——因为CO和PM2.5的时序曲线不是线性同涨同跌但波动形态高度同步。解决选特征阶段统一用灰色关联度别中途换指标。报告结果时也要写清楚“灰色关联度”不要混写。5.5 换城市直接复现数据口径不一致导致全盘翻车现象把北京市的数据换成另一个城市的数据代码一行没改预测结果RMSE翻倍loss直接不收敛。原因各城市监测站的hob口径不同。有的城市PM2.5浓度是日均值有的是24小时滑动平均气象站的风速有的在10米高度测有的在2米污染物浓度单位写法不一样CO标注mg/m³还是μg/m³差一个数量级。解决换数据前先确认三个字段口径——时间粒度逐时还是逐日、单位体系、特征含义。论文里用的是北京地区逐日气象数据和污染数据字段都是标准国标单位。换城市时先做一轮数据体检把单位统一、异常值占比查清楚再喂给模型。6. 验证一个模型值不值RMSE与相关系数之外的事论文用了两个指标评估模型测试集RMSE为10.46预测曲线与真实曲线的相关系数为0.976。RMSE10.46的单位是μg/m³意味着平均每个样本的预测值与真实值偏差约10.46个单位。考虑到北京PM2.5浓度的年波动范围可能覆盖0到300以上10.46的绝对误差已经相当可观如果用它来判断某天是否达到75μg/m³的轻度污染线存在一定误判概率。所以单看RMSE还不够。相关系数0.976才是这篇论文最亮眼的地方。它说明预测曲线与真实曲线的整体变化趋势高度一致峰值出现的时间节点基本吻合。RMSE告诉你“平均偏了多少”相关系数告诉你“趋势跟得准不准”两个指标要配合着看只看任何一个都会被带偏。复现时我习惯再加一个指标看预测残差的时间分布而不是只看汇总数字。把每个测试样本的预测误差按季节分组统计如果冬天误差明显大于夏天说明模型对高浓度时段的拟合还不够可能是训练集中重污染天样本偏少。这种按条件拆分的验证方式比一个总RMSE更能指导下一步优化。还有一个实用的验证技巧把预测结果画成逐日对比曲线人工观察峰值滞后和谷底跟随情况。数字指标再好曲线上的“相位漂移”也骗不了人。论文里那个对比图就是干这个用的。从那以后我评估任何一个LSTM时间序列预测模型都强制走一遍这个流程先看RMSE量级是否可解释再算相关系数看趋势跟随度最后拆分组看残差。三步走完了才敢说这个模型是能落地的。希望帮到你。本文还有配套的精品资源点击获取