CN108519631B - 降水强度预测方法 - Google Patents

降水强度预测方法 Download PDF

Info

Publication number
CN108519631B
CN108519631B CN201810153350.XA CN201810153350A CN108519631B CN 108519631 B CN108519631 B CN 108519631B CN 201810153350 A CN201810153350 A CN 201810153350A CN 108519631 B CN108519631 B CN 108519631B
Authority
CN
China
Prior art keywords
moment
pixel point
optical flow
image data
pixel
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
CN201810153350.XA
Other languages
English (en)
Other versions
CN108519631A (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.)
Qingdao Xinzhongyoushu Technology Co ltd
Original Assignee
Qingdao Xinzhongyoushu Technology Co ltd
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 Qingdao Xinzhongyoushu Technology Co ltd filed Critical Qingdao Xinzhongyoushu Technology Co ltd
Priority to CN201810153350.XA priority Critical patent/CN108519631B/zh
Publication of CN108519631A publication Critical patent/CN108519631A/zh
Application granted granted Critical
Publication of CN108519631B publication Critical patent/CN108519631B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/10—Devices for predicting weather conditions
    • 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
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30181—Earth observation
    • G06T2207/30192—Weather; Meteorology
    • 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)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Environmental & Geological Engineering (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Atmospheric Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明提供了一种降水强度预测方法，包括：获取天气雷达在T时刻、T‑t时刻和T‑2t时刻测量到的格点数据；将T时刻、T‑t时刻和T‑2t时刻的格点数据分别转化为灰度图像数据；分别对T时刻、T‑t时刻和T‑2t时刻的灰度图像数据进行滤波处理，得到T时刻、T‑t时刻和T‑2t时刻的平滑灰度图像数据；根据T时刻、T‑t时刻和T‑2t时刻的平滑灰度图像数据以及光流场函数，确定T时刻的估计光流场；将T时刻的平滑灰度图像数据和T时刻的估计光流场代入半拉格朗日后向外推法，计算得到多个预测灰度图像数据；根据降水强度计算公式，将多个预测灰度图像数据分别转换为降水强度数据，缓解现有技术中的降水强度预测结果准确性低的问题，达到了提高降水强度预测结果准确性的效果。

Description

降水强度预测方法
技术领域
本发明涉及降水强度预测技术领域，尤其是涉及一种降水强度预测方法。
背景技术
天气雷达是监测强对流天气和估测短时降水强度的主要工具。根据天气雷达测量得到的回波强度的格点数据，可以定量预测短时临近时间段内的降水天气。现有技术是将回波强度的格点数据视为灰度图像，基于计算视觉领域中的光流法，计算雷达回波图的光流场并进行外推，然后将回波强度值转换为降水强度值，进而得到临近时间段内的降水预测结果。现有光流法技术基于亮度恒定的假设，适用于无遮挡且像素点进行连续小幅运动的场景。但是，由于相邻时刻的雷达回波变换复杂，并且实际天气雷达回波图像存在杂波多和部分雷达数据没有按时传输等问题，所以，现有技术在进行降水强度的预测时，往往达不到预期的效果，不能得到准确的降水强度预测结果。
发明内容
有鉴于此，本发明的目的在于提供一种降水强度预测方法，以缓解现有技术中存在的降水强度预测结果准确性低的技术问题。
第一方面，本发明实施例提供了一种降水强度预测方法，包括：
获取天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据；
将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据；
分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据；
根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场；
将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到多个预测灰度图像数据，多个所述预测灰度图像数据是T时刻之后的图像数据；
根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据。
结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据，包括：
利用预设的线性变换公式，分别对所述T时刻、T-t时刻和T-2t时刻的格点数据中每一个格点的回波强度数值进行线性变换，分别得到T时刻、T-t时刻和T-2t时刻的所述灰度图像数据。
结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据，包括：
利用预设的双边滤波器，分别对所述T时刻、T-t时刻和T-2t时刻的灰度图像数据中的每一个像素点进行滤波，分别得到所述T时刻、T-t时刻和T-2t时刻的所述平滑灰度图像数据。
结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场，包括：
将所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第一光流场；
将所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第二光流场；
将所述第一光流场乘以预设的第一权重系数，得到第一权重光流场；
将所述第二光流场乘以预设的第二权重系数，得到第二权重光流场；
将所述第一权重光流场和所述第二权重光流场进行相加，得到所述T时刻的估计光流场。
结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述将所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第一光流场，包括：
确定所述T时刻的平滑灰度图像数据中各回波图像块的外层轮廓线、与所述外层轮廓线对应的内部区域，以及未处于所述外层轮廓线内部的无回波区域；
利用预设的FAST算法，计算所述T时刻的平滑灰度图像数据中的所有的角点；
在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-2t时刻的平滑灰度图像数据相匹配的第一匹配像素点集合中任一像素点的光流；
在所述第一匹配像素点集合中，剔除光流矢量长度大于预设的矢量长度阈值的像素点，得到第一稀疏像素点集合；
将与所述第一稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第一初始光流场；
将所述第一初始光流场、所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入预设的Dual TV-L1光流法，得到所述第一光流场。
结合第一方面，本发明实施例提供了第一方面的第五种可能的实施方式，其中，所述在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-2t时刻的平滑灰度图像数据相匹配的第一匹配像素点集合中任一像素点的光流，包括：
将位于所述外层轮廓线上的像素点和所述角点集合成第一像素点集合；
利用预设的Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-2t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第一匹配像素点集合，以及与所述第一匹配像素点集合中的像素点分别对应的位于所述T-2t时刻的平滑灰度图像数据中的匹配点集合；
计算所述第一匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻的平滑灰度图像数据中的坐标减去所述像素点在T-2t时刻的平滑灰度图像中对应的匹配点的坐标。
结合第一方面，本发明实施例提供了第一方面的第六种可能的实施方式，其中，所述在所述第一匹配像素点集合中，剔除光流矢量长度大于预设的矢量长度阈值的像素点，得到第一稀疏像素点集合，包括：
计算所述第一匹配像素点集合中每一个像素点的光流矢量长度；
在所述第一匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第一大于阈值像素点集合；
在所述第一匹配像素点集合中，剔除所述第一大于阈值像素点集合，得到所述第一稀疏像素点集合。
结合第一方面，本发明实施例提供了第一方面的第七种可能的实施方式，其中，所述将与所述第一稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第一初始光流场，包括：
在所述T时刻的平滑灰度图像数据中任取一个像素点；
判断所述像素点是否位于所述第一稀疏像素点集合中；
当所述像素点位于所述第一稀疏像素点集合中时，所述像素点的光流不变；
当所述像素点未位于所述第一稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中；
当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第一稀疏像素点集合中距离所述像素点最近的像素点的光流；
当所述像素点位于所述内部区域中时，确定既位于所述第一稀疏像素点集合，又位于所述内部区域的第一目标像素点；
在所述第一目标像素点中，挑选距离所述像素点最近的预设数量个像素点；
根据预设的加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第一初始光流场。
结合第一方面，本发明实施例提供了第一方面的第八种可能的实施方式，其中，所述将所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第二光流场，包括：
在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-t时刻的平滑灰度图像数据相匹配的第二匹配像素点集合中任一像素点的光流；
在所述第二匹配像素点集合中，剔除光流矢量长度大于所述矢量长度阈值的像素点，得到第二稀疏像素点集合；
将与所述第二稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第二初始光流场；
将所述第二初始光流场、所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述Dual TV-L1光流法，得到所述第二光流场。
结合第一方面，本发明实施例提供了第一方面的第九种可能的实施方式，其中，所述在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-t时刻的平滑灰度图像数据相匹配的第二匹配像素点集合中任一像素点的光流，包括：
利用所述Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第二匹配像素点集合，以及与所述第二匹配像素点集合中的像素点分别对应的位于所述T-t时刻的平滑灰度图像数据中的匹配点集合；
计算所述第二匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻平滑灰度图像数据中的坐标减去所述像素点在T-t时刻的平滑灰度图像中对应的匹配点的坐标。
结合第一方面，本发明实施例提供了第一方面的第十种可能的实施方式，其中，所述在所述第二匹配像素点集合中，剔除光流矢量长度大于所述矢量长度阈值的像素点，得到第二稀疏像素点集合，包括：
计算所述第二匹配像素点集合中每一个像素点的光流矢量长度；
在所述第二匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第二大于阈值像素点集合；
在所述第二匹配像素点集合中，剔除所述第二大于阈值像素点集合，得到所述第二稀疏像素点集合。
结合第一方面，本发明实施例提供了第一方面的第十一种可能的实施方式，其中，所述将与所述第二稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第二初始光流场，包括：
在所述T时刻平滑灰度图像数据中任取一个像素点；
判断所述像素点是否位于所述第二稀疏像素点集合中；
当所述像素点位于所述第二稀疏像素点集合中时，所述像素点的光流不变；
当所述像素点未位于所述第二稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中；
当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第二稀疏像素点集合中距离所述像素点最近的像素点的光流；
当所述像素点位于所述内部区域中时，确定既位于所述第二稀疏像素点集合，又位于所述内部区域的第二目标像素点；
在所述第二目标像素点中，挑选距离所述像素点最近的预设数量个像素点；
根据所述加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第二初始光流场。
结合第一方面，本发明实施例提供了第一方面的第十二种可能的实施方式，其中，所述根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据，包括：
利用预设的逆线性变换公式，分别对多个所述预测灰度图像数据进行逆线性变换，得到T时刻之后多个时刻的预测回波强度格点数据；
针对每一个所述预测回波强度格点数据，当所述预测回波强度格点数据中的格点的回波反射率数据小于0时，则所述格点的降水强度数据为0；
当所述预测回波强度格点数据中的格点的回波反射率数据大于或者等于0时，则利用所述降水强度计算公式计算得到所述格点的降水强度数据。
本发明实施例带来了以下有益效果：本发明实施例提供的降水强度预测方法包括：获取天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据；将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据；分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据；根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场；将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到T时刻之后多个时刻的预测灰度图像数据；根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据。
所以，当获取到天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据后，首先，将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据，然后，分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据，这样就可以去除实际天气雷达回波图像中存在的杂波，使得根据平滑灰度图像数据计算得到的降水强度预测结果更加准确，根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场；将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到T时刻之后多个时刻的预测灰度图像数据；考虑到天气雷达回波图像在实际业务中存在部分数据不能按时传输的情况，只使用两个相邻时刻的天气雷达回波数据计算光流场会导致部分雷达回波图像块的光流计算结果为0，本发明实施例中应用三个相邻时刻的天气雷达回波数据计算光流场，得到的光流预测结果更符合实际情况，根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据，这样，由于多个所述预测灰度图像数据更准确，更符合实际情况，所以根据多个所述预测灰度图像数据分别转换得到的降水强度数据也更准确，更加符合实际情况，避免由于只使用两个相邻时刻的天气雷达回波数据，且实际天气雷达回波图像存在杂波多和部分雷达数据没有按时传输而导致的不能得到准确的降水强度预测结果的问题，缓解了现有技术中的降水强度预测结果准确性低的技术问题，达到了提高降水强度预测结果准确性的技术效果。
本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的降水强度预测方法的流程图；
图2为图1中步骤S104的流程图；
图3为图2中步骤S201的流程图；
图4为图2中步骤S202的流程图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
目前，天气雷达是监测强对流天气和估测降水强度的主要工具。根据天气雷达测量得到的回波强度的格点数据，可以定量预测短时临近时间段内的降水天气。现有技术是将回波强度的格点数据视为灰度图像，基于计算视觉领域中的光流法，计算雷达回波图的光流场并进行外推，然后将回波强度值转换为降水强度值，进而得到临近时间段内的降水预测结果。现有光流法技术基于亮度恒定的假设，适用于无遮挡且像素点进行连续小幅运动的场景。但是，由于相邻时刻的雷达回波变换复杂(例如：回波运动为非刚性变形运动，存在突然出现或消散的回波等)，并且实际天气雷达回波图像存在杂波多和部分雷达数据没有按时传输等问题，所以，现有技术在进行降水强度的预测时，往往达不到预期的效果，不能得到准确的降水强度预测结果，基于此，本发明实施例提供的一种降水强度预测方法，可以缓解现有技术中的降水强度预测结果准确性低的技术问题，达到了提高降水强度预测结果准确性的技术效果。
为便于对本实施例进行理解，首先对本发明实施例所公开的一种降水强度预测方法进行详细介绍，如图1所示，所述降水强度预测方法可以包括以下步骤。
步骤S101，获取天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据。
示例性的，由于现有的天气雷达测量到的回波强度的格点数据都是以6分钟为间隔，所以三个相邻的时刻可以设置为T时刻、T-6min时刻和T-12min时刻。
步骤S102，将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据。
示例性的，利用预设的线性变换公式，分别对所述T时刻、T-t时刻和T-2t时刻的格点数据中每一个格点的回波强度数值进行线性变换，分别得到T时刻、T-t时刻和T-2t时刻的所述灰度图像数据。
示例性的，所述线性变换公式可以为I＝k*dBZ+b，其中I表示灰度图像数据中每个像素点的灰度值，dBZ表示格点数据中每个格点的反射率数值，k和b是已知常数。
步骤S103，分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据。
示例性的，利用预设的双边滤波器，分别对所述T时刻、T-t时刻和T-2t时刻的灰度图像数据中的每一个像素点进行滤波，分别得到所述T时刻、T-t时刻和T-2t时刻的所述平滑灰度图像数据。
示例性的，所述双边滤波器是对图像进行平滑的一种滤波器，所述双边滤波器的输入是一张图像，输出是一张经过滤波后的图像。所述双边滤波器具有三个预先设定的常数参数：窗口大小、位置平滑参数和像素值平滑参数。T时刻的灰度图像数据是输入图像，T时刻的平滑灰度图像数据是经过滤波后的输出图像。T-t时刻的灰度图像数据是输入图像，T-t时刻的平滑灰度图像数据是经过滤波后的输出图像。T-2t时刻的灰度图像数据是输入图像，T-2t时刻的平滑灰度图像数据是经过滤波后的输出图像。
步骤S104，根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场。
示例性的，如图2所示，步骤S104可以包括以下步骤。
步骤S201，将所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第一光流场。
示例性的，如图3所示，步骤S201可以包括以下步骤。
步骤S301，确定所述T时刻的平滑灰度图像数据中各回波图像块的外层轮廓线、与所述外层轮廓线对应的内部区域，以及未处于所述外层轮廓线内部的无回波区域。
步骤S302，利用预设的FAST算法，计算所述T时刻的平滑灰度图像数据中的所有的角点。
示例性的，所述FAST算法是检测图像特征点的一种常用方法。所述FAST算法的输入是一张图像，所述FAST算法的输出是检测到的图像中的特征点集合。如果在圆周上有N个连续的像素点的亮度都比圆心像素点的亮度加上阈值还要亮，或者比圆心像素点的亮度减去阈值还要暗，则圆心像素点被称为角点。其中，N的值、阈值和圆心像素点的亮度可以预先设置好。
步骤S303，在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-2t时刻的平滑灰度图像数据相匹配的第一匹配像素点集合中任一像素点的光流。
示例性的，步骤S303可以包括以下步骤。
步骤S401，将位于所述外层轮廓线上的像素点和所述角点集合成第一像素点集合。
步骤S402，利用预设的Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-2t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第一匹配像素点集合，以及与所述第一匹配像素点集合中的像素点分别对应的位于所述T-2t时刻的平滑灰度图像数据中的匹配点集合。
示例性的，所述的Lucas-Kanade多层金字塔稀疏特征跟踪算法是一种常用的稀疏光流算法，所述Lucas-Kanade多层金字塔稀疏特征跟踪算法的输入是：第一张图像、第一张图像中的特征点位置和第二张图像，输出是在第二张图像中具有相应匹配点的特征点集合(所述特征点集合位于第一张图像中)，以及与特征点集合中的特征点分别对应的位于第二张图像中的匹配点(匹配点带有坐标)。所述T时刻的平滑灰度图像数据是输入的第一张图像，第一像素点集合是第一张图像中的特征点位置，T-2t时刻的平滑灰度图像数据是第二张图像，输出是位于所述T时刻的平滑灰度图像数据中的所述第一匹配像素点集合，以及与所述第一匹配像素点集合中的像素点分别对应的位于所述T-2t时刻的平滑灰度图像数据中的匹配点集合。
步骤S403，计算所述第一匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻的平滑灰度图像数据中的坐标减去所述像素点在T-2t时刻的平滑灰度图像中对应的匹配点的坐标。
步骤S304，在所述第一匹配像素点集合中，剔除光流矢量长度大于预设的矢量长度阈值的像素点，得到第一稀疏像素点集合。
示例性的，步骤S304可以包括以下步骤。
步骤S501，计算所述第一匹配像素点集合中每一个像素点的光流矢量长度。
示例性的，每个像素点的光流包括：第一方向速度vx和第二方向速度vy。所以，每个像素点的光流矢量长度为
Figure BDA0001580692990000131
步骤S502，在所述第一匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第一大于阈值像素点集合。
示例性的，所述矢量长度阈值是预先设置的。
步骤S503，在所述第一匹配像素点集合中，剔除所述第一大于阈值像素点集合，得到所述第一稀疏像素点集合。
示例性的，在所述第一匹配像素点集合中，将所述光流矢量长度大于所述矢量长度阈值的像素点剔除掉，这样做的目的是对于所述第一匹配像素点集合中像素点的光流做质量控制。
步骤S305，将与所述第一稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第一初始光流场。
示例性的，步骤S305可以包括以下步骤。
步骤S601，在所述T时刻的平滑灰度图像数据中任取一个像素点。
步骤S602，判断所述像素点是否位于所述第一稀疏像素点集合中。
步骤S603，当所述像素点位于所述第一稀疏像素点集合中时，所述像素点的光流不变。
步骤S604，当所述像素点未位于所述第一稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中。
步骤S605，当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第一稀疏像素点集合中距离所述像素点最近的像素点的光流。
步骤S606，当所述像素点位于所述内部区域中时，确定既位于所述第一稀疏像素点集合，又位于所述内部区域的第一目标像素点。
步骤S607，在所述第一目标像素点中，挑选距离所述像素点最近的预设数量个像素点。
示例性的，所述预设数量可以是预先设置的。
步骤S608，根据预设的加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第一初始光流场。
示例性的，以p为位于所述内部区域中的任一像素点，所述预设数量为k为例进行说明。p的光流为从所述第一目标像素点中挑选出的距离p最近的k个像素点的光流的加权平均值。所述加权平均公式为：
Figure BDA0001580692990000151
其中，c为已知常数，flow(pi)为k个像素点中第i个像素点的光流，
Figure BDA0001580692990000152
(xp,yp)为已知的像素点p的像素坐标，(xi,yi)为已知的像素点i的像素坐标，i＝1,2,...,k。这k个像素点的光流的权重系数与这k个像素点和p之间的距离成反比。
示例性的，由于每个像素点的光流包括：第一方向速度vx和第二方向速度vy，所以，对于像素点p的光流的第一个分量：第一方向速度vpx而言，
Figure BDA0001580692990000153
其中，vix为k个像素点中第i个像素点的光流的第一个分量。对于像素点p的光流的第二个分量：第二方向速度vpy而言，
Figure BDA0001580692990000154
其中，viy为k个像素点中第i个像素点的光流的第二个分量。
步骤S306，将所述第一初始光流场、所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入预设的Dual TV-L1光流法，得到所述第一光流场。
示例性的，由于所述第一初始光流场是根据局部信息估计得到的，为了减小估计的全局误差，要将所述第一初始光流场作为输入，估计所述第一光流场。所述Dual TV-L1光流法是一种稠密光流算法。所述Dual TV-L1光流法的输入是第一张图像和第二张图像，所述Dual TV-L1光流法的输出是第一张图像的光流场，所述Dual TV-L1光流法还可以输入一个初始光流场，作为算法的起始估计。所述T时刻的平滑灰度图像数据就是输入的第一张图像，所述T-2t时刻的平滑灰度图像数据就是输入的第二张图像，所述第一初始光流场就是作为起始估计的初始光流场，因此，所述Dual TV-L1光流法的输出为所述T时刻的平滑灰度图像数据的光流场，即所述第一光流场。
本发明实施例中，由于强对流天气情形下降水物质的运动形式十分复杂，属于非刚性运动且存在突然出现或消散的情形，本发明实施例中通过匹配多张雷达回波图像中的角点，先估计图像中能够找到匹配的角点的光流矢量，根据这些点可以获取相对准确的光流矢量，提高降水强度预测过程的准确性。同时，鉴于不同雷达回波图像块具有不同的移动速度，本发明实施例计算所述第一光流场时，对雷达回波图像依照回波图像块的轮廓线进行分组，分别计算不同雷达回波图像块内的光流场，计算得到的光流场比现有技术中的光流算法应用于雷达回波图像的结果更为可靠，保证了光流场计算的准确性。
步骤S202，将所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第二光流场。
示例性的，如图4所示，步骤S202可以包括以下步骤。
步骤S701，在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-t时刻的平滑灰度图像数据相匹配的第二匹配像素点集合中任一像素点的光流。
示例性的，步骤S701可以包括以下步骤。
步骤S801，利用所述Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第二匹配像素点集合，以及与所述第二匹配像素点集合中的像素点分别对应的位于所述T-t时刻的平滑灰度图像数据中的匹配点集合。
示例性的，所述的Lucas-Kanade多层金字塔稀疏特征跟踪算法是一种常用的稀疏光流算法，所述Lucas-Kanade多层金字塔稀疏特征跟踪算法的输入是：第一张图像、第一张图像中的特征点位置和第二张图像，输出是在第二张图像中具有相应匹配点的特征点集合(所述特征点集合位于第一张图像中)，以及与特征点集合中的特征点分别对应的位于第二张图像中的匹配点(匹配点带有坐标)。所述T时刻的平滑灰度图像数据是输入的第一张图像，第一像素点集合是第一张图像中的特征点位置，T-t时刻的平滑灰度图像数据是第二张图像，输出是位于所述T时刻的平滑灰度图像数据中的所述第二匹配像素点集合，以及与所述第二匹配像素点集合中的像素点分别对应的位于所述T-t时刻的平滑灰度图像数据中的匹配点集合。
步骤S802，计算所述第二匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻平滑灰度图像数据中的坐标减去所述像素点在T-t时刻的平滑灰度图像中对应的匹配点的坐标。
步骤S702，在所述第二匹配像素点集合中，剔除光流矢量长度大于所述矢量长度阈值的像素点，得到第二稀疏像素点集合。
示例性的，步骤S702可以包括以下步骤。
步骤S901，计算所述第二匹配像素点集合中每一个像素点的光流矢量长度。
步骤S902，在所述第二匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第二大于阈值像素点集合。
示例性的，所述矢量长度阈值是预先设置的。
步骤S903，在所述第二匹配像素点集合中，剔除所述第二大于阈值像素点集合，得到所述第二稀疏像素点集合。
示例性的，在所述第二匹配像素点集合中，将所述光流矢量长度大于所述矢量长度阈值的像素点剔除掉，这样做的目的是对于所述第二匹配像素点集合中像素点的光流做质量控制。
步骤S703，将与所述第二稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第二初始光流场。
示例性的，步骤S703可以包括以下步骤。
步骤S1001，在所述T时刻平滑灰度图像数据中任取一个像素点。
步骤S1002，判断所述像素点是否位于所述第二稀疏像素点集合中。
步骤S1003，当所述像素点位于所述第二稀疏像素点集合中时，所述像素点的光流不变。
步骤S1004，当所述像素点未位于所述第二稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中。
步骤S1005，当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第二稀疏像素点集合中距离所述像素点最近的像素点的光流。
步骤S1006，当所述像素点位于所述内部区域中时，确定既位于所述第二稀疏像素点集合，又位于所述内部区域的第二目标像素点。
步骤S1007，在所述第二目标像素点中，挑选距离所述像素点最近的预设数量个像素点。
示例性的，所述预设数量可以是预先设置的。
步骤S1008，根据所述加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第二初始光流场。
示例性的，以p为位于所述内部区域中的任一像素点，所述预设数量为k为例进行说明。p的光流为从所述第二目标像素点中挑选出的距离p最近的k个像素点的光流的加权平均值。所述加权平均公式为：
Figure BDA0001580692990000191
其中，c为已知常数，flow(pi)为k个像素点中第i个像素点的光流，
Figure BDA0001580692990000192
(xp,yp)为已知的像素点p的像素坐标，(xi,yi)为已知的像素点i的像素坐标，i＝1,2,...,k。这k个像素点的光流的权重系数与这k个像素点和p之间的距离成反比。
示例性的，由于每个像素点的光流包括：第一方向速度vx和第二方向速度vy，所以，对于像素点p的光流的第一个分量：第一方向速度vpx而言，
Figure BDA0001580692990000193
其中，vix为k个像素点中第i个像素点的光流的第一个分量。对于像素点p的光流的第二个分量：第二方向速度vpy而言，
Figure BDA0001580692990000194
其中，viy为k个像素点中第i个像素点的光流的第二个分量。
步骤S704，将所述第二初始光流场、所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述Dual TV-L1光流法，得到所述第二光流场。
示例性的，所述Dual TV-L1光流法是一种稠密光流算法。所述Dual TV-L1光流法的输入是第一张图像和第二张图像，所述Dual TV-L1光流法的输出是第一张图像的光流场，所述Dual TV-L1光流法还可以输入一个初始光流场，作为算法的起始估计。所述T时刻的平滑灰度图像数据就是输入的第一张图像，所述T-t时刻的平滑灰度图像数据就是输入的第二张图像，所述第二初始光流场就是作为起始估计的初始光流场，因此，所述Dual TV-L1光流法的输出为所述T时刻的平滑灰度图像数据的光流场，即所述第二光流场。
步骤S203，将所述第一光流场乘以预设的第一权重系数，得到第一权重光流场。
示例性的，所述第一权重系数可以是0.25。
步骤S204，将所述第二光流场乘以预设的第二权重系数，得到第二权重光流场。
示例性的，所述第二权重系数可以是0.5。
步骤S205，将所述第一权重光流场和所述第二权重光流场进行相加，得到所述T时刻的估计光流场。
步骤S105，将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到多个预测灰度图像数据，多个所述预测灰度图像数据是T时刻之后的图像数据。
示例性的，所述半拉格朗日后向外推法是给定当前雷达回波图像和雷达回波速度场(速度场等于光流场除以间隔时间)后预测未来的回波图像的方法，是用雷达回波预测降水强度时常用的一种外推方法。所述半拉格朗日后向外推法的输入是当前回波图像和速度场，输出是未来若干个时刻的雷达回波图像。在本发明实施例中，所述半拉格朗日后向外推法的输入是：T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场(光流场除以间隔时间t等于速度场)，输出是T时刻之后若干个时刻的预测灰度图像数据。
示例性的，当时间间隔t＝6min时，将所述T时刻的估计光流场除以6就得到每分钟的速度场，将所述T时刻的所述平滑灰度图像数据和所述每分钟的速度场代入所述半拉格朗日后向外推法，输出为T时刻之后每隔一分钟的预测灰度图像数据。
步骤S106，根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据。
示例性的，步骤S106可以包括以下步骤。
步骤S1101，利用预设的逆线性变换公式，分别对多个所述预测灰度图像数据进行逆线性变换，得到T时刻之后多个时刻的预测回波强度格点数据。
示例性的，所述逆线性变换公式可以是：
Figure BDA0001580692990000211
其中，I是所述预测灰度图像数据中每个像素点的灰度值，k和b是已知常数，和所述线性变换公式中的k和b是一样的。dBZ是所述预测格点数据中每个格点的反射率数值。
步骤S1102，针对每一个所述预测回波强度格点数据，当所述预测回波强度格点数据中的格点的回波反射率数据小于0时，则所述格点的降水强度数据为0。
步骤S1103，当所述预测回波强度格点数据中的格点的回波反射率数据大于或者等于0时，则利用所述降水强度计算公式计算得到所述格点的降水强度数据。
示例性的，所述降水强度计算公式可以为：
Figure BDA0001580692990000212
其中，c1和c2都是已知的常量参数，R为降水强度值。这样，就把多个所述预测灰度图像数据中的每个像素点的灰度值分别转换为降水强度数据，根据每个像素点的降水强度数据，就可以进行相应时间的降水强度预测。
本发明实施例中，本发明实施例提供的降水强度预测方法包括：获取天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据；将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据；分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据；根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场；将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到多个预测灰度图像数据，多个所述预测灰度图像数据是T时刻之后的图像数据；根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据。
所以，当获取到天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据后，首先，将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据，然后，分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据，这样就可以去除实际天气雷达回波图像中存在的杂波，使得根据平滑灰度图像数据计算得到的降水强度预测结果更加准确，根据T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据以及预设的光流场函数，确定T时刻的估计光流场；将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到多个预测灰度图像数据，多个所述预测灰度图像数据是T时刻之后的图像数据；考虑到天气雷达回波图像在实际业务中存在部分数据不能按时传输的情况，只使用两个相邻时刻的天气雷达回波数据计算光流场会导致部分雷达回波图像块的光流计算结果为0，本发明实施例中应用三个相邻时刻的天气雷达回波数据计算光流场，得到的光流预测结果更符合实际情况，根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据，这样，由于多个所述预测灰度图像数据更准确，更符合实际情况，所以根据多个所述预测灰度图像数据分别转换得到的降水强度数据也更准确，更加符合实际情况，避免由于只使用两个相邻时刻的天气雷达回波数据，且实际天气雷达回波图像存在杂波多和部分雷达数据没有按时传输而导致的不能得到准确的降水强度预测结果的问题，缓解了现有技术中的降水强度预测结果准确性低的技术问题，达到了提高降水强度预测结果准确性的技术效果。
除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对步骤、数字表达式和数值并不限制本发明的范围。
本发明实施例所提供的装置，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，装置实施例部分未提及之处，可参考前述方法实施例中相应内容。
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在这里示出和描述的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制，因此，示例性实施例的其他示例可以具有不同的值。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
另外，在本发明实施例的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
本发明实施例所提供的进行降水强度预测方法的计算机程序产品，包括存储了处理器可执行的非易失的程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (12)

1.一种降水强度预测方法，其特征在于，包括：
获取天气雷达在T时刻、T-t时刻和T-2t时刻测量到的回波强度的格点数据；
将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据；
分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据；
将所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入光流场函数，计算得到第一光流场；将所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第二光流场；将所述第一光流场乘以预设的第一权重系数，得到第一权重光流场；将所述第二光流场乘以预设的第二权重系数，得到第二权重光流场；将所述第一权重光流场和所述第二权重光流场进行相加，得到所述T时刻的估计光流场；
将T时刻的所述平滑灰度图像数据和所述T时刻的估计光流场代入预设的半拉格朗日后向外推法，计算得到多个预测灰度图像数据，多个所述预测灰度图像数据是T时刻之后的图像数据；
根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据。
2.根据权利要求1所述的降水强度预测方法，其特征在于，所述将T时刻、T-t时刻和T-2t时刻的格点数据分别转化为灰度图像数据，包括：
利用预设的线性变换公式，分别对所述T时刻、T-t时刻和T-2t时刻的格点数据中每一个格点的回波强度数值进行线性变换，分别得到T时刻、T-t时刻和T-2t时刻的所述灰度图像数据。
3.根据权利要求2所述的降水强度预测方法，其特征在于，所述分别对T时刻、T-t时刻和T-2t时刻的灰度图像数据进行滤波处理，得到T时刻、T-t时刻和T-2t时刻的平滑灰度图像数据，包括：
利用预设的双边滤波器，分别对所述T时刻、T-t时刻和T-2t时刻的灰度图像数据中的每一个像素点进行滤波，分别得到所述T时刻、T-t时刻和T-2t时刻的所述平滑灰度图像数据。
4.根据权利要求1所述的降水强度预测方法，其特征在于，所述将所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第一光流场，包括：
确定所述T时刻的平滑灰度图像数据中各回波图像块的外层轮廓线、与所述外层轮廓线对应的内部区域，以及未处于所述外层轮廓线内部的无回波区域；
利用预设的FAST算法，计算所述T时刻的平滑灰度图像数据中的所有的角点；
在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-2t时刻的平滑灰度图像数据相匹配的第一匹配像素点集合中任一像素点的光流；
在所述第一匹配像素点集合中，剔除光流矢量长度大于预设的矢量长度阈值的像素点，得到第一稀疏像素点集合；
将与所述第一稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第一初始光流场；
将所述第一初始光流场、所述T时刻的平滑灰度图像数据和所述T-2t时刻的平滑灰度图像数据代入预设的Dual TV-L1光流法，得到所述第一光流场。
5.根据权利要求4所述的降水强度预测方法，其特征在于，所述在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-2t时刻的平滑灰度图像数据相匹配的第一匹配像素点集合中任一像素点的光流，包括：
将位于所述外层轮廓线上的像素点和所述角点集合成第一像素点集合；
利用预设的Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-2t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第一匹配像素点集合，以及与所述第一匹配像素点集合中的像素点分别对应的位于所述T-2t时刻的平滑灰度图像数据中的匹配点集合；
计算所述第一匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻的平滑灰度图像数据中的坐标减去所述像素点在T-2t时刻的平滑灰度图像中对应的匹配点的坐标。
6.根据权利要求5所述的降水强度预测方法，其特征在于，所述在所述第一匹配像素点集合中，剔除光流矢量长度大于预设的矢量长度阈值的像素点，得到第一稀疏像素点集合，包括：
计算所述第一匹配像素点集合中每一个像素点的光流矢量长度；
在所述第一匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第一大于阈值像素点集合；
在所述第一匹配像素点集合中，剔除所述第一大于阈值像素点集合，得到所述第一稀疏像素点集合。
7.根据权利要求6所述的降水强度预测方法，其特征在于，所述将与所述第一稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第一初始光流场，包括：
在所述T时刻的平滑灰度图像数据中任取一个像素点；
判断所述像素点是否位于所述第一稀疏像素点集合中；
当所述像素点位于所述第一稀疏像素点集合中时，所述像素点的光流不变；
当所述像素点未位于所述第一稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中；
当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第一稀疏像素点集合中距离所述像素点最近的像素点的光流；
当所述像素点位于所述内部区域中时，确定既位于所述第一稀疏像素点集合，又位于所述内部区域的第一目标像素点；
在所述第一目标像素点中，挑选距离所述像素点最近的预设数量个像素点；
根据预设的加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第一初始光流场。
8.根据权利要求5所述的降水强度预测方法，其特征在于，所述将所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述光流场函数，计算得到第二光流场，包括：
在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-t时刻的平滑灰度图像数据相匹配的第二匹配像素点集合中任一像素点的光流；
在所述第二匹配像素点集合中，剔除光流矢量长度大于所述矢量长度阈值的像素点，得到第二稀疏像素点集合；
将与所述第二稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第二初始光流场；
将所述第二初始光流场、所述T时刻的平滑灰度图像数据和所述T-t时刻的平滑灰度图像数据代入所述Dual TV-L1光流法，得到所述第二光流场。
9.根据权利要求8所述的降水强度预测方法，其特征在于，所述在所述T时刻的平滑灰度图像数据中，根据位于所述外层轮廓线上的像素点和所述角点，计算与所述T-t时刻的平滑灰度图像数据相匹配的第二匹配像素点集合中任一像素点的光流，包括：
利用所述Lucas-Kanade多层金字塔稀疏特征跟踪算法，将所述第一像素点集合中的像素点和所述T-t时刻的平滑灰度图像数据进行匹配，得到位于所述T时刻的平滑灰度图像数据中的所述第二匹配像素点集合，以及与所述第二匹配像素点集合中的像素点分别对应的位于所述T-t时刻的平滑灰度图像数据中的匹配点集合；
计算所述第二匹配像素点集合中任一像素点的光流，所述像素点的光流等于所述像素点在所述T时刻平滑灰度图像数据中的坐标减去所述像素点在T-t时刻的平滑灰度图像中对应的匹配点的坐标。
10.根据权利要求9所述的降水强度预测方法，其特征在于，所述在所述第二匹配像素点集合中，剔除光流矢量长度大于所述矢量长度阈值的像素点，得到第二稀疏像素点集合，包括：
计算所述第二匹配像素点集合中每一个像素点的光流矢量长度；
在所述第二匹配像素点集合中，挑选出所述光流矢量长度大于所述矢量长度阈值的像素点，得到第二大于阈值像素点集合；
在所述第二匹配像素点集合中，剔除所述第二大于阈值像素点集合，得到所述第二稀疏像素点集合。
11.根据权利要求10所述的降水强度预测方法，其特征在于，所述将与所述第二稀疏像素点集合对应的稀疏光流插值为稠密光流，得到第二初始光流场，包括：
在所述T时刻平滑灰度图像数据中任取一个像素点；
判断所述像素点是否位于所述第二稀疏像素点集合中；
当所述像素点位于所述第二稀疏像素点集合中时，所述像素点的光流不变；
当所述像素点未位于所述第二稀疏像素点集合中时，判断所述像素点是否位于所述无回波区域中；
当所述像素点位于所述无回波区域中时，则确定所述像素点的光流为位于所述第二稀疏像素点集合中距离所述像素点最近的像素点的光流；
当所述像素点位于所述内部区域中时，确定既位于所述第二稀疏像素点集合，又位于所述内部区域的第二目标像素点；
在所述第二目标像素点中，挑选距离所述像素点最近的预设数量个像素点；
根据加权平均公式，确定所述像素点的光流为所述预设数量个像素点的光流的加权平均值，得到所述第二初始光流场。
12.根据权利要求11所述的降水强度预测方法，其特征在于，所述根据预设的降水强度计算公式，将多个所述预测灰度图像数据分别转换为降水强度数据，包括：
利用预设的逆线性变换公式，分别对多个所述预测灰度图像数据进行逆线性变换，得到T时刻之后多个时刻的预测回波强度格点数据；
针对每一个所述预测回波强度格点数据，当所述预测回波强度格点数据中的格点的回波反射率数据小于0时，则所述格点的降水强度数据为0；
当所述预测回波强度格点数据中的格点的回波反射率数据大于或者等于0时，则利用所述降水强度计算公式计算得到所述格点的降水强度数据。
CN201810153350.XA 2018-02-22 2018-02-22 降水强度预测方法 Active CN108519631B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810153350.XA CN108519631B (zh) 2018-02-22 2018-02-22 降水强度预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810153350.XA CN108519631B (zh) 2018-02-22 2018-02-22 降水强度预测方法

Publications (2)

Publication Number Publication Date
CN108519631A CN108519631A (zh) 2018-09-11
CN108519631B true CN108519631B (zh) 2020-09-25

Family

ID=63433170

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810153350.XA Active CN108519631B (zh) 2018-02-22 2018-02-22 降水强度预测方法

Country Status (1)

Country Link
CN (1) CN108519631B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110378936B (zh) * 2019-07-30 2021-11-05 北京字节跳动网络技术有限公司 光流计算方法、装置及电子设备
CN110824451A (zh) * 2019-11-20 2020-02-21 上海眼控科技股份有限公司 雷达回波图的处理方法、装置、计算机设备及存储介质
CN111142109A (zh) * 2019-12-30 2020-05-12 上海眼控科技股份有限公司 标记方法、装置、计算机设备和存储介质
CN113296074B (zh) * 2021-07-28 2022-02-22 成都远望探测技术有限公司 一种基于气象雷达多层cappi的光流法外推方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6448326B2 (ja) * 2014-11-20 2019-01-09 日本無線株式会社 気象レーダ表示装置および気象レーダ表示プログラム
CN104657994B (zh) * 2015-02-13 2017-12-19 厦门美图之家科技有限公司 一种基于光流法判断图像一致性的方法和系统

Also Published As

Publication number Publication date
CN108519631A (zh) 2018-09-11

Similar Documents

Publication Publication Date Title
CN108519631B (zh) 降水强度预测方法
KR101837407B1 (ko) 영상 기반 표적 추적 장치 및 추적 방법
EP2858008B1 (en) Target detecting method and system
KR100702663B1 (ko) 입자 필터 프레임 워크에서 전방향 시각센서 기반 위치추정 및 매핑을 위한 방법
US20190065824A1 (en) Spatial data analysis
CN112947419B (zh) 避障方法、装置及设备
CN106203381A (zh) 一种行车中障碍物检测方法与装置
KR101051389B1 (ko) 적응적 배경 기반의 객체 검출 및 추적 장치 및 방법
KR20180027242A (ko) 무인 차량을 위한 주변 환경 매핑 방법 및 장치
CN110942484A (zh) 基于遮挡感知和特征金字塔匹配的相机自运动估计方法
CN110717934A (zh) 一种基于strcf的抗遮挡目标跟踪方法
CN112287824A (zh) 基于双目视觉的三维目标检测方法、装置及系统
CN108010065A (zh) 低空目标快速检测方法及装置、存储介质及电子终端
CN111538013A (zh) 雷达回波外推方法、装置、计算机设备和存储介质
CN111723597B (zh) 跟踪算法的精度检测方法、装置、计算机设备和存储介质
CN117313828A (zh) 一种感知模型的迁移方法、装置、设备及介质
KR100635883B1 (ko) 실시간 물체 추적 장치
JP7199974B2 (ja) パラメータ選定装置、パラメータ選定方法、およびパラメータ選定プログラム
CN114549768B (zh) 一种三维重建效果检测方法、装置、设备及存储介质
CN115825946A (zh) 基于无监督学习的毫米波雷达测距方法及装置
CN113450385A (zh) 一种夜间工作工程机械视觉跟踪方法、装置及存储介质
CN112085685A (zh) 一种可消除砖块效应的基于空间混合分解的时空融合方法
CN115431968B (zh) 车辆控制器、车辆及车辆控制方法
CN111123267A (zh) 合成孔径雷达图像舰船检测方法及装置
RILL Intuitive Estimation of Speed using Motion and Monocular Depth Information

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220819

Address after: Unit 620, No. 37, Lianyungang, Shibei District, Qingdao City, Shandong Province, 266000

Patentee after: China Value (Qingdao) Meteorological Technology Co.,Ltd.

Address before: No. 962, Zhonglao Road, Licang District, Qingdao City, Shandong Province, 266000

Patentee before: QINGDAO XINZHONGYOUSHU TECHNOLOGY Co.,Ltd.

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230321

Address after: Room 710, Human Resources Building, No. 9, Yinchuan East Road, Laoshan District, Qingdao, Shandong 266100

Patentee after: QINGDAO XINZHONGYOUSHU TECHNOLOGY Co.,Ltd.

Address before: Unit 620, No. 37, Lianyungang, Shibei District, Qingdao City, Shandong Province, 266000

Patentee before: China Value (Qingdao) Meteorological Technology Co.,Ltd.