CN104977584B - 一种对流天气的临近预报方法及系统 - Google Patents

Abstract

本发明适用于天气预报预测技术领域，提供了一种对流天气的临近预报方法及系统，包括：获取探测区域内的雷达拼图；基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场；依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。在本发明中，采用了使用局部约束的光流法计算得到雷达回波的光流场，以此来代替交叉相关法中的运动矢量场，从而能够准确地完成对雷达回波的识别、追踪和运动估计，提升了对流天气临近预报系统的预报性能。

Description

一种对流天气的临近预报方法及系统

技术领域

本发明属于天气预报预测技术领域，尤其涉及一种对流天气的临近预报方法及系统。

背景技术

灾害性天气预报业务是预报预测体系发展的重要组成部分，尤其是对于地处沿海，经济高度发达，人口密集的地区，常受锋面低槽、热带气旋、季风云团等天气系统的影响，暴雨、雷暴、雷暴大风、冰雹、台风等灾害性天气常有发生，为灾害性天气多发区。根据权威部门统计，对城市构成影响的各类自然灾害中，超过80％的灾害损失是由气象灾害造成的，气象灾害造成的经济损失约占GDP总量的1％～3％。基于上述原因，气象防灾减灾对灾害性天气和气象灾害的监测、预报和预警能力提出了越来越高的要求，对重点时段、重点区域的对流天气临近预报的需求更为迫切。

随着科学技术的快速发展，多种遥感遥测资料的有效应用、数值天气预报技术及高性能计算机能力不断提高，现代天气业务进入了新的快速发展时期，目前，业务上主要采用基于雷达数据的雷暴识别追踪技术和自动外推预报技术来进行0～2小时的对流天气临近预报。

自动外推预报技术主要包括两类：单体质心法和交叉相关法：单体质心法是将雷暴视为三维单体进行识别、分析、追踪，对雷暴进行拟合外推来做临近预报；交叉相关法利用求雷达回波最优空间相关的方法，在二维区域内通过计算连续时次雷达回波在不同区域的最优空间相关系数，来确定雷达回波过去的移动矢量特征，建立不同时次雷达回波的最佳拟合关系，从而追踪到一定区域内雷达回波过去的移动特征，然后通过这些移动特征来外推雷达回波未来的位置和形状。与单体质心法只能用于对流降水系统相比，交叉相关法既可以跟踪对流降水系统，也可以追踪层状云降水，因此，交叉相关法在气象业务部门得到广泛使用。然而，在交叉相关法的使用过程中，发现对于局地生成的降水回波，以及强度和形状随时间变化很快的降水回波来说，交叉相关法给出的运动矢量场的质量降低，导致跟踪失败的情况显著增加。

发明内容

本发明实施例提供一种对流天气的临近预报方法及系统，旨在解决目前的对流天气临近预报中，交叉相关法给出的运动矢量场的质量降低，导致跟踪失败的情况显著增加的问题。

本发明实施例是这样实现的，一种对流天气的临近预报方法，其特征在于，包括：

获取探测区域内的雷达拼图；

基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场，对于所述光流场的光流(u，v)，所述u，v为所述雷达拼图其中一点的光流在x，y方向上的分量，该点在时刻t的灰度值为I(x，y，t)，相应的光流约束方程为：

I_xu+I_yv+I_t＝0，

其中：同时，引入进一步的约束条件来求解所述光流约束方程：

依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

本发明实施例的另一目的在于提供一种对流天气的临近预报系统，所述系统为基于地理信息系统GIS的可视化平台，包括服务器和客户端，所述客户端用于显示所述服务器的临近预报结果，所述服务器包括：

获取单元，用于获取探测区域内的雷达拼图；

运算单元，用于基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场，对于所述光流场的光流(u，v)，所述u，v为所述雷达拼图其中一点的光流在x，y方向上的分量，该点在时刻t的灰度值为I(x，y，t)，相应的光流约束方程为：

I_xu+I_yv+I_t＝0，

其中：同时，引入进一步的约束条件来求解所述光流约束方程：

预报单元，用于依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

在本发明实施例中，采用了使用局部约束的光流法计算得到雷达回波的光流场，以此来代替交叉相关法中的运动矢量场，用于进行对流天气的临近预报，从而能够准确地完成对雷达回波的识别、追踪和运动估计，提升了对流天气临近预报系统的预报性能。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的对流天气的临近预报方法的实现流程图；

图2是本发明实施例提供的对流天气的临近预报系统的结构框图。

具体实施方式

以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透切理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。

在本发明实施例中，采用了使用局部约束的光流法(Optical Flow)计算得到雷达回波的光流场，以此来代替交叉相关法中的运动矢量场，用于进行对流天气的临近预报，从而能够准确地完成对雷达回波的识别、追踪和运动估计，提升了对流天气临近预报系统的预报性能。

为了说明本发明所述的技术方案，下面通过具体实施例来进行说明。

图1示出了本发明实施例提供的对流天气的临近预报方法的实现流程，详述如下：

在S101中，获取探测区域内的雷达拼图。

例如，在本文实验中所使用的雷达资料为新一代S波段多普勒天气雷达(CINRAD/SA)扫描的强度场(反射率因子)拼图资料。在VCP21探测模式下,每隔6分钟对9个仰角进行扫描，雷达最低仰角为0.5°，强度场的最大探测半径为460公里，同时，采用时间同步软件进行时间同步，确保所有雷达能同步扫描。为计算方便，把极坐标格式的雷达资料用ARPS插值法插值到三维直角坐标系中。同时，由于雷达拼图资料经常存在杂波和亮带等虚假回波，而这些虚假回波经常出现在低层，因此，可以取2.5km高度的雷达CAPPI拼图数据作为雷达回波场，此层雷达回波能代表对流水平分布特征，比较具有代表性，并且在很大程度上避免了虚假回波和地物杂波的影响。

在S102中，基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场。

光流法的基本原理如下：以图像亮度变化作为识别对象，运动目标和观测器之间的相对运动在图像序列中产生瞬时位移，从而能够体现图像亮度模式的表观运动，图像中所有像素点的亮度光流就构成了图像的光流场，而光流法的核心就是从连续的图像序列中计算光流场。

因此，简单说来，光流场就是类刚体物体的速度矢量场。在天气预报领域，速度矢量一般用u，v两个分量来表达。假设笛卡尔坐标平面上有一点(x，z)，它代表三维空间中某一点(x，y，z)在图像平面上的投影，该点在时刻t的灰度值(光亮度值)为I(x，y，t)。假定该点在(t+Δt)时运动到(x+Δx，y+Δy)，在一定时间间隔Δt内灰度值保持不变，即：

I(x+uΔt，y+vΔt，t+Δt)＝I(x，y，t)， (1)

(1)式中，u，v分别是该点的光流在x，y方向上的分量。假设亮度I(x，y)随时间t平滑变化，可以将(1)式按泰勒公式展开，得到：

其中，e为包括Δx，Δy，Δz的二次以上无穷小项，(2)式消去I(x，y，t)，用Δt除等式两边，并取Δt→0的极限后，可求得：

(3)式实际上是的展开式，可以简写为：

I_xu+I_yv+I_t＝0， (4)

其中：(4)式就是光流约束方程，I代表的是像素点(x，y)在时刻t的灰度值，且每个像素处的偏导数都是已知的。

令 为图像点灰度的空间梯度，I_t为灰度随时间的变化率，(u，v)称为光流，则图像中的所有点的光流就构成了光流场。

从光流(u，v)的定义可以看到，光流场与交叉相关法所求得的运动矢量场在物理意义上是等价的，如果把传感器从摄像机换成多普勒雷达，被探测的目标从一般的运动物体换成雷达回波，则I(x，y，t)就是雷达图像上像素点(x，y)在时刻t的反射率因子值。因此,可以用光流法替代交叉相关法得到运动矢量场。

由于光流有两个未知变量u和v，而公式(4)只有一个光流约束方程，如果不对光流特性做任何假设，则无法同时求出u和v，因此需要引入进一步的约束条件，才能完全确定光流(u，v)。为了求得u和v，我们认定在图像平面内足够小的区域内，且在足够短的时间间隔内，两帧图像间的运动可以近似为线性的，即，u＝V_x，v＝V_x，也就是认定在邻域区域内的N个像素点的速度是相同的，而且这N个像素点的速度都为V_x和V_x，将其代入(4)式，得到

该方程对邻域区域内中的N个像素点都成立，这样就可以得到N个(N为邻域区域中像素的个数)方程组成的方程组，用矩阵的形式表示如下：

上式中，第一个矩阵为N*2矩阵，等式右边的矩阵为N*1矩阵，这样，N个方程，两个未知数，分别对I求偏导数，就可以求得光流的速度场如下：

由此可以计算出光流场，且可知，在计算光流场时既考虑了连续时次回波的变化，也考虑了相邻位置回波的变化。

在S103中，依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

在本实施例中，将S102计算得到的雷达回波的光流场来代替交叉相关法得到的运动矢量场，从而根据雷达回波光流场的移动特征来外推雷达回波未来的位置和形状，从而完成指定区域内对流天气的临近预报。

获得平滑真实的雷达回波光流场是追踪雷达回波的移动特征并做好外推，从而实现临近预报的关键，由于所获得的雷达数据除了包含雷达真实信号，还包含了各种误差噪声等成分，噪声污染对回波算法的准确性和精度有很大的影响，因此，利用光流法进行雷达回波移动追踪时，除了要进行上述质量控制外，还需分别进行雷达回波噪声消除的质量控制和光流运动矢量场的优化，以获得更加真实和平滑的雷达回波和光流运动矢量场。

首先，作为本发明的一个实施例，对雷达回波进行噪声消除，以实现质量控制的方案如下：

中值滤波是一种非线性平滑技术，是基于排序统计理论的一种能有效抑制雷达回波噪声的非线性信号处理技术，采用中值滤波方法对雷达回波的测量值损失小，平滑效果好，边沿清晰，因此，采用中值滤波对天气雷达数字化回波图作平滑处理是一种切实可行的方法。中值滤波的基本原理是把雷达图像序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的雷达回波强度值更接近真实值，从而消除孤立的噪声点。

对天气雷达数字化回波图作平滑处理的具体方法，是采用预设结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升(或下降)的二维数据序列。二维中值滤波输出为g(x，y)＝med{f(x-k，y-1),(k，1∈w)}，其中f(x，y)，g(x，y)分别为原始的所述雷达拼图的数据和处理后的所述雷达拼图的数据，w为预设结构的二维滑动模板，通常为2*2，3*3区域，也可以是不同的形状，如线状，圆形，十字形，圆环形等。

接下来，作为本发明的一个实施例，对光流场进行运动矢量场的优化。研究表明，用光流法计算雷达回波的运动矢量时，还必须考虑光流场的噪声情况。如果噪声较强,说明光流场数据本身的置信度较低，需要更多的条件加强光流的约束。根据统计试验,算法对雷达拼图某点的光流矢量与其周围相邻若干个点的光流矢量平均值的偏差超过5m/s和/或夹角超过25°的光流矢量予以剔除。进一步地，为了保证一定范围内光流场的平滑和真实，对于与整个探测区域光流矢量大小的整体平均值相比而言，偏差达到15m/s的个别光流矢量,也给予剔除。

应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本发明实施例的实施过程构成任何限定。

通过对2014年天气过程的实验对比，对于锋面低槽系统带来的降水天气过程，本方案效果优于交叉相关算法，对30分钟的预报评价结果，光流法的击中率平均值比交叉相关法高3个百分点，成功临界指数值平均比交叉相关法高6个百分点，60分钟的预报评价结果和30分钟相似。对热带系统带来的降水天气个例，交叉相关法略胜光流法，但两者的效果差别不是太大。而移动型局地加强的对流降水过程，光流法效果优于交叉相关法，30分钟的预报评价结果看,光流法的击中率平均值比交叉相关法高7个百分点，成功临界指数值平均也比交叉相关法高7个百分点，60分钟的预报评价结果光流法效果也明显优于交叉相关法。移动型的分散多单体和超级单体结果和移动型局地加强的对流降水类似。雷暴变化速度较快时，光流法能计算出回波的运动矢量场。预报结果具有很好的业务指示意义，能提升对流临近预报系统的性能，有很大的应用推广价值。

对应于上文实施例所述的对流天气的临近预报方法，图2示出了本发明实施例提供的对流天气的临近预报系统的结构框图，为了便于说明，仅示出了与本实施例相关的部分。该系统为基于地理信息系统的可视化临近预报平台，基于B/S和C/S混合体系结构实现，C/S结构程序(后台程序)运行于服务器端，多线程并发执行，在得到雷达基数据后，启动基数据的预处理，进行雷达拼图、质量控制程序，算法运行得到实况和预报图像。客户端作为B/S结构程序(前台程序)，为显示平台用于显示所述服务器的临近预报结果。设计存储合理且拥有良好数据接口的数据库存储体系，其设计主要包括表结构(数据字段、数据类型、备注)、表间关系以及存储过程和应用层(如索引)的设计和构建。建设后台程序体系，运行于服务器端，其负责多线程并发执行。鉴于数据的时效性，后台程序多线程执行，并统一纳入后台程序体系，后台程序运行的状态可监控可在前台显示。

整个系统基于JAVA语言设计实现，程序在Linux PC服务器上运行，每个CPU至少有4个核。测试结果表明，每个流程所有程序在接收到雷达资料后2分钟之内全部完成，因此，该算法在时间上完全可以作实时预报工具。

参照图2，在服务器中，包括：

获取单元21，获取探测区域内的雷达拼图；

运算单元22，基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场，对于所述光流场的光流(u，v)，所述u，v为所述雷达拼图其中一点的光流在x，y方向上的分量，该点在时刻t的灰度值为I(x，y，t)，相应的光流约束方程为：

I_xu+I_yv+I_t＝0，

其中：同时，引入进一步的约束条件来求解所述光流约束方程：

预报单元23，依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

可选地，所述服务器还包括：

中值滤波单元，采用中值滤波对所述雷达拼图的数据做平滑处理。

可选地，所述中值滤波单元具体用于：

采用预设结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升或下降的二维数据序列g(x，y)＝med{f(x-k，y-1),(k，1∈w)}，其中，所述f(x，y)和g(x，y)分别为原始的所述雷达拼图的数据和处理后的所述雷达拼图的数据，所述w为所述预设结构的二维滑动模板。

可选地，所述预设结构的二维滑动模板的形状包括线状，圆形，十字形或圆环形。

可选地，所述运算单元22还用于：

若所述雷达拼图某点的光流矢量与其相邻点的光流矢量平均值的偏差或夹角超过预设阈值，则在所述雷达回波的光流场中剔除该点的光流矢量。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。实施例中的各功能单元、模块可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中，上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。另外，各功能单元、模块的具体名称也只是为了便于相互区分，并不用于限制本申请的保护范围。上述系统中单元、模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

在本发明所提供的实施例中，应该理解到，所揭露的装置和方法，可以通过其它的方式实现。例如，以上所描述的系统实施例仅仅是示意性的，例如，所述模块或单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通讯连接可以是通过一些接口，装置或单元的间接耦合或通讯连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。

所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实施例的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)或处理器(processor)执行本发明实施例各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例各实施例技术方案的精神和范围。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种对流天气的临近预报方法，其特征在于，包括：

获取探测区域内的雷达拼图；

基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场，对于所述光流场的光流(u，v)，所述u，v为所述雷达拼图其中一点的光流在x,y方向上的分量，该点在时刻t的灰度值为I(x,y,t)，相应的光流约束方程为：

I_xu+I_yv+I_t＝0，

其中：若该点在(t+Δt)时运动到(x+Δx,y+Δy)，引入进一步的约束条件来求解所述光流约束方程：

依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

2.如权利要求1所述的方法，其特征在于，在所述获取探测区域内的雷达拼图之后，所述基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场之前，所述方法还包括：

采用中值滤波对所述雷达拼图的数据做平滑处理。

3.如权利要求2所述的方法，其特征在于，所述采用中值滤波对所述雷达拼图的数据做平滑处理包括：

采用预设结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升或下降的二维数据序列g(x,y)＝med{f(x-k,y-1),(k,1∈ w)}，其中，所述f(x,y),和g(x,y)分别为原始的所述雷达拼图的数据和处理后的所述雷达拼图的数据，所述w为所述预设结构的二维滑动模板。

4.如权利要求3所述的方法，其特征在于，所述预设结构的二维滑动模板的形状包括线状，圆形，十字形或圆环形。

5.如权利要求1所述的方法，其特征在于，所述基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场还包括：

若所述雷达拼图某点的光流矢量与其相邻点的光流矢量平均值的偏差或夹角超过预设阈值，则在所述雷达回波的光流场中剔除该点的光流矢量。

6.一种对流天气的临近预报系统，其特征在于，所述系统为基于地理信息系统的可视化平台，包括服务器和客户端，所述客户端用于显示所述服务器的临近预报结果，所述服务器包括：

获取单元，用于获取探测区域内的雷达拼图；

运算单元，用于基于所述雷达拼图，采用局部约束的光流法计算得到雷达回波的光流场，对于所述光流场的光流(u，v)，所述u，v为所述雷达拼图其中一点的光流在x,y方向上的分量，该点在时刻t的灰度值为I(x,y,t)，相应的光流约束方程为：

I_xu+I_yv+I_t＝0，

其中：若该点在(t+Δt)时运动到(x+Δx,y+Δy)，引入进一步的约束条件来求解所述光流约束方程：

预报单元，用于依据所述雷达回波的光流场进行雷达回波外推，完成对流天气的临近预报。

7.如权利要求6所述的系统，其特征在于，所述服务器还包括：

中值滤波单元，用于采用中值滤波对所述雷达拼图的数据做平滑处理。

8.如权利要求7所述的系统，其特征在于，所述中值滤波单元具体用于：

采用预设结构的二维滑动模板，将板内像素按照像素值的大小进行排序，生成单调上升或下降的二维数据序列g(x,y)＝med{f(x-k,y-1),(k,1∈w)，其中，所述f(x,y),和g(x,y)分别为原始的所述雷达拼图的数据和处理后的所述雷达拼图的数据，所述w为所述预设结构的二维滑动模板。

9.如权利要求8所述的系统，其特征在于，所述预设结构的二维滑动模板的形状包括线状，圆形，十字形或圆环形。

10.如权利要求6所述的系统，其特征在于，所述运算单元还用于：

若所述雷达拼图某点的光流矢量与其相邻点的光流矢量平均值的偏差或夹角超过预设阈值，则在所述雷达回波的光流场中剔除该点的光流矢量。