逆問題としてのデータ同化

琉球大学伊藤耕介
(2014/5/29)
逆問題としてのデータ同化
本資料について
本資料は，データ同化の結果として得られる解 (一般的に，解析値と呼ばれる1) ) がどのような性質
をもつものであるかについて，逆問題の観点からレビューを試みたものである．データ同化の基礎知
識があることを前提にしている．まだまだ完成度は十分でないと考えられるので，著作としての引用
はご遠慮いただきたい．もし，誤りを見つけられた場合には [email protected] までご連絡いただけ
ると幸いである．
本解説では，簡単のため，全体を通じて以下のような条件を設けている．
• 逆問題のうち簡単な線形離散逆問題2) についてだけ考える．
• 複素数は扱わない．
• 内積は要素同士の積の単純和である自然内積として定義する．
• 背景誤差共分散行列と観測誤差共分散行列は正しく設定されている．
(2014.05.29 追記) 式に誤りがあったので修正しました。
逆問題とは？
1
さて，ある理論を表現するモデル3) G が，入力変数 x と出力変数 yout を以下のように結びつけて
いるとする．
yout = Gx
1)
2)
(1)
正確には，場を表現するものであるから解析値場と呼ぶべきであろう．
状態変数の時間発展について非線形性が残っていてもよいが，誤差の場の時間発展については線形性が仮定できると
する．
3)
ここでいうモデルとは，数値計算に用いる力学モデルも含むが，それだけでなく，より一般的に入力変数と出力変数
を結びつけるものを指している．
1
ここで，出力変数 yout と比べることのできる観測値 y が得られたときに，入力変数 x を推定しよう
とするのが逆問題（ここでは特に線形離散逆問題）である．気象学・海洋学における一般的なデータ
同化は，観測値が y に相当し，格子点上の状態変数が x に相当する逆問題であると言ってもよい．
G が正則で逆行列が陽に求められ，yout と y が完全に対応するならば，両辺に G−1 を作用させる
ことで，x を求めることができる．このような問題を平衡決定問題 (ever-determined problem) と呼
ぶ．その一方，観測の数が状態変数の自由度に比べて少ない，もしくは，観測の数が状態変数の自由
度と同じか多いけれども状態変数すべてを復元できるだけの情報が含まれていない場合を指して劣
決定問題 (under-determined problem) であるという．観測の数が少なければ状態を全て決定できない
のはもちろんであるが，観測の数が状態変数の自由度より多くても，劣決定問題となりうることに注
意が必要である．例えば，観測がモデル領域の一部で極めて大量に得られただけであるとき，観測か
ら離れた地点上の変数の影響が集中観測点に及ばないならば，それらの変数は適切に復元できないの
である（このとき G では正則ではない）．劣決定問題とは反対に，観測の数が状態変数の自由度に比
べて多く，全ての状態変数に関する情報が含まれているとき，これを，優決定問題 (over-determined
problem) という．
一般的な大気・海洋におけるデータ同化では，劣決定問題を解いていることがほとんどである（ま
れに，格子点が粗い海洋長期再解析の初期値を逆問題で推定する場合には，実際に優決定問題とな
ることもあり得る）．平衡決定問題でない限り，逆行列を作用させて逆問題を解くことはできないの
で，別の方法で入力変数を求めることになる．優決定問題の場合には，最小二乗法を用いて解を求め
ることができる．
最小二乗法とは，逆問題として，以下のような残差ベクトル e を導入し，
y = Gx + e
(2)
残差ベクトルの最小化問題の解を x とするという方法である．
劣決定問題の場合には，残差ベクトル e をゼロとするような解 x が無数に存在することになるた
め，先験情報と呼ばれる情報を数学的な条件として追加し，解を探索することが一般的である．これ
について，次節で詳しく述べる．
2
一般化逆行列を用いた最適解
観測の数と推定すべき状態変数の数をそれぞれ r, p とし，観測の数が推定すべき状態変数に比べて
少ない劣決定問題 (r ≪ p) を考える．また，変数を観測値空間に内挿する観測行列を H とする．こ
こでは，H が逆問題における「モデル」となる．本節及び次節では数値モデルの力学を通じた時間
発展を考えないが，これを自然な形で拡張することもできる (6 節)．
2
逆問題として，観測値ベクトル y から状態変数ベクトル x を推定する問題は，以下のように表さ
れる．
y = Hx + e
(3)
あとで行う比較のために，x′ ≡ B− 2 x, y′ ≡ R− 2 y という変数変換をかけておき4) ，第一推定値を x′b
1
1
とすると，
y′ − R−1/2 HB1/2 x′b = R−1/2 HB1/2 (x′ − x′b ) + R−1/2 e
(4)
と書き換えられる．
劣決定問題において解を定める方法としては，式 (4) において，e = 0 を満たしつつ，x′ − x′b が最
小となるような「最小長さ解」を求めることが一般的である．単に e = 0 とするだけでは，式の数が
未知変数の数より多く，解が無数に存在するので，e = 0 を満たすもののうち，x′ − x′b が小さい（解
析値が背景値に近い）ものを解として選択するわけである．ここでは，x′ − x′b が小さいものがよい
ということが「先験情報」として利用されている．
この問題の解は，Moore-Penrose 型一般化逆行列 (·)−I (以下，単に一般化逆行列と呼ぶ5) ) を用いて，
x′a = x′b + (R−1/2 HB1/2 )−I R−1/2 d
(5)
と表される．ただし，イノベーション (y − Hxb ) を d とおいている．R−1/2 HB1/2 のことは可観測行
列 (observability matrix) と呼び本資料では L で表すことにする．なぜこのような名前がついている
のかについては，4 節で説明するが，これを使うと，式は以下のように表せる．
x′a = x′b + L−I R−1/2 d
(6)
A が UDVT と特異値分解されるとき，その一般化逆行列 A−I は VD−I UT となる．D は対角成分以
外はゼロ行列で，対角成分は非負の特異値が収められている．そこで，可観測行列 L を UDVT と特
異値分解したとき，その一般化逆行列は，
(
)−I
L−I = R−1/2 HB1/2
= VD−I UT
(7)
として表される6) ．D−I は対角成分以外はゼロで，対角成分は非負の特異値の逆数からなる行列で，
特異値分解ができさえすればすぐに計算できる．さらに，これを各特異ベクトルの成分ごとの重ね合
4)
B は，x′ x′T の期待値ではなく，xxT の期待値である．R についても yyT の期待値としておく．もし，B や R が正
しく設定されているならば，x′ や y′ は O(1) となっている．
5)
A の Moore-Penrose 型一般化逆行列 A−I とは，A の随伴行列を A∗ と表すことにすると，AA−I A = A，
A−I AA−I = A−I ，(AA−I )∗ = AA−I ，(A−I A)∗ = A−I A を満たすものであり，ただ一つに定まる．A が正則の
ときには，逆行列 A−1 がこの条件を満たしていることから，逆行列の一般化となっていることが分かる．また，一般化
逆行列と呼ばれるものには何種類かあるが，Moore-Penrose 型一般化逆行列が最もよく使われる (?)．
6)
複素数体について考えるときや内積の定義が自然内積でない場合には，L = UDV∗ として，特異値分解に転置行列
でなく随伴行列 V∗ を導入するべきだが，ここでは実数体を考え，内積の定義は自然内積としているので転置行列 VT は
随伴行列に等しくなっている．
3
わせの形に書き直すと，以下のように変形できる．
x′a
=
x′b
+
m
∑
uT R−1/2 d
l
l=1
λl
ただし，
L = R−1/2 HB1/2 = UDVT =
vl
m
∑
(8)
ul λl vTl
(9)
l=1
ここで，U は r 行 r 列の直交行列かつ左特異ベクトル ul を並べたものに等しく，V は p 行 p 列の直
交行列かつ右特異ベクトル vl を並べたものに等しい7) ．ここでは，D の値の大きい順にそれぞれ対
応する特異ベクトルを並べ替えたものとする．また，特異ベクトルの性質として，vl は LT L の固有
ベクトルに等しく，ul は LLT の固有ベクトルに等しい．
数学的な説明が長くなったが，要は，劣決定問題として観測ベクトル y をもとに状態変数 x を推
定するという問題は，式 (3) において e = 0，かつ，背景値からの差が最小となるものを解 (式 (6) も
しくは式 (8)) として選択するというわけである．一見して，合理的な解を選んでいるようにも思え
るが，実は，一般的なデータ同化では，これとは異なる解を選択していることになる．それについて
は次節で述べる．
線形不偏分散推定
3
一般的なデータ同化手法では，先ほどと同じように，観測値ベクトル y の情報を用いて，状態ベ
クトル x を推定するのであるが，
xa = xb + K (y − Hxb )
(10)
という式を考え，解析誤差共分散 Pa の trace を最小化する K を求めている．この最小分散推定問題
は以下のように定式化される．
(
)
Pa ≡ E (xa − xt )(xa − xt )T
(11)
∂
trace(Pa ) = 0
∂K
(12)
ただし，E(·) は期待値を表すものとし，下付きの t, a, b はそれぞれ真値，解析値，第一推定値である
ことを指す．この問題の解は，
(
)−1
xa = xb + BHT R + HBHT
(y − Hxb )
(13)
である．これが一般的な線形のデータ同化手法によって採用される解である．
この解は最小分散推定の解であると同時に，以下の評価関数を最小化する解でもある．
1
1
J = (x − xb )T B−1 (x − xb ) + (y − Hx)T R−1 (y − Hx)
2
2
7)
特異値分解にはこのほかにも，r 行 l 列の U ，l 行 l 列の D，p 行 l 列の V などに分解する流儀などもある．
4
(14)
また，前節と同様の x′ ≡ B− 2 x,y′ ≡ R− 2 y という変数変換を導入しておく．これにより，評価関
1
1
数は，
1
1
J = |x′ − x′b |2 + |R−1/2 HB1/2 (x′ − x′b ) − R−1/2 d|2
2
2
のように書き直すことができ，解は，
(
)−1 (
(
)T
)T
−1/2
1/2
−1/2
1/2
′
′
xa = xb + Ip + R
HB
R
HB
R−1/2 HB1/2 R−1/2 d
(15)
(16)
と変形できる．ただし，Ip は p × p の単位行列である．これを特異ベクトルの性質 (UUT = Ir や
VVT = Ip など) と可観測行列 L = R−1/2 HB1/2 の特異値分解を使って変形すると
(
)−1
x′a = x′b + Ip + VDT DVT
VDT UT R−1/2 d
m
∑
λ2l uTl R−1/2 d
vl
= x′b +
2
1
+
λ
λ
l
l
l=1
(17)
(18)
となる．この式の解釈は 4 節で行う．
一般化逆行列を用いた最適解の何がダメか？
4
2 節のように一般化逆行列を作用させて解を求めることは，y′ = Hx′ を満たすもののうち，x′ − x′b
が最小となるという解 (式 (8)) を選んだことになる．これは，一見合理的な解の選び方であるように
思うが，データ同化で通常使われる解 (式 (18)) とは，
ϕl =
λ2l
1 + λ2l
(19)
の分だけ異なっている．すなわち，線形不偏分散推定の解は，一般化逆行列を作用させた解に対し，
特異値が λl ≪ 1 となる成分をフィルターをかけて取り除いたものということになる．それでは，式
(8) ではなく，式 (18) がデータ同化の解としてよく使われるのはなぜであろうか？8)
これは，逆問題の観点でいうと，
「Picard の条件からくる要請に対して，Tikhonov の正則化を施し
た」ことになる．Picard の条件 (Discrete Picard’s condition) について説明しよう．式 (18) と式 (5)
の両方に含まれる
uTl R−1/2 d
vl
λl
(20)
の項は，特異値 λl が小さいものほど値が大きなベクトルとなる．なぜなら，u や v は単位ベクトル
であり長さは変わらない．そして，観測値には通常ノイズが含まれるので，d はどのようなモードで
あっても観測ノイズの大きさ程度かそれより大きいのが典型的となる．そのため，この項は，特異値
の非常に小さなモードに存在する観測ノイズばかりを拡大するのである．ひるがえってみると，式
8)
1990 年代までのデータ同化に関する論文などでは，式 (5) も解として用いられることが多かった．
5
(8) はこの項の重ね合わせからなっており，特異値の非常に小さなモードに存在する観測ノイズばか
りを反映してしまう可能性が高い．Picard の条件とは，uTl R−1/2 d の典型的なサイズと λl の典型的
なサイズを比べて，前者の方が小さいモード l であれば信頼できないとして，その向きの修正はかけ
ないようにすべきという考えである．
結局，典型的なデータ同化の解 (式 (18)) は，特異値の小さなモードについて観測値から得られる
情報を
λ2l 1
1
→
λl
1 + λ2l λl
(21)
と置き換えてダンプすることによって，安定化を施している．この操作を Tikhonov の正則化と呼び，
λ2l /(1 + λ2l ) を Tikhonov のフィルター係数などと呼ぶ．変分法データ同化手法による解が式 (16) に
等しいことを思い出すと，
1
1
J = (x − xb )T B−1 (x − xb ) + (y − Hx)T R−1 (y − Hx)
2
2
(22)
データ同化で通常使われる解は，実は，一般化逆行列を使って得た解に Tikhonov の正則化を導入し
た結果だったのである．
5
可観測行列の解釈
式 (18) の方がデータ同化の解として選択される理由について述べたところで，いくつか追記して
おこう．結局，L = R−1/2 HB1/2 になぜ「可観測行列」という大仰な名前がついているかと言えば，
式 (18) に見えるように，この行列の特異値分解の結果が，イノベーションにどのような解析インク
リメントを生成するかに深くかかわっているからだ．この行列とデータ同化との関わりについて，基
本的な結果をまとめておくと，
• 可観測行列を特異値分解したとき，1 より十分に大きい特異値に対応するモードはデータ同化
の解に反映されるのであり，1 より十分に小さい特異値に対応するモードはダンプされ反映さ
れない．
• 特異値の非ゼロ成分は最大で m 個，つまり，最大でも観測値の数 r かモデルの自由度 p のうち
小さい方の個数しかとれない．
• 観測誤差が小さい場合には R−1 は大きくなるため，特異値も大きくなる．これにより，解析イ
ンクリメントに関わるモード (最大数 m) のうち，ダンプされないモードの数は多くなる．
6
6
４次元データ同化における解の性質
入力変数 x が以下のような時間発展に従うことを考える．
xt = M (x0 )
(23)
ここでは下付きの文字は時刻を表しており，M は時刻 0 から t までの時間発展を表す非線形演算子
である．そして，観測値と入力変数との間に以下のような関係があることを考える．
yt = Hxt + e
(24)
ただし，観測演算子は線形を仮定している．一つの時刻ステップ t にだけ観測値が存在するものとし
たときの同化について考える．線形なシステムを考えているので，別のタイムステップに観測値があ
る場合には，解は単純な重ね合わせとなる．
δxt = M (x0 + δx0 ) − M (x0 )
∂M
∼
δx0
∂x
(25)
(26)
∂M/∂x は接線形演算子と呼ばれる行列であり，以下，M と記す．
線形離散の逆問題を考えるとき，４次元変分法・カルマンフィルタの解はともに，
(
)−1
(y − Hxb,t )
xa,t = xb,t + MBMT HT R + HMBMT HT
(27)
となる．また，3 節で述べてきた最小分散推定による解に関する議論は，以下のような変数の置き換
えで，そのまま 4 次元変分法とカルマンフィルタの解に関する議論に拡張することができる．
B → MBMT
(
)−1
−1/2
1/2
′
′
x ≡B
x → x ≡ MB
x
(28)
L ≡ R−1/2 HB1/2 → L ≡ R−1/2 HMB1/2
(30)
(29)
これにより，データ同化による場の更新は式 (18) と同じ式を用いることができる．ひとつ注目した
いのは，解析インクリメントの方向を決める特異ベクトル vl は LT L = R−1/2 HMBM1/2 HT R−T /2
の固有ベクトルに等しいという点である．もし，ある時刻において，制御変数に対応する観測値が全
て得られるとし，観測誤差共分散行列が R = σo2 I と表されるならば，vl は MBMT の固有ベクトル
の向きに等しくなる．これは，背景誤差共分散を時間発展させたものであり，時間発展 (M を作用さ
せること) によって成長するモードがデータ同化によって修正の対象となりやすいことを如実に示し
ている．
7
謝辞
本文章の誤りを指摘していただいた東京大学の高野雄紀様に感謝申し上げます．
参考文献
メンケ著，柳谷俊・塚田和彦訳 1997: 離散インバース理論―逆問題とデータ解析. 古今出版.
Johnson et al. 2005: A singular vector perspective of 4D-Var: Filtering and interpolation, QJRMS,
131, 1-19.
8

Download Report