GFDオンラインセミナー
2020/11/9

4次元データ同化による解析インクリメントの構造はどう決まるのか？
伊藤耕介 (琉球大学)

データ同化とは何か？

イノベーション d(観測値 – 背景値)

解析インクリメント z(解析値 – 背景値)

誤差共分散の推定値

例: 台風の中で中心気圧の観測だけが得られた
解析値：更新された気圧
中心気圧と風速の共分散情報を追加★

★

観測:960hPa, 背景値:980hPa→イノベーション:-20hPa

解析値：更新された風2

本解説の執筆の経緯• 2009年：指導教員だった石川洋一先生
から最適内挿法における解析インクリメントと固有値の関係に関するメモを渡される。
• 2012年：データ同化夏の学校＠ Les
HouchesでMarc Bocquet先生(Ecole desPonts ParisTech)が講演でデータ同化を逆問題を対比的に論じた。

⇒”Advanced Data Assimilation for Geosciences”(Oxford press)の中の一編として公刊されている。
• 日本語の資料が無いように感じられたので、データ同化夏の学校＠むつで紹介したのち、改めて藤井さんに協力いただき「ながれ」の解説としてまとめた。

(p175左段)「式(29)と対比すると，E=0が条件と・・・正当化する理由は見当たらない」⇒私の誤解に基づく記述ですので削除願います。

本解説について• 観測値を同化したときに得られるインクリメントの構造について深く理解することで、
原理的に解けない問題を避けられる有効な観測システムの設計にも役に立つ
GFDや逆問題に詳しい方がデータ同化に興味を持つきっかけになると嬉しい• 4次元変分法、カルマンフィルタ系列の手法といったアルゴリズムの詳細には立ち入らず、B, H, M, Rは既に与えられているとする研究というより、教科書的な内容です参考文献• Bocquet, M., 2012: An introduction to inverse modelling andparameter estimation for atmosphere and ocean sciences Advanced
data assimilation for geosciences (Oxford university press, 2012)
461–493.
• Johnson, C., Hoskins, B. J., and Nichols N. K., 2005: A singular
vector perspective of 4D-Var: Filtering and interpolation, Q. J. R.
Meteorol. Soc., 131, 1–19.
• メンケ著, 柳谷俊, 塚田和彦訳, 1997: 離散インバース理論, 古今書院, 295pp.
• Tarantola, A., 2005: Inverse problem theory and methods for model
parameter estimation, Society for Industrial and Applied
Mathematics, 352pp.

データ同化をこれから学ぶ人へ• データ同化が初めて学ぶという方。
⇒【公開版】データ同化夏の学校 基礎講義というYouTubeの動画でゆっくり解説していますので興味がありましたらどうぞ。
• 実際のアルゴリズムに触れてみたいという方⇒私のHPの配布資料に「データ同化実習テキスト＆Fortranプログラム」を置いています。
4DVarやEnKFの気持ちが分かると思います。
• より詳しく専門家と話しながら学びたいという方には、毎年８月のでデータ同化夏の学校＠むつ（今年はオンライン）をお勧めします。

問題設定・簡単化• 時刻tにおける観測値と背景値から時刻0の状態変数を推定。
• 以下のような簡単化を施します。
 逆問題のうち，簡単な線形離散逆問題(状態変数の時間発展について非線形性が残っていてもよいが, 摂動の時間発展については線形性を仮定する) について考える．
 観測誤差，第一推定値の誤差のふるまいはそれぞれ正規分布に従う
 背景誤差共分散行列と観測誤差共分散行列は事前に与えられている．
 背景誤差共分散行列は正定値対称行列または半正定値対称行列，観測誤差共分散行列は正定値対称行列とする．
 観測行列は線形とする．
 状態変数は時刻0 から時刻t までの時間発展を考えるが，観測値は時刻t だけで得られるとする(途中の観測がある場
合でも観測行列の拡張で対応可(Johnson et al. 2005))

解説のエッセンス(1/2)• 最良線形不偏推定量（見慣れたデータ同化の解）==解析インクリメント

イノベーション

• 正規化(行列を特異値分解()して、可観測)(特異値λi, 左特異ベクトルui, 右特異ベクトルvi)

解説のエッセンス(2/2)• 4次元データ同化の解は逆問題の解の一種とみなせる。
• その逆問題とは、可観測行列（モデル）と正規化された時刻tのイノベーション（出力）から、正規化された時刻0のインクリメント（入力）を求める問題である。
• 可観測行列は正規化・摂動の時間発展・モデル変数から観測変数へのマッピングを担う行列である。
• 可観測行列の特異値分解で、観測空間上の左特異ベクトルとイノベーションの内積に応じて、右特異ベクトル方向の成分がインクリメントに加わることがわかる。ただし、特異値が１より小さい場合にはノイズが大きいとみなし、寄与を小さくする。
• この解の形より、以下のことなどが分かる。
 インクリメントの構成に使われる特異ベクトルの数は、せいぜいmin(rank(H), rank(B), rank(R), rank(M))である。
 不安定モードに対応する構造は取り出されやすい。
 RとBそのものの大きさより、両者の比が重要である。

解説の流れ• 平方根行列・特異値分解・逆問題のおさらい• ４次元データ同化の解の導出• ４次元データ同化の解の特異値分解• ４次元データ同化の解の形からわかること

平方根行列・特異値分解・逆問題のおさらい平方根行列• m×mの誤差共分散行列は半正定値対称行列• 半正定値対称行列(全ての固有値が非負)は、その性質から以下の通り固有値分解できる

A=CΛCT• ここで正の固有値r個の平方根を収めた対角行−1/2
ˆ
列とその逆数を収めた対角行列 Λ r , r個の正の固有値に対応する固有ベクトルからなるm×rの行列Crを用いて、以下のようにA1/2とA-1/2を定義できる。

A1/2 =A-1/2

特異値分解• 任意のn×mの行列GをUDVTと分解できる。
• U=(u1,u2,...,un)でuiは左特異ベクトル。
• V=(v1,v2,...,vm)でviは右特異ベクトル。
• Dは対角行列 を用いてと表せる。 の要素は特異値と呼ばれ、本解説では簡単のため値が大きい順に並べる。

特異値分解（続き）• 方程式d=GzについてG=UDVTとする• uiとvi はそれぞれdとzの正規直交基底を成す

説明のため、p番目までのベクトルが張る空間をp空間と呼ぶ。

• 正の特異値がp個(p≦min(m,n))のとき，p+1番目以降のベクトルが張る空間をヌル空間という• uiとvi はGを通じて互いにつながっている順問題・逆問題とは？

• 定義
入力変数ベクトル z
出力変数ベクトル dモデル(物理モデルとは限らない) G
残差 e
• 順問題
• 逆問題
• 後述するが、4次元データ同化は「可観測行列G」と「時刻tの正規化したイノベーションd」から「時刻0の正規化したインクリメントz」を求る逆問題を設定し「ダンプ付き最小二乗解」と呼ばれる劣決定問題の解を与えることに相当する。

優決定・劣決定問題• 逆問題• 優決定問題：出力から入力を完全に決定しうるモデルを用いている逆問題。
• 劣決定問題：優決定問題でない逆問題。
dの長さがzの長さに比べて不足している。
（データ同化ではこちらを想定する）
dの長さはzの長さに比べて十分だが、解の構造を完全には復元できない。

優決定問題の解：最小二乗解• 逆問題• 最小二乗法の最小化→• 補足：優決定問題ならば、は存在する。
(証明)優決定問題ならのviの特異値に0がない。このとき、以下のように書ける。

劣決定問題の解• 劣決定問題：逆問題
においてdからzを決めるのに十分な情報が含まれていない。
大きくわけて二つの解法がある。
最小長さ解|dp-Gzp|=0としてzがzq方向に成分を全く持たない解を選択(下付きpは正の特異値の個数)。
ダンプ付き最小二乗解： eTe +ε2zTzを最小とする解を選択（ε2はあるスカラー定数；
データ同化の解は正規化した変数の空間でε2=1としていることに相当）

劣決定問題の解(1)：最小長さ解• 逆問題• 最小長さ解

dp-Gzp =0 を満たしつつ、zがzq方向に成分を全く持たない解

ここでG-Iは一般化逆行列であり，正の特異値p個が収められた対角行列を として

G-I=• この解は特異ベクトルを使ってまとめると
pz=∑
i =1uTi dλivizは(p空間を張る) vの方向の成分の和で構成されている。

劣決定問題の解(2)：ダンプ付き最小二乗解• 逆問題• ダンプ付き最小二乗解
Jを最小化する解→→

これはε≠0ならば必ず存在する。

• この解は特異ベクトルを使ってまとめると

λi2 uTi dz=∑ 2vi2λii =1 ε i + λi
pzは(p空間を張る) vの方向の成分の和で構成されている。

直観的イメージ(順問題)• n=m=2, p=1として、d=Gzを計算する場合は、
v2z2z(ヌル空間)

v1 (p空間) z1zからdを計算

u2(ヌル空間)d=d1(p空間)(u1方向の長さはλ1倍されている)

zは時刻0における正規化されたモデル変数の世界、
dは時刻tにおける正規化された観測変数の世界、Gは2つの世界をつなぐ演算子。

• 当然，u2方向にはdの成分は存在しない。

直観的イメージ(逆問題; 劣決定問題)• n=m=2, p=1として、d=Gz+eを計算する場合は、dからzを計算

v2(ヌル空間)z=z1v1 (p空間)du2(ヌル空間)d1(p空間)zは時刻0における正規化されたモデル変数の世界、dは時刻tにおける正規化された観測変数の世界、Gは2つの世界をつなぐ演算子。
• dからp空間を張るv1の成分の方向にzが得られる。

4次元データ同化の解の導出

変数の定義• 解析値ベクトル x a , 背景値ベクトル xb , 観測値ベクト
ル y ( x a と xb の長さは格子点数×変数の種類)• 背景誤差共分散行列B , 観測誤差共分散行列 R 。
• 観測行列 H : モデルの値を観測値と比較できるようにするための演算子(格子点内挿・変数変換など)。
• イノベーション y－Hxb ，インクリメント xa－xb

x1,b y x2,b1▲●▲x3,b y2 x4,b▲ ●▲x5,b▲

 x1,b  
 x2,b y1  0.5 x1,b + 0.5 x2,b  0.5 0.5 0 0 0 
xb = x3,b y =   H = Hxb = y x+x0.5
0.5
000.5
0.50
 2
3,b4,b x 4,b x 
 5,b 

摂動の時間発展• モデルの背景値が下記のように時刻0から時刻tまで時間発展することを考える。

xb ,t = M ( xb ,0 )• モデルの背景値の摂動 ε b の時間発展は

xb ,t + ε=
M ( xb ,0 + ε b ,0 )b ,tから引き算して、テイラー展開の1次の近似を適用し∂Mεb=
M ( xb ,0 + εb ,0 ) − M ( xb ,0 ) ≅εb ,0,t
∂xb ,0• M ≡ ∂M / ∂x0を時間発展演算子(接線形演算子)と呼ぶ。
タイムステップごとに分けるとM=Mt-1…M1M0

4次元データ同化の定式化(B0が正則)
時刻tの観測ベクトルytから時刻0のモデル変数x0を推定する

p ( y t | x0 ) p ( x0 )ベイズの定理の適用
p ( x0 | y t ) =＆観測誤差と背景誤差の無相関の仮定
p ( yt )以下の確率分布が正規分布に従うと仮定=p ( x0 )=p ( y t | x0 )

1

T
 1

exp  − ( x0 − xb ,0 ) B0−1 ( x0 − xb ,0 ) m
 2

( 2π ) B01
T
 1

−1
exp  − ( y t − Ht xt ) R t ( y t − Ht xt ) n 2( 2π ) Rt( | |は行列式を表すものとする）
T 11T
−1−1
p ( x0 | y t ) ∝ exp  −  ( x0 − xb ,0 ) B ( x0 − xb ,0 ) + ( y t − Ht xt ) R ( y t − Ht xt )  2
 2

次ページの代入・正規化を使うと
TT




J = z 0 z 0 + ( dt − Gz 0 ) ( d t − Gz 0 ) になる。

p ( x0 | y t ) の最大化問題

T11
T
J ( x0 ) = ( x0 − xb ,0 ) B0−1 ( x0 − xb ,0 ) + ( y t − Ht xt ) R t−1 ( y t − Ht xt ) の最小化問題
22∂J ( x0 )= B0−1 ( x0 − xb ,0 ) − MT HTt R t−1 ( y t − Ht xt )= 0∂x0

整理

(x a ,0 =xb ,0 + B + M H R Ht M(−10y t − Ht xb ,t =R t + Ht MB0M H
TTt

T−1tTt) M H R (y −H x )−1T−1tTtttb ,t)(R + H MB M H ) ( y − H x ) を代入
t0tTTt−1ttb ,t

Mを両辺の左から作用させるとカルマンフィルタ系列の解と一致する。
x a ,0 =xb ,0 + B 0 M HTTt

( R + H MB M H ) ( y − H x )tt0TTt−1ttb ,t

y t − H t x b ,t
イノベーションの定義 d=tz 0 x a ,0 − xb ,0
インクリメントの定義 ==z0

(B + M H R Ht M−10TTt−1t)−1

MT HTt R t−1dt−1/2
d=Rdtイノベーションの正規化 ttd Gz + e の逆問題=−1/2
ダンプ付き最小二乗解は インクリメントの正規化 z 0 = B 0 z 0T=
J ε 2 zT z + ( d − Gz ) ( d − Gz )−1/21/2
可観測行列の定義 G = R t H t MB 0の最小化で得られる−1=z(ε I + G G ) G d2T

T

なので、データ同化は逆問題と見なせる。

z =0(I+G GT)−1G T d t

可観測行列の特異値分解 G = UDVT

逆問題の解と4次元データ同化の解の対比•

4次元データ同化の解はε2=1としたダンプ付き最小二乗解と見なせる。

最小長さ解
逆問題

=d Gz + e

解の導出

dp-Gzp =0 かつzq 方向に成分なし

解
Gの特異値解析で分解ダンプ付き最小二乗解

4次元データ同化の解=d t Gz 0 + e
J=( d − Gz ) ( d − Gz )(TJ=
d t − Gz 0 d t − Gz 0+ ε 2 zT z の最小化+ z T0 z 0 の最小化2TTˆ −1UT d=z
εI+GGGd z ==z G=d Vp D()
0pp−1−Ipz=∑
i =1uTi dλivi)(T)( I + G G ) G dT−1Ttpλi2 uTi dλi2 uTi d tz=∑ 2v i z 0 = ∑vi
22λi
λii =1 1 + λii =1 ε i + λipd t :正規化されたイノベーション(=Rt-1/2(yt-Htxb,t)), z 0 :正規化されたインクリメント(=B0-1/2(xa,0-xb,0))
G:可観測行列(=Rt-1/2HtMB01/2), λi:Gのi番目の特異値, ui&vi: Gのi番目の左・右特異ベクトル線形データ同化の定式化(B0が正則でない場合に拡張)
モデル空間では、観測空間では R t−1/2 を掛けて、正規化を施す。

x 0 = B 0−1/2 x 0y t = R t−1/2 y t

p ( y t | x 0 ) p ( x 0 )
正規化された変数に対するベイズの定理
p ( x 0 | y t ) =＆観測誤差と背景誤差の無相関の仮定
p ( y )t以下の確率分布が正規分布に従うと仮定
T
 1





exp  − ( x0 − xb ,0 ) ( x0 − xb ,0 ) r
 2

( 2π )
T1
 1

−1/2
exp  − ( y t − R t H t xt ) ( y t − R t−1/2 H t xt ) n
 2

( 2π )=p ( x 0 )p ( y=t | x0 )1

最尤値を取る
可観測行列の定義の最小化

このままだと規格化された変数の世界のインクリ
x a ,0 − xb ,0 =B1/2
メントなので、0 z 0 でモデル変数の世界のインクリメントに直す

可観測行列Gとは何か？

• G=Rt-1/2HtMB01/2は正規化された時刻0のインクリメントを正規化された時刻tのイノベーションと対比するための行列。
 0 R t−1/2 H t MB1/2d t Gz=• いま考えているのは =0 z 0 という状況。
正規化された時刻0のインクリメント z 0B1/2
0 を左から作用
元の時刻0のインクリメント B 0 z 0 = z 01/2

M を左から作用元の時刻t のインクリメント MB1/20 z0 = zt
H t を左から作用
元の時刻t のインクリメントを観測空間に対応させた値 HMB1/20 z0

R t−1/2 を左から作用元の時刻t のインクリメントを観測空間−1/21/2
=GzR=HMBzd0t00tに対応させたのち正規化した値

特異値解析を用いた成分の分解• 逆問題の解を最小長さ解として求めると
pz 0 = ∑uTi d tviλi
ダンプ付き最小二乗解として求めると
pλi2 uTi d tz 0 = ∑vi
2λii =1 1 + λii =1• ここに出てくる変数がどう決まるか？
モデル変数と観測変数のマッピング⇒H背景誤差共分散と時間発展演算子⇒ui, vi, λi
観測値と背景値⇒dtT• 共通している部分 ui dt / λi v i はどういうもの？

()

( u d / λ ) v の解釈
Titii• G=UDVTは z 0を d t に結び付ける行列。
• ( uTi d t ) ui ：イノベーションのうちui方向の成分
TTu• ( ui d t / λi ) v i ： ( i dt ) ui を説明するインクリメント( u d ) u にTivj u d t 
 vi
 λi Tit～dti

対応するvi方向の成分 ujvi()

uTi d t ui

モデル変数空間

観測変数空間正規化された時刻0の正規化された時刻tの()u d t / λi v i の解釈（物理空間上で）
Ti• パターン1の d t ⇒ uTi d t > 0 ⇒ z0にB1/2方向の成分有0 viT• パターン2の dt ⇒ ui dt = 0 ⇒ z0に B1/20 v i 方向の成分無B vi
1/20MB v i
1/20d t (パターン1)

ui = R t−1/2 H t MB1/20 vi

d t (パターン2)

最小長さ解とダンプ付き最小二乗解の違い• 最小長さ解はpz 0 = ∑uTi d tλii =1vi

• ダンプ付き最小二乗解は
pz 0 = ∑i =1λ2i

1+ λ2iu d tTiλiviλi2• 両者は2 だけ異なる。
1 + λi

これを作用させることを
Tikhonovの正則化と呼ぶ• これは特異値λiが1に比べて十分小さいとき、対応する右特異ベクトルvi成分がダンプ付き最小二乗解にほとんど寄与しないことを意味する。

なぜλiが小さいときに解に反映させない？

• d tR t−1/2 ( y t − H t xb ,t ) にはO(1)のノイズがある。

vjuTi d t ui に観測ノイズが無い～
対応する
真のdtvi方向の成分 uj⇒観測ノイズを拡大してz0に反映するおそれ()vi( u d ) u に
Ti

vjti

実観測に～
基づいたdt

対応するvi方向の成分 uj≈viz 0と d tの関係に関する補足• データ同化により解を得たとする。
 のある２つ成分が同じ大きさのとき、• dtz 0 の空間では特異値が1に近いほど大きくなる。
 0のi番目の成分に• 一方、逆問題により得られた zついて0の特異値が1より十分に大きい場合には、その z
 のi番目の成分を良く説明す
成分を成長させると dt
る 0 の成分が特異値が1より十分に小さい場合には、zが極めて小さく生成され、時間発展で成長させて のi番目の成分の大きさにはならない。
も、 dt4次元データ同化の解の形からわかること

4次元データ同化の解の性質(1)解の構成に利用できる特異ベクトルの数は背景誤差分
散・観測行列・Mのランクの小さいものに制約される

G=R-1/2t

H t MB1/20• 解の構成に使われるのはせいぜい G = R t-1/2 H t MB1/20のランクまで（更にλiの小さい方向には修正効かず）。
• 観測1個:インクリメントとGの唯一のvの向きが一致• さらに B 0 がフルランクでない場合には、そのランクよりも小さくなる。特に、アンサンブル摂動を使って単純に B 0 = XXT とする場合、アンサンブルメンバー数にも制約を受ける。
(ただし、LETKFにおいては、別途考察する必要がある)4次元データ同化の解の性質(2)

時間発展演算子Mを通じて成長する摂動はダンプされにくく解の構成に使われやすい
G=R-1/2tH t MB1/20• 簡単のため、全てのモデル変数に対応する観測が得られたと仮定する。また、観測誤差共分散と背景誤差分散は単位行列の定数倍であるとする。

22=
Rt σ =I,Bσ=
Iob I, H0• この場合、 G = (σ b / σ o ) M となり、成長モード（λiが大きいMの特異ベクトル）はダンプされにくい。
• 全球大気モデルの例：傾圧不安定波(温帯低気圧)

4次元データ同化の解の性質(3)観測誤差共分散が小さい、背景誤差分散が大きい⇒より多くの特異ベクトルが解に反映される
G = R t-1/2 H t MB1/20• 例えば、R の要素が全て1/100になったとする。
⇒ G の定義より、特異値 λi はすべて10倍になる。
⇒ λi がO(1)より大きな特異ベクトル方向に修正がかかるので、さまざまな方向に修正がかかる• 同様に、B の要素が大きくなった時もさまざまな特異ベクトル方向への修正がかかりやすくなる。
• RとBの比が大事。

Lorenz63モデルを用いた4次元変分法• t=0.5の観測(1番目の要素のみ)からt=0の状態推定
真値 xtrue,0 = (-0.890, 0.496, 21.596)
B=I, R=10-4I背景値 xb,0 = (-0.940, 0.446, 21.096)
時間発展した背景値 xb,t = (-2.509, -4.982, 6.108)
観測値 yt = (-0.487) (1番目の要素のみを観測)イノベーション yt-Htxb,t = (-0.487-(-2.509)) = (2.021)
同化で得られた時刻0のインクリメント
正規化後 202.1z0=(0.0440, 0.0724, 0.0032)
可観測行列のλi=2.37×103, vi=(0.519,0.854,0.041)
z0とviがほぼ平行で値は202.1/2.37×103=0.085倍

カルマンフィルタ系列の手法との対比• ここで議論したインクリメントに時間発展演算子を作用させたztは時刻tにおけるカルマンフィルタ系列のインクリメントに一致。
• また、カルマンフィルタ系列のインクリメントについて、出力dtとモデルG=Rt-1/2HtMB01/2から入力ztを求める逆問題を解いていると考えてもよい。
(z0をzt, B0をBt, HtMをHtで置き換えて変形)⇒この場合でも、可観測行列の特異値解析でインクリメントを整理することができる。

可観測行列とその他の諸量の関係• 解析誤差共分散行列• 自己感度• DFS

データ同化による誤差共分散の縮小は特異値λiが大きい方向には有効だが、小さい方向には効いていない。
⇒EnKFによるアンサンブル摂動は成長モードの方向を潰した傾向にあることを示唆する

まとめ(伊藤・藤井，ながれ，2020)• 4次元データ同化の解は逆問題の解の一種とみなせる。
• その逆問題とは、可観測行列（モデル）と正規化された時刻tのイノベーション（出力）から、正規化された時刻0のインクリメント（入力）を求める問題である。
• 可観測行列は正規化・摂動の時間発展・モデル変数から観測変数へのマッピングを担う行列である。
• 可観測行列の特異値分解で、観測空間上の左特異ベクトルとイノベーションの内積に応じて、右特異ベクトル方向の成分がインクリメントに加わることがわかる。ただし、特異値が１より小さい場合にはノイズが大きいとみなし、寄与を小さくする。
• この解の形より、以下のことなどが分かる。
 インクリメントの構成に使われる特異ベクトルの数は、せいぜいmin(rank(H), rank(B), rank(R), rank(M))である。
 不安定モードに対応する構造は取り出されやすい。
 RとBそのものの大きさより、両者の比が重要である。