基本情報(sensitivity_toolについて)

• アンサンブル手法に基づく簡易予報感度解析(Enomoto et al. 2015; 榎本ほか 2014; Matsueda et al. 2011)
• 新しいTIGGEデータについて

基本となる考え方

場の状態ベクトルxの時間発展が非線形モデルM(x)で記述されたとする。初期時刻t=0における摂動y(i=1:m)を与えた、メンバー数mのアンサンブル予報において、時刻tにおける擾乱は次のように表せる。

z(i=1:m,t) = M(x+y(i=1:m),t) - M(x,t)

初期摂動の線形発展を仮定すると、"感度解析"は「検証領域において検証時刻に成長を表すメンバーの線形組み合わせの最適な係数p(i=1:m)を求めること」である。

z(optimum) = p(i=1)z(i=1) + p(i=２)z(i=2) +　... + p(i=m)z(i=m)
p(i=1:m)   = (p(i=1), p(i=2), ..., p(i=m)).T

同じ係数を用いて、発達する摂動に対応する初期摂動が求まる。

y(optimum) = p(i=1)y(i=1) + p(i=２)y(i=2) +　... + p(i=m)y(i=m)
p(i=1:m)   = (p(i=1), p(i=2), ..., p(i=m)).T

コードの紹介

各コードの後半に着く名称は次のことを表している。

• TIGGE: 松枝先生が取得しているTIGGEデータから作成する
• RISH_SERVER: 京都大学のサーバーからアンサンブルデータを取得して作成する

スプレッド解析

• anl_spread.py

アンサンブル随伴感度解析(Ensemble adjoint sensitivity analysis; EnASA)

• anl_EnASA_rate.py

z(i=1:m)を用いてトータルエネルギーノルムnorm(i=1:m) を計算する。各メンバーから求めたnorm(i=1:m)を用いて最適な係数p(i=1:m)を求める。

for i in range(m):
  p(i) = norm(i)/np.sum(norm(1:m))

アンサンブル特異ベクトル感度解析(Ensemble singular vector analysis; EnSVSA)

• anl_EnSVSA_mode_eigen.py: 固有値問題のアプローチ
• ./test/anl_EnSVSA_mode_svds.py, ./test/anl_EnSVSA_part_make_svd.py:特異値分解のアプローチ(上記の固有値問題で解けたので図化まで作成していない)

アンサンブル特異ベクトル法では、共分散(y.T)G(y) = (p.T Y.T)G(Y p) = 1の条件のもとで検証時刻における検証領域における擾乱(p.T Z.T) H (Z)を最大化するp(i=1:m)を求める。この問題ではラグランジュ関数からの微分から固有値問題を得る。

ラグランジュ 関数を次のように定義する。

pについて微分すると、

したがって、下記の固有値問題に置き換えることができる。

はじめに、Y.T G Yは正規直交関数で対角行列になるので、Y.T G Yは考えなくて良い。Z.T G Zについては2つの解法が考えられる。

1. Z.T G Zの固有値問題として解く。 Z.T G Zの行列サイズを確認してみると、(m, dims) (dims, dims) (dims, m) = (m, m)とメンバー数m[~O(10)]と等しくなり、簡単に固有値問題を解くことができる。

2. Zの特異値問題として解く。
   Z.T G Zの固有値問題を解く方法には、Zの特異値問題と置き換えることができる。Zは次のように分解することができる。

左特異値ベクトルUが共分散Z Z.Tの固有ベクトル、右特異値ベクトルV.Tが共分散Z.T Zの正規化された主成分の固有値ベクトルを表す。行列sigmaの対角成分は特異値である。この時、共分散Z.T Zの正規化された主成分の固有値ベクトルがp(i=1:m)に相当する。

1,2で作成したp(i=1:m)を下記のように初期場にかけて感度領域を作成する。固有値ベクトルの行列サイズが(m, m)であることに注意する。

摂動が正負の組で作成されている場合の対処方法

今回はモードを10までに指定して作成する手法にした。※正負に分けるのは下処理を加える必要があったため

• メンバー数を半減して、正負の組を片方のみにする。 (Matsueda et al. 2011)
• モード数をメンバー数以下にして、固有ベクトルを作成する。Z.T Zの固有ベクトルの行列サイズ(m, m)はなので、モード抽出では(m, mode)とする。

(参考)一般的な特異値抽出 モード別の一般的な抽出方法は、次のように行う。(参考サイト)

検証領域を考慮する場合

1. 検証領域における検証時刻の摂動z(i=1:m)から抽出したextraction_z（dims=1:ndims,i=1:m）を作成する。
2. extraction_z（dims=1:ndims,i=1:m）を用いて、固有値・特異値問題を解き、p(i=1:m)を取得する。

計算時の注意点

• 摂動はアンサンブル平均からの差ではなく、コントロールランからの差で構成される。
• ノルムの計算は検証領域によって変化することが知られている。

データ取得※

data_get.sh : grib形式のデータを取得するコード(※1)

京都大学生存圏研究所（RISH: Research Institute for Sustainable Humanosphere）から取得することができる。詳しい格子情報や配信時間等は気象行支援センターに記載されている。

data_grib2bin.sh, data_encode.py : データ形式をGrads形式に変更するコード(※2)

データにおける注意点

※1. データ取得サイト利用時は、下記の点の注意が必要。

教育研究機関向けにデータを提供しています。企業活動等のためにデータを頻繁に必要とされる方は、気象業務支援センターからデータを直接購入し、データ提供スキーム全体の維持発展にご協力ください。（サイトより引用）

※2. 今回使用しているデータは週間予報アンサンブルGPVである。

2006030100以前はgrib形式、それ以降はgrib2形式で提供されている。これに伴いメンバー数が(25→27), 提供データの要素が一部変更されているので注意が必要である。 なお、 2020(令和2)年3月24日（火）以降のデータ提供は終了している。

2020(令和2)年3月24日（火）を以て提供を終了しました。 ※最終提供プロダクトは、2020(令和2)年3月24日（火）12UTC初期値です。高分解能全球域GPV、高分解能日本域GPVをご利用ください。(サイトより引用)

作成情報

• 制作開始日　2020年7月5日
• TIGGEでの制作開始日　2020年8月17日

TIGGE増長分について

下記のスクリプトを増長した。

(1)data_grib2bin_TIGGE.sh
Grib形式のTIGGEデータをbig endianのGrads形式に変更する。
(2)data_encode_TIGGE.py
メンバーごとのディレクトリを作成して、実行しやすいフォーマットにデータを変更する。
(3)anl_(spread/EnASA/EnSVSA)_TIGGE.py
各種手法の場合を増長した。TIGGEデータは乾燥エネルギーノルム/湿潤エネルギーノルムどちらも対応している。

今後の計画

1. 参考文献の再現実験を行う。 -> OK
2. 近年の大気場でも同様の実験で行ってみる。 -> OK
3. 湿潤エネルギーノルムの導入を考える。 -> OK

参考文献

Enomoto, T., S. Yamane, and W. Ohfuchi, 2015: Simple sensitivity analysis using ensemble forecasts. J. Meteor. Soc. Japan, 93, 199- 213.
榎本剛, 山根省三, 大淵済, 2014: アンサンブル手法に基づく簡易予報感度解析. 京都大学防災研究所年報, 57(B), 163-168.
Matsueda, M., M. Kyouda, Z. Toth, H. L. Tanaka, and T. Tsuyuki, 2011: Predictability of an atmospheric blocking event that occurred on 15 December 2005. Mon. Wea. Rev., 139, 957-974.