一种基于时空分布特征的空气质量指数预测方法

摘要

本发明涉及一种空气质量指数预测方法，尤其涉及一种基于时空分布特征的空气质量指数预测方法。本方法采用时间序列预测方法预测时间方向，即未来某个时间的空气质量指数；再采用克里格插值法，利用已知的监测站点经纬度坐标，结合时间序列预测结果，插值估计整个区域任意地点的空气质量指数。因此，本发明具有如下优点：1.降低了模型的复杂度，并且缩短了整体的计算时间，同时又保证了模型的准确性。2.可以更准确地预测区域内每个地方、一天内多个时段的空气质量指数。

说明书

技术领域

本发明涉及一种空气质量指数预测方法，尤其涉及一种基于时空分布 特征的空气质量指数预测方法。

背景技术

空气是地球上的生物赖以生存的物质，是必不可少的一种物质。环境 空气质量与人们的日常生活息息相关，同时也在城市环境综合评价中占有 重要地位。但是随着人类文明和经济的发展，空气污染越来越严重，如何 改善空气质量、合理进行大气环境质量预测预警变得越来越重要，根据当 地的空气质量预测结果，人们可以采取相应措施如带口罩，尽量避免外出 等，保护自己免受空气污染物的侵害。

空气质量指数用AQI表示，是一个用来定量描述空气质量水平的数值。 它由空气中二氧化硫(SO2)、二氧化氮(NO2)、一氧化碳(CO)、臭氧(O3)、 悬浮颗粒物PM10和PM2.5，共6中污染物的含量计算出来的。AQI的取值 范围位于0–500之间。由于空气质量情况与气候因素、地理位置有关， 在数据角度上对应的就是空气质量指数与时间、空间分布都相关。因此需 要有一定的技术方法，利用区域内多个监测站观测并计算得到的空气质量 指数数据，预测整个地区未来的空气质量指数。

但是，传统的空气质量预测方法存在以下问题：

(1)不能准确的对没有观测站的区域进行空气指数预测：一般利用各 地的多个空气质量监测站历史的检测数据，来预测该站点(或附近)未来 一段时间的空气质量指数，而不能预测没有观测站的地方未来的空气质量 情况；也有一种时空克里格插值方法可以仅根据离散分布的观测点的历史 数据资料，预测整个区域范围内某个时间的空气质量指数，但是时空克里 格插值法是内插法，即当要预测的未来时间点超出统计资料样本数据的最 大时间时，预测误差将会变得很大。因此研究出新的方法可以比较准确的 预测整个区域内每个地方在未来某个时间的空气质量指数是有必要的。

(2)不能进行细粒度时间段的空气指数预测：传统的空气质量预测通 常是对接下来一整天的空气质量进行预测，这种粗粒度时间段的空气质量 预测存在一个弊端，假设明天的空气质量预测结果是优，表明对空气污染 物较敏感的人群可以进行户外活动，而实际上某些时间段(如上午 8:00-9:00，下午5:00-6:00)的空气质量有可能是中度污染或者重度污染， 对空气污染较敏感人群在这些时间段不宜进行户外活动。因此，对接下来 一天的空气质量进行细粒度时间段的预测更为合理。

发明内容

本发明主要是解决现有技术所存在的不能准确的对没有观测站的区域 进行空气指数预测以及不能进行细粒度时间段的空气指数预测的技术问题； 提供了一种基于时空分布特征的空气质量指数预测方法。该方法将空间克 里格插值法与时间序列预测方法相结合，把空间插值估计过程和时间序列 预测过程分开进行，降低了模型的复杂度，并且缩短了整体的计算时间， 同时又保证了模型的准确性；并且本方法可以更准确地预测区域内每个地 方、一天内多个时段的空气质量指数，满足人们挑选空气质量较好的时间 段进行外出活动的需要。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种基于时空分布特征的空气质量指数预测方法，包括以下步骤：

样本获取步骤，用于获取待预测区域中的空气质量监测站的经纬度坐 标，每隔一定时间段记录区域中的各个空气质量监测站在同一时刻监测得 到的空气质量指数，得到空气质量样本数据；

时间序列预测步骤，用于根据样本获取步骤中得到空气质量样本数据 建立各个空气质量监测站的时间序列函数，选择置信度较高的时间序列函 数预测该函数所对应的空气质量检测站在未来某个时间采样点的第一监测 点空气指数预测集；

克里格插值步骤，用于根据第一监测点空气指数据预测集拟合得到待 预测区域中的空间变异函数，根据所述空间变异函数以及第一监测点空气 指数预测集进行克里格插值估计。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述时间序列预测步骤进一步包括以下子步骤：

时间序列划分子步骤，用于将空气质量样本数据划分为m个时间序列， 每个时间序列包括来自同一监测站点的n个空气质量样本，其中：m是空气 质量监测站的个数，n是空气质量监测站的数据记录次数；

模型训练子步骤，对于每个时间序列，根据已有样本中n个时间，选 择最后的ι个时间对应的ι×m个样本作为验证样本，即实际值real_ij(i＝ 1,2,…，m；j＝1,2,…，ι)；选择前n-ι个数据作为时间序列训练样本，利 用训练样本训练得到各个站点的移动平均自回归模型，利用移动平均自回 归模型预测验证样本中的ι个时间对应的空气质量指数,记为z_ij(i＝ 1,2,…，m；j＝1,2,…，ι)，分别计算各个z_ij值的80％置信区间和90％的置 信区间；

模型验证子步骤，对m站点的时间序列m_i(i＝1,2,…，m)，若其ι个 实际值分别落在z_ij的80％置信区间和90％的置信区间内部的数量比例分别 达到预设的第一阈值和预设的第二阈值，则将该时间序列m_i对应的移动平 均自回归模型保留，共得到s个观测点的移动平均自回归模型；

时间序列预测子步骤，利用保留的s个观测点的移动平均自回归模型 预测未来t时刻该s个站点的空气质量指数Z(x_i)(i＝1,2,…,s)，将Z(x_i)作 为第一监测点空气指数预测集。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述第一阈值为85％，所述第二阈值为75％。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述克里格插值步骤进一步包括以下子步骤：

区域划分步骤，用于将待预测区域划分成不相交的正方形网格；

变异函数拟合子步骤，利用第一监测点空气指数预测集以及相应的监 测点坐标，并基于以下公式计算得到待预测区域内的变异函数离散值：

${\gamma }^{*}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\Sigma }}{[Z\left(x_i\right)-Z(x_i+h)]}^2$

式中，h为两监测点空间分隔距离，N(h)表示距离为h的监测点的对 数，Z(xi)和Z(xi+h)分别监测点xi和监测点xi+h处第一监测点空气指 数预测值，其中：i∈[1,2,...,N(h)]；选择适当的圆滑曲线来对变异函数 离散值进行拟合，将样本均方差值最小的函数作为待预测区域内的变异函 数；

范围搜索子步骤，利用拟合得到的变异函数的变程a，搜索对网格中心 处x0有影响的监测点，x0处经纬度坐标为(x₀,y₀),计算第一监测点空气指 数预测集中的s个站点与x0的距离di(i＝1,2,…,s),取di小于等于a的所 有的有影响的观测点作为有效观测点构成有效观测点集合s0；

加权系数计算子步骤，利用变异函数、观测点集合s0并基于以下公式 计算得到有效范围内的有效观测点集s0的空气质量指数在对某个网格X0 进行克里格估值计算时的权值λ_i：

$\left(\begin{array}{c}\underset{j=1}{\overset{s_0}{\Sigma }}{\lambda }_j\gamma (x_i,x_j)+\mu =\gamma (x_i,x_0),i=\mathrm{1,2},...,s_0\\ \underset{i=1}{\overset{s_0}{\Sigma }}{\lambda }_i=1\end{array}\right)$

式中:,

γ(x_i,x_j)为有效观测点x_i与x_j之间的变异函数值,γ(x_i,x₀)为待估点x₀与 有效观测点x_i之间的变异函数值；μ为拉格朗日乘数；

空间克里格估值子步骤，利用有效观测集合s0对应的第一监测点空气 质量指数预测值Z(x_i)(i＝1,2,…,s0)，并基于以下公式对网格X₀进行预测：

$Z\left(x_0\right)=\underset{i=1}{\overset{s_0}{\Sigma }}{\lambda }_{i·}Z\left(x_i\right)$

式中，Z(xi)是网格X0的有效观测点集s0对应的第一监测点空气质量 指数预测值，λi是s0个有效观测点对网格x0权值；对区域内其余的每个 网格，重复范围搜索子步骤、加权系数计算子步骤和空间克里格估值子步 骤，最后获得区域内各个网格的空气质量指数预测值。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述变异函数拟合子步骤中变异函数的拟合基于以下公式：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(\frac{3h}{2a}-\frac{h^3}{2a^3})& 0<h\le a\\ C_0+C& h>a\end{array}\right).$

式中，C0为块金常数，(C0+C)为基台值，C为拱高，a为变程。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述变异函数拟合子步骤中变异函数的拟合基于以下公式：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(1-e^{-\frac{h}{a}})& h>0\end{array}\right)$

式中，C0为块金常数，C为拱高，但a为三分之一的变程。

优化的，上述的一种基于时空分布特征的空气质量指数预测方法，所 述变异函数拟合子步骤中变异函数的拟合基于以下公式：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(1-e^{-\frac{h^2}{a^2}})& h>0\end{array}\right)$

式中，C0为块金常数，C为拱高，a为三分之一的变程。

因此，本发明具有如下优点：

1.本方法将空间克里格插值法与时间序列预测方法相结合，把空间插 值估计过程和时间序列预测过程分开进行，降低了模型的复杂度，并且缩 短了整体的计算时间，同时又保证了模型的准确性。

2.本方法可以更准确地预测区域内每个地方、一天内多个时段的空气 质量指数，满足人们挑选空气质量较好的时间段进行外出活动的需要。

附图说明

附图1是本发明的流程图；

附图2是本发明的时间序列的移动平均自回归模型预测图；

附图3是本发明拟合的球形变异函数曲线。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的 说明。

实施例：

图1是本发明的基于时空分布特征的空气质量指数预测方法。该方法 的主要包括三个大的步骤，分别是样本获取步骤、时间序列预测步骤、克 里格插值步骤，下面分别进行介绍。

1、样本获取

获取待预测区域中的空气质量监测站的经纬度坐标，每隔一定时间段 记录区域中的各个空气质量监测站在同一时刻监测得到的空气质量指数， 得到空气质量样本数据；本实施例中，各采样点之间的时间间隔优选为6 小时，为保证训练样本完整性，监测站点有数据历史数据有缺失情况，填 补缺失记录，获取n×m个样本。

2、时间序列预测

根据样本获取步骤中得到空气质量样本数据建立各个空气质量监测站 的时间序列函数，选择置信度较高的时间序列函数预测该函数所对应的空 气质量检测站在未来某个时间采样点的第一监测点空气指数预测集；该步 骤进一步包括以下子步骤：

2.1时间序列划分

用于将空气质量样本数据划分为m个时间序列，每个时间序列包括来 自同一监测站点的n个空气质量样本，其中：m是空气质量监测站的个数， n是空气质量监测站的数据记录次数；

2.2模型训练

对于每个时间序列，根据已有样本中n个时间，选择最后的ι个时间对 应的ι×m个样本作为验证样本，即实际值real_ij(i＝1,2,…，m；j＝ 1,2,…，ι)；选择前n-ι个数据作为时间序列训练样本，利用训练样本训练 得到各个站点的移动平均自回归模型作为该站点的时间序列函数。

移动平均自回归模型,简称为ARIMA(p,d,q)模型，是目前时间序列的众 多模型中比用通用，且预测效果较佳的模型。该模型含有p,d,q三个参数， 通过调整3个参数可保证时间序列的平稳性。因此对以上m个时间序列， 采用该模型时p,d,q三个参数各不相同。

利用训练得到的移动平均自回归模型预测验证样本中的ι个时间对应的 空气质量指数,记为z_ij(i＝1,2,…，m；j＝1,2,…，ι)，其中，z_ij的80％ 置信区间上下限分别为80_lower_ij和80_upper_ij,95％置信区间上下限分别 为95_lower_ij和95_upper_ij(i＝1,2,…，m；j＝1,2,…，ι)。

2.3模型验证

对站点的时间序列m_i(i＝1,2,…，m)，计算ι个实际值分别落在80％ 置信区间内部的数量比例80_ratio_i和落在95％置信区间内部的数量比例 95_ratio_i：

当80_ratio_i>＝80％，并且95_ratio_i>＝75％,就认为时间序列m_i(i＝ 1,2,…，m)的移动平均自回归模型ARIMA的预测精度较高，并且对应的观测 点予以保留，m个观测点中保留的有s个，并记录这s个观测点的经纬度坐 标为(sx_i,sy_i)i＝1,2,…,s)。

如图2所示，是一个时间序列的ARIMA模型的预测值、置信区间和实 际值的比较图，此ARIMA模型的3个参数分别是p＝1,d＝q＝0，右边阴影部分 中的线条是预测结果，阴影部分分别表示80％置信区间和95％置信区间，圆 圈代表实际值，此图比较直观的展示了ARIMA模型的预测效果。

2.4时间序列预测

利用保留的s个观测点的移动平均自回归模型预测未来t时刻该s个 站点的空气质量指数Z(x_i)(i＝1,2,…,s)，将Z(x_i)作为第一监测点空气指 数预测集。

3、克里格插值

根据第一监测点空气指数据预测集拟合得到待预测区域中的空间变异 函数，根据所述空间变异函数以及第一监测点空气指数预测集进行克里格 插值估计。

在知道区域内离散分布的观测点的同一时刻的空气质量指数时，要估 计整个区域任意地方的空气质量指数。由于空气质量情况是有空间分布特 征的，因此空气质量指数也有空间分布特征，对于有空间分布特征的数据， 采用克里格插值法可以无偏估计出区域内任意点的空气质量指数。主要步 骤如下：

3.1区域划分

将待预测区域划分成不相交的正方形网格，在同一时刻每个网格内的 各处的空气质量指数相等。克里格插值过程中网格单元越小，插值估计经 度也越高，因此在划分网格时，网格的面积很小，使每个网格内部的空气 质量指数近似相等的。

区域G的经纬度的最大最小值记为：经度最小值xmin,经度最大值 xmax，纬度最小值ymin，纬度最大值ymax。每个网格g的经纬跨度记为cell， 那么按照cell为单元来划分，区域G的经度可以划分为 XX＝(xmax-xmin)/cell份，其中XX取整数，经度方向可得到 xmin,xmin+cell,xmin+2*cell,…，xmin+XX*cell序列；区域G的纬度可 以划分为YY＝(ymax-ymin)/cell份，其中YY取整数。

纬度方向可划分为ymin,ymin+cell,ymin+2*cell,…，ymin+YY*cell 序列；将这些经纬度序列进行组合，总共有N＝XX*YY个经纬度坐标，即空 间坐标。这些坐标依次记为(x_i,y_i)(i＝1,2,…，N),每个坐标点是网格g的 中心，最终目标就是要预测每个格网g的中心点处的空气质量指数。

3.2变异函数拟合

空间变异的意思是，变量在空间中随着空间位置的变化而不同；变异函 数是描述这种变异性质的表达式。

利用第一监测点空气指数预测集以及相应的监测点坐标，并基于以下 公式计算得到待预测区域内的变异函数离散值：

${\gamma }^{*}\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\Sigma }}{[Z\left(x_i\right)-Z(x_i+h)]}^2$

式中，h为两监测点空间分隔距离，N(h)表示距离为h的监测点的对 数，Z(xi)和Z(xi+h)分别监测点xi和监测点xi+h处第一监测点空气指 数预测值，其中：i∈[1,2,...,N(h)]；

为了考察空间位置x处的值Z(x)受多大范围内的空间点值的影响，所 以要对上面离散的变异函数γ^*h进行拟合，采用加权多项式回归拟合法进 行拟合，拟合的光滑的变异函数曲线，与空间点位置xi无关，只与空间点 距离h有关。

本实施例中的变异函数的拟合可采用以下几种方法：

(1)球形变异函数模型：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(\frac{3h}{2a}-\frac{h^3}{2a^3})& 0<h\le a\\ C_0+C& h>a\end{array}\right)$

式中：C0为块金常数，(C0+C)为基台值，C为拱高，a为变程。当空间 点的距离h>＝变程a时，变异函数γ(h)值不再变化，即以点x处为中心，半 径为a范围内的点xi处的值Z(xi)对x处的值Z(x)有影响，超出这个范围 的无影响。也就是说克里格差值过程中，只选择与所预测的网格距离小于 等于变程的观测点。

采用球形变异函数模型进行拟合的结果如图3所示。

(2)指数模型：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(1-e^{-\frac{h}{a}})& h>0\end{array}\right)$

式中：C0，C意义同前，但a不是变程，由于1-e-3＝1-0.05＝0.95≈1， 则变程为3a

(3)高斯模型：

$\gamma \left(h\right)=\left(\begin{array}{cc}0& h=0\\ C_0+C(1-e^{-\frac{h^2}{a^2}})& h>0\end{array}\right)$

式中：C0，C意义同前，但a不是变程，由于1-e-3＝1-0.05＝0.95≈1， 则变程为3a。

3.3范围搜索

利用拟合得到的变异函数的变程a，搜索对网格中心处x0有影响的监 测点，x0处经纬度坐标为(x₀,y₀),计算第一监测点空气指数预测集中的s 个站点(经纬度坐标为(sx_i,sy_i))与x0的距离di(i＝1,2,…,s),取di小 于等于a的所有的有影响的观测点作为有效观测点构成有效观测点集合s0。

3.4加权系数计算

利用变异函数、观测点集合s0并基于以下公式计算得到有效范围内的 有效观测点集s0的空气质量指数在对某个网格X0进行克里格估值计算时 的权值λ_i：

$\left(\begin{array}{c}\underset{j=1}{\overset{s_0}{\Sigma }}{\lambda }_j\gamma (x_i,x_j)+\mu =\gamma (x_i,x_0),i=\mathrm{1,2},...,s_0\\ \underset{i=1}{\overset{s_0}{\Sigma }}{\lambda }_i=1\end{array}\right)$

式中:γ(x_i,x_j)为有效观测点x_i与x_j之间的变异函数值,γ(x_i,x₀)为待估点 x₀与有效观测点x_i之间的变异函数值；μ为拉格朗日乘数；此公式仅为置 信度高且在网格的有效范围内的监测站对一个网格的权值计算，对每个网 格的权值需重复计算该步骤，其中拉格朗日常数，也称为朗格朗日乘子， 是一种寻找变量受一个或多个条件所限制的多元函数的极值的方法。这种 方法将一个有n个变量与k个约束条件的最优化问题转换为一个有n+k 个变量的方程组的极值问题，因而可用求解极值方法求解条件方程组。

3.5空间克里格估值

利用有效观测集合s0对应的第一监测点空气质量指数预测值Z(x_i) (i＝1,2,…,s0)，并基于以下公式对网格X₀进行预测：

$Z\left(x_0\right)={\Sigma }_{i=1}^{s_0}{\lambda }_{i·}Z\left(x_i\right),$

式中，Z(xi)是网格X0的有效观测点集s0对应的第一监测点空气质量 指数第一监测点空气质量指数预测值，λi是s0个有效观测点对网格x0权 值；对区域内其余的每个网格，重复范围搜索子步骤、加权系数计算子步 骤和空间克里格估值子步骤，最后获得区域内各个网格的空气质量指数预 测值。

通过以上方案，本方法采用空间克里格插值法与时间序列预测方法相结 合，把空间插值估计过程和时间序列预测过程分开进行，降低了模型的复 杂度，并且缩短了整体的计算时间。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明 所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或 补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权 利要求书所定义的范围。

尽管本文较多地使用了等术语，但并不排除使用其它术语的可能性。 使用这些术语仅仅是为了更方便地描述和解释本发明的本质；把它们解释 成任何一种附加的限制都是与本发明精神相违背的。