平成25平成25年度土木学会全国大会・研究討論会
25年度土木学会全国大会・研究討論会
2013/9/4

データ同化とは…
気象学・海洋学の分野で発達
物理数値シミュレーションモデルと観測データを統合する手法
シミュレーションのみでは適切に物理現象を再現できない（例えば）気象予測には適切な初期値の構成が必要
そもそもシミュレーションモデルには，モデル化されない力学，未知境界条件・パラメータといった不確かさが存在

観測データは物理的・社会的制約により得られる情報に限界がある

データ同化- データを活かしたシミュレーションと予測の方法

観測データを使い数値シミュレーション内の変数を修正
説明力の向上明治大学 総合数理学部現象数理学科/大学院先端数理科学研究科
数値シミュレーション

観測データ

中村和幸 (knaka@isc.meiji.ac.jp)

現実の情報の反映
2

2013/9/4

平成25年度土木学会全国大会・研究討論会

方程式からシステムモデルへ

データ同化の狙い・可能とするところ（日本周辺の簡易化した気象モデルの例を用いて説明）観測できない物理変数の推定
予測のための初期条件の構成∂x= cx 2 + L∂t

予報精度の向上を目指す
現業の天気予報ですでに行われている

平成25年度土木学会全国大会・研究討論会

湿度

ξ i = (Ti , H i ,U i ,Vi )

感度解析
効率のよい計測点・データの設計
シミュレーションモデルの改善
経験的パラメータの推定
境界条件の推定（シミュレーションモデルの選択）
3

各格子点は物理量 ξ i を持つ

温度 風速ベクトルxt = f t ( xt −1 )ii +1vt

ξ i +1xt = f t ( xt −1 , vt )2013/9/4

4xt = [ξ1 , ξ 2 ,..., ξ k ]Tk ：格子点数
平成25年度土木学会全国大会・研究討論会

2013/9/4

シミュレーションモデルとシステムモデル

観測情報と観測モデル

シミュレーションモデルの「誤差」，初期・境界条件などによる状態の誤差が反映されていない
このような誤差まで含めたモデルを考え，システムモ
デルを定式化

ほとんどの場合，観測情報はシミュレーションの情報に比べて圧倒的に不足．ダイナミクスを伴う逆問題．
さらに，時点間で独立な「観測ノイズ」もある
観測情報は，「その時点の全物理変数（＝全シミュレーション変数），および「観測ノイズ」が与えられれば，説明できる」という定式化

形式的にこのように書ける：

xt = f t ( xt −1 )シミュレーションモデル（離散時間・空間）

全シミュレーション変数
全観測変数
コーディング

観測ノイズ

yt = ht ( xt , wt )モデル化の誤差も含める：

xt = f t ( xt −1 , vt )dim( xt ) >> dim( yt )104~109101~106

全シミュレーション変数
モデル化誤差など
5

平成25年度土木学会全国大会・研究討論会

2013/9/4

6

2013/9/4

平成25年度土木学会全国大会・研究討論会

データ同化アルゴリズム一覧

状態空間モデルと推定アルゴリズム（非線形）状態空間モデル
全シミュレーション変数

Kalman filter
Extended Kalman filter
Ensemble Kalman filter (EnKF)

モデル化誤差など

 xt = f t ( xt −1 , vt ) yt = ht ( xt , wt )気象・海洋分野で主流
さまざまな変種あり（EAKF,ETKF,…）

dim( xt ) >> dim( yt )逐次型

Particle filter
ＳＩＲ filter の他，ＳＩＳ filter もある
理論的にいい面がある

全観測変数
観測ノイズ

4DVAR

変分（非逐次）型

気象・海洋分野におけるもう一つの主流（現業で使われている）

この形で表現できると，以下の二つが可能となる
観測 yt を得る毎の xt の推定（オンライン推定）
一定時間幅の全ての yt を用いて，その区間の xt を推定（バッチ）１時点の補間と隠れ変数の推定のみ

3DVAR
Nudging, OI, …

原始的8

平成25年度土木学会全国大会・研究討論会

どちらも，ベイジアンモデリングの観点で捉えられる
7

平成25年度土木学会全国大会・研究討論会

2013/9/4

2013/9/4

一期先予測（EnKF,PF(SIR,SIS)共通）

粒子近似：分布を実現値集合で近似する
状態1 N∑ δ ( xt − xt(|it )−1 )N i =1p ( xt | y1:t −1 ) ≅(i )t |tx

δ 関数を使った近似になっている∞∫ p( x | yt−∞

1:t −1

時刻{x }(i ) Nt |t −1 i =1∞

1∫ N ∑ δ (x − x9{x}xt(|t2−)1xt(|1t )−1xt(−N1|)t −1)dx = 1

2013/9/4

10xt(|tN−1)

t −1

t{x }

( 2)
t |t −1x尤度(1)t |txxt(|1t )−1xt(|t2 )時刻
条件の違うシミュレーションを2013/9/4
複数（Nパターン）繰り返す

平成25年度土木学会全国大会・研究討論会観測 :{x }(i ) Nt |t −1 i =1yt

( 2)
t |t −1x

各サンプルの尤度（データへのあてはまり）

xt(|tN )尤度

一期先予測からのサンプルフィルタ分布からのサンプル{x }(i ) Nt |t i =1

11xt(|1t )−1xt(|t1)

xt(|tN−1)

xt(|tN )フィルタリングt

時刻
尤度に比例して復元抽出
2013/9/4yt p( yt | xt(|it )−1 ) ( j)  ∑ p ( yt | xt|t −1 )  j

フィルタリング

平成25年度土木学会全国大会・研究討論会

観測 :xt(|t2 )各サンプルの尤度（データへのあてはまり） p( yt | xt(|it )−1 ) ( j)  ∑ p ( yt | xt|t −1 )  j
xt(|tN−1)

t

一期先予測

xt

一期先予測からのサンプルフィルタ分布からのサンプル{x }フィルタ分布からのサンプルフィルタリング(PF(SIS))
state(i ) N
t |t i =1(i) Nt |t −1 i =1

一期先予測からのサンプル

N : 実現値の個数

xt(i ) Nt |t −1 i =1{x }xt(−11) |t −1xt(−21)|t −1

フィルタリング(PF(SIR))
statext(|it )−1

使用した観測の最後の時刻

N(i)tt |t −1ti =1−平成25年度土木学会全国大会・研究討論会∞)dxt = 1f t ( xt(−i )1|t −1 , vt(i ) )N(i )t −1|t −1 i =1
シミュレーションi 番目の実現値

1 N
p ( xt | y1:t ) ≅∑ δ ( xt − xt(|it ) )N i =1p ( xt | y1:t −1 )

xt

12

時刻
各サンプルの重みを積で蓄積していく
平成25年度土木学会全国大会・研究討論会

2013/9/4

研究紹介１（続）：透水パラメータの逐次推定

研究紹介１ 神戸空港・地盤沈下の予測データ同化

地盤変形モデル（弾塑性変形）

沈下量データ・直接見ることができない地中の土の状態がわかる
・予測精度の向上で，中途での意思決定が可能に（珠玖隆行（岡大），村上章，藤澤和謙（京大），西村伸一（岡大）各氏との共同研究）

13

2013/9/4

平成25年度土木学会全国大会・研究討論会

研究紹介２ マウスの体内時計モデルのデータ解析

14

平成25年度土木学会全国大会・研究討論会

2013/9/4

研究紹介２（続）：パラメータの統計的推定と並列・ハイ
パフォーマンスコンピューティングの活用
１億サンプルによる近似 10万サンプルによる近似

データ同化

現実の系を表すには不完全
未知パラメータ
ノイズ，欠測など

マウスの体内時計（概日周期）のメカニズムの解明と予測（長崎正朗（東北大），宮野悟（東大），吉田亮，樋口知之（統数研）各氏との共同研究）
15

平成25年度土木学会全国大会・研究討論会

2013/9/4

並列計算を活用することで，分布をより正確に近似・推定ができる
16

予測精度が上がるだけでなく，興味ある現象が起こる確率を正確に出すことができる

平成25年度土木学会全国大会・研究討論会

2013/9/4

なぜ分布なのか？

金融工学における商品価値評価・意思決定との関連確率分布で捉えることで．．．

金融市場は不確実性を持つため，確率的な評価がなされてきたどの範囲にどの確率で入るか考えられるので，リスク評価に用いることができる
さらに，「入ってほしくない」範囲の各階級に損失を定義することにより，トータルコストの計算が可能となる

各種確率解析による金融商品の将来価格の評価
たとえば，ブラックショールズモデルでは「株価は幾何ブラウン運動に従う」と仮定
将来の利益または損失の期待値で現在価値を評価
すなわち，（将来の利益または損失）×（その確率）の和

そもそも，大規模空間では，不確実性は必ず現れる
各パラメータが10通りの値を取るとして，10個あったら1010 (＝100億)通り全ての場合をデータから集めるのは不可能

データ同化により将来状態の「確率」を予測すること
で，期待価値によって「現在の価値」を評価可能

不確実性のある現象を定量的に捉え，意思決定するための方法
17

2013/9/4

平成25年度土木学会全国大会・研究討論会

ベイズ推定と確率的モデリング（１）現象Ｘが発生した条件下でデータＹが得られる確率18

2013/9/4

平成25年度土木学会全国大会・研究討論会

ベイズ推定と確率的モデリング（２）現象Ｘが発生する「もともとの」確率
ベイズの定理p (Y | X ) p ( X )p(Y )データＹが得られた時に∑p( X | Y ) =p (Y ) =生成現象

データ p( X | Y ) =p (Y | X ) p ( X )p (Y )p (Y | S ) p( S )( 必要なのは p(Y | X ) と pより，( X ) ．)
S ∈Ω

データＹの生成確率

現象がＸである確率生成モデル
ベイズの定理

データ生成モデルと現象の発生確率を与えれば，データから現象の説明ができるようになる（データ同化も同じ仕組）

生成現象

データ

：事前知識やシミュレーションp(Y | X ) と事前分布 p( X ) の適切な構成（＝モデリング）を通じて，有用な事後分布

p ( X | Y ) による情報抽出をめざす．

：観測を表す式ベイズ補足

19

平成25年度土木学会全国大会・研究討論会

2013/9/4

20

平成25年度土木学会全国大会・研究討論会

2013/9/4

明治大学における取り組み

大量データ・大規模計算の時代の解析人材育成
数理科学・計算機科学を知識と研究の基盤とし，応用分野に挑戦していく人材はまだ少ない「社会に貢献する数理科学」の観点から，次の三本柱に関する研究・教育を進めている現象から数式を作る「モデリング」
数式を計算機上で解いて再現する「シミュレーション」
上記結果や数式の解析から現象を理解する「数理解析」

数理科学・情報科学側の問題
応用研究の位置づけ
応用分野の広さによる教育の難しさ
研究機会創出の困難性問題を「はっきりさせる」までに時間がかかる

研究：先端数理科学インスティテュート（研究機関）・
GCOE「現象数理学の形成と発展」
各分野研究者との交流
テーマ別の研究集会・大学院集中講義

地盤変形の共同研究には，方向性をはっきりさせるのに半年かかっている

集中的にコミットできる環境が必要だが．．．
プロジェクトベースのポスドク雇用とキャリアパスの関係
21

平成25年度土木学会全国大会・研究討論会

2013/9/4

まとめ
データ同化は大規模シミュレーションとデータを統合
基礎方程式がしっかりしている場合に良く使われる非線形性により，難しさが変わってくる予測のための大規模シミュレーション・大規模データ
解析には，新しい道具が必要
確率的モデリングが必須
ベイズ推論により，不確実なものを不確実なまま捉えてデータと融合できる期待価値を計算することで，意思決定まで考慮にいれた議論が可能となる

23

平成25年度土木学会全国大会・研究討論会

2013/9/4

教育：大学院先端数理科学研究科（2011年度開設）・総合数理学部現象数理学科（2013年度開設）
数理と情報を軸に上記を学ぶカリキュラム
各分野と数理科学分野の相互交流可能な人材を育成
22

平成25年度土木学会全国大会・研究討論会

2013/9/4