1. 引言
雷达回波的外推是对流天气临近预报的重要方法,业务上使用的0~2小时外推预报主要是基于雷达数据风暴识别追踪和自动外推预报技术[1]。自动外推技术主要包括了交叉相关法(TREC)和单体质心法。单体质心法是将对流云团视为三维单体进行识别和追踪,而交叉相关法是通过计算雷达回波在连续时次不同区域的最优相关系数,确定其过去的移动特征,通过一定的拟合关系实现回波形状和移动的外推。相较于单体质心法,交叉相关方法在计算方法上比较简单,但随着外推时间增加,交叉相关法计算的运动矢量场质量显著降低,对外推预报的质量影响较大。
光流法是计算机视觉领域中的重要方法之一,最早由Gibson [2]提出。光流法具有高弹性、高灵活的特点,可以模拟出比较理想的运动矢量场,能够较好的对目标物进行识别和跟踪。国内外学者利用光流法在雷达回波临近外推和降水临近预报上做了大量研究。曹春燕等[3]利用光流法和交叉相关法进行对比实验表明光流法的外推效果优于交叉相关法;张蕾[4]等利用金字塔分层技术改进了HS光流算法,并在考虑回波旋转效应的基础上外推回波获得了更好的稳定性和精度;柳士俊[5]等认为在全局平滑假设下,即便对流云团运动和形态发生较大的变化,光流算法也能够较好的描述这种变化,这是交叉相关法无法相比的;Georgy Ayzel [6]等利用多种光流算法进行降水临近预报试验,并由此研发了降水临近预报的rainymotion光流外推框架。可见利用光流法进行雷达回波的外推,相较于其它外推方法具有先天的优势。
本文拟采用Georgy Ayzel [6]等提出的2种稀疏光流模型和2种稠密光流模型对乌鲁木齐机场单部多普勒雷达的组合反射率进行外推实验,探讨4种光流模型在乌鲁木齐机场雷达回波外推的适用性,为建立基于单部雷达回波外推的临近预报系统做准备。
2. 方法介绍
光流是空间运动物体在观测成像平面上的像素运动的瞬时速度,是利用图像序列中像素点在时间阈上的变化以及相邻帧之间的相关性来找到上一帧与当前帧之间的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法[7]。光流法外推雷达回波数据就是利用至少2张临近的雷达回波图像,在引入附加约束条件的基础上,求解光流约束方程,获得雷达回波运动场,利用此运动场对回波的未来时刻进行外推,从而获得未来时刻雷达回波。本文采用Georgy Ayzel [6]等提出的利用光流法进行基于雷达的降水临近预报的4种方法,利用光流法对大尺度背景下强降水天气回波、迅速发展变化对流回波、强度及形态变化较缓慢的对流回波等3组天气过程进行实验和评估,探讨该方法在乌鲁木齐机场进行单部多普勒雷达回波外推中的适用性。
2.1. 计算方法
Georgy Ayzel [6]等在2019年提出了利用光流法进行基于雷达的降水临近预报的4种方法,分别是稀疏组的SparseSD及Sparse方法和稠密组的Dense及DenseRotation。
SparseSD模型仅使用2张临近的雷达回波图形来进行图形特征的识别、跟踪和外推。假定t表示即将进行雷达回波外推的最近一张图像,首先利用Shi-Tomasi角点检测器[8]识别t-1时刻雷达回波图像中的特征,使用局地Lucas-Kanade光流算法对这些特征在t时刻进行跟踪,在t时刻的雷达回波图像上识别出t-1时刻雷达回波图像上特征的位置,线性外推这些特征的运动,获取每个预报时效n的特征位置,使用最小二乘法,根据已经识别特征在t和t + n时刻的位置计算每个预报时效n的仿射变换矩阵,使用相应的仿射变换矩阵将t时刻的雷达回波图像每个像素位置变换为t + 1、t + 2∙∙∙t + n,并对剩余的不连续点进行线性插值,从而得到预报时效n内的连续雷达回波外推图像。而Sparse模型则采用的是预报时刻前最近的24张连续的雷达回波图像,使用Shi-Tomasi角点检测器进行t-23雷达回波图像特征识别,使用局地Lucas-Kanade光流算法跟踪雷达回波图像从t-22至t时刻已识别的特征,对于每个成功跟踪的图像特征,建立从t-23到t时刻坐标变化的线性回归模型,后续计算与SparseSD一致。
本文中稠密组使用的是Kroeger等在2016年提出的全局光流算法Dense Inverse Search (DIS) [9],这种算法能够根据两张连续的雷达回波图像来估计速度场。Georgy Ayzel [6]等在临近预报实验中认为,该方法比DeepFlow等其他全局光流算法相比有更高的准确性和计算效率。稠密组的Dense和DenseRotation在计算方案上是一致的,只是外推上使用的方法不同,Dense使用的是恒定向量平流方案,该方案不考虑回波的旋转效应,而DenseRotation使用的是考虑了旋转效应的半拉格朗日平流方案。
Dense和DenseRotation模型使用时刻t-1和t时刻的雷达回波图像,利用全局DIS光流算法计算连续速度场,然后使用向前或后向恒定向量方法(Dense)和半拉格朗日方法根据速度场外推n个时效的雷达回波图像像素,最后使用插值算法对外推的回波进行插值,最终得到预报时效n内连续的雷达回波外推图像。
本文进行雷达回波外推实验与Georgy Ayzel等进行降水临近预报有所不同,本实验直接对雷达回波图像进行外推,无需将雷达回波数据转化为降水率,只需通过dBZ矩阵与灰度矩阵的变换来实现外推的inputs和outputs,进行效果评估时也只对dBZ矩阵进行评估。外推实验输入为最近2小时的雷达回波图像,预测未来2小时逐6分钟的回波图像。
2.2. 雷达数据预处理
文中所使用的雷达数据为乌鲁木齐国际机场单部雷达多普勒雷达数据,为了简化数据处理,直接使用已经经过质量控制的组合反射率进行外推实验,这样可省去对雷达数据进行质量控制的步骤。本文实验的组合反射率数据时间分辨率为6 min,空间分辨率为1 km,范围为机场周边150公里半径,在回波外推前需要做以下几个方面的处理:
(1) 将组合反射率基数据中大于70 dBZ和小于等于5 dBZ的回波数据滤除,使得组合反射率的回波边界更加清晰,同时也能消除部分孤立的回波噪声。
(2) 由于光流法所需要的雷达回波图像为灰度亮度图像,因此需要将组合反射率dBZ矩阵变换为灰度矩阵。文中dBZ矩阵变化为灰度矩阵使用的是Python的numpy库,变换关系为:
import numpy as np
gray_data = np.clip(dbz_data, R_MIN, R_MAX)
R_gamma = 1
gray_data = 255 * ((gray_data - R_MIN)/(R_MAX-R_MIN))**(1/R_gamma)
gray_data = np.clip(gray_data, 0, 255)
上述关系中R_MIN和R_MAX为组合反射率回波反射率的最大和最小值常量,乌鲁木齐机场多普勒雷达的最小值为−32 dBZ,最大值为94.5 dBZ。
(3) 对由(2)中得到的灰度矩阵进行滤波来消除矩阵中的孤立噪声,常用的滤波方法包括均值滤波、中值滤波、高斯平滑等。本文实验采用中值滤波,该方法能够有效的消除雷达回波噪声,且对回波的测量值损失小,平滑效果好,边缘清晰[10]。通过中值滤波进一步消除灰度矩阵中的孤立噪声,中值滤波器使用OpenCV的medianBlur函数即可实现。
(4) 将(3)中已经处理好的灰度矩阵数据直接存储为组合反射率灰度图像以备后续计算使用。由于乌鲁木齐机场单部多普勒雷达数据量较小,实验时不对雷达图像的分辨率进行降低,在实际应用过程中,若分辨率过高,会出现图像纹理过于细密,不容易满足光流算法图像恒定的假设,这时需要对图像的分辨率进行降低以减少纹理的出现[7],可获得更好的外推效果。
2.3. 外推回波后处理
按2.1计算方法对乌鲁木齐机场多普勒雷达组合反射率未来2小时逐6分钟进行外推,可得到20个外推回波的灰度矩阵,对这20个灰度矩阵进行处理即可得到外推的组合反射率回波图像,处理方法是对2.2中第(2)步进行灰度至dBZ的反变换,反变换关系为dbz_data = grey_data/255*(R_MAX-R_MIN) + R_MIN,其中R_MIN和R_MAX与1.2中第(2)步的值一致,否则变换后的回波图像会出现错误,对dbz_data进行绘图即可得到彩色的组合反射率因子。同时将计算得到的回波外推灰度矩阵进行直接存储,为后续评估做准备。
2.4. 评估方法
(1) 对4组不同计算方法外推的雷达组合反射率回波在t + 6、t + 30、t + 60、t + 90、t + 120与机场实际观测到的组合反射率回波进行比较,从回波发展和形态变化上进行定性评估。
(2) 采用平均绝对误差(MAE)和预报技巧评分即命中率(POD)、空报率(FAR)及临界成功指数(CSI)等进行客观评估。在客观评估过程中,通过读取外推2小时的20张组合反射率灰度图形和对应雷达实际探测的组合反射率灰度图像,将灰度图像利用1.3中的反变换方法转换为dBZ矩阵,然后对每个格点上的dBZ值进行评估。本次实验评估时对dBZ矩阵中的NAN数据取0值,这样整个dBZ矩阵值域范围转化为0~70 dBZ,评估对此值域范围的数据进行。
平均绝对误差MEA计算方法如下:
(1)
上式中nowi和obsi为外推预报和实际观测的第i个像素dBZ值,n为dBZ矩阵像素总数。对于MAE而言,值越小,表明模型外推的误差越小。
对于POD、FAR、CSI评分采用最低阈值为10 dBZ进行计算,POD、FAR、CSI按以下公式进行计算:
(2)
(3)
(4)
其中N成功、N漏报、N空报分别为外推预报中成功、漏报、空报的格点数。
3. 临近外推实验
为了验证4种光流模型在乌鲁木齐机场多普勒雷达组合反射率回波外推中的适用性,本文挑选了3次天过程进行分析(文中涉及的时间均为世界时UTC),分别是2021年7月31日大尺度强降水天气过程、2022年6月28日快速发展的雷雨大风天气过程、2022年7月24日发展较缓慢的雷雨大风天气过程,对这3个过程的4种光流模型外推回波进行定性分析和客观评分。
3.1. 案例1:2021年7月31日大尺度强降水天气过程
2021年7月31日,受较强冷空气影响,新疆天山北坡出现一次大范围的大到暴雨,局地暴雨天气过程,乌鲁木齐机场从31日上午开始出现了长时间的降水,过程降水量达到15.1 mm。从卫星云图和雷达回波图(图略)上分析,导致此次强降水天气的主要原因是强降水回波的列车效应,即西南气流不断生成的降水回波持续经过该区域,回波东移速度缓慢,对该区域的影响时间较长。
图1为4种光流模型从7月31日03:00 UTC开始对组合反射率进行外推预报和机场雷达观测的对比。如图1所示,4种模型在外推T + 30内,其回波形态及强回波区的预报与实际探测基本一致,至T + 60后,Dense和DenseRotation模型在回波形态上略有变形,但强回波的位置与实际较为接近,30 dBZ以上的强回波区的预报较实际略偏小;而Sparse和SparseSD模型预报的回波形态出现顺时针偏转,导致其相对于实际回波在西段偏北而东段偏南,Sparse比SparseSD偏转得更加明显,两个模型对30 dBZ以上的强回波区的预报较实际偏大。从回波预报的演变来看,稠密光流算法Dense和DenseRotation的外推预报效果优于稀疏光流算法Sparse和SparseSD,而SpareSD优于Sparse。
图2为本次过程针对4种光流模型外推预报中大于等于10 dBZ的回波计算得到的均方根误差、命中率POD、空报率FAR及临界成功指数CSI。从图2(a)平均绝对误差来看,4种模型的平均绝对误差随着预报时间均不断增大,而稠密光流算法Dense和DenseRotation的预报效果明显优于稀疏光流算法Sparse和SparseSD,Dense和DenseRotation算法变化相当,SparseSD优于Sparse。从CSI、POD演变来看,随着外推时间增长,4个模型的两个评分均呈下降趋势,FAR均呈上升趋势,但稠密光流算法均优于稀疏光流算法,尤其在外推30 min后,考虑旋转效应的DenseRotation模型的CSI评分最高,空报率最低。
综合来看,就此次个例而言,由于在大尺度环流约束下,实际探测回波强度和移动演变均比较缓慢,4种外推模型均在2 h内的对回波移动和相态的外推预报均有较好的效果,尤其是稠密光流算法其准确率在2 h内均达到了70%以上,空报率均最高控制在20%左右,命中率达到85%以上,结合平均绝对误差,至少85%的预报回波其误差不超过4 dBZ。可见,稠密光流对此类天气系统回波进行2 h的外推预报是成功的,尤其考虑了旋转效应的DenseRotation模型在此类天气下具有更好的外推预报效果。
Figure 1. (a) Mean Absolute Error, (b) CSI Score, (c) POD Score, and (d) FAR Score of 2-hour extrapolation forecasts from 4 optical flow models starting at 03:00 UTC on July 31, 2021
图1. 2021年7月31日03:00 UTC开始用4种光流方法对组合反射率回波进行6 min、30 min、60 min、90 min、120 min的外推预报与实际观测回波的对比
Figure 2. (a) Mean Absolute Error, (b) CSI Score, (c) POD Score, and (d) FAR Score of 2-hour extrapolation forecasts from 4 optical flow models starting at 03:00 UTC on July 31, 2021
图2. 2021年7月31日03:00 UTC开始4种光流模型2小时外推预报的 (a) 平均绝对误差、(b) CSI评分、(c) POD评分及 (d) FAR评分
3.2. 案例2:2022年6月28日快速发展的雷雨大风天气过程
2022年6月28日,受短波槽影响,在乌鲁木齐机场及其西北方向出现了雷雨大风天气过程。乌鲁木齐机场在28日12:40~14:00 UTC出现雷阵雨,偏西风平均风12~16 m/s,阵风达18 m/s,短时伴有700~900 m沙尘暴,降水量1.1 mm。此次天气在组合反射率上的特点是回波发展变化迅速,对流单体不断生成和消散,移动变化复杂。
图3为4种光流模型从28日11:49 UTC开始对组合反射率进行临近外推预报与实际探测组合反射率回波的对比。如图3所示,从观测回波演变来看,此次过程回波演变复杂迅速,具有强回波不断生成、合并、迅速移动的特点。从4组光流模型外推的回波演变上看,所有模型都无法预报出对流回波生成和合并的过程。稠密光流算法Dense和DenseRotation能够推算出整个回波范围变大的过程,但其变化范围较实际仍有较大的误差,尤其是从T + 60开始,强回波相互合并后快速向东南方向发展的过程,稠密光流算法也无法进行推算,在回波形态和强度范围上均产生了较大的误差。而稀疏光流算法的Sparse及SparseSD对此次回波的外推,除了T + 30前回波变化较为缓慢时有一定效果外,之后的回波以类似刚体移动的方式向东北方向移动,且回波向西北–东南向形变。
就此次案例而言,稠密光流算法表现出比稀疏光流算法更优的结果,但其在1 h后的外推预报就已经与实际有了较大的误差。可见对于此类快速发展的强对流系统,稠密光流算法在1 h后就已经失去了外推的价值。
Figure 3. Comparison between extrapolation forecasts and actual observed echoes of composite reflectivity echoes using 4 optical flow methods for 6min, 30min, 60min, 90min, and 120min starting at 11:49 UTC on June 28, 2022
图3. 2022年6月28日11:49 UTC开始用4种光流方法对组合反射率回波进行6 min、30 min、60 min、90 min、120 min的外推预报与实际观测回波的对比
图4为本次过程针对4种光流模型外推预报中大于等于10 dBZ的回波计算得到的均方根误差、命中率POD、空报率FAR及临界成功指数CSI。从平均绝对误差可见,4个模型的平均绝对误差随着外推时间呈明显的增长趋势,但稠密光流算法的平均绝对误差在100 min前比稀疏光流算法低1~2 dBZ,在65 min前稠密光流算法平均绝对误差维持在3 dBZ以下。从CSI评分和POD评分曲线上看,随着外推时间两者在下降趋势上基本一致,稀疏光流算法的评分呈快速下降的趋势,外推40 min后CSI评分低于了0.6,POD评分Sparse在40 min后低于0.6,SparseSD在60 min后低于0.6;而稠密光流算法的CSI评分在80 min后低于0.6,POD评分则在85~90 min后才低于0.6,在至少45~50 min后不考虑旋转效应的Dense模型略优于考虑旋转效应的DenseRotation模型。对于这次过程,在60 min前的FAR评估有一定的意义,60 min后实际探测回波呈发散状,外推预报即便在形态和强度上外推不准确,但回波仍然落在了实际探测回波中,所以其FAR可能保持在较低的水平。从图4(d)的FAR评分中,50分钟前Dense、DenseRotation、SparseSD的FAR都较低,Sparse相对较高一些。
Figure 4. (a) Mean Absolute Error, (b) CSI Score, (c) POD Score, and (d) FAR Score of 2-hour extrapolation forecasts from 4 optical flow models starting at 11:49 UTC on June 28, 2022
图4. 2022年6月28日11:49 UTC开始4种光流模型2小时外推预报的 (a) 平均绝对误差、(b) CSI评分、(c) POD评分及 (d) FAR评分
综合上述定性、定量评估上看,稠密光流算法在此类过程中的外推效果明显优于稀疏光流算法,且稠密光流算法的Dense模型最优。就本次过程而言,光流算法无法对回波的生成、合并等过程进行预报,1小时后外推模型已经无法描述真实回波的形态,但从客观评分上看,稠密光流算法在60~90 min的预报还有一定的参考性,但回波形态已经与实际有了较大的差异。
3.3. 案例3:2022年7月24日发展较缓慢的雷雨大风天气过程
2022年7月24日,受短波槽影响,乌鲁木齐机场及以西区域受雷雨天气影响,乌鲁木齐机场在24日17:07观测到闪电,17:49~18:30 UTC为雷雨,并伴有23 m/s的西北大风。此次天气过程在组合反射率回波上的特点是回波在前期没有快速触发、合并的过程,回波移动较为规律,但在后期随着雷暴云团的减弱消散,回波也出现明显的破碎和减弱。
图5为24日16:49 UTC开始用4种光流模型对组合反射率回波进行外推预报与实际观测回波的对比分析。如图所示,4光流模型在T + 60前的外推预报均能较好的描述实际回波形态的演变,T + 60后随着对流天气减弱,探测回波破碎消亡,光流法的外推预报已经失败。稠密光流算法的Dense和DenseRotation模型在回波形态和位置上更接近于实际探测,但随着外推时间,这两种模型的回波出现了一定的形变,边缘出现模糊,尤其在T + 90后,Dense模型的形变和模糊程度均比DenseRotation大。相比于稠密光流算法,稀疏光流算法的回波形变和边缘模糊相对较小,但Sparse和SparseSD的外推回波均较实际偏北,移动速度偏慢。
可见,在此类过程中,稠密光流算法比稀疏光流算法在回波位置和速度上更接近与实际,同时在前60min不考虑旋转效应的Dense模型在回波形态及强回波范围上比考虑旋转效应的DenseRotation模型更接近于实况。
Figure 5. Comparison between extrapolation forecasts and actual observed echoes of composite reflectivity echoes using 4 optical flow methods for 6 min, 30 min, 60 min, 90 min, and 120 min starting at 16:49 UTC on July 24, 2022
图5. 2022年7月24日16:49 UTC开始用4种光流方法对组合反射率回波进行6 min、30 min、60 min、90 min、120 min的外推预报与实际观测回波的对比
图6为本次过程针对4种光流模型外推预报中大于等于10 dBZ的回波计算得到的均方根误差、命中率POD、空报率FAR及临界成功指数CSI。从图6(a)中的平均绝对误差可见,稠密光流算法比稀疏光流算法具有更大的均方根误差。稠密光流算法中,在前60 min两种模型的均方根误差基本一致,在4dBZ以下,60 min后随着外推过程中的形变和边缘模糊,Dense模型的均方根误差开始比DenseRotation模型的误差大。但从CSI评分和POD评分上看,如图6(b)和图6(c)所示,Dense模型在两个评分上表现得更优,尤其在命中率POD上比其它3个模型更加优秀。而在空报率上,DenseRotation模型由于在形变和边缘模糊上较Dense模型更小,因此其空报率FAR较Dense模型更小一些。稀疏光流的两个模型在CSI和POD评分上比稠密光流的评分更低,空报率FAR更高。
Figure 6. (a) Mean Absolute Error, (b) CSI Score, (c) POD Score, and (d) FAR Score of 2-hour extrapolation forecasts from 4 optical flow models starting at 16:49 UTC on July 24, 2022
图6. 2022年7月24日16:49 UTC开始4种光流模型2小时外推预报的 (a) 平均绝对误差、(b) CSI评分、(c) POD评分及 (d) FAR评分
综合上述主观、客观分析来看,稠密光流算法的外推效果仍优于稀疏光流算法,且在60 min前Dense模型的外推回波形态和强度较DenseRotation更接近于实际,而在CSI评分和POD评分上,Dense模型均比DenseRotation模型更优。
4. 结论与讨论
本文采用Georgy Ayzel [6]等提出的4种光流模型,分别2021年7月31日强降水过程、2022年6月28日快速变化的雷雨大风过程及2022年7月24日变化较缓慢的雷雨大风过程利用乌鲁木齐国际机场单部多普勒雷达的组合反射率进行未来2小时逐6分钟的临近外推实验,对预报结果进行主观、客观评估,结论如下:
(1) 利用光流法对雷达组合反射率图像进行临近外推是可行的,对单部多普勒雷达而言,稠密光流算法较稀疏光流算法在雷达回波的临近外推上更具优势。
(2) 4种光流模型对大尺度系统约束下演变缓慢的回波外推效果均比中小尺度对流系统好,对稠密光流算法而言,外推时间越短,准确性越高,中小尺度对流系统在60 min后外推预报参考价值较低,而对大尺度系统约束下缓慢演变的回波其外推预报至2小时仍具有较高的准确性,具有较高的参考价值。
(3) 光流法无法对回波的生消、强弱变化进行很好的外推预报,因此当回波在机场、或关注区域附近生成时,光流法外推可能是失败的,对此类系统的外推预报还需要如数值预报等手段进行辅助。
(4) 对于大尺度条件约束下的降水回波,考虑旋转效应的DenseRotation模型比其它模型有更好的外推效果,而对于类似中小尺度系统的变化较快的回波,不考虑旋转效应的Dense模型具有更好的外推效果。
(5) 由于本次外推实验采用的是乌鲁木齐机场单部多普勒雷达的组合反射率图像,回波范围小、变化也相对较快,尤其在雷暴等强对流系统下,光流算法的约束条件不容易满足,随着外推时间产生较大的误差,若以大范围的雷达拼图作为基础数据,降低数据的分辨率,则可能产生较好的外推效果。但从评估结果来看,利用稠密光流算法的Dense和DenseRotation模型建立乌鲁木齐机场单部雷达0~1小时客观外推预报系统是可行的,能够为机场预报、流量管理等航空运行战术决策提供有价值的数据支撑。
基金项目
民航新疆空中交通管理局科技项目《新疆辖区内强回波覆盖率监测平台》。