変分オートエンコーダを用いた気象データ画像の生成
1. はじめに
画像生成の分野などで変分オートエンコーダ(Variational Autoencoder)が用いられる場合があります。
Variational Autoencoder(以下、VAE)については既に多くの解説がなされています。
理論的なところは下記を参考にさせて頂きました。
この手法をいくつかの気象データに適用してみました。
2. 手法
2.1 ネットワーク
下記のサイトのソースをほぼ流用させて頂きました。
【Python】Keras で VAE 入門
ネットワークは入力画像から潜在空間へのマッピングを計算する「エンコーダ」と、潜在空間上の指定した値から出力画像を生成する「デコーダ」からなります。
2.2 エンコーダ
対象とする気象画像は２次元カラー画像（縦、横、チャネル数）です。
エンコーダとしては、strides=(2,2)のコンボリューション層（Conv2D）によってサイズを圧縮したのち、Flatten層により一次元化して、Dense層によって2種類のパラメタに縮約します。
この2種類のパラメタはVAEで仮定する正規分布関数の「平均」と「分散」に相当します。
（下記ソースのz_meanとz_log_var）
パラメタを用いて、正規分布関数から値zを取得するのがエンコーダ部分です。
def sampling(args):
z_mean, z_log_var = args
epsilon = K.random_normal(shape=(K.shape(z_mean)[0], latent_dim) ,
mean=0., stddev=1.)
return z_mean + K.exp(z_log_var) * epsilon
cnn_input = Input(shape=(i_dmlat, i_dmlon, in_ch))
cnv2d1 = Conv2D(num_hidden1, kernel_size=(3,3), activation='relu',
padding='same')(cnn_input)
cnv2d1 = Conv2D(num_hidden2, kernel_size=(3,3), strides=(2,2),
activation='relu', padding='same')(cnv2d1)
cnv2d1 = Conv2D(num_hidden2, kernel_size=(3,3), activation='relu',
padding='same')(cnv2d1)
cnv2d1 = Conv2D(num_hidden3, kernel_size=(3,3), strides=(2,2),
activation='relu', padding='same')(cnv2d1)
cnv2d1 = Conv2D(num_hidden3, kernel_size=(3,3),
activation='relu', padding='same')(cnv2d1)
shape_before_flattening = K.int_shape(cnv2d1)
cnv2d1 = Flatten()(cnv2d1)
cnv2d1 = Dense(num_hidden2, activation='relu')(cnv2d1)
z_mean = Dense(latent_dim)(cnv2d1)
z_log_var = Dense(latent_dim)(cnv2d1)
z = Lambda(sampling)([z_mean , z_log_var ])
encoder = Model(input=cnn_input , output=z )
####2.3 デコーダ
エンコーダが作ったzから、エンコーダが一次元化する直前の２次元形式のデータに変換します。これはDense層とReshape層が行っています。
その上で、Conv2Dtransposeを用いて画像を元のサイズに戻します。
この部分がデコーダネットワークです。
decoder_input = Input(K.int_shape(z)[1:])
x = Dense(np.prod(shape_before_flattening[1:]), activation='relu')(decoder_input)
x = Reshape(shape_before_flattening[1:])(x)
x = Conv2DTranspose(num_hidden2, kernel_size=(3,3), strides=(2,2),
activation='relu', padding='same')(x)
x = Conv2DTranspose(num_hidden1, kernel_size=(3,3), strides=(2,2),
activation='relu', padding='same')(x)
x = Conv2D(3, kernel_size=(3,3), activation='sigmoid' , padding='same')(x)
decoder = Model(input=decoder_input , output=x )
z_decoded = decoder(z)
####2.4 ネットワーク全体、変分下限を求めるロス関数
エンコーダ部分にcnn_inputを渡すとデコーダがz_decodedを出力するというのが全体のネットワークです。
このネットワークについて、尤度最大化を行うために「変分下限」を、カスタム損失関数レイヤCunstomVariationalLayerに入力とデコーダ出力を渡して計算させています。
（call はtensorflow内部処理のために必要なのだそうで、直接使用しません）
学習させるネットワークvaeとしては、入力がcnn_inputで出力を損失レイヤの出力yとして作成します。
y = CustomVariationalLayer()([cnn_input, z_decoded])
vae = Model(cnn_input, y )
class CustomVariationalLayer(Layer):
def vae_loss(self, x , z_decoded):
x = K.flatten(x)
z_decoded = K.flatten(z_decoded)
xent_loss = keras.metrics.binary_crossentropy(x,z_decoded)
kl_loss = -5e-4 * K.mean( 1 + z_log_var - K.square(z_mean)
- K.exp(z_log_var) , axis=-1)
return K.mean(xent_loss + kl_loss)
def call(self, inputs):
x = inputs[0]
z_decoded = inputs[1]
loss = self.vae_loss(x, z_decoded)
self.add_loss(loss, inputs=inputs)
return x
損失関数のネットワークへの組み込みについて
この実装はfchollet/deep-learning-with-python-notebooksが元になっていると思いますが、ネットワークvaeの出力は損失関数の出力となっています。
また、compileとfitは以下のようになっていて、trainingの際に教師信号として何も与えていません。
普通に考えると、compileの時にカスタム損失関数を指定して、fitの時には入力はcnn_input, 出力はz_docededに相当する値を渡しそうなものです。どういう理由でこういう実装としているのか、調べきれませんでしたが、
Machine Learning Explained, Machine Learning Tutorialsには通常の構成方法に見える実装ソースが書かれております。
vae.compile(optimizer='rmsprop' , loss=None) # <- lossにNoneが！
vae.fit(x=np_i_data_train, y=None, epochs=num_itter , callbacks=cbks, batch_size=28 , validation_split=0.1 ) # <- yにNoneが！
3. 気象データの生成
###3.1 入力データ
以下の3種類のデータを対象としてみました。
（１）気象庁全国合成レーダー
（２）気象庁メソ予報モデル(MSM)地上気温データ
（３）気象庁速報天気図（SPAS)
（１）と（２）は京都大学生存圏研究所 生存圏データベースが提供しているGPVデータをダウンロードして可視化した画像を用います。
（３）については気象庁が公開しているカラー版の天気図を使用します。
GPVの可視化にはpythonモジュールであるpygribが便利ですが、（１）の全国合成レーダーのGPVに関してはpygribではオープンできないので一度NetCDF形式に変換するということをしております。（２）についてはpygribを用いてデータを操作しました。このあたりのお話については、拙稿ながらQiitaの以下をご参照いただけると幸いです。
気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（2)
pygribでオープンできないGRIB2形式の気象データをnetCDFに変換して可視化する
####3.1.1 全国合成レーダーの可視化画像例
気象庁全国合成レーダーは下記のような画像に可視化しました（2020年3月31日21時UTCの可視化画像）。
####3.1.2 気象庁メソ予報モデル地上気温データの可視化画像例
気象庁メソ予報モデル地上気温データは下記のような画像に可視化しました（2020年3月27日18時UTCの可視化画像）。
###3.2 学習結果
####3.2.1 全国合成レーダーの結果
#####3.2.1 (1)潜在空間の可視化
全国合成レーダー画像については2016年1月から2020年3月までの、３時間ごとのデータを学習させました。
潜在空間を２次元としたので、入力データそれぞれに該当する値（エンコーダの出力）を平面にプロットすると以下のようになりました。色はデータの「月」を示しています。
分布では青や緑のドットつまり11月から2月の冬のデータが第２象限に集まっていてるようです。夏データを示す赤色系はそのほかの領域にばらまかれているように見えます。
#####3.2.1 (2)潜在空間からの画像生成
逆に、この２次元平面のどこかのポイントを指定すると、デコーダによって対応する画像が生成されます。
入力画像が対応していない点からも画像を生成することが出来ます。
指定するポイントを平面上で連続的に動かして生成された画像を並べると、下のような図が出来ます。
平面の縦、横それぞれについて、-3.5から3.5の間を15分割した値から画像を生成したものです。
雨のエリアが異なるのですが、左上は北海道だけが雨、左下は南西諸島だけが雨、右上は東日本全体が雨、右下は西日本全体が雨という特徴で、それぞれの間が連続的に変化しています。
この図の外周あたりに沿って指定する点を移動させて、その点に対応する合成レーダ画像を生成したものが下記です。
####3.2.2 メソ予報モデルの地上気温データの結果
#####3.2.2 (1)潜在空間の可視化
同様に潜在空間を可視化すると下のような図になりました。
色は「月」を表しています。
気温なので、季節による違いがくっきりと現れます。
#####3.2.2 (2) 潜在空間からの画像生成
同様に潜在空間から画像を連続的に生成させて並べました。
上半分がどちらかというと冬データ、下半分が夏データですが、右上の方は実データが該当していない、「存在しない気象状況」という理解できない画像領域になるわけです。
潜在空間の中で指定したポイントを、散布図の上に示しながら画像を生成したものが下記の図です。
4. 潜在空間へのマッピング
4.1 速報天気図の潜在空間へのマッピング
速報天気図は2018年9月から2020年3月までのものを用いて学習を行いました。
速報天気図は、UTC0時をのぞいて３時間ごとに作成されています。
実はオートエンコーダとしての画像の再現性という意味ではうまくいきません。
地形と経度緯度線だけは再現されますが、細かい前線記号や等圧線は全く再現されませんでした。
画像生成という観点では良くない例なのですが、それでもということで、得られた潜在空間に対して、実際の天気図がどこにマッピングされるか？ということを可視化してみました。
2018年9月から、マッピングを始めます。この動画は、毎日3時UTCのデータのプロット位置を大きめのサークルで、その日までのデータ（3時以外も含む）を小さな点でプロットしています。
2020年1月からのプロットです。散布図の方には2018年、2019年のデータも表示されています。
平均は(0,0)の位置を中心に分布していることがわかります。
5. まとめ
少し前に流行ったVAEを用いていくつかの気象データに対して適用してみました。
VAEによって
・潜在空間に対応した気象場画像の生成ができる
・実データの潜在空間でのマッピングができる
ということがわかりました。
課題は、気象予報にとって何かの役に立つだろうか？ということです（笑）。
VAEは異常検知に用いられたりしていますので、観測された、あるいは予報された気象場の「異常度」のようなものが定量化できると良いのかもしれません。