日本機械学会［No.167-1］北陸信越支部 第 53 期総会・講演会 講演論文集 [2016.3.5 長野県長野市]
［No.167-1］日本機械学会北陸信越支部 第 53 期総会・講演会 講演論文集 [2016.3.5 長野県長野市]
305
カルマンフィルタ有限要素法による
流れ場の推定解析における計測点配置に関する検討
吉荒 太一*1，倉橋 貴彦*2，小林 泰秀*2
Examinations for location of observation points in flow estimation analysis based on Kalman filter finite element method
Taichi YOSHIARA*1，Takahiko KURAHASHI*2 and Yasuhide KOBAYASHI*2
*1 Graduate School of Nagaoka University of Technology Kamitomiokamachi 1603-1, Nagaoka-shi, Niigata, 940-2188 Japan
The Kalman filter developed by R. E. Kalman and R. S. Bucy is states estimation theory in the target domain and has been employed in various field of engineering. In general, observation data in any system include observation noise and the computational model include system noise. These noises can be eliminated by using a Kalman filter. The finite element equation is applied to derive the state transition matrix in the system equation of the Kalman filter. As the fundamental study, the estimation of the flow field in open channel is carried out based on the finite element method and the Kalman filter theory. In addition, as the governing equation, the shallow water equations are employed, and the finite element and the selective lumping methods are applied to discretize the governing equations in space and time, respectively. The open channel model is employed in the numerical experiment, and some examinations are carried out by changing the observation variables, number and position of observation points.
Key Words : Finite element method, Kalman filter theory, Estimation flow field, Number of observation points, Number of variables, Position of observation points
1. 緒 言
カルマンフィルタ(1)(2)は1960年代にKalmanとBusyにより開発された観測値に基づく対象領域における状態推 定理論であり，様々な工学の分野で利用されている．カルマンフィルタの特徴として，誤差を含む計測値からシ ステム方程式を下に，対象領域の状態量を高精度に推定することができる点が挙げられる．また，空間的に広が りを持ち，時間進展するモデルへの適用例も多く，空間における状態量分布の離散化式を導く方法として境界要 素法や有限要素法等が用いられている．本研究では，カルマンフィルタ理論および有限要素法に基づき，ある水 域における流れ場の推定解析を行うことに焦点を当てる．これまでにもカルマンフィルタを用いて，水域内の流 れ場を推定する検討(3)が行われてきたが，計測点の数および位置が推定精度に与える影響について議論されてい るものは見当たらない．そこで推定計算において，使用する計測変数（流速および水位変動量），計測点の数およ び位置によっては良好な推定結果が得られるとは限らないため，これらの条件が推定精度に与える影響について 開流路を対象に数値実験を行い，考察する．数値モデルとして平面2次元場の線形浅水長波方程式を適用し，空 間近似には有限要素法に基づくガラーキン法，時間近似にはセレクティブ・ランピング法を適用した．検討から 得られた計測変数，計測点の数および位置が推定精度に与える影響についての知見は，実際の水域における流れ 場の推定を行う場合，例えば潮流発電や潮汐発電のための適切な地点選定や船舶の航路決定等に利用できると考 えられる．
*1 非会員，長岡技術科学大学大学院（〒940-2402 新潟県長岡市上富岡町1603-1）
*2 正員，長岡技術科学大学大学院 E-mail: [email protected]
2. 支配方程式の離散化とカルマンフィルタ
対象とする流れは浅水長波流れで表されるものとし，式（1）~（3）に示す平面2次元場の線形浅水長波方程式 を導入する．
𝑢̇ + 𝑔𝜂,𝑥 = 0 （1）
𝑣̇ + 𝑔𝜂,𝑦= 0 （2）
𝜂̇ + ℎ(𝑢,𝑥+ 𝑣,𝑦) = 0 （3）
ここに，𝑢，𝑣はx方向y方向流速，𝜂は平均水位からの水位変動量，𝑔は重力加速度，hは水深を表す．支配方程 式を離散化するにあたり，有限要素法を適用する．空間近似には三角形1次要素を用いたガラーキン法，時間近 似にはセレクティブ・ランピング法を用いることにより，以下のように離散化された有限要素方程式が得られる．
式(4)をすべての要素について足し合わせ，式（5）のように簡略的に記述する．
{ {𝑢𝑛+1} {𝑣𝑛+1} {𝜂𝑛+1}
} = [
[𝑀̅] [0] [0]
[0] [𝑀̅] [0]
[0] [0] [𝑀̅]
] [
[𝑀̃] [0] −𝑔Δ𝑡[𝑆𝑥] [0] [𝑀̃] −𝑔Δ𝑡[𝑆𝑦]
−ℎ̅Δ𝑡[𝑆𝑥] −ℎ̅Δ𝑡[𝑆𝑦] [𝑀̃]
] { {𝑢𝑛} {𝑣𝑛} {𝜂𝑛}
} （4）
{𝜙̂𝑛+1} = [𝐴]{𝜙̂𝑛} （5）
ここに，∆𝑡は微小時間間隔，ℎ̅は要素内の平均水深，[𝑆𝑥]，[𝑆𝑦]は圧力項に関するマトリックス，[𝑀]は質量マトリ
ックス，[𝑀̅]は質量マトリックスの各要素を対角項に集中させた行列，[𝑀̃]はランピングパラメータ𝑒(0 ≦ 𝑒 ≦ 1)
によって変形した行列（[𝑀̃] = (1 − 𝑒)[𝑀̅] + 𝑒[𝑀]）である．
また，カルマンフィルタの理論では，駆動行列（単位行列）[Γ]にシステムノイズによるベクトル{𝑞𝑛}を乗じた ベクトルを，式（5）に加えることにより真値を表すシステム方程式となる．また，システムノイズは平均0，分 散𝜎2により表した正規乱数により表されるものとすると，真値の時間進展の方程式は，式（6）のように与えられ る．ここに，行列[𝐴]は状態遷移行列と呼ばれる．また，観測値{𝑧𝑛+1}については真値{𝜙𝑛+1}に，観測ノイズ{𝑟𝑛+1} を加えた値により表されるものとすると，式（7）に示す観測方程式が得られる．観測ノイズもシステムノイズと 同様に，平均0，分散𝜎2により表した正規乱数により表されるものとする．行列[𝐻]は観測行列と呼ばれ，各観測 値に対して対応する真値の成分に1，それ以外には0が入る行列とする．また，カルマンフィルタ理論による状 態推定計算のアルゴリズムを以下に示す．
{𝜙𝑛+1} = [𝐴]{𝜙𝑛} + [Γ]{𝑞𝑛} （6）
{𝑧𝑛+1} = [𝐻]{𝜙𝑛+1} + {𝑟𝑛+1} （7）
1. 初期値の入力：[𝐴], [𝑃(+)0 ], {𝜙(+)0 }, [Γ], [𝑄], [𝑅], 𝜀, {𝑧𝑛+1}(𝑛 = 0~𝑖𝑚𝑎𝑥) 2. 予測誤差共分散行列の計算：[𝑃(−)] = [𝐴][𝑃(+)][𝐴]𝑇+ [Γ][𝑄][Γ]𝑇 3. カルマンゲイン行列の計算：[𝐾1] = [𝑃(−)][𝐻]𝑇([𝐻][𝑃(+)][𝐻]𝑇+ [𝑅])−1 4. 推定誤差共分散行列の計算：[𝑃(+)] = [𝑃(−)] − [𝐾1][𝐻][𝑃(−)]
5. 推定誤差共分散行列のフロベニウスノルムによる収束判定：
もし、 √𝑡𝑟 (([𝑃(+)𝑘+1] − [𝑃(+)𝑘 ])([𝑃(+)𝑘+1] − [𝑃(+)𝑘 ])𝑇) ≤ 𝜀 ならステップ6へ，そうでなければステップ2へ 6.支配方程式の有限要素方程式によるステップの更新：{𝜙̂(−)𝑛+1} = [𝐴]{𝜙̂(+)𝑛 }
7.同化後の最適推定値の計算：{𝜙̂(+)𝑛+1} = {𝜙̂(−)𝑛+1} + [𝐾1]({𝑧𝑛+1} − [𝐻]{𝜙̂(−)𝑛+1})
3. 観測変数の設定の違い・観測点の設置数および位置に関する数値実験
観測点における観測変数の設定を方向流速𝑢, 𝑣および水位変動量𝜂とした場合と, 水位変動量𝜂のみにした場合, また, 観測点を設置する数を変えた場合, 観測点の位置を変えた場合について検討を行う．流入境界から周期2𝜋 秒，振幅1.0mの水位変動を与え，システムノイズを考慮して得られた数値解析結果に観測ノイズを付加したもの を人工的な観測値として取り扱う．ここでシステムノイズ，観測ノイズはそれぞれ平均を0.0，分散を0.0001 と した正規ホワイトノイズ，平均を0.0，分散を0.1000とした正規ホワイトノイズとする．図1に検討条件，表1に 計算条件をそれぞれ整理する．
Fig.1 Numerical test conditions.
ここに，各検討ケースにおいて，観測変数の設定を方向流速𝑢, 𝑣および水位変動量𝜂とした場合を「𝑎」，水位変 動量𝜂のみとした場合を「𝑏」と付記する．駆動行列[Γ]は単位行列とする．また，システム誤差共分散行列[𝑄]は
0.0001，観測誤差共分散行列[𝑅]は0.1000を対角成分に持つ対角行列とする．図2にT=2.0s後における流路中心
線上における水位変動量の分布をそれぞれ示す．Case 1，2では，良好な水位変動量の推定結果が得られている．
Case3 では，計測点の配置されている流入境界側においては良好な推定が行われ，計測点の設置されていない流
出境界側でも傾向を捉えた推定が行われている．Case 4では，流入境界から流路中央までの観測点を設置してい ない区間で水位が一定となっていることがわかる．Case 5では，流路壁面に計測点を設置しても良好に推定が行 われていることが見て取れる．図3にCase 1aにおける流路中央点（𝑥=5.0m，𝑦=0.4m）における水位の時系列変 化を示す．また，表2に各ケースの流路中央点における真値と推定値の差に関するL2ノルムを示す．L2ノルムに よる比較を行うことで，各時間ステップにおける水位変動量の真値と推定値の一致率について検証する．Case 3a，
Case 3b以外は「a」のケースの方が「b」のケースに比べて，L2ノルムは小さくなることがわかる．
Case 1a and 1b
Time increment ∆𝑡, sec. 1.0 × 10−3
Time steps 2000
Number of nodes 153
Number of elements 200
Gravitational acceleration 𝑔, 𝑚/𝑠𝑒𝑐2 9.8
Lumping parameter 𝑒 0.8
Initial of estimated error covariance 𝑃(+)0 1.0 Initial of estimated value 𝜙̂(+)0 0 Convergence determination constant ε 10−2
Table 1 Computational conditions
Case 2a and 2b
Case 3a and 3b
Case 4a and 4b
Case 5a and 5b
Fig.3 Distribution of estimated water elevation at 𝑇 = 2.0s.
Fig.3 Time history of estimated water elevation at center point of channel in Case 1a.
4. 結 語
水域における観測変数の設定，観測点の数および位置を変え，カルマンフィルタ有限要素法に基づき，浅水域 における流れ場の推定精度について検証を行った．対象とする流れは浅水長波方程式より表し，有限要素法に基 づくガラーキン法，セレクティブ・ランピング法による離散化を行った．数値実験では，開流路モデルにおいて，
システムノイズベクトル{𝑞}および観測ノイズベクトル{𝑟}を正規ホワイトノイズと仮定し，シミュレーションに よる結果に付加した値を人工的な観測値として，推定計算に用いた．この結果，以下の所見が得られた．
（1）検討ケースCase 1a～Case 5bの流路中央点における水位の時系列変化において，真値と推定値の差に関す るL2ノルムを計算すると，Case 3a，Case 3bはおおむね等しいが，それ以外は「a」のケース（計測変数 を𝑢, 𝑣, 𝜂とした場合）の方が「b」のケース（計測変数を𝜂のみにした場合）に比べて，小さくなることが わかった．
（2）Case3では，観測点を設定しない流路中央から流出境界までの間の水位分布は定性的に真値に合う結果が
得られることを確認した．また，Case 4では，観測点を設定しない流入境界から流路中央まで間では，お おむね一定の水位となり，定性的にも推定が困難な結果となった．
（3）本研究では，1次元の流れ場を対象としたため，観測点の設置位置を流路壁面に設定した場合(Case 5)にお いても，流路中央部に観測点を設定した場合と同様に結果が得られるか確認を行ったところ，Case 1と同 様の結果を得ることを確認した．
謝 辞
本論文を執筆するにあたり，科学研究費補助金（基盤(C)）15K05786 の援助を受けた．ここに謝意を表す．
文 献
(1) Kalman, R.E., A New Approach to Linear Filtering and Prediction Problems, Transactions for the ASME – Journal of Basic of Engineering, Vol. 82, No.1(1960), pp.35-45.
(2) Kalman, R. E. and Bucy, R. S., New Results in Linear Filtering and Prediction Theory, Transactions of the ASME - Journal of Basic of Engineering, Vol. 83, No. 1 (1961), pp.95–108.
(3) Heemink, A. W., Two-dimensional Shallow Water Flow Identification, Applied Mathematical Modelling, Vol. 12(1988), pp.109-118.
Case L2 norm for difference of true and estimated values
1a 1.899
1b 1.984
2a 2.339
2b 4.050
3a 2.531
3b 2.524
4a 4.161
4b 4.988
5a 2.126
5b 2.206
Table 2 Comparison L2 norm for difference of true and estimated values.
