発表資料

Report
データ同化ワークショップ
2014/01/08
4次元変分法を用いた
初期値とパラメータの同時最適化
伊藤耕介
海洋研究開発機構
■謝辞:石川洋一先生、淡路敏之先生、杉浦望実様、
川畑拓矢様、本田有機様、加藤輝之様
数値予報の誤差要因
• 初期値の誤差
• 境界値の誤差
• モデルの不完全性
(必要な要素の欠如,
パラメタリゼーションスキーム,
固定パラメータの値など)
今日のメッセージ
重要な固定パラメー
タに不確定性がある
のなら、制御すべき。
発表内容
• 4次元変分法における初期値の最適化
• 初期値と固定パラメータの同時最適化
理論的な話
台風状況下の海面交換係数の最適化
Frequently Asked Question
• (余裕があれば)Hybrid EnKF-4DVAR
※簡単のため、今日の話は擾乱の時間発展に線形を仮
定し、観測演算子は線形であるとします。また、
内積が要素の単純積の和である標準内積を考えます。
4次元変分法における
初期値の最適化
(まずは固定パラメータに
誤差がないと仮定する)
典型的な4次元変分法
• 評価関数Jを定義

1
=
2
1
+
2
0 − b

観測値がない所で同化結果が
第一推定値から離れないように
T  −1
0 − b
T −1
 −   
• 解いている問題
 −  
観測値と同化結果が
近づくように
ある固定区間(同化ウィンドウ)で、xが
モデルの時間発展に従うという条件の
もと、Jを最小化する初期値x0を求める
評価関数Jを最小にするx0を探索
• 評価関数J の分布は事前に分からない
• 方法1: x0の値を変えて予報を繰り返す.
⇒単純明快だが変数が多いと実行不可能
• 方法2: ∂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の結果は気象研究所国井勝氏に提供していただいた。

similar documents