[发明专利]基于时空数据模型的城市级PM2.5浓度预测方法

权利要求书

1.一种基于时空数据模型的城市级PM2.5浓度预测方法，其特征在于，步骤如下：

步骤1：选择获取国内城市级PM2.5浓度数据的API接口，并根据选择的API接口建立在线数据采集系统；用在线数据采集系统实时采集所有待测站点每小时的PM2.5浓度值，得到样本数据；从中国气象科学数据共享服务网http://cdc.cma.gov.cn/中获取气象数据集，包含温度、湿度和风速；

步骤2：考虑时空数据模型可解释性和数据易获取，选择适当的预测因子；根据PM2.5浓度值的空间特性，通过分析各地区PM2.5浓度值之间的相关系数，选择城市其他区域的PM2.5浓度值作为空间变量；考虑到温度、湿度和风速对PM2.5浓度值的影响很大，选择温度、湿度和风速作为模型的外部变量；根据PM2.5浓度的时间特性，选择过去时刻的PM2.5浓度值作为时间变量；

步骤3：建立时空数据模型，具体步骤如下：

1)建立时空数据模型

时空数据模型的输入、输出变量如表1所示：

表1 时空数据模型输入、输出变量表

建立公式(1)的时空数据模型：

y_k＝α₀+α^Tx_k+ε_k(1)

x_k＝[ti_k-1,ti_k-2,…,ti_k-d,w1_k-1,w1_k-2,…,w1_k-d,w2_k-1,w2_k-2,…,w2_k-d,w3_k-1,w3_k-2,…,w3_k-d,

s1_k-1,s1_k-2,…,s1_k-d,s2_k-1,s2_k-2,…,s2_k-d,s3_k-1,s3_k-2,…s3_k-d]^T(2)

时空数据模型中的x_k∈R^7m，具体形式如公式(2)，x_k表示k时刻的预测因子，y_k表示k时刻A地区PM2.5浓度值，d表示延时，ε_k表示k时刻的白噪声项，α₀和α表示模型未知参数；当模型中的x_k只含有时间变量ti时，模型称为自回归模型；

2)模型参数拟合

使用最小二乘算法对公式(1)的时空数据模型参数α₀和α进行估计，令β＝(α₀,α^T)^T，最小二乘算法估计参数如公式(3)所示：

β＝(X^TX)^-1X^TY(3)

公式(3)中X＝[1,x]，1表示元素全部为1的列向量，x表示由输入变量的数据构成的矩阵；Y表示输出变量的数据构成的向量；

根据拟合的α₀和α参数，根据公式(1)计算出拟合值

步骤4：最优模型结构选择，具体过程如下：

使用公式(4)的赤池信息量准则AIC选择最优模型结构：

AIC=n×ln(Σi=1n((yi-y^i)2)n)+2×p---(4)]]>

其中：y_i表示i时刻的观测值，表示i时刻的预测值，p表示模型参数个数，n表示样本个数；AIC的表达式分为两部分，第一部分意味着模型的精度，当很小时，第一部分值很小，AIC也小；第二部分意味着模型的复杂度，当p小时，模型结构简单，AIC也小；因此，AIC用于平衡模型精度和模型复杂度，选择使AIC最小的模型结构为最优模型结构；

使用AIC选择最优模型结构的伪代码如下所示：

步骤5：期望最大化算法EM对缺失PM2.5浓度值进行插值，具体步骤如下：

EM算法用于含有隐变量的概率参数模型的最大似然和极大后验概率估计，使用EM算法进行插值计算时，将缺失数据设为隐变量，并根据似然函数最大化，求出模型的参数表达式；最后通过实际数据进行迭代，求出时空数据模型的参数α₀和α，并使用参数α₀和α计算当前的缺失值；EM算法推导出参数α₀、α和方差σ²的过程如下所示：

PM2.5浓度观测值：y_i～N(μ,σ²)，其中i＝1,2,3…m，N(μ,σ²)表示均值为μ，方差为σ²的高斯分布，观测值构成的聚合为Y；

PM2.5浓度缺失值：z_i0，其中i0＝1,2,3…r，缺失值构成的集合为Z；

PM2.5浓度完整数据集：T＝Y∪Z，∪表示并集，样本总数为n＝m+r，t_i1表示i1时刻的PM2.5浓度值；

其中，μ^(k)表示第k次的期望，β^(k)表示第k次的时空数据模型参数，β^T表示β的转置；

对公式(1)进行变形，得即模型转化为y＝β^Tx+ε的形式；

对模型y＝β^Tx+ε取期望，得E(y)＝E(β^Tx)+E(ε)，由于E(ε)＝0，则E(y)＝E(β^Tx)，由E(y)＝μ，得E(y)＝E(β^Tx)＝μ；

观测值y的似然函数L如公式(5)所示:

L=Πi=1n12πσ2e-(ti-μi)22σ2---(5)]]>

对公式(5)两边取对数，得到公式(6)的对数似然函数:

log L=Σi=1nlog12πσ2e-(ti-μi)22σ2=Σi=1n[-12(log2πσ2)-(ti-μi)22σ2]=n2(log2πσ2)-12σ2Σi=1n(ti-μi)2=n2(log2πσ2)-12σ2Σi=1m(yi-μi)2-12σ2Σj=1r(zj-μj)2---(6)]]>

对公式(6)取期望，得到公式(7)的Q函数:

Q=E(log L)=n2(log2πσ2)-12σ2Σi=1m(yi-μi)2-12σ2Σj=1rE((zj-μj)2)---(7)]]>

对公式(7)中的E((z_j-μ_j)²)进行变形和展开，得到公式(8)：

E((zj-μj)2)=E((zj-μ(k)+μ(k)-μj)2)=E((zj-μ(k))2)+E((μ(k)-μj)2)+2E(zj-μ(k))(μ(k)-μj)=E((zj-μ(k))2)+E((μ(k)-μj)2)+2(E(zj)-μ(k))(μ(k)-μj)=E((zj-μ(k))2)+E((μ(k)-μj)2)+2(μ(k)-μ(k))(μ(k)-μj)=E((zj-μ(k))2)+E((μ(k)-μj)2)---(8)]]>

将公式(8)的结果代入公式(7)的Q函数中：

Q=n2(log2πσ2)-12σ2Σi=1m(yi-βTxi)2-12σ2[rσ2(k)+Σj=1r(μ(k)-μj)2]]]>

=n2(log2πσ2)-12σ2[Σi=1m(yi-βTxi)2+rσ2(k)+Σj=1r(β(k)Txj-βTxj)2]---(9)]]>

EM算法的E步是指期望Q函数如公式(9)所示，在E步的基础上，然后推导M步，以Q函数为目标函数，求目标函数的极大值或最大值，即M步转为公式(10)

目标函数：J＝maxQ(10)

公式(9)的Q函数对σ²求偏导：

∂Q∂σ2=-n2σ2+12(σ2)2[Σi=1m(yi-βTxi)2+rσ2(k)+Σj=1r(β(k)Txj-βTxj)2]---(11)]]>

令公式(11)等于零，解得σ²如公式(12)所示：

σ2=[Σi=1m(yi-βTxi)2+rσ2(k)+Σj=1r(β(k)Txj-βTxj)2]/n---(12)]]>

公式(10)的目标函数对β求偏导：

∂Q∂β=-12σ2[Σi=1m-2(yi-βTxi)(xi)+Σj=1r-2(β(k)Txj-βTxj)(xj)]---(13)]]>

令公式(13)等于零，解得β如公式(14)所示：

β=[Σi=1myixi+Σj=1rβ(k)xjTxj]/[Σi=1mxiTxi+Σj=1rxjTxj]---(14)]]>

综上所述：σ²和β的值如公式(15)所示：

β=[Σi=1myixi+Σj=1rβ(k)xjTxj]/[Σi=1mxiTxi+Σj=1rxjTxj]σ2=[Σi=1m(yi-βTxi)2+rσ2(k)+Σj=1r(β(k)Txj-βTxj)2]/n---(15)]]>

使用EM算法进行参数求解的伪代码如下所示：

为了对插值的效果进行分析，使用均方根误差RMSE作为分析的指标，均方根误差的公式如公式(19)所示：

RMSE=Σi=1n(yi-y^i)2n---(19)]]>

其中n表示样本数，为i时刻预测值，y_i为i时刻观测值；

步骤6：利用遗忘因子递推辨识算法对实时采集的数据进行在线预测，并使用决定系数R²对计算结果进行分析，具体步骤如下：

1)RFF算法

RFF算法适用于在线递推计算，公式为(20)，(21)，(22)：

h(k+1)=Pkxk+1α+xk+1TPkxk+1---(20)]]>

β^(k+1)=β^(k)+h(k+1)[yk+1-xk+1Tβ^(k)]---(21)]]>

Pk+1=1α[I-h(k+1)xk+1T]Pk---(22)]]>

其中：(x_k+1,y_k+1)为新数据；P为数据协方差矩阵；h为增益矩阵；α(0＜α＜1)为遗忘因子，α一般由人工确定，取值在0.95～0.99；

2)R²介绍

R²用于对模型预测效果进行拟合优度分析，R²的具体公式如下所示：

R2=1-Σ(yi-y^i)2Σ(yi-y‾)2---(23)]]>

其中，y_i表示实际观测值，表示i时刻的预测值，表示所有观测样本y的均值；从公式知：R²取值在0到1之间，对于确定的观测样本，的值不变，当时，R²越大，模型的预测效果好。