２．データ同化法の変遷	

データ同化の歴史（１）	
データ同化は数値天気予報(NWP)から始まった。
 1922 Richarsonが手計算で数値天気予報を試みたが失敗（145hPa/6hの気圧変化を予測）	
 1950 Charneyらが電子計算機ENIACにより数値予報実験に成功（1日予報に1日要した）
初期値：主観（手書き）解析から客観（計算機）解析へ 1954 スウェーデン気象水理局が数値予報業務開始
 1959 気象庁が数値予報業務開始（世界で3番目）
客観解析法：関数あてはめ法、逐次修正法、最適内挿法（後二者では、モデル予測値を第一推定値に採用）	

データ同化の歴史（２）	
非定時観測データへの対応• 1969 Charneyらがデータ同化を提唱	
• 1979 全球大気研究計画第1回全球実験（GARP/FGGE)⇒ 全球大気観測システムの確立
データ同化の発展• 1980年代〜 変分法やカルマンフィルタの研究• 1990年代〜 全球大気長期再解析の実施
他の分野では• 1980年代〜
海洋データ同化の研究• 1990年代後半〜 大気微量成分データ同化の研究	Richardsonの数値予報実験（1910年5月20日07UTC）	

（Richardson，1922）	

Charneyらによる最初の数値予報	

最初の数値予報に使われた計算機、ENIAC。
Platzman et al. (1979)より 	

最初の数値予報に使われたモデルの格子配置。
Charney et al. (1950) より 	

ECMWFにおける500hPa高度の予報精度の推移	

・衛星輝度温度データ同化
・3次元変分法導入
・4次元変分法導入
・モデル物理過程改良など
・同化システムの誤差統計量改訂など	

(Simmons and Hollingsworth, 2002)	

全球数値予報の精度の推移	

3DVar	

4DVar	
3DVar	
3DVar	 4DVar	
4DVar	

関数あてはめ法	
Weierstrassの近似定理
閉区間の連続関数は、多項式によって任意の精度で近似できる。	
Panofsky (1949), Gilchrist and Cressman (1954)	
気象要素の空間分布を有限個の基底関数の線形結合で表し，その係数を観測データから最小二乗法で求める。解析領域を小領域に分割し、基底関数として多項式を用いるなどの方法がある。地衡風バ
ランスや傾度風バランスを考慮することもできる。
ただし、予測値を用いないのでデータ同化とはいえない。 	

問題点• 観測データの分布に大きな地域的偏りがあると、基底関数の個数を増やして観測データに近づけるほど、データ空白域の解析場が不自然になる。
•基底関数の個数があまりに少ないと、観測データから得られる小さなスケールの現象が表現されない。
• 基底関数の選択と観測データの分布によっては，係数を求めるための逆行列計算に計算不安定が生じる。	

逐次修正法	Bergthorsson and Doos（1955）	
ある範囲内の観測データに経験的な重みを与えて、各格子点の予測値 を修正する。範囲を縮小させながらこの修正を反復し、解析値とする。	

xia ( 0)xia (k )= xi

f= xia ( k −1)J∑w+( k −1)(o(rij ) x j − H j (xj =1

Ja ( k −1))

( k −1)2(rij ) + λ

j =1○(k = 1, 2, 3, )w(r ) ：重み（ r ：格子点と観測点の距離）
H ：格子点値を観測点に内挿する演算子
λ ：観測値と予測値の誤差標準偏差の比	

○

○

○

○

○
×○

○

○

×○

∑w

○

○

×
格子点値 i○●○○rij
R
×観測データ j
○
○
○
○
××
○
○
○
○

Barnes（1964）	

Cressman（1959）	
⎧ R 2 − r 2⎪w(r ) = ⎨ R 2 + r 2⎪⎩ 0)○(r ≤ R)(r > R)⎛ r 2 ⎞
w(r ) = exp⎜⎜ − 2 ⎟⎟⎝ 2 R ⎠

最適内挿法（１）	
Eliassen (1954), Gandin (1963)	
観測データの重みを経験的に与えるのではなく、解析誤差標準偏差が最小になるように統計的に計算する。	
aJf(oxi = xi + ∑ wij x j − H j (x f )j =1)

解析誤差	
Jafof
Δxi = Δxi + ∑ wij Δx j − ΔH j
aj =1
affoo(○○○○○×○○○

×
格子点値 i●○○×○)○

Δxi := xi − xiΔxi := xi − xi
Δx j := x j − H j (x)f○f

ΔH j := H j (x ) − H j (x)

○

○

○○○×○

○

×観測データ j○○×○○

解析誤差分散	J(σ ) = (σ ) − 2∑ w Δx ΔHa 2if 2ifij

j =1ifjJJ(o+ ∑∑ wij wik Δx j Δxk
j =1 k =1of+ ΔH j ΔH kf)

最適内挿法（２）	
解析誤差分散の期待値が最小になるように，重み wij を決める。	
J∑ ( Δx Δxojokf+ ΔH j ΔH kk =1f) w = Δx ΔHfikifj( j = 1, 2,, J )oσ
観測誤差の間には相関がないものとし，それらの標準偏差を j とする。
f
予測値の誤差標準偏差は水平面内で局所的に一様でσ とし，2点間の予測誤差の間の相関係数をρ ij とする。	
J∑ (λ λ δ + ρ )w = ρjkjkjkikij( j = 1, 2,, J )k =1σ jo
λ j := fσ

最適内挿法は当初は水平2次元内挿であったが，その後，空間3次元内挿に拡張された。また，地衡風平衡に基づいて，風とジオポテンシャル高度を同時に解析する多変量解析が導入され，1970年代にはいると次第に主流の客観解析法になった。Gustafsson（1981）の第5.1表によれば，当時の10の現業数値予報センターのうち，7つのセンターで最適内挿法が採用されていた。	

多変量解析（１）	
地衡風平衡などを利用して、風の観測データからジオポテンシャル高度の推定、あるいはその逆の推定などを行う。例えば、格子点jにおける東西風の観測データから、格子点iの高度を推定すると	
afzi = zi +σ ( z ) f σ (u ) f(( zu )ofρu−uijjj(u ) f 2(u ) o 2(σ ) + (σ ))必要な誤差統計量は、地衡風平衡が成り立つとして高度の誤差統計量から計算する。	

σ ( z ) f σ (u ) f ρ12 ( zu ) = Δz f ( x1 , y1 )Δu f ( x2 , y2 )g
∂Δz ff=−Δz ( x1 , y1 )( x2 , y 2 )f0∂y2=−g ∂Δz f ( x1 , y1 )Δz f ( x2 , y2 )f 0 ∂y2=−g (z) f 2 ∂σρ12 ( zz )
f0
∂y2()

多変量解析（２）	

逆符号	

逆符号	

（多田，1997）	

観測データの代表性誤差	
観測誤差の分解	o(ot) (tΔx j = x j − x j + x j − H j (x)測定誤差状態変数	)（ j ：観測点の位置）	

代表性誤差 	観測値	現実大気の状態	

数値モデルの格子点値	

現実大気に対応するモデル大気の状態	 
（格子点値の内挿）	
空間	
数値モデルの格子点	誤差統計量の推定（１）	
D値（Departure）	
d := xo − H(x f ) = Δxo − ΔH fddToo Tf T( ) + ΔH (ΔH )= Δx Δxf

仮定• 予測誤差は水平面内で局所的に一様である。
• 予測誤差の相関係数が2点間の距離のみの関数である。
• 観測誤差が互いに無相関である。
このとき、必要な誤差統計量はD値統計から次のように算出できる。	

(σ ) + (σ ) = do 2f 2j2jf 2(σ ) ρ (r ) = d d(σ ) = lim d djkjf 2r jk →0jkk( j ≠ k)誤差統計量の推定（２）	
同じ種類の観測データだけでD値統計を取れば、観測誤差標準偏差は一定とみなせるので、次の式が成り立つ。	

d j dkdj2

f 2(σ ) ρ (r )=(σ ) + (σ )
o 2

f 2jk

( j ≠ k)（多田，1997）	

最適内挿法の問題点	
• 空間内挿に基づいた方法であるため、数値モデルの状態変数と同じ物理量の観測データしか同化できない。したがって、衛星による輝度温度データなどは、逆解析によって気温や湿度にあらかじめ変換する必要があるが、それによってさらなる仮定が導入される。
• 予測誤差の統計量には長期間の平均的な値が使われるため、それが循環場の状態や観測データの分布によって変動することが考慮されない。
• 観測データの分布が解析時刻ごとに大きく変動しないことを仮定しているため、衛星観測などの非定時の観測データを適切に扱えない。
• 計算上の問題から、各格子点ごとに同化する観測データの範囲が狭く、そのために解析値が力学的にバランスすることが保証されない。
これらの問題点を解消するものとして、4次元変分法やアンサンブルカルマンフィルタなどの高度なデータ同化法が導入された。	

変分法の最適内挿法に対する利点	
３次元変分法	・数値モデルの状態変数から算出可能な任意の物理量の観測データを同化できる	
・すべての解析値を一括して算出するため、力学的なバランスに優れた解析値が得られる	

４次元変分法（上記に加えて）	
・数値モデルによって表現される物理法則をデータ同化に直接利用できる	
・非定時観測データをその時刻のデータとして同化できる	解析値の力学的バランス	• 現実の大気や海洋の運動を散逸カオス系と見なせるとすると、外力の急激な変動がない限り、状態変数で張られる相空間において、それらの状態はアトラクターのごく近傍を運動していると考えられる。これを力学的にバランスした状態という。	
• それらの推定値である解析値が力学的にバランスしていないと、解析値を初期値にして数値モデルを時間積分したときに、現実には見られない過渡的な運動が発生する。	
• 数値モデルが大きなバイアスを持っているためにそのアトラクターが現実のアトラクターからずれていると、解析値の力学的バランスを保証するうえで困難が生じる。

初期値化（Initialization）	
数値天気予報は、当初は順圧モデルや準地衡風モデルなど、重力波を含まない数値モデルを用いていたが、1960年代からより精度の高いプリミティブ方程式モデルを採用し始めた。このため，解析値に含まれる重力波ノイズを効果的に取り除く方法が必要とされた。これを初期値化とよぶ。	

気象庁領域モデルのある格子点における地上気圧の初期値からの時間変化。非線形ノーマルモード・イ
ニシャリゼーションを行わなかった場合（細実線）と行った場合（太実線：カットオフ周期12時間，破線：6時間）。横軸と縦軸の単位は，それぞれhourとhPa。
（Takano et al., 1990）	
松野・島崎 (1981）	

非線形 Normal Mode Initialization	
Machenhauer (1977)
数値モデルの方程式を高周波重力波成分に射影する。	
dxG= Lx G (t ) + N(xG (t ), xG (t ), t )dt
ここで、 x G と x G はそれぞれ高周波重力波成分とそれ以外の成分、L は正則な定数行列、右辺第2項は非線形項など残りの成分を表す。	
高周波重力波の発生を抑えるために、解析値に含まれる xG (0) を次*式を満たす x G (0) で置き換える。	
Lx*G (0) + N(x*G (0), x G (0), 0) = 0

この非線形方程式を次のように反復法で解く。	
( 0)

x*G (0) = 0(k )x*G (0) = −L−1N(x*G( k −1)(0), xG (0))(k = 1, 2, )デジタルフィルタと IAU	
デジタルフィルタ	
Lynch and Huang (1992), Lynch (1997)
デジタルフィルタを利用すると、高周波成分をより効果的に除去できる。	
*Nx (0) = ∑ hn x(nΔt )n=− NΔt ：数値モデルのタイムステップhn ：適当なローパスフィルタの重み係数

データ同化では、Lanczos filter やChebyshev filterがよく利用される。	

IAU (Intermittent Analysis Update）	
Bloom et al (1996)	
観測データによる予測値の修正を、解析時刻の予測値にいきなり加えるのではなく、時間積分しながら少しずつ加えることによって力学的なバランスを保証する。	

ナッジング（Nudging）（１）	
Hoke and Anthes (1976)	
衛星やレーダーなどによる時間的に連続な観測データを同化する経験的な方法。観測データのある格子点で、観測データに近づけるための強制項を数値モデルに追加して時間積分する。	
dxdxˆ
= F(x, t ) ⇒= F(xˆ , t ) + γ x o − xˆ
dtdt()(γ > 0)数値モデルが線形で１変数の場合に、最適な γ の値を求めてみる。	
dx=α x(α < 0)数値モデル：	
dtdx
ε t = 0, ε t ε t ʹ′ = Qδ (t − t ʹ′)
現実システム：	 = α x + ε t ,dt
観測データ：	
ナッジング：	

xo = x + δ t ,δ t = 0, δ tδ tʹ′ = Rδ (t − t ʹ′)dxˆ
= α xˆ + γ (x o − xˆ ),xˆ (0) = x(0) + ηdt

ナッジング（２）	
真値	
tx(t ) = x(0)eαt + ∫ ετ eα (t −τ ) dτ0

解析値	t(α −γ ) txˆ (t ) = xˆ (0)e+ γ ∫ x o (τ )e (α −γ )( t −τ ) dτ0= x(t ) + η e(α −γ ) tt+ ∫ (γδτ − ε τ )e (α −γ )( t −τ ) dτ0

解析誤差分散	22

σ (t ) := (xˆ (t ) − x(t ) )γ 2 R + Q ⎡ 2 γ 2 R + Q ⎤ −2(γ −α )t=+ ⎢ η −⎥ e2(γ − α ) ⎣2(γ − α )⎦

γ 2R + Q((γ − α )t → ∞ )→2(γ − α )これを最小にする γ の値と、そのときの解析誤差分散	

γm = α2 +σ m2 = γ m R

Q+αRナッジング（３）	
連続時間カルマンフィルタ： 線形システムに対する最適推定を与える。	
dx
T= F(t )x + ε t ,ε t = 0, ε t (ε t ʹ′ ) = Q(t )δ (t − t ʹ′)
dtTx o = x + δt ,
δt = 0, δt (δt ʹ′ ) = R (t )δ (t − t ʹ′)
解析値と解析誤差共分散行列	
dxˆ
= F(t )xˆ + PR (t ) −1 (x o − xˆ )dt
dP= F(t )P + PF (t ) T + Q(t ) − PR (t ) −1 P
dt
例題の１変数の場合はF, Q, R が定数であることから、定常解では	⎛ 2 Q⎞σ4220 = 2ασ + Q −∴σ = ⎜⎜ α + + α ⎟⎟ RRR⎝⎠
最適な γ を用いたナッジングによる解析値と一致する。	
多変数の場合には、ナッジングを連続時間カルマンフィルタのように最適化することは一般に困難である。