利用机器学习进行降水临近预报

ML、DL、CNN、天气预报、UNet
利用机器学习进行降水临近预报 cover image

如今，气象学家估计 5 天内的天气预报有 90% 是正确的。做出的预测通常基于两种不同的方法：

  1. 基于物理的方法：这些方法使用集成可测量量的模型，例如压力、云的运动、天空条件……此类模型擅长预测未来几天或几周的天气。

  2. 不受物理影响（基于数据）的方法：这些方法使用历史数据来创建可以进行预测的模型。此类模型在预测长达 6 小时的天气或所谓的临近天气预报方面显示出良好的效果。

在本文中，我们将讨论第二类方法。我们将讨论天气数据的不同格式、如何利用机器学习来进行天气临近预报，以及开发的架构如何有利于解决类似问题。

气候预测数据

由于无物理方法使用历史数据，因此我们首先研究可用的数据。

我们将使用两个主要数据源：

  1. 图像数据：该数据采用指定地理区域的雷达或卫星图像的形式。它用于预测降水、风向或湿度。

Satellite Image

  1. 表格数据：该数据采用可测量量记录的形式，例如温度、湿度或风速。

Tabular Data

虽然这两个数据源对于构建强大的模型都很重要，但出于简单原因，我们将重点关注前者（从雷达或卫星收集的图像数据）。最常用的图像数据模型是卷积神经网络 (CNN ）。

在此工作之后，我们将使用 U-Net 架构来构建我们自己的即时预报模型。

＃＃ 建筑学

从现有架构开始有很多帮助，其中包括：

架构作为创建新模型的指导方针。创建新架构的人们采用反复试验的方法来获得最终结果。通过重用他们的最终结果，我们可以节省大量时间。

预训练模型通常可以立即使用。当研究人员发布他们的新架构时，他们通常也会发布经过训练的参数，这样用户就不必经历从头开始训练/优化的麻烦。这对于非常大的、资源匮乏的模型特别有用。

著名的视觉架构示例包括：

  • LeNet（60k 参数）

LeNet

  • AlexNet（60m参数）

AlexNet

  • VGG-16（138m参数）

VGG-16

U网

U-Net是一种基于全卷积网络的架构，这意味着它没有任何全连接层。它首次被引入用于医学图像分割任务。这些结果启发研究人员将其扩展到计算机视觉的其他任务。

2019年，谷歌使用基于U-Net架构创建了降水预报模型。

“U-Net”这个名字来源于其架构的“U”形。

U-net Architecture

我们确定了三个主要组成部分：

  1. 收缩/编码器：连续的卷积/池化层，将输入图像压缩为较小尺寸的表示。

  2. 桥/瓶颈：“U”的底部连接编码器和解码器。它是由一系列卷积运算形成的。

  3. 解压缩/解码器：连续的上卷积层和卷积层，“解压缩”瓶颈的输出。

到目前为止，U-Net 的架构看起来就像一个自动编码器。然而，差异在于编码器和解码器之间传递的信息。该信息是通过将编码器的卷积结果与解码器的上卷积结果连接起来传递的。此修改提高了分辨率，并允许模型输出空间上更精确的输出（输出中像素的位置将更精确，解决了不包含串联的模型中出现的问题）。连接是对称完成的。

构建临近降水预报的 U-Net 模型

在本节中，我们将构建一个 U-Net 模型来根据卫星图像预测天气。我们将使用一个名为 Rain-Net 的预训练模型。

以下代码可在此 colab 中找到。

我们首先安装wradlib，一个用于天气雷达数据处理的开源库

!pip install wradlib
import wradlib as wrl

然后，我们编写两个实用函数来从 DWD 的开放数据服务器 下载卫星数据


import urllib.request
import io
import numpy as np
import datetime

def download_and_read_RY(RY_timestamp):
   url = f"https://opendata.dwd.de/weather/radar/radolan/ry/raa01-ry_10000-{RY_timestamp}-dwd---bin"
   data_binary = urllib.request.urlopen(url).read()
   data, attr = wrl.io.read_radolan_composite(io.BytesIO(data_binary), missing = 0)
   data = data.astype("float32")
   return data,attr

def download_data():
   latest_scan, latest_attr = download_and_read_RY("latest")
   latest_datetime = latest_attr["datetime"]
   list_for_downloading = [ts.strftime("%y%m%d%H%M") for ts in
       [latest_datetime - datetime.timedelta(minutes=t) for t in [15, 10, 5]]]
   previous_scans = np.array([download_and_read_RY(ts)[0] for ts in list_for_downloading])
   print(list_for_downloading)
   scans = np.concatenate([previous_scans, latest_scan[np.newaxis, ::, ::]], axis = 0)

   return scans, latest_datetime

这些实用程序允许我们下载最新的 4 个卫星图像，这是我们使用预训练模型进行预测所需的卫星图像数量。

然后我们可以使用创建的函数来获取最新的 4 张图像


RY_latest, RY_latest_timestep = download_data()

获取图像后，我们使用wradlib的vis.plot_ppi方法来绘制数据

for i in range(RY_latest.shape[0]):
   wrl.vis.plot_ppi(RY_latest[i])

VIS Radar

现在我们已经加载了数据。接下来让我们加载模型。

我们首先导入相关的类。我们将在本文中使用 TensorFlow。

从tensorflow.keras.layers导入输入、Conv2D、激活、连接、Conv2DTranspose、MaxPool2D

从tensorflow.keras.models导入模型

让我们构建 3 个原始构建块。根据此实现，这些“构建块”将用于创建整个架构。

第一个块对应于连续的卷积层，我们称之为“conv_block”

def conv_block(input, num_filters):
   x = Conv2D(num_filters, 3, padding ="same", kernel_initializer="he_normal"(input))
   x = Activation("relu")(x)
   x = Conv2D(num_filters, 3, padding ="same", kernel_initializer="he_normal"(x))
   x = Activation("relu")(x)
   return x

第二个块用于构建编码器部分（卷积块+最大池化）。我们称之为“encoder_block”

def encoder_block(input, num_filters):
   x = conv_block(input, num_filters)
   p = MaxPool2D((2,2))(x)
   return x,p

第三个也是最后一个块是“decoder_block”（上卷积+串联+卷积）。

def decoder_block(input, skip_features, num_filters):
   x = Conv2DTranspose(num_filters, (2,2), strides=2, padding="same", kernel_initializer="he_normal")
   x = Concatenate()([x, skip_features])
   x = conv_block(x, num_filters)
   return x

我们结合这些构建块来构建 U-Net 模型

def build_unet(input_shape = (928, 928, 4)):

   inputs = Input(input_shape)
   e1, p1 = encoder_bock(inputs, 64)
   e2, p2 = encoder_bock(p1, 128)
   e3, p3 = encoder_bock(p2, 256)
   e4, p4 = encoder_bock(p3, 512)

   b1 = conv_block(p4, 1024)

   d1 = decoder_block(b1, e4, 512)
   d2 = decoder_block(d1, e3, 256)
   d3 = decoder_block(d2, e2, 128)
   d4 = decoder_block(d3, e1, 64)

   outputs = Conv2D(1, 1, padding="same", activation="linear")(d4)

   unet = Model(inputs, outputs, name="U-Net")

   return unet

如果您想将其用于其他用途，请随意尝试实施、改进它或根据您的需求进行调整。

训练这个模型可能需要很多时间。幸运的是，有一个名为 Rain-Net 的模型，它是基于 U-Net 架构创建的，专门用于降水临近预报。

让我们克隆它的 GitHub 存储库

! git clone https://github.com/hydrogo/rainnet.git

然后我们下载该模型的预训练权重

!wget -O /content/rainnet/rainnet_weights.h5 https://zenodo.org/record/3630429/files/rainnet_weights.h5

下一步是根据存储库中找到的架构创建一个模型，然后加载下载到该模型的权重

import sys
from rainnet import rainnet
model = rainnet.rainnet()
model.load_weights('/content/rainnet/rainnet_weights.h5')
model.summary()

我们下载的图像大小为 900*900 像素。我们将重塑这些图像以匹配 Rain-Net 的预期输入

def Scaler(array):
   return np.log(array+0.01)


def invScaler(array):
   return np.exp(array) - 0.01


def pad_to_shape(array, from_shape=900, to_shape=928, how="mirror"):
   # calculate how much to pad in respect with native resolution
   padding = int( (to_shape - from_shape) / 2)
   # for input shape as (batch, W, H, channels)
   if how == "zero":
       array_padded = np.pad(array, ((0,0),(padding,padding),(padding,padding),(0,0)), mode="constant", constant_values=0)
   elif how == "mirror":
       array_padded = np.pad(array, ((0,0),(padding,padding),(padding,padding),(0,0)), mode="reflect")
   return array_padded


def pred_to_rad(pred, from_shape=928, to_shape=900):
   # pred shape 12,928,928
   padding = int( (from_shape - to_shape) / 2)
   return pred[::, padding:padding+to_shape, padding:padding+to_shape].copy()

'''
the spatial extent of input data has to be a multiple of 2n+1
where n is the number of max pooling layers
'''

def data_preprocessing(X):

   # 0. Right shape for batch
   X = np.moveaxis(X, 0, -1)
   X = X[np.newaxis, ::, ::, ::]
   # 1. To log scale
   '''
   We use a log scale to respond to skewness towards large values
   '''
   X = Scaler(X)
   # 2. from 900x900 to 928x928
   X = pad_to_shape(X)

   return X


def data_postprocessing(nwcst):

   # 0. Squeeze empty dimensions
   nwcst = np.squeeze(np.array(nwcst))
   # 1. Convert back to rainfall depth
   nwcst = invScaler(nwcst)
   # 2. Convert from 928x928 back to 900x900
   nwcst = pred_to_rad(nwcst)
   # 3. Return only positive values
   nwcst = np.where(nwcst>0, nwcst, 0)
   return nwcst

然后我们创建一个进行预测的函数。

def prediction(model_instance, input_data, lead_time=24):

   input_data = data_preprocessing(input_data)
   nwcst = []

   for _ in range(lead_time):
       pred = model_instance.predict(input_data)
       nwcst.append(pred)
       input_data = np.concatenate([input_data[::, ::, ::, 1:], pred], axis = -1)

   nwcst = data_postprocessing(nwcst)

   return nwcst

然后我们对之前下载的数据调用这个函数

Y_pred = prediction(model, RY_latest)

我们可以绘制预测，并保存它们以使用保存的结果创建 gif 图像，使我们能够可视化预测

import matplotlib.pyplot as plt
names = []
for i in range(0, 19):
   im = wrl.vis.plot_ppi(Y_pred[i])
   name = '/content/image'+str(i)+'.png'
   names.append(name)
   plt.savefig(name)

import imageio
from google.colab import files

imgs = []

for name in names:
   imgs.append(imageio.imread(name))

imageio.mimsave('/content/pred.gif', imgs, fps=4)
files.download('/content/pred.gif')

Resuts

恭喜您已经走到这一步了！您现在可以使用 Rain-Net 进行预测并将其可视化。

＃＃ 结论

在本文中，我们使用机器学习模型（Rain-Net）进行降水临近预报。我们使用了使用 TensorFlow 构建的 U-Net 架构。我们加载了一个预训练模型来预测卫星图像。

该实施可以通过多种方式进行改进。例如：

  1. 在数据集上微调模型

  2. 在架构中使用基于注意力的模块，例如 CBAM（Convolutional Block Attention Module）。

＃＃ 参考

来参加我们的免费研讨会之一！

通过我们的免费研讨会 开始您作为数据科学家的职业生涯，这些研讨会基于适应性强的课程并由行业专家指导。

Career Services background pattern

职业服务

Contact Section background image

让我们保持联系

概述

关于

课程

法律

联系我们

跟着我们:

Code Labs Academy © 2024 版权所有.