本文详细介绍了基于TrajGRU的短临预报模型的代码实现，涵盖雷达数据预处理、PyTorch模型构建、训练参数设置、模型测试验证和部署推理模型的步骤。重点讨论了数据解码、噪声剔除和TrajGRU在网络中的作用，特别是其在捕捉光流运动中的关键作用。
前言
接上篇一、基于TrajGRU的短临预报实现：模型结构，本篇文章主要是关于代码实现以及训练时的一些调优技巧和经验型的错误。
本文分为以下几个部分：
一、雷达数据预处理
二、模型代码pytorch实现
三、训练参数设置
四、模型测试验证
五、部署推理模型
本文依次根据上述章节进行展开。
chapter1 雷达数据预处理
雷达基数据以二进制存储，需要按位读取并解码成十进制数据。目前使用的数据是S波段多普勒雷达数据，雷达扫描方式是VCP11或者VCP21，二者的区别是仰角不同。美国雷达的标准是VCP-11有14个仰角，仰角范围是0.5～19.5度，VCP-21有9个仰角，仰角范围是0.5～19.5度。我国S波段多普勒雷达使用CINRAD 标准，
11：降水模式，16层仰角
21：降水模式，14层仰角
31：晴空模式，8层仰角
32：晴空模式，7层仰角
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()多普勒雷达每个仰角的扫描方式是PPI，每个仰角扫描得到的数据长度都是固定的，存储长度是2432个字节，数据格式如下：
字节顺序 双字节顺序 数据类型 说明
1-14 1-7 保留 雷达信息头
（28字节）
15-16 8 2字节 1-表示雷达数据
17-28 9-14 保留
29-32 15-16 4字节 径向数据收集时间(毫秒,自00:00开始)
33-34 17 2字节 儒略日（Julian）表示，自1970年1月1日开始
35-36 18 2字节 不模糊距离（表示：数值/10.=千米）
37-38 19 2字节 方位角（编码方式：[数值/8.]*[180./4096.]=度）
39-40 20 2字节 当前仰角内径向数据序号
41-42 21 2字节 径向数据状态 0：该仰角的第一条径向数据
1：该仰角中间的径向数据
2：该仰角的最后一条径向数据
3：体扫开始的第一条径向数据
4：体扫结束的最后一条径向数据
43-44 22 2字节 仰角 （编码方式：[数值/8.]*[180./4096.]=度）
45-46 23 2字节 体扫内的仰角数
47-48 24 2字节 反射率数据的第一个距离库的实际距离(单位:米)
49-50 25 2字节 多普勒数据的第一个距离库的实际距离(单位:米)
51-52 26 2字节 反射率数据的距离库长（单位：米）
53-54 27 2字节 多普勒数据的距离库长（单位：米）
55-56 28 2字节 反射率的距离库数
57-58 29 2字节 多普勒的距离库数
59-60 30 2字节 扇区号
61-64 31-32 4字节 系统订正常数
65-66 33 2字节 反射率数据指针（偏离雷达数据信息头的字节数）
表示第一个反射率数据的位置
67-68 34 2字节 速度数据指针（偏离雷达数据信息头的字节数）
表示第一个速度数据的位置
69-70 35 2字节 谱宽数据指针（偏离雷达数据信息头的字节数）
表示第一个谱宽数据的位置
71-72 36 2字节 多普勒速度分辨率。 2：表示0.5米/秒
4：表示1.0米/秒
73-74 37 2字节 体扫（VCP）模式 11：降水模式，16层仰角
21：降水模式，14层仰角
31：晴空模式，8层仰角
32：晴空模式，7层仰角
75-82 38-41 保留
83-84 42 2字节 用于回放的反射率数据指针，同33
85-86 43 2字节 用于回放的速度数据指针，同34
87-88 44 2字节 用于回放的谱宽数据指针，同35
89-90 45 2字节 Nyquist速度（表示：数值/100. = 米/秒）
91-128 46-64 保留
129-588 65-294 1字节 反射率
距离库数：0-460
编码方式：（数值-2）/2.-32 = DBZ
当数值为0时，表示无回波数据（低于信噪比阀值）
当数值为1时，表示距离模糊 基数据
部分
(2300字节)
129-1508 65-754 1字节 速度
距离库数：0-920
编码方式：
分辨率为0.5米/秒时
（数值-2）/2.-63.5 = 米/秒
分辨率为1.0米/秒时
（数值-2）-127 = 米/秒
当数值为0或1时，意义同上
129-2428 65-1214 1字节 谱宽
距离库数：0-920
编码方式：
（数值-2）/2.-63.5 = 米/秒
当数值为0或1时，意义同上
2429-2432 1215-1216 保留
说明：
1． 数据的存储方式
每个体扫存储为一个单独的文件
2． 数据的排列方式
按照径向数据的方式顺序排列，对于CINRAD SA/SB雷达，体扫数据排列自低仰角开始到高仰角结束。
3． 径向数据的长度
径向数据的长度固定，为2432字节。
4． 距离库长和库数
反射率距离库长为1000米，最大距离库数为460；
速度和谱宽距离库长为250米，最大距离库数为920。经过上述方法解析，可以得到雷达在各个仰角的基本反射率数据。解析后的雷达反射率数据经过变换可以得到0-255范围内的可见光图像。
常用的转化公式
得到连续的雷达反射率拼图后，还需要对数据做进一步的处理。在晴天云量较少时，雷达反射率较小甚至没有，这一部分数据对于模型来说等同于噪声，所以需要剔除。
另外，在剔除数据值一定要注意保持数据时空连续性，比如十分钟一次的雷达反射率组合拼图，剔除噪声数据之后需要保证剩下的图像仍然是以十分钟为间隔的。
chapter2 模型代码pytorch实现
模型采用的是Encoder-Decoder架构，论文中叫Encoder-Forecaster，是相同的意思。前文已经提到过这个架构，这里主要是代码实现。
在论文中，Encoder阶段是通过卷积+trajGRU实现图像降采样和光流以及时序特征提取，Forecaster阶段相当于Encoder的镜像，即通过trajGRU+卷机实现光流以及时序特征提取和图像上采样，最终输出预测结果。
Encoder
Encoder由三层Conv+trajGRU堆叠而成。
conv -> rnn -> conv -> rnn -> conv -> rnn ->states在卷积阶段通过设置卷积核尺寸步长和padding大小进行特征提取和图像降采样。
三层卷机的设置分别如下：
stage1
1,8,7,5,1 -> inputchannels,outputchannels,kernel-size,stride,padding
(N + 2*1 - 7)/5 + 1 = N'
N' = N/5
stage2
64,192,5,3,1
(N + 2*1 - 5)/3 + 1 = N'
N' = N/3
stage3
192,192,3,2,1
(N + 2*1 - 3)/2 + 1 = N'
N' = N/2可以看到，经过第一次卷积之后，输出尺寸是输入的1/5，经过第二次卷积之后，输出尺寸是输入的1/3，经过第三次卷积之后，输出尺寸是输入的1/2。
基于torch的实现如下：
stage1
conv2d = nn.Conv2d(1,8,7,5,1)
relu = nn.LeakyReLU(negative_slope=0.2, inplace=True)
stage2
conv2d = nn.Conv2d(64,192,5,3,1)
relu = nn.LeakyReLU(negative_slope=0.2, inplace=True)
stage3
conv2d = nn.Conv2d(192,192,3,2,1)
relu = nn.LeakyReLU(negative_slope=0.2, inplace=True)
rnn结构（TrajGRU）
这一部分可以说是整个模型架构的核心部分，也是能够成功实现外推的关键。
上篇已经分析过TrajGRU的网络结构，本文在这里结合代码实现进一步说明。
TrajGRU和GRU的区别主要有两点：一、前者使用卷积结构实现状态转移，后者使用全连接结构；二、前者比后者多了一个光流运动捕捉结构，即traj结构部分。
先来说一下traj结构，traj结构通过引入一层卷积结构，同时限制输出的维度，再经过不断的迭代调整卷积参数达到捕捉复杂光流运动的效果。需要注意的是，作者通过限制输出的维度是2维其实是模拟了经典光流算法的输出设置，所以traj结构的输出其实是前一帧加上光流矢量之后的效果，也即是某点在下一时刻图像中的位置。由于位置通常是离散不可微分，所以利用了一个插值方法取所在位置索引对应的特征值作为反馈学习的对象。
结构组成
1、input2hidden
2、inputs2flow
3、hidden2flow
4、flows_conv -->> wrap function -->> pixel value
5、ret（reset gate）trajGRU网络结构由上述5个部分组成，实现从输入到光流运动学习以及状态转移等。最终利用3个trajGRU的堆叠实现雷达回波外推。
需要注意的是，trajGRU结构不会改变输入输出的尺寸，即特征的宽和高不会被改变，只在卷积层改变特征的尺寸。一种合理的解释是类似LSTM和GRU这种含有记忆单元的结构，只是产生了状态转移，没有发生尺寸变换。
trajGRU代码实现
class TrajGRU(BaseConvRNN):
# b_h_w: input feature map size
def __init__(self, input_channel, num_filter,b_h_w,zoneout=0.2,L=5,
i2h_kernel=(3,3), i2h_stride=(1,1), i2h_pad=(1,1),
h2h_kernel=(5,5), h2h_dilate=(1,1), prefix='BaseConvRNN'):
super(TrajGRU,self).__init__(num_filter=num_filter,
b_h_w=b_h_w,
h2h_kernel=h2h_kernel,
h2h_dilate=h2h_dilate,
i2h_kernel=i2h_kernel,
i2h_stride=i2h_stride,
i2h_pad=i2h_pad,
prefix=prefix)
self._L = L
self._zoneout = zoneout
#self._act_type = F.leaky_relu()
# *3 according to 3 hidden states. t-1,t,t+1
## self.i2h -->> input2hidden
self.i2h = nn.Conv2d(in_channels=input_channel,
out_channels=self._num_filter*3,
kernel_size=self._i2h_kernel,
stride=self._i2h_stride,
padding=self._i2h_pad,
dilation=self._i2h_dilate)
# inputs to flow
self.i2f_conv1 = nn.Conv2d(in_channels=input_channel,
out_channels=32,
kernel_size=(5,5),
stride=1,
padding=(2,2),
dilation=(1,1))
# hidden to flow
self.h2f_conv1 = nn.Conv2d(in_channels=self._num_filter,
out_channels=32,
kernel_size=(5,5),
stride=1,
padding=(2,2),
dilation=(1,1))
# generate flow
self.flows_conv = nn.Conv2d(in_channels=32,
out_channels=self._L*2,
kernel_size=(5,5),
stride=1,
padding=(2,2))
# reset gate hh,hz,hr,1*1 ks
self.ret = nn.Conv2d(in_channels=self._num_filter*self._L,
out_channels=self._num_filter*3,
kernel_size=(1,1),
stride=1)
# inputs: B C H W
# flow comes from current inputs and forward states
def _flow_generator(self,inputs,states):
if inputs is not None:
i2f_conv1 = self.i2f_conv1(inputs)
else:
i2f_conv1 = None
h2f_conv1 = self.h2f_conv1(states)
f_conv1 = i2f_conv1 + h2f_conv1 if i2f_conv1 is not None else h2f_conv1
f_conv1 = F.leaky_relu(f_conv1,0.2,inplace=True)
flows = self.flows_conv(f_conv1)
# channels L*2 ,split according to 2
# get L flow maps each have 2 channels
flows = torch.split(flows,2,dim=1)
return flows
# inputs states
# inputs: S B C H W
def forward(self,inputs=None,states=None,seq_len=5):
if states is None:
states = torch.zeros((inputs.size(1),self._num_filter,self._state_height,
self._state_width),dtype= torch.float).cuda()
if inputs is not None:
S,B,C,H,W = inputs.size()
i2h = self.i2h(torch.reshape(inputs,(-1,C,H,W)))
i2h = torch.reshape(i2h,(S,B,i2h.size(1),i2h.size(2),i2h.size(3)))
i2h_slice = torch.split(i2h,self._num_filter,dim=2)
else:
i2h_slice = None
prev_h = states
outputs = []
for i in range(seq_len):
if inputs is not None:
flows = self._flow_generator(inputs[i,...],prev_h)
else:
flows = self._flow_generator(None,prev_h)
wrapped_data = []
for j in range(len(flows)):
flow = flows[j]
wrapped_data.append(wrap(prev_h,-flow))
wrapped_data = torch.cat(wrapped_data,dim=1)
h2h = self.ret(wrapped_data)
h2h_slice = torch.split(h2h,self._num_filter,dim=1)
if i2h_slice is not None:
reset_gate = torch.sigmoid(i2h_slice[0][i,...]+h2h_slice[0])
update_gate = torch.sigmoid(i2h_slice[1][i,...]+h2h_slice[1])
new_mem = F.leaky_relu(i2h_slice[2][i, ...] + reset_gate * h2h_slice[2],0.2,inplace=True)
else:
reset_gate = torch.sigmoid(h2h_slice[0])
update_gate = torch.sigmoid(h2h_slice[1])
new_mem = F.leaky_relu(reset_gate * h2h_slice[2],0.2,inplace=True)
next_h = update_gate * prev_h + (1-update_gate) * new_mem
if self._zoneout > 0.0:
mask = F.dropout2d(torch.zeros_like(prev_h),p=self._zoneout)
next_h = torch.where(mask,next_h,prev_h)
outputs.append(next_h)
prev_h = next_h
return torch.stack(outputs),next_h整体的实现还是严格按照网络结构进行的，继承了一个BaseConvRNN类，用来计算网络特征转移中的特征尺寸以及padding设置等。
下一篇按照trajGRU的5个结构拆解开逐个分析。
扫一扫
1516
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
