CN109100722A - 基于雷达回波图像扇区分量分析的风暴趋势预测方法 - Google Patents

基于雷达回波图像扇区分量分析的风暴趋势预测方法 Download PDF

Info

Publication number
CN109100722A
CN109100722A CN201810826357.3A CN201810826357A CN109100722A CN 109100722 A CN109100722 A CN 109100722A CN 201810826357 A CN201810826357 A CN 201810826357A CN 109100722 A CN109100722 A CN 109100722A
Authority
CN
China
Prior art keywords
sector
value
radar
storm
key
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.)
Granted
Application number
CN201810826357.3A
Other languages
English (en)
Other versions
CN109100722B (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.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
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 Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN201810826357.3A priority Critical patent/CN109100722B/zh
Publication of CN109100722A publication Critical patent/CN109100722A/zh
Application granted granted Critical
Publication of CN109100722B publication Critical patent/CN109100722B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/10—Devices for predicting weather conditions
    • 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)
  • Environmental & Geological Engineering (AREA)
  • General Physics & Mathematics (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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提出一种基于雷达回波图像扇区分量分析的风暴趋势预测方法，该方法与当前行业广泛采用的TREC、TITAN等风暴识别追踪算法不同，本发明方法不依赖于对雷达强回波中心的识别，特别对于弱回波或没有明显风暴核的回波，也能够较好地分析出风暴过去若干时刻的位移矢量，进而为准确预测奠定了基础。本发明方法直接对极坐标形式的基本反射率数据进行计算，没有采用空间插值等高空间复杂度的算法进行平面直角坐标系的转换，减少了相关计算的运算量，进而提升了短临天气预报业务的实时性和效能。

Description

基于雷达回波图像扇区分量分析的风暴趋势预测方法
技术领域
本发明属于大气科学领域，尤其涉及一种利用雷达探测资料进行风暴移动发展趋势预测的方法。
背景技术
气象上的风暴一般是指强烈天气系统过境时出现的天气过程，通常伴随有大风、雷暴、冰雹或短时强降水等天气现象。强热带风暴还可能引发雷暴大风、龙卷、下击暴流和极端强降水等灾害性天气，严重影响人民生命和财产安全。多普勒天气雷达是一种用来监测和预报风暴等强对流天气的重要工具，在大气科学和大气探测领域得到广泛的应用。利用雷达探测资料进行风暴识别和追踪的研究，最早可溯源到上世纪60年代，经过半个多世纪的发展，如今已积累了多种行业普遍认可的理论方法，并在实际业务中得以应用。典型的方法有如：TITAN(Thunder Storm Identification,Tracking,Analysis andNowcasting)、TREC(Tracking Radar Echoes by Correlation)和SCIT(Storm CellIdentification and Tracking Algorithm)等。这些方法的主要思想是利用雷达基本反射率因子对风暴主体进行识别，通过滑动示踪窗分析相邻时刻风暴主体之间的相关性，预判风暴移动发展的趋势，再利用线性力学关系函数进行计算，最终达到预测的目的。
理论上，上述基于风暴的三维结构及其生消演变过程的外推预报，由于其算法模型能够反映大气中水汽等粒子运动的物理过程，因此，随着研究和算法模型的不断完善，其外推预报的效果必然也会越来越理想。但遗憾的是，这类方法存在一个共性的问题：风暴主体能否被准确识别，对雷达回波的形态和强度都有一定的要求，特别是对于非对流性的连续性降水系统，这类方法往往无法适用。
发明内容
本发明提出一种基于雷达回波图像扇区分量分析的风暴趋势预测方法。
本发明所采用的技术方案是：
一种基于雷达回波图像扇区分量分析的风暴趋势预测方法，包括以下步骤：
步骤1：设时间T为进行风暴趋势预测的起报时刻，需要预测在未来T+ε时刻的风暴信息；读取某部雷达T时刻及之前最接近T时刻的N次雷达探测资料，设T时刻的雷达探测资料为RD(T)，则任一时刻的雷达资料可表示为RD(T-t+1)，t∈[1,N]；N表示所读取的雷达资料的数量，N取值为3～10的整数；
步骤2：对T时刻的雷达探测资料RD(T)，定义一个数据集存储这些基本反射率数据，其中，表示仰角，λ表示距离库长，λ∈[1,MaxDis]，MaxDis表示雷达探测的最大有效距离，ω表示方位角，ω∈[0,360)；
步骤3：定义一个阈值α，α∈[1,10]的整数，用于将值域[1,MaxDis]按α为单位划分为若干份；定义一个阈值β，β∈[2,20]的整数，用于将值域面上，按α划分的距离库长和β划分的方位角划分成若干个扇区面；
定义表示T时刻的雷达探测资料RD(T)在仰角面上的某一个扇区，该扇区中心位置的距离库长和方位角分别为λ和ω；定义一个阈值δ，δ∈[2,20]的整数，将雷达可探测到的基本反射率划分为若干个值域区间，即Thd＝{[0,δ),[δ,2δ),[2δ,3δ),…,[(m-2)δ,MaxdBZ)}，其中，MaxdBZ表示多普勒天气雷达可能探测到的基本反射率值最大值，m的取值可由MaxdBZ和δ来确定；
步骤4：定义一个包含若干组键值对的数据集其中，key为用于区分和标识Thd的阈值区间，value用于计数，初始值均为0；每个key对应一个value，每一组key和value构成一个键值对；对雷达探测资料RD(T)的任一个仰角面，根据步骤3扇区的划分，依次对每一个扇区做如下操作：
第一，对于雷达探测资料RD(T)中任一扇区查找落入该扇区的所有探测点；对这些探测点逐一分析其基本反射率在域值Thd中的位置，即找到相应的key；再在相应key的键值对中，将value自增1；
第二，将计数value转换为其占所在扇区value累加和的比例，计算方法为：
其中，表示数据集中，键为key时的值，即key所对应的value；表示数据集中，所有key对应的value值的累加和；由此，计算出雷达探测资料RD(T)的各个仰角面每一个扇区中所有探测点的基本反射率在域值Thd中的分布情况；
步骤5：对于雷达探测资料RD(T)中任一扇区依次做如下操作：
第一，读取T时刻前一时刻的雷达探测资料RD(T-1)，生成数据集对雷达探测资料RD(T-1)划定一个与同一仰角面且位置和大小完全相同的扇区，即与步骤3中所述相同，此处的λ和ω分别表示该扇区中心位置的距离库长和方位角；
第二，保持扇区的形状和大小不变，将扇区的中心位置由(λ,ω)平移到(λ-L,ω-M)；此时，扇区变为其中，L表示距离库长的偏移量，L∈[2,20]的整数，M表示方位角的偏移量，L∈[2,20]的整数；按照步骤4的方法，计算该扇区的DSC和DSP；
第三，依次计算扇区 的DSC和DSP；
第四，定义表示扇区的DSP值，其中，x∈[-L,L]，μ∈[-M,M]；计算与的近似指数，计算公式如下：
其中，x∈[-L,L]且以1递增，μ∈[-M,M]且以1°递增；k∈[0,m-1]，m为Thd被划分的区间数量；表示key＝k时，键值对对应的pret值，同理，表示key＝k时，键值对对应的pret值；
由此，计算出雷达资料RD(T-1)中，扇区中心点从至范围内，每个扇区与RD(T)中扇区的基本反射率分布近似指数；从这些指数中找出数值最大的一个，记录下此扇区的位置，即中的x和μ；
进一步地，得到雷达资料RD(T)在扇区处的偏移量为
步骤6：上述步骤得到RD(T)相对于前一时刻RD(T-1)雷达探测资料各仰角面上各个扇区的偏移量，记为采用同样方法，分别计算N次雷达探测资料中，各个时刻相对于前一时刻雷达探测资料各仰角面上各扇区的偏移量；
步骤7：设定N次雷达探测资料划分后的任一扇区为定义Sec[1]、Sec[2]、…、Sec[8]为与为紧相邻的8个扇区；再定义(T-t+1,d,e,f)为扇区中任一探测点，探测点(T-t+1,d,e,f)与Sec[g]中心位置的距离为Dist(T-t+1,d,e,f,g)，其中，g∈[1,2,…,8]；其中，d表示仰角，e表示距离库长，f表示方位角；如图3所示，(T-t+1,d,e,f,g)和Sec[g]都是同一圆形体扫面中的点，因此，Dist(T-t+1,d,e,f,g)可通过几何方法计算得到；
步骤8：定义如下函数，用于预测风暴移动发展的趋势：
其中，表示雷达探测资料RD(T)中任一探测点在未来T+ε时刻的位移矢量；假设T与T-1时刻的时间步长为1，则ε是该时间步长的自然数倍或正小数倍；h(t)是一个权重函数，该函数的构造始终满足单调递减的原则，且表示与雷达资料RD(T-t+1)中仰角λ距离库长ω方位角的风暴所在的扇区紧相邻的某一扇区，则表示扇区的风暴在T-t+1与T-t时刻的单位偏移量；指数p为正数，经验取值范围是[1.0,3.0]；
由此完成雷达探测资料RD(T-t+1)中任一探测点在未来T+ε时刻的位移矢量进一步完成该雷达N组雷达探测资料中所有探测点的移动发展趋势，即完成雷达探测基本反射率所表征的风暴体在T时刻移动趋势预测。
本发明的进一步设计在于：
其中ε表示未来某一时间与时间T的增量，是该时间步长的自然数倍或正小数倍；
步骤3)中，按α划分的距离库长和β划分的方位角划分扇区面时，数据集中每一个探测点都落在且仅落在一个扇区中，对于落在交织线上的探测点，可规定其归属于λ较小或ω较小的一个扇区。
步骤3)中，阈值α将值域[1,MaxDis]划分为若干份，各份为等份。
步骤3)中阈值β将值域[0,360]划分为若干扇区，各份为等份。
该方法还包括步骤9，在步骤8的基础上，利用空间插值算法可计算出风暴在经度、纬度和高度三个维度的移动矢量，由此可做出风暴综合移动发展趋势的预测
本发明相比现有技术具有的优点如下：
1、本发明的预测方法与当前行业广泛采用的TREC、TITAN等风暴识别追踪算法不同，本发明方法不依赖于对雷达强回波中心的识别，特别对于弱回波或没有明显风暴核的回波，也能够较好地分析出风暴过去若干时刻的位移矢量，进而为准确预测奠定了基础。
2、本发明的预测方法与大多数雷达外推预报方法不同，本方法直接对极坐标形式的基本反射率数据进行计算，没有采用空间插值等高空间复杂度的算法进行平面直角坐标系的转换，减少了相关计算的运算量，进而提升了短临天气预报业务的实时性和效能。
附图说明
图1是本发明实施例一的主控流程图；
图2是本发明实施例一示意图之一；
图2展示了在雷达体扫的一个仰角面上，由按α划分的距离库长和由按β划分的方位角交织成的若干个扇区面。
图3是本发明实施例一示意图之二；
图3展示了任一扇区和与其紧相邻的8个扇区Sec[1]、Sec[2]、…、Sec[8]的分布情况；
表1展示了实施例一中键值对中的key与Thd之间的关系。
表2展示了应用实例一中某一扇区的键值对的初始取值情况。
表3展示了应用实例一中某一扇区的键值对的计算结果。
表4展示了应用实例一中某一扇区的键值对的计算结果。
具体实施方式
下面结合附图对本发明的技术方案进一步说明。
实施例一：
如图1所示，本发明的基于雷达回波图像扇区分量分析的风暴趋势预测方法，其具体步骤如下：
步骤1：设时间T为进行风暴趋势预测的起报时刻，需要预测在未来T+ε时刻的风暴信息。读取某部雷达T时刻及之前最接近T时刻的N次雷达探测资料，设T时刻的雷达探测资料为RD(T)，则任一时刻的雷达资料可表示为RD(T-t+1)，t∈[1,N]；N表示所读取的雷达资料的数量，N取值为3～10的整数。其中ε表示未来某一时间与时间T的增量，是该时间步长的自然数倍或正小数倍。
步骤2：步骤1读取的T时刻(即t＝1时)的雷达探测资料RD(T)，该资料RD(T)是由多个仰角，多个距离库长和多个方位角的若干探测点的基本反射率数据构成，这里用“探测点”来描述雷达探测的有效分辨率下的基本单位，探测点的值为基本反射率。定义一个数据集DS用于存储这些基本反射率数据，(对于某种型号DS中包含多个(如9个)仰角、多个距离库长(如460，距离库长是探测点与雷达中心点之间的距离，1个距离库长近似等于1km)、多个方位角(雷达每探测1次，方位角近似转1度，共360度；并且每探测一次，探测一整个距离库长，即460个探测点，)对于任一探测点的基本反射率数据表示为其中，表示仰角，它是雷达体扫时抬升的高度角，通常以平行于地面的仰角为0°，仰角的数量以及每个仰角的取值由雷达体扫所采用的模式决定；λ表示距离库长，λ∈[1,MaxDis]，MaxDis表示雷达探测的最大有效距离，该距离也是由雷达体扫所采用的模式以及雷达的技术参数决定；ω表示方位角，ω∈[0,360)，通常以正北方向的取值为0，按顺时针方向度数逐步增加。
步骤3：定义一个阈值α，α∈[1,10]的整数，用于将值域[1,MaxDis]按α为单位划分为若干等份；定义一个阈值β，β∈[2,20]的整数，用于将值域[0,360]按β为单位划分为若干等份。对雷达探测资料RD(T)的每一仰角面上，由按α划分的距离库长和由按β划分的方位角“交织”成若干个扇区面，扇区面的形状如图2所示。很显然，数据集中每一个探测点都落在且仅落在一个扇区中，对于落在“交织线”上的探测点，可规定其归属于λ较小或ω较小的一个扇区。
为了便于表述，不妨定义表示T时刻的雷达探测资料RD(T)在仰角面上的某一个扇区，该扇区中心位置的距离库长和方位角分别为λ和ω。需要补充说明的是，扇区是由空间上相邻的多个探测点构成的面，中λ和ω分别是扇区面上中心位置的距离库长和方位角，即表征的是一个点，而表征的是一个面。
定义一个阈值δ，δ∈[2,20]的整数，将雷达可探测到的基本反射率划分为若干个值域区间，即Thd＝{[0,δ),[δ,2δ),[2δ,3δ),…,[(m-2)δ,MaxdBZ)}，其中，MaxdBZ是一个经验值，表示多普勒天气雷达可能探测到的基本反射率值最大值，这里MaxdBZ取值为90dBZ。举例来说，假设设定δ＝5，则Thd＝{[0,5),[5,10),…,[85,90)}。m的取值可由MaxdBZ和δ来确定。
步骤4：定义一个包含若干组“键值对”的数据集其中，key为用于区分和标识Thd中不同的阈值区间，如[0,δ)、[δ,2δ)和[2δ,3δ)对应的key取值可依次为1、2和3或者a、b和c；value用于计数，初始值均为0。每个key对应一个value，每一组key和value构成一个“键值对”，DSC中包含了若干组“键值对”。用于标识某一扇区。对雷达探测资料RD(T)的任一个仰角面，根据步骤3扇区的划分，依次对每一个扇区做如下操作：
首先，对于雷达探测资料RD(T)中任一扇区查找落入该扇区的所有探测点。对这些探测点逐一分析其基本反射率在域值Thd中的位置，即找到相应的key；再在相应key的键值对DSC中，将value自增1。举例来说，假设Thd＝{[0,5),[5,10),[10,15),…,[85,90)}，键值对中的key与Thd的关系如表1所示。表1中key的各个值，如0、1、2，…，都只是用于标识和区分不同的“键值对”，0、1、2这样的数字也可以用a、b、c或其他字符来表示。
表1
进一步举例，假设某扇区中包含有100个探测点，其中，第1个探测点的基本反射率值为3，由于3∈[0,5)，则键值对中的value值自增1；类似地，假设第2个探测点的基本反射率值为12，由于12∈[10,15)，则键值对中的value值自增1；以此类推。
接着，将计数value转换为其占所在扇区value累加和的比例，计算方法为：
其中，表示数据集中，“键”为key时的“值”，即key所对应的value。表示数据集中，所有key对应的value值的累加和。很显然，也是一个由若干组“键值对”构成的数据值。
由此，计算出雷达探测资料RD(T)的各个仰角面每一个扇区中所有探测点的基本反射率在域值Thd中的分布情况。
步骤5：对于雷达探测资料RD(T)中任一扇区依次做如下操作：
首先，读取T时刻前一时刻(即t＝2时)的雷达探测资料RD(T-1)，生成数据集对雷达探测资料RD(T-1)划定一个与同一仰角面且位置和大小完全相同的扇区，即与步骤3中所述相同，此处的λ和ω分别表示该扇区中心位置的距离库长和方位角。
然后，保持扇区的形状和大小不变，将扇区的中心位置由(λ,ω)平移到(λ-L,ω-M)；此时，扇区变为其中，L表示距离库长的偏移量，L∈[2,20]的整数，M表示方位角的偏移量，L∈[2,20]的整数。按照步骤4的方法，计算该扇区的DSC和DSP。
接着，依次计算扇区 的DSC和DSP。
为了便于表述，不妨定义表示扇区 的DSP值，其中，x∈[-L,L]，μ∈[-M,M]。
再接着，计算与的近似指数，计算方法为：
其中，x∈[-L,L]且以1递增，μ∈[-M,M]且以1°递增。k∈[0,m-1]，m为Thd被划分的区间数量。表示key＝k时，键值对对应的pret值，同理，表示key＝k时，键值对对应的pret值。
由此，可计算出雷达资料RD(T-1)中，扇区中心点从至范围内，每个扇区与RD(T)中扇区的基本反射率分布近似指数。从这些指数中找出数值最大的一个，记录下此扇区的位置，即中的x和μ。
进一步地，得到雷达资料RD(T)在扇区处的偏移量为
步骤6：上述步骤只计算了RD(T)相对于前一时刻RD(T-1)雷达探测资料各仰角面上各个扇区的偏移量，记为采用相同方法，分别计算N次雷达探测资料中，各个时刻相对于前一时刻雷达探测资料各仰角面上各扇区的偏移量，具体如：由RD(T-1)与RD(T-2)计算得到由RD(T-2)与RD(T-3)计算得到以次类推，可得N-1组扇区的偏移量。为了统一以下步骤与上述公式中的变量符号，上述RD(T)、RD(T-1)、可扩展表示为RD(T-t+1)、RD(T-t)、其中，t的取值为步骤1相同，即t∈[1,N]。
步骤7：设定N次雷达探测资料划分后的任一扇区为为了便于表述，定义Sec[1]、Sec[2]、…、Sec[8]为与紧相邻的8个扇区，如图3所示。
再定义(T-t+1,d,e,f)为某一扇区中任一探测点(如上述步骤所述，一个扇区中包含有多个探测点)。点(T-t+1,d,e,f)与Sec[g]中心位置的距离为Dist(T-t+1,d,e,f,g)，其中，g∈[1,2,…,8]，Dist可通过数学几何方法计算得到。其中，d表示仰角，e表示距离库长，f表示方位角。
步骤8：定义如下函数，用于预测风暴移动发展的趋势：
其中，表示雷达探测资料RD(T)中任一探测点在未来T+ε时刻的位移矢量。这里之所以称之为“矢量”，是因为包括2个维度的信息，即距离库长和方位角。假设T与T-1时刻的时间步长为1，则ε是该时间步长的自然数倍或正小数倍。h(t)是一个权重函数，该函数的构造始终满足单调递减的原则，且 表示与雷达资料RD(T-t+1)中仰角λ距离库长ω方位角的风暴所在的扇区紧相邻的某一扇区，则表示扇区的风暴在T-t+1与T-t时刻的单位偏移量。指数p为正数，一般经验取值范围是[1.0,3.0]，取值越大，相邻扇区sec[g]对的影响越小。
由此完成雷达探测资料RD(T-t+1)中任一探测点在未来T+ε时刻的位移矢量进一步完成该雷达N组雷达探测资料中所有探测点的移动发展趋势，即完成雷达探测基本反射率所表征的风暴体在T时刻移动趋势预测。
步骤9：在步骤8的基础上，利用空间插值算法可计算出风暴在经度、纬度和高度三个维度的移动矢量，由此可做出风暴综合移动发展趋势的预测。
应用实例一：
本实施例选用南京地区的多普勒天气雷达资料，以2018年5月1日16点整作为起报时间，预测未来10分钟、30分钟和60分钟风暴的移动趋势。
步骤1：读取距离起报时刻最近的6个多普勒雷达探测资料，这些资料的真实探测时间分别为当日的07:29、07:35、07:41、07:47、07:53和07:59(均为世界时，与北京时间相差8小时，该部雷达探测的体扫周期约为6分钟，实际体扫周期存在±1分钟的误差，其中07:59时刻的资料最接近北京时间的16点)，为了便于表述，将这6个资料分别命名为RD(T-5)、RD(T-4)、RD(T-3)、RD(T-2)、RD(T-1)和RD(T)。
步骤2：读取T时刻的雷达探测资料RD(T)，该资料是由多个仰角，多个距离库长和多个方位角的基本反射率数据构成，这里用“探测点”来描述雷达探测的有效分辨率下的基本单位，探测点的值为基本反射率。定义一个数据集DS用于存储这些基本反射率数据，对于任一探测点的基本反射率数据表示为其中，是雷达体扫时抬升的高度角，通常以平行于地面的仰角为0°，仰角的数量以及每个仰角的取值由雷达体扫所采用的模式决定；λ∈[1,MaxDis]，根据该雷达的技术参数和体扫方式，此处MaxDis＝460，单位为Km；ω∈[0,360)，通常以正北方向的取值为0，按顺时针方向度数逐步增加。
步骤3：定义阈值α＝10Km，用于将雷达探测的距离库数划分为46等份；定义阈值β＝5°，用于将[0,360°]的仰角划分为72等份。定义表示雷达资料RD(T)在仰角面上的某个扇区，该扇区中心位置的距离库长和方位角分别为λ和ω。定义一个阈值δ＝5，将雷达可探测到的基本反射率划分为若干个值域区间，即Thd＝{[0,5),[5,10),[10,15),…,[85,90)}。
步骤4：定义一个包含若干组“键值对”的数据集其中，key为键，用于区分和标识Thd中不同的阈值范围。value为值，用于计数，初始值均为0。每个key对应一个value，每一组key和value构成一个“键值对”，DSC中包含了若干组“键值对”。DSCSec(T,0.5,37.5,187.5)(key,value)的初始取值情况如表2所示。
表2
对雷达探测资料RD(T)的任一个仰角面，根据步骤3扇区的划分，依次对每一个扇区做如下相同操作。为了更清楚地表述，此处随机选择一个具体的扇区，该扇区中心位置的λ＝37.5Km，ω＝187.5°，即扇区Sec(T,0.5,37.5,187.5)，步骤4的以下内容是针对该扇区的计算过程展开描述：
找到扇区Sec(T,0.5,37.5,187.5)所包含的所有探测点，对这些探测点逐一分析其基本反射率在域值Thd中的位置，即找到相应的key；再在相应key的键值对DSC中，将value自增1。当前扇区的键值对DSCSec(T,0.5,37.5,187.5)(key,value)的计算结果如表3所示。
表3
接着，将计数value转换为其占所在扇区value累加和的比例，计算方法为：
其中，DSCSec(T,0.5,37.5,187.5)(key)表示数据集DSC(T,0.5,37.5,187.5)(key,value)中，“键”为key时的“值”，即key所对应的value。∑DSCSec(T,0.5,37.5,187.5)表示数据集DSC(T,0.5,37.5,187.5)(key,value)中，所有key对应的value值的累加和。
DSPSec(T,0.5,37.5,187.5)的计算结果如表4所示。
表4
由此，计算出扇区Sec(T,0.5,37.5,187.5)中所有探测点的基本反射率在域值Thd中的分布情况。上述雷达探测资料RD(T)每个仰角面上各个扇区的键值对DSP计算方法相同。
步骤5：对于RD(T)中任一扇区依次做如下操作：
首先，读取T时刻前一时刻的雷达探测资料，即RD(T-1)，生成数据集 对雷达探测资料RD(T-1)划定一个与同一仰角面且位置和大小完全相同的扇区，即与步骤3中所述相同，此处的λ和ω分别表示该扇区中心位置的距离库长和方位角。
然后，保持扇区的形状和大小不变，将扇区的中心位置由(λ,ω)平移到(λ-L,ω-M)，此时，扇区变为其中，L表示距离库长的偏移量，L＝10Km，M表示方位角的偏移量，L∈5°。按照步骤4的方法，计算该扇区的DSC和DSP。
接着，依次计算扇区 的DSC和DSP。
为了便于表述，不妨定义表示扇区 的DSP值，其中，x∈[-10,10]，μ∈[-5,5]。
再接着，计算与的近似指数，计算方法为：
其中，x∈[-10,10]且以1递增，μ∈[-5,5]且以1°递增。k∈[0,17]。表示key＝k时，键值对对应的pret值，同理，表示key＝k时，键值对对应的pret值。
由此，可计算出雷达资料RD(T-1)中，扇区中心点从至范围内，每个扇区与RD(T)中扇区的基本反射率分布近似指数。从这些指数中找出数值最大的一个，记录下此扇区的位置，即中的x和μ。
仍以步骤4中Sec(T,0.5,37.5,187.5)的计算为例，在这些近似指数中，数值最大的一个IDX＝0.91，此扇区为Sec(T-1,0.5,34.5,200.5)，通过简单计算，得到雷达资料RD(T)在扇区Sec(T,0.5,37.5,187.5)处的偏移量为OffsetSec(T,0.5,37.5,187.5)＝(3,-3)。
步骤6：上述步骤只计算了RD(T)与RD(T-1)两个时刻的雷达探测资料各仰角面上各个扇区的偏移量，记为采用相同方法，分别计算RD(T-1)与RD(T-2)、RD(T-2)与RD(T-3)、…、RD(T-4)与RD(T-5)各仰角面上各扇区的偏移量，记为
步骤7：为了便于表述，定义任一扇区名称为Sec[0]，则与该扇区紧相邻的8个扇区编号为Sec[1]、Sec[2]、…、Sec[8]，如图3所示。再定义扇区Sec[0]中任一探测点与Sec[g]中心位置的距离为其中，g∈[1,2,…,8]，Dist可通过数学几何方法计算得到。
步骤8：定义如下函数，用于预测风暴移动发展的趋势：
其中，表示雷达资料RD(T)中仰角λ距离库长ω方位角的风暴在未来T+ε时刻的位移矢量。由于本实施例的雷达资料时间间隔为6分钟，要预测未来10、30和60分钟，因此ε取值分别为：1.67、5、10。h(t)是一个权重函数，按照函数单调递减且的原则，这里h(t)采用一组权重系数[0.57,0.20,0.11,0.07,0.05]，t＝[1,2,3,4,5]。表示与雷达资料RD(T-t+1)中仰角λ距离库长ω方位角的风暴所在的扇区 紧相邻的某一扇区，则表示 扇区的风暴在T-t+1与T-t时刻的单位偏移量。指数p为正数，一般经验取值范围是[1.0,3.0]，取值越大，相邻扇区sec[g]对OffsetSec[g]的影响越小。
由此完成雷达探测资料RD(T-t+1)中任一探测点在未来T+ε时刻的位移矢量进一步完成该雷达N组雷达探测资料中所有探测点的移动发展趋势，即完成雷达探测基本反射率所表征的风暴体在T时刻移动趋势预测。
步骤9：在步骤8的基础上，利用空间插值算法可计算出风暴在经度、纬度和高度三个维度的移动矢量，由此可做出风暴综合移动发展趋势的预测。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种基于雷达回波图像扇区分量分析的风暴趋势预测方法，包括以下步骤：
步骤1：设时间T为进行风暴趋势预测的起报时刻，需要预测在未来T+ε时刻的风暴信息；读取某部雷达T时刻及之前最接近T时刻的N次雷达探测资料，设T时刻的雷达探测资料为RD(T)，则任一时刻的雷达资料可表示为RD(T-t+1)，t∈[1,N]；N表示所读取的雷达资料的数量，N取值为3～10的整数；
步骤2：对T时刻的雷达探测资料RD(T)，定义一个数据集存储这些基本反射率数据，其中，表示仰角，λ表示距离库长，λ∈[1,MaxDis]，MaxDis表示雷达探测的最大有效距离，ω表示方位角，ω∈[0,360)；
步骤3：定义一个阈值α，α∈[1,10]的整数，用于将值域[1,MaxDis]按α为单位划分为若干份；定义一个阈值β，β∈[2,20]的整数，用于将值域面上，按α划分的距离库长和β划分的方位角划分成若干个扇区面；
定义表示T时刻的雷达探测资料RD(T)在仰角面上的某一个扇区，该扇区中心位置的距离库长和方位角分别为λ和ω；定义一个阈值δ，δ∈[2,20]的整数，将雷达可探测到的基本反射率划分为若干个值域区间，即Thd＝{[0,δ),[δ,2δ),[2δ,3δ),…,[(m-2)δ,MaxdBZ)}，其中，MaxdBZ表示多普勒天气雷达可能探测到的基本反射率值最大值，m的取值可由MaxdBZ和δ来确定；
步骤4：定义一个包含若干组键值对的数据集其中，key为用于区分和标识Thd的阈值区间，value用于计数，初始值均为0；每个key对应一个value，每一组key和value构成一个键值对；对雷达探测资料RD(T)的任一个仰角面，根据步骤3扇区的划分，依次对每一个扇区做如下操作：
第一，对于雷达探测资料RD(T)中任一扇区查找落入该扇区的所有探测点；对这些探测点逐一分析其基本反射率在域值Thd中的位置，即找到相应的key；再在相应key的键值对中，将value自增1；
第二，将计数value转换为其占所在扇区value累加和的比例，计算方法为：
其中，表示数据集中，键为key时的值，即key所对应的value；表示数据集中，所有key对应的value值的累加和；由此，计算出雷达探测资料RD(T)的各个仰角面每一个扇区中所有探测点的基本反射率在域值Thd中的分布情况；
步骤5：对于雷达探测资料RD(T)中任一扇区依次做如下操作：
第一，读取T时刻前一时刻的雷达探测资料RD(T-1)，生成数据集对雷达探测资料RD(T-1)划定一个与同一仰角面且位置和大小完全相同的扇区，即与步骤3中所述相同，此处的λ和ω分别表示该扇区中心位置的距离库长和方位角；
第二，保持扇区的形状和大小不变，将扇区的中心位置由(λ,ω)平移到此时，扇区变为其中，L表示距离库长的偏移量，L∈[2,20]的整数，M表示方位角的偏移量，L∈[2,20]的整数；按照步骤4的方法，计算该扇区的DSC和DSP；
第三，依次计算扇区 的DSC和DSP；
第四，定义表示扇区的DSP值，其中，x∈[-L,L]，μ∈[-M,M]；计算与的近似指数，计算公式如下：
其中，x∈[-L,L]且以1递增，μ∈[-M,M]且以1°递增；k∈[0,m-1]，m为Thd被划分的区间数量；表示key＝k时，键值对对应的pret值，同理，表示key＝k时，键值对对应的pret值；
由此，计算出雷达资料RD(T-1)中，扇区中心点从至范围内，每个扇区与RD(T)中扇区的基本反射率分布近似指数；从这些指数中找出数值最大的一个，记录下此扇区的位置，即中的x和μ；
进一步地，得到雷达资料RD(T)在扇区处的偏移量为
步骤6：上述步骤得到RD(T)相对于前一时刻RD(T-1)雷达探测资料各仰角面上各个扇区的偏移量，记为采用同样方法，分别计算N次雷达探测资料中，各个时刻相对于前一时刻雷达探测资料各仰角面上各扇区的偏移量；
步骤7：设定N次雷达探测资料划分后的任一扇区为定义Sec[1]、Sec[2]、…、Sec[8]为与为紧相邻的8个扇区；再定义(T-t+1,d,e,f)为扇区中任一探测点，探测点(T-t+1,d,e,f)与Sec[g]中心位置的距离为Dist(T-t+1,d,e,f,g)，其中，g∈[1,2,…,8]；其中，d表示仰角，e表示距离库长，f表示方位角；求得Dist(T-t+1,d,e,f,g)；
步骤8：定义如下函数，用于预测风暴移动发展的趋势：
其中，表示雷达探测资料RD(T)中任一探测点在未来T+ε时刻的位移矢量；假设T与T-1时刻的时间步长为1，则ε是该时间步长的自然数倍或正小数倍；h(t)是一个权重函数，该函数的构造始终满足单调递减的原则，且 表示与雷达资料RD(T-t+1)中仰角λ距离库长ω方位角的风暴所在的扇区紧相邻的某一扇区，则表示扇区的风暴在T-t+1与T-t时刻的单位偏移量；指数p为正数，经验取值范围是[1.0,3.0]；
由此完成雷达探测资料RD(T-t+1)中任一探测点在未来T+ε时刻的位移矢量进一步完成该雷达N组雷达探测资料中所有探测点的移动发展趋势，即完成雷达探测基本反射率所表征的风暴体在T时刻移动趋势预测。
2.根据权利要求1所述风暴趋势预测方法，其其中ε表示未来某一时间与时间T的增量，是该时间步长的自然数倍或正小数倍。
3.根据权利要求1所述风暴趋势预测方法，其中步骤3)中，按α划分的距离库长和β划分的方位角划分扇区面时，数据集中每一个探测点都落在且仅落在一个扇区中，对于落在交织线上的探测点，可规定其归属于λ较小或ω较小的一个扇区。
4.根据权利要求3所述风暴趋势预测方法，其中步骤3)中，阈值α将值域[1,MaxDis]划分为若干份，各份为等份。
5.根据权利要求3所述风暴趋势预测方法，其中步骤3)中阈值β将值域[0,360]划分为若干扇区，各份为等份。
6.根据权利要求1所述风暴趋势预测方法，还包括步骤9，在步骤8的基础上，利用空间插值算法可计算出风暴在经度、纬度和高度三个维度的移动矢量，由此可做出风暴综合移动发展趋势的预测。
CN201810826357.3A 2018-07-25 2018-07-25 基于雷达回波图像扇区分量分析的风暴趋势预测方法 Active CN109100722B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810826357.3A CN109100722B (zh) 2018-07-25 2018-07-25 基于雷达回波图像扇区分量分析的风暴趋势预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810826357.3A CN109100722B (zh) 2018-07-25 2018-07-25 基于雷达回波图像扇区分量分析的风暴趋势预测方法

Publications (2)

Publication Number Publication Date
CN109100722A true CN109100722A (zh) 2018-12-28
CN109100722B CN109100722B (zh) 2022-08-12

Family

ID=64847429

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810826357.3A Active CN109100722B (zh) 2018-07-25 2018-07-25 基于雷达回波图像扇区分量分析的风暴趋势预测方法

Country Status (1)

Country Link
CN (1) CN109100722B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109917394A (zh) * 2019-03-13 2019-06-21 南京信息工程大学 一种基于天气雷达的短临智能外推方法
CN110515081A (zh) * 2019-06-26 2019-11-29 南京信息工程大学 一种雷达回波零度层亮带智能识别预警方法
CN110574614A (zh) * 2019-08-07 2019-12-17 大连市人工影响天气办公室 一种人工防御冰雹科学应用安全射界的有效作业方法及系统
CN110673146A (zh) * 2019-10-12 2020-01-10 上海眼控科技股份有限公司 气象预测图检测方法、装置、计算机设备和可读存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6654689B1 (en) * 2000-11-06 2003-11-25 Weather Central, Inc. System and method for providing personalized storm warnings
CN102721987A (zh) * 2012-06-12 2012-10-10 中国海洋大学 多普勒雷达遥感强风暴的预警方法
CN103337133A (zh) * 2013-06-14 2013-10-02 广东电网公司中山供电局 基于识别预报技术的电网雷暴灾害预警系统及方法
CN106570598A (zh) * 2016-11-15 2017-04-19 中国海洋大学 一种基于扩展卡尔曼滤波的风暴潮灾害预警方法
CN107193060A (zh) * 2017-06-20 2017-09-22 厦门大学 一种多路径台风风暴潮快速预测方法及系统
CN107526083A (zh) * 2017-10-18 2017-12-29 国网新疆电力公司电力科学研究院 一种基于天气雷达数据的强对流风力等级预测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6654689B1 (en) * 2000-11-06 2003-11-25 Weather Central, Inc. System and method for providing personalized storm warnings
CN102721987A (zh) * 2012-06-12 2012-10-10 中国海洋大学 多普勒雷达遥感强风暴的预警方法
CN103337133A (zh) * 2013-06-14 2013-10-02 广东电网公司中山供电局 基于识别预报技术的电网雷暴灾害预警系统及方法
CN106570598A (zh) * 2016-11-15 2017-04-19 中国海洋大学 一种基于扩展卡尔曼滤波的风暴潮灾害预警方法
CN107193060A (zh) * 2017-06-20 2017-09-22 厦门大学 一种多路径台风风暴潮快速预测方法及系统
CN107526083A (zh) * 2017-10-18 2017-12-29 国网新疆电力公司电力科学研究院 一种基于天气雷达数据的强对流风力等级预测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
热苏里 阿布拉等: "基于多普勒天气雷达的冰雹云早期识别与预警方法研究", 《冰川冻土》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109917394A (zh) * 2019-03-13 2019-06-21 南京信息工程大学 一种基于天气雷达的短临智能外推方法
CN109917394B (zh) * 2019-03-13 2022-12-23 南京信息工程大学 一种基于天气雷达的短临智能外推方法
CN110515081A (zh) * 2019-06-26 2019-11-29 南京信息工程大学 一种雷达回波零度层亮带智能识别预警方法
CN110574614A (zh) * 2019-08-07 2019-12-17 大连市人工影响天气办公室 一种人工防御冰雹科学应用安全射界的有效作业方法及系统
CN110574614B (zh) * 2019-08-07 2021-05-18 大连市人工影响天气办公室 一种人工防雹安全射界图科学有效应用的技术方法及系统
CN110673146A (zh) * 2019-10-12 2020-01-10 上海眼控科技股份有限公司 气象预测图检测方法、装置、计算机设备和可读存储介质

Also Published As

Publication number Publication date
CN109100722B (zh) 2022-08-12

Similar Documents

Publication Publication Date Title
CN109100722A (zh) 基于雷达回波图像扇区分量分析的风暴趋势预测方法
Krajewski et al. Towards better utilization of NEXRAD data in hydrology: An overview of Hydro-NEXRAD
Smith et al. Analyses of a long‐term, high‐resolution radar rainfall data set for the Baltimore metropolitan region
Fausto et al. A snow density dataset for improving surface boundary conditions in Greenland ice sheet firn modeling
Nanding et al. Comparison of different radar-raingauge rainfall merging techniques
Biggs et al. A comparison of gauge and radar precipitation data for simulating an extreme hydrological event in the Severn Uplands, UK
CN111401602B (zh) 基于神经网络的卫星以及地面降水测量值同化方法
CN109523066B (zh) 一种基于克里金插值的pm2.5新增移动站点选址方法
CN105891833A (zh) 基于多普勒雷达信息识别暖云降水率的方法
McRoberts et al. Detecting beam blockage in radar-based precipitation estimates
CN108876172B (zh) 一种基于改进型modis植被供水指数的地表土壤含水量评估方法
McCabe et al. Glacier variability in the conterminous United States during the twentieth century
CN111044110A (zh) 一种基于相似度分析的气体超声波流量计信号处理方法
Damant et al. Errors in the Thiessen technique for estimating areal rain amounts using weather radar data
Tian et al. Accuracy assessment and error cause analysis of GPM (V06) in Xiangjiang river catchment
Xia et al. A model to interpolate monthly mean climatological data at Bavarian forest climate stations
Bakış et al. Analysis and comparison of spatial rainfall distribution applying different interpolation methods in Porsuk river basin, Turkey
Wang et al. Ice thickness, volume and subglacial topography of Urumqi Glacier No. 1, Tianshan Mountains, central Asia, by ground penetrating radar survey
CN113390471B (zh) 一种基于gnss反射信号的河流流量估算方法
Shadeed et al. Comparative analysis of interpolation methods for rainfall mapping in the Faria catchment, Palestine
Bebbington et al. Modelling of weather radar echoes from anomalous propagation using a hybrid parabolic equation method and NWP model data
Berndtsson et al. Some Eulerian and Lagrangian statistical properties of rainfall at small space-time scales
Borup et al. Comparing the impact of time displaced and biased precipitation estimates for online updated urban runoff models
CN107807388A (zh) 一种基于多普勒效应的地震断层滑动速度计算方法
Mazari et al. Validation of the NEXRAD DSP product with a dense rain gauge network

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