2. 美国国家海洋和大气管理局 合作遥感科学与技术中心, 纽约 NY 10031, 美国
2. Cooperative Remote Sensing Science and Technology Center, National Oceanic and Atmospheric Administration, New York NY 10031, USA
长江三峡区间是地形复杂的山区流域,流域内支流多、坡降陡、汇流迅速; 同时位于大巴山南麓及鄂西北暴雨区内,区间产生的暴雨洪水将直接迅速注入三峡水库,对水库的防洪安全及运行调度影响很大。研究表明,现有地面雨量观测站网仍无法很好地描述降雨的时空分布特征[1, 2],基于雨量站观测的短期洪水预报还存在较大不确定性[3],未来天气雷达提供的高时空分辨率降雨定量估测信息则有可能用于提高洪水预报精度[4, 5]。另一方面,由于区间暴雨洪水汇入迅速,适当增加预见期可大幅提升洪水预报的价值,因此需要进一步发展有效的降雨定量预报方法。
短时降雨定量预报(简称临近预报)是水文气象学中的重要课题之一[6],它是指对未来短时间内 (0~3 h)降雨情况的定量预报。概括而言,临近预报方法主要分为2类: 数值天气预测模式和基于遥感观测的外推方法[7]。数值天气预测模式,如Weather Research and Forecasting Model基于微云物理方案进行预报,能够预报降雨的形成、发展、消亡,具有很强的物理基础,但其结果通常对初始条件、空间精度、同化算法等因素非常敏感[8]。近年来,大量研究致力于改进数值天气预测模式的预报精度,但计算速度慢且短时预报精度低等问题仍难以有效克服[9, 10]。基于遥感观测的临近预报方法则是根据之前时刻的观测结果(如雷达、卫星观测)找出降雨运动、变化的规律,并以此规律推断出下一时刻降雨分布的方法。基于观测外推的临近预报算法简单有效,短时预报效果甚至好于基于微云物理方案的数值天气预测模式,因此在水文气象研究中有着广泛地应用。
传统上,基于观测的临近预报外推算法如 TITAN算法、Eulerian-Persistence算法、Lagrangian -Persistence算法[11, 12]等,多是通过直接的空间相关来计算降雨的位移,这类算法能够预报降雨的运动,但无法很好地刻画降雨的旋转与变形。山区复杂地形往往会导致降雨时空分布不均匀,局部地形效应使降雨向高海拔区域集中。针对山区临近预报的这一困难,Zahraei等提出了一种新的基于网格追踪的临近预报外推(pixel-based nowcasting,PBN)算法[8]。该算法可以利用高时空精度的雷达数据,能够有效追踪小尺度降雨的平移与变形,尤其适用于有雷达观测区域的短时降雨定量预报。
本研究将PBN算法应用于长江三峡区间的降雨短时预报,利用万县S波段雷达的观测到的2010年汛期的主要降雨过程进行预报检验,定量评估该算法在三峡区间的应用效果,并进一步讨论区分该算法对不同降雨类型的预报精度。
1 PBN算法PBN算法是一种基于网格追踪的临近预报算法,主要适用于雷达、卫星等遥感栅格数据。其基本操作分为2个阶段:
1) 识别与追踪: 根据当前时刻(t)以及之前2个观测时刻(t-Δt,t-2Δt;Δt为相邻时间间隔)的雷达图像识别出降雨变动的动态信息;
2) 外推: 基于此信息再通过定量演绎外推得到下一观测时刻(t+Δt)的降雨分布。具体来说,第一阶段的识别与追踪是通过寻找本时刻每一个网格上的降雨在上一时刻的位置,确立对应关系; 根据这个对应关系,可以获得每一个网格对应的位移向量以及降雨强度变化。 第二阶段是基于对应关系进行外推,推求下一时刻降雨的分布。
1.1 识别与追踪目前多数临近预报的算法均使用模块匹配(template-matching algorithms)的方法识别及追踪降雨。这种方法是寻找目前图像中指定大小的窗口在上一幅图像中可能的位置,通过定义并计算相似性评价指标(如相关系数)来确定最佳匹配位置。这种算法对空间精度较高的观测图像应用效果较差,经常会出现弥散或多个最优解的现象[13, 14],而且这种仅依靠窗口最优相关的方法无法模拟观测图像(降雨遥感影像)的局部旋转与变形[15]。上述问题可采用位移的分层表示方法来解决,分层法将时序图像的运动特征分解为一系列较粗网格的变化趋势叠加上较细网格的小幅局地修正[16]。目前已有若干基于网格数据和分层技术的追踪算法[17],PBN采用了此类追踪算法并将降雨的变化视为一种矩形的拓扑变换,主要包括以下5个步骤[18]:
1) 基于相关系数,寻找基础图像(如t-Δt时刻图像)各网格的中心点在参考图像(如t时刻图像)上最相关的对应位置;
2) 将基础图像与参考图像的网格节点分别替换为步骤一中基础图像各网格的网格中心点和其在参考图像中的对应点;
3) 检查新的基础图像中相邻格点的一致性,不一致的格点由交互相关匹配替代;
4) 校正边缘效应,将基础图像和参考图像同时线性插值为原来空间精度的2倍;
5) 去除插值中产生的凹四边形网格。
1.2 外推外推过程如公式(1)—公式(5)所示。
Pt+Δt(xt+Δt,yt+Δt)= min{Pt(xt,yt)+ΔP,T}, |
(1) |
(xt+Δt,yt+Δt)=(xt,yt)+Δ(xt,yt). | (2) |
其中: Pt+Δt(xt+Δt,yt+Δt)是预报的下一时刻t+Δt时刻(xt+Δt,yt+Δt)网格处的降雨强度,它与当前t时刻(xt,yt)网格处的降雨强度Pt(xt,yt)的关系如式(1)和式(2)所示。
追踪算法找到每一个降雨网格在上一时刻的对应位置,例如Pt(xt,yt)的位置是(xt,yt),在上一时刻t-Δt时刻对应的是Pt-Δt(xt-Δt,yt-Δt),其位置是(xt-Δt,yt-Δt)。Vp表示的是Pt(xt,yt)的运动向量,它和当前时刻(t)、 上一时刻(t-Δt)以及上上时刻(t-2Δt)降雨的位置有关。如果Vp-2表示 t-2Δt时刻的Pt-2Δt(xt-2Δt,yt-2Δt)与t-Δt时刻的Pt-Δt(xt-Δt,yt-Δt)之间的位移,Vp-1表示t-Δt时刻的Pt-Δt(xt-Δt,yt-Δt)与t时刻的Pt(xt,yt)之间的位移,那么,如式(3)和式(4)所示,在PBN算法中,t时刻Pt(xt,yt)的位移向量 V p是这2个向量的线性叠加。
Vp=Δ(xt,yt)= f{(xt,yt),(xt-Δt,yt-Δt),(xt-2Δt,yt-2Δt)} = $\frac{1}{2}$(Vp-2+Vp-1). |
(3) |
V p=mean(Dist1,Dist2). |
其中:
Dist1= ((xt-xt-Δt)2+ (yt-yt-Δt)2)0.5, Dist2= ((xt-Δt-xt-2Δt)2+(yt-Δt-yt-2Δt)2)0.5. | (4) |
其中: | V p|表示的是Pt(xt,yt)位置变化的大小。
ΔP=g{Pt(xt,yt),Pt-Δt(xt-Δt,yt-Δt), Pt-2Δt(xt-2Δt, yt-2Δt)}= $\frac{1}{2}\left[\begin{array}{l} ({P_t}({x_t},{y_t}) - {P_{t - \Delta t}}({x_{t - \Delta t}},{y_{t - \Delta t}})) + \\ ({P_{t - \Delta t}}({x_{t - \Delta t}},{y_{t - \Delta t}}) - {P_{t - 2\Delta t}}({x_{t - 2\Delta t}},{y_{t - 2\Delta t}})) \end{array} \right],$ |
(5) |
为了避免观测噪声的污染以及相邻格点间的不连续问题,算法进一步在3×3个网格范围内进行了低通滤波处理。同理,降雨强度的变化ΔP也是由当前时刻(t)、 及前一时刻(t-Δt)、 前2个时刻 (t-2Δt)的降雨强度所决定。如式(5)所示,它是相邻两时刻变化值的平均值,当预报的降雨强度超过阈值T时,取该阈值作为该网格的雨强,当降雨强度小于0时,该网格降雨强度为0。
2 数据来源三峡库区内的万县雷达是国家气象局布设的S波段Doppler天气雷达。该雷达采用VCP21体扫模式: 即在6 min内完成9个仰角上半径460 km的扫描, 扫描范围覆盖整个三峡区间。万县雷达沿角度方向的空间分辨率为1°,沿径向反射率因子(Z)的空间分辨率为1 000 m,其主要参数见表 1,本文采用的是0.5°仰角的观测结果。
采用万县雷达2010年汛期观测到的11场典型降雨过程对PBN算法进行应用与验证,同时考虑降雨特征的差别,重点分析观测样本中4场代表性降雨过程(图 1)。6月7日降雨是面积稍小、形状复杂且移动的层状降雨; 6月24日降雨是包含多个快速变动雨团的对流降雨; 7月3日降雨是含有多个对流雨团移动变形并伴随分离融合的降雨; 9月9日降雨是相对稳定的大面积层状降雨。
![]() |
图 1 4场代表性降雨过程初始时刻雷达反射率空间分布(2010年) |
研究表明,临近预报所选取的时间步长以 10 min为宜[19]。结合三峡雷达数据观测特点,本研究分别采用6 min与12 min进行测试,结果显示,采用 6 min时间步长预报并未显著提高半小时预见期内的预报精度,整体预报精度也显示12 min设置下预报效果更好。因此,本研究设定每12 min预报一次,直至180 min。
4 评价指标对于一个二元事件,观测与预报结果都会有发生、不发生这2种可能。观测发生且预报发生的事件称为正确预报,观测未发生而预报发生的事件称为错误预报,观测发生而预报未发生的事件称为漏报,观测未发生且预报未发生的事件称为合理拒绝。对于降雨事件而言,通常会设定一个阈值(TR),降雨强度大于阈值认为发生,小于则认为不发生。评估指标计算方法如表 2所示,其中,nh、 nfa、 nf、 ncn分别表示正确预报(hit)、错误预报(false alarm)、漏报(failure)、合理拒绝(correct rejection)的数量。
预报值F | 观测值O | 总计 | |
>TR | <TR | ||
>TR | nh | nfa | nh+nfa |
<TR | nf | ncn | nf+ncn |
总计 | nh+nf | nfa+ncn | nh+nf+nfa+ncn |
本研究采用4种评价指标来验证PBN算法的准确性。除相关系数外,还采用水文气象学中常用的探测率、误报率、临界成功指数这3个指标,它们的定义如下:
1) 探测率(probability of detection,POD): 描述观测到阈值以上的降雨时间可以被正确预报的比例,即:
$POD = \frac{{{n_h}}}{{{n_i} + {n_h}}}$ | (6) |
2) 误报率(false to alarm,FAR): 描述预报的阈值以上的降雨事件中被误报的比例,即:
$FAR = \frac{{{n_{f\alpha }}}}{{{n_{f\alpha }} + {n_h}}}$ | (7) |
3) 临界成功指数(critical success index,CSI): 区别于探测率,该指标描述的是正确预报占所有可能发生降雨的情况的比例,即:
$CSI = \frac{{{n_{f\alpha }}}}{{{n_f} + {n_{f\alpha }} + {n_h}}}$ | (8) |
由于所有可能发生降雨的情况包括正确预报、错误预报和漏报,所以这个值通常会小于探测率。
5 验证结果表 3展示了全部11场降雨过程预报结果的相关系数。对于全部11场降雨,临近预报12 min的相关系数平均值可达0.80,36 min可达0.67,60 min 可达0.59。可以看到,在36 min预见期时,11场降雨中有8场降雨预报结果的相关系数大于 0.60,相关性较强; 在1小时预见期内,仍有7场降雨预报结果的相关系数大于0.60。相关系数低于0.60的结果出现在6月20日、 6月24日、 7月3日、 8月3日这4场降雨中,这4场降雨均具有对流型降雨的典型特征,包含多个对流雨团的生成与消亡。研究对象的尺度是影响临近预报准确性的重要因素之一[13],一般小尺度降雨的生命周期较短,而PBN追踪算法基于拓扑关系,短期内能够较好 地模拟降雨移动和变形,但没有考虑降雨的生成与消退,因此对局地强对流降雨的预报效果不佳。故当观测区域中存在显著的对流雨生成与衰退过程时,该算法精度将大幅下降。
日期 | 时间/min | |||
12 | 36 | 60 | 180 | |
20100607 | 0.82 | 0.69 | 0.60 | 0.23 |
20100618 | 0.84 | 0.76 | 0.70 | 0.50 |
20100620 | 0.72 | 0.50 | 0.47 | 0.41 |
20100624 | 0.77 | 0.65 | 0.59 | — |
20100703 | 0.71 | 0.43 | 0.37 | 0.12 |
20100710 | 0.86 | 0.78 | 0.72 | 0.49 |
20100715 | 0.87 | 0.72 | 0.61 | 0.21 |
20100803 | 0.59 | 0.34 | 0.22 | 0.08 |
20100814 | 0.83 | 0.74 | 0.63 | 0.42 |
20100905 | 0.90 | 0.85 | 0.82 | 0.73 |
20100909 | 0.90 | 0.86 | 0.80 | 0.41 |
平均 | 0.80 | 0.67 | 0.59 | 0.36 |
图 2表示了这4场代表性降雨过程预报值与观测值的相关系数随预见期变化的曲线。结果显示,对于着重关注的4场降雨,9月9日的相关系数最高,1小时预见期时相关系数可达0.80。6月7日的相关系数其次、 6月24日的相关系数再次,这2场降雨1小时预见期时的相关系数均接近0.6。7月3日的相关系数最低,1小时预见期的相关系数只有0.37。
![]() |
图 2 预报值与观测值的相关系数随预见期变化图 |
考虑到降雨事件二元分类评价时需要设定阈值,所以本文先采用了一种降雨评估算法将反射率转换为降雨,之后再对结果进行评价[4]。
图 3表示了在不同阈值下,4场代表性降雨过程的探测率随预见期变化的曲线。阈值分别设为 1 mm/h、 5 mm/h、 10 mm/h和20 mm/h。当阈值不断增大时,同一预报时刻探测率减小,这是因为阈值增大,观测值与预报值中超过阈值的点的数量减少,即正确预报的数量减少,虽然错误预报的数量也减小,但总体来说,探测率减小。对于同一阈 值下不同类型的降雨,可以看出,在低阈值(1 mm/h 和5 mm/h)时,9月9日的探测率最高,6月7日 的探测率次之,6月24日的探测率再次之,7月3日的探测率最低。这说明算法对独立且相对稳定、移动缓慢的大面积层状降雨(9月9日)探测率最高,2个阈值下1小时预见期探测率分别可以达到0.76和0.62。对多个对流降雨在短时间内移动、变形、融合成大的降雨 中心(7月3日)预报效果探测率 最低,2个阈值下1小时预见期探测率分别只有 0.22 和0.03。在高阈值时(10 mm/h和20 mm/h),由于9月9日降雨本身为层状降雨,降雨强度较少超过阈值,所以正确预报的数量较少,导致探测率也下降。这种情况下反而是6月24日和7月3日这2场对流型降雨初始时刻的探测率高。
![]() |
图 3 降雨过程的探测率随预见期变化图 |
图 4表示了在不同阈值下,4场代表性降雨过程的误报率随预见期变化的曲线。同样针对4个阈值1 mm/h、 5 mm/h、 10 mm/h、 20 mm/h分别计算。当阈值增大时,与探测率的结果相反,同一预报时刻误报率增大。这是因为阈值增大时,预报的大于阈值的降雨明显减少,这使得错误预报在预报的大于阈值的降雨情况中所占的比例增大。对于不同类型的降雨,在低阈值(1 mm/h和5 mm/h)时,9月9日的误报率最低,6月7日的误报率次之,6月24日的误报率再次之,7月3日的误报率最高。算法对独立且相对稳定、移动缓慢的大面积层状降雨(9月9日)误报率最低,2个阈值下1小时预见期错误预报率分别只有0.19和0.53。对多个对流降雨在短时间内移动、变形、并融合成大的降雨中心(7月3日)预报效果误报率最高,2个阈值下1小时预见期错误预报率分别达到0.79和0.89。在高阈值时(10 mm/h和20 mm/h),同样由于9月9日降雨本身降雨强度较少超过阈值,所以误报率与6月24日和7月3日的误报率较为接近。
![]() |
图 4 降雨过程的误报率随预见期变化图 |
图 5表示了4场代表性降雨过程的临界成功指数,临界成功指数与探测率变化的规律相似,9月9日的临界成功指数最高,7月3日的临界成功指数最低。
![]() |
图 5 降雨过程的临界成功指数随预见期变化图 |
本文基于长江三峡区间内万县S波段雷达2010年汛期的观测资料,应用外推临近预报PBN算法,定量评估了对三峡地区不同类型降雨的临近预报效果。
算法对三峡区间1 h内的降雨预报效果良好,平均相关系数达0.59。对于不同类型的降雨,在低阈值时,探测率、误报率、临界成功指数这3个评价指标的评价结果均显示: 算法对独立且相对稳定的大面积层状降雨预报效果最好,对包含多个对流型雨团生成与消亡的降雨预报效果最差; 在高阈值时,由于层状降雨本身降雨强度达不到设定阈值,所以初始时刻结果对流型降雨的预报效果较好。
如上所述,因原理所限,算法对包含对流型雨团生成与消亡的降雨预报效果不佳,所以在未来的研究中,需进一步发展具有物理基础且能同化其它气象观测资料的外推临近预报算法,使之能够在较短时间内更加准确的预报局地对流降雨。同时将临近预报结果与区域分布式水文算法结合,检验其延长水文预报预见期的效果。
[1] | 卢江涛. 利用天气雷达估测和预报降雨分布的研究 [D]. 清华大学, 2011.LU Jiangtao. Study on Precipitation Estimation and Nowcasting Based on Weather Radar [D]. Beijing: Tsinghua University, 2011. (in Chinese) |
[2] | 袁晓清, 倪广恒, 潘安君, 等. 基于最优化算法的北京市新一代天气雷达Z-R 关系研究 [J]. 水文, 2010, 30 (1): 1-6. YUAN Xiaoqing, NI Guangheng, PAN Anjun, et al. NEXRAD Z-R power relationship in Beijing based on optimization algorithm [J]. Journal of China Hydrology, 2010, 30(1): 1-6. (in Chinese) |
[3] | 潘安君, 卢江涛, 田富强, 等. 雷达测雨的最优化方法在分布式城市洪水算法中的应用 [J]. 中国水利, 2011(23): 54-55. PAN Anjun, LU Jiangtao, TIAN Fuqiang, et al. Radar rainfall optimization method in a distributed model of urban flood [J]. China Water Recourses, 2011 (23): 54-55.(in Chinese) |
[4] | 李哲, 杨大文, 洪阳, 等. 基于天气雷达的长江三峡区间降雨定量估测方法 [J]. 水力发电学报, 2014, 33(3): 5.LI Zhe, YANG Dawen, HONG Yang, et al. Radar-based quantitative precipitation estimate in the Three Gorges Region of Yangtze River [J].Journal of Hydroelectric Engineering, 2014, 33(3): 5 (in Chinese) |
[5] | 李哲, 杨大文, 田富强, 等. 基于地面雨情信息的长江三峡区间洪水预报研究 [J]. 水力发电学报, 2013, 32(1): 42-49. LI Zhe, YANG Dawen, TIAN Fuqiang, et al. Flood forecast for Three Gorges region of the Yangtze based on ground-observed rainfall [J]. Journal of Hydroelectric Engineering, 2013, 32(1): 42-49.(in Chinese). |
[6] | Berenguer F M, Sempere T D. Radar-based rainfall nowcasting at European scale: Long-term evaluation and performance assessment [C]//Proceedings of the 36th AMS Conference on Radar Meteorology. Breckenridge, USA: American Meteorological Society, 2013. |
[7] | Wilson J W, Ebert E E, Saxen T R, et al. Sydney 2000 forecast demonstration project: Convective storm nowcasting [J]. Weather & Forecasting, 2004, 19: 131-150. |
[8] | Zahraei A, Hsu K, Sorooshian S, et al. Quantitative precipitation nowcasting: A Lagrangianpixel-based approach [J]. Atmospheric Research, 2012, 118: 418-434. |
[9] | Golding B W. Nimrod: A system for generating automated very short range forecasts [J].Meteorological Applications, 1998, 5(1), 1-16. |
[10] | Seed A W. A dynamic and spatial scaling approach to advection forecasting [J].Journal of Applied Meteorology and Climatology, 2003, 42: 381-388. |
[11] | Dixon M, Wiener G. TITAN: Thunderstorm identification, tracking, analysis, and nowcasting: A radar-based methodology [J]. Journal of Atmospheric and Oceanic Technology, 1993, 10(6):785-797. |
[12] | Montanari L, Montanari A, Toth E A. A comparison and uncertainty assessment of system analysis techniques for short-term quantitative precipitation nowcasting based on radar images [Z/OL]. (2014-06-20), http://onlinelibrary.wiley.com/doi/10.1029/2005JD006729/. |
[13] | Evans A N. Cloud tracking using ordinal measures and relaxation labelling [C]//Proceedings of International Geoscience and Remote Sensing Symposium. Hamburg, Germany: IEEE, 1999: 1259-1261. |
[14] | Evans A N. On the use of ordinal measures for cloud tracking [J]. International Journal of Remote Sensing, 2000, 21(9):1939-1944. |
[15] | Kambhamettu C, Palaniappan K, Hasler A F. Hierarchical motion decomposition for cloud-tracking [C]//Proceedings of AMS 17th Conference Interactive Information and Processing Systems (IIPS) for Meteorology, Oceanography, and Hydrology. Albuquerque, USA: American Meteorological Society, 2001: 318-323. |
[16] | Bergen J R, Anandan P, Hanna K J, et al. Hierarchical model-based motion estimation [C]//Proceedings of Computer Vision-ECCV. Berlin, Germany: Springer, 1992: 237-252. |
[17] | Toklu C, Erdem A T, Sezan M I, et al. Tracking motion and intensity variations using hierarchical 2-D mesh modeling for synthetic object transfiguration [J]. Graphical Models and Image Processing, 1996, 58(6): 553-573 |
[18] | Bellerby T J. High-resolution 2-D cloud-top advection from geostationary satellite imagery [J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(12): 3639-3648. |
[19] | 王改利, 刘黎平. 多普勒雷达资料在暴雨临近预报中的应用 [J]. 气象, 2006, 31(10): 12-15.WANG Gaili, LIU Liping. Application of Doppler radar data to nowcasting of heavy rainfall [J]. Meteorological Monthly, 2006, 31(10): 12-15. (in Chinese) |