期刊检索:
  暴雨灾害   2021, Vol. 40 Issue (4): 401-409.  DOI: 10.3969/j.issn.1004-9045.2021.04.008

论文

DOI

10.3969/j.issn.1004-9045.2021.04.008

资助项目

国家重点研发计划专项项目(2017YFC1501800);中国铁路总公司科技研发项目(K2018T007);中国民用航空中南地区空中交通管理局科技项目(ZNLX202031)

第一作者

罗义, 主要从事短时临近天气预报和航空气象预报。E-mail: luoyiromon@163.com.

文章历史

收稿日期:2021-02-05
定稿日期:2021-07-25
雷达反演的多尺度风场在临近预报中的应用研究
罗义1,2 , 梁旭东1 , 王刚2 , 曹正2 , 陈春元2     
1. 中国气象科学研究院, 北京 100081;
2. 中国民用航空中南地区空中交通管理局气象中心, 广州 510405
摘要:基于多尺度IVAP(MIVAP,the Multi-scale Integrating Velocity Azimuth Process)方法反演风场进行临近外推的预报方法原理可以概括为:首先采用大尺度的分析单元通过IVAP(Integrating Velocity Azimuth Process)方法得到反演风场作为引导气流,表示雷暴的系统性移动,再通过小尺度的分析单元来反演风场,作为雷暴内部移动信息的表征,最后将不同尺度的风场进行合成得到外推的运动矢量场以进行临近预报。相比传统的多尺度TREC(MTREC,the Multi-scale Tracking Radar Echoes by Correlation)方法依赖于追踪回波过去的变化,需要相邻一个或者几个时次的雷达回波来跟踪雷达回波的移动,使用MIVAP方法外推预报只需要最新时刻的雷达反射率因子和径向速度。通过个例分析和统计分析,并且与MTREC方法进行对比试验,结果表明:MIVAP方法反演风场用于临近外推预报具有可行性,MIVAP方法具有实际应用的意义;MIVAP方法得到的外推运动矢量场方向要比MTREC方法方向总体整体偏右,且外推的平均速度要比MTREC方法快;MIVAP方法在30 min以内前几个时次的预报效果略低于MTREC方法,但是30 min以后的预报评分明显优于MTREC方法,MIVAP方法的检验评分整体要优于MTREC方法。
关键词多普勒雷达    MIVAP方法    MTREC方法    临近预报    
A study of the multi-scale winds from radar retrieval for nowcasting
LUO Yi1,2 , LIANG Xudong1 , WANG Gang2 , CAO Zheng2 , CHEN Chunyuan2     
1. State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Science, Beijing, 100081;
2. Meteorological Center of Central South Air Traffic Management Bureau of CAAC, Guangzhou, 510405
Abstract: A new method to obtain multiple scale wind retrievals is proposed by using the Multi-scale Integrating Velocity Azimuth Process (MIVAP) method with different scale analysis sectors for nowcasting. The retrieved winds from the IVAP method with the large-scale analysis sectors are considered as the steering winds involving the systematic movement, while the small-scale internal motion vectors are estimated by the IVAP method with small scale analysis sectors. Then the retrieved winds with different scales are synthesized as the motion vectors for extrapolation nowcasting. Compared to the traditional Multi-scale Tracking Radar Echoes by Correlation (MTREC) methods, which rely on tracking radar echoes of a few time slices, the MIVAP method only needs the radar echo and radial velocity observations at the latest time. Comparison experiments based on a series of convective rainfall cases demonstrate that it is practical to apply the retrieved winds from the MIVAP method for nowcasting. The statistical verification tests show that extrapolated motion from the MIVAP method appears at the right of that from the MTREC method, and the speed of the former is faster than the latter. Nowcasting using the MTREC method performs better than that applying the MIVAP method in the few periods of the first 30 minutes, while the MIVAP method achieves better evaluation scores than the MTREC method after 30 minutes of forecasting periods, so the MIVAP method has a better performance than the MTREC on the whole.
Key words: Doppler radar    MIVAP method    MTREC method    nowcasting    
引言

多普勒天气雷达是探测中小尺度天气系统的有力工具,能够提供高时空分辨率的反射率、径向速度和径向速度谱宽的资料。在目前业务系统中,特别对于暴雨强对流天气的临近预报业务,外推预报依然是对流临近预报技术的重要基础(张家国等,2008苏华英等,2020)。国内外不少学者对于对流系统的临近预报也做了非常详细的研究和总结(Wilson等,1998杨洪平等,2007)。

目前,对于雷达反射率因子资料在临近外推预报中的应用主要以追踪方法为基础,可以分为两类:第一类方法是质心法,主要通过雷暴单体质心的识别与跟踪技术,基于2D或者3D反射率资料来识别、跟踪、预报雷暴的位置和降水分布,相关的反射率资料包括最大反射率、回波顶高、垂直液态水含量(VIL),特别是质心在单体追踪方法扮演着重要的角色。早期典型的质心法WSR-88D风暴序列算法(Richard et al,1998) 在单体追踪和回波合并分裂的基础上发展而来,对于对流单体识别跟踪预报非常有效,但是对于多单体雷暴和雷暴带效果不佳。Johnson等(1998)在WSR-88D风暴序列算法的基础上,提出了改进的风暴识别追踪(SCIT,the Storm Cell Identification and Tracking)算法,通过预设多个不同的反射率阈值而不是像WSR-88D算法那样采用单一反射率阈值来进行雷暴识别,不仅能够识别孤立的雷暴,而且能够识别包含多个单体核的簇状或线状雷暴。Handwerker (2002)综合了单体识别、合并和分裂的算法,提出了雷达体扫数据跟踪对流单体(TRACE-3D)的方法,将质心法从2D平面二维拓展为3D空间的追踪算法。此外,Dixon和Wiener (1993)结合TREC (the Tracking Radar Echoes by Correlation)方法、质心法提出了雷暴识别跟踪分析和预报(TITAN,the Thunderstorm Identification,Tracking,Analysis and Nowcasting)算法(Han et al., 2009),不仅能够识别对流单体,处理对流单体的合并和分离,还能实时跟踪对流单体,同时能基于过去的趋势分析对流单体的发展和消亡。

另一类方法是交叉相关法TREC (Rinehart和Carvey,1978)。其原理是假定空间散射粒子完全由风力驱动,通过对两个连续时次体扫描的雷达回波进行相关分析,估计出运动矢量场来对降水回波进行移动预报。TREC方法提出以后,在实际应用中得到了极大发展,主要改进方法包括:COTREC (the Continuity of TREC)和MTREC (the Multi-scale Tracking Radar Echoes by Correlation)。Li等(1995)提出的COTREC方法以二维连续性方程作为强约束条件,利用观测距离和计算位移构造价值函数,通过变分技术平滑速度场,相当于对传统交叉相关法反演得到的风场进行水平无辐散处理,既改进了与地形相关的降水的临近预报效果,也一定程度上改进了强风暴的外推预报效果。但是这种方法也存在缺陷,由于采用水平无辐散约束限制,使得外推反演得到的风场在不同程度上受到削弱,从而外推得到的回波略慢于实况观测的回波(陈雷等,2009)。Wang等(2013)基于交叉相关原理提出的多尺度雷达回波跟踪和预报(MTREC)方法则提供了另一种新的思路:通过分析不同尺度的回波移动速度,将大尺度的系统性的移动速度和小尺度的雷暴内部的移动速度进行合成,综合得到外推运动场。此外,光流法(Horn和Schunck,1981Lucas和Kanade,1981Bowler et al., 2004)通过计算雷达回波的光流场来得到回波的运动矢量场,从而进行外推临近预报,其效果跟TREC方法有一定的类似性(曹春燕等,2015)。

与此同时,对于多普勒雷达径向速度资料的研究主要集中风场信息的反演方面。目前根据单多普勒雷达径向速度进行反演风场的主要方法有:VAD(Velocity Azimuth Display)方法(Lhermitte,1961),VARD(Velocity Area Display)方法(Easterbrook,1975),VVP(Velocity Volume Processing)方法(Waldteufel and Corbin, 1979Koscielny et al., 1982),UW(Uniform Winds)方法(Persson and Anderson, 1987),VAP (Velocity Azimuth Processing)方法(陶祖钰,1992),VPP (Velocity Plan Processing)方法(郎需兴等,2001),IVAP (Integrating Velocity Azimuth Process)方法(Liang,2007王欣颖和梁旭东,2009)等。其中IVAP方法对于常用的UW方法、VAD方法、VAP方法具有兼容性,具有过滤短波的特性,且分析单元的大小可以根据精度或分辨率进行灵活调整(Liang,2007)。此外,基于双多普勒或多多普勒雷达技术(Armijo,1969Brandes,1977Kessinger et al., 1987Shapiro and Mewes, 1999)或者采用变分技术结合物理约束(Shapiro et al., 2003Snyder and Zhang, 2003Tong and Xue, 2005李红莉和王叶红,2007Xiao et al., 2008Shapiro et al., 2009),进行反演二维或者三维风场的方法也得到大量研究。这些方法采用变分的方法,通过利用反射率守恒、运动方程甚至数值模式作为约束条件,以得到更加精确的风场。但是这些方法计算过于复杂且计算量大,目前还难以应用于临近预报。

总的来说,目前雷达临近外推预报主要利用反射率资料进行追踪算法来获得外推的运动矢量场,而如何利用径向速度资料获得反演风场作为外推运动矢量场的研究则较少。早期有研究对于雷暴的移动与环境平均风场的关系进行了统计分析,为利用环境风场外推预报提供了有力的可预报性依据。例如Wilson (1966)发现雷暴的回波移动都与对流的尺度相关,较小尺度的对流系统随着平均风移动,而尺度较大的对流系统移动比平均风速慢且偏向平均风向的右侧。Andersson和Ivarsson (1991)在概率临近预报模式中采用850 hPa风场作为外推风场,并考虑了利用多普勒雷达,通过VAD、UW等反演风场或者类似TREC反演风场作为引导风场进行外推预报的可能性。Turner等(2004)开发的临近预报系统(MAPLE, the McGill algorithm for precipitation nowcasting by Lagrangian extrapolation)采用结合交叉相关和变分分析的反演风场作为外推运动矢量。刘红艳等(2015)曾利用VAD方法得到不同高度的反演风场作为雷达回波的引导风场,对不同高度上的回波进行了外推实验。

为了进一步研究多尺度的雷达反演风场在临近外推预报中的应用,本文在IVAP方法的基础上,通过选用不同大小的分析单元,得到不同尺度的反演风场,再将多尺度风场进行合成得到外推运动矢量场,从而外推临近预报。通过实验表明,基于多尺度IVAP方法(MIVAP,the Multi-scale Integrating Velocity Azimuth Process)反演风场进行临近外推的预报方法具有可行性,提供了一种新的临近预报方法和思路。

1 数据和质量控制

本文所使用的雷达资料来自广州市气象局的S波段雷达,采样扫描频率为6 min,体扫模式为VCP21,包含9个仰角(仰角由低到高,即0.5°, 1.5°, 2.4°, 3.4°,4.3°,6.0°, 9.9°, 14.6°, 19.5°)。本文采用2016年4月9日的飑线过程作为个例分析,来展示MIVAP方法的原理和预报效果。最后通过统计2016年前汛期4-6月一共15个例来评估MIVAP方法的总体效果和实际应用能力。

雷达反射率因子资料的预处理和质量控制主要有:(1) 去除地物回波;(2) 中值滤波;(3) 线性插值将雷达坐标系数据插值到直角坐标系。雷达径向速度资料预处理和质量控制主要有:(1) 去除地物回波;(2) 中值滤波;(3) 距离去折叠;(4) 速度退模糊;(5) 线性插值将雷达坐标系数据插值到直角坐标系。经过预处理和质量控制后得到水平分辨率为1 km × 1 km,垂直分辨率为1 km的雷达反射率因子和径向速度资料。本文采用3 km高度的反射率作为对流系统的表征,考虑到仰角的影响,通常3 km高度的数据在230 km半径范围内有效,所以本文分析的天气个例发生在雷达中心230 km以内。

2 方法和方案设计 2.1 MTREC方法介绍

TREC方法的基本原理是将反射率图像划分为一系列固定大小的区域,通过计算连续时次雷达回波不同区域的空间相关,搜索与目标区域具有最优相关系数的区域,从而得到回波的运动矢量场。其相关系数可以表示为

$ R(m, n)=\frac{\sum\limits_{i} \sum\limits_{j}[Z(i+n, j+m, t+\Delta t)-\bar{Z}(t+\Delta t)][Z(i, j, t)-\bar{Z}(t)]}{\sqrt{\sum\limits_{i} \sum\limits_{j}[Z(i+n, j+m, t+\Delta t)-\bar{Z}(t+\Delta t)]^{2} \sum\limits_{i} \sum\limits_{j}[Z(i, j, t)-\bar{Z}(t)]^{2}}} $ (1)

其中Z(i, j, t) 表示t0时刻横坐标i纵坐标为j的反射率值,m表示从t0t1时间间隔ΔtX轴方向移动的格点数,n表示从t0t1时间间隔ΔtY轴方向移动的格点数。通过在搜索半径内改变mn的值,求得最大相关系数Rmax对应的位置,作为目标回波在t0t1间隔内所移动的位置,从而得到移动速度。

MTREC方法(Wang et al., 2013)是在TREC方法的基础上,先采用跟踪区域相对较大的TREC方法来获得大尺度的系统性的移动信息,作为雷暴系统移动的表征,然后再采用跟踪区域相对较小的TREC方法来估计小尺度的雷暴内部的运动信息。主要步骤分为4步:(1) 采用跟踪区域较大的TREC方法得到t0时刻到t1时刻的空间分辨率较低的系统性运动矢量场;(2) 通过已知的系统性运动矢量场,将t0时刻的雷达回波外推t1时刻,得到t1时刻预报场;(3) 通过对比t1时刻预报场和t1时刻观测场,采用跟踪区域尺度较小的TREC方法再次进行交叉相关计算,得到高分辨率的雷暴内部运动矢量场;(4) 将第一步得到的系统性运动矢量场和第三步得到的雷暴内部运动矢量场进行合成,从而得到用于外推预报的运动矢量场。

MTREC方法将回波运动分为大尺度和小尺度两种不同尺度的移动信息,对应取不同大小的跟踪区域。本文参考陈明轩等(2007)对于TREC方法计算参数的分析和统计,并结合当地经验,将MTREC方法设置如下:采用3 km高度的反射率作为对流回波的表征,MTREC方法的大尺度跟踪区域的大小为50 km × 50 km,跟踪区域间隔为25 km,小尺度的跟踪区域取25 km × 25 km,跟踪区域间隔为10 km。为有效降低噪声等低阈值回波对TREC算法的影响,通常需要设定最低阈值,本文最低阈值采用15 dBz;同时,每个跟踪区域中有效的数据格点数占跟踪区域内总格点数的百分率低于一定阈值时,则该“区域”将不进行TREC方法计算,本文将有效数据格点占百分率的阈值设置为60%。

2.2 MIVAP方法介绍

根据IVAP方法原理(Liang,2007),在均匀风场条件假设下,在分析单元[θ1, θ2],[r1, r2] 内,局地的平均反演风场公式可以表达为

$ u = \frac{{\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} V_R^*\cos \theta {\rm{d}}\theta {\rm{d}}r\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\sin }^2}\theta {\rm{d}}\theta {\rm{d}}r - \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} V_R^*\sin \theta {\rm{d}}\theta {\rm{d}}r\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \sin \theta \cos \theta {\rm{d}}\theta {\rm{d}}r} } } } }}{{\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\cos }^2}\theta {\rm{d}}\theta {\rm{d}}r} \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\sin }^2}\theta {\rm{d}}\theta {\rm{d}}r} - \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \sin \theta {\rm{d}}\theta {\rm{d}}r} \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \sin \theta \cos \theta {\rm{d}}\theta {\rm{d}}r} }} $ (2)
$ v = \frac{{\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} V_R^*\sin \theta {\rm{d}}\theta {\rm{d}}r\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\cos }^2}\theta {\rm{d}}\theta {\rm{d}}r - \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} V_R^*\cos \theta {\rm{d}}\theta {\rm{d}}r\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \cos \theta \sin \theta {\rm{d}}\theta {\rm{d}}r} } } } }}{{\int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\sin }^2}\theta {\rm{d}}\theta {\rm{d}}r} \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} {{\cos }^2}\theta {\rm{d}}\theta {\rm{d}}r} - \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \sin \theta \cos \theta {\rm{d}}\theta {\rm{d}}r} \int\limits_{{r_1}}^{{r_2}} {\int\limits_{{\theta _1}}^{{\theta _2}} {} \sin \theta \cos \theta {\rm{d}}\theta {\rm{d}}r} }} $ (3)

公式(2)、(3)中,VR*表示雷达观测到的径向速度,θ表示方位角,r表示距离,uv表示在分析单元[θ1, θ2],[r1, r2] 中心点处的反演风场。

MIVAP方法利用多尺度合成原理,在IVAP方法的基础上,先采用分析单元相对较大的IVAP方法来获得大尺度的系统性移动信息,然后通过采用分析单元相对较小的IVAP方法来反演小尺度的雷暴内部运动信息。MIVAP方法的主要特点是:不依赖跟踪雷达回波过去的移动,只需要最新时刻的雷达回波和径向速度。为保证实验对比的一致性,在直角坐标系下,MIVAP方法采用与MTREC方法相同的分析单元大小,即大尺度分析单元为50 km × 50 km,小尺度的分析单元取25 km × 25 km。

MIVAP方法具体步骤主要分为4步:

(1) 采用相对较大的分析单元,由原始的径向速度场VR*通过IVAP方法反演公式得到空间分辨率较低的反演风场u0v0作为大尺度运动场。本文以2016年4月9日21∶00 (世界时,下同)飑线个例的探测资料为例,图 1a为该时刻飑线的雷达反射率,整体呈西南-东北走向分布,偏东方向移动。图 1b中对应填色图为该时刻的原始径向速度场VR*,而图 1b中对应矢量图为采用大尺度的分析单元反演得到的大尺度运动场u0v0。由于采用比较大的分析单元,所以大尺度运动场u0v0比较平滑,且与飑线的移动比较一致。

图 1 2016年4月9日21:00 3 km高度回波(a, 阴影, 单位: dBz)、3 km原始径向速度(阴影, 单位: m·s-1)叠加大尺度反演运动场(b, 矢量, 单位: m·s-1)、3 km小尺度径向速度(阴影, 单位: m·s-1)叠加小尺度反演运动场(c, 矢量, 单位: m·s-1)、3 km高度回波(阴影, 单位: dBz)叠加合成运动场(d, 矢量, 单位: m·s-1) Fig. 1 (a) Radar reflectivity at 3 km level (shaded, unit: dBz), (b) retrieved motion vectors of large scale (vector, unit: m·s-1) overlaid original radial velocity (shaded, unit: m·s-1) at the 3 km level, (c) retrieved motion vectors of small scale (vector, unit: m·s-1) overlaid small scale radial velocity (shaded, unit: m·s-1) at the 3 km level, and (d) synthesizing motion vectors (vector, unit: m·s-1) overlaid radar reflectivity at the 3 km level (shaded, unit: dBz) at 21:00 UTC 9 April, 2016.

(2) 通过反演风场u0v0得到反演风场的径向速度VR0*,再由VR1* = VR*-VR0*得到新的径向速度场VR1*,定义VR1*为小尺度的径向速度,如图 1c中填色图为消除掉大尺度反演径向速度VR0*得到的新径向速度场VR1*VR1*包含着将大尺度运动场信息进行分离后所得到的小尺度运动信息。

(3) 采用分析单元相对较小的IVAP方法再次进行计算,得到高分辨率的反演风场u1v1作为小尺度运动场,如图 1c中矢量图,小尺度运动场u1v1作为表征小尺度的雷达内部运动。包含小尺度信息的运动场u1v1相比大尺度的运动场u0v0量级由50 m·s-1变为5 m·s-1,说明将大尺度运动场与小尺度运动场的平均速度不一致,同时小尺度运动场u1v1的整体方向相对大尺度运动场u0v0更加偏南,说明大尺度的系统性运动与小尺度的雷暴内部运动存在一定的差异。

(4) 将第(1)步得到的系统性运动矢量场u0v0和第(3)步得到的雷暴内部的运动矢量场u1v1进行合成,如图 1d为对应大尺度运动场和小尺度运动场合成得到的矢量图。最后得到的合成运动场在大尺度的系统性运动基础上兼顾小尺度的雷暴内部运动,将作为雷达回波移动的表征,用于外推临近预报。

为了保证外推运动矢量场的连续性,需要对IVAP反演风场进行质量控制,主要设置包括:(1) 在反演风场时,当每个分析单元内径向速度的数据格点数占分析单元内总格点数的百分率低于阈值60 %时,将不进行IVAP反演计算;(2) 当一个格点的IVAP反演风场的风向与周围平均风场的风向角度偏差超过20°时,该点风场使用平均风场进行替代;(3) 在对大尺度反演风场和小尺度反演风场进行合成时,当某一格点合成矢量的方向与该点的大尺度反演风场的方向相差超过20°时,小尺度风场将被忽略,该点采用大尺度反演风场替代合成矢量场;(4) 采用Cressman权重插值对反演风场进行连续性平滑(王改利等,2007)。

此外,应用雷达反演风场外推必须考虑的问题之一是,确定具体哪一高度上的风场最合适进行外推预报。早期基于个例统计的研究表明,对流系统的移动和平均环境风具有一定的相关性,根据以往关于对流天气移动方向与环境风场相关性的研究,一般参考700 hPa的气流作为引导气流,风场高度大致在2.5 - 3.0 km (杨凡等,2009)。Andersson和Ivarsson (1991)考虑850 hPa (1.5 km)的模式预报风场做为引导气流外推。刘红艳等(2015)利用VAD方法反演风场进行外推预报研究时发现,预报效果最好的高度层与实际天气过程有关,一般对流系统的强度越强,引导气流的高度越高。本文根据Bunkers等(2014)的研究,采用0.5-6 km内各层高度上的风场按照不同权重进行合成,即通常低层风场赋予更高权重,高层给予较小的权重。按照不同的高度给予不同的权重,其权重公式为

$ w=\left\{\begin{array}{l} \frac{R^{2}-d_{m}^{2}}{R^{2}+d_{m}^{2}}, d_{m}^{2} <R^{2} \\ 0, \quad d_{m}^{2} <R^{2} \end{array}\right. $ (4)

其中R为影响半径取6 km,d为垂直距离,这里取任意雷达资料网格点到0.5 km水平高度的距离。

3 个例分析

本节选取2016年4月9日的飑线过程进行个例分析,同时使用MTREC方法作为对比试验,来展示MIVAP方法的特性和实际预报能力。

参照第二部分MIVAP方法和MTREC方法的参数设置,对应得到两种方法的外推运动矢量场如图 2。对比图 2ab可以看出:MTREC方法运动场整体偏东移动,具有一定偏北的分量,MIVAP方法运动场则整体偏东南运动,具有明显的偏南分量。由该个例的过程分析表明,飑线在实际偏东移动过程中,不断向南发展,所以MIVAP方法运动场一定程度上更接近飑线的实际移动。

图 2 2016年4月9日21:00 3 km高度回波(阴影, 单位: dBz)叠加不同反演方法得到运动场(矢量, 单位: m·s-1): (a) MTREC方法, (b) MIVAP方法 Fig. 2 Retrieved motion vectors (vector, unit: m·s-1) overlaid radar reflectivity (shaded, unit: dBz) at the 3 km level at 21:00 UTC 9 April, 2016: (a) MTREC, and (b) MIVAP.

60 min的预报结果如图 3所示,图 3a为初始场,图 3b为60 min后的实况场,MIVAP方法(图 3d)和MTREC方法(图 3c)预报结果都能够保持飑线的整体形状和强度,但是MTREC方法60 min后的预报场对比实况场整体偏慢(图 3c),而MIVAP方法的外推结果相比MTREC方法的预报结果明显移动更快,且其预报结果跟实况匹配更好,特别是在40 dBz以上的区域,MIVAP方法得到的预报场与实况场更为接近(图 3d)。

图 3 2016年4月9日21:00雷达实况回波(a, 阴影, 单位: dBz), 22:00雷达实况回波(b, 阴影, 单位: dBz), MTREC方法外推60 min预报(c, 阴影, 单位: dBz), MIVAP方法外推60 min预报(d, 阴影, 单位: dBz) (反射率40 dBz以上回波落区, 蓝色线为初始场, 黑色线为实况场, 红色线为预报场) Fig. 3 (a) Observed reflectivity at 21:00 UTC (shaded, unit: dBz), (b) observed reflectivity at 22:00 UTC (shaded, unit: dBz), (c) 60-minute nowcasts of MTREC (shaded, unit: dBz), and (d) 60-minute nowcasts of MIVAP (shaded, unit: dBz) at 21:00 UTC 9 April, 2016. The contours enveloped the area with reflectivity over 40 dBz, the blue contour represents initial field, the black contour represents observed field, and the red contour represents forecasting field.

为了定量评估两种方法外推的预报效果,采用Dixon和Wiener (1993)所使用列联表法,利用命中率(POD)、虚假警报率(FAR)和临界成功指数(CSI)来对两种方法进行定量检验评估。将预报数据和雷达实况数据逐个格点对比:设定检验评分阈值为35 dBz,如果实况格点数据和预报格点数据都大于该阈值,则判定该格点预报成功(格点数记为n1);若实况格点数据大于阈值而预报格点数据小于阈值,则判定该格点空报失败(格点数记为n2);若实况格点数据阈值小于阈值而预报的格点数据大于阈值,则该格点是虚假警报(格点数记为n3)。命中率(POD)、虚假警报率(FAR)和临界成功指数(CSI)的计算公式如下

$ \begin{aligned} P O D=\frac{n_{1}}{n_{1}+n_{2}} \end{aligned} $ (5)
$ F A R=\frac{n_{2}}{n_{1}+n_{2}} $ (6)
$ C S I=\frac{n_{1}}{n_{1}+n_{2}+n_{3}} $ (7)

预报每隔6 min进行一次输出,以便于实况和预报进行比对,最后得到检验评分,如图 4所示。从预报评分来看:两者的预报效果都是随着预报时间延长而递减,MIVAP方法的命中率(POD)和临界成功指数(CSI) 随预报时效降低,虚假警报率(FAR)随着预报时效升高,但是MIVAP方法的检验评分明显整体优于MTREC方法。

图 4 检验评分阈值为35 dBz时外推预报60 min评分(黑色线为POD, 红色线为FAR, 蓝色线为CSI, 实线为MTREC方法, 虚线为MIVAP方法) Fig. 4 Evaluation scores with threshold value 35 dBz: POD (black line), FAR (red line) and CSI (blue) within 60-minute forecast for MTREC (solid line) and MIVAP (dashed line).
4 统计分析

为了进一步对MIVAP方法的预报效果进行评估,本文选取了2016年华南地区广州雷达在4-6月份观测到的15个个例进行检验评分的统计,以探讨MIVAP方法应用于实际业务系统的可行性。

15个雷达个例的类型包括:飑线(SQL)、中尺度对流复合体(MCS)、多单体对流系统(MC)、混合降水(Mixed)等。按照第三部分个例研究的步骤,检验评分阈值依然设置为35 dBz,15个个例的发生日期,雷暴特征和60 min的检验评分如表 1所示。

表 1 个例发生日期、雷暴类型和不同预报方法的检验评分 Table 1 List of the cases with date, storm type and mean evaluation using the MTREC and MIVAP methods for 60-minute nowcasting.

15个例按照PODFARCSI三个指标对MTREC方法和MIVAP方法进行评估。当MTREC方法或MIVAP方法在三个指标中取获得2个或以上更好的评分时,则评估认定该方法的预报效果更好。经过统计如表 1,15个个例当中,有10个个例MIVAP方法检验评分占优,5个个例MTREC方法表现更好。

60 min内预报评分随预报时效变化的平均值统计如图 5所示,从图中可知,MIVAP方法在30 min以内前几个时次的预报效果略低于MTREC方法,但是30 min以后的预报评分明显优于MTREC方法,总体来看,MIVAP方法的检验评分整体要优于MTREC方法。

图 5 检验评分阈值为35 dBz时外推预报60 min共15个例的评分平均值(黑色线为POD, 红色线为FAR, 蓝色线为CSI, 实线为MTREC方法, 虚线为MIVAP方法) Fig. 5 Mean evaluation scores with threshold value 35 dBz based on 15 cases: POD (black line), FAR (red line) and CSI (blue) within 60-minute forecast for MTREC (solid line) and MIVAP (dashed line).

将两种方法使用的外推运动矢量场进行平均,再转化为方位角和速度,得到15个个例的外推运动平均值如表 2所示。由表 2可见,MIVAP方法外推运动场的平均方位角要比MTREC方法的方位角大,说明MIVAP方法外推的总体方向相对MTREC方法要整体偏右,而根据实际经验,通常情况下雷暴在移动过程中会向右前方法发展,所以MIVAP方法更加契合雷暴向右偏运动的经典概念模型;再对比两者的平均速度可知,MIVAP方法总体外推的速度也要比MTREC方法快。

表 2 个例发生日期、雷暴类型和不同预报方法得到的外推运动平均值 Table 2 List of the cases with date, storm type and mean motion from the MTREC and MIVAP methods.
5 结论和讨论

MIVAP方法在IVAP方法的基础上,通过径向速度反演不同尺度的风场,再通过多尺度合成反演风场得到外推运动矢量场,从而实现临近外推预报。MIVAP方法先采用大尺度分析单元得到反演风场作为雷暴的系统性移动,再通过小尺度分析单元来分析反演风场,作为雷暴内部的移动信息表征,最后将不同尺度的风场进行合成得到外推的运动矢量场。通过个例分析和统计分析,得到如下结论:

(1) MIVAP方法具有实际的应用意义,也证明了多尺度合成IVAP方法反演风场用于临近外推预报的可行性。

(2) MIVAP方法外推运动场的方向要比MTREC方法运动场方向总体上偏右,且外推的平均速度要比MTREC方法快。

(3) MIVAP方法在30 min以内前几个时次的预报效果略低于MTREC方法,但是30 min以后的预报评分明显优于MTREC方法,MIVAP方法的检验评分整体上要优于MTREC方法。

相比传统的多尺度方法MTREC方法需要相邻一个或者几个时次的雷达回波来跟踪雷达回波的移动,MIVAP方法外推预报只需要最新时刻的雷达回波和径向速度。当雷达回波剧烈变化时,MTREC方法得到的外推运动矢量场准确性和平滑性会受到一定程度的影响,所以MTREC方法这类跟踪算法对于剧烈变化的雷达回波具有一定的局限性。而MIVAP方法不依赖雷达回波过去的变化,所以能够一定程度上避免MTREC方法这类跟踪算法的局限性问题。

参考文献
曹春燕, 陈元昭, 刘东华, 等. 2015. 光流法及其在临近预报中的应用[J]. 气象学报, 73(3): 471-480.
陈雷, 戴建华, 陶岚. 2009. 一种改进后的交叉相关法(COTREC)在降水临近预报中的应用[J]. 热带气象学报, 25(1): 117-122. DOI:10.3969/j.issn.1004-4965.2009.01.015
陈明轩, 王迎春, 俞小鼎. 2007. 交叉相关外推算法的改进及其在对流临近预报中的应用[J]. 应用气象学报, 18(5): 690-701. DOI:10.3969/j.issn.1001-7313.2007.05.014
郎需兴, 魏鸣, 葛文忠, 等. 2001. 一种新的单多普勒雷达风场反演方法[J]. 气象科学, 21(4): 417-424. DOI:10.3969/j.issn.1009-0827.2001.04.005
李红莉, 王叶红. 2007. 单多普勒雷达资料在伴随模式同化系统中的应用研究[J]. 暴雨灾害, 26(3): 21-26.
刘红艳, 魏鸣, 管理. 2015. 多普勒雷达风场资料在临近预报中的应用[J]. 大气科学学报, 38(4): 483-491.
苏华英, 曾莉萍, 罗乃兴, 等. 2020. 基于雷达资料的定量降水预报外推算法试验[J]. 沙漠与绿洲气象, 14(5): 27-35.
陶祖钰. 1992. 从单Doppler速度场反演风矢量场的VAP方法[J]. 气象学报, 50(1): 81-90. DOI:10.3321/j.issn:0577-6619.1992.01.009
王改利, 刘黎平, 阮征. 2007. 多普勒雷达资料在暴雨临近预报中的应用[J]. 应用气象学报, 18(3): 388-395. DOI:10.3969/j.issn.1001-7313.2007.03.016
王欣颖, 梁旭东. 2009. 基于风向线性变化假定的单Doppler雷达风场IVAP反演技术[J]. 暴雨灾害, 28(1): 73-93.
杨凡, 孙贞, 孙琪, 等. 2009. 改进的TREC算法在奥帆赛测试期的应用[J]. 气象科学, 29(3): 408-413. DOI:10.3969/j.issn.1009-0827.2009.03.022
杨洪平, 张沛源, 程明虎, 等. 2007. 多普勒天气雷达短时预报技术研究进展[J]. 暴雨灾害, 26(2): 90-96.
张家国, 王珏, 周金莲, 等. 2008. 暴雨多普勒天气雷达回波特征分析及临近预警[J]. 暴雨灾害, 27(4): 40-43.
Andersson T, Ivarsson K I. 1991. A model for probability nowcasts of accumulated precipitation using radar[J]. Journal of Applied Meteorology, 30(1): 135-141. DOI:10.1175/1520-0450(1991)030<0135:AMFPNO>2.0.CO;2
Armijo L. 1969. A theory for the determination of wind and precipitation velocities with Doppler radars[J]. Journal of the Atmospheric Sciences, 26(3): 570-573. DOI:10.1175/1520-0469(1969)026<0570:ATFTDO>2.0.CO;2
Bowler N E, Pierce C E, Seed A. 2004. Development of a precipitation nowcasting algorithm based on optical flow techniques[J]. Journal of Hydrology, 288(1-2): 74-91. DOI:10.1016/j.jhydrol.2003.11.011
Brandes E A. 1977. Flow in severe thunderstorms observed by dual-doppler radar[J]. Monthly Weather Review, 105(1): 113-120. DOI:10.1175/1520-0493(1977)105<0113:FISTOB>2.0.CO;2
Bunkers M J, Barber D A, Thompson R L, et al. 2014. Choosing a Universal Mean Wind for Supercell Motion Prediction[J]. Journal of Operational Meteorology, 2(11): 115-129. DOI:10.15191/nwajom.2014.0211
Dixon M, Wiener G. 1993. TITAN: Thunderstorm identification, tracking, analysis, and nowcasting-A radar-based methodology[J]. Journal of Atmospheric and Oceanic Technology, 10(6): 785-797. DOI:10.1175/1520-0426(1993)010<0785:TTITAA>2.0.CO;2
Easterbrook C C. 1975. Estimating horizontal wind fields by two-dimensional curve fitting of single Doppler radar measurements[C]//Radar Meteorology Conference, 16 th, Houston, Tex, 214-219
Fulton R A, Breidenbach J P, Seo D J, et al. 1998. The WSR-88D rainfall algorithm[J]. Weather and Forecasting, 13(2): 377-395. DOI:10.1175/1520-0434(1998)013<0377:TWRA>2.0.CO;2
Han Lei, Fu Shengxue, Zhao Lifeng, et al. 2009. 3D convective storm identification, tracking, and forecasting-An enhanced TITAN algorithm[J]. Journal of Atmospheric and Oceanic Technology, 26(4): 719-732. DOI:10.1175/2008JTECHA1084.1
Handwerker J. 2002. Cell tracking with TRACE3D-A new algorithm[J]. Atmospheric Research, 61(1): 15-34. DOI:10.1016/S0169-8095(01)00100-4
Horn B K, Schunck B G. 1981. Determining optical flow[J]. Artificial Intelligence, 17(1-3): 185-203. DOI:10.1016/0004-3702(81)90024-2
Johnson J T, Mackeen P L, Arthur Witt, et al. 1998. The storm cell identification and tracking algorithm: An enhanced WSR-88D algorithm[J]. Weather and Forecasting, 13(2): 263-276. DOI:10.1175/1520-0434(1998)013<0263:TSCIAT>2.0.CO;2
Kessinger C J, Ray P S, Hane C E. 1987. The Oklahoma squall line of 19 May 1977. Part Ⅰ: A multiple Doppler analysis of convective and stratiform structure[J]. Journal of the Atmospheric Sciences, 44(19): 2840-2865. DOI:10.1175/1520-0469(1987)044<2840:TOSLOM>2.0.CO;2
Koscielny A J, Doviak R, Rabin R. 1982. Statistical considerations in the estimation of divergence from single-Doppler radar and application to prestorm boundary-layer observations[J]. Journal of applied meteorology, 21(2): 197-210. DOI:10.1175/1520-0450(1982)021<0197:SCITEO>2.0.CO;2
Lhermitte R M. 1961. Precipitation motion by pulse Doppler[C]//Proc. 9th Weather Radar Conf, Amer Meteor Soc, 218-223
Li L, Schmid W, Joss J. 1995. Nowcasting of motion and growth of precipitation with radar over a complex orography[J]. Journal of Applied Meteorology and Climatology, 34(6): 1286-1300. DOI:10.1175/1520-0450(1995)034<1286:NOMAGO>2.0.CO;2
Liang Xudong. 2007. An integrating velocity-azimuth process single-Doppler radar wind retrieval method[J]. Journal of Atmospheric and Oceanic Technology, 24(4): 658-665. DOI:10.1175/JTECH2047.1
Lucas B D, Kanade T. 1981. An iterative image registration technique with an application to stereo vision[C]//Proc. Of 7th International Joint Conference on Artificial Intelligence (IJCAI), 674-679
Persson P, Anderson T. 1987. A real-time system for automatic single-Doppler wind field analysis[C]//Mesoscale Analysis and Forecasting, ESA, 61-66
Rinehart R, Garvey E. 1978. Three-dimensional storm motion detection by conventional weather radar[J]. Nature, 273(5660): 287-289. DOI:10.1038/273287a0
Shapiro A, Mewes J J. 1999. New formulations of dual-Doppler wind analysis[J]. Journal of Atmospheric and Oceanic Technology, 16(6): 782-792. DOI:10.1175/1520-0426(1999)016<0782:NFODDW>2.0.CO;2
Shapiro A, Potvin C K, Gao J. 2009. Use of a vertical vorticity equation in variational dual-Doppler wind analysis[J]. Journal of Atmospheric and Oceanic Technology, 26(10): 2089-2016. DOI:10.1175/2009JTECHA1256.1
Shapiro A, Robinson P, Wurman J, et al. 2003. Single-Doppler velocity retrieval with rapid-scan radar data[J]. Journal of Atmospheric and Oceanic Technology, 20(12): 1758-1775. DOI:10.1175/1520-0426(2003)020<1758:SVRWRR>2.0.CO;2
Snyder C, Zhang F Q. 2003. Assimilation of simulated Doppler radar observations with an ensemble Kalman filter[J]. Monthly Weather Review, 131(8): 1663-1677. DOI:10.1175//2555.1
Tong M J, Xue M. 2005. Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model: OSS experiments[J]. Monthly Weather Review, 133(7): 1789-1807. DOI:10.1175/MWR2898.1
Turner B, Zawadzki I, Germann U. 2004. Predictability of precipitation from continental radar images. Part Ⅲ: operational nowcasting implementation (MAPLE)[J]. Journal of Applied Meteorology, 43(2): 231-248. DOI:10.1175/1520-0450(2004)043<0231:POPFCR>2.0.CO;2
Waldteufel P, Corbin H. 1979. On the analysis of single-Doppler radar data[J]. Journal of Applied Meteorology, 18(4): 532-542. DOI:10.1175/1520-0450(1979)018<0532:OTAOSD>2.0.CO;2
Wang G, Wong W, Liu L, et al. 2013. Application of multi-scale tracking radar echoes scheme in quantitative precipitation nowcasting[J]. Advances in Atmospheric Sciences, 30(2): 448-460. DOI:10.1007/s00376-012-2026-7
Wilson J W. 1966. Movement and Predictability of Radar Echoes[C]//Tech Memo ERTM-NSSL-28, National Severe Storms Laboratory, 30
Wilson J W, Crook N A, Mueller C K, et al. 1998. Nowcasting thunderstorms: a status report[J]. Bulletin of the American Meteorological Society, 79(10): 2079-2100. DOI:10.1175/1520-0477(1998)079<2079:NTASR>2.0.CO;2
Xiao Q N, Eunha Lim, Duke Jin Won, et al. 2008. Doppler radar data assimilation in KMA's operational forecasting[J]. Bulletin of the American Meteorological Society, 89(1): 39-43. DOI:10.1175/BAMS-89-1-39