技術資料
Feel&Think
第5回 有限要素法とカルマンフィルタ
今回からは、いよいよカルマンフィルタを使った有限要素法による逆解析に挑戦します。この場合は、状態空間モデルやフィルタ方程式が多次元となるため、もう一度だけカルマンフィルタの方法論を説明します。
回目の観測において、観測値
(m次ベクトル)が状態
(n次ベクトル)と次のような関係があり、また状態は次のような変動をするものとします(状態空間モデル)。
(観測方程式)

(状態方程式)

ここに、は観測行列(m×n行列)、
は観測ノイズ(m次ベクトル)、
はシステムノイズ(r次ベクトル)とします。また、
は観測ノイズの共分散行列(m×m行列)、
はシステムノイズの共分散行列(r×r行列)、
は状態の予測誤差共分散(n×n行列)とします。
これらの条件の下で、状態の推定量を算出するカルマンフィルタは次のようになります。
(フィルタ方程式)

(初期条件)

ただし、は状態の初期値
の共分散とします。
一方、変形問題では有限要素方程式が次のような形をとります。

ここに、Fは荷重ベクトル、uは変位ベクトル、Kは剛性行列です。
変位は有限要素方程式に従い、これが状態(例えば、ヤング率やポアソン比)の関数とします。

観測値は、有限要素方程式によって求まる変位の近辺にあるものとし、ノイズを加えて観測方程式を作ります。

上式は非線形であるため、そのままではカルマンフィルタにあてはめることができないため線形化を施します。
まず、非線形関数まわりでTaylor展開し2次項以降を無視し、

新しい観測値として次式を定義すれば、

最終的に観測方程式ができあがります。

ちなみに、

を求めてみると、有限要素方程式を状態で偏微分し、


が得られます。荷重は状態の影響を受けないと考えれば次式が得られます。

この方程式を解くことで、

が得られます。
計算の手順は次のようになります。
(1)初期状態を設定します。

(2)予測誤差共分散行列の初期値と状態および観測値の分散
、
と状態および観測値の分散

(3)有限要素方程式を解き、変位を求めます。

(4)観測行列を得るため、次の方程式を解きます。

(5)予測誤差共分散行列を求めます。

(6)カルマンゲインを求めます。

(7)観測値をもとに状態の推定値を求めます。

(本来はですが、
を考慮)
(8)予測誤差共分散行列を補正します。

(9)として(3)に戻って繰り返します。
これらの過程を実現するには、複雑な数学的手法とプログラミング技術が必要です。地層科学研究所では、カルマンフィルタを用いた有限要素法による逆解析を広く使っていただくため、Geo-Inverse*1と名付けたソフトウェアを開発しました。次回以降は、このソフトを使っていろいろな逆解析に挑戦します。
*1:受託解析に使用するソフトウェアです。業務実績についてはお問い合わせください。
Geo-Inverseに関するテクニカルレポート(TR47 1-4)も併せてご参照ください。