1 引　言

GRAPES全球四维变分资料同化（4DVar）系统已实现业务运行（Zhang，et al，2019）。4DVar使用气候背景误差（简称为背景误差），该方法的基本原理表明，虽然在其同化时间窗内能通过切线性与伴随模式积分来隐式计算流依赖背景误差，但在下一次同化分析时4DVar对流依赖背景误差信息失去“记忆”能力。要克服这一缺陷需要在4DVar中引入集合样本估计的集合预报误差，发展集合四维变分资料同化（En4DVar）方法，从而实现在同化时间窗起始时刻就引入流依赖背景误差。为了构建和运行En4DVar，首要问题是生成适当的集合样本并在同化过程中合理应用。如何生成适合GRAPES同化的集合样本并在GRAPES全球4DVar中合理使用集合预报误差，在GRAPES 全球En4DVar的研发中还没有成熟经验可循。

利用集合卡尔曼滤波（EnKF，ensemble Kalman filtering）类方法或基于4DVar的集合资料同化（EDA）方法生成集合样本，并将其用于变分资料同化，动态估计随天气流型变化的集合预报误差，是国际资料同化研究的重点和业务发展主流（Clayton，et al，2013；Kleist，et al，2015a，2015b；Buehner，et al，2013；Bonavita，et al，2012；Raynaud，et al，2011）。集合卡尔曼滤波（Evensen，1994；Houtekamer，et al，2005）、集合方根滤波（EnSRF，Whitaker，et al，2002，2008；Tippett，et al，2003）、集合卡尔曼变换（ETKF，Bishop，et al，2001）等分别被加拿大气象局、美国环境预报中心、英国气象局采用。以上方法中，EnKF方法是基础性方法，其他方法是EnKF方法的变形或简化（Tippett，et al，2003）。其中，EnKF与集合方根滤波都属于资料同化方法，通过集合样本估计卡尔曼增益矩阵，并通过使用观测资料和观测误差计算分析场。这两种方法的差异在于如何在不同样本中引入观测扰动的影响，EnKF通过一定的误差模型直接扰动观测资料（Burgers，et al，1998），但这种显式扰动会造成抽样误差增大；集合方根滤波方法通过重新定义卡尔曼增益矩阵的表达式来隐式表述观测扰动的影响，但其只能采用顺序同化方式每次使用一个观测资料，无法像EnKF那样同时使用多个观测资料，需要特殊的并行计算处理。集合卡尔曼变换方法不是一种资料同化方法，而是一种集合预报样本产生方法，只使用观测误差而不使用观测值，并不求解卡尔曼增益矩阵。由于集合卡尔曼变换是在观测空间计算变换矩阵，其优点是计算量小、生成集合样本快速；其不足是集合预报误差在观测空间的变化与模式空间地理位置的分布特征不对应。

基于EnKF类方法生成集合样本，可以平行更新每个集合样本的分析场（Evensen，1994），具有平行计算程度高可扩展性好、可以同时生成大量集合样本的优点，但该方法要求在业务变分资料同化系统之外独立发展一套集合卡尔曼滤波系统，开发和维护成本较高。基于4DVar的EDA方法利用扰动观测资料和表述模式不确定性的模式扰动技术来生成集合样本（Bonavita，et al，2012），为欧洲中期天气预报中心（ECMWF）、法国气象局所采用。与基于EnKF类方法的方案相比，EDA直接使用业务已有4DVar系统，只需对观测资料、边界条件、模式物理过程倾向叠加随机扰动并调整运行流程即可以计算分析样本，更加方便业务数值预报系统的维护。而EDA方法的不足在于4DVar的计算量非常大，限制了集合样本数的增加，并且为了减少计算消耗系统一般在较低空间分辨率下运行。目前在世界各主要业务中心的混合同化系统发展中采用接近或超过100个集合样本数并对集合预报误差给予较大的权重是主要趋势，因而提高EDA方法的计算效率非常关键。

Houtekamer等（2005）研究指出，估计集合预报误差时会存在由于集合样本有限造成的抽样误差、观测误差估计不当、正演观测算子存在误差、预报与观测误差出现非高斯分布特征、预报模式中存在未能考虑的误差等问题，这些因素的综合作用会造成集合预报误差低估。为减缓低估产生的不利影响研究者们发展了不同的方差膨胀方法，如分析误差方差膨胀法（Anderson，et al，1999）、误差额外膨胀法（Houtekamer，et al，2005）、分析扰动幅度张驰到预报扰动幅度法（Zhang，et al，2004），以及分析误差幅度张弛到先验估计的集合预报误差幅度方法（Whitaker，et al，2012）等。此外，Wang等（2003）用新息向量来对分析扰动进行膨胀处理。ECMWF在减缓集合预报误差存在低估现象时，一方面在观测扰动基础上，进一步引入针对多种物理过程的随机扰动来表征已有参数化方案的不确定性，另一方面通过比较集合离散度与集合预报误差来定量诊断集合预报误差低估幅度，并将全球划分为南、北半球、赤道3个区域，对每个区域采用不同的膨胀系数（Bonavita，et al，2012）。

中国在集合卡尔曼滤波方法上进行了多项研究探索，其中基于GRAPES区域模式的研究工作有EnKF方法（庄照荣等，2011），以及将ETKF估计的集合预报误差用于区域三维变分资料同化的混合资料同化方法（Chen，et al，2015）。在GRAPES全球模式中也开展了ETKF初值扰动方法研究（马旭林等，2008），利用模拟探空观测资料进行试验，还没有考虑卫星等非局地观测资料的局地化技术。刘柏年等（2015）开展了集合预报误差的方差滤波方法研究，采用谱滤波方法根据信号和噪声尺度的统计特征构造低通滤波器来滤除集合预报方差估计值中的大部分随机取样噪声。

考虑到GRAPES全球模式还没有形成业务可用的EnKF系统，本研究借助现有GRAPES全球4DVar系统，通过扰动观测资料和扰动海温分析场，利用EDA方法来生成集合样本。研究的第一部分（龚建东，2021）讨论了如何有效应用集合预报误差的科学方案，确定了集合预报误差在GRAPES全球4DVar中应用的分析框架。然而要将集合预报误差实际应用于GRAPES全球4DVar中，还需要解决如何高效生成集合样本的问题，以及确定与系统相匹配的同化关键参数。本研究重点是将实际估计的集合预报误差用于GRAPES_GFS全球4DVar中，通过批量数值试验合理给定多项关键参数，并通过数值试验来验证En4DVar相对于4DVar的优势。

2 EDA集合样本生成方案 2.1 集合样本生成方法

EDA方法利用不同扰动观测资料产生多组4DVar分析，并在模式积分时扰动海温分析场来考虑下垫面的不确定性，以及通过随机物理过程表示模式不确定性，最后获得多组集合样本。扰动观测和海温分析场可以表示为分析和预报步骤

$\begin{split}J\left({{\rm{\delta }}{\tilde X}} \right) =& \frac{\,1\,}{\,2\,}{\rm{\delta }}{{{\tilde X}}^{\rm{T}}}{{B}^{ - 1}}{\rm{\delta }}{\tilde X} + \frac{\,1\,}{\,2\,}\sum\limits_{{{i}} = 0}^{{L}} {{{\Big[\left({{{{\tilde d}}_{{i}}} - {{H}_{{i}}}{{M}_{0,{{i}}}}{\rm{\delta }}{\tilde X}} \right)}^{\rm{T}}}}{\text •} \\& {{{R}}^{ - 1}}\left({{{{\tilde d}}_{{i}}} - {{H}_{{i}}}{{M}_{0,{{i}}}}{\rm{\delta }}{\tilde X}} \right)\Big] + {J_{\rm{c}}}\\[-10pt]\end{split}$

(1)

${{{{\tilde X}}^{\rm{b}}} = {{\cal M}_{0,{{L}}}}\left({{{{\tilde X}}^{\rm{b}}} + {\rm{\delta }}{\tilde X}} \right) + {\tilde Q}}$

(2)

式（1）为分析步骤，右端第1项为背景项，第2项为观测项，J_c是基于数字滤波的弱约束项，用来抑制切线性模式时间积分中因分析场动力不平衡产生的噪音。B是背景误差协方差，R是观测误差协方差， ${\cal H}$ 为非线性观测算子。 $\widetilde {\left(\cdot \right)}$ 表示对一组集合样本中的任意一个样本进行分析或预报，形成分析或预报场。 ${\rm{\delta }}{\tilde X}$ 为集合样本中的任意一个样本的分析增量。M_0,i是非线性预报模式 ${\cal M}$ 的切线性模式，下标i表示时间t_i。 ${{\tilde d}_{{i}}} = {{y}^{\rm{o}}} + {{\eta }_{{i}}} - {{\cal H}_{{i}}}{{\cal M}_{{{0,i}}}}\left({{{{\tilde X}}^{\rm{b}}}} \right)$ 为在观测场y^o上叠加观测误差扰动扰动 ${{\eta }_{{i}}}$ 后的新息向量， ${{\eta }_{{i}}}$ 服从期望值为0、方差为R的观测误差分布。式（2）为预报步骤，其中 ${{\tilde Q}}$ 为在模式预报中叠加的海温分析场随机扰动，它服从期望值为0、方差为海温分析误差的随机扰动。EDA方法在GRAPES模式中具体实现时，考虑了海温分析场误差存在水平相关的影响（龚建东等，2020）。式（1）—（2）所表示的集合样本生成的具体内容如表1所示。

上述EDA方法基于GRAPES全球4DVar系统，由于单个4DVar计算消耗较大限制了集合样本数的增加。Courtier等（1994）研究表明，在4DVar中采用增量分析方案并适当降低内循环水平分辨率，可以在基本不影响分析结果的情形下显著减小计算消耗。GRAPES全球4DVar就采用了这样的分析方案，模式轨迹及新息向量计算在高分辨率上进行，而切线性与伴随模式时间积分、目标函数极小化迭代计算在低分辨率上进行。本研究中GRAPES全球4DVar的高、低分辨率配置分别为0.25°与1.0°，由于EDA方法只需提供与全球4DVar内循环分辨率匹配的集合样本，因而其高、低分辨率可以设置为1.0°与1.5°，这显著减少了EDA方法所需的计算量。在此基础上顺序进行EDA分析样本和预报样本的生成计算，所有样本累计所耗时间仍然明显大于单次高分辨率GRAPES全球4DVar的计算时间，通过采用平行作业方法生成集合样本的计算时间还将继续缩短。

2.2 极小化预调节信息

分析EDA方法的计算消耗，其计算量主要集中在目标函数极小化迭代求解过程中切线性与伴随模式的时间积分。张林等（2017）研究表明，间隔相近的两次4DVar分析对应的海森矩阵（目标泛函相对于控制变量的二阶偏导数构成的矩阵）变化不大，前一次极小化迭代过程中产生的信息可提供给下一次极小化做预调节，经预调节后的极小化算法的目标函数收敛效率能得到明显提高。对于同一时刻不同集合样本中的4DVar分析而言，由于它们的背景场和观测场均较为接近，对应的海森矩阵变化也不大。因而当第一个样本的4DVar完成计算后其在极小化迭代过程中产生的信息也可以用于后续其他样本。对2018年12月1—30日逐6 h同化循环的总计120个样本的统计表明，在引入极小化预调节的情形下，目标函数经过20次极小化迭代后的数值与未做预调节情形下40次迭代的结果相当，也即计算效率提高了1倍（图1a）。进一步比较极小化迭代后目标函数梯度的下降比例，未做预调节情形下经40次迭代梯度下降为初始值的7.5%，而使用预调节情形下经20次迭代梯度下降为初始值的19.5%，梯度下降幅度相差近2倍（图1b）。由于EDA方法中的4DVar主要用于集合样本生成，并不涉及单次分析的精度，梯度下降程度上的差异可以接受。近1个月EDA方法的试验表明，使用极小化预调节后目标函数的减小程度和目标函数梯度下降程度逐日变化不大，表现稳定。以上是在EDA方法中顺序进行集合样本生成时使用前次预调节信息的方式。当集合样本采用并发作业生成时，借鉴ECMWF的经验可以使用高分辨率控制预报的预调节信息用于EDA方法中集合样本的计算。

3 集合预报误差应用 3.1 高频噪声滤除与异常扰动抑制

使用集合样本来估计集合预报误差时存在抽样误差，抽样误差主要与集合样本数有关，对N个集合样本其对应的蒙特卡洛抽样误差是 $\dfrac{1}{{\sqrt {{N}} }}$ 。抽样误差在不同空间尺度上表现不同（Bonavita，et al，2011），长波对应的抽样误差小而短波对应的抽样误差大。对较短的波动而言其误差构成主要是无用的噪声，需要滤除。Raynaud等（2009）、Bonavita等（2012）采用基于三角函数的低通滤波器，此处采用同样的滤波函数

$\rho \left({{n}} \right) = {\rm{co}}{{\rm{s}}^2}\left({\frac{{\rm{\,{\text{π}}\, }}}{\,2\,}{\rm{min}}\left({{{n}},{{{N}}_{{\rm{trunc}}}}} \right)/{{{N}}_{{\rm{trunc}}}}} \right)$

(3)

式中，n为球谐谱函数对应的二维总波数， ${{{{N}}_{{\rm{trunc}}}}}$ 为最大谱截断波数。在试验中按照Bonavita等（2012）经验将最大谱截断波数取为106，滤除随机噪声同时保留天气到中尺度主要扰动波数的扰动振幅。垂直方向的噪声直接来源于扰动观测资料对分析的影响。GRAPES全球4DVar中由于位温扰动是通过气压扰动的垂直梯度计算获得，而垂直梯度计算对气压扰动非常敏感，对气压扰动进行噪音滤除能显著减小位温扰动存在“之”字形虚假振荡现象。采用简单的五点平滑方法进行气压扰动的噪声滤除就能收到较好效果。

通过EDA方法获得集合样本来估计集合预报误差时，在扣除集合平均后所生成的集合扰动样本之间的数值差异较大，当样本数较小时直接用这些扰动样本来估计集合预报误差时容易受到数值异常偏大样本的干扰而出现误差高估现象，且这种高估不能通过高频噪声滤除的方式消除，需要对数值异常偏大值进行抑制。GRAPES全球4DVar中，已经有纬向平均的气候背景误差，以此作为参考，当扰动样本数值明显大于气候背景误差时，采用衰减扰动幅度的方式控制数值。表示为

$\left\{ {\begin{array}{*{20}{l}} {v = v}&{{}\;\!\!v {\text{≤}} {{a}}{{\rm{\sigma }}_{\rm{c}}}}\\ {v = \dfrac{v}{{\left| v \right|}}\left({{{a}}{{\rm{\sigma }}_{\rm{c}}} + \left({{{b}} - {{a}}} \right){{\rm{\sigma }}_{\rm{c}}}{\rm{\gamma }}} \right)}& {{v{\text{＞}}{a}{{\rm{\sigma }}_{\rm{c}}}}} \\ {{{\gamma = a \tan}}\left({\dfrac{{\rm{\,{\text{π}}\, }}}{\,2\,}\dfrac{{{\rm{min}}\left({\left| v \right|/{{\rm{\sigma }}_{\rm{c}}} - {{a}},{{b}}} \right)}}{{{b}}}} \right)}&{} \end{array}} \right.$

(4)

式中，v为任一预报变量， ${{{\rm{\sigma }}_{\rm{c}}}}$ 为背景误差，系数a为扰动开始衰减的阈值系数，b为衰减后的最大值阈值系数。式（4）限制了预报扰动的幅度，当扰动幅度大于 ${{{a}}{{\rm{\sigma }}_{\rm{c}}}}$ 时，开始采用三角函数进行衰减，最大扰动被限制在 ${{{b}}{{\rm{\sigma }}_{\rm{c}}}}$ 以下。试验中风与气压扰动取 ${{a}} = 3{\text{，}} {{b =}} {\rm{4.6}}$ ，对比湿 ${{a}} = 1.8{\text{，}}{{b = 3}}{\rm{.0}}$ ，能获得较好的衰减效果。

3.2 集合预报误差与背景误差

利用2018年6月数据，比较纬向平均集合预报误差与背景误差。在南、北半球对流层1000—100 hPa（对应模式1—48层）风场集合预报误差存在低估现象，普遍低50%—60%，赤道地区集合预报误差低估情况要小于南、北半球，低估30%—40%（图2a）。同样气压集合预报误差也存在低估，其低估程度比风场要大，南、北半球低估60%—70%，赤道地区低估50%—60%（图2b）。值得注意的是，背景误差自身存在一定的偏差，并且这种偏差特征在赤道以及南、北半球表现并不一致。赤道地区的大气混沌特征明显，模式预报误差不易增长，用同一检验时刻不同预报时效预报的差异估计背景误差时容易在赤道地区造成背景误差低估（龚建东等，2006），而在中高纬度地区容易出现背景误差高估。纬向平均的分布特征表明，相对于背景误差，集合预报误差明显小，需要进行膨胀处理。考虑到背景误差在赤道地区低估、南北半球存在高估的问题，不能仅参照背景误差来膨胀集合预报误差。

4 En4DVar参数确定 4.1 权重系数与相关尺度

集合预报误差主要用于对流层大气。在对流层顶附近（约100 hPa，模式48层），随模式层次增加使用衰减函数将集合预报误差权重系数（ ${{\rm{\beta }}_{\rm{e}}}$ ）的数值衰减为0，背景误差权重系数（ ${{\rm{\beta }}_{\rm{c}}}$ ）相应增大，两者之和为1。图3a给出 ${{\rm{\beta }}_{\rm{e}}} = 0.7$ 时权重系数廓线的垂直分布情况。对湿度分析批量试验结果表明，引入集合预报湿度误差对改进4DVar没有明显效果，其原因和现有4DVar中湿度观测数据量偏少（Zhang，et al，2018）、集合预报误差偏小有关。在后续试验中暂不考虑集合预报湿度误差。

水平局地化的相关尺度需要考虑4DVar中背景误差水平相关尺度随高度的变化。GRAPES全球4DVar中背景误差水平相关尺度随高度的变化如图3b所示，对流层不同分析变量之间差异显著，其中比湿的相关尺度最小（约175 km），流函数约600 km，非平衡势函数最大（700—1200 km），非平衡无量纲气压在500 km左右。由于GRAPES全球4DVar的背景误差中流函数是主导控制变量，其他变量的平衡部分由其计算得到（Zhang，et al，2018）。因此，在选择水平局地化相关尺度（L_loc）时，参照流函数背景误差水平相关尺度（ ${L_{\rm{\psi }}}$ ）的垂直分布特征来给出。批量试验的结果表明，当 ${L_{{\rm{loc}}}} = 1.4{L_{\rm{\psi }}}$ 时（图3b），也即在对流层 ${L_{{\rm{loc}}}}$ 取值在850 km左右时，能获得较好的改进结果。

4.2 集合样本数影响

当集合样本足够大时能减少抽样误差影响，可以合理地估计集合预报误差协方差。在增加集合样本数的方法中，Xu等（2008a，2008b）提出在分析时刻前后利用集合样本在时间序列上扩展来增加集合预报数，可以明显改进区域EnKF的分析与预报效果。Huang等（2018）提出利用时间错位扰动（VTSP，valid time shifting perturbation method）方法，通过增加样本在时间上的输出频次，由时间错位来增加集合样本数。具体做法是，对需要集合样本的时间T，使用由同一分析场起报、分别预报到 $T - \tau{\text{、}}T{\text{、}}T + \tau$ 这3个时次的集合样本，其中 $\tau $ 为时间错位的间隔。为了避免不同时间的大尺度场存在的相位误差，采用扣除对应自身时间集合平均的扰动样本，并采用 $T - \tau $ 与 $T + \tau $ 对称采样方式。研究表明，这种方法对台风等系统扰动样本基本不存在相位误差。在GRAPES全球EDA方法生成集合样本时，没有使用同一起报时间不同预报时长的样本，而是采用当前起报时间的 $T{\text{、}}T + \tau$ 集合样本，和6 h之前起报、时间对应 $T - \tau $ 时刻的集合样本。由式（1）可知，同化时间窗起始和终止时间的6 h，是有观测约束的时段，代表了分析场。本研究所使用集合样本，所有样本的预报时长都大于6 h，代表了预报场。采用以上方法对每个集合预报通过增加3 h预报时效的较小计算代价，实现集合样本增加到3倍，最终产生60个集合样本。

对EDA方法产生的20个集合样本，以及通过VTSP方法进行样本扩展后产生的60个集合样本，比较其引入4DVar后对分析与预报的影响。取背景误差权重系数 ${{\rm{\beta }}_{\rm{c}}} = 0.8$ ，集合预报误差权重系数 ${{\rm{\beta }}_{\rm{e}}} = 0.2$ ，集合预报误差未做膨胀处理，结果表明20个集合样本引入后，GRAPES全球En4DVar与4DVar的分析与预报效果相当，集合预报误差引入4DVar后的作用不明显。当集合样本增加到60个后，对南、北半球第1天的风、温度与位势高度预报效果均有改善，更长时效的改进不明显，对赤道地区的改善略大一些，特别是风场和位势高度场的改进较明显，且正影响可以达到5 d，En4DVar的效果略优于4DVar（图略）。其他多组对照试验表明，只引入20个集合样本，通过集合预报方差膨胀、集合预报误差权重系数加大，均不能产生明显改进效果，表明集合样本数对GRAPES全球En4DVar起到关键作用。

4.3 集合预报方差膨胀

集合预报误差与背景误差的比较（图2）表明，在利用集合样本估计集合预报误差时出现低估问题，且南、北半球比赤道地区的低估更加严重。Bonavita等（2012）在比较ECMWF EDA法生成集合样本的离散度和均方根误差时，发现集合样本的离散度明显低于均方根误差，且在南、北半球要高于赤道地区。一种基于集合样本离散度和均方根误差的膨胀方法被发展出来，利用前5天样本序列的滑动平均来计算方差膨胀因子，北半球1和8月的平均膨胀幅度在1.5左右。由于EDA法生成的每个集合样本相当于单个4DVar分析预报循环，并不依赖于集合预报方差，集合预报方差主要用于En4DVar。在本研究中采用了静态膨胀处理法，具体做法是将EDA法生成的1个月集合样本的离散度与预先估计的气候误差进行比较来计算膨胀因子。预先估计的气候误差是将GRAPES全球4DVar背景误差均方根误差和全球集合预报计算的离散度进行加权平均得到，其中两者加权系数分别是0.6和0.4。计算得到的集合预报误差膨胀系数如图4所示。对风场南北半球的膨胀为1.3—1.4倍，而赤道地区基本没有膨胀。在垂直层次上风场膨胀主要在48层以下的对流层，48层以上系数逐步减小。结合图3a，在平流层权重主要放在背景误差上，集合预报误差的作用减少。对于无量纲气压膨胀发生在所有层次，其中48层以下膨胀幅度基本在1.6左右，在48层以上膨胀系数逐步减小，南北半球与赤道的差异不明显。可以发现风场与无量纲气压的膨胀特点在数值和分布上呈现出明显的不同。

取背景误差权重系数 ${{\rm{\beta }}_{\rm{c}}} = 0.2$ ，集合预报误差权重系数 ${{\rm{\beta }}_{\rm{e}}} = 0.8$ ，采用60个集合样本，比较集合预报误差未做方差膨胀和进行方差膨胀的影响。在进行方差膨胀时分别考虑风速与无量纲气压变量均膨胀1.7倍，以及按照图4所示的系数进行膨胀，并进行1个月的批量试验（图略）。结果表明集合预报误差膨胀1.7倍对南、北半球有正影响，而在赤道地区为负影响，En4DVar不如4DVar。而对风速与无量纲气压合理膨胀后北半球不同等压面层次的风场、位势高度场与温度场均有明显改进，在南半球也有所改进但不及北半球显著。赤道地区的改善不明显。这表明对不同分析变量采用符合各自误差特点的膨胀系数要优于对分析变量采用均值膨胀系数。后续试验将采用图4所示的系数进行膨胀。

4.4 误差权重系数影响

合理给定背景误差与集合预报误差权重系数对改善En4DVar的分析与预报质量至关重要。在全球混合资料同化系统发展中，世界上不同数值预报中心采用不同的权重系数。分析各业务数值预报中心的背景误差与集合预报误差的权重系数，对20个左右的集合样本，其抽样误差为20%，集合预报误差权重系数在0.2左右，背景误差在0.8左右（Clayton，et al，2013），而当集合样本数达到80个时，抽样误差减小到11%，集合预报误差权重系数达到0.75（Kleist，et al，2015a，2015b）。文中使用60个集合样本，其中20个通过EDA方法生成，其余40个集合样本通过VTSP方法生成，并不是完全独立的集合样本，需要进行不同权重系数的试验来寻找最佳参数。对集合预报误差权重系数 ${{\rm{\beta }}_{\rm{e}}}$ 分别取权重0.5、0.7和0.9，对2018年6月10日至7月9日，以及2018年11月21日—2019年1月3日，进行夏、冬两个代表时段的数值试验。由图5、6所示北半球夏季结果可见，随着 ${{\rm{\beta }}_{\rm{e}}}$ 取值增大，温度场在100 hPa以上层次改进增大，但在对流层250 hPa附近，当 ${{\rm{\beta }}_{\rm{e}}} = 0.9$ 时预报误差明显增大；同时，对流层风场第1天预报的误差明显减小，但在96 h以后的预报时效，当 ${{\rm{\beta }}_{\rm{e}}} = 0.7$ 时的效果最好。赤道和南半球也有类似结果（图略），表明权重系数取值0.7较为合理。对于北半球冬季试验结果类似。为此，文中后续试验中设定 ${{\rm{\beta }}_{\rm{e}}} = 0.7$ ，这与Kleist等（2015a，2015b）选定的权重系数接近。

5 试验效果分析

通过数值试验确定集合预报误差与背景误差权重系数、水平局地化相关尺度、风速与无量纲气压的集合预报误差膨胀系数后，还需要通过批量试验验证其对分析和预报的影响。分别进行北半球夏、冬两季各52 d的En4DVar相对于4DVar的批量试验，重点关注对南、北半球及赤道地区的影响。

5.1 夏季预报试验

夏季是北半球大范围强降雨、极端高温、强对流等高影响天气易发且影响范围较广的时段。文中北半球夏季试验的时段设定为2018年6月10日至7月31日。试验中GRAPES全球模式采用0.25°水平分辨率，垂直60层，模式顶高36354 m（4—5 hPa）。相应的4DVar外循环水平分辨率为0.25°，内循环为1.0°，下降算法采用带预调节的L-BFGS（Limited memory Broyden Fletcher Goldfarb Shanno）方法。4DVar中使用的观测资料包括探空、地面、船舶、飞机等常规观测资料，卫星微波温度探测资料，红外高光谱资料，全球导航系统掩星折射率资料，卫星散射计测风资料，以及极轨与静止气象卫星云导风资料（Zhang，et al，2019）。这些资料在4DVar中按对应时间分配在同化时间窗的30 min间隔的时间槽内。En4DVar的系统配置与4DVar相同，背景误差与集合预报误差的权重系数分别为 ${{\rm{\beta }}_{\rm{c}}} = 0.3$ 与 ${{\rm{\beta }}_{\rm{e}}} = 0.7$ ，采用通过VTSP方法将EDA方法产生的20个集合样本扩展为60个集合样本。

5.1.1 对分析的影响

首先分析集合预报误差引入4DVar后对分析场的影响，取ECMWF再分析场作为“实况”，夏季52 d的统计结果如图7。南、北半球，集合预报误差引入4DVar后，从对流层700 hPa一直到平流层50 hPa的位势高度场均有改善，且南半球更加显著，主要表现在均方根误差减小、偏差略有减小。赤道地区，En4DVar对风场，尤其是经向风场改善明显。

En4DVar相对于4DVar在赤道地区的改善也可以从均方根误差的时间序列上表现出来，以赤道地区300 hPa为例，位势高度场、温度与经向风场误差的演变情况如图8所示。可以看到在7月20日之前，En4DVar的误差虽然略小于4DVar，但差异并不大，而在7月20日之后的时段内4DVar误差出现了突然增大的现象，而En4DVar则抑制了这一现象。有必要对此做进一步分析，以理解En4DVar表现出正效果的内在原因。

选择7月22日00时，分析集合预报误差的特征及对分析造成的影响。60个集合样本估计的集合预报误差在大西洋赤道东部地区及非洲大陆中南部，出现风场集合预报误差大值区。此外，在太平洋中东部和南美洲北部地区也存在较大的风场集合预报误差（图9a）。当集合预报误差被引入4DVar后对应的分析增量出现明显的变化（图9b），在大西洋赤道东部地区及非洲大陆中南部出现分析值差异的大值区，在太平洋中东部地区也出现明显的分析差异。集合预报误差和分析值差异存在较好的对应关系，这是En4DVar分析误差和6 h预报误差减小的内在原因。

5.1.2 对北半球夏季预报的影响

En4DVar相对于4DVar的预报，对北半球24 h预报有明显改进，48 h有改善，更长时效的预报技巧与4DVar相当；对南半球48 h预报有明显改进，72 h有改善；对赤道地区72 h预报有明显改善，120 h有改善。北半球预报改进的时效短与北半球夏季的气候背景有关。夏季中纬度天气尺度到小尺度天气系统较为活跃，集合预报误差的数值偏小，造成集合预报误差的影响较小。在后文中进行了冬季预报试验，来比较季节变化对北半球的影响。从预报效果来看，改进最显著的是赤道地区。图10给出的是赤道地区250 hPa温度与风场的预报误差，可以看到En4DVar对应的误差更小，通过显著性检验的有改进时效可以达到7 d。

5.2 北半球冬季预报试验 5.2.1 对分析影响

北半球冬季试验的配置与夏季一致，具体时段为2018年11月20日至2019年1月22日，其中预报时段是2018年12月1日至2019年1月22日。冬季试验主要侧重讨论北半球在季节变换时集合预报误差的影响。首先分析集合预报误差在冷、暖季节的变化，结果如图11所示。对于风场误差而言，南、北半球冷暖季的数值相当，冷季误差略小一些。赤道地区低层风场误差大值区出现随季节南北移动的现象，这与赤道对流辐合带的季节变动相对应。此外，北半球冷季风场的集合预报误差在边界层增大明显。对于无量纲气压冷、暖两季的变化特征明显，北半球冷季的集合预报误差明显大于暖季，而南半球冷季的集合预报误差要明显小于暖季，这与风场存在明显的不同。而在赤道地区无量纲气压的误差特征同样体现出对流辐合带随季节南北移动的影响。

北半球冬季52 d的统计结果和夏季预报试验接近，主要差异是对北半球位势高度的均方根误差，En4DVar较4DVar的改进主要集中在700—30 hPa，较夏季试验结果更加深厚，且在200—50 hPa的改进幅度较大（图12a），而夏季试验主要改进在400—150 hPa。北半球位势高度的均方根误差变化与集合预报无量纲气压误差有明显的对应关系。北半球冬季试验时南半球处于暖季，En4DVar主要改进集中在700—150 hPa，而对应的南半球冷季如图7b所示，改进主要位于700—30 hPa。这可能与冷、暖季节主导的天气系统不同有关，冷季以大尺度锋面天气系统为主，天气系统比较深厚，而暖季以中小尺度天气系统为主，其影响主要是强降雨、强对流天气等过程的高度。对于北半球风场的改进，同样是冷季（图12b）大于暖季（图12c），且改进的垂直高度更加深厚。而对于南半球风场冷、暖季的改进相近。

以300 hPa位势高度的均方根误差为例来讨论集合预报误差对分析的影响。2018年12月25—30日东亚发生了一次寒潮天气过程，影响范围广，中国中东部地区普遍降温6—10℃。从形势场来看北半球被5波大尺度系统控制（图14b）。从分析误差来看位势高度的分析误差与背景误差在寒潮天气过程期间，En4DVar要明显小于4DVar（图13），对风场和温度场表现不明显。检查25—30日所使用的观测资料的数量并不能发现明显的差异。改进的原因与对照试验4DVar采用背景误差，而En4DVar可以使用集合预报误差，在资料相同的情况下可以更好地使用观测资料有关。由图14a可见，在北极冷涡附近、亚洲等与形势场槽线对应的位置，气压的集合预报误差明显增大，表明背景场的不确定性较大，这一预报误差在En4DVar中使用后，与集合预报误差大值区对应的位置，如在北极冷涡附近En4DVar相对于4DVar的分析增量差异明显增大，其他区域的分析增量也有明显变化（图14b），相应的分析误差也有所减小（图略）。这表明集合预报误差能反映天气形势变化的特点，可以改善分析质量。

5.2.2 对冬季预报的影响

北半球冬季预报试验效果要明显好于夏季预报试验效果，特别是对北半球。从距平相关系数和均方根误差上看，En4DVar相对于4DVar在北半球48 h预报有明显改进，风场改进可以达到72 h，南半球72 h也有明显改进。改进最显著的是赤道地区，前120 h预报评分明显提高（图略），这些改进在700—250 hPa整个层次表现非常一致。

6 总结与讨论

在GRAPES全球业务4DVar系统中引入集合样本估计的集合预报误差，发展En4DVar同化方法，进一步改进系统的分析和预报质量。第一部分（龚建东，2021）重点讨论了集合预报误差引入方式、水平与垂直局地化对分析平衡性的影响，以及En4DVar的计算效率。第二部分重点讨论集合预报误差与背景误差在GRAPES全球4DVar中合理使用问题，并通过批量试验验证改进效果。

考虑到En4DVar系统包括4DVar和集合样本生成两个系统，研发面向业务应用，集合样本生成的计算效率是非常重要的一个方面。集合样本生成方法上选择基于GRAPES全球4DVar的EDA方法。EDA方法生成的每个集合样本都需要通过4DVar的分析和预报来实现，计算量较大。除了采用低分辨率轨迹计算和更低分辨率的内循环计算外，本研究重复使用第一个样本极小化迭代过程中产生的预调节信息，并用于其他样本极小化做预调节，使得20次极小化迭代后目标函数的减小程度和40次相当，将计算效率提高了1倍。此外，借鉴Xu等（2008a，2008b）的连续时间采样方法和Huang等（2018）时间错位扰动方法，通过增加样本在时间上的输出频次，由时间错位来增加集合样本数，在对每个集合预报通过增加3 h预报的少许计算代价，实现集合样本额外增加2倍，最终产生60个集合样本。比较EDA方法产生的20个集合样本，以及通过VTSP方法进行样本扩展后产生的60个集合样本，多组对照试验表明，20个集合样本引入后，通过集合预报方差膨胀、集合预报误差权重系数加大，均不能产生明显改进效果，而当集合样本增加到60个后，对分析与预报有改进效果。

在此基础上，通过高频噪声滤除与异常扰动抑制，使得集合样本估计的集合预报误差分布特征更加合理。最后采用了静态膨胀处理方法，将EDA方法生成的1个月集合样本的离散度与预先估计的气候误差进行比较，计算膨胀因子，由该膨胀因子对每个集合样本进行放大。处理后的集合样本进入到4DVar中，用于估计集合预报误差。文中主要考虑在4DVar中应用对流层大气的集合预报误差。在对流层顶附近（约100 hPa的位置，模式48层），随模式层次增加使用衰减函数将集合预报误差权重系数衰减为0。通过批量试验，选择水平局地化相关尺度（L_loc）为流函数背景误差水平相关（ ${L_{\rm{\psi }}}$ ）的1.4倍，在对流层约为850 km。此外，通过批量数值试验方法来确定背景误差与集合预报误差的权重系数，当 ${{\rm{\beta }}_{\rm{e}}} = 0.7$ 时的预报效果最好。

文中进一步通过北半球夏、冬季各52 d的批量试验，验证En4DVar相对于4DVar的改进效果。结果表明，北半球En4DVar对位势高度误差的改进冷季更加深厚，暖季试验主要改进在400—150 hPa，而冷季的改进主要集中在700—30 hPa，且在200—50 hPa的改进幅度较大。北半球位势高度的均方根误差变化与集合预报无量纲气压误差有明显的对应关系。南半球En4DVar在暖季的改进主要集中在700—150 hPa，在冷季改进主要位于700—30 hPa。赤道地区受季节影响较小，En4DVar对位势高度、风场与温度的改进都较为明显，其中最显著的是经向风场。预报试验表明，引入En4DVar后南、北半球2—3 d的预报有明显改进，赤道地区5 d有明显改进。

En4DVar还有一些方面需要继续优化。首先，集合样本估计的湿度预报误差还没有在4DVar中使用。GRAPES全球模式在赤道边界层顶附近存在较大的湿度偏差，围绕湿度偏差的改进工作正在开展，后续将考虑使用湿度集合预报误差。再者，目前集合预报误差采用静态膨胀处理方法，区分了不同变量、不同高度和纬度的膨胀系数，但还不能实现随季节变化的动态变化。ECMWF发展了一种基于集合样本离散度和均方根误差的膨胀方法，利用前5 d样本序列的滑动平均来计算方差膨胀因子（Bonavita，et al，2011），后续工作将考虑动态的膨胀处理技术，更好的利用集合预报误差。

此外，将集合预报误差用于4DVar目前还没有针对台风进行专门的工作，如台风涡旋位置和强度扰动、台风涡旋重定位等。在夏季试验期间的台风路径预报没有发现改进效果，针对台风还需要进行研究工作。

最后，使用基于4DVar的EDA方法生成集合样本，其最大优势是实现比较简单，但计算量较大，样本需要逐一生成。从业务应用上看，发展基于集合的同化方法，如四维集合变分同化（4DEnVar，four dimensional ensemble variational data assimilation）方法、DRP-4DVar（dimension reduced projection 4DVar）方法（Wang，et al，2010），来并行生成大量的集合样本是未来发展的一个趋势，这方面的探索工作已经在开展中。

致　谢：感谢中国气象局数值预报中心陆慧娟博士、刘艳博士、王金成博士在全球四维变分同化基础版本上所进行的基础性工作，以及支持完成本研究提供的帮助。