引言

近年来，我国常用的0 ~ 3 h临近预报技术主要包括基于雷达资料的雷暴识别追踪、外推预报技术和数值预报技术以及概念模型预报技术等（陈明轩等，2004）。一方面，尽管通过同化多普勒天气雷达资料等改善数值模式初始场来预报对流天气系统的生消变化已取得一定进展（Benjamin et al, 2009），但在实际预报业务应用中还不够完善，预报起转延迟(spin-up)现象在开始一段时间内仍会存在；此外，在通常的同化技术中一般是先反演风场，再反演热力场的“两步走”方式，这样会造成较大的反演误差，导致预报误差（顾建峰，2006）。另一方面，在基于雷达资料的识别追踪和外推技术中，交叉相关外推（Rinehart et al, 1978；Bjerkaas et al, 1980；Li et al,1995；Li et al, 2000；张亚萍等，2006；陈明轩等，2007；王明筠等，2010；符式红等，2012）和回波特征追踪识别外推（Dixon et al, 1993；Einfalt et al, 1990；陈明轩等，2006）应用虽然较为成熟，但仍具有预报时效较短，预报雷暴的发展演变能力较弱等缺陷（陈明轩等，2004；俞小鼎等，2012）。因此，将快速更新循环的高时空分辨率数值预报和基于雷达资料的外推预报相融合的预报技术是改善两者不足的重要方法，也是近年来国内外临近预报方法研究的重要方向（郑永光等，2010）。

目前，国内外的融合方法主要可分为三类（郑永光等，2010）：第一类是随时效分别计算外推预报和数值模式预报结果的权重系数，然后通过加权平均进行融合预报。比如由Golding（1998）研发的Nim-rod强降水预报系统，Pierce等（2000）研发的Gan-dolf临近预报系统和美国NCAR的NIWOT定量降水融合预报系统（Wilson et al,（2006）），以及Tera-da et al,（2004）、颜琼丹等（2010） 、杨丹丹等（2010）的工作中均用到类似的物理空间融合技术，融合过程中均认为数值模式预报所占权重随时效增加，而外推预报权重随时效减小。第二类是趋势融合法，比如NIWOT系统（Wilson et al,（2006）），该系统的一部分是利用模式降水落区和强度预报的趋势变化，对雷达外推预报结果进行修正，从而进行定量降水预报融合。第三类是通过计算当前时刻模式降水预报的落区或强度误差，并估计误差的时间变化趋势，利用估计的误差趋势特征，对未来相应时段的模式降水预报结果进行修正。比如，Wong 等,2009利用香港天文台研发的第二代“小涡旋”临近预报系统（Short-range Warning of Intense Rainstorms in Localized Systems，SWIRLS-2；Yeung et al,2009）中的多尺度光流变分法 (Multi-scale Optical flow by Variational Analysis，MOVA)识别模式定量降水预报的落区误差并进行相位修正，同时利用韦伯累积分布函数识别降水强度误差，根据雷达定量降雨估测值调整模式降水强度；程从兰等（2011）则借鉴香港天文台的融合技术分别应用快速傅里叶变换FFT和MOVA计算模式降水的相位偏差并进行修正，另外程从兰等（2013）还对中尺度数值模式输出的定量降水预报在谱空间进行相位校正；此外，Du-Fran 等2009利用VET方法（the Variational ra-dar Echo Tracking algorithm）计算模式预报的降水与雷达观测的相位误差，并导出误差的时间变化趋势来修正降水相位误差。

本文设计的融合方法属于上述第三类融合法，与上述第三类方法相比，其主要不同是尝试利用相关分析方法计算当前相同时刻雷达估测降水与模式反射率因子估测降水（简称为模式估测降水，下同）之间的位置偏差，并估计位置偏差的时间变化趋势，对未来相应时段的模式估测降水预报结果进行订正。文中主要探索0 ~ 2 h的模式估测降水的订正效果，并检验其能否好于利用雷达反射率因子外推的定量降水预报效果，进而为预报员发布0 ~ 2 h临近预警提供重要的科学依据。

1 资料选取及预处理 1.1 资料的选取

由于直接观测的降水在时、空分辨率上较难与数值预报降水产品匹配，所以选用与降水预报密切相关的模式预报的反射率因子(dBz)作为研究对象，便于与雷达观测反射率因子比较。雷达在3 km的探测效果最佳，故选用的雷达反射率因子和模式反射率因子的高度均为3 km。但由于反射率因子(dBz)呈对数量级增长，在订正结果检验时，不便于统计量的计算。因此本文的处理方式是将3 km高度的雷达反射率因子和模式反射率因子均转换成降水后再进行融合试验。

1.2 资料的预处理

(1) 本文以1°×1°的GFS预报场资料作为模式背景场，模式采用双重单向嵌套网格，外层格距为12 km，格点数为283×223，内层格距为4 km，格点数为423×363，内外层中心位置均为(29.88°N、107.73°E），均采用相同的垂直层，共37层，垂直分辨率为625 m。外层的输出结果作为内层的背景场，雷达资料在内层中加入，用ARPS的复杂云分析方案（Hu et al, 2006a）及ARP-3DVAR系统（Hu et al, 2006b）整点单次同化重庆雷达资料（CINRAD/SA）的反射率和径向风数据更新初始场，将同化后的初始场应用中尺度数值模式WRF进行数值预报。预报结果输出间隔为12 min，预报时效为3 h，水平分辨率为4 km×4 km。将3 km高度的数值预报反射率因子数据，通过Z-R关系转换为12 min的降水数据作为模式估测降水。本文所用的Z-R关系（Fulton et al, 1998）为：

$Z=300{{R}^{1.4}}$

(1)

式中，Z为反射率因子，单位: mm⁶·h^-1，R为雨强，单位: mm·h^-1。

(2) 雷达反射率因子拼图数据（重庆雷达、宜宾雷达，宜宾雷达型号为CINRAD/SC）和反射率因子拼图的外推预报数据由中国气象局灾害天气短时临近预报业务系统（Severe Weather Automatic Nowcast System，SWAN；胡胜等，2011；李俊等，2012）提供。产品输出的时间分辨率为6 min，水平空间分辨率均为0.01°×0.01°，外推预报数据的预报时效为3 h。将3 km高度间隔12 min的反射率因子拼图数据和外推预报数据通过Z-R关系转换为12 min降水数据，然后用双线性方法插值成水平分辨率为4 km×4 km的降水数据，下文分别简称为雷达估测降水和雷达外推定量降水预报。

2 融合方法介绍

本文设计的融合方法主要包括三部分：(1) 查找匹配网格阵列块并计算位移偏差矢量场；(2) 将矢量场分区和平滑处理，计算出各分区的平均位移偏差矢量；(3) 估计各分区的平均位移偏差矢量随时间的变化趋势特征，根据此特征，对相应时段和分区的模式估测降水位置偏差进行订正。

2.1 位移偏差矢量场的计算

参考TREC跟踪方法中划分网格阵列求相关的思路(王改利等，2007)，取3 km高度、同一时刻t的模式估测降水(图 1a)和雷达估测降水(图 1b)作为试验对象。图 1示意地给出了模式估测降水和雷达估测降水之间位移偏差矢量的计算过程。首先将模式估测降水和雷达估测降水均划分成若干个大小相同的网格阵列块，称之为阵列块，每个阵列块包含相同的网格点数，且各个阵列块几何中心有一定间隔距离。然后将其中一个模式估测降水阵列块与一定搜索半径r内的所有雷达估测降水阵列块求相关，相关系数(R)的计算公式（王改利等，2007）为：

$R=\frac{\sum\limits_{k}{{{Z}_{r}}(k){{Z}_{m}}(k)}-\frac{1}{N}\sum\limits_{k}{{{Z}_{r}}(k)}\sum\limits_{k}{{{Z}_{m}}(k)}}{{{\left[ (\sum\limits_{k}{Z_{r}^{2}(k)-N\overline{Z}_{r}^{2}})(\sum\limits_{k}{Z_{m}^{2}(k)-N\overline{Z}_{m}^{2}}) \right]}^{\frac{1}{2}}}}$

(2)

式中Z_r、Z_m分别为t时刻雷达估测降水阵列块和模式估测降水阵列块内各个网格点上的值，N是每个阵列块的格点数（k=1，2，…，N）。

将相关系数最大的雷达估测降水阵列块称之为与该模式估测降水阵列块最为匹配的阵列块。以该模式估测降水阵列块中心为起点，与之最匹配的雷达估测降水阵列中心为终点，构成模式估测降水阵列块与雷达估测降水阵列块的位移偏差矢量C。位移偏差矢量反映了模式估测降水阵列块与其匹配的雷达估测降水阵列块之间的位置偏差，位移偏差矢量的长度为位置偏差的大小，指向为位置偏差的方向。同理，求得所有模式估测降水阵列块的位移偏差矢量，构成位移偏差矢量场。

需要说明的是：在初始时刻，对于某一个模式估测降水阵列块而言，由于模式同化雷达数据，认为该阵列块与其最匹配的雷达估测降水阵列块的位置是无偏差的，但经过较短的时间 Δ t后，两个阵列块之间位置逐渐出现偏差，这种偏差具有时空连续性。由于 Δ t较短，这两个阵列块之间的距离应相距不远，为避免在较远处找到错误的最匹配阵列块，导致得到错误的位移偏差矢量，因此针对位移偏差，需指定搜索半径r来加以约束。

图 2给出了搜索半径的确定方法示意图。如图 2所示，这里做几点假设：在间隔较短时间 Δ t内，3 km 高度的模式估测降水阵列块移动由该高度上模式风场引导，移动距离≤V_模式× Δ t；同样地，3 km 高度的雷达估测降水阵列块移动由该高度上实际风场引导，移动距离≤V_实际× Δ t；3 km高度上，当模式预报的风场与实际风场差距不大时，图 2中α较小。经这样假设后，位移偏差矢量C_偏差的大小不应超过V_模式× Δ t，也不超过V_实际× Δ t。为计算方便，将搜索半径定为：

$r={{V}_{\max }}\times Δ t$

(3)

式中，V_max为模式起报时刻（即初始同化时刻）所选研究区域，3 km高度的模式预报的最大风速 Δ t为时间间隔。

2.2 分区的平均位移偏差矢量的计算

由于夏季降水往往分布不均，具有成片区的分布特点，因此，位移偏差矢量场同样具有这样的特征。为此，本文参考NCAR的Auto-nowcast system(Mueller et al, 2003)人机交互操作的处理方式，把位移偏差矢量较一致的区域划为一个分区。此外，为提高每个分区位移偏差矢量场的连续性和一致性，对矢量场进行了平滑处理。即在每个分区内，如果一个位移偏差矢量与其周围矢量的平均方向相差超过40°，则用其周围一定扫描半径内非零矢量（不少于10个）的平均值代替；然后再计算每个分区的位移偏差矢量的平均值，简称为平均位移偏差。用同样的方法，对未来连续几个时次（间隔 Δ t不宜过长）的位移偏差矢量场做分区、平滑处理，计算出各分区的平均位移偏差。

2.3 模式估测降水的订正

由于各分区的平均位移偏差随时间在变化，因此需计算平均位移偏差在未来时段方向、大小的变化趋势特征。文中采用最小二乘法对每个分区连续几个时次的平均位移偏差分别进行线性拟合，得到每个分区平均位移偏差随时间的变化趋势特征。将未来某一时刻各分区的平均位移偏差作为模式估测降水订正的依据，即模式估测降水沿着对应时刻、对应分区的平均位移偏差的方向、大小移动，由此订正模式估测降水的位置偏差，从而进行融合预报。

3 试验及检验 3.1 试验设计

为了验证上述融合方法在临近预报中的效果，选取2012和2013年夏季发生在重庆西部、四川东部的3次强降水天气过程来进行预报试验，试验参数见表 1，所用时间均为世界时。

下文以试验1为例给出了详细的融合试验过程及分析，试验2和3类似，不再赘述。

图 3给出了2012年7月21日13:12、13:24和13:36 UTC，3 km高度的模式估测降水和雷达估测降水。由图 3可见，尽管与同时刻雷达估测降水相比，模式估测降水能大致反映出降水分布的基本特点，即大致以29.5°N为界，南、北部各有一个呈带状分布的降水中心；但南部的模式估测降水中心位置偏向雷达估测降水中心位置的西南方，北部的模式估测降水中心位置在13:36 UTC也略偏向雷达估测降水位置的西南方，在空间分布上存在着一定的位置偏差。利用第二节所述方法对2012年7月21日的强降水过程进行融合试验。如表 1中试验1，根据此次过程降水分布情况，选取矩形区域(26.938° ~ 31.0°N、104.288° ~ 107.0°E)作为试验区域；计算位移偏差矢量场的时次为13:12、13:24和13:36 UTC；由于V_max约为9 m·s^-1，因此3个时次搜索半径分别定为6、12和18 km。Tuttle等(1990)研究认为，虽然计算最大相关系数的阵列块大小可以任意选择，但若阵列块太小，不能包含足够的资料点计算一个相对稳定的相关系数，若太大，则不能获得较精细的特征；而且在相关系数的计算过程中对噪声等低阈值回波非常敏感。因此，此个例以29.5°N以北的降雨中心范围大小为选取阵列块的标准，对于水平分辨率为4 km×4 km的网格，选用长方形阵列块大小11×17，即范围40 km×64 km，基本将北部的降雨中心包括在内（如图 3a红色实线方框所示），阵列块移动间隔定为1个格距，即4 km。此外，本文将计算相关系数的12 min反射率因子估测降水的最低阈值均设为0.472 mm（相当于反射率因子最低阈值为30 dBz）；为有效减少边界错误的矢量，在选取阵列块时，如果阵列块内≥最低阈值的格点数之和＜阵列块内总点数的3/5，则该阵列块被视为无效，不计算位移偏差矢量。

3.2 位移偏差矢量场

图 4a、4c和4e分别是2012年7月21日13:12、13:24和13:36 UTC未经过分区、平滑处理的位移偏差矢量场，从图中可以看出两个特征，一是某些成片的区域内位移偏差的大小、方向较一致，但与周边区域的位移偏差的方向、大小有较大差别；二是在成片的位移偏差的大小、方向较一致的区域内，有少数位移偏差与附近主要位移偏差差别较大、不连续。于是，采用2.2节所述方式进行分区、平滑处理，得到4个分区的平均位移偏差场（图 4b、4d和4f）。由图 4b、4d和4f与图 4a、4c和4e 对比可见，经过分区、平滑处理后，每个分区的平均位移偏差场与对应位置未经过分区、平滑处理的位移偏差场的整体方向保持一致，并且对每个分区内少数差别较大的位移偏差得以一定改善。由图 4b、4d和4f可见，平均位移偏差在1区和2区基本指向东北方，3区基本指向偏北方，且1 ~ 3区平均位移偏差的大小随时间逐渐增大，表明模式估测降水与雷达估测降水之间的位置偏差在平均位移偏差方向上随时间逐渐增大，4区由指向西南方逐渐转为指向偏北方，表明4区平均位移偏差随时间有指向东北方的增长趋势。对每个分区的平均位移偏差进行最小二乘线性拟合，得到每个分区平均位移偏差随时间趋势变化特征，对13:48—16:00 UTC的模式估测降水进行订正。

3.3 预报结果及分析

图 5给出了2012年7月21日13:48—16:00 UTC时段3 km高度的模式估测降水预报和融合估测降水预报结果。由图可见，在13:48—14:12 UTC时段，雷达估测降水基本稳定在重庆西部与四川接壤的区域，呈准南北向；融合前，模式估测降水位置与雷达估测降水位置较一致，但模式估测降水范围较雷达估测降水范围窄，其东西两侧有明显的“空白”区（图 5a ~ 5c），融合后，东侧边缘的“空白”区显著减少，但是西部边缘的“空白”区却扩大（图 5a₁ ~ 5c₁），可见，在这一时段的融合预报效果较模式预报略有改善但优势不显著。在14:24—14:48 UTC时段，北部（29.5°N以北，下同）的雷达估测降水位置基本稳定不动，南部（29.5°N以南，下同）的雷达估测降水逐渐向东北方移动，逐渐进入重庆地区；融合前，模式估测降水位置偏向雷达估测降水位置的西南方，南部模式估测降水位置偏离较明显（图 5d ~ 5f）；融合后，与雷达估测降水位置比较，虽然北部的融合估测降水位置略微出现偏东现象，但总体略胜于模式估测降水预报，而南部的融合估测降水位置与雷达估测降水位置十分接近，模式估测降水位置偏西南的情况得以改善（图 5d₁ ~ 5f₁），可以看出，此时段融合预报效果更优。在15:00—16:00 UTC时段，北部的雷达估测降水向东北方缓慢移动，南部仍继续向东北方移动；融合前，模式估测降水位置明显偏向雷达估测降水位置的西南方，形状有了明显不同，两者间的偏差已经非常明显（图 5g ~ 51），融合后，融合估测降水位置与雷达估测降水位置基本重合（图 5g₁ ~ 51_l），由此表明，此时段融合预报效果明显优于模式预报效果。综上所述，对于0 ~ 2 h左右的临近预报，融合预报效果总体优于模式预报效果。用同样的方法，进行试验2和试验3，所得的融合预报效果与试验1类似（图略）。

3.4 定量检验分析

为进一步检验上述融合方法对临近预报效果的改进，将对3次试验的融合估测降水预报、模式估测降水预报和雷达外推定量降水预报的预报效果进行定量对比分析。选用的检验方法是将格点上对应时刻的估测降水预报数据与雷达估测降水数据逐一比较，如果一个格点上雷达估测降水数据与预报数据都大于预先设定的阈值，则认为该格点击中，如果格点上雷达估测降水数据大于阈值而预报数据小于阈值，则认为该格点为漏报，如果格点上雷达估测降水数据小于阈值而预报数据大于阈值，则认为该格点为空报。采用以下3个指标来进行客观评价（胡胜等，2012；傅娜等，2013）：

$POD=\frac{{{N}_{\text{击中}}}}{{{N}_{\text{击中}}}+{{N}_{\text{漏报}}}}$

(4)

$FAR=\frac{{{N}_{\text{空报}}}}{{{N}_{\text{击中}}}+{{N}_{\text{空报}}}}$

(5)

$CSI=\frac{{{N}_{\text{击中}}}}{{{N}_{\text{击中}}}+{{N}_{\text{漏报}}}+{{N}_{\text{空报}}}}$

(6)

式(4) ~ (6)中，POD为命中率、FAR为空报率、CSI为临界成功指数，N_击中、N_漏报、N_空报分别为预报中击中、漏报、空报的格点数目。

图 6是3次强降水个例试验的模式估测降水预报、融合估测降水预报和雷达外推定量降水预报各自的CSI时间变化曲线及总的CSI、POD和FAR时间变化曲线。由图 6a ~ 6c可见，在144 min（约2 h）的整个订正过程中，总体而言，3次试验融合估测降水预报的CSI均高于模式估测降水预报的CSI；雷达外推定量降水预报的CSI均呈下降趋势。由试验1的CSI可见，在0 ~ 60 min订正过程中，雷达外推定量降水预报的CSI最高；60 ~ 144 min是融合估测降水预报的CSI最高，雷达外推定量降水预报的CSI高于模式估测降水预报预报的CSI（图 6a）。由试验2的CSI可见，在0 ~ 48 min(约1 h)订正过程中，雷达外推定量降水预报的CSI最高；48 ~ 144 min是融合估测降水预报的CSI最高，雷达外推定量降水预报的CSI低于模式估测降水预报预报的CSI（图 6b）。由试验3的CSI可见，在0 ~ 72 min(约1 h)订正过程中，雷达外推定量降水预报的CSI最高，72 ~ 144 min是融合估测降水预报的CSI最高（图 6c）。综上可见0 ~ 1 h，3次试验的雷达外推定量降水预报CSI最高，预报效果最优；1 ~ 2 h，融合估测降水预报CSI最高，预报效果最优，雷达外推定量降水预报效果和模式估测降水预报效果各有优劣。

如图 6d所示，在144 min（约2 h）的整个订正过程中，3次试验融合估测降水预报总的CSI高于模式估测降水预报总的CSI；雷达外推定量降水预报的总的CSI呈下降趋势，0 ~ 60 min是雷达外推定量降水预报的CSI最高，60 ~ 144 min为融合估测降水预报总的CSI最高。图 6e与图 6d比较，3次试验的3种预报方法总的POD与总的CSI趋势均相似。如图 6f所示，在144 min的整个订正过程中，3次试验的3种预报方法总的FAR整体均呈上升趋势；融合估测降水预报总的FAR低于模式估测降水预报总的FAR。0 ~ 72 min(约1 h)订正过程中，雷达外推定量降水预报总的FAR最低，72 ~ 144 min为融合估测降水预报总的FAR最低。由此可见，0 ~ 1 h，3次试验雷达外推定量降水预报总的CSI和总的POD都最高、总的FAR最低，雷达外推定量降水预报效果最优；1 ~ 2 h，融合估测降水预报总的CSI和总的POD都最高、总的FAR最低，融合估测降水预报效果最优。

综上所述，基于所分析的3次夏季强降水个例，对0 ~ 2 h左右的临近预报，融合预报效果总体上优于模式预报效果；0 ~ 1 h，雷达外推定量降水预报效果最优，雷达外推定量降水预报效果总体呈下降趋势，1 ~ 2 h，融合预报效果最优。

4 结论和讨论

本文发展了一种对多普勒天气雷达资料与中尺度数值预报产品进行相关分析的融合方法，并利用该方法对夏季发生在重庆西部、四川东部的3次强降水天气个例进行临近预报试验及检验评价，得出以下结论：

(1) 对0 ~ 2 h左右的临近预报，融合预报效果总体上优于模式预报效果。

(2) 与雷达外推定量降水预报效果相比，对0 ~ 2 h左右的临近预报，0 ~ 1 h 雷达外推定量降水预报效果优于融合预报效果；1 ~ 2 h，融合预报效果优于雷达外推定量降水预报效果，可见，在临近预报1 ~ 2 h 时效内，本文设计的融合技术对夏季强降水临近预报有一定应用价值。

另外，由于该融合方法是通过对当前几个时次各分区平均位移偏差进行线性拟合，并根据线性拟合趋势特征订正模式估测降水预报结果，因此主要应用于0 ~ 2 h临近预报中。如果订正时间较长，这种线性拟合趋势特征可能与实际情况之间存在较大误差，所以文中没有研究该融合方法应用于2 ~ 6 h的数值预报结果。今后可以尝试非线性的方式进行拟合并探讨该融合方法对2 ~ 6 h数值预报降水订正的可行性。

同时，该融合方法亦存在一些局限：首先，该融合方法依赖于模式对反射率因子的预报效果，本方法不适用于雷达观测的反射率因子和模式预报的反射率因子结构有较大差别的情况。其次，本文对位移偏差矢量场进行分区主要是采用人机交互操作的处理方式，把位移偏差矢量较一致的区域划为一个分区，在实际应用时，还需要发展对位移偏差进行自动分区的方法。另外，本文选择求取相关的阵列范围对于局地或尺度小的降水过程不一定适合，阵列大小的选择对不同尺度降水过程订正效果的影响等方面还需要进一步探讨。