雲画像予測をチャレンジしてみた

気象も・画像処理も0知識だけど、勉強としてSignateでの雲画像予測SOTAをやってみました。1ヶ月間やったことを書いてみます。

データ確認

大量の画像ファイルの処理

Colab上で実行したいけど、画像データが相当多くて、ドライブにアップロードしたくても難しくて、zipファイルだけアップロードして、Colab上で解凍になりました。

!unzip "/content/drive/MyDrive/met_data_2016_01.zip" -d .
!unzip "/content/drive/MyDrive/met_data_2016_02.zip" -d .
!unzip "/content/drive/MyDrive/met_data_2017_01.zip" -d .
!unzip "/content/drive/MyDrive/met_data_2017_02.zip" -d .
!unzip "/content/drive/MyDrive/met_data_2018_01.zip" -d .
!unzip "/content/drive/MyDrive/met_data_2018_02.zip" -d .

!unzip "/content/drive/MyDrive/sat_image_2016_01.zip" -d .
!unzip "/content/drive/MyDrive/sat_image_2016_02.zip" -d .
!unzip "/content/drive/MyDrive/sat_image_2017_01.zip" -d .
!unzip "/content/drive/MyDrive/sat_image_2017_02.zip" -d .
!unzip "/content/drive/MyDrive/sat_image_2018.zip" -d .

これで約10分かかるけど、学習時間と比べるとそこまで長いではないです。

衛星画像ファイルの処理

pngファイルをcv2で読み込むと、numpy arrayというものが生成されます。

img = cv2.imread(file, 0)

今回の衛星画像はグレースケールなので、チャネルも一つだけで問題ないです。
またMatplotlibのimshowを使うと、画像の表示もできます。

fig, axes = plt.subplots(1, 2, figsize=(12, 8))
axes[0].imshow(img1, cmap="gray")
axes[0].set_title("Satellite Image 1")
axes[1].imshow(img2, cmap="gray")
axes[1].set_title("Satellite Image 2")
plt.show()

今回の提出ファイルは日本とその周辺だけになるので、画像のトリミングが必要です。numpy arrayなので簡単にできます。カットしたものの画像表示もできます。

img_cut = img[40:40+420, 130:130+340]
fig, axes = plt.subplots(1, 1, figsize=(12, 8))
axes.imshow(img_cut, cmap="gray")
axes.set_title("Satellite Image")
plt.show()

気象データの処理

次はgzファイルの処理、これは公式の方でも処理のfunctionが提供されたのでそのまま使います。np.fromfileを使ってnumpy arrayが生成されるので、衛星画像と同じく簡単に処理と画像の表示ができます。

気象データは少し範囲が衛星データより狭くて、アルゴリズムで補完することができるけど、実際今回の予測範囲は真ん中の部分なので、補完しなくても大きな影響がないと考えます。

気象データの特徴を確認

提供されたデータは地面から200hPaまでで、各レイヤーで提供されるデータの種類が少し違います。検証は後ほどですが、今回の衛星画像予測において、解像度は5kmで、つまりより大規模の変動を予測になるので、300hPa・200hPaの方が役に立つと推測します。

またデータの種類は、海面気圧・高度・気温・東西風・南北風・鉛直流・湿度があって、推測としては、気温・湿度・高度の辺は一番影響しそうなものだと思って、この部分は合わせて後ほど検証します。

タスクの確認と設計

インプットとアウトプット

チャレンジ自体は2年間の学習期間と1年間の評価期間になって、評価は24時間後の6時間ごとの衛星画像です。予測の時に最大インプットできる画像は96時間前の96枚画像だけど、インプットの衛星画像が多いほど、学習が難しくて、学習できる画像も少なくなるので、元コンペの解法を見ると、インプットは1画像で、プラス気象データというやり方もあります。

気象データの追加にはモデルの編集が必要なので、一旦アウトプットと同じフォーマットで24時間前・18時間前・12時間前・6時間前のインプットで、公式提供のモデルを使って精度を検証します。

モデル

まずは公式提供のConvLSTMをそのまま使います。Colab上のメモリ及び処理時間を考慮して、画像のインプットとアウトプットは1/4にします。Lossの計算はMSELossで、OptimizerはAdamWで、バッチサイズは16にして、1回の学習は10分ぐらいかかります。

データ処理

現段階はあまり工夫していなくて、

1. 欠損データはnp.zeros((1, 168, 128))で空白1チャネル画像を使う
2. 画像の数字は0-255なので、正規化で割る255で処理を行う

ちなみにオーグメンテーションはしていないです。いろんな研究結果を見ると、オーグメンテーションによる精度の改善は僅かなので、やった方が良い程度だと思うけど、モデルの改善にまず工夫したいと思います。

訓練と大規模改修

1回目の提出

画像のリサイズ・トリミングもあって、いろいろ前後処理の部分も工夫しました。正しく処理できたかを検証したいため、軽く5 epochぐらい学習して提出しました。結果としては50.66で、108位の中の82位です。ベンチマークの35.50より低いけど、方向は間違っていないことを確認できました。

2回目の提出

今度は真面目に1回学習して提出したくて、少し修正を入れました。まず提供されているソースコードの中のスケジューラーはCosineAnnealingLRなので、こちらをより長い学習の適するLambdaLRに変更します。

ソースコードの修正の時に見つけたけど、optimizer.stepとscheduler.stepがlossの計算functionにありました。怪しいと思っていて調べたところ確かに怪しいです。

1. optimizer.stepはミニバッチごとに減衰
2. scheduler.stepはepochごとに減衰

とのことで、scheduler.stepをlossの計算functionに入れると、ミニバッチごとに減衰するので、相当な減衰速度になります。

2回目の提出結果は、43.47で、まだベンチマークより低いけど、確実に改善がありました。

モデルの改修

フォーラムの解法を見て、EfficientNetを採用してみました。これに伴ってネットワークも大きくなるので、インプットの衛星画像をさらに24時間前と直前の2枚だけを使います。モデルを大きく修正して、EfficientNet + デコーダ + Conv2dみたいな構成にしました。加えて、インプットの気象データも64枚追加します。

時系列のデータだけど、モデルがすでに膨大なので、LSTMを使わず、単純に各インプットを連結しただけです。ただしそのまま連結すると衛星画像の特徴量と合わないので、1回Batch Normを行います。

結局デコーダの部分は2048層まで増えて、ColabでもT4 GPUギリギリまで動けます。1 epochの訓練でも1時間かかるので、訓練と提出もせずに、1 epoch程度の動作テストだけ実行しました。

気象データの検証

気象データを追加すると設計したけど、気象データは6層で合計34種類があるので、全部入れるつもりもなくて、一番影響しそうな種類を検証してみました。

比較対象は2017年7月以降半年間のデータで、衛星画像と気象データの関係性を計算します。行列の類似度を比較するために、RV係数を使います。

import numpy as np
from scipy.stats import special_ortho_group

def rv_coefficient(X, Y):
  X_centered = X - X.mean(axis=0)
  Y_centered = Y - Y.mean(axis=0)

  X_cov = X_centered.T @ X_centered
  Y_cov = Y_centered.T @ Y_centered

  numerator = np.trace((X_cov @ Y_cov))
  denominator = np.sqrt(np.trace((X_cov @ X_cov)) * np.trace((Y_cov @ Y_cov)))

  if denominator == 0:
    return 0
  else:
    return numerator / denominator

予想通り、やや上空の300hPa・500hPaの方が関係性が高いです。要素別で見ると、VGRD（南北風速）・VVEL（鉛直流）は明らかに関係性が弱くて、一方でRH（湿度）・UGRD（東西風速）は関係性が強いです。

これで、モデルに追加する気象データを絞ることができます。さらに、各気象データの数値が違うので、数値の標準化も適用します。

モデルの再改修

EfficientNetを諦める

上記で作成したモデルがやはり大きすぎて、訓練できないほど時間がかかるので、本来のConvLSTM構成に戻ります。そして他に試した手法としては、

アウトプットの画像枚数を変更、3回目の訓練は、アウトプット画像が1枚目と4枚目だけで、外挿で2枚目と3枚目を生成します。アウトプット画像を4枚にすると、本当に精度が落ちないかを検証します。（結果は一気に4枚生成の方が良い）

1時間ごとの衛星画像は本来それに対応する1時間ごとの気象データが必要だけど、精度を少し犠牲して、00hr ~ 02hrも同じ00hrの気象データを使うと、訓練データ量が3倍に戻るので、全体精度が改善できると推測します。（結果は少しだけ精度UP）

欠損の画像について、前後の画像を使って補完します。欠損の比率がそもそも高くないので、どこまで影響するのかはわからないけど、悪くなることがないと思います。

半精度に変更して速度を早くしたいけど、T4のFP32処理速度が予想より遅くて、処理時間がほぼ改善できなかったです。

最後のスコアと感想

最後に何回か回して提出した結果、Publicスコアは28.97で、29位となりました。無事ブロンズメダル圏内に到達したけど、そこまで良い成績でもないです。ただし画像データの訓練において、GPUリソースがないと大型モデルの構築と訓練ができないので、Colabを使う時点でリソースが限られます。あとは引き続き画像処理の知識を強化しないといけなくて、また他の関連チャレンジを探してやってみます。