CN114624715A - 一种基于自注意力时空神经网络模型的雷达回波外推方法 - Google Patents

一种基于自注意力时空神经网络模型的雷达回波外推方法 Download PDF

Info

Publication number
CN114624715A
CN114624715A CN202210193798.0A CN202210193798A CN114624715A CN 114624715 A CN114624715 A CN 114624715A CN 202210193798 A CN202210193798 A CN 202210193798A CN 114624715 A CN114624715 A CN 114624715A
Authority
CN
China
Prior art keywords
echo
layer
conv2d
input
size
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.)
Pending
Application number
CN202210193798.0A
Other languages
English (en)
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.)
National University of Defense Technology
Original Assignee
National University of Defense 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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202210193798.0A priority Critical patent/CN114624715A/zh
Publication of CN114624715A publication Critical patent/CN114624715A/zh
Pending legal-status Critical Current

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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/047—Probabilistic or stochastic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • 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)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Image Processing (AREA)

Abstract

本发明提供了一种基于自注意力时空神经网络模型的雷达回波外推方法，包括：数据预处理：对经相关系数质量控制的原始PPI组合反射率雷达回波数据依次进行降噪、单帧分类、事件分类、归一化和历史未来段划分处理；训练神经网络模型：将训练集传入SAB‑STN中进行训练，训练完成后保存神经网络模型参数；测试神经网络模型：将测试集中的未标记网格数据输入到神经网络模型中，最终输出未来回波序列预测结果。

Description

一种基于自注意力时空神经网络模型的雷达回波外推方法
技术领域
本发明属于临近预报的自动化分析领域，具体涉及一种基于自注意力时空神经网络模型的雷达回波外推方法。
背景技术
临近预报指对当前至未来0-2小时的天气预报，是非常重要的天气预报业务之一，而临近预报最主要的预报对象则包括雷暴、降水和强对流，并常借助多普勒天气雷达产品与卫星云图资料进行，因为不仅在大多数情况下，这些预报对象的预警时效上限不大于2小时，而且应用数值模式来预报未来0-2小时的这些天气，误差通常极大。由于使用的是图像资料，临近预报仍以识别、追踪及预测雷达回波和云图的图像特征为基本点，但随着相关业务的发展成熟，临近预报也通过结合其他技术，其预报效果得以不断优化完善。作为临近预报的一种方式，雷达回波外推技术已在业务中被广泛应用，其所依赖的回波产品对应不同的业务需求具有多样性。其中单站预报常综合考虑平显(PPI，Plan Position Indicator)和体扫(RHI，Range Height Indicator)基本反射率回波，而较广区域的临近预报常利用由多个站点各自PPI基本反射率产品所生成的组合反射率图像。
当前雷达回波外推方法分为传统方法与机器学习方法，而传统方法又包括人工经验外推法，质心法，交叉相关法和光流法。具体而言，人工经验外推法仍是应用最广泛的方法，主要依赖预报员的知识与经验积累，能在某些预报任务，针对某些天气目标，结合某些非雷达数据综合分析，得出较为概括且具有一定可靠度的预测结果。人工经验外推法的局限在于不仅需要足够的专业知识，而且无法更直观、更准确地获得外推结果；质心法关注回波单体并追踪其运动，通常先从回波强度、形态、大小和分布等多种因素分析单体特征，再计算出单体质心在各时刻位置，并线性拟合出单体运动矢量与路径，进而外推出质心未来位置；交叉相关法先将整个图像划分多个子图，后在相邻回波帧的子图之间计算相关系数，再以最大相关系数获取不同时刻子图关联的方法确定作为外推依据的回波运动矢量；光流法基于光流模型分析回波运动，利用光流获得整个观测域内较为准确的运动矢量。以上传统方法的局限在于：(1)受限于各类人工规则、经验认知以及实地情况，无法更精确把握混沌系统；(2)基于运动矢量的线性外推，在回波生消演化方面缺少更多的预测细节，甚至对于强对流等变化急剧的天气系统，其仍缺乏令人满意的预测能力。相较于传统外推方法，基于机器学习的外推法善于从给定的数据分布中挖掘出隐式回波运动特征，并具有较好的非线性拟合能力，包括卷积神经网络CNN，循环神经网络RNN，条件生成对抗网络cGAN等等，但当前影响机器学习方法效果的最主要因素仍然是足量且高质量的数据，其次，机器学习算法的表现也因应用场景不同而不同。
发明内容
发明目的：本发明所要解决的技术问题是针对现有雷达回波外推传统方法依靠经验和运动矢量线性外推导致预报精度不足，以及已有机器学习方法对回波数据利用和回波运动特征提取两方面能力不足等问题，提出了一种基于数据驱动的方法设计出自注意力时空网络来实现雷达回波外推，以及一种与之匹配的训练数据集清洗方法，利用通过该数据清洗方法所得的大量历史回波序列训练神经网络模型来预测未来回波状态，可减少基于运动矢量线性外推所产生的误差和不稳定性，增强对混沌系统的感知与预测能力，提高雷达回波外推的智能化水平。包括以下步骤：
步骤1，数据预处理：对经相关系数质量控制的PPI组合反射率雷达回波产品进行数据预处理，将每个定长回波序列划分为历史段和未来段，得到训练集和测试集；
步骤2，训练神经网络模型：将训练集中的回波序列数据输入到基于自注意力时空神经网络模型(Self-Attention-Based Spatiotemporal Net)SAB-STN神经网络模型中进行训练，训练完成后保存神经网络模型参数；
步骤3，测试神经网络模型：将测试集中的历史回波序列数据输入到步骤2训练好的神经网络模型中进行神经网络模型测试，输出从历史回波向未来外推的预测结果。
步骤1包括以下步骤：
步骤1-1：回波图像预处理：对数据集中所有回波图像进行预处理，其中数据集来自国家气象数据网，包含近5年雷达组合反射率；
步骤1-2：回波序列预处理：根据回波图像的预处理结果，对所有回波序列进行预处理，包括序列类别标签标注、序列分类、序列类别均衡采样与序列的历史未来段划分；
步骤1-3：划分数据集：将预处理后的回波序列按年份将最近一年的数据划分为测试集，其余部分化为训练集。
步骤1-1包括以下步骤：
步骤1-1-1：图像降噪。去除回波图像中可疑非降水回波，包括以下步骤：
步骤1-1-1-1：将原始单通道(通道数channel＝1)灰度值回波图像xraw的像素值pix转化为回波强度值I，计算公式如下：
Figure BDA0003526079590000031
步骤1-1-1-2，根据数据来源设置降水回波强度最低阈值Ith，根据Ith与步骤1-1-1-1得到的回波强度值I对回波图像进行一次降噪，一次降噪规定如下：对于回波强度小于最低阈值Ith的像素和缺测像素NaN，将它们的值归零，其他情况则像素值不变，公式如下：
Figure BDA0003526079590000032
一次降噪完成之后，得到去除较低回波强度的单通道回波图像xraw；
步骤1-1-1-3：对一次降噪所得回波图像xraw进行连通域分析，统计得到每幅图像中各连通域面积ccA，连通域面积ccA，设定连通域面积最低阈值ccAth＝13，根据ccAth与ccA对回波图像X进行二次降噪，二次降噪规定如下：对于面积小于ccAth的连通域，连通域内像素值归零，否则像素值不变，公式如下：
Figure BDA0003526079590000033
二次降噪完成之后，得到不含较小面积杂波的单通道回波图像x；
步骤1-1-2：图像分类，包括以下步骤：
步骤1-1-2-1：统计每幅降噪后回波图像x中像素值不为0的像素个数numpix＞0，计算其占整个观测域像素个数numob的百分比，即有回波区与可观测区的面积比P1，公式如下：
Figure BDA0003526079590000041
根据P1是否为0首先将x二分类，其中，
Figure BDA0003526079590000042
表示降噪后无回波的回波图像，
Figure BDA0003526079590000043
表示降噪后有回波的回波图像，随后设置面积比区间集合
Figure BDA0003526079590000044
(集合内数值为百分数)：
Figure BDA0003526079590000045
并根据
Figure BDA0003526079590000046
的P1所属面积比区间，将
Figure BDA0003526079590000047
划分为4类；
步骤1-1-2-2：设置回波强度等级集合ΘI(集合内数值单位为dBZ)：
ΘI＝{(0，20]，(20，30]，(30，40]，(40，70]}，
统计每幅
Figure BDA0003526079590000048
和像素值在各回波等级的像素个数
Figure BDA0003526079590000049
计算与numpix＞0之间的百分比，即某回波等级区与有回波区的面积比Pi，公式如下：
Figure BDA00035260795900000410
根据按序排列后最后不为0的Pi所对应的i，将
Figure BDA00035260795900000411
划分为4类，随后结合步骤1-1-2-1中对
Figure BDA00035260795900000412
划分的4类，将
Figure BDA00035260795900000413
最终划分为16类，加上
Figure BDA00035260795900000414
则x共分为17类；
另外
Figure BDA00035260795900000415
所对应标签为(0，0)；
步骤1-1-3：进行图像标准化，公式如下：
Figure BDA00035260795900000416
其中，pix是原始像素值pixnorm是标准化像素值。
步骤1-2包括以下步骤：
步骤1-2-1：标注序列类别标签，制作对应回波序列的回波类别标签序列：
根据步骤1-1-2所得回波图像类别标签，将一个回波序列X中所有图像的类别标签按回波序列顺序排列，得到回波类别标签序列Xcl；
步骤1-2-2：序列分类：利用K-Means聚类算法，分别对回波类别标签序列Xcl的两个类别维度进行序列聚类，在两个类别维度上得到各自回波序列聚类结果，随后综合这两类聚类结果得到一个完整的聚类结果。其中，在K-Means聚类算法中所使用的序列样本距离为欧氏距离；
步骤1-2-3：对序列聚类结果进行类别均衡采样，回波序列采样公式如下：
Figure BDA0003526079590000051
其中，N为回波序列数据集的类别数，nj为类别j的样本总数，pi为从类别i中采样一个回波序列的概率；
步骤1-2-4：划分序列的历史未来段，设置历史帧数与未来帧数分别作为SAB-STN神经网络模型的输入长度Tinput与外推目标长度Ttarget，则Ttinput+Ttarget表示整个序列的长度，并以此将回波序列X划分为历史输入段Xinput与未来外推段Xtarget前后两部分：
Figure BDA0003526079590000052
其中xi表示序列中第i帧图像(i＝1，2，...，Tinput+Ttarget)；
将回波图像大小重塑：(height，width)重塑至(448，448)，height表示高度，width表示宽度，此时一个回波序列的数据形状(length，channel，height，width)为(Tinput，1，448，448)，表示为回波序列
Figure BDA0003526079590000053
表示长度L，通道数C，高度H，宽度W四个维度构成的高维实数空间，length表示长度，channel表示通道数。
步骤1-2-2包括以下步骤：
步骤1-2-2-1：选择K-Means聚类中心，设置类别数k，选择初始化的k个回波类别标签序列样本Xcl作为初始聚类中心cs：
cs＝cs1，cs2，...，csk；
其中csk指第k个聚类中心；
步骤1-2-2-2：归类每个样本。针对每个样本Xcl计算它到k个聚类中心的欧式距离dist并将其分到距离最小的聚类中心所对应的类中：
Figure BDA0003526079590000061
其中arg min函数指返回使被作用函数达到最小值的变量取值，c(i)表示第i个样例Xcl，i与k个类中距离最近的那个聚类中心csj，其中j＝1，2，..k，函数dist计算公式如下：
Figure BDA0003526079590000062
其中，a，b表示长度同为length的不同回波类别标签序列，ai，bi分别表示b，b在序列位置i处的回波标签向量；
步骤1-2-2-3：更新类的质心：针对每个类别，重新计算属于该类别所有样本的质心即聚类中心：
Figure BDA0003526079590000063
其中ci表示csi所对应的类，i＝1，2，..k；
步骤1-2-2-4：获得聚类结果。重复步骤1-2-2-3，直到收敛。此时确定k个类中各自包含的样本。
步骤2包括以下步骤：
步骤2-1：建立神经网络模型：搭建模型各模块及初始化模型参数；
步骤2-2：初始化训练超参数：采用Adam优化器训练神经网络模型，设置学习率lr、每次输入一组样本所含的样本数batchsize和训练集遍历数epoch；
步骤2-3：神经网络模型参数训练及更新：采用平均平方误差MSE作为损失函数，每次取一个batchsize的数据输入神经网络模型，以训练和优化神经网络模型参数，通过反向传播算法更新参数，模型参数初始化采用模型预训练和随机初始化相结合的方法，损失函数MSE计算公式如下：
Figure BDA0003526079590000064
其中pixi表示一幅真实回波图像内第i个像素，
Figure BDA0003526079590000071
表示真实回波图像对应生成回波图像内第i个像素，height×width表示图像内像素点个数。
步骤2-1包括以下步骤：
步骤2-1-1：构造编码器En：编码器为卷积神经网络，用于将原始回波图像转化为较高级的图像表征；
步骤2-1-2：构造回波序列处理器ESP(Echo Series Processor)：序列处理器包括ESP编码器ESP-En和ESP解码器ESP-De两部分，以自注意力模块SAB为基本单元，用于学习一个回波序列内由编码器生成的各时刻图像表征之间的时间关联，以及各时刻图像表征自身的空间关联；
步骤2-1-3：构造解码器De：解码器为卷积神经网络，用于将序列处理器所生成的特征转化成外推结果。
步骤2-1-4：初始化模型参数：对于SAB-STN神经网络模型中的卷积核，其权重采用kaiming正态分布初始化方法，其偏执采用0值初始化方法；对于模型中的标准化层，其权重采用1值初始化方法，其偏执采用0值初始化方法。
步骤2-1-1包括以下步骤：
步骤2-1-1-1：编码器En包括三个卷积模块Conv2dEn，1，Conv2dEn，2，Conv2dEn，3，两层最大池化层MaxPool2dEn，1，MaxPool2dEn，2。
各模块以串行方式排列，顺序依次为Conv2dEn，1，MaxPool2dEn，1，Conv2dEn，2，MaxPool2dEn，2，Conv2dEn，3，其中，所有卷积模块均依次包含一层卷积层，一层非线性激活函数层LeakyReLU，而Cony2dEn，2与Conv2dEn，3各额外包含一层标准化层GroupNorm，并排在卷积层之前。Conv2dEn，1的输入数据尺寸为(448，448)，卷积层的卷积核大小为(3，3)，输入通道数为1，输出通道数为4，步长为1，补零为1，LeakyReLU的α参数为0.02，GroupNorm的组数为1；MaxPool2dEn，1为核大小为2、步长为2、补零为0、能够将特征图尺寸分别缩小1/2的最大池化层；Conv2dEn，2的输入数据尺寸为(224，224)，卷积层输入通道数为4，输出通道数为8，其他参数设置与Conv2dEn，1一致；MaxPool2dEn，2与MaxPoolEdEn，1一致；CvnvEdEn，3的输入数据尺寸为(112，112)，卷积层输入通道数为8，输出通道数为16，其他参数设置与Conv2dEn，1一致。此外，输入数据在进入Conv2dEn，1前由5维形状(batchsize，length，channel，height，width)重塑至4维(batchsize*length，channel，height，width)，表示为
Figure BDA0003526079590000081
此时形状为(batchsize，length，channel，height，width)＝(batchsize，Tinput，1,448,448)的输入数据经编码器En处理，得到形状为(batchsize*Tinput，16，112,112)的回波图像表征。
步骤2-1-2包括以下步骤：
步骤2-1-2-1：自注意力模块SAB包括一层自注意力层SA，一层非线性激活函数层GeLU，以及一层标准化层GroupNorm，其中SAB的核心在于自注意力层SA；
在自注意力层SA中，输入端同时传入查询项源与键值项源，查询项源通过查询卷积Conv2dQ生成查询项Q，键值项源通过键值卷积Conv2dKV生成键值项KV，Conv2dQ，Conv2dKV均为逐通道channel-wise二维卷积层，其中Conv2dQ的输入数据尺寸为(112,112)，卷积核大小为(3，3)，输入通道数与输出通道数相等且为16，步长为1，补零为1，组数为16，Conv2dKV的输入数据尺寸为(112，112)，卷积核大小为(3，3)，输入通通道数为16，输出通道数为32，步长为1，补零为1，组数为16，图块边长值regionsize，Q与KV在图像高宽维度分别切分成若干不重叠且边长值为regionsize的图块，并且KV沿channel维度均分为两份得到键项K和值项V，此时Q，K，V三项的形状为(batchsize，length，channel，height*width，patch)，patch表示图块数，即表示为
Figure BDA0003526079590000082
根据自注意力机制，利用Q，K计算注意力权重score：
score＝softmax(similarity(Q，K)))，
其中函数similarity用于计算Q，K之间的匹配度mr，此处为Q，K先分别重塑，后各维度对应数值相乘
Figure BDA0003526079590000083
并以尺度因子scale对结果进行尺度控制，此时匹配度mr的形状为(batchsiz，，lengthq，lengthkv，chann，l，height*wid，h，patch)，计算公式为：
Figure BDA0003526079590000091
函数softmax则用于将匹配度mr在lengthkv维度转化为和为1的权重值score，表示为：
Figure BDA0003526079590000092
求得注意力权重score后，重塑V，使V与score进行各维度对应数值相乘，并沿lengthkv维度求和实现注意力加权求和，得到自注意力计算结果，以Q，K，V为输入的自注意力计算Self-Attention(Q，K，V))表示为：
Figure BDA0003526079590000093
自注意力计算结果需将图块还原与重塑，并再经一层深度可分离卷积模块SepConv处理。SepConv内网络串行排列，依次为一层逐通道卷积层Conv2dSC1和一层逐点卷积层Conv2dSC2，其中，Conv2dSC1的输入数据尺寸为(112，112)，卷积核大小为(3，3)，输入通道数为16，输出通道数为64，步长为1，补零为1，组数为16，Conv2dSC2的输入数据尺寸为(112，112)，卷积核大小为(1，1)，输入通道数为64，输出通道数为16，步长为1，补零为0。
整个自注意力模块SAB内各模块串行排列，依次为GroupNorm，SA和GeLU，另外，SAB的输出为GroupNorm的输出与GeLU的输出之和，即两者之间存在一个残差连接。
步骤2-1-2-2：ESP编码器ESP-En内部各模块串行拼接，依次包括一个时间嵌入层TE(Time Embedding)和一个多层自注意力网络ESP-En-M；
时间嵌入层TE不含可训练参数，时间嵌入te基于正、余弦函数生成，公式如下：
Figure BDA0003526079590000094
其中，t表示第t时刻，i表示第i通道。由该公式生成的时间嵌入te，其形状为(batchsize，length，channel)，即
Figure BDA0003526079590000101
还需要将时间嵌入te在高宽维度复制扩展与重塑，使其形状变为(batchsize*length，channel，height,width)即
Figure BDA0003526079590000102
时间嵌入层TE的输入为步骤2-1-1中编码器En的输出，而时间嵌入层TE的输出为时间嵌入层TE的输入与te之和，并作为多层自注意力网络ESP-En-M的输入，其形状为(batchsize*Tinput，16，112，112)；
多层自注意力网络ESP-En-M由numen_layer个SA模块级联拼接而成，其输入数据只有来自时间嵌入层TE的输出数据而且同时作为网络的查询项源与键值项源，其输出数据的形状为(batchsize*Tinput，16，112,112)；
步骤2-1-2-3：ESP解码器ESP-De内部各模块串行拼接，依次包括一个时间嵌入层TE和一个多层自注意力网络ESP-De-M。
ESP解码器ESP-De的时间嵌入层TE结构与ESP-En的TE一致，ESP解码器ESP-De的输入数据是形状为(batchsize*Ttarget，16，112,112)的全0张量。
多层自注意力网络ESP-De-M由numde_layer个SA模块级联拼接而成，但与ESP-En-M的不同之处在于：第一，numde_layer必须为偶数；第二，在满足第一个条件下，numde_layer个SA模块前后两两成对，在一个模块对中，较低层SA的查询项源与键值项源均来自比该层低一层的模块，数据形状为(batchsize*Ttarget，16，112，112)，而较高层SA的查询项源与键值项源分别为较低层SA的输出和ESP-En的输出，对应的数据形状分别为(batchsize*Ttarget，16，112,112)和(batchsizn*Tinput，16，112,112)；
步骤2-1-3包括以下步骤：
步骤2-1-3-1：解码器De内各模块串联，共有三个卷积模块与两层反池化层，依次为Conv2dDe，1，MaxUnpool2dDe，1，Conv2dDe，2，MaxUnpool2dDe，2，Conv2dDe，3。其中，所有卷积模块均依次包含一层标准化层GroupNorm，一层卷积层，而Conv2dDe，1与Conv2dDe，2各额外包含一层非线性激活函数层LeakyReLU，并排在卷积层之后。Conv2dDe，1的输入数据尺寸为(112，112)，GroupNorm的组数为1，卷积层的卷积核大小为(3，3)，输入通道数为16，输出通道数为8，步长为1，补零为1，LeakyReLU的α参数为0.02；MaxUnpool2dDe，1为核大小为2、步长为2、补零为0、能够将特征图尺寸分别放大至2倍的最大反池化层；Conv2dDe，2的输入数据尺寸为(224，224)，卷积层输入通道数为8，输出通道数为4，其他参数设置与Conv2dDe，1一致；MaxUnpool2dDe，2与MaxUnpool2dDe，1一致；Conv2dDe，3的输入数据尺寸为(448，448)，卷积层的卷积核大小为(1，1)，输入通道数为4，输出通道数为1，步长为1，补零为0。此时输出数据的形状为4维(batchsizeTtarget，1,448,448)，整个自注意力时空网络SAB-STN的生成结果需要经重塑操作使其形状转为(batchsize，Ttarget，448,448，1)，随后以截断操作将生成结果限制在区间[-1，1]。
步骤3包括以下步骤：
步骤3-1：读取测试样本：将步骤1生成的测试集输入到步骤2生成的神经网络模型中；
步骤3-2：神经网络模型前向传播，将测试集输入神经网络模型后，生成外推结果。
本发明具有如下有益效果：
(1)数据预处理过程实现把数据中代表不同降水过程的回波序列区分出来，并进一步将回波序列做类别均衡采样，使得神经网络模型有效识别不同降水天气并学习出不同降水系统下的回波运动特征，从而细化回波外推结果；
(2)神经网络模型的ESP模块具有自注意力机制，通过自注意力机制，ESP能从一个回波序列整个时间尺度学习回波运动的时间相关性，同时ESP依靠二维卷积来获取每帧回波图像内的局部空间相关性，而注意力机制与二维卷积的共同作用，使得神经网络模型具备学习全时域下局部二维空间内回波运动特征的能力。
附图说明
下面结合附图和具体实施方式对本发明做更进一步的具体说明，本发明的上述和/或其他方面的优点将会变得更加清楚。
图1为本发明流程图。
图2为网络架构图。
图3为编码器构造图。
图4为自注意力模块结构图。
图5为回波序列处理器的编码器构造图。
图6为回波序列处理器的解码器构造图。
图7为解码器构造图。
图8是一幅回波图像样例的降噪可视化效果图。
图9是一个针对混合云降水的回波事件的预测图。
具体实施方式
如图1所示，本发明公开了本发明提供了一种基于自注意力时空网络的雷达回波外推方法，并以外推2021年上海地区雷达回波为例演示该方法，包括以下步骤：
步骤1，数据预处理：从国家气象数据网收集2017～2021年上海市的PPI组合反射率作为数据集，对其进行数据预处理，将每个定长回波序列划分为历史段和未来段，得到包含2017～2020年数据的训练集和包含2021年数据的测试集；
步骤2，训练网络：将训练集中的回波序列数据输入到SAB-STN网络中进行训练，训练完成后保存网络参数；
步骤3，测试网络：将测试集中的历史回波序列数据输入到步骤2训练好的神经网络模型中进行网络测试，输出从历史回波向未来外推的预测结果。
步骤1包括以下步骤：
步骤1-1：回波图像预处理：对数据集中所有回波图像进行预处理；
步骤1-2：回波序列预处理：根据回波图像的预处理结果，对所有回波序列进行预处理，包括序列类别标签标注、序列分类、序列类别均衡采样与序列的历史未来段划分；
步骤1-3：划分数据集：将预处理后的回波序列按年份将2021年的数据划分为测试集，其余部分化为训练集，其中测试集包括1405个序列，训练集包括40000个序列。
步骤1-1包括以下步骤：
步骤1-1-1：图像降噪。去除回波图像中可疑非降水回波，包括以下步骤：
步骤1-1-1-1：将原始单通道(通道数channel＝1)灰度值回波图像xraw的像素值pix转化为回波强度值I，计算公式如下：
Figure BDA0003526079590000131
步骤1-1-1-2，根据数据来源设置降水回波强度最低阈值Ith，根据Ith与步骤1-1-1-1得到的回波强度值I对回波图像进行一次降噪，一次降噪规定如下：对于回波强度小于最低阈值Ith的像素和缺测像素(NaN)，其值归零，其他情况则像素值不变，公式如下：
Figure BDA0003526079590000132
一次降噪完成之后，得到去除较低回波强度的单通道回波图像xraw；
步骤1-1-1-3：对一次降噪所得回波图像xraw进行连通域分析，统计得到每幅图像中各连通域面积ccA，连通域面积ccA，取连通域面积最低阈值ccAth＝13，根据ccAth与ccA对回波图像X进行二次降噪，二次降噪规定如下：对于面积小于ccAth的连通域，连通域内像素值归零，反之则像素值不变，公式如下：
Figure BDA0003526079590000133
二次降噪完成之后，得到不含较小面积杂波的单通道回波图像x，一幅回波图像样例的降噪可视化效果如图8所示；
步骤1-1-2：图像分类，包括以下步骤：
步骤1-1-2-1：统计每幅降噪后回波图像x中像素值不为0的像素个数numpix＞0，计算其占整个观测域像素个数numob的百分比，即有回波区与可观测区的面积比P1，公式如下：
Figure BDA0003526079590000134
根据P1是否为0首先将x二分类，其中，
Figure BDA0003526079590000135
表示降噪后无回波的回波图像，
Figure BDA0003526079590000136
表示降噪后有回波的回波图像，随后设置面积比区间集合
Figure BDA0003526079590000137
(集合内数值为百分数)：
Figure BDA0003526079590000141
并根据
Figure BDA0003526079590000142
的P1所属面积比区间，将
Figure BDA0003526079590000143
划分为4类；
步骤1-1-2-2：设置回波强度等级集合ΘI(集合内数值单位为dBZ)：
ΘI＝{(0，20]，(20，30]，(30，40]，(40，70]}，
统计每幅
Figure BDA0003526079590000144
中像素值在各回波等级的像素个数
Figure BDA0003526079590000145
计算其与numpix＞0之间的百分比，即某回波等级区与有回波区的面积比Pi，公式如下：
Figure BDA0003526079590000146
根据按序排列后最后不为0的Pi所对应的i，将
Figure BDA0003526079590000147
划分为4类，随后结合步骤1-1-2-1中对
Figure BDA0003526079590000148
划分的4类，将
Figure BDA0003526079590000149
最终划分为16类，加上
Figure BDA00035260795900001410
则x共分为17类，设置各类别对应标签如下表1：
表1
Figure BDA00035260795900001411
另外
Figure BDA00035260795900001412
所对应标签为(0，0)；
步骤1-1-3：图像标准化，标准化公式如下：
Figure BDA00035260795900001413
其中，pix是原始像素值pixnorm是标准化像素值。
步骤1-2包括以下步骤：
步骤1-2-1：序列类别标签标注，制作对应回波序列的回波类别标签序列：
根据步骤1-1-2所得回波图像类别标签，将一个回波序列X中所有图像的类别标签按回波序列顺序排列，得到回波类别标签序列Xcl；
步骤1-2-2：序列分类，利用K-Means聚类算法，分别对回波类别标签序列Xcl的两个类别维度进行序列聚类，在两个类别维度上得到各自回波序列聚类结果，随后综合这两类聚类结果得到一个完整的聚类结果。其中，在K-Means聚类算法中所使用的序列样本距离为欧氏距离。在所用40000个训练数据中，无监督聚类结果表示如下：强度在30dBZ以下的有7类，包含12117个样本；强度在30～45dBZ的有6类，包含17250个样本；强度在45dBZ以上的有6类，包含10663个样本；
步骤1-2-3：序列类别均衡采样，对序列聚类结果进行类别均衡采样，回波序列采样公式如下：
Figure BDA0003526079590000151
其中，N为回波序列数据集的类别数，nj为类别j的样本总数pi为从类别i中采样一个回波序列的概率；
步骤1-2-4：序列的历史未来段划分，设置历史帧数与未来帧数分别作为SAB-STN神经网络模型的输入长度Tinput＝10与外推目标长度Ttarget＝20，Tinput+Ttarget表示整个序列的长度为30，并以此将回波序列X划分为历史输入段Xinput与未来外推段Xtarget前后两部分：
Figure BDA0003526079590000152
其中xi表示序列中第i帧图像(i＝1，2，...，Tinput+Ttarget)；
随后将回波图像大小重塑(高度height，宽度width)至(448，448)，此时一个回波序列的数据形状(length，channel，height，width)为(Tinput，1，448，448)，表示为
Figure BDA0003526079590000153
即回波序列X属于由长度L，通道数C，高度H，宽度W四个维度构成的高维实数空间。
步骤1-2-2包括以下步骤：
步骤1-2-2-1：K-Means聚类中心选择，设置类别数k＝5，选择初始化的k个回波类别标签序列样本Xcl作为初始聚类中心cs：
cs＝cs1，cs2，...，csk；
其中cs1指第1个聚类中心，其他依此类推。
步骤1-2-2-2：归类每个样本。针对每个样本Xcl计算它到k个聚类中心的欧式距离dist并将其分到距离最小的聚类中心所对应的类中：
Figure BDA0003526079590000161
其中arg min函数指返回使被作用函数达到最小值的变量取值，c(i)表示第i个样例Xcl，i与k个类中距离最近的那个类心csj，其中j＝1，2，..k，dist计算公式如下：
Figure BDA0003526079590000162
其中，a，b表示长度同为length的不同回波类别标签序列，ai，bi分别表示a，b在序列位置i处的回波标签向量；
步骤1-2-2-3：更新类的质心。针对每个类别，重新计算属于该类所有样本的质心即聚类中心：
Figure BDA0003526079590000163
其中ci表示csi所对应的类(i＝1，2，..k)；
步骤1-2-2-4：获得聚类结果。重复步骤1-2-2-3和步骤1-2-2-4，直到收敛。此时确定k个类中各自包含的样本。
步骤2包括以下步骤：
步骤2-1：神经网络模型建立：如图2所示，搭建SAB-STN模型各模块及初始化模型参数；
步骤2-2：训练超参数初始化：采用Adam优化器训练网络，设置学习率lr＝0.001、每次输入一组样本所含的样本数batchsize＝4和训练集遍历数epoch＝20；
步骤2-3：网络参数训练及更新：采用平均平方误差MSE作为损失函数，每次取一个batchsize的数据输入网络，以训练和优化网络参数，通过反向传播算法更新参数，模型参数初始化采用模型预训练和随机初始化相结合的方法，MSE损失函数公式如下：
Figure BDA0003526079590000171
其中pixi表示一幅真实回波图像内第i个像素，
Figure BDA0003526079590000172
表示该真实回波图像对应生成回波图像内第i个像素，448×448表示图像内像素点个数。
步骤2-1包括以下步骤：
步骤2-1-1：编码器En构造：编码器为卷积网络，用于将原始回波图像转化为较高级的图像表征；
步骤2-1-2：回波序列处理器ESP(Echo Series Processor)构造：序列处理器包括ESP编码器ESP-En和ESP解码器ESP-De两部分，自注意力模块SAB为基本单元，用于学习一个回波序列内由编码器生成的各时刻图像表征之间的时间关联，以及各时刻图像表征自身的空间关联；
步骤2-1-3：解码器De构造：解码器为卷积网络，用于将序列处理器所生成的特征转化成外推结果。
步骤2-1-4：初始化模型参数：对于SAB-STN神经网络模型中的卷积核，其权重采用kaiming正态分布初始化方法，其偏执采用0值初始化方法；对于模型中的标准化层，其权重采用1值初始化方法，其偏执采用0值初始化方法。
步骤2-1-1包括以下步骤：
步骤2-1-1-1：编码器En包括三个卷积模块Conv2dEn，1，Conv2EEn，2，Conv2EEn，3，两层最大池化层MaxPool2dEn，1，MaxPool2dEn，2。
如图3所示，编码器En各模块以串行方式排列，顺序依次为Conv2dEn，1，MaxPool2dEn，1，Conv2dEn，2，MaxPool2dEn，2，Conv2dEn，3，其中，所有卷积模块均依次包含一层卷积层，一层非线性激活函数层LeakyReLU，而Conv2dEn，2与Conv2dEn，3各额外包含一层标准化层GroupNorm，并排在卷积层之前。Conv2dEn，1的输入数据尺寸为(448，448)，卷积层的卷积核大小为(3，3)，输入通道数为1，输出通道数为4，步长为1，补零为1，LeakyReLU的α参数为0.02，GroupNorm的组数为1；MaxPool2dEn，1为核大小为2、步长为2、补零为0、可将特征图尺寸分别缩小1/2的最大池化层；Conv2dEn，2的输入数据尺寸为(224，224)，卷积层输入通道数为4，输出通道数为8，其他参数设置与Conv2dEn，1一致；MaxPool2dEn，2与MaxPool2dEn，1一致；Conv2dEn，3的输入数据尺寸为(112，112)，卷积层输入通道数为8，输出通道数为16，其他参数设置与Conv2dEn，1一致。此外，输入数据在进入Conv2dEn，1前由5维形状(batchsize＝4，length＝10，channel＝16，height＝11，，width＝112)重塑至4维(batchsize*length＝4，，channel＝1，，height＝11，，width＝112)，表示为
Figure BDA0003526079590000181
此时形状为(batchsize，length，channel，height，width)＝(4，10，1，448，448)的输入数据经En处理，得到形状为(40，16，112，112)的回波图像表征。
步骤2-1-2包括以下步骤：
步骤2-1-2-1：如图4所示，自注意力模块SAB包括一层自注意力层SA，一层非线性激活函数层GeLU，以及一层标准化层GroupNorm，其中SAB的核心在于自注意力层SA；
在SA层中，输入端同时传入查询项源与键值项源，查询项源通过查询卷积Conv2dQ生成查询项Q，键值项源通过键值卷积Conv2dKV生成键值项KV，Conv2dQ，Conv2dKV均为逐通道channel-wise二维卷积层，其中Conv2dQ的输入数据尺寸为(112，112)，卷积核大小为(3，3)，输入通道数与输出通道数相等且为16，步长为1，补零为1，组数为16，Conv2dKV的输入数据尺寸为(112，112)，卷积核大小为(3，3)，输入通通道数为16，输出通道数为32，步长为1，补零为1，组数为16，图块边长值regionsize＝14，Q与KV在图像高宽维度分别切分成若干不重叠且边长值为regionsize的图块，并且KV沿channel维度均分为两份得到键项K和值项V，此时Q，K，V三项的形状为(batchsize，length，channel，height*width，patch)而patch表示图块数，即表示为
Figure BDA0003526079590000191
根据自注意力机制，利用Q，K计算注意力权重score，即：
score＝softmax(similarity(Q，K))，
其中函数similarity用于计算Q，K之间的匹配度mr，此处为Q，K先分别重塑，后各维度对应数值相乘
Figure BDA0003526079590000192
并以尺度因子scale对结果进行尺度控制，此时匹配度mr的形状为(batchsiz，，lengthq，lengthkv，chann，l，height*width，patch)，计算公式为：
Figure BDA0003526079590000193
函数softmax则用于将匹配度mr在lengthkv维度转化为和为1的权重值score，表示为：
Figure BDA0003526079590000194
求得注意力权重score后，重塑V使之与score进行各维度对应数值相乘，并沿lengthkv维度求和实现注意力加权求和，得到自注意力计算结果，以Q，K，V为输入的自注意力计算Self-Attention(Q，K，V)表示为：
Figure BDA0003526079590000195
自注意力计算结果需将图块还原与重塑，并再经一层深度可分离卷积模块SepConv处理。SepConv内网络串行排列，依次为一层逐通道卷积层Conv2dSC1和一层逐点卷积层Conv2dSC2，其中，Conv2dSC1的输入数据尺寸为(112，112)，卷积核大小为(3，3)，输入通道数为16，输出通道数为64，步长为1，补零为1，组数为16，Conv2dSC2的输入数据尺寸为(112，112)，卷积核大小为(1，1)，输入通道数为64，输出通道数为16，步长为1，补零为0。
整个SAB内各模块串行排列，依次为GroupNorm，SA和GeLU，另外，SAB的输出为GroupNorm的输出与GeLU的输出之和，即两者之间存在一个残差连接。
步骤2-1-2-2：如图5所示，ESP编码器ESP-En内部各模块串行拼接，依次包括一个时间嵌入层TE(Time Embedding)和一个多层自注意力网络ESP-En-M；
时间嵌入层TE不含可训练参数，时间嵌入te基于正、余弦函数生成，公式如下：
Figure BDA0003526079590000201
其中，t表示第t时刻，i表示第i通道。由该公式生成的时间嵌入te，其形状为(batchsize＝4，length＝10，channel＝16)，即
Figure BDA0003526079590000202
还需要将其在高宽维度复制扩展与重塑，使其形状变为(batchsize*length＝40,channel＝16，height＝112，width＝112)，即
Figure BDA0003526079590000203
TE的输入为步骤2-1-1中En的输出，而TE的输出为TE的输入与te之和，并作为ESP-En-M的输入，其形状为(40，16，112，112)；
多层自注意力网络ESP-En-M由numen_layer＝4个SA模块级联拼接而成，其输入数据只有来自TE的输出数据而且同时作为网络的查询项源与键值项源，其输出数据的形状为(40，16，112，112)；
步骤2-1-2-3：如图6所示，ESP解码器ESP-De内部各模块串行拼接，依次包括一个时间嵌入层TE和一个多层自注意力网络ESP-De-M。
ESP-De的TE结构与ESP-En的TE一致，但就输入数据而言，前者的数据是形状为(80，16，112，112)的全0张量。
多层自注意力网络ESP-De-M由numde_layer＝4个SA模块级联拼接而成，但与ESP-En-M的不同之处在于：第一，numde_layer必须为偶数；第二，在满足第一个条件下，numde_layer个SA模块前后两两成对，在一个模块对中，较低层SA的查询项源与键值项源均来自比该层低一层的模块，数据形状为(80，16，112，112)，而较高层SA的查询项源与键值项源分别为较低层SA的输出和ESP-En的输出，对应的数据形状分别为(80，16，112，112)和(80，16，112，112)；
步骤2-1-3包括以下步骤：
步骤2-1-3-1：如图7所示，解码器De内各模块串联，共有三个卷积模块与两层反池化层，依次为Conv2dDe，1，MaxUnpool2dDe，1，Conv2dDe，2，MaxUnpool2dDe，2，Conv2dDe，3。其中，所有卷积模块均依次包含一层标准化层GroupNorm，一层卷积层，而Conv2dDe，1与Conv2dDe，2各额外包含一层非线性激活函数层LeakyReLU，并排在卷积层之后。Conv2dDe，1的输入数据尺寸为(112，112)，GroupNorm的组数为1，卷积层的卷积核大小为(3，3)，输入通道数为16，输出通道数为8，步长为1，补零为1，LeakyReLU的α参数为0.02；MaxUnpool2dDe，1为核大小为2、步长为2、补零为0、可将特征图尺寸分别放大至2倍的最大反池化层；Conv2dDe，2的输入数据尺寸为(224，224)，卷积层输入通道数为8，输出通道数为4，其他参数设置与Conv2dDe，1一致；MaxUnpool2dDe，2与MaxUnpool2dDe，1一致；Conv2dDe，3的输入数据尺寸为(448，448)，卷积层的卷积核大小为(1，1)，输入通道数为4，输出通道数为1，步长为1，补零为0。此时输出数据的形状为4维(80，1，448，448)，整个自注意力时空网络SAB-STN的生成结果需要经重塑操作使其形状转为(4，20，448，448，1)，随后以截断操作将生成结果限制在区间[-1，1]。
步骤3包括以下步骤：
步骤3-1：读取测试样本：将步骤1生成的测试集输入到步骤2生成的神经网络模型中；
步骤3-2：神经网络模型前向传播，将测试集输入神经网络模型后，生成外推结果。以一个针对混合云降水的回波事件的预测为例，其结果如图9所示，其中前10帧为输入数据，后20帧为外推数据；
雷达回波外推中最常用的数据是2小时内单通道PPI组合反射率，图上像素格点表示现实地理水平空间按一定分辨率的投影，像素值则表示回波强度的映射；通过一个回波序列，能较为精确地观察到局部回波所代表天气现象的运动变化过程及可能的运动趋势。如果收集到足量某类天气现象全过程的回波序列，则可以利用统计方法先识别出该天气现象的早期阶段，然后以此段回波序列作为外推依据，分析出潜在的运动变化，进而得到未来时间内可能的回波状态；当经过规定时间间隔，需要获取新的回波序列以判断当前天气的变化阶段，再做外推。研究表明深度神经网络具有优秀的数据拟合能力，使之作为用于外推的统计方法，具体来说，深度神经网络不仅能拟合出潜在空间内回波基本的水平移动以及回波强度变化，还能拟合出不同天气现象所独有的回波运动，换句话说，深度神经网络能根据已知回波运动过程识别出该过程最有可能属于什么天气现象的什么变化阶段，再进行外推预测。
当前已出现基于CNN、RNN和cGAN等神经网络的回波外推方法，但这些神经网络在原理上缺乏从全时域学习回波序列时空相关性的能力。与其他网络相比，自注意力网络的能从整个回波序列时间段学习回波图像上任何一块水平空间的回波运动特征，因此本发明围绕自注意力机制，同时结合回波序列数据特点，提出了适用于回波外推的SAB-STN网络。
本发明提供了一种基于自注意力时空神经网络模型的雷达回波外推方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (10)

1.一种基于自注意力时空神经网络模型的雷达回波外推方法，其特征在于，包括以下步骤：
步骤1，数据预处理：对经相关系数质量控制的PPI组合反射率雷达回波产品进行数据预处理，将每个定长回波序列划分为历史段和未来段，得到训练集和测试集；
步骤2，训练神经网络模型：将训练集中的回波序列数据输入到SAB-STN神经网络模型中进行训练，训练完成后保存神经网络模型参数；
步骤3，测试神经网络模型：将测试集中的历史回波序列数据输入到步骤2训练好的神经网络模型中进行神经网络模型测试，输出从历史回波向未来外推的预测结果。
2.根据权利要求1所述的方法，其特征在于，步骤1包括以下步骤：
步骤1-1：回波图像预处理：对数据集中所有回波图像进行预处理；
步骤1-2：回波序列预处理；
步骤1-3：划分数据集：将预处理后的回波序列划分为测试集和训练集。
3.根据权利要求2所述的方法，其特征在于，步骤1-1包括以下步骤：
步骤1-1-1：图像降噪：去除回波图像中可疑非降水回波，包括以下步骤：
步骤1-1-1-1：将原始单通道灰度值回波图像xraw的像素值pix转化为回波强度值I，计算公式如下：
Figure FDA0003526079580000011
步骤1-1-1-2，根据数据来源设置降水回波强度最低阈值Ith，根据最低阈值Ith与步骤1-1-1-1得到的回波强度值I对回波图像进行一次降噪，一次降噪规定如下：对于回波强度小于最低阈值Ith的像素和缺测像素NaN，将它们的值归零，其他情况则像素值不变，公式如下：
Figure FDA0003526079580000012
一次降噪完成之后，得到去除较低回波强度的单通道回波图像xraw；
步骤1-1-1-3：对一次降噪所得回波图像xraw进行连通域分析，统计得到每幅图像中各连通域面积ccA，设定连通域面积最低阈值ccAth，根据ccAth与ccA对回波图像X进行二次降噪，二次降噪规定如下：对于面积小于ccAth的连通域，连通域内像素值归零，否则像素值不变，公式如下：
Figure FDA0003526079580000021
二次降噪完成之后，得到不含较小面积杂波的单通道回波图像x；
步骤1-1-2：图像分类，包括以下步骤：
步骤1-1-2-1：统计每幅降噪后回波图像x中像素值不为0的像素个数numpix＞0，计算其占整个观测域像素个数numob的百分比，即有回波区与可观测区的面积比P1，公式如下：
Figure FDA0003526079580000022
根据P1是否为0首先将x二分类，其中，
Figure FDA0003526079580000023
表示降噪后无回波的回波图像，
Figure FDA0003526079580000024
表示降噪后有回波的回波图像，设置面积比区间集合
Figure FDA0003526079580000025
Figure FDA0003526079580000026
并根据
Figure FDA0003526079580000027
的P1所属面积比区间，将
Figure FDA0003526079580000028
划分为4类；
步骤1-1-2-2：设置回波强度等级集合ΘI：
ΘI＝{(0，20]，(20，30]，(30，40]，(40，70]}，
统计每幅
Figure FDA0003526079580000029
中像素值在各回波等级的像素个数
Figure FDA00035260795800000215
计算与numpix＞0之间的百分比，即回波等级区与有回波区的面积比Pi，公式如下：
Figure FDA00035260795800000210
根据按序排列后最后不为0的Pi所对应的i，将
Figure FDA00035260795800000211
划分为4类，结合步骤1-1-2-1中对
Figure FDA00035260795800000212
划分的4类，将
Figure FDA00035260795800000213
最终划分为16类，加上
Figure FDA00035260795800000214
则x共分为17类；
Figure FDA0003526079580000031
所对应标签为(0,0)；
步骤1-1-3：进行图像标准化，公式如下：
Figure FDA0003526079580000032
其中，pix是原始像素值，pixnorm是标准化像素值。
4.根据权利要求3所述的方法，其特征在于，步骤1-2包括以下步骤：
步骤1-2-1：标注序列类别标签，制作对应回波序列的回波类别标签序列：
根据步骤1-1-2所得回波图像类别标签，将一个回波序列X中所有图像的类别标签按回波序列顺序排列，得到回波类别标签序列Xcl；
步骤1-2-2：利用K-Means聚类算法，分别对回波类别标签序列Xcl的两个类别维度进行序列聚类，在两个类别维度上得到各自回波序列聚类结果，随后综合这两类聚类结果得到一个完整的聚类结果，其中，在K-Means聚类算法中所使用的序列样本距离为欧氏距离；
步骤1-2-3：对序列聚类结果进行类别均衡采样，回波序列采样公式如下：
Figure FDA0003526079580000033
其中，N为回波序列数据集的类别数，nj为类别j的样本总数，pi为从类别i中采样一个回波序列的概率；
步骤1-2-4：划分序列的历史未来段，设置历史帧数与未来帧数分别作为SAB-STN神经网络模型的输入长度Tinput与外推目标长度Ttarget，则Tinput+Ttarget表示整个序列的长度，将回波序列X划分为历史输入段Xinput与未来外推段Xtarget前后两部分：
Figure FDA0003526079580000034
其中xi表示序列中第i帧图像；
将回波图像大小重塑：(height，width)重塑至(448，448)，height表示高度，width表示宽度，此时一个回波序列的数据形状(length，channel，height，width)为(Tinput，1，448，448)，表示为回波序列
Figure FDA0003526079580000041
Figure FDA0003526079580000042
表示长度L，通道数C，高度H，宽度W四个维度构成的高维实数空间，length表示长度，channel表示通道数。
5.根据权利要求4所述的方法，其特征在于，步骤1-2-2包括以下步骤：
步骤1-2-2-1：选择K-Means聚类中心，设置类别数k，选择初始化的k个回波类别标签序列样本Xcl作为初始聚类中心cs：
cs＝cs1，cs2，…，csk；
其中csk指第k个聚类中心；
步骤1-2-2-2：归类每个样本：针对每个样本Xcl计算它到k个聚类中心的欧式距离dist并将其分到距离最小的聚类中心所对应的类中：
Figure FDA0003526079580000043
其中arg min函数指返回使被作用函数达到最小值的变量取值，c(i)表示第i个样例Xcl，i与k个类中距离最近的那个聚类中心csj，其中j＝1,2,..k，函数dist计算公式如下：
Figure FDA0003526079580000044
其中，a，b表示长度同为length的不同回波类别标签序列，ai，bi分别表示a，b在序列位置i处的回波标签向量；
步骤1-2-2-3：更新类的质心：针对每个类别，重新计算属于该类别所有样本的质心即聚类中心：
Figure FDA0003526079580000045
其中ci表示csi所对应的类，i＝1,2,..k；
步骤1-2-2-4：获得聚类结果：重复步骤1-2-2-3，直到收敛，此时确定k个类中各自包含的样本。
6.根据权利要求5所述的方法，其特征在于，步骤2包括以下步骤：
步骤2-1：建立神经网络模型：搭建神经网络模型各模块及初始化模型参数；
步骤2-2：初始化训练超参数：采用Adam优化器训练神经网络模型，设置学习率lr、每次输入一组样本所含的样本数batchsize和训练集遍历数epoch；
步骤2-3：神经网络模型参数训练及更新：采用平均平方误差MSE作为损失函数，每次取一个batchsize的数据输入神经网络模型，以训练和优化神经网络模型参数，通过反向传播算法更新参数，模型参数初始化采用模型预训练和随机初始化相结合的方法，损失函数MSE计算公式如下：
Figure FDA0003526079580000051
其中pixi表示一幅真实回波图像内第i个像素，
Figure FDA0003526079580000052
表示真实回波图像对应生成回波图像内第i个像素，height×width表示图像内像素点个数。
7.根据权利要求6所述的方法，其特征在于，步骤2-1包括以下步骤：
步骤2-1-1：构造编码器En：编码器为卷积神经网络模型，用于将原始回波图像转化为抽象图像表征；
步骤2-1-2：构造回波序列处理器ESP：序列处理器包括ESP编码器ESP-En和ESP解码器ESP-De两部分，以自注意力模块SAB为基本单元，用于学习一个回波序列内由编码器生成的各时刻图像表征之间的时间关联，以及各时刻图像表征自身的空间关联；
步骤2-1-3：构造解码器De：解码器为卷积神经网络模型，用于将序列处理器所生成的特征转化成外推结果；
步骤2-1-4：初始化模型参数：对于SAB-STN神经网络模型中的卷积核，权重采用kaiming正态分布初始化方法，偏执采用0值初始化方法；对于模型中的标准化层，权重采用1值初始化方法，偏执采用0值初始化方法。
8.根据权利要求7所述的方法，其特征在于，步骤2-1-1包括以下步骤：
步骤2-1-1-1：编码器En包括三个卷积模块Conv2dEn，1，Conv2dEn，2，Conv2dEn，3，两层最大池化层MaxPool2dEn，1，MaxPool2dEn，2；
各模块以串行方式排列，顺序依次为Conv2dEn，1，MaxPool2dEn，1，Conv2dEn，2，MaxPool2dEn，2，Conv2dEn，3，其中，所有卷积模块均依次包含一层卷积层，一层非线性激活函数层LeakyReLU，而Conv2dEn，3与Conv2dEn，3各额外包含一层标准化层GroupNorm，并排在卷积层之前；Conv2dEn，1的输入数据尺寸为(448,448)，卷积层的卷积核大小为(3,3)，输入通道数为1，输出通道数为4，步长为1，补零为1，LeakyReLU的α参数为0.02，GroupNorm的组数为1；MaxPool2dEn，1为核大小为2、步长为2、补零为0、能够将特征图尺寸分别缩小1/2的最大池化层；Conv2dEn，2的输入数据尺寸为(224,224)，卷积层输入通道数为4，输出通道数为8；MaxPool2dEn，2与MaxPool2dEn，1一致；Conv2dEn，3的输入数据尺寸为(112,112)，卷积层输入通道数为8，输出通道数为16；输入数据在进入Conv2dEn，1前由5维形状(batchsize,length,channel,height,width)重塑至4维(batchsize*length,channel,height,width)，表示为
Figure FDA0003526079580000061
此时形状为(batchsize,length,channel,height,width)＝(batchsize,Tinput,1,448,448)的输入数据经编码器En处理，得到形状为(batchsize*Tinput,16,112,112)的回波图像表征。
9.根据权利要求8所述的方法，其特征在于，步骤2-1-2包括以下步骤：
步骤2-1-2-1：自注意力模块SAB包括一层自注意力层SA，一层非线性激活函数层GeLU，以及一层标准化层GroupNorm；
在自注意力层SA中，输入端同时传入查询项源与键值项源，查询项源通过查询卷积Conv2dQ生成查询项Q，键值项源通过键值卷积Conv2dKV生成键值项KV，Conv2dQ，Conv2dKV均为逐通道channel-wise二维卷积层，其中Conv2dQ的输入数据尺寸为(112,112)，卷积核大小为(3,3)，输入通道数与输出通道数相等且为16，步长为1，补零为1，组数为16，Con2dKV的输入数据尺寸为(112,112)，卷积核大小为(3,3)，输入通通道数为16，输出通道数为32，步长为1，补零为1，组数为16，图块边长值regionsize，Q与KV在图像高宽维度分别切分成不重叠且边长值为regionsize的图块，并且KV沿channel维度均分为两份得到键项K和值项V，此时Q，K，V三项的形状为(batchsize,length,channel,height*width,patch)，patch表示图块数，即表示为
Figure FDA0003526079580000071
根据自注意力机制，利用Q，K计算注意力权重score：
score＝softmax(similarity(Q，K))，
其中函数similarity用于计算Q，K之间的匹配度mr，此处为Q，K先分别重塑，后各维度对应数值相乘，并以尺度因子scale对结果进行尺度控制，此时匹配度mr的形状为(batchsize,lengthq,lengthkv,channel,height*width,patch)，计算公式为：
Figure FDA0003526079580000072
函数softmax用于将匹配度mr在lengthkv维度转化为和为1的权重值score，表示为：
Figure FDA0003526079580000073
求得注意力权重score后，重塑V，使V与score进行各维度对应数值相乘，并沿lengthkv维度求和实现注意力加权求和，得到自注意力计算结果，以Q，K，V为输入的自注意力计算Self-Attention(Q，K，V)表示为：
Figure FDA0003526079580000074
自注意力计算结果需将图块还原与重塑，并再经一层深度可分离卷积模块SepConv处理；SepConv内网络串行排列，依次为一层逐通道卷积层Conv2dSC1和一层逐点卷积层Conv2dSC2，其中，Conv2dSC1的输入数据尺寸为(112,112)，卷积核大小为(3,3)，输入通道数为16，输出通道数为64，步长为1，补零为1，组数为16，Conv2dSC2的输入数据尺寸为(112,112)，卷积核大小为(1,1)，输入通道数为64，输出通道数为16，步长为1，补零为0；
整个自注意力模块SAB内各模块串行排列，依次为GroupNorm，SA和GeLU，另外，SAB的输出为GroupNorm的输出与GeLU的输出之和；
步骤2-1-2-2：ESP编码器ESP-En内部各模块串行拼接，依次包括一个时间嵌入层TE和一个多层自注意力网络ESP-En-M；
时间嵌入层TE不含可训练参数，时间嵌入te基于正、余弦函数生成，公式如下：
Figure FDA0003526079580000081
其中，t表示第t时刻，i表示第i通道；时间嵌入te形状为(batchsize,length,channel)，即
Figure FDA0003526079580000082
还需要将时间嵌入te在高宽维度复制扩展与重塑，使其形状变为(batchsize*length,channel,height,width)，即
Figure FDA0003526079580000083
时间嵌入层TE的输入为步骤2-1-1中编码器En的输出，而时间嵌入层TE的输出为时间嵌入层TE的输入与te之和，并作为多层自注意力网络ESP-En-M的输入，形状为(batchsize*Tinput,16,112,112)；
多层自注意力网络ESP-En-M由numen_layer个SA模块级联拼接而成，其输入数据只有来自时间嵌入层TE的输出数据而且同时作为网络的查询项源与键值项源，输出数据的形状为(batchsize*Tinput,16,112,112)；
步骤2-1-2-3：ESP解码器ESP-De内部各模块串行拼接，依次包括一个时间嵌入层TE和一个多层自注意力网络ESP-De-M；
ESP解码器ESP-De的时间嵌入层TE结构与ESP-En的TE一致，ESP解码器ESP-De的输入数据是形状为(batchsize*Ttarget,16,112,112)的全0张量；
多层自注意力网络ESP-De-M由numde_layer个SA模块级联拼接而成，numde_layer为偶数，并且numde_layer个SA模块前后两两成对，在一个模块对中，较低层SA的查询项源与键值项源均来自比该层低一层的模块，数据形状为(batchsize*Ttarget,16,112,112)，而较高层SA的查询项源与键值项源分别为较低层SA的输出和ESP-En的输出，对应的数据形状分别为(batchsize*Ttarget,16,112,112)和(batchsize*Tinput,16,112,112)。
10.根据权利要求9所述的方法，其特征在于，步骤2-1-3包括以下步骤：
步骤2-1-3-1：解码器De内各模块串联，共有三个卷积模块与两层反池化层，依次为Conv2dDe，1，MaxUnpool2dDe，1，Conv2dDe，2，MaxUnpool2dDe，2，Conv2dDe，3；其中，所有卷积模块均依次包含一层标准化层GroupNorm，一层卷积层，而Conv2dDe，1与Conv2dDe，2各额外包含一层非线性激活函数层LeakyReLU，并排在卷积层之后；Conv2dDe，1的输入数据尺寸为(112,112)，GroupNorm的组数为1，卷积层的卷积核大小为(3,3)，输入通道数为16，输出通道数为8，步长为1，补零为1，LeakyReLU的α参数为0.02；MaxUnpool2dDe，1为核大小为2、步长为2、补零为0、能够将特征图尺寸分别放大至2倍的最大反池化层；Conv2dDe，2的输入数据尺寸为(224,224)，卷积层输入通道数为8，输出通道数为4；MaxUnpool2dDe，2与MaxUnpool2dDe，1一致；Conv2dDe，3的输入数据尺寸为(448,448)，卷积层的卷积核大小为(1,1)，输入通道数为4，输出通道数为1，步长为1，补零为0；此时输出数据的形状为4维(batchsize*Ttarget,1,448,448)，整个自注意力时空网络SAB-STN的生成结果需要经重塑操作使其形状转为(batchsize,Ttarget,448,448,1)，随后以截断操作将生成结果限制在区间[-1,1]；
步骤3包括以下步骤：
步骤3-1：读取测试样本：将步骤1生成的测试集输入到步骤2生成的神经网络模型中；
步骤3-2：神经网络模型前向传播，将测试集输入神经网络模型后，生成外推结果。
CN202210193798.0A 2022-03-01 2022-03-01 一种基于自注意力时空神经网络模型的雷达回波外推方法 Pending CN114624715A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210193798.0A CN114624715A (zh) 2022-03-01 2022-03-01 一种基于自注意力时空神经网络模型的雷达回波外推方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210193798.0A CN114624715A (zh) 2022-03-01 2022-03-01 一种基于自注意力时空神经网络模型的雷达回波外推方法

Publications (1)

Publication Number Publication Date
CN114624715A true CN114624715A (zh) 2022-06-14

Family

ID=81899533

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210193798.0A Pending CN114624715A (zh) 2022-03-01 2022-03-01 一种基于自注意力时空神经网络模型的雷达回波外推方法

Country Status (1)

Country Link
CN (1) CN114624715A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116106856A (zh) * 2023-04-13 2023-05-12 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 雷暴大风的识别模型建立方法、识别方法及计算设备

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116106856A (zh) * 2023-04-13 2023-05-12 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 雷暴大风的识别模型建立方法、识别方法及计算设备
CN116106856B (zh) * 2023-04-13 2023-08-18 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 雷暴大风的识别模型建立方法、识别方法及计算设备

Similar Documents

Publication Publication Date Title
CN110210463B (zh) 基于Precise ROI-Faster R-CNN的雷达目标图像检测方法
CN110135267B (zh) 一种大场景sar图像细微目标检测方法
CN110414377B (zh) 一种基于尺度注意力网络的遥感图像场景分类方法
CN107016677B (zh) 一种基于fcn和cnn的云图分割方法
CN112036335B (zh) 一种反卷积引导的半监督植物叶部病害识别与分割方法
CN110033473B (zh) 基于模板匹配和深度分类网络的运动目标跟踪方法
CN109190491B (zh) 残差卷积神经网络sar影像海冰分类方法
CN108629370B (zh) 一种基于深度置信网络的分类识别算法及装置
CN105913081B (zh) 基于改进的PCAnet的SAR图像分类方法
CN108537102A (zh) 基于稀疏特征与条件随机场的高分辨sar图像分类方法
CN111611972B (zh) 基于多视图多任务集成学习的作物叶片种类识别方法
CN108345856B (zh) 基于异构卷积神经网络集成的sar自动目标识别方法
CN111273288B (zh) 一种基于长短期记忆网络的雷达未知目标识别方法
Afework et al. Detection of bacterial wilt on enset crop using deep learning approach
CN110853070A (zh) 基于显著性与Grabcut的水下海参图像分割方法
Panati et al. Feature relevance evaluation using grad-CAM, LIME and SHAP for deep learning SAR data classification
CN110490894B (zh) 基于改进的低秩稀疏分解的视频前背景分离方法
Pinto et al. Crop disease classification using texture analysis
Liu et al. Deep fuzzy graph convolutional networks for PolSAR imagery pixelwise classification
CN112749675A (zh) 一种基于卷积神经网络的马铃薯病害识别方法
CN114359631A (zh) 基于编码-译码弱监督网络模型的目标分类与定位方法
CN114581434A (zh) 基于深度学习分割模型的病理图像处理方法、电子设备
Periyasamy et al. How to get the most out of u-net for glacier calving front segmentation
CN114624715A (zh) 一种基于自注意力时空神经网络模型的雷达回波外推方法
CN107766792A (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