「天気図っぽい前線」は気象データのどこを見て描かれているのか？(1) Grad-CAMのセマンティックセグメンテーションへの適用
0.はじめに
天気図でよく目にする温暖前線、寒冷前線といった「前線」は気象庁の予報官の方が気象データを元に解析して描画しています。
どこに前線を描くのかは機械的に決まるものでもなく、さらに前線描画のスタイルは各国の気象機関によっても異なっております。
私は数年前から機械学習の勉強題材として、この前線描画をニューラルネットワークで行おうということをはじめました。
格子点気象データ(GPV)を可視化した画像を入力として、セマンティックセグメンテーションによって気象庁が作成した「速報天気図」の前線と同じような「天気図っぽい前線」を描画させるということです。この過程で前線の描き方を調べたりしているうちに、気象予報士の資格まで取ってしまいました。Domain Knowledgeは重要です。
前線の自動描画の話は5回＋番外編1回に分けて書きました（気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（1）(2)(3)(4)(5)[「天気図っぽい」気象前線の自動描画を日本から世界に拡張してみる〜領域拡張と地図投影]
(https://qiita.com/m-taque/items/b4ea490667d88cd10a55)）
それなりに前線を描画してくれるようになったのですが、ニューラルネットワークがいったい気象データ画像のどこを見て前線を作成しているのかを知りたいところです。
そこで、画像の推論根拠の可視化でよく使われるGrad-CAMおよびGrad-CAM++を使ってみたので、記録として纏めます。
Grad-CAM++適用の話はこちら：
天気図っぽい前線」は気象データのどこを見て描かれているのか？(2) Grad-CAM++のセマンティックセグメンテーションへの適用と実装
1.Seg-Grad-CAM
1.1 Seg-Grad-CAMとは
ニューラルネットワークが画像のどこを注視しているのか、を可視化する技法としてGrad-CAMが有名です。このGrad-CAMをセマンティックセグメンテーションに適用したSeg-Grad-CAMの論文があります。
オリジナルのGrad-CAMは画像の分類結果のOne-hotベクトルを元に最後のConvolution層への勾配を計算していきます。これに対して、Seg-Grad-CAMでは関心のあるクラスに分類されたピクセル群を元に勾配を計算していきます。
CNNのある層での、クラスCのヒートマップ（$L^c_{lm}$）は、その層のForward Propagationの出力（Activation Map）を、その層のチャネル方向に重み付き平均をとって作成されます。そのチャネルごとの重みは、Cに分類されたピクセル群の数値の、Activation Mapの値に対する勾配を、チャネルの画像全体で平均することで得られます。
言葉ではわかりづらくなってしまいますが、これを数式で表すと次のようになります。
\begin{align}
L^c_{lm} &= Relu( \sum_k \alpha^k_c A^k_{lm} ) \\
\alpha^k_c &= \frac{1}{N}\sum_{u,v} \frac{\partial Y^c}{\partial A^k_{uv}} \\
Y^c &= \sum_{i,j \in \mathcal M} y^c_{ij}
\end{align}
最後の式はクラスに所属するピクセルの値の和を取る表式ですが、論文やそれを実装したコードでは、実際には和を取っているわけではなくそのまま２次元のマップとして勾配計算を行います。
重みはピクセルが属するクラスを示す数値$y^c_{ij}$のActivation Mapでの勾配を、そのActivation Mapで平均した値を使っています。
Activation MapはCNNのForward Propagationを実行すると計算されて、Kerasにもその値を取り出す機能があります。また、勾配$\frac{\partial Y^c}{\partial A^k_{uv}}$もKerasのバックエンド関数を用いることで、
Forward Propagationさせたときの各層の情報と、$y^c_{ij}$から計算することができます。
1.2 写真画像のsemantic segmentationでの例
Kira Vinogradovaらは、ネットワークとしてU-Net(Ronnebeger, Fischer, and Brox et.al 2015)、データセットとしてCityspcapes(Cordts et al. 2016)を用いてSeg-Grad-CAMの実験を行っています。
この例はセマンティックセグメンテーションにより道路（上図）および空（下図）と分類されたエリアのヒートマップ（各図右上）です。空や道路と分類される時に、写真のどういう部分に着目しているかを表しているのですが、道路はともかく、空については木を見ているようにも思えます。
つまり、木があると空があると、少しバイアスがかかっている可能性があります。
いずれにせよ、こういった判定根拠を示すことができるわけです。
##2. 気象前線検出のネットワークへの適用
写真のセマンティックセグメンテーションについては、人間の目には空や車に属するピクセルは見えています。しかし、「天気図っぽい前線」の入力とするような気象図には、明確な正解エリアが見えているわけではありません。また１枚の写真だけではなく、複数の物理量に対応して、入力画像も数種類を使用しています。写真の処理と、気象図から前線を描画することの間にはこのような違いがあります。
とはいうものの、学習それ自体は、入力画像のどのピクセルをどういうクラス分けをするべきなのか、正解データに基づいて学習するということで、どちらもほぼ同じです。ということで、画像どのエリアに着目してクラス分けをしているのか、なんらかの可視化がなされるのでは？という
期待です。
###2.1 ネットワーク構造
前線検出に用いるニューラルネットワークの概念図です。
詳細は気象データをもとに「天気図っぽい前線」を機械学習で描いてみる（5）を参照頂ければ幸いです。
ピクセルサイズが半分になるとチャネルを倍のサイズにするのが一般的ですが、このネットワークではチャネルはおおよそ1.5倍くらいにしか増やさない構造です。これによってパラメタの数がかなり抑えられています。
図ではConv2Dを実施している箇所が青い直方体で示されています。Seg-Grad-CAMのヒートマップもこれらの箇所で計算していきます。
###2.2 入力データ
気象データは7種類、サイズは256x256ピクセル、それを128x128ピクセルに縮小したもの、さらに64x64に縮小したものを用います。各種別ごとにRGBの3チャンネルを持っています。
U-Netと似たようなスキップコネクションを有するCNNにより、最終的に4クラスのセマンティックセグメンテーションを行います。
入力データの可視化の詳細は、気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(2)に記載しています。
具体的な入力画像は下記の7種類です。
| # | 気象データ種別 | 内容 |
|---|---|---|
| 1 | 地表 風・気圧・気温 | 地表面での風ベクトル、気圧(等高線)、気温(等高線) |
| 2 | 850hPa 風・高度・気温 | 850hPa面での風ベクトル、気圧高度(等高度線)、気温(等高線） |
| 3 | 850hPa 相当温位 | 850hPa面での相当温位(等高線）、空気の暖湿度を示す |
| 4 | 700hPa 上昇流 | 700hPa面での鉛直方向の風速。対流的な空気の流れの強さを表す |
| 5 | 700hPa 湿数 | 700hPa面の空気の露点温度との温度差。降水活動と関連がある |
| 6 | 500hPa 風・高度・気温 | 500hPa面での風ベクトル、気圧高度(等高度線)、気温(等高線） |
| 7 | 300hPa 風速(等値線) | 300hPa面での風速分布。強風軸の所在を示す |
###2.3 描画データ
ニューラルネットワークの最終層は各ピクセルごとに４つのクラスへの所属確率を算出しています。
これによって前線が描画されます。
| # | 色 | 内容 |
|---|---|---|
| 0 | 白 | 前線要素に該当しない |
| 1 | 赤 | 温暖前線および停滞前線の一部 |
| 2 | 青 | 寒冷前線および停滞前線の一部 |
| 3 | ピンク | 閉塞前線 |
この他、正解画像にこれらと同じ色を使用した次のマークが含まれているため、ネットワークはそれも込みで学習してしまいます。
赤： 低気圧を示す「低」の文字、および「熱低」の文字
青： 高気圧を示す「高」の文字
ピンク：台風番号を示す文字列
正解画像を作る際に、気象庁のカラー版SPASのラスタ形式版から色を元に正解データを抽出したため、同じ色を持つこれらの文字列も抽出されてしまいます。
すべての図に出てくるわけではありませんが、たまたまそれぞれの色のマークが合わせて出てきている具体例を下に示します。
停滞前線の他に、低気圧・高気圧・熱低・台風番号が表示されています。
なお、気象庁の天気図からカラーを元に抜き出す話は下記にも記載しています。
気象データをもとに「天気図っぽい前線」を機械学習で描いてみる(3)
##3. Seg-Grad-CAMの実装
論文著者らによる実装である、Github:kiraving/SegGradCAMを参考にさせていただいています。
(1) Kerasのバックエンド関数
Kerasのバックエンド関数を使用しますので、それらをImportします。
import keras.backend.tensorflow_backend as KTF
(2) Seg-Grad-CAM計算部分
ネットワークの重みや定義は、以前に学習済みのものを呼び出しています。
Grad-CAM適用にあたり、あらためて学習し直す必要があるわけではありません。
"""
以前に学習したパラメタと重みを読み込む
"""
NN_1 = load_model( paramfile , \
custom_objects={'softMaxAxis': softMaxAxis })
NN_1.load_weights(weightfile)
NN_1.compile(optimizer='adam', loss='categorical_crossentropy' , \
metrics=['accuracy'])
"""
Seg-Grad-CAMの計算のためにForward Propagationを実行する
"""
prepro_input = [ \ # ネットワークに必要な3種類の解像度の入力
np.expand_dims(np_i_data_tr[c_datalist], 0),\
np.expand_dims(np_i_128data_tr[c_datalist],0),\
np.expand_dims(np_i_64data_tr[c_datalist] ,0) ]
np_p_data = NN_1.predict(prepro_input) # Forward Prop.実行
np_p_data = np.transpose(np_p_data,(0,2,3,1))
# このネットワーク定義は'Channels_first'であるので後の処理のために置換
w_array = np_p_data[0] # Forward Prop.の結果（予測結果）を取得
"""
勾配を計算する起点となるConvolution層:prop_from_layer
勾配を計算する対象となるConvolution層:out_layers_list
"""
prop_from_layer = "conv2d_23"
output_layers_list = [ 'conv2d_2' , 'conv2d_3' , 'conv2d_7' ,\
'conv2d_8' , 'conv2d_12' , 'conv2d_13' ,\
'conv2d_15' , 'conv2d_16' , 'conv2d_17' ,\
'conv2d_18' , 'conv2d_19' , 'conv2d_20' ,\
'conv2d_21', 'conv2d_22' \
]
勾配を計算する対象となるConvolution層は、ネットワーク定義の際に明示的に指定していればその値が使えます。指定していない場合は、NN_1.summary()で表示すると内部でどういう名称となっているか、表示されます。
max_preds = w_array.argmax(axis=-1) # 各ピクセルの所属クラスを取得
# 対象とするconvolution層の数だけ以下を実行
for prop_to_layer in output_layers_list :
row_cnt = 0
# 可視化の準備
fig, ax = plt.subplots(4, 7, figsize=(20,13.75), \
sharex="all", sharey="all" , tight_layout=True )
"""
最終層の出力取得の定義
"""
y_c_o = NN_1.get_layer(prop_from_layer).output
conv_output = NN_1.get_layer(prop_to_layer ).output
# 各クラスごとにヒートマップ（デフォルトの白は除くので1からスタート）
for cnt_cls in range(1, 4):
cls = cnt_cls
# クラスに所属すると算定されたピクセルだけ値が入るようにするマスクを
# 計算する( roi )
# さらに、その中で指定した領域のみ有効にするためのマスクも
# 作成する( roib )
roi = np.zeros((1,w_array.shape[0],w_array.shape[1]))
roi = np.round( \
w_array[:,:,cls] * ( max_preds == cls ) \
).reshape(w_array.shape[0], w_array.shape[1])
if _flag == True :
roib = np.zeros((1,w_array.shape[0],w_array.shape[1]))
roib[:, _ruc_y:_llc_y+1, _llc_x:_ruc_x+1] = 1.0
# 左下(_llc_x, _llc_y) 右上 (_ruc_x, _ruc_y)で囲まれた
# 矩形領域のみ1.0
else:
roib = np.ones((1,w_array.shape[0],w_array.shape[1]))
"""
最終層の結果についてクラスに所属する部分、かつ指定領域に所属する部分を
roi, roib を掛けることで抽出する定義（テンソル計算）
"""
y_c = y_c_o[:,:,:,:] * roi * roib
"""
conv_outputが指定しているconvolution層の勾配を計算する定義
"""
grads = KTF.gradients( y_c , conv_output )[0]
"""
ここまで定義してきたテンソルの計算を関数として集約定義
"""
gradient_function = \
KTF.function( \
[NN_1.input[0], NN_1.input[1], NN_1.input[2]], \
[conv_output, grads] \
)
"""
集約定義した関数に実際の値(prepro_input)を代入して計算実行
"""
output , grads_val = gradient_function( prepro_input )
# channels_firstのため配置転換
output_r, grads_val_r = np.transpose(output,(0,2,3,1)),\
np.transpose(grads_val, (0,2,3,1))
"""
Activation map(A_m)と勾配(grads_m)を取り出す
画像横サイズ×縦サイズ×チャンネルの次元
"""
A_m, grads_m = output_r[0,:] , grads_val_r[0,:,:,:]
"""
勾配を計算した層のそれぞれのチャネルで、勾配を平均計算
---> Activation mapを足し合わせる重みとなる
alpha_c:チャンネル数の次元
"""
alpha_c = np.mean(grads_m, axis=(0,1))
"""
重み付き足し合わせをテンソルとベクトルの内積計算で実行
---> 画像横サイズ×縦サイズの次元（1枚の図）
"""
cam_o = np.dot(A_m, alpha_c)
# Relu及びNormalizeの実行
cam_o = np.maximum(cam_o, 0)
cam_o = cam_o / cam_o.max()
# 配列からPIL画像に変換
cam = array_to_img(cam_o.reshape( \
cam_o.shape[0], cam_o.shape[1], 1)\
)
# 入力画像のサイズに拡大
cam_img = cam.resize(\
(np_i_data_tr.shape[2], np_i_data_tr.shape[3]),\
Image.LANCZOS)
こうして得られたcam_imgを入力画像にオーバーレイして表示すると、求めるヒートマップが得られます。
##4.表示例
###4.1 自動描画された前線例
実際のデータ例を見てみます。
このデータは2017年3月22日UTC0時のデータで前線描画したものです。
地上気圧図を見ると、関東の東、東経149度北緯35度付近に閉塞前線を伴う低気圧が見えており、閉塞点から南東方向に温暖前線、南から南西方向へ弧状に寒冷前線が検出されています。比較的うまく前線が検出されている例です。
また、この低気圧のすぐ上には赤い「低」のような文字が出現しています。
さらに北海道付近にひとつ、他にもいくつかの「低」マークと、朝鮮半島付近に青い「高」の文字が見られます。
この時の実際の気象図（SPAS）は下のようでした。
さきほどの自動検出結果と比べると、本物の方は寒冷前線が南西方向により長く引かれています。また閉塞点の位置もややずれているかと思われます。
低気圧マークは低気圧本体中心付近と北海道付近のものは合格圏内？その他は実際には解析されておらず、CNNは何を見て書きたくなったのか謎です。朝鮮半島付近の高気圧マークはまあ良いのではないかと思われます。
###4.2 Seg-Grad-CAMによる可視化例
この前線を描画するためにニューラルネットワークが入力データのどういったところに着目したのかをSeg-Grad-CAMにより見ていきます。
####(1)ニューラルネットワークの構造
まず、再びニューラルネットワーク概念図を示します。
各convolution層の上に書かれている数字が、conv2d_NNという層のNNを表しています。
以後の説明では、conv2d_NNで得られたヒートマップというような言い方をします。
####(2)着目したconvolution層のヒートマップ全体像
まず比較的浅い層であるconv2d_7のヒートマップを見てみます。この層は128x128のサイズで入力画像256x256に対して縦横ともに半分になっています。チャネル数は36です。
図：conv2d_7で得られたヒートマップ
上記の絵は、上段から下に向かってクラス「1(赤)」，「2(青)」，「3(ピンク)」のヒートマップを、左から右方向に7種類の入力データのそれぞれにオーバレイして表示させたものです。3✕7＝21種類の図になっています。
例えば上段の一番左側ですと、クラス「赤」のヒートマップを850hPa相当温位の図に重ねている状態です。
最下段には計算された前線と入力画像のオーバレイも表示さいています。
ヒートマップ自体のカラーは、青色から黄色を経て赤色になるほど大きな値を示しています。
####(3)入力データ種別との対応
(2)に示した図のいくつかを拡大してみます。
クラス「赤」とクラス「ピンク」では、関東の東海上の低気圧領域、そこから延びる前線を描画すべき領域を注視しているように見えます。
クラス「赤」は850hPa相当温位の気象図との対応が最もよさそうです。
低気圧中心の周辺と低気圧中心の南東側の等相当温位線が集中した領域を注視しているようです。
図：850hPa相当温位に重ね合わせた、クラス「赤」のヒートマップ
クラス「ピンク」は低気圧本体の構造を注視しているように見え、地上気圧図との対応が良さそうです。
図：地上気圧・温度・風速に重ね合わせた、クラス「ピンク」のヒートマップ
一方、クラス「青」については対応のよい種別がどれか、少し微妙です。また、赤やピンクに比べると少し全体的な構造を見ているようにも思えます。あえて言えば、等相当温位線の集中帯を挟む南北の領域を見ているようにも思えます。
図：850hPa相当温位に重ね合わせた、クラス「青」のヒートマップ
####(4)深い層でのヒートマップ
次に、もう少しネットワークの深い段階に進んだconv2d_15でのヒートマップです。
この層では画像サイズは32x32、80チャンネルになっています。
図：conv2d_15で得られたヒートマップ
クラス「青」とクラス「ピンク」はやはり低気圧を見ているようですが、クラス「赤」については何やら注視領域が増えているようです。赤い低気圧マーク「低」をどこに描くべきかを見ているのでしょうか？
さらに進んだconv2d_17およびconv2d_18のヒートマップです。ここは画像サイズが16x16x120(ch)と最も抽象度が進んでいると考えられる、U-Netライクな構造のボトルネック部分になります。
図：conv2d_17で得られたヒートマップ
図：conv2d_18で得られたヒートマップ
前線を描くための大域的な構造というよりは、どうやって「高」や「低」を描くための領域を見ているように思えます。特に、conv2d_18でのクラス「青」は日本の東にある低気圧と大陸の華中から華南にかけての図を見て、**「それらに挟まれた中間あたりに『高』を書こう」**と考えているように見えます。
4.3 出力画像の一部領域に絞ったヒートマップ
勾配計算させる際に最終結果のピクセル範囲を、興味のある範囲に限定することで、影響の分離を試みます。
先に示したソースで、roibを有効にする_llc_x, _llc_y, _ruc_x, _ruc_yを指定します。
(1)前線と高気圧記号
クラス「青」について、関東の東にある低気圧の寒冷前線の部分、朝鮮半島付近にある高気圧記号、それぞれのみを使用してヒートマップを作成しました。指定した部分は緑色の矩形で示しています。矩形領域に含まれないクラス「赤」やクラス「ピンク」はヒートマップに表現されなくなります。
図：寒冷前線部分に絞ったヒートマップ（conv2d_15で得られたヒートマップ）
図：高気圧記号部分に絞ったヒートマップ（conv2d_15で得られたヒートマップ）
それぞれの850hPaの相当温位にオーバレイさせたヒートマップです。
図：寒冷前線に絞ったヒートマップ（conv2d_15, 850hPa相当温位にオーバレイ）
図：高気圧記号に絞ったヒートマップ（conv2d_15, 850hPa相当温位にオーバレイ）
高気圧記号に絞った方は低気圧から西側の気象場の構造を全体的に捉えようとしているようにも見えます。
これに対して前線記号に絞った側は、低気圧周辺の等値線の集中しているところに注視しているようです。
####(2)低気圧に関わる前線に絞ったヒートマップ
高気圧や低気圧マークの影響を排除するため、関東東海上の低気圧の前線に絞ったヒートマップを作成してみます。
図：低気圧に関わる前線に絞ったヒートマップ(conv2d_15で得られたヒートマップ)
これを先ほどの、すべてのピクセル値を使用した前線に絞らないヒートマップと比較してみます。
図：範囲限定しないでconv2d_15で得られたヒートマップ（低気圧マークなども含む）
クラス「青」と「ピンク」はそれほど差異はありませんが、クラス「赤」の日本列島に西側から大陸にかけての広い範囲に見られる注目点がなくなっています。なくなった部分は、低気圧マークの影響であったと考えられます。
図：クラス赤：conv2d_15で得られたヒートマップ（地上気圧にオーバレイ。低気圧マークなども含む）
図：クラス赤：conv2d_15で得られたヒートマップ（地上気圧にオーバレイ。前線記号に絞った）
5.まとめ
天気図っぽい前線を描くCNNに対して、Grad-CAMをセマンティックセグメンテーションに拡張したSeg-Grad-CAMを適用しました。
もともと、写真画像のセグメンテーションのように対象がズバリ見えているわけではありません。また、クラスに分類された複数のピクセル群の影響が合成されていることもあります。こういったことから、ズバリここを見て描いているというところまで言い切るのは難しいところです。
気象的に見た場合、このようなところが言えそうです。
・前線描画は、850hPa相当温位と前線の対応がよいように思われ、入力画像と前線の対応を見てこれまで思っていたことがある程度裏付けられました。
気象予報士試験での前線描画の問題は多くの場合に850hPa相当温位図を用います。
・特に深い層で、前線というより高気圧・低気圧のマークの位置を決めるための注視が行われていているように思えます。
機械学習の観点からは、本来の目的よりも高気圧、低気圧といった記号を特定するのにかなりの層を費やしているのではないかと思われ、前線検出をメインタスクとするニューラルネットワークとしては前線のみの教師データを作る必要性があるように感じられました。
今回はGrad-CAMを用いましたが、次回はGrad-CAM++を同様にセマンティックセグメンテーションに拡張してみた話です。
天気図っぽい前線」は気象データのどこを見て描かれているのか？(2) Grad-CAM++のセマンティックセグメンテーションへの適用と実装