引言

随着天气雷达技术的进一步发展, 强对流天气短时临近预报技术的研究取得了长足进步^[1]。目前使用的短时临近预报方法主要有以雷达资料为基础的外推预报技术^[2]、基于数值模式的短时临近预报技术以及概念模型预报技术。以雷达资料为基础的外推预报技术中，交叉相关追踪法和回波特征追踪法是现在发展比较成熟的技术^[3]，已应用于许多短时临近预报系统中。但是这两种方法时效比较短，不能预报出对流系统的增强和减弱趋势。随着精细数值天气预报技术和计算机技术的发展, 利用多普勒天气雷达资料和其他中小尺度观测资料进行数值模式初始化, 来预报雷暴的发生、发展和消亡已经成为一个研究的热点，但该技术发展还不够成熟。通过分析和借鉴目前使用的短时临近预报方法^{[4-8, 12-16]}，发现仅依赖于雷达外推预报而不经过强度和形状的修正，对于未来变化趋势明显的对流系统，其预报误差将随时间快速增大。对于中尺度数值预报的降水产品来说，数值模式算法虽然还不能准确预报每个时间节点的降水量，但还是能够成功地做出整个预报周期内的降水频率和降水趋势预报。

本文尝试利用雷达资料外推，配合中尺度数值模式产品作为环境场，充分利用两者的优点，提出动态权重融合法和趋势演变叠加法这两种方法将雷达资料与数值模式产品融合进行强对流天气短时临近预报，提高强对流天气短时临近预报的准确率。

1 融合资料及资料预处理

融合方法采用了新一代天气雷达的风暴跟踪信息图形产品(STI)进行雷达外推预报，3 km的中尺度数值模式降水产品作为系统发展演变趋势预报的背景场，上海、浙江、江苏的自动站雨量作为实况对比资料。

1.1 雷达外推 1.1.1 风暴跟踪信息图形产品(STI)

风暴跟踪信息图形产品(STI)是风暴单体识别和跟踪(SCIT)算法的图形方式输出产品，SCIT算法由四个部分组成：风暴段识别、风暴中心定位、风暴单体跟踪和风暴位置预报。应用该产品时首先对1小时降水估测产品OHP和回波强度进行资料解码，提取文件头中实际需要的信息，并根据预报需要读取回波强度值或者经过订正的降水量估测值。之后通过滤波功能块，将弱降水数据滤除，获取边界清晰的对流回波区域图。在滤波功能块中提供两种滤波方法，一是固定过滤法，即选择过滤阈值(30 dBz或35 dBz)进行过滤。另一种方法是搜索最大降水量值，并要求相邻达3个像素点，根据最大降水量值减去一个阈值(可选)，获得最佳滤波效果，这种方法较适合灾害天气的早期识别。线性外推回波的位置是直接利用STI产品中对应的回波移动速度和方位对选定的回波区域进行位置的配对和外推预报。当新单体生成没有历史记录时，采用所有已被识别单体的平均移动。如果没有其他单体被确定，则采用UCP中预设的缺省速度和方向^[9]。

1.1.2 雷达外推预报

由于对较小的单体，尤其相互距离较近的多个单体，交叉相关法识别和追踪的效果较差，而单体质心法能较好地识别和追踪较小的孤立单体^[10]，因此在此直接利用STI产品中对应的雷暴移动速度和方位数值，对选定回波区域的反射率因子进行外推，再根据Z-R关系得到降水率，对降水率进行累加，得到1小时雨量预报。研究表明，上海地区对流性降水的Z-R关系为Z=200R^1.46[11]。当没有单体被确定时，外推算法停止运行。当新单体生成时没有历史记录，采用所有已被识别单体的平均移动。没有实时的STI文件时，可用前一时次的STI文件，向前沿用一个过程，即直到没有单体被确定生成时。该过程一直没有STI文件时，即没有其他单体被确定，则采用预设的缺省速度和方向。上海地区风暴单体移动的缺省值为移向260°，移速40 km·h^-1。

1.2 数值模式产品

使用的数值模式产品是上海台风所的快速同化更新系统STI-WARR的输出结果。WARR于2008年底实现业务化，主要针对短时临近预报。该预报系统逢整点同化一次观测资料，然后进行预报，当起报时次为02 UTC、08 UTC、14 UTC、20 UTC时预报时效为12小时，其余起报时次的预报时效为6小时。WARR在每日的14 UTC初始启动一次，初猜场为GFS 08 UTC起报的6小时预报场，其余时次的初猜场为WARR系统的1小时预报场，经过ADAS-3DVar同化观测数据后得到模式初始场，边界场都由GFS 08 UTC预报场给出，每3小时更新一次。自动站资料在经过时间和空间一致性检验后，也进入了STI-WARR的同化系统，而且模式使用了DFI(Digital Filter Initialization)技术以滤除频繁同化资料容易导致的高频扰动。目前，DFI往前往后各积分了20分钟。该系统在IBM/AIX系统下运行，每次预报滞后70分钟左右出结果。

1.3 自动站雨量资料

选取上海、浙江、江苏自动站在研究范围内(28.7°~33.34°N、119.585°~122°E)的自动站，将分钟雨量资料累加为1小时累积雨量，并使用Cressman插值方法将其统一插值为1 km×1 km分辨率，和雷达外推产品分辨率一致。Cressman插值方法具有较大的灵活性，可用于不同要素的分析。其误差订正原理独具特色，主要根据实测站点资料，对分析场进行连续订正，适当的平滑处理过程也包含其中。

2 融合技术及方法实现

提出两种实现雷达外推与数值模式融合的方法：动态权重融合法和趋势演变叠加法。融合方法既考虑了雷达外推预报在预报降水系统位置方面的准确性，又考虑了数值模式预报反映降水系统变化的能力，避免了单个方法的局限性。

动态权重融合法：雷达外推预测值和数值模式的输出预报值的相对权重随着时间的改变需要调整，在较短的预报时间内，雷达外推方法预测值取最大权重，但是随着预报时效的延长，雷达外推预报的误差增加，因此在预报时段较长时给数值预报结果一个比较大的权重以延长预报的时效和准确性。使用三种方法：正弦权重、双曲正切权重、实时滚动权重计算数值模式预报的权重，最后，根据得到的权重，将雷达外推与数值模式预报融合。

趋势演变叠加法：计算每个像素点对应的外推时间内数值模式预报的对流系统增强减弱的定量化趋势，并与雷达外推叠加，使得预报结果能够更好地反映对流系统的发展趋势。

2.1 动态权重融合法 2.1.1 正弦权重

不同的对流系统具有不同的生命史，大量的观测资料研究表明, 对流风暴单体的生命史一般小于半小时, 而多单体和超级单体风暴包括飑线等生命史则大于半小时，研究发现，对于局地的强对流系统，半小时内雷达外推的效果比较好，可以直接作为预报值使用，而数值模式在两小时以上预报误差相对较小^[1]。因此在前半小时内假设雷达外推的权重为1，数值模式的权重为0，两小时后雷达外推权重为0，数值模式权重为1，即前半小时全部用雷达外推；两小时后完全使用数值模式的预报，0.5小时到2小时之间雷达外推的权重W_r(t)从1逐渐下降到0；而数值预报的权重W_m(t)则从0增加到1。W_r(t)与W_m(t)必须满足：

$ \left\{ {\begin{array}{*{20}{c}} {{W_m}(0.5) = 0}\\ {{W_r}(0.5) = 1}\\ {{W_r}(2) = 0}\\ {{W_m}(2) = 1} \end{array}} \right. $

取适合此条件的函数曲线sin²(at+b)和cos²(at+b)分别作为数值模式以及雷达外推的权重曲线。根据以上列出的条件，有：

$ \left\{\begin{aligned} {\sin ^2}(0.5{a} + {b}) = 0\\ {\sin ^2}(2{a} + {b}) = 1 \end{aligned}\right. $

则雷达外推和数值模式的权重W_r(t)与W_m(t)分别为：

$ {W_m}({t}) = {\sin ^2}( - \frac{\pi }{3}t + \frac{{7\pi }}{6})\\ {W_r}({t}) = {\cos ^2}( - \frac{\pi }{3}t + \frac{{7\pi }}{6}) $

W_m(t)的权重曲线如图 1所示。

从权重曲线看出，用sin²(at+b)作为数值模式的权重曲线，在融合时间段的前段和后段时间里权重随时间的增加上升得比较缓慢，在中间段时间数值模式权重增长很快。

2.1.2 双曲正切权重

借鉴香港RAPIDS暴雨预报系统，在计算数值模式权重时取经验方程：

$ {{W}_m}({t}) = \alpha + (\frac{{\beta - \alpha }}{2}) \times \{ 1 + tan{h}{\rm{[}}\gamma {\rm{(}}{t}{\rm{-3)]}}\} $

式中，α和β分别是t=0及t=6的权重，γ代表在融合时段中间部分W_m(t)的斜率，根据该方程，前半小时的预报主要基于雷达外推预报的结果，而在预报的后段，数值模式的比重相应地增加。为了令W_m(t)于2小时至4小时的变化较为平滑，γ值被设定为等于1。

根据图 1b画出的权重曲线可以看出，在融合前段时间数值模式权重上升较缓慢，融合时段中间部分上升比较快，斜率为1，在融合时段最后部分数值模式权重增加率又趋于平缓并接近于1。

最后，根据用3种方法得到的权重，将雷达外推与数值模式预报根据以下公式进行融合，得到融合后的降水预报：

$ R({t}) = {W_r}({t}){R_r}({t}) + {W_m}({t}){R_m}({t}) $

式中, R(t)为融合后的雨量，R_r(t)为雷达外推预报得到的雨量，R_m(t)为数值模式预报的雨量。

2.1.3 实时滚动权重

使用预报时次前1小时数值模式预报误差滚动计算数值模式权重，分别用数值模式预报结果以及雷达外推预报结果与自动站实况之差，计算数值模式和雷达外推在预报时次前1小时的误差：

$ \Delta {R_r} = \left| {{{R}_r}({t - 1}) - {{R}_z}({t - 1})} \right|\\ \Delta {R_m} = \left| {{{R}_m}({t - 1}) - {{R}_z}({t - 1})} \right| $

式中，ΔR_r和ΔR_m分别为雷达外推以及数值模式的误差，R_r(t-1) 和R_m(t-1) 分别为雷达外推和数值模式在预报时次前1小时(t-1)时刻的预报值，R_z(t-1) 为自动站雨量计资料在(t-1) 时刻的实测值。

根据预报时次前1小时的误差计算出数值模式的权重W_m(t)：

$ {W_m}({t}) = \frac{{\frac{1}{{\Delta {R_m}}}}}{{\frac{1}{{\Delta {R_m}}} + \frac{1}{{\Delta {R_r}}}}} $

使用预报时次前1小时数值模式预报误差动态调整数值模式权重对于对流系统增强减弱趋势的修正更为准确。

2.2 趋势演变叠加法

在一定区域内，找到每个格点对应的外推时间内数值模式预报雨量的发展趋势，将其定量化，再将这种定量化趋势叠加到雷达外推预报结果中。在选定区域内，使用实际t时刻起报的数值模式预报对(t+Δt)和(t+Δt-1) 预测时刻的数值预报，可以得到每个像素点对应的外推时间内数值模式的发展趋势：

$ {R_t}({t} + \Delta {t}) = {R_m}({t} + \Delta {t}) - {R_m}({t} + \Delta {t - 1}) $

式中，R_t(t+Δt)为每个格点对应的数值模式的发展趋势，R_m(t+Δt)为数值模式预报的(t+Δt)时刻的雨量，R_m(t+Δt-1) 为数值模式预报的(t+Δt-1) 时刻的雨量。

再将数值模式预报的发展趋势与已插值处理后的雷达外推预报叠加，得到融合值：

$ R(t + \Delta t) = {R_r}(t + \Delta t) + {R_t}(t + \Delta t) $

式中，R(t+Δt)为融合后(t+Δt)时刻的雨量，R_r(t+Δt)为雷达外推雨量。

与风暴单体路径外推预报相比，这样得到的产品虽然对强度较弱的风暴单体发展趋势不能准确预报，弱降水区变形严重，但对于强降水的强度和位置预报更加准确。因此，在实际应用中，在强对流天气条件下也具有一定的参考价值。

3 个例分析与误差检验 3.1 个例分析

对6月21日09时至13时一次降水过程中雷达回波变化及数值模式1小时降水量预报结果进行比较。从图 2中看出，数值预报较准确地预报出了对流系统随时间增强减弱的趋势。09时，强回波区主要位于浙江北部及上海北部；10时雨区范围扩大，09时位于浙江北部及上海北部的回波连成一片，且在南通附近回波增强；数值预报结果预报出了这个增强的趋势。11时至13时，位于上海西北方的强回波逐渐减弱，且在南通以东回波加强，而位于上海及上海以西地区的雷达回波则逐渐减弱，这种加强与减弱的趋势也能从数值模式的预报结果中体现出来。

根据雷达探测的范围以及上海所处的地理位置，选取上海附近区域(28.7°~33.34°N、119.585°~122°E)作为研究区域，分辨率为1 km×1 km，如不特殊说明，文中所用时间均为世界时，降水量单位为mm。为了更直观地对各种方法进行对比，将各种方法预报出的降水量均按照统一的降水分级画出。

下面以2009年6月21日10时及12时的1小时降水量预报为例分析雷达外推预报、数值模式预报以及将它们融合后的预报效果。

图 3a和图 3b显示了2009年6月21日10时实况1小时雨量，从图中可以看出，10时主要雨区在上海西北部，太湖东南以及江苏南部地区。下面将雷达外推预报，数值模式预报以及融合后的1小时雨量预报与实况进行对比。

将2009年6月21日上午09时外推到10时的雷达外推预报结果(图 3c)与数值模式预报结果(图 3d)和10时实况的1小时降水OHP(图 3a)以及自动站1小时雨量(图 3b)进行对比可以看出，雷达外推预报的降水区域位置基本准确，但是区域大小以及对降雨量的估计不是很准确，没有预报出雨区的形变以及雨量大小的变化而数值模式对雨量的估计普遍偏大，且对对流单体中心位置的预报与实况有偏差。根据图 4可见，融合了雷达外推和数值模式的预报均修正了雷达外推预报的雨区，并且对雨量也进行了修正，反映出对流系统增强减弱的趋势。其中，正弦权重修正后的结果(图 4a)保留了雷达外推预报的位置，并且因为融合了数值模式的预报，也加入了对雨区形状及雨量大小的修正，改进效果比较明显。双曲正切权重(图 4b)以及实时滚动权重(图 4c)体现了雨量最大区域且位置准确，但后者预报的雨量量级普遍偏大。经过趋势演变叠加法融合后预报的雨区范围比前几种方法有较大的改进(图 4d)。

图 5a为2009年6月21日12时实况的雷达1小时雨量产品，这时主要的雨区已移至长江以北，10时在上海西南方的雨区已明显减弱。图 5c中，雷达外推对整体雨区的预报较为准确，但是并没有预报出上海西南方雨区减弱以及长江以北的雨量增强的趋势，而数值模式则预报出长江以北雨量的增强趋势，但对雨区位置的预报并不十分准确(图 5d)。图 6显示经过融合后，在雷达外推预报的基础上融合了数值模式预报的雨量增强趋势，修正了雷达外推预报。其中，用正弦权重(图 6a)以及双曲正切权重(图 6b)实现融合，既将数值模式预报出的雨量增强趋势加入了雷达外推预报中，并且对于雷达外推预报的雨区扩大范围较小，与其他两种方法相比，更接近实况。

3.2 误差检验

使用RMSF检验方法对这几种融合方法进行检验。

RMSF定义为：

$ RMSF = \exp {\left\{ {\frac{1}{N}\sum\limits_{i = 1}^N {\ln {{(\frac{{{F_i}}}{{{G_i}}})}^2}} } \right\}^{\frac{1}{2}}} $

其中F_i为预报值，G_i为分析值。为了避免赋值于微小预报，只计算F_i或G_i在一定范围内的量。即α＜F_i＜β或α＜G_i＜β，其中对于降水累积，α=0.2，β未赋值。

RMSF值越高，说明预报准确率越差。

表 1列出了雷达外推预报，数值模式预报以及用两种方法将雷达外推与数值模式进行融合的预报误差。

根据表 1列出的结果可以看出，数值模式准确性最差，动态权重融合法中的正弦权重以及双曲正切权重对两个时次的预报与实况误差较小，明显优于雷达外推和数值模式预报。动态权重融合法中的实时滚动权重以及趋势演变叠加法在对两个时次的预报中与实况误差都大于雷达外推，需要做进一步的改进。

4 结论与讨论

在分析讨论了目前常用的两种强对流短时临近预报方法：雷达外推和数值模式预报的基础上，参考了国内外最新发展的短时临近预报系统，提出两种将雷达外推预报与数值模式预报融合的方法：动态权重融合法以及趋势演变叠加法。在动态权重融合法提出3种方案计算权重：正弦权重、双曲正切权重以及实时滚动权重。

(1) 通过分析个例以及误差检验证明利用动态权重融合法和趋势演变叠加法进行融合后的预报结果对雷达外推和数值模式预报均有所改进，动态权重融合法中的正弦权重和双曲正切权重改进效果最为显著，对今后的短时临近预报技术发展研究有一定的参考价值。

(2) 融合后的预报结果与雷达外推和数值模式预报相比，预报雨区的范围更接近实况，改进了数值模式预报雨强普遍偏大的问题，修正了雨强中心位置的偏差。

(3) 在分析融合结果中发现，融合后的降水预报扩大了雨区，而且融合时使用固定的数值模式权重曲线也会使融合的误差增大，因此要进一步优化和调整雷达资料和数值模式融合技术，改进预报效果。

(4) 在动态权重融合法中对于数值模式权重的选取，需要统计更多的强对流个例来优化和改善，让雨量融合方案借助经验积累而得以更具代表性和更能够反映雨区的变化趋势。在以后的研究中也可以对降水类型及降水过程的不同阶段进行分类讨论。