基于拉格朗日力学的雷达图像追踪预测方法

王兴+罗永月+苗春生+张皓

摘 要：提出一种多普勒天气雷达回波图像追踪预测方法。该方法以拉格朗日力学作为动力学基础，利用雷达探测到的基本反射率因子作为输入场，结合分层抽稀技术捕捉回波在空间和时间上的位移矢量，并由此预测雷达回波的位置和强度，进而为短临天气预报业务的实施提供必要支撑。实验表明，该方法能够较准确地预测局地天气演变过程，1小时预测结果与实况的相关性达到70%以上，相关性的个体差异在20%以下。同时，表明该方法对于辅助预报未来1小时内的强对流事件是行之有效的。

关键词：拉格朗日力学； 图像预测； 图像追踪； 天气雷达； 临近预报

DOIDOI：10.11907/rjdk.162240

中图分类号：TP317.4

文献标识码：A文章编号：1672-7800（2016）012-0001-04

0 引言

多普勒天气雷达（以下简称雷达）是现代气象业务研究及应用中不可或缺的重要工具，基于雷达回波及图像的追踪与外推是短时临近预报的关键性技术，也是长期以来的难点问题。准确而及时的雷达外推预报，可以为局地龙卷风、雷暴、短时强降水等极端灾害性天气提供预警，以便做好防御措施，最大限度保障人民及生命财产安全。

近半个世纪以来，很多学者在基于雷达资料的预报领域作出了大量卓有成效的贡献，并逐步形成了一系列较为通行的方法。例如，TREC算法通过逐区域寻求相邻时刻雷达反射率的最大相关，跟踪整个回波区域的移动，并且假设回波具有一致的移动方向[1]。CTREC算法则利用交叉相关分析，跟踪反射率因子大于一定阈值区域的移动，进而推算回波的发展[2]。TITAN是由美国国家大气研究中心（NCAR） 研发的一套风暴识别、跟踪、分析和预报系统，其利用雷达一次完整体扫所构成的三维结构数据对强回波中心进行识别追踪[3-4]。SCIT算法则更加侧重对雷暴单体的有效追踪和预测[5-7]。近年来，又有一些学者试从图形图像学中的光流技术入手，通过分析雷达回波时序图像中的光流场特征进行回波强度和位置的外推预测[8-11]。此外，还有基于神经网络、模式匹配等多种方法[12-14]。

总体来说，这些方法的共同之处是需要分析相邻时刻雷达图像（或基数据）中区域的相似性，而该相似性的度量，如最大相关法，虽然可以计算出最匹配的位置，但匹配和外推结果往往表现出发散性或多个最优解。并且，基于窗口平移的模板匹配算法无法适应区域图像的旋转和变形等情况。

为了克服相关性度量算法所遇到的问题，本文考虑对位移场的分析加以分层，也就是每个特征运动被认为是确定在相对粗糙的空间分辨率下的平稳变化趋势的总和，然后衍生为更高的空间分辨率下的小幅度局部修正，并对此过程进行多次迭代。对于位移场的分析，将充分考虑雷达回波所指示风暴的动力学特征，运用拉格朗日力学相关理论构建预测模型。该方法的研究意义在于进一步提高中小尺度、强对流天气事件的预测能力，且相对于光流等一些大运算量算法，本方法运算规模更小，进而能够更好地满足短临预报业务高时效性的要求。

1 理论与方法基础

1.1 雷达回波外推预测

大量研究表明，合理的外推预测算法可以为降水、雷暴、冰雹等对流天气的预报提供重要支撑[15-17]。基于外推预报的一般性描述为：

其中，pt（x，y）表示任一位置的回波强度，U和V分别表示回波在水平和垂直方向上的偏移量，由U和V共同组成回波移动的速度矢量。g表示一个函数，用来计算单位时间间隔后回波强度的变化。根据式（1），Δpt（x，y）反映了任一点（x， y） 在t 时刻回波强度的变化情况，U和V反映了回波移动的方向和速度。此外，函数g代表一个拉格朗日动力学过程，在此过程中雷达回波的强度是由其沿回波路径移动时在拉格朗日坐标系统中的历史变化推导出的，也就是用当前回波演变的趋势预测回波未来的位置和强度。究其趋势预测的方法，多年来诸多学者作出很多研究，本文主要从拉格朗日力学角度进行分析，提出一种雷达回波图像追踪预测的方法。

1.2 拉格朗日力学

拉格朗日力学是由Joseph Lagrange[18-19]最早提出的一种力学分析方法。由于该方法引用了广义坐标的概念，使得对力学相关问题的研究更具普适性。

在雷达回波图像预测研究中，如何准确得到回波运动矢量是预测需要解决的关键问题。在不考虑天气系统的非线性变化时，拉格朗日力学模型能够满足构建回波发展演变过程的算法要求，式（1） 可以改写为：

有研究表明，对整个回波图像采用统一的U和V所构成的位移矢量，可适用于对大尺度天气系统的预测分析，如对卫星图像的外推预测和云导风的分析[20]。但对于局地强对流天气系统，预测结果往往与实际偏差较大。因此，这也是本文将重点阐述解决的问题。

2 雷达图像追踪预测

2.1 基于拉格朗日力学的追踪算法

在上述理论基础上，根据大气运动演变发展的规律及其在雷达回波图像上的表征特点，构建基于拉格朗日力学的追踪算法模型，如式（3） 所示。

该模型假定所预测的回波图像是当前和过去若干个回波图像以固定时间间隔而变化的函数。f2表征一个用于估测单位时间间隔前后回波图像各相应网格点回波强度值变化率的函数，即回波的演变趋势。在不断生消、发展的对流系统中，ΔP可以为正，也可以为负数。式（4）～式（6） 进一步表明了函数f1和f2，即回波移动矢量及回波强度的计算方法。

2.2 中心极值滤波

为了降低雷达杂波对位移矢量计算的不良影响，本节提出采用一种滤波器对雷达基数据进行滤波处理。其基本思想是：逐网格分析回波强度特征，如果某格点的值大于周边最相邻一圈（共8个网格）的最大值，或者该值小于周边最相邻一圈的最小值，则将当前网格点回波强度值用上述8个网格的最大值或最小值替代。

如图1所示，位于当前中心点的数值39大于其最邻近一圈8个网格的最大值。因此，使用数值23替换当前网格的39。

从图像上看，该滤波方法可显著降低图像中的椒盐噪声，从实际效果上看，该方法可以有效过滤单点的杂波奇异值，较传统的均值滤波和中值滤波更好地保留了回波细节[21-22]，特别是回波中梯度变化较大的边缘区域。

2.3 分层外推预测算法

由于天气系统复杂多变，特别是尺度较小的局地强对流，其生命周期短的只有几分钟到几十分钟，由于其空间尺度小，生消速度快，因此，包括基于拉格朗日力学在内的各种线性关系外推算法，其预测准确性都存在一定的局限性。为了改善这一问题，本节提出采用分层的外推预测方法。该思想最早由Bellerby等 [20]提出，并研究应用于卫星图像的云顶平流场分析中。

该算法的关键流程为：先将当前雷达回波图像逐级抽稀，降低图像的分辨率，以模糊回波细节，由此粗略估算出回波主体的移动趋势；然后再反向逐级提高图像分辨率，在较粗的移动趋势基础上，细化和订正位移矢量的细节。从而产生一个在空间上连续和平滑的且不受模板边界不连续性影响的矢量场。计算方法如式（7）所示：

对于每一级抽稀计算，都是将当前一级各网格点的回波强度值经由公式（7）计算，并往复迭代。其中P表示某一点的回波强度，L和L-1代表抽稀的层级，在本文下述实验中，采用的最高层级为4。

在计算两个相邻时刻图像中回波的位移时，可以借鉴交叉相关法，计算方法如式（8），在每个选定的匹配窗口遍历出最大相关矩阵的位置，从而输出位移矢量。

式（8） 中，P和P 分别表示相邻两个时刻（如t-Δt与t）的回波，（x， y） 表示图像中的某一点，X和Y表示匹配窗口的大小。再将两幅回波图像之间的网格还原或内插到其先前空间分辨率的两倍，重复上述匹配。该迭代过程还考虑到了由非矩形网格代表的局部扭曲，结合这些局部扭曲，使外推预测算法能够适应旋转、扩展、缩小等回波图像形态上的变化。如此插值和匹配计算，迭代直到网格分辨率达到原始雷达图像分辨率。

3 实验与结果分析

3.1 实验数据说明

为检验所述方法的预测效果，本实验数据使用2016年6月南京地区多普勒天气雷达的基数据文件。该雷达使用VCP-21体扫模式，探测周期为6分钟。实验选用1.5°和2.4°仰角的基本折射率数值。为方便计算，实验前将原始数据由极坐标系统转换为平面直接坐标系统，数据的图像分辨率为920×920。为减少样本数量，从全部7199个基数据文件中筛选出以230库长为半径，其覆盖区域内具有大面积强回波的数据文件，共计880个。

3.2 实验结果分析

为检验雷达回波图像预测的准确性，使用与预测同一时刻的雷达实际探测数据作比对分析，计算过程采用交叉相关检验法。

首先以自然日为单位，统计逐日样本数据中每份预测结果与实况交叉检验的相关系数的平均值，如图2所示。

图2中3种图案标记分别表示预测6分钟、30分钟和60分钟的检验结果，每个值代表当日所有样本检验结果的平均值。横坐标为2016年6月的逐个日期，纵坐标为相关系数，其中横坐标4、5、9、10等日期没有标记图案，原因是这些日期的当日为晴天或少云，体现在雷达上没有强的大面积回波，因此没有列入样本数据进行分析。从图2中还可以看出，本方法预测未来6分钟的结果与实况相比，相关系数超过87%，平均达到93%以上，而随着预测时效的延长，预测准确率逐步下降，在未来60分钟的预测中，全月平均相关系数为70%左右。

进一步统计分析每次预测准确率的稳定性。以6月19日全天样本数据为例，统计每批样本所预测6、12、18至60分钟结果分别与实况交叉相关检验的情况，如图3所示。

图3中，每个柱状条的顶端和底端分别表示检验的相关系数的最大值和最小值，柱状条中间的黑色方形表示相关系数的均值。可以看出，随着预测时效的增长，其预测准确率的个体差异也随之增大。在前6分钟的预测中，该差异约为3%，30分钟时约为8%，而到预测60分钟时，差异进一步增大到20%。结果与强对流天气系统具有生命史短、突发性强，水气生消发展变化快的特点是相一致的。

4 结语

由于天气系统复杂多变，特别是对于中小尺度的对流系统，其生消、发展时间短、变化快，如何进行准确、有效的预报是提升当今短时临近预报的关键环节之一。考虑到大气中水气等物质的移动变化应遵循一般力学规律，而拉格朗日力学正是表征和计算动力学问题的普适性方法，因此，本文的预测动力模型建立在拉格朗日力学关系基础之上。又由于天气系统的变化表现在雷达图像上，其回波图形具有相当的不确定性，因此，本文提出采用分层的位移场分析方法，先假定位移矢量是在相对粗糙的空间分辨率下的平稳变化趋势的总和，然后在更高空间分辨率下作小幅度局部修正，并如此进行多次迭代。为了减少低仰角杂波对实验结果的影响，提出采用中心极值滤波对实验数据进行处理。结合上述理论构建起基于拉格朗日力学的追踪预测模型及算法流程，以雷达基本反射率因子作为输入场，追踪和预测回波在空间和时间上的位移矢量，并由此预测未来一段时间雷达回波的位置和强度。

实验部分采用1个月样本数据对本算法模型进行检验，通过预测结果与同时刻实况的比对分析，得出两者的相关性和个例稳定性等评价指标。结果表明，该方法能够较好地预测局地天气系统的演变过程，在未来30分钟的预测中准确率平均超过80%，且对于辅助预报未来60分钟内的局地龙卷风、强降水、雷暴等灾害性事件具有实践应用的价值。

参考文献：

[1] 刘红艳，魏鸣.多普勒雷达风场资料在临近预报中的应用[J].大气科学学报，2015（4）：483-491.

[2] 郑永光，林隐静，朱文剑，等.强对流天气综合监测业务系统建设[J].气象，2013（2）：234-240.

[3] 周康辉，郑永光，蓝渝.基于闪电数据的雷暴识别、追踪与外推方法[J].应用气象学报，2016，（2）：173-181.

[4] DIXON M， WIENER G.TITAN：thunderstorm identification，tracking，analysis，and nowcasting—a radar-based methodology[J].Journal of Atmospheric & Oceanic Technology， 1993， 10（6）：785-797.

[5] SHAH S， NOTARPIETRO R， BRANCA M.Storm identification，tracking and forecasting using high-resolution images of short-range X-band radar[J].Atmosphere， 2015， 6（5）：579-606.

[6] 庄旭东，胡胜，陈荣，等.“雨燕”中风暴算法与新一代雷达SCIT产品的对比分析[J].热带气象学报，2011（3）：299-306.

[7] DUAN Y， XU Y， ZHI S.Application analysis of the hail suppression operation based on the improved SCIT Algorithm[J].Meteorology & Disaster Reduction Research， 2014（15）：23-29.

[8] 曹春燕，陈元昭，刘东华，等.光流法及其在临近预报中的应用[J].气象学报，2015（3）：471-480.

[9] GARCIA F， CERRI P， BROGGI A， et al.Data fusion for overtaking vehicle detection based on radar and optical flow[J].2012， 7（2272）：494-499.

[10] 王兴，王新，苗春生，等.基于GPU加速的雷暴追踪外推方法研究[J].南京师范大学学报：工程技术版，2015（1）：35-42.

[11] STAINVAS OLSHANSKY I， BILIK I， BIALER O.Doppler-Based Segmentation and Optical Flow in Radar Images： US20160084953[P].2016.

[12] 盛仲飙.BP神经网络在数据预测中的应用[J].软件导刊，2016（1）：147-148.

[13] 王利卿，黄松杰.基于多尺度卷积神经网络的图像检索算法[J].软件导刊，2016（2）：38-40.

[14] WANG X， GU Y H， MIAO C S， et al.Parallelization and performance optimization of radar extrapolation algorithm with OpenCL[J].Journal of Internet Technology， 2016（17）：323-330.

[15] 王丹.雷达外推预报与暴雨数值模式融合预报降水方法研究[D].北京：中国气象科学研究院，2013.

[16] 张蕾.多普勒雷达回波演变的动力学分析及临近预报算法改进[D].南京：南京信息工程大学，2015.

[17] FOX N I， WEBB R， BALLY J， et al.The impact of advanced nowcasting systems on severe weather warning during the sydney 2000 forecast demonstration project：3 November 2000[J].Weather & Forecasting， 2004， 19（1）：97-114.

[18] 李艳艳.相似空间中不变的欧拉-拉格朗日方程[J].河南大学学报：自然科学版，2014，03：273-276.

[19] ERICKSEN R E， GUITERAS J J， LARRIVEE J A， et al.A parachute recovery system dynamic analysis [J].Journal of Spacecraft & Rockets， 1967， 4（3）：321-326.

[20] BELLERBY T J.High-resolution 2-D cloud-top advection from geostationary satellite imagery[J].IEEE Transactions on Geoscience & Remote Sensing， 2006， 44（12）：3639-3648.

[21] 莫晓齐，何爱.基于自适应开关中值滤波算法的工程图像处理[J].软件导刊，2014（3）：55-59.

[22] 李佐勇，汤可宗，胡锦美，等.椒盐图像的方向加权均值滤波算法[J].中国图象图形学报，2013（11）：1407-1415.

（责任编辑：陈福时）