气象学报  2020, Vol. 78 Issue (3): 451-476   PDF    
http://dx.doi.org/10.11676/qxxb2020.030
中国气象学会主办。
0

文章信息

沈学顺, 王建捷, 李泽椿, 陈德辉, 龚建东. 2020.
SHEN Xueshun, WANG Jianjie, LI Zechun, CHEN Dehui, GONG Jiandong. 2020.
中国数值天气预报的自主创新发展
China's independent and innovative development of numerical weather prediction
气象学报, 78(3): 451-476.
Acta Meteorologica Sinica, 78(3): 451-476.
http://dx.doi.org/10.11676/qxxb2020.030

文章历史

2019-09-02 收稿
2019-12-31 改回
中国数值天气预报的自主创新发展
沈学顺1,2,3 , 王建捷1,2 , 李泽椿1 , 陈德辉1,2 , 龚建东1,2     
1. 国家气象中心,北京,100081;
2. 中国气象局数值预报中心,北京,100081;
3. 中国气象科学研究院灾害天气国家重点实验室,北京,100081
摘要: 数值天气预报是天气预报业务和防灾、减灾的核心科技。中国数值天气预报研究和业务应用一直受到高度重视,在理论、方法和数值模式研究方面取得了有广泛国际影响的研究成果。在回顾新中国数值天气预报自主创新研究成果的基础上,重点对GRAPES(Global Regional Assimilation and PrEdiction System)半隐式半拉格朗日格点模式与物理过程的研发和业务应用的状况以及所取得的重要科学进展进行了综述。近年来,通过自主研发建立了中国数值天气预报业务体系—GRAPES体系。首次以自主技术实现了从区域3—10 km到全球25—50 km分辨率的确定性预报和集合预报系统,并在模式动力框架、四维变分同化和卫星资料同化技术等方面有所突破,建立了大气化学数值天气预报、台风数值预报和海浪预报等系统。自主研发的数值天气预报体系的建立是长期坚持既定科学技术方向以及研究和业务紧密结合、经验不断积累的结果,是中国自主发展数值天气预报技术的重要起点。
关键词: 数值天气预报    GRAPES    半隐式半拉格朗日格点模式    物理过程    卫星资料同化    四维变分同化    
China's independent and innovative development of numerical weather prediction
SHEN Xueshun1,2,3 , WANG Jianjie1,2 , LI Zechun1 , CHEN Dehui1,2 , GONG Jiandong1,2     
1. National Meteorological Center,Beijing 100081,China;
2. Numerical Weather Prediction Center of CMA,Beijing 100081,China;
3. State Key Laboratory for Severe Weather,Chinese Academy of Meteorological Sciences,Beijing 100081,China
Abstract: Numerical weather prediction (NWP) is the core technology of weather forecast and disaster prevention and mitigation. China's NWP research and operational applications have been attached great importance by the meteorological society. Great achievements have been made in the theories, methods and numerical model developments in China, which have a wide range of international impacts. The scientific and technological progresses of NWP since 1949 was summarized. The main part of this paper then devotes to presenting the current status and recent progresses of the domestically developed NWP system, i.e., GRAPES (Global Regional Assimilation and PrEdiction System). In the recent 10 years, a new operational NWP system has been established through independent research and development. This system includes both regional and global deterministic and ensemble prediction models with the resolutions ranging from 3 to 10 km for regional and 25 to 50 km for global forecasts. New progresses have been made in dynamic core, 4-dimensional variational assimilation and satellite data assimilation. The GRAPES system also is extended to drive the atmospheric chemistry model and ocean wave model.
Key words: Numerical weather prediction    GRAPES    Semi-implicit semi-Lagrangian grid point model    Physical processes    Satellite data assimilation    Four dimensional assimilation    
1 引 言 1.1 中国数值天气预报自主创新历史的回顾

数值天气预报是基于数学物理方法客观定量计算未来天气演变的科学。数值天气预报从20世纪50年代试验成功以来,经过70 a的发展,已经成为一个跨学科的复杂而严格的系统性工程,使得天气预报从传统的以统计和经验为主的天气图方法转变成为客观定量的科学(曾庆存,2013Benjamin,et al,2019)。数值天气预报从100年前概念的提出到70 a的发展应用历史表明,数值天气预报的进步是建立在多年稳步持续的科学认知和技术进步的积累之上,被认为是物理学科各领域中最具影响的成就之一(Bauer,et al,2015纪立人,2011)。可以说,数值天气预报是气象科技水平的综合体现,是保障国计民生和国防安全的核心科技的重要组成部分。

中国的数值天气预报研究始于1954年,是国际上较早开展数值天气预报的国家之一,20世纪50年代到60年代末是中国在该领域研究非常活跃的时期,并取得了在国际上有广泛影响的成果(顾震潮,1959Blumen,et al,1973)。20世纪50年代的研究多集中于准地转正压模式。顾震潮(1959)做了系统性的总结,并指出数值天气预报中最根本的还是对大气过程本质的了解、物理因子的重要作用以及在模式的建立、选择和解的方法上必须做大量工作。早期的研究是用图解法求解两层准地转方程,并在寒潮爆发的短期预报、500 hPa高度场的48 h预报方面取得了较好的效果(顾震潮等,1957廖洞贤,1956),值得指出的是,陈雄山等(1957)应用位温坐标系的两层模式研究了锋面的预报问题,这应该是国际上较早将位温用作垂直坐标的研究。1960年2月,中央气象局数值预报组与中国科学院地球物理研究所和计算技术研究所合作,用准地转正压模式制作亚欧范围内24和48 h 500 hPa形势预报(中央气象局数值预报组,1965)。中央气象局数值预报组在《我国的数值预报业务工作》一文中做了很好的总结,指出了差分格式的设计、初始风场和地形影响的重要性,并指出斜压模式是发展的方向(中央气象局数值预报组,1965)。沈如桂等(1965)对该正压模式1964—1965年500 hPa 48 h数值预报图在中央气象局的应用情况进行了系统的总结,指出认识数值天气预报误差和能力、充分利用正确的一面对有效应用数值天气预报的重要帮助作用,这些结论至今在应用数值天气预报结果时仍然具有指导意义。

在开展准地转正压模式研究及在天气预报业务中应用的同时,中国科学家也开展了在准地转模式中考虑斜压过程及原始方程模式的研究。朱永禔(19611962)在苏联科学院应用地球物理研究所访问期间(1959—1960年)发展了一个球坐标3层非线性模式(1000、500和300 hPa),研究了地形对24 h预报的影响。1961年,曾庆存首次用原始方程组做出了实际天气预报,并在莫斯科世界气象中心实现了业务应用,比美国1966年开始用原始方程做预报早了5 a(Zeng,1961曾庆存,2013)。特别是在求解原始方程模式时,曾庆存以大气中运动的适应过程理论为基础,提出“半隐式(或称半显式)差分格式”,对影响时间步长和计算稳定性的快波取隐式格式,增加了计算稳定性且大幅度减少了计算量(Zeng,1961曾庆存,1963a1963c)。Robert(1969)Robert等(1985)进一步丰富和发展了半隐式格式,并与半拉格朗日方法结合,形成了半隐式半拉格朗日方法,这是迄今仍为世界业务数值天气预报模式广泛采用的计算格式。针对数值模式中地形带来的计算复杂性问题,曾庆存(1963b)提出了静力扣除(或标准层结扣除)方法,使得在气压梯度力计算中不再出现基本量的计算误差,避免了大量小差问题。这种静力扣除法后来在谱模式以及半拉格朗日平流计算减缓地形影响中发挥了重要作用(Simmons,et al,1991Temperton,et al,2001)。刘瑞芝等(1965)发展了一个北半球的正压原始方程,该模式采用时间空间中央差分和Lilly(1961)提出的跳点网格,空间网格距为600 km,对4个实例进行了模拟研究,并与地转正压模式的结果进行了对比,得到了比较好的结果。

20世纪50年代及60年代初期,以赵九章、顾震潮、叶笃正、谢义炳、朱抱真、黄仕松等为代表的老一辈科学家,关于动力气象、大气环流、青藏高原动力热力作用、大气运动的适应理论和大气频散理论、风压场的平衡关系等方面的研究成果(徐尔灏,1959赵九章,1959叶笃正,1963),直接或间接地指导并推动了中国数值天气预报的发展(纪立人等,2005)。以周秀骥、巢纪平、周晓平、苏从先、赵柏林等为代表的科学家,在云微物理、积云动力学、中小尺度动力学、边界层湍流等方面的研究成果不仅在理论上当时处于国际前列,而且对在数值天气预报中考虑非绝热过程、理解物理机理等方面起到了重要的理论支撑作用(苏从先,19581959周秀骥,1963巢纪平等,1964)。中国学者在这个时期提出了将初值问题的数值天气预报提升为基于近期历史演变资料的创新性思想(顾震潮,1958a1958b)。之后,以丑纪范(1974)为代表对数值天气预报中使用历史资料的问题进一步深化,将数值天气预报这一微分方程的定解问题变成使用多时刻历史观测资料的泛函极值问题—变分问题。顾震潮的思想和丑纪范引入的泛函极值处理方法,正是今天四维变分同化的主体思想。

20世纪70年代是国际数值天气预报真正实现业务应用,并成为逐日天气预报不可或缺的科技支撑的重要时期(Kalnay,2003)。1975年成立的欧洲中期天气预报中心(ECMWF,European Centre for Medium range Weather Forecast)是这期间的标志性事件。在此期间,中国科学家虽然没有条件开展系统性的研究,但在数值天气预报基础理论研究方面做了大量的工作(陶诗言等,2003);而且,中国科学院大气物理研究所数值预报研究组及时跟踪国际进展,将国际上先进的数值天气预报成果介绍到国内(数值预报研究组,1975)。这期间具有代表性的工作是《数值天气预报的数学物理基础(第一卷)》的出版(曾庆存,1979)。在该著作中,曾庆存针对数值天气预报模式发展中的数学物理基础问题开展了系列阐述,将重点关注的大气流体力学的计算问题进一步升华为科学理论,形成了数值天气预报的数学物理基础。

改革开放之后,中国数值天气预报研究和业务应用迎来了新的契机。在20世纪70—80年代,中国科学家有关数值天气预报基础理论和方法的研究主要集中于与计算格式相关的领域,包括计算稳定性的理论问题和差分格式的构造等(曾庆存,1978曾庆存等,1981季仲贞等,1982)。20世纪90年代后,在可预报性、大气数值模式新算法、资料同化新方法等领域也涌现了在国际上有影响或领先的成果。同时,改革开放之后也是中国数值天气预报模式在自主研发方面活跃的时期,并在中国现代化数值天气预报业务建立的过程中发挥了重要作用(陶诗言等,2003)。

设计适当的计算格式使得离散化的模式方程组满足一定的物理属性或约束是发展数值天气预报模式的核心问题之一,如总能量守恒、总质量守恒和总拟能守恒等。这些约束一方面保证模式具有良好的长期积分性能,另一方面也有助于克服计算不稳定问题。Arakawa(19661972)提出的通过构造格点间合理相互作用的平流雅可比差分格式保持平均动能和涡度拟能的守恒,是研究计算格式的经典之作。中国学者研究指出了Arakawa提出的瞬时守恒格式存在的问题(季仲贞,1981曾庆存等,1981)。曾庆存等(1981)提出了总能量守恒隐式平流项的差分格式,在时间离散和任何初值情况下均计算稳定,具有能量守恒、广义能量守恒和平均尺度守恒。尽管隐式完全能量守恒格式计算是稳定的,但实际求解困难,需费巨大的工作量。王斌等(1990)季仲贞等(1991)在此基础上提出了显式完全平方守恒格式。钟青等(19921993)在显式完全平方守恒格式上进一步发展了完全平方和三阶守恒格式。也有学者探索了全隐式差分格式的可能应用(陈长胜等,2007)。这些与计算格式有关的研究显示了中国科学家对数值计算基本方法研究方面的一贯重视。

中国在可预报性研究方面也做出了国际上有影响的研究成果。Mu等(2003)将线性奇异向量(linear singular vector,LSV)推广到非线性领域,提出条件非线性最优扰动(Conditional Nonlinear Optimal Perturbation,CNOP)。线性奇异向量代表了切线性模式中具有最大增长率的一类初始扰动;条件非线性最优扰动代表了满足一定物理约束条件下,在预报时刻具有最大非线性发展的一类初始扰动。Mu等(2010)进一步将条件非线性最优扰动拓展到初始误差和模式参数误差同时存在的情形,将仅与初始扰动有关的条件非线性最优扰动称为CNOP-I;与模式参数扰动有关的条件非线性最优扰动称为CNOP-P。条件非线性最优扰动在各类时间尺度的可预报性研究、目标观测和集合预报中有广泛的应用前景。

现有世界主要数值天气预报业务中心的谱模式和格点模式均采用球面高斯或经纬度网格以及半隐式半拉格朗日时间积分方案。这些模式尤其是全球模式在应用众核高性能计算机实现千米尺度全球数值天气预报时,网格和时、空离散化算法存在可扩展性的计算瓶颈问题,不能满足高效及时、高可扩展性、高分辨率高精度的数值天气预报业务发展需求。因此,世界主要数值天气预报业务中心均将可扩展性作为发展下一代模式考虑的重要问题。研发高精度、守恒、适用百米到数十千米分辨率的数值天气预报模式动力框架和相应的尺度自适应物理过程是发展下一代数值天气预报模式的重点(Mengaldo,et al,2019)。沈学顺及其研究团队研究并发展了基于多矩约束的有限体积离散化算法求解全可压大气方程组的新方法,该方法严格保证数值守恒,与传统的有限体积法相比具有局地自由度高、各类网格适应性好、无全局通信、算法局地性强等特点,具有高的可扩展性(Li,et al,2013Chen,et al,2014)。该方法得到了国际发展同类下一代模式同行专家的认可(Smolarkiewicz,et al,2016)。在尺度自适应物理过程这一前沿研究领域,Huang等(2018)基于Arakawa提出的积云对流尺度自适应概念发展了一个新的表征积云对流参数化尺度自适应的表述方法;Zhang X等(2018)发展了一个三维湍流尺度自适应参数化方案,在百米级分辨率和千米尺度数值模拟中取得了物理上更为合理的结果。

先进的资料同化技术被认为是数值天气预报系统业务应用效果大幅度提高的重要原因之一(Bannister,2017)。中国科学家有关资料同化的系统性研究相对于计算格式等的研究起始较晚,但发展很快。结合数值天气预报业务需求,在20世纪80年代至90年代初,曾经发展了采用最优插值客观分析与非线性正规模初始化相结合的初始化方案(屠伟铭等,1995朱宗申等,1992薛纪善等,1992)。特别是,进入21世纪以来,中国科学家研究发展的资料同化系统已经实现了先进资料同化系统的基本功能(Zhang L,et al,2019)。GRAPES变分同化系统是21世纪中国变分类同化方法研究的一个代表性综合成果(薛纪善等,2008)。

中国科学家也开展了资料同化新方法的研究。Wang B等(2010)提出的降维投影四维变分同化(DRP-4DVar,4DVar based on a Dimension-Reduced Projection)以及Tian等(2018)提出的集合四维变分同化方法(NLS-En4DVar,Non-linear Least Squares enhanced proper orthogonal decomposition based Ensemble 4DVar algorithm)是代表性的研究工作。DRP-4DVar方法的基本思路是利用集合方法替代伴随方法实现代价函数梯度的计算,将模式变量空间和观测变量空间降维投影到样本空间,并结合使用局地化相关函数的主模态对样本空间进行扩展,最终在扩展样本空间完成四维变分分析。DRP-4DVar在计算目标函数梯度中回避了切线性与伴随模式,计算量大幅度减少。NLS-En4DVar实现了4DVar (four-dimensional variational data assimilation)与EnKF (Ensemble Kalman Filter)的优势互补,该方法通过严格的理论推导将集合四维变分同化问题转换为非线性最小二乘最优化问题,实现非线性最小二乘最优化问题的理论体系与集合数据同化方法的融合。这两种四维变分同化新方法,显示出较强的业务应用潜力。

20世纪70年代后期开始,也是中国自主研制数值天气预报模式的热点时期。国家气象局与中国科学院大气物理研究所(以后简称大气所)合作研制的亚欧区域3层原始方程绝热模式(当时称为A模式,1980年开始发布48 h形势预报)、中国科学院大气物理研究所-北京大学-国家气象中心联合研制的北半球5层格点模式(当时称为B模式,1982年业务运行;陶诗言等,2003)和一个有限区域的5层格点模式(当时称为B小模式,1982年业务运行)、国家气象中心发展的LAFS(Limited Area Forecast System;后逐步发展为High-resolution LAFS—HLAFS;该模式是以北京大学张玉玲教授提出的基本模式架构建立起来的)有限区域模式(郭肖容等,1995),代表了20世纪80年代中后期中国业务数值天气预报模式的自主创新成就,并支撑了中国20世纪80年代到90年代初的逐日天气预报和防灾、减灾工作(陈德辉等,2004),其中,HLAFS一直业务运行至2006年才被GRAPES_Meso取代。同时,上海、广州、兰州等地都独立研制了各自的模式,如上海的台风模式(朱永禔等,1987)、广州的热带模式(Xue,et al,1988)、兰州的复杂地形模式(颜宏,1987)等,并且都曾经应用于预报业务。针对中国的复杂地形和暴雨预报问题,宇如聪参考Mesinger(1984)的工作于1986年开始发展一个通过引入η坐标考虑陡峭地形的预报模式(后称为REM,Regional Eta-coordinate Model;宇如聪等,1994)。经过多年对模式物理过程的更新、完善,REM进一步发展为AREM(Advanced Regional Eta-coordinate Model),在武汉暴雨研究所和军队等单位的科研和业务中曾得到应用(宇如聪,2004)。这时期,除了面向业务的数值天气预报模式自主创新研究外,中国科学家自主研制的模式在数值天气预报研究中也发挥了重要作用。郑庆林(19801989)基于完全形式的涡度方程和散度方程,采用球函数为基底,构造了7层半球谱模式,并成功给出了一批数小时的实例预报结果。他的工作大致和国际上谱模式的研究同步,其模式对非线性项的计算采用经典的“相互作用系数法”,而不是现在流行的“变换方法”。中国也比较早地开展了非静力平衡模式的研究,但多限于研究之用,如非静力平衡研究模式(王东海等,1996)、三维非静力平衡研究模式(胡志晋等,1991)等。

需要指出的是,文中重点关注对中国数值天气预报方面自主创新成果的综述,不包括长期积分的大气环流模式和气候模式。有关大气环流模式自主创新发展的成果可以参见曾庆存等(1985;IAP AGCM-I,Institute of Atmospheric Physics Atmospheric General Circulation Model-I)和Wang等(2004;GAMIL,Grid-point Atmospheric General Circulation Model of IAP/LASG)以及沈学顺等(2010)的综述。

1.2 中国现代数值天气预报业务体系的建立和创新发展

如上所述,中国科学家虽然在20世纪50—60年代早期探讨性地开展了数值天气预报业务试验,但限于当时的条件很难形成业务能力。真正形成实际业务能力则始于1980年之后A模式、B模式的业务运行,建成了短期数值天气预报业务系统。中国全球中期数值天气预报业务系统的研发始于20世纪80年代中期,1985年被国家科学技术委员会批准为“七五”国家重点科技攻关项目(李泽椿,2010)。鉴于当时中国数值天气预报系统的研发能力和支撑条件与先进国家存在很大差距,确定了预报模式从先进国家整体引进的技术路线。“七五”期间,李泽椿主持了17个单位参加的近300人的跨部门联合攻关小组,以从ECMWF引进的全球谱模式为核心,在国产银河大型机上建立了T42L9全球预报系统,并于1991年投入业务运行。在20世纪90年代与21世纪初几经升级(T63L16、T106L19、T213L31、TL639L60;之后称为T系列),在中国的气象预报业务与服务中发挥了重要作用,为之后开展自主研发工作奠定了重要基础。20世纪90年代该系统的建立标志着中国具备了从全球中期到区域短期的现代化数值天气预报业务和服务能力。

20世纪末至21世纪初,中国气象局做出了未来数值天气预报业务发展技术路线由引进为主转为自主开发为主的重大决策,开启了中国新一代数值天气预报系统的自主开发工作。2001年开始,在科学技术部“十五”国家重点科技攻关项目“中国气象数值预报系统创新研究”支持下,中国气象局联合多家单位,自主研究并初步建立了新一代多尺度通用资料同化与数值天气预报系统—GRAPES(薛纪善等,2008)。之后,在科学技术部“十一五”“十二五”科技支撑计划项目以及中国气象局GRAPES专项的连续支持下,经过十几年持续不懈的努力,以自主技术研发建成了从区域3—10 km到全球25—50 km分辨率的确定性与集合预报的完整数值天气预报业务技术体系,并锻炼培养出一支数值天气预报业务全链条(即从观测资料前处理和质量控制、资料同化、模式动力物理、模式并行计算到模式系统集成、数值天气预报产品后处理等)研发队伍。通过GRAPES的研发,中国在非静力全可压格点动力框架、四维变分同化、云降水物理、高精度大气模式数值算法、卫星资料偏差订正和卫星资料同化技术等方面取得了创新性的成果,为今后中国数值预报的全面深入发展奠定了重要科研和技术积累以及人才储备。

伴随数值天气预报业务的发展,一些学者对不同阶段的数值天气预报业务研发进展也做了总结和综述(陈德辉等,2004;李泽椿,2010沈学顺等,2010)。文中在回顾中国数值天气预报自主创新发展历史的基础上,重点总结近10年GRAPES的发展成果,包括GRAPES科学技术进展、业务系统研发和应用状况以及未来发展技术路线。

2 GRAPES研发历程和业务体系的形成

GRAPES由区域和全球一体化模式及同化系统构成,其研发酝酿于2000年初。发展非静力格点大气模式以克服高分辨率谱模式在大规模并行计算时的瓶颈问题,是当时国际数值天气预报界的共识。应该说,在设计GRAPES系统时做出的格点模式技术方向选择是科学合理的。同时,发展同化系统的策略则采取先三维再逐步走向四维变分同化的路线,也是符合科学发展规律的。

2.1 研发和业务应用的重要节点

至“十五”科技攻关项目结束(2001—2005年),完成了半隐式半拉格朗日非静力动力框架、适用于中尺度数值天气预报的物理过程、区域等压面三维变分同化及并行计算的研发。基于这些成果,形成了GRAPES中尺度数值天气预报系统(GRAPES_Meso 2.0),并于2006年7月实现业务应用(30 km水平分辨率),代替HLAFS成为中央气象台业务中尺度预报系统。

GRAPES_Meso 2.0是GRAPES的第一个业务应用成果。2007年,经过水平分辨率从30 km到15 km的升级,形成了GRAPES_Meso 2.5版。在业务应用中,发现预报雨带散乱、降水强度偏低等系统性问题。针对此状况,研发了基于分段有理函数的高精度守恒标量平流方案PRM (Piecewise Rational Method)、考虑陡峭地形的单调水平扩散等技术,并依据模式有效分辨率在模式中考虑了有效地形。这些改进大幅度提高了模式对水汽等水物质的模拟精度,显著改进了降水预报效果(沈学顺等,20112013)。高精度守恒标量平流PRM的应用也显著改进了GRAPES区域台风预报模式(GRAPES_TYM)对台风强度的预报效果。这些成果使得GRAPES中尺度数值天气预报系统有了第一次重要升级(2010年2月GRAPES_Meso 2.5版升级到3.0版)。

GRAPES_Meso业务应用中影响预报效果的另外一个关键问题是中尺度资料同化系统。之前,GRAPES_Meso所用三维变分同化是以等压面上的位势高度、风及水汽(比湿或相对湿度)作为基本大气状态变量,与GRAPES_Meso预报模式并不匹配。无论预报模式给同化系统提供背景场或同化系统向模式提供分析结果作为初始场都要经过附加的中间过程,不可避免会引入额外的误差(薛纪善等,2012)。这个问题在研发GRAPES全球中期预报系统的初期同样存在。此外,还缺乏能够为中尺度模式提供云水等云信息初始场的同化功能或者云分析系统。为此,研发了模式空间三维变分同化系统,形成了与GRAPES模式匹配的区域/全球一体化的三维变分同化系统(薛纪善等,2012),这也成为后续GRAPES全球四维变分同化系统顺利研发成功的重要一步。在区域模式空间三维变分同化系统的基础上,Wang R C等(2018)进一步将控制变量更新为纬向和经向风,且实现了温度和地表气压的直接同化,并通过背景误差协方差、平衡约束关系、水平相关模型等的优化与研发以及大尺度约束的引入,使之能够有效应用在千米尺度分辨率的GRAPES_Meso系统中。同时,以ARPS(Advanced Regional Prediction System)中的云分析方案为基础,通过优化改进,研发了有效应用分钟级探空、地面观测、风云二号静止气象卫星云产品及多普勒天气雷达三维组网拼图产品等资料的云分析系统,通过松弛逼近方法在GRAPES_Meso中实现了对云微物理初始信息的应用(朱立娟等,2017)。以上研发的这些核心技术成为继GRAPES_Meso 3.0之后的重要科学技术进步。当然,包括积云对流、云物理及动力物理耦合计算的优化工作也对GRAPES_Meso的不断进步起到了积极作用。

以上工作促进了GRAPES_Meso业务应用效果的不断提升,并使得GRAPES_Meso系统得到多样化应用。基于GRAPES_Meso的快速循环系统(GRAPES_RAFS)和区域台风预报系统(GRAPES_TYM)于2012年实现了业务应用。2014年,GRAPES_Meso分辨率从15 km升级为10 km(黄丽萍等,2017),并在垂直方向进行了加密(由原来的30层增加到50层)。2016年建立了中国东部区域3 km分辨率的GRAPES_Meso准业务系统,该系统在2019年进一步扩展为覆盖全中国的3 km分辨率预报系统,在逐日天气预报业务中发挥了重要作用。

2007年7月开始系统地研发GRAPES全球中期数值天气预报系统(GRAPES_GFS)。在半隐式半拉格朗日非静力动力框架和扩展至全球的等压面三维变分同化系统的基础上,研发了适用于全球预报的动力框架、物理过程以及全球卫星资料的系列同化应用技术。2009年3月完成了GRAPES_GFS的前期试验,并在业务环境中建立了运行试验平台,确定了GRAPES_GFS的准业务版本—GRAPES_GFS 1.0,实现准业务化运行(沈学顺等,2009)。这期间,不仅研发了GRAPES全球预报模式、全球同化和卫星资料应用等核心技术,而且建立了全球多卫星平台资料的业务接收和处理流程,数值天气预报自主创新发展过程也带动了气象业务全链条的进步发展。

但是,不可否认的是,GRAPES研发时采用了传统两时间层半隐式半拉格朗日算法(Temperton,et al,2001Gospodinov,et al,2001),为保持积分稳定,模式中取较大的半隐式权重系数使得时间离散化精度仅为1阶,且计算上游点和方程右端非线性项时需要时间外插,造成模式计算噪音并导致个别情况下计算不稳定;同时,采用了非守恒方程组和半拉格朗日算法,模式的质量守恒问题在算法设计上没有很好考虑。另外一个严重影响全球预报效果的问题是资料同化,如前所述,GRAPES_GFS 1.0的同化系统是一个在标准等压面上的三维变分分析系统,分析变量为位势高度、风和相对湿度,在17层标准等压面和Arakawa-A网格上分析后,通过静力关系导出模式预报变量,再插值到模式格点上作为初值。由于插值误差,等压面3DVar存在精度上的局限,对预报初值的精度带来严重影响。这些影响预报效果的问题涉及到模式动力框架、同化框架和数值算法的根本问题,也是GRAPES在2005年初步研发成功之后走向业务应用所面临的重要挑战。

在GRAPES_GFS 1.0准业务化之后,开展了3个阶段的攻关。第一阶段以完善模式物理过程为代表,包括辐射引入RRTMG(The Rapid Radiative Transfer Model for GCMs)方案提高辐射计算精度、以中国自主研发的CoLM(Common Land Model)陆面模式代替简单的5层热扩散模型(SLAB)、引入次网格地形重力波参数化、对积云对流参数化进行优化,这些工作使得GRAPES全球模式具备了适用于全球预报的一套物理过程。第二阶段,重点攻克模式空间GRAPES三维变分同化系统以避免同化后插值到模式空间带来的误差,并在改进卫星资料同化应用效果的同时增加更多卫星资料的同化应用(薛纪善等,2012)。第三阶段以同化和模式协同攻关为代表,重点解决了动力框架中部分影响计算精度和稳定性的问题、以云-降水为重点的物理过程协调性改进和优化、同化系统中风压平衡问题、背景误差协方差随季节变化问题、同化系统中背景廓线的精确插值和求解、观测资料的质量控制和偏差订正以及更多卫星资料的同化应用等。这些成果不仅在GRAPES_GFS中取得很好的效果,同时也在GRAPES_Meso系统中得到了应用(具体的改进内容参见沈学顺等(2017))。

在上述研究成果的基础上,建立了水平0.25°×0.25°分辨率垂直方向60层的全球预报系统新版本—GRAPES_GFS 2.0,并于2015年年底通过业务化验收(沈学顺等,2015)。到目前为止的业务运行结果表明,GRAPES_GFS 2.0预报形势场统计检验指标虽与ECMWF、NCEP相比尚存在些许差距,但全面超越了国家气象中心引进的业务模式TL639,且雨带预报能力已经接近ECMWF的预报水平,不同阈值下的模式空报问题得到了抑制,中国区域的降水预报形势,特别是中国东南部的主降水区的形势与实况更为接近。

在GRAPES_GFS 2.0基础上,进一步解决了GRAPES-4DVar面向业务应用存在的问题。重点解决了计算效率、观测剖分、多次外循环更新等技术细节问题,并发展了包括垂直扩散、深积云对流、次网格地形阻塞流拖曳、大尺度凝结等线性化物理过程,研发了卫星资料动态偏差订正技术,引入数字滤波弱约束,并对同化框架中的风压平衡、背景误差协方差等进行了进一步优化。有关GRAPES-4DVar的具体科学技术细节可参考Zhang等(2019)。2018年7月带有GRAPES-4DVar的GRAPES_GFS 2.4实现业务运行。GRAPES-4DVar的业务运行标志着中国业务数值天气预报同化技术迈入国际前列,成为国际上少数具有自主研发和业务应用四维变分同化系统的国家级预报中心之一。

集合预报是业务数值天气预报系统的重要组成部分。2010年,中国气象局数值预报中心开始研发GRAPES中尺度集合预报(GRAPES_REPS,GRAPES Regional Ensemble Predication System)和全球集合预报系统(GRAPES_GEPS,GRAPES Global Ensemble Predication System)。对于中尺度集合预报系统,先后研发了考虑初值不确定性的ETKF(Ensemble Transform Kalman Filter)方法,以及考虑模式不确定性的随机物理过程倾向扰动SPPT(Stochastically Perturbed Parameterization Tendencies)和随机动能后向散射SKEB(Stochastic Kinetic-energy Backscatter)方案。GRAPES_REPS于2015年实现业务应用,分辨率为15 km×15 km,2019年GRAPES_REPS进一步升级为水平分辨率10 km×10 km和垂直50层。对于全球集合预报,基于GRAPES全球切线伴随模式,研发了基于奇异向量(SV)构造初值扰动的方法,并研发了随机物理过程倾向扰动(SPPT)和随机物理参数扰动(SPP)。全球集合预报系统GRAPES_GEPS于2018年实现业务应用,水平分辨率为50 km×50 km,垂直方向60层。

伴随GRAPES模式和同化技术的不断深化发展,至2018年实现了从区域3—10 km到全球25—50 km分辨率的确定性与集合预报的完整数值天气预报体系,并建立了观测资料前处理和质量控制、数值天气预报产品后处理、预报检验和产品解释应用、观测和预报数据库、试验分析平台、监控系统全链条支撑系统。回溯数值天气预报自主研发历程,可以看到,与GRAPES体系之前的T系列相比,中国业务数值天气预报系统又向前跨了一大步,实现了自主研发、应用、改进、发展的良性循环。尤为重要的是,在坚持不懈的自主发展过程中,培养了一支稳定、全技术链条的研发和业务应用队伍。

2.2 GRAPES业务体系概述及应用效果 2.2.1 主要科学技术特点

引进的T系列自1991年开始服务中国气象预报和防灾、减灾逾20年,在国家、省、市、县各级预报员中有着深刻的影响。基于自主研发的GRAPES成果构建全新的业务体系满足预报业务的需要并不是一件简单的事。从上一节对GRAPES研发重要节点的叙述中可以看到,研究成果的叠加与反复迭代、在应用中研发业务系统并不断改进需要团队长期的研究和知识经验的积累。沿着既定的技术路线、注重研究成果的积累、在应用实践中解决科学技术难题谋发展,对于自主发展数值天气预报业务体系并扎实掌握气象预报核心科技而言至关重要。

自2005年GRAPES研发成功之后,历经13 a的深化改进,到2018年建成了完整数值天气预报体系。与早期相比,GRAPES的动力框架、模式物理过程、同化框架以及观测资料应用等诸方面有了长足的进步。有关GRAPES的科学技术设计和初期研发成果参见《数值预报系统GRAPES的科学设计与应用》(薛纪善等,2008)。

表1可以看出,模式动力框架在科学技术方案上有了长足的进步,与目前英国气象局、加拿大气象局数值天气预报模式采用的预估修正半隐式半拉格朗日方案类似,与传统半隐式半拉格朗日相比时间离散化精度更高、稳定性更好,突破了传统半隐式半拉格朗日中非线性项和拉格朗日轨迹中间点风速的时间外推带来的计算噪音问题。GRAPES业务体系区域/全球预报系统采用一体化框架的非静力格点半隐式半拉格朗日模式,在当前国际业务数值天气预报单位中也只有英国气象局和加拿大气象局;另外,全球预报采用非静力模式的目前仅有英国(ENDGame/UKMO)、美国(FV3_GFS/NCEP)、德国(ICON/DWD)和中国(GRAPES/CMA)(Mengaldo,et al,2019)。模式物理过程的研发与应用也有了长足的进步,从早期借用WRF的物理过程,到自主研发双参数云物理过程,并在此基础上研发了云宏物理和预报云方案(Ma,et al,2018),同时研发了次网格地形重力波拖曳和小尺度地形湍流拖曳等过程,使得GRAPES模式的物理过程更加完备。GRAPES模式物理过程在应用中得到了不断的改进与发展,成为提高区域/全球预报技巧的重要因素之一。

表 1  GRAPES模式现状与研发初期的对比 Table 1  Comparison of present GRAPES model with the version at early development stage
动力框架 物理过程
GRAPES 2.0 GRAPES业务体系 GRAPES 2.0 GRAPES业务体系
1. 非静力全可压方程组
 (浅大气近似)
2. 等温参考大气
3. 传统两时间层半隐式半拉 格朗日(SISL)算法
4. 动量方程三维矢量离散化
5. 传统 SISL 求解位温预报 方程
6. 准单调正定标量平流
  QMSL
1. 非静力全可压方程组(浅 大气近似)
2. 三维参考大气
3. 预估修正半隐式半拉格朗 日算法
4. 动量方程三维矢量离散
 化+ W damping
5. 垂直方向非插值 SISL 求 解位温预报方程
6. 高精度守恒标量平流
  PRM
1. NCEP-3/5 种水物质的微物
 理方案
2. Betts-Miller 积云对流参数化
3. Dudhia 短波、RRTM 长波辐
 射方案
4. MRF 边界层参数化
5. SLAB 陆面过程
6. 动力物理反馈:线性插值
1. 双参数云物理方案或 WSM-6
 (中尺度模式)
2. 宏云物理与预报云方案(全球模式)
3. NSAS 积云对流参数化或 Meso-SAS
 (区域模式)
4. RRTMG 长短波辐射
5. CoLM (全球模式)或 Noah(中尺度模
 式)陆面过程
6. NMRF 或 YSU (区域模式)边界层参数化
7. 次网格地形重力波拖曳
8. 小尺度地形湍流拖曳
9. 动力物理反馈:三次样条、准单调三次
 样条插值

表2可以看到,初期是与预报模式完全不匹配的等压面三维变分同化,历经10多年发展,开发了基于GRAPES非静力切线和伴随模式的全球四维变分同化、适用于千米尺度高分辨率数值天气预报的区域三维变分同化,与模式进步一样展现出了显著的进步,具体的技术进步在后文进一步阐述。全球GRAPES-4DVar的业务应用是中国自主发展数值天气预报的重要节点,中国成为继ECMWF、日本、法国、英国和加拿大之后业务运行四维变分同化的国家。

表 2  GRAPES变分同化系统现状与研发初期的对比 Table 2  Comparison of present GRAPES data assimilation system with the system at its early stage
GRAPES 2.0 
区域等压面三维变分同化 
业务体系区域  
GRAPES-3DVar  
业务体系全球  
GRAPES-4DVar  
基本公式 增量分析方案 增量分析方案 增量分析方案
网格 水平 A 网格、垂直标准等压面 水平 C 网格,垂直高度地形追随坐标 水平 C 网格,垂直高度地形追随坐标
分析变量 φu、v、RH/q Tpsuv、q/RH/RH* πu、v、q/RH/RH*
背景误差协方差 NMC 方法 NMC 方法 NMC 方法
平衡关系 线性平衡 线性平衡+统计 线性平衡+统计
最优化方法 有限记忆变尺度方法(LBFGS) 有限记忆变尺度方法(LBFGS) Lanczos-CG
初始化 全量数字滤波 增量数字滤波 数字滤波弱约束
1. 非静力切线伴随模式
2. 切线性物理过程(垂直扩散、次网 格地形阻塞流拖曳、深积云对流、 云物理)

观测资料的精细化处理与同化是影响数值天气预报效果的重要因素(Magnusson,et al,2013)。与早期GRAPES_Meso业务系统仅使用常规资料相比,目前的GRAPES系统具备了同化应用更多非常规资料的能力,包括全中国雷达反射率和径向风、地基GPS可降水量、中外静止与极轨卫星多平台多通道辐射率、GPS掩星资料等。各种类型观测资料处理和同化应用能力的提升是自主发展数值天气预报技术的重要成果之一。具体可参考有关文献(Han,et al,2010Liu,et al,2014Li,et al,2014王金成等,2016万晓敏等,2017)。表34分别给出了目前GRAPES全球和区域预报系统中的资料应用情况。

表 3  GRAPES_GFS所用观测资料统计(2019年7月18日03时—27日09时) Table 3  Utilization of observations in GRAPES_GFS based on statistics during 03:00 BT 18—09:00 BT 27 July 2019
观测类型 仪器 平台 变量 占比(%) 分类占比(%)
常规资料 TEMP uvpq 3.4 30.1
SYNOP p 3.1
SHIP p 0.2
BUOY uv 0.1
AIREP uvT 23.3
卫星资料 AMSUA NOAA-15,-18,-19,Metop-A,-B Radiance 9.6 69.9
MWHS-2 FY-3C,-3D Radiance 0.2
GIIRS FY-4A Radiance 0.7
ATMS NPP Radiance 4.1
IASI Metop-A,-B Radiance 39.8
AIRS EOS-2 Radiance 2.3
GNSSRO COSMIC-1,Metop-A,B,FY-3C/D Refractivity 2.1
AMVs-IR/VIS/WVs FY-2E,GOES-13,GOES-15,MTSAT-2,METEOSAT_10 uv 10.4
ASCAT Metop-A,-B uv 0.7
表 4  GRAPES_Meso所用观测资料统计(2018年7月1日06时—7日00时) Table 4  Utilization of observations in GRAPES_Meso based on statistics during 06:00 BT 1—00:00 BT 7 July 2018
观测类型 仪器   平台 变量 占比(%) 分类占比(%)
常规资料 TEMP uvT、RH 4.22 6.24
SYNOP p、RH 0.74
SHIP p、RH 0.04
AIREP uvT 1.24
非常规资料 Doppler radar SA/SB/CB/SC/CD/CC VAD wind
Radial wind
Refractivity
0.09
53.35
20.85
93.76
WPR Wind 0.57
GPSPW Precipitable Water 0.14
GNSSRO COSMIC-1,Metop-A,B,FY-3C/D Refractivity 0.11
AMVs-IR/VIS/WVs FY-2E,GOES-13,GOES-15,MTSAT-2,METEOSAT_10 uv 0.36
Rain Rain(nudging) 18.29

表3中可以看到,全球预报系统中卫星资料应用取得了大幅度进步,这在后文中有进一步叙述。就表4中中尺度预报系统的资料应用情况而言,虽然雷达等高时、空分辨率资料的应用有了长足进步,但在卫星资料应用方面仍然需要加强。

2.2.2 业务体系概况和应用效果

到2018年为止,基于GRAPES成果建成了高、低分辨率搭配适当、确定性和集合预报有机结合的完整数值天气预报体系。与基于引进技术的体系相比,自主业务体系的建立不仅需要模式和同化研发成果向业务转化的大量工作,而且涉及从观测资料检索、前处理和监控,直至预报数据的后处理、检验、可视化等各个环节的重构或者研发。在GRAPES业务体系建设过程中,除模式、同化和观测资料应用方面的显著进步外,在大数据库、试验平台、分析诊断工具和检验平台等支撑系统的研发和应用方面也取得了丰硕成果。

GRAPES业务体系由0.25°×0.25°水平分辨率垂直方向60层和0.5°×0.5°水平分辨率垂直方向60层的全球确定性和集合预报、亚太区域0.09°×0.09°垂直方向68层的中尺度确定性预报系统和覆盖中国区域0.03°×0.03°垂直方向50层的高分辨率数值天气预报系统以及覆盖中国区域0.1°×0.1°垂直方向50层的中尺度集合预报系统构成。同时,基于GRAPES驱动的全球/区域海浪预报系统、东亚沙尘暴预报系统和核污染应急预报系统作为专业数值天气预报系统也是GRAPES体系的重要组成部分。表56分别给出了GRAPES确定性预报系统和集合预报系统的配置情况。

表 5  GRAPES确定性业务预报系统配置 Table 5  Configuration of GRAPES deterministic forecast system
GRAPES_GFS  
全球预报系统  
亚太区域中尺度预报系统  
GRAPES_TYM  
中国千米尺度高分辨率预报系统 GRAPES_Meso 3 km
分辨率 0.25°×0.25°L60 0.09°×0.09°L68 0.03°×0.03°L50
模式区域 全球,模式顶 3 hPa 40°E−180°,15°S−60°N,模式顶 10 hPa 70°−145°E,10°−65°N,模式顶 10 hPa
预报时次 00:00、06:00、12:00、18:00 UTC 00:00、06:00、12:00、18:00 UTC 00:00、06:00、12:00、18:00 UTC
预报时长 240 h 120 h 36 h
同化系统 4DVar 全球分析降尺度+涡旋初始化 全球分析降尺度+云分析
物理过程 1. 双参数云物理方案
2. 宏云物理与预报云方案
3. NSAS 积云对流参数化
4. RRTMG 长短波辐射
5. CoLM 陆面过程
6. NMRF 边界层参数化
7. 次网格地形重力波拖曳
8. 小尺度地形湍流拖曳
1. Meso-SAS 积云对流参数化
2. WSM-6 云物理方案
3. RRTM 长波/Goddard 短辐射
4. NOAH 陆面过程
5. YSU 边界层参数化和垂直扩散
1. WSM-6 云物理方案
2. RRTMG 长短波辐射
3. NOAH 陆面过程
4. NMRF 边界层参数化和垂直扩散
表 6  GRAPES集合预报系统配置 Table 6  Configuration of GRAPES ensemble forecast system
GRAPES_GEPS 全球集合预报系统 区域中尺度集合预报系统 GRAPES_REPS
模式区域 全球,模式顶3 hPa 70°−145°E,10°−65°N,模式顶 10 hPa
预报时次 00:00、06:00、12:00、18:00 UTC 00:00、06:00、12:00、18:00 UTC
预报时长 360 h 84 h
同化系统 全球 4DVar 升尺度 全球分析降尺度+云分析
初值扰动 奇异向量 ETKF
模式扰动 SPPT、SKEB SPPT
侧边界条件 GRAPES_GEPS
集合成员数 30 15

应该指出,过去的20 a是国际上数值天气预报持续高速发展,且各国纷纷加强数值天气预报投入的时期。中国气象局数值预报系统尽管取得了很大的成绩,但与国际先进业务系统相比,GRAPES业务体系仍存在一定的差距,包括GRAPES全球预报系统的分辨率较低、所用资料同化系统尚未考虑集合信息、缺乏云雨区卫星资料同化和卫星资料变分偏差订正等先进技术、尚未建立对流尺度集合预报系统、缺乏对流尺度同化系统和千米尺度分辨率快速循环系统以及集合预报的成员数较少且分辨率较低等。这些问题是一个全新的数值天气预报业务体系在形成之初不可避免但必须认真对待的问题,也是今后加快补缺的重点。

从业务GRAPES_Meso系统降水预报评分的逐年提高情况(图1)可以看出,GRAPES_Meso在各量级降水预报方面均有长足的进步,成为中央气象台逐日天气预报不可或缺的重要支撑。GRAPES_Meso的进步包含了分辨率的提高、模式动力框架和物理过程改进以及有效应用雷达资料的云分析系统的引入。2019年6月,水平分辨率3 km、覆盖全中国的GRAPES_Meso系统正式上线,该系统对强降水的预报能力将GRAPES业务应用水平推上了新的台阶。图2给出2019年6—8月平均的GRAPES_Meso 3 km预报的逐3 h暴雨预报的ETS评分,为了比较也给出了上海SMS-WARMS(9 km×9 km分辨率)以及ECMWF(水平分辨率9 km×9 km)的预报ETS评分。从图2可以看到,GRAPES_Meso 3 km在预报前20 h表现出了明显的优势。需要指出的是,不仅是前20 h的暴雨预报,对于其他时段降水的预报ETS评分也表现出明显的优势。

图 1  GRAPES_Meso从2006年7月业务化以来24 h预报逐月TS演变及变化趋势 Fig. 1  Monthly-averaged TS and its linear trend of GRAPES_Meso 24 h forecasts since its operational application in July 2006
图 2  GRAPES_Meso 3 km暴雨预报(≥50 mm) ETS评分 (2019年6月10日—8月18日的平均,红色柱:GRAPES_Meso 3 km,蓝色柱:上海9 km WRF模式,绿色柱:ECMWF的预报;数据来自中国气象局业务内网idata.cma/areaHighResolution) Fig. 2  ETS skill scores of 3 h heavy rainfall forecast by GRAPES_Meso 3 km version (red color). Also shown are the forecasts from Shanghai WRF and ECMWF model. This figure is based on data from the CMA internal network: idata.cma/areaHighResolution

图3给出了区域模式GRAPES_TYM 2019年业务版本(V3.0)对台风路径和强度回报的误差随时间的演变。图3abc分别为对2016—2018年所有西太平洋台风回报的平均路径误差、中心最低气压平均误差和最大风速平均误差5 d预报的演变情况。为了比较,在图中也给出了预报员常参考的ECMWF和NCEP全球模式的预报误差情况。从路径预报来看,目前仍是ECMWF全球预报占优,GRAPES_TYM与NCEP全球预报的结果相当。从强度预报来看,无论是中心最低气压平均预报误差或者是最大风速平均预报误差,GRAPES_TYM均表现出了比较明显的优势。

图 3  GRAPES_TYM 2019年业务版本(V3.0)对台风路径和强度回报的误差随时间的演变 (a、b、c分别为对2016—2018年所有西太平洋台风回报的平均路径误差、中心最低气压平均误差和最大风速平均误差5 d预报的演变情况;蓝线:GRAPES_TYM、绿线:NCEP、红线:ECMWF) Fig. 3  Time evolutions of forecast errors for TC track,center pressure and maximum wind by the 2019 operational version of GRAPES_TYM (a,b,c are for forecasts of all the typhoons over the western Pacific during 2016—2018;blue line: GRAPES_TYM,green line: NCEP,red line: ECMWF)

图4是2010年以来GRAPES_GFS预报第5天500 hPa高度场距平相关系数(Anomaly correlation coefficient,ACC)的时间序列,同时在图中也给出了ECMWF、NCEP的预报结果。可以看到,GRAPES_GFS对形势场的预报技巧10年来逐步提高,虽与ECMWF、NCEP还有差距,但逐年预报技巧的不断提高反映了GRAPES_GFS模式、同化以及观测资料应用等诸方面的巨大进步,这将在第4节中进一步说明。

图 4  2010年以来GRAPES_GFS预报第5天500 hPa高度场距平相关系数时间序列 Fig. 4  Time series of ACC of 500 hPa height on the 5th day forecast by GRAPES_GFS since 2010

特别是,2019年汛期以来,GRAPES_GFS对中国强降水的预报效果多次超过ECMWF,对江南、江淮、江汉大型暖切变等导致的系统性降水(2019年5月26日),GRAPES_GFS(TS评分0.29)性能明显优于ECMWF(TS评分0.16)。另外,在江淮、江南、江汉、华南冷切变强降水(2019年6月6—7日、12日)GRAPES_GFS的性能也明显优于ECMWF,为预报员提供了重要的参考。图5是2019年5月24日20时(北京时)起报的36 h预报降水与观测的散点分布。可以清楚地看出,此次强降水过程的预报GRAPES_GFS明显好于ECMWF。

图 5  降水(单位:mm)预报与观测的散点分布 (a. GRAPES_GFS,b. ECMWF) Fig. 5  Scatter plots of precipitation (unit:mm) forecasts and observations (a. GRAPES_GFS,b. ECMWF)

目前,对于确定性预报而言,GRAPES_GFS与GRAPES_Meso的预报产品形成了很好的搭配,GRAPES_GFS把握形势场和雨带的中期变化,3 km分辨率的GRAPES_Meso则重点为预报员提供强降水雨量和落区的预报。这两个系统目前已经成为中央气象台逐日预报的核心支撑。

同样地,GRAPES全球和区域集合预报的业务应用为GRAPES“大家庭”增加了丰富的概率预报产品,在提前一周预报大范围强降水概率及局地强天气定量、定点概率预报中发挥了重要作用。

3 GRAPES主要科学技术进展

在GRAPES模式和同化技术不断深化发展和业务应用过程中,取得了自主创新的成果。与GRAPES设计和研发之初相比(薛纪善等,2008),在模式动力框架的核心算法、变分同化框架、卫星资料同化技术、适合千米尺度分辨率的变分同化技术、集合预报初值和模式扰动方法等方面有了显著的进步。

3.1 GRAPES模式动力框架实现了从经典半隐式半拉格朗日算法到预估修正半隐式半拉格朗日算法的转变,模式计算精度和运算效率大幅度提高

GRAPES模式动力框架在设计和研发之初,采用了当时国际业务数值天气预报模式主流采用的两时间层半隐式半拉格朗日算法解决离散化后平流和快波计算的稳定性问题。两时间层半隐式半拉格朗日算法虽然在节省计算内存提高计算效率方面有明显优势,但需要解决线性化时采用等温参考大气、半隐式处理时权重系数选取和非线性项时间外推及上游点计算时中间点风速外推带来的计算精度和稳定性问题(Simmons,et al,1997Ritchie,et al,1996Temperton,et al,2001)。这些方面在GRAPES模式开始研发时没有很好地考虑。为此,借鉴UKMO等的成功经验,进一步在GRAPES模式中发展了预估修正半隐式半拉格朗日算法,并在线性化时引入三维参考大气(Côté,et al,1998)。在动力框架中将三维参考大气与预估-修正算法结合,减小非线性项数量级的同时,通过迭代算法避免时间外推带来的影响,也可以使半隐式系数取值接近0.5,时间积分算法接近二阶精度(苏勇等,2018)。这些新的改进大幅度提高了模式的稳定性、计算效率和计算精度。

半拉格朗日模式均面临的另一个重要问题是对诸如水汽等标量平流计算的正定、守恒。半拉格朗日算法通过寻找上游点进而通过上游点处物理量的插值获得下一时间步的预报。插值过程必然带来不守恒和一定程度的耗散。Bermejo等(1992)提出了准单调半拉格朗日算法(QMSL,Quasi-Monotone Semi-Lagrangian),简捷且易于实现,但在水汽的强梯度、不连续区计算精度较低,且不能做到守恒。Bermejo等(2002)随后在准单调半拉格朗日算法基础上提出了订正的守恒算法。这些方法由于简捷和计算高效,在半拉格朗日模式中应用较多,ECMWF的半隐式半拉格朗日谱模式中至今仍在应用(Diamantakis,et al,2014)。不少研究者也发展了精度高且严格守恒的“体元映射”型半拉格朗日算法,如Cell-integrated semi-Lagrangian advection scheme(CSLAM,Lauritzen,et al,2010)、Semi-Lagrangian Inherently Conserving and Effi-cient scheme(SLICE,Zerroukat et al,20042009)等。但是,由于这些方法计算代价过高,在业务数值天气预报模式中较少采用。针对这种状况,研发了基于分段有理函数的标量平流方案(PRM,Piecewise Rational Method)。PRM方案通过将拉格朗日形式的标量预报方程重新写成易于构造守恒格式的通量形式,结合单元格界面通量更新的半拉格朗日算法,实现了高精度守恒且高效的水汽平流计算(苏勇等,2013)。

3.2 研发了GRAPES全球四维变分同化系统,实现业务应用

自1997年ECMWF将四维变分同化成功业务应用并显示出在提高数值天气预报技巧方面的巨大潜力后,4DVar一直在研究和业务应用中被认为是先进的资料同化技术(Rabier,et al,2000Bannister,2017Bonavita,et al,2017Kwon,et al,2018)。由于4DVar中所需切线伴随模式开发的巨大工作量和难度,研发4DVar是一项挑战性的工作。GRAPES全球4DVar系统的研发始于2008年,至2018年7月才实现业务应用,历经10 a。在全球非静力GRAPES模式基础上,发展了非静力的切线和伴随模式。为保持在同化时间窗内切线性模式与非线性模式在演变轨迹上较高的近似程度,研发了基于简化物理过程的切线性物理过程,包括积云对流、云物理、边界层垂直扩散、次网格地形阻塞流拖曳。这些研发工作成为GRAPES全球4DVar系统的重要基础。如前所述,GRAPES全球4DVar采用增量分析方案,在代价函数中引入数字滤波作为弱约束以抑制高频重力波,在系统中设计了多重外循环(目前在业务中采用单重外循环),极小化算法采用Lanczos-CG算法(Liu,et al,2018)。4DVar分析框架采用了水平和垂直不可分离的背景场误差协方差模型,水平相关模型采用二阶自回归模型,相关尺度随高度变化,垂直相关模型则直接由集合样本统计得到。平衡约束考虑了旋转风和散度风、旋转风和质量场、非平衡散度风和质量场的平衡,采用动力和统计结合方案来实现。

GRAPES全球4DVar系统每天进行4次6 h的分析,同化时间窗的起始时刻分别是03:00、09:00、15:00和21:00 UTC。全球4DVar分析结束后生成同化时间窗起始时刻的分析场,然后做6 h全球模式预报,为下个时次的全球4DVar分析提供背景场,同时输出3 h的模式预报结果,也就是06:00、12:00、18:00和00:00 UTC四个标准时刻的分析场,作为10 d全球预报的初值。GRAPES全球4DVar业务同化系统使用的观测资料包括探空报的温度、风场、相对湿度,地面报的气压,船舶报的气压,飞机报的风场,云导风,洋面散射风,GNSSRO掩星折射率和NOAA15 AMSUA、NOAA18 AMSUA、NOAA19 AMSUA、METOP-A AMSUA、METOP-B AMSUA、ATMS AMSUA、AIRS、FY4A GIIRS、FY-3C MWHS等卫星辐射率资料(表3)。相对于3DVar来说,GRAPES全球4DVar增加了近50%的观测资料使用量。

在分析质量方面,GRAPES全球4DVar全面超过3DVar的结果,尤其是质量场和风场。从均方根误差的角度来看,GRAPES全球4DVar背景场(6 h预报场)甚至优于同时刻的全球3DVar分析场。从分区域统计结果来看,南半球位势高度场和热带高层风场的改进更为明显。在不同时效和不同变量的预报方面,GRAPES全球4DVar都有改进。从图6 500 hPa位势高度的距平相关系数来看,平均预报时效有望提高5 h,其中南半球能够达到7—8 h。对于降水预报来说,GRAPES全球4DVar能够明显提高中雨以上量级的预报技巧。在台风预报方面,路径和强度预报均有明显改进,其中路径预报误差能够减小15%左右。

图 6  500 hPa预报位势高度场的距平相关 (a. 北半球,b. 南半球) Fig. 6  ACC of 500 hPa geopotential height (a. Northern Hemisphere,b. Southern Hemisphere)
3.3 形成了系列卫星资料同化应用技术

自从NCEP和ECMWF分别于1995和1996年在业务三维变分同化中实现卫星辐射率资料(TOVS,The Television and Infrared Observation Satellite(TIROS)Operational Vertical Sounder)的直接同化后,卫星资料作为全球数值天气预报的重要资料源及其在提高数值天气预报技巧方面的重要作用一直受到重视(Eyre,2007)。ECMWF在此领域处于领先地位,其同化系统中所用观测资料超过90%来自卫星探测。不同于引进和本地化国外系统,GRAPES卫星资料同化技术的系统性研发始于2004年,在早期的GRAPES等压面三维变分同化系统中加入了NOAA-16微波温度计AMSU-A辐射率资料的直接同化(Zhang,et al,2004)。包括NOAA系列(NOAA 15、16、17、18)ATOVS(Advanced TOVS)、GPS掩星资料、静止和极轨卫星云导风等批量卫星资料同化技术的系统研发和业务应用始于全球中期预报系统GRAPES_GFS的开发(沈学顺等,2009)。在GRAPES_GFS准业务化、业务化过程中,陆续开发了中国风云系列、欧洲Metop等微波温度计AMSU-A、微波湿度计MHS、高光谱大气红外等辐射率(IASI、AIRS)资料,多平台卫星资料的成功同化是GRAPES_GFS预报能力持续提高的重要因素(沈学顺等,2015)。

图7给出的GRAPES系统中同化卫星资料的种类和平台数目逐年变化情况可以看到,近10年卫星资料的同化应用进展非常迅速,尤其是GRAPES_GFS 2.0版本2016年正式业务化之后(沈学顺等,2017)。目前GRAPES系统同化应用的卫星遥感资料占观测资料总量的约70%(表3),在提高GRAPES_GFS的预报能力方面显示出了重要作用。

图 7  GRAPES全球模式可同化的卫星资料种类逐年演变 (图的最右侧给出的是先进数值天气预报中心的卫星资料应用情况) Fig. 7  Changes of satellite data types and number of platforms assimilated by GRAPES_GFS (The right side of the figure shows the status at advanced NWP center)

上述卫星资料同化应用的进步得益于观测数据质量分析、偏差订正、云检测、通道选择、质量控制等关键技术的逐步掌握和成熟运用,当然,GRAPES同化框架本身的不断进步也至关重要。而且,根据GRAPES_GFS模式的特点,还发展了国际先进的有约束卫星资料偏差订正方案。传统的观测资料偏差订正主要依赖于观测和背景的差异(O−B),在业务应用中存在当模式有显著偏差时观测信息会被错误订正,以及当质量控制不完善时会造成不合理的观测偏差估计和订正等问题。这两个问题在GRAPES_GFS发展过程中表现得尤为凸出,主要表现在同化预报循环中,预报偏差与观测偏差订正会产生相互作用,从而导致观测偏差订正向模式偏差的漂移。借鉴数学物理反问题中正则化思想中的“极小模解”方法,把卫星辐射率资料定标和辐射传输模式不确定性估计通过正则化参数把不等式约束引入目标泛函,发展了有约束的卫星资料偏差订正技术—CBC(Constrained Bias Correction)(Han,2014Han,et al,2016),减小了模式的背景偏差对卫星辐射率资料偏差订正的影响,较好地去除了资料本身的系统性偏差,更好地利用了观测信息。

有约束卫星资料偏差订正技术是在发展GRAPES数值天气预报系统中提出的原创技术,在2016年被引入到ECMWF的同化系统(Han,et al,2016),发展为有约束变分偏差订正(Constrained Variational Bias Correction,CVarBC),克服了卫星观测偏差订正向模式平流层偏差的漂移,提高了平流层初值分析精度,并改进了全球中期预报效果。CVarBC在ECMWF的2018年业务版本(IFS CY46R1)中开始业务化应用(ECMWF,2018a),此技术还将用于ECMWF第6代大气再分析(ECMWF,2018b)。

随着中国风云气象卫星的不断进步,风云卫星观测资料在GRAPES业务系统中同化应用的数量和效益不断提高。同时,通过在GRAPES中的定量评估和同化应用,也不断发现不足和问题并反馈,促使风云卫星观测资料质量和精度不断提升,中国气象局业务中逐步形成了卫星资料同化应用与改进的全链条业务协作。中国第一代静止气象卫星FY-2的卫星云导风产品十几年来精度持续提高,在GRAPES全球和区域同化预报中发挥了重要作用(Han,et al,2006万晓敏等,20172018)。可见光和红外自旋辐射扫描仪(VISSR,Visible and Infrared Spin Scan Radiometer)是中国自主研发的第一代静止卫星FY-2上搭载的有效载荷,VISSR的晴空水汽通道辐射率资料实现了在GRAPES全球4DVar中的同化应用(王皓,2017),充分发挥了其高时间分辨率观测的优势。

为充分发挥中国第二代极轨气象卫星FY-3多平台多载荷观测的综合优势,弥补单一仪器通道信息不足的缺点,在定量应用中发展了综合云和降水检测算法,实现了FY-3卫星微波温度计(MWTS)、湿度计(MWHS)、成像仪(MWRI)和掩星观测(GNOS)在GRAPES业务系统中的同化应用(Li,et al,2016Wang,et al,2019)。

2016年发射成功的中国新一代静止气象卫星FY-4A,搭载了地球静止轨道干涉式探测仪(GIIRS,Geosynchronous Interferometric Infrared Sounder),在国际上首次实现了静止轨道高光谱探测(Yang,et al,2017)。在GRAPES同化系统中发展了GIIRS观测算子(Di,et al,2018)、基于观测区域背景误差协方差的通道选择技术(尹若莹等,2019),以及适合大阵型探测器的偏差在线估计和订正算法等关键技术(Yin,et al,2020),2018年12月实现了GIIRS辐射率资料在业务GRAPES全球4DVar中的同化。同时,充分利用FY-4A探测仪机动灵活观测能力和数值天气预报系统的云预报和敏感区识别技术,在国际上首次实现了“智能快速晴空业务观测模式”和“高影响天气目标观测模式”。2018年在实时业务环境下,针对台风“玛丽亚”“安比”和“山竹”,应用GRAPES奇异矢量识别的敏感区,对卫星观测区域、观测频次等参数进行最优确定,通过观测和预报的协同互动,实现了面向预报对象的静止卫星探测模式,并同化到GRAPES全球4DVar,改进了对台风路径的预报。

3.4 研发了适合千米尺度分辨率数值天气预报的三维变分同化技术和云信息初始化技术

发展千米尺度资料同化系统,改善模式初值中的中小尺度信息,对于提高极端天气和近地面气象要素的预报能力十分重要(Gustafsson,et al,2018)。GRAPES千米尺度资料同化系统基于GRAPES全球区域一体化变分同化框架(薛纪善等,2008)开发,采用三维变分同化(3DVar)选项(以下简称为km-scale 3DVar),框架主体、观测算子、资料流程等与区域3DVar系统保持一致。

在同化分析框架上,km-scale 3DVar面向千米尺度数值天气预报的需求做了针对性开发(Wang R C,et al,2018)。首先,km-scale 3DVar将极小化目标泛函中的动力学控制变量由流函数(ψ)、势函数(χ)、以及π转换为u风、v风以及温度(T)和地面气压(ps),并且不再考虑控制变量间的地转平衡关系。新控制变量对应的背景误差空间相关尺度更小,且为观测系统的直接测量要素。如图8所示,与业务10 km 3DVar相比,km-scale 3DVar(试验水平分辨率为3 km)的分析增量要更加局地,适合中小尺度系统的同化分析。其次,为了能在分析好中小尺度系统的同时兼顾天气尺度环境气流的正确描述,还研发了多尺度分析方案,包括采用多尺度水平相关模型(吴洋等,2018)、取千米尺度分析场中的中小尺度信息与全球分析场中的大尺度信息进行重组的直接混合,以及在km-scale 3DVar的目标泛函中引入大尺度信息作为弱约束的隐式混合(Yang M J,et al,2019)等。此外,为适应中小尺度信息发展迅速、生命史短的特征,也为弥补3DVar时间维的缺失,km-scale 3DVar采用快速同化预报循环更新的方式运行,同化分析间隔为1—3 h,采用数字滤波抑制快速循环过程中的噪音问题。

图 8  单点探空温度观测 (位置:31.93°N,118.9°E,850 hPa) 分别在GRAPES区域10 (a) 和3 (b) km变分同化系统中强迫出的温度分析增量 (沿32°N的垂直剖面,等值线单位:K) Fig. 8  Temperature analysis increments forced by single-point sounding temperature observations (location: 31.93°N,118.9°E,850 hPa) by 10 (a) and 3 (b) km variational assimilation systems in GRAPES_Meso (vertical profiles along 32°N,isoline unit: K)

在观测资料应用方面,km-scale 3DVar除了同化现有业务系统所用观测资料,还积极发展时、空稠密资料的同化应用。首先,在雷达资料的应用上,径向风和风廓线雷达资料已在km-scale 3DVar中实现直接同化应用,改进了流场分析;反射率资料则在云分析系统中用于诊断水凝物和潜热倾向信息,并通过张驰逼近的形式引入到模式轨迹中。其次,对中国新一代地球静止卫星FY-4A提供的高时、空分辨率观测资料,km-scale 3DVar