2. 南京信息工程大学 计算机与软件学院, 南京 210044;
3. 南京信息工程大学 江苏省网络监控中心, 南京 210044
2. School of Computer and Software, Nanjing University of Information Science and Technology, Nanjing Jiangsu 210044, China;
3. Jiangsu Engineering Center of Network Monitoring, Nanjing University of Information and Technology, Nanjing Jiangsu 210044, China
霾是悬浮在大气中的颗粒物积累到一定程度,使水平能见度小于10 km的一种天气现象[1]。近年来,中国中东部霾事件频发,霾过程中低能见度会引发交通事故,高气溶胶浓度会影响空气质量和气候,也会引发多种呼吸道疾病,危害人类健康,受到各级政府、部门及社会各界的广泛关注。因此,急需建立和发展高精度的霾预报系统,进而加强霾的防御和防治。
目前霾预报主要有天气学预报法、数值预报法和统计学预报法。天气学预报法通过分析和诊断霾形成的天气学条件,结合预报员的经验和实况,用外推方法预报霾过程。廖碧婷等[2]利用1998—2008年观测和模式资料计算了垂直交换系数,并利用垂直交换系数对典型霾过程进行预报。但天气学预报法的预报水平和精细化程度远远满足不了用户对预报服务的需求。
数值预报法主要通过空气质量数值预报模式预报能见度和主要污染物浓度,从而预报霾过程。中国气象局自主研发的雾霾数值预报系统(CUACE/haze-fog)[3]对主要雾霾过程预报效果较好,但对过程中污染物峰值浓度及能见度预报存在一定偏差。珠江三角洲空气质量数值预报系统能基本预报出霾过程中污染物的产生、聚集和消亡过程和时间节点,但低估了污染物的浓度[4]。Wang等[5]建立的长三角空气质量预报系统预报的上海和南京空气污染指数的准确率分别为50%~83%和80%。北京区域环境气象数值预报系统(Beijing Regional Environmental Meteorology Prediction System, BREMPS)对2014年京津冀PM2.5浓度的预报效果较好,相关系数在0.6以上,但预报值总体低于观测值,且24 h之后预报效果略有下降[6]。由于现阶段模式中污染物排放、传输和沉降过程中的各参数均存在较大不确定性,导致预报的污染物浓度也存在较大不确定性[5]。
统计学预报法是通过分析霾天气出现时的大尺度天气系统活动规律,建立中低空大气温度、相对湿度和风速等气象要素和能见度、污染物浓度的统计关系,最终得到霾的统计预报模型。已有研究利用最小二乘法、回归分析、时间序列分析、遗传算法、神经网络、卡尔曼滤波法、小波分析、支持向量机等方法预报霾。邹乐强[7]基于气象要素、天气形势、能见度以及污染物浓度资料,建立了1~3 d霾集成预报系统,对有无霾预报准确率较高,但对重度霾的预报准确率偏低。高辉[8]基于K最近邻(K Nearest Neighbor, KNN)数据挖掘算法构建了北京地区霾等级预报模型,对轻度霾、中度霾和重度霾的空报率仅为4.7%、1.4%和2.6%。这些预测方法各有各的特点和适用场合, 但由于污染物浓度受排放量、气象条件、谱分布以及化学成分和性质等的共同影响,用单一统计方法建立准确度较高的统计预报模型难度较大。将各预报模型进行合理的组合,可以使各个模型的特点相互补充,达到更加准确的预报效果[9]。已有文献利用时间序列分析结合卡尔曼滤波应用在导航算法[10]、交通路段行程预测[11]和微机械陀螺仪漂移误差补偿[12]中,尚没有研究将时间序列分析结合卡尔曼滤波法应用至霾预报中。因此,本文选择时间序列分析和卡尔曼滤波进行混合建模,有利于降低预报误差。
综上所述,基于霾预报的重要性和难度,本文将采用时间序列分析方法,建立能见度预报模型,再结合卡尔曼滤波订正方法建立能见度预报订正方程,并排除降水、沙尘暴、扬沙、浮尘、吹雪、雪暴等天气现象,进而得到霾预报结果。本文旨在建立具有更高实时性和精度的霾预报模型,这将有助于提高霾预报能力,为政府防控提供科学参考。
1 基于时间序列和卡尔曼滤波霾预报模型 1.1 时间序列分析Box等[13]在20世纪70年代初提出了差分自回归滑动平均(AutoRegressive Integrated Moving Average, ARIMA)模型。它是指在将非平稳化的时间序列经过差分运算转化为平稳的时间序列,然后将预报量仅对它的滞后值以及随机误差的现值和滞后值进行差分自回归滑动平均的建模。其基本思想是用一定的数学模型来描述预报量随时间变化而形成的一组随机有序时间序列,通过时间序列建立的模型结合预报量的过去值和现在值预报未来值。ARIMA(p, d, q)模型由AR(p)模型、平稳化时间序列所做的差分次数d和MA(q)模型组成。AR(p)模型为自回归模型,其中p为自回归项;差分次数d是将时间序列由非平稳序时间序列通过差分运算转化成平稳时间序列所需的次数;MA(q)为滑动平均模型,其中q为滑动平均项数。
ARIMA模型建模的基本过程如下:
1) 分别求取该样本序列的自相关系数(Auto Correlation Function, ACF)和偏相关函数(Partial Auto Correlation Function, PACF),通过图检验法初步判断样本序列的平稳性。然后通过单位根检验法ADF来验证序列的平稳性。
2) 通过图检验法和ADF检验后,如果样本序列是非平稳的,那么需要对非平稳样本序列进行差分运算,直到非平稳序列转化称平稳序列为止。
3) 对经过平稳化后的样本序列使用贝叶斯信息量(Bayesian Information Criterion, BIC)对模型定阶。
4) 模型确定后,对模型参数进行估计,并且检验模型是否具有统计意义。
5) 检验残差序列是否为白噪声,进行假设性检验。
6) 使用通过检验的模型进行预报。
1.2 卡尔曼滤波的原理和应用卡尔曼滤波方法采用量测方程和状态方程组成的线性随机系统的状态空间模型来描述滤波器,并且利用状态方程的递推性,按线性无偏最小均方差估计准则,采用递推算法对滤波器的状态变量作出最佳估计,从而求得滤掉噪声后有用信号的最佳估计。大气环境预报领域中的卡尔曼滤波基本公式如下:
$ {\mathit{\boldsymbol{Y}}_t} = {\mathit{\boldsymbol{X}}_t}{\mathit{\boldsymbol{\beta }}_t} + {\mathit{\boldsymbol{v}}_t} $ | (1) |
$ {\mathit{\boldsymbol{\beta }}_t} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{t - 1}}{\mathit{\boldsymbol{\beta }}_{t - 1}} + {\mathit{\boldsymbol{w}}_{t - 1}} $ | (2) |
式(1)称为量测方程,式(2)称为状态方程。式中,βt=[β1, β2, …, βm]tT为卡尔曼滤波系统状态向量,也称回归系数;Φ为状态转移矩阵;Yt为测量向量即预报量,Yt=[Y1, Y2, …, Ym]tT;Xt为预报因子矩阵。假定Φ为单位矩阵,经过简化后的卡尔曼滤波状态方程为:
$ {\mathit{\boldsymbol{\beta }}_t} = {\mathit{\boldsymbol{\beta }}_{t - 1}} + {\mathit{\boldsymbol{v}}_{t - 1}} $ | (3) |
假定w和v是两个互不相关的随机向量,应用最小二乘法,即可得到卡尔曼滤波的递推方程组,即:
$ {\mathit{\boldsymbol{\hat Y}}_t} = {\mathit{\boldsymbol{X}}_t}{\mathit{\boldsymbol{\hat \beta }}_{t-1}} $ | (4) |
$ {\mathit{\boldsymbol{R}}_t} = {\mathit{\boldsymbol{C}}_{t-1}} + \mathit{\boldsymbol{W}} $ | (5) |
$ {\mathit{\boldsymbol{\sigma }}_t} = {\mathit{\boldsymbol{X}}_t}{\mathit{\boldsymbol{R}}_t}{\mathit{\boldsymbol{X}}}_{t}^{\rm{T}} + \mathit{\boldsymbol{V}} $ | (6) |
$ {\mathit{\boldsymbol{A}}_t} = {\mathit{\boldsymbol{R}}_t}{\mathit{\boldsymbol{X}}}_{t}^{\text{T}}{\mathit{\boldsymbol{\sigma }}}_{t}^{-1} $ | (7) |
$ \mathit{\boldsymbol{\hat \beta }} = {\mathit{\boldsymbol{\hat \beta }}_{t - 1}} + {\mathit{\boldsymbol{A}}_t}({\mathit{\boldsymbol{Y}}_t} - {\mathit{\boldsymbol{\hat Y}}_t}) $ | (8) |
$ {\mathit{\boldsymbol{C}}_t} = {\mathit{\boldsymbol{R}}_t} - {\mathit{\boldsymbol{A}}_t}{\mathit{\boldsymbol{\sigma }}_t}{\mathit{\boldsymbol{A}}}_{t}^{\rm{T}} $ | (9) |
式(4)~(9)组成的递推系统中:W和V分别表示动态噪声和量测噪声的方差阵,
卡尔曼滤波订正方法通过不断递推,不断更新实测信息,采用更新实测的数据对状态估计进行修正,从而对预报产品进行订正,提高预报准确率。
1.3 时序差分自回归滑动平均与卡尔曼滤波模型本方案中采用的资料包括:气象观测资料、欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts, ECMWF)细网格数值预报产品。在实验之前,采用双线性插值法,将ECMWF格点数据转化成站点数据作为实验数据。
本方案首先选用时序差分自回归滑动平均模型得到最优的能见度预报模型,然后结合使用卡尔曼滤波订正方法建立能见度客观预报订正模型,再排除降水、沙尘暴、扬沙、浮尘、吹雪、雪暴等视程障碍天气,得到霾预报结果。模型流程如图 1所示,具体步骤如下:
1) 选取4个典型代表站(北京、南京、广州、上海),根据对霾天气的分析,利用中国气象局实况资料,构建时序差分自回归滑动平均模型。
2) 将各代表站观测值作为样本序列,选取700个数据进行建模,使用图检验法和根检验法判断序列是否平稳,如果样本序列为非平稳序列,则对非平稳序列进行差分运算,直至样本序列平稳为止。
3) 分别计算经过平稳化的各代表站样本序列的自相关函数(ACF)和偏自相关函数(PACF)。根据BIC定阶并确定ARIMA模型,然后检查模型的有效性。确定混合模型流程(如图 1所示)。
4) 利用ECMWF模式产品,将ARIMA模型确定的方程模型作为卡尔曼滤波的初始量测方程,然后根据ECMWF模式产品确定W0、V0、C0、β0。由于本文使用时间序列分析建立能见度最优预报方程,所以将时间序列分析预报方程中的回归系数作为卡尔曼滤波过程中回归系数的初始值。由于β0是从样本资料中通过严格的数学计算得到的,因此理论上认为它是十分精确的,即与理论值相等,所以C0是m阶零方阵,即C0=[0]m×m。经过实验证明状态量测噪声V0和输入噪声W0的初值对回归系数β的订正并不敏感,因此,本文将W0、V0都设为零方阵。最后根据卡尔曼滤波递推公式,对时间序列建立的能见度预报方程动态修改权重。将卡尔曼滤波订正能见度预报方程的过程及误差订正过程设计成自动批处理程序。
5) 根据霾天气发生时能见度、相对湿度、降水等气象要素进一步判断,得到霾预报结果。最后利用观测实况进行预报检验。选用气象上衡量预报准确率的校正得分率(True Shooting percentage, TS)评分算法,对基于时间序列分析和卡尔曼滤波的霾预报结果进行评分。TS评分的定义如下:
$ TS = \frac{{{t_1}}}{{{t_1} + {t_2} + {t_3}}} \times 100\% $ | (10) |
其中:t1为正确预报霾的次数,t2为漏报霾的次数,t3为空报霾的次数。由定义看出,TS越大,预报效果就越好,预报准确率就越高。
2 实验结果 2.1 实验资料和平台本文选取北京、南京、杭州和广州4个典型站点的观测值作为实验数据。将观测值(逐小时能见度实况资料)作为样本序列{Zt},运用时间序列分析分别对4个典型站点进行建模。
选取北京站连续700 h的数据作为能见度样本时间序列{Zt}(如图 2(a)所示),然后求取北京站样本时间序列的前50个自相关系数和偏相关函数(图 2(b)、(c))。北京站能见度样本资料自相关函数大致按指数衰减至零,可以初步推断原始能见度时间序列是平稳序列。同时本文使用单位根检验法(ADF检验)来检验北京站能见度时间序列,原始序列平稳性检验结果:h=1(这里的h为平稳性检验标志,非零代表平稳,零代表非平稳)。检验结果与之前通过分析自相关系数图得出的结论是一致的,即北京站能见度样本时间序列是平稳数据,可以用来作为预报序列。
根据BIC确定ARIMA模型,最小的BIC值为3.561×103,所以p=1,q=2。至此p、q定阶完成。由此北京站能见度样本序列的时间序列模型为ARIMA(1,0,2)模型。由表 1模型参数计算获得该模型方程为:
$ X(t) = 0.933345X(t - 1) + {\alpha _t} + 0.0330181{\alpha _{t - 1}} - 0.191094{\alpha _{t - 2}} $ | (11) |
式中:αt为残差;0.033 018 1αt-1对预报方程影响较少,可以省去,所以公式如式(12):
$ X(t) = 0.933345X(t - 1) + {\alpha _t} - 0.191094{\alpha _{t - 2}} $ | (12) |
同理,选取南京站连续700 h的数据作为能见度样本时间序列{Zt}(如图 3(a)所示),求取的前50个自相关系数和偏相关函数(图 3(b)、(c))。南京站能见度样本资料自相关函数也很快衰减至零,初步推断原始能见度时间序列是平稳序列,进一步使用单位根检验法(ADF检验)检验出南京站能见度时间序列是平稳数据,与之前通过分析自相关系数图得出的结论是一致的,可以作为预报序列。
根据BIC确定ARIMA模型,最小的BIC值为2.661 4×103,所以p=5,q=4。至此p、q定阶完成。由此南京站能见度样本序列的时间序列模型为ARIMA(5,0,4)模型。由表 2模型参数计算获得该模型方程为式(13):
$ X(t) = 1.51105X(t - 1) + 0.388007X(t - 2){\rm{ - }} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;0.92432X(t - 3) - 0.596349X(t - 4) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;0.606433X(t - 5) + {\alpha _t} + 0.661273{\alpha _{t - 1}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;0.872461{\alpha _{t - 2}} - 0.228326{\alpha _{t - 3}} - 0.475871{\alpha _{t - 4}} \\ $ | (13) |
同理,选取杭州站连续700个小时的数据作为能见度样本时间序列{Zt},如图 4(a)所示,求取的前50个自相关系数和偏相关函数(图 4(b)、(c))。杭州站能见度样本资料自相关函数也很快衰减至零,初步推断原始能见度时间序列是平稳序列,进一步使用单位根检验法(ADF检验)检验出南京站能见度时间序列是平稳数据,可以作为预报序列。根据BIC准则确定ARIMA模型,最小的BIC值为2.496 7×103,所以p=2,q=1。所以至此p、q定阶完成。由此杭州站能见度样本序列的时间序列模型为ARIMA(2,0,1)模型。由表 3模型参数计算获得该模型方程为式(14):
$ X(t) = 1.62102X(t - 1){\rm{ - }}0.670331X(t - 2) + {\alpha _t} + 0.607047{\alpha _{t - 1}} $ | (14) |
同理,选取广州站连续700 h的数据作为样本时间序列{Zt},如图 5(a)所示,求取的前50个自相关系数和偏相关函数(图 5(b)、(c))。广州站能见度样本资料自相关函数也很快衰减至零,初步推断原始能见度时间序列是平稳序列,进一步使用单位根检验法(ADF检验)检验出广州站能见度时间序列是平稳数据,可以作为预报序列。根据BIC确定ARIMA模型,最小的BIC值为3.562×103,所以p=6,q=5。所以至此p, q定阶完成。由此广州站能见度样本序列的时间序列模型为ARIMA(6,0,5)模型。
由表 4模型参数计算获得该模型方程为式(15):
$ X(t) = 0.91651X(t - 1) + 1.00275X(t - 2) - \\ \;\;\;\;\;\;\;\;\;\;\;\;0.29851X(t - 3) - 0.756201X(t - 4) - \\ \;\;\;\;\;\;\;\;\;\;\;\;0.632804X(t - 5) + 0.738651X(t - 6) + \\ \;\;\;\;\;\;\;\;\;\;\;\;{\alpha _t} + 0.432253{\alpha _{t - 1}} + 1.25093{\alpha _{t - 2}} - 0.404755{\alpha _{t - 3}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;0.356343{\alpha _{t - 4}} - 0.123638{\alpha _{t - 5}} \\ $ | (15) |
随后将以上通过时序模型得到的北京站、南京站、杭州站和广州站预报方程用于卡尔曼滤波,通过卡尔曼滤波订正时序模型求得结果。
2.2 实验结果分析 2.2.1 本方法实验结果在确定了时间序列差分自回归滑动平均和卡尔曼滤波的霾客观预报订正模型后,做北京、南京、杭州、广州预报实验。实验日期:2014年10月— 2015年2月。实验结果如图 6所示。结果显示,北京、南京、杭州和广州4个站点的实验数据散点图中能见度客观预报订正结果与实况的波形图基本上一致,说明时间序列分析和卡尔曼滤波的霾客观预报模型能很好地反映实况能见度的变化趋势。将时间序列差分自回归滑动平均建立的能见度预报结果与经历过卡尔曼滤波订正之后的结果进行对比发现,经过卡尔曼滤波订正后的结果与能见度实况拟合情况更好。在实验过程中发现,如果只用能见度预报霾是否发生,则误差会很大,所以将能见度、风速、相对湿度和温度露点差结合建立霾的预报模型。结果显示:利用ECMWF模式预报因子数据和霾预报模型对2014年— 2015年2月进行预报实验来检验基于时间序列差分自回归滑动平均和卡尔曼滤波的霾预报模型的预报准确率。预报准确率较CUACE有明显提高。
利用本文的模型输出结果与中国气象局第一代雾霾数值预报业务系统(CUACE/haze-fog_54km)输出结果进行对比分析。选取2013年10月— 2014年2月各代表站点的霾天气进行预报实验,利用每日8:00的资料进行对比检验,检验结果见表 5。
由检验的对比结果可看出,利用时间序列差分自回归滑动平均与卡尔曼滤波建立的预报模型得到的霾预报结果与CUACE模式在北京、南京、杭州、广州的准确率相比分别提高了11%、13%、11%、1%(平均提高了10个百分点),而空报率和漏报率则有不同程度减低,平均值分别降低了3个百分点和7个百分点。进一步验证了本文的方法在霾预报上有着更好的效果。
3 结语为提高霾预报准确率,本文构建了时间序列差分自回归滑动平均和卡尔曼滤波混合的霾客观预报订正模型:先通过时间序列分析推导出能见度预报方程,然后利用卡尔曼滤波递推性动态订正预报过程中的权重,避免了时间序列建立低阶模型预报准确率低的情形,也降低了卡尔曼滤波在建立初始状态方程和测量方程的难度。该混合算法与CUACE模式模拟结果相比,能有效提高霾预报准确率。
目前本文研究内容仅限于北京、南京、杭州、广州几个典型站点,且客观预报准确率还有待提高,后续研究将增加实验站点、延长实验时间序列、改进算法,进一步提高客观预报准确率。
[1] | 吴兑. 关于霾与雾的区别和灰霾天气预警的讨论[J]. 气象, 2005, 31(4): 3-7. (WU D. Discussion about the differences between haze and fog and the gray smog weather warning[J]. Meteorological Monthly, 2005, 31(4): 3-7. DOI:10.7519/j.issn.1000-0526.2005.04.001) |
[2] | 廖碧婷, 吴兑, 陈静, 等. 灰霾天气变化特征及垂直交换系数的预报应用[J]. 热带气象学报, 2012, 28(3): 417-424. (LIAO B T, WU D, CHEN J, et al. The forecast application of grey haze weather variations and vertical exchange coefficients[J]. Journal of Tropical Meteorology, 2012, 28(3): 417-424.) |
[3] | ZHOU C H, GONG S, ZHANG X Y, et al. Towards the improvements of simulating the chemical and optical properties of Chinese aerosols using an online coupled model:CUACE/Aero[J]. Tellus Series B:Chemical & Physical Meteorology, 2012, 64(1): 91-102. |
[4] | 邓雪娇, 邓涛, 吴兑, 等. 珠江三角洲空气质量与能见度数值预报模式系统[J]. 广东气象, 2010, 32(4): 18-22. (DENG X J, DENG T, WU D, et al. The numerical forecast model system of air quality and visibility over the Pearl River Delta region[J]. Guangdong Meteorological, 2010, 32(4): 18-22.) |
[5] | WANG T J, JIANG F, DENG J J, et al. Urban air quality and regional haze weather forecast for Yangtze River Delta region[J]. Atmospheric Environment, 2012, 58(15): 70-83. |
[6] | 赵秀娟, 徐敬, 张自银, 等. 北京区域环境气象数值预报系统及PM2.5预报检验[J]. 应用气象学报, 2016, 27(2): 160-172. (ZHAO X J, XU J, ZHANG Z Y, et al. Beijing regional environmental meteorological numerical forecasting system and PM2. 5 forecast test[J]. Journal of Applied Meteorology, 2016, 27(2): 160-172. DOI:10.11898/1001-7313.20160204) |
[7] | 邹乐强. 最小二乘法原理及其简单应用[J]. 科技信息, 2010, 2(23): 282-283. (ZOU L Q. The principle of the least square method and its simple application[J]. Technology Information, 2010, 2(23): 282-283. DOI:10.3969/j.issn.1673-1328.2010.23.278) |
[8] | 高辉. 几类常用非线性回归分析中最优模型的构建与SAS智能化实现[D]. 北京: 中国人民解放军军事医学科学院, 2012. (GAO H. The construction of optimal model in several types of nonlinear regression analysis and the intelligent realization of SAS[D]. Beijing:The PLA Academy of Military Medicine, 2012.) http://cdmd.cnki.com.cn/Article/CDMD-90106-1012421845.htm |
[9] | 原继东, 王志海. 时间序列的表示与分类算法综述[J]. 计算机科学, 2015, 43(3): 1-7. (YUAN J D, WANG Z H. The representation and classification algorithm of the time series[J]. Computer Science, 2015, 43(3): 1-7. DOI:10.11896/j.issn.1002-137X.2015.03.001) |
[10] | 周俊, 张鹏, 刘成良. 基于时间序列分析的卡尔曼滤波组合导航算法[J]. 农业工程学报, 2010, 26(12): 254-258. (ZHOU J, ZHANG P, LIU C L. Based on time series analysis of integrated navigation Kalman filtering algorithm[J]. Journal of Agricultural Engineering, 2010, 26(12): 254-258. DOI:10.3969/j.issn.1002-6819.2010.12.043) |
[11] | 田甜, 王秀玲, 吕芳. 基于Kalman和ARIMA组合模型的路段行程时间预测[J]. 信息技术, 2016(3): 148-150, 155. (TIAN T, WANG X L, LYU F. Urban travel time prediction based on Kalman and ARIMA model[J]. Information Technology, 2016(3): 148-150, 155.) |
[12] | 李杰, 张文栋, 刘俊. 基于时间序列分析的Kalman滤波方法在MEMS陀螺仪随机漂移误差补偿中的应用研究[J]. 传感技术学报, 2006, 19(5): 2215-2219. (LI J, ZHANG W D, LIU J. Kalman filtering method based on time series analysis in the application of MEMS gyroscope random drift error compensation study[J]. Journal of Sensing Technology, 2006, 19(5): 2215-2219.) |
[13] | BOX G E P, JENKINS G M. Time Series Analysis:Forecasting and Control[M]. San Francisco: Holden Day, 1976: 303. |
[14] | 周国强, 韩晟. 基于ARIMA模型的中国医院药品费用增长趋势分析[J]. 药品评价, 2014(14): 15-19. (ZHOU G Q, HAN S. Predictive analysis of China's hospital drug expenditure based on ARIMA model[J]. Drug Evaluation, 2014(14): 15-19. DOI:10.3969/j.issn.1672-2809.2014.14.002) |