引言

对雷达回波进行外推，即对回波的运动进行跟踪。早期的回波运动跟踪都是基于雷达回波反射率因子图像，使用常用的TREC^[1-9]技术，它既适合于对流云也适合于层状云的跟踪，但不能给出变化趋势。

Hilst等(1960)^[1-4]先后将交叉相关技术应用于跟踪风暴簇的整体移动。Rinehar t等(1978)^[5]首先发展了用相关法跟踪回波运动的TREC(Tracking Radar Echo by Correla tion)技术，该技术可以较好地反演风暴内部运动。Rinehart(1981)^[6]又提出建议在应用TREC技术时，时间间隔最好小于5分钟，否则由于回波型的变化将导致无序矢量的产生。Smythe等(1983)^[7]和Tuttle等(1990)^[8]沿用同一概念反演出边界层内的气流。Li等(1995)^[9]在TREC技术的基础上发展了一种TREC矢量的连贯性技术，即COTREC(Continuity of TREC vectors)技术，使所得矢量场较为连续。之后，TREC技术在估计回波运动场中得到广泛应用^[10-12]。乔春贵等(2006)^[13]使用质心法对整幅雷达回波图进行线性外推，所得结果对稳定性的层状云回波有较好的外推能力，而对其它类型的云的回波外推能力较弱。张亚萍等(2006)^[14]在常用的TREC技术基础上，发展了DITREC(Difference Image-based Tracking Radar Echo byCorrelations)技术，DITREC矢量场的时间连续性和空间连续性好于TREC矢量场，但在外推中没有考虑回波生消的变化，需要三张间隔5分钟左右的反射率因子图像序列才能生成一个DITREC矢量场, 并且DITREC矢量数少于TREC矢量数。

Wang Z H等(1977)^[15]提出相关技术追踪云块的位移时存在“亚像素位移”问题，用傅立叶相位分析技术能避免该问题，且能够计算出处于发展或消散阶段的云块的运动矢量。王振会等(2005)^[16]提出TCFM(Tracking cloud with combined Fourier phase analys is and Maximum correlation)技术，该技术能在最大相关法整数倍像素位移的基础上进行亚像素尺度修正，进一步提高目标云块的追踪精度。

当地C波段雷达每次体扫用时约7分钟，TREC外推中必然有无序矢量的产生，本文在COTREC和质心跟踪技术的基础上，再使用傅立叶相位分析技术，它能计算出1分钟间隔云块的移动速度，处理处于生消阶段的回波的一定程度的变化，使无序矢量的产生得到更多地抑制，并且通过计算“亚像元”速度以提高追踪精度。将外推回波强度及1小时降水量图像应用到预报实践中，为当地短时临近降水预报提供参考。

1 前提条件、方法说明

回波图像块在5分钟之内不能发生很大变形或翻转，其变化必须满足平稳随机过程，这是任何使用相关技术、质心跟踪以及傅立叶相位分析技术处理雷达回波的前提。

为提高计算速度和减小一维傅立叶相位分析的误差，采用二维快速傅立叶相位分析法^[17](简称FFT2D)。外推回波的时间步长为7分钟。首先根据相邻两个时次(t₀, t₁)的雷达强度回波数据，在M×M像素的回波图像内，以N×N像素为外推计算单位，将t₁时次作为目标图像块的所在回波图像，用COTREC法与FFT2D得到回波块从t₀运动到t₁的移动速度，再用质心法得到对应回波块的质心位置，并由强度数据得到质心强度，从而得到质心强度变化率，作为对应外推回波块的强度变化率，在t₁时次回波图像数据的基础上，根据各回波块的移动速度和强度变化率得到外推的t₂时次的各回波块的位置和强度数据。计算过程见图 1的说明。

图 1中，M×M表示强度回波图像的尺寸(500×500)，N×N表示回波的外推计算单位(10×10像素的回波块)，P₀、P₁分别表示t₀、t₁时次的对应匹配回波图像块的质心，P₂所在回波块是由P₀和P₁所在回波块外推所得的t₂时次回波块，Z_P0、Z_P1、Z_P2分别表示质心处的强度，u、v分别表示回波块从t₀时刻运动到t₁时刻的x、y、方向的速度，令(i_t2，j_t2)表示外推的t₂时刻回波块的像素坐标，(i_t1，j_t1)表示t₁时刻回波块的像素坐标，则

$ \left\{ \begin{array}{l} {i_{t2}} = {i_{t1}} + u \times \left( {{t_2} - {t_1}} \right)\\ {j_{t2}} = {j_{t1}} + \upsilon \times \left( {{t_2} - {t_1}} \right) \end{array} \right. $

令r_at表示回波块从t₀到t₁时刻的质心强度变化率，则

$ {r_{at}} = \frac{{{Z_{p1}} - {Z_{p{\rm{o}}}}}}{{{t_1} - {t_{\rm{o}}}}} $

作为P₂所在回波块的强度变化率，则在t₁时刻回波块的基础上，外推的t₂ 时次回波图像上任意像元的回波强度为:

$ \begin{array}{*{20}{l}} {\sum\limits_{}^n {} \sum\limits_{}^N {} \sum\limits_{}^N {{Z_{k(i,j)t2}} = \sum\limits_{k = 1}^n {} } \sum\limits_{j = 1}^N {} \sum\limits_{i = 1}^N {{Z_{k(i,j)t1}} + } }\\ {\sum\limits_{k = 1}^n {} \sum\limits_{j = 1}^N {} \sum\limits_{i = 1}^N {{r_{atk(i,j)}} \times ({t_2} - {t_1})} } \end{array} $

其中，k=1, 2, 3, …n表示10×10回波图像块的顺序号，n表示回波图像内共有50×50个回波块，(i, j)表示第k个回波图像块内的像元坐标。

最后计算外推t₂时次与实况t₂时次回波的强度相关系数(K)，若K大于域值0.6，则认为此次外推有效，可以显示出外推的t₂时次回波图像。同理，在t₁的基础上外推t₃、t₄…时次的回波。将有效外推回波结合当地早期统计所得的本地化Z-R关系参数，得到雨强分布图像，并与WSR_88D的PUP产品(1小时降水量：OHP)，以及自动雨量站的记录实况做比较分析。

2 数据来源

本文所采用的数据是西宁C波段天气雷达体扫基数据，通过读取基数据并对每4个连续的库取平均回波功率得到反射率因子(强度)数据，对9个仰角上的极坐标数据进行双线性插值，得到1~18km高度上的CAPPI数据(并且本文垂直剖面数据由此得到)。另外，以1km×1km为底面积、从地面到云顶的柱体中找到回波强度最大值的资料点，从而得到组合反射率数据，以0. 1°为插值单位对每个高度层上的360个方位进行资料填补等，将极坐标系下的回波数据转换为直角坐标系下的回波数据，回波图像的像素分辨率是1km/pixel，在500×500像素的回波图像区域内用上述方法进行外推。本文使用的CAPPI高度缺省是3km，显示范围缺省是25 0km。以2007年8月25日以对流为主的混合性降水、2007年10月31日的层云降水过程为例。

3 结果比较与分析 3.1 回波图像的对比分析

图 2(见彩页)是2007年8月25日的外推回波和实况。从15:19开始，对流在西宁的西面和北面产生并发展(如图 2a所示)，到15:40时还处于发展阶段(图 2b)；20:08开始成熟(图 2c)，可见在西宁周围的回波形成了强中心带，20:33成熟的对流从乐都伸展到民和(图 2d)，大通站附近的回波局地顶高达到9km，强中心高度为1km，位置偏下(图 2dt)，20:52回波形状、强度分布无明显变化(图 2e)，回波顶高维持不变，强中心高度上升(图略)，20:58强中心高度上升到2km(图 2g)，此时回波发展到最成熟，之后回波强中心带的强度和顶高及强中心高度开始下降；从23:50开始回波正式处于衰亡阶段(图 2f)。

外推的回波图像与实况基本一致。15:40发展阶段的外推回波图像如图 2bx所示，是由15 :19(图 2a)和相邻前一时次的体扫回波数据外推而得，外推时效是21分钟，即3个外推时间步长。可见回波形状与图 2a相似，回波位置与其实况图 2b基本一致。20:33成熟阶段的外推回波图像如图 2dx所示，是由20:08(图 2c)和相邻前一时次的体扫回波数据外推而得，可见除了民和站所在的线状回波(图 2dx)不具备外，其余大部分地区的回波形状和强度分布均与实况图 2d几乎一致；20:52图 2ex所示回波是由20:33(图 2d)和相邻前一时次的体扫回波数据外推而得，回波形状和强度分布与其实况图 2e基本一致。23:50衰亡阶段的外推回波图像如图 2fx所示，可见与其实况图 2f也是接近一致的。对于对流降水回波外推视其回波演变情况而定，通常情况下，最大外推时效为5个时间步长，发展阶段的回波外推时效为2个外推时间步长，成熟阶段的外推时效为4个时间步长，衰亡阶段的外推时效为4个时间步长。

图 3(见彩页)是2007年10月31日的外推回波和实况。20:51开始层云在大通、门源、互助三县产生并发展(图 3a)，21:24此三县层云已向西宁方向移动(图 3b)，22:30此三县层云回波的强度增强，尺度扩大，且湟源、湟中境内有新生回波，此时层云开始成熟，23:04达到最佳成熟状态，之后回波逐渐减弱直到消亡。发展阶段的外推回波如21:24(图 3bx)所示，由于回波形状、强度分布及移动速度的变化较大，外推所得图 3bx与实况图 3b有较小差异存在，例如图 3bx在门源境内还有回波而实况已经加速移动到了大通县内。23:04的外推回波(图 3dx)与实况(图 3d)位置基本一致，但强度分布有差别。这些差别的存在是由于TREC不能给出回波变化趋势的特性决定的。此次过程中的外推最大时效是4个外推时间步长，最小为2个步长。对于层状云降水回波的外推仍然视其回波演变程度而定，通常情况下，对较稳的层云回波的外推时效最大能达到8个外推时间步长，即约60分钟。

3.2 平均绝对差异分析

图 4(见彩页)和图 5(见彩页)分别是2007年8月25日和10月31日由COTREC外推所得回波的图像，在此只需分别列出一张图像来说明COTREC外推的结果，可见图 4和图 2(dx)、图 5和(图 3(dx)都是有区别的。为了得到具体差别程度，以及说明计算“亚像元”速度的必要性，下面计算外推回波和实况的各质心位置和强度，将计算“亚像元”速度(S)和不计算“亚像元”速度(S)的外推结果，与实况(S₀)进行比较，分别得到质心位置及强度的平均绝对差。

如表 1所示，质心的比较是用外推回波与实况相同位置的回波块的质心位置及强度的比较。表 1中的序号Ⅰ、Ⅱ、Ⅲ分别代表发展、成熟、衰亡阶段，ⅰ、ⅱ、ⅲ同理。ΔX表示外推回波与实况的质心位置在水平方向的平均绝对差，ΔY表示垂直方向位置的平均绝对差，ΔZ表示外推回波与实况回波对应质心处的强度平均绝对差，单位是dBz。从表 1中可以看出，Ⅰ、Ⅱ和ⅰ、ⅲ中用S法外推所得ΔX、ΔY、ΔZ的值均大于用S法的值，Ⅲ中S法所得ΔY、ΔZ的值分别大于S法的值，ⅱ中S的ΔZ大于S 的值6dBz。可见其一般规律是，外推时计算“亚像元”速度要比不计算“亚像元”速度的结果更接近实况，回波演变越快，外推回波与实况的差异越大，并且不计算“亚像元”速度的ΔZ比计算“亚像元”速度的值越大。

3.3 改善效果分析

令S_x、S_x、S_0x分别表示S、S、S₀所得质心在水平方向的位置，IMD表示计算“亚像元”速度的外推比不计算的改善程度

$ IMD = \frac{{({{\bar S}_x} - {S_{0x}}) - ({S_x} - {S_{0x}})}}{{{S_x} - {S_{0x}}}} $

若IMD>0，则表示S更准确；若IMD＜0，则S更准确；若IMD=0，则两种方法同效。在垂直方向的质心位置及质心处的强度的改善程度可同理计算。

表 1中Ⅱ时次的X方向的S_x-S_0x=3、S_x-S_0x=2，所以，IMD>0，用S更准确。Ⅰ到Ⅲ及ⅰ到ⅲ的IMD均不小于0，可见用S比用S的准确概率更大。由表 1可得，从发展到衰亡的回波变化过程的平均IMD，2007年8月25日外推位置改善x方向约20%、y方向约30%、相同位置处的强度改善约43%，2007年10月31日分别为17%、15%、50%。因此，改善程度同样视回波演变速度而定，一般情况下，外推位置改善至少15%、强度改善至少40%。这对于气象业务中要求更准确地预报降水落区甚至落点、降水强度、影响范围等将起到一定作用。另外，在计算时间方面，此方法比COTREC用时约多1500ms，但这并不影响降水的预报。

4 外推雨强分布及对比分析

下面将WSR_88D的PUP所做的1小时降水量产品(OHP)、由本文外推数据所得雨强分布(OHP)、自动雨量站的记录相互比较。例如，2007年8月25日20:52的地面雨量站记录是21:00—22:00，西宁城西6.4mm、湟中5.0mm、其余地区无降水；22:00—23:00西宁城西7.6mm、湟中10.3mm、其余地区无降水。但实际在大通站到西宁北部之间有强降水出现，而自动雨量站无记录。PUP的OHP产品的估计结果如图 6(见彩页)(a0和b0)所示，20:52—21:52无降水；21:17—22:17湟中和大通有最大降水量91.4mm、其余地区空报。由20:33的CAPPI外推所得OHP如图 6(a1)和图 (6b1-b3)所示，20:52西宁城西和湟中北部的降水在6.35mm附近，与实况接近，但其余地区空报(图a1)；21:17在1km高度层上西宁城西、湟中的降水平均在6.35~12.70mm之间，大通到西宁之间的降水强度则更大(图 6b1)，与实况更接近而PUP则差距太大，特别是湟中北部PUP估计值比实况大约平均30mm(图6b0)。在2km高度层上主要降水带仍然在西宁城西、湟中和大通，但雨强略有减弱(图6b2)，在3km高度层上该三地依然是主要降水带，且雨强比2km高度稍微减弱(图6b3)；从20:33到20:58回波强中心带的位置仍然分布在西宁城西、湟中和大通、最大回波强度稍微减弱但变化不明显图 7(见彩页)和图 8(见彩页)所示，并且，由前几个时次的CAPPI图像(图 2)、以图 2dt)和(图 2gt)为代表的前几时次的垂直剖面等，可以得到最大回波强度约50dBz，回波顶高约8km，强中心高度偏下，回波已经触地，强度和高度变化不明显，回波无明显外型特征等PUP指标，据此指标，本文所用系统的自动预报结果是：未来2小时内西宁城西、湟中和大通的降水仍然维持，且降水量可达到中雨量级。此结果与雨量自动站出现的实况降水位置一致且降水量级别相同。可见，本文外推的雨强与实况更加接近，但与PUP一样易空报，而PUP更易漏报。外推OHP与PUP的OHP同样能做到约每7分钟更新一次，而且外推OHP可直接将回波做最大时效的外推，使预报时间更长。

5 小结与讨论

(1) 外推技术对于对流云、混合云、层状云降水均适用，并且，将计算“亚像元”速度应用到外推算法中，能够将某些处于生消阶段的对流纳入计算，对于以往的TREC算法能起到一定改进作用，从而提高外推运算的准确率。

(2) 外推的时效，对处于发展或消亡阶段的对流云视其变化速度而定，一般在30分钟以内；对相对稳定的降水回波通常约35分钟有时甚至1小时。

(3) 在短时临近降水预报中，用本文算法做最大时效的外推，经初步检验，不会存在漏报现象，可以做到较准确地预报降水发生的落区甚至落点，但空报现象与PUP产品同样存在，并且外推降水量的准确度受Z-R关系的影响。

(4) 当地天气雷达每做一次体扫所用时间约7分钟，回波在此时间段内极易发生形状和强度分布的剧烈变化，而不满足本文算法的前提条件；对于生命史不到7分钟的对流，不能做外推运算。雷达体扫所用时间和云块的生命史是限制本文外推算法的客观因素。

经当地多次短时临近预报检验，结果证明本文的外推技术是适用的，能够在一定程度上提高预报准确率和相应地提高当地气象部门短时临近预报评分。对于该算法存在的主客观问题的解决，将在今后的工作中进一步深入研究。