引言

暖季期间，我国南方暴雨频发，易引发洪涝灾害，如何提高降水预报落区和强度的准确率是天气预报员一直关注的难点问题。中国气象局2013年3月公布的天气研究计划(2013—2020年)将提高暴雨等重大灾害天气预报能力作为核心目标之一，认为目前我国的研究水平在暴雨等灾害性天气的精细化预报、延伸期预报、预报准确率等方面，还存在明显的预报难点和技术方法瓶颈。并指出大量常规观测资料和高时空分辨率的雷达等观测信息，在天气预报中的综合应用和定量应用技术非常欠缺，资料同化和精细模式技术与国际先进水平的差距呈加大趋势。天气预报中的这些难点问题对科学研究工作提出了更高的要求：缩小差距必须依靠科学技术进步，以先进的气象科研成果推动业务发展。因此，中国气象局近年着重推动气象预报预测业务与科研结合试点工作，使科研工作能够围绕业务难点和重点开展，科研成果能支撑和指导业务，实现理论与实践的紧密结合。为了促进气象科研与业务结合，2013年起中国气象科学研究院与国家气象中心通过暖季试验、中小尺度灾害性天气研究和学术交流等方面开展合作，并将暴雨个例预报和分析作为重点，将科研成果转化到业务工作中，从而提高预报效率。

多种观测资料不能得到充分利用，资料同化方法的不完善，均制约着预报水平的提高。作为业务科研结合的试点工作，本研究针对2013年暖季试验中的一次强降水过程，使用国外成熟的WRF-EnKF实时预报系统，同化常规探空观测资料，进行确定性降水预报，检验EnKF方法同化探空资料对模式初值的改进能力，以及该系统对我国南方降水的预报效果。

目前，数值预报已成为天气预报的核心工具，其准确程度是一个国家气象现代化水平的重要标志。数值天气预报本质是一个初值问题，即依靠大气当前的状态，用模式模拟将来的演变和发展。提高模式初始条件最有效的方法就是尽可能多地、合理地吸收观测信息，数据同化技术就是获取模式准确初值的有效途径。近年来，多种资料同化方法逐步发展，广泛应用，成为了现代数值天气预报的核心(邹晓蕾，2009；熊春晖等，2013)。从最早期的多项式函数拟合方法(Panofsky，1949)、逐步订正法(Bergthorsson et al, 1955)、到最优插值法(Gandin, 1965)，三维变分(Lorenc, 1986)和四维变分方法(Rabier et al, 2000)，以及集合卡尔曼滤波资料同化方法(Evensen, 2003)，观测资料的处理技术和质量都获得了很大的提高，也促进了模式预报水平的进步。同化方法中，背景误差协方差的估计是关键。集合卡尔曼滤波方法(ensemble Kalman filter，EnKF)将卡尔曼滤波与集合预报结合在一起，通过预报集合的样本方差估计得到随气流变化的背景误差协方差，成为当前数据同化的主流方法之一(张涵斌等，2014)。EnKF方法的应用领域非常广泛，从理想数值模式到实际数值模式，从对流尺度系统的预报到全球尺度系统的预报(Snyder et al, 2003；Tong et al, 2005；Whitaker et al, 2008；Zhang et al, 2009)。

同化技术对数值天气预报的关键作用在中国也受到了越来越多的重视(龚建东，2013)。近年来，国内的相关研究在使用不同的同化方法、同化多种观测资料方面取得了很大进展：WRF-3DVAR同化系统同化常规和地面观测资料(莫毅等，2008；王益柏等，2014)；多尺度/分块逐批资料同化技术同化海洋气象观测资料(万齐林等，2012)；WRF集合均方根滤波(EnSRF)同化系统同化多普勒雷达资料(闵锦忠等，2013)等。基于EnKF流依赖背景误差协方差的优势，众多研究将不同模式与EnKF方法耦合，以改进模式初值，提高预报准确率。兰伟仁等(2010a；2010b)利用雷暴尺度ARPS模式及Tong等(2005)发展的EnKF系统，分别探讨了不考虑、考虑模式误差的情形下同化多普勒雷达资料的效果问题。庄照荣等(2011a；2011b)基于GRAPES模式建立了EnKF同化系统，能够成批同化常规观测资料，与三维变分相比其降水预报更接近实况。郝民等(2014)在GRAPES同化系统中同化了L波段秒级探空资料，发现高分辨率资料更好地改进了模式预报初值，降水预报评分提高，预报偏差明显变小。秦琰琰等(2012)将EnKF同化技术应用到对流系统中，进行了基于WRF模式的多普勒雷达径向风和反射率因子同化试验。张晶等(2014)利用LAPS系统同化GPS(Global Positioning System)/PWV(Precipitable Water Vapor)资料，发现同化GPS/PWV资料对LAPS湿度场有显著改善，对随后的确定性预报产生正影响。

美国宾州州立大学张福青和其团队发展的WRF-EnKF(Real-time Penn State WRF-EnKF System)实时预报系统经过多年改进，已经非常成熟，应用在不同地区的台风预报和敏感性试验中，体现了其优越的性能(Xie et al, 2013；Munsell et al, 2014)。早在21世纪初，张福青和其合作团队就通过理想模式试验验证了EnKF对于对流尺度雷达资料同化的可行性(Snyder et al, 2003；Zhang et al, 2004)；还将EnKF方法与WRF模式结合，发现参数化集合可以部分消除模式误差的影响(Meng et al, 2007)；并采用合理的雷达资料稀疏化方法，证实真实雷达径向风速度的EnKF同化能明显改善模式的初值(Zhang et al, 2009)；随后通过大量深入系统的研究工作改进WRF-EnKF系统的性能，将飞机雷达观测资料同化进入初始场，并首次把该对流尺度集合分析系统应用于飓风实时预报，证实了其在飓风路径和强度预报准确率方面的优势，其性能明显优于其他预报系统(Zhang et al, 2011；Weng et al, 2012)；除了常规资料和雷达资料，卫星反演的云迹风资料也被同化到该系统之中(Zhang et al, 2012)。总之，该WRF-EnKF资料同化系统能够同化地面观测资料、探空资料、雷达径向风等多种资料，并且能够改善模式初值，大幅度地提高飓风的预报准确率。

目前，该WRF-EnKF系统已经成功移植到中国气象局新一代IBM高性能计算机上，并针对我国南方(关键区域为广东)暖季降水设计了实时确定性预报方案(图 1a)。为了检验该系统对于我国南方降水的预报能力，本文使用这一系统同化探空观测资料，对2013年5月14—16日我国南方一次暴雨过程进行确定性预报试验。

1 资料同化方法和数值模拟设计 1.1 资料同化方法

WRF-EnKF系统基于中尺度数值模式WRF-ARW(Skamarock et al, 2005)建立，通过WRF三维变分资料同化系统(WRF-3DVAR)产生扰动而得到初始集合成员，扰动变量包括水平风要素、位势温度、位势高度和水汽混合比，其他的预测变量如空气垂直速度、云水、雨水、雪、霰混合比等不被扰动；使用集合卡尔曼滤波(EnKF)方法进行资料同化。

集合卡尔曼滤波(EnKF)资料同化方法是1994年Evensen(1994)提出的对经典Kalman滤波(KF)资料同化方法的近似，是集合预报和资料同化的结合。EnKF的最大优点是可以利用蒙特卡罗方法给出背景误差统计特征的预报，这样可以方便地得到依流型而变(flow-dependent)的背景误差协方差，避免了静态背景误差协方差的缺陷，并且总体简单化，容易实现并行计算。Zhang等(2004)在经典理论的基础上将Kalman filter方程改进如下：

$ {X^a} = {X^f} + \mathit{\boldsymbol{K}}\left({y - H{X^f}} \right) $

(1)

式中，X^f代表没有进行资料同化的模式场(预报场)，X^a是加入资料同化之后的模式场(分析场)，y是大气的真实状态(观测场)，H是由模式量到观测量之间转换的观测算子，K是Kalman增益矩阵，定义如下：

$ \mathit{\boldsymbol{K}} = {P^f}{H^{\rm{T}}}{(H{P^f}{H^{\rm{T}}} + R)^{ - 1}} $

(2)

式中，P^f和R在这里分别代表背景和观测误差协方差。集合卡尔曼滤波同化方法的关键优势是利用一组预报集合的样本方差来统计出随天气形势变化的背景误差协方差P^f：

$ {P^f} = \frac{1}{{N - 1}}\sum\limits_{i = 1}^N {\left[ {x_i^f - {{\overline x }^f}} \right]{{\left[ {x_i^f - {{\overline x }^f}} \right]}^{\rm{T}}}} $

(3)

式中，N表示集合样本数，x_i^f是第i个集合成员的预报场，X^f为集合平均值，T表示共轭矩阵运算。简单而言，同化分析场X^a为观测场y与预报场X^f之差对预报场的加权修正，权重系数里P^f为背景误差协方差，EnKF的背景误差协方差通过预报集合的样本方差来估计，它是随气流形态变化的，同时也与模式误差密切相关，因而能够合理估计模式误差。

1.2 数值模拟设计

本文采用WRFV3.5.1模式进行两组试验：一组为不同化任何观测资料的控制试验(NODA)；另一组为同化探空资料的同化试验(DA)。两组试验均选用2013年5月14日08时(北京时，下同)起报的6 h间隔的1°×1° NCEP GFS(globle forecast model)60 h预报数据(预报至5月16日20时)作为初始和边界条件。模拟采用三重双向嵌套(图 1a)，水平分辨率分别为40.5、13.5和4.5 km，垂直层43层。云微物理过程选用Morrison双参数方案，辐射方案选用RRTM长波辐射方案和Dudhia短波辐射方案，边界层方案选用YSU方案，陆面过程选用thermal diffusion方案。第一、第二层嵌套采用Kain-Fritsch积云对流参数化方案，第三层嵌套关闭了对流参数化。控制试验模拟时段为2013年5月14日08时至5月16日20时，共积分60 h。同化试验的物理方案配置与控制试验相同，同化的探空资料包含11个等压面上的温度、湿度和风场资料。同化试验在5月14日08时通过WRF-3DVAR扰动产生30个集合成员，并分别向前积分12 h，然后在5月14日20时、15日02时和15日08时三个时次(每6 h一次)通过EnKF方法循环同化探空资料，之后使用15日08时的集合平均分析场作为模式的初始条件进行确定性预报，向前预报36 h(到16日20时)。由于此次过程降水范围较大，本文使用网格D2的模式输出结果对定量降水预报(QPF)进行分析，图 1b给出了D2中被同化进入系统的探空站位置分布，其中14日20时同化了53个站(实心圆+空心圆)，15日02时同化了4个站(空心圆)，15日08时同化了55个站(实心圆+空心圆+菱形)。

2 天气过程及背景场

2013年5月14—16日，江汉南部、江淮南部、江南大部、华南北部及东部，自西向东，自北向南出现大范围暴雨、局地大暴雨天气。14日08时起500 hPa上有短波槽从西南地区东部向东移动(图略)，15日08时加深为切断低涡(图 2a)；850 hPa上，14日开始四川盆地附近有低涡活动(图略)，15日08时两湖盆地附近又出现新的低涡，低涡东部至长江中下游地区存在一条暖式切变线，并伴随西南急流的显著加强(图 2b)。总体而言，15日08时，长江中游地区各层次都有低涡活动，系统相对深厚，低涡后部冷空气南压，与华南地区的西南暖湿气流相互对峙，形成有利于降水的天气形势。15日08时至16日08时，降水主要发生在江南和华南地区(图 3a)，本文主要检验WRF-EnKF系统对这段时期降水的预报效果。

3 定量降水预报效果检验

为检验WRF-EnKF系统确定性降水预报效果，本节将通过降水落区对比和TS评分客观检验方法对控制试验和同化试验的定量降水预报进行对比检验分析。参考公颖(2007)的划分方法，雨量等级划分为以下5个等级：小雨(0.01~9.99 mm)，中雨(10~24.99 mm)，大雨(25~49.99 mm)，暴雨(50~99.99 mm)，大暴雨(100~249.99 mm)和特大暴雨(＞250 mm)。由于此次天气过程的定量降水预报主要检验时段中(5月15日08时至5月16日08时)地面观测站点的24 h累积降水量没有出现＞250 mm的特大暴雨，因此本文暂不考虑特大暴雨的TS评分情况。

图 3a为5月15日08时至16日08时的24 h累积降水量实况分布图，雨带主要位于江南和华南地区，呈东北—西南向分布，存在南、北两个超过50 mm的暴雨区，南部暴雨区主要位于广西、广东、湖南和江西交界处，影响范围较广，且强度大，多数地区降水超过100 mm，广东中北部的降水中心强度接近250 mm；北部暴雨区主要位于江西和福建北部，范围较小，强度较弱，仅局部地区降水达到大暴雨等级。图 3b和3c分别为控制试验(NODA)和同化试验(DA)预报的24 h累积降水量分布图，由图可见，控制试验预报出了雨带的大体走势，但对于南部暴雨区预报范围偏小，量级略偏大；北部雨区范围偏大，强度显著偏强，大部分地区预报降水量超过100 mm，达到大暴雨级别。同化探空资料后(图 3c)，降水预报在强度和落区上均有较大改进。其中，南部暴雨区范围比控制试验有所增大，暴雨中心强度略减弱，更接近于实况；北部暴雨区的预报改进更为显著，强度和范围都大幅变弱、变小，尽管对于大暴雨落区的预报仍略大于实况，但位置和强度已和实况很近似。总体而言，对于24 h定量降水确定性预报同化试验与控制试验相比，在暴雨范围和强度预报上均有明显的改善。

为进一步定量化评估同化试验对定量降水预报的改进效果，图 4给出了控制试验和同化试验2013年5月15日08时至16日08时的24 h降水预报TS评分。TS评分值在0~1变化，其值越大表明对降水的预报准确率越高。从24 h降水预报TS评分结果看，控制试验和同化试验的TS评分均随着降水量级的增大而降低，但同化试验的TS评分值，除了在小雨量级上比控制试验略低之外，其他降水量极的评分均高于控制试验的结果。其中，大暴雨量级的TS评分提高幅度最大，从同化前的约0.03提高到同化后的约0.11，其次是中雨和大雨。暴雨TS评分的增幅略小，大约提高了10%。

4 EnKF方法同化探空资料对模式初值的改进

同化试验通过EnKF方法同化探空观测资料后得到了集合平均的分析场，并将此分析场作为初始条件进行确定性预报。为检验EnKF方法同化探空资料后对模式初值的改变情况，借鉴Bao等(2013)的方法，分别统计了2013年5月15日08时同化试验的分析场(即确定性预报的初始场)和同一时刻控制试验的预报场中各要素的平均值、相对于观测资料的偏差和均方根误差的垂直变化(图 5)。检验区域为区域D2，检验要素包括纬向风(U)、经向风(V)、温度(T)和相对湿度(RH)。具体方法如下：将区域D2内55个探空站点的观测资料作为参照(图 1b)，把控制试验和同化试验的模式格点场插值到探空站点上，分别计算观测(OBS)、未同化资料的控制试验(NODA)和同化试验(DA)中各要素的站点平均值(图 5a、5b、5c、5d)；NODA和DA中各要素与观测值之间的偏差站点平均值(图 5e、5f、5g、5h)；以及NODA和DA相对于观测的均方根误差站点平均值(图 5i、5j、5k、5l)。

从各要素的站点平均值及其与观测的偏差来看，NODA和DA试验均模拟出了与观测近似的垂直分布特征，且同化后的各要素偏差与未同化试验相比总体上呈现减小的趋势。

具体来看，观测纬向风整层为西风，由低层到高层逐渐增大，150~200 hPa达到峰值(大约20 m·s^-1)。对流层中低层(400 hPa以下)，DA试验纬向风和观测偏差与NODA相比差别不大，在近地层附近(1000~925 hPa)都存在大值(大约-2 m·s^-1), 在对流层中层(700~400 hPa)均在-1 m·s^-1~0间变化；对流层高层(300~100 hPa)两者差别变得显著，除200 hPa附近，DA偏差均小于NODA。与纬向风峰值对应，在150~200 hPa，两个模式与观测的偏差也达到最大，且DA试验的偏差(1.5 m·s^-1)与NODA(2.3 m·s^-1)试验相比减小幅度明显(图 5a和5e)。

由探空观测平均可知，经向风低层为弱南风(5 m·s^-1)；高层盛行北风，150 hPa附近达到峰值(-8 m·s^-1)。并且，除了平均经向风峰值(150 hPa)附近之外，DA试验的经向风垂直曲线在各层都比NODA更接近OBS。与之对应，DA的偏差比NODA显著减小，其中，400~200 hPa间改善最显著，250 hPa附近偏差大小从NODA的3 m·s^-1减小到DA的1.5 m·s^-1(图 5b和5f)。

对于温度和湿度的站点平均垂直廓线，两个试验均模拟出了其随高度降低的变化趋势; 从偏差大小上看，DA试验的偏差比NODA明显减小(图 5c、5d、5g、5h)。其中，NODA与观测的温度偏差在-1~1℃之间波动，同化探空后整层的温度偏差均明显小于NODA，减小幅度高达50%。而NODA的湿度偏差200 hPa以下为正偏差，近地层和250 hPa附近最大，约为12%，200 hPa以上由正转负，150 hPa为-10%。同化探空后，200 hPa以下DA试验的湿度偏差均比NODA有不同程度的改善，其中250 hPa偏差峰值附近减小幅度最大：偏差由NODA的12%减小到DA的6.5%。

与偏差的整体减小相对应，同化了探空资料后，纬向风、经向风、温度和相对湿度相对于观测的均方根误差均有所减小，并且其改进效果更明显(图 5i~5l)。两个试验的四个要素在均方根误差的垂直分布特征上均较相似，但从减小幅度看，风场在中高层的减小幅度较大，其中，纬向风减小幅度在250 hPa附近达到峰值(约2.8 m·s^-1)，DA试验的均方根误差比NODA减小了一半；经向风最大减小幅度在200 hPa附近，由NODA的6.7 m·s^-1减小到DA的4 m·s^-1。与风场不同，温度均方根误差的最显著改善位于低层和高层(最大改善幅度为0.7℃)，中层的改善较小(0.15~0.4℃)。相对湿度均方根误差在高层和低层的减小幅度都不大，其最显著改善位于中层400 hPa左右，最大减小幅度为12%。

总体而言，通过EnKF方法同化了探空资料后，DA试验初始条件的各气象要素相对观测的偏差和均方根误差均较NODA试验相同时刻的预报场有不同程度的减小，风场和相对湿度在对流层中高层的改善最明显；温度场整层都有改善，低层和高层改善更显著。

为了更加直观地分析WRF-EnKF系统同化探空资料后对大气环境场的改变情况，图 6和图 7分别给出了两组试验2013年5月15日08时的850和500 hPa温度场、相对湿度场、风场水平空间分布，及相应的探空观测值，以及DA试验相对于NODA试验的增值(DA减去NODA)。而为了更清晰地对比NODA和DA试验对风场的预报和分析效果，将两次模拟试验的格点资料分别插值到站点上，并绘制矢量风，直接与观测进行对比。从图中可以看出，NODA试验基本上预报出了探空观测各要素在两个等压面上的总体分布形势，但在具体数值上均与探空观测值存在较大偏差(图 6a、6b、6c和图 7a、7b、7c)。同化探空资料的DA试验分析场在减小偏差方面有明显效果，其各要素的分析值均比NODA更接近观测(图 6d、6e、6f和图 7d、7e、7f)。

具体来看，850 hPa上NODA试验的温度值在32°N以北以及东南沿海地区均低于观测，其中江浙附近低估最显著；而在32°N以南的内陆地区，NODA的温度值总体高于探空观测，且高估在重庆、湖北、广西、湖南的探空站点附近表现最为明显(图 6a)。与NODA相比，DA分析场的温度在D2区域北部边缘和东部沿海显著增高，云南贵州附近也较观测略偏大，其余的内陆地区则呈明显偏低趋势(图 6g)。可见，DA试验通过EnKF方法同化探空资料对温度场进行了调整，减小了与探空观测之间的偏差(图 6d)。DA试验的分析场相对于NODA相同时刻的预报场的正增量和NODA与观测的负偏差呈现了很好的对应关系，反之亦然(图 6a和6g)。NODA预报场的850 hPa湿度与观测的偏差大值区主要位于湿度梯度较大的区域，例如：苏浙皖沪地区有干湿空气交汇，NODA的高湿度区的相对湿度值比观测值偏低而低湿度区的相对湿度值偏高。在广西附近，其空气较干，NODA预报场的湿度值要远低于观测(图 6b)。相对于NODA试验，DA分析场的湿度值较大幅度的降低主要发生在NODA与观测正偏差显著的地区(苏浙皖沪地区的低湿度区)，而大幅度升高主要发生在NODA与观测负偏差的地区(苏浙皖交界处和广西附近)(图 6e)。因此，从850 hPa湿度来看，DA分析场与观测之间的偏差幅度较NODA试验有显著减小，即DA分析场比NODA预报场更接近实况(图 6b)。850 hPa的风场，在低涡附近NODA试验的风矢量和观测之间在风向(箭头夹角)和风速(箭头长短)上都有较大偏差，而同化探空后，DA试验的风矢量与观测风矢量之间的夹角减小、大小更接近。可见DA试验的分析场对风场也有很好的调整，其更接近观测。同样，同化探空资料对位势高度场也有明显改进，DA试验的位势高度值比NODA试验普遍增大，在低涡中心及其西北部增幅最明显，北部略有减小，导致低涡强度减弱，更接近于探空观测(图 6c、6f、6i)。低涡减弱改进了降水预报：NODA试验对该地区降水预报偏强，DA试验预报的降水强度降低。

对于500 hPa温度、相对湿度和风场，DA试验相对于NODA试验具有同样显著的改善：DA与观测的偏差减小，更接近观测场(图 7)。例如：与NODA相比，DA的温度值在区域西北部显著降低，南北温差增加；对于湿度场，由于NODA相对于观测的500 hPa相对湿度偏差比850 hPa层明显大很多，因此DA试验对500 hPa湿度的调整幅度也较850 hPa大很多，区域东北部西北-东南走向的高湿度带经过调整后宽度变窄、强度加强，两侧空气更干，湿度梯度增加；此外，DA试验分析场调整后的风矢量的风向和风速均更接近观测，与观测的夹角减小、风速大小偏差减小；对于位势高度场，NODA中北部的横槽在DA中加深为一低涡，其他区域高度场略有增加，系统强度增强。

5 结论和讨论

为了检验WRF-EnKF系统对我国南方降水确定性预报的应用效果，本文使用2013年5月14日08时起报的每6 h间隔的NCEP GFS 60 h预报场作为数值预报的初、边值条件，进行了一组对比试验：控制试验(NODA)和同化试验(DA)。控制试验不同化任何资料；同化试验通过EnKF方法循环同化三次常规探空资料，并使用最后一次同化时刻(15日08时)30个成员的集合平均分析场作为初始场进行确定性预报。将两组试验的降水预报和实况进行对比，并分析了同化试验末次同化的分析场(即确定性预报的初始场)相对于控制试验同一时刻的预报场的改进情况。结论如下：

(1) WRF-EnKF系统同化探空观测资料可以改善暴雨落区和强度的预报，尤其是0~24 h预报的改善效果较好。

(2) 同化试验的分析场与控制试验的预报场相比，各物理量在不同层次上与观测的偏差和均方根误差均有不同程度的减小，风场和相对湿度在对流层中高层减小幅度较大；温度整层都有改善，在低层和高层改善更显著。

(3) 通过对比控制试验预报场和同化试验分析场各变量在对流层低层(850 hPa)和中层(500 hPa)的水平分布可见，EnKF方法同化探空资料对温度、湿度、风场和位势高度场在这两个层次均有大幅度的调整，调整后的DA分析场比NODA预报场更接近实际观测，有效改进了模式初值。

综上，WRF-EnKF系统同化探空资料后，有效减小了模式初始条件与实况的偏差，使预报的气象要素场更加接近实况，24 h定量降水预报在强度和落区上均较未同化资料的数值试验与实况更相符。但是，同化试验的36 h时效的24 h累积降水落区预报并不理想(图略)，可见WRF-EnKF同化系统对于QPF的改善具有一定的时效性。但是，由于本研究仅针对一个降水个例进行试验，检验结果不具有广泛的代表性，还需要进行更多的样本试验来对比验证这一系统的性能。此外，本研究中的个例试验同化的资料只有探空观测，并且分析所用的网格分辨率(13.5 km)对于中小尺度系统而言仍然较粗，因此，今后的工作要针对更多的降水个例进行试验，并提高模式分辨率，同化更多的观测资料(如地面观测、雷达径向风资料等)，进一步检验该系统对不同天气背景下降水的预报能力。