PM2.5时空长序列预测研究

### PM2.5时空长序列预测研究 #### 研究方法 ##### 研究区域与数据 研究区域选定为中国京津冀地区和长三角地区的184个城市，这些地区是中国空气污染最严重的区域。研究选取了2015 - 2018年的城市PM2.5数据以及气象数据，每3小时记录一次。其中，PM2.5数据来自国家控制站点，气象数据来自欧洲中期天气预报中心（ECMWF）的ERA5气候再分析大型数据集，该数据集由卫星数据反演并经过多源数据再分析校正。 ##### 问题定义 PM2.5预测属于典型的时空序列预测问题，与时间序列预测不同，它不仅要考虑数据时间的推移和连接，还要考虑节点间的空间相关性。接下来将PM2.5预测问题转化为数学公式： 假设在时间$t$共有$N$个城市节点，对于城市节点$i \in N$，其PM2.5浓度为$x_t^i$，其周边关联城市构成集合$M = \{l_1, l_2, l_3, \cdots, l_m\}$，$M \in N$。在时刻$t$，城市节点$i \in N$存在节点属性矩阵$P_t^i \in R^p$，其中$p$表示气象属性的数量。同时，通过经验模态分解（EMD）对城市节点的时间序列进行分解，得到$s$个子序列，那么每个节点的子序列在时刻$t$的值为$^1u_t^i, \cdots, ^su_t^i$，它们构成了EMD矩阵$U_t^i$。目标节点PM2.5预测问题可定义为以下形式： \[ \left[x_t^i, U_t^i, P_t^i, x_t^{l_1}, P_t^{l_1}, \cdots, x_t^{l_m}, P_t^{l_m}, x_{t + 1}^i, U_{t + 1}^i, P_{t + 1}^i, \cdots, x_{t + 1}^{l_m}, P_{t + 1}^{l_m}, \cdots, x_{t + \tau - 1}^{l_m}, P_{t + \tau - 1}^{l_m}\right] \to f(\cdot) \left[x_{t + \tau}^{l_i}, x_{t + \tau + 1}^{l_i}, \cdots, x_{t + T}^{l_I}\right] \] 其中$\tau \in [1, \cdots, T]$，$f(\cdot)$通过模型迭代$T$步得到： \[ f(\cdot) = g(\cdots g(g(\cdot))) \quad (T\text{次}) \] 训练时，使用均方根误差作为损失函数，假设预测值为$[\hat{x}_1^i, \cdots, \hat{x}_T^i]$，真实观测值为$[x_1^i, \cdots, x_T^i]$，则： \[ MSE_{loss} = \frac{1}{T} \sum_{t = 1}^{T} (\hat{x}_T^i - x_T^i)^2 \] ##### 图构建 由于PM2.5受多种因素影响，且这种影响通常具有空间相关性，为提高预测准确性，基于数据以及PM2.5自身的迁移和扩散规律构建图，将天气变量作为节点和边的属性，然后在图的基础上学习PM2.5的迁移和扩散过程，以获取所需的空间PM2.5污染相关性。 - **节点属性**：城市的气象条件会影响PM2.5污染的扩散传输，本文选取的气象属性如下表所示，这些气象属性可作为节点的特征。 | 变量 | 单位 | | --- | --- | | 边界层高度 | m | | K指数 | K | | 风速分量u | m/s | | 风速分量v | m/s | | 地表温度 | K | | 相对湿度 | % | | 总降水量 | m | | 地表压力 | Pa | - **邻接矩阵**：两个城市之间的距离以及城市间山脉的高度决定了它们之间PM2.5的关联程度，邻接矩阵计算公式如下： \[ A_{ij} = H(d_{\theta} - d_{ij}) \cdot H(m_{\theta} - m_{ij}) \] 其中$d_{ij} = \rho_i - \rho_j$，$m_{ij} = \sup_{\lambda \in (0, 1)} \left[h(\lambda\rho_i + (1 - \lambda)\rho_j) - \max\left(h(\rho_i), h(\rho_j)\right)\right]$，$\rho_i$表示节点$i$的位置（经纬度），$d_{ij}$表示两个节点之间的距离，$H(\cdot)$是阶跃函数，当且仅当$x > 0$时，$H(x) = 1$。$d_{\theta}$和$m_{\theta}$分别为距离和高度的阈值，此处分别为300 km和1200 m。在上述约束下，只有当两个城市节点之间的距离小于300 km且它们之间的山脉高度小于1200 m时，PM2.5才能进行传输和扩散。 ##### 模型结构 模型主要由三个模块组成： - **GAT模块**：基于2.3节构建的图提取PM2.5的空间信息。 - **EMD模块**：对原始PM2.5时间序列进行经验模态分解，得到多个子序列。 - **GRU模块**：整合GAT模块和EMD模块获得的空间和时间信息，在天气条件的影响下进行PM2.5预测。 下面详细解释这三个模块： ```mermaid graph LR classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px; classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; A([输入数据]):::startend --> B(GAT模块):::process B --> C(EMD模块):::process C --> D(GRU模块):::process D --> E([输出预测结果]):::startend ``` - **GAT模块**：由于注意力机制在各种深度学习模型中的成功应用，一些学者开始将注意力机制与图神经网络相结合。图注意力网络模型由velickovic等人首次提出，该模型在图神经网络的基础上引入注意力机制，大大提高了模型的预测效率。在PM2.5预测问题中，图注意力网络通过注意力机制计算图中每个节点的隐藏状态，关注节点的相邻节点，为相邻节点分配不同的权重信息，并根据权重信息聚合相邻节点的特征信息。具体步骤如下： 1. 定义图注意力层，在时间节点$t$，图注意力层的输入是$N$个节点的特征矩阵$P_t = [P_t^1, \cdots, P_t^i, \cdots, P_t^N]$，$P_t^i \in R^p$，$p$是节点特征的数量，通过图注意力层生成新的节点特征矩阵$P'_t = [P'_t^1, \cdots, P'_t^i, \cdots, P'_t^N]$，$P'_t^i \in R^{p'}$。 2. 输入PM2.5的城市网络图形结构并训练所有节点，得到图的权重矩阵$W \in R^{p \times p'}$，该矩阵用于将低维特征转换为高维特征。 3. 对每个节点执行自注意力机制$a: R^{p'} \times R^{p'} \to R$，通过自注意力机制计算节点之间的注意力系数： \[ e_{ij} = a(WP_t^i, WP_t^j) \] 4. 使用SoftMax函数对注意力系数进行归一化： \[ \alpha_{ij} = \text{softmax}(e_{ij}) = \frac{\exp(e_{ij})}{\sum_{k \in M} \exp(e_{ik})} \] 其中$M$是节点$i$的相邻节点集合。 5. 注意力机制$a(\cdot)$是一个单向前馈神经网络，其大小由权重向量$\vec{a} \in R^{2p'}$确定，并使用非线性激活函数LeakyReLU(·)激活，斜率取0.2，代入可得： \[ \alpha_{ij} = \frac{\exp(\text{LeakyReLU}(\vec{a}^T [WP_t^i \| WP_t^j]))}{\sum_{k \in M} \exp(\text{LeakyReLU}(\vec{a}^T [WP_t^i \| WP_t^k]))} \] 6. 每个节点的输出为： \[ P'_t^i = \sigma\left(\sum_{j \in M} \alpha_{ij}WP_t^i\right) \] 7. 为使预测结果更准确稳定，使用多头注意力机制，对$K$个独立的注意力机制进行并行操作，然后将结果拼接，最终空间特征输出结果为： \[ P'_t^i = \|_{k = 1}^{K} \sigma\left(\sum_{j \in M} \alpha_{ij}^k W^k P_t^i\right) \] 其中$\alpha_{ij}^k$表示第$k$个注意力机制的注意力系数归一化结果，最终输出特征矩阵$P'_t$有$Kp'$个特征。 - **EMD模块**：经验模态分解（EMD）作为一种自适应时间序列分解技术，主要思想是使用希尔伯特 - 黄变换（HHT）对非线性和非平稳时间序列数据进行滤波，直到得到平滑的序列数据。EMD的关键创新是引入了本征模态函数（IMF），IMF可定义为数据序列中的隐藏振荡模式，它可以是非平滑的，并且可以进行幅度或频率调制。IMF需满足以下两个条件： - 在整个数据序列中，极值点的数量（最大值和最小值的数量之和）与过零点的数量必须相等或最多相差1。 - 在序列的任何一点上，由局部最大值定义的包络和由局部最小值定义的包络的平均值为0。 基于以上定义，对单个节点$i$的PM2.5时间序列$x_i(t)$进行分解的步骤如下： 1. 确定序列$x_i(t)$的所有极值点，拟合上下极值点包络（上下包络应包含所有数据）。 2. 求上下包络的平均值$m_i(t)$，从原始值中减去平均值得到新序列$h_i(t)$，判断该序列是否满足IMF函数的定义。 3. 若满足则输出IMF，若不满足则对$h_i(t)$重复上述操作，直到满足IMF定义，此时$h_i(t)$即为所需的IMF函数$^ku_i(t)$。 4. 每次得到一个阶数的IMF后，将其从原始值中减去，直到最后剩余部分$r_n$是单调序列或常数序列。 经过上述方法分解后，单个节点的原始PM2.5时空序列$x_i(t)$可分解为一系列IMF函数和$r_n$的线性叠加： \[ x_i(t) = \sum_{k = 1}^{s} ^ku_i(t) + r_n(t) \] 其中$s$表示分解完成后IMF函数的数量。之后，将城市节点$i$在时刻$t$获得的IMF子序列拼接，得到城市节点$i$在时刻$t$的EMD矩阵$U_t^i$。 - **GRU模块**：门控循环单元（GRU）是一种用于处理时间序列的RNN模型，与传统RNN相比，其网络模型架构经过优化，缓解了梯度消失等问题。与LSTM相比，GRU所需的计算资源更少，效果也不逊色。GRU可以通过门控单元选择性地遗忘和保留过去时刻和当前时刻的信息，从而实现对未来时刻的预测。在时刻$t$，获取包含周边节点的空间特征信息$P'_t^i$和EMD矩阵$U_t^i$，将它们拼接得到时空GRU组件的输入。具体计算过程如下： \[ y_t^i = [P'_t^i, U_t^i] \] \[ z_t^i = \sigma(W_z \cdot [h_{t - 1}^i, y_t^i]) \] \[ r_t^i = \sigma(W_r \cdot [h_{t - 1}^i, y_t^i]) \] \[ \tilde{h}_t^i = \tanh(W \cdot [r_t^i \odot h_{t - 1}^i, y_t^i]) \] \[ h_t^i = (1 - z_t^i) \odot h_{t - 1}^i + z_t^i \odot \tilde{h}_t^i \] 其中$W_z$、$W_r$和$W$是在训练过程中可学习的权重矩阵，$y_t^i$是输入，$z_t^i$是更新门，$r_t^i$是重置门，$\tilde{h}_t^i$表示当前时刻的时空信息，通过更新门保留所需的时空信息，丢弃不需要的信息，得到最终的时空信息$h_t^i$。最后，通过全连接层得到最终预测值： \[ \hat{X}_t^i = \phi(h_t^i) \] 其中$\phi$表示全连接层。 ### PM2.5时空长序列预测研究 #### 实验部分 ##### 数据集 基于前文所述信息，构建了一个涵盖长达4年（2015年1月1日至2018年12月31日）、包含184个城市以及8个气象变量的数据集，数据粒度为每3小时记录一个数据点。为确定与PM2.5浓度变化相关的这些影响因素之间的相关性，首先对数据进行了相关性分析，并计算了每个气象变量与PM2.5序列之间的皮尔逊相关系数。计算结果如下表所示： | 变量 | 皮尔逊相关系数 | | --- | --- | | 边界层高度 | -0.34 | | K指数 | -0.64 | | 风速 | -0.53 | | 地表温度 | 0.49 | | 相对湿度 | 0.15 | | 总降水量 | -0.44 | | 地表压力 | 0.69 | 从表中可以看出，PM2.5序列与边界层高度、K指数、风速和年降雨量呈负相关，与地表温度、相对湿度和地表压力呈正相关。 随后，对数据集进行了预处理，包括两个部分： 1. **数据归一化**：在深度学习算法中，数据归一化非常必要。因为各个气象变量的量级不同，会对模型的收敛和最终数据的准确性产生较大影响。这里使用零均值归一化。 2. **数据集划分**：将经过零均值处理后的数据按照2:1:1的比例划分为训练集、验证集和测试集，以检验模型在一般远程设置下的预测能力。 ##### 实验设置 实验任务是基于起始的PM2.5浓度和未来72小时的气象变量数据，预测未来3、12、24、48、72和96小时的PM2.5浓度。为验证所提出模型的有效性，选择了以下模型进行比较： - **nodesFC - GRU**：这是所提出的GAT - EGRU模型的简化版本。通过将模型的前两个模块替换为全连接层来实现该模型的预测能力。目的是了解前两个模块所提取的时空信息对模型预测能力的实际提升程度。该模型也可视为类似CNN模型的替代方案。 - **GAT**：即图注意力网络。作为模型的简化版本，用于确定GRU和EMD模块对该模型的影响。所使用的GAT版本基于Pytorch - Geometric框架，在图数据处理和图神经网络操作方面具有优势。仅对最后一层网络进行了修改，以使模型适用于此类回归问题。 - **GC - LSTM**：该模型是基于图卷积神经网络（GCN）和长短期记忆网络（LSTM）的混合PM2.5时空预测模型。其中LSTM模块与本文的GRU模块差异不大，该模型旨在比较GAT在提取空间信息方面的优势。 - **HighAir**：这是一种基于分层图神经网络的空气质量预测方法。它考虑了空间质量影响因素对天气质量恶化的影响，并从分层角度构建图数据结构，使用LSTM作为解码器。选择该模块是为了测试EMD模块的性能，确保其在GAT - EGRU中的有效性。 **评价指标**：为更好地评估模型性能，选择了平均绝对误差（MAE）和均方误差（MSE）作为评价指标。MAE是观测值与真实值之间误差绝对值的平均值，用于描述预测值与真实值的误差；MSE是观测值与真实值之间差值平方的求和平均值，通常用于检测模型预测值与真实值之间的偏差，反映预测的准确性。这两个指标的值越小，模型的预测能力越好，越符合预期。 **超参数优化**：对模型的超参数进行贝叶斯优化，以获得模型的最优超参数值，并设计了获取函数以避免局部最优。贝叶斯优化基于高斯过程，通过不断添加新的采样点来更新先验分布并考虑先前的参数信息，因此比随机搜索和网格搜索更可靠、更高效。所提出的方法使用Adam优化器进行优化，批量大小设置为32，学习率设置为1e - 3，并在验证集上使用了早停策略。 **实验平台**：所有模型均在单个Tesla - V100 - SXM2 - 32GB GPU上运行。使用Pytorch实现模型，并使用Pytorch Geometric进行图形计算。数据开发平台是大连理工大学的大数据AIStation平台，本文的计算工作得到了大连理工大学大数据与智能决策研究中心的支持。 ##### 实验结果 - **成功构建图数据结构**：基于PM2.5和关联规则，得到了PM2.5的空间关联图，其中边表示城市节点之间的关联，城市节点的大小和颜色深浅表示污染的严重程度。选取了2018年11月23日18:00污染严重的时刻进行分析。从图中可以看出，大部分污染集中在京津冀和长三角地区，这些地区的PM2.5具有很强的相关性。西南地区的污染与其他地区的关联性较小，考虑到其地形，这种污染分布是合理的。同样，北方地区的污染扩散受到山脉的影响而受阻。污染区域的形成说明了相关性判断的重要性，因为无法从缺乏相关性的区域判断其他区域的污染状况。南北污染差异明显，东部中部地区的污染区域基本有污染中心，向外扩散形成一个或多个大小不一的污染区域；而南方地区各城市节点的污染程度相似，没有明显的污染中心，但由于交通便利，该地区的污染关联程度远大于北方。此外，图中还展示了沿海城市形成的污染关联。 ```mermaid graph LR classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px; classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px; A([获取PM2.5数据和关联规则]):::startend --> B(构建空间关联图):::process B --> C(分析污染分布特征):::process C --> D([得出相关性结论]):::startend ``` - **模型性能与评估**：在6个预测步骤和两个评价指标的12次比较中，GAT - EGRU模型获胜11次，这表明该模型在PM2.5预测方面优于其他模型。随着预测步长的增加，其他模型的误差开始大幅增加。当预测步长为96小时时，经过零均值处理的数据中，GC - LSTM和HighAir的MSE和MAE均上升到1以上，这意味着这两个模型在该步长下基本失去了预测能力，而GAT - EGRU模型仍保持在较低水平。以人们在现实世界中常关注的48小时预测步长为例，这两个模型的MSE分别为0.895和0.794，与GAT - EGRU模型的0.367相比，偏差过大，无法为人们的日常出行和生活提供有效的预警和指导。GAT - EGRU与nodesFC - GRU的得分比为12:1，从数值上看，这两个模型在两个评价指标上得分相近，但随着预测步长的增加，nodesFC - GRU的指标值开始大幅增大，例如从72小时到96小时，其MSE和MAE分别增加了0.221和0.205，而GAT - EGRU仅增加了0.198和0.182。这表明随着预测步长的增加，GAT - EGRU模型提取的时空信息在长期预测中发挥了作用，成功捕捉到了扩散过程的长期依赖性。同样，与GAT的比较也说明了EMD和GRU在捕捉时间依赖性方面的出色性能，这在PM2.5预测过程中是必不可少的能力。具体实验结果如下表所示： | 模型 | GAT - EGRU | nodesFC - GRU | GAT | GC - LSTM | HighAir | | --- | --- | --- | --- | --- | --- | | 指标(h) | MSE MAE | MSE MAE | MSE MAE | MSE MAE | MSE MAE | | 3 | 0.067 0.193 | 0.114 0.205 | 0.157 0.199 | 0.194 0.218 | 0.225 0.348 | | 12 | 0.315 0.442 | 0.353 0.349 | 0.370 0.397 | 0.445 0.517 | 0.578 0.696 | | 24 | 0.354 0.419 | 0.401 0.433 | 0.428 0.454 | 0.543 0.589 | 0.688 0.726 | | 48 | 0.367 0.425 | 0.474 0.452 | 0.496 0.517 | 0.895 0.648 | 0.794 0.819 | | 72 | 0.503 0.514 | 0.591 0.539 | 0.578 0.594 | 0.984 0.801 | 1.092 1.114 | | 96 | 0.701 0.696 | 0.812 0.744 | 0.766 0.751 | 1.088 0.991 | 1.605 1.828 | - **模型细粒度分析与真实数据比较**：作为模型性能评估的补充，将模型的预测结果与测试集中某些真实的PM2.5浓度进行了比较。目的是展示GAT - EGRU模型在数值上表现良好的原因以及其在哪些方面表现出色。选择了大连作为目标城市，比较预测值和真实值。由于大连的PM2.5时间和数值分布相对均匀，能够很好地反映模型在预测各个数值段的性能。从图中可以大致看出，PM2.5的预测趋势与观测值一致，验证了模型捕捉空间和时间变化的可行性。对于极值点，模型往往会低估或高估数值，这在一定程度上是合理的，因为模型没有专门处理异常值的模块。相反，模型在预测位于[-1, 1]数值区间的趋势时更为准确。总体而言，模型对PM2.5浓度值做出了正确的预测。 - **鲁棒性实验结果分析**：为验证模型各模块的有效性，进行了蒸馏实验，旨在验证模型的鲁棒性并确定各模块的作用。蒸馏实验还可以验证气象变量对GAT - EGRU模型性能的影响程度，从而确定此类领域知识在时空长序列预测中的重要性。实验结果如下表所示： | 模型 | GAT - EGRU | GAT - GRU | NodesFC - EGRU | GAT - EMD | GAT - EGRU - O | | --- | --- | --- | --- | --- | --- | | 指标(h) | MSE MAE | MSE MAE | MSE MAE | MSE MAE | MSE MAE | | 48 | 0.367 0.425 | 0.432 0.460 | 0.533 0.481 | 0.662 0.698 | 0.419 0.437 | | 72 | 0.503 0.514 | 0.569 0.581 | 0.672 0.588 | 0.981 1.113 | 0.563 0.571 | | 96 | 0.701 0.696 | 0.793 0.770 | 0.795 0.782 | 1.347 1.568 | 0.778 0.726 | 在进行实验结果分析之前，先对鲁棒性实验中的模型进行说明： - **GAT - GRU**：从模型中移除EMD模块，仅使用原始的PM2.5序列与空间信息拼接，用于验证EMD模块在时间信息提取方面的能力。 - **NodesFC - EGRU**：与NodesFC - GRU不同，保留了EMD模块，与原始模型相比仅移除了GAT模块。这意味着模型的空间特征提取部分被完全消除，仅剩下经过EMD模块处理后输入到GRU的信息，用于验证空间特征提取在模型中的有效性。 从实验结果可以看出，移除相应模块后，模型的性能有所下降，这表明EMD模块在提取时间信息和GAT模块在提取空间信息方面都发挥了重要作用，进一步验证了GAT - EGRU模型各模块的有效性和合理性。