"システムを観測して観測値を各々出力する複数のセンサと、前記システムの状態を表すための複数の状態変数を要素とする状態ベクトルと、前記複数のセンサから各々出力される観測値を要素とする観測値ベクトルとを用いて、前記状態ベクトルを更新することで、前記システムの状態を推定する複数のカルマンフィルタと、互いに異なる複数の初期ベクトルを生成し、複数の前記初期ベクトルのそれぞれを、前記状態ベクトルの初期値として前記複数のカルマンフィルタの各々に供給する初期値生成部と、前記複数のカルマンフィルタの推定精度を各々評価し、前記複数のカルマンフィルタのうち、推定精度が最も高いカルマンフィルタを特定し、特定したカルマンフィルタを除く他のカルマンフィルタの動作を停止させるカルマンフィルタ制御部と、を備えることを特徴とする状態推定装置。前記カルマンフィルタは、前記状態ベクトルと、前記観測値ベクトルとに基づいて観測残差を算出するとともに、前記状態ベクトルの推定誤差の共分散を算出し、前記カルマンフィルタ制御部は、前記複数のカルマンフィルタの各々が出力する前記観測残差に基づいて、前記状態ベクトルと前記システムの状態との誤差を評価する誤差評価値が最小となるカルマンフィルタを選択する観測残差比較部と、前記観測残差比較部によって選択されたカルマンフィルタが出力する前記状態ベクトルの推定誤差の共分散に基づいて、前記状態ベクトルの収束の程度を表す収束評価値を算出し、前記収束評価値が閾値以下であるか否かを判定する収束度判定部と、前記収束度判定部の判定結果が肯定である場合、前記観測残差比較部によって選択されたカルマンフィルタを、推定精度が最も高いカルマンフィルタであると特定し、特定したカルマンフィルタを除く他のカルマンフィルタの動作を停止させる、停止制御部と、を備えることを特徴とする、請求項1に記載の状態推定装置。前記観測残差比較部は、前記観測残差のうち少なくとも一部の要素のノルムの、現在時刻から所定時間だけ過去の時刻までの所定期間における最大値を、前記複数のカルマンフィルタごとに算出し、算出した複数の前記最大値が最小となるカルマンフィルタを、前記誤差評価値が最小となるカルマンフィルタとして選択する、ことを特徴とする、請求項2に記載の状態推定装置。前記収束評価値は、前記状態ベクトルの推定誤差の共分散を表す行列のトレースである、ことを特徴とする、請求項2または3に記載の状態推定装置。前記カルマンフィルタ制御部が、前記他のカルマンフィルタを停止させる前は、前記複数のカルマンフィルタは、第1のデータ型を用いて演算を行い、前記カルマンフィルタ制御部が、前記他のカルマンフィルタを停止させた後は、前記カルマンフィルタ制御部によって特定されたカルマンフィルタは、前記第1のデータ型よりもデータサイズの大きな第2のデータ型を用いて演算を行うことを特徴とする、請求項1乃至4のうちいずれか1項に記載の状態推定装置。ところで、観測残差にはノイズが重畳しているため、ある時刻で選択したカルマンフィルタが、次の時刻において最も観測残差が小さくなるとは限らない。従って、従来の技術のように、単に観測残差が最小になるカルマンフィルタを選択すると、不適切なものを選択してしまう可能性があるといった問題が存在する。さらに、従来の技術のように、複数のカルマンフィルタを並列で動作させる場合、状態ベクトルを収束させる時間を短縮することができるが、処理負荷が増加する。すなわち、状態ベクトルが収束するまでの時間と処理負荷とはトレードオフの関係にあるといった問題があった。本発明は、上述した点に鑑みてなされたものであり、状態ベクトルを適切な値に高速に収束させるとともに、非線形カルマンフィルタの演算に係る処理負荷を低減させることを解決課題とする。上述した課題を解決するため、本発明に係る状態推定装置は、システムを観測して観測値を各々出力する複数のセンサと、前記システムの状態を表すための複数の状態変数を要素とする状態ベクトルと、前記複数のセンサから各々出力される観測値を要素とする観測値ベクトルとを用いて、前記状態ベクトルを更新することで、前記システムの状態を推定する複数のカルマンフィルタと、互いに異なる複数の初期ベクトルを生成し、複数の前記初期ベクトルのそれぞれを、前記状態ベクトルの初期値として前記複数のカルマンフィルタの各々に供給する初期値生成部と、前記複数のカルマンフィルタの推定精度を各々評価し、前記複数のカルマンフィルタのうち、推定精度が最も高いカルマンフィルタを特定し、特定したカルマンフィルタを除く他のカルマンフィルタの動作を停止させるカルマンフィルタ制御部と、を備えることを特徴とする。カルマンフィルタは、状態ベクトルと観測値ベクトルとを用いて、状態ベクトルを更新する演算を繰り返すことにより、状態ベクトルの各要素を真値に近い値に収束させ、システムの状態を推定する。従って、状態ベクトルの初期値すなわち、初期ベクトルの要素が真値から大きく乖離する場合には、状態ベクトルの各要素が真値に近い値に収束するまでに長時間を要することがあり、更には、状態ベクトルの各要素が真値とは異なる不正確な値に収束することも起こり得る。このような初期ベクトルが供給されるカルマンフィルタは、推定精度が低く、システムの状態を正確に推定することはできない。この発明によれば、状態推定装置は、複数のカルマンフィルタに対して、互いに異なる複数の初期ベクトルをそれぞれ供給し、これら複数のカルマンフィルタを並列で動作させる。状態推定装置が互いに異なる複数の初期ベクトルを生成する場合、1つの初期ベクトルを生成する場合に比べて、真値に近い値を要素とする初期ベクトルが含まれる可能性が高い。真値に近い値を要素とする初期ベクトルが供給されるカルマンフィルタは、状態ベクトルの各要素が正確且つ高速に真値に収束させることができるため、システムの状態の推定精度も高い。そして、状態推定装置は、複数のカルマンフィルタの中から、推定精度が最も高いカルマンフィルタを特定し、特定したカルマンフィルタにより状態推定を行うとともに、特定したカルマンフィルタ以外のカルマンフィルタを停止する。従って、本発明に係る状態推定装置によれば、正確且つ高速な状態推定が可能であるとともに、状態推定に係る処理負荷の低減が可能となる。また、上述した状態推定装置において、前記カルマンフィルタは、前記状態ベクトルと、前記観測値ベクトルとに基づいて観測残差を算出するとともに、前記状態ベクトルの推定誤差の共分散を算出し、前記カルマンフィルタ制御部は、前記複数のカルマンフィルタの各々が出力する前記観測残差に基づいて、前記状態ベクトルと前記システムの状態との誤差を評価する誤差評価値が最小となるカルマンフィルタを選択する観測残差比較部と、前記観測残差比較部によって選択されたカルマンフィルタが出力する前記状態ベクトルの推定誤差の共分散に基づいて、前記状態ベクトルの収束の程度を表す収束評価値を算出し、前記収束評価値が閾値以下であるか否かを判定する収束度判定部と、前記収束度判定部の判定結果が肯定である場合、前記観測残差比較部によって選択されたカルマンフィルタを、推定精度が最も高いカルマンフィルタであると特定し、特定したカルマンフィルタを除く他のカルマンフィルタの動作を停止させる、停止制御部と、を更に備えることが好ましい。この発明によれば、状態推定装置は、誤差評価値と収束評価値とを用いて、推定精度が最も高いカルマンフィルタを特定する。誤差評価値は、観測残差に基づいて定められ、状態ベクトルの示す値と真値との誤差を表す値である。すなわち、誤差評価値が小さい場合、状態ベクトルの各要素は真値に近い値を示す。従って、誤差評価値が最小となるカルマンフィルタは、システムの状態を最も正確に推定するカルマンフィルタであると看做すことができる。一方、収束評価値は、状態ベクトルの推定誤差の共分散に基づいて定められるものであり、状態ベクトルの収束の程度を表す値である。状態ベクトルが収束している場合、ある時刻において状態ベクトルの各要素が真値に近い値を示していれば、当該ある時刻以降においても状態ベクトルの各要素は真値に近い値を示す可能性が高い。すなわち、収束評価値が小さい場合、カルマンフィルタは安定的にシステムの状態を推定することができる。このような、誤差評価値と、収束評価値とを用いることにより、ある時刻においてシステムの状態を最も正確に推定し、且つ、当該ある時刻以降も安定的にシステムの状態を推定することができるカルマンフィルタを特定することができる。すなわち、本発明に係る状態推定装置は、推定精度の最も高いカルマンフィルタを正確に特定し、正確且つ安定的な状態推定を行うことが可能になる。また、上述した状態推定装置において、前記観測残差比較部は、前記観測残差のうち少なくとも一部の要素のノルムの、現在時刻から所定時間だけ過去の時刻までの所定期間における最大値を、前記複数のカルマンフィルタごとに算出し、算出した複数の前記最大値が最小となるカルマンフィルタを、前記誤差評価値が最小となるカルマンフィルタとして選択する、ことを特徴とすることが好ましい。誤差評価値は、観測残差に基づいて算出される値である。しかし、観測残差は、複数のセンサから各々出力される観測値を要素とする観測値ベクトルに基づいて算出されるベクトルであるため、ノイズが重畳する。この発明によれば、観測残差のうち少なくとも一部の要素のノルムについて、所定期間における最大値を誤差評価値とする。従って、観測残差が観測ノイズ等により変動する場合であっても、ノイズによる影響を小さく抑えて誤差評価値を算出することができ、状態ベクトルの示す値と真値との誤差の大きさを正確に把握することが可能となる。また、上述した状態推定装置において、前記収束評価値は、前記状態ベクトルの推定誤差の共分散を表す行列のトレースである、ことを特徴とすることが好ましい。一般的に共分散行列の固有値は、当該固有値に対応する固有ベクトル方向の分散の大きさを表す値である。従って、状態ベクトルの推定誤差の共分散のトレースすなわち、共分散行列が有する複数の固有値の和は、状態ベクトルの推定誤差の大きさを表す。また、カルマンフィルタは、状態ベクトルを、初期状態から現在時刻に至る全ての状態ベクトルの値を考慮して更新する。従って、状態ベクトルの推定誤差の共分散も、初期状態から現在時刻に至る全ての共分散の値を考慮して更新される。よって、状態ベクトルの推定誤差の共分散のトレースが小さな値となる場合、状態ベクトルの推定誤差が小さな値で安定している状態すなわち、状態ベクトルが収束している状態であると看做すことができる。本発明は、収束評価値として、状態ベクトルの推定誤差の共分散のトレースを採用したため、状態ベクトルが集束している状態であるか否かを正確に判定することが可能となる。また、上述した状態推定装置において、前記カルマンフィルタ制御部が、前記他のカルマンフィルタを停止させる前は、前記複数のカルマンフィルタは、第1のデータ型を用いて演算を行い、前記カルマンフィルタ制御部が、前記他のカルマンフィルタを停止させた後は、前記カルマンフィルタ制御部によって特定されたカルマンフィルタは、前記第1のデータ型よりもデータサイズの大きな第2のデータ型を用いて演算を行うことを特徴とすることが好ましい。この発明によれば、複数のカルマンフィルタが並列で動作するときは、これら複数のカルマンフィルタは、データサイズの小さな第1のデータ型を用いた演算を行う。従って、本発明に係る状態推定装置は、複数のカルマンフィルタを並列で動作させる場合の計算負荷を、低減することができる。一方、この発明によれば、特定されたカルマンフィルタが単独で動作するときは、当該特定されたカルマンフィルタは、データサイズの大きな第2のデータ型を用いた演算を行う。従って、本発明に係る状態推定装置は、特定された1つのカルマンフィルタにより演算を行う場合に、高い精度でシステムの状態を推定することが可能となる。<A.実施形態>本発明の実施の形態について図面を参照して説明する。[1.機器構成及びソフトウェア構成]図1は、本発明の実施形態に係る携帯機器1のブロック図であり、図2は携帯機器1の外観を示す斜視図である。携帯機器1は、携帯機器1の姿勢の変化に応じて、携帯機器1が備える画面に表示される地図などの画像を回転させ、画像によって表される方位を現実空間の方位に追従させる機能を備える。この機能は、携帯機器1が備える各種センサの出力に基づいてカルマンフィルタの演算を行い、携帯機器1の姿勢または携帯機器1から見た地磁気の向き等を推定することによって実現される。携帯機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM蓄積部20、状態推定プログラム100等の各種のプログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及びGPS部60を備える。また、携帯機器1は、地磁気等の磁気を検出して磁気データqを出力する3次元磁気センサ70、加速度を検出して加速度データaを出力する3次元加速度センサ80、及び角速度を検出して角速度データgを出力する3次元角速度センサ90を備える。表示部50は、CPU10が状態推定プログラム100を実行することにより推定したシステムの状態に基づいて算出された地磁気の向きや携帯機器1の姿勢μ等を、矢印等の画像により表示する。GPS部60は、GPS衛星からの信号を受信して携帯機器1の位置情報緯度、経度を生成する。3次元磁気センサ70は、地磁気を検出するためのセンサであり、X軸磁気センサ71、Y軸磁気センサ72、及びZ軸磁気センサ73を備える。各センサは、MI素子磁気インピーダンス素子、MR素子磁気抵抗効果素子などを用いて構成することができる。磁気センサI/F74は、各センサからの出力信号をAD変換して磁気データqを出力する。この磁気データqは、x軸、y軸およびz軸の3成分を有するベクトルとして表されるデータである。より具体的には、磁気データqは、3次元磁気センサ70と一体となって運動する携帯機器1に固定された座標系において、X軸磁気センサ71からの出力値がx軸成分として表され、Y軸磁気センサ72からの出力値がy軸成分として表され、Z軸磁気センサ73からの出力値がz軸成分として表されるベクトルデータである。なお、3次元磁気センサ70は、携帯機器1に設けられるため、地磁気の他に、携帯機器1の部品が発する内部磁界と、携帯機器1の外部にある外部物が発する外部磁界とを検出する。但し、本実施形態では、外部磁界は無視できる程度に小さいものと想定する。従って、磁気データqを用いて地磁気の向き及び大きさを算出するためには、磁気データqから、内部磁界を表す成分を減算する補正処理を行う必要がある。この補正処理において、減算される値つまり、内部磁界を表す成分を、3次元磁気センサ70のオフセットと呼ぶ。3次元加速度センサ80は、X軸加速度センサ81、Y軸加速度センサ82、及びZ軸加速度センサ83を備える。各センサは、ピエゾ抵抗型、静電容量型、熱検知型などどのような検知方式であってもよい。加速度センサI/F84は、各センサからの出力信号をAD変換して加速度データaを出力する。この加速度データaは、3次元加速度センサ80と一体となって運動する携帯機器1に固定された座標系における慣性力と重力との合力を、x軸、y軸及びz軸の3成分を有するベクトルとして表すデータである。従って、携帯機器1が静止状態または等速直線運動状態にあれば、加速度データaは、携帯機器1に固定された座標系から見た重力加速度の大きさと方向とを示すベクトルデータとなる。3次元角速度センサ90は、X軸角速度センサ91、Y軸角速度センサ92、及びZ軸角速度センサ93を備える。角速度センサI/F94は、各センサからの出力信号をAD変換して角速度データgを出力する。角速度データgは、各軸の回りの角速度を示すベクトルデータである。CPU10は、ROM30に格納されている状態推定プログラム100を実行することによって、携帯機器1の姿勢や磁気センサのオフセット等のシステムの状態を表す複数の物理量を推定する。すなわち、携帯機器1は、CPU10が状態推定プログラム100を実行することで、状態推定装置として機能する。図3は、状態推定装置のうち、CPU10が状態推定プログラム100を実行することにより実現される機能を表す、機能ブロック図である。図3に示すように、状態推定装置は、互いに異なる複数の初期ベクトルを出力する初期値生成部200、複数のカルマンフィルタを備えカルマンフィルタの演算を行うカルマンフィルタ部300、カルマンフィルタ部300の動作を制御するカルマンフィルタ制御部400、及び、カルマンフィルタ部300からの出力に基づいて出力情報を生成する出力情報生成部500を備える。カルマンフィルタ部300は、10個のカルマンフィルタKFを備え、各々のカルマンフィルタKFは、非線形カルマンフィルタの演算を行う。なお、説明の便宜上、10個のカルマンフィルタKFを区別する場合、10個のカルマンフィルタを、それぞれ、KF[1]、KF[2]、...、KF[10]と表記する。また、10個のカルマンフィルタの中で、第m番目mは、1≦m≦10を満たす自然数のカルマンフィルタを、KF[m]と表記する。10個のカルマンフィルタKF[1]~KF[10]のそれぞれは、第1演算モードまたは第2演算モードのいずれかの演算モードで、非線形カルマンフィルタの演算を行う。ここで、第1演算モードとは、CPU10及び状態推定プログラム100が対応可能な実数を表すデータ型の中で、データサイズの最も小さな第1のデータ型を用いて、非線形カルマンフィルタの演算を行うモードである。また、第2演算モードとは、第1のデータ型よりもデータサイズの大きな第2のデータ型を用いて、非線形カルマンフィルタの演算を行うモードである。すなわち、第2演算モードは、第1演算モードと比べ、演算の精度が高く、処理負荷が大きい演算モードである。非線形カルマンフィルタの演算とは、以下で説明する状態ベクトルx及び観測値ベクトルyを用いて算出される観測残差zに基づいて、状態ベクトルxを周期的に更新することにより、システムの状態を推定する演算である。なお、詳細は後述するが、観測残差zとは、観測値ベクトルyと、推定観測値ベクトルy-との差分である。また、非線形カルマンフィルタの演算は、状態ベクトルの更新と同じ周期で、状態ベクトルxの推定誤差の共分散Pも更新する。状態ベクトルxとは、複数の状態変数を要素とするベクトルである。状態ベクトルxの要素である複数の状態変数は、それぞれ、非線形カルマンフィルタが推定するシステムの状態を示す複数の物理量を表す。すなわち、状態ベクトルxは、非線形カルマンフィルタが推定するシステムの状態を表すベクトルである。本実施形態は、状態ベクトルxの要素状態変数として、携帯機器1の姿勢μ、地磁気の強さr、地磁気の伏角φ、携帯機器1の角速度ω、角速度センサ91~93のオフセット推定値bO、及び磁気センサ71~73のオフセット推定値cOを採用し、これらを非線形カルマンフィルタの演算における推定の対象とする。時刻T=kにおける状態ベクトルxkは、以下の式1で表される。なお、各状態変数の右下に付された添え字kは、当該状態変数が時刻T=kにおける値であることを表す。ここで、地磁気の強さrはスカラーであり、地磁気の伏角φはスカラーであり、角速度ωは3次元ベクトルであり、角速度センサのオフセット推定値bOは3次元ベクトルであり、磁気センサのオフセット推定値cOは3次元ベクトルである。また、姿勢μは、クォータニオンを用いて表現される。クォータニオンとは、物体の姿勢回転状態を表す4次元数である。例えば、携帯機器1に固定された座標系の各軸と、地上に固定された座標系の各軸とが一致する姿勢μを基準姿勢と定義し、基準姿勢をμ=0、0、0、1と表現する。このとき、携帯機器1の任意の姿勢μは、基準姿勢に対して単位ベクトルρの方向を回転軸として角度ψだけ回転した姿勢として表現できる。回転後の姿勢μは、クォータニオンを用いて、式2で与えられる。一方、観測値ベクトルyとは、3次元磁気センサ70から出力される磁気データq、3次元加速度センサ80から出力される加速度データa、及び、3次元角速度センサ90から出力される角速度データgを要素とするベクトルである。時刻T=kにおける観測値ベクトルykは式3で与えられる。各々のカルマンフィルタKFは、非線形カルマンフィルタの演算の結果、周期的に算出される、状態ベクトルxk、観測残差zk、及び状態ベクトルの推定誤差の共分散Pkを、出力する。なお、説明の便宜上、第m番目のカルマンフィルタKF[m]からの出力は、状態ベクトルxk[m]、観測残差zk[m]、状態ベクトルの推定誤差の共分散Pk[m]等のように、符号[m]を付して表記する。また、詳細は後述するが、本実施形態では、各々のカルマンフィルタKFが出力する状態ベクトルxkとは、厳密には、更新後の状態ベクトルxkであり、各々のカルマンフィルタKFが出力する状態ベクトルの推定誤差の共分散Pkとは、厳密には、更新後の状態ベクトルの推定誤差の共分散Pkである。初期値生成部200は、ROM30に格納された各種設定情報等に基づいて、互いに異なる10個の初期ベクトルINI[1]~INI[10]を生成し、これらを出力する。10個の初期ベクトルINI[1]~INI[10]は、それぞれ、10個のカルマンフィルタKF[1]~KF[10]に供給され、状態ベクトルxの初期値すなわち、時刻T=0における状態ベクトルx0として適用される。第m番目の初期ベクトルINI[m]は、姿勢μの初期値μ0[m]、地磁気の強さrの初期値r0[m]、地磁気の伏角φの初期値φ0[m]、角速度ωの初期値ω0[m]、角速度センサのオフセット推定値bOの初期値bO0[m]、及び磁気センサのオフセット推定値cOの初期値cO0[m]を要素として含み、以下の式4により示される。なお、初期ベクトルINI[m]の各要素の生成方法については、後述する。カルマンフィルタ制御部400は、10個のカルマンフィルタKF[1]~KF[10]からの出力値に基づいて、10個のカルマンフィルタKF[1]~KF[10]の中で、推定精度が最も高いカルマンフィルタKF[Sel]を選択し、選択したカルマンフィルタKF[Sel]以外の9個のカルマンフィルタKFの動作を停止する制御を行う。なお、カルマンフィルタKFの推定精度とは、カルマンフィルタKFがシステムの状態をどの程度正確且つ安定的に推定しているかについて表す指標であり、具体的には、状態ベクトルxkの示す値と状態ベクトルの真値との誤差の大きさを評価する値誤差評価値、及び、状態ベクトルxkの収束の程度を表す値収束評価値により表現される指標である。そして、カルマンフィルタKFの推定精度が高い状態とは、状態ベクトルxkの示す値と状態ベクトルの真値との誤差が小さくつまり、誤差評価値が小さく、且つ、状態ベクトルxkが収束しているつまり、収束評価値が小さい状態のことを意味する。ここで、真値とは、状態ベクトルxkの各要素に対応する実際の物理量を表す値である。また、状態ベクトルの真値とは、状態ベクトルxkの各要素に対応する実際の物理量を表す値真値を要素とするベクトルである。状態ベクトルxkの示す値と状態ベクトルの真値との誤差が小さい場合つまり、状態ベクトルxkの各要素が真値に近い正確な値を示す場合、カルマンフィルタKFは、時刻T=kにおいて、システムを正確に推定していると看做すことができる。また、状態ベクトルxkが収束している場合、時刻T=kにおいて状態ベクトルxkの各要素が真値に近い正確な値を示していれば、時刻T=k以降においても状態ベクトルxの各要素が真値に近い正確な値を示す可能性が高い。つまり、状態ベクトルxkが収束している場合、カルマンフィルタKFは、システムの状態を安定的に推定していると看做すことができる。以下、カルマンフィルタ制御部400の詳細について説明する。カルマンフィルタ制御部400は、状態ベクトルxkの示す値と状態ベクトルの真値との誤差誤差評価値が最小のカルマンフィルタを選択する観測残差比較部420、状態ベクトルxkが収束していると看做すことができるか否かを判定する収束度判定部440、及び、カルマンフィルタKFの動作を停止させる制御を行う停止制御部460を備える。観測残差比較部420は、10個のピークモニタPM[1]~PM[10]と、比較部422とを備える。第m番目のピークモニタPM[m]は、以下の式5に示すように、第m番目のカルマンフィルタKF[m]が出力する観測残差zk[m]を、時刻T=k-τ0から時刻T=kまでの所定期間モニターし、当該期間における観測残差zk[m]のノルムの最大値を、誤差評価値MAXzk[m]として出力する。比較部422は、10個のピークモニタPM[1]~PM[10]がそれぞれ出力する誤差評価値MAXzk[1]~MAXzk[10]のうち、最小となる誤差評価値MAXzk[Sel]に対応するカルマンフィルタKF[Sel]を示す選択値Selを出力する但し、Selは、1≦Sel≦10を満たす自然数。なお、観測残差比較部420は、カルマンフィルタ部300が観測残差zk[1]~zk[10]の全てを出力する場合のみ、処理を実行する。例えば、CPU10は、カルマンフィルタ部300が、観測残差zk[1]~zk[10]のうち1つの観測残差zk[m]のみを出力する場合には、観測残差比較部420の動作を停止させる。前述のとおり、観測残差zkは、状態ベクトルxと観測値ベクトルyとを用いて算出される値であり、観測残差zkに基づいて、状態ベクトルxの示す値と状態ベクトルの真値との誤差の大きさを表すことができる。しかし、観測値ベクトルyには、ノイズが重畳し、時刻Tの経過に伴い変動する値であるため、観測残差zkも時刻Tの経過に伴い変動する。そこで、観測残差比較部420は、観測残差zk[m]を所定期間モニターし、所定期間における観測残差zk[m]のノルムの最大値を誤差評価値MAXzk[m]として採用し、誤差評価値MAXzk[m]を用いて、状態ベクトルxkの示す値と状態ベクトルの真値との誤差の大きさを評価する。これにより、観測残差比較部420は、状態ベクトルxkの示す値と真値との誤差の大きさを評価する際の、各種センサの測定誤差やノイズ等の影響を低減することができ、状態ベクトルxkの示す値と状態ベクトルの真値との誤差が最小となるカルマンフィルタKFを、正しく選択することが可能となる。なお、本実施形態では、一定期間における観測残差zk[m]のノルムの最大値を誤差評価値MAXzk[m]として採用するが、本発明はこのような形態に限定されるものではなく、観測残差zk[m]の一部の要素のノルムの所定期間における最大値を誤差評価値として採用してもよい。例えば、観測残差比較部420は、観測残差zk[m]を構成する要素のうち、観測値ベクトルykの磁気データqkに対応する要素について一定期間モニターし、当該要素のノルムの最大値を誤差評価値として採用するものであってもよい。また、本実施形態では、状態ベクトルxkの示す値と状態ベクトルの真値との誤差の大きさを評価する値として、観測残差zk[m]に基づいて定められる誤差評価値MAXzk[m]を用いるが、本発明はこれに限定するものではなく、誤差評価値MAXzk[m]以外であっても、状態ベクトルxkの示す値と状態ベクトルの真値との誤差の大きさを評価することのできる値であれば、いかなる値を用いてもよい。収束度判定部440は、10個のカルマンフィルタKF[1]~KF[10]のそれぞれが出力する状態ベクトルの推定誤差の共分散Pk[1]~Pk[10]のうち、選択値Selが示すカルマンフィルタKF[Sel]が出力する状態ベクトルの推定誤差の共分散Pk[Sel]のトレースTrPk[Sel]を、収束評価値として採用する。そして、収束度判定部440は、収束評価値TrPk[Sel]が、閾値PTH以下であるか否かを判定し、判定結果を出力する。一般的に、共分散行列の固有値は、当該固有値に対応する固有ベクトル方向の分散の大きさを表す。また、ある行列のトレース対角成分は、当該ある行列の有する複数の固有値の和を表す。つまり、共分散行列のトレースは、当該共分散行列の全ての固有ベクトルの方向の分散の大きさの和を表す値である。従って、収束評価値TrPk[Sel]によって、状態ベクトルxkの推定誤差の大きさ、つまり状態ベクトルxkの収束の程度を表すことができる。一方、状態ベクトルxkは、初期状態T=0から時刻T=k-1に至る全ての状態ベクトルx0~xk-1を考慮して、更新される値であるため、状態ベクトルの推定誤差の共分散Pkも、初期状態T=0から時刻T=k-1に至る全ての状態ベクトルの推定誤差の共分散P0~Pk-1を考慮した値となる。よって、収束評価値TrPk[Sel]は、初期状態T=0から時刻T=k-1に至るまでの期間における、状態ベクトルx0~xk-1の推定誤差を考慮した値となる。つまり、収束評価値TrPk[Sel]が小さな値となる場合、状態ベクトルxkの推定誤差が小さな値で安定していると看做すことができる。例えば、時刻T=kにおいて、状態ベクトルxkの各要素が真値に近い値を示していれば、時刻T=k以降においても、状態ベクトルxの各要素が真値に近い値を示す可能性が高い。このような収束評価値TrPk[Sel]を用いることにより、カルマンフィルタKF[Sel]の有する状態ベクトルxkの収束の程度を評価することができる。そして、収束評価値TrPk[Sel]が、閾値PTH以下である場合すなわち、状態ベクトルxkが収束していると看做すことができる場合、カルマンフィルタKF[Sel]は、システムの状態の正確な推定を安定的に行っていると看做すことができる。なお、本実施形態では、状態ベクトルxkの収束の程度を表す収束評価値として、状態ベクトルの推定誤差の共分散Pk[Sel]のトレースTrPk[Sel]を採用するが、本発明はこれに限定するものではなく、収束評価値TrPk[Sel]以外であっても、状態ベクトルxkの収束の程度を表す値であれば、いかなる値を用いてもよい。例えば、収束度判定部440は、状態ベクトルの推定誤差の共分散Pk[Sel]の最大固有値を収束評価値として採用してもよい。すなわち、収束度判定部440は、状態ベクトルの推定誤差の共分散Pk[Sel]の最大固有値が、予め定めた閾値以下であるか否かを判定するものであってもよい。この場合、状態ベクトルの推定誤差の共分散Pk[Sel]の最大固有値が、閾値以下であれば、状態ベクトルxkの推定誤差が小さな値で安定しているつまり、状態ベクトルxkが収束していると看做すことができる。なお、収束度判定部440は、カルマンフィルタ部300が共分散Pk[1]~Pk[10]の全てを出力する場合のみ、処理を実行する。例えば、CPU10は、カルマンフィルタ部300が共分散Pk[1]~Pk[10]のうち1個の共分散Pk[m]のみを出力する場合、収束度判定部440の動作を停止させる。停止制御部460は、収束度判定部440が出力する判定結果が閾値PTH以下である場合、比較部422が出力する選択値Selが示すカルマンフィルタKF[Sel]以外の9個のカルマンフィルタKFの動作を停止させる制御信号Ctrを出力する。出力情報生成部500は、カルマンフィルタ制御部400が選択したカルマンフィルタKF[Sel]から出力される状態ベクトルxk[Sel]に基づいて、携帯機器1の姿勢や、携帯機器1から見た地磁気の向きなどの出力情報を生成する。なお、出力情報生成部500は、カルマンフィルタ部300から、複数の状態ベクトルxk[1]~xk[10]が出力されている場合には、出力情報の生成は行わず、カルマンフィルタ部300から、単一の状態ベクトルxk[Sel]が出力される場合のみ、出力情報の生成を行う。前述のとおり、非線形カルマンフィルタの演算は、時刻T=kにおける状態ベクトルxkを、観測残差zkと時刻T=k-1における状態ベクトルxk-1とに基づいて算出する。すなわち、状態ベクトルxkは、初期状態時刻T=0から時刻T=k-1に至る全ての状態ベクトルx0~xk-1を考慮して、更新される。従って、初期ベクトルINIが、時刻T=0における状態ベクトルの真値から大きく乖離した値に設定された場合、カルマンフィルタ演算処理を行っても、状態ベクトルxkの示す値が状態ベクトルの真値に近い値に収束するまでに長い時間を要し、更には状態ベクトルxkの示す値が状態ベクトルの真値に近づかないこともある。これに対して、本実施形態に係る状態推定装置は、10個のカルマンフィルタKF[1]~KF[10]に、互いに異なる10個の初期ベクトルINI[1]~INI[10]を供給し、これら10個のカルマンフィルタKF[1]~KF[10]を並列で動作させる。そして、状態推定装置は、推定精度が最も高いカルマンフィルタKF[Sel]を特定し、特定された1個のカルマンフィルタKF[Sel]から出力される状態ベクトルxk[Sel]に基づいて、出力情報を生成する。このように、本実施形態に係る状態推定装置は、10個の初期ベクトルINI[1]~INI[10]のうち、時刻T=0における状態ベクトルの真値に近い値に設定された初期ベクトルINI[Sel]に基づいて動作するカルマンフィルタKF[Sel]により状態推定を行うため、状態ベクトルxkの正確且つ高速な収束が可能となった。すなわち、本実施形態に係る状態推定装置によれば、正確且つ高速な状態推定を行うことが可能となる。[2.状態推定装置の処理フロー]図4は、状態推定装置の動作を説明するためのフローチャートである。このフローチャートに示される処理は、CPU10が、状態推定プログラム100を実行することにより実施される。ステップS101において、CPU10は、初期ベクトル生成処理を実行する。具体的には、CPU10は、複数の初期ベクトルINI[1]~INI[10]を生成する。すなわち、CPU10は、ステップS101の初期ベクトル生成処理を実行することにより、初期値生成部200として機能する。ステップS102において、CPU10は、カルマンフィルタ並列演算処理を実行する。具体的には、CPU10は、複数のカルマンフィルタKF[1]~KF[10]の全てを用いて非線形カルマンフィルタの演算を実行することで、観測残差zk[1]~zk[10]と状態ベクトルの推定誤差の共分散Pk[1]~Pk[10]とを生成する。なお、カルマンフィルタ並列演算処理において、CPU10は、第1のデータ型を用いて第1演算モードにより非線形カルマンフィルタの演算を行う。カルマンフィルタ並列演算処理は、10個のカルマンフィルタKF[1]~KF[10]を並列に動作させる処理であるため、処理負荷が大きい。しかし、本実施形態に係る状態推定装置は、データサイズの最も小さな第1のデータ型を用いた第1演算モードによりカルマンフィルタ並列演算処理を行うため、処理負荷を大きく増加させることなく、並列処理を実行することができる。ステップS103において、CPU10は、観測残差比較処理を実行する。具体的には、CPU10は、ステップS102のカルマンフィルタ並列演算処理において生成された観測残差zk[1]~zk[10]に基づいて、誤差評価値MAXzk[1]~MAXzk[10]を算出することにより、状態ベクトルxkの示す値と状態ベクトルの真値との誤差が最も小さなカルマンフィルタKF[Sel]を選択する。すなわち、CPU10は、ステップS103の観測残差比較処理を実行することにより、観測残差比較部420として機能する。ステップS104において、CPU10は、収束度判定処理を実行する。具体的には、CPU10は、ステップS102のカルマンフィルタ並列演算処理において生成された状態ベクトルの推定誤差の共分散Pk[Sel]に基づいて、収束評価値TrPk[Sel]を算出し、収束評価値TrPk[Sel]が閾値PTH以下であるか否かを判定する。判定結果が肯定である場合つまり、収束評価値TrPk[Sel]が閾値PTH以下である場合、CPU10は処理をステップS105に進める。一方、判定結果が否定である場合、CPU10は処理をステップS102に戻す。すなわち、CPU10は、ステップS104の収束度判定処理を実行することにより、収束度判定部440として機能する。ステップS105において、CPU10は、カルマンフィルタ停止処理を実行する。具体的には、CPU10は、カルマンフィルタKF[Sel]以外の9個のカルマンフィルタの動作を停止させる。すなわち、CPU10は、ステップS105のカルマンフィルタ停止処理を実行することにより、停止制御部460として機能する。なお、CPU10は、観測残差比較処理ステップS103及び収束度判定処理ステップS104を実行する間も、カルマンフィルタ並列演算処理を実施する。つまり、ステップS105のカルマンフィルタ停止処理において、CPU10が9個のカルマンフィルタの動作を停止させるまでは、CPU10は、10個のカルマンフィルタKF[1]~KF[10]の全てを用いた非線形カルマンフィルタの並列演算を継続する。ステップS106において、CPU10は、カルマンフィルタ単独演算処理を実行する。具体的には、CPU10は、1個のカルマンフィルタKF[Sel]のみを用いて非線形カルマンフィルタの演算を行う。カルマンフィルタ単独演算処理において、CPU10は、第2のデータ型を用いる第2演算モードにより、カルマンフィルタKF[Sel]を動作させることで、非線形カルマンフィルタの演算を行う。また、CPU10は、ステップS106のカルマンフィルタ単独演算処理において、カルマンフィルタKF[Sel]が周期的に更新する状態ベクトルxkを逐次出力する。すなわち、CPU10は、ステップS102のカルマンフィルタ並列演算処理と、ステップS106のカルマンフィルタ単独演算処理とを実行することにより、カルマンフィルタ部300として機能する。ステップS107において、CPU10は、出力情報生成処理を実行する。具体的には、CPU10は、ステップS106のカルマンフィルタ単独演算処理において出力された状態ベクトルxkに基づいて、出力情報を生成する。本実施形態に係る状態推定装置の出力情報は、携帯機器1の姿勢μである。CPU10は、ステップS107の出力情報生成処理を実行することにより、出力情報生成部500として機能する。なお、本実施形態では、ステップS107の出力情報生成処理において、CPU10は、姿勢μを出力情報として生成するが、本発明はこれに限定されるものではなく、姿勢μ以外の値例えば、3次元磁気センサ70のオフセット推定値cOを出力情報として生成してもよい。また、本実施形態では、CPU10は、状態ベクトルxkのみに基づいて出力情報を生成するが、状態ベクトルxk以外の情報例えば、観測値ベクトルykを用いて出力情報を生成しても構わない。例えば、CPU10は、磁気センサのオフセット推定値cO及び磁気データqに基づいて、地上に固定された座標系における地磁気の向きを算出してもよい。ステップS106のカルマンフィルタ単独演算処理は、1個のカルマンフィルタKF[Sel]のみを動作させる処理であるため、処理負荷は、ステップS102のカルマンフィルタ並列演算処理に比べて小さい。一方、ステップS107の出力情報生成処理において生成される出力情報は、カルマンフィルタ単独演算処理において出力される状態ベクトルxkに基づいて生成されるため、カルマンフィルタ単独演算処理の演算の精度は高いことが好ましい。本実施形態では、CPU10は、第2のデータ型を用いる第2演算モードによりカルマンフィルタ単独演算処理を行うため、精度の高い状態ベクトルxkを算出することができる。これにより、状態推定装置は、正確な出力情報を算出することが可能となる。次に、CPU10は、状態推定プログラム100を終了する終了条件を充足するか否かを判定するステップS108。終了条件は、携帯機器1の仕様等に基づいて適宜定めればよい。例えば、携帯機器1の電源がオフ状態になることを終了条件としてもよい。終了条件を充足する場合、CPU10は、当該フローチャートに示す状態推定処理を終了する。一方、終了条件を充足しない場合、CPU10は、処理をステップS106に進める。[3.初期ベクトルの生成]前述のとおり、初期値生成部200は、初期ベクトルINI[1]~INI[10]を生成する。この初期ベクトルINI[1]~INI[10]は、時刻T=0における状態ベクトルの真値に近い値を示すものが少なくとも1つ含まれるように生成すれば、どのような方法で生成しても構わない。以下において、初期ベクトルINI[1]~INI[10]の各要素の生成方法の一例を示す。まず姿勢μの初期値μ0[m]は、例えば、以下の式6のように設定することが可能である。ここで、式6に示した姿勢の初期値μ0[m]の意味を、図5乃至図7を用いて説明する。なお、図5乃至図7は、理解を容易にするために、初期姿勢を、ティーポットの図形を用いて例示する。まず、図5に示すように、初期値μ0[1]を定める。初期値μ0[1]は、どのような値に設定してもよいが、本実施形態では、式6に示すように基本姿勢つまり、携帯機器1に固定された座標系の各軸と、地上に固定された座標系の各軸とが一致する姿勢を、初期値μ0[1]とする。次に、図5に示すように、地上に固定された座標系におけるz軸を回転軸として、初期値μ0[1]をそれぞれ90度、180度、270度回転させて3つの姿勢を生成し、これら3つの姿勢を、それぞれ、初期値μ0[2]、μ0[3]、μ0[4]とする。また、図6に示すように、地上に固定された座標系におけるy軸を回転軸として、初期値μ0[1]をそれぞれ90度、180度、270度回転させて3つの姿勢を生成し、これら3つの姿勢を、それぞれ、初期値μ0[5]、μ0[6]、μ0[7]とする。同様に、図7に示すように、地上に固定された座標系におけるx軸を回転軸として、初期値μ0[1]をそれぞれ90度、180度、270度回転させて3つの姿勢を生成し、これら3つの姿勢を、それぞれ、初期値μ0[8]、μ0[9]、μ0[10]とする。初期値生成部200は、以上に示した互いに異なる10個の初期値μ0[1]~μ0[10]を生成し、それぞれ、初期ベクトルINI[1]~INI[10]の要素とする。仮に、姿勢μの初期値μ0を、事前に定められた1つの固定値例えば、μ0=[0,0,0,1]に定める場合、初期値μ0は真の姿勢真値から最大で180度ずれることになる。これに対して、本実施形態に係る初期値生成部200は、基本姿勢を90度ずつ回転させることで、10個の初期姿勢を定める。この場合、10個の初期姿勢のうち、真の姿勢に最も近い初期姿勢と、真の姿勢とのなす角度は、90度以下となる。従って、本実施形態に係る初期値生成部200は、携帯機器1がどのような姿勢であっても、真値に近い値となる初期値μ0[m]を要素とする初期ベクトルINI[m]を生成することができる。次に、3次元磁気センサ70のオフセット推定値cOの初期値cO0について、生成方法の一例を述べる。前述のとおり、磁気センサのオフセットは、携帯機器1に搭載された部品が発する内部磁界の方向及び大きさを、携帯機器1に固定された座標系において表現したベクトルである。内部磁界は、3次元磁気センサ70の検出対象である地磁気に比べて大きな磁界となる場合があり、内部磁界を無視した場合には、正確な地磁気の向きを算出することはできない。従って、地磁気の向きを正確に算出するためには、オフセットを正確に把握し、3次元磁気センサ70が出力する磁気データqから、オフセットを減算することが必要となる。しかし、3次元磁気センサ70のオフセット推定値cOの初期値cO0に対して、事前に定められた固定値例えば、原点[0,0,0]T等を適用する場合、磁気センサのオフセット推定値cOの初期値cO0が、磁気センサの真のオフセット真値から大きく乖離した値となる可能性が高い。本実施形態に係る初期値生成部200は、3次元磁気センサ70のオフセット推定値cOの初期値cO0[m]を、姿勢μの初期値μ0[m]に基づいて算出する。具体的には、本実施形態に係る初期値生成部200は、3次元磁気センサ70のオフセット推定値cOの初期値cO0[m]を、時刻T=0に3次元磁気センサ70が観測する磁気データq0、姿勢の初期値μ0[m]、式9に示すベクトルBg、及び、式8に示す行列Bμを用いて、以下に示す式7で与えられる値に設定する。ここで、行列Bμは、携帯機器1が姿勢μである場合に、地上に固定された座標系において表現されたベクトルを、携帯機器1に固定された座標系において表現するための座標変換を行う行列である。また、式9に示すベクトルBgは、地上に固定された座標系において携帯機器1を使用する地域の地磁気の向き及び大きさを表すベクトルである。なお、本実施形態では、地上に固定された座標系は、y軸が磁極北を向き、z軸が鉛直上を向くように定められることとする。なお、ROM30には、位置情報と地磁気ベクトルBgとを対応づけて記憶したルックアップテーブルLUT2が記憶されている。初期値生成部200は、GPS部60で生成される位置情報に基づいてルックアップテーブルLUT2を参照して地磁気ベクトルBgを取得し、式7を演算することによって、磁気センサのオフセット推定値cOの初期値cO0[m]を得る。前述のとおり、10個の初期値μ0[1]~μ0[10]の中には、真の姿勢に近い値となるものが含まれるため、10個の初期値cO0[1]~cO0[10]の中には、磁気センサの真のオフセット真値に近い値となるものが含まれる。従って、本実施形態に係る初期値生成部200は、真のオフセットに近い値を有する3次元磁気センサ70のオフセット推定値cOの初期値cO0[m]を要素とする初期ベクトルINI[m]を生成することができる。地磁気の強さrの初期値r0[1]~r0[10]、及び地磁気の伏角φの初期値φ0[1]~φ0[10]は、例えば、GPS部60から供給される位置情報に基づいて生成することができる。これは、地球上の位置が特定できれば、その位置における地磁気の強さ及び伏角を知ることができるからである。具体的には、ROM30には、位置情報と地磁気の強さ及び伏角とを対応づけて記憶したルックアップテーブルLUT1が格納される。初期値生成部200は、ルックアップテーブルLUT1を参照して、位置情報に対応する地磁気の強さ及び伏角を、地磁気の強さrの初期値r0[1]~r0[10]及び地磁気の伏角φの初期値φ0[1]~φ0[10]として設定する。なお、携帯機器1が衛星からの電波の届かない場所例えば、地下街にある場合には、GPS部60から位置情報が得られない。そのような場合には、携帯機器1が使われ得る地域の典型的な値を採用すればよい。その値もルックアップテーブルLUT1に格納されている。初期値生成部200は、位置情報の取得が不能な場合には、典型的な値をルックアップテーブルLUT1から読み出す。例えば、日本で販売された携帯機器1については、明石市の地磁気の強さ及び伏角に基づいて、地磁気の強さrの初期値r0[1]~r0[10]、及び地磁気の伏角φの初期値φ0[1]~φ0[10]を算出する。角速度ωの初期値ω0[1]~ω0[10]は、例えば、携帯機器1が静止していることを仮定して、0に設定する。また、角速度センサのオフセット推定値bOの初期値bO0[1]~bO0[10]は、通常0近辺に調整されているため、0に設定する。このように、本実施形態に係る状態推定装置は、真値に近い値を示すものが少なくとも1個は含まれるように、互いに異なる10個の初期ベクトルINI[1]~INI[10]を生成し、10個のカルマンフィルタKF[1]~KF[10]を並列動作させる。従って、本実施形態に係る状態推定装置は、時刻T=0における状態ベクトルの真値に近い値を示す初期ベクトルINI[m]が供給されて動作するカルマンフィルタKF[m]により、状態ベクトルxkを正確且つ高速に収束させることが可能となる。[4.非線形カルマンフィルタによる演算]次に、非線形カルマンフィルタの演算について説明する。一般的に、カルマンフィルタは、動的システムの状態の経時的な変化を表す状態遷移モデルを用いて、ある時刻時刻T=k-1のシステムの状態を示す状態ベクトルxk-1から、単位時間が経過した後時刻T=kの状態ベクトルxkを推定する。この推定値を、推定状態ベクトルx-kと称する。そして、カルマンフィルタは、各種センサからの出力を要素とする観測値ベクトルykと状態ベクトルxkとの関係を表す観測モデルを用いて、推定状態ベクトルx-kから、観測値ベクトルykを推定する。この推定値を、推定観測値ベクトルy-kと称する。次に、カルマンフィルタは、推定観測値ベクトルy-kと、実際の観測値を要素とする観測値ベクトルykとの差分を観測残差zkとして算出するとともに、観測残差zkに基づいてカルマンゲインKkを算出する。その後、カルマンフィルタは、観測残差zk、カルマンゲインKk、及び、推定状態ベクトルx-kを用いて、状態ベクトルxkを更新し、更新後の状態ベクトルxkを算出する。カルマンフィルタは、以上のような状態ベクトルxkの更新を繰り返す演算を行うことにより、状態ベクトルxkの各要素を、実際の物理量を正確に表した値真値に近い正確な値へと近付ける。なお、詳細は後述するが、本実施形態ではカルマンフィルタとして、非線形カルマンフィルタであるシグマポイントカルマンフィルタを用いる。シグマポイントカルマンフィルタとは、状態ベクトルxk-1の周囲に複数のシグマポイントχk-1を生成し、これらの複数のシグマポイントχk-1を状態遷移モデルに適用することで、単位時間経過後のシグマポイントχ-kを推定し、推定された複数のシグマポイントχ-kの平均を算出することにより、推定状態ベクトルx-kを求める演算である。従って、推定観測値ベクトルy-kは、厳密には、推定状態ベクトルx-kの周辺に存在する複数のシグマポイントχ-kに基づいて算出される。本実施形態における、非線形カルマンフィルタの状態遷移モデルは以下に示す式10で与えられ、観測モデルは以下に示す式11で与えられる。なお、本実施形態では、式10に現れる関数f及び式11に現れる関数hは、非線形の関数である。ここで、状態ベクトルxkの次元をdimxと表し、観測値ベクトルykの次元をdimyと表す。本実施形態では、dimx=15であり、dimy=9である。また、式10に現れるプロセスノイズwk及び式11に現れる観測ノイズvkは、0を中心とするガウスノイズである。式10は、時刻T=kにおける状態ベクトルxkが、時刻T=k-1における状態ベクトルxk-1を関数fに代入して得られた値と、時刻T=k-1におけるプロセスノイズwk-1とを加算することにより推定されることを示している。また、式11は、時刻T=kにおける観測値ベクトルykが、時刻T=kにおける状態ベクトルxkを関数hに代入して得られる値と、時刻T=kにおける観測ノイズvkとを加算することにより推定されることを示している。なお、プロセスノイズwkの共分散をQk、観測ノイズvkの共分散をRkと表す。時刻T=kにおける観測残差zkは、観測値ベクトルykと、推定観測値ベクトルy-kとに基づいて定められるベクトルであり、以下に示す式12で表される。カルマンフィルタは、以下の式14に示すように、観測残差zk、推定状態ベクトルx-k、及び、式13に示すカルマンゲインKkを用いて、更新後の状態ベクトルxkを算出する。また、カルマンフィルタは、以下の式15に示すように、状態ベクトルxkの推定誤差の共分散Pkを更新する。ここで、P-kは、推定状態ベクトルx-kの推定誤差の共分散であり、Pkは、更新後の状態ベクトルxkの推定誤差の共分散である。また、Pyykは、観測残差zkの共分散であり、Pxykは、推定状態ベクトルx-kと、推定観測値ベクトルy-kとの相互共分散行列である。なお、推定観測値ベクトルy-kは、式11に示すように、状態遷移モデルを用いて算出された推定状態ベクトルx-kを、観測モデルに適用することで算出される値である。よって、推定観測値ベクトルy-kと実際のセンサ出力に基づく観測値ベクトルykとの差分である観測残差zkに基づいて、推定状態ベクトルx-kと実際の物理量を正確に表した値状態ベクトルの真値との誤差の大きさを評価することが可能である。式14に示すように、観測残差zkを用いて、推定状態ベクトルx-kを、更新後の状態ベクトルxkへと更新することにより、状態ベクトルxkの各要素を真値に近い値へと近付けることができる。本実施形態において、10個のカルマンフィルタKF[1]~KF[10]は、アンセンテッド変換を用いたシグマポイントカルマンフィルタの演算を実行する。以下では、シグマポイントカルマンフィルタの演算について、具体的に説明する。図8は、状態推定装置がカルマンフィルタKF[m]を実行することで実現される機能を表す機能ブロック図である。以下、図8を参照しつつ、カルマンフィルタKF[m]の動作について説明する。なお、10個のカルマンフィルタKF[1]~KF[10]の構成は同一であり、初期値生成部200からそれぞれ供給される初期ベクトルINI[1]~INI[10]のみが相違する。遅延部301は、加算器311から出力される更新後の状態ベクトルxkを、単位時間時刻T=k-1から時刻T=kに相当する時間だけ遅延させることで、状態ベクトルxk-1を生成し、これを、シグマポイント生成部302に対して出力する。なお、初回の演算時刻T=0では、遅延部301は、初期値生成部200が出力する初期値INI[m]を、状態ベクトルxk-1として採用する。また、遅延部301は、減算器312から出力される更新後の状態ベクトルxkの推定誤差の共分散Pkを単位時間遅延させることで、状態ベクトルxk-1の推定誤差の共分散Pk-1を生成し、これを、シグマポイント生成部302に対して出力する。シグマポイント生成部302は、dimx行dimx列の共分散行列Pk-1及びQk-1を用いて、2dimx1個のシグマポイントを生成する。具体的には、まず、複数のシグマポイントの平均に対して、平均のまわりのシグマポイントの広がりを表すスケーリングパラメータλを用いて、式16及び式17に示すベクトルσkを定義する。このとき、シグマポイント生成部302は、ベクトルσk-1と状態ベクトルxk-1とに基づいて、式18及び式19で示される2dimx1個のシグマポイントχk-1を生成する。状態遷移モデル部303は、式20に示すように、時刻T=k-1における2dimx1個のシグマポイントχk-1jの各々を、状態遷移モデルに適用することで、時刻T=kにおける2dimx1個のシグマポイントの推定値χ-kjを算出する。次に、平均算出部304は、式21に示すように、時刻T=kにおける2dimx1個のシグマポイントの推定値χ-kの平均値を計算することで、推定状態ベクトルx-kを算出する。また、平均算出部304は、式22に示す、推定状態ベクトルx-kの推定誤差の共分散P-kを算出する。一方、観測モデル部305は、式23に示すように、時刻T=kにおける2dimx1個のシグマポイントχkjの各々を、観測モデルに適用することにより、2dimx1個の推定観測値γkjを算出する。そして、平均処理部306は、式24に示すように、2dimx1個の推定観測値γkjの平均を演算することにより、推定観測値ベクトルy-kを算出する。次に、平均処理部306は、式25に示す観測残差の共分散Pyykを算出する。減算器307は、式12に示したように、観測値ベクトルykと、推定観測値ベクトルy-kとの差分として、観測残差zkを算出する。カルマンフィルタKF[m]は、この観測残差zkを出力する。カルマンゲイン付与部308は、式26に示す相互共分散行列Pxykを算出する。そして、カルマンゲイン付与部308は、式13に示したように、観測残差zkの共分散Pyykと、相互共分散行列Pxykとに基づいて、カルマンゲインKkを算出し、式14の右辺第2項Kkzkの演算を実行する。また、カルマンゲイン付与部308は、式15の右辺第2項KkPyykKTkの演算を実行する。加算器311は、式14に示したように、推定状態ベクトルx-kと、カルマンゲイン付与部308から出力される式14の右辺第2項Kkzkとを加算することにより、更新後の状態ベクトルxkを算出する。カルマンフィルタKF[m]は、更新後の状態ベクトルxkを出力する。減算器312は、式15に示したように、推定状態ベクトルx-kの推定誤差の共分散P-kと、カルマンゲイン付与部308から出力される式15の右辺第2項KkPyykKTkとの差分として、更新後の状態ベクトルxkの推定誤差の共分散Pkを算出するすなわち、状態ベクトルの推定誤差の共分散Pkを、P-kからPkに更新する。カルマンフィルタKF[m]は、更新後の状態ベクトルxkの推定誤差の共分散Pkを出力する。次に、状態遷移モデル部303が行う演算で用いられる状態遷移モデルについて説明する。本実施形態に係る状態遷移モデルにおいて、状態ベクトルxkを構成する状態変数のうち、姿勢μについての状態遷移は、式27により定義される。式27は、時刻T=k-1における姿勢μk-1から、単位時間経過後の時刻T=kにおける姿勢μkを推定する演算を表す。ここで、μ-kは、推定状態ベクトルx-kのうち、姿勢μを表す状態変数に対応する要素である。なお、式27の右辺の演算子Ωは、式28により定義される。ここで、I3×3は3行3列の単位行列を表す。3次元ベクトルl=l1,l2,l3に対して、演算子[l×]は、式29で定義される。また、測定時間間隔時刻T=k-1から時刻T=kまでの間隔をΔtで表し、時刻T=kにおける更新後の状態ベクトルxkのうち角速度を表す状態変数に対応する要素をωkで表したとき、演算子Ωを構成する成分Ψkは、式30で定義される。式2に示すように、姿勢μは、クォータニオンで表現され、正規化条件μ=1が満たされる必要がある。しかし、シグマポイントカルマンフィルタを用いて状態ベクトルxkを更新する場合、更新後の状態ベクトルxkは、式21に示すように、シグマポイントの推定値χ-kiの平均として算出されるため、更新前の状態ベクトルxk-1のノルムと、更新後の状態ベクトルxkのノルムとは、異なる値となることがある。従って、シグマポイントカルマンフィルタの演算を行う場合、更新前の状態ベクトルxk-1の要素である姿勢μk-1のノルムと、更新後の状態ベクトルxkの要素である姿勢μkのノルムとは、異なる値となる可能性が存在する。つまり、シグマポイントカルマンフィルタの演算を行う場合、姿勢μについての正規化条件が満たされなくなる可能性が存在する。そこで、姿勢μに対して何らかの演算が行われるときには、演算後の結果をそのベクトル自身の大きさで正規化する。なお、より厳密に正規化条件を保つためには、状態ベクトルxkを構成する要素のうち姿勢μkについてはMRPsを用いて前時刻との差分情報だけに限定し、カルマンフィルタの外部にある姿勢情報をカルマンフィルタから得られる差分情報に基づいて更新すればよい。地磁気の強さr及び地磁気の伏角φは、変化を予測することが難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける地磁気の強さrk及び伏角φkと、時刻T=k-1における地磁気の強さrk-1及び伏角φk-1とは等しい値であると仮定する。同様に、角速度センサのオフセット推定値bOは、変化を予測することが難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける角速度センサのオフセット推定値bOKと、時刻T=k-1における角速度センサのオフセット推定値bOK-1とは等しい値であると仮定する。携帯機器1の角速度ωは、携帯機器1の利用者による携帯機器1の動かし方に依存して変化するため、時刻T=k-1の角速度ωk-1を用いて、時刻T=kにおける角速度ωkを定式化することは難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける角速度ωkと、時刻T=k-1における角速度ωk-1とは等しいと仮定する。前述のとおり、磁気センサのオフセットは、携帯機器1の部品が発する内部磁界の方向及び大きさを携帯機器1に固定された座標系において表現したベクトルである。従って、携帯機器1の内部状態が一定である間は、磁気センサのオフセットも変化しない。一方、携帯機器1の内部状態が変化した場合、例えば、携帯機器1に搭載された部品を流れる電流の大きさが変化した場合や、携帯機器1に搭載された部品の着磁状況が変化した場合には、磁気センサのオフセットも変化する。すなわち、携帯機器1の内部状態は、携帯機器1の利用者による携帯機器1の操作や、携帯機器1の外部の環境等に依存して変化する。従って、磁気センサのオフセットが変化するタイミングや磁気センサのオフセットの変化量を予測することは困難であり、時刻T=k-1における磁気センサのオフセット推定値cOK-1を用いて、時刻T=kにおける磁気センサのオフセット推定値cOKを定式化することは難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける磁気センサのオフセット推定値cOKと、時刻T=k-1における磁気センサのオフセット推定値cOK-1とは等しいと仮定する。このように、本実施形態に係る状態遷移モデルでは、以下に示す式31のように、状態ベクトルxkを構成する複数の状態変数のうち、姿勢μを表す状態変数以外は、前の時刻から変化しないという前提を置く。なお、カルマンフィルタの演算においては、状態ベクトルxkの各要素は、観測残差zkに基づいて真値に近づくように更新される。従って、カルマンフィルタの演算において、式31に示された姿勢μ以外の値が全く変化しないわけではない。次に、観測モデル部305が行う演算で用いられる観測モデルについて説明する。3次元角速度センサ90から出力される角速度データgの推定値γgyroは、角速度ωと、角速度センサのオフセット推定値bOとを用いて、式32で与えられる。また、地上に固定された座標系において地磁気を表すベクトルBgは式33で与えられる。従って、3次元磁気センサ70から出力される磁気データqの推定値γmagは、磁気センサのオフセット推定値cOと、式8で示した行列Bμとを用いて、式34で与えられる。また、地上に固定された座標系において、重力加速度を表したベクトルを、重力加速度の大きさで正規化した3次元のベクトルGRVは、以下の式35で表される。従って、3次元加速度センサ80から出力される加速度データaの推定値γaccは、式8で示した行列BμとベクトルGRVとを用いて、式36で与えられる。式3を式12の右辺第1項に代入し、式32、式34、及び式36を用いて式12の右辺第2項を表すことにより、式12を、以下の式37に変形することができる。このとき、観測残差zkは、式37により算出される。[5.結論]以上に示したように、本実施形態に係る状態推定装置は、10個のカルマンフィルタKF[1]~KF[10]に、互いに異なる10個の初期ベクトルINI[1]~INI[10]を供給し、これら10個のカルマンフィルタKF[1]~KF[10]を並列で動作させる。そして、状態推定装置は、推定精度が最も高いカルマンフィルタKF[Sel]を特定し、特定されたカルマンフィルタKF[Sel]により状態推定を行う。従って、本実施形態に係る状態推定装置は、状態ベクトルxkを正確且つ高速に収束させることが可能であり、正確且つ高速な状態推定が可能となった。また、本実施形態に係る状態推定装置は、カルマンフィルタKFの推定精度を、誤差評価値MAXzk[m]及び収束評価値TrPk[Sel]を用いて表現した。誤差評価値MAXzk[m]は、状態ベクトルxkの示す値と状態ベクトルの真値との誤差の大きさを評価する。従って、本実施形態に係る状態推定装置は、真値に近い正確な値を要素とする状態ベクトルxkを有し、真値に短時間で状態ベクトルxkを収束させることができるカルマンフィルタKF[Sel]を選択することが可能となった。また、収束評価値TrPk[Sel]は、状態ベクトルxkの収束の程度を表す。従って、本実施形態に係る状態推定装置は、システムの状態の正確な推定を安定的に行うことのできるカルマンフィルタKF[Sel]を特定することが可能となった。このように、本実施形態に係る状態推定装置は、カルマンフィルタKFの推定精度に基づいて特定された1個のカルマンフィルタKF[Sel]を用いて、システムの状態を推定するため、正確且つ高速な状態推定を行うことが可能となった。また、本実施形態に係る状態推定装置は、第1演算モードによりカルマンフィルタ並列演算処理を行うとともに、第2演算モードによりカルマンフィルタ単独演算処理を行う。これにより、本実施形態に係る状態推定装置は、処理負荷を低減しつつ、正確且つ精度の高い状態推定を行うことが可能となった。<B.変形例>本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。1変形例1上述した実施形態では、カルマンフィルタKFは、シグマポイントカルマンフィルタの演算を実行したが、本発明はこれに限定するものではなく、カルマンフィルタKFは、拡張カルマンフィルタ等、シグマポイントカルマンフィルタ以外の非線形カルマンフィルタの演算を実行するものであってもよい。本発明に係る状態推定装置は、複数のカルマンフィルタKFの中で、推定精度が最も高いカルマンフィルタKFを特定し、当該特定されたカルマンフィルタKFを用いてシステムの状態を推定する。そして、カルマンフィルタKFが、どのような非線形カルマンフィルタの演算を実行する場合であっても、最も推定精度の高いカルマンフィルタKFの状態ベクトルは、正確な値に高速に収束する。すなわち、本発明に係る状態推定装置は、カルマンフィルタKFが、どのような非線形カルマンフィルタの演算を実行するものであっても、正確且つ高速な状態推定を行うことができる。2変形例2上述した実施形態及び変形例では、カルマンフィルタ部300は、10個のカルマンフィルタKF[1]~KF[10]を備えるが、本発明はこのような形態に限定されるものではなく、カルマンフィルタ部300は、2個以上のカルマンフィルタKFを備えるものであればよい。この場合、初期値生成部200は、カルマンフィルタKFの個数に対応して、互いに異なる複数の初期ベクトルINIを生成するものであればよい。なお、カルマンフィルタ部300が備えるカルマンフィルタKFの個数は、携帯機器1の処理性能を考慮して、適宜定めればよい。例えば、携帯機器1の処理性能が高い場合には、カルマンフィルタ部300は、10個よりも多くのカルマンフィルタKFを備えてもよい。この場合、状態ベクトルxを高速且つ正確に収束させることが可能となる。また、携帯機器1の処理性能が低い場合には、カルマンフィルタ部300は、10個よりも少ないカルマンフィルタを備えてもよい。この場合、カルマンフィルタ並列演算処理に係る処理負荷を軽減することができる。3変形例3上述した実施形態及び変形例では、状態推定装置は、第1演算モードによりカルマンフィルタ並列演算処理を行うとともに、第2演算モードによりカルマンフィルタ単独演算処理を行うが、本発明は、状態推定装置の動作をこのような態様に限定するものではない。例えば、状態推定装置は、カルマンフィルタ並列演算処理及びカルマンフィルタ単独演算処理の双方を、第1演算モードにより実行するものであってもよい。この場合、カルマンフィルタ単独演算処理に係る処理負荷を軽減することができる。また、状態推定装置は、カルマンフィルタ並列演算処理及びカルマンフィルタ単独演算処理の双方を、第2演算モードにより実行するものであってもよい。この場合、状態推定装置は、精度の高い推定を行うことが可能となる。4変形例4上述した実施形態及び変形例では、初期値生成部200は、初期ベクトルINIを構成する複数の状態変数のうち、姿勢μの初期値μ0及び3次元磁気センサ70のオフセット推定値cOの初期値cO0を互いに異なる値に設定することで、互いに異なる複数の初期ベクトルINIを生成したが、本発明はこのような初期ベクトルINIの生成方法に限定するものではない。例えば、初期値生成部200は、初期ベクトルINIを構成する状態変数の全てについて、互いに異なる値に設定するものであってもよい。5変形例5上述した実施形態及び変形例では、状態推定装置は出力情報生成部500を備えるが、本発明はこのような構成に限定されるものではなく、状態推定装置は、出力情報生成部500を備えずに構成されてもよい。この場合、状態推定装置は、状態ベクトルxkを出力するものであってもよい。また、この場合、状態推定プログラム100と連携する各種アプリケーションが、姿勢や地磁気の向き等の出力情報を生成してもよい。6変形例6上述した実施形態及び変形例では、携帯機器1は、3次元磁気センサ70、3次元加速度センサ80、及び、3次元角速度センサ90を備えたが、本発明はこのような形態に限定するものではなく、携帯機器1は、どのようなセンサを備えるものであってもよい。なお、観測値ベクトルyは、携帯機器1が備える複数のセンサからの出力値の全部または一部を要素とするように定めればよい。本発明は、異種の物理量を測定する複数のセンサの出力を統合し、動的システムの状態を推定する非線形カルマンフィルタの演算において、正確且つ高速な状態推定を行うものである。従って、この発明によれば、状態推定装置がどのようなセンサを備える場合であっても、正確且つ高速な状態推定を行うことができる。7変形例7上述した実施形態及び変形例では、状態ベクトルxを構成する状態変数として、式1に示した変数、すなわち、携帯機器1の姿勢μ、地磁気の強さr、地磁気の伏角φ、携帯機器1の角速度ω、角速度センサ91~93のオフセット推定値bO、及び磁気センサ71~73のオフセット推定値cOを採用し、これらを推定の対象としたが、本発明は、状態ベクトルxを式1に限定するものではない。例えば、状態変数として、式1に示した変数のうちの一部を採用するものであってもよい。また、状態変数として、式1に示した変数以外の変数を採用するものであってもよい。"