"磁界を発生させる部品を備えた機器に組み込まれる状態推定装置であって、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサと、加速度を直交する3軸で検出して、3軸の座標系において表現される3次元のベクトルデータとして加速度データを順次出力する3次元加速度センサと、前記3次元磁気センサから順次出力される磁気データを蓄積する蓄積手段と、前記蓄積手段に蓄積された複数の前記磁気データに基づいて、前記部品の発する磁界を表す前記3次元磁気センサのオフセットの近似値を示す中心点の3軸の座標を算出する中心点算出手段と、複数の前記加速度データに基づいて、前記機器の動きが安定しているか否かを判定する安定性判定手段と、前記安定性判定手段の判定結果が肯定である場合、前記中心点の3軸の座標、前記磁気データ、及び、前記加速度データに基づいて、前記オフセットの推定値及び前記機器の姿勢を推定するベクトルを要素として含む初期ベクトルを生成する初期ベクトル生成手段と、前記初期ベクトルを初期値とする状態ベクトルを、前記3次元磁気センサ及び前記3次元加速度センサを含む前記機器に組み込まれた複数のセンサの出力を要素とする観測値ベクトルを用いて更新するカルマンフィルタ演算手段と、を備えることを特徴とする、状態推定装置。前記中心点算出手段が前記中心点の算出に用いる複数の前記磁気データの取得を開始する前に、利用者に前記機器の姿勢を変化させるように促す第1の報知手段を、更に備えることを特徴とする、請求項1に記載の状態推定装置。前記中心点算出手段が前記中心点を算出した後に、利用者に前記機器の姿勢を変化させないように促す第2の報知手段を、更に備えることを特徴とする、請求項2に記載の状態推定装置。前記安定性判定手段は、複数の前記加速度データに基づいて基準加速度データを算出するローパスフィルタを備え、前記加速度データの示す座標と前記基準加速度データの示す座標との誤差を表す安定性指標を、所定個の前記加速度データの各々について算出し、前記安定性指標の全てが第1閾値以下である場合、前記機器の動きが安定していると判定することを特徴とする、請求項1乃至3のうちいずれか1項に記載の状態推定装置。前記安定性判定手段は、所定個の前記加速度データの示す3軸の座標の分散を、安定性指標として算出し、前記安定性指標が、第2閾値以下である場合、前記機器の動きが安定していると判定することを特徴とする、請求項1乃至3のうちいずれか1項に記載の状態推定装置。しかし、非線形カルマンフィルタは、状態ベクトルと、状態ベクトルに基づいて算出される観測残差とを用いて、状態ベクトルを逐次更新する演算であるため、状態ベクトルの初期値の要素が真値から乖離した値に設定された場合、状態ベクトルの各要素が真値に近い値に収束するまでに長時間を要することがあり、更には、状態ベクトルの各要素が真値とは異なる不正確な値に収束することも起こり得る。本発明は、上述した点に鑑みてなされたものであり、真値に近い値を要素とする状態ベクトルの初期値を算出して非線形カルマンフィルタの演算を行うことで、地磁気の向きや物体の姿勢を高速且つ正確に推定することを解決課題とする。上述した課題を解決するため、本発明に係る状態推定装置は、磁界を発生させる部品を備えた機器に組み込まれる状態推定装置であって、互いに直交する3方向の磁気成分をそれぞれ検出して、3軸の座標系において表現される3次元のベクトルデータとして磁気データを順次出力する3次元磁気センサと、加速度を直交する3軸で検出して、3軸の座標系において表現される3次元のベクトルデータとして加速度データを順次出力する3次元加速度センサと、前記3次元磁気センサから順次出力される磁気データを蓄積する蓄積手段と、前記蓄積手段に蓄積された複数の前記磁気データに基づいて、前記部品の発する磁界を表す前記3次元磁気センサのオフセットの近似値を示す中心点の3軸の座標を算出する中心点算出手段と、複数の前記加速度データに基づいて、前記機器の動きが安定しているか否かを判定する安定性判定手段と、前記安定性判定手段の判定結果が肯定である場合、前記中心点の3軸の座標、前記磁気データ、及び、前記加速度データに基づいて、前記オフセットの推定値及び前記機器の姿勢を推定するベクトルを要素として含む初期ベクトルを生成する初期ベクトル生成手段と、前記初期ベクトルを初期値とする状態ベクトルを、前記3次元磁気センサ及び前記3次元加速度センサを含む前記機器に組み込まれた複数のセンサの出力を要素とする観測値ベクトルを用いて更新するカルマンフィルタ演算手段と、を備えることを特徴とする。カルマンフィルタの演算において、状態ベクトルの初期値初期ベクトルの各要素が真値から乖離した値である場合、状態ベクトルの各要素が真値に近い値に収束するまでに長い時間を必要とする。この発明によれば、中心点算出手段が算出する中心点の座標に基づいて定められる3次元磁気センサのオフセット推定値と、機器の動きが安定していると判定されたときの磁気データ及び加速度データとに基づいて算出される姿勢の推定値とを、初期ベクトルの要素とする。中心点算出手段が算出する中心点は、3次元磁気センサのオフセットの近似値を示す座標なので、オフセット推定値を真値に近い値に設定することができる。そして、3次元磁気センサのオフセット推定値及び磁気データを用いて、正確な地磁気の向きを得ることができる。また、機器の動きが安定しているときの加速度データより重力加速度の向きを得ることができるため、重力加速度の向き及び地磁気の向きに基づいて、真値に近い姿勢の推定値を算出することができる。すなわち、この発明によれば、初期ベクトルの各要素を真値に近い値に設定することができ、カルマンフィルタの演算において状態ベクトルの各要素を高速に真値に近い値に収束させることが可能となる。また、上述した状態推定装置において、前記中心点算出手段が前記中心点の算出に用いる複数の前記磁気データの取得を開始する前に、利用者に前記機器の姿勢を変化させるように促す第1の報知手段を、更に備えることが好ましい。3次元磁気センサは、地磁気と、部品の発する磁界とを検出する。3次元磁気センサの姿勢を変化させる場合、地磁気は、大きさは一定であるが向きが変化する磁界として検出される。一方、部品の発する磁界は、3次元磁気センサの姿勢をどのように変化させても、一定の向き及び大ききを有する磁界として検出される。従って、3次元磁気センサの姿勢を3次元的に変化させつつ検出された複数の磁気データが、3次元磁気センサに固定された座標系において示す座標は、部品の発する磁界を表すベクトルの示す座標を中心点とし地磁気の大きさを半径とする球面近傍に分布する。そして、当該球面の中心点の座標が、3次元磁気センサのオフセットの近似値を表す。よって、複数の磁気データに基づいてオフセットの近似値を算出するためには、複数の磁気データの示す座標を近傍に有する球面を一意に特定することが必要であり、そのためには、複数の磁気データの示す座標が、2次元的に分布せず、3次元的な広がりを有して分布していることが必要である。この発明によれば、利用者に機器の姿勢を変化させるように促す第1の報知手段を備えるため、複数の磁気データを取得する際に利用者が機器の姿勢を3次元的に変化させることが期待される。機器の姿勢を3次元的に変化させつつ取得した複数の磁気データの示す座標は、3次元的な広がりを有して分布する。従って、複数の磁気データの示す座標を近傍に有する球面を一意に特定することが可能となり、3次元磁気センサのオフセットの近似値を示す中心点の座標を算出することが可能となる。また、上述した状態推定装置において、前記中心点算出手段が前記中心点を算出した後に、利用者に前記機器の姿勢を変化させないように促す第2の報知手段を、更に備えることが好ましい。この発明によれば、利用者に機器の姿勢を変化させないように促す第2の報知手段を備えるため、利用者が機器を動かさず安定的な状態に置くことが期待される。機器の動きが安定している場合、加速度データより重力加速度の向きを得ることができるため、重力加速度の向き及び地磁気の向きに基づいて、真値に近い姿勢の推定値を得ることが可能となる。また、上述した状態推定装置において、前記安定性判定手段は、複数の前記加速度データに基づいて基準加速度データを算出するローパスフィルタを備え、前記加速度データの示す座標と前記基準加速度データの示す座標との誤差を表す安定性指標を、所定個の前記加速度データの各々について算出し、前記安定性指標の全てが第1閾値以下である場合、前記機器の動きが安定していると判定することを特徴とすることが好ましい。この発明によれば、加速度データをローパスフィルタに通過させることで得られる基準加速度データと、所定個の加速度データとの誤差を表す所定個の安定性指標を算出することで、所定個の加速度データがほぼ一定の値を示すものであるか否かを判定する。複数の加速度データの示す値の変動が小さくほぼ一定の値を示すとき、例えば機器が等加速度運動をする等の特殊な場合を除き、加速度データは重力加速度を示すものと看做すことができる。従って、本発明によれば、加速度データから重力加速度の向きを算出することができるか否かを判定し、判定結果が肯定である場合には、加速度データから算出された重力加速度の向きと地磁気の向きとに基づいて、真値に近い姿勢の推定値を得ることが可能となる。また、上述した状態推定装置において、前記安定性判定手段は、所定個の前記加速度データの示す3軸の座標の分散を、安定性指標として算出し、前記安定性指標が、第2閾値以下である場合、前記機器の動きが安定していると判定することを特徴とすることが好ましい。この発明によれば、所定個の前記加速度データの示す座標の分散を表す安定性指標を算出し、所定個の加速度データがほぼ一定の値を示すものであるか否かを判定する。従って、本発明によれば、加速度データから重力加速度の向きを算出することができるか否かを判定し、判定結果が肯定である場合には、加速度データから算出された重力加速度の向きと地磁気の向きとに基づいて、真値に近い姿勢の推定値を得ることが可能となる。<A.実施形態>本発明の実施の形態について図面を参照して説明する。[1.機器構成及びソフトウェア構成]図1は、本発明の実施形態に係る携帯機器のブロック図であり、図2は携帯機器の外観を示す斜視図である。携帯機器1は、画面の姿勢に応じて地図などの画像を回転させることによって、画像によって表されている方位を現実空間の方位に追従させる機能を備える。この機能は、各種センサの出力に基づいてカルマンフィルタの演算を行い、携帯機器1の姿勢と、携帯機器1から見た地磁気Bgの向きとを推定することによって実現される。携帯機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM蓄積部20、状態推定プログラム100等の各種のプログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及びGPS部60を備える。また、携帯機器1は、地磁気Bg等の磁気を検出して磁気データを出力する3次元磁気センサ70、加速度を検出して加速度データを出力する3次元加速度センサ80、及び角速度を検出して角速度データを出力する3次元角速度センサ90を備える。表示部50は、CPU10が状態推定プログラム100を実行することにより推定したシステムの状態に基づいて算出された地磁気Bgの向きや携帯機器1の姿勢μ等を、矢印等の画像により表示する。GPS部60は、GPS衛星からの信号を受信して携帯機器1の位置情報緯度、経度を生成する。3次元磁気センサ70は、X軸磁気センサ71、Y軸磁気センサ72、及びZ軸磁気センサ73を備える。各センサは、MI素子磁気インピーダンス素子、MR素子磁気抵抗効果素子などを用いて構成することができる。磁気センサI/F74は、各センサからの出力信号をAD変換して磁気データqを出力する。この磁気データqは、x軸、y軸およびz軸の3成分によって表されるベクトルデータである。ところで、3次元磁気センサ70が搭載される携帯機器1は、着磁性を有する各種金属や、電気回路等、磁界を発生させる部品が備えられることが多い。このため、3次元磁気センサ70が出力するベクトルデータ磁気データqは、地磁気Bgを表すベクトルの他に、携帯機器1に搭載された部品が発する磁界内部磁界Biを表すベクトルも含む値となる。従って、地磁気Bgの向きを正確に知るためには、3次元磁気センサが出力するベクトルデータ磁気データqから、携帯機器1の部品が発する内部磁界Biを表すベクトルを取り除く補正処理が必要となる。補正処理を説明するための前提として、まず、地磁気Bg及び内部磁界Biの性質について説明する。図3Aに示すように、地磁気Bgは、磁極北に向かう水平成分と伏角方向の鉛直成分を有する磁界であり、地上に固定された座標系である地上座標系ΣGにおいて、一定の向き及び大きさを有するベクトルGBgとして表される。一方、内部磁界Biは、携帯機器1の部品が発する磁界であり、携帯機器1に対して一定の方向を向いた、一定の大きさの磁界である。つまり、内部磁界Biは、地上座標系ΣGにおいて、携帯機器1の姿勢μの変化に連動して向きを変える、一定の大きさのベクトルGBiμとして表される。ここで、携帯機器1厳密には3次元磁気センサ70に固定された座標系であるセンサ座標系ΣSを導入し、磁気データqを、センサ座標系ΣS上の座標として表現する。より具体的には、磁気データqが示すセンサ座標系ΣS上の座標のうち、x軸成分は、X軸磁気センサ71からの出力値であり、y軸成分は、Y軸磁気センサ72からの出力値であり、z軸成分は、Z軸磁気センサ73からの出力値である。このようなセンサ座標系ΣSから見た場合、地磁気Bgは、携帯機器1の姿勢μに連動して向きを変える、一定の大きさのベクトルSBgμとして表され、内部磁界Biは、一定の向き及び大きさを有するベクトルSBiとして表される。従って、地上座標系ΣGから見た携帯機器1の姿勢μをμ1~μNと変化させつつ取得した複数の磁気データq1~qNは、センサ座標系ΣSにおいて、図3Bに示すように、内部磁界Biを表すベクトルSBiの示す座標を中心とし、地磁気Bgの大きさを半径とする、球面S上に分布する。そして、磁気データqの示す座標から、内部磁界Biを表すベクトルSBiの示す座標を減算することにより、センサ座標系ΣSにおける地磁気Bgの向きが特定される。このように、検出対象である地磁気Bgの正確な向きを特定するために、補正処理において、3次元磁気センサ70から出力される磁気データqから取り除かれる内部磁界Biを表すベクトルを磁気センサのオフセットcOFFと呼ぶ。なお、複数の磁気データq1~qNの示す座標は、3次元磁気センサ70の測定誤差を考慮すると、球面S上には分布せず、球面Sの近傍に分布する。従って、複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの座標は、内部磁界Biを表すベクトルSBiの示す座標の近似値すなわち、オフセットcOFFの近似値となる。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を実行することにより、状態推定装置として機能する。状態推定プログラム100は、バッファ管理モジュール200、中心点算出モジュール300、安定性判定モジュール400、初期ベクトル生成モジュール500、及びカルマンフィルタモジュール600等のモジュール群で構成される。バッファ管理モジュール200は、3次元磁気センサ70から順次出力される磁気データqをバッファに蓄積する。これら複数の磁気データの蓄積手段としては、RAM20を用いる。中心点算出モジュール300は、RAM20に蓄積された複数の磁気データq1~qNに基づいて、センサ座標系ΣSにおいて複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの座標を算出する。安定性判定モジュール400は、3次元加速度センサ80から出力される加速度データaに基づいて、携帯機器1の動きが安定しているか否かを判定する。初期ベクトル生成モジュール500は、安定性判定モジュール400による判定結果が肯定である場合に、初期ベクトルINIを生成する。カルマンフィルタモジュール600は、初期ベクトルINIを初期値とする状態ベクトルxを、観測値ベクトルyに基づいて更新する演算非線形カルマンフィルタの演算を実行することで、システムの状態を推定する。ここで、状態ベクトル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次元数である。例えば、センサ座標系ΣSの各軸と、地上座標系ΣGの各軸とが一致する姿勢μを基準姿勢と定義し、基準姿勢をμ=0、0、0、1と表現する。このとき、携帯機器1の任意の姿勢μは、基準姿勢に対して単位ベクトルρの方向を回転軸として角度ψだけ回転した姿勢として表現できる。回転後の姿勢μは、クォータニオンを用いて、式2で与えられる。観測値ベクトルyとは、3次元磁気センサ70から出力される磁気データq、3次元加速度センサ80から出力される加速度データa、及び、3次元角速度センサ90から出力される角速度データgを要素とするベクトルである。時刻T=kにおける観測値ベクトルykは式3で与えられる。ところで、非線形カルマンフィルタの演算は、状態ベクトルと観測値ベクトルとに基づいて算出される観測残差を用いて、状態ベクトルを周期的に更新する演算であるため、状態ベクトルの初期値の各要素各状態変数の示す値が、真値から乖離した値である場合、状態ベクトルの各要素を真値に近い値に収束させるのに長い時間を要することがある。また、非線形カルマンフィルタの演算は、非線形演算であるため、状態ベクトルの初期値の各要素が真値から大きく乖離した値である場合、状態ベクトルの更新を繰り返しても、状態ベクトルの各要素が真値または、真値に近い値に収束せずに、真値とは異なる値に収束することがある。従って、非線形カルマンフィルタの演算において、システムの状態を正確且つ高速に推定するためには、状態ベクトルxの初期値つまり初期ベクトルINIの各要素を、真値に近い値に設定する必要がある。ここで、真値とは、状態ベクトルxkの各要素に対応する実際の物理量を表す値である。また、状態ベクトルの真値とは、状態ベクトルxkの各要素を真値としたときのベクトルである。[2.状態推定装置の処理フロー]図4は、CPU10が実行する状態推定プログラム100の処理の流れを表すフローチャートである。以下では、CPU10が状態推定プログラム100を実行したときの処理の流れを、図4に加えて、当該フローチャートの各ステップが実行された時に実現される機能の詳細を説明する図5乃至図6を参照しつつ、説明する。まず、ステップS101において、CPU10は、初期化処理を実行する。具体的には、CPU10は、バッファ管理モジュール200を呼び出すことで、RAM20に記憶した磁気データqの全部を廃棄する。なお、本実施形態では、CPU10は、磁気データqの全部を廃棄するが、RAM20に蓄積されたデータのうち古い方から一定割合のデータのみを廃棄しても良い。ステップS102において、CPU10は、第1の報知を行う。具体的には、CPU10は、表示部50に画像や動画等を表示させることで、携帯機器1の利用者に対して、携帯機器1の姿勢を3次元的に変化させるように指示する。すなわち、CPU10は、ステップS102の第1の報知を行うことにより、第1の報知手段として機能する。なお、本実施形態では、CPU10は、第1の報知を、表示部50を用いて画像や動画等を表示することにより行うことを想定しているが、音声等を用いて行ってもよい。ステップS103において、CPU10は、磁気データ蓄積処理を実行する。具体的には、CPU10は、バッファ管理モジュール200を呼び出すことにより、3次元磁気センサ70から順次出力されるN個の磁気データq1~qNを、RAM20に格納するここで、Nは、後述する中心点算出処理において、精度よく中心点cSを導出するために必要な磁気データの規定測定回数を表す5以上の自然数。すなわち、CPU10は、ステップS103の磁気データ蓄積処理を実行することにより、蓄積手段として機能する。次に、ステップS104において、CPU10は、中心点算出処理を実行する。具体的には、CPU10は、中心点算出モジュール300を呼び出すことにより、RAM20に蓄積された複数の磁気データq1~qNの示す座標を近傍に有する球面Sの中心点cSの座標を算出する。すなわち、CPU10は、ステップS104の中心点算出処理を実行することにより、中心点算出手段として機能する。なお、前述の通り、中心点cSの座標を示すベクトルは、オフセットcOFFの近似値である。中心点cSの座標の算出は、公知の方法を適宜適用して行えばよい。例えば、本件出願人が、特開2007-240270号公報で開示している方法等を適用してもよい。次に、ステップS105において、CPU10は、第2の報知を行う。具体的には、CPU10は、表示部50に画像や動画等を表示させることで、携帯機器1の利用者に対して、携帯機器1を動かさないように指示する。すなわち、CPU10は、ステップS105の第2の報知を行うことにより、第2の報知手段として機能する。なお、本実施形態では、CPU10は、第2の報知を、表示部50を用いて画像や動画等を表示することにより行うことを想定しているが、音声等を用いて行ってもよい。その後、ステップS106において、CPU10は、安定性判定処理を実行する。具体的には、CPU10は、安定性判定モジュール400を呼び出すことにより、3次元加速度センサ80から出力される加速度データaに基づいて、携帯機器1の動きが安定しているか否かを判定する。すなわち、CPU10は、ステップS106の安定性判定処理を実行することにより、安定性判定手段として機能する。図5は、CPU10が、ステップS106の安定性判定処理を実行することで実現される安定性判定手段の機能を表す機能ブロック図である。図5に示すように、安定性判定手段には、平滑化手段410、安定性指標算出手段420、及び条件判定手段430が含まれる。平滑化手段410は、3次元加速度センサ80から順次出力される加速度データaの各成分を、ローパスフィルタを通して平滑化することで、基準加速度データaCを算出する。ここで、ローパスフィルタには、公知のフィルタを適宜適用すればよい。例えば、FIRFiniteimpulseresponseフィルタ、または、IIRInfiniteimpulseresponseフィルタを用いてもよい。なお、本実施形態では、平滑化手段410への加速度データaの入力の開始は、CPU10がステップS106の安定性判定処理を開始するタイミングであるが、それ以外のタイミングで入力を開始してもよい。例えば、CPU10がステップS101の初期化処理を開始するタイミングで、平滑化手段410への加速度データaの入力を開始してもよい。この場合、CPU10は、ステップS101の初期化処理を開始するタイミングで、ステップS106の安定性判定処理を開始して、平滑化手段410を機能させる。次に、安定性指標算出手段420は、3次元加速度センサ80から出力される加速度データaと、平滑化手段410から出力される基準加速度データaCとに基づいて、安定性指標αを算出する。時刻tにおいて算出される安定性指標αtは、以下の式4aで表される。なお、式4bに示される加速度データatは、時刻tにおける加速度データであり、式4cに示される基準加速度データaCtは、時刻tにおける基準加速度データである。安定性指標算出手段420は、時刻t=1から時刻t=Mまでの期間に、M個の安定性指標αtt=1、...、Mを算出する。ここで、Mは、1以上の任意の自然数である。なお、携帯機器1の動きが安定しているか否かを、より厳密に判定することが必要な場合には、Mを、できるだけ大きな自然数に設定すればよい。条件判定手段430は、以下の式5に示すように、M個の安定性指標αtt=1、...、Mの全てが、第1閾値αTH以下であるか否かを判定する。このように、CPU10は、ステップS106の安定性判定処理を実行することにより、まず、M個の安定性指標αtを算出し、次に、M個の安定性指標αtの全てが第1閾値αTH以下であるか否かを判定する。そして、CPU10は、安定性判定処理における判定結果が肯定である場合、処理をステップS107に進める。一方、CPU10は、安定性判定処理における判定結果が否定である場合、処理をステップS105に戻す。次に、ステップS107において、CPU10は、初期ベクトル生成処理を実行する。具体的には、CPU10は、初期ベクトル生成モジュール500を呼び出すことにより、時刻t=Mつまり、条件判定手段430の判定結果が肯定となったタイミングに3次元磁気センサ70から出力される磁気データq、時刻t=Mに3次元加速度センサ80から出力される加速度データa、及び、ステップS104において算出される中心点cSの座標に基づいて、初期ベクトルINIを生成する。すなわち、CPU10は、ステップS107の初期ベクトル生成処理を実行することにより、初期ベクトル生成手段として機能する。前述の通り、初期ベクトルINIは、状態ベクトルxkの初期値、つまり時刻T=0における状態ベクトルx0であり、姿勢μの初期値μ0、地磁気の強さrの初期値r0、地磁気の伏角φの初期値φ0、角速度ωの初期値ω0、角速度センサのオフセット推定値bOの初期値bO0、及び磁気センサのオフセット推定値cOの初期値cO0を要素として含む、以下の式6により示されるベクトルである。なお、以下では、CPU10は、初期ベクトル生成処理を実行することにより初期ベクトルINIを生成するタイミングを、時刻T=0と表現する。時刻T=0は、時刻t=Mと等しい時刻である。以下、初期ベクトルINIの生成方法を、具体的に説明する。まず、時刻T=0に3次元磁気センサ70から出力される磁気データq0は、地上座標系ΣGから地磁気Bgを表したベクトルGBgと、3次元磁気センサ70のオフセット推定値cOの初期値cO0とを用いて、以下の式7で表される。ここで、行列Bμは、携帯機器1が姿勢μである場合に、地上座標系ΣGにおいて表現されたベクトルを、センサ座標系ΣSで表現するための座標変換行列であり、以下の式8で表される。また、時刻T=0における地磁気Bgを地上座標系ΣGから表したベクトルGBgは、地磁気の強さrの初期値r0と、地磁気の伏角φの初期値φ0とを用いて、以下の式9で表される。なお、本実施形態では、地上座標系ΣGを、y軸が磁極北を向き、z軸が鉛直上を向くように定める。また、3次元磁気センサ70のオフセット推定値cOの初期値cO0には、以下の式10に示すように、中心点cSの座標が設定される。一方、時刻T=0に3次元加速度センサ80から出力される加速度データa0は、携帯機器1の動きが安定していれば例えば、携帯機器1が静止状態または等速直線運動状態にあれば、以下の式11で表される。なお、ベクトルGGRVは、地上座標系ΣGにおいて重力加速度を表したベクトルを、重力加速度の大きさで正規化した3次元のベクトルであり、以下の式12で表される。センサ座標系ΣSにおいて地磁気Bgの向き及び大きさを表すベクトルSBgは、磁気データqと3次元磁気センサ70のオフセット推定値cOとを用いて、以下の式13aによって表すことができる。携帯機器1の姿勢は、以下の式13aに示すベクトルSBgと、加速度データaとを用いて、以下の式13bに示す3行3列の行列Aとして表現することができる。そして、行列Aにより表された携帯機器1の姿勢を、クォータニオンを用いて表現し直すことにより、携帯機器1の姿勢μを算出することができる。つまり、姿勢μの初期値μ0は、以下の式13a及び式13bに対して、時刻T=0に3次元加速度センサ80から出力される加速度データa0と、磁気センサのオフセット推定値cOの初期値cO0とを適用することにより算出される行列Aに基づいて、求めることができる。なお、3行3列の行列Aからクォータニオンへの変換は、公知の方法を適宜用いて行えば良い。行列Aは、式8に示した行列Bμの逆行列であるので、行列Aの各成分と行列Bμの逆行列の各成分とを比較することで、クォータニオンを用いて表される姿勢μの各要素を計算すればよい。例えば、以下の公知文献に記載された方法を用いて、行列Aから姿勢μを算出することもできる。MalcolmDShusterandGregoryANatansonQuaternionComputationfromaGeometricPointofViewTheJournaloftheAstronauticalSciencesVol41No4OctoberDecember1993pp545556地磁気の強さrの初期値r0は、以下の式14aに基づいて算出することができる。具体的には、地磁気の強さrの初期値r0は、まず、式13aに対して磁気データqと磁気センサのオフセット推定値cOの初期値cO0とを代入してベクトルSBgを求め、次に、ベクトルSBgを以下の式14aに代入することで算出される。また、地磁気の伏角φの初期値φ0は、ベクトルSBgの各要素を以下の式14bで表した場合、以下の式14c及び式14dに基づいて算出することができる。角速度ωの初期値ω0は、例えば、携帯機器1が静止していることを仮定して、0に設定する。また、角速度センサのオフセット推定値bOの初期値bO0は、通常0近辺に調整されているため、0に設定する。このように、CPU10は、ステップS107の初期ベクトル生成処理を実行することにより、姿勢μの初期値μ0、地磁気の強さrの初期値r0、地磁気の伏角φの初期値φ0、角速度ωの初期値ω0、角速度センサのオフセット推定値bOの初期値bO0、及び磁気センサのオフセット推定値cOの初期値cO0を要素とするベクトルである初期値INIを生成し、これを出力する。次に、ステップS108において、CPU10は、カルマンフィルタ演算処理を実行する。具体的には、CPU10は、カルマンフィルタモジュール600を呼び出すことにより、非線形カルマンフィルタの演算を行い、システムの状態を推定する。すなわち、CPU10は、ステップS108のカルマンフィルタ演算処理を実行することにより、カルマンフィルタ演算手段として機能する。以下、非線形カルマンフィルタの演算の詳細な内容について説明する。[3.非線形カルマンフィルタによる演算]一般的にカルマンフィルタは、動的システムの状態の経時的な変化を表す状態遷移モデルを用いて、ある時刻時刻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-1を更新することにより、更新後の状態ベクトルxkを算出する。以上のような状態ベクトルxkの更新を繰り返すカルマンフィルタの演算により、状態ベクトルxkの各要素を実際の物理量を正確に表した値真値に近い正確な値へと近付ける。なお、詳細は後述するが、本実施形態ではカルマンフィルタとして、非線形カルマンフィルタであるシグマポイントカルマンフィルタを用いる。シグマポイントカルマンフィルタは、まず、状態ベクトルxk-1の周囲に複数のシグマポイントχk-1を設定し、これらの複数のシグマポイントχk-1を状態遷移モデルに適用することで、単位時間経過後の複数のシグマポイントを推定する。この推定された複数のシグマポイントを、推定シグマポイントχ-kと称する。次に、シグマポイントカルマンフィルタは、複数の推定シグマポイントχ-kの平均を算出することにより、推定状態ベクトルx-kを求める。従って、推定観測値ベクトルy-kは、厳密には、推定状態ベクトルx-kの周辺に存在する複数の推定シグマポイントχ-kに基づいて算出される。本実施形態における、非線形カルマンフィルタの状態遷移モデルは以下に示す式15で与えられ、観測モデルは以下に示す式16で与えられる。なお、本実施形態では、式15に現れる関数f及び式16に現れる関数hは、非線形の関数である。ここで、状態ベクトルxkはdimx次元のベクトルであり、観測値ベクトルykはdimy次元のベクトルである。なお、本実施形態では、dimx=15であり、dimy=9である。また、式15に現れるプロセスノイズwk及び式16に現れる観測ノイズvkは0を中心とするガウスノイズである。式15は、時刻T=kにおける状態ベクトルxkが、時刻T=k-1における状態ベクトルxk-1を関数fに代入して得られた値と、時刻T=k-1におけるプロセスノイズwk-1とを加算することにより推定されることを示している。また、式16は、時刻T=kにおける観測値ベクトルykが、時刻T=kにおける状態ベクトルxkを関数hに代入して得られる値と、時刻T=kにおける観測ノイズvkとを加算することにより推定されることを示している。なお、プロセスノイズwkの共分散をQk、観測ノイズvkの共分散をRkと表す。時刻T=kにおける観測残差zkは、実際の観測値ベクトルykと、推定観測値ベクトルy-kとに基づいて定められるベクトルであり、以下に示す式17で表される。カルマンフィルタは、以下の式19に示すように、観測残差zk、推定状態ベクトルx-k、及び、式18に示すカルマンゲインKkを用いて、更新後の状態ベクトルxkを算出する。また、カルマンフィルタは、以下の式20に示すように、状態ベクトルxkの推定誤差の共分散Pkを更新する。ここで、P-kは、推定状態ベクトルx-kの推定誤差の共分散であり、Pkは、更新後の状態ベクトルxkの推定誤差の共分散である。また、Pyykは、観測残差zkの共分散であり、Pxykは、推定状態ベクトルx-kと、推定観測値ベクトルy-kとの相互共分散行列である。なお、推定観測値ベクトルy-kは、式16に示すように、状態遷移モデルを用いて算出された推定状態ベクトルx-k厳密には、推定シグマポイントχ-kを、観測モデルに適用することで算出される値である。よって、推定観測値ベクトルy-kと実際のセンサ出力に基づく観測値ベクトルykとの差分である観測残差zkは、推定状態ベクトルx-kと実際の物理量を正確に表した値真値との近似度を示す値である。式19に示すように、観測残差zkを用いて、推定状態ベクトルx-kを、更新後の状態ベクトルxkへと更新することにより、状態ベクトルxkの各要素を真値に近い値へと近付けることができる。本実施形態に係る状態推定装置は、アンセンテッド変換を用いたシグマポイントカルマンフィルタを用いて、カルマンフィルタ演算処理を実行する。図6は、CPU10がカルマンフィルタモジュール600を呼び出すことで実現されるカルマンフィルタ演算手段の機能を表す機能ブロック図である。図6に示すように、遅延部610は、加算器690から出力される更新後の状態ベクトルxkを、単位時間時刻T=k-1から時刻T=kに相当する時間だけ遅延させることで、状態ベクトルxk-1を生成し、これを、シグマポイント生成部620に対して出力する。なお、初回の演算時刻T=0では、遅延部610は、初期ベクトル生成モジュール500が出力する初期値INIを、状態ベクトルxk-1に適用する。シグマポイント生成部620は、dimx行dimx列の共分散行列Pk-1及びQk-1を用いて、2dimx1個のシグマポイントを生成する。具体的には、まず、複数のシグマポイントの平均に対して、平均のまわりのシグマポイントの広がりを表すスケーリングパラメータλを用いて、式21及び式22に示すベクトルσkを定義する。このとき、シグマポイント生成部620は、ベクトルσk-1と状態ベクトルxk-1とに基づいて、式23及び式24で示される2dimx1個のシグマポイントχk-1を生成する。状態遷移モデル部630は、式25に示すように、時刻T=k-1における2dimx1個のシグマポイントχk-1jの各々を、状態遷移モデルに適用することにより、時刻T=kにおける2dimx1個の推定シグマポイントχ-kjを算出する。次に、平均算出部640は、式26に示すように、時刻T=kにおける2dimx1個の推定シグマポイントχ-kの平均値を計算することで、推定状態ベクトルx-kを算出する。また、平均算出部640は、式27に示す、推定状態ベクトルx-kの誤差の共分散P-kを算出する。一方、観測モデル部650は、式28に示すように、時刻T=kにおける2dimx1個の推定シグマポイントχ-kjの各々を、観測モデルに適用することにより、2dimx1個の推定観測値γkjを算出する。そして、平均処理部660は、式29に示すように、2dimx1個の推定観測値γkjの平均を演算することにより、推定観測値ベクトルy-kを算出する。次に、平均処理部660は、式30に示す観測残差の共分散Pyykを算出する。減算器670は、式17に示したように、観測値ベクトルykと、推定観測値ベクトルy-kとの差分として、観測残差zkを算出する。カルマンゲイン付与部680は、式31に示す相互共分散行列Pxykを算出する。そして、カルマンゲイン付与部680は、式18に示したように、残差の共分散Pyykと、相互共分散行列Pxykとに基づいて、カルマンゲインKkを算出し、式19の右辺第2項Kkzkの演算を実行する。また、カルマンゲイン付与部680は、式20に示したように、状態ベクトルxkの推定誤差の共分散Pkを、P-kからPkに更新する。加算器690は、式19に示したように、推定状態ベクトルx-kと、カルマンゲイン付与部680から出力される式19の右辺第2項Kkzkとを加算することにより、更新後の状態ベクトルxkを算出する。次に、状態遷移モデル部630における演算で用いる状態遷移モデルについて説明する。本実施形態に係る状態遷移モデルにおいて、状態ベクトルxkを構成する状態変数のうち、姿勢μについての状態遷移は、式32により定義される。式32は、時刻T=k-1における姿勢μk-1から、単位時間経過後の時刻T=kにおける姿勢μ-kを推定する演算を表す。ここで、μ-kは、時刻kにおける推定状態ベクトルx-kのうち、姿勢μを表す状態変数に対応する要素である。なお、式32の右辺の演算子Ωは、式33により定義される。ここで、I3×3は3行3列の単位行列を表す。3次元ベクトルl=l1,l2,l3に対して、演算子[l×]は、式34で定義される。また、測定時間間隔時刻T=k-1から時刻T=kまでの間隔をΔtで表し、時刻T=kにおける更新後の状態ベクトルxkのうち角速度を表す状態変数に対応する要素をωkで表したとき、演算子Ωを構成する成分Ψkは、式35で定義される。姿勢μは、クォータニオンで表現され、式2に示すように、正規化条件μ=1が満たされる必要がある。しかし、シグマポイントカルマンフィルタを用いて状態ベクトルxkを更新する場合、更新後の状態ベクトルxkは、式26に示すように、推定シグマポイントχ-kjの平均として算出されるため、更新前の状態ベクトル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とは等しいと仮定する。前述の通り、磁気センサのオフセットcOFFは、携帯機器1の部品が発する内部磁界Biの方向及び大きさをセンサ座標系ΣSにおいて表現したベクトルである。従って、携帯機器1の内部状態が一定である間は、磁気センサのオフセットcOFFも変化しない。一方、携帯機器1の内部状態が変化した場合、例えば、携帯機器1に搭載された部品を流れる電流の大きさが変化した場合や、携帯機器1に搭載された部品の着磁状況が変化した場合には、磁気センサのオフセットcOFFも変化する。すなわち、携帯機器1の内部状態は、携帯機器1の利用者による携帯機器1の操作や、携帯機器1の外部の環境等に依存して変化する。従って、磁気センサのオフセットcOFFが変化するタイミングや磁気センサのオフセットcOFFの変化量を予測することは困難であり、時刻T=k-1における磁気センサのオフセット推定値cOK-1を用いて、時刻T=kにおける磁気センサのオフセット推定値cOKを定式化することは難しい。そこで、本実施形態に係る状態遷移モデルでは、便宜上、時刻T=kにおける磁気センサのオフセット推定値cOKと、時刻T=k-1における磁気センサのオフセット推定値cOK-1とは等しいと仮定する。このように、本実施形態に係る状態遷移モデルでは、以下に示す式36のように、状態ベクトルxkを構成する複数の状態変数のうち、姿勢μを表す状態変数以外は、前の時刻から変化しないという前提を置く。なお、カルマンフィルタの演算においては、状態ベクトルxkは、観測残差zkに基づいて更新される。従って、式36に示された姿勢μ以外の値が、全く変化しないわけではない。具体的には、状態ベクトルxkの各要素は、加算器690において、観測残差zkに基づいて真値に近づくように更新される。次に、観測モデル部650における演算で用いられる観測モデルについて説明する。3次元角速度センサ90から出力される角速度データgの推定値γgyroは、角速度ωと、角速度センサのオフセット推定値bOとを用いて、式37で与えられる。また、地上座標系ΣGにおいて地磁気Bgを表すベクトルGBgは式38で与えられる。従って、3次元磁気センサ70から出力される磁気データqの推定値γmagは、磁気センサのオフセット推定値cOと、式8で示した行列Bμとを用いて、式39で与えられる。また、観測モデルにおいて、3次元加速度センサ80から出力される加速度データaの推定値γaccは、式8で示した行列Bμと、式12で示したベクトルGGRVとを用いて、式40で与えられる。式3を式17の右辺第1項に代入し、式37、式39、及び式40を用いて式17の右辺第2項を表すことにより、式17を、以下の式41に変形することができる。このとき、観測残差zkは、式41により算出される。[4.結論]以上に示したように、CPU10は、カルマンフィルタモジュール600を呼び出すことで、ステップS108のカルマンフィルタ演算処理を実行することにより、状態ベクトルxkを更新する。具体的には、時刻T=kにおける状態ベクトルxkは、時刻T=k-1における状態ベクトルxk-1を、状態ベクトルxk-1を用いて算出される観測残差zkに基づいて更新することで算出される。同様に、時刻T=k-1における状態ベクトルxk-1は、時刻T=k-2における状態ベクトルxk-2を、状態ベクトルxk-2を用いて算出される観測残差zk-1に基づいて更新することで算出される。すなわち、カルマンフィルタ演算処理において、状態ベクトルxkは、初期状態時刻T=0から時刻T=k-1に至る全ての状態ベクトルx0~xk-1の値を考慮して、更新される。従って、初期ベクトルINIの各要素が真値から大きく乖離した値に設定された場合、カルマンフィルタ演算処理を行っても、状態ベクトルxkの各要素は真値に近づかなくなることがあり、また、仮に状態ベクトルxkの各要素が真値に近づくとしても、真値に近い値に収束するまでに長い時間を要することになる。前述の通り、磁気センサのオフセットcOFFは、携帯機器1に搭載された部品が発する内部磁界Biを表すベクトルである。従って、携帯機器1に搭載された部品が大きな内部磁界Biを発する場合、磁気センサのオフセット推定値cOの初期値cO0に対して、事前に定められた固定値例えば、原点[0,0,0]T等を適用すると、磁気センサのオフセット推定値cOの初期値cO0が、磁気センサのオフセットcOFFから大きく乖離した値となることがある。また、姿勢μについても、予め定めた固定値を姿勢μの初期値μ0に設定する場合、真の姿勢から最大で180度ずれることになる。このように、初期ベクトルINIを構成する要素のうち、磁気センサのオフセット推定値cOの初期値cO0と、姿勢μの初期値μ0とは、真値から大きく乖離した値となることがある。従って、初期ベクトルINIの各要素が不適切な真値から乖離した値に設定されることを防止するためには、磁気センサのオフセット推定値cOの初期値cO0と、姿勢μの初期値μ0とを、真値に近い正確な値に設定することが、特に重要となる。本実施形態では、まず、中心点算出処理において、オフセットcOFFの近似値である中心点cSの座標を算出し、次に、初期ベクトル生成処理において、磁気センサのオフセット推定値cOの初期値cO0に、中心点cSの座標を適用することで、初期ベクトルINIを算出した。これにより、磁気センサのオフセット推定値cOの初期値cO0が、真値から大きく乖離した値に設定されることを防止した。また、本実施形態では、安定性判定処理において、携帯機器1の動きが安定していることを判定したうえで、姿勢μの初期値μ0を算出することで、姿勢μの初期値μ0が、真値から大きく乖離した値に設定されることを防止した。このように、本実施形態は、初期ベクトルINIの各要素を、真値に近い正確な値に設定するため、カルマンフィルタ演算処理において、正確且つ高速な状態推定が可能となった。<B.変形例>本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。1変形例1上述した実施形態では、安定性判定処理において、ローパスフィルタを用いて算出された基準加速度データaCと、加速度データaとに基づいて安定性指標αを算出したが、本発明はこのような形態に限定されるものでは無く、ローパスフィルタを用いずに安定性指標を算出してもよい。例えば、以下の式43で示される安定性指標βが、以下の式42で示される条件を満たすか否かにより、携帯機器1の動きが安定しているか否かを判定してもよい。ここで、安定性指標βは、以下の式43に示すように、時刻t=1から時刻t=Mまでの期間に3次元加速度センサ80が出力する複数の加速度データa1~aMの分散を表す値である。また、式42に示すように、安定性指標βが第2閾値βTH以下となり十分に小さな値となる場合、加速度データa1~aMの示す値の変動が小さいため、携帯機器1の動きが安定していると看做すことができるMは、加速度データaの分散を算出するために必要な加速度データの規定測定回数を表す自然数。なお、式43に現れる値axAVE、ayAVE、及びazAVEは、複数の加速度データa1~aMの各要素の平均値であり、それぞれ、式44、式45、及び式46により表される。2変形例2上述した実施形態及び変形例に係る状態推定装置は、CPU10が、ステップS101~S108の各処理を実行するものであったが、本発明はこのような形態に限定されるものでは無く、CPU10が、ステップS101~S108のうち一部の処理のみを実行するものであってもよい。例えば、状態推定装置は、第2の報知手段を備えず、CPU10が、ステップS105の第2の報知を実行しないものであってもよい。携帯機器1が、長期間継続的に動かされることは想定し難く、携帯機器1は、通常安定した状態に置かれる。そして、ステップS106の安定性判定処理で、携帯機器1が安定していることを確認したうえで、初期ベクトルINIが算出されるため、状態推定プログラム100が第2の報知手段を備えなくても、姿勢μの初期値μ0が、真値から大きく乖離した値に設定されることを防止することができる。また、例えば、状態推定装置は、第1の報知手段を備えず、CPU10が、ステップS102の第1の報知を実行しないものであってもよい。この場合、CPU10は、ステップS103の磁気データ蓄積処理において、蓄積された複数の磁気データq1~qNが、球面Sの中心点cSを算出するために適切であるか否かを判定し、適切であると判定した場合に、処理をステップS104へと進めてもよい。具体的には、CPU10は、ステップS103の磁気データ蓄積処理において、複数の磁気データq1~qNの各々がセンサ座標系ΣSにおいて示す座標の広がりが3次元的であるか否かを判定し、座標の広がりが3次元的であると判定した場合に、処理をステップS104へと進めてもよい。また、例えば、状態推定装置は、第1の報知手段及び第2の報知手段を備えず、CPU10が、ステップS102の第1の報知、及び、ステップS105の第2の報知を実行しないものであってもよい。3変形例3上述した実施形態及び変形例に係る状態推定装置において、CPU10は、カルマンフィルタ演算処理を、シグマポイントカルマンフィルタを用いて実行したが、本発明はこのような形態に限定されるものでは無く、公知の非線形カルマンフィルタを適宜適用して演算を行ってもよい。例えば、拡張カルマンフィルタを用いて、カルマンフィルタ演算処理を実行してもよい。"