CN108732570B - 基于粒子滤波融合算法的灾害性对流天气的临近预报方法 - Google Patents

基于粒子滤波融合算法的灾害性对流天气的临近预报方法 Download PDF

Info

Publication number
CN108732570B
CN108732570B CN201710259932.1A CN201710259932A CN108732570B CN 108732570 B CN108732570 B CN 108732570B CN 201710259932 A CN201710259932 A CN 201710259932A CN 108732570 B CN108732570 B CN 108732570B
Authority
CN
China
Prior art keywords
echo
radar echo
motion vector
vector field
particle
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710259932.1A
Other languages
English (en)
Other versions
CN108732570A (zh
Inventor
陈元昭
兰红平
曹春燕
陈训来
张文海
王明洁
徐婷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
METEOROLOGICAL BUREAU OF SHENZHEN MUNICIPALITY
Original Assignee
METEOROLOGICAL BUREAU OF SHENZHEN MUNICIPALITY
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by METEOROLOGICAL BUREAU OF SHENZHEN MUNICIPALITY filed Critical METEOROLOGICAL BUREAU OF SHENZHEN MUNICIPALITY
Priority to CN201710259932.1A priority Critical patent/CN108732570B/zh
Publication of CN108732570A publication Critical patent/CN108732570A/zh
Application granted granted Critical
Publication of CN108732570B publication Critical patent/CN108732570B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/95—Radar or analogous systems specially adapted for specific applications for meteorological use
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开一种基于粒子滤波融合算法的灾害性对流天气临近预报方法，包括如下步骤：对所述雷达回波资料进行去除噪音的质量控制预处理；将去噪后的雷达回波基数据分别采用基于Lucas–Kanade约束法的光流法和基于Harris角点算法追踪计算得到雷达回波运动矢量场，再采用粒子滤波算法对所述雷达回波运动矢量场进行融合得到逼近雷达回波真实运动轨迹的最优雷达回波运动矢量估计；利用所述最优雷达回波运动矢量估计采用半拉格朗日外推方法进行雷达回波的外推预报，最后获得预定区域整个雷达回波的预报图像。初步结果显示本发明的方法能够获得更接近实况、更真实的、更精细的雷达回波运动矢量场，从而提高灾害性对流天气临近预报的准确率。

Description

基于粒子滤波融合算法的灾害性对流天气的临近预报方法
技术领域
本发明涉及气象雷达资料数据处理及临近预报领域，尤其涉及一种基于粒子滤波融合算法的临近预报方法。
背景技术
临近预报一般是指0-6h对雷暴及其产生的灾害性对流天气的预报。目前业务上使用的0-1h临近预报主要是基于多普勒雷达的雷暴识别和自动外推临近预报技术。客观外推预报30分钟、60分钟或更长时间的雷暴单体位置往往比主观外推预报更有效。20世纪50年代，Ligda首先进行了雷达回波的外推工作。此后雷达自动外推技术得到不断发展。对流临近预报外推方法主要有单体质心法、交叉相关法和光流法。单体质心法是将雷暴看作一个三维单体进行识别、追踪、外推做临近预报。交叉相关法是计算连续时次雷达回波的相关系数，建立雷达回波的最佳拟合关系，来确定回波在过去时次的移动矢量，以此追踪该回波的移动特征，通过得到的回波移动矢量来对回波的位置和形状进行外推。单体质心外推法主要用于对流降水的外推预报。而交叉相关法既可以对对流降水进行预报，也可以追踪层状云降水，在气象业务部门得到比较广泛的使用。但交叉相关法对局地生成的回波及强度和形状变化较快的回波外推失败的情况会显著增加。近年来光流法逐渐应用在短临预报中，取得较好的应用效果。
对华南地区降水检验表明，光流法对移动型局地生成的回波及强度和形状随时间变化很快的雷暴效果优于交叉相关预报方法，而对热带系统降水，交叉相关法预报效果优于光流法，表明光流法仍存在一定的局限性。在未来几年把高时空分辨率的数值预报产品应用在临近预报的实际业务中的可能性依然较低，2020年之前临近预报仍然主要以雷达外推技术为主。如何将各种临近预报方法结合起来，取长补短，得到最优化的预报结果，是目前临近预报的一个难题。
发明内容
本发明要解决的技术问题是提供一种能够得到更逼近回波真实运动矢量场的方法，并利用获得的运动矢量场进行外推，以提高0-1h客观临近预报能力。
本发明采用如下技术方案解决上述技术问题：
设计一种基于粒子滤波融合算法的灾害性对流天气临近预报方法，该方法以分布在预报区域多个站点的气象雷达回波资料为临近预报的初始资料，包括如下步骤：
第一步：对预报区域内所述雷达回波资料进行质量控制预处理，包括：
S1：将雷达资料从极坐标格式利用最近邻居法和垂直方向的线性内插法相结合的插值法插值到三维直角坐标系中，采用“膨胀-侵蚀”法解决雷达回波资料所存在的回波空洞、数据缺失的问题，得到雷达回波基数据；
S2：采用双边滤波对所述雷达回波基数据进行去除噪音处理；
第二步：将去噪后的雷达回波基数据分别采用基于Lucas–Kanade约束法的光流法和基于Harris角点算法追踪计算得到第一和第二回波运动矢量场，再采用粒子滤波算法对所述第一和第二回波运动矢量场进行融合得到逼近雷达回波真实运动轨迹的最优雷达回波运动矢量估计；
第三步：利用所述最优雷达回波运动矢量估计采用半拉格朗日外推方法进行雷达回波的外推预报，最后获得预报区域整个雷达回波的预报图像；
所述半拉格朗日外推方法的平流公式如下：
Figure GDA0003179052010000031
公式中的
Figure GDA0003179052010000032
F分别为雷达回波的预报场和观测场，t0为预报开始时间，τ为预报时效，α为位移矢量；
采用半拉格朗日平流方法将预报时效τ分为N步进行外推，外推的时间步长为△t，采用前向外推，每一步的位移量α由下式迭代得到：
Figure GDA0003179052010000033
其中
Figure GDA0003179052010000034
为在格点
Figure GDA0003179052010000035
处的速度；α的初始值为0，总位移是N步位移的总和。
所述第一步S2中的双边滤波是基于中心像素与中心周围像素的亮度差值的加权，对相似的像素赋予较高权重，不相似的像素赋予较小的权重，具体算法如下：
所述双边滤波的具体算法为:
Figure GDA0003179052010000036
其中：
Figure GDA0003179052010000037
Figure GDA0003179052010000038
上式中，I、FB(I)分别为输入图像和滤波后图像，高斯核函数
Figure GDA0003179052010000039
代表以(x,y)为中心，以w为半径的矩形内点的空间相似度，σs为方差参数，取值3；
Figure GDA00031790520100000310
为以(x,y)为中心，以w为半径的矩形内点的像素相似度，σr为方差参数，取值0.1。
所述第二步包括如下步骤：
S1：采用Lucas–Kanade约束的光流法得到第一回波运动矢量场f(x)；
S2：采用基于Harris角点算法追踪回波得到第二回波运动矢量场h(x)；
所述Harris角点定义的基础是雷达回波图像灰度强度的二阶导数矩阵，对原始雷达回波图像全部像素求导，得到一个新的二阶导数图像矩阵，即Hessian矩阵：
Figure GDA0003179052010000041
其中，I为输入图像，x,y分别为直角坐标系中x方向和y方向坐标，对于Harris角点，使用每个像素邻域内的二阶导数图像的自相关矩阵来定义，其表达式如下：
Figure GDA0003179052010000042
上式中，Wi,j为权重比例，I(x+i,y+j)为二阶导数图像像素点；
S3：采用粒子滤波算法对第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行融合得到逼近雷达回波真实运动的最优运动矢量场估计。
所述第二步粒子滤波的具体算法如下：
回波矢量场可表示为：
Figure GDA0003179052010000043
上式中，f(x)、h(x)分别为状态方程和观测方程，xt、yt、mt、nt分别表示系统的状态、观测值、过程噪声和观测噪声；
在得到观测值yt后，估计状态变量xt的条件概率密度函数Ω(xt|yt)，在进行雷达回波状态跟踪时，假定回波当前时刻的状态xt只与上一时刻的状态xt-1有关，同时假定观测值yt只与t时刻的状态xt有关，即观测值之间互相独立；
在进行求解时引入了优化算法对粒子滤波参数进行寻优，使算法对整个参数空间进行高效搜索以获得最优解，通过对各参数时间更新与量测更新的反复迭代估计出Ω(xt|yt)的最优解；粒子滤波利用一系列带权值的空间随机采样，来逼近后验概率密度函数；后验概率密度函数采用一组加权的粒子
Figure GDA0003179052010000051
近似表示：
Figure GDA0003179052010000052
式中
Figure GDA0003179052010000053
为归一化权重，满足
Figure GDA0003179052010000054
其中，xk为一组加权的粒子样本，
Figure GDA0003179052010000055
为单个粒子样本，N为粒子数。
所述第二步中的S1：采用Lucas–Kanade约束的光流法获得第一回波运动矢量场f(x)具体包括如下步骤：
S1-1：假定某点的雷达回波在(t+△t)时段运动到(x+△x,y+△y)，在设定的时间间隔△t内，雷达回波的灰度值保持不变，可得到如下光流约束方程：
Ixu+Iyv+It＝0
上式中：
Figure GDA0003179052010000056
S1-2：采用Lucas-Kanade局部约束法作为计算雷达回波光流场的约束条件，可得到如下光流场的计算公式：
Figure GDA0003179052010000057
Figure GDA0003179052010000058
Figure GDA0003179052010000059
上述Ix,Iy,It分别为雷达回波数值在x,y,t方向相邻格点限制条件下的约束值；Ii，j，k为雷达回波在坐标(i,j)和时间k对应的灰度值，x，y分别为雷达回波格点在X方向和Y方向的坐标,t为时间坐标。
所述第二步的S3：采用粒子滤波算法对第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行融合得到逼近雷达回波真实运动的最优运动矢量场估计的具体算法包括如下步骤：
S3-1：初始化：k＝0，
在先验密度函数Ω(x0)中采集样本
Figure GDA0003179052010000061
该粒子权重取为
Figure GDA0003179052010000062
S3-2：k≥1时，按照如下步骤迭代计算：
S3-2-1：粒子重要性采样
粒子
Figure GDA0003179052010000063
被采样后，根据
Figure GDA0003179052010000064
分别计算相应的权重，并对其进行归一化处理
Figure GDA0003179052010000065
S3-2-2：重采样
如果有效采样粒子数N小于阈值Nth，那么对粒子样本
Figure GDA0003179052010000066
进行重新采样；重采样的样本形成新的粒子集
Figure GDA0003179052010000067
如果新的粒子集满足
Figure GDA0003179052010000068
则权重值重新定义为
Figure GDA0003179052010000069
否则由步骤S3-2-1跳转到步骤S3-2-3；
这里，
Figure GDA00031790520100000610
Nth为事前给定的阈值，取值为3N/4，P(·)为概率；
S3-2-3：状态更新
Figure GDA00031790520100000611
这里，
Figure GDA00031790520100000612
为相邻的单个粒子样本，q为粒子密度，
Figure GDA00031790520100000613
为单个粒子的权重值，p为粒子概率函数，zk为一组加权的粒子样本；
所述粒子滤波以递推的形式对所述第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行寻优融合，给出后验概率密度函数的最优解，从而计算出预报区域内回波运动的最优轨迹。
所述第三步采用半拉格朗日平流方法将预报时效τ分为N步进行外推的时间步长△t取值6分钟。
与现有技术相比，本发明基于粒子滤波融合算法的临近预报方法具有如下优点：结合粒子滤波法是基于Lucas–Kanade约束法的光流法和基于Harris角点算法追踪计算得到回波运动矢量场进行融合计算，可得到更逼近回波真实运动的最优矢量估计，能够获得更接近实况、更真实的、更精细的回波运动矢量场。
附图说明
图1是本发明基于粒子滤波融合算法的灾害性对流天气临近预报方法的流程图；
图2是2016年4月13日4：30时雷达回波双边滤波前后对比效果图，其中，(a)为滤波前的雷达回波效果图，(b)为滤波后的雷达回波效果图；
图3是对2016年4月13日6:30时广东范围的雷达回波采用多种算法获得的回波运动矢量场的对比效果图；其中，(a)为采用粒子滤波融合法获得的回波运动矢量场效果图，(b)为采用交叉相关法获得的回波运动矢量场效果图，(c)为采用光流法获得的回波运动矢量场效果图；
图4是对2016年5月21日广东范围降水过程采用多种算法的未来30分钟雷达回波图的外推预报的对比效果图；其中，(a)为起报时间为2:00的雷达回波实况，(b)为30分钟后2:30的雷达回波实况，(c)为采用粒子滤波融合法预报30分钟后2:30的回波，d为采用交叉相关法预报30分钟后2:30的回波，e为采用光流法预报30分钟后2:30的回波；
图5是2016年8月2日妮妲台风广东雷达拼图采用多种算法的未来30分钟外推预报个例对比效果图；(a)为起报时间为04:30的雷达实况，(b)为30分钟后04:30的雷达实况，(c)为采用粒子滤波融合法预报未来30分钟后05:00的回波，(d)为采用交叉相关法预报未来30分钟后05:00的回波，(e)为采用光流法预报未来30分钟后05:00的回波。
具体实施方式
为了使本发明的目的、技术方案、有益效果更加清楚明白，以下结合附图及实施例，对本发明做进一步详细说明。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。
本发明的基于粒子滤波融合算法的灾害性对流天气临近预报方法以分布在预定预报区域多个站点的气象雷达回波的拼图资料为临近预报的初始资料。例如，本发明的实施例选用广东12部S波段多普勒天气雷达2.5公里高度的拼图资料，本发明的的灾害性对流天气临近预报方法包括如下步骤：
第一步：首先对雷达资料进行质量控制处理，该质量控制处理的具体步骤和算法如下：
S1：将雷达资料从极坐标格式利用最近邻居法和垂直方向的线性内插法相结合的插值法插值到三维直角坐标系中，采用“膨胀-侵蚀”法解决回波空洞、缺失缺陷问题，得到雷达回波基数据。所述“膨胀-侵蚀”法属于现有技术，不赘述；
S2：再采用双边滤波对三维直角坐标系的雷达回波基数据进行质量控制，双边滤波的具体算法为:
Figure GDA0003179052010000081
其中：
Figure GDA0003179052010000082
Figure GDA0003179052010000083
公式(1)中，I、FB(I)分别为输入图像和滤波后图像，高斯核函数
Figure GDA0003179052010000084
代表以(x,y)为中心，以w为半径的矩形内点的空间相似度，σs为方差参数，取值3；
Figure GDA0003179052010000085
为以(x,y)为中心，以w为半径的矩形内点的像素相似度，σr为方差参数，取值0.1。
由算法可知，双边滤波是基于中心像素与中心周围像素的亮度差值的加权，对相似的像素赋予较高权重，不相似的像素赋予较小的权重。双边滤波信号保持优异，降噪能力强，并且能够完整保持回波的边缘。图2是2016年4月13日6:30雷达回波双边滤波前后对比效果图，图2a为滤波前的雷达回波效果图，图2b为滤波后的雷达回波效果图。由图2中不难看出，经双边滤波处理后，雷达回波没有出现明显失真，回波变得更平滑、回波边缘完整清晰，证明双边滤波的去噪效果良好，达到了较理想的滤波效果。
第二步：将去噪后的雷达回波基数据分别采用基于Lucas–Kanade约束法的光流法和基于Harris角点算法追踪计算得到第一回波运动矢量场和第二回波运动矢量场，再采用粒子滤波算法对所述第一回波运动矢量场和第二回波运动矢量场进行融合得到逼近回波真实运动的最优运动矢量估计，具体算法如下：
S1:采用Lucas–Kanade约束的光流法得到第一回波运动矢量场f(x)，具体步骤包括：
S1-1：假定某点在(t+△t)时运动到(x+△x,y+△y)，在一定时间间隔△t内灰度值保持不变，可得到光流约束方程：
Ixu+Iyv+It＝0 (4)
公式(4)中：
Figure GDA0003179052010000091
S1-2：再采用Lucas-Kanade局部约束法作为计算光流场的约束条件，可得到光流场的计算方法：
Figure GDA0003179052010000092
S2：采用基于Harris角点算法追踪回波得到第二回波运动矢量场h(x)。
在回波运动估计领域中，角点是一类信息含量足够丰富并且可以从前后时次的图像中识别出来的点。Harris角点定义的基础是图像灰度强度的二阶导数矩阵，对原始图像全部像素求导，得到一幅新的二阶导数图像，即Hessian矩阵：
Figure GDA0003179052010000093
对于Harris角点，使用每个像素邻域内的二阶导数图像的自相关矩阵来定义，其表达如下：
Figure GDA0003179052010000101
公式(7)中，Wi,j为权重比例，I(x+i,y+j)为二阶导数图像像素点。
S3：采用粒子滤波算法对所述第一回波运动矢量场和第二回波运动矢量场进行融合得到逼近回波真实运动的最优运动矢量估计。
粒子滤波方法为采用一组加权的粒子来近似表示状态变量的概率分布，然后再根据观测方程群迭代更新每一个粒子对应权重的大小来获取实际分布的样本，最后再对粒子进行重要性采样以获得均匀权重分布。粒子滤波的具体算法如下：
在粒子滤波算法中，回波矢量场可表示为：
Figure GDA0003179052010000102
公式(8)中，f(x)、h(x)分别为状态方程和观测方程。xt、yt、mt、nt分别表示系统的状态、观测值、过程噪声和观测噪声。
在得到观测值yt后，估计状态变量xt的条件概率密度函数Ω(xt|yt)。在进行雷达回波状态跟踪时，假定回波当前时刻的状态xt只与上一时刻的状态xt-1有关，同时假定观测值yt只与t时刻的状态xt有关，即观测值之间互相独立。
在进行求解时引入了优化算法对粒子滤波参数进行寻优，使算法对整个参数空间进行高效搜索以获得最优解。通过对各参数时间更新与量测更新的反复迭代估计出Ω(xt|yt)的最优解。粒子滤波利用一系列带权值的空间随机采样，来逼近后验概率密度函数，后验概率密度函数采用一组加权的粒子
Figure GDA0003179052010000103
近似表示：
Figure GDA0003179052010000104
公式(9)中，
Figure GDA0003179052010000105
为归一化权重，满足
Figure GDA0003179052010000106
其中，，xk为一组加权的粒子样本，
Figure GDA0003179052010000107
为单个粒子样本，N为粒子数。
采用粒子滤波算法估计后验概率密度函数主要分3步：首先从先验分布中抽取粒子进行采样，再计算粒子的权重进行定权，最后为避免粒子退化进行重采样，计算方法如下：
S3-1：初始化：k＝0，
在先验密度函数Ω(x0)中采集样本
Figure GDA0003179052010000111
该粒子权重取为
Figure GDA0003179052010000112
S3-2：k≥1时，按照如下步骤迭代计算：
S3-2-1：粒子重要性采样
粒子
Figure GDA0003179052010000113
被采样后，根据
Figure GDA0003179052010000114
分别计算相应的权重，并对其进行归一化处理
Figure GDA0003179052010000115
粒子滤波算法中的关键问题是粒子重要性分布函数的选择和再采集，本方法中选择先验密度作为重要性分布函数。
S3-2-2：重采样
如果有效采样粒子数N小于阈值Nth(其中
Figure GDA0003179052010000116
Nth为事前给定的阈值，本方法取值为3N/4)，那么对粒子样本
Figure GDA0003179052010000117
进行重新采样。重采样的样本形成新的粒子集
Figure GDA0003179052010000118
新的粒子集满足
Figure GDA0003179052010000119
(P(·)为概率)。权重值重新定义为
Figure GDA00031790520100001110
否则由步骤S3-2-1跳转到步骤S3-2-3。
S3-2-3：状态更新
Figure GDA00031790520100001111
粒子滤波以递推的形式对两种回波运动矢量场进行寻优融合，给出后验概率密度函数的最优解，从而可计算出回波运动的最优轨迹。图3是2016年4月13日6:30时回波运动矢量场的对比图，其中图3a为采用粒子滤波融合法获得的回波运动矢量场效果图，图3b为采用交叉相关法获得的回波运动矢量场效果图，图3c为采用光流法获得的回波运动矢量场效果图。结合图2和图3不难看出，粒子滤波融合法获得了更接近实况、更真实的、更精细的回波运动矢量。粒子滤波融合采用寻优的方法，把不同的运动矢量场结合起来，取长补短，得到最优化的预报结果，所以能够获得更好的回波运动矢量。
第三步：利用寻优融合后的回波运动矢量场采用半拉格朗日外推方法进行雷达回波的外推预报。
通过粒子滤波融合算法获得雷达回波的运动矢量场后就可以对回波进行临近外推预报。半拉格朗日平流公式为：
Figure GDA0003179052010000121
公式(11)中，
Figure GDA0003179052010000122
F分别为雷达回波的预报场和观测场，t0为预报开始时间，τ为预报时效，α为位移矢量。
采用半拉格朗日平流方法将预报时效τ分为N步进行外推，时间步长△t取值6分钟，采用前向外推，每一步的位移量α由下式迭代得到：
Figure GDA0003179052010000123
公式(12)中，
Figure GDA0003179052010000124
为在格点
Figure GDA0003179052010000125
处的速度。α的初始值为0，总位移是N步位移的总和。
利用粒子滤波融合算法得到的回波的运动矢量场，以此来进行外推预报，最终可获得整个回波的预报图像。
图4是针对2016年4月13日广东的一次强飑线过程，采用多种算法的未来60分钟外推预报的对比效果图，图4a为起始时间6:30的雷达实况，图4b为60分钟后7:30的回波实况，图4c为采用粒子滤波融合法预报60分钟后7:30的回波，图4d为采用交叉相关法预报60分钟后7:30回波，图4e为采用光流法预报60分钟后7:30回波；
图5为对2016年8月2日妮妲台风雷达拼图，基于多种算法的未来30分钟外推预报对比效果图，图5a为起报时间04:30的雷达实况，图5b为30分钟后04:30的雷达实况，图5c为采用粒子滤波融合法预报30分钟后05:00的回波，图5d为采用交叉相关法预报30分钟后05:00回波，图5e为采用光流法预报30分钟后05:00回波。
从图4和图5的对比效果图不难看出，粒子滤波融合法给出的30分钟、60分钟雷达回波的外推预报效果在大多数情况下好于采用交叉相关法和光流法的外推预报效果，粒子滤波融合法对对流天气过程雷达回波的外推预报具有更好的效果，可用于回波的临近预报，对业务有指导意义。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于粒子滤波融合算法的灾害性对流天气临近预报方法，该方法以分布在预报区域多个站点的气象雷达回波资料为临近预报的初始资料，包括如下步骤：
第一步：对预报区域内所述雷达回波资料进行质量控制预处理，包括：
S1：将雷达资料从极坐标格式利用最近邻居法和垂直方向的线性内插法相结合的插值法插值到三维直角坐标系中，采用“膨胀-侵蚀”法解决雷达回波资料所存在的回波空洞、数据缺失的问题，得到雷达回波基数据；
S2：采用双边滤波对所述雷达回波基数据进行去除噪音处理；
第二步：将去噪后的雷达回波基数据分别采用基于Lucas–Kanade约束法的光流法和基于Harris角点算法追踪计算得到第一和第二回波运动矢量场，再采用粒子滤波算法对所述第一和第二回波运动矢量场进行融合得到逼近雷达回波真实运动轨迹的最优雷达回波运动矢量估计；
第三步：利用所述最优雷达回波运动矢量估计采用半拉格朗日外推方法进行雷达回波的外推预报，最后获得预报区域整个雷达回波的预报图像；
所述半拉格朗日外推方法的平流公式如下：
Figure FDA0003179052000000011
公式中的
Figure FDA0003179052000000012
F分别为雷达回波的预报场和观测场，t0为预报开始时间，τ为预报时效，α为位移矢量；
采用半拉格朗日平流方法将预报时效τ分为N步进行外推，外推的时间步长为△t，采用前向外推，每一步的位移量α由下式迭代得到：
Figure FDA0003179052000000013
其中
Figure FDA0003179052000000014
为在格点
Figure FDA0003179052000000015
处的速度；α的初始值为0，总位移是N步位移的总和。
2.根据权利要求1所述的灾害性对流天气临近预报方法，其特征在于：所述第一步S2中的双边滤波是基于中心像素与中心周围像素的亮度差值的加权，对相似的像素赋予较高权重，不相似的像素赋予较小的权重，具体算法如下：
所述双边滤波的具体算法为:
Figure FDA0003179052000000021
其中：
Figure FDA0003179052000000022
Figure FDA0003179052000000023
上式中，I、FB(I)分别为输入图像和滤波后图像，高斯核函数
Figure FDA0003179052000000024
代表以(x,y)为中心，以w为半径的矩形内点的空间相似度，σs为方差参数，取值3；
Figure FDA0003179052000000025
为以(x,y)为中心，以w为半径的矩形内点的像素相似度，σr为方差参数，取值0.1。
3.根据权利要求1所述的灾害性对流天气临近预报方法，其特征在于：所述第二步包括如下步骤：
S1：采用Lucas–Kanade约束的光流法得到第一回波运动矢量场f(x)；
S2：采用基于Harris角点算法追踪回波得到第二回波运动矢量场h(x)；
所述Harris角点定义的基础是雷达回波图像灰度强度的二阶导数矩阵，对原始雷达回波图像全部像素求导，得到一个新的二阶导数图像矩阵，即Hessian矩阵：
Figure FDA0003179052000000026
其中，I为输入图像，x,y分别为直角坐标系中x方向和y方向坐标，对于Harris角点，使用每个像素邻域内的二阶导数图像的自相关矩阵来定义，其表达式如下：
Figure FDA0003179052000000031
上式中，Wi,j为权重比例，I(x+i,y+j)为二阶导数图像像素点；
S3：采用粒子滤波算法对第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行融合得到逼近雷达回波真实运动的最优运动矢量场估计。
4.根据权利要求1所述的灾害性对流天气临近预报方法，其特征在于：所述第二步粒子滤波的具体算法如下：
回波矢量场可表示为：
Figure FDA0003179052000000032
上式中，f(x)、h(x)分别为状态方程和观测方程，xt、yt、mt、nt分别表示系统的状态、观测值、过程噪声和观测噪声；
在得到观测值yt后，估计状态变量xt的条件概率密度函数Ω(xt|yt)，在进行雷达回波状态跟踪时，假定回波当前时刻的状态xt只与上一时刻的状态xt-1有关，同时假定观测值yt只与t时刻的状态xt有关，即观测值之间互相独立；
在进行求解时引入了优化算法对粒子滤波参数进行寻优，使算法对整个参数空间进行高效搜索以获得最优解，通过对各参数时间更新与量测更新的反复迭代估计出Ω(xt|yt)的最优解；粒子滤波利用一系列带权值的空间随机采样，来逼近后验概率密度函数；后验概率密度函数采用一组加权的粒子
Figure FDA0003179052000000033
近似表示：
Figure FDA0003179052000000034
式中
Figure FDA0003179052000000035
为归一化权重，满足
Figure FDA0003179052000000036
其中，xk为一组加权的粒子样本，
Figure FDA0003179052000000037
为单个粒子样本，N为粒子数。
5.根据权利要求3所述的灾害性对流天气临近预报方法，其特征在于：所述S1中采用Lucas–Kanade约束的光流法获得第一回波运动矢量场f(x)具体包括如下步骤：
S1-1：假定某点的雷达回波在(t+△t)时段运动到(x+△x,y+△y)，在设定的时间间隔△t内，雷达回波的灰度值保持不变，可得到如下光流约束方程：
Ixu+Iyv+It＝0
上式中：
Figure FDA0003179052000000041
S1-2：采用Lucas-Kanade局部约束法作为计算雷达回波光流场的约束条件，可得到如下光流场的计算公式：
Figure FDA0003179052000000042
Figure FDA0003179052000000043
Figure FDA0003179052000000044
上述Ix,Iy,It分别为雷达回波数值在x,y,t方向相邻格点限制条件下的约束值；Ii，j，k为雷达回波在坐标(i,j)和时间k对应的灰度值，x，y分别为雷达回波格点在X方向和Y方向的坐标,t为时间坐标。
6.根据权利要求3所述的灾害性对流天气临近预报方法，其特征在于：所述步骤S3采用粒子滤波算法对第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行融合得到逼近雷达回波真实运动的最优运动矢量场估计的具体算法包括如下步骤：
S3-1：初始化：k＝0，
在先验密度函数Ω(x0)中采集样本
Figure FDA0003179052000000045
该粒子权重取为
Figure FDA0003179052000000046
S3-2：k≥1时，按照如下步骤迭代计算：
S3-2-1：粒子重要性采样
粒子
Figure FDA0003179052000000047
被采样后，根据
Figure FDA0003179052000000048
分别计算相应的权重，并对其进行归一化处理
Figure FDA0003179052000000049
S3-2-2：重采样
如果有效采样粒子数N小于阈值Nth，那么对粒子样本
Figure FDA0003179052000000051
进行重新采样；重采样的样本形成新的粒子集
Figure FDA0003179052000000052
如果新的粒子集满足
Figure FDA0003179052000000053
则权重值重新定义为
Figure FDA0003179052000000054
否则由步骤S3-2-1跳转到步骤S3-2-3；
这里，
Figure FDA0003179052000000055
Nth为事前给定的阈值，取值为3N/4；
S3-2-3：状态更新
Figure FDA0003179052000000056
这里，
Figure FDA0003179052000000057
为相邻的单个粒子样本，q为粒子密度，
Figure FDA0003179052000000058
为单个粒子的权重值，p为粒子概率函数，zk为一组加权的粒子样本；
所述粒子滤波以递推的形式对所述第一回波运动矢量场f(x)和第二回波运动矢量场h(x)进行寻优融合，给出后验概率密度函数的最优解，从而计算出预报区域内回波运动的最优轨迹。
7.根据权利要求1所述的灾害性对流天气临近预报方法，其特征在于：所述第三步采用半拉格朗日平流方法将预报时效τ分为N步进行外推的时间步长△t取值6分钟。
CN201710259932.1A 2017-04-20 2017-04-20 基于粒子滤波融合算法的灾害性对流天气的临近预报方法 Active CN108732570B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710259932.1A CN108732570B (zh) 2017-04-20 2017-04-20 基于粒子滤波融合算法的灾害性对流天气的临近预报方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710259932.1A CN108732570B (zh) 2017-04-20 2017-04-20 基于粒子滤波融合算法的灾害性对流天气的临近预报方法

Publications (2)

Publication Number Publication Date
CN108732570A CN108732570A (zh) 2018-11-02
CN108732570B true CN108732570B (zh) 2021-10-19

Family

ID=63924313

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710259932.1A Active CN108732570B (zh) 2017-04-20 2017-04-20 基于粒子滤波融合算法的灾害性对流天气的临近预报方法

Country Status (1)

Country Link
CN (1) CN108732570B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110379207B (zh) * 2019-07-19 2020-11-06 知天(珠海横琴)气象科技有限公司 一种自动船舶气象信息发布系统及其方法
CN111538013A (zh) * 2020-05-12 2020-08-14 上海眼控科技股份有限公司 雷达回波外推方法、装置、计算机设备和存储介质
CN112363140B (zh) * 2020-11-05 2024-04-05 南京叁云科技有限公司 一种基于循环神经网络的热力约束外推客观订正方法
CN114333241B (zh) * 2021-12-08 2022-10-11 电子科技大学 基于事件触发的滑坡灾害点大数据采集与样本库更新方法
CN117630946B (zh) * 2024-01-26 2024-04-12 中国人民解放军国防科技大学 基于双偏振雷达的强对流联合观测指挥方法、系统及设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101399969A (zh) * 2007-09-28 2009-04-01 三星电子株式会社 基于运动相机的运动目标检测与跟踪的系统、设备和方法
CN104869387A (zh) * 2015-04-19 2015-08-26 中国传媒大学 基于光流法的双目图像最大视差获取方法
JP5862023B2 (ja) * 2011-03-04 2016-02-16 日本電気株式会社 目標追跡システム及び目標追跡方法
CN105741324A (zh) * 2016-03-11 2016-07-06 江苏物联网研究发展中心 移动平台上的运动目标检测识别与跟踪方法
WO2016156292A1 (en) * 2015-03-31 2016-10-06 Koninklijke Philips N.V. Filter assembly and airway pressure support system employing same

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101399969A (zh) * 2007-09-28 2009-04-01 三星电子株式会社 基于运动相机的运动目标检测与跟踪的系统、设备和方法
JP5862023B2 (ja) * 2011-03-04 2016-02-16 日本電気株式会社 目標追跡システム及び目標追跡方法
WO2016156292A1 (en) * 2015-03-31 2016-10-06 Koninklijke Philips N.V. Filter assembly and airway pressure support system employing same
CN104869387A (zh) * 2015-04-19 2015-08-26 中国传媒大学 基于光流法的双目图像最大视差获取方法
CN105741324A (zh) * 2016-03-11 2016-07-06 江苏物联网研究发展中心 移动平台上的运动目标检测识别与跟踪方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking;M. Sanjeev Arulampalam;《IEEE TRANSACTIONS ON SIGNAL PROCESSING,》;20020228;全文 *
粒子滤波理论、方法及其在多目标跟踪中的应用;李天成;《自动化学报》;20151231;全文 *

Also Published As

Publication number Publication date
CN108732570A (zh) 2018-11-02

Similar Documents

Publication Publication Date Title
CN108732570B (zh) 基于粒子滤波融合算法的灾害性对流天气的临近预报方法
CN108508505B (zh) 基于多尺度卷积神经网络的强降雨及雷暴预报方法和系统
CN106875415B (zh) 一种动态背景中弱小动目标的连续稳定跟踪方法
CN108460764B (zh) 基于自动上下文和数据增强的超声图像智能分割方法
CN107256225B (zh) 一种基于视频分析的热度图生成方法及装置
CN107146239B (zh) 卫星视频运动目标检测方法及系统
CN108761461B (zh) 基于天气雷达回波时序图像的降水预报方法
CN110084201B (zh) 一种监控场景下基于特定目标跟踪的卷积神经网络的人体动作识别方法
CN102142085B (zh) 一种林区监控视频中运动火焰目标的鲁棒跟踪方法
CN111275740B (zh) 一种基于高分辨率孪生网络的卫星视频目标跟踪方法
CN104156984A (zh) 一种不均匀杂波环境下多目标跟踪的概率假设密度方法
CN110717934B (zh) 一种基于strcf的抗遮挡目标跟踪方法
CN105374049B (zh) 一种基于稀疏光流法的多角点跟踪方法及装置
CN112180375A (zh) 一种基于改进的TrajGRU网络的气象雷达回波外推方法
CN107622507B (zh) 一种基于深度学习的空中目标跟踪方法
CN110376582B (zh) 自适应gm-phd的机动目标跟踪方法
CN111652790A (zh) 一种亚像素图像配准方法
CN114648547B (zh) 用于反无人机红外探测系统的弱小目标检测方法和装置
Yao et al. Prediction of weather radar images via a deep lstm for nowcasting
CN113191427B (zh) 一种多目标车辆跟踪方法及相关装置
CN108876807B (zh) 一种实时星载卫星图像运动对象检测跟踪方法
CN112668615B (zh) 一种基于深度跨尺度外推融合的卫星云图预测方法
CN103985139A (zh) 基于颜色模型与预测向量簇模型信息融合的粒子滤波目标跟踪方法
CN110751671B (zh) 一种基于核相关滤波与运动估计的目标跟踪方法
CN116012421A (zh) 目标跟踪方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant