1 引言

通过资料同化的方法改善数值模式初始场的精度是提高数值预报准确性的主要手段之一。目前，三维变分方法^[1-5]是应用较广泛的资料同化方法，但其静态的、各向同性的背景场误差协方差，不符合真实大气“流依赖”的特征。张旭斌等^[6-7]针对三维变分方法的这一缺陷，将时间依赖信息及多尺度信息引入背景场误差协方差，从而有效提取观测资料的多尺度信息，对降水预报改善明显。陈耀登等^[8]构造了具有纬度依赖特征的背景场误差协方差，使模式的动力平衡关系更为合理。四维变分方法^[9-12]虽然可以隐性提供“流依赖”的背景场误差协方差，但其求取切线性模式和伴随模式的过程较为复杂。集合卡尔曼滤波同化方法中背景场误差协方差不用直接表达，且具有依天气形势而变的流依赖特征^[13-16]。但是集合卡尔曼滤波法存在着有限的集合样本较难合理估计背景场误差、矩阵不满秩，以及分析场中的变量难以平衡等问题。为结合变分方法和集合方法的优点，同时回避各自的不足，研究人员提出和采用了将变分方案（3DVar或4DVar）与集合方案结合的混合（Hybrid）同化方案^[17-23]。混合同化方案改善了变分方案背景场误差方差模型化、无法依天气形势而变的问题，也一定程度上缓解了集合方案不满秩、变量不协调问题。

然而，集合-变分混合同化方法中如果采用的集合成员数量过少，依然会存在采样误差、流依赖特征不一定合理等问题，为了获得足够数量的集合样本，集合预报过程同样需要消耗较大的计算量^{[17, 24]}。为缓解集合-变分混合同化方案中集合预报所需大量计算资源的问题，杨毅等^[25]尝试将时间扩展取样的预报样本引入集合卡尔曼滤波同化系统，效果好于样本数较少的控制试验。Wang等^[26]将时间滞后的预报场应用于集合-变分混合同化系统，通过不同起报时间获得同一时刻的预报场，构建了流依赖的时间滞后预报样本，从而避免了实际集合预报带来的计算问题，然而由于较长预报时效的时间滞后预报样本的预报偏离等问题，时间滞后集合有效样本数量依然较为有限。Chen等^[27]尝试利用NMC（National Meteorological Center）^[28]方法，将不同预报时效的历史样本引入集合-变分混合同化系统中，节省计算量的同时获得了大量样本，构建了“各向异性”的背景场误差协方差，然而该背景场误差协方差具有较为明显的气候态特征，并不能很好地表示“依流而变”的真实大气状态。

为充分利用时间滞后预报样本流依赖的特性，和历史预报样本数量大的优点，本文尝试在集合-变分混合同化中，利用优选的历史预报样本来扩展时间滞后预报样本，既能够避免时间滞后预报样本数不足的问题，又能够通过优选历史预报样本引入更多的流依赖样本信息，同时继续保持计算资源消耗较少的优势。

2 集合-变分混合同化方法及样本组合 2.1 集合变分混合同化方法

对于集合-变分混合同化方法，其代价函数定义如下^[18-19]：

$ J\left( {\delta {{\bf{x}}_1}, \alpha } \right) = {\beta _1}\frac{1}{2}\delta {\bf{x}}_1^{\rm{T}}{{\bf B}^{ - 1}}\delta {{\bf{x}}_{\rm{1}}}{\beta _2} + \frac{1}{2}{\alpha ^{\rm{T}}}{{\bf{{\rm A}}}^{ - 1}}\alpha + \frac{1}{2}{\left( {{\bf{d}} - {\bf H}\delta {\bf{x}}} \right)^{\rm{T}}}{{\bf{R}}^{ - 1}}\left( {{\bf{d}} - {\bf H}\delta {\bf{x}}} \right) $

(1)

在代价函数中，右边第一项表示为静态背景场误差协方差相关项，其中B为静态的背景场误差协方差矩阵。第二项表示集合误差协方差相关项，α为集合扩展控制变量，A为空间协方差矩阵。β₁和β₂分别为前两项的权重系数，且满足关系式$ \frac{1}{{{\beta _1}}} + \frac{1}{{{\beta _2}}}{\rm{ = }}1 $。第三项为观测项，H表示线性观测算子，R表示观测误差协方差。式中，分析增量δx被定义为：

$ \delta {\bf{x = }}\delta {{\bf{x}}_1} + \sum\limits_{k = 1}^N {\left( {{\alpha _k} \circ {\bf{x}}_k^e} \right)} $

(2)

式中，δx₁表示静态的分析增量，$ {\bf{x}}_k^e $表示集合预报误差的无偏估计，○表示舒尔积（Schur product）。

2.2 样本组合方法

如前所述，为保持计算资源消耗较少的优势，又能充分利用时间滞后预报样本流依赖的特性和历史预报样本数量大的优点，本文尝试在集合-变分混合同化中，利用优选的历史预报样本来扩展时间滞后预报样本。两类样本组合后引入集合-变分混合同化系统的流程如图 1所示。时间滞后样本，历史预报样本的构建及组合方法具体如下。

（1）时间滞后样本。

时间滞后样本由预报过程中不同起报时间、同一时刻的预报场构成，因此时间滞后样本所生成的集合预报误差可被定义为^[23]：

$ {\bf{x}}_{{\rm{T}}, k}^e = \frac{{\left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}}{{\sqrt {N - 1} }}, 0 < i < j \le N $

(3)

式中，x_i、x_j表示不同时刻初始场经过不同预报时长积分后得到的同一时刻的预报场，N表示时间滞后样本数，第k个差异场生成的预报误差可表示为$ {\bf{x}}_{{\rm{T}}, k}^e $。以间隔3 h同化并向前预报24 h的时间滞后方案为例，若每3 h保存一次预报场结果，则共有8个时间滞后样本。由8个时间滞后样本两两求差异，可得到x₈ - x₇, x₈ - x₆, ……, x₈ - x₁, ……, x₂ - x₁，共计28个差异场，再依据公式（3）可得到28个时间滞后样本估计的集合预报误差。

（2）历史预报样本。

历史预报样本由过去一段时间（如一个月）模拟的历史预报数据，每6 h保存的12 h和24 h预报场构建，并将同一时刻24 h预报场与12 h预报场的差异，作为历史预报误差的计算样本。历史预报样本所生成的集合预报误差可由如下公式定义^[27]：

$ {\bf{x}}_{{\rm{H}}, l}^e = \frac{{\left( {{\bf{x}}_{{\rm{H}}, l}^{{\rm{diff}}} - \overline {{{\bf{x}}_{\rm{H}}}} } \right)}}{{\sqrt {M - 1} }}, 0 < l \le M $

(4)

$ {\bf{x}}_{{\rm{H}}, l}^{{\rm{diff}}} = {\bf{x}}_{{\rm{H}}, l}^{{\rm{T}}1} - {\bf{x}}_{{\rm{H}}, l}^{{\rm{T2}}} $

(5)

其中，$ {\bf{x}}_{{\rm{H}}, l}^e $表示由优选历史预报样本估计的集合预报误差，M表示优选历史预报样本个数，$ {\bf{x}}_{{\rm{H}}, l}^{{\rm{diff}}} $表示预报时效分别为T₁、T₂的历史预报差异场，$ \overline {{{\bf{x}}_{\rm{H}}}} $表示该差异场的平均。

历史预报样本数量较多，且包含气候态信息，若直接全部引入，则与大气的“流依赖性”不符。因此，在同化阶段之前，本研究需优选出与预报时刻天气形势更相符的历史预报样本。具体做法是：首先，对于同化前的背景场x^b，与历史预报样本x_m (0 < m ≤ M)，其中M表示历史预报样本总数，提取两者中低层（500 hPa、700 hPa及850 hPa）的纬向风场U、经向风场V及温度场T，一一计算各个变量皮尔森相关系数，并以其加权平均$ Co{r_m} = \overline {\left| {Co{r_{u, m}}} \right| + \left| {Co{r_{V, m}}} \right| + \left| {Co{r_{T, m}}} \right|} \left( {0 < m \le M} \right) $, 为两者的相关。其次，将历史预报样本按相关系数从高到低进行排序，优选出相关系数更高的历史预报样本进行后续计算。

（3）样本组合。

综合两类样本的分析增量可由以下公式表示：

$ \delta {\rm{x = }}\delta {{\rm{x}}_1} + \sum\limits_{k = 1}^{\rm{N}} {\left( {{\alpha _k} \circ {\bf{x}}_{{\rm{T}}, k}^e} \right)} + \sum\limits_{l = 1}^{\rm{M}} {\left( {{\alpha _l} \circ {\bf{x}}_{{\rm{H}}, l}^e} \right)} $

(6)

其中，右侧第一项为静态的分析增量，第二项中的$ {\bf{x}}_{{\rm{T}}, k}^e $表示时间滞后样本估计的集合预报误差，第三项中的$ {\bf{x}}_{{\rm{H}}, l}^e $表示由优选历史预报样本估计的集合预报误差。N表示时间滞后样本数，M表示历史预报样本数。式（6）为式（2）基于扩展控制变量方法^[29]的多样本拓展形式，参考Wang等^[30]中的式（6-7），所对应的代价函数为：

$ \begin{array}{l} J\left( {\delta {{\bf{x}}_1}, \alpha } \right) = {\beta _1}\frac{1}{2}\delta {\bf{x}}_1^{\rm{T}}{{\bf B}^{ - 1}}\delta {{\bf{x}}_{\rm{1}}} + {\beta _2}\frac{1}{2} < {\alpha _l}, {\alpha _k}{ > ^{\rm{T}}}{{\bf{{\bf A}}}^{ - 1}} < {\alpha _l}, {\alpha _k} > + \frac{1}{2}{\left( {{\bf{d}} - {\bf H}\delta {\rm{x}}} \right)^{\rm{T}}}{{\bf{R}}^{ - 1}}\left( {{\bf{d}} - {\bf H}\delta {\bf{x}}} \right) \end{array} $

(7)

其中α_l和α_k分别表示历史预报样本和时间滞后样本对应的集合扩展控制变量的向量。

3 模式参数设置及试验方案设计

研究所用的模式为WRF（Weather Research and Forecasting Model）V 3.9.1以及WRFDA（Weather Research and Forecasting Model Data Assimilation System）V 3.9.1。区域分辨率9 km，水平格点数为271×271，垂直分41层，模式顶气压50 hPa。试验的初边界条件来自0.5 ° ×0.5 °的GFS（Global Forecast System）资料，模式参数采用：WSM -5微物理方案^[31]、YSU边界层方案^[32]，Kaine-Fritsch积云参数化方案^[33]、Dudhia短波辐射方案^[34]以及RRTM（Rapid Radiative Transfer Model）长波辐射方案^[35]。观测资料来自GTS（Global Telecommunications System）观测数据集（图 2）。

研究中静态的背景场误差协方差由NMC（National Meteorological Center）^[28]方法估计得到。统计样本为2015年7月1日00时—31日12时（世界时，下同）连续一个月每隔12 h WRF模拟的12 h预报场与24 h预报场。同化系统所用的控制变量为CV5方案^[36]，包括流函数ψ、非平衡速度势χ_u、非平衡温度T_u、假相对湿度RH_s及非平衡地表气压P_{s, u}。

根据同化试验中不同样本来源及样本估计的集合预报误差数，设计了四组对比试验（表 1），时间段为2016年7月1日00时—7日00时，每隔3 h同化一次，整个模拟过程共循环同化49次。分别是基于三维变分同化方法的控制试验（3DVar）；使用7个时间滞后样本得到28个时间滞后预报误差的混合同化试验（TL28）；使用12个时间滞后样本，得到66个时间滞后预报误差的混合同化试验（TL66）；以及结合了时间滞后样本与优选历史预报样本的混合同化试验（TH66）。在循环同化过程中，TH66试验在每次同化前均从所有历史预报样本中优选与当前天气形势相近的历史预报样本，引入同化系统。对比样本较为有限的TL66与TL28试验，TH66试验扩展了样本个数，且几乎不需要增加额外的计算资源。试验期间，3DVar、TL28及TL66预报时长为24 h，而在同化间隔均为3 h前提下，为获得更多的时间滞后样本，TL66的预报时长为36 h。

所有历史预报样本来自2015年6月7日00时—7月12日00时的预报数据，每隔6 h得到一组预报时效分别为12 h与24 h的同一时刻预报场，共进行了为期35天的样本构建，生成139组历史预报样本及对应的历史预报误差。所用的初边界条件均为GFS资料，不同化观测资料，且其余参数设置与试验模拟过程相同。

历史预报样本与同化时刻背景场的相关系数随时间的变化如图 3所示，黑色点线表示最低相关系数的变化，黑色细竖线表示各时刻38个优选样本的相关系数数值范围，蓝色点线表示各时刻引入的38个优选样本的平均相关系数。从图 3可以看出，优选样本的平均相关系数（蓝色点线）在7月1—4日的平均相关系数较为稳定，始终保持在0.7左右，7月4—7日平均相关系数有明显下降，最小值约为0.5，而最低相关系数数值（黑色实线）则在0.4左右波动。造成该现象的主要原因是本研究中可供优选的历史预报样本总数较为有限，当天气形势发生明显变化时难以始终选出高相关系数的历史预报样本。另外，由于不同时段构建的历史预报样本所包含的天气形势也不相同，这种差异会对同化结果产生一定影响，实际业务中若构建足够庞大的历史预报样本数据库，理论上可持续引入相关系数更高的优选样本。

4 单点观测理想试验

为探究不同预报样本对于同化系统的影响，本研究首先进行单点观测理想试验的分析。在模式层第9层（约800 hPa）的112 °E，31 °N处放置一个新息向量为1.0 m/s的纬向风U理想观测，同化时间为2016年7月2日00时，集合误差协方差权重设置为0.75，局地化尺度为200 km。风场的分析增量如图 4所示，可以看出，3DVAR单点试验的风场增量为各向同性（图 4a），而基于集合-变分混合同化方法的风场分析增量均表现为各向异性的特征：由于TL28试验（图 4b）的时间滞后样本过少，带来了较大的采样误差，分析增量分布较为杂乱；在TL66试验中当时间滞后样本数增加后，采样误差的问题被一定程度改善（图 4c）；TH66试验，通过时间滞后样本结合优选历史预报样本，分析增量表现出了与当前天气形势较好的关联性（图 4d）。

在模式第16层（约500 hPa）的118 °E，31 °N位置处放置一个新息向量为1 K的温度T理想观测，时间为2017年7月2日00时，集合权重及局地化尺度的设置与纬向风U的单点理想观测试验一致。分析水汽场的增量特征可知，3DVAR中CV5方案的背景场误差协方差不存在温度T与水汽q之间的变量相关（图 5a），混合同化试验存在T、q之间的变量相关性，且水汽增量的分布形态与背景场湿度的等值线较为一致（图 5b~5d），均能够引入“流依赖”的特征。

进一步进行了TH66方案不同时刻的温度T的单点理想观测试验（7月2—5日每隔24 h进行一次），集合误差协方差权重为0.75，局地化尺度为200 km。从图 6可以看出，TH66方案的水汽分析增量随时间不断变化，且分析增量能够随着背景湿度场发生改变（图 6a~6d），其“依流而变”的特征来源于不断更新的时间滞后预报样本及与当前天气形势更相近的历史预报样本，说明本方案能够仅依靠少量计算资源，达到引入合理的“流依赖”的背景场误差协方差的目的。

5 梅雨过程循环同化及预报试验 5.1 梅雨过程介绍

2016年长江区梅雨期从6月19日入梅，7月21日出梅，雨量为584.3 mm，较常年平均偏多107.9%，位于1951年以来第3位。7月1—6日的总降水量达到270.4 mm，已接近该区整个梅雨期的平均降水量（280.4 mm），其中前三天的降水量分别为60.4 mm、64.3 mm、52.5 mm^[37]。

本研究的模拟时间段为2016年7月1日00时—7日00时，每隔3 h同化一次，整个模拟过程共循环同化49次。

5.2 均方根误差

图 7为四组试验连续一周的分析场及预报场中各变量相对于ERA5（5^th generation climate reanalysis dataset from ECMWF）的平均均方根误差。从分析场看（图 7实线），蓝色实线所表示的TH66的各变量分析结果，明显优于仅依赖于时间滞后样本的混合同化试验TL28、TL66；从12 h预报场看（图 7虚线），蓝色虚线所表示的TH66试验各变量仍然能够保持较小的均方根误差，而仅采用时间滞后样本的TL28、TL66相较于3DVAR的改进稍小。

四组试验分析场与12 h预报场较探空观测资料平均一周的均方根误差如图 8所示。对于U、V变量，蓝线所表示的样本组合的混合同化试验TH66整体表现最优，其分析场和12 h的预报场均与实况探空观测更接近。时间滞后样本数较多的TL66表现略好于TL28，3DVAR与TL28试验整体表现稍差。从水汽q变量的分析场来看（图 8d），混合同化试验没有显示出优势，在12 h预报场中，基于样本组合的TH66试验（图 8d蓝色虚线）表现略好于其它三组试验。

5.3 降水循环同化预报FSS评分

为对比不同试验对降水预报的效果，本研究对四组试验（3DVAR、TL28、TL66、TH66），按不同阈值进行了逐6 h的平均累计降水FSS（Fraction Skill Score）评分（图 9），降水观测为中国地面与CMORPH融合逐小时降水产品（CHMPA）^[38]。在小雨量级中，3DVAR表现稍差，TL28与TH66表现稍优（图 9a）。在中雨量级中，TH66在整个预报时效内均表现最优，尤其在预报时效较长的24 h，TH66能够保持较好的表现，展现出较好的鲁棒性（图 9b）。在大雨量级中，预报初期TH66的评分表现稍低于TL28与TL66试验，但在24 h评分与TL66试验接近（图 9c）。在暴雨量级中，TH66的评分表现明显优于其它方案（图 9d）。综合来看，3DVAR在各降水阈值内的降水评分表现最差，由优选历史预报样本与时间滞后样本结合的混合同化方案TH66表现最佳。

5.4 降水分布及诊断

为更好地理解各方案对强降水过程的模拟，本研究挑选了在此次梅雨过程中降水强度较大时间段（2016年7月1日00时—2日00时）的24 h累积降水进行分析（图 10）。如图所示观测降水呈西南-东北向，累计24 h的降水量超过150 mm的区域主要集中在安徽中部地区（图 10a）。3DVAR试验模拟出了整体的雨带走势，但对于强降水中心的模拟表现较差（图 10b）。相较于3DVAR，TL28对于降水强度较大区域的模拟稍有改善，但在江苏南部地区存在过报的问题（图 10c）。增加时间滞后样本数的TL66试验对于安徽中部强降水的模拟有一定改进，但其强降水中心略微偏北，且在江苏南部仍然存在过报的问题（图 10d）。由优选历史预报样本组合时间滞后样本的混合同化方案TH66对于强降水中心的模拟在位置和强度上均最接近实况，并且改善了江苏南部降水过报的问题（图 10e）。

为进一步探究样本组合的混合同化方案较三维变分方案的降水模拟改进的原因，两组试验在2016年7月1日12时（图 11a、11b）和18时（图 11c、11d）850 hPa的水汽通量和水汽通量散度如图 11所示。参照累计降水分布图可知（图 10），两组试验在强降水中心附近均能模拟出较强的水汽通量以及水汽通量散度。在12时的850 hPa层上，3DVAR试验的水汽通量的切变线更弱且稍偏北（图 11a红线），而TH66试验具有更为明显的水汽通量切变（图 11b红线）。对比两组试验在18时的水汽通量散度分布可以看出（图 11c、11d），3DVAR试验的水汽辐合区明显向东偏移（图 11c红框），从而导致3DVar试验没有较好地模拟出强降水中心，而TH66试验的水汽辐合区仍能够较好地分布在降水中心（图 11d红框），这与实况降水分布特征更为接近。

6 总结

集合-变分混合同化方法是资料同化领域的研究热点之一。一方面，该方法能为同化系统提供“流依赖”的背景场误差协方差，另一方面，该同化方法中集合预报需消耗较多的计算资源。本文探究讨论了时间滞后样本的构建方案与历史预报样本的优选方案，通过计算风场与温度场的相关系数，引入与同化时刻更为相近的历史预报样本。研究进行了一系列单点观测理想试验和针对2016年梅雨期强降水过程（2016年7月1日00时—7日00时）的连续循环同化及预报试验。

单点观测理想试验的结果表明，时间滞后样本能够为集合-变分混合同化系统引入“流依赖”的背景场误差协方差；通过优选历史预报样本扩展有限数量的时间滞后样本后，风场增量的表现更为合理，缓解了时间滞后样本量不足所带来的问题。单点试验也表明，不断更新的时间滞后预报样本及优选得到的与当前天气形势更相近的历史预报样本，能够为同化带来“依流而变”的分析增量，说明本方案能够仅依靠少量计算资源，便达到引入合理的“流依赖”的背景场误差协方差的目的。

对于连续循环同化及预报试验，整体而言，采用时间滞后的集合-变分混合同化方案分析场与预报场均方根误差比三维变分更小；通过优选历史预报样本扩展的时间滞后样本后，其分析场和预报场的均方根误差相对于单纯采用时间滞后样本进一步减小。从降水预报效果上看，优选历史预报样本组合时间滞后样本的混合同化方案效果最好，仅使用时间滞后样本的混合同化方案效果次之，三维变分表现最差。通过分析降水个例过程中的水汽通量条件可知，样本组合方案较三维变分方案能够更好地模拟出水汽通量散度分布，使水汽辐合区与实况强降水区域位置更为接近，从而改善降水预报的表现。

本研究进行了将时间滞后样本与历史预报样本共同引入集合-变分混合同化系统的尝试，并取得了较好的研究结果，但对于该方法的集合协方差权重系数及局地化尺度的参数敏感性试验尚未进行，仍需进一步展开研究。此外，本研究主要构建了与梅雨过程相关的历史预报样本，并在此基础上进行了针对梅雨期强降水的研究，因此下一步拟构建包含更多天气过程、更长时间序列的历史预报样本，以改进不同类型天气过程的预报效果。