流体工学分野のデータ同化に関する最近の取組み
三坂 孝志（東北大学 学際科学フロンティア研究所）
大林 茂（東北大学 流体科学研究所）

第1回設計情報駆動研究会
2016年3月5日
1目次流体工学分野におけるデータ同化
 計算コストの克服と適切な問題設定

 近似モデルによるデータ同化の高速化
 RANS乱流モデルのパラメータ推定
 次元縮約モデルによる高速流動予測（今回は省略）

 データ同化の問題設定（DAシステムの評価）
 可観測性と可制御性（主に可観測性に関して）2

流体工学分野におけるデータ同化

3

データ同化とは? ・・・ 気象予測に不可欠なツール
 気象場の再現性向上のために，観測値をモデル予測に融合し，統計的に尤もらしい初期値を生成する
 気象予測のための４つ目の道具として利用されている
• 理論
• 観測
• 数値シミュレーション• データ同化 気象再解析データセットの作成
気象庁ウェブサイトより(http://www.jma.go.jp/jma/kishou/know/whitep/1-3-6.html)

データ同化を利用した
数値気象予測雲画像（観測）4データ同化とは? ・・・ 統計的データ解析及びデータ同化との違い
統計を利用したデータ解析
𝑥𝑥 少数のパラメータ（a, b）を持つモデルを計測データに当てはめる
線型モデルを最小二乗法で当てはめる場合
𝑁𝑁

𝐽𝐽 𝑎𝑎, 𝑏𝑏 = � 𝑦𝑦𝑛𝑛 − 𝑎𝑎𝑡𝑡𝑛𝑛 − 𝑏𝑏 2 の最小化 ・・・
𝑛𝑛=1

線形モデルの例
計測値𝑦𝑦𝑛𝑛
線型モデル
𝑥𝑥 = 𝑎𝑎𝑎𝑎 + 𝑏𝑏

𝜕𝜕𝐽𝐽
𝜕𝜕𝐽𝐽
= = 0から𝜕𝜕𝑎𝑎𝜕𝜕𝑏𝑏

データ同化

cf. 逆問題，スパースモデリング

状態𝑥𝑥
初期値の修正 データに大自由度のシミュレーションモデル（ナビエ・ストークス方程式など）を当てはめることで実現象を推定
 このとき変数となるのはモデルの初期・境界条件・パラメータ
 疎な観測をヒントに観測のない領域の状態（または観測されていない量）を物理モデルに基づき推定𝑡𝑡

変分型初期値推定の例
モデル予測

計測値𝑦𝑦𝑛𝑛

時間𝑡𝑡5

データ同化による設計・開発技術の革新（目指す姿）
CAE解析

CAE解析

模型試験

センサーネットワークを擁した先進計測

データ同化を生かした解析プロセス
性能評価

性能評価

得られた結果の統合（人為的プロセス）
再設計/
最適設計

データ同化の導入

開発プロセスへのCAE導入により試験コストの低減を実現してきた統一的視点での性能評価

ビッグデータ
統合解析

CAEモデル高度化・最適化

試験システ
ムの高度化

再設計/
最適設計

データ同化による
CAE解析と模型試験の統合的アプローチ

データ同化による試験・CAEの信頼性・精度向上
 時間コストの低減，製品の高度化による競争力強化
http://www.ifs.tohoku.ac.jp/edge/DAE/

CAE: Computer-Aided Engineering6流体工学分野におけるデータ同化
気象分野のデータ同化
 設定の難しい初期・境界条件を観測データに基づき決めるというデータ同化の目的が明確
 予測結果を（信頼性の高い）事前分布として用いることができる
 観測データ・数値モデル・問題設定は（ほぼ）固定で運用される

流体工学分野のデータ同化
 解析精度の向上という目的はあるものの問題設定は多種多様…
データ同化への要求
 数値解析単体よりも良い（実験・実現象に近い）結果が得られる
 もしくは，数値解析と実験・実現象との相関を定量化できる
 計算コストが増加しない（既存数値解析に要する時間と同程度）

7流体工学分野におけるデータ同化の方向性
データ同化の方向性① 次元縮約モデル（ROM）を用いた計算コスト削減（シミュレーションモデルの自由度削減は状態推定を容易にする）
 流体解析モデルは数百万～数億自由度（格子点数×変数）② モデルパラメータの最適化（乱流モデルなど）
 モデルはほぼ正しい．工学問題の初期・境界条件はほぼ決まり（モデル構築自体が困難な問題もあるが，モデルパラメータ推定に落とし込むことができればデータ同化を適用可能）③ 計測方法の工夫（流れ場に応じた計測，効率的に大自由度を制御） 通常は与えられている．計測の事前検討をデータ同化で行う

8

データ同化のコスト削減を実現する方法
Kriging

高次元（非線形）数値モデル(Full-Order Model, FOM)

定常問題
応答曲面法
(Response surface)

RBF

モデルパラメータ推定
境界条件推定（少数）ガレルキン射影，RBF（POD基底，DMD基底）

4DVar, EnKF

KF, EnKF

背景場による正則化

PF, MCMC

インクリメンタル法

非定常問題

非定常場の推定次元縮約モデル(Reduced-Order Model, ROM)

応答曲面法
 Kriging法 放射基底関数(RBF)法
 多項式
 など

KF, EnKF
PF, 4DVar

場のモード分解 固有直交分解(POD) 動的モード分解(DMD) など

効果的なデータ同化

データ同化・統計手法
 アンサンブルカルマンフィルタ(EnKF)
 マルコフ連鎖モンテカルロ法(MCMC)
 粒子フィルタ(PF) 4次元変分法(4DVar)

9データ同化の問題設定の検討方法
 感度解析や制御理論の枠組みを利用

計測手法・計測
データの検討

可観測性
計測感度

初期問題設定状態推定する
変数の検討

問題設定が良かった場合

可制御性
感度解析

経験的可観測性グラム行列

経験的可制御性グラム行列

アジョイント感度解析

統計手法（モンテカルロ法,
Stochastic Collocation法）

アンサンブル感度解析

効果的なデータ同化

10近似モデルによるデータ同化の高速化
RANS乱流モデルのパラメータ推定（株式会社IHIとの共同研究）

11

タービン翼のフィルム冷却

本ページの図はWebから適宜引用させて頂いております 環境適合性の高いジェットエンジン
‐‐‐‐

騒音の低減
NOx排出量の削減
省エネルギー
CO2排出量の削減

熱効率の向上  タービン入口温度の上昇
 部材への熱負荷の増大‐ 熱負荷に耐え得る材料の開発‐ 冷却技術の開発

タービン入口温度の変遷

実用化された冷却技術
‐‐‐‐

対流冷却
インピンジメント冷却
フィルム冷却全面フィルム冷却

冷却技術の進歩12冷却効率の数値シミュレーションによる予測
フィルム冷却タービン翼のRANSシミュレーション
冷却流の流路G. M. Laskowski, et al., HEATTRANSFER PREDICTIONS OF
FILM COOLED STATIONARY
TURBINE AIRFOILS, GT200727497, Proceedings of ASME
Turbo Expo 2007, 2007.

フィルム冷却孔冷却効率の予測性能はRANS乱流モデルに大きく依存する
K. L. Harrison, et al., Comparison of RANS Turbulence Modelsfor Prediction of Film Cooling Performance, GT2008-51423,
Proceedings of ASME Turbo Expo 2008, 2008.

目的：データ同化によってフィルム冷却効率の予測精度を向上する
13

データ駆動型の乱流モデリング
 データ同化 はく離流れでパラメータを推定し，最適化したパラメータを用いて衝撃波を伴う流れのはく離を予測
H. Kato, et al., Optimization of Parameter Values in
the Turbulence Model Aided by Data Assimilation,
AIAA Journal, Vol. 54, No. 5, 2016.

 統計手法 サンプル計算で事前分布を作成し，応答曲面法を用いてベイズ推定
J. Ray, et al., Bayesian Parameter Estimation of a k-ϵ
Model for Accurate Jet-in-Crossflow Simulations,
AIAA Journal, Vol. 54, No. 8, pp. 2432-2448, 2016.

 機械学習 基準流れでの学習結果を用いて，異なる流れのRANS予測の不確性を評価
J. Ling, et al., Evaluation of Machine Learning
Algorithms for Prediction of Regions of High Reynolds
Averaged Navier Stokes Uncertainty, Physics of Fluids,
Vol. 27, 085103, 2015.

近年，乱流モデルのパラメータを修正することへの抵抗感が薄らいできた?

14

データ駆動型の乱流モデリング手法の概要
解析の流れ

誤差評価関数
計測方法の評価

乱流モデルのパラメータ選定（上下限値の調査）
実験計画法により乱流モデルパラメータを変化させて計算計測誤差分布を得る

例）SST乱流モデル，𝜔𝜔方程式の崩壊項の係数（内層𝛽𝛽1 , 外層𝛽𝛽2 ）計算-計測誤差のKriging応答曲面

誤差に関する応答曲面の生成
マルコフ連鎖モンテカルロ法，粒子フィルタによるパラメータおよび確率密度推定
分散分析・固有直交分解
パラメータ推定モデルパラメータと温度場の関連づけ（POD）

分散分析によるモデルパラメータの寄与度

粒子フィルタによるモデルパラメータの推定
𝛽𝛽1 = 0.0389

b1b2𝛽𝛽2 = 0.0424

15

Kriging応答曲面法
 サンプル点を用いてモデル出力の近似関数を生成 パラメータ推定をKriging応答
曲面の出力で行うことで計算コストを大幅に低減yˆ (x) = µ + ε (x)局所モデル

大域的定数モデル

大域的定数モデル
 定数（サンプル点の平均）
1T R −1 fµ (x) = T −1
1 R 1関数値

サンプル点真の関数値
Kriging モデル
標準誤差

局所モデル
 ガウス関数の重ね合わせ
Additional
sample point

ε (x) = r T R −1 ( f − μ)Corr [ε (x ), ε (x )] = ∏ exp − θ x − x mijk

設計変数[ikj 2kk

Rij = Corr ε (x i ), ε (x j )] r = Corr [ε (x), ε (x )]ji16マルコフ連鎖モンテカルロ法
 モンテカルロ法(MC)：ランダムに値を発生
 マルコフ連鎖モンテカルロ法(MCMC)：相関を持つ乱数を順次発生
Metropolis-Hastings法
①②
③
④
⑤12

𝐿𝐿: 尤度（誤差）初期位置（パラメータ）を設定: 𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛
𝑃𝑃𝑃𝑃: 事前分布
新しい位置の計算: 𝒙𝒙𝑛𝑛𝑒𝑒𝑒𝑒 = 𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 + 𝝈𝝈 rand()
遷移確率の計算: 𝛼𝛼 = min 1, 𝐿𝐿 𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 𝑃𝑃𝑃𝑃 / 𝐿𝐿 𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 𝑃𝑃𝑃𝑃
𝛼𝛼 < rand()なら𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 = 𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 ，そうでなければ𝒙𝒙𝑛𝑛𝑛𝑛𝑛𝑛 を更新しない②に戻って繰り返す

10推定値（真値）

86Y420-2-4

初期位置-6-6-4-202X4681017粒子フィルター 非ガウス・非線形フィルター① 初期アンサンブルを作る：𝑥𝑥0|0 (𝑖𝑖) , 𝑖𝑖 = 1, 𝑁𝑁② 各個体の時間発展(非線形)： 𝑥𝑥𝑡𝑡|𝑡𝑡−1 (𝑖𝑖) = 𝑓𝑓𝑡𝑡 𝑥𝑥𝑡𝑡−1|𝑡𝑡−1 (𝑖𝑖) , 𝑣𝑣𝑡𝑡 , 𝑖𝑖 = 1, 𝑁𝑁③ 各粒子の尤度を計算する(非線形観測演算子も可)(𝑖𝑖)
𝜆𝜆𝑡𝑡 = 𝑝𝑝 𝒚𝒚𝑡𝑡 |𝒙𝒙𝑡𝑡|𝑡𝑡−1 (𝑖𝑖)=11
�𝑡𝑡|𝑡𝑡−1 (𝑖𝑖)exp − 𝒚𝒚𝑡𝑡 − ℎ𝑡𝑡 𝒙𝒙
2(2𝜋𝜋)𝑙𝑙 |𝑅𝑅𝑡𝑡 |𝑇𝑇

�𝑡𝑡|𝑡𝑡−1 (𝑖𝑖)
𝑅𝑅𝑡𝑡 −1 𝒚𝒚𝑡𝑡 − ℎ𝑡𝑡 𝒙𝒙④ 尤度の比に応じて復元抽出をして新しいアンサンブルを得る（尤度の高い個体は複製し，尤度の低い個体は消去する）⑤ ➁に戻って繰り返す
アンサンブルカルマンフィルタ
尤度予測分布

粒子フィルタ観測
分布

フィルタ分布

尤度予測
分布

観測
分布

フィルタ分布18

フィルム冷却（問題設定）
 単純化した流れ場を設定（1つのフィルム冷却孔に注目）

 壁面温度分布からコスト関数（尤度関数）を定義
コスト関数の例

𝐽𝐽𝑇𝑇 = �

12𝑇𝑇𝑐𝑐𝑐𝑐𝑐𝑐 − 𝑇𝑇𝑒𝑒𝑒𝑒𝑒𝑒219

SST乱流モデル(ANSYS Fluent)Turbulent kinetic energy

Specific dissipation rate

パラメータ変更時の壁面温度の収束
元の流れ場

パラメータ修正後

𝛽𝛽1 = 0.075 𝛽𝛽2 = 0.0828
𝜎𝜎𝑘𝑘𝑘 = 1.176 𝜎𝜎𝑘𝑘𝑘 = 1.0𝜎𝜎𝜔𝜔1 = 2.0 𝜎𝜎𝜔𝜔2 = 1.168
𝑎𝑎1 = 0.31 𝛽𝛽 ∗ = 0.09𝛽𝛽1 = 0.1397 𝛽𝛽2 = 0.1196
𝜎𝜎𝑘𝑘𝑘 = 1.2808 𝜎𝜎𝑘𝑘𝑘 = 0.8134
𝜎𝜎𝜔𝜔1 = 3.1768 𝜎𝜎𝜔𝜔2 = 2.0533
𝑎𝑎1 = 0.2885 𝛽𝛽 ∗ = 0.0934

20中間のまとめ
フィルム冷却流の乱流モデルパラメータ推定を行った

 少数の（時間不変）パラメータ推定には最適化手法
でよいはず パラメータの不確実性とモデル構造の不確実性を考える必要あり（本スライドでは結果省略）
 修正されたパラメータと物理現象を関連づけなければならない

21

データ同化の問題設定（DAシステムの評価）
可観測性と可制御性（一部，日立製作所との共同研究）

22能動的なデータ取得（適応型計測） 気象分野では台風予測などではAdaptive Observation（航空機でドロップゾンデをばらまく）
 UAV (Unmanned Aerial Vehicle, 自律飛行するセンサー)の利用可能性の増大（低価格化，汎用品化）
 “Trillion Sensors Universe”センサーリッチの社会へ(iPhone)http://mainichi.jp/articles/20160901/ddm/016/040/003000c

 小型センサーの爆発的普及が進んでいる
 数十億人が毎日数十個のセンサーを持ち歩くようになる（なっている）
 データ同化にとってもチャンス?
ビッグデータ（解析）：各種センサーで得た大量のデータから情報を探し当てる“ビッグ”データ同化：物理モデルを各種センサーで得た大量のデータに当てはめて現象を推定・予測する
23

データ同化に基づく適応型計測システムの開発

実証実験
Supercomputer
計測とトラバース制御

Labview PC
低乱熱伝達風洞@流体研
三坂孝志，鵜飼孝博，小西康郁，大林茂，データ同化に基づく計測
位置の最適化, 日本機械学会 第28回計算力学講演会, 2015.

24データ同化に基づく適応型計測の応用可能性
自律飛行ドローンを用いた計測データ駆動型の環境計測システム
自律飛行試験
Xbeeデータリンク

センサー試作

MissionPlanner(PC) ドローンの素早い展開を活かし，データ同化に基づく効率的な計測経路設定を行うアルゴリズ
ムを開発する（データ同化の応用として）
25

データ同化における計測・制御変数の設定
データ同化問題において，計測位置や制御量（コントロールパラメータ）の選択が状態推定結果にどのように影響するか?
二次元熱伝達問題の例
2

𝜕𝜕𝑇𝑇/𝜕𝜕𝑦𝑦 = 02

𝜕𝜕 𝑇𝑇 𝜕𝜕 𝑇𝑇+=0
𝜕𝜕𝑥𝑥 2 𝜕𝜕𝑦𝑦 2<数値計算法>- 二階微分  二次精度中心差分- SOR法のよる反復解法<データ同化手法>- アンサンブルカルマンフィルタ- アンサンブル数: 200<データ同化で推定する量>- 左側面の境界条件

温度分布を推定する
𝑇𝑇 = 𝑇𝑇(𝑦𝑦)𝑇𝑇 = 0𝜕𝜕𝑇𝑇/𝜕𝜕𝑦𝑦 = 026

二次元熱伝達のデータ同化例
赤丸：計測点
推定値（Case1）

疑似計測値
（真値）

推定値（Case2）

1.6
初期値
1.4

計測値
推定値（Case1）

温度

1.2

推定値（Case2）

1.0

計測点が左側面に近い方が温度を正確に推定できる

0.80.60.412345678910格子点

目的：計測（計測量，位置）や被制御量（状態推定変数）の善し悪しを系統的に検討する手法を構築する
27

線形システムの可制御性と可観測性
𝒙𝒙̇ = 𝐴𝐴𝒙𝒙 + 𝐵𝐵𝒖𝒖
𝒚𝒚 = 𝐶𝐶𝒙𝒙

𝒙𝒙: 状態ベクトル

𝒚𝒚: 観測ベクトル

𝒖𝒖: 外部入力ベクトル

𝐴𝐴, 𝐵𝐵, 𝐶𝐶: それぞれ行列

可制御性：有限時間で任意の状態から希望の状態へ移すことが可能な制御入力が存在するかどうかを示す指標
 可制御行列がフルランク

rank 𝐵𝐵, 𝐴𝐴𝐴𝐴, 𝐴𝐴2 𝐵𝐵, ⋯ , 𝐴𝐴𝑛𝑛−1 𝐵𝐵 = 𝑛𝑛𝑡𝑡

 可制御性グラム行列が正則 𝑊𝑊𝐶𝐶 = � 𝑒𝑒 −𝐴𝐴𝜏𝜏 𝐵𝐵 𝑒𝑒 −𝐴𝐴𝜏𝜏 𝐵𝐵 𝑇𝑇 𝑑𝑑𝜏𝜏
0

可観測性：有限時間だけ観測される値により状態を把握することができるかどうかを示す指標
 可観測行列がフルランク

rank 𝐶𝐶, 𝐶𝐶𝐶𝐶, ⋯ , 𝐶𝐶𝐴𝐴𝑛𝑛−1 𝑇𝑇 = 𝑛𝑛𝑡𝑡

 可観測性グラム行列が正則 𝑊𝑊𝑜𝑜 = � 𝑒𝑒 −𝐴𝐴𝜏𝜏 𝐶𝐶 𝑒𝑒 −𝐴𝐴𝜏𝜏 𝐶𝐶 𝑇𝑇 𝑑𝑑𝜏𝜏028

グラム行列の近似計算法
 線形システムのグラム行列
 前ページのように行列のランクから判定（YES or NO）
 非線形（大規模）システムのグラム行列
 経験的グラム行列（複数回の解析プログラム実行で評価）
経験的可制御性グラム行列 ・・・ 外部入力のインパルス応答（外部入力の数のrun）
経験的可観測性グラム行列 ・・・ 状態変数のインパルス応答（変数の数のrun）
アジョイント方程式の計測のインパルス応答（計測の数のアジョイントrun）
29

可観測性の検討手法
 システムの行列𝐴𝐴, 𝐵𝐵, 𝐶𝐶が直接わからない場合にも適用できる必要
 ブラックボックスコードに対応
 計算コストは可能な限り小さく 経験的可観測性グラム行列の近似評価∆𝒚𝒚𝑖𝑖 𝑡𝑡𝑗𝑗

King, S., Kang, W., Xu, K.: Observability foroptimal sensor locations in data assimilation,
Int. J. Dynam. Control, 3, 416-424, 2015.

1=𝒚𝒚 𝑡𝑡𝑗𝑗 , 𝜆𝜆, 𝒙𝒙0 + 𝜌𝜌𝒘𝒘𝑖𝑖 − 𝒚𝒚 𝑡𝑡𝑗𝑗 , 𝜆𝜆, 𝒙𝒙0 − 𝜌𝜌𝒘𝒘𝑖𝑖
2𝜌𝜌

∆𝒚𝒚𝑖𝑖 𝑡𝑡𝑗𝑗 : ベクトル𝒘𝒘𝑖𝑖 （直交系など，ここではPODモード）で状態変数ベクトルを変化させたときの観測ベクトルの差，𝜆𝜆は観測に関するパラメータ（位置など）

𝑌𝑌 = [∆𝒚𝒚1 𝑡𝑡1 , ⋯ , ∆𝒚𝒚1 𝑡𝑡𝑚𝑚 ⋯ ∆𝒚𝒚𝑝𝑝 𝑡𝑡1 , ⋯ , ∆𝒚𝒚𝑝𝑝 𝑡𝑡𝑚𝑚 �𝑡𝑡𝑗𝑗 = 𝑡𝑡1 ~𝑡𝑡𝑚𝑚

𝑡𝑡𝑗𝑗 = 𝑡𝑡1 ~𝑡𝑡𝑚𝑚

mステップ分のスナップショットをモード毎に並べたもの
𝐺𝐺𝑜𝑜 = 𝑌𝑌 𝑇𝑇 𝑌𝑌  固有値を調べる
 固有値の最大化することで計測位置を最適化

30経験的可観測性グラム行列の固有値
x方向に三カ所（x2, x4, x8）で，y方向に異なる点数の計測（○，□）を行った

計測点○0.7x20.6x8x40.5
固有値

𝜕𝜕𝑇𝑇/𝜕𝜕𝑦𝑦 = 00.4
0.3
0.20.1𝑇𝑇 = 𝑇𝑇(𝑦𝑦)0𝑇𝑇 = 001234モード

5678

計測点□

0.7x2

0.6x4x8

𝜕𝜕𝑇𝑇/𝜕𝜕𝑦𝑦 = 0 物理量の変動が大きいところを計測した方が固有値が大きくなる
 計測点数が多い方が固有値が大きい

固有値

0.5
0.4
0.3
0.2
0.1
0.0
01234
モード

567831

可観測性グラム行列の固有値とデータ同化結果
可観測性グラム行列の固有値

計測点の配置

0.7
Case3

0.6

Case3

Case40.5

固有値

Case40.40.30.20.1
0.001234
モード

5678

推定された壁面温度分布
1.6
初期値
1.4

計測値
推定値（Case3）

温度1.2

推定値（Case4）

RMSE

Case 3

Case 4（参考）Case 1

2.31e-1

1.35e-1

1.12e-2

1.00.8 固有値の大小と壁面温度の推定精度の相関

0.6
0.4
124

5678

910

格子点

32

後半のまとめ
適切なデータ同化の問題設定に関して，特に計測位置の最適化を簡単な事例で検討した
 経験的可観測性は評価が簡単なので実システムにも適用可能であり，観測方法の相対評価が可能
 勾配法を用いたためセンサー位置の最適化では初期値依存性がある（本スライドでは結果省略）
 計測点の変化に対応する必要あり

33