气象学报  2021, Vol. 79 Issue (3): 443-457   PDF    
http://dx.doi.org/10.11676/qxxb2021.021
中国气象学会主办。
0

文章信息

庄照荣, 李兴良, 刘艳, 陈春刚. 2021.
ZHUANG Zhaorong, LI Xingliang, LIU Yan, CHEN Chungang. 2021.
数字滤波初始化方案在高分辨率GRAPES区域模式中的应用研究
A study on digital filter initialization in high-resolution GRAPES regional model
气象学报, 79(3): 443-457.
Acta Meteorologica Sinica, 79(3): 443-457.
http://dx.doi.org/10.11676/qxxb2021.021

文章历史

2020-06-16 收稿
2021-01-07 改回
数字滤波初始化方案在高分辨率GRAPES区域模式中的应用研究
庄照荣1,2 , 李兴良1,2 , 刘艳1,2 , 陈春刚3     
1. 国家气象中心,北京,100081;
2. 中国气象局数值预报中心,北京,100081;
3. 西安交通大学,西安,710049
摘要: 高分辨率同化系统中引入的高频率高密度观测信息会造成分析变量之间的不平衡,因而在模式积分中产生的虚假重力波会严重影响模式预报质量和运行稳定。为了抑制虚假重力波在模式中快速增长,文中通过在高分辨率GRAPES区域模式中应用非绝热数字滤波和增量数字滤波初始化方案,研究不同滤波截断周期的初始化方案对分析预报的影响。试验结果表明,不同初始化方案的分析预报效果依赖于截断周期,截断周期越长滤除的高频信息越多,对于3 km分辨率模式,30 min与15 min截断周期的初始化方案既滤除了高频噪音,也基本保留了有意义的分析信息。不同截断周期的数字滤波方案在GRAPES-MESO系统中对气象要素预报场的影响基本可忽略;而在GRAPES-RAFS系统中,初始化方案对分析和预报的影响会逐渐累积,15 min截断周期的初始化方案比2 h截断周期的初始化方案的风场预报和降水预报质量更优。
关键词: 非绝热数字滤波    增量数字滤波    初始化方案    GRAPES区域模式    
A study on digital filter initialization in high-resolution GRAPES regional model
ZHUANG Zhaorong1,2 , LI Xingliang1,2 , LIU Yan1,2 , CHEN Chungang3     
1. National Meteorological Centre,Beijing 100081,China;
2. Numerical Weather Prediction Center of CMA,Beijing 100081,China;
3. Xi'an Jiaotong University,Xi'an 710049,China
Abstract: The imbalance of the analysis variables created by introducing irregular high frequency and high density observations into the high-resolution analysis gives rise to unphysical gravity-inertia waves in the model integration, and thus affects the forecast quality and stability of operational models. In this paper, the diabatic digital-filtering initialization and incremental digital filter initialization are applied to 3 km resolution GRAPES regional model to suppress noises in the early forecast hours, and the impact on analysis and forecast with different cutoff periods of initializations are evaluated. The analysis and forecast experiments show that the filtering effect relies on the cutoff period. If the cutoff period is longer, more high frequency signals are removed. In the 3 km resolution modelling, both 30 min and 15 min cutoff periods initialization can remove high frequency noises and retain meaningful analysis information. The impact produced by different cutoff period initialization schemes in the GRAPES-MESO system on the model variables can be neglected, whereas the accumulated effects of initialization schemes are gradually observed in the GRAPES-RAFS system. The results show that wind and precipitation forecast qualities are obviously improved with a 15 min cutoff period initialization compared to that with the 2 h cutoff period. The impact of digital filter initialization scheme on high-resolution model prediction depends on the cutoff period. The cutoff period of 30 min or 15 min is suitable for 3 km resolution model, and the initialization scheme of incremental digital filter demonstrates obvious advantages.
Key words: Diabatic digital filtering    Incremental digital filtering    Initialization scheme    GRAPES regional model    
1 引 言

高分辨率同化系统不仅同化常规观测资料,还同化高频率高密度雷达和卫星等非常规观测资料。高分辨率分析场作为数值模式的初值对千米尺度的模式预报质量非常重要。分析场中引入非规则分布的观测信息后,分析要素之间不满足准平衡关系,而分析变量之间的不平衡噪音在模式积分过程中会激发出没有气象意义的高频振荡重力波,严重时会影响模式预报质量和稳定性。为了抑制虚假重力波在模式中快速增长,在分析要素进入模式预报系统之前需要进行初始化过程。

重力波的控制可以通过调整初值中的重力波分量来实现,方法主要有牛顿松弛(Bao,et al,1997)、非线性正规模(Daley,1979Williamson,et al,1981)、数字滤波和增量分析更新初始化等。目前牛顿松弛和非线性正规模初始化方案在数值预报模式中很少应用,而数字滤波和增量分析更新方案因有很好的初始化效果,在许多数值预报模式中被广泛使用。最初的数字滤波为绝热数字滤波(ADFI:Adiabatic Digital Filter Initialization)初始化方案,该方案分别向前和向后模式积分,而且都关闭物理过程(Lynch,et al,1992),对于当时6 min时间步长的模式,采用3和6 h的截断周期都可以滤除高频噪音,并获得相似的结果,但12 h的截断周期会阻尼气象模态的发展。虽然绝热数字滤波初始化增量改动较小,但是由于采用的绝热模式和后来模式预报使用的非绝热模式不一致,造成该初始化方案不能控制非绝热模式产生的噪音,因而又提出了非绝热数字滤波方案(DDFI:Diabatic Digital Filter Initialization)(Huang,et al,19931994Lynch,et al,1997)。Huang等(1993)研究指出,非绝热滤波方案采用向前非绝热积分的时间序列场进行滤波,其预报产生的噪音程度比绝热数字滤波方案低,同时考虑到较长的截断周期影响计算时效,通过最优滤波设计对于6 min时间步长可采用最小滤波周期如2—3 h的截断周期就可以达到可接受的初始化结果。由于数字滤波方案简单有效,在许多数值预报中心的四维变分同化系统中也将数字滤波作为目标函数的约束项进行初始化过程(Polavarapu,et al,2000Gauthier,et al,2001Wee,et al,2004王舒畅等,2011刘艳等,2019)。由于非绝热数字滤波初始化增量比较大,并且模式绝热向后积分过程中地面变量是保持不变的,但在非绝热模式向前积分时,随着积分步长变化的高空变量和保持不变的地面变量并不匹配。进而又提出了增量数字滤波方案(IDFI:Incremental Digital Filter Initialization)。该方案只有很小的初始化增量,不仅能减少动力热力调整时间,还能保持背景场快速增长的模态(Lynch,et al,1994)。但是增量数字滤波的计算量比非绝热数字滤波翻倍,而且若背景场有噪音,对模式预报也有隐患。增量分析更新初始化方案(IAU:Incremental Analysis Updating)通过在一段时间里把瞬时分析增量平滑地插值到预报模式系统来缓解分析场引入模式产生的噪音(Bloom,et al,1996Polavarapu,et al,2004)。在快速更新分析预报循环系统中增量分析更新初始化方案提高了湿度场和动力场的平衡,因而能加快极小化过程和减少模式预报的动力、热力调整时间(Lee,et al,2006)。由于数字滤波的初始化效果较好并容易实施,因此,全球/区域分析和模式预报系统(GRAPES:Global/Regional Assimilation and Prediction System)中都采用非绝热数字滤波方案(薛纪善等,20082012)。在GRAPES-3DVar分析与预报系统中,资料同化和初始化是两个独立实施的过程,会增加计算量。当同化系统升级到全球四维变分同化时进一步发展为弱约束的数字滤波初始化方法,约束强加在分析增量上与极小化迭代过程同步进行,数字滤波初始化方案是变分同化系统的一部分,不会对四维变分产生额外的计算资源消耗(刘艳等,2019)。

目前GRAPES-RAFS(RAFS:Rapid Analysis and Forecast System)分析预报循环系统分辨率从10 km升级为3 km,原数字滤波方案是否适用?非绝热数字滤波方案是否能滤除高频噪音而又不影响预报质量等问题都有待解决。有研究表明,在高水平分辨率模式系统中有天气意义的波动与高频噪音的界限更加模糊,数字滤波初始化能对特定切断频率的高频振荡进行滤波,但无法区分滤除的信号是具有天气意义的高频信号还是虚假的高频噪音。而要消除积分过程中的高频振荡,必须延长滤波时间窗,但这是以增加额外积分时间以及有可能损害气象上有意义的波动为代价的(陈敏等,2012)。

文中研究不同截断频率的非绝热数字滤波和增量数字滤波初始化方案在3 km分辨率GRAPES区域模式中的应用情况,分析不同数字滤波方案对初始场和预报场的滤波效果,以及通过批量GRAPES-MESO和GRAPES-RAFS试验考察不同初始化方案对气象要素预报质量的影响。

2 数字滤波方案 2.1 非绝热数字滤波

按照数字滤波的理论(Lynch,et al,1992Huang,et al,1993),通过对短时间内的模式预报场进行时间滤波可以滤除高频部分噪音,在GRAPES区域模式系统中的数字滤波方案采用非绝热数字滤波(薛纪善等,20082012),即

$ {f_0}^* = \sum\limits_{k = - N}^N {{H_k}{f_k}} $ (1)

若非绝热数字滤波的截断周期为T $\Delta t$ 为时间步长,则在周期内时刻序列 $k$ 取值为 ${\rm{ - }}N$ $N$ $N = T/(2\Delta t)$ ${f_k}$ 为截断周期 $ - T/2$ $T/2$ 时段,离散时刻 $k$ 的计算值, ${f_0}^*$ 为滤波后的初始场。非绝热数字滤波方案先把分析场通过模式绝热向后积分到 $ - T/2$ 处,然后在 $ - T/2$ 处模式非绝热向前积分到 $T/2$ 处,此时 ${f_k}$ 为从 $ - T/2$ $T/2$ 时段内模式非绝热向前积分的预报场,如图1所示。 ${H_k}$ 为数字滤波系数,即

图 1  非绝热数字滤波初始化方案示意 Fig. 1  Schematic diagram of diabatic digital filtering initialization scheme
$ {H_k} = \frac{{\sin (k{\theta _{\rm c}})}}{{k{\text π} }}\frac{{\sin [k{\text π} /(N + 1)]}}{{k{\text π} /(N + 1)}} $ (2)

式中,定义 ${h_k} = \dfrac{{\sin (k{\theta _{\rm c}})}}{{k{\text π} }}$ ,其中 ${\theta _{\rm c}}$ 为截断频率,在数字滤波中通过 ${h_k}$ 的作用会把频率高于 ${\theta _{\rm c}}$ 的高频信息滤除, ${w_k} = \dfrac{{\sin [k{\text π} /(N + 1)]}}{{k{\text π} /(N + 1)}}$ 为窗函数,由于傅里叶系数的截断会产生吉布斯振荡,因而窗函数的作用是滤除截断函数在高频部分的不稳定噪音以减少吉布斯振荡(Lynch,et al,1992)。若 $\omega $ 为截断圆频率, $\Delta t$ 为时间步长,则有

${\theta _{\rm c}} = \omega \Delta t = \frac{{2{\text π} \Delta t}}{{{T_{\rm c}}}}$ (3)

式中, ${T_{\rm c}}$ 为截断周期。从式(3)可以看出,当模式的时间步长确定,截断频率随着截断周期而变化,

${T_{\rm c}} = 2{\text π} /\omega $ (4)

从以上公式可知,非绝热数字滤波方案为通过对一段时间的模式积分预报场进行时间滤波,滤波效果与截断周期有关。定义转换函数 $T(\theta)$ 为一次滤波过程,若 ${f_k}{\rm{ = }}{{\rm e}^{{\rm i}k\theta }}$ ${f_0}^* = T(\theta) \cdot {f_0}$ ,其中i为虚数单位,θ为频率,f0k=0时的气象场,由于对称系数Hk=Hk,代入式(1)可得(Lynch,et al,1992),

$T(\theta) = \sum\limits_{k = - N}^N {{H_k}} {{\rm e}^{ - {\rm i}k\theta }} = \left({H_0} + 2\sum\limits_{k = 1}^N {{H_k}} \cos (k\theta)\right)$ (5)

$T(\theta)$ 函数可以看出滤波过程中不同截断频率下高频信息滤除的程度。图2为模式时间步长30 s情况下不同截断周期下的转换函数变化情况,可以看出理想截断频率为0.1。当选用不同的截断周期,转换函数分布大不相同。截断周期为2 h时,0.05以上频率的信息滤除过多;截断周期为30 min时,频率0.25以上的波动都被滤除,截断周期为15 min,频率0.45以上的信息大部分被滤除。也就是说截断周期越短,保留的高频信息越多。从图2中也可以看出,当不使用窗函数时,同样截断周期的转换函数在高频部分会有不稳定波动,采用窗函数后,不稳定波动噪音也被滤除。

图 2  不同截断周期下的转换函数 Fig. 2  Transfer functions for the filter with different cutoff periods Tc

数字滤波为时间滤波,滤波系数 ${H_k}$ 相当于在滤波区间内对各预报场所取的权重,即滤波区间内各预报场对初始场的影响程度。时间步长为30 s的水平分辨率3 km模式中,滤波系数与截断周期的关系如图3所示。可以看出,当截断周期越长,滤波系数在初始时刻的权重越小;截断周期越短,滤波权重在初始时刻越大。意味着截断周期越长,分析时刻模式预报场占的比重越小,通过滤波区间内所有模式场的影响,滤除了分析场中可能引起虚假重力波的噪音。

图 3  不同截断周期下的滤波系数 Fig. 3  Filtering coefficients with different cutoff periods Tc
2.2 增量数字滤波

在非绝热数字滤波方案中,对分析全量滤除了噪音,但实际上背景场是模式预报场,变量之间通常是协调的,一般假定背景场没有噪音,因而只对加入不规则观测资料信息后的分析增量进行滤波。文中采用增量数字滤波方案,即

${x_{{\rm{ini}}}} = {x_{\rm{b}}} + (x_{\rm{a}}^{{\rm{df}}} - x_{\rm{b}}^{{\rm{df}}})$ (6)

式中, ${x_{\rm{b}}}$ 为背景场, $x_{\rm{b}}^{{\rm{df}}}$ 为对背景场进行非绝热数字滤波, $x_{\rm{a}}^{{\rm{df}}}$ 为对分析场进行非绝热数字滤波,增量数字滤波的初始场 ${x_{{\rm{ini}}}}$ 是在原背景场上叠加滤波后的分析增量。增量数字滤波的优点为可以保留背景场中快速增长的模态,也可以加速动力、热力调整过程,获得更小的初始化增量。但由式(6)可见,增量数字滤波需要分别对背景场和分析场进行非绝热数字滤波,因而滤波过程比非绝热数字滤波初始化方案增加了1倍计算量。

3 GRAPES-MESO数值预报试验 3.1 试验设置

为了测试数字滤波方案的影响,文中对GRAPES区域模式系统进行5组数字滤波试验,试验时段为2018年6月12日00时—28日12时(世界时,下同),试验范围为中国东部(17º—50ºN,102º—135ºE),模式分辨率为3 km,垂直为51层,模式层顶高33 km。试验在GRAPES-3DVar分析和云分析的基础上进行不同数字滤波初始化方案,预报时效为24 h,GRAPES区域模式时间步长30 s,各组试验采用相同的GRAPES区域模式物理过程参数设置。不同数字滤波初始化方案的试验设置见表1

表 1  试验设置 Table 1  Experiments design
试验名称 试验说明 数字滤波参数
Ana 不采用初始化方案(No DFI)
T2h 非绝热数字滤波(DDFI)        ${T_{\rm c}}=2\;{\rm h}{\text{,}}{\theta _{\rm c}}{\rm{ = }}\dfrac{{\text π} }{{{\rm{1}}2{\rm{0}}}}=0.026$
T30 非绝热数字滤波(DDFI)        ${T_{\rm c}}=30\;{\rm{min}}{\text{,}}{\theta _{\rm c}}{\rm{ = }}\dfrac{{\text π} }{{{\rm{30}}}}=0.104$
T30inct 在T30基础上进行增量数字滤波(IDFI)        ${T_{\rm c}}=30\;{\rm{min}}{\text{,}}{\theta _{\rm c}}{\rm{ = }}\dfrac{{\text π} }{{{\rm{30}}}}=0.104$
T15 非绝热数字滤波(DDFI)        ${T_{\rm c}}=15\;{\rm{min}}{\text{,}}{\theta _{\rm c}}{\rm{ = }}\dfrac{{\text π} }{{{\rm{15}}}}=0.209$
T15inct 在T15基础上进行增量数字滤波(IDFI)        ${T_{\rm c}}=15\;{\rm{min}}{\text{,}}{\theta _{\rm c}}{\rm{ = }}\dfrac{{\text π} }{{{\rm{15}}}}=0.209$

Ana为对照试验,直接采用分析结果进行预报,不做数字滤波初始化处理。T2h试验按照业务水平10 km分辨率 GRAPES-RAFS的数字滤波设置,滤波截断周期为2 h。T30/T15试验的数字滤波截断周期分别为30 min和15 min,T30inct/T15inct试验为在T30/T15试验基础上进行的增量数字滤波,保留了背景场中所有高频信息。

数字滤波初始化方案主要用来抑制初始时刻动力场和质量场不协调激发出来的高频振荡噪音(Lynch,et al,1992薛纪善等,2008),因此文中仅对模式层上的无量纲气压、位温、UVW风场进行数字滤波初始化过程,未调整湿度场,与区域GRAPES-3DVar和全球GRAPES-4DVar业务中设置相同(刘艳等,2019)。

3.2 滤波效果分析

文中在批量试验中选取2018年6月15日00时的初始场和预报场进行分析,研究采用数字滤波初始化方案后对初始场和预报场的滤波效果。

3.2.1 初始场影响

从2018年6月15日00时模式面第24层的初始场U风场(图4)可以看出,经过三维变分分析后,U分析风场有更多细微信息(图4a),当滤波截断周期从15 min(图4d)、30 min(图4c)到2 h(图4b),U风场逐渐滤除较多的中小尺度信息,T2h试验的初始U风场已经非常平滑。如果采用增量数字滤波,T30inct与T15inct的初始U风场都非常接近Ana试验,T15inct与Ana的初始风场差别最小。

图 4  2018年6月15日00时第24层初始U风场 (a. 无数字滤波,b. 2 h数字滤波,c. 30 min数字滤波,d. 15 min数字滤波,e. 30 min增量数字滤波,f. 15 min增量数字滤波;单位:m/s) Fig. 4  Initial fields of U-component at the 24th model level at 00:00 UTC 15 June 2018 for different experiments (a. No DFI,b. DDFI with Tc=2 h,c. DDFI with Tc=30 min,d. DDFI with Tc=15 min,e. IDFI with Tc=30 min,f. IDFI with Tc=15 min;unit:m/s)

初始化方案目的是在尽可能小的改动分析场基础上,使初始场更协调,减小分析噪音。因而初始化方案的改动要相对小于分析增量场,不影响分析质量。表2统计2018年6月15日00时不同数字滤波方案的初始化增量(初始场减去分析场),以及Ana试验的分析增量(分析场减去背景场)全场平均的均方根误差(RMSE)以及绝对最大值(MAX)情况。

表 2  2018年6月15日00时不同试验的分析增量与初始化增量 Table 2  Analysis increments and initialization increments at 00:00 UTC 15 June 2018
试验名称 U(m/s) V(m/s) 位温(K) 无量纲气压
RMSE MAX RMSE MAX RMSE MAX RMSE MAX
Ana 1.33 −8.28 1.44 −8.75 0.72 5.56 1.46×10−4 5.9×10−4
T2h 0.67 19.65 0.72 −20.27 0.56 −13.49 1.28×10−4 1.77×10−3
T30 0.26 10.59 0.27 −11.37 0.22 −8.67 6.08×10−5 1.2×10−3
T15 0.15 9.53 0.16 8.52 0.14 −7.0 4.24×10−5 1.0×10−3
T30inct 0.12 −3.94 0.11 −6.08 0.09 −2.36 3.69×10−5 −3.2×10−4
T15inct 0.07 4.00 0.06 5.19 0.07 1.79 1.96×10−5 −1.8×10−4

表2可以看出,数字滤波方案各个变量的初始化增量的均方根误差都比分析增量的均方根误差小,而且滤波周期越短,对分析场的改动越小,初始化增量的均方根误差越小。对于同样的滤波周期,增量非绝热数字滤波方案明显比全量非绝热数字滤波方案对分析场改动的少。5组数字滤波试验中的T15inct初始化增量的均方根误差最小,初始场最接近Ana试验,与图4的结论一致。由于没有对湿度变量进行数字滤波初始化,湿度初始场即为分析场。T2h方案采用的滤波周期太长,最大初始化增量比分析增量大得多,从图4也可以看出T2h方案对分析场改动过大,初始场太平滑;T30与T15数字滤波方案的最大初始化增量与最大分析增量基本相当;T30inct与T5inct方案的最大初始化增量明显小于最大分析增量。

3.2.2 功能谱分析

上一节主要从时空域分析不同初始化方案中初始场的差别,本节从频谱空间研究不同数字滤波初始化方案对初始场的影响。文中采用二维离散余弦转换(2D-DCT:Discrete Cosine Transform)方法对有限区域二维气象场进行谱分解(郑永骏等,2008庄照荣等,20182020),比较不同数字滤波方案的初始场(xinit)与分析场(xa)的功率谱差别(即初始化增量功率谱)。为了便于比较,文中对初始化增量功率谱的绝对值取对数(即:lg|xinitxa|)。从2018年6月15日00时模式面第24层不同初始化方案的初始化增量功率谱(图5)可以看出,对于无量纲气压变量,由于其本身随时间的变率很小,不同初始化方案的初始化增量功率谱差别很小。对于位温变量,不同初始化方案的初始化增量功率谱在339 km以下部分差别比较明显,其中T2h试验的初始化增量功率谱最大,其次是T30试验,T15试验和两组增量数字滤波方案的初始化增量功率谱差别不大。对于风场变量,不同初始化方案的初始化增量功率谱差别很大。T2h试验的初始场和分析场的功率谱相差最大,并且不仅在中小尺度部分有较大差别,在大尺度部分差别也非常明显。说明T2h试验的初始化方案滤除了有气象意义的波动信息。T30试验的初始化增量功率谱进一步缩小,T15试验的初始化增量功率谱又进一步下降。同样,两组增量数字滤波方案的风场初始化增量功率谱相当,并且在几组试验中最小,说明增量数字滤波方案能保留更多的中小尺度信息。从图5还可以看出,在不同波段上,不同初始化方案的初始化增量功率谱差别也有所不同,例如T2h试验在所有波段上的风场增量功率谱和其他初始化方案的差别都很大,T30、T15试验及同一截断周期的增量数字滤波试验的风场初始化增量功率谱的差别主要在641 km以下尺度范围内。这也说明T2h试验的初始风场与分析场差别在大尺度上都比较明显,T30、T15试验的初始风场与分析场的差别主要在中小尺度信息上,与图4结果一致。

图 5  2018年6月15日00时第24层分析场与数字滤波后的初始场功率谱差别 (a. 无量纲气压,b. 位温,c. U风场,d. V风场;单位:m3/s2 Fig. 5  Power spectra of analysis minus initial field at the 24th model level at 00:00 UTC 15 June 2018 (a. non-dimensional pressure,b. potential temperature,c. U-component,d. V-component;unit:m3/s2
3.2.3 预报场影响

分别选取海上和平原两点来考察采用不同数字滤波方案后预报的地面气压变化情况。从3 h地面气压预报(图6)可以看出,数字滤波周期越长,信息滤除越多,地面气压随时间演变越平滑。而增量数字滤波保留了背景场中所有信息,因而地面气压变化略比同样滤波周期的非绝热数字滤波波动大,其中周期为2 h的数字滤波方案使地面气压变化过于平滑。比较海上和平原两点的地面气压变化也可以看出,在预报前1 h内平原上地面气压变化幅度明显大于海上,随后平原上地面气压随时间积分增加,其变化幅度有所减缓。

图 6  前3 h预报的地面气压变化 (单位:hPa;a. 海上 (30ºN,128ºE),b. 平原 (30ºN,108ºE)) Fig. 6  Surface pressure variations (unit:hPa) in the first 3 h at two model grids (a. ocean point (30ºN,128ºE),b. land point (30ºN,108ºE))

从全场平均的地面气压倾向(图7)可以看出,初始时刻的地面气压倾向最大,随着预报步数的增加,变量间进行调整,气压倾向逐渐减小。在前3 h预报中,滤波周期为15 min(T15试验)的地面气压倾向略低于对照试验,滤波周期30 min的地面气压倾向比T15试验进一步减小,但预报3 h后,Ana、T15、T30的地面气压倾向的差别越来越小。当滤波周期为2 h时,初始地面气压倾向就很小,约0.14 hPa/(10 min),随着模式积分,前3 h的地面气压倾向略微减小且趋于稳定,说明T2h试验虽然消除了分析噪音,但只保留了大尺度部分信息,因而变量之间很协调,气压倾向非常小。T15inct,T30inct的增量数字滤波方案的地面气压倾向几乎与同样滤波周期的非绝热数字滤波方案重合,可以达到非绝热数字滤波滤除噪音的效果,同时也保留背景场中协调的气象场信息。

图 7  2018年6月15日00时起报的地面气压倾向 (单位:hPa/(10 min)) Fig. 7  Surface pressure tendency starting from 00:00 UTC 15 June 2018 (unit:hPa/(10 min))
3.3 批量试验结果

为了研究不同数字滤波初始化方案对预报质量的影响,本节将不同初始化方案的数值预报试验结果与观测要素进行比较,考察初始化方案对高分辨率天气预报的影响。

3.3.1 气温及风预报场

比较17 d地面气象要素24 h预报与观测的平均偏差(图8),对于前12 h预报,不同初始化方案的2 m气温预报略有不同,T2h试验与Ana试验无初始化方案的预报偏差差别最大,达0.1 K;T30和T15试验与Ana试验预报偏差差别小于0.05 K;两个增量数字滤波方案的2 m气温偏差与Ana试验几乎相同,意味着增量数字滤波方案相对于同样滤波周期的全量数字滤波方案对预报场影响更小;对于10 m风场,不同初始化方案与Ana试验的预报偏差几乎没有区别。对于12 h后的预报,不同初始化方案对地面气象要素预报几乎没有影响。

图 8  

无论对2 m气温还是10 m风场,不同初始化方案的预报均方根误差也几乎一样。对于高空气象要素场12 h预报(图略),几组初始化方案对要素预报的全场平均预报偏差和均方根误差也几乎一样,表明初始化方案对12 h高空气温和风预报场的影响可以忽略。

好的数字滤波方案的应用应使分析场改变很小且在可接受范围内,并能使预报质量不下降(Lynch,et al,1992)。由上述分析可知,文中不同数字滤波初始化方案对气温和风预报的影响很小,只有T2h试验的2 m气温3 h预报偏差与Ana试验差别略大。

3.3.2 降水预报

这里采用常用的公平风险评分(ETS,Equitable Threat Score)来客观评估降水预报准确度。比较6组试验各时段平均累计降水预报ETS评分(图9)可以看出,不同初始化方案的降水预报质量略有不同,但差别并不显著。相对于Ana试验,T2h试验0—6 h除暴雨量级外ETS评分都略低于Ana试验;6—12 h不同量级降水的ETS评分都明显低于Ana试验;12—18 h大雨ETS评分显著低于Ana试验,中雨和暴雨量级降水预报优于Ana试验;18—24 h从小雨到大雨ETS评分都低于Ana试验。总之,T2h试验滤波截断周期过长,滤除了一些有气象意义的信息造成24 h内整体降水预报质量有所下降。

图 9  不同时段降水预报ETS (a. 0—6 h,b. 6—12 h,c. 12—18 h,d. 18—24 h) Fig. 9  ETS of accumulated precipitation forecast at various forecast lead times (a. 0—6 h,b. 6—12 h,c. 12—18 h,d. 18—24 h)

T30试验的降水ETS评分0—6 h大雨和暴雨时明显高于Ana试验,12—18 h除去大雨的其他量级降水都略高于Ana试验;其他时段和量级降水都与Ana试验评分相当或略低。当在滤波周期30 min基础上采用增量滤波时,相对于T30试验,6—24 h中雨降水ETS评分有所提高,6—18 h大雨降水ETS评分也有所提高。

T15试验的降水ETS评分0—6 h、12—18 h大雨以上量级和Ana试验的差别略大,其他时段和量级降水ETS评分接近Ana试验结果。T15inct试验在各个时段和各量级的降水ETS评分与T30inct基本相当,二者的降水预报ETS评分相对最接近Ana试验。

4 GRAPES-RAFS数值预报试验

在快速分析预报循环系统中,多次分析会造成初始场变量之间更加不协调,严重影响模式的稳定性,而初始化方案在循环中多次应用在降低分析噪音的同时,对模式预报质量也会有影响。本节在GRAPES-RAFS分析预报循环中应用不同的数字滤波方案检验初始化方案对预报质量的影响。

4.1 试验设置

采用GRAPES-RAFS进行一天两次(00/12时起始)每3 h循环的12 h间歇分析预报试验,试验时段为2018年6月11—29日。每个循环过程开始时采用全球模式降尺度的预报作为背景场,循环过程中采用区域模式自身预报作为背景场,文中每次三维变分分析后都采用云分析,在分析场进入模式预报系统前进行数字滤波初始化,模式预报24 h。试验范围为中国东部(17º—50ºN,102º—135ºE),模式分辨率为3 km,垂直为51层,模式层顶33 km。进行了两组数字滤波初始化试验,第一组采用2 h截断周期(即T2h),第二组采用15 min截断周期(即T15),进行分析预报循环试验。为了测试不同初始化方案对快速分析预报循环试验的影响,本节主要比较分析预报循环12 h后,分析00时两组数字滤波初始化方案对分析场和预报场的影响。

4.2 气温、风及相对湿度检验

文中比较经过快速分析预报循环12 h后,试验时段00时的背景场、分析场与探空观测的17 d全场平均偏差与标准差(图10),可以看出,采用不同数字滤波方案,在经过多次分析和模式预报后,T15试验的背景风场偏差比T2h试验更接近0线,同时T15试验的背景风场标准差从低层到高层也明显小于T2h试验,分析风场结果类似(图10ab)。不同数字滤波方案温度场(图10c)的背景场和分析场偏差相差均不大。T15试验的背景场温度标准差略小于T2h试验。不同数字滤波方案湿度场(图10d)的背景场和分析场偏差基本没有区别,标准差也相差不大。

图 10  00时分析场/背景场与探空观测的偏差 (a1—d1) 和标准差 (a2—d2) (a. U风场,单位:m/s;b. V风场,单位:m/s;c. 温度,单位:K;d. 相对湿度,单位:%) Fig. 10  Biases (a1—d1) and standard deviations (a2—d2) of analyses /background fields with radiosonde observations at 00:00 UTC (a. U-component,unit:m/s;b. V-component,unit:m/s;c. temperature,unit:K;d. relative humidity,unit:%)
图 10   Fig. 10  Continued

综上可知,在快速分析预报循环中,经过多次数字滤波初始化,采用15 min截断周期的数字滤波方案比2 h截断周期方案对风场分析场和风场的3 h预报场(背景场)有明显正贡献。

经过12 h预报后,与T2h试验相比,T15试验的初始化方案对预报风场还有明显正贡献(11),预报风场与探空观测的偏差从低层到高层比T2h试验更接近0线,风场标准差也略小于T2h试验。但12 h预报的温度场、湿度场的偏差与标准差在两组试验中差别并不明显。

图 11  12 h预报场与探空观测的偏差 (a1—d1) 和标准差 (a2—d2) (a. U风场,单位:m/s;b. V风场,单位:m/s;c. 温度,单位:K;d. 相对湿度,单位:%) Fig. 11  Biases (a1—d1) and standard deviations (a2—d2) of 12 h forecasts with radiosonde observations (a. U-component,unit:m/s;b. V-component,unit:m/s;c. temperature,unit:K;d. relative humidity,unit:%)
图 11   Fig. 11  Continued
4.3 降水预报

2018年6月12—28日12 h循环后00时起报的T15和T2h试验的平均每6 h预报降水ETS如图12所示。可以看出,T15试验的0—6 h降水预报ETS除中雨量级外其他各个降水量级都比T2h试验显著提高,但T15试验的6—12 h降水预报ETS在中雨和暴雨之间低于T2h试验。12—18 h,T15试验在小雨、中雨和大雨量级的ETS显著高于T2h试验;18—24 h,T15试验的降水预报ETS在小雨和暴雨量级高于T2h试验。整体来看,T15试验比T2h试验对24 h内降水预报有优势。

图 12  不同预报时段降水ETS (a. 0—6 h,b. 6—12 h,c. 12—18 h,d. 18—24 h) Fig. 12  ETS of accumulated precipitation forecast at various forecast lead times (a. 0—6 h,b. 6—12 h,c. 12—18 h,d. 18—24 h)
5 结论和讨论

对3 km分辨率的GRAPES区域模式采用非绝热数字滤波和增量数字滤波初始化方案抑制虚假重力波在模式中的增长,研究不同截断周期的数字滤波方案对分析和预报的影响。通过滤波效果分析及数值预报批量试验,得到以下结论:

(1)对水平3 km分辨率,模式时间步长30 s的GRAPES区域模式,2 h截断周期的初始场过于平滑,滤除了有气象意义的信息。30 min与15 min截断周期的初始场既滤除了高频噪音,也基本保留了有用的分析信息。而增量数字滤波方案可以保留背景场中快速增长的模态,只滤除分析增量不协调的噪音。几组方案中T15inct试验的初始化增量最小,初始场最接近分析场。从功率谱分析也可以看出不同数字滤波方案对初始动力场影响最大,两组增量数字滤波方案保留了更多的动力中小尺度信息,与分析场的功率谱差别最小。

(2)GRAPES-MESO的批量数值预报试验表明,不同数字滤波方案对12 h内预报的地面2 m温度场影响略有不同,几组初始化方案中初始场与分析场的平均偏差最大可达0.1 K,而对地面10 m风场预报影响不大。不同数字滤波方案对高空气象要素12 h预报场的影响也不大。不同数字滤波方案对降水预报的影响略有不同,T2h试验整体降低24 h内降水预报的ETS,T30inct与T15inct试验的降水预报ETS相对更接近Ana试验。

(3)GRAPES-RAFS的批量数值预报试验表明,经过多次分析预报循环后,初始化方案对分析和预报的影响会逐渐累积,因而不同初始化方案的分析预报质量有明显的区别。T15试验对风场分析和12 h预报比T2h试验有明显正贡献,同时T15试验24 h内各时段的降水预报质量也整体优于T2h试验。

对于高分辨率模式,由于高频噪音与具有天气意义的高频信号很难区分,增量数字滤波能避免破坏初猜场的气象信息,因而在计算效率和资源有保障的情况下,增量数字滤波初始化方案是更优选择。另外,在能消除积分过程中的高频振荡基础上,较小滤波截断周期可以节约一定的计算时间,也能最大限度地保留有气象意义的高频信息。因而需要在高分辨率模式中选择合适的数字滤波参数来保证中小尺度的预报质量。

文中只对动力场和质量场进行数字滤波,而对于湿度变量以及相关的水成物变量的噪音没有进行处理,云分析后水成物变量是通过松弛逼近(nudging)方式吸收到模式中(朱立娟等,2017)。虽然初始化方案主要用来抑制初始动力场和质量场不平衡而激发的虚假重力波,但一般在数字滤波方案中湿度变量也进行了初始化过程(Lynch,et al,1992),而且在间歇循环中对云水的初始化过程能缓解云水和其他变量的不平衡,同时缩短云和降水产生的调整适应时间(Huang,1996)。因而下一步工作将对GRAPES-RAFS中的水成物变量和其他变量的协调性进行研究。由于非绝热数字滤波需要绝热向后积分,然后再向前积分,在高频同化及高分辨率系统中会占用计算时间,而增量数字滤波方案不仅对分析场还需要对背景场进行非绝热数字滤波,因而在实施过程中又增加了计算时间。此外,在快速分析预报循环中,若背景场的预报时效不长,变量之间难以达到平衡,背景场也可能存在一些噪音,因而增量数字滤波初始化方案是否能有效滤除噪音而保持模式稳定运行还需要进一步研究。同时,在高分辨率模式下,很难区分高频噪音和有气象意义的小尺度波动,因而采用的数字滤波方案也有可能影响中小尺度系统的预报质量。今后针对千米尺度的高频快速分析预报循环系统,也需要寻求其他能满足业务运行时效要求的初始化方案,例如增量分析更新初始化方案。

参考文献
陈敏, Huang X Y, Wang W. 2012. 数字滤波初始化(DFI)在高水平分辨率模式中的应用. 气象学报, 70(1): 109-118. Chen M, Huang X Y, Wang W. 2012. Difficulties in the implementation of the digital filter initialization in a high resolution numerical weather forecast model. Acta Meteor Sinica, 70(1): 109-118. (in Chinese)
刘艳, 薛纪善. 2019. GRAPES的新初始化方案. 气象学报, 77(2): 165-179. Liu Y, Xue J S. 2019. The new initialization scheme of the GRAPES. Acta Meteor Sinica, 77(2): 165-179. (in Chinese)
王舒畅, 李毅, 张卫民等. 2011. 资料同化中的数字滤波弱约束试验及分析. 物理学报, 60(9): 825-832. Wang S C, Li Y, Zhang W M, et al. 2011. Tests and analysis of digital filter weak constrain in data assimilation. Acta Phys Sinica, 60(9): 825-832. (in Chinese)
薛纪善, 陈德辉. 2008. 数值预报系统GRAPES的科学设计与应用. 北京: 科学出版社, 40-41. Xue J S, Chen D H. 2008. Scientific Design and Application of GRAPES. Beijing: Science Press, 40-41 (in Chinese)
薛纪善, 刘艳, 张林等. 2012. GRAPES全球三维变分同化系统模式变量分析版文档. 中国气象局数值预报中心技术手册. 北京: 中国气象局, 92-96. Xue J S, Liu Y, Zhang L, et al. 2012. Scientific Documentation of GRAPES-3DVar Version for Global Model. Numerical Weather Prediction Center. Beijing: China Meteorological Administration, 92-96 (in Chinese)
郑永骏, 金之雁, 陈德辉. 2008. 半隐式半拉格朗日动力框架的动能谱分析. 气象学报, 66(2): 143-157. Zheng Y J, Jin Z Y, Chen D H. 2008. Kinetic energy spectrum analysis in a semi-implicit semi-Lagrangian dynamical framework. Acta Meteor Sinica, 66(2): 143-157. (in Chinese)
朱立娟, 龚建东, 黄丽萍等. 2017. GRAPES三维云初始场形成及在短临预报中的应用. 应用气象学报, 28(1): 38-51. Zhu L J, Gong J D, Huang L P, et al. 2017. Three-dimensional cloud initial field created and applied to GRAPES numerical weather prediction nowcasting. J Appl Meteor Sci, 28(1): 38-51. DOI:10.11898/1001-7313.20170104 (in Chinese)
庄照荣, 陈静, 黄丽萍等. 2018. 全球和区域分析的混合方案对区域预报的影响试验. 气象, 44(12): 1509-1517. Zhuang Z R, Chen J, Huang L P, et al. 2018. Impact experiments for regional forecast using blending method of global and regional analyses. Meteor Mon, 44(12): 1509-1517. DOI:10.7519/j.issn.10000526.2018.12.001 (in Chinese)
庄照荣, 王瑞春, 李兴良. 2020. 全球大尺度信息在3 km GRAPES-RAFS系统中的应用. 气象学报, 78(1): 33-47. Zhuang Z R, Wang R C, Li X L. 2020. Application of global large scale information to GRAEPS RAFS system. Acta Meteor Sinica, 78(1): 33-47. (in Chinese)
Bao J W, Errico R M. 1997. An adjoint examination of a nudging method for data assimilation. Mon Wea Rev, 125(6): 1355-1373. DOI:10.1175/1520-0493(1997)125<1355:AAEOAN>2.0.CO;2
Bloom S C, Takacs L L, da Silva A M, et al. 1996. Data assimilation using incremental analysis updates. Mon Wea Rev, 124(6): 1256-1271. DOI:10.1175/1520-0493(1996)124<1256:DAUIAU>2.0.CO;2
Daley R. 1979. The application of non-linear normal mode initialization to an operational forecast model. Atmos Ocean, 17(2): 97-124. DOI:10.1080/07055900.1979.9649054
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
Huang X Y, Lynch P. 1993. Diabatic digital-filtering initialization: Application to the HIRLAM Model. Mon Wea Rev, 121(2): 589-603. DOI:10.1175/1520-0493(1993)121<0589:DDFIAT>2.0.CO;2
Huang X Y, Cederskov A, Källén E. 1994. A comparison between digital filtering initialization and nonlinear normal-mode initialization in a data assimilation system. Mon Wea Rev, 122(5): 1001-1015. DOI:10.1175/1520-0493(1994)122<1001:ACBDFI>2.0.CO;2
Huang X Y. 1996. Initialization of cloud water content in a data assimilation system. Mon Wea Rev, 124(3): 478-486. DOI:10.1175/1520-0493(1996)124<0478:IOCWCI>2.0.CO;2
Lee M S, Kuo Y H, Barker D M, et al. 2006. Incremental analysis updates initialization technique applied to 10 km MM5 and MM5 3DVAR. Mon Wea Rev, 134(5): 1389-1404. DOI:10.1175/MWR3129.1
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. 1994. Diabatic initialization using recursive filters. Tellus, 46A: 583-597.
Lynch P, Giard D, Ivanovici V. 1997. Improving the efficiency of a digital filtering scheme for diabatic initialization. Mon Wea Rev, 125(8): 1976-1982. DOI:10.1175/1520-0493(1997)125<1976:ITEOAD>2.0.CO;2
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
Polavarapu S, Ren S Z, Clayton A M, et al. 2004. On the relationship between incremental analysis updating and incremental digital filtering. Mon Wea Rev, 132(10): 2495-2502. DOI:10.1175/1520-0493(2004)132<2495:OTRBIA>2.0.CO;2
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 II: Nonlinear aspects. Mon Wea Rev, 109(4): 744-757. DOI:10.1175/1520-0493(1981)109<0744:NMIFAM>2.0.CO;2