CN113296074B - 一种基于气象雷达多层cappi的光流法外推方法 - Google Patents

一种基于气象雷达多层cappi的光流法外推方法 Download PDF

Info

Publication number
CN113296074B
CN113296074B CN202110853509.0A CN202110853509A CN113296074B CN 113296074 B CN113296074 B CN 113296074B CN 202110853509 A CN202110853509 A CN 202110853509A CN 113296074 B CN113296074 B CN 113296074B
Authority
CN
China
Prior art keywords
data
layer
height
cappi
extrapolation
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
CN202110853509.0A
Other languages
English (en)
Other versions
CN113296074A (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.)
Chengdu Yuanwang Detection Technology Co ltd
Original Assignee
Chengdu Yuanwang Detection 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 Chengdu Yuanwang Detection Technology Co ltd filed Critical Chengdu Yuanwang Detection Technology Co ltd
Priority to CN202110853509.0A priority Critical patent/CN113296074B/zh
Publication of CN113296074A publication Critical patent/CN113296074A/zh
Application granted granted Critical
Publication of CN113296074B publication Critical patent/CN113296074B/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
    • 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/411—Identification of targets based on measurements of radar reflectivity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Operations Research (AREA)
  • Computing Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种基于气象雷达多层CAPPI的光流法外推方法，所述外推方法包括：S1、解析雷达体积扫描文件中的反射率数据；S2、采用多层CAPPI算法计算出各个高度层的CAPPI数据；S3、采用DARTS算法基于各个高度层CAPPI数据计算出光流场；S4、采用半拉格朗日算法计算出各个高度层的光流场外推数据；S5、将不同高度的外推数据进行融合处理。本发明的优点在于：考虑到空气受到气压梯度力而形成的梯度风，不同高度层的风运动速度和方向不同，从而影响雷达反射率外推结果，通过计算不同高度层的反射率，将每层反射率进行外推计算，最终融合每层反射率作为最终的反射率外推结果。

Description

一种基于气象雷达多层CAPPI的光流法外推方法
技术领域
本发明涉及气象雷达技术领域，尤其涉及一种基于气象雷达多层CAPPI的光流法外推方法。
背景技术
气象雷达的外推应用中，常常使用气象雷达的反射率数据和光流法生成预报区域的外推流场；由于空气受到气压梯度力而形成的梯度风，不同高度层的风运动速度和方向有所不同，从而影响气象雷达的反射率外推结果；因此，如何通过计算不同高度层的反射率，进而得到最终融合每层反射率作为最终反射率外推场结果的，是现阶段需要解决的问题。
发明内容
本发明的目的在于克服现有技术的缺点，提供了一种基于气象雷达多层CAPPI的光流法外推方法，解决了现有技术存在的不足。
本发明的目的通过以下技术方案来实现：一种基于气象雷达多层CAPPI的光流法外推方法，所述外推方法包括：
S1、解析雷达体积扫描文件中的反射率数据；
S2、采用多层CAPPI算法计算出各个高度层的CAPPI数据；
S3、采用DARTS算法基于各个高度层CAPPI数据计算出光流场；
S4、采用半拉格朗日算法计算出各个高度层的光流场外推数据；
S5、将不同高度的外推数据进行融合处理。
所述采用多层CAPPI算法计算出各个高度层的CAPPI数据包括：
S21、以等高平面上的距离为单位点，从某方位角上第一个距离库开始计算其对应的仰角，得到该仰角对应的上下体扫描仰角；
S22、判断该仰角与上下体扫描仰角的关系，用线性插值方法得到该高度上的数据；
S23、对每个方位角径向上的每一点都进行步骤S21和步骤S22的操作，直到平面全部遍历完成。
所述采用DARTS算法基于各个高度层CAPPI数据计算出光流场包括：
Figure 171179DEST_PATH_IMAGE002
来表示雷达反射率场的时间序列的降水模式和演变，其中，F(x,y,t)表示雷达图像序列，U(x,y)表示东西方向速度场，V(x,y)为南北方向速度场，S(x,y,t)表示加性演化机制的序列，其中x是图像的横坐标，y是图像的纵坐标，t是外推时间；
将上述公式进行离散化得到离散公式
Figure 405238DEST_PATH_IMAGE004
；
选定
Figure 210383DEST_PATH_IMAGE006
的最大谐波数的整数集为N={Nx，Ny，Nt}、
Figure 732631DEST_PATH_IMAGE008
和
Figure 493783DEST_PATH_IMAGE010
的最大谐波数的整数集为M={Mx，My}以及
Figure 708864DEST_PATH_IMAGE012
的最大谐波数的整数集为L={Lx，Ly，Lt}，其中Nx，Ny，Nt，Mx，My，Lx，Ly，Lt分别表示对应数据集的谐波数，其中DFT表示离散傅里叶变换；
通过快速傅里叶变换对离散公式求解偏导数并设置索引变量和表
Figure 922676DEST_PATH_IMAGE014
，
Figure 514194DEST_PATH_IMAGE016
，
Figure 146164DEST_PATH_IMAGE018
，得到块矩阵Hm×n =[A,B,C]，其中，
Figure 656780DEST_PATH_IMAGE019
、
Figure 880258DEST_PATH_IMAGE020
、
Figure 744309DEST_PATH_IMAGE021
、
Figure 480052DEST_PATH_IMAGE022
、
Figure 302515DEST_PATH_IMAGE023
、
Figure 241652DEST_PATH_IMAGE024
、
Figure 424241DEST_PATH_IMAGE025
和
Figure 296382DEST_PATH_IMAGE026
皆表示索引变量，
Figure 227429DEST_PATH_IMAGE028
；
Figure 965446DEST_PATH_IMAGE030
；
Figure 374562DEST_PATH_IMAGE032
；
设置向量元素集
Figure 632368DEST_PATH_IMAGE034
、
Figure 986513DEST_PATH_IMAGE036
、
Figure 696981DEST_PATH_IMAGE038
、
Figure 159055DEST_PATH_IMAGE040
，将离散公式变换为
Figure 536946DEST_PATH_IMAGE042
，进而得到矩阵点积形式
Figure 544217DEST_PATH_IMAGE043
，其中，
Figure 256827DEST_PATH_IMAGE045
为伪逆矩阵；
采用最小二乘解得到反向得到伪逆矩阵的分割子矩阵的表示：
Figure 804483DEST_PATH_IMAGE047
，进而通过分割子矩阵最终得到速度场水平方向的东西方向速度场U和垂直方向的南北方向速度场V。
采用半拉格朗日算法计算出各个高度层的光流场随时间变化的结果，外推出之后数小时或数分钟的雷达回波数据包括：
采用二维守恒方程的微分形式
Figure 240143DEST_PATH_IMAGE049
表示平流，对其进行变换得到
Figure 464320DEST_PATH_IMAGE051
；
设置局部变化率∂Ψ/∂t为零得到
Figure 352642DEST_PATH_IMAGE053
，根据平流跟随降水包的运动得到
Figure 969568DEST_PATH_IMAGE055
，设置源汇项
Figure 246353DEST_PATH_IMAGE057
来表示降水的增长和消散，进而得到
Figure 923322DEST_PATH_IMAGE059
，表示在
Figure DEST_PATH_IMAGE061
时刻和位置x时预测到的变化率，其中，Ψ为观测到的降水场，t 0为外推的开始时间；
将平流分成N个时间步Δt来表示间隔时间 τ ，即τ =NΔt，对每个时间步的α 通过迭代的方法确定位移矢量
Figure DEST_PATH_IMAGE063
；
从α =0开始，最终位移矢量是单个时间步的N个矢量的矢量和，因此，在半拉格朗日格式中通过假设速度平稳性，即
Figure 954732DEST_PATH_IMAGE065
，沿着向前的时间线或向后的时间线来确定降水包的轨迹；
根据平流层光流场DARTS的计算方法获得的回波运动场在公式
Figure DEST_PATH_IMAGE066
中进行迭代收敛，循环多次后中止迭代得到相应数据，其中，t表示外推时间，u是x方向上的速度场，v表示y方向的速度场；u，v都是上面的光流法得到的结果，
Figure 31141DEST_PATH_IMAGE068
是网格点
Figure 425082DEST_PATH_IMAGE070
处的回波运动。
所述将不同高度的外推数据进行融合处理步骤如下：
根据得到的每个层的雷达回波预报数据，每个层的数据用二维网格的图像表示，每个格点的值为-1到127的反射率值，最终的数据张量的形状为(c, height, width)，其中c表示层数、height表示图像的高度、width表示图像的宽度；
每个外推时刻t 0+Δt的结果是通过在N的方向上取极大值，得到不同时刻的网格点面数据；最终的数据张量的形状为(1, height, width)，就是只剩下一层图像方便气象员可视化雷达图；
总的外推时刻为t 0+NΔt，根据迭代过程就会得到N个形状为(1, height, width)的雷达外推数据。
本发明具有以下优点：一种基于气象雷达多层CAPPI的光流法外推方法，考虑到空气受到气压梯度力而形成的梯度风，不同高度层的风运动速度和方向不同，从而影响雷达反射率外推结果，通过计算不同高度层的反射率，将每层反射率进行外推计算，最终融合每层反射率作为最终的反射率外推结果。
附图说明
图1 为本发明的流程示意图；
图2 为CAPPI算法原理示意图；
图3 为半拉格朗日向量示意图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。因此，以下结合附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的保护范围，而是仅仅表示本申请的选定实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。下面结合附图对本发明做进一步的描述。
如图1所示，本发明涉及一种基于气象雷达多层CAPPI的光流法外推方法，气象雷达的体积扫描反射率作为外推原始数据，通过算法生成不同高度的等高平面位置显示器（Constant Altitude Plan Position Indicato，CAPPI）数据，The Dynamic and AdaptiveRadar Tracking of Storms（DARTS）方法计算各个高度层的CAPPI数据在空间上的运动矢量估计能力，结合半拉格朗日方法得出各个高度层的外推结果，最终将各个高度层数据融合得出外推结果，图中，H为不同海拔的高度，t为外推时间；具体包括以下内容：
S1、解析雷达体积扫描文件中的反射率数据；
S2、采用多层CAPPI算法计算出各个高度层的CAPPI数据；
具体为，多层CAPPI产品是雷达体积扫描中的反射率数据进行计算。多层由不同高度的单层CAPPI组成，高度分为：500，1000，2000，3000，3500，4000，4500，5000，6000，7000米共10层高度。以下介绍单层CAPPI的计算方法，CAPPI是一个投影面，该投影面上的数据是来自于切割平面与各个扫描层相交处的数据。
按照设定的高度,应用测高公式，选取临近该高度平面上的上下两个仰角相应雷达测距上的数据,然后用内插方法得到该高度上的数据。为提高数据精度,常采用双线性插值。
S21、以等高平面上的距离库为单位点，从某方位角上第一个距离库开始计算其对应的仰角，找到该点上下两个体扫仰角；
S22、判断该仰角与上下体扫仰角的关系，用线性插值方法得到该高度上的数据；
S23、对每个方位角径向上的每一点都进行步骤S21和步骤S22的计算，知道平面全部遍历完成。
如图2所示，图中，地球平均半径为R，切割高度为h，天线海拔高度为H；举例如下：
第一步：根据测高公式计算切割面上某径向上点A的仰角，得到测高公式：
Figure 741794DEST_PATH_IMAGE072
，其中，α即点A所在的仰角，SlatRan为点A的径向距离；
第二步：判断点A的仰角α与体扫各层仰角的关系；
第三步：根据仰角关系进行内插取值；如果α刚好等于某仰角值，则直接取相应仰角PPI上的值为该点的CAPPI值；如果α小于最低体扫仰角值，则取最低仰角PPI上的值为该点的CAPPI值；如果α大于最高体扫仰角值，则视该点无回波值；如果α在两个体扫仰角之间，则进行线性插值；
第四步：对每个径向上每个距离库都进行第一到第三步的计算，如果α在两个体扫仰角之间，则对点A进行垂直方向上线性插值。首先根据点A、B、E所对应的水平距离相同，计算出B、E在对应仰角层对应的径向距离和到水平面的高度，再判断点B、E对应的回波值是否有效，最后根据高度权重进行插值。
S3、采用光流场DARTS算法基于各个高度层CAPPI数据计算出光流场；
The Dynamic and Adaptive Radar Tracking of Storms（DARTS）基于CAPPI数据在单层空间上的运动矢量计算。
通过描述由下列公式给出的雷达反射率场的时间序列表示的降水模式的通量和演变来实现的，
其中，F(x, y, t)是雷达图像序列（2张t=2），U(x, y)是东西方向速度场，V(x, y)为南北方向速度场，S(x, y, t) 可以解释为加性演化机制的序列，如强度的增长和衰减；其中x是图像的横坐标，y是图像的纵坐标，t是外推时间；
Figure 401445DEST_PATH_IMAGE002
离散形式公式为：
Figure 550055DEST_PATH_IMAGE073
；
此公式可以转换为线性形式，允许有效的程序实现，其中，“+”表示左边的像素，“-”表示右边的像素。
设参数：
N={Nx,Ny,Nt}，M={Mx,My}，L={Lx,Ly,Lt}为
Figure 611551DEST_PATH_IMAGE074
、
Figure 568006DEST_PATH_IMAGE075
和
Figure 964221DEST_PATH_IMAGE076
、
Figure 398745DEST_PATH_IMAGE077
选定的最大谐波数的整数集。
本算法选取了N={20，20，1}，M={2，2}，L={0,0,0}；
算法中使用快速傅里叶变换来求解偏导数，Nx默认值20，Ny=20， Nt=1，Mx=2，My=2，L={ 0，0，0}；
并且定义以下索引变量和表：
Figure DEST_PATH_IMAGE078
，
Figure 236120DEST_PATH_IMAGE079
，
Figure DEST_PATH_IMAGE080
其中
Figure DEST_PATH_IMAGE082
，
Figure DEST_PATH_IMAGE084
，
Figure DEST_PATH_IMAGE086
，
Figure DEST_PATH_IMAGE088
，
Figure DEST_PATH_IMAGE090
，
Figure DEST_PATH_IMAGE092
，
Figure DEST_PATH_IMAGE094
，
Figure DEST_PATH_IMAGE096
；
这里只考虑了DFT的正整数部分，DFT是快速傅里叶变换。
Figure DEST_PATH_IMAGE098
表示矩阵相乘操作即笛卡儿积。由离散形式公式和上面定义的索引变量得到下面的块矩阵：
Figure 142238DEST_PATH_IMAGE099
；
Figure 25750DEST_PATH_IMAGE100
；
Figure 263964DEST_PATH_IMAGE032
；
Hm×n =[A,B,C] (A,B,C就是上面的三个块矩阵拼接成H矩阵)；
其中
Figure 34474DEST_PATH_IMAGE102
满足条件：
Figure 584929DEST_PATH_IMAGE104
；同时我们定义以下向量元素集：
Figure DEST_PATH_IMAGE105
、
Figure 627840DEST_PATH_IMAGE106
、
Figure DEST_PATH_IMAGE107
、
Figure 138587DEST_PATH_IMAGE108
离散公式就可以写成如下的形式：
Figure 278450DEST_PATH_IMAGE110
进而得到线性系统的观测和响应向量以及设计矩阵（伪逆矩阵的结果的矩阵）：
Figure 544346DEST_PATH_IMAGE112
，
Figure DEST_PATH_IMAGE113
这里
Figure DEST_PATH_IMAGE115
，
Figure DEST_PATH_IMAGE117
，
Figure DEST_PATH_IMAGE119
，
Figure DEST_PATH_IMAGE121
；
其中Re是取复数的实数部分，Im是取复数的虚数部分。
Figure DEST_PATH_IMAGE123
，
Figure DEST_PATH_IMAGE125
，
Figure DEST_PATH_IMAGE127
，
Figure DEST_PATH_IMAGE129
，
Figure DEST_PATH_IMAGE131
，
Figure DEST_PATH_IMAGE133
，
Figure DEST_PATH_IMAGE135
，
Figure 310439DEST_PATH_IMAGE137
，
Figure 156035DEST_PATH_IMAGE139
，
Figure DEST_PATH_IMAGE141
；
进而得到如下矩阵点积形式：
Figure DEST_PATH_IMAGE143
矩阵伪逆
Figure 559859DEST_PATH_IMAGE144
，
Figure DEST_PATH_IMAGE145
由最小二乘法得到解，
Figure 183608DEST_PATH_IMAGE146
用于估计运动矢量场的DFT快速傅里叶变换的系数和演化项S。
应用最小二乘解反向求得：
Figure DEST_PATH_IMAGE147
矩阵最终得到想要的速度场水平方向的U和垂直方向的V；
Figure 935532DEST_PATH_IMAGE149
，
Figure 850398DEST_PATH_IMAGE151
，
Figure DEST_PATH_IMAGE153
，
Figure 840220DEST_PATH_IMAGE155
；
通过反向傅里叶法的实数部分得到雷达图在横向、纵向上的速度场U、V结果的形状跟观测数据的形状一样都是(1, height, width)，物理意义就是对应像素的速度。
S4、采用半拉格朗日算法计算出各个高度层的光流场外推数据；
采用半拉格朗日算法计算出各个高度层的光流场随时间变化的结果，外推出之后数小时或数分钟的雷达回波数据。
因为可以将相同高度的大气层上的雷达回波想象成水平运动，所以CAPPI的每一层就相当于相同高度上的水平运动观测到的雷达回波数据，每层单独使用半拉格朗日算法就可以算出每个高度上的雷达回波情况，最终就是输入的有多少高度层，就会产生多少层的外推结果了。
每个层算法相同，应用简单的平流外推semilagrangian半拉格朗日算法生成一个预报步骤如下：
用符号Ψ 表示平流，可以用二维守恒方程的微分形式来表示Ψ，具体内容如下：
Figure DEST_PATH_IMAGE156
忽略压缩性项Ψ(∂u/∂x+∂v/∂y）我们得到：
Figure 182340DEST_PATH_IMAGE157
局部变化率∂Ψ/∂t设置为零。
Figure DEST_PATH_IMAGE158
平流Ψ 跟随降水包的运动（拉格朗日持久性）；
Figure 412771DEST_PATH_IMAGE159
这里的α 是位移向量表示（u，v）的场。我们可以计算解相关时间，从而得到另一种可预测性度量，即通过拉格朗日持久性得到的可预测性。作为拉格朗日持久性的一种改进，我们首先引入了源汇项
Figure 583858DEST_PATH_IMAGE160
，它表示降水的增长和消散，然后替换上面的等式右边的Ψ，得到：
Figure DEST_PATH_IMAGE161
在这里，
Figure DEST_PATH_IMAGE163
是从Ψ 通过任何类型的变换，如光谱滤波或时间和/或空间上的平均，或者是傅里叶变换。公式
Figure DEST_PATH_IMAGE164
是本发明所用拉格朗日预报程序的一般方程。
其中的位移向量α的常向量形式为：
Figure DEST_PATH_IMAGE166
。
Figure DEST_PATH_IMAGE168
是网格点P处的回波运动。对于每个网格点，常量向量方法使用一个常量平移向量，并且不允许旋转。为了克服这个缺点，本发明使用半拉格朗日（semi-Lagrangian）格式：平流被分成N个时间步为Δt来表示 τ ，即τ =NΔt；
对于每个时间步的α 通过迭代的方法确定：
Figure DEST_PATH_IMAGE169
（半拉格朗日形式）。
从α = 0开始，最终位移矢量是单个时间步的N个矢量的矢量和。因此，在半拉格朗日形式中通过假设速度平稳性，即
Figure DEST_PATH_IMAGE170
，沿着向前的时间线或向后的时间线来确定降水包的轨迹。
由上面的平流层光流场DARTS的计算方法获得的回波运动场在公式
Figure 149225DEST_PATH_IMAGE066
的迭代过程中收敛速度很快，在2或3次循环后可以中止迭代。
如图3所示，显示了半拉格朗日向前 semi-Lagrangian forward（slf）、半拉格朗日向后 semi-Lagrangian backward（slb）、常量向量向前 constant-vector forward（cvf）和常量向量向后constant-vector backward（cvb）的向量。
如果旋转不可忽略，如天气尺度，半拉格朗日方法显然是首选。定性地表明，即使在强自转的情况下，该方案也能保持降水系统的形状。
S5、融合不同高度层的外推数据；
将不同高度的外推数据融合处理，处理方法是：遍历每个网格点所对应不同高度反射率外推数据取最大值，最终得出不同时刻的网格点面数据。
根据得到的每个层的雷达回波预报数据，每个层的数据用二维网格的图像表示，每个格点的值为-1到127的反射率值，最终的数据张量的形状为(c, height, width)，其中c表示层数、height表示图像的高度、width表示图像的宽度；
每个外推时刻t 0+Δt的结果是通过在N的方向上取极大值，得到不同时刻的网格点面数据；最终的数据张量的形状为(1, height, width)，就是只剩下一层图像方便气象员可视化雷达图；
总的外推时刻为t 0+NΔt，根据迭代过程就会得到N个形状为(1, height, width)的雷达外推数据。
以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (3)

1.一种基于气象雷达多层CAPPI的光流法外推方法，其特征在于：所述外推方法包括：
S1、解析雷达体积扫描文件中的反射率数据；
S2、采用多层CAPPI算法计算出各个高度层的CAPPI数据；
S3、采用DARTS算法基于各个高度层CAPPI数据计算出光流场；
S4、采用半拉格朗日算法计算出各个高度层的光流场外推数据；
S5、将不同高度的外推数据进行融合处理；
所述采用多层CAPPI算法计算出各个高度层的CAPPI数据包括：
S21、以等高平面上的距离为单位点，从某方位角上第一个距离库开始计算其对应的仰角，得到该仰角对应的上下体扫描仰角；
S22、判断该仰角与上下体扫描仰角的关系，用线性插值方法得到该高度上的数据；
S23、对每个方位角径向上的每一点都进行步骤S21和步骤S22的操作，直到平面全部遍历完成；
所述采用DARTS算法基于各个高度层CAPPI数据计算出光流场包括：
通过公式
Figure 93792DEST_PATH_IMAGE001
来表示雷达反射率场的时间序列的降水模式和演变，其中，F(x,y,t)表示雷达图像序列，U(x,y)表示东西方向速度场，V(x,y)为南北方向速度场，S(x,y,t)表示加性演化机制的序列，其中x是图像的横坐标，y是图像的纵坐标，t是外推时间；
将上述公式进行离散化得到离散公式
Figure 884025DEST_PATH_IMAGE002
；
选定
Figure 952475DEST_PATH_IMAGE003
的最大谐波数的整数集为N={Nx，Ny，Nt}、
Figure 153649DEST_PATH_IMAGE004
和
Figure 550127DEST_PATH_IMAGE005
的最大谐波数的整数集为M={Mx，My}以及
Figure 127738DEST_PATH_IMAGE006
的最大谐波数的整数集为L={Lx，Ly，Lt}，其中Nx，Ny，Nt，Mx，My，Lx，Ly，Lt分别表示对应数据集的谐波数，其中DFT表示离散傅里叶变换；
通过快速傅里叶变换对离散公式求解偏导数并设置索引变量和表
Figure 188711DEST_PATH_IMAGE007
，
Figure 826366DEST_PATH_IMAGE008
，
Figure 100352DEST_PATH_IMAGE009
，其中，
Figure 232387DEST_PATH_IMAGE010
、
Figure 400063DEST_PATH_IMAGE011
、
Figure 959352DEST_PATH_IMAGE012
、
Figure 455055DEST_PATH_IMAGE013
、
Figure 640049DEST_PATH_IMAGE014
、
Figure 412964DEST_PATH_IMAGE015
、
Figure 533367DEST_PATH_IMAGE016
和
Figure 906579DEST_PATH_IMAGE017
皆表示索引变量，得到块矩阵
Figure 129486DEST_PATH_IMAGE018
；
Figure 6175DEST_PATH_IMAGE019
；
Figure 563059DEST_PATH_IMAGE020
，其中，Hm×n =[A,B,C] 是A、B、C三个块矩阵拼接成的H矩阵；
设置向量元素集
Figure 908721DEST_PATH_IMAGE021
、
Figure 435517DEST_PATH_IMAGE022
、
Figure 917445DEST_PATH_IMAGE023
、
Figure 379650DEST_PATH_IMAGE024
，将离散公式变换为
Figure 461876DEST_PATH_IMAGE025
，进而得到矩阵点积形式
Figure 277516DEST_PATH_IMAGE026
，其中，
Figure 128797DEST_PATH_IMAGE027
为伪逆矩阵；
采用最小二乘解得到反向得到伪逆矩阵的分割子矩阵的表示：
Figure 761904DEST_PATH_IMAGE028
，进而通过分割子矩阵最终得到速度场水平方向的东西方向速度场U和垂直方向的南北方向速度场V，其中
Figure 813649DEST_PATH_IMAGE029
，
Figure 682248DEST_PATH_IMAGE030
，
Figure 404348DEST_PATH_IMAGE031
，
Figure 942776DEST_PATH_IMAGE032
，
Figure 999594DEST_PATH_IMAGE033
，
Figure 891458DEST_PATH_IMAGE034
，
Figure 717331DEST_PATH_IMAGE035
，
Figure 302028DEST_PATH_IMAGE036
，
Figure 721508DEST_PATH_IMAGE037
，
Figure 931909DEST_PATH_IMAGE038
，Re是取复数的实数部分，Im是取复数的虚数部分。
2.根据权利要求1所述的一种基于气象雷达多层CAPPI的光流法外推方法，其特征在于：采用半拉格朗日算法计算出各个高度层的光流场随时间变化的结果，外推出之后数小时或数分钟的雷达回波数据包括：
采用二维守恒方程的微分形式
Figure 365951DEST_PATH_IMAGE039
表示平流，对其进行变换得到
Figure 636396DEST_PATH_IMAGE040
；
设置局部变化率∂Ψ/∂t为零得到
Figure 277592DEST_PATH_IMAGE041
，根据平流跟随降水包的运动得到
Figure 776838DEST_PATH_IMAGE042
，设置源汇项
Figure 577304DEST_PATH_IMAGE043
来表示降水的增长和消散，进而得到
Figure 503803DEST_PATH_IMAGE044
，表示在
Figure 632296DEST_PATH_IMAGE045
时刻和位置x时预测到的变化率，其中，Ψ为观测到的降水场，t 0为外推的开始时间；
将平流分成N个时间步Δt来表示间隔时间 τ ，即τ =NΔt，对每个时间步的α 通过迭代的方法确定位移矢量
Figure 184500DEST_PATH_IMAGE046
；
从α = 0开始，最终位移矢量是单个时间步的N个矢量的矢量和，因此，在半拉格朗日格式中通过假设速度平稳性，即
Figure 590204DEST_PATH_IMAGE047
，沿着向前的时间线或向后的时间线来确定降水包的轨迹；
根据平流层光流场DARTS的计算方法获得的回波运动场在公式
Figure 202451DEST_PATH_IMAGE048
中进行迭代收敛，循环多次后中止迭代得到相应数据，其中，t表示外推时间，u是x方向上的速度场，v表示y方向的速度场；u，v都是上面的光流法得到的结果，
Figure 818240DEST_PATH_IMAGE049
是网格点
Figure 656359DEST_PATH_IMAGE050
处的回波运动。
3.根据权利要求2所述的一种基于气象雷达多层CAPPI的光流法外推方法，其特征在于：所述将不同高度的外推数据进行融合处理步骤如下：
根据得到的每个层的雷达回波预报数据，每个层的数据用二维网格的图像表示，每个格点的值为-1到127的反射率值，最终的数据张量的形状为(c, height, width)，其中c表示层数、height表示图像的高度、width表示图像的宽度；
每个外推时刻t 0+Δt的结果是通过在N的方向上取极大值，得到不同时刻的网格点面数据；最终的数据张量的形状为(1, height, width)，就是只剩下一层图像方便气象员可视化雷达图；
总的外推时刻为 t 0+NΔt，根据迭代过程就会得到N个形状为(1, height, width)的雷达外推数据。
CN202110853509.0A 2021-07-28 2021-07-28 一种基于气象雷达多层cappi的光流法外推方法 Active CN113296074B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110853509.0A CN113296074B (zh) 2021-07-28 2021-07-28 一种基于气象雷达多层cappi的光流法外推方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110853509.0A CN113296074B (zh) 2021-07-28 2021-07-28 一种基于气象雷达多层cappi的光流法外推方法

Publications (2)

Publication Number Publication Date
CN113296074A CN113296074A (zh) 2021-08-24
CN113296074B true CN113296074B (zh) 2022-02-22

Family

ID=77331142

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110853509.0A Active CN113296074B (zh) 2021-07-28 2021-07-28 一种基于气象雷达多层cappi的光流法外推方法

Country Status (1)

Country Link
CN (1) CN113296074B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111208517B (zh) * 2020-01-15 2023-08-15 成都信息工程大学 基于多普勒天气雷达的短临外推预报流场构建方法
CN115542431B (zh) * 2022-11-25 2023-03-10 成都远望探测技术有限公司 一种基于地基云雷达和卫星数据的对流初生监测方法
CN116953653B (zh) * 2023-09-19 2023-12-26 成都远望科技有限责任公司 一种基于多波段天气雷达组网回波外推方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6035057A (en) * 1997-03-10 2000-03-07 Hoffman; Efrem H. Hierarchical data matrix pattern recognition and identification system
WO2007022376A2 (en) * 2005-08-18 2007-02-22 Honeywell International Inc. Constant altitude plan position indicator display for multiple radars
CN108519631A (zh) * 2018-02-22 2018-09-11 青岛心中有数科技有限公司 降水强度预测方法
CN111521990A (zh) * 2020-05-11 2020-08-11 沈阳工业大学 一种基于多层雷达回波数据的降雨量分析方法
CN212160076U (zh) * 2020-04-21 2020-12-15 成都远望科技有限责任公司 一种激光雷达性能监控装置
CN112232232A (zh) * 2020-10-20 2021-01-15 城云科技（中国）有限公司 一种目标检测方法
CN112505707A (zh) * 2021-01-29 2021-03-16 成都远望探测技术有限公司 一种x波段双极化快速扫描相控阵天气雷达

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004037366A (ja) * 2002-07-05 2004-02-05 Mitsubishi Electric Corp データ処理方法及びこれを用いた気象レーダ装置
US7843378B2 (en) * 2008-03-04 2010-11-30 Colorado State University Research Foundation Dynamic and adaptive radar tracking of storms (DARTS)
CN104950186B (zh) * 2014-03-31 2018-06-12 乌托巴斯洞察公司 雷电预测的方法和装置
CN106932766B (zh) * 2017-04-27 2019-07-09 中国人民解放军海军航空大学 基于变参数广义结构的距离扩展目标自适应检测方法
CN108535731B (zh) * 2018-04-18 2020-12-29 青岛心中有数科技有限公司 短临降水预报方法与装置
US11442164B2 (en) * 2019-06-07 2022-09-13 Honeywell International Inc. Systems and methods for determining convective cell growth from weather radar reflectivity data
CN112379345B (zh) * 2020-10-23 2024-04-19 吴海英 融合数值模式的雷达短临外推预报方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6035057A (en) * 1997-03-10 2000-03-07 Hoffman; Efrem H. Hierarchical data matrix pattern recognition and identification system
WO2007022376A2 (en) * 2005-08-18 2007-02-22 Honeywell International Inc. Constant altitude plan position indicator display for multiple radars
CN108519631A (zh) * 2018-02-22 2018-09-11 青岛心中有数科技有限公司 降水强度预测方法
CN212160076U (zh) * 2020-04-21 2020-12-15 成都远望科技有限责任公司 一种激光雷达性能监控装置
CN111521990A (zh) * 2020-05-11 2020-08-11 沈阳工业大学 一种基于多层雷达回波数据的降雨量分析方法
CN112232232A (zh) * 2020-10-20 2021-01-15 城云科技（中国）有限公司 一种目标检测方法
CN112505707A (zh) * 2021-01-29 2021-03-16 成都远望探测技术有限公司 一种x波段双极化快速扫描相控阵天气雷达

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
Fully Spectral Method for Radar-Based Precipitation Nowcasting;Seppo Pulkkinen等;《 IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing》;20190531;第12卷(第5期);全文 *
Short-term predictability of weather radar quantities and lightning activity;Evan Ruzanski等;《 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)》;20151112;全文 *
Stochastic Spectral Method for Radar-Based Probabilistic Precipitation Nowcasting;Seppo Pulkkinen等;《Journal of Atmospheric and Oceanic Technology》;20190101;全文 *
光流法及其在临近预报中的应用;曹春燕等;《气象学报》;20151231;第73卷(第3期);全文 *
半拉格朗日方法及其在数值模拟和数值预报中的应用;董敏;《应用气象学报》;19970228;第8卷(第1期);全文 *
基于天气雷达反演的参量和光流法在回波外推中的应用研究;刘燕斐;《中国优秀硕士学位论文全文数据库 基础科学辑》;20200315;全文 *
基于改进光流法的雷达图像运动估计;王志斌等;《计算机技术与发展》;20170801(第12期);全文 *

Also Published As

Publication number Publication date
CN113296074A (zh) 2021-08-24

Similar Documents

Publication Publication Date Title
CN113296074B (zh) 一种基于气象雷达多层cappi的光流法外推方法
Brusch et al. Underwater bottom topography in coastal areas from TerraSAR-X data
CN104950306B (zh) 一种海杂波背景下前视海面目标角超分辨成像方法
CN102982555B (zh) 基于自适应流形粒子滤波的制导红外小目标跟踪方法
CN108828691A (zh) 短临降水预报方法与装置
Naaijen et al. Phase resolved wave prediction from synthetic radar images
CN108107434B (zh) 基于双多普勒雷达反演的区域三维风场拼图方法
CN108594230A (zh) 一种海船场景的合成孔径雷达图像仿真方法
Gerstoft et al. Refractivity estimation using multiple elevation angles
CN109100723A (zh) 基于多普勒天气雷达数据的高空风反演方法
CN116953653B (zh) 一种基于多波段天气雷达组网回波外推方法
Damiani et al. A high-resolution dual-Doppler technique for fixed multiantenna airborne radar
CN115267770A (zh) 一种sar图像海洋涡旋检测方法及系统
CN115755043A (zh) 一种基于x波段非相参雷达的波浪场重构及预测方法
Boag et al. A fast physical optics (FPO) algorithm for double-bounce scattering
KR101781760B1 (ko) 해양파 스펙트럼을 이용한 해양 표면의 사실적 시뮬레이션 방법 및 장치
CN112731564B (zh) 一种基于多普勒天气雷达数据的雷电智能预报方法
Fowdur et al. Tracking targets with known spatial extent using experimental marine radar data
CN111487621A (zh) 一种基于雷达图像的海表流场反演方法及电子设备
CN116451465A (zh) 一种星载sar中尺度涡成像仿真方法及系统
CN116500648A (zh) 一种地基激光雷达目标区风廓线反演方法
Kim et al. Real-time inverse estimation of multi-directional random waves from vessel-motion sensors using Kalman filter
Pfaff et al. Filtering on the unit sphere using spherical harmonics
CN114236522A (zh) 前向散射雷达网目标三维空间位置估计方法及存储介质
JP3589186B2 (ja) 機上海洋予察装置

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