2. 中国气象局数值预报中心,北京 100081
2. Center of Numerical Weather Prediction of CMA, Beijing 100081
数值天气预报是当代大气科学技术进步的最伟大的成就之一,已成为当代天气预报最主要的科学途径,并在极端天气事件的预报方面显示了经验预报方法所不具备的独特优势[1]。然而数值预报所采用的模式初值、模式本身均有误差,加上大气是一个非线性系统,具有混沌特性,这使得数值预报具有不确定性,显著地影响着模式大气的可预报性[1-3]。如果数值预报只提供给用户单一的确定性预报,那么它是不完备的,需要将预报的不确定性定量地提供给用户,从单一值的确定性预报向多个预报值的概率预报转变,集合预报是实现这种转变的关键技术。1992年12月美国(NCEP)和欧洲中期天气预报中心(ECMWF)的中期集合预报系统先后投入业务应用,标志着数值天气预报进入了一个新纪元,随后,其他国家如加拿大、英国、日本、澳大利亚等国的集合预报系统也相继投入业务应用。集合预报作为一个国家数值天气预报系统的重要组成部分,被WMO列为数值预报领域的四个发展方向之一,显示出强大的生命力,能够成功处理数值预报的不确定性问题[2]。
我国全球数值集合预报业务自1998年开始发展。1998年6月国家气象中心在国产神威巨型计算机上建立了基于T106L19全球模式的中期集合预报系统,2001年3月实现业务运行,2006年底建成了基于T213L31的全球集合预报系统。初值扰动方法采用增长模繁殖法(Breeding Vector),集合预报成员为14个,提供集合预报概率,集合平均和离散度的数据、图形、MICAPS格式等的集合预报产品。该系统的业务流程如图 1所示,包括观测资料预处理系统、客观分析系统、多初值产生系统、模式预报系统、模式后处理系统和产品生成系统[4]。
目前,数值预报的不确定性不仅仅是由模式初值的不确定性造成的,模式本身的不确定性也是不容忽视的。目前关于集合预报初值扰动方法的研究比较多,集合预报的最新研究结果表明:同时采用初值扰动和模式扰动的集合预报方法能够更好地表现数值预报的不确定性[5-6]。国外一些集合预报系统中都引入了模式扰动, 使得集合预报的离散度增加了,同时,预报技巧也有所提高[7-8]。我国T213集合预报系统采用的是初值扰动方法,还未对模式进行扰动,构建的集合预报系统离散度偏小,也很难表现模式误差所带来的预报不确定性[4]。国外的研究也表明仅采用初值扰动方法构建的集合预报系统普遍存在离散度偏小的问题[7-9]。实际上,天气预报系统对模式误差很敏感,尤其是次网格尺度物理过程参数化所带来的不确定性对系统的可预报性起着重要作用[10-11]。但是到目前为止,还没有一种理论意义完备且具有实践意义的公式能够将物理过程参数化的不确定性引入到集合预报系统中。目前国际上采用的模式扰动主要有三种方法:(1) 在物理过程参数化方案中加入随机扰动项[7, 12];(2) 使用不同的参数化方案进行随机组合来模拟参数化误差[9];(3) 采用多模式方法[13-14]。结果表明扰动模式物理过程可以增加集合预报成员之间的发散度,提高小尺度天气事件的预报概率。对于全球集合预报系统,各大中心主要采用物理过程随机扰动方案。表 1是国际上全球集合预报系统的技术特征现状。英国气象局对统一模式中的夹卷率系数、临界相对湿度、对流有效位能(CAPE)释放时间、临界弗罗德(Froude)数等进行随机扰动以构造集合预报成员[15]。ECMWF则是在每一次积分时间步长后,在非绝热强迫项中加入随机增倍噪音,当所有参数化过程得到的强迫项的值都加上后,净强迫项值乘上0.5~1.5之间的一个随机数,以此反应模式次网格尺度物理过程参数化中存在的不确定性[7]。ECMWF从1998年10月21日将该方案应用于业务运行。
目前我国T213全球集合预报系统只采用初值扰动方案,没有考虑模式扰动方法,在技术方法上滞后于国际主流的集合预报系统。本文针对这个问题,参考ECMWF的模式扰动方法,设计了我国T213全球中期集合预报系统的物理过程随机扰动方案,并进行集合预报试验,分析评估了T213模式对物理过程随机扰动的敏感性以及集合预报效果,为发展我国T213多初值多模式集合预报系统提供科学依据。
1 T213模式扰动试验方案设计 1.1 T213模式简介国家气象中心T213全球中期数值预报模式是从ECMWF引进的全球谱模式,经过移植改造和自行开发与其匹配的资料分析同化方案(SSI)、模式后处理方案、大规模并行机环境下的自动化运行流程及作业监控方案等,形成我国第四代全球中期数值天气预报系统[16]。T213模式采用地形追随-等压面混合坐标,在垂直方向上有31层,模式层顶到达10 hPa,水平分辨率为60 km,空间格点数为640×320,模式积分时间步长为15 min。模式系统包含一整套比较先进的物理过程, 包括长波辐射方案, 短波辐射方案,对可分辨山脉的描述采用平均地形方法,次网格地形拖曳参数化方案, 湍流扩散方案,云方案和陆面过程,采用质量通量方案描述了各种类型的对流[17]。模式还采用高效的计算方法,如半隐式-半拉格朗日方法、精简格点、分布并行算法等现代技术。T213全球中期数值预报模式是我国T213集合预报系统的基础模式。
1.2 模式随机扰动方案设计本文参照ECMWF的方法,重点研究在辐射传输、湍流混合、次网格尺度地形拖曳、湿对流、积云对流过程产生的综合非绝热强迫项中加入随机噪音,随机扰动方案包括随机数产生方法和物理过程随机扰动方案,详细介绍如下。
1.2.1 随机扰动数产生函数为了获得在给定区间[a, b]内均匀分布的一个随机整数,首先利用公式(1) 产生[0, S]内均匀分布的随机整数。
$\left\{ \begin{array}{l} {R_i} = \bmod \left( {5{R_{i - 1}},4M} \right)\\ RN{D_i} = {\rm{INT}}\left( {{R_i}/4} \right) \end{array} \right.$ | (1) |
式中R0≥1的奇数是初值(随机数种子),S=b-a+1, M=2k,k=[log2S]+1,i=1, 2, 3, …。
然后将每个随机数加a, 即实际需要的随机数为Ri=a+RNDi。本次试验中a=-50, b=50。R0是借助物理过程产生的综合温度趋向项得到随机种子数,然后将此随机整数利用式(2) 获得0.5~1.5区间内的随机数PERTi。
$PER{T_i} = 1 + 0.01{R_i}$ | (2) |
扰动值完全是随机确定的。当给定不同的初值R0后,所获得的随机整数序列不同,对同一格点,当随机数发生器的初值R0不同时,扰动误差型态也完全不同,以此可以获得由物理过程参数化方案产生的温度倾向项、比湿倾向项、风速倾向项的随机扰动值。
1.2.2 物理过程随机扰动方案目前,我国的T213集合预报系统中每个成员可看做是从初始条件开始对模式方程进行数值积分计算。
${e_j}\left( t \right)\int_{t = 0}^t {\left\{ {A\left( {{e_j},t} \right) + P\left( {{e_j},t} \right)} \right\}{\rm{d}}t} $ | (3) |
$\frac{{\partial {e_j}}}{{\partial t}} = A\left( {{e_j},t} \right) + P\left( {{e_j},t} \right)$ | (4) |
其中A和P分别表示大尺度和参数化过程,初始条件可表示为
$\begin{array}{l} {e_j}\left( {t = 0} \right) = {e_0}\left( {t = 0} \right) + \delta {e_j}\left( {t = 0} \right)\\ \quad \quad \quad \quad j = 1,2, \cdots ,m \end{array}$ | (5) |
其中e0(t=0) 表示在t=0时的分析场。j=1表示控制预报,m=15表示包括控制预报在内的共计15个集合预报成员,初始扰动δej(t=0) 是由增长模培育法产生的。
本文尝试在以上集合预报的基础上,通过对辐射传输、湍流混合、次网格尺度地形拖曳、湿对流、积云对流过程这些物理过程参数化方案计算结束后的综合温度倾向项、比湿倾向项、风速倾向项进行随机扰动来模拟物理过程产生的模式不确定性。每个集合预报成员可看做从初始条件式(5) 开始的对扰动模式方程的时间积分,如式(6)。
${e_j}\left( t \right) = \int_{t = 0}^t {\left\{ {A\left( {{e_j},t} \right) + P'\left( {{e_j},t} \right)} \right\}{\rm{d}}t} $ | (6) |
对于每一个格点x=(λ, ϕ, σ)(纬度、经度和垂直混合坐标),扰动了的物理过程可表示为:
${{P'}_j}\left( {{e_j},t} \right) \equiv PER{T_i} \times {P_j}\left( {{e_j},t} \right)$ | (7) |
为了分析评估T213全球中期数值预报模式对本文设计的物理过程随机扰动的敏感性以及该物理过程随机扰动方案对集合预报系统产生的影响,本文设计了以下两个试验方案(表 2)。
(1) 多初值扰动试验方案。采用目前业务增长模培育法产生的7对初值,构造集合预报试验1。
(2) 多初值加上物理过程随机扰动方案。在第一个试验方案的基础上,利用前面介绍的物理过程随机扰动方案,对T213全球中期数值预报模式进行扰动,产生多初值+模式物理过程随机扰动的集合预报试验2。
试验时段为2008年7月20日至7月31日,总计12天。积分开始时间为试验时段内每天的12时(UTC),每天积分240小时,共获得12天的集合预报结果,以下对试验结果进行详细分析。
2 T213模式对物理过程随机扰动的敏感性大气是多种天气系统相互作用的结果,物理过程随机扰动对各预报变量产生的影响特征是我们关心的问题。为研究T213模式对物理过程随机扰动的敏感性,选取2008年7月20—31日集合预报试验1中Member0和集合预报试验2中Member0(表 2)的24~240小时积分结果,定义某一积分时效(如24小时)两者的平均绝对离差Svar(x, y, z),以衡量物理过程随机扰动对预报变量的影响。
$\begin{array}{l} {S_{{\mathop{\rm var}} }}\left( {x,y,z} \right) = \frac{1}{n}\sum\limits_{t = 1}^n {\left| {{{{\mathop{\rm var}} }_1}\left( {x,y,z,t} \right)} \right.} - \\ \quad \quad \quad \quad \quad \left. {{{{\mathop{\rm var}} }_2}\left( {x,y,z,t} \right)} \right| \end{array}$ | (8) |
其中var1(x, y, z, t)和var2(x, y, z, t)分别为集合预报试验1中Member0和集合预报试验2中Member0的预报变量,x, y, z, t分别为经向格点、纬向格点、垂直层次、积分日期,n=12为最大积分天数。
2.1 物理过程随机扰动对预报变量水平空间分布的影响图 2是500 hPa位势高度场24小时至240小时预报的12天平均绝对离差随积分时间的演变,有如下特征:
(1) 500 hPa位势高度场12天平均绝对离差随积分时间快速增长,并且在中高纬度地区的平均绝对离差大于热带地区。24小时预报的500 hPa位势高度场12天平均绝对离差在热带地区小于1 gpm,在北半球中高纬度地区小于2 gpm,在南半球中高纬度地区小于5 gpm。至积分240小时,500 hPa位势高度场12天平均绝对离差在热带大部分地区小于5 gpm,在南北半球中高纬度地区大于50 gpm,在个别地区达到90 gpm以上。这表明:由于大气具有非线性动力系统的不稳定性,对物理过程随机扰动是非常敏感的,模式的微小变动,就能对预报结构产生较大影响,也说明对T213模式进行模式扰动是非常必要的,可以反映由模式误差所带来的预报不确定性。
(2) 物理过程随机扰动对南半球500 hPa位势高度场的影响较北半球大。对于各个预报时效,南半球的离差大小和范围总是大于北半球。造成这种现象可能有两点原因:一是T213模式对南半球的预报技巧较对北半球的预报技巧偏低;二是由于南半球海洋面积较陆地面积大,所以表现出比北半球更大的预报不确定性和对物理过程随机扰动的高敏感性。
(3) 积分后期,在北半球,物理过程随机扰动主要影响高纬度至极地区域,对极地区域的影响最大;在南半球,主要影响中高纬度地区。至积分192小时,500 hPa位势高度场12天平均绝对离差在30°S以南地区大于40 gpm,在60°S附近出现3个大值中心,中心值大于90 gpm;在45°N以北地区大于30 gpm,最大值位于80°N附近,达到了62 gpm。至积分216小时,在30°S以南地区大于45 gpm,在60°S附近出现了多个达90 gpm以上的高值中心;在30°N以北,大于45 gpm的区域明显扩大,在北极地区出现了65 gpm以上的高值中心。至积分240小时,在30°S以南地区大于50 gpm,大于90 gpm的区域主要分布在60°S附近,范围显著增大;在北半球极地区域的最大值达到90 gpm以上。
另外,我们对250 hPa风速、850 hPa温度的12天平均绝对离差分布和随积分时间的演变也进行了分析(图略),发现存在与500 hPa高度场类似的分布和演变特征,不同的是分布状况较500 hPa位势高度场略均匀,在此不再赘述。
2.2 物理过程随机扰动对预报变量垂直分布的影响图 3是全球位势高度场12天平均绝对离差经向平均剖面图随积分时间的演变,可以发现:
(1) 随着积分时间增加,平均绝对离差迅速增大。在24小时预报场中,南半球300 hPa高度上存在一个3 gpm的大值中心;北半球的高值中心位于60°N附近的300 hPa高度上,中心值为1.5 gpm;在240小时预报场中,从南极到50°S附近的300 hPa高度上出现三个并列的高值中心,三个中心的强度分别为70、60和70 gpm;北半球的高值中心位于极地地区,强度为80 gpm左右。
(2) 各层次上的扰动结构具有较高的一致性。对于各预报时次,离差的分布结构很相似,主要表现为从低层到高层,在南北半球中高纬地区为大值区,在赤道地区为小值区。
(3) 从低层到高层,在赤道地区,离差随积分时间增长缓慢;在南北半球中高纬度地区增长速度较快。对于赤道地区的24小时预报场,整层的离差值小于1.5 gpm,之后缓慢增长,至积分120小时,整层的离差值小于6 gpm,此后增长并不明显,几乎达到饱和状态。从积分24小时至240小时,南半球的最大离差值从2.5 gpm增长到70 gpm,高值中心位置变化不大,位于60°S附近的300 hPa高度上;北半球的最大离差值从1.5 gpm增长到80 gpm,高值中心位置从60°N附近逐渐移到90°N附近。
(4) 南北半球中高纬度地区的最大离差值从低层向高层逐渐增加,至300 hPa达到最大,在300 hPa以上变小。
由此可以看出,进行物理过程随机扰动后,随积分时间增加,位势高度场离差在热带地区容易达到饱和,在南北半球中高纬度地区增长较快。其中,南半球(60°S,300 hPa)和北半球(60°~90°N,300 hPa)地区对物理过程随机扰动最敏感。
对绝对涡度、经向风速、纬向风速和温度进行了类似分析(图略),发现12天平均绝对离差经向平均剖面图与位势高度场的分布特征相似,这里不再赘述。
图 4是全球垂直速度场的12天平均绝对离差经向平均剖面图随积分时间的演变。与图 3相比,最大的差别就是在赤道地区,从积分24小时至240小时,存在从低层到高层的大值中心。散度、相对湿度与垂直速度的分布特征相似(图略)。由此可见,与上升运动有关的物理量(如垂直速度、散度等)在赤道地区对物理过程随机扰动很敏感;表征大尺度运动特征的物理量(如位势高度、温度等)在赤道地区对物理过程随机扰动不是很敏感,而在南北半球中高纬度地区对物理过程随机扰动很敏感。
从上述分析可以看出,物理过程随机扰动对预报效果产生了影响,随着积分时间的增加,各预报变量平均绝对离差均迅速增大,表明大气是不稳定动力系统,不仅对初值误差敏感,对模式误差也非常敏感,说明发展T213集合预报模式扰动是非常必要的。
3 T213集合预报效果检验与分析前面讨论了T213全球中期数值预报模式对物理过程随机扰动的敏感性,下面将对表 2构造的两个集合预报试验的预报效果进行检验评估,以深入地分析物理过程随机扰动对集合预报的影响。
3.1 形势场检验集合离散度表征集合预报成员与集合平均的总体偏离程度,也就是集合预报成员之间的发散程度。离散度值小,表示可信度高,同时也说明出现未来实况不在集合预报成员预报范畴内的可能性更高,这也是目前世界上集合预报普遍存在的问题;离散度值大,可信度就会偏低。因此,集合离散度要保持在合适的范围内[2]。由此,我们计算了850 hPa温度的集合预报离散度、均方根误差。
图 5是2008年7月20日至7月31日集合预报试验1和集合预报试验2的850 hPa温度离散度、集合平均的均方根误差的12天平均值的时间演变图。由图可见,集合预报试验1与集合预报试验2的离散度差别不大。积分48至144小时,集合预报试验2的均方根误差大于集合预报试验1,从积分168小时开始,集合预报试验2的均方根误差开始小于集合预报试验1的。这说明加入物理过程随机扰动后,积分后期的均方根误差得到改善。500 hPa位势高度、250 hPa风速也有类似的分布特征(图略),这里不再赘述。
以2008年7月20日的暴雨过程为试验个例,对集合预报试验1和集合预报试验2的降水预报结果进行了对比分析。图 6是7月23日20时24小时累计降水量(实况图),可以看出:主降水带位于山东、江苏、河南南部、安徽北部地区,最大降水中心位于山东和江苏两省交界处,最大降水量在100 mm以上。
图 7a和7b分别是集合预报试验1和试验2的72小时降水量集合预报平均。由图可以看出:两者之间存在细微的差别,例如从湖南地区的降水范围来看,集合预报试验2的降水范围更加接近实况,而集合预报试验1在湖南北部地区出现漏报。由图 8集合平均中雨的TS评分来看,集合预报试验2的TS值在24、96、216小时接近或略低于集合预报试验1,在其他预报时刻都明显高于集合预报试验1。
本文尝试设计T213模式物理过程随机扰动方法,分析模式扰动对T213集合预报的影响,得到以下结论:
(1) 模式物理过程随机扰动方案是在每一次积分时间步长后,通过在辐射传输、湍流混合、次网格尺度地形拖曳、湿对流、积云对流这些物理过程计算结束后的综合温度倾向项、比湿倾向项、风速倾向项乘上0.5~1.5区间内均匀分布的随机数,进行随机扰动来模拟物理过程参数化产生的模式不确定性。
(2) T213全球中期数值预报模式对物理过程随机扰动非常敏感,加入物理过程随机扰动后,各物理量的预报情况发生了明显改变,并且这种变化随积分时间增长而迅速扩大。在水平方向上主要表现为中高纬度地区较赤道地区更敏感;在垂直方向上,表征大尺度运动特征的物理量(如位势高度、温度、风速等)在南北半球中高纬度地区的低层到高层都很敏感,尤以300 hPa最为明显,垂直速度、散度等物理量除了在以上地区敏感外,在赤道地区的底层到高层也非常敏感。
(3) 在多初值集合预报中加入物理过程随机扰动后,对850 hPa温度、500 hPa位势高度、250 hPa风速的离散度没有明显影响,对集合平均均方根误差在积分后期略有改善,对降水预报水平有较为明显的提高,表明模式物理过程随机扰动具有较好的业务应用前景。
本文所做的研究工作仅是初步的,在以后的研究中,对随机扰动参数的选取、集合预报效果的评估等,均需进行更多的试验和评估,以便尽早在我国业务集合预报中加入模式物理过程随机扰动方案,提高集合预报系统的预报技巧,缩短我国与国际先进数值中心在集合预报技术方面的差距。
[1] |
矫梅燕. 天气业务的现代化发展[J]. 气象, 2010, 36(7): 1-4. DOI:10.7519/j.issn.1000-0526.2010.07.002 |
[2] |
陈静, 陈德辉, 颜宏. 集合数值预报发展与研究进展[J]. 应用气象学报, 2002, 13(4): 497-506. |
[3] |
杜钧, 陈静. 单一值预报向概率预报转变的基础:谈谈集合预报及其带来的变革[J]. 气象, 2010, 36(11): 1-11. DOI:10.7519/j.issn.1000-0526.2010.11.001 |
[4] |
田华, 邓国, 胡江凯. 全球T213数值集合预报业务系统简介[C]. 2007年中国气象学会年会, 2007.
|
[5] |
David J Stensrud, Bao Jian-Wen. Using initial condition and model physics perturbations in short-range[J]. Mon Wea Rev, 2000, 128: 2077-2107. DOI:10.1175/1520-0493(2000)128<2077:UICAMP>2.0.CO;2 |
[6] |
Joāo Teieira, Carolyn A Reynolds. Stochastic nature of physical parameterizations in ensemble prediction a stochastic[J]. Mon Wea Rev, 2008, 136: 483-496. DOI:10.1175/2007MWR1870.1 |
[7] |
Buizza R, Miller M, Palmer T N. Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[J]. Quart J Roy Meteor Soc, 1999b, 125: 2887-2908. DOI:10.1002/qj.49712556006 |
[8] |
Alberto Arribas.Results of an initial stochastic physics scheme for the Met Office Unified Model[R]. Ex-eter: Forecasting Research Technical Report No.452, 2004.
|
[9] |
Houtekamer P L, Lefaivre L, Derome J. A System simulation approach to ensemble prediction[J]. Mon Wea Rev, 1996, 124: 1225-1242. DOI:10.1175/1520-0493(1996)124<1225:ASSATE>2.0.CO;2 |
[10] |
谭燕, 陈德辉. 河南"75.8"大暴雨的中尺度集合预报试验[J]. 气象, 2008, 34(9): 11-20. |
[11] |
陈静, 薛纪善, 颜宏. 物理过程参数化方案对中尺度暴雨数值模拟影响的研究[J]. 气象学报, 2003, 61(2): 203-218. DOI:10.11676/qxxb2003.019 |
[12] |
Berner J, Shutts G J, Leutbecher M, et al. A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system[J]. J Atmos Sci, 2009, 66(3): 603-626. DOI:10.1175/2008JAS2677.1 |
[13] |
Harrison M S. Analysis and model dependencies in medium-range forecast: two transplant case studies[J]. Q J R Meteorol Soc, 1999, 125: 2487-2515. DOI:10.1002/(ISSN)1477-870X |
[14] |
严明良, 缪启龙, 沈树勤. 基于超级集合思想的数值预报产品变权集成方法探讨[J]. 气象, 2009, 35(6): 19-25. DOI:10.7519/j.issn.1000-0526.2009.06.003 |
[15] |
Vannitsem S, Toth Z. Short-term dynamics of modelerrors[J]. Atmos Sci, 2002, 59: 2594-2604. DOI:10.1175/1520-0469(2002)059<2594:STDOME>2.0.CO;2 |
[16] |
胡江凯. 国家气象中心T213L31数值预报运行监控方案及预报效果评估[J]. 应用气象学报, 2005, 16(2): 249-259. DOI:10.11898/1001-7313.20050231 |
[17] |
陈起英, 姚明明, 王雨. 国家气象中心新一代业务中期预报模式T213L31的主要特点[J]. 气象, 2004, 30(10): 16-21. DOI:10.3969/j.issn.1000-0526.2004.10.004 |