toyo-daichi / da_basic Goto Github PK

最も有名なトイモデルの一つであるLorenz96モデルを実装する。

 wikiの引用

Lorenz96 modelは、エドワード・ローレンツ(1996)により提唱された力学系方程式である。

Lorenz96での問題。

#36 にてLorenz96が観測の情報がフル条件であれば、
EKF・EnKFが非常に精度よく回ることが確認できた。

が、観測の条件を緩めてしまうと計算が発散する。

プログラムに何かしらのミスがあるはず。。。
アルゴリズムのミスではないので。。。

チェックする！！！

カルマンフィルタ系列の値、予測値、真値の値に整合性がない。
明らかにカルマンフィルター 系列の値が発散している。
※おそらく、状態遷移行列が悪さをしている。

予測精度評価指標のライブラリを作成する

• Root Mean Squared Error (RMSE)
• 相関係数など

SRF(Squared root Filter: 平方根フィルタ)では、アンサンブル平均を最も確からしい状態として、解析方程式を解く。

x_anl = x_sim + matmul(K, hdfx)

参考: 気象学におけるデータ同化、実装してある技術ブログ

完成したKFとEnKFでは、filter divergenceが起きずにうまく回っている。
この原因は何であるかを突き止める。

前回の作成より理解度が進んだ。
適宜、コメント等を追加してリファクタリングを行う。

【特に留意する点】

• 命名の厳密化(有識者が見れば、何の変数か識別できるように。)
• python描画コードの詳細化

lorenz63モデルがある程度の時間がかかる。
下記のスクリプトを並列化してみること。

def lorenz_3ddraw(true_data: list, sim_data: list, da_data: list, *, timestep :int='None'):

for i_num in tqdm(range(0, len(time_list[0]))):
    lorenz_3ddraw(true_list, sim_list, da_list, timestep=i_num)

前回のissue(カルマンフィルターの逆行列計算が間違っていた件)を、
行列サイズが2x2のため、公式を用いて行い計算エラーを修正した。

ただ、本来の目的はLAPACKを使用しての逆行列計算だった点ができていない。
サブルーチンでLU分解して逆行列計算する経験を行う。

#33 でのプルリクエストのマージ後、現在の状態ではバグが出ないことを確認している。

その一方で、Lorenz63/Lorenz96ともに出力形式がでたらめになっているので修正する。

誤差共分散行列の対角成分とカイ２乗分布を用いれば、きれいにかけると思われる。

【参考リンク】

• https://myenigma.hatenablog.com/entry/20140621/1403336277

カルマンフィルタ(Kalman Filter: KF)を非線形モデルに適用した場合、仮にモデルが完全でランダム過程を含まない理想的な場合でも、KFは発散することが知られている。

解析値が真の値から離れていき、データ同化が上手く働かないことを言う。filter divergence とも言う。

この一因は、誤差共分散の予報方程式の近似において紛れ込む誤差により、KFが見積もる誤差が過小となって第一推定値を信頼しすぎるためである。これを防ぐために共分散膨張(covariance inflation)を行い、より観測データを取り組むことでKFを安定化することができる。

現段階(tag: assimilation_tool_ver2.0)では、KFに共分散膨張は実装していない。結果を確認してみたところ、データ精度が著しく下がっていたため検討する。

風上側に感度がある場合、インクリメント相関は風上側にインクリメント相関のピークが寄ると思われる？

現在は、下記の対応で行っている。
...コメントアウトなので、根本的な解決ではない

oscillation_compiler.sh

boolen='true' # write output
./${prg} > ./log.${prg} << EOF
&output
    output_file  = '${outputfile}',
    !opt_veach   = '.${boolen}.'
/
EOF

最終的にはLETKFを実装したい。その過程で、KFとEnKFの実装を加える。
LETKFでは、局所化関数に風上領域の活用が加えられる方法を模索する。

プルリクエストを出すのは、KFとEnKFが実装できた段階で。。。

commonでは、SFMT(SIMD-oriented Fast Mersenne Twister)ライブラリを活用している。
ライブラリの概要説明と引用元の明記を忘れないようにREADME.mdを作成すること。

引用元(https://github.com/takemasa-miyoshi/letkf)

Rungeの手法が間違っていると思われる。