气象学报  2019, Vol. 77 Issue (2): 165-179   PDF    
http://dx.doi.org/10.11676/qxxb2019.029
中国气象学会主办。
0

文章信息

刘艳, 薛纪善. 2019.
LIU Yan, XUE Jishan. 2019.
GRAPES的新初始化方案
The new initialization scheme of the GRAPES
气象学报, 77(2): 165-179.
Acta Meteorologica Sinica, 77(2): 165-179.
http://dx.doi.org/10.11676/qxxb2019.029

文章历史

2017-11-21 收稿
2018-11-05 改回
GRAPES的新初始化方案
刘艳1,2, 薛纪善3     
1. 中国气象局数值预报中心, 北京, 100081;
2. 国家气象中心, 北京, 100081;
3. 中国气象科学研究院, 北京, 100081
摘要: 四维变分同化由于引入预报模式作为一项约束,理论上它的分析场已经具有较好的平衡性,但实施时还会有诸多因重力波导致的高频振荡过程,因此,四维变分同化(4DVar)分析仍需要初始化。文中描述了GRAPES全球四维变分同化系统(GRAPES-4DVar)的新初始化方案的科学设计、公式演绎以及试验结果。GRAPES-4DVar的新初始化方案采用数字滤波方案作为代价函数的一项约束控制重力波引发的不平衡结构,约束强加在分析增量上与极小化迭代过程同步进行。新的初始化方案是变分同化系统的一部分,数字滤波的积分时间与4DVar的同化时间窗一致,不会对4DVar产生额外的计算资源消耗;并能适应长时间窗的同化,不会因为时间窗的延长而削弱慢波过程。新初始化方案中,模式轨迹的光滑程度可在变分同化中通过重力波控制项的权重系数方便控制。GRAPES全球四维变分同化的理想和循环同化批量试验都表明,在四维变分同化中,重力波的控制依然非常重要,具有初始化的GRAPES试验,无论分析还是预报技巧都较无初始化的有明显优势。与以前分析和滤波独立实施的旧初始化方案相比,新方案的分析和预报效果略优,同时有效地节省循环同化系统的运行时间,这对四维变分同化来说非常重要。
关键词: 初始化     数字滤波     四维变分同化     GRAPES    
The new initialization scheme of the GRAPES
LIU Yan1,2, XUE Jishan3     
1. Numerical Weather Prediction Center, China Meteorology Administration, Beijing 100081, China;
2. National Meteorological Center, China Meteorology Administration, Beijing 100081, China;
3. Chinese Academy of Meteorological Sciences, Beijing 100081, China
Abstract: Theoretically, the analysis of a four-dimensional variational (4DVar) data assimilation system is in better balance as the forecast model is its constraint. In actual implementation of 4DVar, however, there are many processes that can lead to high-frequency gravity waves. Thus, the initialization process is still necessary for 4DVar. In the present paper, the scientific design and formula derivation of the new initialization scheme of GRAPES (Global and Regional Assimilation Prediction System) numerical weather prediction system and experiments using the new scheme are presented. The four-dimensional variational data assimilation system of GRAPES (hereafter GRAPES-4DVar) adopts the digital filter scheme for initialization. The digital filter is added to the penalty function as a weak constraint to efficiently diminish the unbalanced structure associated with gravity-inertia waves. The constraint is only imposed on the incremental analysis and run synchronously with the minimization. The weak constraining term does not add extra computation cost since the integration period of the digital filter is the same as the 4DVar's time window. Furthermore, the new initialization scheme is suitable for long time window assimilation and will not weaken slow wave processes as the time window extends. It is easy to control the smoothness of model trajectory by tuning digital filter coefficient in the variational assimilation procedure. Results of GRAPES-4DVar assimilation cycling experiments show that the initialization is necessary in the 4DVar analysis, and the analysis and forecast skills have been improved significantly in the experiments with initialization compared to those without initialization. Compared to the old initialization scheme, in which the analysis and filter are two independent procedures, the new scheme efficiently saves the running time of assimilation cycling, although its impact on analysis and forecast is similar to that of the old scheme.
Key words: Initialization     Digital filter     Four-dimensional variation assimilation     GRAPES    
1 引言

资料同化过程所产生的数值预报初始场,不仅应与观测信息及模式提供的背景场充分接近,还应反映所关注的天气过程所具有的动力特性,包括维持气象要素间应有的准平衡关系。如不顾及后者,单纯考虑观测与背景场的信息融合所得到的客观分析场在时间积分过程中往往会产生具有比所感兴趣的天气过程更高时间频率的“噪音”,这些噪音带来预报误差,甚至妨碍预报模式积分的正常进行(纪立人,2011)。为了过滤掉这些“噪音”,需要对客观分析结果做初始化处理,这是当代资料同化的一个重要环节,即初始化(Daley, 1997; Lynch, et al, 2010)。实质是通过对初始场的调整,控制模式积分过程中重力波的快速增长。

早期的初始化是资料同化中的一个独立阶段,它与从观测资料得到关于模式预报变量的分析场(即通常所谓的客观分析)是分别实施的。随着现代变分资料同化方法的发展与广泛应用,特别是在四维变分同化中,预报模式本身作为同化过程一个约束条件,同化系统产生的分析场已经具有较高的平衡特性(Courtier, et al, 1994; Rabier, 2005; Rawlins, et al, 2007)。从理论上讲,四维变分同化所寻求的是相空间中与一段时间(即整个同化时间窗)的观测最接近的模式运行的轨迹,在理想的情况下,这条轨迹不再包含高频的噪音。但实践表明,以四维变分同化的结果作初始场做出的预报并不能完全排除预报过程出现高频噪音。这是因为数值预报模式本身允许各种时、空尺度的运动共存,而在四维变分同化的实施过程中,还存在众多会引入各种噪音与误差的环节,如通行的内、外循环的粗、细网格间的反复插值与转换(Courtier, et al, 1994)等。因此,即便对于采用四维变分方案的同化系统,初始化或重力波的控制依然是必要的。

重力波的控制可以通过调整初值中的重力波分量来实现,但不能采用对初值作正规模态的分解后消除其中高频重力波分量的简单做法,而必须考虑快波与慢波间的非线性作用,这就是非线性正规模态的初始化方法(Williamson, et al, 1981; Temperton, 1988, 1989; Ballish, et al, 1992)。这种方法的基础是对初始场作正规模态分解,涉及模式的正规模态计算,不仅计算量较大,并且某些情况下,正规模态的计算很困难。另一种初始化方法是对模式积分的时间序列作滤波处理,消除其中的高频部分而保留低频部分,并将初始时刻的低频分量作为模式的初值,这就是数字滤波初始化(Digital Filter Initialization, DFI)方法(Lynch, et al, 1992)。这种方法的动力学理论依据与非线性正规模态是相同的,但更易实施,因而被广泛采用。

考虑到数字滤波作初始化的良好效果,以及容易实施等优点,GRAPES业务三维变分同化系统将数字滤波选为初始化方案(薛纪善等,2008b;薛纪善等)。尽管理论上重力波的控制是通过在变分的代价函数中增加反映重力波振幅项来实现,初始化已成为变分同化的一部分,但由于三维变分同化并不涉及模式积分,而数字滤波初始化则是对模式积分所产生的预报变量时间序列进行滤波,即前者只涉及初始时刻的要素场的处理,后者则属于纯粹的时间序列处理而不考虑空间点间的联系,三维变分分析与数字滤波初始化在计算上并不交叉,因此在实际实施时依然可以将三维变分分析与初始化作为两个阶段。随着GRAPES同化系统的发展(张华等,2004庄世宇等,2005薛纪善等,2008a马旭林等,2009;薛纪善等刘艳等,2016),GRAPES三维变分分析(3DVar)被四维变分分析(4DVar)所替代(薛纪善等,2008b张林等,2017刘永柱等,2017;Liu et al, 2019)。四维变分同化可以考虑观测资料的时间分布对分析场的影响,并且使背景误差在一定程度上呈现流依赖特征,因而具有更好的效果。四维变分同化引入预报模式作为一项约束,如前所述,它所寻求的是相空间中与整个同化时间窗的观测最接近的模式运行的轨迹,而轨迹的光滑程度自然影响到与观测的距离,即数字滤波初始化与四维变分分析是互相影响或渗透的。这意味着初始化不应再作为一个跟随变分同化之后的单纯的时间序列的处理过程,需要建立一个新的四维变分分析与数字滤波相结合的初始化方案。

薛纪善,刘艳,张林等. 2012. GRAPES全球三维变分同化系统模式变量分析版科学文档.中国气象局数值预报中心内部技术手册.北京:中国气象局,1-105

四维变分同化与数字滤波初始化的结合有两种方式,其一是以数字滤波得到的光滑的轨迹作为四维变分分析的强约束条件,其二是数字滤波只作为四维变分分析的弱约束(Gauthier, et al, 2001; Polavarapu, et al, 2000; Wee, et al, 2004)。在前一种情况下,四维变分分析所求得的是一个初始化的结果,它与同化时间窗内的观测整体上最靠近。在后一种情况下,四维变分分析的模式运行轨迹与观测的差异以及经过数字滤波的光滑轨迹的差异总体最小。参照国际上一些主要数值预报业务中心的做法(Gauthier, et al, 2001; Polavarapu, et al, 2000),GRAPES选择了后者。本研究介绍了这一新的初始化方案的科学设计与试验结果。

2 GRAPES-4DVar与数字滤波初始化的回顾 2.1 GRAPES-4DVar

四维变分同化使如下的目标泛函达到极小

(1)

式中,δx为分析增量,MH分别为预报模式与观测算子的切线性算子,BR分别为背景场与观测误差的协方差,d为新息向量,向量或矩阵的上标“T”表示其转置,而右下标表示时间(文中所有公式的变量符号的更完整定义见薛纪善等(2008b))。引入基于背景误差协方差矩阵的预条件变换等变量变换

(2)

式中,P为物理变换算子,σu为由格点上均方根构成的对角矩阵,U为预条件变换算子,w为目标函数极小化的控制变量,变换矩阵按B=PBPT=uUUTσuPT来定义(薛纪善等,2008b)。

再将预报模式展开

(3)

并定义

M0, 0为单位矩阵。

(4)

得到

(5)
(6)

在导出式(6)时,已对相应项做了重新归并。由于目前的模式初值依然采用静力平衡假定,x0并未包括所有的模式预报变量,对于x0中未包括的模式预报变量则通过静力关系,由已在x0中指定的预报变量获取。以xm表示模式的完整预报变量,则

(7)

式中,D是根据δx0诊断δx0m的算子,D-1并不表示严格的矩阵逆。显然模式积分须以δx0m为初值,故式(4)和(6)表示模式积分的部分均应按(7)插入算子D

(4')
(6')
2.2 数字滤波初始化

数字滤波初始化可以简单表示为

(8)

式中,上标“F”表示经初始化的值,αk为数字滤波系数。关于αk的取值,见Lynch等(2010)、薛纪善等,这里不再重复。

薛纪善,刘艳,张林等. 2012. GRAPES全球三维变分同化系统模式变量分析版科学文档.中国气象局数值预报中心内部技术手册.北京:中国气象局,1-105

以预报模式代入

(9)

式中,表示预报模式,以与前面的切线性预报模式Mk-1, k相区别

(10)

数字滤波可以以强或弱约束的方式施加于四维变分同化,文中采取弱约束,即在四维变分同化的目标泛函中增加重力波控制项

(11)

式中,(*,*)表示内积, 也可按照向量运算习惯写成

(12)

式中,Q是一个正定矩阵。根据研究,这一约束只需在同化时间窗的一个特定时刻考虑(Gauthier, et al, 2001),文中将其取在同化时间窗的中间时刻

(13)

以数字滤波的公式代入

(14)

(15)

由于

(16)

故数字滤波公式可展开为

(17)

代入上面关于bI/2的表达式

(18)

式中,

(19)

在四维变分同化中需要计算Jcx0的梯度。考虑对一个已经经过滤波的参考态(如背景场)的扰动量δx0所产生的bI/2的变化(经线性化)为

(20)

Jc的一阶变分

(21)

(22)

注意,式中δx0是相对已经经过滤波的参考态,即背景场的扰动量δx0。需要说明的是,本节中至今并未考虑式(7)中δx0mδx0的区别,为了与模式一致,须将式(20)和(22)修正为

(20')
(22')
3 四维变分同化与数字滤波初始化的结合

JvarJc相加得到四维变分同化与数字滤波初始化相结合的目标泛函

(23)

式中,βvβc分别是观测(含背景场)项与重力波控制项的权重。在本节推导中暂将其取为1。假定背景场不包含噪音,相当于已经过数字滤波,则数字滤波的代价函数中只有分析增量的贡献,故

(24)

式中,

(25)

前面已经推导出Jvar对控制变量及Jc对分析增量的梯度分别为式(6)与(22)。为了合并计算,需要将变量统一,故根据δx0=uUw,将后一梯度转化为

(26)

式中,包含多个预报进程,不便于编程计算,参照四维变分同化的推导重新排列为

(27)

合并后目标泛函关于同化控制变量的梯度为

(28)

此式与四维变分同化的计算相比,在伴随模式反向积分的起始(强迫)项增加了数字滤波贡献项γi-I/2QbI/2

4 实施中的问题 4.1 内、外循环问题

与三维变分同化不同,四维变分同化因计算量过大,一般都需要采用两套或多套分辨率的增量分析策略。简言之,由分析估算的观测相当量在高分辨循环中计算,最优化则在低分辨循环中计算。前者称为外循环,后者为内循环。关于内、外循环问题,在GRAPES-4DVar的相关文献中已有叙述(薛纪善等,2008b),在此不再重复。这里只就Jc、∇Jc的计算做简单讨论,它们只在内循环涉及。

按照增量分析的公式,在内循环过程中对背景场的增量

(29)
(30)

注意,式中δx0是内循环中对本次外循环猜值的增量,它是相对于本次外循环的猜值内插到低分辨率后的值的增量。简单起见,略去表示完整模式变量的上标“m”。切线性模式的积分是以δx0为初值。因此

(31)
(32)

式(31)右端第二项在最优化(内循环计算)中是不变的。而实际上如进一步假定外循环的猜值也是经过滤波的,则这一项可以忽略,此时数字滤波只在内循环进行,但为了控制内、外循环的分辨率转换带来的噪音,在最后一次外循环后又进行了一次数字滤波的初始化(Gauthier, et al, 2001)。

外循环的主要计算包括:(1)积分非线性预报模式,得到各个时刻的背景场(首次迭代)或猜值场(其后的迭代),以及各时刻猜值场与背景场的差(xig-xib);(2)利用非线性观测算子计算各时刻背景场(猜值场)的观测相当量;(3)计算各时刻的新息量与猜值余差。

内循环的计算是在猜值固定的条件下求使目标函数达到极小的分析增量近似值(注意:这里所说的分析增量都是相对背景场,而非相对内循环的猜值,但内循环的切线性模式的积分是以相对内循环的猜值的扰动量作为初值的),所涉及的计算主要是目标函数及其梯度值。为此要先积分切线模式求得各观测时刻与w相对应的分析的增量的猜值(后面直接称分析增量),并计算相对于猜值的增量对应的观测增量,计算余差,计算目标函数,再反向积分伴随模式,每一步都要加上重力波控制的贡献,而遇到有观测的时刻,还要加上该时刻的观测项贡献。根据返向积分计算目标函数的梯度,做最优化计算求w的极小值。由于整个计算流程与四维变分同化完全相同,只是在正、反向积分中增加了重力波控制的贡献项,因此初始化所增加的计算机代价是不大的。

4.2 内积的计算

式(11)与(12)没有给出内积(*, *)与矩阵Q的具体定义,它们反映了各个模式变量在重力波控制中的贡献,根据对初始场平滑性的不同要求,可以有不同的选择,如希望各个要素都具有相似的光滑性,可选用“总能量范数”(Ehrendorfer, et al,1999; Errico,2000), 其定义为

(33)

式中,ρ为密度,g为重力加速度,δuδvδθδp分别表示纬向风速、经向风速、位温、气压的扰动变量,φ为纬度,λ为经度,z为模式高度,A为总区域的面积,下标dom表示区域求和。这里带“~”的变量均指参考大气的值,其中

(34)

式中,R为干空气气体常数,cp为定压比热,cv为定容比热,它们的值分别为287.04 J/(kg·K)、1005.7 J/(kg·K)和718.0 J/(kg·K)。

另一种选择是

(35)

式中,σu2σv2σθ2σπ2分别是纬向风速(u)、经向风速(v)、位温(θ)、无量纲气压(π)的背景场误差方差(Wee, et al, 2004)。

当前GRAPES-4DVar中矩阵Q定义为背景场误差方差。

5 试验结果 5.1 有无初始化的试验

图 1给出了变分同化分析中有、无初始化弱约束项时,观测、背景以及弱约束这3个目标泛函随迭代过程的演变。本节的试验是个例试验,采用了2013年5月1日03时(世界时,下同)的背景场,GRAPES-4DVar同化系统的外循环分辨率是0.5°,内循环水平分辨率1.0°,垂直分辨率60层,只同化了探空、船舶、地面、云导风、飞机观测等常规观测。根据最佳线性无偏估计理论,资料同化应使分析场与观测及背景场的平均距离最小,而当加入数字滤波初始化弱约束项时,分析的结果使得分析场与观测、背景场以及约束条件三者的平均距离达到最小。可以看出,在没有任何约束时,观测项Jo的目标泛函与有约束项时的目标泛函相当,背景场项Jb的目标泛函比有约束时略有增加,但约束项Jc的目标泛函比有约束项时的目标泛函明显增加。图 1说明了数字滤波初始化约束条件在保证分析对观测的拟合程度的情况下,对分析场做了调整和适应,既控制了四维变分同化的分析噪音,又改善了极小化的收敛性能。

图 1 极小化过程中无(a)和有(b)数字滤波弱约束项时的观测、背景以及约束项的代价函数的变化 Fig. 1 Evolutions of cost functions in the minimization for experiments (a) without and (b) with the digital filter constraint

利用图 1试验的分析场提供给模式作为初值积分,进一步了解有无初始化处理的模式初值对模式预报的影响。图 2是随机给出的格点(30°N,124°E)的地表气压随时间的变化曲线,其中黑线表示未做初始化的地表气压的变化,红线表示做了弱约束数字滤波初始化的地表气压的变化。可见在模式积分过程中,尤其在前几小时,模式激发的不合理的重力波扰动可以被初始化过程过滤掉,地表气压的预报曲线比较平稳和光滑,而没有初始化时,地表气压的预报振荡比较明显。这说明该初始化方案对抑制高频振荡是有效的。

图 2 格点(30°N,124°E)地表气压的时间演变,黑线(红线)表示无(有)初始化的模式积分结果 Fig. 2 Variations of surface pressure with the time step in the grid point of 30°N, 124°E for two forecasts. The black line indicates the forecast using the uninitialized analysis, and the red line shows the forecast after the digital filter constraint initialization

高频振荡表现为大气快速的辐散、辐合,地表气压倾向则代表了整层气压的辐散、辐合,因而地表气压倾向可以作为衡量高频振荡的一个重要指标(Lynch, et al, 1992)

(36)

式中,IJ分别表示区域的格点数。根据式(36)计算上述试验的地表气压倾向。图 3给出了第一步模式积分时地表气压倾向的地理分布,即每10 min地表气压倾向的变化。可见如果分析场没有经过初始化,会产生许多高频重力波,地表气压倾向为-2—3 hPa/(10 min),表现为等值线散乱,尺度很小,且密集(图 3a),而初始化后,地表气压倾向的最大、最小值分别为0.9和-1.2 hPa/(10 min),闭合等值线的区域和范围均明显减少,有效抑制了积分初期变量不协调产生的高频重力波(图 3b)。

图 3 无(a)和有(b)初始化过程的试验在2013年5月1日03时模式第一步积分时的地表气压倾向水平分布(单位:hPa/(10 min)) Fig. 3 Distributions of surface pressure tendency for experiment (a) without and (b) with Jc initialization in the first integration time step at 03:00 UTC 1 May 2013 (unit: hPa/(10 min))

初始化通过改变初值减少在积分初期的噪音,提高积分初期的预报质量。但随着积分时间延长,无初始化的模式变量开始协调,这时无初始化和有初始化的预报结果应该趋于一致。因而初始场的改变不能太大,否则必然会严重影响预报结果。图 4给出了有、无初始化的第40层模式纬向风及其差异的水平分布,图 4a中等值线表示无初始化时的情况,色阶分布表示做了初始化的情况。由图可见,两者的分布形态一致,等值线和色阶几乎重合,而图 4b上两个初值是有差异的,但各格点的偏差不是很大,偏差大的地方主要在地形陡峭,海、陆交界,或有资料的地区,也正是这些因素容易在模式积分初期激发重力波和高频噪音,因此,虚假振荡比别的地方大。

图 4 有、无初始化的模式纬向风(单位:m/s)初值在2013年5月1日03时的模式第40层的水平分布及其差异的水平分布 (a.等值线为无初始化,色阶为有数字滤波弱约束的初始化,b.两者的差异) Fig. 4 Distributions of zonal wind (unit: m/s) at model level 40 simulated by two experiments at 03:00 UTC 1 May 2013 (The two experiments are the same as those shown in Fig. 3. a, the contours show results of the uninitialized experiment and the shadings show results of the initialized experiment. b displays the difference between these two experiments)

理想试验表明,四维变分同化虽然寻求的是整个同化时间窗的观测最接近的模式运行的轨迹,在理想的情况下,这条轨迹不再包含高频的噪音,但实际上,以四维变分同化的结果作初始场做出的预报并不能完全排除预报过程出现高频噪音。因此,在四维变分同化系统中,初始化依然是必要和重要的。

5.2 数字滤波初始化权重的影响

初始化在变分同化的式(23)中可以取不同的重力波控制系数(βc),该系数反映了希望被过滤掉的高频噪音的程度。继续利用5.1节的试验研究重力波控制系数对分析和预报的影响(图 5)。由图可见,βc越大,模式起转时间则越短,但大的重力波约束系数同时也意味着增加了四维变分同化的模式运行轨迹与观测的差异,即分析场远离观测,而太小的βc又会使得分析场仍有较多的高频噪音残余。因此,合理的βc值对改进模式的起转有重要意义,尤其在模式积分的前2—3 h。

图 5 根据不同重力波控制系数得到的分析场作为模式初值积分得到的地表气压倾向绝对值随时间的演变 Fig. 5 Evolutions of the absolute value of surface pressure tendency as a function of integration time in experiments initialized with analysis fields with different gravity control coefficients
5.3 数字滤波初始化与四维变分同化结合和分离的两种初始化方案的对比

前面已经提及,在三维变分同化系统中由于不涉及模式的积分,分析与初始化是两个独立实施的部分,数字滤波初始化必须放在分析完成之后进行。而四维变分同化系统带有模式积分,因此,模式初值可以通过初始化与四维变分同化相结合的方式得到,也可以参照三维变分循环同化系统的流程得到。在上节试验基础上,图 6分别给出了数字滤波初始化与四维变分同化结合与分离时的两种初始化方案对循环同化的影响,图中曲线表示全球平均的地面气压倾向的绝对值随时间的变化,其中模式积分步长取600 s,每一间隔为10 min,积分6 h。由图可见,数字滤波初始化与四维变分同化结合(红线)与分离(蓝线)这两种初始化方式得到的模式初值,对分析和预报的作用相当,它们都可以抑制初始时刻风场和质量场不协调激发出来的重力波,以及观测资料带来的噪音,气压场扰动约由原来的0.68降至0.3 hPa/(10 min)左右,这里数字滤波初始化与四维变分同化结合的方式中重力波控制系数βc取值为1.0。作为对比,图 6同时给出了背景场以及不做初始化的分析场作为模式初值的6 h预报结果,可以发现背景场(绿线)是光滑的,但做完分析后(黑线),会带进观测资料的噪音,以及激发出高频重力波。一旦加入初始化后,能够有效减少高频噪音。一般说来,模式经过6 h左右的积分时间,即使不做初始化,高频噪音也基本消除。

图 6 不同初始场积分模式得到的地表气压倾向绝对值的时间演变 (红线表示数字滤波初始化与四维变分同化结合的分析场作为模式初值,蓝线表示数字滤波初始化与四维变分同化分离的分析场作为模式初值,黑线表示未做初始化的分析场做模式初值,绿线表示背景场作为模式初值) Fig. 6 Evolutions of the absolute value of surface pressure tendency as a function of model integration time using different initial fields (The blue line indicates the initialization field in which the analysis and filter are two independent procedures, the red line indicates the initialization field in which the analysis and initialization are coupled, the black line denotes the uninitialized analysis field, and the green line is the background)
5.4 对循环同化的分析和预报的影响

为进一步说明初始化对四维变分循环同化系统的分析和预报的影响,利用GRAPES全球四维变分循环同化系统,按照GRAPES业务系统的设置分别做了3组试验,第1组是没有初始化的试验,称之为NODF,另两组分别是数字滤波初始化与四维变分同化结合与分离的初始化试验,称之为JCDF和OLDDF。同化的观测资料有常规与遥感观测,包括探空,船舶、地面、云导风,飞机观测,掩星资料,美国国家海洋大气局和欧洲气象卫星组织的极轨卫星的微波与红外资料,以及红外高光谱资料。同化的时间窗长度为6 h,同化系统外循环水平分辨率为0.5°,内循环水平分辨率为1.5°。GRAPES全球预报模式为2015年底集成的试验版本,垂直方向60层,模式顶高35 km,水平分辨率为0.5°,时间积分步长600 s,所选取的物理过程参数化方案包括:云与降水的微物理方案(Hong, et al,2006)、长、短波辐射方案(Mlawer, et al,1997)、边界层参数化方案(Troen, et al,1986; Liu, et al, 2015)、积云对流参数化方案(Arakawa, et al,1974Liu, et al, 2015)等。试验时间从2013年5月1日03时到2013年5月31日21时,即一个月的循环同化试验。

图 7为这3组试验得到的模式初值与欧洲中期数值预报中心的全球再分析资料ERA-Interim(Dee, et al,2011)的比较,这里给出的是位势高度、温度、纬向风和经向风4个要素与ERA-Interim分析场的差的平均值与均方根的全球平均垂直分布廓线。仅就单次分析来说,有、无初始化并不会造成分析效果有太大的差异,但连续同化循环的累积作用,会使有初始化处理的同化效果明显优于无初始化的分析,减小了GRAPES变分同化系统与再分析资料的差异,这在高层更为明显。数字滤波初始化与四维变分同化结合和分离的两种初始化方式,对模式初值的影响效果相当。但以弱约束结合的初始化方式与最优化迭代过程同步进行,它的积分时间与四维变分的同化时间窗一致,因此,这种初始化方式没有额外的计算机资源耗散,不像数字滤波独立的初始化过程那样,先经过模式绝热过程的反向积分,再到非绝热过程的正向积分,来获得模式要素场的时间序列进行虚假惯性重力波的过滤。并且,随着同化时间窗的延长,数字滤波初始化更倾向于为一理想化的滤波器,它不会削弱慢波(Gauthier, et al, 2001)。因此,多数情况下,数字滤波初始化与四维变分同化结合的初始化方式分析效果略好于数字滤波初始化与四维变分同化分离。

图 7 GRAPES-4DVar的初值与ERA-Interim分析场均方根误差的全球平均垂直分布廓线 (a.位势高度,b.温度,c.纬向风,d.经向风;蓝线表示无初始化的初值,黑线表示数字滤波初始化与四维变分同化分离的初值,红线表示数字滤波初始化与四维变分同化结合的初值) Fig. 7 Root mean square errors of GRAPES analysis of (a) geopotential height, (b) temperature, (c) zonal wind and (d) meridional wind compared to ERA_Interim reanalysis (GRAPES minus ERA-Interim) (The blue curve is uninitialized; the black curve shows the initialization in which the digital filter initialization and 4DVAR analysis are independent; the red curves shows initialization with a weak constraint digital filter)

衡量初始化效果的好坏,不仅要求初始化有效地减少初始积分时产生的噪音,对初始场的改变程度在合理范围内,还要求初始化后的预报场质量不能降低。分别以上述3组试验的分析场作为模式初值进行全球中期预报试验,来判断初始化是否改进了初值质量,并采用目前通用的全球预报检验指标,即500 hPa预报与分析的高度场距平相关系数来衡量。图 8是这3组预报试验的距平相关系数随着预报时效的变化曲线。可见,无论在南半球还是北半球,初始化对改进数值预报都有正的贡献,各地区的预报时效均不少于6.5 d。经过数字滤波初始化,在南半球5 d以后的距平相关系数明显提高。

图 8 北(实线)、南(虚线)半球500 hPa的距平相关系数随预报时效的变化曲线 (蓝线表示无初始化的初值,黑线表示数字滤波初始化与4DVar分离的初值,红线表示数字滤波初始化与4DVar结合的初值) Fig. 8 Anomaly correlation coefficients (ACCs) of 500 hPa height in the 8 d forecast in the Northern and Southern Hemisphere (Colors indicate the same meanings as in Fig. 7)

初始化对降水预报有重要的影响。不经过初始化,在模式积分初期会产生虚假的重力波,导致虚假的大气辐合、辐散,从而影响垂直速度,产生不合理的降水。模式积分6 h后,有、无初始化的垂直速度已开始趋向一致,这表明此时初值作用已不明显,其后的降水差异主要是预报模式本身系统性的误差造成的。因此,考察了初始化对6 h累计降水的影响。从2013年5月1日03时不同初值预报的6 h累计降水的水平分布(图 9)看,不同初始化方案预报的降水落区差别不大,但预报的降水强度有差异。做了初始化的初值比无初始化的初值预报的降水强度要大,原因可能是初始化使模式起转时间更短,更早地产生合理的辐合、辐散场。但实际上,在这一个月的循环同化试验中,不同方案的初值在各时次预报的6 h累计降水量的全球平均值差别不大(图略)。这可能因为降水是一个复杂的物理过程,目前初始化方案仅调整了质量场和动力场,即对无量纲气压、位温以及风场调整,未调整湿度场。因此,关于初始化对降水预报的影响还需做更多细致的研究。

图 9 不同初始化方案在2013年5月1日03时预报的6 h累计降水(单位:mm) (a.无初始化,b.数字滤波弱约束初始化,c.数字滤波和同化分析相分离的初始化,d.弱约束初始化与无初始化时的降水差异) Fig. 9 Distributions of 6 h accumulated rainfall in three experiments (unit: mm) (a. no initialization, b. Jc initialization, c. initialization in which DFI and 4DVar are independent, d. difference between (b) and (a) ((b) minus (a)))
6 小结与讨论

初值的形成是数值预报的基础,只有在合理的、协调的初值条件下才能提高模式预报的准确性和稳定性。在数值预报中,如果客观分析的结果不满足动力学平衡而被用来当作模式积分的初值,会激发无气象意义的高频振荡,严重时可以破坏模式的预报。因此,在数值模式积分之前需要把初始条件中不协调和不平衡的噪音过滤掉。四维变分同化虽然引入了预报模式作为同化过程的一个约束条件,理论上,它产生的分析场已具有较好的平衡特性,但由于诸多环节,如内、外循环的粗、细网格间的反复插值,流依赖背景误差协方差的引进,更多的物理过程参数化引进等,仍会导致噪音,造成初始场的动力不平衡。因此,在四维变分同化系统的分析中,仍需要有初始化的处理。

数字滤波是一种控制重力波的有效手段,它通过对模式积分的时间序列作滤波处理,可消除高频部分而保留低频部分,具有容易实施等优点。在变分同化系统中,重力波的控制可以通过在代价函数中增加反映重力波振幅的项来实现。三维变分同化不涉及到模式积分,其初始化是一个跟随在变分分析之后的单纯的模式积分时间序列的处理过程。四维变分同化引入了预报模式作为一项约束,寻求的是相空间中与整个同化时间窗的观测最接近的模式运行的轨迹,轨迹的光滑程度自然影响到与观测的距离,因此,四维变分与初始化不再是两个独立的阶段,而是相互影响,相互渗透。

全球GRAPES-4DVar的新初始化方案采用了数字滤波初始化与四维变分同化弱约束结合的方式控制高频重力波。该方案的优点为:数字滤波初始化的积分时间与四维变分同化时间窗一致,对重力波的处理与极小化过程同步进行,不会产生额外的计算时间;弱约束的数字滤波初始化更适合于长时间窗同化,不会因同化时间窗的延长而削弱慢波;可在同化过程中,通过调节重力波控制项权重系数方便地控制模式轨迹的光滑程度。相对于数字滤波初始化与四维变分同化独立处理的旧初始化方案,新方案的分析和预报效果比旧方案略好或相当,但新初始化方案有效节省了全球GRAPES-4DVar循环同化系统每次分析的运行时间,这对需要大量计算时间的四维变分来说尤为重要。

GRAPES-4DVar的新初始化方案可进一步优化。如考虑使用“总能量范数”来定义重力波控制项的归一化矩阵Q,将使各个模式要素具有相似的光滑程度。在最后一次的高分辨率外循环中再做一次数字滤波,可以减少从低分辨率的分析增量到高分辨率的分析带来的插值振荡。此外,还可尝试水汽的滤波。

参考文献
纪立人. 2011. 数值天气预报发展进程中若干亮点的回顾及其启迪. 气象科技进展, 1(1): 40–43. Ji L R. 2011. Some highlights and their implication in the early progress of numerical weather prediction:A review. Adv Meteor Sci Technol, 1(1): 40–43. (in Chinese)
刘艳, 薛纪善, 张林, 等. 2016. GRAPES全球三维变分同化系统的检验与诊断. 应用气象学报, 27(1): 1–15. Liu Y, Xue J S, Zhang L, et al. 2016. Verification and diagnostics for data assimilation system of global GRAPES. J Appl Meteor Sci, 27(1): 1–15. (in Chinese)
马旭林, 庄照荣, 薛纪善, 等. 2009. GRAPES非静力数值预报模式的三维变分资料同化系统的发展. 气象学报, 67(1): 50–60. Ma X L, Zhuang Z R, Xue J S, et al. 2009. Development of 3-D variational data assimilation system forthe nonhydrostatic numerical weather prediction model-GRAPES. Acta Meteor Sinica, 67(1): 50–60. DOI:10.3321/j.issn:0577-6619.2009.01.006 (in Chinese)
张林, 刘永柱. 2017. GRAPES全球四维变分同化系统极小化算法预调节. 应用气象学报, 28(2): 168–176. Zhang L, Liu Y Z. 2017. The Preconditioning of minimization algorithm in GRAPES global four-dimensional variational data assimilation system. J Appl Meteor Sci, 28(2): 168–176. (in Chinese)
刘永柱, 张林, 金之雁. 2017. GRAPES全球切线性和伴随模式的调优. 应用气象学报, 28(1): 62–71. Liu Y Z, Zhang L, Jin Z Y. 2017. The optimization of GRAPES global tangent linear model and adjoint model. J Appl Meteor Sci, 28(1): 62–71. (in Chinese)
薛纪善, 庄世宇, 朱国富, 等. 2008. GRAPES新一代全球/区域变分同化系统研究. 科学通报, 53(20): 2408–2417. Xue J S, Zhuang S Y, Zhu G F, et al. 2008. Scientific design and preliminary results of three-dimensional variational data assimilation system of GRAPES. Chinese Sci Bull, 53(20): 2408–2417. DOI:10.3321/j.issn:0023-074X.2008.20.003 (in Chinese)
薛纪善, 陈德辉. 2008. 数值预报系统GRAPES的科学设计与应用. 北京: 科学出版社: 383pp. Xue J S, Chen D H. 2008. Scientific Design and Application of GRAPES. Beijing: Science Press: 383pp. (in Chinese)
张华, 薛纪善, 庄世宇, 等. 2004. GRAPES三维变分同化系统的理想试验. 气象学报, 62(1): 31–41. Zhang H, Xue J S, Zhuang S Y, et al. 2004. Idea experiments of GRAPES three-dimensional variational data assimilation system. Acta Meteor Sinica, 62(1): 31–41. DOI:10.3321/j.issn:0577-6619.2004.01.004 (in Chinese)
庄世宇, 薛纪善, 朱国富, 等. 2005. GRAPES全球三维变分同化系统:基本设计方案与理想试验. 大气科学, 29(6): 872–884. Zhuang S Y, Xue J S, Zhu G F, et al. 2005. GRAPES global 3D-Var System:Basic scheme design and single observation test. Chinese J Atmos Sci, 29(6): 872–884. DOI:10.3878/j.issn.1006-9895.2005.06.04 (in Chinese)
Arakawa A, Schubert W H. 1974. Interaction of a cumulus cloud ensemble with the large-scale environment, Part Ⅰ. J Atmos Sci, 31(3): 674–701. DOI:10.1175/1520-0469(1974)031<0674:IOACCE>2.0.CO;2
Ballish B, Cao X H, Kalnay E, et al. 1992. Incremental nonlinear normal-mode initialization. Mon Wea Rev, 120(8): 1723–1734. DOI:10.1175/1520-0493(1992)120<1723:INNMI>2.0.CO;2
Courtier P, Thépaut J N, Hollingsworth A. 1994. A strategy for operational implementation of 4D-Var, using an incremental approach. Quart J Roy Meteor Soc, 120(519): 1367–1387. DOI:10.1002/(ISSN)1477-870X
Daley R. 1997. Atmospheric Data Analysis. New York: Cambridge University Press: 457pp.
Dee D P, Uppala S M, Simmons A J, et al. 2011. The ERA-Interim reanalysis:Configuration and performance of the data assimilation system. Quart J Roy Meteor Soc, 137(656): 553–597. DOI:10.1002/qj.v137.656
Ehrendorfer M, Errico R M, Raeder K D. 1999. Singular-vector perturbation growth in a primitive equation model with moist physics. J Atmos Sci, 56(11): 1627–1648. DOI:10.1175/1520-0469(1999)056<1627:SVPGIA>2.0.CO;2
Errico R M. 2000. Interpretations of the total energy and rotational energy norms applied to determination of singular vectors. Quart J Roy Meteor Soc, 126(566): 1581–1599. DOI:10.1256/smsqj.56602
Gauthier P, Thépaut J N. 2001. Impact of the digital filter as a weak constraint in the preoperational 4DVAR assimilation system of Météo-France. Mon Wea Rev, 129(8): 2089–2102. DOI:10.1175/1520-0493(2001)129<2089:IOTDFA>2.0.CO;2
Hong S Y, Lim J O J. 2006. The WRF single-moment 6-class microphysics scheme (WSM6). J Korean Meteor Soc, 42: 129–151.
Liu K, Chen Q Y, Sun J. 2015. Modification of cumulus convection and planetary boundary layer schemes in the GRAPES global model. J Meteor Res, 29(5): 806–822. DOI:10.1007/s13351-015-5043-5
Liu Y Z, Zhang L, Lian ZH. 2018. Conjugate gradient algorithm in the four-dimensional variational data assimilation system in GRAPES. J Meteor Res, 32(6): 974–984. DOI:10.1007/s13351-018-8053-2
Lynch P, Huang X Y. 1992. Initialization of the HIRLAM model using a digital filter. Mon Wea Rev, 120(6): 1019–1034. DOI:10.1175/1520-0493(1992)120<1019:IOTHMU>2.0.CO;2
Lynch P, Huang X Y. 2010. Initialization//Lahoz W, Khattatov B, Menard R. Data Assimilation: Making Sense of Observations. Berlin, Heidelberg: Springer, 241-259
Mlawer E J, Taubman S J, Brown P D, et al. 1997. Radiative transfer for inhomogeneous atmospheres:RRTM, a validated correlated-k model for the longwave. J Geophys Res, 102(D14): 16663–16682. DOI:10.1029/97JD00237
Polavarapu S, Tanguay M, Fillion L. 2000. Four-dimensional variational data assimilation with digital filter initialization. Mon Wea Rev, 128(7): 2491–2510. DOI:10.1175/1520-0493(2000)128<2491:FDVDAW>2.0.CO;2
Rabier F. 2005. Overview of global data assimilation developments in numerical weather-prediction centres. Quart J Roy Meteor Soc, 131(613): 3215–3233. DOI:10.1256/qj.05.129
Rawlins F, Ballard S P, Bovis K J, et al. 2007. The Met Office global four-dimensional variational data assimilation scheme. Quart J Roy Meteor Soc, 133(623): 347–362. DOI:10.1002/(ISSN)1477-870X
Temperton C. 1988. Implicit normal mode initialization. Mon Wea Rev, 116(5): 1013–1031. DOI:10.1175/1520-0493(1988)116<1013:INMI>2.0.CO;2
Temperton C. 1989. Implicit normal mode initialization for spectral models. Mon Wea Rev, 117(2): 436–451. DOI:10.1175/1520-0493(1989)117<0436:INMIFS>2.0.CO;2
Troen I B, Mahrt L. 1986. A simple model of the atmospheric boundary layer; sensitivity to surface evaporation. Bound-Layer Meteor, 37(1-2): 129–148. DOI:10.1007/BF00122760
Wee T K, Kuo Y H. 2004. Impact of a digital filter as a weak constraint in MM5 4DVAR:An observing system simulation experiment. Mon Wea Rev, 132(2): 543–559. DOI:10.1175/1520-0493(2004)132<0543:IOADFA>2.0.CO;2
Williamson D L, Temperton C. 1981. Normal mode initialization for a multilevel grid-point model. Part Ⅱ:Nonlinear aspects. Mon Wea Rev, 109(4): 744–757. DOI:10.1175/1520-0493(1981)109<0744:NMIFAM>2.0.CO;2