基于时空融合的缺失值填补算法

2020-02-08苏思凡朱文军赵枫朝

计算机工程与设计 2020年1期

苏思凡，竹翠，朱文军，赵枫朝

(北京工业大学信息学部，北京 100022)

0 引言

空气污染物预报的方法可分为机理模型和机器学习模型。前者基于空气动力学知识对污染物的产生和扩散建模。代表算法有CMAQ和WRF-Chem。其理论性强，能对污染成因进行诊断，但需要使用者有丰富的专业知识，时间与计算资源开销大，还需详细的污染源资料。为寻找更加简便的方式，自回归移动平均[1]、线性回归[2]、支持向量机[3]、神经网络[4,5]等机器学习模型被广泛尝试。长短时记忆神经网络[6]、深度置信网络[7]等深度学习模型也用于该领域。此类方法不用详述特征与目标间的关系，预测结果可精确到点位，但无法对污染成因给出原理上的解释。

从数值角度看，PM2.5浓度是一组时间序列。此外，其它地区的污染物也在不断扩散，任何地区的PM2.5浓度都不能“独善其身”，势必会受到相邻地区的影响[8]。本文分别从时、空两个维度进行预测，再将预测结果融合，对不同气象条件下的PM2.5有针对性地回归。

由于传感器故障，数据集往往存在大量缺失数据。通过填补缺失值弥补其对预测带来损失也是本文研究重点之一。尽管缺失值填补已有一些常用的方法[9,10]，但这些方法很少能充分利用数据的时空相关性。因此，本文结合时空背景提出了一种缺失值填补算法，辅助PM2.5浓度预测模型进一步精提升精度。

1 模型原理简述

1.1 人工神经网络(ANN)原理简述

ANN一般由3类结构组成：输入层、隐藏层和输出层。每一层包含一个或多个神经元。神经元将前一层的输出作为当前层的输入，将这些输入乘以不同的权值，可以实现对输入信号的放大与抑制。然后利用激活函数将计算后的值压缩到某一区间后传递给下一层的神经元。数据经过层与层之间连接最终传导到输出层，由输出层神经元将结果送出。

对ANN的训练就是根据输出结果与真实值计算误差，然后根据误差从后向前不断调整每层神经元之间权值的过程。这个过程被称为“误差的反向传播”。该过程利用梯度下降的方法对参数进行调整，使模型的输出值与真实值之间的误差尽可能减小，从而不断优化模型，获得理想的输出结果。

1.2 长短时记忆神经网络(LSTM)原理简述

LSTM[11]是循环神经网络的变体。与ANN不同的是，LSTM中的神经元除了可以把输出传递到下层神经元，还引入了时序的概念，可以把当前的输出作为自身下一时刻的输入，能够让信息在网络中循环。LSTM与ANN相比，最大的优势在于它能将网络之前获得的信息应用在当前神经元的计算中。LSTM的神经元细胞中包含“输入门”、“忘记门”和“输出门”，这种结构能有选择性地记住和忘记曾经获得的信息，并解决长期依赖问题。正是由于这些特点，LSTM可以从以前的数据中提取知识，从而为当前的预测提供帮助。

1.3 模型树(MT)原理简述

MT[12]是将决策树和线性回归算法进行融合的一种模型。它首先利用树构建算法寻找当前最佳的切分属性与最佳的切分点，不断将数据集进行二元切分。与普通的决策树不同的是，MT中的每个叶子结点并不是所在分支对应的标签，而是基于划分后的子样本集所建立的线性回归模型。与线性回归算法不同的是，MT的结果是分段线性的，而线性回归算法是全局线性的。当数据集较大且问题的影响机制复杂时，使用普通线性回归算法显然不够灵巧，不能很好地对全体数据进行拟合。而MT将全体数据切分成多组易建模的数据后再回归，能够更好地处理非线性问题，同时提升了算法的可解释性和预测准确性。

2 基于时空融合的预测模型

一个地区的气象数据和空气污染物数据通常由分布在该地区的多个气象和空气质量监测站获取并发布的。本模型也以站点为粒度对PM2.5浓度进行预测。针对每一个站点，分别建立一套预测模型。每套模型的结构相同，但分别依赖各自站点的数据进行训练与预测。

该模型由3个子模型构成：基于时间维度的LSTM预测器、基于空间维度的ANN预测器和基于时空融合的MT预测器，如图1所示。其中LSTM负责从时间序列中获取PM2.5的变化规律，从时间维度进行预测；ANN负责学习站点间PM2.5的影响机制，从空间维度进行预测；而MT负责将时、空间维度的结果整合起来，获取最终的预测结果。

图1 基于时空融合的PM2.5浓度预测总体架构

2.1 基于时间维度的LSTM预测器

LSTM可以从历史数据中获取变量在时间维度上呈现出的趋势性和周期性规律。本文利用LSTM通过输入某一站点前12小时的逐小时PM2.5浓度以及气象数据(包含6个变量：天气、湿度、气压、风向、风速、温度)，对该站点下一小时的PM2.5浓度进行多因子的时间序列预测。

在LSTM模型的结构设置上，本文使用了一层LSTM叠加dropout，再叠加一层全连接层的结构来进行训练与测试，如图2所示。模型使用mae作为损失函数，adam为优化算法。在训练前，首先要将每12个相邻小时的气象和PM2.5浓度数据处理为一组时间序列，作为训练数据输入模型。输入数据首先通过LSTM层，提取出气象和PM2.5随时间变化的时序信息与规律。由于LSTM的学习能力十分强大，为防止过拟合而出现模型泛化性不好的问题，本文在LSTM层结束后加入了dropout，利用以一定比例随机断开神经元间连接的策略对网络进行正则化。经过多次实验，本文选择以0.01作为dropout的参数，即在模型训练过程中，会有1%的神经元的连接会被随机断开。经过dropout处理过的数据将通过全连接层进行特征提取与特征映射，根据输入的气象和PM2.5数据计算出下一小时的PM2.5浓度值。

图2 基于时间维度的LSTM预测器结构

2.2 基于空间维度的ANN预测器

某地区的当前的PM2.5浓度值除了会受到本地气象扩散条件和前期PM2.5水平的影响，周围地区的PM2.5也会通过大气运动传输到本地。因此，对某一地区PM2.5浓度的预测还要考虑空间维度上存在的相关性。由于PM2.5所呈现的空间相关性除了和距离有关，还会受到空气对流情况的影响，度量起来较为复杂。为了降低特征工程的难度并获得一种可靠的空间预测模型，本文利用ANN来模拟周围地区对本地PM2.5浓度的影响。

基于空间维度的ANN预测器使用了与“基于时间维度的LSTM预测器”类似的序贯模型。该网络由3个全连接层和dropout组成。dropout的参数设为0.01，全连接层的激活函数选用relu。该模型的输入数据由本地与其邻居站点前一小时的PM2.5浓度、本地站点前一小时的风速与风向数据组成。

对于图3而言，星型代表当前研究的站点，我们以它作为中心，向四周发散出8条射线，每条线间呈45度角。这8条射线将研究站点的周围划分成8个区域。对于每一块区域，我们取离研究站点最近的站点为邻居站点。按照上述规则，每一个站点最多有8个邻居站点，分别来自于不同的方向。图中圆圈代表邻居站点，方块代表非邻居站点。当研究站点前一小时的风向为东北风时，我们通常认为邻居站点A和邻居站点B对当前时刻研究站点的PM2.5浓度水平有着较为直接的影响。我们把当前研究站点及其所有邻居站点前一小时的PM2.5浓度数据以及当前研究站点前一小时的风速风向数据输入ANN中，试图让网络学习到类似以上风对站点之间PM2.5浓度值的影响机制，从而从空间角度对PM2.5水平进行回归。如果直接将所有站点的PM2.5浓度数据全部输入到神经网络中，输入的数据量过大，为了得到相对稳定和准确的效果，势必要增大神经网络的规模，这将会提升调整参数的难度，也会导致训练的时间过长，甚至模型震荡难以收敛。

图3 邻居站点

2.3 基于时空融合的MT预测器

有文献指出，如果各个神经网络的误差互不相关，将这些神经网络作为子预测器进行集成后，总体的预测效果与泛化能力将获得明显提升。为了尽量提升子预测器的独立性，本文分别利用不同的特征构造不同的样本，并利用不同类型的神经网络从不同的角度进行预测。现在我们已分别获得了从时间和空间角度的PM2.5浓度预测结果。当对这两个结果进行融合时，需要考虑如何进行权重分配的问题。从问题本身来看，当某一站点所处区域的空气流动性较差时，邻居站点的PM2.5很难通过大气运动传输到本地。此时，利用ANN所得的结果可能不如利用LSTM所得的结果影响力大。相反，当该站点所处区域的风速较大时，LSTM所得的局部预测结果就应该弱化，而ANN所得的预测结果应该适当得到加强。因此对于这种较为复杂的问题，时、空预测结果的权值不应一概而论，应由当时的气象扩散条件决定。时空融合预测器所使用的MT算法可以很好地解决这个问题，其示意图如图4所示。MT中每个线性回归方程所确定的参数都适应于当前的气象条件，比对整体样本统一进行一次回归有更好的针对性与适应性。

图4 基于时空融合的MT预测器结构

3 基于时空融合的缺失值填补算法

缺失数据会使数据中蕴含的规律更难把握，阻碍模型发现与学习特征之间的影响机制。在实际建模中，缺失值的危害主要体现在数据量不足而导致的模型的泛化性差，输出结果不可靠等方面。此外，由于传感器硬件故障或通信失败等情况，数据集经常会存在缺失数据。且本文所使用的数据集是精确到小时粒度的，如果以上故障排除不能及时排除将导致连续缺失值的出现，严重降低了数据质量。为了解决传感器网络缺失的问题，并进一步提高模型的预测性能，本文提出了一种基于时空融合的缺失值填补算法。该算法的整体思路是利用站点自身的时间关系和站点间的空间关系找到与缺失值最相近的4个参考值，然后将相关系数作为这4个参考值对应的权重。对于连续性变量，将4个参考值进行加权平均后的结果作为填充值。对于非连续性变量，将4个参考值中权重最大的值作为填充值。本方法是对热卡填充法进行优化后的一种算法，与普通热卡填充法相比可以有更多的参考值，并从时间和空间两个维度寻找相似的样本，提高了填补算法的鲁棒性与准确性。

3.1 数据集结构说明

对于n个数据采集站点在相同时间粒度下获得的观测值，可以以图5的方式进行组织排列。其中，每一列代表一个数据采集站点获得的观测值的时间序列，每一行代表某一时刻所有站点的观测量的取值情况。对于记录中的缺失值，用符号“na”表示。

图5 数据集结构

对于图5，我们以3号站点t时刻的缺失值作为研究对象，对其进行填补。

3.2 时间维度

对于时间序列，变量历史的情况会对当前和未来时刻的情况产生影响，当前时刻的变量是在前一时刻的基础上发展而来的。特别是对于本文使用的小时粒度的气温、气压等气象数据和PM2.5浓度数据，这些变量在一小时内很少会出现大幅度跳变的情况。所以对于时间间隔较小的时间序列，可以简单地认为与当前时刻与前后两个时刻的指标水平比较相似。因此，对于所研究的t时刻站点的取值序列，我们认为与t-1时刻和t+1时刻的取值序列最为相近，其相似程度选用皮尔逊相关系数进行量化，公式如式(1)所示。R>0表示变量X和Y间存在正相关关系。R越大，表明两组变量的正相关程度越强。当R=1时，X和Y呈完全正相关线性关系。反之亦然。R=0表示两者间无线性关系。本文利用皮尔逊相关系数分别计算t-1时刻序列与t时刻序列的相似度cor1、t时刻序列与t+1时刻序列的相似度cor2

(1)

val1和val2可以看作是与当前缺失值最相似的两个参考值。本算法利用序列的相似度cor1和cor2分别代表val1与当前缺失值、val2与当前缺失值的相似度，并以此分别作为val1、val2的权值，参与3.4部分时空融合填补的算法。

3.3 空间维度

在空间维度，由于数据采集站点散布在研究区域的不同位置，站点间的距离与位置关系比较复杂，每个站点所处的环境也各不相同，因此很难通过量化这些外部因素的方式来判断空间中与哪些站点与当前站点最相似。为了找到衡量空间相似度的标准，与时间维度不同，本文直接从各站点采集到的观测值序列的角度出发，通过计算序列间的相关系数矩阵，从中找到与当前站点在观测值表现上最相似的其它站点。

对于图5，通过计算n列两两之间的相关系数，我们可得到一个大小为n*n的相关系数矩阵。易知该矩阵是对角线为1的对称矩阵。由于皮尔逊相关系数的绝对值越大表明两组变量相关程度越强。因此，遍历相关系数矩阵中的每一列，排除对角线元素，取其中绝对值最大的前两个元素，这两个元素对应的行标签就是与当前站点相似度最高的两个站点。对于图5中当前研究的缺失值，假设通过上述方法找到与其所在的3号站点在相似度最高的两个站点分别为1和n-1号站点，那么与当前缺失值同一时刻下的val3和val4被认为是空间维度下与当前缺失值相似度最高的两个参考值。通过相关系数矩阵找到的当前列绝对值最大的前两个元素cor3和cor4分别对应为val3和val4的权值，参与3.4部分时空融合填补的算法。

3.4 时空融合

经过3.2和3.3两个部分，对于当前缺失值，我们已从时间和空间维度找到了与其相似度较高的4个参考值，并确定了它们各自的权重。本算法从时间和空间维度共选取4个参考值的原因在于，本课题所用数据中的缺失值较多，如果在时、空维度分别仅选取一个参考值，经常会出现两个参考值同时都是缺失值的情况，导致参考值不能发挥作用，缺失值依然不能被填补的情况发生。为了在保证填补精度的前提下尽量提高填补效率，本算法分别从时空维度各选取两个参考值，以降低参考值同时为缺失值的概率。

在对缺失值进行填补时，需要区分两种情况：一是待填补值是连续型变量，二是待填补值是非连续性变量。对于前者，可使用加权平均的思想。对于后者，由于待填补值非数值，且只能是该特征属性集合中的某一个，因此加权平均的方式并不适用，这里选择权值最大的参考值作为填补值。需要说明的是，当某一个或多个参考值为缺失值时，需要将它和它所对应的权值分别从参考值集合和权重集合中剔除，使其不参与时空融合缺失值填补的计算。此时，参考值的总个数从4个退化为不足4个。当且仅当4个参考值全缺失时，当前缺失值无法填补。

填补值fill_val的计算方式如式(2)所示。并利用同样的方法对整个数据集缺失值进行填补

(2)

4 实验及分析

本文共设计了3组实验：时空融合与单维度预测的效果对比、不同时空融合预测器的效果对比、两种缺失值填补方法的效果对比。第1组实验是为了验证本文提出的基于时空融合的预测框架比从时间或空间某一维度进行独立预测的优越性。第2组实验是为了验证在时空融合预测框架的前提下，MT比其它主流回归算法更适合作为当前问题下的时空融合预测器。由于本文所使用的原始数据集是精确到小时粒度的，且数据缺失程度较高，我们也未能找到其它可以准确参考的小时级数据源，因此本文利用天级粒度的气象后报与PM2.5浓度后报对所用数据集进行手工近似填补，试图尽量排除缺失值带来的不良影响。前两组实验均是在手工近似填补后的数据集的基础上进行的。第3组实验利用相同的预测模型，分别基于手工近似填补缺失值的数据集和时空融合填补缺失值的数据集进行预测，旨在验证本文所提出的基于时空融合填补缺失值的有效性。

4.1 实验数据说明

本文以北京地区为研究对象对PM2.5浓度进行预测。北京共有16个行政区县和35个空气质量监测站。本文所使用的数据源可以分为3类：各区县气象数据、各空气质量监测站PM2.5浓度数据、各空气质量监测站经纬度数据。根据第3类数据可以得到站点间的邻居关系。前两组数据的时间跨度为2014年5月1日至2015年4月30日，共计12个月，均以小时为单位进行记录。但由于缺失数据的存在，前两组数据的样本量无法达到小时级粒度下应有的数量。对于以下实验，我们把前11个月(2014年5月至2015年3月)的数据作为训练集，把第12个月(2015年4月)的数据作为测试集，对35个站点未来一小时的PM2.5浓度进行预测。

4.2 评价指标

本课题属于回归类问题。因此，本文选用了回归问题常用的两种评价指标：均方根误差(RMSE)和决定系数(R2)，对实验效果进行衡量。这两个评价指标的公式如式(3)、式(4)所示

(3)

(4)

4.3 时空融合与单维度预测的效果对比

本文第2部分提出了基于时空融合的PM2.5预测模型。该模型分别使用LSTM和ANN从时空两个角度独立预测，再将预测结果融合起来用模型树回归。为论证该方法的有效性，本实验将时空融合预测(LSTM+ANN+MT)的效果与时间维度单模型预测(LSTM)和空间维度单模型预测(ANN)的效果进行对比。

表1对时空融合预测模型的整体表现进行了统计。可以看出经过时空融合预测后35个站点的平均指标均好于时空单独预测。且利用LSTM+ANN+MT进行预测后，大部分站点比时间或空间单维度预测的最好效果还有提升。可见，本文提出的时空融合预测模型是有效的。

4.4 不同时空融合预测器的效果对比

本文2.3部分介绍了基于时空融合的MT预测器。为证明MT作为时空融合预测器的有效性，对于相同的数据集，本文选择了回归算法中较为常用的线性回归(LR)和支持向量回归(SVR)分别替代MT作为时空融合预测器。

表1 融合预测与单独预测效果统计

表2对3个模型的整体表现进行了统计。可以看出MT作为时空融合预测模型时，35个站点的平均RMSE和R2均好于LR和SVR。可见，MT在回归前先对数据集分组，然后分别在每个类中进行回归，这种方式可以更好地适应样本间差异较大的数据集，使每次回归考虑的样本量缩小，提升回归算法的针对性和拟合效果。因此利用MT作为时空融合预测模型是有效的。

表2 不同时空融合预测器效果统计

另外，即使是使用LR或SVR作为时空融合预测器，与表1对比，35个站点的平均指标也都优于时空单独预测，进一步验证了本文提出的基于时空融合的预测框架的有效性。

4.5 两种缺失值填补方法的效果对比

本文利用第3部分提出的算法，对数据集中各属性分别进行缺失值填补。填补前后的数据集缺失程度对比可见表3。从中可以看出，缺失值填补后，缺失值大量减少，可为模型训练提供更加充足的样本。其中，填补后数据集依然不完整是由于3.4中所提到的4个参考值均为缺失值的原因所导致的。

此外，对于相同的原始数据集，本文选用两种方式对缺失值进行填补：利用天级后报数据进行手工近似填补(下称“手工填补”)、本文提出的基于时空融合进行缺失值填补(下称“时空填补”)。对于这两个填补后的数据集，本文利用相同的时空融合预测模型(LSTM+ANN+MT)进行预测，并从RMSE和R2的角度对预测效果进行对比，如图6、图7所示。其中，灰色标记代表手工近似填补，黑色标记代表时空融合填补。并以不同形状对不同预测方式进行区分。