CN107133686A - 基于时空数据模型的城市级pm2.5浓度预测方法 - Google Patents
基于时空数据模型的城市级pm2.5浓度预测方法 Download PDFInfo
- Publication number
- CN107133686A CN107133686A CN201710198680.6A CN201710198680A CN107133686A CN 107133686 A CN107133686 A CN 107133686A CN 201710198680 A CN201710198680 A CN 201710198680A CN 107133686 A CN107133686 A CN 107133686A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msup
- msub
- sigma
- munderover
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000013499 data model Methods 0.000 title claims abstract description 36
- 238000000034 method Methods 0.000 title claims abstract description 19
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 50
- 241001269238 Data Species 0.000 claims abstract description 19
- 230000006870 function Effects 0.000 claims description 24
- 238000004364 calculation method Methods 0.000 claims description 19
- 230000000694 effects Effects 0.000 claims description 18
- 238000004458 analytical method Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 238000007476 Maximum Likelihood Methods 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 abstract description 8
- 238000012217 deletion Methods 0.000 abstract 1
- 230000037430 deletion Effects 0.000 abstract 1
- 238000011160 research Methods 0.000 description 3
- 239000000443 aerosol Substances 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 229910002651 NO3 Inorganic materials 0.000 description 1
- NHNBFGGVMKEFGY-UHFFFAOYSA-N Nitrate Chemical compound [O-][N+]([O-])=O NHNBFGGVMKEFGY-UHFFFAOYSA-N 0.000 description 1
- QAOWNCQODCNURD-UHFFFAOYSA-L Sulfate Chemical compound [O-]S([O-])(=O)=O QAOWNCQODCNURD-UHFFFAOYSA-L 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 239000003570 air Substances 0.000 description 1
- 208000030303 breathing problems Diseases 0.000 description 1
- 235000009508 confectionery Nutrition 0.000 description 1
- 238000007596 consolidation process Methods 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 239000010419 fine particle Substances 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 239000013618 particulate matter Substances 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 210000002345 respiratory system Anatomy 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
Landscapes
- Business, Economics & Management (AREA)
- Engineering & Computer Science (AREA)
- Strategic Management (AREA)
- Economics (AREA)
- Tourism & Hospitality (AREA)
- Human Resources & Organizations (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Marketing (AREA)
- General Physics & Mathematics (AREA)
- General Business, Economics & Management (AREA)
- Quality & Reliability (AREA)
- Game Theory and Decision Science (AREA)
- Entrepreneurship & Innovation (AREA)
- Operations Research (AREA)
- Development Economics (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
- Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
Abstract
本发明属于多学科和技术相互交叉结合的应用领域,具体涉及一种基于时空数据模型的城市级PM2.5浓度预测方法。首先利用开发的在线数据采集系统采集城市级PM2.5浓度数据;综合考虑模型可解释性、数据易获得和模型复杂度,确定最佳的模型预测因子,进而确定模型类型和模型结构;再根据实际采集PM2.5浓度数据的缺失情况和确定的模型结构,确定对缺失PM2.5浓度数据的插值算法;最后根据模型预测精度,确定能进行实时在线预测的算法。本发明能够充分利用随手可得的数据,确定模型结构,进而建立有效的城市级PM2.5浓度预测模型。
Description
技术领域
本发明属于多学科和技术相互交叉结合的应用领域,具体涉及一种基于时空数据模型的城市级PM2.5浓度预测方法。
背景技术
PM2.5是对空气中直径小于或等于2.5微米的固体颗粒或液滴的总称。PM2.5能降低空气的能见度,导致灰霾天的出现;PM2.5通过呼吸道进入人的肺部,增加呼吸道疾病和心血管疾病的发病率。为了控制PM2.5对环境和人类健康造成的危害,需要实时掌握现在和将来的PM2.5浓度值。因此,对PM2.5浓度值进行实时预测具有重要的现实意义。
由于无法获得PM2.5浓度值的历史数据,开发了基于互联网的在线PM2.5浓度数据采集系统。在线PM2.5浓度数据采集系统获取的国内城市级PM2.5浓度值存在部分缺失,缺失数据无法用于PM2.5浓度的数据分析与建模研究。因此,需要对缺失的PM2.5浓度数据进行插值计算,有研究者使用了反向距离加权插值(Sree Dhevi,A.T."Imputing missingvalues using Inverse Distance Weighted Interpolation for time series data."Advanced Computing(ICoAC),2014Sixth International Conference on.IEEE,2014),文中使用了缺失数据的采样时刻到周围数据采样时刻的时间距离作为权值进行插值计算,反向距离加权插值只能对单个缺失数据进行插值,忽略了数据连续缺失的情况;也有使用土地使用回归模型和克里格方法进行插值(Zou B,Luo Y,Wan N,et al.Performancecomparison of LUR and OK in PM2.5concentration mapping:a multidimensionalperspective[J].Scientific reports,2015,5),对土地使用回归模型和克里格插值的插值效果进行了比较,说明土地使用回归模型的插值效果好,土地使用回归模型忽略了PM2.5浓度的时间和空间特性。
上述对缺失PM2.5浓度的插值算法存在不足:反向距离加权插值适合对单个缺失数据进行插值计算,而现实采集的数据存在连续缺失的情况,使得反向距离加权插值的应用受到限制;土地使用回归模型使用土地使用类型和人口密度作为预测因子,忽略了PM2.5浓度的时间和空间特性,使得插值的效果不好。对缺失PM2.5浓度数据的插值计算是预测PM2.5浓度的前提,是保证PM2.5浓度预测准确度的重要环节。
过去对PM2.5浓度预测的研究中,大多使用了线性模型,把研究的重点放在了预测因子的选择上,有学者使用湿度作为模型预测因子(Barman S C,Singh R,Negi M P S,etal.Fine particles(PM2.5)in residential areas of Lucknow city and factorsinfluencing the concentration[J].CLEAN–Soil,Air,Water,2008,36(1):111-117.),建立湿度与PM2.5浓度的线性模型,模型最终的拟合度只有0.53;有使用土地类型、人口密度作为PM2.5的预测因子(Zou B,Luo Y,WanN,et al.Performance comparison of LUR andOK in PM2.5concentration mapping:a multidimensional perspective[J].Scientificreports,2015,5),建立土地使用回归模型对PM2.5浓度进行预测,预测效果有所提升,拟合度为0.69;Song W和Jia H等人的文章中,在温度、湿度、风速的基础上,加入了边界层高度和卫星获取的气溶胶光学厚度作为PM2.5的预测因子(Song W,Jia H,Huang J,et al.Asatellite-based geographically weighted regression model for regional PM2.5estimation over the Pearl River Delta region in China[J].Remote Sensing ofEnvironment,2014,154:1-7),建立地理加权回归模型,模型对测试集的拟合度为0.74。因此,在城市级PM2.5浓度数据集获取非常困难、数据采集站点少以及遥感技术获得的数据集精确度不高的条件下,如何选择适当的预测因子成为预测PM2.5浓度的关键。
上述预测PM2.5浓度的模型存在不足:Barman S C,Singh R等人最早使用气象变量中的湿度作为模型的预测因子,模型的预测因子过于单一;考虑到PM2.5浓度受到PM2.5采集站点的地理环境的影响,Zou B,Luo Y等人使用了土地使用类型和人口密度作为模型的预测因子,该模型的预测精度不高,以及需要的数据比较难获取;近年来,随着技术不断进步,遥感技术得到广泛使用,Song W和Jia H等人使用卫星获取的气溶胶光学厚度对PM2.5进行预测,预测效果得到了明显的提升,但由于遥感技术对图像采集技术要求很高,使得获取的数据准确度有待提升。因此,在选择模型预测因子的过程中,要考虑获取数据的难易程度和数据的准确性、PM2.5浓度自身的时间和空间特性、以及预测因子影响PM2.5浓度的原因。
发明内容
本发明要解决的技术问题有两个方面:一是如何选择适当的算法对缺失数据进行插值计算;另一方面是如何选择预测因子,建立模型实现PM2.5浓度的高精度在线预测。为了解决这两个问题,应该考虑如何充分利用随手可得的数据实现插值计算和预测。针对模型的预测因子选择问题,首先考虑到温度会导致硝酸盐、硫酸盐等化学物质发生复杂的化学反应,形成二次PM2.5;湿度会导致微小的颗粒物凝聚起来,从而使PM2.5浓度升高;风速会导致PM2.5的扩散,影响PM2.5浓度,因此,温度、湿度和风速对PM2.5浓度的影响很大;其次,考虑到PM2.5浓度具有时间特性,即当前时刻的PM2.5浓度受到过去时刻PM2.5浓度的影响;PM2.5浓度具有空间特性,即PM2.5浓度容易由高浓度区向低浓度区扩散。考虑到本发明的适用性,在选择预测因子时,要充分考虑数据获取的难易程度,考虑到模型的可解释性,尽量选择能解释影响PM2.5浓度原因的预测因子。确定模型的输入变量为过去PM2.5浓度值,温度,湿度,风速和空间变量,模型的输出为将来时刻的PM2.5浓度值。通过输入变量和输出变量,建立时空数据模型,并根据AIC对模型结构进行选择。在选择的模型结构基础上,选择机器学习的EM算法对缺失数据进行插值计算,首先根据模型结构推导出时空数据模型的参数,然后使用模型参数计算出缺失数据。为了实现对PM2.5浓度的实时预测,选择可以根据新数据随时更新模型参数的RFF算法。利用该发明可以对缺失的数据进行插值计算和对将来的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
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时,模型称为自回归模型(Auto regression,AR模型)。
2)模型参数拟合
使用最小二乘算法对公式(1)的时空数据模型参数α0和α进行估计,令β=(α0,αT)T,最小二乘算法估计参数如公式(3)所示:
β=(XTX)-1XTY (3)
公式(3)中X=[1,x],1表示元素全部为1的列向量,x表示由输入变量的数据构成的矩阵;Y表示输出变量的数据构成的向量;
根据拟合的α0和α参数,根据公式(1)计算出拟合值
步骤4:最优模型结构选择,具体过程如下:
使用公式(4)的赤池信息量准则(Akaike information criterion,AIC)选择最优模型结构:
其中:yi表示i时刻的观测值,表示i时刻的预测值,p表示模型参数个数,n表示样本个数;AIC的表达式分为两部分,第一部分意味着模型的精度,当很小时,第一部分值很小,AIC也小;第二部分意味着模型的复杂度,当p小时,模型结构简单,AIC也小。因此,AIC可以平衡模型精度和模型复杂度,选择使AIC最小的模型结构为最优模型结构。
使用AIC选择最优模型结构的伪代码如下所示:
步骤5:期望最大化算法(Expectation Maximization Algorithm,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)所示:
对公式(5)两边取对数,得到公式(6)的对数似然函数:
对公式(6)取期望,得到公式(7)的Q函数:
对公式(7)中的E((zj-μj)2)进行变形和展开,得到公式(8):
E((zj-μj)2)=E((zj-μ(k)+μ(k)-μj)2) (8)
=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)的结果代入公式(7)的Q函数中:
EM算法的E步是指期望Q函数如公式(9)所示,在E步的基础上,然后推导M步,以Q函数为目标函数,求目标函数的极大值或最大值,即M步转为公式(10)
目标函数:J=maxQ (10)
公式(9)的Q函数对σ2求偏导:
令公式(11)等于零,解得σ2如公式(12)所示:
公式(10)的目标函数对β求偏导:
令公式(13)等于零,解得β如公式(14)所示:
综上所述:σ2和β的值如公式(15)所示:
使用EM算法进行参数求解的伪代码如下所示:
为了对插值的效果进行分析,使用均方根误差(root-mean-square error,RMSE)作为分析的指标,均方根误差的公式如公式(19)所示:
其中n表示样本数,为i时刻预测值,yi为i时刻观测值。由于均方根误差对数据中的特大或特小误差反映非常敏感,可以很好地反映插值方法的插值效果。
步骤6:利用遗忘因子递推辨识算法(Recursive forgetting factor,RFF算法)对实时采集的数据进行在线预测,并使用决定系数(R2)对计算结果进行分析,具体步骤如下:
1)RFF算法
RFF算法由最小二乘算法演变而来,主要是为了克服数据饱和现象,基本思想是对旧数据加遗忘因子,降低旧数据信息的占有量。RFF算法适用于在线递推计算,公式为(20),(21),(22):
其中:(xk+1,yk+1)为新数据;P为数据协方差矩阵;h为增益矩阵;α(0<α<1)为遗忘因子,α一般由人工确定,取值在0.95~0.99之间为宜。
2)R2介绍
R2用于对模型预测效果进行拟合优度分析。R2的具体公式如下所示:
其中,yi表示实际观测值,表示i时刻的预测值,表示所有观测样本y的均值;从公式可知:R2取值在0到1之间,对于确定的观测样本,的值不变,当时,R2越大,模型的预测效果好。
本发明的效果和益处是:
本发明在建立在线数据采集系统时,充分考虑到当前国内公布的PM2.5浓度数据是以城市为单位,没有考虑城市各区域的PM2.5浓度值,在线数据采集系统可以实时采集城市各站点的PM2.5浓度值。在选择预测因子时,考虑了数据获取的简便性和模型的可解释性,结合预测因子影响PM2.5浓度的原因和PM2.5浓度自身具有的特性。在选择模型结构时,考虑了模型精度和模型复杂度之间的关系,即模型结构过于简单和过于复杂,都会降低模型的预测精度。考虑到缺失PM2.5浓度数据受过去PM2.5浓度数据和过去气象数据的影响,使用选择的最优模型结构,并结合可以进行批量数据处理的算法进行插值计算。考虑到实时预测的精度,使用实变的算法进行在线预测,从而便于模型参数的更新。
本发明可以充分利用采集的城市级PM2.5浓度和气象数据,对建立的时空数据模型进行模型结构选择和插值计算,最终根据完整数据集实现PM2.5浓度的实时在线预测,为雾霾天气的治理和人类活动提供可靠的依据。
附图说明
图1为模型结构选择图。
图2(a)为EM算法的缺失数据插值效果图。
图2(b)为回归法的缺失数据插值效果图。
图3为RFF算法实时在线预测效果图。
具体实施方式
下面结合附图对本发明的实施实例进行说明,另注意,此处描述的优选实施实例仅仅用于解释本发明,并不限定。
为了更好地理解本发明的技术方案,以大连市甘井子区的PM2.5浓度预测为例,对本发明的实施方式作进一步描述。大连市分布的PM2.5浓度采集站点有甘井子、周水子、开发区、青泥洼桥、傅家庄、七贤岭、星海三站和旅顺。
本发明的具体实施步骤如下:
步骤1:使用在线数据采集系统实时采集大连市八个站点每小时的PM2.5浓度值,并获取足够多的样本数据;从中国气象科学数据共享服务网中获取大连市温度、湿度和风速;
步骤2:对采集的大连市八个站点PM2.5浓度数据进行相关性分析,从而获得时空数据模型的空间变量。使用Matlab软件计算出的相关系数如表2所示。
表2大连市八个站点的PM2.5浓度相关系数表
从表可知,周水子PM2.5浓度、开发区PM2.5浓度、青泥洼桥PM2.5浓度和甘井子PM2.5浓度的相关系数值最大,分别为0.8643、0.8222、0.8711。因此,选择周水子、开发区和青泥洼桥三个区域的PM2.5浓度数据作为模型的空间变量。
步骤3:建立时空数据模型,具体步骤如下:
3)建立时空数据模型
时空数据模型的输入、输出变量如表3所示:
表3时空数据模型输入、输出变量表
根据输入输出变量,建立公式(24)的时空数据模型:
yk=α0+αTxk+εk (24)
使用公式(3)的最小二乘算法对公式(24)的时空数据模型参数α0和α进行估计,根据拟合的α0和α参数,根据公式(24)计算出拟合值
步骤4:最优模型结构选择,具体过程如下:
使用公式(4)的AIC选择时空数据模型的最优模型结构,伪代码如下所示:
选取没有缺失的PM2.5浓度数据进行模型选择,使用Matlab进行仿真,模型选择的结果如图1所示。图1中的AIC基本趋势是先减小后上升,最后趋于平稳;在整个过程中,AR模型的AIC都要高于时空数据模型的AIC,说明空间变量和气象变量是选择的最优变量。在时空数据的AIC中,延时为4时,AIC取得最小值。因此,选择的最优模型为时空数据模型,模型的延时为4。
步骤5:EM算法插值,具体步骤过程如下:
EM算法是由Dempster、Rubin等人于1977年提出的一种简单实用的迭代算法,用于含有隐变量的概率参数模型的最大似然和极大后验概率估计。使用EM算法进行插值计算时,将缺失数据设为隐变量,并根据似然函数最大化,求出模型的参数表达式;最后通过实际数据进行迭代,求出时空数据模型的参数α0和α,并使用参数α0和α计算当前的缺失值。使用公式(25)计算参数σ2和β:
使用EM算法进行参数求解的伪代码如下所示:
使用公式(29)的均方根误差分析EM算法的插值效果:
为了对插值算法进行验证,使用人为缺失的数据集作为样本数据,样本数据的缺失情况包含了单个数据缺失和连续数据缺失两种情况,使用EM算法、以气象数据和过去PM2.5浓度数据作为外部变量的回归模型对样本数据进行插值计算,插值效果如图2(a)和2(b)所示,两种插值方法的RMSE如表4所示。
表4插值效果统计表
对比两种插值方法的插值效果,可以得出EM插值算法的插值效果好。回归模型使用气象变量和过去时刻的PM2.5浓度作为预测因子,使用建立的模型对缺失数据进行预测,而模型的参数与训练集有关,且随着数据的不断增多,数据带的新信息无法在模型中得到积累,使得模型的预测效果逐渐变弱,从图2(b)的回归方法插值效果可知,回归模型对前面的缺失数据插值效果好,随着数据的增多,插值效果逐渐变差,并出现了负值;EM算法使用回归模型作为基础,并根据批量数据对模型参数进行了调整,加入了新信息的内容,使得EM算法的插值效果优于回归模型法。
步骤6:利用RFF算法对实时采集的数据进行在线预测,并使用R2对计算结果进行分析,具体步骤如下:
1)使用公式(30),(31),(32)的RFF算法对甘井子区PM2.5浓度数据进行在线预测:
遗传因子α取值为0.98。
2)使用公式(33)的R2对RFF算法的预测效果进行拟合优度分析:
3)在线计算和结果分析
使用Matlab对采集的数据进行实时在线预测,结果如图3所示。图3中的样本数为219,在线计算的R2为0.8403,与具有高影响度和代表性的论文的预测效果(表5)进行对比,模型拟合度高;为了进一步说明在线计算的效果好,使用气象数据作为预测因子,对同一数据集进行了在线计算,两个模型的计算效果对比情况如表6所示,表6中的Err表示实际值和真实值的误差平方和,1.0e4表示Err的单位是10的4次方。通过与表5和表6进行对比,说明本发明的实时在线预测效果精度高。
表5代表论文预测效果统计表
表6模型效果统计表
Claims (1)
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)的时空数据模型:
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选择最优模型结构:
<mrow>
<mi>A</mi>
<mi>I</mi>
<mi>C</mi>
<mo>=</mo>
<mi>n</mi>
<mo>&times;</mo>
<mi>l</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>y</mi>
<mo>^</mo>
</mover>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
</mrow>
<mi>n</mi>
</mfrac>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>2</mn>
<mo>&times;</mo>
<mi>p</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
其中: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)所示:
<mrow>
<mi>L</mi>
<mo>=</mo>
<munderover>
<mo>&Pi;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mfrac>
<mn>1</mn>
<msqrt>
<mrow>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mfrac>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>t</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
</msup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
对公式(5)两边取对数,得到公式(6)的对数似然函数:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>log</mi>
<mi> </mi>
<mi>L</mi>
<mo>=</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mi>log</mi>
<mfrac>
<mn>1</mn>
<msqrt>
<mrow>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mfrac>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>t</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mrow>
<mi>log</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>t</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfrac>
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mrow>
<mi>log</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>t</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfrac>
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mrow>
<mi>log</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
对公式(6)取期望,得到公式(7)的Q函数:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>Q</mi>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<mrow>
<mi>log</mi>
<mi> </mi>
<mi>L</mi>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfrac>
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
对公式(7)中的E((zj-μj)2)进行变形和展开,得到公式(8):
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>+</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>2</mn>
<mi>E</mi>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mrow>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>z</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
将公式(8)的结果代入公式(7)的Q函数中:
<mrow>
<mi>Q</mi>
<mo>=</mo>
<mfrac>
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&lsqb;</mo>
<msup>
<mi>r&sigma;</mi>
<mrow>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&mu;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<mo>-</mo>
<msub>
<mi>&mu;</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>=</mo>
<mfrac>
<mi>n</mi>
<mn>2</mn>
</mfrac>
<mrow>
<mo>(</mo>
<mi>l</mi>
<mi>o</mi>
<mi>g</mi>
<mn>2</mn>
<msup>
<mi>&pi;&sigma;</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&lsqb;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>r&sigma;</mi>
<mrow>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mi>T</mi>
</mrow>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
</mrow>
EM算法的E步是指期望Q函数如公式(9)所示,在E步的基础上,然后推导M步,以Q函数为目标函数,求目标函数的极大值或最大值,即M步转为公式(10)
目标函数:J=maxQ (10)
公式(9)的Q函数对σ2求偏导:
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>Q</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mi>n</mi>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>+</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>&lsqb;</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>r&sigma;</mi>
<mrow>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mi>T</mi>
</mrow>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
令公式(11)等于零,解得σ2如公式(12)所示:
<mrow>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mo>&lsqb;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>r&sigma;</mi>
<mrow>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
<mi>T</mi>
</mrow>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>&rsqb;</mo>
<mo>/</mo>
<mi>n</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>12</mn>
<mo>)</mo>
</mrow>
</mrow>
公式(10)的目标函数对β求偏导:
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>Q</mi>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&beta;</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<mo>-</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<mo>-</mo>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&beta;</mi>
<mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mi>T</mi>
</mrow>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>13</mn>
<mo>)</mo>
</mrow>
</mrow>
令公式(13)等于零,解得β如公式(14)所示:
<mrow>
<mi>&beta;</mi>
<mo>=</mo>
<mo>&lsqb;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>&rsqb;</mo>
<mo>/</mo>
<mo>&lsqb;</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>14</mn>
<mo>)</mo>
</mrow>
</mrow>
综上所述:σ2和β的值如公式(15)所示:
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>&beta;</mi>
<mo>=</mo>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mi>&beta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msup>
<msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mo>/</mo>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
<mo>+</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mi>&sigma;</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mrow>
<mo>&lsqb;</mo>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>m</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>i</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mi>r&sigma;</mi>
<mrow>
<mn>2</mn>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
<mo>+</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>r</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<msup>
<mi>&beta;</mi>
<mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mi>T</mi>
</mrow>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>x</mi>
<mi>j</mi>
</msub>
</mrow>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<mo>/</mo>
<mi>n</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
使用EM算法进行参数求解的伪代码如下所示:
为了对插值的效果进行分析,使用均方根误差RMSE作为分析的指标,均方根误差的公式如公式(19)所示:
<mrow>
<mi>R</mi>
<mi>M</mi>
<mi>S</mi>
<mi>E</mi>
<mo>=</mo>
<msqrt>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>y</mi>
<mo>^</mo>
</mover>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
<mi>n</mi>
</mfrac>
</msqrt>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>19</mn>
<mo>)</mo>
</mrow>
</mrow>
其中n表示样本数,为i时刻预测值,yi为i时刻观测值;
步骤6:利用遗忘因子递推辨识算法对实时采集的数据进行在线预测,并使用决定系数R2对计算结果进行分析,具体步骤如下:
1)RFF算法
RFF算法适用于在线递推计算,公式为(20),(21),(22):
<mrow>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>P</mi>
<mi>k</mi>
</msub>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
<mrow>
<mi>&alpha;</mi>
<mo>+</mo>
<msup>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mi>T</mi>
</msup>
<msub>
<mi>P</mi>
<mi>k</mi>
</msub>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>20</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mover>
<mi>&beta;</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mover>
<mi>&beta;</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&lsqb;</mo>
<msub>
<mi>y</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>-</mo>
<msup>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mi>T</mi>
</msup>
<mover>
<mi>&beta;</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>21</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msub>
<mi>P</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>&alpha;</mi>
</mfrac>
<mo>&lsqb;</mo>
<mi>I</mi>
<mo>-</mo>
<mi>h</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<msup>
<msub>
<mi>x</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mi>T</mi>
</msup>
<mo>&rsqb;</mo>
<msub>
<mi>P</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>22</mn>
<mo>)</mo>
</mrow>
</mrow>
其中:(xk+1,yk+1)为新数据;P为数据协方差矩阵;h为增益矩阵;α(0<α<1)为遗忘因子,α一般由人工确定,取值在0.95~0.99;
2)R2介绍
R2用于对模型预测效果进行拟合优度分析,R2的具体公式如下所示:
<mrow>
<msup>
<mi>R</mi>
<mn>2</mn>
</msup>
<mo>=</mo>
<mn>1</mn>
<mo>-</mo>
<mfrac>
<mrow>
<mi>&Sigma;</mi>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>y</mi>
<mo>^</mo>
</mover>
<mi>i</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mi>&Sigma;</mi>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>y</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<mover>
<mi>y</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>23</mn>
<mo>)</mo>
</mrow>
</mrow>
其中,yi表示实际观测值,表示i时刻的预测值,表示所有观测样本y的均值;从公式知:R2取值在0到1之间,对于确定的观测样本,的值不变,当时,R2越大,模型的预测效果好。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710198680.6A CN107133686A (zh) | 2017-03-30 | 2017-03-30 | 基于时空数据模型的城市级pm2.5浓度预测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710198680.6A CN107133686A (zh) | 2017-03-30 | 2017-03-30 | 基于时空数据模型的城市级pm2.5浓度预测方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN107133686A true CN107133686A (zh) | 2017-09-05 |
Family
ID=59715910
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710198680.6A Pending CN107133686A (zh) | 2017-03-30 | 2017-03-30 | 基于时空数据模型的城市级pm2.5浓度预测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107133686A (zh) |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108120661A (zh) * | 2017-12-19 | 2018-06-05 | 北京理工大学 | 一种城市空气中颗粒物含量时空分布测定方法 |
CN109472366A (zh) * | 2018-11-01 | 2019-03-15 | 郑州云海信息技术有限公司 | 一种机器学习模型的编码解码方法与装置 |
CN109541730A (zh) * | 2018-11-23 | 2019-03-29 | 长三角环境气象预报预警中心(上海市环境气象中心) | 一种大气污染物浓度预测的方法及设备 |
CN110046771A (zh) * | 2019-04-25 | 2019-07-23 | 河南工业大学 | 一种pm2.5浓度预测方法与装置 |
CN110346420A (zh) * | 2019-06-09 | 2019-10-18 | 重庆工商大学融智学院 | 一种时空数据智能聚合方法 |
CN110738354A (zh) * | 2019-09-18 | 2020-01-31 | 北京建筑大学 | 预测颗粒物浓度的方法、装置、存储介质及电子设备 |
CN111401605A (zh) * | 2020-02-17 | 2020-07-10 | 北京石油化工学院 | 大气污染的可解释预测方法 |
CN111859054A (zh) * | 2020-07-23 | 2020-10-30 | 中国科学院计算机网络信息中心 | 气象卫星数据的处理方法及装置 |
CN115876964A (zh) * | 2023-01-31 | 2023-03-31 | 北方工业大学 | 一种城市街区气候环境与碳排放移动监测预警方法及系统 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103201754A (zh) * | 2010-11-18 | 2013-07-10 | 索尼公司 | 数据处理设备、数据处理方法以及程序 |
-
2017
- 2017-03-30 CN CN201710198680.6A patent/CN107133686A/zh active Pending
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103201754A (zh) * | 2010-11-18 | 2013-07-10 | 索尼公司 | 数据处理设备、数据处理方法以及程序 |
Non-Patent Citations (2)
Title |
---|
胡隽等: "随机缺失数据下的时间序列分析建模", 《数学的实践与认识》 * |
陈丽等: "基于时空数据模型的PM2.5预测", 《仪器仪表学报》 * |
Cited By (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108120661A (zh) * | 2017-12-19 | 2018-06-05 | 北京理工大学 | 一种城市空气中颗粒物含量时空分布测定方法 |
CN109472366B (zh) * | 2018-11-01 | 2020-07-24 | 苏州浪潮智能科技有限公司 | 一种机器学习模型的编码解码方法与装置 |
CN109472366A (zh) * | 2018-11-01 | 2019-03-15 | 郑州云海信息技术有限公司 | 一种机器学习模型的编码解码方法与装置 |
CN109541730A (zh) * | 2018-11-23 | 2019-03-29 | 长三角环境气象预报预警中心(上海市环境气象中心) | 一种大气污染物浓度预测的方法及设备 |
CN110046771A (zh) * | 2019-04-25 | 2019-07-23 | 河南工业大学 | 一种pm2.5浓度预测方法与装置 |
CN110346420A (zh) * | 2019-06-09 | 2019-10-18 | 重庆工商大学融智学院 | 一种时空数据智能聚合方法 |
CN110346420B (zh) * | 2019-06-09 | 2022-03-01 | 重庆工商大学融智学院 | 一种时空数据智能聚合方法 |
CN110738354B (zh) * | 2019-09-18 | 2021-02-05 | 北京建筑大学 | 预测颗粒物浓度的方法、装置、存储介质及电子设备 |
CN110738354A (zh) * | 2019-09-18 | 2020-01-31 | 北京建筑大学 | 预测颗粒物浓度的方法、装置、存储介质及电子设备 |
CN111401605A (zh) * | 2020-02-17 | 2020-07-10 | 北京石油化工学院 | 大气污染的可解释预测方法 |
CN111401605B (zh) * | 2020-02-17 | 2023-05-02 | 北京石油化工学院 | 大气污染的可解释预测方法 |
CN111859054A (zh) * | 2020-07-23 | 2020-10-30 | 中国科学院计算机网络信息中心 | 气象卫星数据的处理方法及装置 |
CN111859054B (zh) * | 2020-07-23 | 2023-12-26 | 中国科学院计算机网络信息中心 | 气象卫星数据的处理方法及装置 |
CN115876964A (zh) * | 2023-01-31 | 2023-03-31 | 北方工业大学 | 一种城市街区气候环境与碳排放移动监测预警方法及系统 |
CN115876964B (zh) * | 2023-01-31 | 2024-01-23 | 北方工业大学 | 一种城市街区气候环境与碳排放移动监测预警方法及系统 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107133686A (zh) | 基于时空数据模型的城市级pm2.5浓度预测方法 | |
Kannan et al. | A nonparametric kernel regression model for downscaling multisite daily precipitation in the Mahanadi basin | |
Kuligowski et al. | Experiments in short-term precipitation forecasting using artificial neural networks | |
Heine et al. | Development and comparison of approaches for automated mapping of stream channel networks | |
CN108009674A (zh) | 基于cnn和lstm融合神经网络的空气pm2.5浓度预测方法 | |
CN110232471B (zh) | 一种降水传感网节点布局优化方法及装置 | |
Kohail et al. | Implementation of data mining techniques for meteorological data analysis | |
Huang et al. | An analytical comparison of four approaches to modelling the daily variability of solar irradiance using meteorological records | |
CN112506990B (zh) | 一种基于时空信息的水文数据异常检测方法 | |
CN107292098A (zh) | 基于前期气象因子与数据挖掘技术的中长期径流预报方法 | |
Solaiman et al. | Development of probability based intensity-duration-frequency curves under climate change | |
Sartini et al. | Comparing different extreme wave analysis models for wave climate assessment along the Italian coast | |
CN106231609A (zh) | 一种基于重点目标区域的水下传感器网络优化部署方法 | |
CN107798425A (zh) | 一种基于大数据的时空混淆暴露度评估系统及方法 | |
CN112287294B (zh) | 一种基于深度学习的时空双向土壤含水量插值方法 | |
Olsson et al. | Statistical atmospheric downscaling of short-term extreme rainfall by neural networks | |
CN105869100A (zh) | 一种基于大数据思维的滑坡多场监测数据的融合及预测方法 | |
CN115689293B (zh) | 一种基于压力-状态-响应框架的城市内涝韧性评估方法 | |
CN115220133B (zh) | 一种多气象要素降雨预测方法、装置、设备及存储介质 | |
Huang et al. | Research on urban modern architectural art based on artificial intelligence and GIS image recognition system | |
CN113836808A (zh) | 一种基于重污染特征约束的pm2.5深度学习预测方法 | |
CN114186491A (zh) | 基于改进lur模型的细颗粒物浓度时空特征分布方法 | |
CN114049545A (zh) | 一种基于点云体素的台风定强方法、系统、设备及介质 | |
CN110414088B (zh) | 结合水动力模型的涉禽栖息地适宜性空间模糊评价方法 | |
CN117332909B (zh) | 基于智能体的多尺度城市内涝道路交通暴露性预测方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
WD01 | Invention patent application deemed withdrawn after publication | ||
WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20170905 |