CN116953653B - 一种基于多波段天气雷达组网回波外推方法 - Google Patents

Abstract

本发明涉及天气雷达组网领域，具体涉及一种基于多波段天气雷达组网回波外推方法，避免了组网时各部雷达的数据存在时间差异，降低了组网数据与真实回波存在的误差。本发明基于多波段天气雷达组网回波外推方法，包括：对多波段雷达数据进行衰减订正；计算雷达各个高度层数据；对雷达数据进行光流外推至同一时刻；时间同步后，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据；对组网后的数据采用RMP‑Net网络进行外推。本发明适用于对多波段天气雷达组网回波外推。

Description

一种基于多波段天气雷达组网回波外推方法

技术领域

本发明涉及天气雷达组网领域，具体涉及一种基于多波段天气雷达组网回波外推方法。

背景技术

天气雷达回波及相关产品可有效监测中小尺度的强对流天气。目前全国已建成上百部雷达，包含S、C、X等不同波段的雷达。然而，由于雷达的体积扫描模式、地形遮挡和地球曲率等因素，远距离雷达存在近地面盲区和近高度锥形区的盲区。因此，采用多部天气雷达组网，相互补盲以减少探测盲区，可以提高雷达监测效果。

现有技术如CN114660591A公开的一种基于多部天气雷达三维组网生成方法，其三维组网方法包括S1、分析组网雷达的数据探测范围，计算组网总网格范围及网格点数，最大探测高度及高度层数；S2、计算单部雷达在不同海拔高度层的数据；S3、计算单部雷达各个海拔高度层的数据与总网格的对应关系；S4、计算多部雷达在各个海拔高度层的融合数据，得到多部天气雷达的三维组网数据。

该方法具有根据组网雷达个数动态计算组网网格大小，所有组网雷达独立处理各层数据及映射关系加快计算速度，快速处理多雷达重复覆盖区域的数据。

但该方法中雷达数据时间不同步，导致组网时各部雷达的数据存在时间差异，导致组网数据与真实回波存在较大误差，甚至出现错误，影响对天气的掌握。

该方法也无法实现对组网后的数据进行外推，获得短期外推的回波。因此该方法对天气监测准确性不高。

发明内容

本发明的目的在于克服现有技术的缺点，提供一种基于多波段天气雷达组网回波外推方法，将超出时间范围的雷达数据外推至同一时刻，避免了组网时各部雷达的数据存在时间差异，导致组网数据与真实回波存在较大误差，并且实现了对组网后的数据进行外推，获得短期外推的回波，提高了对天气监测的准确性。

本发明采取如下技术方案实现上述目的，一种基于多波段天气雷达组网回波外推方法，包括：

S1、对多波段雷达数据进行衰减订正；

S2、计算雷达各个高度层数据；

S3、对雷达数据进行光流外推至同一时刻；

S4、时间同步后，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据；

S5、对组网后的数据采用RMP-Net网络进行外推。

进一步的是，计算雷达各个高度层数据具体包括：

S21、根据测高公式计算CAPPI面上某径向上点A的仰角a；

S22、判断点A的仰角a与体扫各层仰角的关系；

S23、根据仰角关系进行内插取值；

S24、对每个径向上的点都执行S21-S23的计算。

进一步的是，根据仰角关系进行内插取值具体包括：

若a刚好等于某仰角值，则直接取相应仰角PPI上的值为点A的CAPPI值；若a小于最低体扫仰角值，则取最低仰角PPI上的值为点A的CAPPI值；若a大于最高体扫仰角值，则视点A无回波值；若a在两个体扫仰角之间，则进行线性插值，a为点A所在的仰角。

进一步的是，对雷达数据进行光流外推至同一时刻具体包括：

S31、通过光流法获得外推风场；

S32、对获得的外推风场进行无辐散约束，获得无辐散风场；

S33、通过半拉格朗日外推方法将原始雷达数据外推至同一时刻。

进一步的是，通过半拉格朗日外推方法将原始雷达数据外推至同一时刻具体包括：

根据不同时刻的回波图获得雷达回波的移动矢量，根据移动矢量进行半拉格朗日外推，方式如下：；

其中为外推回波强度，F为实际观测值，t₀为开始时间，τ为外推时效，x为空间位置，α为像素点在时间τ内的位移矢量，即当前位置未来时刻的回波是由当前时候位置上的回波平移得到；

通过将预报时效分为多步进行外推方式将雷达数据进行光流外推至同一时刻。

进一步的是，多步外推时，每一步的时长为Δt，位移矢量α由如下公式迭代得到：a⁽ⁿ⁺¹⁾=ΔtU(t₀, x - a⁽ⁿ⁾)，其中U(t₀, x - a⁽ⁿ⁾)为云图在x-a处的速度，像素点总的位移矢量是多步位移的矢量和；

若是二维的情况，则：

a⁽ⁿ⁺¹⁾=ΔtU[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]；

β⁽ⁿ⁺¹⁾=ΔtV[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]；

初值a⁽⁰⁾和β⁽⁰⁾设置为0，经过多步迭代后求得α和β。

进一步的是，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据具体包括：

获取三维组网网格体系：

根据动态输入的雷达站点的体扫数据，获取每个雷达站的基本信息，包括站点经纬度和最大探测距离；判断每个雷达站点经纬度值大小，遍历所有雷达站点的最大探测距离，得到所有雷达站点中的最大探测距离的最大值；将站点经纬度的范围外扩最大探测距离的最大值，将外扩后的站点经纬度的范围作为组网网格的经纬度范围；

根据组网网格的经纬度范围以及网格分辨率计算二维平面上的横向网格数和纵向网格数，结合雷达海拔高度和探测距离，得到以海平面为基准的三维组网的最大探测高度；

根据已计算的最大探测高度和预设的高度分辨率，计算三维网格的高度层数；将所有高度层的二维平面融合，形成三维组网网格体系，每个三维网格的格点用于保存雷达的一个回波数据；

单部雷达组网融合：

根据雷达最大探测距离，以外接正方形形式，计算当前雷达探测的经纬度范围；根据雷达探测的经纬度范围对应组网网格位置关系，所述位置关系包括起始、结束网格序号、横向以及纵向网格数；根据每个扫描点的方位和距离雷达中心的距离，计算每一个扫描点对应的经纬度值，获取单部雷达不同高度层上的网格数据，再根据单部雷达网格与组网网格对应关系将单部雷达数据融合到组网网格中；

多部雷达组网融合：

对三维组网网格不同的高度层进行分别处理，各个高度层的间隔为预设的高度分辨率；获取各个高度层二维网格内的所有雷达站点对应的CAPPI结果值，对获取的结果值进行融合处理，融合方式如下：

，/>是单元网格中来自第n个雷达的分析值，w_n为对应给与的权重，N_rad为网格单元处有分析值的雷达总个数。

进一步的是，对组网后的数据采用RMP-Net网络进行外推具体包括：

对组网后的数据进行归一化处理，得到对应的雷达回波灰度图像，并生成为原始雷达回波图像的时间；

设置输入雷达回波灰度图像的序列长度、总时长以及相邻两张雷达回波图像的间隔时间；

设置目标雷达回波灰度图像数据的序列长度、总时长以及相邻两张雷达回波图像的间隔时间；

设置完成后，通过输入雷达回波灰度图像对RMP-Net网络进行训练，训练完成后，对组网后的数据采用RMP-Net网络进行外推。

进一步的是，所述RMP-Net网络包含LERM 模块，所述LERM 模块包含用于存储回声演化规律记忆的记忆池以及用于记忆回忆的感知注意机制；记忆池具有张量形式，表示为，其中K为记忆项的数量，每个记忆项的特征为m_i；

外推具体包括：给定输入回波序列，提取其时空动态特征，并使用两层感知结构将提取的时空动态特征投影到查询向量中q中，，其中u_m为差分输入回波序列上应用2D-CNNs得到的回波运动特征，u₁为直接应用3D-CNNs得到的回波生命周期特征，fc₁和fc_m都是用于仿射变换的全联通层，W_m和W₁为线性变换矩阵，b_q为偏置项；

使用查询向量q从记忆池中召回记忆特征项，并与感知注意机制保持一致，然后对记忆项加权求和：，m^r为回忆回声-进化记忆特征，a_i是对每个记忆项m_i的注意权重；

记忆对齐后，召回的m^r是与输入回波动态相对应的最相关的记忆特征，然后将该最相关的记忆特征嵌入到外推主干中，指导外推。

进一步的是，步骤S5还包括对RMP-Net网络进行评估，具体为：选择检测概率、误报率、关键成功指数和结构相似性指数度量进行定量评估，检测概率、误报率、关键成功指数用于降水临近预报性能评估，结构相似性指数用于评估生成的图像的感知质量；

检测概率POD、误报率FAR以及关键成功指数CSI计算方式如下：

，

其中，TP、FP 和 FN分别表示真阳性、假阳性以及假阴性，真阳性表示外推=真值=1，假阳性表示外推=1，真值=0，假阴性表示外推=0，真值=1；检测概率、关键成功指数和结构相似性指数度量得分越高表示RMP-Net网络性能越好，而误报率得分越低表示RMP-Net网络性能越好。

本发明的有益效果为：

本发明将超出时间范围的雷达数据外推至同一时刻，避免了组网时各部雷达数据存在的时间差异，降低了组网数据与真实回波存在的误差。

本发明实现了对组网后的数据进行外推，获得了短期外推的回波，提高了对天气监测的准确性。

附图说明

图1是本发明实施例提供的一种基于多波段天气雷达组网回波外推方法流程图；

图2是本发明实施例提供的CAPPI计算高度示意图；

图3是本发明实施例提供的VCP21扫描模式下的9个仰角及在标准大气下的波束传播路径示意图；

图4是本发明实施例提供的两部雷达在VCP21扫描模式下组网覆盖情况示意图；

图5是本发明实施例提供的多部雷达组网经纬度范围计算关系示意图；

图6是本发明实施例提供的经纬度距离计算示意图；

图7是本发明实施例提供的多部雷达组网高度范围计算关系示意图；

图8是本发明实施例提供的单雷达经纬度范围映射网格数据示意图；

图9是本发明实施例提供的三维组网网格数据示意图；

图10是本发明实施例提供的RMP-Net的整体架构示意图；

图11是本发明实施例提供的RMP-Unit内部结构；

图12是本发明实施例提供的LERM模块结构图示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。

如图1所示，本发明提供一种基于多波段天气雷达组网回波外推方法，包括：

S1、对多波段雷达数据进行衰减订正

具体采用KDP综合分类法对X波段雷达的反射率数据进行衰减订正。

S2、计算雷达各个高度层数据

具体包括：根据三维组网中各个高度层对应海拔高度，应用测高公式，选取临近该高度平面上的上下两个仰角相应雷达测距上的数据，然后用内插方法得到该高度上的数据（即CAPPI，等高平面位置显示，CAPPI是等高平面位置显示的英文名称，是基于气象雷达计算出的等高平面数据，CAPPI是一个投影面，该投影面上的数据是来自于切割平面与各个扫描层相交处的数据）。

第一步：根据测高公式计算CAPPI面上某径向上点A的仰角，测高公式如下：；其中：h为CAPPI高度，a为点A所在的仰角，SlatRan为点A到雷达中心的径向距离，H为站点天线海拔高度，R为地球平均半径；

第二步：判断点A的仰角a与体扫各层仰角的关系；

第三步：根据仰角关系进行内插取值，如果a刚好等于某仰角值，则直接取相应仰角PPI上的值为该点的CAPPI值；如果a小于最低体扫仰角值，则取最低仰角PPI（PlainPosition Indicator，平面位置显示）上的值为该点的CAPPI值；如果a大于最高体扫仰角值，则视该点无回波值；如果a在两个体扫仰角之间，则进行线性插值；

第四步：对每个径向每个距离库都进行第一至三步计算。

线性插值：如图2所示，如果a在两个体扫仰角（a1，a2）之间，则对点A进行垂直方向上线性插值。首先根据点A、B、E所对应的水平距离相同，计算出B、E在对应仰角层对应的径向距离和到水平面的高度，再判断点B、E对应的回波值是否有效，最后根据高度权重进行插值。

S3、对雷达数据进行光流外推至同一时刻

第一步：通过Farneback光流法获得外推风场。

假设图像梯度恒定且局部光流恒定，将图像视为二维信号的函数，利用二次多项式进行近似建模：

；

采用全局位移d构建一个新信号：

，

假设A₁非奇异，由f₂(X)求解可得位移d：；在真实情况下，通过局部多项式近似代替全局多项式，可将误差控制在足够小的范围。

首先对两个图像进行多项式展开，为第一图像提供展开系数A₁(X)、b₁(X)和c₁(X)，为第二图像提供A₂(X)、b₂(X)和c₂(X)。理想情况下，应该A₁=A₂，但在实际实践中我们必须满足于近似：

在X的邻域I上尽可能找到满足A(X)d(X)的最小化函数：

其中ω(∆xX)是邻域中点的权重函数，获得的最小值为：，该公式中，只有Δb未知，其余项可直接或间接求得，为便于求解Δb，引入先验位移估计的/>，令/>，则：

计算时，将d(X)设置为0，联立上述公式，经过多次迭代以后，可计算得到d(X)的精确值。为兼顾计算效率和光流估计精度，设定金字塔层数为3层。

第二步：对获得的外推风场进行无辐散约束，获得无辐散风场。

对光流法反演出的风场进行无辐散约束可以减少风矢量杂乱无序的现象，并且可以对等速进行平滑处理，使风速空间梯度更平滑。

假设光流法获得的风场在x，y方向上的分量为u^T (x，y)和v^T (x，y)，无辐散处理后的风场为u^c (x，y) 和v^c (x，y)，使其既要逼近原始风场u^T (x，y), v^T (x，y)，又要满足二维连续方程(约束条件)：

，构造泛函J (u^c, v^c)，则

；

公式中∑为风场的空间范围，可通过求解泛函J (u^c, v^c)的最小值实现对u^c (x，y)和v^c (x，y)的求解。在约束条件下问题转化为求解如下泛函的最小值：

上式中λ为拉格朗日乘子，F的极小值可通过欧拉-拉格朗日方程求解：

，将u^c、v^c代入构造泛函J (u^c, v^c)可得：

，该函数为泊松方程，写为差分形式后可用高斯迭代法进行求解，解出λ后带入u^c、v^c方程即可求得无辐散风场u^c (x，y)和v^c (x，y)。

第三步：通过半拉格朗日外推方法将原始雷达数据外推至同一时刻。

光流法可以根据不同时刻的回波图想获得雷达回波的移动矢量，根据移动矢量进行外推预报，结合无辐散风场，利用半拉格朗日外推方法，能够保持螺旋状回波的旋转状态。半拉格朗日外推方法如下：

其中为外推回波强度，F为实际观测值，t₀为开始时间，τ为外推时效，x为空间位置，a为像素点在时间τ内的位移矢量，即当前位置未来时刻的回波是由当前时候特定位置上的回波平移得到；半拉格朗日平流方案将预报时效分为N步进行外推（这样即可得到外推至任意步长时刻），每一步的时长为Δt，位移矢量α由下面的公式迭代得到：

a⁽ⁿ⁺¹⁾=ΔtU(t₀, x - a⁽ⁿ⁾)，其中，U(t₀, x - a⁽ⁿ⁾)为云图在x-a处的速度，像素点总的位移矢量是N步位移的矢量和，因此半拉格朗日外推法获得的运动轨迹是流线型，具有旋转效应。若是二维的情况则可改写为：

a⁽ⁿ⁺¹⁾=ΔtU[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]，

β⁽ⁿ⁺¹⁾=ΔtV[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]，

上式的初值a⁽⁰⁾和β⁽⁰⁾可设置为0，经过多步迭代后可求得a和β，β也表示位移矢量。

S4、时间同步后，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据；

单部天气雷达扫描方式及空间分辨率：

单部天气雷达基数据大多都是VCP21体扫数据，如图3，由于扫描方式的局限性，在最高扫描19.5°仰角以上和最低扫描仰角0.5°以下有资料空白区（图3中灰色区域），仰角之间也有间隙，还受分辨率、探测范围和地形阻挡等影响，导致资料不完整，而多部雷达组网拼图则可以弥补单部雷达的资料空白区（如图4所示），红色和黑色线条各为一部雷达的VCP21的扫描范围，相对于图3明显改善单部雷达的资料空白区域。

计算组网网格范围及网格数：

由于参与拼图的雷达个数根据实时扫描情况动态输入，导致组网拼图经纬度范围和最大探测高度可能不一致，所以需要动态计算组网经纬度网格范围及网格数，最大探测高度及高度层数。

根据动态输入的雷达站点（R1、R2、R3、R4等）的体扫数据，获取每个雷达站的基本信息，包括站点经纬度和最大测距；然后判断站点经纬度值大小，获取最小经度(minLon)、最大经度(maxLon)、最小纬度(minLat)、最大纬度(maxLat)，如图5所示的虚线所表示得到站点经纬度的范围；遍历所有雷达站点的最大探测测距，得到最大探测距离的最大值(maxD)。多部雷达组网经纬度范围计算关系示意图如图5所示。

然后，根据站点经纬度的范围外扩最大测距，获取组网网格的经纬度范围，即最小经度(minLonPZ)、最大经度(maxLonPZ)、最小纬度(minLatPZ)、最大纬度(maxLatPZ)。假设小范围内，经纬度线近似于直线，故可采用左上角点和右下角点外扩最大测距，得到最大最小经纬度值作为组网网格的经纬度范围。

本发明使用已知一点的经纬度（long1，lat1），距离（d）和方位角（α），计算另一点的经纬度（long2，lat2），如图6所示，具体计算如下：

1、计算第二点的经度，即水平平移的距度（d*sinα）除以当前纬度切面周长（2π*arc），再每乘以360度，就能够获得水平横向平移了多少度，再加上long1，就得到long2的值了。

2、计算第二点的纬度，即垂直平移的距离d（d*cosα）除以地球纵向周长，再乘上360度，就能够获取纵向平移了多少度，再加上lat1，就得到lat2的值了。

最后，根据组网网格经纬度范围，以及网格分辨率（LenofWin），计算二维平面上的横向网格数（MaxWidth=（maxLon-minLon）/ LenofWin）和纵向网格数（MaxHeight=（maxLat-minLat）/ LenofWin）。

另外，根据动态输入的雷达站点数据（R1、R2、R3、R4等），获取每个雷达站的基本信息，包括站点经纬度和最大探测测距，如图7所示，例如D1、D2、D3和D4为雷达站点R1、R2、R3和R4的最大探测距离，结合雷达海拔高度和探测距离，得到以海平面为基准的三维组网的最大探测高度(maxH)。

根据已计算的最大探测高度（maxH）和预设的高度分辨率（LenofHei），计算三维网格的高度层数（HeightNum= maxH/ LenofHei）。

最后，所有高度层的二维平面融合，形成三维组网网格体系，三维组网网格总数为：MaxWidth* MaxHeight* HeightNum，每个三维网格的格点用于保存雷达的一个回波数据。

单站径向数据与网格的对应关系：

由于雷达体积扫描形式（如VCP21模式）决定雷达数据为三维锥面的极坐标形式，导致单部雷达计算得到的CAPPI结果也是极坐标形式的（投影到平面，仰角为0），与三维组网网格的笛卡尔坐标不一致，需要进行坐标转换。

已知单部雷达的站点经纬度，首先根据雷达最大探测距离，以外接正方形形式，计算当前雷达探测的经纬度范围，计算方法同组网网格范围计算，形成单部雷达的网格。然后根据单部雷达探测经纬度范围，对应组网网格位置关系，包括起始、结束网格序号，横向、纵向网格数。

根据每个距离库（扫描点）的方位、离雷达中心的距离，可利用long2、lat2的计算方式计算每一个距离库对应的经纬度值，计算获取单部雷达不同高度层上的网格数据，再根据单部雷达网格与组网网格对应关系（如图8所示），将单部雷达数据融合到组网网格中。

图8中，R为示例雷达站点，A、B、C、D为雷达在高度层H（指任意高度）所计算出的反射率数据，A1、B1、C1、D1为雷达映射到高度层H的组网网格位置。

由于多个高度层上的CAPPI对应于网格的位置是固定的，所以考虑将每一部雷达的极坐标与组网网格的笛卡尔坐标对应关系保存到缓存，不再重复计算，加快计算效率。

多雷达同一高度层数据融合：

在三维组网网格的很多区域，特别是在对流层中高层，有来自多个雷达的资料重叠区，在每个高度层上的组网网格中的每个单元点的值需要进行合成：

，式中，/>是单元网格中来自第n个雷达的分析值，w_n为对应给与的权重，N_rad为网格单元处有分析值的雷达总个数。

为了避免噪声的干扰，反射率小于0dBZ的格点认为无回波点。如果N_rad=0，那么该网格点不被任何一个雷达覆盖，直接赋值为无效值；如果N_rad=1，那么网格单元值直接等于对应雷达分析值；如果N_rad＞1，需要对多个雷达的分析值的权重平均进行融合处理。目前采用最大值法，即覆盖同一网格单元的多个雷达反射率分析值中的最大值的权重赋值为1，其他的权重全赋值为0，也就是把覆盖同一网格单元的多个雷达反射率分析值中的最大值赋值给网格单元。

针对三维组网网格不同的高度层（H1，H2…Hn）进行分别处理，各个层的间隔为预设的高度分辨率，即对每一个高度层上的二维网格（MaxWidth* MaxHeight）内的数据进行识别结果值融合分析，首先获取可落在当前二维平面网格单元内的所有雷达站点对应的CAPPI结果值，然后对这些结果值进行上述数据融合处理，得到当前网格单元的最终结果值。当所有层的网格单元处理完成，即可得到完整的多部天气雷达三维组网网格数据结果，如图9所示。

S5、对组网后的数据采用RMP-Net网络进行外推。

组网后的多普勒气象雷达数据归一化处理得到的雷达回波灰度图像，图像名称为雷达回波探测并生成为原始雷达回波图像的时间。模型输入雷达回波灰度图像数据的序列长度为10张，总时长为1小时，相邻两张雷达回波图像的间隔时间为6分钟；目标雷达回波灰度图像数据的序列长度为20张，总时长为2小时，相邻两张雷达回波图像的间隔时间为6分钟。总雷达回波图像序列长度(模型输入雷达回波图像序列长度加上目标雷达回波图像序列长度)为30，覆盖时长为3小时。设置完成后开始对RMP-Net网络模型进行训练。训练完成后，对组网后的数据采用RMP-Net网络进行外推。

整体架构：

本发明的RMP-Net的整体框架如图10所示。它利用PredRNN（引用）的结构并堆叠多个RMP-unit来生成从时间步长2~T的预测。在任何时间t，模型都会预测下一个时间戳t+1雷达图，根据当前雷达图Xt和历史雷达序列X0:t。红色和蓝色箭头分别表示空间记忆M和时间记忆C的传递方向。这两个记忆分别保存空间和时间表征。与PredRNN不同，RMP-Net利用不同的数据流来传输长期时空信息/>保留所有历史表征。

RMP-Unit模块：

感知记忆力预测单元(RMP-Unit)的内部结构如图11所示。输入包括当前输入，之前的隐藏状态/>，时间记忆/>、空间记忆/>和长期记忆/>。根据图10可知，RMP-Net由四个堆叠的RMP-Unit组成。在底层，/>表示所有历史输入X0: t，而在其他层，/>是上一层的输出。RAP-Unit 的输出是当前隐藏状态/>、空间记忆/>，时间记忆/>，以及新的长期记忆/>。详细计算公式如下：

由于雷达回波形状不规则且分布不同，RAB可以更好地捕获不同雷达回波之间的相关性。因此，该块的引入可以提高模型的时空能力，特别是因为高反射率的雷达回波的信息更容易提取，因为它们具有更稳定的外观和形状。

长期演化规律性记忆模块：

由于复杂的时空回波演化模式，根据历史观测推断未来回波是一项具有挑战性的任务。现有外推方法存在的一个固有范式缺陷是它们仅关注基于有限观测动力学的回波运动外推，而未能全面模拟回波演化过程。

由于专业预报员通常可以在实践中掌握当地天气条件的独特演变规律，然后用它们来指导未来的预报，因此本发明在网络中加入一种新颖的 LERM 模块来模拟这种行为，以补充缺陷回波运动外推范式。LERM模块由两个部分组成，一个是用于存储回声演化规律记忆的记忆池，另一个是用于记忆回忆的感知注意机制（如图12）所示。

记忆池具有张量形式表示为，其中K为记忆项的数量，每个记忆项的特征为m_i。它可以看作是多种代表性回波演化模式的密集表示的集合，任何其他复杂模式都可以通过加权张量加法运算来表示，这对应于自然大气系统演化模式通常是各种典型回波演化模式的组合子模式。一旦内存池构建完毕，给定输入回波序列，我们首先提取其时空动态特征，并使用两层感知结构将它们投影到查询向量中q中：

其中u_m为差分输入回波序列上应用2D-CNNs得到的回波运动特征，u₁为直接应用3D-CNNs得到的回波生命周期特征。fc₁和fc_m都是用于仿射变换的全联通层。W_m和W₁为线性变换矩阵，b_q为偏置项。然后使用形成的查询向量q从记忆池中召回记忆特征项，这与感知注意机制Att保持一致。回忆回声-进化记忆特征m^r是由注意权重(即进化模式组合)对记忆项目的加权求和：

，其中a_i是对每个记忆项m_i的注意权重：/>

其中g表示沿通道维度的全局平均池化，d[·，·]计算查询向量与每个通道池化内存项之间的余弦相似度，代表 softmax功能。记忆对齐后，召回的m^r是与输入回波动态相对应的最相关的记忆特征，提供了回波未来如何发展的信息，然后嵌入到外推主干中，指导外推。

由于LERM模块的所有计算都是完全可微的，因此基于历史观测的地面真相，在训练过程中可以通过梯度反向传播优化记忆参数以自动收敛，并且所提出的网络可以学习隐式识别和存储回声进化模式。因此，在网络中采用正态分布初始化内存参数。

网络模型评估：

选择检测概率（POD）、误报率（FAR）、关键成功指数（CSI）和结构相似性指数度量（SSIM）进行定量评估，其中常用的是POD、FAR和CSI用于降水临近预报性能评估。使用 SSIM来评估生成的图像的感知质量。本发明根据列联表得出 POD、FAR 和 CSI 的技能得分。给定对数雷达反射率因子阈值τ，使用 TP、FP 和 FN 来表示真阳性（外推 = 真值 = 1）、假阳性（外推 = 1，真值 = 0）和假阴性（外推 = 0，真值 = 1）。POD、FAR 和 CSI 计算如下：

其中，x和y是待比较的两幅图像，μ_x和μ_y分别是x和y的均值，σ_x和σ_y分别是x和y的标准差，σ_xy是x和y的协方差，C₁和C₂是两个常数，用于避免分母为0。

所有指标均在 [0, 1] 范围内。POD、CSI 和 SSIM 得分越高表示模型性能越好，而 FAR 得分越低越好。此模型评估得出了CSI 在三个雷达反射率阈值上的技能得分，τ{10 , 25 , 35 }，分别对应轻、中、重级别。

0-2小时指标均值图下表：

各指标12min间隔变化情况如下表：

以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。

Claims (8)

1.一种基于多波段天气雷达组网回波外推方法，其特征在于，包括：

S1、对多波段雷达数据进行衰减订正；

S2、计算雷达各个高度层数据；

S3、对雷达数据进行光流外推至同一时刻；

S4、时间同步后，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据；

S5、对组网后的数据采用RMP-Net网络进行外推，具体包括对组网后的数据进行归一化处理，得到对应的雷达回波灰度图像，并生成为原始雷达回波图像的时间；

设置输入雷达回波灰度图像的序列长度、总时长以及相邻两张雷达回波图像的间隔时间；

设置目标雷达回波灰度图像数据的序列长度、总时长以及相邻两张雷达回波图像的间隔时间；

设置完成后，通过输入雷达回波灰度图像对RMP-Net网络进行训练，训练完成后，对组网后的数据采用RMP-Net网络进行外推；

所述RMP-Net网络包含LERM 模块，所述LERM 模块包含用于存储回声演化规律记忆的记忆池以及用于记忆回忆的感知注意机制；记忆池具有张量形式，表示为，其中K为记忆项的数量，每个记忆项的特征为m_i；

外推具体包括：给定输入回波序列，提取其时空动态特征，并使用两层感知结构将提取的时空动态特征投影到查询向量中q中，，其中u_m为差分输入回波序列上应用2D-CNNs得到的回波运动特征，u₁为直接应用3D-CNNs得到的回波生命周期特征，fc₁和fc_m都是用于仿射变换的全联通层，W_m和W₁为线性变换矩阵，b_q为偏置项；

使用查询向量q从记忆池中召回记忆特征项，并与感知注意机制保持一致，然后对记忆项加权求和：，m^r为回忆回声-进化记忆特征，a_i是对每个记忆项m_i的注意权重，/>其中g表示沿通道维度的全局平均池化，d[·，·]计算查询向量与每个通道池化内存项之间的余弦相似度；

记忆对齐后，召回的m^r是与输入回波动态相对应的最相关的记忆特征，然后将该最相关的记忆特征嵌入到外推主干中，指导外推。

2.根据权利要求1所述的基于多波段天气雷达组网回波外推方法，其特征在于，计算雷达各个高度层数据具体包括：

S21、根据测高公式计算CAPPI面上某径向上点A的仰角a；

S22、判断点A的仰角a与体扫各层仰角的关系；

S23、根据仰角关系进行内插取值；

S24、对每个径向上的点都执行S21-S23的计算。

3.根据权利要求2所述的基于多波段天气雷达组网回波外推方法，其特征在于，根据仰角关系进行内插取值具体包括：

若a刚好等于某仰角值，则直接取相应仰角PPI上的值为点A的CAPPI值；若a小于最低体扫仰角值，则取最低仰角PPI上的值为点A的CAPPI值；若a大于最高体扫仰角值，则视点A无回波值；若a在两个体扫仰角之间，则进行线性插值，a为点A所在的仰角。

4.根据权利要求1所述的基于多波段天气雷达组网回波外推方法，其特征在于，对雷达数据进行光流外推至同一时刻具体包括：

S31、通过光流法获得外推风场；

S32、对获得的外推风场进行无辐散约束，获得无辐散风场；

S33、通过半拉格朗日外推方法将原始雷达数据外推至同一时刻。

5.根据权利要求4所述的基于多波段天气雷达组网回波外推方法，其特征在于，通过半拉格朗日外推方法将原始雷达数据外推至同一时刻具体包括：

根据不同时刻的回波图获得雷达回波的移动矢量，根据移动矢量进行半拉格朗日外推，方式如下：；

其中为外推回波强度，F为实际观测值，t₀为开始时间，τ为外推时效，x为空间位置，α为像素点在时间τ内的位移矢量，即当前位置未来时刻的回波是由当前时候位置上的回波平移得到；

通过将预报时效分为多步进行外推方式将雷达数据进行光流外推至同一时刻。

6. 根据权利要求5所述的基于多波段天气雷达组网回波外推方法，其特征在于，多步外推时，每一步的时长为Δt，位移矢量α由如下公式迭代得到：a⁽ⁿ⁺¹⁾=ΔtU(t₀, x - a⁽ⁿ⁾)，其中U(t₀, x - a⁽ⁿ⁾)为云图在x-a处的速度，像素点总的位移矢量是多步位移的矢量和；

若是二维的情况，则：

a⁽ⁿ⁺¹⁾=ΔtU[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]；

β⁽ⁿ⁺¹⁾=ΔtV[x-a⁽ⁿ⁾,y-β⁽ⁿ⁾, t₀]；

初值a⁽⁰⁾和β⁽⁰⁾设置为0，经过多步迭代后求得α和β，β也表示位移矢量。

7.根据权利要求1所述的基于多波段天气雷达组网回波外推方法，其特征在于，将不同波段雷达数据进行融合，获得多个天气雷达的三维组网数据具体包括：

获取三维组网网格体系：

根据动态输入的雷达站点的体扫数据，获取每个雷达站的基本信息，包括站点经纬度和最大探测距离；判断每个雷达站点经纬度值大小，遍历所有雷达站点的最大探测距离，得到所有雷达站点中的最大探测距离的最大值；将站点经纬度的范围外扩最大探测距离的最大值，将外扩后的站点经纬度的范围作为组网网格的经纬度范围；

根据组网网格的经纬度范围以及网格分辨率计算二维平面上的横向网格数和纵向网格数，结合雷达海拔高度和探测距离，得到以海平面为基准的三维组网的最大探测高度；

根据已计算的最大探测高度和预设的高度分辨率，计算三维网格的高度层数；将所有高度层的二维平面融合，形成三维组网网格体系，每个三维网格的格点用于保存雷达的一个回波数据；

单部雷达组网融合：

根据雷达最大探测距离，以外接正方形形式，计算当前雷达探测的经纬度范围；根据雷达探测的经纬度范围对应组网网格位置关系，所述位置关系包括起始、结束网格序号、横向以及纵向网格数；根据每个扫描点的方位和距离雷达中心的距离，计算每一个扫描点对应的经纬度值，获取单部雷达不同高度层上的网格数据，再根据单部雷达网格与组网网格对应关系将单部雷达数据融合到组网网格中；

多部雷达组网融合：

对三维组网网格不同的高度层进行分别处理，各个高度层的间隔为预设的高度分辨率；获取各个高度层二维网格内的所有雷达站点对应的CAPPI结果值，对获取的结果值进行融合处理，融合方式如下：

，/>是单元网格中来自第n个雷达的分析值，w_n为对应给与的权重，N_rad为网格单元处有分析值的雷达总个数。

8.根据权利要求1所述的基于多波段天气雷达组网回波外推方法，其特征在于，步骤S5还包括对RMP-Net网络进行评估，具体为：选择检测概率、误报率、关键成功指数和结构相似性指数度量进行定量评估，检测概率、误报率、关键成功指数用于降水临近预报性能评估，结构相似性指数用于评估生成的图像的感知质量；

检测概率POD、误报率FAR以及关键成功指数CSI计算方式如下：

，

其中，TP、FP 和 FN分别表示真阳性、假阳性以及假阴性，真阳性表示外推=真值=1，假阳性表示外推=1，真值=0，假阴性表示外推=0，真值=1；检测概率、关键成功指数和结构相似性指数度量得分越高表示RMP-Net网络性能越好，而误报率得分越低表示RMP-Net网络性能越好。