引言

已有实践表明，确定性的临近外推预报对于1~2 h的临近预报有较好的预报能力，利用现代化的高分辨率实况探测手段，结合相对简单的外推模型就能够取得较为满意的高时空分辨率的降水预报结果。而当时效增加至1~5 h，基于实况的直接外推无法预测天气系统的生消演变，此时结合高分辨率数值预报产品对中小尺度天气的发展进行分析判断，将外推与数值预报相互融合，是将预报时效提高至2 h以上的直接途径(郑永光等，2010；Sun et al, 2014；杨璐等，2019)。

当前外推算法主要是利用前期回波的连续观测计算得到其移动矢量场用于之后的预报，常见的有TITAN法(thunderstorm identification, tracking, analysis, and nowcasting; Dixon and Wiener, 1993)、TREC法(tracking radar echoes by correlation; Rinehart and Garvey, 1978；Wilson et al, 1998);光流法(Gibson, 1979；Yilmaz et al, 2006)等。传统的单站雷达短时外推预报时效往往在1~2 h，若要用于2 h以上的融合预报则需对其进行扩展，如奥地利气象局的INCA(Integrated Nowcasting through Comprehensive Analysis)系统中通过在外推矢量计算中加入数值预报风场约束进行6 h时效的外推预报(Haiden et al, 2011)。而对于区域高分辨率数值模式结果，当其直接应用于短时临近预报时，需要考虑几个方面问题：一是时效性，数值预报计算相对消耗时间，因而往往其积分的前几小时结果不具备实际使用时效，而当时效增加时，其预报误差难免产生增长，因而常采用快速循环更新同化模式作为短时预报参考(陈敏等，2010；陈葆德等，2013)；二是分辨率，对于当前常用区域数值模式，其业务产品时间分辨率一般为1 h，空间分辨率在3~10 km，虽已足够反映中小尺度天气系统，但相比于实际雷达外推(10 min/1 km)的分辨率仍显得粗糙，要将其与外推预报进行融合得到平稳变化的结果，对数值预报产品进行降尺度及偏差订正显得十分必要。

因此，对于短时定量降水预报，如何将常规的临近外推预报与短时数值预报相互订正并融合是融合算法的关键问题。英国的“猎人”(Nimrod)系统(Golding, 1998)通过给予外推与数值预报不同时效上不同权重来进行短时临近强降水预报；美国国家大气研究中心的“尼沃特”(Niwot)系统(Wilson et al, 2010)则通过数值预报调整外推回波范围；我国香港天文台的SWIRLS(Short-Range Warning of Intense Rainstorms in Localized Systems)系统中，通过位相修正和强度修正技术订正数值预报，并由一个双曲函数确定外推与数值预报的融合权重因子，有效提高了0~6 h的降水预报效果(Li et al, 2000)。由于SWIRLS系统相对较为成熟，并有较好的开放性，北京、湖北等一些省级气象部门也在其基础上开展了相关的本地化工作(程丛兰等，2013；薛峰等，2017；崔春光等，2021)。然而，前期融合预报一般直接采用降水这一累积预报量进行，并且订正多针对外推或数值预报某一方面进行，程丛兰等(2019)指出，对于短时临近融合预报方法，将结合数值预报改进传统的雷达外推和结合雷达实况改进数值模式预报两者相结合是未来的发展方向。本研究在主流融合框架基础上，基于长三角地区雷达观测资料及高分辨率数值模式预报，对雷达外推及模式预报进行订正融合，建立了一种精细化的长三角地区短时临近融合定量降水预报技术。

1 预报融合技术

本文实现了一种长三角地区的0~6 h短时临近融合定量降水预报技术，空间分辨率为3 km，时间分辨率为10 min，采用与SWIRLS(Li et al, 2000)及程丛兰等(2013)相似的融合预报总体框架，但选用组合反射率进行前期的订正及融合，采用时间降尺度、目标识别订正、分时分类权重系数融合等技术，对降水回波进行直观的订正，得到精细化预报产品，技术路线如图 1：(1)首先对实况雷达资料进行组网扩展及时效延长，得到预报时效6 h的外推预报产品；(2)然后对数值模式预报反射率通过与当前实况雷达拼图对比进行实时的强度和位相修正，得到预报时效6 h的模式订正预报产品；(3)对0~6 h的外推和模式产品进行分时效的权重融合，在临近时效内更多地反映外推预报信息，随时效增加逐渐提高数值预报部分所占权重；(4)最后通过实时分类Z-R关系选择转换得到定量降水预报产品。

1.1 基于0~2 h雷达外推延长的0~6 h雷达回波外推预报

当前常规业务单站雷达外推产品时效为2 h，为将外推时效延长至6 h以便与数值预报进行融合，首先需要增大用于外推的雷达资料范围，将单站雷达资料扩展至周边区域雷达组网资料。本文利用长三角周边地区11部多普勒雷达反射率因子实时观测，对其进行综合组网质量控制及外推(表 1)。采用中国气象科学研究院研发的天气雷达三维组网软件(RADARMOSAIC，肖艳姣和刘黎平，2006；张志强等，2007)对各雷达站基数据进行质量控制及组网拼图，该软件通过将雷达观测的多仰角锥面数据转换并拼接成以经纬度、海拔高度为坐标的大范围三维网格数据，数据产品使用便捷，能够提高对降水系统监测的完整性。特别是对于回波演变较慢的系统性天气过程，通过雷达拼图的方式扩展观测数据覆盖范围进行整体外推，能够有效地增加雷达外推预报时效。

本文基于一种改进的COTREC方法(continuity of TREC；Li et al, 1995；陈雷等，2009)生成外推预报产品。该技术首先采用交叉相关法(TREC)求出TREC风场，对求出的风场先进行9点平滑处理，去掉明显失真的风(极大或为零)，再经过水平无辐散处理后得到COTREC风场。由于TREC法反演的风场只在有回波的区域有值，经处理后得到的COTREC风场也会受到一定程度的削弱，特别是孤立的块状回波或线状回波受到的削弱更加明显。采用引入数值预报平均风场作为TREC风场引导流场的方法，较好地解决了TREC风场经过水平无辐散处理所受到的削弱问题。通过将数值模式预报中低层风场(500/700/850/925 hPa)取矢量平均，并插值到TREC风场格点上，代替TREC风场中的缺值点，经过这一步骤后得到新TREC风场[nu⁰(i, j)，nv⁰(i, j)]，对新的风场进行水平无辐散限制求出COTREC风场[u(i, j)，v(i, j)]。最后，用该风场将t₂时刻的回波外推到t₂+Δt时刻，完成回波的外推(图 2)。

由于COTREC方法的平移矢量场是由当前相邻帧雷达观测资料计算得到，在短时效内可作为未来回波移动的近似预报，而当时效增加时，由于大气流场演变，COTREC方法必然丧失其有效性，因而本文在陈雷等(2009)工作基础上进一步提高数值模式风场的融合应用，当时效达到2 h以上时，逐渐使用该时次数值预报引导风场代替COTREC风场，提高对回波整体移动趋势的预报，权重系数采用一个双曲函数确定：

$ \begin{gathered} \boldsymbol{V}=W_{\mathrm{C}} \boldsymbol{V}_{\mathrm{C}}+W_{\mathrm{S}} \boldsymbol{V}_{\mathrm{S}} \\ W_{\mathrm{S}}=0.5\left[1+\tanh \left(\frac{t-t_0}{h}\right)\right] \\ W_{\mathrm{C}}=1-W_{\mathrm{S}} \end{gathered} $

式中：V表示外推方法的回波移动矢量，V_C表示COTREC方法计算移动矢量，V_S表示数值预报平均引导风场。W_C, W_S分别表示COTREC和模式引导风的权重系数，t₀表示引导风场权重0.5时的预报时效，本文中取2 h。

1.2 数值预报时空分布处理

为提高对更长时效的短时降水预报效果，需要进一步加入数值预报结果与雷达外推产品进行有效融合。数值模式预报部分基于华东区域快速更新同化系统WARR(陈葆德等，2013)的数值预报组合反射率资料订正生成，该模式水平分辨率为3 km，逐小时更新未来12 h华东地区的预报。然而，当前业务中使用的区域数值预报系统的产品输出间隔往往在1 h左右，与雷达外推产品的分辨率(10 min)存在较大差距。此外，由于数值预报的时效性相比外推预报仍有2 h左右的延迟，并且与实况观测相比，模式降水系统仍不可避免地存在着强度及位置上的误差，因此在与外推结果进行融合之前必须先对模式的反射率预报进行时空分布处理。

(1) 数值预报时间降尺度。对于降水这类在时空上不断移动变化的物理量进行时间降尺度处理，难以使用传统的插值方法，为将数值模式预报组合反射率产品时间分辨率降尺度到10 min的雷达外推产品上，本文采用Bao et al(2019)开发的基于深度感知的视频帧插值技术(depth-aware video frame interpolation, DAIN)，该方法采用机器学习算法，结合光流投影对上下层不同距离的对象进行采样，其结果在各种数据集上均取得了最优的测试结果，利用该视频帧内插技术对逐小时输出的模式预报降水产品内插值到10 min分辨率上，不仅能够对回波强度演变进行内插，同时能够很好地估计回波的中间运动及形态演变。

(2) 数值预报反射率强度修正。通过与临近时次实况观测的反射率概率分布对比，计算表征模拟与实况的分布态特征量，通过调整量级、周期等特征量，使其分布特征接近实际观测，达到强度修正的目的(图 3)。采用Weibull函数分别拟合当前观测与模式预报雷达组合反射率的频率分布曲线，拟合函数形如：

$ f(x ; \lambda, k)= \begin{cases}\frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1} \mathrm{e}^{-(x / \lambda)^k} & x \geqslant 0 \\ 0 & x <0\end{cases} $

式中：x为反射率值，k为形状参数，λ为尺度参数。

将模式预报反射率强度F修正至F_cal：

$ F_{\text {cal }}= \begin{cases}\left(\frac{O_{\max }}{F_{\max }}\right) F, & F \leqslant F_{\max } \\ \left(\frac{O_{10 p}}{F_{10 p}-O_{\max }}\right)\left(F-F_{\max }\right)+O_{\max } & F>F_{\max }\end{cases} $

式中：F_max/F_10p，O_max/O_10p分别为预报及观测的Weibull拟合曲线的峰值对应x坐标值及位于峰值右侧大小为峰值10%处位置所对应的x坐标值。

(3) 数值预报反射率位相修正。对于数值模式降水预报的落区误差判断常常采用基于目标对象的检验方法(徐同等，2019；盛春岩等，2021)，本文采用基于多尺度目标识别的技术(茅懋等，2016)，进行预报与实况回波目标识别与匹配。针对雷达组合反射率因子，首先进行设定阈值的空间连续对象单体目标识别，筛选出大于一定强度量级(20 dBz)的观测目标和预报目标，随后计算预报目标集和单个观测目标之间对应的面积、形状、位置参数评分(Score-Area，ScoreGC，ScoreShape)，通过三者加权计算得到匹配分数(茅懋等，2016)，根据匹配分数最高的预报目标与观测目标，确定预报和观测的目标匹配对应关系。

在实际应用中，对于高分辨率的雷达观测和数值预报，回波形态和分布复杂、变化较快，直接进行目标识别的结果往往较为零乱，对于本文中针对与外推预报融合的模式反射率位相修正来说，主要是调整预报和观测之间系统性的雨带、主要回波落区的位置偏差，以使短时效内数值预报能够更好地与外推进行融合，因此，在进行目标匹配之前首先对观测与预报的雷达组合反射率进行基于快速傅里叶变化的低通滤波处理，去除小尺度回波的影响，仅考虑整体雨区间的目标对象差异。如图 4所示，首先经过对低通滤波后的雷达组合反射率进行目标识别，观测上识别出两块回波目标(图 4a：O1，O2)，而预报上识别出三块目标(图 4b：F1, F2, F3)；随后经目标匹配计算得到预报目标与观测目标的对应关系(F1-O1，F3-O2，F2无匹配)；再计算匹配目标间的重心位置差异矢量V_i(图 4b：F1, F3目标红色虚线内的均一矢量场)，可见F1相对O1位置偏西，F3相对O2位置偏东北，因此需分别向东、西南方向调整；对于目标范围以外以及未匹配到观测目标的预报区域回波调整，则根据已匹配目标的矢量场进行插值得到：

$ \boldsymbol{V}(x, y)= \begin{cases}\boldsymbol{V}_i & (x, y) \in F_i, i=1, N \\ \frac{\sum\limits_{i=1}^N \boldsymbol{V}_i S_i / D_i}{\sum\limits_{i=1}^N S_i / D_i} & (x, y) \notin F_i, i=1, N\end{cases} $

式中：V(x, y)为(x, y)处的位相平移调整矢量，V_i为处于目标匹配对象F_i范围内的调整矢量，S_i为目标对象F_i的面积，D_i为(x, y)距离对象F_i重心的距离，因此，在进行目标范围外的平移矢量插值时，综合考虑了距离不同目标的距离以及目标的尺度，距离目标越近、目标尺度越大则越接近该目标的平移矢量(图 4b，F1，F3目标红色虚线外的矢量场)。对模式原始预报应用矢量场V(x, y)进行平移，得到调整后的预报场(图 4c)。

1.3 融合降水预报

在将雷达外推与订正后的数值预报进行融合时，采用不同的权重系数，使得前期(0~3 h)以雷达外推占比为主，后期(4~6 h)以数值预报为主。权重系数参考香港天文台的经验方程，同样采用一个双曲函数确定：

$ \begin{gathered} W_{\mathrm{NWP}}=0.5\left[1+\tanh \left(\frac{t-t_0}{h}\right)\right] \\ W_{\mathrm{EX}}=1-W_{\mathrm{NWP}} \end{gathered} $

式中：t₀表示模式预报权重0.5时的预报时效。权重系数参数按照实际使用进行调整，对于尺度小、变化快的系统t₀一般取为3 h，而对于较大尺度的系统其生命史也较长，本文采用根据系统尺度调整t₀取值的方式提高外推部分权重，前期检验发现，外推预报效果低于数值预报的节点往往出现在预报时效3~4 h区间，因而本文以当前实况回波20 dBz 以上面积作为判据，取4×10⁴(9×10⁴) km²为下(上)阈值，回波尺度位于阈值区间内时对应t₀取值在3~4 h线性变化，而对小于(大于)下(上)阈值的t₀取为3(4) h。

此外，在进行实际融合的过程中，时常仍会遇到外推及模式预报位置不一致的情况(如图 5b, 5c中江苏浙江交界处回波预报，而实况图 5a中回波基本位于两者之间)，特别是对于融合时效3~4 h的结果，此时外推及模式的权重系数一般都在0.5左右，由于两个融合成员的位置偏差造成该时效上融合结果明显偏弱(如图 5d，同时弱于图 5b外推及图 5c模式预报)。因此对于模式在外推预报晴空区的降水，则需适当降低t₀取值(本文中取1.5 h)提高模式预报部分的权重，使融合结果在中后期时效的强度变化更为合理(图 5e)。

最后，将融合后的雷达反射率预报通过Z-R关系转换为降水率，采用预报区域前3 h的地面自动站降水与对应的雷达回波作为统计样本，结合上海中心气象台雷达降水分类Z-R关系历史分类结果(表 2)，选取样本平均绝对误差最小的Z-R关系，供当前融合预报采用(图 6)。

2 预报个例及2021年梅雨期预报检验

本部分通过两个预报个例试验说明融合预报方法效果，预报个例分别选取一次天气尺度的台风影响过程以及一次天气尺度雨带与中尺度局地对流共存的梅雨期降水过程，之后对2021年长三角梅雨期间批量融合预报试验结果进行定量检验。

2.1 2019年8月10日台风利奇马过程

2019年第9号超强台风利奇马于8月10日凌晨在浙江省温岭市登陆，10日白天经江苏省苏州市北上，受其影响，江苏、浙江一带普降暴雨到大暴雨，局地产生特大暴雨。从8月10日13时起报的融合预报来看，时效2 h后的外推与模式预报出现了不同的误差，实况40 dBz以上的强回波区主要位于台风中心北侧的江苏及上海北部一带(图 7a)，外推预报对于台风及降水雨带的整体位置预报较好，但对于安徽东南部台风西北象限螺旋雨带的回波预报偏强(图 7b)，数值模式预报的原始反射率强度分布与实际较为接近，因而强度调整前后基本没有变化(图 7c，7d)，主要是通过目标识别匹配对模式台风的位置进行订正，原始数值模式预报台风整体位置较实况偏南，经过初始时刻目标位置匹配订正后，模式未来预报台风整体北调了约80 km，与实况更接近(图 7e)，通过与外推预报融合之后，虽然整体回波范围有所扩大，但由于外推及模式订正结果对于台风中心北侧35 dBz以上的强回波带均有体现，融合产品对于此处的强回波预报带得到了保留，而对于其西侧安徽东北部的外推预报过强的回波预报由于外推与模式存在分歧，融合进行了一定程度的减弱，整体作出了较为准确的调整(图 7f)。

图 8给出此次融合预报14—15时的小时累计降水对比，其中观测由国家气象信息中心研发的中国多源降水融合系统(CMPA)的逐小时降水实况分析产品提供(图 8a)。外推预报对于降水落区范围整体预报偏小(图 8b)，而对于江苏南部到上海北部的降水预报偏强，数值模式对此次过程降水整体预报偏强，但经位相调整后，在上海北部的50 mm以上降水落区预报与外推预报较为一致(图 8c)，经两者融合后，0.1 mm以上降水整体落区范围、位于上海北部50 mm以上的强降水预报均与实况更为一致(图 8d)。

2.2 2021年7月5日梅雨期局地强对流过程

2021年7月5日，受中低层切变线及副热带高压边缘影响，长三角地区出现午后局地强对流天气。5日14时之前，雨带主要位于长三角北部的安徽、江苏北部地区，降水强度不大，观测雷达反射率基本都在35 dBz以下(图略)，14时之后，随着地面气温上升，江苏、浙江交界一带开始有分散性局地对流出现，16时雷达观测图上，北部天气尺度的梅雨锋雨带与南部的新生局地对流共同存在(图 9a)，对于此类过程，外推预报及数值模式具有不同的预报能力，从16时起报的融合预报初始时刻各部分对比来看，COTREC外推预报(图 9b)与观测回波(图 9a)一致，而模式预报由于产品到达时间滞后(融合预报初始时刻的数值模式产品起报时间一般为2 h前)存在初始偏差，由图 9c可见，模式预报上此时也存在北侧梅雨锋及江苏南部的局地对流两部分降水回波，但北侧回波带已开始减弱消散，而南侧局地对流部分也与实况存在强度及位置上的偏差。经过回波强度订正之后，模式预报的北侧回波带得到了较明显的增强(图 9d)，特别是原本强度20 dBz以上的分散回波块经订正后连接呈带状，与观测回波形态较为一致，由于本方法中采用20 dBz作为目标识别的判断阈值，此时再对经过强度调整的模式预报进行与观测回波的低通滤波后的目标识别匹配订正，能够较好地对其偏东北方向的预报误差做出判断并调整(图 9e)。对于初始时刻的融合预报结果(图 9f)，在雷达观测资料覆盖范围内采用雷达观测，而在雷达观测范围之外利用目标识别订正后的模式预报结果进行补充，可以看到对位于图 9f东北角海上的回波，融合方法较平滑的将该处的模式预报与西侧陆地的外推回波进行了连接。到5日20时，江苏南部的对流发展增强(图 10a)，外推预报由于缺少对于回波发展演变的预报能力，此时回波与初始时刻在形态强度上均没有明显变化，仅向东南方向移入浙江(图 10b)，而数值模式对于该处对流系统的发展作出了较好的预报，但强度上略偏弱，位置也略偏东北(图 10c)，经过强度调整之后，数值预报回波的强度得到了改进，其中心部分反射率达到50 dBz以上(图 10d)，经过位相调整之后，对回波的位置向西南方向进行了偏移(图 10e)，也与实况更为接近。需要指出的是，对于南侧对流部分模式预报回波的位相调整，由于该处回波初始时刻尺度较小，未能与实况该处的回波达到目标识别匹配，其平移矢量实际是由北侧梅雨锋雨带位置偏移矢量插值得到，在此个例中取得了较好的结果。此外，数值模式预报北侧梅雨锋雨带回波消散速度过快，20时回波已基本消散，而融合结果保留了部分外推预报的信息，对此处回波作出了一定程度的补充(图 10f)。

2.3 2021年6月10日至7月10日长三角梅雨期间融合预报检验

为客观说明融合预报效果，对2021年6月10日至7月10日长三角梅雨期间的逐1 h起报更新的融合预报进行批量预报检验(共计744次起报的融合反射率预报结果)。2021年江淮梅雨期间中高纬环流以经向型为主，副热带高压不稳定，雨带在长江一带南北摆动，期间既有受梅雨锋影响的天气尺度降水过程，也有局地发展的中小尺度对流过程。本节采用与实况区域组网组合反射率观测20 dBz阈值的临界成功指数评分(CSI)对融合预报各成员进行客观检验。由图 11可见，各方法CSI评分总体随时效增加而降低，其中COTREC外推预报在预报时效1 h内评分最高，但随时间下降最显著，而数值模式预报效果在6 h时效内整体下降较平缓。经过组网扩展改进的外推预报在3 h时效内评分均高于模式预报，到3 h之后模式预报逐渐开始优于外推预报。值得注意的是，对于数值预报的订正检验上看，整体强度订正和位置订正均有正效果，但模式强度订正(NWP_AMP)相对于原始模式(NWP)的改进相比模式位置订正(NWP_OBJ)相对模式强度订正(NWP_AMP)的改进更为显著。除在1 h时效内评分略低于外推预报外，融合预报产品在2~6 h的预报效果基本都优于其他各产品。

为进一步评估算法对定量降水的预报效果，采用CMPA逐小时降水实况分析产品对融合预报的逐小时累计降水进行分阈值的TS评分检验。由图 12可见，对于小雨量级(0.1 mm以上)的降水，除最初1 h时效上略低于外推外，融合预报在其他时次的TS评分均为最高，特别是在4~6 h时效上较明显地优于原外推及数值模式预报，说明融合方法对雨带的整体落区作出了较好的预报(图 12a)。而对于较大量级(5 mm以上)的降水，融合预报在1 h、3 h时效上略低于外推，而在其他时效上整体优于外推及数值模式预报，说明融合预报对于0~6 h时效内较大量级的降水预报整体也优于单一的外推或数值模式预报(图 12b)。

3 结论及讨论

本文采用当前主流的外推和数值预报融合方案，结合区域雷达外推、数值预报、目标识别订正、机器学习降尺度等技术进行应用，完成了时效6 h、空间/时间分辨率为3 km/10 min的长三角地区短时临近融合定量降水技术及业务试运行。主要结论如下。

(1) 通过对长三角地区11部雷达反射率的组网拼图和COTREC外推，同时结合在时效增加后逐渐使用数值预报引导风替代外推矢量场的方法，将原COTREC外推预报2 h的预报时效延长至6 h。

(2) 通过对数值预报反射率进行时空分布处理以使其能与外推预报进行融合。采用视频帧插值技术对数值预报1 h分辨率雷达反射率产品进行时间降尺度插值，避免了传统插值方法对于空间上移动变化的物体难以处理的困难。采用基于韦伯变换的强度误差修正和基于目标识别匹配的位相修正技术，根据实况雷达观测先调整预报回波整体强度，再通过选定阈值的目标位置误差计算调整回波空间分布，使数值预报回波的时间分辨率、强度和位置与实况外推更为接近。

(3) 采用分时效的权重融合系数对外推和数值预报进行融合，使短时时效内以外推结果为主，随后逐渐变为以数值预报为主。同时，在融合时需结合不同的系统尺度、外推和数值预报差异情况调整融合系数，最后根据实时Z-R关系判断将融合反射率转换得到定量降水预报。

(4) 通过两个不同尺度和类型的个例预报和2021年梅雨期间批量融合预报试验对融合方法进行了分析检验，结果表明本文中采用融合外推和数值模式预报的方法能较有效地结合两者优势，提高0~6 h的降水预报效果。

对于2 h以上时效的短时预报，融合外推和数值预报产品的做法是主流的共识，除了对外推和数值预报本身进行改进外[例如当前一些深度学习算法(Shi et al, 2015;Chen et al, 2020)基于历史雷达观测强度、范围、生消等信息，对回波的发展变化有一定的预报能力]，融合预报的难点在于对外推和数值预报产品的可信度进行判断。然而，当前这方面的各类判断技术往往较为线性单一，如根据模式中低层平均的引导气流计算，根据时效的权重系数变化和根据当前误差对后期预报的直接订正等，没有对天气系统的不同类型和发生发展阶段进行细致的判断区分。如何应用一些大数据处理、深度学习技术加强实况和数值预报产品的融合是短时预报的未来发展方向。