气象学报  2019, Vol. 77 Issue (6): 999-1014   PDF    
http://dx.doi.org/10.11676/qxxb2019.062
中国气象学会主办。
0

文章信息

吴剑坤, 陈明轩, 秦睿, 高峰, 张宇, 闫雪瑾, 宋林烨. 2019.
WU Jiankun, CHEN Mingxuan, QIN Rui, GAO Feng, ZHANG Yu, YAN Xuejin, SONG Linye. 2019.
变分回波跟踪算法及其在对流临近预报中的应用试验
The variational echo tracking method and its application in convective storm nowcasting
气象学报, 77(6): 999-1014.
Acta Meteorologica Sinica, 77(6): 999-1014.
http://dx.doi.org/10.11676/qxxb2019.062

文章历史

2019-03-21 收稿
2019-07-05 改回
变分回波跟踪算法及其在对流临近预报中的应用试验
吴剑坤1 , 陈明轩1 , 秦睿1 , 高峰1 , 张宇2 , 闫雪瑾3 , 宋林烨1     
1. 北京城市气象研究院, 北京, 100089;
2. 广东海洋大学海洋与气象学院, 湛江, 524088;
3. 河北省气象台, 石家庄, 050021
摘要: 目前业务上0—1 h对流天气临近预报仍旧以客观外推为主,采用不同外推算法,得到雷达回波以及降水的外推临近预报。以业务应用为目标,开展了变分回波跟踪算法在强对流天气临近预报中的应用研究。利用京津冀地区的8部新一代多普勒天气雷达逐6 min雷达组网拼图资料,选取2016—2018年夏季发生在京津冀地区的18个典型对流个例,开展变分回波跟踪算法和交叉相关法的0—1 h临近预报对比试验及检验评估。与传统的交叉相关法相比,变分回波跟踪算法采用变分技术求解雷达回波运动矢量场,在计算中使用两个严格的约束条件,运用迭代法进行求解,其得到的运动矢量场更为准确。结果表明,变分回波跟踪算法优于传统的交叉相关法,得到的30、60 min内雷达回波的形状、位置及强度的外推预报和实况更接近,定量检验评分更高:(1)京津冀地区4次典型对流天气过程临近预报对比试验表明,和交叉相关法相比,变分回波跟踪算法可以更好地预报出未来1 h内雷达回波的位置、形态和强度。(2)通过对18个典型对流个例定量检验,发现当雷达回波强度阈值为35和45 dBz时,无论是30或是60 min外推预报,变分回波跟踪算法的命中率(POD)和临界成功指数(CSI)都明显高于交叉相关法,且虚警率(FAR)更低;分天气类型定量检验发现,绝大多数天气类型,变分回波跟踪算法外推预报效果优于交叉相关法。
关键词: 临近预报    中值滤波    雷达回波    变分回波跟踪算法    交叉相关法    
The variational echo tracking method and its application in convective storm nowcasting
WU Jiankun1 , CHEN Mingxuan1 , QIN Rui1 , GAO Feng1 , ZHANG Yu2 , YAN Xuejin3 , SONG Linye1     
1. Institute of Urban Meteorology, China Meteorological Administration, Beijing 100089, China;
2. College of Ocean and Meteorology, Guangdong Ocean University, Zhanjiang 524088, China;
3. Hebei Meteorological Observatory, Shijiazhuang 050021, China
Abstract: The objective extrapolation forecast is the main method for 0-1 h convective storm nowcasting. Radar echo extrapolation is performed by using radar mosaics at 6 min intervals obtained from the radar images provided by eight radars in Beijing-Tianjin-Hebei region. A comparative study of two suites of extrapolated forecasts of eighteen typical convective precipitation events occurred in Beijing-Tianjin-Hebei region from 2016 to 2018 was conducted. Compared with the method of tracking radar echoes by correlation method, the variational echo tracking method utilizes variational technique to compute the motion vector field, and uses two strict constraints to get a better motion vector field. The results indicate that the variational echo tracking method performs better in prediction of the radar echo pattern, echo location, and echo intensity at 30 and 60 min forecast lead times. (1) A comparative study of two suites of forecasts of four precipitation events in Beijing-Tianjin-Hebei region has been conducted. The results indicate that the radar echo location, the echo pattern and echo intensity produced by the variational echo tracking method are closer to observations within 1 h. (2) Quantitative evaluation of the two suites of forecasts of eighteen typical convective precipitation events has been conducted. Compared with the correlation method of tracking radar echoes, the variational echo tracking method yields higher probability of detection and the critical success index of the 30 or 60 min extrapolated forecast of echoes larger than 35 and 45 dBz, while the false alarm rate is lower. Also, a quantitative evaluation classified by the weather type indicates that the variational echo tracking method performs better than the correlation method of tracking radar echoes for most weather types.
Key words: Nowcasting    Median filter    Radar echo    Variational echo tracking method    Tracking radar echoes by correlation method    
1 引言

临近预报指的是短时间(0—6 h)内雷暴、强对流等灾害性天气的预报(俞小鼎等,2012)。临近预报方法主要有外推预报、数值预报以及概念模型预报等(Wilson,et al, 1998)。中尺度数值模式预报受到模式起转(Spin-up)过程限制,最初1 h的预报无效,所以目前业务上0—1 h对流天气临近预报仍以外推预报为主(Wilson,et al,2010俞小鼎等,2012Wang,et al,2016)。

外推预报方法是20世纪60—70年代在外推雷达回波的基础上发展起来的,主要有两类:一类是单体质心跟踪算法,另一类是区域跟踪算法。单体质心跟踪算法(Dixon,et al,1993Johnson,et al,1998)是基于三维雷暴跟踪技术,通过识别和分析雷达回波,获得雷暴单体的特征(如雷暴中心、反射率因子权重中心、雷暴顶等),然后外推这些雷暴特征的运动来做出临近预报。研究表明,这类技术仅适用于强对流风暴的跟踪和临近预报。目前中国区域跟踪算法的代表是交叉相关法(Tracking Radar Echoes by Correlation method,TREC)和光流法(Optical Flow method,OF)(郑永光等,2015吴剑坤等,2018)。交叉相关法是用来追踪雷达回波移动的一种传统算法,其核心是通过计算连续时次雷达回波不同区域的最优空间相关来确定雷达回波在过去的移动矢量特征(陈明轩等,2007王丹等,2014)。研究发现,强对流过程尤其是局地生成、强度和形状随时间变化较快的雷暴,只通过简单的计算相关系数的方法难以保证跟踪的准确性,比如运动矢量可能因为相关系数太小而被剔除,还会出现相邻运动矢量对吹的情况等,这种跟踪失败会最终影响外推预报效果(韩雷等,2008曹春燕等,2015)。近几年来,关于光流法在临近预报中的应用研究逐渐增多(张蕾等,2014曹春燕等,2015Bechini,et al,2017Chen,et al,2017),如曹春燕等(2015)应用光流法开展了华南强对流天气临近预报研究,发现光流法给出的30、60 min内雷达回波的位置、形状的外推预报和实况较为接近,其外推预报具有很好的业务指示意义。光流法通过从连续时次雷达回波图像中计算光流场,用光流场代替雷达回波运动矢量场,并基于该运动矢量场对雷达回波进行外推从而达到临近预报的目的。大量学者开展了光流法和交叉相关法的0—1 h外推预报效果的定量对比评价。结果表明,对移动型局地生成及强度和形状随时间变化较快的雷暴及锋面低槽降水,光流法优于交叉相关法;对热带系统降水尤其是台风降水,光流法基于图像间运动近似为线性的这一假设而给出的运动矢量场,并没有考虑到雷达回波的旋转性,和实况回波运动有较大的误差,而交叉相关法只考虑相关系数,受此影响较小,从而导致光流法预报效果较差,不如交叉相关法(曹春燕等,2015Chen,et al,2017)。鉴于上述两种外推预报方法存在的不足,本研究尝试一种新的外推算法——变分回波跟踪算法(Variational Echo Tracking method,VET),进行强对流天气临近预报。近年来,国际上许多学者基于连续时次的雷达组网拼图资料,采用变分回波跟踪算法得到雷达回波在过去的移动矢量特征,广泛地开展了变分回波跟踪算法在降水临近预报中的应用研究(Germann,et al,2002Bellon,et al,2010Lee,et al,2010Mandapaka,et al,2012)。Germann等(2002)基于美国雷达组网资料分析了4次不同尺度降水事件,得出基于变分回波跟踪算法的临近预报有效时效(TS评分在0.3以上)至少2.5 h;Lee等(2010)基于韩国雷达组网资料,利用变分回波跟踪算法开展了2008年夏季降水试验,同样得出基于变分回波跟踪算法的临近预报有效时效为2.5 h;Mandapaka等(2012)基于瑞士雷达组网资料,分析了2005—2010年20个夏季降水事件,得出基于变分回波跟踪算法的0—3 h降水临近预报效果要优于欧拉算法和高分辨率数值天气预报模式COSMO2(the Consortium for Small-scale Modeling model),可见基于变分回波跟踪算法的降水临近预报有很好的可预报性。

本研究选取2016—2018年京津冀地区发生的18个典型对流个例,采用变分回波跟踪算法计算得到的雷达回波运动矢量场,进行对流天气的0—1 h临近预报试验,并对试验结果进行分析和检验评估,揭示该方法相比北京市气象局现有的对流风暴自动临近预报系统(Beijing-Auto-NowCasting,BJ-ANC)外推算法(交叉相关法)的优势,评判其对对流天气0—1 h临近预报的性能。

2 变分回波跟踪算法、资料和质量控制 2.1 变分回波跟踪算法基本原理

变分回波跟踪算法最早由Laroche等(1994, 1995)提出,采用变分技术,从连续时次雷达回波资料中反演出运动矢量场。

尽管相邻两个时次的雷达反射率因子会随时间变化而发生变化,但在短时间(6 min)内,这种变化可以认为比较小,所以假设雷达回波的运动基本满足拉格朗日守恒,即雷达回波满足反射率因子在短时间内守恒。变分回波跟踪算法的核心就是采用变分法,通过极小化某一包含了反射率因子守恒项和平滑约束项的代价函数,得到其最优解即最终的雷达回波运动矢量场(uv)。文中采用的代价函数(Bellon,et al,2010)为

(1)

式中,FZ为反射率因子守恒约束项,可以用全区域(Ω)所有点(xy)的连续时次的反射率因子残差的平方总和代替,其求解表达式为

(2)

式中,ΨZ为反射率因子守恒约束项的权重系数,代表了雷达资料质量,如果雷达资料可信度低,ΨZ取较小值,反之,取较大值,为了计算方便,一般整个雷达回波区域取一个常量;uv分别是雷达回波在XY方向上的运动速度;Z(xyt0)为t0时刻的反射率因子,Z(x-uΔty-vΔtt0t)为t0t时刻的反射率因子。

式(1)中,FV为引入的平滑约束条件项,其作用是限制(uv)在空间上的可变性,防止得到的(uv)与周边值差异过大,其表达式为

(3)

式中,ΨV为平滑约束条件项的权重系数,是常量。

综上可知,代价函数(式(1))中包含两个控制变量(uv),其求解转化为寻找最优的(uv)使该代价函数值最小,即转化为全局最小化的无约束优化问题,可以采用迭代方法求其最优解

(4)

式中,XN是包含了控制变量(uv)的第N次迭代点,dN为第N次的搜索方向,由代价函数相对于控制变量(uv)的梯度计算得到,αN为第N次的步长因子,N为迭代次数。不同的步长因子αN和搜索方向dN构成了不同的迭代方法。文中的迭代选用应用较广的拟牛顿法(L-BFGS法,Liu,et al,1989),该方法具有计算稳定、收敛快、节省计算资源等优点,特别适合求解大规模无约束优化问题。具体步骤如下:(1)设定收敛条件,并给出控制变量(uv)的一个初猜值(可以设为0);(2)计算代价函数值以及代价函数相对于控制变量(uv)的梯度,按照一定的方法求得搜索方向dN;(3)确定步长因子αN,使目标值有某种意义的下降;(4)应用式(4)计算得到XN+1,如果XN+1满足收敛条件,则停止迭代,得到最优解XN+1;(5)如果不收敛,则N=N+1,转到第2步,直到得到最优解。

综上,通过变分回波跟踪算法和L-BFGS法,求解得到最终的雷达回波运动矢量场。

2.2 资料说明

资料包括京津冀地区的8部新一代多普勒天气雷达(图 1)扫描的强度场拼图资料和雷达四维变分分析系统(VDRAS)的三维风场资料。其中北京、天津、石家庄、秦皇岛、沧州和邯郸是S波段多普勒天气雷达(CINRAD/SA),张北和承德是C波段多普勒天气雷达(CINRAD/CB),均在VCP21探测模式下运行,逐6 min对9个仰角进行扫描。大量研究指出,雷达强度资料在低层存在地物杂波和由于大气超折射造成的异常传播(AP)等虚假回波,在4—5 km高度可能存在“0℃层亮带”回波,这些回波会对真实回波运动矢量的分析和跟踪产生影响,而对流活动常常出现在中、低层,所以,文中选用2.5 km高度的雷达等高平面位置显示(CAPPI)拼图资料作为研究的雷达回波场。雷达四维变分分析系统是在Sun等(1997, 1998)研究基础上,通过改进雷达资料快速更新四维变分同化(RR4DVar)技术和数值云模式进一步发展的一个可实时运行系统。该系统使用了雷达探测资料、自动气象站观测资料和中尺度数值模式结果,通过雷达四维变分分析系统,最终可快速更新模拟分析得到18 min间隔的低层大气三维动力、热动力和微物理场的精细结构(陈明轩等, 2016a, 2016b)。文中将雷达四维变分分析系统的临近预报风场资料作为变分回波跟踪算法的初猜场。

图 1  京津冀地区地形高度(色阶,单位:m),北京、天津、石家庄、秦皇岛、沧州和邯郸S波段雷达及张北和承德C波段雷达的位置(·表示),雷达四维变分分析系统模拟范围(黑色虚线方框表示) Fig. 1  Topography (shaded, unit: m) in Beijing-Tianjin-Hebei region, S-band radar sites in Beijing, Tianjin, Shijiazhuang, Qinhuangdao, Cangzhou and Handan, and C-band radar sites in Zhangbei and Chengde (denoted by black dots). The numerical simulation domain of VDRAS is denoted by black dotted box
2.3 质量控制 2.3.1 资料质量控制——中值滤波

中值滤波是常用的一种非线性平滑滤波方法,能够有效抑制脉冲噪声,经常被用于去除图像或者其他信号中的噪声(赵悦等,2007)。中值滤波的主要原理是把数字图像或数字序列中某点的值用该点一个邻域中各点值的中值代替,这里的邻域通常被称为窗口。通过中值滤波算法可以很好地对图像进行平滑处理,从而消除孤立的噪声点。文中对雷达回波采用正方形窗口滤波,具体为:用点(xy)周围(2n+1)×(2n+1)大小子窗口内的中值代替(xy)点的值f(xy)

(5)

式中,M表示取数值序列f(ij)的中值,f(ij)是点(ij)的反射率因子,文中n取2,即取(5×5)大小的窗口。

受蒙古低涡影响,2017年7月7日傍晚至夜间,京津冀地区出现了一次强飑线天气过程。选取该次过程的21时36分(北京时,下同)京津冀地区反射率因子拼图资料做中值滤波对比试验。可见,中值滤波后(图 2b),较好地保留了飑线原有的形状、强回波中心以及层状云回波特征;较中值滤波前(图 2a)去噪效果好,回波变得更加平滑,边沿更加清晰。

图 2  2017年7月7日21时36分京津冀地区反射率因子拼图中值滤波前(a)后(b)对比 Fig. 2  Comparison of radar mosaics over Beijing-Tianjin-Hebei region at 21:36 BT 7 July 2017 before (a) and after (b) the median filetering
2.3.2 算法主要参数设置

大量研究和统计分析(Germann,et al,2002Bellon,et al,2010)指出,对雷达回波中值滤波后,还有一些因素会影响变分回波跟踪算法的效果,主要包括初猜场、权重系数ΨZΨV、时间平滑等。

变分回波跟踪算法在使用迭代方法求最优解时,需要给出控制变量(uv)的初猜场。大量试验研究表明,一个较为准确的初猜场对于最终雷达回波运动矢量的精确快速求解至关重要。起初,研究人员开发了一个“尺度猜测方法”(Laroche,et al,1994),即通过逐步提高初猜场格点的分辨率,得到较为准确的初猜场,分3个步骤进行(图 3)。但由于该方法求解初猜场用时较多,不能满足业务实时运行的要求。若采用上一次的运动矢量场或者数值模式风场代替尺度猜测方法得到的初猜场,可以大幅度缩短程序运行时间(Germann,et al,2002)。雷暴的运动包括平流和传播,当环境为强气流控制时,其运动主要取决于平流。而雷暴单体主要随着承载层的平均气流方向移动,所以本研究采用雷达四维变分分析系统的临近预报三维风场资料,取其3—6 km高度的平均水平风场作为初猜场,利用变分回波跟踪算法求解雷达回波运动矢量场(图 4)。从图 4可以看出,雷达四维变分分析系统的3—6 km高度平均水平风场不仅较好地描绘了台风“安比”的逆时针螺旋状旋转的风场结构,还给出了台风向西北方向移动的大趋势,大致符合台风“安比”的环境风场特征,有较高的参考价值。

图 3  尺度猜测方法求解最终初猜场示意 (a.以(0,0)为初猜值得到1×1雷达回波运动矢量场,b.以(a)为初猜场计算得到的5×5雷达回波运动矢量场,c.以(b)为初猜场计算得到的25×25雷达回波运动矢量场)(Germann,et al,2002) Fig. 3  Scaling-guess procedure: the echo motion field is retrieved in three runs with increasing resolution (The retrieval starts with (a) a uniform field, which is used as the first guess to retrieve (b) the field on 5×5 grids, which in turn is used as the first guess for (c) the final minimization with 25×25 grids) (Germann, et al, 2002)
图 4  2018年7月24日台风“安比”过程的08时48分反射率因子拼图叠加雷达四维变分分析系统的3—6 km平均水平风场(a)和09时18分反射率因子拼图(b) Fig. 4  Radar mosaics of Typhoon Ampil on 24 July 2018 (a) Radar mosaic at 08:48 BT, and the average horizontal wind field from 3 km to 6 km by VDRAS, (b) radar mosaic at 09:18 BT

代价函数(式(1))中包含了两个权重系数ΨZΨV,其取值不同会影响L-BFGS迭代求解时的收敛性以及求得的运动矢量场是否最优(Bellon,et al,2010)。很显然,ΨV为0时,仅使用反射率因子项得到的运动矢量场会导致运动矢量场空间分布不规则。如果ΨV太大,则会削弱反射率因子项作用,导致不收敛或收敛不合理。经过大量试验,最终获得最优权重系数ΨZ=0.5和ΨV=10。

虽然变分回波跟踪算法通过加入平滑约束项限制了运动矢量场在空间上的可变性,但从时间上来看,仍旧存在一些小波动。运动矢量场在时间上的差异反映了平均运动微小的不稳定性,可能导致一些效果较差的外推预报(Bellon,et al,2010)。为了改进外推预报效果,对得到的运动矢量场做时间平滑处理,即给予当前的运动矢量场(Vc)和过去的时间平滑后运动矢量场(Vp)不同权重,从而得到当前最终运动矢量场(Vf)

(6)

式中,取Ψc=0.6,Ψp=0.4,Vp是上一个时次经过时间平滑后的雷达回波运动矢量场。

3 临近预报对比试验

对0—1 h外推预报来说,雷达回波运动矢量场的准确描述至关重要。基于变分回波跟踪算法得到的运动矢量场,不仅考虑了雷达回波的移动,还考虑了不同区域回波的形变特征。目前北京市气象局业务实时运行的自动临近预报系统以京津冀地区8部多普勒天气雷达组网拼图资料为基础,采用传统的交叉相关法得到雷达回波运动矢量,利用拉格朗日持续性预报,滚动进行60 min内间隔6 min的对流临近预报。为了揭示变分回波跟踪算法和交叉相关法的优、缺点,采用变分回波跟踪算法得到的雷达回波运动矢量场,利用同样的预报方式,得到0—1 h雷达回波外推预报,并开展这两种外推算法的临近预报对比试验。

3.1 2017年7月7日强飑线过程

受蒙古低涡影响,2017年7月7日傍晚至夜间,京津冀地区出现了一次强飑线天气过程。20时30分前后,北京西部地区开始形成东北—西南向的飑线,该飑线自西北向东南移动,8日03时前后,强回波减弱东移入海,过程结束,历时超过6 h,横扫京津冀大部分地区,造成多地区出现短时强降水、雷暴大风和冰雹等灾害性天气。京津冀地区8部多普勒天气雷达完整地探测到此次强飑线演变过程,于是基于京津冀地区2.5 km高度反射率因子拼图资料,开展变分回波跟踪算法和交叉相关法的临近预报对比试验。

对比两种外推算法对这次强飑线过程的20时54分起报的1 h外推预报(图 5cd)和21时54分的反射率因子拼图实况(图 5b)可知,利用变分回波跟踪算法的1 h外推预报得到的飑线位置、形状与雷达观测到的实况比较接近。该方法较好预报出了1 h后飑线的强回波(≥35 dBz)位置、形态和强度等特征。实况飑线上的强回波区域呈“两头小中间大”带状分布,东北部最强回波约50 dBz,西南部最强回波约55 dBz,中间最大,最强回波超过65 dBz。变分回波跟踪算法的1 h外推预报和实况相比,东北部范围略偏小,西南部略偏强,中间相当。交叉相关法也进行了较好的外推预报,和变分回波跟踪算法相比,其飑线东北部强回波强度偏弱,范围小,且落区偏西,这和交叉相关法给出的运动矢量场在该区域较弱有关(图 5a)。同时也发现,21时54分的反射率因子拼图实况显示,飑线西南部开始分裂成两个强回波中心,两种外推算法均没有预报出分裂。这是因为目前外推预报无法从物理机制上预报雷暴的生成、发展、减弱和消亡。

图 5  2017年7月7日强飑线过程反射率因子拼图 20时54分实况叠加变分回波跟踪算法(蓝色箭头表示)和交叉相关法(黑色箭头表示)的运动矢量场(a),21时54分实况(b),变分回波跟踪算法的20时54分的1 h外推预报(c)及交叉相关法的20时54分的1 h外推预报(d) Fig. 5  Radar mosaics of a strong squall line on 7 July 2017: (a) observations at the initial time of the forecast at 20:54 BT, and the motion vector fields obtained by the methods of VET (denoted by blue arrows) and TREC (denoted by black arrows), (b) observations at 21:54 BT, (c) forecast of radar mosaic at 21:54 BT predicted by the VET method, (d) forecast of radar mosaic at 21:54 BT predicted by the TREC method

从本次强飑线过程的两种外推算法的临近预报对比试验来看,和交叉相关法相比,变分回波跟踪算法更好地预报出了1 h后飑线的回波形态、位置以及强回波特征。表明利用变分回波跟踪算法对类似飑线这种组织性比较强的对流系统的1 h外推预报有较高的参考价值。

3.2 2016年8月12日副热带高压边缘型强降水过程

受冷空气与副热带高压外围偏南暖湿气流共同影响,2016年8月12日京津冀地区出现一次强降水天气过程。06时开始,北京北部地区不断有强回波生成,并向西传播,而河北西北部,有一条近似南北向对流回波以较快的速度东移,09时前后,两处强回波在北京怀柔区合并,强度和范围迅速增大,继续东移,最后发展成“逗点状”强回波(图 6a)。这次对流过程给京津冀地区造成了局地的强降水和雷暴大风天气。

图 6  2016年8月12日强降水过程反射率因子拼图 12时18分实况叠加变分回波跟踪算法(蓝色箭头表示)和交叉相关法(黑色箭头表示)的运动矢量场(a),12时48分实况(b),变分回波跟踪算法的12时18分的30 min外推预报(c)及交叉相关法的12时18分的30 min外推预报(d) Fig. 6  Radar mosaics of a heavy rainfall event on 12 August 2016 (a) observations at the initial time of the forecast at 12:18 BT, and the motion vector fields obtained by the methods of VET (denoted by blue arrows) and TREC (denoted by black arrows), (b) observations at 12:48 BT, (c) forecast of radar mosaic at 12:48 BT predicted by the VET method, (d) forecast of radar mosaic at 12:48 BT predicted by the TREC method

从2016年8月12日12时18分强降水过程反射率因子拼图(图 6a)和12时48分实况(图 6b)以及两种外推算法的12时18分起报的30 min外推预报(图 6cd)可见,雷达回波实况的半小时演变显示,“逗点状”强回波在向东移动过程中,北部向东北偏东方向移动,南部向东南偏东方向移动。利用变分回波跟踪算法得到的运动矢量场与雷达回波实况移动趋势较为一致,且预报的30 min后雷达回波的位置和强度与实况基本相符,也预报出了“逗点状”强回波形态。利用交叉相关法得到的运动矢量场显示南部强回波向东北方向移动,导致预报的30 min后南部的雷达回波位置略偏西偏北。相对而言,变分回波跟踪算法对此次强降水过程的30 min外推预报效果更优。

3.3 2018年8月11日强对流天气过程

受高空西来槽以及低层切变线共同影响,2018年8月11日下午至夜间,北京地区出现一次明显的强对流天气过程。从雷达回波演变来看,12时前后,在河北西部有弱对流回波生成,最强回波为35 dBz,以较快速度向东移动。进入北京地区后,迅速发展加强,最后发展成为较为明显的线状多单体风暴,最强回波达50 dBz(图略)。此次强对流天气过程给北京地区造成了雷暴大风以及东北部(顺义、密云和平谷等地)强降水。

对比两种外推算法17时30分起报的30 min外推预报(图 7cd)和此次强对流过程18时的反射率因子拼图实况(图 7b)发现,利用变分回波跟踪算法的30 min外推预报得到的线状多单体风暴的位置、形状与雷达观测到的实况基本相符,也基本预报出南、北两段强回波区,最大反射率因子强度均达50 dBz,只是南段强回波范围略比实况小。从两种外推算法17时30分的运动矢量场(图 7a)可以看出,交叉相关法的运动矢量场略小,方向偏北,这也是交叉相关法30 min外推预报雷达回波整体偏西偏北的原因。总体来说,变分回波跟踪算法对此次强对流过程的30 min外推预报效果略优,表明该算法对类似这种线状多单体风暴30 min外推预报有一定的参考价值。

图 7  2018年8月11日强对流天气过程反射率因子拼图 17时30分实况叠加变分回波跟踪算法(蓝色箭头表示)和交叉相关法(黑色箭头表示)的运动矢量场(a),18时实况(b),变分回波跟踪算法的17时30分的30 min外推预报(c)及交叉相关法的17时30分的30 min外推预报(d) Fig. 7  Radar mosaics of a severe convective weather process on 11 August 2018: (a) observations at the initial time of the forecast at 17:30 BT, and the motion vector fields obtained by the methods of VET (denoted by blue arrows) and TREC (denoted by black arrows), (b) observations at 18:00 BT, (c) forecast of radar mosaic at 18:00 BT predicted by the VET method, (d) forecast of radar mosaic at 18:00 BT predicted by the TREC method
3.4 2018年7月24日台风“安比”过程

台风“安比”于2018年7月18日在西北太平洋洋面生成,22日在上海登陆后一路北上,24日03时由沧州进入河北,17时由承德移出河北并继续北上减弱。受台风“安比”影响,京津冀东部地区(包括北京东部地区、天津以及河北东部)出现大范围暴雨,局地大暴雨天气。

从24日09时48分雷达观测到的反射率因子拼图实况(图 8b)上可以看出,台风眼区回波未完全闭合,眼区外是大范围的螺旋雨带。回波大值区主要集中在台风中心西侧和东北侧,其中西侧的强回波区覆盖了第二象限和部分第三象限,东北侧强回波覆盖了第一象限,这两个区域的最强回波超过40 dBz。而台风中心的南侧(第四象限和部分第三象限)有零星的弱回波。

图 8  2018年7月24日台风“安比”过程的反射率因子拼图 09时48分实况叠加变分回波跟踪算法(蓝色箭头表示)和交叉相关法(黑色箭头表示)的运动矢量场(a),09时48分实况(b),变分回波跟踪算法的08时48分的1 h外推预报(c)及交叉相关法的08时48分的1 h外推预报(d) Fig. 8  Radar mosaics of Typhoon Ampil on 24 July 2018 (a) observations at the initial time of the forecast at 08:48 BT, and the motion vector fields obtained by the methods of VET (denoted by blue arrows) and TREC (denoted by black arrows), (b) observations at 09:48 BT, (c) forecast of radar mosaic at 09:48 BT predicted by the VET method, (d) forecast of radar mosaic at 09:48 BT predicted by the TREC method

基于08时48分的实况雷达回波和运动矢量场(图 8a),开展两种外推算法的1 h外推预报(图 8cd)。可以看出,两种外推算法的运动矢量场都具有逆时针螺旋状旋转的风场结构,眼区的运动矢量小,螺旋雨带区的运动矢量大等特征;且均预报1 h后台风眼区回波未完全闭合,眼区外是大范围的螺旋雨带。预报的回波大值区也主要集中在台风中心西侧和东北侧,但西侧强回波区覆盖了整个第二和第三象限,和实况相比,范围稍大,而东北侧回波强度略小。还可以看到,两种外推算法均预报螺旋雨带最外侧有一条西南—东北向带状回波,而实况上只有很小的一块,出现较为明显的空报,这主要是因为外推预报无法预报雷暴的生消过程。总体来说,两种外推算法预报的台风螺旋雨带的整体形态、强度和位置与1 h后的实况比较接近。可见,变分回波跟踪算法用来追踪台风降水回波是可行的。

综上,选取京津冀地区4次典型对流天气过程,开展了变分回波跟踪算法和交叉相关法临近预报对比试验及评估。从结果来看,大部分对流天气过程,两种外推算法可以较好地预报出未来30、60 min内雷达回波的位置、形态和强度。总体来说,变分回波跟踪算法的外推预报效果优于交叉相关法,其外推预报结果对业务有更好的指导意义。

4 定量检验

4次典型对流天气过程的变分回波跟踪算法和交叉相关法临近预报对比试验结果表明,3次对流天气过程(2017年7月7日强飑线、2016年8月12日强降水和2018年8月11日强对流天气),变分回波跟踪算法较交叉相关法更好地预报出雷达回波的位置、形态和强度。而2018年7月24日台风“安比”过程,两种外推算法的1 h外推效果相似,这可能和台风系统降水尺度大、比较稳定,变分回波跟踪算法和交叉相关法得到的运动矢量场比较接近有关。

为了对两种外推算法的预报结果进行定量评价,引入气象领域中应用较广的三个评价指标:命中率(probability of detection,POD)、虚警率(false alarm rate,FAR)和临界成功指数(critical success index,CSI)(曹春燕等,2015Chen, et al,2017)。这些指标可以通过列联表评价方法得到。假设检验的雷达回波阈值为M,则NH表示实况不小于M且预报也不小于M的格点个数,即命中次数;NM表示实况不小于M而预报小于M的格点个数,即漏报次数;NF表示实况小于M而预报不小于M的格点个数,即空报次数。本试验使用的雷达回波阈值分别为15、25、35、和45 dBz。则评价指标计算为

(7)
(8)
(9)

选取2016—2018年京津冀地区发生的18个典型对流个例(表 1),分别为蒙古低涡型降水4个、东北冷涡型降水4个、副热带高压边缘型降水5个、西风槽型降水4个和台风型降水1个,基本上涵盖了京津冀地区夏季所有降水天气类型。对18个典型对流个例分别使用变分回波跟踪算法和交叉相关法进行了30和60 min的外推预报,并进行不同阈值下的检验评估(表 23)。该结果为18个典型个例统计时段内逐6 min检验的平均值,样本总数为1120个,网格范围为550×550。从检验结果可以看出,当雷达回波阈值为35和45 dBz时,无论是30还是60 min外推预报,变分回波跟踪算法的命中率和临界成功指数都明显高于交叉相关法,且虚警率更低。但也可以看出,当雷达回波阈值为15 dBz时,两种外推算法的效果无明显差异;25 dBz时,变分回波跟踪算法的命中率和临界成功指数略高于交叉相关法。可见,相比交叉相关法,在强回波区的外推预报方面,变分回波跟踪算法更优。

表 1  试验选取的18个典型对流个例 Table 1  18 typical convective cases selected for experimental study
个例起止时间 天气类型
2016年6月27日19时—28日06时 蒙古低涡型降水
2017年7月7日19时—8日00时 蒙古低涡型降水
2017年8月11日19时—12日02时 蒙古低涡型降水
2018年6月12日13—18时 蒙古低涡型降水
2016年6月6日16时—7日01时 东北冷涡型降水
2017年6月25日13—19时 东北冷涡型降水
2017年7月13日18—22时 东北冷涡型降水
2017年8月8日19时—9日00时 东北冷涡型降水
2016年7月24日19时—25日00时 副热带高压边缘型降水
2016年8月12日09—16时 副热带高压边缘型降水
2017年7月20日21时—21日03时 副热带高压边缘型降水
2018年7月9日09—15时 副热带高压边缘型降水
2018年7月16日14—20时 副热带高压边缘型降水
2016年8月6日20时—7日04时 西风槽型降水
2016年9月11日00—04时 西风槽型降水
2018年6月30日15—20时 西风槽型降水
2018年8月11日14—20时 西风槽型降水
2018年7月24日08—15时 台风型降水
表 2  两种外推算法30 min预测评分对比 Table 2  30 min prediction scores of two extrapolation algorithms
雷达回波阈值
(dBz)
命中率 虚警率 临界成功指数
VET TREC VET TREC VET TREC
15 0.66 0.66 0.31 0.30 0.52 0.52
25 0.61 0.60 0.36 0.36 0.46 0.45
35 0.42 0.38 0.56 0.58 0.27 0.25
45 0.27 0.21 0.72 0.77 0.16 0.12
表 3  两种外推算法60 min预测评分对比 Table 3  60 min prediction scores of two extrapolation algorithms
雷达回波阈值
(dBz)
命中率 虚警率 临界成功指数
VET TREC VET TREC VET TREC
15 0.55 0.55 0.41 0.41 0.40 0.40
25 0.48 0.45 0.49 0.50 0.33 0.31
35 0.24 0.19 0.75 0.78 0.14 0.11
45 0.09 0.05 0.90 0.94 0.05 0.03

图 9给出了雷达回波阈值为35 dBz时的变分回波跟踪算法和交叉相关法外推预报定量检验结果,图 9a是5种不同降水类型的平均,可以看出变分回波跟踪算法优于交叉相关法,30 min(60 min)外推预报的命中率是42%(24%),高交叉相关法4个百分点(5个百分点),其临界成功指数也高交叉相关法2个百分点(3个百分点),而虚警率均低于交叉相关法。5类不同降水类型下两种外推算法评价对比可以看出,对于蒙古低涡型和东北冷涡型降水天气过程(图 9bc),变分回波跟踪算法的30 min (60 min)外推预报的命中率高交叉相关法5—6个百分点(5—7个百分点),且临界成功指数高4个百分点(3—4个百分点),可见变分回波跟踪算法具有明显优势。对于副热带高压边缘型和西风槽型降水天气过程(图 9de),变分回波跟踪算法的30 min(60 min)外推预报命中率高交叉相关法2个百分点(4—5个百分点),且临界成功指数高1—2个百分点(2—3个百分点),可见变分回波跟踪算法略优于交叉相关法。对于台风型降水天气过程(图 9f),30 min外推预报,变分回波跟踪算法和交叉相关法的各项评价指标接近;60 min的外推预报,变分回波跟踪算法略优于交叉相关法,由于京津冀地区台风型降水个例较少,需要更多的个例参与检验,才能给出更准确的结果。

图 9  雷达回波阈值为35 dBz时的变分回波跟踪算法和交叉相关法临近预报试验结果的命中率、虚警率和临界成功指数对比 (a.总个例平均,b.蒙古低涡型降水,c.东北低涡型降水,d.副热带高压边缘型降水,e.西风槽型降水,f.台风型降水) Fig. 9  Comparison of the POD, FAR and CSI for forecasts of echoes greater than or equal to 35 dBz produced by VET and TREC methods (a. averages of five types of precipitation, b. precipitation induced by Mongolia Cold Vortex, c. precipitation induced by Northeast Cold Vortex, d. precipitation induced by subtropical high, e. precipitation induced by Westerly trough, f. typhoon precipitation)
5 结论与讨论

利用京津冀地区的8部新一代多普勒天气雷达逐6 min组网拼图资料,选取2016—2018年夏季发生在京津冀地区的18个典型对流个例,进行变分回波跟踪算法和交叉相关法的0—1 h临近预报对比试验及检验评估。与传统的交叉相关法相比,变分回波跟踪算法采用变分技术,使用两个严格的约束条件(反射率因子守恒约束和平滑约束)计算运动矢量场,即使雷暴的移动和外形变化比较剧烈,也可以准确地得到其整体移动趋势。试验及检验结果表明,变分回波跟踪算法优于交叉相关法,给出的30、60 min内雷达回波的形状、位置及强度的外推预报和实况更为接近,定量检验评分更高。

(1) 通过对雷达资料进行中值滤波处理,去除噪声的影响,得到更加平滑、清晰和真实的雷达回波。通过对变分回波跟踪算法的主要参数进行合理设置,得到更平滑、更准确的雷达回波运动矢量场。

(2) 通过对京津冀地区4次典型对流天气过程详细分析,发现变分回波跟踪算法可以较好地预报未来30、60 min内雷达回波的位置、形态和强度。总体来说,变分回波跟踪算法的外推预报效果要优于交叉相关法,其外推预报结果对业务有更好的指导意义。

(3) 通过对18个典型对流个例定量检验发现,当雷达回波阈值为35和45 dBz时,无论是30还是60 min外推预报,变分回波跟踪算法的命中率和临界成功指数都明显高于交叉相关法,而虚警率更低。但也可以看出,当雷达回波阈值为15 dBz时,两种外推算法的效果无明显差异;25 dBz时,变分回波跟踪算法的命中率和临界成功指数略高于交叉相关法。可见,相比交叉相关法,在强回波区的外推预报方面,变分回波跟踪算法更优。雷达回波阈值为35 dBz时两种算法外推预报定量检验结果表明,对于蒙古低涡型和东北冷涡型降水天气过程,变分回波跟踪算法优于交叉相关法,这可能是因为变分回波跟踪算法能更好地捕捉到强回波区域的整体运动趋势,而交叉相关法做不到;对于副热带高压边缘型和西风槽型降水天气过程,变分回波跟踪算法略优于交叉相关法;对于台风型降水天气过程,两种外推算法效果接近,这可能和台风系统降水尺度大、比较稳定,变分回波跟踪算法和交叉相关法得到的运动矢量场比较接近有关,且该类型降水个例较少,后期需要更多的个例参与检验才能给出更准确的结果。

需要指出的是,同一雷达回波阈值下,无论是变分回波跟踪算法还是交叉相关法,命中率和临界成功指数都随着预报时效的延长而降低,虚警率均随着时效的延长而增大,说明两种外推算法的可用性是随着时间的延长而降低的。但也可以看出,相比交叉相关法,变分回波跟踪算法给出的运动矢量场更接近实况雷达回波移动趋势,其可以弥补交叉相关法的不足,提升0—1 h对流天气临近预报、预警能力。

目前,北京市气象局已有交叉相关法和变分回波跟踪算法两种外推算法,正在研发应用较为广泛的多尺度光流法。未来的工作重点是开展三种外推算法性能对比和融合技术研发。

参考文献
曹春燕, 陈元昭, 刘东华, 等. 2015. 光流法及其在临近预报中的应用. 气象学报, 73(3): 471-480. Cao C Y, Chen Y Z, Liu D H, et al. 2015. The optical flow method and its application to nowcasting. Acta Meteor Sinica, 73(3): 471-480. (in Chinese)
陈明轩, 王迎春, 俞小鼎. 2007. 交叉相关外推算法的改进及其在对流临近预报中的应用. 应用气象学报, 18(5): 690-701. Chen M X, Wang Y C, Yu X D. 2007. Improvement and application test of TREC algorithm for convective storm nowcast. J Appl Meteor Sci, 18(5): 690-701. DOI:10.3969/j.issn.1001-7313.2007.05.014 (in Chinese)
陈明轩, 高峰, 孙娟珍, 等. 2016a. 基于VDRAS的快速更新雷达四维变分分析系统. 应用气象学报, 27(3): 257-272. Chen M X, Gao F, Sun J Z, et al. 2016a. An analysis system using rapid-updating 4-D variational radar data assimilation based on VDRAS. J Appl Meteor Sci, 27(3): 257-272. (in Chinese)
陈明轩, 肖现, 高峰, 等. 2016b. 基于雷达四维变分分析系统的强对流高分辨率模拟个例分析和批量检验. 气象学报, 74(3): 421-441. Chen M X, Xiao X, Gao F, et al. 2016b. A case study and batch verification on high resolution numerical simulations of severe convective events using an analysis system based on rapid-refresh 4-D variational radar data assimilation. Acta Meteor Sinica, 74(3): 421-441. (in Chinese)
韩雷, 王洪庆, 林隐静. 2008. 光流法在强对流天气临近预报中的应用. 北京大学学报(自然科学版), 44(5): 751-755. Han L, Wang H Q, Lin Y J. 2008. Application of optical flow method to nowcasting convective weather. Acta Scient Natur Univ Pekin, 44(5): 751-755. DOI:10.3321/j.issn:0479-8023.2008.05.014 (in Chinese)
王丹, 王改利, 刘黎平, 等. 2014. 基于雷达回波外推和中尺度模式预报的短时降水对比分析. 高原气象, 33(3): 811-822. Wang D, Wang G L, Liu L P, et al. 2014. Comparisons analysis on short-term precipitation between the radar-based extrapolation and the meso-scale numerical model weather prediction. Plateau Meteor, 33(3): 811-822. (in Chinese)
吴剑坤, 陈明轩. 2018. 基于雷达回波区域跟踪算法的临近预报技术进展. 气象科技, 46(5): 899-909. Wu J K, Chen M X. 2018. Advances in nowcasting techniques based on radar echo area tracking algorithm. Meteor Sci Technol, 46(5): 899-909. (in Chinese)
俞小鼎, 周小刚, 王秀明. 2012. 雷暴与强对流临近天气预报技术进展. 气象学报, 70(3): 311-337. Yu X D, Zhou X G, Wang X M. 2012. The advances in the nowcasting techniques on thunderstorms and severe convection. Acta Meteor Sinica, 70(3): 311-337. DOI:10.3969/j.issn.1004-4965.2012.03.003 (in Chinese)
张蕾, 魏鸣, 李南, 等. 2014. 改进的光流法在回波外推预报中的应用. 科学技术与工程, 14(32): 133-137, 148. Zhang L, Wei M, Li N, et al. 2014. Improved optical flow method application to extrapolate radar echo. Sci Technol Eng, 14(32): 133-137, 148. DOI:10.3969/j.issn.1671-1815.2014.32.027 (in Chinese)
赵悦, 陈家华, 章建军, 等. 2007. 基于中值滤波和小波变换的天气雷达回波图像处理. 气象科学, 27(1): 63-68. Zhao Y, Chen J H, Zhang J J, et al. 2007. Weather radar echo image processing based on median filter and wavelet transform. Scientia Meteor Sinica, 27(1): 63-68. DOI:10.3969/j.issn.1009-0827.2007.01.009 (in Chinese)
郑永光, 周康辉, 盛杰, 等. 2015. 强对流天气监测预报预警技术进展. 应用气象学报, 26(6): 641-657. Zheng Y G, Zhou K H, Sheng J, et al. 2015. Advances in techniques of monitoring, forecasting and warning of severe convective weather. J Appl Meteor Sci, 26(6): 641-657. (in Chinese)
Bechini R, Chandrasekar V. 2017. An enhanced optical flow technique for radar nowcasting of precipitation and winds. J Atmos Ocean Technol, 34(12): 2637-2658. DOI:10.1175/JTECH-D-17-0110.1
Bellon A, Zawadzki I, Kilambi A, et al. 2010. McGill Algorithm for precipitation nowcasting by Lagrangian extrapolation (MAPLE) applied to the South Korean radar network. Part Ⅰ:Sensitivity studies of the variational echo tracking (VET) technique. Asia-Pacific J Atmos Sci, 46(3): 369-381.
Chen Y Z, Lan H P, Chen X L, et al. 2017. A nowcasting technique based on application of the particle filter blending algorithm. J Meteor Res, 31(5): 931-945. DOI:10.1007/s13351-017-6557-9
Dixon M, Wiener G. 1993. TITAN:Thunderstorm identification, tracking, analysis, and nowcasting:A radar-based methodology. J Atmos Ocean Technol, 10(6): 785-797. DOI:10.1175/1520-0426(1993)010<0785:TTITAA>2.0.CO;2
Germann U, Zawadzki I. 2002. Scale-dependence of the predictability of precipitation from continental radar images. Part Ⅰ:Description of the methodology. Mon Wea Rev, 130(12): 2859-2873.
Johnson J T, MacKeen P L, Witt A, et al. 1998. The storm cell identification and tracking algorithm:An enhanced WSR-88D algorithm. Wea Forecasting, 13(2): 263-276. DOI:10.1175/1520-0434(1998)013<0263:TSCIAT>2.0.CO;2
Laroche S, Zawadzki I. 1994. A variational analysis method for retrieval of three-dimensional wind field from single-Doppler radar data. J Atmos Sci, 51(18): 2664-2682. DOI:10.1175/1520-0469(1994)051<2664:AVAMFR>2.0.CO;2
Laroche S, Zawadzki I. 1995. Retrievals of horizontal winds from single-Doppler clear-air data by methods of cross correlation and variational analysis. J Atmos Ocean Technol, 12(4): 721-738. DOI:10.1175/1520-0426(1995)012<0721:ROHWFS>2.0.CO;2
Lee H C, Lee Y H, Ha J C, et al. 2010. McGill algorithm for precipitation nowcasting by Lagrangian extrapolation (MAPLE) applied to the South Korean radar network. Part Ⅱ:Real-time verification for the summer season. Asia-Pacific J Atmos Sci, 46(3): 383-391.
Liu D C, Nocedal J. 1989. On the limited memory BFGS method for large scale optimization. Math Program, 45(1-3): 503-528. DOI:10.1007/BF01589116
Mandapaka P V, Germann U, Panziera L, et al. 2012. Can Lagrangian extrapolation of radar fields be used for precipitation nowcasting over complex alpine orography?. Wea Forecasting, 27(1): 28-49. DOI:10.1175/WAF-D-11-00050.1
Sun J Z, Crook N A. 1997. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part Ⅰ:Model development and simulated data experiments. J Atmos Sci, 54(12): 1642-1661.
Sun J Z, Crook N A. 1998. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part Ⅱ:Retrieval experiments of an observed Florida convective storm. J Atmos Sci, 55(5): 835-852.
Wang G L, Yang J, Wang D, et al. 2016. A quantitative comparison of precipitation forecasts between the storm-scale numerical weather prediction model and auto-nowcast system in Jiangsu, China. Atmos Res, 181: 1-11. DOI:10.1016/j.atmosres.2016.06.004
Wilson J W, Crook N A, Mueller C K, et al. 1998. Nowcasting thunderstorms:A status report. Bull Amer Meteor Soc, 79(10): 2079-2100. DOI:10.1175/1520-0477(1998)079<2079:NTASR>2.0.CO;2
Wilson J W, Feng Y R, Chen M, et al. 2010. Nowcasting challenges during the Beijing Olympics:Successes, failures, and implications for future nowcasting systems. Wea Forecasting, 25(6): 1691-1714. DOI:10.1175/2010WAF2222417.1