引言

临近预报是指对天气未来几小时的预报。目前，很多业务上的预报系统采用数值预报模式，但由于数值预报模式有一个spin-up过程，对起始的几个小时的降水预报效果不是很理想，为了弥补数值预报模式的spin-up问题所造成的降水延迟现象的缺陷，我们这里采用多普勒雷达外推技术。

由于新一代多普勒天气雷达的探测数据具有很高的时间和空间分辨率，这些数据对预报天气系统的发生、发展具有很重要的意义，因此多普勒天气雷达已经成为改进临近预报最有效的工具之一，它具有其他气象探测设备无法取代的优势，张沛源等^[1]把新一代天气雷达应用在临近预报和灾害性天气警报中。张家国等^[2-3]利用多普勒雷达资料对2008年7月22日湖北省襄樊市出现的一次长历时的特大暴雨的雷达回波结构特征、中尺度系统活动和地形作用进行了分析，并研究了几类区域性暴雨雷达回波模型。唐仁茂等^[4]根据雷达回波自动选取对比进行人工增雨的方法，该方法实用性较强，能够快速识别出对比云，在一定程度上消除人为判别的误差，提高对流云人工增雨作业效果分析的科学性。刘晓阳等^[5]利用多种雨量计校准雷达定量降水，研究表明：校准雨量计数量和估测精度有明显正相关，校准雨量数越多，降水估测精度越高。目前，我国在很多地方已经建立了新一代天气多普勒雷达观测网，其探测范围和能力都迅速地增加。

目前，采用最普遍的外推方法有：质心跟踪法和线性交叉相关外推法。Dixon等^[6]指出质心跟踪法也另称为对流单体追踪或雷暴质心跟踪，它仅仅适用于强对流风暴的跟踪。而线性交叉相关外推算法是计算出回波最优空间相关的方法，追踪一定区域内雷达回波在过去时间段内的变化特征，通过这些变化特征来确定未来回波的移动位置和形状^[7-8]。线性交叉相关外推方法是Rinehart等^[9]于1978年首先提出的，王改利等^[10]把线性交叉相关外推算法用在暴雨临近预报中，陈明轩等^[11]利用线性交叉相关外推方法外推飑线移动过程。交叉相关外推技术也用在了追踪热带气旋以及台风风场和螺旋雨带的移动过程中^[12-13]，赵放等^[13]也采用改进的交叉相关技术追踪台风移动位置与降水的落区。NCAR的临近预报系统“Auto-Nowcast”^[8]也是采用了交叉相关外推方法构造回波运动矢量场。在这些系统中都是通过多普勒雷达回波的反射率因子来外推未来回波的移动方向与演变。

线性交叉相关算法只是考虑两个时刻的雷达回波变化，求这两个雷达回波之间的空间最优相关，来寻找回波的移动速度和移动方向，再以这样的移动矢量外推未来回波的移动情况，它是假设回波的演变过程是线性的来计算的，但是在实际情况下，回波的变化表现为非线性。为了解决这个问题，我们采用集合交叉相关外推技术，它考虑了多时刻回波的演变过程，对这些时刻的雷达回波求空间最优相关，求得几组回波在过去一段时间内的移动矢量，综合这些移动矢量，最后得出回波在未来的移动位置和形状，因此它考虑了回波演变的非线性特性，具有一定的物理意义。

1 资料和方法 1.1 资料说明和处理

采用新一代多普勒天气雷达WSR-98D扫描强度场，该雷达体扫为每隔6 min一次。首先，将WSR-98D格式的雷达资料转化为WSR-88D格式；其次，为了方便计算，把极坐标格式的雷达资料插值到三维直角坐标系中。采用ARPS对雷达资料进行插值，插值后的资料共有53层，但由于雷达资料经常存在杂波和亮带等虚假回波，其中由于大气超折射异常的传播等虚假回波和地物杂波经常出现在低层，在此，我们选第7层(高度为2.2 km)的雷达资料进行外推，此层为能代表对流水平分布特征的一层雷达回波，因此，它比较具有代表性，并且在很大程度上能够避免由于大气超折射的异常传播等虚假回波和地物杂波的影响。由于噪声等具有低值回波特性，这些低值回波对算法有很大的影响，为了避免由噪声带来的外推误差，我们对外推回波进行最低值限制。为了能够更好地追踪降水系统(包括一般性对流降水系统，对流云和层状云混合降水系统)，本文把外推雷达回波的最低值定在10 dBz^[11]。

1.2 方法

线性交叉相关外推算法是目前临近预报的主要算法之一^[15]。它主要是通过计算雷达回波等资料在连续时次的空间最优相关，得到对流系统不同位置的移动矢量特征，并基于这些获得的移动矢量对雷达回波等进行外推，从而达到预报的目的。交叉相关外推算法考虑了回波移动矢量大小和方向的变化，而集合交叉相关外推算法考虑到多时刻回波的演变过程，考虑到回波的移速、移向以及回波在移动过程中的形变具有非线性，还考虑了整个回波在移动过程中的演变过程，所以，集合交叉相关外推算法在临近预报中具有一定的物理意义。

我们取通过插值后两个时刻(t₁, t₂)的第7层雷达资料，将这两个时刻的雷达回波分成若干个“小区域”，将t₁时刻每个“小区域”的雷达回波在搜索半径范围内分别与t₂ (t₂=t₁+Δt)时刻每个“小区域”的雷达回波进行空间交叉相关计算，相关系数R表示为

$R = \frac{{\sum\limits_k {\left[ {{Z_1}\left( k \right){Z_2}\left( k \right)} \right] - \frac{1}{N}\sum\limits_k {{Z_1}\left( k \right)\sum\limits_k {{Z_2}\left( k \right)} } } }}{{{{\left[ {\left( {\sum\limits_k {Z_1^2\left( k \right) - N\bar Z_1^2} } \right)\left( {\sum\limits_k {Z_2^2\left( k \right) - N\bar Z_2^2} } \right)} \right]}^{1/2}}}}$

(1)

其中Z₁、Z₂分别为t₁、t₂时刻这些“小区域”各个格点上的反射率因子, N是“小区域”内的格点总数。将t₁时刻的每个“小区域”与t₂时刻的搜索半径内所有“小区域”求相关以找出相关系数最大的“小区域”, 则t₂时刻该最大相关系数的“小区域”中心即是经过时间步长为Δt内该“小区域”内所有回波移动矢量的终点(见图 1交叉相关示意图)。

搜索半径太小，第一时刻“小区域”与第二时刻最大相关的“小区域”有可能不在搜索半径范围之内，这样找到的就不是最大相关，就会把最大相关的“小区域”漏掉。因此，所求的回波移动矢量是错误的。我们以低空急流的速度作为参考来决定搜索半径的大小(R_搜索半径=V_{低空急流平均速度}×T_{回波间隔时间})。

交叉相关外推的结果受到“小区域”大小的影响，对于不同的系统，我们选取“小区域”的大小也不一样。“小区域”太小，包含的数据信息太少，就得不到很稳定的相关关系，回波的移动方向比较多，比较凌乱；“小区域”选得合适，就会看到回波的合并与演变过程。“小区域”太大，就会把一些中小尺度的扰动平滑掉，上面所说的回波的合并与演变会被平滑掉，对于不同的天气系统和不同的分辨率，我们选取“小区域”的大小也会不同。为了避免噪声低值回波的影响，在检查到每个“小区域”的回波值小于10 dBz的格点总数超过该“小区域”格点总数的40%时^[11]，该“小区域”将不被作为我们追踪的对象。我们选的两个“小区域”的矩形中心的距离要小于“小区域”的大小，这样，相邻矩形“小区域”之间就有重叠的网格。

线性交叉相关外推的算法比较简单，计算量不大，很容易实现，由于它仅仅需要两个时次的雷达回波资料，根据这两个时次的雷达回波资料线性地外推第三时刻的雷达回波资料，这样仅考虑了水平方向上回波移动的大小和方向的变化。在此，我们在基于线性交叉相关算法的基础上，采用集合的方法，利用多时刻的雷达回波的资料，对这些雷达回波资料(两个时刻一组)求交叉相关，得出它们的移向和移速，再根据权重系数，最后得出未来1个小时内雷达回波的变化图像。此文采用集合的方法，如图 2所示，t₂保持不变，t₁分别取00：00、00：06、00：12、00：18、00：24和00：30这6个时次(世界时，下同)，然后将这6个t₁分别与t₂组合外推出t₃时刻的雷达回波，再将得到的6个t₃时刻的雷达回波集合，本文采取平均的方法，最后得到t₃时刻的雷达回波。这个方法不仅考虑了雷达回波的移向和移速，同时还考虑了雷达回波在过去1小时之内的演变过程。

雨强是将回波强度经过Z-I(Z=AI^b)关系转化得到的。Z为反射率因子，I为雨强，A和b为参数。

参数A和b的确定方法有最优化处理法^[15]，概率配对法等。我们这里采用最优化处理法。由雨量计测量每小时降水量G_i，同时由雷达测量这些雨量计上方对应的放射率因子Z_i，假设一个Z-I关系，把这些Z_i转化成雨强I_i，再把I_i值对时间进行累加，得到由雷达资料估算的每小时累积降水量A_i，最后通过判别函数CTF=min[∑(A_i-G_i)²+∣A_i-G_i∣]达到最小值来判断A和b，此时参数A和b所确定的Z-I关系为最优化的。

2 预报的评分标准

按照北京市气象局1小时降雨量级小雨、中雨、大雨、暴雨的分级标准分别为0.1、2.6、8.1和16 mm做了POD、FAR、ETS、CIS、BIAS评分。两类预报取值见表 1。

$POD = \frac{{{N_{11}}}}{{{N_{11}} + {N_{21}}}}$

(2)

$FAR = \frac{{{N_{12}}}}{{{N_{11}} + {N_{12}}}}$

(3)

$CSI = \frac{{{N_{11}}}}{{{N_{11}} + {N_{21}} + {N_{12}}}}$

(4)

$\begin{array}{l} ETS = \frac{{{N_{11}} - {N_{{\rm{11random}}}}}}{{{N_{11}} + {N_{21}} + {N_{12}} - {N_{{\rm{11random}}}}}}{N_{{\rm{11random}}}}\\ \quad \quad = \frac{{({N_{11}} + {N_{21}})({N_{11}} + {N_{12}})}}{{{N_{11}} + {N_{12}} + {N_{21}} + {N_{22}}}} \end{array}$

(5)

$BIAS = \frac{{{N_{11}} + {N_{12}}}}{{{N_{11}} + {N_{21}}}}$

(6)

POD反映命中的预报站数占总降水站数的概率，值越大越好；FAR用于衡量空报的概率，值越小越好；CSI为临界成功指数，值越大越好；ETS用于衡量模式去掉随机命中后的模拟命中概率，当ETS大于0时，说明模式优于随机命中，值越大越好，当ETS小于0时，说明模式劣于随机命中；BIAS用于描述模式预报降水的偏干、偏湿程度，越接近于1越好。

3 两个个例分析 3.1 山东一次暴雨过程的试验

2009年8月17日，山东地区出现一次自西向东移动的暴雨过程。我们选取17日05—06时(世界时，下同)的雷达资料为例，利用集合交叉相关追踪的方法进行试验。从回波强度CAPPI以集合交叉相关的矢量场回波的移向和移速，预报未来1小时内(逐6 min外推)雷达回波的位置与强度(见图 3，因篇幅所限，图 3只给出了每隔18 min的回波图)。

如图 3所示，利用集合交叉相关外推算法，获得这次降水过程未来1小时内的雷达回波的位置和发展形势变化与未来1小时雷达观测到的实况很相似。整个回波都呈西南—东北走向，回波的强中心与实际观测对应比较好，特别是在回波的左下方，整个外推演变过程回波与实况一致。雷达所在的中心位置左右两侧的单体回波在整个外推过程中与实况相似。从图中可以看出，回波反射率在30 dBz以上所覆盖范围要在山东以北的小单体预报中，预报的雷达反射率与实际观测的雷达反射率相比，强度要偏弱，因为05—06时，山东以北的小单体的回波在这段时间内是逐渐变弱的，而06—07时，这些小单体在增强。这正是集合交叉相关外推技术所存在的缺陷，它可以把回波的移动位置与形变预报得很好，虽然它利用集合的思想，考虑单体的演变过程，也考虑到回波强度的变化，但它不能够预报出雷达回波是否增强，而在山东以北的回波单体预报中，预报前1小时回波单体是减弱的，但是在此时，回波强度又开始加强，所以预报的回波强度要比实况偏弱。图 4给出了17日07时北京雷达实况和集合交叉相关外推及最小二乘拟合外推的对比，由图 4可见，最小二乘拟合外推法外推的雷达回波与实况相比，回波比较零碎，在济南东部和南部，有很多散乱的回波中心，而在实况上是一条强的带状回波中心，集合算法得到的回波比较集中，与实况比较一致。

预报的1小时累计降水量与观测的1小时累计降水量(图略)对比可以看出，1小时累计降水量的实际观测的图像与预报的图像还是很相似的，雨区和降水量级还是比较相同的，尤其是山东西南角的降水预报与实况对应非常好。预报的降水范围比实况偏大，因为预报的1小时累积降水量是由外推的雷达回波的反射率通过Z-I关系反演出降水强度，再对降水强度进行时间积分得到的，因此预报的降水覆盖是与雷达反射率所对应的范围一致，这也是利用Z-I关系预报降水的缺陷。

表 2为2009年8月17日06—07时山东地区雷达预报的1小时累计降水各个量级与实况的评分结果。POD总体来说是比较高的，就连大雨都达到0.4以上，其中小雨很高，达到了0.9599，这说明小雨几乎都报出来了，中雨和大雨大多数也都能够报出来，都超过了0.7。由于预报的降水范围比实况的降水偏大，所以空报率略偏高，但是也只是暴雨偏高一点，其他的还是比较正常的。小雨、中雨、大雨的临界成功指数比较高，都达到0.5以上。ETS均高于0.2，说明模式的预报概率远远优于随机模拟命中率，具有一定的预报水平。BIAS均大于1，表明模式预报降水实际偏湿，这也是由于预报降水的范围偏大所导致的，是由于Z-R关系反演降水存在的不足。

3.2 北京一次暴雨过程的试验

2007年8月1日，北京地区出现一次由西北向东南移动的暴雨过程。我们选取1日13—14时的雷达资料为例，利用交叉相关追踪的方法进行试验。从回波强度CAPPI以集合交叉相关的矢量场回波的移向和移速，预报未来1小时内(逐6 min外推)雷达回波的位置与强度(见图 5)。

如图 5可示，采用集合交叉相关外推算法，得到这次降水过程未来1小时内的雷达回波的反射率的移动方向与未来1小时雷达观测到的实况基本一致，都是由西北向东南移动。前18 min，预报的回波的移动位置与发展趋势与实况相比，两者基本上是重合的，说明预报与实际观测是一致的，北京南部的强回波带预报很好，与实况的位置是重合的，强度也与实况相当；北京东部的回波单体的位置与强度也是一致的。但是从24 min开始，预报的移动速度与实况相比是慢慢变快的，这是由于在实况图上可以看到，外推之前的1小时以内雷达回波(图略)的位置移动速度要远远大于外推时刻的未来1小时以后的回波移动速度。随着外推时间步长的增加，回波强度与实况相比慢慢偏弱，因为由集合交叉相关技术得到的回波的移动矢量的一致性随着外推时间步长的增加而变弱，外推时间步长越长，外推的误差也就变得越大，并且集合方法，外推的误差增大时，加上权重系数以后的结果也会使回波强度偏弱。在强回波带的西段和东南端，都有新的单体生成，这点外推技术就没有办法预报出来，外推技术只是参考已经出现的回波前一段时间的移动矢量，以这个移动矢量来计算回波未来的移动矢量，而后来出现的单体式没有办法计算到得，这也是外推技术存在的局限性。如图 6对比可见，最小二乘拟合外推法外推的雷达回波与实况相比，在北京西部，回波中心偏弱，范围偏小，而集合算法得到回波的范围和强度和实况比较一致。与最小二乘拟合方法相比，集合方法得到的回波更接近与实况。

统计表明(表 3)，对小雨和中雨的预报命中率还是很高的，其中小雨全部都能够报出来，总的来说，小雨和中雨的预报比较好，大雨和暴雨的预报效果不是很理想，这也是对回波的外推偏弱所导致的。

4 结论

本文基于线性交叉相关外推的算法，采用集合的思想，分别对山东和北京的两次暴雨过程的雷达回波与1小时累计降水量进行分析，得到以下结论。

(1) 集合交叉相关外推技术能够跟踪回波的移动矢量，能够反映降水系统所对应回波的移动方向和速度，集合方法考虑到多时刻对流单体的演变过程，采用非线性特征对回波的位置和形状的预报与实况还是比较一致的，所以该方法可以在临近预报业务中应用，本文外推1小时内回波演变过程，这符合张杰^[16]所指出的，理想的外推时间应选在1小时之内。

(2) 采用Z-I关系对未来1小时累积降水的预报，预报雨带走向与实况还是基本上一致的，整个雨带的形状与实况基本相同，对降水进行TS评分，评分结果比较高，总体效果比较可信。

(3) 采用集合交叉相关外推技术预报准确率比较高，与最小二乘拟合的方法外推相比，预报效果要好，回波移动的位置与实况较为一致。集合方法，采用多时刻雷达回波资料进行外推，这样就考虑了雷达回波在过去1小时之内的演变过程，考虑了雷达回波在过去1小时内演变是非线性的，与实际情况相符合。

需要指出的是集合交叉相关算法也存在局限性：它是参考前1小时雷达回波的移动矢量和发展情况，对未来1小时雷达回波进行计算，这就可能会出现计算回波的移动速度会偏快或偏慢，回波强度会偏强或偏弱。集合的组合有很多种，不同的集合组合对预报效果是否会有很大的影响？希望在今后的研究中能够深入探讨这个问题。从理论上讲，集合交叉相关外推技术比传统的交叉相关外推算法效果要好，在下一步工作中，将会试验大量的个例来验证这种方法的有效性。本文只是对外推1小时的回波进行比较，详细的比较将在今后的工作中进一步展开。