蒙特卡洛与滞后时间集合预报的比较

蒙特卡洛与滞后时间集合预报的比较王丽婧南京信息工程大学大气科学学院，南京210044摘要：本文首先介绍了集合预报中生成初始扰动的两种方法—蒙特卡洛法和滞后时间法，然后利用Lorenz1963年模式，着重研究蒙特卡洛法和滞后时间法的比较，分别在不同的观测误差大小，不同扰动个数的情况下进行了数值试验研究。数值试验结果表明，在观测误差较小或是在不同扰动个数情况下，在预报初期控制预报和蒙特卡洛法的预报效果要比滞后时间法的预报效果好,而在预报中期滞后时间的预报效果要比蒙特卡洛和控制预报的效果好，预报后期滞后时间和蒙特卡洛预报的效果一样；观测误差较大的时候，滞后时间的预报效果在预报前期和预报中期都比控制预报和蒙特卡洛的效果好，后期两种方法的预报效果一致。关键词：蒙特卡洛；滞后时间；观测误差；扰动个数1引言集合预报在国际上被认为是最具发展空间的，用来解决确定性数值预报“不确定性”问题的新一代随机动力预报理论和方法。集合预报的概念是由Lorenz在1965年首先提出，就是针对初值的不确定性，尽可能给出能最好的反映初始时刻误差分布的“误差集”（初始扰动），将不同的初始扰动叠加于预报模式的初始场然后分别制作预报，这样，扰动预报与控制预报的全体就称为集合预报[1]。从集合预报的概念我们可以看出它包含有三项技术问题，分别是初始扰动的生成、数值模式的运用以及预报集合中有用信息的提取，其中主要的一个难点就是初始扰动的生成，因为扰动初值的质量直接影响着模式的预报质量[2]。目前关于生成初始扰动的方法主要有4种：蒙特卡洛法、增长模繁殖法、时间滞后法和奇异向量法。经典的蒙特卡洛（MonteCarlo）法是Leith[3]在1974年通过总结Lorenz和Epstein的集合预报观点，然后在他们的基础上发明的，它是在初始观测场上叠加随机扰动，产生多个初始状态的集合，然后通过对集合成员积分得到预报；增长模繁殖(BreedingofGrowingModels)法是以NCEP为代表，在已有的数值预报模式的基础上，通过模式的繁殖循环求取最快增长模，把最快的增长模作为集合预报的初始扰动；滞后时间（LaggedAverageForecasting）法[4]是1983年Hoffman和Kalnay提出的，它是用相距6小时的数据循环产生相继的分析场作为集合预报的系列初始场，分别作预报然后把相同时刻的求平均；奇异向量（SingularVector）法是在一定原则下选取切线性模式中若干个奇异向量，其中最大的奇异值对应相应的奇异向量，最大的奇异值就代表增长最快的扰动，再按照一定的要求对这些特征值进行线性组合得到初始扰动[5]。当今世界上最为有效的集合预报技术分别为,欧洲中心的基于奇异向量分解的集合预报,以及美国NOAA的基于繁殖模培育法的集合预报技术，但作为最简洁的使用方法还是蒙特卡洛和滞后时间法，国内业务上最常使用的动力延伸预报方法也是蒙特卡洛和滞后时间集合预报，但这两种方法的比较还没有系统性的研究。本文的主要目的就是熟悉集合预报的基本概念和基本方法，主要熟悉蒙特卡洛法集合预报和滞后时间集合预报；通过设计试验研究给出两种集合预报方法的定性比较结论，为今后的实际应用和理论研究打下良好的基础。具体做法是利用Lorenz’63[6]年模式利用吸引子提取变量作为初始场，并利用计算机产生随机数来模拟观测误差，在不考虑模式误差情况下，分别在不同观测误差大小和不同扰动个数下用两种方法分别作预报，计算预报均方根，通过预报均方根对两种方法进行比较。12蒙特卡洛法和时间滞后法介绍经典的蒙特卡洛法认为真实大气中的初始误差是随机分布的。具体做法是：在初始时刻的观测分析场上叠加和扣除不同的随机扰动，从这些具有不同随机扰动的初始场出发分别作预报，然后把这些不同的预报结果求统计平均。集合预报成员的选取和成员个数的确定直接关系到集合预报结果的好坏。按Leith的理论，8个成员足以通过集合预报达到提高预报质量的目的，另外考虑到本文只限理论研究两种方法的比较，所以在试验中成员个数取的较大。滞后时间法的基本思路是取相距6h或12h（取样时间可以自定）的数据同化循环相继产生的分析场作为集合预报的一组初始场，分别对模式积分，最后取相同时刻预报结果的平均值[5]。其具体做法是在过去（-mp）时刻上叠加以随机扰动，预报到0时刻，再从0时刻预报，然后将所有预报中时间相同的预报求统计平均。3数值模型及数值试验设计Lorenz’1963年系统[6]如下：dx/dt=-ax+aydy/dt=rx-y-xz（1）dz/dt=-bz+xy.a,r,b分别为普朗特数、瑞利数、外形比。其中当a=10，b=2.67，r=24.74时表示的是精确的数值模式。当参数r取不同的值时，会对lorenz系统长期积分产生影响，本文重点在研究蒙特卡洛法和滞后时间法的比较，所以在本次的试验设计中不考虑模式参数误差，将参数都取为固定值a=10，b=2.67，r=24.74。取定X，Y，Z，时间积分方案采用简单的一阶欧拉算法，如下：)x(hfxxn1n)y(hfyyn1n（2）)z(hfzzn1nnnnzyx,,分别是方程积分到第n步的值，h是积分步长，本次试验中取为0.01。另外，在本次试验中，我们要考虑不同的观测误差下两种方法的效果对比，而观测误差又具有一定的随机性，所以我们利用计算机产生随机数来模拟观测误差：nnaindhxxnnaindhyy（3）nnaindhzz其中nnnzyx,,表示无观测误差值，nnnzyx,,表示有观测误差值，ain是计算机产生的随机数，它的样本平均值为0,均方差为1。dh为观测误差参数，通过改变dh的大小可以控制观测误差的大小，特别2当dh=0，即观测误差为零时，nnnzyx,,与nnnzyx,,对应相等。接下来进行试验设计：取初始值x、y、z，对初值运行klost步进行适应期计算（klost=5000）得出相应值。接下来继续运行mp步，得到mp个无观测误差的大气状态值nnnzyx,,，接着通过（3）式给nnnzyx,,引入观测误差，得到mp个有观测误差的大气状态值nnnzyx,,。然后用第mp个无观测误差值nnnzyx,,进行理想试验，取预报步数nstep=3000，所得到的数据作为精确值；方案一：控制预报。用第mp个有误差观测值nnnzyx,,作为初值，由式（2）预报3000步，作为控制实验，它相对于理想试验得到的精确值可以求出有观测误差情况下各时刻的均方根误差（rootmeansquarederror,RMSE），就得到控制预报的RMSE曲线。方案二：蒙特卡洛预报。用第mp个有误差观测值nnnzyx,,作为初值，模式预报3000步，得到初始的有误差的大气状态值。然后利用（3）式，给有误差初始场引入mp个随机扰动，最终我们得到（mp+1）个预报成员，将每一个预报结果累加求平均，就得到蒙特卡洛法预报结果，再将其与理想试验得到的精确值求有观测误差下蒙特卡洛预报的各时刻均方根误差，得到蒙特卡洛预报的RMSE曲线。方案三：滞后时间预报。取过去–mp时刻的有误差观测值nnnzyx,,作为初值，利用（2）式对过去mp时刻的观测值积分到0时刻，然后再以积分后的值作为初始值，模式积分3000步。这样循环mp次，得到（mp+1）个预报值，对这（mp+1）个预报值求统计平均，再将其与理想试验得到的精确值求有观测误差下时间滞后法预报的各时刻均方根误差，得到时间滞后预报的RMSE曲线。在以上三个方案中分别取不同的观测误差大小（改变dh）和不同扰动个数作预报，将所得到的均方根误差值画图，就可以得到不同情况下两种方法效果的比较。4不同观测误差大小下两种方法效果的比较为了研究不同的观测误差大小下，蒙特卡洛法和滞后时间法的比较，我们取定mp=10（mp为扰动个数），观测误差参数dh定为变量，取值范围为0.001至6，分别进行控制实验、蒙特卡洛试验和滞后时间试验，通过计算后得到的均方根误差（RMSE）如图1所示：30.001000.002000.003000.000.004.008.0012.00dh=0.001mp=10r=24.74nstep=3000预报步数均方根误差0.001000.002000.003000.000.004.008.0012.00dh=0.01mp=10r=24.74nstep=3000预报步数均方根误差0.001000.002000.003000.000.004.008.0012.00dh=0.1mp=10r=24.74nstep=3000预报步数均方根误差0.001000.002000.003000.000.004.008.0012.00dh=1mp=10r=24.74nstep=3000预报步数均方根误差0.001000.002000.003000.000.004.008.0012.00dh=3mp=10r=24.74nstep=3000预报步数均方根误差0.001000.002000.003000.000.004.008.0012.00dh=6mp=10r=24.74nstep=3000预报步数均方根误差4图1不同大小观测误差下控制预报、MonteCarlo预报、LAF预报RMSE曲线对比图控制实验的RMSEMonteCarlo试验的RMSELAF试验的RMSE分析图1可以得到，在3000步的预报中，两种方法以及控制预报的预报均方根误差都随着预报步数的增加而增大。在观测误差很小（1.0dh）的时候，前500步内三条RMSE值都比较小基本处在6以下，说明预报刚开始的时候的预报效果还是比较好的。到大约500步以后，控制预报的RMSE值迅速增大到10左右，而MonteCarlo和LAF的RMSE曲线只增大到了7附近。也就是从这时开始可以看出MonteCarlo和LAF的优势。另外当005.0001.0和dh时，预报500步以前LAF的RMSE值大于MonteCarlo的RMSE值，从500到1000步之间，LAF的RMSE值又低于MonteCarlo的RMSE值。为了方便看出预报开始时的比较，下面取dh=0.001、0.005、0.01和0.05，mp取10，只画出预报前1000步的RMSE曲线。0.00200.00400.00600.00800.001000.000.002.004.006.008.00预报步数均方根误差dh=0.0010.00200.00400.00600.00800.001000.000.002.004.006.008.0010.00预报步数均方根误差dh=0.0050.00200.00400.00600.00800.001000.000.002.004.006.008.0010.00预报步数均方根误差dh=0.010.00200.00400.00600.00800.001000.000.002.004.006.008.0010.00预报步数均方根误差dh=0.05图2不同大小观测误差下控制预报、MonteCarlo预报、LAF预报RMSE曲线对比图（取预报前1000）控制实验的RMSEMonteCarlo试验的RMSE5LAF试验的RMSE由图2可以看出，预报刚开始的时候(0-第500步左右），LAF的RMSE值最大，也就是它的预报效果最不好，预报一段时间后，LAF的RMSE达到值最小，即它的预报效果达到最好，而MonteCarlo和控制预报的RMSE曲线在前700步时几乎是重合的，到700步以后，MonteCarlo的RMSE值逐渐小于控制预报的RMSE值。说明在预报刚起步时，滞后时间的预报效果不好，而蒙特卡洛的优势也不明显，而到预报中期时，滞后时间的预报效果要优于控制预报和蒙特卡洛预报的。另外可以明显看出，随着观测误差的增大，滞后时间的优势在逐渐提前，当dh取0.1时，在大约预报的第100步滞后时间的均方根误差就已经达到最小了。再看图1，随着观测误差的增大，三种RMSE曲线的值都增大，也就是随观测误差的增大，预报的效果越不好。自dh取0.5往后，