CN116384565A - 一种基于缺失数据填补的层次式大气臭氧浓度预测方法 - Google Patents

Abstract

本发明公开一种基于缺失数据填补的层次式大气臭氧浓度预测方法，首先通过季节性时间序列分解算法，将含有缺失数据的大气污染物数据以及气象数据分解为平滑趋势项、季节性周期项与短期波动项。然后使用基于时空序列建模的数据填补算法，对缺失数据进行有效的填补，获得连续数据。最后基于多层次视图臭氧预测方法，以及多次采样预测方法，构建大气臭氧浓度预测模型。本发明综合考虑了大气污染物数据以及气象数据的季节性特征，分解后的序列能够更好反映数据的情况；时空序列数据填补方法能够考虑数据的时空相关性；多层次视图臭氧浓度预测方法，能够表征污染物多尺度扩散传输情况；多次采样方法，能够获得更加稳定而有效的预测结果。

Description

一种基于缺失数据填补的层次式大气臭氧浓度预测方法

技术领域

本发明涉及一种大气臭氧浓度预测方法，其中包括基于季节性时间序列分解方法与基于时空序列建模的数据填补方法，以及基于多层次视图的臭氧浓度预测方法，属于软件应用技术领域以及数据预测技术领域。

背景技术

空气污染六项首要污染物分别为PM_2.5、PM₁₀、O₃、SO₂、NO₂以及CO。臭氧是一种具有强烈刺激性的气体，高浓度的臭氧会刺激和损害深部呼吸道，并可损害中枢神经系统，会严重影响生产生活的开展。因此对臭氧形成控制开展更好的认识，进而控制臭氧浓度是非常有必要的。同时构建空气质量臭氧浓度预测平台，提前预警与控制，也很有必要。

大气臭氧浓度预测技术分为数值方法和统计方法。数值方法基于主流的光化学烟雾反应理论，通过偏微分方程求解臭氧输送、扩散、迁移与转化的物理化学过程。目前数值方法已经发展到第三代，代表性模型包括Model3、CMAQ以及WRF-Chem等。由于数值方法基于物理化学过程驱动，机制非常复杂并且需要设置众多参数，同时对数据精度以及边界要求都比较高，因此其使用受到比较大的限制。统计方法基于统计学理论，根据与臭氧相关的空气污染物以及气象数据建立模型，其不考虑反应机理而是通过分析自变量与因变量之间的统计关系，进行预测。常见的统计模型有广义线性回归模型(GLM)、广义加性模型(GAM)等。随着深度学习方法的发展，以及算力的提高，深度学习方法在各统计方法中得到了广泛运用。深度学习方法基于深度神经网络，拥有很强的时空非线性拟合能力，在众多方法中展现出最优的效果。

目前基于深度学习的大气臭氧浓度预测方法，使用深度神经网络提取目标区域以及时间段的数据时空特征，通过解码网络进行预测得到未来目标区域以及时间的臭氧浓度。这些方法主要使用了同为时空序列预测的交通流预测模型，没有考虑空气质量情景的具体情况，基于固定的站点图提取空间特征，没有考虑大气污染物传输多层次的情况；没有结合未来的气象预测信息，仅仅使用过去的大气污染物以及气象信息。

深度学习网络输入数据包括大气污染物以及气象数据，这些数据通常来自监测站点。由于若干原因，这些监测数据通常会出现空缺，空缺长度可达数周之久。深度学习网络训练需要连续的数据，因此需要对这些空缺进行填补。目前填补方法有临近填补法，前向填补法等。这些方法只是基于临近数据对空缺数据进行填补，没有考虑污染物以及气象数据的特征，包括季节性周期起伏、长期变化趋势以及短期波动情况。

发明内容

发明目的：现有大气臭氧浓度预测方法，没有考虑大气污染场景下的特有情况。同时由于监测数据中存在大量的数据缺失情况，需要对数据进行填补，现有填补方法没有针对大气污染物数据以及气象数据的特征专门建模。为了解决这些问题，首先是数据填补问题。本发明设计了基于季节性时间序列分解以及时空序列建模的数据填补方法。季节性时间序列分解方法能够考虑污染物以及气象数据的特征，包括季节性周期起伏、长期变化趋势以及短期波动情况。时空序列建模的数据填补方法，能够综合建模单站点数据的时间相关性，以及多站点数据的空间相关性，获得更加有效的填补结果。对于预测问题，本发明设计了基于图池化-反池化的多层次视图神经网络模型。大气污染扩散传输是其浓度变化的重要环节。扩散传输需要考虑不同尺度情况，基于图池化-反池化的多层次视图神经网络模型能够考虑这种多尺度扩散传输情况，完成站点-城市-区域的多层级传输模式匹配。同时考虑到气象情况的混沌性，采取了多次采样方法，能够使得预测结果更加稳定且有效。综合上述内容，本发明提出了一种基于缺失数据填补的层次是大气臭氧浓度预测方法，简称STGUNet。

技术方案：一种基于缺失数据填补的层次式大气臭氧浓度预测方法，包括基于季节性时间序列分解、时空序列数据填补以及多层次视图臭氧浓度预测三个步骤。利用季节性时间序列分解方法，将含有缺失数据的大气污染物数据以及气象数据分解为平滑趋势项、季节性周期项与短期波动项，得到污染物以及气象数据季节性周期起伏、长期变化趋势以及短期波动情况。利用分解结果，通过基于时空序列建模的数据填补方法，对缺失数据进行有效的填补，得到连续数据。使用基于多层次视图的臭氧浓度预测方法，提取多层次视图，完成站点-城市-区域的多层级传输模式匹配，并设计多次采样方法，构建大气臭氧浓度预测模型，提高预测结果稳定性与有效性。

(1)基于季节性的时间序列分解

本发明研究数据包括大气污染数据以及气象数据。这些数据通常呈现季节性长周期变化。比如臭氧在夏季一般浓度会比较高，而在冬季浓度会比较低。同时包含长时间的趋势变化，这来源于当地长期工业生产水平、人口密度变化以及环保政策的实施。此外由于一些突发性变化，数据会出现短期的波动。基于上述分析，设计了一套季节性时间序列分解方法。

1.1局部加权回归算法

局部加权回归算法用于拟合有趋势或季节性的数据。对于时间序列S＝(T，Y)∈R^{L *2}，其中L代表时间序列的长度，Y表示序列的值，T表示序列记录的时间点。给定特征点向量X∈R^L，算法使用S＝(T，Y)拟合生成目标序列(X，Y^*)。对于每一个需要拟合的特征点x∈X，截取前后长l_p的序列

作为训练集，每个点对(t_i，y_i)表示记录的时间点以及对应的序列值，满足/>

并且截取的序列以x为中心，使得|t_k-x|与|x-t_k+lp|最小，即t_k满足：

argmin(·)表示使目标函数取最小值时的变量值。对于训练集中第i个样本点，其权值为：

σ为训练集方差。这样距离特征点x越近的样本点，其权值越大，对最终结果的贡献度越高。算法通过最小化损失函数进行训练，损失函数如下：

该函数使用矩阵形式描述为：

其中权值矩阵W＝diag(w₁，w₂，...，w_lp)。

训练过程使用梯度下降法，对损失函数求导：

令导数为0可得该函数极小值点：

θ＝(X^TWX)^-1X^TWY

参数θ即为所需。预测目标y^*＝h_θ(x)＝xθ，由此得到目标点(x，y^*)。对整个特征点向量X进行相同操作，即可得目标序列Y^*。

1.2季节性时间序列分解方法

对于时间序列Y∈R^L，其中L代表时间序列的长度，Y中包含若干缺失点。本方法将序列值向量Y表示为Y＝T+S+R，T，S，R∈R^L。其中T为平滑趋势项，S为季节性周期项，R为短期波动项。

本方法使用迭代运算。记Tⁱ，Sⁱ，Rⁱ分别为第i次迭代产生的平滑趋势项、季节性周期项与短期波动项。本方法算法如下所示：

1.2.1初始化T⁰≡0，即T⁰为全0序列。

1.2.2描述迭代算法。

对于第k次迭代得到的平滑趋势项T^k，计算去趋势化序列A^k＝Y-T^k。由于Y中包含缺失点，故A^k中也包含缺失点。

1.2.3循环子序列平滑操作。

给定周期项长度l_p，将A^k划分为长为l_p的子序列。在这些子序列上分解进行局部加权回归算法。得到回归参数θ。使用特征向量X与回归参数α进行计算，得到周期循环序列C^k。

1.2.4低通滤波操作。

对周期循环序列C^k进行低通滤波操作。一次低通滤波操作可表示为：

是一次低通滤波操作的原数据点，k表示第k次迭代，i表示在周期循环序列C^k的第i个值取/>

前后各l_da项，l_da表示第a次低通滤波操作采用的窗口长度量；低通滤波序列共2l_da+1项，得到平均值/>

对C^k每一个点进行低通滤波操作，得到周期序列/>

为一次滤波的结果。

在本阶段，总共需要进行三次滤波操作。首先取

进行第一次滤波操作，得到/>

再取/>

进行第二次滤波操作，得到/>

最后取l_d3＝3进行第三次滤波操作，得到/>

三次滤波操作使得数据达到充分平滑。第三次滤波操作后得到的/>

即为本阶段最终得到的平滑序列L^k。

1.2.5平滑循环子序列的去趋势化。

第k+1层的季节性周期项S^k+1＝C^k-L^k。这是为了防止低频能量分量进入季节性周期项。

1.2.6去季节项与趋势项平滑。

去季节项D^k+1＝Y-S^k+1。对于Y中包含的缺失值，D^k+1中也会出现相应的缺失。在整个D^k+1上使用局部加权回归算法，可以得到第k+1层的连续的平滑趋势项T^k+1。

1.2.7上文1.2.2-1.2.6的迭代算法描述了从第k次序列T^k、S^k得到T^k+1、S^k+1的过程。重复迭代算法直到前后两次季节项相差平均值小于预设的误差ε，即当mean(S^n′+1-S^n′)＜ε时，选取S^n′+1作为季节性周期项S，以及T^n′+1作为平滑趋势项T。

1.2.8获取短期波动项R＝Y-S-T。对于Y中包含的缺失值，R中也会出现相应的缺失。

(2)基于时空序列建模的数据填补

从季节性时间序列分解得到单个站点数据的平滑趋势项T，季节性周期项S以及短期波动项R。这三项包含时间序列的季节性特征，使用它们用于空缺序列填补可以获得较好的效果。其中T与S不包含空缺值，R包含空缺值。同时考虑各监测站点分布，在空间上存在关联性。因此本部分设计了基于深度神经网络的时空空缺序列填补方法。主要分为以下部分：

2.1构建训练用子图数据。

设一共有N个站点。这些站点空间上构成有向图G。整段序列长度为L，输入数据X规模为N*L，包含不同程度的空缺值。为了使得模型能够训练，需要构建训练用子数据。训练用子图数据构建具体流程如下：

2.1.1设定超参数，移动窗口长度h、样本迭代次数S以及批迭代次数I_max。

2.1.2生成随机数N_o与N_m，N_o表示可观测到的站点数，N_m表示未观测到的站点数，即数据缺失的站点。满足N_o+N_m＜N。

2.1.3从所有站点中随机选取N_o+N_m个站点，得到样本空间下标集

随机选取时间点j满足1≤j≤L-h。得到样本时间下标集J_sample＝[j，j+h]。得到子样本/>

2.1.4构建邻接矩阵

邻接矩阵每个元素W_i，j运算方法为：

其中σ是正规化距离因子。v_i，、v_j，表示对应站点的地理坐标。

2.1.5构建掩码矩阵

其值为：

2.1.6上文2.1.2-2.1.5描述的迭代步骤，一次迭代可得到子样本X＝X_sample，邻接矩阵W＝W_sample，掩码矩阵M＝M_sample。重复迭代步骤共S次，得到一次训练集合{X_1：S}，{M_1：S}，{W_1：S}，用于后续基于GNN的深度神经网络模型的填补训练。

2.1.7迭代步骤2.1.6共I_max次，共得到I_max批数据，用于训练。

2.2基于GNN的深度神经网络模型填补流程

深度神经网络训练需要多批数据，每批数据的训练方法一致。下面描述一批数据的训练过程，基于步骤2.1得到一批训练数据集合{X_sample[1：S]}，{M_sample[1：S]}，{W_sample[1：S]}。

使用基于GNN的深度神经网络模型进行训练，具体流程如下：

2.2.1对训练集合中的每一个数据，构建掩码后训练数据

其中N^*为本批训练数据站点数目。

2.2.2深度神经网络模型单层描述。

记H_l为神经网络第l层隐状态，

神经网络基本计算流程为：

其中

分别为前向转移矩阵与后向转移矩阵。其中W_sample即为上文提到的邻接矩阵，/>

表示该矩阵的转置，rowsum(·)为对矩阵的每行进行求和。k是转移次数。T_k(X)为递归定义函数，满足T_k(X)＝2XT_k-1(X)-T_k-2(X)，且T₀(X)＝I，T₁(X)＝X。/>

与/>

是第l层训练参数，用于控制每个节点中接收到信息的流动。H_l+1是第l层的输出。由于每批数据选取的节点集不同，因此邻接矩阵W_sample在每批数据能够获取不同的子图信息，并且不同批之间信息流动方法也具有差异性。

2.2.3整体模型

使用具有三层的神经网络模型。其中第0层输入

H₁使用步骤2.2.2计算流程得到。由于被掩盖的节点在第一层只能向邻居节点传播0，因此设计第二层计算流程为：

其中σ(x)＝1/(1+e^-x)。在H₂的更新算法中加入H₁是因为H₁包含含有缺失数据的节点。最后第三层使用步骤2.2.2基本计算流程得到

2.2.4损失函数

使用生成的

以及未加掩码的原始X_sample，评价两者间的总重构误差，损失函数如下：

使用该损失函数进行神经网络训练，通过梯度下降的方法优化模型参数。

2.2.5多批训练

如2.1.7所述，有I_max批数据，用于训练。对每批数据均使用2.2.1-2.2.4所述方法优化模型参数，最终得到训练好的模型。

2.2.6空缺填补

使用训练好的模型，对空缺数据进行预测即可完成填补。输入数据为N个站点，整段序列长度L，总数据X∈R^N*L。将X切分为若干长度为h的序列，X＝[X₁，X₂，...，X_k]，序列中任一子项X_i∈R^N*h。构建掩码矩阵M_i∈R^N*h，其中X_i中含有数据的点在M_i中值为1，X_i中不含数据的点在M_i中值为0。构建邻接矩阵W_i∈R^N*N。经过模型预测得到输出Y_i∈R^N*h，使用Y_i的值作为X_i中相应空缺点的值。对序列X＝[X₁，X₂，...，X_k]依次填补即可完成对输入数据的全部填补。

(3)基于多层次视图的臭氧浓度预测

将包含空缺的气象数据以及污染物数据完成填补后，便可以使用这些完整的数据构建大气臭氧浓度预测模型。

监测站点组成的网状结构可以被表示成一个有向图G＝(V，E，W)，其中V为图顶点集合，|V|＝N为站点总数。E是有向图的边集合。W∈R^N*N是带权邻接矩阵，表示了对应两个站点间的空间关系。

每个站点包含与大气臭氧污染预测的相关气象数据与污染物数据。预测模型任务是使用历史气象预计污染物数据预测未来各站点臭氧浓度。大气臭氧浓度预测模型各部分如下：

3.1图池化与反池化算法

各类气象因子以及污染物的扩散传输过程是影响臭氧生成的重要因素。监测站点的分布存在空间不均匀性。具体表现为：在城市中站点数目多，特别是重点城市，如上海有10个国控站点，南京有9个国控站点。而城市间地带，监测站点数目少，零星分布。监测站点呈现城市内聚集，密度过高；城市间又密度过低的情况。

为了表征这种分布情况，需要将扩散传输过程分层，包括城市内部站点，近距离扩散传输过程；城市间，中等距离扩散传输过程；城市间，长距离扩散传输过程。针对此需要，本发明提出了一种基于图池化的动态分层生成算法，来构建多层次网络，以表征多层扩散传输过程。图池化算法分为压缩，TopK选取、梯度传递，以及后续预测使用的反池化操作。

3.1.1压缩操作

定义X^l∈R^N*C为第l层顶点数据矩阵，W^l∈R^N*N为相应的带权邻接矩阵。压缩操作如下所示：

其中p^l∈R^C是可训练的投影向量。

表示的是各个节点经过投影后的权重值。

3.1.2 TopK选取

TopK选取操作如下：

idx＝rank(y^l，k)

W^l+1＝W^l[idx，idx]

其中rank(·)操作对权重向量y^l进行排序，并得到前k大的值的下标集合idx。k＝β·N，β∈[0，1]为池化缩小系数。记idx中这k个下标为i₁，i₂，...，i_k，其满足i_m＜i_n且1≤m＜n≤k，选取的下标顺序与原图顺序一致。通过X^l[idx，：]可以得到特征矩阵

包含X^l中前k大权重值的点，其他点则不包含在/>

中；通过W^l[idx，idx]可以得到X^l中前k大权重值的点构成的邻接矩阵W^l+1。这样原本节点数为N的图，就被池化到了节点数为k的子图。并且该生成过程是动态的，能够表征动态的扩散传输过程。

3.1.3梯度传递

由于rank(·)是离散操作，为了能让投影向量p^l可被训练，需要使用门控操作让梯度能够在反向传播中存在。其过程如下：

使用σ(·)得到门控向量

其中σ(x)＝1/(1+e^-x)。然后对矩阵/>

与/>

进行内积操作，得到下一层特征矩阵X^l+1∈R^k*C。其中X^l+1的第i行向量，是X^l第i行向量与/>

第i行值的乘积。

3.1.4反池化操作

反池化操作需要将rank(·)产生的子图恢复为原图。在图池化操作中，记录下标集合idx。借助此集合，使用逐层传播规则(layer-wise propagation rule)：

X^l＝distribute(0_N*C，X^l+1，idx)

其中idx∈Z^k包含了对应图池化算法中，从上层N节点图取出k节点的下标集合。X^{l +1}∈R^k*C是当前层特征矩阵，0_N*C是与原始图大小相等的全零特征矩阵。distribute(0_N*C，X^{l +1}，idx)根据idx下标集合将X^l+1的值分配至0_N*C，特征矩阵中。在Xi中，如果行向量下标包含在idx中，那么使用X^l+1中对应的值赋值，否则该行向量保持为0。

3.2多层次编码器-解码器网络

深度神经网络算法预测空气质量臭氧浓度的过程在于提取输入数据的时空特征，并使用时空特征进行预测，得到未来的臭氧浓度。编码器-解码器网络(Encoder-Decoder)是符合该需求的合适网络。编码器用于提取数据的时空特征，解码器用于将特征解耦至输出空间，并进行预测。Unet是一种优秀的网络结构，含有长跳跃连接、对称结构以及多尺度的网络特征，其中长跳跃连接能够有效利用不同尺度与精度的特征信息，同时能够防止网络过深出现梯度消失的问题。本发明提出一种Graph Unet，以表征气象污染物多层次扩散传输过程。其具体结构如下：

3.2.1Encoder端

对于输入X∈R^N*T*C，包含历史气象数据与污染物数据。N为站点个数，T为时间序列长度，C为数据总维度。X首先经过扩散卷积门控网络进行时空特征提取：

H_t＝u_t⊙H_t-1+(1-u_t)⊙C_t

其中σ(x)＝1/(1+e^-x)，tanh(x)＝(e^x-e^-x)/(e^x+e^-x)，均为激活函数。X_t与H_t分别为t时刻的输入与输出，H_t为隐藏矩阵，用于记录特征信息。r_t与u_t为t时刻的重置门与更新门。DGCN_W为扩散图卷积操作，W_r，W_u，W_c则是相应的可训练参数。DGCN算法每层具体形式如下：

其中U^l∈R^N*C是第l层特征矩阵，U^l+1∈R^N*C为第l+1层特征矩阵，为扩散图卷积算法在该层输出。D_o∈R^N*N＝diag(W1)，D_I∈R^N*N＝diag(W^T1)。其中D_O是邻接矩阵W的出度矩阵，D_I是邻接矩阵W的入度矩阵。1∈R^N是元素均为1的向量。θ_k，1与θ_k，2分别对应第k次扩散中的两个可训练参数矩阵。在扩散卷积中，

与/>

分别代表前向转移矩阵与后向转移矩阵。

扩散卷积使用的邻接矩阵W依赖于原图中各站点的地理位置，基于阈值高斯核构造：

其中W_i，j表示顶点v_i与v_j间的边权重，dist(v_i，v_j)表示点v_i与v_j间的欧氏距离。σ是距离的标准差，κ是阈值。

经过Encoder端，最终得到时空特征H_t∈R^N*C，并再经过自注意力操作：

K_t＝H_tW_k，V_t＝H_tW_v，Q_t＝H_tW_q

其中

分别是键向量、值向量和查询向量。

分别是可训练参数。

将经过自注意力得到的特征

在进行上文图池化操作，可得到不同层级的不同尺寸的特征图/>

3.2.2Decoder端

Decoder端输入包括三部分：1.来自Encoder端的分层时空特征

2.预测目标时间短，经其他摸型模拟得到的未来气象矩阵Wea_future；3.随机采样的种子向量。这是由于气象系统本身是个混沌系统，包含大量不确定性。加入随机向量进行多次采样，能够减小模型方差，提高预测准确性。随机向量生成器从高斯分布中随机采样向量，得到用于输入的随机向量Z_sp。

Decoder主干网络采取与Encoder端一致的扩散卷积门控网络。将随机向量与气象矩阵一并作为输入X_decoder＝[Z_sp，Wea_future]。分层时空特征

则作为隐藏矩阵输入到各层中，下一层的输出作为上一层的输入，一直到顶层。

Deocder主干网络输出

为N个站点在预测时长T′长度数据。/>

还需要经过最后的LayerNorm-FC层：

其中FC为全连接层，layernorm(·)为逐层正则化操作，对于输入x，输出y：

其中E[x]、Var[x]分别是输入x的期望与方差，∈是为数值稳定性添加到分母的值，β、γ则为可训练参数。最终得到Y＝[y₁，y₂，...，y_predlen]R^N*T′，即为单次采样下的预测。

3.2.3多次采样与训练预测

由于气象系统具有混沌性，现实情况中也有若干短期波动情况。因此需要多次采样随机向量生成器，进行多次预测，以平衡波动情况。记一次预测结果为Y_i，共采样k次，得到若干预测结果[Y₁，Y₂，...，Y_k]。对这些预测结果进行平均得到最终预测结果

训练采用专门设计的损失函数。对于时空序列预测任务，一般采用均方误差(MSE)作为损失函数以及评价指标。该函数评价预测结果与真实结果间的平均均方误差。而空气质量臭氧预测任务中，更加关注高值。高值往往意味着臭氧浓度超标，是需要预警、处理的情况。为此，设计了评价指标HVBE(High Value Balanced Error，高值平衡误差)：

HVBE(Y_true，Y_pred)＝Mean(Abs(Y_true-Y_pred)*weights)

weights＝Clip(Y_true，lowbound，upbound)

其中Mean(·)是平均值函数，Abs(·)是绝对值函数。采用Clip(·)对惩罚项weights进行裁剪。对于真实值低于下限lowbound的值，统一设置为lowbound；对于真实值高于上限upbound的值，统一设置为upbound；对于介于上下间的值，保持不变。由于weights在高值时更大，因此HVBE能够更好关注高值情况，使得模型在高值预测上更精准。但如果不加限制，模型会过于关注高值，导致整体预测偏高，因此设计了upbound。同时为了裁减掉一些过于低的值，设计了lowbound。

附图说明

图1位本发明实施例的整体流程示意图；

图2为本发明实施例的季节性空缺时间序列分解与填补方法示意图；

图3为本发明实施例的多层次时空视图空气质量臭氧预测模型示意图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图2所示，季节性空缺时间序列分解与填补方法：

(1)基于季节性的时间序列分解

给定含有空缺数据的时间序列Y∈R^L，给定特征点向量X∈R^L，将其分解为Y＝T+S+R。其中T为平滑趋势项，S为季节性周期项，R为短期波动项。使用迭代算法。对于第k次迭代得到的平滑趋势项T^k，计算去趋势化序列A^k＝Y-T^k。由于Y中包含缺失点，故A^k中也包含缺失点。给定周期项长度l_p，将A^k划分为长为l_p的子序列。在这些子序列上分解进行局部加权回归算法。得到回归参数θ。使用特征向量X与回归参数α进行计算，得到周期循环序列C^k。对周期循环序列C^k进行低通滤波操作，一次滤波操作如下：

是一次低通滤波操作的原数据点，k表示第k次迭代，i表示在周期循环序列C^k的第i个值/>

前后各l_da项，l_da表示第a次低通滤波操作采用的窗口长度量；低通滤波序列共2l_da+1项的平均值/>

对Ck每一个点进行低通滤波操作，得到周期序列/>

为一次滤波的结果。

在本阶段，总共需要进行三次滤波操作。首先取

进行第一次滤波操作，得到/>

再取/>

进行第二次滤波操作，得到/>

最后取l_d3＝3进行第三次滤波操作。以使得数据达到充分平滑。第三次滤波操作后得到的/>

即为本阶段最终得到的平滑序列L^k。。第k+1层的季节性周期项S^k+1＝C^k-L^k。这是为了防止低频能量分量进入季节性周期项。去季节项D^k+1＝Y-S^k+1。对于Y中包含的缺失值，D^k+1中也会出现相应的缺失。在整个D^k+1上使用局部加权回归算法，可以得到第k+1层的连续的平滑趋势项T^k+1。重复上述迭代过程，直到前后两次季节项相差平均值小于预设的误差ε，即当mean(S^n′+1-S^n′)＜ε时，选取S^n′+1作为季节性周期项S，以及T^n′+1作为平滑趋势项T。获取短期波动项R＝Y-S-T。对于Y中包含的缺失值，R中也会出现相应的缺失。

(2)基于时空序列建模的数据填补

输入数据一共包含N个站点，这些站点空间上构成有向图G。整段序列长度为L，总数据记为x∈R^N*L。X包含不同程度的空缺值。

设定超参数，移动窗口长度h、样本迭代次数S。生成随机数N_o与N_m，N_o表示可观测到的站点数，N_m表示未观测到的站点数，即数据缺失的站点。满足N_o+N_m＜N。从所有站点中随机选取N_o+N_m个站点，得到样本空间下标集

随机选取时间点j满足1≤j≤L-h。得到样本时间下标集J_sample＝[j，j+h]。得到子样本/>

构建邻接矩阵/>

邻接矩阵每个元素W_i，j运算方法为：

构建掩码矩阵

其值为：

重复上述迭代步骤S次，得到一批训练集合{X_sample[1：S]}，{M_sample[1：S]}，{W_sample[1：S]}，用于后续基于GNN的深度神经网络模型的填补训练。

迭代上述步骤共I_max次，共得到I_max批数据。

深度神经网络训练需要多批数据，每批数据的训练方法一致。下面描述一批数据的训练过程。对训练集合中的每一个数据，构建掩码后训练数据

其中N^*为本批训练数据站点数目。

深度神经网络模型单层描述：记H_l为神经网络第l层隐状态，

神经网络基本计算流程为：

其中

分别为前向转移矩阵与后向转移矩阵。其中W_sample即为上文提到的邻接矩阵，/>

表示该矩阵的转置，rowsum(·)为对矩阵的每行进行求和。k是转移次数。T_k(X)为递归定义函数，满足T_k(X)＝2XT_k-1(X)-T_k-2(X)，且T₀(X)＝I，T₁(X)＝X。/>

与/>

是第l层训练参数，用于控制每个节点中接收到信息的流动。H_l+1是第l层的输出。由于每批数据选取的节点集不同，因此邻接矩阵W_sample在每批数据能够获取不同的子图信息，并且不同批之间信息流动方法也具有差异性。

深度神经网络整体具有三层。其中第0层输入

H₁使用上文计算流程得到。由于被掩盖的节点在第一层只能向邻居节点传播0，因此设计第二层计算流程为：

其中σ(x)＝1/(1+e^-x)。在H₂的更新算法中加入H₁是因为H₁包含含有缺失数据的节点。最后第三层使用同H₁基本计算流程得到

使用生成的/>

以及未加掩码的原始X_sample，评价两者间的总重构误差，损失函数如下：

使用该损失函数进行神经网络训练，通过梯度下降的方法优化模型参数。

对I_max批数据均采用相同方法训练，优化模型参数，最终得到训练好的模型。

使用训练好的模型，对空缺数据进行预测即可完成填补。输入数据为N个站点，整段序列长度L，总数据X∈R^N*L。将X切分为若干长度为h的序列，X＝[X₁，X₂，...，X_k]，序列中任一子项X_i∈R^N*h。构建掩码矩阵M_i∈R^N*h，其中X_i中含有数据的点在M_i中值为1，X_i中不含数据的点在M_i中值为0。构建邻接矩阵W_i∈R^N*N。经过模型预测得到输出Y_i∈R^N*h，使用Y_i的值作为X_i中相应空缺点的值。对序列X＝[X₁，X₂，...，X_k]依次填补即可完成对输入数据的全部填补。

如图3所示，基于多层次视图的臭氧浓度预测模型：

(3)图池化与反池化过程

定义X^l∈R^N*C为第l层顶点数据矩阵，W^l∈R^N*N为相应的带权邻接矩阵。首先使用压缩操作：

其中p^l∈R^C是可训练的投影向量。

表示的是各个节点经过投影后的权重值。

接着使用TopK选取操作：

idx＝rank(y^l，k)

W^l+1＝W^l[idx，idx]

其中rank(·)操作对权重向量y^l进行排序，并得到前k大的值的下标集合idx。k＝β·N，β∈[0，1]为池化缩小系数。记idx中这k个下标为i₁，i₂，...，i_k，其满足i_m＜i_n且1≤m＜n≤k。由于选取的下标顺序与原图顺序一致通过X^l[idx，：]可以得到特征矩阵

包含X^l中前k大权重值的点，其他点则不包含在/>

中；通过W^l[idx，idx]可以得到X^l中前k大权重值的点构成的邻接矩阵W^l+1。这样原本节点数为N的图，就被池化到了节点数为k的子图。并且该生成过程是动态的，能够表征动态的扩散传输过程。

为了能让投影向量p^l可被训练，需要使用门控操作让梯度能够在反向传播中存在。其过程如下：

使用σ(·)得到门控向量

其中σ(x)＝1/(1+e^-x)。然后对矩阵/>

与/>

进行内积操作，得到下一层特征矩阵X^l+1∈R^k*C。其中X^l+1的第i行向量，是X^l第i行向量与/>

第i行值的乘积。

反池化操作需要将rank(·)产生的子图恢复为原图。在图池化操作中，记录下标集合idx。借助此集合，使用逐层传播规则(layer-wise propagation rule)：

X^l＝distribute(0_N*C，X^l+1，idx)

其中idx∈Z^k包含了对应图池化算法中，从上层N节点图取出k节点的下标集合。X^{l +1}∈R^k*C是当前层特征矩阵，0_N*C是与原始图大小相等的全零特征矩阵。distribute(0_N*C，X^{l +1}，idx)根据idx下标集合将X^l+1的值分配至0_N*C，特征矩阵中。在Xl中，如果行向量下标包含在idx中，那么使用X^l+1中对应的值赋值，否则该行向量保持为0.

(4)多层次编码器-解码器网络(Encoder-Decoder)

对于输入X∈R^N*T*C，包含历史气象数据与污染物数据。N为站点个数，T为时间序列长度，C为数据总维度。

在Encoder端，X首先经过扩散卷积门控网络进行时空特征提取：

H_t＝u_t⊙H_t-₁+(1-u_t)⊙C_t

其中DGCN(·)是基于扩散卷积的图卷积算法，算法每层具体形式如下：

其中σ(x)＝1/(1+e^-x)，tanh(x)＝(e^x-e^-x)/(e^x+e^-x)，均为激活函数。X_t与H_t分别为t时刻的输入与输出，H_t为隐藏矩阵，用于记录特征信息。r_t与u_t为t时刻的重置门与更新门。DGCN_W为扩散图卷积操作，W_r，W_u，W_c则是相应的可训练参数。DGCN算法每层具体形式如下：

其中U^l∈R^N*C是第l层特征矩阵，U^l+1∈R^N*C为第l+1层特征矩阵，为扩散图卷积算法在该层输出。D_o∈R^N*N＝diag(W1)，D_I∈R^N*N＝diag(W^T1)。其中D_o是邻接矩阵W的出度矩阵，D_I是邻接矩阵W的入度矩阵。1∈R^N是元素均为1的向量。θ_k，1与θ_k，2分别对应第k次扩散中的两个可训练参数矩阵。在扩散卷积中，

与/>

分别代表前向转移矩阵与后向转移矩阵。

扩散卷积使用的邻接矩阵W依赖于原图中各站点的地理位置，基于阈值高斯核构造：

其中W_i，j表示顶点v_i与v_j间的边权重，dist(v_i，v_j)表示点v_i与v_j间的欧氏距离。σ是距离的标准差，κ是阈值。

经过该网络，最终得到时空特征H_t∈R^N*C，并再经过自注意力操作：

K_t＝H_tW_k，V_t＝H_tW_v，Q_t＝H_tW_q

将经过自注意力得到的特征

在进行上文图池化操作，可得到不同层级的不同尺寸的特征图/>

Decoder端输入包括三部分：1.来自Encoder端的分层时空特征

2.预测目标时间短，经其他摸型模拟得到的未来气象矩阵Wea_future；3.随机采样的种子向量。这是由于气象系统本身是个混沌系统，包含大量不确定性。加入随机向量进行多次采样，能够减小模型方差，提高预测准确性。随机向量生成器从高斯分布中随机采样向量，得到用于输入的随机向量Z_sp。

Decoder主干网络采取与Encoder端一致的扩散卷积门控网络。将随机向量与气象矩阵一并作为输入X_decoder＝[Z_sp，Wea_future]。分层时空特征

则作为隐藏矩阵输入到各层中，下一层的输出作为上一层的输入，一直到顶层。

Deocder主干网络输出

为N个站点在预测时长T′长度数据。/>

还需要经过最后的LayerNorm-FC层，即为单次下的预测：

由于气象系统具有混沌性，现实情况中也有若干短期波动情况。因此需要多次采样随机向量生成器，进行多次预测，以平衡波动情况。记一次预测结果为Y_i，共采样k次，得到若干预测结果[Y₁，Y₂，...，Y_k]。对这些预测结果进行平均得到最终预测结果

训练采用专门设计的损失函数。对于时空序列预测任务，一般采用均方误差(MSE)作为损失函数以及评价指标。该函数评价预测结果与真实结果间的平均均方误差。而空气质量臭氧预测任务中，更加关注高值。高值往往意味着臭氧浓度超标，是需要预警、处理的情况。为此，设计了评价指标HVBE(High Value Balanced Error，高值平衡误差)：

HVBE(Y_true，Y_pred)＝Mean(Abs(Y_true-Y_pred)*weights)

weights＝Clip(Y_true，lowbound,upbound)

实例一：基于本发明的STGUNet的填补有效性与预测准确性评测

1、实验数据

`使用了长三角地区国控站点污染物数据。选取了纬度29°N～33°N、经度115°E～124°E间的站点，共140个。污染物种类包括PM_2.5、PM₁₀、O₃、SO₂、NO₂和CO。气象数据选取来自欧洲中期天气预报中心(ECMWF)的欧洲气象中心资料每日数据(ERA5，Daily fields)，包含云覆盖比率、空气湿度、温度、U向风场、V向风场、地势。污染物数据与气象数据单步长度1小时。选取2018年至2020年数据。其中2018年至2019年数据作为训练集，2020年数据作为验证集。

2、评价指标

将从两个方面来评价本算法的有效程度：时间序列填补方法的有效性、空气质量臭氧预测方法的准确性。同样的算法模型将会在不同对比方法消歧后的训练数据上训练，并在测试集上验证其有效性与准确性。

时间序列填补有效性评价指标为均方误差(MSE)，公式如下所示：

MSE(Y_true，Y_pred)＝(Y_true-Y_pred)²

空气质量臭氧预测方法准确性评价指标包括HVBE(High Value Balanced Error，高值平衡误差)与F1-score，公式如下所示：

HVBE(Y_true，Y_pred)＝Mean(Abs(Y_true-Y_pred)*weights)

weights＝Clip(Y_true，lowbound,upbound)

F1-score为评价模型对臭氧超标情况预测的准确性。选取《中华人民共和国环境空气质量标准》中规定的臭氧1小时平均一级浓度限值160g/μm³与1小时平均二级浓度限值200g/μm³作为分界线。其中TP为预测超标，实际也超标的样本数；FP为预测超标，实际未超标的样本数；FN为预测未超标，实际超标的样本数。precision为预测超标的数据中，实际也超标样本的比例；recall为所有实际超标的样本中，有多少能被正确预测的比例。由于一般情况下precision与recall难以兼顾，因此使用F1-score进行综合考量。HVBE(·)越小表示预测结果越好，F1-score越大表示预测结果越好。

3、比较方法

为了更好地验证本发明提出的STGUNet效果，对于时间序列填补有效性，设置以下两个对比方法进行比较：

(1)Linear-Inter，使用空缺数据前后有值点线性填补。

(2)STGUNet-NoDec，不进行季节性时间序列分解方法，直接对源数据使用本发明中基于深度神经网络的时空空缺序列填补方法。

对于空气质量臭氧预测方法的准确性，设置以下三个对比方法进行比较：

(1)STGUNet-NoPool，不使用图池化与反池化算法，Graph Unet的每层特征图均为原图。

(2)STGUNet-NoSample，不使用随机向量生成器进行多次采样方法，将单次预测结果直接作为最终结果。

(3)STGUNet-MSE，采用MSE作为损失函数训练。

(4)CLRNN，使用基于球面卷积算法的时空序列预测模型。

4、实验设置

对于时间序列填补实验，季节性时间序列分解部分设置移动窗口长度为4380，样本迭代次数100，批迭代次数128。时空序列建模的数据填补部分，神经网络模型批大小128，中间层维度均为128，学习率0.001，训练轮次100轮。

对于空气质量臭氧预测实验，数据全部经过归一化处理，设置历史数据长度48、预测数据长度24。中间层维度均为128,。图池化层数2，每层池化缩小系数均为0.8。神经网络模型批大小128，学习率0.001，训练轮次100轮。

5、实验结果

首先是时间序列填补实验，表1展示了具体实验结果。可以看到本发明提出的基于季节性周期时间序列分解与深度神经网络的时空空缺序列填补方法，在填补有效性上显著优于Linear-Inter和STGUNet-NoDec方法。

表1时间序列填补算法有效性对比

方法	MSE
		STGUNet	54.877
Linear-Inter	113.80
		STGUNet-NoDec	79.771

对于空气质量臭氧预测实验，表2展示了具体实验结果。在预测准确度上显著优于STGUNet-NoPool，STGUNet-NoSample，STGUnet-MSE，CLRNN。

表2空气质量臭氧预测准确性对比

方法	HVBE	L1 F1-Score	L2 F1-Score
				STGUNet	3.457	0.6263	0.4801
STGUNet-NoPool	3.993	0.5149	0.3335
				STGUNet-NoSample	3.623	0.5317	0.3849
STGUNet-MSE	4.307	0.4576	0.2628
				CLRNN	4.723	0.3638	0.1388

通过对实验结果的分析，使用本发明提供的方法，通过季节性周期序列分解，以及使用深度神经网络能够有效完成空缺时间序列的填补，获得高质量的连续数据。同时使用图池化与反池化的多层次视图深度神经网络模型，能够有效表征多层次扩散传输过程，获得了最优的预测结果。综上所述，证明了本发明的有效性。

Claims (7)

1.一种基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，将大气污染物数据以及气象数据作为输入数据，利用基于季节性的时间序列分解方法，分解含有空缺的输入数据，得到污染物以及气象数据季节性周期起伏、长期变化趋势以及短期波动情况的分解结果；对分解结果使用基于时空序列建模的数据填补模型进行填补；使用基于多层次视图的臭氧浓度预测模型进行预测，并采用多次采样预测方法，生成未来区域稳定的臭氧预测结果，方法包括以下步骤：

(1)基于季节性的时间序列分解；

(2)基于时空序列建模的数据填补；

(3)基于多层次视图的臭氧浓度预测。

2.根据权利要求1所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，基于季节性的时间序列分解，包括如下内容：

1.1局部加权回归算法

首先定义局部加权回归算法；对于时间序列S＝(T，Y)∈R^L*2，其中L代表时间序列的长度，Y表示序列的值，T表示序列记录的时间点；给定特征点向量X∈R^L，使用S＝(T，Y)拟合生成目标序列(X，Y^*)；对于每一个需要拟合的特征点x∈X，截取前后长l_p的序列

作为训练集，每个点对(t_i，y_i)表示记录的时间点以及对应的序列值，满足t_k＜t_k+1＜…＜t_k+lp、并且截取的序列以x为中心，使得|t_k-x|与

最小，即t_k满足：

argmin(·)表示使目标函数取最小值时的变量值；对于训练集中第i个样本点，其权值为：

σ为训练集方差；距离特征点x越近的样本点，其权值越大，对最终结果的贡献度越高；算法通过最小化损失函数进行训练，损失函数如下：

该函数使用矩阵形式描述为：

其中权值矩阵W＝diag(w₁，w₂，...，w_lp)；

求解方法使用梯度下降法，对损失函数求导：

令导数为0可得该函数极小值点：

θ＝(X^TWX)^-1X^TWY

参数θ即为所需，预测目标y^*＝h_θ(x)＝xθ；由此得到目标点(x，y^*)，对整个特征点向量X进行相同操作，即可得由预测目标y^*组成的目标序列Y^*；

1.2季节性时间序列分解方法

对于包含缺失点的时间序列Y，将其分解为Y＝T+S+R；其中T为平滑趋势项，S为季节性周期项，R为短期波动项；使用迭代运算，记Tⁱ，Sⁱ，Rⁱ分别为第i次迭代产生的平滑趋势项、季节性周期项与短期波动项。

3.根据权利要求1所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，季节性时间序列分解方法包括：

1.2.1，初始化T⁰≡0，即T⁰为全0序列；

1.2.2，描述迭代算法

对于第k次迭代得到的平滑趋势项T^k，计算去趋势化序列A^k＝Y-T^k，由于Y中包含缺失点，故A^k中也包含缺失点；

1.2.3，循环子序列平滑操作

给定周期项长度l_p，将A^k划分为长为l_p的子序列，在这些子序列上分解进行局部加权回归算法，得到回归参数θ，使用特征向量X与回归参数θ进行计算，得到周期循环序列C^k；

1.2.4，低通滤波操作

对周期循环序列C^k进行低通滤波操作，一次低通滤波操作可表示为：

是一次低通滤波操作的原数据点，k表示第k次迭代，i表示在周期循环序列C^k的第i个值取/>

前后各l_da项，l_da表示第a次低通滤波操作采用的窗口长度量；低通滤波序列共2l_da+1项，得到平均值/>

对C^k每一个点/>

进行低通滤波操作，得到周期序列/>

为一次滤波的结果；

在本阶段，总共需要进行三次滤波操作；首先取

进行第一次滤波操作，得到

再取/>

进行第二次滤波操作，得到/>

最后取l_d3＝3进行第三次滤波操作，得到

三次滤波操作使得数据达到充分平滑；第三次滤波操作后得到的/>

即为本阶段最终得到的平滑序列L^k；

1.2.5，平滑循环子序列的去趋势化

第k+1层的季节性周期项S^k+1＝C^k-L^k，这是为了防止低频能量分量进入季节性周期项；

1.2.6，去季节项与趋势项平滑

去季节项D^k+1＝Y-S^k+1，对于Y中包含的缺失值，D^k+1中也会出现相应的缺失，在整个D^k+1上使用局部加权回归算法，能得到第k+1层的连续的平滑趋势项T^k+1；

1.2.7，上文1.2.2-1.2.6的迭代算法描述了从第k次序列T^k、S^k得到T^k+1、S^k+1的过程；重复迭代算法直到前后两次季节项相差平均值小于预设的误差ε，即当mean(S^n′+1-S^n′)＜ε时，选取S^n′+1作为季节性周期项S以及T^n′+1作为与平滑趋势项T；

1.2.8，获取短期波动项R＝Y-S-T；对于Y中包含的缺失值，R中也会出现相应的缺失。

4.根据权利要求1所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，所述基于时空序列建模的数据填补，包括如下内容：

从季节性时间序列分解得到单个站点数据的平滑趋势项T，季节性周期项S以及短期波动项R；其中T与S不包含空缺值，R包含空缺值；同时考虑各监测站点分布，在空间上存在关联性；因此设计了基于深度神经网络的时空空缺序列填补方法，分为以下部分：

2.1，构建训练用子图数据

设一共有N个站点，这些站点空间上构成有向图G，站点间存在关联性，整段序列长度为L，输入数据X规模为N*L，包含不同程度的空缺值；训练用子图数据构建具体流程如下：

2.1.1，设定超参数，移动窗口长度h、样本迭代次数S以及批迭代次数I_max；

2.1.2，生成随机数N_o与N_m，N_o表示可观测到的站点数，N_m表示未观测到的站点数，即数据缺失的站点，满足N_o+N_m＜N；

2.1.3，从所有站点中随机选取N_o+N_m个站点，得到样本空间下标集

随机选取时间点j满足1≤j≤L-h；得到样本时间下标集J_sample＝[j，j+h]；得到子样本/>

2.1.4，构建邻接矩阵

邻接矩阵每个元素W_i，j运算方法为：

其中σ是正规化距离因子；v_i，、v_j，表示对应站点的地理坐标；

2.1.5，构建掩码矩阵

其值为：

2.1.6，2.1.2至2.1.5描述的迭代步骤，重复迭代步骤共S次，得到一次训练集合{X_sample[1：S]}，{M_sample[1：S]}，{W_sample[1：S]}，用于后续基于GNN的深度神经网络模型的填补训练；

2.1.7，迭代步骤2.1.6共I_max次，共得到I_max批数据，用于训练；

2.2，基于GNN的深度神经网络模型填补流程

深度神经网络训练需要多批数据，每批数据的训练方法一致；下面描述一批数据的训练过程，基于步骤2.1中得到训练集合{X_sample[1：S]}，{M_sample[1：S]}，{W_sample[1：S]}；使用基于GNN的深度神经网络模型进行训练，具体流程如下：

2.2.1，对训练集合中的每一个数据，构建掩码后训练数据

其中N^*为本批训练数据站点数目；

2.2.2，深度神经网络模型单层描述

记H_l为神经网络第l层隐状态，

神经网络基本计算流程为：

其中

分别为前向转移矩阵与后向转移矩阵；其中W_sample即为上文提到的邻接矩阵，/>

表示该矩阵的转置，rowsum(·)为对矩阵的每行进行求和；k是转移次数；T_k(X)为递归定义函数，满足T_k(X)＝2XT_k-1(X)-T_k-2(X)，且T₀(X)＝I，T₁(X)＝X；/>

与/>

是第l层训练参数，用于控制每个节点中接收到信息的流动；H_l+1是第l层的输出；

2.2.3，整体模型

使用具有三层的神经网络模型；其中第0层输入

H₁使用步骤2.2.2计算流程得到；由于被掩盖的节点在第一层只能向邻居节点传播0，因此设计第二层计算流程为：

其中σ(x)＝1/(1+e^-x)。在H₂的更新算法中加入H₁是因为H₁包含含有缺失数据的节点。最后第三层使用步骤2.2.2基本计算流程得到

2.2.4，损失函数

使用生成的

以及未加掩码的原始X_sample，评价两者间的总重构误差，损失函数如下：

使用该损失函数进行神经网络训练，通过梯度下降的方法优化模型参数；

2.2.5，多批训练

如步骤217所述，有I_max批数据，用于训练；对每批数据均使用步骤221-22.4所述方法优化模型参数，最终得到训练好的模型；

2.2.6，空缺填补

使用训练好的模型，对空缺数据进行预测即可完成填补；输入数据为N个站点，整段序列长度L，总数据X∈R^N*L；将X切分为若干长度为h的序列，X＝[X₁，X₂，...，X_k]，序列中任一子项X_i∈R^N*h；构建掩码矩阵M_i∈R^N*h，其中X_i中含有数据的点在M_i中值为1，X_i中不含数据的点在M_i中值为0；构建邻接矩阵W_i∈R^N*N；经过模型预测得到输出Y_i∈R^N*h，使用Y_i的值作为X_i中相应空缺点的值；对序列X＝[X₁，X₂，...，X_k]依次填补即可完成对输入数据的全部填补。

5.根据权利要求1所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，基于多层次视图的臭氧浓度预测，使用图池化与反池化算法获得多层次视图；图池化与反池化算法分为压缩，TopK选取、梯度传递，以及后续预测使用的反池化操作；

压缩操作中，定义X^l∈R^N*C为第l层顶点数据矩阵，W^l∈R^N*N为相应的带权邻接矩阵，压缩操作如下所示：

其中p^l∈R^C是可训练的投影向量，

表示的是各个节点经过投影后的权重值；

TopK选取操作如下：

idx＝rank(y^l，k)

W^l+1＝W^l[idx，idx]

其中rank(·)操作对权重向量y^l进行排序，并得到前k大的值的下标集合idx，k＝β·N，β∈[0，1]为池化缩小系数，记idx中这k个下标为i₁，i₂，...，i_k，其满足i_m＜i_n且1≤m＜n≤k，由于选取的下标顺序与原图顺序一致，通过X^l[idx，：]得到特征矩阵

包含X^l中前k大权重值的点，其他点则不包含在/>

中；通过W^l[idx，idx]，得到X^l中前k大权重值的点构成的邻接矩阵W^l+1，这样原本节点数为N的图，就被池化到了节点数为k的子图，并且该生成过程是动态的；

梯度传递中，由于rank(·)是离散操作，为了能让投影向量p^l可被训练，需要使用门控操作让梯度能够在反向传播中存在，其过程如下：

使用σ(·)得到门控向量

其中σ(x)＝1/(1+e^-x)，然后对矩阵/>

与/>

进行内积操作，得到下一层特征矩阵X^l+1∈R^k*C，其中X^l+1的第i行向量，是X^l第i行向量与/>

第i行值的乘积；

反池化操作需要将rank(·)产生的子图恢复为原图，在图池化操作中，记录下标集合idx，借助此集合，使用逐层传播规则：

X^l＝distribute(0_N*C，X^l+1，idx)

其中idx∈Z^k包含了对应图池化算法中，从上层N节点图取出k节点的下标集合；X^l+1∈R^k*C是当前层特征矩阵，0_N*C是与原始图大小相等的全零特征矩阵；distribute(0_N*C，X^l+1，idx)根据idx下标集合将X^l+1的值分配至0_N*C，特征矩阵中；在X^l中，如果行向量下标包含在idx中，那么使用X^l+1中对应的值赋值，否则该行向量保持为0。

6.根据权利要求1所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，基于多层次视图的臭氧浓度预测，预测模型主体是多层次编码器-解码器网络，对于输入X∈R^N*T*C，包含历史气象数据与污染物数据，N为站点个数，T为时间序列长度，C为数据总维度，在Encoder端，X首先经过扩散卷积门控网络进行时空特征提取：

H_t＝u_t⊙H_t-1+(1-u_t)⊙C_t

其中σ(x)＝1/(1+e^-x)，tanh(x)＝(e^x-e^-x)/(ex+e^-x)，均为激活函数，X_t与H_t分别为t时刻的输入与输出，H_t为隐藏矩阵，用于记录特征信息，r_t与u_t为t时刻的重置门与更新门，DGCN(·)表示基于扩散卷积的图卷积算法，W_r，W_u，W_c则是相应的可训练参数，DGCN算法每层具体形式如下：

其中U^l∈R^N*C是第l层特征矩阵，U^l+1∈R^N*C为第l+1层特征矩阵，为扩散图卷积算法在该层输出，D_o∈R^N*N＝diag(W1)，D_I∈R^N*N＝diag(W^T1)，其中D_O是邻接矩阵W的出度矩阵，D_I是邻接矩阵W的入度矩阵，1∈R^N是元素均为1的向量，θ_k，1与θ_k，2分别对应第k次扩散中的两个可训练参数矩阵，在扩散卷积中，

与/>

分别代表前向转移矩阵与后向转移矩阵；

扩散卷积使用的邻接矩阵W依赖于原图中各站点的地理位置，基于阈值高斯核构造：

其中W_i，j表示顶点v_i与v_j间的边权重，dist(v_i，v_j)表示点v_i与v_j间的欧氏距离；σ是距离的标准差，κ是阈值；

经过Encoder端，最终得到时空特征H_t∈R^N*C，并再经过自注意力操作：

K_t＝H_tW_k，V_t＝H_tW_v，Q_t＝H_tW_q

将经过自注意力得到的特征

在进行上文图池化操作，可得到不同层级的不同尺寸的特征图/>

Decoder端输入包括三部分：(1)来自Encoder端的分层时空特征

预测目标时间段，经其他摸型模拟得到的未来气象矩阵Wea_future；(3)随机采样的种子向量；加入随机向量进行多次采样，随机向量生成器从高斯分布中随机采样向量，得到用于输入的随机向量Z_sp。

7.根据权利要求6所述的基于缺失数据填补的层次式大气臭氧浓度预测方法，其特征在于，Decoder主干网络采取与Encoder端一致的扩散卷积门控网络，将随机向量与气象矩阵一并作为输入X_decoder＝[Z_sp，Wea_future]，分层时空特征

则作为隐藏矩阵输入到各层中，下一层的输出作为上一层的输入，一直到顶层；

Deocder主干网络输出

为N个站点在预测时长T′长度数据，/>

还需要经过最后的LayerNorm-FC层，即为单次下的预测：

多次采样随机向量生成器，进行多次预测，以平衡波动情况，记一次预测结果为Y_i，共采样k次，得到若干预测结果[Y₁，Y₂，...，Y_k]，对这些预测结果进行平均得到最终预测结果

对于时空序列预测任务，采用均方误差作为损失函数以及评价指标，该函数评价预测结果与真实结果间的平均均方误差；而空气质量臭氧预测任务中，更加关注高值，高值意味着臭氧浓度超标，是需要预警、处理的情况；为此，设计了评价指标HVBE(High ValueBalanced Error，高值平衡误差)：

HVBE(Y_true，Y_pred)＝Mean(Abs(Y_true-Y_pred)*weights)

weights＝Clip(Y_true，lowbound，upbound)

Images