引言

近年来数值天气预报技术有了飞速发展，数值模式的分辨率也越来越高，其中对流可分辨模式分辨率可达1~4 km，取消积云对流参数化方案，对流过程可以显式表达(Yano et al，2018)，业务单位也越来越多地依赖对流可分辨数值模式来对局地高影响天气进行预报。集合预报技术(Epstein，1969；Leith，1974)作为一种数值天气预报中的概率预报手段，也逐步从全球大尺度天气预报拓展到有限区域以及对流可分辨尺度数值天气预报当中(陈静等，2005)。对于暴雨、强对流等局地强天气现象，其发生发展的动力和热力机制较为复杂，导致它们具有很强的预报不确定性，而对流尺度集合预报技术为解决此类强天气的概率预报问题提供了一个有效途径。

在集合预报研究中，如何获得合理的初值扰动是一个研究难点。集合预报初值扰动方法自全球集合预报开始已发展得较为成熟，如增长模繁殖法(BGM)(Toth and Kalney, 1993；1997)、奇异向量法(SVs)(Molteni et al，1996)以及集合变换卡尔曼滤波(ensemble transform Kalman filter，ETKF)(Wang and Bishop, 2003；Bowler et al，2008) 方法等，此类方法也沿用到了有限区域集合预报中(Stensrud and Fritsch, 1994；Stensrud et al，1999；Du et al，2003；Walser et al，2006；Bishop et al，2009；Bojarova et al，2011；张涵斌等，2014)，在初始时刻即可为区域集合预报产生足够的中尺度扰动信息；另外一些有限区域集合预报初值扰动采用大尺度集合直接动力降尺度(Marsigli et al，2005；Frogner et al，2006；Bowler et al，2008；张涵斌等，2017)，也能为区域集合产生一定的效果。此外，为解决有限区域集合预报中的多尺度不确定性问题，有学者开始尝试尺度混合初值扰动方法的研究，即将全球集合预报动力降尺度初值扰动和区域版本的BGM、ETKF等方法相结合，来为有限区域集合预报生成初值扰动，此种方法可以获得充分的大尺度和小尺度扰动信息(Caron，2013；Wang et al，2014；Zhang et al，2015；Keresturi et al, 2019)。

对流尺度集合预报是目前集合预报研究的热点和难点，其模式分辨率可达到1~4 km，关闭积云对流参数化方案，国内外关于对流尺度集合预报技术的发展尚处于相对初步的阶段(Johnson and Wang, 2016; Raynaud and Bouttier, 2016)，如何为对流尺度集合预报产生有效的初值扰动，是一个值得继续深入研究的科学问题，与中尺度区域集合预报采用的方法类似，对于对流尺度区域集合预报而言，目前国际上较为流行的方法有动力降尺度，如德国气象局2.8 km分辨率的对流尺度集合预报系统COSMO-DE-EPS(Kühnlein et al, 2014)，初值扰动来自四个业务全球模式动力降尺度生成，并做了中心化处理；英国气象局发展了2.2 km分辨率的集合预报系统MOGREPS-UK，初值、侧边界扰动均来自全球集合预报系统MOGREPS-G的动力降尺度扰动场，其中初值扰动做了中心化处理(Hagelin et al, 2017)；法国气象局的对流尺度集合预报系统AROME-EPS在2016年业务运行，系统水平分辨率为2.5 km，初值扰动来自短期集合预报系统PEARP的动力降尺度场(Nuissier et al, 2016)。动力降尺度方法较为简便易行，且能够产生较好的效果，但该方法的局限性是无法在初始时刻产生充足的小尺度扰动结构。

有学者也将中尺度区域集合预报中的初值扰动方法沿用到了对流尺度区域集合预报中，如沿用区域BGM方法(高峰等，2010；马申佳等，2018)、区域ETKF方法(庄潇然等，2016)以及混合扰动方法(庄潇然等，2017)等为对流尺度集合预报产生初值，也均取得了一定的效果。但是鉴于对流可分辨模式主要瞄准短预报时效(1~3 d)、局地天气的发生发展，因此观测资料的同化必不可少。集合资料同化(ensemble data assimilation，EDA)将资料同化与观测扰动相结合，是一种有效的初值扰动方法，目前在ECMWF全球集合预报系统中得到了应用(Buizza et al，2005)，而目前国内将EDA引入对流尺度集合预报系统中还未见报道。

目前华北区域中心业务运行确定性快速循环同化系统，为华北地区强对流天气的预报提供有效的业务支撑，但目前对华北地区强对流天气的发生发展不确定性的研究还较为有限，华北地区强对流天气概率预报的研究基础还较为匮乏，针对华北地区暴雨等强对流天气开展强对流集合预报研究势在必行。本文利用华北地区丰富的实时观测资料，探索华北对流尺度集合预报的EDA初值扰动方法。本研究不仅对华北对流尺度集合预报的发展具有重要意义，也可为强对流天气可预报性研究和业务对流尺度集合预报提供新方法新思路，具有较好的应用前景。

1 资料与方法 1.1 华北对流尺度区域集合预报系统

北京城市气象研究院2018年发展的对流尺度区域集合预报系统属于华北“睿图”模式体系的一个组成部分，该系统采用WRF模式(Weather Research Forecasting Model)V4.1.2版，模式区域设置为水平分辨率3 km，模拟区域范围为35.5°~46.3°N、105.2°~122.4°E(图 1)，共550×424个格点，覆盖华北大部分区域，垂直层次为51层模式面，模式层顶为50 hPa。该系统包括1个控制预报和20个扰动成员预报，共21个集合成员，系统每天从00 UTC和12 UTC起报两次，预报时效为24 h，预报输出间隔为6 h。所有成员物理过程参数化方案设置与业务确定性预报3 km区域相一致，即Thompson微物理方案、Mellor-Yamada-Janjic(MYJ) 边界层方案及RRTMG长短波辐射方案, 关闭了积云对流参数化方案。

1.2 集合资料同化初值扰动构建方案

传统的动力降尺度方法，将全球集合预报直接通过区域模式进行初始化来得到区域集合的初值场，并进行模式积分，此过程并未有资料同化过程，少了区域模式分析场的参与，因此传统的动力降尺度方法虽然简便，但具有一定局限性，首先是初值扰动尺度较大，与区域模式分辨率不匹配，其次是缺少观测信息，初值不够准确。

目前欧洲中期天气预报中心(ECMWF)集合预报系统采用奇异向量与集合资料同化EDA相结合产生初值扰动，EDA方法通过对所有成员进行资料同化形成多组分析初值，从而改善所有成员的初值质量。本文构建的集合资料同化方法，首先将NCEP全球集合预报系统(Global Ensemble Forecast System，GEFS)资料初值场进行动力降尺度，获得区域模式降尺度的初值场，并基于WRF三维变分模块(WRF data assimilation，WRFDA)，与多组观测进行资料同化，获得分析初值。目前华北业务数值预报系统中同化的观测包括常规地面站和探空、飞机报、地基GPS以及京津冀6部雷达径向风速度和反射率因子资料，考虑到非常规资料的影响较为复杂，本文中EDA的观测资料主要为华北地区常规地面站和探空观测。图 2给出了对流尺度集合预报系统框架: GEFS成员1动力降尺度得到WRF的背景场，利用WRFDA同化观测成员1，得到模式的分析场成员1，进行预报，以此类推，得到21个集合预报。

EDA方法中的变分同化采用传统的三维变分(3D-Var)。在3D-Var框架中，为了获得最优分析场x_a, 引入目标函数J来描述向量x与背景场x_b以及观测y的距离, 即目标函数可表示为:

$J(\boldsymbol{x})=J_{\mathrm{b}}(\boldsymbol{x})+J_{\mathrm{o}}(\boldsymbol{x}) $

(1)

对于基于3D-Var的一组集合预报的资料同化，为获得每个成员的分析场，第i个集合成员(包括控制预报ctl)的目标函数可表示为:

$\begin{aligned} J\left(\boldsymbol{x}_{i}\right)=& J_{\mathrm{b}}\left(\boldsymbol{x}_{i}\right)+J_{\mathrm{o}}\left(\boldsymbol{x}_{i}\right) \\ =& 1 / 2\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{\mathrm{b} i}\right)^{\mathrm{T}} \boldsymbol{B}^{-1}\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{\mathrm{b} i}\right)+\\ & 1 / 2\left(\boldsymbol{y}_{i}-\boldsymbol{H} \boldsymbol{x}_{i}\right)^{\mathrm{T}} \boldsymbol{R}^{-1}\left(\boldsymbol{y}_{i}-\boldsymbol{H} \boldsymbol{x}_{i}\right) \\ &(i=c t l, 1, 2, \cdots, 20) \end{aligned} $

(2)

式中: B和R分别为背景误差协方差矩阵和观测误差协方差矩阵，H为线性观测算子。本文中背景误差协方差矩阵的估计采用NMC方法获得，变分同化的控制变量为流函数、势函数、不平衡温度、假相对湿度以及不平衡表面气压。

EDA流程中，不同成员的同化的观测是基于常规观测资料进行扰动形成的多组观测。观测扰动构造过程如下:

第一步，获取WRFDA中的三维变分常规资料观测文件中的观测变量的值以及观测误差信息，其中观测变量包含气压、风速、风向、温度以及相对湿度；观测误差信息来自WRFDA中自带的观测误差统计量，这也是目前华北业务确定性数值预报系统中采用的观测误差。

第二步，利用随机数发生器生成20个正态分布的随机数R_k(k=1, 2, …, 20)，数值范围在-2~2。

第三步，原有的观测文件中包含多个变量在不同高度以及地面的观测，包括相对湿度、风向、气压湿度、风速，构成了所有观测样本，针对每一个观测样本，用对应的观测误差ε乘以20个随机数R_k(k=1, 2, …, 20)得到20个扰动，如对于500 hPa温度T，其观测误差为ε=1.0 K，如果随机生成的第一个随机数值为R₁=-1.2，则第一个观测扰动值为R_k×ε=-1.2 K。

第四步，原观测值O_ctl加上20组观测扰动值可获得20组新的观测，即扰动观测O_k:

$ O_{k}=O_{a l}+R_{k} \times \varepsilon \quad k=1, 2, \cdots, 20 $

(3)

20组经过扰动的观测与原观测共同构成21组观测文件，且扰动后的多组观测是相对独立的，用于每个集合成员的资料同化。

图 3是2018年7月15日00 UTC时次探空观测扰动的绝对值|R_k×ε|以及观测误差ε的分布，对于探空资料来说，在整层大气中，无论是风速、温度还是气压，观测扰动值基本上分布在探空观测误差附近，如对于风速(图 3a)，250 hPa左右的观测误差为3.4 m·s^-1，而观测扰动绝对值则在2.2~4.2 m·s^-1范围内分布；对于地面站观测来说情况也类似，地面温度的观测误差为1.1 K，而地面站点的观测扰动绝对值在0.7~1.4 K分布(图略)。

观测扰动R_k×ε是否正态分布也在一定程度上反映了观测扰动的质量。通过计算一段时间内(2018年7月1—15日)不同量级的扰动样本占总样本比率可获得观测扰动的概率密度分布，图 4给出了探空观测的850 hPa温度(T₈₅₀)、500 hPa纬向风速(U₈₅₀)扰动概率密度分布，由图 4a可以看出，T₈₅₀的观测扰动主要分布在-2~2 K附近，这主要是由于T₈₅₀观测误差统计为εT₈₅₀=1 K，而随机数R值在-2~2。此外可以看出观测扰动的概率密度分布与标准正态分布(均值为0，标准差为0.5)较接近，即扰动值越大的样本占比越小，而超过-2~2 K区间的扰动样本占比较小，几乎为0；对于U₈₅₀观测扰动(图 4b)结果也类似，表明构建的观测扰动基本合理。

至此21组扰动观测数据已构建完成，包括原有的观测数据(用于控制预报同化), 以及20组扰动观测数据(用于扰动成员同化)。

1.3 EDA初值扰动方案试验设计

为了验证集合资料同化方法在对流尺度集合预报系统中的作用，本文基于EDA初值扰动方法开展了对流尺度集合预报试验，并与动力降尺度初值扰动方法进行对比。系统模式分辨率及范围设置均与1.1节中介绍的一致。

设计了三种试验方案: 方案一基于GEFS多成员预报场的动力降尺度构建初值扰动并作为参照试验，称为DOWN方案；方案二为采用GEFS多成员动力降尺度初值与单一未扰动观测进行集合资料同化的方案，称为EDA-NOPO方案；方案三为GEFS多成员动力降尺度初值与扰动观测进行集合资料同化的方案，称为EDA-PO方案。从方案二与方案一的对比可以看出资料同化的作用，而方案三与方案二的对比可以看出观测扰动的作用。三组试验设置见表 1，试验时段选取2018年7月11日至8月11日，连续一个月。其中2018年7月16日在北京城区发生了对流强降水过程，作为批量试验中的重点分析个例。

1.4 试验数据

本文构建的对流尺度集合预报背景场和侧边界条件来自GEFS全球集合分析场和预报场，水平分辨率为1°×1°，该资料每天00 UTC、12 UTC获取2次，预报间隔为6 h。本文中高空和地面要素检验采用全国探空和地面观测站资料，这两种资料在整点时刻(00 UTC和12 UTC)站点较密，其中探空站在华北范围内有24个，地面观测站在华北范围内有420个，可较好地用于检验，而06 UTC和18 UTC时刻的探空和地面测站数量较少，因此仅采用整点观测进行检验，即检验间隔为12 h；降水检验采用全国基本降水观测站(约2 507个)，包括逐3 h和逐6 h累计降水。所有针对站点的检验均通过双线性插值方法将模式格点预报插值到观测站点上进行。

2 试验结果分析 2.1 扰动特征分析

本节具体分析EDA的引入会对集合初值场产生怎样的效果。图 5给出了三种方案500 hPa U分量风的集合初值场离散度分布以及500 hPa观测站点分布，从DOWN方案可以看出其离散度较大，尤其是内蒙古东部和宁夏地区，全场离散度可达1.26 m·s^-1；对于EDA-NOPO方案，各成员同化相同观测之后，初值离散度相对于DOWN方案呈现出显著减小特征，如内蒙古西部以及宁夏等地区，与该地区观测位置对应较好，全场整体的离散度为0.91 m·s^-1；而EDA-PO方案同化了扰动观测之后，虽然在内蒙古西部以及宁夏等地区离散度有所减小，但是某些位置有所增加，如大连东部地区、河北南部、山西和山东南部观测密集区域，体现了观测扰动的作用，使得EDA-PO方案全场离散度达到了1.2 m·s^-1，虽然不及DOWN方案，但是显著优于EDA-NOPO方案。以上结果说明资料同化减小了背景场的不确定性，而观测扰动有效地表达了观测的不确定性。

为进一步研究EDA方法对初值扰动的影响，本节对三种方案扰动增长特征进行分析。对于一个M_i×M_j的二维格点场引入平均绝对扰动:

$ P=\frac{1}{M_{i} \times M_{j}} \sum\limits_{i=1}^{M_{i}} \sum\limits_{j=1}^{M_{j}}\left(\frac{1}{N} \sum\limits_{k=1}^{N}\left|I C_{k}^{i j}+I C_{c t l}^{i j}\right|\right) $

(4)

式中: IC_k^ij为第k个成员在格点(i, j)的初值, IC_ctl^ij为对应得控制预报的初值, N为扰动的集合成员数(本研究中为20个)。

从逐6 h的平均绝对扰动演变特征(图 6)来看，在初始时刻DOWN方案具有很大的扰动幅度，且扰动增长较为明显，12 h扰动最大值可达2 m·s^-1；EDA-NOPO方案因为引入了观测资料同化，使得扰动能量在各个层次相对于DOWN方案有所减小，12 h扰动最大值也未超过2 m·s^-1；而EDA-PO方案由于进行了观测扰动，初值扰动幅度相对于EDA-NOPO方案有显著提高，在随后时刻虽然其扰动有所增长速率较快，12 h预报时效扰动最大值可达2 m·s^-1左右，与DOWN方案相当。从24 h预报时效来看，DOWN和EDA-PO的扰动最大值均可达3.2 m·s^-1，而EDA-NOPO仅达到3 m·s^-1，说明同化单一观测会对集合扰动显著削弱，而同化扰动观测会最大限度上保持扰动的振幅。

2.2 集合预报检验结果

对批量试验时段内的结果进行评分统计，进一步定量分析EDA方法的预报效果。由于EDA-NOPO方案离散度较低，且不作为主要分析对象，这里分析的EDA方案特指EDA-PO方案。采用的集合预报检验方法为集合平均的均方根误差(RMSE)、集合离散度以及连续等级概率评分(CRPS)。批量试验检验主要包括等压面要素、地面要素以及降水三个部分。

RMSE和集合离散度(Spread)是集合预报最常用的检验方法。本文首先对低空和近地面温度、纬向风风场进行检验。图 7给出了试验时段内统计的850 hPa温度和纬向风(T₈₅₀, U₈₅₀)、2 m温度(T_{2 m})以及10 m纬向风(U_{10 m})的检验结果。从图 7可知EDA方案的RMSE在0~24 h时效内均小于DOWN方案，尤其是短预报时效最为明显，如12 h预报时效，EDA集合的RMSE为1.3 K，而DOWN集合的RMSE为1.55 K，说明EDA方法有效地降低了集合预报短预报时效的误差；从集合离散度来看，DOWN方案的离散度在短时效内略优于EDA方案，如DOWN方案T₈₅₀的12 h集合离散度为0.61 K，而EDA方案的12 h集合离散度为0.52 K，对于其他层次和要素的检验结果也较为类似，这里不再赘述。以上结果说明EDA方法可以有效提高集合预报的准确性，减小预报误差，离散度会略有降低，与图 5中得出的结论类似，但是整体集合预报的可靠性会有所提高。

CRPS评分(Hersbach，2000)是能够有效衡量集合预报中多种事件预报准确性的综合指标，是对集合预报质量的整体评价，评分越小表示集合预报质量越好。从CRPS评分来看(图 8)，不管是等压面要素还是近地面要素，EDA方案的CRPS值均小于DOWN方案，说明EDA方案能够显著提高对流尺度集合预报的准确性。

2.3 降水预报检验结果

为了检验EDA初值扰动方法的预报效果，研究了试验时段内一次典型的强降水个例。图 9a给出了观测的2018年7月15日18 UTC至16日00 UTC的6 h累计降水。可以看出该个例在华北地区表现为明显的东北—西南向雨带分布，尤其是在北京地区具有明显的大值区，某些站点的6 h累计降水可达60 mm以上。图 9b, 9c给出了DOWN和EDA两种集合方法从2018年7月15日12 UTC起报、对该降水个例时段6 h累计降水大于4 mm量级概率预报。DOWN集合(图 9b)由于没有资料同化，对实况的位置和范围均存在较明显的漏报；对于EDA集合(图 9c)，大概率区域较好地给出了实况观测的强降水中心位置，北京地区大于4 mm降水概率局地可达70%以上，有效地改善了DOWN集合的漏报现象。此外两种方案对河北、山西、内蒙古交界区域均表现出了较大的降水概率，相对于实况均略偏北，其中业务确定性模式强降水也在该地区存在一定的偏北(图略)，主要由于WRF模式对本次系统性降水过程的环流形势预报略有偏差，尤其是西南气流位置偏北所致，这里不再做深入分析。

为了探索对流尺度集合预报对局地降水的量级以及时段的识别能力，研究了2018年7月15—16日降水过程中北京地区逐3 h累计降水变化。图 10给出了2018年7月15日15 UTC至16日12 UTC的北京区域内平均的逐3 h观测降水量及成员预报降水量。图中的观测降水量演变可以看出该降水过程在北京地区具有明显的强降水时段，一个是15日21 UTC左右，一个是16日12 UTC左右，这两个时段北京地区平均3 h累计降水可达8 mm以上。从集合成员预报可以看出，DOWN集合各成员预报不够发散，且对这两个时段的降水变化描述不够准确，9 h预报显著低估了实况降水，18 h预报又明显高估了实况；EDA成员来看，首先各个成员预报较为发散，9 h预报有若干成员非常接近实况量级，对于18~24 h第二次强降水时段，各个成员预报相对于DOWN集合要更加接近。

我们采用ROC面积(AROC)和Brier评分(BS)来评估两种集合扰动方案在试验时段内对所有降水个例预报的统计检验结果，其中AROC代表了某一量级降水命中率和假警报率之间的关系，AROC越大说明概率预报的命中率越高；BS评分衡量了预报概率和事件发生频率之间的距离，值越小说明概率预报表现越好。图 11给出了2018年7月11日至8月11日这连续一个月内所有降水个例统计计算的0.1~13 mm三个降水量级的逐6 h降水AROC和BS评分演变。可以看出，不管是小雨、中雨还是大雨，EDA集合的AROC评分均优于DOWN集合，尤其是12~18 h预报时效的优势更加明显，24 h预报时效两种集合的AROC表现趋于一致；BS评分情况也类似，EDA集合方案在所有量级均显示出相对于DOWN集合的优势，其中12~18 h预报时效的改进最为明显，说明EDA方法能够显著改善DOWN的降水预报效果。

3 结论与讨论

本文基于华北对流尺度集合预报系统探索了集合资料同化EDA的构建方法，该方法利用多成员变分同化，将全球集合动力降尺度初值场与观测扰动相结合，为对流尺度集合预报获得更优的初值扰动场。设计了三组对比试验，即动力降尺度初值扰动，多成员同化单一观测的EDA，以及多成员同化不同观测的EDA。从扰动特征、概率预报检验以及降水预报检验等方面进行了综合评估。得出以下结论:

(1) 设计了针对华北对流尺度集合预报的EDA初值扰动方案，该方案以GEFS全球集合分析场作为背景场，集合成员分别同化常规观测资料构建对流尺度集合预报的EDA初值场。对观测资料设计了随机扰动方案，观测扰动幅度能够与观测误差相一致，且呈现正态分布特征。

(2) EDA方法会改变初值场中要素分布特征，进而改变初始离散度分布，使得初始扰动可以有效减少初值场中背景场的不确定性，同时体现观测的不确定性。相对而言，观测扰动会最大程度上保持充分的初始离散度且扰动具有较好的增长能力，而不对观测进行扰动则会限制初始集合离散度以及其增长效果。

(3) 集合预报检验结果表明相对于动力降尺度，EDA初值扰动方法可以大幅减少短预报时效的预报误差，离散度略有减少；对于对流尺度集合预报的概率技巧评分具有较明显提升。

(4) 降水试验结果表明EDA集合能够显著改善动力降尺度集合对强降水的漏报现象，对于局地降水的量级和时段具有更准确的预报能力；试验时段内的统计降水评分结果也表明不管是小雨、中雨还是大雨量级，EDA方案均能够比动力降尺度获得更好的概率预报评分效果。

以上结论表明EDA初值扰动方法在华北对流尺度集合预报中是行之有效的，该方法充分利用华北地区丰富的观测资料，可获得能够反映真实分析场中背景场和实际观测不确定性, 而且扰动振幅更合理的初值扰动场，因此具有较好的实际应用价值。本文开展试验基于常规观测资料进行集合资料同化，考虑到业务应用需要与雷达资料同化相结合，未来需要进一步开展试验。