CN116307068A

CN116307068A - 基于四维有向gcn-lstm模型的多城市多种大气污染物预测方法

Info

Publication number: CN116307068A
Application number: CN202310035696.0A
Authority: CN
Inventors: 王立; 唐千慧; 王小艺; 许继平; 赵峙尧; 张慧妍; 于家斌; 孙茜; 白玉廷; 王昭洋; 崔晓玉
Original assignee: Beijing Technology and Business University
Current assignee: Beijing Technology and Business University
Priority date: 2023-01-10
Filing date: 2023-01-10
Publication date: 2023-06-23

Abstract

本发明是一种基于四维有向GCN‑LSTM模型的多城市多种大气污染物预测方法，解决目前不能准确预测多种大气污染物浓度的问题。本发明方法包括：对多城市多种大气污染物浓度数据进行相关性分析，确定数据间具有相关性；根据城市的地理位置建立四维有向GCN模型的图；进行四维有向GCN模型的谱分解及张量运算，获取图傅里叶系数和图傅里叶基；改进四维有向GCN模型的图滤波器；引入LSTM网络架构，构建四维有向GCN‑LSTM模型进行大气污染预测。本发明构建了一种能够包含二维空间有向性信息、捕捉图节点内多个因素一维相关性信息、提取不同时刻时间一维相关性信息的四维有向图模型，提高了对多城市多种大气污染物的预测精度。

Description

基于四维有向GCN-LSTM模型的多城市多种大气污染物预测方法

技术领域

本发明涉及一种大气污染物浓度预测方法，属于人工智能技术与环境学科交叉领域，具体涉及一种四维有向GCN-LSTM模型的多城市多种大气污染物浓度大数据预测方法。

背景技术

近些年，伴随全国工业化和城市化水平越来越高，能源消耗也急剧增加，以煤炭、石油和天然气为主的一次能源消耗占比高达84.7％，这些化石能源的大量使用直接致使我国经济发达地区空气质量较差，大气污染物浓度指数超标。根据国家颁布的《环境空气质量标准》，明确规定了六个基本的大气污染物，包括了大气颗粒物和气态污染物：PM2.5、PM10、NO₂、SO₂、CO、O₃。这些大气污染物被世界卫生组织评定为世界最大的环境健康风险诱因。据统计，近十年来，我国多个省市出现不同程度的极端低能见度和重度空气污染事件，大气污染物种类也由工业烟尘转为PM2.5、PM10这类直径更细小的颗粒物。大气污染除了会直接对人体健康造成危害和影响人类社会活动，还会给自然生态系统带来负面影响。如：SO₂致使农耕土壤酸化森林衰亡、PM2.5、PM10致使大气能见度降低、O₃大量分解全球气候变暖、大气污染物携带氮磷等元素通过干湿沉降致使湖泊河流富营养化。因此对大气污染物浓度预测非常有必要。

目前，大气污染物浓度预测建模方法主要分为两类：基于机理驱动的大气污染预测建模方法和基于数理统计的大气污染预测建模方法。机理驱动的大气污染预测建模方法是以大气动力学、大气环境化学为基础，根据大气污染排放源数据、气象数据，用方程组构建数学模型来分析污染物时空分布和扩散迁移，再计算求解实现预测未来污染物的浓度变化和分布。由于大气污染物的产生和浓度大小受污染源的空间地理位置、区域气象条件等诸多因素影响，使得系统结构难以精准建立，因此仅仅依赖机理驱动的建模方法并不适合大气污染物浓度的预测。大气污染物数理统计预测模型是建立在数据分析与数据挖掘理论之上，此类模型不依赖污染物产生和扩散机理，不依赖相关的物理、化学作用及生物过程，以统计学为基础，通过大气污染物及相关影响因素数据，分析污染物浓度变化规律，给出未来一段时间内大气污染物浓度的预测结果。相对于机理驱动方法，数理统计的预测成本较低，尤其是多频次的短期预测中具备极大优势。而在现有的大气污染物浓度预测方法中，大多采用单一的数据驱动模型，如回归模型、灰色理论模型、进化算法、机器学习等，它们都存在预测结果不够准确的弊端。

多城市多种大气污染物浓度数据具有海量、高维、时空演变、多因素等特征，符合大数据特点，对其进行预测应采用适用于大数据分析的深度学习方法。传统的循环神经网络和长短期记忆网络在处理时间序列数据时，任务中的数据样本是由多个因素的时间序列值组成，捕捉时间序列中元素之间的短记忆和长记忆依赖关系，使得序列中排序靠后的元素为排序靠前的元素的计算提供指导信息。然而，多个城市大气污染物的扩散运动受动力因素和热力因素影响，具有空间相关性。因此，研究城市的大气污染物浓度则需要考虑其周围城市的大气污染物浓度之间的时间-空间相关性，采用一种适用于研究多因素的时空序列建模方法。

图卷积神经网络是以图论为基础的一种利用深度学习对图论中的图数据搭建起的建模框架，图数据包括图的节点和边，主要描述了节点之间的两两关系，可将各城市作为图的节点，城市间关联作为图的边，大气污染物即大气污染的表征因素，作为图节点的属性。传统的图卷积神经网络在处理时空序列预测问题时，一般仅考虑建立一维空间距离图来表示多城市多种大气污染物在空间上的相关性，且大多数研究仅采用无向图来反映多种大气污染物在多个城市间的扩散运动影响。但是多城市间的大气污染物实际包含经度和纬度的二维空间关系，并且大气污染物扩散具有一定的有向性。应采用包含二维空间信息的有向图描述多城市间的大气污染物空间相关性，同时，需要考虑多种大气污染物间的多因素相关性以及时间相关性。

发明内容

针对上述现有技术在预测城市多种大气污染物浓度数据不够准确的问题，本发明提供了一种基于四维有向GCN-LSTM模型的多城市多种大气污染物预测方法，解决了传统图卷积神经网络在预测多城市大气污染物浓度时，建立的图模型未考虑多城市间的二维空间关系和大气污染物扩散有向性的问题，通过构建一种能够包含二维空间有向性信息、捕捉图节点内多个因素一维相关性信息、提取不同时刻时间一维相关性信息的四维有向图模型，为具有时空特性的多城市多种大气污染物浓度预测提供一种新思路。

本发明提供的基于四维有向GCN-LSTM模型的多城市多种大气污染物预测方法，包括以下五个步骤：

步骤一，对多城市的多种大气污染物浓度数据进行相关性分析，确定不同城市的多种大气污染物浓度数据间具有相关性；

步骤二，根据所要研究的城市的地理位置建立四维有向GCN(图卷积神经网络)模型的图；

该图为二维空间有向图，图中将每座城市作为一个节点，节点位置为城市中心的经纬度二维坐标，图中节点间都存在有向边，有向边表示为一节点指向另一节点的平面向量；将每种大气污染物作为节点的一个属性，属性值随时间变化；构建二维空间有向图的邻接矩阵、度矩阵和拉普拉斯矩阵；其中邻接矩阵用于描述图中节点和节点间边的信息，度矩阵用于描述图中节点与相连节点间边的权重之和；可以将多城市多种大气污染物浓度时间序列表示为上述二维空间有向图信息序列；

步骤三，进行四维有向GCN模型的谱分解及张量运算；

四维有向GCN模型的谱分解是指：由于二维空间有向图的拉普拉斯矩阵为反对称矩阵，不能做频域卷积，因此将二维空间有向图的邻接矩阵通过线性代数变换构建一个新的邻接矩阵，新的邻接矩阵为对称矩阵，包含旧的度矩阵和旧的邻接矩阵的信息；根据新的邻接矩阵更新度矩阵，获取新的对称的拉普拉斯矩阵，对新的拉普拉斯矩阵进行谱分解，获取特征值矩阵和特征向量矩阵，分别作为图傅里叶系数和图傅里叶基，能够进行谱域上的图卷积运算；

四维有向GCN模型的张量运算是指：传统GCN模型的输入数据为二维矩阵，经谱分解获取的特征值矩阵和特征向量矩阵均为三维矩阵，存在张量维数不匹配的问题，因此，采用广播机制将输入矩阵进行维数扩展和元素复制，变换为与所述特征向量矩阵维数相同的三维矩阵，将变换后的输入矩阵与特征向量矩阵进行张量运算，进行二维空间有向性的信息整合；所述输入矩阵是根据多城市多种大气污染物浓度数据获取的矩阵，在下面步骤五中说明；

步骤四，改进四维有向GCN模型的图滤波器；

图信号在进行傅里叶变换转到频域后，得到图傅里叶基和图傅里叶系数，通过图滤波器对图傅里叶系数进行调制，传统的图滤波器采用切比雪夫滤波器，但在通带频率的末端部分会存在较强的波动，导致信号处理不稳定，因此本发明利用Hermite滤波器作为图滤波器，对变换后的输入矩阵进行滤波，以解决末端部分信号出现受噪声干扰的问题；

步骤五，构建四维有向GCN-LSTM模型进行大气污染预测；

所述的四维有向GCN-LSTM模型是将四维有向GCN模型嵌入到LSTM(长短期记忆网络)网络架构里，该模型首先将当前时刻的多城市多种大气污染物浓度X_t与四维有向GCN模型上一时刻的输出数据

一起输入LSTM模型，利用LSTM模型提取不同时刻间多因素的时序特征，输出四维有向GCN模型的输入矩阵H_t，将输入矩阵H_t采用广播机制变换为与特征向量矩阵维数相同的三维矩阵后，输入四维有向GCN模型；四维有向GCN模型利用图滤波器对输入矩阵H_t进行滤波后，输入网络层进行卷积操作，输入矩阵H_t与网络层的输出一起送入输入门、遗忘门和输出门，最后四维有向GCN模型输出t时刻的输出数据/>

输出数据/>

记录下一时刻的各城市每种大气污染物的浓度预测值。训练所述四维有向GCN-LSTM模型，利用训练好的模型进行大气污染预测。

相比现有技术，本发明的优点在于：

(1)本发明方法首先采用典型相关分析方法分析并验证了具有一定空间地理位置关联的不同城市不同时刻多种大气污染物浓度存在相关性，说明了多种大气污染物浓度不仅在时间维度上具有相关性而且在空间维度上也具有一定相关性，为后续预测模型的提出奠定基础。

(2)本发明方法通过建立二维空间有向图来拓扑多城市间大气污染物的空间相关性，将传统GCN模型的图由传统包含一维空间距离信息的三维无向图改为包含二维空间有向信息的四维有向图，通过此方式构建的图邻接矩阵和度矩阵能够充分挖掘城市间的二维空间有向信息，从而提出能够全面提取二维空间相关性、时间相关性、多因素相关性特征的四维有向GCN模型，提升GCN模型特征提取的有效性，提高模型预测精度。

(3)本发明方法提出四维有向GCN模型的谱分解和张量运算算法，利用反对称矩阵性质，将有向图的反对称的图邻接矩阵和对称的图度矩阵通过线性代数变换成对称的图邻接矩阵，解决了无法对有向图的拉普拉斯矩阵进行谱分解的问题，并采用广播机制对四维有向GCN模型的输入信号和有向图拉普拉斯矩阵的特征向量组成的矩阵进行张量运算，解决了四维有向GCN模型的输入信号与二维空间信息的图拉普拉斯矩阵的特征向量矩阵存在维度不匹配的问题。

(4)本发明方法对四维有向GCN模型的滤波器进行了改进，传统GCN模型中使用的是切比雪夫滤波器，即Chebyshev多项式滤波器，而Chebyshev多项式在拟合逼近信号时对通带频率末端这部分信号的处理过程容易受噪声干扰不稳定，因此本发明使用Hermite多项式滤波器替换Chebyshev多项式滤波器，Hermite滤波器具有良好的带通特性可以有效地减少频率末端这部分信号受噪声干扰的现象出现，提高模型运算的准确性和速度。

(5)本发明方法提出将四维有向GCN模型嵌入LSTM网络构成四维有向GCN-LSTM模型，利用四维有向GCN模型具有刻画空间分布能力来提取不同城市的多种大气污染物浓度的空间相关性，同时再通过LSTM网络循环递归来捕捉记录多种大气污染物浓度变化的多因素相关性和时序相关性，有效解决多种大气污染物浓度这类多因素时空序列数据进行多因素相关性、时间相关性和空间相关性的统一建模的问题，经试验验证，采用本发明方法的预测模型具有更高的预测准确度。

附图说明

图1是本发明的多城市多种大气污染物预测方法流程示意图；

图2是本发明方法构建的四维有向图的示意图；

图3是本发明四维有向图中的二维空间有向图的示意图；

图4是本发明四维有向GCN模型(4D-DGCN)的结构图；

图5是本发明四维有向GCN-LSTM模型(4D-DGCN-LSTM)的结构图；

图6是本发明实施例使用典型相关分析方法对典型变量x_r和X组变量分析所得载荷系数三维可视曲线图；

图7是本发明实施例使用典型相关分析方法对典型变量x_r和X组变量分析所得载荷系数三维可视曲线图的左视图；

图8是本发明实施例使用典型相关分析方法对典型变量y_r和Y组变量分析所得载荷系数三维可视曲线图；

图9是本发明实施例使用典型相关分析方法对典型变量y_r和Y组变量分析所得载荷系数三维可视曲线图的左视图；

图10是本发明实施例中使用本发明方法、GCN及LSTM方法对苏州市PM2.5浓度的预测值和真实值曲线图；

图11是本发明实施例中使用三种方法对苏州市PM2.5浓度在第21～40个时间节点的预测值和真实值曲线图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明方法提出一种四维有向图卷积神经网络模型，基于该模型实现了多城市多种大气污染物时空大数据预测方法。本发明在构建四维有向图卷积神经网络模型中，还提出相应的二维空间有向图谱分解和张量计算算法，在此基础上，更新四维有向图卷积神经网络的图滤波器以便更好地调制输入信号的频率，改进四维有向图卷积神经网络和长短期记忆网络联结方式，以便实现同步提取多因素数据的时间和空间特征，为多城市多种大气污染物的更加准确的预测提供一种新方案。

如图1所示，本发明实现的基于四维有向GCN-LSTM模型的多城市多种大气污染物预测方法，通过下面五个步骤来说明。

步骤一、基于典型相关分析方法对大气污染物时空数据的相关性进行分析。

大气污染物浓度数据不仅在时间维度上具有一定相关性，同时，不同种类的大气污染物PM2.5、PM10、NO₂、SO₂、CO、O₃之间、不同城市间的大气污染物浓度数据也具有一定相关性，因此，多城市多种大气污染物的浓度数据是一种多因素相关且时间-空间相关的序列。因此，本步骤使用适用于分析和量化两组变量间相关性大小问题的典型相关分析方法，来衡量多个城市的多种大气污染物待预测数值和多个城市的多种大气污染物历史数值间相关性强弱。

选取N个城市M种大气污染物当前时刻t和p个历史时刻数值作为X组数据集，X组变量个数共N*M*(p+1)个。选取N个城市M种大气污染物的t+1时刻数值作为Y组，Y组变量个数共N*M个，构成的两组变量如公式(1)(2)所示：

X＝(Var₁₁(t),…,Var_NM(t),…,Var_ij(t-1),…,Var₁₁(t-p),…,Var_NM(t-p)) (1)

Y＝(Var₁₁(t+1),…,Var_ij(t+1),…,Var_NM(t+1)) (2)

其中，t,t-1,…,t-p表示当前时刻t和p个历史时刻，t+1表示下一时刻，i＝1,2,…,N，j＝1,2,…,M，Var表示大气污染物浓度，Var_ij(t-1)表示t-1时刻第i个城市第j种大气污染物浓度。

对X和Y这两组变量进行线性组合得到第一对典型变量x₁和y₁，具体如公式(3)所示：

其中{a0₁₁,…,a1_ij,…ap_NM}表示X组的典型系数，{b₁₁,…,b_ij…,b_NM}表示Y组的典型系数。

利用Pearson相关系数计算公式得出第一对典型变量x₁和y₁的典型相关系数P₁如下：

其中，Cov为求取协方差，Var为求取标准差。

重复公式(3)、(4)直至进行到r步，X和Y两组变量的相关性被提取完，其中r＝min(X组变量个数,Y组变量个数)。得到r对典型变量，分析这r对典型变量和X、Y两组变量的典型载荷系数和交叉载荷系数，典型载荷系数表征典型变量x_r和X组变量、典型变量y_r和Y组变量的简单相关系数，交叉载荷系数表征典型变量x_r和Y组变量、典型变量y_r和X组变量的简单相关系数，载荷系数绝对值越大说明该项与典型变量之间的相关关系越强。

进行相关性分析时，可以看典型载荷系数也可以看交叉载荷系数。例如，x₁和X产生的N*M*(p+1)个典型载荷系数中，第1个和第N个系数绝对值最高；y₁和Y产生的N*M个典型载荷系数中，第2个和第M个系数绝对值最高；那么就认为X组里的第1个和第N个变量和Y组里的第2个和第M个变量存在很强的相关性。

经过试验数据验证，可以确定N个城市M种大气污染物在待预测时刻的浓度和N个城市M种大气污染物的p个历史时刻浓度值存在相关性，即多城市多种大气污染物的浓度数据具有时空相关性。

步骤二、根据所要研究的城市的地理位置建立四维有向GCN(图卷积神经网络)模型的图。

由图2所示，四维有向GCN模型的图所包含的空间关系图是一个二维空间有向图，因此对四维有向GCN模型的图定义关键是针对二维空间有向图的定义。如图2所示，将多城市多种大气污染物的时空序列构建成了图数据形式。本发明实施例中根据经度和纬度建立一个二维平面坐标系，如图3所示。

(2.1)构建二维空间有向图，如下：

定义二维空间有向图为G＝(v,ε)，其中：

ν＝{v₁:"城市1",v₂:"城市2",…,v_N:"城市N"}为图的节点集合，v_i表示第i个城市节点，i＝1,2,…,N；节点的位置用城市中心的地理经度和维度的二维空间信息表示；

为节点之间边的集合，/>

表示节点v_i指向节点v_j的平面向量，i,j＝1,2,…,N，如图3所示。设节点v_i、v_j的经纬度坐标分别为(x_vi,y_vi)、(x_vj,y_vj)，则/>

所有节点间都存在有向边。

(2.2)构建二维空间有向图邻接矩阵，如下：

给定二维空间有向图G＝(v,ε)，定义邻接矩阵A来描述该图的节点和节点间边的信息，邻接矩阵A的第i行第j列的元素A_ij表示节点v_i和v_j的连接关系，i,j＝1,2,…,N。考虑城市大气污染浓度变化受周围城市大气污染物扩散影响，即不同城市节点间的空间地理位置关系影响城市大气污染物浓度变化，因此建立以经纬度坐标构成的平面向量来描述节点间边信息的邻接矩阵，每个城市节点都具有唯一的经纬度坐标，两两城市间边关系为对应经纬度坐标之差构成的平面向量，具体邻接矩阵A如下：

其中x_vi和y_vi表示城市i的纬度坐标和经度坐标。

(2.3)构建二维空间有向图的度矩阵，如下：

定义度矩阵D来描述图G＝(ν,ε)中节点与相连节点间边的权重之和，与邻接矩阵一一对应。度矩阵D是对角矩阵，即除了第i行第i列的元素D_ii表示与节点i相关的边的权重之和，其余位置元素均为0，i＝1,2,…,N。度矩阵D具体如下：

其中：

(2.4)构建二维空间有向图的拉普拉斯矩阵。

定义拉普拉斯矩阵L为度矩阵D与邻接矩阵A之差。拉普拉斯矩阵L计算公式如下：

步骤三、四维有向GCN模型的谱分解及张量运算。

(3.1)四维有向GCN模型的谱分解。

传统的基于频域图卷积神经网络中，对图的拉普拉斯矩阵做频域卷积条件是图为无向图，即图的度矩阵和邻接矩阵均为对称矩阵。图的拉普拉斯矩阵L进行谱分解，需要满足前提该矩阵为对称矩阵，考虑式(5)中度矩阵D为对角矩阵即是对称矩阵，因此若拉普拉斯矩阵L为对称矩阵则需要邻接矩阵A一定为对称矩阵。无向图的邻接矩阵A为对称矩阵，而有向图的邻接矩阵A为反对称矩阵。因此，有向图无法直接对其拉普拉斯矩阵L进行谱分解来获得特征值Λ和特征向量U，进而无法直接对有向图进行图卷积处理。

根据反对称矩阵性质3，若D为对称矩阵，A为反对称矩阵，则存在如下关系：

D＝D^Τ,A＝-A^Τ

(DA-AD)＝(DA-AD)^Τ (7)

上角标T表示矩阵转置。

即由A、D构成的新邻接矩阵A′＝DA-AD为对称矩阵，则更新度矩阵D′如下：

更新拉普拉斯矩阵L′如下：

L′＝D′-A′＝D′-(A-AD) (8)

更新后的拉普拉斯矩阵L′为对称矩阵可进行谱分解，谱分解后得到特征值矩阵Λ和特征向量矩阵U如下：

其中，

表示矩阵L′的第i个特征值，μ_i表示第i个特征值对应的特征向量，

表示第j个特征向量内的第i个元素。

(3.2)四维有向GCN模型的张量运算。

传统GCN模型的输入数据为二维矩阵H，如下：

其中，矩阵H为N*M的二维矩阵，h_ij表示第i个城市第j种大气污染物对传统GCN模型的输入数据，i＝1,2,…,N，j＝1,2,…,M。传统输入数据H是大气污染物浓度形成的矩阵，本发明对其进行改进，具体在下面步骤五中说明。

对于四维有向GCN模型而言，其矩阵Λ和U均为N*N*2的三维矩阵，在矩阵Λ和U与矩阵H进行张量运算时存在张量维数不匹配的问题，因此使用Numerical Python库中广播机制对H进行维数扩展和元素复制来解决该问题，将H变换为四维有向GCN模型的输入数据，如下：

变换后的矩阵H为N*M*2的三维矩阵，(h_NM,h_NM)^Τ表示元素h_NM在第三维度上的广播拉伸并进行一次转置。变换后的H为三维矩阵，与U进行点乘的张量运算的公式如下：

步骤四、改进四维有向GCN模型的图滤波器。

图信号在进行傅里叶变换转到频域后，得到图傅里叶基和图傅里叶系数，通过图滤波器对图傅里叶系数进行调制，使得其中一些频率分量被保留或放大，另外一些频率分量被移除或减小。传统的图滤波器采用的是切比雪夫滤波器，可以实现信号的快速衰减但在通带频率的末端部分会存在较强的波动，导致信号处理不稳定。本发明利用Hermite滤波器具有良好的带通特性可以有效地去除干扰噪声，从而提高了频率计算的准确性和速度。

定义图滤波器多项式为g_θ(Λ)，具体如下：

其中θ_κ表示滤波器的系数，K表示图滤波器多项式截断参数，用来调制图输入数据H，通过对数据H和图滤波器g_θ(Λ)进行卷积操作来实现特征变换，具体公式如下：

然而图滤波器多项式g_θ(Λ)的基{Λ⁰,Λ¹,Λ²,…,Λ^K}彼此为非正交，即系数之间是相互依赖的，模型训练时各部分基容易相互干扰。因此，采用具有正交性的Hermite多项式来更新图滤波器，Hermite多项式递归关系如下：

Her_k(Λ)＝2Λ·Her_k-1(Λ)-2(k-1)Her_k-2(Λ) (13)

其中k≥2，Her₀(Λ)＝Λ⁰，Her₁(Λ)＝2Λ。

将式(13)带入式(11)中得：

将式(14)代入式(12)中，对数据H使用Hermite多项式滤波器，具体计算如公式(15)所示：

其中，ST表示Hermite多项式滤波器的输出数据。

步骤五、构建四维有向GCN-LSTM模型，用于大气污染预测。

多城市多种大气污染物浓度数据不仅具有空间相关性，即二维空间有向性，同时也具有时间相关性和多因素相关性，而四维有向GCN模型通过二维空间有向图结构可以拓扑出多因素数据间的空间相关性，还需要对数据间的多因素相关性和时间相关性进行提取。考虑长短期记忆网络LSTM可以有效捕捉数据的多因素时序特征，因此将四维有向GCN模型嵌入到长短期记忆网络LSTM中构成四维有向GCN-LSTM模型，实现在提取数据空间性的同时能够捕捉其多因素时序性，提高模型预测精度。

(5.1)首先说明四维有向GCN模型的计算公式。

多城市多种大气污染物浓度矩阵经LSTM模型输出四维有向GCN模型的输入矩阵，对各时间点输入矩阵H_t经步骤三维度变换和步骤四滤波，获取滤波后数据，设滤波后时刻t的输入数据为ST_t。将Hermite多项式滤波器的输出数据ST_t传入四维有向GCN模型网络层进行权重和偏置的训练，其中四维有向GCN模型网络结构如图4所示，四维有向GCN模型的计算公式如下：

ST′_t＝Sigmoid(ST_t*W_G+b_G) (17)

i′_t＝Sigmoid(W′_i*[ST′_t,H_t]+b′_i) (18)

g′_t＝tanh(W′_g*[ST′_t,H_t]+b′_g) (19)

G_t＝(1+i′_t)·G_t-1+i′_t·g′_t (20)

o′_t＝Sigmoid(W′_o·[ST′_t,H_t]+b′_o) (21)

其中H_t表示t时刻四维有向GCN模型的输入数据，ST_t表示t时刻Hermite多项式滤波器的输出数据，W_G、b_G分别表示四维有向GCN模型中卷积操作的权重和偏置，ST′_t表示t时刻四维有向GCN模型网络层训练输出，W′_i、b′_i分别表示四维有向GCN模型中输入门的权重和偏置，i′_t表示t时刻四维有向GCN模型中输入门的输出，W′_g、b′_g分别表示四维有向GCN模型中遗忘门的权重和偏置，g′_t表示t时刻四维有向GCN模型中遗忘门的输出，W′_o、b′_o分别表示四维有向GCN模型中输出门的权重和偏置，O′_t表示t时刻四维有向GCN模型中输出门的输出，G_t-1表示t-1时刻空间特征的短期记忆和长期记忆整合的输出，G_t表示t时刻空间特征的短期记忆和长期记忆整合的输出，

表示t时刻四维有向GCN模型的输出数据。

如图4所示，四维有向GCN模型的结构和实现功能主要包括：将四维有向GCN模型的二维空间有向图分别在经度和维度提取信息，获取邻接矩阵、度矩阵和拉普拉斯矩阵，对拉普拉斯矩阵进行谱分解，获得特征值矩阵和特征向量矩阵，分别对应为图傅里叶系数和图傅里叶基，图滤波器对图傅里叶系数进行调制，对输入数据H_t进行滤波，滤波后数据ST_t经四维有向GCN模型的网络层训练输出ST′_t，与输入数据H_t一起送入输入门、遗忘门和输出门，结合输入门的输出i′_t与遗忘门的输出g′_t获取t时刻空间特征的短期记忆和长期记忆整合的输出G_t，由输出门的输出O′_t与G_t得到t时刻四维有向GCN模型的输出数据

(5.2)构建四维有向GCN-LSTM模型。

再将四维有向GCN模型嵌入到LSTM网络架构里，利用LSTM的长短期记忆循环结构提取不同时刻间多因素的时序特征，即在对图数据进行空间特征提取的同时保留多因素数据的时序关系，进而实现对多因素数据的空间相关性和时间相关性的同时建模，其中同步提取多因素时空特征的四维有向GCN-LSTM模型具体结构如图5所示，四维有向GCN-LSTM模型中LSTM结构图的计算公式如下：

C_t＝f_t·C_t-1+i_t·g_t (27)

H_t＝O_t·_tanh(C_t) (28)

其中

分别表示t-1时刻四维有向GCN-LSTM模型的输出，X_t表示t时刻大气污染物浓度矩阵，f_t表示t时刻LSTM网络中遗忘门输出，W_f、b_f分别表示LSTM网络中遗忘门的权重和偏置，i_t表示t时刻LSTM网络中输入门输出，W_i和b_i分别表示LSTM网络中输入门、遗忘门的权重和偏置，g_t表示一个中间写入信息，W_g和b_g分别表示中间写入信息网络的权重和偏置，o_t表示t时刻LSTM网络中输出门输出，W_o、b_o分别表示LSTM网络中输出门的权重和偏置，C_t-1和C_t分别表示LSTM网络一个记忆单元的t-1时刻状态和t时刻状态，Sigmoid和tanh分别为模型网络训练使用的两种激活函数。

如图5所示，获取多城市多种大气污染物浓度的时间序列，将t时刻的浓度数据X_t与四维有向GCN模型上一时刻的输出数据

一起输入LSTM模型，然后LSTM模型输出H_t上给四维有向GCN模型。

(5.3)利用构建好的四维有向GCN-LSTM模型进行大气污染预测。

预测多城市多种大气污染物浓度的工作过程如下：

使用构建好的四维有向GCN-LSTM模型对t+1时刻的大气污染物浓度进行预测，首先确定模型的输入为{X_t-p,…,X_t-1,X_t}是由N个城市M种大气污染物的当前时刻t和p个历史时刻浓度数值组成，具体表示为：

/>

其中，Var_ij(t)表示t时刻第i个城市第j个大气污染物浓度，i＝1,2,…,N，j＝1,2,…,M，使用四维有向GCN-LSTM模型来预测t+1时刻的大气污染物浓度数值，则t时刻四维有向GCN-LSTM模型的输出

是N个城市M种大气污染物的t+1时刻浓度预测值，即：

实现四维有向GCN-LSTM模型对多城市多种大气污染物浓度的预测。

实施例：

以江苏省苏州市及环苏州市城市群的六项大气污染物PM2.5、PM10、NO₂、SO₂、CO、O₃数据为例，采用本发明方法进行六项大气污染物的预测，其中城市群包括：苏州、南通、常州、宜兴、湖州、嘉兴和上海市。以2009年城市气象监测数据为例，经过数据筛选处理后，共选取使用的数据集包括七个城市的六个大气参数，时间跨度为2009年1月1日至2009年12月31日，共372天，每天采集8次数据，数据纬度(2976,7,6)，分别表示时间长度为2976，图节点数为7个，对应7座城市，数据特征为6个，对应6种大气污染物浓度。实验过程选取前330天数据作为模型训练集，余下天数中在选取前30天作为测试集。具体获取的数据集如表1所示，采用本发明方法处理如下步骤一至五。

表1数据集

属性	数值	内容
			节点个数	7个	苏州、南通、常州、宜兴、湖州、嘉兴、上海
节点特征	6个	PM2.5、PM10、CO、O₃、NO₂、SO₂
			时序长度	2009年	372天×8个/天＝2972个时间节点
数据分割	(330,30)	第1-330天为训练集，第331-360天为测试集

步骤一、基于典型相关分析方法进行大气污染物时空数据相关性分析。

首先提取典型变量对，计算成对典型变量之间的典型相关系数，使得具有最大相关性；其次分别计算集合X和集合Y的典型变量的系数，得到典型变量的组成公式以及具体典型相关系数。如表2所示，为本发明实施例计算的典型相关系数。

表2 X组和Y组的典型相关系数

典型变量对	典型相关系数	典型变量对	典型相关系数	典型变量对	典型相关系数
						第1对	0.973	第15对	0.809	第29对	0.674
第2对	0.963	第16对	0.798	第30对	0.659
						第3对	0.932	第17对	0.796	第31对	0.654
第4对	0.927	第18对	0.786	第32对	0.647
						第5对	0.905	第19对	0.762	第33对	0.634
第6对	0.9	第20对	0.761	第34对	0.623
						第7对	0.889	第21对	0.756	第35对	0.608
第8对	0.88	第22对	0.731	第36对	0.601
						第9对	0.871	第23对	0.715	第37对	0.587
第10对	0.854	第24对	0.709	第38对	0.57
						第11对	0.847	第25对	0.705	第39对	0.563
第12对	0.833	第26对	0.701	第40对	0.538
						第13对	0.829	第27对	0.688	第41对	0.508
第14对	0.813	第28对	0.676	第42对	0.48

根据典型载荷系数进行典型结果分析，由典型载荷系数得到典型变量与本组所有变量的相关性，即典型变量x_r和X组变量、典型变量y_r和Y组变量简单相关系数，交叉载荷系数得到典型变量与另一组所有变量的相关性，即典型变量x_r和Y组变量、典型变量y_r和X组变量简单相关系数，绘制成三维曲线。如图6和图7所示，为典型变量x_r和X组变量分析所得载荷系数三维可视曲线图、典型变量x_r和X组变量分析所得载荷系数三维可视曲线图左视图，如图8和图9所示，为对典型变量y_r和Y组变量分析所得载荷系数三维可视曲线图以及典型变量y_r和Y组变量分析所得载荷系数三维可视曲线图左视图。

由上面典型相关分析方法分析结果可以看出，第一对典型变量和两组变量的典型载荷系数和交叉载荷系数分别大于0.5、0.6，相关性很高，这说明不同城市的多种大气污染物浓度数据具有时空相关性。

步骤二、建立四维有向GCN模型的图。

首先根据七个城市经纬度信息建立一个二维平面坐标轴如图3所示，每个城市作为图G的节点ν＝{v₁,…,v₇}，图G的边ε使用节点间向量来表达。城市v₁对城市v₂空间关联定义为边

城市v₂对城市v₁空间关联定义为边

同理可得图G的邻接矩阵A为：

对应的度矩阵D具体如下：

根据图拉普拉斯矩阵的定义，可得到拉普拉斯矩阵L具体如下：

步骤三、四维有向GCN模型的谱分解及张量运算。

在进行谱域上的图卷积计算时，需要对图的拉普拉斯矩阵进行谱分解来获得构成图傅里叶变换的基，在此前提下需要确保图的拉普拉斯为对称矩阵才能进行矩阵的谱分解，因此，利用对称矩阵和反对称矩阵的性质3即公式(7)，变换得到新形式的图邻接矩阵A′，即通过公式变换将原先的有向图转换为新无向图并保留了有向图内节点的空间方向信息。公式变换后的邻接矩阵A′具体如下：

经过变换得到新的邻接矩阵为对称矩阵，相应更新得到其度矩阵D′，具体如下：

再根据公式(8)可得到拉普拉斯矩阵L′，并进行谱分解计算相应的特征值Λ和特征向量U。

通过维数扩展和元素复制将输入的二维矩阵H变换为四维有向GCN模型的输入数据。

步骤四、四维有向GCN模型的图滤波器改进。

根据公式(13)设置图滤波器的参数k＝3，可得如下公式：

Her₃(Λ)＝2Λ·Her₂(Λ)-4Her₁(Λ)

＝2Λ·(2Λ·Her₁(Λ)-2Her₀(Λ))-4Her₁(Λ)

＝8Λ³-12Λ (29)

再将公式(29)带入公式(14)并令参数θ_k均为1，可得如下公式：

其中E为单位矩阵，最后将公式(30)带入公式(15)得到输入信号：

H*g′_θ(Λ)

＝U(θ₀M₀(Λ)+θ₁M₁(Λ)+θ₂M₂(Λ)+θ₃M₃(Λ))U^Τ·H

＝U(E+2Λ+2(2Λ²-Λ)+4(2Λ³-3Λ))U^Τ·H

＝U(8Λ³+4Λ²-12Λ+E)U^Τ·H

步骤五、四维有向GCN-LSTM模型构建及大气污染预测。

将2009年苏州、南通、常州、宜兴、湖州、嘉兴和上海市七个城市的PM2.5、PM10、CO、O₃、NO₂、SO₂六个大气污染物从第1天至第330天的浓度数据作为四维有向GCN-LSTM模型输入数据，七个城市第331天至360天的六个大气污染物浓度作为输出数据，得到苏州市第331天至360天的PM2.5浓度预测结果。对比LSTM、GCN和四维有向GCN-LSTM模型(4D-DGCN-LSTM)三种模型测试集的第1至第80个时间节点苏州市PM2.5浓度预测值和真实值曲线，如图10所示，其中第21至第40个时间节点预测值和真实值曲线图如图11所示，各模型的评估指标如表3所示：

表1.3各模型评估指标

如图10和图11所示，其中4D-DGCN-LSTM表示本发明方法，target表示目标曲线，由图中可以看出，采用本发明方法关于苏州市PM2.5浓度预测结果曲线相比于GCN模型、LSTM模型的预测结果曲线更贴近苏州市PM2.5浓度真实目标曲线，由表3的各评估指标可以看出，采用本发明方法的预测值和真实值间的误差最小，准确度最高。

除说明书所述的技术特征外，均为本专业技术人员的已知技术。本发明省略了对公知组件和公知技术的描述，以避免赘述和不必要地限制本发明。上述实施例中所描述的实施方式也并不代表与本申请相一致的所有实施方式，在本发明技术方案的基础上，本领域技术人员不需要付出创造性的劳动即可做出的各种修改或变形仍在本发明的保护范围内。

Claims

1.一种基于四维有向GCN-LSTM模型的多城市多种大气污染物预测方法，其特征在于，包括如下步骤：

步骤一，对多城市的多种大气污染物浓度数据进行相关性分析，确定多城市多种大气污染物的浓度数据具有时空相关性；

步骤二，根据所要研究的城市的地理位置建立四维有向GCN模型的图；GCN表示图卷积神经网络；

所述图为二维空间有向图，将每座城市作为图中的节点，节点位置用城市中心的经纬度表示，图中节点间都存在有向边，有向边是一节点指向另一节点的平面向量；将每种大气污染物作为节点的一个属性，属性值随时间变化；将某一时刻多城市多种大气污染物浓度数据表示为二维空间有向图，构建二维空间有向图的邻接矩阵、度矩阵和拉普拉斯矩阵；其中邻接矩阵用于描述图中节点和节点间边的信息，度矩阵用于描述图中节点与相连节点间边的权重之和；

步骤三，进行四维有向GCN模型的谱分解及张量运算；

四维有向GCN模型的谱分解是指：将二维空间有向图的反对称的邻接矩阵通过线性代数变换构建一个新的邻接矩阵，新的邻接矩阵为对称矩阵，包含旧的度矩阵和旧的邻接矩阵的信息；根据新的邻接矩阵更新度矩阵，获取新的对称的拉普拉斯矩阵，对新的拉普拉斯矩阵进行谱分解，获取特征值矩阵和特征向量矩阵；

四维有向GCN模型的张量运算是指：四维有向GCN模型的输入矩阵与所述特征向量矩阵维数不同，采用广播机制将输入矩阵进行维数扩展和元素复制，变换为与所述特征向量矩阵维数相同的三维矩阵，将变换后的输入矩阵与特征向量矩阵进行张量运算，进行二维空间有向性的信息整合；

步骤四，改进四维有向GCN模型的图滤波器；

采用Hermite滤波器作为图滤波器，对变换后的输入矩阵进行滤波；

步骤五，构建四维有向GCN-LSTM模型进行大气污染预测；LSTM表示长短期记忆网络；

所述四维有向GCN-LSTM模型是将四维有向GCN模型嵌入到LSTM网络架构里，该模型首先将当前时刻的多城市多种大气污染物浓度X_t与四维有向GCN模型上一时刻的输出数据

输出数据/>

记录下一时刻的各城市每种大气污染物的浓度预测值；

训练所述四维有向GCN-LSTM模型，利用训练好的模型进行大气污染预测。

2.根据权利要求1所述的方法，其特征在于，所述的步骤一，采用典型相关分析方法进行相关性分析，包括：

设研究对象为N个城市M种大气污染物浓度变化，首先选取N个城市M种大气污染物的当前时刻t和p个历史时刻的浓度值作为X组数据集，X组变量共N*M*(p+1)个，选取N个城市M种大气污染物的t+1时刻浓度值作为Y组，Y组变量共N*M个，如下所示：

x＝(yar₁₁(t)，，yar_NM(t)，，yar_ij(t-1)，，yar₁₁(t-p)，…，yar_NM(t-p))

Y＝(yar₁₁(t+1)，…，yar_ij(t+1)，…，Var_NM(t+1))

其中，Var表示大气污染物浓度，Var_ij(t-1)表示t-1时刻第i个城市第j种大气污染物浓度，i＝1,2,…,N，j＝1,2,…,M；

(1.1)对两组变量X和Y进行线性组合，得到第一对典型变量x₁和y₁如下：

其中，{a0₁₁,…,a1_ij,…ap_NM}表示X组的典型系数，{b₁₁,…,b_ij…,b_NM}表示Y组的典型系数；

(1.2)利用Pearson相关系数计算第一对典型变量x₁和y₁的典型相关系数P₁如下：

其中，Cov为求取协方差，Var为求取标准差；

(1.3)重复(1.1)和(1.2)，提取r对典型变量并计算典型相关系数，r＝min(X组变量个数,Y组变量个数)；

(1.4)分析r对典型变量和两组变量X、Y的典型载荷系数和交叉载荷系数；典型载荷系数表征典型变量x_r和X组变量、典型变量y_r和Y组变量的简单相关系数；交叉载荷系数表征典型变量x_r和Y组变量、典型变量y_r和X组变量的简单相关系数；载荷系数绝对值越大说明该项与典型变量之间的相关关系越强。

3.根据权利要求1所述的方法，其特征在于，所述的步骤二，包括如下步骤：

(2.1)构建二维空间有向图G＝(v,ε)，其中：v为图的节点集合，v_i表示第i个城市节点，i＝1,2,…,N；N为城市数量；ε为节点之间边的集合，节点v_i指向节点v_j的有向边e_ij为

设节点v_i、v_j的经纬度坐标分别为(x_vi,y_vi)、(x_vj,y_vj)，则/>

(2.2)构建二维空间有向图G的邻接矩阵A，邻接矩阵A的第i行第j列的元素A_ij表示节点v_i和v_j的连接关系，

(2.3)构建二维空间有向图G的度矩阵D，度矩阵D是对角矩阵，第i行第i列的元素D_ii表示与节点i相关的边的权重之和，为邻接矩阵A中第i行元素的和，度矩阵D中除了D_ii外的其余位置元素均为0；

(2.4)构建二维空间有向图G的拉普拉斯矩阵L，L为度矩阵D与邻接矩阵A之差。

4.根据权利要求1或3所述的方法，其特征在于，所述的步骤三中，四维有向GCN模型的谱分解包括：

设由步骤二获取的二维空间有向图的邻接矩阵为A，度矩阵D，拉普拉斯矩阵为L，A和L均为反对称矩阵；则由A和D构成的新的邻接矩阵A′＝DA-AD，A′为对称矩阵，则根据A′更新得到度矩阵D′，以及拉普拉斯矩阵L′，L′计算如下：

L′＝D′-A′＝D′-(DA-AD)

更新后的拉普拉斯矩阵L′为对称矩阵，对L′进行谱分解，获取特征值矩阵Λ和特征向量矩阵U如下：

其中，

表示第j个特征向量内的第k个元素，k＝1,2,…N，N为城市数量。

5.根据权利要求4所述的方法，其特征在于，所述的步骤三中，进行张量运算包括：

设四维有向GCN模型的输入矩阵H，如下：

其中，矩阵H为N*M的二维矩阵，h_ij表示根据第i个城市第j种大气污染物浓度得到的输入数据，N为城市数量，M为大气污染物种类；

特征值矩阵Λ和特征向量矩阵U均为N*N*2的三维矩阵，因此使用Numerical Python库中广播机制对H进行维数扩展和元素复制，将H变换为N*M*2的三维矩阵，如下：

变换后的矩阵H(h_NM,h_NM)^Τ表示元素h_NM在第三维度上的广播拉伸并进行一次转置；

将变换后的H与特征向量矩阵U进行张量运算。

6.根据权利要求4所述的方法，其特征在于，所述的步骤四中，采用的Hermite滤波器中Hermite多项式递归关系表示如下：

Her_k(Λ)＝2Λ·Her_k-1(Λ)-2(k-1)Her_k-2(Λ)

其中，k≥2，Her₀(Λ)＝Λ⁰，Her₁(Λ)＝2Λ，Λ为获取的特征值矩阵；

图滤波器多项式

θ_κ表示滤波器的系数，K表示图滤波器多项式截断参数；

将变换后的输入矩阵H使用Hermite滤波器，如下：

其中，ST表示H经Hermite滤波器滤波后的输出数据，U为特征向量矩阵，上角标T表示转置。

7.根据权利要求1所述的方法，其特征在于，所述的步骤五中，建立的四维有向GCN-LSTM模型中，将t时刻大气污染物浓度矩阵X_t和四维有向GCN模型的t-1时刻的输出数据

一起输入LSTM模型，在LSTM模型中进行如下计算：

C_t＝f_t·C_t-1+i_t·g_t

H_t＝O_t·tanh(C_t)

其中，f_t、i_t和O_t分别表示t时刻LSTM网络中遗忘门输出、输入门输出和输出门输出，W_f、b_f分别表示LSTM网络中遗忘门的权重和偏置，W_i、b_i分别表示LSTM网络中输入门的权重和偏置，W_o、b_o分别表示LSTM网络中输出门的权重和偏置，g_t表示一个中间写入信息，W_g和b_g分别表示中间写入信息网络的权重和偏置，C_t-1和C_t分别表示LSTM网络一个记忆单元的t-1时刻状态和t时刻状态，Sigmoid和tanh为两种激活函数；H_t为LSTM模型根据X_t和

输出的矩阵，将作为四维有向GCN模型的输入矩阵。

8.根据权利要求1或7所述的方法，其特征在于，所述的步骤五中，建立的四维有向GCN模型的计算公式如下：

ST′_t＝Sigmoid(ST_t*W_G+b_G)

i′_t＝Sigmoid(W′_i*[ST′_t,H_t]+b′_i)

g′_t＝tanh(W′_g*[ST′_t,H_t]+b′_g)

G_t＝(1+i′_t)·G_t-1+i′_t·g′_t

O′_t＝Sigmoid(W′_o·[ST′_t,H_t]+b′_o)

其中，H_t表示采用广播机制变换后的t时刻的输入矩阵；ST_t表示采用图滤波器对H_t滤波后的输出数据，Λ为特征值矩阵，U为特征向量矩阵，θ_k表示滤波器的系数，K表示图滤波器多项式截断参数；ST′_t表示t时刻四维有向GCN模型网络层的训练输出，W_G、b_G分别表示四维有向GCN模型中卷积操作的权重和偏置；i′_t表示t时刻四维有向GCN模型中输入门的输出，W′_i、b′_i分别表示四维有向GCN模型中输入门的权重和偏置；g′_t表示t时刻四维有向GCN模型中遗忘门的输出，W′_g、b′_g分别表示四维有向GCN模型中遗忘门的权重和偏置；O′_t表示t时刻四维有向GCN模型中输出门的输出，W′_o、b′_o分别表示四维有向GCN模型中输出门的权重和偏置；G_t、G_t-1分别表示t时刻、t-1时刻的空间特征的短期记忆和长期记忆整合的输出，

表示t时刻四维有向GCN模型的输出数据。