基于convLSTM模型的雷达图像外推算法

最新推荐文章于 2024-09-14 16:30:14 发布

nobrody

最新推荐文章于 2024-09-14 16:30:14 发布

阅读量5k

点赞数 10

CC 4.0 BY-SA版权

分类专栏：短临降水预报文章标签：深度学习短临预报

本文链接：https://blog.youkuaiyun.com/zshluckydogs/article/details/113701857

短临降水预报专栏收录该内容

4 篇文章

订阅专栏

本文探讨了基于convLSTM模型的雷达图像外推算法，旨在进行短时预报。通过将卷积网络与LSTM结合，提取时空特征。然而，实验结果显示模型在捕捉非平移运动时面临挑战，未来可能需要引入local-variant网络以提高预测准确性。作者提供了网络架构代码，并分享了模型训练与评估指标，包括acc、loss、psnr、ssim、POD和FAR。

前言

现在的数值预报模式是一种长时预报，预报未来一天或者一周的降水或者天气情况。这种预报密度已无法满足现今的生产生活需要，临近预报或者短时预报对于人们的出行，民航飞行计划制定，工农业生产有重要的指导意义。

云与雨水的关系

不去谈降水形成背后复杂的大气物理运动，我们从生活常识来理解云与降水的关系。我们常用万里晴空或者万里无云来描述晴天，常用阴云密布来形容雨天，显然云层的状态直接影响降雨。

云层的状态包含云层的高度、厚度、形状、运动趋势等。为了能够利用云层的状况或者说是特征来建立云与降水的关系，首先要获取云层的特征描述也叫观测数据。目前成熟的方法是用雷达获取云层的反射率数据，基于雷达反射率数据做降水分析等。

雷达基数据解析

关于多普勒雷达基数据解析，可以在网上找到很多专业的解释，如果对这一块不是很了解的可以查阅资料快速的理解一下，或者直接找经过解析的雷达反射率图像。简单解释一下，雷达基数据是以二进制形式存储的数据，一般是2位连在一起解析出实际的反射率数值。

可以参考我的实现（同样是参考前人的实现）

class Radar(object):
    ## path + file_name
    ## 
    def __init__(self,*args):
        path,name = args
        self.AbsolutePath = os.path.join(path,name)
        file = open(self.AbsolutePath,'rb')
        file.seek(0)
        self.RawData = np.array([int(i) for i in file.read()])
        self.Name = name[:-4]
        self.Count = int(len(self.RawData)/2432) ## 总个数
        self.RawArray = self.RawData.reshape(self.Count,2432)
        _,_,_,self.Station,self.time = self.Name.split('_')[:5]
        
        
        self.time_ms = read_data(self.AbsolutePath)
        logger.info(f'ms {self.time_ms}')
        #self.Julian = [self.RawArray[i][32] + self.RawArray[i][33]*256 for i in range(0,self.Count)]#儒略日
        self.NumberofElevation = [self.RawArray[i][44] + self.RawArray[i][45]*256 for i in range(0,self.Count)]#层数
        self.StartOfReflectivity = [self.RawArray[i][46]+self.RawArray[i][47]*256 for i in range(0,self.Count)]#起始距离
        self.StartOfSpeed = [self.RawArray[i][48] + self.RawArray[i][49]*256 for i in range(0,self.Count)]
        self.StepOfReflectivity = [self.RawArray[i][50] + self.RawArray[i][51]*256 for i in range(0,self.Count)]#库长
        self.StepOfSpeed = [self.RawArray[i][52]+self.RawArray[i][53]*256 for i in range(0,self.Count)]
        self.NumberOfReflectivity = [self.RawArray[i][54] + self.RawArray[i][55]*256 for i in range(0,self.Count)]#反射率的距离库数
        self.NumberOfSpeed = [self.RawArray[i][56] + self.RawArray[i][57]*256 for i in range(0,self.Count)]
        self.PointerOfReflectivity = [self.RawArray[i][64]+self.RawArray[i][65]*256 for i in range(0,self.Count)]#数据位置指针
        self.PointerOfSpeed = [self.RawArray[i][66]+self.RawArray[i][67]*256 for i in range(0,self.Count)]
        self.PointerOfSpectralWidth = [self.RawArray[i][68]+self.RawArray[i][69]*256 for i in range(0,self.Count)] # 谱宽指针
        self.ResolutionOfSpeed = [self.RawArray[i][70]+self.RawArray[i][71]for i in range(0,self.Count)]#多普勒速度分辨率
        self.vcp = [self.RawArray[i][72]+self.RawArray[i][73]*256 for i in range(0,self.Count)]# 11:降水，16；21：降水，14；31：晴空，8；32：晴空，7.
        self.Elevation = [(self.RawArray[i][42]+256*self.RawArray[i][43])/8*180/4096 for i in range(0,self.Count)]#仰角
        self.Azimuth = [(self.RawArray[i][36]+256*self.RawArray[i][37])/8*180/4096 for i in range(0,self.Count)]#方位角
        self.Storage = self.getStorage()
        #self.AllInfo = self.getAllInfo()
        #self.x,self.y,self.z,self.r = self.getXyzr()
        #self.space_info = self.getAllInfo()
        #self.elevation_list = self.get_elevation_list()
        #self.grid_data = self.grid()
        file.close()

需要注意的是从雷达基数据解析得到的是反射率数据dbz，对于深度学习模型，常见的图像数据类型是值域在0～255的可见光图像，所以需要把dbz映射到0～255，可以参考下面的公式进行实现

$pixel = 255 * \frac{dbz+10}{70} + 0.5$

依照这个公式计算得到的pixel是浮点类型的数据，还需要进一步转换成无符号整形数据（uint8)。

算法模型

雷达图像外推与时空相关，即受到时间影响也跟空间运动有关系。所以模型应该同时能够捕捉连续时间图像之间的时序特征以及图像中的空间运动特征。卷积网络在提取图像中静态特征，描述图像方面性能优异，LSTM网络则是对时序特征有较好的感受度，所以对这两种模型进行融合用于提取空间特征和时序特征是一个较为理想的选择。基于这个思路，Xingjian Shi 团队提出了convLSTM模型用于解决雷达外推问题。

与FC-LSTM的不同点是状态之间的转换以及门的状态从全连接计算变成了卷积运算。

下面是用作实验的convLSTM网络架构，使用了5层ConvLSTM2D，在层与层之间进行BatchNormalization，最后使用一个Conv3D卷积层输出推理结果。

Model: "model"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
content_input (InputLayer)   [(None, None, 64, 64, 1)] 0         
_________________________________________________________________
layer1 (ConvLSTM2D)          (None, None, 64, 64, 64)  150016    
_________________________________________________________________
batch_normalization (BatchNo (None, None, 64, 64, 64)  256       
_________________________________________________________________
layer2 (ConvLSTM2D)          (None, None, 64, 64, 64)  295168    
_________________________________________________________________
batch_normalization_1 (Batch (None, None, 64, 64, 64)  256       
_________________________________________________________________
layer3 (ConvLSTM2D)          (None, None, 64, 64, 64)  295168    
_________________________________________________________________
batch_normalization_2 (Batch (None, None, 64, 64, 64)  256       
_________________________________________________________________
layer4 (ConvLSTM2D)          (None, None, 64, 64, 64)  295168    
_________________________________________________________________
batch_normalization_3 (Batch (None, None, 64, 64, 64)  256       
_________________________________________________________________
layer5 (ConvLSTM2D)          (None, None, 64, 64, 64)  295168    
_________________________________________________________________
batch_normalization_4 (Batch (None, None, 64, 64, 64)  256       
_________________________________________________________________
conv3d (Conv3D)              (None, None, 64, 64, 1)   1729      
=================================================================
Total params: 1,333,697
Trainable params: 1,333,057
Non-trainable params: 640
_________________________________________________________________

模型采用的是类似encoder-forecater架构。

基于keras实现的网络架构代码

def Net():
    contentInput = Input(shape=(None,Config.WIDTH,Config.HEIGHT,1),name='content_input')

    x = ConvLSTM2D(64,(3,3),padding='same',dropout=0.2,return_sequences=True,kernel_initializer='he_normal',name='layer1')(contentInput)
    x = BatchNormalization()(x)
    x = ConvLSTM2D(64,(3,3),padding='same',dropout=0.2,return_sequences=True,kernel_initializer='he_normal',name='layer2')(x)
    x = BatchNormalization()(x)
    x = ConvLSTM2D(64,(3,3),padding='same',dropout=0.2,return_sequences=True,kernel_initializer='he_normal',name='layer3')(x)
    x = BatchNormalization()(x)
    x = ConvLSTM2D(64,(3,3),padding='same',dropout=0.2,return_sequences=True,kernel_initializer='he_normal',name='layer4')(x)
    x = BatchNormalization()(x)
    x = ConvLSTM2D(64,(3,3),padding='same',dropout=0.2,return_sequences=True,kernel_initializer='he_normal',name='layer5')(x)
    x = BatchNormalization()(x)

    predictions = Conv3D(filters=1,kernel_size=(3,3,3),activation='linear',padding='same',data_format='channels_last')(x)
    model = Model(inputs=contentInput,outputs=predictions)
    return model

为了避免L1、L2正则化引入的模糊，使用logcosh损失函数。鉴于数据特殊性，使用多种评价准则，实验发现acc和loss无法准确的表征模型的性能。

在模型训练时使用acc、loss、psnr、ssim、POD、FAR等要素综合更新网络参数。

acc：准确率

loss：logcosh值

psnr：峰值信噪比

ssim：图像相似度

POD：命中率

FAR：误警率

模型表现

模型外推结果

输入8张当前时序雷达反射率图像，向外推16次得到16*10分钟外推结果，如果使用的训练数据是6分钟/次，则是16*6分钟的外推结果。

可以看到，外推结果中存在大量噪声，并且没有捕捉到原图像中的非平移运动。这是因为对于卷积网络，每一个位置的局部连接都是固定的，也即是local-invariant，局部不变性，但是对于自然界中的运动，很多都是local-variant，如旋转和消散都是动态的，一团云在某一帧的位置对应到下一帧的位置不一样，所以需要建立一个动态的链接。这是当前ConvLSTM无法做到的，在下一篇文章中会讲解local-variant网络TrajGRU。

整理数据、调研资料、敲代码不易，所以下载资源会收取少量费用作为更新和研究的动力，对于付费下载的朋友，博主会一直提供长期技术支持，希望大家一起进步，在各自的领域各有所成。

完整代码和模型下载连接

有需要别的材料支持的，可以加QQ群：897283503