Transcript 発表資料
データ同化ワークショップ 2014/01/08 ４次元変分法を用いた 初期値とパラメータの同時最適化 伊藤耕介 海洋研究開発機構 ■謝辞：石川洋一先生、淡路敏之先生、杉浦望実様、 川畑拓矢様、本田有機様、加藤輝之様 数値予報の誤差要因 • 初期値の誤差 • 境界値の誤差 • モデルの不完全性 (必要な要素の欠如， パラメタリゼーションスキーム, 固定パラメータの値など) 今日のメッセージ 重要な固定パラメー タに不確定性がある のなら、制御すべき。 発表内容 • 4次元変分法における初期値の最適化 • 初期値と固定パラメータの同時最適化 理論的な話 台風状況下の海面交換係数の最適化 Frequently Asked Question • (余裕があれば)Hybrid EnKF-4DVAR ※簡単のため、今日の話は擾乱の時間発展に線形を仮 定し、観測演算子は線形であるとします。また、 内積が要素の単純積の和である標準内積を考えます。 4次元変分法における 初期値の最適化 (まずは固定パラメータに 誤差がないと仮定する) 典型的な４次元変分法 • 評価関数Jを定義 𝐽 1 = 2 1 + 2 𝐱0 − 𝐱b 𝑡 観測値がない所で同化結果が 第一推定値から離れないように T 𝐁 −1 𝐱0 − 𝐱b T −1 𝐲𝑡 − 𝐇𝐱 𝑡 𝐑 • 解いている問題 𝐲𝑡 − 𝐇𝐱 𝑡 観測値と同化結果が 近づくように ある固定区間(同化ウィンドウ)で、xが モデルの時間発展に従うという条件の もと、Jを最小化する初期値x0を求める 評価関数Jを最小にするx0を探索 • 評価関数J の分布は事前に分からない • 方法１: x0の値を変えて予報を繰り返す. ⇒単純明快だが変数が多いと実行不可能 • 方法２: ∂J/∂x0を求める ⇒効率的に最適なx0にたどり着く J ★★ ★ ★★ ★ ★★ ★★ ★★ ★★ ★★ ★ ★★ ★★★ ★★ ★★ ★ ★★ ★ xb xupdate xoptimal x0 ∂J/∂x0の計算方法 • xの摂動の時間発展をdδx/dt=Mδxと書く • 勾配に関し以下の式が得られる(導出は略) 𝑑 𝜕𝐽 𝑇 𝜕𝐽 =𝐌 +𝐅 𝑑(−𝑡) 𝜕𝐱 𝜕𝐱 アジョイント行列 (転置行列に等しい） 観測値とモデル結果 のミスフィット等 • 同化ウィンドウをt=0からt=T(>0)まで, ∂J/∂x=0(t=T)として、逆向きに時間積分 • これで∂J/∂x0が計算できる！ 非 線 形 モ デ ル 𝑑𝑥 = −𝜎(𝑥 − 𝑦) 𝑑𝑡 𝑑𝑦 = −𝑦 − 𝑥𝑧 + 𝑟𝑥 𝑑𝑡 𝑑𝑧 = 𝑥𝑦 − 𝑏𝑧 𝑑𝑡 x 𝑑𝐱 = 𝑀(𝐱) 𝑑𝑡 𝐱 ≡ (𝑥, 𝑦, 𝑧)𝑇 接 𝑑𝛿𝑥 = −𝜎(𝛿𝑥 − 𝛿𝑦) 線 𝑑𝑡 形 𝑑𝛿𝑦 モ 𝑑𝑡 = −𝛿𝑦 − 𝑥𝛿𝑧 − 𝑧𝛿𝑥 + 𝑟𝛿𝑥 デ 𝑑𝛿𝑧 ル 𝑑𝑡 = 𝑥𝛿𝑦 + 𝑦𝛿𝑥 − 𝑏𝛿𝑧 δxo 𝑑𝛿𝐱 = 𝐌𝛿𝐱 𝑑𝑡 −𝜎 𝐌 = −𝑧 + 𝑟 𝑦 ア 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 𝜕𝐽 ジ = −𝜎 − 𝑧 − 𝑟 + 𝑦 + 𝐹𝑥 ョ 𝑑(−𝑡) 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧 イ 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 𝜕𝐽 ン =𝜎 − + 𝑥 + 𝐹𝑦 ト 𝑑(−𝑡) 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧 モ 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 デ = −𝑥 − 𝑏 + 𝐹𝑧 ル 𝑑(−𝑡) 𝜕𝑧 𝜕𝑦 𝜕𝑧 t 𝜎 −1 𝑥 0 −𝑥 −𝑏 𝑑 𝜕𝐽 𝜕𝐽 𝑇 =𝐌 +𝐅 𝑑(−𝑡) 𝜕𝐱 𝜕𝐱 −𝜎 𝐌𝑇 = 𝜎 0 −𝑧 + 𝑟 −1 −𝑥 𝑦 𝑥 −𝑏 第一推定値からちゃんと真値に近づく 真値・同化結果 第一推定値 〇 真の初期値 〇 第一推定値の初期値 〇 同化結果の初期値 同化ウィンドウ：t=0-1, ⊿tobs=0.1 σ=10.0, r=28.0, b=8/3 x0t=10.0, y0t=14.0, z0t=24.0, x0f=11.0, y0f=13.0, z0f=25.0 初期値と固定パラ メータの同時最適化 固定パラメータにも誤差 • ここまでは、固定パラメー タに誤差がないとしていた。 • Lorenz63モデルで、固定パ ラメータrにも誤差があるの に、初期値だけを最適化す るとどうなる？ (他の固定パラメータは真値とする) 𝑟𝑡𝑟𝑢𝑒 = 28.0 𝑟𝑏 = 𝑟𝐷𝐴 = 26.0 真値 同化結果 第一推定値 〇 真の初期値 〇 第一推定値の初期値 〇 同化結果の初期値 うまくいかなかった • 初期値と固定パラメータの両方に 誤差があるのに、初期値だけを制 御してもうまくいかない。 • 固定パラメータにも誤差があるこ とをちゃんと考慮しよう。 • J を定義し直して、∂J/∂rを計算。 𝐽= 1 2 1 + 2 T 𝐱0 − 𝐱b 𝐁 𝑡 −1 𝐱0 − 𝐱b + T −1 𝐲𝑡 − 𝐇𝐱 𝑡 𝐑 𝐲𝑡 − 𝐇𝐱 𝑡 1 (𝑟−𝑟𝑏 )2 2 𝛼2 接 線 形 モ デ ル ア ジ ョ イ ン ト モ デ ル 𝑑𝛿𝑥 = −𝜎(𝛿𝑥 − 𝛿𝑦) 𝑑𝑡 𝑑𝛿𝑦 = −𝛿𝑦 − 𝑥𝛿𝑧 − 𝑧𝛿𝑥 + 𝑟𝛿𝑥 + 𝑥𝛿𝑟 𝑑𝑡 𝑑𝛿𝑧 = 𝑥𝛿𝑦 + 𝑦𝛿𝑥 − 𝑏𝛿𝑧 𝑑𝑡 𝑑𝛿𝑟 =0 𝑑𝑡 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 𝜕𝐽 = −𝜎 − 𝑧−𝑟 + 𝑦 + 𝐹𝑥 𝑑(−𝑡) 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 𝜕𝐽 =𝜎 − + 𝑥 + 𝐹𝑦 𝑑(−𝑡) 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝑑 𝜕𝐽 𝜕𝐽 𝜕𝐽 = −𝑥 − 𝑏 + 𝐹𝑧 𝑑(−𝑡) 𝜕𝑧 𝜕𝑦 𝜕𝑧 𝑑 𝜕𝐽 𝜕𝐽 =𝑥 + 𝐹𝑟 𝑑(−𝑡) 𝜕𝑟 𝜕𝑦 𝑑 𝛿𝐱 𝐌 = 𝟎 𝑑𝑡 𝛿𝐜 −𝜎 𝐌 = −𝑧 + 𝑟 𝑦 𝐒 = (0, 𝑥, 0)𝑇 𝐒 𝟎 𝜎 0 −1 −𝑥 𝑥 −𝑏 𝐜 ≡ (𝑟)𝑇 𝐱 ≡ (𝑥, 𝑦, 𝑧)𝑇 𝜕𝐽 𝑑 𝜕𝐱 𝑑(−𝑡) 𝜕𝐽 𝜕𝐜 𝑇 𝐌 = 𝑇 𝐒 𝛿𝐱 𝛿𝐜 𝟎 𝟎 −𝜎 −𝑧 + 𝑟 𝐌𝑇 = 𝜎 −1 0 −𝑥 𝐒 𝑇 = 0, 𝑥, 0 𝜕𝐽 𝜕𝐱 + 𝐅 𝜕𝐽 𝜕𝐜 𝑦 𝑥 −𝑏 ∂J/∂x, ∂J/∂y, ∂J/∂zは変更なし＆∂J/∂rは既知の∂J/∂yから計算可 今度はうまくいった！ 〇 真の初期値 〇 第一推定値の初期値 〇 同化結果の初期値 𝑟𝑡𝑟𝑢𝑒 = 28.0 𝑟𝑏 = 26.0 𝑟𝐷𝐴 = 27.7 位相空間での概念図 Jの極小値 ★ ★ ★ ★ 初期値 現実問題における固定パラメータ • 固定パラメータと言ってもいろいろ 重力加速度のように、精密に値が観測されてお り不確定性が小さいもの。 パラメタリゼーションを行うにあたり、素過程 のモデル化に現れる係数 • 後者の場合、観測と照らし合わせることになるが、 そのような観測に基づいた値は往々にして、ユニ バーサルに適用可能でなかったり、観測の困難さ から不確定性が高くなっている。 • もちろん、良いパラメタリゼーションと十分な観 測が望ましいが、それが完璧でないときには不確 定性があることも考慮に入れておいたほうがよい。 台風状況下の海面交換 係数の最適化への応用 National Hurricane Centerの予報結果 120時間予報 96時間予報 72時間予報 48時間予報 24時間予報 24時間予報 1990年 120時間予報 96時間予報 72時間予報 48時間予報 2008年 1990年 2008年 RSMC Tokyo – Typhoon Centerの予報結果 72時間予報 48時間予報 24時間予報 1982年 (沢田雅洋博士提供） 2010年 WISHE: 台風の最大風速の理論式 (Emanuel, 1986) k ～ ■ Ck: 熱交換係数 (～CE: 潜熱交換係数), CD: 摩擦係数 眼 CD 摩擦係数 波浪や波しぶきに 依存する係数 壁雲 CE 潜熱交換係数 2003年以降に提案された 摩擦係数と水蒸気交換係数の風速依存性 ・依然として、強風速域の海面交換係数は 観測が難しく、不確定性が大きい 軸対称台風モデルを用いた最大風速の 強風速域におけるCD, CEに対する依存性 100 0.0 20 0.0 ⊿CD (伊藤, 2011) 水蒸気交換係数に対する勾配 • モデル最下層の比湿の時間発展 𝑑𝑞𝑙 = ⋯ − |𝑉𝑙 |(𝑞𝑙 − 𝑞𝑠𝑒𝑎 )𝐶𝐸 𝑑𝑡 モデル最下層の比湿 水蒸気交換係数 • 水蒸気交換係数に対する勾配 𝑑 𝜕𝐽 𝜕𝐽 = − 𝑉𝑙 𝑞𝑙 − 𝑞𝑠𝑒𝑎 +𝐹 𝑑(−𝑡) 𝜕𝐶𝐸 𝜕𝑞𝑙 物理空間における模式図 4次元変分法による海面交換係数の最適化 （簡単のため，定常な移流過程だけを考えています） 軸対称台風モデルを用いた双子実験 ■ある数値実験の結果を「真値」とみなし，「誤った初 期値・交換係数」から開始した数値計算に対して観測 値を同化し、真値に近い値を復元できるか確認する． ■ 本研究で行った双子実験の構成要素 True Observation NoAsm 「真の」CD, CEと初期値を与えた実験 Trueの計算結果 + ガウシアンノイズ 「誤った」CD, CEと初期値を与えた実験 Asm_NoCoef NoAsmから出発。初期値のみを最適化 Asm_Coef NoAsmから出発。CD, CEと初期値を最適 化 双子実験の主な設定 同化ウィンドウ Day 6.0からDay 10.0 観測物理量・時間間隔 12 hごと (動径風速, 接線風速, 水蒸気量, 混合層の運動量) 観測量 r<100 km, z<5 kmの全グリッド 強風時の交換係数 [誤った値] CD, CE: Large and Pond (1981) (V≦30m/s), V>30 m/sで は一定．(誤差は近年提案された値のばらつき程度の大きさ) 背景誤差共分散行列 B 既知とする: 典型的影響半径をフーリエ解析から決め、大きさ は”TRUE”の時間変動の分散のスケール．同じ物理量のみ 観測誤差共分散行列 O 既知とする: 対角成分のみ。大きさはBの対角成分の1/20 制御変数 ◆Asm_NoCoef: 初期値 ◆Asm_Coef: 初期値 & CD & CE 最適値の探索手法 最急降下法（探索回数は40回） 4 NOAA/HRDによる台風 の集中観測の一例 0 0 r(km) 60 (Montgomery et al., 2006) CDとCEの最適化の結果 NoAsm:「誤った」 CD, CEによる実験 データ同化 ～ 擬似的に生成した ドロップゾンデ観測 CDとCEの復元値 「真の」CDとCE 誤ったC D 誤ったCD データ同化後のCD 真のCE 真のC E 真のCD 真のC D データ同化後のCE 誤ったCE 誤ったC E (Ito et al. 2010) 最大接線風速の再現性 観測値 真値からの平均二乗誤差 “NoAsm”: 7.3 m/s “Asm_NoCoef” 6.3 m/s “Asm_Coef” 2.1 m/s “NoAsm”や”Asm_NoCoef” では，台風が系統的に弱い。 “Asm_Coef” では真値に近 い。潜熱交換係数CEが大き くなり，摩擦係数CDが小さ い場合に，台風が強くなる というWISHEメカニズムと (Ito et al. 2010) 整合的． 風速・温位・水蒸気場の誤差（同化期間平均） NoAsm 陰影⊿v, ベクトル(⊿u,⊿w) AsmNoCoef 陰影⊿v, ベクトル(⊿u,⊿w) ⊿θ<-4K 壁雲域での 水蒸気不足 ⊿θ<-2K 陰影⊿v, ベクトル(⊿u,⊿w) 接線風速 の弱まり 接線風速 の弱まり コンター⊿θ, 陰影⊿q AsmCoef コンター⊿θ, 陰影⊿q ⊿θ<-4K 壁雲域での 水蒸気不足 ⊿θ<-2K コンター⊿θ, 陰影⊿q ⊿θ<-2K 非静力学メソ4次元変分法(JNoVA)への適用 • 2010年台風第14号(Chaba) • 計算機資源の制約から最盛期の1日間のみ 海面交換係数の補正量(第一推定値=1) • 摩擦係数は強風速域で小さく、水蒸気交 換係数で大きくするような修正。 摩擦係数 水蒸気交換係数 Ito et al. 2013 海面交換係数の風速依存性 • NoCoef：初期値のみを最適化 • Coef：初期値及び海面交換係数を最適化 Jobs: 第一推定値と観測値のミスフィット （主要項で全体の90%以上を占める） 評価関数 1 T J (x0 xb ) B (x0 xb ) y Hx R 1 y Hx J p J pc 2 T 1 5 0 ミスフィット項 -5 の変化率(%) Coef NoCoef -10 Jobs Jobs NoCoef -15 Jobs -20 -25 バルク係数を制御変数に加えることで， ミスフィットは5-20%程度低下した． 1 2 3 4 5 6 7 サイクル数 TC Chabaの構造の変化：軸対称モデルと整合的 Wind field [Coef] Wind field [Coef - NoCoef] 海面付近での 接線風速に大 きな差が出た Psrf, (Usrf, Vsrf) [Coef] ⊿Psrf, (⊿Usrf, ⊿Vsrf) [Coef - NoCoef] 気圧が低下し、 中心位置はベス トトラックに近く 200 km 200 km 強度と中心位置のずれの改善 • 解析場における台風強度の再現性が向上した ほか、解析時刻の中心位置のずれも改善。 50.0 ベストトラック 45.0 40.0 35.0 初期値のみ最適化 30.0 初期値と海面交換 係数の同時最適化 25.0 20.0 1 2 3 4 5 同化サイクル 6 7 ベストトラックからの位置ずれ(km) 55.0 最大風速 (m/s) 中心位置のずれ(km) 最大風速(m/s) 50.0 45.0 初期値のみ最適化 初期値と海面交換 係数の同時最適化 40.0 35.0 30.0 25.0 20.0 15.0 10.0 5.0 0.0 1 2 3 4 5 同化サイクル 6 7 台風の進路予報精度向上 初期値と海面交換係数の同時最適化を施 した結果、進路予報の精度も向上した。 FAQ FAQ 1:チューニングと同じ？ • 非常に効率的なチューニングと言える。 →5個の固定パラメータを10種類ずつ試 すと10^5回の実験が必要。本手法は既 存の4DVARシステムとほぼ同じ計算量 J ★★ ★ ★★ ★ ★★ ★★★ ★ ★★ ★★ ★ ★★ ★★★ ★★ ★★ ★ ★★ ★ xb xupdate xa x FAQ 2: 固定パラメータの不確定性を どのようにして与えればよいか？ • 物理的にあり得ない固定パラメータ にならないように、評価関数の定義 で制約をかけておくべきである。 • 台風状況下の海面交換係数について は、2003年以降に提案された値を 集め、そのばらつきを不確定性を代 表する値として用いた。 FAQ 3: 固定パラメータの値に その他の誤差を押し付ける可能性は？ • 現実的なデータ同化において、その 可能性は十分にある。 • 提案されたモデルと手に入る観測値 の範囲内で最適な固定パラメータを 探そうとしている。 • ただし、固定パラメータに誤差があ ることを全く想定せずに初期値に押 し付けてしまうよりは、一歩前進し ている。 FAQ 4: 初期値と固定パラメータ どちらの制御が大事か？ • 問題によるが、一般的に ☆相対的に長い時間スケールを扱う ☆固定パラメータの不確定性が大 ☆固定パラメータの場に対する影響 が大きい ときには固定パラメータの制御が重 要になるだろう。 FAQ 5: アジョイントモデルは難しそう • 示した通り、固定パラメータに対す る勾配は、アジョイントモデルMTの 出力結果だけを使って計算できる。 • グリーン関数法やEnKFを通した固定 パラメータ推定法など、アジョイン トモデルを使わないパラメータ推定 もある。 (Menemenlis, 2005；Ruiz et al. 2013) まとめ (Ito et al. 2010, 2013) • 気象学におけるデータ同化では、初 期値を制御するのが一般的であっ た。 • しかし、不確定性が高く、場に大き な影響を与える固定パラメータにつ いては、最適化を図るのが自然。 • その例として、台風状況下の海面交 換係数を制御変数に加えることによ り、同化システムの性能が改善する 4DVARでは、第一推定値の誤差共分散Bを 事前に仮定しておく必要がある。 第一推定値と同化結果のずれ 観測値と同化結果のずれ • 通常は、気候値的なBが使われるが、アンサンブル計算の出 力を使ってBを作ることで、よりもっともらしくなる。 It is used for B. Model Analysis Obs. EnKF Cycle 4D-Var 2011年台風12号における1点疑似観測実験 ・同化ウィンドウ：t=0-1h, Obs: U (t=0,z=20m) ・ベクトル：水平風の修正量 (⊿U, ⊿V) (t=0h) NMC-based B (Typical) 基本場 東風 Hybrid B U, Vの修正量 基本場 西風 渦の構造とあまり関係 のない修正 渦を南偏させる 水平風の修正 ■EnKFの結果は気象研究所国井勝氏に提供していただいた。