近年来,随着计算机能力的迅速提高、气象观测技术的快速发展和数值预报技术的不断创新,数值模式所能考虑的大气物理过程越来越细,模式的水平和垂直分辨率越来越高,相应的数值预报产品的种类和预报准确率也在不断增加和提高.如何开发和应用这些产品,为首都的经济建设和社会发展提供更准确的定时、定点、定量预报,是我们面临的一项新课题.目前对数值预报产品进行释用的方法主要有:统计学释用(MOS 、PP 和卡尔曼滤波等)、天气学释用、动力学释用和人工智能释用等[1] .应用一定的释用方法对数值天气预报模式直接输出产品进行处理,得到本地区天气要素的预报或某种特殊专业预报是目前国际上通行的做法.国内的数值天气预报产品释用技术,与发达国家相比,还有一定的差距,尤其是高分辨率中尺度数值预报产品的释用更为欠缺.
北京地区中尺度数值业务系统的建立[2] 至今已4 年有余.在中尺度数值预报业务和应用方面积累了较丰富的经验,并存储了大量资料,为开展各种释用方法研究提供了较好的基础.利用统计或动力释用方法,对中尺度数值业务预报产品进行后处理,进而得到本地区天气要素预报和空气污染预报,对城市气象预报服务具有重大意义.这是本文要讨论的主要内容.此外,还讨论采用卡尔曼滤波方法制作北京地区9 个自动观测站每隔3 h的温度、风和每天最高、最低气温的预报试验;介绍了如何将中尺度模式产品应用到一个简单的空气污染预报动力模式中,进行北京市区空气污染浓度预报及其预报效果.
1 卡尔曼滤波方法在高分辨率数值预报产品中的应用卡尔曼滤波方法是一种统计估算算法,通过处理带有误差的实际测量数据而得到所需物理参数的最佳估算值.它的主要优点在于能够根据前一时刻的预报误差大小及其它统计量的变化来调整预报方程的系数.它不仅利用了样本所提供的信息,同时也吸收了前一时刻预报方程的反馈信息,从而有利于提高预报精度[3] .利用卡尔曼滤波方法可对数值预报产品进行统计释用,主要适用于制作连续性的天气要素温度、湿度和风速等的预报.该方法有两大特点:一是所需历史资料少,便于建立方程,并且能够适应不断更新的数值预报模式;二是所建立方程的通用性好,使用期限长,便于实际业务应用.
1.1 卡尔曼滤波基本原理和递推公式[4]统计预报方法的最通常的一种做法是在预报量和预报因子间建立回归方程.根据这些统计关系(即回归方程式)可做出相应气象要素或天气现象的预报.
回归方程的一般形式为:
![]() |
(1) |
一般而言,要得到稳定可靠的统计关系,要求样本数量至少为数百个.而卡尔曼滤波采用一种新的思路:只用少量样本(少于100 个)建立回归方程,然后根据前一时刻的预报误差来订正回归(预报)方程中的回归系数.
卡尔曼滤波方法要求滤波对象是离散时间线性动态系统,假定某些气象预报对象是具有这种特征的动态系统,可用下列两组方程来描述:
![]() |
(2) |
![]() |
(3) |
式(2)对应于回归方程(1),在卡尔曼滤波中称为量测方程,其中yt 为预报变量在t时刻的实际观测值,x =(1 ,x 1 ,x 2 ,… ,xm)′(这里上面加一撇表示转置,下同)为预报因子向量,bt =(b0t ,b1t ,b2t ,…,bmt)′为t 时刻的回归系数向量,et 为量测(预报)误差,是一个随机量,其方差为υ.式(3)在卡尔曼滤波中称为状态方程.将回归系数向量bt 作为状态向量,它是变化的,用该状态方程来描述其变化.φt 为转移矩阵,在数值预报产品的统计释用的问题中其值难以确定,一般简单地将其取为单位矩阵.尽管这是一个缺陷,但目前没有更好的办法.εt 为动态噪声向量,与随机量et 相互独立.两者为均值为零的白噪声.εt 的误差方差矩阵为W .
由式(2)、式(3)和上述关于et 、εt 的假定,运用广义最小二乘法,可得到(推导过程繁琐,这里省略)下面一组卡尔曼滤波技术应用于数值预报产品释用时的递推公式:
(i)在t 时刻作对t +1 时刻的回归系数估计,进而给出对预报量在t +1 时刻的估计值(预报值):
![]() |
(4) |
(ii)在t 时刻给出对t +1 时刻的回归系数估计误差矩阵的估计值Pt +1/ t ,进而给出增益向量Kt +1 :
![]() |
(5) |
(iii)得到t +1 时刻预报量的实测值yt +1后,对该时刻回归系数的估计值进行订正:
![]() |
(6) |
并对回归系数估计误差矩阵的预估值Pt +1/ t进行订正:
![]() |
(7) |
(iv)将值t +1 赋给t ,即可进行下一个递推段的计算.
为了启动上述递推过程,需要给定回归系数向量及其估计误差矩阵的初始值
选取9 个有代表性的预报站点,分别是天安门、东直门、西直门、永定门、丽泽桥、怀柔、门头沟、房山、霞云岭.预报对象为各站点的当日17 :00 到次日20 :00 ,时间间隔为3h 的风和气温以及日最低和最高气温.
考虑到近地面气温和风主要是受低层大气影响,故选取了500 hPa 以下和地面的14个数值天气预报产品(见表 1)作为候选预报因子.将模式格点预报值插值到预报站点上,用经验和最优回归方法进行相关性分析统计,挑选出最佳因子组合.其中,最低气温和最高气温的预报因子分别用05 :00和08 :00 及14 :00 和17 :00 的气温平均值代替,而修正回归系数时使用的观测值是实际观测到的最低和最高气温.
![]() |
表 1 14 个数值预报产品候选因子 |
最后确定的预报因子为:用于地面气温预报的因子T 4 、T2 、V 8 和RH13 ;用于地面风U 分量预报的因子U5 、T 2 、T4 和RH12 ;地面风的V 分量预报用V 8 、T1 、T2 、RH12 .
利用2000 年3 、4 和5 月3 个月的中尺度数值预报产品(预报时效为36 h)和北京地区自动站观测资料建立初始回归方程,进而得到
图 1 为由卡尔曼滤波方法预报的天安门站2000 年10 月份最高、最低气温与实测值的逐日对比.从气温的变化起伏可以看出,预报的趋势与实况是基本吻合的,且10 月季节的变化特征及气温的转折也可以反映出来,这说明卡尔曼滤波方法对数值预报产品的解释应用是有效的,具有较好实用价值.图 2 为天安门站2000 年10 月份最低气温三种预报模型的预报结果与实况的对比.图中4 因子指的是使用4 个预报因子的本方法确定的方程,单因子是指仅用2 m 的气温作为预报因子的方程计算所得结果,2 m 温度就是指MM5 直接输出的2 m 温度.从中可以看出,本方法所计算的预报结果好于模式直接输出的2 m 温度,可见通过使用卡尔曼滤波方法对MM5 的预报产品进行解释应用,地面温度的预报有明显的提高.
![]() |
|
图 1. 2000 年10 月份天安门最低和最高气温36 h 预报与实况对比
(-最高气温实况 |
![]() |
|
图 2. 三种不同预报模型所做的2000 年10 月份天安门逐日最高气温36 h 预报与实况的对比 |
表 2 给出了2000 年10 月所选的9 个站12 个预报时次气温预报的月平均绝对误差.可以看出多数时次该误差小于2 .5 ℃,基本在预报误差允许的范围之内.因此预报结果对预报员是有参考价值的.
![]() |
表 2 2000 年10 月9 个站12 个预报时次气温36 h预报的月绝对平均误差 |
图 3 和图 4 分别是2000 年6~10 月9 个站卡尔曼滤波方法与会商室预报员的最低气温、最高气温预报准确率的对比.会商室的预报是指每天15 :30 预报会商后对外发布的当天夜间最低气温和第二天白天的最高气温预报准确率.我们把预报与实况绝对误差小于2 .0 ℃视为准确,准确率为报对的次数除以预报的总次数.总体来看,该系统的预报准确率接近预报员水平,具有很好的参考价值.
![]() |
|
图 3. 2000 年6~10 月9 个站的卡尔曼滤波与预报员的最低温度预报准确率的比较 (1 :天安门,2 :东直门,3 :西直门,4 :永定门,5 :丽泽桥,6 :怀柔,7 :门头沟,8 :房山,9 :霞云岭,黑色为会商室预报员预报值) |
![]() |
|
图 4. 2000 年6~10 月9 个站的卡尔曼滤波与预报员的最高温度预报准确率的比较(说明同图 3) |
表 3 和表 4 分别是2000 年9 个站6~10 月17 :00 卡尔曼滤波方法的风向预报准确率和2000 年10 月份9 个站10 个预报时次风速预报的月绝对平均误差分析.由于北京地方性风为南风、北风,所以在分析风向的预报准确率时按东西南北4 个方向来判断,由于风向的准确率较高,所以风速的月绝对平均误差多数小于1 .6 m/s 也是可信的.因此,如果按级别来预报风速,该系统是有指导意义的.
![]() |
表 3 2000 年9 个站6~10 月逐日17 :00 风向预报月平均准确率 |
![]() |
表 4 2000 年10 月份9 个站10 个预报时次风速预报的月绝对平均误差 |
2 运用中尺度模式预报产品制作空气污染预报
考虑全面的城市空气污染数值预报模式包括数百个化学反应方程并需要掌握详细的污染源和汇的资料,用于业务比较困难.本文采用的方法是在我局中尺度数值天气预报模式的输出产品的基础上,利用引进的CAPPS (City Air Pollution Prediction System)系统[5] ,将其中的大气污染物平流扩散的箱格模式连接到每日业务运行的高分辨率中尺度数值天气预报业务系统,进行每日的污染浓度预报.
2.1 一个简单的污染预报动力模式CAPPS 模式假设在有限的时间积分步长内,大气对污染物清除的总能力和污染源强的变化可以忽略不计,采用有限体积积分法求解大气平流扩散方程[5] ,得到了污染浓度和污染源强之间的明确函数关系:
![]() |
(15) |
其中c 为污染物浓度;c0 为初始污染物浓度;Q 为污染源强;Vc 为大气对污染物清除的总能力,包括污染物的平流和湍流扩散、干沉降和湿沉降;τ为箱格的体积;δT 为时间步长.在已知污染物实际监测浓度的情况下,就可以反算出污染物排放的虚拟源强.因此CAPPS 模式不需要污染源强资料就可以预报出城市空气污染物的浓度.CAPPS 模式每3 h 输出一次SO2 、NO2 、RSP 、CO 的污染指数和污染等级,最后给出24 h 平均的SO2 、NO2 、RSP 、CO 的污染指数和污染等级.
2.2 预报流程及效果在高分辨率中尺度数值天气预报业务系统基础上,运行污染预报动力模式,每日进行污染浓度预报.其业务流程如图 5 所示.
![]() |
|
图 5. 污染预报流程图 |
本文把1999 年12 月13 日~2000 年1 月13 日、2000 年6 月11 日~8 月10 日和9月22 日到11 月25 日三个时段分别代表北京市冬、夏和秋季,利用北京市环境监测中心提供的浓度资料作为污染初始场,对SO2 、NO2 、RSP 、CO 进行了冬季、夏季和秋季24 h 日均值试预报.预报准确率定义为:
![]() |
相关系数为预报污染物指数与监测得到的污染物指数之间的相关系数.为节省篇幅文中只给出了中秋到初冬的预报结果与实测值的比较(图 6).表 5 为2000 年夏秋冬三季SO2 、NO2 、RSP 、CO 日均值浓度等级预报准确率的统计结果,表 6 为CAPPS 污染数值模式中用MM5 和MM4 气象背景场预报结果.综合分析图 6 ,表 5 ,表 6 得出,经修改使用中尺度数值模式MM5 作为污染数值预报的气象背景场,SO2 、NO2 、CO 的预报准确率在夏秋冬三季中均可达到65 %以上.而原有的以MM4 中尺度模式为基础CAPPS 系统,预报的准确率较低(见表 6).RSP(可吸入颗粒物)秋冬季预报准确率在65 %以上,夏季较低,在50 %左右(图略).由于颗粒物污染形成的原因较为复杂,夏季太阳辐射强、气温高、空气相对湿度大,空气中污染物容易发生气-粒转化过程[6~7] ,而基于大气平流扩散方程的污染数值模式没有考虑大气的光化学反应,因而夏季可吸入颗粒物的预报准确率较低.
![]() |
|
图 6. 2000 年中秋到初冬的空气污染指数的24 h 预报与实况的比较(实线:实况,直方条:预报) |
![]() |
表 5 冬夏秋三季SO2 、NO2 、RSP 、CO 4 种污染物日均值浓度等级预报准确率 |
![]() |
表 6 CAPPS污染数值模式中二种不同气象背景场(MM5、MM4)预报结果比较 |
3 小结
(1) 利用卡尔曼滤波方法制作北京地区9 个地面自动站处的每3 h 的风向、风速、温度和日最高/最低温度预报.从预报效果比较来看,温度和风的预报效果高于中尺度数值天气预报模式的直接输出结果.气温的变化趋势和转折都能较好地反映和预报出来.9个站12 个预报时次温度的月平均绝对误差在大多数情况下小于2 .5 ℃,基本在预报误差允许的范围之内,对预报员有较好的参考价值;风向的预报准确率很高,风速的月平均绝对误差多数小于1 .6 m/s ,按级别预报风速,该系统是有指导意义的.
(2) 使用中尺度数值模式(MM5)作为一个简单的污染动力数值预报模式的气象背景场,与原来的以MM4 中尺度模式为基础CAPPS 系统相比预报的准确率有一定的提高.SO2 、NO2 、CO 的预报准确率在夏、秋、冬三季中均可达到65 %以上,RSP 秋冬季预报准确率在65 %以上,夏季较低(为50 %),这与污染模式中没有考虑夏季的光化学反应有关.
[1] | 朱盛明, 曲学实. 数值预报产品统计解释技术的进展. 北京: 气象出版社, 1988. |
[2] | Yingchun Wang, Jianjie Wang, Bo Cui, et al. The Beijing area mesoscale NWP system and its application. Acta Meteor. Sinica, 2000, 14: 233–246. |
[3] | 陆如华, 徐传玉, 张玲, 等. 卡尔曼滤波的初值计算方法及其应用. 应用气象学报, 1997, 8, (1): 34–43. |
[4] | 章国材, 邓北胜, 于玉斌, 等. 卫星气象数据广播接受系统培训教材. 北京: 气象出版社, 2001: 191-198. |
[5] | 徐大海, 朱蓉. 大气平流扩散的箱格预报模型与污染潜势指数预报. 应用气象学报, 2000, 11, (1): 1–12. |
[6] | 王明星. 大气化学(第二版). 北京: 气象出版社, 1999: 186-190. |
[7] | 德利克*埃尔森著. 烟雾警报. 北京: 科学出版社, 1999: 50-53. |