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

摘要

本发明属于多学科和技术相互交叉结合的应用领域，具体涉及一种基于时空数据模型的城市级PM2.5浓度预测方法。首先利用开发的在线数据采集系统采集城市级PM2.5浓度数据；综合考虑模型可解释性、数据易获得和模型复杂度，确定最佳的模型预测因子，进而确定模型类型和模型结构；再根据实际采集PM2.5浓度数据的缺失情况和确定的模型结构，确定对缺失PM2.5浓度数据的插值算法；最后根据模型预测精度，确定能进行实时在线预测的算法。本发明能够充分利用随手可得的数据，确定模型结构，进而建立有效的城市级PM2.5浓度预测模型。

主权项

一种基于时空数据模型的城市级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)的时空数据模型：yk＝α0+αTxk+εk  (1)xk＝[tik‑1,tik‑2,…,tik‑d,w1k‑1,w1k‑2,…,w1k‑d,w2k‑1,w2k‑2,…,w2k‑d,w3k‑1,w3k‑2,…,w3k‑d,s1k‑1,s1k‑2,…,s1k‑d,s2k‑1,s2k‑2,…,s2k‑d,s3k‑1,s3k‑2,…s3k‑d]T  (2)时空数据模型中的xk∈R7m，具体形式如公式(2)，xk表示k时刻的预测因子，yk表示k时刻A地区PM2.5浓度值，d表示延时，εk表示k时刻的白噪声项，α0和α表示模型未知参数；当模型中的xk只含有时间变量ti时，模型称为自回归模型；2)模型参数拟合使用最小二乘算法对公式(1)的时空数据模型参数α0和α进行估计，令β＝(α0,αT)T，最小二乘算法估计参数如公式(3)所示：β＝(XTX)‑1XTY  (3)公式(3)中X＝[1,x]，1表示元素全部为1的列向量，x表示由输入变量的数据构成的矩阵；Y表示输出变量的数据构成的向量；根据拟合的α0和α参数，根据公式(1)计算出拟合值步骤4：最优模型结构选择，具体过程如下：使用公式(4)的赤池信息量准则AIC选择最优模型结构：AIC=n×ln(Σi=1n((yi-y^i)2)n)+2×p---(4)]]>其中：yi表示i时刻的观测值，表示i时刻的预测值，p表示模型参数个数，n表示样本个数；AIC的表达式分为两部分，第一部分意味着模型的精度，当很小时，第一部分值很小，AIC也小；第二部分意味着模型的复杂度，当p小时，模型结构简单，AIC也小；因此，AIC用于平衡模型精度和模型复杂度，选择使AIC最小的模型结构为最优模型结构；使用AIC选择最优模型结构的伪代码如下所示：步骤5：期望最大化算法EM对缺失PM2.5浓度值进行插值，具体步骤如下：EM算法用于含有隐变量的概率参数模型的最大似然和极大后验概率估计，使用EM算法进行插值计算时，将缺失数据设为隐变量，并根据似然函数最大化，求出模型的参数表达式；最后通过实际数据进行迭代，求出时空数据模型的参数α0和α，并使用参数α0和α计算当前的缺失值；EM算法推导出参数α0、α和方差σ2的过程如下所示：PM2.5浓度观测值：yi～N(μ,σ2)，其中i＝1,2,3…m，N(μ,σ2)表示均值为μ，方差为σ2的高斯分布，观测值构成的聚合为Y；PM2.5浓度缺失值：zi0，其中i0＝1,2,3…r，缺失值构成的集合为Z；PM2.5浓度完整数据集：T＝Y∪Z，∪表示并集，样本总数为n＝m+r，ti1表示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)的对数似然函数:logL=Σ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(logL)=n2(log2πσ2)-12σ2Σi=1m(yi-μi)2-12σ2Σj=1rE((zj-μj)2)---(7)]]>对公式(7)中的E((zj‑μj)2)进行变形和展开，得到公式(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函数对σ2求偏导：∂Q∂σ2=-n2σ2+12(σ2)2[Σi=1m(yi-βTxi)2+rσ2(k)+Σj=1r(β(k)Txj-βTxj)2]---(11)]]>令公式(11)等于零，解得σ2如公式(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)]]>综上所述：σ2和β的值如公式(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时刻预测值，yi为i时刻观测值；步骤6：利用遗忘因子递推辨识算法对实时采集的数据进行在线预测，并使用决定系数R2对计算结果进行分析，具体步骤如下：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)]]>其中：(xk+1,yk+1)为新数据；P为数据协方差矩阵；h为增益矩阵；α(0＜α＜1)为遗忘因子，α一般由人工确定，取值在0.95～0.99；2)R2介绍R2用于对模型预测效果进行拟合优度分析，R2的具体公式如下所示：R2=1-Σ(yi-y^i)2Σ(yi-y‾)2---(23)]]>其中，yi表示实际观测值，表示i时刻的预测值，表示所有观测样本y的均值；从公式知：R2取值在0到1之间，对于确定的观测样本，的值不变，当时，R2越大，模型的预测效果好。