基于马尔科夫模型的太原市空气质量等级预测分析
2021-12-16夏桂梅邵升清
夏桂梅,邵升清,张 绚
(太原科技大学 应用科学学院,山西 太原 030024)
0 引言
空气质量是环境监测的一项重要指标,其质量好坏直接或间接地影响到了人们的身体健康、动植物的生长和发育以及生产生活中的各个领域。如果可以较为准确地预测出某地区的空气质量等级情况,就可以提前做好防护及相关应对措施,减少其带来的危害。目前,空气质量预测的统计模型有马尔科夫模型[1]、神经网络(BP)模型[2]、灰色系统GM(1,1)模型[3]等。有时也使用组合模型进行预测,如:Diaz-Robles等[4]将ARIMA(p,d,q)模型与神经网络模型组合来提高预测精度,Zheng 等[5]将神经网络和线性回归组合进行预测。本文基于太原市两年内的空气质量等级历史数据,应用马尔科夫模型的相关理论知识,构建模型并使用Matlab编程实现了对太原市的空气质量等级的预测。
1 原理与方法
1.1 马尔科夫过程
设X=(X1,X2,…,Xt)是随机变量序列,其中每个随机变量的取值在有限集S={s1,s2,…,sn},称为状态空间。Markov特征是:
1)有限历史假设
P(Xt+1=sk|X1,…,Xt)=P(Xt+1=sk|Xt)
(1)
2)时间不变性[6]假设
∀i∈{1,2,3,…,T}∀x,y∈S,P(Xi=y|Xi-1=x)=p(y|x)
(2)
如果X具有以上的这些特征,则这个X序列称为一个马尔科夫过程(链)。
(3)
(4)
P(n)=P(n-1)P
(5)
1.2 马尔科夫分析法
马尔科夫分析法的基本模型为:
X(K+1)=X(K)×P
(6)
式中:X(K)表示趋势分析与预测对象在T=K时刻的状态向量;P表示一步转移概率矩阵;X(K+1)表示趋势分析与预测对象在T=K+1时刻的状态向量。
1.3 构建马尔科夫模型的基本过程
1.3.1构造一步转移概率矩阵
设pij=P(Xi=j|X0=i),i,j∈I,则称pij为转移概率,它表示由状态i到状态j的转移概率,以pij为元素的转移概率矩阵[7]如下所示:
(7)
1.3.2马尔科夫性检验
数列具有马尔科夫性[8]是用马尔科夫模型进行预测的重要前提。所以要先对数据进行检验,使用χ2来对检验结果进行判断。假设指标数列分为m个状态,将状态i向状态j转化的频数记为fij,转移概率为pij,并将状态转移频数矩阵的条件概率当作边际概率。根据fij(i,j∈I)计算出矩阵每列和,再计算出所有元素的和,最后进行除法运算。计算公式如下:
(8)
计算χ2统计量的公式如下:
(9)
这里需要说明的是:式中的log是因为在编程中lnx常写为logx.
1.3.3构造多步转移概率矩阵
多步转移概率矩阵的计算需要用到C-K方程(chapman-colmogorov方程),计算公式如下:
P(u+v)=P(u)P(v)
(10)
在公式(10)中,令u=1,v=n-1,可以得到递推关系:P(n)=P(1)P(n-1)=PP(n-1)=…=Pn.故对齐次马氏链而言,n步转移概率矩阵就等于一步转移概率的n次方。
1.3.4计算不变概率测度
2 建模与检验
2.1 数据整理分类
选取太原市2019年4月1日至2021年3月31日的730条空气质量等级数据,数据来源于 “2345天气预报”网站(http://waptianqi.2345.com).根据环境空气质量标准,将我国空气质量分为6个等级。在预测过程中将其分别替换为表1的这6种状态,划分情况见表1.
表1 空气质量状态分类表Table 1 Table of air quality status classification
2.2 计算一步转移概率矩阵
对于730个空气质量等级状态数据,表2给出了太原市空气质量状态的转移情况。
表2 空气等级状态转移次数表Table 2 The number of air class status transitions
基于空气质量状态的转移情况,再用matlab计算,可得到一步转移概率矩阵P(1):
2.3 检验马尔科夫性
对太原市空气质量等级所选取的历史数据进行统计,可计算出转移频数矩阵(fij)6×6和它所对应的概率矩阵(pij)6×6,根据公式(8)可求得边际概率见表3,统计量χ2的值见表4.
表3 边际概率表Table 3 Table of marginal probability
表4 统计量χ2计算表Table 4 Calculation table of statistic χ2
2.4 计算转移概率矩阵的极限分布
根据C-K方程的公式知P(2)=P(1)×P(1)=P(1)2,P(3)=P(1)3,…,P(n)=P(1)n,计算多步转移概率矩阵,发现最终转移概率矩阵的极限分布如下:
(0.153 3,0.597 9,0.177 2,0.050 9,0.013 8,0.006 9)
2.5 计算不变概率测度
由不变概率测度的概念和公式可知πP=π,其中π=(π1,π2,π3,π4,π5,π6),P为一步转移概率矩阵P(1),从而我们可得到如下方程组:
即P的不变测度为(0.153 3,0.597 9,0.177 2,0.050 9,0.013 8,0.006 9).
2.6 空气质量等级预测与检验
以2021年4月15日的空气质量等级为初始状态,空气质量等级为中度污染,状态数为4,则初始状态向量设定为P(0)=(0,0,0,1,0,0),对后5 d的空气质量等级进行预测。预测结果见表5.
表5 空气质量等级状态预测Table 5 Prediction of air quality level status
2.7 预测结果检验
由表5知,4月15日后5 d的空气质量预测状态数分别为3、2、2、2、2,即空气质量等级为轻度污染、良、良、良、良的概率最大,而实际空气质量等级为良、良、良、中度污染、良,预测准确率为80%.
3 结果
利用P的不变测度概率(0.153 3,0.597 9,0.177 2,0.050 9,0.013 8,0.006 9)计算2020年这6种空气质量出现的天数,并与实际出现天数进行比较,检验模型的预测效果,结果见表6.
表6 空气质量等级状态预测Table 6 Prediction of air quality level status
2020年度空气质量这6个等级的预测天数分别为56 d、218 d、65 d、18 d、5 d和3 d,实际天数分别为54 d、223 d、65 d、15 d、8 d和0 d,相对误差为4.38%,模型预测效果较好。
4 结语
本文将样本数据(太原市2019年4月1日至2021年3月31日的730条空气污染指数)划分为6种状态,后利用样本数据建立马尔科夫模型对空气等级进行预测,并根据2020年空气质量的历史数据进行模型检验。结果表明,马尔科夫模型对空气等级进行预测的效果较好,人们能够利用该模型及时掌握空气质量变化情况,并根据预测数据采取相应措施应对环境变化[10]。