雷达反演的多尺度风场在临近预报中的应用研究

罗义，梁旭东，王刚，曹正，陈春元

(1.中国气象科学研究院，北京100081；2.中国民用航空中南地区空中交通管理局气象中心，广州510405)

引 言

多普勒天气雷达是探测中小尺度天气系统的有力工具，能够提供高时空分辨率的反射率、径向速度和径向速度谱宽的资料。在目前业务系统中，特别对于暴雨强对流天气的临近预报业务，外推预报依然是对流临近预报技术的重要基础(张家国等，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，1981；Lucas和Kanade，1981；Bowler et al.，2004)通过计算雷达回波的光流场来得到回波的运动矢量场，从而进行外推临近预报，其效果跟TREC方法有一定的类似性(曹春燕等，2015)。

与此同时，对于多普勒雷达径向速度资料的研究主要集中风场信息的反演方面。目前根据单多普勒雷达径向速度进行反演风场的主要方法有：VAD(Velocity Azimuth Display)方法(Lhermitte，1961)，VARD(Velocity Area Display)方法(Easterbrook，1975)，VVP(Velocity Volume Processing)方法(Waldteufel and Corbin，1979；Koscielny 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，1969；Brandes，1977；Kessinger et al.，1987；Shapiro and Mewes，1999)或者采用变分技术结合物理约束(Shapiro et al.，2003；Snyder and Zhang，2003；Tongand Xue，2005；李红莉和王叶红，2007；Xiaoet al.，2008；Shapiro et al.，2009)，进行反演二维或者三维风场的方法也得到大量研究。这些方法采用变分的方法，通过利用反射率守恒、运动方程甚至数值模式作为约束条件，以得到更加精确的风场。但是这些方法计算过于复杂且计算量大，目前还难以应用于临近预报。

总的来说，目前雷达临近外推预报主要利用反射率资料进行追踪算法来获得外推的运动矢量场，而如何利用径向速度资料获得反演风场作为外推运动矢量场的研究则较少。早期有研究对于雷暴的移动与环境平均风场的关系进行了统计分析，为利用环境风场外推预报提供了有力的可预报性依据。例如Wilson(1966)发现雷暴的回波移动都与对流的尺度相关，较小尺度的对流系统随着平均风移动，而尺度较大的对流系统移动比平均风速慢且偏向平均风向的右侧。Andersson和Ivarsson(1991)在概率临近预报模式中采用850 hPa风场作为外推风场，并考虑了利用多普勒雷达，通过VAD、UW等反演风场或者类似TREC反演风场作为引导风场进行外推预报的可能性。Turner等(2004)开发的临近预报系统(MAPLE,the Mc-Gill 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方法的基本原理是将反射率图像划分为一系列固定大小的区域，通过计算连续时次雷达回波不同区域的空间相关，搜索与目标区域具有最优相关系数的区域，从而得到回波的运动矢量场。其相关系数可以表示为

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

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)，在均匀风场条件假设下，在分析单元内，局地的平均反演风场公式可以表达为公式(2)、(3)中，表示雷达观测到的径向速度，θ表示方位角，r表示距离，u和v表示在分析单元中心点处的反演风场。

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

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

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

图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 vectorsof largescale(vector,unit:m·s-1)overlaid original radial velocity(shaded,unit:m·s-1)at the3 km level,(c)retrieved motion vectorsof small scale(vector,unit:m·s-1)overlaid small scaleradial velocity(shaded,unit:m·s-1)at the3 kmlevel,and(d)synthesizingmotion vectors(vector,unit:m·s-1)overlaid radar reflectivity at the3 kmlevel(shaded,unit:dBz)at 21:00 UTC 9 April,2016.

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

(4)将第(1)步得到的系统性运动矢量场u0和v0和第(3)步得到的雷暴内部的运动矢量场u1和v1进行合成，如图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内各层高度上的风场按照不同权重进行合成，即通常低层风场赋予更高权重，高层给予较小的权重。按照不同的高度给予不同的权重，其权重公式为

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

3 个例分析

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

参照第二部分MIVAP方法和MTREC方法的参数设置，对应得到两种方法的外推运动矢量场如图2。对比图2a、b可以看出：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 kmlevel 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 reflectivityover 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)的计算公式如下

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

图4 检验评分阈值为35 dBz时外推预报60 min评分(黑色线为POD,红色线为FAR,蓝色线为CSI,实线为MTREC方法,虚线为MIVAP方法)Fig.4 Evaluation scoreswith threshold value35 dBz:POD(black line),FAR(red line)and CSI(blue)within 60-minuteforecast 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 thecaseswith date,storm typeand mean evaluation using the MTRECand MIVAPmethods for 60-minutenowcasting.

15个例按照POD、FAR、CSI三个指标对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 scoreswith threshold value35 dBzbased 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 个例发生日期、雷暴类型和不同预报方法得到的外推运动平均值Table2 List of thecaseswith date,stormtypeand mean motion fromthe MTRECand MIVAPmethods.

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方法这类跟踪算法的局限性问题。