JP2014215822A - 状態推定装置、方法、及びプログラム - Google Patents

Abstract

【課題】精度よく推定対象の状態を推定することができる。
【解決手段】カルマンフィルタ部３２により、状態空間モデルで表わされるカルマンフィルタに従って、入力された観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、状態推定誤差を算出し、ロバストエンハンスメント部３４により、状態推定値に基づく観測値の残差の二乗の期待値を算出し、算出された期待値が予め定められた値以上であるときに、算出された期待値に基づいて、状態推定手段によって算出された状態推定値及び状態推定誤差の各々を補正し、出力部５０により、時系列データの観測値の各々について、カルマンフィルタ部３２による算出及びロバストエンハンスメント部３４による補正を繰り返したときに得られる状態推定値の各々を出力する。
【選択図】図２

Description

本発明は、状態推定装置、方法、及びプログラムに係り、特に、カルマンフィルタを用いて推定対象の状態を観測値から推定する状態推定装置、方法、及びプログラムに関する。

従来技術として、カルマンフィルタのロバスト性を確保するため、システムノイズ及び状態の推定誤差共分散行列を人工的に拡大させる方法がある（非特許文献１）。

D. Simon, "Optimal State Estimation: Kalman, H∞, and Nonlinear Approaches", John Wiley & Sons,Inc., Hoboken, New Jersey, 2006.

非特許文献１の方法では、いつ、どれ位、システムノイズ及び状態の推定誤差共分散行列を拡大させれば良いかが解析者の経験に基づいて決定されることが多く、フィルタの初心者にとってはパラメータ設定の敷居が高く扱いづらいという問題がある。

本発明では、上記問題点を解決するために成されたものであり、精度よく推定対象の状態を推定することができる状態推定装置、方法、及びプログラムを提供することを目的とする。

上記目的を達成するために、第１の発明に係る状態推定方法は、カルマンフィルタ手段と、ロバストエンハンスメント手段と、出力手段とを含み、入力された推定対象の観測値の時系列データに基づいて、前記推定対象の状態を推定する状態推定装置における状態推定方法において、前記カルマンフィルタ手段は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出し、前記ロバストエンハンスメント手段は、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正し、前記出力手段は、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する。

第２の発明に係る状態推定装置は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段と、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段と、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段と、を含んで構成されている。

第１の発明及び第２の発明によれば、カルマンフィルタ手段により、予め定められた状態空間モデルであらわされるカルマンフィルタに従って、入力された観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、状態推定誤差を算出し、ロバストエンハンスメント手段により、状態推定値に基づく観測地の残差の二乗の期待値を算出し、算出された期待値が予め定められた値以上であるときに、算出された期待値に基づいて、状態推定値及び状態推定誤差の各々を補正する。

そして、出力手段により、カルマンフィルタ手段による算出及びロバストエンハンスメント手段による補正を繰り返したときに得られる状態推定値の各々を出力する。

このように、第１の発明及び第２の発明によれば、精度よく推定対象の状態を推定することができる。

また、本発明のプログラムは、コンピュータを、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段、及び前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段、として機能させるためのプログラムである。

以上説明したように、本発明の状態推定装置、方法、及びプログラムによれば、精度よく推定対象の状態を推定することができる。

観測データと状態推定結果のデータ構造の例を示す図である。 本発明の実施の形態に係る状態推定装置の機能的構成を示すブロック図である。 本発明の実施の形態に係る状態推定装置における状態推定処理ルーチンを示すフローチャート図である。 カルマンフィルタの結果とロバストカルマンフィルタの推定結果の例を示す図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

＜発明の概要＞
本実施の形態では、入力されるデータは、観測機器によって観測された推定対象（現象）の観測データと、推定対象の背後にある支配方程式を解析者がモデル化した状態空間モデルである。そして、入力された観測データと状態空間モデルとに基づいて、カルマンフィルタ及びロバストエンハンスメントからなるロバストカルマンフィルタによる処理を行い、推定対象の状態を推定する結果を出力する。

本実施の形態において、フィルタが異常を検知した際、カルマンフィルタを通常の二乗誤差規範に対して準最適化する。これは、最小化する損失関数に罰則項を追加することで実現され、結果としてフィルタのモデリング誤差に対するロバスト性が向上する。一方、フィルタが正常な区間においては、通常の二乗誤差規範に対して最適であるカルマンフィルタと等価なアルゴリズムを適用する。

観測データは時系列データとなっており、観測データの次元をＭとする。状態空間モデルは解析対象の背後にある支配方程式を解析者がモデル化したものである。具体的には、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて定められ、下記（１）式及び（２）式に表される線形な状態空間モデルである。

ここで、確率変数Ｘ_ｔを時刻ｔにおけるＬ次元の状態（推定対象）、確率変数Ｙ_ｔを時刻ｔにおける観測値と定義すると、上記（１）式及び（２）式は、上記（１）式が状態Ｘ_ｔの時間発展となり、上記（２）式が状態Ｘ_ｔの観測機構を表す連立方程式となる。

上記（１）式及び（２）式のＦ、Ｈはそれぞれ状態遷移行列、観測出力行列と定義する。μ_ｔ−１、ω_ｔはシステムノイズ、及び観測ノイズと定義し、それぞれ平均０、共分散行列Ｑ及びＲの正規分布に従うノイズとする。状態空間モデルにおけるＦ，Ｈ、Ｑ、Ｒの値は解析者により設定する。

フィルタの問題は、観測値Ｙ_ｔ：１＝Ｙ_ｔ，Ｙ_ｔ−１，・・・，Ｙ_１を使用した状態Ｘ_ｔの推定であり、カルマンフィルタはその推定誤差の二乗を最小化する最適な推定値＾Ｘ_ｔを出力する。つまり、カルマンフィルタは下記（３）式の損失関数を最小化する推定値を出力する。

しかし、カルマンフィルタには、フィルタの回数が増加すると推定誤差の二乗が小さくなり過ぎ、推定誤差の共分散行列が小さくなり過ぎるために状態の時間変化に対する素早い追従性を失うという問題がある。推定対象は一般的に複雑な時間変化を示すので、状態空間モデルで事前に想定されたモデルの誤差もそれに伴い時間変化するが、フィルタの回数が大きくなるとその誤差に対してフィルタが敏感になってしまい、状態の時間変化に適切に追従することができなくなる。

＾Ｖ_ｔは、観測値Ｙ_ｔ：１が与えられた後の状態の推定誤差共分散行列であり、上記（３）式の二乗誤差はこの値により正規化されている。本実施の形態におけるロバストカルマンフィルタは下記（４）式及び（５）式の損失関数を最小化するように設計される。

は、予測残差の二乗の期待値（観測値の残差の二乗の期待値）であり、上記（４）式及び（５）式の損失関数は、状態の推定誤差は小さいが、フィルタの予測残差が大きくなっているときに通常の二乗誤差項（第一項）に罰則を与える形になっている。

従って、フィルタが異常を検知した場合には、状態の二乗誤差を最小化することを止め、第二項の付加によって準最適な推定値を出力するようにする。準最適な推定値の推定誤差共分散行列は、最適な推定値の推定誤差共分散行列より大きくなるため、結果として次のフィルタリングにおけるフィルタのロバスト性が自動的に回復する。

上記（５）式の予測残差の二乗は、￣Ｖ［ｖ_ｔ］（観測値Ｙ_{ｔ−１：１}が与えられた後の残差の予測共分散行列）で正規化されており、第一項の二乗誤差の正規化と合わせてある。

上記（４）式の第一項、第二項共にカイ二乗分布に従う確率変数の期待値になっており、第一項の期待値は状態の次元のＬ、第二項の期待値は観測値の次元のＭである。従って、両項のバランスをとるため、кを第二項にかける。

上記（４）式及び（５）式の損失関数を最小化するロバストカルマンフィルタの状態推定アルゴリズムは下記（６）〜（１３）式に示すようになる。

＜カルマンフィルタの状態推定アルゴリズム＞

＜準最適化によるロバストエンハンスメントの状態推定アルゴリズム＞

なお、上記の予測残差の二乗の期待値が上記（１３）式のγ_tに対応する。

カルマンフィルタには、状態空間モデルで設定された上記（１）式及び（２）式と時刻ｔにおける観測データの観測データＹ_ｔが入力される。この情報に基づいて、上記（６）〜（１０）式を用いて、状態推定値の初期値＾Ｘ´_０及び状態推定誤差の初期値＾Ｖ´_０又は、ロバストエンハンスメントにおいて前回算出された状態推定値＾Ｘ´_ｔ及び状態推定誤差＾Ｖ´_ｔを時間更新、及び観測更新させ、カルマンフィルタにおける状態推定値＾Ｘ_ｔと、カルマンフィルタにおける状態推定誤差＾Ｖ_ｔを算出する。なお、状態推定値及び状態推定誤差の初期値は予め適切な値を解析者が設定しておく。

また、ロバストエンハンスメントにはカルマンフィルタの状態推定結果である＾Ｘ_ｔと、＾Ｖ_ｔとが入力され、上記（１１）式〜（１３）式を用いて、準最適化によるロバストエンハンスメントで新たに計算される状態推定値＾Ｘ´_ｔと、状態推定誤差＾Ｖ´_ｔが状態推定結果として出力される。

具体的には、上記（１３）式においてγ_ｔを算出し、γ_ｔ≧τの場合には、カルマンフィルタの状態推定結果である＾Ｘ_ｔと、＾Ｖ_ｔとに基づいて、上記（１１）式及び（１２）式を用いて、状態推定値＾Ｘ_ｔを補正した状態推定値＾Ｘ´_ｔと、状態推定誤差＾Ｖ_ｔを補正した状態推定誤差＾Ｖ´_ｔとを求める。

また、γ_ｔ≦＜τの場合には、状態推定値＾Ｘ_ｔを＾Ｘ´_ｔとし、状態推定誤差＾Ｖ_ｔを状態推定誤差＾Ｖ´_ｔとする。

ここで、

であり、ｖ_ｉは時刻ｉにおける観測残差（観測値の予測誤差）のベクトルを表す。β_{ｔ，ｍｉｎ}は、＾Ｖ_ｔの最小の固有値であり、Ｍ´_{ｔ，ｍｉｎ}はと＾Ｖ_ｔの最小の固有値に対応する固有ベクトルを長さが１となるように正規化したベクトルである。τ及びＮの標準的な設定値（推奨値）は観測値が１次元の場合、

、Ｎ＝５〜１０である。観測値が２次元以上の場合もカイ二乗分布の有意水準５％の点を採用する。

＾Ｘ´_ｔの計算における±はどちらでもよいが、＋と−を交互に使うことで、ロバストカルマンフィルタの推定量のバイアス性を軽減する効果を得る事が出来る。

カルマンフィルタとロバストエンハンスメントは繰返しの処理関係になっており、時刻ｔの観測値Ｙ_ｔを得る度にカルマンフィルタ→ロバストエンハンストメントを繰り返し、ロバストエンハンスメントは状態の推定結果を出力する。この繰返しの処理体系全体をロバストカルマンフィルタと定義する。図１に観測データと状態推定結果のデータ構造を示す。

＜状態推定装置の構成＞
次に、本発明の実施の形態に係る状態推定装置の構成ついて説明する。図２に示すように、本発明の実施の形態に係る状態推定装置１００は、ＣＰＵとＲＡＭと後述する状態推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この状態推定装置１００は、機能的には図２に示すように入力部１０と、演算部２０と、出力部５０とを備えている。

入力部１０は、観測機器（図示省略）によって得られた、推定対象の次元がＭの時系列データである時刻ｔの観測データＹ_ｔの入力を受け付ける。

入力部１０は、解析者により入力される上記（１）式及び（２）式からなる状態空間モデルを受け付ける。また、入力部１０は、状態遷移行列Ｆ、観測出力行列Ｈ、共分散行列Ｑ及びＲ、の値の入力を受け付ける。

演算部２０は、カルマンフィルタ部３２と、ロバストエンハンスメント部３４とを備えている。

カルマンフィルタ部３２は、入力部１０により受け付けた各時刻ｔにおける観測データＹ_ｔと、状態空間モデルと、状態推定値の初期値＾Ｘ´_０及び状態推定誤差の初期値＾Ｖ´_０又は、ロバストエンハンスメント部３４において前回算出された状態推定値＾Ｘ´_ｔ及び状態推定誤差＾Ｖ´_ｔとに基づいて、カルマンフィルタにおける状態推定値＾Ｘ´_ｔと、カルマンフィルタにおける状態推定誤差＾Ｖ_ｔを算出する。

具体的には、上記（６）式と、（７）式と、（１０）式とにより，￣Ｘ_ｔと、￣Ｖ_ｔと、Ｋ_ｔとを求めた後に、上記（８）式及び（９）式により、＾Ｘ_ｔと、＾Ｖ_ｔとを求め、ロバストエンハンスメント部３４に出力する。

ロバストエンハンスメント部３４は、カルマンフィルタ部３２で算出された￣Ｖ_ｔ、￣Ｘ_ｔ、及び観測データＹ_ｔに基づいて、上記（１３）式に従って、γ_ｔを算出し、γ_ｔとτの値を比較する。γ_ｔ≧τの場合には、カルマンフィルタ部３２から入力された状態推定値＾Ｘ_ｔと、状態推定誤差＾Ｖ_ｔと、算出したγ_ｔとに基づいて、上記（１１）式、及び（１２）式に従って、補正後の状態推定値＾Ｘ´_ｔと、状態推定誤差＾Ｖ´_ｔを求めカルマンフィルタ部３２及び出力部５０に出力する。なお、β_{ｔ，ｍｉｎ}及びＭ´_ｔは、＾Ｖ_ｔの値から求められる。

また、γ_ｔ＜τの場合には、カルマンフィルタ部３２から入力された状態推定値＾Ｘ_ｔを補正せずに状態推定値＾Ｘ´_ｔとし、状態推定誤差＾Ｖ_ｔを補正せずに状態推定誤差＾Ｖ´_ｔとしてカルマンフィルタ部３２及び出力部５０に出力する。

＜状態推定装置の作用＞
次に、本発明の実施の形態に係る状態推定装置１００の作用について説明する。図３に示す状態推定処理ルーチンを実行する前に、入力部１０により状態空間モデルを受け付ける。そして、入力部１０により、各時刻ｔの観測データＹ_ｔを受け付けているときに、状態推定装置１００は、図３に示す状態推定処理ルーチンを実行する。

まず、ステップＳ１００では、入力部１０において受け付けた上記（１）式及び（２）式からなる状態空間モデルを読み込む。なお、同時に状態遷移行列Ｆ、観測出力行列Ｈ、共分散行列Ｑ及びＲの値も読み込む。

次に、ステップＳ１０２では、入力部１０において受け付けた現時刻ｔの観測データＹ_ｔを読み込む。

次に、ステップＳ１０４では、ステップＳ１００において読み込んだ状態空間モデルと、ステップＳ１０２において読み込んだ現時刻ｔにおける観測データＹ_ｔと、状態推定値の初期値＾Ｘ´_０及び状態推定誤差の初期値＾Ｖ´_０又は、ステップＳ１１０又はステップＳ１１２において前回算出された状態推定値＾Ｘ´_ｔ及び状態推定誤差＾Ｖ´_ｔとに基づいて、上記（６）式〜（１０）式に従って、カルマンフィルタにおける状態推定値＾Ｘ_ｔ及び状態推定誤差＾Ｖ_ｔを算出する。

次に、ステップＳ１０６では、ステップＳ１０４において、上記（６）式及び（７）式に従って算出された￣Ｖ_ｔ、￣Ｘ_ｔ、及びステップＳ１０２において取得した観測データＹ_ｔに基づいて、上記（１３）式を用いてγ_ｔを算出する。

次に、ステップＳ１０８では、ステップＳ１０６において算出したγ_ｔの値がτの値以上か否かの判定を行う。γ_ｔの値がτの値以上である場合には、ステップＳ１１０に移行し、γ_ｔの値がτの値未満の場合には、ステップＳ１１２に移行する。

ステップＳ１１０では、ステップＳ１０４において算出されたカルマンフィルタにおける状態推定値＾Ｘ_ｔ、状態推定誤差＾Ｖ_ｔ、及びステップＳ１０６において算出されたγ_ｔに基づいて、上記（１１）式及び（１２）式に従って状態推定値＾Ｘ_ｔを＾Ｘ´ｔに補正し、状態推定誤差＾Ｖ_ｔを＾Ｖ´_ｔに補正する。

ステップＳ１１２では、ステップＳ１０４において算出されたカルマンフィルタにおける状態推定値＾Ｘ_ｔを補正をせずに＾Ｘ´_ｔとし、状態推定誤差＾Ｖ_ｔを補正をせずに＾Ｖ´_ｔとする。

次に、ステップＳ１１４では、ステップＳ１１０又はステップＳ１１２において求められた状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔを推定結果としてカルマンフィルタ部３２及び出力部５０に出力し、ステップＳ１０２に移行し、次の時刻ｔの観測データＹ_ｔを読み込みステップＳ１０２〜ステップＳ１０８の処理を繰り返す。

＜実験例＞
本発明の状態推定装置におけるロバストカルマンフィルタの有効性を、以下（１４）式に示すＡＲ（２）過程の係数推定の実験で検証した。

上記（１４）式で示されるＡＲ（２）過程はｔ≦１０００の区間は定常で、１０００＜ｔ≦２０００の区間は独立過程になり、その後また元の定常過程に戻るといった構造の変化を有するものである。このＡＲ（２）過程を状態空間モデルで表すと下記（１５）式になる。

問題はこの２次元の状態ベクトルＸ_ｔのフィルタ推定であり、通常のカルマンフィルタの結果と共に、ロバストカルマンフィルタの推定結果を以下に説明する。

図４の左上が観測値の時系列で、右上が時変ＡＲ係数の真値をプロットしたものである。左下がカルマンフィルタの推定結果で、最初ｎ区間（ｔ≦１０００）の状態は正しく推定されているが、その後の係数の変化は追従できていないことがわかる。これは最初の１０００回のフィルタ後、状態の推定誤差共分散行列が極端に小さくなり、フィルタのロバスト性がほぼ無くなったためである。右下がロバストカルマンフィルタの結果で、ＡＲ係数の時間変化を素早く追従していることがわかる。推定値の準最適化が正しく動作し、必要な時にフィルタのロバスト性が適切にエンハンスされている。

上述した数値実験の結果もシステムモデルの時間変化に対する優れた追従性を見せていることから、本アルゴリズムの有効性が確認できた。

以上説明したように、本発明の実施の形態に係る状態推定装置によれば、システムモデルの時間変化に対する追従性を向上させて、精度よく推定対象の状態を推定することができる。

また、推定対象（現象）の状態を観測値から推定するフィルタリングの問題において、モデリング誤差に対してロバストなカルマンフィルタを提案することができる。

また、フィルタが異常を感知すると自動的にロバスト性を向上させる自己回復機能を持ち、状態の推定値を観測値に基づき素早く修正することができる。

また、必要な時にフィルタのロバスト化が適切に行われるため、解析者の熟練した経験を必要とせず、結果として安定した状態の推定結果を出力することができる。

なお、本発明は、上記の実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

上記の実施の形態では、出力部５０において、状態推定値＾Ｘ´_ｔ及び状態推定誤差＾Ｖ´_ｔを推定結果として出力する場合について説明したが、これに限定されるものではない。例えば、出力部５０において、状態推定値＾Ｘ´_ｔのみを出力するようにしてもよい。

１０ 入力部
２０ 演算部
２２ 入力部
３２ カルマンフィルタ部
３４ ロバストエンハンスメント部
５０ 出力部
１００ 状態推定装置

Claims (4)

1. カルマンフィルタ手段と、ロバストエンハンスメント手段と、出力手段とを含み、入力された推定対象の観測値の時系列データに基づいて、前記推定対象の状態を推定する状態推定装置における状態推定方法において、
   前記カルマンフィルタ手段は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出し、
   前記ロバストエンハンスメント手段は、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正し、
   前記出力手段は、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する
   状態推定方法。
2. 前記カルマンフィルタ手段は、以下の（１）式〜（５）式に従って、前記状態推定値^Ｘ_tを算出すると共に、前記状態推定誤差^Ｖ_tを算出し、
   前記ロバストエンハンスメント手段は、以下の（８）式に従って、前記観測値の残差の二乗の期待値γ_tを算出し、前記算出された前記期待値γ_tが予め定められた値τ以上であるときに、以下の（６）式、（７）式に従って、前記状態推定手段によって算出された前記状態推定値^Ｘ_t及び前記状態推定誤差^Ｖ_tの各々を状態推定値＾Ｘ´_ｔ及び状態推定誤差＾Ｖ´_ｔに補正する請求項１記載の状態推定方法。
   
   ただし、Ｆは状態遷移行列であり、Ｈは観測出力行列であり、Ｒ及びＱは共分散行列であり、ｖ_ｉは時刻ｉにおける前記観測値の残差を示すベクトルであり、β_{ｔ，ｍｉｎ}は＾Ｖ_ｔの最小の固有値であり、Ｍ´_ｔは＾Ｖ_ｔの最小の固有値に対応する固有ベクトルを長さが１となるように正規化したベクトルであり、Ｎは定数である。
3. システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段と、
   前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段と、
   前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段と、
   を含む、状態推定装置。
4. コンピュータを、
   システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段、
   前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段、及び
   前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段、
   として機能させるためのプログラム。